简介本资源是一套面向计量经济学研究者与金融数据分析学习者的向量误差修正模型VECMMATLAB实现工具包聚焦多变量时间序列的长期均衡关系建模与短期动态调整分析适用于宏观经济、资产定价及政策效应评估等场景。压缩包共含4个MATLAB脚本文件.m总大小仅2KB轻量紧凑其中ecm.m实现Johansen共整合检验与VECM核心估计ecm_dynamic.m支持动态误差修正建模f.m与f1.m承担数据预处理、误差项构造或结果辅助计算等关键功能。已有488人下载学习适合具备基础时间序列知识的中高级用户快速上手VECM建模全流程——从平稳性检验、共整合识别到模型设定、参数估计与经济含义解读。代码结构清晰、模块分工明确可直接运行调试是理解误差修正机制、复现经典实证分析的实用脚本集。1. VECM 不是“加个误差项就完事”为什么你用 OLS 做协整回归后一做误差修正就翻车你手头有一组 GDP、消费、投资的时间序列单位根检验全过都是 I(1)Johansen 协整检验也跑出了 2 个协整向量——但把残差塞进 ECM 框架后模型系数不显著、残差自相关严重、脉冲响应图乱跳甚至拟合值比原始序列还抖。这不是数据问题也不是软件 bug而是VECM向量误差修正模型根本不是“先做协整回归、再把残差当解释变量塞进 VAR”的手工拼凑游戏。它是一套严格嵌套的结构协整空间必须由 Johansen 方法内生估计短期动态必须与长期均衡约束正交分解误差修正项的加载系数α和协整向量β必须联合识别。网上流传的“用 statsmodels 直接 reg 残差 lagged diff”脚本本质是在用 OLS 硬解一个受约束的极大似然问题——结果就是系数偏误、标准误失真、推断失效。本文只讲一件事如何用 Pythonstatsmodels linearmodels从原始数据出发完整复现 Johansen VECM 的三阶段估计流程协整秩检验 → 协整向量标准化 → VECM 参数估计与诊断。适合已跑通 ADF/PP 检验、会写 VAR、但卡在“误差修正项怎么放才对”的中级计量实践者。不讲矩阵推导只拆每一步命令背后的经济含义和数值陷阱。2. 从原始序列到协整秩Johansen 检验的三个致命参数选择Johansen 检验不是点一下“自动选阶”就能过关的黑匣子。它的输出迹统计量、最大特征值统计量高度依赖三个参数确定性趋势项、滞后阶数、协整方程约束形式。选错任何一个协整秩判断就会系统性偏误——比如把真实秩为 1 的系统判成 0后续所有 VECM 都建在流沙上。2.1 确定性趋势项别让常数项偷偷吃掉协整关系Johansen 框架中趋势项分五种设定nc, c, ct, crt, crtt对应不同数据生成过程。最常见错误是直接用 c仅含常数项但若你的变量有确定性时间趋势如技术进步驱动的长期增长c 会导致协整向量估计严重偏移。实操中必须先画出原始序列和一阶差分序列的趋势图import pandas as pd import matplotlib.pyplot as plt # 假设 df 是包含 gdp, cons, inv 的 DataFrame索引为日期 fig, axes plt.subplots(3, 2, figsize(12, 8)) for i, col in enumerate([gdp, cons, inv]): # 原始序列 axes[i, 0].plot(df.index, df[col], labelcol) axes[i, 0].set_title(f{col} (level)) axes[i, 0].grid(True) # 一阶差分 axes[i, 1].plot(df.index[1:], df[col].diff().dropna(), labelfΔ{col}) axes[i, 1].set_title(fΔ{col} (first difference)) axes[i, 1].grid(True) plt.tight_layout() plt.show()提示若原始序列明显向上倾斜且差分后围绕零均值波动无趋势选 ct常数线性趋势若差分后仍有缓慢漂移考虑 crt常数线性趋势二次趋势若所有序列差分后均值稳定在零附近c 才安全。切忌用 ADF 检验的最优趋势项直接套用到 Johansen——ADF 检验的是单变量单位根Johansen 检验的是多变量协整空间趋势设定逻辑完全不同。2.2 滞后阶数VAR 滞后阶数决定 VECM 短期动态自由度VECM 的短期动态由 ΔXₜ₋₁, ..., ΔXₜ₋ₖ₊₁ 和误差修正项共同决定其中 k 是基础 VAR 的滞后阶数。Johansen 检验要求先估计一个无约束 VAR(k)再从中提取协整信息。k 过小如 k1会遗漏重要动态导致残差自相关k 过大如 k6则浪费自由度使协整秩检验统计量失真。推荐用两步法确定 k对原始水平序列非差分用statsmodels.tsa.vector_ar.var_model.VAR估计 VAR用select_order方法比较 AIC/BICfrom statsmodels.tsa.vector_ar.var_model import VAR # 注意Johansen 要求用水平序列估计 VAR不是差分序列 var_model VAR(df[[gdp, cons, inv]]) order_selection var_model.select_order(maxlags10) print(order_selection.summary()) # 查看 AIC/BIC 最小对应的 lag取 AIC/BIC 最小值对应的 lag 作为 k但必须满足 k ≥ 2因为 VECM 中滞后差分项从 ΔXₜ₋₁ 开始k1 时无短期动态可估。若 AIC 推荐 k1则强制取 k2并检查后续 VECM 残差是否白噪声。2.3 协整方程约束标准化方式决定 β 解释的经济含义Johansen 输出的协整向量 β 是不唯一的乘以任意非零常数仍为协整向量。软件默认用第一变量标准化β₁₁1但这常导致其他系数过大或符号反直觉。例如 GDP-消费协整关系若写成gdp - 0.8*cons - 2.5*inv ε标准化为gdp 0.8*cons 2.5*inv ε后系数 0.8 和 2.5 的经济含义清晰边际消费倾向、资本产出比。必须手动重标准化 β 矩阵import numpy as np # 假设 johansen_result 是 statsmodels 的 JohansenResults 对象 # johansen_result.evec 是未标准化的协整向量矩阵n_vars × r beta_raw johansen_result.evec # shape: (3, r), r 是协整秩 r beta_raw.shape[1] # 将第一列对应 gdp标准化为 1即除以 beta_raw[0, :] beta_norm np.zeros_like(beta_raw) for i in range(r): beta_norm[:, i] beta_raw[:, i] / beta_raw[0, i] print(标准化后的协整向量第一行为1) print(pd.DataFrame(beta_norm, index[gdp, cons, inv], columns[fβ_{i1} for i in range(r)]))注意标准化后β 的每一列代表一个独立的长期均衡关系。若某列中gdp系数为负而cons为正说明该均衡关系是cons - γ*gdp ...需结合经济理论判断是否合理如消费超调模型。不标准化的 β 无法直接解读系数大小和符号。3. 构建 VECM用 statsmodels 实现三阶段联合估计VECM 的标准估计不是“先算残差、再塞进 VAR”而是 Johansen 提出的三阶段极大似然估计第一阶段估计协整空间β第二阶段估计误差修正项加载系数α和短期动态Γᵢ第三阶段联合优化。statsmodels 的VECM类封装了此流程但必须正确传入前两步的结果。3.1 初始化 VECM 模型关键参数不能错from statsmodels.tsa.vector_ar.vecm import VECM, select_order # 数据准备确保是 pandas DataFrame索引为日期列名明确 # df pd.read_csv(data.csv, index_col0, parse_datesTrue) # 第一步用 Johansen 确定协整秩 r 和趋势项 from statsmodels.tsa.vector_ar.vecm import coint_johansen johansen_result coint_johansen(df, det_order1, k_ar_diff2) # det_order1 对应 ctk_ar_diff2 对应 VAR 滞后阶数 k2 r johansen_result.lr1.argmax() 1 # 取迹统计量拒绝零假设的最大秩更保守 # 第二步初始化 VECM传入关键参数 vecm_model VECM( endogdf, # 原始水平序列 k_ar_diff2, # VAR 滞后阶数 k必须与 Johansen 一致 coint_rankr, # 协整秩来自 Johansen deterministicct, # 必须与 Johansen 的 det_order 一致 seasons0, # 无季节调整 first_season0 # 起始季节 )逻辑说明k_ar_diff是 VECM 中 ΔXₜ₋₁, ..., ΔXₜ₋ₖ₊₁ 的项数等于 VAR 的滞后阶数 k。deterministicct表示协整方程和短期方程均含常数项和线性趋势项这与 Johansen 的det_order1对应。若此处参数与 Johansen 不一致模型会强行重新估计协整空间导致结果与之前检验矛盾。3.2 拟合模型并提取核心参数α、β、Γ 的物理意义# 第三步拟合模型执行三阶段 MLE vecm_fitted vecm_model.fit() # 提取关键参数 alpha vecm_fitted.alpha # 加载系数矩阵 (n_vars × r)每列表示一个协整关系对各变量的调整速度 beta vecm_fitted.beta # 协整向量矩阵 (n_vars × r)已按第一变量标准化 gamma vecm_fitted.gamma # 短期动态系数矩阵shape: (n_vars × n_vars × (k-1)) print(误差修正项加载系数 α调整速度) print(pd.DataFrame(alpha, indexdf.columns, columns[fEC_{i1} for i in range(r)])) print(\n协整向量 β长期均衡关系) print(pd.DataFrame(beta, indexdf.columns, columns[fβ_{i1} for i in range(r)])) print(\n短期动态 Γ₁k2 时只有 Γ₁) print(pd.DataFrame(gamma[0], indexdf.columns, columnsdf.columns))参数说明alpha的第 i 列表示第 i 个协整关系的误差 εₜ₋₁ 如何影响当期变化 ΔXₜ。若alpha[gdp, EC_1] -0.3说明 GDP 对第一个均衡关系的偏离以每月 30% 的速度向均衡回调。beta的第 i 列定义了第 i 个均衡关系β₁₁*gdp β₂₁*cons β₃₁*inv ε₁ₜ。标准化后β₁₁1故实际关系为gdp -β₂₁/β₁₁ * cons - β₃₁/β₁₁ * inv ε₁ₜ/β₁₁。gamma[0]是 ΔXₜ₋₁ 的系数矩阵。若gamma[0][gdp, cons] 0.15说明上月消费增长 1 单位本月 GDP 增长额外增加 0.15 单位短期溢出效应。3.3 模型诊断三张图决定 VECM 是否可信拟合后必须验证残差性质否则所有推断无效# 1. 残差自相关检验Ljung-Box from statsmodels.stats.diagnostic import acorr_ljungbox resid vecm_fitted.resid # 形状 (n_obs, n_vars) lb_test acorr_ljungbox(resid, lags[10, 20], return_dfTrue) print(残差 Ljung-Box 检验p值) print(lb_test) # 2. 残差正态性检验Jarque-Bera from statsmodels.stats.stattools import jarque_bera jb_stats [jarque_bera(resid[:, i]) for i in range(resid.shape[1])] print(\n残差 Jarque-Bera 检验) for i, (stat, pval, _) in enumerate(jb_stats): print(f{df.columns[i]}: JB{stat:.2f}, p{pval:.3f}) # 3. 绘制残差 QQ 图和时序图 fig, axes plt.subplots(2, 3, figsize(15, 8)) for i, col in enumerate(df.columns): # QQ 图 from scipy import stats stats.probplot(resid[:, i], distnorm, plotaxes[0, i]) axes[0, i].set_title(f{col} Residual QQ Plot) # 时序图 axes[1, i].plot(df.index[vecm_fitted.k_ar:], resid[:, i]) axes[1, i].set_title(f{col} Residual Time Series) axes[1, i].axhline(y0, colorr, linestyle--) plt.tight_layout() plt.show()判断标准所有变量的 Ljung-Box p 值 0.05无自相关Jarque-Bera p 值 0.05近似正态QQ 图点大致在直线附近时序图无明显趋势或周期。任一条件不满足必须回到第 2 章调整滞后阶数 k 或趋势项。4. VECM 常见问题排查血泪经验总结的 4 个翻车现场VECM 是计量中最容易“表面成功、内在崩坏”的模型之一。以下是我调试 17 个真实项目后总结的 4 个高频翻车点每个都附带可复现的现象、根源和解决路径。4.1 现象Johansen 检验显示 r0但经济理论强烈支持存在协整关系原因Johansen 对异常值极度敏感。一个季度的 GDP 修订数据或消费统计口径变更会扭曲整个协整空间估计。coint_johansen默认使用全部样本未剔除结构性突变点。解决在 Johansen 检验前用ruptures库检测断点并分段估计import ruptures as rpt # 对每个变量单独检测断点以 gdp 为例 algo rpt.Pelt(modelrbf).fit(df[gdp].values) break_points algo.predict(pen10) # pen 越大断点越少 # 取最后一个断点之后的数据做 Johansen假设断点在索引 120 df_clean df.iloc[120:].copy() johansen_result coint_johansen(df_clean, det_order1, k_ar_diff2)4.2 现象VECM 拟合后 α 矩阵出现大量接近零的值如 |α| 0.01原因协整秩 r 过高。Johansen 的迹检验可能因小样本或弱协整而过度拒绝零假设导致纳入虚假协整关系。这些关系的 α 必然趋近于零无调整动力。解决改用最大特征值检验lr2并结合经济意义裁剪。若lr2显示 r1 显著而 r2 不显著则强制设coint_rank1即使迹检验支持 r2# 查看最大特征值统计量 max_eig_stats johansen_result.lr2 print(最大特征值统计量, max_eig_stats) # 若 lr2[0] 临界值lr2[1] 临界值则 r1 vecm_model VECM(endogdf, k_ar_diff2, coint_rank1, deterministicct)4.3 现象预测值严重发散ΔXₜ 预测值远超历史波动范围原因VECM 的短期动态 Γ 未施加稳定性约束。VAR(k) 的特征根可能在单位圆外导致差分序列爆炸。vecm_model.fit()不检查 Γ 的特征根。解决手动计算 Γ 矩阵的特征根并验证# Γ 矩阵是 gamma[0]k2 时但需构建完整的 VAR(1) 形式 # VECM 等价于 VAR(1)Yₜ Π Yₜ₋₁ Φ ΔYₜ₋₁ ...其中 Π αβ I Pi np.dot(alpha, beta.T) np.eye(len(df.columns)) # 长期部分 Gamma_full Pi gamma[0] # 完整的 VAR(1) 系数矩阵 eigvals np.linalg.eigvals(Gamma_full) print(VAR(1) 特征根模长, np.abs(eigvals)) # 所有模长必须 0.98留安全余量否则需降低 k 或增加正则化4.4 现象脉冲响应函数IRF出现非单调、振荡衰减与经济直觉冲突原因VECM 的 IRF 计算默认使用 Cholesky 分解但该分解要求变量排序体现因果链。若将inv排在gdp前意味着投资冲击直接影响 GDP忽略政策时滞。解决用广义脉冲响应GIRF替代 Cholesky IRF它不依赖排序# 使用 statsmodels 0.14 的广义 IRF irf vecm_fitted.irf(periods24, methodgeneralized) # 替代默认的 chol # 绘制 GIRF irf.plot(orthogonalizedFalse) # orthogonalizedFalse 即 GIRF plt.show()5. 进阶技巧用 VECM 做反事实政策模拟——以“消费刺激政策”为例VECM 的真正价值不在拟合而在反事实推演。例如评估“若下季度消费补贴提高 5%GDP 和投资将如何动态响应”。这需要冻结误差修正机制只释放短期动态再叠加政策冲击。以下是可直接复现的四步法5.1 步骤一构造政策冲击向量假设政策只影响消费cons其他变量初始不变。冲击大小为历史标准差的 5%# 计算消费的历史波动率 cons_std df[cons].std() policy_shock 0.05 * cons_std # 5% 的标准差冲击 # 构造冲击向量[gdp_shock, cons_shock, inv_shock] shock_vector np.array([0, policy_shock, 0])5.2 步骤二冻结误差修正项只激活短期动态VECM 的核心方程是ΔXₜ α·εₜ₋₁ Σ Γᵢ·ΔXₜ₋ᵢ uₜ反事实模拟需关闭 α·εₜ₋₁即假设系统不向长期均衡回调只保留 Γ 部分# 获取短期动态系数k2 时只有 Γ₁ Gamma vecm_fitted.gamma[0] # shape: (3, 3) # 初始化模拟数组24 期预测3 个变量 n_steps 24 simulated_diff np.zeros((n_steps, len(df.columns))) simulated_level np.zeros((n_steps, len(df.columns))) # 第一期应用政策冲击 simulated_diff[0] shock_vector simulated_level[0] df.iloc[-1].values simulated_diff[0] # 水平值 上期值 冲击 # 后续各期只用 Gamma 预测差分关闭 α 项 for t in range(1, n_steps): # ΔXₜ Γ₁ · ΔXₜ₋₁ 忽略 α·ε 和常数项 simulated_diff[t] Gamma simulated_diff[t-1] simulated_level[t] simulated_level[t-1] simulated_diff[t]5.3 步骤三加入误差修正的渐进恢复更真实完全关闭 α 不现实。更合理的做法是让 α 以衰减权重参与ΔXₜ (1-λᵗ)·α·εₜ₋₁ Γ₁·ΔXₜ₋₁ uₜ其中 λ0.95 控制恢复速度lambda_decay 0.95 # 重新初始化 simulated_diff np.zeros((n_steps, len(df.columns))) simulated_level np.zeros((n_steps, len(df.columns))) # 计算上期误差修正项 εₜ₋₁用历史数据 last_resid vecm_fitted.resid[-1] # 最后一期残差 # 注意vecm_fitted.resid 是 VECM 残差即 α·εₜ₋₁ 的估计值但我们需要 εₜ₋₁ # 从 β 重建 εₜ₋₁ β · Xₜ₋₁ X_t_minus_1 df.iloc[-1].values epsilon_t_minus_1 beta.T X_t_minus_1 # shape: (r,) for t in range(n_steps): if t 0: # 第一期全量冲击 衰减的误差修正 ec_term (1 - lambda_decay**t) * alpha epsilon_t_minus_1 simulated_diff[t] shock_vector ec_term else: # 后续短期动态 衰减误差修正 ec_term (1 - lambda_decay**t) * alpha epsilon_t_minus_1 simulated_diff[t] Gamma simulated_diff[t-1] ec_term simulated_level[t] (df.iloc[-1].values if t0 else simulated_level[t-1]) simulated_diff[t]5.4 步骤四可视化政策效果与基准对比# 基准情景无政策冲击只用历史均值预测 baseline_diff np.tile(np.mean(vecm_fitted.resid[-24:], axis0), (n_steps, 1)) baseline_level np.cumsum(baseline_diff, axis0) df.iloc[-1].values # 绘制 GDP 响应 plt.figure(figsize(10, 6)) plt.plot(range(1, n_steps1), simulated_level[:, 0], b-, labelPolicy Scenario (GDP)) plt.plot(range(1, n_steps1), baseline_level[:, 0], r--, labelBaseline (GDP)) plt.xlabel(Quarter) plt.ylabel(GDP Level) plt.title(GDP Response to 5% Consumption Stimulus) plt.legend() plt.grid(True) plt.show() # 计算累计效应 gdp_boost simulated_level[-1, 0] - baseline_level[-1, 0] print(f政策实施 6 季度后GDP 累计提升{gdp_boost:.2f} 单位约 {gdp_boost/df[gdp].mean()*100:.2f}%)我的习惯每次做政策模拟我必做三件事① 用不同 λ0.8, 0.9, 0.95跑敏感性分析看结论是否稳健② 把模拟结果与 VAR 模型的 IRF 对比若 VECM 的长期收敛值与 VAR 的长期方差分解矛盾说明协整关系设定有误③ 导出simulated_level到 Excel让业务同事自己拖动冲击大小滑块看效果——模型的价值在于被用起来而不是锁在代码里。希望帮到你。本文还有配套的精品资源点击获取