Pythonによる状態空間モデル¶
In [1]:
# 基本のライブラリを読み込む
import numpy as np
import pandas as pd
from scipy import stats
# グラフ描画
from matplotlib import pylab as plt
import seaborn as sns
%matplotlib inline
# グラフを横長にする
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
# 統計モデル
import statsmodels.api as sm
データの読み込み¶
In [2]:
# 日付形式で読み込む
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
data = pd.read_csv('AirPassengers.csv', index_col='Month', date_parser=dateparse, dtype='float')
data.head()
Out[2]:
In [3]:
# 日付形式にする
ts = data['#Passengers']
ts.head()
Out[3]:
In [4]:
# プロット
plt.plot(ts)
Out[4]:
ローカルレベルモデルの推定¶
In [5]:
# ローカルレベルモデルの推定
mod_local_level = sm.tsa.UnobservedComponents(ts, 'local level')
# 最尤法によるパラメタの推定
res_local_level = mod_local_level.fit()
# 推定されたパラメタ一覧
print(res_local_level.summary())
# 推定された状態・トレンドの描画
rcParams['figure.figsize'] = 15, 15
fig = res_local_level.plot_components()
ローカル線形トレンドモデルの推定¶
In [6]:
# ローカル線形トレンドモデル
mod_trend = sm.tsa.UnobservedComponents(
ts,
'local linear trend'
)
# 最尤法によるパラメタの推定
# ワーニングが出たのでBFGS法で最適化する
res_trend = mod_trend.fit(method='bfgs')
# 推定されたパラメタ一覧
print(res_trend.summary())
# 推定された状態・トレンドの描画
rcParams['figure.figsize'] = 15, 20
fig = res_trend.plot_components()
季節変動の取り込み¶
In [7]:
# 季節変動ありのローカルレベルモデル
mod_season_local_level = sm.tsa.UnobservedComponents(
ts,
'local level',
seasonal=12)
# まずはNelder-Meadでパラメタを推定し、その結果を初期値としてまた最適化する。2回目はBFGSを使用。
res_season_local_level = mod_season_local_level.fit(
method='bfgs',
maxiter=500,
start_params=mod_season_local_level.fit(method='nm', maxiter=500).params,
)
# 推定されたパラメタ一覧
print(res_season_local_level.summary())
# 推定された状態・トレンド・季節の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = res_season_local_level.plot_components()
In [8]:
# 季節変動ありのローカル線形トレンドモデル
mod_season_trend = sm.tsa.UnobservedComponents(
ts,
'local linear trend',
seasonal=12)
# まずはNelder-Meadでパラメタを推定し、その結果を初期値としてまた最適化する。2回目はBFGSを使用。
res_season_trend = mod_season_trend.fit(
method='bfgs',
maxiter=500,
start_params=mod_season_trend.fit(method='nm', maxiter=500).params,
)
# 推定されたパラメタ一覧
print(res_season_trend.summary())
# 推定された状態・トレンド・季節の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = res_season_trend.plot_components()