基于AIRMA模型对订单总额未来七天的预测

Posted 2020-12-30

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_function
import pandas as pd
import numpy as np
from scipy import  stats
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics.api import qqplot
from 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‘] = False
data.plot()

订单走势如下图：

3.2 时间序列的差分dd

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

#时间序列的差分d
dta=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.9672550548284564
3174.789227240986 3183.698667139714 3178.4095208845374 0.7514460553259659
3005.171853576973 3034.8699865727326 3017.239499055478 1.945755370815453
3001.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)

预测结果如下：

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