時系列データ解析 05 ～移動平均過程

www.asakura.co.jp

　沖元先生の本の続きです。ここでは、移動平均モデルについての話からです。

　まず、1次MA過程（以下、MA(1)過程）はこのようにモデル化されます。

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

　移動平均過程は、ホワイトノイズの線形和になります。

　ホワイトノイズはすべての時点において期待値が $0$ で、分散が一定の確率過程です。また、自己相関をもたないので、弱定常過程となります。

　 $y_{t}$ が、移動平均過程に従うことを、 $y_{t} \sim MA(1)$ とあらわします。

　もっとも基本的な弱定常過程は

$\displaystyle y_{t} = \mu + \epsilon_{t}, \epsilon \sim W.N. (\sigma^{2} )$

ですから、移動平均過程はこの過程に、 $\theta_{1} \epsilon_{t - 1}$ という項が追加されたものになります。

　また、 $t$ 期のMA(1)過程と $t - 1$ 期のMA(1)を比較してみます。

$\displaystyle y_{t} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t - 1}$
$\displaystyle y_{t - 1} = \mu + \epsilon_{t - 1} + \theta_{1} \epsilon_{t - 2}$

　この二つのモデルは、 $\epsilon_{t - 1}$ という項を共通にもちます。そのため、 $y_{t}$ と $y_{t - 1}$ に相関が生じます。

　そして、MA(1)過程は、弱定常過程と比べて、 $\theta_{1}$ というパラメータが多くなっています。このパラメータが、自己相関の強さを決定づけます。

　MA(1)過程の確率変動は、ホワイトノイズ $\epsilon_{t}$ によって決まります。まず、 $\epsilon$ が決まり、続いて $y$ が決まっていきます。沖本先生の本では、具体例として、

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

というMA(1)過程で、 $\epsilon_{0} = 0.5$ 、 $\epsilon_{1} = - 2.1$ 、 $\epsilon_{2} = 0.9$ とすれば、

$\displaystyle y_{1} = 1 + \epsilon_{1} + 0.5 \epsilon_{0} = 1 - 2.1 + 0.5 \cdot 0.5 = - 0.85$

$\displaystyle y_{2} = 1 + \epsilon_{2} + 0.5 \epsilon_{1} = 1 + 0.9 + 0.5 \cdot ( - 2.1 ) = - 0.85$

という順番で $y$ が決められることが示されています。

　ここで、 $\epsilon$ を正規ホワイトノイズとして、それぞれのパラメータの組み合わせたMA(1)をグラフ化してみます。

　Rを使って、シミュレーションをしてみます。Rにはさまざまなツールが用意されていまして、ここではARIMAを例にMA(1)過程をみていきます。ちなみに、ARIMA過程はこのようなコードで記述します。ARIMA過程そのものについては次回以降に説明します。

d <- arima.sim(n=100, model=list(order=c(2, 0, 4), ar=c(0.1, 5), ma=c(-0.2, 0.3, 0.4, -0.5)), sd=1)

　arima.simの引数について説明します。n=100がシミュレーションの数です。今回の話の流れでは、時系列のデータを $T = 100$ まで実験してみるということです。

　第2引数は、model=list()となっていますが、中身をみていきます。order=()の引数は、自己回帰過程、和分過程、移動平均過程の次数になります。自己回帰過程、和分過程は後述します。ar=()の引数は、係数で、ma=()も同様です。

　第3引数については、分散の大きさといいますか、上の例では分散も標準偏差（標準誤差）も $1$ です。

　それでは、シミュレーションしてみます。沖本先生の例にならって６パターンにします。

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

　ついでにコードも記述します。

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

　ar=()については省略しています。plot()はグラフに表示するコマンドです。main=""はグラフのタイトルのようなものです。set.seet(1)は、乱数種を指定して乱数を発生させる関数です。すべて同じ乱数を発生させてシミュレーションしています。

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

　グラフから確認できることは、系列が $\mu$ の周りを変動していることです。つまり、MA(1)過程の期待値は $\mu$ になります。

$\begin{align} E(y_{t}) &= E(\mu + \epsilon_{t} + \theta_{1} \epsilon_{t - 1}) \\ &= E(\mu) + E(\epsilon_{t}) + E(\theta_{1} \epsilon_{t - 1}) \\ &= E(\mu) + E(\epsilon_{t}) + \theta_{1} E(\epsilon_{t - 1}) \\ &= \mu \because E(\epsilon_{t}) = 0, E(\epsilon_{t - 1}) = 0)\end{align}$

　次にMA(1)過程の分散は、撹乱項 $\epsilon_{t}$ の分散 $\sigma^{2}$ に等しいと既に説明しました。しかし、グラフを見ると分かるとおり、撹乱項の分散が同じでも若干変動幅が異なることが分かると思います。

　つまり、MA(1)過程分散は、撹乱項の分散よりも大きくなるのです。

$\displaystyle y_{t} - \mu = \epsilon_{t} - \theta_{1} \epsilon_{t - 1}$

$\begin{align} (y_{t} - \mu)^{2} &= (\epsilon_{t} - \theta_{1} \epsilon_{t - 1})^{2} \\ &= \epsilon_{t}^{2} - 2 \theta_{1} \epsilon_{t} \epsilon_{t - 1} + \theta_{1}^{2} \epsilon_{t - 1}^{2} \end{align}$

$\begin{align} \gamma_{0} &= E \left[ \epsilon_{t}^{2} \right] - 2 E \left[ \theta_{1} \epsilon_{t} \epsilon_{t - 1} \right] + \theta_{1}^{2} E \left[ \epsilon_{t - 1}^{2} \right] \\ &= (1 + \theta_{1}^{2} ) \sigma^{2} \end{align}$

　したがって、MA(1)過程の分散は $\theta_{1}^{2} \sigma^{2}$ だけ、撹乱項の分散よりも大きくなります。

　また、グラフの滑らかさは $\theta_{1}$ が大きくなるにしたがって、より滑らかになります。(a)が一番滑らかであり、(f)が一番ギザギザしているます。これは、MA(1)過程の自己相関が $0$ でないためです。MA(1)過程をもう一度見てみましょう。

$\displaystyle y_{t} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t - 1}, \epsilon$
$\displaystyle y_{t - 1} = \mu + \epsilon_{t - 1} + \theta_{1} \epsilon_{t - 2}$

　 $\theta_{1}$ が正の場合、 $y_{t}$ の $\theta_{1} \epsilon_{t-1}$ と $y_{t-1}$ の $\epsilon_{t-1}$ が同じ正の符号になるため、 $y_{t}$ と $y_{t-1}$ が同じ方向に動く傾向をもちます。

　その結果、１次の正の自己相関が生じて、その自己相関は $\theta_{1}$ の値が $1$ に近づくほど強くなります。

　逆に、 $\theta_{1}$ が負の場合は、 $y_{t}$ と $y_{t-1}$ は逆に動く傾向があり、１次の負の自己相関が生まれます。そのため、グラフはよりギザギザになります。

　次に、MA(1)過程の自己共分散を求めます。

$\displaystyle \begin{align} \gamma_{1} &= Cov( y_{t}, y_{t-1}) \\ &= Cov(\mu + \epsilon_{t} + \theta_{1} \epsilon_{t-1}, \mu + \epsilon_{t-1} + \theta_{1} \epsilon_{t-2} \\ &= Cov(\epsilon_{t}, \epsilon_{t-1}) + Cov(\epsilon_{t}, \theta_{1} \epsilon_{t-2} ) + Cov(\theta_{1} \epsilon_{t-1} , \epsilon_{t-1}) + Cov(\theta_{1} \epsilon_{t-1}, \theta_{t-1} \epsilon_{t-2}) \\ &= \theta_{1} Cov(\epsilon_{t-1}, \epsilon_{t-1}) \\ &= \theta_{1} \sigma_{2} \end{align}$

　これは、ホワイトノイズ $\epsilon$ の $k$ 次の自己共分散が $0$ であり、分散が $0$ であることから導けます。

　よって、MA(1)過程の１次自己相関を求めると、

$\displaystyle \rho_{1} = \frac{\gamma_{1}}{\gamma_{0}} = \frac{\theta_{1}}{1 + \theta_{1}^{2}}$

となります。MA(1)過程の自己相関が $\mp 1$ のとき、 $\frac{1}{2}$ になります。したがいまして、１次の自己相関の絶対値が $\frac{1}{2}$ より大きいMA(1)過程はモデル化できません。

　MA(1)過程の２次以降の自己共分散は、

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

となります。つまり、MA(1)過程の２次以降の自己相関は $0$ ということです。そのため、MA(1)過程は１次の自己相関をモデル化できますが、２次以降の自己相関を記述することはできません。

　以上のとおり、期待値も分散も $t$ によらず、定常性の性質を満たすことがわかりました。

　最後に、MA(1)過程のコレログラムを確認してみます。まずRのコードを記載しておきます。

set.seed(1)
d1 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.2, 0)), sd=1)
acf(d1, main="(a)  θ1=0.8, θ2=0")
set.seed(1)
d2 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(-0.8, 0)), sd=1)
acf(d2, main="(b) θ1=-0.8, θ2=0")
set.seed(1)
d3 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.8, 0.3)), sd=1)
acf(d3, main="(c) θ1=0.8, θ2=0.3")
set.seed(1)
d4 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.3, 0.8)), sd=1)
acf(d4, main="(d) θ1=0.3, θ2=0.8")
set.seed(1)
d5 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(0.8, -0.3)), sd=1)
acf(d5, main="(e) θ1=0.8, θ2=-0.3")
set.seed(1)
d6 <- arima.sim(n=100, model=list(order=c(0, 0, 2), ma=c(-0.8, -0.3)), sd=1)
acf(d6, main="(f)  θ1=-0.8, θ2=-0.3")

　acf()は自己相関関数、つまりコレログラムを表しています。先ほどのarima.sim()の中身ですが、今回はMA(1)過程のオーダーを2にしていますので、ma=c()中に指定する係数も2つになっています。

　まず、それぞれの係数 [tex \theta] を見るとコレログラムが係数に応じて変化しているのがわかると思います。(a)、(b)については、1次のみ有意な自己相関があり、2次以降は有意な自己相関はありません。

　一般化した $q$ 次移動平均過程は、

$y_{t} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t-1} + \theta_{2} \epsilon_{t-2} + \cdots + \theta_{q} \epsilon_{t-q}, \epsilon_{t} \sim W.N.(\sigma^{2})$

と表現でき、MA(q)過程と書きます。MA(q)過程は、現在とq期間の過去のホワイトノイズの線形和に定数 $\mu$ を加えたものになります。

　MA(q)過程の性質は、

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

$\displaystyle \gamma_{0} = Var(y_{t}) = (1 + \theta_{1}^{2} + \theta_{2}^{2} + \cdots + \theta_{q}^{2}) \sigma^{2}$

$\displaystyle \gamma_{k} = \begin{cases}{} (\theta_{k} + \theta_{1} \theta_{k+1} + \cdots + \theta_{q-k} \theta_{q}) \sigma^{2}, &1 \leq k \leq q \\ 0, & k \geq q + 1 \end{cases}$

$\displaystyle \rho_{k} = \begin{cases}{} \frac{\theta_{k} + \theta_{1} \theta_{k+1} + \cdots + \theta_{q-k} \theta_{q}}{1 + \theta_{1}^{2} + \theta_{2}^{2} + \cdots + \theta_{q}^{2}}, & 1 \leq k \leq q \\ 0, & k \geq q + 1 \end{cases}$

であり、重要なことはMA過程は定常過程であり、MA(q)過程のq+1次以降の自己相関は $0$ ということです。したがいまして、q次の自己相関構造をモデル化するには、q個のパラメータが必要になります。長期間の自己相関のモデル化はあまりに多くのパラメータが必要になるため、望ましいことではありません。

　次回は自己相関過程の話に移ります。

マツリさんの日記

androidプログラミング初心者の奮闘日記です。たまに統計学もしてます。

時系列データ解析 05 ～移動平均過程