算法_ARIMA时间序列分析总结

简介

ARIMA算法流程步骤（算法数学推导自行查阅相关论文），本文只讲工程技术和方法。

参考文章
标题：ARIMA模型
地址：http://wiki.mbalib.com/wiki/ARIMA%E6%A8%A1%E5%9E%8B
标题：[python] 时间序列分析之ARIMA
地址：http://blog.csdn.net/u010414589/article/details/49622625

01，什么是ARIMA模型

ARIMA模型全称为自回归移动平均模型(Autoregressive Integrated Moving Average Model,简记ARIMA)，是由博克思(Box)和詹金斯(Jenkins)于70年代初提出的一著名时间序列预测方法，所以又称为box-jenkins模型、博克思-詹金斯法。其中ARIMA（p，d，q）称为差分自回归移动平均模型，AR是自回归, p为自回归项; MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数。

02，ARIMA模型的基本思想

ARIMA模型的基本思想是：将预测对象随时间推移而形成的数据序列视为一个随机序列，用一定的数学模型来近似描述这个序列。这个模型一旦被识别后就可以从时间序列的过去值及现在值来预测未来值。现代统计方法、计量经济模型在某种程度上已经能够帮助企业对未来进行预测。

03，ARIMA模型预测的基本程序

 （一）根据时间序列的散点图、自相关函数和偏自相关函数图以ADF单位根检验其方差、趋势及其季节性变化规律，对序列的平稳性进行识别。一般来讲，经济运行的时间序列都不是平稳序列。
 （二）对非平稳序列进行平稳化处理。如果数据序列是非平稳的，并存在一定的增长或下降趋势，则需要对数据进行差分处理，如果数据存在异方差，则需对数据进行技术处理，直到处理后的数据的自相关函数值和偏相关函数值无显著地异于零。
 （三）根据时间序列模型的识别规则，建立相应的模型。若平稳序列的偏相关函数是截尾的，而自相关函数是拖尾的，可断定序列适合AR模型；若平稳序列的偏相关函数是拖尾的，而自相关函数是截尾的，则可断定序列适合MA模型；若平稳序列的偏相关函数和自相关函数均是拖尾的，则序列适合ARMA模型。
 （四）进行参数估计，检验是否具有统计意义。
 （五）进行假设检验，诊断残差序列是否为白噪声。
 （六）利用已通过检验的模型进行预测分析。

小结：看完上面明白ARIMA干嘛的就行了，简单来说“基于时间序列的回归(数值型预测都看作回归)”。ARIMA已有现成方法了(python,R)，只关注参数怎么填就ok（其实也有填参方法^_^，先不说这个），3参数p，d，q，分别讲解。

04，平稳化

由于ARIMA 模型对时间序列的要求是平稳型。因此对一个非平稳的时间序列时，需先转为平稳序列。

定义平稳序列
恒定的平均数
恒定的方差
不随时间变化的自协方差

不平稳如何转为平稳
a，预测&消除趋势
i:消除趋势的第一个方法是转换。
若曲线有一个显著的趋势。可以通过变换，惩罚较高值而不是较小值。这可以采用平方根,立方跟等等。（其实就是把较大的数值压缩）
ii:提取主趋势
聚合-取一段时间的平均值（月/周平均值）
平滑-取滚动平均数（或指数加权移动平均法）
多项式回归分析-适合的回归模型
iii:差值得新序列。
通过原序列-主趋势的差值得到新序列（主趋势可以看作能独立预测的变量了） 。此时新序列是消除主趋势后的，较原序列更平稳。
b，差分–采用一个特定时间差的差值
c，分解——建立有关趋势和季节性的模型和从模型中删除它们。
分解法有现成的季节代码块调用，见 参考资料《时间序列预测全攻略（附带Python代码）》

如何确定已经平稳
1、绘制滚动统计：我们可以绘制移动平均数和移动方差，观察它是否随着时间变化。随着移动平均数和方差的变化，我认为在任何“t”瞬间，我们都可以获得去年的移动平均数和方差。如：上一个12个月份。但是，这更多的是一种视觉技术。（看图法，均匀的绕0轴震动）
2、DF检验：这是一种检查数据稳定性的统计测试。无效假设：时间序列是不稳定的。测试结果由测试统计量和一些置信区间的临界值组成。如果“测试统计量”少于“临界值”，我们可以拒绝无效假设，并认为序列是稳定的。

05，确定参数d

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

1
2
3
4
5
6
7
8
9
10

dta=[10930,10318,10595,10972,7706,6756,9092,10551,9722,10913,11151,8186,6422,  
6337,11649,11652,10310,12043,7937,6476,9662,9570,9981,9331,9449,6773,6304,9355,  
10477,10148,10395,11261,8713,7299,10424,10795,11069,11602,11427,9095,7707,10767,  
12136,12812,12006,12528,10329,7818,11719,11683,12603,11495,13670,11337,10232,  
13261,13230,15535,16837,19598,14823,11622,19391,18177,19994,14723,15694,13248,  
9543,12872,13101,15053,12619,13749,10228,9725,14729,12518,14564,15085,14722,  
11999,9390,13481,14795,15845,15271,14686,11054,10395]  
dta=pd.Series(dta)  
dta.index = pd.Index(sm.tsa.datetools.dates_from_range('2001','2100'))  
dta.plot(figsize=(12,8))

图1
看起来不错，保险起见，做下二阶的

1
2
3
4

fig = plt.figure(figsize=(12,8))  
ax1= fig.add_subplot(111)  
diff1 = dta.diff(1)  
diff1.plot(ax=ax1)

图2
可以看出二阶差分后的时间序列与一阶差分相差不大，并且二者随着时间推移，时间序列的均值和方差保持不变。这里所谓的一阶差分和二阶差分就是需要的参数d了，此处一阶已经满足，故d=1（此时应做df检验，暂忽略）。

06，确定参数p和q

利用前面得到的差分后的平稳序列来确定p和q。这里先明确两个概念。
a.自相关函数（ACF）：这是时间序列和它自身滞后版本之间的相关性的测试。比如在自相关函数可以比较时间的瞬间‘t1’…’t2’以及序列的瞬间‘t1-5’…’t2-5’ (t1-5和t2 是结束点)。
b.部分自相关函数(PACF):这是时间序列和它自身滞后版本之间的相关性测试，但是是在预测（已经通过比较干预得到解释）的变量后。如：滞后值为5，它将检查相关性，但是会删除从滞后值1到4得到的结果。

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

1
2
3
4
5
6

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至7阶（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法则。我们知道：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。不仅仅包括AIC准则，目前选择模型常用如下准则：

• AIC=-2 ln(L) + 2 k 中文名字：赤池信息量 akaike information criterion
• BIC=-2 ln(L) + ln(n)*k 中文名字：贝叶斯信息量 bayesian information criterion
• HQ=-2 ln(L) + ln(ln(n))*k hannan-quinn criterion
  构造这些统计量所遵循的统计思想是一致的，就是在考虑拟合残差的同时，依自变量个数施加“惩罚”。但要注意的是，这些准则不能说明某一个模型的精确度，也即是说，对于三个模型Ａ，Ｂ，Ｃ，我们能够判断出Ｃ模型是最好的，但不能保证Ｃ模型能够很好地刻画数据，因为有可能三个模型都是糟糕的。

  1 2 3 4 5 6 7 8

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

  
  可以看到ARMA(7,0)的aic，bic，hqic均最小，因此是最佳模型。（todo:aic,bic,hqic含义）

如何通过ACF和PACF观察pq取值
自相关和偏自相关图作为时间序列判断阶数的重要方法，很多童鞋在刚接触的时候都会在如何判断拖尾截尾上有疑问。今天楼主就通过图片来告诉你怎么看！

AR模型：自相关系数拖尾，偏自相关系数截尾；
MA模型：自相关系数截尾，偏自相关函数拖尾；
ARMA模型：自相关函数和偏自相关函数均拖尾。
根据输出结果，自相关函数图拖尾，偏自相关函数图截尾，且n从2或3开始控制在置信区间之内，因而可判定为AR(2)模型或者AR(3)模型。

这张图可以看到，很明显的自相关和偏自相关都是拖尾，因为数据到后面还有增大的情况，没有明显的收敛趋势。

如果你的图片成这样，估计十有八九是一个ARMA模型了。自相关7阶拖尾（n从7开始缩至置信区间），偏自相关2阶拖尾。

07，模型检验

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

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

1
2
3
4
5

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)

07.2 做D-W检验

德宾-沃森（Durbin-Watson）检验。德宾-沃森检验,简称D-W检验，是目前检验自相关性最常用的方法，但它只使用于检验一阶自相关性。因为自相关系数ρ的值介于-1和1之间，所以 0≤DW≤４。并且DW＝O＝＞ρ＝１　　 即存在正自相关性
DW＝４＜＝＞ρ＝－１　即存在负自相关性
DW＝２＜＝＞ρ＝０　　即不存在（一阶）自相关性
因此，当DW值显著的接近于O或４时，则存在自相关性，而接近于２时，则不存在（一阶）自相关性。这样只要知道ＤＷ统计量的概率分布，在给定的显著水平下，根据临界值的位置就可以对原假设Ｈ０进行检验。

1

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

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

07.3 观察是否符合正态分布

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

1
2
3
4

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

07.4 Ljung-Box检验

Ljung-Box test是对randomness的检验,或者说是对时间序列是否存在滞后相关的一种统计检验。对于滞后相关的检验，我们常常采用的方法还包括计算ACF和PCAF并观察其图像，但是无论是ACF还是PACF都仅仅考虑是否存在某一特定滞后阶数的相关。LB检验则是基于一系列滞后阶数，判断序列总体的相关性或者说随机性是否存在。
时间序列中一个最基本的模型就是高斯白噪声序列。而对于ARIMA模型，其残差被假定为高斯白噪声序列，所以当我们用ARIMA模型去拟合数据时，拟合后我们要对残差的估计序列进行LB检验，判断其是否是高斯白噪声，如果不是，那么就说明ARIMA模型也许并不是一个适合样本的模型。

1
2
3
4

r,q,p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
data = np.c_[range(1,41), r[1:], q, p]
table = pd.DataFrame(data, columns=['lag', "AC", "Q", "Prob(>Q)"])
print(table.set_index('lag'))

检验的结果就是看最后一列前十二行的检验概率（一般观察滞后1~12阶），如果检验概率小于给定的显著性水平，比如0.05、0.10等就拒绝原假设，其原假设是相关系数为零。就结果来看，如果取显著性水平为0.05，那么相关系数与零没有显著差异，即为白噪声序列。

08，模型预测

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

1
2
3
4
5

predict_sunspots = arma_mod20.predict('2090', '2100', dynamic=True)
print(predict_sunspots)
fig, ax = plt.subplots(figsize=(12, 8))
ax = dta.ix['2001':].plot(ax=ax)
predict_sunspots.plot(ax=ax)

前面90个数据为测试数据，最后10个为预测数据；从图形来，预测结果较为合理。

参考资料

R语言时间序列分析：http://blog.csdn.net/jiabiao1602/article/details/43153139
ARIMA模型：http://wiki.mbalib.com/wiki/ARIMA%E6%A8%A1%E5%9E%8B
时间序列预测全攻略（附带Python代码）：http://www.36dsj.com/archives/44065
[python] 时间序列分析之ARIMA：http://blog.csdn.net/u010414589/article/details/49622625
（侧重数据原理）ARIMA模型：http://blog.csdn.net/aspirinvagrant/article/details/46323271
prophet:时间序列预测原理：http://blog.csdn.net/a358463121/article/details/70194279
时间序列完全教程（R）：http://blog.csdn.net/Earl211/article/details/50957029
图片告诉你怎么看拖尾截尾！：http://bbs.pinggu.org/thread-3241517-1-1.html
时间序列分析(一) 如何判断序列是否平稳：http://blog.csdn.net/bi_hu_man_wu/article/details/64918870
时间序列分析之ARIMA模型预测__R篇：http://www.cnblogs.com/bicoffee/p/3838049.html
时间序列分析—(ARIMA模型)：http://www.17bigdata.com/%E6%97%B6%E9%97%B4%E5%BA%8F%E5%88%97%E5%88%86%E6%9E%90-arima%E6%A8%A1%E5%9E%8B.html
R学习日记——时间序列分析之ARIMA模型预测：http://www.17bigdata.com/r%E5%AD%A6%E4%B9%A0%E6%97%A5%E8%AE%B0-%E6%97%B6%E9%97%B4%E5%BA%8F%E5%88%97%E5%88%86%E6%9E%90%E4%B9%8Barima%E6%A8%A1%E5%9E%8B%E9%A2%84%E6%B5%8B.html