Python实战：格兰杰因果检验全流程详解-北京尧图网络

做时序预测、特征筛选或者归因分析的朋友大概率都绕不过“格兰杰因果关系”这个词。之前两篇把理论部分铺开了这篇开始进入真正的Python实战环节。这一篇会用一组实际数据从平稳性检验、滞后阶数确定到格兰杰因果检验跑通全流程并详细解释每个步骤背后的逻辑以及如何解读结果。不管你是刚接触时序分析的新手还是已经能熟练调用statsmodels但想知道“为什么这么做”的进阶玩家这篇都适合花十几分钟认真读一遍。1. 动手之前先把格兰杰因果检验的底摸清1.1 格兰杰因果检验到底在检验什么先说人话格兰杰因果检验的本质是验证一个时间序列的历史信息是否对另一个时间序列的当前值有统计上显著的预测能力。举个例子我现在有两个序列A是某商场每天的人流量B是该商场的销售额。我想知道“人流量是不是销售额的格兰杰原因”实际上是在做这样一个对比模型一只用销售额的滞后值历史数据来预测今天的销售额模型二在模型一的基础上加入人流量的滞后值再预测今天的销售额。如果模型二的预测效果显著优于模型一那就说明人流量的历史信息确实能提供关于销售额未来变化的额外信息此时我们就说“人流量是销售额的格兰杰原因”。这里有一个特别容易踩的认知误区格兰杰因果不等于真正意义上的业务因果。它只是从统计预测能力的角度说明“先发生的事有助于预测后发生的事”。比如公鸡打鸣和太阳升起之间存在格兰杰因果关系但公鸡并不能“导致”太阳升起。所以在实际业务里格兰杰因果检验通常作为特征筛选、变量间关联探索的参考工具而不是直接下因果结论的判决书。1.2 为什么先要“平稳”这个前提条件理论上格兰杰因果检验的基础是向量自回归模型VAR而VAR模型的参数估计、F检验等统计推断都建立在时间序列平稳的假设之上。这里的“平稳”指的是宽平稳序列的均值和方差不随时间的推移发生系统性变化且任意两个时间点之间的协方差只依赖时间间隔而不依赖具体的时间点。如果不平稳会出现什么问题最典型的场景就是两个都带趋势的序列比如逐月增长的销售金额和逐月递增的广告投放费用。这两个变量即使毫无业务关系也会因为“都随时间上涨”而在回归中出现显著相关。格兰杰因果检验在这种情况下很容易产生伪因果——明明没有关联检验结果却显示有显著因果关系。所以正式检验之前必须先做平稳性判断。实践中判断平稳性的主流方法是单位根检验最常用的是ADF检验。后面我会用代码演示如何解读ADF的输出结果。1.3 Python里做格兰杰检验需要哪些库这一篇的操作是基于statsmodels这个库完成它把ADF检验、VAR模型、格兰杰因果检验都封装得很完善不用自己造轮子。statsmodels.tsa.stattools.adfuller进行ADF单位根检验statsmodels.tsa.stattools.grangercausalitytests进行格兰杰因果检验pandas用于数据处理与结构整理matplotlib用于可视化虽然不参与运算但图表能帮你在检验前快速感知序列特征。安装环境不多说了pip install pandas matplotlib statsmodels一条命令搞定。值得注意的是不同版本的statsmodels在格兰杰检验的输出格式上有一些细微差异这部分我在“常见问题”里具体讲。2. 准备数据从业务场景到DataFrame2.1 先构造一组适合演示的数据用真实案例、真实数据演示永远是最直观的。这里我构造一组模拟数据模拟的是某电商平台“每日广告曝光量”和“每日成交订单量”这两个序列时间跨度是180天。广告曝光量这部分我给它设计了一个小幅上升的趋势加上一定的周期波动和随机噪声模拟真实投放中素材生命周期、周末效应等带来的波动。成交订单量则设计为由一部分“广告曝光带来的拉动”和一部分“自身惯性”构成同时叠加随机噪声。这样设计的意图很明确构造出一组理论上存在格兰杰因果关系的数据方便后面验证检验结果是否与设计的业务逻辑吻合。import numpy as np import pandas as pd np.random.seed(42) dates pd.date_range(start2023-01-01, periods180, freqD) # 广告曝光量趋势周期性噪声 t np.arange(180) exposure 5000 15 * t 200 * np.sin(2 * np.pi * t / 7) np.random.normal(0, 80, 180) # 成交订单量受曝光滞后影响自身惯性噪声 orders np.zeros(180) orders[0] 200 for i in range(1, 180): # 用前一天的曝光量影响今天的订单量 orders[i] 0.6 * orders[i - 1] 0.02 * exposure[i - 1] np.random.normal(0, 15) df pd.DataFrame({exposure: exposure, orders: orders}, indexdates)这里注意一个细节构造数据时我对订单量的处理用的是**“前一天曝光影响今日订单”**的结构设计也就是真实滞后阶数为1。这样跑检验的时候心里有底知道理论上的滞后阶数应该在哪里再去看statsmodels给出来的结果就能一目了然判断工具是否用对了。2.2 数据清洗与结构性处理时序数据最怕的就是脏数据。订单量偶尔会出现负值或者异常缺失这在电商场景里非常常见比如退款、系统漏单、节假日流量暴涨等。为了演示我先故意在数据里加入两个异常值再演示标准的处理流程。# 人为制造两个异常值 df.loc[df.index[50], orders] -30 df.loc[df.index[90], exposure] 300 # 处理负值用前后均值替代异常低值同理 df[orders] df[orders].mask(df[orders] 0) df[orders] df[orders].fillna((df[orders].shift(1) df[orders].shift(-1)) / 2) df[exposure] df[exposure].mask(df[exposure] 1000) df[exposure] df[exposure].fillna((df[exposure].shift(1) df[exposure].shift(-1)) / 2)在处理异常值的时候有一个经验不要轻易删除整条数据。时间序列最大的价值在于顺序信息删除一个点等于丢掉了这个点的顺序关系还会让后续的滞后运算出现断裂。使用前后均值替代通常是比较稳妥的做法虽然会略微降低数据本身的波动信息但不会对整体结构造成破坏。2.3 画图看趋势心里有底再检验任何时序分析的第一步都应该是可视化。眼睛先看一遍数据判断有没有明显趋势、周期序列之间的大致走势关系这样后续分析的方向心里有底。特别是做格兰杰检验之前画图能提前预判“会不会有伪回归风险”。import matplotlib.pyplot as plt fig, axes plt.subplots(2, 1, figsize(12, 7), sharexTrue) axes[0].plot(df.index, df[exposure], color#4C72B0) axes[0].set_title(Daily Ad Exposure) axes[1].plot(df.index, df[orders], color#DD8452) axes[1].set_title(Daily Orders) plt.tight_layout() plt.show()从图上观察如果“exposure”和“orders”都体现出明显的上升趋势那就意味着两者的水平值很可能不平稳直接拿去做格兰杰检验大概率得到伪因果。这也是为什么我在文章里反复强调画图不只是为了展示更是为了提前预判后续检验中可能遇到的坑。3. 平稳性检验ADF检验实操3.1 ADF检验的核心原理与判定逻辑ADF检验全称是Augmented Dickey-Fuller检验本质上是回归检验方程$$\Delta y_t \alpha \beta t \gamma y_{t-1} \sum_{i1}^{p} \phi_i \Delta y_{t-i} \varepsilon_t$$核心要看的参数是γ方程中y_{t-1}的系数。如果γ显著为0说明序列含有单位根即序列非平稳。如果γ显著小于0说明序列是平稳的。对于ADF检验的输出结果重点看两样东西P值如果P值小于显著性水平通常取0.05拒绝“存在单位根”的原假设说明序列平稳检验统计量把它和1%、5%、10%显著性水平下的临界值对比如果统计量小于临界值同样说明平稳。3.2 用adfuller对两个序列做单位根检验实际操作很简单statsmodels里的adfuller函数封装得很完善from statsmodels.tsa.stattools import adfuller for col in [exposure, orders]: result adfuller(df[col]) print(f{col}: ADF Statistic {result[0]:.4f}, p-value {result[1]:.4f}) print(Critical Values:, {k: round(v, 4) for k, v in result[4].items()})像我构造的这组数据“exposure”带了明显的上升趋势ADF检验大概率不会拒绝单位根假设P值会明显大于0.05。“orders”本身受自身惯性影响往往也呈现非平稳特征。这就是典型的不平稳序列检验结果。在实际工作中有时P值在0.05边界附近徘徊这时不要急着下结论。我的建议是多做几个检验交叉验证比如KPSS检验它的原假设与ADF相反——原假设是序列平稳。如果你做完ADF和KPSS得到的结果互相矛盾那说明序列很可能是差分平稳的需要通过差分转换后再检验。3.3 不平稳怎么办差分处理的完整细节对于非平稳序列最常用的手段是差分。一阶差分含义直观每个时刻的数值与上一时刻数值的差即每天的“增量”。df[exposure_diff] df[exposure].diff() df[orders_diff] df[orders].diff() df df.dropna() for col in [exposure_diff, orders_diff]: result adfuller(df[col]) print(f{col}: ADF Statistic {result[0]:.4f}, p-value {result[1]:.4f})绝大多数经济类、业务类的时序数据一阶差分后都能变得平稳。如果一阶差分后仍然不平稳那就需要尝试二阶差分但这种情况越少越好——差分阶数越高原始数据的信息损耗越严重模型解释起来也越困难。差分处理之后检验的就是“增量与增量之间”的格兰杰因果关系了。比如原来我们问“广告曝光量能否预测订单量”差分后的问题就变成了“广告曝光的变化量能否预测订单量的变化量”。这两个问题的业务含义有本质区别解读结果时千万不能搞混。3.4 关于差分解读的一个高频问题很多人问我差分后做格兰杰因果检验是不是就已经偏离了最初的业务问题我的看法是格兰杰因果检验的本质是预测能力的检验它依赖于“平稳序列”这个统计前提。如果你研究的对象本身是存量、价格这类带有累积性质的序列用原始值做检验得到的结果往往不可信。差分后再检验是在“增量层面”上验证预测能力这在金融、宏观经济分析中都是非常标准的做法。如果你实在不希望差分还有一种思路是协整分析。如果两个非平稳序列之间存在协整关系——即它们的线性组合是平稳的——那么即使各自不平稳也可以建立误差修正模型来考察长期均衡关系和短期动态调整。不过那是另一个层面的方法了本篇先锚定在格兰杰检验框架内。4. 滞后阶数选择与格兰杰检验实操4.1 滞后阶数是格兰杰检验里最容易出问题的参数格兰杰因果检验对滞后阶数的敏感程度远超大多数人的预期。选择不同的滞后阶数有时得到完全相反的结论。要理解这一点需要回顾格兰杰检验的数学本质其实它就是在做受约束的回归F检验用y的滞后项回归y得到残差平方和RSS_restricted用y的滞后项加上x的滞后项回归y得到残差平方和RSS_unrestricted构建F统计量F ((RSS_restricted - RSS_unrestricted) / p) / (RSS_unrestricted / (n - k))其中p是约束个数n是样本量k是无约束模型参数个数。滞后阶数太小比如只取1阶可能会遗漏掉变量之间更长滞后的影响滞后阶数太大又会消耗大量自由度尤其样本量不大的时候统计功效明显下降。更麻烦的是不同滞后阶数下x的历史信息对y的解释能力不一样很可能出现2阶上显著、3阶上不显著的情况。4.2 用信息准则自动选择最优滞后阶数实践中选择滞后阶数通常参考信息准则最常用的是AIC赤池信息准则和BIC贝叶斯信息准则。from statsmodels.tsa.api import VAR # 先拟合不同阶数的VAR模型选出最优滞后阶数 model_order_selection VAR(df[[exposure_diff, orders_diff]]) lag_order_results model_order_selection.select_order(maxlags10) print(lag_order_results.summary())AIC和BIC的本质都是在“拟合优度”与“模型复杂度”之间寻找平衡点。滞后阶数增加时残差通常会减小但参数个数也在增加信息准则通过惩罚项来抑制过度拟合。AIC的惩罚项比BIC小所以AIC选出来的滞后阶数往往更多模型拟合更充分BIC选出来的阶数更少模型更简洁。实际操作中不同准则选出的阶数不一致很正常。我的经验是如果AIC和BIC选的阶数一致直接采用如果不一致优先用BIC的结果。因为BIC对复杂度的惩罚更严格在小样本场景下更不容易发生过拟合。但这也不是铁律最终还是要结合业务场景判断。4.3 跑通grangercausalitytests的完整流程选定滞后阶数后就可以正式调用grangercausalitytests做检验。from statsmodels.tsa.stattools import grangercausalitytests # 检验 exposure_diff 是否是 orders_diff 的格兰杰原因 gc_result grangercausalitytests(df[[orders_diff, exposure_diff]], maxlag2, verboseTrue)调用时有一个非常容易搞混的坑一定要仔细看grangercausalitytests接收的是二维数组第一个变量被视为“结果变量”被解释变量第二个变量被视为“原因变量”解释变量。这个顺序和日常直觉正好相反。我第一次用的时候就在这里栽过跟头把顺序传反了检验出来的结果完全和预期相反。输出结果中每个滞后阶数对应两行检验ssr_ftest基于残差平方和的F检验这是最常用的检验结果ssr_chi2test基于卡方分布的检验大样本下渐进等价于F检验lrtest似然比检验params_ftest基于参数约束的F检验。实际操作中绝大多数场景下看ssr_ftest这一行的P值就够了。4.4 完整解读怎么判断“谁是因谁是果”格兰杰因果检验是方向敏感的。A是B的格兰杰原因不代表B也是A的格兰杰原因。想确认两个方向就要做两次检验把两个变量的位置对调。# 反向检验orders_diff 是否是 exposure_diff 的格兰杰原因 gc_result_reverse grangercausalitytests(df[[exposure_diff, orders_diff]], maxlag2, verboseTrue)假设正向检验曝光变化量 → 订单变化量在1%或5%水平上显著反向检验订单变化量 → 曝光变化量不显著就有比较充分的证据说明广告曝光的变化情况确实对订单量的变化具有预测能力而反向的预测关系不成立。这种“单向显著”的结果在业务上很有价值说明广告曝光可能是订单量的一个领先指标可以用来做销售预测的特征输入。如果出现双向显著业务含义就更微妙了可能是两个变量相互影响也可能是遗漏了某个同时驱动两个变量的共同因素。这时建议引入更多变量做进一步分析。4.5 多变量场景下怎么处理上面的演示是经典的双变量场景。但实际业务里往往有多个潜在原因变量比如广告曝光、搜索热度、竞品价格等都可能对订单量有预测作用。多变量场景下的标准做法是向量自回归VAR模型然后在VAR框架下做格兰杰因果检验也就是块外生性检验Granger causality / block exogeneity test。statsmodels里也有相应的接口用VARModel.test_causality()来实现。这类多变量场景涉及VAR模型的完整建模流程本篇暂不展开后面会单独用一篇文章来写多变量场景下的实操。5. 常见问题与排查技巧实录5.1 结果里P值全是NaN发生了什么这是新手最容易碰到的报错之一。出现NaN最常见的原因是输入数据里包含NaN值。差分操作会引入NaN如果差分后忘记dropna()grangercausalitytests直接报错或者输出结果全部变成NaN。我的建议是在调用检验之前养成一个强制习惯——先用print(df.isnull().sum())检查数据再用df df.dropna().reset_index(dropTrue)清理干净。5.2 statsmodels版本不同输出格式不一样grangercausalitytests在不同版本的statsmodels中返回的结构有细微差别。旧版本返回的是字典每个滞后阶数的value包含ssr_ftest、lrtest等键新版本在部分场景下也会调整输出结构。如果你在写自动化脚本、循环读取检验结果建议不要依赖verboseTrue的打印输出而是手动解析返回的数据结构gc_result grangercausalitytests(df[[orders_diff, exposure_diff]], maxlag2, verboseFalse) for lag, result in gc_result.items(): p_value result[0][ssr_ftest][1] print(fLag {lag}: SSR F-test p-value {p_value:.4f})5.3 检验结果不显著但业务上明明有关系这种情况很常见遇到时先别慌按照顺序排查第一样本量是否足够。格兰杰检验本质上是F检验样本量太小时统计功效非常有限。一般经验是数据量少于50个时间点很难检验出显著结果。第二滞后阶数是否合适。这个前面提过不同滞后阶数直接影响结论。多尝试几个阶数看看结果的稳定性。如果一个变量在多个滞后阶数下都对另一个变量显著那结论就可靠得多。第三数据是否被过度处理了。有些时候为了处理异常值用了大量平滑结果把序列里的有效波动信息也抹掉了。判断标准很简单看差分后序列的方差如果方差趋近于0说明数据已经“死”了。第四很可能你的业务逻辑里影响不是线性的。格兰杰因果检验本质上是线性框架下的预测能力检验。如果变量之间存在明显的非线性关系它会变得很不敏感。这种情况下可以考虑用非线性方法做进一步验证。5.4 稳定输出结果建议用自定义函数封装写代码次数多了我发现把完整的检验流程封装成一个函数能节省大量重复劳动。这里提供一个可以直接复制的参考版本def granger_test(df, x_col, y_col, maxlag5, methodssr_ftest): 对两个时间序列做格兰杰因果检验 df: DataFrame包含两列时间序列数据 x_col: 原因变量预测变量的列名 y_col: 结果变量被预测变量的列名 from statsmodels.tsa.stattools import grangercausalitytests data df[[y_col, x_col]].dropna() results grangercausalitytests(data, maxlagmaxlag, verboseFalse) output [] for lag, res in results.items(): p_value res[0][method][1] output.append({lag: lag, p_value: p_value, significant: p_value 0.05}) return pd.DataFrame(output)注意函数内部传入grangercausalitytests时列顺序依然是先结果变量、后原因变量这和x_col、y_col的参数顺序正好相反。用这个封装的时候按业务逻辑传参即可函数内部会把顺序调整正确。5.5 一个很容易被忽视的环节业务解释要有场景边界文章开头说过格兰杰因果检验的结论不能直接等同于业务因果。但实际工作里很多业务方看到“P值显著”就直接认定“A导致B”后面的沟通就会出问题。我的处理方法是在向业务方汇报的时候明确说明检验结论是“A的历史数据对B有显著的统计预测能力”并且在业务上下文中补充判断如果两者在时间先后顺序上符合逻辑业务传导链路也说得通那这个结果才比较有可能是真正的因果信号。6. 写在最后实战里的个人体会格兰杰因果检验做久了最大的感受是工具的调用很简单但每一步决策都不简单。我在实际项目中见过太多“拿着grangercausalitytests一顿跑输出一堆P值然后挑显著的结果写进报告”的操作这本质上和撞大运没什么区别。真正的价值在于每一步都清楚自己在做什么——差分如何改变了业务问题的含义为什么选这个滞后阶数结果在不同参数下是否稳健。如果让我给刚上手的读者一个建议那就是拿到一组新数据先画图、再查平稳性、再定阶数、再跑检验每一步都要能给自己一个说得通的理由。这句话听起来简单做到了它你的时序分析水平就已经超过了绝大多数人。下一篇会进入格兰杰因果关系的更进阶场景围绕多变量VAR框架下的因果检验和结果解释来展开。这一篇里的代码和思路只要消化透了下一篇衔接起来会很顺手。

Python实战：格兰杰因果检验全流程详解

相关推荐

MR25H40CDF+PIC32：工业现场存储选型与SPI MRAM落地

Android音频问题五层分析法：从App到Codec的系统化诊断

微信小程序在线阅读毕设全解析：Java后端+MySQL+避坑实战

更多资讯

如何实时掌握用户健康数据：Open Wearables Webhooks 完整配置与调试教程

LT9211深度解析：MIPI DSI重定时器与双路分路核心技术

一文拆解MySQL索引：B+树、回表、覆盖索引与最左匹配

PIC18F86K90使用MRAM替代EEPROM：工业仪表高频写入与掉电保存方案

AI获客系统技术选型指南：从GEO优化到数字员工架构的落地路径

机械臂抓取入门指南：从硬件选型到算法实现的关键路径

今日更新

CMOS设计三大物理效应：体效应、沟道长度调制与亚阈值导电

微信小程序在线阅读毕设全解析：Java后端+MySQL+避坑实战

Android音频问题五层分析法：从App到Codec的系统化诊断

本周精选

CMOS设计三大物理效应：体效应、沟道长度调制与亚阈值导电

微信小程序在线阅读毕设全解析：Java后端+MySQL+避坑实战

Android音频问题五层分析法：从App到Codec的系统化诊断

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？