GluonTS - 概率时间序列建模(Probabilistic Time Series Modeling)

最近在研究时间序列预测模型的的研究。关于时间序列的更多介绍，知乎已经有大佬进行详细系统的分类介绍了。有兴趣的可以直接去这里看一下。

这里是关于GluonTS官方API中Quick Start Tutorial部分的源码学习，通过阅读教程案例进行相关的翻译和一些自己的心得总结，如有错误，欢迎指正。官方API案例地址

1. 快速开始向导

GluonTS工具箱包含用于使用MXNet构建时间序列模型的组件和工具。 当前包含的模型是预测模型，但组件还支持其他时间序列用例，例如分类或异常检测。
该工具包并非旨在作为企业或最终用户的预测解决方案，而是针对想要调整算法或构建和试验自己模型的科学家和工程师。
内容包括：

• 用于构建新模型的组件（释然函数，特征处理的pipelines，日期特征，等）
• 数据加载和处理
• 多种预设模型
• 绘图和评估指标
• 人工数据集和真实数据集

导入相关库：

# Third-party imports
%matplotlib inline
import mxnet as mx
from mxnet import gluon
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import json

2. 数据集datasets

GluonTS自带了许多公开的数据集，可以直接导入

from gluonts.dataset.repository.datasets import get_dataset, dataset_recipes
from gluonts.dataset.util import to_pandas

要下载其中一个内置数据集，只需使用上述名称之一调用get_dataset。 GluonTS可以重新使用保存的数据集，因此无需再次下载：只需设置regenerate = False。

dataset = get_dataset("m4_hourly", regenerate=True)

通常，GluonTS提供的数据集是由三个主要成员组成的对象：

• dataset.train 是用于训练的数据条目的可迭代集合。 每个条目对应一个时间序列
• dataset.test 是用于推理的数据条目的可迭代集合。 测试数据集是火车数据集的扩展版本，在每个时间序列的末尾包含一个在训练期间未看到的窗口。 该窗口的长度等于建议的预测长度。
• dataset.metadata 包含数据集的元数据，例如时间序列的频率，建议的预测范围，相关特征等。

# 绘制训练集数据图像
entry = next(iter(dataset.train))
train_series = to_pandas(entry)
train_series.plot()
plt.grid(which="both")
plt.legend(["train series"], loc="upper left")
plt.show()

# 绘制测试集数据图像
entry = next(iter(dataset.test))
test_series = to_pandas(entry)
test_series.plot()
plt.axvline(train_series.index[-1], color='r') # end of train dataset
plt.grid(which="both")
plt.legend(["test series", "end of train series"], loc="upper left")
plt.show()

print(f"Length of forecasting window in test dataset: {len(test_series) - len(train_series)}")
print(f"Recommended prediction horizon: {dataset.metadata.prediction_length}")
print(f"Frequency of the time series: {dataset.metadata.freq}")

3. 自定义数据集

在这一点上，必须强调的是，对于用户可能拥有的自定义数据集，GluonTS不需要这种特定格式。 自定义数据集的唯一要求是可迭代的，并具有“目标”和“开始”字段。 为了更清楚地说明这一点，请假设常见的情况是数据集采用numpy.array的形式，并且时间序列的索引以pandas.Timestamp表示（每个时间序列可能不同）：

N = 10  # number of time series
T = 100  # number of timesteps
prediction_length = 24
freq = "1H"
custom_dataset = np.random.normal(size=(N, T))
start = pd.Timestamp("01-01-2019", freq=freq)  # can be different for each time series

现在，您可以拆分数据集，并仅需两行代码即可将其转换为GluonTS适当的格式：

from gluonts.dataset.common import ListDataset

# train dataset: cut the last window of length "prediction_length", add "target" and "start" fields
train_ds = ListDataset([{'target': x, 'start': start}
                        for x in custom_dataset[:, :-prediction_length]],
                       freq=freq)
# test dataset: use the whole dataset, add "target" and "start" fields
test_ds = ListDataset([{'target': x, 'start': start}
                       for x in custom_dataset],
                      freq=freq)

4. 训练一个现有的模型(Estimator)

GluonTS带有许多预先构建的模型。 用户所需要做的就是配置一些超参数。 现有模型专注于（但不限于）概率预测。 概率预测是以概率分布的形式进行的预测，而不是简单的单点估计。
我们将从GulonTS的预先构建的前馈神经网络估计器开始，这是一个简单但功能强大的预测模型。 我们将使用此模型来演示训练模型，生成预测和评估结果的过程。
GluonTS的内置前馈神经网络（SimpleFeedForwardEstimator）接受长度为context_length的输入窗口，并预测随后的projection_length 值的分布。 用GluonTS的话来说，前馈神经网络模型是Estimator的一个示例。 在GluonTS中，Estimator对象代表一种预测模型以及诸如其系数，权重等详细信息。
通常，每个估计器（预先构建或自定义）由许多超参数配置，这些超参数可以是所有估计器之间（例如，prediction_length）相同（但不具有约束力），也可以特定于特定估计器（例如，层数） 用于神经网络或CNN中的步幅）。
最后，每个估计器都由一个训练器配置Trainer，该训练器定义如何训练模型，即时期数，学习率等。

from gluonts.model.simple_feedforward import SimpleFeedForwardEstimator
from gluonts.trainer import Trainer

estimator = SimpleFeedForwardEstimator(
    num_hidden_dimensions=[10],
    prediction_length=dataset.metadata.prediction_length,
    context_length=100,
    freq=dataset.metadata.freq,
    trainer=Trainer(ctx="cpu",
                    epochs=5,
                    learning_rate=1e-3,
                    num_batches_per_epoch=100
                   )
)

在为我们的估算器指定了所有必要的超参数之后，我们可以通过调用估算器的训练方法，使用训练数据集data.train对其进行训练。 训练算法返回可用于构建预测的拟合模型（或GluonTS术语中的Predictor)。

predictor = estimator.train(dataset.train)

现在有了预测器，我们就可以预测数据集dataset.test的最后一个窗口。测试并评估模型的性能。
GluonTS带有make_evaluation_predictions函数，该函数可自动进行预测和模型评估过程。 大致来说，此功能执行以下步骤：

• 删除我们要预测的dataset.test的最后一个length_precision_length的最终窗口
• 估计器使用剩余的数据来预测（以样本路径的形式）刚刚删除的“未来”窗口
• 模块输出预测样本路径和dataset.test（作为python生成器对象）

from gluonts.evaluation.backtest import make_evaluation_predictions

forecast_it, ts_it = make_evaluation_predictions(
    dataset=dataset.test,  # test dataset
    predictor=predictor,  # predictor
    num_samples=100,  # number of sample paths we want for evaluation
)

首先，我们可以将这些生成器转换为列表，以简化后续计算。

(ps:这里有些个人看法，个人觉得llist()函数的操作会开辟大量的内存空间，因为forecast_it本身是一个迭代器，个人感觉直接使用next(iter(forecast_it)))操作可以节省大量的内存空间，当然List可以直接选择索引号更加灵活。

forecasts = list(forecast_it)
tss = list(ts_it)

我们可以检查这些列表的第一个元素（与数据集的第一个时间序列相对应）。 让我们从包含时间序列的列表开始，即tss。 我们希望tss的第一个条目包含dataset.test的第一个时间序列（目标）。

# first entry of the time series list
ts_entry = tss[0]

# first 5 values of the time series (convert from pandas to numpy)
np.array(ts_entry[:5]).reshape(-1,)

# first entry of dataset.test
dataset_test_entry = next(iter(dataset.test))

# first 5 values
dataset_test_entry['target'][:5]

预测forecast列表中的条目要复杂一些。 它们是包含所有样本路径的对象，这些样本路径的形式为numpy.ndarray，其维度为（num_samples，prediction_length），预测的开始日期，时间序列的频率等。我们可以通过简单地调用 预测对象的相应属性。

# first entry of the forecast list
forecast_entry = forecasts[0]

print(f"Number of sample paths: {forecast_entry.num_samples}")
print(f"Dimension of samples: {forecast_entry.samples.shape}")
print(f"Start date of the forecast window: {forecast_entry.start_date}")
print(f"Frequency of the time series: {forecast_entry.freq}")

我们还可以进行计算以总结样本路径，例如计算预测窗口中48个时间步长中的每个时间步长的均值或分位数。

print(f"Mean of the future window:\n {forecast_entry.mean}")
print(f"0.5-quantile (median) of the future window:\n {forecast_entry.quantile(0.5)}")

Forecast对象具有一种plot方法，可以将预测路径总结为平均值，预测间隔等。预测间隔以不同的颜色阴影显示为“扇形图”。

def plot_prob_forecasts(ts_entry, forecast_entry):
    plot_length = 150
    prediction_intervals = (50.0, 90.0)
    legend = ["observations", "median prediction"] + [f"{k}% prediction interval" for k in prediction_intervals][::-1]

    fig, ax = plt.subplots(1, 1, figsize=(10, 7))
    ts_entry[-plot_length:].plot(ax=ax)  # plot the time series
    forecast_entry.plot(prediction_intervals=prediction_intervals, color='g')
    plt.grid(which="both")
    plt.legend(legend, loc="upper left")
    plt.show()

plot_prob_forecasts(ts_entry, forecast_entry)

我们还可以通过数字评估我们的预测质量。 在GluonTS中，Evaluator类可以计算综合性能指标以及每个时间序列的指标（这对于分析异构时间序列的性能很有用）。

from gluonts.evaluation import Evaluator

evaluator = Evaluator(quantiles=[0.1, 0.5, 0.9])
agg_metrics, item_metrics = evaluator(iter(tss), iter(forecasts), num_series=len(dataset.test))

汇总指标跨时间步长和跨时间序列汇总。

print(json.dumps(agg_metrics, indent=4))

###输出展示###
{
    "MSE": 8443908.678894352,
    "abs_error": 8877123.154870987,
    "abs_target_sum": 145558863.59960938,
    "abs_target_mean": 7324.822041043146,
    "seasonal_error": 336.9046924038305,
    "MASE": 3.246015285156499,
    "sMAPE": 0.18506043295168975,
    "MSIS": 38.217292786814824,
    "QuantileLoss[0.1]": 5679215.008755971,
    "Coverage[0.1]": 0.09983896940418678,
    "QuantileLoss[0.5]": 8877123.216818333,
    "Coverage[0.5]": 0.5041767310789049,
    "QuantileLoss[0.9]": 7424048.535117529,
    "Coverage[0.9]": 0.892109500805153,
    "RMSE": 2905.8404427797395,
    "NRMSE": 0.39671140493208645,
    "ND": 0.0609864829618992,
    "wQuantileLoss[0.1]": 0.0390166209622099,
    "wQuantileLoss[0.5]": 0.06098648338748198,
    "wQuantileLoss[0.9]": 0.051003754436685866,
    "mean_wQuantileLoss": 0.05033561959545924,
    "MAE_Coverage": 0.004076086956521706
}

5. 创建你自己的预测模型

ps : 创建自己的预测模型，需要自己构件神经网络层次结构，以及设置自己的损失函数等，灵活的构造一个自己的Estimator，由于个人实际情况，并没有进行深入的研究，仅仅是调用现有的模型接口API，对官方模型有兴趣的可以直接去官方模型API – estimator 阅读具体的参数和方法。如果有兴趣自己构建模型的也欢迎和我一起探讨研究。

下面仅仅放出官方代码，不过多进行翻译解读。

class MyTrainNetwork(gluon.HybridBlock):
    def __init__(self, prediction_length, **kwargs):
        super().__init__(**kwargs)
        self.prediction_length = prediction_length

        with self.name_scope():
            # Set up a 3 layer neural network that directly predicts the target values
            self.nn = mx.gluon.nn.HybridSequential()
            self.nn.add(mx.gluon.nn.Dense(units=40, activation='relu'))
            self.nn.add(mx.gluon.nn.Dense(units=40, activation='relu'))
            self.nn.add(mx.gluon.nn.Dense(units=self.prediction_length, activation='softrelu'))

    def hybrid_forward(self, F, past_target, future_target):
        prediction = self.nn(past_target)
        # calculate L1 loss with the future_target to learn the median
        return (prediction - future_target).abs().mean(axis=-1)


class MyPredNetwork(MyTrainNetwork):
    # The prediction network only receives past_target and returns predictions
    def hybrid_forward(self, F, past_target):
        prediction = self.nn(past_target)
        return prediction.expand_dims(axis=1)

from gluonts.model.estimator import GluonEstimator
from gluonts.model.predictor import Predictor, RepresentableBlockPredictor
from gluonts.core.component import validated
from gluonts.support.util import copy_parameters
from gluonts.transform import ExpectedNumInstanceSampler, Transformation, InstanceSplitter
from gluonts.dataset.field_names import FieldName
from mxnet.gluon import HybridBlock

class MyEstimator(GluonEstimator):
    @validated()
    def __init__(
        self,
        freq: str,
        context_length: int,
        prediction_length: int,
        trainer: Trainer = Trainer()
    ) -> None:
        super().__init__(trainer=trainer)
        self.context_length = context_length
        self.prediction_length = prediction_length
        self.freq = freq


    def create_transformation(self):
        # Feature transformation that the model uses for input.
        # Here we use a transformation that randomly select training samples from all time series.
        return InstanceSplitter(
                    target_field=FieldName.TARGET,
                    is_pad_field=FieldName.IS_PAD,
                    start_field=FieldName.START,
                    forecast_start_field=FieldName.FORECAST_START,
                    train_sampler=ExpectedNumInstanceSampler(num_instances=1),
                    past_length=self.context_length,
                    future_length=self.prediction_length,
                )

    def create_training_network(self) -> MyTrainNetwork:
        return MyTrainNetwork(
            prediction_length=self.prediction_length
        )

    def create_predictor(
        self, transformation: Transformation, trained_network: HybridBlock
    ) -> Predictor:
        prediction_network = MyPredNetwork(
            prediction_length=self.prediction_length
        )

        copy_parameters(trained_network, prediction_network)

        return RepresentableBlockPredictor(
            input_transform=transformation,
            prediction_net=prediction_network,
            batch_size=self.trainer.batch_size,
            freq=self.freq,
            prediction_length=self.prediction_length,
            ctx=self.trainer.ctx,
        )