【挖掘模型】： Python 主要时序模式算法—StatsMod

数据源：arima_data.xls

arima_data.xls

数据结果

Paste_Image.png

ARIMA模型包

    # 使用该模型需要一些列判别操作
        # 平稳性检测 -》 白噪声检验 ——》 是否差分 -》 AIC和BIC指标 --> 模型定阶 --》预测
        # acf()  计算自相关系数
            # autocorr = acf(data,unbiased = False,nlags = 40,qstat = False,fft = False ,alpha = None)
            # data_观测值序列(时间序列,dataFrame or Series)
            # autocorr_观测值序列自相关函数
            # 其他为可选参数，如 qstat = True,同事返回Q统计量和对应的p值
        # plot_acf()/plot_pacf()   画自相关系数图和偏相关系数图
            # p = plot_acf(data) 
            # 返回一个Matplotlib对象，可以用show()方法显示图像
        # adfuller()  对观测值序列进行单位根检验(ADF test)
            # h = adffuller(Series,maxlag = None,regression = 'c',autolag = 'AIC',store = False,regresults = False)
            # 输入参数Series为一维观测值序列，依次返回adf,pvalue,usedlag,nobs,critical,values,icbest,regresults,resstore
        # diff()  对观测值序列进行差分计算
            # D。diff()_D为dataframe 或者 Series
        # arima  设置时序模式的建模参数，创建ARIMA时序模型
            # arima = ARIMA(data,(p,d,q)).fit()
            # data参数为输入时间序列，d为差分次数
        # summary()/summary2()   生成已有模型报告
            # 包含模型系数，标准误差，p值，AIC和BIC值等
        # aic/bic/hqic  计算模型的AIC、BIC、HQIC值
            # arima.aic/arima.bic/arima.hqic
        # forecast()   预测
            # a,b,c = arima.forecase(num)
            # num为预测的天数,a_返回的预测值,b_预测的误差,c_预测值置信区间
        # acorr_ljungbox()  检测是否为白噪声序列
            # acorr_ljungbox(data,lags=1)
            # lags 为滞后数,返回统计量和p值

源代码

import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

#参数初始化
discfile = 'F:/python 数据挖掘分析实战/Data/arima_data.xls'
forecastnum = 5

#读取数据，指定日期列为指标，Pandas自动将“日期”列识别为Datetime格式
data = pd.read_excel(discfile, index_col = u'日期')

#时序图
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()
# 从该图说明自相关系数长期大于0，序列间有很强的长期相关性

#平稳性检测
from statsmodels.tsa.stattools import adfuller as ADF
print(u'原始序列的ADF检验结果为：', ADF(data[u'销量']))
#返回值依次为adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore
# P值显著大于0.05，该序列是非平稳的时间序列

#差分后的结果，对原始数据进行一阶拆分，并进行平稳性和白噪声检验
D_data = data.diff().dropna()
D_data.columns = [u'销量差分']
D_data.plot() #时序图


plt.show()
plot_acf(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，所以一阶差分之后的序列是平稳非白噪声序列

from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(D_data).show() #偏自相关图，显示出拖尾性，本来从自相关图和偏自相关图可以看出p,q...看不懂

"""
#第二步、确定p、q，---》 0,1
res = sm.tsa.arma_order_select_ic(
    D_data.dropna(), 
    max_ar=8, 
    max_ma=8, 
    ic=['aic', 'bic', 'hqic'], 
    trend='nc'
)
"""



from statsmodels.tsa.arima_model import ARIMA

data[u'销量'] = data[u'销量'].astype(float)

#定阶
pmax = int(len(D_data)/10) #一般阶数不超过length/10
qmax = int(len(D_data)/10) #一般阶数不超过length/10
bic_matrix = [] #bic矩阵
for p in range(pmax+1):
  tmp = []
  for q in range(qmax+1):
    try: #存在部分报错，所以用try来跳过报错。
      tmp.append(ARIMA(data, (p,1,q)).fit().bic)
    except:
      tmp.append(None)
  bic_matrix.append(tmp)

bic_matrix = pd.DataFrame(bic_matrix) #从中可以找出最小值

p,q = bic_matrix.stack().idxmin() #先用stack展平，然后用idxmin找出最小值位置。
print(u'BIC最小的p值和q值为：%s、%s' %(p,q)) 

model = ARIMA(data, (p,1,q)).fit() #建立ARIMA(0, 1, 1)模型
model.summary2() #给出一份模型报告
model.forecast(5) #作为期5天的预测，返回预测结果、标准误差、置信区间。

参考资料：《Python数据分析与挖掘实战》