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)

mv <- c(mean(Yt),sd(Yt))
mv #平均と分散 
## [1] 1.906485 1.242185

平均は変動の中心、標準偏差は変動幅としての意味を持っていることが分かります。
青の点線の間が 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 2.093269 1.165108 1.712376 2.625682 3.021953
head(tail(Yt,-1))
## [1] 2.093269 1.165108 1.712376 2.625682 3.021953 1.521073

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

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

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

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

整理すると、

  • 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 の自己相関係数といいます。

公式.

　\(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=10)
## 
## Autocorrelations of series 'Yt', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.489  0.272 -0.133 -0.001  0.067 -0.062  0.034  0.031  0.006  0.107
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}\)

定常性の仮定をおくことにより、推定に \(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 ユール・ウォーカー方程式

定常性の仮定のもとでは、自己共分散は t に依存せず、
　　　　\(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}=c+\phi y_{t-1}+\varepsilon _{t}\) を用いると、
　　　\(Cov\left( y_{t},y_{t-k}\right)=Cov\left( c+\phi y_{t-1}+\varepsilon _{t},y_{t-k}\right)\)

\(Cov(\varepsilon_t, y_{t-k})=0\)　　\(Cov(c, y_{t-k})=0\)　に注意して、
　　　\(Cov\left( y_{t},y_{t-k}\right)=\phi Cov\left( y_{t-1},y_{t-k}\right)=\phi \gamma_{k-1}\)

よって、\(\gamma_{k}=\phi \gamma_{k-1}\)　（自己共分散の漸化式）

両辺を \(\gamma_0\)で割ると、
　　　　\(\rho_{k}=\phi \rho_{k-1}\)　（自己相関係数の漸化式）

これより、　\(\rho_{k}=\phi^k\)

自己相関係数の漸化式を、ユール・ウォーカー方程式といいます。
ユール・ウォーカー方程式は、実データから時系列モデルの係数を求めるために用いられます。

p過去までの影響を受ける AR(p) モデルの漸化式は、

　　　\(\displaystyle{ y_t = c+\phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_p y_{t-p}+\varepsilon_t}\)

　　　\(\displaystyle{ \rho_k = \phi_1\rho_{k-1} + \phi_2\rho_{k-2} + \dots + \phi_p\rho_{k-p}~~~( k=1,2,\cdots,p)}\)

このp個の連立方程式を解くことで、\(\phi_1,\phi_2 \cdots,\phi_p\)が得られます。

AR(2) を例にとれば、
\(\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\) で割ると、
　　　\(\rho_1=1\cdot \phi_{1} + \phi_{2}\rho_0\)
　　　\(\rho_2=\phi_{1}\rho_1 + 1\cdot \phi_{2}\)

3.2 定常 AR(1).

AR(1)モデル.

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

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

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

3.2.1 \(\phi=0.5\) のとき.

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2+0.5*yt+rnorm(1,0,1)
  Yt <- c(Yt,yt)
}
head(Yt)
## [1] 3.000000 3.988726 3.729293 4.404064 4.556515 4.038689
length(Yt)
## [1] 100
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),var(Yt))#平均と標準偏差
MV <- c(2/(1-0.5),1/(1-0.5^2))  
mv
## [1] 3.658915 1.072341
MV
## [1] 4.000000 1.333333

3.3 定常の意味.

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

  • \(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}\) .

上記の例では、平均は理論値とほぼ一致しています。

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

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

ということが言えます。

参考までに、誤差項を含まないシミュレーションを行います。

yt=3
Yt=c(3)
for (i in 1:(100-1)) {
  yt <- 2+0.5*yt
  Yt <- c(Yt,yt)
}
head(Yt)
## [1] 3.00000 3.50000 3.75000 3.87500 3.93750 3.96875
length(Yt)
## [1] 100
plot(Yt,type="l",xlab="time")
abline(h=4,col="red")

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

\(\phi=-0.5\) として同様のシミュレーションをすると、

次に、定常でない場合を考えます。

\(y_{t}=c+\phi y_{t-1}+\varepsilon _{t}\) の定常性に関わらず、

  • \(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.4 非定常 AR(1).

AR(1)モデル.

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

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

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

yt=2
Y=c(2)
for (i in 1:100) {
  yt <- 1.5*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
}
plot(Y,type="l",xlab="time")

mean(Y)
## [1] 2.773448e+16

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

3.4.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)
for (i in 1:100) {
  yt <- 1*yt+rnorm(1,0,1)
  Y <- c(Y,yt)
}
plot(Y,type="l",xlab="time")

mean(Y)
## [1] 11.92183

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

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

3.5 まとめ.

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

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

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

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

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\) の絶対値が１より小さければ収束するわけであり、それに誤差項が加わって複雑化しているだけなのです。

3.6 ユールとウォーカー.

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

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

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

参考サイト.