适应性学习与行为预期
Bounded Rationality / Adaptive Learning

"""
E-stability 数值检验: 对简化 reduced form y_t = mu + alpha E*y_{t+1} + delta w_t
PLM: y_t = a + b w_t,  T(a)=mu+alpha*a,  T(b)=delta+alpha*rho*b
REE: a*=mu/(1-alpha),  b*=delta/(1-alpha*rho)
E-stability: DT = diag(alpha, alpha*rho) 的特征值实部 < 1
"""
import numpy as np

def check_estability(alpha, rho, verbose=True):
    DT = np.array([[alpha, 0.0],
                   [0.0,    alpha*rho]])
    eig = np.linalg.eigvals(DT)
    stable = np.all(np.real(eig) < 1.0)
    if verbose:
        print(f"alpha={alpha:.3f}, rho={rho:.3f}")
        print(f"  DT 特征值: {eig}")
        print(f"  最大实部 {np.max(np.real(eig)):.4f}  -> E-stable? {stable}")
    return stable, eig

# 教科书 NK 校准: Taylor 原则 phi_pi>1 对应 alpha 较小
for alpha in [0.5, 0.9, 1.1]:
    check_estability(alpha, rho=0.8)
# 输出: alpha=0.5 与 0.9 时 E-stable; alpha=1.1 时 E-unstable
# 经济含义: 央行反应越温和(phi_pi 接近 1), 前向系数 alpha 越大,
#           学习动态越可能发散 => 通胀脱锚风险.

"""
适应性学习（Adaptive Learning）最小可运行模拟
模型: y_t = mu + alpha * E* y_{t+1} + delta * w_t + eps_t,  w_t = rho*w_{t-1} + eta_t
PLM : y_t = a + b * w_t,  即 x_t = [1, w_t]', theta = [a, b]
更新: RLS (gamma_t = 1/(t+1)) 或 Constant Gain (gamma = gbar)
"""
import numpy as np
import matplotlib.pyplot as plt

# ---------- 1. 参数 ----------
alpha   = 0.90     # 前向系数 (对应 NK 模型中 phi_pi 较强时的值)
rho     = 0.80     # 外生冲击持续性
delta   = 0.50     # 外生冲击传导
mu      = 1.0      # 常数项
sig_e   = 0.10     # 内生冲击标准差
sig_w   = 0.10     # 外生冲击标准差
T       = 2000     # 模拟期数
gbar    = 0.02     # 恒定增益

# 解析的 REE（理性预期均衡）
a_star = mu / (1 - alpha)
b_star = delta / (1 - alpha * rho)
print(f"REE: a*={a_star:.4f}, b*={b_star:.4f}")

def run_learning(gain_mode="decreasing", seed=42):
    """跑一次学习模拟, 返回 (theta 路径, y 路径, w 路径)"""
    rng = np.random.default_rng(seed)
    theta = np.array([0.0, 0.0])     # 初始信念: 全部为 0 (远离 REE)
    R     = np.eye(2) * 10.0         # 初始信息矩阵 (伪逆保护)
    w     = 0.0
    a_path = np.zeros(T); b_path = np.zeros(T)
    y_path = np.zeros(T)
    for t in range(1, T):
        # --- (a) 外生状态 ---
        w = rho * w + sig_w * rng.standard_normal()
        # --- (b) 代理人基于 PLM 形成 E* y_{t+1} ---
        # PLM: y_t = theta' x_{t-1},  E*[y_{t+1}|I_t] = theta0 + theta1 * E*[w_{t+1}]
        # 代理人知道 w 的 AR(1) 系数 rho
        Et_y_next = theta[0] + theta[1] * (rho * w)
        # --- (c) ALM: 真实 y_t 由结构方程决定 ---
        eps = sig_e * rng.standard_normal()
        y   = mu + alpha * Et_y_next + delta * w + eps
        y_path[t] = y
        # --- (d) RLS / CG  更新 theta 与 R ---
        x_lag = np.array([1.0, w])          # 回归元 [1, w_t]
        gamma = 1.0/(t+1) if gain_mode == "decreasing" else gbar
        R     = R + gamma * (np.outer(x_lag, x_lag) - R)
        err   = y - theta @ x_lag           # 预测误差
        theta = theta + gamma * np.linalg.solve(R, x_lag) * err
        a_path[t] = theta[0]; b_path[t] = theta[1]
    return a_path, b_path, y_path

# ---------- 2. 跑两种增益 ----------
a_dec, b_dec, y_dec = run_learning("decreasing")
a_cg,  b_cg,  y_cg  = run_learning("constant")

# ---------- 3. 画收敛路径 ----------
fig, axes = plt.subplots(2, 2, figsize=(11, 7))
for ax, path, star, name in zip(
    axes[0],
    [a_dec, b_dec], [a_star, b_star], ["a (intercept)", "b (slope on w)"]):
    ax.plot(path, lw=1.2, label="decreasing gain (RLS)")
    ax.axhline(star, color="k", ls="--", lw=1, label=f"REE = {star:.3f}")
    ax.set_title(f"PLM coefficient: {name}"); ax.legend(fontsize=8)
for ax, path, star, name in zip(
    axes[1],
    [a_cg, b_cg], [a_star, b_star], ["a (intercept)", "b (slope on w)"]):
    ax.plot(path, lw=1.2, label="constant gain (g=0.02)")
    ax.axhline(star, color="k", ls="--", lw=1, label=f"REE = {star:.3f}")
    ax.set_title(f"PLM coefficient: {name}"); ax.legend(fontsize=8)
plt.tight_layout(); plt.savefig("learning_convergence.png", dpi=130); plt.show()

# 输出渐近样本均值, 看 CG 是否围绕 REE 波动
print(f"RLS  末段均值: a={a_dec[-200:].mean():.3f}, b={b_dec[-200:].mean():.3f}")
print(f"CG   末段均值: a={a_cg[-200:].mean():.3f}, b={b_cg[-200:].mean():.3f}")
print(f"CG   末段标准差: a={a_cg[-200:].std():.3f}, b={b_cg[-200:].std():.3f}  <-- 内生波动")

"""
NK 三方程下, 一次 1 单位成本推动冲击 u_t 的 IRF 对比
RE: 直接用解析的 RE 决策规则 pi_t = phi_pi^RE * u_t / (1 - beta*rho_u) 简化近似
Learning: 代理人用 RLS 学习 PLM pi_t = a + b*pi_{t-1} + c*x_{t-1}
"""
import numpy as np
import matplotlib.pyplot as plt

# NK 参数
beta, sigma, kappa, phi_pi, rho_u = 0.99, 1.0, 0.10, 1.5, 0.8
H = 40                       # IRF 长度
gbar = 0.03                  # 恒定增益

# ---- (A) 理性预期 IRF (解析近似: 把 IS 代入 Phillip 后) ----
# 在 Taylor 规则与 rho_u 下, RE 通胀对 u_t 的解析系数
pi_re = np.zeros(H); x_re = np.zeros(H)
coef_pi = kappa / (1 - beta*rho_u - sigma**(-1)*(phi_pi - 1)*kappa)
# 注: 这是教科书线性化解析解的简化形式, 教学用
u = np.zeros(H); u[0] = 1.0
for t in range(1, H):
    u[t] = rho_u * u[t-1]
    pi_re[t] = coef_pi * u[t]
    x_re[t] = (pi_re[t] - beta*pi_re[t-1]) / kappa

# ---- (B) 学习 IRF: 先做一段 burn-in 让信念接近 REE, 再在 t=0 注入冲击 ----
rng = np.random.default_rng(0)
theta = np.array([0.0, 0.9, 0.1])   # 初始信念: 高度外推
R = np.eye(3) * 5.0
pi_lrn = np.zeros(H); x_lrn = np.zeros(H)
pi_prev, x_prev = 0.0, 0.0
for t in range(H):
    # 代理人主观预期: E* pi_{t+1} = theta0 + theta1*pi_t + theta2*x_t
    Ep_next = theta[0] + theta[1]*pi_prev + theta[2]*x_prev
    Ex_next = Ep_next * 0.5          # 简化: 产出预期与通胀预期同向
    # 结构方程 (把 IS 与 Phillips 联立)
    u_t = u[t]
    # 简化求解: x_t 由 IS 给出, pi_t 由 Phillips 给出
    i_t = phi_pi * pi_prev
    x_t = Ex_next - sigma**(-1)*(i_t - Ep_next)
    pi_t = beta*Ep_next + kappa*x_t + u_t
    pi_lrn[t] = pi_t; x_lrn[t] = x_t
    # RLS 更新 (用 [1, pi_{t-1}, x_{t-1}] 回归 pi_t)
    x_lag = np.array([1.0, pi_prev, x_prev])
    R = R + gbar*(np.outer(x_lag, x_lag) - R)
    theta = theta + gbar*np.linalg.solve(R, x_lag)*(pi_t - theta@x_lag)
    pi_prev, x_prev = pi_t, x_t

fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].plot(pi_re,  "-o", label="Rational Expectations", lw=2)
ax[0].plot(pi_lrn, "-s", label="Adaptive Learning (CG=0.03)", lw=2)
ax[0].set_title("Inflation IRF to cost-push shock"); ax[0].legend()
ax[1].plot(x_re,  "-o", label="RE", lw=2)
ax[1].plot(x_lrn, "-s", label="Learning", lw=2)
ax[1].set_title("Output gap IRF"); ax[1].legend()
plt.tight_layout(); plt.savefig("irf_re_vs_learning.png", dpi=130)

"""
Coibion-Gorodnichenko 回归教学版
输入: 实际通胀序列 actual[t] 与调查预测序列 forecast[t]
输出: 误差对滞后误差的回归系数 beta_CG
"""
import numpy as np
from numpy.linalg import lstsq

def coibion_gorodnichenko_reg(actual, forecast):
    """actual[t]: t 期实际通胀; forecast[t]: t-1 期对 t 期通胀的预测
       返回 (alpha, beta_CG), 并打印 RE 是否被拒绝"""
    err = actual[1:] - forecast[1]      # 实际预测误差
    err_lag = err[:-1]                  # 滞后误差
    y = err[1:]
    X = np.column_stack([np.ones(len(y)), err_lag[:len(y)]])
    coef, *_ = lstsq(X, y, rcond=None)
    alpha, beta_CG = coef
    # 简单 t 检验 (OLS 标准误, 教学用)
    resid = y - X @ coef
    sigma2 = (resid @ resid) / (len(y) - 2)
    cov = sigma2 * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(cov))
    t_stat = beta_CG / se[1]
    print(f"alpha={alpha:.4f}, beta_CG={beta_CG:.4f} (t={t_stat:.2f})")
    print("RE 被拒绝 (beta_CG>0 显著)?" , abs(t_stat) > 1.96 and beta_CG > 0)
    return alpha, beta_CG

# 模拟: 恒定增益学习生成的预测 => CG 回归应显著为正
rng = np.random.default_rng(0)
T = 400
pi = np.zeros(T); pi[0] = 2.0
theta = 2.0  # 代理人对稳态通胀的初始信念
gbar = 0.03
for t in range(1, T):
    pi[t] = 0.9*pi[t-1] + 0.5*(theta - 0.9*pi[t-1]) + 0.2*rng.standard_normal()
    forecast_t = theta
    theta = theta + gbar*(pi[t] - theta)   # 简化 RLS
coibion_gorodnichenko_reg(pi[10:], np.roll(pi, 1)[10:])