Zenn
Open8

ARIMAを復習するのであります！（時系列解析）

統計学
計量経済学
arima
時系列解析
畳屋民也畳屋民也

AR(1) 過程と期待値・分散・自己共分散

設定

y_{t+1} = c + \phi y_t + \varepsilon_{t+1}

ただし、|\phi|<1 であり、\varepsilon_t は下記の条件を満たす分散が\sigma^2のホワイトノイズとする：

E[\varepsilon_t] = 0, \quad E[\varepsilon_t^2] = \sigma^2\\ E[\varepsilon_{t_1} \varepsilon_{t_2}] = 0 \quad (t_1 \ne t_2)

統計量（結果）

以下、定常状態を仮定する。

期待値

E[y_t] = \frac{c}{1-\phi}

分散

V[y_t] = \frac{\sigma^2}{1 - \phi^2}

自己共分散

Cov[y_t, y_{t+k}] = \phi^k \frac{\sigma^2}{1 - \phi^2}

自己相関係数

\gamma_k = \frac{Cov[y_t, y_{t+k}]}{\sqrt{V[y_t]} \sqrt{V[y_{t+k}]}} = \phi^k

統計量（導出）

期待値

\begin{aligned} E[y_t] &= E[c + \phi y_{t-1} + \varepsilon_t] = c + \phi E[y_{t-1}] + 0 \\ &= c + \phi (c + \phi E[y_{t-2}]) = c(1 + \phi) + \phi^2 E[y_{t-2}]\\ &= ... \\ &= c(1 + \phi + \phi^2 + ... + \phi^{n-1}) + \phi^n E[y_{t-n}] = c \frac{1-\phi^n}{1-\phi} + \phi^n E[y_{t-n}] \\ &\xrightarrow{n\to\infty} \frac{c}{1-\phi} \quad (\because |\phi|<1 \Rightarrow \lim_{n\to\infty} \phi^n = 0) \end{aligned}

分散

E[y_t] = \mu とおくと、c = \mu (1 - \phi) より、

y_{t+1} - \mu = \phi (y_t - \mu) + \varepsilon_{t+1}

したがって、

\begin{aligned} V[y_t] &= E[(y_t - \mu)^2] \\ &= E[\{ \phi (y_{t-1} - \mu) + \varepsilon_t \}^2] = \phi^2 E[(y_{t-1} - \mu)^2] + \sigma^2 \quad (\because E[\varepsilon_t^2]=\sigma^2)\\ &= \phi^4 E[(y_{t-2} - \mu)^2] + \phi^2 \sigma^2 + \sigma^2\\ &= ...\\ &= \phi^{2n} E[(y_{t-n} - \mu)^2] + \sigma^2 (1 + \phi^2 + \phi^4 + ... + \phi^{2(n-1)})\\ &= \phi^{2n} E[(y_{t-n} - \mu)^2] + \sigma^2 \frac{1 - \phi^{2n}}{1 - \phi^2}\\ &\xrightarrow{n\to\infty} \frac{\sigma^2}{1 - \phi^2} \end{aligned}

自己共分散

\begin{aligned} Cov[y_t, y_{t+k}] &= E[(y_t-\mu)(y_{t+k}-\mu)]\\ &= E[(y_t-\mu)\{ \phi (y_{t+k-1} - \mu) + \varepsilon_{t+k} \}] = \phi E[(y_t-\mu)(y_{t+k-1}-\mu)] + 0\\ &= \phi^2 E[(y_t-\mu)(y_{t+k-2}-\mu)]\\ &= ... \\ &= \phi^k E[(y_t-\mu)(y_t-\mu)] = \phi^k V[y_t]\\ &= \phi^k \frac{\sigma^2}{1 - \phi^2} \end{aligned}

自己相関係数

\begin{aligned} \gamma_k &= \frac{Cov[y_t, y_{t+k}]}{\sqrt{V[y_t]} \sqrt{V[y_{t+k}]}} \\ &= \frac{\phi^k V[y_t]}{V[y_t]}\\ &= \phi^k \end{aligned}
畳屋民也畳屋民也

AR(1) 過程の予測値

結果

予測期待値

\begin{aligned} E[y_{t+k}\vert y_t ] &= \mu + \phi^k (y_t - \mu) \\ &\xrightarrow{k\to\infty} \mu \end{aligned}

ただし、\mu = c / (1 - \phi)。

予測分散

\begin{aligned} V[y_{t+k}\vert y_t] &= \sigma^2 \frac{1-\phi^{2k}}{1-\phi^2}\\ &\xrightarrow{k\to\infty} \frac{\sigma^2}{1-\phi^2} \end{aligned}

導出

予測期待値

\begin{aligned} E[y_{t+k}\vert y_t ] &= E[c + \phi y_{t+k-1} + \varepsilon_{t+k}\vert y_t] = c + \phi E[y_{t+k-1}\vert y_t] \quad (\because E[\varepsilon_{t+k} \vert y_t] = 0)\\ &= c + \phi c + \phi^2 E[y_{t+k-2}\vert y_t]\\ &= ... \\ &= c(1 + \phi + \phi^2 + ... \phi^{k-1}) + \phi^k E[y_t\vert y_t]\\ &= c \frac{1-\phi^k}{1-\phi} + \phi^k y_t\\ &= \mu + \phi^k (y_t - \mu) \quad \left(\because \frac{c}{1 - \phi} = \mu \right) \end{aligned}

予測分散

\begin{aligned} V[y_{t+k}\vert y_t] &= E[(y_{t+k} - E[y_{t+k}\vert y_t])^2 \vert y_t]\\ &= E[\{ y_{t+k} - \mu - \phi^k (y_t - \mu) \}^2 \vert y_t]\\ &=E[\{ \phi (y_{t+k-1} - \mu) + \varepsilon_{t+k} - \phi^k (y_t - \mu) \}^2 \vert y_t]\\ &= E[ \{ \phi (y_{t+k-1} - \mu - \phi^{k-1}(y_t - \mu) ) + \varepsilon_{t+k} \}^2 \vert y_t]\\ &= E[\{ \phi (y_{t+k-1} - E[y_{t+k-1}\vert y_t]) \}^2 \vert y_t] + E[\varepsilon_{t+k}^2\vert y_t]\\ &= \phi^2 V[y_{t+k-1}\vert y_t] + \sigma^2\\ &= \phi^4 V[y_{t+k-2}\vert y_t] + \phi^2 \sigma^2 + \sigma^2\\ &= ... \\ &= \phi^{2k} V[y_t \vert y_t] + \sigma^2 (1 + \phi^2 + \phi^4 + ... + \phi^{2(k-1)})\\ &=\sigma^2 \frac{1 - \phi^{2k}}{1-\phi^2} \quad (\because V[y_t \vert y_t] = 0) \end{aligned}
畳屋民也畳屋民也

図で見てわかる AR(1) の予測値の挙動

https://zenn.dev/link/comments/c1e650535540b8

前回求めたAR(1)過程の予測期待値・予測分散が何を意味しているかを視覚的に理解するために、Pythonでシミュレーションを行なってみた。

以下では、t=0における初期値y_0=-1.0を出発点として、パラメータ c=1.2, \phi=0.8, \sigma^2=1.0 の AR(1) 過程よる予測値をシミュレートしたものである。

import numpy as np
import matplotlib.pyplot as plt

t = np.arange(0, 51, 1)

y0 = -1.0

c = 1.2
phi = 0.8
sigma = 1.0

fig, ax = plt.subplots(1,1)
for i in range(0,10):
  list_ar1 = []
  y_prev = y0
  list_ar1.append(y0)
  eps = np.random.normal(0, sigma, 50)
  for i in range(0,50):
    yt = c + phi * y_prev + eps[i]
    list_ar1.append(yt)
    y_prev = yt

  ax.plot(t, list_ar1)

mu = c / (1 - phi)
muk = mu + phi ** t * (y0 - mu)
vk = sigma**2 * (1 - phi ** (2 * t)) / (1 - phi**2)
upper_bound = muk + 1.96 * np.sqrt(vk)
lower_bound = muk - 1.96 * np.sqrt(vk)

ax.plot(t, muk, ls="-", c="gray")
ax.plot(upper_bound, ls="--", c="gray")
ax.plot(lower_bound, ls="--", c="gray")

ax.set_ylabel("$y_t$")
ax.set_xlabel("$t$")

色のついた実線は、異なる乱数を用いて描いた10本のシミュレーション軌道である。

灰色の実線は予測期待値

\begin{aligned} E[y_t\vert y_0 ] &= \mu + \phi^t (y_0 - \mu) \\ &\xrightarrow{t\to\infty} \mu \end{aligned}

である（ただし、\mu = c / (1 - \phi)）。

t=0でy_0から出発し、時間の経過と共に値 \mu = 6.0に収束しているのがわかる。

灰色の点線は、予測分散

\begin{aligned} V[y_t \vert y_0] &= \sigma^2 \frac{1-\phi^{2t}}{1-\phi^2}\\ &\xrightarrow{t\to\infty} \frac{\sigma^2}{1-\phi^2} \end{aligned}

に基づいて計算した95%信頼区間である。
予測分散はt=1 では\sigma^2=1.0 だが、時間の経過と共に \sigma^2 / (1 - \phi^2) = 2.777...に収束する。
それに合わせて、信頼区間も時間と共に幅が広がりながら最終的には一定の長さに落ち着くことが見てとれる。

シミュレーション軌道は、予測期待値を中心に、おおよそこの点線の範囲内で変動している。

畳屋民也畳屋民也

AR(1) 過程のパラメータ推定

設定

以下のようなモデルを考える：

y_{t+1} = c + \phi y_t + \varepsilon_{t+1}

ただし、\varepsilon_t は下記の条件を満たす分散が\sigma^2のホワイトノイズとする：

E[\varepsilon_t] = 0, \quad E[\varepsilon_t^2] = \sigma^2\\ E[\varepsilon_{t_1} \varepsilon_{t_2}] = 0 \quad (t_1 \ne t_2)

パラメータ推定値（結果）

T+1 個の観測値 \{ y_1, y_2, ..., y_T, y_{T+1}\}が与えられたとき、以下はそれぞれパラメータ \boldsymbol{\theta} = \{c, \phi, \sigma^2\} の一致推定量 になる：

\begin{aligned} \hat{c} & = \bar{y}_{T+1} - \hat{\phi} \bar{y}_T \\ \hat{\phi} &= \frac{\sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1})(y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \\ s^2 &= \frac{1}{T} \sum_{t=1}^T (y_{t+1} - \hat{c} - \hat{\phi} y_t)^2 \end{aligned}

ここで、

\bar{y}_{T+1} = \frac{1}{T} \sum_{t=1}^T y_{t+1}, \quad \bar{y}_{T} = \frac{1}{T} \sum_{t=1}^T y_{t}

とした。

なお、いずれも、不偏推定量ではない。

パラメータ推定方法

係数 c, \phi

最小二乗法による推定を行う。

残差平方和を以下のように定義する：

RSS(c, \phi) = \sum_{t=1}^T (y_{t+1} - c - \phi y_t)^2.

これを最小化するような c, \phi を求めればよい。

そのためには、残差平方和を偏微分して値が0になるc=\hat{c}, \phi=\hat{\phi} の値を探せば良い：

\frac{\partial RSS(c, \phi)}{\partial c} \bigg |_{\substack{c=\hat{c}\\ \phi=\hat{\phi}}} = 0, \quad \frac{\partial RSS(c, \phi)}{\partial \phi} \bigg |_{\substack{c=\hat{c}\\ \phi=\hat{\phi}}} = 0.

これを解くと、

\begin{aligned} \hat{c} & = \bar{y}_{T+1} - \hat{\phi} \bar{y}_T \\ \hat{\phi} &= \frac{\sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1})(y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \end{aligned}

のように求まる。

誤差項分散 \sigma^2

いったん天下りに以下のように定義する：

s^2 = \frac{1}{T} \sum_{t=1}^T (y_{t+1} - \hat{c} - \hat{\phi} y_t)^2

この時、この推定量は \sigma^2 の一致推定量になる：

\plim_{T\to \infty} s^2 = \sigma^2.

補足： 最尤推定による方法

以下、誤差項が正規分布であると仮定を置く：

\varepsilon_t \sim N(0, \sigma^2).

このとき、\{c, \phi\} の推定量 \{\hat{c}, \hat{\phi}\} が上述の最小二乗法によるものと一致するほか、
\sigma^2 の推定量 s^2 についても自然に求まることを示す：

係数 c, \phi

Y_tで条件づけられたときのY_{t+1} の尤度関数は以下のように表せる：

f_{Y_{t+1} \vert Y_t }(y_{t+1} \vert y_t ), \boldsymbol{\theta}) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left \{ - \frac{(y_{t+1} - c - \phi y_t )^2}{2\sigma^2} \right\}

したがって、Y_1 で条件づけられた \{ Y_1, Y_2, ..., Y_{T+1}, Y_T\} の条件付き尤度関数は、以下のように表される：

\begin{aligned} f_{Y_{T+1}, Y_T, ..., Y_3, Y_2 \vert Y_1} (y_{T+1}, y_T, ..., y_3, y_2 \vert y_1; \boldsymbol{\theta}) &= \prod_{t=1}^T f_{Y_{t+1} \vert Y_t}(y_{t+1} \vert y_t, \boldsymbol{\theta}) \end{aligned}

対数尤度は、

\begin{aligned} \log f_{Y_{T+1}, Y_T, ..., Y_3, Y_2 \vert Y_1} (y_{T+1}, y_T, ..., y_3, y_2 \vert y_1; \boldsymbol{\theta}) \\ = - \frac{T}{2} \log (2\pi) - \frac{T}{2} \log \sigma^2 - \sum_{t=1}^T \frac{(y_{t+1} - c - \phi y_t)^2}{2\sigma^2} \tag{1} \end{aligned}

この対数尤度関数が最大になる \{c, \phi\} を求めることは、残差平方和

\sum_{t=1}^T (y_{t+1} - c - \phi y_t)^2

を最小化する \{c, \phi\} を求めることと同値である。

したがって、先ほどと同様に、

\begin{aligned} \hat{c} & = \bar{y}_{T+1} - \hat{\phi} \bar{y}_T \\ \hat{\phi} &= \frac{\sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1})(y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \end{aligned}

が得られる。

誤差項分散 \sigma^2

\sigma^2=\eta とおき、対数尤度関数を最大化する \eta=s^2 を求める。
そのためには、式(1) を \eta で偏微分すればよい。

すると、以下のように求まる：

s^2 = \frac{1}{T} \sum_{t=1}^T (y_{t+1} - \hat{c} - \hat{\phi} y_t)^2
畳屋民也畳屋民也

AR(1) 過程において係数の最小二乗法による推定値が不偏推定量にならないことの説明

以下の記事でAR(1) 過程のパラメータ \{c, \phi\} の最小二乗法による推定方法を記述した：
https://zenn.dev/link/comments/8a28e47602c1cf

ここでは、この時の係数\phi の最小二乗推定量

\hat{\phi} = \frac{\sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1})(y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2}

が不偏推定量にならないことを示す：

E[\hat{\phi}] \ne \phi

（時系列でない）単回帰ではなぜ最小二乗推定が不偏だったか？

まずシンプルな単回帰に立ち戻り、単回帰ではなぜ最小二乗推定による結果が不偏推定量になったのかを復習する。

モデルは、

y_i = \beta_0 + \beta_1 x_i + \varepsilon_i

ただし、\varepsilon_i は以下のように分散が\sigma^2 の独立同時なホワイトノイズである：

E[\varepsilon_i]=0, \quad E[\varepsilon_{i}\varepsilon_{j}] = \left\{ \begin{array}{ll} \sigma^2 & i=j \\ 0 & i \ne j \end{array} \right. .

また、\{x_i\} は独立同時である。

N個の観測値\{(x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\}が与えられたとき、係数\beta_1 の最小二乗推定量 \hat{\beta}_1 は以下のようになる：

\begin{aligned} \hat{\beta}_1 &= \frac{\sum_{i=1}^N (y_i - \bar{y})(x_i - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2}\\ &= \beta_1 + \frac{\sum_{i=1}^N \varepsilon_i (x_i - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2} \end{aligned}

ただし、

\bar{x} = \frac{1}{N} \sum_{i=1}^N x_i, \quad \bar{y} = \frac{1}{N} \sum_{i=1}^N y_i

このとき \hat{\beta}_1 の期待値を求めると、

\begin{aligned} E[\hat{\beta}_1] &= \beta_1 + E \left[ \frac{\sum_{i=1}^N \varepsilon_i (x_i - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2} \right]\\ &= \beta_1 + \sum_{i=1}^N E \left[ \frac{ (x_i - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2} \varepsilon_i \right]\\ &= \beta_1 + \sum_{i=1}^N \frac{(x_i - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2} E \left[ \varepsilon_i \right] \\ &= \beta_1 + 0 \end{aligned}

のようになり、\hat{\beta}_1 が \beta_1 の不偏推定量になっていることがわかる。

AR(1) モデルの場合

AR(1) モデルでも、\hat{\phi} を以下のように表すことができる：

\hat{\phi} = \phi + \frac{\sum_{t=1}^T \varepsilon_{t+1} (y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2}.

しかし、前述の単回帰の時と異なり、上式の右辺第二項の期待値が0にならない：

E\left[ \frac{\sum_{t=1}^T \varepsilon_{t+1} (y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \right] \ne 0.

なぜかというと、\varepsilon_t と y_{t+k} \, (k\ge 0) が独立でないからである（補足参照）。

実際、

E\left[ \frac{\sum_{t=1}^T \varepsilon_{t+1} (y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \right] = \sum_{t=1}^T E\left[ \frac{ (y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \varepsilon_{t+1} \right].

と表したとき、0\le t < T では \varepsilon_{t+1} と (y_t - \bar{y}_T) / \sum_{t=1}^T (y_t - \bar{y}_T)^2 が独立にならないため、

E\left[ \frac{ (y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \varepsilon_{t+1} \right] \ne E\left[　\frac{ (y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \right] E\left[ \varepsilon_{t+1} \right]

のように \varepsilon_{t+1} についての期待値を外に出すことができない。

以上から、E[\hat{\phi}] \ne \phi であり、\hat{\phi} は \phi の不偏推定量ではない。

まとめ

以上から、AR(1) モデルでは、観測値 y_t が誤差項 \varepsilon_t, \varepsilon_{t-1}, \varepsilon_{t-2}, ... と独立でないため、最小二乗法で求めた推定量が不偏推定量にならない。

補足： y_{t+k} と \varepsilon_t が独立でないことの説明

共分散を求めてみると、k\ge 0では0にならない：

\begin{aligned} Cov[y_{t+k}, \varepsilon_t] &= E[(y_{t+k} - \mu) \varepsilon_t] \\ &= E[(y_{t+k} - \mu) \{ (y_t - \mu) - \phi (y_{t-1} - \mu) \}]\\ &=Cov[y_{t+k}, y_t] - \phi Cov[y_{t+k}, y_{t-1}]\\ &= \frac{\sigma^2}{1 - \phi^2} \phi^{|k|} - \frac{\sigma^2}{1 - \phi^2} \phi \phi^{|k+1|}\\ &= \left\{ \begin{array}{ll} \sigma^2 \phi^k & k\ge 0 \\ 0 & k < 0 \end{array} \right. . \end{aligned}
畳屋民也畳屋民也

AR(1)過程でパラメータ推定量が一致推定量であることを示す

https://zenn.dev/link/comments/8a28e47602c1cf
上記の投稿で述べたように、AR(1)過程

y_{t+1} = c + \phi y_t + \varepsilon_{t+1}

（ただし、\varepsilon_t は分散が\sigma^2のホワイトノイズ）

のもと、T+1 個の観測値 \{ y_1, y_2, ..., y_T, y_{T+1}\}が与えられたときのパラメータ \{c, \phi, \sigma^2\} の推定量は、以下のように表せた：

\begin{aligned} \hat{c} & = \bar{y}_{T+1} - \hat{\phi} \bar{y}_T \\ \hat{\phi} &= \frac{\sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1})(y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \\ s^2 &= \frac{1}{T} \sum_{t=1}^T (y_{t+1} - \hat{c} - \hat{\phi} y_t)^2 \end{aligned}

ただし、

\bar{y}_{T+1} = \frac{1}{T} \sum_{t=1}^T y_{t+1}, \quad \bar{y}_{T} = \frac{1}{T} \sum_{t=1}^T y_{t}

とした。

今回は、これらがパラメータ \{c, \phi, \sigma^2\} の一致推定量となることを示す：

\plim_{T\to\infty}\hat{\phi}=\phi, \quad \plim_{T\to\infty}\hat{c}=c, \quad \plim_{T\to\infty}s^2=\sigma^2

回帰係数推定量 \hat{\phi} の一致性

https://zenn.dev/link/comments/fcd4d35340701b
上記投稿でも述べたとおり、\hat{\phi} を以下のように書き換えることができる：

\begin{aligned} \hat{\phi} &= \phi + \frac{\sum_{t=1}^T \varepsilon_{t+1} (y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2}\\ &= \phi + \frac{\frac{1}{T}\sum_{t=1}^T \varepsilon_{t+1} y_t - \frac{1}{T}\sum_{t=1}^T \varepsilon_{t+1} \bar{y}_T}{\frac{1}{T}\sum_{t=1}^T (y_t - \bar{y}_T)^2}. \end{aligned}

このとき、

\frac{1}{T}\sum_{t=1}^T \varepsilon_{t+1} y_t \xrightarrow{p} 0 \tag{1}
\frac{1}{T}\sum_{t=1}^T \varepsilon_{t+1} \bar{y}_T \xrightarrow{p} E[\varepsilon_{t+1}] \cdot \mu = 0 \tag{2}
\frac{1}{T}\sum_{t=1}^T (y_t - \bar{y}_T)^2 \xrightarrow{p} V[y_t] \tag{3}

が成り立つので、

\hat{\phi}\xrightarrow{p} \phi + \frac{0 - 0}{V[y_t]}=\phi

となる。

定数項推定量 \hat{c} の一致性

上述のように

\hat{c} = \bar{y}_{T+1} - \hat{\phi} \bar{y}_T

であるが、ここで

\bar{y}_{T+1} \xrightarrow{p} \mu, \quad \bar{y}_T \xrightarrow{p} \mu, \quad \hat{\phi} \xrightarrow{p} \phi

であるので、

\begin{aligned} \hat{c} \xrightarrow{p} \mu - \phi \mu = \mu (1 - \phi) = c \end{aligned}

が成り立つ。

なお、

\mu = \frac{c}{1-\phi}

を用いた。

誤差項分散推定量 s^2 の一致性

AR(1)モデルの定義より \varepsilon_{t+1} = y_{t+1} - \phi y_t - c と表すことができるので、

\begin{aligned} \sum_{t=1}^T \varepsilon_{t+1}^2 &= \sum_{t=1}^T (y_{t+1} - \phi y_t - c)^2\\ &= \sum_{t=1}^T \left\{ (y_{t+1} - \hat{\phi} y_t - \hat{c}) + (\hat{c} - c) + (\hat{\phi} - \phi) y_t \right\}^2\\ &= \sum_{t=1}^T \hat{\varepsilon}_{t+1}^2 + T(\hat{c}-c)^2 + (\hat{\phi}-\phi)^2 \sum_{t=1}^T y_t^2\\ &\quad + 2 (\hat{c}-c) \sum_{t=1}^T \hat{\varepsilon}_{t+1} + 2 (\hat{\phi}-\phi) \sum_{t=1}^T \hat{\varepsilon}_{t+1} y_t\\ &\quad + 2 (\hat{c}-c)(\hat{\phi}-\phi) \sum_{t=1}^T y_t \end{aligned}

が成立する。

なお、

\hat{\varepsilon}_{t+1} = y_{t+1} - \hat{\phi} y_t - \hat{c}

とした。

ここで、

\sum_{t=1}^T \hat{\varepsilon}_{t+1} = 0 \tag{4}
\sum_{t=1}^T \hat{\varepsilon}_{t+1}y_t = 0 \tag{5}

が成り立つことから、

\begin{aligned} \sum_{t=1}^T \varepsilon_{t+1}^2 &= \sum_{t=1}^T \hat{\varepsilon}_{t+1}^2 + T(\hat{c}-c)^2 + (\hat{\phi}-\phi)^2 \sum_{t=1}^T y_t^2 + 2 (\hat{c}-c)(\hat{\phi}-\phi) \sum_{t=1}^T y_t \end{aligned}

となる。
従って、移項して両辺を T でわることで、

\begin{aligned} s^2 &= \frac{1}{T} \sum_{t=1}^T \hat{\varepsilon}_{t+1}^2\\ &= \frac{1}{T} \sum_{t=1}^T \varepsilon_{t+1}^2 - (\hat{c}-c)^2 - (\hat{\phi}-\phi)^2 \frac{1}{T} \sum_{t=1}^T y_t^2\\ &\quad - 2 (\hat{c}-c)(\hat{\phi}-\phi) \frac{1}{T} \sum_{t=1}^T y_t \end{aligned}

が得られる。

ここで、

\begin{aligned} \hat{c}\xrightarrow{p} c, \quad \hat{\phi} \xrightarrow{p} \phi\\ \frac{1}{T} \sum_{t=1}^T \varepsilon_{t+1}^2 \xrightarrow{p} \sigma^2 \end{aligned}

に着目すると、

\begin{aligned} s^2\xrightarrow{p} \sigma^2 - 0 - 0\cdot E[y_t^2] - 2 \cdot 0 \cdot 0 \cdot E[y_t] = \sigma^2 \end{aligned}

となることが示される。

補足

(1) \frac{1}{T}\sum_{t=1}^T \varepsilon_{t+1} y_t \xrightarrow{p} 0

\{\varepsilon_{t+1} y_t\} が Martingale 差分系列であることを利用する。

詳細は下記書籍の §7.2 を参照のこと：
J. D. Hamilton, Time Series Analysis (1994, Princeton Univ. Press)
https://amzn.asia/d/eCnCBm3

(4) \sum_{t=1}^T \hat{\varepsilon}_{t+1} = 0, (5) \sum_{t=1}^T \hat{\varepsilon}_{t+1}y_t = 0

直接的な示し方

\begin{aligned} \sum_{t=1}^T \hat{\varepsilon}_{t+1} &= \sum_{t=1}^T (y_{t+1} - \hat{\phi} y_t - \hat{c})\\ &= \sum_{t=1}^T \left\{y_{t+1} - \hat{\phi} y_t - (\bar{y}_{T+1} - \hat{\phi} \bar{y}_T) \right\}\\ &= T \bar{y}_{T+1} - \hat{\phi} T \bar{y}_T - T (\bar{y}_{T+1} - \hat{\phi} \bar{y}_T)\\ &= 0 \end{aligned}
\begin{aligned} \sum_{t=1}^T \hat{\varepsilon}_{t+1}y_t &= \sum_{t=1}^T (y_{t+1} - \hat{\phi} y_t - \hat{c}) y_t\\ &= \sum_{t=1}^T \left\{ (y_{t+1} - \bar{y}_{T+1}) - \hat{\phi} (y_t - \bar{y}_T ) \right\} y_t\\ &= \sum_{t=1}^T \left\{ (y_{t+1} - \bar{y}_{T+1}) - \hat{\phi} (y_t - \bar{y}_T ) \right\} (y_t - \bar{y}_T) \quad ... (*)\\ &= \sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1})(y_t - \bar{y}_T) - \hat{\phi} \sum_{t=1}^T (y_t - \bar{y}_T )^2\\ &=0 \end{aligned}

(*) では、

\sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1}) =0, \quad \sum_{t=1}^T (y_t - \bar{y}_T ) =0

であることから

\sum_{t=1}^T \left\{ (y_{t+1} - \bar{y}_{T+1}) - \hat{\phi} (y_t - \bar{y}_T ) \right\} \bar{y}_T =0

が成り立つことを利用した。

一般的な示し方

下記の投稿のように計画行列 (design matrix) \boldsymbol{X} を用いて考える。
https://zenn.dev/link/comments/32302a7a779d30
https://zenn.dev/link/comments/db9e98ef7a0301

上記の投稿でも記述したように、回帰式

\boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}

が与えられたとき、最小二乗法による回帰係数推定量\hat{\boldsymbol{\beta}}は以下のように求められる：

\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}

従って、以下の関係が成り立つ：

\boldsymbol{X}^\top (\boldsymbol{Y} - \boldsymbol{X}\hat{\boldsymbol{\beta}}) = \boldsymbol{0}

ここで、\hat{\boldsymbol{\varepsilon}} = \boldsymbol{Y} - \boldsymbol{X}\hat{\boldsymbol{\beta}} と置くと、上記の式は以下のように表すことができる：

\begin{aligned} \boldsymbol{X}^\top \hat{\boldsymbol{\varepsilon}} &= \sum_{t=1}^T \boldsymbol{x}_t \hat{\varepsilon}_t\\ &= \boldsymbol{0}. \end{aligned}

今回の場合、

\boldsymbol{X} = \begin{pmatrix} \boldsymbol{x}^\top_1\\ \boldsymbol{x}^\top_2\\ \vdots \\ \boldsymbol{x}^\top_T \end{pmatrix} = \begin{pmatrix} 1 & y_1\\ 1 & y_2\\ \vdots & \vdots\\ 1 & y_T \end{pmatrix}, \quad \boldsymbol{Y} = \begin{pmatrix} y_2\\ y_3\\ \vdots \\ y_{T+1} \end{pmatrix}, \quad \hat{\boldsymbol{\varepsilon}} = \begin{pmatrix} \hat{\varepsilon}_2\\ \hat{\varepsilon}_3\\ \vdots \\ \hat{\varepsilon}_{T+1} \end{pmatrix}
\hat{\boldsymbol{\beta}} = \begin{pmatrix} \hat{c}\\ \hat{\phi} \end{pmatrix}

に対応するので、

\begin{aligned} \boldsymbol{X}^\top \hat{\boldsymbol{\varepsilon}} &= \sum_{t=1}^T \begin{pmatrix} 1\\ y_t \end{pmatrix} \hat{\varepsilon}_{t+1}\\ &= \begin{pmatrix} 0\\ 0 \end{pmatrix}. \end{aligned}

が成り立つことから、

\sum_{t=1}^T \hat{\varepsilon}_{t+1} = 0, \quad \sum_{t=1}^T \hat{\varepsilon}_{t+1}y_t = 0

が示される。

畳屋民也畳屋民也

AR(1)過程におけるパラメータ推定量の漸近分布

https://zenn.dev/link/comments/8a28e47602c1cf
上記の投稿で述べたように、AR(1)過程

y_{t+1} = c + \phi y_t + \varepsilon_{t+1}

では、T+1 個の観測値 \{ y_1, y_2, ..., y_T, y_{T+1}\}が与えられたときのパラメータ \{c, \phi, \sigma^2\} の推定量は、以下のように表せた（\varepsilon_t は分散が\sigma^2のホワイトノイズ）：

\begin{aligned} \hat{c} & = \bar{y}_{T+1} - \hat{\phi} \bar{y}_T \\ \hat{\phi} &= \frac{\sum_{t=1}^T (y_{t+1} - \bar{y}_{T+1})(y_t - \bar{y}_T)}{\sum_{t=1}^T (y_t - \bar{y}_T)^2} \\ s^2 &= \frac{1}{T} \sum_{t=1}^T (y_{t+1} - \hat{c} - \hat{\phi} y_t)^2 \end{aligned}

ただし、

\bar{y}_{T+1} = \frac{1}{T} \sum_{t=1}^T y_{t+1}, \quad \bar{y}_{T} = \frac{1}{T} \sum_{t=1}^T y_{t}

とした。

さらに、下記の投稿ではこれらの推定量がT\to\inftyで真の値 \{c, \phi, \sigma^2\} に確率収束する一致推定量であることを示した：
https://zenn.dev/link/comments/c43ebfbee4a69f

今回は、推定量 \{\hat{c}, \hat{\phi}, s^2\} の漸近分布について扱う。

結果

ホワイトノイズ \{\varepsilon_t \} の4次モーメント \mu_4 = E[\varepsilon_t^4] が有限の大きさをもつとする。

このとき、推定量 \{\hat{c}, \hat{\phi}, s^2\} は、漸近的に以下のような分布に従う（\xrightarrow{L}は、T\to\infty における分布収束を表す）：

定数項 c、係数 \phi

\begin{aligned} \sqrt{T} \left\{ \begin{pmatrix} \hat{c}\\ \hat{\phi} \end{pmatrix} - \begin{pmatrix} c\\ \phi \end{pmatrix} \right\} &\xrightarrow{L} N\left( \begin{pmatrix} 0\\ 0 \end{pmatrix}, \begin{pmatrix} \sigma^2 + \mu^2 (1 - \phi^2) & -\mu (1-\phi^2)\\ -\mu (1-\phi^2) & 1-\phi^2 \end{pmatrix} \right) \end{aligned}

特に、\hat{\phi}, \hat{c} 個別に見ると、

\begin{aligned} \sqrt{T}(\hat{c} - c) &\xrightarrow{L} N(0, \sigma^2 + \mu^2 (1 - \phi^2))\\ \sqrt{T}(\hat{\phi} - \phi) &\xrightarrow{L} N(0, 1-\phi^2) \end{aligned}

誤差項分散 \sigma^2

\begin{aligned} \sqrt{T} (s^2 - \sigma^2) &\xrightarrow{L} N(0, \mu_4 - \sigma^4) \end{aligned}

特に、誤差項 \{ \varepsilon_t \} が正規分布 \varepsilon_t \sim N(0, \sigma^2) に従うとき、

\begin{aligned} \sqrt{T} (s^2 - \sigma^2) &\xrightarrow{L} N(0, 2\sigma^4) \end{aligned}

導出

以下にならってベクトル・行列表記をとる：
https://zenn.dev/link/comments/32302a7a779d30

\begin{aligned} \boldsymbol{\beta} = \begin{pmatrix} c\\ \phi \end{pmatrix}, \quad \hat{\boldsymbol{\beta}} = \begin{pmatrix} \hat{c}\\ \hat{\phi} \end{pmatrix} \end{aligned}
\boldsymbol{x}_t = \begin{pmatrix} 1\\ y_t \end{pmatrix}, \quad \boldsymbol{X} = \begin{pmatrix} \boldsymbol{x}_1^\top\\ \boldsymbol{x}_2^\top\\ \vdots \\ \boldsymbol{x}_T^\top \end{pmatrix} = \begin{pmatrix} 1 & y_1\\ 1 & y_2\\ \vdots & \vdots \\ 1 & y_T \end{pmatrix}
\boldsymbol{Y} = \begin{pmatrix} y_2\\ y_3\\ \vdots \\ y_{T+1} \end{pmatrix}, \quad \boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_2\\ \varepsilon_3\\ \vdots \\ \varepsilon_{T+1} \end{pmatrix}

定数項 c、係数 \phi

以下の投稿と同様に、
https://zenn.dev/link/comments/db9e98ef7a0301

\begin{aligned} \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}\\ &= \boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon} \end{aligned}

と表すことができる。

従って、

\begin{aligned} \sqrt{T} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) &= \sqrt{T} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon}\\ &= \left( \frac{1}{T}\boldsymbol{X}^\top \boldsymbol{X} \right)^{-1} \sqrt{T} \left(\frac{1}{T}\boldsymbol{X}^\top \boldsymbol{\varepsilon} - \boldsymbol{0}\right) \end{aligned}

が成り立つ。

ここで、

E[\boldsymbol{x}_t \boldsymbol{x}_t^\top] = \boldsymbol{Q}

と置いて、

\sqrt{T} \left(\frac{1}{T}\boldsymbol{X}^\top \boldsymbol{\varepsilon} - \boldsymbol{0}\right) \xrightarrow{L} N(\boldsymbol{0}, \sigma^2 \boldsymbol{Q})
\frac{1}{T}\boldsymbol{X}^\top \boldsymbol{X} = \frac{1}{T} \sum_{t=1}^T \boldsymbol{x}_t \boldsymbol{x}_t^\top \xrightarrow{p} \boldsymbol{Q}

が成り立つ（詳細は省略）。

従って、

\begin{aligned} \left( \frac{1}{T}\boldsymbol{X}^\top \boldsymbol{X} \right)^{-1} \sqrt{T} \left(\frac{1}{T}\boldsymbol{X}^\top \boldsymbol{\varepsilon} - \boldsymbol{0}\right) \xrightarrow{L}& N(\boldsymbol{0}, \sigma^2 \boldsymbol{Q}^{-1}\boldsymbol{Q} (\boldsymbol{Q}^{-1})^\top)\\ &= N(\boldsymbol{0}, \sigma^2 \boldsymbol{Q}^{-1}) \end{aligned}

が成り立ち、さらに今回のケースでは

\boldsymbol{Q} = \begin{pmatrix} 1 & E[y_t]\\ E[y_t] & E[y_t^2] \end{pmatrix} = \begin{pmatrix} 1 & \mu \\ \mu & \frac{\sigma^2}{1 - \phi^2} + \mu^2 \end{pmatrix}

であることから

\begin{aligned} \sqrt{T} \left\{ \begin{pmatrix} \hat{c}\\ \hat{\phi} \end{pmatrix} - \begin{pmatrix} c\\ \phi \end{pmatrix} \right\} &\xrightarrow{L} N\left( \begin{pmatrix} 0\\ 0 \end{pmatrix}, \begin{pmatrix} \sigma^2 + \mu^2 (1 - \phi^2) & -\mu (1-\phi^2)\\ -\mu (1-\phi^2) & 1-\phi^2 \end{pmatrix} \right) \end{aligned}

が成り立つ。

誤差項分散 \sigma^2

下記の誤差項分散推定量の一致性の箇所で議論したように、
https://zenn.dev/link/comments/c43ebfbee4a69f

\begin{aligned} s^2 &= \frac{1}{T} \sum_{t=1}^T (y_{t+1} - \hat{\phi} y_t - \hat{c})^2\\ &= \frac{1}{T} \sum_{t=1}^T \varepsilon_{t+1}^2 - (\hat{c}-c)^2 - (\hat{\phi}-\phi)^2 \frac{1}{T} \sum_{t=1}^T y_t^2\\ &\quad - 2 (\hat{c}-c)(\hat{\phi}-\phi) \frac{1}{T} \sum_{t=1}^T y_t \end{aligned}

であるので、

\begin{aligned} \sqrt{T} (s^2 - \sigma^2) &= \sqrt{T} \left( \frac{1}{T} \sum_{t=1}^T \varepsilon_{t+1}^2 - \sigma^2 \right) \\ &\quad - \sqrt{T}(\hat{c}-c) \cdot (\hat{c}-c)\\ &\quad - \sqrt{T} (\hat{\phi}-\phi) \cdot (\hat{\phi}-\phi) \frac{1}{T} \sum_{t=1}^T y_t^2\\ &\quad - \sqrt{T} (\hat{c}-c) \cdot 2 (\hat{\phi}-\phi) \frac{1}{T} \sum_{t=1}^T y_t \end{aligned}

のように表すことができる。

このうち右辺第一項については、E[\varepsilon_{t+1}^4]=\mu_4<\inftyより

\begin{aligned} E[(\varepsilon_{t+1}^2 - \sigma^2)^2] &= E[\varepsilon_{t+1}^4] - \sigma^4\\ &= \mu_4 - \sigma^4 \end{aligned}

が成り立つことから

\sqrt{T} \left( \frac{1}{T} \sum_{t=1}^T \varepsilon_{t+1}^2 - \sigma^2 \right) \xrightarrow{L} N(0, \mu_4 - \sigma^4)

となる（詳細は省略）。

さらに右辺第二項以降については、

\hat{c} - c \xrightarrow{p} 0, \quad \hat{\phi} - \phi \xrightarrow{p} 0

かつ

\frac{1}{T} \sum_{t=1}^T y_t \xrightarrow{p} E[y_t] < \infty , \quad \frac{1}{T} \sum_{t=1}^T y_t^2 \xrightarrow{p} E[y_t^2] < \infty

が成り立つので、いずれもT\to \infty で消滅する。

以上から、

\sqrt{T} (s^2 - \sigma^2) \xrightarrow{L} N(0, \mu_4 - \sigma^4)

が示された。

なお、誤差項\{\varepsilon_{t+1}\}が正規分布 \varepsilon_{t+1} \sim N(0, \sigma^2)に従う場合、
E[\varepsilon_{t+1}^4] = 3\sigma^4 となるので、この場合は

\sqrt{T} (s^2 - \sigma^2) \xrightarrow{L} N(0, 2 \sigma^4)

となる。

補足

ベクトル・行列表現を用いた誤差分散推定量s^2の漸近分布の導出

上記ではs^2の漸近分布を求める際にスカラー値による表記を用いたが、こちらもベクトル・行列を用いた表現で導出してみる。

まず、推定量s^2は、以下のように表すことができる：

\begin{aligned} s^2 &= \frac{1}{T}\sum_{t=1}^T (y_{t+1} - \hat{\phi} y_t - \hat{c})^2\\ &= \frac{1}{T}\sum_{t=1}^T (y_{t+1} - \boldsymbol{x}_t^T \hat{\boldsymbol{\beta}})^2\\ &= \frac{1}{T} (\boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{\beta}})^\top (\boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{\beta}}) \end{aligned}

ここで、

\boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{\beta}} = (\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}) - \boldsymbol{X} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})

に着目してさらに式変形をすると、

\begin{aligned} s^2 &= \frac{1}{T} \left\{(\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}) - \boldsymbol{X} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \right\}^\top \left\{ (\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}) - \boldsymbol{X} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \right\}\\ &= \frac{1}{T} \left\{ (\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta})^\top (\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}) \right.\\ &\qquad - 2 (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \boldsymbol{X}^\top (\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta})\\ &\qquad \left. + (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \boldsymbol{X}^\top \boldsymbol{X}(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \right\}\\ &= \frac{1}{T} \boldsymbol{\varepsilon}^\top \boldsymbol{\varepsilon} - \frac{1}{T} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \boldsymbol{X}^\top \boldsymbol{X} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \end{aligned}

のようになる。

なお、二行目から三行目の変換の際に、

\boldsymbol{\varepsilon} = \boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}
\begin{aligned} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \boldsymbol{X}^\top (\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}) &= (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \boldsymbol{X}^\top \left\{ (\boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{\beta}}) + \boldsymbol{X} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \right\}\\ &= (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \boldsymbol{X}^\top \boldsymbol{X}(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})\\ & \qquad \because \boldsymbol{X}^\top (\boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{\beta}}) = \boldsymbol{0} \end{aligned}

という関係を用いた。

よって、

\sqrt{T}(s^2 - \sigma^2) = \sqrt{T} \left( \frac{1}{T} \boldsymbol{\varepsilon}^\top \boldsymbol{\varepsilon} - \sigma^2 \right) - \sqrt{T} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \left( \frac{1}{T} \boldsymbol{X}^\top \boldsymbol{X} \right) (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})

となる。
この式の右辺第一項は上述の通り

\sqrt{T} \left( \frac{1}{T} \boldsymbol{\varepsilon}^\top \boldsymbol{\varepsilon} - \sigma^2 \right) \xrightarrow{L} N(0, \mu_4 - \sigma^4)

であり、右辺第二項については

\begin{aligned} \sqrt{T} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) &\xrightarrow{L} N(\boldsymbol{0}, \sigma^2 \boldsymbol{Q}^{-1})\\ \frac{1}{T} \boldsymbol{X}^\top \boldsymbol{X} &\xrightarrow{p} \boldsymbol{Q}\\ (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) &\xrightarrow{p} \boldsymbol{0} \end{aligned}

より

\sqrt{T} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top \left( \frac{1}{T} \boldsymbol{X}^\top \boldsymbol{X} \right) (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \xrightarrow{p} \boldsymbol{0}

となる。

以上から、

\sqrt{T}(s^2 - \sigma^2) \xrightarrow{L} N(0, \mu_4 - \sigma^4)
畳屋民也畳屋民也

statsmodels で AR(1) モデル

以下のように AR(1) モデル

y_{t+1} = c + \phi y_t + \varepsilon_{t+1}, \quad E[\varepsilon_t]=0, \, E[\varepsilon_{t_1}\varepsilon_{t_2}]= \begin{cases} \sigma^2 \quad (t_1 = t_2)\\ 0 \quad (t_1 \ne t_2) \end{cases}\\
\begin{aligned} c = 1.2, \, \phi = 0.8\\ \sigma^2 = 1.0 \end{aligned}

に従う時系列データ\{y_t\}_{t=0}^Tを人工的に生成し、それをもとに statsmodels でモデルを作成した。

そのうえで、出力された値が理論をもとにどのように計算されるかを見ていく。

なお、今回はデータ長T=100、初期値y_0=-1.0 とした。

import numpy as np
import matplotlib.pyplot as plt

y0 = -1.0

c = 1.2
phi = 0.8
sigma = 1.0

nsample = 100

t = np.arange(0, nsample+1, 1)

fig, ax = plt.subplots(1,1)

list_ar1 = []
y_prev = y0
list_ar1.append(y0)
eps = np.random.normal(0, sigma, nsample)
for i in range(0,nsample):
  yt = c + phi * y_prev + eps[i]
  list_ar1.append(yt)
  y_prev = yt

ax.plot(t, list_ar1)

モデル作成

AR モデルの作成には、AutoReg クラスを使用する：

https://www.statsmodels.org/stable/generated/statsmodels.tsa.ar_model.AutoReg.html

AutoReg では、最小二乗法によりパラメータ \{c, \phi\} の推定を行う：

https://zenn.dev/link/comments/8a28e47602c1cf

from statsmodels.tsa.ar_model import AutoReg
from statsmodels.graphics.tsaplots import plot_predict

# モデル作成
## AR(p) モデルの次数 p は `lags` で指定
## 今回は定数項 c も含めるため、`trend="c"` のように指定した。
ar_mod = AutoReg(list_ar1, lags=1, trend="c")
ar_res = ar_mod.fit()

# 結果のプロット
fig, ax = plt.subplots()
ax.plot(list_ar1)  # 入力データのプロット
plot_predict(ar_res, ax=ax);  # 推測値と信頼区間（default では 95%）の表示

青線は観測値 \{y_t\}_{t=0}^T、オレンジ線はモデルによる予測値 \{\hat{y}_t\}_{t=1}^T である。

予測値は、\hat{y}_{t+1} = E[y_{t+1}\vert y_t] として定義され、今回の場合はパラメータ推定値 \hat{\phi}, \hat{c} を用いて以下のように求められる：

\hat{y}_{t+1} = \hat{c} + \hat{\phi} y_{t}

信頼区間の下限・上限は、誤差項分散推定値 s^2 を用いて以下のように表される：

\hat{y}_{t+1} \pm s z_{\alpha/2}

ここで、1-\alpha は信頼係数、z_{\alpha/2} は標準正規分布N(0, 1)の両側\alpha点である。今回はデフォルトの95%信頼区間を表示するので \alpha = 0.05 となる。

予測

観測値を取得した期間より未来の値を out-of-sample 予測することもできる。

以下では、20期先までの予測を行い、その結果を観測値とともにプロットしている。

fig, ax = plt.subplots()
plot_predict(ar_res, start=0, end=120, ax=ax)

ax.vlines(101, *ax.get_ylim(), color="red", linestyles="dashed")

ここでの青線は、区間1\le t \le T では観測値\{y_t\}を、t>T では予測値 \{\hat{y}_t\} を示している。

t = T + \tau \quad (\tau > 0) における予測値は、以下のように表される（下記投稿を参照）：

https://zenn.dev/link/comments/c1e650535540b8

\begin{aligned} \hat{y}_{T+\tau} = E[y_{T+\tau}\vert y_T ] &= \hat{\mu} + \hat{\phi}^k (y_T - \hat{\mu}) \\ &\xrightarrow{\tau\to\infty} \hat{\mu} \end{aligned}

ただし、\hat{\mu} = \hat{c} / (1 - \hat{\phi}) である。

一方で灰色で表示されているのは95%信頼区間である。t = T + \tau \quad (\tau > 0) においては、こちらは予測分散\hat{V}_tを用いて以下のように表される：

\hat{y}_{T + \tau} \pm \sqrt{\hat{V}_{T+\tau}} z_{\alpha/2}
\begin{aligned} \hat{V}_{T+\tau} = V[y_{T + \tau}\vert y_T] &= s^2 \frac{1-\hat{\phi}^{2\tau}}{1-\hat{\phi}^2}\\ &\xrightarrow{\tau\to\infty} \frac{s^2}{1-\hat{\phi}^2} \end{aligned}

モデル Summary

summary() メソッドを使うことで、モデルに関する統計量を表示できる。

print(ar_res.summary())

        AutoReg Model Results                             
==============================================================================
Dep. Variable:                      y   No. Observations:                  101
Model:                     AutoReg(1)   Log Likelihood                -137.147
Method:               Conditional MLE   S.D. of innovations              0.954
Date:                Tue, 01 Aug 2023   AIC                            280.294
Time:                        13:35:22   BIC                            288.110
Sample:                             1   HQIC                           283.457
                                  101                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.7493      0.348      5.025      0.000       1.067       2.432
y.L1           0.7151      0.057     12.467      0.000       0.603       0.828
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            1.3983           +0.0000j            1.3983            0.0000
-----------------------------------------------------------------------------

上記で3段に分かれたうちの2段目が、予測パラメータについての統計である。

ここでは、const は定数項 c を、y.L1 は係数 \phi を示す。

パラメータ推定値

coef 列はパラメータ推定値である。ここでは、

\begin{aligned} \hat{c} & = \bar{y}_{T} - \hat{\phi} \bar{y}_{T-1} = 1.7493\\ \hat{\phi} &= \frac{\sum_{t=0}^{T-1} (y_{t+1} - \bar{y}_{T})(y_t - \bar{y}_{T-1})}{\sum_{t=0}^{T-1} (y_t - \bar{y}_{T-1})^2} =0.7151 \end{aligned}

であったことがわかる。

ただし、

\bar{y}_{T-1} = \frac{1}{T}\sum_{t=0}^{T-1}y_t, \quad \bar{y}_{T} = \frac{1}{T}\sum_{t=1}^{T}y_t

なお、誤差項標準偏差 \sigma の推定値 \sqrt{s^2} は1段目の S.D. of innovations に書かれており、

\sqrt{s^2} = \sqrt{\frac{1}{T} \sum_{t=0}^{T-1} (y_{t+1} - \hat{c} - \hat{\phi} y_t)^2 } = 0.954

である。

パラメータ標準誤差

std err 列は標準誤差 \sqrt{\hat{V}[\hat{c}]} = 0.348, \sqrt{\hat{V}[\hat{\phi}]}=0.057 を示している。

それぞれ以下のように計算されている：

\begin{aligned} \hat{V}[\hat{\phi}] &= \frac{s^2}{\sum_{t=0}^{T-1} (y_t - \bar{y}_{T-1})^2}\\ \hat{V}[\hat{c}] &= \frac{s^2 \frac{1}{T}\sum_{t=0}^{T-1} y_t^2}{\sum_{t=0}^{T-1} (y_t - \bar{y}_{T-1})^2} \end{aligned}

なぜこうなるかを簡単に示す。下記の投稿より、
https://zenn.dev/link/comments/3dca8f9b5b6735

\begin{aligned} \sqrt{T} \left\{ \begin{pmatrix} \hat{c}\\ \hat{\phi} \end{pmatrix} - \begin{pmatrix} c\\ \phi \end{pmatrix} \right\} &\xrightarrow{L} N\left( \begin{pmatrix} 0\\ 0 \end{pmatrix}, \sigma^2 \boldsymbol{Q}^{-1} \right)\\ \boldsymbol{Q} &= E[\boldsymbol{x}_t \boldsymbol{x}_t^\top] = \begin{pmatrix} 1 & E[y_t]\\ E[y_t] & E[y_t^2] \end{pmatrix} \end{aligned}

が成り立つ。さらに、

\begin{aligned} \frac{1}{T}\boldsymbol{X}^\top \boldsymbol{X} &= \begin{pmatrix} 1 & \frac{1}{T} \sum_{t=0}^{T-1} y_t\\ \frac{1}{T} \sum_{t=0}^{T-1} y_t & \frac{1}{T} \sum_{t=0}^{T-1} y_t^2 \end{pmatrix}\\ &\xrightarrow{p} \boldsymbol{Q}\\ s^2 &\xrightarrow{p} \sigma^2 \end{aligned}

が成り立つので、T が十分に大きいと仮定して

P\left( \begin{pmatrix} \hat{c}\\ \hat{\phi} \end{pmatrix} \right) \approx N \left( \begin{pmatrix} c\\ \phi \end{pmatrix}, \frac{1}{T} s^2 \left( \frac{1}{T}\boldsymbol{X}^\top \boldsymbol{X} \right)^{-1} \right)

と見做すと、

\begin{aligned} \boldsymbol{X}^\top \boldsymbol{X} &= \begin{pmatrix} T & \sum_{t=0}^{T-1}y_t\\ \sum_{t=0}^{T-1}y_t & \sum_{t=0}^{T-1}y_t^2 \end{pmatrix}\\ \hat{V}[\hat{\phi}] &= s^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1}_{2,2}\\ &= \frac{s^2}{\sum_{t=0}^{T-1} (y_t - \bar{y}_{T-1})^2}\\ \hat{V}[\hat{c}] &= s^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1}_{1,1}\\ &= \frac{s^2 \frac{1}{T}\sum_{t=0}^{T-1} y_t^2}{\sum_{t=0}^{T-1} (y_t - \bar{y}_{T-1})^2} \end{aligned}

となる。

補足： ARIMAクラスとの違い

今回使用した AutoReg クラスのほかに、ARIMA クラスを用いて AR モデルを作成することもできる。

両者の違いは、モデルの推定方法である。

AutoReg は最小二乗法を用いるのに対して、ARIMA はカルマンフィルタを使用する。

https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima.model.ARIMA.html