from statsmodels.tsa.stattools import acf, pacf
import statsmodels.api as sm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pdimport scipy.stats as statsreceive=pd.read_excel(r'E:\Ering\data\HF_5.xlsx')

#设置一下时间索引并进行一阶差分
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

# 利用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]))

#建立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')

# 差分数据转化为原始值
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)

# 模型残差项的白噪声检验：
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