python学习日志3--ARIMA时间序列模型预测_section_3_arima_modeling_tutorial

前言

这篇文章主要讲述如何使用python实现时间序列ARIMA预测算法

一、代码

代码如下（示例）：

#跟着视频学习的代码，记录一下。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import itertools

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

from matplotlib.pylab import style #自定义图表风格
style.use('ggplot')

# 解决中文乱码问题
plt.rcParams['font.sans-serif'] = ['Simhei']
# 解决坐标轴刻度负号乱码
plt.rcParams['axes.unicode_minus'] = False

#pip install statsmodels   导入python统计分析库
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf  #自相关图、偏自相关图
from statsmodels.tsa.stattools import adfuller as ADF #平稳性检验
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
import statsmodels.api as sm #D-W检验,一阶自相关检验
from statsmodels.graphics.api import qqplot #画QQ图,检验一组数据是否服从正态分布
from statsmodels.tsa.arima_model import ARIMA

#读取往年数据
sale=pd.read_excel('C:/Users/john/Desktop/arima_data.xls',index_col='日期')

print(sale.head())   #打印数据头五行
print(sale.tail())   #打印数据尾五行

print(sale.info())     #打印数据整体信息
print('-----')
sale.销量=sale.销量.astype('float')   #完成数据的类型转换
print(sale.info())

# sale.plot()                #将当前曲线绘制出来
# plt.show()

# plot_acf(sale,lags=35)  #绘制自相关函数 并且停留五秒钟  这是第一种平稳性检验的方法
# #plt.pause(5)

print('原始序列的ADF检验结果为：',ADF(sale.销量))  #第二项大于0.05非平稳(接受原假设)，小于0.05平稳
d1_sale=sale.diff(periods=1, axis=0).dropna()   #计算一阶差分  看是否平稳

#时序图
# plt.figure(figsize=(10,5))   #进行一次差分后绘图
# d1_sale.plot()
# plt.show()
#解读：在均值附件比较平稳波动

#自相关图
#plot_acf(d1_sale,lags=34).show()   #acf自相关图
#plt.pause(10)
#解读：有短期相关性，但趋向于零。

#平稳性检验
# print('一阶差分的ADF检验结果为：',ADF(d1_sale.销量))   #小于0.05   平稳性得到确定
#
# #解读：P值小于显著性水平α（0.05），拒绝原假设（非平稳序列），说明一阶差分序列是平稳序列。
# #
# print('一阶差分序列的白噪声检验结果为：',acorr_ljungbox(d1_sale,lags=1))#返回统计量、P值

#解读：p值小于0.05，拒绝原假设（纯随机序列），说明一阶差分序列是非白噪声。
#第二个值:p的值小于0.05  说明不是白噪声  需要大于0.05

#下面需要确定参数  p  q  d
d1_sale=sale.diff(periods=1, axis=0).dropna()
# # print('1阶序列的ADF检验结果为：',ADF(d1_sale.销量))  #第二项大于0.05非平稳(接受原假设)，小于0.05平稳
# # #自相关图
# plot_acf(d1_sale,lags=15).show()
#
# #解读：有短期相关性，但趋向于零。
#
# #偏自相关图
# plot_pacf(d1_sale,lags=15).show()

# plt.pause(5)
#偏自相关图
# plot_pacf(d1_sale,lags=15).show()
# plt.pause(15)
#解读：自相关图，1阶截尾；偏自相关图，拖尾。则ARIMA(p,d,q)=ARIMA(0,1,1)

# pmax=4 #一般阶数不超过length/10
# qmax=4#一般阶数不超过length/10
# print('pmax',pmax)
# print('qmax',qmax)

# 当多组值都不符合时，遍历多组值，得出最好的值
# bic_matrix=[]
# for p in range(pmax+1):
#     tmp=[]
#     for q in range(qmax+1):
#         try:
#             # tmp.append(ARIMA(sale,(p,1,q)).fit().bic)
#             tmp.append(sm.tsa.ARIMAX(d1_sale, order=(p, 1, q)).fit().bic)
#             print(sm.tsa.ARIMAX(d1_sale, order=(p, 1, q)).fit().bic)
#         except:
#             tmp.append(None)
#     bic_matrix.append(tmp)
#
# bic_matrix=pd.DataFrame(bic_matrix)
# print(bic_matrix)


# 当多组值都不符合时，遍历多组值，得出最好的值
p_min = 0
d_min = 0
q_min = 0
p_max = 6
d_max = 0
q_max = 4


# results_bic = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min, p_max + 1)],
#                            columns=['MA{}'.format(i) for i in range(q_min, q_max + 1)])
#
# for p, d, q in itertools.product(range(p_min, p_max + 1),
#                                  range(d_min, d_max + 1),
#                                  range(q_min, q_max + 1)):
#     if p == 0 and d == 0 and q == 0:
#         # results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
#         continue
#
#     try:
#
#         model = sm.tsa.ARIMAX(d1_sale, order=(p, d, q),
#                                # enforce_stationarity=False,
#                                # enforce_invertibility=False,
#                                )
#         results = model.fit()
#         # results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.bic
#     except:
#         continue
# results_bic = results_bic[results_bic.columns].astype(float)


# fig, ax = plt.subplots(figsize=(10, 8))
# ax = sns.heatmap(results_bic,
#                  mask=results_bic.isnull(),
#                  ax=ax,
#                  annot=True,
#                  fmt='.2f',
#                  );
# ax.set_title('BIC');
# plt.show()
#
train_results = sm.tsa.arma_order_select_ic(d1_sale, ic=['aic', 'bic'], max_ar=6, max_ma=4)

aic_results=sm.tsa.arma_order_select_ic(d1_sale,max_ar=6,max_ma=4,ic='aic')['aic_min_order']
bic_results=sm.tsa.arma_order_select_ic(d1_sale,max_ar=6,max_ma=4,ic='bic')['bic_min_order']
print(aic_results)
print(bic_results)
# print('AIC', train_results.aic_min_order)
# print('BIC', train_results.bic_min_order)
print('-----')
p=bic_results[0];
q=bic_results[1];
print(p)
print(q)
# model=sm.tsa.arima.ARIMA(sale, order=(p, 1, q)).fit()
model=sm.tsa.arima.ARIMA(sale, order=(p, 1, q)).fit()
print(model.summary())
#
resid=model.resid

#自相关图  判断残差
# plot_acf(resid,lags=35).show()
# #
# #解读：有短期相关性，但趋向于零。
#
# # #偏自相关图
# plot_pacf(resid,lags=17).show()
# #
# # #偏自相关图
# plot_pacf(resid,lags=17).show()
# plt.pause(50)
#
# qqplot(resid, line='q', fit=True).show()
# plt.pause(100)


# qqplot(resid, line='q', fit=True).show()
# 并且DW＝O ＜ ＝＞ ρ＝１　　即存在正自相关性
# DW＝４ ＜＝＞ ρ＝－１　即存在负自相关性
# DW＝２ ＜＝＞ ρ＝０　　即不存在（一阶）自相关性
# plt.pause(10)
print('D-W检验的结果为：',sm.stats.durbin_watson(resid.values))
# 方法一
print('残差序列的白噪声检验结果为：',acorr_ljungbox(resid,lags=1))#返回统计量、P值

#解读：残差是白噪声
#预测
print('未来7天的销量预测：')
model.forecast(7) #预测、标准差、置信区间
forecast=pd.Series(model.forecast(7)[0],index=pd.date_range('2015-2-7',periods=7,freq='D'))
print(forecast)
data=pd.concat((sale,forecast),axis=0)
data.columns=['销量','未来7天销量']
plt.figure(figsize=(10,5))
data.plot()
plt.show()
plt.pause(100)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205

总结

代码可以直接运行，如有不懂请留言哦，