标准 Theta 模型

​

目录

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

​

引言

Theta 方法 (Assimakopoulos & Nikolopoulos, 2000, 以下简称 A&N) 应用于非季节性或去季节化的时间序列，其中去季节化通常通过乘法经典分解来执行。该方法通过所谓的 theta 系数（用 ${\theta}_1$ 和 ${\theta}_2$ 表示，其中 ${\theta}_1, {\theta}_2 \in \mathbb{R}$），将原始时间序列分解为两条新线，这些系数应用于数据的二阶差分。当 ${\theta}<1$ 时，二阶差分会减小，从而更好地逼近序列的长期行为（Assimakopoulos, 1995）。如果 ${\theta}$ 等于零，则新线是一条直线。当 ${\theta}>1$ 时，局部曲率增加，放大了时间序列的短期变动（A&N）。生成的新线称为 theta 线，此处用 $\text{Z}(\theta_1)$ 和 $\text{Z}(\theta_2)$ 表示。这些线与原始数据具有相同的平均值和斜率，但局部曲率根据 $\theta$ 系数的值被过滤或增强。

换句话说，分解过程的优势在于可以利用数据中通常无法通过原始时间序列外推完全捕获和建模的信息。theta 线可以被视为新的时间序列，并使用适当的预测方法单独进行外推。一旦每条 theta 线的外推完成，通过组合方案进行重组，以计算原始时间序列的点预测。组合长期以来被认为是预测文献中的一种有用的实践（例如，Clemen, 1989, Makridakis and Winkler, 1983, Petropoulos et al., 2014），因此将其应用于 Theta 方法有望产生更准确和稳健的预测。

Theta 方法在选择 theta 线的数量、theta 系数、外推方法以及组合这些方法以获得稳健预测方面具有很强的灵活性。然而，A&N 提出了一种简化版本，仅使用两条具有预设 $\theta$ 系数的 theta 线，其中 ${\theta}_1 =0$ 对应的 theta 线使用线性回归 (LR) 模型进行外推，而 ${\theta}_2 =2$ 对应的 theta 线使用简单指数平滑 (SES) 进行外推。最终预测通过将两条 theta 线的外推结果以相等权重组合而产生。

Theta 方法的性能也得到了其他实证研究的证实（例如 Nikolopoulos et al., 2012, Petropoulos and Nikolopoulos, 2013）。此外，Hyndman 和 Billah (2003)，以下简称 H&B，表明带漂移的简单指数平滑模型 (SES-d) 是 Theta 方法简化版的统计模型。最近，Thomakos 和 Nikolopoulos (2014) 提供了额外的理论见解，而 Thomakos 和 Nikolopoulos (2015) 推导出了将该方法应用于多元时间序列的新理论公式，并研究了双变量 Theta 方法预计优于单变量方法的条件。尽管取得了这些进展，鉴于其简洁性和优越的预测性能，我们认为 Theta 方法值得预测界给予更多关注。

Theta 方法的一个关键方面是，根据定义，它是动态的。可以选择不同的 theta 线，并使用相等或不等的权重组合生成的预测结果。然而，AN 通过将 theta 系数固定为预定义值来限制这一重要属性。

​

标准 Theta 模型

Assimakopoulos 和 Nikolopoulo 在标准 theta 模型中提出 Theta 线是以下方程的解：

其中 $Y_1, \cdots , Y_T$ 表示原始时间序列数据，$DX_t = (X_t − X_{t−1})$。初始值 $\zeta_1$ 和 $\zeta_2$ 通过最小化 $\sum_{i=1}^{T} [Y_t - \zeta_t (\theta) ]^2$ 获得。然而，方程 (1) 的解析解由下式给出：

其中 $A_T$ 和 $B_T$ 是对 $Y_1, \cdots,Y_T$ 进行简单线性回归的最小二乘系数，它们仅取决于原始数据，并由下式给出：

从这个角度来看，可以将 Theta 线理解为直接应用于数据的线性回归模型的函数。实际上，Theta 方法对未来 h 步的预测是 $\zeta(0)$ 和 $\zeta(2)$ 的线性外推的特别组合（50% - 50%）。

• 当 $\theta < 1$ 应用于数据的二阶差分时，分解过程由一个 theta 系数定义，该系数会减小二阶差分并改善序列行为的逼近。

• 如果 $\theta = 0$，分解后的线将变成一条恒定直线。（见图）

• 如果 $\theta > 1$，则分析序列的短期波动显示出更多的局部曲率（见图）

我们将上述设置称为标准 Theta 方法。构建 theta 方法的步骤如下：

1. 去季节化：首先，对时间序列数据进行统计显著季节性行为测试。如果时间序列满足以下条件，则具有季节性：

$|\rho_m| > q_{1- \frac{\alpha}{2} } \sqrt{\frac{1+2 \sum_{i=1}^{m-1} \rho_{i}^{2} }{T} }$

其中 ρk 表示滞后 $k$ 自相关函数，$m$ 是季节周期内的周期数（例如，月度数据为 12），$T$ 是样本大小，$q$ 是标准正态分布的分位数函数，$(1 − a)\%$ 是置信水平。Assimakopoulos 和 Nikolopoulo [标准 Theta 模型] 选择了 90% 的置信水平。如果时间序列被识别为季节性，则通过经典分解方法进行去季节化，假设季节成分具有乘法关系。

1. 分解： 第二步是将季节调整后的时间序列分解为两条 Theta 线：线性回归线 $\zeta(0)$ 和 theta 线 $\zeta(2)$。

2. 外推： $\zeta(2)$ 使用简单指数平滑 (SES) 进行外推，而 $\zeta(0)$ 则作为正常的线性回归线进行外推。

3. 组合：最终预测是使用相等权重组合两条 $\theta$ 线预测结果。

4. 季节性恢复： 如果第一步存在季节性，则最终预测结果乘以相应的季节指数。

​

加载库和数据

提示

需要安装 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/milk_production.csv", usecols=[1,2])
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)

​

自相关图

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();

​

时间序列分解

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

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

• 水平：这是随时间变化的平均主要值。
• 趋势：趋势是指导致时间序列出现增长或下降模式的值。
• 季节性：这是时间序列中短期发生的周期性事件，导致时间序列出现短期增长或下降模式。
• 残差/噪声：这是时间序列中的随机变动。

将这些组成部分随时间组合起来就形成了时间序列。大多数时间序列由水平和噪声/残差组成，而趋势或季节性是可选的。

如果季节性和趋势是时间序列的一部分，则会影响预测值。因为预测时间序列的模式可能与之前的时间序列不同。

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

​

加法时间序列

如果时间序列的组成部分相加形成时间序列，则该时间序列称为加法时间序列。通过可视化，如果时间序列的增长或下降模式在整个序列中相似，则可以判断该时间序列是加法的。任何加法时间序列的数学函数可以表示为： $y(t) = level + Trend + seasonality + noise$

​

乘法时间序列

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

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

​

加性

from statsmodels.tsa.seasonal import seasonal_decompose 
a = seasonal_decompose(df["y"], model = "additive", period=12)
a.plot();

​

乘性

from statsmodels.tsa.seasonal import seasonal_decompose 
a = seasonal_decompose(df["y"], model = "Multiplicative", period=12)
a.plot();

​

将数据分成训练集和测试集

我们将数据分成两部分：1. 用于训练 Theta 模型的数据 2. 用于测试模型的数据

对于测试数据，我们将使用最近 12 个月的数据来测试和评估模型的性能。

train = df[df.ds<='1974-12-01'] 
test = df[df.ds>'1974-12-01']

train.shape, test.shape

((156, 3), (12, 3))

​

使用 StatsForecast 实现 StandardTheta

​

加载库

from statsforecast import StatsForecast
from statsforecast.models import Theta

​

实例化模型

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

season_length = 12 # Monthly data 
horizon = len(test) # number of predictions

models = [Theta(season_length=season_length, 
                decomposition_type="additive")] # multiplicative   additive

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

models: 模型列表。从可用模型中选择所需的模型并导入它们。

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

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

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

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

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

​

拟合模型

sf.fit(df=train)

StatsForecast(models=[Theta])

我们来看看 Theta 模型的结果。可以使用以下指令查看

result=sf.fitted_[0,0].model_
print(result.keys())
print(result['fit'])

dict_keys(['mse', 'amse', 'fit', 'residuals', 'm', 'states', 'par', 'n', 'modeltype', 'mean_y', 'decompose', 'decomposition_type', 'seas_forecast', 'fitted'])
results(x=array([225.82002697,   0.76015625]), fn=10.638733596938778, nit=19, simplex=array([[241.83142594,   0.76274414],
       [225.82002697,   0.76015625],
       [212.41789302,   0.76391602]]))

现在我们来可视化模型的残差。

如我们所见，上面获得的结果是一个字典输出，要从字典中提取每个元素，我们将使用 .get() 函数来提取元素，然后将其保存到 pd.DataFrame() 中。

residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual

import scipy.stats as stats

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

residual.plot(ax=axs[0,0])
axs[0,0].set_title("Residuals");

sns.distplot(residual, ax=axs[0,1]);
axs[0,1].set_title("Density plot - Residual");

stats.probplot(residual["residual Model"], dist="norm", plot=axs[1,0])
axs[1,0].set_title('Plot Q-Q')

plot_acf(residual,  lags=35, ax=axs[1,1],color="fuchsia")
axs[1,1].set_title("Autocorrelation");

plt.show();

​

预测方法

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

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

forecast 方法接受两个参数：预测未来 h (时间范围) 和 level。

• h (int): 表示未来 h 个时间步的预测。在本例中，即未来 12 个月。

• level (浮点数列表): 这个可选参数用于概率预测。设置预测区间的水平（或置信百分比）。例如，level=[90] 意味着模型预期真实值有 90% 的几率落在此区间内。

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

# Prediction
Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat

values=sf.forecast_fitted_values()
values.head()

使用 forecast 方法添加 95% 置信区间

sf.forecast(df=train, h=horizon, level=[95])

sf.plot(train, Y_hat)

​

带有置信区间的 Predict 方法

使用 predict 方法生成预测。

predict 方法接受两个参数：预测未来 h (时间范围) 和 level。

• h (int): 表示未来 h 个时间步的预测。在本例中，即未来 12 个月。

• level (浮点数列表): 这个可选参数用于概率预测。设置预测区间的水平（或置信百分比）。例如，level=[95] 意味着模型预期真实值有 95% 的几率落在此区间内。

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

此步骤应少于 1 秒。

sf.predict(h=horizon)

forecast_df = sf.predict(h=horizon, level=[80,95]) 
forecast_df

sf.plot(train, test.merge(forecast_df), level=[80, 95])

​

交叉验证

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

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

下图展示了这种交叉验证策略

​

执行时间序列交叉验证

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

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

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

• df: 训练数据框

• h (int): 表示正在预测的未来 h 个时间步。在本例中，即未来 12 个月。

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

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

crossvalidation_df = sf.cross_validation(df=train,
                                         h=horizon,
                                         step_size=12,
                                         n_windows=3)

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

• unique_id: 序列标识符
• ds: 日期时间戳或时间索引
• cutoff: n_windows 的最后一个日期时间戳或时间索引。
• y: 真实值
• "model": 包含模型名称和拟合值的列。

​

模型评估

现在我们将根据预测结果评估模型，我们将使用不同类型的评估指标 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,
)

​

致谢

我们要感谢 Naren Castellon 撰写本教程。

​

参考文献

1. Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016)。“优化 Theta 方法及其与状态空间模型关系的模型”。International Journal of Forecasting.
2. V. Assimakopoulos, K. Nikolopoulos，“The theta model: a decomposition approach to forecasting”（Theta 模型：一种用于预测的分解方法）
3. Nixtla 参数.
4. Pandas 可用频率.
5. Rob J. Hyndman 和 George Athanasopoulos (2018)。“预测原则与实践 (第三版)”.
6. 季节周期 - Rob J Hyndman.