TVP-VAR 与风险溯源
Time-Varying Parameter VAR

% ==============================================================
% TVP-VAR 估计核心骨架（Primiceri 2005 简化版）
% 生产环境请直接使用:
%   - Gary Koop & Dimitris Korobilis 的 TVP-VAR 工具箱
%     (https://sites.google.com/site/garykoop/)
%   - Primiceri 个人主页提供的 replication code
% 本文件只演示"如何组织数据 + 调用 FFBS + 保存后验"
% ==============================================================
clear; clc; rng(2024);

% ---------- 0. 读数据：n x T ----------
% y 三列: [inflation, unemployment, interest_rate]
load us_macro.mat;
y = us_macro;
[T, n] = size(y);
p = 2;                       % VAR 滞后阶数
m = n*(n*p + 1);             % beta_t 的维度

% ---------- 1. 堆叠 VAR 右侧 ----------
Y = y(p+1:end, :);
[T2, ~] = size(Y);
X = zeros(T2, m);
for t = 1:T2
    xt = 1;
    for lag = 1:p
        xt = [xt, y(p+t-lag, :)];
    end
    X(t,:) = kron(eye(n), xt);
end

% ---------- 2. 先验（Primiceri 2005 第 4 节） ----------
B0_hat  = (X'*X)\X'*Y(:);        % 用全样本 OLS 做初值
beta_0  = zeros(m,1);
Vbeta_0= 4*eye(m);               % 扩散先验
s2      = 4;                     % 超参数尺度
Q       = (s2+1)*eye(m)*0.01;    % 状态协方差初值

% ---------- 3. MCMC 主循环（示意 5000 次） ----------
nrep = 5000; nburn = 2000;
beta_draws = zeros(nrep, T2, m);
for irep = 1:nrep
    % (a) FFBS 抽 beta_t（给定 A, Sigma）
    %   beta_t 状态方程: beta_t = beta_{t-1} + u_t, Var=Q
    %   观测方程: Y = X*beta + A^-1*Sigma*eps
    beta_draws(irep,:,:) = FFBS_linear(Y, X, Q, beta_0, Vbeta_0);

    % (b) 抽 A_t（下三角同期结构）
    % A_draws(irep,:,:) = FFBS_A(...);

    % (c) 抽 Sigma_t（随机波动率，KSC 1998 mixture）
    % Sigma_draws(irep,:,:) = KSC_sv(Y, beta_draws(irep,:,:));

    % (d) 抽超参数 Q（Inverse-Wishart 后验）
    % Q = iwishrnd(...)

    if mod(irep,500)==0
        fprintf('MCMC iteration %d / %d\n', irep, nrep);
    end
end

% ---------- 4. 后处理：剔除 burn-in ----------
post_beta = beta_draws(nburn+1:end,:,:);
% 平均时变系数、后验分位数
beta_med = median(post_beta, [1 3]);

% ---------- 5. 计算时变脉冲响应 ----------
% 在每个时点 t 用后验 beta_t, A_t, Sigma_t 构造 VAR companion form,
% 向前迭代 h=0..20, 得到 IRF(t, h, shock)
% plot(0:20, squeeze(IRF(2008,:,3)))  % 2008 年货币政策冲击对产出
% plot(0:20, squeeze(IRF(2020,:,3)))  % 2020 年同一冲击

# ==============================================================
# 简化 TVP-VAR：时变系数 + 固定协方差
# 教学用途。完整 Primiceri 估计见上节 MATLAB 工具箱。
# 依赖: pip install numpy pandas statsmodels matplotlib
# ==============================================================
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)

# ---------- 0. 模拟一段数据：真实 beta_t 缓慢漂移 ----------
T, n, p = 400, 2, 2
# 真实系数随时间随机游走（演示"结构变化"）
beta_true = np.zeros((T, n*(n*p+1)))
beta_true[0] = np.tile([0.5, -0.2, 0.1, 0.3, 0.0, 0.2], 1)
for t in range(1, T):
    beta_true[t] = beta_true[t-1] + np.random.normal(0, 0.02, beta_true.shape[1])

y = np.zeros((T, n))
y[:p] = np.random.normal(0, 1, (p, n))
for t in range(p, T):
    xt = np.concatenate([[1.0], y[t-1], y[t-2]])      # 1 + 2*2 = 5 维
    # 注意: 这里演示用 n=2, m=5; 与上面 beta_true 维度示意略作简化
    y[t] = xt @ beta_true[t][:5].reshape(n, 5).T + np.random.normal(0, 0.3, n)

# ---------- 1. 构造观测矩阵 X ----------
m = n*(n*p + 1)
X = np.zeros((T-p, m))
Y = np.zeros((T-p, n))
for i, t in enumerate(range(p, T)):
    xt = np.concatenate([[1.0], y[t-1], y[t-2]])
    X[i] = np.tile(xt, n)
    Y[i] = y[t]

# ---------- 2. 卡尔曼滤波 (线性高斯状态空间) ----------
# 观测: Y_t = X_t * beta_t + eps,  eps ~ N(0, R)
# 状态: beta_t = beta_{t-1} + eta, eta ~ N(0, Q)
Q = np.eye(m) * 0.02**2          # 状态扰动方差（关键超参）
R = np.eye(n) * 0.3**2           # 观测噪声
beta_filt = np.zeros((T-p, m))
P = np.eye(m) * 10               # 初始状态方差
beta_tt = np.zeros(m)

for i in range(T-p):
    # predict
    beta_pred = beta_tt
    P_pred = P + Q
    # update
    S = X[i].reshape(1,-1) @ P_pred @ X[i].reshape(-1,1) + R
    K = P_pred @ X[i].reshape(-1,1) / S
    innovation = Y[i] - X[i]*0   # 简化为单变量情形（演示）
    beta_tt = beta_pred + (K.flatten() * (Y[i,0] - X[i] @ beta_pred))
    P = P_pred - K @ X[i].reshape(1,-1) @ P_pred
    beta_filt[i] = beta_tt

# ---------- 3. 画出"估计出的时变系数" ----------
plt.figure(figsize=(10,5))
plt.plot(beta_filt[:,0], label='filtered beta(1,1)')
plt.plot(beta_true[p:,0], 'k--', label='true beta(1,1)')
plt.legend(); plt.title('TVP-VAR filtered coefficient')
plt.savefig('tvpvar_filtered.png', dpi=120, bbox_inches='tight')
print("已保存 tvpvar_filtered.png")

# ---------- 4. 后续：在每个时点做时变 IRF ----------
# 对每个 t, 用 beta_filt[t] 构造 companion matrix,
# 迭代 h 步, 得到 IRF(t, h, shock_j)
# 比较 2008 vs 2020 的 IRF 曲线（见下一节）