時系列データを扱う際に、ある変数が別の変数の予測に役立つかどうかを統計的に検証したいことがあります。このような「予測的因果関係」を検定する方法として、グレンジャー因果検定（Granger Causality Test）があります。
本記事では、グレンジャー因果検定の理論と、Pythonでの実装方法を解説します。
本記事の内容
- グレンジャー因果の定義と直感的な理解
- VARモデルに基づく検定の数学的定式化
- Pythonのstatsmodelsによる実装
グレンジャー因果とは
グレンジャー因果（Granger Causality）は、1969年にClive Grangerによって提案された因果関係の概念です。
直感的に言うと、「$X$ が $Y$ のグレンジャー因果である」とは、$Y$ の過去の値だけで予測するよりも、$X$ の過去の値も加えて予測した方が $Y$ の予測精度が向上する、ということを意味します。
重要な注意点として、グレンジャー因果は「予測的因果」であり、真の因果関係（介入に基づく因果）とは異なります。$X$ が $Y$ のグレンジャー原因であっても、$X$ が $Y$ の真の原因であるとは限りません。
数学的定式化
2つの定常時系列 $X_t$ と $Y_t$ を考えます。以下の2つのVARモデル（ベクトル自己回帰モデル）を比較します。
制限モデル（Restricted Model）: $Y$ の過去のみ使用
$$ Y_t = \alpha_0 + \sum_{j=1}^{p} \alpha_j Y_{t-j} + \epsilon_t $$
非制限モデル（Unrestricted Model）: $X$ の過去も使用
$$ Y_t = \alpha_0 + \sum_{j=1}^{p} \alpha_j Y_{t-j} + \sum_{j=1}^{p} \beta_j X_{t-j} + \eta_t $$
ここで $p$ はラグ次数です。
グレンジャー因果検定の帰無仮説と対立仮説は以下のようになります。
$$ H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0 \quad (\text{$X$ は $Y$ のグレンジャー原因でない}) $$
$$ H_1: \exists j, \; \beta_j \neq 0 \quad (\text{$X$ は $Y$ のグレンジャー原因である}) $$
検定統計量はF検定で構成されます。
$$ F = \frac{(\text{RSS}_R – \text{RSS}_U) / p}{\text{RSS}_U / (n – 2p – 1)} $$
ここで、$\text{RSS}_R$ は制限モデルの残差二乗和、$\text{RSS}_U$ は非制限モデルの残差二乗和、$n$ はデータ数です。帰無仮説のもとで、この統計量はF分布 $F(p, n – 2p – 1)$ に従います。
Pythonでの実装
合成データでの検証
まず、因果関係を持つ合成データを生成して、グレンジャー因果検定を行ってみましょう。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import grangercausalitytests
np.random.seed(42)
n = 500
# X → Y の因果関係を持つデータを生成
X = np.zeros(n)
Y = np.zeros(n)
for t in range(2, n):
X[t] = 0.7 * 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) # XがYに影響
# 可視化
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
axes[0].plot(X[:200], 'b-', linewidth=0.8)
axes[0].set_title("X (Cause)")
axes[0].set_xlabel("Time")
axes[0].grid(True, alpha=0.3)
axes[1].plot(Y[:200], 'r-', linewidth=0.8)
axes[1].set_title("Y (Effect)")
axes[1].set_xlabel("Time")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
statsmodelsでのグレンジャー因果検定
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import grangercausalitytests
np.random.seed(42)
n = 500
X = np.zeros(n)
Y = np.zeros(n)
for t in range(2, n):
X[t] = 0.7 * 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)
# DataFrameに格納（第1列がY、第2列がX）
# grangercausalitytests は「第2列 → 第1列」の因果を検定する
data = pd.DataFrame({'Y': Y, 'X': X})
print("=== X → Y のグレンジャー因果検定 ===")
result_xy = grangercausalitytests(data[['Y', 'X']], maxlag=5, verbose=True)
print("\n=== Y → X のグレンジャー因果検定 ===")
result_yx = grangercausalitytests(data[['X', 'Y']], maxlag=5, verbose=True)
合成データでは $X \to Y$ の因果関係を組み込んであるため、$X \to Y$ の検定ではp値が有意水準（通常0.05）より小さくなり、帰無仮説が棄却されます。一方、$Y \to X$ の検定ではp値が大きく、帰無仮説が棄却されないことが期待されます。
因果関係がない場合の確認
独立な2つの時系列でも検定してみましょう。
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import grangercausalitytests
np.random.seed(123)
n = 500
# 互いに独立な2つのAR(1)過程
A = np.zeros(n)
B = np.zeros(n)
for t in range(1, n):
A[t] = 0.6 * A[t-1] + np.random.normal(0, 1)
B[t] = 0.5 * B[t-1] + np.random.normal(0, 1)
data_ind = pd.DataFrame({'A': A, 'B': B})
print("=== B → A のグレンジャー因果検定（独立な系列）===")
grangercausalitytests(data_ind[['A', 'B']], maxlag=5, verbose=True)
独立な系列では、全てのラグでp値が有意水準を上回り、グレンジャー因果関係がないと判定されます。
結果の可視化
複数のラグに対するp値をグラフにまとめると、結果を見やすく整理できます。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import grangercausalitytests
np.random.seed(42)
n = 500
X = np.zeros(n)
Y = np.zeros(n)
for t in range(2, n):
X[t] = 0.7 * 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})
# p値を収集
max_lag = 10
p_values_xy = []
p_values_yx = []
for lag in range(1, max_lag + 1):
res_xy = grangercausalitytests(data[['Y', 'X']], maxlag=lag, verbose=False)
res_yx = grangercausalitytests(data[['X', 'Y']], maxlag=lag, verbose=False)
p_values_xy.append(res_xy[lag][0]['ssr_ftest'][1])
p_values_yx.append(res_yx[lag][0]['ssr_ftest'][1])
lags = list(range(1, max_lag + 1))
plt.figure(figsize=(10, 5))
plt.plot(lags, p_values_xy, 'bo-', linewidth=2, markersize=8, label='X -> Y')
plt.plot(lags, p_values_yx, 'rs--', linewidth=2, markersize=8, label='Y -> X')
plt.axhline(y=0.05, color='gray', linestyle=':', linewidth=1.5, label='alpha=0.05')
plt.xlabel("Lag")
plt.ylabel("p-value")
plt.title("Granger Causality Test: p-values by Lag")
plt.legend()
plt.grid(True, alpha=0.3)
plt.yscale('log')
plt.tight_layout()
plt.show()
$X \to Y$ のp値が全ラグで有意水準（0.05）を大きく下回り、$Y \to X$ のp値は有意水準を上回っていることが視覚的に確認できます。
グレンジャー因果検定の注意点
グレンジャー因果検定を使う際には、以下の点に注意が必要です。
- 定常性の仮定: 入力データは定常でなければなりません。非定常な場合は差分を取るなどの前処理が必要です
- ラグ次数の選択: 適切なラグ次数をAICやBICなどの情報量規準で選択する必要があります
- 見せかけの因果: 第三の変数（交絡因子）の影響で、見せかけのグレンジャー因果が検出されることがあります
- 因果の方向: グレンジャー因果は「予測的因果」であり、真の因果メカニズムを示すものではありません
まとめ
本記事では、グレンジャー因果検定について解説しました。
- グレンジャー因果とは「変数 $X$ の過去が変数 $Y$ の予測に有用かどうか」を検定する手法
- VARモデルの制限モデルと非制限モデルを比較するF検定で実現される
- statsmodelsの
grangercausalitytestsで簡単に実行できる - 定常性の仮定やラグ次数の選択に注意が必要
- 「予測的因果」であり、真の因果関係を保証するものではない