引言

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

在本教程中,我们将致力于对时间序列进行可靠的预测。 我们将首先介绍和讨论自相关、平稳性和季节性的概念,然后应用最常用的时间序列预测方法之一,即 ARIMA。

在 Python 中可用于建模和预测时间序列未来点的方法之一是 SARIMAX,它代表带有外生回归变量的季节性自回归综合移动平均值。 这里,我们将主要关注 ARIMA 分量,它用于拟合时间序列数据,以更好地理解和预测时间序列中的未来点。

先决条件

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

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

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

第一步ー安装软件包

为了建立时间序列预测环境,让我们首先进入本地编程环境或基于服务器的编程环境:

  • cd environments
  • . my_env/bin/activate

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

  • mkdir ARIMA
  • cd ARIMA

本教程将要求使用警告、迭代工具、 pandas、 numpy、 matplotlib 和 statsmodels 库。 警告和迭代工具库包含在标准的 Python 库集中,所以您不需要安装它们。

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

  • pip install pandas numpy statsmodels matplotlib

现在,我们已经准备好开始使用已安装的包了。

步骤2ー导入包和载入数据

为了开始处理我们的数据,我们将启动 Jupyter Notebook:

  • jupyter notebook

要创建一个新的笔记本文件,从右上角的下拉菜单中选择 New python3:

这会打开一个笔记本。

作为最佳实践,从导入你笔记本顶部需要的库开始:

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 风格的538。

我们将使用一个名为“来自美国夏威夷莫纳罗亚火山天文台的连续空气样本的大气二氧化碳”的数据集,该数据集收集了1958年3月至2001年12月的二氧化碳样本。 我们可以提供以下数据:

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

在继续之前,让我们先对数据进行一些预处理。 周数据可能很难处理,因为它的时间更短,所以让我们使用月平均值代替。 我们将使用重采样函数进行转换。 为简单起见,我们还可以使用 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)
Output
co2 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()

在对数据进行绘图时,会出现一些可区分的模式。 时间序列具有明显的季节性特征,总体呈上升趋势。

要了解更多关于时间序列预处理的信息,请参阅“使用 Python 3进行时间序列可视化的指南” ,其中更详细地描述了上述步骤。

现在,我们已经转换和研究了我们的数据,让我们进入时间序列预测与 ARIMA。

第三步ー ARIMA 时间序列模型

时间序列预测中最常用的方法之一是 ARIMA 模型,意思是 a utoreg r 过剩综合 m 通过平均。 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 时间序列模型的最优参数集的过程。

第四步: ARIMA 时间序列模型的参数选择

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

我们将使用“网格搜索”迭代地探索参数的不同组合。 对于每个参数组合,我们使用状态模型中的 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]))
Output
Examples 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 (赤池信息量准则)值,这个值可以方便地返回,并使用状态空间模型拟合 ARIMA 模型。 Aic 衡量模型与数据的匹配程度,同时考虑模型的总体复杂性。 与使用较少特征以获得相同拟合优度的模型相比,在使用大量特征时非常适合数据的模型将获得较大的 AIC 得分。 因此,我们有兴趣寻找产生最低 AIC 值的模型。

下面的代码块迭代参数组合,并使用 statsmodels 中的 SARIMAX 函数来拟合相应的季节 ARIMA 模型。 这里,顺序参数指定(p,d,q)参数,而季节顺序参数指定季节 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

由于某些参数组合可能导致数值错误规格,因此我们显式禁用警告消息,以避免警告消息过载。 这些错误规范还可能导致错误并抛出异常,因此我们确保捕获这些异常并忽略导致这些问题的参数组合。

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

Output
SARIMAX(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 | 列告诉我们每个特征权重的重要性。 在这里,每个权重的 p 值都低于或接近0.05,因此在我们的模型中保留所有权重是合理的。

在拟合季节 ARIMA 模型(以及与此相关的任何其他模型)时,运行模型诊断以确保模型所做的假设没有被违反,这一点非常重要。 绘图诊断对象允许我们快速生成模型诊断并调查任何不寻常的行为。

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

我们主要关心的是确保残差的我们的模型是不相关的和正态分布的零均值。 如果季节 ARIMA 模型不满足这些性质,这是一个很好的迹象,它可以进一步改进。

在这种情况下,我们的模型诊断表明,模型残差是基于以下内容正态分布的:

  • 在右上角的图中,我们可以看到红色 KDE 线紧跟 n (0,1)线(其中 n (0,1)是平均值为0,标准差为1的正态分布的标准符号)。 这很好地表明残差是正态分布的。
  • 左下角的 qq 图显示残差(蓝点)的有序分布服从 n (0,1)标准正态分布样本的线性趋势。 再次,这是一个强烈的迹象表明残差是正态分布的。
  • 剩余时间(左上角的图)没有显示任何明显的季节性,似乎是白噪音。 右下角的自相关图(即相关图)证实了这一点,它表明时间序列的残差与其本身的滞后版本之间的相关性很低。

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

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

第六步ー验证预测

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

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

上述代码要求天气预报从1998年1月开始。

动态的 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)))
Output
The 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)))
Output
The Mean Squared Error of our forecasts is 1.01

从动态预报中获得的预测值的最小相对误差为1.01。 考虑到我们对时间序列历史数据的依赖程度较低,这个数字略高于预期的一步。

提前一步预测和动态预测都证实了该时间序列模型的有效性。 然而,围绕时间序列预测的许多兴趣是能够提前预测未来的价值。

第七步ー制作及可视化预测

在本教程的最后一步,我们描述了如何利用季节 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 模型。 我们广泛使用了熊猫和国家图书馆,并展示了如何运行模型诊断,以及如何产生二氧化碳时间序列的预测。

以下是一些你可以尝试的方法:

  • 更改动态预测的开始日期,以查看这将如何影响预测的整体质量。
  • 尝试更多的参数组合,看看是否可以改进模型的拟合优度。
  • 选择一个不同的度量标准来选择最佳模型。 例如,我们使用 AIC 测度来寻找最佳模型,但你可以寻求优化的样本外均方误差代替。

为了获得更多的实践,您还可以尝试加载另一个时间序列数据集来生成您自己的预测。