一、 前言

時間序列分析是統計學科的一個重要分支。它主要是通過研究隨著時間的推移事物發展變化過程中的規律，來進行事物未來發展情況的預測。在我們的日常生活中，股票的價格走勢，奶茶店每天的銷售額，一年的降雨量分布，河水四季的漲落情況等都屬于時間序列。時間序列的分析深入諸多行業。

時間序列的分類：

圖1

• 根據指標的平穩性，分為平穩時間序列和非平穩時間序列；
• 根據指標的性質分類，分為總量指標時間序列，相對指標和平均指標時間序列；
• 根據指標的時間屬性分類，分為時期指標時間序列，時點指標時間序列；

時期指標時間序列是可以相加的，并且相加是有意義的，比如每天的訂單量，一個月的訂單量直接將這個月對應的每天的訂單量相加即可。時點指標時間序列是不可以相加的，反映的是某一時間點達到的水平，比如每天庫存量，庫存量相加是沒有統計意義的，每月總庫存量不等于每天庫存量加和。

對于互聯網公司而言，業務量是公司經營關注的重要指標之一。實際情況的復雜性給業務量的分析預測帶來了許多挑戰：

• 具有業務特征的周期性影響
• 節假日等特定時序節點的變異
• 地域差異，空間的相互作用
• 受到庫存、實際市場容量的影響
• 其他外生變量，不可控自然或社會因素

對于時間序列的分析，例如訂單量，話務量，庫存管理等，實現的方式有ANN，RNN，LR，ARIMA，Prophet等，這里我們重點關注ARIMA分析方法。

二、 時間序列分析實踐

2.1 ARIMA模型簡介

ARMA模型的全稱是自回歸移動平均模型，可以說是目前最常用的擬合平穩序列的模型。

ARMA模型由兩部分組成：

p階自回歸模型AR(p)

當φ0=0時，自回歸模型又稱為中心化AR(p)模型。非中心化的AR(p)序列也可以轉化(通過平移)為中心化的AR(p)模型。

AR模型將某時刻t的值用過去若干時刻t-1到t-p的值通過線性組合以及噪聲來表示。

q階移動平均模型MA(q)

當μ=0時，模型MA(q)稱為中心化MA(q)模型，對于非中心化的MA(q)模型只要做簡單的位移就可以轉化為中心化的MA(q)模型。

MA模型是通過歷史點的噪聲線性組合來表示當前時刻的值。

ARMA模型其實就是AR(P)和MA(q)的組合：

同樣的，當φ0=0時該模型稱為中心化的ARMA(p,q)模型。他結合了兩個模型的特點，AR模型處理當前數據與后期數據之間的關系，MA則處理隨機變動的影響。

對于平穩時間序列可以采用ARMA模型直接進行擬合，但是實際場景中，我們的時間序列都是有趨勢的，即一般時序為非平穩的，所以需要做平穩處理，其中最常用的是差分處理，使得時序平穩后進行ARMA分析。

這一過程其實就是ARIMA，在原始非平穩時間序列基礎上做一階或二階差分處理后的平穩時間序列上應用ARMA模型。ARIMA(p,d,q)模型是在ARMA(p,q)兩元組階數基礎上增加差分d的三元組階數模型。

2.2 ARIMA模型實踐分析步驟

圖2

具體實現以python為例。

Step1、讀取時間序列

df = pd.read_csv('testdata.csv', encoding='gbk', index_col='ddate')#時間序列索引轉換為日期格式df.index = pd.to_datetime(df.index)#指標量轉為float類型df['cnt'] = df['cnt'].astype(float)plt.figure(facecolor='white',figsize=(20,8))plt.plot(df.index,df['cnt'],label='Time Series')plt.legend(loc='best')plt.show()?

圖3

Step2、時間序列平穩性檢驗

什么是平穩？

平穩分為嚴平穩和寬平穩，嚴平穩保證時間序列的任何有限維分布對于時間的平移是不變的，比如高斯白噪聲就是嚴平穩序列；寬平穩則要求協方差結構隨時間的平移而不變，或均值和方差是不變的。

為什么需要平穩？

ARIMA包含了AR模型，AR模型的實質是用歷史時間點數據預測當前時間點對應的值。這就要求序列的相關性不會隨著時間變化而變化。

from statsmodels.tsa.stattools import adfullerdef test_stationarity(timeseries): dftest = adfuller(timeseries, autolag='AIC') return dftest[1]

原始時間序列平穩性檢驗未通過(0.94)。從圖3也可以看到，時間序列具有明顯的上升趨勢，所以需要嘗試對時間序列進行差分處理，再次檢驗其平穩性。

Step3、差分處理后檢驗平穩性

pred_day = 7train_start = datetime(2017,3,1)train_end = datetime(2019,8,16)pred_start = train_end+timedelta(1)pred_end = train_end+timedelta(pred_day)train_diff=df[train_start:train_end]train_diff['cnt']=train_diff.diff()print(test_stationarity(train_diff['cnt'][train_start+timedelta(1):train_end]))plt.figure(facecolor='white',figsize=(20,8))plt.plot(train_diff.index,train_diff['cnt'],label='Time Series after diff')plt.legend(loc='best')plt.show()

圖4

差分后的時序平穩性檢驗值9.51*e(-15)，說明差分后時序已經是平穩時間序列了，可以應用ARIMA模型。

Step4、畫出ACF和PACF圖

自相關函數ACF，反映了兩個點之間的相關性。

偏自相關函數PACF則是排除了兩點之間其他點的影響，反應兩點之間的相關性。比如：在AR(2)中，即使y(t-3)沒有直接出現在模型中，但是y(t)和y(t-3)之間也相關。

import statsmodels.api as smfig = plt.figure(figsize=(12,8))ax1 = fig.add_subplot(211)fig = sm.graphics.tsa.plot_acf(train_diff['cnt'][1:], lags=20, ax=ax1)ax1.xaxis.set_ticks_position('bottom')fig.tight_layout()ax2 = fig.add_subplot(212)fig = sm.graphics.tsa.plot_pacf(train_diff['cnt'][1:], lags=20, ax=ax2)ax2.xaxis.set_ticks_position('bottom')fig.tight_layout()plt.show()

圖5

嚴格來看，ACF和PACF顯示存在一定程度的拖尾和振蕩。但是，ACF和PACF在3階后有驟降和平穩的趨勢，考慮到是短期預測的場景，可進一步結合預測效果和模型檢驗來進行判斷。

Step5、ARIMA模型定階

雖然ACF和PACF為我們提供了選擇模型參數的參考依據，但是一般實際情況中，我們總會需要通過模型訓練效果確定最終采用的參數值。在ARMA模型中，我們通常采用AIC法則(赤池信息準則，AIC=2k-2ln(L)，k為模型參數個數，n為樣本數量，L為似然函數)。AIC鼓勵數據擬合的有良性但是盡量避免出現過擬合的情況。所以實際操作中，我們會選擇模型AIC值最小的那組參數。

#定階warnings.filterwarnings("ignore") # specify to ignore warning messagespmax = 8qmax = 8aic_matrix = [] #aic矩陣for p in range(1,pmax+1): tmp = [] for q in range(1,qmax+1): try: #存在部分報錯，所以用try來跳過報錯。 model = ARIMA(endog=df['cnt'],order=(p,1,q)) results = model.fit(disp=-1) tmp.append(results.aic) print('ARIMA p:{} q:{} - AIC:{}'.format(p, q, results.aic)) except: tmp.append(None) aic_matrix.append(tmp)aic_matrix = pd.DataFrame(aic_matrix) #從中可以找出最小值p,q = aic_matrix.stack().idxmin() #先用stack展平，然后用idxmin找出最小值位置。print(u'AIC最小的p值和q值為：%s、%s' %(p+1,q+1))

圖6

由于時間序列為1階差分平穩時間序列，所以模型參數d=1，根據AIC最小原則得到p=7，q=7。

Step6、模型測試與優化

將訓練得到的參數加入模型，分析模型效果。

model = ARIMA(endog=df['cnt'], order=(p,1,q)) #建立ARIMA(7, 1,7)模型result_ARIMA = model.fit(disp=-1,method='css')predict_diff=result_ARIMA.predict()#一階差分還原df_shift=df['cnt'].shift(1)predict=predict_diff+df_shiftplt.figure(figsize=(18,5),facecolor='white')predict[train_start+timedelta(p+1):train_end].plot(color='blue', label='Predict')df['cnt'][train_start+timedelta(p+1):train_end].plot(color='red', label='Original')err=sum(np.sqrt((predict[train_start+timedelta(p+1):train_end]-df['cnt'][train_start+timedelta(p+1):train_end])**2)/df['cnt'][train_start+timedelta(p+1):train_end])/df['cnt'][train_start+timedelta(p+1):train_end].sizeplt.legend(loc='best')plt.title('Error: %.4f'%err) plt.show()

圖7

用訓練好的模型進行未來預測。

y_forecasted =result_ARIMA.forecast(steps=pred_day, alpha=0.01)[0] #作為期7天的預測y_truth = df[pred_start:pred_end]['cnt']# 均方根誤差 #平均錯誤率mse = np.sqrt( ((y_forecasted - y_truth) ** 2) ).mean()error_rate = ( abs(y_forecasted - y_truth)/y_truth ).mean()print('The Mean Error rate of our forecasts is {}'.format(round(error_rate, 4)))?

模型預測誤差8.58%(【偏差/真實值】的均值)，結果并不是太理想，所以我們需要對模型進行優化，考慮是因為指標受到了節假日和周的影響，所以在模型的外生變量里面我們加入節假日和周的識別參數。

加入exog外生變量后，需要重新定階，重新訓練模型，步驟與上類似。優化后的預測誤差1.77%，相比之前有了很大程度的提升。

圖8

Step7、模型檢驗

用模型殘差來檢驗模型的合理性。

resid = result_ARIMA_improve.resid #賦值plt.figure(figsize=(12,8))qqplot(resid,line='q',fit=True)#利用D-W檢驗,檢驗殘差的自相關性print('D-W檢驗值為{}'.format(durbin_watson(resid.values)))

圖9

通過圖9的qq圖可以看出，殘差基本滿足了正態分布。D-W檢驗結果值為1.99，接近于2，說明殘差序列不存在自相關性，即模型較好。

三、總結與展望

• 對于時間序列的分析一定做好前期評估工作，直觀的圖表分析會助力我們的決策。多探索優秀的開源工具庫，往往會使我們事半功倍。

• 模型選擇至關重要，明確模型的適用場景，根據自身的時序選擇適合的模型分析。

• ARIMA模型在短時間內的預期效果還算可以，但是長時間比如未來一年的預測不太適用，因為偏差會逐漸增大。

• 現實中的復雜場景，單一模型很難解決，需要考慮多模型結合的方式實現分析預測。

【作者簡介】June，攜程數據分析經理，對數倉搭建，數據治理，數據分析等方面有較濃厚的興趣。

總結

以上是生活随笔為你收集整理的二阶差分预测后数据还原公式_携程如何基于ARIMA时序分析做业务量的预测的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。