【Rで時系列解析】グレンジャー因果 (Granger causality)の周波数領域表現を数値実験で検証

グレンジャー因果 (Granger causality)のGewekeのアプローチを少しずつ理解していきます．前回の説明の続きで，今回は数値実験をします．
chaos-kiyono.hatenablog.com

1. 問題設定

　単純な例から考えはじめないと，私には理解できません．一般的な表現を見ても，私にはすぐに理解できません．簡単な例で試行錯誤してやっと理解できることがほとんどです．ということで，２変量の1次自己回帰過程を考えます．

　 $\begin{equation} \left\{\begin{array}{l} x_t = a_{11} \, x_{t-1} + a_{12} \, y_{t-1} + \varepsilon^{(x)}_{t} \\ y_t = a_{21} \, x_{t-1} + a_{22} \, y_{t-1} + \varepsilon^{(y)}_{t} \\ \end{array}\right. \end{equation}$

ここで， $\varepsilon^{(x)}_{t}$ ， $\varepsilon^{(x)}_{t}$ は，それぞれ，分散 $\sigma_x^2$ ， $\sigma_y^2$ の白色ノイズで，お互いに相関がないとします．

　この過程で， $a_{12} \neq 0$ であれば， $x$ の予測に，過去に観測した $y$ が役立つということなので，「 $y$ から $x$ へのグレンジャー因果性あり」です．また， $a_{21} \neq 0$ であれば，「 $x$ から $y$ へのグレンジャー因果性あり」です．

　実際の時系列解析で，係数 $\{a_{ij}\}$ をどうやって決めるの？と思うかもしれません．実際の時系列解析では，最小２乗法などで係数を決められるそうです．

2. 計算内容

　前回，詳しく説明しようと行列の要素を書いてみましたが，結局，式が見にくく，計算の見通しが悪くなっただけかもしれません．今回の計算の流れを式で書いておきます．

をラグオペレータ $L$ を使って，行列で表現すると，

　 $\begin{eqnarray} \begin{bmatrix} 1 - a_{11} \, L & - a_{12} \, L \\ {-} a_{21} \, L & 1- a_{22} \, L \\ \end{bmatrix} \begin{bmatrix} x_t \\ y_t \\ \end{bmatrix} &=& \begin{bmatrix} \varepsilon^{(x)}_{t} \\ \varepsilon^{(y)}_{t} \\ \end{bmatrix} \end{eqnarray}$

です．この式の左辺の行列を

　 $\begin{eqnarray} A &=& \begin{bmatrix} 1 - a_{11} \, L & - a_{12} \, L \\ {-} a_{21} \, L & 1- a_{22} \, L \\ \end{bmatrix} \end{eqnarray}$

として，逆行列 $A^{-1}$ を計算します (下のRスクリプトではBです)．

　逆行列の成分を， $\{\bar{a}_{ij}(L)\}$ とすれば， $x_t$ のパワースペクトルは ( $L$ に $e^{-i 2 \pi f}$ を代入して)，

　 $\displaystyle S_x(f) = \bar{a}_{11}\!\left(e^{-i 2 \pi f} \right) \, \bar{a}_{11}\! \left(e^{i 2 \pi f} \right) \sigma_x^2 + \bar{a}_{12}\! \left(e^{-i 2 \pi f} \right) \, \bar{a}_{12}\! \left(e^{i 2 \pi f} \right) \sigma_y^2$

となります．この式の第２項が，過去の $y$ が $x$ の予測に役立つ程度を表しています．

　ということで， $y$ から $x$ へのグレンジャー因果の指標は，

　 $\displaystyle \phi_{y \to x} (f) = \log \frac{S_x(f)}{\bar{a}_{11}\!\left(e^{-i 2 \pi f} \right) \, \bar{a}_{11}\! \left(e^{i 2 \pi f} \right) \sigma_x^2 }$

です (たぶん)．分子分母とも実数になるので，絶対値は付けてません．Rスクリプトでは，複素数表現 $a+0 i$ を実数表現 $a$ に変えるため，absとか，Reとかを使っています．

　 $x$ から $y$ へのグレンジャー因果の指標は，

　 $\displaystyle \phi_{x \to y} (f) = \log \frac{\bar{a}_{21}\!\left(e^{-i 2 \pi f} \right) \, \bar{a}_{21}\! \left(e^{i 2 \pi f} \right) \sigma_x^2 + \bar{a}_{22}\! \left(e^{-i 2 \pi f} \right) \, \bar{a}_{22}\! \left(e^{i 2 \pi f} \right) \sigma_y^2 }{\bar{a}_{22}\!\left(e^{-i 2 \pi f} \right) \, \bar{a}_{22}\! \left(e^{i 2 \pi f} \right) \sigma_y^2 }$

です (おそらく)．

3. 数値実験結果 (時系列は生成してません)

　係数 $\{a_{ij}\}$ を与えて計算した結果の図を示していきます．使ったRスクリプトは下の方にあります．ここでは，「コヒーレンス」と呼ばれる周波数成分の相関の程度を表す指標も計算しました．コヒーレンスについて，今回は説明しません．

yからxへのグレンジャー因果性ありの例
　 $\begin{equation} \left\{\begin{array}{l} x_t = 0.8 \, x_{t-1} + 0.4 \, y_{t-1} + \varepsilon^{(x)}_{t} \\ y_t = 0 \, x_{t-1} +0.6 a_{22} \, y_{t-1} + \varepsilon^{(y)}_{t} \\ \end{array}\right. \end{equation}$

ここで， $\varepsilon^{(x)}_{t}$ ， $\varepsilon^{(x)}_{t}$ は，分散１の白色ノイズで，お互いに相関がないとします．この例では，yからxへのグレンジャー因果性はありますが，xからyへのグレンジャー因果性はありません．

　下の図では，左の $\phi_{y \to x} (f)$ は，0より大きな値が見られますが，真ん中の $\phi_{x \to y} (f)$ は，ずっと0です．0なら，グレンジャー因果性なしということです．

xからyへのグレンジャー因果性ありの例
　 $\begin{equation} \left\{\begin{array}{l} x_t = 0.8 \, x_{t-1} + 0 \, y_{t-1} + \varepsilon^{(x)}_{t} \\ y_t = 0.2 \, x_{t-1} +0.6 a_{22} \, y_{t-1} + \varepsilon^{(y)}_{t} \\ \end{array}\right. \end{equation}$

ここで， $\varepsilon^{(x)}_{t}$ ， $\varepsilon^{(x)}_{t}$ は，分散１の白色ノイズで，お互いに相関がないとします．この例では，xからyへのグレンジャー因果性はありますが，yからxへのグレンジャー因果性はありません．

　下の図では，左の $\phi_{y \to x} (f)$ は，ずっと0ですが，真ん中の $\phi_{x \to y} (f)$ は，0よりも大きな値になっています．コヒーレンスは，因果性の向きに関係なく，相関があることを示しています．

xからyへも，yからxへもグレンジャー因果性ありの例
　 $\begin{equation} \left\{\begin{array}{l} x_t = 0.8 \, x_{t-1} + 0.4 \, y_{t-1} + \varepsilon^{(x)}_{t} \\ y_t = 0.2 \, x_{t-1} +0.6 a_{22} \, y_{t-1} + \varepsilon^{(y)}_{t} \\ \end{array}\right. \end{equation}$

ここで， $\varepsilon^{(x)}_{t}$ ， $\varepsilon^{(x)}_{t}$ は，分散１の白色ノイズで，お互いに相関がないとします．

　結果は下の図です．どちらの向きへもグレンジャー因果性があるということです．

Rスクリプト

　自己回帰過程 (AR過程)の次数が上がったときにも計算できるように改良してみてください．

##################
# ２変量(x,y)の1次ARの係数
# x[t] = a.11*x[t-1]+a.12*y[t-1] + eps.1[t]
# y[t] = a.21*x[t-1]+a.22*y[t-1] + eps.2[t]
#####
a.11 <- 0.8
a.12 <- 0.4
a.21 <- 0.0
a.22 <- 0.6
#####
# 白色ノイズeps.1, eps.2の分散
sig2.1 <- 1
sig2.2 <- 1
##################
# グレンジャー因果 (Granger causality)の周波数領域表現の計算
# AR係数の行列
A <- matrix(c(a.11,a.12,a.21,a.22),nrow=2,ncol=2,byrow=TRUE)
# 単位行列
I <- diag(2)
###
# 周波数の点数(描画のためのfを何点とるか)
n <- 500
f.list <- seq(0,0.5,length.out=n)
###
Sx <- c()
Sy <- c()
Sx0 <- c()
Sy0 <- c()
Sxy <- c()
for(j in 1:n){
 f <- f.list[j]
 z <- exp(2*pi*f*1i)
 # 逆行列
 B <- solve(I-A*z)
 # スペクトルの計算
 Sx[j] <- B[1,1]*Conj(B[1,1])*sig2.1+B[1,2]*Conj(B[1,2])*sig2.2
 Sx0[j] <- B[1,1]*Conj(B[1,1])*sig2.1
 Sy[j] <- B[2,1]*Conj(B[2,1])*sig2.1+B[2,2]*Conj(B[2,2])*sig2.2
 Sy0[j] <- B[2,2]*Conj(B[2,2])*sig2.2
 Sxy[j] <- B[1,1]*Conj(B[2,1])*sig2.1+B[1,2]*Conj(B[2,2])*sig2.2
}
################
# 線形因果の指標
fy2x <- log(abs(Sx)/abs(Sx0))
fx2y <- log(abs(Sy)/abs(Sy0))
# 描画
y.max <- max(1,fy2x,fx2y)
par(mfrow=c(1,3),cex.axis=1.2,cex.lab=1.5,las=1)
plot(f.list,fy2x,"l",col=2,lwd=3,xaxs="i",ylim=c(0,y.max),xlab="f [Hz]",ylab="measure of linear causality",main=expression(y %->% x))
abline(h=0,col=8,lty=2)
plot(f.list,fx2y,"l",col=2,lwd=3,xaxs="i",ylim=c(0,y.max),xlab="f [Hz]",ylab="measure of linear causality",main=expression(x %->% y))
abline(h=0,col=8,lty=2)
##################
# 比較用の２乗コヒーレンス
plot(f.list,abs(Sxy)^2/(Re(Sx)*Re(Sy)),"l",col=2,lwd=3,xaxs="i",ylim=c(0,1),xlab="f [Hz]",ylab="coherence",main=expression(x %<->% y))
abline(h=c(0,1),col=8,lty=2)

まだ，十分に検証できていませんが (間違っているかも)，自分のメモとして高次の自己回帰過程の場合のスクリプトも載せておきます．

##################
# ２変量ARの係数
a.xx <- c(1.6,-0.9)
a.xy <- c(0)
a.yx <- c(0,0,0.1)
a.yy <- c(0.6,-0.2,0.1)
########################
n.xx <- length(a.xx)
n.xy <- length(a.xy)
n.yx <- length(a.yx)
n.yy <- length(a.yy)
#####
# 白色ノイズeps.1, eps.2の分散
sig2.1 <- 1
sig2.2 <- 1
##################
# グレンジャー因果 (Granger causality)の周波数領域表現の計算
# 単位行列
I <- diag(2)
###
# 周波数の点数(描画のためのfを何点とるか)
n <- 500
f.list <- seq(0,0.5,length.out=n)
###
Sx <- c()
Sy <- c()
Sx0 <- c()
Sy0 <- c()
Sxy <- c()
for(j in 1:n){
 f <- f.list[j]
 z <- exp(2*pi*f*1i)
 a.11 <- a.xx %*% z^(1:n.xx)
 a.12 <- a.xy %*% z^(1:n.xy)
 a.21 <- a.yx %*% z^(1:n.yx)
 a.22 <- a.yy %*% z^(1:n.yy)
 A <- matrix(c(a.11,a.12,a.21,a.22),nrow=2,ncol=2,byrow=TRUE)
 # 逆行列
 B <- solve(I-A)
 # スペクトルの計算
 Sx[j] <- B[1,1]*Conj(B[1,1])*sig2.1+B[1,2]*Conj(B[1,2])*sig2.2
 Sx0[j] <- B[1,1]*Conj(B[1,1])*sig2.1
 Sy[j] <- B[2,1]*Conj(B[2,1])*sig2.1+B[2,2]*Conj(B[2,2])*sig2.2
 Sy0[j] <- B[2,2]*Conj(B[2,2])*sig2.2
 Sxy[j] <- B[1,1]*Conj(B[2,1])*sig2.1+B[1,2]*Conj(B[2,2])*sig2.2
}
################
# 線形因果の指標
fy2x <- log(abs(Sx)/abs(Sx0))
fx2y <- log(abs(Sy)/abs(Sy0))
# 描画
y.max <- max(1,fy2x,fx2y)
par(mfrow=c(2,3),cex.axis=1.2,cex.lab=1.5,las=1)
plot(f.list,Re(Sx),"l",col=2,lwd=3,xaxs="i",xlab="f [Hz]",ylab="Sx(f)",main="Power spectrum of x")
plot(f.list,Re(Sy),"l",col=2,lwd=3,xaxs="i",xlab="f [Hz]",ylab="Sy(f)",main="Power spectrum of y")
plot(f.list,abs(Sxy),"l",col=2,lwd=3,xaxs="i",xlab="f [Hz]",ylab="|Sxy(f)|",main="|Cross spectrum|")
###
plot(f.list,fy2x,"l",col=2,lwd=3,xaxs="i",ylim=c(0,y.max),xlab="f [Hz]",ylab="measure of linear causality",main=expression(y %->% x))
abline(h=0,col=8,lty=2)
plot(f.list,fx2y,"l",col=2,lwd=3,xaxs="i",ylim=c(0,y.max),xlab="f [Hz]",ylab="measure of linear causality",main=expression(x %->% y))
abline(h=0,col=8,lty=2)
##################
# 比較用の２乗コヒーレンス
plot(f.list,abs(Sxy)^2/(Re(Sx)*Re(Sy)),"l",col=2,lwd=3,xaxs="i",ylim=c(0,1),xlab="f [Hz]",ylab="coherence",main=expression(x %<->% y))
abline(h=c(0,1),col=8,lty=2)