基於AIRMA模型對訂單總額未來七天的預測
? 從時間上看,訂單量時間序列有兩個明顯的特征:
? 1)周期性。每天訂單量的變化趨勢都大致相同,午高峰和晚高峰訂單量集中;
? 2)實時性。當天的訂單量可能會受天氣等因素影響,呈現整體的上漲或下降;
? 預測可以反映未來司機成單的情況,能給運營部門及時調整有效的運營策略。預測又有好幾種方向:基於訂單總額的預測,基於乘客目的地預測,基於上車地點的供需預測等,這裏闡述訂單總額未來七天的預測。
2 數據建模基本步驟
-
獲取被觀測系統時間序列數據;
-
對數據繪圖,觀測是否為平穩時間序列;對於非平穩時間序列要先進行d階差分運算,化為平穩時間序列;
-
經過第二步處理,已經得到平穩時間序列。要對平穩時間序列分別求得其自相關系數ACF 和偏自相關系數PACF ,通過對自相關圖和偏自相關圖的分析,得到最佳的階層 p 和階數 q;
- 由以上得到的d、q、pd、q、p ,得到ARIMA模型。然後開始對得到的模型進行模型檢驗
3 ARIMA實戰解剖
? 原理大概清楚,實踐卻還是會有諸多問題。相比較R語言,Python在做時間序列分析的資料相對少很多。下面就通過Python語言詳細解析後三個步驟的實現過程。
文中使用到這些基礎庫:
pandas,numpy,scipy,matplotlib,statsmodels,pandas, 對其調用如下:
from __future__ import print_function import pandas as pd import numpy as np from scipy import stats import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.graphics.api import qqplot from statsmodels.tsa.stattools import adfuller as ADF
3.1 獲取數據
? 對乘用車日運營報表訂單總額數據,進行分析。
數據如下:
日期 | 訂單總額 |
---|---|
2018/1/1 | 22093.02 |
2018/1/2 | 18917.3 |
2018/1/3 | 6539.3 |
2018/1/4 | 19442.9 |
2018/1/5 | 24425.74 |
2018/1/6 | 29198.1 |
2018/1/7 | 26611.8 |
2018/1/8 | 25389.75 |
2018/1/9 | 22802.62 |
2018/1/10 | 31000.7 |
2018/1/11 | 28794.53 |
2018/1/12 | 28746.15 |
2018/1/13 | 30382 |
2018/1/14 | 31726.1 |
2018/1/15 | 32041.8 |
2018/1/16 | 30514.5 |
2018/1/17 | 31783.9 |
2018/1/18 | 29976.9 |
2018/1/19 | 33292.82 |
2018/1/20 | 35055.3 |
2018/1/21 | 33054.2 |
2018/1/22 | 33000.27 |
2018/1/23 | 34889.85 |
2018/1/24 | 33495.7 |
2018/1/25 | 33960.66 |
2018/1/26 | 35901.61 |
2018/1/27 | 2276.5 |
2018/1/28 | 1268.7 |
2018/1/29 | 32588.14 |
2018/1/30 | 33881 |
2018/1/31 | 29342.02 |
2018/2/1 | 24991.15 |
2018/2/2 | 30995.08 |
2018/2/3 | 31219.9 |
2018/2/4 | 31040.2 |
2018/2/5 | 31152.5 |
2018/2/6 | 31190.76 |
2018/2/7 | 32778.69 |
2018/2/8 | 31938.7 |
2018/2/9 | 31857.8 |
2018/2/10 | 31969.73 |
2018/2/11 | 29406.12 |
2018/2/12 | 26877.6 |
2018/2/13 | 24472.7 |
2018/2/14 | 9663.4 |
2018/2/15 | 142.2 |
2018/2/16 | 108.6 |
2018/2/17 | 321.9 |
2018/2/18 | 16647.8 |
2018/2/19 | 18353.85 |
2018/2/20 | 21131.1 |
2018/2/21 | 20511.5 |
2018/2/22 | 25521.6 |
2018/2/23 | 25701 |
2018/2/24 | 28140.4 |
2018/2/25 | 31535 |
2018/2/26 | 27903.9 |
2018/2/27 | 29730.48 |
2018/2/28 | 28317.9 |
2018/3/1 | 12704.4 |
2018/3/2 | 13736.45 |
2018/3/3 | 13093.1 |
2018/3/4 | 11186.6 |
2018/3/5 | 7279.68 |
2018/3/6 | 9042.7 |
2018/3/7 | 8524.87 |
2018/3/8 | 10346.3 |
2018/3/9 | 9219.9 |
2018/3/10 | 9931.8 |
2018/3/11 | 9772.8 |
2018/3/12 | 9100.1 |
2018/3/13 | 7822.2 |
2018/3/14 | 8863.92 |
2018/3/15 | 20644.8 |
2018/3/16 | 29468.6 |
2018/3/17 | 31152.9 |
2018/3/18 | 29558.9 |
2018/3/19 | 28728.9 |
2018/3/20 | 26359.25 |
2018/3/21 | 27735.05 |
2018/3/22 | 31034.86 |
2018/3/23 | 33682.55 |
2018/3/24 | 38766.6 |
2018/3/25 | 36706.35 |
2018/3/26 | 29303.9 |
2018/3/27 | 26164.2 |
2018/3/28 | 28287.7 |
2018/3/29 | 26742.22 |
2018/3/30 | 30283.25 |
2018/3/31 | 32708.4 |
2018/4/1 | 39168.62 |
2018/4/2 | 40197.3 |
2018/4/3 | 39577.86 |
2018/4/4 | 38930.96 |
2018/4/5 | 33691.44 |
2018/4/6 | 31831.3 |
2018/4/7 | 36942.8 |
2018/4/8 | 35314.52 |
2018/4/9 | 33843.3 |
2018/4/10 | 35063.18 |
2018/4/11 | 35706.36 |
2018/4/12 | 40346.94 |
2018/4/13 | 42388.1 |
2018/4/14 | 41396.29 |
2018/4/15 | 42374.3 |
2018/4/16 | 16489.5 |
2018/4/17 | 31132.9 |
2018/4/18 | 32511 |
2018/4/19 | 35183.22 |
2018/4/20 | 65664.7 |
2018/4/21 | 73029.91 |
2018/4/22 | 73738.58 |
2018/4/23 | 60052.2 |
2018/4/24 | 66994.34 |
2018/4/25 | 70996.55 |
2018/4/26 | 66477.66 |
2018/4/27 | 82832.23 |
2018/4/28 | 88265.1 |
2018/4/29 | 85795.8 |
2018/4/30 | 79327.25 |
2018/5/1 | 87192.7 |
2018/5/2 | 81075.1 |
2018/5/3 | 84212.59 |
2018/5/4 | 98093.35 |
2018/5/5 | 84365.13 |
2018/5/6 | 82469.8 |
2018/5/7 | 83517.96 |
2018/5/8 | 85643.2 |
2018/5/9 | 93718.3 |
2018/5/10 | 94353.3 |
2018/5/11 | 103593.85 |
2018/5/12 | 106824.98 |
2018/5/13 | 107588.58 |
2018/5/14 | 98374.04 |
2018/5/15 | 104670.4 |
2018/5/16 | 118165.44 |
2018/5/17 | 60852.3 |
2018/5/18 | 48933.27 |
2018/5/19 | 49102.2 |
2018/5/20 | 44945.4 |
2018/5/21 | 34141.62 |
2018/5/22 | 32317.6 |
2018/5/23 | 29779.3 |
2018/5/24 | 25971 |
discfile = ‘D:/PycharmProjects/OrderForecast/datafile/data.xls‘
forecastnum = 7
#讀取數據,指定日期列為指標,Pandas自動將“日期”列識別為Datetime格式
data = pd.read_excel(discfile, index_col=u‘日期‘)
#用來正常顯示中文標簽
plt.rcParams[‘font.sans-serif‘] = [‘SimHei‘]
#用來正常顯示負號
plt.rcParams[‘axes.unicode_minus‘] = False
data.plot()
訂單走勢如下圖:
3.2 時間序列的差分dd
? ARIMA 模型對時間序列的要求是平穩型。因此,當你得到一個非平穩的時間序列時,首先要做的即是做時間序列的差分,直到得到一個平穩時間序列。如果你對時間序列做d次差分才能得到一個平穩序列,那麽可以使用ARIMA(p,d,q)模型,其中d是差分次數。
#時間序列的差分d
dta=data[u‘訂單總額‘]
fig = plt.figure(figsize=(12,8))
ax1= fig.add_subplot(111)
diff1 = dta.diff(2)
diff1.plot(ax=ax1)
print(u‘原始序列的ADF檢驗結果為:‘, ADF(dta))
一階差分的時間序列的均值和方差已經基本平穩,不過我們還是可以比較一下二階差分的效果:
fig = plt.figure(figsize=(12,8))
ax2= fig.add_subplot(111)
diff2 = dta.diff(2)
diff2.plot(ax=ax2)
? 可以看出二階差分後的時間序列與一階差分相差不大,並且二者隨著時間推移,時間序列的均值和方差保持不變。因此可以將差分次數d設置為1。
3.3 合適的p,qp,q
得到一個平穩的時間序列後,接來下就是選擇合適的ARIMA模型,即ARIMA模型中合適的p,qp,q。
第一步要先檢查平穩時間序列的自相關圖和偏自相關圖。
dta= dta.diff(1)#已經知道要使用一階差分的時間序列,之前判斷差分的程序可以註釋掉
fig = plt.figure(figsize=(12,8))
ax1=fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(dta,lags=40,ax=ax2)
其中lags 表示滯後的階數,以上分別得到acf 圖和pacf 圖
通過兩圖觀察得到:
- 自相關圖顯示滯後有三個階超出了置信邊界;
- 偏相關圖顯示在滯後1至2階(lags 1,2...7)時的偏自相關系數超出了置信邊界,從lag 7之後偏自相關系數值縮小至0
則有以下模型可以供選擇:- ARMA(0,1)模型:即自相關圖在滯後1階之後縮小為0,且偏自相關縮小至0,則是一個階數q=1的移動平均模型;
- ARMA(7,0)模型:即偏自相關圖在滯後7階之後縮小為0,且自相關縮小至0,則是一個階層p=3的自回歸模型;
- ARMA(7,1)模型:即使得自相關和偏自相關都縮小至零。則是一個混合模型。
4.還可以有其他供選擇的模型
現在有以上這麽多可供選擇的模型,我們通常采用ARMA模型的AIC法則。
arma_mod20 = sm.tsa.ARMA(dta,(7,0)).fit()
print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)
print(sm.stats.durbin_watson(arma_mod20.resid.values))
arma_mod30 = sm.tsa.ARMA(dta,(0,1)).fit()
print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)
print(sm.stats.durbin_watson(arma_mod30.resid.values))
arma_mod40 = sm.tsa.ARMA(dta,(7,1)).fit()
print(arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic)
print(sm.stats.durbin_watson(arma_mod40.resid.values))
arma_mod50 = sm.tsa.ARMA(dta,(8,0)).fit()
print(arma_mod50.aic,arma_mod50.bic,arma_mod50.hqic)
print(sm.stats.durbin_watson(arma_mod50.resid.values))
3001.733127872179 3029.461447568363 3013.5940088028337 1.9672550548284564
3174.789227240986 3183.698667139714 3178.4095208845374 0.7514460553259659
3005.171853576973 3034.8699865727326 3017.239499055478 1.945755370815453
3001.8331059050256 3031.5312389007854 3013.900751383531 1.9401970059406308
可以看到ARMA(7,0)的aic,bic,hqic值最小,因此是最佳模型。
3.4 模型檢驗
在指數平滑模型下,觀察ARIMA模型的殘差是否是平均值為0且方差為常數的正態分布(服從零均值、方差不變的正態分布),同時也要觀察連續殘差是否(自)相關。
3.4.1 對ARMA(7,0)模型所產生的殘差做自相關圖
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
3.4.2 做D-W檢驗
德賓-沃森(Durbin-Watson)檢驗。德賓-沃森檢驗,簡稱D-W檢驗,是目前檢驗自相關性最常用的方法,但它只使用於檢驗一階自相關性。因為自相關系數ρ的值介於-1和1之間,所以 0≤DW≤4。
並且DW=0=>ρ=1 即存在正自相關性
DW=4<=>ρ=-1 即存在負自相關性
DW=2<=>ρ=0 即不存在(一階)自相關性
因此,當DW值顯著的接近於O或4時,則存在自相關性,而接近於2時,則不存在(一階)自相關性。這樣只要知道DW統計量的概率分布,在給定的顯著水平下,根據臨界值的位置就可以對原假設H0進行檢驗。
print(sm.stats.durbin_watson(arma_mod20.resid.values))
檢驗結果是1.9672550548284564,說明不存在自相關性。
3.4.3 觀察是否符合正態分布
這裏使用QQ圖,它用於直觀驗證一組數據是否來自某個分布,或者驗證某兩組數據是否來自同一(族)分布。在教學和軟件中常用的是檢驗數據是否來自於正態分布。
resid = arma_mod20.resid#殘差
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line=‘q‘, ax=ax, fit=True)
3.5 模型預測
模型確定之後,就可以開始進行預測了,對未來七天的數據進行預測。
predict_sunspots = arma_mod20.predict(‘2018-05-23‘, ‘2018-05-29, dynamic=True)
print(predict_sunspots)
fig=plt.figure(figsize=(12, 8))
ax3=fig.add_subplot(713)
predict_sunspots.plot(ax=ax3)
預測結果如下:
? 從圖形來,預測結果較為合理。
基於AIRMA模型對訂單總額未來七天的預測