広告費を増やしたら売上は本当に伸びるのでしょうか？金利を引き上げるとインフレ率は下がるのでしょうか？あるいは、脳の前頭葉の活動が運動野の活動を「引き起こしている」と言えるのでしょうか？
こうした問いに共通するのは、2つの時系列のあいだに因果関係があるかどうかを判定したいという欲求です。しかし、観測データから真の因果関係を証明することは、一般的には不可能に近い難題です。ランダム化比較試験（RCT）を行えない状況で、私たちにできるのは「因果の手がかり」を見つけることくらいです。
1969年、経済学者クライブ・グレンジャーは逆転の発想を提案しました。「原因とは何か」という哲学的な問いを脇に置き、「ある変数の過去の情報が、別の変数の予測精度を改善するか」という予測の問題に帰着させたのです。この「予測的因果」を統計的に検定する方法がグレンジャー因果検定（Granger Causality Test）です。
グレンジャー因果検定は、以下のような幅広い分野で活用されています。
- マクロ経済学 — GDP成長率とマネーサプライ、為替レートと貿易収支の予測的因果を調べ、経済政策の効果を評価する
- 神経科学 — fMRIやEEGの信号から脳領域間の情報伝達の方向性を推定し、認知機能のネットワーク構造を解明する
- 気候科学 — 海面水温とエルニーニョ指数の時間的先行関係を分析し、気候変動のメカニズムを探る
- 金融工学 — 市場間のリード・ラグ関係を特定し、リスク伝播やポートフォリオ構築に活用する
本記事の内容
- グレンジャー因果の概念 —「予測的因果」と真の因果の違い
- VARモデルに基づく数学的定式化（制限モデル vs 非制限モデル）
- F検定統計量の導出と解釈
- 多変量システムへの拡張
- ラグ次数の選択方法（AIC, BIC）
- Pythonでのスクラッチ実装
- statsmodelsを用いた実装
- 合成データでの検証実験
- グレンジャー因果の限界と注意点
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
特にVARモデルの基本的な構造（定式化、OLS推定、安定性条件）を理解していると、本記事の数学的な議論がスムーズに読めます。
グレンジャー因果の概念 — 「予測的因果」とは何か
日常の例から考える
朝、空が赤く染まると、その日の午後に雨が降ることが多いという経験則があります。では、「朝焼けが雨を引き起こしている」と言えるでしょうか？もちろん、朝焼けそのものが雨雲を発生させるわけではありません。しかし、朝焼けの情報を知ることで、雨の予測精度が向上するのは事実です。
グレンジャー因果はまさにこの「予測に役立つかどうか」を検定します。もう少し正確に述べると、時系列 $X$ が時系列 $Y$ に対してグレンジャー因果を持つとは、以下の条件を意味します。
$Y$ の過去の値だけで予測するよりも、$X$ の過去の値も加えて予測した方が、$Y$ の予測精度が統計的に有意に向上する。
イメージとしては、友人と天気予報対決をしているようなものです。あなたは過去の気温データだけで明日の天気を予測します。一方、友人は気温データに加えて気圧データも使います。もし友人の予測が統計的に有意にあなたの予測よりも正確であれば、「気圧は気温に対してグレンジャー因果を持つ」と言います。
予測的因果と真の因果の違い
ここで注意が必要です。グレンジャー因果は「予測的因果（predictive causality）」であって、真の因果関係（true causality）とは別物です。この区別は、検定結果を解釈する上で極めて重要です。
具体例で考えてみましょう。クリスマスの時期になると、以下の2つの現象が同時に観察されます。
- デパートの売上が増加する
- 街のイルミネーション点灯数が増加する
イルミネーション点灯数の過去の値が売上の予測精度を改善するかもしれません（両者とも12月に向けて増加するため）。しかし、イルミネーションを増やしたからといって売上が増えるわけではありません。真の原因は「クリスマスシーズンの到来」という共通の要因です。このように、背後に交絡因子（confounding variable）が存在するとき、グレンジャー因果は因果関係がなくても検出されることがあります。
グレンジャー自身もこの限界を認識しており、この検定はあくまで「予測的因果」を調べるものであって、真の因果メカニズムの存在を証明するものではないと明言しています。それでもなお、この検定が広く使われているのは、観測データだけから変数間の時間的先行関係を定量的に評価できるという実用的な価値が大きいためです。
さて、グレンジャー因果の概念的な意味がわかったところで、これを数学的にどう定式化するのかを見ていきましょう。鍵となるのはVARモデルです。
数学的定式化 — 制限モデルと非制限モデル
直感的な枠組み
グレンジャー因果検定の数学的な枠組みは、実にシンプルな発想に基づいています。「$X$ の過去の情報が $Y$ の予測に役立つか」を調べたいのですから、$X$ の情報を使うモデルと使わないモデルを作り、予測精度を比較すれば良いのです。
予測精度の比較には、回帰分析でおなじみの「残差二乗和（RSS）」を使います。$X$ の情報を加えることでRSSが有意に小さくなれば、$X$ は $Y$ の予測に貢献していると言えます。
2変量のケース
2つの定常時系列 $\{X_t\}$ と $\{Y_t\}$ を考えます。「$X$ が $Y$ にグレンジャー因果を持つか」を検定するために、2つのモデルを立てます。
制限モデル（Restricted Model）: $Y$ の過去の値だけで $Y_t$ を予測するモデルです。$X$ の情報は一切使いません。
$$ Y_t = \alpha_0 + \sum_{j=1}^{p} \alpha_j Y_{t-j} + \varepsilon_t^{(R)} $$
これは単純なAR(p)モデルにほかなりません。$\alpha_0$ は定数項、$\alpha_j$ は $Y$ の $j$ 期前の値に対する係数、$\varepsilon_t^{(R)}$ は誤差項です。
非制限モデル（Unrestricted Model）: $Y$ の過去の値に加えて、$X$ の過去の値も使って $Y_t$ を予測するモデルです。
$$ Y_t = \beta_0 + \sum_{j=1}^{p} \beta_j Y_{t-j} + \sum_{j=1}^{p} \gamma_j X_{t-j} + \varepsilon_t^{(U)} $$
ここで $\gamma_j$ が $X$ の $j$ 期前の値に対する係数です。
2つのモデルの違いは、$X$ の過去の値を説明変数に含めるかどうかだけです。もし $X$ が $Y$ の予測にまったく貢献しないのであれば、$\gamma_1 = \gamma_2 = \cdots = \gamma_p = 0$ となるはずです。つまり、グレンジャー因果検定の帰無仮説と対立仮説は次のように書けます。
$$ \begin{align} H_0 &: \gamma_1 = \gamma_2 = \cdots = \gamma_p = 0 \quad (\text{$X$ は $Y$ にグレンジャー因果を持たない}) \\ H_1 &: \exists\, j \in \{1, \ldots, p\} \text{ such that } \gamma_j \neq 0 \quad (\text{$X$ は $Y$ にグレンジャー因果を持つ}) \end{align} $$
帰無仮説は「$X$ の過去の値がすべて無関係」と述べており、対立仮説は「少なくとも1つのラグで $X$ の過去の値が有意」と述べています。
VARモデルとの関係
上の非制限モデルは、実はVAR(p)モデルの $Y$ に関する方程式そのものです。VAR(p)モデルの完全な表現は以下の通りです。
$$ \begin{pmatrix} Y_t \\ X_t \end{pmatrix} = \begin{pmatrix} \beta_0 \\ \delta_0 \end{pmatrix} + \sum_{j=1}^{p} \begin{pmatrix} \beta_j & \gamma_j \\ \phi_j & \psi_j \end{pmatrix} \begin{pmatrix} Y_{t-j} \\ X_{t-j} \end{pmatrix} + \begin{pmatrix} \varepsilon_{Y,t} \\ \varepsilon_{X,t} \end{pmatrix} $$
この行列表現の1行目を取り出すと、先ほどの非制限モデルが得られます。VARモデルの係数行列の中で、$\gamma_j$（$X$ から $Y$ への影響を表す非対角要素）がすべてゼロかどうかを検定するのがグレンジャー因果検定です。
このようにVARモデルの構造を使って帰無仮説を立てたところで、次に「$\gamma_j$ がすべてゼロかどうか」を実際にどう検定するかを見ていきましょう。ここで登場するのがF検定です。
F検定統計量の導出
F検定の直感
F検定の考え方を直感的に理解しましょう。2つのモデルの残差二乗和を比較するのですが、単純に差を取るだけでは不十分です。なぜなら、パラメータの数が増えれば残差二乗和は必ず（少なくとも悪くならない方向に）減少するからです。これは、大学受験で必須科目と選択科目を増やせば合計点が上がりやすいのと同じ理屈です。
F検定は「パラメータを追加したことによるRSSの改善幅」を「追加したパラメータの数」で正規化し、さらに「非制限モデルの1パラメータあたりのRSS」で割ることで、モデルの改善が偶然を超えているかを判定します。
残差二乗和の定義
制限モデルの残差二乗和を $\text{RSS}_R$、非制限モデルの残差二乗和を $\text{RSS}_U$ とします。
$$ \text{RSS}_R = \sum_{t=p+1}^{T} \left(\hat{\varepsilon}_t^{(R)}\right)^2, \quad \text{RSS}_U = \sum_{t=p+1}^{T} \left(\hat{\varepsilon}_t^{(U)}\right)^2 $$
ここで $T$ はサンプルサイズ、$p$ はラグ次数です。$\hat{\varepsilon}_t^{(R)}$ と $\hat{\varepsilon}_t^{(U)}$ は、それぞれのモデルをOLSで推定したときの残差です。
非制限モデルは制限モデルの説明変数をすべて含んでいる（ネストされている）ため、$\text{RSS}_U \leq \text{RSS}_R$ が常に成り立ちます。問題は、この差が統計的に有意かどうかです。
F統計量の構成
制限モデルのパラメータ数（定数項を含む）は $p + 1$、非制限モデルのパラメータ数は $2p + 1$ です。追加されたパラメータの数は $q = p$ です。
F統計量は以下のように定義されます。
$$ F = \frac{(\text{RSS}_R – \text{RSS}_U) / q}{\text{RSS}_U / (T – 2p – 1)} $$
この式の意味を丁寧に読み解きましょう。分子の $(\text{RSS}_R – \text{RSS}_U) / q$ は、$X$ のラグ変数を $q = p$ 個追加したことによる1パラメータあたりのRSSの改善量です。
分母の $\text{RSS}_U / (T – 2p – 1)$ は、非制限モデルの1自由度あたりの残差分散の推定値（平均二乗誤差）です。$T – 2p – 1$ は非制限モデルの残差自由度で、サンプルサイズ $T$ から推定したパラメータ数 $2p + 1$ を引いたものです。
帰無仮説の下でのF統計量の分布
帰無仮説 $H_0: \gamma_1 = \cdots = \gamma_p = 0$ の下で、誤差項が正規分布に従うと仮定すると、F統計量はF分布に従います。
$$ F \sim F(q, \, T – 2p – 1) = F(p, \, T – 2p – 1) $$
第1自由度は追加パラメータ数 $q = p$、第2自由度は非制限モデルの残差自由度 $T – 2p – 1$ です。
F値が大きいほど、$X$ のラグ変数を追加したことによるモデルの改善が大きいことを意味します。有意水準 $\alpha$（例えば0.05）に対して、$F > F_{\alpha}(p, T – 2p – 1)$ であれば帰無仮説を棄却し、「$X$ は $Y$ にグレンジャー因果を持つ」と結論します。あるいは、p値を計算して $\alpha$ と比較する方法が一般的です。
ワルド検定との関係
実は、グレンジャー因果検定はF検定の代わりにワルド検定（カイ二乗検定）で実施することもできます。ワルド検定統計量 $W$ は以下で定義されます。
$$ W = T \cdot \frac{\text{RSS}_R – \text{RSS}_U}{\text{RSS}_U} $$
帰無仮説の下で、$W$ は漸近的にカイ二乗分布 $\chi^2(p)$ に従います。サンプルサイズが十分に大きいとき、$W = q \cdot F$ の関係があるため、F検定とワルド検定は本質的に同等です。statsmodelsの grangercausalitytests 関数では、F検定とカイ二乗検定の両方の結果が出力されます。
ここまでで2変量のケースにおけるF検定の仕組みが理解できました。しかし、現実の時系列データは3つ以上の変数を含むことが多いです。次に、多変量への拡張を見ていきましょう。
多変量への拡張
なぜ多変量が必要か
2変量のグレンジャー因果検定には大きな弱点があります。第3の変数（交絡因子）を考慮できないのです。
例えば、アイスクリームの売上 $X$ と水難事故の発生件数 $Y$ の2変量だけを見ると、$X$ が $Y$ にグレンジャー因果を持つように見えるかもしれません。しかし、真の原因は「気温」$Z$ であり、$Z$ を考慮すれば $X \to Y$ の見かけの因果関係は消えるはずです。
多変量のグレンジャー因果検定では、他の変数の影響を「コントロール」した上で、特定の変数ペア間の予測的因果を検定します。
多変量VARモデルにおける定式化
$K$ 個の定常時系列 $\bm{z}_t = (z_{1,t}, z_{2,t}, \ldots, z_{K,t})^\top$ を考えます。VAR(p)モデルは以下のように書けます。
$$ \bm{z}_t = \bm{c} + \bm{A}_1 \bm{z}_{t-1} + \bm{A}_2 \bm{z}_{t-2} + \cdots + \bm{A}_p \bm{z}_{t-p} + \bm{u}_t $$
ここで $\bm{c}$ は $K \times 1$ の定数ベクトル、$\bm{A}_j$ は $K \times K$ の係数行列、$\bm{u}_t$ は $K \times 1$ の誤差ベクトルで $\bm{u}_t \sim \mathcal{N}(\bm{0}, \bm{\Sigma}_u)$ です。
変数 $z_m$ が変数 $z_k$ にグレンジャー因果を持つかを検定するには、$\bm{A}_j$ の $(k, m)$ 成分を調べます。帰無仮説は以下の通りです。
$$ H_0: [\bm{A}_1]_{k,m} = [\bm{A}_2]_{k,m} = \cdots = [\bm{A}_p]_{k,m} = 0 $$
つまり、すべてのラグにおいて、$z_m$ から $z_k$ への影響係数がゼロであることを検定します。
この場合の制限モデルは、$z_k$ の方程式から $z_m$ の全ラグ変数を除外したものです。非制限モデルはVARモデルの $z_k$ に関する方程式そのものです。F統計量の構成は2変量の場合と同じですが、非制限モデルの自由度が変わります。
非制限モデルのパラメータ数は $Kp + 1$（$K$ 変数の各 $p$ ラグ + 定数項）であるため、残差自由度は $T – Kp – 1$ となります。追加パラメータ数は依然として $q = p$ です。
$$ F = \frac{(\text{RSS}_R – \text{RSS}_U) / p}{\text{RSS}_U / (T – Kp – 1)} $$
多変量のグレンジャー因果検定では、$z_m$ 以外の変数（$z_1, \ldots, z_{m-1}, z_{m+1}, \ldots, z_K$）の影響をすべて制御した上で $z_m$ の予測的因果を検定するため、2変量の検定よりも交絡因子の問題に対してロバストです。ただし、モデルに含まれていない交絡因子の影響は依然として除去できません。
多変量への拡張を理解したところで、次に実務上極めて重要な問題に移りましょう。ラグ次数 $p$ をどう決めるかという問題です。
ラグ次数の選択
ラグ次数が結果に与える影響
グレンジャー因果検定の結果は、ラグ次数 $p$ の選択に強く依存します。これは検定の信頼性に直結する問題です。
ラグ次数が小さすぎると、重要なラグの影響を見落としてしまい、本来存在するグレンジャー因果を検出できない（第2種の過誤）可能性が高まります。例えば、金融政策の効果が6か月後に表れるのにラグを3か月に設定すれば、因果関係を見逃してしまいます。
一方、ラグ次数が大きすぎると、推定すべきパラメータが増えて各推定値の精度が悪化し、検出力が低下します。さらに、データが限られている場合はモデルが過適合し、サンプル外での予測精度が悪化します。
情報量規準によるラグ選択
最も一般的なアプローチは、情報量規準（Information Criterion） を用いてVARモデルの最適なラグ次数を選択する方法です。代表的な規準は以下の2つです。
赤池情報量規準（AIC）:
$$ \text{AIC}(p) = \ln |\hat{\bm{\Sigma}}_u(p)| + \frac{2}{T} \cdot K^2 p $$
ベイズ情報量規準（BIC / SIC）:
$$ \text{BIC}(p) = \ln |\hat{\bm{\Sigma}}_u(p)| + \frac{\ln T}{T} \cdot K^2 p $$
ここで $\hat{\bm{\Sigma}}_u(p)$ はラグ次数 $p$ のVARモデルから推定した残差の共分散行列、$|\cdot|$ はその行列式です。$K$ は変数の数、$T$ はサンプルサイズです。
第1項 $\ln |\hat{\bm{\Sigma}}_u(p)|$ はモデルの当てはまりの良さを表し、ラグを増やすほど小さくなります。第2項はパラメータ数に対するペナルティで、ラグを増やすほど大きくなります。この2項のバランスが最適なラグを決定します。
BICはAICよりもペナルティが強い（$\ln T > 2$ for $T \geq 8$）ため、一般にBICはより少ないラグを選びます。実務的には、AICとBICの両方を計算し、結果が大きく異なる場合はロバスト性を確認するために複数のラグ次数で検定を行うことが推奨されます。
実務的な手順
グレンジャー因果検定における典型的なラグ選択の手順をまとめます。
- $p = 1, 2, \ldots, p_{\max}$ の各ラグ次数でVARモデルを推定する
- AICまたはBICを計算し、最小値を与える $p^*$ を選ぶ
- 選ばれた $p^*$ でグレンジャー因果検定を実施する
- ロバスト性チェックとして、$p^* \pm 1$ や $p^* \pm 2$ でも検定を行い、結論が変わらないか確認する
ラグ次数の上限 $p_{\max}$ は、データの長さに応じて決めます。目安としては $p_{\max} = \lfloor 12(T/100)^{1/4} \rfloor$（Schwert, 1989）や、季節データなら周期（12か月、4四半期）の2倍程度がよく使われます。
ここまでで、グレンジャー因果検定の理論的な枠組みがすべて揃いました。次はいよいよPythonで実装していきましょう。まずは理論の理解を確認するために、スクラッチで実装します。
Pythonでのスクラッチ実装
実装方針
ここでは、2変量のグレンジャー因果検定をNumPyだけで実装します。ライブラリに頼らず自分で実装することで、F統計量の計算過程が具体的に理解できます。実装は以下のステップで進めます。
- ラグ行列の作成
- 制限モデルと非制限モデルのOLS推定
- 残差二乗和の計算
- F統計量とp値の算出
まず、ラグ行列を作成するヘルパー関数を定義します。
import numpy as np
from scipy import stats
def create_lag_matrix(data, max_lag):
"""時系列データからラグ行列を作成する"""
T = len(data)
n_cols = data.shape[1] if data.ndim > 1 else 1
if data.ndim == 1:
data = data.reshape(-1, 1)
# 各ラグの値を横に並べる
lag_matrices = []
for lag in range(1, max_lag + 1):
lag_matrices.append(data[max_lag - lag : T - lag, :])
# 定数項(1の列)を先頭に追加
n_samples = T - max_lag
X = np.column_stack([np.ones(n_samples)] + lag_matrices)
y = data[max_lag:, :]
return X, y
このコードは、時系列データを受け取り、各ラグの値を列方向に並べた行列を返します。定数項（切片）も含めています。例えばラグ2の場合、$t$ 時点の行には $[1, z_{t-1}, z_{t-2}]$ が並びます。
次に、グレンジャー因果検定本体を実装します。
def granger_causality_test_scratch(x, y, max_lag):
"""
xがyにグレンジャー因果を持つかを検定する（スクラッチ実装）
Parameters
----------
x : array-like, shape (T,)
原因系列（の候補）
y : array-like, shape (T,)
結果系列
max_lag : int
ラグ次数
Returns
-------
dict : F統計量、p値、各モデルのRSS
"""
x = np.asarray(x, dtype=float)
y = np.asarray(y, dtype=float)
T = len(y)
# --- 制限モデル: y ~ y のラグのみ ---
X_r, y_r = create_lag_matrix(y, max_lag)
# X_r の形状: (T-p, p+1) [定数項 + yのpラグ]
# OLS推定: beta = (X'X)^{-1} X'y
beta_r = np.linalg.lstsq(X_r, y_r, rcond=None)[0]
residuals_r = y_r - X_r @ beta_r
rss_r = float(np.sum(residuals_r ** 2))
# --- 非制限モデル: y ~ y のラグ + x のラグ ---
data_both = np.column_stack([y, x])
X_u_full, y_u_full = create_lag_matrix(data_both, max_lag)
# X_u_full の形状: (T-p, 2p+1) [定数項 + yのpラグ + xのpラグ]
# y_u_full の列0 がyの目的変数
y_u = y_u_full[:, 0:1]
X_u = X_u_full # 全列を説明変数に使う
beta_u = np.linalg.lstsq(X_u, y_u, rcond=None)[0]
residuals_u = y_u - X_u @ beta_u
rss_u = float(np.sum(residuals_u ** 2))
# --- F統計量の計算 ---
n = len(y_u) # 有効サンプルサイズ (T - p)
q = max_lag # 追加パラメータ数
df_u = n - 2 * max_lag - 1 # 非制限モデルの残差自由度
f_stat = ((rss_r - rss_u) / q) / (rss_u / df_u)
p_value = 1.0 - stats.f.cdf(f_stat, q, df_u)
return {
"f_stat": f_stat,
"p_value": p_value,
"rss_restricted": rss_r,
"rss_unrestricted": rss_u,
"df_num": q,
"df_denom": df_u,
"n_obs": n
}
この関数では、まず制限モデル（$Y$ の自己ラグのみ）と非制限モデル（$Y$ の自己ラグ + $X$ のラグ）をOLSで推定し、それぞれの残差二乗和を求めています。F統計量は先ほど導出した式をそのまま計算しています。p値はF分布の上側確率として求めます。
np.linalg.lstsq を使っているのは、$(\bm{X}^\top \bm{X})^{-1}$ を直接計算するよりも数値的に安定だからです。正規方程式 $\bm{X}^\top \bm{X} \hat{\bm{\beta}} = \bm{X}^\top \bm{y}$ を解くのと数学的には等価ですが、特異行列に近い場合でもロバストに動作します。
それでは、スクラッチ実装の動作確認をしましょう。まず、因果関係が存在するデータで検定が正しく棄却できるか確認します。
np.random.seed(42)
T = 500
# x -> y の因果関係を持つデータを生成
x = np.zeros(T)
y = np.zeros(T)
for t in range(2, T):
x[t] = 0.5 * x[t-1] + np.random.normal(0, 1)
y[t] = 0.3 * y[t-1] + 0.6 * x[t-1] + np.random.normal(0, 1)
# x -> y のグレンジャー因果検定
result_xy = granger_causality_test_scratch(x, y, max_lag=2)
print("=== x -> y の検定 ===")
print(f"F統計量: {result_xy['f_stat']:.4f}")
print(f"p値: {result_xy['p_value']:.6f}")
print(f"RSS(制限): {result_xy['rss_restricted']:.2f}")
print(f"RSS(非制限): {result_xy['rss_unrestricted']:.2f}")
# y -> x のグレンジャー因果検定（こちらは棄却されないはず）
result_yx = granger_causality_test_scratch(y, x, max_lag=2)
print("\n=== y -> x の検定 ===")
print(f"F統計量: {result_yx['f_stat']:.4f}")
print(f"p値: {result_yx['p_value']:.6f}")
上のコードでは、$x_t = 0.5 x_{t-1} + \varepsilon_{x,t}$、$y_t = 0.3 y_{t-1} + 0.6 x_{t-1} + \varepsilon_{y,t}$ というデータ生成過程を設定しています。設計上、$x \to y$ の因果は存在しますが、$y \to x$ の因果は存在しません。
実行すると、$x \to y$ の検定ではF統計量が非常に大きくなり、p値はほぼ0になります。これは $x$ の過去の値が $y$ の予測に強く貢献していることを正しく検出しています。一方、$y \to x$ の検定ではF統計量が小さく、p値は有意水準0.05を大きく上回ります。$y$ の過去の値は $x$ の予測に貢献しておらず、検定は正しく棄却を回避しています。制限モデルと非制限モデルのRSSの差に注目すると、$x \to y$ の方では非制限モデルのRSSが大幅に小さくなっているのに対し、$y \to x$ の方では両モデルのRSSがほぼ同じであることが確認できます。
スクラッチ実装で検定の仕組みが理解できたところで、次に実務で使うことが多いstatsmodelsのAPIを用いた実装方法を見ていきましょう。
statsmodelsでの実装
grangercausalitytests の使い方
statsmodelsには grangercausalitytests 関数が用意されており、わずか数行でグレンジャー因果検定を実行できます。この関数は内部でF検定とカイ二乗検定の両方を計算し、複数のラグ次数についてまとめて結果を出力します。
import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.stattools import grangercausalitytests
np.random.seed(42)
T = 500
# 先ほどと同じデータを生成
x = np.zeros(T)
y = np.zeros(T)
for t in range(2, T):
x[t] = 0.5 * x[t-1] + np.random.normal(0, 1)
y[t] = 0.3 * y[t-1] + 0.6 * x[t-1] + np.random.normal(0, 1)
# statsmodelsのgrangercausalitytests
# 注意: 第1列が結果変数(y), 第2列が原因変数(x)
data = np.column_stack([y, x])
print("=== statsmodels: x -> y の検定 ===")
results = grangercausalitytests(data, maxlag=4, verbose=True)
grangercausalitytests を使う際に注意すべき点がいくつかあります。まず、入力データの列の順序です。第1列が結果変数（$Y$）、第2列が原因変数（$X$）です。「$X$ が $Y$ にグレンジャー因果を持つか」を検定する場合、直感とは逆に $Y$ を先に置く必要があります。
また、verbose=True にすると各ラグについてF検定とカイ二乗検定の結果が表示されます。出力される4つの検定は、ssr_ftest（F検定）、ssr_chi2test（カイ二乗検定）、lrtest（尤度比検定）、params_ftest（パラメータF検定）です。通常はssr_ftestの結果を使います。
実行すると、ラグ1から4のすべてにおいてp値が非常に小さく（0.05を大幅に下回り）、帰無仮説が棄却されます。これはスクラッチ実装の結果と整合しており、$x$ が $y$ に対して強い予測的因果を持つことを示しています。ラグ1のF統計量が最も大きいのは、データ生成過程で $x_{t-1}$ が $y_t$ に直接影響しているためです。
VARモデル経由でのグレンジャー因果検定
statsmodelsのVARモデルクラスにも test_causality メソッドが実装されています。こちらは多変量の場合に便利です。
from statsmodels.tsa.api import VAR
# VARモデルの推定
data_var = np.column_stack([y, x])
model = VAR(data_var, names=["y", "x"])
# AICで最適ラグを選択
lag_order_results = model.select_order(maxlags=8)
print("=== ラグ次数の選択 ===")
print(f"AIC: {lag_order_results.aic}")
print(f"BIC: {lag_order_results.bic}")
print(f"選択されたラグ (AIC): {lag_order_results.aic}")
print(f"選択されたラグ (BIC): {lag_order_results.bic}")
print(lag_order_results.summary())
# 最適ラグでVARを推定
fitted = model.fit(maxlags=lag_order_results.aic, ic=None)
# グレンジャー因果検定
print("\n=== VARモデルによる因果検定: x -> y ===")
gc_result = fitted.test_causality("y", causing="x", kind="f")
print(gc_result.summary())
print("\n=== VARモデルによる因果検定: y -> x ===")
gc_result_rev = fitted.test_causality("x", causing="y", kind="f")
print(gc_result_rev.summary())
VARモデル経由で検定する利点は、最適ラグの選択（select_order）と因果検定（test_causality）がシームレスにつながることです。select_order はAIC、BIC、HQIC、FPEの各規準を一覧表示するため、ラグの選択根拠を透明に示すことができます。
出力を見ると、AICおよびBICで選択されたラグ次数でVARモデルを推定した上で、$x \to y$ の因果が有意に検出され、$y \to x$ の因果は検出されないことが確認できます。VARモデルの test_causality メソッドは、ワルド検定（デフォルト）またはF検定（kind="f"）を選択できます。
2つの実装方法（grangercausalitytests とVARモデル経由）の使い分けとして、2変量の簡単な分析には前者が手軽で、3変数以上の多変量分析や他の分析（インパルス応答など）と組み合わせたい場合には後者が適しています。
理論と実装の基本が揃ったところで、次はより本格的な合成データを使って検定の性質を実験的に確認しましょう。
実験: 合成データでの検証
実験1: 検出力の評価
因果の強さ（係数の大きさ）と検定の検出力（帰無仮説を正しく棄却できる確率）の関係を調べましょう。直感的には、因果が強いほど検出しやすく、弱いほど見逃しやすいはずです。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
def granger_f_test(x, y, max_lag):
"""簡略版グレンジャーF検定（p値を返す）"""
T = len(y)
p = max_lag
# 制限モデルのラグ行列
Y_lags_r = np.column_stack([
np.ones(T - p),
*[y[p - j : T - j] for j in range(1, p + 1)]
])
y_target = y[p:]
beta_r = np.linalg.lstsq(Y_lags_r, y_target, rcond=None)[0]
rss_r = np.sum((y_target - Y_lags_r @ beta_r) ** 2)
# 非制限モデルのラグ行列
Y_lags_u = np.column_stack([
Y_lags_r,
*[x[p - j : T - j] for j in range(1, p + 1)]
])
beta_u = np.linalg.lstsq(Y_lags_u, y_target, rcond=None)[0]
rss_u = np.sum((y_target - Y_lags_u @ beta_u) ** 2)
n = T - p
df_denom = n - 2 * p - 1
f_stat = ((rss_r - rss_u) / p) / (rss_u / df_denom)
p_value = 1.0 - stats.f.cdf(f_stat, p, df_denom)
return p_value
# 因果の強さを変化させて検出力を測定
np.random.seed(123)
causal_strengths = np.arange(0.0, 0.81, 0.05)
n_simulations = 500
T = 300
significance_level = 0.05
power_curve = []
for gamma in causal_strengths:
rejections = 0
for _ in range(n_simulations):
x = np.zeros(T)
y = np.zeros(T)
for t in range(2, T):
x[t] = 0.5 * x[t-1] + np.random.normal(0, 1)
y[t] = 0.4 * y[t-1] + gamma * x[t-1] + np.random.normal(0, 1)
p_val = granger_f_test(x, y, max_lag=2)
if p_val < significance_level:
rejections += 1
power_curve.append(rejections / n_simulations)
plt.figure(figsize=(9, 5))
plt.plot(causal_strengths, power_curve, "o-", color="#00bcd4", linewidth=2, markersize=5)
plt.axhline(y=significance_level, color="#ff9800", linestyle="--", label=f"有意水準 α = {significance_level}")
plt.axhline(y=0.80, color="#4caf50", linestyle=":", alpha=0.7, label="検出力 80%")
plt.xlabel("Causal strength (γ)", fontsize=12)
plt.ylabel("Power (rejection rate)", fontsize=12)
plt.title("Granger Causality Test: Power vs Causal Strength", fontsize=13)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.ylim(-0.02, 1.02)
plt.tight_layout()
plt.show()
このグラフは因果の強さ $\gamma$ と検出力の関係を示しています。$\gamma = 0$（因果なし）のとき、棄却率は有意水準 $\alpha = 0.05$ にほぼ一致します。これは検定の第1種の過誤率が正しく制御されていることを意味します。$\gamma$ が増加するにつれて検出力は急速に上昇し、$\gamma \approx 0.2$ で検出力80%に達します。$\gamma \geq 0.4$ ではほぼ確実（検出力99%以上）に因果を検出できます。この結果は、適切なサンプルサイズ（$T=300$）があれば、中程度以上の因果関係は十分に検出可能であることを示しています。
実験2: サンプルサイズの影響
次に、サンプルサイズが検出力に与える影響を調べます。
np.random.seed(456)
sample_sizes = [50, 100, 200, 300, 500, 1000]
gamma_fixed = 0.3 # 中程度の因果
n_simulations = 500
power_by_T = []
for T in sample_sizes:
rejections = 0
for _ in range(n_simulations):
x = np.zeros(T)
y = np.zeros(T)
for t in range(2, T):
x[t] = 0.5 * x[t-1] + np.random.normal(0, 1)
y[t] = 0.4 * y[t-1] + gamma_fixed * x[t-1] + np.random.normal(0, 1)
p_val = granger_f_test(x, y, max_lag=2)
if p_val < 0.05:
rejections += 1
power_by_T.append(rejections / n_simulations)
plt.figure(figsize=(9, 5))
plt.plot(sample_sizes, power_by_T, "s-", color="#e91e63", linewidth=2, markersize=7)
plt.axhline(y=0.80, color="#4caf50", linestyle=":", alpha=0.7, label="検出力 80%")
plt.xlabel("Sample size (T)", fontsize=12)
plt.ylabel("Power (rejection rate)", fontsize=12)
plt.title(f"Granger Causality Test: Power vs Sample Size (γ = {gamma_fixed})", fontsize=13)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.ylim(-0.02, 1.02)
plt.tight_layout()
plt.show()
サンプルサイズが50の場合、検出力は50%程度にとどまり、因果関係を見逃す可能性が半分もあります。サンプルサイズを200に増やすと検出力は80%を超え、実用的に十分な水準になります。1000では検出力はほぼ100%です。この実験から、グレンジャー因果検定で信頼できる結果を得るためには、最低でも200以上のサンプルサイズが望ましいことがわかります。もちろん、因果が弱い場合にはさらに多くのデータが必要です。
実験3: 双方向因果と見せかけの因果
最後に、2つの重要なシナリオを検証します。双方向因果（フィードバック）のケースと、共通因子による見せかけの因果（疑似因果）のケースです。
np.random.seed(789)
T = 500
# --- シナリオA: 双方向因果 ---
x_bi = np.zeros(T)
y_bi = np.zeros(T)
for t in range(2, T):
x_bi[t] = 0.3 * x_bi[t-1] + 0.4 * y_bi[t-1] + np.random.normal(0, 1)
y_bi[t] = 0.3 * y_bi[t-1] + 0.5 * x_bi[t-1] + np.random.normal(0, 1)
p_xy_bi = granger_f_test(x_bi, y_bi, max_lag=2)
p_yx_bi = granger_f_test(y_bi, x_bi, max_lag=2)
print("=== シナリオA: 双方向因果 ===")
print(f"x -> y: p値 = {p_xy_bi:.6f} {'***' if p_xy_bi < 0.001 else '**' if p_xy_bi < 0.01 else '*' if p_xy_bi < 0.05 else 'ns'}")
print(f"y -> x: p値 = {p_yx_bi:.6f} {'***' if p_yx_bi < 0.001 else '**' if p_yx_bi < 0.01 else '*' if p_yx_bi < 0.05 else 'ns'}")
# --- シナリオB: 共通因子による見せかけの因果 ---
z = np.zeros(T) # 共通因子
x_sp = np.zeros(T)
y_sp = np.zeros(T)
for t in range(2, T):
z[t] = 0.7 * z[t-1] + np.random.normal(0, 1)
x_sp[t] = 0.3 * x_sp[t-1] + 0.8 * z[t-1] + np.random.normal(0, 0.5)
y_sp[t] = 0.3 * y_sp[t-1] + 0.8 * z[t-2] + np.random.normal(0, 0.5)
# 2変量検定（zを無視）
p_xy_sp = granger_f_test(x_sp, y_sp, max_lag=3)
p_yx_sp = granger_f_test(y_sp, x_sp, max_lag=3)
print("\n=== シナリオB: 共通因子による見せかけの因果（2変量検定）===")
print(f"x -> y: p値 = {p_xy_sp:.6f} {'***' if p_xy_sp < 0.001 else '**' if p_xy_sp < 0.01 else '*' if p_xy_sp < 0.05 else 'ns'}")
print(f"y -> x: p値 = {p_yx_sp:.6f} {'***' if p_yx_sp < 0.001 else '**' if p_yx_sp < 0.01 else '*' if p_yx_sp < 0.05 else 'ns'}")
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].plot(x_bi[:200], label="x", alpha=0.8, color="#00bcd4")
axes[0].plot(y_bi[:200], label="y", alpha=0.8, color="#ff9800")
axes[0].set_title("Scenario A: Bidirectional Causality", fontsize=12)
axes[0].set_xlabel("Time")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].plot(x_sp[:200], label="x", alpha=0.8, color="#00bcd4")
axes[1].plot(y_sp[:200], label="y", alpha=0.8, color="#ff9800")
axes[1].plot(z[:200], label="z (hidden)", alpha=0.5, linestyle="--", color="#9c27b0")
axes[1].set_title("Scenario B: Spurious Causality (common factor z)", fontsize=12)
axes[1].set_xlabel("Time")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
シナリオAでは、$x$ と $y$ が互いに影響し合うフィードバックシステムを構成しています。実行結果を見ると、$x \to y$ と $y \to x$ の両方向で帰無仮説が棄却されます。これは設計通りです。双方向因果はグレンジャー因果検定で適切に検出できることが確認できました。
シナリオBがより興味深いケースです。$x$ と $y$ は直接的な因果関係を持たず、共通因子 $z$ が $x$ にはラグ1で、$y$ にはラグ2で影響しています。つまり $z$ の影響が $x$ に先に現れるため、$z$ を無視した2変量検定では「$x$ が $y$ を予測する」ように見えてしまいます。実行結果では、$x \to y$ のp値が有意になり、見せかけのグレンジャー因果が検出されます。グラフを見ても、$x$ と $y$ が類似した動きをしていますが、隠れた因子 $z$（紫の破線）が両者を駆動していることがわかります。このような見せかけの因果を回避するためには、$z$ を含めた多変量の検定が必要です。
この実験結果は、次のセクションで議論するグレンジャー因果の限界と直結しています。
グレンジャー因果の限界と注意点
グレンジャー因果検定は強力なツールですが、いくつかの重要な限界があります。これらを理解しないまま結果を鵜呑みにすると、誤った結論を導いてしまいます。
1. 定常性の仮定
グレンジャー因果検定は定常な時系列を前提としています。トレンドや単位根を持つ非定常時系列にそのまま適用すると、見せかけの回帰（spurious regression）が発生し、本来存在しない因果関係が検出されてしまいます。
例えば、2つの独立なランダムウォーク $X_t = X_{t-1} + \varepsilon_{X,t}$、$Y_t = Y_{t-1} + \varepsilon_{Y,t}$ に対してグレンジャー因果検定を行うと、高い確率で有意な結果が得られます。これは両系列がたまたま同じ方向にドリフトしているだけで、真の因果関係はありません。
対策: 検定の前にADF検定（拡張ディッキー-フラー検定）やKPSS検定で定常性を確認し、非定常であれば差分を取るか、共和分検定（Johansenの検定）を先に行います。共和分関係がある場合は、ベクトル誤差修正モデル（VECM）の枠組みでグレンジャー因果を検定します。
2. 省略変数バイアス（交絡因子）
先ほどの実験3のシナリオBで見たように、モデルに含まれていない第3の変数が $X$ と $Y$ の両方に影響している場合、見せかけのグレンジャー因果が検出されます。
これはグレンジャー因果検定に限らず、観測データに基づく因果推論全般の問題です。しかし、グレンジャー因果検定は「因果」という名前がついているために、この限界が見落とされやすいです。
対策: 分析に影響を与える可能性のある変数をできるだけ多くモデルに含め、多変量の検定を行います。ただし、すべての交絡因子を特定することは現実には困難であるため、結果の解釈は常に慎重に行う必要があります。
3. ラグ次数の選択への依存
ラグ次数の選択が検定結果に大きく影響することは先に述べました。最適なラグ次数は情報量規準で選べますが、AICとBICで異なるラグが選択されることがあります。
対策: 複数の情報量規準でラグを選び、さらにラグ次数を前後に変えても結果が安定するか（ロバスト性）を確認します。結果がラグ次数に敏感な場合は、その因果関係の証拠は弱いと判断すべきです。
4. 瞬時因果と同時刻の相関
標準的なグレンジャー因果検定は、$X$ の過去の値が $Y$ の予測に役立つかを検定します。$X_t$ と $Y_t$ の同時刻の相関は考慮しません。そのため、データのサンプリング間隔が因果の伝搬速度に比べて粗い場合、因果関係を見逃す可能性があります。
例えば、月次データでは1日以内に伝搬する因果関係は検出できません。秒単位の高頻度データでは検出可能でも、日次データに集約すると見えなくなることがあります。
対策: データのサンプリング頻度が分析対象のダイナミクスに対して十分に細かいかを事前に検討します。必要に応じてVARモデルの瞬時因果（instantaneous causality）の検定を追加で行います。
5. 非線形因果の見落とし
グレンジャー因果検定は線形モデル（VARモデル）に基づいているため、$X$ と $Y$ の間に非線形の因果関係がある場合は検出できません。例えば、$Y_t = f(X_{t-1}) + \varepsilon_t$ で $f$ が強い非線形関数のとき、線形モデルでは $X$ の寄与が過小評価されます。
対策: 非線形グレンジャー因果検定（Hiemstra-Jones検定、Transfer Entropy、ニューラルネットワークベースの方法など）を検討します。ただし、これらの方法はより多くのデータを必要とし、解釈も複雑になります。
6. 多重検定の問題
多変量システムで全ペアのグレンジャー因果を検定すると、$K$ 変数で $K(K-1)$ 回の検定を行うことになります。有意水準5%で検定すれば、因果がなくても約5%は偽陽性になります。10変数なら90回の検定で約4〜5個の偽陽性が期待されます。
対策: ボンフェローニ補正や偽発見率（FDR）制御（Benjamini-Hochberg法）を適用して、多重比較の問題に対処します。
これらの限界を踏まえた上で、グレンジャー因果検定を正しく使いこなすためのポイントをまとめましょう。
実践的なワークフロー
グレンジャー因果検定を実務で使う際の標準的な手順を、コードとともに整理します。以下のワークフローは、ここまでの議論を踏まえたベストプラクティスです。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller, grangercausalitytests
from statsmodels.tsa.api import VAR
def granger_workflow(y, x, names=("y", "x"), max_lag_search=12, alpha=0.05):
"""
グレンジャー因果検定の実践的ワークフロー
Step 1: 定常性の確認（ADF検定）
Step 2: 最適ラグの選択（AIC/BIC）
Step 3: グレンジャー因果検定
Step 4: ロバスト性の確認
"""
print("=" * 60)
print("グレンジャー因果検定ワークフロー")
print("=" * 60)
# Step 1: 定常性の確認
print("\n--- Step 1: 定常性の確認 (ADF検定) ---")
for name, series in zip(names, [y, x]):
adf_result = adfuller(series, autolag="AIC")
stationary = "定常" if adf_result[1] < alpha else "非定常 (差分が必要)"
print(f" {name}: ADF統計量 = {adf_result[0]:.4f}, "
f"p値 = {adf_result[1]:.4f} → {stationary}")
# Step 2: 最適ラグの選択
print(f"\n--- Step 2: 最適ラグの選択 (maxlag={max_lag_search}) ---")
data = np.column_stack([y, x])
model = VAR(data, names=list(names))
lag_results = model.select_order(maxlags=max_lag_search)
aic_lag = lag_results.aic
bic_lag = lag_results.bic
print(f" AIC最適ラグ: {aic_lag}")
print(f" BIC最適ラグ: {bic_lag}")
# Step 3: グレンジャー因果検定
selected_lag = aic_lag
print(f"\n--- Step 3: グレンジャー因果検定 (ラグ = {selected_lag}) ---")
print(f"\n {names[1]} → {names[0]}:")
gc_result = grangercausalitytests(data, maxlag=selected_lag, verbose=False)
f_stat = gc_result[selected_lag][0]["ssr_ftest"][0]
p_value = gc_result[selected_lag][0]["ssr_ftest"][1]
sig = "***" if p_value < 0.001 else "**" if p_value < 0.01 else "*" if p_value < alpha else "ns"
print(f" F統計量 = {f_stat:.4f}, p値 = {p_value:.6f} {sig}")
data_rev = np.column_stack([x, y])
print(f"\n {names[0]} → {names[1]}:")
gc_result_rev = grangercausalitytests(data_rev, maxlag=selected_lag, verbose=False)
f_stat_rev = gc_result_rev[selected_lag][0]["ssr_ftest"][0]
p_value_rev = gc_result_rev[selected_lag][0]["ssr_ftest"][1]
sig_rev = "***" if p_value_rev < 0.001 else "**" if p_value_rev < 0.01 else "*" if p_value_rev < alpha else "ns"
print(f" F統計量 = {f_stat_rev:.4f}, p値 = {p_value_rev:.6f} {sig_rev}")
# Step 4: ロバスト性の確認
print(f"\n--- Step 4: ロバスト性の確認 ---")
print(f" {names[1]} → {names[0]} の検定結果 (各ラグ):")
lag_min = max(1, selected_lag - 2)
lag_max = selected_lag + 2
gc_robust = grangercausalitytests(data, maxlag=lag_max, verbose=False)
for lag in range(lag_min, lag_max + 1):
f_val = gc_robust[lag][0]["ssr_ftest"][0]
p_val = gc_robust[lag][0]["ssr_ftest"][1]
marker = " ← selected" if lag == selected_lag else ""
print(f" ラグ {lag}: F = {f_val:.4f}, p = {p_val:.6f}{marker}")
return {"selected_lag": selected_lag, "p_value_forward": p_value, "p_value_reverse": p_value_rev}
このワークフロー関数は、定常性の確認から検定の実施、ロバスト性の確認までを一貫して行います。各ステップの出力を確認しながら、信頼できる結論を導くことができます。
以下のコードで動作確認をしてみましょう。
# 使用例: 因果関係を持つ合成データ
np.random.seed(42)
T = 500
x_demo = np.zeros(T)
y_demo = np.zeros(T)
for t in range(2, T):
x_demo[t] = 0.5 * x_demo[t-1] + np.random.normal(0, 1)
y_demo[t] = 0.3 * y_demo[t-1] + 0.4 * x_demo[t-1] + np.random.normal(0, 1)
results = granger_workflow(y_demo, x_demo, names=("y", "x"))
ワークフロー関数の出力から、以下の情報が一目で確認できます。まず、ADF検定で両系列が定常であること。次に、AICとBICで選択された最適ラグ。そして、双方向の検定結果（$x \to y$ は有意、$y \to x$ は非有意）。最後に、前後のラグでも結果が安定していること（ロバスト性）。このように体系的に分析を進めることで、見落としのない信頼性の高い検定が実現します。
まとめ
本記事では、グレンジャー因果検定の理論と実装について解説しました。
- グレンジャー因果の概念: 真の因果関係ではなく「予測的因果」を検定する手法です。「$X$ の過去の情報が $Y$ の予測精度を有意に改善するか」を、制限モデルと非制限モデルの比較で判定します
- 数学的定式化: VARモデルのフレームワークにおいて、非制限モデル（$X$ のラグを含む）と制限モデル（$X$ のラグを含まない）の残差二乗和をF検定で比較します。帰無仮説は「$X$ のラグ係数がすべてゼロ」です
- F統計量の導出: $F = [(\text{RSS}_R – \text{RSS}_U)/p] / [\text{RSS}_U/(T-2p-1)]$ で計算し、帰無仮説の下で $F(p, T-2p-1)$ に従います
- 多変量への拡張: 他の変数をコントロールした上で因果を検定でき、交絡因子の影響を軽減できます
- ラグ次数の選択: AIC・BICなどの情報量規準で選び、ロバスト性の確認を行うことが重要です
- 限界と注意点: 定常性の仮定、省略変数バイアス、非線形因果の見落とし、多重検定問題など、結果の解釈には慎重さが求められます
グレンジャー因果検定は、時系列データから変数間の予測的依存関係を抽出するための第一歩です。しかし、「グレンジャー因果がある」ことは「真の因果関係がある」ことを意味しません。この検定の結果は、より精密な因果分析（構造VAR、局所射影法、自然実験など）や、ドメイン知識に基づく因果的推論の出発点として位置づけるのが適切です。
次のステップとして、以下の記事も参考にしてください。
- VARモデル(ベクトル自己回帰)の理論と実装 — グレンジャー因果検定の基盤であるVARモデルのインパルス応答分析や分散分解も合わせて学ぶことで、時系列間の動的関係をより深く理解できます
- 【時系列解析】ARモデルを分かりやすく解説して実装する — 単変量の自己回帰モデルの理論と実装を振り返ることで、制限モデルの構造をより明確に理解できます