​

目录

• 引言
• Croston Classic 模型
• 加载库和数据
• 使用绘图方法探索数据
• 将数据分为训练集和测试集
• 使用 StatsForecast 实现 CrostonClassic
• 交叉验证
• 模型评估
• 参考资料

​

引言

Croston 模型是一种用于时间序列分析的方法，用于预测间歇性数据或频繁出现零值情况下的需求。它由 J.D. Croston 于 1972 年开发，在库存管理、零售销售以及低销售频率产品的需求预测等行业中特别有用。

Croston 模型基于两个主要组成部分：

1. 间歇性需求率：计算发生销售或事件期间的需求率，忽略没有销售的期间。此率用于估计将来发生需求的概率。

2. 需求间隔：计算销售或事件发生之间的时间间隔，同样忽略非销售期间。此间隔用于估计下一期间发生需求的概率。

Croston 模型结合了这两个估计值来生成加权预测，该预测同时考虑了间歇性需求率和需求之间的间隔。这种方法有助于解决时间序列具有许多零值或缺失值情况下的需求预测挑战。

需要注意的是，Croston 模型是一种简化模型，并未考虑需求数据中其他可能的变异源或模式。因此，在存在外部因素或需求行为发生变化的情况下，其准确性可能会受到影响。

​

Croston Classic 模型

​

什么是间歇性需求？

间歇性需求是一种需求模式，其特征是事件或销售的发生不规律和零星。换句话说，它指的是产品或服务需求间歇性发生的情况，其中存在没有销售或重大事件的时期。

间歇性需求与持续或规律性需求不同，后者销售随时间以可预测且一致的方式发生。相反，在间歇性需求中，没有销售的时期可能很长，并且可能没有规律的事件序列。

这种类型的需求可能发生在不同的行业和背景下，例如低消费品、季节性产品、高变异性产品、生命周期短的产品，或者需求取决于特定事件或外部因素的情况。

间歇性需求可能给预测和库存管理带来挑战，因为难以预测销售何时发生以及数量多少。像我前面提到的 Croston 模型这样的方法被用来处理间歇性需求，并为此类需求模式生成更准确和适当的预测。

​

间歇性需求的问题

间歇性需求可能在库存管理和需求预测中带来各种挑战和问题。与间歇性需求相关的一些常见问题如下：

1. 不可预测的变异性：间歇性需求可能具有不可预测的变异性，使得规划和预测困难。需求模式可能不规律，并在有销售期间和无销售期间之间大幅波动。

2. 低销售频率：间歇性需求的特点是没有销售的时期很长。这可能导致库存管理困难，因为需要持有足够的库存以在需求发生时满足需求，同时避免非销售期间的过剩库存。

3. 预测误差：预测间歇性需求比预测持续性需求更难确定。传统的预测模型可能不足以捕捉间歇性需求中的变异性和缺乏模式，这可能导致对未来需求的估计出现重大误差。

4. 对供应链的影响：间歇性需求可能影响供应链的效率，并给生产计划、供应商管理和物流带来困难。必须调整提前期和库存水平以满足不可预测的需求。

5. 运营成本：在间歇性需求情况下管理库存会增加运营成本。在非销售期间维持充足的库存以及管理库存水平可能需要额外的仓储和物流投资。

为了解决这些问题，采用了一些针对间歇性需求的特定方法，例如专门的预测模型、产品分类技术和定制的库存策略。这些解决方案旨在最小化间歇性需求中变异性和缺乏模式的影响，优化库存管理并提高供应链效率。

​

Croston 方法 (CR)

Croston 方法 (CR) 是一种专门处理间歇性需求的经典方法，它是在简单指数平滑法的基础上开发的。当 Croston 处理间歇性需求时，他发现通过使用 SES，每个期间需求的预测水平通常高于其实际值，这导致准确性非常低。经过一段时间的研究，他提出了一种优化间歇性需求预测结果的方法。

该方法基本上将间歇性需求分解为两部分：非零需求的规模和这些需求发生的时间间隔，然后对两部分应用简单指数平滑法。其公式如下：

如果 $Z_t=0$ 则

$Z'_t= Z'_{t-1}$

$P'_t= P'_{t-1}$

否则

$Z'_t=\alpha Z_t +(1-\alpha) Z'_{t-1}$

$P'_t=\alpha P_t +(1-\alpha) P'_{t-1}$

其中 $0< \alpha < 1$

最后通过结合这些预测

${Y'}_t = \frac{{Z'}_t}{{P'}_t}$

其中

• ${Y'}_t:$ 每期间平均需求。
• $Z_t:$ 期间 $t$ 的实际需求。
• $Z'_t:$ 两次正需求之间的时间。
• $P:$ 下一个期间的需求规模预测。
• $P_t:$ 需求间隔预测。
• $\alpha :$ 平滑常数。

Croston 方法将间歇性需求时间序列转换为非零需求时间序列和需求间隔时间序列，许多案例表明该方法效果很好，但在应用 Croston 方法之前，应进行三个假设：

• 非零需求是独立的且服从正态分布；
• 需求间隔是独立的且服从几何分布；
• 需求规模和需求间隔之间相互独立。

根据许多实际案例表明，Croston 方法适用于提前期服从正态分布的情况，对于包含大量零值的需求序列，Croston 方法并未表现出突出性能，有时甚至不如 SES 方法。

此外，Croston 方法只能提供每个期间的平均需求，它无法预测每个期间的需求规模，无法预测哪个期间会发生需求，也无法给出某个期间是否会发生需求的概率。

总而言之，虽然 Croston 方法是一种非常经典且广泛使用的方法，但它仍然存在许多局限性，不过经过多年统计学家和学者们的研究，也提出了一些 Croston 方法的变体。

​

Croston 的变体

Croston 方法是需求预测领域中使用的主要模型，大多数工作都基于该模型。然而，Syntetos 和 Boylan 在 2001 年提出 Croston 方法并非无偏方法，同时一些经验证据也表明使用 Croston 方法会导致性能损失（Sani 和 Kingsman，1997）。在改进 Croston 方法方面进行了大量进一步研究。Syntetos 和 Boylan (2005) 提出了一种近似无偏程序，该程序能提供更少的估计结果方差，该程序被称为 SBA（Syntetos 和 Boylan 近似）。最近，Teunter 等人 (2011) 也提出了一种可以处理过时商品的间歇性预测方法，该方法基于 Croston 方法，称为 TSB 方法（Teunter、Syntetos 和 Babai）。

​

Croston 方法的应用领域

Croston 方法通常应用于间歇性需求情况下的库存管理和需求预测领域。Croston 模型可以应用的一些具体领域有：

1. 库存管理：Croston 模型用于预测销售零星或间歇性产品的需求。有助于确定最佳库存水平和补货策略，最大程度地降低库存成本并确保充足的可供性以满足间歇性需求。

2. 零售销售：在零售业中，特别是在销售频率低或销售不规律的产品中，Croston 模型对于预测需求和优化商店或仓库的库存计划可能很有用。

3. 需求预测：一般来说，当时间序列缺乏清晰模式或具有高变异性时，Croston 模型应用于需求预测。它可用于各种行业，例如制药业、汽车业、易腐品行业以及其他间歇性需求常见的行业。

4. 供应链计划：Croston 模型可用于供应链计划和管理，以提高间歇性需求预测的准确性。这有助于精简生产、库存管理、供应商订单调度以及供应链的其他方面。

需要注意的是，Croston 模型只是解决间歇性需求的众多可用方法之一。根据时间序列的背景和具体特征，可能存在其他更适合的方法和技术。

​

用于平稳时间序列的 Croston 方法

不，Croston 方法中的时间序列不必是平稳的。Croston 方法是一种有效的间歇性时间序列预测方法，即使它们不是平稳的。但是，如果时间序列是平稳的，Croston 方法可能更准确。

Croston 方法基于这样一个想法：间歇性时间序列可以分解为两个组成部分：需求组成部分和需求之间的时间组成部分。需求组成部分使用标准时间序列预测方法（例如单指数平滑或双指数平滑）进行预测。需求之间的时间组成部分使用概率分布函数（例如泊松分布或 Weibull 分布）进行预测。

然后，Croston 方法结合这两个组成部分的预测以获得下一期间的总需求预测。

如果时间序列是平稳的，则时间序列的两个组成部分也将是平稳的。这意味着 Croston 方法将能够更准确地预测这两个组成部分。

然而，即使时间序列不平稳，Croston 方法仍然是一种有效的预测方法。Croston 方法是一种鲁棒的方法，可以处理具有不规则需求模式的时间序列。

如果您使用 Croston 方法预测非平稳的间歇性时间序列，重要的是选择一种对非平稳时间序列有效的使用标准时间序列预测方法。双指数平滑法是一种对非平稳时间序列有效的预测方法。

​

加载库和数据

提示

需要 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"])

​

使用绘图方法探索数据

使用 StatsForecast 类中的 plot 方法绘制一些序列。此方法从数据集中随机打印一个序列，对于基本探索性数据分析非常有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

自相关图

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 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. 用于训练我们的 Croston Classic 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))

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

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

​

使用 StatsForecast 实现 CrostonClassic

要了解有关 CrostonClassic Model 函数参数的更多信息，请参阅以下列表。更多信息，请访问文档。

alias : str
    Custom name of the model.

​

加载库

from statsforecast import StatsForecast
from statsforecast.models import CrostonClassic

​

实例化模型

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

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

models = [CrostonClassic()]

我们通过实例化一个新的 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=[CrostonClassic])

让我们看看我们的 Croston Classic Model 的结果。我们可以通过以下指令观察：

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

{'mean': array([27.41841685]),
 'fitted': array([     nan,  0.     ,  5.     , ..., 30.61961, 30.61961, 30.61961],
       dtype=float32),
 'sigma': np.float32(49.5709)}

​

预测方法

如果您希望在具有多个序列或模型的生产环境中提高速度，我们建议使用 StatsForecast.forecast 方法代替 .fit 和 .predict。

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

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

• h (int)： 表示预测未来 h 步。在此例中，未来 25 周。

这里的 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 方法接受两个参数：预测未来 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=) 中的性能，每隔一小时 (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)。《应用时间序列分析和 Python 预测》。
2. Ivan Svetunkov。使用增强型动态自适应模型 (ADAM) 进行预测和分析
3. James D. Hamilton。《时间序列分析》普林斯顿大学出版社，普林斯顿，新泽西州，第 1 版，1994 年。
4. Nixtla 参数.
5. Pandas 可用的频率.
6. Rob J. Hyndman 和 George Athanasopoulos (2018)。《预测原理与实践（第 3 版）》.
7. 季节周期 - Rob J Hyndman.