時系列データを予測するために、簡単な確率モデルを考えます。モデルを考える上では、確率過程を使います。その中で、定常確率過程が大事です。今回の記事では、定常確率過程にも、非定常確率過程にもなるARモデルを紹介し、どんな時に非定常になるのか計算してみます。
参考文献は一応3冊あります。
記事で使っているソースコードはgithub に置いてあります。
https://github.com/msamunetogetoge
ARモデル
定常確率過程を作れるMA(q)モデルを紹介する記事を書きました。 MA(q)モデルでは、q次までの相関係数の絶対値がパラメーター次第で それなりにコントロール出来ました。しかし、MAモデルの問題点は、高次(次数N)の相関係数を持つデータを生成するにはMA(N) モデルを考える必要がある事でした。その問題点を解決できるARモデル( autoregressive model ) を紹介します。
AR(1)や、AR(p)は、以下の式で表されるモデルです。\( \epsilon _t \sim { \rm W.N. } (\sigma ^2 ) \)で\( \{ \epsilon _t \} \)が分散\( \sigma ^2 \)のホワイトノイズに従う事を示します。
[AR(1)モデル]
$$\begin{eqnarray}
y_t&=& c+ \phi _1 y_{t-1} +\epsilon _t \\
\epsilon _t &\sim & { \rm W.N. } (\sigma ^2 )
\end{eqnarray}$$
[AR(p)モデル]
$$\begin{eqnarray}
y_t&=& c+ \sum_i ^p\phi _i y_{t-i} +\epsilon _t \\
\epsilon _t &\sim & { \rm W.N. } (\sigma ^2 )
\end{eqnarray}$$
AR(1) モデルを見ると、\(y_t \)の式に直接\(y_{t-1} \)が入っているので、高次まで自己相関がある事が分かると思います。同じ理由から、平均を取った時に無条件で\(t\)に寄らないかは怪しい所です。
MA(1)モデルは以下の式でした。
$$\begin{eqnarray}
y_t&=& c+ \theta _1 \epsilon _{t} +\epsilon _{t-1} \\
\epsilon _t &\sim & { \rm W.N. } (\sigma ^2 )
\end{eqnarray}$$
MA(1)では、\(y_t \)を生成するのは基本的にホワイトノイズなので、高次の自己相関がありませんでした。
AR(1)モデルの平均や分散を求めましょう。 計算をするために、AR(1)モデルが定常だと仮定します。1定常性の仮定の下2
$$\begin{eqnarray}
E[y_t]=\mu &=& c+ \phi _1 \mu \\
\mu &=& \frac{c}{1-\phi _1} \\
{\rm Var } [y_t] = \gamma _0 &=& \phi _1 ^2 \gamma _0 + \sigma ^2 \\
\gamma _0 &=& \frac{\sigma ^2 }{1-\phi _1 ^2} \\
\gamma _k &=& {\rm Cov } (y_t , y_{t-k} ) = {\rm Cov } ( \phi _1 y_{t-1} , y_{t-k} ) \\
\gamma _k &=& \phi _1 \gamma _{k-1}
\end{eqnarray}$$
最後の式は、\( \gamma _0 \)で割って、自己相関係数の方程式に書き直すことが出来ます。
$$\begin{eqnarray}
\rho _k = \phi _1 \rho _{k-1}
\end{eqnarray}$$
この差分方程式をユール・ウォーカー方程式( Yule–Walker equation )と呼びます。
AR(1)の場合は簡単に解けます。
$$\begin{eqnarray}
\rho _k = \phi _1 ^k
\end{eqnarray}$$
この式から、自己相関係数が求まります。3
次に、AR(p)モデルが定常だと仮定して、平均や自己相関係数を計算しましょう。定常と仮定してしまえば、期待値や分散の線形性で計算出来ます。
[AR(p)モデルの平均や分散]
$$\begin{eqnarray}
E[y_t] = \mu &=& c/(1- \sum_{i=1}^p \phi _i) \\
\gamma _0 &=& \sigma ^2 /(1- \sum_{i=1}^p \phi _i ^2 ) \\
\gamma _k &=& \sum _{i=1} ^p \phi _i \gamma _{k-i} \\
\rho _k &=& \sum _{i=1} ^p \phi _i \rho _{k-i}
\end{eqnarray}$$
最後の、自己相関係数に関する差分方程式も、 ユール・ウォーカー方程式( Yule–Walker equation )と呼びます。 \(\rho _k \)は、kがp未満か、p以上かの場合に分けて解くことが出来ます。
\(p \)未満の時は、以下のように解きます。 4 AR(3)の場合を考えておけば、pが3以上の場合でも同じ感じです。ユール・ウォーカー方程式から、以下の式が成り立ちます。
\(
\left(
\begin{array}{c}
\rho _1 \\
\rho _2
\end{array} \right)
=
\left(
\begin{array}{c}
\phi _1 \\
\phi _2
\end{array} \right)
+
\left(
\begin{array}{cc}
\phi _2 & \phi _3 \\
\phi _2 + \phi _3 & 0
\end{array} \right)
\left(
\begin{array}{c}
\rho _1 \\
\rho _2
\end{array} \right)
\)
\( \left(
\begin{array}{cc}
1- \phi _2 & -\phi _3 \\
-\phi _2 – \phi _3 & 1
\end{array} \right)
\left(
\begin{array}{c}
\rho _1 \\
\rho _2
\end{array} \right)
=
\left(
\begin{array}{c}
\phi _1 \\
\phi _2
\end{array} \right)
\)
この連立方程式を解けば、自己相関係数が出てきます。
左辺の行列が逆行列が持つかは怪しい訳ですが、定常性の仮定から逆行列を持つと思って計算すると、\(\rho _1 , \rho _2 \)は以下のようになります。
\(
\left(
\begin{array}{c}
\rho _1 \\
\rho _2
\end{array} \right)
=
\frac{1}{1-\phi _2 – \phi _3 (\phi _2 + \phi _3 )}
\left(
\begin{array}{c}
\phi _1 + \phi _2 \phi _3 \\
\phi _1 (\phi _2 + \phi _3 )+ \phi _2 (1 – \phi _3 )
\end{array} \right)
\)
\( \rho _3 \)以降は、ユール・ウォーカー方程式に逐一代入することで、求まります。
また、ユール・ウォーカー方程式から、φの値を求める事が出来ます。行列を書くのが大変なので、AR(2) の場合で考えます。\( \rho _0 =1 \)を使って、以下の方程式が成り立つ事に注目します。
\( \left(
\begin{array}{c}
\rho _1 \\
\rho _2 \\
\end{array}
\right) =
\left(
\begin{array}{cc}
1 & \rho _1\\
\rho _1 & 1 \\
\end{array} \right)
\left(
\begin{array}{c}
\phi _1 \\
\phi _2
\end{array} \right)
\)
この式を
$$\begin{eqnarray}
r=R\phi
\end{eqnarray}$$
と書くと、
$$\begin{eqnarray}
\phi =R^{-1} r
\end{eqnarray}$$
です。
定常モデルと非定常モデル
AR(1) モデルで、\( \phi _1 \)の値を変えて時系列データを生成すると、定常な確率過程と非定常な確率過程が表れる事が分かります。
どのような条件で定常と非定常が入れ替わるのか考えます。
AR(p)モデルを、AR(1)モデルみたいな見た目に書き換える事が出来ます。
\(
\left(
\begin{array}{c}
y_t \\
y_{t-1}\\
\cdots \\
y_{t-p+1}
\end{array}
\right)
=
\left(
\begin{array}{cccc}
\phi _1 & \phi _2 & \cdots & \phi _p \\
1& 0 & \cdots &0 \\
0& 1 & \cdots &0 \\
0 & \cdots &0 &1
\end{array}
\right)
\left(
\begin{array}{c}
y_{t-1} \\
y_{t-2}\\
\cdots \\
y_{t-p}
\end{array}
\right)
+ \left(
\begin{array}{c}
\epsilon _t \\
0\\
\vdots \\
0
\end{array}
\right)
+ \left(
\begin{array}{c}
c \\
0\\
\vdots \\
0
\end{array}
\right)
\)
この、ベクトル方程式5を
$$\begin{eqnarray}
\xi _t = F \xi _{t-1} + \epsilon _t +c
\end{eqnarray}$$
と書き直せば、AR(1)モデルみたいになります。右辺の\( \xi _{t-1} \)に対しても同じような式が成り立つので、
$$\begin{eqnarray}
\xi _t = F^t \xi_0 +F \epsilon _{t-1} \cdots + F^{t-1} \epsilon _{1} + F^t \epsilon _0 +c
\end{eqnarray}$$
が成り立ちます。ここで、Fという行列について、少し思いをはせましょう。
Fの固有値を計算して、対角行列に出来る事が分かれば、
\(F= \Gamma ^{-1} \Lambda \Gamma \)と書けます。ただし、\( \Lambda = {\rm diag}[\lambda _1 , \cdots , \lambda _p ]\)は対角行列で、固有値が並んでいます。この時、
$$\begin{eqnarray}
F^t =\Gamma ^{-1} \Lambda ^t \Gamma
\end{eqnarray}$$
となり、\( \Lambda \)の対角成分に 1以上の成分があると、 \( t\rightarrow \infty \)の極限で、\( \xi _t \)の全ての成分が発散します。
この事から、\(\{ y_t \} \)の平均が時間によって変化するようになる事が分かります。
つまり、AR過程が定常か、非定常かは\(F \)の固有値の値にかかっています。
という訳で、\(F-\lambda I \)の行列式を考えましょう。
適当なpで計算してみてほしいのですが、1行目で余因子展開して行列式を計算していくと、p-1列目までは、下三角行列、p-1列目と、p列目では上三角行列が出てきます。そしてi列目の計算では下/上三角行列の対角成分には、p-i個の\( -\lambda \)と、1だけが出てきます。
$$\begin{eqnarray}
{\rm det }(F- \lambda I ) &=& ( \phi _1 -\lambda ) (- \lambda )^p -\phi _2 (-\lambda )^{p-1} +\cdots +(-1)^{1+p-1} \phi _{p-1} (-\lambda ) + (-1)^p \phi _p \\
&=& (-1)^{p+1} ( \lambda ^p -\phi _1 \lambda ^{p-1} – \cdots – \phi _{p-1} \lambda – \phi _p )
\end{eqnarray}$$
\( {\rm det }(F- \lambda I ) =0 \)の事を、AR特性多項式と呼んだりします。ARモデルが定常か否かは、AR多項式の根\( \lambda _1 , \cdots , \lambda _p \)の絶対値が1より大きいかどうかで決まります。
議論をまとめておきます。
[AR過程の定常性]
AR特性多項式
$$\begin{eqnarray}
\lambda ^p -\phi _1 \lambda ^{p-1} – \cdots – \phi _{p-1} \lambda – \phi _p =0
\end{eqnarray}$$
の解\( \lambda _1 , \cdots , \lambda _p \)の中に、\( | \lambda _i | \geq 1\)を満たすものがある時、AR(p)モデル
$$\begin{eqnarray}
y_t = c+ \sum_i ^p\phi _i y_{t-i} +\epsilon _t
\end{eqnarray}$$
は非定常確率過程を生成する。
逆に、任意の\(i \)について、 \( | \lambda _i | < 1\) が成り立つとき、AR(p)は定常確率過程を生成する。
人によっては、AR特性多項式を\(\lambda ^p \)で割って、
$$\begin{eqnarray}
1-\phi _1 \lambda ^{-1} – \cdots – \phi _{p-1} \lambda ^{-p + 1}- \phi _p \lambda ^{-p } =0
\end{eqnarray}$$
の事をAR多項式と呼ば場合があります。この場合、 \( \lambda \)を\(z\neq 0 \)の複素数だと思うと、AR(p)が定常確率過程であるための条件は、
AR多項式の解が全て単位円の外にあること
と言い換える事が出来ます。
別のアプローチ
\(\lambda ^{-1} \)を係数にした多項式は、 Z変換を考えると出てきます。
Z変換は、フーリエ変換とかラプラス変換の一般化です。
数列\( \{x_n \} \)のZ変換\( F_x (z) \)は以下の式6で定義されます。
$$\begin{eqnarray}
F_x (z) = \sum_{n=-\infty} ^{\infty} x_n z^{-n}
\end{eqnarray}$$
Z変換は線形性と、\( \{x_n \} \)の添え字をずらすと、\(z^{-1} \)が出てくる性質があります。
$$\begin{eqnarray}
F_{ax+by} (z) &=& aF_x (z) + bF_y (z) \\
F_{y_{t} }(z) &=& \sum y_t z^{-t} = z^{-1} \sum y_t z^{-t+1} = z^{-1} F_{y_{t+1}} (z)
\end{eqnarray}$$
また、この変換には逆変換があり、Z変換した後何かした後に、元の関数の空間に戻ってくることが出来ます。\(c=0 \)のときのAR(p)モデルを書き換えてZ変換すると、AR特性多項式が出てきます。
$$\begin{eqnarray}
y_t – \sum \phi _i y_{t-i} =\epsilon _t \\
F_{y_t – \sum \phi _i y_{t-i} }(z) = F_{ \epsilon _t } (z)\\
(1 – \sum_i ^p \phi _i z^{-i} )F_{y_t} (z) = F_{ \epsilon _t } (z)
\end{eqnarray}$$
最後の式で、\( F_{y_t} (z) \) の係数が、AR特性多項式になっています。両辺をAR特性多項式で割った時に変な事が起きない為に、 AR多項式の解が全て単位円の外にあること という条件が出てきます。7
一番簡単なAR(1) で確認してみましょう。AR(1)の特性方程式は以下のようになります。
$$\begin{eqnarray}
\lambda -\phi _1 =0
\end{eqnarray}$$
従って、\( |\phi _1 |<1 \)が定常な確率過程になる条件です。
また、ユール・ウォーカー方程式から、\( \rho _k = \phi _1 ^k \)なので、定常なAR(1)モデルの自己相関係数は、絶対値が指数関数的に小さくなります。8
python による実装
定常/非定常な時系列データを生成してみます。また、 定常過程の時、自己相関係数が指数的に減少することを確かめます。
初めにAR(1)で、\( \phi _1 \)を調整して 確率過程を生成してグラムを描いてみます。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
phi = np.round(np.linspace(-1.01, 1, 4),2)
epsilon = np.random.randn(4,10000)
T=100
n=len(phi)
Y=np.zeros((n,T))
Y[:,0]=0
for i in range(1,T):
Y[:,i]=phi*Y[:, i-1] +epsilon[:,i]
plt.figure(figsize=(10,10))
plt.subplot(2,2,1)
plt.plot(Y[0,:])
plt.title("AR(1)φ = {}".format(phi[0]))
plt.subplot(2,2,2)
plt.plot(Y[1,:])
plt.title("AR(1)φ = {}".format(phi[1]))
plt.subplot(2,2,3)
plt.plot(Y[2,:])
plt.title("AR(1)φ = {}".format(phi[2]))
plt.subplot(2,2,4)
plt.plot(Y[3,:])
plt.title("AR(1)φ = {}".format(phi[3]))AR(1)が定常か否かは、
$$\begin{eqnarray}
y_t&=& c+ \phi _1 y_{t-1} +\epsilon _t \\
\epsilon _t &\sim & { \rm W.N. } (\sigma ^2 )
\end{eqnarray}$$
において、\( | \phi _1 |<1 \)か否かでした、\( \phi _1 =1 \)(右下のグラフ)はありがちな時間によって大きくなるデータのように見えます。一方で、\( \phi _1 =-1.01 \)(左上のグラフ)は、振幅が大きくな振動のグラフです。定義式から、前のデータに\(\phi _1 \)を掛けるのが基本なので納得出来ると行くと思います。
\( |\phi _1 |<1 \)を満たすデータは、平均値があまり変わらないように見えるので、定常っぽいですね。
次に、生成したデータを使ってコレログラムを描いてみます。定常なAR過程では、自己相関係数の絶対値が指数関数的に減少していきます。
def autcor(X, k):
X_bar = np.mean(X)
X_bar_k = X_bar * np.ones(k)
X_k = np.append(X_bar_k ,X)
k_X = np.append(X, X_bar_k )
r= np.dot(X_k -X_bar , k_X - X_bar )
r/= np.linalg.norm(X-X_bar)**2
return r
n=30
Aut_r =np.zeros((4,n))
for j in range(4):
for i in range(n) :
Aut_r[j-1,i] = autcor(Y[j,:], i)
plt.figure(figsize=(6,10))
plt.subplot(2,1,1)
plt.stem(Aut_r[1,1:])
plt.title("Autocorrelation of AR(1) φ={}".format(phi[1]))
plt.subplot(2,1,2)
plt.stem(Aut_r[2,1:])
plt.title("Autcorrelation of AR(1) φ={}".format(phi[2]))
\( \phi =-0.4 \)では良く分かりませんが、\( \phi =0.3 \)はまさに、という感じです。
まとめ
- ARモデルの定義を確認した
- AR(p)モデルの平均や分散を計算した
- AR(p)モデルの定常/非定常の条件を計算した
- グラフを描いて計算結果や、定常/非定常確率過程のグラフを確認した
- 定常のための条件は後で計算します。
- \(E[y_t ] =\mu , {\rm Cov } (y_t , y_s) = \gamma _{t-s} \)という事です。
- 定常である事と、\( | \phi _1 | <1 \)が同値なので、高次になるにつれてAR(1)モデルの自己相関係数は、指数関数的に小さくなります。
- 計算の際には、\( {\rm Cov} (x,y) = {\rm Cov} (y,x) \)より、 \(\gamma _k = \gamma _{-k} \)に注意しましょう。
- 1行目が\(y_t \)の式を表しており、2行目以降は\( y_{t-i} = y_{t-i} \)を表しています。
- \( z\neq 0 \)上で定義された複素関数を生成します
- Ar(p)の特性方程式の根の形が分からないので、解き方があれば教えてください
- AR(p)でも同じことが起こります。