状態空間モデル
時系列を解析してますか。
何故かRの情報が多い気がする状態空間モデルですが、Pythonでも解析できます。
statsmodelsというモジュールを使います。ARIMAなどの時系列分析の有名なモデルも使うことができるのですが、今回はUnobservedComponentsという状態空間モデルの一番簡単なモデルを利用したいと思います。
UnobservedComponentsという名前は構造時系列モデルの書籍を書いたHarvey A. C. (1989). Forecasting, Structural Time Series Models and the Kalman Filterの表現から取られている気がします。
環境
- Google Colaboratory（最近敷居が高くなった）
下準備
kaggleのデータをお借りすることにします。
2 Years Restaurant Sale With Multiple External Var
2018年から2020年におけるブラジル南部のレストランでのランチの売上個数だそうです。2019年末までのデータを使おうと思います（下記は一部）。
| DATE | SALES | IS_WEEKEND | IS_HOLIDAY | IS_FESTIVE_DATE | AMOUNT_OTHER_OPENED_RESTAURANTS |
|---|---|---|---|---|---|
| 20180214 | 132 | 0 | 0 | 0 | 7 |
| 20180215 | 149 | 0 | 0 | 0 | 7 |
| 20180216 | 130 | 0 | 0 | 0 | 7 |
| 20180217 | 174 | 1 | 0 | 0 | 7 |
| 20180218 | 185 | 1 | 0 | 0 | 7 |
データフレームのインデックスを日付 (DatetimeIndex) として、頻度 (freq) を設定すると、学習データから期間外の予測するときにスムーズに進んだりして便利なので、そうしたいと思います。
日付と日付の間に空きがあると頻度を日次 (D) に設定できないので、目的変数の中身はNaNでいいので、連続した日付でデータフレームを作ってマージし、これをインデックスにします。
time_df = pd.DataFrame(pd.date_range(data.DATE.min(), datetime.strptime("20191231", "%Y%m%d"), freq='d'),
columns=["DATE"])
data = pd.merge(time_df, data, how="left", on="DATE")
data = data.set_index("DATE").sort_index()
data.index.freq = "d"
お祭りの日（ブラジルだからかな？）、周囲の開いているレストランの数などのデータがあったので、あとで外生変数として加えようと思います。目的変数の中身はNaNでいいんですが外生変数は怒られるので前後の値で補完しておきます。
data.IS_FESTIVE_DATE = data.IS_FESTIVE_DATE.interpolate('ffill')
data.AMOUNT_OTHER_OPENED_RESTAURANTS = data.AMOUNT_OTHER_OPENED_RESTAURANTS.interpolate('ffill')
data.IS_FESTIVE_DATE = data.IS_FESTIVE_DATE.interpolate('bfill')
data.AMOUNT_OTHER_OPENED_RESTAURANTS = data.AMOUNT_OTHER_OPENED_RESTAURANTS.interpolate('bfill')
モデル作成
モデル概形
level=”Local linear trend” の部分で、どこに時間変化を入れるかなどのモデルの概形が指定できます。今回はローカル線形トレンドモデルに決め打ちします。
model = sm.tsa.UnobservedComponents(
endog=data[["SALES"]],
level='local linear trend',
seasonal=7,
freq_seasonal=[{'period':365, "harmonics": 10}],
exog=data[["IS_FESTIVE_DATE", "AMOUNT_OTHER_OPENED_RESTAURANTS"]]
)
季節と曜日の効果
季節や曜日のような周期的なデータは、ダミー変数で入れる方法と三角関数で入れる方法があります。
参考書籍：カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2) 野村 俊一
ダミー変数として入れる場合は、seasonalに周期数を入力します。今回は曜日によって売上が異なると仮定して、seasonal=7を設定します。例えばこれが月次のデータなら、seasonal=12とすると月をダミー変数として扱えます。
三角関数として入れる場合は、freq_easonalに {”period”: 周期数, “harmonics”: 周波数成分の個数} とした辞書のリストを指定します。公式に丁寧な説明のページがあります。周波数成分の数は周期の半分まで設定できますが、増やすと推定にかかる時間は伸びます。
参考公式ページ：Seasonality in time series data
外生変数
複数指定できます。
ところでこの外生変数は時間によって変化せず、目的変数に対して切片のように加わるだけです。これを時間で変化させるにはUnobservedComponentsでなく、MLEmodelというモデルが用意されています。
MLEmodelでは他にも目的変数を複数用意するなど多彩なモデルづくりが可能です。
参考公式ページ：Custom statespace models
モデル推定
Nelder-Meadでパラメーター推定してからBFGSで推定しておきます（参考URL）。
max-iterを大きくしておいて推定が収束するように祈ります（祈るだけではダメです）。
Unobserved Components Results
===============================================================================================
Dep. Variable: SALES No. Observations: 728
Model: local linear trend Log Likelihood -3162.793
+ stochastic seasonal(7) AIC 6339.586
+ stochastic freq_seasonal(365(10)) BIC 6371.444
Date: Wed, 12 Oct 2022 HQIC 6351.901
Time: 00:58:00
Sample: 01-03-2018
- 12-31-2019
Covariance Type: opg
========================================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------------------------------
sigma2.irregular 542.2542 26.031 20.831 0.000 491.234 593.275
sigma2.level 1.834e-13 0.363 5.05e-13 1.000 -0.712 0.712
sigma2.trend 1.85e-18 8.32e-06 2.22e-13 1.000 -1.63e-05 1.63e-05
sigma2.seasonal 3.353e-11 0.147 2.29e-10 1.000 -0.287 0.287
sigma2.freq_seasonal_365(10) 0.0024 0.010 0.227 0.820 -0.018 0.023
beta.IS_FESTIVE_DATE 42.5521 14.901 2.856 0.004 13.347 71.757
beta.AMOUNT_OTHER_OPENED_RESTAURANTS -5.6407 1.467 -3.844 0.000 -8.517 -2.764
===================================================================================
Ljung-Box (L1) (Q): 20.72 Jarque-Bera (JB): 31.20
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.95 Skew: -0.04
Prob(H) (two-sided): 0.71 Kurtosis: 4.03
===================================================================================
できました。時系列の最初の方が推定ができていないために全体の図示がおかしいですが許してください。
実際にはこの結果を見て、AICを確認しつつモデルの変数を増やしたり減らしたり、ローカル線形トレンド以外のモデルに変更したりするかと思いますが今回はこのまま使います。
モデル利用
forecast, predict, get_forecast, get_predictionなどのメソッドがあるので（それぞれやる仕事や返り値が微妙に違います）、これで将来の予測を行います。未来の外生変数は用意しておかなくてはいけません。
steps = 60
fes_future = [0]*steps
restaurants_future = [0]*steps
forecast = result.forecast(steps=steps,
exog=[fes_future, restaurants_future])
forecast.head()
| 2020-01-01 | 197.202209 |
| 2020-01-02 | 191.954819 |
| 2020-01-03 | 185.610514 |
| 2020-01-04 | 188.698018 |
| 2020-01-05 | 189.558425 |
また、「お祭りが存在しなかったらどうなっていたか」など、外生変数のみ異なる値を利用して推定を行いたい場合は、applyメソッドで推定結果を出すことができます。
data_apply = data.copy()
data_apply[["IS_FESTIVE_DATE"]] = 0
result_apply = result.apply(endog=data_apply[["SALES"]],
exog=data_apply[["IS_FESTIVE_DATE", "AMOUNT_OTHER_OPENED_RESTAURANTS"]]
)
prediction = result.predict(start=datetime.strptime("20190101", "%Y%m%d"),
end=datetime.strptime("20191231", "%Y%m%d"))
prediction_apply = result_apply.predict(start=datetime.strptime("20190101", "%Y%m%d"),
end=datetime.strptime("20191231", "%Y%m%d"))
print(prediction.sum() - prediction_apply.sum())
# 42.287497166733374
2019年には11月と12月の2回お祭りがあったようですが、これによって2019年のランチは42個多く売れたみたいですね。毎日お祭りをやった方がいいです。
