脉冲响应与模拟：IRF、历史分解、方差分解与条件预测

% ======================================================================
% plot_irf.m  从 Dynare oo_.irf 画论文级 IRF 图
% 前提: 已运行 dynare nk_baseline.mod, 工作区有 oo_, M_
% ======================================================================
clear; clc; close all;

% --- 1) 取出 IRF: 假设冲击为 ea(技术), 变量为 y, pi, r ---
% Dynare 把 IRF 存在 oo_.irf 下, 字段名 = 变量名_冲击名
horizon = 1:40;                 % 横轴: 40 期
y_irf   = oo_.irf.yea;          % 技术冲击 ea 对 y 的 IRF
pi_irf  = oo_.irf.piea;         % 技术冲击 ea 对 pi 的 IRF
r_irf   = oo_.irf.rea;          % 技术冲击 ea 对 r 的 IRF

% --- 2) 贝叶斯后验: 取 16/50/84 分位数画置信带 ---
% 若用估计命令, 后验 IRF 存在 oo_.PosteriorIRF.metropolis.HPDinf / HPDsup
try
  y_lo = oo_.PosteriorIRF.metropolis.HPDinf.yea;
  y_hi = oo_.PosteriorIRF.metropolis.HPDsup.yea;
catch
  y_lo = y_irf*0.8;  y_hi = y_irf*1.2;   % 校准模型无置信带, 用 ±20% 示意
end

% --- 3) 画三行一列子图 ---
figure('Position',[100 100 700 750]);

subplot(3,1,1);
  fill([horizon fliplr(horizon)], [y_lo fliplr(y_hi)], ...
       [0.85 0.9 0.95],'EdgeColor','none'); hold on;
  plot(horizon, y_irf, 'b-','LineWidth',1.8);
  yline(0,'k:');  title('Output \hat{y}_t');
  ylabel('% dev. from SS');  grid on;

subplot(3,1,2);
  plot(horizon, pi_irf, 'r-','LineWidth',1.8);
  yline(0,'k:');  title('Inflation \hat{\pi}_t');
  ylabel('annualized pp');  grid on;

subplot(3,1,3);
  plot(horizon, r_irf, 'g-','LineWidth',1.8);
  yline(0,'k:');  title('Nominal Rate \hat{r}_t');
  xlabel('Quarters after shock'); ylabel('annualized pp'); grid on;

% --- 4) 导出为论文用 PDF/PNG ---
exportgraphics(gcf,'irf_tech_shock.pdf','ContentType','vector');
exportgraphics(gcf,'irf_tech_shock.png','Resolution',300);
disp('IRF 图已导出 irf_tech_shock.pdf / .png');

% ======================================================================
% 进阶: 把两个模型(无摩擦 vs BGG)叠在一张图对比
% plot(horizon, y_irf_noFric,'b-', horizon, y_irf_BGG,'r--','LineWidth',1.6);
% legend('No friction','BGG accelerator','Location','best');
% ======================================================================