时间序列分析之Python实现

为啥使用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()