IMAPA 模型

​

目录

• 简介
• IMAPA 模型
• 加载库和数据
• 使用 plot 方法探索数据
• 将数据拆分为训练集和测试集
• 使用 StatsForecast 实现 IMAPA
• 交叉验证
• 模型评估
• 参考资料

​

简介

IMAPA 是一种算法，它使用多个模型来预测间歇性时间序列的未来值。该算法首先以固定的时间间隔对时间序列值进行聚合。然后使用预测模型来预测聚合值。

IMAPA 是间歇性时间序列的良好选择，因为它对缺失值具有鲁棒性且计算效率高。IMAPA 也易于实现。

IMAPA 已在各种间歇性时间序列上进行测试，并被证明在预测未来值方面是有效的。

​

IMAPA 方法

间歇性多重聚合预测算法 (IMAPA) 模型是一种时间序列模型，用于预测间歇性时间序列的未来值。IMAPA 模型基于在固定时间间隔内聚合时间序列值，然后使用预测模型来预测聚合值的思想。可以使用任何预测模型来预测聚合值。它使用优化的 SES 在新的水平上生成预测，然后使用简单平均值将它们组合起来。

IMAPA 模型可以用数学形式定义如下

$\hat{y}_{t+1} = f(\hat{y}_{t-\tau}, \hat{y}_{t-2\tau}, ..., \hat{ y}_{t-m\tau})$

其中 $\hat{y}_{t+1}$ 是预测时间值 $t+1$，$f$ 是预测模型，$\hat{y}_{t-\tau} , \hat{y}_{t-2\tau}, ..., \hat{y}_{t-m\tau}$ 是在时间点 $t-\tau, t-2 \tau, ..., t-m\tau$ 聚合值的预测，而 $\tau$ 是对时间序列值进行聚合的时间间隔。

IMAPA 是间歇性时间序列的良好选择，因为它对缺失值具有鲁棒性且计算效率高。IMAPA 也易于实现。

IMAPA 已在各种间歇性时间序列上进行测试，并被证明在预测未来值方面是有效的。

​

IMAPA 一般属性

• 多重聚合：IMAPA 使用多重聚合级别来分析和预测间歇性时间序列。这涉及到将原始序列分解为不同时间尺度的组成部分。

• 间歇性：IMAPA 专注于处理间歇性时间序列，这类序列表现出不规则和非平稳模式，有活动期和非活动期。

• 自适应预测：IMAPA 采用自适应方法根据收集到的新数据调整预测模型。这使得算法能够随时间变化适应时间序列行为的变化。

• 对缺失值具有鲁棒性：IMAPA 可以处理数据中的缺失值而不会牺牲准确性。这对于经常有缺失值的间歇性时间序列很重要。

• 计算效率高：IMAPA 计算效率高，这意味着它可以快速预测未来值。这对于大型时间序列很重要，使用其他方法可能需要很长时间进行预测。

• 分解属性：时间序列可以分解为趋势、季节性和残差等组成部分。

​

加载库和数据

提示

需要 Statsforecast。要安装，请参阅说明。

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

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.style.use('grayscale') # fivethirtyeight  grayscale  classic
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#008080',  # #212946
    'axes.facecolor': '#008080',
    'savefig.facecolor': '#008080',
    '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': '#000000',  #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)

import pandas as pd

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

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             int64
unique_id    object
dtype: object

我们可以看到时间变量 (ds) 是对象格式，需要将其转换为日期格式

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

​

使用 plot 方法探索数据

使用 StatsForecast 类中的 plot 方法绘制一些时间序列。此方法会打印数据集中的一个随机时间序列，对基本 EDA（探索性数据分析）很有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

自相关图

自相关 (ACF) 和偏自相关 (PACF) 图是用于分析时间序列的统计工具。ACF 图显示时间序列值与其滞后值之间的相关性，而 PACF 图显示时间序列值与其滞后值之间的相关性，同时去除了先前滞后值的影响。

ACF 和 PACF 图可用于识别时间序列的结构，这有助于选择合适的时间序列模型。例如，如果 ACF 图显示重复的峰谷模式，则表示时间序列是平稳的，这意味着其统计属性随时间保持不变。如果 PACF 图显示快速下降的尖峰模式，则表示时间序列是可逆的，这意味着可以对其进行逆转以获得平稳时间序列。

ACF 和 PACF 图的重要性在于它们可以帮助分析师更好地理解时间序列的结构。这种理解有助于选择合适的时间序列模型，从而提高预测时间序列未来值的能力。

分析 ACF 和 PACF 图

• 查找图中的模式。常见模式包括重复的峰谷、锯齿模式和平台模式。
• 比较 ACF 和 PACF 图。PACF 图通常比 ACF 图的尖峰少。
• 考虑时间序列的长度。较长时间序列的 ACF 和 PACF 图将有更多的尖峰。
• 使用置信区间。ACF 和 PACF 图还显示自相关值的置信区间。如果自相关值超出置信区间，则可能具有统计显著性。

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

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

# Grafico
plot_pacf(df["y"],  lags=30, 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
from plotly.subplots import make_subplots
import plotly.graph_objects as go

def plotSeasonalDecompose(
    x,
    model='additive',
    filt=None,
    period=None,
    two_sided=True,
    extrapolate_trend=0,
    title="Seasonal Decomposition"):

    result = seasonal_decompose(
            x, model=model, filt=filt, period=period,
            two_sided=two_sided, extrapolate_trend=extrapolate_trend)
    fig = make_subplots(
            rows=4, cols=1,
            subplot_titles=["Observed", "Trend", "Seasonal", "Residuals"])
    for idx, col in enumerate(['observed', 'trend', 'seasonal', 'resid']):
        fig.add_trace(
            go.Scatter(x=result.observed.index, y=getattr(result, col), mode='lines'),
                row=idx+1, col=1,
            )
    return fig

plotSeasonalDecompose(
    df["y"],
    model="additive",
    period=24,
    title="Seasonal Decomposition")

​

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

让我们将数据分成集合

1. 用于训练 IMAPA 模型 的数据。
2. 用于测试模型的数据

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

train = df[df.ds<='2023-01-31 19:00:00'] 
test = df[df.ds>'2023-01-31 19:00:00']

train.shape, test.shape

((9500, 3), (500, 3))

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

sns.lineplot(train,x="ds", y="y", label="Train", linestyle="--",linewidth=2)
sns.lineplot(test, x="ds", y="y", label="Test", linewidth=2, color="yellow")
plt.title("Store visit");
plt.xlabel("Hours")
plt.show()

​

使用 StatsForecast 实现 IMAPA 方法

​

加载库

from statsforecast import StatsForecast
from statsforecast.models import IMAPA

​

实例化模型

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

season_length = 24 # Hourly data 
horizon = len(test) # number of predictions

models = [IMAPA()]

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

models：模型列表。从 models 中选择所需的模型并导入。

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

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

• fallback_model: 如果模型失败，将使用的备用模型。

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

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

​

拟合模型

sf.fit(df=train)

StatsForecast(models=[IMAPA])

让我们看看 IMAPA 模型 的结果。我们可以通过以下说明进行观察

result=sf.fitted_[0,0].model_
result

{'mean': array([28.579695], dtype=float32)}

​

预测方法

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

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

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

• h (int): 表示预测未来 h 步长。在本例中，即未来 500 小时。

这里的 forecast 对象是一个新的数据帧，包含一列模型名称和 y hat 值，以及不确定性区间列。根据您的计算机性能，此步骤大约需要 1 分钟。

Y_hat = sf.forecast(df=train, h=horizon)
Y_hat

sf.plot(train, Y_hat)

​

带有置信区间的 predict 方法

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

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

• h (int): 表示预测未来 h 步长。在本例中，即未来 500 小时。

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

此步骤应少于 1 秒。

forecast_df = sf.predict(h=horizon) 
forecast_df

​

交叉验证

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

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

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

​

执行时间序列交叉验证

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

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

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

• df: 训练数据帧

• h (int): 表示预测未来 h 步长。在本例中，即未来 500 小时。

• step_size (int): 每个窗口之间的步长。换句话说：您希望多久运行一次预测过程。

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

crossvalidation_df = sf.cross_validation(df=df,
                                         h=horizon,
                                         step_size=50,
                                         n_windows=5)

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

• unique_id: 索引。如果您不喜欢使用索引，只需运行 crossvalidation_df.resetindex()。
• ds: 日期戳或时间索引
• cutoff: n_windows 的最后一个日期戳或时间索引。
• y: 真实值
• model: 包含模型名称和拟合值的列。

crossvalidation_df

​

模型评估

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

from functools import partial

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

evaluate(
    test.merge(Y_hat),
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=season_length), ufl.rmse, ufl.smape],
    train_df=train,
)

​

参考资料

1. Changquan Huang • Alla Petukhina. Springer series (2022). Applied Time Series Analysis and Forecasting with Python.
2. Ivan Svetunkov. Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM)
3. James D. Hamilton. Time Series Analysis Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.
4. Nixtla 参数.
5. Pandas 可用频率.
6. Rob J. Hyndman and George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”.
7. Seasonal periods - Rob J Hyndman.