导入数据,并转化为时间序列
#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()
平稳性检测
- 方法一:时序图
from pylab import *
plt.plot(ts)
plt.title(u'当天出库')
show()
输出:
序列始终在一个常数值附近随机波动,且波动范围有界,且没有明显的趋势性或周期性,所以可认为是平稳序列。
- 方法二:自相关图
from statsmodels.graphics.tsaplots import plot_pacf,plot_acf
plot_acf(ts)
show()
输出:
自相关系数会很快衰减向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
输出:
单位根检验统计量对应的P值远小于0.05,故该序列可确认为平稳序列。
纯随机性检验(白噪声检验)
from statsmodels.stats.diagnostic import acorr_ljungbox
print u'序列的纯随机性检测结果为:',acorr_ljungbox(ts,lags = 1)
输出:
序列的纯随机性检测结果为: (array([ 9.10802245]), array([ 0.00254491]))
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()
输出:
可以看出,模型的阶次应该为(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
输出:
(3, 2)
sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='bic')['bic_min_order'] # BIC
输出:
(1, 0)
sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='hqic')['hqic_min_order'] # HQIC
输出:
(3, 2)
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()给出一份模型报告
输出:
接下来预测最后10天的数据:
tempModel.forecast(10)
输出:
(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]]))
最后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
输出:
0.0353600467617
拟合效果远小于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()
效果图:
至此,整个流程结束。但是,拟合效果并不好。
总结
导致拟合效果欠佳的原因可能有:
- 使用数据为原始数据,未加任何预处理(主要原因)。原始数据中存在着异常值、不一致、缺失值,严重影响了建模的执行效率,造成较大偏差。;
- 在模型定阶过程中,为了控制计算量,限制AR最大阶不超过6,MA最大阶不超过4,从了影响了参数的确定,导致局部最优。
接下来会从这两个方面考虑,改进并完善结果。