グレンジャー因果検定について - 機械学習 Memo φ(・ω・ )

概要

仕事でグレンジャー因果検定を使う機会があったので、グレンジャー因果検定について勉強したことを記載します。また、statsmodelsのAPIを使い、株式データを使って簡単な検定をしてみました。

概要
グレンジャー因果検定
株式データを使って検定
- statsmodelsのグレンジャー因果検定
おわり
- 参考

間違いがありましたら、コメントいただけたら嬉しいです。

グレンジャー因果検定

考え方

時系列データ $x$ , $y$ において、 $x$ が増減すると $y$ も同じように増減するという関係なのかを検証したい。

未来の $y$ の値の予測に、現在と過去の $y$ の値を使って予測した時より、 $x$ の値も加えて予測したほうが精度が改善される時、 $x$ から $y$ にGrengerの意味で因果があるという。

モデルにVARを利用し、F検定で統計的有意性の検定します。

f:id:YukoIshizaki:20200802155701p:plain:w400
Granger causality - Wikipedia

確かに $x$ が増減したら $y$ も同じように増減する関係なら、 $y$ の予測に $x$ の値も知っているモデルの方が精度が高く、その逆も言えると考えるのは感覚的も理解しやすいです。しかし、一般的に変数が多ければモデルの精度は高くなりやすいので、それが統計的有意性があるかを調べる必要があり、それを２つのモデルの誤差の比が十分に大きいかF検定で調べます。

ARモデル(自己回帰モデル)

グレンジャー因果検定で使うVARモデル(ベクトル自己回帰モデル)を理解するために、ARモデル(自己回帰モデル)を理解する。

AR(p)モデルは、 $t$ 時点での $y$ の値を $p$ 期内の $y$ の過去値の加重和で表現したモデル。 $c$ は定数項、 $\phi$ は自己回帰係数、 $\varepsilon_t$ はホワイトノイズ(期待値0, 分散 $\sigma^2$ を持ち、自己共分散0の誤差項)。

$y_t = c + \phi_1 y_{t-1} + \cdots +\phi_p y_{t-p} + \varepsilon_t$

( $y_t \in \mathbb{R}, \phi_i \in \mathbb{R}$ )
(正確には $y$ は実数に値をもつ確率変数)

VARモデル(ベクトル自己回帰モデル)

VAR(p)モデルは、AR(p)モデルをベクトルに拡張したもの。 $\textbf{c}$ と $\varepsilon_t$ は $n \times 1$ の列ベクトル、 $\Phi_i$ は $n \times n$ の行列。

$\textbf{y}_t = \textbf{c} + \Phi_1 \textbf{y}_{t-1} + \cdots +\Phi_p \textbf{y}_{t-p} + \varepsilon_t$

( $\textbf{y}_t \in \mathbb{R}^n, \Phi_i \in M(n, \mathbb{R})$ )
(正確には $\textbf{y}$ は実n次元数ベクトルに値をもつ確率変数ベクトル)

例えば、下記のようなVAR(2)モデルは

$\left( \begin{array}{c} x_t \\ y_t \\ \end{array} \right) = \left( \begin{array}{c} c_1 \\ c_2 \\ \end{array} \right) + \left( \begin{array}{cc} \phi^{(1)}_{11} & \phi^{(1)}_{12} \\ \phi^{(1)}_{21} & \phi^{(1)}_{22} \\ \end{array} \right) \left( \begin{array}{c} x_{t-1} \\ y_{t-1} \\ \end{array} \right)+ \left( \begin{array}{cc} \phi^{(2)}_{11} & \phi^{(2)}_{12} \\ \phi^{(2)}_{21} & \phi^{(2)}_{22} \\ \end{array} \right) \left( \begin{array}{c} x_{t-2} \\ y_{t-2} \\ \end{array} \right) +\left( \begin{array}{c} \varepsilon_{xt} \\ \varepsilon_{yt} \\ \end{array} \right)$

下記のように書き換えられる

$\begin{cases} x_t = c_1 + \phi^{(1)}_{11} x_{t-1} + \phi^{(1)}_{12} y_{t-1} + \phi^{(2)}_{11} x_{t-2} + \phi^{(2)}_{12} y_{t-2} + \varepsilon_{xt} \\ y_t = c_2 + \phi^{(1)}_{21} x_{t-1} + \phi^{(1)}_{22} y_{t-1} + \phi^{(2)}_{21} x_{t-2} + \phi^{(2)}_{22} y_{t-2} + \varepsilon_{yt} \end{cases}$

検定方法

考え方で記載したとおり、 $y_t$ を予測するとき、 $x$ の値を含めた方が精度が高い、つまり

$y_t = c_2 + \phi^{(1)}_{21} x_{t-1} + \phi^{(1)}_{22} y_{t-1} + \phi^{(2)}_{21} x_{t-2} + \phi^{(2)}_{22} y_{t-2} + \varepsilon_{yt} \:\:\:\:\:\:\:\:\:\: \cdots(a)$

のほうが

$y_t = c_2 + \phi^{(1)}_{22} y_{t-1} + \phi^{(2)}_{22} y_{t-2} + \varepsilon_y \:\:\:\:\:\:\:\:\:\:\:\:\:\:\: \cdots(b)$

より誤差が小さいことを示したい。
$x$ から $y$ にGrengerの意味で因果がないことを帰無仮説とし、２つのモデルの誤差を使いF検定で帰無仮説を棄却する流れになる。

上記の式だと、 $H_0: \phi^{(1)}_{21} = \phi^{(2)}_{21} =0$ ( $x$ の係数が $0$ )と同値になる。

F検定統計量は、式 (a) をOLSで推定した時の誤差(残差平方和) $SSR_1$ と、式 (b) をOLSで推定した時の誤差(残差平方和) $SSR_0$ を使って以下のように表され、 $2F$ は漸近的にカイ2乗分布 ( $\chi^2(r)$ ) に従うので、 $2F$ の値が $\chi^2(r)$ の95%点より大きければ帰無仮説を棄却できる。

$F = \displaystyle \frac{(SSR_0 - SSR_1 )\: \: / \: 2}{SSR_1 \: / \: (T-5)}$

ここで $r$ は制約(xの係数を0とする制約)の数、 $T$ はサンプル数。

F検定

F検定では２つの分散の比がF分布に従うかを調べることで、帰無仮説を母分散が等しい(２つの分散に差がない、今回の場合は２つのモデルの誤差に差がない)、対立仮説を母分散が等しくない(２つの分散に差がある、今回の場合は２つのモデルの誤差に差がある)として、検定を行います。

上記の説明ではF検定統計量を使用していますが、以下の例に記載するstatsmodelsの例ではp値を使用し、p値が0.05未満であれば有意水準5%で有意であり帰無仮説は棄却するとして、検定を行います。

注意点

サンプル数が多い時

単純にサンプル数が多ければ多いほどp値は小さくなりやすいので、サンプル数が多い時はp値に惑わされないよう注意する必要があるようです。

交絡因子について

交絡因子(間接的に影響する変数)が存在していても「Grengerの意味で因果がある」となる可能性があります。因果関係を調べるには、十分ではないことに注意する必要があるようです。

株式データを使って検定

statsmodelsのグレンジャー因果検定

statsmodelsに、grangercausalitytestsが用意されているのでそれを使い、データはkaggleのdatasetのS&P500のデータを使います。

statsmodels ドキュメント: statsmodels.tsa.stattools.grangercausalitytests — statsmodels
Kaggle Dataset : S&P 500 stock data | Kaggle

株に関しては完全に素人なのですが、単純に前日の市場が開いてる間に値上がりすれば(前日終値 - 前日始値 > 0)買いな株と判断され、前場の寄り付き前に値上がりして、当日の始値は前日の終値より高いじゃないか(当日始値 - 前日終値 > 0)、逆に前日に値が下がったら、寄り付き前に値下がるのじゃないかと考え検証してみます。

値上がり金額を $x$ とし、始値と前日終値との差を $y$ とした時、 $x$ から $y$ へGrengerの意味で因果があるかを調べます。

grangercausalitytestsの第一引数にはarray_likeなデータを指定しますが、

The data for test whether the time series in the second column Granger causes the time series in the first column.

と書かれている通り、１番目のカラムには結果側のデータ $y$ 、２番目のカラムには起因側のデータ $x$ となっている配列を指定します。今回は前日の値上がり金額に対して検証したいのでmaxlagを1にします。

import pandas as pd 
from statsmodels.tsa.stattools import grangercausalitytests
import matplotlib.pyplot as plt
pd.options.plotting.backend = "plotly"

PATH = '/kaggle/input/sandp500/individual_stocks_5yr/individual_stocks_5yr/'

# Exxon Mobil Corporation のデータ
XOM = pd.read_csv(f'{PATH}/XOM_data.csv', index_col='date', parse_dates=True)
# 期間を限定
XOM = XOM[(XOM.index >= pd.to_datetime('2017-07-01')) 
        & (XOM.index <= pd.to_datetime('2017-12-31'))]
# 終値と始値の差
XOM['price_increase'] = XOM['close'] - XOM['open']
# 始値と前日終値との差
XOM['before_opening'] = XOM['open'] - XOM['close'].shift()

# データ確認
XOM.head()
XOM[['price_increase', 'before_opening']].plot()

# グレンジャー因果検定の結果取得
x = XOM[['before_opening', 'price_increase']][1:].values
results = grangercausalitytests(x=x, maxlag=[1],  verbose=False)
print(results)

f:id:YukoIshizaki:20200803002323p:plain:w450
f:id:YukoIshizaki:20200803002419p:plain:w700

>> {1: ({'ssr_ftest': (2.137422914716022, 0.14633552782619325, 121.0, 1),
   'ssr_chi2test': (2.190416871279229, 0.1388717255785229, 1),
   'lrtest': (2.171295251305196, 0.1406077610113857, 1),
   'params_ftest': (2.137422914716034, 0.14633552782619272, 121.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x7f86b1e70c10>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x7f86b1e87d90>,
   array([[0., 1., 0.]])])}

４つのタイプで検定された結果が返ってきます。lagのkeyに対して、各検定の結果が格納されています。検定名がkeyで、valueは１番目の要素に検定統計量、２番目の要素にp値が入ってます。params_ftestとssr_ftestはF分布に従うとした場合の検定、ssr_chi2testとlrtestはカイ2乗分布に従うとした場合の検定です。

全ての検定でp値が0.05を大きく上回り帰無仮説が棄却されませんでした...やはりそんな簡単な話じゃなさそうですね😇

例えば、前日に値上げされたら、当日の始値と前日の始値の差はプラスになりそうです。前日終値が前日始値よりすでに高いので。値上げ金額を $x$ とし、始値と前日始値との差を $y$ とした時であれば、 $x$ から $y$ へGrengerの意味で因果がある、となりそうなので、検定してみます。

# 値上がり金額
XOM['price_increase'] = XOM['close'] - XOM['open']
# 始値と前日始値の差
XOM['diff_open'] = XOM['open'] - XOM['open'].shift()

XOM[['price_increase', 'diff_open']].plot()

x = XOM[['diff_open', 'price_increase']][1:].values
results = grangercausalitytests(x=x, maxlag=[1],  verbose=False)
results

f:id:YukoIshizaki:20200803005611p:plain:w700

>> {1: ({'ssr_ftest': (168.13795446547945, 1.220451046664481e-24, 121.0, 1),
   'ssr_chi2test': (172.30666408032607, 2.3194627990098294e-39, 1),
   'lrtest': (108.01805908422944, 2.663402872664279e-25, 1),
   'params_ftest': (168.13795446547954, 1.2204510466644514e-24, 121.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x7f86b1e87a90>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x7f86b1f10e90>,
   array([[0., 1., 0.]])])}

今回は４つの検定でp値が0.05を大きく下回ったため帰無仮説を棄却し、Grengerの意味で因果があると言えそうです。