1 前言

​ 从时间上看，订单量时间序列有两个明显的特征：

​ 1）周期性。每天订单量的变化趋势都大致相同，午高峰和晚高峰订单量集中；

​ 2）实时性。当天的订单量可能会受天气等因素影响，呈现整体的上涨或下降；

​ 预测可以反映未来司机成单的情况，能给运营部门及时调整有效的运营策略。预测又有好几种方向：基于订单总额的预测，基于乘客目的地预测，基于上车地点的供需预测等，这里阐述订单总额未来七天的预测。

2 数据建模基本步骤

1. 获取被观测系统时间序列数据；

2. 对数据绘图，观测是否为平稳时间序列；对于非平稳时间序列要先进行d阶差分运算，化为平稳时间序列；

3. 经过第二步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF ，通过对自相关图和偏自相关图的分析，得到最佳的阶层 p 和阶数 q；

4. 由以上得到的d、q、pd、q、p ，得到ARIMA模型。然后开始对得到的模型进行模型检验

3 ARIMA实战解剖

​ 原理大概清楚，实践却还是会有诸多问题。相比较R语言，Python在做时间序列分析的资料相对少很多。下面就通过Python语言详细解析后三个步骤的实现过程。

文中使用到这些基础库：

pandas,numpy,scipy,matplotlib,statsmodels,pandas, 对其调用如下：

from __future__ import print_functionimport pandas as pdimport numpy as npfrom scipy import  statsimport matplotlib.pyplot as pltimport statsmodels.api as smfrom statsmodels.graphics.api import qqplotfrom statsmodels.tsa.stattools import adfuller as ADF

3.1 获取数据

​ 对乘用车日运营报表订单总额数据，进行分析。

数据如下：

日期	订单总额
2018/1/1	22093.02
2018/1/2	18917.3
2018/1/3	6539.3
2018/1/4	19442.9
2018/1/5	24425.74
2018/1/6	29198.1
2018/1/7	26611.8
2018/1/8	25389.75
2018/1/9	22802.62
2018/1/10	31000.7
2018/1/11	28794.53
2018/1/12	28746.15
2018/1/13	30382
2018/1/14	31726.1
2018/1/15	32041.8
2018/1/16	30514.5
2018/1/17	31783.9
2018/1/18	29976.9
2018/1/19	33292.82
2018/1/20	35055.3
2018/1/21	33054.2
2018/1/22	33000.27
2018/1/23	34889.85
2018/1/24	33495.7
2018/1/25	33960.66
2018/1/26	35901.61
2018/1/27	2276.5
2018/1/28	1268.7
2018/1/29	32588.14
2018/1/30	33881
2018/1/31	29342.02
2018/2/1	24991.15
2018/2/2	30995.08
2018/2/3	31219.9
2018/2/4	31040.2
2018/2/5	31152.5
2018/2/6	31190.76
2018/2/7	32778.69
2018/2/8	31938.7
2018/2/9	31857.8
2018/2/10	31969.73
2018/2/11	29406.12
2018/2/12	26877.6
2018/2/13	24472.7
2018/2/14	9663.4
2018/2/15	142.2
2018/2/16	108.6
2018/2/17	321.9
2018/2/18	16647.8
2018/2/19	18353.85
2018/2/20	21131.1
2018/2/21	20511.5
2018/2/22	25521.6
2018/2/23	25701
2018/2/24	28140.4
2018/2/25	31535
2018/2/26	27903.9
2018/2/27	29730.48
2018/2/28	28317.9
2018/3/1	12704.4
2018/3/2	13736.45
2018/3/3	13093.1
2018/3/4	11186.6
2018/3/5	7279.68
2018/3/6	9042.7
2018/3/7	8524.87
2018/3/8	10346.3
2018/3/9	9219.9
2018/3/10	9931.8
2018/3/11	9772.8
2018/3/12	9100.1
2018/3/13	7822.2
2018/3/14	8863.92
2018/3/15	20644.8
2018/3/16	29468.6
2018/3/17	31152.9
2018/3/18	29558.9
2018/3/19	28728.9
2018/3/20	26359.25
2018/3/21	27735.05
2018/3/22	31034.86
2018/3/23	33682.55
2018/3/24	38766.6
2018/3/25	36706.35
2018/3/26	29303.9
2018/3/27	26164.2
2018/3/28	28287.7
2018/3/29	26742.22
2018/3/30	30283.25
2018/3/31	32708.4
2018/4/1	39168.62
2018/4/2	40197.3
2018/4/3	39577.86
2018/4/4	38930.96
2018/4/5	33691.44
2018/4/6	31831.3
2018/4/7	36942.8
2018/4/8	35314.52
2018/4/9	33843.3
2018/4/10	35063.18
2018/4/11	35706.36
2018/4/12	40346.94
2018/4/13	42388.1
2018/4/14	41396.29
2018/4/15	42374.3
2018/4/16	16489.5
2018/4/17	31132.9
2018/4/18	32511
2018/4/19	35183.22
2018/4/20	65664.7
2018/4/21	73029.91
2018/4/22	73738.58
2018/4/23	60052.2
2018/4/24	66994.34
2018/4/25	70996.55
2018/4/26	66477.66
2018/4/27	82832.23
2018/4/28	88265.1
2018/4/29	85795.8
2018/4/30	79327.25
2018/5/1	87192.7
2018/5/2	81075.1
2018/5/3	84212.59
2018/5/4	98093.35
2018/5/5	84365.13
2018/5/6	82469.8
2018/5/7	83517.96
2018/5/8	85643.2
2018/5/9	93718.3
2018/5/10	94353.3
2018/5/11	103593.85
2018/5/12	106824.98
2018/5/13	107588.58
2018/5/14	98374.04
2018/5/15	104670.4
2018/5/16	118165.44
2018/5/17	60852.3
2018/5/18	48933.27
2018/5/19	49102.2
2018/5/20	44945.4
2018/5/21	34141.62
2018/5/22	32317.6
2018/5/23	29779.3
2018/5/24	25971

discfile = 'D:/PycharmProjects/OrderForecast/datafile/data.xls'forecastnum = 7#读取数据，指定日期列为指标，Pandas自动将“日期”列识别为Datetime格式data = pd.read_excel(discfile, index_col=u'日期')#用来正常显示中文标签plt.rcParams['font.sans-serif'] = ['SimHei']#用来正常显示负号plt.rcParams['axes.unicode_minus'] = Falsedata.plot()

订单走势如下图：

3.2 时间序列的差分dd

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

#时间序列的差分ddta=data[u'订单总额']fig = plt.figure(figsize=(12,8))ax1= fig.add_subplot(111)diff1 = dta.diff(2)diff1.plot(ax=ax1)print(u'原始序列的ADF检验结果为：', ADF(dta))

一阶差分的时间序列的均值和方差已经基本平稳，不过我们还是可以比较一下二阶差分的效果:

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

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

3.3 合适的p,qp,q

得到一个平稳的时间序列后，接来下就是选择合适的ARIMA模型，即ARIMA模型中合适的p,qp,q。

第一步要先检查平稳时间序列的自相关图和偏自相关图。

dta= dta.diff(1)#已经知道要使用一阶差分的时间序列，之前判断差分的程序可以注释掉fig = plt.figure(figsize=(12,8))ax1=fig.add_subplot(211)fig = sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax1)ax2 = fig.add_subplot(212)fig = sm.graphics.tsa.plot_pacf(dta,lags=40,ax=ax2)

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

通过两图观察得到：

• 自相关图显示滞后有三个阶超出了置信边界；
• 偏相关图显示在滞后1至2阶（lags 1,2...7）时的偏自相关系数超出了置信边界，从lag 7之后偏自相关系数值缩小至0
  则有以下模型可以供选择：
  1. ARMA(0,1)模型：即自相关图在滞后1阶之后缩小为0，且偏自相关缩小至0，则是一个阶数q=1的移动平均模型；
  2. ARMA(7,0)模型：即偏自相关图在滞后7阶之后缩小为0，且自相关缩小至0，则是一个阶层p=3的自回归模型；
  3. ARMA(7,1)模型：即使得自相关和偏自相关都缩小至零。则是一个混合模型。
     4.还可以有其他供选择的模型
     现在有以上这么多可供选择的模型，我们通常采用ARMA模型的AIC法则。

arma_mod20 = sm.tsa.ARMA(dta,(7,0)).fit()print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)print(sm.stats.durbin_watson(arma_mod20.resid.values))arma_mod30 = sm.tsa.ARMA(dta,(0,1)).fit()print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)print(sm.stats.durbin_watson(arma_mod30.resid.values))arma_mod40 = sm.tsa.ARMA(dta,(7,1)).fit()print(arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic)print(sm.stats.durbin_watson(arma_mod40.resid.values))arma_mod50 = sm.tsa.ARMA(dta,(8,0)).fit()print(arma_mod50.aic,arma_mod50.bic,arma_mod50.hqic)print(sm.stats.durbin_watson(arma_mod50.resid.values))

3001.733127872179 3029.461447568363 3013.5940088028337 1.96725505482845643174.789227240986 3183.698667139714 3178.4095208845374 0.75144605532596593005.171853576973 3034.8699865727326 3017.239499055478 1.9457553708154533001.8331059050256 3031.5312389007854 3013.900751383531 1.9401970059406308

可以看到ARMA(7,0)的aic，bic，hqic值最小，因此是最佳模型。

3.4 模型检验

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

3.4.1 对ARMA(7,0)模型所产生的残差做自相关图

fig = plt.figure(figsize=(12,8))ax1 = fig.add_subplot(211)fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)ax2 = fig.add_subplot(212)fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)

3.4.2 做D-W检验

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

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

DW＝4＜＝＞ρ＝-1　即存在负自相关性 DW＝2＜＝＞ρ＝0　　即不存在（一阶）自相关性 因此，当DW值显著的接近于O或４时，则存在自相关性，而接近于２时，则不存在（一阶）自相关性。这样只要知道ＤＷ统计量的概率分布，在给定的显著水平下，根据临界值的位置就可以对原假设Ｈ０进行检验。

print(sm.stats.durbin_watson(arma_mod20.resid.values))

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

3.4.3 观察是否符合正态分布

这里使用QQ图，它用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布。

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

3.5 模型预测

模型确定之后，就可以开始进行预测了，对未来七天的数据进行预测。

predict_sunspots = arma_mod20.predict('2018-05-23', '2018-05-29, dynamic=True)print(predict_sunspots)fig=plt.figure(figsize=(12, 8))ax3=fig.add_subplot(713)predict_sunspots.plot(ax=ax3)

预测结果如下：

​ 从图形来，预测结果较为合理。