Pythonで状態空間モデルを使う（StatsModels）

今回は、代表的な時系列モデルである状態空間モデルをPythonで使う方法を書いていきます。先日、『時系列データ分析とPython』という記事を書きましたが、今回はその内容の実装部分にあたります。（状態空間モデルって？という方はぜひ前回の記事を見て下さい！）

blog.tsurubee.tech

用いるデータは、時系列のデータセットとして大変有名な『各年ごとのナイル川流量』のデータです。

f:id:hirotsuru314:20170702224936p:plain

このデータは今回用いる統計ライブラリであるStatsModelsの中のデータセットにも入っています。

statsmodels/nile.csv at master · statsmodels/statsmodels · GitHub

さっそくですが、このデータに状態空間モデルをあてはめてみます。

import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt
%matplotlib inline
import statsmodels.api as sm

#日付形式で読み込む
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y') 
#CSV読み込み
data = pd.read_csv('nile.csv', index_col='year', date_parser=dateparse)

# ローカルレベルモデルによる状態推定
model = sm.tsa.UnobservedComponents(data, 'local level')
result = model.fit()
fig = result.plot_components()

f:id:hirotsuru314:20170813224510p:plain

　状態空間モデルの実質的なコードは最後のたった3行だけです。
　黒線のObserved（観測データ）に対して、青線が推定した状態（潜在変数）になります。
　今回あてはめたモデルは状態空間モデルの中でも一番単純なローカルレベルモデルと呼ばれるモデルです。このモデルはランダムウォークプラスノイズモデルとも呼ばれ、ランダムウォークで推移している潜在変数（真の状態）に、何かしらのノイズ（不規則要素）がのっかって、観測データが得られるといったシンプルなモデルです。
　このようにStatsModelsを使うと、とても短いコードで状態空間モデルを使うことができます。下記のリンクがStatsModelsの状態空間モデル部分のリファレンスです。

statsmodels.tsa.statespace.structural.UnobservedComponents — statsmodels 0.9.0 documentation

今回の例では

model = sm.tsa.UnobservedComponents(data, 'local level')

のように第2引数として'local level'を渡していますが、上記リンクのStatsModelsリファレンスに記載のとおり、ここが実際にどのようなモデルを用いるかを決めます。したがって、例えば引数に、‘local linear trend’を渡せばローカルレベルモデルに線形のトレンドを加えたモデルに拡張することができます。また、『seasonal=12』といった引数を渡せば月の季節性成分を入れることも容易にできます。
　実際のモデリングのプロセスでは、データをよく再現するモデルを見つけるため、様々な要素（トレンド成分、季節性成分など）を入れながら、モデルの試行錯誤を行います。その上で参考となる指標も、下記のように簡単にsummaryとして出すことができます。

print(result.summary())

f:id:hirotsuru314:20170813230712p:plain

　推定したパラメータだけでなくLjung-Box検定などの検定の結果もまとめて出すことができます。
　今回の状態推定はStatsModelsのUnobservedComponentsの中のfit関数を呼び出しただけです。公式リファレンスによると、『Fits the model by maximum likelihood via Kalman filter.』と書いてあり、裏側ではカルマンフィルタによる最尤推定を行なっているようです。