# 此单元格将不会被渲染，但用于隐藏警告并限制显示的行数。
import warnings
warnings.filterwarnings("ignore")
import logging
logging.getLogger('statsforecast').setLevel(logging.ERROR)
import pandas as pd
pd.set_option('display.max_rows', 6)优化的 Theta 模型
使用
Statsforecast的优化 Theta 模型的逐步指南。
在本指南中，我们将熟悉主要的 StatsForecast 类以及一些相关的方法，例如 StatsForecast.plot、StatsForecast.forecast 和 StatsForecast.cross_validation 等。
目录
介绍
优化的Theta模型是一种基于将时间序列分解为三个组成部分：趋势、季节性和噪声的时间序列预测方法。该模型首先预测长期趋势和季节性，然后利用噪声来调整短期预测。研究表明，优化的Theta模型比其他时间序列预测方法更准确，尤其是在处理具有复杂趋势和季节性的时间序列时。
优化的Theta模型由Athanasios N. Antoniadis和Nikolaos D. Tsonis于2013年开发。该模型基于Theta预测方法，而Theta方法是由George E. P. Box和Gwilym M. Jenkins于1976年开发的。Theta方法是基于将时间序列分解为三种成分：趋势、季节性和噪声的时间序列预测方法。Theta模型随后预测长期趋势和季节性，并利用噪声来调整短期预测。
Theta优化模型通过使用优化算法来寻找模型的最佳参数，从而改进了Theta方法。优化算法基于赤池信息量（AIC）损失函数，这是一个衡量模型与数据拟合优度的标准。优化算法寻找能够最小化AIC函数的参数。
研究表明，优化的Theta模型比其他时间序列预测方法更准确，特别是在处理具有复杂趋势和季节性的时间序列时。该模型已被用于预测各种时间序列，包括销售、生产、价格和天气。
以下是优化的Theta模型的一些优点：
- 比其他时间序列预测方法更准确。
- 使用简单。
- 可用于预测各种时间序列。
- 具有灵活性，可以适应不同的场景。
如果您在寻找一种易于使用且准确的时间序列预测方法，优化的Theta模型是一个不错的选择。
优化后的Theta模型可以应用于多个领域，包括：
- 销售： 优化后的Theta模型可以用于预测产品或服务的销售。这可以帮助公司做出有关生产、库存和营销的决策。
- 生产： 优化后的Theta模型可以用于预测商品或服务的生产。这可以帮助公司确保他们有能力满足需求，并避免过度生产。
- 价格： 优化后的Theta模型可以用于预测商品或服务的价格。这可以帮助公司做出有关定价和营销策略的决策。
- 天气： 优化后的Theta模型可以用于预测天气。这可以帮助公司做出有关农业生产、旅行计划和风险管理的决策。
- 其他： 优化后的Theta模型还可以用于预测其他类型的时间序列，包括交通、能源需求和人口。
优化后的Theta模型是一种强大的工具，可以用于提高时间序列预测的准确性。它易于使用，可应用于多个领域。如果您正在寻找一种改善时间序列预测的工具，优化后的Theta模型是一个不错的选择。
优化的 Theta 模型 (OTM)
假设时间序列 \(Y_1, \cdots Y_n\) 是非季节性的，或者已通过乘法经典分解方法进行了季节性调整。
设 \(X_t\) 为两个 theta 线的线性组合，
\[ \begin{equation} X_t=\omega \text{Z}_t (\theta_1) +(1-\omega) \text{Z}_t (\theta_2) \tag 1 \end{equation} \]
其中 \(\omega \in [0,1]\) 是权重参数。假设 \(\theta_1 <1\) 且 \(\theta_2 \geq 1\)，权重 \(\omega\) 可以通过以下公式推导得到
\[ \begin{equation} \omega:=\omega(\theta_1, \theta_2)=\frac{\theta_2 -1}{\theta_2 -\theta_1} \tag 2 \end{equation} \]
从公式 (1) 和 (2) 可以清楚地看出 \(X_t=Y_t, \ t=1, \cdots n\)，即权重的计算方式能够正确重现原始序列。
定理 1: 设 \(\theta_1 <1\) 且 \(\theta_2 \geq 1\)。我们将证明
- 线性系统 \(X_t=Y_t\) 对所有 \(t=1, \cdots, n\)，其中 \(X_t\) 由公式 (4) 给出，具有唯一解
\[\omega= (\theta_2 -1)/(\theta_2 - \theta_1)\]
- 选择非最优权重 \(\omega_{\delta} =\omega + \delta\) 的误差与简单线性回归模型的误差成正比。
在定理 1 中，我们证明了解是唯一的，并且不选择最优权重 (\(\omega\) 和 \(1-\omega\)) 的误差与线性回归模型的误差成正比。因此，STheta 方法仅需设置 \(\theta_1=0\) 和 \(\theta_2=2\)，从公式 (2) 我们得到 \(\omega=0.5\)。因此，公式 (1) 和 (2) 使我们能够构建一个 Theta 模型的推广，保持原始时间序列在任何 theta 线 \(\text{Z}_t (\theta_1)\) 和 \(\text{Z}_t (\theta_2)\) 的重组合性质。
为了维持长期成分的建模，并与 STheta 方法保持公正的比较，在本研究中我们固定 \(\theta_1=0\)，并专注于短期成分的优化，\(\theta_2=0\) 且 \(\theta \geq 1\)。因此，\(\theta\) 是目前唯一需要估计的参数。现在，theta 分解表示为
\[Y_t=(1-\frac{1}{\theta}) (\text{A}_n+\text{B}_n t)+ \frac{1}{\theta} \text{Z}_t (\theta), \ t=1, \cdots , n\]
在原点计算的 \(h\) 步预测值为
\[ \begin{equation} \hat Y_{n+h|n} = (1-\frac{1}{\theta}) [\text{A}_n+\text{B}_n (n+h)]+ \frac{1}{\theta} \tilde {\text{Z}}_{n+h|n} (\theta) \tag 3 \end{equation} \]
其中 \(\tilde {\text{Z}}_{n+h|n} (\theta)=\tilde {\text{Z}}_{n+1|n} (\theta)=\alpha \sum_{i=0}^{n-1}(1-\alpha)^i \text{Z}_{n-i}(\theta)+(1-\alpha)^n \ell_{0}^{*}\) 是通过 SES 模型对 \(\text{Z}_t(\theta)\) 的外推，其中 \(\ell_{0}^{*} \in \mathbb{R}\) 是初始水平参数，\(\alpha \in (0,1)\) 是平滑参数。注意，对于 \(\theta=2\)，公式 (3) 对应于 STheta 算法的步骤 4。经过一些代数运算，我们可以写成
\[ \begin{equation} \tilde {\text{Z}}_{n+1|n} (\theta)=\theta \ell{n}+(1-\theta) \{ \text{A}_n [1-(1-\alpha)^n] + \text{B}_n [n+(1-\frac{1}{\alpha}) [1-(1-\alpha)^n] ] \} \tag 4 \end{equation} \]
其中 \(\ell_{t}=\alpha Y_t +(1-\alpha) \ell_{t-1}\)，对于 \(t=1, \cdots, n\)，且 \(\ell_{0}=\ell_{0}^{*}/\theta\)。
根据公式 (3)、(4)，我们建议四种随机方法。这些方法因参数 \(\theta\) 而异，\(\theta\) 可以固定为二或优化，且系数 \(\text{A}_n\) 和 \(\text{B}_n\)，可以是固定的或动态函数。为了建立状态空间模型，有助于采用 \(\mu_{t}\) 作为原点 \(t-1\) 的一步预测值，\(\varepsilon_{t}\) 作为相应的加性误差，即如果 \(\mu_{t}= \hat Y_{t|t-1}\)，则 \(\varepsilon_{t}=Y_t - \mu_{t}\)。我们假设 \(\{ \varepsilon_{t} \}\) 是一个均值为零、方差为 \(\sigma^2\) 的高斯白噪声过程。
关于优化的Theta模型
设 \(\text{A}_n\) 和 \(\text{B}_n\) 为所有 \(t=1, \cdots, n\) 的固定系数，使得方程 (3)、(4) 形成由以下公式给出的状态空间模型：
\[ \begin{equation} Y_t=\mu_{t}+\varepsilon_{t} \tag 5 \end{equation} \]
\[ \begin{equation} \mu_{t}=\ell_{t-1}+(1-\frac{1}{\theta}) \{(1-\alpha)^{t-1} \text{A}_n +[\frac{1-(1-\alpha)^t}{\alpha} \text{B}_n] \tag 6 \end{equation} \]
\[ \begin{equation} \ell_{t}=\alpha Y_t +(1-\alpha)\ell_{t-1} \tag 7 \end{equation} \]
参数 \(\ell_{0} \in \mathbb{R}\), \(\alpha \in (0,1)\) 和 \(\theta \in [1,\infty)\) 。参数 \(\theta\) 将与 \(\alpha\) 和 \(\ell_{0}\) 一同进行估计。我们称之为优化的Theta模型（OTM）。
在原点 \(n\) 处的 \(h\) 步预测为
\[\hat Y_{n+h|n}=E[Y_{n+h}|Y_1,\cdots, Y_n]=\ell_{n}+(1-\frac{1}{\theta}) \{(1-\alpha)^n \text{A}_n +[(h-1) + \frac{1-(1-\alpha)^{n+1}}{\alpha}] \text{B}_n \}\]
这等同于方程 (3)。条件方差 \(\text{Var}[Y_{n+h}|Y_1, \cdots, Y_n]=[1+(h-1)\alpha^2]\sigma^2\) 可以从状态空间模型中轻松计算得出。因此，\(Y_{n+h}\) 的 \((1-\alpha)\%\) 预测区间为 \[\hat Y_{n+h|n} \ \pm \ q_{1-\alpha/2} \sqrt{[1+(h-1)\alpha^2 ]\sigma^2 }\]
对于 \(\theta=2\)，OTM 复制了 STheta 方法的预测；此后我们将这个特定案例称为标准Theta模型（STM）。
定理 2： SES-d \((\ell_{0}^{**}, \alpha, b)\) 模型，其中 \(\ell_{0}^{**} \in \mathbb{R}, \alpha \in (0,1)\) 和 \(b \in \mathbb{R}\) 等价于 \(\text{OTM} (\ell_{0}, \alpha, \theta )\)，其中 \(\ell_{0} \in \mathbb{R}\) 和 \(\theta \geq 1\)，如果
\[\ell_{0}^{**} = \ell_{0} + (1- \frac{1}{\theta} )A_n \ \ and \ \ b=(1-\frac{1}{\theta} )B_n\]
在定理 2 中，我们表明 OTM 在数学上等价于 SES-d 模型。作为定理 2 的一个推论，STM 在数学上等价于 SES-d，且 \(b=\frac{1}{2} \text{B}_n\)。因此，对于 \(\theta=2\)，该推论也再次确认了 H&B 关于 STheta 与 SES-d 模型之间关系的结果。
加载库和数据
Tip
需要使用 Statsforecast。有关安装的方法，请参见 说明。
接下来，我们导入绘图库并配置绘图样式。
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.style.use('grayscale') # 五三八 灰度 经典
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()| month | production | |
|---|---|---|
| 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, engine="matplotlib")
自相关图
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. 用于训练我们的优化的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))现在让我们绘制训练数据和测试数据。
sns.lineplot(train,x="ds", y="y", label="Train", linestyle="--")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.title("")
plt.ylabel("Monthly Milk Production")
plt.xlabel("Monthly")
plt.show()
使用StatsForecast实现OptimizedTheta
为了更好地了解OptimizedTheta Model函数的参数，以下是它们的列表。有关更多信息，请访问 文档。
season_length : int
每单位时间的观察次数。例如：24小时数据。
decomposition_type : str
季节分解类型，'multiplicative'（默认）或'additive'。
alias : str
模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
用于计算符合预测区间的信息。
默认情况下，模型将计算本地预测
区间。加载库
from statsforecast import StatsForecast
from statsforecast.models import OptimizedTheta 实例化模型
导入并实例化模型。设置参数有时会很棘手。大师 Rob Hyndmann 在季节性周期上的文章对 season_length 可能会有所帮助。
season_length = 12 # 月度数据
horizon = len(test) # 预测数量
models = [OptimizedTheta(season_length=season_length,
decomposition_type="additive")] # 乘法 加法我们通过实例化一个新的 StatsForecast 对象来拟合模型，使用以下参数：
models: 模型列表。从模型中选择您想要的模型并导入它们。
freq:字符串，表示数据的频率。（请参阅 pandas 可用的频率。）n_jobs:n_jobs: int， 并行处理时使用的作业数，使用 -1 表示所有核心。fallback_model:如果某个模型失败，则使用的模型。
任何设置都将被传递到构造函数中。然后您调用它的 fit 方法，并传入历史数据框。
sf = StatsForecast(df=train,
models=models,
freq='MS',
n_jobs=-1)拟合模型
sf.fit()StatsForecast(models=[OptimizedTheta])让我们看看我们的 优化的 Theta 模型 (OTM) 的结果。我们可以通过以下指令进行观察：
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([-83.14191626, 0.73681394, 12.45013763]), fn=10.448217519858634, nit=47, simplex=array([[-58.73988124, 0.7441127 , 11.69842922],
[-49.97233449, 0.73580297, 11.41787513],
[-83.14191626, 0.73681394, 12.45013763],
[-77.04867427, 0.73498431, 11.99254037]]))现在让我们可视化我们模型的残差。
如我们所见，上述结果输出为字典，要从字典中提取每个元素，我们将使用 .get() 函数提取元素，然后将其保存到 pd.DataFrame() 中。
residual=pd.DataFrame(result.get("residuals"), columns=["residual Model"])
residual| residual Model | |
|---|---|
| 0 | -271.899414 |
| 1 | -114.671692 |
| 2 | 4.768066 |
| ... | ... |
| 153 | -60.233887 |
| 154 | -92.472839 |
| 155 | -44.143982 |
156 rows × 1 columns
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 不会存储拟合值，并且在分布式环境中具有高度可扩展性。
预测方法接受两个参数：预测下一个 h（时间范围）和 level。
h (整数):代表预测未来 h 步。此案例中，预测 12 个月后。level (浮点数列表):这是一个可选参数，用于概率预测。设置您的预测区间的水平（或置信百分位）。例如，level=[90]意味着该模型预计真实值在该区间内的概率为 90%。
这里的预测对象是一个新的数据框，包含模型名称和 y hat 值的列，以及不确定性区间的列。根据您的计算机，这一步预计耗时约 1 分钟。（如果您想将时间缩短到几秒钟，可以移除像 ARIMA 和 Theta 的自动模型）
Y_hat = sf.forecast(horizon, fitted=True)
Y_hat| ds | OptimizedTheta | |
|---|---|---|
| unique_id | ||
| 1 | 1975-01-01 | 839.682800 |
| 1 | 1975-02-01 | 802.071838 |
| 1 | 1975-03-01 | 896.117126 |
| ... | ... | ... |
| 1 | 1975-10-01 | 824.135498 |
| 1 | 1975-11-01 | 795.691162 |
| 1 | 1975-12-01 | 833.316345 |
12 rows × 2 columns
让我们可视化拟合值
values=sf.forecast_fitted_values()
values.head()| ds | y | OptimizedTheta | |
|---|---|---|---|
| unique_id | |||
| 1 | 1962-01-01 | 589.0 | 860.899414 |
| 1 | 1962-02-01 | 561.0 | 675.671692 |
| 1 | 1962-03-01 | 640.0 | 635.231934 |
| 1 | 1962-04-01 | 656.0 | 614.731323 |
| 1 | 1962-05-01 | 727.0 | 609.770752 |
StatsForecast.plot(values)
添加95%的置信区间与预测方法
sf.forecast(h=horizon, level=[95])| ds | OptimizedTheta | OptimizedTheta-lo-95 | OptimizedTheta-hi-95 | |
|---|---|---|---|---|
| unique_id | ||||
| 1 | 1975-01-01 | 839.682800 | 742.509583 | 955.414307 |
| 1 | 1975-02-01 | 802.071838 | 643.581360 | 945.119263 |
| 1 | 1975-03-01 | 896.117126 | 710.785217 | 1065.057495 |
| ... | ... | ... | ... | ... |
| 1 | 1975-10-01 | 824.135498 | 555.948975 | 1084.320312 |
| 1 | 1975-11-01 | 795.691162 | 503.148010 | 1036.519531 |
| 1 | 1975-12-01 | 833.316345 | 530.259949 | 1106.636963 |
12 rows × 4 columns
Y_hat=Y_hat.reset_index()
Y_hat| unique_id | ds | OptimizedTheta | |
|---|---|---|---|
| 0 | 1 | 1975-01-01 | 839.682800 |
| 1 | 1 | 1975-02-01 | 802.071838 |
| 2 | 1 | 1975-03-01 | 896.117126 |
| ... | ... | ... | ... |
| 9 | 1 | 1975-10-01 | 824.135498 |
| 10 | 1 | 1975-11-01 | 795.691162 |
| 11 | 1 | 1975-12-01 | 833.316345 |
12 rows × 3 columns
# 将预测结果与真实值合并
test['unique_id'] = test['unique_id'].astype(int)
Y_hat1 = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat1| ds | y | unique_id | OptimizedTheta | |
|---|---|---|---|---|
| 0 | 1975-01-01 | 834 | 1 | 839.682800 |
| 1 | 1975-02-01 | 782 | 1 | 802.071838 |
| 2 | 1975-03-01 | 892 | 1 | 896.117126 |
| ... | ... | ... | ... | ... |
| 9 | 1975-10-01 | 827 | 1 | 824.135498 |
| 10 | 1975-11-01 | 797 | 1 | 795.691162 |
| 11 | 1975-12-01 | 843 | 1 | 833.316345 |
12 rows × 4 columns
fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat1]).set_index('ds')
plot_df[['y', "OptimizedTheta"]].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Monthly Milk Production', fontsize=20)
ax.set_xlabel('Month [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)
带置信区间的预测方法
要生成预测，请使用预测方法。
预测方法接受两个参数：预测下一个 h（代表时间范围）和 level。
h (int):表示未来 h 步的预测。在这种情况下，表示未来 12 个月。level (list of floats):此可选参数用于概率预测。设置预测区间的置信水平（或置信百分位）。例如，level=[95]意味着模型预计真实值在该区间内的概率为 95%。
此处的预测对象是一个新的数据框，其中包含模型名称和 y hat 值的列，以及不确定性区间的列。
此步骤应少于 1 秒。
sf.predict(h=horizon) | ds | OptimizedTheta | |
|---|---|---|
| unique_id | ||
| 1 | 1975-01-01 | 839.682800 |
| 1 | 1975-02-01 | 802.071838 |
| 1 | 1975-03-01 | 896.117126 |
| ... | ... | ... |
| 1 | 1975-10-01 | 824.135498 |
| 1 | 1975-11-01 | 795.691162 |
| 1 | 1975-12-01 | 833.316345 |
12 rows × 2 columns
forecast_df = sf.predict(h=horizon, level=[80,95])
forecast_df| ds | OptimizedTheta | OptimizedTheta-lo-80 | OptimizedTheta-hi-80 | OptimizedTheta-lo-95 | OptimizedTheta-hi-95 | |
|---|---|---|---|---|---|---|
| unique_id | ||||||
| 1 | 1975-01-01 | 839.682800 | 766.665955 | 928.326233 | 742.509583 | 955.414307 |
| 1 | 1975-02-01 | 802.071838 | 704.290100 | 899.335876 | 643.581360 | 945.119263 |
| 1 | 1975-03-01 | 896.117126 | 761.334717 | 1007.408630 | 710.785217 | 1065.057495 |
| ... | ... | ... | ... | ... | ... | ... |
| 1 | 1975-10-01 | 824.135498 | 623.904114 | 996.567322 | 555.948975 | 1084.320312 |
| 1 | 1975-11-01 | 795.691162 | 576.546753 | 975.490967 | 503.148010 | 1036.519531 |
| 1 | 1975-12-01 | 833.316345 | 606.713989 | 1033.886230 | 530.259949 | 1106.636963 |
12 rows × 6 columns
我们可以使用pandas函数pd.concat()将预测结果与历史数据合并，然后能够利用这个结果进行图表绘制。
pd.concat([df, forecast_df]).set_index('ds')| y | unique_id | OptimizedTheta | OptimizedTheta-lo-80 | OptimizedTheta-hi-80 | OptimizedTheta-lo-95 | OptimizedTheta-hi-95 | |
|---|---|---|---|---|---|---|---|
| ds | |||||||
| 1962-01-01 | 589.0 | 1 | NaN | NaN | NaN | NaN | NaN |
| 1962-02-01 | 561.0 | 1 | NaN | NaN | NaN | NaN | NaN |
| 1962-03-01 | 640.0 | 1 | NaN | NaN | NaN | NaN | NaN |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 1975-10-01 | NaN | NaN | 824.135498 | 623.904114 | 996.567322 | 555.948975 | 1084.320312 |
| 1975-11-01 | NaN | NaN | 795.691162 | 576.546753 | 975.490967 | 503.148010 | 1036.519531 |
| 1975-12-01 | NaN | NaN | 833.316345 | 606.713989 | 1033.886230 | 530.259949 | 1106.636963 |
180 rows × 7 columns
现在，让我们可视化我们的预测结果和时间序列的历史数据，同时绘制我们在以95%置信度进行预测时获得的置信区间。
def plot_forecasts(y_hist, y_true, y_pred, models):
_, ax = plt.subplots(1, 1, figsize = (20, 7))
y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
df_plot = pd.concat([y_hist, y_true]).set_index('ds').tail(12*10)
df_plot[['y'] + models].plot(ax=ax, linewidth=3 , )
colors = ['green', "lime"]
ax.fill_between(df_plot.index,
df_plot['OptimizedTheta-lo-80'],
df_plot['OptimizedTheta-lo-80'],
alpha=.20,
color='orange',
label='OptimizedTheta_level_80')
ax.fill_between(df_plot.index,
df_plot['OptimizedTheta-lo-95'],
df_plot['OptimizedTheta-hi-95'],
alpha=.3,
color='lime',
label='OptimizedTheta_level_95')
ax.set_title('', fontsize=22)
ax.set_ylabel("Montly Mil Production", fontsize=20)
ax.set_xlabel('Month', fontsize=20)
ax.legend(prop={'size': 20})
ax.grid(True)
plt.show()plot_forecasts(train, test, forecast_df, models=['OptimizedTheta'])
让我们使用Statsforecast中提供的plot函数绘制相同的图表，如下所示。
sf.plot(df, forecast_df, level=[95])
交叉验证
在之前的步骤中，我们利用历史数据来预测未来。然而，为了评估其准确性，我们还希望了解模型在过去的表现。为了评估模型在数据上的准确性和稳健性，执行交叉验证。
对于时间序列数据，交叉验证是通过在历史数据上定义滑动窗口并预测随后的时间段来完成的。这种交叉验证形式使我们能够更好地估计模型在更广泛的时间实例上的预测能力，同时保持训练集中的数据连续，这是我们的模型所要求的。
下图描绘了这种交叉验证策略：
执行时间序列交叉验证
时间序列模型的交叉验证被视为最佳实践，但大多数实现都非常缓慢。statsforecast库将交叉验证实现为分布式操作，使得执行过程不那么耗时。如果您有大型数据集，也可以使用Ray、Dask或Spark在分布式集群中执行交叉验证。
在这种情况下，我们希望评估每个模型在过去5个月的表现(n_windows=5)，每隔两个月进行一次预测(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| ds | cutoff | y | OptimizedTheta | |
|---|---|---|---|---|
| unique_id | ||||
| 1 | 1972-01-01 | 1971-12-01 | 826.0 | 828.836365 |
| 1 | 1972-02-01 | 1971-12-01 | 799.0 | 792.592346 |
| 1 | 1972-03-01 | 1971-12-01 | 890.0 | 883.269592 |
| ... | ... | ... | ... | ... |
| 1 | 1974-10-01 | 1973-12-01 | 812.0 | 812.183838 |
| 1 | 1974-11-01 | 1973-12-01 | 773.0 | 783.898376 |
| 1 | 1974-12-01 | 1973-12-01 | 813.0 | 821.124329 |
36 rows × 4 columns
模型评估
现在我们可以使用适当的准确性指标来计算预测的准确性。在这里，我们将使用均方根误差 (RMSE)。为此，我们首先需要 install datasetsforecast，这是一个由 Nixtla 开发的 Python 库，其中包含计算 RMSE 的函数。
%%capture
!pip install datasetsforecastfrom datasetsforecast.losses import rmse计算RMSE的函数需要两个参数：
- 实际值。
- 预测值，这里指的是“优化的Theta模型（OTM）”。
rmse = rmse(crossvalidation_df['y'], crossvalidation_df["OptimizedTheta"])
print("RMSE using cross-validation: ", rmse)RMSE using cross-validation: 14.504839正如您所注意到的，我们使用交叉验证结果来评估我们的模型。
现在我们将使用预测结果来评估我们的模型，我们将使用不同类型的指标MAE, MAPE, MASE, RMSE, SMAPE来评估精度。
from datasetsforecast.losses import (mae, mape, mase, rmse, smape)def evaluate_performace(y_hist, y_true, y_pred, model):
y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
evaluation = {}
evaluation[model] = {}
for metric in [mase, mae, mape, rmse, smape]:
metric_name = metric.__name__
if metric_name == 'mase':
evaluation[model][metric_name] = metric(y_true['y'].values,
y_true[model].values,
y_hist['y'].values, seasonality=12)
else:
evaluation[model][metric_name] = metric(y_true['y'].values, y_true[model].values)
return pd.DataFrame(evaluation).Tevaluate_performace(train, test, Y_hat, model="OptimizedTheta")| mae | mape | mase | rmse | smape | |
|---|---|---|---|---|---|
| OptimizedTheta | 6.740209 | 0.782753 | 0.30312 | 8.701501 | 0.778689 |
鸣谢
我们要感谢 Naren Castellon 撰写本教程。
参考文献
- Kostas I. Nikolopoulos, Dimitrios D. Thomakos. 使用Theta方法进行预测-理论与应用。2019年，约翰·威利父子公司。
- Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). “优化Theta方法的模型及其与状态空间模型的关系”。《国际预测杂志》。
- Nixtla参数。
- Pandas可用频率。
- Rob J. Hyndman和George Athanasopoulos (2018). “预测原则与实践，时间序列交叉验证”。。
- 季节性周期- Rob J. Hyndman。
Give us a ⭐ on Github