ARモデルの定常性の判定方法

時系列データを予測するために、簡単な確率モデルを考えます。モデルを考える上では、確率過程を使います。その中で、定常確率過程が大事です。今回の記事では、定常確率過程にも、非定常確率過程にもなるARモデルを紹介し、どんな時に非定常になるのか計算してみます。

参考文献は一応3冊あります。

https://amzn.to/2T39nSL

https://amzn.to/2umqXXt

https://amzn.to/2N7378V

記事で使っているソースコードはgithub に置いてあります。
https://github.com/msamunetogetoge

ARモデル
定常モデルと非定常モデル
別のアプローチ
python による実装
まとめ

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)モデルが定常だと仮定します。¹定常性の仮定の下²
$$\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}$$
この式から、自己相関係数が求まります。³
次に、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 $未満の時は、以下のように解きます。 ⁴ 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)
$
この、ベクトル方程式⁵を
$$\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) $は以下の式⁶で定義されます。
$$\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多項式の解が全て単位円の外にあること という条件が出てきます。⁷
一番簡単なAR(1) で確認してみましょう。AR(1)の特性方程式は以下のようになります。
$$\begin{eqnarray}
\lambda -\phi _1 =0
\end{eqnarray}$$
従って、$ |\phi _1 |<1 $が定常な確率過程になる条件です。
また、ユール・ウォーカー方程式から、$ \rho _k = \phi _1 ^k $なので、定常なAR(1)モデルの自己相関係数は、絶対値が指数関数的に小さくなります。⁸

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)でも同じことが起こります。