文章来源于网络收集而来，版权归原创者所有，如有侵权请及时联系！

11.1 ARIMA 模型

发布于 2024-01-21 22:13:24 字数 11512 浏览 0 评论 0 收藏 0

下面应用Python语言建模步骤，对表11-2中2013年1月到2016年1月某餐厅的营业数据进行建模。

表11-2 某餐厅的销量数据

*数据详见：示例程序/data/arima_data.csv

1.时间序列对象

加载基础库：pandas，numpy，scipy，matplotlib，statsmodels对其调用如下：

import pandas as pd

2.获取数据

从xls文件中读取数据，代码见代码清单11-1：

代码清单11-1 获取数据

# 参数初始化
discfile = '../data/arima_data.xls'
# 读取数据，指定日期列为指标，
Pandas自动将
"日期
"列识别为
Datetime格式
data = pd.read_excel(discfile,index_col=0)
print(data.head())
print('\n Data Types:')
print(data.dtypes)

3.绘制时间序列图

对读取的数据绘制时间序列图，观察图形的特征，代码见代码清单11-2，所得图形见图11-1。

代码清单11-2 绘制时间序列图

# 时序图
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']     #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False     #用来正常显示负号
data.plot()
plt.show()

4.自相关

对时间序列做自相关图（见图11-2），判断序列是否自相关，代码见代码清单11-3：

代码清单11-3 自相关图

#自相关图
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data).show()

图11-1 时间序列图

图11-2 自相关图

由自相关图可以看出，在4阶后才落入区间内，并且自相关系数长期大于零，显示出很强的自相关性。

5.平稳性检验

还需要对时间序列做平稳性检验，代码见代码清单11-4：

代码清单11-4 平稳性检验

#平稳性检测
from statsmodels.tsa.stattools import adfuller as ADF
print(u'原始序列的
ADF检验结果为：
', ADF(data[u'销量
']))
#返回值依次为
adf、
pvalue、
usedlag、
nobs、
critical values、
icbest、
regresults、
resstore
#result:原始序列的
ADF检验结果为：
 (1.8137710150945272, 0.99837594215142644, 10L, 
26L, {'5%': -2.9812468047337282, '1%': -3.7112123008648155, '10%': -2.6300945562130176}, 
299.46989866024177)

从返回的结果可以看出检验结果的pvalue即p值显著大于0.05，判断该序列为非平稳序列。

6.时间序列的差分d

ARIMA模型对时间序列的要求是平稳型。因此，当你得到一个非平稳的时间序列时，首先要做的是时间序列的差分，直到得到一个平稳时间序列。如果你对时间序列做d次差分才能得到一个平稳序列，那么可以使用ARIMA（p，d，q）模型，其中d是差分次数。

首先对时间序列做差分，并观察差分后的时序图，见图11-3，代码如代码清单11-5：

代码清单11-5 做差分并绘制时序图

D_data = data.diff().dropna()
D_data.columns = [u'销量差分
']
D_data.plot() #时序图
plt.show()

图11-3 差分后的时序图

对差分后的序列做自相关检验，见图11-4，观察是否自相关，代码如代码清单11-6：

代码清单11-6 对序列做自相关检验

plot_acf(D_data).show() #自相关图

由图11-4可以看出，差分后的序列迅速落入区间内，并呈现出向0靠拢的趋势，序列没有自相关性。

对差分后的序列做偏自相关检验，见图11-5，观察是否偏自相关，代码如代码清单11-7：

代码清单11-7 对序列做偏自相关检验

from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(D_data).show() #偏自相关图

图11-4 差分后的自相关图

图11-5 差分后的偏自相关图

由偏自相关图可以看出，差分后的序列也没有显示出偏自相关性。

再对差分后的序列做平稳性检测，代码见代码清单11-8：

代码清单11-8 平稳性检测

#平稳性检测
print(u'差分序列的
ADF检验结果为：
', ADF(D_data[u'销量差分
'])) 
#result:差分序列的
ADF检验结果为：
 (-3.1560562366723537, 0.022673435440048798, 0L, 
35L, {'5%': -2.9485102040816327, '1%': -3.6327426647230316, '10%': -2.6130173469387756}, 287.59090907803341)

从返回的ADF检验结果得到，p值为0.022673435440048798，小于0.05。

还需要对差分后的序列做白噪声检验，见代码清单11-9：

代码清单11-9 白噪声检验

#白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
print(u'差分序列的白噪声检验结果为：
', acorr_ljungbox(D_data, lags=1))
#返回统计量和
p值
#result:差分序列的白噪声检验结果为：
 (array([ 11.30402222]), array([ 0.00077339]))

从得到的白噪声检验结果可以看出，检验的p值为0.00077339，小于0.05，通过白噪声检验，序列为白噪声序列。

接下来我们比较下一阶差分后的序列和二阶差分后的序列，一阶差分序列的代码见代码清单11-10：

代码清单11-10 一阶差分序列

# 一阶差分
fig = plt.figure(figsize=(12,8))
ax1= fig.add_subplot(111)
diff1 = data.diff(1)
diff1.plot(ax=ax1)

图11-6 一阶差分后的序列图

一阶差分的时间序列的均值和方差已经基本平稳，不过我们还是可以比较一下二阶差分的效果，二阶差分的代码见代码清单11-11：

代码清单11-11 二阶差分序列

# 二阶差分
fig = plt.figure(figsize=(12,8))
ax2= fig.add_subplot(111)
diff2 = dta.diff(2)
diff2.plot(ax=ax2)

图11-7 二阶差分后的时序图

可以看出二阶差分后的时间序列与一阶差分相差不大，并且二者随着时间推移，时间序列的均值和方差保持不变。因此可以将差分次数d设置为1。

7.合适的p，q

现在我们已经得到一个平稳的时间序列，接下来就是选择合适的ARIMA模型，即ARIMA模型中合适的p，q。

第一步我们要先检查平稳时间序列的自相关图和偏自相关图，代码见代码清单11-12。

代码清单11-12 模型定阶

# 合适的
p,q
dta = data.diff(1)[1:]
fig = plt.figure(figsize=(12,8))
ax1=fig.add_subplot(211)
fig1 = sm.graphics.tsa.plot_acf(dta[u'销量
'],lags=10,ax=ax1)
ax2 = fig.add_subplot(212)
fig2 = sm.graphics.tsa.plot_pacf(dta[u'销量
'],lags=10,ax=ax2)

其中lags表示滞后的阶数，以上分别得到acf图和pacf图，见图11-8：

图11-8 自相关和偏自相关图

通过两图观察得到：

1）自相关图显示滞后有两个阶超出了置信边界。

2）偏相关图显示在滞后1阶时的偏自相关系数超出了置信边界，从lag 1之后偏自相关系数值缩小至0。

则有以下模型可以供选择：

1）ARMA（0，2）模型：即自相关图在滞后2阶之后缩小为0，且偏自相关缩小至0，则是一个阶数q=2的移动平均模型。

2）ARMA（1，0）模型：即偏自相关图在滞后1阶之后缩小为0，且自相关缩小至0，则是一个阶数p=1的自回归模型。

3）ARMA（0，1）模型：即自相关图在滞后1阶之后缩小为0，且偏自相关缩小至0，则是一个阶数q=1的自回归模型。

现在有以上这么多可供选择的模型，我们通常采用ARMA模型的AIC法则。我们知道：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性，但是应尽量避免出现过度拟合（Overfitting）的情况。所以优先考虑的模型应是AIC值最小的那一个。赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。不仅仅包括AIC准则，目前选择模型常用如下准则：

1）AIC=-2 ln（L）+2 k中文名字：赤池信息量Akaike Information Criterion

2）BIC=-2 ln（L）+ln（n）*k中文名字：贝叶斯信息量Bayesian Information Criterion

3）HQ=-2 ln（L）+ln（ln（n））*k Hannan-Quinn Criterion

构造这些统计量所遵循的统计思想是一致的，就是在考虑拟合残差的同时，依自变量个数施加“惩罚”。但要注意的是，这些准则不能说明某一个模型的精确度，也就是说，对于三个模型A、B、C，我们能够判断出C模型是最好的，但不能保证C模型能够很好地刻画数据，因为有可能三个模型都是糟糕的。

对三个模型分别做AIC、BIC、HQ统计量检验，代码见代码清单11-13。

代码清单11-13 AIC、BIC、HQ、统计量检验

#模型
arma_mod20 = sm.tsa.ARMA(dta,(2,0)).fit()
print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)
arma_mod01 = sm.tsa.ARMA(dta,(0,1)).fit()
print(arma_mod01.aic,arma_mod01.bic,arma_mod01.hqic)
arma_mod10 = sm.tsa.ARMA(dta,(1,0)).fit()
print(arma_mod10.aic,arma_mod10.bic,arma_mod10.hqic)
#result:
print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)
420.440748036 426.77482379 422.651510127
print(arma_mod01.aic,arma_mod01.bic,arma_mod01.hqic)
417.759525387 422.510082203 419.417596956
print(arma_mod10.aic,arma_mod10.bic,arma_mod10.hqic)
418.877719334 423.628276149 420.535790902

对比3个模型，可以看到ARMA（0，1）的aic、bic、hqic均最小，因此是最佳模型。

8.模型检验

对于选择的模型，观察ARIMA模型的残差是否是平均值为0且方差为常数（服从零均值、方差不变的正态分布），同时也要观察连续残差是否（自）相关。

首先使用QQ图，它用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布，代码见代码清单11-14，其效果图如图11-9所示。

代码清单11-14 残差QQ图

#残差
QQ图
resid = arma_mod01.resid
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)

图11-9 残差QQ图

我们对ARMA（0，1）模型所产生的残差做自相关图，代码见代码清单11-15，其效果图如图11-10所示。

代码清单11-15 残差自相关检验

#残差自相关检验
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(arma_mod01.resid.values.squeeze(), lags=10, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(arma_mod01.resid, lags=10, ax=ax2)

图11-10 残差自相关检验效果图

还需要对残差做D-W检验。

德宾–沃森（Durbin-Watson）检验，简称D-W检验，是目前检验自相关性最常用的方法，但它只适用于检验一阶自相关性。因为自相关系数ρ的值介于-1和1之间，所以0≤DW≤4。并且

DW＝0＝＞ρ＝1 即存在正自相关性

DW＝4＜＝＞ρ＝-1 即存在负自相关性

DW＝2＜＝＞ρ＝0 即不存在（一阶）自相关性

因此，当DW值显著的接近于0或4时，则存在自相关性，而接近于2时，则不存在（一阶）自相关性。这样只要知道DW统计量的概率分布，在给定的显著水平下，根据临界值的位置就可以对原假设进行检验。代码见代码清单11-16。

代码清单11-16 D-W检验

#D-W检验
print(sm.stats.durbin_watson(arma_mod01.resid.values))
#result:1.95414900233

检验结果是1.95414900233，说明不存在自相关性。

最后还需要对残差做Ljung-Box检验。

Ljung-Box检验是对随机性的检验，或者说是对时间序列是否存在滞后相关的一种统计检验。对于滞后相关的检验，我们常常采用的方法还包括计算ACF和PCAF并观察其图像，但是无论是ACF还是PACF都仅仅考虑是否存在某一特定滞后阶数的相关。LB检验则是基于一系列滞后阶数，判断序列总体的相关性或者随机性是否存在。

时间序列中一个最基本的模型就是高斯白噪声序列。而对于ARIMA模型，其残差被假定为高斯白噪声序列，所以当我们用ARIMA模型去拟合数据时，拟合后我们要对残差的估计序列进行LB检验，判断其是否是高斯白噪声，如果不是，那么就说明ARIMA模型也许并不是一个适合样本的模型。

对残差做Ljung-Box检验的代码见代码清单11-17。

代码清单11-17 Ljung-Box检验

# Ljung-Box检验
import numpy as np
r,q,p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
datap = np.c_[range(1,36), r[1:], q, p]
table = pd.DataFrame(datap, columns=['lag', "AC", "Q", "Prob(>Q)"])
print(table.set_index('lag'))
#result:
            AC          Q  Prob(>Q)
lag                                
1.0   0.009994   0.003904  0.950179
2.0   0.151097   0.922489  0.630498
3.0   0.119392   1.513404  0.679180
4.0  -0.212564   3.445002  0.486290
5.0   0.034075   3.496239  0.623957
6.0  -0.053349   3.626021  0.727135
7.0  -0.157088   4.790085  0.685562
8.0   0.082868   5.125590  0.744072
9.0   0.180436   6.775153  0.660516
10.0 -0.119683   7.528822  0.674754
11.0  0.051306   7.672864  0.742274
12.0 -0.062678   7.896792  0.793143
13.0 -0.020659   7.922177  0.848633
14.0 -0.078650   8.306822  0.872737
15.0 -0.024755   8.346742  0.909130
16.0  0.001821   8.346969  0.937859
17.0  0.081164   8.821276  0.945702
18.0  0.181184  11.316173  0.880461
19.0 -0.036607  11.424010  0.908751
20.0  0.049095  11.630091  0.928220
21.0  0.095998  12.470556  0.926035
22.0 -0.186408  15.865930  0.822484
23.0 -0.066136  16.326203  0.840953
24.0 -0.160985  19.280651  0.736861
25.0 -0.218461  25.215923  0.450326
26.0  0.054818  25.627015  0.483747
27.0 -0.067221  26.313862  0.501239
28.0 -0.073881  27.247238  0.504816
29.0  0.025513  27.374445  0.551511
30.0 -0.068816  28.454178  0.546382
31.0 -0.001377  28.454697  0.597610
32.0  0.008883  28.481683  0.645347
33.0 -0.015105  28.585729  0.686711
34.0 -0.001370  28.587014  0.730003
35.0 -0.001850  28.591694  0.769544

检验的结果就是看最后一列前十二行的检验概率（一般观察滞后1~12阶），如果检验概率小于给定的显著性水平，比如0.05、0.10等就拒绝原假设，其原假设是相关系数为零。就结果来看，如果取显著性水平为0.05，那么相关系数与零没有显著差异，即为白噪声序列。

9.模型预测

模型确定之后，就可以开始进行预测了，我们对未来9日的数据进行预测，代码见代码清单11-18。

代码清单11-18 模型预测

#预测
predict_sunspots = arma_mod01.predict('2015-2-07', '2015-2-15', dynamic=True)
fig, ax = plt.subplots(figsize=(12, 8))
print(predict_sunspots)
predict_sunspots[0] += data['2015-02-06':][u'销量
']
data=pd.DataFrame(data)
for i in range(len(predict_sunspots)-1):
    predict_sunspots[i+1]=predict_sunspots[i]+predict_sunspots[i+1]
print(predict_sunspots)
ax = data.ix['2015':].plot(ax=ax)
predict_sunspots.plot(ax=ax)
plt.show()

*代码详见：示例程序/code/11-1.py

所得的时序图如图11-11所示。