1.やりたいこと

日経平均株価(NIKKEI225)はS&P500と相関があることはよく知られた話である。今回、多変量の時系列データ分析を扱う際に使用されるVAR（ベクトル自己回帰）モデルを使用してNIKKEI225とS&P500の関係性（因果探索）の分析を行なってみる。

2.VAR（ベクトル自己回帰）モデル

ベクトル自己回帰(VAR)モデルは、自己回帰(AR)モデルを多変量に拡張したものである。互いに相関する複数の時系列データを分析する際に使用する。
VARモデルは以下の式で示めされるように過去のp期前までの各系列データの値を使用して現時点の各系列データの値を予測しようとするモデルである。

$y_t=c+\Phi_1 y_{t-1}+ \dots +\Phi_1 y_{t-1}+\epsilon_t, \epsilon_t \sim W.N.(\Sigma)$

ここで、 $W.N.(\Sigma)$ はベクトルホワイトノイズ（撹乱項）であり、以下の式を満たす。

$E(\epsilon_t)=0$

$E(\epsilon_t \epsilon_{t-k}^t) = \left\{ \begin{array}{ll} \Sigma & (k = 0) \\ 0 & (k \neq 0) \end{array} \right.$

つまり、自己相関を持たず、異時点で系列データ間でも相関を持たない。ただ、同一時点では相関を許容する。( $\Sigma$ を分散分散行列とする。)
VARモデルの推定は、撹乱項を多変量正規分布に従うとするなどして、各系列データに関する方程式を個別にOLSすることによって推定していく。
VARモデル解析においては、グレンジャー因果性、インパルス応答解析、分散分解という３つのツールが重要な役割を果たす。ここではグレンジャー因果性とインパルス応答解析を扱う。グレンジャー因果性はある変数(群)が他の変数(群)の予測の向上に役立つかどうか判定するものであり、インパルス応答解析はある変数に対するショックがその変数やその他の変数の値に与える影響を分析するものである。グレンジャー因果性やインパルス応答解析を行うことで系列データ間の因果関係の探索、予測精度向上に重要なラグや変数の特定および変数の影響度を把握することが可能になります。例えば、目的変数yに異常が生じた場合に系列データ間の関連性(因果関係)を把握できていれば、原因の特定につながる。また、影響にラグがある場合にどの程度影響が残存してしまうのかは多変量の系列データを制御する場合に重要な情報となる。

3.グレンジャー因果性

前説で簡単に説明したが、ここではもう少しだけ詳しく説明する。
基本的な考え方は、時系列 $y_t$ の予測において、時系列 $x_t$ を用いると予測精度が向上した場合、 $x_t$ から $y_t$ への因果性があるという考え方である。ここで、２変量VAR(2)モデルについて考える。

$\left\{ \begin{array} ~y_{1t}=c_1+\phi_{11}^{(1)} y_{1,t-1}+\phi_{12}^{(1)} y_{2,t-1}+\phi_{11}^{(2)} y_{1,t-2}+\phi_{12}^{(2)} y_{2,t-2}+\epsilon_{1t} \\ y_{2t}=c_2+\phi_{21}^{(1)} y_{1,t-1}+\phi_{22}^{(1)} y_{2,t-1}+\phi_{21}^{(2)} y_{1,t-2}+\phi_{22}^{(2)} y_{2,t-2}+\epsilon_{2t} \end{array} \right.$

この時、 $y _ {2t}$ から $y _ {1t}$ へのグレンジャー因果性が存在している場合、 $\phi _ {12}^{(1)}=\phi _ {12}^{(2)}=0$ と同値になる。逆にグレンジャー因果性がない場合は、係数が全て０であると同値になる。そのため、グレンジャー因果性を検定するためには,F検定を用いて $\phi _ {12}^{(1)}=\phi _ {12}^{(2)}=0$ を検定し求めていく。

手順としては、はじめに、VARモデルにおける $y _ kt$ のモデルをOLSで推定し、その残差平方を $SSR _ 1$ とする。次に、VARモデルにおける $y _ kt$ のモデルに制約を貸したモデルをOLSで推定し、その残差平方を $SSR _ 0$ とする。 F統計量を次のように算出します。

$F=\frac{(SSR_0-SSR_1)/r}{SSR_1/(T-np-1)}$

ここで, $r$ はグレンジャー因果性検定に必要な制約数である。そして、 $rF$ を $\chi ^ 2(r)$ の95%点と比較し、 $rF$ の方が大きければ、ある変数群からの $y _ {kt}$ へのグレンジャー因果性は存在し、小さければグレンジャー因果性は存在しないと結論付ける。グレンジャー因果性はグレンジャー因果性は通常の因果性とは異なることに注意したい。グレンジャー因果性は通常の因果性が存在する必要条件ではあるが、十分条件ではない。また、因果の方向が実際と真逆になることもある。そのため因果の方向がはっきりしている場合をのぞいて、定義通りに予測に優位かどうかという観点で解釈するのが良いとされている。

4.インパルス応答解析

ある変数時系列の変化がその変数時系列やその他の変数時系列に与える影響を分析することができるのが、インパルス応答関数である。変化の識別の仕方によって、複数のインパルス応答関数が存在し、直行化インパルス応答関数と非直行化インパルス応答関数とがある。

非直行化インパルス応答解析
$y _ {jt}$ の撹乱項 $\epsilon _ {jt}$ だけに1単位（または1標準偏差）のショックを与えたときの $y _ {i,t+k}$ の変化は、 $y _ j$ のショックに対する $y _ i$ の $k$ 期後の非直交化インパルス応答と呼ばれる。また、それを $k$ の関数としてみたものは、 $y_j$ のショックに対する $y_i$ の非直交化インパルス応答関数といわれる。
直行化インパルス応答解析
撹乱項の分散共分散行列の三角分解を用いて、撹乱項を互いに相関する部分と無相関な部分に分解したとき、無相関な部分は直交化撹乱項といわれる。また、 $y_j$ の直交化撹乱項に1単位または1標準偏差のショックを与えたときの $y_i$ の直交化インパルス応答を時間の関数としてみたものは、 $y_j$ に対する $y_i$ の直交化インパルス応答関数といわれる。

5.NIKKEI225とS&P500の因果探索(Pythonでの実装)

それでは、Pythonのstatsmodelライブラリを使用してVARモデルの構築からグレンジャー因果性の統計的仮説検定、インパルス応答解析をおこなってみる。今回は冒頭に書いたようにNIKKEI225とS&P500の関係性を分析してみる。

5.1.データの用意

pandas_datareaderモジュールを使用してNIKKEI225、S&P500のデータを取得した。

import pandas_datareader.data as web
import datetime

#2013年10月10日を収集開始時点にする
start =datetime.date(2013,10,10)
#2023年9月30日を収集終了時点にする
end = datetime.date(2023,9, 30)

df_nikkei = web.DataReader('NIKKEI225', 'fred', start, end)
df_sp500 = web.DataReader('SP500', 'fred', start, end)

#データがない日付を特定
drop_ind_nikkei = df_nikkei[df_nikkei['NIKKEI225'].isnull() == True].index
drop_ind_sp500 = df_sp500[df_sp500['SP500'].isnull() == True].index
drop_ind = drop_ind_nikkei.append(drop_ind_sp500)
#データがない日付を削除
df_nikkei = df_nikkei.drop(drop_ind)
df_sp500 = df_sp500.drop(drop_ind)

取得したデータを描画してみる。

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(16, 8))

ax1 = fig.add_subplot(2, 1, 1)
ax1.plot(df_nikkei.index, df_nikkei['NIKKEI225'])
ax1.set_title('NIKKEI225')

ax2 = fig.add_subplot(2, 1, 2)
ax2.plot(df_sp500.index, df_sp500['SP500'])
ax2.set_title('S&P500')
plt.show()

相関がありそうなことが確認できる。

5.2.定常性の確認

上図から明らかなように両系列データともに定常性を満たしていない。一応ADF検定により定常性を確認してみる。

import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller

#ADF検定(単位根検定、元データ)
c_results = adfuller(df_nikkei['NIKKEI225'], regression='c')
print('NIKKEI225:', c_results[1])
c_results = adfuller(df_sp500['SP500'], regression='c')
print('SP500:', c_results[1])

仮に有意水準5%と仮定した場合、両系列データの検定値は0.05より大きく定常性が認められない。

NIKKEI225: 0.7150331280648695
SP500: 0.873103776458502

上記結果から定常性が認められなかったので、そのままのデータでは定常性を前提とするVARモデルによるモデル化ができない。
この問題に対処するため、対数差分処理を行う。一般に下記に示すような対数差分処理や差分処理を行うことで定常になるように変換する。

$log(y_t) - log(y_{t-1}) = log( \dfrac{y_{t}}{y_{t-1}}) = log(1 + \dfrac{y_{t} - y_{t-1}}{y_{t-1}}) = \dfrac{y_{t} - y_{t-1}}{y_{t-1}}$

対数差分処理を行うことで前の時刻からの変化率に変換できる。変化分が小さいと仮定すると、一次のテイラー展開近似により変化率に変換ができる。今回はこれに100をかけることで割合にしておく。

#対数差分をとって変化率に変換する処理
import numpy as np

#時系列データにおける対数差分は、前期からの比率を表すので、このあとの結果がわかりやすくなるように100をかけています。
df_nikkei_log = np.log(df_nikkei['NIKKEI225']).diff(1)*100
df_nikkei_log = df_nikkei_log[1:]
df_sp500_log = np.log(df_sp500['SP500']).diff(1)*100
df_sp500_log = df_sp500_log[1:]

対数差分処理後のデータを描画してみる。

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(16, 8))

ax1 = fig.add_subplot(2, 1, 1)
ax1.plot(df_nikkei.index[1:], df_nikkei_log)
ax1.set_title('NIKKEI225')
ax1.set_ylim(-10, 10)

ax2 = fig.add_subplot(2, 1, 2)
ax2.plot(df_sp500.index[1:], df_sp500_log)
ax2.set_title('S&P500')
ax2.set_ylim(-10, 10)
plt.show()

定常データになっていそうなことが確認できる。

再度ADF検定を行い、定常性の確認を行う。

#ADF検定(単位根検定、対数差分)
c_results = adfuller(df_nikkei_log, regression='c')
print('NIKKEI225:', c_results[1])
c_results = adfuller(df_sp500_log, regression='c')
print('SP500:', c_results[1])

NIKKEI225: 0.0
SP500: 6.1407766905440326e-21

定常な系列データに変換することができた。一度データフレーム形式で両系列データをまとめておく。

df_nikkei_log = pd.DataFrame(df_nikkei_log, columns=['NIKKEI225'])
df_sp500_log = pd.DataFrame(df_sp500_log, columns=['SP500'])
df = pd.concat([df_nikkei_log, df_sp500_log], axis=1)

いよいよVAR(ベクトル自己回帰モデル)を構築し、グレンジャー因果性に関する統計的仮説検定、インパルス応答解析を行なっていく。その前に相互相関を確認しておく。

#相互相関コレログラム
fig = plt.figure(figsize=(16, 2))

plt.xcorr(df['NIKKEI225'], df['SP500'])
plt.title('NIKKEI225*S&P500')
plt.show()

相互相関コレログラムは以下のようになり、0期(同一時点)と1期ずれた点で相関が高いことが確認できる。

5.3.VARモデルの構築

statsmodelsライブラリのtsa.vector_ar.var_model を使用することでVARモデルの推定が行える。今回は最大となるラグを指定し、ベイズ情報量(BIC)を基準にラグを最適化し、最適化したラグでVARモデルの推定を行った。

# VARモデルの構築
from statsmodels.tsa.vector_ar.var_model import VAR

# 最大のラグ数
maxlags = 10

# モデルのインスタンス生成
var_model = VAR(df)

# 最適なラグの探索
# ベイズ情報量基準で最適なラグを決定する
lag = var_model.select_order(maxlags).selected_orders
print('最適なラグ：',lag['bic'],'\n')

# モデルの学習
# 最適なラグでVARモデルを学習させる
results = var_model.fit(lag['bic'])

# 結果出力
print(results.summary())

ベイズ情報量基準(BIC)による最適なラグ探索結果と最適なラグで学習させた時のVARモデルのサマリーが出力される。また最後に量系列データの相関行列が表示される。

最適なラグ： 2 

  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Mon, 16, Oct, 2023
Time:                     17:44:44
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                   0.427853
Nobs:                     2350.00    HQIC:                  0.412263
Log likelihood:          -7132.93    FPE:                    1.49680
AIC:                     0.403333    Det(Omega_mle):         1.49046
--------------------------------------------------------------------
Results for equation NIKKEI225
===============================================================================
                  coefficient       std. error           t-stat            prob
-------------------------------------------------------------------------------
const                0.011922         0.023616            0.505           0.614
L1.NIKKEI225        -0.203864         0.021595           -9.440           0.000
L1.SP500             0.543548         0.022090           24.606           0.000
L2.NIKKEI225        -0.002393         0.019161           -0.125           0.901
L2.SP500             0.180303         0.024491            7.362           0.000
===============================================================================

Results for equation SP500
===============================================================================
                  coefficient       std. error           t-stat            prob
-------------------------------------------------------------------------------
const                0.040638         0.023308            1.744           0.081
L1.NIKKEI225         0.027808         0.021313            1.305           0.192
L1.SP500            -0.144626         0.021802           -6.634           0.000
L2.NIKKEI225        -0.019569         0.018911           -1.035           0.301
L2.SP500             0.083522         0.024172            3.455           0.001
===============================================================================

Correlation matrix of residuals
             NIKKEI225     SP500
NIKKEI225     1.000000  0.322193
SP500         0.322193  1.000000

モデルサマリーをみると、NIKKEI225の予測に関しては、１日前のNIKKEI225と１日前、２日前のS&P500の値を使用することは優位であることがわかる。（有意水準5%の場合）

Results for equation NIKKEI225
===============================================================================
                  coefficient       std. error           t-stat            prob
-------------------------------------------------------------------------------
const                0.011922         0.023616            0.505           0.614
L1.NIKKEI225        -0.203864         0.021595           -9.440           0.000
L1.SP500             0.543548         0.022090           24.606           0.000
L2.NIKKEI225        -0.002393         0.019161           -0.125           0.901
L2.SP500             0.180303         0.024491            7.362           0.000
===============================================================================

一方、S&P500の予測に関しては、1日前と２日前のS&P500の値を使用することは統計的に優位であるが、NIKKEI225の値は優位でないと分かる。つまり、NIKKEI225はS&P500の影響を受けるが、S&P500はNIKKEI225の影響をほぼ受けないことが分かる。

Results for equation SP500
===============================================================================
                  coefficient       std. error           t-stat            prob
-------------------------------------------------------------------------------
const                0.040638         0.023308            1.744           0.081
L1.NIKKEI225         0.027808         0.021313            1.305           0.192
L1.SP500            -0.144626         0.021802           -6.634           0.000
L2.NIKKEI225        -0.019569         0.018911           -1.035           0.301
L2.SP500             0.083522         0.024172            3.455           0.001
===============================================================================

そのため、グレンジャー因果性的には、S%P500→NIKKEI225の向きに因果があると想定されます。

5.4.グレンジャー因果性の統計的仮説検定

続いて推定したVARモデルでグレンジャー因果性の統計的仮説検定を行う。以下のコードで実行できる。

# グレンジャー因果性の検定
# 帰無仮説：グレンジャー因果なし
# 対立仮説：グレンジャー因果あり

for i in range(len(df.columns)):
    for j in range(len(df.columns)):
        if i != j :
            #因果の検定
            test_results = results.test_causality(causing=i, caused=j) 
            #p値
            test_results.pvalue
            #検定結果の出力
            print(df.columns[i],'->',df.columns[j],'\t',test_results.pvalue)

検定結果は以下の通りである。SP500 -> NIKKEI225の方向の因果は統計的に優位であることが確認でき、NIKKEI250を予測するのにS&P500は有効であることが確認できる。

NIKKEI225 -> SP500     0.19007514052849195
SP500 -> NIKKEI225    2.227002166652272e-124

5.5.インパルス応答解析

インパルス応答関数は、statsmodelsモジュールを使えば簡単に実行できる。インパルス応答関数は、非直交化と直交化の２種類あり、plotのorth引数で切り替えることができる。orthをTrueにすると直交化された結果が表示される。まず非直行化インパルス応答解析を行なってみる。

5.5.1 非直行化インパルス応答解析

５期先の影響までを描画する。

# インパルス応答関数(直交化)

irf = results.irf(5)
irf.plot(orth=False)
plt.show()

上記のインパルス応答解析結果から、S&P500に１単位のショックを与えた1日後に大きくNIKKEI225の値が増えていることが分かる。一方でNIKKEI225に１単位分のショックを与えた場合に正の方向への変化は見られない。インパルス応答解析結果は、VARモデルサマリーやグレンジャー因果性の検定結果と同等の分析結果を示唆した。

非直行化インパルス応答解析を行う場合、1単位分(今回は1%変化分)のショックを撹乱項に与えた場合の結果になる。外乱項の１偏差分のショックを与えた場合の影響を定量的に把握したい場合は、以下で求められる外乱項の偏差を応答値に書けることで確認できる。

print("ε_NIKKEI225:",np.sqrt(results.sigma_u['NIKKEI225']['NIKKEI225']))
print("ε_S&P500:",np.sqrt(results.sigma_u['SP500']['SP500']))

それぞれの撹乱項の偏差は以下である。

ε_NIKKEI225 1.1430889085971325
ε_S&P500 1.1281818499705119

ちなみにインパルス応答解析の生データはirf.irfsで取得できるので、１偏差のショックを与えた時のインパルス応答解析結果のプロットも可能である。

fig = plt.figure(figsize=(12, 12))

#NIKKEI225→NIKKEI225
ax1 = fig.add_subplot(2,2,1)
ax1.plot(irf.irfs[:, 0, 0]*np.sqrt(results.sigma_u['NIKKEI225']['NIKKEI225']))
ax1.set_title('NIKKEI225→NIKKEI225')

#SP500→NIKKEI225
ax2 = fig.add_subplot(2,2,2)
ax2.plot(irf.irfs[:, 0, 1]*np.sqrt(results.sigma_u['SP500']['SP500']))
ax2.set_title('SP500→NIKKEI225')

#NIKKEI225→SP500
ax3 = fig.add_subplot(2,2,3)
ax3.plot(irf.irfs[:, 1, 0]*np.sqrt(results.sigma_u['NIKKEI225']['NIKKEI225']))
ax3.set_title('NIKKEI225→SP500')

#SP500→SP500
ax4 = fig.add_subplot(2,2,4)
ax4.plot(irf.irfs[:, 1, 1]*np.sqrt(results.sigma_u['SP500']['SP500']))
ax4.set_title('SP500→SP500')

plt.show()

5.5.2 直行化インパルス応答解析

続いて直行化インパルス応答解析を行う。

# インパルス応答関数(直交化)

irf = results.irf(5)
irf.plot(orth=True)
plt.show()

非直行化インパルス応答解析では現れていなかったNIKKEI225→S&P500の0期での影響が大きく出ていることが確認できる。直行化することでNIKKEI225にショックを与えたことでそれと相関する撹乱項が同じ時点のS&P500にも与えられるために0期の値が大きくなったと考えられる。また、直行化インパルス応答解析では、１偏差分のショックが与えられた場合の結果がでてくる。この際の標準偏差は下記で導出できる。

sigma = results.sigma_u

L = np.linalg.cholesky(sigma)
L_inv = np.linalg.inv(L)
L_t_inv = np.linalg.inv(L.T)
D = L_inv * sigma * L_t_inv

new_sigma_sqrt = np.sqrt(np.diag(D))
new_sigma_sqrt

今回のケースでは、下記のようになり、ショックとしてNIKKEI225に1%増加を与えた場合、S&P500に1.05%の増加を与えた場合の結果がインパルス応答解析結果として描画されていることになる。

array([1.        , 1.05632992])

ここでもう少しだけインパルス応答解析結果に関して、考察しておく。S&P500→NIKKEI225の0期の影響が０になっているが、これは日本と米国の株式市場の開いてる時間に起因する。米国市場は日本市場が閉まった後に開くため、同一期(同じ日)にS&P500→NIKKEI225に影響（インパルス）があるはずがない。これが反映されているために0になっていると考えられる。
また、インパルス応答解析では系列データの並びが重要である。今回は市場が開く順番(日本→米国)でデータを並べたが、並びが変わるとインパルス応答解析結果も変わってしまう。仮にS&P500→NIKKEI225に並びを変えてインパルス応答解析を行なってみると、以下のようになる。

#VARモデル作成時にNIKKEI225とS&P500の並びを変えるとインパルス応答解析の結果が変わる。
df = pd.concat([df_sp500_log, df_nikkei_log], axis=1)

# VARモデルの構築
from statsmodels.tsa.vector_ar.var_model import VAR

# 最大のラグ数
maxlags = 10

# モデルのインスタンス生成
var_model = VAR(df)

# 最適なラグの探索
# ベイズ情報量基準で最適なラグを決定する
lag = var_model.select_order(maxlags).selected_orders
print('最適なラグ：',lag['bic'],'\n')

# モデルの学習
# 最適なラグでVARモデルを学習させる
results = var_model.fit(lag['bic'])

# インパルス応答関数(直交化)

irf = results.irf(5)
irf.plot(orth=True)
plt.show()