モチベーション
時系列データは動画、FXの取引額の遷移、気温の遷移など幅広い分野で使用されるデータです。そのデータを学習し予測できればそれに合わせた施作を打つことができ、ちょっとした未来予測につながります。
そのような魅力的な時系列データを簡単に予測することができるライブラリPyFluxについて紹介します。
pyfluxの特徴
- scikit learnライクに記述可能
- 描画のメソッドが標準で実装済み
- 学習後のパラメータの確認が容易
時系列データ独特な特徴
- 連続
- 時間、空間、場所
- 潜在的な依存構造
- トレンド、シーズン、周期性など
- 動的な挙動
- 突然の変化、ゆるやかな変化
時系列データの予測に対する2つのアプローチ
- 予測
- 深層学習や線形回帰などのアプローチ
- 確率分布
- 確率分布を用いたアプローチ
PyFluxでは確率分布を用いたアプローチをしています。
ワークフロー
下記のワークフローで考えます。
- 'Build a time series model'でモデルを作成します。
- 'Inference'でモデルを学習します。
- 'Validation'で'Data'との整合性を取ります。
- 整合性が取れるまで1-3の処理を繰り返します。
- 整合性が十分になれば'Forecasting'で実際の予測を行います。
引用元
具体的な手法
今回の解析に使用するデータ
年間の太陽の黒点の数の遷移のデータです。
引用
import numpy as np
import pandas as pd
import pyflux as pf
from datetime import datetime
import matplotlib.pyplot as plt
%matplotlib inline
data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/sunspot.year.csv')
data.index = data['time'].values
plt.figure(figsize=(15, 5))
plt.plot(data.index, data['sunspot.year'])
plt.ylabel('Sunspots')
plt.title('Yearly Sunspot Data')
上記のコードで可視化しました。
データの特徴として
- 周期性がある
- 振幅は一定でない
ARIMA
ガウシアンモデルにより安定してスムーズな予測が可能
- 線形 + ガウシアンモデル
- AR + MA
- 定常データ
- 次数はAICで決定
モデルの定義
下記で行います。引数の意味は下記です。
- data: 学習に使用するデータ
- ar: ARIMAモデルの自己回帰の次数
- ma: 移動平均の次数
- family: 確率分布の指定
model = pf.ARIMA(data=data, ar=4, ma=4, target='sunspot.year', family=pf.Normal())
学習
学習モデルをデータに学習させる手法として'MLE'(最尤推定)を使用しています。
その他の予測手法について気になる方は下記をご覧下さい。
x = model.fit('MLE')
性能評価
モデルの情報をsummaryを使用すると確認できます。
- Log Likelihood: 確率を対数にしているため負の値になっています。このケースの場合はマイナスなので絶対値が小さいほど良いです。
- AIC： モデルのパラメータ数と対数尤度を考慮しており、低いほど良いとされています。
- BIC: AICの基準に加えて標本数も考慮しており、低いほど良いとされています。
- Std Error: エラーの標準偏差
- z: 標準偏差
- 95% C.I. : 各値が95%の範囲でどのレンジにあるかを示しています。
- P>|z|: P値とZスコアの関係
x.summary()
下記から分かることは定数の範囲が安定しておらず他の係数に比べて異常に大きいということです。
このような場合は他の学習手法も試した方が良いです。
Normal ARIMA(4,0,4)
======================================================= ==================================================
Dependent Variable: sunspot.year Method: MLE
Start Date: 1704 Log Likelihood: -1191.0131
End Date: 1988 AIC: 2402.0263
Number of observations: 285 BIC: 2438.5512
==========================================================================================================
Latent Variable Estimate Std Error z P>|z| 95% C.I.
======================================== ========== ========== ======== ======== =========================
Constant 11.5047 3.7053 3.1049 0.0019 (4.2422 | 18.7672)
AR(1) 1.4764 0.0738 20.0147 0.0 (1.3318 | 1.621)
AR(2) -0.2544 0.2076 -1.2257 0.2203 (-0.6613 | 0.1524)
AR(3) -0.9033 0.1831 -4.932 0.0 (-1.2622 | -0.5443)
AR(4) 0.452 0.0761 5.9379 0.0 (0.3028 | 0.6013)
MA(1) -0.2985 0.0727 -4.1047 0.0 (-0.441 | -0.156)
MA(2) -0.497 0.1008 -4.9307 0.0 (-0.6946 | -0.2995)
MA(3) 0.1824 0.0988 1.8468 0.0648 (-0.0112 | 0.376)
MA(4) 0.329 0.0837 3.9315 0.0001 (0.165 | 0.4931)
Normal Scale 15.8645
==========================================================================================================
下記は潜在変数の確率分布です。
この時系列は偏りのある分布で構成されていることが確認できます。
- figsize: 描画する図のサイズです。
model.plot_z(figsize=(15, 5))
学習データとのマッチ度合いを下記で確認できます。
model.plot_fit(figsize=(15, 10))
下記のように予測のサイズを調整して細かく見ることも可能です。
- h： どれだけ過去の時系列データを使用するかを意味しています。
他の引数については下記をご覧下さい
model.plot_predict_is(h=50, figsize=(15, 10))
実際に予測するとどれくらいの範囲のデータになるかが確認できます。色の濃さでどの程度の予測範囲か確認できます。時間が先になればなるほど予測が難しくなり、範囲が広くなっていることが確認できます。
- h: 予測するタイムステップ数を指定
- past_values
model.plot_predict(h=20, past_values=20, figsize=(15,5))
ARIMAX
時系列データにはない外部要因を入れて評価したい場合があります。そのような場合はARIMAXを使用します。
外部要素が入っているデータとして自動車における死亡事故のデータを使用します。
下記のデータからオイルショックとシートベルトの有無が死亡率に寄与していることを含めたい場合に使用してみます。下記のようなデータ形式です。
Unnamed: 0 time drivers seat_belt oil_crisis
time
1969.000000 1 1969.000000 1687 0.0 0.0
1969.083333 2 1969.083333 1508 0.0 0.0
1969.166667 3 1969.166667 1507 0.0 0.0
通常の処理と同様にデータを設定
data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/MASS/drivers.csv')
data.index = data['time']
plt.figure(figsize=(15, 5))
plt.plot(data.index, data['drivers'])
plt.ylabel('Driver Deaths')
plt.title('Deaths of Car Drivers in Great Britain 1969-84')
plt.plot()
data_d = data['drivers']
条件を加えて再度描画。データには別のパラメータが加えられているだけで変わらない。
data.loc[(data['time']>=1983.05), 'seat_belt'] = 1;
data.loc[(data['time']<1983.05), 'seat_belt'] = 0;
data.loc[(data['time']>=1974.00), 'oil_crisis'] = 1;
data.loc[(data['time']<1974.00), 'oil_crisis'] = 0;
plt.figure(figsize=(15, 5))
plt.plot(data.index, data['drivers'])
plt.ylabel('Driver Deaths')
plt.title('Deaths of Car Drivers in Great Britain 1969-84')
plt.plot()
外部要因の変数を変更して性能がどのような影響を受けるか見てみます。
formulaの部分がデータの指定と外部要因で何を有効にするかを指定しています。
- シートベルトのみ
model = pf.ARIMAX(data=data, formula='drivers~1+seat_belt',
ar=1, ma=1, family=pf.Normal())
- オイルショックのみ
model = pf.ARIMAX(data=data, formula='drivers~1+oil_crisis',
ar=1, ma=1, family=pf.Normal())
- オイルショックとシートベルトの両者を考慮
model = pf.ARIMAX(data=data, formula='drivers~1+seat_belt+oil_crisis',
ar=1, ma=1, family=pf.Normal())
評価指標であるAICとBIC、Log Likelihoodに着目してまとめるとオイルショックが外部要因として有用そうであることが分かります。
| Data | AIC | BIC | Log Likelihood |
|---|---|---|---|
| drivers~1+seat_belt+oil_crisis | 2606 | 2626 | -1297 |
| drivers~1+oil_crisis | 2571 | 2588 | -1282 |
| drivers~1+seat_belt | 2574 | 2590 | -1282 |
ここまではARIMAとARIMAXを見てきました。安定的かつ周期性のあるデータには有用ですがこれらのモデルには欠点が存在します。
欠点
- データから潜在的なプロセス（データの要因）は分解不可能
- 正規分布でなく非線形でもない
次に紹介するモデルはデータをノイズと予測データに分けて表す状態空間モデルを紹介します。
状態空間モデル
状態空間モデルではノイズの分布も模倣するため、実データとノイズの分解が可能なモデルになります。
観測データに対する数式
y_{t} = Z_t\alpha_t+\epsilon_t
潜在空間に対する数式
\alpha_{t} = T_t\alpha_{t-1}+\eta_t
観測データのノイズに関する数式
\epsilon_{t} \sim N(0, \Sigma^{2}_{\epsilon})
潜在空間のノイズに関する数式
\eta_{t} \sim N(0, \Sigma_{\eta})
- フィルター：過去のデータで傾向をフィルター
- 予測：過去のデータから予測
- スムージング：過去のデータでスムーズに
カルマンフィルター
パラメータαの更新方法がカルマンフィルターです。予測した値と実測値の差を考慮して時系列ごとに更新しているため、ノイズに強い特徴があります。
- $y_t - Z_t\alpha_t$は線形予測のエラー
- $K_t$はエラーに対してどれだけ反応するかをカルマンゲインで決定
\alpha_{t+1} = \alpha_{t} + K_t(y_t - Z_t\alpha_t)
残念ながらPyFluxでは状態空間モデルはカバーされていないのでもっとシンプルなローカルレベルモデルを使用します。
ローカルレベルモデル
観測データに対する数式
y_{t} = \mu_t+\epsilon_t
潜在空間に対する数式
\mu_{t} = \mu_{t-1}+\eta_t
観測データのノイズに関する数式
\epsilon_{t} \sim N(0, \sigma^{2}_{\epsilon})
潜在空間のノイズに関する数式
\eta_{t} \sim N(0, \sigma^{2}_{\eta})
状態空間モデルと過去の値を考慮する点は同一となります。
ナイル川のデータを見てみます。
nile = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/Nile.csv')
nile.index = pd.to_datetime(nile['time'].values,format='%Y')
plt.figure(figsize=(15,5))
plt.plot(nile.index,nile['Nile'])
plt.ylabel('Discharge Volume')
plt.title('Nile River Discharge');
plt.show()
下記で実装を行います。引数はARIMAと同一です。
model = pf.LocalLevel(data=nile, target='Nile', family=pf.Normal())
x = model.fit()
x.summary()
ローカルレベルモデルでは推定するのが観測データと異常データにおけるノイズの標準偏差のみなので下記のような情報になります。ノイズの標準偏差は実データとの差が取れないので推定値のみになります。
LLEV
======================================================= ==================================================
Dependent Variable: Nile Method: MLE
Start Date: 1871-01-01 00:00:00 Log Likelihood: -641.5238
End Date: 1970-01-01 00:00:00 AIC: 1287.0476
Number of observations: 100 BIC: 1292.258
==========================================================================================================
Latent Variable Estimate Std Error z P>|z| 95% C.I.
======================================== ========== ========== ======== ======== =========================
Sigma^2 irregular 15098.5410
Sigma^2 level 1469.13678
==========================================================================================================
model.plot_fit(figsize=(15,10))
下記がモデルが学習した結果を描画したものです。
予測値がスムージングされているので全体の傾向を見るのに適しています。
図は上から
- 予測値と予測値の95%で取りうる値と実データ
- 予測値と予測値の95%で取りうる値
- ノイズ
になります。
Gaussian local level model notebook
動的自己回帰モデル（ランダムウォーク）
先ほどのローカルレベルモデルではスムーズになりすぎる問題があります。ノイズではなく必要なスパイクを残したい場合は動的回帰モデルを使用する方がベターです。
y_{t} = \sum^{P}_{i=1}\phi_{i,t}y_{t-i}+\epsilon_t
潜在空間に対する数式
\phi_{i,t} = \phi_{i,t-1}+\eta_t
観測データのノイズに関する数式
\epsilon_{t} \sim N(0, \sigma^{2}_{\epsilon})
潜在空間のノイズに関する数式
\eta_{t} \sim N(0, \sigma^{2}_{\eta})
具体的な例として下記のようなケースになります。
上記のようなスムーズにしすぎる問題を解決
自己回帰モデルと同様のデータを使用して性能を比較してみます。
data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/sunspot.year.csv')
model= pf.DAR(data=data, ar=9, integ=0, target='sunspot.year')
x = model.fit('MLE')
x.summary()
下記のように動的自己回帰モデルの方が若干、性能が良いことが確認できます。
| Method | AIC | BIC | Log Likelihood |
|---|---|---|---|
| DAR | 2380 | 2420 | -1179 |
| AR | 2402 | 2438 | -1191 |
また各変数ごとの波形の再現率が確認できます。
Dynamic Regression Model notebook
動的自己回帰モデルの欠点
- ガウシアンモデルのみ
- ガウシアンのフィルター
- 線形の要素を外していない
スコア駆動時系列モデル
ガウシアンモデル以外のモデルを使用して問題に適した分布を使用したい場合にこのモデルを使用します。
ここで着目するのが更新部分です。
ARIMAの予測式
y_{t+1} = \phi y_t + \theta(y_t - \mu_t)
上記の更新式の場合は予測値から平均値を引いた部分が更新式となっています。これは予測値が分布の平均から大きくずれることを防ぐことを意味しています。
スコア駆動時系列モデルの導出
上記の式の更新部分の導出ができればその部分を置き換えることで分布の置き換えが可能になります。ではここでは正規分布の場合を用いて
(y_t - \mu_t)
を導出します。
予測データが正規分布に基づいて生成されると仮定
p(y_t) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(y_t-\mu)^2}{2\sigma^2})
上記を対数に置き換え
\log p(y_t) = \frac{1}{2}\log(2\pi\sigma^2)-\frac{1}{2}\frac{(y_t-\mu)^2}{\sigma^2}
平均値で微分
\nabla_{\mu_t}\log p(y_t) = \frac{(y_t-\mu)^2}{\sigma^2}
平均値で2回微分
\nabla^{2}_{\mu_t}\log p(y_t) = \frac{1}{\sigma^2}
比を導出（平均の変化度合い
\frac{\nabla_{\mu_t}\log p(y_t)}{-\nabla^{2}_{\mu_t}\log p(y_t)} = (y_t-\mu)
ここまでのプロセスを整理すると
- 予測関数の分布を定義
- 予測関数の分布を対数化
- 2.を対数を微分
- 2.を２回微分
- 3.を4.で割りマイナスをかける
上記の方法を使用すれば更新部分の分布の置き換えが可能になります。
ポアソン分布によるスコア駆動時系列モデル
最後の表現を綺麗にするためλを置き換え
\lambda_t = \exp(\theta_t)
ポアソン分布
p(y_t) = \frac{\exp(\theta_t)^{y_t}\exp^{-\exp(\theta_t)}}{y_t!}
対数をとる
\log p(y_t) = y_t\theta_t-\exp(\theta_t)-\log(y_t!)
θで微分
\nabla_{\theta_t}\log p(y_t) = y_t-\exp(\theta_t)
θで2回微分
\nabla^2_{\theta_t}\log p(y_t) = -\exp(\theta_t)
比を導出
\frac{\nabla_{\mu_t}\log p(y_t)}{-\nabla^{2}_{\mu_t}\log p(y_t)} = \frac{1}{\exp(\theta_t)}(y_t-\exp(\theta_t))
expをスケーリング係数として予測値が平均値（ポアソン分布の場合はλが平均値のため）より離れすぎないように更新される
実データへの適用
では実際にポアソン分布を用いて、時系列データを予測します。ポアソン分布は離散データを扱える分布なので今まで使用した連続データとは異なる離散データに対して適用します。
使用するのはゴールスコアの時系列データで今回のモデルに適した離散データです。
leicester = pd.read_csv('http://www.pyflux.com/notebooks/leicester_goals_scored.csv')
leicester.columns= ["Time","Goals","Season2"]
plt.figure(figsize=(15,5))
plt.plot(leicester["Goals"])
plt.ylabel('Goals Scored')
plt.title('Leicester Goals Since Joining EPL');
plt.show()
ガウシアン以外のモデルをローカルレベルモデルに適用します。
model = pf.LocalLevel(data=leicester, target='Goals', family=pf.Poisson())
モデルを学習します。
x = model.fit(iterations=5000)
x.summary()
出力結果は下記のようになります。
モデルの学習には"Black-Box Variational Inference"を使用しており、先ほど導出したスコア関数を勾配としてデータの分布とモデルの分布のKL距離を最小化するように学習します。ELBOは'Evidence Lower Bound'の略称で大きければ大きいほど分布間の距離が小さくなるのでこの値がおおきくなるように最適化されます。
ELBOの詳細が気になる方は下記をご覧ください。
BBVIを使用している理由は最尤推定に使用するMCMCだと時間的なコストが大きいためです。
10% done : ELBO is -251.57569112887904, p(y,z) is -166.3745900900033, q(z) is 85.20110103887575
20% done : ELBO is -251.0404488657637, p(y,z) is -163.08148512110355, q(z) is 87.95896374466014
30% done : ELBO is -248.3826857773143, p(y,z) is -160.95781417531765, q(z) is 87.42487160199666
40% done : ELBO is -243.78188227888018, p(y,z) is -158.15537328287218, q(z) is 85.62650899600798
50% done : ELBO is -240.84398713641178, p(y,z) is -156.96138472006595, q(z) is 83.88260241634582
60% done : ELBO is -238.7879147438848, p(y,z) is -156.14707335243622, q(z) is 82.6408413914486
70% done : ELBO is -237.27773593445806, p(y,z) is -155.37226587768643, q(z) is 81.90547005677165
80% done : ELBO is -234.66469485120308, p(y,z) is -153.9445088318954, q(z) is 80.72018601930766
90% done : ELBO is -234.21050603276024, p(y,z) is -153.86353460095933, q(z) is 80.34697143180091
100% done : ELBO is -232.41596530920148, p(y,z) is -153.02493554086158, q(z) is 79.3910297683399
Final model ELBO is -232.6451770092436
Poisson Local Level Model
======================================================= ==================================================
Dependent Variable: Goals Method: BBVI
Start Date: 0 Unnormalized Log Posterior: -153.3266
End Date: 74 AIC: 308.65314788761225
Number of observations: 75 BIC: 310.97063600114853
==========================================================================================================
Latent Variable Median Mean 95% Credibility Interval
======================================== ================== ================== =========================
Sigma^2 level 0.0404 0.0404 (0.0353 | 0.0462)
==========================================================================================================
学習データに対する予測分布を確認してみます。
model.plot_fit(figsize=(15,10))
Non-Gaussian Local Level Model notebook
今回使用したnotebook
参考
Ross Taylor | Time Series for Python with PyFlux