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

今回は，時系列に対して自己回帰過程をフィットし，最適な次数を決定します．自己回帰過程の次数の推定には，赤池情報量規準 (AIC: Akaike Information Criterion)を使います．AICの導出をここで説明することは，私には難しいのでやりません．でも，計算過程については見て分かるように，Rスクリプトを書きました．つまり，Rのarなどの便利コマンドは使いません (ちょっと説明しますが)．

　ここでは，推定されたパラメタが正しいか確認したいので，時系列は自分で生成します．使ったRスクリプトは，最後に掲載しました．

赤池情報量規準 (AIC: Akaike Information Criterion)とは
数値実験で自己回帰過程の次数を決定
- 時系列の生成
- 自己回帰過程のパラメタの推定とAICの計算
予告と目標
使ったRスクリプト
- 時系列の生成に使ったRスクリプト
- 自己回帰過程のフィットとAICの計算に使ったRスクリプト

赤池情報量規準 (AIC: Akaike Information Criterion)とは

　時系列に対して何かモデルを当てはめる場合，複数のモデル候補から，どれが一番ふさわしいかを評価する必要があります．その評価基準の一つがAICです．

　既に観測済みの時系列にモデルをフィットするとき，その時系列を良く再現する (予測誤差が小さい)モデルが良いもののような気がします．しかし，未知の時系列に対しても，推定したモデルの予測が良い (予測誤差が小さい)かどうかは分かりません．既に観測済みのデータの再現にこだわっていると，当てはめすぎ (overfitting)という恐ろしい罠にはまっているかもしれません．AICでは，モデルのパラメタ数の増加は当てはめすぎにつながり，むしろ，未知のデータに対する予測を悪くする，という影響も考慮して，モデル予測の良さを評価しています．

　ということで，今回は $m$ 次の自己回帰過程

　 $x[t]=a_0+a_1 x[t-1] + \cdots + a_m x[t-m] + \varepsilon[t]$

をサンプル時系列 $\{x[t]\}$ に当てはめ，AICが最小になるという基準で，自己回帰過程の次数を決定します．推定するパラメタは， $\{a_0, a_1, \cdots, a_m\}$ と， $\varepsilon[t]$ の分散 $\sigma^2$ です ( $\varepsilon[t]$ の平均は0です)．そして， $m$ も推定します．

　このとき，AICは，

　 ${\rm AIC} =$ (用いるデータ数) $\cdot \log \big\{ 2 \pi\cdot$ (残差の分散) $+1\big\}+2$ (パラメタ数)

で与えられます．田村先生が書かれた記事［田村義保. (2011). 統計的データ解析入門第 2 回: 情報量規準と時系列解析. 知能と情報, 23(1), 107-113］を参考に，勉強させていただきました．田村先生には，昔，統計数理研究所の共同研究集会で10年ほどお世話になりました．ありがとうございました．

数値実験で自己回帰過程の次数を決定

時系列の生成

　下に掲載したRスクリプトにある

# 自己回帰 (AR)過程のパラメタ
a <- c(0.8,-0.2)
sig2 <- 1
#########################
# 時系列の長さ
N <- 1000

の部分で，自己回帰過程のパラメタと，生成する時系列の長さを指定してください．上の例は， $a_1=0.8$ ， $a_2=-0.2$ ， $\sigma^2=1$ ということです．自己回帰過程の次数 $m$ は，a <- c(0.8,-0.2)の長さで自動的に設定されます．

　Rスクリプトを実行すれば，x.simに時系列が入っています．平均は0にしてあります．

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

　下の掲載したRスクリプトを実行すれば，自己回帰過程の次数とパラメタが推定されます．AICのグラフが描かれて，推定結果はRコンソールに出力されます．AICは，最小のAICを引いた結果がプロットされます．つまり，0になるのが最小の次数です．

　Rスクリプトは，そのままで動くと思いますが，必要に応じて評価する自己回帰過程の最大次数を変更してください．これは，ord.max <- 10の部分の値を変更するだけです．

　うまく推定できた例を掲載すると，

> # AICが最小になる次数
> print(ord.est <- which.min(aic)-1)
[1] 2
> # 自己回帰過程の係数 a0, a1, ...
> a.est[[ord.est+1]]
[1]  0.00127796  0.78739197 -0.19826019
> # 白色ノイズの分散
> sig2.est[ord.est+1]
[1] 0.9860711

のような結果がでます．次数は2で正しく評価できています．注意点は，本来，0のはずの $a_0$ が，ここでは推定されているということです．つまり，パラメタの推定結果は，

　 $a_0=0.00127796, a_1=0.78739197, a_2=-0.19826019$

となってます．これは，伝統に従ったもので，田村先生の記事でも，Rのarコマンドでも， $a_0$ を推定しています．

　同じことを，Rのarコマンドで実行すれば，

> ar(x.sim, method="ols")

Call:
ar(x = x.sim, method = "ols")

Coefficients:
      1        2  
 0.7874  -0.1983  

Intercept: 0.001278 (0.03143) 

Order selected 2  sigma^2 estimated as  0.9861

となります．Intercept: 0.001278が， $a_0$ の推定結果です．

　繰り返し数値実験を試みてください．必ずしも正しい次数にならないことがわかります．しかし，次数が間違っていても，元々0のモデルパラメタの推定結果は，0に近くなっていると思います．

予告と目標

　多変量の自己回帰過程でもAICを計算して次数を決めてみたいと思います．目標は，長時間相互相関過程のグレンジャー因果を評価する方法を開発することです．

使ったRスクリプト

時系列の生成に使ったRスクリプト

# 自己回帰 (AR)過程のパラメタ
a <- c(0.8,-0.2)
sig2 <- 1
#########################
# 時系列の長さ
N <- 1000
#########################
ord.ar <- length(a)
#########################
# 白色ノイズの生成
WN <- rnorm(N)
#　白色ノイズのフーリエ変換
fft.WN <- fft(WN)
#########################
# 周波数を定義
freq <- (0:(N-1))/N
#########################
#　ARモデルのパワースペクトル
Sf <- c()
for(k in 1:N){
 Sf[k] <- 1/abs((c(1,-a) %*% exp(2*pi*freq[k]*(0:ord.ar)*1i)))^2
}
fft.sim <-  sqrt(Sf)*fft.WN
#########################
# フーリエ逆変換で時系列を生成
x.sim <- Re(fft(fft.sim,inverse=TRUE))/N
# 平均0に変換
x.sim <- x.sim - mean(x.sim)

自己回帰過程のフィットとAICの計算に使ったRスクリプト

# パラメタの推定（ord.maxで推定する最大次数を設定）
ord.max <- 10
a.est <- 0
sig2.est <- mean(x.sim^2)
aic <- N*(log(2*pi*sig2.est)+1)+2
for(ord.est in 1:ord.max){
 Y <- matrix(x.sim[(1+ord.est):N],ncol=1)　# ncolで列数のみ指定
 tmp <- rep(1,N-ord.est)
 for(k in 1:ord.est){
  tmp <- c(tmp,x.sim[(ord.est+1-k):(N-k)])
 }
 X <- matrix(tmp,ncol=ord.est+1)　# ncolで列数のみ指定
 #####################
 # 係数の推定
 a.tmp <- solve(t(X) %*% X) %*% (t(X) %*% Y)
 a.est <- c(a.est,list(as.numeric(a.tmp)))
 # 残差の２乗平均
 sig2.tmp <- sum((Y - X %*% a.tmp)^2)/(N-ord.est)
 sig2.est <- c(sig2.est,sig2.tmp)
 aic <- c(aic,N*(log(2*pi*sig2.tmp)+1)+2*(ord.est+1))
 ######################
}
plot(0:ord.max,aic-min(aic),"o",pch=16,col=2,lwd=2,cex=2,xlab="AR order")
abline(h=0,lty=2,col=8)
####################################
# AICが最小になる次数
print(ord.est <- which.min(aic)-1)
# 自己回帰過程の係数 a0, a1, ...
a.est[[ord.est+1]]
# 白色ノイズの分散
sig2.est[ord.est+1]