简介这是一份围绕TVP-VAR时变参数向量自回归模型估计的MATLAB代码包适合经济学、金融学等领域的研究者、教师及高年级学生使用帮助解决时序数据中参数漂移、结构突变以及非线性动态传导等经典VAR模型难以处理的问题。与静态VAR相比TVP-VAR允许系数随时期平滑演变能够更细致地识别政策冲击、外部冲击在不同时点上的传导效应因而广泛应用于宏观实证与金融市场波动研究。代码以中岛上智教授2011年的经典实现为基础经二次开发后显著增强可用性增加时间标签便于将估计结果与具体年份或事件一一对应新增三维脉冲响应图能够从多个维度呈现冲击响应的动态演变路径补充sa2参数的统计信息为时变波动特征提供更完整的诊断依据。压缩包采用zip格式整体大小约2MB便于快速获取和部署目前已有92人关注学习在论文中引用时可参考Nakajima2011的标注格式。对于需要快速开展实证研究的用户这份代码可直接运行省去从零搭建与调试的精力也可作为理解模型内部机制的入门样例。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(p1:end, :); X X(p1:end, :);k里的1是截距项由于有 n 个方程β_t 的实际维度是 k×nm来自组合数 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编辑器里按CtrlShiftF全局搜索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的关系β_tT×k×nDK仿真平滑无直接关系Qk×k逆Wishart抽样无直接关系α_tT×mDK仿真平滑若sa2指S的尺度则有关系h_tT×nMH逐点抽样候选方差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拆开再拼回去避免依赖金融工具箱。不同频率的数据按下面方式处理数据频率标签格式生成方式季度1990Q1sprintf 拼年季度月度1990-01datestr(d, yyyy-mm)日度1990-01-01datestr(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以上增加迭代次数轨迹图无明显分段迁移增大nburnac 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候选方差重复使用了一次。本文还有配套的精品资源点击获取 TVP-VAR模型MATLAB复现：sa2参数定位与三维脉冲响应图绘制
简介这是一份围绕TVP-VAR时变参数向量自回归模型估计的MATLAB代码包适合经济学、金融学等领域的研究者、教师及高年级学生使用帮助解决时序数据中参数漂移、结构突变以及非线性动态传导等经典VAR模型难以处理的问题。与静态VAR相比TVP-VAR允许系数随时期平滑演变能够更细致地识别政策冲击、外部冲击在不同时点上的传导效应因而广泛应用于宏观实证与金融市场波动研究。代码以中岛上智教授2011年的经典实现为基础经二次开发后显著增强可用性增加时间标签便于将估计结果与具体年份或事件一一对应新增三维脉冲响应图能够从多个维度呈现冲击响应的动态演变路径补充sa2参数的统计信息为时变波动特征提供更完整的诊断依据。压缩包采用zip格式整体大小约2MB便于快速获取和部署目前已有92人关注学习在论文中引用时可参考Nakajima2011的标注格式。对于需要快速开展实证研究的用户这份代码可直接运行省去从零搭建与调试的精力也可作为理解模型内部机制的入门样例。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(p1:end, :); X X(p1:end, :);k里的1是截距项由于有 n 个方程β_t 的实际维度是 k×nm来自组合数 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编辑器里按CtrlShiftF全局搜索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的关系β_tT×k×nDK仿真平滑无直接关系Qk×k逆Wishart抽样无直接关系α_tT×mDK仿真平滑若sa2指S的尺度则有关系h_tT×nMH逐点抽样候选方差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拆开再拼回去避免依赖金融工具箱。不同频率的数据按下面方式处理数据频率标签格式生成方式季度1990Q1sprintf 拼年季度月度1990-01datestr(d, yyyy-mm)日度1990-01-01datestr(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以上增加迭代次数轨迹图无明显分段迁移增大nburnac 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候选方差重复使用了一次。本文还有配套的精品资源点击获取