時系列データ解析 06 ～自己回帰過程

www.asakura.co.jp

　沖本先生の本の続きです。MA過程の次はAR過程です。

　自己回帰(AR)過程は、過程が自己の過去に回帰されたモデルで表現されるものです。１次のAR過程、AR(1)は、

$\displaystyle y_{t} = c + \phi_{1} y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{2})$

で定義され、 $y_{t} \sim AR(1)$ と表記されます。

　AR過程は代表的な弱定常過程である

$\displaystyle y_{t} = \mu + \epsilon_{t}$

に、 $\phi_{1} y_{t-1}$ が追加された形をしていて、そのため $y_{t}$ と $y_{t-1}$ が相関をもつことになります。

　AR過程の確率的変動は、MA過程と同様、撹乱項のホワイトノイズ $\epsilon_{t}$ によって決まります。 $\epsilon$ が決まり、続いて $y$ が決まっていきます。

　AR過程では初期値の決定が問題になります。 $y$ が分布が定まっている場合は、その分布に従う確率変数とし、 $y$ の分布が定まっていない場合は、ある定数とすることが一般的です。ただ、定常過程では、初期値の影響は時間とともになくなるので、大きな問題にはなりません。

　沖本先生の本のとおり、具体例を示します。次のようなAR(1)を考えます。

$y_{t} = 1 + 0.5 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N.(1)$

　初期値を $y_{0} = 0$ として、 $\epsilon_{1} = - 2.1$ 、 $\epsilon_{2} = 0.9$ とします。このときに順番に解いていくと、

$y_{1} = 1 + 0.5 y_{0} \epsilon_{1} = 1 + 0.5 \dot - 2.1 = - 1.1$
$y_{2} = 1 + 0.5 y_{1} \epsilon_{2} = 1 + 0.5 \dot (-1.1) + 0.9 = 1.35$

となります。

　ここで、MA(1) 過程と同様に、パラメータを調整したAR(1)過程をグラフ化してみます。RのコマンドについてはMA(1)で多少説明していますので、必要な部分のみ解説します。

　以下の６つのパラメータのAR(1)過程をグラフ化します。

$\displaystyle \begin{align} &(a) y_{t} = 2 + 0.8 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \\ &(b) y_{t} = - 2 + 0.3 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{0.5}) \\ &(c) y_{t} = 0 - 0.3 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{2}) \\ &(d) y_{t} = - 2 - 0.8 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \\ &(e) y_{t} = 2 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \\ &(f) y_{t} = 2 + 1.1 y_{t-1} + \epsilon_{t}, \epsilon_{t} \sim W.N. (\sigma^{1}) \end{align}$

　Rのコードです。arima.sim関数の第2引数ですが、今回はAR過程なので、order=c(1, 0, 0) となっています。ちなみに、MA過程の時はorder=c(0, 0, 1) でした。次にAR過程の係数はar=c(0.8)として指定します。

set.seed(1)
d1 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(0.8)), sd=1) +2
plot(d1, main="(a) c=2, Φ=0.8, σ=1")
set.seed(1)
d2 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(0.3)), sd=0.5) - 2
plot(d2, main="(b) c=-2, Φ=0.3, σ=0.5")
set.seed(1)
d3 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(-0.3)), sd=2)
plot(d3, main="(c) c=0, Φ-=0.3, σ=2")
set.seed(1)
d4 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(-0.8)), sd=1) - 2
plot(d4, main="(d) c=-2, Φ=-0.8, σ=1")
set.seed(1)
d5 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(1)), sd=1) + 2
plot(d5, main="(e) c=2, Φ=1, σ=1")
set.seed(1)
d6 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(1.1)), sd=1) + 2
plot(d6, main="(f) c=2, Φ=1.1, σ=1")

　グラフ化するとこのようになります。

　まず、(e)と(f)ですが、係数が $1$ 以上の場合、AR過程は定常ではなくなるため、エラーで出力されません。そのため、グラフも出力されません。

　AR過程はMA過程と違い、パラメータの値によって定常か非定常かが変わります。先ほど、係数が $1$ 以上と書きましたが、正確には係数 $|\phi_{1}| \ge 1$ の時に非定常となります。
　ここではAR(1)過程が定常という場合のみを解説していきます。

　AR(1)過程の期待値についてです。グラフでは概ね期待値と定数項が一致しているようにも見えますが、実際は期待値と定数項は一致しません。(a)あたりでは定数項は $2$ ですが、やや下の方に過程の期待値があるようです。AR(1)過程の期待値をみています。

$\displaystyle y_{t} = c + \phi_{1} y_{t-1} + \epsilon_{t}$
$\displaystyle E(y_{t}) = E(c + \phi_{1} y_{t-1} + \epsilon_{t}) = c + \phi_{1} E(y_{t-1})$

ここで、ホワイトノイズの性質により $E(\epsilon_{t}) = 0$ です。また、 $y$ は定常過程ですから、

$\displaystyle E(y_{t}) = E(y_{t-1}) = \mu$

です。よって、AR(1)過程の期待値

$\displaystyle \mu = c + \phi_{1} \mu$

$\displaystyle \mu = \frac{c}{(1 - \phi_{1})}$

が得られます。

　次に分散ですが、MA過程と同様に過程の分散は撹乱項の分散 $\sigma^{2}$ と異なります。MA(1)過程で使用した(a)のグラフと今回の(a)のグラフを比べてみます。ちなみに左側はMA(1)過程、右側がAR(1)過程です。

　撹乱項の分散は $1$ ですが、過程の変動幅は違います。AR(1)過程の方がやや大きいです。つまり、AR(1)過程の分散は撹乱項のそれより大きくなります。では分散を見ていきます。

$\displaystyle y_{t} = c + \phi_{1} y_{t-1} + \epsilon_{t}$
$\displaystyle \begin{align} Var(y_{t}) &= Var(c + \phi_{1} y_{t-1} + \epsilon_{t}) \\ &= \phi_{1}^{2} y_{t-1}^{2} + 2 \phi_{1} y_{t-1} (c + \epsilon_{t} - \mu) + (c + \epsilon_{t} - \mu)^{2} \\ &= \phi_{1}^{2} Var(y_{t-1}) + 2 \phi_{1} Cov(y_{t-1} , \epsilon_{t}) + Var(\epsilon_{t}) \\ &= \phi_{1}^{2} Var(y_{t-1}) + \sigma^{2} \end{align}$

　ここで、 $y_{t}$ が定常であれば、 $\displaystyle \gamma_{0} = Var(y_{t}) = Var(y_{t-1})$ ですので、AR(1)過程の分散は、
$\displaystyle \gamma_{0} = \frac{\sigma^{2}}{(1 - \phi_{1}^{2})}$
となります。

　最後に自己相関をみてみましょう。先ほどのMA(1)過程のグラフ(a)とAR(1)過程のグラフを比べると、係数( $\theta_{1}$ と $\phi_{1}$ はともに $0.8$ )は同じでもAR(1)過程のグラフがより滑らかであることがわかります。

　これは、この場合のように係数が同じとき、１次自己相関の絶対値はAR(1)過程の方が大きくなることと、AR(1)過程は２次以降の相関係数もすべて正になることの２点によるからです。

　AR過程のk次の自己共分散はこのようになります。

$\displaystyle \begin{align} \gamma_{k} &= Cov(y_{t} , y_{t-k}) \\ &= Cov(\phi_{1} y_{t-1} + \epsilon_{t} , y_{t-k}) \\ &= Cov(\phi_{1} y_{t-k} , y_{t-k}) + Cov(\epsilon_{t} , y_{t-k}) \\ &= \phi_{1} \gamma_{k-1} \end{align}$

　ここで、両辺を $\gamma_{0}$ で割ると、

$\displaystyle \gamma_{k} = \phi_{1} \rho_{k-1}$

となり、これはユール・ウォーカー方程式と呼ばれます。ユール・ウォーカー方程式は、AR過程の自己相関が、 $y_{t}$ が従うAR過程と同一の係数 $\phi_{1}$ をもつ差分方程式に従うということを表しています。

　ユール・ウォーカー方程式と $\rho_{0} = 1$ を用いて、AR過程の自己相関は順番に求めることができます。また、ユール・ウォーカー方程式により、AR(1)過程の一般的な自己相関は、

$\displaystyle \rho_{k} = \phi_{1}^{2} \rho_{k-2} = \phi_{1}^{3} \rho_{k-3} = \cdots = \phi_{1}^{k} \rho_{0} = \phi_{1}^{k}$

となります。また、[ $| \phi_{1} | < 1$ より、AR(1)過程の自己相関の絶対値は指数的に減衰していきます。コレログラムを描いてみます。

set.seed(1)
d1 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(0.8)), sd=1)
acf(d1, main="(a)  θ1=0.8, θ2=0")
set.seed(1)
d2 <- arima.sim(n=100, model=list(order=c(1, 0, 0), ar=c(-0.8)), sd=1)
acf(d2, main="(a)  θ1=0.8, θ2=0")

　コレログラムを描くacf()関数とarima.sim()関数は以前に説明していますのでここでは割愛します。AR(1)過程の自己相関の絶対値は指数的に減衰していくことが、(a)と(b)からよくわかると思います。AR(1)過程のコレログラムは $\phi_{1}$ の符号によって、減衰していくことが確認できます。

　AR(1)過程を一般化したp次のAR過程AR(p)はこのように記述されます。

$\displaystyle y_{t} = c + \phi_{1} y_{t-1} + \phi_{2} y_{t-2} + \cdots + \phi_{p} y_{t-p} + \epsilon_{t}$

　AR(p)過程はy_{t} を定数と自身のp期過去の値に回帰したモデルになります。AR(1)同様にAR(q)もパラメータにより定常となるか非定常となるか変わってきます。定常AR(p)過程の性質は、

$\displaystyle \mu = E(y_{t}) = \frac{c}{1 - \phi_{1} - \phi_{2} - \cdots - \phi_{p}}$
$\displaystyle \gamma_{0} = Var(y_{t}) = \frac{\sigma^{2}}{1 - \phi_{1} \rho_{1} - \phi_{2} \rho_{2} - \cdots - \phi_{p} \rho_{p}}$
$\displaystyle \gamma_{k} = \phi_{1} \gamma_{k-1} + \phi_{2} \gamma_{k-2} + \cdots + \phi_{p} \gamma_{k-p}, k \ge 1$
$\displaystyle \rho_{k} = \phi_{1} \rho_{k-1} + \phi_{2} \rho_{k-2} + \cdots + \phi_{p} \rho_{k-p}, k \ge 1$

であり、特に $\gamma_{k}$ はユール・ウォーカー方程式と呼ばれます。そして、自己相関が指数的に減衰していきます。

　ここで、AR(2)過程のユール・ウォーカー方程式は、このようになります。

$\displaystyle \rho_{k} = \phi_{1} \rho_{k-1} + \phi_{2} \rho_{k-2}$

これに、 $k=1$ を代入し、 $\rho_{0} = 1$ と $\rho_{1} = \rho_{-1}$ から

$\displaystyle \rho_{1} = \phi_{1} \rho_{1} + \phi_{2} \rho_{-1} = \phi_{1} + \phi_{2} \rho_{1}$

が得られます。これを変形すると、

$\displaystyle \rho_{1} = \frac{\phi_{1}}{(1 - \phi_{2})}$

となります。次に、 $k=2$ を代入すると、

$\displaystyle \rho_{2} = \phi_{1} \rho_{1} + \phi_{2} \rho_{0} = \frac{\phi_{1}^{2}}{1 - \phi_{2}} + \phi_{2} = \frac{\phi_{1}^{2} + \phi_{2} - \phi_{2}^{2}}{1 - \phi_{2}}$

となり、 $\rho_{2}$ を求めることができます。AR(2)過程の場合は、特に $\rho_{k}$ の一般的な

$\displaystyle 1 - \phi_{1} z - \phi_{2} z^{2} = 0$

という２次方程式が異なる２つの解をもつとき、その解の逆数を $\lambda_{1}$ 、 $\lambda_{2}$ とすれば

$\displaystyle \rho_{k} = \frac{(1 - \lambda_{2}^{2}) \lambda_{1}^{k+1} - (1 - \lambda_{1}^{2}) \lambda_{2}^{k+1}}{(\lambda_{1} - \lambda_{2})(1 + \lambda_{1} \lambda_{2})}$

で自己相関係数が得られます。

　AR(p)の場合、 $\rho_{1}$ から $\rho_{k-1}$ までは、ユール・ウォーカー方程式と $\rho_{0} = 1$ 、 $\rho_{k} = \rho_{-k}$ から $p-1$ 次の連立方程式を解いて求めることができます。

　AR(2)過程のコレログラムを図示します。

set.seed(1)
d3 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.5, 0.35)), sd=1)
acf(d3, main="(a)  θ1=0.5, θ2=0.35")
set.seed(1)
d4 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.1, 0.5)), sd=1)
acf(d4, main="(a)  θ1=0.1, θ2=0.5")
set.seed(1)
d5 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.5, -0.8)), sd=1)
acf(d5, main="(a)  θ1=0.5, θ2=-0.8")
set.seed(1)
d6 <- arima.sim(n=100, model=list(order=c(2, 0, 0), ar=c(0.9, -0.8)), sd=1)
acf(d6, main="(a)  θ1=0.9, θ2=-0.8")