渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

偏Granger因果で「第三者効果」を排除しつつ因果性検定してみる

遥か古の時代、まだ自分が研究者だった頃にデータ分析に使っていた手法のひとつに偏Granger因果 (partial Granger causality) というものがありました。これはGuo et al. (2008)で提唱されたもので、当時は著者グループ提供のオリジナルMatlabツールボックスを使っていたのですが、仕事も変わりRやPythonをメインに使うようになってからは触る機会は全くなくなっていたのでした。


ところが、先日偶然その偏Granger因果について触れる機会があったので、もしかしてと思ってググってみたらR実装があることに気付きました。こちらの{FIAR}パッケージです。

とりあえずちょっと触ってみた感じでは、当時のMatlabツールボックスよりもGUIという面では弱いものの、一方でMatlabツールボックスでは未整備でユーザーが自分で書かなきゃいけなかったようなところは逆に整備されているということで、一旦は実用に耐えそうだという印象です。


ということで、今回の記事では偏Granger因果について説明しつつ、その{FIAR}パッケージを用いた実践方法について簡単に紹介してみることにします。なお、ここでは自分で延々とTeXで式を書くのは極めてだるい(どれも割と冗長な表現がされていて式の長さの割に大したことは言ってない)ので、個々の式に関しては最小限TeXで書くに留め、代わりに参考文献・資料へのポインタを示しておきます。詳しい数式展開が見たいという方はそれらをご自分でご覧ください。


偏相関(partial correlation)について


これは僕がまだ駆け出しだった頃の記事で取り上げたことがあります。

要は、多変量のデータセットの中で普通の単相関即ち1-by-1のピアソン積率相関係数を計算すると、相関を計算したA & Bの相関だけでなくその他の例えばCやDからの相関の影響を受けて歪んだ結果になるかもしれない。それを補正するためには偏相関という考え方を用いるべき、というお話です。その定義や中身については赤本こと『統計学入門』p.52-53の説明やWikipedia英語版の記事が分かりやすいかと思います。

相関係数とは、いま、変数1から変数3まで三つの変数があるとき、変数3の影響を除いたあとの変数1と変数2の間の相関係数のことで、一般に r_{12 \cdot 3}と書き


(3. 3)  r_{12 \cdot 3} = \frac{r_{12} - r_{13} r_{12}}{\sqrt{1 - r_{13}^2} \sqrt{1 - r_{23}^2}}


と定義される。

上記の駆け出しの頃の記事では「重回帰分析(多変量線形回帰モデル)の方がもっとエレガントに第3の変数からの影響を排除したパラメータを計算できる」という趣旨のことを書きましたが、全てが全てそのような計算で済むわけでもないので、偏相関のような考え方は今でも有意義だと思っています。


偏相関のアナロジーでGranger因果も扱う=偏Granger因果(partial Granger causality)


ここからが本題です。そもそもGranger因果というのは、以前にも書いたように沖本本p.81の2次VARモデルを用いた例によれば以下のように定義されます。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

Granger因果性検定の手順


 y_{1t} = c_1 + \phi^{(1)}_{11} y_{1,t-1} + \phi^{(1)}_{12} y_{2,t-1} + \phi^{(2)}_{11} y_{1,t-2} + \phi^{(2)}_{12} y_{2,t-2} + \epsilon_{1t}


なる2変量VAR(2)モデルをOLSで推定し、その残差平方和をSSR_1とする。次に、 y_{2t}がない(Granger因果を与え得る時系列を伴わない)ARモデル


 y_{1t} = c_1 + \phi^{(1)}_{11} y_{1,t-1} + \phi^{(2)}_{11} y_{1,t-2} + \epsilon_{1t}


をOLSで推定し、その残差平方和をSSR_0とする。この時F統計量は


F \equiv \frac{(SSR_0 - SSR_1) / 2}{SSR_1/(T-5)}


で定義され、2Fは漸近的に\chi^2(2)に従うことが知られている。従って2Fの値を\chi^2(2)の95%点を比較して2Fの方が大きければ、y_{2t}からy_{1t}へのGranger因果性が存在しないという帰無仮説を棄却し、y_{2t}y_{1t}の将来を予測するのに有用であると結論することになる。(※Tはサンプルサイズ)

Granger因果は基本的に二変量(bivariate)を前提としたものです*1。一方これを多変量(multivariate)に拡張させることを考えると、先に述べたようにA & Bの因果関係に対してその他の既知のCやDからの因果的影響(第三者効果)がある場合、それ単体では本質的には第三者効果を排除もしくは均衡させることができません。そこで考え出されたのが偏Granger因果です。


最初に提唱したGuo et. al (2008)のpp. 81-83を見るとその導出過程が出ていますが、全てTeXで書き出すのは壮絶に面倒臭い&沖本 / Hamiltonの数式の書き方と異なり読みにくいので、本文中の(8)式と(9)式だけを並べておきます。興味のある方はGoogle Scholarで検索して本文をお読みください*2。(8)が偏Granger因果のF統計量、(9)が従来のGranger因果のF統計量をGuo論文中の定義に合わせて表したものです。

 F_1 = ln \left( \frac{|R^{(1)}_{XX|Z}|}{|R^{(2)}_{XX|Z}|}  \right) = ln \left( \frac{ S_{11} - S_{12} S_{22}^{-1} S_{21} }{ \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21} }  \right) (8)

 F_2 = ln \left(  \frac{|S_{11}|}{|\Sigma_{11}|} \right) (9)

添字だけ見ても分かりにくいかもですが、一応(8)式は3つの時系列 X_t, Y_t, Z_t同士の関係性を表しています( Zが外生的影響)。これらは S \Sigmaが出てきているのを見れば分かるように共分散行列を用いて個々のVARモデルの予測誤差分散を表現している式ですが、その詳細はTeXで長々と書き写すのが面倒なのでここでは割愛します。その代わり、この2つの式の直前にある説明を見てみましょう。

The value of  R^{(1)}_{XX|Z} measures the accuracy of the autoregressive prediction of X based on its previous values conditioned on Z by eliminating the influence of the common exogenous input and latent variables, whereas the value of  R^{(2)}_{XX|Z} represents the accuracy of predicting present value of X based on the previous history of both X and Y conditioned on Z by eliminating the influence of the exogenous input and latent variables. According to the causality definition of Granger, if the prediction of one process is improved by incorporating the information of the second process, then the second process causes the first process.

つまり、あらかじめ Z Xへの影響を考慮した Xの自己予測(予測誤差分散)に対して、 Z X Yへの影響を考慮した Xがその自己予測を改善する(予測誤差分散を小さくする)のであれば、偏Granger因果の意味で X Yとの間には因果性があると認める(なので前者と後者の比を取ってその比が大きくなればなるほど因果性が強いと言える)、というロジックになっています。これが(8)式で、 Zの影響を入れなければ(9)式になるというわけです。


一方偏Granger因果のF値は通常のGranger因果のそれとは異なり、経験的にもどのような確率分布に従うかが分かっていません。このため、有意性を判定するためにはブートストラップ法を用いる必要がある点に注意が必要です。


なお本文中にはカーネル自己回帰を使った非線形VARへの拡張という話題も出ていましたが、そもそも非線形VARという考え方自体がchallengingというかcrazyな気がするのでここでは割愛します。それに、例えばトレンド過程でやりたければ普通に差分系列でやるのが筋だと思いますので。。。他にも共和分の時は差分系列VAR自体が成立しないとか様々な制約条件があるわけですが、論文でもこの後紹介するパッケージでもそこまで考慮されているとは到底思えない*3のでこちらも割愛します。


{FIAR}パッケージで偏Granger因果検定をやってみる


では試しにやってみましょう。沖本本のサンプルデータであるMSCIを使います。

因果関係が明確な場合

手始めに、ベースラインとして{vars}パッケージでVARモデルを推定し、試しにUSの各国への影響をインパルス応答関数推定で見ておきます。ここでは沖本本p. 91の3番目の推奨に従い、us - uk - jpの順に並べることとします。

> library(vars)
> d <- read.csv('msci.txt', sep = '\t')
> d <- d[,-1]
> d1 <- d[, c(7, 6, 5)]
> names(d1)
[1] "us" "uk" "jp"
> VARselect(d1)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      3      2      3 

$criteria
                  1            2            3            4            5            6
AIC(n) 1.652867e+01 1.612393e+01 1.609777e+01 1.609881e+01 1.610355e+01 1.610417e+01
HQ(n)  1.654568e+01 1.615369e+01 1.614028e+01 1.615407e+01 1.617156e+01 1.618493e+01
SC(n)  1.657412e+01 1.620347e+01 1.621140e+01 1.624652e+01 1.628535e+01 1.632006e+01
FPE(n) 1.507690e+07 1.005856e+07 9.798843e+06 9.809018e+06 9.855645e+06 9.861771e+06
                  7            8            9           10
AIC(n) 1.610701e+01 1.611700e+01 1.612183e+01 1.612373e+01
HQ(n)  1.620052e+01 1.622327e+01 1.624084e+01 1.625550e+01
SC(n)  1.635698e+01 1.640107e+01 1.643998e+01 1.647597e+01
FPE(n) 9.889822e+06 9.989214e+06 1.003759e+07 1.005674e+07

> d1.var <- VAR(d1, p = 3)
> d1.irf <- irf(d1.var, impulse = 'us', n.ahead = 30, ci = 0.95, boot = T)
> plot(d1.irf, lwd = 2)

f:id:TJO:20180911150136p:plain

大体こんな感じになるはずです。次に、causality関数で通常のGranger因果を見てみましょう。

> causality(d1.var, cause = 'us', boot = T, boot.runs = 1000)
$Granger

	Granger causality H0: us do not Granger-cause uk jp

data:  VAR object d1.var
F-Test = 55.94, boot.runs = 1000, p-value < 2.2e-16


$Instant

	H0: No instantaneous causality between: us and uk jp

data:  VAR object d1.var
Chi-squared = 271.69, df = 2, p-value < 2.2e-16

インパルス応答で見たのと同じように、USの影響は各国に対して強く出ているということが見て取れます。では、今度は{FIAR}パッケージの出番です。やることはシンプルで、partGranger関数で偏Granger因果検定を行うだけです。Matlabツールボックス時代と違って、ブートストラップ法によるp値算出機能が初めからついているので楽です。

> partGranger(d1, nx = 1, ny = 2, order = 3, perm = T, prob = T, bs = 1000)
$orig
[1] 0.1131896

$prob
[1] 0

当たり前ですが、USの因果的影響が有意に強いという結果になりました。

因果関係が不明確な場合

では、これを例えば外生性という点で見ると逆になるjp - uk - usという順に並び替えたらどうなるでしょうか?

> d2 <- d[, 5:7]
> names(d2)
[1] "jp" "uk" "us"
> VARselect(d2)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      3      2      3 

$criteria
                  1            2            3            4            5            6
AIC(n) 1.652867e+01 1.612393e+01 1.609777e+01 1.609881e+01 1.610355e+01 1.610417e+01
HQ(n)  1.654568e+01 1.615369e+01 1.614028e+01 1.615407e+01 1.617156e+01 1.618493e+01
SC(n)  1.657412e+01 1.620347e+01 1.621140e+01 1.624652e+01 1.628535e+01 1.632006e+01
FPE(n) 1.507690e+07 1.005856e+07 9.798843e+06 9.809018e+06 9.855645e+06 9.861771e+06
                  7            8            9           10
AIC(n) 1.610701e+01 1.611700e+01 1.612183e+01 1.612373e+01
HQ(n)  1.620052e+01 1.622327e+01 1.624084e+01 1.625550e+01
SC(n)  1.635698e+01 1.640107e+01 1.643998e+01 1.647597e+01
FPE(n) 9.889822e+06 9.989214e+06 1.003759e+07 1.005674e+07

> d2.var <- VAR(d2, p = 3)
> d2.irf <- irf(d2.var, impulse = 'jp', n.ahead = 30, ci = 0.95, boot = T)
> plot(d2.irf, lwd = 2)
> causality(d2.var, cause = 'jp', boot = T, boot.runs = 1000)
$Granger

	Granger causality H0: jp do not Granger-cause uk us

data:  VAR object d2.var
F-Test = 0.76183, boot.runs = 1000, p-value = 0.759


$Instant

	H0: No instantaneous causality between: jp and uk us

data:  VAR object d2.var
Chi-squared = 37.704, df = 2, p-value = 6.496e-09


> partGranger(d2, nx = 1, ny = 2, order = 3, perm = T, prob = T, bs = 1000)
$orig
[1] 0.001248189

$prob
[1] 0.098

f:id:TJO:20180911150847p:plain

インパルス応答関数で見ても、通常のGranger因果でも、偏Granger因果でも、JPからの因果的影響は有意ではないという結果になりました。が、偏Granger因果のp値が大幅に低くなっているのがちょっと気になりますね。USからUKへの第三者効果としての因果的影響を引いた結果として、JPからUKへの予測誤差分散が小さくなったということなのかもしれません。。。


なお、この手法を探索的(F値に正負の別をつけて因果の「向き」も推定する)に適用する"difference Granger causality"と称する方法論もあり、{FIAR}パッケージではそれぞれdiffGranger, pdiffGrangerという関数で実装されていますが、最初にインパルス応答関数などで当たりをつける方が妥当だと個人的には考えているのでここでは取り上げません。


最後に


{FIAR}パッケージは"Functional Integration Analysis in R"というフルネームが示すように、僕にとっては昔取った杵柄でもある機能的MRIデータの時系列分析に用いられるパッケージなのですが、偏Granger因果以外の手法*4については正直ビジネス系の実務で使う場面が思いつかない&その信頼性については個人的に懸念があるので、今回の記事では全面的に割愛しました。興味のある方はご自身でお調べください。ということで、お後がよろしいようで。

*1:ただし沖本本pp. 81-82にあるようにn変量VARモデルでは検定統計量をrFと変更し(rは因果性検定に必要な制約の数)、このrFが自由度rのカイ二乗分布に従うことから有意性を判定する方法がある

*2:PDFがどこかで公開されています笑

*3:実際論文本文にはunit rootやcointegrationという単語が一語も出てこない

*4:例えばDynamic Causal Modelingとか