自己回帰モデルのソースを表示

'''自己回帰モデル'''（じこかいきモデル、{{lang-en-short|autoregressive model}}）は時点 t におけるモデル出力が時点 t 以前のモデル出力に依存する[[確率過程]]である。'''ARモデル'''とも呼ばれる。

自己回帰モデルは、例えば[[自然科学]]や[[経済学]]において、時間について変動する過程を描写している。（古典的な）自己回帰モデルは実現値となる変数がその変数の過去の値と'''確率項'''（確率、つまりその値を完全には予測できない項）に'''線形に'''依存している。ゆえに自己回帰モデルは一種の確率差分方程式の形状を取る。

自己回帰モデルはより一般的な時系列の[[自己回帰移動平均モデル]]（ARMAモデル）の特別なケースである。また、一つ以上の確率差分方程式からなる{{仮リンク|ベクトル自己回帰モデル|en|vector autoregressive model}}（VARモデル）の特別ケースでもある。[[推計統計学]]・[[機械学習]]における生成モデルとしても自己回帰モデルは表現でき、古典的な（線形）自己回帰生成モデルを拡張した'''非線形自己回帰生成モデル'''も盛んに研究されている。

==定義==
<math>AR(p)</math> という記法はオーダー ''p'' の自己回帰モデルを意味している。AR(''p'') モデルは以下のように定義される。
:<math> X_t = c + \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t \,</math>
ここで <math>\varphi_1, \ldots, \varphi_p</math> はモデルの'''''パラメーター'''''であり、<math>c</math> は定数項、 <math>\varepsilon_t</math> は[[ホワイトノイズ]]である。この式は{{仮リンク|後退オペレーター|en|backshift operator}} ''B'' を用いることで以下のような同値である表現で書き表すことが出来る。
:<math> X_t = c + \sum_{i=1}^p \varphi_i B^i X_t + \varepsilon_t </math>
よって、左辺の総和を移項し[[多項式]]表現を用いれば、
:<math>\phi (B)X_t= c + \varepsilon_t \, .</math>
と表せる。ゆえに自己回帰モデルは、ホワイトノイズを入力値とする、全ての[[特異点|極]]における[[無限インパルス応答]]の出力値として見なすことも出来る。

自己回帰モデルが[[定常過程|弱定常]]であるためにはいくつかのパラメーター制約が必要になる。例えば、<math>|\varphi_1 | \geq 1</math> であるAR(1) モデルで表現される過程は定常ではない。より一般的に、AR(''p'') モデルが弱定常であるためには、多項式 <math>\textstyle z^p - \sum_{i=1}^p \varphi_i z^{p-i}</math> の根が[[単位円]]の内側になくてはならない。つまり全ての根 <math>z_i</math> が <math>|z_i|<1</math> を満たさなくてはならない。

==ショックの異時点間における影響==
自己回帰モデルにおいて、一時点でのショックは将来の更新変数の値に恒久的に影響を与える。例えば、AR(1) モデル <math> X_t = c + \varphi_1 X_{t-1} + \varepsilon_t</math> を考えてみよう。''t''=1 時点での <math>\varepsilon_t</math> の値がゼロでなければ、<math>\varepsilon_1</math> の量だけ <math>X_1</math> に影響がある。この時、<math>X_1</math> から見た <math>X_2</math> についてのAR方程式により、<math>\varepsilon_1</math> は <math>\varphi_1 \varepsilon_1</math> の量だけ <math>X_2</math> に影響を与える。さらに、<math>X_2</math> から見た <math>X_3</math> についてのAR方程式により、<math>\varepsilon_1</math> は <math>\varphi_1^2 \varepsilon_1</math> の量だけ <math>X_3</math> に影響を与える。これを繰り返すことで <math>\varepsilon_1</math> の効果は永久に波及することが分かる。しかしながら、過程が[[定常過程]]ならば、この効果は極限において0となる。

全てのショックが、それが起こった時点から ''X'' に恒久的に影響を与えるため、任意の与えられた ''X''<sub>''t''</sub> の値は過去に起こったショック全てから影響を受ける。これは自己回帰方程式
:<math>\phi (B)X_t=  \varepsilon_t \,</math>
（ここで定数項は変数が平均からの逸脱として測定されると仮定することで無視できる）が以下のように書き直せることからもまた分かる。
:<math>X_t= \frac{1}{\phi (B)}\varepsilon_t \, .</math>
右辺における[[多項式]]の除算が可能なのであれば、<math>\varepsilon_t</math> に適用される後退オペレーターによる多項式は無限次元のオーダーを持つ。つまり、<math>\varepsilon_t</math> のラグ値が方程式の右辺において無限個現れる。

==特性多項式==
AR(''p'') 過程の[[自己相関]]関数は以下のように表すことが出来る<ref>{{Citation
|title = Estimation of parameters for a linear difference equation with application to EEG analysis
|journal = Mathematical Biosciences
|volume = 5
|number = 3
|pages = 227--275
|year = 1969
|issn = 0025-5564
|doi = 10.1016/0025-5564(69)90044-3
|first = Lars H.
|last = Zetterberg
}}</ref>。
:<math>\rho(\tau) = \sum_{k=1}^p a_k y_k^{-|\tau|} ,</math>
ここで <math>y_k</math> は以下の多項式の根である。
:<math>\phi(B) = 1- \sum_{k=1}^p \varphi_k B^k </math>
ここで ''B'' は{{仮リンク|後退オペレーター|en|backshift operator}}であり、<math>\phi(.)</math> は自己回帰を定義する関数、<math>\varphi_k</math> は自己回帰における係数である。

AR(''p'') 過程の自己相関関数は指数減衰する部分の和となっている。
* 全ての実数根は指数減衰する自己相関関数の構成要素として寄与する。
* 同様にすべての複素数の共役根の組は指数的に減衰する循環として寄与する。

==AR(''p'') 過程のグラフ==

[[File:ArTimeSeries.svg|thumb|right|alt="Figure has 5 plots of AR processes. AR(0) and AR(0.3) are white noise or look like white noise. AR(0.9) has some large scale oscillating structure."|AR(0); AR(1) with AR parameter 0.3; AR(1) with AR parameter 0.9; AR(2) with AR parameters 0.3 and 0.3; and AR(2) with AR parameters 0.9 and −0.8]]
最も単純なARモデルは AR(0) であり、項の間に依存関係がない。誤差/イノベーション/ノイズ項のみが過程の出力に寄与し、ゆえに図で示されているように AR(0) はホワイトノイズに対応する。

<math>\varphi</math> の値が正である AR(1) 過程について、その過程の以前の項とノイズ項のみが出力に寄与する。もし <math>\varphi</math> が0に近ければ、その過程は依然としてホワイトノイズのように見える。しかし、<math>\varphi</math> が1に近いならば、出力はノイズに比べて現在の項に大きな影響を受ける。結果として出力の"スムージング"もしくは和分が起こり、[[ローパスフィルタ]]と似たものとなる。

AR(2) 過程について、以前の二つの項とノイズ項が出力に寄与する。<math>\varphi_1</math> と <math>\varphi_2</math> が共に正ならば、出力はノイズの高周波数領域が減衰するローパスフィルタに似通ったものとなる。もし <math>\varphi_1</math> が正である一方で <math>\varphi_2</math> が負であれば、過程はその項の間で符号が変わりやすくなる。出力は循環的となる。これは方向におけるエッジ検出もしくは変化検出と結びつけることが出来る。

==例: AR(1) 過程==
AR(1) 過程は以下で与えられる。
:<math>X_t = c + \varphi X_{t-1}+\varepsilon_t\,</math>
ここで <math>\varepsilon_t</math> は平均0のホワイトノイズ過程でありその分散は定数 <math>\sigma_\varepsilon^2</math> である（注記: <math>\varphi_1</math> の下につく添え字をなくしている）。もし <math>|\varphi|<1</math> ならば、この確率過程は[[定常過程|弱定常]]である。というのもこの過程はホワイトノイズを入力とする定常フィルターの出力として得られるからである（もし <math>\varphi=1</math> ならば <math>X_t</math> の分散は無限大となり、ゆえに弱定常ではなくなる）。結果として、<math>|\varphi|<1</math> を仮定すれば、平均 <math>\operatorname{E} (X_t)</math> は全ての ''t'' の値について同じとなる。もし平均を <math>\mu</math> と書くのであれば、以下の式
:<math>\operatorname{E} (X_t)=\operatorname{E} (c)+\varphi\operatorname{E} (X_{t-1})+\operatorname{E}(\varepsilon_t),</math>
より次の式
:<math> \mu=c+\varphi\mu+0,</math>
が成り立ち、ゆえに以下が得られる。
:<math>\mu=\frac{c}{1-\varphi}.</math>
特に、 <math>c = 0</math> ならば、平均は0である。

[[分散 (確率論)|分散]]は以下のように定まる。
:<math>\operatorname{var}(X_t)=\operatorname{E}(X_t^2)-\mu^2=\frac{\sigma_\varepsilon^2}{1-\varphi^2},</math>
ここで <math>\sigma_\varepsilon</math> は <math>\varepsilon_t</math> の標準偏差である。これは以下の式、
:<math>\operatorname{var}(X_t) = \varphi^2\operatorname{var}(X_{t-1}) + \sigma_\varepsilon^2,</math>
と上の量は安定な不動点となることから示される。

[[自己共分散]]は以下で与えられる。
:<math>B_n=\operatorname{E}(X_{t+n}X_t)-\mu^2=\frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\,\varphi^{|n|}.</math>
自己共分散関数は <math>\tau=-1/\ln(\varphi)</math> の減衰時間（または[[時定数]]）で減衰していくことが分かる[これを見るために <math>B_n=K\phi^{|n|}</math> と書く。ここで <math>K</math> は <math>n</math> に依存しない。この時、<math>\phi^{|n|}=e^{|n|\ln\phi}</math> であり、指数減衰法則 <math>e^{-n/\tau}</math> と一致する]。

[[スペクトル密度]]とは自己共分散関数の[[フーリエ変換]]である。離散時間の場合、フーリエ変換は離散時間フーリエ変換に対応する。
:<math>\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).
</math>
この表現は <math>X_j</math> の離散的性質により周期的となり、それは分母におけるコサイン項によって明らかとなっている。もしサンプリング時間 (<math>\Delta t=1</math>) が減衰時間 (<math>\tau</math>) より非常に小さいと仮定するならば、<math>B_n</math> の連続体近似を用いることが出来る。
:<math>B(t)\approx \frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\,\varphi^{|t|}</math>
これにより、[[コーシー分布]]（ローレンツ型プロファイル）のスペクトル密度が得られる。
:<math>\Phi(\omega)=\frac{1}{\sqrt{2\pi}}\,\frac{\sigma_\varepsilon^2}{1-\varphi^2}\,\frac{\gamma}{\pi(\gamma^2+\omega^2)}</math>
ここで <math>\gamma=1/\tau</math> は減衰時間 <math>\tau</math> に対応した角周波数である。

<math>X_{t-1}</math> についての <math>c+\varphi X_{t-2}+\varepsilon_{t-1}</math> を定義式にまず代入することで、<math>X_t</math> の別表現が得られる。これを ''N'' 回繰り返せば
:<math>X_t=c\sum_{k=0}^{N-1}\varphi^k+\varphi^NX_{t-N}+\sum_{k=0}^{N-1}\varphi^k\varepsilon_{t-k}.</math>
となる。''N'' を無限大まで発散させれば、<math>\varphi^N</math> は0に近づき、
:<math>X_t=\frac{c}{1-\varphi}+\sum_{k=0}^\infty\varphi^k\varepsilon_{t-k}.</math>
となる。<math>X_t</math> は <math>\varphi^k</math> の核で畳み込まれたホワイトノイズに定数の平均を足したものとなることが分かる。もしホワイトノイズ <math>\varepsilon_t</math> が[[ガウス過程]]ならば <math>X_t</math> もまたガウス過程である。他の場合として、[[中心極限定理]]により <math>\varphi</math> が1に近づけば <math>X_t</math> は正規分布に近似的に近づくことが分かる。

=== AR(1) 過程の解析的な平均と差分の形式 ===
AR(1) 過程は連続時間における[[オルンシュタイン＝ウーレンベック過程]]の離散時間のアナロジーである。ゆえに AR(1) モデルの性質を理解するために同様の形式に変換することが時として有用になる。この形式において AR(1) モデルは以下で与えられる。
:<math>X_{t+1} = X_t + \theta(\mu - X_t) + \epsilon_{t+1}\,</math>
ここで <math>|\theta| < 1 \,</math> であり <math>\mu</math> はモデルの平均である。これを <math> X_{t+1} = c + \phi X_t\, </math> の式に当てはめ <math>X_{t+n}\,</math> についての系列に展開することで次が示される。
:<math> \operatorname{E}(X_{t+n} | X_t) = \mu\left[1-\left(1-\theta\right)^n\right] + X_t(1-\theta)^n</math>, and
:<math> \operatorname{Var} (X_{t+n} | X_t) = \sigma^2 \frac{ \left[ 1 - (1 - \theta)^{2n} \right] }{1 - (1-\theta)^2}</math>.

==最大ラグの選択==
AR(''p'') 過程の偏自己相関はラグが p + 1 より大きい時にゼロとなり、結果として適切な最大ラグはそのラグより大きいラグでの偏自己相関が全てゼロになるものである。

==ARパラメーターの計算==
ARモデルの係数の推定には多数の方法があり、例えば[[最小二乗法]]の手続きや、もしくは（ユール–ウォーカー方程式（{{lang-en-short|Yule–Walker equation}}）を通した）[[モーメント法]]がある。

AR(''p'') モデルは以下の方程式で与えられる。
:<math> X_t = \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t.\,</math>
この方程式はパラメーター <math>\varphi_i</math> ''i'' = 1, ..., ''p'' に基いている。これらのパラメーターと過程の共分散関数の間には直接的な対応が存在し、その対応は（共分散から得られる）自己相関関数からパラメーターを決定する為に裏返すことができる。これはユール–ウォーカー方程式を用いて行われる。

===ユール–ウォーカー方程式===
<!-- this heading is linked from other articles -->
ユール–ウォーカー方程式は、{{仮リンク|ウドニー・ユール|en|Udny Yule}}と[[ギルバート・ウォーカー]]にちなんで名づけられたもので<ref>{{Citation
|last = Yule
|first = G. Udny
|year = 1927
|url = http://visualiseur.bnf.fr/Visualiseur?Destination=Gallica&O=NUMM-56031
|title = On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers
|journal = Philosophical Transactions of the Royal Society of London, Ser. A
|volume = 226
|pages = 267–298}}</ref><ref>{{Citation
|last = Walker
|first = Gilbert
|year = 1931
|url = http://visualiseur.bnf.fr/Visualiseur?Destination=Gallica&O=NUMM-56224
|title = On Periodicity in Series of Related Terms
|journal = Proceedings of the Royal Society of London, Ser. A
|volume = 131
|pages = 518–532}}</ref>、以下の方程式からなる<ref name="Hamilton59">{{Harvnb|Hamilton|(1994)|page=59|ref=Hamilton1994}}</ref>。
:<math>\gamma_m = \sum_{k=1}^p \varphi_k \gamma_{m-k} + \sigma_\varepsilon^2\delta_{m,0},</math>
ここで {{nowrap|''m'' &#61; 0, ..., ''p''}} であり、{{nowrap|''p'' + 1}} 個の方程式からなる。さらに <math>\gamma_m</math> は X<sub>t</sub> の自己共分散関数、<math>\sigma_\varepsilon</math> は入力ホワイトノイズの標準偏差、<math>\delta_{m,0}</math> は[[クロネッカーのデルタ]]である。

各方程式の最後の部分がゼロとならないのは {{nowrap|''m'' &#61; 0}} の時に限られるので、この方程式は {{nowrap|''m'' > 0}} の方程式を行列形式に表すことで解くことが出来る。よって次の方程式が得られる。
:<math>\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}

</math>
これは全ての <math>\{\varphi_m; m=1,2, \cdots ,p\}</math> について解くことが出来る。残りの ''m'' = 0 についての方程式は
:<math>
\gamma_0 = \sum_{k=1}^p \varphi_k \gamma_{-k} + \sigma_\varepsilon^2 ,
</math>
となり、一度 <math>\{\varphi_m ; m=1,2, \cdots ,p \}</math> の値を知ってしまえば、この方程式を <math>\sigma_\varepsilon^2 </math> について解くことが出来る。

他の定式化として[[自己相関]]についてのものがある。ARパラメーターは自己相関 <math>\rho(\tau)</math> の最初の p + 1 個の要素で決定する。完全な自己相関関数はこの時、再帰的な計算によって導出できる<ref name="Hamilton59" /><ref name=Storch>{{Cite book
| publisher = Cambridge Univ Pr
| isbn = 0-521-01230-9
| last = Von Storch
| first = H.
|author2=F. W Zwiers
| title = Statistical analysis in climate research
| year = 2001
}}{{Page needed|date=March 2011}}</ref>。
: <math>\rho(\tau) = \sum_{k=1}^p \varphi_k \rho(k-\tau)</math>
幾つかの低い次数の AR(''p'') 過程についての例は
* p=1
** <math>\gamma_1 = \varphi_1 \gamma_0</math>
** ゆえに <math>\rho_1 = \gamma_1 / \gamma_0 = \varphi_1</math>
* p=2
** AR(2) 過程のユール–ウォーカー方程式は
**: <math>\gamma_1 = \varphi_1 \gamma_0 + \varphi_2 \gamma_{-1}</math>
**: <math>\gamma_2 = \varphi_1 \gamma_1 + \varphi_2 \gamma_0</math>
*** <math>\gamma_{-k} = \gamma_k</math> であることを思い出せば、
*** 第一の方程式を用いることで <math>\rho_1 = \gamma_1 / \gamma_0 = \frac{\varphi_1}{1-\varphi_2}</math> となり、
*** 第二の方程式を用いることで <math>\rho_2 = \gamma_2 / \gamma_0 = \frac{\varphi_1^2 - \varphi_2^2 + \varphi_2}{1-\varphi_2}</math> となる。

===ARパラメーターの推定===
上の方程式（ユール–ウォーカー方程式）は、理論的な共分散を推定値に置き換えることで、 AR(''p'') モデルのパラメーターを推定する為にいくつかの方法を提供する<ref>{{Harvnb|Hamilton|(1994)|ref=Hamilton1994}}, Chapter 3 and 5</ref>。下記のような方法が考えられる。
* 自己共分散もしくは自己相関の推定。便利な推定法を用いて自己共分散もしくは自己相関の項のそれぞれを分割して推定したものとする。推定の方法は多様であり、どれを選択するかは推定のスキームが持つ性質に影響を与える。例えば、ある方法では分散の負の推定量が生じうる。
* ''X''<sub>''t''</sub> の予測値を同じ系列の過去の ''p'' 個の値として基礎づける最小二乗予測問題を構築する上での[[最小二乗法|最小二乗回帰]]問題としての定式化。これは前方予測スキームとして考えられる。この問題についての{{仮リンク|正規方程式|en|normal equation}}は同じラグで現れる自己共分散を少し違った推定値で置き換えたユール–ウォーカー方程式の行列形式の近似と対応するように見える。
* 最小二乗予測問題の拡張形式としての定式化。ここで二つの予測方程式のセットを一つの推定スキームと単一の正規方程式に結合する。一つのセットは前方予測方程式のセットとなっており、もう片方は対応する後方予測方程式のセットとなっている。これはARモデルの後退表現と関連している。
::<math> X_t = c + \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon^*_t \,.</math>
:ここで予測値 ''X''<sub>''t''</sub> は同じ系列の ''p'' 個の将来の値に基づいている。このARパラメーターの推定方法はジョン・バーグ（John P. Burg）<ref name=Burg/> によるものでバーグの方法（{{lang-en-short|the Burg method}}）と呼ばれる<ref name=Brockwell/>。バーグや後続の研修者はこの特別な推定値を"最大エントロピー推定量"と呼ぶが<ref name=Burg1/>、この背後にある理論は推定パラメータ―のどのようなセットについても適用できる。前進予測方程式のみを用いた推定スキームと比べると、異なる自己共分散の推定値が得られ、推定量は異なる安定性の性質を持つ。バーグ推定量は特に{{仮リンク|最大エントロピースペクトル推定|en|maximum entropy spectral estimation}}と関連している<ref name=Bos/>。

他の考えられる方法として[[最尤法]]がある。異なる二つの最尤法が利用できる。一つは（大体は最小二乗スキームによる前方予測と同じだが）考慮する尤度関数を系列における当初の ''p'' 此の値を所与とした系列の後の値の条件つき分布に対応させるものである。もう一つは考慮する尤度関数を観測された系列の全ての値の無条件の同時分布に対応させるものである。これらの方法の結果における本質的な違いは観測系列が短い、もしくは過程が非定常に近い時に現れる。

==スペクトル==
[[File:AutocorrTimeAr.svg|thumb|right]]
[[File:AutoCorrAR.svg|thumb|right]]
ノイズ分散が <math>\mathrm{Var}(Z_t) = \sigma_Z^2</math> である AR(''p'') 過程の[[スペクトル密度#パワースペクトル密度|パワースペクトル密度]]は以下のようになる<ref>{{Harvnb|Hamilton|(1994)|page=155|ref=Hamilton1994}}</ref><ref name=Storch />。
: <math>S(f) = \frac{\sigma_Z^2}{| 1-\sum_{k=1}^p \varphi_k e^{-2 \pi i k f} |^2}.</math>

===AR(0)===
ホワイトノイズ (AR(0)) については以下のようになる。
: <math>S(f) = \sigma_Z^2.</math>

===AR(1)===
AR(1) については以下のようになる。
: <math>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 }</math>
* もし <math>\varphi_1 > 0</math> ならば、スペクトルは f = 0 において単峰で、[[レッドノイズ]]と呼ばれる。<math>\varphi_1 </math> が1に近ければ低周波においてパワーが強くなる。つまり時間のラグが大きくなる。これはローパスフィルタであり、フルスペクトル光に適用された時、赤の波長を除いてすべてがフィルタリングされる。
* もし <math>\varphi_1 < 0</math> ならば、スペクトルは f = 0 において最小値を取り、[[カラードノイズ#ブルーノイズ|ブルーノイズ]]と呼ばれる。これはハイパスフィルタのように振る舞い、青の波長を除いてすべてがフィルタリングされる。

===AR(2)===
AR(2) 過程は特性方程式の根に依存する3つのグループに分割される。
:<math>z_1,z_2 = \frac{1}{2}\left(\varphi_1 \pm \sqrt{\varphi_1^2 + 4\varphi_2}\right)</math>
* <math>\varphi_1^2 + 4\varphi_2 < 0</math> の時、過程は複素共役根のペアを一つ持ち、中周波でピークを作る。
:<math>f^* = \frac{1}{2\pi}\cos^{-1}\left(\frac{\varphi_1(\varphi_2-1)}{4\varphi_2}\right)</math>
そうでなければ実数根を持ち、
* <math>\varphi_1 > 0</math> の時、<math>f=0</math> で頂点を持つホワイトノイズに対するローパスフィルタのように振る舞い、
* <math>\varphi_1 < 0</math> の時、<math>f=1/2</math> で頂点を持つホワイトノイズに対するハイパスフィルタのように振る舞う。
根が単位円の外側にある時、この過程は定常である。
根が単位円の内側にある、もしくは同じことだが係数が三角形 <math>-1 \le \varphi_2 \le 1 - |\varphi_1|</math> の内部にある時、安定である。
完全なパワースペクトル密度関数は以下のように表される。
:<math>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)}</math>

==統計パッケージにおける実装==
* [[R言語|R]], ''stats'' パッケージに ''ar'' 関数が含まれている<ref>[http://finzi.psych.upenn.edu/R/library/stats/html/ar.html "Fit Autoregressive Models to Time Series"] (in R)</ref>。
* [[MATLAB]], Econometric Toolbox<ref>[http://www.mathworks.com/products/econometrics/ Econometrics Toolbox Overview]</ref> と System Identification Toolbox<ref>[http://www.mathworks.com/products/sysid/ System Identification Toolbox overview]</ref> に自己回帰モデルが含まれている<ref>[http://www.mathworks.com/help/econ/autoregressive-model.html "Autoregressive modeling in MATLAB"]</ref>。
* [[MATLAB]]と[[GNU Octave|Octave]]: ''TSA toolbox'' に単一変数、複数変数、適応自己回帰モデルについてのいくつかの推定関数が含まれている<ref>[http://pub.ist.ac.at/~schloegl/matlab/tsa/ "Time Series Analysis toolbox for Matlab and Octave"]</ref>。

==''n'' 期先予測==
自己回帰
:<math> X_t = c + \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t \,</math>
のパラメーターが一度推定されてしまえば、この自己回帰は将来の任意の時点での予測に用いることが出来る。まず、''t'' をデータが使えない最初の時点とする。既知の値''X''<sub>''t-i''</sub> for ''i=''1, ..., ''p'' を自己回帰方程式に代入し、誤差項 <math>\varepsilon_t</math> をゼロと置く（なぜならば ''X''<sub>''t''</sub> をその条件つき期待値と一致させるように予測し、観測されない誤差項の期待値は0であるから）ことで予測ができる。自己回帰方程式の出力は最初のデータが観測されない時点についての予測となる。次に、''t'' をデータが使えない ''次の'' 時点とする。もう一度自己回帰方程式を予測を作るために使うことができる。ただし一つ異なる点がある。''X'' の今予測している時点より一期前の値は未知である。よってその期待値、つまり前の予測ステップでの予測値を代わりに用いる。この時、将来の時点において同じ手続きが用いられ、''p'' 回の予測の後に、全ての ''p'' 個の右辺の値が事前のステップによる予測値となるまで、予測方程式の右辺における予測値を用いる。

この方法で得られた予測値について四つの不確実性のソースがある。(1) 自己回帰モデルが正しいモデルかどうかという不確実性、(2) 自己回帰方程式の右辺においてラグ値として用いられる予測値の正しさについての不確実性、(3) 自己回帰係数の真の値についての不確実性、(4) 予測機関における誤差項 <math>\varepsilon_t \,</math> の値についての不確実性である。最後の三つは定量化可能で ''n'' ステップ後の予測についての[[信頼区間]]として与えられる。右辺の変数についての推定値が増えるため信頼区間は ''n'' が増えれば広くなる。

==予測の質の評価==
自己回帰モデルの予測性能は、[[クロス・バリデーション]]が行われるならば、推定の後に即座に評価できる。この方法においては、最初の方の利用可能なデータはパラメーターの推定の為に用いられ、データセットにおける後の方のデータはアウトオブサンプルのテストとして残しておく。他には、パラメーター推定が行われた後にしばらくしたあと、より多くのデータが利用可能になり予測性能を新しいデータを使うことで評価できる。

どちらのケースも、評価可能な予測性能には2つの側面がある。1期先予測の性能と ''n'' 期先予測の性能である。1期先予測の性能について、推定パラメーターは予測を行った期以前の全ての期における ''X'' の観測値と共に自己回帰方程式が用いられ、方程式の出力は1期先予測となる。この手続きはアウトオブサンプルの観測値についての予測を得るために用いられる。''n'' 期先予測の質を評価する為には、予測を得るために前の節での予測手続きが用いられる。

予測値のセットと対応する様々な期間の ''X'' の本当の値のセットが与えられたとして一般的な評価のテクニックは{{仮リンク|平均二乗予測誤差|en|mean squared prediction error}}を用いることである。他の尺度もまた用いられる。

ここで測定された予測の正しさをどのように解釈するのかという問題が持ち上がる。例えば平均二乗予測誤差が"高い"（悪い）もしくは"低い"（良い）とはどういう事なのだろうか。比較の上で二つのポイントがある。第一に他のモデルの仮定もしくは推定手法の下で推定された代替モデルの予測の正しさは比較目的に使用できる。第二にアウトオブサンプルの正確さの尺度は十分に前のデータを用いることが出来るならば、つまり最初の ''p'' 個のデータポイントを落として ''p'' 期以前のデータを使わないならば（パラメーター推定に用いられた）インサンプルのデータポイントでの同じ尺度と比較できる。モデルはインサンプルのデータポイントに出来るだけ適合するように特定化されて推定されるので、普通はアウトオブサンプルの予測性能はインサンプルの予測性能より悪い。しかし予測の質がアウトオブサンプルで（正確には定義できないが）"そう悪くない"のであれば、予測値は十分なパフォーマンスを見せていると言える。

== 統計モデル・生成モデルとしての表現 ==
上記のように自己回帰モデルは、決定論的/deterministicな線形変換に、確率的/probabilisticなブレ（確率項）が線形に追加されるモデルである。別の表現として、自己回帰モデルは統計モデル（生成モデル）で表すことができる<ref name=":0">亀岡 (2019) 深層生成モデルを用いた音声音響信号処理. http://www.kecl.ntt.co.jp/people/kameoka.hirokazu/publications/Kameoka2019SICE03_published.pdf</ref>。

自己回帰モデル AR(n) を考えるとし、因果関係を持つ時点 t から時刻 t-n までの値の組を

<math>X_{n}=(x_{t}, x_{t-1}, ..., x_{t-n})</math>

とする。X<sub>n</sub> の確率分布すなわち xt ~ xt-nの同時確率 p<sub>AR(n)</sub> は

<math>p_{AR(n)}(X_{n})=p_{AR(n)}(x_{t}, x_{t-1}, ..., x_{t-n})</math>

であり、条件付き確率の定義を用いる（{{仮リンク|確率の連鎖律|en|Chain rule (probability)|label=}}）と

<math>p_{AR(n)}(X_{n})=p(x_{t}|x_{t-1}, ..., x_{t-n})\cdot p(x_{t-1}, ..., x_{t-n})=...=\prod_{i=0}^n p(x_{t-i}|\bigcap_{j=i+1}^{n} x_{t-j})</math>

となる（条件付確率分布の積、因子分解の積）。モデルのn次自己回帰性（== 過去方向にn次までしか依存しない）より x<sub>t-n</sub> はそれ単体で分布が定まる。

p(x<sub>t-k</sub>|x<sub>t-k-1</sub>, x<sub>t-k-2</sub>, ..., x<sub>t-n</sub>)を考えると（線形変換・ガウス確率項）ARモデルは確率項がガウス分布に従い、その平均値統計量は決定論的な線形変換で決まるため、

<math>p(x_{t-k}|x_{t-k-1}, ..., x_{t-n})=N(x_{t-k}|w_{t-k-1}\cdot x_{t-k-1}, \sigma )</math>

であるといえる。x<sub>t-k-1</sub>より過去の系列の情報はすべてx<sub>t-k-1</sub>の実現値として集約されている。

以上をまとめると、n次自己回帰モデルは統計モデル/生成モデルとして以下のように定式化できる。

<math>X_{n}:=(x_{t}, x_{t-1}, ..., x_{t-n})</math>

<math>x_{t-k} \thicksim p(x_{t-k}|x_{t-k-1}, ..., x_{t-n})=N(x_{t-k}|w_{t-k-1}\cdot x_{t-k-1}, \sigma )</math>・・・(1)

<math>X_{n}\thicksim p_{AR(n)}(X_{n})
=(\prod_{i=0}^{n-1} p(x_{t-i}|\bigcap_{j=i+1}^{n} x_{t-j}))\cdot p(x_{t-n})
=(\prod_{i=0}^{n-1}N(x_{t-i}|w_{t-i-1}\cdot x_{t-i-1}, \sigma))\cdot p(x_{t-n})</math>

確率分布が計算可能なため、データ（標本）が与えられた際の母集団パラメータを[[最尤推定]]を用いて推定することができる。また生成モデルであるため、モデルに従う系列の生成（例: 現時刻までのデータに基づいた将来系列の生成）が可能である。もし音声を時系列とみなせば音声合成をおこなうことが可能になる。

== 非線形自己回帰生成モデル ==
古典的な自己回帰モデルは系列要素間の関係を'''線形'''と仮定してきたが、近年では'''非線形'''自己回帰モデルも提唱されている。[[人工ニューラルネットワーク]]と[[深層学習]]の発達により発達した自己回帰生成ネットワーク（Autoregressive Generative Networks）がその代表例である。

自己回帰モデルを生成モデルとして表現したとき、x<sub>t</sub> はそれ以前の値（x<sub>t-1</sub>, x<sub>t-2</sub> ,...）で条件づけられた確率分布からサンプリングされる。線形の自己回帰モデルではガウス分布の平均値が前要素の線形変換になるとモデル化するが、この条件は緩和することができる。上記の（１）式、すなわち過去値に条件づけられた確率分布関数を非線形関数を含む任意の関数とすることでこれが達成できる。

確率分布を[[人工ニューラルネットワーク]]によって表現することで非線形性を導入し[[深層学習]]によって実データに基づく分布の推定/学習（最尤推定など）をおこなったものが自己回帰生成ネットワークである<ref name=":0" />。[[DeepMind]]社が開発したWaveNet<ref>The model is fully probabilistic and autoregressive, with the predictive distribution for each audio sample conditioned on all previous ones...

Aaron van den Oord, et al.. (2016) WaveNet: A Generative Model for Raw Audio</ref>はAutoregressive Generative Networksの代表例であり、音声波形を系列とみなして自己回帰モデル化・学習することにより、人の声と区別がつかない音声の合成に成功している。

=== バリエーション ===
非線形自己回帰生成モデルは制約を緩めた分、いくつかの変種（モデル上の違い）がある。

==== 確率分布 ====

* カテゴリカル分布（[[ソフトマックス関数]])

==== サンプリング ====

* ランダムサンプリング
* 最大確率（ArgMax）
* ビームサーチ

=== 統計的推論（学習）・最適化 ===
非線形自己回帰生成モデルの難点の1つは、非線形性からくるパラメータ推定の難しさにある。統計的推論には[[最尤推定]]がしばしば用いられるが、古典的な（線形）ARモデルと比較してパラメータ推定の難易度が高い。

==== teacher forcing ====
teacher forcing は自己回帰入力に教師信号を用いる、自己回帰モデル学習技法の1つである<ref>"to replace the actual output <math>y_k(t)</math> of a unit by the teacher signal <math>d_k(t)</math> in subsequent computation of the behavior of the network, whenever such a value exists. We call this technique 'teacher forcing.' " Williams & Zipser. (1989). ''A Learning Algorithm for Continually Running Fully Recurrent Neural Networks''. doi: 10.1162/neco.1989.1.2.270</ref>。自己回帰モデルは系列長が長くなるほど誤差を蓄積する特性をもつ。teacher forcingは学習時に教師信号を自己回帰入力することで、誤差のない（理想的な）入力に基づいた学習を可能にする。

===== exposure bias =====
teacher forcingで学習したモデルに関して、推論時に与えられる自己回帰入力は文字通り「自己回帰」であり、学習時に得られるようなノイズの無い理想信号とは限らない。ゆえに学習が不十分なモデルでは推論時の自己回帰入力が学習時と乖離してしまうため、その振る舞いは予期できないものになる。この問題を '''exposure bias''' という。

==脚注==
{{Reflist | refs=

<ref name=Brockwell>{{cite journal| first1=Peter J.|last1=Brockwell |first2= Rainer|last2= Dahlhaus| first3=A. Alexandre|last3= Trindade
 |journal=Statistica Sinica | volume= 15 |year=2005 |pages=197–213
 |title=Modified Burg Algorithms for Multivariate Subset Autoregression
 |url=http://www3.stat.sinica.edu.tw/statistica/oldpdf/A15n112.pdf}}</ref>

<ref name=Bos>{{Cite journal | last1 = Bos | first1 = R. | last2 = De Waele | first2 = S. | last3 = Broersen | first3 = P. M. T. | doi = 10.1109/TIM.2002.808031 | title = Autoregressive spectral estimation by application of the burg algorithm to irregularly sampled data | journal = IEEE Transactions on Instrumentation and Measurement | volume = 51 | issue = 6 | pages = 1289 | year = 2002 | pmid =  | pmc = }}</ref>

<ref name=Burg>{{Citation
|last = Burg
|first = John P.
|year = 1968
|contribution = A new analysis technique for time series data
|title = Modern Spectrum Analysis
|editor = D. G. Childers
|series = NATO Advanced Study Institute of Signal Processing with emphasis on Underwater Acoustics
|publisher = IEEE Press
|location = New York}}</ref>

<ref name=Burg1>{{Citation
|last = Burg
|first = John P. 
|year = 1967
|title = Maximum Entropy Spectral Analysis
|journal = Proceedings of the 37th Meeting of the Society of Exploration Geophysicists
|publisher = Oklahoma City
|location = Oklahoma}}</ref>

}}

==参考文献==
* {{cite book |last=Mills |first=Terence C. |year=1990 |title=Time Series Techniques for Economists |location= |publisher=Cambridge University Press |isbn= }}
* {{cite book |last=Percival |first=Donald B. |first2=Andrew T. |last2=Walden |year=1993 |title=Spectral Analysis for Physical Applications |location= |publisher=Cambridge University Press |isbn= }}
* {{cite book |last=Pandit |first=Sudhakar M. |last2=Wu |first2=Shien-Ming |year=1983 |title=Time Series and System Analysis with Applications |location= |publisher=John Wiley & Sons }}
* {{Citation 
|last=Hamilton
|first=James D.
|authorlink = ジェームズ・ハミルトン (計量経済学者)
|year=1994 
|title=Time Series Analysis
|isbn = 0691042896
|publisher = [[Princeton University Press]]
|ref=Hamilton1994}}

==関連項目==
* [[移動平均モデル]]
* [[予測分析]]
* [[線形予測符号|線形予測コーディング]]
* [[共鳴]]

==外部リンク==
*[http://paulbourke.net/miscellaneous/ar/ AutoRegression Analysis (AR)] by Paul Bourke
*{{YouTube|id=b8yslhlIsA0&list=PLD15D38DC7AA3B737&index=8#t=34m25s|title=Econometrics lecture (topic: Autoregressive models)}} by {{仮リンク|マーク・ソーマ|en|Mark Thoma}}

{{統計学}}
{{デフォルトソート:しこかいきもてる}}
[[Category:統計学]]
[[Category:時系列分析]]
[[Category:信号処理]]