python的ARIMA Python的arima模型

近来，一个项目需要预测数据，想到了利用ARIMA算法来解决这个问题，遂将最近一段时间关于ARIMA算法的研究内容做以总结。

【ARIMA算法介绍】

【ARMA与ARIMA】

“ARIMA”实际上并不是一整个单词，而是一个缩写。其全称是：Autoregressive Integrated Moving Average Model，即自回归移动平均模型。它属于统计模型中最常见的一种，用于进行时间序列的预测。其原理在于：在将非平稳时间序列转化为平稳时间序列的过程中，将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。 ARIMA模型与ARMA模型的区别：ARMA模型是针对平稳时间序列建立的模型。ARIMA模型是针对非平稳时间序列建模，需要对时间序列先做差分，平稳后再建立ARMA模型。换句话说，非平稳时间序列要建立ARMA模型，首先需要经过差分转化为平稳时间序列，然后建立ARMA模型。

【Python实现】
#先说明一下哈，本文不选用项目数据，示例数据源是一个通信运营商的销售数据，代码都相同，换一个数据源即可。
#示例问题描述：某通信运营商为了研究其某款产品的销售情况，收集了从2001年1月到2008年8月的销售额月度数据。请根据此数据作时间序列分析和预测。
###################先确定差分阶次d#####################

#加载数据
import pandas as pd
sales=pd.read_sas('sales_monthly.sas7bdat')
sales.head(12)
#设置一下时间索引，这一步很有必要，否则后面程序会报错
sales.index=pd.Index(pd.date_range('1/2001','9/2008',freq='1M'))
sales['Sales'].plot()
sales['Sales'].diff(1).plot() #可以看出1阶差分后是一个零均值化的平稳序列,第一个数是NaN

从上图可以看出，原数据序列是不平稳的，再一阶差分后数据平稳了，ARIMA算法的一个重要应用前提是保证算法入口的数据是平稳的。差分阶次对应ARIMA（p,d,q）中的d数值，因此本文中的d可以设置为1。

###################确定p和q#####################

#ARIMA模型分析：
#先把ACF图和PACF图画出来看看：
fig=plt.figure(figsize=(10,6))
ax1=fig.add_subplot(211)
fig=sm.graphics.tsa.plot_acf(sales['Sales'].diff(1).iloc[1:92].dropna(),lags=24,ax=ax1) # 注意：要去掉第1个空值                             
ax2=fig.add_subplot(212)
fig=sm.graphics.tsa.plot_pacf(sales['Sales'].diff(1).iloc[1:92].dropna(),  lags=24,ax=ax2)# 注意：要去掉第1个空值
#判断：ACF图在2之后截尾，而PACF拖尾。模型可以由MA(2)→ARIMA(0,1,2).

关于截尾与拖尾：

拖尾：始终有非零取值，不会在k大于某个常数后就恒等于零(或在0附近随机波动)
截尾：在大于某个常数k后快速趋于0为k阶截尾

AR模型：自相关系数拖尾，偏自相关系数截尾；
MA模型：自相关系数截尾，偏自相关函数拖尾；
ARMA模型：自相关函数和偏自相关函数均拖尾。

AR模型对应p，根据偏自相关图的截尾位置确定p，ACF与PACF图的横轴均从0开始，如果PACF在2时还在置信区间以外，从3开始之后均在置信区间以内，则定p为2。MA模型对应q，根据自相关图的截尾位置确定q，确定方法同p。

###################模型预测#####################

from statsmodels.tsa.stattools import acf, pacf
import statsmodels.api as sm
import matplotlib.pyplot as plt
import numpy as np
#建立ARIMA模型
model=sm.tsa.ARMA(sales['Sales'].diff(1).iloc[1:92].dropna(),(0,2)).fit(method='css') #使用最小二乘，‘mle’是极大似然估计
#画图比较一下预测值和真实观测值之间的关系
fig=plt.figure(figsize=(8,6))
ax=fig.add_subplot(111)
ax.plot(sales['Sales'].diff(1).iloc[1:92],color='blue',label='Sales')
ax.plot(model.fittedvalues,color='red',label='Predicted Sales')
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 Sales':prediction})
    return prediction[1:]  #第一个值是原序列最后一个值，故第二个值是预测值。
forecast(4,sales['Sales'],model)

至此，完成了ARIMA算法的主体部分。如果想锦上添花，把工作做得细致一些，在(p,d,q)参数的选择以及模型效果的检验等方面还有其它工作可做。

首先，是自相关系数的白噪声检验，这个是用来判断差分后，数据是否平稳。

# 还可以再进行自相关系数的白噪声检验：
r,q,p=sm.tsa.acf(sales['Sales'].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，表明白噪声检验不显著，所以数据经过一阶差分后平稳。

上面这个截图是示例数据结果，P值均大于0.05，表明白噪声检验不显著，所以数据经过一阶差分后平稳。但是，项目这一步结果的P均小于0.05，但是也不影响后续的预测效果，先忽略吧。

其次，是单位根检验，这一步的目的同样是判断差分后的数据是否平稳。上代码！

# 单位根检验(原假设H0：时间序列具有单位根（不平稳）)：
# http://www.statsmodels.org/stable/generated/statsmodels.tsa.stattools.adfuller.html?
# highlight=adfuller#statsmodels.tsa.stattools.adfuller
from statsmodels.tsa.stattools import adfuller
def DFTest(sales,regression, maxlag,autolag='AIC'):
    print("ADF-Test Result:")
    dftest=adfuller(sales,regression=regression,maxlag=maxlag, autolag=autolag)

    dfoutput=pd.Series(dftest[0:4],
                   index=['Test Statistic','P-value',
                          'Lags Used','nobs'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value at %s'%key]=value
    print(dfoutput)

# regression可以根据模型形式指定为：
# ‘c'：仅常数，默认；’ct‘：常数和长期趋势；’ctt’：常数、线性和二次曲线趋势；‘nc’：无常数趋势。

DFTest(sales['Sales'],regression='nc', maxlag=6,autolag='AIC')  #对原始数据进行单位根检验

print(37*'-')
DFTest(sales['Sales'].diff(1).iloc[1:92], regression='nc', maxlag=5,autolag='AIC')
# 结论：一阶差分后平稳。

------上面一个是原始数据，差分前的单位根检验结果，下面一个是差分后的检验结果。单位根检验结果的分析主要看两个方面：第一，ADF-Test result同时小于1%、5%、10%即说明非常好地拒绝该假设(原假设H0：时间序列具有单位根（不平稳）)。第二，P-value是否非常接近0，本示例中，P-value 为 9e-10,接近0。

综上，差分后，本示例数据平稳了。

然后，在确定(p,d,q)参数方面，还有一种方法—BIC最小原则。

# 利用BIC最小的模型作为识别的依据：
order_p,order_q,bic=[],[],[]
model_order=pd.DataFrame()
for p in range(4):
    for q in range(4):
        arma_model=sm.tsa.ARMA(sales['Sales'].diff(1).iloc[1:92].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('\n最好的模型是ARMA(%s,%s)' %(P[0],Q[0]))

下面是运行结果：

分析程序发现，将p与q均设置一个遍历范围[0,3]，步长为1，最终返回BIC值最小对应的p与q。

继而，估计并检验模型参数。

# 模型参数估计及检验：
   model=sm.tsa.ARMA(sales['Sales'].diff(1).iloc[1:92].dropna(),
                                   (0,2)).fit(method='css') #使用最小二乘，‘mle’是极大似然估计
   params=model.params
   tvalues=model.tvalues
   pvalues=model.pvalues
   result_mat=pd.DataFrame({'Estimate':params,'t-values':tvalues,'pvalues':pvalues})
   result_mat  # 注意：这里得到的是-θ的结果（Xt=at-θ1a(t-1)-θ2a(t-2)）。从t值来看是显著的。
    
   model.aic
   print('方差估计值是：%f' %model.sigma2)

进一步，

# 如果ARMA模型估计的好，应当使得估计值后的残差项是白噪声。下面对模型的残差项进行白噪声检验：
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
# 原假设为残差是白噪声。

# 再看一下正态性：
import scipy.stats as stats
sm.ProbPlot(resid, stats.t,fit=True).ppplot(line='45')
sm.ProbPlot(resid, stats.t,fit=True).qqplot(line='45')
plt.show()

plt.figure()
x=pd.Series(resid)
p1=x.plot(kind='kde')   #可以是’line’, ‘bar’, ‘barh’, ‘kde’.kde是密度图。
p2=x.hist(normed=True)
plt.grid(True)
plt.show()

下图是预测结果的残差的白噪声检验结果，分析可知Prob值均较大，查阅资料显示Prob值较大时，接受原假设-残差是白噪声。

Prob值接近于0时拒绝原假设；接近于1时接受原假设；Prob值为10%时，表示10%置信区间下通过。

下图是预测结果残差的pp图与qq图，这两个图通常用来分析数据是否符合正态分布，这两个图中，数据点均沿着红色斜线分布，因此认为预测结果的残差满足正态分布。关于pp图与qq图，可以看一下参考文献中的两篇文章就会明白，不做赘述。

下面这个是pp图。

这个是qq图。

下图是预测结果残差的概率密度图，呈现出正态分布的特性。

至此，完成了ARIMA算法的建模，平稳性、白噪声检验、残差正态性分析、BIC分析等工作。简单的总结一下，ARIMA算法建模需要先看一下原数据是否平稳，如果平稳直接建立ARMA模型，如果不平稳则可能需要建立ARIMA模型，两者的区别就是是否需要对原数据进行差分处理，差分的目的是让数据变得平稳。之后，可以通过ACF、PACF图或者BIC最小原则确定p与q。最后，根据(p,d,q)建模，进而预测。文末提及的各种检验工作只是对中间结果的辅助分析，用以评判当前的参数是否适用。

项目开发过程中，发现一个有趣的现象，如果要预测60步结果，则前20步结果还是有一定的波动的，21步以后，预测结果便呈现明显的线性变化关系，很直的一条线哦，初步来看这样的预测效果不会太好！这也让我想起之前文献中的一句话，“ARIMA可以做短期预测”，这其中的深层原因还待后续研究！