1 AR(2)モデル.

さて少し欲張って、AR(2) モデルについてシミュレーションしましょう。
高校で学ぶ隣接３項間漸化式の特性方程式と関連が深く興味深いところです。
将来の学習のために覗いてみましょう。
尚、数理の難しいところは公式の確認をして進んでください。

AR(2)モデルは、２過去までのデータで表されるモデルで、隣接３項間の漸化式に誤差項を加えた式で表されます。

　　\(y_{t}=m+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\varepsilon _{t}\)

コードは、ゼロ数列（全ての項が０）を作っておき、漸化式に従って数値を入れ替えます。

例　\(y_1=0,~y_2=1,~~~y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}+\varepsilon _{t}\)

  • AR(1)では、定常であるためには係数の絶対値が１より小さいことが必要でした。
  • ここでは、\(\phi_1=1.2>1\) ですが、定常性はどうなるでしょうか？
y <- rep(0, 200)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 2 + 1.2 * y[i - 1] -  0.8 * y[i - 2]+rnorm(1,0,1) 
}
plot(y,type="l",xlab="time")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

\(\phi_1=1.2>1\) の影響は、\(\phi_2=-0.8<0\) で相殺され、定常となるようです。
　

＜平均と分散＞

\(y_{t}=m+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\varepsilon _{t}\)　が定常性を満たせば、

  • \(E(y_{t})=m+\phi_1 E(y_{t-1})+\phi_1 E(y_{t-2})+E(\varepsilon _{t})\) から、

　　　　　\(\mu=m+\phi_1 \mu+\phi_2 \mu~~~\therefore~\mu=\dfrac{m}{1-\phi_1-\phi_2}\)

一方、分散については、\(y_{t-1},~y_{t-2}\) が独立でないので簡単にはいきません。
ユール・ウォーカー方程式.
　　　\(\begin{split} \gamma _{0}&=\phi_1 \gamma_{1}+\phi_2~\gamma_{2}+\sigma^2\\ \gamma_1 &=\phi_{1}\gamma_0 + \phi_{2}\gamma_1\\ \gamma_2 &=\phi_{1}\gamma_1 + \phi_{2}\gamma_0 \end{split}\)

を \(\gamma_0,\gamma_1,\gamma_2\) の連立方程式として解くと、\(\gamma_0\) が分散です。

自己相関係数 \(\rho_1,\rho_2\) を用いる場合は、ユール・ウォーカーの式を \(\gamma_0\) で割った.

　　　\(1=\phi_{1}\rho_1 + \phi_{2}\rho_2+\sigma^2/\gamma_0\)

　　\(\begin{pmatrix} \rho_{1} \\ \rho_{2} \end{pmatrix}= \begin{pmatrix} 1& \rho_{1} \\ \rho_{1} &1 \end{pmatrix} \begin{pmatrix} \phi _{1} \\ \phi _{2} \end{pmatrix}\) 　 から 　

　　　\(\gamma_0=\dfrac{\sigma^2}{1-\phi_{1}\rho_1 - \phi_{2}\rho_2}\)　　によって分散を求めることができます。
まとめると、

  • \(V(y_{t})=\dfrac{\sigma^2}{1-\phi_1 \rho_1-\phi_2 \rho_2}=\dfrac{1-\phi_{2}}{1+\phi _{2}}\cdot \dfrac{\sigma^2}{(1-\phi_{2})^2-\phi _{1}^2}\)
  • \(\rho_1=\dfrac{\phi_1}{1-\phi_2}\)

標本値と理論値を比べてみましょう。

mv <- c(mean(y),sd(y))
mv #平均と標準偏差の標本値
## [1] 3.383714 2.070898
MV <- c(2/(1-1.2+0.8),sqrt((1+0.8)/(1-0.8)*1^2/((1+0.8)^2-1.2^2)))  
MV #理論値
## [1] 3.333333 2.236068

標本値は理論値と近い値になっています。
グラフを見ると、平均は変動の中心、標準偏差は変動幅としての意味を持っていることが分かります。
青の点線の間が 95% 信頼区間であり、変動の 95% がこの範囲内に入ります。

1.0.1 自己相関・偏自己相関

\(Y_{t}=\phi_{0}+\phi_{1} Y_{t-1}+ \phi_{2} Y_{t-2}+\varepsilon _{t}\) に対し、
　　　\(\begin{split} \gamma _{h}&=cov\left( Y_{t},Y_{t-h}\right)\\ &=cov\left(\phi_{1} Y_{t-1}+\phi_{2} Y_{t-2}+\varepsilon _{t},Y_{t-h}\right)\\ &=\phi_{1} cov\left( Y_{t-1},Y_{t-h}\right)+\phi_{2} cov\left( Y_{t-2},Y_{t-h}\right)+cov\left( \epsilon_{t},Y_{t-h}\right)\\ &=\phi_{1} \gamma_{h-1}+\phi_{2}\gamma_{h-2} \end{split}\)
　自己共分散は、この隣接３項間の漸化式によって順次計算できます。
　また、この式において h=1,2 としたものがユール・ウォーカー方程式です。
　（\(\gamma_{h}=\gamma_{-h}\)　に注意）.

両辺を \(\gamma_0\)で割れば、
　　　\(\begin{split} \rho _{h} &=\phi_{1} \rho_{h-1}+\phi_{2}\rho_{h-2} \end{split}\)　　（\(\rho_{h}=\rho_{-h}\)、\(\rho_{0}=1\) に注意）

漸化式の考えにより逐次 \(\rho_h\) を求めることができます。
（この式は、時系列モデルの式から定数項と誤差項を除いたものと同等です）

\(y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}+\varepsilon _{t}\)　のときの理論値は、

  • \(\rho_1=\dfrac{\phi_1}{1-\phi_2}=\dfrac{1.2}{1+0.8}=0.667\)
  • \(\rho_2=\phi_{1}\dfrac{\phi_1}{1-\phi_2} + \phi_{2}=\dfrac{\phi_2+\phi_1^2-\phi_2^2}{1-\phi_2}=\dfrac{-0.8+1.2^2-(-0.8)^2}{1+0.8}=0\)
acf(y,plot=F,lag.max=7)
## 
## Autocorrelations of series 'y', by lag
## 
##      0      1      2      3      4      5      6      7 
##  1.000  0.648 -0.004 -0.481 -0.550 -0.292  0.060  0.272
par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(mfrow=c(1,2))
acf(y,plot=T,main="自己相関")
acf(y,plot=T,type="p",main="偏自己相関")

自己相関係数の標本値は理論値とほぼ一致しています。
ラグ１では偏自己相関係数と自己相関係数は同じになります。
間接的影響を与える項が間にないためです。

コレログラムでは、ラグ20までしか表示がありませんが、ラグ200まで表示してみましょう。

par(family = "HiraKakuProN-W3") #日本語フォントの指定
acf(y,plot=T,main="自己相関係数",lag.max=200)

時系列が定常であれば、自己相関係数は 0 に収束していることが分かります。
（逆は成り立ちません）

時系列が定常であれば、誤差項を除いた漸化式は収束するはずです。

　　　\(y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}+\varepsilon _{t}\)

から誤差項を除いた漸化式をシミュレーションします。

＜誤差項のない時系列＞

　　　\(y_1=0,~y_2=1,~~~y_{t}=2+1.2 y_{t-1}-0.8y_{t-2}\).

y <- rep(0, 200)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 2 + 1.2 * y[i - 1] -  0.8 * y[i - 2]
}

plot(y,type="l",xlab="time")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

mv <- c(mean(y),sd(y))
mv #平均と標準偏差 
## [1] 3.3194444 0.4033081

誤差項のない漸化式で与えられる数列は 0 に収束しています。
変動の様子はコレログラムと似ています。
AR(2) でも、誤差項のない漸化式が収束するとき、時系列は定常になります。

2 AR(2) の定常性.

ARモデルが定常となる条件を整理しておきます。

AR(1)モデル　\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\)　では、

  • 漸化式　\(y_{t}=c+\phi y_{t-1}\)　が収束するとき定常になりました.
    公比 \(|\phi|<1\) が収束の条件です。

AR(2)モデル　\(y_{t}=m+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\varepsilon _{t}\)　では、

  • 漸化式　\(y_{t}=m+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}\)　が収束するとき定常になります。
    漸化式の解法では、２次方程式 \(x^2=\phi_{1}x+\phi_{2}\) の解が公比として現れます。
    （高校数学の隣接３項間漸化式の範囲で説明できます）

従って、 \(x^2-\phi_{1}x-\phi_{2}=0\)　 の解の絶対値が１より小さいことが定常の条件となります。

２次方程式の解の理論として考えると、

　　\(1-\phi_{1}-\phi_{2}>0\)、\(1+\phi_{1}-\phi_{2}>0\)、\(\phi _{2} >-1\)

\(\left( \phi _{1},\phi _{2}\right)\) が \((0,1),(-2,-1),(2,-1)\) を頂点とする三角形の内部の領域にあるとき定常となります。
（高校数学の２次方程式の範囲で説明できます）

課題.　上記のことを説明しなさい。
（\(x^2=\phi_{1}x+\phi_{2}\) の解が実数解の場合と虚数解の場合に分けて考える必要があります）

ひょっこり統計「時系列分析の基礎」参照.

ここでは証明には立ち入らず、シミュレーションで確かめることにしましょう。

2.1 ＜発散する例＞.

　　　\(y_1=0,~y_2=1,~~~y_{t}=3 y_{t-1}+10 y_{t-2}+\varepsilon _{t}\).

パラメーター \((\phi_1,~\phi_2)=(3,~10)\) は三角形領域の外部にあり、 \(y_t\) は発散します。

y <- rep(0, 120)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 3 * y[i - 1] +  10 * y[i - 2] +rnorm(1,0,1)
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="Time series")
abline(h=mean(y),col="red")
par(family = "HiraKakuProN-W3") #日本語フォントの指定
acf(y,plot=T,main="ACF",lag.max=120)

時系列は発散しているので定常ではありません。
コレログラムは 0 に収束しているように見えます。
しかし、定常でなければ自己相関係数はラグだけでなく時間にも依存します。
自己相関係数は、時間の初期値に応じて不安定な挙動をする可能性があり、コレログラムから何かを判断する訳にはいきません。

2.2 ＜収束する例＞.

　　　\(y_1=0,~y_2=1,~~~y_{t}=0.4 y_{t-1}+0.3 y_{t-2}+\varepsilon _{t}\).

パラメーター \((\phi_1,~\phi_2)=(0.4,~0.3)\) は三角形領域の内部にあり、 \(y_t\) は収束します。

y <- rep(0, 120)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 0.4 * y[i - 1] +  0.3 * y[i - 2]+rnorm(1,0,1)
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="Time series")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
par(family = "HiraKakuProN-W3") #日本語フォントの指定
acf(y,plot=T,main="ACF",lag.max=120)

時系列は定常であり、コレログラムは減衰しながら０に収束しているようです。

誤差項を含まない漸化式のグラフと並置すると、

par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(y,type="l",xlab="time",main="誤差項あり")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

y <- rep(0, 120)
y[1] <- 0
y[2] <- 1
for (i in 3:length(y)) {
  y[i] <- 0.4 * y[i - 1] +  0.3 * y[i - 2]
}
plot(y,type="l",xlab="time",main="誤差項なし")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

誤差項を除いた漸化式は急速に収束します。
誤差項のある時系列データは、誤差項が影響し大きく変動している訳です。

参考サイト.
無粋な日々に様「時系列分析：ARモデルの性質と定常性（1/2）」.

2.3 非定常である例.

例　\(y_1=9,~y_2=5,~~~y_{t}=2+y_{t-1}+0.3y_{t-2}+\varepsilon _{t}\).

パラメーター \((\phi_1,~\phi_2)=(1,~0.3)\) は三角形領域の外部にあり、 \(y_t\) は発散します。

y <- rep(0, 120)
y[1] <- 9
y[2] <- 5
for (i in 3:length(y)) {
  y[i] <- 2 + 1 * y[i - 1] +  0.3 * y[i - 2]+rnorm(1,0,1) 
}
par(mfrow=c(1,2))
plot(y,type="l",xlab="time",main="Time series")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)
acf(y,plot=T,main="ACF",lag.max=200)

時系列は発散し、コレログラムは 0 に収束していません。

誤差項を含まない漸化式のグラフと並置すると、

par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(y,type="l",xlab="time",main="誤差項あり")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

y <- rep(0, 120)
y[1] <- 9
y[2] <- 5
for (i in 3:length(y)) {
  y[i] <- 2 + 1 * y[i - 1] +  0.3 * y[i - 2]
}
plot(y,type="l",xlab="time",main="誤差項なし")
abline(h=mean(y),col="red")
abline(h=mean(y)+1.96*sd(y),col="blue", lty=2)
abline(h=mean(y)-1.96*sd(y),col="blue", lty=2)

2.4 漸化式との関係.

高校数学Bの数列分野では隣接３項間の漸化式を学びます。
大学教養の線形代数では数列の項をずらす変換を行列で表し、行列の固有値を用いて漸化式を解きます。
斎藤正彦先生の線形代数入門に、一次変換の例として紹介されていて衝撃を受けたものです。
少し踏み込んで紹介します。

例.
\(y_1=0,~y_2=1,~~~y_{t}=3 y_{t-1}+10y_{t-2}\).

　　　\(\begin{pmatrix} y_t \\ y_{t-1} \end{pmatrix}= \begin{pmatrix} 3 & 10 \\ 1 & 0 \end{pmatrix}\begin{pmatrix} y_{t-1} \\ y_{t-2} \end{pmatrix}\)

行列をかけることで数列の項が１つずれると考えれば、
　　　\(\begin{pmatrix} y_t \\ y_{t-1} \end{pmatrix}= \begin{pmatrix} 3 & 10 \\ 1 & 0 \end{pmatrix}^{t-2} \begin{pmatrix} y_{2} \\ y_{1} \end{pmatrix}= A^{t-2}\begin{pmatrix} y_{2} \\ y_{1} \end{pmatrix}\)

行列 A の固有方程式　\(\lambda ^{2}-3 \lambda -10 =0\)　から.
固有値 \(\lambda_1=-2 ,~\lambda_2=5\)　とその固有ベクトルによって行列は対角化され、
　　　\(\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix}^{-1} A\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix} =\begin{pmatrix} \lambda_{1} & 0 \\ 0 & \lambda_{2} \end{pmatrix}\).

　　　\(\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix}^{-1} A^{t-2}\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix} =\begin{pmatrix} \lambda_{1}^{t-2} & 0 \\ 0 & \lambda_{2}^{t-2} \end{pmatrix}\).

　　　\(\begin{pmatrix} y_t \\ y_{t-1} \end{pmatrix}= \begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix} \begin{pmatrix} \lambda_{1}^{t-2} & 0 \\ 0 & \lambda_{2}^{t-2} \end{pmatrix}\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix}^{-1} \begin{pmatrix} y_{2} \\ y_{1} \end{pmatrix}\)

よって、\(y_t=A\lambda_{1}^{t-2}+B\lambda_{2}^{t-2}\) の形となり、

　　　　\(\lambda_1=-2,~\lambda_2=5\)　なので、\(y_t\) は発散します。

重要な点は、

  • 行列の固有方程式が漸化式の特性方程式と同じであること.

  • \(y_t\) は特性方程式の解の指数乗で表され、

  • \(y_t\) の収束は、特性方程式の解の絶対値が１より小さいか否かで決まる.

ということです。

＜まとめ＞.
https://mathlog.info/articles/3816 Mathlog

3 ラグオペレーター.

前節で、隣接３項間の漸化式を行列で表現しました。
そのアイディアは、漸化式を「数列の項をずらす変換」と考えることです。
時系列分析の理論では時間tをずらす変換としてラグオペレーターが用いられます。

ラグオペレーターを、 　\(Ly_{t}=y_{t-1}\)　　によって定めます。

ラグオペレーターは、時系列 \(y_t\) の時間 \(t\) を過去に遡る変換です。
数学では新しい記号が大きな威力を発揮します。
L は数列を１つずらす変換ですが、ずらさない変換を 1 と表し、
　　　\(1~y_{t}=y_{t}\)　と表すことにします。

隣接２項間の漸化式.
\(y_{t}=3-0.5 y_{t-1}\)　は、

　　　　\(1~y_{t}+0.5~L y_{t}=3\)　　　　\(\therefore~~(1+0.5~L) y_{t}=3\).

\(\phi=1+0.5~L\)　は数列をラグ変換する写像を表しラグ多項式と呼ばれます。.

一般項は漸化式を解くことにより、
　　　\(y_{t}=2+(y_0-2)(-0.5)^t\).
ラグ方程式 \(\phi=1+0.5~L=0\) の解は、
　　　\(L=\frac{1}{-0.5}\)
一般項に現れる公比の逆数であり、 ラグ方程式の解の絶対値が１より大きいとき数列は収束します。

隣接３項間の漸化式.
\(y_{t}=\phi_{1}y_{t-1}+\phi _{2}y_{t-2}\) は、
　 \(1y_{t}=y_{t}\) 、\(Ly_{t}=y_{t-1}\) 、\(L^{2}y_{t}=y_{t-2}\)

ラグ多項式を \(\phi \left( L\right) =1-\phi _{1}L-\phi _{2}L^{2}\)　と定義すると、

　　　\(\phi\left( L\right) y_{t}=y_{t}-\phi _{1}y_{t-1}-\phi _{2}y_{t-2}\).

\(y_{t}=2+0.4 y_{t-1}+0.3y_{t-2}\) は、

\(\phi=1-0.4~L-0.3~L^2\). 　

　　　\(\phi \left( L\right) y_{t}=3\)　と表されます。

　　　\(y_{t}-\phi_{1}y_{t-1}-\phi _{2}y_{t-2}=m+\varepsilon _{t}\). 　　　 \(y_{t}=m+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\varepsilon _{t}\).

AR(2) モデルはラグ多項式によって表され、

前後２季節までのラグ相関係数に着目し、気候変動を説明しようとしたのがウォーカーでした。このことは ユール・ウォーカー方程式を用いた時系列分析につながっていきます。時間をずらすことの重要性に注目したのが ウォーカーだと言えるでしょう。

そこで、時間をずらす操作をラグオペレーターとして定義し、定常となる条件を考えます。
まず、AR(1) から始めて、AR(2)へと拡張するのが自然でしょう。

時刻 t を k ずらした \(Y_{t-k}\) と、\(Y_{t}\) との相関係数を ラグ k の自己相関係数 といい、
　　　\(\rho_k=Cor(Y_t,~Y_{t-k})\) 　と書きます。

4 展望.

これまで、ARモデル、MAモデルについてみてきました。
実際のデータに適応し、予測性能を上げるためには更なるモデルの構築が必要です。

  • ARモデルとMAモデルを合わせたSARMAモデル.
  • 周期性を考慮したSARIMAモデル.

気象の予測などに活用されています。 これらは将来の学習に委ねます。