前言：我用的是python3.6版本。这篇博客参考了以下两篇博文，根据自己的理解重新整合，稍作修改代码，希望能够帮到大家。如果在过程中发现错误或者有更好的方式实现，欢迎大家发表评论。
参考1: https://blog.csdn.net/u010414589/article/details/49622625
参考2: https://blog.csdn.net/hal_sakai/article/details/51965657
正文：
一、时间序列建模基本步骤
1、获取被观测系统时间序列数据；
2、对数据绘图，观测是否为平稳时间序列；对于非平稳时间序列要先进行d阶差分运算，化为平稳时间序列；
3、经过第二步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF ，通过对自相关图和偏自相关图的分析，得到最佳的阶层 p 和阶数 q
4、由以上得到的d、q、p得到ARIMA模型。然后开始对得到的模型进行模型检验。
二、实战解析
基础库： pandas；numpy；scipy；matplotlib；statsmodels ：
#-*- coding:utf-8 -*- 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 plt.rcParams['font.sans-serif']=['SimHei'] #用来画图时正常显示中文标签 plt.rcParams['axes.unicode_minus']=False #用来画图时正常显示坐标轴负数 from statsmodels.graphics.tsaplots import plot_acf from statsmodels.graphics.tsaplots import plot_pacf
2.1获取数据
这里我们使用一个具有周期性的测试数据，进行分析。
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=np.array(dta,dtype=np.float) #转换数据类型 dta=pd.Series(dta) #print(type(dta)) dta.index=pd.Index(sm.tsa.datetools.dates_from_range('2001','2090')) plt.plot(dta,label='原始数据') plt.legend() #让label生效
#plt.show() #PyCharm IDE要输入这个命令才能显示图片，由于下面还要画图，
因此我只要在程序最后附上这句话即可
2.2 差分d
ARIMA 模型对时间序列的要求是平稳型。因此，当你得到一个非平稳的时间序列时，首先要做的即是做时间序列的差分，直到得到一个平稳时间序列。如果你对时间序列做d次差分才能得到一个平稳序列，那么可以使用ARIMA(p,d,q)模型，其中d是差分次数。plt.figure() #新开一个图窗 diff1=dta.diff(1) plt.plot(diff1,label='一阶') diff2=dta.diff(2) plt.plot(diff2,label='二阶') plt.legend()
一阶差分的时间序列的均值和方差已经基本平稳，不过我们还是可以比较一下二阶差分的效果。可以看出二阶差分后的时间序列与一阶差分相差不大，并且二者随着时间推移，时间序列的均值和方差保持不变。因此可以将差分次数d设置为1。
其实还有针对平稳的检验，叫“ADF单位根平稳型检验”(还没有学习这部分，等我学习了会更这部分内容的)。
2.3 合适的p,q
现在我们已经得到一个平稳的时间序列，接来下就是选择合适的ARIMA模型，即ARIMA模型中合适的p,qp,q。
2.3.1 第一步我们要先检查平稳时间序列的自相关图和偏自相关图。
其中lags
表示滞后的阶数，以下分别得到差分后的平稳序列的acf 图和pacf 图。
因为差分后的序列第一个是NAN，因此要去掉它，直接对diff1运行图形会运行不出来，有两种方法，#后是第二种
fig=plt.figure() ax1=fig.add_subplot(211) # 画子图 plot_acf(diff1[1:],lags=40,ax=ax1) #diff11=np.diff(dta) #plot_acf(diff11,lags=40,ax=ax1) ax2=fig.add_subplot(212) plot_pacf(diff1[1:],lags=40,ax=ax2)
通过两图观察得到：
* 自相关图显示滞后有三个阶超出了置信边界；
* 偏相关图显示在滞后1至7阶（lags 1,2,…，7）时的偏自相关系数超出了置信边界，从lag 7之后偏自相关系数值缩小至0
2.3.2 模型选择（p,q选择）
这部分文字来自参考的第一篇博文，根据acf和pacf图得到以下模型，但是我并不理解这部分，我不明白这些模型是如何通过观察得到的，如果有知道的童鞋希望可以告知。
根据上图，猜测有以下模型可以供选择：
1. ARMA(0,1)模型：即自相关图在滞后1阶之后缩小为0，且偏自相关缩小至0，则是一个阶数q=1的移动平均模型；
2. ARMA(7,0)模型：即偏自相关图在滞后7阶之后缩小为0，且自相关缩小至0，则是一个阶层p=7的自回归模型；
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
构造这些统计量所遵循的统计思想是一致的，就是在考虑拟合残差的同时，依自变量个数施加“惩罚”。但要注意的是，这些准则不能说明某一个模型的精确度，也即是说，对于三个模型Ａ，Ｂ，Ｃ，我们能够判断出Ｃ模型是最好的，但不能保证Ｃ模型能够很好地刻画数据，因为有可能三个模型都是糟糕的。
arma_mod70 = sm.tsa.ARMA(dta,(7,0)).fit() arma_mod01 = sm.tsa.ARMA(dta,(0,1)).fit() arma_mod71 = sm.tsa.ARMA(dta,(7,1)).fit() arma_mod80 = sm.tsa.ARMA(dta,(8,0)).fit() print(arma_mod70.aic,arma_mod70.bic,arma_mod70.hqic) print(arma_mod01.aic,arma_mod01.bic,arma_mod01.hqic) print(arma_mod71.aic,arma_mod71.bic,arma_mod71.hqic) print(arma_mod80.aic,arma_mod80.bic,arma_mod80.hqic)1619.1917731179842 1641.69006015 1628.26440493
1657.2172636424214 1664.71669265 1660.24147424
1605.686576696957 1630.6846734 1615.7672787
1597.935980998883 1622.9340777 1608.01668301
可以看到ARMA(8,0)的aic，bic，hqic均最小，因此是最佳模型。
2.4 模型检验
在指数平滑模型下，观察ARIMA模型的残差是否是平均值为0且方差为常数的正态分布（服从零均值、方差不变的正态分布），同时也要观察连续残差是否（自）相关。
2.4.1 对ARMA(8,0)模型所产生的残差做自相关图
残差的ACF和PACF图，可以看到序列残差基本为白噪声。
2.4.2 做D-W检验
德宾-沃森（Durbin-Watson）检验。德宾-沃森检验,简称D-W检验，是目前检验自相关性最常用的方法，但它只使用于检验一阶自相关性。因为自相关系数ρ的值介于-1和1之间，所以
0≤DW≤４。并且DW＝O＝＞ρ＝１ 即存在正自相关性
DW＝４＜＝＞ρ＝－１ 即存在负自相关性
DW＝２＜＝＞ρ＝０ 即不存在（一阶）自相关性
因此，当DW值显著的接近于0或４时，则存在自相关性，而接近于２时，则不存在（一阶）自相关性。这样只要知道ＤＷ统计量的概率分布，在给定的显著水平下，根据临界值的位置就可以对原假设Ｈ０进行检验。
print(sm.stats.durbin_watson(arma_mod80.resid.values))
结果是 2.02323317584，说明不存在自相关性。
2.4.3 是否符合正态分布
这里使用QQ图，它用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布。（QQ图还未学习，待更）
resid = arma_mod80.resid#残差 qqplot(resid, line='q',fit=True)
2.4.4 Ljung-Box 检验
Ljung-Box
test是对randomness的检验,或者说是对时间序列是否存在滞后相关的一种统计检验。对于滞后相关的检验，我们常常采用的方法还包括计算ACF和PCAF并观察其图像，但是无论是ACF还是PACF都仅仅考虑是否存在某一特定滞后阶数的相关。LB检验则是基于一系列滞后阶数，判断序列总体的相关性或者说随机性是否存在。
时间序列中一个最基本的模型就是高斯白噪声序列。而对于ARIMA模型，其残差被假定为高斯白噪声序列，所以当我们用ARIMA模型去拟合数据时，拟合后我们要对残差的估计序列进行LB检验，判断其是否是高斯白噪声，如果不是，那么就说明ARIMA模型也许并不是一个适合样本的模型。
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'))结果：
AC Q Prob(>Q)
lag
1.0 -0.014699 0.020101 0.887255
2.0 -0.045553 0.215348 0.897920
3.0 0.101611 1.197985 0.753488
4.0 0.063415 1.585168 0.811455
5.0 0.176177 4.608680 0.465475
6.0 0.004990 4.611135 0.594563
7.0 -0.209066 8.971443 0.254713
8.0 0.115697 10.323082 0.243078
9.0 0.020907 10.367765 0.321541
10.0 -0.220087 15.381123 0.118772
11.0 -0.050643 15.649936 0.154632
12.0 -0.031516 15.755374 0.202699
13.0 -0.055326 16.084532 0.244595
14.0 0.195279 20.239100 0.122782
15.0 -0.204407 24.851871 0.051968
16.0 0.034586 24.985716 0.070078
17.0 0.181442 28.719772 0.037199
18.0 -0.043290 28.935286 0.049176
19.0 -0.045298 29.174578 0.063288
20.0 0.044176 29.405417 0.080078
21.0 0.141981 31.824456 0.060989
22.0 0.118182 33.525145 0.054859
23.0 0.004506 33.527654 0.072303
24.0 -0.133570 35.765894 0.057824
25.0 0.061832 36.252912 0.067849
26.0 -0.021823 36.314524 0.086035
27.0 0.047269 36.608186 0.102616
28.0 0.131266 38.909330 0.082392
29.0 -0.080900 39.797706 0.087299
30.0 -0.026185 39.892329 0.106967
31.0 0.011452 39.910736 0.131066
32.0 -0.015213 39.943774 0.157970
33.0 -0.042655 40.208071 0.181254
34.0 0.006051 40.213485 0.214256
35.0 -0.015970 40.251881 0.248995
36.0 0.000407 40.251907 0.287553
37.0 -0.083808 41.349205 0.286411
38.0 -0.091955 42.695609 0.276349
39.0 0.002442 42.696577 0.315252
40.0 -0.071594 43.545399 0.322996
检验的结果就是看最后一列前十二行的检验概率（一般观察滞后1~12阶），如果检验概率小于给定的显著性水平，比如0.05、0.10等就拒绝原假设，其原假设是相关系数为零。就结果来看，如果取显著性水平为0.05，那么相关系数与零没有显著差异，即为白噪声序列。即：prob值均大于0.05，所以残差序列不存在自相关性。
2.5 预测结果
模型确定之后，就可以开始进行预测了，我们对未来十年的数据进行预测。
predict_sunspots = arma_mod80.predict('2090', '2100', dynamic=True) print(predict_sunspots) fig, ax = plt.subplots() ax = dta.ix['2001':].plot(ax=ax) predict_sunspots.plot(ax=ax)
plt.show()
2090-12-31 9542.423714
2091-12-31
12907.339818
2092-12-31 13980.500757
2093-12-31 14500.625514
2094-12-31 13893.455315
2095-12-31 13248.652060
2096-12-31 10960.214775
2097-12-31 10072.055950
2098-12-31 12682.506171
2099-12-31 13475.109548
2100-12-31 13614.243276
Freq: A-DEC, dtype: float64
前面90个数据为测试数据，最后10个为预测数据。（还没有预测模型好坏的标准）
三、
在模型选择方面我用的是自动选择BIC值最小的值：
from statsmodels.tsa.arima_model import ARIMA
#定阶:导入ARIMA pmax=int(len(diff1)/10) #一般阶数不超过length/10 qmax=int(len(diff1)/10) bic_matrix=[] #BIC for p in range(pmax+1): tmp=[] for q in range(qmax+1): try: #处理报错 tmp.append(ARIMA(dta,(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("BIC最小的p值和q值为：%s、%s"%(p,q)) model=ARIMA(dta,(p,1,q)).fit()
predict_sunspots = model.predict('2090', '2100', dynamic=True) print(predict_sunspots) fig, ax = plt.subplots() ax = dta.ix['2001':].plot(ax=ax) predict_sunspots.plot(ax=ax) plt.show()
结果：BIC最小的p值和q值为：6、6
2090-12-31
-867.898294
2091-12-31 4643.628545
2092-12-31 775.860782
2093-12-31 357.981031
2094-12-31 -1137.189339
2095-12-31 -406.305865
2096-12-31 -2943.157268
2097-12-31 -654.484380
2098-12-31 4519.613997
2099-12-31 727.625300
2100-12-31 119.724101
Freq: A-DEC, dtype: float64
问题：这样做出的模型 （6,6）效果还不如（8,0）的效果好，而且预测结果还出现了负数，希望有人能够解答这个问题。