前回の記事でグレンジャー因果性という用語を紹介したので、今回はPythonで実際に計算してみましょう。
参考: グレンジャー因果性
データは、VARモデルのパラメーター推定を行った時のように、乱数で生成します。
参考: Python(statsmodels)でVARモデルのパラメーター推定
コードの重複する部分は省略して、変数resultに、推定したモデルが格納されているとします。
(今回の記事用に実行し直したので、パラメーターは前回と違います。)
この推定結果のオブジェクトが、 test_causality というメソッドを持っており、
それを使って、グレンジャー因果性の検証ができます。
ドキュメント: statsmodels.tsa.vector_ar.var_model.VARResults.test_causality
statsmodelsにはこの他にも、二つの過程のグレンジャー因果性を確認するツールとして、
grangercausalitytests というのが用意されているようなのですが、今回はこちらは使いません。
それでは早速やってみましょう。
# 推定結果のVARモデルの確認。(以前の記事のコードで生成されている前提)
print(result.summary())
"""
Summary of Regression Results
==================================
Model: VAR
Method: OLS
Date: Sun, 09, Feb, 2020
Time: 16:37:19
--------------------------------------------------------------------
No. of Equations: 2.00000 BIC: 1.41978
Nobs: 99.0000 HQIC: 1.32614
Log likelihood: -337.444 FPE: 3.53438
AIC: 1.26250 Det(Omega_mle): 3.32953
--------------------------------------------------------------------
Results for equation y1
========================================================================
coefficient std. error t-stat prob
------------------------------------------------------------------------
const 1.263791 0.705621 1.791 0.073
L1.y1 0.604660 0.124431 4.859 0.000
L1.y2 0.267836 0.145104 1.846 0.065
========================================================================
Results for equation y2
========================================================================
coefficient std. error t-stat prob
------------------------------------------------------------------------
const -3.289405 0.292191 -11.258 0.000
L1.y1 0.601292 0.051526 11.670 0.000
L1.y2 0.263057 0.060086 4.378 0.000
========================================================================
Correlation matrix of residuals
y1 y2
y1 1.000000 0.596470
y2 0.596470 1.000000
"""
# y2 から y1 へのグレンジャー因果性の確認
test_result = result.test_causality(caused="y1", causing="y2")
print(test_result.summary())
"""
Granger causality F-test. H_0: y2 does not Granger-cause y1. Conclusion: fail to reject H_0 at 5% significance level.
==============================================
Test statistic Critical value p-value df
----------------------------------------------
3.407 3.890 0.066 (1, 192)
----------------------------------------------
"""
# y1 から y2 へのグレンジャー因果性の確認
test_result = result.test_causality(caused="y2", causing="y1")
print(test_result.summary())
"""
Granger causality F-test. H_0: y1 does not Granger-cause y2. Conclusion: reject H_0 at 5% significance level.
==============================================
Test statistic Critical value p-value df
----------------------------------------------
136.2 3.890 0.000 (1, 192)
----------------------------------------------
"""
当然データによって結果が変わる(乱数で生成しているので、モデルが同じでも実行するたびに変わります。)のですが、
y1からy2へのグレンジャー因果性がないという過程は棄却されているので、y1はy2に影響しているようです。
y2からy1の向きは、そうはならなかったですね。
さて、上のコード例のように、 test_causality というメソッドに、
caused と causing の引数で変数名を渡してあげるだけで、テスト結果を返してくれました。
例では変数名を指定しましたが、 0 や 1など、何番目の変数なのかを指定してもOKです。
これが２変数ではなくもっと多くの変数のVARモデルで、多変数から多変数へのグレンジャー因果性を調べる場合は、
causedやcausingをリストで指定すれば、同じようにテストできます。
さて、 test_causality は CausalityTestResults というオブジェクトを返してくれますが、
上の例のように、 summary() というメソッドを持っているので、人間の目で見る場合は、これをprintするのが良いでしょう。
p値などのプログラムで使いたい場合は、それぞれ変数に持っているのでそちらを使えます。
# 帰無仮説
print(test_result.h0)
# H_0: y1 does not Granger-cause y2
# 統計量
print(test_result.test_statistic)
# 136.18235785962648
# p値
print(test_result.pvalue)
# 3.9579973656016236e-24
# 結論
print(test_result.conclusion_str)
# Conclusion: reject H_0