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

tmtk

Data Science

2019.2.22

Topics

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

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

ARモデルとは

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

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

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

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

$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)$ となります。これは多変数の線形回帰モデル

$\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)$ の尤度関数と全く同じです。したがって、ARモデルのパラメータの最尤推定は、多変数の線形回帰の場合の最尤推定と同じように最小二乗法で

$\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$ と

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

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)

今回は $m=20$ 次のARモデルを構築することにします。線形回帰を実行するため、データを整形します。 $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)$ と $Y = \left(\begin{array}{c}y_m \\ \vdots \\ y_{N-1} \end{array}\right)$ を作成します。

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モデルの式

$y_t = \sum_{i = 0}^{m-1} a_{i}y_{t-i-1}$ を使うと、今後の値の予測ができます。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:]

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