TSB 模型

​

目录

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

​

引言

Teunter-Syntetos-Babai (TSB) 模型是一种用于时间序列中的库存管理和需求预测领域的模型。它由 Teunter、Syntetos 和 Babai 于 2001 年提出，是 Croston 需求预测模型的扩展。

TSB 模型专门用于预测具有间歇性需求特征的产品需求，即需求期之后接着无需求期的产品。它旨在处理具有大量零值和非零观测值之间间隔变化的时间序列数据。

TSB 模型基于两个主要组成部分：水平模型和间隔模型。水平模型估计发生需求时的需求水平，而间隔模型估计需求发生之间的间隔。这两个组成部分结合起来生成对未来需求的准确预测。

TSB 模型已被证明在间歇性需求预测方面是有效的，并已广泛应用于各种工业部门。然而，需要注意的是，需求预测还有其他可用模型和方法，选择合适的模型将取决于数据的具体特征及其应用的环境。

​

TSB 模型

TSB (Teunter, Syntetos 和 Babai) 是 2011 年提出的一种新方法，该方法用每个周期更新的需求概率取代了需求间隔。原因是 Croston 方法只在需求发生时更新需求，但在现实生活中，有很多情况存在大量的零需求，因此，预测结果不适合估计因信息过时而导致的陈旧风险。

在 TSB 方法中，$D_t$ 表示周期 $t$ 的需求发生指示符，因此

如果 $D_t=0$，则

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

$D_t=D'_{t-1}+\beta (0- D'_{t-1})$

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

$D'_t=D'_{t-1}+\beta(1-D'_{t-1})$

因此，预测由下式给出

$Y'_t=D'_t \cdot Z'_t$

其中

• $Y'_t:$ 每期平均需求
• $Z_t:$ 周期 $t$ 的实际需求
• $Z'_t:$ 两次正需求之间的时间
• $D'_t:$ 估计周期 $t$ 结束时发生需求的概率
• $\alpha, \beta:$ 平滑常数，$0 \leq \alpha, \beta \leq 1$

​

TSB 一般特性

用于时间序列的 Teunter-Syntetos-Babai (TSB) 模型具有以下特性

1. 间歇性需求建模：TSB 模型专门设计用于预测间歇性需求，其特点是无需求期后接着需求期。该模型有效地解决了需求的这一特性。

2. 水平和间隔组成部分：TSB 模型基于两个主要组成部分：水平模型和间隔模型。水平模型估计发生需求时的需求水平，而间隔模型估计需求发生之间的间隔。

3. 处理包含大量零值的数据：TSB 模型可以有效地处理包含大量零值的时间序列数据，这在间歇性需求中很常见。该模型在预测过程中会适当考虑这些零值。

4. 指数平滑：TSB 模型使用指数平滑方法来估计需求水平和发生间隔。指数平滑是时间序列预测中广泛使用的技术。

5. 置信区间估计：TSB 模型为生成的预测提供置信区间估计。这使得对预测相关的不确定性有所了解，并有助于决策制定。

6. 简单易于实现：与其他更复杂的方法相比，TSB 模型相对简单且易于实现。它不需要对需求分布做出复杂的假设，并且可以以实用方式应用。

这些是 Teunter-Syntetos-Babai 模型在时间序列和间歇性需求预测背景下的一些基本特性。

​

加载库和数据

提示

需要 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 plotly.subplots import make_subplots

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. 用于训练 TSB 模型的数据。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 实现 TSB 模型

​

加载库

from statsforecast import StatsForecast
from statsforecast.models import TSB

​

构建模型

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

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

models = [TSB(alpha_d=0.8, alpha_p=0.9)]

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

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

• freq: 表示数据频率的字符串。（参见 pandas 可用的频率别名。）

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

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

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

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

​

拟合模型

sf.fit(df=train)

StatsForecast(models=[TSB])

让我们看看 TSB 模型的结果。我们可以通过以下指令观察：

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

{'mean': array([65.58645721]),
 'fitted': array([        nan,  0.        ,  9.        , ..., 14.937817  ,
         1.4937816 ,  0.14937817], dtype=float32),
 'sigma': np.float32(63.87893)}

​

预测方法

如果您想在具有多个时间序列或模型的生产环境中提高速度，我们建议使用 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=) 中的性能，每隔两个月 (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 series (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 and George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”.
7. 季节周期 - Rob J Hyndman.