基礎からはじめる時系列解析入門

tmtk

その他

2019.2.22

Topics

ARモデルとは
ARモデルのパラメータ推定
Pythonでモデリングしてみる
まとめ
参考文献

こんにちは。データサイエンスチームのtmtkです。
この記事では、時系列解析で使われるARモデルについて、基礎から解説します。

ARモデルとは

ARモデルというのは、時系列のデータ[latex] \{y_t\}_{t=1}^N \subset \mathbb{R}[/latex] に対して、

[latex]y_t = \sum_{i = 1}^m a_{i}y_{t-i} + v_t \, (t = m + 1, \ldots, N)[/latex]という生成過程を仮定したものです。ただし、[latex]m \in \mathbb{N}[/latex]はARモデルの次数とよばれる自然数、[latex] a_1, \ldots, a_m \in \mathbb{R}[/latex]は自己回帰係数とよばれる係数、[latex] v_t [/latex] は平均[latex]0[/latex]、分散[latex]\sigma^2[/latex]の正規分布[latex]N(0, \sigma^2)[/latex]に独立にしたがう確率変数です。

ARモデルのパラメータ推定

ARモデルのパラメータ[latex] m, (a_1, \ldots, a_m), \sigma^2[/latex]のうち、[latex](a_1, \ldots, a_m), \sigma^2[/latex]は最尤推定により推定することができます。残りのパラメータ[latex]m[/latex]の選択にはAICが使われますが、ここでは解説しません。
尤度関数は

[latex] L(a_1, \ldots, a_m, \sigma^2) = \prod_{t=m+1}^N p(v_{t} \mid a_1, \ldots, a_m, \sigma^2) = \prod_{t = m + 1}^N \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{\left(y_t – \sum_{i=1}^m a_{i}y_{t-i}\right)^2}{2\sigma^2}\right)[/latex]となります。これは多変数の線形回帰モデル

[latex]\left(\begin{array}{c}y_{m+1} \\ \vdots \\ y_N \end{array}\right) = \left(\begin{array}{ccc}y_m & \dots & y_1 \\ \vdots & \ddots & \vdots \\ y_{N-1} & \dots & y_{N-m}\end{array}\right)\left(\begin{array}{c} a_1 \\ \vdots \\ a_m\end{array}\right) + \left(\begin{array}{c}v_{m+1} \\ \vdots \\ v_N \end{array}\right)[/latex]の尤度関数と全く同じです。したがって、ARモデルのパラメータの最尤推定は、多変数の線形回帰の場合の最尤推定と同じように最小二乗法で

[latex] \hat{a}_i = \mathop{\rm argmin}_{a_i} \sum_{t=m+1}^N \left(y_t – \sum_{j=1}^m a_jy_{t-j}\right)^2 [/latex]と

[latex] \hat{\sigma}^2 = \frac{\sum_{t = m+1}^N \left(y_t – \sum_{i=1}^m \hat{a}_iy_{t-i}\right)^2}{N}[/latex]として推定することができます。

Pythonでモデリングしてみる

それでは、実際にPython 3でARモデルを使ってみます。
Pythonでは、StatsModelsというパッケージでARモデルが提供されていますが、ここではARモデルが多変数の線形回帰モデルと同一視できることを活かして、scikit-learnの線形回帰モデルを使ってみることにします。なお、上の説明で添え字が1-originだったのとは異なり、Pythonでは添え字が0-originであることに注意します。
まずは必要なライブラリをインポートします。

import numpy as np
from sklearn.linear_model import LinearRegression
np.random.seed(1234)

次に、コサインカーブに正規分布によるノイズを加えた、サンプルサイズ200の時系列データセットを作成します。

N = 200
std = 0.3
x = np.arange(0, N)
y = np.cos(x*np.pi/10) + np.random.normal(0, std, N)

今回は[latex]m=20[/latex]次のARモデルを構築することにします。線形回帰を実行するため、データを整形します。[latex]X = \left(\begin{array}{ccc}y_{m-1} & \dots & y_0 \\ \vdots & \ddots & \vdots \\ y_{N-2} & \dots & y_{N-m-1}\end{array}\right)[/latex]と[latex] Y = \left(\begin{array}{c}y_m \\ \vdots \\ y_{N-1} \end{array}\right) [/latex]を作成します。

m = 20
X = np.zeros((N - m, m))
for i in range(N - m):
    X[i] = y[i: i + m][::-1]
Y = y[m:]

最後に、パラメータの推定を行います。普通の線形回帰モデルと異なり、切片がモデルに含まれていないことに注意します。

lr = LinearRegression(fit_intercept=False)
lr.fit(X, Y)

これでARモデルのパラメータの学習が完了します。ARモデルの式

[latex] y_t = \sum_{i = 0}^{m-1} a_{i}y_{t-i-1} [/latex]を使うと、今後の値の予測ができます。ARモデルによる100期先までの予測値をy_predictedとして計算します。

prediction_len = 100
y_predicted = y
for i in range(prediction_len):
    t = [y_predicted[-m:][::-1]]
    y_predicted = np.append(y_predicted, [lr.predict(t)[0]])
y_predicted = y_predicted[N:]

すると、以下のように未来の値を予測ができていることがわかります。