記事について
以前SARIMAXで予測したTeamsのレスポンスデータを利用し、状態空間モデルの構築練習をしてみました。。
状態空間モデルの理論は学習中のため、ここでは理論を一旦置いておき、分析パラメーターを変更しながら、それっぽい予測ができることを一先ずのゴールとしています。
以前のSARIMAXの記事
https://qiita.com/ramutarafarm/items/bc6ed907a0a44177672e
とっても参考になった記事
https://logics-of-blue.com/python-state-space-models/
環境、利用データについてはSARIMAX記事を参照ください
ローカルレベルモデル
statsmodelsのサイトには13のモデルがありますが、まずはよく紹介されているローカルレベルモデルで分析。
実質2行、すごく簡単！ですが、、この裏でカルマンフィルタやら最尤推定やら、私がまだ理解できてないことが行われているんですね。。
mod_local_level = sm.tsa.UnobservedComponents(df_30T[["teams"]],
'local level')
res_local_level = mod_local_level.fit()
# 各コンポーネント（local level）の描画
fig = plt.figure(figsize=[9,12])
res_local_level.plot_components(fig=fig)
測定値と予測値のグラフ
ローカルレベルで予測した部分のグラフ
（今回はローカルレベルだけで分析しているので、ローカルレベルで全部分析しました、、と言うような結果になっている）
はじめ上記の結果を見て、「あれ、ローカルレベルだけで完全に分析できるのか。状態空間モデルすごいな」とおもってましたが、、
以下で実施した9/19～22の予測線を見てみてそういう事か、、と理解しました。
pred_l = res_local_level.get_prediction('2021-09-19', '2021-09-22')
pred_mean = pred_l.predicted_mean
pred_ci = pred_l.conf_int(alpha = .05)
plt.figure(figsize=(12, 6))
plt.plot(df_30T[df_30T.index > '2021-09-18']['teams'], label="train_data")
plt.plot(pred_mean, c="r", label="predict")
# 信頼区間を描画しても発散するだけなのでコメントアウト
#plt.fill_between(pred_ci.index,
# pred_ci.iloc[:, 0],
# pred_ci.iloc[:, 1], color='g', alpha=.2)
plt.legend()
plt.show()
予測結果
きっちり1コマ遅れの相関になってます。
また予測も80辺りで一本調子です。（実際にはノイズが乗っているので微妙に変動してます）
ただこれはローカルモデルの数式を見れば当たり前で、ローカルモデルでは一つ前の状態$\mu_{t-1}$にノイズを足して$\mu_t$を算出し、当該$\mu_t$にさらにノイズを足しているだけなので、このような結果になります。
（statsmodelsのサイト説明では、$\mu$はtrendとなってますが、このデータではイメージ定数項って感じですね）
Y_t = \mu_t + \varepsilon_t\\
\mu_t = \mu_{t-1} + \eta_t
ローカルレベル＋季節成分モデル
それでは次に季節成分を考慮して分析してみます。
季節成分はSARIMAX時と同様に48となります。（30分間隔データなので24時間周期になります）
mod_seasonal_local_level = sm.tsa.UnobservedComponents(df_30T["teams"],
'local level',
seasonal=48,
)
res_seasonal_local_level = mod_seasonal_local_level.fit()
# 各コンポーネント（local level,seasonal）の描画
fig = plt.figure(figsize=[9,12])
res_seasonal_local_level.plot_components(fig=fig)
測定値と予測値のグラフ
前述のローカルレベルだけで分析した方がよい結果に見えますね。。
ローカルレベルで予測した部分のグラフ
季節成分で予測した部分のグラフ
今回はローカルレベルと季節成分で分析しているので、コンポーネントグラフもそれぞれ描画されてます。
ちなみに、statusmodelsの状態空間モデルでは、各コンポーネントの総和が予測結果となるので、各コンポーネントの分析結果を見ながら要因を検討しやすいですね。
数式はローカルレベルに$\gamma_t$（季節成分）が足されているだけですが、、季節成分の計算方法は分からなかったです。48周期で階差取ってるんでしょうか？？
Y_t = \mu_t + \varepsilon_t + \gamma_t\\
\mu_t = \mu_{t-1} + \eta_t\\
\gamma_t = ???
測定値と予測値のグラフからあまり良い分析になってないのは明白ですが、いったんこのモデルで9/19～22の予測線を見てみてみます。
pred_s = res_seasonal_local_level.get_prediction('2021-09-19', '2021-09-22')
pred_mean = pred_s.predicted_mean
pred_ci = pred_s.conf_int(alpha = .05)
plt.figure(figsize=(12, 6))
plt.plot(df_30T[df_30T.index > '2021-09-18']['teams'], label="train_data")
plt.plot(pred_mean, c="r", label="predict")
#plt.fill_between(pred_ci.index,
# pred_ci.iloc[:, 0],
# pred_ci.iloc[:, 1], color='g', alpha=.2)
plt.legend()
plt.show()
季節成分を入れたおかげで予測線はそれっぽい変動をしてますが、、
本データのドメイン知識から「60を下回るのはおかしい」と分かるので、、単純に季節成分のせるだけではダメなようです。。
ローカルレベル＋季節成分＋自己回帰モデル
単純に季節節分を載せるだけではダメと言うのは上で分かりましたが、、
ではどうすればよいか、、は分からなかったので、パラメーターを総当たりしていろいろやって見た結果、以下のようになりました。
（コメントアウトはいろいろ検討した残骸です）
mod_seasonal_local_level = sm.tsa.UnobservedComponents(df_30T["teams"],
'local level',
#'local linear deterministic trend',
seasonal=48,
autoregressive=2,
#cycle=True,
#method='newton',
)
res_seasonal_local_level = mod_seasonal_local_level.fit()
# 各コンポーネント（local level,seasonal,autoregressive）の描画
fig = plt.figure(figsize=[9,18])
res_seasonal_local_level.plot_components(fig=fig)
測定値と予測値のグラフ
ローカルレベルで予測した部分のグラフ
完全に残差＋定数項のようになってます。
季節成分と自己回帰で変動要因をきれいに抜けた、、と言う事になるのでしょうか。
季節成分で予測した部分のグラフ
自己回帰で予測した部分のグラフ
分析結果サマリ
分析結果サマリのcoef（係数）を見ると、
sigma2.level（ローカルレベル）、sigma2.seasonal（季節成分）、ar.L2（t-2の回帰成分）がかなり低く、重回帰分析の知識からすると不要に思えますが、、
どれを消してもおかしな予測結果になりました。
P値、検定結果もよろしくはないですが、、現在の知識ではこの辺りが限界です。
最後に予測値をのせておきます。
いままでで一番それっぽくはなってますが、いまはまだこの結果の解釈ができないので、まだまだ勉強が必要です。。
