Pythonのstatsmodels==0.10.0を用いて、主に因果推論のやり方をメモ程度に書きました。
それぞれの説明は全くないですがご容赦下さい。
単位根検定
Augmented Dickey-Fuller 単位根検定
from statsmodels.tsa.stattools import adfuller
results = adfuller(x) # x: array_like, 1d
stat = results[0] # 統計量
pvalue = results[1] # p値
pvalueが基準(0.05など)より小さい場合、単位根過程であるという帰無仮説が棄却
参考:
https://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.adfuller.html
共和分検定
Augmented Engle-Granger two-step 共和分検定
from statsmodels.tsa.stattools import coint
results = coint(x, y) # # x, y: array_like, 1d
stat = results[0] # 統計量
pvalue = results[1] # p値
pvalueが基準(0.05など)より小さい場合、共和分の関係にないという帰無仮説が棄却
参考:
https://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.coint.html
グレンジャー因果性検定
前提
下のブログにあるように、以下のようなステップを踏んでから因果性検定を行います
- VARモデルを推定する前に、まず個々の原系列に対して単位根検定を実施する
- そもそもただの回帰関係が見たいだけなら、気にせずVARモデルを計算して良い
- 単位根過程が混ざっていなければ、そのままVARモデルを計算して因果性分析まで行って良い
- 単位根過程が混ざっていたら、まず共和分関係（ランク）を推定する
- 共和分関係がなければ、差分系列に対してVARモデルを推定して因果性分析まで行って良い
- 共和分関係があったら、VECMを推定してからVARモデルに変換して初めて因果性分析を行える
単位根検定、共和分検定は上で説明したように行ってください。
VECMモデルに関してはたぶんいつか追記します。
以下では、上のステップを踏んだ後のVARモデルを推定して因果性分析をするところについて書きます。
その１
####VARモデルで近似して、icがもっとも低くなる次数を見つけてから因果性検定
from statsmodels.tsa.vector_ar.var_model import VAR
maxlags = 100
model = VAR(x) # x: ndarray (data数, column数)
selected_orders = model.select_order(maxlags).selected_orders
# e.g. selected_orders = {'aic': 91, 'bic': 7, 'fpe': 91, 'hqic': 10}
results = model.fit(selected_orders['aic'])
test_results = results.test_causality(causing=0, caused=1) # x0=>x1への因果を検定
pvalue = test_results.pvalue
pvalueが基準(0.05など)より小さい場合、グレンジャー因果がないという帰無仮説が棄却
その２
####VARモデルの次数を決めうちで因果性検定
次数が決められるときは、もっと簡単にできるっぽいので一応書いておきますが、その１でやればいいと思います。
from statsmodels.tsa.stattools import grangercausalitytests
maxlag = 30
results = grangercausalitytests(x, maxlag) # x: 2darray (data数, 2), x[:,0] => x[:,1]への因果を検定
stat = results[0] # 統計量
p_value = results[1] # p値
p_valueが基準(0.05など)より小さい場合、x1=>x2のグレンジャー因果がないという帰無仮説が棄却
参考:
https://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.grangercausalitytests.html
インパルス応答
VARモデルを推定するところまでは、グレンジャー因果性のその１と同じです。
maxlags = 100
model = VAR(x) # x: ndarray (data数, column数)
selected_orders = model.select_order(maxlags).selected_orders
# e.g. selected_orders = {'aic': 91, 'bic': 7, 'fpe': 91, 'hqic': 10}
results = model.fit(selected_orders['aic'])
period = 25
irf = results.irf(period)
irf.plot()
こんな感じの図が出るはず