【Rで時系列解析】AICで多変量自己回帰過程の次数を決定 - ケィオスの時系列解析メモランダム

今回は，2変量の自己回帰過程

　 $\begin{eqnarray} \left\{\begin{array}{rcl} y [ t ]\!\!\!\! &=& \!\!\!\! a_{11} y [t-1] + a_{12} y [t-1] + \cdots + a_{1m} y [t-m] \\ && + b_{11} z [t-1] + b_{12} z [t-2] + \cdots + b_{1m} z [t-m] + \varepsilon_1 [t] \\ z [ t ] \!\!\!\! &=& \!\!\!\! a_{21} y [t-1] + a_{22} y [t-1] + \cdots + a_{2m} y [t-m] \\ && + b_{21} z [t-1] + b_{22} z [t-2] + \cdots + b_{2m} z [t-m] + \varepsilon_2 [t] \\ \end{array} \right. \end{eqnarray}$

について最小2乗法でパラメタを推定し，AICで次数 $m$ を決定します．お気楽な"vars"パッケージのVARなどは使いません．計算過程や式を確認できるように，Rスクリプトを書きました．Rスクリプトは一番下に掲載してあります．AICについては，前回の記事を参考にしてください．
chaos-kiyono.hatenablog.com

Rスクリプトの使い方
AICの計算式
数値実験の例
まとめ
varsパッケージのVARselectの不思議
Rスクリプト
- 時系列の生成
- 2変量自己回帰過程のパラメタ推定とAICの計算

Rスクリプトの使い方

　まずは，サンプル時系列を数値的に生成してください．下にある「時系列の生成」で示したRスクリプトを実行すれば，yとzに時系列が入ります．

　時系列が，yとzに収まっている状態で，「2変量自己回帰過程のパラメタ推定とAICの計算」のRスクリプトを実行すれば，AICが最小になるパラメタが推定されます．

AICの計算式

　AICの表式は，Helmut Luetkepohl著の"New Introduction to Multiple Time Series Analysis" pp. 147の式(4.3.2)を参考にしました．この本は，先週買いました．
　
　私が使った， $m$ 次の2変量自己回帰過程のAICは，

　 $\displaystyle {\rm AIC}(m) = \log \left( {\rm det }\, \tilde{\Sigma}_{\varepsilon} (m) \right) + \frac{2}{N-m}$ (モデルのパラメタ数)

ここで，
　 $\tilde{\Sigma}_{\varepsilon} (m) = \begin{bmatrix} \displaystyle \frac{1}{N-m} \sum_{t=m+1}^N \hat{\varepsilon_1} [t]^2 & \displaystyle\frac{1}{N-m} \sum_{t=m+1}^N \hat{\varepsilon_1} [t] \hat{\varepsilon_2} [t] \\ \displaystyle \frac{1}{N-m} \sum_{t=m+1}^N \hat{\varepsilon_2} [t] \hat{\varepsilon_1} [t] & \displaystyle\frac{1}{N-m} \sum_{t=m+1}^N \hat{\varepsilon_2} [t]^2 \end{bmatrix}$

です ( $N-m$ で割るのか， $N$ で割るのか，といったことに自信がありません)．"vars"パッケージの結果と一致するので， $N-m$ で割るでいいと判断しました．もしかすると第2項を， $N$ で割る流儀もあるかもしれません．

数値実験の例

　下に掲載した「2変量自己回帰過程のパラメタ推定とAICの計算」のRスクリプトを実行すると，以下のような結果がRコンソールに出力されます．ここでは、時系列の生成に，

　 $\begin{eqnarray} \left\{\begin{array}{l} y [ t ]= 0.7 y [t-1] + 0.3 z [t-2] + \varepsilon_1 [t] \\ z [ t ]= 0 y [t-1] - 0.4 z [t-1] + \varepsilon_2 [t] \end{array} \right. \end{eqnarray}$

を使いました．

> # 計算したAIC
> data.frame(AIC)
         AIC
1 0.06288797
2 0.06639792
3 0.06505095
4 0.06923661
5 0.07457428
> # AICが最小になる次数
> print(i.min <- which.min(AIC))
[1] 1
> ###
> # y[t] = a.11*y[t-1] + ... + b.11*z[t-1] + ... + const
> a1_.est[[i.min]]
[1]  0.71541602  0.30857341 -0.01492857
> # z[t] = a.21*y[t-1] + ... + b.21*z[t-1] + ... + const
> a2_.est[[i.min]]
[1] -0.03231752 -0.44302379  0.06051948
> ###
> diag(sig2.est[[i.min]])
[1] 1.0884622 0.9668476

　AICは， $m=1$ のときに最小になっています．
　
　 $y [ t ]= 0.7 y [t-1] + 0.3 z [t-2] + \varepsilon_1 [t]$ のパラメタについては，

[1]  0.71541602  0.30857341 -0.01492857

と推定されています．正解が0.7のパラメタは，0.71541602，正解が0.3のパラメタは，0.30857341と推定されています．3番目の数値は，モデルの式にはない，定数項を推定しています．

　また， $z [ t ]= 0 y [t-1] - 0.4 z [t-1] + \varepsilon_2 [t]$ のパラメタについては，

[1] -0.03231752 -0.44302379  0.06051948

と推定されています．正解が0のパラメタは，-0.03231752，正解が-0.4のパラメタは，-0.44302379と推定されています．ここでも，3番目の数値は，モデルの式にはない，定数項を推定しています．

　より高次のモデルの場合も正しく推定できるか，数値実験してみてください．

まとめ

　これは、私の勉強メモですので，間違いがあるかもしれません．間違いを見つけたら正しい答えを教えてもらえると大変ありがたいです．タダで教えてもらおうなどと，甘えていますが，なにとぞよろしくお願いします．

　私のRスクリプトでAICが正しく計算できているのか心配だったので，varsパッケージのVARselectと比べてみました．VARselectで計算されたAICと，私のRスクリプト計算したAICでは，

VARselect(... , lag.max=m.max, type="const")

のように，m.maxでの値は一致するのに，m.max未満では一致しませんでした．奇妙なことに，VARselectでは，lag.maxの設定値m.maxを増やすと，m.max未満のAICの値がm.maxに依存して変化しました．ちょっとおかしいと感じましたが，私が理解できていないAICの調整法が存在するのかもしれません．秘密を知っている方は教えてください．

varsパッケージのVARselectの不思議

　VARselectコマンドを使って，同じデータを評価しても，lag.maxの設定を変えるとAICの値が変ってしまいます．arコマンドで計算されるAICは，order.maxの設定を変えても，値は同じままです．理由が分かりませんが，理由を教えてもらいたいので例を示しておきます．

lag.max=1のとき
> VARselect(data.frame(y,z), lag.max=1, type="const")
（省略）
$criteria
1
AIC(n) -0.05334486

lag.max=2のとき　AIC(1)の値が大きくなります．
> VARselect(data.frame(y,z), lag.max=2, type="const")
（省略）
$criteria
1 2
AIC(n) -0.05152919 -0.044359411

lag.max=3のとき　AIC(1)とAIC(2)の値が大きくなります．
> VARselect(data.frame(y,z), lag.max=3, type="const")
（省略）
$criteria
1 2 3
AIC(n) -0.05095212 -0.043739762 -0.04107069

私のRスクリプトでAICを計算した結果
AIC
1 -0.05334486　(lag.max=1のときの値と一致)
2 -0.04435941　(lag.max=2のときの値と一致)
3 -0.04107069　(lag.max=3のときの値と一致)

　どこかの論文にAICの調整法が提案されているのか，VARselectにバグがあるのか，分かる人は理由を教えてください．

Rスクリプト

時系列の生成

　数値的にサンプル時系列を生成するスクリプトです．長さNの時系列が生成されます．係数のベクトルは，0が入る場合も，すべて長さをそろえてください．

1次の例
　下の例では，

$\begin{eqnarray} \left\{\begin{array}{l} y [ t ]= 0.7 y [t-1] + 0.3 z [t-2] + \varepsilon_1 [t] \\ z [ t ]= 0 y [t-1] - 0.4 z [t-1] + \varepsilon_2 [t] \end{array} \right. \end{eqnarray}$

の時系列を生成しています． $\varepsilon_1$ ， $\varepsilon_2$ の分散sig2.1，sig2.2はどちらも1にしてあります．

# データ長
N <- 1000
##################
# AR係数
# y[t] = a.11*y[t-1] + a.12*y[t-2] + ... + b.11*z[t-1] + b.12*z[t-2] + ... + eps1[t] 
# z[t] = a.21*y[t-1] + a.22*y[t-2] + ... + b.21*z[t-1] + b.22*z[t-2] + ... + eps2[t] 
a.1_ <- c(　0.7) # y[t-1] の係数
b.1_ <- c(　0.3) # z[t-1] の係数
a.2_ <- c(　  0) # y[t-1] の係数
b.2_ <- c(-0.4) # z[t-1] の係数
# eps1, eps2の分散
sig2.1 <- 1
sig2.2 <- 1
##################
n.a1 <- length(a.1_)
n.b1 <- length(b.1_)
n.a2 <- length(a.2_)
n.b2 <- length(b.2_)
##################
# サンプル時系列の生成
y <- c()
z <- c()
##################
# 過渡状態として余分にまわす回数
N.add <- 10
N.add <- max(N.add,n.a1*4)# 念のため
##################
w1 <- rnorm(N+N.add,0,sqrt(sig2.1))
w2 <- rnorm(N+N.add,0,sqrt(sig2.2))
y[1:n.a1] <- w1[1:n.a1]
z[1:n.a1] <- w2[1:n.a1]
for(i in (n.a1+1):(N+N.add)){
 y[i] <- a.1_ %*% y[(i-1):(i-n.a1)]+ b.1_ %*% z[(i-1):(i-n.b1)] + w1[i]
 z[i] <- a.2_ %*% y[(i-1):(i-n.a2)]+ b.2_ %*% z[(i-1):(i-n.b2)] + w2[i] 
}
y <- y[-(1:N.add)]
z <- z[-(1:N.add)]

2次の例
　下の例では，

$\begin{eqnarray} \left\{\begin{array}{l} y [ t ]= 0.7 y [t-1] + 0.1 y [t-2] + 0 z [t-1] + 0.3 z [t-2] + \varepsilon_1 [t] \\ z [ t ]= 0 y [t-1] + 0 y [t-2] + 0.8 z [t-1] - 0.4 z [t-2] + \varepsilon_2 [t] \end{array} \right. \end{eqnarray}$

の時系列を生成しています． $\varepsilon_1$ ， $\varepsilon_2$ の分散sig2.1，sig2.2はどちらも1にしてあります．

# データ長
N <- 1000
##################
# AR係数
# y[t] = a.11*y[t-1] + a.12*y[t-2] + ... + b.11*z[t-1] + b.12*z[t-2] + ... + eps1[t] 
# z[t] = a.21*y[t-1] + a.22*y[t-2] + ... + b.21*z[t-1] + b.22*z[t-2] + ... + eps2[t] 
a.1_ <- c(0.7, 0.1) # y[t-1] y[t-2] の係数
b.1_ <- c(  0, 0.3) # z[t-1] z[t-2] の係数
a.2_ <- c(  0,   0) # y[t-1] y[t-2] の係数
b.2_ <- c(0.8,-0.4) # z[t-1] z[t-2] の係数
# eps1, eps2の分散
sig2.1 <- 1
sig2.2 <- 1
##################
n.a1 <- length(a.1_)
n.b1 <- length(b.1_)
n.a2 <- length(a.2_)
n.b2 <- length(b.2_)
##################
# サンプル時系列の生成
y <- c()
z <- c()
##################
# 過渡状態として余分にまわす回数
N.add <- 10
N.add <- max(N.add,n.a1*4)# 念のため
##################
w1 <- rnorm(N+N.add,0,sqrt(sig2.1))
w2 <- rnorm(N+N.add,0,sqrt(sig2.2))
y[1:n.a1] <- w1[1:n.a1]
z[1:n.a1] <- w2[1:n.a1]
for(i in (n.a1+1):(N+N.add)){
 y[i] <- a.1_ %*% y[(i-1):(i-n.a1)]+ b.1_ %*% z[(i-1):(i-n.b1)] + w1[i]
 z[i] <- a.2_ %*% y[(i-1):(i-n.a2)]+ b.2_ %*% z[(i-1):(i-n.b2)] + w2[i] 
}
y <- y[-(1:N.add)]
z <- z[-(1:N.add)]

2変量自己回帰過程のパラメタ推定とAICの計算

　2変量時系列に自己回帰過程を当てはめるRスクリプトです．lag.maxで，推定する最大次数を設定してください．下のスクリプトでは，lag.max <- 5です．

# 推定する最大次数
lag.max <- 5
###############
AIC <- c()
a1_.est <- c()
a2_.est <- c()
sig2.est <- c()
###############
for(ord.est in 1:lag.max){
　# 最小2乗フィット
　Y <- matrix(y[(1+ord.est):N],ncol=1)
　Z <- matrix(z[(1+ord.est):N],ncol=1)
　tmp <- c()
 for(k in 1:ord.est) tmp <- c(tmp,y[(ord.est+1-k):(N-k)])
 for(k in 1:ord.est) tmp <- c(tmp,z[(ord.est+1-k):(N-k)])
 tmp <- c(tmp,rep(1,N-ord.est))
　X <- matrix(tmp,ncol=2*ord.est+1)
　#####################
　# 係数の推定
　a1_ <- solve(t(X) %*% X) %*% (t(X) %*% Y)
　a2_ <- solve(t(X) %*% X) %*% (t(X) %*% Z)

　a1_.est <- c(a1_.est,list(as.numeric(a1_)))
　a2_.est <- c(a2_.est,list(as.numeric(a2_)))
　# 残差の２乗平均
　eps1 <- Y - X %*% a1_
　eps2 <- Z - X %*% a2_
　sig2.11 <- mean(eps1^2)
　sig2.12 <- sig2.21 <- mean(eps1*eps2)
　sig2.22 <- mean(eps2^2)
 tmp <- matrix(c(sig2.11,sig2.21,sig2.12,sig2.22),nrow=2)
　det.sig <- det(tmp)
 sig2.est <- c(sig2.est,list(tmp))
　##############################################
　AIC <- c(AIC,log(det.sig)+2*(ord.est*2^2+2)/(N-ord.est))
}
###############################################
# 計算したAIC
data.frame(AIC)
# AICが最小になる次数
print(i.min <- which.min(AIC))
###
# y[t] = a.11*y[t-1] + ... + b.11*z[t-1] + ... + const
a1_.est[[i.min]]
# z[t] = a.21*y[t-1] + ... + b.21*z[t-1] + ... + const
a2_.est[[i.min]]
###
diag(sig2.est[[i.min]])