动态优化 Theta 模型

引言

**动态优化 Theta 模型 (DOTM)** 是 StatsForecast 中经典 Theta 模型的一个变体。它结合了另外两个扩展的关键特性：**优化 Theta 模型 (OTM)** 和 **动态标准 Theta 模型 (DSTM)**。

DOTM 在标准 Theta 模型的基础上引入了两个主要改进：theta 参数的**优化**以及模型组件随时间进行的**动态更新**。

**优化**：与 OTM 一样，此版本会根据数据自动搜索最佳的 theta 值，而不是依赖于固定参数。这种灵活性使模型能够更好地适应具有复杂季节性或趋势模式的序列。
**动态更新**：与 DSTM 一样，DOTM 会随着新数据的到来不断更新其内部组件。这使得它非常适合非平稳序列，其中潜在的数据结构随时间演变。

DOTM 还支持**季节性分解**，由 decomposition_type 参数控制。您可以选择：- 'multiplicative'（默认），它假设季节性影响与序列的水平成比例缩放，或 - 'additive'，它假设季节性影响在绝对量级上保持不变。

动态优化 Theta 模型是 Theta 系列中最灵活的模型，在预测具有不断变化的趋势和季节性的序列时特别有效。

加载库和数据

提示

将需要 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()

	月份	产量
0	1962-01-01	589
1	1962-02-01	561
2	1962-03-01	640
3	1962-04-01	656
4	1962-05-01	727

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

	ds	y	unique_id
0	1962-01-01	589	1
1	1962-02-01	561	1
2	1962-03-01	640	1
3	1962-04-01	656	1
4	1962-05-01	727	1

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

时间序列分解

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

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

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

随着时间的推移组合这些组件会形成时间序列。大多数时间序列包含水平和噪声/残差，趋势或季节性是可选值。

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

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

加法时间序列

如果时间序列的组件相加构成时间序列。则该时间序列称为加法时间序列。通过可视化，我们可以说，如果时间序列的增加或减少模式在整个序列中相似，则该时间序列是加法的。任何加法时间序列的数学函数可以表示为： $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();

将数据拆分为训练集和测试集

让我们将数据划分为数据集

用于训练我们的 动态优化 Theta 模型(DOTM) 的数据。
用于测试我们的模型的数据

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

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

train.shape, test.shape

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

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

sns.lineplot(train,x="ds", y="y", label="Train", linestyle="--")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.title("Monthly Milk Production");
plt.show()

在 StatsForecast 中实现 DynamicOptimizedTheta

加载库

from statsforecast import StatsForecast
from statsforecast.models import DynamicOptimizedTheta

实例化模型

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

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

# We call the model that we are going to use
models = [DynamicOptimizedTheta(season_length=season_length, 
                decomposition_type="additive")] # multiplicative   additive

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

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

freq: 一个字符串，指示数据的频率。（参见 pandas 可用频率。）
n_jobs: n_jobs: int，并行处理中使用的作业数，使用 -1 表示所有核心。
fallback_model: 如果模型失败时使用的模型。

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

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

拟合模型

sf.fit(df=train)

StatsForecast(models=[DynamicOptimizedTheta])

让我们看看我们的 动态优化 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([250.83206219,   0.75624902,   4.67964777]), fn=10.697554045462667, nit=55, simplex=array([[237.42074763,   0.75306547,   4.46023813],
       [250.83206219,   0.75624902,   4.67964777],
       [257.16444246,   0.75229688,   4.42377059],
       [256.90853867,   0.75757957,   4.43171897]]))

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

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

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

	残差模型
0	-18.247106
1	-75.757706
2	6.001494
…	…
153	-59.747044
154	-91.901521
155	-43.503294

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 (float 列表): 此可选参数用于概率预测。设置预测区间的置信水平（或百分位数）。例如，level=[90] 表示模型期望实际值有 90% 的时间落在此区间内。

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

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

	unique_id	ds	DynamicOptimizedTheta
0	1	1975-01-01	839.259705
1	1	1975-02-01	801.399170
2	1	1975-03-01	895.189148
…	…	…	…
9	1	1975-10-01	821.271240
10	1	1975-11-01	792.530518
11	1	1975-12-01	829.854553

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

	unique_id	ds	y	DynamicOptimizedTheta
0	1	1962-01-01	589.0	607.247131
1	1	1962-02-01	561.0	636.757690
2	1	1962-03-01	640.0	633.998535
3	1	1962-04-01	656.0	608.461243
4	1	1962-05-01	727.0	604.808899

StatsForecast.plot(values)

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

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

	unique_id	ds	DynamicOptimizedTheta	DynamicOptimizedTheta-lo-95	DynamicOptimizedTheta-hi-95
0	1	1975-01-01	839.259705	741.952332	955.151001
1	1	1975-02-01	801.399170	641.867920	946.045776
2	1	1975-03-01	895.189148	707.189087	1066.356812
…	…	…	…	…	…
9	1	1975-10-01	821.271240	546.081726	1088.193481
10	1	1975-11-01	792.530518	494.623718	1037.459839
11	1	1975-12-01	829.854553	519.661133	1108.213867

带有置信区间的 predict 方法

使用 predict 方法生成预测。

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

h (int): 表示预测未来 h 步。在本例中，提前 12 个月。
level (float 列表): 此可选参数用于概率预测。设置预测区间的置信水平（或百分位数）。例如，level=[95] 表示模型期望实际值有 95% 的时间落在此区间内。

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

此步骤应少于 1 秒。

sf.predict(h=horizon)

	unique_id	ds	DynamicOptimizedTheta
0	1	1975-01-01	839.259705
1	1	1975-02-01	801.399170
2	1	1975-03-01	895.189148
…	…	…	…
9	1	1975-10-01	821.271240
10	1	1975-11-01	792.530518
11	1	1975-12-01	829.854553

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

	unique_id	ds	DynamicOptimizedTheta	DynamicOptimizedTheta-lo-80	DynamicOptimizedTheta-hi-80	DynamicOptimizedTheta-lo-95	DynamicOptimizedTheta-hi-95
0	1	1975-01-01	839.259705	766.142090	928.025513	741.952332	955.151001
1	1	1975-02-01	801.399170	702.981262	899.884216	641.867920	946.045776
2	1	1975-03-01	895.189148	760.125916	1008.335022	707.189087	1066.356812
…	…	…	…	…	…	…	…
9	1	1975-10-01	821.271240	617.391724	996.698364	546.081726	1088.193481
10	1	1975-11-01	792.530518	568.303162	975.070312	494.623718	1037.459839
11	1	1975-12-01	829.854553	598.098267	1035.476196	519.661133	1108.213867

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": 包含模型名称和拟合值的列。

crossvalidation_df

	unique_id	ds	cutoff	y	DynamicOptimizedTheta
0	1	1972-01-01	1971-12-01	826.0	828.692017
1	1	1972-02-01	1971-12-01	799.0	792.444092
2	1	1972-03-01	1971-12-01	890.0	883.122620
…	…	…	…	…	…
33	1	1974-10-01	1973-12-01	812.0	810.304688
34	1	1974-11-01	1973-12-01	773.0	781.804688
35	1	1974-12-01	1973-12-01	813.0	818.811096

模型评估

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

	unique_id	指标	DynamicOptimizedTheta
0	1	mae	6.861949
1	1	mape	0.008045
2	1	mase	0.308595
3	1	rmse	8.647459
4	1	smape	0.004010

参考

CrostonSBA 模型动态标准 Theta 模型

在本页

目录
引言
加载库和数据
读取数据
使用 plot 方法探索数据
自相关图
时间序列分解
加法时间序列
乘法时间序列
加法
乘法
将数据拆分为训练集和测试集
在 StatsForecast 中实现 DynamicOptimizedTheta
加载库
实例化模型
拟合模型
预测方法
带有置信区间的 predict 方法
交叉验证
执行时间序列交叉验证
模型评估
参考

入门

教程

操作指南

分布式

实验

模型参考

API 参考

贡献

目录

引言

加载库和数据

读取数据

使用 plot 方法探索数据

自相关图

时间序列分解

加法时间序列

乘法时间序列

加法

乘法

将数据拆分为训练集和测试集

在 StatsForecast 中实现 DynamicOptimizedTheta

加载库

实例化模型

拟合模型

预测方法

带有置信区间的 predict 方法

交叉验证

执行时间序列交叉验证

模型评估

参考

入门

教程

操作指南

分布式

实验

模型参考

API 参考

贡献

​目录

​引言

​加载库和数据

​读取数据

​使用 plot 方法探索数据

​自相关图

​时间序列分解

​加法时间序列

​乘法时间序列

​加法

​乘法

​将数据拆分为训练集和测试集

​在 StatsForecast 中实现 DynamicOptimizedTheta

​加载库

​实例化模型

​拟合模型

​预测方法

​带有置信区间的 predict 方法

​交叉验证

​执行时间序列交叉验证

​模型评估

​参考

目录

引言

加载库和数据

读取数据

使用 plot 方法探索数据

自相关图

时间序列分解

加法时间序列

乘法时间序列

加法

乘法

将数据拆分为训练集和测试集

在 StatsForecast 中实现 DynamicOptimizedTheta

加载库

实例化模型

拟合模型

预测方法

带有置信区间的 predict 方法

交叉验证

执行时间序列交叉验证

模型评估

参考