時間序列預測（一）：趨勢分解法+ARIMA
背景介紹：最近剛接手的項目是時間序列預測，我分別用3中方法進行了實現，分別是
時間序列預測方法一：趨勢分解法+ARIMA
後續還會更新：
時間序列預測方法二：Facebook Prophet模型方法
時間序列預測方法三：基於LTSM神經網路搭建時間預測模型
並會對其三種方法進行性能優略的詳細比較。。。。敬請期待
項目描述
項目是一個典型的時間序列預測。
我的數據樣本是近兩年的活躍數據，其中分為三列，分別是登錄量、行為量和註冊量。三類數據都和時間密切相關，有著明顯的week周期性、year周期特性和明顯的節假日帶來的拐點數據。目標就是通過過去的數據對未來數據的預測。
而且要準確預測出假期的特性、week周期性和year周期性。
工具：Anaconda
圖表 1 用戶行為數據
Load和準備數據
數據準備的過程中導入必要的庫函數，以及此種時間序列預測需要的包statsmodels模型中用於進行平穩性檢測（AD檢測）的adfuller、用於進行模型分解的seasonal_decompose，以及用於分解後的部分進行建模的ARIMA模型，以及確認ARIMA模型參數的自相關函數和偏自相關函數acf,pacf。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import os
import warnings
warnings.filterwarnings(ignore)
from datetime import datetime
from matplotlib.pyplot import rcParams
rcParams[figure.figsize]=15,6
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.arima_model import ARIMA
from datetime import timedelta
from statsmodels.tsa.stattools import acf,pacf
接下來就是導入業務數據，經排查，2016年9月23日，2016年9月24日，2016年9月25是周五、周六、周日，後台系統在上線維護，所以將數據進行恢復處理ts3註冊數據在2016年9月23日，數據為0點，所以選擇使用2016年9月22日數據進行填充ts2投遞數據和ts1登錄數據，2016年9月24日，2016年9月25選擇使用一周前20160917、20160918的數據進行替代。到此數據異常點處理完畢。
#import the data
dataparse=lambda dates : datetime.strptime(dates,%Y-%m-%d)
data=pd.read_excel(C:/Users/xxxxx/Desktop/時間序列預測/acitivedata.xlsx,pares_dates=[Day],index_col=Day,data_parser=dataparse)
ts1=data[login]
ts2=data[delivery]
ts3=data[registration]
ts1[2016-09-24]=ts1[2016-09-17]
ts1[2016-09-25]=ts1[2016-09-18]
ts2[2016-09-24]=ts2[2016-09-17]
ts2[2016-09-25]=ts2[2016-09-18]
演算法原理簡介
時間序列趨勢分解法原理
時間序列分解法（Time-series Decomposition）是數年來一直非常有用的方法,這種方法包括譜分析、時間序列分析和傅立葉級數分析等。
時間序列分解模型
時間序列y可以表示為以上四個因素的函數，即：
時間序列分解的方法有很多，較常用的模型有加法模型和乘法模型。本文的演算法模型是加法模型。
加法模型為：
其中，公式中代表的含義，是將時間序列分別分解為長期時間趨勢T、季節性時間趨勢S，和周期性時間趨勢C以及剩餘的殘差項。
本項目中，趨勢分解出趨勢部分、周期部分和殘差項部分
其中趨勢部分經過預處理（取對數和差分的方法），再次通過ARIMA模型進行建模模擬，最後預測數據由趨勢部分的建模結果累加周期部分。
ARIMA模型原理
ARIMA (Auto Regressive Integrated Moving Average) 可以用來對時間序列進行預測，常被用於需求預測和規劃中。
可以用來對付 『隨機過程的特徵隨著時間變化而非固定』 且 『導致時間序列非平穩的原因是隨機而非確定』 的問題。不過，如果是從一個非平穩的時間序列開始， 首先需要做差分，直到得到一個平穩的序列。
模型的思想就是從歷史的數據中學習到隨時間變化的模式，學到了就用這個規律去預測未來。
ARIMA(p,d,q)模型，其中 d 是差分的階數，用來得到平穩序列。
AR是自回歸, p為相應的自回歸項。
MA為移動平均，q為相應的移動平均項數
ARIMA（p，d，q）模型是ARMA（p，q）模型的擴展。
ARIMA（p，d，q）模型可以表示為：
AR：當前值只是過去值的加權求和。
MA：過去的白噪音的移動平均。
ARMA：AR和MA的綜合。
ARIMA：和ARMA的區別，就是公式左邊的x變成差分運算元，保證數據的穩定性。
差分運算元就是：
令 wt 為：
則 ARIMA 就可以寫成：
網上關於ARIMA原理眾多，不在累述，作者在這裡只想強調，最後ARIMA是時間序列的預測精準的前提是首先時間序列是平穩序列，本模型中使用常規的ADfuller檢測以及移動均值和方差進行數據平穩化的驗證，將非平穩序列進行平穩性的轉化通常是對序列進行取對數、差分、移動平均或其一或組合使用。
模型模擬
# take aDFuller test the statioinarity of the timeseries
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean=pd.Series.rolling(timeseries,window=12).mean()
#rolmean=pd.rolling_mean(timeseries,window=12) 因為後續pandas因為安裝了最新的版本
#rolstd=pd.rolling_std(timeseries,window=12)
rolstd=pd.Series.rolling(timeseries,window=12).std()
#plot rolling statistics
orig=plt.plot(timeseries,color=blue,label=Original)
mean=plt.plot(rolmean,color=red,label=Rolling mean)
std=plt.plot(rolstd,color=black,label=Rolling std)
plt.legend(loc=best)
plt.title(Rolling mean & Standard Deviatiom)
plt.show(block=False)
#perfprmance Dickey_Fuller Test
print(Resluts of Dickey-Fuller Test:)
dftest=adfuller(timeseries,autolag=AIC)
dfountput=pd.Series(dftest[0:4],index=[Test Statistic,p-value,#lags used,Number of obervations used])
for key,value in dftest[4].items():
dfountput[Critical Value(%s)%key]=value
print(dfountput)
#找出ARIMA(p,I,q)參數
def show_acf_pacf(data):
lag_acf=acf(data,nlags=20)
laf_pacf=pacf(data,nlags=20,method=ols)
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestylex=--,color=gray)
plt.axhline(y=-1.96/np.sqrt(len(data)),linestylex=--,color=gray)#上下邊間是屬於置信區間
plt.axhline(y=1.96/np.sqrt(len(data)),linestylex=--,color=gray)#上下邊間是屬於置信區間
plt.title(Autocorrelation Function)
#Plot PACF:
plt.subplot(122)
plt.plot(laf_pacf)
plt.axhline(y=0,linestylex=--,color=gray)
plt.axhline(y=-1.96/np.sqrt(len(data)),linestylex=--,color=gray)
plt.axhline(y=1.96/np.sqrt(len(data)),linestylex=--,color=gray)
plt.title(Partial Autocorrelation Function)
plt.tight_layout()
#這個是ARIMA模型
from statsmodels.tsa.arima_model import ARIMA
def arima_performance(data,order1):
model=ARIMA(data,order=order1)
results_arima=model.fit(disp=-1)
results_arima_value=results_arima.fittedvalues
result_future=results_arima.forecast(7)
return results_arima_value,result_future
#請注意，這裡面的數據如果是差分輸入，那麼請輸入差分序列
def arima_plot(data,results_arima_value):
plt.plot(data)
plt.plot(results_arima_value,color=red)
plt.title(RSS:%.4f%sum((results_arima_value-data)**2))
def add_season(ts_recover_trend,startdate):
ts2_season=ts2_seasonal
values=[]
low_conf_values=[]
high_conf_values=[]
for i,t in enumerate(ts_recover_trend):
trend_part=ts_recover_trend[i]
#begin=datetime(2016,8,9)
begin=startdate
day=begin+timedelta(days=i)
#添加周期部分
season_part=ts2_season[ts2_season.index.dayofweek==day.weekday()].mean()
#趨勢+周期+誤差界限
#ts2_low_error,ts2_high_error
predict=trend_part+season_part
low_bound=trend_part+season_part+ts2_low_error
high_bound=trend_part+season_part+ts2_high_error
values.append(predict)
low_conf_values.append(low_bound)
high_conf_values.append(high_bound)
final_pred=pd.Series(values,index=ts_recover_trend.index,name=predict)
low_conf=pd.Series(low_conf_values,index=ts_recover_trend.index,name=low_conf)
high_conf=pd.Series(high_conf_values,index=ts_recover_trend.index,name=high_conf)
return final_pred,low_conf,high_conf
數據樣本的單位在幾十萬級別，首先對數據進行取對數
"""
先用對數變換，對數據進行取log變換
"""
ts1= np.log(ts1)
ts2= np.log(ts2)
ts3= np.log(ts3)
def model_decompose(ts):
decomposition=seasonal_decompose(ts,freq=7)
ts_trend=decomposition.trend
ts_seasonal=decomposition.seasonal
ts_resid=decomposition.resid
"""
使用分位數進行預測
"""
d=ts_resid.describe()
delta=d[75%]-d[25%]
ts_low_error,ts_high_error= (d[25%] - 1 * delta, d[75%] + 1 * delta)
"""
將分解出來的趨勢數據單獨建模,通過差分的方法對ts2_trend轉化為平穩序列，對ts2_trend部分進行ARIMA建模
"""
#因為分解出來的trend前3項和後3項會有nan,1階查分後第一個數也會是nan值，所以要去nan
ts_trend.count()
ts_trend_diff=ts_trend.diff()
ts_trend_diff.dropna(inplace=True)
return ts_trend_diff,ts_trend,ts_seasonal, ts_low_error,ts_high_error
def data_recovery(ts_results_arima_value,result_future,ts_trend,ts,ts_seasonal, ts_low_error,ts_high_error):
"""
注意：ts2_trend再進行預測的時候，由於前面自身的前3個和尾部3個數據是nan,所以要恢復的數據和要預測的數據化成一致後再進行比對,在恢復的時候應該填充為0
千萬別忘記了。
"""
ts_results_arima_value_cunsum=ts_results_arima_value.cumsum()
ts_recover =pd.Series(ts_trend.iloc[3],index=ts_trend.index)
ts_recover=ts_recover.add(ts_results_arima_value_cunsum,fill_value=0)
ts_trend.fillna(12.617903,inplace=True)
"""
缺少一步，將預測的趨勢數據前6項和後6項進行填充回0
"""
ts_recover_trend=ts_recover
for i in range(0,3):
ts_recover_trend[i]=ts_recover_trend[3]
ts_recover_trend[-i-1]=ts_recover_trend[-4]
"""
在這裡就生成預測將來的數據
"""
predict_future=result_future[0]
predict_future_cumsum=predict_future.cumsum()
#生成預測將來的index
n=7
begin=ts.index[-1]+timedelta(1)
predict_time_index=pd.date_range(start=begin,periods=n,freq=D)
ts_trend_future_recover =pd.Series(ts_recover_trend.iloc[-4],index=predict_time_index)
ts_trend_future__recover1=ts_trend_future_recover.add(predict_future_cumsum)
print("預測將來部分的trend部分：
",ts_trend_future__recover1)
"""
預測將來的數據
ts_trend_future__recover1
"""
startdate=datetime(2018,8,9)
final_pred_fut,low_conf_fut,high_conf_fut = add_season(ts_trend_future__recover1,startdate)
"""
預測原來的數據
"""
startdate=datetime(2016,8,9)
final_pred,low_conf,high_conf= add_season(ts_recover_trend,startdate)
print("RMSE :%.4f"% np.sqrt((sum((ts-final_pred)**2))/len(ts)))
return final_pred,final_pred_fut
ts1_trend_diff,ts1_trend,ts1_seasonal, ts1_low_error,ts1_high_error = model_decompose(ts1)
ts2_trend_diff,ts2_trend,ts2_seasonal, ts2_low_error,ts2_high_error = model_decompose(ts2)
ts3_trend_diff,ts3_trend,ts3_seasonal, ts3_low_error,ts3_high_error = model_decompose(ts3)
test_stationarity(ts1_trend_diff)
以其中的數據ts1登陸數據為例：
test_stationarity(ts1_trend_diff)
show_acf_pacf(ts1_trend_diff)
性能評估
最後也僅以ts1登陸數據為例展示：
ts1_results_arima_value,result_future1=arima_performance(ts1_trend_diff,(1,0,0))
arima_plot(ts1_trend_diff,ts1_results_arima_value)
final_pred1,final_pred_fut1 = data_recovery(ts1_results_arima_value,result_future1,ts1_trend,ts1,ts1_seasonal, ts1_low_error,ts1_high_error)
"""
最後數據的預測輸出
"""
def future_output(final_pred,final_pred_fut,ts):
ts_backlog=np.exp(ts)
final_pred_backlog=np.exp(final_pred)
plt.plot(final_pred_backlog,color=red,label=predict)
plt.plot(ts_backlog,color=blue,label=Orignal)
plt.legend(loc=best)
plt.title(" RMSE decompose+ARIMA model on the test data is :%.4f"% np.sqrt((sum((ts_backlog-final_pred_backlog)**2))/len(ts_backlog)))
plt.savefig("decompose_best.png")
print(" last 145 days decompose+ARIMA model RMSE :%.4f"% np.sqrt((sum((ts_backlog.tail(145)-final_pred_backlog.tail(145))**2))/len(ts_backlog.tail(145))))
plt.savefig(last_145_decompose_best.png)
"""
添加預測出最新一周的數據進行可視化
"""
def final_pred_fut_plot(ts,final_pred,final_pred_fut):
ts_backlog=np.exp(ts)
final_pred_backlog=np.exp(final_pred)
final_pred_fut_log=np.exp(final_pred_fut)
final_pred_backlog1=pd.concat([final_pred_backlog,final_pred_fut_log])
plt.plot(final_pred_backlog1.tail(152),color=red,label=predict)
plt.plot(ts_backlog.tail(145),color=blue,label=Orignal)
plt.legend(loc=best)
RMSE= np.sqrt((sum((ts_backlog.tail(145)-final_pred_backlog.tail(145))**2))/len(ts_backlog.tail(145)))
plt.title(" last 145 days decompose+ARIMA model RMSE :%.4f"% np.sqrt((sum((ts_backlog.tail(145)-final_pred_backlog.tail(145))**2))/len(ts_backlog.tail(145))))
print("最後145天的誤差百分比是：%.4f"% (RMSE/np.mean(ts_backlog.tail(145))))
print("最新一周預測數據是：",final_pred_fut_log)
plt.savefig(oneweek_future.png)
future_output(final_pred1,final_pred_fut1,ts1)
final_pred_fut_plot(ts1,final_pred1,final_pred_fut1)
至此，趨勢分解法預測時間序列完畢。
總結
最後145天的誤差在12%，統計最後100天的誤差就已經小於10%，這個模型也是能夠符合項目要求的，這個模型實現相對簡單，誤差在可接受範圍內，也可以預測長期的趨勢。
推薦閱讀：
