时序分析写到现在终于到了动手跑代码的阶段。如果你是从系列上篇一路看下来的应该已经理解了格兰杰因果关系的基本思想——它探讨的不是哲学意义上“谁导致了谁”而是“一个变量的历史值能否显著提升对另一个变量未来值的预测精度”。这套理论框架搭起来了剩下的就是把它落到Python里真正跑通一组数据看懂每行输出再把这些结果翻译成业务上能听懂的结论。这篇“实践1”我准备用模拟数据做完整的流程演示原因后面会细说。整个流程会覆盖环境准备、数据构造、平稳性检验、滞后阶数确定、格兰杰因果检验、结果解读、可视化报告以及我在实际项目中踩过的一些坑。内容偏长但每一步都是串起来的建议你跟着敲一遍。1. 动笔前先把检验逻辑捋一遍H0假设与F统计量决定了代码怎么写1.1 格兰杰因果的“果”不是我们日常说的“果”在写代码之前我强烈建议先把检验逻辑再确认一遍。很多人跑完grangercausalitytests函数看到p值小于0.05就兴奋地写进报告“X导致Y”。这个表述在学术上是不严谨的。格兰杰因果检验的原假设是H0X的历史值对预测Y的未来值没有显著帮助。当p值小于显著性水平比如0.05时我们拒绝原假设结论是X的滞后项对Y具有“格兰杰意义上的预测能力”。注意这里说的是预测能力而不是因果机制。我经常拿天气预报打比方。燕子低飞和下雨之间统计上很可能存在格兰杰因果关系因为燕子低飞确实对预测下雨有帮助。但你要是说“燕子导致了下雨”就闹笑话了。格兰杰检验识别的是统计意义上的“先导-滞后”关系它为因果关系提供了必要但不充分的证据。1.2 原假设与统计量F检验和卡方检验到底在比什么statsmodels里grangercausalitytests函数返回的结果包括ssr_ftest、ssr_chi2test、lrtest、params_ftest这四组统计量很多初学者看到这四个就懵了。它们其实是从不同角度去检验同一个问题。核心逻辑是构造受约束模型仅含Y自身滞后项和不受约束模型加入X的滞后项比较两者的残差平方和是否有显著差异。ssr_ftest基于残差平方和的F检验最常用小样本下表现更稳健。ssr_chi2test基于残差平方和的卡方检验大样本下F分布趋近卡方分布。lrtest似然比检验基于极大似然函数值的差异。params_ftest直接对X的滞后项系数做F检验检验这些系数是否联合为零。实际应用中我的习惯是以ssr_ftest的结果为主其他统计量作为参考。当样本量足够大比如超过500个观测各组统计量给出的结论通常一致如果出现矛盾优先信任ssr_ftest同时检查数据是否满足前提假设。理解了这一点再回头看代码你就知道为什么函数返回的是一个包含多项统计量的字典而不是简单打印一个p值。这提醒我们统计检验不是拿一个数字交差而是要对结果的多维度信息做交叉验证。2. 环境准备与数据构造为什么我坚持先用模拟数据跑通全流程2.1 Python工具链statsmodels是主力pandas和numpy是底座格兰杰因果检验的Python实现核心依赖是statsmodels库。加上pandas做数据处理、numpy做数值计算、matplotlib做可视化这套组合足够覆盖整个分析流程。import numpy as np import pandas as pd import statsmodels.api as sm from statsmodels.tsa.stattools import adfuller, grangercausalitytests from statsmodels.tsa.api import VAR import matplotlib.pyplot as plt import seaborn as sns安装可以直接用pippip install statsmodels pandas numpy matplotlib seaborn如果你用的是Anaconda发行版statsmodels通常已经预装。检查版本的时候注意一下statsmodels 0.12以上的版本接口都比较稳定grangercausalitytests函数从0.11开始就在tsa.stattools模块下位置没变过。2.2 模拟数据设计让“真相”掌握在自己手里为什么要用模拟数据因为在真实数据上我们永远不知道真正的因果关系是什么。用模拟数据则不同数据生成过程是我自己写的我知道x和y之间到底有没有关系、方向是什么、滞后几期。这样跑完代码拿到结果后可以和“标准答案”对比验证自己的流程和解读是否正确。我构造一个VAR(1)过程np.random.seed(42) n 500 x np.zeros(n) y np.zeros(n) # x_t 0.6 * x_{t-1} epsilon_x # y_t 0.5 * y_{t-1} 0.4 * x_{t-1} epsilon_y for t in range(1, n): x[t] 0.6 * x[t-1] np.random.normal(0, 1) y[t] 0.5 * y[t-1] 0.4 * x[t-1] np.random.normal(0, 1) data pd.DataFrame({y: y, x: x})这个设计里x_{t-1}的系数是0.4也就是说x对y存在真实的格兰杰因果关系而y对x则没有影响因为x的生成过程里不含y。这就是我们预设的“正确答案”。用真实业务数据做分析时我强烈建议也先做一个模拟数据跑通全流程再去分析真实数据。这个习惯帮我排查过很多代码层面的低级错误比如列顺序放反、数据类型不对、缺失值没处理等。2.3 数据形态要求DataFrame列顺序的隐藏陷阱grangercausalitytests函数接收的输入是一个二维数组或DataFrame要求第一列是被检验的响应变量第二列是解释变量。比如grangercausalitytests(data[[y, x]], maxlag2)检验的是“x是否格兰杰引起y”。很多人第一次用的时候把列顺序弄反了检验的方向就完全反了。我在一个项目里就犯过这个错误。当时有两列数据一列是广告投放量一列是销售额。我本来想验证“广告投放是否能预测销售额”结果列顺序写反了跑出来p值极低差点得出“销售额能预测广告投放”的荒谬结论。幸好做稳健性检查时发现方向不对回去一看才发现列顺序的问题。另外需要注意时间是隐式的——函数要求数据按时间顺序排列但不会去验证索引是否为时间类型。如果你传入的数据是无序的检验结果就没有任何意义。2.4 缺失值处理不要随意插补先看缺失机制时间序列数据里缺失值很常见。grangercausalitytests函数遇到空值会直接报错所以必须提前处理。我体验过的错误做法是一上来就均值填充。如果你的数据是趋势性的均值填充会引入大量虚假信息直接改变滞后结构的计算结果。正确的步骤是先看缺失值占比和位置。如果缺失值很少比如1%-2%可以先用前向填充或线性插值。如果缺失集中在某一段比如疫情封控期间门店关闭导致销售数据为零那这个缺失本身包含业务信息直接插值反而掩盖了问题。如果数据量足够大缺失比例较高优先考虑删除缺失时段而不是插值。我自己处理金融日频数据时遇到过节假日导致的结构性缺失。这时候直接用工作日日历对齐或者使用last observation carried forward前向填充都比均值填充合理。3. 平稳性检验与预处理ADF检验决定的不仅是“能不能用”还有“怎么解读”3.1 为什么非平稳数据会让格兰杰检验“翻车”格兰杰因果检验的理论推导建立在平稳随机过程的基础上。如果序列含有单位根即非平稳回归分析会面临伪回归问题——两个完全无关的非平稳序列可能表现出非常高的相关性和显著的回归系数格兰杰检验结果也会失真。我见过有人直接拿原始的GDP和人口数据跑格兰杰检验跑出来双向因果显著。这就是典型的伪回归两个序列都有明显的时间趋势各自的自相关性导致假象。处理思路有两种第一如果序列非平稳先差分使其平稳后再检验第二如果序列之间存在协整关系即使非平稳它们之间也存在长期均衡关系此时可以用误差修正模型VECM做格兰杰检验。在“实践1”里我先把差分和平稳性流程走一遍协整的情况后面单独开篇讲。3.2 ADF检验的代码与判读规则我用adfuller函数对模拟数据做平稳性检验def adf_test(series, name): result adfuller(series, autolagAIC) print(f{name}: ADF统计量{result[0]:.4f}, p值{result[1]:.4f}) print(f临界值: 1%{result[4][1%]:.4f}, 5%{result[4][5%]:.4f}, 10%{result[4][10%]:.4f}) return result adf_test(data[x], x) adf_test(data[y], y)ADF检验的原假设是序列存在单位根非平稳。p值小于0.05时拒绝原假设认为序列平稳。因为我构造的VAR(1)过程系数都小于1模拟数据大概率是平稳的。但真实数据往往不是这样这时候需要对非平稳序列做一阶差分data[x_diff] data[x].diff().dropna() data[y_diff] data[y].diff().dropna()差分后的序列通常变得平稳可以继续做格兰杰因果检验。有一个细节需要注意差分会损失一个观测值而且改变了变量的经济含义。原始变量是水平值差分后是变化量。业务解读时要相应调整比如“广告投放量的一阶差分能预测销售额的一阶差分”而不是“广告投放量能预测销售额”。3.3 如果一步差分还不行怎么办二阶差分与对数变换的选择有些强趋势序列一步差分后仍非平稳比如带有指数增长特征的序列例如某些增长速度极快的科技指标。这时可以考虑二阶差分。但我个人不太喜欢盲目二阶差分因为每次差分都会损失信息、放大噪声。我在处理宏观经济数据时的经验是先看数据是否有明显的季节性如果有先用季节差分如果没有季节性再考虑对数变换加一阶差分。取对数能把乘法关系变成加法关系同时压缩方差对很多经济金融序列很有效。判断标准始终是变换后的序列是否平稳、变换后的结果是否还能对应到有实际业务含义的解释。有个容易被忽略的点是经过差分的序列在做格兰杰因果检验时需要对滞后阶数重新进行选择不能再沿用原序列的滞后阶数。差分改变了数据的时间依赖结构原来滞后2期显著差分后可能滞后1期就够用了。这一步很多人会忘导致检验结果不稳定。4. 滞后阶数选择信息准则比拍脑袋可靠但也别完全迷信4.1 为什么滞后阶数会直接改变检验结论格兰杰因果检验的核心是比较“含X滞后项的模型”和“不含X滞后项的模型”的预测效果。滞后阶数决定了两件事如果阶数选得太小可能会遗漏真实存在于更早时期的预测关系发生“遗漏变量偏误”。如果阶数选得太大会增加待估参数数量消耗自由度降低检验效能还可能引入噪声。滞后期不同p值可能是天壤之别。我之前处理过一组传感器数据数据本身存在每12小时一次的周期性波动。如果滞后阶数选到5检验不显著选到12结果高度显著。原因是序列的自相关结构集中在12期只有滞后阶数覆盖到12时X的预测能力才能体现出来。4.2 用VAR模型的select_order选择滞后阶数一种是最简单的方法用statsmodels的VAR.fit方法配合select_order函数。model VAR(data[[x, y]]) lag_order model.select_order(maxlags12) lag_order.summary()输出会显示AIC、BIC、HQIC、FPE各个准则在不同滞后阶数下的取值。通常选择带星号的最小值。不过在实际使用中不同准则给出的最优阶数可能不一致。我的经验是AIC倾向于选择更多的滞后阶数因为它对参数数量的惩罚较弱。适合样本量较大的情况。BIC和HQIC惩罚更强选择的阶数通常更小。适合样本量有限的情况。在金融高频数据场景我倾向于看AIC选出来的阶数再结合业务上已知的周期长度做微调。4.3 结合业务周期调整滞后阶数这是一种“有根据的手动修正”纯粹按信息准则选出来的阶数并不总是合理的。比如在零售行业做销售预测促销活动的影响可能持续两周而AIC可能只选出滞后3期。这种情况下我一般会在信息准则建议值附近多试几个阶数检查结果的稳健性。操作上可以做一个stability check对滞后阶数从1到某个上限分别运行grangercausalitytests观测p值的变化趋势。如果p值在某个区间持续显著、在另一个区间持续不显著那结论相对稳健如果p值在相邻阶数间剧烈波动比如0.01跳到0.4那就要警惕结论是否取决于阶数选择。grangercausalitytests函数本身可以接受一个maxlag参数一次性返回多个阶数的检验结果为这种稳健性检查提供了便利。后面核心实战部分会展示具体用法。grangercausalitytests(data[[y, x]], maxlag4, verboseTrue)每个滞后阶数的检验结果都会打印出来。 ### 5.3 输出结果逐项解读到底看哪一行、哪个数值 grangercausalitytests的输出分为多个块每个滞后阶数一个块。以滞后阶数1为例输出的核心信息是Granger Causality number of lags (no zero) 1 ssr based F test: F83.2950 , p0.0000 , df_denom496, df_num1 ssr based chi2 test: chi283.7911 , p0.0000 , df1 likelihood ratio test: chi279.2679 , p0.0000 , df1 parameter F test: F83.2950 , p0.0000 , df_denom496, df_num1解读的优先级结构如下 1. **看ssr based F test**这是最常用的结果。p0.0000小于0.01说明在1%的显著性水平下x对y存在格兰杰因果关系。 2. **看F统计量的大小**F83.2950数值越大说明x的滞后项联合解释力越强。 3. **df_num和df_denom自由度**df_num是x滞后项的个数即滞后阶数df_denom是样本量减去参数数量。 另外需要注意ssr based F test和parameter F test的F统计量在单变量场景下数值基本一致但一个基于残差平方和的比较一个基于参数联合检验数学上等价不用困惑。 ### 5.4 逆向检验别忘了验证“反方向”是否也成立 格兰杰因果检验的一个关键特征是非对称性——x可能是y的格兰杰因但y不一定是x的格兰杰因。实际操作中两个方向都要跑一下 python # 检验 x - y result_x_to_y grangercausalitytests(data[[y, x]], maxlag2, verboseFalse) # 检验 y - x result_y_to_x grangercausalitytests(data[[x, y]], maxlag2, verboseFalse)在我们的模拟数据里x - y方向应该显著y - x方向应该不显著。这就是我们预设的“标准答案”。如果两个方向都显著有几种可能双向因果关系真实存在比如经济周期和股市指标互相驱动存在第三个共同因素在同时推动x和y或者数据存在伪回归问题需要回头检查平稳性。很多实务报告只报告显著的那个方向隐藏不显著的结果这种做法很危险。格兰杰检验的价值恰恰在于它能够帮你发现哪些“想当然的因果关系”在统计上并不成立。5.5 加入控制变量的扩展场景partial Granger causality现实问题比模拟数据复杂得多。你检验广告投入x和销售额y之间的关系时季节性因素、节假日因素、宏观经济环境都可能同时影响两者。这时需要加入控制变量做“偏格兰杰因果检验”partial Granger causality。statsmodels没有直接提供现成的partial granger函数但实现思路很直接先建立控制变量z对x和y的回归提取残差。对残差做格兰杰因果检验。from statsmodels.api import add_constant # 假设z是控制变量 # 第一步剔除z的影响 reg_x sm.OLS(data[x], add_constant(data[z])).fit() resid_x reg_x.resid reg_y sm.OLS(data[y], add_constant(data[z])).fit() resid_y reg_y.resid # 第二步对残差做格兰杰因果检验 grangercausalitytests(pd.DataFrame({y_resid: resid_y, x_resid: resid_x}), maxlag2, verboseTrue)这个思路用在多变量场景下特别实用。比如在电商场景里“促销活动”“搜索热度”“成交额”三者互相纠缠控制搜索热度之后再看促销活动对成交额是否有额外的预测能力这个结论在业务上才有指导意义。6. 进阶实践VAR框架下的多变量检验与结果可视化6.1 从二维到多维用test_causality做多变量格兰杰检验当你有三个或更多个变量时逐对做grangercausalitytests的效率会很低而且会忽略变量间的交互影响。更好的方式是构建一个VAR模型然后通过test_causality方法一次性完成多变量的因果检验。# 假设有三个变量 data_3var pd.DataFrame({y: y, x: x, w: np.random.normal(0, 1, n)}) var_model VAR(data_3var) var_result var_model.fit(maxlags2) # 检验x、w共同是否格兰杰引起y causality_all var_result.test_causality(y, [x, w], kindf) print(causality_all.summary()) # 检验仅x是否格兰杰引起y控制w的前提下 causality_x var_result.test_causality(y, x, kindf) print(causality_x.summary())test_causality方法的参数里第一个是受影响的变量effect第二个是原因的变量或变量列表cause。kind参数选f表示F检验也可以选wald做Wald检验。这个方法的优势是它在一个统一的VAR框架内能够控制所有其他变量的影响。当你说“x是y的格兰杰因”时其实是说“在控制了w的情况下x依然对y有增量预测能力”。这种表述在学术和业务报告里都更有说服力。6.2 用热力图一口气看完所有变量对之间的格兰杰关系变量一多跑出来的结果会是一堆p值。直接看数字很累我习惯把结果整理成一个矩阵然后画热力图一目了然。def granger_heatmap(data, maxlag2, significance_level0.05): n_vars data.shape[1] pval_matrix np.zeros((n_vars, n_vars)) for i in range(n_vars): col_i data.columns[i] for j in range(n_vars): if i j: pval_matrix[i, j] 1.0 continue col_j data.columns[j] # 检验列j是否格兰杰引起列i test_result grangercausalitytests(data[[col_i, col_j]], maxlagmaxlag, verboseFalse) pval_matrix[i, j] test_result[maxlag][0][ssr_ftest][1] # 画热力图 plt.figure(figsize(8, 6)) mask (pval_matrix significance_level) sns.heatmap(pval_matrix, annotTrue, fmt.3f, cmapRdYlGn_r, xticklabelsdata.columns, yticklabelsdata.columns, maskmask, cbar_kws{label: p-value}) plt.title(Granger Causality p-values) plt.show() granger_heatmap(data[[y, x]])热力图的配色逻辑是p值越小颜色越红表示格兰杰因果越显著。被mask覆盖的格子里p值大于显著性水平显示为灰色一眼扫过去就能看出哪些变量对之间存在显著的预测关系。这个方法在处理十个以内变量的时候非常高效我在项目汇报里经常直接贴这张图。6.3 样本外预测验证给格兰杰结论上“双保险”很多严谨的研究者会在格兰杰检验之后加一个样本外预测验证用不含x的模型预测y对比包含x的模型预测y看预测误差是否有显著下降。格兰杰检验是“样本内”的拟合比较而样本外验证是“预测”层面的比较两者结合结论更可靠。一个简单的实现思路from sklearn.metrics import mean_squared_error from statsmodels.tsa.ar_model import AutoReg train_len 400 y_train, y_test data[y][:train_len], data[y][train_len:] # 基准模型AR(2)只用y自身滞后项 ar_model AutoReg(y_train, lags2).fit() ar_forecast ar_model.predict(starttrain_len, endn-1) # 扩展模型VAR(2)包含x的滞后项 var_model VAR(data[[y, x]][:train_len]).fit(maxlags2) var_forecast var_model.forecast(data[[y, x]][:train_len].values[-2:], stepsn - train_len)[:, 0] mse_ar mean_squared_error(y_test, ar_forecast) mse_var mean_squared_error(y_test, var_forecast) print(fAR基准模型MSE: {mse_ar:.4f}) print(fVAR扩展模型MSE: {mse_var:.4f})如果VAR模型的MSE显著低于AR模型说明x的滞后项确实提升了y的预测精度这与格兰杰检验的结论相互印证。需要注意的是预测误差的比较需要做Diebold-Mariano检验来判断差异是否统计显著否则仅凭MSE数值大小下结论还是不够严谨。7. 我在真实项目中踩过的几个坑从数据陷阱到结果误读7.1 采样频率不一致导致的全盘皆错有次分析会员活跃度和购买金额之间的格兰杰关系会员活跃度是日度数据购买金额是周度数据。我没有对齐频率直接把两组数据丢进grangercausalitytests结果跑出了非常显著的双向因果。后来检查数据才发现周度数据被pandas自动重采样填充后产生了大量的人为平滑完全扭曲了滞后关系。正确的做法是先统一采样频率。低频数据可以用ffill或bfill对齐到高频但更好的方案是用高频数据聚合成低频聚合成周度、月度数据# 将日度数据聚合成周度 weekly_activity daily_activity.resample(W).mean() weekly_sales daily_sales.resample(W).sum()聚合时均值还是求和取决于业务含义活跃度用均值销售额用总和。这里没有标准答案业务理解永远是第一位的。7.2 异常值引发虚假的显著性格兰杰检验对异常值相当敏感。一次有个变量的某个时间点出现了极大值导致之后几期都偏离正常水平。这个异常值恰好和另一个变量的滞后项高度相关最终跑出了显著的格兰杰因果。加上一个简单的中位数滤波或winsorize处理后结论就消失了。对真实数据我建议在跑格兰杰检验之前先画出时序图肉眼扫描一遍是否有明显的异常值。如果发现异常点先确认是数据错误还是真实的极端事件数据错误就修正或剔除真实事件就考虑是否要单独建哑变量控制而不是让一个异常值主导整个回归结果。7.3 结构性突变会导致结论分段改变时间序列最常见的坑之一是结构性突变。某个政策出台、某个重大事件发生都可能让变量之间的关系发生根本性变化。整个样本期内跑出的“平均效应”可能掩盖了前后两段截然不同的关系模式。我在分析某行业数据时就遇到过2019年之前广告投入对销售额有显著预测作用2020年之后渠道结构变化导致关系弱化到不显著。整个样本一起跑格兰杰检验给出的是“不显著”。但如果把样本分段前后结论完全不同。实操上可以先画滚动窗口的格兰杰检验结果图滚动窗口回归或者用Chow检验确认变点位置再分段建模。很多工具包支持递归窗口估计但原理并不复杂取固定长度窗口逐期滑动计算每个窗口内的p值观察p值是否随时间出现结构性变化。p值在某个时点突然大幅上升往往就是结构突变的位置。7.4 结果表述这句话能不能写进报告取决于你怎么说最后强调一下结果表述的问题。格兰杰检验的结果在报告里应该这样写正确“在控制了自身滞后项和其他变量后x对y具有统计上显著的增量预测能力。”错误“x是y的原因。”较稳妥“x领先y可以作为y的前瞻指标使用。”很多业务部门看到“因果”两个字就兴奋以为找到了杠杆解——控制x就能改变y。这是对格兰杰因果的滥用。它只能告诉你x的历史值携带着关于y未来值的信息至于x是否真的“导致”了y需要经济学逻辑、实验设计或其他因果推断方法如断点回归、工具变量等来进一步确认。从我自己的使用体验来说格兰杰因果检验在实务场景里最有价值的应用其实是筛选变量、验证直觉、辅助构建预测模型。比如你有很多个候选指标想知道哪些指标对目标变量有增量预测能力跑一遍多变量格兰杰检验把不显著的剔除掉再用剩下的变量构建预测模型这个用法既合理又高效。它像是用统计方法把不必要的变量过滤掉让你有更多精力去深入分析那些真正重要的驱动因素。以上是实践1的核心内容。下一篇实践2我准备写协整与误差修正模型下的格兰杰因果检验那才是处理非平稳经济金融数据的重头戏。先用好这篇里的工具和思路处理日常80%的时序预测问题都够了。 
长河君 · 2026/10/4 1:27:11 · 阅读 1.2k Python实现格兰杰因果检验：从平稳性检验到结果解读的完整实战
时序分析写到现在终于到了动手跑代码的阶段。如果你是从系列上篇一路看下来的应该已经理解了格兰杰因果关系的基本思想——它探讨的不是哲学意义上“谁导致了谁”而是“一个变量的历史值能否显著提升对另一个变量未来值的预测精度”。这套理论框架搭起来了剩下的就是把它落到Python里真正跑通一组数据看懂每行输出再把这些结果翻译成业务上能听懂的结论。这篇“实践1”我准备用模拟数据做完整的流程演示原因后面会细说。整个流程会覆盖环境准备、数据构造、平稳性检验、滞后阶数确定、格兰杰因果检验、结果解读、可视化报告以及我在实际项目中踩过的一些坑。内容偏长但每一步都是串起来的建议你跟着敲一遍。1. 动笔前先把检验逻辑捋一遍H0假设与F统计量决定了代码怎么写1.1 格兰杰因果的“果”不是我们日常说的“果”在写代码之前我强烈建议先把检验逻辑再确认一遍。很多人跑完grangercausalitytests函数看到p值小于0.05就兴奋地写进报告“X导致Y”。这个表述在学术上是不严谨的。格兰杰因果检验的原假设是H0X的历史值对预测Y的未来值没有显著帮助。当p值小于显著性水平比如0.05时我们拒绝原假设结论是X的滞后项对Y具有“格兰杰意义上的预测能力”。注意这里说的是预测能力而不是因果机制。我经常拿天气预报打比方。燕子低飞和下雨之间统计上很可能存在格兰杰因果关系因为燕子低飞确实对预测下雨有帮助。但你要是说“燕子导致了下雨”就闹笑话了。格兰杰检验识别的是统计意义上的“先导-滞后”关系它为因果关系提供了必要但不充分的证据。1.2 原假设与统计量F检验和卡方检验到底在比什么statsmodels里grangercausalitytests函数返回的结果包括ssr_ftest、ssr_chi2test、lrtest、params_ftest这四组统计量很多初学者看到这四个就懵了。它们其实是从不同角度去检验同一个问题。核心逻辑是构造受约束模型仅含Y自身滞后项和不受约束模型加入X的滞后项比较两者的残差平方和是否有显著差异。ssr_ftest基于残差平方和的F检验最常用小样本下表现更稳健。ssr_chi2test基于残差平方和的卡方检验大样本下F分布趋近卡方分布。lrtest似然比检验基于极大似然函数值的差异。params_ftest直接对X的滞后项系数做F检验检验这些系数是否联合为零。实际应用中我的习惯是以ssr_ftest的结果为主其他统计量作为参考。当样本量足够大比如超过500个观测各组统计量给出的结论通常一致如果出现矛盾优先信任ssr_ftest同时检查数据是否满足前提假设。理解了这一点再回头看代码你就知道为什么函数返回的是一个包含多项统计量的字典而不是简单打印一个p值。这提醒我们统计检验不是拿一个数字交差而是要对结果的多维度信息做交叉验证。2. 环境准备与数据构造为什么我坚持先用模拟数据跑通全流程2.1 Python工具链statsmodels是主力pandas和numpy是底座格兰杰因果检验的Python实现核心依赖是statsmodels库。加上pandas做数据处理、numpy做数值计算、matplotlib做可视化这套组合足够覆盖整个分析流程。import numpy as np import pandas as pd import statsmodels.api as sm from statsmodels.tsa.stattools import adfuller, grangercausalitytests from statsmodels.tsa.api import VAR import matplotlib.pyplot as plt import seaborn as sns安装可以直接用pippip install statsmodels pandas numpy matplotlib seaborn如果你用的是Anaconda发行版statsmodels通常已经预装。检查版本的时候注意一下statsmodels 0.12以上的版本接口都比较稳定grangercausalitytests函数从0.11开始就在tsa.stattools模块下位置没变过。2.2 模拟数据设计让“真相”掌握在自己手里为什么要用模拟数据因为在真实数据上我们永远不知道真正的因果关系是什么。用模拟数据则不同数据生成过程是我自己写的我知道x和y之间到底有没有关系、方向是什么、滞后几期。这样跑完代码拿到结果后可以和“标准答案”对比验证自己的流程和解读是否正确。我构造一个VAR(1)过程np.random.seed(42) n 500 x np.zeros(n) y np.zeros(n) # x_t 0.6 * x_{t-1} epsilon_x # y_t 0.5 * y_{t-1} 0.4 * x_{t-1} epsilon_y for t in range(1, n): x[t] 0.6 * x[t-1] np.random.normal(0, 1) y[t] 0.5 * y[t-1] 0.4 * x[t-1] np.random.normal(0, 1) data pd.DataFrame({y: y, x: x})这个设计里x_{t-1}的系数是0.4也就是说x对y存在真实的格兰杰因果关系而y对x则没有影响因为x的生成过程里不含y。这就是我们预设的“正确答案”。用真实业务数据做分析时我强烈建议也先做一个模拟数据跑通全流程再去分析真实数据。这个习惯帮我排查过很多代码层面的低级错误比如列顺序放反、数据类型不对、缺失值没处理等。2.3 数据形态要求DataFrame列顺序的隐藏陷阱grangercausalitytests函数接收的输入是一个二维数组或DataFrame要求第一列是被检验的响应变量第二列是解释变量。比如grangercausalitytests(data[[y, x]], maxlag2)检验的是“x是否格兰杰引起y”。很多人第一次用的时候把列顺序弄反了检验的方向就完全反了。我在一个项目里就犯过这个错误。当时有两列数据一列是广告投放量一列是销售额。我本来想验证“广告投放是否能预测销售额”结果列顺序写反了跑出来p值极低差点得出“销售额能预测广告投放”的荒谬结论。幸好做稳健性检查时发现方向不对回去一看才发现列顺序的问题。另外需要注意时间是隐式的——函数要求数据按时间顺序排列但不会去验证索引是否为时间类型。如果你传入的数据是无序的检验结果就没有任何意义。2.4 缺失值处理不要随意插补先看缺失机制时间序列数据里缺失值很常见。grangercausalitytests函数遇到空值会直接报错所以必须提前处理。我体验过的错误做法是一上来就均值填充。如果你的数据是趋势性的均值填充会引入大量虚假信息直接改变滞后结构的计算结果。正确的步骤是先看缺失值占比和位置。如果缺失值很少比如1%-2%可以先用前向填充或线性插值。如果缺失集中在某一段比如疫情封控期间门店关闭导致销售数据为零那这个缺失本身包含业务信息直接插值反而掩盖了问题。如果数据量足够大缺失比例较高优先考虑删除缺失时段而不是插值。我自己处理金融日频数据时遇到过节假日导致的结构性缺失。这时候直接用工作日日历对齐或者使用last observation carried forward前向填充都比均值填充合理。3. 平稳性检验与预处理ADF检验决定的不仅是“能不能用”还有“怎么解读”3.1 为什么非平稳数据会让格兰杰检验“翻车”格兰杰因果检验的理论推导建立在平稳随机过程的基础上。如果序列含有单位根即非平稳回归分析会面临伪回归问题——两个完全无关的非平稳序列可能表现出非常高的相关性和显著的回归系数格兰杰检验结果也会失真。我见过有人直接拿原始的GDP和人口数据跑格兰杰检验跑出来双向因果显著。这就是典型的伪回归两个序列都有明显的时间趋势各自的自相关性导致假象。处理思路有两种第一如果序列非平稳先差分使其平稳后再检验第二如果序列之间存在协整关系即使非平稳它们之间也存在长期均衡关系此时可以用误差修正模型VECM做格兰杰检验。在“实践1”里我先把差分和平稳性流程走一遍协整的情况后面单独开篇讲。3.2 ADF检验的代码与判读规则我用adfuller函数对模拟数据做平稳性检验def adf_test(series, name): result adfuller(series, autolagAIC) print(f{name}: ADF统计量{result[0]:.4f}, p值{result[1]:.4f}) print(f临界值: 1%{result[4][1%]:.4f}, 5%{result[4][5%]:.4f}, 10%{result[4][10%]:.4f}) return result adf_test(data[x], x) adf_test(data[y], y)ADF检验的原假设是序列存在单位根非平稳。p值小于0.05时拒绝原假设认为序列平稳。因为我构造的VAR(1)过程系数都小于1模拟数据大概率是平稳的。但真实数据往往不是这样这时候需要对非平稳序列做一阶差分data[x_diff] data[x].diff().dropna() data[y_diff] data[y].diff().dropna()差分后的序列通常变得平稳可以继续做格兰杰因果检验。有一个细节需要注意差分会损失一个观测值而且改变了变量的经济含义。原始变量是水平值差分后是变化量。业务解读时要相应调整比如“广告投放量的一阶差分能预测销售额的一阶差分”而不是“广告投放量能预测销售额”。3.3 如果一步差分还不行怎么办二阶差分与对数变换的选择有些强趋势序列一步差分后仍非平稳比如带有指数增长特征的序列例如某些增长速度极快的科技指标。这时可以考虑二阶差分。但我个人不太喜欢盲目二阶差分因为每次差分都会损失信息、放大噪声。我在处理宏观经济数据时的经验是先看数据是否有明显的季节性如果有先用季节差分如果没有季节性再考虑对数变换加一阶差分。取对数能把乘法关系变成加法关系同时压缩方差对很多经济金融序列很有效。判断标准始终是变换后的序列是否平稳、变换后的结果是否还能对应到有实际业务含义的解释。有个容易被忽略的点是经过差分的序列在做格兰杰因果检验时需要对滞后阶数重新进行选择不能再沿用原序列的滞后阶数。差分改变了数据的时间依赖结构原来滞后2期显著差分后可能滞后1期就够用了。这一步很多人会忘导致检验结果不稳定。4. 滞后阶数选择信息准则比拍脑袋可靠但也别完全迷信4.1 为什么滞后阶数会直接改变检验结论格兰杰因果检验的核心是比较“含X滞后项的模型”和“不含X滞后项的模型”的预测效果。滞后阶数决定了两件事如果阶数选得太小可能会遗漏真实存在于更早时期的预测关系发生“遗漏变量偏误”。如果阶数选得太大会增加待估参数数量消耗自由度降低检验效能还可能引入噪声。滞后期不同p值可能是天壤之别。我之前处理过一组传感器数据数据本身存在每12小时一次的周期性波动。如果滞后阶数选到5检验不显著选到12结果高度显著。原因是序列的自相关结构集中在12期只有滞后阶数覆盖到12时X的预测能力才能体现出来。4.2 用VAR模型的select_order选择滞后阶数一种是最简单的方法用statsmodels的VAR.fit方法配合select_order函数。model VAR(data[[x, y]]) lag_order model.select_order(maxlags12) lag_order.summary()输出会显示AIC、BIC、HQIC、FPE各个准则在不同滞后阶数下的取值。通常选择带星号的最小值。不过在实际使用中不同准则给出的最优阶数可能不一致。我的经验是AIC倾向于选择更多的滞后阶数因为它对参数数量的惩罚较弱。适合样本量较大的情况。BIC和HQIC惩罚更强选择的阶数通常更小。适合样本量有限的情况。在金融高频数据场景我倾向于看AIC选出来的阶数再结合业务上已知的周期长度做微调。4.3 结合业务周期调整滞后阶数这是一种“有根据的手动修正”纯粹按信息准则选出来的阶数并不总是合理的。比如在零售行业做销售预测促销活动的影响可能持续两周而AIC可能只选出滞后3期。这种情况下我一般会在信息准则建议值附近多试几个阶数检查结果的稳健性。操作上可以做一个stability check对滞后阶数从1到某个上限分别运行grangercausalitytests观测p值的变化趋势。如果p值在某个区间持续显著、在另一个区间持续不显著那结论相对稳健如果p值在相邻阶数间剧烈波动比如0.01跳到0.4那就要警惕结论是否取决于阶数选择。grangercausalitytests函数本身可以接受一个maxlag参数一次性返回多个阶数的检验结果为这种稳健性检查提供了便利。后面核心实战部分会展示具体用法。grangercausalitytests(data[[y, x]], maxlag4, verboseTrue)每个滞后阶数的检验结果都会打印出来。 ### 5.3 输出结果逐项解读到底看哪一行、哪个数值 grangercausalitytests的输出分为多个块每个滞后阶数一个块。以滞后阶数1为例输出的核心信息是Granger Causality number of lags (no zero) 1 ssr based F test: F83.2950 , p0.0000 , df_denom496, df_num1 ssr based chi2 test: chi283.7911 , p0.0000 , df1 likelihood ratio test: chi279.2679 , p0.0000 , df1 parameter F test: F83.2950 , p0.0000 , df_denom496, df_num1解读的优先级结构如下 1. **看ssr based F test**这是最常用的结果。p0.0000小于0.01说明在1%的显著性水平下x对y存在格兰杰因果关系。 2. **看F统计量的大小**F83.2950数值越大说明x的滞后项联合解释力越强。 3. **df_num和df_denom自由度**df_num是x滞后项的个数即滞后阶数df_denom是样本量减去参数数量。 另外需要注意ssr based F test和parameter F test的F统计量在单变量场景下数值基本一致但一个基于残差平方和的比较一个基于参数联合检验数学上等价不用困惑。 ### 5.4 逆向检验别忘了验证“反方向”是否也成立 格兰杰因果检验的一个关键特征是非对称性——x可能是y的格兰杰因但y不一定是x的格兰杰因。实际操作中两个方向都要跑一下 python # 检验 x - y result_x_to_y grangercausalitytests(data[[y, x]], maxlag2, verboseFalse) # 检验 y - x result_y_to_x grangercausalitytests(data[[x, y]], maxlag2, verboseFalse)在我们的模拟数据里x - y方向应该显著y - x方向应该不显著。这就是我们预设的“标准答案”。如果两个方向都显著有几种可能双向因果关系真实存在比如经济周期和股市指标互相驱动存在第三个共同因素在同时推动x和y或者数据存在伪回归问题需要回头检查平稳性。很多实务报告只报告显著的那个方向隐藏不显著的结果这种做法很危险。格兰杰检验的价值恰恰在于它能够帮你发现哪些“想当然的因果关系”在统计上并不成立。5.5 加入控制变量的扩展场景partial Granger causality现实问题比模拟数据复杂得多。你检验广告投入x和销售额y之间的关系时季节性因素、节假日因素、宏观经济环境都可能同时影响两者。这时需要加入控制变量做“偏格兰杰因果检验”partial Granger causality。statsmodels没有直接提供现成的partial granger函数但实现思路很直接先建立控制变量z对x和y的回归提取残差。对残差做格兰杰因果检验。from statsmodels.api import add_constant # 假设z是控制变量 # 第一步剔除z的影响 reg_x sm.OLS(data[x], add_constant(data[z])).fit() resid_x reg_x.resid reg_y sm.OLS(data[y], add_constant(data[z])).fit() resid_y reg_y.resid # 第二步对残差做格兰杰因果检验 grangercausalitytests(pd.DataFrame({y_resid: resid_y, x_resid: resid_x}), maxlag2, verboseTrue)这个思路用在多变量场景下特别实用。比如在电商场景里“促销活动”“搜索热度”“成交额”三者互相纠缠控制搜索热度之后再看促销活动对成交额是否有额外的预测能力这个结论在业务上才有指导意义。6. 进阶实践VAR框架下的多变量检验与结果可视化6.1 从二维到多维用test_causality做多变量格兰杰检验当你有三个或更多个变量时逐对做grangercausalitytests的效率会很低而且会忽略变量间的交互影响。更好的方式是构建一个VAR模型然后通过test_causality方法一次性完成多变量的因果检验。# 假设有三个变量 data_3var pd.DataFrame({y: y, x: x, w: np.random.normal(0, 1, n)}) var_model VAR(data_3var) var_result var_model.fit(maxlags2) # 检验x、w共同是否格兰杰引起y causality_all var_result.test_causality(y, [x, w], kindf) print(causality_all.summary()) # 检验仅x是否格兰杰引起y控制w的前提下 causality_x var_result.test_causality(y, x, kindf) print(causality_x.summary())test_causality方法的参数里第一个是受影响的变量effect第二个是原因的变量或变量列表cause。kind参数选f表示F检验也可以选wald做Wald检验。这个方法的优势是它在一个统一的VAR框架内能够控制所有其他变量的影响。当你说“x是y的格兰杰因”时其实是说“在控制了w的情况下x依然对y有增量预测能力”。这种表述在学术和业务报告里都更有说服力。6.2 用热力图一口气看完所有变量对之间的格兰杰关系变量一多跑出来的结果会是一堆p值。直接看数字很累我习惯把结果整理成一个矩阵然后画热力图一目了然。def granger_heatmap(data, maxlag2, significance_level0.05): n_vars data.shape[1] pval_matrix np.zeros((n_vars, n_vars)) for i in range(n_vars): col_i data.columns[i] for j in range(n_vars): if i j: pval_matrix[i, j] 1.0 continue col_j data.columns[j] # 检验列j是否格兰杰引起列i test_result grangercausalitytests(data[[col_i, col_j]], maxlagmaxlag, verboseFalse) pval_matrix[i, j] test_result[maxlag][0][ssr_ftest][1] # 画热力图 plt.figure(figsize(8, 6)) mask (pval_matrix significance_level) sns.heatmap(pval_matrix, annotTrue, fmt.3f, cmapRdYlGn_r, xticklabelsdata.columns, yticklabelsdata.columns, maskmask, cbar_kws{label: p-value}) plt.title(Granger Causality p-values) plt.show() granger_heatmap(data[[y, x]])热力图的配色逻辑是p值越小颜色越红表示格兰杰因果越显著。被mask覆盖的格子里p值大于显著性水平显示为灰色一眼扫过去就能看出哪些变量对之间存在显著的预测关系。这个方法在处理十个以内变量的时候非常高效我在项目汇报里经常直接贴这张图。6.3 样本外预测验证给格兰杰结论上“双保险”很多严谨的研究者会在格兰杰检验之后加一个样本外预测验证用不含x的模型预测y对比包含x的模型预测y看预测误差是否有显著下降。格兰杰检验是“样本内”的拟合比较而样本外验证是“预测”层面的比较两者结合结论更可靠。一个简单的实现思路from sklearn.metrics import mean_squared_error from statsmodels.tsa.ar_model import AutoReg train_len 400 y_train, y_test data[y][:train_len], data[y][train_len:] # 基准模型AR(2)只用y自身滞后项 ar_model AutoReg(y_train, lags2).fit() ar_forecast ar_model.predict(starttrain_len, endn-1) # 扩展模型VAR(2)包含x的滞后项 var_model VAR(data[[y, x]][:train_len]).fit(maxlags2) var_forecast var_model.forecast(data[[y, x]][:train_len].values[-2:], stepsn - train_len)[:, 0] mse_ar mean_squared_error(y_test, ar_forecast) mse_var mean_squared_error(y_test, var_forecast) print(fAR基准模型MSE: {mse_ar:.4f}) print(fVAR扩展模型MSE: {mse_var:.4f})如果VAR模型的MSE显著低于AR模型说明x的滞后项确实提升了y的预测精度这与格兰杰检验的结论相互印证。需要注意的是预测误差的比较需要做Diebold-Mariano检验来判断差异是否统计显著否则仅凭MSE数值大小下结论还是不够严谨。7. 我在真实项目中踩过的几个坑从数据陷阱到结果误读7.1 采样频率不一致导致的全盘皆错有次分析会员活跃度和购买金额之间的格兰杰关系会员活跃度是日度数据购买金额是周度数据。我没有对齐频率直接把两组数据丢进grangercausalitytests结果跑出了非常显著的双向因果。后来检查数据才发现周度数据被pandas自动重采样填充后产生了大量的人为平滑完全扭曲了滞后关系。正确的做法是先统一采样频率。低频数据可以用ffill或bfill对齐到高频但更好的方案是用高频数据聚合成低频聚合成周度、月度数据# 将日度数据聚合成周度 weekly_activity daily_activity.resample(W).mean() weekly_sales daily_sales.resample(W).sum()聚合时均值还是求和取决于业务含义活跃度用均值销售额用总和。这里没有标准答案业务理解永远是第一位的。7.2 异常值引发虚假的显著性格兰杰检验对异常值相当敏感。一次有个变量的某个时间点出现了极大值导致之后几期都偏离正常水平。这个异常值恰好和另一个变量的滞后项高度相关最终跑出了显著的格兰杰因果。加上一个简单的中位数滤波或winsorize处理后结论就消失了。对真实数据我建议在跑格兰杰检验之前先画出时序图肉眼扫描一遍是否有明显的异常值。如果发现异常点先确认是数据错误还是真实的极端事件数据错误就修正或剔除真实事件就考虑是否要单独建哑变量控制而不是让一个异常值主导整个回归结果。7.3 结构性突变会导致结论分段改变时间序列最常见的坑之一是结构性突变。某个政策出台、某个重大事件发生都可能让变量之间的关系发生根本性变化。整个样本期内跑出的“平均效应”可能掩盖了前后两段截然不同的关系模式。我在分析某行业数据时就遇到过2019年之前广告投入对销售额有显著预测作用2020年之后渠道结构变化导致关系弱化到不显著。整个样本一起跑格兰杰检验给出的是“不显著”。但如果把样本分段前后结论完全不同。实操上可以先画滚动窗口的格兰杰检验结果图滚动窗口回归或者用Chow检验确认变点位置再分段建模。很多工具包支持递归窗口估计但原理并不复杂取固定长度窗口逐期滑动计算每个窗口内的p值观察p值是否随时间出现结构性变化。p值在某个时点突然大幅上升往往就是结构突变的位置。7.4 结果表述这句话能不能写进报告取决于你怎么说最后强调一下结果表述的问题。格兰杰检验的结果在报告里应该这样写正确“在控制了自身滞后项和其他变量后x对y具有统计上显著的增量预测能力。”错误“x是y的原因。”较稳妥“x领先y可以作为y的前瞻指标使用。”很多业务部门看到“因果”两个字就兴奋以为找到了杠杆解——控制x就能改变y。这是对格兰杰因果的滥用。它只能告诉你x的历史值携带着关于y未来值的信息至于x是否真的“导致”了y需要经济学逻辑、实验设计或其他因果推断方法如断点回归、工具变量等来进一步确认。从我自己的使用体验来说格兰杰因果检验在实务场景里最有价值的应用其实是筛选变量、验证直觉、辅助构建预测模型。比如你有很多个候选指标想知道哪些指标对目标变量有增量预测能力跑一遍多变量格兰杰检验把不显著的剔除掉再用剩下的变量构建预测模型这个用法既合理又高效。它像是用统计方法把不必要的变量过滤掉让你有更多精力去深入分析那些真正重要的驱动因素。以上是实践1的核心内容。下一篇实践2我准备写协整与误差修正模型下的格兰杰因果检验那才是处理非平稳经济金融数据的重头戏。先用好这篇里的工具和思路处理日常80%的时序预测问题都够了。