当前位置:   article > 正文

时间序列实战(一)_sm.tsa.arma_order_select_ic

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,从了影响了参数的确定,导致局部最优。

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

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/2023面试高手/article/detail/91850
推荐阅读
相关标签
  

闽ICP备14008679号