自己回帰モデル

自己回帰モデル（じこかいきモデル、英: 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 に恒久的に影響を与えるため、任意の与えられた X_t の値は過去に起こったショック全てから影響を受ける。これは自己回帰方程式

\phi (B)X_{t}=\varepsilon _{t}\,

（ここで定数項は変数が平均からの逸脱として測定されると仮定することで無視できる）が以下のように書き直せることからもまた分かる。

X_{t}={\frac {1}{\phi (B)}}\varepsilon _{t}\,.

右辺における多項式の除算が可能なのであれば、 $\varepsilon _{t}$ に適用される後退オペレーターによる多項式は無限次元のオーダーを持つ。つまり、 $\varepsilon _{t}$ のラグ値が方程式の右辺において無限個現れる。

特性多項式

AR(p) 過程の自己相関関数は以下のように表すことが出来る^[1]。

\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(p) 過程のグラフ

"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) はホワイトノイズに対応する。

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

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

例: AR(1) 過程

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 ^{|n|}.

自己共分散関数は $\tau =-1/\ln(\varphi )$ の減衰時間（または時定数）で減衰していくことが分かる[これを見るために $B_{n}=K\phi ^{|n|}$ と書く。ここで $K$ は $n$ に依存しない。この時、 $\phi ^{|n|}=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 ^{|t|}

これにより、コーシー分布（ローレンツ型プロファイル）のスペクトル密度が得られる。

\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 ^{N}X_{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) モデルの性質を理解するために同様の形式に変換することが時として有用になる。この形式において 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パラメーターの計算

ARモデルの係数の推定には多数の方法があり、例えば最小二乗法の手続きや、もしくは（ユール–ウォーカー方程式（英: Yule–Walker equation）を通した）モーメント法がある。

AR(p) モデルは以下の方程式で与えられる。

X_{t}=\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}.\,

この方程式はパラメーター $\varphi _{i}$ i = 1, ..., p に基いている。これらのパラメーターと過程の共分散関数の間には直接的な対応が存在し、その対応は（共分散から得られる）自己相関関数からパラメーターを決定する為に裏返すことができる。これはユール–ウォーカー方程式を用いて行われる。

ユール–ウォーカー方程式

ユール–ウォーカー方程式は、ウドニー・ユール（英語版）とギルバート・ウォーカーにちなんで名づけられたもので^[2]^[3]、以下の方程式からなる^[4]。

\gamma _{m}=\sum _{k=1}^{p}\varphi _{k}\gamma _{m-k}+\sigma _{\varepsilon }^{2}\delta _{m,0},

ここで m = 0, ..., p であり、p + 1 個の方程式からなる。さらに $\gamma _{m}$ は X_t の自己共分散関数、 $\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 個の要素で決定する。完全な自己相関関数はこの時、再帰的な計算によって導出できる^[4]^[5]。

\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パラメーターの推定

上の方程式（ユール–ウォーカー方程式）は、理論的な共分散を推定値に置き換えることで、 AR(p) モデルのパラメーターを推定する為にいくつかの方法を提供する^[6]。下記のような方法が考えられる。

自己共分散もしくは自己相関の推定。便利な推定法を用いて自己共分散もしくは自己相関の項のそれぞれを分割して推定したものとする。推定の方法は多様であり、どれを選択するかは推定のスキームが持つ性質に影響を与える。例えば、ある方法では分散の負の推定量が生じうる。
X_t の予測値を同じ系列の過去の p 個の値として基礎づける最小二乗予測問題を構築する上での最小二乗回帰問題としての定式化。これは前方予測スキームとして考えられる。この問題についての正規方程式（英語版）は同じラグで現れる自己共分散を少し違った推定値で置き換えたユール–ウォーカー方程式の行列形式の近似と対応するように見える。
最小二乗予測問題の拡張形式としての定式化。ここで二つの予測方程式のセットを一つの推定スキームと単一の正規方程式に結合する。一つのセットは前方予測方程式のセットとなっており、もう片方は対応する後方予測方程式のセットとなっている。これはARモデルの後退表現と関連している。

X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}^{*}\,.

ここで予測値 X_t は同じ系列の p 個の将来の値に基づいている。このARパラメーターの推定方法はジョン・バーグ（John P. Burg）^[7] によるものでバーグの方法（英: the Burg method）と呼ばれる^[8]。バーグや後続の研修者はこの特別な推定値を"最大エントロピー推定量"と呼ぶが^[9]、この背後にある理論は推定パラメータ―のどのようなセットについても適用できる。前進予測方程式のみを用いた推定スキームと比べると、異なる自己共分散の推定値が得られ、推定量は異なる安定性の性質を持つ。バーグ推定量は特に最大エントロピースペクトル推定（英語版）と関連している^[10]。

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

スペクトル

ノイズ分散が $\mathrm {Var} (Z_{t})=\sigma _{Z}^{2}$ である AR(p) 過程のパワースペクトル密度は以下のようになる^[11]^[5]。

S(f)={\frac {\sigma _{Z}^{2}}{|1-\sum _{k=1}^{p}\varphi _{k}e^{-2\pi ikf}|^{2}}}.

AR(0)

ホワイトノイズ (AR(0)) については以下のようになる。

S(f)=\sigma _{Z}^{2}.

AR(1)

AR(1) については以下のようになる。

S(f)={\frac {\sigma _{Z}^{2}}{|1-\varphi _{1}e^{-2\pi if}|^{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)

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\leq \varphi _{2}\leq 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 関数が含まれている^[12]。
MATLAB, Econometric Toolbox^[13] と System Identification Toolbox^[14] に自己回帰モデルが含まれている^[15]。
MATLABとOctave: TSA toolbox に単一変数、複数変数、適応自己回帰モデルについてのいくつかの推定関数が含まれている^[16]。

n 期先予測

自己回帰

X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}\,

のパラメーターが一度推定されてしまえば、この自己回帰は将来の任意の時点での予測に用いることが出来る。まず、t をデータが使えない最初の時点とする。既知の値X_t-i for i=1, ..., p を自己回帰方程式に代入し、誤差項 $\varepsilon _{t}$ をゼロと置く（なぜならば X_t をその条件つき期待値と一致させるように予測し、観測されない誤差項の期待値は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なブレ（確率項）が線形に追加されるモデルである。別の表現として、自己回帰モデルは統計モデル（生成モデル）で表すことができる^[17]。

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

$X_{n}=(x_{t},x_{t-1},...,x_{t-n})$

とする。X_n の確率分布すなわち xt ~ xt-nの同時確率 p_AR(n) は

$p_{AR(n)}(X_{n})=p_{AR(n)}(x_{t},x_{t-1},...,x_{t-n})$

であり、条件付き確率の定義を用いる（確率の連鎖律（英語版））と

$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})$

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

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

$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 )$

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

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

$X_{n}:=(x_{t},x_{t-1},...,x_{t-n})$

$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 )$ ・・・(1)

$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})$

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

非線形自己回帰生成モデル

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

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

確率分布を人工ニューラルネットワークによって表現することで非線形性を導入し深層学習によって実データに基づく分布の推定/学習（最尤推定など）をおこなったものが自己回帰生成ネットワークである^[17]。DeepMind社が開発したWaveNet^[18]はAutoregressive Generative Networksの代表例であり、音声波形を系列とみなして自己回帰モデル化・学習することにより、人の声と区別がつかない音声の合成に成功している。

バリエーション

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

確率分布

カテゴリカル分布（softmax)

サンプリング

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

統計的推論（学習）・最適化

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

脚注

^ Zetterberg, Lars H. (1969), “Estimation of parameters for a linear difference equation with application to EEG analysis”, Mathematical Biosciences 5 (3): 227--275, doi:10.1016/0025-5564(69)90044-3, ISSN 0025-5564
^ Yule, G. Udny (1927), “On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers”, Philosophical Transactions of the Royal Society of London, Ser. A 226: 267–298
^ Walker, Gilbert (1931), “On Periodicity in Series of Related Terms”, Proceedings of the Royal Society of London, Ser. A 131: 518–532
^ ^a ^b Hamilton & (1994), p. 59
^ ^a ^b Von Storch, H.; F. W Zwiers (2001). Statistical analysis in climate research. Cambridge Univ Pr. ISBN 0-521-01230-9 ^{[要ページ番号]}
^ Hamilton & (1994), Chapter 3 and 5
^ Burg, John P. (1968), “A new analysis technique for time series data”, in D. G. Childers, Modern Spectrum Analysis, NATO Advanced Study Institute of Signal Processing with emphasis on Underwater Acoustics, New York: IEEE Press
^ Brockwell, Peter J.; Dahlhaus, Rainer; Trindade, A. Alexandre (2005). “Modified Burg Algorithms for Multivariate Subset Autoregression”. Statistica Sinica 15: 197–213.
^ Burg, John P. (1967), “Maximum Entropy Spectral Analysis”, Proceedings of the 37th Meeting of the Society of Exploration Geophysicists (Oklahoma: Oklahoma City)
^ Bos, R.; De Waele, S.; Broersen, P. M. T. (2002). “Autoregressive spectral estimation by application of the burg algorithm to irregularly sampled data”. IEEE Transactions on Instrumentation and Measurement 51 (6): 1289. doi:10.1109/TIM.2002.808031.
^ Hamilton & (1994), p. 155
^ "Fit Autoregressive Models to Time Series" (in R)
^ Econometrics Toolbox Overview
^ System Identification Toolbox overview
^ "Autoregressive modeling in MATLAB"
^ "Time Series Analysis toolbox for Matlab and Octave"
^ ^a ^b 亀岡 (2019) 深層生成モデルを用いた音声音響信号処理. http://www.kecl.ntt.co.jp/people/kameoka.hirokazu/publications/Kameoka2019SICE03_published.pdf
^ 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

参考文献

Mills, Terence C. (1990). Time Series Techniques for Economists. Cambridge University Press
Percival, Donald B.; Walden, Andrew T. (1993). Spectral Analysis for Physical Applications. Cambridge University Press
Pandit, Sudhakar M.; Wu, Shien-Ming (1983). Time Series and System Analysis with Applications. John Wiley & Sons
Hamilton, James D. (1994), Time Series Analysis, Princeton University Press, ISBN 0691042896

外部リンク

AutoRegression Analysis (AR) by Paul Bourke
Econometrics lecture (topic: Autoregressive models) - YouTube by マーク・ソーマ（英語版）

[1] Zetterberg, Lars H. (1969), “Estimation of parameters for a linear difference equation with application to EEG analysis”, Mathematical Biosciences 5 (3): 227--275, doi:10.1016/0025-5564(69)90044-3, ISSN 0025-5564

[2] Yule, G. Udny (1927), “On a Method of Investigating Periodicities in Disturbed Series, with Special Reference to Wolfer's Sunspot Numbers”, Philosophical Transactions of the Royal Society of London, Ser. A 226: 267–298

[3] Walker, Gilbert (1931), “On Periodicity in Series of Related Terms”, Proceedings of the Royal Society of London, Ser. A 131: 518–532

[Hamilton59-4] Hamilton & (1994), p. 59

[Storch-5] Von Storch, H.; F. W Zwiers (2001). Statistical analysis in climate research. Cambridge Univ Pr. ISBN 0-521-01230-9 ^{[要ページ番号]}

[6] Hamilton & (1994), Chapter 3 and 5

[Burg-7] Burg, John P. (1968), “A new analysis technique for time series data”, in D. G. Childers, Modern Spectrum Analysis, NATO Advanced Study Institute of Signal Processing with emphasis on Underwater Acoustics, New York: IEEE Press

[Brockwell-8] Brockwell, Peter J.; Dahlhaus, Rainer; Trindade, A. Alexandre (2005). “Modified Burg Algorithms for Multivariate Subset Autoregression”. Statistica Sinica 15: 197–213.

[Burg1-9] Burg, John P. (1967), “Maximum Entropy Spectral Analysis”, Proceedings of the 37th Meeting of the Society of Exploration Geophysicists (Oklahoma: Oklahoma City)

[Bos-10] Bos, R.; De Waele, S.; Broersen, P. M. T. (2002). “Autoregressive spectral estimation by application of the burg algorithm to irregularly sampled data”. IEEE Transactions on Instrumentation and Measurement 51 (6): 1289. doi:10.1109/TIM.2002.808031.

[11] Hamilton & (1994), p. 155

[12] "Fit Autoregressive Models to Time Series" (in R)

[13] Econometrics Toolbox Overview

[14] System Identification Toolbox overview

[15] "Autoregressive modeling in MATLAB"

[16] "Time Series Analysis toolbox for Matlab and Octave"

[:0-17] 亀岡 (2019) 深層生成モデルを用いた音声音響信号処理. http://www.kecl.ntt.co.jp/people/kameoka.hirokazu/publications/Kameoka2019SICE03_published.pdf

[18] 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

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

定義