VARそして時系列因果性分析の復習 - 渋谷駅前で働くデータサイエンティストのブログ

f:id:TJO:20201123172036p:plain

「新型コロナウイルス感染症における治療の進展（令和2年10月29日に開催された第13回新型コロナウイルス感染症対策分科会事務局提出資料を基に内閣官房・内閣府作成）」という資料が世間で物議を醸しているようです。ただ、これを見ていて僕が個人的に気になったのは、その議論の内容や結論ではなく、「グレンジャー(Granger)因果」が使われているという点でした。

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

作者:竜義, 沖本
発売日: 2010/02/01
メディア: 単行本

Time Series Analysis

作者:Hamilton, James D.
発売日: 1994/01/11
メディア: ハードカバー

以前このブログでも一通り計量時系列分析を取り上げて一生懸命沖本本やHamiltonで勉強しながらシリーズ記事を書いたものですが、その時の記憶から言えば「Granger因果ってそんなに軽い気持ちで使って大丈夫な代物だったっけ？」という印象を強く受けたものです。そこで、今回の記事では以前のブログ記事から要所要所を抜粋しながら、Granger因果を含む時系列因果性分析がどんなものであるかを復習し、その上で冒頭に紹介した資料で行われている分析の妥当性について考察してみようと思います。

いつもながらですが、僕の計量時系列分析に関する知識は中途半端な点が多く記事中には誤りや理解不足が含まれている可能性がありますので*1、お気付きの方はどしどしコメント欄などでご指摘下さると有難いですm(_ _)m

VAR及び因果性分析の復習
問題の資料におけるGranger因果の適用について

VAR及び因果性分析の復習

この辺の話は過去にこのブログで連続シリーズとして散々取り上げた内容ですので、あくまでも復習ということで必要な箇所だけ適宜引用して紹介するに留めたいと思います。なお、これらの議論において本来参照すべき沖本本もHamiltonも勤務先オフィスに置いたままであり、コロナ禍で出社が禁じられている現状ではどちらも手に入りませんので、どちらの内容であれこのブログに掲載した際の情報を元に孫引きする形になる点予めご了承ください。

また、Rによる実践例は過去記事に全て収録してありますので、今回の記事では割愛します。大した手間ではありませんし、{vars}パッケージをインストールの上皆さんのお手元で是非実践されることをお薦めします。

VARモデル

まず、VAR（ベクトル自己回帰）モデルについて。要は多変量時系列同士の関係性をまとめた自己回帰モデルですが、以前のブログ記事から説明を以下に抜粋しておきます。

VAR(p)モデルを $\mathbf{y}_t$ を定数と自身のp期の過去の値に回帰したモデルとすると、
$\mathbf{y_t}=\mathbf{c}+\mathbf{\Phi_1}\mathbf{y_{t-1}}+\cdots+\mathbf{\Phi_p}\mathbf{y_{t-p}}+\mathbf{\epsilon_t}, \mathbf{\epsilon_t} \sim W.N.(\mathbf{\Sigma})$
（ただし $\mathbf{c}$ は $n \times 1$ 定数ベクトル、 $\mathbf{\Phi_i}$ は $n \times n$ 係数行列）
これでは分かりにくいので、例えば2変量VAR(1)モデルに書き下すとこんな感じになります。
$\left\{ \begin{array} ~y_{1t}=c_1+\phi_{11} y_{1,t-1}+\phi_{12} y_{2,t-1}+\epsilon_{1t} \\ y_{2t}=c_2+\phi_{21} y_{1,t-1}+\phi_{22} y_{2,t-1}+\epsilon_{2t} \end{array} \right.$
$\left ( \begin{array} ~ \epsilon_{1t} \\ \epsilon_{2t} \end{array} \right ) \sim W.N.(\Sigma)$
n変量VAR(p)モデルはn本の回帰式からなり、それぞれの回帰式は各変数を定数と全変数のp期間の過去の値に回帰した形になっています。

VARモデルの各方程式は同時点のその他の変数は含まないので、同時方程式モデル(simultaneous equation model)ではない*2とされます。ただし、各方程式は誤差項の相関を通じて関係しており、見かけ上無関係な回帰(SUR)モデル(seemingly unrelated regression model)と呼ばれるようです。

一般にSURモデルを推定するためには誤差項同士の相関を考慮に入れなければならず、全ての回帰式を同時に推定する必要があるそうですが、VARモデルは全ての回帰式が同一の説明変数を持つという特殊性があるため、実は各方程式を個別にOLSによって推定するだけで良いということが分かっています。しかも撹乱項 $\epsilon_t$ が多変量正規分布に従うと仮定した場合、OLS係数推定量は最尤推定量とも一致します。

ということで、VARモデルというのは基本的には多時系列自己回帰モデルです。なおARモデルに対してMAモデルやARMA / ARIMAモデルがあるように、VARモデルにもVARMA / VARIMAといったモデルがありますし、さらには季節調整を加味したモデルもありますが、ここでは割愛します。

Granger因果

さて、問題のGranger因果ですが。沖本本ではこのように定義が紹介されています。

定義4.1（グレンジャー因果性）現在と過去の $x$ の値だけに基づいた将来の $x$ の予測と、現在と過去の $x$ と $y$ の値に基づいた将来の $x$ の予測を比較して、後者のMSEの方が小さくなる場合、 $y_t$ から $x_t$ へのグレンジャー因果性(Granger causality)が存在するといわれる。

定義4.2（一般的なグレンジャー因果性） $\mathbf{x}_t$ と $\mathbf{y}_t$ をベクトル過程とする。また、 $x$ と $y$ の現在と過去の値を含む、時点 $t$ において利用可能な情報の集合を $\Omega_t$ とし、 $\Omega_t$ から現在と過去の $\mathbf{y}$ を取り除いたものを $\tilde{\Omega}$ とする。このとき、 $\tilde{\Omega}_t$ に基づいた将来の $\mathbf{x}$ の予測と、 $\Omega_t$ に基づいた将来の $\mathbf{x}$ の予測を比較して、後者のMSEのほうが小さくなる場合、 $\mathbf{y}_t$ から $\mathbf{x}_t$ へのグレンジャー因果性が存在するといわれる。ここで、MSEの大小は行列の意味での大小であることに注意されたい。

この定義に基づいて、実際にGranger因果性検定を行う流れが以下です。

n変量VAR(p)モデルにおけるグレンジャー因果性検定の手順
(1) VARモデルにおける $y_{kt}$ のモデルをOLSで推定し、その残差平方和を $SSR_1$ とする。
(2) VARモデルにおける $y_{kt}$ のモデルに制約を課したモデルをOLSで推定し、その残差平方和を $SSR_0$ とする。
(3) F統計量を
$F \equiv \frac{(SSR_0 - SSR_1)/r}{SSR_1 / (T-np-1)}$
で計算する。ここで、rはグレンジャー因果性検定に必要な制約の数である。
(4) $rF$ を $\chi^2(r)$ の95%点と比較し、 $rF$ のほうが大きければ、ある変数（群）から $y_{kt}$ へのグレンジャー因果性は存在し、小さければグレンジャー因果性は存在しないと結論する。

この説明を何となくそれっぽく模式図にしたのが、5年前の因果フェスで登壇した際に作った以下の図です。

f:id:TJO:20201126192037p:plain

こんな感じなのですが、過去にこのブログで書いた統計的因果推論のシリーズ記事をご記憶の方や、はたまた畏友・安井君が著した名著『効果検証入門』でモダンな統計的因果推論についてある程度知識をお持ちの方であれば、「これは因果推論とは言えない」ということがお分かりいただけるかと思います。

効果検証入門〜正しい比較のための因果推論／計量経済学の基礎

作者:安井翔太
発売日: 2019/12/27
メディア: Kindle版

というのは、モダンな統計的因果推論であれば一般には「反実仮想」(counterfactual)即ち「原因Xがなかった場合に結果Yはどうなるべきだったか」ということを「原因Xがあった場合に結果Yは実際にどうなったか」と比べるというプロセスを程度問題ながら伴いますが、Granger因果にはそれはありません。そこにあるのは「時間的ラグを置いて先行するXがあった場合にYの予測を改善するかどうか」という、モデルの予測性能に関する議論です。Wikipedia記事にも載っている例で言うと「雷が光ると雷鳴が轟くが、雷がピカッと光ること自体は雷鳴の原因ではない」わけで、当該記事の冒頭でも「グレンジャー因果性は、XがYを引き起こすかどうかを検定するのではなく、XによってYを予測できるかどうかを検定するものである」とはっきり断っています。

修論でVAR使ったとき、グレンジャー検定で政策の効果が統計的に確認されたと言ったら、「グレンジャーコーザリティは政策の因果効果を示すものじゃないから、◯◯がXXの正のインジケーターになってることが示唆された、くらいに言っておけ」と指導された思い出。
— 江口允崇 (@maseguchi) 2020年11月22日

この辺の議論に関連する、駒澤大学の江口先生のコメントがこちら。これは全くもってその通りだと個人的には思います。Granger因果そのものは何かのアクションの因果性を導くものではないはずなのです。

Granger因果と単位根過程＆見せかけの回帰との関係

そしてもう一点。Granger因果には、実は様々な制約事項があります。それは、「見せかけの回帰」が起きてしまうような単位根過程同士から成るVARモデルに対しては、そのままでは適用できないのです。単位根過程については以下に定義を孫引きしておきますが、要はランダムウォークなどと同じで「差分系列を取ると定常（平均回帰的）だがそれを累積系列として原系列に戻すと非定常（トレンドなど平均回帰しない）になる」時系列です。

定義5.1（単位根過程）　原系列 $y_t$ が非定常過程であり、差分系列 $\Delta y_t = y_t - y_{t-1}$ が定常過程であるとき、過程は単位根過程(unit root process)といわれる。

で、この単位根過程同士で単なる線形回帰を行うと、何故か偏回帰係数のp値が極端に下がって（t値が極端に上がって）、本来全く互いに無関係な時系列同士だったとして統計的有意であることになってしまうのが見せかけの回帰です。定義は以下の通り。

定義6.1（見せかけの回帰）　単位根過程 $y_t$ を定数と $y_t$ と関係のない単位根過程 $x_t$ に回帰すると、 $x_t$ と $y_t$ の間に有意な関係があり、回帰の説明力が高いように見える現象は見せかけの回帰(spurious regression)といわれる。

このような状況では、そのままVARモデルを推定してGranger因果を適用しようとしても各時系列が平均回帰的でないなどの理由で標準的な漸近理論が使えない＝F分布が使えない*3ということで、上記の通りF分布に従うことを前提とするGranger因果は使えなくなってしまう*4ことが知られています。これでは困ってしまいますね。

単位根過程同士で時系列因果性分析をやりたければどうするべきか

解決策は2つあって、1つは差分系列に直すこと。ところが、ここにも大きな落とし穴があります。単位根過程同士で2つ以上並べた際に時々見られる関係性として「共和分」というものがあります。

定義6.2（共和分）　 $x_t$ と $y_t$ を単位根( $I(1)$ )過程とする。 $ax_t+by_t$ が定常( $I(0)$ )過程となるようなaとbが存在するとき、 $x_t$ と $y_t$ との間には共和分(cointegration)の関係がある、もしくは $x_t$ と $y_t$ は共和分している(cointegrated)といわれる。より一般的には、 $\mathbf{y_t}$ を $I(1)$ とする。 $\mathbf{a'y_t}$ が $I(0)$ 過程となるような $\mathbf{a}$ が存在するとき、 $\mathbf{y_t}$ には共和分の関係がある、もしくは $\mathbf{y_t}$ は共和分しているといわれる。また、このとき、 $(a,b)'$ や $\mathbf{a}$ は共和分ベクトル(cointegrating vector)と呼ばれる。

これは平たく言うと「2つの時系列が互いに和を取ると一定になるなど均衡関係がある」という状況で、この場合は実は通常のVARモデルが使えません。この点に関するGrangerの表現定理を簡潔に僕がまとめたのが以下です。

共和分関係にある2つの単位根過程 $y_{1t}, y_{2t}$ を仮定した場合（これを共和分システムと仮に呼ぶ）、この二者を表すVARモデルも単位根を持つが、共和分システムの場合その差分系列のVMAモデルは単位根を持つため、反転不可能となる*5。

反転不可能な(V)MAモデルは(V)ARモデルでは表現不可能であるため、仮に共和分システムの差分系列に対してVARモデルを推定したとしても、それは本来ならば成立し得ない誤ったモデルとなる。

平たく言えば、共和分関係にある2つの時系列を差分系列VARモデルで正確に表現することは不可能だということです。つまり、共和分関係にある単位根過程同士ではGranger因果を適用すること自体がそもそも不可能なのです。ちなみに、共和分関係による影響を何とかして修正してVARモデルの代わりに多変量時系列同士の関係性をうまく表現するモデルがないか？ということで登場したのが、ベクトル誤差修正モデル(VECM: vector error correction model)です。

定理6.1（VAR(p)表現を持つ共和分システム $\mathbf{y}_t$ の性質）　VAR(p)表現を持つ共和分システム $\mathbf{y}_t$ は以下の性質をもつ。
（※ただしここで $\mathbf{y}_t = (y_{1t},y_{2t},\cdots,y_{kt})'$ とする。 $y_{kt}$ は個々の時系列を表す）

(1) VAR(p)表現は単位根を持つ。

(2) $\Delta \mathbf{y}_t$ のVMA表現は反転不可能である。

(3) $\Delta \mathbf{y}_t$ のVAR表現は存在しない。

(4) $\mathbf{y}_t$ は以下のVECM(p-1)で表現できる。

$\Delta \mathbf{y}_t = \zeta_1 \delta \mathbf{y}_{t-1} + \zeta_2 \delta \mathbf{y}_{t-2} + \cdots + \zeta_{p-1} \delta \mathbf{y}_{t-p+1} + \mathbf{\alpha} + \zeta_0 \delta \mathbf{y}_{t-1} + \epsilon_t$
$= \zeta_1 \delta \mathbf{y}_{t-1} + \zeta_2 \delta \mathbf{y}_{t-2} + \cdots + \zeta_{p-1} \delta \mathbf{y}_{t-p+1} + \mathbf{\alpha} - \mathbf{BA'} \mathbf{y}_{t-1} + \epsilon_t$
ここで、 $\mathbf{A, B}$ はn×h行列であり、 $\mathbf{A'y}_t \sim I(0)$ はh個の共和分関係を表す。また、共和分ランクhは行列 $\zeta_0$ のランクによって決まる。

(5) VECMには定常な変数しか含まれない。また、
$|\mathbf{I}_n - \zeta_1 z - \zeta_2 z^2 - \cdots - \zeta_{p-1} z^{p-1}| = 0$
のすべての解の絶対値は1より大きい。

このVECMであれば、共和分関係にある単位根過程同士の関係性を適切にモデリングすることができます。が、ここまで来るとそもそもGranger因果をどうやるかという問題はどこかに吹っ飛んでしまっています。実際、Rでcausality::varsやgrangertest::lmtestといった関数でGranger因果性検定をやりたくとも、VECMは適用外となっていてどうにもなりません*6。またしても困ってしまいました。

VAR/VECMにおけるインパルス応答推定

そこで出てくる解決策のもう1つが、インパルス応答の推定というアプローチです。これは過去記事からそのまま孫引きも含めて以下に抜粋しておきます。

まず、VARモデルの定義式を眺めてみましょう。
$\mathbf{y_t}=\mathbf{c}+\mathbf{\Phi_1}\mathbf{y_{t-1}}+\cdots+\mathbf{\Phi_p}\mathbf{y_{t-p}}+\mathbf{\epsilon_t}, \mathbf{\epsilon_t} \sim W.N.(\mathbf{\Sigma})$
（ただし $\mathbf{c}$ は $n \times 1$ 定数ベクトル、 $\mathbf{\Phi_i}$ は $n \times n$ 係数行列）
ここで $\mathbf{\Sigma}$ は対角行列ではないとします。この時、インパルス応答関数は以下のように定義されます。

定義4.3（非直交化インパルス応答関数） $y_{jt}$ の撹乱項 $\epsilon_{jt}$ だけに1単位（または1標準偏差）のショックを与えたときの $y_{i,t+k}$ の変化は、 $y_j$ のショックに対する $y_i$ のk期後の非直交化インパルス応答と呼ばれる。また、それをkの関数としてみたものは、 $y_j$ のショックに対する $y_i$ の非直交化インパルス応答関数といわれる。

定義4.4（直交化インパルス応答関数）撹乱項の分散共分散行列の三角分解を用いて、撹乱項を互いに相関する部分と無相関な部分に分解したとき、無相関な部分は直交化撹乱項といわれる。また、 $y_j$ の直交化撹乱項に1単位または1標準偏差のショックを与えたときの $y_i$ の直交化インパルス応答を時間の関数としてみたものは、 $y_j$ に対する $y_i$ の直交化インパルス応答関数といわれる。

理論的な定義を全部書いていくと死んでしまうので、要点だけ抜粋しておきます。Rで主に用いることになる直交化インパルス応答関数を求める際に重要になるのは、撹乱項 $\mathbf{\epsilon}_t$ の分散共分散行列 $\mathbf{\Sigma}$ です。この $\mathbf{\Sigma}$ は、
$\mathbf{\Sigma} = \mathbf{ADA'}$
（ただし $\mathbf{A}$ は対角成分が1に等しい下三角行列、 $\mathbf{D}$ は対角行列）
と三角分解することができ、このとき直交化撹乱項 $\mathbf{u}_t$ を
$\mathbf{u}_t = \mathbf{A}^{-1} \epsilon_t$
と定義し、各変数固有の変動を表すとみなすことができます。なお、 $\mathbf{u}_t$ は互いに無相関です（ただし $Var(\mathbf{u}_t)=\mathbf{D}$ ：証明は沖本本を参照のこと）。

$y_j$ のショックに対する $y_i$ のインパルス応答関数は、 $u_{jt}$ に1単位のショックを与えたときの $y_i$ の反応を時間の関数として見たものであり、
$IRF_{ij}(k) = \frac{\partial y_{i,t+k}}{\partial u_{jt}}$
$k = 0,1,2,\cdots$
と表されます。1単位ではなく1標準偏差のショックを与えたいときは、上の式にそのまま $\sqrt{Var(u_{jt})}$ を乗じれば良いだけです。

なのですが、Rでは三角分解の代わりにコレスキー分解を用いて撹乱項を分解し、その分解した撹乱項に1単位のショックを与えてインパルス応答関数を計算するという方法を用いています。これについてちょっと書いておきましょう。 $\mathbf{D}^{1/2}$ を(j,j)成分が $u_{jt}$ の標準偏差に等しい対角行列であるとすると、
$\mathbf{\Sigma} = \mathbf{ADA'}$
は、 $\mathbf{P} \equiv \mathbf{AD}^{1/2}$ として
$\mathbf{\Sigma} = \mathbf{AD}^{1/2} \mathbf{D}^{1/2} \mathbf{A'} = \mathbf{PP'}$
と書けます。この式は正定値行列 $\mathbf{\Sigma}$ のコレスキー分解になっていて、 $\mathbf{P}$ はコレスキー因子と呼ばれます。この $\mathbf{P}$ を用いれば、撹乱項を
$\mathbf{v}_t = \mathbf{P}^{-1}\mathbf{\epsilon}_t = \mathbf{D}^{-1/2}\mathbf{A}^{-1}\mathbf{\epsilon}_t = \mathbf{D}^{-1/2}\mathbf{u}_t$
と分解することもできます。なお $Var(\mathbf{u}_t)=\mathbf{D}$ であることから、 $v_{jt}$ は $u_{jt}$ をその標準偏差 $\sqrt{Var(u_{jt})}$ で割ったものになっています。これは言い換えると $v_{jt}$ における1単位の増加が $u_{jt}$ における1標準偏差の増加に等しいということであり、同じ要領で言えば $v_{jt}$ に1単位のショックを与えて計算したインパルス応答関数は、 $u_{jt}$ に1標準偏差のショックを与えて計算したインパルス応答関数に等しくなるということです。

Rであればirf::varsで非常に簡単に実践することができます。またVECMはインパルス応答を推定する上ではVARモデルと全く同じように扱うことができるため、irf::varsはVECMにも対応しています。よって、インパルス応答を推定するというのがGranger因果の代替策としてはベストであろうと個人的には考えます。

問題の資料におけるGranger因果の適用について

ここまで書けばもうお分かりかと思いますが、要は

そもそもGranger因果はモダンな意味での因果関係の有無を検証するものではない
単位根過程（トレンドのある時系列）同士でGranger因果を検定する際は非常に多くの制約がある

ということなので、個人的には冒頭で紹介した資料に掲載されているような「うねりがありトレンドが強く乗っていて単位根過程のように見える」時系列データ同士に対してGranger因果を適用するのも不適切であれば、そこから「両者には因果関係はなかった」と結論づけるのも不適切だ、と感じた次第です。せめてインパルス応答推定ぐらいすれば良かったのに、というのが正直な感想です。勿論、インパルス応答とて基本的には「時系列ラグの前後関係に基づく見かけ上の因果性」を可視化するものに過ぎませんが、そもそも理論的に成立しない可能性のあるGranger因果を使うよりはだいぶマシなのではないかと思います。

もっと言えば、VAR/VECMを問わずインパルス応答推定は比較的容易に計算できる上に、過去に見かけた学術資料の中には「単位根過程の原系列同士のVARモデルであっても推定されたインパルス応答の精度は悪くない（つまり必ずしもVECMや差分VARモデルなどをわざわざ検討する必要はない）」と指摘するものもあったので、制約だらけのGranger因果なんて最初からやらずにいきなりインパルス応答推定でも良かったのではないでしょうか。

……と色々書きましたが、生データが入手できないのでここに書かれていることの大半は冒頭の資料1枚を見ただけで思いついた感想に過ぎないことを最後に付記しておきます。もっともGranger因果に対する制約の相対的な緩さを考えたら「インパルス応答で良いのでは」という結論は変わらないと思いますが（笑）。そして、そんなことをするくらいなら、どれほど緩くても良いので政策的介入の有無によるDIDテストを行ってCausalImpactで決着をつけた方がより確実だったのではないかと思います。

*1:そもそも一生懸命勉強していたのがもう7年も前なので色々と記憶が曖昧になってきている

*2:同時方程式については例えば http://upo-net.ouj.ac.jp/tokei/contents/sub_contents/c01_05_00.xml 辺りを参照のこと

*3:ここは確率過程論に僕が詳しくないので完全に沖本本の受け売りになっています、ごめんなさい

*4:正しい統計量を返さない

*5:特性方程式の全ての解の絶対値が1より大きいという条件を満たさない（単位根過程なので絶対値が1に等しい解を持つ）

*6:VARにしか対応していない