世の中には色んな種類のデータがあり、売上の推移であったり、勉強へのモチベーションの移り変わりであったりといった、数字の並び順自体に意味があるデータがあります。この数字の並び順自体に意味があるデータは時系列データと呼ばれます。時系列データは、並び順（前後関係）に意味があり、機械学習でよく使うクロスセクションデータとは異りなます。データの時系列特有の要素を加味するため、時系列に特化したモデルがあります。
時系列に特化したモデルを扱うことで、時系列特有の情報を扱うことができます。今回は、時系列に特化したモデルの中でも、多変量に対する予測モデルである自己回帰モデル（VAR）を構築します。
このブログでは、とにもかくにも実践することに重きをおいています。体感することが実用で活かす一番の近道です。なので、まずPythonで動かしてみて体感してから、理論を見ています。今回のコードはこちら（Github）にあります。
Pythonを使ってVARを実践
多変量時系列データの予測はVARにより行うことができます。ただし、どんな時系列データでもよいわけではなく、定常な時系列データである必要があります。なので、定常的ではない（非定常的）な時系列データは、VARのモデルでは扱えません。定常性の具体的な説明は他に譲るとして、簡単に述べると、平均・自己共分散が時間 \( t \)に対して一定であり、トレンドをもたない（だんだん増加したり減少したりしない）データのことです。年々増えて行く体重の時系列データは非定常であり、勉強へのモチベーションが周期的に変動する時系列データは定常でしょうかね。
それでは、さっそく実践していきます。まず、多変量時系列データを用意し、データの中身を確認したあと、VARモデルを構築します。
データの準備
必要なモジュールのインポート
今回は『statsmodels』モジュールを使ってVARを構築します。statsmodelsはpipを使うことで、簡単にインストールできます。
from pathlib import Path import pandas as pd import numpy as np import statsmodels.api as sm from statsmodels.tsa import stattools from statsmodels.tsa.api import VAR import seaborn as sns import matplotlib.pyplot as plt result_dir_path = Path('result') if not result_dir_path.exists(): result_dir_path.mkdir(parents=True)
データの取得
多変量時系列データを取得します。今回扱うデータは、statsmodelsモジュールに組み込まれているデータであり、簡単に呼び出せます。このデータはアメリカの経済データであり、1959年Q1〜2009年Q3までの経済指標がクオーターごとに集計されたデータです。各変数の意味はこちら（statsmodels-macrodata.html）を御覧ください
macro_data = sm.datasets.macrodata.load_pandas().data macro_data.head()
データの加工
取得したデータのデータフレームのインデックスにクオーターを代入し、今回扱う変数である、「realgdp」「realcons」
「realinv」だけを取り出しています。
- realgdb：実質国内総生産
- realcons：実質個人消費支出
- realinv：実質国内総投資
macro_data['year'] = macro_data['year'].astype(int) macro_data.index = pd.date_range( str(macro_data['year'].min()), periods=len(macro_data), freq='Q' ) target_data = macro_data[['realgdp', 'realcons', 'realinv']]
データの確認
データの準備が終わったので、それぞれの時系列データを見てみます。
plot_data = target_data.stack() plot_data = plot_data.reset_index() plot_data.columns = ['quarter', 'variable', 'value'] g = sns.relplot( data=plot_data, x='quarter', y='value', col='variable', kind='line' ) g.set_axis_labels('quarter', '').set_titles(size=20).tight_layout() g.savefig(result_dir_path.joinpath('macro_data_plot.png')) g
上のグラフを見て分かる通り、３つの変数とも上がり幅は異なりますが、徐々に上昇しています。まさしくトレンドのあるデータであり、非定常的なデータです。ということは、このデータのままではVARモデルでは扱えません。そのことを定量的に測定するためにADF検定を行います。
単位根過程の確認
トレンドがあるような非定常的なデータでは、VARを扱えません。そこで、単位根過程であることを帰無仮説とするADF検定を行い、非定常的なデータかどうかを判断します。単位根過程という言葉が急に出できましたが、単位根過程とは、原系列 \(y_t\) が非定常過程であり、差分系列 \(\Delta y_t=y_t-y_{t+1}\) が定常過程である過程のことです。単位根過程であることの帰無仮説が棄却されないと、原系列は非定常過程である可能性が高いです。
それでは、今回扱うデータのそれぞれの変数において、ADF検定を実践していきます。
adf = stattools.adfuller(target_data['realgdp'], regression='ctt') print('t値 : {:.2f}, p値 : {:.1f}%'.format(adf[0], adf[1]*100)) print('データ数 : {}, 使用されたラグ数 : {}'.format(adf[3], adf[2])) print('検定統計量における棄却値 : 1%={:.2f}, 5%={:.2f}, 10%={:.2f}'.format( adf[4]['1%'], adf[4]['5%'], adf[4]['10%'] ))
adf = stattools.adfuller(target_data['realcons'], regression='ctt') print('t値 : {:.2f}, p値 : {:.1f}%'.format(adf[0], adf[1]*100)) print('データ数 : {}, 使用されたラグ数 : {}'.format(adf[3], adf[2])) print('検定統計量における棄却値 : 1%={:.2f}, 5%={:.2f}, 10%={:.2f}'.format( adf[4]['1%'], adf[4]['5%'], adf[4]['10%'] ))
adf = stattools.adfuller(target_data['realinv'], regression='ctt') print('t値 : {:.2f}, p値 : {:.1f}%'.format(adf[0], adf[1]*100)) print('データ数 : {}, 使用されたラグ数 : {}'.format(adf[3], adf[2])) print('検定統計量における棄却値 : 1%={:.2f}, 5%={:.2f}, 10%={:.2f}'.format( adf[4]['1%'], adf[4]['5%'], adf[4]['10%'] ))
３つの変数のすべてにおいて、p値が5%以上であり、単位根過程であるという帰無仮説を棄却できない事がわかります。つまり、非定常過程である可能性が高く、このままではVARモデルを構築できません。そこで、１つ前の値との差分を取ることで、データを変形させてみます。
１次の階差
これまでの結果で、取得したデータは非定常的であり、このままではVARモデルを構築できないことがわかりました。そこで、各変数１つ前の値との差分をとった形にデータを変形させてみます。
target_data_diff = target_data.diff().dropna() plot_data = target_data_diff.stack() plot_data = plot_data.reset_index() plot_data.columns = ['quarter', 'variable', 'value'] sns.relplot( data=plot_data, x='quarter', y='value', col='variable', kind='line' )
上のグラフを見ると、１次の階差を取ることで、平均が一定であり、トレンドが消失しているように見えます。同じくADF検定でも確認してみます。
adf = stattools.adfuller(target_data_diff['realgdp'], regression='ctt') print('t値 : {:.2f}, p値 : {:.1f}%'.format(adf[0], adf[1]*100)) print('データ数 : {}, 使用されたラグ数 : {}'.format(adf[3], adf[2])) print('検定統計量における棄却値 : 1%={:.2f}, 5%={:.2f}, 10%={:.2f}'.format( adf[4]['1%'], adf[4]['5%'], adf[4]['10%'] ))
adf = stattools.adfuller(target_data_diff['realcons'], regression='ctt') print('t値 : {:.2f}, p値 : {:.1f}%'.format(adf[0], adf[1]*100)) print('データ数 : {}, 使用されたラグ数 : {}'.format(adf[3], adf[2])) print('検定統計量における棄却値 : 1%={:.2f}, 5%={:.2f}, 10%={:.2f}'.format( adf[4]['1%'], adf[4]['5%'], adf[4]['10%'] ))
adf = stattools.adfuller(target_data_diff['realinv'], regression='ctt') print('t値 : {:.2f}, p値 : {:.1f}%'.format(adf[0], adf[1]*100)) print('データ数 : {}, 使用されたラグ数 : {}'.format(adf[3], adf[2])) print('検定統計量における棄却値 : 1%={:.2f}, 5%={:.2f}, 10%={:.2f}'.format( adf[4]['1%'], adf[4]['5%'], adf[4]['10%'] ))
上記の結果を見る通り、３つともp値が1%を切っており、単位根過程であるという帰無仮説は棄却されます。つまり、３つの変数とも単位根過程ではないです。
VARモデルの構築
ここまでの話で、今回扱うデータはそのままだと非定常過程である可能性が高く、１次の階差を取ることで定常過程にかわる可能性が高いことがわかりました。そこで、１次の階差をとった時系列データでVARモデルを構築します。
learning_data = target_data.diff().dropna() model = VAR(learning_data)
なんと、モデルの構築は、１行で終わりです。ありがたいですね。ただ、VARモデルは、過去さかのぼって何時点前まで（ラグと呼ぶ）をモデルに組み込むかで、精度が変わってきます。そこで、ラグを変えた結果を求めてみます。
model.select_order(15).summary()
上記では、15ラグまでのモデルの結果を表示しており、AICで見ると、ラグ3が最適であることがわかります。statsmodelsのVARは利便性が高く、指定した指標に基づいた最適なラグのモデルを自動で作成してくれる関数があります。
model_result = model.fit(maxlags=15, ic='aic') model_result.summary()
未来の予測
VARモデルを構築できたので、このモデルを使ってデータにおける将来の値を予測してみます。
steps = 20 result = model_result.forecast(model_result.endog, steps=steps) fig = plt.figure(i, figsize=(8, 8)) for i, variable_name in enumerate(['realgdp', 'realcons', 'realinv']): ax = fig.add_subplot(3, 1, i+1) y = model_result.endog[:, i] ax.plot( np.arange(len(y)), y ) ax.plot( np.arange(len(y), len(y)+steps), result[:, i] ) ax.set_title(variable_name, size=15) fig.tight_layout() plt.savefig(result_dir_path.joinpath('micro_data_predict.png'))
予測結果を見ると、早い段階で一定値に収束してしまっています。時系列データ予測モデルでよくあることですが、長期間の予測には向いてなさそうですね。短期間では、細かい変動が見られます。VARモデルの利点は、変数間の関係性を分析できることであり、次回そのことを書こうと思います。
VARの定義
１変数の時系列データに特化したモデルとしてARIMAなどがありますが、これを多変量に拡張したのがVARです。多変量に拡張することで、複数の変数による予測精度の向上や変数間の動学的関係の分析を行うことができます。変数間の分析として有名なのは、グレンジャー因果性、インパルス応答、分散分解です。これらに関しては、次回書こうと思います。
まず、VARの定義を見てみます。定義に使う各変数の定義は以下のとおりです。
- \( n \) : 変数の数
- \(p\) : ラグ
- \( \mathbf{y}_t \) : 時刻 \(t\) のときの \(n\) 個の時系列データのベクトル
- \( \mathbf{\phi_t} \) : 時刻 \( t \) の変数に対する \( n \times n \) の係数行列
- \( \mathbf{c} \) : \(n \times 1 \) 定数ベクトル
- \( \mathbf{\epsilon}_t \) : 時刻 \( t \) のときの \( n \) 個の誤差ベクトル
- \( \mathbf{\Sigma} \) : \( n \times n \) の共分散行列
数式で書いてみると、なんだかスッキリしていて簡単な式なように見えます。そこで、実際にはどのような式になるのかを \( n=2 \)、 \( p=1 \) という条件を与えて見てみます。
\(
\begin{cases}
\ y_{1t} = c_1 + \phi_{11} y_{1, t-1} + \phi_{12} y_{2, t-1} + \epsilon_{1t} \\\
y_{2t} = c_2 + \phi_{21} y_{1, t-1} + \phi_{22} y_{2, t-1} + \epsilon_{2t}
\end{cases} ,
\left(
\begin{array}{c}
\epsilon_{1t} \\
\epsilon_{2t} \\
\end{array}
\right)
\sim \mathbf{W.H.( \mathbf{\Sigma} )} \\
\mathbf{ \Sigma } =
\begin{pmatrix}
\sigma^2_{1} & \rho\sigma_1\sigma_2 \\
\rho\sigma_1\sigma_2 & \sigma^2_2
\end{pmatrix}
\)
ここで \( \rho = \mathrm{Corr} ( \epsilon_{1t}, \epsilon_{2t}) \) です。たかだか \( n=2 \)、 \( p=1 \) という簡易な条件でも、パラメータがそこそこあるのがわかります。実際、\( n \) 変量 \( \mathrm{VAR}(p)\) は、 \( n(np+1) + n(n+1)/2\) のパラメータを持つため、\( O(n^2p)\) です。つまり、計算時間はラグ数より変数の数の方が影響します。
終わりに
時系列データはLSTMを扱う際に以前勉強しましたが、実データを扱ってモデルを確認してみるということはしていなかったです。実際にやってみることで理解がはかどりました。そのために、こちらの本をとても参考にさせてもらいました。
