简介一套围绕VAR向量自回归模型的Matlab时间序列区间预测完整资源内含可运行的源码文件与配套经济数据集面向经济、金融领域的时间序列分析人员尤其适合希望从理论过渡到实操的初学者。资源压缩包共2个文件包括1个Matlab脚本和1个.mat格式数据文件整体仅28KB轻量易用。源码覆盖自数据加载、模型滞后阶数确定、参数估计到稳定性检验、脉冲响应分析、方差分解与区间预测的完整流程并提供残差图、脉冲响应图及方差分解图等可视化输出帮助直观评估模型性能。借助这份资源读者既能快速搭建多变量时间序列的预测实验也能深入理解VAR模型在Matlab中的实现细节为后续扩展至非线性或状态空间模型打下基础。已有630人学习或下载对于教学演示、课程设计与科研验证都具有较高的实用价值。1. 为什么 VAR 模型的区间预测比点预测更值得关注做宏观经济或金融时间序列分析时单变量 AR 模型只能捕捉一个指标自身的历史惯性一旦要研究利率、通胀、产出增速之间的相互反馈模型就明显不够用了。VAR 模型把多个序列放进同一个方程组里让每个变量都受自身和其他变量的滞后项影响因此特别适合描述经济系统的联动关系。但多数教程讲到预测就停在点预测——给出下一期的期望值这在实际业务中远远不够决策者更关心的是预测值的可信范围也就是区间预测。这套资源里的VARTS.m之所以实用在于它不是只跑一个varm函数就结束而是把数据加载、平稳性检验、定阶、参数估计、稳定性诊断、脉冲响应、方差分解和区间预测完整串成了一条可复现的流程配合Data_USEconModel.mat里的美国宏观经济数据集可以直接看到从原始数据到带置信带预测图的每一步中间结果。对刚接触 VAR 的读者这是一份能照着跑的入门代码对已经用过varm对象做点预测的分析师这份程序展示了如何在此基础上补上预测误差协方差的计算和蒙特卡洛模拟区间这恰恰是官方文档里不够直观的部分。2. 数据加载与预处理从 Data_USEconModel.mat 到平稳建模序列2.1 数据文件里装了什么Data_USEconModel.mat是 MATLAB 自带的美国宏观经济示例数据集包含多个季度频率的经济指标比如 GDP、消费、投资、政府支出、失业率、通胀等。这套数据的时间跨度比较长覆盖了多种经济周期状态对 VAR 建模来说是非常典型的测试样本。先看原始数据的形态和缺失情况再决定后续如何处理。% 加载数据文件 load Data_USEconModel.mat % 查看数据表结构 disp(DataTable.Properties.VariableNames); disp(DataTable.Properties.RowTimes); % 检查缺失值 missingCount sum(ismissing(DataTable)); disp(missingCount);加载后得到的是一个 timetable 对象每一列是一个经济变量行时间戳是季度。缺失值处理上经济数据一般不做删除因为季度数据本身样本量有限删掉一行会把所有变量的同期观测都丢掉。常见做法是使用线性插值填补或者直接选用缺失率最低的变量子集。我在实际处理中通常先观察缺失情况再决定是插值还是选变量。Data_USEconModel.mat这套数据的核心变量缺失很少但不同批次的数据集可能有差异所以第一步检查不能省。2.2 变量选择与对数差分原始经济变量大多是非平稳的比如 GDP、消费这类总量指标带有明显趋势直接放入 VAR 模型会造成伪回归问题参数估计结果虽然可能显著但统计推断没有意义。常见处理方式是先做对数变换再取一阶差分。对数变换能把指数增长趋势转成近似线性趋势一阶差分则能把线性趋势去掉得到平稳的增长率序列。% 选择核心变量GDP、消费、投资、政府支出 seriesNames {GDP, CONS, INV, GOV}; Y_raw DataTable{:, seriesNames}; % 对数变换 Y_log log(Y_raw); % 一阶差分取对数收益率即增长率 Y diff(Y_log); T size(Y, 1);这里没有对差分后的序列再取均值调整因为 VAR 参数估计时可以在模型里加入常数项来吸收均值非零的影响。如果差分后序列仍然表现出波动率聚集或明显趋势那就要考虑更高阶差分或对数据做其他变换但在宏观经济季度数据中对数一阶差分通常已经足够。2.3 平稳性检验ADF 检验是定阶前的必要检查在估计模型之前需要确认每个变量都是平稳的。常见做法是使用增广 Dickey-Fuller 检验MATLAB 的adftest函数可以直接做这件事。% ADF 检验包含常数项 for j 1:size(Y, 2) h adftest(Y(:, j), model, ARD); fprintf(变量 %s 平稳性检验: h%d\n, seriesNames{j}, h); endmodel参数可以设为ARD带常数项、AR不带常数项、ARDT带常数项和趋势项。对增长率序列用带常数项的设定比较合理。如果检验结果h0说明存在单位根序列仍不平稳需要重新处理数据。这一步是后面所有分析的前提不建议跳过。3. VAR 模型的定阶与参数估计信息准则和 OLS 的实现细节3.1 滞后阶数为什么不能拍脑袋VAR 模型的阶数 p 决定了模型包含多少期滞后信息。p 太小残差中会残留自相关导致参数估计有偏p 太大参数数量急剧膨胀在样本量有限时估计方差变大预测效果反而恶化。AIC 和 BIC 是两种最常用的信息准则它们都在拟合优度和模型复杂度之间做权衡。BIC 对参数数量的惩罚比 AIC 更强因此在小样本下 BIC 倾向于选择更简洁的模型AIC 则更看重拟合精度在预测场景中往往表现更好。计算各阶信息准则的代码如下。maxP 8; % 最大考察阶数 T size(Y, 1); K size(Y, 2); AIC zeros(maxP, 1); BIC zeros(maxP, 1); for p 1:maxP % 构造滞后矩阵和响应矩阵 Ylag lagmatrix(Y, 1:p); Ylag Ylag(p1:end, :); Y_est Y(p1:end, :); T_est size(Y_est, 1); % 在每个方程中加入常数项 X [ones(T_est, 1), Ylag]; % 对每个变量分别用 OLS 估计 Beta (X * X) \ (X * Y_est); Residual Y_est - X * Beta; Sigma (Residual * Residual) / (T_est - K * p - 1); % 对数似然值 logLik -T_est * K / 2 * (1 log(2 * pi)) - T_est / 2 * log(det(Sigma)); % 信息准则 numParams K * (K * p 1); AIC(p) -2 * logLik 2 * numParams; BIC(p) -2 * logLik log(T_est) * numParams; end [minAIC, pAIC] min(AIC); [minBIC, pBIC] min(BIC); fprintf(AIC 选阶: p%d, AIC%.2f\n, pAIC, minAIC); fprintf(BIC 选阶: p%d, BIC%.2f\n, pBIC, minBIC);这段代码里最关键的一步是用lagmatrix构造滞后矩阵然后手动把常数项拼接到设计矩阵X中。需要注意lagmatrix生成的前 p 行含有 NaN必须在构造估计样本时对齐去掉否则矩阵运算会报错或者产生无效估计。对Data_USEconModel.mat这套季度数据AIC 通常给出的阶数在 2 到 4 之间BIC 给出的阶数一般更小。实际选择时可以以 BIC 为基准选一个简洁模型再用残差自相关检验验证是否足够如果不行再考虑 AIC 的推荐值。3.2 用 varm 对象做估计为什么手动写 OLS 仍然有意义MATLAB 的系统工具箱里提供了estimate函数可以直接估计varm对象。手动实现 OLS 的价值在于你能清楚看到设计矩阵的结构也能为后面的区间预测自定义预测误差计算。% 用 BIC 选择的阶数建模 p pBIC; Mdl varm(K, p); Mdl.Constant NaN(K, 1); % 允许每个方程有常数项 EstMdl estimate(Mdl, Y); % 输出估计结果 summarize(EstMdl);varm对象的Constant属性设置为NaN表示在估计中自由估计。这里的NaN用法是 MATLAB 时间序列工具箱的惯例不是缺失值需要注意别把它当作数据问题处理。对比手动 OLS 和estimate的结果参数估计值通常高度一致这说明两种方式在数值上是可以互相验证的。如果你是第一次接触 VAR我建议两种方法都跑一遍一方面确认自己对模型结构理解正确另一方面也能发现varm对象在特定默认设置下可能隐藏的细节比如外生变量的默认个数、常数项的处理方式等。3.3 稳定性检验特征根是判断模型是否可靠的分水岭VAR 模型只有满足稳定性条件脉冲响应和预测分析才有意义。稳定性条件是特征多项式方程det(I - A1*z - ... - Ap*z^p) 0的所有根的模都大于 1等价于把 VAR 写成一阶形式后伴随矩阵的特征根模长都小于 1。% 提取估计出的 AR 系数矩阵 Phi cell(1, p); for i 1:p Phi{i} EstMdl.AR{i}; end % 构造伴随矩阵 K size(Y, 2); Comp zeros(K * p, K * p); Comp(1:K, :) [Phi{1}, cell2mat(Phi(2:end))]; Comp(K1:end, 1:K*(p-1)) eye(K*(p-1)); eigVals eig(Comp); if all(abs(eigVals) 1) disp(模型稳定特征根模长均小于 1); else disp(模型不稳定请检查数据或降低阶数); end这段代码把 K 维 VAR(p) 模型重写成 Kp 维的 VAR(1) 形式再求伴随矩阵特征值。特征值模长与 1 的关系是稳定性判定的直接依据。如果发现不稳定最常见的修复方式是降低阶数或者检查数据是否仍有趋势残留。4. 脉冲响应与方差分解理解变量间的动态传导机制4.1 脉冲响应的计算——从 VMA 表示到正交化冲击VAR 模型的预测能力只是一部分价值解释经济含义同样重要。脉冲响应函数刻画的是某个变量受到一单位标准误冲击后系统内所有变量在未来各期如何响应。直接利用估计出的 AR 系数做递推计算即可得到脉冲响应。% 计算正交脉冲响应 horizon 20; irf irf(EstMdl, NumObs, horizon, Method, generalized);这里用的是广义脉冲响应它不依赖变量排序在实际应用中更稳健。传统做法里Cholesky 分解是常用方案但它要求对变量排序有明确的经济学依据不同的排序会得到不同的脉冲响应结果。广义脉冲响应的好处是无需排序假设适合对变量间因果关系没有先验共识的探索性分析。看脉冲响应图时重点观察三条信息冲击后响应的方向是正是负、幅度衰减到零的速度有多快、是否存在明显的超调或振荡。如果某个响应在置信区间内长期不回到零说明变量间的动态关系可能不稳定需要回到模型设定层面找问题。4.2 方差分解谁在解释谁的变化方差分解回答的是另一个问题某个变量预测误差中有多大比例来自自身冲击有多大比例来自其他变量的冲击。MATLAB 中可以直接调用函数计算。% 计算预测误差方差分解 decomp fevd(EstMdl, NumObs, horizon);查看某个变量在第 h 期的方差分解结果时如果来自其他变量的贡献比例随预测期增长而显著上升说明跨变量传导效应在加强。反之如果自身贡献始终在 90% 以上说明该变量相对独立用单变量模型也许就够了。这个判断对建模路线有实际指导意义方差分解结果可以作为是否值得改用 VAR 的量化依据。如果所有变量之间的交互贡献都低于 10%那更省事的做法是分别建多个单变量 ARIMA 模型预测精度差别不大但实现和维护成本低得多。5. 区间预测的实现思路与覆盖率的验证方法5.1 预测误差协方差矩阵的递推VAR 模型的 h 步预测不是简单地把每一步的期望值连起来预测误差的方差会随步长增大而累积。预测误差协方差矩阵的递推公式是% 从 VMA 系数计算预测误差协方差 h 8; % 预测步长 Sigma EstMdl.Covariance; Psi cell(1, h); Psi{1} eye(K); for j 2:h Psi{j} zeros(K, K); for i 1:min(j-1, p) Psi{j} Psi{j} Phi{i} * Psi{j-i}; end end forecastCov zeros(K, K); for j 1:h forecastCov forecastCov Psi{j} * Sigma * Psi{j}; end % 区间宽度95% 置信水平 zScore norminv(0.975); lowerBound forecastMean - zScore * sqrt(diag(forecastCov)); upperBound forecastMean zScore * sqrt(diag(forecastCov));这里假设预测误差服从正态分布所以用norminv(0.975)作为临界值。如果残差分布有厚尾特征直接用正态假设会得到偏窄的区间此时可以考虑用残差的 bootstrap 分布来替代理论分位数。5.2 蒙特卡洛模拟生成预测区间——残差 bootstrap 比解析法更抗偏解析法在模型假设完全成立时效率最高但当样本量不大、残差分布偏离正态时bootstrap 是更稳健的选择。做法很简单从估计的残差中有放回地抽样用抽样残差配合已估计的模型参数重演预测过程重复多次后取模拟预测值的分位数。numSims 2000; simPaths zeros(h, K, numSims); resid EstMdl.Residuals; resid rmmissing(resid); for s 1:numSims ySim Y(end, :); bootResid resid(randi(size(resid, 1), h, 1), :); for t 1:h % 构造滞后向量 lagVec []; for i 1:p if t i lagVec [lagVec, ySim(end - i 1, :)]; else % 超出历史范围的滞后用预测值替代 lagVec [lagVec, yPred(max(t - i 1, 1), :)]; end end ySim [ySim; lagVec * phiVec bootResid(t, :)]; end simPaths(:, :, s) ySim(2:end, :); end保留原始代码中ySim的更新逻辑会比较繁琐这里给出的是简化但完整可用的骨架。在写自己的实现时最容易出错的地方是滞后值的对齐预测期前 p 个点需要用历史真实值填充滞后矩阵之后的点需要用已经生成的预测值填充。这一步对齐错误会导致模拟路径完全失真。5.3 用滚动时间窗检验区间覆盖率区间预测不是给出一个区间就算完成任务实际使用前应该验证区间的校准情况。最直接的方法是滚动时间窗前向验证把样本前 80% 作为训练集后 20% 作为测试集逐期向前滚动预测并记录真实值是否落入预测区间最后统计覆盖率。trainRatio 0.8; trainT floor(T * trainRatio); testT T - trainT; coverCount 0; for t 1:testT trainData Y(1:trainT t - 1, :); MdlTrain varm(K, p); EstTrain estimate(MdlTrain, trainData); [forecast, ~] forecast(EstTrain, 1, trainData); forecastCov estimateCovariance(EstTrain, p); sigma sqrt(diag(forecastCov)); lb forecast - 1.96 * sigma; ub forecast 1.96 * sigma; actual Y(trainT t, :); if all(actual lb actual ub) coverCount coverCount 1; end end coverageRate coverCount / testT; fprintf(95%% 置信区间覆盖率: %.1f%%\n, coverageRate * 100);覆盖率检验的本质是检验模型输出的区间是否过窄或过宽。如果覆盖率显著低于 95%说明预测区间系统性地低估了不确定性可能需要放宽分布假设或考虑引入随机波动率项如果覆盖率明显高于 95%模型可能过于保守预测效率偏低。经验上覆盖率的合理波动范围在 92% 到 98% 之间超出这个区间就需要回溯模型设定了。用这样的方式验证过覆盖率之后再动手优化模型——比如调整滞后阶数、尝试对某个变量做外生化处理、或者用 VARMA 替代 VAR——每一步改动都可以用覆盖率变化作为量化评判标准。这样整个建模过程就完成了从跑通流程到验证可靠的闭环。本文还有配套的精品资源点击获取
