用 Python 轻松实现时间序列预测:Darts 快速入门 Quickstart

文中内容仅限技术学习与代码实践参考,市场存在不确定性,技术分析需谨慎验证,不构成任何投资建议。

Darts

Darts 是一个 Python 库,用于对时间序列进行用户友好型预测和异常检测。它包含多种模型,从 ARIMA 等经典模型到深度神经网络。所有预测模型都能以类似 scikit-learn 的方式使用 fit()predict() 函数。该库还可以轻松地对模型进行回溯测试,将多个模型的预测结果结合起来,并将外部数据考虑在内。Darts 支持单变量和多变量时间序列和模型。基于 ML 的模型可以在包含多个时间序列的潜在大型数据集上进行训练,其中一些模型还为概率预测提供了丰富的支持。

快速入门

在本笔记本中,我们将介绍该库的主要功能:

  • 安装 Darts
  • 构建和操作时间序列
  • 训练预测模型并生成预测
  • 回测
  • 机器学习和全局模型
  • 协变量:使用外部数据
  • 回归预测模型
  • 训练样本权重
  • 预测起点偏移
  • 概率预测
  • 模型集成
  • 滤波模型
  • 异常检测

这里我们仅展示最简化的"入门"示例。更深入的信息请参考我们的用户指南示例笔记本

安装 Darts

建议使用虚拟环境。主要有两种安装方式:

使用 pip:

pip install darts

使用 conda:

conda install -c conda-forge -c pytorch u8darts-all

如果遇到问题或需要安装不同版本(避免某些依赖),请查阅详细安装指南

首先导入必要的模块:

## 本地工作时修复 Python 路径
from utils import fix_pythonpath_if_working_locally

fix_pythonpath_if_working_locally()

import warnings

warnings.filterwarnings("ignore", category=FutureWarning)
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

构建和操作 TimeSeries

TimeSeries 是 Darts 中的核心数据类。一个 TimeSeries 表示单变量或多变量时间序列,带有时间索引。时间索引可以是 pandas.DatetimeIndex(包含日期时间)或 pandas.RangeIndex(包含整数,适用于表示无特定时间戳的顺序数据)。在某些情况下,TimeSeries 甚至可以表示概率序列,例如用于获取置信区间。Darts 中的所有模型都使用 TimeSeries 并生成 TimeSeries

读取数据并构建 TimeSeries

可通过工厂方法轻松构建 TimeSeries

  • DataFrame 构建,使用 TimeSeries.from_dataframe()文档)。支持多种框架后端如 pandas.DataFramepolars.DataFramepyarrow.table 等(参见 (a))。
  • 从时间索引和对应值数组构建,使用 TimeSeries.from_times_and_values()文档)。
  • 从 NumPy 值数组构建,使用 TimeSeries.from_values()文档)。
  • Series 构建,使用 TimeSeries.from_series()文档)。同样支持多种序列后端(参见 (a))。
  • 从 Pandas DataFrame 按组提取多个 TimeSeries,使用 TimeSeries.from_group_dataframe()文档)。
  • xarray.DataArray 构建,使用 TimeSeries.from_xarray()文档)。
  • 从 CSV 文件构建,使用 TimeSeries.from_csv()文档)。
  • 从 JSON 文件构建,使用 TimeSeries.from_json()文档)。
  • 或使用内置的 TimeSeries 生成函数(文档)。

(a) 我们使用 narwhals 作为 DataFrame 和 Series 库之间的兼容层。所有支持的后端参见 narwhals 文档

让我们从 DataFrame 创建一个多元(两列)TimeSeries。该数据框包含正弦和余弦序列,时间索引为日频率。

from darts import TimeSeries
from darts.utils.utils import generate_index

## 生成正弦波和余弦波
x = np.linspace(0, 2 * np.pi, 100)
sine_vals = np.sin(x)
cosine_vals = np.cos(x)

## 生成日频率的日期时间索引
dates = generate_index("2020-01-01", length=len(x), freq="D")

## 创建 DataFrame(这里使用 pandas,但也可用其他后端如 polars、pyarrow 等)
df = pd.DataFrame({"sine": sine_vals, "cosine": cosine_vals, "time": dates})

series = TimeSeries.from_dataframe(df, time_col="time")
series.plot();

img

加载内置数据集

也可以加载我们现有的数据集(文档)。

下面直接加载航空乘客序列数据集获取 TimeSeries

from darts.datasets import AirPassengersDataset

series = AirPassengersDataset().load()
series.plot();

导出 TimeSeries

TimeSeries 可通过多种方式导出:

  • 导出到 DataFrame,使用 TimeSeries.to_dataframe()文档)。支持多种框架后端(参见上面的 (a))。
  • 导出到 NumPy 值数组,使用 TimeSeries.values()文档)或 TimeSeries.all_values()文档)。
  • 导出到 Series,使用 TimeSeries.to_series()文档)。同样支持多种序列后端(参见上面的 (a))。
  • 导出到 xarray.DataArray,使用 TimeSeries.data_array()文档)。
  • 导出到 CSV 文件,使用 TimeSeries.to_csv()文档)。
  • 导出到 JSON 文件,使用 TimeSeries.to_json()文档)。
  • 以及更多!

以下是将数据导出到 polars.DataFrame 的示例。Polars 是可选依赖项,请先安装(安装指南在此):

series[-5:].to_dataframe(backend="polars", time_as_index=False)
Month#Passengers
datetime[ns]f64
1960-08-01 00:00:00606.0
1960-09-01 00:00:00508.0
1960-10-01 00:00:00461.0
1960-11-01 00:00:00390.0
1960-12-01 00:00:00432.0

部分 TimeSeries 操作

TimeSeries 支持多种操作 - 以下是几个示例。

分割
可在序列的分数位置、pandas Timestamp 或整数索引值处分割:

series1, series2 = series.split_after(0.75)
series1.plot()
series2.plot();

img

切片:

series1, series2 = series[:-36], series[-36:]
series1.plot()
series2.plot();

img

算术运算:

series_noise = TimeSeries.from_times_and_values(
    series.time_index, np.random.randn(len(series))
)
(series / 2 + 20 * series_noise - 10).plot();

img

堆叠
沿新维度拼接以生成新的多元序列:

(series / 50).stack(series_noise).plot();

img

映射:

series.map(np.log).plot();

img

同时映射时间戳和值:

series.map(lambda ts, x: x / ts.days_in_month).plot();

img

添加日期时间属性作为额外维度(生成多元序列):

(series / 20).add_datetime_attribute("month").plot();

img

添加二元节假日成分:

(series / 200).add_holidays("US").plot();

img

差分:

series.diff().plot();

img

填充缺失值(使用 utils 函数)
缺失值用 np.nan 表示。

from darts.utils.missing_values import fill_missing_values

values = np.arange(50, step=0.5)
values[10:30] = np.nan
values[60:95] = np.nan
series_ = TimeSeries.from_values(values)

(series_ - 10).plot(label="with missing values (shifted below)")
fill_missing_values(series_).plot(label="without missing values")

img

创建训练集和验证集

接下来我们将 TimeSeries 拆分为训练集和验证集。注意:通常还应保留测试集直到流程结束才使用。为简化起见,这里仅构建训练集和验证集。

训练集包含 1958 年 1 月(不含)之前的值,验证集包含剩余部分:

train, val = series.split_before(pd.Timestamp("19580101"))
train.plot(label="training")
val.plot(label="validation")

img

训练预测模型并生成预测

使用基础模型

Darts 包含一组"朴素"基线模型,可用于了解预期的最低准确度。例如,NaiveSeasonal(K) 模型总是"重复" K 个时间步之前的值。

K=1 时,该模型仅重复训练序列的最后一个值:

from darts.models import NaiveSeasonal

naive_model = NaiveSeasonal(K=1)
naive_model.fit(train)
naive_forecast = naive_model.predict(36)

series.plot(label="actual")
naive_forecast.plot(label="naive forecast (K=1)")

img

在时间序列上拟合模型并进行预测非常容易。所有模型都包含 fit()predict() 函数。这两个函数与 Scikit-learn 类似,只是它专门用于时间序列。fit() 函数的参数是要拟合模型的训练时间序列,而 predict() 函数的参数是要进行预测的时间步数(训练序列结束后)。

检查季节性

上面的模型可能过于简单。我们可以利用数据的季节性改进模型。数据似乎存在年度季节性,这可以通过自相关函数(ACF)确认,并突出显示滞后 m=12

from darts.utils.statistics import check_seasonality, plot_acf

plot_acf(train, m=12, alpha=0.05, max_lag=24)

img

ACF 在 x=12 处出现峰值,表明存在年度季节性趋势(红色突出)。蓝色区域表示 α=5% 置信水平下的统计显著性。我们还可以对每个候选周期 m 进行季节性统计检查:

for m in range(2, 25):
    is_seasonal, period = check_seasonality(train, m=m, alpha=0.05)
    if is_seasonal:
        print(f"There is seasonality of order {period}.")
There is seasonality of order 12.
改进的模型

让我们使用季节性周期 12 再次尝试 NaiveSeasonal 模型:

seasonal_model = NaiveSeasonal(K=12)
seasonal_model.fit(train)
seasonal_forecast = seasonal_model.predict(36)

series.plot(label="actual")
seasonal_forecast.plot(label="naive forecast (K=12)")

img

这有所改进,但我们仍缺失趋势。幸运的是,还有另一个捕捉趋势的朴素基线模型 NaiveDrift。该模型仅生成线性预测,斜率由训练集的首尾值确定:

from darts.models import NaiveDrift

drift_model = NaiveDrift()
drift_model.fit(train)
drift_forecast = drift_model.predict(36)

combined_forecast = drift_forecast + seasonal_forecast - train.last_value()

series.plot()
combined_forecast.plot(label="combined")
drift_forecast.plot(label="drift")

img

发生了什么?我们只是拟合了一个朴素的漂移模型,并将其预测添加到我们之前的季节性预测中。我们还从结果中减去训练集的最后一个值,以便最终的组合预测从正确的偏移量开始。

计算误差指标

这看起来已经相当不错,而我们尚未使用任何非朴素模型。实际上,任何模型都应能超越此结果。

那么我们需要超越的误差是多少?我们将使用平均绝对百分比误差 (MAPE)(注意实践中通常有充分理由使用 MAPE - 这里因其便捷性和尺度无关性而使用)。在 Darts 中只需简单调用函数:

from darts.metrics import mape

print(
    f"Mean absolute percentage error for the combined naive drift + seasonal: {mape(series, combined_forecast):.2f}%."
)
Mean absolute percentage error for the combined naive drift + seasonal: 5.66%.

darts.metrics 包含更多用于比较时间序列的指标。当两个序列未对齐时,指标仅比较共同的切片,并可并行计算大量序列对 - 但让我们先不深入讨论。

快速尝试多个模型

Darts 的设计使得可以统一方式轻松训练和验证多个模型。让我们再训练几个模型,并计算它们在验证集上的 MAPE:

from darts.models import AutoARIMA, ExponentialSmoothing, Theta


def eval_model(model):
    model.fit(train)
    forecast = model.predict(len(val))
    print(f"model {model} obtains MAPE: {mape(val, forecast):.2f}%")


eval_model(ExponentialSmoothing())
eval_model(AutoARIMA())
eval_model(Theta())
model ExponentialSmoothing() obtains MAPE: 5.11%
model AutoARIMA() obtains MAPE: 12.70%
model Theta() obtains MAPE: 8.15%
使用 Theta 方法搜索超参数

Theta 模型实现了 Assimakopoulos 和 Nikolopoulos 的 Theta 方法。该方法取得了一定成功,尤其在 M3 竞赛中。

虽然 Theta 参数值通常设为 0,但我们的实现支持可变值以进行参数调优。让我们尝试寻找最佳 Theta 值:

## 尝试 50 个不同值搜索最佳 theta 参数
thetas = 2 - np.linspace(-10, 10, 50)

best_mape = float("inf")
best_theta = 0

for theta in thetas:
    model = Theta(theta)
    model.fit(train)
    pred_theta = model.predict(len(val))
    res = mape(val, pred_theta)

    if res < best_mape:
        best_mape = res
        best_theta = theta
best_theta_model = Theta(best_theta)
best_theta_model.fit(train)
pred_best_theta = best_theta_model.predict(len(val))

print(f"Lowest MAPE is: {mape(val, pred_best_theta):.2f}, with theta = {best_theta}.")
Lowest MAPE is: 4.40, with theta = -3.5102040816326543.
train.plot(label="train")
val.plot(label="true")
pred_best_theta.plot(label="prediction")

img

我们可以观察到,就 MAPE 而言,具有 best_theta 的模型是迄今为止最好的。

回测:模拟历史预测

此时我们有一个在验证集表现良好的模型,这很好。但我们如何知道如果历史上使用该模型会获得什么性能?

回测模拟使用给定模型在历史上会获得的预测。生成过程可能需要时间,因为(默认情况下)每次模拟预测时间推进时都会重新训练模型。

此类模拟预测始终相对于预测范围定义,即预测时间与预测时间点之间的时间步数。在下面的示例中,我们模拟对未来 3 个月(与预测时间相比)的预测。调用 historical_forecasts() 的结果(默认)是一个 TimeSeries,仅包含每个 3 个月提前预测的最后一个预测值:

hfc_params = {
    "series": series,
    "start": pd.Timestamp("1956-01-01"),  ## 也可以是序列起始比例的浮点数
    "forecast_horizon": 3,
    "verbose": True,
}
historical_fcast_theta = best_theta_model.historical_forecasts(
    last_points_only=True, **hfc_params
)

series.plot(label="data")
historical_fcast_theta.plot(label="backtest 3-months ahead forecast (Theta)")
print(f"MAPE = {mape(series, historical_fcast_theta):.2f}%")
historical forecasts: 100% 58/58 [00:00<00:00, 194.38it/s] MAPE = 7.99%

img

通过设置 last_points_only=False,我们还可以检索每个历史预测的所有预测值。使用 stride 参数定义连续预测之间的步长。我们将其设为 3 个月,以便将预测连接成单个 TimeSeries

historical_fcast_theta_all = best_theta_model.historical_forecasts(
    last_points_only=False, stride=3, **hfc_params
)

series.plot(label="data")
for idx, hfc in enumerate(historical_fcast_theta_all):
    hfc.plot(label=f"forecast {idx}")

from darts import concatenate

historical_fcast_theta_all = concatenate(historical_fcast_theta_all, axis=0)
print(f"MAPE = {mape(series, historical_fcast_theta_all):.2f}%")
historical forecasts: 100% 20/20 [00:00<00:00, 97.92it/s] MAPE = 5.61%

img

所以,当我们对验证集上的最佳模型进行回测时,它似乎表现得不那么好了(我是不是听到过拟合了 😄)

为了更仔细地观察误差,我们还可以使用 backtest() 方法获取模型本来应该获得的所有原始误差(比如 MAPE 误差):

best_theta_model = Theta(best_theta)

raw_errors = best_theta_model.backtest(
    metric=mape, reduction=None, last_points_only=False, stride=1, **hfc_params
)

from darts.utils.statistics import plot_hist

plot_hist(
    raw_errors,
    bins=np.arange(0, max(raw_errors), 1),
    title="Individual backtest error scores (histogram)",
)
historical forecasts: 100% 58/58 [00:00<00:00, 202.54it/s]

img

最后,使用 backtest() 我们还可以更简单地查看历史预测的平均误差:

average_error = best_theta_model.backtest(
    metric=mape,
    reduction=np.mean,  ## 这实际上是默认的
    **hfc_params,
)

print(f"Average error (MAPE) over all historical forecasts: {average_error:.2f}")
historical forecasts: 100% 58/58 [00:00<00:00, 197.04it/s]
Average error (MAPE) over all historical forecasts: 6.30

例如,我们也可以指定参数 reduction=np.median 来获取中值 MAPE。

上面,backtest() 每次调用时都会重新计算历史预测值。我们也可以使用一些预先计算的预测值来更快地获得结果!

hfc_precomputed = best_theta_model.historical_forecasts(
    last_points_only=False, stride=1, **hfc_params
)
new_error = best_theta_model.backtest(
    historical_forecasts=hfc_precomputed, last_points_only=False, stride=1, **hfc_params
)

print(f"Average error (MAPE) over all historical forecasts: {new_error:.2f}")
historical forecasts: 100% 58/58 [00:00<00:00, 213.58it/s]
Average error (MAPE) over all historical forecasts: 6.30

让我们看一下当前 Theta 模型的拟合值残差,即通过对所有先前点进行模型拟合获得的每个时间点的 1 步预测值与实际观测值之间的差异:

hfc_precomputed = best_theta_model.historical_forecasts(
    last_points_only=False, stride=1, **hfc_params
)
new_error = best_theta_model.backtest(
    historical_forecasts=hfc_precomputed, last_points_only=False, stride=1, **hfc_params
)

print(f"Average error (MAPE) over all historical forecasts: {new_error:.2f}")
historical forecasts: 100% 58/58 [00:00<00:00, 213.58it/s]
Average error (MAPE) over all historical forecasts: 6.30

让我们看一下当前 Theta 模型的拟合值残差,即通过对所有先前点进行模型拟合获得的每个时间点的 1 步预测值与实际观测值之间的差异:

from darts.utils.statistics import plot_residuals_analysis

plot_residuals_analysis(best_theta_model.residuals(series))

img

我们可以看到,分布并非以 0 为中心,这意味着我们的 Theta 模型存在偏差。我们还可以在滞后 12 时观察到较大的 ACF 值,这表明残差包含模型未使用的信息。

我们的残差方法实际上更加强大!它可以用于计算 Darts 中的任何每一步指标(请参阅此处的列表),即使是多步预测也是如此。它还支持类似于回测的预先计算的历史预测。

现在,让我们检查一下 3 个月预测中每一步的绝对误差分布。

from darts.metrics import ae

residuals = best_theta_model.residuals(
    historical_forecasts=hfc_precomputed,
    metric=ae,  ## 每个时间步的绝对误差
    last_points_only=False,
    values_only=True,  ## 返回一个 NumPy 数组列表
    **hfc_params,
)
residuals = np.concatenate(residuals, axis=1)[:, :, 0]

fig, ax = plt.subplots()
for forecast_step in range(len(residuals)):
    ax.hist(residuals[forecast_step], label=f"step {forecast_step}", alpha=0.5)
ax.legend()
ax.set_title("Absolute errors per forecast step")

img

我们可以清楚地看到,预测的误差会随着未来的推移而增大。

我们能否使用一个简单的指数平滑模型来做得更好?

model_es = ExponentialSmoothing(seasonal_periods=12)
historical_fcast_es = model_es.historical_forecasts(**hfc_params)

series.plot(label="data")
historical_fcast_es.plot(label="backtest 3-months ahead forecast (Exp. Smoothing)")
print(f"MAPE = {mape(historical_fcast_es, series):.2f}%")
historical forecasts: 100% 58/58 [00:01<00:00, 31.05it/s]
MAPE = 4.42%

img

这好多了!在这种情况下,当我们以 3 个月的预测期进行回测时,我们得到的平均绝对百分比误差约为 4-5%。

plot_residuals_analysis(model_es.residuals(series, verbose=True))
historical forecasts: 100% 120/120 [00:03<00:00, 31.33it/s]

img

残差分析也反映出性能的改善,因为现在我们的残差分布以值 0 为中心,并且 ACF 值虽然不小,但幅度较低。

机器学习和全局模型

Darts 对机器学习和深度学习预测模型有丰富支持;例如:

  • SKLearnModel 可包装任何 sklearn 兼容的回归模型以生成预测(在下文有专门章节)。
  • RNNModel 是灵活的 RNN 实现,可像 DeepAR 一样使用。
  • NBEATSModel 实现 N-BEATS 模型。
  • TFTModel 实现时序融合变换器模型。
  • TCNModel 实现时序卷积网络。

除了支持与其他模型相同的基础 fit()/predict() 接口外,这些模型也是全局模型,因为它们支持在多个时间序列上训练(有时称为元学习)。

这是使用基于 ML 的模型进行预测的关键点:通常,ML 模型(尤其是深度学习模型)需要大量数据训练,这通常意味着大量相关但独立的时间序列。

在 Darts 中,指定多个 TimeSeries 的基本方式是使用 TimeSeriesSequence(例如,简单的 TimeSeries 列表)。

包含两个序列的示例

这些模型可在数千个序列上训练。为简单说明,这里加载两个不同的序列 - 航空乘客数量和每月每头奶牛产奶磅数。为加速训练,我们将序列转换为 np.float32

from darts.datasets import AirPassengersDataset, MonthlyMilkDataset

series_air = AirPassengersDataset().load().astype(np.float32)
series_milk = MonthlyMilkDataset().load().astype(np.float32)

## 将每个序列的最后 36 个月留作验证集:
train_air, val_air = series_air[:-36], series_air[-36:]
train_milk, val_milk = series_milk[:-36], series_milk[-36:]

train_air.plot()
val_air.plot()
train_milk.plot()
val_milk.plot();

img

首先,将这两个序列缩放到 0 到 1 之间,这对大多数 ML 模型有益。我们使用 Scaler

from darts.dataprocessing.transformers import Scaler

scaler = Scaler()
train_air_scaled, train_milk_scaled = scaler.fit_transform([train_air, train_milk])

train_air_scaled.plot()
train_milk_scaled.plot();

img

注意我们如何一次性扩展多个系列。我们还可以通过指定 n_jobs 在多个处理器上并行执行此类操作。

使用深度学习:N-BEATS 示例

注意:神经网络模型 (TorchForecastingModels) 的详细用户指南见此处

接下来构建 N-BEATS 模型。该模型可通过许多超参数调优(如堆栈数、层数等)。为简单起见,这里使用默认超参数。我们必须提供的两个超参数是:

  • input_chunk_length:模型的"回看窗口" - 即神经网络在前向传递中作为输入的历史时间步数。
  • output_chunk_length:模型的"前向窗口" - 即神经网络在前向传递中输出的未来值时间步数。

random_state 参数仅用于确保结果可复现。

Darts 中的大多数神经网络都需要这两个参数。这里我们使用季节性的倍数。现在准备在序列上拟合模型(通过向 fit() 提供包含两个序列的列表):

from darts.models import NBEATSModel

model = NBEATSModel(input_chunk_length=24, output_chunk_length=12, random_state=42)

model.fit([train_air_scaled, train_milk_scaled], epochs=50, verbose=True);
GPU available: True (mps), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs

  | Name            | Type             | Params | Mode
-------------------------------------------------------------
0 | criterion       | MSELoss          | 0      | train
1 | train_criterion | MSELoss          | 0      | train
2 | val_criterion   | MSELoss          | 0      | train
3 | train_metrics   | MetricCollection | 0      | train
4 | val_metrics     | MetricCollection | 0      | train
5 | stacks          | ModuleList       | 6.2 M  | train
-------------------------------------------------------------
6.2 M     Trainable params
1.4 K     Non-trainable params
6.2 M     Total params
24.787    Total estimated model params size (MB)
396       Modules in train mode
0         Modules in eval mode
Epoch 49: 100% 7/7 [00:00<00:00,  7.81it/s, train_loss=0.00179]
`Trainer.fit` stopped: `max_epochs=50` reached.

现在预测未来 36 个月。只需使用 predict()series 参数告知模型要预测的序列。重要的是,output_chunk_length 不直接约束 predict() 的预测范围 n。这里我们使用 output_chunk_length=12 训练模型,并预测 n=36 个月;这只需在后台以自回归方式完成(网络递归消耗其先前的输出)。

pred_air, pred_milk = model.predict(series=[train_air_scaled, train_milk_scaled], n=36)

## 逆缩放:
pred_air, pred_milk = scaler.inverse_transform([pred_air, pred_milk])

plt.figure(figsize=(10, 6))
series_air.plot(label="actual (air)")
series_milk.plot(label="actual (milk)")
pred_air.plot(label="forecast (air)")
pred_milk.plot(label="forecast (milk)")
GPU available: True (mps), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
Predicting DataLoader 0: 100% 1/1 [00:00<00:00, 16.92it/s]

img

我们的预测实际上并不差,考虑到我们使用一个具有默认超参数的模型来捕捉航空乘客和牛奶产量!

该模型似乎能很好地捕捉年度季节性,但错过了航空序列的趋势。下一节中,我们将尝试使用外部数据(协变量)解决此问题。

协变量:使用外部数据

除了目标序列(我们想要预测的序列)之外,Darts 中的许多模型还接受输入中的协变量序列。协变量是我们不想预测但可为模型提供有用额外信息的序列。目标和协变量都可以是多元或单变量。

Darts 中有两种协变量时间序列:

  • past_covariates 是不一定在预测时间之前已知的序列。例如,它们可表示必须测量且无法提前获知的事物。模型在预测时不使用 past_covariates 的未来值(对于全局模型,当 n > output_chunk_length 时因自回归而使用)。有关过去/未来协变量的更多信息,请查看此用户指南
  • future_covariates 是提前已知的序列,直至预测范围。它们可表示日历信息、节假日、天气预报等。接受 future_covariates 的模型在预测时将查看未来值(直至预测范围)。
  • static_covariates 是目标序列不随时间变化的特征。它们直接嵌入目标序列中。可表示产品类别、国家信息等。有关静态协变量的更多信息,请查看此用户指南

img

每个协变量都可能是多元的。如果有多个协变量序列(如月份和年份值),应 stack()concatenate() 它们以获得多元序列。

提供的协变量可以比必要长度更长。Darts 的模型会自动处理相关时间范围的提取! 但如果协变量的时间跨度不足,会收到错误提示。

并非所有模型都支持每种协变量类型。可在此处查看模型列表,了解它们支持的类型。

现在为航空和牛奶序列构建包含月份和年份值的外部协变量。
在下面的单元格中,我们使用 darts.utils.timeseries_generation.datetime_attribute_timeseries() 函数生成包含月份和年份值的序列,并沿 "component" 轴(同 axis=1concatenate() 这些序列以获得每个目标序列一个包含两个分量(月份和年份)的协变量序列。为简单起见,我们直接将月份和年份值缩放到 0 到 1 之间:

from darts import concatenate
from darts.utils.timeseries_generation import datetime_attribute_timeseries as dt_attr

air_covs = concatenate(
    [
        dt_attr(series_air, "month", dtype=np.float32),
        dt_attr(series_air, "year", dtype=np.float32),
    ],
    axis="component",
)

milk_covs = concatenate(
    [
        dt_attr(series_milk, "month", dtype=np.float32),
        dt_attr(series_milk, "year", dtype=np.float32),
    ],
    axis="component",
)

air_covs_scaled, milk_covs_scaled = Scaler().fit_transform([air_covs, milk_covs])
air_covs_scaled.plot()
milk_covs_scaled.plot()
plt.title(
    "one multivariate time series of 2 dimensions, containing covariates for the air series:"
)

img

NBEATSModel 仅支持 past_covariates。因此,即使协变量表示日历信息且提前已知,我们仍将其作为 past_covariates 与 N-BEATS 一起使用。训练时,只需将它们作为 past_covariates 提供给 fit() 函数,顺序与目标相同:

model = NBEATSModel(input_chunk_length=24, output_chunk_length=12, random_state=42)

model.fit(
    [train_air_scaled, train_milk_scaled],
    past_covariates=[air_covs_scaled, milk_covs_scaled],
    epochs=50,
    verbose=True,
);
GPU available: True (mps), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs

  | Name            | Type             | Params | Mode
-------------------------------------------------------------
0 | criterion       | MSELoss          | 0      | train
1 | train_criterion | MSELoss          | 0      | train
2 | val_criterion   | MSELoss          | 0      | train
3 | train_metrics   | MetricCollection | 0      | train
4 | val_metrics     | MetricCollection | 0      | train
5 | stacks          | ModuleList       | 6.6 M  | train
-------------------------------------------------------------
6.6 M     Trainable params
1.7 K     Non-trainable params
6.6 M     Total params
26.314    Total estimated model params size (MB)
396       Modules in train mode
0         Modules in eval mode
Epoch 49: 100% 7/7 [00:00<00:00,  7.67it/s, train_loss=0.00264]
`Trainer.fit` stopped: `max_epochs=50` reached.

然后生成预测时,必须再次将协变量作为 past_covariates 提供给 predict() 函数。

preds = model.predict(
    series=[train_air_scaled, train_milk_scaled],
    past_covariates=[air_covs_scaled, milk_covs_scaled],
    n=36,
)

## 逆缩放:
pred_air, pred_milk = scaler.inverse_transform(preds)

plt.figure(figsize=(10, 6))
series_air.plot(label="actual (air)")
series_milk.plot(label="actual (milk)")
pred_air.plot(label="forecast (air)")
pred_milk.plot(label="forecast (milk)")
`predict()` was called with `n > output_chunk_length`: using auto-regression to forecast the values after `output_chunk_length` points. The model will access `(n - output_chunk_length)` future values of your `past_covariates` (relative to the first predicted time step). To hide this warning, set `show_warnings=False`.
GPU available: True (mps), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
Predicting DataLoader 0: 100% 1/1 [00:00<00:00, 19.90it/s]

img

现在看来,该模型可以更好地捕捉空气系列的趋势(这也对牛奶系列的预测产生了一些干扰)。

模型警告非常重要。它告知我们预测时使用了协变量的未来值,因为我们选择了范围 n > output_chunk_length 从而激活了自回归。然后模型将消耗自身的预测作为下一个预测的输入。由于预测点前移,模型需要来自过去协变量的新值。然而,这些值仅从预测点之前的过去提取(永远不会来自预测范围本身)!

在我们的案例中这没问题,因为我们仅使用提前已知的日历信息。如果使用仅在过去已知的特征,则应仅使用 n <= output_chunk_length 进行预测。

编码器:免费使用协变量

使用与日历或时间轴相关的协变量(如上述示例中的月份和年份)非常频繁,因此所有带协变量的 Darts 模型都内置了自动生成它们的功能。

为轻松将此类协变量集成到模型中,只需在模型创建时指定 add_encoders 参数。该参数必须是包含关于应编码为额外协变量的信息的字典。以下是此类字典的示例,适用于同时支持过去和未来协变量的模型:

def extract_year(idx):
    """提取每个时间索引条目的年份并归一化。"""
    return (idx.year - 1950) / 50


encoders = {
    "cyclic": {"future": ["month"]},
    "datetime_attribute": {"future": ["hour", "dayofweek"]},
    "position": {"past": ["absolute"], "future": ["relative"]},
    "custom": {"past": [extract_year]},
    "transformer": Scaler(),
}

在上面的字典中,指定了以下内容:

  • 月份应用作未来协变量,使用循环(sin/cos)编码。
  • 小时和星期几应用作未来协变量。
  • 绝对位置(序列中的时间步)应用作过去协变量。
  • 相对位置(相对于预测时间)应用作未来协变量。
  • 年份的额外自定义函数应用作过去协变量。
  • 所有上述协变量应使用 Scaler 缩放,该缩放器将在调用模型 fit() 函数时拟合,并在之后用于变换协变量。

有关如何使用编码器的更多信息,请参阅 API 文档。注意不能使用 lambda 函数,因为它们不可 pickle。

要在 N-BEATS 中使用月份和年份作为过去协变量复制示例,可以如下使用一些编码器:

encoders = {"datetime_attribute": {"past": ["month", "year"]}, "transformer": Scaler()}

现在,使用这些协变量训练 N-BEATS 模型的整个过程如下:

model = NBEATSModel(
    input_chunk_length=24,
    output_chunk_length=12,
    add_encoders=encoders,
    random_state=42,
)

model.fit([train_air_scaled, train_milk_scaled], epochs=50, verbose=True);
GPU available: True (mps), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs

  | Name            | Type             | Params | Mode
-------------------------------------------------------------
0 | criterion       | MSELoss          | 0      | train
1 | train_criterion | MSELoss          | 0      | train
2 | val_criterion   | MSELoss          | 0      | train
3 | train_metrics   | MetricCollection | 0      | train
4 | val_metrics     | MetricCollection | 0      | train
5 | stacks          | ModuleList       | 6.6 M  | train
-------------------------------------------------------------
6.6 M     Trainable params
1.7 K     Non-trainable params
6.6 M     Total params
26.314    Total estimated model params size (MB)
396       Modules in train mode
0         Modules in eval mode
Epoch 49: 100% 7/7 [00:00<00:00,  7.63it/s, train_loss=0.000928]
`Trainer.fit` stopped: `max_epochs=50` reached.

获取航空乘客序列的一些预测:

preds = model.predict(
    series=[train_air_scaled, train_milk_scaled], n=36, show_warnings=False
)

## 逆缩放:
pred_air, pred_milk = scaler.inverse_transform(preds)

plt.figure(figsize=(10, 6))
series_air.plot(label="actual (air)")
series_milk.plot(label="actual (milk)")
pred_air.plot(label="forecast (air)")
pred_milk.plot(label="forecast (milk)")

img

回归预测模型

注意:SKLearnModel 的详细示例笔记本见此处

SKLearnModel 是包装 sklearn 兼容回归模型的预测模型。内部回归模型用于预测目标序列的未来值,作为目标、过去和未来协变量的某些滞后的函数。在后台,时间序列被表格化以构建格式正确的训练数据集。

默认情况下,SKLearnModel 执行线性回归。通过指定 model 参数可以轻松使用任何所需的 sklearn 兼容回归模型,但为方便起见,Darts 还提供了一些开箱即用的现成模型:

例如,以下是将贝叶斯岭回归拟合到玩具双序列问题的示例:

from sklearn.linear_model import BayesianRidge

from darts.models import SKLearnModel

model = SKLearnModel(lags=72, lags_future_covariates=[-6, 0], model=BayesianRidge())

model.fit(
    [train_air_scaled, train_milk_scaled],
    future_covariates=[air_covs_scaled, milk_covs_scaled],
);

上述发生了以下几件事:

  • lags=72 告知 SKLearnModel 查看目标的过去 72 个滞后/月份。
  • 此外,lags_future_covariates=[-6, 0] 意味着模型还将查看提供的 future_covariates 的滞后。这里我们指定模型应考虑的确切滞后;“-6th” 和 “0th” 滞后。“0th” 滞后表示"当前"滞后(即预测的时间步);显然,了解此滞后需要提前知道数据(因此我们使用 future_covariates)。类似地,-6 表示我们还查看预测时间步之前 6 个月的协变量值(如果预测范围超过 6 步,这也需要提前知道协变量)。
  • model=BayesianRidge() 提供实际的内部回归模型。

现在获取一些预测:

preds = model.predict(
    series=[train_air_scaled, train_milk_scaled],
    future_covariates=[air_covs_scaled, milk_covs_scaled],
    n=36,
)

## 逆缩放:
pred_air, pred_milk = scaler.inverse_transform(preds)

plt.figure(figsize=(10, 6))
series_air.plot(label="actual (air)")
series_milk.plot(label="actual (milk)")
pred_air.plot(label="forecast (air)")
pred_milk.plot(label="forecast (milk)")

img

类似地,我们还可以在序列序列上获取一些指标:

mape([series_air, series_milk], [pred_air, pred_milk])
[3.435389, 5.129032]

或"所有"序列的平均指标:

mape([series_air, series_milk], [pred_air, pred_milk], series_reduction=np.mean)
4.2822104

该模型在航空交通序列上表现良好,当在此序列上回测时表现如何?

bayes_ridge_model = SKLearnModel(
    lags=72, lags_future_covariates=[0], model=BayesianRidge()
)

backtest = bayes_ridge_model.historical_forecasts(
    future_covariates=[air_covs_scaled, milk_covs_scaled],
    **hfc_params,
)

print(f"MAPE = {mape(series_air, backtest):.2f}%")
series_air.plot()
backtest.plot();
`enable_optimization=True` is ignored because `retrain` is not `False` or `0`. To hide this warning, set `show_warnings=False` or `enable_optimization=False`.
`enable_optimization=True` is ignored because `forecast_horizon > model.output_chunk_length`. To hide this warning, set `show_warnings=False` or `enable_optimization=False`.
historical forecasts: 100% 58/58 [00:00<00:00, 210.78it/s]
MAPE = 3.76%

img

我们迄今为止最好的模型!

为避免通过目标和未来协变量缩放器导致数据泄漏,我们将其传递给 historical_forecasts,以便在每次迭代时拟合并应用于目标和未来协变量。

backtest_auto_scaling = bayes_ridge_model.historical_forecasts(
    future_covariates=air_covs,
    data_transformers={"series": Scaler(), "future_covariates": Scaler()},
    retrain=True,  ## 确保每次迭代时缩放器被拟合
    **hfc_params,
)

print(f"MAPE = {mape(series_air, backtest_auto_scaling):.2f}%")
series_air.plot()
backtest_auto_scaling.plot();
`enable_optimization=True` is ignored because `retrain` is not `False` or `0`. To hide this warning, set `show_warnings=False` or `enable_optimization=False`.
`enable_optimization=True` is ignored because `forecast_horizon > model.output_chunk_length`. To hide this warning, set `show_warnings=False` or `enable_optimization=False`.
historical forecasts: 100% 58/58 [00:00<00:00, 114.27it/s]
MAPE = 3.92%

img

MAPE 分数略差,这是预期的,因为信息不再从未来泄漏,但回测条件更可靠。

样本权重

所有全局模型(回归、集成和神经网络模型)都支持训练样本权重。它们按观察值、标签(output_chunk_length 中的每一步)和分量应用。这非常有用:

  • 加权某些时间框架(例如,过去越远的点权重衰减,业务价值更高的时间权重更高)
  • 减少/消除异常值或缺失值对模型的影响

让我们在航空乘客序列中引入一些异常值,看看具有 12 个月回看窗口(lags=12)的线性回归模型在此数据上的表现。

## 转换为 pandas DataFrame 并添加异常值
air_df = series_air.to_dataframe()
outlier_mask = (air_df.index.year >= 1950) & (air_df.index.year <= 1951)
air_df.loc[outlier_mask] = 600.0
air_outlier = TimeSeries.from_dataframe(air_df)

from darts.models import LinearRegressionModel

model_lin = LinearRegressionModel(lags=12, lags_future_covariates=[0])
hist_fc_kwargs = {
    "series": air_outlier,
    "future_covariates": air_covs_scaled,
    "start": 0.6,
    "forecast_horizon": 3,
    "verbose": True,
    "show_warnings": False,
}
backtest = model_lin.historical_forecasts(**hist_fc_kwargs)

print(f"MAPE = {mape(air_outlier, backtest):.2f}%")
air_outlier.plot(label="air series with outliers")
backtest.plot(label="non-weighted predictions")
historical forecasts: 100% 57/57 [00:00<00:00, 233.79it/s]
MAPE = 26.40%

img

如果能以某种方式忽略异常值就太好了。这就是样本权重的用武之地!

在 Darts 中,样本权重被视为协变量 - 它们定义为 TimeSeries 本身。因此我们的模型可以自动为我们提取相关时间范围!

让我们创建一个权重序列,将所有异常时间设为 0.,其余时间设为 1.。我们还将最后一个异常值之后的 12 个月设为 0.,因为这些时间属于模型回看窗口。否则模型中仍会有一些异常值。

weight_mask = (air_df.index.year >= 1950) & (air_df.index.year <= 1952)
sample_weight = np.ones((len(air_df), 1))
sample_weight[weight_mask, 0] = 0.0
sample_weight = air_outlier.with_values(sample_weight)

## 绘制结果
air_outlier.plot(label="air series with outliers")
(sample_weight * 300).plot(label="sample weight * 300")

img

现在我们可以使用 fit()historical_forecastsbacktest() 等方法训练带样本权重的模型。

model_lin = LinearRegressionModel(lags=12, lags_future_covariates=[0])
backtest_weighted = model_lin.historical_forecasts(
    sample_weight=sample_weight, **hist_fc_kwargs
)

print(f"MAPE = {mape(series_air, backtest_weighted):.2f}%")
air_outlier.plot(label="air series with outliers")
backtest_weighted.plot(label="weighted predictions")
historical forecasts: 100% 57/57 [00:00<00:00, 321.96it/s]
MAPE = 5.19%

img

我们成功消除了异常值对模型的负面影响!

注意: 我们的样本权重还支持多步预测、多个时间序列、多元序列以及评估集的权重(如果模型支持)。

预测起点偏移

我们可能还对在目标序列结束后有偏移的预测感兴趣。例如:

  • 在(日)前市场中,(每天)我们需要为未来某个时间点(次日)出价 - 我们对当前到该未来时间点之间发生的事情并不真正感兴趣。通过仅关注感兴趣的时间可以减少模型复杂性。
  • 当协变量(或目标序列)延迟报告时

所有全局模型都支持通过模型创建参数 output_chunk_shift 进行此类带偏移的预测 - 将第一个预测步长偏移到未来的步数。

带输出偏移时,模型无法再执行自回归。因此让我们创建一个线性模型,直接预测未来 12 个月,偏移 12 个月。

model_shifted = LinearRegressionModel(
    lags=12,
    lags_future_covariates=(0, 12),
    output_chunk_length=12,
    output_chunk_shift=12,
)

model_shifted.fit(series_air[:-24], future_covariates=air_covs)
preds = model_shifted.predict(n=12)

series_air[:-24].plot(label="train series")
series_air[-24:].plot(label="val_series")
preds.plot(label="shifted prediction")

img

我们看到预测在训练序列结束后 12 个月开始。

概率预测

Darts 中的大多数模型都支持概率预测(一些局部模型及所有回归、集成和神经网络模型)。完整支持列表见 Darts README 页面

使用局部模型

对于局部模型(ARIMAExponentialSmoothing 等),调用 predict() 时只需指定 num_samples 参数。返回的 TimeSeries 将包含 num_samples 个蒙特卡洛样本,描述时间序列值的分布。依赖蒙特卡洛样本(与显式置信区间相比)的优势在于它们可用于描述分量上的任何参数或非参数联合分布,并计算任意分位数。

model_es = ExponentialSmoothing()
model_es.fit(train)
probabilistic_forecast = model_es.predict(len(val), num_samples=500)

series.plot(label="actual")
probabilistic_forecast.plot(label="probabilistic forecast")
plt.legend()
plt.show()

img

使用神经网络模型 (TorchForecastingModel)

使用神经网络时,必须在模型创建时使用 Likelihood 对象。似然指定模型将尝试拟合的分布,以及分布参数的潜在先验值。可用似然的完整列表见文档

除了分布外,我们还支持 QuantileRegression 作为似然,直接估计目标序列的未来分位数值。

使用似然很简单。例如,训练 TCNModel 以拟合拉普拉斯似然如下:

from darts.models import TCNModel
from darts.utils.likelihood_models.torch import LaplaceLikelihood

model = TCNModel(
    input_chunk_length=24,
    output_chunk_length=12,
    random_state=42,
    likelihood=LaplaceLikelihood(),
    pl_trainer_kwargs={"accelerator": "cpu"},
)

model.fit(train_air_scaled, epochs=400, verbose=True);
GPU available: True (mps), used: False
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
/Users/dennisbader/miniconda3/envs/darts/lib/python3.10/site-packages/pytorch_lightning/trainer/setup.py:177: GPU available but not used. You can set it by doing `Trainer(accelerator='gpu')`.

  | Name            | Type             | Params | Mode
-------------------------------------------------------------
0 | criterion       | MSELoss          | 0      | train
1 | train_criterion | MSELoss          | 0      | train
2 | val_criterion   | MSELoss          | 0      | train
3 | train_metrics   | MetricCollection | 0      | train
4 | val_metrics     | MetricCollection | 0      | train
5 | res_blocks      | ModuleList       | 166    | train
-------------------------------------------------------------
166       Trainable params
0         Non-trainable params
166       Total params
0.001     Total estimated model params size (MB)
23        Modules in train mode
0         Modules in eval mode
Epoch 399: 100% 3/3 [00:00<00:00, 188.73it/s, train_loss=-1.71]
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.

Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
概率样本预测

然后要获取概率预测,我们只需指定一些 num_samples >> 1。这将从预测分布中采样。

pred = model.predict(n=36, num_samples=500)

## 逆缩放:
pred = scaler.inverse_transform(pred, series_idx=0)

series_air.plot()
pred.plot();
GPU available: True (mps), used: False
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
Predicting DataLoader 0: 100% 1/1 [00:00<00:00, 12.61it/s]

img

直接参数预测

我们似然的强大之处在于,我们可以直接预测分布/分位数参数,而不是从分布/分位数中采样。只需设置 predict_likelihood_parameters=True

下面我们获取拉普拉斯分布的位置(mu)和尺度(b)预测(在 0 到 1 的缩放空间中)。

pred = model.predict(n=12, predict_likelihood_parameters=True)

train_air_scaled.plot()
pred.plot(label="laplace_dist")
GPU available: True (mps), used: False
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
/Users/dennisbader/miniconda3/envs/darts/lib/python3.10/site-packages/pytorch_lightning/trainer/setup.py:177: GPU available but not used. You can set it by doing `Trainer(accelerator='gpu')`.
Predicting DataLoader 0: 100% 1/1 [00:00<00:00, 171.41it/s]

img

此外,例如,我们可以指定先验信念:分布的尺度约为 0.1 0.1 0.1(在变换后的域中),同时仍捕捉分布的某些时间依赖性,通过指定 prior_b=.1

在后台,这将使用 Kullback-Leibler 散度项正则化训练损失。

model = TCNModel(
    input_chunk_length=24,
    output_chunk_length=12,
    random_state=42,
    likelihood=LaplaceLikelihood(prior_b=0.1),
    pl_trainer_kwargs={"accelerator": "cpu"},
)

model.fit(train_air_scaled, epochs=400, verbose=True);
GPU available: True (mps), used: False
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
/Users/dennisbader/miniconda3/envs/darts/lib/python3.10/site-packages/pytorch_lightning/trainer/setup.py:177: GPU available but not used. You can set it by doing `Trainer(accelerator='gpu')`.

  | Name            | Type             | Params | Mode
-------------------------------------------------------------
0 | criterion       | MSELoss          | 0      | train
1 | train_criterion | MSELoss          | 0      | train
2 | val_criterion   | MSELoss          | 0      | train
3 | train_metrics   | MetricCollection | 0      | train
4 | val_metrics     | MetricCollection | 0      | train
5 | res_blocks      | ModuleList       | 166    | train
-------------------------------------------------------------
166       Trainable params
0         Non-trainable params
166       Total params
0.001     Total estimated model params size (MB)
23        Modules in train mode
0         Modules in eval mode
Epoch 399: 100% 3/3 [00:00<00:00, 167.12it/s, train_loss=-1.38]
`Trainer.fit` stopped: `max_epochs=400` reached.
pred = model.predict(n=36, num_samples=500)

## 逆缩放:
pred = scaler.inverse_transform(pred, series_idx=0)

series_air.plot()
pred.plot();
GPU available: True (mps), used: False
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
Predicting DataLoader 0: 100% 1/1 [00:00<00:00, 12.68it/s]

img

默认情况下,TimeSeries.plot() 显示中位数以及第 5 和第 95 百分位数(如果是多元 TimeSeries,则为边际分布)。可以控制此行为:

pred.plot(low_quantile=0.01, high_quantile=0.99, label="1-99th percentiles")
pred.plot(low_quantile=0.2, high_quantile=0.8, label="20-80th percentiles")

img

分布类型

似然必须与时间序列值的域兼容。例如,PoissonLikelihood 可用于离散正值,ExponentialLikelihood 可用于实正值,BetaLikelihood 用于 ( 0 , 1 ) (0,1) (0,1) 中的实值。

也可以使用 QuantileRegression 应用分位数损失并直接拟合所需分位数。

评估概率预测

如何评估概率预测的质量?默认情况下,大多数指标函数(如 mape())仅查看中位数预测。也可以使用平均分位数损失指标 mql(),它量化每个预测分位数的误差。对于分位数=0.5(中位数),它与平均绝对误差(MAE)相同:

from darts.metrics import mae, mql

print(f"MAPE of median forecast: {mape(series_air, pred):.2f}")
print(f"MAE of median forecast: {mae(series_air, pred):.2f}")
for q in [0.05, 0.1, 0.5, 0.9, 0.95]:
    q_loss = mql(series_air, pred, q=q)
    print(f"quantile loss at quantile {q:.2f}: {q_loss:.2f}")
MAPE of median forecast: 12.14
MAE of median forecast: 51.64
quantile loss at quantile 0.05: 5.33
quantile loss at quantile 0.10: 13.14
quantile loss at quantile 0.50: 51.64
quantile loss at quantile 0.90: 20.74
quantile loss at quantile 0.95: 12.20
使用分位数损失

通过直接拟合这些分位数能做得更好吗?我们可以使用 QuantileRegression 似然:

from darts.utils.likelihood_models.torch import QuantileRegression

model = TCNModel(
    input_chunk_length=24,
    output_chunk_length=12,
    random_state=42,
    likelihood=QuantileRegression([0.05, 0.1, 0.5, 0.9, 0.95]),
)

model.fit(train_air_scaled, epochs=400, verbose=True);
GPU available: True (mps), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs

  | Name            | Type             | Params | Mode
-------------------------------------------------------------
0 | criterion       | MSELoss          | 0      | train
1 | train_criterion | MSELoss          | 0      | train
2 | val_criterion   | MSELoss          | 0      | train
3 | train_metrics   | MetricCollection | 0      | train
4 | val_metrics     | MetricCollection | 0      | train
5 | res_blocks      | ModuleList       | 208    | train
-------------------------------------------------------------
208       Trainable params
0         Non-trainable params
208       Total params
0.001     Total estimated model params size (MB)
23        Modules in train mode
0         Modules in eval mode
Epoch 399: 100% 3/3 [00:00<00:00, 102.18it/s, train_loss=0.0563]
`Trainer.fit` stopped: `max_epochs=400` reached.
pred = model.predict(n=36, num_samples=500)

## 逆缩放:
pred = scaler.inverse_transform(pred, series_idx=0)

series_air.plot()
pred.plot()

print(f"MAPE of median forecast: {mape(series_air, pred):.2f}")
for q in [0.05, 0.1, 0.5, 0.9, 0.95]:
    q_loss = mql(series_air, pred, q=q)
    print(f"quantile loss at quantile {q:.2f}: {q_loss:.2f}")
GPU available: True (mps), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
Predicting DataLoader 0: 100% 1/1 [00:00<00:00,  3.28it/s]
MAPE of median forecast: 4.91
quantile loss at quantile 0.05: 9.05
quantile loss at quantile 0.10: 13.50
quantile loss at quantile 0.50: 20.07
quantile loss at quantile 0.90: 11.99
quantile loss at quantile 0.95: 7.39

img

使用回归模型

SKLearnModels 的概率支持与神经网络类似。必须在模型创建时指定 likelihood
我们可以简单选择 "quantile"(带一些 quantiles)或 "poisson",而不是给出似然对象。

model = LinearRegressionModel(
    lags=24,
    lags_future_covariates=[0],
    likelihood="quantile",
    quantiles=[0.05, 0.1, 0.5, 0.9, 0.95],
)
model.fit(train_air, future_covariates=air_covs)
pred = model.predict(n=36, num_samples=500)

series_air.plot()
pred.plot()

print(f"MAPE of median forecast: {mape(series_air, pred):.2f}")
for q in [0.05, 0.1, 0.5, 0.9, 0.95]:
    q_loss = mql(series_air, pred, q=q)
    print(f"quantile loss at quantile {q:.2f}: {q_loss:.2f}")
MAPE of median forecast: 8.52
quantile loss at quantile 0.05: 20.76
quantile loss at quantile 0.10: 25.44
quantile loss at quantile 0.50: 35.33
quantile loss at quantile 0.90: 11.15
quantile loss at quantile 0.95: 6.19

img

模型集成

集成是将多个模型产生的预测组合以获得最终(且希望更好的)预测。

例如,在上面的改进模型示例中,我们手动组合了朴素季节性模型和朴素漂移模型。这里,我们将展示如何自动组合模型预测 - 朴素地使用 NaiveEnsembleModel,或学习使用 RegressionEnsembleModel

当然也可以使用 past 和/或 future_covariates 与集成模型,但它们仅传递给支持它们的预测模型。

朴素集成

朴素集成仅对多个模型的预测取平均。Darts 的 NaiveEnsembleModel 正是这样做的,并通过与预测模型相同的 API(拟合、预测、历史预测、回测等)实现。

from darts.models import NaiveEnsembleModel

models = [NaiveDrift(), NaiveSeasonal(12)]

ensemble_model = NaiveEnsembleModel(forecasting_models=models)

backtest = ensemble_model.historical_forecasts(**hfc_params)

print(f"MAPE = {mape(backtest, series_air):.2f}%")
series_air.plot()
backtest.plot();
historical forecasts: 100% 58/58 [00:00<00:00, 194.96it/s]
MAPE = 11.93%

img

学习集成

在这种情况下,朴素集成效果不佳(但在某些情况下可能不错!)

如果我们将集成视为监督回归问题,有时可以做得更好:给定一组预测(特征),找到一个模型将它们组合以最小化目标误差。
这就是 RegressionEnsembleModel 所做的。它接受参数:

  • forecasting_models 是我们想要集成其预测的预测模型列表。
  • regression_train_n_points 是用于拟合"集成回归"模型(即组合预测的内部模型)的时间步数。
  • regression_model 是可选的,用于集成回归的 sklearn 兼容回归模型或 Darts SKLearnModel。如果未指定,则使用线性回归。使用 sklearn 模型很容易,但使用 SKLearnModel 允许将各个预测的任意滞后作为回归模型的输入。
  • 更多信息请阅读我们的集成模型用户指南

一旦这些元素就位,RegressionEnsembleModel 可以像常规预测模型一样使用:

from darts.models import RegressionEnsembleModel

ensemble_model = RegressionEnsembleModel(
    forecasting_models=models, regression_train_n_points=12
)

backtest = ensemble_model.historical_forecasts(**hfc_params)

print(f"MAPE = {mape(backtest, series_air):.2f}")
series_air.plot()
backtest.plot();
historical forecasts: 100% 58/58 [00:00<00:00, 137.11it/s]
MAPE = 4.77

img

我们还可以检查线性组合中两个内部模型的权重系数:

ensemble_model.fit(series_air)
ensemble_model.regression_model.model.coef_
array([0.01368814, 1.0980107 ], dtype=float32)

集成模型本身也可以是概率的!可在集成模型指南中阅读相关内容。

RegressionEnsembleModel 使用堆叠技术训练和组合 forecasting_models:每个模型独立训练,然后使用它们的预测作为 future_covariates 训练 regression_model

滤波模型

除了能够预测序列未来值的预测模型外,Darts 还包含几个有用的滤波模型,可对"样本内"序列值分布建模。

拟合卡尔曼滤波器

KalmanFilter 实现卡尔曼滤波器。该实现依赖于 nfoursid,因此可以提供一个包含转移矩阵、过程噪声协方差、观测噪声协方差等的 nfoursid.kalman.Kalman 对象。

也可以通过调用 fit() 进行系统辨识,使用 N4SID 系统辨识算法"训练"卡尔曼滤波器:

from darts.models import KalmanFilter

kf = KalmanFilter(dim_x=3)
kf.fit(train_air_scaled)
filtered_series = kf.filter(train_air_scaled, num_samples=100)

train_air_scaled.plot()
filtered_series.plot();

img

使用高斯过程推断缺失值

Darts 还包含 GaussianProcessFilter,可用于序列的概率建模:

from sklearn.gaussian_process.kernels import RBF

from darts.models import GaussianProcessFilter

## 创建带缺失值的序列:
values = train_air_scaled.values()
values[20:22] = np.nan
values[28:32] = np.nan
values[55:59] = np.nan
values[72:80] = np.nan
series_holes = TimeSeries.from_times_and_values(train_air_scaled.time_index, values)
series_holes.plot()

kernel = RBF()

gpf = GaussianProcessFilter(kernel=kernel, alpha=0.1, normalize_y=True)
filtered_series = gpf.filter(series_holes, num_samples=100)

filtered_series.plot();

img

注意事项

那么 N-BEATS、指数平滑还是在牛奶产量上训练的贝叶斯岭回归是预测未来航空乘客数量的最佳方法吗?嗯,此时很难确切地说出哪个最好。我们的时间序列很小,验证集更小。在这种情况下,很容易将整个预测过程过度拟合到如此小的验证集上。尤其是当可用模型及其自由度很高时(如深度学习模型),或者如果我们在单个测试集上尝试了许多模型(如本笔记本中所做)。

作为数据科学家,我们有责任了解模型的可信程度。因此,请谨慎对待结果,尤其是在小型数据集上,并在进行任何预测前应用科学方法 😃 祝建模愉快!

异常检测

Darts 还有一个完整的时间序列异常检测模块。更多信息请阅读我们的异常检测用户指南

风险提示与免责声明
本文内容基于公开信息研究整理,不构成任何形式的投资建议。历史表现不应作为未来收益保证,市场存在不可预见的波动风险。投资者需结合自身财务状况及风险承受能力独立决策,并自行承担交易结果。作者及发布方不对任何依据本文操作导致的损失承担法律责任。市场有风险,投资须谨慎。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

船长Q

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值