自己回帰モデル
自己回帰モデル(じこかいきモデル、autoregressive model)は時点 t におけるモデル出力が時点 t 以前のモデル出力に依存する確率過程である。ARモデルとも呼ばれる。
自己回帰モデルは、例えば自然科学や経済学において、時間について変動する過程を描写している。(古典的な)自己回帰モデルは実現値となる変数がその変数の過去の値と確率項(確率、つまりその値を完全には予測できない項)に線形に依存している。ゆえに自己回帰モデルは一種の確率差分方程式の形状を取る。
自己回帰モデルはより一般的な時系列の自己回帰移動平均モデル(ARMAモデル)の特別なケースである。また、一つ以上の確率差分方程式からなるベクトル自己回帰モデル(VARモデル)の特別ケースでもある。推計統計学・機械学習における生成モデルとしても自己回帰モデルは表現でき、古典的な(線形)自己回帰生成モデルを拡張した非線形自己回帰生成モデルも盛んに研究されている。
AR(p) という記法はオーダー p の自己回帰モデルを意味している。AR(p) モデルは以下のように定義される。
X_t = c + \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t \,
ここで \varphi_1, \ldots, \varphi_p はモデルのパラメーターであり、c は定数項、 \varepsilon_t はホワイトノイズである。この式は後退オペレーター B を用いることで以下のような同値である表現で書き表すことが出来る。
X_t = c + \sum_{i=1}^p \varphi_i B^i X_t + \varepsilon_t
よって、左辺の総和を移項し多項式表現を用いれば、
\phi (B)X_t= c + \varepsilon_t \, .
と表せる。ゆえに自己回帰モデルは、ホワイトノイズを入力値とする、全ての極における無限インパルス応答の出力値として見なすことも出来る。
自己回帰モデルが弱定常であるためにはいくつかのパラメーター制約が必要になる。例えば、|\varphi_1 | \geq 1 であるAR(1) モデルで表現される過程は定常ではない。より一般的に、AR(p) モデルが弱定常であるためには、多項式 \textstyle z^p - \sum_{i=1}^p \varphi_i z^{p-i} の根が単位円の内側になくてはならない。つまり全ての根 z_i が |z_i|<1 を満たさなくてはならない。
自己回帰モデルにおいて、一時点でのショックは将来の更新変数の値に恒久的に影響を与える。例えば、AR(1) モデル  X_t = c + \varphi_1 X_{t-1} + \varepsilon_t を考えてみよう。t=1 時点での \varepsilon_t の値がゼロでなければ、\varepsilon_1 の量だけ X_1 に影響がある。この時、X_1 から見た X_2 についてのAR方程式により、\varepsilon_1 は \varphi_1 \varepsilon_1 の量だけ X_2 に影響を与える。さらに、X_2 から見た X_3 についてのAR方程式により、\varepsilon_1 は \varphi_1^2 \varepsilon_1 の量だけ X_3 に影響を与える。これを繰り返すことで \varepsilon_1 の効果は永久に波及することが分かる。しかしながら、過程が定常過程ならば、この効果は極限において0となる。
全てのショックが、それが起こった時点から X に恒久的に影響を与えるため、任意の与えられた Xt の値は過去に起こったショック全てから影響を受ける。これは自己回帰方程式
\phi (B)X_t=  \varepsilon_t \,
(ここで定数項は変数が平均からの逸脱として測定されると仮定することで無視できる)が以下のように書き直せることからもまた分かる。
X_t= \frac{1}{\phi (B)}\varepsilon_t \, .
右辺における多項式の除算が可能なのであれば、\varepsilon_t に適用される後退オペレーターによる多項式は無限次元のオーダーを持つ。つまり、\varepsilon_t のラグ値が方程式の右辺において無限個現れる。
AR(p) 過程の自己相関関数は以下のように表すことが出来る。
\rho(\tau) = \sum_{k=1}^p a_k y_k^{-|\tau〕 ,
ここで y_k は以下の多項式の根である。
\phi(B) = 1- \sum_{k=1}^p \varphi_k B^k
ここで B は後退オペレーターであり、\phi(.) は自己回帰を定義する関数、\varphi_k は自己回帰における係数である。
AR(p) 過程の自己相関関数は指数減衰する部分の和となっている。
全ての実数根は指数減衰する自己相関関数の構成要素として寄与する。
同様にすべての複素数の共役根の組は指数的に減衰する循環として寄与する。
最も単純なARモデルは AR(0) であり、項の間に依存関係がない。誤差/イノベーション/ノイズ項のみが過程の出力に寄与し、ゆえに図で示されているように AR(0) はホワイトノイズに対応する。
\varphi の値が正である AR(1) 過程について、その過程の以前の項とノイズ項のみが出力に寄与する。もし \varphi が0に近ければ、その過程は依然としてホワイトノイズのように見える。しかし、\varphi が1に近いならば、出力はノイズに比べて現在の項に大きな影響を受ける。結果として出力の"スムージング"もしくは和分が起こり、ローパスフィルタと似たものとなる。
AR(2) 過程について、以前の二つの項とノイズ項が出力に寄与する。\varphi_1 と \varphi_2 が共に正ならば、出力はノイズの高周波数領域が減衰するローパスフィルタに似通ったものとなる。もし \varphi_1 が正である一方で \varphi_2 が負であれば、過程はその項の間で符号が変わりやすくなる。出力は循環的となる。これは方向におけるエッジ検出もしくは変化検出と結びつけることが出来る。
AR(1) 過程は以下で与えられる。
X_t = c + \varphi X_{t-1}+\varepsilon_t\,
ここで \varepsilon_t は平均0のホワイトノイズ過程でありその分散は定数 \sigma_\varepsilon^2 である(注記: \varphi_1 の下につく添え字をなくしている)。もし |\varphi|<1 ならば、この確率過程は弱定常である。というのもこの過程はホワイトノイズを入力とする定常フィルターの出力として得られるからである(もし \varphi=1 ならば X_t の分散は無限大となり、ゆえに弱定常ではなくなる)。結果として、|\varphi|<1 を仮定すれば、平均 \operatorname{E} (X_t) は全ての t の値について同じとなる。もし平均を \mu と書くのであれば、以下の式
\operatorname{E} (X_t)=\operatorname{E} (c)+\varphi\operatorname{E} (X_{t-1})+\operatorname{E}(\varepsilon_t),
\mu=c+\varphi\mu+0,
が成り立ち、ゆえに以下が得られる。
\mu=\frac{c}{1-\varphi}.
特に、 c = 0 ならば、平均は0である。
分散は以下のように定まる。
\operatorname{var}(X_t)=\operatorname{E}(X_t^2)-\mu^2=\frac{\sigma_\varepsilon^2}{1-\varphi^2},
ここで \sigma_\varepsilon は \varepsilon_t の標準偏差である。これは以下の式、
\operatorname{var}(X_t) = \varphi^2\operatorname{var}(X_{t-1}) + \sigma_\varepsilon^2,
と上の量は安定な不動点となることから示される。
自己共分散は以下で与えられる。
B_n=\operatorname{E}(X_{t+n}X_t)-\mu^2=\frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\,\varphi^.
自己共分散関数は \tau=-1/\ln(\varphi) の減衰時間(または時定数)で減衰していくことが分かる[これを見るために B_n=K\phi^ と書く。ここで K は n に依存しない。この時、\phi^=e^〔n|\ln\phi} であり、指数減衰法則 e^{-n/\tau} と一致する]。
スペクトル密度とは自己共分散関数のフーリエ変換である。離散時間の場合、フーリエ変換は離散時間フーリエ変換に対応する。
\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_\varepsilon^2}{1+\varphi^2-2\varphi\cos(\omega)}\right).
この表現は X_j の離散的性質により周期的となり、それは分母におけるコサイン項によって明らかとなっている。もしサンプリング時間 (\Delta t=1) が減衰時間 (\tau) より非常に小さいと仮定するならば、B_n の連続体近似を用いることが出来る。
B(t)\approx \frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\,\varphi^
これにより、コーシー分布(ローレンツ型プロファイル)のスペクトル密度が得られる。
\Phi(\omega)=\frac{1}{\sqrt{2\pi}\,\frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\frac{\gamma}{\pi(\gamma^2+\omega^2)}
ここで \gamma=1/\tau は減衰時間 \tau に対応した角周波数である。
X_{t-1} についての c+\varphi X_{t-2}+\varepsilon_{t-1} を定義式にまず代入することで、X_t の別表現が得られる。これを N 回繰り返せば
X_t=c\sum_{k=0}^{N-1}\varphi^k+\varphi^NX_{t-N}+\sum_{k=0}^{N-1}\varphi^k\varepsilon_{t-k}.
となる。N を無限大まで発散させれば、\varphi^N は0に近づき、
X_t=\frac{c}{1-\varphi}+\sum_{k=0}^\infty\varphi^k\varepsilon_{t-k}.
となる。X_t は \varphi^k の核で畳み込まれたホワイトノイズに定数の平均を足したものとなることが分かる。もしホワイトノイズ \varepsilon_t がガウス過程ならば X_t もまたガウス過程である。他の場合として、中心極限定理により \varphi が1に近づけば X_t は正規分布に近似的に近づくことが分かる。
AR(1) 過程は連続時間におけるオルンシュタイン=ウーレンベック過程の離散時間のアナロジーである。ゆえに AR(1) モデルの性質を理解するために同様の形式に変換することが時として有用になる。この形式において AR(1) モデルは以下で与えられる。
X_{t+1} = X_t + \theta(\mu - X_t) + \epsilon_{t+1}\,
ここで |\theta| < 1 \, であり \mu はモデルの平均である。これを  X_{t+1} = c + \phi X_t\,  の式に当てはめ X_{t+n}\, についての系列に展開することで次が示される。
\operatorname{E}(X_{t+n} | X_t) = \mu\left[1-\left(1-\theta\right)^n\right] + X_t(1-\theta)^n, and
\operatorname{Var} (X_{t+n} | X_t) = \sigma^2 \frac{ \left[ 1 - (1 - \theta)^{2n} \right] }{1 - (1-\theta)^2}.
AR(p) 過程の偏自己相関はラグが p + 1 より大きい時にゼロとなり、結果として適切な最大ラグはそのラグより大きいラグでの偏自己相関が全てゼロになるものである。
ARモデルの係数の推定には多数の方法があり、例えば最小二乗法の手続きや、もしくは(ユール–ウォーカー方程式(Yule–Walker equation)を通した)モーメント法がある。
AR(p) モデルは以下の方程式で与えられる。
X_t = \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t.\,
この方程式はパラメーター \varphi_i i = 1, ..., p に基いている。これらのパラメーターと過程の共分散関数の間には直接的な対応が存在し、その対応は(共分散から得られる)自己相関関数からパラメーターを決定する為に裏返すことができる。これはユール–ウォーカー方程式を用いて行われる。
ユール–ウォーカー方程式は、ウドニー・ユールとギルバート・ウォーカーにちなんで名づけられたもので、以下の方程式からなる。
\gamma_m = \sum_{k=1}^p \varphi_k \gamma_{m-k} + \sigma_\varepsilon^2\delta_{m,0},
ここで m = 0, ..., p であり、p + 1 個の方程式からなる。さらに \gamma_m は Xt の自己共分散関数、\sigma_\varepsilon は入力ホワイトノイズの標準偏差、\delta_{m,0} はクロネッカーのデルタである。
各方程式の最後の部分がゼロとならないのは m = 0 の時に限られるので、この方程式は m > 0 の方程式を行列形式に表すことで解くことが出来る。よって次の方程式が得られる。
\begin{bmatrix}
\gamma_1 \\
\gamma_2 \\
\gamma_3 \\
\vdots \\
\gamma_p \\
\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 \\
\vdots      & \vdots         & \vdots       & \ddots \\
\gamma_{p-1} & \gamma_{p-2} & \gamma_{p-3} & \dots \\
\end{bmatrix}
\begin{bmatrix}
\varphi_{1} \\
\varphi_{2} \\
\varphi_{3} \\
\vdots \\
\varphi_{p} \\
\end{bmatrix}
これは全ての \{\varphi_m; m=1,2, \cdots ,p\} について解くことが出来る。残りの m = 0 についての方程式は
\gamma_0 = \sum_{k=1}^p \varphi_k \gamma_{-k} + \sigma_\varepsilon^2 ,
となり、一度 \{\varphi_m ; m=1,2, \cdots ,p \} の値を知ってしまえば、この方程式を \sigma_\varepsilon^2  について解くことが出来る。
他の定式化として自己相関についてのものがある。ARパラメーターは自己相関 \rho(\tau) の最初の p + 1 個の要素で決定する。完全な自己相関関数はこの時、再帰的な計算によって導出できる。
\rho(\tau) = \sum_{k=1}^p \varphi_k \rho(k-\tau)
幾つかの低い次数の AR(p) 過程についての例は
p=1
\gamma_1 = \varphi_1 \gamma_0
ゆえに \rho_1 = \gamma_1 / \gamma_0 = \varphi_1
p=2
AR(2) 過程のユール–ウォーカー方程式は
\gamma_1 = \varphi_1 \gamma_0 + \varphi_2 \gamma_{-1}
\gamma_2 = \varphi_1 \gamma_1 + \varphi_2 \gamma_0
\gamma_{-k} = \gamma_k であることを思い出せば、
第一の方程式を用いることで \rho_1 = \gamma_1 / \gamma_0 = \frac{\varphi_1}{1-\varphi_2} となり、
第二の方程式を用いることで \rho_2 = \gamma_2 / \gamma_0 = \frac{\varphi_1^2 - \varphi_2^2 + \varphi_2}{1-\varphi_2} となる。
上の方程式(ユール–ウォーカー方程式)は、理論的な共分散を推定値に置き換えることで、 AR(p) モデルのパラメーターを推定する為にいくつかの方法を提供する。下記のような方法が考えられる。
自己共分散もしくは自己相関の推定。便利な推定法を用いて自己共分散もしくは自己相関の項のそれぞれを分割して推定したものとする。推定の方法は多様であり、どれを選択するかは推定のスキームが持つ性質に影響を与える。例えば、ある方法では分散の負の推定量が生じうる。
Xt の予測値を同じ系列の過去の p 個の値として基礎づける最小二乗予測問題を構築する上での最小二乗回帰問題としての定式化。これは前方予測スキームとして考えられる。この問題についての正規方程式は同じラグで現れる自己共分散を少し違った推定値で置き換えたユール–ウォーカー方程式の行列形式の近似と対応するように見える。
最小二乗予測問題の拡張形式としての定式化。ここで二つの予測方程式のセットを一つの推定スキームと単一の正規方程式に結合する。一つのセットは前方予測方程式のセットとなっており、もう片方は対応する後方予測方程式のセットとなっている。これはARモデルの後退表現と関連している。
X_t = c + \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon^*_t \,.
ここで予測値 Xt は同じ系列の p 個の将来の値に基づいている。このARパラメーターの推定方法はジョン・バーグ(John P. Burg) によるものでバーグの方法(the Burg method)と呼ばれる。バーグや後続の研修者はこの特別な推定値を"最大エントロピー推定量"と呼ぶが、この背後にある理論は推定パラメータ―のどのようなセットについても適用できる。前進予測方程式のみを用いた推定スキームと比べると、異なる自己共分散の推定値が得られ、推定量は異なる安定性の性質を持つ。バーグ推定量は特に最大エントロピースペクトル推定と関連している。
他の考えられる方法として最尤法がある。異なる二つの最尤法が利用できる。一つは(大体は最小二乗スキームによる前方予測と同じだが)考慮する尤度関数を系列における当初の p 此の値を所与とした系列の後の値の条件つき分布に対応させるものである。もう一つは考慮する尤度関数を観測された系列の全ての値の無条件の同時分布に対応させるものである。これらの方法の結果における本質的な違いは観測系列が短い、もしくは過程が非定常に近い時に現れる。
ノイズ分散が \mathrm{Var}(Z_t) = \sigma_Z^2 である AR(p) 過程のパワースペクトル密度は以下のようになる。
S(f) = \frac{\sigma_Z^2}〔 1-\sum_{k=1}^p \varphi_k e^{-2 \pi i k f} |^2}.
ホワイトノイズ (AR(0)) については以下のようになる。
S(f) = \sigma_Z^2.
AR(1) については以下のようになる。
S(f) = \frac{\sigma_Z^2}〔 1- \varphi_1 e^{-2 \pi i f} |^2}
= \frac{\sigma_Z^2}{ 1 + \varphi_1^2 - 2 \varphi_1 \cos 2 \pi f }
もし \varphi_1 > 0 ならば、スペクトルは f = 0 において単峰で、レッドノイズと呼ばれる。\varphi_1  が1に近ければ低周波においてパワーが強くなる。つまり時間のラグが大きくなる。これはローパスフィルタであり、フルスペクトル光に適用された時、赤の波長を除いてすべてがフィルタリングされる。
もし \varphi_1 < 0 ならば、スペクトルは f = 0 において最小値を取り、ブルーノイズと呼ばれる。これはハイパスフィルタのように振る舞い、青の波長を除いてすべてがフィルタリングされる。
AR(2) 過程は特性方程式の根に依存する3つのグループに分割される。
z_1,z_2 = \frac{1}{2}\left(\varphi_1 \pm \sqrt{\varphi_1^2 + 4\varphi_2}\right)
\varphi_1^2 + 4\varphi_2 < 0 の時、過程は複素共役根のペアを一つ持ち、中周波でピークを作る。
f^* = \frac{1}{2\pi}\cos^{-1}\left(\frac{\varphi_1(\varphi_2-1)}{4\varphi_2}\right)
そうでなければ実数根を持ち、
\varphi_1 > 0 の時、f=0 で頂点を持つホワイトノイズに対するローパスフィルタのように振る舞い、
\varphi_1 < 0 の時、f=1/2 で頂点を持つホワイトノイズに対するハイパスフィルタのように振る舞う。
根が単位円の外側にある時、この過程は定常である。
根が単位円の内側にある、もしくは同じことだが係数が三角形 -1 \le \varphi_2 \le 1 - |\varphi_1| の内部にある時、安定である。
完全なパワースペクトル密度関数は以下のように表される。
S(f) = \frac{\sigma_Z^2}{1 + \varphi_1^2 + \varphi_2^2 - 2\varphi_1(1-\varphi_2)\cos(2\pi f) - 2\varphi_2\cos(4\pi f)}
R, stats パッケージに ar 関数が含まれている。
MATLAB, Econometric Toolbox と System Identification Toolbox に自己回帰モデルが含まれている。
MATLABとOctave: TSA toolbox に単一変数、複数変数、適応自己回帰モデルについてのいくつかの推定関数が含まれている。
X_t = c + \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t \,
のパラメーターが一度推定されてしまえば、この自己回帰は将来の任意の時点での予測に用いることが出来る。まず、t をデータが使えない最初の時点とする。既知の値Xt-i for i=1, ..., p を自己回帰方程式に代入し、誤差項 \varepsilon_t をゼロと置く(なぜならば Xt をその条件つき期待値と一致させるように予測し、観測されない誤差項の期待値は0であるから)ことで予測ができる。自己回帰方程式の出力は最初のデータが観測されない時点についての予測となる。次に、t をデータが使えない 次の 時点とする。もう一度自己回帰方程式を予測を作るために使うことができる。ただし一つ異なる点がある。X の今予測している時点より一期前の値は未知である。よってその期待値、つまり前の予測ステップでの予測値を代わりに用いる。この時、将来の時点において同じ手続きが用いられ、p 回の予測の後に、全ての p 個の右辺の値が事前のステップによる予測値となるまで、予測方程式の右辺における予測値を用いる。
この方法で得られた予測値について四つの不確実性のソースがある。(1) 自己回帰モデルが正しいモデルかどうかという不確実性、(2) 自己回帰方程式の右辺においてラグ値として用いられる予測値の正しさについての不確実性、(3) 自己回帰係数の真の値についての不確実性、(4) 予測機関における誤差項 \varepsilon_t \, の値についての不確実性である。最後の三つは定量化可能で n ステップ後の予測についての信頼区間として与えられる。右辺の変数についての推定値が増えるため信頼区間は n が増えれば広くなる。
自己回帰モデルの予測性能は、クロス・バリデーションが行われるならば、推定の後に即座に評価できる。この方法においては、最初の方の利用可能なデータはパラメーターの推定の為に用いられ、データセットにおける後の方のデータはアウトオブサンプルのテストとして残しておく。他には、パラメーター推定が行われた後にしばらくしたあと、より多くのデータが利用可能になり予測性能を新しいデータを使うことで評価できる。
どちらのケースも、評価可能な予測性能には2つの側面がある。1期先予測の性能と n 期先予測の性能である。1期先予測の性能について、推定パラメーターは予測を行った期以前の全ての期における X の観測値と共に自己回帰方程式が用いられ、方程式の出力は1期先予測となる。この手続きはアウトオブサンプルの観測値についての予測を得るために用いられる。n 期先予測の質を評価する為には、予測を得るために前の節での予測手続きが用いられる。
予測値のセットと対応する様々な期間の X の本当の値のセットが与えられたとして一般的な評価のテクニックは平均二乗予測誤差を用いることである。他の尺度もまた用いられる。
ここで測定された予測の正しさをどのように解釈するのかという問題が持ち上がる。例えば平均二乗予測誤差が"高い"(悪い)もしくは"低い"(良い)とはどういう事なのだろうか。比較の上で二つのポイントがある。第一に他のモデルの仮定もしくは推定手法の下で推定された代替モデルの予測の正しさは比較目的に使用できる。第二にアウトオブサンプルの正確さの尺度は十分に前のデータを用いることが出来るならば、つまり最初の p 個のデータポイントを落として p 期以前のデータを使わないならば(パラメーター推定に用いられた)インサンプルのデータポイントでの同じ尺度と比較できる。モデルはインサンプルのデータポイントに出来るだけ適合するように特定化されて推定されるので、普通はアウトオブサンプルの予測性能はインサンプルの予測性能より悪い。しかし予測の質がアウトオブサンプルで(正確には定義できないが)"そう悪くない"のであれば、予測値は十分なパフォーマンスを見せていると言える。
上記のように自己回帰モデルは、決定論的/deterministicな線形変換に、確率的/probabilisticなブレ(確率項)が線形に追加されるモデルである。別の表現として、自己回帰モデルは統計モデル(生成モデル)で表すことができる。
自己回帰モデル AR(n) を考えるとし、因果関係を持つ時点 t から時刻 t-n までの値の組を
X_{n}=(x_{t}, x_{t-1}, ..., x_{t-n})
とする。Xn の確率分布すなわち xt ~ xt-nの同時確率 pAR(n) は
p_{AR(n)}(X_{n})=p_{AR(n)}(x_{t}, x_{t-1}, ..., x_{t-n})
であり、条件付き確率の定義を用いる(確率の連鎖律)と