前回は、トレンド + 周期性 + その他のその他部分から祝日効果を取り出しました。
今回は、ここからさらに他の時系列データから推定し、残りをノイズとすることにします。
コード
他の時系列データは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()
思っていた形とかなりちがいますね。
元のデータが上昇する形になっているため、トレンドに吸収されているのかもしれません。
pred['trend'].plot() (pred['trend'] + (pred['x'] * pred['trend'])).plot()
xの補正を組み合わせるとよりシグモイドに近い形になっています。
今回は直接元のyは推定できず、トレンド項の補正にとどまった模様です。他の変数からの影響は、やはり、回帰を使ったり、ラグがとれるVARなどを使った方が無難そうです。
その他の残りのノイズを確認する
最後に、ノイズに傾向がないか把握します。
ノイズは想定しているモデルの値と実測の差なので、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)
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
祝日やググっと上がるところにノイズの傾向が残っているので、このあたりもまだ改善できそうです。
それ以外の日付は実務では、そのあたりのデータを掘ることになりそうです。