自己回帰移動平均モデル

自己回帰移動平均モデルとは：Autoregressive moving average model、ARMAモデル †

統計学において時系列データに適用されるモデルである。George Box と G. M. Jenkins の名をとって "ボックス・ジェンキンスモデル" とも呼ばれる。

モデルは自己回帰(AR)部分と移動平均(MA)部分からなる。一般に ARMA(p,q)モデルと表記され、p は自己回帰部分の次数、q は移動平均部分の次数を表す。

ARMA(p, q)の定義 †

ARMA(p, q)という表記は、p次の自己回帰とq次の移動平均を組合わせたモデルを指す。

yt=a1yt-1+a2yt-2+...+apyt-p + b1ut-1+b2ut-2+...+bqut-q

ただし、誤差項 ut は一般に「独立かつ同一の分布に従う」(i.i.d.)無作為変数であり、ゼロを平均値とする正規分布に従う。すなわち ut ~ N(0,σ2) で、σ2 は分散である。

実データに適用する場合、ARMAモデルの p と q を選択後、誤差項を最小化するパラメータを探るため最小二乗法を使うのが普通である。また、実データに適合する最小の p および q を見つけることでよい結果が得られることが知られている。
純粋なARモデルでは、これに Yule-Walker 方程式を利用することができる

↑

自己回帰モデルと最小二乗法 †

AR(p)モデルは次の式で表される。

yt=c+a1yt-1+a2yt-2+...+apqut-p+ut
C:定数項
誤差項utは、σ2の分散をもつホワイトノイズ

モデルとして定常的であるために、パラメータの値には何らかの制約が必要である。例えば、|a1| > 1 となる AR(1)モデルは定常的ではない。
時系列データが系列相関をもつかか否かの検定方法として、ダービンワトソン比は有名である。

↑

AR(1)の場合 †

yt=C+ayt-1+ut

この過程は |a|<1 であれば、共分散定常性を有する。 a=1であれば、ランダムウォークと見なされ、共分散定常性を有しない

期待値は、両辺の期待値をとると
```
E(yt)=E(C+ayt-1+ut)
     =C+aE(yt-1)+E(ut)
     =C+aμy
```
従って
```
μy=E(yt)=c/(1-a)
```
となる。c = 0 なら、平均も 0 になる。
c=0ならば、分散は次のようになる。
```
Var(yt)=E(yt^2)-μy^2=σ2/(1-a^2)
```
c=0ならば、自己共分散は次の式で表される。
```
V(yt+iyt)=E(yt+iyi)-μy^2=[σ2/(1-a^2)]a^|i|
```
すなわち、|a|<1なので、i個ずれれば、自己相関が減衰して消えいく。

↑

パラメータ推定：予測誤差二乗和の最小化：The Yule Walker Equations for the AR Coefficients の方程式 †

観測値の自己相関係数は、信号の特徴抽出にも用いられます。この自己相関係数を用いる代わりに自己回帰モデルを作って、分析したり予測することも考えられます。パラメータの推定方法に、観測値の自己相関係数を用いたYule Walker Equations(ユール・ウオーカーの方程式があります。

自己相関係数の定義
```
Rk=E(yt-kyt)
```
ARMA(p,q)の時系列信号yt の相関関数Rk は、AR パラメータのみに関係する！！。まずARMA(p.q)の式両辺にyt-k をかけると
```
yt-kyt=a1yt-kyt-1+..+apyt-kyt-p + b1yt-kut-1+..+bqyt-kut-q
```
さらに、両辺の時間平均をとると
```
Rk+a1Rk-1+a2Rk-2+...+ apRk-p=b1E(yt-kut-1)+...+bqE(yt-kqut-q)
```
一方、システムのインパルス応答をhjとすれば
```
yt=Σhjut-j　ΣはJ=0~無限大の和
```
ut は平均値0、分散σ の白色雑音系列であるので、上の期待値は
```
Rk+a1Rk-1+a2Rk-2+...+ apRk-p=b1E(yt-kut-1)+...+bqE(yt-kqut-q)　k<qまたはK=q
Rk+a1Rk-1+a2Rk-2+...+ apRk-p=0　　　　　　　　　　　　　ｋ>q+1またはK=q+1
```
kが q+1 以上のときにARMA(p, q) モデルの時系列信号yt の相関関数Rk がARパラメータにのみ関係していることを表している。 kが q+1 以上のときに
```
Rk+a1Rk-1+a2Rk-2+...+ apRk-p=0
```
であるので、この方程式を解けば、自己回帰項のパラメータが求められることが判る。この方程式をYule-Walker の方程式とよぶ。
- 移動平均項のない自己回帰モデルの場合は、ｋ>q+1またはK=q+1の制約なしに、パラメータ推定が出来ることになる。

↑

自己回帰モデルの係数を求める正規方程式（Yule-Walker 方程式）の導出 †

ytを過去のP 個の観測値yt-1,yt-2,....yt-pの線形結合で表現してみよう。

y*t=a1yt-1+a2yt-2+...+apyt-p

このとき、t期までの観測データの予測誤差は、

et=yt-y*t　t=0,1,2,...t

で与えられるが、この予測誤差の二乗期待値σ^2は　J=E( et^2) であたえられる。

J=E｛[yt-a1yt-1-a2yt-2-...-apyt-p]^2｝
=ΣΣ　akaj　E(yt-lyt-j) 前のΣはk=0~p,後ろのΣはj=0~pの和
=Σak　Σaj　Rj-k

ノイズ関数の値は互いに独立であり、ゼロより大きいj についてyt-jは誤差項utに独立であるという性質を用いている。

この予測誤差の二乗和を最小にする係数を求める式がYule-Walker 方程式であることを示そう。なお、係数は

a = (1,a1,a2,...ap)

とする（実現値yt に乗じるa0 は1 と考える）。つまり、式をa0 = 1 を除くakで偏微分して0 とおけば、次のp個の線形式が得られる。

ΣRj-kaj=0　k=1~p　Σはj=1~pの合計　式（１）

また、a0 についてはそのまま書き下すと、

σ*^2=ao(a1R1+a2R2+...+apRp)
    =a1R1+a2R2+...+apRp=ΣRj　　式(2)

が成立する。なお、一般に、式(1) と式(2) をあわせてYule-Walker 方程式と呼ばれている。

↑

Yule-Walker の方程式の逐次解法 †

逐次解法が存在する。 Levinson アルゴリズムとも呼ばれる。詳細は、下記を参考のこと。

↑

自己相関係数のパワースペクトラム（フーリエ変換）による計算 †

自己相関関数Rk とパワースペクトルE(n) = ｜X(n)|^2 は、フーリエ変換によって結ばれている。この関係はウィーナーヒンチンの定理とよばれ、

E(n) = F[Rk]　：F(・)はフーリエ変換
Rk = F-1[E(n)]：F-1(・)は逆フーリエ変換

と表される。この関係を用いると、次のように自己相関関数が計算できる。

1. 測定データ対してFFT を実施し、測定データのスペクトルX(n) (n = 0,1,...) を求める。
2. 測定データのパワースペクトル｜X(n)|^2 を求める。すなわち
```
Re[ ｜X(n)|^2 ] = Re[｜X(n)|]^2  + Im[｜X(n)|] ^2
Im[ ｜X(n)|^2 ] = 0
```
を計算する。
3. 測定データのパワースペクトル｜X(n)|^2 (k = 0,1,...) に対して逆FFT を実施し、測定データの自己相関関数Rk (k = 0,1,...) を求める。
4. 以下のように自己相関関数の規格化を行う。
```
Rk = Rk/R0 (k = 0,1,...)
```
- PSD - MATLABでの説明：Yule-walker 法では、指定した次数の AR 予測モデルを信号に近似することで、スペクトル密度が計算されます.まず、指定した次数の AR (全極) モデルから信号を生成します。関数 freqz を使用して、AR フィルターの周波数応答の振幅をチェックできます。これにより、関数 pyulear を使用してパワースペクトラム密度(PSD) を推定する場合、どのような結果が得られるかを予測できます。

↑

人気の30件

最新の20件