python时间减法_时间序列分析

大家晚上好，我是阿涛。

本期继续炒冷饭，时间序列分析，也称时序模型，将用python实现arima模型，当然还有很多丰富的模型ARCH、GARCH就不过多展示咯。

Reference:《python数据挖掘分析与实战》——代码、数据都来源于这里，需要可联系我

ok，开始之前先放个本文小目录：

一.时间序列模型

二.时间序列的预处理

2.1 平稳性检验

2.2 纯随机性检验

三.平稳时间序列分析

3.1 AR模型

3.2 MA模型

3.3 ARMA模型

3.4 平稳时间序列建模

四.非平稳时间序列分析

4.1 差分运算

4.2 ARIMA模型

4.3 实例：某餐厅销量预测

【1】.时间序列模型

基于菜品历史销售数据做好餐饮销售预测，以便减少菜脱销，或者避免因备料不够而导致的生产延误。                             

 【2】.时间序列的预处理

拿到一个观察值序列后，首先要对它的纯随机性和平稳性进行检验，这两个重要的检验称为序列的预处理。根据检验结果可以将序列分为不同的类型，对不同类型的序列会采取不同的分析方法。

纯随机序列又叫白噪声序列，序列的各项之间没有任何相关关系，序列在进行完全无序的随机波动，可以终止对该序列的分析。白噪声序列是没有信息可提取的平稳序列。

对于平稳非白噪声序列，它的均值和方差是常数，现已有一套非常成熟的平稳序列的建模方法。通常是建立一个线性模型来拟合该序列的发展，借此提取该序列的有用信息。所以我们只分析非白噪声的序列。才有价值。

ARMA模型是最常用的平稳序列拟合模型。对于非平稳序列，由于它的均值和方差不稳定，处理方法一般是将其转变为平稳序列，这样就可以应用有关平稳时间序列的分析方法，如建立ARMA模型来进行相应的研究。如果一个时间序列经差分运算后具有平稳性，则称该序列为差分平稳序列，可以使用ARIMA模型进行分析。

2.1 平稳性检验

2.1.1 平稳时间序列的定义

如果时间序列{Xt,t∈T}在某一常数附近波动且波动范围有限，即有常数均值和常数方差，并且延迟k期的序列变量的自协方差和自相关系数是相等的，或者说延迟k期的序列变量之间的影响程度是一样的，则称{Xt,t∈T}为平稳序列。

2.1.2 平稳性的检验：

对序列平稳性的检验有两种方法：一种是根据时序图和自相关图的特征做出判断的图检验，该方法操作简单、应用广泛，缺点是带有主观性；另一种是构造检验统计量进行检验的方法，目前最常用的方法是单位根检验。

①时序图检验

根据平稳时间序列的均值和方差都为常数的性质，平稳序列的时序图显示该序列值始终在一个常数附近随机波动，而且波动的范围有界；如果表现出明显的趋势性或者周期性，那它通常不是平稳序列。

②自相关图检验

平稳序列具有短期相关性，这表明对平稳序列而言通常只有近期的序列值对现时值的影响比较明显，间隔越远的过去值对现时值的影响越小。随着延迟期数k的增加，平稳序列的自相关系数ρk(延迟k期)会比较快地衰减趋向于零，并在零附近随机波动。而非平稳序列的自相关系数衰减的速度比较慢，这就是利用自相关图进行平稳性检验的标准。

③单位根检验

单位根检验是指检验序列中是否存在单位根，存在单位根就是非平稳时间序列了。

2.2 纯随机性检验

如果一个序列是纯随机序列，那么它的序列值之间应该没有任何关系，即满足γ(k)=0，k≠C，这是一种理论上才会出现的理想状态，实际上纯随机序列的样本自相关系数不会绝对为零，但是很接近零，并在零附近随机波动。

纯随机性检验也称白噪声检验，一般是构造检验统计量来检验序列的纯随机性。常用的检验统计量有Q统计量和LB统计量，由样本各延迟期数的自相关系数可以计算得到检验统计量，然后计算出对应的p值，如果p值明显大于显著性水平α，则表示该序列不能拒绝纯随机的原假设，可以停止对该序列的分析。

【3】.平稳时间序列分析

自回归移动平均模型(Autoreg Ressive Moving Average Model，简称ARMA模型)，是目前最常用的拟合平稳序列的模型。它又可以细分为AR模型、MA模型和ARMA这3大类，都可以看作是多元线性回归模型。

1.AR模型

(1)AR(p)

具有如下结构的模型称为p阶自回归模型，简称为AR(p)

Xt受p期序列值影响。

(2)平稳AR模型的性质：

均值：常数均值

方差：常数方差

自相关系数(ACF)：拖尾

偏自相关系数(PACF)：p阶截尾

2.MA模型

(1) MA(q)

具有如下结构的模型称为q阶自回归模型，简称为MA(q)

 

Xt受过去q期误差项的影响。

(2)平稳MA模型的性质：

均值：常数均值

方差：常数方差

自相关系数(ACF)：q阶截尾

偏自相关系数(PACF)：拖尾

3.ARMA模型

(1)ARMA(p,q)模型

Xt受过去p期序列值和q期误差项的共同影响。

(2)平稳ARMA模型的性质：

均值：常数均值

方差：常数方差

自相关系数(ACF)：拖尾

偏自相关系数(PACF)：拖尾

4.平稳时间序列建模

1)计算ACF和PACF。先计算非平稳白噪声序列的自相关系数(ACF)和偏自相关系数(PACF)。

2)ARMA模型识别。也称为模型定阶，根据AR(p)、MA(q)和ARMA(p，q)模型的自相关系数和偏自相关系数的性质选择合适的模型。流程图如下：

【4】.非平稳时间序列分析

前面介绍了平稳时间序列的分析方法。实际上，在自然界中绝大部分序列都是非平稳的。因而非平稳时间序列的分析更普遍、更重要，创造出来的分析方法也更多。非平稳时间序列的分析方法可以分为确定性因素分解的时序分析和随机时序分析两大类。(考虑Y=T+S+C+I / Y=T*S*C*I 对应下列四个因素)

确定性因素分解的方法把所有序列的变化都归结为4个因素(长期趋势、季节变动、循环变动和随机波动)的综合影响，其中长期趋势和季节变动的规律性信息通常比较容易提取，而由随机因素导致的波动则非常难以确定和分析，对随机信息浪费严重，会导致模型拟合精度不够理想。

随机时序分析法的发展就是为了弥补确定性因素分解方法的不足。根据时间序列的不同特点，随机时序分析可以建立的模型有ARIMA模型、残差自回归模型、季节模型(考虑季节因素分解)、异方差模型等。接下来重点介绍ARIMA模型对非平稳时间序列进行建模。

4.1 差分运算

(1)p阶差分相距一期的两个序列值之间的减法运算称为1阶差分运算。

(2)k步差分相距k期的两个序列值之间的减法运算称为k步差分运算。

4.2 ARIMA模型

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

4.3 实例：某餐厅销量预测

4.3.1 原始数据：

4.3.2 检验序列的平稳性：

明显单调递增>说明是非平稳

4.3.3 自相关图：

自相关图长期大于0>说明是具有很强的长期相关性

4.3.4 单位根检验：

P>0.05 说明>非平稳序列(非平稳序列一定不是白噪声序列)

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

一阶差分后序列的时序图在100附近波动，自相关图有很强的短期相关性，单位根P值小于0.05，说明一阶差分后符合平稳序列了。

4.3.6 白噪声检验：

P值<0.05说明是非白噪声

4.3.7 下面开始对模型进行定阶了ARMA(p,q)

第一种方法：人为定阶

由差分后的自相关图可以看出是1阶截尾，由上图(偏自相关图)可以看出是拖尾(仔细看这个图和前面那个差分后的自相关图)，有统计学基础的同学应该是可以知晓的。

所以认为是MA(1)模型拟合1阶差分后的序列，即建立ARIMA(0,1,1)模型。

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

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

由bic矩阵得到最小的值是422.510082>>得到p=0，q=1。

但模型并非唯一，ARIMA(1,1,0)和ARIMA(1,1,1)也能通过检验。通常来说Numpy和Pandas结合已经很强了，只要要得到较为深入的统计模型才会用到statsmodels.

 补充：python主要时序模型算法

 最后代码见文末.

下期预告：

深度学习推荐系统~

 感谢大家看到这里，码字不易，觉得有用的话留下小心心再走~❤

往期推荐阅读 ▽爬虫之scrapy框架极速可视化BI——Tableau用python实现微信自动回复五种聚类算法一览与python实现分类算法一览与python实现

Lasso回归、Ridge回归与ElasticNet回归

End

作者：涛

半壶水全栈工程师，好读书，甚喜之

code:

# -*- coding: utf-8 -*-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 pltplt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号data.plot()plt.show()# 自相关图from statsmodels.graphics.tsaplots import plot_acfplot_acf(data).show()# 平稳性检测from statsmodels.tsa.stattools import adfuller as ADFprint(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_pacfplot_pacf(D_data).show()  # 偏自相关图print(u'差分序列的ADF检验结果为：', ADF(D_data[u'销量差分']))  # 平稳性检测# 白噪声检验from statsmodels.stats.diagnostic import acorr_ljungboxprint(u'差分序列的白噪声检验结果为：', acorr_ljungbox(D_data, lags=1))  # 返回统计量和p值from statsmodels.tsa.arima_model import ARIMA# 定阶data[u'销量'] = data[u'销量'].astype(float) pmax = int(len(D_data)/10)  # 一般阶数不超过length/10qmax = int(len(D_data)/10)  # 一般阶数不超过length/10bic_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)模型print('模型报告为：\n', model.summary2())print('预测未来5天，其预测结果、标准误差、置信区间如下：\n', model.forecast(5))