あれもPython,これもPython

Pythonメモ※本サイトはアフィリエイトを利用しています

pythonでProphetを使いたい(6:ノイズと別の時系列データ)

前回は、トレンド + 周期性 + その他のその他部分から祝日効果を取り出しました。
今回は、ここからさらに他の時系列データから推定し、残りをノイズとすることにします。

コード

他の時系列データはadd_regressorでモデルに追加できます。これもmodeを足すことができます。
観測値自体は、モデルに食わせるデータに直接いれます。 実務では、予測に使うにはx自体が予測できないといけないので注意が必要です。

mdl = Prophet(
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=False,
    growth = 'logistic',
    changepoint_prior_scale=0.3,
    holidays = dts[dts['holiday'].notna()],
    holidays_prior_scale=0.05 
)
mdl_df['cap'] = 6
# 観測された別のデータ
mdl_df['x'] = x
#fitの前
mdl.add_seasonality(name = '_25',period=25,fourier_order=10)
#説明変数として追加
mdl.add_regressor('x',mode='multiplicative')
mdl.fit(
    mdl_df
)
future = mdl.make_future_dataframe(periods=0,freq='D')
future['cap'] = 6
future['holiday'] = hols.values
future['x'] = x
pred  = mdl.predict(future)

確認してみる

#単体
pred['x'].plot()

#トレンドを乗せたものを、実際の値と比較
(pred['x']*pred['trend']).plot()
df['y'].plot()

f:id:esu-ko:20200913172806p:plain

f:id:esu-ko:20200913172922p:plain 思っていた形とかなりちがいますね。

元のデータが上昇する形になっているため、トレンドに吸収されているのかもしれません。

pred['trend'].plot()
(pred['trend'] + (pred['x'] * pred['trend'])).plot()

f:id:esu-ko:20200913173014p:plain

xの補正を組み合わせるとよりシグモイドに近い形になっています。

今回は直接元のyは推定できず、トレンド項の補正にとどまった模様です。他の変数からの影響は、やはり、回帰を使ったり、ラグがとれるVARなどを使った方が無難そうです。

esu-ko.hatenablog.com

その他の残りのノイズを確認する

最後に、ノイズに傾向がないか把握します。

ノイズは想定しているモデルの値と実測の差なので、yhat-yで算出可能です。

ノイズの自己相関の確認

import statsmodels.api as sm
sm.graphics.tsa.plot_acf((ts - pred['yhat']), lags=50)

sm.graphics.tsa.plot_pacf((ts - pred['yhat']), lags=50)

f:id:esu-ko:20200913173454p:plain f:id:esu-ko:20200913173507p:plain

7の倍数が飛び出しているので、まだ曜日周期は改善できるようです。

ノイズの多い日とノイズにトレンドが残っていないか確認

tmp = pred[['ds','yhat']].assign(y=ts).set_index('ds')
(tmp.y-tmp.yhat).sort_values(ascending=False).head(10)

(ts - pred['yhat']).rolling(10).mean().plot()
2019-11-23    4.004385
2019-02-11    3.834714
2019-05-06    2.668610
2019-10-14    2.261944
2019-09-23    1.919278
2019-12-28    1.918154
2019-11-24    1.900130
2019-12-01    1.792548
2019-05-18    1.761975
2019-09-14    1.730068

f:id:esu-ko:20200913173635p:plain

祝日やググっと上がるところにノイズの傾向が残っているので、このあたりもまだ改善できそうです。
それ以外の日付は実務では、そのあたりのデータを掘ることになりそうです。