はじめに
時系列分析は、regression・classification・clusteringに同じくらいに使用頻度が高いものであるが、結構扱いが難しく敷居が高い。私も触りを学んだ程度で避けていたところがあった。しかしながら、今般、いろいろあって時系列分析に向き合う必要が出てきた。そこで知識のおさらいも兼ねて時系列解析を扱うライブラリーを使ってみたので、結果を報告したい。
実際に使ったモジュールは以下の通りである。
１．statsmodels
python・時系列解析で検索するとかなりの確率でヒットするほど有名なモジュールである。まさに王道といってもいいようにさえ思える。このstatsmodelsであるが、使ってみてRを意識した設計になっていると感じた。実装されている関数はおおよそRのものをカバーしている。また関数のパラメータの仕様・出力結果ともにRに類似している。
今回、教本として「時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装」を参考にした。掲載されるコードはRであるが、記載されるコードの解説とstatsmodelsのドキュメントを見ながらコード変換していくことで同じ結果を得ることができた。
２．fbprophet
fbprophetは、facebookが開発・公開する時系列分析用モジュールである。fbprophetの売りは、「柔軟性を保持しつつ、統計知識がそれほどないユーザでも使うことができる」ということである。実際に拍子抜けするほど、簡単にモデルを構築することができた。公開されてからそれほど時間を経ていないが、ネットで検索すると構築事例がヒットするようになっており、今後、大いに流行る予感がする。
３．tensorflow
tensorflowは時系列解析専門のモジュールでなく、ニューラルネットワーク開発用の汎用モジュールである。しかしながら、RNNやLSTMなどの手法に適合した関数が用意されることにより時系列分析にも応用できるようになった。また、tensorflowの解説本もRNNやLSTMをカバーしているものも多くなっており、敷居は確実に低くなっているように思う。
今回は、上記の１および２について報告する。tensorflowについては、別に報告する。
Dataの準備
今回使用するデータは東京電力の公開する電力使用実績データを使用した。（urlは以下）
http://www.tepco.co.jp/forecast/html/download-j.html
データは2016年からの3年分である。できれば5年分ほどあるといいのであるが、30分値の公開が2016年の電力自由化以降のこともあり、やむを得ないところである。
このデータは30分単位になっているもので、raw data状態で推移を示すと次のようになる。
これをそのままモデル化することも不可能ではないかもしれないが、問題を簡単にするため、1週間単位の平均値にリサンプリングしたものをデータとして使用する。
statsmodelsによるモデル化
statsmodelsによってモデルを構築するためにパラメータの設定に至るいくつかの手順をこなす必要がある。まず、リサンプリングしたデータと移動平均化したデータをプロットする。
これより周期的な変動を繰り返していることがうかがえる。そこで、この周期的変動がどこから来ているのかを自己相関・偏自己相関をプロットすることで可視化する。この自己相関・偏自己相関の可視化は、statsmodelの関数を用いることで実現できる。コードと結果を以下に示す。
import statsmodels.api as sm
fig = plt.figure(figsize=(20,10))
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
fig = sm.graphics.tsa.plot_acf(dat_w['実績'],ax=ax1)
fig = sm.graphics.tsa.plot_pacf(dat_w['実績'],ax=ax2)
これからも周期性の存在がわかる。また、この周期性は偏自己相関より1期前のデータとの相関に起因していることを読み取ることができる。なお、この自己相関・偏自己相関は数値の計算のみ実施し、別に可視化することも可能である。
続いて、1期前のデータの差分を可視化する。これよって定常性の有無が確認できる。
最後に長期的な季節性の有無について可視化する。リサンプリングしたデータから該当する月の平均値を減算し、月別の平均と標準偏差をプロットする。
左の月別平均値の差分からわかるようにほぼゼロという結果となった。右の標準偏差の大きさが気になるところであるが、おおよそ、予測対象時期に月が決定することで電力需要を推測できることが分かる。
いよいよモデル化に着手する。モデル化には①モデルの特定、②パラメータの特定、のステップを踏む必要がある。今回のモデルは、「自己相関が見られること」・「階差を取ることで相関がほぼなくなること（定常性）」・「長期的な季節周期が見られること」よりSARIMAを適用する。次にSARIMAに設定するパラメータの目安を決定する。
まず、階差を決定する。先のグラフよりおおよそ1階差分でよさそうであるが、statsmodelsの関数を使って検定する。検定にはkpss検定とADF検定がある。ここではkpss検定と用いる。コードはシンプルで以下のように書く。lagsオプションが検定したい階差を示す。
stats,p_value,lags,crit = sm.tsa.kpss(dat_w['実績'],lags=1)
print(P_value)
このp_value（p値）が十分に小さければ、検定対象であるlags=1を差分を取るべきということになる。今回の場合、p_value=0.0236ということなので1階差分をパラメータとして設定する。
次に自己相関と移動平均の次元数を決定する。statsmodelにはこれらを決定する関数が提供されているのでこれを使う。
res = sm.tsa.arma_order_select_ic(dat_w['実績'], ic='aic', trend='nc')
print(res['aic_min_order'])
オプションパラメータであるtrend='c'は、トレンドをどのようなものとして想定するかを示す。今回は見た目からトレンドなしとした。このaic_min_orderから(2,0)を次数として活用する。
最後に季節周期であるが、1年間を周期としてもいいのだが、自己相関の周期性から6か月に相当する26週間が1サイクルであることが見込まれることからこの値を用いることにする。つまり、冬と夏は電力需要が大きくなり、春と秋にかけて減少していくというモデルを表現することにした。
以上をモデルパラメータに設定して学習させれば完了である。以下に学習から予測までのコードとその結果を示す。
model = sm.tsa.SARIMAX(dat_w,order=(2,1,0),seasonal_order=(1,1,1,26))
result_sarima = model.fit()
fig = plt.figure(figsize=(20,5))
plt.plot(result_sarima.predict(start=1,end=150),label='予測')
plt.plot(dat_w['実績'],label='実績')
plt.legend()
plt.show()
order=(2,1,0)は、先のarma_order_select_icの結果とkpss検定の結果からパラメータを設定している。また、seasonal_order=(1,1,1,26)は、周期性から半年周期であるという仮定のもと設定している。今回はモデル化の方向性を確認するため、与えられた全データを学習データに用い、再現性を評価することとした。予測期間は、start=1とend=150より1期目から150期目までを予測範囲とした。0期目を除外したのは、1階差分をとっていることから0期目は予測されないものと解釈したためである。結果からすると、実績値をうまく捉えたモデルになっていることが見て取れる。実績データのない先のものについても、これまでの周期性をうまく捉えることができた。1点問題があるとすると、予測値が1期ずれているように感じることである。実際に予測値を1期ずらしたものを以下に示す。
これよりわかるように、1期ずらしたほうがフィットしているように見える。実際に二乗誤差平均を計算したところ、1期ずらしたほうがスコアが良かった。このズレは、過去の実装例において指摘はなかったことからすると、本モデル化にあたって考慮漏れ・パラメータ設定漏れがあったように思われる。この点については今後の調査事項としたい。
なお、学習データであるが、予測モデルのパラメータによって必要なデータ件数が制約を受ける。今回のケースでいうと、106件のデータがないと学習できない。総データ件数が107件であることを考えると学習データと試験データに分割する意味合いがほとんどない。sarimaxを適用する際の重要な考慮事項といえよう。
最後にstatsmodelsを用いたモデル化についてを総括する。
実際にモデル化を行ってみてわかったこととして難易度が高い。そもそも時系列分析自体が難易度が高いこともあるが、パラメータ設定に至る作業も時系列モデルについてそれなりの知識などがないと苦労しそうである。実は、今回のモデルであるSARIMAXに至るまでARやARIMAを試している。想定としては、複雑なモデルになるにつれて精度が向上するとみていたが、実際には訓練データに対する際限度はどのモデルもそれほどの差がなかった。一方、未知の期間の予測は、AR・ARIMAとものうまくいかなかった。つまり、モデルおよびパラメータの特定は、あらかじめ入念な分析の上で行うものであって試行錯誤により徐々に結論に近づけるという方法は難しいようである。ただ、うまく使いこなすことができれば、外因性をモデル化して組み込むなどといったことも実現できそうなので、使い道は大いにありそうと思われる。
fbprophetによるモデル化
続いてfbprophetを用いたモデル化について言及する。
まず、簡単にインストール方法について説明する。インストールはpipとanacondaと2パターンで可能である。私はcondaコマンド（conda install -c conda-forge fbprophet）でインストールを行った。なお、用意されているバージョンは、0.21で、linux,Mac,Win(32bit),Win(64bit)に対応している。
続いてモデル化である。fbprophetのモデル化の関数はひとつだけであり、statsmodelsのようにモデルの特定を必要としない。さらに設定すべきパラメータも重要なものは、yearly_seasonality,daily_seasonality,weekly_seasonalityの3つであり、True/Falseを設定するだけなので単純である。これによってfbprophet側が必要なパラメータを自動設定してモデルを構築する。仮にこれらのパラメータの設定も難しいようであれば、デフォルトのままにすることでこのTrue/Falseも自動設定してくれる。したがって、比較的単純なデータであればデフォルトのまま学習させればそれなりの精度のモデルが構築でき、その後に他のパラメータをチューニングしていくという方法で進めることができる。
今回の場合、データは週間単位でまとまっていること、statsmodelsの分析において季節周期の存在が確認できたことから、yearly_seasonality=Trueとしてそれ以外はFalseを設定する。モデル化のコードは以下のようになる。
from fbprophet import Prophet
model = Prophet(yearly_seasonality=True,
weekly_seasonality=False,
daily_seasonality=False)
続いて学習と予測のコードを示す。
tmp = dat_w.reset_index()
tmp.columns = ['ds','y']
result = model.fit(tmp)
predict = model.predict(pd.DataFrame(tmp['ds']))
predict.index = predict['ds']
fbprophetは、学習・予測時のデータに指定がある。予測データは時間と予測値（今回の場合は電力需要実績）を列名dsとyとしたDataFrameにする必要がある。これ以外の列名であるとエラーとなるので注意が必要である。予測については予測対象に日付を列名dsとして設定する必要がある。とりあえず、再現性を確認するため、予測期間は学習に用いた日付データをそのまま使用する。結果をプロットしたものが以下の通りである。
fbprophetには可視化において便利なメソッドがいくつか存在する。まず、plot()を使うと以下の通りに予測値・区間推定・実績値をセットで可視化することができる。
加えて、plot_components()を使うとtrend・yearly_seasonality・weekly_seasonalityなどに分解した推移を可視化することができる。
これによって「なぜ、このような予測結果になるのか」について説明することができる。fbprophetはstatsmodelsと異なり、パラメータ設定の多くがモジュールの裏側で行われているため、モデルの説明が難しい。しかしながら、このグラフがあることで、相応の説明ができる。
最後にデータをトレーニングデータとテストデータに分けて未知の期間の予測精度を評価する。トレーニングデータは前半の80件とし、テストデータは後半の27件とする。結果は以下の通りである。
赤線の左側がトレーニングデータの状況であり右側が未知のデータである。これから、未知のデータに対する予測精度も十分に確保できていると言えよう。
以上を総括すると、fbprephetはその売り文句に違わない柔軟性と使い勝手が確保できている上、予測精度も十分である。モジュールが裏側で行っていることでブラックボックス化している部分についても事後に出力されるグラフを使ってうまく説明できるであろう。
ここまでのまとめ
以上が、statsmodelsとfbprphetの実装に関する報告である。両者の比較を実装した感触からいうと、statsmodelsはプロ向き、fbprophetは一般ユーザー向きということになろう。今回は紹介していないが、statsmodelsは相当に充実した機能が提供されており、単純に予測モデルの構築ができるということに加えて、統計的にそのモデルが妥当であることをチェックする様々な検定方法が関数として提供されている。例えば、summary()というメソッドを使うとモデルに関する統計量がコンパクトにまとめられており、各係数のt検定値や予測結果の評価値であるAIC値などを確認することができる。
一方、fbprophetは、こうした統計的な評価は基本的にモジュールの裏側で行い、その結果として最適と見込まれるものをモデルとして自動的に設定する。パラメータの設定についてもドメイン知識から設定可能なものを中心に整理されており、時系列分析特有の概念を意識する必要がない。その代わりに「なぜその結果となるのか」については、fbprophetの公式ドキュメントに基づく理解が必要とされる。