动态标准Theta模型

动态标准Theta模型 (DSTM)

动态标准Theta模型是优化动态Theta模型的一种特殊情况变体。

此外，对于 $\theta=2$ ，我们采用了一种Theta的随机方法，下文将其称为动态标准Theta模型 (DSTM)。

加载库和数据

提示

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

使用绘图方法探索数据

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

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

让我们将数据分为两组：1. 用于训练`Dynamic Standard Theta Model`的数据 2. 用于测试模型的数据

对于测试数据，我们将使用最近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中实现DynamicStandardTheta

加载库

from statsforecast import StatsForecast
from statsforecast.models import DynamicTheta

实例化模型

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

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

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

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

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

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

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

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

拟合模型

sf.fit(df=train)

StatsForecast(models=[DynamicTheta])

让我们看看我们的`Dynamic Standard Theta model`的结果。我们可以通过以下指令观察它

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([393.28739991,   0.76875   ]), fn=10.787112115489622, nit=20, simplex=array([[399.92916541,   0.771875  ],
       [393.28739991,   0.76875   ],
       [384.74798713,   0.771875  ]]))

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

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

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

	残差模型
0	-18.247131
1	-46.247131
2	17.140198
…	…
153	-58.941711
154	-91.055420
155	-42.624939

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 (list of floats): 这个可选参数用于概率预测。设置预测区间的置信水平（或置信百分位）。例如，level=[90]意味着模型预期真实值有90%的时间会落在此区间内。

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

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

	unique_id	ds	DynamicTheta
0	1	1975-01-01	838.531555
1	1	1975-02-01	800.154968
2	1	1975-03-01	893.430786
…	…	…	…
9	1	1975-10-01	815.959351
10	1	1975-11-01	786.716431
11	1	1975-12-01	823.539368

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

	unique_id	ds	y	DynamicTheta
0	1	1962-01-01	589.0	607.247131
1	1	1962-02-01	561.0	607.247131
2	1	1962-03-01	640.0	622.859802
3	1	1962-04-01	656.0	606.987793
4	1	1962-05-01	727.0	605.021179

StatsForecast.plot(values)

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

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

	unique_id	ds	DynamicTheta	DynamicTheta-lo-95	DynamicTheta-hi-95
0	1	1975-01-01	838.531555	741.237366	954.407166
1	1	1975-02-01	800.154968	640.697205	945.673096
2	1	1975-03-01	893.430786	703.900635	1065.418701
…	…	…	…	…	…
9	1	1975-10-01	815.959351	536.422791	1086.643433
10	1	1975-11-01	786.716431	484.476593	1033.687134
11	1	1975-12-01	823.539368	509.187256	1104.107788

使用predict方法生成置信区间

使用predict方法生成预测。

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

h (int): 表示预测未来h步。在这种情况下，提前12个月。
level (list of floats): 这个可选参数用于概率预测。设置预测区间的置信水平（或置信百分位）。例如，level=[95]意味着模型预期真实值有95%的时间会落在此区间内。

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

此步骤应少于1秒。

sf.predict(h=horizon)

	unique_id	ds	DynamicTheta
0	1	1975-01-01	838.531555
1	1	1975-02-01	800.154968
2	1	1975-03-01	893.430786
…	…	…	…
9	1	1975-10-01	815.959351
10	1	1975-11-01	786.716431
11	1	1975-12-01	823.539368

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

	unique_id	ds	DynamicTheta	DynamicTheta-lo-80	DynamicTheta-hi-80	DynamicTheta-lo-95	DynamicTheta-hi-95
0	1	1975-01-01	838.531555	765.423828	927.285339	741.237366	954.407166
1	1	1975-02-01	800.154968	701.099854	899.316162	640.697205	945.673096
2	1	1975-03-01	893.430786	758.326416	1007.631165	703.900635	1065.418701
…	…	…	…	…	…	…	…
9	1	1975-10-01	815.959351	608.699463	992.552673	536.422791	1086.643433
10	1	1975-11-01	786.716431	558.429810	970.648376	484.476593	1033.687134
11	1	1975-12-01	823.539368	588.706787	1031.564941	509.187256	1104.107788

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: 索引。如果您不喜欢使用索引，只需运行crossvalidation_df.resetindex()。
ds: 日期戳或时间索引
cutoff: n_windows的最后一个日期戳或时间索引。
y: 真实值
"model": 包含模型名称和拟合值的列。

crossvalidation_df

	unique_id	ds	cutoff	y	DynamicTheta
0	1	1972-01-01	1971-12-01	826.0	827.107239
1	1	1972-02-01	1971-12-01	799.0	789.924194
2	1	1972-03-01	1971-12-01	890.0	879.664429
…	…	…	…	…	…
33	1	1974-10-01	1973-12-01	812.0	804.398560
34	1	1974-11-01	1973-12-01	773.0	775.329285
35	1	1974-12-01	1973-12-01	813.0	811.767639

模型评估

现在我们将根据预测结果评估我们的模型，我们将使用不同类型的指标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	指标	DynamicTheta
0	1	MAE	8.182119
1	1	MAPE	0.009736
2	1	MASE	0.367965
3	1	RMSE	9.817624
4	1	SMAPE	0.004874

参考资料

动态优化Theta模型 GARCH模型

在本页上

目录
动态标准Theta模型 (DSTM)
加载库和数据
读取数据
使用绘图方法探索数据
自相关图
时间序列分解
加法时间序列
乘法时间序列
加法
乘法
将数据分割为训练集和测试集
在StatsForecast中实现DynamicStandardTheta
加载库
实例化模型
拟合模型
预测方法
使用predict方法生成置信区间
交叉验证
执行时间序列交叉验证
模型评估
参考资料

入门

教程

操作指南

分布式

实验

模型参考

API参考

贡献

目录