自己回帰移動平均モデル

自己回帰移動平均モデル（じこかいきいどうへいきんモデル、英: autoregressive moving average model、ARMAモデル）は自己回帰モデルによる線形フィードバックと移動平均モデルによる線形フィードフォワードによりシステムを表現するモデルである。George Box と G. M. Jenkins の名をとって "ボックス・ジェンキンスモデル" とも呼ばれる。 ARMAモデルは時系列データの将来値を予測するツールとして機能する。 p {\displaystyle p} 次の自己回帰 (AR) および q {\displaystyle q} 次の移動平均 (MA) からなる自己回帰移動平均モデル ARMA(p, q) {\displaystyle {\text{ARMA(p, q)}}} は以下のように定義される。 X t = c + ∑ i = 1 p φ i X t − i + ∑ i = 0 q θ i ε t − i {\displaystyle X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\sum _{i=0}^{q}\theta _{i}\varepsilon _{t-i}} ここで c {\displaystyle c} は定数、 φ k {\displaystyle \varphi _{k}} は自己回帰パラメータ、 θ k {\displaystyle \theta _{k}} は移動平均パラメータ ( θ 0 = 1 {\displaystyle \theta _{0}=1} )、 ε t {\displaystyle \varepsilon _{t}} は時刻 t {\displaystyle t} におけるホワイトノイズである。 すなわちARMAモデルでは、各時刻でサンプリングされたホワイトノイズが過去時刻 q {\displaystyle q} まで重み付け和でフィードフォワードされ、また過去時刻 p {\displaystyle p} まで出力が線形フィードバックされ、定数に足しこまれることで現在値が得られる。 AR(p) という表記は次数 p の自己回帰モデルを表す。AR(p)モデルは次の式で表される。 X t = c + ∑ i = 1 p φ i X t − i + ε t . {\displaystyle X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}.\,} ここで φ 1 , … , φ p {\displaystyle \varphi _{1},\ldots ,\varphi _{p}} はモデルのパラメータ、 c {\displaystyle c} は定数項、 ε t {\displaystyle \varepsilon _{t}} は誤差項（後述）である。定数項は単純化するために省かれることが多い。 自己回帰モデルは基本的に無限インパルス応答フィルタに一種の変形を加えたものである。 モデルとして定常的であるために、パラメータの値には何らかの制約が必要である。例えば、|φ1| > 1 となる AR(1)モデルは定常的ではない。 AR(1)過程は次の式で表される。 X t = c + φ X t − 1 + ε t , {\displaystyle X_{t}=c+\varphi X_{t-1}+\varepsilon _{t},\,} ここで、 ε t {\displaystyle \varepsilon _{t}} は、 σ 2 {\displaystyle \sigma ^{2}} の分散に従うホワイトノイズである（ φ 1 {\displaystyle \varphi _{1}} のような添え字は省いてある）。この過程は | φ | < 1 {\displaystyle |\varphi |<1} であれば、共分散定常性を有する。 φ = 1 {\displaystyle \varphi =1} であれば、 X t {\displaystyle X_{t}} は単位根を表し、ランダムウォークと見なされ、共分散定常性を有しない。そうでない場合、 X t {\displaystyle X_{t}} の期待値の計算は単純である。ここで共分散定常性を以下のように定式化する。 E ( X t ) = E ( c ) + φ E ( X t − 1 ) + E ( ε t ) ⇒ μ = c + φ μ + 0. {\displaystyle {\mbox{E}}(X_{t})={\mbox{E}}(c)+\varphi {\mbox{E}}(X_{t-1})+{\mbox{E}}(\varepsilon _{t})\Rightarrow \mu =c+\varphi \mu +0.} 従って、次のようになる。 μ = c 1 − φ , {\displaystyle \mu ={\frac {c}{1-\varphi }},} ここで μ {\displaystyle \mu } は平均である。c = 0 なら、平均も 0 になり、分散は次のようになる。 var ( X t ) = E ( X t 2 ) − μ 2 = σ 2 1 − φ 2 . {\displaystyle {\textrm {var}}(X_{t})=E(X_{t}^{2})-\mu ^{2}={\frac {\sigma ^{2}}{1-\varphi ^{2}}}.} 自己共分散は次の式で表される。 B n = E ( X t + n X t ) − μ 2 = σ 2 1 − φ 2 φ | n | . {\displaystyle B_{n}=E(X_{t+n}X_{t})-\mu ^{2}={\frac {\sigma ^{2}}{1-\varphi ^{2}}}\,\,\varphi ^{|n|}.} この自己共分散関数は減衰時間 τ = − 1 / ln ⁡ ( φ ) {\displaystyle \tau =-1/\ln(\varphi )} で減衰する（これを確かめるには、 B n = K ϕ | n | {\displaystyle B_{n}=K\phi ^{|n|}} で K {\displaystyle K} が n {\displaystyle n} に独立な場合を考えればよい。 ϕ | n | = e | n | ln ⁡ ϕ {\displaystyle \phi ^{|n|}=e^{|n|\ln \phi }} であり、指数関数的減衰の法則 e − n / τ {\displaystyle e^{-n/\tau }} に適合することに注意されたい）。スペクトル密度関数は自己共分散関数の逆フーリエ変換である。離散系では、離散時間逆フーリエ変換が適用される。 Φ ( ω ) = 1 2 π ∑ n = − ∞ ∞ B n e − i ω n = 1 2 π ( σ 2 1 + φ 2 − 2 φ cos ⁡ ( ω ) ) . {\displaystyle \Phi (\omega )={\frac {1}{\sqrt {2\pi }}}\,\sum _{n=-\infty }^{\infty }B_{n}e^{-i\omega n}={\frac {1}{\sqrt {2\pi }}}\,\left({\frac {\sigma ^{2}}{1+\varphi ^{2}-2\varphi \cos(\omega )}}\right).} X j {\displaystyle X_{j}} が離散的であるため、この式の分母にあるコサインの項が折り返し雑音（エイリアス）を表している。標本化間隔( Δ t = 1 {\displaystyle \Delta t=1} )が減衰時間( τ {\displaystyle \tau } )より十分に小さいと仮定すると、 B n {\displaystyle B_{n}} に連続体近似を適用できる。 B ( t ) ≈ σ 2 1 − φ 2 φ | t | {\displaystyle B(t)\approx {\frac {\sigma ^{2}}{1-\varphi ^{2}}}\,\,\varphi ^{|t|}} この場合、スペクトル密度はローレンツ分布に従う。 Φ ( ω ) == 1 2 π σ 2 1 − φ 2 γ π ( γ 2 + ω 2 ) {\displaystyle \Phi (\omega )=={\frac {1}{\sqrt {2\pi }}}\,{\frac {\sigma ^{2}}{1-\varphi ^{2}}}\,{\frac {\gamma }{\pi (\gamma ^{2}+\omega ^{2})}}} ここで γ = 1 / τ {\displaystyle \gamma =1/\tau } は減衰時間 τ {\displaystyle \tau } に関する角周波数である。 X t {\displaystyle X_{t}} の別の表現方法として、最初の式で X t − 1 {\displaystyle X_{t-1}} を c + φ X t − 2 + ε t − 1 {\displaystyle c+\varphi X_{t-2}+\varepsilon _{t-1}} に置き換える方法がある。これを再帰的に N回繰り返すと次の式になる。 X t = c ∑ k = 0 N − 1 φ k + φ N X φ − N + ∑ k = 0 N − 1 φ k ε t − k . {\displaystyle X_{t}=c\sum _{k=0}^{N-1}\varphi ^{k}+\varphi ^{N}X_{\varphi -N}+\sum _{k=0}^{N-1}\varphi ^{k}\varepsilon _{t-k}.} N が無限大に近づくと、 φ N {\displaystyle \varphi ^{N}} はゼロに近づき、最終的に次の式が得られる。 X t = c 1 − φ + ∑ k = 0 ∞ φ k ε t − k {\displaystyle X_{t}={\frac {c}{1-\varphi }}+\sum _{k=0}^{\infty }\varphi ^{k}\varepsilon _{t-k}} AR(p)モデルは次の方程式で与えられる。 X t = ∑ i = 1 p φ i X t − i + ε t . {\displaystyle X_{t}=\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}.\,} これはパラメータ φ i {\displaystyle \varphi _{i}} （i = 1, ..., p）に基づいている。これらパラメータは以下の Yule-Walker方程式で計算できる可能性がある。 γ m = ∑ k = 1 p φ k γ m − k + σ ε 2 δ m {\displaystyle \gamma _{m}=\sum _{k=1}^{p}\varphi _{k}\gamma _{m-k}+\sigma _{\varepsilon }^{2}\delta _{m}} ここで m = 0, ... , p であり、p + 1 個の方程式となる。 γ m {\displaystyle \gamma _{m}} は X の自己共分散関数、 σ ε {\displaystyle \sigma _{\varepsilon }} は入力ノイズ過程の標準偏差、δm はクロネッカーのデルタである。 この式の最後の部分は m = 0 のときだけ 0 でない値となるので、この方程式は一般に m > 0 のときの行列式で表すことで解ける。 [ γ 1 γ 2 γ 3 ⋮ ] = [ γ 0 γ − 1 γ − 2 … γ 1 γ 0 γ − 1 … γ 2 γ 1 γ 0 … … … … … ] [ φ 1 φ 2 φ 3 ⋮ ] {\displaystyle {\begin{bmatrix}\gamma _{1}\\\gamma _{2}\\\gamma _{3}\\\vdots \\\end{bmatrix}}={\begin{bmatrix}\gamma _{0}&\gamma _{-1}&\gamma _{-2}&\dots \\\gamma _{1}&\gamma _{0}&\gamma _{-1}&\dots \\\gamma _{2}&\gamma _{1}&\gamma _{0}&\dots \\\dots &\dots &\dots &\dots \\\end{bmatrix}}{\begin{bmatrix}\varphi _{1}\\\varphi _{2}\\\varphi _{3}\\\vdots \\\end{bmatrix}}} これにより φ {\displaystyle \varphi } が全て求められる。また、m = 0 のときは次のようになる。 γ 0 = ∑ k = 1 p φ k γ − k + σ ε 2 {\displaystyle \gamma _{0}=\sum _{k=1}^{p}\varphi _{k}\gamma _{-k}+\sigma _{\varepsilon }^{2}} これにより σ ε 2 {\displaystyle \sigma _{\varepsilon }^{2}} が求められる。 AR過程を定義する方程式は次の通りである。 X t = ∑ i = 1 p φ i X t − i + ε t . {\displaystyle X_{t}=\sum _{i=1}^{p}\varphi _{i}\,X_{t-i}+\varepsilon _{t}.\,} 両辺に Xt-m をかけて、期待値を求めるとしたとき、次のようになる。 E [ X t X t − m ] = E [ ∑ i = 1 p φ i X t − i X t − m ] + E [ ε t X t − m ] . {\displaystyle E[X_{t}X_{t-m}]=E\left[\sum _{i=1}^{p}\varphi _{i}\,X_{t-i}X_{t-m}\right]+E[\varepsilon _{t}X_{t-m}].} 自己共分散関数の定義から、 E [ X t X t − m ] = γ m {\displaystyle E[X_{t}X_{t-m}]=\gamma _{m}} である。ノイズ関数の値は互いに独立であり、ゼロより大きい m について Xt − m は εt に独立である。m ≠ 0 の場合、 E [ ε t X t − m ] = 0 {\displaystyle E[\varepsilon _{t}X_{t-m}]=0} となる。m = 0 の場合、次のようになる。 E [ ε t X t ] = E [ ε t ( ∑ i = 1 p φ i X t − i + ε t ) ] = ∑ i = 1 p φ i E [ ε t X t − i ] + E [ ε t 2 ] = 0 + σ ε 2 , {\displaystyle E[\varepsilon _{t}X_{t}]=E\left[\varepsilon _{t}(\sum _{i=1}^{p}\varphi _{i}\,X_{t-i}+\varepsilon _{t})\right]=\sum _{i=1}^{p}\varphi _{i}\,E[\varepsilon _{t}\,X_{t-i}]+E[\varepsilon _{t}^{2}]=0+\sigma _{\varepsilon }^{2},} 従って、次が得られる。 γ m = E [ ∑ i = 1 p φ i X t − i X t − m ] + σ ε 2 δ m . {\displaystyle \gamma _{m}=E\left[\sum _{i=1}^{p}\varphi _{i}\,X_{t-i}X_{t-m}\right]+\sigma _{\varepsilon }^{2}\delta _{m}.} さらに E [ ∑ i = 1 p φ i X t − i X t − m ] = ∑ i = 1 p φ i E [ X t X t − m + i ] = ∑ i = 1 p φ i γ m − i , {\displaystyle E\left[\sum _{i=1}^{p}\varphi _{i}\,X_{t-i}X_{t-m}\right]=\sum _{i=1}^{p}\varphi _{i}\,E[X_{t}X_{t-m+i}]=\sum _{i=1}^{p}\varphi _{i}\,\gamma _{m-i},} これにより次の Yule-Walker方程式が導かれる。 γ m = ∑ i = 1 p φ i γ m − i + σ ε 2 δ m . {\displaystyle \gamma _{m}=\sum _{i=1}^{p}\varphi _{i}\gamma _{m-i}+\sigma _{\varepsilon }^{2}\delta _{m}.} 誤差項 εt は一般に「独立かつ同一の分布に従う」(i.i.d.)無作為変数であり、ゼロを平均値とする正規分布に従う。すなわち εt ~ N(0,σ2) で、σ2 は分散である。このような仮定を弱めることもあるが、そうするとモデルとしての性質が変化する。特に、i.i.d. という仮定を変更すると根本的な性質が変化する。 ARMAモデルをラグ作用素（遅れ作用素） L を使って表す場合もある。この場合、AR(p)モデルは次のように表される。 ε t = ( 1 − ∑ i = 1 p φ i L i ) X t = φ X t {\displaystyle \varepsilon _{t}=\left(1-\sum _{i=1}^{p}\varphi _{i}L^{i}\right)X_{t}=\varphi X_{t}\,} ここで、φ は次の多項式で表される。 φ = 1 − ∑ i = 1 p φ i L i . {\displaystyle \varphi =1-\sum _{i=1}^{p}\varphi _{i}L^{i}.\,} また、MA(q)モデルは次のように表される。 X t = ( 1 + ∑ i = 1 q θ i L i ) ε t = θ ε t {\displaystyle X_{t}=\left(1+\sum _{i=1}^{q}\theta _{i}L^{i}\right)\varepsilon _{t}=\theta \varepsilon _{t}\,} ここで θ は次の多項式で表される。 θ = 1 + ∑ i = 1 q θ i L i . {\displaystyle \theta =1+\sum _{i=1}^{q}\theta _{i}L^{i}.\,} 以上から、ARMA(p, q)モデルは次のように表される。 ( 1 − ∑ i = 1 p φ i L i ) X t = ( 1 + ∑ i = 1 q θ i L i ) ε t {\displaystyle \left(1-\sum _{i=1}^{p}\varphi _{i}L^{i}\right)X_{t}=\left(1+\sum _{i=1}^{q}\theta _{i}L^{i}\right)\varepsilon _{t}\,} あるいは、もっと簡潔に記せば、次のようになる。 φ X t = θ ε t . {\displaystyle \varphi X_{t}=\theta \varepsilon _{t}.\,} ラグ作用素とは、時系列データのある時点のデータで他の時点のデータを表すように係数化したもの。上記の式はいずれも Xt しか出現しない（他の時点のデータが出てこない）ことに注意されたい。他の時点のデータは全てラグ作用素によって表されている。 実データに適用する場合、ARMAモデルの p と q を選択後、誤差項を最小化するパラメータを探るため最小二乗法を使うのが普通である。また、実データに適合する最小の p および q を見つけることでよい結果が得られることが知られている。純粋なARモデルでは、これに Yule-Walker 方程式を利用することができる。 ARMAモデルの一般化として次が挙げられる。 非線型自己回帰移動平均モデル (NARMA): Xt の過去の値や誤差項 εt との依存関係を線形に限定しない 自己回帰条件付き分散変動モデル (ARCH) 自己回帰和分移動平均モデル (ARIMA) ベクトルARIMAモデル 季節ARIMAモデル (SARIMA): 季節変動効果の考慮 多変量自己回帰モデル (MAR) 自己回帰モデル (AR) 移動平均モデル (MA) 自己回帰和分移動平均モデル (ARIMA) 予測分析 放射基底関数