自己回帰移動平均モデルとは: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             k>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 の方程式とよぶ。
    • 移動平均項のない自己回帰モデルの場合は、k>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の合計 式(1)

また、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) を推定する場合、どのような結果が得られるかを予測できます。

参考


トップ   差分 バックアップ リロード   一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-10-18 (月) 11:21:00 (3312d)