简介：这是一份围绕TVP-VAR（时变参数向量自回归）模型估计的MATLAB代码包，适合经济学、金融学等领域的研究者、教师及高年级学生使用，帮助解决时序数据中参数漂移、结构突变以及非线性动态传导等经典VAR模型难以处理的问题。与静态VAR相比，TVP-VAR允许系数随时期平滑演变，能够更细致地识别政策冲击、外部冲击在不同时点上的传导效应，因而广泛应用于宏观实证与金融市场波动研究。代码以中岛上智教授2011年的经典实现为基础，经二次开发后显著增强可用性：增加时间标签，便于将估计结果与具体年份或事件一一对应；新增三维脉冲响应图，能够从多个维度呈现冲击响应的动态演变路径；补充sa2参数的统计信息，为时变波动特征提供更完整的诊断依据。压缩包采用zip格式，整体大小约2MB，便于快速获取和部署；目前已有92人关注学习，在论文中引用时可参考Nakajima（2011）的标注格式。对于需要快速开展实证研究的用户，这份代码可直接运行，省去从零搭建与调试的精力，也可作为理解模型内部机制的入门样例。
1. TVP-VAR模型的MATLAB复现，先别急着跑主程序
拿到一份TVP-VAR模型的MATLAB代码，先花十分钟把sa2这个变量找出来，比先跑回归更重要。TVP-VAR把VAR的系数、同期矩阵和波动率都放宽为随时间演化的状态量，估计依赖MCMC，代码里到处是先验方差和候选方差，sa2正是最容易被当作常数误写的那一处。
sa2设定错了，三维脉冲响应图照样画得出来，但曲线会安静地失真。这套代码要做三件事：跑通MCMC主循环、把脉冲响应按时间维度展成三维图、在图上给出可读的季度时间标签。
适合做宏观或金融时间序列研究、需要报告时变脉冲响应结果的人。下面的做法以Nakajima的经典设定为蓝本，只用到MATLAB基础函数，不依赖深度学习之类额外工具箱，按常规的matlab安装流程装好环境就能跑。
2. TVP-VAR模型的MCMC估计主框架：先验、状态方程与sa2的位置
2.1 从状态空间表达式到MATLAB的变量布局
TVP-VAR的状态空间形式写成y_t = X_t β_t + A_t^{-1} Σ_t ε_t，其中 β_t 是含截距和滞后项的共同系数向量，A_t 是下三角化的同期矩阵，Σ_t = diag(exp(h_t/2)) 由随机波动率驱动。MCMC要抽的后验对象就是 {β_t, α_t, h_t} 三条路径，以及它们各自状态方程里的超参数。动手写码之前，先把维度算清楚：
[T, n] = size(Y); % 样本期数、变量个数 p = 2; % 滞后阶数 k = n * p + 1; % 每个方程右侧变量个数（含截距） m = n * (n - 1) / 2; % 下三角矩阵 A_t 的自由参数个数 X = lagmatrix(Y, 1:p); % 全部方程共用的回归元 Y = Y(p+1:end, :); X = X(p+1:end, :);k里的+1是截距项，由于有 n 个方程，β_t 的实际维度是 k×n；m来自组合数 C(n,2)，是 A_t 中非对角线元素的个数。lagmatrix生成的前 p 行是 NaN，所以数据窗口要对齐。
需要注意 X 是全方程共用的，这和“每个方程右侧变量各自不同”的贝叶斯VAR写法有区别。混用会直接导致维度报错，而且这种报错常被误判成数据问题，实际是模型设定问题。
2.2 sa2在估计框架里的两个常见位置
在流传的TVP-VAR代码版本里，sa2这个变量名通常出现在两处。第一处是随机波动率状态方程h_t = h_{t-1} + η_t的创新方差，固定为常数，不参与抽样；第二处是同期参数块先验逆Wishart分布的尺度参数。两种位置的后验含义完全不同，拿到代码第一步是定位它。
% sa2 位置一：SV状态方程创新方差（固定常数） sa2 = 0.02 * ones(n, 1); % sa2 位置二：alpha块先验（逆Wishart）的尺度参数 S_scale = sa2 * eye(m); S_df = m + 1;2.2.1 用全局搜索定位sa2的真实口径
在MATLAB编辑器里按Ctrl+Shift+F全局搜索sa2。如果它出现在生成候选随机数的语句里，比如h_cand = h_old + sqrt(sa2) * randn，那它是MH步的调节方差；如果它出现在iwishrnd或invwishrnd调用的参数位置，那它是先验尺度。搜索时配合编辑器自带的代码补全逐个点开引用位置，比肉眼扫更快，也能顺带看清它有没有被第二个函数改写过。两处一旦混用，后验会悄悄收缩或发散，三维脉冲响应图的形态会变得异常平滑或异常尖锐。
2.3 Gibbs主循环的骨架与后验结果保存
主循环按固定顺序更新六个块：β_t、Q、α_t、S、h_t。下面是最小可运行骨架，四个子函数按Nakajima经典设定实现：
nsave = 2000; nburn = 500; beta_save = zeros(nsave, T, k*n); alpha_save = zeros(nsave, T, m); h_save = zeros(nsave, T, n); for iter = 1:(nburn + nsave) beta = sim_smoother_beta(Y, X, beta, Q, h); Q = iwishrnd(inv(SSR_beta + inv(prior_Q)), df_Q + T); alpha = sim_smoother_alpha(Y, X, beta, alpha, S, h); S = iwishrnd(inv(SSR_alpha + inv(prior_S)), df_S + T); h = mh_sampler_h(Y, alpha, Q, sa2, h); if iter > nburn idx = iter - nburn; beta_save(idx, :, :) = beta; alpha_save(idx, :, :) = alpha; h_save(idx, :, :) = h; end end参数说明：nburn是预烧期，这段抽样不保存，用于让马尔可夫链进入平稳分布；nsave是正式抽样次数，后面所有后验统计量都只用这一段。iwishrnd的第一个参数要传协方差矩阵的逆，第二个是自由度，MATLAB的记号习惯和教科书公式差在这里，最容易写错。mh_sampler_h内部对每个时间点逐点做Metropolis-Hastings更新，候选方差就是sa2，它决定整条h链的移动效率，也是后文调参的主要对象。
| 变量 | 维度 | 更新方式 | 与sa2的关系 |
|---|---|---|---|
| β_t | T×k×n | DK仿真平滑 | 无直接关系 |
| Q | k×k | 逆Wishart抽样 | 无直接关系 |
| α_t | T×m | DK仿真平滑 | 若sa2指S的尺度则有关系 |
| h_t | T×n | MH逐点抽样 | 候选方差=sa2 |
| σ_h² | n×1 | 固定常数或逆Gamma | 若固定，即sa2本身 |
提示：如果下载的代码里sim_smoother_beta和mh_sampler_h写在同一个脚本中，先确认主程序调用顺序与这里一致。顺序颠倒不会报错，但后验会错得非常隐蔽。
3. 三维脉冲响应图的MATLAB实现：数据排成T×H后再贴时间标签
3.1 从MCMC抽样里计算出每个时间点的脉冲响应
三维图的三根轴分别是冲击后的响应期数、样本时间、响应值。因此画图前必须先把后验抽样整理成“每一期、每一响应期数”一个数。常见做法是保留每次MCMC抽样的结果，再对某一维度取均值或分位数：
H = 16; % 冲击后看16期 irf_draw = zeros(nsave, T, H, n); % 抽样×时间×响应期×变量 for ii = 1:nsave for tt = 1:T Bt = reshape(beta_save(ii, tt, :), k, n); At = build_A(squeeze(alpha_save(ii, tt, :)), n); St = diag(exp(squeeze(h_save(ii, tt, :)))); irf_draw(ii, tt, :, :) = calc_irf(Bt, At, St, H); end end irf_mean = squeeze(mean(irf_draw, 1)); % T×H×nbuild_A把堆成m维的α_t还原成下三角矩阵，calc_irf做完Cholesky分解后按VAR(1)递归累积。这里的关键是“累积”二字：如果不累积，三维图看到的是逐期增量，和经济含义中“某一期冲击对之后各期的累计影响”不是一回事。贴报告前先确认这个函数最后一步做了cumsum。
提示：当 T 很大时，irf_draw会占掉几百MB内存，可以只保留16%和84%分位抽样，或把nsave拆成若干块循环累加。
3.2 用surf画三维脉冲响应图：网格、着色与视角
MATLAB里最顺手的是surf，数据矩阵行是时间、列是响应期数：
i_var = 2; % 画第二个变量的响应 figure('Color', 'w', 'Position', [80 80 780 440]); hs = surf(1:H, 1:T, irf_mean(:, :, i_var)); hs.EdgeAlpha = 0.35; % 网格线半透明 shading interp; colormap(parula); colorbar; xlabel('冲击后的季度数'); ylabel('样本时间（截面）'); zlabel('累积脉冲响应'); view(38, 42);hs.EdgeAlpha设成0.35让网格线半透明，颜色带全部用shading interp平滑。很多三维图画出来黑糊糊一团，就是默认edge颜色权重太重。view(38,42)是常用视角，改成view(0,90)相当于俯视图，会直接退回二维热力图形态，适合快速自查数值范围。
3.3 时间标签怎么贴：从截面序号映射到季度字符串
surf的纵轴默认是行索引1:T，需要替换成真实日期。季度数据先造好标签数组，再通过set(gca,'YTickLabel',...)覆盖：
date_labels = strings(T, 1); for i = 1:T d = dates(i + p); % 对齐去掉的 p 期 date_labels(i) = sprintf('%04dQ%d', year(d), ceil(month(d) / 3)); end step = max(1, floor(T / 8)); % 纵向最多放 8 个刻度 yt = 1:step:T; set(gca, 'YTick', yt, 'YTickLabel', date_labels(yt));年份和季度用year、month拆开再拼回去，避免依赖金融工具箱。不同频率的数据按下面方式处理：
| 数据频率 | 标签格式 | 生成方式 |
|---|---|---|
| 季度 | 1990Q1 | sprintf 拼年+季度 |
| 月度 | 1990-01 | datestr(d, 'yyyy-mm') |
| 日度 | 1990-01-01 | datestr(d, 'yyyy-mm-dd') |
时间标签决定你能不能直接回答“样本后期脉冲响应有没有变化”，这也是“增加时间标签”要解决的最后一个环节。贴完标签把图旋转一圈，确认纵轴刻度没有重叠，再去考虑画多个变量的对比子图。
4. sa2的取值逻辑与收敛检查：先看接受率再谈调参
4.1 sa2充当SV创新方差时的平滑度判断
如果sa2停在h的状态方程里，它控制波动率路径的平滑程度。sa2越小，h_t变化越慢，后验越接近恒定波动率；sa2越大，h_t会追踪残差的短促波动，也容易把估计拉向极端。经验窗口是0.01到0.05，样本越长、波动越剧烈，越取中上端。
判定办法是输出h的后验均值路径，看相邻期差分的量级：
h_mean = squeeze(mean(h_save, 1)); % T×n dh = abs(diff(h_mean, 1, 1)); % 相邻期变化 fprintf('max %.4f, median %.4f\n', max(dh(:)), median(dh(:)));如果h路径像锯齿一样在相邻期之间大幅往复，说明sa2偏大，按量级缩小后重跑；如果h几乎是一条平线，说明sa2偏小，波动率部分没被识别出来，同样要重跑。这个判断比单纯看后验均值图更量化。
4.2 sa2充当MH候选方差时，用接受率反馈调整
MH候选方差过大，候选点被持续拒绝，h链停在原地形成阶梯状轨迹；方差过小，接受率接近1但遍历速度极慢。标准做法是让接受率落在25%到50%区间，自适应调整放在预烧期：
acc_rate = acc_count / window_len; if acc_rate < 0.25 sa2 = sa2 * 0.8; elseif acc_rate > 0.50 sa2 = sa2 * 1.2; endsa2的初值决定第一次尝试能否被接受，初值不要拍脑袋填0.1，用之前已跑通数据集的结果做起点更稳。调整动作集中在nburn阶段完成，正式抽样期固定sa2不再改动，否则后验样本不满足同分布假设。
4.3 收敛检查：自相关、有效样本量与轨迹图
调完sa2之后输出三类诊断，任何一类不过都回到前两节修改参数：
| 诊断项 | 通过标准 | 失败时的处理方向 |
|---|---|---|
| h_t自相关 | 10阶滞后内降到0.1以下 | 增nsave或调MH候选方差 |
| 有效样本量 | 核心参数在500以上 | 增加迭代次数 |
| 轨迹图 | 无明显分段迁移 | 增大nburn |
ac = autocorr(h_mean(:, 1), 20); bar(0:20, ac); xlabel('滞后'); ylabel('自相关');h的自相关收敛慢，优先回去调sa2；β_t的自相关大，多半是Q先验太紧，和sa2没有关系。把这两个问题的排查方向分开，能省掉大量试错时间。
5. 三维脉冲响应图的验证技巧：时间截面切片与sa2敏感性对拍
三维图画完不要直接贴进报告。第一步，抽一个时间截面，把三维数据的这一行拉出来画普通二维脉冲响应，并叠上后验分位带：
t_c = round(T * 0.7); % 样本后期某截面 irf_q = squeeze(quantile(irf_draw(:, t_c, :, i_var), [0.16 0.50 0.84], 1)); figure; plot(1:H, irf_q(2, :), 'b-', 'LineWidth', 1.5); hold on; fill([1:H fliplr(1:H)], [irf_q(1, :) fliplr(irf_q(3, :))], ... [0.8 0.8 0.8], 'FaceAlpha', 0.4); xlabel('响应期数'); ylabel('脉冲响应'); grid on;16%和84%分位构成68%后验带。对照你熟悉的经济直觉，确认响应方向、峰值出现的期数是否合理，再贴三维图。这一步能同时暴露两类问题：数据排列错位，以及结构识别矩阵方向错误。
第二步是sa2敏感性检查。把sa2初值分别设为0.005、0.02、0.05三档重跑MCMC，比较同一时间截面的峰值响应期是否移动超过一期。移动小于一期说明结论稳健；移动大说明三维图上看到的“时变”可能是SV采样噪声，而不是真实的经济结构变化。三维图与二维图对拍时若量级差两倍以上，先排查sa2是否在同一份代码里既被当作SV创新方差、又被当作MH候选方差重复使用了一次。
本文还有配套的精品资源，点击获取
