以ARIMA模型为例介绍时间序列算法在python中是如何实现的，一下是应用Python语言建模步骤：

–– coding: utf-8 ––

“””

Created on Mon Apr 2 16:45:36 2018

@author: houy

“”“

arima模型对时间序列的要求是平稳型

import pandas as pd

参数初始化7ku9

discfile = ‘arima_data.xlsx’

读取数据，制定日期列为指标，Pandas自动将“日期”列识别为Datatime格式

data = pd.read_excel(discfile,index_col = 0)

print(data.head())

print(‘\n Data Types:’)

print(data.dtypes)

时序图

import matplotlib.pyplot as plt

plt.rcParams[‘font.sans-serif’] = [‘SimHei’] #用来正常显示中文标签

plt.rcParams[‘axes.unicode_minus’] = False #用来正常显示负号

data.plot()

plt.show()

自相关图

from statsmodels.graphics.tsaplots import plot_acf

plot_acf(data).show() #显示出很强的自相关性

平稳性检测

from statsmodels.tsa.stattools import adfuller as ADF

print(u’原始序列的ADF检验结果为：’,ADF(data[u’销量’]))

返回值依次为：adf,pvalue,usedlag,nobs,critical value,icbest,regresults,resstore

返回结果中，pvalue的值显著大于0.05.该序列为非平稳序列

##########通过时间序列的差分将非平稳的时间序列变为平稳时间序列，d次差分后平稳，使用ARIMA(p,d,q)模型

D_data = data.diff().dropna()

D_data.columns= [u’销量差分’]

D_data.plot() #时序图

plt.show()

plot_acf(D_data).show() #自相关图 差分后的序列迅速落入区间内，并向0靠拢，序列没有自相关性

对差分后的序列做偏自相关检验

from statsmodels.graphics.tsaplots import plot_pacf

plot_pacf(D_data).show() #偏自相关图 差分后的序列没有显示出骗子相关性

print(u’差分序列的ADF检验结果为：’,ADF(D_data[u’销量差分’])) #p值小于0.05

对差分后的序列做白噪声检验

from statsmodels.stats.diagnostic import acorr_ljungbox

print(u’差分序列的白噪声检验结果为：’,acorr_ljungbox(D_data,lags=1)) #返回统计量和p值

通过检验，p值小于0.05，该序列为白噪声序列

##########比较一阶差分后的序列和二阶差分后的序列

一阶差分

fig = plt.figure(figsize=(12,8))

ax1 = fig.add_subplot(111)

diff1 = data.diff(1)

diff1.plot(ax = ax1)

一阶差分的序列的均值和方差已基本平稳

二阶差分

fig = plt.figure(figsize = (12,8))

ax2 = fig.add_subplot(112)

diff2 = data.diff(2)

diff2.plot(ax = ax2)

可以看出二阶差分后的时间序列与一阶差分相差不大，并且二者随时间推移，序列的均值和方差保持不变，可设置差分次数d为1

##########已找到稳定的时间序列，选择合适的ARIMA模型，即确定p,q

合适的pq,第一步检验平稳序列的自相关图和偏自相关图

import statsmodels.api as sm

dta = data.diff(1)[1:]

fig = plt.figure(figsize = (12,8))

ax1 = fig.add_subplot(211)

fig1 = sm.graphics.tsa.plot_acf(dta[u’销量’],lags=10,ax=ax1)

ax2 = fig.add_subplot(212)

fig2 = sm.graphics.tsa.plot_pacf(dta[u’销量’],lags=10,ax=ax2)

从两图观察到：可有ARIMA(0,2)、ARIMA(1,0)、ARIMA(0,1)模型作为选择

q =2 的移动平均模型， p=1的自回归模型，q = 1的自回归模型

”’

使用AIC准则进行模型选择，有限考虑AIC值小的那个模型

(其中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

”’

对三个模型分别做三个统计量检验

模型

arma_mod20 = sm.tsa.ARMA(data,(2,0)).fit()

print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)

arma_mod01 = sm.tsa.ARMA(data,(0,1)).fit()

print(arma_mod01.aic,arma_mod01.bic,arma_mod01.hqic)

arma_mod10 = sm.tsa.ARMA(data,(1,0)).fit()

print(arma_mod10.aic,arma_mod10.bic,arma_mod10.hqic)

结果中，ARMA(0，1)的三个值均最小，是最佳模型

###########模型检验

残差QQ图

from statsmodels.graphics.api import qqplot

resid = arma_mod01.resid

fig = plt.figure(figsize=(12,8))

ax = fig.add_subplot(111)

fig = qqplot(resid,line =’q’,ax=ax,fit = True)

残差自相关检验

fig = plt.figure(figsize=(12,8))

ax1 = fig.add_subplot(211)

fig = sm.graphics.tsa.plot_pacf(arma_mod01.resid.values.squeeze(),lags = 10,ax = ax1)

ax2 = fig.add_subplot(212)

fig = sm.graphics.tsa.plot_pacf(arma_mod01.resid,lags = 10,ax = ax2)

德宾-沃森D-W检验

”’

它只使用于检验一阶自相关性。因为自相关系数ρ的值介于-1和1之间，

所以 0≤DW≤４。并且DW＝O＝＞ρ＝１　　 即存在正自相关性

DW＝４＜＝＞ρ＝－１　即存在负自相关性

DW＝２＜＝＞ρ＝０　　即不存在(一阶)自相关性

因此，当DW值显著的接近于O或４时，则存在自相关性，而接近于２时，则不存在(一阶)自相关性。

”’

print(sm.stats.durbin_watson(arma_mod01.resid.values))

检验结果为1.95，说明不存在自相关性

Ljung-Box检验

import numpy as np

”’

输出有点问题，还未解决

r,q,p = sm.tsa.acf(resid.values.squeeze(), qstat=True)

datap = np.c_[range(1,36), r[1:], q, p]

table = pd.DataFrame(datap, columns=[‘lag’, “AC”, “Q”, “Prob(>Q)”])

print(table.set_index(‘lag’))

”’

###########模型预测

模型确定后，对未来9日的数据进行预测

predict_sunspots = arma_mod01.predict(‘2015-2-07’,’2015-2-15’,dynamic = True)

fig, ax = plt.subplots(figsize = (12,8))

print(predict_sunspots)

predict_sunspots[0] += data[‘2015-02-06’:][u’销量’]

data = pd.DataFrame(data)

for i in range(len(predict_sunspots)-1):

predict_sunspots[i+1] = predict_sunspots[i] + predict_sunspots[i+1]

print(predict_sunspots)

ax = data.ix[‘2015’:].plot(ax=ax)

predict_sunspots.plot(ax=ax)

plt.show()