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

こんにちは。データサイエンスチームの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:]

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

まとめ

この記事では、時系列解析の入門として、ARモデルを紹介しました。ARモデルのパラメータ推定の方法や線形回帰との関係について説明し、Python 3で実際にARモデルの学習と予測を行いました。
次回の記事では、AICによるパラメータmの選択について解説します。

参考文献

  • 樋口知之『予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで』
  • 北川源四郎『時系列解析入門』
  • 線形回帰 – Wikipedia
AWS EC2リザーブドインスタンス 10%OFFキャンペーン

あなたにおすすめの記事