​

目录

• StatsForecast 中的 AutoArima 是什么？
• Arima 模型的定义
• 使用 AutoArima 的优势
• 加载库和数据
• 使用 plot 方法探索数据
• 将数据分割为训练集和测试集
• 使用 StatsForecast 实现 AutoARIMA
• 交叉验证
• 模型评估
• 参考资料

​

StatsForecast 中的 AutoArima 是什么？

AutoARIMA 是一种时间序列模型，它使用自动化过程为给定的时间序列选择最佳 ARIMA（自回归积分滑动平均）模型参数。ARIMA 是一种广泛用于时间序列建模和预测的统计模型。

AutoARIMA 模型中的自动参数选择过程采用统计和优化技术（例如 Akaike 信息准则 (AIC) 和交叉验证），以确定 ARIMA 模型的自回归、积分和滑动平均参数的最佳值。

自动参数选择很有用，因为如果没有对生成时间序列的基础随机过程的深入理解，很难确定给定时间序列的 ARIMA 模型最佳参数。autoARIMA 模型自动化了参数选择过程，可以为时间序列建模和预测提供快速有效的解决方案。

statsforecast.models 库提供了 Python 的 AutoARIMA 函数实现，可以自动为给定时间序列选择 ARIMA 模型的最佳参数。

​

Arima 模型的定义

ARIMA 模型（自回归积分滑动平均）过程是自回归过程 AR(p)、积分 I(d) 和滑动平均过程 MA(q) 的组合。

与 ARMA 过程一样，ARIMA 过程指出，当前值取决于过去的值（来自 AR(p) 部分）和过去的误差（来自 MA(q) 部分）。然而，ARIMA 过程不使用表示为 yt 的原序列，而是使用表示为$y'_{t}$ 的差分序列。请注意，$y'_{t}$ 可以表示经过多次差分的序列。

因此，ARIMA(p,d,q) 过程的数学表达式指出，差分序列$y'_{t}$ 的当前值等于一个常数 $C$、差分序列的过去值 $\phi_{p}y'_{t-p}$、差分序列的均值 $\mu$、过去的误差项 $\theta_{q}\varepsilon_{t-q}$ 以及当前的误差项 $\varepsilon_{t}$，如方程所示

其中 $y'_{t}$ 是差分序列（可能经过多次差分）。等号右侧的“预测变量”包括 $y_{t}$ 的滞后值和滞后误差。我们称之为 ARIMA(p,d,q) 模型，其中

用于自回归和滑动平均模型的平稳性和可逆性条件同样适用于 ARIMA 模型。

我们已经讨论过的许多模型都是 ARIMA 模型的特例，如表所示

一旦我们开始以这种方式组合组件以形成更复杂的模型，使用滞后算子（backshift notation）会容易得多。例如，公式 (1) 可以用滞后算子表示为

选择合适的 p、d 和 q 值可能很困难。然而，statsforecast 库中的 AutoARIMA() 函数将自动为您完成此操作。

更多信息请参阅此处

​

加载库和数据

使用 AutoARIMA() 模型进行时间序列建模和预测有几个优点，包括

1. 参数选择过程自动化：AutoARIMA() 函数自动化了 ARIMA 模型参数选择过程，通过消除手动尝试不同参数组合的需要，为用户节省时间和精力。

2. 降低预测误差：通过自动选择最优参数，与手动选择的 ARIMA 模型相比，ARIMA 模型可以提高预测的准确性。

3. 识别复杂模式：AutoARIMA() 函数可以识别数据中可能难以通过可视化或其他时间序列建模技术检测到的复杂模式。

4. 参数选择方法的灵活性：ARIMA 模型可以使用不同的方法来选择最优参数，例如赤池信息准则（AIC）、交叉验证等，这使用户能够选择最适合其需求的方法。

通常，使用 AutoARIMA() 函数有助于提高时间序列建模和预测的效率和准确性，特别是对于不熟悉手动选择 ARIMA 模型参数的用户。

​

主要结果

我们比较了与 pmdarima、Rob Hyndman 的 forecast 包以及 Facebook 的 Prophet 在准确性和速度上的差异。我们使用了 M4 竞赛中的 Daily（日度）、Hourly（小时）和 Weekly（周度）数据。

下表总结了结果。可以看出，我们的 auto_arima 在准确性（由 MASE 损失衡量）和时间方面是最好的模型，甚至优于 R 中的原始实现。

[1] pmdarima 的 auto_arima 模型在小时数据上存在问题。已开启一个 issue。

下表总结了数据详情。

​

加载库和数据

提示

需要安装 Statsforecast。安装方法请参阅说明。

接下来，我们导入绘图库并配置绘图样式。

import numpy as np
import pandas as pd

import scipy.stats as stats

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#212946',
    'axes.facecolor': '#212946',
    'savefig.facecolor':'#212946',
    'axes.grid': True,
    'axes.grid.which': 'both',
    'axes.spines.left': False,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.spines.bottom': False,
    'grid.color': '#2A3459',
    'grid.linewidth': '1',
    'text.color': '0.9',
    'axes.labelcolor': '0.9',
    'xtick.color': '0.9',
    'ytick.color': '0.9',
    'font.size': 12 }
plt.rcParams.update(dark_style)

from pylab import rcParams
rcParams['figure.figsize'] = (18,7)

​

加载数据

df = pd.read_csv("https://raw.githubusercontent.com/Naren8520/Serie-de-tiempo-con-Machine-Learning/main/Data/candy_production.csv")
df.head()

StatsForecast 的输入始终是一个长格式的数据帧，包含三列：unique_id、ds 和 y

• unique_id（字符串、整数或类别）代表时间序列的标识符。

• ds（日期戳）列应采用 Pandas 期望的格式，日期通常为 YYYY-MM-DD，时间戳通常为 YYYY-MM-DD HH:MM:SS。

• y（数值）代表我们希望预测的度量值。

df["unique_id"]="1"
df.columns=["ds", "y", "unique_id"]
df.head()

print(df.dtypes)

ds            object
y            float64
unique_id     object
dtype: object

我们需要将 ds 从 object 类型转换为 datetime 类型。

df["ds"] = pd.to_datetime(df["ds"])

​

使用 plot 方法探索数据

使用 StatsForecast 类中的 plot 方法绘制一个时间序列。此方法会从数据集中打印一个随机系列，对进行基本探索性数据分析 (EDA) 非常有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

自相关图

fig, axs = plt.subplots(nrows=1, ncols=2)

plot_acf(df["y"],  lags=60, ax=axs[0],color="fuchsia")
axs[0].set_title("Autocorrelation");

plot_pacf(df["y"],  lags=60, ax=axs[1],color="lime")
axs[1].set_title('Partial Autocorrelation')

plt.show();

​

时间序列分解

如何以及为何分解时间序列？

在时间序列分析中，为了预测新值，了解过去的数据非常重要。更正式地说，了解数值随时间变化的模式非常重要。许多原因可能导致我们的预测值偏离正确的方向。基本上，时间序列由四个组成部分组成。这些组成部分的变化导致时间序列模式的改变。这些组成部分是：

• 水平项： 这是随时间平均的主要值。
• 趋势项： 趋势是导致时间序列呈现上升或下降模式的值。
• 季节项： 这是时间序列中周期性发生的事件，持续时间较短，导致时间序列呈现短期上升或下降模式。
• 残差/噪声项： 这些是时间序列中的随机变化。

随着时间的推移，这些组成部分的组合形成了时间序列。大多数时间序列由水平项和噪声/残差项组成，而趋势项或季节项是可选的。

如果季节项和趋势项是时间序列的一部分，那么会对预测值产生影响。因为预测时间序列的模式可能与之前的时间序列不同。

时间序列中组成部分的组合可以是两种类型： * 加法模型 * 乘法模型

加法时间序列模型

如果时间序列的组成部分相加形成时间序列。则称该时间序列为加法时间序列。通过可视化，我们可以说如果时间序列的上升或下降模式在整个序列中相似，则它是加法模型。任何加法时间序列的数学函数可以表示为：$y(t) = level + Trend + seasonality + noise$

​

乘法时间序列模型

如果时间序列的组成部分相乘形成时间序列，则称该时间序列为乘法时间序列。从可视化角度看，如果时间序列随时间呈现指数增长或下降，则可以将其视为乘法时间序列。乘法时间序列的数学函数可以表示为。

$y(t) = Level * Trend * seasonality * Noise$

from statsmodels.tsa.seasonal import seasonal_decompose 
a = seasonal_decompose(df["y"], model = "add", period=12)
a.plot();

​

将数据分割为训练集和测试集

我们将数据分为两部分：1. 用于训练 AutoArima 模型的数据 2. 用于测试模型的数据

对于测试数据，我们将使用最后 12 个月的数据来测试和评估模型的性能。

Y_train_df = df[df.ds<='2016-08-01'] 
Y_test_df = df[df.ds>'2016-08-01']

Y_train_df.shape, Y_test_df.shape

((536, 3), (12, 3))

现在让我们绘制训练数据和测试数据。

sns.lineplot(Y_train_df,x="ds", y="y", label="Train")
sns.lineplot(Y_test_df, x="ds", y="y", label="Test")
plt.show()

​

使用 StatsForecast 实现 AutoArima

​

加载库

from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA
from statsforecast.arima import arima_string

​

实例化模型

导入并实例化模型。设置 season_length 参数有时很棘手。大师 Rob Hyndmann 关于季节周期的文章可能很有帮助。

season_length = 12 # Monthly data 
horizon = len(Y_test_df) # number of predictions

models = [AutoARIMA(season_length=season_length)]

我们通过实例化一个新的 StatsForecast 对象并使用以下参数来拟合模型

models：模型列表。从模型中选择您想要的模型并导入它们。

• freq： 一个字符串，表示数据的频率。（请参阅 panda 可用的频率。）

• n_jobs： n_jobs：整数，并行处理中使用的作业数，使用 -1 表示所有核心。

• fallback_model： 如果某个模型失败，则使用的备用模型。

任何设置都通过构造函数传递。然后调用其 fit 方法并传入历史数据帧。

sf = StatsForecast(models=models, freq='MS')

​

拟合模型

sf.fit(df=Y_train_df)

StatsForecast(models=[AutoARIMA])

输入模型后，我们可以使用arima_string 函数查看模型找到的参数。

arima_string(sf.fitted_[0,0].model_)

'ARIMA(4,0,3)(0,1,1)[12]                   '

自动化过程告诉我们找到的最佳模型形式为 ARIMA(4,0,3)(0,1,1)[12]，这意味着我们的模型包含 $p=4$，即它包含一个非季节性自回归成分；另一方面，我们的模型包含一个季节性部分，其阶数为 $D=1$，即它包含一个季节性差分；以及 $q=3$，它包含 3 个移动平均成分。

为了了解模型项的值，我们可以使用以下语句来查看模型的全部结果。

result=sf.fitted_[0,0].model_
print(result.keys())
print(result['arma'])

dict_keys(['coef', 'sigma2', 'var_coef', 'mask', 'loglik', 'aic', 'arma', 'residuals', 'code', 'n_cond', 'nobs', 'model', 'bic', 'aicc', 'ic', 'xreg', 'x', 'lambda'])
(4, 3, 0, 1, 12, 0, 1)

现在让我们可视化模型的残差。

如我们所见，上面得到的结果是一个字典输出，要从字典中提取每个元素，我们将使用 .get() 函数提取元素，然后将其保存在 pd.DataFrame() 中。

residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual

fig, axs = plt.subplots(nrows=2, ncols=2)

# plot[1,1]
residual.plot(ax=axs[0,0])
axs[0,0].set_title("Residuals");

# plot
sns.distplot(residual, ax=axs[0,1]);
axs[0,1].set_title("Density plot - Residual");

# plot
stats.probplot(residual["residual Model"], dist="norm", plot=axs[1,0])
axs[1,0].set_title('Plot Q-Q')

# plot
plot_acf(residual,  lags=35, ax=axs[1,1],color="fuchsia")
axs[1,1].set_title("Autocorrelation");

plt.show();

要生成预测，我们只需使用 predict 方法并指定预测范围 (h)。此外，为了计算与预测相关的预测区间，我们可以包含 level 参数，该参数接收一个列表，其中包含我们想要构建的预测区间的置信水平。在本例中，我们将只计算 90% 的预测区间 (level=[90])。

​

预测方法

如果您想在具有多个时间序列或模型的生产环境中提高速度，我们建议使用 StatsForecast.forecast 方法，而不是 .fit 和 .predict。

主要区别在于 .forecast 不存储拟合值，并且在分布式环境中具有高度可扩展性。

forecast 方法接受两个参数：预测未来 h（范围）和 level。

• h（整数）： 表示未来 h 步的预测。在本例中，即未来 12 个月。

• level（浮点数列表）： 此可选参数用于概率预测。设置预测区间的置信水平（或置信百分位数）。例如，level=[90] 表示模型期望真实值有 90% 的时间落在此区间内。

这里的 forecast 对象是一个新的数据帧，包含模型名称和 y 帽值（预测值）的列，以及不确定性区间的列。根据您的计算机性能，此步骤大约需要 1 分钟。（如果您想将速度提高到几秒钟，请移除 AutoModels，例如 ARIMA 和 Theta）

Y_hat_df = sf.forecast(df=Y_train_df, h=horizon, fitted=True)
Y_hat_df.head()

values=sf.forecast_fitted_values()
values

使用 forecast 方法添加 95% 置信区间

sf.forecast(df=Y_train_df, h=12, level=[95])

Y_hat_df = Y_test_df.merge(Y_hat_df, how='left', on=['unique_id', 'ds'])

fig, ax = plt.subplots(1, 1, figsize = (18, 7))
plot_df = pd.concat([Y_train_df, Y_hat_df]).set_index('ds')
plot_df[['y', 'AutoARIMA']].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Monthly ', fontsize=20)
ax.set_xlabel('Timestamp [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid()

​

带有置信区间的 Predict 方法

要生成预测，请使用 predict 方法。

predict 方法接受两个参数：预测未来 h（范围）和 level。

• h（整数）： 表示未来 h 步的预测。在本例中，即未来 12 个月。

• level（浮点数列表）： 此可选参数用于概率预测。设置预测区间的置信水平（或置信百分位数）。例如，level=[95] 表示模型期望真实值有 95% 的时间落在此区间内。

这里的 forecast 对象是一个新的数据帧，包含模型名称和 y 帽值（预测值）的列，以及不确定性区间的列。

此步骤应少于 1 秒。

sf.predict(h=12)

forecast_df = sf.predict(h=12, level = [80, 95]) 
forecast_df

我们可以使用 pandas 函数 pd.concat() 将预测结果与历史数据连接起来，然后将此结果用于绘图。

df_plot=pd.concat([df, forecast_df]).set_index('ds').tail(220)
df_plot

现在让我们可视化预测结果以及时间序列的历史数据，同时绘制我们在进行 95% 置信度预测时获得的置信区间。

sf.plot(df, forecast_df, level=[95], max_insample_length=12 * 5)

​

交叉验证

在之前的步骤中，我们使用了历史数据来预测未来。然而，为了评估其准确性，我们还想了解模型在过去会如何表现。要评估您的模型在数据上的准确性和稳健性，请执行交叉验证。

对于时间序列数据，交叉验证是通过在历史数据上定义一个滑动窗口并预测其后的一段时间来完成的。这种形式的交叉验证使我们能够更准确地估计模型在更广泛的时间点上的预测能力，同时保持训练集中的数据连续性，这是模型所要求的。

下图描述了这种交叉验证策略

​

执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现都非常慢。Statsforecast 库将交叉验证实现为分布式操作，从而减少执行过程所需的时间。如果数据集很大，您还可以使用 Ray、Dask 或 Spark 在分布式集群中执行交叉验证。

在本例中，我们希望评估每个模型在过去 5 个月中的性能 (n_windows=5)，每隔 12 个月预测一次 (step_size=12)。根据您的计算机性能，此步骤大约需要 1 分钟。

StatsForecast 类中的 cross_validation 方法接受以下参数。

• df： 训练数据帧

• h（整数）： 表示预测未来的 h 步。在本例中，即未来 12 个月。

• step_size（整数）： 每个窗口之间的步长。换句话说：您希望多久运行一次预测过程。

• n_windows（整数）： 用于交叉验证的窗口数量。换句话说：您希望评估过去多少个预测过程。

crossvalidation_df = sf.cross_validation(df=Y_train_df,
                                         h=12,
                                         step_size=12,
                                         n_windows=5)

crossvaldation_df 对象是一个新的数据帧，包含以下列

• unique_id： 时间序列标识符
• ds： 日期戳或时间索引
• cutoff： n_windows 的最后一个日期戳或时间索引。
• y： 真实值
• "model"： 包含模型名称和拟合值的列。

crossvalidation_df.head()

​

模型评估

现在我们将使用预测结果来评估我们的模型，我们将使用不同类型的指标，如 MAE、MAPE、MASE、RMSE、SMAPE 来评估准确性。

from functools import partial

import utilsforecast.losses as ufl
from utilsforecast.evaluation import evaluate

evaluate(
    Y_test_df.merge(Y_hat_df),
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=season_length), ufl.rmse, ufl.smape],
    train_df=Y_train_df,
)

​

参考资料

1. Nixtla-Arima
2. Rob J. Hyndman and George Athanasopoulos (2018)。《预测原则与实践（第三版）》.