從時間序列ariam模型到Prophet模型的應用

所遇問題

前兩個星期，我的一個朋友(任職於上海的一家數據公司的機器學習崗位)在工作中有一個華為手機各價位銷售趨勢預測的項目，需要用到時間序列，找到我一起做一下協助。數據是2016年12月-2019年12月的脫敏數據，拿到的數據很簡單，我們的思路是按照價格段來分組，根據不同的價格段用前三年的數據劃分訓練集和測試集，預測2020年前6個月的趨勢走向。但是遇到的問題是數據中2017年整個6月份的數據是缺失的，各年的3月，5月，11月分別有零散缺失。其中每年份中6月18日，11月11日，12月12日等特殊日期的峰值特別的高。

了解完數據的特點，我首先想到的是用時間序列裡的arima模型，把2017年1月，2018年1月的數據作為訓練集，2019年1月的數據作為測試集，跑入模型，預測2020年1月的數據，同理用2017年2月，2018年2月的數據作為訓練集，2019年2月的數據作為測試集，預測2020年2月的數據，用這個思路分別來預測2020年1月，2月，3月，4月，5月，6月的數據，拋開缺失不談，這種思路預測效果偏差很大，直接原因是因為訓練數據量太少。第二種做法是直接把三年的數據全放進去，按七三比例拆分成訓練集和測試集，利用步長差分消除季節性去跑模型，一起預測出2020年前六個月的數據。但是arima模型不能很好的解決2017年6月份缺失的數據，針對節假日數據也沒辦法很好的控制，我們採取的方法是把屬於6月18日，11月11日，12月12日等購物節假日的數據剔除，零散缺失值用-1填充，但是最後結果也是不盡人意。

探索新方法

通過查找資料和交流，最終確定了處理時間序列的另外一種非常實用的模型Prophet模型。Prophet模型可以很好的解決異常的節假日級別的數據，在這裡使數據和代碼簡單的介紹一下Prophet模型的使用。

Prophet 遵循 sklearn 庫建模的應用程式接口。我們創建了一個 Prophet 類的實例，其中使用了「擬合模型」 fit 和「預測」 predict 方法。

# Python import pandas as pdimport numpy as npfrom fbprophet import Prophetimport matplotlib.pyplot as pltimport warningswarnings.filterwarnings("ignore")data = pd.read_csv('原版市場價格段.csv',encoding='gbk')data['y'] = np.log(data['1200元~2400元'])data['ds'] = data['時間']df = data[['ds','y']]df# 擬合模型m = Prophet()m.fit(df)# 構建待預測日期數據框，periods = 365 代表除歷史數據的日期外再往後推 365 天future = m.make_future_dataframe(periods=365)future.tail()# 預測數據集forecast = m.predict(future)forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()# 展示預測結果m.plot(forecast)# 預測的成分分析繪圖，展示預測中的趨勢、周效應和年度效應m.plot_components(forecast)

forecast[['ds','yhat_1']]forecast['y1'] = data['1200元~2400元']forecast[['ds','yhat_1','y1']].describe(percentiles= [0.1,0.9,0.95,0.99,0.999])forecast.loc[forecast['y1'] > 3.492611e+07,'y1'] = 8.811842e+06forecast[['ds', 'yhat']] forecast['yhat_1'] = np.exp(forecast[ 'yhat'])forecast[['ds', 'yhat','yhat_1']]plt.figure(figsize=(15,6))plt.plot(forecast['ds'],forecast[['yhat_1','y1']])# plt.xticks(np.arange(57)[0::6], train.時間[0::6])#plt.xticks(np.arange(57)[0::12], train.時間[0::12],rotation=90)plt.xlabel('時間')plt.ylabel('1200元~2400元')plt.title('手機1200元~2400元價格段')plt.show()

飽和預測

data = pd.read_csv('原版市場價格段.csv',encoding='gbk')data['y'] = np.log(data['1200元~2400元'])data['ds'] = data['時間']df = data[['ds','y']]df['cap'] = 1.5dfm = Prophet(growth='logistic')m.fit(df)future = m.make_future_dataframe(periods=365)future['cap'] = 1.5fcst = m.predict(future)fig = m.plot(fcst)

df['y'] = 10 - df['y']df['cap'] = 6df['floor'] = 1.5future['cap'] = 6future['floor'] = 1.5m = Prophet(growth='logistic')m.fit(df)fcst = m.predict(future)fig = m.plot(fcst)

趨勢突變點

在之前的部分，我們可以發現真實的時間序列數據往往在趨勢中存在一些突變點。默認情況下， Prophet 將自動監測到這些點，並對趨勢做適當地調整。不過，要是對趨勢建模時發生了一些問題，例如：Prophet 不小心忽略了一個趨勢速率的變化或者對歷史數據趨勢變化存在過擬合現象。如果我們希望對趨勢的調整過程做更好地控制的話，那麼下面將會介紹幾種可以使用的方法。

Prophet 中的自動監測突變點

data = pd.read_csv('原版市場價格段.csv',encoding='gbk')data['y'] = np.log(data['1200元~2400元'])data['ds'] = data['時間']df = data[['ds','y']]# 擬合模型m = Prophet()m.fit(df)# 構建待預測日期數據框，periods = 365 代表除歷史數據的日期外再往後推 365 天future = m.make_future_dataframe(periods=365)future.tail()# 預測數據集forecast = m.predict(future)forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()from fbprophet.plot import add_changepoints_to_plotfig = m.plot(forecast)a = add_changepoints_to_plot(fig.gca(), m, forecast)

調整趨勢的靈活性

如果趨勢的變化被過度擬合（即過於靈活）或者擬合不足（即靈活性不夠），可以利用輸入參數 changepoint_prior_scale 來調整稀疏先驗的程度。默認下，這個參數被指定為 0.05 。增加這個值，會導致趨勢擬合得更加靈活。代碼和圖如下所示：

data = pd.read_csv('原版市場價格段.csv',encoding='gbk')

data['y'] = np.log(data['1200元~2400元'])

data['ds'] = data['時間']df = data[['ds','y']]# 擬合模型

m = Prophet(changepoint_prior_scale=0.5)m.fit(df)# 構建待預測日期數據框，

periods = 365 代表除歷史數據的日期外再往後推 365 天

future = m.make_future_dataframe(periods=365)

forecast = m.predict(future)fig = m.plot(forecast)

指定突變點的位置

如果你希望手動指定潛在突變點的位置而不是利用自動的突變點監測，可以使用 changepoints 參數。

m = Prophet(changepoints=['2017-06-18', '2017-11-11', '2017-12-12' ])m.fit(df)future = m.make_future_dataframe(periods=365)forecast = m.predict(future)m.plot(forecast)

季節性，假期效果和回歸量

對假期和特徵事件建模

playoffs = pd.DataFrame({ 'holiday': 'playoff', 'ds': pd.to_datetime(['2016-12-12','2017-06-18', '2017-11-11', '2017-12-12', '2018-06-18', '2018-11-11', '2018-12-12', '2019-06-18', '2019-11-11', '2019-12-12', '2020-06-18', '2020-11-11', '2020-12-12']), 'lower_window': 0, 'upper_window': 1,})superbowls = pd.DataFrame({ 'holiday': 'superbowl', 'ds': pd.to_datetime(['2017-11-11', '2018-11-11', '2019-11-11', '2020-11-11']), 'lower_window': 0, 'upper_window': 1,})holidays = pd.concat((playoffs, superbowls))m = Prophet(changepoint_prior_scale=0.005,holidays=holidays,yearly_seasonality=80)m.fit(df)future = m.make_future_dataframe(periods=365)forecast = m.predict(future)# 可通過 forecast 數據框，來展示節假日效應：# 看一下假期的最後10行數據forecast[(forecast['playoff'] + forecast['superbowl']).abs() > 0][ ['ds', 'playoff', 'superbowl']][-10:]# 在成分分析的圖中，如下所示，也可以看到節假日效應。我們可以發現，在節日日期附近有一個穿透，而在超級碗日期時穿透則更為明顯。 fig = m.plot_components(forecast)m.plot(forecast)forecast['y'] = df['y']forecast[['ds', 'yhat']]forecast[['ds', 'yhat','y']]plt.figure(figsize=(15,6))plt.plot(forecast['ds'],forecast[['yhat','y']])# plt.xticks(np.arange(57)[0::6], train.時間[0::6])#plt.xticks(np.arange(57)[0::12], train.時間[0::12],rotation=90)plt.xlabel('時間')plt.ylabel('1200元~2400元')plt.title('手機1200元~2400元價格段')plt.show()

季節性的傅立葉級數

季節性是用部分傅立葉和估計的。有關完整的細節，請參閱論文，以及維基百科上的這個圖，以說明部分傅立葉和如何近似於一個線性周期信號。部分和(order)中的項數是一個參數，它決定了季節性的變化有多快。為了說明這一點， 我們仍似乎用第一部分中佩頓 · 曼寧的數據。每年季節性的默認傅立葉級數是10，這就產生了這樣的擬合:

from fbprophet.plot import plot_yearlym = Prophet().fit(df)a = plot_yearly(m)

默認值10通常是合適的，但是當季節性需要適應更高頻率的變化時，它們可以增加，並且通常不那麼平滑。在實例化模型時，可以為每個內置季節性指定傅立葉級數，這裡增加到20:

from fbprophet.plot import plot_yearlym = Prophet(yearly_seasonality=20).fit(df)a = plot_yearly(m)

可以看到，曲線更加的多變了。增加傅立葉項的數量可以使季節性適應更快的變化周期，但也可能導致過度擬合:N個傅立葉項對應於用於建模周期的2N個變量。