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

グレンジャー因果 (Granger causality)の評価についてのGewekeのアプローチを理解するために，２変量時系列でその方法をまとめます．２よりも多い変量だと私には難しいので，まず，２変量で考えます．これは，自分の理解を深めるためのメモです．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}$ のように時間を１つ戻す線形演算子です． $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|}$