5.5 时间序列预测

ARIMA模型

时间序列模型 ARIMA

ARIMA模型运用的流程

根据时间序列的散点图、自相关函数和偏自相关函数图识别其平稳性。
对非平稳的时间序列数据进行平稳化处理。直到处理后的自相关函数和偏自相关函数的数值非显著非零。
根据所识别出来的特征建立相应的时间序列模型。平稳化处理后，若偏自相关函数是截尾的，而自相关函数是拖尾的，则建立AR模型；若偏自相关函数是拖尾的，而自相关函数是截尾的，则建立MA模型；若偏自相关函数和自相关函数均是拖尾的，则序列适合ARMA模型。
参数估计，检验是否具有统计意义。
假设检验，判断（诊断）残差序列是否为白噪声序列。
利用已通过检验的模型进行预测。

## 加载包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
## 图像在jupyter notebook中显示
%matplotlib inline
## 显示的图片格式（mac中的高清格式），还可以设置为"bmp"等格式
%config InlineBackend.figure_format = "retina"
## 输出图显示中文
from matplotlib.font_manager import FontProperties
fonts = FontProperties(fname = "D:\Desktop\python在机器学习中的应用\方正粗黑宋简体.ttf",size=14)
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn import metrics
from sklearn.model_selection import train_test_split
## 忽略提醒
import warnings
warnings.filterwarnings("ignore")

## 读取飞机乘客数量的数据集
Airp = pd.read_csv("D:\Desktop\python在机器学习中的应用\AirPassengers.csv",index_col="Month")
Airp = Airp.astype("float64",copy=False) ## 调整数据类型
## 调整数据索引
Airp.index = pd.date_range(start=Airp.index.values[0],periods=Airp.shape[0],freq="MS")
Airp.head(5)

指定了索引列index_col=“Month”，为了方便使用statamodel库，将数据转换为*“float64”，使用astype*方法，并且重新制定数据的索引时间格式，*pd.date_range（）*函数用来生成一个时间序列，start指定开始时间，periods指定向后生成序列的长度，*freq=“MS”*表示时间方式，为每个月的月初，绘制折线图

## 该数据集为1949年到1960年的数据集
Airp.plot(kind = "line",figsize=(12,6))
plt.grid("on")
plt.title("航空乘客数量走势",FontProperties = fonts)
plt.show()

数据呈现波形逐渐上升的趋势，说明数据集具有很强的季节性规律。所谓季节性规律就是数据每隔一段时间其变化规律和之前一样，例如一年四季气温变化

时间序列的白噪声检验

该检验用来检查序列是否为随机序列，如果是随机序列，那它们的值之间没有任何关系

## 使用LB检验来检验序列是否为白噪声，原假设为在延迟期数内序列之间相互独立。
lags = [4,8,12,16,20,24,28,32,36,64,128]
LB = sm.stats.diagnostic.acorr_ljungbox(Airp,lags = lags)
LB = pd.DataFrame(data = np.array(LB).T,columns=["LBvalue","P-value"])
LB["lags"] = lags
LB

延迟阶数 [4,8,12,16,20,24,28,32,36,64,128]
LB检验的P值小于0.05，则认为该数据不是随机数据，数据是有规律的
从上面的LB检验的P值，都远小于0.05，说明我们可以拒绝原假设，不能认为序列为白噪声。
白噪声检验后要对数据平稳性进行判断

##  序列的单位根检验，即检验序列的平稳性
dftest = sm.tsa.adfuller(Airp.Passengers,autolag='BIC')
print("adf:",dftest[0])
print("pvalue:",dftest[1])
print("usedlag:",dftest[2])
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
dfoutput

P值远远大于0.05，说明该序列不是平稳的，需要对数据集进行一阶差分，一阶差分后重新进行白噪声检验和单位根检验

## 对数据进行一阶差分
datadiff = Airp.diff().dropna()
## 将时间序列可视化
datadiff.plot(kind = "line",figsize=(12,6))
plt.grid("on")
plt.title('一阶差分',FontProperties = fonts)
plt.show()

## 一阶差分后序列的白噪声检验
## 该检验用来检查序列是否为随机序列，如果是随机序列，那它们的值之间没有任何关系
## 使用LB检验来检验序列是否为白噪声，原假设为在延迟期数内序列之间相互独立。
lags = [4,8,12,16,20,24,28,32,36,64,128]
LB = sm.stats.diagnostic.acorr_ljungbox(datadiff,lags = lags)
LB = pd.DataFrame(data = np.array(LB).T,columns=["LBvalue","P-value"])
LB["lags"] = lags
print("差分后序列的白噪声检验:\n",LB)
print("-------------------")
##  差分后序列的单位根检验，即检验序列的平稳性
dftest = sm.tsa.adfuller(datadiff.Passengers,autolag='AIC')
print("adf:",dftest[0])
print("pvalue:",dftest[1])
print("usedlag:",dftest[2])

一阶差分后：进行白噪声检验，可以排除序列为白噪声；进行单位根检验时，P－value=0.054,非常接近于0.05，在置信度为90%的情况下，认为改序列是平稳的。可以确定参数 I＝1

## 差分后序列的自相关和偏相关图
fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(datadiff, lags=80, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(datadiff, lags=80, ax=ax2)
plt.subplots_adjust(hspace = 0.3)
plt.ylim(-0.5,1.5)
plt.show()

根据自相关系数图可以确定p，根据偏相关系数图可以确定q，图中可以确定p=2;q=10

并且从图中我们也能看出相关系数的周期性的变化，通过观察自相关图，周期在12左右波动，我们可以将周期定位［12，24］等12的倍数计算多个模型，并选出最优的。seasonal_order=(1,0,8,12)中的另外三个参数的确定

ARIMA 模型

## 我们可以通过循环计算多个模型，通过比较bic值对模型定阶
#一般阶数不超过length/10
pmax = 12
#一般阶数不超过length/10
qmax = 15
#bic矩阵
bic_matrix = []
for p in range(pmax+1):tmp = []for q in range(qmax+1):# 存在部分报错，所以用try来跳过报错。try: tmp.append(sm.tsa.ARIMA(Airp, (p,1,q)).fit().bic)except:tmp.append(None)bic_matrix.append(tmp)
print("模型迭代结束")#从中可以找出最小值
bic_matrix = pd.DataFrame(bic_matrix)
# print("bic矩阵:\n",bic_matrix)
#先用stack展平，然后用idxmin找出最小值位置。
p,q = bic_matrix.stack().idxmin()
print("比较合适的p:",p)
print("比较合适的q:",q)
bic_matrixARIMAmol = sm.tsa.ARIMA(Airp,order=(4,1,9)).fit()
print(ARIMAmol.summary2())fig, ax = plt.subplots(figsize=(12, 8))
ax = Airp.ix['1949':].plot(ax=ax,marker = "*")
fig = ARIMAmol.plot_predict(start=143,end=168,dynamic=True,ax=ax, plot_insample=False)
plt.xlabel("Year")
plt.ylabel("Number")
plt.title("ARIMA(4,1,9)")
plt.legend(loc = "upper left")
plt.show()

ARIMAX模型

## 我们可以通过循环计算多个模型，通过比较bic值对模型定阶
#一般阶数不超过length/10 （我了减少运行时间）
pmax = 10
#一般阶数不超过length/10
qmax = 10
#bic矩阵bic_matrix = []
for p in range(pmax+1):tmp = []for q in range(qmax+1):# 存在部分报错，所以用try来跳过报错。try:
#             tmp.append(sm.tsa.ARIMA(Airp, (p,1,q)).fit().bic)tmp.append(SARIMAX(Airp,order=(p,1,q),seasonal_order=(0,0,0,12)).fit().bic)except:tmp.append(None)bic_matrix.append(tmp)
print("模型迭代结束")#从中可以找出最小值
bic_matrix = pd.DataFrame(bic_matrix)
print("bic矩阵:\n",bic_matrix)
#先用stack展平，然后用idxmin找出最小值位置。
p,q = bic_matrix.stack().idxmin()
print("比较合适的p:",p)
print("比较合适的q:",q)modl = SARIMAX(Airp,order=(10,1,0),seasonal_order=(0,0,0,12)).fit()
print(modl.summary())## 使用模型进行预测
datapre = modl.predict(start=120,end= 160)## 绘制原始数据和训练数据
Airp.plot(kind = "line",figsize=(12,6))
datapre.plot(kind = "line",color = "r")
# FontProperties = fonts
plt.title("ARIMAX((10,1,0),(0,0,0,12))")
plt.show()