R Time Series Analysis 時系列解析(7) (線形定常過程2) 2次ARモデル

こっちでAR1次モデルを扱ったが、2次以降も項が増えるだけで、理屈は似ている。

待ち時間にやってみた。今度はACFの結果も載せる。

AR2 <- function(n,b1,b2,a=0,sd=1,y0=0) {
  c <- rnorm(n,sd=sd)
  y <- rep(0,n)
  y[1] <- y0
  
  for(j in 1:(n-2)) {
    y[j+2] <- a + b1*y[j+1] + b2*y[j] + c[j+2] 
  }
  
  m <- round(mean(y),digits=2)
  e_y <- round(m/(1-b1-b2),digits=2)
  
  title <- paste("AR Model y_t=a + b1*y_t-1 + b2*y_t-2 + c_t, 
                 a=",a,",b1=",b1,",b2=",b2,",sd=",sd)
  yl <- c(min(y),max(y))
  plot(y,type="l",main=title)
  abline(h=e_y,col=2)
  
  y_ts <- ts(y,start=1,freq=1)
  acf(y_ts,main="Autocorrelation coefficient")
}

以下、実行コード。

par(mfrow=c(7,2))
AR2(n=200,b1=0.3,b2=0.3)    #1番目
AR2(n=200,b1=0.3,b2=0.6)　  #2番目
AR2(n=200,b1=0.3,b2=0.9)    #3番目
AR2(n=200,b1=0.6,b2=0.3)　  #4番目
AR2(n=200,b1=0.9,b2=0.3)　  #5番目
AR2(n=200,b1=-0.3,b2=-0.3)  #6番目
AR2(n=200,b1=-0.3,b2=0.3)   #7番目

b1, b2の値を様々に変えることで、様々な系列を表現出来る。1, 2, 4番目は株価、6,7番目は音声データのような形になっている気がする。こうやって波形の特徴を係数の条件に表せれば便利ですね。ちなみに、2次ARモデルの収束条件は、

　　　　　　　　b2 - b1 < 1,  b1 + b2 < 1, -1 < b1 < 1

となる(漸化式を解き、等比級数部分の公比部分から導ける)。これはb1b2平面上で3角形となる。これを拡張すれば、n次ARモデルの収束条件はn次元空間のn+1個の頂点を持つ立体となるのか？？暇があったら考えることにする。

以下グラフ。赤い線は平均の理論値。青い線は標本自己相関係数が0と有意に異なることを示す95%信頼区間。係数が大きい方が遠く先まで影響を及ぼしていることがわかります。