为啥使用Python?,虽然经典的Python机器学习库scikit-learn中并没有内置时间序列分析的API,但是Python中也有时间序列模型,主要是统计分析库statsModels中的tsa模块,内置了ARMA,ARIMA等算法模型;另外,在使用模型建模的时候,需要进行的一系列判别操作及检验也都提供了方法和函数,比如:平稳性检验,白噪声检验,是否差分(pandas库带有差分方法),AIC,BIC模型定阶等等。
导入模块
import pandas as pd import numpy as np import matplotlib.pyplot as plt import statsmodels.api as smfrom statsmodels.graphics.api import qqplot#读取数据,数据是2015/1/1----2015/2/6某个餐厅的销售数据data = pd.read_excel('timelist.xlsx',index_col='日期')#数据时序图data.plot()plt.show()
平稳性检测:
从图上看时间序列是不平稳的,我们可以通过adf检验来判断是否平稳,也可使用Kpss检验
from statsmodels.tsa.stattools import adfuller,kpssadf = adfuller(data['销量']) #print(dftest)adfout = pd.Series(adf[0:4],index=['Test statistic','p-value','lag Users',\ 'Number of Observations Used'])print('ADF检测:')print(adfout)kpsstest=kpss(data['销量'])print('KPSS检测:',kpsstest)
结果如下:adf检测的p值远远大于0.05,而KPSS检测的p值0.055与0.05相差不大,可以判定这个时间序列是不平稳的。ADF检测:
Test statistic 1.813771
p-value 0.998376
lag Users 10.000000
KPSS检测: (0.4503223095873123, 0.05546452172960678, 10, {'10%': 0.347, '5%': 0.463, '2.5%': 0.574, '1%': 0.739})
差分处理
针对不平稳的数据进行处理,差分处理
diff1=data.diff(1) #一阶差分处理diff1=diff1.dropna()diff1.columns=['差分销量']
差分后在进行平稳性检测
print('差分后再ADF检测平稳性')adf2=adfuller(diff1['差分销量'])print(adf2)
检测结果如下:p值为0.022,小于0.05,可以判定,差分后的序列已经平稳(-3.1560562366723537, 0.022673435440048798, 0, 35, {'1%': -3.6327426647230316, '5%': -2.9485102040816327, '10%': -2.6130173469387756}, 287.5909090780334)
白噪声检测,返回结果:[统计量,pvalue]
from statsmodels.stats.diagnostic import acorr_ljungboxprint('白噪声检测结果:',acorr_ljungbox(diff1,lags=1))
白噪声检测结果,p值0.00077远远小于0.05,所以差分后的序列是平稳非白噪声序列白噪声检测结果: (array([11.30402222]), array([0.00077339]))
模型定阶
1,通过序列的自相关图和偏自相关图的性质进行定阶
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf#一阶差分的时间序列的均值和方差已经基本平稳,因此可以将差分次数d设置为1。 #自相关图和偏相关图plot_acf(diff1).show()plot_pacf(diff1).show()
一阶差分后的序列自相关图
一阶差分后序列的偏自相关图
2,通过相对最优模型识别,取bic信息量最小的模型阶数
from statsmodels.tsa.arima_model import ARIMA#定阶,一般p,q最大值都不会超过数据长度的十分之一pmax=int(len(data)/10)qmax=int(len(data)/10)bic_matrix=[]for p in range(pmax+1): tmp=[] for q in range(qmax+1): try: tmp.append(ARIMA(data,(p,1,q)).fit().bic) except: tmp.append(None)bic_matrix.append(tmp)bic_data = pd.DataFrame(bic_matrix)print(bic_data)p,q = bic_data.stack().idxmin()print('BIC最小的p和q值为:%s,%s'%(p,q))
计算的BIC数据如下:432.068472 422.510082 426.088911 426.595507
423.628276 426.073601 NaN NaN
426.774824 427.395847 NaN NaN
430.317524 NaN NaN 436.478109
BIC最小的p和q值为:0,1
创建ARIMA(0,1,1)模型
model = ARIMA(data,(0,1,1)).fit()
模型检验
需要对模型进行残差检测是否为白噪声序列,满足才可以接着进行预测
resids = model.resid #残差,#要求残差为纯随机序列,也就是白噪声序列print('残差的白噪声检测结果:',acorr_ljungbox(resids,lags=1))#残差独立性检测,durbin_watson 在2附近说明满足残差独立性print('durbin_watson检测结果:',sm.stats.durbin_watson(resids))# qq图,检测是否符合正态分布qqplot(resids)
检验结果如下:白噪声检测p值0.95,残差是白噪声,durbin_watson检测结果1.973,在2附近,说明满足残差独立,qq图如下,符合正态分布残差的白噪声检测结果: (array([0.003904]), array([0.95017903]))
durbin_watson检测结果: 1.973487626529122
预测,现在我们可以使用训练的模型进行预测了
print(model.forecast(5)) #预测5天的,返回结果:预测结果,标准误差,置信区间(array([4873.96625288, 4923.92173955, 4973.87722621, 5023.83271288,
5073.78819955]),
array([ 73.08574135, 142.32683622, 187.54287785, 223.8028904 ,
254.95712673]),
array([[4730.72083205, 5017.2116737 ],[4644.96626651, 5202.87721258],
[4606.29994008, 5341.45451235],[4585.18710806, 5462.4783177 ],
[4574.08141355, 5573.49498555]]))
pred = model.predict('2015/2/7','2015/2/15',dynamic=True,typ='levels') #预测的时间段2015-02-07 4873.966253
2015-02-08 4923.921740
2015-02-09 4973.877226
2015-02-10 5023.832713
2015-02-11 5073.788200
2015-02-12 5123.743686
2015-02-13 5173.699173
2015-02-14 5223.654660
2015-02-15 5273.610146
Freq: D, dtype: float64
预测趋势图
fig, ax = plt.subplots()ax = data.loc['2015/1/20':].plot(ax=ax)fig = model.plot_predict('2015/2/7', '2015/2/15', dynamic=True, ax=ax,plot_insample=False)plt.show()