Autoregressive Denoising Diffusion Models for Multivariate Probabilistic Time Series Forecasting

Introduction

Classical time series forecasting methods such as those in (Hyndman & Athanasopoulos, 2018) typically provide univariate point forecasts, require hand-tuned features to model seasonality, and are trained individually on each time series. Deep learning based time series models (Benidis et al., 2020) are popular alternatives due to their end-to-end training of a global model, ease of incorporating exogenous covariates, and automatic feature extraction abilities. The task of modeling uncertainties is of vital importance for downstream problems that use these forecasts for (business) decision making. More often the individual time series for a problem data set are statistically dependent on each other. Ideally, deep learning models need to incorporate this inductive bias in the form of multivariate (Tsay, 2014) probabilistic methods to provide accurate forecasts.

To model the full predictive distribution, methods typically resort to tractable distribution classes or some type of low-rank approximations, regardless of the true data distribution. To model the distribution in a general fashion, one needs probabilistic methods with tractable likelihoods. Till now several deep learning methods have been proposed for this purpose such as autoregressive (van den Oord et al., 2016c) or generative ones based on normalizing flows (Papamakarios et al., 2019) which can learn flexible models of high dimensional multivariate time series. Even if the full likelihood is not be tractable, one can often optimize a tractable lower bound to the likelihood. But still, these methods require a certain structure in the functional approximators, for example on the determinant of the Jacobian (Dinh et al., 2017) for normalizing flows. Energy-based models (EBM) (Hinton, 2002; LeCun et al., 2006) on the other hand have a much less restrictive functional form. They approximate the unnormalized log-probability so that density estimation reduces to a non-linear regression problem. EBMs have been shown to perform well in learning high dimensional data distributions at the cost of being difficult to train (Song & Kingma, 2021).

In this work, we propose autoregressive EBMs to solve the multivariate probabilistic time series forecasting problem via a model we call TimeGrad and show that not only are we able to train such a model with all the inductive biases of probabilistic time series forecasting, but this model performs exceptionally well when compared to other modern methods. This autoregressive-EBM combination retains the power of autoregressive models, such as good performance in extrapolation into the future, with the flexibility of EBMs as a general purpose high-dimensional distribution model, while remaining computationally tractable.

The paper is organized as follows. In Section 2 we first set up the notation and detail the EBM of (Ho et al., 2020) which forms the basis of our per time-step distribution model. Section 3 introduces the multivariate probabilistic time series problem and we detail the TimeGrad model. The experiments with extensive results are detailed in Section 4. We cover related work in Section 5 and conclude with some discussion in Section 6.

Diffusion Probabilistic Model

is not trainable but fixed to a Markov chain (called the forward process) that gradually adds Gaussian noise to the signal:

The forward process uses an increasing variance schedule $\beta_{1},\ldots,\beta_{N}$ with $\beta_{n}\in(0,1)$. The joint distribution $p_{\theta}(\mathbf{x}^{0:N})$ is called the reverse process, and is defined as a Markov chain with learned Gaussian transitions starting with $p(\mathbf{x}^{N})=\mathcal{N}(\mathbf{x}^{N};\mathbf{0},\mathbf{I})$, where each subsequent transition of

is given by a parametrization of our choosing denoted by

This upper bound can be shown to be equal to

As shown by (Ho et al., 2020), a property of the forward process is that it admits sampling $\mathbf{x}^{n}$ at any arbitrary noise level $n$ in closed form, since if $\alpha_{n}:=1-\beta_{n}$ and $\bar{\alpha}_{n}:=\Pi_{i=1}^{n}\alpha_{i}$ its cumulative product, we have:

By using the fact that these processes are Markov chains, the objective in (2) can be written as the KL-divergence between Gaussian distributions:

and (Ho et al., 2020) shows that by the property (3) the forward process posterior in these KL divergences when conditioned on $\mathbf{x}^{0}$, i.e. $q(\mathbf{x}^{n-1}|\mathbf{x}^{n},\mathbf{x}^{0})$ are tractable given by

Further, (Ho et al., 2020) shows that the KL-divergence between Gaussians can be written as:

where $\mathbf{\epsilon}_{\theta}$ is a network which predicts $\mathbf{\epsilon}\sim{\cal{N}}(\mathbf{0},\mathbf{I})$ from $\mathbf{x}^{n}$, so that the objective simplifies to:

resembling the loss in Noise Conditional Score Networks (Song & Ermon, 2019, 2020) using score matching. Once trained, to sample from the reverse process $\mathbf{x}^{n-1}\sim p_{\theta}(\mathbf{x}^{n-1}|\mathbf{x}^{n})$ (1) we can compute

where $\mathbf{z}\sim{\cal{N}}(\mathbf{0},\mathbf{I})$ for $n=N,\ldots,2$ and $\mathbf{z}=\mathbf{0}$ when $n=1$. The full sampling procedure for $\mathbf{x}^{0}$, starting from white noise sample $\mathbf{x}^{N}$, resembles Langevin dynamics where we sample from the most noise-perturbed distribution and reduce the magnitude of the noise scale until we reach the smallest one.

TimeGrad Method

In the univariate probabilistic DeepAR model (Salinas et al., 2019b), the log-likelihood of each entity $x^{0}_{i,t}$ at a time step $t\in[t_{0},T]$ is maximized over an individual time series’ prediction window. This is done with respect to the parameters of some chosen distributional model via the state of an RNN derived from its previous time step $x^{0}_{i,t-1}$ and its corresponding covariates $\mathbf{c}_{i,t-1}$. The emission distribution model, which is typically Gaussian for real-valued data or negative binomial for count data, is selected to best match the statistics of the time series and the network incorporates activation functions that satisfy the constraints of the distribution’s parameters, e.g. a softplus() for the scale parameter of the Gaussian.

A straightforward time series model for multivariate real-valued data could use a factorizing output distribution instead. Shared parameters can then learn patterns across the individual time series entities through the temporal component — but the model falls short of capturing dependencies in the emissions of the model. For this, a full joint distribution at each time step has to be modeled, for example by using a multivariate Gaussian. However, modeling the full covariance matrix not only increases the number of parameters of the neural network by $O(D^{2})$, making learning difficult but computing the loss is $O(D^{3})$ making it impractical. Furthermore, statistical dependencies for such distributions would be limited to second-order effects. Approximating Gaussians with low-rank covariance matrices do work however and these models are referred to as Vec-LSTM in (Salinas et al., 2019a).

Instead, in this work we propose TimeGrad which aims to learn a model of the conditional distribution of the future time steps of a multivariate time series given its past and covariates as:

were we assume that the covariates are known for all the time points and each factor is learned via a conditional denoising diffusion model introduced above. To model the temporal dynamics we employ the autoregressive recurrent neural network (RNN) architecture from (Graves, 2013; Sutskever et al., 2014) which utilizes the LSTM (Hochreiter & Schmidhuber, 1997) or GRU (Chung et al., 2014) to encode the time series sequence up to time point $t$, given the covariates $\mathbf{c}_{t}$, via the updated hidden state $\mathbf{h}_{t}$:

where now $\theta$ comprises the weights of the RNN as well as denoising diffusion model. This model is autoregressive as it consumes the observations at the time step $t-1$ as input to learn the distribution of, or sample, the next time step as shown in Figure 1.

Training is performed by randomly sampling context and adjoining prediction sized windows from the training time series data and optimizing the parameters $\theta$ that minimize the negative log-likelihood of the model (10):

starting with the hidden state $\mathbf{h}_{t_{0}-1}$ obtained by running the RNN on the chosen context window. Via a similar derivation as in the previous section, we have that the conditional variant of the objective (4) for time step $t$ and noise index $n$ is given by the following simplification of (7) (Ho et al., 2020):

2 Inference

After training, we wish to predict for each time series in our data set some prediction steps into the future and compare with the corresponding test set time series. As in training, we run the RNN over the last context sized window of the training set to obtain the hidden state $\mathbf{h}_{T}$ via (9). Then we follow the sampling procedure in Algorithm 2 to obtain a sample $\mathbf{x}_{T+1}^{0}$ of the next time step, which we can pass autoregressively to the RNN together with the covariates $\mathbf{c}_{T+1}$ to obtain the next hidden state $\mathbf{h}_{T+1}$ and repeat until the desired forecast horizon has been reached. This process of sampling trajectories from the “warm-up” state $\mathbf{h}_{T}$ can be repeated many times (e.g. $S=100$) to obtain empirical quantiles of the uncertainty of our predictions.

3 Scaling

In real-world data, the magnitudes of different time series entities can vary drastically. To normalize scales, we divide each time series entity by their context window mean (or $1$ if it’s zero) before feeding it into the model. At inference, the samples are then multiplied by the same mean values to match the original scale. This rescaling technique simplifies the problem for the model, which is reflected in significantly improved empirical performance as shown in (Salinas et al., 2019b). The other method of a short-cut connection from the input to the output of the function approximator, as done in the multivariate point forecasting method LSTNet (Lai et al., 2018), is not applicable here.

4 Covariates

We employ embeddings for categorical features (Charrington, 2018), that allows for relationships within a category, or its context, to be captured when training time series models. Combining these embeddings as features for forecasting yields powerful models like the first place winner of the Kaggle Taxi Trajectory Predictionhttps://www.kaggle.com/c/pkdd-15-predict-taxi-service-trajectory-i challenge (De Brébisson et al., 2015). The covariates $\mathbf{c}_{t}$ we use are composed of time-dependent (e.g. day of week, hour of day) and time-independent embeddings, if applicable, as well as lag features depending on the time frequency of the data set we are training on. All covariates are thus known for the periods we wish to forecast.

Experiments

We benchmark TimeGrad on six real-world data sets and evaluate against several competitive baselines. The source code of the model will be made available after the review process.

For our experiments we use Exchange (Lai et al., 2018), Solar (Lai et al., 2018), Electricityhttps://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014, Traffichttps://archive.ics.uci.edu/ml/datasets/PEMS-SF, Taxihttps://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page and Wikipediahttps://github.com/mbohlkeschneider/gluon-ts/tree/mv_release/datasets open data sets, preprocessed exactly as in (Salinas et al., 2019a), with their properties listed in Table 1. As can be noted in the table, we do not need to normalize scales for Traffic.

2 Model Architecture

We train TimeGrad via SGD using Adam (Kingma & Ba, 2015) with learning rate of $1\text{\times}{10}^{-3}$ on the training split of each data set with $N=100$ diffusion steps using a linear variance schedule starting from $\beta_{1}=$1\text{\times}{10}^{-4}$$till$\beta_{N}=0.1$. We construct batches of size$64$by taking random windows (with possible overlaps), with the context size set to the number of prediction steps, from the total time steps of each data set (see Table 1). For testing we use a rolling windows prediction starting from the last context window history before the start of the prediction and compare it to the ground-truth in the test set by sampling$S=100$ trajectories.

All experiments run on a single Nvidia V100 GPU with $16$GB of memory.

3 Results

VAR (Lütkepohl, 2007) a mutlivariate linear vector auto-regressive model with lags corresponding to the periodicity of the data,

GARCH (van der Weide, 2002) a multivariate conditional heteroskedastic model and

VES a innovation state space model (Hyndman et al., 2008);

as well as deep learning based methods namely:

KVAE (Fraccaro et al., 2017) a variational autoencoder to represent the data on top of a linear state space model which describes the dynamics,

Vec-LSTM-ind-scaling (Salinas et al., 2019a) which models the dynamics via an RNN and outputs the parameters of an independent Gaussian distribution with mean-scaling,

Vec-LSTM-lowrank-Copula (Salinas et al., 2019a) which instead parametrizes a low-rank plus diagonal covariance via Copula process,

GP-scaling (Salinas et al., 2019a) which unrolls an LSTM with scaling on each individual time series before reconstructing the joint distribution via a low-rank Gaussian,

GP-Copula (Salinas et al., 2019a) which unrolls an LSTM on each individual time series and then the joint emission distribution is given by a low-rank plus diagonal covariance Gaussian copula and

Transformer-MAF (Rasul et al., 2021) which uses Transformer (Vaswani et al., 2017) to model the temporal conditioning and Masked Autoregressive Flow (Papamakarios et al., 2017) for the distribution emission model.

4 Ablation

To highlight the predictions of TimeGrad we show in Figure 4 the predicted median, $50\%$ and $90\%$ distribution intervals of the first $6$ dimensions of the full $963$ dimensional multivariate forecast of the Traffic benchmark.

Related Work

The EBM of (Ho et al., 2020) that we adapt is based on methods that learn the gradient of the log-density with respect to the inputs, called Stein Score function (Hyvärinen, 2005; Vincent, 2011), and at inference time use this gradient estimate via Langevin dynamics to sample from the model of this complicated data distribution (Song & Ermon, 2019). These models achieve impressive results for image generation (Ho et al., 2020; Song & Ermon, 2020) when trained in an unsupervised fashion without requiring adversarial optimization. By perturbing the data using multiple noise scales, the learnt Score network captures both coarse and fine-grained data features.

The closest related work to TimeGrad is in the recent non-autoregressive conditional methods for high fidelity waveform generation (Chen et al., 2021; Kong et al., 2021). Although these methods learn the distribution of vector valued data via denoising diffusion methods, as done here, they do not consider its temporal development. Also neighboring dimensions of waveform data are highly correlated and have a uniform scale, which is not necessarily true for multivariate time series problems where neighboring entities occur arbitrarily (but in a fixed order) and can have different scales. (Du & Mordatch, 2019) also use EBMs to model one and multiple steps for a trajectory modeling task in an non-autoregressive fashion.

2 Time Series Forecasting

Neural time series methods have recently become popular ways of solving the prediction problem via univariate point forecasting methods (Oreshkin et al., 2020; Smyl, 2020) or univariate probabilistic methods (Salinas et al., 2019b). In the multivariate setting we also have point forecasting methods (Lai et al., 2018; Li et al., 2019) as well as probabilistic methods, like this method, which explicitly model the data distribution using Gaussian copulas (Salinas et al., 2019a), GANs (Yoon et al., 2019), or normalizing flows (de Bézenac et al., 2020; Rasul et al., 2021). Bayesian neural networks can also be used to provide epistemic uncertainty in forecasts as well as detect distributional shifts (Zhu & Laptev, 2018), although these methods often do not perform as well empirically (Wenzel et al., 2020).

Conclusion and Future Work

We have presented TimeGrad, a versatile multivariate probabilistic time series forecasting method that leverages the exceptional performance of EBMs to learn and sample from the distribution of the next time step, autoregressivly. Analysis of TimeGrad on six commonly used time series benchmarks establishes the new state-of-the-art against competitive methods.

We note that while training TimeGrad we do not need to loop over the EBM function approximator $\epsilon_{\theta}$, unlike in the normalizing flow setting where we have multiple stacks of bijections. However while sampling we do loop $N$ times over $\epsilon_{\theta}$. A possible strategy to improve sampling times introduced in (Chen et al., 2021) uses a combination of improved variance schedule and an $L_{1}$ loss to allow sampling with fewer steps at the cost of a small reduction in quality if such a trade-off is required. A recent paper (Song et al., 2021) generalize the diffusion processes via a class of non-Markovian processes which also allows for faster sampling.

The use of normalizing flows for discrete valued data dictates that one dequantizes it (Theis et al., 2016), by adding uniform noise to the data, before using the flows to learn. Dequantization is not needed in the EBM setting and future work could explore methods of explicitly modeling discrete distributions.

As noted in (Du & Mordatch, 2019) EBMs exhibit better out-of-distribution (OOD) detection than other likelihood models. Such a task requires models to have a high likelihood on the data manifold and low at all other locations. Surprisingly (Nalisnick et al., 2019) showed that likelihood models, including flows, were assigning higher likelihoods to OOD data whereas EBMs do not suffer from this issue since they penalize high probability under the model but low probability under the data distribution explicitly. Future work could evaluate the usage of TimeGrad for anomaly detection tasks.

For long time sequences, one could replace the RNN with a Transformer architecture (Rasul et al., 2021) to provide better conditioning for the EBM emission head. Concurrently, since EBMs are not constrained by the form of their functional approximators, one natural way to improve the model would be to incorporate architectural choices that best encode the inductive bias of the problem being tackled, for example with graph neural networks (Niu et al., 2020) when the relationships between entities are known.

References