1 はじめに.

高校情報Ⅰでは時系列分析が取り上げられています。時系列データは身近に見られる統計データであり、時系列解析は未来を予測する手段として重要です。
時系列解析の背景には確率過程の理論があり、その数理が高校で扱われることはありません。しかし、時系列モデルは数列の漸化式の応用として扱うことができます。
ここでは、漸化式を用いた時系列モデルでデータを生成し、その性質に触れることを目的とします。

数理的には「定常性：平均と自己共分散が t に依存しない」という説明から入るべきですが、ここでは具体的な時系列モデルから入り、そのモデルを通して定常性の理解を深めます。

2 AR(1)モデル.

自己回帰過程 (auto regressive process) は過去のデータから現在のデータを説明するもので、ARモデルと呼ばれます。１過去（１つ前）データの一次式で表されるモデルを AR(1) モデルといいます。AR(1)モデルは、隣接２項間の漸化式に誤差項 \(\varepsilon _{t}\) を加えた式で表されます。

　　　\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\).

\(\varepsilon _{t}\) は誤差項であり、 \(N(0,\sigma^2)\)に従います。

初期値を定めれば、繰り返し構文を用いて時系列データを生成することができます。
簡単のため、誤差項は \(N(0,1)\) に従うとします。

例　\(y_1=3,~~~y_{t}=3-0.5 y_{t-1}+\varepsilon _{t}\)

yt=3 #初期値
Yt=c(3)
for (i in 1:200) {
  yt <- 3-0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
plot(Yt,type="l",xlab="time")
abline(h=mean(Yt),col="red")
abline(h=mean(Yt)+1.96*sd(Yt),col="blue", lty=2)
abline(h=mean(Yt)-1.96*sd(Yt),col="blue", lty=2)

＜平均と分散＞

AR(1) モデル：\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\) に対し、

  • \(E(y_{t})=c+\phi E(y_{t-1})+E(\varepsilon _{t})\).
  • \(V(y_{t})=\phi^2 V(y_{t-1})+\sigma^2\).

もし、\(E(y_{t})=E(y_{t-1})=\mu\)、\(V(y_{t})=V(y_{t-1})=V\)　であれば、
（この条件は、変動が安定していることを意味します）

  • \(\mu=c+\mu~~~~~~~~~~\therefore~\mu=\dfrac{c}{1-\phi}\)
  • \(V=\phi^2 V+\sigma^2~~~\therefore~V=\dfrac{\sigma^2}{1-\phi^2}\) .

上記の例では、

  • \(\mu=\dfrac{c}{1-\phi}=\dfrac{3}{1+0.5}=2\)
  • \(V=\dfrac{\sigma^2}{1-\phi^2}=\dfrac{1}{1-(-0.5)^2}\fallingdotseq1.33~~~\therefore sd\fallingdotseq1.15\) 　　　
mv <- c(mean(Yt),sd(Yt))
mv #平均と分散 
## [1] 1.929185 1.118101

実データの平均と標準偏差は理論値とほぼ一致しています。
また、変動のグラフから、「平均は変動の中心、標準偏差は変動幅」としての意味を持っていることが分かります。
青の点線の間が 95% 信頼区間であり、変動の 95% がこの範囲内に入ります。

2.0.1 自己相関係数.

自己相関係数は、過去データが現在に及ぼす影響の度合いを表します。

１過去（１つ前）データとの自己相関係数 \(\rho_1\) は、
　　　\(\rho_1=Cor(Y_t,~Y_{t-1})\)　ですが、
　　　\(\rho_1=Cor(Y_t,~Y_{t+1})\)　でも計算することができます。

上記の例で、 \(Y_t\) から最初の１項を除けば、\(Y_{t+1}\) になります。
最初の１項を除く操作は、tail(Yt,-1) です。 はじめの６項を head()関数で表示すると、

head(Yt)
## [1] 3.000000 1.067379 4.302638 1.952458 2.858613 1.005472
head(tail(Yt,-1))
## [1] 1.067379 4.302638 1.952458 2.858613 1.005472 2.712663

tail(Yt,-1)によって、項が１つスライドしていることが分かります。

データの個数を揃えるため head(Yt,-1) によってYt の最後の１項を除き、相関係数を計算します。

r1 <- cor(head(Yt,-1),tail(Yt,-1))
r1
## [1] -0.3830295

これは、ラグ１（ズレ１）の標本自己相関係数で、漸化式の係数に近い値となっています。
（理論上、ラグ１の自己相関係数は漸化式の係数に一致します）.

整理すると、

  • head(Yt,-k) はYtから最後のk項を除いたもの.
  • tail(Yt,-k) はYtから最初のk項を除いたもの.

ラグ k（ズレ k）の標本自己相関係数は、
　　　cor(head(Yt,-k), tail(Yt,-k))　によって計算されます。

課題.

　ラグ２、ラグ３の自己相関係数 \(\rho_2,\rho_3\) を計算しましょう。
　［注意］求めた値は標本相関係数であり真値ではありません。

自己相関係数の定義.

　k過去（k個前）データとの自己相関係数 \(\rho_k\) は、
　　　\(\rho_k=Cor(Y_t,~Y_{t-k})\)　であり、ラグ k の自己相関係数といいます。

AR(1) の自己相関係数.

　\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\)　のとき、　\(\rho_1=\phi、~\rho_k=\phi^k\).

これを数理的に示すには定常という性質が必要になりますが、これは以下で説明します。

　　　 .
　参考サイト　AVILEN様：1次ARモデルの特徴や統計量について .

2.0.2 コレログラム.

Rにはacf関数(autocorrelation function) が用意されており、自己相関係数を計算し図示することができます。

  • plot=F で数値を表示し、
  • plot=T でグラフを表示します.

グラフをコレログラムと呼びます。

#自己相関係数
acf(Yt,plot=F,lag.max=7)
## 
## Autocorrelations of series 'Yt', by lag
## 
##      0      1      2      3      4      5      6      7 
##  1.000 -0.381  0.141 -0.047 -0.032 -0.028  0.021 -0.106
acf(Yt,plot=T)

青の点線で挟まれた範囲は、無相関という帰無仮説が棄却されない範囲です。この帯の間にグラフがあれば、そのラグで相関はないと判断されます。
上図では、はっきりした影響を及ぼすのは１過去までですが、２過去も若干影響が残っているようです。

2過去のデータは１過去のデータを通じて現在に影響を及ぼします。この間接的な影響を取り除いた自己相関係数を偏自己相関係数といいます。
acf()に、type=“p”（partialの意味）を付け加えれば、偏自己相関係数がプロットされます。

acf(Yt,plot=T,type="p")

間接的影響を除くと、直接影響を与えるのは１過去のデータのみであることが分かります。

2.0.3 漸化式との関係.

誤差項のない隣接２項間の漸化式は、
　　　\(y_{t}=3-0.5 y_{t-1}\).
一般項は漸化式を解くことにより、
　　　\(y_{t}=2+(y_0-2)(-0.5)^t\).
変化の様子は、

t=1:10
y=2+(3-2)*(-0.5)^t
y
##  [1] 1.500000 2.250000 1.875000 2.062500 1.968750 2.015625 1.992188 2.003906
##  [9] 1.998047 2.000977
plot(y,type="l",ylim=c(1.8,2.2),xlab="time")
abline(h=2,col="lightblue")

漸化式の一般項は２に収束し、変動もシンプルです。
時系列では誤差項が加わるため、変動が複雑化していることが分かります。

3 定常性.

一般の時系列データでは、平均と自己共分散は時刻 t に依存します。
　　　\(E\left( y_{t}\right) =\mu_t\)　　\(Cov\left( y_{t},y_{t-k}\right)=\gamma _{kt}\)

時系列データの定常性とは、平均と自己共分散が t に依存しないことをいいます。
　　　\(E\left( y_{t}\right) =\mu\)　　\(Cov\left( y_{t},y_{t-k}\right)=\gamma _{k}\)

\(k=0\) のとき、\(Cov\left( y_{t},y_{t}\right)=\gamma _{0}\) は分散を意味します。
分散は時刻に依存しないので、
　　　\(V(y_t)=Cov\left( y_{t},y_{t}\right)=\gamma _{0}\)
　　　\(V(y_{t-k})=Cov\left( y_{t-k},y_{t-k}\right)=\gamma _{0}\)

ラグ k の自己相関係数を \(\rho_k\) とすると、
　　　\(\rho_k=\dfrac{Cov\left( y_{t},y_{t-k}\right)}{\sqrt{V(y_t)}\sqrt{V(y_{t-k})}}=\dfrac{\gamma_k}{\sqrt{\gamma_0}\sqrt{\gamma_0}}=\dfrac{\gamma_k}{\gamma_0}\)

＜標本統計量＞
定常性の仮定をおくことにより、推定に \(y_{t},y_{t-1},\)… を用いる事が可能になります。
即ち、標本値から平均や自己共分散などの統計量を求めることが妥当性を持つようになります。

　時点 \(t=1,2,\cdots,T\) までのデータの集合 \(\left(y_1,y_2,\cdots,y_T\right)\) に対し、
　標本平均　　　　\(\displaystyle \overline{y}=\dfrac{1}{T}\sum ^{T}_{t=1}y_{t}\)

　標本自己共分散 　\(\widehat{\gamma }_{k}=\dfrac{1}{T}\sum ^{T}_{t=k+1}\left( y_{t}-\overline{y}\right) \times \left( y_{t-k}-\overline{y}\right)\)

　標本自己相関係数 　　\(\widehat{\rho}_{k}=\dfrac{\widehat{\gamma}_{k}}{\widehat{\gamma}_{0}}\)

　（\(\gamma _{0}=Cov\left( y_{t},y_{t}\right)=V\left( y_{t}\right)\)　なので、\(\widehat{\gamma}_{0}\) は標本分散を表します）

　これらの標本値は T が大きくなれば真値に近づきます。

前節で説明した自己相関係数も標本値を用いており、定常性を前提としています。
定常性は、データからモデルを推定する際に前提となる条件です。

トレンドや季節変動は、統計量が時刻 t に依存するので定常ではありません。

　　（時系列データ）＝（トレンド）＋（季節変動）＋（定常変動）

この分解例のように、一般の時系列データは定常でない成分を含んでおり、定常ではありません。
分析の際は、トレンドや季節変動を除去しデータを定常化する必要があります。
（非定常な時系列データの分析については、ここでは立ち入りません）.

https://qiita.com/nori0__/items/ed279f0f82616f9665e2.
https://otexts.com/fppjp/regression-intro.html

セールスアナリティックス：時系列データの「定常性」と「3つの非定常性」. SARIMAモデル. 　

3.1 定常 AR(1).

　　\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\)　、誤差変動：\(\varepsilon _{t}\sim N(0,\sigma^2)\).

は、\(|\phi|<1\) のとき定常となります。

定常性を満たせば、

  • \(E(y_{t})=c+\phi E(y_{t-1})+E(\varepsilon _{t})\) から、　\(\mu=c+\mu~~~\therefore~\mu=\dfrac{c}{1-\phi}\)
  • \(V(y_{t})=\phi^2 V(y_{t-1})+\sigma^2\)から、　\(V=\phi^2 V+\sigma^2~~~\therefore~V=\dfrac{\sigma^2}{1-\phi^2}\) .

このことをシミュレーションで確かめましょう。
簡単のため、誤差項は \(N(0,1)\) に従うとします。

3.1.1 \(c=2,~~\phi=0.5\) のとき.

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

yt=3
Yt=c(3)
for (i in 1:(200-1)) {
  yt <- 2+0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
head(Yt)
## [1] 3.000000 2.892613 3.035954 3.357312 3.437594 5.791026
length(Yt)
## [1] 200
plot(Yt,type="l",xlab="time")
abline(h=4,col="red")
abline(h=mean(Yt)+1.96*sd(Yt),col="blue", lty=2)
abline(h=mean(Yt)-1.96*sd(Yt),col="blue", lty=2)

＜平均と標準偏差＞.
　実データによる標本値と理論値を比較します。

mv <- c(mean(Yt),sd(Yt))
mv #標本値
## [1] 4.018626 1.223233
MV <- c(2/(1-0.5),sqrt(1/(1-0.5^2)))  
MV #理論値
## [1] 4.000000 1.154701

標本値と理論値はほぼ一致しています。

3.2 定常の意味.

定常であることは、時系列データが収束することを意味しません。
分散が一定の値になるということは、平均を中心として一定の幅で変動することを意味します。

  • 定常かどうかは、誤差項のない漸化式で定義される数列に現れる公比によって決まり、
  • 誤差項の分散は、平均を中心とした変動の幅を決定する.

ということが言え、
　　　平均：\(\mu=\dfrac{c}{1-\phi}\)（変動の中心）、　　分散：\(V=\dfrac{\sigma^2}{1-\phi^2}\)（変動の幅）
となります。

誤差項を含まないシミュレーションを並置して比べます。

3.2.1 \(c=2,~~\phi=0.5\) のとき.

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

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2+0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Yt,type="l",xlab="time",main="誤差項あり")
abline(h=4,col="red")

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2+0.5*yt
  Yt <- c(Yt,yt)
}
plot(Yt,type="l",xlab="time",main="誤差項なし")
abline(h=4,col="red")

時系列データは単調かつ急速に収束します。
誤差項を加えることで、時系列データは平均を挟んだ変動となっています。

3.2.2 \(c=2,~~\phi=-0.5\) のとき.

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

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2-0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Yt,type="l",xlab="time",main="誤差項あり")
abline(h=2/1.5,col="red")

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2-0.5*yt
  Yt <- c(Yt,yt)
}
plot(Yt,type="l",xlab="time",main="誤差項なし")
abline(h=2/1.5,col="red")

係数がマイナスだと振動しますが、急速に収束します。
誤差項を加えるた時系列データは平均を挟んで激しく振動しています。

3.3 非定常 AR(1).

　　\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\)　、誤差変動：\(\varepsilon _{t}\sim N(0,\sigma^2)\).

は、\(|\phi|\ge 1\) のとき非定常となります。

  • \(E(y_{t})=c+\phi E(y_{t-1})\) から、　\(E(y_t)=\alpha+\phi^{t-1}\left(E(y_0)-\alpha\right)\)
    （\(\alpha\) は \(x=c+\phi x\) の解）

  • \(V(y_{t})=\phi^2 V(y_{t-1})+\sigma^2\)から、　\(V(y_{t})=\beta+\phi^{2(n-1)}\left(V(y_{0})-\beta)\right)\) .
    （\(\beta\) は \(x=\phi^2 x+\sigma^2\) の解）

平均も分散も、\(|\phi|\ge 1\) のとき発散します。

シミュレーションで確かめましょう。

3.3.1 \(c=2,~~~\phi=1.5\) のとき.

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

yt=3
Y=c(2)
for (i in 1:20) {
  yt <- 2+1.5*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Y,type="l",xlab="time",main="誤差項あり")

yt=3
Y=c(2)
for (i in 1:20) {
  yt <- 2+1.5*yt
  Y <- c(Y,yt)
}
plot(Y,type="l",xlab="time",main="誤差項なし")

この場合、時系列データは発散します。

yt=3
Y=c(2)
for (i in 1:200) {
  yt <- 2+1.5*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
}

3.3.2 \(c=0,~~~\phi=1\) のとき.

　　\(y_{t}=y_{t-1}+\varepsilon _{t}\)　、誤差変動：\(\varepsilon _{t}\sim N(0,\sigma^2)\)　のとき.

　　\(y_{t}=\varepsilon _{1}+\varepsilon _{2}+\cdots +\varepsilon _{t}\)　（誤差項の和）

これは、ランダムウォークと呼ばれます。平均値は収束しません。
時系列データと平均変動を並置します。

yt=0
Y=c(0)
M=c()
for (i in 1:50) {
  yt <- 1*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
  M <- c(M,mean(Y))
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Y,type="l",xlab="time",main="時系列データ")

plot(M,type="l",xlab="time",main="平均変動")

平均もランダムに変動し、一定した平均値を持ちません 。

3.3.3 \(c=2,~~~\phi=1\) のとき.

　　\(y_{t}=2+y_{t-1}+\varepsilon _{t}\)　、誤差変動：\(\varepsilon _{t}\sim N(0,\sigma^2)\)　のとき.

　　\(y_{t}=2t+\varepsilon _{1}+\varepsilon _{2}+\cdots +\varepsilon _{t}\)　（誤差項の和）

この場合、増加トレンドになります。

yt=2
Y=c(2)
M=c()
for (i in 1:50) {
  yt <- 2+1*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
  M <- c(M,mean(Y))
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(Y,type="l",xlab="time",main="時系列データ")
plot(M,type="l",xlab="time",main="平均変動")

課題. c=-2 のとき、どのようになるかシミュレーションしましょう。

（c の値によって、増加トレンドになるか減少トレンドになるかが決まります）

3.4 まとめ.

非定常を見れば、定常の意味がわかります。

  • 非定常は変動が発散することを意味します.
  • 定常は一定の幅に変動が収まり続けることを意味します.

定常であっても誤差項のため変動幅は０にはなりません。
ある変動幅の中で変動し続けることになります。

  • 定常であれば、その事象は変動しながらも存続し続け、
  • 定常でなければ、その事象は破綻し終わりを迎えます。

AR(1)モデルの定常性の判定について、

  • \(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\)　が定常であるかどうかは、
  • 漸化式　\(y_{t}=c+\phi y_{t-1}\)　が収束するかどうかによって決まります.
  • 即ち、　\(y_{t}=\phi y_{t-1}\)　が収束するかどうかによって決まります.

漸化式の公比 \(\phi\) の絶対値が１より小さければ収束するわけであり、それに誤差項が加わって複雑化しているだけなのです。

4 ユール&ウォーカー.

イギリスの統計学者 ウドニー・ユール（Udny Yule）は、太陽黒点の変動周期について、2次の自己回帰モデル AR(2) を用い、太陽黒点の変動周期がおよそ11年であることを示しました。

イギリスの数学者・気象学者 ギルバート・ウォーカー（Gilbert Walker） は、モンスーン予報やエルニーニョ・南方振動など地球規模の気象予測に自己相関係数が深く関わっていることに注目し、気候変動の予兆を探ろうとしました。熱波を予測し飢饉に備えることは重要なことだからです。

時系列分析は株価や物価など短期的な変動に用いられている印象ですが、気象変動や太陽の黒点活動など地球物理学的な長期の変動の分析から生まれていることは興味深いことです。気象変動や太陽の活動がある一定の幅の範囲内で続いてくれなければ（定常でなければ）地球の存続は危ういでしょう。

参考サイト.

4.1 時系列のモデリング.

ユールは、どのようにして太陽の黒点変動の自己回帰モデルを発見したのでしょうか？
具体的な資料はネットでは見当たりません。
ユールが、11年間、22年間、黒点を観察し続けたとは思えません。
そもそも、数学の発想は短期間の観察から遠い未来を予測することにあります。

時系列分析の一般的な流れは次のようなものです。

＜時系列分析の流れ＞.

  • あるモデルを仮定する.
  • 観測によって得られたデータからモデルの係数を定め、モデルの式を確定する。
  • モデルの当てはまりを検証する。

モデルの係数を定めるのに考えられたのが、ユール・ウォーカー方程式です。
ユール・ウォーカー方程式は、自己共分散を用いた連立方程式です。

4.2 ユール・ウォーカー方程式

実際のデータから自己共分散の標本値を求めることができます。
AR(p) モデルは、ラグに応じた自己共分散、
　　　\(\gamma _{0},\gamma _{1},\cdots,\gamma _{p}\)　　を持ちます。
これらの関係式に自己共分散の標本値を代入し、連立方程式としたものがユール・ウォーカー方程式です。

式の構成には共分散の性質を用います。少し難しいですが挑戦してみてください。

＜共分散の性質＞.

  • \(Cov(kX,Y)=k~Cov(X,Y)\).
  • \(Cov(X+Y,Z)=Cov(x,Z)+Cov(Y,Z)\).
  • X,Y が独立のとき、\(Cov(X,Y)=0\).
    特に定数c に対し、\(Cov(X,c)=0\).

まず AR(1) モデルについて考え、次に AR(2) モデルについて考えます。

　　\(Cov(c,~y_{t-k})=0\)　、\(k\geq 1\) のとき \(Cov(\varepsilon_t,~y_{t-k})=0\)　に注意して計算します。

＜AR(1) モデル＞.

　　　\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\) を用いると、

\(\begin{split} \gamma _{0}&=Cov\left( y_{t},y_{t}\right) =Cov\left( c+\phi y_{t-1}+\varepsilon _{t},y_{t}\right)\\ &=\phi Cov\left( y_{t-1},y_{t}\right)+Cov\left( \epsilon_{t},y_{t}\right)\\ &=\phi \gamma_{1}+Cov\left(\varepsilon _{t},c+\phi y_{t-1}+\varepsilon _{t}\right)\\ &=\phi \gamma_{1}+\sigma^2 \end{split}\).

\(\begin{split} \gamma _{h} &=Cov\left( y_{t},y_{t-h}\right) =Cov\left( c+\phi y_{t-1}+\varepsilon _{t},y_{t-h}\right)\\ &=\phi Cov\left( y_{t-1},y_{t-h}\right)+Cov\left( \epsilon_{t},y_{t-h}\right)\\ &=\phi \gamma_{h-1} \end{split}\).
\(\therefore~~~\gamma_{1}=\phi \gamma_{0}\)　（ユール・ウォーカー方程式）

連立方程式

　　　\(\gamma _{0}=\phi \gamma_{1}+\sigma^2\)、　　\(\gamma_{1}=\phi \gamma_{0}\)

に標本から求めた \(\gamma _{0},\gamma _{1}\) を代入すれば、\(\phi,~\sigma^2\) が求まり、モデルが確定します。

自己相関係数から求める方法もあります。
　　　\(\gamma_{h}=\phi \gamma_{h-1}\) の両辺を \(\gamma_0\)で割ると、 \(\rho_{h}=\phi \rho_{h-1},~~\rho_1=\phi\)

（自己相関係数の標本値からパラメーター\(\phi\) を求めることができます）
　

＜AR(2) モデル＞

　　　\(y_{t}=c+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\varepsilon _{t}\)　　を用いると、 .

\(\begin{split} \gamma _{0}&=Cov\left( y_{t},y_{t}\right) =Cov\left( c+\phi_1 y_{t-1}+\phi_{2}y_{t-2}+\varepsilon _{t},y_{t}\right)\\ &=\phi_1 Cov\left( y_{t-1},y_{t}\right)+\phi_2 Cov\left( y_{t-2},y_{t}\right)+Cov\left( \epsilon_{t},y_{t}\right)\\ &=\phi_1 \gamma_{1}+\phi_2~\gamma_{2}+\sigma^2\\ \end{split}\).

\(\begin{split} \gamma_1 &= {Cov}(y_{t}, y_{t-1})\\ &= {Cov}(c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2}+ \varepsilon_t, y_{t-1})\\ &=\phi_{1}{Cov}(y_{t-1},y_{t-1}) + \phi_{2}{Cov}(y_{t-2}, y_{t-1})\\ &=\phi_{1}\gamma_0 + \phi_{2}\gamma_1 \end{split}\)

\(\begin{split} \gamma_2 &= {Cov}(y_{t}, y_{t-2})\\ &= {Cov}(c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2}+ \varepsilon_t, y_{t-2})\\ &=\phi_{1}{Cov}(y_{t-1},y_{t-2}) + \phi_{2}{Cov}(y_{t-2}, y_{t-2})\\ &=\phi_{1}\gamma_1 + \phi_{2}\gamma_0 \end{split}\)

標本から求めた \(\gamma _{0},\gamma _{1},\gamma _{2}\) を代入すれば、 \(\phi_1,~\phi_2,~\sigma^2\) が求まり、モデルが確定します。

自己相関係数を用いた連立方程式をユール・ウォーカー方程式ということもあります。
上の３つの式の両辺を \(\gamma_0\) で割ると、
　　　\(1=\phi_{1}\rho_1 + \phi_{2}\rho_2+ \dfrac{\sigma^2}{\gamma_0}\)
　　　\(\rho_1=1\cdot \phi_{1} + \phi_{2}\rho_1\)
　　　\(\rho_2=\phi_{1}\rho_1 + 1\cdot \phi_{2}\)

観測によって得られた自己相関係数の標本値を代入してパラメーターと誤差分散を求めればモデルが確定し、シミュレーションや予測が可能になります。即ち、

　　\(\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}\) 　　からパラメータを求め、
　
　　　\(\sigma^2=\gamma_0\cdot (1-\phi_{1}\rho_1 - \phi_{2}\rho_2)\)　　から誤差分散を求めます。

AR(2) については、後の節でも扱います。

参考サイト.
　アベリオシステム数理技術チーム様 「時系列解析入門 2. ARモデル」.