如何在Python中实现Toda-Yamamoto式Granger因果检验？

1. 残差自相关LM检验与模型动态稳定性检查 #

残差LM检验 #

拟合VAR模型后，直接用VARResults.test_serial_correlation()执行LM检验，覆盖1到12阶滞后：

# 假设已得到拟合好的VAR模型结果var_results
for lag in range(1, 13):
    lm_test = var_results.test_serial_correlation(lag)
    print(f"滞后{lag}阶LM检验:")
    print(f"统计量: {lm_test.statistic:.4f}, p值: {lm_test.pvalue:.4f}")

若所有滞后阶数的p值均大于0.05（或你设定的显著性水平），则无法拒绝“残差无自相关”原假设，模型满足序列独立性要求。

动态稳定性检查 #

VAR模型稳定的核心要求是所有特征根的模小于1，可通过以下代码验证：

# 检查稳定性
is_stable = var_results.is_stable()
print(f"模型是否稳定: {is_stable}")

# 绘制特征根分布图
var_results.plot_eigen()

若所有特征根都落在单位圆内（图中虚线圆），则模型动态稳定。

2. Johansen协整检验 #

使用statsmodels的coint_johansen函数直接实现迹检验和最大特征值检验，注意输入为水平时间序列，滞后阶数设为VAR最优滞后阶数p减1：

from statsmodels.tsa.vector_ar.vecm import coint_johansen

# 假设data为双变量水平时间序列DataFrame
p = var_results.k_ar  # 之前确定的VAR最优滞后阶数
johansen_result = coint_johansen(data, det_order=0, k_ar_diff=p-1)

# 输出迹检验结果
print("Johansen迹检验:")
print(f"原假设协整秩r | 迹统计量 | 5%临界值")
for r, stat, crit in zip(johansen_result.r0, johansen_result.lr1, johansen_result.cvt[:,1]):
    print(f"{r}            | {stat:.4f}    | {crit:.4f}")

# 输出最大特征值检验结果
print("\nJohansen最大特征值检验:")
print(f"原假设协整秩r | 最大特征值统计量 | 5%临界值")
for r, stat, crit in zip(johansen_result.r0, johansen_result.lr2, johansen_result.cvm[:,1]):
    print(f"{r}            | {stat:.4f}          | {crit:.4f}")

若统计量大于对应临界值，则拒绝原假设，认为变量间存在协整关系。

3. Toda-Yamamoto扩展VAR与Wald因果检验 #

核心逻辑：在最优滞后阶数p基础上，加入最大单整阶数d_max的滞后项拟合扩展VAR，仅对前p个滞后项的系数做Wald检验（后d_max个滞后仅用于保证渐近卡方分布）。

步骤1：确定最大单整阶数d_max #

通过ADF/KPSS检验得到两个变量的单整阶数，取最大值作为d_max。

步骤2：拟合扩展VAR模型 #

# p为最优滞后阶数，d_max为最大单整阶数
extended_lag = p + d_max

var_extended = VAR(data)
var_extended_results = var_extended.fit(extended_lag)

步骤3：执行Wald因果检验 #

用VARResults.test_causality()指定仅检验前p个滞后项：

# 检验变量X是否Granger引起Y
causality_test = var_extended_results.test_causality(
    endog='Y',  # 被解释变量
    exog='X',   # 解释变量
    kind='wald',
    lag_order=list(range(1, p+1))  # 仅检验前p个滞后
)

print("Toda-Yamamoto Granger因果检验结果:")
print(causality_test.summary())

若p值小于显著性水平，则拒绝原假设，认为X是Y的Granger原因。

完整流程示例代码 #

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tsa.vector_ar.var_model import VAR
from statsmodels.tsa.vector_ar.vecm import coint_johansen

# 1. 加载双变量时间序列数据（替换为你的数据路径）
data = pd.read_csv('your_data.csv', index_col=0, parse_dates=True)
data = data.dropna()

# 2. 检验单整阶数
def test_unit_root(series):
    adf_p = adfuller(series)[1]
    kpss_p = kpss(series)[1]
    print(f"ADF检验p值: {adf_p:.4f}")
    print(f"KPSS检验p值: {kpss_p:.4f}")
    # 根据检验结果返回单整阶数，此处示例假设d=1
    return 1

d1 = test_unit_root(data['X'])
d2 = test_unit_root(data['Y'])
d_max = max(d1, d2)

# 3. 选择VAR最优滞后阶数
var_model = VAR(data)
lag_order = var_model.select_order(maxlags=12)
print(lag_order.summary())
p = lag_order.selected_orders['aic']  # 可选择FPE/AIC/SC/HQ中的一种

# 4. 拟合VAR模型
var_results = var_model.fit(p)
print(var_results.summary())

# 5. 残差LM检验与稳定性检查
for lag in range(1, 13):
    lm_test = var_results.test_serial_correlation(lag)
    print(f"滞后{lag}阶LM检验: 统计量={lm_test.statistic:.4f}, p值={lm_test.pvalue:.4f}")

print(f"模型稳定: {var_results.is_stable()}")
var_results.plot_eigen()
plt.show()

# 6. Johansen协整检验
johansen_result = coint_johansen(data, det_order=0, k_ar_diff=p-1)
print("\nJohansen迹检验:")
print(f"协整秩r | 迹统计量 | 5%临界值")
for r, stat, crit in zip(johansen_result.r0, johansen_result.lr1, johansen_result.cvt[:,1]):
    print(f"{r}      | {stat:.4f}    | {crit:.4f}")

print("\nJohansen最大特征值检验:")
print(f"协整秩r | 最大特征值统计量 | 5%临界值")
for r, stat, crit in zip(johansen_result.r0, johansen_result.lr2, johansen_result.cvm[:,1]):
    print(f"{r}      | {stat:.4f}          | {crit:.4f}")

# 7. Toda-Yamamoto因果检验
extended_lag = p + d_max
var_extended = VAR(data)
var_extended_results = var_extended.fit(extended_lag)

# X→Y的因果检验
test_xy = var_extended_results.test_causality('Y', 'X', kind='wald', lag_order=list(range(1, p+1)))
print("\nX→Y的Toda-Yamamoto检验结果:")
print(test_xy.summary())

# Y→X的因果检验
test_yx = var_extended_results.test_causality('X', 'Y', kind='wald', lag_order=list(range(1, p+1)))
print("\nY→X的Toda-Yamamoto检验结果:")
print(test_yx.summary())

内容的提问来源于stack exchange，提问作者chaserone10