ARIMA

# -*- coding: utf-8 -*-
# 用 ARIMA 进行时间序列预测

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.arima.model import ARIMA
#我用的新版本没有ARMA模块了
from statsmodels.graphics.tsaplots import acf,pacf,plot_acf,plot_pacf
from statsmodels.graphics.api import qqplot

# 1.创建数据
data = [5922, 5308, 5546, 5975, 2704, 1767, 4111, 5542, 4726, 5866, 6183, 3199, 1471, 1325, 6618, 6644, 5337, 7064, 2912, 1456, 4705, 4579, 4990, 4331, 4481, 1813, 1258, 4383, 5451, 5169, 5362, 6259, 3743, 2268, 5397, 5821, 6115, 6631, 6474, 4134, 2728, 5753, 7130, 7860, 6991, 7499, 5301, 2808, 6755, 6658, 7644, 6472, 8680, 6366, 5252, 8223, 8181, 10548, 11823, 14640, 9873, 6613, 14415, 13204, 14982, 9690, 10693, 8276, 4519, 7865, 8137, 10022, 7646, 8749, 5246, 4736, 9705, 7501, 9587, 10078, 9732, 6986, 4385, 8451, 9815, 10894, 10287, 9666, 6072, 5418]

data = pd.Series(data)#seies的作用是生成索引
'''例如my_list=[1,2,3]
series=pd.series(my_list)
print(series)结果是
0 1  2  
1 2  3'''

data.index = pd.Index(sm.tsa.datetools.dates_from_range('1901','1990'))

'''这段代码的功能是将数据集的索引值设置为时间序列索引。在这个示例中，
使用了statsmodels（sm）中的tsa（时间序列分析）模块的datetools对象来根据时间范围生成日期序列。
具体来说，它使用了`dates_from_range()`方法来生成从"1901"到"1990"之间的日期序列，然后使用`pd.Index()`方法将其转换为Pandas中的索引对象（Index），再将其赋值给数据集的索引。
时间序列索引的好处是，可以方便地对时间序列数据进行基于时间的切片和索引，
例如窗口统计和滚动平均。此外，它还可以与许多时间序列分析工具（例如ARIMA模型）进行良好的兼性。'''

data.plot(figsize=(12,8))#padas中的plot方法绘制长12，宽8的图形。
#绘制时序的数据图
plt.show()

#2.下面我们先对非平稳时间序列进行时间序列的差分，找出适合的差分次数d的值：
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(111)
'''创建子图'''
diff1 = data.diff(1)
'''pandas 中的diff()函数是指做差分，(1)表示一阶差分，(2)表示2阶差分'''
diff1.plot(ax=ax1)
#这里是做了1阶差分，可以看出时间序列的均值和方差基本平稳，不过还是可以比较一下二阶差分的效果：
plt.show()
#这里进行二阶差分
#fig = plt.figure(figsize=(12, 8))#创建一个大小为12*8的图
#ax2 = fig.add_subplot(111)
#diff2 = data.diff(2)#diff函数进行两次差分
#diff2.plot(ax=ax2)#ax=ax2表示将一个已经存在的subplot对象，将这个图赋给变量ax2，并将其
#指定
#plt.show()

#由下图可以看出来一阶跟二阶的差分差别不是很大，所以可以把差分次数d设置为1，上面的一阶和二阶程序我们注释掉

#这里我们使用一阶差分的时间序列
#3.接下来我们要找到ARIMA模型中合适的p和q值：
data1 = data.diff(1)
data1.dropna(inplace=True)#数据预处理，删除有缺失值的行，最后生成data1，inplace=True表示是否替换，true表示替换，直接在原数据文件修改，不产生副本
#加上这一步，不然后面画出ueacf和pacf图会是一条直线

#第一步：先检查平稳序列的自相关图和偏自相关图
fig = plt.figure(figsize=(12, 8))#figure()函数表示：创建一个新的绘图窗口，规定图的大小
ax1 = fig.add_subplot(211)#add.subpolt()函数创建组合窗口，2两行，1，一列，1表示第一个图
fig = sm.graphics.tsa.plot_acf(data1,lags=40,ax=ax1)#
#lags 表示滞后的阶数
#第二步：下面分别得到acf 图和pacf 图
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(data1, lags=40,method='ywm',ax=ax2)
#plt.show()'''重点记住plot_acf()函数，用来画自相关图，lags参数默认为40，data1表示用来显示自相关性的数据'''
#由上图可知，我们可以分别用ARMA(7,0)模型、ARMA(7,1)模型、ARMA(8,0)模型等来拟合找出最佳模型：
#第三步:找出最佳模型ARMA

'''源代码有问题，版本不一致
arma_mod1 = sm.tsa.ARMA(data1,(7,0)).fit()
print(arma_mod1.aic, arma_mod1.bic, arma_mod1.hqic)
arma_mod2 = sm.tsa.ARMA(data1,(7,1)).fit()
print(arma_mod2.aic, arma_mod2.bic, arma_mod2.hqic)
arma_mod3 = sm.tsa.ARMA(data1,(8,0)).fit()
print(arma_mod3.aic, arma_mod3.bic, arma_mod3.hqic)'''


#from statsmodels.tsa.arima_model import ARIMA
#from statsmodels.tsa.arima.model import ARIMA

arma_mod1 = sm.tsa.ARIMA(data1,order=(7,0,0)).fit()
print(arma_mod1.aic, arma_mod1.bic, arma_mod1.hqic)
arma_mod2 = sm.tsa.ARIMA(data1,order=(7,0,1)).fit()
print(arma_mod2.aic, arma_mod2.bic, arma_mod2.hqic)
arma_mod3 = sm.tsa.ARIMA(data1,order=(8,0,0)).fit()
print(arma_mod3.aic, arma_mod3.bic, arma_mod3.hqic)


#由上面可以看出ARMA(7,0)模型最佳

#第四步：进行模型检验
#首先对ARMA(7,0)模型所产生的残差做自相关图

resid = arma_mod1.resid#获得模型的残差
#一定要加上这个变量赋值语句，不然会报错resid is not defined 
'''  这行代码用于获取 ARMA 模型 `arma_mod1` 残差（residual）序列，`resid` 是该序列的名称，可以根据需要进行修改。

在时间序列分析中，残差序列是指利用模型预测值与真实值之差得到的一组残差（即预测误差）数据。这些数据可以用于检查模型拟合的好坏，
如果模型拟合良好，则残差序列中不应该存在显著的自相关性和异方差性（即残差的方差不应该与时间有关），否则可能说明模型还需要调整
或者存在未建模的变量。在建立好模型后，我们通常会对残差序列进行进一步的检验，对残差序列进行自相关性和噪声方差的分析，以确定
拟合模型的有效性和优化程度。 '''

fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(),lags=40,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40,ax=ax2)
plt.show()
#接着做德宾-沃森（D-W）检验
print(sm.stats.durbin_watson(arma_mod1.resid.values))
#得出来结果是不存在自相关性的
#2.033219223457169表示不存在或极少程度的自相关性

#再观察是否符合正态分布,这里用qq图
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q',ax=ax, fit=True)
'''这三行代码绘制QQ图，检验残差是否在一条直线上，在就符合正态分布，通过残差检验
fig=表示绘制一张图
ax=表示在组合图111中添加一个子图像ax
fig=qqplot()表示绘制qq图，其中的参数resid是残差序列，line='q'是绘制qq图时显示的
参考线，ax=ax表示将QQ图添加到ax中， fit=True表示在图像上绘制拟合线，fig是将qqplot产生的
图形对象赋值给变量fig,fig有助于后续对图形进行修改或保存,一般不需要这样做，因为直接调用plt.show()
即可显示图形'''
plt.show()


#最后用Ljung-Box检验:检验的结果就是看最后一列前十二行的检验概率（一般观察滞后1~12阶），
#如果检验概率小于给定的显著性水平，比如0.05、0.10等就拒绝原假设，其原假设是相关系数为零。
#就结果来看，前12阶的P值都是大于0.05，所以在0.05的显著性水平下，不拒绝原假设，即残差为白噪声序列。
r,q,p = sm.tsa.acf(resid.values.squeeze(),qstat=True)
data2 = np.c_[range(1,41), r[1:], q, p]
table= pd.DataFrame(data2, columns=[ 'lag','AC','Q','Prob(>Q)'])
print(table.set_index('lag'))

#第五步：平稳模型预测,对未来十年进行预测
predict_y =arma_mod1.predict('1990', '2000', dynamic=True)
print(predict_y)
fig, ax = plt.subplots(figsize=(12,8))
ax = data1.loc['1901':].plot(ax=ax)
predict_y.plot(ax=ax)

#第六步：使用ARIMA模型对原始序列进行预测
model = ARIMA(data,order=(7,0,0)) #导入ARIMA模型
result = model.fit(disp=-1)
#print(result.summary())
result.conf_int()#模型诊断，可以发现所有的系数置信区间都不为0；即在5%的置信水平下，所有的系数都是显著的，即模型通过检验。

#画出时序图
fig, ax = plt.subplots(figsize=(12, 10))
ax = data.loc['1901':].plot(ax=ax)   #注意起点是从1901开始
fig = result.plot_predict(5, 100)  #因为前面是90个数，所以加上预测的10个就是100
plt.show()   #数据预测并画图

#预测原始序列的未来10年数据
pred = result.predict(start = 90, end = 99, dynamic = True)
pred

本节的学习链接是：(39条消息) 利用ARIMA模型对时间序列进行分析的经典案例（详细代码）_时间序列模型经典案例_小白掌柜的博客-CSDN博客