零利率下限（ZLB）与偶发约束：OccBin 分段线性解与全局非线性求解

// ============================================================
// zlb_occbin.mod  ——  NK 三方程 + 名义利率零下限（OccBin）
// 运行：dynare zlb_occbin.mod
// 思想：两套 model 块分别对应 "约束松弛" 与 "约束绑定"
//       occbin 块声明 i >= i_lb，绑定时把 i 钉到下界
// ============================================================

var  y pi i rn ;
varexo eps_rn ;

parameters beta sigma varphi epsilon theta kappa ...
           phi_pi phi_y rho_rn i_lb;

// ---------- 深参数校准（与 NK 三方程页一致） ----------
beta    = 0.99;
sigma   = 1.0;
varphi  = 1.0;
epsilon = 6.0;
theta   = 0.75;
kappa   = (1-theta)*(1-beta*theta)/theta*(sigma+varphi)/(1+epsilon*varphi);
phi_pi  = 1.5;
phi_y   = 0.125;
rho_rn  = 0.90;
i_lb    = 0.0;          // 零下界（对数线性模型里取 0）

// ============================================================
// 【model 块 1】约束【松弛】：Taylor 规则完全生效
// ============================================================
model;
// (1) IS 曲线
y = y(+1) - (1/sigma)*( i - pi(+1) - rn );
// (2) NKPC
pi = beta*pi(+1) + kappa*y;
// (3) Taylor 规则（无下界的"影子"意愿利率）
i = phi_pi*pi + phi_y*y;
// (4) 自然利率 AR(1)
rn = rho_rn*rn(-1) + eps_rn;
end;

// ============================================================
// 【occbin 块】声明偶发约束：名义利率不能低于 i_lb
// ============================================================
occbin;
var i;                          // 被约束的变量
parameters i_lb;                // 下界参数
constraints;
  i >= i_lb;                    // 不等式约束
restrictions;
  i = i_lb;                     // 一旦绑定，把 i 钉到下界（替换 Taylor 规则）
end;

// ============================================================
// 【model 块 2】约束【绑定】：利率钉死在 i_lb
// 只把第三式换成 i = i_lb，其余与块1相同
// ============================================================
model(linear);
y = y(+1) - (1/sigma)*( i - pi(+1) - rn );
pi = beta*pi(+1) + kappa*y;
i = i_lb;                       // 关键区别：不再按 Taylor 规则调
rn = rho_rn*rn(-1) + eps_rn;
end;

// ---------- 初始值（正常稳态原点） ----------
initval;
 y = 0; pi = 0; i = 0; rn = 0;
end;
steady;

// ---------- 自然利率冲击（一次性大负冲击推入 ZLB） ----------
shocks;
var eps_rn; stderr = 0.01;
end;

// 给一个大的负自然利率冲击（-3%），观察经济如何撞入并走出 ZLB
// periods=80：模拟 80 个季度；replic：确定性模拟路径
shocks;
var eps_rn = (-0.03);   // 第 1 期一个 -3% 的自然利率冲击
end;

occbin_setup(display_occbin);
simul(periods=80);

// 想对比"无约束"版本：把 occbin 块与第二个 model 块删掉，
// 直接 stoch_simul 跑块1即可，会看到利率被允许降到负值（错误！）

# ============================================================
# zlb_global_solve.py  ——  NK三方程 + ZLB 的全局非线性求解
# 方法：在状态 rn 网格上，把每个网格点的 IS/NKPC/ZLB 残差
#       全部堆叠成一个大的非线性方程组，用 fsolve 一次性求根。
# 说明：纯前望 NK 在"深 ZLB 绑定区"局部特征根>1，朴素的
#       时间迭代会发散；直接联立求根能稳定收敛到全局政策函数。
# 状态：自然利率 rn（AR(1)，用 Rouwenhorst 离散）
# 控制：产出缺口 y, 通胀 pi, 名义利率 i
# 依赖：numpy, scipy
# ============================================================
import numpy as np
from scipy.optimize import fsolve

# ---------- 1. 深参数（与 Dynare 一致） ----------
beta    = 0.99      # 贴现因子
sigma   = 1.0       # 风险厌恶
kappa   = 0.0245    # NKPC 斜率（由深参数推出）
phi_pi  = 1.5       # Taylor 通胀系数
phi_y   = 0.125     # Taylor 产出系数
rho_rn  = 0.90      # 自然利率持续性
sd_rn   = 0.01      # 自然利率 innovation 标准差
i_lb    = 0.0       # 名义利率零下界（对数线性口径）

# ---------- 2. 用 Rouwenhorst 把 AR(1) 离散成状态网格 ----------
N = 21                              # 自然利率网格点数
p = (1 + rho_rn) / 2.0              # Rouwenhorst 相邻转移概率
def rouwenhorst(N, p, q):
    # 递归生成 N 点马尔可夫链转移矩阵 Pi
    Pi = np.array([[p, 1-p], [1-q, q]])
    for _ in range(N-2):
        a = np.zeros((Pi.shape[0]+1, Pi.shape[1]+1))
        a[:-1, :-1] += p*Pi
        a[:-1, 1:]  += (1-p)*Pi
        a[1:, :-1]  += (1-q)*Pi
        a[1:, 1:]   += q*Pi
        a[1:-1, :] /= 2
        Pi = a
    return Pi
Pi = rouwenhorst(N, p, p)
# 状态网格：对称中心 0，跨约 3.5 个无条件标准差
g = 3.5*sd_rn/np.sqrt(1-rho_rn**2)
rn_grid = np.linspace(-g, g, N)

# ---------- 3. 定义残差方程组（把 N 个网格点全部堆叠） ----------
# 未知向量 x = [y(0..N-1), pi(0..N-1)]，共 2N 个未知量
def residuals(x):
    y  = x[:N]          # 各网格点的产出缺口政策值
    pi = x[N:]          # 各网格点的通胀政策值
    # 期望项：E[y']=Pi@y, E[pi']=Pi@pi（对明天政策函数求期望，
    #         天然包含"将来可能撞 ZLB"的前向预期）
    Ey  = Pi @ y
    Epi = Pi @ pi
    r = np.zeros(2*N)
    for j in range(N):
        # 意愿（影子）利率与实际利率：ZLB 取 max(i_lb, 意愿)
        i_star = phi_pi*pi[j] + phi_y*y[j]
        i = i_star if i_star >= i_lb else i_lb   # 关键的分段 max
        # IS 曲线残差：y_j = Ey_j - (1/sigma)(i_j - Epi_j - rn_j)
        r[j]   = y[j] - Ey[j] + (1/sigma)*(i - Epi[j] - rn_grid[j])
        # NKPC 残差：pi_j = beta*Epi_j + kappa*y_j
        r[N+j] = pi[j] - beta*Epi[j] - kappa*y[j]
    return r

# ---------- 4. 初始猜测：无约束线性解附近 ----------
# 无约束时 y 对 rn 的导数约为 2.78（由深参数推出），pi 约 0.625*y
x0 = np.concatenate([2.78*rn_grid, 0.625*2.78*rn_grid])

# ---------- 5. 直接联立求根（全局政策函数） ----------
sol, info, ier, msg = fsolve(residuals, x0, full_output=True, xtol=1e-11)
y_policy  = sol[:N]
pi_policy  = sol[N:]
print(f"[求解] fsolve 收敛状态 ier={ier}")

# ---------- 6. 后处理：实际利率政策函数与绑定标记 ----------
i_star  = phi_pi*pi_policy + phi_y*y_policy
i_policy = np.maximum(i_lb, i_star)            # 实际利率（含 ZLB）
binding  = i_policy <= i_lb + 1e-9             # 哪些状态 ZLB 绑定

# ---------- 7. 打印关键结果 ----------
print("\n=== 全局政策函数（自然利率由低到高） ===")
print(f"{'rn':>8}{'y':>10}{'pi':>10}{'i(实际)':>10}{'ZLB绑定?':>10}")
for j in [0, N//4, N//2, 3*N//4, N-1]:
    print(f"{rn_grid[j]:8.4f}{y_policy[j]:10.4f}{pi_policy[j]:10.4f}"
          f"{i_policy[j]:10.4f}{str(binding[j]):>10}")
print("\n经济直觉：")
print("· rn 很负（深衰退）时 y 大幅为负、i 被钉在 0（ZLB 绑定）；")
print("· 因为央行降不动息，通缩+衰退被前向预期放大；")
print("· rn 回到 0 附近，约束松弛，i 随 Taylor 规则正常变动。")