Question

差分运算具有强大的确定性信息提取能力，许多非平稳序列差分后会显示出平稳序列的性质，这时称这个非平稳序列为差分平稳序列。对差分平稳序列可以使用ARMA模型进行拟合。ARIMA模型的实质就是差分运算与ARMA模型的组合，掌握了ARMA模型的建模方法和步骤以后，对序列建立ARIMA模型是比较简单的。

差分平稳时间序列建模步骤如图5-18所示。

图5-18　差分平稳时间序列建模步骤

下面应用以上的理论知识，对表5-25中2015/1/1~2015/2/6某餐厅的销售数据进行建模。

表5-25　某餐厅的销量数据

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

（1）检验序列的平稳性

图5-19　原始序列的时序图

图5-20　原始序列的自相关图

表5-26　原始序列的单位根检验

图5-19时序图显示该序列具有明显的单调递增趋势，可以判断为是非平稳序列；图5-20的自相关图显示自相关系数长期大于零，说明序列间具有很强的长期相关性；表5-26单位根检验统计量对应的p值显著大于0.05，最终将该序列判断为非平稳序列（非平稳序列一定不是白噪声序列）。

（2）对原始序列进行一阶差分，并进行平稳性和白噪声检验

1）对一阶差分后的序列再次做平稳性判断。

过程同上。

图5-21　一阶差分之后序列的时序图

图5-22　一阶差分之后序列的自相关图

一阶差分之后序列的单位根检验如表5-27所示。

表5-27　一阶差分之后序列的单位根检验

结果显示，一阶差分之后的序列的时序图在均值附近比较平稳的波动、自相关图有很强的短期相关性、单位根检验p值小于0.05，所以一阶差分之后的序列是平稳序列。

2）对一阶差分后的序列做白噪声检验（结果见表5-28）。

表5-28　一阶差分之后序列的白噪声检验

输出的p值远小于0.05，所以一阶差分之后的序列是平稳非白噪声序列。

（3）对一阶差分之后的平稳非白噪声序列拟合ARMA模型

下面进行模型定阶。模型定阶就是确定p和q。

第一种方法：人为识别的方法，根据表5-24进行模型定阶。

一阶差分后自相关图（见图5-23）显示出1阶截尾，偏自相关图显示出拖尾性，所以可以考虑用MA（1）模型拟合1阶差分后的序列，即对原始序列建立ARIMA（0，1，1）模型。

图5-23　一阶差分后序列的偏自相关图

第二种方法：相对最优模型识别。

计算ARMA（p，q）。当p和q均小于等于3的所有组合的BIC信息量，取其中BIC信息量达到最小的模型阶数。

计算完成BIC矩阵如下。

p值为0、q值为1时最小BIC值为：430.1374。p、q定阶完成！

用AR（1）模型拟合一阶差分后的序列，即对原始序列建立ARIMA（0，1，1）模型。虽然两种方法建立的模型是一样的，但模型是非唯一的，可以检验ARIMA（1，1，0）和ARIMA（1，1，1），这两个模型也能通过检验。

下面对一阶差分后的序列拟合AR（1）模型进行分析。

1）模型检验。残差为白噪声序列，p值为：0.627016。

2）参数检验和参数估计见表5-29。

表5-29　模型参数

（4）ARIMA模型预测

应用ARIMA（0，1，1）对表530中2015/1/1~2015/2/6某餐厅的销售数据做为期5天的预测，结果如下。

需要说明的是，利用模型向前预测的时期越长，预测误差将会越大，这是时间预测的典型特点。

在Python中实现ARIMA模型建模过程的代码如代码清单5-7所示。可以看到，我们使用了StatsModels，读者或许记得，我们在第2章介绍了它，在这里才真正用上它。这表明对于通常的数据探索任务来说，Numpy与Pandas的结合已经相当强大了，只有到较为深入的统计模型之时，才用到StatsModels。

代码清单5-7　ARIMA模型实现代码

#-*- coding: utf-8 -*-
#arima时序模型
import pandas as pd
#参数初始化
discfile = '../data/arima_data.xls'
forecastnum = 5
#读取数据，指定日期列为指标，Pandas自动将“日期”列识别为Datetime格式
data = pd.read_excel(discfile, index_col = u'日期')
#时序图
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号
data.plot()
plt.show()
#自相关图
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data).show()
#平稳性检测
from statsmodels.tsa.stattools import adfuller as ADF
print(u'原始序列的ADF检验结果为：', ADF(data[u'销量']))
#返回值依次为adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore
#差分后的结果
D_data = data.diff().dropna()
D_data.columns = [u'销量差分']
D_data.plot() #时序图
plt.show()
plot_acf(D_data).show() #自相关图
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(D_data).show() #偏自相关图
print(u'差分序列的ADF检验结果为：', ADF(D_data[u'销量差分'])) #平稳性检测
#白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
print(u'差分序列的白噪声检验结果为：', acorr_ljungbox(D_data, lags=1)) #返回统计量和p值
from statsmodels.tsa.arima_model import ARIMA
#定阶
pmax = int(len(D_data)/10) #一般阶数不超过length/10
qmax = int(len(D_data)/10) #一般阶数不超过length/10
bic_matrix = [] #bic矩阵
for p in range(pmax+1):
  tmp = []
  for q in range(qmax+1):
    try: #存在部分报错，所以用try来跳过报错。
      tmp.append(ARIMA(data, (p,1,q)).fit().bic)
    except:
      tmp.append(None)
  bic_matrix.append(tmp)
bic_matrix = pd.DataFrame(bic_matrix) #从中可以找出最小值
p,q = bic_matrix.stack().idxmin() #先用stack展平，然后用idxmin找出最小值位置。
print(u'BIC最小的p值和q值为：%s、%s' %(p,q))
model = ARIMA(data, (p,1,q)).fit() #建立ARIMA(0, 1, 1)模型
model.summary2() #给出一份模型报告
model.forecast(5) #作为期5天的预测，返回预测结果、标准误差、置信区间。

代码详见：示例程序/code/arima_test.py

运行代码清单5-7可以得到输出结果如下。

原始序列的ADF检验结果为： (1.8137710150945285, 0.99837594215142644, 10, 26, {'5%': 
   -2.9812468047
337282, '10%': -2.6300945562130176, '1%': -3.7112123008648155}, 299.46989866024177)差分序列的ADF检验结果为： (-3.1560562366723537, 0.022673435440048798, 0, 35, {'5%': 
   -2.948510204
0816327, '10%': -2.6130173469387756, '1%': -3.6327426647230316}, 287.59090907803341)差分序列的白噪声检验结果为： (array([ 11.30402222]), array([ 0.00077339]))
BIC最小的p值和q值为：0、1
                           Results: ARIMA
====================================================================
Model:              ARIMA            Log-Likelihood:      -205.88
Dependent Variable: D.销量             Scale:               1.0000
Date:               2015-08-06 10:24 Method:              css-mle
No. Observations:   36               Sample:              01-02-2015
Df Model:           2                                     02-06-2015
Df Residuals:       34               S.D. of innovations: 73.086
AIC:                417.7595         HQIC:                419.418
BIC:                422.5101
----------------------------------------------------------------------
               Coef.    Std.Err.     t      P>|t|     [0.025    0.975]
----------------------------------------------------------------------
const         49.9560    20.1390   2.4806   0.0182   10.4843   89.4277
ma.L1.D.销量     0.6710     0.1648   4.0712   0.0003    0.3480    0.9941
-----------------------------------------------------------------------------
                 Real           Imaginary          Modulus          Frequency
-----------------------------------------------------------------------------
MA.1           -1.4902             0.0000           1.4902             0.5000
(array([ 4873.96649322,  4923.9224731 ,  4973.87845298,  5023.83443286,
        5073.79041274]), array([  73.08574266,  142.32680564,  187.54283091,  
            223.80283119,
        254.95705731]), array([[ 4730.72106983,  5017.21191661],
       [ 4644.96706002,  5202.87788618],
       [ 4606.30125884,  5341.45564712],
       [ 4585.18894409,  5462.47992162],
       [ 4574.08376281,  5573.49706266]]))

5.4.5　Python主要时序模式算法

Python实现时序模式的主要库是StatsModels（当然，如果Pandas能做的，就可以利用Pandas先做），算法主要是ARIMA模型，在使用该模型进行建模时，需要进行一系列判别操作，主要包含平稳性检验、白噪声检验、是否差分、AIC和BIC指标值、模型定阶，最后再做预测。与其相关的函数见表5-30。

表5-30　时序模式算法函数列表

（1）acf（）

功能：计算自相关系数

使用格式：

autocorr=acf（data，unbiased=False，nlags=40，qstat=False，fft=False，alpha=None）

输入参数data为观测值序列（即为时间序列，可以是DataFrame或Series），返回参数autocorr为观测值序列自相关函数。其余为可选参数，如qstat=True时同时返回Q统计量和对应p值。

（2）plot_acf（）

功能：画自相关系数图

使用格式：

p=plot_acf（data）

返回一个Matplotlib对象，可以用.show（）方法显示图像。

（3）acf（）/plot_acf（）

功能：计算偏相关系数/画偏相关系数图

使用格式：使用跟acf（）/plot_acf（）类似，不再赘述。

（4）adfuller（）

功能：对观测值序列进行单位根检验（ADF test）

使用格式：

h=adffuller（Series，maxlag=None，regression='c'，autolag='AIC'，store=False，regresults=False）

输入参数Series为一维观测值序列，返回值依次为adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore。

（5）diff（）

功能：对观测值序列进行差分计算

使用格式：

D.diff（）D为Pandas的DataFrame或Series。

（6）arima

功能：设置时序模式的建模参数，创建ARIMA时序模型

使用格式：

arima=ARIMA（data，（p，1，q））.fit（）

data参数为输入的时间序列，p、q为对应的阶，d为差分次数。

（7）summary（）/summary2（）

功能：生成已有模型的报告

使用格式：

arima.summary（）/arima.summary2（）

其中，arima为已经建立好的ARIMA模型，返回一份格式化的模型报告，包含模型的系数、标准误差、p值、AIC和BIC等详细指标。

（8）aic/bic/hqic

功能：计算ARIMA模型的AIC、BIC、HQIC指标值

使用格式：

arima.aic/arima.bic/arima.hqic

其中，arima为已经建立好的ARIMA模型，返回值是Model时序模型得到的AIC、BIC和HQIC指标值。

（9）forecast（）

功能：用得到的时序模型进行预测

使用格式：

a，b，c=arima.forecast（num）

输入参数num为要预测的天数，arima为已经建立好的ARIMA模型。a为返回的预测值，b为预测的误差，c为预测置信区间。

（10）acorr_ljungbox（）

功能：检测是否为白噪声序列

使用格式：

acorr_ljungbox（data，lags=1），

输入参数data为时间序列数据，lags为滞后数，返回统计量和p值。