​

目录

• 引言
• Croston SBA 模型
• 加载库和数据
• 使用 plot 方法探索数据
• 将数据分割为训练集和测试集
• 在 StatsForecast 中实现 CrostonSBA
• 交叉验证
• 模型评估
• 参考资料

​

引言

Croston 模型是一种用于预测间歇性需求时间序列数据的方法，即在许多时期需求为零，只有少数时期存在非零需求的数据。Croston 的方法最初由 J.D. Croston 于 1972 年提出。随后，Syntetos 和 Boylan 在 2001 年对原始模型进行了改进，称为 Croston-SBA（Syntetos 和 Boylan 近似）。

Croston-SBA 模型基于间歇性需求遵循二项式过程的假设。模型的重点不是直接对需求进行建模，而是对需求期之间的间隔进行建模。该模型有两个主要组成部分：一个用于对需求期之间的间隔进行建模（假设遵循泊松分布），另一个用于对发生需求时的需求量进行建模。

值得注意的是，Croston-SBA 模型假设非零需求期之间的间隔是独立的，并遵循泊松分布。然而，该模型是一种近似方法，可能不适用于所有情况。建议在实际使用前评估其在历史数据上的表现。

​

Croston SBA 模型

SBA 的公式与原始 Croston 方法非常相似，但它应用了一个校正因子，减少了最终估计结果中的误差。

如果 $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}\ where \ 0<\alpha < 1$

$Y'_t=(1-\frac{\alpha}{2}) \frac{Z'_t}{P'_t}$

其中

• $Y'_t:$ 平均每期需求量
• $Z_t:$ 第 $t$ 期的实际需求量
• $Z'_t:$ 两次正需求之间的时间
• $P:$ 下一期的需求量预测
• $P'_t:$ 需求间隔预测
• $\alpha:$ 平滑常数

注意：在 Croston 方法中，结果通常会呈现相当大的正偏差，而在 SBA 中，偏差会减少，有时甚至会出现轻微的负偏差。

​

Croston SBA 方法的原理

Croston SBA（Syntetos 和 Boylan 近似）方法是一种用于预测具有间歇性或零星数据的时间序列的技术。该方法基于原始 Croston 方法，该方法旨在预测数据稀疏或未按固定间隔提供时的库存需求。

Croston SBA 方法的主要特性如下

1. 适用于间歇性数据：当数据呈现间歇性模式时，即需求期后紧跟着非需求期，Croston SBA 方法特别有用。Croston SBA 方法不将非需求期的数据视为零，而是估计需求发生率和条件需求率。

2. 频率和水平分离：Croston SBA 方法的一个关键特性是它将需求数据中的频率和水平信息分离开来。这使得可以分别对这两个组件进行建模和预测，从而获得更好的预测结果。

3. 发生率和需求率估计：Croston SBA 方法使用简单的指数平滑技术来估计条件发生率和需求率。然后利用这些比率来预测未来需求。

4. 不假设数据分布：与一些假定特定数据分布的预测技术不同，Croston SBA 方法不对需求分布做任何假设。这使其更加灵活，适用于各种情况。

5. 不需要完整的历史数据：即使历史数据稀疏或未按固定间隔提供，Croston SBA 方法也能奏效。这使得它成为使用有限数据预测间歇性需求的一个有吸引力的选择。

需要注意的是，Croston SBA 方法是一种近似方法，可能不适用于所有情况。建议将其性能与其他预测技术一起进行评估，并根据数据的具体特征和问题的背景进行调整。

在 Croston SBA 方法中，数据系列无需平稳。Croston SBA 方法适用于预测具有间歇性数据的时间序列，其中需求期穿插着非需求期。

Croston SBA 方法基于使用简单的指数平滑技术估计发生率和条件需求率。这些比率用于预测未来需求。

在时间序列的上下文中，平稳性指序列的统计特性（如均值和方差）随时间保持不变的特性。然而，对于间歇性数据，序列通常不满足平稳性假设，因为需求量在不同时间段可能差异很大。

Croston SBA 方法不基于数据序列平稳性的假设。相反，它侧重于使用简单的指数平滑技术分别对间歇性需求的频率和水平进行建模。这使得无需序列平稳性即可捕获需求发生模式并估计条件需求率。

​

加载库和数据

提示

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

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 SBA 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 中实现 CrostonSBA

​

加载库

from statsforecast import StatsForecast
from statsforecast.models import CrostonSBA

​

实例化模型

导入并实例化模型。设置参数有时很棘手。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 = [CrostonSBA()]

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

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

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

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

• fallback_model：如果模型失败时使用的备用模型。

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

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

​

拟合模型

sf.fit(df=train)

StatsForecast(models=[CrostonSBA])

让我们看看 Croston SBA Model 的结果。我们可以通过以下指令来观察它

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

{'mean': array([26.04749601]),
 'fitted': array([      nan,  0.      ,  4.75    , ..., 29.088629, 29.088629,
        29.088629], dtype=float32),
 'sigma': np.float32(49.512943)}

​

预测方法

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

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

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

• h（整数）：表示未来 h 步的预测。在此例中，向前预测 500 小时。

此处的预测对象是一个新的数据框，其中包含模型名称列和 y hat 值，以及不确定性区间列。根据您的计算机，此步骤应花费大约 1 分钟。（如果想将速度提高到几秒，请移除 ARIMA 和 Theta 等自动模型。）

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（整数）：表示未来 h 步的预测。在此例中，向前预测 500 小时。

此处的预测对象是一个新的数据框，其中包含模型名称列和 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（整数）：表示预测未来 h 步。在此例中，向前预测 500 小时。

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

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

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.