【確率過程・時系列解析】自己回帰過程の特性方程式の根 (数値解)

前回は，自己回帰過程の特性方程式の話をしました．
chaos-kiyono.hatenablog.com
つまり， $M$ 次自己回帰過程

　 $\begin{eqnarray}　y[n] &=& a_1 y[n-1] + a_2 y[n-2] + \cdots + a_m y[n-m] + w[n] \\ &=& \sum_{m=1}^M a_m y[n-m] + w[n] \end{eqnarray}$

について，

特性方程式
　 $\displaystyle 1- a_1 z- a_2 z^{2} - \cdots - a_M z^{M} = 0$
の解 (根)の絶対値がすべて1より大きければ，自己回帰過程は発散しない．

ということを説明しました．今回は，Rをつかって実際に特性方程式の根を求めてみます．

　使うスクリプトは，以下です．

# 自己回帰過程の係数を与える
# a <- c(a1,a2,a3, ... )
a <- c(0.9*sqrt(3),-0.81)
##############################
# 根を計算
roots <- polyroot(c(-1,a))
##############################
# 根を複素数平面にプロット
r <- abs(roots)
r.max <- max(1,r)
x <- Re(roots)
y <- Im(roots)
par(pty="s")
plot(x,y,pch=4,col=2,xlim=c(-r.max,r.max),ylim=c(-r.max,r.max),xlab="Re(z)",ylab="Im(z)")
# 単位円の描画
curve( sqrt(1-x^2),xlim=c(-1,1),col=8,add=TRUE)
curve(-sqrt(1-x^2),xlim=c(-1,1),col=8,add=TRUE)
# 根
roots

　このスクリプトでは，例として，２次自己回帰過程で， $a_1=0.9 \sqrt{3}$ ， $a_2=-0.81$ の場合を計算しています．実行すると以下の図が描かれます．円は半径１の単位円 $|z|=1$ です．赤いxが根ですので，xが円の外側にあれば，この過程は発散しません．

　以前の記事　
chaos-kiyono.hatenablog.com
で例として使った，a <- c(2.373596,-2.528090,1.404494)で試してみると，

となります．すいません，発散してしまいます．安定になる条件を，間違えました．

# ３次自己回帰過程の係数
a <- c(2.373596,-2.528090,1.404494)
# 繰り返し回数
N <- 30
########################################
y <- numeric(N)
w <- rnorm(N)
y[1:3] <- w[1:3]
for(n in 4:N){
 y[n] <- a[1]*y[n-1]+a[2]*y[n-2]+a[3]*y[n-3]+w[n]
}
plot(1:n,y,"l",las=1)

数値実験してみると，下の図のようになります．この例では，マイナス側に発散しています (常にプラスかマイナスに発散します)．