今回は、代表的な時系列モデルである状態空間モデルをPythonで使う方法を書いていきます。 先日、『時系列データ分析とPython』という記事を書きましたが、今回はその内容の実装部分にあたります。(状態空間モデルって?という方はぜひ前回の記事を見て下さい!)
blog.tsurubee.tech
用いるデータは、時系列のデータセットとして大変有名な『各年ごとのナイル川流量』のデータです。
このデータは今回用いる統計ライブラリである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()
状態空間モデルの実質的なコードは最後のたった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())
推定したパラメータだけでなくLjung-Box検定などの検定の結果もまとめて出すことができます。
今回の状態推定はStatsModelsのUnobservedComponentsの中のfit関数を呼び出しただけです。公式リファレンスによると、 『Fits the model by maximum likelihood via Kalman filter.』と書いてあり、裏側ではカルマンフィルタによる最尤推定を行なっているようです。
参考

カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)
- 作者: 野村俊一
- 出版社/メーカー: 共立出版
- 発売日: 2016/09/08
- メディア: 単行本
- この商品を含むブログを見る
本のタイトルは『カルマンフィルタ』とピンポイントですが、状態空間モデルの基礎からカルマンフィルタ、粒子フィルタによる非線形非ガウス状態空間モデルまで幅広く学べる良書です。