时间序列实战（一）_sm.tsa.arma_order_select_ic

导入数据，并转化为时间序列

#coding:utf-8
import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pylab as plt
plt.rcParams['font.sans-serif']=['SimHei']
from matplotlib.pylab import rcParams
from statsmodels.tsa.stattools import adfuller

dateparse = lambda dates: pd.datetime.strptime(dates, '%Y/%m/%d')
data = pd.read_csv('mian.csv', parse_dates='date', index_col='date',date_parser=dateparse)

rcParams['figure.figsize'] = 10, 5
ts = data['out'] 
ts.tail()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

平稳性检测

• 方法一：时序图

from pylab import *
plt.plot(ts)
plt.title(u'当天出库')
show()
1
2
3
4

输出：

序列始终在一个常数值附近随机波动，且波动范围有界，且没有明显的趋势性或周期性，所以可认为是平稳序列。

• 方法二：自相关图

from statsmodels.graphics.tsaplots import plot_pacf,plot_acf
plot_acf(ts)
show()
1
2
3

输出：

自相关系数会很快衰减向0，所以可认为是平稳序列。

• 方法三：ADF单位根检验（精确判断）

temp = np.array(ts)
t = sm.tsa.stattools.adfuller(temp)  # ADF检验
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
output
1
2
3
4
5
6
7
8
9
10
11

输出：

单位根检验统计量对应的P值远小于0.05，故该序列可确认为平稳序列。

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

from statsmodels.stats.diagnostic import acorr_ljungbox
print u'序列的纯随机性检测结果为：',acorr_ljungbox(ts,lags = 1)
1
2

输出：

序列的纯随机性检测结果为： (array([ 9.10802245]), array([ 0.00254491]))
1

P=0.00254491，统计量的P值小于显著性水平0.05，则可以以95%的置信水平拒绝原假设，认为序列为非白噪声序列（否则，接受原假设，认为序列为纯随机序列。）

综上：原序列为平稳非白噪声序列，适用于ARMA模型。

识别ARMA模型阶次

• 方法一：ACF、PACF 判断模型阶次

from statsmodels.graphics.tsaplots import plot_pacf,plot_acf
plot_acf(ts)
plot_pacf(ts)
show()
1
2
3
4

输出：

可以看出，模型的阶次应该为(200,400)，阶数高，计算量过大。采用另外一种方法确定阶数。

• 方法二：信息准则定阶
  目前选择模型常用如下准则： （其中L为似然函数，k为参数数量，n为观察数）
  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
  我们常用的是AIC准则，同时需要尽量避免出现过拟合的情况。所以优先考虑的模型应是AIC值最小的那一个模型。

为了控制计算量，在此限制AR最大阶不超过6，MA最大阶不超过4。 但是这样带来的坏处是可能为局部最优。

import statsmodels.api as sm
sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='aic')['aic_min_order']  # AIC
1
2

输出：

(3, 2)
1

sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='bic')['bic_min_order']  # BIC
1

输出：

(1, 0)
1

sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='hqic')['hqic_min_order'] # HQIC
1

输出：

(3, 2)
1

AIC求解的模型阶次为(3,2)
BIC求解的模型阶次为（1,0）
HQIC求解的模型阶次为（3,2）
这里就以AIC准则为准，选择（3,2），也可依次尝试每一种准则，选择最优。

模型的建立及预测

上一步骤已确定了ARMA模型的阶数为（3,2），接下来进行模型的建立和预测工作。将原数据分为训练集和测试集，选择最后10个数据用于预测。

order = (3,2)
train = ts[:-10]
test = ts[-10:]
tempModel = sm.tsa.ARMA(train,order).fit()
#tempModel.summary2()给出一份模型报告
1
2
3
4
5

输出：

接下来预测最后10天的数据：

tempModel.forecast(10)
1

输出：

(array([  5.29077389,   3.45200299,   4.61117218,   6.3501017 ,
          8.20994055,   9.98513142,  11.51878938,  12.68732223,
         13.40622711,  13.6351102 ]),
 array([ 23.03022024,  23.10739399,  23.12377736,  23.15402121,
         23.17973782,  23.19614362,  23.20345806,  23.20486727,
         23.20499493,  23.20820176]),
 array([[-39.84762834, 50.42917612],
        [-41.837657  , 48.74166298],
        [-40.71059863, 49.93294299],
        [-39.03094596, 51.73114937],
        [-37.22151076, 53.64139185],
        [-35.47847465, 55.4487375 ],
        [-33.95915273, 56.99673149],
        [-32.79338188, 58.16802634],
        [-32.07472721, 58.88718143],
        [-31.85212939, 59.12234979]]))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

最后10天的预测数据为：
5.29077389, 3.45200299, 4.61117218, 6.3501017 ,8.20994055,
9.98513142, 11.51878938, 12.68732223,13.40622711, 13.6351102

拟合效果：

delta = tempModel.fittedvalues - train
score = 1 - delta.var()/train.var()
print score
1
2
3

输出：

0.0353600467617
1

拟合效果远小于1，可见效果不好。。。

predicts = tempModel.predict('2016/4/21', '2016/4/30', dynamic=True)
print len(predicts)
comp = pd.DataFrame()
comp['original'] = test
comp['predict'] = predicts
comp.plot()
1
2
3
4
5
6

效果图：

至此，整个流程结束。但是，拟合效果并不好。

总结

导致拟合效果欠佳的原因可能有：

• 使用数据为原始数据，未加任何预处理（主要原因）。原始数据中存在着异常值、不一致、缺失值，严重影响了建模的执行效率，造成较大偏差。；
• 在模型定阶过程中，为了控制计算量，限制AR最大阶不超过6，MA最大阶不超过4，从了影响了参数的确定，导致局部最优。

接下来会从这两个方面考虑，改进并完善结果。