Pythonで時系列解析がしたい(多変量に対してVARを使う)

時系列間の関係を確認したいとき、どちらも定常であればVARモデルやVMAモデルが使えます。 特にVARは時系列間のlagを細かくこちらで仮定しきらなくて良いので便利です。

この書籍のRのVARをpythonで書き直してみます。

データの用意

書籍で用いるデータをRから取得します。
Rとのインターフェイスとして、pyperを使いました。 esu-ko.hatenablog.com

import pandas as pd

import pyper
r = pyper.R(use_pandas='True')
code="""
library(FinTS)
data(m.ibmspln)
"""
res = r.get("m.ibmspln")
df = pd.DataFrame(res,columns = ["ibm","sp"])

単位根過程かの検定を行う

定常の前提が欲しいので、adf検定をかけます。 esu-ko.hatenablog.com

from statsmodels.tsa.stattools import adfuller

adfuller(x=df.ibm,regression="c")
adfuller(x=df.sp,regression="c")
(-27.572016610758457,
 0.0,
 0,
 887,
 {'1%': -3.437743827988169,
  '5%': -2.864803950716061,
  '10%': -2.5685079558676054},
 5759.594725186838)

(-7.085102028822751,
 4.559991788038839e-10,
 20,
 867,
 {'1%': -3.437914898140353,
  '5%': -2.864879373440662,
  '10%': -2.5685481313814624},
 5448.609024607928)

若干、spで用いるlagが書籍と異なりましたが、単位根過程であるという帰無仮説を棄却できました。

varモデルをつくる

適切なパラメータを探します。

from statsmodels.tsa.vector_ar.var_model import VAR

maxlags = 5
model = VAR(df)
selected_orders = model.select_order(maxlags,trend='c').selected_orders
{'aic': 3, 'bic': 0, 'hqic': 0, 'fpe': 3}

bic,hqicの結果が書籍とは異なりますね。(書籍だとどちらも1)
また、maxlagsの値によってaicなども変わっていたので、この辺のいい感じの値を探す方法は宿題としておきます。

というわけでモデリングします。

res = model.fit()
res.summary()

とこの時に、select_orderと同じようにmaxlags,ic,trendなどを渡したり、渡さなかったりで適切なlagを探してくれます。 上の結果が書籍と違うこともあり、ひとまず、勝手にフィッティングしてもらいます。

結果、lag=1となり各パラメータ含めて書籍と同じ結果になりました。

グレンジャー因果/インパルス応答なんかをみる

varを使うメリットとして、上記のものなんかもみておきます。

#インパルス応答
ip = res.irf()
ip.plot()

#グレンジャー因果
#causingがRのcauseの意味なので注意
c1 = res.test_causality(caused='sp',causing='ibm')
c2 = res.test_causality(caused='ibm',causing='sp')

print(c1)
print(c2)

インパルス応答はデフォルトで全指標間のショックを可視化してくれます。

グレンジャー因果は書籍と同様ibm->spは棄却されず、sp->ibmは棄却され、グレンジャー因果性が存在しました。
なお、causedだけ渡すと、自身で自身の因果性を検証しに行くので注意が必要です。