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()

思っていた形とかなりちがいますね。

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

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

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)

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

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

pythonでProphetを使いたい(5:祝日効果を追加したい)

周期性の設定まで完了したため、その他(ノイズなど）の分析を行っていきます。 まず、あまりフィットしていない日付を見てみます。

future = mdl.make_future_dataframe(periods=0,freq='D')
future['cap'] = 6
pred  = mdl.predict(future)

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

#はずれが大きい順に日をならべる
pd.Series(not_fitted_dt).apply(jpholiday.is_holiday_name)

0         憲法記念日
1          春分の日
2         みどりの日
3           海の日
4        勤労感謝の日
5          敬老の日
6          成人の日
7    こどもの日 振替休日
8         こどもの日
9       即位礼正殿の儀

きちんと祝日たちでした。
単純にこの差を祝日効果にしてもいいかもしれませんが、せっかくなのでモデルに組み込んでみます。

コード

まずはProphetの祝日設定用のデータフレームを作成します。
属性として祝日名と日付を持つ必要があります。

dt = pd.Series(pd.date_range(start='2019-01-01',periods=N))
hols = dt.apply(jpholiday.is_holiday_name )

dts = pd.DataFrame(
    {
      'holiday':dt.apply(jpholiday.is_holiday_name),
      'ds':dt.values,
      'lower_window': 0,
      'upper_window': 1,  
    }
)

あとはモデルと予測用のデータにholidaysにこのデータを渡します。

フィッティング

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=1.0 #祝日効果の強さをいじるパラメータ
)
mdl_df['cap'] = 6
#fitの前
mdl.add_seasonality(name = '_25',period=25,fourier_order=10)
mdl.fit(
    mdl_df
)
future = mdl.make_future_dataframe(periods=0,freq='D')
future['cap'] = 6
future['holiday'] = hols.values
pred  = mdl.predict(future)

#可視化
mdl_df.y.plot()
pred.yhat.plot()

祝日時に同じように上がる用になりました。 各祝日の影響を見ておきます。

import japanize_matplotlib

#祝日情報
mdl.holidays
hol_name = mdl.holidays['holiday']
hol_e = pred[hol_name.to_list()]
hol_e.max().plot.barh()

※なお、曜日効果や祝日効果をすべて合わせたものがadditive_terms,multiplicative_termsになります。

PythonでProphetをつかいたい(4:周期性に独自のものを追加)

時系列=トレンド + 周期性 + その他のうちの周期性に独自の周期性を追加していきます。

▼こちらで作ったデータ、モデルをさらにチューニングしていきます。 esu-ko.hatenablog.com

前回は、曜日の影響を確認しました。

現状のモデルでは、長いトレンドと、短期七日の周期性はとらえていますが、その他の周期性が把握できていません。

新しい季節性の追加

まず、acf,pacfを見て周期にあたりをつけます。

import statsmodels.api as sm

sm.graphics.tsa.plot_acf(mdl_df['y'], lags=50)
sm.graphics.tsa.plot_pacf(mdl_df['y'], lags=50)

ざっくりいってしまうと、青枠から出ているところが自己相関がありそうなところです。

• 10以下のところはありそう => トレンドと曜日(7日)
• 25前後と、50あたりがありそう=>25の倍数

後者を追加します。

周期性の追加方法

Prophetではadd_seasonalityで周期性を追加できます。

mdl = Prophet(
    yearly_seasonality=False,
    weekly_seasonality=True,
    daily_seasonality=False,
    growth = 'logistic',
    changepoint_prior_scale=0.3,
    seasonality_mode='multiplicative'
)
#fitの前に記載
#名前、周期、大体の範囲、トレンドの影響を受けるか(今回はうけない)
mdl.add_seasonality(name = '_25',period=25,fourier_order=10,mode='additive')
mdl.fit(
    mdl_df
)
future = mdl.make_future_dataframe(periods=0,freq='D')
future['cap'] = 6
pred  = mdl.predict(future)

#可視化
ts.plot()
pred.yhat.plot()

だいぶ形が近くなってきました。

今回設定した周期性の確認

周期性は、nameで指定したカラム名でpredの中に格納されています。

red._25.plot()
df['trend2'].plot()

かなり、近い形で周期性を推定しています。

現在のモデルで予測をしてみる

future = mdl.make_future_dataframe(periods=365,freq='D')
future['cap'] = 6
pred  = mdl.predict(future)

ts.plot()
pred.yhat.plot()

最初のモデルと比べると、かなり納得しやすい形になりつつあります。