Multivariate Probabilistic Time Series Forecasting via Conditioned Normalizing Flows
Authors: Kashif Rasul
Abdul-Saboor Sheikh
Ingmar Schuster
Urs M Bergmann
Roland Vollgraf
Venue: ICLR 2021 Spotlight (2021)
Keywords: time series
normalizing flows
attention
probabilistic multivariate forecasting

Abstract: Time series forecasting is often fundamental to scientific and engineering problems and enables decision making. With ever increasing data set sizes, a trivial solution to scale up predictions is to assume independence between interacting time series. However, modeling statistical dependencies can improve accuracy and enable analysis of interaction effects. Deep learning methods are well suited for this problem, but multi-variate models often assume a simple parametric distribution and do not scale to high dimensions. In this work we model the multi-variate temporal dynamics of time series via an autoregressive deep learning model, where the data distribution  is represented by a conditioned normalizing flow. This combination retains the power of autoregressive models, such as good performance in extrapolation into the future, with the flexibility of flows as a general purpose high-dimensional distribution model, while remaining computationally tractable. We show that it improves over the state-of-the-art for standard metrics on many real-world data sets with several thousand interacting time-series.

Published as a conference paper at ICLR 2021
MULTIVARIATE PROBABILISTIC TIME SERIES FORE- CASTING VIA CONDITIONED NORMALIZING FLOWS
Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann & Roland Vollgraf
Zalando Research
MühlenstraBe 25 10243 Berlin
Germany
Skashif.rasue symmaar.ccuster,,
roland. voland..olllgraf)eeaalaaaaa de
ABSTRACT
Time series forecasting is often fundamental to scientific and engineering prob- lems and enables decision making. With ever increasing data set sizes, trivial solution to scale up predictions to assume independence between interacting time series. However, modeling statistical dependencies can improve accuracy and enable analysis of interaction effeets. Deep lcarning methods are well suited for this problem, but multivariate models often assume simple parametric distribution and do not scale to high dimensions. In this work we model the multivariate temporal dynamics of time series via autoregressive deep learning model, where the data distribution is represented by conditioned normalizing flow. This com- bination retains the power of autoregressive models, such good performance in extrapolation into the future, with the flexibility of flows as general purpose high- dineesiieniio distribution model, while remaining computationally tractable. We show that it improves over the state-of-the-at -herrrr for standard metrics on many real-world data sets with several thousand interacting time-series.
INTRODUCTION
Classical time series forecasting methods such those in Hyndman & Athanasopoulos (2018) typically provide univariate forecasts and require hand-tuned features to model seasonality and other parameters. Time series models based recurrent neural networks (RNN), like LSTM (Hochreiter & Schmidhuber) 1997). have become popular methods due to their end-to-end training, the ease of incorporating exogenous covariates, and their automatic feature extraction abilities, which the hallmarks deep learning. Forecasting outputs can either be points or probability distributions, in which case the forecasts typically come with uncertainty bounds.
The problem of modeling uncertainties in time series forecasting is of vital importance for assessing how much trust the predictions for downstream tasks, such anomaly detection (business) decision making. Without probabilistic modeling, the importance of the forecast in regions of low noise (small variance around mean value) versus scenario with high noise cannot be distinguished. Hence, point estimation models ignore risk stemming from this noise, which would be of particular importance in some contexts such as making (business) decisions.
Finally, individual time series, in many cases, are statistically dependent on each other, and models need the capacity to adapt this in order to improve forecast accuracy (Tsay 2014). For example, to model the demand for a retail article, it is important to not only model its sales dependent on its own past sales, but also to take into account the effect of interacting articles, which can lead to cannibalization effects in the case of article competition. As another example, consider traffic flow in a network of streets as measured by occupancy sensors. A disruption on one particular street will also ripple to occupancy sensors of nearby streets univariate model would arguably not be able to account for these effects.
In this work, propose end-to-end trainable autoregressive deep learning architectures for proba- bilistic forecasting that explicitly models multivariate time series and their temporal dynamics by employing normalizing flow, like the Masked Autoregressive Flow (Papamakarios al. 2017)

---

Published as a conference paper at ICLR 2021
Real NVP (Dinh et al. 2017). These models are able to scale to thousands of interacting time series, we show that they are able to learn ground dependency structure on toy data and we establish new state-of-the-an results diverse real world data sets by comparing to competitive baselines. Additionally, these methods adapt to a broad class of underlying data distribution on account using normalizing flow and our Transformer based model is highly efficient due the parallel nature of attention layers while training.
The paper first provides some background context in Section 2 We cover related work in Section Section | introduces our model and the experiments are detailed in Section [ We conclude with some discussion Section 6 The Appendix contains details of the datasets, additional metries and exploratory plots of forecast intervals as well as details of our model.
2
BACKGROUND
2.1 DENSITY ESTIMATION VIA NORMALIZING FLOWS
Normalizing flows Tabak Turner 2013 Papamakarios et al. 2019) are mappings from ID to R. such that densities px the input space X' = IR" transformed into some simple distribution pz (e.g. an isoroopic Gaussian) on the space 2 RRD. These mappings, f r 2, are composed of a sequence of bijections invertible functions. Due the change of variables formula can express px (x) by
px(x)) = pz( det (21(x)
0x
where à 0f(x)/Ox is the Jacobian of f at x. Normalizing flows have the property that the inverse (z) is easy evaluate and computing the Jacobian determinant takes O(D) time.
The bijection introduced by Real NVP (Dinh et al. (017) called the coupling layer satisfies the above two properties. It leaves part of its imputs unchanged and transforms the other part via functions of the un nn-transformed variables (with superscript denoting the coordinate indices)
y Iid x [yd+l.D x+1I @ exp(s(x<lid) - t(x1:d).
where O is element wise product, s() is scaling and t() translation function from Rd D-d given by neural networks. To model a nonlinear density map f fxx), number of coupling layers which + VK-1 + Z are composed together all the while alternating the dimensions
map X' + 3/1 which are unchanged and transformed. Via the change variables formula the probability density function (PDF) of the flow given a data point be written
log (x) log pz (z) log (ōt(0x))x) log
Σ log dot((yy/0y (1)
Note that the Jacobian for the Real NVP is block-triangular matrix and thus the log-determinant of each map simply becomes
log det(ōy, = log exp(sum(s; ))1,
(2)
where sum() is the sum over all the vector elements. This model, parameterized by the weights of the scaling and translation neural networks is then trained via stochastic gradient descent (SGD) on training data points where for each batch D we maximize the average log likelihood (I) given by
1
L
1x;px(x;0).
D
xED
In practice, Batch Normalization (Ifff & Szegedy 2015) is applied as a bijection to outputs of successive coupling layers to stabilize the training of normalizing flows. This bijection implements the normalization procedure using a weighted moving average of the layer's mean and standard deviation values, which has to be adapted to either training inference regimes.
2

---

Published as a conference paper at ICLR 2021
The Real NVP approach can be generalized, resulting in Masked Autoregressive Flows (Papamakarios et a 2017) (MAF) where the transformation layer is built as an autoregressive neural network in the sense that it takes in some input E R2 and outputs y (y'..... . ?²) with the requirement that this transformation is invertible and any output y cannot depend on input with dimension indices 2 i, x2e. The Jacobian of this transformation is triangular and thus the Jacobian determinant is tractable. Instead of using RNN to share parameters across the D dimensions of x one avoids this sequential computation by using masking, giving the method its name. The inverse however, needed for generating samples, is sequential.
By realizing that the scaling and translation function approximators don't need to be invertible, is straight-forward to implement conditioning of the PDF px (x(h) on some additional information E " we concatenate h to the inputs of the scaling and translation function approximators of the coupling layers, i.e. s(concat(x!* 1:d h)) and ((concat(x" 1:d h)) which are modified to map Rd+H - R D-d Another approach is to add bias computed from h to every layer inside the and t networks as proposed by Korshunova et al. (2018). This does not change the log-determinant of the coupling layers given by (2). More importantly for us, for sequential data, indexed by t, we can share parameters across the different conditioners h, by using RNNs Attention in an autoregressive fashion.
For discrete data the distribution has differential entropy of negative infinity, which leads arbitrary high likelihood when training normalizing flow models, even on test data. To avoid this one can dequantize the data, often by adding Uniform(0 (0, 1) noise to integer- data. The log- likeliidl of the resulting continuous model is then lower-bounded by the log-likelihood of the discrete one as shown in Theis et al. 6
2.2 SELF-ATTENTION
The self-attention based Transformer (Vaswani et al., 2017) model has been used for sequence modeling with great success. The multi-head self-attention mechanism enables it to capture both long- and short-term dependencies in time series data. Essentially, the Transformer takes in sequence X = [1...... x E TxD, and the multi-head self-attention transforms this into H distinct query Q2 = XW
2 A key K, = XWW h. and value V XWX matrices, where the W w and W are learnable parameters. After these linear projections the scaled dot-product attention computes sequence of vector outputs via:
A
(Q Vdx
Oj = Attentionn Qi. KR. Vv. = sof tmax
M
V
where a mask M can be applied to filter out right-ward attention (or future information leakage) by setting its upper-triangular elements to -o and normalize by dæ the dimension of the W matrices. Afterwards, all H outputs O₁ are concatenated and linearly projected again.
One typically uses the Transformer in an encoder-decoder setup, where some warm-up time series is passed through the encoder and the decoder can be used to learn and autoregressively generate outputs.
RELATED WORK
Related this work models that combine normalizing flows with sequential modeling in some way. Transformation Autoregressive Networks (Oliva et al., 2018) which model the density of a multi-variate variable x R. as D conditional distributions P="'''=' ........"" ....,.").... where the conditioning is given by mixture model coming from the state of a RNN, and is then transformed via a bijection. The PixeISNAIL. (Chen et al., 2018) method also models the joint as a product of conditional distributions, optionally with some global conditioning, via causal convolutions and self-attention (Vaswani et al. 2017) to capture long-term temporal dependencies. These methods are well suited to modeling high dimensional data like images, however their use in modeling the temporal development of data has only recently been explored for example in VideoFlow (Kumar et al., 2019) in which they model the distribution of the next video frame via flow where the model outputs the parameters of the flow's, 's base distribution via a ConvNet, whereas our approach will be based on conditioning of the PDF as described above.
3

---

Published as a conference paper at ICLR 2021
Using RNNs for modeling either multivariate or temporal dynamics introduces sequential com- putational dependencies that are not amenable to parallelization. Despite this, RNNs have been shown to be very effective in modeling sequential dynamics. A recent work in this direction (Hwang et al. 2019) employs bipartite flows with RNNs for temporal conditioning to develop a conditional generative model of multivariate sequential data. The authors use a bidirectional training procedure to learn a generative model of observations that together with the temporal conditioning through a RNN, can also be conditioned on (observed) covariates that are modeled as additional conditioning variables in the latent space, which adds extra padding dimensions to the normalizing flow.
The other aspect of related works deals with multivariate probabilistic time series methods which are able to model high dimensional data. The Gaussian Copula Process method (Salinas et al. 2019a> a RNN- time series method with a Gaussian copula process output modeled using a low-rank covariance structure to reduce computational complexity and handle non-Gaussian marginal distributions. By using low-rank approximation of the covariance matrix they obtain a computationally tractable method and are able to scale to multivariate dimensions in the thousands with ofate-of-the-art results. We will compare our model to this method in what follows.
TEMPORAL CONDITIONED NORMALIZING FLOWS
We denote the entities of multivariate time series by x R for E (1...., D) where t is the time index. Thus the multivariate vector at time is given by Xx E RD, We will in what follows consider time series with & E [1, sampled from the complete time series history of our data, where for training we will split this time series by some context window to) and prediction window [to. T7.
In the DeepAR model (Salinas et al. 2019b). the log- likeliiddhood of each entity æi time step t E [8o, 2" is maximized given an individual time series' prediction window. This is done with respect to the parameters of the chosen distributional model (e.g. negative binomal for count data) via the state of a RNN derived from its previous time step -;- and current covariates c₁. 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 these distribution parameters, e.g. a softplus () for the scale parameter of the Gaussian.
A simple model for multivariate real-valued data could use a factorizing distribution in the emissions. Shared parameters can then learn patterns across the individual time series through the temporal component- but the model falls short of capturing dependencies in the emissions of the model. For this, full joint distribution at each time step must be modeled, for example by using multivariate Gaussian model. However, modeling the full covariance matrix not only increases the number of parameters of the neural network by O(D²), making learning difficult, but computing the loss becomes expensive when D is large. Furthermore, statistical dependencies in the emissions would be limited second-order effects. These models are referred to as Vec-LSTT in Salinas al. (2019a).
We wish to have a scalable model of D interacting time-series x,, and further to use a flexible distribution model on the emissions that allows for capturing and representing higher order moments. To this end, we model the conditional joint distribution at time t of all time series px (x,) |h;;) f) with a flow, e.g. Real NVP, conditioned on either the hidden state of RNN time or an embedding of the time series up to / 1 from an attention module. In the case of an autoregressive RNN (either LSTM or GRU (Chung et al., 2014)), its hidden state h, is updated given the previous time step observation x and the current time step's covariates C₁ (as in Figure
h
(oONNcccttxxx1,....0 h,-1)
This model is autoregressive since it consumes the observation of the last time step x,-1 as well as the recurrent state h to produce the state h, on which we condition the current observation.
To get a powerful and general emission distribution model, we stack K layers of a conditional flow module (Real NVP or MAF) and together with the RNN, we arrive our model of the conditional distribution of the future of all time series, given its past E [1, to) and all the covariates in E [1,7]. As the model is autoregressive can be written as a product of factors
( X1:to-
0) = II Pxx (xx bh; 8),
(4)
where 0 denotes the set of all parameters of both the flow and the RNN.
4

---

Published as a conference paper at ICLR 2021
For modeling the time evolution, we also investi- gate an encoder-decoder Transformer (Vaswami et al. 2017) architecture where the encoder em- beds X1 11tto= and the decoder outputs the con- ditioning for the flow over xt,T via masked attention module. See Figure for a schematic of the overall model in this case. While training, care has to be taken to prevent using information from future time points as well as to preserve the autoregressive property by utilizing mask that reflects the causal direction of the progressing time, i.e. to mask out future time points. The Transformer allows the model to access any part of the historic time series regardless of temporal distance (Li et al. 2019) and thus is potentially able to generate better conditioning for the nor- malizing flow head.
Batcch
Goupling Layer
Batchh Norm
Coupling Layer
x,
Norm
%
RNN
RNN Conditioned Real NVP model schematic at time 1, consisting of K blocks of coupling layers and Batch Normalization, where in each coupling layer we condition X and its trans- formations on the state of a shared RNN from the previous time step x and current time covari- ates C₁ which are typically time dependent and time independent features.
Figure 1:
In real-world data the magnitudes of different
time series can vary drastically. To normalize scales, we divide each individual time series by their training window means before feeding it into the model. At inference the distributions are then correspondingly transformed with the same mean values to match the original scale. This rescaling technique simplifies the problem for the model, which reflected in significantly improved empirical performance as noted in Salinas et al. 2019b).
4.1 TRAINING
Given D, a batch of time series, where for each time series and each time step have X, IR and their associated covariates C₁ we maximize the log-likelihood given by and (3, i.e.
Coupling Loyer
Coupling Layer
T
L
Σ
px (x:
Muls Head Alm
Massled
Muli Head
DIT
XLTED I=
Altin
via SGD using Adam (Kingma & Ba 2013) with respect to the parameters O of the conditional flow and the RNN or Transformer. In practice, the time series X1:T in a batch D are selected from a random time window of size within training data, and the relative time steps kept constant. This allows the model to learn to cold-start given only the covariates. This also increases the size of our training data when the training data has small time history and allows us to trade off computation time with memory consumption especially when D T are large. Note that information about absolute time is only available to the RNN or Transformer via the covariates and not the relative position of x, in the training data.
Sagel
Cooy
Xu7-1
Transformer Conditioned Real NVP model schematic consisting of an encoder-decoder stack where the encoder takes in some context length of time series and then uses it to generate conditioning for the prediction length portion of the time series via a causally masked decoder stack. Note that the positional encodings are part of the covariates and unlike the RNN model, here all x time points trained in parallel.
Figure 2:
The Transformer has computational complexity D) compared to a RNN which is O(TD²), where is the time series length and the assumption that the dimension of the hidden states are proportional to the number of simultaneous time-series modeled. This means for large multivariate time series, i.e. D > T. the Transformer flow model has smaller computational complexity and unlike the RNN, all computation while training, over the time dimension happens in parallel.
4.2 COVARIATES
We employ embeddings for categorical features (Charrington) 2018), which allows for relationships within a category, or its context, to be captured while training models. Combining these embeddings
5

---

Published as a conference paper at ICLR 2021
as features for time series forecasting yields powerful models like the first place winner of the Kaggle Taxi Trajectory Prediction challenge (De Brébisson et al. 2015). The covariates с 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 time periods wish to forecast.
4.3 INFERENCE
For inference we either obtain the hidden state h4s by passing a "warm up" time series X1:f1 through the RNN or use the cold-start hidden state, i.c. we set h2, = h1 = Ô, and then by sampling a noise vector Zt,- Rº from isotropic Gaussian, go backward through the flow obtain sample of our time series for the next time step, X₁2 (z,, [h,,). conditioned on this starting state. We then use this sample and its covariates to obtain the next conditioning state h via the RNN and repeat till our inference horizon. This process of sampling trajectories from some initial state can be repeated many times to obtain empirical quantiles of the uncertainty of our prediction for arbitrary long forecast horizons.
The attention model similarly uses a warm-up time series X1:4-11 and covariates and passes them through the encoder and then uses the decoder output the conditioning for sampling from the flow. This sample then used again in the decoder to iteratively sample the next conditioning state, similar to the inference procedure in seq-to-seq models.
Note that we do not sample from a peduuedd-temperatee model, e.g. by scaling the variance of the isotropic Gaussian, unlike what is done in likelihood-based generative models (Parmar al. 2018 to obtain higher quality samples.
EXPERIMENTS
Here we discuss a toy experiment for sanity-checking our model and evaluate probabilistic forecasting results six real-world data sets with competitive baselines. The source code of the model, as well as other time series models, is available at lattps://github: http://cc//naaanoeeeeeeccco pytorca-ts
5.1 SIMULATED FLOW IN A SYSTEM OF PIPES
In this toy experiment, check if the inductive bias of incorporating relations between time series is learnt in our model by simulating flow of a liquid in system of pipes with valves. See Figure 3 for depiction of the system.
Liquid flows from left to right, where pressure the first sensor in the system is given by So = X + 3, X Gamma(1, 0.2) in the shape/scale parameterization of the Gamma dis- tribution. The valves are given by Vi, V2 erid Beta(0.5, 0.5). and we have
V1
S1
SO
S3
S2
Figure System of pipes with liquid flowing from left to right with sensors (S;) and valves (V;).
V
S; = V +
So +
1/2
for i E (1, and finally S3 S1 + S2 E3 with e, (00,, 0.1). With this simulation we check whether our model captures correlations in space and time. The correlation between S₁ and S2 results from both having the same source, measured by Sg. This is reflected by Cov(S1, S2) > which is captured by our model as shown in Figure 4 left.
The cross- covaricve structure between consecutive time points in the ground truth and as captured by our trained model is depicted in Figure 4 right. It reflects the true flow of liquid in the system from So at time l S₁ and S2 at time &. on to S3 time 2.

---

Published as a conference paper at ICLR 2021
Ground truth
Fitted model
Ground truth
Fitted model
v1
V1
V1
measured at
V2
at
51
51
52
52
s3
v2 v1
50
s1
s2 23
v1 v2 s1 52 measured
Figure 4: Estimated (cross-ccovariance matrices. Darker means higher positive values. left: Co- variance matrix for a fixed time point capturing the correlation between S1 and S2. right: Cross- covariance matrix between consecutive time points capturing true flow of liquid in the pipe system.
Table 1: Test set CRPS comparison (lower is better) of models from Salinas al. (019a) and our models LSSM-RRel--- ReV.-NVP LSTM-MAF and Trans f Temer-MAF. The two best methods are in bold and the mean and standard errors of our methods are obtained by rerunning them 20 times.
Data set
Veg-LSTM
Vec-LSIM
SP
GP
LSTX
LSTX
Trans Eormer
ind-sca
lawrenk-Copula
Copula
Real- I-VF
MAF
MAP
Exchange
0.008
0.007=0.000
0.009=0.cmm
0.007
0.0064
0.005:00008
0.005
Solar
0.319_0.011
0.368_0.012
0.337a0.024
0.331
0.315.0.om
40.044
Electricity
0.025
0.064=0.00s
0.022=0.000
0.024:0.002
0.024
0.0208
0.0207
Traffic
0.087.000m
0. 0..0000.00
0.079_0.00m
0.078.2.002
0.078
0.069
0.05610000
Tax1
0.506
0.326
0.183=0.00
0.208
0.175500.01
0.161:00
0.179
Wikipedia
0.133=00002
0.211=0.033
1.483=1.034
0.056=20004
0.078
0.067=00001
0 +..63
5.2 REAL WORLD DATA SETS
For evaluation we compute the Continuous Ranked Probability Score (CRPS) (Matheson Winkler 1976) each individual time series, as well the sum of all time series (the latter denoted by CRPSsum). CRPS measures the compatibility of a cumulative distribution function with an observation x:
CRPS(F,2) (F(=) If(2) 2
(5)
where I(x < 2} is the indicator function which is one x: S 2 and zero otherwise. CRPS is a proper scoring function, hence CRPS attains its minimum when the predictive distribution F and the data distribution are equal. Employing the empirical CDF of F. i.e. F(2) Ei 1(X; =} with samples X; F' as a natural approximation of the predictive CDF, CRPS can be directly computed from simulated samples of the conditional distribution 4 each time point (Jordan al. 2019). We take 100 samples to estimate the empirical CDF in practice. Finally, CRPS is obtained by first summing across the D time-series both for the ground- data, and sampled data (yielding F (t) for each time point). The results are then averaged over the prediction horizon, i.e. formally
CRPS E, CRPS
Arum (t),
Σ
Our model trained on the training split each data set, and for testing we use a rolling windows prediction starting from the last point seen in the training data set and compare it the test set.
We train on Exchange (Lai et al. 2018), Solar (Lai et al. 2018 Electricio siciitiii Traffid Tax and wi i widi open data sets, preprocessed exactl in Salinas et al. 2019a). with their

---

Published as a conference paper at ICLR 2021
properties listed in Table of the appendix. Both Taxi and Wikipedia consist of count data and are thus dequantized before being fed to the flow (and -canedc
We compare method using LSTM and two different normalizing flows (LSTM-Real-VV)
and LSTM-MAF based on Real NVP and MAF, respectively) as well as a Transformer with the most competitive baseline probabilistic models from Salinas et al. (2019a) on the six data sets and report the results in Table [ Vec-LSTM-ind scal scalinggee outputs the parameters of an independent Gaussian distribution with mean-scaling, Vec via Copula process. GP-scaling unrolls LSTM with scaling on each individual time series before reconstructing the joint distribution via a low-rank Gaussian. Similarly, GP GP-Copula unrolls a LSTT on each individual time series and then the joint emission distribution given by a low-rank plus diagonal covariance Gaussian copula.
model with MAF (Trans
Cooplle a parametrizes low-rank plus diagonal covariance
Data
Sa maples
Absolute Difference
200
0.5
400
400
400
0.2
800
800
800
0.1
400
600
800
600
200
600
Figure 5: Visual analysis of the dependency structure extrapolation of the model. Left: Cross- covariance matrix computed from the test split Traffic benchmark. Middle: Cross-covariance matrix computed from the mean of 100 sample trajectories drawn from the Trans Transformer former-Mar model's 's extrapolation into the future (test split). Right: The absolute difference of the two matrices mostly shows small deviations between ground-truth and extrapolation.
In Table [ we observe that MAF with either RNN or self-attention mechanism for temporal condi- tioning achieves the state-of-the-an (to the best our knowledge) CRPSsum on all benchmarks. Moreover, bipartite flows with RNN either also outperform or are found to be competitive w.r.t. the previous state-of-tha -art---rttt--tttt- results as listed in the first four columns of Table Further analyses with other metrics (e.g. CRPS and MSE) are reported in Section B of the appendix.
To showcase how well our model captures dependencies in extrapolating the time series into the future versus real data, plot in Figure cross-covariance matrix of observations (plotted left) as well as the mean of 100 sample trajectories (middle plot) drawn from Trans formerr-r-rrrr --MAF model for the test split of Traffic fc data set. As can be seen, most of the covariance structure especially in the top-left region of highly correlated sensors is very well reflected in the samples drawn from the model.
6 CONCLUSION
We have presented general method to model high-dimensional probabilistic multivariate time series by combining conditional normalizing flows with an autoregressive model, such as a recurrent neural network or an attention module. Autoregressive models have long-standing reputation for working very well for time series forecasting, as they show good performance in extrapolation into the future. The flow model, on the other hand, does not assume any simple fixed distribution class, but instead can adapt to broad range of high-dimensional data distributions. The combination hence combines the extrapolation power of the autoregressive model class with the density estimation flexibility of flows. Furthermore, it is computationally efficient, without the need resorting to approximations (e.g. low-rank approximations of covariance structure in Gaussian copula methods) and is robust compared to Deep Kernel learning methods especially for large D. Analysis on six commonly used time series benchmarks establishes the state-of-the performance against competitive methods.
natural way to improve our method is to incorporate a better underlying flow model. For example, Table shows that swapping the Real NVP flow with MAF improved the performance, which is
8

---

Published as a conference paper at ICLR 2021
consequence of Real NVP lacking in density modeling performance compared MAF Likewise, we would expect other design choices of the flow model to improve performance, e.g. changes to the dequantization method, the specific affine coupling layer more expressive conditioning, say via another Transformer. Recent improvements to flows, e.g. as proposed in the Flow++ (Ho et al. 2019), to obtain expressive bipartite flow models, models to handle discrete categorical data (Tran et al. 2019), are left as future work to assess their usefulness. To our knowledge, it is however still an open problem how to model discrete ordinal data via flows-which would best capture the nature of some data sets (e.g. sales data).
ACKNOWLEDGMENTS
K.R.: would like to thank Rob Hyndman and Zaeem Burq for the helpful discussions and suggestions. I would like to acknowledge the traditional owners of the land on which I have lived and worked, the Wurundjeri people of the Kulin nation who have been custodians of their land for thousands of years. pay my respects to their elders, past and present as well as past and present aboriginal elders of other communities.
We wish to acknowledge and thank the authors and contributors of the following open source libraries that were used in this work: GluonTS (Alexandrov et al.] (020), NumPy (Harris et al., (020). Pandas (Pandas development team 2020), Matplotlib (Hunter) 2007) and PyTocc (Paszke et al. 2019). We would also like to thank and acknowledge the hard work of the reviewers whose comments and suggestions have without doubt help improve this paper.
REFERENCES
Alexander Alexandrov. Konstantinos Benidis, Michael Bohlke-Schneiderr Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Danielle C. Maddix, Syama Rangapuram, David Salinas, Jasper Schulz, Lorenzo Stella, Ali Caner Türkmen, and Yuyang Wang. Gluon'TS: Probabilistic and Neural Time Series Modeling in Python. Journal of Machine Learning Research, 21(116):1-6, 2020. URL http::/5mlx org wwwpar/g/pppprreeeeeee www.rr//vv111/1 820. html
Sam Charrington. TWiML Al Podcast: Systems and software for machine learning scale with LmGg
Jeff Dean, 2018. URL
XI Chen, Nikhil Mishra, Mostafa Rohaninejad, and Pieter Abbeel. PixeISNAIL: An improved autoregressive generative model. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 872, Stockholmsmässan, Stockholm Sweden, 2018. PMLR. URL http::/prroceednnge mlr. ttees/vv0/cheeee
html
Junyoung Chung, Caglar Gulcehre, Kyunghyun Cho, and Yoshua Bengio. Empirical evaluation of gated recurrent neural networks on sequence modeling. In NIPS 2014 Workshop on Deep Learning, December 2014, 2014.
Djork-Aréé Clevert Thomas Unterthiner, and Sepp Hochreiter. Fast and accurate deep network learning by exponential linear units (elus). In Yoshua Bengio and Yann LeCun (eds.), 4th Inter- national Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico May 2-4, 2016, Conference Track Proceedings, 2016. URL
/tt.ab//1111... 07289
Emmanuel de Bézenac, Syama Sundar Rangapuram. Konstantinos Benidis, Michael Bohlke- Schneider, Richard Kurle, Lorenzo Stella, Hilaf Hasson, Patrick Gallinari. and Tim Januschowski. Normalizing Kalman Filters for Multivariate Time series Analysis. In Advances in Neural Infor- mation Processing Systems, volume 33, Curran Associates, Inc., 2020.
Alexandre De Brébisson, Étienne Simon, Alex Auvolat, Pascal Vincent, and Yoshua Bengio. Artificial neural networks applied to taxi destination prediction. In Proceedings of the 201 International Conference ECML PKDD Discovery Challenge- Volume 1526, ECMLPKDDDC 15, pp. 40-51. Aachen, Germany, Germany, 2015. BEUR-WS.org htLhttp::// //dt: acm. wwtp///ggtatioo.o 01210056112333366177 3056178

---

Published as a conference paper at ICLR 2021
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP. In International Conference on Learning Representations 2017 (Conference Track), 2017. URL
Charles R. Harris, K. Jarrod Millman, St' Sf'n J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk. Matthew Brett, Allan Haldane, Jaime Fern't andez del R R'1o, Mark Wiebe, Pearu Peterson, Pierre G' Braard-Marrchnt, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 3010:um268887772555353----2 September 2020. doi: 10. 10.10388/441588-00---666
URL /attps://bii.org//////
020 2649-2
Jonathan Ho, Xi Chen, Aravind Srinivas, Yan Duan, and Pieter Abbeel. Flow++: Improving flow- based generative models with variational dequantization and architecture design. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 2722-2730, Long Beach, California, USA, 2019. PMLR. URL http: sttp://proceeddings mlr .tt.sress..pr./....me.me. press/ v97/hol 1ahol9..hai. html
S. Hochreiter and J. Schmidhuber. Long Short-Term Memory. Neural Computation, 98)):1735- November 1997. ISSN 0899-76677 doi: 10.1 162/neco. 1997.9.8. 19715.55..
J. D. Hunter. Matplotlib: A 2d graphics environment. Computing in Science & Engineering, 9()): 90-95, 2007. doi: 10. 109/
Scong Jae Hwang, Zirui Tao, Won Hwa Kim, and Vikas Singh. Conditional recurrent flow: Con- ditional generation of longitudinal samples with applications to neuroimaging. In The IEEE International Conference on Computer Vision (ICCV), October 2019.
R.J. Hyndman and G. Athanasopoulos. Forecasting: Principles and pructice. OTexts, 2018. ISBN 9780987507112.
Rob Hyndman, Anne Koehler, Keith Ord, and Ralph Snyder. Forecasting with exponential smoothing. The state space approach, chapter 17, pp. 287-300. Springer- Springer-Verlag, 2008. doi: 10.1007/
978-3-540-71918-22
Sergey loffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings of the 32Nd International Conference on International Conference on Machine Learning- Volume 37, ICML' pp. 448 456. JMLR.org. ctt.stiioonccm....mm 3000=3105111... 3045167
2015. http: acm.
Alexander Jordan, Fabian Krüger, and Sebastian Lerch. Evaluating probabilistic forecasts with scoringRules. Journal of Statistical Software, Articles, 90(12) 2019. ISSN 1548-7660. doi: 10.18637/js..0090.12. v09...... URL htt /ww::://wwwwww.... http::///wwwwsstssssststteee wwtp::/ssattttt.oroooooooo
Diederick P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), 2015.
Iryna Korshunova, Yarin Gal, Arthur Gretton, and Joni Dambre. Conditional BRUNO: A Deep Recurrent Process for Exchangeable Labelled Data. In Bayesian Deep Learning workshop, NIPS, 2018.
Rahul G Krishnan, Uri Shalit, and David Sontag. Structured inference networks for nonlinear state space models. In AAAI, 2017.
Manoj Kumar, Mohammad Babaeizadeh, Dumitru Erhan, Chelsea Finn. Sergey Levine, Laurent Dinh, and Durk Kingma. VideoFlow: A Flow-Based Generative Model for Video. In Workshop on Invertible Neural Nets and Normalizing Flows, JCMI., 2019.
Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu. Modeling long- and short-term temporal patterns with deep neural networks. In The 41st International ACM SIGIR Conference Research & Development in Information Retrieval, SIGIR 18, pp. 95-104, New York, NY, USA. 2018. ACM. ISBN 978- 4503--5500--557222 doi: 10.1145/3209978.22100 320000 URL http:
7/doi acm. www.oo.g/10000... 14/ 321000
10

---

Curran Published as a conference paper at ICLR 2021
Shiyang Li, Xiaoyong Jin, Yao Xuan, Xiyou Zhou, Wenhu Chen, Yu-Xiang Wang, and Xifeng Yan. Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 5254. Curran Associates, Inc., 2019.
H. Lütkepohl. New Introduction to Multiple Time Series Analysis. Springer Berlin Heidelberg, 2007. ISBN 9783540262398. URL /tt.ttt////ooos..ooolllee google de/ midokksddddmmodeee EC
James E. Matheson and Robert I.. Winkler. Scoring rules for continuous probability distributions. Management Science, 22(10):1087-1096, 1976.
Junier Oliva, Avinava Dubey, Manzil Zaheer, Barnabas Poczos, Ruslan Salakhutdinov, Eric Xing, and Jeff Schneider. Transformation autoregressive networks. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 3898-3907, Stockholmsmässann, Stockholm Sweden, 201 8. PMLR. http:://ooccccinn....rrsrs/////p
htm
Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. N-BEATS: Neural basis ex- pansion analysis for interpretable time series forecasting. In International Conference on Learning Representations, 2020. URL htRItttts///:ppereeeet http://pperreeeeeeeeeeeerrorrrra tttt://forrmmmdddreeeeeip/m
The Pandas development team. pandas-dev/pandas: Pandas, February 2020. URL https://dot orc/10..521/zenodb.. 3509134
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. Advances in Neural Information Processing Systems 30, 2017.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference, 2019.
Niki Parmar, Ashish Vaswani, Jakob Uszkoreit, Lukasz Kaiser, Noam Shazeer, Alexander Ku, and Dustin Tran. Image transformer. Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 4055-4064, Stockholmsmässaan Stockholm Sweden, 10-15 Jul 2018, PMLR. URL http:://proceedings. mlr. wwwps /vB0 pa rmar 8a. .html
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy. Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Py PyTorch: An imperative style, perforrefree deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché Buc. E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. Associates, Inc., 2019.
David Salinas, Michael
Laurent Callot, Roberto Medico, and Jan Gasthaus. High- dimensional multivariate forecasting with low-rank Gaussian copula processes. In H. Wallach, H. Larochelle, A. Beygelzimer, F. Alché Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 6824 6834. Curran Associates, Inc., 2019a.
David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. DeepAR: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 2019b. ISSN 0169-2070. URL http: //wp::://wwww.... S01 5016920701930188
com/ wttp://wwtt...ooooeeeceeeceeeeee wttp//ccieecce//rrrcce/e///////e
E.G. Tabak and C.V. Turner. A family of nonparametric density estimation algorithms. Communica- fions on Pure and Applied Mathematics, 145-164,, 2013.
L. Theis, A. van den Oord, and M. Bethge. A note the evaluation of generative models. In International Conference on Learning Representations, 2016. URL
abs/ has/15111 01844 1111v:1511.018..
11

---

Published as a conference paper at ICLR 2021
Dustin Tran, Keyon Vafa, Kumar Agrawal, Laurent Dinh, and Ben Poole. Discrete flows: Invertible generative models of discrete data. In H. Wallach, H. Larochelle, A. Beygelzimer, d' Alché Buc. E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 14692-14701. Curran Associates, Inc., 2019.
Ruey S. Tsay. Multivariate Time Series Analysis: With R and Financial Applications. Wiley Series in Probability and Statisties. Wiley, 2014. ISBN 97811 18617908.
Roy van der Weide. GO-GARCH: a multivariate generalized orthogonal GARCH model. Journal of Applied Econometrics, 17(5):549-564, 2002. doi: 10. 1002/jae.688.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, E ukasz Kaiser, and Illia Polosukhin. Attention is all you need. In I. Guyon, U.V. Luxburg, S. Bengio, H. Wallach, R. Fergus, Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems 30, pp. 5998 6008. Curran Associates, Inc., 2017. URL http: 7/papers ppapers.niss.oom www.piisscccpppppppppptpp 7781 att-ttteeettttt inn--issalll -oo---allll---------- need. pdf
A
DATA SET DETAILS
Table Properties of the data sets used in experiments.
DIMENSION
TOTAL TIME STEPS
PREDICTION
DATA SET
DOMAIN
FREQ.
D
LENGTH
EXCHANCE
8
R
DAILY
6,071
30
SOLAR
R+
HOURLY
7,009
ELECTRICITY
370
R+
HOURLY
5,833
24
TRAFFIC
963
(0,
HOURLY
4,001
24
1,214
N
30-MIN
1,488
24
WIXTPROTA
2,000
N
DAILY
792
30
B ADDITIONAL METRICS
We used exactly the same open source code to evaluate our metrics as provided by the authors of Salinas al. (2019a).
B. COMPARISON AGAINST CLASSICAL BASELINES
We report test set CRPS results VAR ( Lütkepohl 2007) mutlivariate linear vector auto- regressive model with lags corresponding to the periodicity of the data, VAR-Lasso a Lasso regularized VAR, GARCH (van der Weide 2002) multivariate conditional heteroskedastic model, GP Gaussian process model, KVAE (Krishnan et all 2017) variational autoencoder on top of linear state space model and VES a innovation state space model (Hyndman al. (0088) Table Note that VAR-Lasso, KVAE and VES metrics are from (de Bézenac et al.) 2020).
Table 3: Test set CRPS (lower is better) of classical methods and our Trans former-MA. model, where the mean and standard errors of our model are obtained over a mean of 20 runs.
Transformas
Data set
VAR
VAR Lasao
GP
GARCH
VES
KVBE
HAF
Exchange
0.010=0.30
0.012+0.000
0.011+0.001
0.020#0.com
0.005=a.00
0.01-22.002
0. 0.005=00.00
Solar
0.524:n.00.
0.51:0.008
0.828+0.0.
0.869:0.c0
0.949.0c
0.34=50c25
0.301=9.014
Electricity
0.031=0.00
0.025e0.00
0.941=00.0.
0.278s0.co
0.88 oca
0.051+0.019
0. 0.007=00009
Traffic
0.144:n.30
2.198
0.368
0.35
0.140.005
Taxi
0.292:0..0
0.425
0.179.00.02
Nikipadia padia
3.420.005
0.9
0.095.0.012
0.0633.9.0
12

---

Published as a conference paper at ICLR 2021
B.2 CONTINUOUS RANKED PROBABILITY SCORE (CRPS)
The average marginal CRPS over dimensions D and over the predicted time steps compared to the test interval is given in Table 4
Table 4: Test set CRPS comparison (lower is better) of models from Salinas et al. 2019a and our models L.STM-Ree Rea NVVP, LSTM-MAF and Trans ormer-MAF. The mean and standard errors obtained by re-running each method three times.
Veo 1SIM
LSIM
LSTM
LSTM
Tranaforme
Data set
Ind-scal1-
lowrank-Copola
scal inn
Copula
Real NNVP
MAP
MAP
Exchange
0.013:0.000
0.009
0.017.0.000
0.008 o0e
0.010.20001
0.01220.005
0.01210.00
Solar
0.434:00
0.384=0.010
0.415=0.00m
0.371 20.022
0.365
0.378+0.032
0.368
Electricity
0.039=0.0m 001
0.081=0.0m%
0.053=0.000
0.056 0.002
0.039=3.001
0.051=0.000
0.052+0.000
Traftic
0.168
0.165 oot
0.140
0.133 oo1
0.172
P..124...
0.134.1.001
Taxi
0,586 004
0.116=a0.00
0.316=00.44
0.360 201
0.327=3.001
0.314=0.00m ona
0.377
Mikipedia
0.379:0.00
0.247=0.001
1.519=11017
0.236i=0.000
0.333+00.00
0,282+0.002
0.274:00007
B.3 MEAN SQUARED ERROR (MSE)
The MSE is defined as the mean squared error over all the time series dimensions D and over the whole prediction range with respect to the test data. Table shows the MSE results for the the marginal MSE.
Table 5: Test set MSE comparison (lower is better) of models from Salinas et al. 2019a and our models LSTM Real-elllleeeeeeell -VV,, 4STM-MAF and Trans for Former -aeemmee------
Data set
Vea-1.STM Ind acaling
Vec-1.8TM lowrank Copula
GP acaling
T.STM
1.5TM MAF
Trans former
Copula
Exchange
1.6 10
1.9 x
2.
1.7 10
2.4 10
3.8 10
3.4 10
Solar
9.3 10²
2.9 x 10
1.1 10
9.8 10"
9.1 x 102
9.8 102
9.3 16"
Electricity
2.1
5.5 10"
1.8 x 10a
24x
1.8
2.0
Tractic
6.3 10-1
1.5 10-3
5.2 10
6.9 10
6.9 10-1
x 10-4
5.0 10-1
Taxi
7.3 10"
5.1 10
2.7 x 10a
3.1xx00 10"
2.6 10a
2.4 10
4.5 10*
Rikipedia
7.2
3.8 107
5.5 x 107
4.0
3.8
3.1 x
C UNIVARIATE AND POINT FORECASTS
Univariate methods typically give better forecasts than multivariate ones, which is counter-intuitive. the reason being the difficulty in estimating the cross-series correlations. The additional variance that multivariate methods add often ends up harming the forecast, even when one knows that individual time series are related. Thus as additional sanity check, that this method is good enough to improve the forecast and not make worse, we report the metrics with respect to a modern univariate point forecasting method as well as a multivariate point forecasting method for the Tra f f data set.
Figure 6 reports the metrics from LSTNet (Lai et al. 2018 multivariate point forecasting method and Figure 7 reports the metrics from I--BEATS (Oreshkin et al. 2020) univariate model. As can be seen, our methods improve on the metrics for the Traffi i data set this pattern holds for other data sets in experiments, As visual comparison, have also plotted the prediction
intervals using our models Figures g 9 and
D EXPERIMENT DETAILS
D. FEATURES
For hourly data sets used hour of day, day of week, day of month features which are normalized. For daily data sets we use the day of week features. For data sets with minute granularity we use minute of hour, hour of day and day of week features. The normalized features concatenated to
13

---

Published as a conference paper at ICLR 2021
0:2
0:20
4.00
0.08
022
2.02
0.64
0.64
1527
304
1527
0.22
01.5
0.23
0.0
0.22
01000
0.000
505V
12-6
5050
1375
18:30
05.00
3040
1009
0800
2003
3008
2003
0.37
6.8
0.23
0.8
0.06
0.80
0.06
0:0
0.0
0.00
5000 1270 1050 5100
10.70
01
00.20
23.00
0000
12
50:00
2108
00.00
52.00 5000 1278
2050
Figure 6: Point forecast and test set ground-truth from LSTNet multivariate model for Traffic data of the first 16 of 963 time series. CRPS 0.125, CRPS: 0.202 and MSF: 7.4 10-4
the RNN or Transformer input at each time step. also concatenate lag values as inputs according to the data set's 's time frequency; [1, 24, 168] for hourly data, [1, 14] for daily and [1, 4, 12, 24, 48] the half-hourly data.
D.2 HYPERPARAMETERS
We batch sizes of 64, with 100 batches per epoch and train for maximum of 40 epochs with learning rate of le-3. The LSTM hyperparameters were the ones from Salinas et al, (2019a) and we used K stacks of normalizing flow bijections layers. The components of the normalizing flows (f and g) are linear feed forward layers (with fixed input and final output sizes because we model bijections) with hidden dimensions of 100 and ELU (Clevert al. 2016) activation functions. We sample 100 times to report the metrics the test set. The Transformer uses H = 8 heads and 3 encoding and m decoding layers and dropout rate of 0. All experiments run on a single Nvidia V-100 GPU and the code to reproduce the results will be made available after the review process.
14

---

Published as a conference paper at ICLR 2021
4.300
0.2
0.20
8:00
0.66
0.0
0.9
600
4.500
0.2
2.00
-007
0.35
0..9
0100-
0.35
0.22
0.47
0.0
21.0
0800
1308
1800
1400 1900
1808
(076
(000
1900
18300
0600
0608
4000
2008
1308
1829
0.24
0.28
0:2-
0.37
0.1
2007
18/8
02.50
1858
4.25
0.44
0.04
0:00
0.25
13.0
6360
1358 1287
1800
1800 1400
R000
1500
18:00
54:0
13:00
0000
0600 13.30
00
5400
1568
0618
00000 1358 1829 RISH 0000
Figure 7: Point forecast and test set ground-truth from N-BEATS univariate model for Traffic data of the first 16 of 963 time series. CRPS 0.174, CRPS: 0.228 and MSE: 8.4 10-4
15

---

Published as a conference paper at ICLR 2021
$.44-
5:40
000
3.00
680
-007
2.00
150
0.14
0.35
0.35
0.42
03.00
1308
1800
9400 1800
1808
(076
9820
13:00
1869
1900 13.30
18300
0600
1129
2008
1308
1829
0.35
0.28
0.25
0.28
0.80
0.00
1207
1838
02.50
0.3
0:0
000
33.0 6000 1800
1358 1800 R090 1400 1500
18:00
1400 13:00
0900 0600 13.30
1400
3568
0618
13.0 0050
1358 1869 RISH 0600
1287
Figure 8: Prediction intervals and test set ground-truth from LSTM -STT-MEAL NVP model for Traffic data of the first 16 of 963 time series.
16

---

Published as a conference paper at ICLR 2021
0:28
pradition
0:24
A04
0.6
000
D00P
12.00
52:00
0.21
0.35
1258
(460
1068
1258
12:50
1.57
0.0
50.00
1206
2000
10.00
3008
120
0.18
0.30
4.25
0.04
0:2
0.04
0.5
0.25
3.50
2000
0600 1308 180
1400 1000
5440 23:00
1800 13:00
1440
3008
0638
550 20:00
2000
1308
1827
RE0 1548
Figure 9: Prediction intervals and test set ground-trun from Traffic c data of the first 963 time series.
model for
17

---

Published as a conference paper at ICLR 2021
308
1256 BOU
SATA 20Ud
30 1ESS 19T8 2oud
1246 (60 18:00 5000
12:0
DOOD
205
354
20.0 1208
00:00 1308 1800
1300 1800
3700
3450
1150
1208
1923
1992
3:00
Figure 10: Prediction intervals and test set ground-truth from LSTM-REAL-NVZ model for Electricity data of the first 16 of 370 time series.
18

---

Published as a conference paper at ICLR 2021
1208
Wii LNve 1865.
2318 OWe La 1892.
9136 LO0D
12:6
12700
1248
3000
2002
2016
200
10.00
02900
02.00
Figure 11: Prediction intervals and set ground-truth from Transformer MAF model for Electricity data of the first 16 of 370 time series.
19


Reviews:
[{"review_id": "o6_WuIukKHK", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes a method to provide probabilistic forecasts of multivariate time-series taking dependencies between series into account even for large dimensions. The approach consists in using a normalizing flow to model the distribution of observations at a time-step condition on a state that can be obtained either with a RNN or a transformer. The motivation of using the normalizing flow is to be able to model various type of distributions without having to make specific hypothesis on the distribution which could hinder accuracy. Experiments are performed both on a synthetic task and on real-world datasets where accuracy is shown to outperform previous methods.\n\n### Strengths:\n+ empirical improvements obtained over real-world datasets when measuring CRPS-sum accuracy metric for a problem that has significant applications\n+ the clarity of the model description: models are well described and introduced. Generally speaking the paper is well written and reads well.\n\n### Weaknesses:\n- the improvement cannot be traced directly to the author contribution (the use of flow to model the multivariate distribution) given that they change the underlying model of the baseline (a LSTM) to a GRU. Given that the difference are small, it is unclear using GRU is not the source of improvement. \n- the stability regarding the flow architecture is not analyzed whereas it is its core contribution. No detail is given on how the flow architecture was found, two architectures are used with 4 or 5 layers and no detail is given how they are chosen (the methods compared use only a single architecture on all datasets).\n- the paper claims that one contribution is that their method is efficient but this is not backed-up by either theoretical complexity nor empirical data. \n\t- The time complexity given for RNN (O(TD^2)) is misleading given that this is a naive approach, for instance the GP copula methods have O(TD) runtime with a RNN. \n\t- The transformer method may still be faster due to parallelism but to make an efficiency claim, the runtime should be measured given that the theoretical complexity is higher\n\nDespite the accuracy improvement and the fact that the paper is well written and motivated, I am currently not inclined to accept it given some ambiguities raised by the experiments but I believe those points could be clarified in the author rebuttal.\n\nIn particular, I think the following points needs clarification:\n- running their method with an LSTM would be very important to have an apple-to-apple comparison and discard that the improvement comes from GRU rather than normalizing flows\n- How did you find K = 4 or K = 5 parameters? by looking at a single dataset? How do you choose between 4 and 5 in your test-set? Given that the methods you compare with use a single configuration, the comparison should also be made with a single configuration.\n- a performance claim is made but is not backed-up at the moment. The complexity of the models are O(TD^2) or O(T^2D) but the GP copula baseline is O(TD). Rather than complexity, perhaps runtime would be a better advocate for this point? (otherwise I believe the efficiency claim should be dropped as it is not backed up by theoretical complexity nor runtime measurement)\n\n# Additional feedback (no part of the decision)\n\nDetails:\n* section 4: \"In the case of an autoregressive\".\n* the artificial experiments is ill-suited to show-case the benefit of the method, couldnt one model directly a full-covariate matrix in this case given the low-dimensionality? How many time-steps were considered for the sequence of pipes?\n\nFinally, I would like to note that I reviewed a previous version of this paper and that I updated my review with the changes made by the author.", "extra_scores": {}, "extra_text": {}}, {"review_id": "aEYYRrR2A2t", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "# Update after author discussion\nMy original review mentioned an apparent similarity to VideoFlow and a weak experimental setup as the main reasons for my score. Through discussions with the authors, I am now convinced that VideoFlow is indeed a quite different model, which significantly adds to the novelty of the proposed model. \n\nAs for the experiments, I had misunderstood the setup. In terms of probabilistic multivariate methods, my experience lies in the multi-output Gaussian process domain where it is common to consider tasks where only some outputs are considered missing (i.e., part of the test set), thus testing the models' ability to make use of partial information at a given input. This was the setup I had in mind, but the proposed model is not easily adapted for this. This is not a limitation of the proposed model as much as an inherent difficulty of handling missing data with neural networks. The experimental setup used in this paper instead considers a training set and a test set entirely separated in time, which is not a setup where GPs usually perform well, so the additional experiments I had requested do not make much sense. Indeed, the baseline methods chosen for comparison in the paper are both reasonable and strong.\n\nIn conclusion, the authors have satisfactorily addressed my concerns with the submission. I still think the proposed idea is good (although not groundbreakingly novel), the paper is well-written, and the model thoroughly discussed and tested. I view the submission as a significant contribution to time-series modelling and is therefore recommending acceptance.\n\n\n----------------------------------------------------------\n\n\n# Summary\nThe paper proposes to use normalising flows for probabilistic modelling of multivariate time-series. By taking advantage of normalising flows' ability to model complex interactions between a large number of outputs with arbitrary output distributions, the proposed model naturally scales to many concurrent time-series. The temporal component is modelled with an autoregressive model, such as an RNN or a Transformer, which outputs a latent representation of previous times and covariates, which is then used as the condition in the flow transformations.\nCombined, the two parts of the model allow for probabilistic trajectories to be sampled for predictions into the future. The model is tested on six datasets and compared to four other models, showing very good performance.\n\n# Evaluation\n## Strong points\n* Using normalising flows for multivariate time-series modelling is a good idea.\n* Two variations of the proposed model are tested, showing how design decision change the performance.\n\n## Weak points\n* Limited novelty over previous models.\n* Experimental section a bit weak.\n\n## Recommendation\nBorderline reject.\n\nThe model in itself is interesting, as is the application to time-series modelling. However, the model does seem very related to VideoFlow, which is mentioned briefly in the paper but never discussed properly, so the novelty seems limited. Furthermore, the experimental section only considers competing models from one other paper.\n\n## Detailed feedback\nThe paper is well-written, well-structured, and the background section is comprehensive and detailed. The writing is clear and of high quality.\n\nThe related work section discusses some relevant work, but given that the paper concerns probabilistic time-series modelling, I would have expected to see some discussions of multi-output Gaussian processes (MO-GPs) and state-space models. I am also quite surprised how little discussion there is of Kumar et al. (2019) given that VideoFlow is very related to the proposed model. The way the two models introduce temporal dependencies differ, but they are otherwise similar. A comparison of the two models seems in order.\n\nThe proposed model is certainly interesting, and the idea of using flows to correlate time-series is very good. Again, it is not entirely new, but it is a very underexplored path. There is not much novelty from a theoretical perspective either, as the model is mostly a combination of two sub-models; a normalising flow architecture and an autoregressive component. I see this as a strength of a kind, as it makes the model quite modular, which should make it interesting to practitioners - it's just that there is a trade-off in novelty.\n\nTo me, the main shortcoming of the paper is experimental evaluation. The model is evaluated on six datasets, which are quite different, so this is definitely good. The experiments are run several times and results are reported with standard deviations - excellent! The model, however, is compared against two baselines and two competing models, which seems thorough, but the two baselines are very related to each other, as are the two competing models. All four models are taken from the same paper (Salinas et al. 2019a). Given how similar this model is to VideoFlow, I do not understand why that model was not tested. Furthermore, I think it would have made sense to include independent Gaussian processes (GPs) as baselines and to test additional multi-output GP models, e.g., a semi-parametric latent factor model and the Gaussian process autoregressive regression model (GPAR, Requeima et al., 2019).\nIn addition, I find the CRPS quite hard to interpret, and I think it would have made sense to include more metrics than just those based on CRPS. MAE is included in the supplementary, but reporting something like the standardised MSE and the (negative) log probability of the data under the model would have been useful too. In particular, the log probability would have been interesting to see.\n\nLastly, I think it would be interesting to see plots of the predictive densities by the different models for a couple of datasets. This is, of course, only a qualitative assessment, but I am very interested in seeing how the probabilistic trajectory of the model looks like in terms of smoothness and predictive density.\n\n# Questions\nComments that will affect my score:\n1. Could you please include experiments with other models? In particular, I'd be interested in seeing the performance of VideoFlow and GPAR.\n2. Could you please include additional metrics? In particular, I'd be interested in seeing the negative predictive log probability.\n3. Could you please include some plots of the predictions of the models on future data so we can assess the quality of the predictions? (Including these in the supplementary would be fine.)\n\nMinor questions not likely to affect my score:\n1. You say that \"Shared parameters can then learn patterns across the individual time series through the temporal component - but the model falls short of capturing dependencies in the emissions of the model.\" Why can't shared parameters capture dependencies between outputs?\n2. In section 4.3 you mention a cold start, but I am confused about when you would ever prefer this to a warm-start. Can you elaborate?\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "ZcbCWZMHWi6", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 9, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper presents a new approach to multivariate probabilistic time series forecasting. The authors propose to combine recurrent neural networks with conditioned normalizing flows to model the output distribution. They explain a few variations of the method and evaluate them against the baselines. \n\n*Quality*\nThe choice of the output distribution is always a hard step in the probabilistic time series forecasting using RNNs. The authors propose an interesting workaround by trying to learn the shape of the output distribution directly with conditioned normalizing flows. In combination with a Transformer model, this also has the potential to reduce training/inference time for high-dimensional datasets. The authors explain and motivate the proposed method really well and I think the paper is a useful contribution to the field of probabilistic time series forecasting.\n\n*Clarity*\nThe paper is very well-written. \n\n*Originality*\nTo my best knowledge, the proposed approach is new.\n\n*Significance*\nRather significant. The multivariate forecasting models suffer from the problems with scaling and the paper’s proposed combination of transformers with conditioned normalizing flows has the potential to overcome this problem.\n\nPros\n* Conditioned normalizing flows can ease the design choice of the output distribution\n* The method scales well in terms of number of time series to forecast\n* Great empirical results\n\nCons\n* Since the scaling is a big issue for multivariate models, it would be useful to compare training/inference times between the evaluated methods.", "extra_scores": {}, "extra_text": {}}, {"review_id": "1Bf7v1Ix5tI", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This work explores combining an RNN and a neural density estimator for forecasting in multivariate time series. RNN is stacked with a density estimator, MAF for best results, to forecast density of a multivariate time series at future time steps. In addition, variations of the architecture with attention and other density estimators are examined. The architecture, RNN+MAF and variants, is evaluated by CRPS score on several datasets.\n\nThe work may have  a value but I have many issues reading it. I have a feeling that the motivation and the  related work are incomplete. Too many variants of the architectures are discussed and unnecessary details are provided while important questions are not discussed.  Reconstructing the 'best' architecture from the paper is not a trivial task. Empirical evaluation is somewhat sketchy.\n\nIn detail, \n\n* while deep learning is a powerful tool for forecasting in multivariate timeseries, including settings with uncertainty, it is not the only way. 'Classical' approaches like Kalman filter and multi-task gaussian processes work well in many settings. \n\n* there was work on forecasting time series with uncertainty using end-to-end deep learning architectures (rather than involving GP GCP etc) which the paper does not cite. For example, https://deepai.org/publication/sequential-neural-models-with-stochastic-layers and other variants of stochastic RNNs, as well as neural processes, https://arxiv.org/abs/1807.01622. \n\n* As I read the paper, RNN+MAF gives the best results empirically, and MAF is indeed a sato density estimator, so this is quite expected. However, if the purpose is *forecasting* that is, obtaining the future states, I am wondering why isn't IAF used instead, which is more suited for sampling.\n\n* Empirical evaluation --- on the toy example, I am wondering why the learned cross-variance matrices are visually different in Figure 4. I would expect them to be visually indistinguishable on such a simple model. Even using simpler multi-variate time series models. In Table 1, \"the mean and standard errors are obtained by re-running ... three times\". I see that the same method of \"obtaining the mean and standard error\" appears in the cited work, but it does not make the method more credible. You cannot \"obtain\" or even reliably estimate the mean and the standard error from three runs. Either a statistically credible number of runs is needed, or there isn't much sense in reporting the standard error.\n\n~~~~~\n\nBased on the discussion and the updated manuscript, I am happy to say that most my concerns were addressed and at least partially resolved.  I have updated my score accordingly.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an approach to probabilistic time series forecasting based on combining autoregressive deep learning models with normalizing flows. In terms of strengths, time series forecasting is a fundamental problem. The proposed approach is a reasonable combination of existing model components that provides a flexible, end-to-end trainable framework for multivariate probabilistic forecasting. The experiments are well-conducted and the results outperform recently published methods. While the reviewers raised a number of questions, all of the reviewers agree that their questions have be answered satisfactorily by the authors during the discussion and the paper should be accepted. The authors should be sure to incorporate the reviewer suggestions and author responses into the final paper. 
Decision: Accept (Spotlight)