​

目录

• 引言
• Croston 优化模型
• 加载库和数据
• 使用 plot 方法探索数据
• 将数据分为训练集和测试集
• 使用 StatsForecast 实现 CrostonOptimized
• 交叉验证
• 模型评估
• 参考文献

​

引言

Croston 优化模型是一种专为间歇性需求时间序列数据设计的预测方法。它是 Croston 方法的扩展，后者最初是为了预测零星需求模式而开发的。

间歇性需求时间序列的特点是，非零需求值出现不规律且零星，通常伴随着长时间的零需求期。传统的预测方法可能难以有效地处理此类模式。

Croston 优化模型通过整合两个关键组成部分来解决这一挑战：指数平滑和间歇性需求估计。

1. 指数平滑：Croston 优化模型使用指数平滑来捕捉间歇性需求数据中的趋势和季节性。这有助于识别潜在模式并进行更准确的预测。

2. 间歇性需求估计：由于间歇性需求数据通常包含长时间的零需求期，Croston 优化模型对非零需求值的出现和大小采用单独的估计过程。它估计出现概率和非零需求区间的平均大小，从而能够更好地预测间歇性需求。

Croston 优化模型旨在在间歇性需求的过度预测和预测不足之间取得平衡，这是传统预测方法中的常见挑战。通过显式建模间歇性需求模式，它可以为间歇性需求时间序列数据提供更准确的预测。

值得注意的是，Croston 优化模型存在多种变体和适应版本，为适应特定的预测场景进行了不同的修改和增强。这些变体可能包含额外的特征或算法，以进一步提高预测的准确性。

​

Croston 优化方法

Croston 优化模型可以数学定义如下

1. 初始化
   • 令 $(y_t)$ 表示时间 $t$ 的间歇性需求时间序列数据。
   • 初始化两组变量：$(p_t)$ 用于表示出现概率，$(q_t)$ 用于表示非零需求区间的平均大小。
   • 将预测值 $(F_t)$ 和预测误差 $(E_t)$ 变量初始化为零。
2. 计算 $(p_t)$ 和 $(q_t)$
   • 使用指数平滑计算间歇性需求出现概率 $(p_t)$：$[p_t = \alpha + (1 - \alpha)(p_{t-1}),]$ 其中 $(\alpha)$ 是平滑参数（通常设置在 0.1 到 0.3 之间）。
   • 使用指数平滑计算非零需求区间的平均大小 $(q_t)$：$[q_t = \beta \cdot y_t + (1 - \beta)(q_{t-1}),]$ 其中 $(\beta)$ 是平滑参数（通常设置在 0.1 到 0.3 之间）。
3. 预测
   • 如果 $(y_t > 0)$ （非零需求发生）
     • 计算预测值 $(F_t)$ 为前一个预测值 $(F_{t-1})$ 除以非零需求区间的平均大小 $(q_{t-1})$：$[F_t = \frac{{F_{t-1}}}{{q_{t-1}}}]$
     • 计算预测误差 $(E_t)$ 为实际需求 $(y_t)$ 与预测值 $(F_t)$ 之间的差值：$[E_t = y_t - F_t]$
   • 如果 $(y_t = 0)$ （零需求发生）
     • 将预测值 $(F_t)$ 和预测误差 $(E_t)$ 设置为零。
4. 更新模型
   • 使用指数平滑更新间歇性需求出现概率 $(p_t)$ 和非零需求区间的平均大小 $(q_t)$，如步骤 2 所述。
5. 对时间序列中的每个时间点重复步骤 3 和 4。

Croston 优化模型利用指数平滑捕捉间歇性需求数据中的趋势和季节性，并分别估计非零需求区间的出现概率和平均大小，以有效处理间歇性需求模式。通过根据观测数据更新模型参数，它为间歇性需求时间序列提供了改进的预测。

​

优化 Croston 模型的一些特性

优化 Croston 模型是经典 Croston 模型的修改版本，用于预测间歇性需求。经典 Croston 模型使用历史订单的加权平均值和订单之间的平均间隔来预测需求。优化 Croston 模型使用概率函数来预测订单之间的平均间隔。

优化 Croston 模型已被证明对于具有不规律需求的时间序列比经典 Croston 模型更准确。优化 Croston 模型也更能适应不同类型的间歇性时间序列。

优化 Croston 模型具有以下特性

• 即使对于具有不规律需求的时间序列，它也很准确。
• 它能适应不同类型的间歇性时间序列。
• 易于实现和理解。
• 对异常值具有鲁棒性。

优化 Croston 模型已成功用于预测各种间歇性时间序列，包括产品需求、服务需求和资源需求。

以下是优化 Croston 模型的一些特性

• 精度：优化 Croston 模型已被证明对于具有不规律需求的时间序列比经典 Croston 模型更准确。这是因为优化 Croston 模型使用概率函数来预测订单之间的平均间隔，这比历史订单的加权平均值更准确。
• 适应性：优化 Croston 模型也更能适应不同类型的间歇性时间序列。这是因为优化 Croston 模型使用概率函数来预测订单之间的平均间隔，使其能够适应不同的需求模式。
• 易于实现和理解：优化 Croston 模型易于实现和理解。这是因为优化 Croston 模型是经典 Croston 模型的修改版本，而后者是一个众所周知且易于理解的模型。
• 鲁棒性：优化 Croston 模型对异常值也具有鲁棒性。这是因为优化 Croston 模型使用概率函数来预测订单之间的平均间隔，这使得它能够忽略异常值。

​

加载库和数据

提示

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

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

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import plotly.graph_objects as go
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");

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

plt.show();

​

时间序列分解

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

在时间序列分析中预测新值时，了解过去的数据非常重要。更正式地说，我们可以说了解值随时间变化的模式非常重要。预测值出现偏差可能有很多原因。基本上，时间序列由四个组成部分构成。这些组成部分的变动导致时间序列模式的变化。这些组成部分是

• 水平项 (Level)： 这是随时间变化的平均主要值。
• 趋势项 (Trend)： 趋势是导致时间序列出现递增或递减模式的值。
• 季节性项 (Seasonality)： 这是时间序列中短暂出现的周期性事件，导致时间序列的短期递增或递减模式。
• 残差/噪声 (Residual/Noise)： 这些是时间序列中的随机变动。

这些组成部分随时间的结合形成时间序列。大多数时间序列包含水平项和噪声/残差，而趋势项或季节性项是可选值。

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

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

​

加法时间序列

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

​

乘法时间序列

如果时间序列的组成部分相乘在一起，则该时间序列称为乘法时间序列。通过可视化，如果时间序列随时间呈现指数增长或下降，则该时间序列可视为乘法时间序列。乘法时间序列的数学函数可以表示为。

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

from plotly.subplots import make_subplots

from statsmodels.tsa.seasonal import seasonal_decompose

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. 用于训练 Croston Optimized Model 的数据。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))

​

使用 StatsForecast 实现 CrostonOptimized

​

加载库

from statsforecast import StatsForecast
from statsforecast.models import CrostonOptimized

​

实例化模型

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

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

# We call the model that we are going to use
models = [CrostonOptimized()]

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

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

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

• n_jobs: n_jobs: int，并行处理中使用的任务数，使用 -1 表示所有核心。

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

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

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

​

拟合模型

# fit the models
sf.fit(df=train)

StatsForecast(models=[CrostonOptimized])

让我们看看 Croston optimized Model 的结果。我们可以使用以下说明来观察它

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

{'mean': array([27.41841685])}

​

预测方法

如果您想在具有多个序列或模型的生产环境中提高速度，建议使用 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, max_insample_length=500)

​

带置信区间的 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: 序列标识符
• 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 系列 (2022)。Applied Time Series Analysis and Forecasting with Python。
2. Ivan Svetunkov. 使用增广动态自适应模型 (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 和 George Athanasopoulos (2018)。“预测原理与实践（第 3 版）”.
7. 季节周期 - Rob J Hyndman.