Python3中使用ARIMA进行时间序列预测的指南

介绍

时间序列提供了预测未来价值的机会。基于以前的值，时间序列可用于预测经济、天气和容量规划等方面的趋势。时间序列数据的特定属性意味着通常需要专门的统计方法。

在本教程中，我们的目标是生成可靠的时间序列预测。我们将首先介绍和讨论自相关、平稳性和季节性的概念，然后应用最常用的时间序列预测方法之一，即 ARIMA。

Python 中可用于对时间序列的未来点进行建模和预测的方法之一称为 SARIMAX，它代表 Seasonal AutoRegressive Integrated Moving Averages with eXogenous regressors。在这里，我们将主要关注 ARIMA 组件，它用于拟合时间序列数据，以更好地理解和预测时间序列中的未来点。

先决条件

本指南将介绍如何在本地桌面或远程服务器上进行时间序列分析。处理大型数据集可能会占用大量内存，因此无论哪种情况，计算机都需要至少 2GB 内存 来执行本指南中的一些计算。

为了充分利用本教程，熟悉时间序列和统计数据可能会有所帮助。

在本教程中，我们将使用 Jupyter Notebook 来处理数据。如果您还没有它，您应该按照我们的教程安装和设置适用于 Python 3 的 Jupyter Notebook。

第 1 步 — 安装软件包

为了设置我们的时间序列预测环境，让我们首先进入我们的本地编程环境或基于服务器的编程环境：

cd environments

. my_env/bin/activate

从这里，让我们为我们的项目创建一个新目录。我们将它命名为ARIMA，然后进入目录。如果您给项目起一个不同的名称，请务必在整个指南中用您的名称替换 ARIMA

mkdir ARIMA
cd ARIMA

本教程将需要 warnings、itertools、pandas、numpy、matplotlib 和 statsmodels 库。 warnings 和 itertools 库包含在标准 Python 库集中，因此您不需要安装它们。

与其他 Python 包一样，我们可以使用 pip 安装这些要求。我们现在可以安装 pandas、statsmodels 和数据绘图包 matplotlib。它们的依赖项也将被安装：

pip install pandas numpy statsmodels matplotlib

至此，我们现在已设置好开始使用已安装的软件包。

第 2 步 — 导入包和加载数据

要开始处理我们的数据，我们将启动 Jupyter Notebook：

jupyter notebook

要创建新的笔记本文件，请从右上角的下拉菜单中选择 New > Python 3：

这将打开一个笔记本。

作为最佳实践，首先在笔记本顶部导入您需要的库：

import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

我们还为我们的绘图定义了一个 matplotlib 样式为五三八。

我们将使用一个名为“来自美国夏威夷莫纳罗亚天文台连续空气样本的大气二氧化碳”的数据集，该数据集收集了 1958 年 3 月至 2001 年 12 月的二氧化碳样本。我们可以按如下方式引入这些数据：

data = sm.datasets.co2.load_pandas()
y = data.data

在继续之前，让我们对数据进行一些预处理。每周数据可能很难处理，因为它的时间较短，所以让我们改用每月平均值。我们将使用 resample 函数进行转换。为简单起见，我们还可以使用 fillna() 函数来确保我们的时间序列中没有缺失值。

# The 'MS' string groups the data in buckets by start of the month
y = y['co2'].resample('MS').mean()

# The term bfill means that we use the value before filling in missing values
y = y.fillna(y.bfill())

print(y)

Outputco2
1958-03-01  316.100000
1958-04-01  317.200000
1958-05-01  317.433333
...
2001-11-01  369.375000
2001-12-01  371.020000

让我们探索这个时间序列 e 作为数据可视化：

y.plot(figsize=(15, 6))
plt.show()

当我们绘制数据时，会出现一些可区分的模式。时间序列具有明显的季节性规律，整体呈上升趋势。

要了解有关时间序列预处理的更多信息，请参阅“A Guide to Time Series Visualization with Python 3”，其中对上述步骤进行了更详细的描述。

现在我们已经转换并探索了我们的数据，让我们继续使用 ARIMA 进行时间序列预测。

第 3 步 — ARIMA 时间序列模型

时间序列预测中最常用的方法之一被称为 ARIMA 模型，它代表 AutoregRessive Iintegrated M[ X167X]oving A平均。 ARIMA 是一个可以拟合时间序列数据的模型，以便更好地理解或预测序列中的未来点。

有三个不同的整数（p、d、q）用于参数化 ARIMA 模型。因此，ARIMA 模型用符号 ARIMA(p, d, q) 表示。这三个参数一起说明了数据集中的季节性、趋势和噪声：

p 是模型的 auto-regressive 部分。它允许我们将过去值的影响合并到我们的模型中。直观地说，这类似于说如果过去 3 天一直很暖，明天可能会暖和。
d是模型的集成部分。这包括模型中包含差异量的项（即要从当前值中减去的过去时间点的数量）以应用于时间序列。直观地说，这类似于说如果过去三天的温度差异非常小，明天可能是相同的温度。
q是模型的移动平均线部分。这允许我们将模型的误差设置为在过去的先前时间点观察到的误差值的线性组合。

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

由于涉及多个调整参数，季节性 ARIMA 方法可能看起来令人生畏。在下一节中，我们将描述如何自动化识别季节性 ARIMA 时间序列模型的最佳参数集的过程。

第 4 步 — ARIMA 时间序列模型的参数选择

在寻找使用季节性 ARIMA 模型拟合时间序列数据时，我们的第一个目标是找到优化感兴趣指标的 ARIMA(p,d,q)(P,D,Q)s 的值。实现这一目标有许多指导方针和最佳实践，但 ARIMA 模型的正确参数化可能是一个艰苦的手动过程，需要领域专业知识和时间。 R 等其他统计编程语言提供了自动解决此问题的方法，但这些尚未移植到 Python。在本节中，我们将通过编写 Python 代码以编程方式为我们的 ARIMA(p,d,q)(P,D,Q)s 时间序列模型选择最佳参数值来解决此问题。

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

# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)

# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))

# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

print('Examples of parameter combinations for Seasonal ARIMA...')
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))

OutputExamples of parameter combinations for Seasonal ARIMA...
SARIMAX: (0, 0, 1) x (0, 0, 1, 12)
SARIMAX: (0, 0, 1) x (0, 1, 0, 12)
SARIMAX: (0, 1, 0) x (0, 1, 1, 12)
SARIMAX: (0, 1, 0) x (1, 0, 0, 12)

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

在评估和比较配备不同参数的统计模型时，可以根据其与数据的拟合程度或其准确预测未来数据点的能力，对每个模型进行排名。我们将使用 AIC（Akaike 信息标准）值，使用 statsmodels 拟合的 ARIMA 模型可以方便地返回该值。 AIC 衡量模型与数据的拟合程度，同时考虑模型的整体复杂性。与使用较少特征以实现相同拟合优度的模型相比，在使用大量特征时非常适合数据的模型将被分配更高的 AIC 分数。因此，我们有兴趣找到产生最低 AIC 值的模型。

下面的代码块遍历参数组合，并使用 statsmodels 中的 SARIMAX 函数来拟合相应的季节性 ARIMA 模型。此处，order 参数指定 (p, d, q) 参数，而 seasonal_order 参数指定季节性 ARIMA 模型的 (P, D, Q, S) 季节性分量。在拟合每个 SARIMAX() 模型后，代码会打印出其各自的 AIC 分数。

warnings.filterwarnings("ignore") # specify to ignore warning messages

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

由于某些参数组合可能导致数值错误指定，我们明确禁用警告消息以避免警告消息过载。这些错误指定也可能导致错误并引发异常，因此我们确保捕获这些异常并忽略导致这些问题的参数组合。

上面的代码应该会产生以下结果，这可能需要一些时间：

OutputSARIMAX(0, 0, 0)x(0, 0, 1, 12) - AIC:6787.3436240402125
SARIMAX(0, 0, 0)x(0, 1, 1, 12) - AIC:1596.711172764114
SARIMAX(0, 0, 0)x(1, 0, 0, 12) - AIC:1058.9388921320026
SARIMAX(0, 0, 0)x(1, 0, 1, 12) - AIC:1056.2878315690562
SARIMAX(0, 0, 0)x(1, 1, 0, 12) - AIC:1361.6578978064144
SARIMAX(0, 0, 0)x(1, 1, 1, 12) - AIC:1044.7647912940095
...
...
...
SARIMAX(1, 1, 1)x(1, 0, 0, 12) - AIC:576.8647112294245
SARIMAX(1, 1, 1)x(1, 0, 1, 12) - AIC:327.9049123596742
SARIMAX(1, 1, 1)x(1, 1, 0, 12) - AIC:444.12436865161305
SARIMAX(1, 1, 1)x(1, 1, 1, 12) - AIC:277.7801413828764

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

第 5 步 — 拟合 ARIMA 时间序列模型

使用网格搜索，我们已经确定了为我们的时间序列数据生成最佳拟合模型的参数集。我们可以继续更深入地分析这个特定的模型。

我们将从将最优参数值插入新的 SARIMAX 模型开始：

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])

Output==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.3182      0.092      3.443      0.001       0.137       0.499
ma.L1         -0.6255      0.077     -8.165      0.000      -0.776      -0.475
ar.S.L12       0.0010      0.001      1.732      0.083      -0.000       0.002
ma.S.L12      -0.8769      0.026    -33.811      0.000      -0.928      -0.826
sigma2         0.0972      0.004     22.634      0.000       0.089       0.106
==============================================================================

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

在拟合季节性 ARIMA 模型（以及与此相关的任何其他模型）时，运行模型诊断以确保不违反模型所做的任何假设非常重要。 plot_diagnostics 对象允许我们快速生成模型诊断并调查任何异常行为。

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

我们主要关心的是确保我们模型的残差是不相关的，并且正态分布为零均值。如果季节性 ARIMA 模型不满足这些属性，则很好地表明它可以进一步改进。

在这种情况下，我们的模型诊断表明模型残差正态分布基于以下几点：

在右上角的图中，我们看到红色的 KDE 线紧跟 N(0,1) 线（其中 N(0,1)）是均值为 0 和 1 的标准差）。这很好地表明残差是正态分布的。
左下角的 qq-plot 显示残差的有序分布（蓝点）遵循从具有 N(0, 1) 的标准正态分布中提取的样本的线性趋势。同样，这强烈表明残差是正态分布的。
随时间变化的残差（左上图）没有显示任何明显的季节性，并且似乎是白噪声。这由自相关证实（即 correlogram) 图在右下角，这表明时间序列残差与自身滞后版本的相关性较低。

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

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

第 6 步 - 验证预测

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

pred = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=False)
pred_ci = pred.conf_int()

上面的代码要求预测从 1998 年 1 月开始。

dynamic=False 参数确保我们产生提前一步的预测，这意味着每个点的预测都是使用该点之前的完整历史生成的。

我们可以绘制 CO2 时间序列的真实值和预测值，以评估我们的表现。请注意我们如何通过切片日期索引来放大时间序列的末尾。

ax = y['1990':].plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)

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()

总体而言，我们的预测与真实值非常吻合，呈现整体增长趋势。

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

y_forecasted = pred.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)))

OutputThe Mean Squared Error of our forecasts is 0.07

我们提前一步预测的 MSE 产生的值为 0.07，该值非常低，因为它接近于 0。 MSE 为 0 表示估计器以完美的准确度预测参数的观测值，这将是一个理想的场景，但通常不可能。

但是，可以使用动态预测更好地表示我们的真实预测能力。在这种情况下，我们只使用从时间序列到某个时间点的信息，然后使用之前预测时间点的值生成预测。

在下面的代码块中，我们指定从 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()

绘制时间序列的观察值和预测值，我们看到即使使用动态预测，整体预测也是准确的。所有预测值（红线）与基本事实（蓝线）非常匹配，并且完全在我们预测的置信区间内。

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()

再次，我们通过计算 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)))

OutputThe Mean Squared Error of our forecasts is 1.01

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

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

第 7 步 - 生成和可视化预测

在本教程的最后一步，我们将描述如何利用我们的季节性 ARIMA 时间序列模型来预测未来值。我们的时间序列对象的 get_forecast() 属性可以计算提前指定步数的预测值。

# Get forecast 500 steps ahead in future
pred_uc = results.get_forecast(steps=500)

# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

我们可以使用此代码的输出来绘制时间序列和对其未来值的预测。

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')
ax.set_ylabel('CO2 Levels')

plt.legend()
plt.show()

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

随着我们对未来的进一步预测，我们对自己的价值观变得不那么自信是很自然的。这反映在我们的模型生成的置信区间上，随着我们进一步走向未来，置信区间会变得更大。

结论

在本教程中，我们描述了如何在 Python 中实现季节性 ARIMA 模型。我们广泛使用了 pandas 和 statsmodels 库，并展示了如何运行模型诊断，以及如何生成 CO2 时间序列的预测。

以下是您可以尝试的其他一些事情：

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

如需更多练习，您还可以尝试加载另一个时间序列数据集以生成您自己的预测。