关注前沿技术,分享热点话题,QCon全球软件开发大会三站同启,重磅回归!立即查看 了解详情

携程如何基于ARIMA时序分析做业务量的预测

2020 年 10 月 04 日

携程如何基于ARIMA时序分析做业务量的预测

一、 前言

时间序列分析是统计学科的一个重要分支。它主要是通过研究随着时间的推移事物发展变化过程中的规律,来进行事物未来发展情况的预测。在我们的日常生活中,股票的价格走势,奶茶店每天的销售额,一年的降雨量分布,河水四季的涨落情况等都属于时间序列。时间序列的分析深入诸多行业。

时间序列的分类:

图 1
  • 根据指标的平稳性,分为平稳时间序列和非平稳时间序列;
  • 根据指标的性质分类,分为总量指标时间序列,相对指标和平均指标时间序列;
  • 根据指标的时间属性分类,分为时期指标时间序列,时点指标时间序列;

时期指标时间序列是可以相加的,并且相加是有意义的,比如每天的订单量,一个月的订单量直接将这个月对应的每天的订单量相加即可。时点指标时间序列是不可以相加的,反映的是某一时间点达到的水平,比如每天库存量,库存量相加是没有统计意义的,每月总库存量不等于每天库存量加和。

对于互联网公司而言,业务量是公司经营关注的重要指标之一。实际情况的复杂性给业务量的分析预测带来了许多挑战:

  • 具有业务特征的周期性影响
  • 节假日等特定时序节点的变异
  • 地域差异,空间的相互作用
  • 受到库存、实际市场容量的影响
  • 其他外生变量,不可控自然或社会因素

对于时间序列的分析,例如订单量,话务量,库存管理等,实现的方式有 ANN,RNN,LR,ARIMA,Prophet 等,这里我们重点关注 ARIMA 分析方法。

二、 时间序列分析实践

2.1 ARIMA 模型简介

ARMA 模型的全称是自回归移动平均模型,可以说是目前最常用的拟合平稳序列的模型。

ARMA 模型由两部分组成:

p 阶自回归模型 AR§

φ0=0φ0=0 时,自回归模型又称为中心化 AR§模型。非中心化的 AR§序列也可以转化(通过平移)为中心化的 AR§模型。

AR 模型将某时刻 t 的值用过去若干时刻 t-1 到 t-p 的值通过线性组合以及噪声来表示。

q 阶移动平均模型 MA(q)

μ=0μ=0 时,模型 MA(q) 称为中心化 MA(q) 模型,对于非中心化的 MA(q) 模型只要做简单的位移就可以转化为中心化的 MA(q) 模型。

MA 模型是通过历史点的噪声线性组合来表示当前时刻的值。

ARMA 模型其实就是 AR§和 MA(q) 的组合:

同样的,当 φ0=0φ0=0 时该模型称为中心化的 ARMA(p,q) 模型。他结合了两个模型的特点,AR 模型处理当前数据与后期数据之间的关系,MA 则处理随机变动的影响。

对于平稳时间序列可以采用 ARMA 模型直接进行拟合,但是实际场景中,我们的时间序列都是有趋势的,即一般时序为非平稳的,所以需要做平稳处理,其中最常用的是差分处理,使得时序平稳后进行 ARMA 分析。

这一过程其实就是 ARIMA,在原始非平稳时间序列基础上做一阶或二阶差分处理后的平稳时间序列上应用 ARMA 模型。ARIMA(p,d,q) 模型是在 ARMA(p,q) 两元组阶数基础上增加差分 d 的三元组阶数模型。

2.2 ARIMA 模型实践分析步骤

图 2

具体实现以 python 为例。

Step1、读取时间序列

复制代码
df = pd.read_csv('testdata.csv', encoding='gbk', index_col='ddate')
#时间序列索引转换为日期格式
df.index = pd.to_datetime(df.index)
#指标量转为 float 类型
df['cnt'] = df['cnt'].astype(float)
plt.figure(facecolor='white',figsize=(20,8))
plt.plot(df.index,df['cnt'],label='Time Series')
plt.legend(loc='best')
plt.show()

图 3

Step2、时间序列平稳性检验

什么是平稳?

平稳分为严平稳和宽平稳,严平稳保证时间序列的任何有限维分布对于时间的平移是不变的,比如高斯白噪声就是严平稳序列;宽平稳则要求协方差结构随时间的平移而不变,或均值和方差是不变的。

为什么需要平稳?

ARIMA 包含了 AR 模型,AR 模型的实质是用历史时间点数据预测当前时间点对应的值。这就要求序列的相关性不会随着时间变化而变化。

复制代码
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
dftest = adfuller(timeseries, autolag='AIC')
return dftest[1]

原始时间序列平稳性检验未通过 (0.94)。从图 3 也可以看到,时间序列具有明显的上升趋势,所以需要尝试对时间序列进行差分处理,再次检验其平稳性。

Step3、差分处理后检验平稳性

复制代码
pred_day = 7
train_start = datetime(2017,3,1)
train_end = datetime(2019,8,16)
pred_start = train_end+timedelta(1)
pred_end = train_end+timedelta(pred_day)
train_diff=df[train_start:train_end]
train_diff['cnt']=train_diff.diff()
print(test_stationarity(train_diff['cnt'][train_start+timedelta(1):train_end]))
plt.figure(facecolor='white',figsize=(20,8))
plt.plot(train_diff.index,train_diff['cnt'],label='Time Series after diff')
plt.legend(loc='best')
plt.show()

图 4

差分后的时序平稳性检验值 9.51*e(-15),说明差分后时序已经是平稳时间序列了,可以应用 ARIMA 模型。

Step4、画出 ACF 和 PACF 图

自相关函数 ACF,反映了两个点之间的相关性。

偏自相关函数 PACF 则是排除了两点之间其他点的影响,反应两点之间的相关性。比如:在 AR(2) 中,即使 y(t-3) 没有直接出现在模型中,但是 y(t) 和 y(t-3) 之间也相关。

复制代码
import statsmodels.api as sm
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(train_diff['cnt'][1:], lags=20, ax=ax1)
ax1.xaxis.set_ticks_position('bottom')
fig.tight_layout()
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(train_diff['cnt'][1:], lags=20, ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
fig.tight_layout()
plt.show()

图 5

严格来看,ACF 和 PACF 显示存在一定程度的拖尾和振荡。但是,ACF 和 PACF 在 3 阶后有骤降和平稳的趋势,考虑到是短期预测的场景,可进一步结合预测效果和模型检验来进行判断。

Step5、ARIMA 模型定阶

虽然 ACF 和 PACF 为我们提供了选择模型参数的参考依据,但是一般实际情况中,我们总会需要通过模型训练效果确定最终采用的参数值。在 ARMA 模型中,我们通常采用 AIC 法则(赤池信息准则,AIC=2k-2ln(L),k 为模型参数个数,n 为样本数量,L 为似然函数)。AIC 鼓励数据拟合的有良性但是尽量避免出现过拟合的情况。所以实际操作中,我们会选择模型 AIC 值最小的那组参数。

复制代码
#定阶
warnings.filterwarnings("ignore") # specify to ignore warning messages
pmax = 8
qmax = 8
aic_matrix = [] #aic 矩阵
for p in range(1,pmax+1):
tmp = []
for q in range(1,qmax+1):
try: #存在部分报错,所以用 try 来跳过报错。
model = ARIMA(endog=df['cnt'],order=(p,1,q))
results = model.fit(disp=-1)
tmp.append(results.aic)
print('ARIMA p:{} q:{} - AIC:{}'.format(p, q, results.aic))
except:
tmp.append(None)
aic_matrix.append(tmp)
aic_matrix = pd.DataFrame(aic_matrix) #从中可以找出最小值
p,q = aic_matrix.stack().idxmin() #先用 stack 展平,然后用 idxmin 找出最小值位置。
print(u'AIC 最小的 p 值和 q 值为:%s、%s' %(p+1,q+1))

图 6

由于时间序列为 1 阶差分平稳时间序列,所以模型参数 d=1,根据 AIC 最小原则得到 p=7,q=7。

Step6、模型测试与优化

将训练得到的参数加入模型,分析模型效果。

复制代码
model = ARIMA(endog=df['cnt'], order=(p,1,q)) #建立 ARIMA(7, 1,7) 模型
result_ARIMA = model.fit(disp=-1,method='css')
predict_diff=result_ARIMA.predict()
#一阶差分还原
df_shift=df['cnt'].shift(1)
predict=predict_diff+df_shift
plt.figure(figsize=(18,5),facecolor='white')
predict[train_start+timedelta(p+1):train_end].plot(color='blue', label='Predict')
df['cnt'][train_start+timedelta(p+1):train_end].plot(color='red', label='Original')
err=sum(np.sqrt((predict[train_start+timedelta(p+1):train_end]-df['cnt'][train_start+timedelta(p+1):train_end])**2)/df['cnt'][train_start+timedelta(p+1):train_end])/df['cnt'][train_start+timedelta(p+1):train_end].size
plt.legend(loc='best')
plt.title('Error: %.4f'%err)
plt.show()

图 7

用训练好的模型进行未来预测。

复制代码
y_forecasted =result_ARIMA.forecast(steps=pred_day, alpha=0.01)[0] #作为期 7 天的预测
y_truth = df[pred_start:pred_end]['cnt']
# 均方根误差 #平均错误率
mse = np.sqrt( ((y_forecasted - y_truth) ** 2) ).mean()
error_rate = ( abs(y_forecasted - y_truth)/y_truth ).mean()
print('\nThe Mean Error rate of our forecasts is {}'.format(round(error_rate, 4)))

模型预测误差 8.58%(【偏差 / 真实值】的均值),结果并不是太理想,所以我们需要对模型进行优化,考虑是因为指标受到了节假日和周的影响,所以在模型的外生变量里面我们加入节假日和周的识别参数。

加入 exog 外生变量后,需要重新定阶,重新训练模型,步骤与上类似。优化后的预测误差 1.77%,相比之前有了很大程度的提升。

图 8

Step7、模型检验

用模型残差来检验模型的合理性。

复制代码
resid = result_ARIMA_improve.resid #赋值
plt.figure(figsize=(12,8))
qqplot(resid,line='q',fit=True)
#利用 D-W 检验, 检验残差的自相关性
print('D-W 检验值为{}'.format(durbin_watson(resid.values)))

图 9

通过图 9 的 qq 图可以看出,残差基本满足了正态分布。D-W 检验结果值为 1.99,接近于 2,说明残差序列不存在自相关性,即模型较好。

三、总结与展望

  • 对于时间序列的分析一定做好前期评估工作,直观的图表分析会助力我们的决策。多探索优秀的开源工具库,往往会使我们事半功倍。
  • 模型选择至关重要,明确模型的适用场景,根据自身的时序选择适合的模型分析。
  • ARIMA 模型在短时间内的预期效果还算可以,但是长时间比如未来一年的预测不太适用,因为偏差会逐渐增大。
  • 现实中的复杂场景,单一模型很难解决,需要考虑多模型结合的方式实现分析预测。

作者介绍

June,携程数据分析经理,对数仓搭建,数据治理,数据分析等方面有较浓厚的兴趣。

本文转载自公众号携程技术(ID:ctriptech)。

原文链接

携程如何基于 ARIMA 时序分析做业务量的预测

2020 年 10 月 04 日 10:00 1668

评论

发布
暂无评论
发现更多内容

向 "忙忙碌碌泯然众人" say no

Sicolas Flamel

时间分配 时间管理

配置微软Azure大数据HDInsight云集群

虚拟世界的懒猫

microsoft 大数据 微软 azure hdinsght

跬步贴|5分钟搞定缓存击穿问题

架构师跬步营

张小龙 的 22 年和微信的 8 年

池建强

微信 张小龙

消息队列Kafka - Kafka中的选举

Java收录阁

kafka

管理信息系统课程基础知识

Sicolas Flamel

程序员,你最重要的选择是和谁结婚,你最重要的能力是赚钱,钱和女友两手抓

陆陆通通

程序员 赚钱 女朋友 找对象

idea+spring4+springmvc+mybatis+maven实现简单增删改查CRUD

虚拟世界的懒猫

spring mybatis ssm springmvc java8

容器日志采集利器:Filebeat深度剖析与实践

傅轶

Kubernetes 容器 云原生 日志 Filebeat

当 Redis 发生高延迟时,到底发生了什么

程序员历小冰

redis Linux 延迟

Hadoop集群搭建-02安装配置Zookeeper

虚拟世界的懒猫

hadoop zookeeper centos7

【终于解决】ubuntu19安装nvidia驱动后屏幕亮度默认最亮不可调节

虚拟世界的懒猫

ubuntu 英伟达

思考如何节省时间,节省出时间进行思考

伯薇

思考 时间管理 思考力 工作效率 提升效率

万字破解云原生可观测性

谭建

云原生 APM 可观测性 链路追踪 Skywalking

Hadoop集群搭建-03编译安装hadoop

虚拟世界的懒猫

hadoop centos7

安装VMware16兼容Hyper-v+WSL2+Docker+解决0x80370102报错

虚拟世界的懒猫

Docker vmware vm hyper-v WSL2

你真的理解 Java 的基础数据类型吗

Rayjun

Java

为什么开源是基础软件的未来

顾钧

开源 基础软件

机器学习中常用的处理手段

子夜

深度学习

centos6搭建NEXUSphp pt私人种子站

虚拟世界的懒猫

centos nexusphp pt bt

从“成为作者”到“立即创作”:开启你的“写作极客”生活

岛乾坤

写作

Hadoop集群搭建-05安装配置YARN

虚拟世界的懒猫

hadoop

让你写出来的代码像诗一样优美!《Java开发手册》PDF下载

Kareza

Java 阿里巴巴 Java规范 Java开发手册

Hadoop集群搭建-04安装配置HDFS

虚拟世界的懒猫

hadoop

利用Translate ToolKit 2.5.0 API构建Flask web app

虚拟世界的懒猫

Python nginx flask uwgsi translate

选赵敏还是选小昭,这可真是个问题 | Decision Tree

张利东

Python 机器学习 算法 决策树

程序员陪娃漫画系列——喂药

孙苏勇

程序员 生活 程序员人生 陪伴 漫画

如何无缝的将Flutter引入现有应用?

稻子

flutter ios android 开源 移动应用

Zookeeper选举机制

tunsuy

zookeeper 开源 源码分析 分布式协同

做好仓储控制系统(WCS)的关键

阿喜伯

仓储控制系统 WCS

Hadoop集群搭建-01前期准备

虚拟世界的懒猫

hadoop hdfs mapreduce zookeeper centos

携程如何基于ARIMA时序分析做业务量的预测-InfoQ