【数据分析】Granger因果关系分析-EW帮帮网

微信交流群：https://732903.hlcode.cn/?id=NK1fWUR

文章目录

一. Granger 因果关系的基本概念
二. Granger 因果关系的数学表述
三. Granger 因果关系的检验步骤
四. Granger 因果关系的局限性
五. 实际应用
六、Python 示例
七、完整代码和数据获取

Granger（格兰杰）因果关系分析（Granger Causality Analysis）是一种用于检验时间序列变量之间因果关系的统计方法。该方法由诺贝尔经济学奖得主 Clive Granger 提出，被广泛应用于经济学、金融学、神经科学等领域，用于揭示两个或多个时间序列变量之间的因果关系。

一. Granger 因果关系的基本概念

Granger 的因果关系定义基于预测能力，而非传统意义上的“原因-结果”关系。

在时间序列情形下，两个经济变量X、Y之间的格兰杰因果关系定义为：若在包含了变量X、Y的过去信息的条件下，对变量Y的预测效果要优于只单独由Y的过去信息对Y进行的预测效果，即变量X有助于解释变量Y的将来变化，则认为变量X是引致变量Y的格兰杰原因。

举个例子：

要预测商品的销售额，使用历史销售额 + 历史销售量来预测，效果比只用历史销售额来预测未来销售额效果好，则说销售量Granger导致了销售额。

二. Granger 因果关系的数学表述

Granger 因果关系通常通过以下两个回归模型来检验：

$X_t = \sum_{i=1}^{p} a_i X_{t-i} + \epsilon_{1t}$

$X_t = \sum_{i=1}^{p} a_i X_{t-i} + \sum_{j=1}^{q} b_j Y_{t-j} + \epsilon_{2t}$

其中， $p$ 和 $q$ 为滞后阶数， $\epsilon_{1t}$ 和 $\epsilon_{2t}$ 为白噪声项（也就是实际值和预测值之间的误差，随机分布的，没有任何趋势）。

若第二个模型显著优于第一个模型，表明 $Y$ 的历史值能显著提高对 $X$ 的预测能力，则认为 $Y$ Granger 导致 $X$ 。

那么，2个时间序列之间，就有4种可能的关系，X Granger 导致 Y、Y Granger 导致 X、互为因果、相互独立。

三. Granger 因果关系的检验步骤

3.1 时间序列平稳性检验

进行格兰杰因果关系检验的一个前提条件是时间序列必须具有平稳性，否则可能会出现虚假回归问题。因此在进行格兰杰因果关系检验之前首先应对各指标时间序列的平稳性进行单位根检验(unit root test)。常用增广的迪基—富勒检验(ADF检验)来分别对各指标序列的平稳性进行单位根检验。

时间序列的平稳性（Stationarity）是指时间序列的统计特性（如均值、方差和自相关性）不随时间变化而变化。

换句话说，如果一个时间序列是平稳的，那么它的这些统计特性在时间序列的任意子区间内都是稳定的。平稳性是时间序列分析中的一个重要概念，因为它允许使用某些统计方法和模型来分析和预测数据。

并不是说时间序列的方差、均值是一个固定值，这几乎不可能，可以有小幅度的变化，但不会随时间的增加而显著地增加或减少。

常用来检测时间序列平稳性的方法是ADF检验。

ADF检验，即Augmented Dickey-Fuller检验，是一种用于检测时间序列数据是否存在单位根的统计方法，也就是判断时间序列是否平稳。单位根的存在意味着时间序列数据是非平稳的，这将直接影响许多时间序列分析方法的有效性。ADF检验通过自回归模型来检验时间序列数据的残差（即误差）是否存在单位根，如果残差存在单位根，则时间序列数据是非平稳的；否则，数据是平稳的。

单位根指的是在自回归模型（AR模型）中，如果滞后一期的系数（即自回归系数）等于1，即存在一个根为1的多项式，这时我们称该时间序列具有单位根。

ADF检验的原假设是存在单位根，即时间序列是非平稳的。如果检验结果拒绝了原假设，即在给定的显著性水平下，可以认为时间序列是平稳的。在实际应用中，ADF检验的结果包括统计量（如t值）、p值和临界值。如果p值小于显著性水平（如0.05），则可以拒绝原假设，认为时间序列数据是平稳的。

如果时间序列不平稳，需要进行一阶差分，即将序列中每一项与其前一项相减。差分得到的结果就是原始数据的波动情况，当然原始数据的许多特征也会随之消失。

当我们对时间序列进行差分后，分析的是差分序列之间的因果关系。这种因果关系反映了原始数据的增量或变化之间的关系，而不是原始水平数据之间的关系。

短期因果关系：差分后的 Granger 因果关系可以揭示变量之间的短期动态关系，即一个变量的变化是否会引发另一个变量的变化。
长期关系的丧失：差分处理可能会丧失原始数据中的长期趋势（上升、下降）信息，因此差分后的 Granger 因果关系分析更多反映的是变量之间的短期因果关系，而非长期因果关系。

在某些情况下，两个或多个不平稳的时间序列可能在某种线性组合下是平稳的，这种情况称为协整。如果存在协整关系，那么即使原始数据不平稳，也可以通过协整检验结合误差修正模型（ECM）来分析数据间的长期因果关系。误差修正模型可以捕捉变量之间的长期均衡关系，同时考虑短期动态调整。

如果你还有研究序列之间的长期因果关系，可以用此法；

如果一阶差分后不平稳，亦可考虑分析长期因果关系。

进行Granger分析的时间序列，如果有一个不平稳，就要对所有的时间序列进行差分，你不能拿一个时间序列的波动情况和另一个时间序列进行因果关系分析。

3.2 确定滞后阶数

首先，需要确定模型的滞后阶数 $p$ 和 $q$ （即使用过去的多少个值来预测t时刻的值，好像一个滑动的窗口一样，窗口右侧第一个值根据窗口内的几个值来预测，窗口逐渐向右滑动，就可以以此预测未来的多个值，窗口的长度就是滞后阶数）。

常用的方法包括赤池信息准则（AIC）和贝叶斯信息准则（BIC）。选择使信息准则最小化的滞后阶数，即为最优滞后阶数。

AIC（Akaike Information Criterion，赤池信息准则）和 BIC（Bayesian Information Criterion，贝叶斯信息准则）是两种用于模型选择的统计量，它们衡量模型在拟合数据的精度和复杂度之间的平衡。在时间序列分析和其他统计模型中，AIC 和 BIC 经常用于选择最优模型的滞后阶数或最优模型。

（1）AIC（赤池信息准则）

AIC 是由日本统计学家赤池弘次提出的，用于评估模型的优劣。AIC 的计算公式为：

$\text{AIC} = 2k - 2\ln(L)$

其中：

$k$ 是模型中的参数数量（即模型的自由度）。
$L$ 是模型的似然函数值， $\ln(L)$ 是该值的对数。

（极大）似然估计是根据观测值（即已知的数据），来找到模型的最优参数（模型使用这些参数能够最大概率地预测出观测值）。在进行Granger因果关系分析时，使用VAR（向量自回归）模型来确定滞后阶数。VAR 模型允许一个变量不仅依赖于自身的滞后值，还依赖于其他变量的滞后值。 Granger 模型本身并不是一种独立的模型，而是基于时间序列模型（如 VAR 模型）进行的因果关系分析。所以将使得VAR模型的AIC或BIC取得最小值时对应的滞后阶数作为Granger因果分析的滞后阶数。

（2）BIC（贝叶斯信息准则）

BIC 是基于贝叶斯推断理论，由 Schwarz 提出的。它与 AIC 类似，但在模型复杂度的惩罚项中对参数数量进行了更严格的惩罚。BIC 的计算公式为：

$\text{BIC} = \ln(N)k - 2\ln(L)$

其中：

$N$ 是样本的数量。
$k$ 和 $L$ 的定义同 AIC。

AIC 和 BIC 的比较：

惩罚力度：BIC 对参数数量的惩罚力度比 AIC 更大（BIC 的惩罚项为 $\ln(N)k$ ，而 AIC 为 $2 k$ ）。因此，在样本量较大时，BIC 更倾向于选择更简单的模型。
模型选择：在模型选择中，AIC 和 BIC 都倾向于选择 AIC 或 BIC 值最小的模型。由于 BIC 更偏向于惩罚复杂模型，因此它通常会选择更简单的模型。
应用场景：
- AIC：更适用于探索模型时，尤其是在样本量相对较小时，可以容忍较多参数以更好地拟合数据。
- BIC：更适用于实际应用中，尤其是在样本量较大时，倾向于选择更简单、更稳健的模型。

3.3 构建回归模型

在确定滞后阶数后，构建包含和不包含因变量 $Y$ 的回归模型，并计算残差平方和。

3.4 F 检验

通过 F 检验比较两个模型的拟合优度。F 统计量的计算公式为：

$\frac{(RSS_1 - RSS_2) / q}{RSS_2 / (N - p - q)}$

其中， $RSS_1$ 和 $RSS_2$ 分别为两个模型的残差平方和， $N$ 为样本大小。根据 F 统计量和临界值判断是否存在 Granger 因果关系。

四. Granger 因果关系的局限性

尽管 Granger 因果关系分析在时间序列分析中非常有用，但其存在以下局限性：

不对称性：Granger 因果关系是基于预测能力的方向性分析，无法处理同步事件。
平稳性假设：Granger 因果关系假设时间序列是平稳的，但现实中许多时间序列并不平稳。
间接关系：Granger 因果关系无法区分直接和间接因果关系。

五. 实际应用

Granger 因果关系分析常用于经济学变量分析，例如检验股市之间的相互影响，或货币政策对经济增长的因果影响。在神经科学中，该方法用于分析不同脑区之间的相互影响，以理解大脑的功能连接。

六、Python 示例

主要模块：statsmodels

import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.api import VAR

以上证指数为例(一年数据)，假设我要分析它的成交量和成交额之间的关系。（你也可以分析不同指数之间的关系）
在这里插入图片描述

主要步骤：

（1）数据读取

df = pd.read_csv('上证指数.txt', encoding='ansi', sep=r'\s+', engine='python')

（2）ADF平稳性检验

result = adfuller(df['成交量'])

返回结果：

(-4.97912147654726, 2.440447767294336e-05, 1, 301, {'1%': -3.452263435801039, '5%': -2.871190526189069, '10%': -2.571911967527952}, 10846.097320776977)

分别是：ADF 统计量、P值、滞后阶数、样本数量、不同显著性水平下的P值的上限、信息准则最佳值。

通常只需要验证p值是否小于0.05即可。

（3）差分（不平稳时）

df= df.diff().dropna()  # 对所有列进行差分，去掉差分后的NaN

（4）选择滞后阶数

# 初始化VAR类
model = VAR(df)
# 模型拟合
lag_order = model.select_order()

select_order 方法用于基于给定的数据和指定的信息准则（如 AIC、BIC 等）来计算最佳滞后阶数。

参数：

maxlags (int):
- 指定允许的最大滞后阶数。如果不指定（即设置为 None），默认值为 $12 \times \left(\frac{nobs}{100}\right)^{\frac{1}{4}}$ ，其中 nobs 是样本数量。这是一个经验法则，用于在数据量较大时设置合理的最大滞后阶数。（比如300个样本默认最大滞后阶数就是15）
trend (str):
- 用于指定模型中包含的确定性项（例如常数项、线性趋势等）。可以选择以下选项：
  - "n": 无确定性项（不包括常数项、线性或二次趋势）。
  - "c": 包含常数项（默认选项）。
  - "ct": 包含常数项和线性趋势。
  - "ctt": 包含常数项、线性趋势和二次趋势。

返回值：

selections (LagOrderResults):
- 返回一个 LagOrderResults 对象，该对象包含基于不同信息准则（如 AIC、BIC、FPE、HQIC）计算的最佳滞后阶数。

结果打印示例：

<statsmodels.tsa.vector_ar.var_model.LagOrderResults object. Selected orders are: AIC -> 4, BIC -> 1, FPE -> 4, HQIC ->  2>

显示了不同信息准则最优时对应的最佳滞后阶数。

比如你要选择里面最小的阶数：

lag = min(lag_order.aic, lag_order.bic, lag_order.hqic, lag_order.fpe)

（5）Granger因果关系分析

result = sm.tsa.stattools.grangercausalitytests(df, maxlag=lag)

grangercausalitytests 可以测试两个时间序列之间是否存在因果关系。

参数：

df (array-like or pandas DataFrame):
- 输入的数据。通常是一个包含两个时间序列的数组或 DataFrame。第一个列通常被认为是被解释变量，第二列是解释变量。
maxlag (int):
- 指定最大滞后阶数。grangercausalitytests 函数将对 1 到 maxlag 的每个滞后阶数进行检验。
addconst (bool, optional, default=True):
- 是否在回归模型中添加常数项（截距项）。
verbose (bool, optional, default=True):
- 如果设置为 True，函数将打印详细的检验结果，包括每个滞后阶数的统计量和 p 值。如果设置为 False，则不打印结果。

返回值:

result (dict):
- 返回一个字典，键为滞后阶数，值为每个滞后阶数对应的检验结果。每个滞后阶数的结果也是一个字典，包含以下内容：
- params_ftest:
  - params_ftest[0]: F 检验的统计量值。
  - params_ftest[1]: F 检验的 p 值。
  - params_ftest[2]: F 检验的自由度。
  - params_ftest[3]: F 检验的回归残差平方和。
- ssr_ftest:
  - ssr_ftest[0]: 基于残差平方和的 F 检验统计量值。
  - ssr_ftest[1]: F 检验的 p 值。
  - ssr_ftest[2]: F 检验的自由度。
  - ssr_ftest[3]: F 检验的回归残差平方和。
- ssr_chi2test:
  - ssr_chi2test[0]: 基于残差平方和的卡方检验统计量值。
  - ssr_chi2test[1]: 卡方检验的 p 值。
  - ssr_chi2test[2]: 卡方检验的自由度。
  - ssr_chi2test[3]: 卡方检验的回归残差平方和。
- lrtest:
  - lrtest[0]: 似然比检验的统计量值。
  - lrtest[1]: 似然比检验的 p 值。
  - lrtest[2]: 似然比检验的自由度。
- params_ftest:
  - params_ftest[0]: 基于参数的 F 检验统计量值。
  - params_ftest[1]: F 检验的 p 值。
  - params_ftest[2]: F 检验的自由度。
  - params_ftest[3]: F 检验的回归残差平方和。
每个滞后阶数的结果都会返回多个统计量和 p 值，用于判断两个时间序列之间的 Granger 因果关系。
如果某个滞后阶数的 p 值小于显著性水平（如 0.05），则可以认为解释变量在该滞后阶数下对被解释变量具有 Granger 因果关系。

示例：

        # 屏蔽输出
        with contextlib.redirect_stdout(io.StringIO()):
            result = sm.tsa.stattools.grangercausalitytests(data_to_analysis, maxlag=self.lag_order)
        # 检查每个滞后阶数的 p 值
        for lag in range(1, self.lag_order + 1):
            test_result = result[lag][0]['ssr_ftest']  # 获取 F 检验结果
            p_value = test_result[1]  # p 值
            if p_value < 0.05:
                significant_lags.append(lag)
        return len(significant_lags) > 0

只要返回值为1，就代表df中第二列是第一列的因。

通常你需要，互换两列时间序列，判断他们是否互为Granger因果。

比如：
在这里插入图片描述

七、完整代码和数据获取

前面已经说的很明白了，如果你需要完整代码和数据示例，可以获取：

https://mbd.pub/o/bread/ZpmamZxt

￥2

在这里插入图片描述

【数据分析】Granger因果关系分析