稳态求解：解析推导、fsolve 数值解与 Dynare steady_state_model

% ============================================================
% solve_ss.m —— 用 fsolve 数值求解 baseline NK 模型稳态
% 内生变量 x = [Y; C; N; w; mc; r; i; pi]
% ============================================================
clear; clc;

% --- 1. 校准参数 ---
p.beta   = 0.99;
p.sigma  = 1.0;
p.phi    = 5.0;
p.alpha  = 0.25;
p.eps    = 9.0;
p.chi    = (1-p.alpha)*(p.eps-1)/p.eps;   % 归一化系数 -> N=1

% --- 2. 定义稳态方程组 ---
ss_eq = @(x) [
    x(1) - x(2);                          % (1) Y = C
    x(1) - x(3)^(1-p.alpha);              % (2) Y = N^(1-alpha), A=1
    p.chi*x(3)^p.phi * x(2)^p.sigma - x(4);   % (3) 劳动供给: chi*N^phi*C^sigma = w
    x(4) - (1-p.alpha)*x(5)*x(3)^(-p.alpha);  % (4) mc = w / ((1-alpha)*N^(-alpha))
    x(5) - (p.eps-1)/p.eps;               % (5) 稳态加成 mc = (eps-1)/eps
    1/p.beta - (1+x(6));                  % (6) 欧拉: 1 = beta*(1+r)
    x(7) - x(6) - x(8);                   % (7) Fisher: i = r + pi
    x(8)                                   % (8) 零通胀稳态 pi = 0
];

% --- 3. 初值 guess ---
x0 = [1; 1; 1; 1; (p.eps-1)/p.eps; 1/p.beta-1; 1/p.beta-1; 0];

% --- 4. fsolve 求解 ---
opts = optimoptions('fsolve','Display','iter','TolFun',1e-10);
[x, fval, exitflag] = fsolve(ss_eq, x0, opts);

% --- 5. 输出 ---
fprintf('\n===== 稳态解 =====\n');
fprintf('Y   = %.6f\n', x(1));
fprintf('C   = %.6f\n', x(2));
fprintf('N   = %.6f\n', x(3));
fprintf('w   = %.6f\n', x(4));
fprintf('mc  = %.6f\n', x(5));
fprintf('r   = %.6f (年化 %.4f%%)\n', x(6), 4*x(6)*100);
fprintf('i   = %.6f\n', x(7));
fprintf('pi  = %.6f\n', x(8));
fprintf('残差范数 = %.2e\n', norm(fval));

% ============================================================
% baseline_nk_ss.mod —— 完整可运行的 baseline NK 模型
% ============================================================

// ---------- 变量 ----------
var Y C N w mc r i pi y_gap;
varexo eps_a eps_m;

// ---------- 参数 ----------
parameters beta sigma phi alpha eps_p chi kappa rho_a rho_i phi_pi phi_y;

beta=0.99; sigma=1.0; phi=5.0; alpha=0.25;
eps_p=9.0;
chi = (1-alpha)*(eps_p-1)/eps_p;     % 归一化 -> N_ss = 1
kappa = (1-0.75)*(1-beta*0.75)/0.75 * (sigma+phi);  % theta=0.75
rho_a=0.9; rho_i=0.8; phi_pi=1.5; phi_y=0.125;

// ---------- 模型方程（对数线性化后）----------
model;
  y_gap = y_gap(+1) - (1/sigma)*r;
  r     = i - pi(+1);
  pi    = beta*pi(+1) + kappa*y_gap;
  i     = rho_i*i(-1) + (1-rho_i)*(phi_pi*pi + phi_y*y_gap) + eps_m;
end;

// ---------- 稳态（关键块）----------
steady_state_model;
  % --- 零通胀归一化稳态 ---
  pi   = 0;
  r    = 1/beta - 1;
  i    = r;
  N    = 1;
  Y    = N^(1-alpha);          % = 1
  C    = Y;
  w    = N^phi * C^sigma / chi;   % 注意归一化 chi
  mc   = (eps_p-1)/eps_p;
  y_gap = 0;
end;

// ---------- 冲击 ----------
shocks;
var eps_a; stderr 0.01;
var eps_m; stderr 0.0025;
end;

// ---------- 求解 ----------
steady;
check;
stoch_simul(irf=20, nograph);