利用ARIMA算法建立短期預測模型

　　周五福利日活動是電信為回饋老用戶而做的活動，其主要回饋老用戶的方式是讓用戶免費領取對應的優惠券，意在提升老用戶的忠誠度和活躍度。今日，為保證倉庫備貨優惠券資源充足，特別是5元話費券等，需要對該類優惠券領取效果進行預測，從而指導備貨。經研究選用ARIMA算法建立預測模型，對5元話費券進行日領取量的短期預測。數據集收集了2019年1月到2019年2月5元話費券的日領取量數據，並根據此數據做時間序列分析並建立預測模型。

　　1、進行數據的加載

from statsmodels.tsa.stattools import acf, pacf
import statsmodels.api as sm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as stats

receive=pd.read_excel(r'E:\Ering\data\HF_5.xlsx')

　　2、進行一階差分和檢驗一階差分的效果

#設置一下時間索引並進行一階差分
receive.index=pd.Index(pd.date_range('1/1/2019','31/3/2019',freq='1D'))
receive['number'].plot()
receive['number'].diff(1).plot() 

　　利用自相關系數的白噪聲檢驗差分效果

# 利用自相關系數的白噪聲檢驗差分后的數據是否是平穩序列：
r,q,p=sm.tsa.acf(receive['number'].diff(1).iloc[1:92].values.squeeze(),qstat=True) #squeeze: 除去size為1的維度
mat=np.c_[range(1,41),r[1:],q,p] #np.c_是按行連接兩個矩陣，把兩矩陣左右相加，要求行數相等，類似pandas的merge()
table=pd.DataFrame(mat,columns=['lag','AC','Q','Prob(>Q)'])
LB_result=table.iloc[[5,11,17]]
LB_result.set_index('lag',inplace=True)
LB_result

　　上圖是自相關系數白噪聲檢驗的結果，P值均大於0.05，表明白噪聲檢驗不顯著，所以數據經過一階差分后平穩，故此知道一階差分后所形成的序列是平穩序列。ARIMA算法的一個重要應用前提是保證算法入口的數據是平穩的。差分階次對應ARIMA（p,d,q）中的d數值，因此本文中的d可以設置為1。

　　3、確定系數p和q

# 利用BIC最小的模型作為識別的依據,確定參數p和參數q：
order_p,order_q,bic=[],[],[]
model_order=pd.DataFrame()
for p in range(4):
    for q in range(4):
        arma_model=sm.tsa.ARMA(receive['number'].diff(1).iloc[1:89].dropna(),
                               (p,q)).fit()
        order_p.append(p)
        order_q.append(q)
        bic.append(arma_model.bic)
        print('The BIC of ARMA(%s,%s) is %s'%(p,q,arma_model.bic))

model_order['p']=order_p
model_order['q']=order_q
model_order['BIC']=bic
P=list(model_order['p'][model_order['BIC']==model_order['BIC'].min()])
Q=list(model_order['q'][model_order['BIC']==model_order['BIC'].min()])
print('最好的模型是ARMA(%s,%s)' %(P[0],Q[0]))

　　

　　根據BIC法則確定可知，當p值和q值分別為1和0的時候，可以取到模型的效果最好。

　　4、建立ARIMA預測模型

#建立ARIMA模型
model=sm.tsa.ARMA(receive['number'].diff(1).iloc[1:89].dropna(),(1,0)).fit(method='css') #使用最小二乘，‘mle’是極大似然估計
#畫圖比較一下預測值和真實觀測值之間的關系
fig=plt.figure(figsize=(8,6))
ax=fig.add_subplot(111)
ax.plot(receive['number'].diff(1).iloc[1:89],color='blue',label='number')
ax.plot(model.fittedvalues,color='red',label='Predicted number')
plt.legend(loc='lower right')

　　

　　5、差分值轉化為原始值

# 差分數據轉化為原始值
def forecast(step,var,modelname):
    diff=list(modelname.predict(len(var)-1,len(var)-1+step,dynamic=True))
    prediction=[]
    prediction.append(var[len(var)-1])
    seq=[]
    seq.append(var[len(var)-1])
    seq.extend(diff)
    for i in range(step):
        v=prediction[i]+seq[i+1]
        prediction.append(v)

    prediction=pd.DataFrame({'Predicted number':prediction})
    return prediction[1:]  #第一個值是原序列最后一個值，故第二個值是預測值。
forecast(15,receive['number'][1:89],model)

　　如圖是4月1-5號的預測值。

　　6、模型殘差項的白噪聲檢驗及正態性檢驗

# 模型殘差項的白噪聲檢驗：
resid=model.resid
r,q,p=sm.tsa.acf(resid.values.squeeze(),qstat=True)
mat_res=np.c_[range(1,41),r[1:],q,p] #np.c_是按行連接兩個矩陣，把兩矩陣左右相加，要求行數相等，類似pandas的merge()
table=pd.DataFrame(mat_res,columns=['to lag','AC','Q','Prob(>Q)'])
LB_result_res=table.iloc[[5,11,17,23]]
LB_result_res.set_index('to lag',inplace=True)
LB_result_res

　　

　　如果ARMA模型估計的好，應當使得估計值后的殘差項是白噪聲。上圖是預測結果的殘差的白噪聲檢驗結果，分析可知Prob值均較大，查閱資料顯示Prob值較大時，接受原假設-殘差是白噪聲Prob值接近於0時拒絕原假設；接近於1時接受原假設；Prob值為10%時，表示10%置信區間下通過。

　　至此模型建立完畢。