ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【時系列解析】グレンジャー因果 (Granger causality)の周波数領域表現

グレンジャー因果 (Granger causality)の評価についてのGewekeのアプローチを理解するために,2変量時系列でその方法をまとめます.2よりも多い変量だと私には難しいので,まず,2変量で考えます.これは,自分の理解を深めるためのメモです.Gewekeの論文はこれです.
https://www.tandfonline.com/doi/abs/10.1080/01621459.1982.10477803
さらに,この論文も参考にしました.
https://doi.org/10.1016/j.jneumeth.2005.06.011

線形因果の周波数領域での評価

 2変量の弱定常な確率過程 \left\{(x_t, y_t)\right\}を考えます.

 [解説]弱定常」というのは,アンサンブルの期待値の意味で,平均,分散,自己共分散 (自己相関)関数が,時間シフトに対して不変ということで,自己共分散 (自己相関)関数は,時間差 (ラグ)の関数で与えられます.

 \left\{(x_t, y_t)\right\}が無限次の移動平均過程で記述できる線形過程とします.つまり,ラグオペレータL多項式

 \begin{eqnarray} A^{(11)}(L) &=& a_1^{(11)} L + a_2^{(11)} L ^2 + a_3^{(11)} L ^3 + \cdots \\
A^{(12)}(L) &=& a_1^{(12)} L + a_2^{(12)} L ^2 + a_3^{(12)} L ^3 + \cdots \\
A^{(21)}(L) &=& a_1^{(21)} L + a_2^{(21)} L ^2 + a_3^{(21)} L ^3 + \cdots \\
A^{(22)}(L) &=& a_1^{(22)} L + a_2^{(22)} L ^2 + a_3^{(22)} L ^3 + \cdots \\
\end{eqnarray}

を使って.

  \begin{eqnarray} x_t &=& A^{(11)}(L)\, \varepsilon^{(x)}_{t} + A^{(12)}(L)\, \varepsilon^{(y)}_{t}\\
y_t &=& A^{(21)}(L)\, \varepsilon^{(x)}_{t} + A^{(22)}(L)\, \varepsilon^{(y)}_{t} \end{eqnarray}
 
と書けると言うことです.

 [解説] ラグオペレータは,L \varepsilon_{t} = \varepsilon_{t-1}のように時間を1つ戻す線形演算子です.L^j \varepsilon_{t} = \varepsilon_{t-j}です.Lは値を持つ変数ではありませんが,変数と同ように四則演算や指数法則が使えます.

 元の論文との対応が分かるように行列を使って書けば,

 \begin{bmatrix}
x_t \\
y_t \\
\end{bmatrix}  = \begin{bmatrix}
A^{(11)}(L)& A^{(12)}(L) \\
A^{(21)}(L)& A^{(22)}(L) \\
\end{bmatrix} \begin{bmatrix}
\varepsilon^{(x)}_{t} \\
\varepsilon^{(y)}_{t} \\
\end{bmatrix}

です.

 元の論文の流れに従い,移動平均過程から入りましたが,計算に使うのは自己回帰過程としての表現です.ということで,\left\{(x_t, y_t)\right\}は自己回帰過程として

  \begin{eqnarray} B^{(11)}(L) x_t + B^{(12)}(L) y_t &=&  \varepsilon^{(x)}_{t}\\
B^{(11)}(L) x_t + B^{(12)}(L) y_t &=& \varepsilon^{(y)}_{t} \end{eqnarray}

と書けるとします.ここで,B^{(11)}(L)B^{(12)}(L) B^{(21)}(L)B^{(22)}(L) は,L多項式です (自己回帰過程は無限次の移動平均過程で書けます).行列で表せば,

 \begin{bmatrix}
B^{(11)}(L)& B^{(12)}(L) \\
B^{(21)}(L)& B^{(22)}(L) \\
\end{bmatrix} \begin{bmatrix}
x_t \\
y_t \\
\end{bmatrix}  =  \begin{bmatrix}
\varepsilon^{(x)}_{t} \\
\varepsilon^{(y)}_{t} \\
\end{bmatrix}

です.

 パワースペクトルを計算するのに,\varepsilon^{(x)}_{t}\varepsilon^{(y)}_{t}が無相関のほうが見通しがいいので,

 \sigma_x^2 = {\rm var} \left( \varepsilon^{(x)}_{t} \right) = {\rm E} \left( \left(\varepsilon^{(x)}_{t}\right)^2 \right)
C = {\rm cov} \left(\varepsilon^{(x)}_{t}, \varepsilon^{(y)}_{t}\right) = {\rm E} \left(\varepsilon^{(x)}_{t} \varepsilon^{(y)}_{t}\right)

を使って,

 \begin{eqnarray} \begin{bmatrix}1 & 0 \\ {-} C /\sigma_x^2 & 1 \\
\end{bmatrix} \begin{bmatrix}
B^{(11)}(L)& B^{(12)}(L) \\
B^{(21)}(L)& B^{(22)}(L) \\
\end{bmatrix} \begin{bmatrix}
x_t \\
y_t \\
\end{bmatrix}  &=&  \begin{bmatrix} 1 & 0 \\ {-} C / \sigma_x^2 & 1 \\
\end{bmatrix} \begin{bmatrix}
\varepsilon^{(x)}_{t} \\
\varepsilon^{(y)}_{t} \\
\end{bmatrix} \\ 
\begin{bmatrix}
B^{(11)}(L)& B^{(12)}(L) \\
{-} \frac{C \, B^{(11)}(L)}{\sigma_x^2} + B^{(21)}(L)& {-} \frac{C\, B^{(12)}(L)}{\sigma_x^2}+ B^{(22)}(L) \\
\end{bmatrix} \begin{bmatrix}
x_t \\
y_t \\
\end{bmatrix}  &=&  \begin{bmatrix}
\varepsilon^{(x)}_{t} \\
{-} \frac{C\, \varepsilon^{(x)}_{t}}{\sigma_x^2} + \varepsilon^{(y)}_{t} \\
\end{bmatrix} \quad (*)
\end{eqnarray}

と変形します.

 [解説] \varepsilon^{(x)}_{t} {-} \frac{C\, \varepsilon^{(x)}_{t}}{\sigma_x^2} + \varepsilon^{(y)}_{t}の相関を計算すると,
 {\rm E} \left(\varepsilon^{(x)}_{t} \cdot \left( {-} \frac{C\, \varepsilon^{(x)}_{t}}{\sigma_x^2} + \varepsilon^{(y)}_{t}\right) \right)= {-} \frac{C\, {\rm E}\left(\left(\varepsilon^{(x)}_{t}\right)^2\right) }{\sigma_x^2}+{\rm E}\left(\varepsilon^{(x)}_{t} \varepsilon^{(y)}_{t}\right)={-} \frac{C \sigma_x^2}{\sigma_x^2}+C = 0
となり,無相関になってます.

 式の表現を簡単にするために,上の式 (*)

  \begin{eqnarray} \begin{bmatrix}
B^{(11)}(L)& B^{(12)}(L) \\
\tilde{B}^{(21)}(L)& \tilde{B}^{(22)}(L) \\
\end{bmatrix} \begin{bmatrix}
x_t \\
y_t \\
\end{bmatrix}  &=&  \begin{bmatrix}
\varepsilon^{(x)}_{t} \\
\tilde{ \varepsilon}^{(y)}_{t} \\
\end{bmatrix} 
\end{eqnarray}
と表します.

 この式を使って,x_tパワースペクトルを計算すると,

 \begin{eqnarray}
S_x (f) &=& \frac{1}{\Delta \Delta^{*}} \left(\tilde{B}^{(22)}(e^{i 2\pi f})\, \sigma_x^2 \, \tilde{B}^{(22)}(e^{- i 2\pi f})  +  {B}^{(12)}(e^{i 2\pi f})\, \tilde{\sigma}_y^2 \, {B}^{(12)}(e^{-i 2\pi f}) \right)
 \end{eqnarray}

になると思います (行列を使ったパワースペクトルの計算はまだ説明してません).ここで,

 \displaystyle \Delta = \frac{1}{B^{(11)}(e^{i 2\pi f}) \tilde{B}^{(22)}(e^{i 2\pi f}) - {B}^{(12)}(e^{i 2\pi f}) B^{(12)}(e^{i 2\pi f}) }

です.

 [解説] ラグオペレータLフーリエ変換すると,e^{-i 2\pi f}になります.

 ここで計算したS_x (f)では,後ろの項が\tilde{ \varepsilon}^{(y)}_{t}の影響,つまり,y_tの影響を表しています.

 ということで,Gewekeが導入した,y_tが原因としてx_tへ影響を与える線形因果の指標は

 \displaystyle \phi_{y \to x} (f) = \log \frac{\left|\tilde{B}^{(22)}(e^{i 2\pi f})\, \sigma_x^2 \, \tilde{B}^{(22)}(e^{-i 2\pi f})  +  {B}^{(12)}(e^{i 2\pi f})\, \tilde{\sigma}_y^2 \, {B}^{(12)}(e^{-i 2\pi f})\right|}{\left|\tilde{B}^{(22)}(e^{i 2\pi f})\, \sigma_x^2 \, \tilde{B}^{(22)}(e^{-i 2\pi f})\right|}

です.具体的に書こうとした分,式の変形や,タイプ結果にちょっと自信がないので,近いうちに数値実験で検証します.これは,十分に確認された結果はないので,信じないでください.間違っている部分は指摘してもらえると助かります.