基于Python3 使用ARIMA进行时间序列预测

灰太狼 2022-05-16 03:36 1860阅读 0赞

## **引言：** ##

在本文章中，我们将提供可靠的时间序列预测。我们将首先介绍和讨论自相关，平稳性和季节性的概念，并继续应用最常用的时间序列预测方法之一，称为ARIMA。

## 简介： ##

时间序列提供了预测未来价值的机会。 基于以前的价值观，可以使用时间序列来预测经济，天气和能力规划的趋势，其中仅举几例。 时间序列数据的具体属性意味着通常需要专门的统计方法。

在本教程中，我们将针对时间序列产生可靠的预测。 我们将首先介绍和讨论自相关，平稳性和季节性的概念，并继续应用最常用的时间序列预测方法之一，称为ARIMA。

用于建模和预测时间序列未来点的Python中的一种方法被称为**SARIMAX** ，其代表**具有eXogenous回归的季节性自动反馈集成移动平均值** 。 在这里，我们将主要关注ARIMA组件，该组件用于适应时间序列数据，以更好地了解和预测时间序列中的未来点。

## 1 实验条件 ##

本将介绍如何在本地桌面或远程服务器上进行时间序列分析。 使用大型数据集可能是内存密集型的，所以在这两种情况下，计算机将至少需要**2GB的内存**来执行本指南中的一些计算。

我使用的IDE: PyCharm ，实验Python版本：3.6

实验所用的库： pandas(操作数据)    statsmodels(统计分析库)    matplotlib(可视化)   itertools(迭代工具)   warnings(设置警告)

引入包：

import pandas as pd
    import statsmodels.api as sm
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import itertools
    import warnings

## 2 数据集： ##

我们将使用一个名为“来自美国夏威夷Mauna Loa天文台的连续空气样本的大气二氧化碳”的数据集，该数据集从1958年3月至2001年12月期间收集了二氧化碳样本。这些数据可以通过statsmodels.api 引进来；

#引入数据
    data=sm.datasets.co2.load()
    #对引入的数据进行变形
    index = pd.DatetimeIndex(start=data.data['date'][0].decode('utf-8'),
                                 periods=len(data.data),freq='W-SAT')
    co2=pd.DataFrame(data.data['co2'], index=index, columns=['co2'])
    print(co2.index)#检查co2的索引

运行结果：

![70][]

我们稍后再进行一些数据处理。从上面的的结果可以看出，每周数据可能很棘手，因为它是一个很短的时间，所以让我们使用每月平均值。 我们将使用resample函数进行转换。 另外当数据量比较大的时候，容易存在缺失值，我们需要先检测是否存在缺失值，如果缺失值过多，则说明数据不可用，如果数据缺失的不太多，可以用一些方法进行拟合补全。为了简单起见，这里 我们使用[`fillna()`函数][fillna]来确保我们的时间序列中没有缺少值。当然可以自己训练模型根据预测结果，进行自动补全。

y=co2["co2"].resample("MS").mean()#获得每个月的平均值
    print(y.isnull().sum)#5个 检测空白值
    #处理数据中的缺失项
    y=y.fillna(y.bfill())#填充缺失值

经过上部，我们已经补全了所有的空白值。

## 3 初始数据可视化： ##

我们可以先看看原始的数据，可视化的展示：

plt.figure(figsize=(15,6))
     plt.title("原始数据",loc="center",fontsize=20)
     plt.plot(y)

实验结果：

![70 1][]

由实验结果可以看出：当我们绘制数据时，会出现一些可区分的模式。 时间序列具有明显的季节性格局，总体呈上升趋势。

现在我们已经转换和探索了原始的数据，接下来我们继续使用ARIMA进行时间序列预测。

## 4 ARIMA时间序列模型 ##

时间序列预测中最常用的方法之一就是被称为ARIMA模型，它代表了**A** utoreg **R** essive综合**M** oving **A版本** 。 ARIMA是可以适应时间序列数据的模型，以便更好地了解或预测系列中的未来点。

有三个不同的整数（ `p` ， `d` ， `q` ）用于参数化ARIMA模型。 因此，ARIMA模型用符号`ARIMA(p, d, q)` 。 这三个参数共计数据集中的季节性，趋势和噪音：

*  `p`是模型的*自回归*部分。 它允许我们将过去价值观的影响纳入我们的模型。 直观地说，这将是类似的，表示如果过去3天已经变暖，明天可能会变暖。
 *  `d`是模型的*集成*部分。 这包括模型中包含差异量（即从当前值减去的过去时间点的数量）以适用于时间序列的术语。 直观地说，这将类似于说如果过去三天的温度差异非常小，明天可能会有相同的温度。
 *  `q`是模型的*移动平均*部分。 这允许我们将模型的误差设置为过去以前时间点观察到的误差值的线性组合。

在处理季节性影响时，我们利用*季节性* ARIMA，表示为**`ARIMA(p,d,q)(P,D,Q)s`** 。 这里， `(p, d, q)`是上述非季节性参数，而`(P, D, Q)`遵循相同的定义，但适用于时间序列的季节分量。 术语`s`是时间序列的周期（季度为`4` ，年度为`12` ，等等）。

由于所涉及的多个调整参数，季节性ARIMA方法可能会令人望而生畏。所以这里我们采用年度的预测。 在下一节中，我们将介绍如何自动化识别季节性ARIMA时间序列模型的最优参数集的过程。

## 5  ARIMA时间序列模型的参数选择 ##

当考虑使用季节性ARIMA模型拟合时间序列数据时，我们的第一个目标是找到优化感兴趣度量的`ARIMA(p,d,q)(P,D,Q)s`的值。 实现这一目标有许多指导方针和最佳实践，但ARIMA模型的正确参数化可能是一个需要领域专长和时间的艰苦的手工过程。

这里我们使用“网格搜索”来迭代地探索参数的不同组合。 对于参数的每个组合，我们使用`statsmodels`模块的`SARIMAX()`函数拟合一个新的季节性ARIMA模型，并评估其整体质量。 一旦我们探索了参数的整个范围，我们的最佳参数集将是我们感兴趣的标准产生最佳性能的参数。 我们开始生成我们希望评估的各种参数组合：

#找合适的p d q
    #初始化 p d q
    p=d=q=range(0,2)
    print("p=",p,"d=",d,"q=",q)
    #产生不同的pdq元组,得到 p d q 全排列
    pdq=list(itertools.product(p,d,q))
    print("pdq:\n",pdq)
    seasonal_pdq=[(x[0],x[1],x[2],12) for x in pdq]
    print('SQRIMAX:{} x {}'.format(pdq[1],seasonal_pdq[1]))

![70 2][]

我们现在可以使用上面定义的参数三元组来自动化不同组合对ARIMA模型进行培训和评估的过程。 在统计和机器学习中，这个过程被称为模型选择的网格搜索（或超参数优化）。

在评估和比较配备不同参数的统计模型时，可以根据数据的适合性或准确预测未来数据点的能力，对每个参数进行排序。 我们将使用`AIC` （Akaike信息标准）值，该值通过使用`statsmodels`安装的ARIMA型号方便地返回。 `AIC`衡量模型如何适应数据，同时考虑到模型的整体复杂性。 在使用大量功能的情况下，适合数据的模型将被赋予比使用较少特征以获得相同的适合度的模型更大的AIC得分。 因此，我们有兴趣找到产生最低`AIC`值的模型。

下面的代码块通过参数的组合来迭代，并使用`SARIMAX`函数来适应相应的季节性ARIMA模型。 这里， `order`参数指定`(p, d, q)`参数，而`seasonal_order`参数指定季节性ARIMA模型的`(P, D, Q, S)`季节分量。 在安装每个`SARIMAX()`模型后，代码打印出其各自的`AIC`得分。

for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(y,
                                                order=param,
                                                seasonal_order=param_seasonal,
                                                enforce_stationarity=False,
                                                enforce_invertibility=False)
    
                results = mod.fit()
    
                print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
            except:
                continue

运行结果：

![70 3][]

【注：实验结果会有很多的数据，这里只截取了AIC的结果，且需要执行一会儿】

我们的代码的输出表明， `SARIMAX(1, 1, 1)x(1, 1, 1, 12)`产生最低的`AIC`值为277.78。 因此，我们认为这是我们考虑过的所有模型中的最佳选择。

## 6 配置ARIMA时间序列模型 ##

##    利用以上的估计，使该模型为 (p,d,q)=(1,1,1) ##

mod = sm.tsa.statespace.SARIMAX(y,
                                    order=(1, 1, 1),
                                    seasonal_order=(1, 1, 1, 12),
                                    enforce_stationarity=False,
                                    enforce_invertibility=False)
    
    results = mod.fit()
    
    print(results.summary().tables[1])

结果：

![70 4][]

由SARIMAX的输出产生的`SARIMAX`返回大量的信息，但是我们将注意力集中在系数表上。 `coef`列显示每个特征的重量（即重要性）以及每个特征如何影响时间序列。 `P>|z|` 列通知我们每个特征重量的意义。 这里，每个重量的p值都低于或接近`0.05` ，所以在我们的模型中保留所有权重是合理的。

在适合季节性ARIMA模型（以及任何其他模型）的情况下，运行模型诊断是非常重要的，以确保没有违反模型的假设。 `plot_diagnostics`对象允许我们快速生成模型诊断并调查任何异常行为。

results.plot_diagnostics(figsize=(15, 12))
    plt.show()

运行结果：

![70 5][]

我们的主要关切是确保我们的模型的残差是不相关的，并且平均分布为零。 如果季节性ARIMA模型不能满足这些特性，这是一个很好的迹象，可以进一步改善。

在这种情况下，我们的模型诊断表明，模型残差正常分布如下：

*  在右上图中，我们看到红色`KDE`线与`N(0,1)`行（其中`N(0,1)` ）是正态分布的标准符号，平均值`0` ，标准偏差为`1` ） 。 这是残留物正常分布的良好指示。
 *  左下角的[qq图][qq]显示，残差（蓝点）的有序分布遵循采用`N(0, 1)`的标准正态分布采样的线性趋势。 同样，这是残留物正常分布的强烈指示。
 *  随着时间的推移（左上图）的残差不会显示任何明显的季节性，似乎是白噪声。 这通过右下角的自相关（即相关图）来证实，这表明时间序列残差与其本身的滞后版本具有低相关性。

这些观察结果使我们得出结论，我们的模型产生了令人满意的合适性，可以帮助我们了解我们的时间序列数据和预测未来价值。

虽然我们有一个令人满意的结果，我们的季节性ARIMA模型的一些参数可以改变，以改善我们的模型拟合。 例如，我们的网格搜索只考虑了一组受限制的参数组合，所以如果我们拓宽网格搜索，我们可能会找到更好的模型。

## 7 验证预测 ##

我们已经获得了我们时间序列的模型，现在可以用来产生预测。 我们首先将预测值与时间序列的实际值进行比较，这将有助于我们了解我们的预测的准确性。 `get_prediction()`和`conf_int()`属性允许我们获得时间序列预测的值和相关的置信区间。

# #进行验证预测
    pred=results.get_prediction(start=pd.to_datetime('1998-01-01'),dynamic=False)
    pred_ci=pred.conf_int()
    print("pred ci:\n",pred_ci)#获得的是一个预测范围，置信区间
    print("pred:\n",pred)#为一个预测对象
    print("pred mean:\n",pred.predicted_mean)#为预测的平均值

上述规定需要从1998年1月开始进行预测。

`dynamic=False`参数确保我们产生一步前进的预测，这意味着每个点的预测都将使用到此为止的完整历史生成。

![70 6][]

我们可以绘制二氧化碳时间序列的实际值和预测值，以评估我们做得如何。 注意我们如何在时间序列的末尾放大日期索引。

#进行绘制预测图像
    ax=y['1990':].plot(label="observed")
    pred.predicted_mean.plot(ax=ax,label="static ForCast",alpha=.7,color='red',linewidth=5)
    #在某个范围内进行填充
    ax.fill_between(pred_ci.index,
                    pred_ci.iloc[:, 0],
                    pred_ci.iloc[:, 1], color='k', alpha=.2)
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
    plt.legend()
    plt.show()

![70 7][]

总体而言，我们的预测与真实价值保持一致，呈现总体增长趋势。

量化我们的预测的准确性也是有用的。 我们将使用MSE（均方误差），它总结了我们的预测的平均误差。 对于每个预测值，我们计算其到真实值的距离并对结果求平方。 结果需要平方，以便当我们计算总体平均值时，正/负差异不会相互抵消。

# 求取 MSE（均方误差）
    y_forecasted=pred.predicted_mean
    y_truth=y['1998-01-01':]
    mse=((y_forecasted-y_truth)**2).mean()
    print("MSE:",mse)

结果：

![70 8][]

我们前进一步预测的MSE值为`0.07` ，这是接近0的非常低的值。0的MSE是估计器将以完美的精度预测参数的观测值，这将是一个理想的场景但通常不可能。

然而，使用动态预测可以获得更好地表达我们的真实预测能力。 在这种情况下，我们只使用时间序列中的信息到某一点，之后，使用先前预测时间点的值生成预测。

### 在下面的代码块中，我们指定从1998年1月起开始计算动态预测和置信区间 ###

pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
    pred_dynamic_ci = pred_dynamic.conf_int()

绘制时间序列的观测值和预测值，我们看到即使使用动态预测，总体预测也是准确的。 所有预测值（红线）与地面真相（蓝线）相当吻合，并且在我们预测的置信区间内。

# #使用动态预测
    pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
    pred_dynamic_ci = pred_dynamic.conf_int()
    ax = y['1990':].plot(label='observed', figsize=(20, 15))
    pred_dynamic.predicted_mean.plot(label='Dynamic Forecast', ax=ax)
    
    ax.fill_between(pred_dynamic_ci.index,
                    pred_dynamic_ci.iloc[:, 0],
                    pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)
    
    ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('1998-01-01'), y.index[-1],
                     alpha=.1, zorder=-1)
    
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
    
    plt.legend()
    plt.show()

![70 9][]

我们再次通过计算MSE量化我们预测的预测性能：

# Extract the predicted and true values of our time-series
    y_forecasted = pred_dynamic.predicted_mean
    y_truth = y['1998-01-01':]
    
    # Compute the mean square error
    mse = ((y_forecasted - y_truth) ** 2).mean()
    print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

![70 10][]

从动态预测获得的预测值产生1.01的MSE。 这比前进一步略高，这是预期的，因为我们依赖于时间序列中较少的历史数据。

前进一步和动态预测都证实了这个时间序列模型是有效的。 然而，关于时间序列预测的大部分兴趣是能够及时预测未来价值观。

## 8 生成和可视化预测 ##

我们将介绍如何利用季节性ARIMA时间序列模型来预测未来的价值。 我们的时间序列对象的`get_forecast()`属性可以计算预先指定数量的步骤的预测值。

# Get forecast 500 steps ahead in future
    pred_uc = results.get_forecast(steps=200)#steps 可以代表（200/12）年左右
    
    # Get confidence intervals of forecasts
    pred_ci = pred_uc.conf_int()
    plt.title("预测",fontsize=15,color="red")
    ax = y.plot(label='observed', figsize=(20, 15))
    pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
    ax.fill_between(pred_ci.index,
                    pred_ci.iloc[:, 0],
                    pred_ci.iloc[:, 1], color='k', alpha=.25)
    ax.set_xlabel('Date',fontsize=15)
    ax.set_ylabel('CO2 Levels',fontsize=15)
    
    plt.legend()
    plt.show()

![70 11][]

现在可以使用我们生成的预测和相关的置信区间来进一步了解时间序列并预见预期。 我们的预测显示，时间序列预计将继续稳步增长。

随着我们对未来的进一步预测，我们对自己的价值观信心愈来愈自然。 这反映在我们的模型产生的置信区间，随着我们进一步走向未来，这个模型越来越大。

## 9 结论 ##

在本教程中，我们描述了如何在Python中实现季节性ARIMA模型。 我们广泛使用了`pandas`和`statsmodels`图书馆，并展示了如何运行模型诊断，以及如何产生二氧化碳时间序列的预测。

这里还有一些其他可以尝试的事情：

*  更改动态预测的开始日期，以了解其如何影响预测的整体质量。
 *  尝试更多的参数组合，看看是否可以提高模型的适合度。
 *  选择不同的指标以选择最佳模型。 例如，我们使用`AIC`测量来找到最佳模型，但是您可以寻求优化采样均方误差。

对于更多的实践，您还可以尝试加载另一个时间序列数据集来生成您自己的预测。

【本篇文章源自前辈的实验，因为前辈描述的较为详细，这里只是做了细微的调整。原文地址：[https://www.howtoing.com/a-guide-to-time-series-forecasting-with-arima-in-python-3/][https_www.howtoing.com_a-guide-to-time-series-forecasting-with-arima-in-python-3]】

[70]: /images/20220516/84f76b63b8e04afd8cb49251b744491e.png
[fillna]: https://www.howtoing.com/an-introduction-to-the-pandas-package-and-its-data-structures-in-python-3/#handling-missing-values
[70 1]: /images/20220516/4726d5b68c2a4036a3569cee0a291c39.png
[70 2]: /images/20220516/bf9fce5035b7489993a8119cbedb0345.png
[70 3]: /images/20220516/7692b59f226c431d9b579c88fc051e92.png
[70 4]: /images/20220516/189679c063924827910e8a8bfd7869d2.png
[70 5]: /images/20220516/4cfa7ab0a50a4244bcca6e5d2b99e3f2.png
[qq]: https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot
[70 6]: /images/20220516/59c6c192105e44d295b9e493d72c2bcf.png
[70 7]: /images/20220516/f2b7998d8cdb4ff1a2d134ae3c856a1e.png
[70 8]: /images/20220516/4e000c884941474481de3f9ab109cf62.png
[70 9]: /images/20220516/4da2b7e61d014ae4bc96a26c284ed0f9.png
[70 10]: /images/20220516/75719d531ad54ad9b01110b498f4775c.png
[70 11]: /images/20220516/f6e3c14a6b7443c5bc32ee0690f9046a.png
[https_www.howtoing.com_a-guide-to-time-series-forecasting-with-arima-in-python-3]: https://www.howtoing.com/a-guide-to-time-series-forecasting-with-arima-in-python-3/