やりたいこと
VAR過程はAR過程を多変数に拡張したものです。変数が増えることによってモデルの説明力が向上することが期待されるほか、変数間の動学的関係の分析手法を提供します。「グレンジャー因果性」、「インパルス応答関数」、「分散分解」です。
- グレンジャー因果性→ある変数がほかの変数の将来の予測に役に立つか
- インパルス応答関数→ある変数に与えたショックがその変数やほかの変数にどの程度影響を与えるか
- 分散分解→各変数の不確実性に対して各変数がどれだけ寄与しているか
今回の目的はグレンジャー因果性検定の実装です。但しプログラムはVARから実装したのでコードが少し長くなっています。VARの詳しい実装は私のGitHubにおいてあるのでそちらを見てください。卒論の付録部分がそれにあたります。なお、本記事を書くにあたってプログラムを少し書き換えました。それに伴いGitHubと異なる部分があると思いますがご了承ください。
グレンジャー因果性とは
まずは定義します。(1)
定義4.2(一般的なグレンジャー因果性) $x_{t}$ と $y_{t}$ をベクトル過程とする. また, $x$ と $y$ の現在と過去の値を含む, 時点 $t$ において利用可能な情報の集合を $\Omega_{t}$ とし, $\Omega_{t}$ から現在と過去の $y$ を取り除いたものを $\tilde{\Omega_{t}}$ とする. このとき, $\tilde{\Omega_{t}}$ に基づいた将来の $x$ の予測と, $\Omega_{t}$ に基づいた将来の $x$ の予測を比較して,後者のMSEのほうが小さくなる場合, $y_{t}$ から $x_{t}$ へのグレンジャー因果性が存在するといわれる. ここで, MSEの大小は行列の意味での大小であることに注意されたい.
具体例でみるとわかりやすいので、2変数VAR(1)の場合を考えてみましょう。2変数VAR(1)は次のようになります。
y_{1t} = c_{1} + \phi_{11} y_{1, t-1} + \phi_{21} y_{2, t-1} + noise_{1} \tag{1} \\
y_{2t} = c_{2} + \phi_{12} y_{1, t-1} + \phi_{22} y_{2, t-1} + noise_{2}
ここで、$y_{2}$ から $y_{1}$ のグレンジャー因果性が存在しないというのはモデルでいうと次のようになります。
\begin{align}
y_{1t} &= c_{1} + \phi_{11} y_{1, t-1} + noise_{1} \tag{2} \\
y_{2t} &= c_{2} + \phi_{12} y_{1, t-1} + \phi_{22} y_{2, t-1} + noise_{2}
\end{align}
つまり、$y_{1t}$ の予測に関して$y_{2, t-1}$がなんら役に立たないことを意味します。ただし逆は言えないので注意しましょう。
よくありがちなのが間違いとして、通常の因果関係に関しての間違いです。グレンジャー因果性は通常の因果関係があることの必要条件ではあるが、十分条件ではありません。
加えて、グレンジャー因果性の方向と通常の因果関係の方向が必ず一致するとも限らないです。$y$ から $x$ にグレンジャー因果性があるとしても、因果関係が $x$ から $y$ の方向にあることもあります。
また、グレンジャー因果性があるからと言って相関関係があると言えませんし、その逆も言えません。
グレンジャー因果性は、ある変数がほかの変数の予測精度を向上させることができるのかという点について議論していることに注意しましょう。
単位根VAR過程では通常のF検定統計量を用いたグレンジャー因果性検定を行ってはいけない
まず単位根VAR過程ってなんじゃいって話ですが、簡単に言うと過程の一つが単位根であるようなVAR過程のことです。
もう少し詳しくAR過程を例に考えます。AR過程は特性方程式の解のすべてに対して絶対値を取った値がが1より大きくなければならない、という定常条件があります。この方程式の解の一つが1になるとき単位根過程と呼びます。これを多変数に拡張したものが単位根VAR過程になります。もう少し詳しく見たいよって方は過去記事を見てください。
グレンジャー因果性検定ですが通常はF分布を用いることで検定を行うことができます（正確にF分布に従うわけではなく漸近的にのみ正当化されることに注意しましょう）。ここでモデルが単位根過程になると推定量の漸近分布に対して大数の法則や中心極限定理が適用できないので、用いる漸近分布が変わります。詳しくは沖本本の5章単位根過程をご覧ください。以上のことから単位根VAR過程にグレンジャー因果性検定を適用させることができません。
じゃあ、単位根になったらどうするのって話になるのですが、私としては、データを加工してモデルを作り直すというのが一番簡単だと思います。基本的に定常の系列データに対してモデルが定常条件を満たすはず、、、経験的に。ほかにも検定統計量を自分で設定してP値を計算してってすれば検定はできるんじゃないかと思っていますが、そこらへんは勉強不足でよくわからないのでここら辺でこの話を終わりたいと思います。
繰り返しになりますが単位根VAR過程ではグレンジャー因果性検定を単純に適用させることができないということです。
グレンジャー因果性検定の実装
具体的な実装は計量時系列分析4.3のP81～P82にのっとって作成しています。走らせるメソッドはgrangerだけです。先にこのメソッドの説明を軽くしたいと思います。全体のソースコードは最後にまとめて載せておきます。
まずは帰無仮説と対立仮説に対応するモデルを生成します。式(2)が帰無仮説、式(1)が対立仮説になります。次にfit!メソッドでモデルにデータを学習させます。余談ですがJuliaでは引数の値を変更するような破壊的メソッドには名前の最後に"!"を付けることになっています。次に検定統計量を計算するためにそれぞれの残差平方和を計算します。最終的な検定の判定はcrit_val < r*Fで行っています。rは制約で、モデルの係数を0にした数です。
引数は下にまとめました。
- data：時系列データ
- caused：上の例でいうと $y_{1}$
- causing：上の例でいうと $y_{2}$
- lags：モデルにどれくらい過去の自分のデータを含めるかを決める
function granger( data, caused, causing ; lags::Int=1 )
N, D = size(data)
# 帰無仮説のモデル--制約あり
# 2変数の場合、制約を課すとARモデルを推定することと同じになる
if D == 2
h_0 = ar( lags, rand( lags + 1 ), randn(1) )
else
h_0 = var( lags, rand( ( D - 1 )*lags + 1, D - 1 ), rand( D - 1, D - 1 ) )
end
# 対立仮設のモデル--制約なし
h_1 = var( lags, rand( D*lags + 1, D ), rand( D, D ) )
# 影響を与えると考えられる変数の列を削除する
data_pan = data[:, setdiff(begin:end, causing)]
fit!( h_0, data_pan )
fit!( h_1, data )
# F検定統計量
ssr0 = ssr( h_0, data_pan, caused )
ssr1 = ssr( h_1, data, caused )
# 制約の数
r = lags
tmp_up = ( ssr0 - ssr1 ) / r
tmp_down = ssr1 / ( N - D*lags - 1 )
F = tmp_up / tmp_down
# カイ二条分布
x² = Chisq( r )
# 95%点
crit_val = quantile( x², 0.95 )
println( "critical value：", crit_val )
println( "statistic value：", F )
print( "conclusion: " )
if crit_val < r*F
println( "exists granger causality" )
else
println( "not exists granger causality" )
end
end
ソースコード
using Statistics
using Random
using LinearAlgebra
using Distributions
mutable struct var
p # ラグ
params # パラメータ
Σ # 共分散行列
end
mutable struct ar
p
params
σ # 分散
end
# パラメータの学習
function fit!( m::var, data )
N, D = size(data)
# 計画行列
global Φ = zeros( N-m.p, D*m.p+1 )
# 定数項部分を1にする
Φ[:,1] .+= 1
# 計画行列を作る。詳しくは付録を参照
for i in 0 : m.p-1
Φ[ :, i*D+2 : i*D+D+1 ] += data[ m.p-i : end-1-i, : ]
end
Y = data[ m.p+1 : end, : ]
m.params[:,:] = inv(Φ'Φ)Φ'Y
tmp = ( Y'Y - Y'Φ*m.params - m.params'Φ'Y + m.params'Φ'Φ*m.params ) ./ ( length(data) - m.p*( D + 1 ) - 2 )
m.Σ[:,:] = Symmetric( tmp )
return
end
function fit!( m::ar, data )
N = size(data)[1]
Φ = zeros( N-m.p, m.p+1 )
# 定数項
Φ[:,1] .+= 1
for i in 0 : m.p-1
Φ[ :,i+2 ] += data[ m.p-i : end-1-i ]
end
Y = data[ m.p+1 : end ]
m.params[:] = inv(Φ'Φ)Φ'Y
m.σ = mean( ( Y .- (m.params'Φ')' ).^2 )
return
end
function predict( m::var, data ; length::Int=1 )
D = size(data)[2]
# 時系列データは直近のデータを用いて予測を行う。なのでデータを反転させるためのコピーが必要になる
tmp_data = data
for i in 1:length
rev_data = reverse( tmp_data', dims=2 )
basis_func = [1 ; rev_data[begin:m.p*D]]
pre = m.params'basis_func + rand( MvNormal( zeros(D), Symmetric( m.Σ ) ), 1 )
# 観測値と予測値をつなげて返す
tmp_data = cat( tmp_data, pre', dims=1 )
end
return tmp_data
end
# 条件付き期待値
function conditional_expectation( m::var, data ; length_pre::Int=1 )
D = size(data)[2]
tmp_data = data
for i in 1:length_pre
rev_data = reverse( tmp_data', dims=2 )
basis_func = [1 ; rev_data[begin:m.p*D] ]
pre = m.params'basis_func
tmp_data = cat( tmp_data, pre', dims=1 )
end
# 条件付き期待値だけが欲しい
return tmp_data[end,:]
end
function conditional_expectation( m::ar, data ; length_pre::Int=1 )
tmp_data = data
for i in 1:length_pre
rev_data = tmp_data[end:-1:1]
basis_func =[ 1 ; rev_data[begin:m.p] ]
pre = dot( m.params, basis_func )
push!( tmp_data, pre )
end
return tmp_data[end]
end
function granger( data, caused, causing ; lags::Int=1 )
N, D = size(data)
# 帰無仮説のモデル--制約あり
# 2変数の場合、制約を課すとARモデルを推定することと同じになる
if D == 2
h_0 = ar( lags, rand( lags + 1 ), randn(1) )
else
h_0 = var( lags, rand( ( D - 1 )*lags + 1, D - 1 ), rand( D - 1, D - 1 ) )
end
# 対立仮設のモデル--制約なし
h_1 = var( lags, rand( D*lags + 1, D ), rand( D, D ) )
#=
e.g. P = 2, dims = 2, w11 = w12 = 1
params => [ w11 w12
w21 w22
w31 w32
w41 w42
w51 w52 ]
=#
# 影響を与えると考えられる変数の列を削除する
data_pan = data[:, setdiff(begin:end, causing)]
fit!( h_0, data_pan )
fit!( h_1, data )
# F検定統計量
ssr0 = ssr( h_0, data_pan, caused )
ssr1 = ssr( h_1, data, caused )
# 制約の数
r = lags
tmp_up = ( ssr0 - ssr1 ) / r
tmp_down = ssr1 / ( N - D*lags - 1 )
F = tmp_up / tmp_down
# カイ二条分布
x² = Chisq( r )
# 95%点
crit_val = quantile( x², 0.95 )
println( "h_0：", h_0.params )
println( "h_1：", h_1.params )
println( "critical value：", crit_val )
println( "statistic value：", F )
print( "conclusion: " )
if crit_val < r*F
println( "exists granger causality" )
else
println( "not exists granger causality" )
end
end
# 残差平方和
# sum of squared residuals
function ssr( m::var, data, caused )
N = size(data)[1]
ssr = 0
for i in 1 + m.p : N
ŷ = conditional_expectation( m, data[begin : i - 1, :], length_pre=1 )[caused]
y = data[i, caused]
ssr += ( y - ŷ ) ^2
end
ssr
end
function ssr( m::ar, data, caused )
N = size(data)[1]
ssr = 0
for i in 1 + m.p : N
ŷ = conditional_expectation( m, data[begin : i - 1, caused], length_pre=1 )
y = data[i, caused]
ssr += ( y - ŷ ) ^2
end
ssr
end
使い方
サンプルを用意しました。Xには適当な時系列データを入れてください。
# size( X ) = ( 100, 3 )
# 3変数var(2)
params = var( 2, rand(2, 3), rand( 2, 2 ) )
# 学習
fit!( params, X )
# 2期先まで予測
predict( params, X, length=2 )
# 2->1 のグレンジャー因果性検定
granger( X, 1, 2, lags=1 )