时间序列：ARIMA

1 ARIMA模型

1.1 ARIMA模型是什么：

全称差分自回归移动平均模型 （Autoregressive Integrated Moving Average Model）

• 它其实可以拆分为：AR自回归模型、I差分、MA移动平均模型

• 自回归模型（AR模型） ——》 自己的数据跟自己的数据进行一个相关性的回归分析！

• 描述当前值与历史值之间的关系，用变量自身的历史事件数据对自身进行预测

• 自回归模型必须满足平稳性的要求

• 这里有个模型的参数p
  

• 移动平均模型（MA模型）

• 移动平均模型关注的是自回归模型中的误差项的累加

• 移动平均法能有效地消除预测中的随机波动

• 这里有个模型的参数q
  

• 差分法

• 时间序列在t与t-1时刻的差值

• 差分参数d

• 能够使时间序列变得平稳

• 原时间序列
  

• 一阶差分
  

• 二阶差分

• 总结
• ARIMA（p,d,q），原理：将非平稳序列转化为平稳序列，然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。

1.2 ARIMA模型建模流程

• 白噪声检验（纯随机性检验）

• 什么是纯随机序列？

• 答：如果一个序列是纯随机序列，那么序列值之间没有任何关系，且序列没有规律可循。

• 用途：

• 建模之前，检验数据是否满足白噪声检验，非白噪声才能进一步建模

• 建模后，检验残差是否满足白噪声检验，通过检验，建模才能成立

2 ARIMA模型的应用

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

from matplotlib.pylab import style #自定义图表风格
style.use('ggplot')

# 解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['Simhei']
# 解决坐标轴刻度负号乱码
plt.rcParams['axes.unicode_minus'] = False

#pip install statsmodels
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf  #自相关图、偏自相关图
from statsmodels.tsa.stattools import adfuller as ADF #平稳性检验
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
## 白噪声检验的
# import statsmodels.api as sm #D-W检验,一阶自相关检验
# from statsmodels.graphics.api import qqplot #画QQ图,检验一组数据是否服从正态分布
from statsmodels.tsa.arima_model import ARIMA

2.1导入数据

sale=pd.read_excel('./arima_data.xls',index_col='日期')
sale.head()

2.2画时序图

• 解读：具有单调递增趋势，则是非平稳序列。

plt.figure(figsize=(10,5))
sale.plot()
plt.show()

2.3平稳性检验

• 方法：单位根检验（这里采用ADF检验）
• 时间序列单位根的检验就是对时间序列平稳性的检验，时间序列如果存在单位根，则说明它为非平稳序列。
• 作出的假设检验为：
• HO：时间序列为非平稳序列
• H1：时间序列为平稳序列

#方法：单位根检验

print('原始序列的ADF检验结果为：',ADF(sale.销量))

#0.9983759421514264表示的是p值（个人理解是犯错的概率）
#解读：P值大于显著性水平α（0.05），没有足够的证据去拒绝原假设（非平稳序列），说明原始序列是非平稳序列。

原始序列的ADF检验结果为： (1.813771015094526, 0.9983759421514264, 10, 26, {'1%': -3.7112123008648155, '5%': -2.981246804733728, '10%': -2.6300945562130176}, 299.4698986602418)

2.4差分法

• 差分法能够使时间序列变为平稳序列
• ARIMA的d = 1

#差分法
d1_sale=sale.diff(periods=1, axis=0).dropna()

#时序图
plt.figure(figsize=(10,5))
d1_sale.plot()
plt.show()
#解读：在均值附件比较平稳波动

sale.plot(color='b')

# #自相关图
# plot_acf(d1_sale,lags=34).show()
# #解读：有短期相关性，但趋向于零。

#平稳性检验
print('原始序列的ADF检验结果为：',ADF(d1_sale.销量))

#解读：P值小于显著性水平α（0.05），拒绝原假设（非平稳序列），说明一阶差分序列是平稳序列。

2.5白噪声检验

• Ljung-Box Test（LB法）

print('一阶差分序列的白噪声检验结果为：',acorr_ljungbox(d1_sale,lags=1))#返回统计量、P值

# 参数lags表示的是置后值[1,2,3,4,5,6],一般1就可以了，不放心可以多填几个
# 假设P值<0.05的、那么该模型是白噪声，也就是有规律的，可以用于建模的

#解读：p值小于0.05，拒绝原假设（纯随机序列），说明一阶差分序列是非白噪声。

一阶差分序列的白噪声检验结果为： (array([11.30402222]), array([0.00077339]))

2.6模型构建

2.6.1 图检验（偏主观）

• 模型ARIMA(p,d,q)，三个参数，我们需要进行定阶

• 自相关函数ACF：有序的随机变量与其自身比较，自相关函数反映了同一序列在不同时序的取值之间的相关性

• 偏自相关函数（PACF）：对于一个平稳的AR模型，求出滞后k个自相关函数ACF时，实际上不是得到x(t)与x(t-k)之间单纯的相关关系，它中间还会受到k-1个随机变量x(t-1),x(t-2),…x(t-k+1)的影响。为剔除这些影响统计学家研究出PACF偏自相关函数。

• PACF是严格两个变量之间的相关性。

• 截尾：落在置信区间内

d1_sale=sale.diff(periods=1, axis=0).dropna()

#自相关图
plot_acf(d1_sale,lags=20).show()

#解读：有短期相关性，但趋向于零。

#偏自相关图
plot_pacf(d1_sale,lags=20).show()


#解读：自相关图，0阶截尾；偏自相关图，1阶截尾。则ARIMA(p,d,q)=ARIMA(0,1,1)

2.6.2 相关评估指标检验：BIC和AIC法

• 此处我们介绍一下常用的两个模型选择方法——赤池信息准则（Akaike Information Criterion，AIC）和贝叶斯信息准则（Bayesian Information Criterion，BIC）。

• AIC是衡量统计模型拟合优良性的一种标准，由日本统计学家赤池弘次在1974年提出，它建立在熵的概念上，提供了权衡估计模型复杂度和拟合数据优良性的标准。

• BIC（Bayesian InformationCriterion）贝叶斯信息准则与AIC相似，用于模型选择，1978年由Schwarz提出。训练模型时，增加参数数量，也就是增加模型复杂度，会增大似然函数，但是也会导致过拟合现象，针对该问题，AIC和BIC均引入了与模型参数个数相关的惩罚项，BIC的惩罚项比AIC的大，考虑了样本数量，样本数量过多时，可有效防止模型精度过高造成的模型复杂度过高。

• 值越小越好

• 参数调优的方法非常多，用不同方法得出的结论可能不同

pmax=int(len(d1_sale)/10) #一般阶数不超过length/10
qmax=int(len(d1_sale)/10) #一般阶数不超过length/10
pmax
qmax

• BIC法

bic_matrix=[]
for p in range(pmax+1):#0，1，2，3
    tmp=[]
    for q in range(qmax+1):#0，1，2，3
        try:
            tmp.append(ARIMA(sale,(p,1,q)).fit().bic)
        except:
            tmp.append(None)
    bic_matrix.append(tmp)

bic_matrix

[[432.0684724517514, 422.51008220208723, 426.0889106675776, 426.595507416633],
 [423.62827614854785, 426.07360131952146, None, None],
 [426.77482379014015, 427.3958203685418, 431.00352092490533, None],
 [430.3175243609237,
  431.92488894885025,
  434.76171244384716,
  436.4781092333494]]

#行是p，列是q
df_bic

df_bic = pd.DataFrame(bic_matrix)

df_bic.min().min()

422.51008220208723

df_bic==df_bic.min().min()

• 总结，用BIC得到的最优参数为p=0,q=1

• AIC法

pmax=int(len(d1_sale)/10) #一般阶数不超过length/10
qmax=int(len(d1_sale)/10) #一般阶数不超过length/10

aic_matrix=[]
for p in range(pmax+1):
    tmp=[]
    for q in range(qmax+1):
        try:
            tmp.append(ARIMA(sale,(p,1,q)).fit().aic)
        except:
            tmp.append(None)
    aic_matrix.append(tmp)
df_aic=pd.DataFrame(aic_matrix)

df_aic==df_aic.min().min()

• 总结，用AIC得到的最优参数为p=0,q=1

2.7 建模与预测

• 建模

model=ARIMA(sale,(0,1,1)).fit()

• 残差检验
• 残差是什么？其实就是y_true - y_predict
• 解读：残差是白噪声

resid=model.resid

print('残差序列的白噪声检验结果为：',acorr_ljungbox(resid,lags=1))#返回统计量、P值

残差序列的白噪声检验结果为： (array([0.00390451]), array([0.95017574]))

【补充说明】

如果残差不是白噪声，可以考虑取另外一组p、q值；或者该数据不适应ARIMA模型。

• 预测

print('未来7天的销量预测：')
model.forecast(7) #预测、标准差、置信区间

未来7天的销量预测：

(array([4873.9667493 , 4923.92317644, 4973.87960359, 5023.83603073,
        5073.79245787, 5123.74888501, 5173.70531215]),
 array([ 73.08574293, 142.32679918, 187.542821  , 223.80281869,
        254.95704265, 282.69857718, 307.95109593]),
 array([[4730.72132537, 5017.21217324],
        [4644.96777602, 5202.87857687],
        [4606.30242887, 5341.4567783 ],
        [4585.19056646, 5462.48149499],
        [4574.08583666, 5573.49907907],
        [4569.66985525, 5677.82791477],
        [4570.13225512, 5777.27836918]]))

forecast=pd.Series(model.forecast(7)[0],index=pd.date_range('2015-2-7',periods=7,freq='D'))
forecast

2015-02-07    4873.966749
2015-02-08    4923.923176
2015-02-09    4973.879604
2015-02-10    5023.836031
2015-02-11    5073.792458
2015-02-12    5123.748885
2015-02-13    5173.705312
Freq: D, dtype: float64

data=pd.concat((sale,forecast),axis=0)
data.columns=['销量','未来7天销量']
plt.figure(figsize=(10,5))
data.plot()
plt.show()

MSE，sum(yture - y_pre)**0.5)