# DynaConF: Dynamic Forecasting of Non-Stationary Time-Series

Siqi Liu
Borealis AI &Andreas Lehrmann
Borealis AI

###### Abstract

Deep learning models have shown impressive results in a variety of time series
forecasting tasks, where modeling the conditional distribution of the future
given the past is the essence. However, when this conditional distribution is
non-stationary, it poses challenges for these models to learn consistently and
to predict accurately. In this work, we propose a new method to model non-
stationary conditional distributions over time by clearly decoupling
stationary conditional distribution modeling from non-stationary dynamics
modeling. Our method is based on a Bayesian dynamic model that can adapt to
conditional distribution changes and a deep conditional distribution model
that can handle large multivariate time series using a factorized output
space. Our experimental results on synthetic and popular public datasets show
that our model can adapt to non-stationary time series better than state-of-
the-art deep learning solutions.

##  1 Introduction

Time series forecasting is a cornerstone of modern machine learning and has
applications in a broad range of domains, such as operational processes [1],
energy [2], and transportation [3]. In recent years, models based on deep
neural networks have shown particularly impressive results [4, 5, 3, 1] and
demonstrated the effectiveness of deep feature and latent state
representations.

Despite this exciting progress, current time series forecasting methods often
make the implicit assumption that training and test data follow the same
distribution. In real-world applications this assumption is often violated,
which is known as _non-stationarity_ and poses serious practical challenges to
a model’s robustness and predictive power. The statistics literature contains
several related concepts of (non-)stationarity for time series, with weak and
strong stationarity being the most widely used ones [6, 7]. Common to these
variants of (non-)stationarity is that they are defined in terms of a
stochastic process’ joint or marginal distribution. For example, given a
discrete time series $\\{y_{t}\in\mathbb{R}\\}_{t\in\mathbb{Z}}$ and any
subset of time points $\\{t_{1},t_{2},\ldots,t_{k}\\}$, we call the time
series _strongly stationary_ if
$\forall\tau\in\mathbb{Z}:\mathrm{p}(y_{t_{1}},y_{t_{2}},\ldots,y_{t_{k}})=\mathrm{p}(y_{t_{1}+\tau},y_{t_{2}+\tau},\ldots,y_{t_{k}+\tau})$.

While non-stationarity in a stochastic process’ joint or marginal distribution
is important and has been widely studied [7, 8, 9, 10], we argue that temporal
_forecasting_ relies more heavily on the properties of a time-series’
_conditional_ distribution $\mathrm{p}(y_{t}|y_{t-B:t-1},\bm{x}_{t-B:t})$,
where $y_{t-B:t-1}=(y_{t-B},y_{2},\ldots,y_{t-1})$, $B\in\mathbb{Z}_{>0}$ can
be arbitrarily large, and $\bm{x}_{t}\in R^{Q}$ contains auxiliary
information. Most forecasting methods, from traditional approaches (e.g.,
Autoregressive Integrated Moving Average (ARIMA; [11]), Generalized
Autoregressive Conditional Heteroskedasticity (GARCH; [12, 13]), state-space
models (SSMs; [14])) to more recent models (e.g., recurrent neural networks
(RNNs; [15, 1]), temporal convolutional networks (TCNs; [16]), Transformers
[17, 4], Neural Processes (NPs; [18, 19])), rely on this conditional
distribution for predictions, but many of them implicitly assume its
stationarity by using time-invariant model parameters. A model with a
stationary conditional distribution can still handle non-stationarity in the
joint or marginal distribution, such as seasonality and trend, by conditioning
on extra features in $\bm{x}_{t}$, such as day of the week, but the
conditional distribution itself may also change due to (1) unobserved causes
and (2) new causes. For example, the daily number of posts from each user on a
social media platform is unlikely to be robustly predictable from historical
data, even with input features like day of the week, because (1) user activity
is often affected by events that are not reflected in the observable data
(e.g., illness); and (2) events that have not been seen before, such as a new
functionality being added to the platform, may occur and change the functional
pattern of the input-output relation between the conditioning and target
variables in an unpredictable way. How to deal with these changes in the
model’s _conditional_ distribution, which are based on a dynamic cause-effect
structure, is the main focus of this work.

Autoregressive (AR) models, TCNs, and Transformers model
$\mathrm{p}(y_{t}|y_{t-B:t-1},\bm{x}_{t-B:t};\bm{\psi})$ with time-invariant
parameters $\bm{\psi}$ and therefore assume stationarity in
$\mathrm{p}(y_{t}|y_{t-B:t-1},\bm{x}_{t-B:t})$. In contrast, SSMs, RNNs, and
NPs model $\mathrm{p}(y_{t}|y_{1:t-1},\bm{x}_{1:t};\bm{\psi})$, which has a
growing number of conditioning variables (note the time range $1\\!:\\!t$) and
therefore can potentially model different conditional distributions
$\mathrm{p}(y_{t}|y_{t-B:t-1},\bm{x}_{t-B:t})$ at different time points $t$.
However, these models need to achieve two goals using the same time-invariant
structure: (1) modeling $\mathrm{p}(y_{t}|y_{t-B:t-1},\bm{x}_{t-B:t})$; and
(2) modeling its changes over time. Because they do not incorporate explicit
inductive biases for changes in
$\mathrm{p}(y_{t}|y_{t-B:t-1},\bm{x}_{t-B:t})$, they either cannot learn
different (seemingly inconsistent) relations between the conditioning and
target variables (if the model capacity is limited) or tend to memorize the
training data and are not able to generalize to new changes at test time.

In this work we take a different approach to dealing with non-stationary
conditional distributions in time series. The core of our model, called
DynaConF, is a clean decoupling of the time-variant (non-stationary) part and
the time-invariant (stationary) part of the distribution. The time-invariant
part models a stationary conditional distribution, given some control
variables, while the time-variant part focuses on modeling the changes in this
conditional distribution over time through those control variables. Using this
separation, we build a flexible time-invariant conditional model and make
efficient inferences about how the model changes over time. At test time, our
model takes both the uncertainty of the conditional distribution and non-
stationarity into account when making predictions and can adapt to new changes
in the conditional distribution over time in an online manner.

##  2 Related Work

Time-series forecasting models have a rich history [7, 11], with many recent
advances due to deep neural networks [20]. Here we discuss relations of
related works to non-stationarity and our work.

Non-Stationary Marginal Distribution. There are three common ways of dealing
with non-stationarity in the marginal distribution, such as seasonality and
trend: (1) Data transformation. In ARIMA, taking the difference of the time
series over time can remove trend and seasonality. More advanced approaches
based on exponential smoothing [21] or seasonal-trend decomposition [22] have
also been combined with deep neural networks and achieve good performance [23,
24]. More recently, Kim et al. [10] propose to use reversible
normalization/denormalization on the input/output of the time series model to
account for (marginal) distribution shifts over time. (2) Inductive bias. As
an alternative to data transformations, the underlying ideas of exponential
smoothing and decomposition can also be incorporated into deep learning models
as inductive biases to deal with seasonality/trend end-to-end [25, 26, 27, 2].
Similarly, for models based on Gaussian processes, the inductive biases can be
added as specific (e.g. periodic/linear) kernel choices [28]. (3) Conditioning
information. Adding features such as relative time (e.g, day of the week) and
absolute time to the model input as conditional information is commonly used
in deep probabilistic forecasting models [1, 3, 5, 29, 4]. In our work, we
focus on proper handling of changes in the _conditional_ distribution. To deal
with marginal distribution shifts, we simply add conditional information as in
approach (3), although we could potentially utilize approach (1) and (2) as
well.

Non-Stationary Conditional Distribution. State-space models [30, 31, 32, 5,
33] and recurrent neural networks [1, 3, 4, 29] are among the most popular
choices to model temporal dependencies in time series. When these models are
applied to, and therefore conditioned on, the entire history of the time
series, they can theoretically “memorize” different conditional distributions
at different time points. However, for these models to generalize and adapt to
new changes in the future, it is critical to have appropriate inductive biases
built into the model. A common approach is to allow the state space model to
switch between a discrete set of dynamics, which can be learned from training
data [34, 35]. However, the model cannot adapt to continuous changes or
generalize to new changes that have not been observed in the training data. In
contrast, our model has explicit inductive biases to account for both
continuous and discontinuous changes, and to adapt to new changes.

Observation Model. The expressivity and flexibility of the observation model
is a topic that is especially relevant in case of multivariate time series.
Different observational models have been employed in time series models,
including low-rank covariance structures [3], auto-encoders [30, 36],
normalizing flows [31, 4], determinantal point processes [37], denoising
diffusion models [29], and probabilistic circuits [38]. In this work, to deal
with multivariate time series, we take a simple approach by assuming
conditional independence in the output dimensions for scalability and consider
more expressive observation models as orthogonal to our work.

Online Approach. Continual learning [39, 40, 41, 42, 43, 44] also addresses
(conditional) distribution changes in an online manner, but usually in a
multi-task supervised learning setting and not in time series. In this work,
our focus is on conditional distribution changes in time series, and the
conditional distribution can change either continuously or discontinuously
over time.

##  3 Method

We study the problem of modeling and forecasting time series with changes in
the conditional distribution
$\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$ over time $t$, where
$\bm{y}_{t}\in\mathbb{R}^{P}$ is the target time series, and
$\bm{x}_{t}\in\mathbb{R}^{Q}$ is the input containing contextual information,
under the following assumption:

###### Assumption 1

$\bm{y}_{t}$ only depends on a bounded history of $\bm{y}$ and $\bm{x}$ for
all $t$. That is, there exists $B\in\mathbb{Z}_{>0}$ such that for all $t$,
$\mathrm{p}(\bm{y}_{t}|\bm{y}_{<t},\bm{x}_{\leq
t})=\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$.

Although we assume $\bm{y}_{t}$ can only depend on the history up to $B$ time
steps, its conditional distribution can change over time based on information
beyond $B$ steps. Notice that this is not particularly restrictive in
practice, since we usually have a finite amount of training data while $B$ can
be arbitrarily large, although usually not needed. Specifically, given the
historical data $\\{(\bm{x}_{t},\bm{y}_{t})\\}_{t=1}^{T}$, the task is to fit
a model to the data and use it to forecast
$\\{\bm{y}_{t}\\}_{t=T+1}^{T+H},\\{\bm{y}_{t}\\}_{t=T+H+1}^{T+H+H},\ldots,$
continually with a horizon and step size $H$. For notational simplicity, we
assume that $\bm{x}_{t}$ only contains information known in advance at time
$t$, such as day of the week.

###  3.1 Decoupled Model

We assume the distribution
$\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$ to be parameterized
by $\bm{\theta}_{t}\in\mathbb{R}^{M}$:

$\bm{\theta}_{t}=f(\bm{y}_{t-B:t-1},\bm{x}_{t-B:t};\bm{\phi}_{t},\bm{\psi}).$
(1)

$f:\mathbb{R}^{B\times P}\times\mathbb{R}^{(B+1)\times Q}\to\mathbb{R}^{M}$
models the conditional distribution, with its own static parameters, denoted
collectively as $\bm{\psi}$, and is modulated by a dynamic control variable
$\bm{\phi}_{t}\in\mathbb{R}^{F}$, which can change over time according to a
dynamic process we define later. A property we try to guarantee in our model
is that if $\bm{\phi}_{t}$ stays the same at different time points, then the
conditional distribution
$\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$ stays the same as
well. Figure 1 shows an overview of our model.

Our key idea is to separate the time-variant part of the model from the time-
invariant part, instead of allowing all components to vary across time. This
simplifies the probabilistic inference and improves the generalization
capabilities of the learned model by allowing the time-invariant part to learn
from time-invariant input-output relations with time-variant modulations.

###  3.2 Conditional Distribution at One Time Point

First we describe how we model the conditional distribution
$\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$ at each time point
$t$ without accounting for non-stationary effects (Figure 1, bottom). We use a
neural network $g$ to summarize the historical and contextual information into
a vector $\bm{h}_{t}\in\mathbb{R}^{D}$ as
$\bm{h}_{t}=g(\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$. For example, $g$ could be a
multi-layer perceptron (MLP) or a recurrent neural network (RNN). The
parameters of $g$ are time-invariant and included in $\bm{\psi}$ (Eq. 1).

We note that a key distinction between our model’s use of an RNN and a typical
deep time series model using an RNN is that the latter keeps unrolling the RNN
over time to model the dynamics of the time series. In contrast, we unroll the
RNN for $B+1$ steps to summarize $(\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$ in the
exact same way at each time point $t$, i.e., we apply it in a _time-invariant_
manner.

We construct the distribution of $\bm{y}_{t}$ such that each dimension $i$ of
$\bm{y}_{t}$, denoted as $y_{t,i}$, is conditionally independent from the
others given $\bm{h}_{t}$, as this helps the learning and inference algorithms
to scale better with the dimensionality of $\bm{y}_{t}$. First we transform
$\bm{h}_{t}$ into $P$ vectors of dimension $E$,
$\bm{z}_{t,i}\in\mathbb{R}^{E}$, $E<D$, as

$\bm{z}_{t,i}=\tanh(\bm{W}_{z,i}\bm{h}_{t}+\bm{b}_{z,i}),\quad\forall
i=1,\ldots,P$ (2)

where $\bm{W}_{z,i}\in\mathbb{R}^{E\times D}$ and
$\bm{b}_{z,i}\in\mathbb{R}^{E}$, so $\bm{z}_{t,i}$ corresponds to $y_{t,i}$.
Then, from $\bm{z}_{t,i}$, we obtain the distribution parameters
$\bm{\theta}_{t,i}$ of $y_{t,i}$. Specifically, we assume a normal
distribution with a diagonal covariance for $\bm{y}_{t}$, so
$y_{t,i}\sim\mathcal{N}(\mu_{t,i},\sigma^{2}_{t,i})$ and
$\bm{\theta}_{t,i}:=(\mu_{t,i},\sigma^{2}_{t,i})$, with

$\mu_{t,i}=\bm{w}_{\mu,i}^{T}\bm{z}_{t,i}+b_{\mu,i},\quad\sigma_{t,i}=s(\bm{w}_{\sigma,i}^{T}\bm{z}_{t,i}+b_{\sigma,i}),$
(3)

where $\bm{w}_{\mu,i}\in\mathbb{R}^{E}$, $\bm{w}_{\sigma,i}\in\mathbb{R}^{E}$,
and $s:\mathbb{R}\to\mathbb{R}_{>0}$ is the soft-plus function. We use
$\bm{W}_{\mu}\in\mathbb{R}^{P\times E}$ and $\bm{b}_{\mu}\in\mathbb{R}^{P}$ to
denote the result of stacking $\bm{w}_{\mu,i}^{T}$ and $b_{\mu,i}$ along $i$,
and similarly $\bm{W}_{\sigma},\bm{b}_{\sigma}$, $\bm{z}_{t}$, $\bm{\mu}_{t}$,
$\bm{\sigma}_{t}$. As we will see in Section 3.5, the use of a Normal
distribution enables more efficient inference.

###  3.3 Conditional Distributions Across Time Points

We have explained how we model the conditional distribution
$\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$ at each time $t$. To
model changes in the conditional distribution over time, we first specify
which parameters to include in the control variable
$\bm{\phi}_{t}\in\mathbb{R}^{F}$, which changes over time and modulates the
conditional distribution (Figure 1, top).

We choose $\bm{\phi}_{t}:=\operatorname{vec}(\bm{W}_{\mu})$, i.e.,
$\bm{\phi}_{t}$ is the vectorization of $\bm{W}_{\mu}$. Recall that
$\bm{W}_{\mu}$ transforms $\bm{z}_{t}$ into the mean $\bm{\mu}_{t}$ of
$\bm{y}_{t}$, where $\bm{z}_{t}$ is essentially a summary of the information
in the conditioning variables $(\bm{y}_{t-B:t-1},\bm{x}_{t-B:t})$. By allowing
$\bm{W}_{\mu}$ to be different at each time point $t$, the conditional mean of
$\bm{y}_{t}$, $\mathrm{E}[\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t}]$, can
change as well. We could allow $\bm{W}_{\sigma}$ to change over time to model
a time-variant conditional variance as well, but focusing on $\bm{W}_{\mu}$
reduces the dimensionality of $\bm{\phi}_{t}$ and enables more efficient
inference utilizing Rao-Blackwellization (see Section 3.5).

We propose to decompose $\bm{\phi}_{t}$ into a dynamic stochastic process
$\bm{\chi}_{t}\in\mathbb{R}^{F}$ and a static vector
$\bm{b}_{\phi}\in\mathbb{R}^{F}$ as

$\bm{\phi}_{t}=\bm{\chi}_{t}+\bm{b}_{\phi}.$ (4)

The intuition is that $\bm{b}_{\phi}$ captures the global information of
$\bm{\phi}_{t}$ and acts as a baseline, while $\bm{\chi}_{t}$ captures the
time-variant changes of $\bm{\phi}_{t}$ relative to $\bm{b}_{\phi}$.

We assume that $\bm{\chi}_{t}$ follows the following generative process

$\displaystyle\pi_{t}\sim\mathcal{B}(\lambda);$
$\displaystyle\bm{\chi}_{t}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{0}),\text{ if
}\pi_{t}=0;$ (5)
$\displaystyle\bm{\epsilon}_{t}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{d});$
$\displaystyle\bm{\chi}_{t}=\bm{\chi}_{t-1}+\bm{\epsilon}_{t},\text{ if
}\pi_{t}=1.$

$\mathcal{B}$ and $\mathcal{N}$ denote the Bernoulli and normal distributions,
respectively.

$\pi_{t}\in\\{0,1\\}$ is a binary random variable that either generates the
current $\bm{\chi}_{t}$ as a new sample drawn from a global distribution
$\mathcal{N}(\bm{0},\bm{\Sigma}_{0})$, or as a continuation from the previous
$\bm{\chi}_{t-1}$ following a simple stochastic process in the form of a
random walk. The intention is to allow $\bm{\chi}_{t}$ to change both
continuously (when $\pi_{t}=1$) through the random walk, and discontinuously
(when $\pi_{t}=0$) through the global distribution, which captures the variety
of possible changes of $\bm{\chi}_{t}$ in its parameter $\bm{\Sigma}_{0}$.

We assume $\bm{\chi}_{t}$ to start at $t=B$, since it controls the
_conditional_ distribution, whose first observation occurs at $t=B+1$. For the
initial $\bm{\chi}_{B}$, we assume generation from
$\mathcal{N}(\bm{0},\bm{\Sigma}_{0})$ as well. Our intention is that
$\mathcal{N}(\bm{0},\bm{\Sigma}_{0})$ should be the distribution to generate
new $\bm{\chi}_{t}$ whenever there is a drastic change in the conditional
distribution of $\bm{y}_{t}$, so at $t=B$ it is natural to use that
distribution.

Recall that $\bm{\chi}_{t}\in\mathbb{R}^{F}$, with $F=P\times E$. We propose
to separate $\bm{\chi}_{t}$ along the $P$ dimensions of $\bm{y}_{t}$ into $P$
groups. For each $i=1,\ldots,P$, we define $\bm{\chi}_{t,i}\in\mathbb{R}^{E}$
as in Eq. 5. The final $\bm{\chi}_{t}$ is the concatenation of
$\bm{\chi}_{t,i}$ for all $i$. The intuition is to allow the group of
components of $\bm{\chi}_{t}$ modulating each dimension of $\bm{y}_{t}$ to
change independently from the others, corresponding to the conditional
independence assumption we made in Section 3.2, so we can sample a subset of
dimensions of $\bm{y}$ in each iteration during training to scale to high-
dimensional time series.

Figure 1: Overview. DynaConF is built on the principle of a clean decoupling
of stationary conditional distribution modeling (red) and non-stationary
dynamics modeling (blue). We predict the parameters ($\bm{\theta}$) of the
conditional distribution (green) by aggregating time-invariant local context
($\bm{z}$) and modulating this context with time-variant global dynamics
($\bm{\phi}$) driven by a random walk ($\bm{\chi}$) with Bernoulli-restarts
($\bm{\pi}$).

###  3.4 Learning

All parameters of the conditional distribution model and the prior are learned
by fitting the model to the historical training data
$\\{(\bm{y}_{t},\bm{x}_{t})\\}_{t=1}^{T}$. We train our model by maximizing
the marginal log-likelihood, where the latent variables $\bm{\chi}_{t}$ are
marginalized out. Given a trajectory of $\bm{\chi}_{B:T}$, the conditional
log-likelihood is

$\log\mathrm{p}(\bm{y}_{B+1:T}|\bm{y}_{1:B},\bm{x}_{1:T},\bm{\chi}_{B:T})=\sum_{t=B+1}^{T}\log\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t},\bm{\chi}_{t}).$
(6)

Marginalizing out $\bm{\chi}_{t}$ gives us the log-likelihood objective

$\log\mathrm{p}(\bm{y}_{B+1:T}|\bm{y}_{1:B},\bm{x}_{1:T})=\log\int\mathrm{p}(\bm{y}_{B+1:T}|\bm{y}_{1:B},\bm{x}_{1:T},\bm{\chi}_{B:T})\mathrm{p}(\bm{\chi}_{B:T})\mathrm{d}\bm{\chi}_{B:T}.$
(7)

Since the integral is intractable, we instead introduce a variational
distribution $\mathrm{q}(\bm{\chi}_{B:T})$ and maximize the following
variational lower-bound $\mathcal{L}$ of the log-likelihood in Eq. 7:

$\mathcal{L}:=\mathrm{E}_{q}[\log\mathrm{p}(\bm{y}_{B+1:T}|\bm{y}_{1:B},\bm{x}_{1:T},\bm{\chi}_{B:T})]+\mathrm{E}_{q}\left[\log\frac{p(\bm{\chi}_{B:T})}{\mathrm{q}(\bm{\chi}_{B:T})}\right]\leq\log\mathrm{p}(\bm{y}_{B+1:T}|\bm{y}_{1:B},\bm{x}_{1:T}).$
(8)

Based on the conditional independence structure of the prior process, we
construct the variational distribution $\mathrm{q}(\bm{\chi}_{B:T})$ similarly
as an autoregressive process, but we assume a simple Normal distribution at
each time step for efficient sampling and back-propagation. At the beginning,
we assume $\mathrm{q}(\bm{\chi}_{B})=\mathrm{p}(\bm{\chi}_{B})$. Then,
conditioning on the previous $\bm{\chi}_{t-1}$, we recursively define

$\mathrm{q}(\bm{\chi}_{t}|\bm{\chi}_{t-1})=\mathcal{N}(\bm{a}_{t}\odot\bm{\chi}_{t-1}+(1-\bm{a}_{t})\odot\bm{m}_{t},\mathrm{diag}(\bm{s}_{t}^{2})),\quad\forall
t=1,\ldots,T,$ (9)

where $\bm{a}_{t}$, $\bm{m}_{t},\bm{s}_{t}\in\mathbb{R}^{F}$ are variational
parameters. Intuitively, $\bm{a}_{t}$ is a gate that chooses between
continuing from the previous $\bm{\chi}_{t-1}$, with noise
$\mathcal{N}(0,\mathrm{diag}(\bm{s}^{2}_{t})$, and using a new distribution
$\mathcal{N}(\bm{m}_{t},\mathrm{diag}(\bm{s}_{t}^{2}))$.

We note that both terms in Eq. 8 factorize over time $t=B+1,\ldots,T$ as
follows:

$\displaystyle\mathcal{L}$
$\displaystyle=\mathrm{E}_{q(\bm{\chi}_{B:T})}\left[\sum_{t=B+1}^{T}\log\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t},\bm{\chi}_{t})\right]+\mathrm{E}_{q(\bm{\chi}_{B:T})}\left[\sum_{t=B+1}^{T}\log\frac{p(\bm{\chi}_{t}|\bm{\chi}_{t-1})}{\mathrm{q}(\bm{\chi}_{t}|\bm{\chi}_{t-1})}\right]$
(10)
$\displaystyle=\sum_{t=B+1}^{T}\mathrm{E}_{q(\bm{\chi}_{t})}\left[\log\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t},\bm{\chi}_{t})\right]+\sum_{t=B+1}^{T}\mathrm{E}_{q(\bm{\chi}_{t-1:t})}\left[\log\frac{p(\bm{\chi}_{t}|\bm{\chi}_{t-1})}{\mathrm{q}(\bm{\chi}_{t}|\bm{\chi}_{t-1})}\right].$

In the above derivation, $\mathrm{p}(\bm{\chi}_{B})=\mathrm{q}(\bm{\chi}_{B})$
were canceled out due to the definition of $\mathrm{q}(\bm{\chi}_{B})$. The
expectations in this equation can be evaluated by Monte-Carlo sampling from
$\mathrm{q}(\bm{\chi})$ with the reparameterization trick [45] for back-
propagation.

In cases where training efficiency is critical, sequential sampling in the
autoregressive posterior might not be feasible, and we leverage Inverse
Autoregressive Flows (IAFs; [46]) for parallel sampling. While the latter does
not reflect the structure of the prior and true posterior as closely as the
autoregressive posterior, it is a viable option in cases where parallel
sampling is necessary. In the IAF-based posterior,
$\mathrm{q}(\bm{\chi}_{B:T})$ can be sampled over $t=B:T$ in parallel by
sampling from a standard Normal distribution of dimension $T-B+1$ and
transforming the sample through several MADE layers [47]. In our case, each
MADE layer not only takes the output from the previous layer as input but is
also conditioned on a learnable embedding representing each dimension of
$\bm{\chi}_{t}$. Details of this posterior model can be found in Appendix A.

We develop optimization procedures based on stochastic gradient descent (SGD)
that work in practice on large datasets utilizing our modeling assumptions
from Section 3.2 and Section 3.3. Specifically, we alternate between
optimizing the conditional distribution model and the prior and posterior of
the dynamic model with different sampling strategies to accommodate high
dimensionalities and long time spans. More details of our optimization
procedures are in Appendix B.

###  3.5 Forecasting

At test time we are given the past observations $\bm{y}_{1:T}$ as well as the
input features $\bm{x}_{1:T+H}$, including $H$ future steps, and infer the
conditional distribution
$\mathrm{p}(\bm{y}_{T+1:T+H}|\bm{y}_{1:T},\bm{x}_{1:T+H})$. Based on our
modeling assumptions, the latter can be computed as

$\displaystyle\mathrm{p}(\bm{y}_{T+1:T+H}|\bm{y}_{1:T},\bm{x}_{1:T+H})$ (11)
$\displaystyle=$
$\displaystyle\int\mathrm{p}(\bm{y}_{T+1:T+H}|\bm{y}_{T+1-B:T},\bm{x}_{T+1-B:T+H},\bm{\chi}_{T+1:T+H})\mathrm{p}(\bm{\chi}_{T+1:T+H}|\bm{y}_{1:T},\bm{x}_{1:T})\mathrm{d}\bm{\chi}_{T+1:T+H}.$

The first factor in the integrand can be computed recursively by step-by-step
predictions based on our conditional distribution model given
$\bm{\chi}_{T+1:T+H}$,

$\mathrm{p}(\bm{y}_{T+1:T+H}|\bm{y}_{T+1-B:T},\bm{x}_{T+1-B:T+H},\bm{\chi}_{T+1:T+H})=\prod_{t=T+1}^{T+H}\mathrm{p}(\bm{y}_{t}|\bm{y}_{t-B:t-1},\bm{x}_{t-B:t},\bm{\chi}_{t}).$
(12)

The second factor in the integrand can be further factorized into

$\mathrm{p}(\bm{\chi}_{T+1:T+H}|\bm{y}_{1:T},\bm{x}_{1:T})=\int\mathrm{p}(\bm{\chi}_{T+1:T+H}|\bm{\chi}_{T})\mathrm{p}(\bm{\chi}_{T}|\bm{y}_{1:T},\bm{x}_{1:T})\mathrm{d}\bm{\chi}_{T}.$
(13)

We use Rao-Blackwellized particle filters [48] to infer
$\mathrm{p}(\bm{\chi}_{T}|\bm{y}_{1:T},\bm{x}_{1:T})$, so our model can keep
adapting to new changes in an online manner. Specifically, we jointly infer
$\bm{\pi}_{t}$ and $\bm{\chi}_{t}$ with particles representing $\bm{\pi}_{t}$
and closed-form inference of $\bm{\chi}_{t}$. Once we have the posterior
samples of $\mathrm{p}(\bm{\chi}_{T}|\bm{y}_{1:T},\bm{x}_{1:T})$, we use the
prior model to sample trajectories of $\bm{\chi}_{T+1:T+H}$ conditioned on the
samples of $\bm{\chi}_{T}$. With the sample trajectories of
$\mathrm{p}(\bm{\chi}_{T+1:T+H}|\bm{y}_{1:T},\bm{x}_{1:T})$, we sample the
trajectories of
$\mathrm{p}(\bm{y}_{T+1:T+H}|\bm{y}_{T+1-B:T},\bm{x}_{T+1-B:T+H},\bm{\chi}_{T+1:T+H})$
using the aforementioned step-by-step predictions with our conditional
distribution model.

##  4 Experiments

Table 1: Baselines Models.

Method | S | R | MV | Implementation
---|---|---|---|---
DeepAR [1] | ✓ |  |  | GluonTS[49, 50]
DeepSSM [5] | ✓ |  |  | GluonTS
TransformerMAF [4] | ✓ | ✓ | ✓ | PyTorchTS[51]
DeepVAR (I)nd. [3] | ✓ | ✓ | ✓ | GluonTS
DeepVAR (C)opula [3] |  | ✓ | ✓ | GluonTS
GP-Scaling [3] |  | ✓ | ✓ | GluonTS
GP-Copula [3] |  | ✓ | ✓ | GluonTS
LSTM-NVP [4] |  | ✓ | ✓ | PyTorchTS
LSTM-MAF [4] |  | ✓ | ✓ | PyTorchTS
TimeGrad [29] |  | ✓ | ✓ | PyTorchTS

[S = Synthetic; R = Real-World; MV = Multivariate]

We compare our approach with $2$ univariate and $8$ multivariate time series
models, both on synthetic (Section 4.1) and real-world (Section 4.2) datasets;
see Table 1 for an overview, including references to the relevant literature
and implementations. We note that DeepVAR is also called Vec-LSTM in previous
works [3, 4]. For our own model we consider two variants: an ablated model
without the dynamic updates to the conditional distribution described in
Section 3.3 (StatiConF); and our full model including those contributions
(DynaConF). In both cases we experiment with different encoder architectures.
For synthetic data, we use either a two-layer MLP with $32$ hidden units (* –
MLP) or a point-wise linear + $\tanh$ mapping (* – PP) as the encoder. For
real-world data, we use an LSTM as the encoder.

###  4.1 Experiments on Synthetic Data

#### Datasets

For our experiments on synthetic data we simulate four conditionally non-
stationary stochastic processes for $T=2500$ time steps. We use the first
$1000$ steps as training data, the next $500$ steps as validation data, and
the remaining $1000$ steps as test data:

  * •

(AR(1) – Flip) We simulate an AR(1) process,
$y_{t}=w_{t}y_{t-1}+\epsilon,\epsilon\sim\mathcal{N}(0,1)$, but resample its
coefficient $w_{t}$ from a uniform categorical distribution over
$\\{-0.5,+0.5\\}$ every 100 steps to introduce non-stationarity.

  * •

(AR(1) – Dynamic) We simulate the same process as above but now resample
$w_{t}$ from a continuous uniform distribution $\mathcal{U}(-1,1)$ every 100
steps.

  * •

(AR(1) – Sin) We simulate the same process as above but now resample $w_{t}$
according to $w_{t}=\sin(2\pi t/T)$. Different from the two processes above,
this process has a continuously changing non-stationary conditional
distribution with its own time-dependent dynamics.

  * •

(VAR(1) – Dynamic) This process can be viewed as a multivariate generalization
of AR(1) – Dynamic and is used in our comparisons with multivariate baselines.
We use a four-dimensional VAR process with an identity noise matrix. Similar
to the univariate case, we resample the entries of the coefficient matrix from
a continuous uniform distribution $\mathcal{U}(-0.8,0.8)$ every 250 steps. In
addition, we ensure the stability of the resulting process by computing the
eigenvalues of the coefficient matrix and discard it if its largest absolute
eigenvalue is greater than $1$.

#### Experimental Setup

For univariate data (AR(1) – Flip/Sin/Dynamic), we compare our approach with
the two univariate baselines DeepAR and DeepSSM. For multivariate data (VAR(1)
– Dynamic), most baselines are redundant because their focus is on better
observation models, while the underlying temporal backbone is similar. Since
our synthetic observation distributions are simple, we compare with two model
families that differ in their temporal backbone: DeepVAR (RNN backbone) and
TransformerMAF (Transformer backbone). We use a context window size of $200$
to give the model access to the information needed to infer the current
parameter of the true conditional distribution. We tried increasing the window
size to $500$ on VAR(1) – Dynamic but did not see performance improvements. We
also removed the unnecessary default input features of these models to prevent
overfitting. For details of the setup, please see Appendix D. DynaConF is
trained with the autoregressive posterior.

For evaluation we use a rolling-window approach with a window size of $10$
steps. The final evaluation metrics are the aggregated results from all $100$
test windows. For univariate data we report the continuous ranked probability
score (CRPS) [52], a commonly used score to measure how close the predicted
distribution is to the true distribution. For multivariate data we also report
$\textrm{CRPS}_{\Sigma}$, a popular variant of CRPS for multivariate
distributions [3, 29, 4]. Additional results, including an evaluation based on
mean squared error (MSE), can be found in the Appendix E. In all cases lower
values indicate better performance.

#### Results

The results on synthetic data are shown in Table 2. For univariate data, our
full model (DynaConF) outperforms its ablated counterpart (StatiConF) by 2.0%
– 12.3%, validating the importance of our dynamic adaptation to non-starionary
effects. DynaConF – PP is also superior to all univariate baselines, with its
closest competitor DeepAR – 10 behind by an average of 6.1%. Furthermore, we
note that our model with the pointwise encoder tends to outperform the MLP
encoder, both for the ablated and full model. For multivariate data we observe
similar trends. Here, our full model (DynaConF – PP) performs 24.3% (CRPS) /
33.8% ($\textrm{CRPS}_{\Sigma}$) better than the ablated model with static
conditional distribution (StatiConF – PP) and 22.6% (CRPS) / 31.0%
($\textrm{CRPS}_{\Sigma}$) better than the best-performing baseline (DeepVAR –
160). As before, pointwise encoders tend to perform better than MLP encoders.
Since for synthetic data we also have access to the ground-truth models, we
include the corresponding scores as a reference and upper bound in terms of
performance.

Figure 2 shows qualitative results of our model on AR(1) – Flip/Sin/Dynamic.
Note that because the encoder in our model is non-linear, the inferred
$\bm{\phi}_{t}$ does not necessarily correspond to the original parameter.
However, we can clearly see how it differs for Sin (continuous changes) vs
Flip/Dynamic (discrete jumps) and how it tracks the ground-truth up to
scale/sign.

(a) Flip

(b) Sin

(c) Dynamic

Figure 2: Qualitative Results. We show one dimension of $\bm{\phi}_{t}$ of
DynaConF – PP inferred with particle filters at test time on (a) AR(1) – Flip,
(b) AR(1) – Sin, and (c) AR(1) – Dynamic. The red dashed lines show the
ground-truth parameters of the conditional distribution varying over time. The
blue curves and bands are the medians and 90% confidence intervals of the
posterior. Table 2: Quantitative Evaluation (Synthetic Data). We compare
StatiConF/DynaConF to (a) $4$ univariate and (b) $6$ multivariate baselines.
For univariate processes (AR(1)-Flip/Sin/Dynamic) we report CRPS. For
multivariate processes (VAR(1)-Dynamic) we report both CRPS and
$\textrm{CRPS}_{\Sigma}$. Values behind models indicate the number of hidden
units.

| Method | CRPS
---|---|---
| AR(1)-F | AR(1)-S | AR(1)-D
| GroundTruth | 0.731$\scriptscriptstyle\pm\text{0.001}$ | 0.710$\scriptscriptstyle\pm\text{0.001}$ | 0.624$\scriptscriptstyle\pm\text{0.001}$
| DeepAR – 10 | 0.741$\scriptscriptstyle\pm\text{0.005}$ | 0.764$\scriptscriptstyle\pm\text{0.004}$ | 0.768$\scriptscriptstyle\pm\text{0.011}$
| DeepAR – 40 | 0.740$\scriptscriptstyle\pm\text{0.003}$ | 0.776$\scriptscriptstyle\pm\text{0.002}$ | 0.820$\scriptscriptstyle\pm\text{0.053}$
| DeepAR – 160 | 0.740$\scriptscriptstyle\pm\text{0.001}$ | 0.774$\scriptscriptstyle\pm\text{0.004}$ | 0.801$\scriptscriptstyle\pm\text{0.047}$
| DeepSSM | 0.755$\scriptscriptstyle\pm\text{0.001}$ | 0.761$\scriptscriptstyle\pm\text{0.001}$ | 0.803$\scriptscriptstyle\pm\text{0.001}$
(ours) | StatiConF – MLP | 0.753$\scriptscriptstyle\pm\text{0.001}$ | 0.784$\scriptscriptstyle\pm\text{0.002}$ | 0.764$\scriptscriptstyle\pm\text{0.003}$
StatiConF – PP | 0.752$\scriptscriptstyle\pm\text{0.002}$ | 0.763$\scriptscriptstyle\pm\text{0.001}$ | 0.783$\scriptscriptstyle\pm\text{0.002}$
DynaConF – MLP | 0.750$\scriptscriptstyle\pm\text{0.001}$ | 0.727$\scriptscriptstyle\pm\text{0.019}$ | 0.691$\scriptscriptstyle\pm\text{0.006}$
DynaConF – PP | 0.737$\scriptscriptstyle\pm\text{0.001}$ | 0.721$\scriptscriptstyle\pm\text{0.005}$ | 0.687$\scriptscriptstyle\pm\text{0.010}$

[AR(1)-*: F = Flip; S = Sin; D = Dynamic]

(a) Univariate

| Method | CRPS | $\textrm{CRPS}_{\Sigma}$
---|---|---|---
| GroundTruth | 0.496$\scriptscriptstyle\pm\text{0.001}$ | $0.437\scriptscriptstyle\pm\text{0.001}$
| DeepVAR – 10 | 0.797$\scriptscriptstyle\pm\text{0.009}$ | $0.824\scriptscriptstyle\pm\text{0.015}$
| DeepVAR – 40 | 0.792$\scriptscriptstyle\pm\text{0.000}$ | $0.821\scriptscriptstyle\pm\text{0.002}$
| DeepVAR – 160 | 0.787$\scriptscriptstyle\pm\text{0.001}$ | $0.816\scriptscriptstyle\pm\text{0.002}$
| TransformerMAF – 8 | 0.800$\scriptscriptstyle\pm\text{0.001}$ | $0.831\scriptscriptstyle\pm\text{0.008}$
| TransformerMAF – 32 | 0.806$\scriptscriptstyle\pm\text{0.008}$ | $0.855\scriptscriptstyle\pm\text{0.011}$
| TransformerMAF – 128 | 0.866$\scriptscriptstyle\pm\text{0.077}$ | $0.873\scriptscriptstyle\pm\text{0.025}$
(ours) | StatiConF – MLP | 0.806$\scriptscriptstyle\pm\text{0.004}$ | $0.843\scriptscriptstyle\pm\text{0.005}$
StatiConF – PP | 0.805$\scriptscriptstyle\pm\text{0.002}$ | $0.850\scriptscriptstyle\pm\text{0.003}$
DynaConF – MLP | 0.762$\scriptscriptstyle\pm\text{0.036}$ | 0.777$\scriptscriptstyle\pm\text{0.047}$
DynaConF – PP | 0.609$\scriptscriptstyle\pm\text{0.012}$ | 0.563$\scriptscriptstyle\pm\text{0.019}$

(b) Multivariate: VAR(1) – Dynamic

###  4.2 Experiments on Real-World Data

#### Datasets

We evaluate the proposed method on six publicly available datasets [53, 2, 3]:
(Exchange) daily exchange rates of 8 different countries from 1990 to 2016;
(Solar) [54] solar power production in 10-minute intervals of 137 PV plants in
2006; (Electricity) [55] hourly electricity consumption of 370 customers from
2012 to 2014; (Traffic) [56] hourly occupancy data at 963 sensor locations in
the San Francisco Bay area; (Taxi) rides taken in 30-minute intervals at 1214
locations in New York City in January 2015/2016; (Wikipedia) daily page views
of 2000 Wikipedia articles.

#### Experimental Setup

We use the same train/test splits and the same input features, such as time of
the day, as previous works with published code and results [3, 4, 29]. For our
method, we first train StatiConF and then reuse its learned encoder in
DynaConF, so the learning of DynaConF is focused on the dynamic model.
DynaConF is trained with the IAF posterior for efficiency. Our models use a
two-layer LSTM with 128 hidden units as the encoder, except for the
8-dimensional Exchange data, where the hidden size is 8. We stress again that,
different from DeepVAR or LSTM-MAF, we use LSTM as an encoder of
$(\bm{y}_{t-B:t-1},\bm{x}_{t-B,t})$ only, so we actually “restart” it at every
time step. More details of our hyperparameters can be found in Appendix C.

#### Results

The results are shown in Table 3. As we can see, for different datasets and
different evaluation metrics, the relative performance of each method can be
different. This shows the diversity of these datasets and that different
models may benefit from dataset-specific structure in different ways. However,
DynaConF achieves the best performance more often than all the other
baselines. Where it does not outperform, its performance is competitive
consistently, unlike the baselines, whose relative performance (compared to
others) varies significantly across datasets. We also note that our full
model, DynaConF, which adapts to changes in the conditional distribution
consistently outperforms our ablated model, StatiConF, which only models a
static conditional distribution, except for electricity, where the performance
is similar. This shows the effectiveness of modeling the dynamic changes in
the conditional distribution. Further results, including CRPSΣ, and standard
deviations of MSE and CRPS are in Appendix E.

Table 3: Quantitative Evaluation (Real-World Data). We compare the CRPS and
MSE scores of StatiConF/DynaConF and 7 baselines on 6 publicly available
datasets. Lower values are better.

Method | Exchange | Solar | Electricity | Traffic | Taxi | Wikipedia
---|---|---|---|---|---|---
CRPS | MSE | CRPS | MSE | CRPS | MSE | CRPS | MSE | CRPS | MSE | CRPS | MSE
| [e-4] |  | [e+2] |  | [e+5] |  | [e-4] |  | [e+1] |  | [e+7]
DeepVAR (I) | 0.013 | 1.6 | 0.434 | 9.3 | 1.059 | 2.1 | 0.168 | 6.3 | 0.586 | 7.3 | 0.379 | 7.2
DeepVAR (C) | 0.009 | 1.9 | 0.384 | 29 | 0.084 | 55 | 0.165 | 15 | 0.416 | 5.1 | 0.247 | 3.8
GP-scaling | 0.017 | 2.9 | 0.415 | 11 | 0.053 | 1.8 | 0.140 | 5.2 | 0.346 | 2.7 | 1.549 | 5.5
GP-Copula | 0.008 | 1.7 | 0.371 | 9.8 | 0.056 | 2.4 | 0.133 | 6.9 | 0.360 | 3.1 | 0.236 | 4.0
LSTM-NVP | 0.010 | 2.4 | 0.365 | 9.1 | 0.059 | 2.5 | 0.172 | 6.9 | 0.327 | 2.6 | 0.333 | 4.7
LSTM-MAF | 0.012 | 3.8 | 0.378 | 9.8 | 0.051 | 1.8 | 0.124 | 4.9 | 0.314 | 2.4 | 0.282 | 3.8
TransformerMAF | 0.012 | 3.4 | 0.368 | 9.3 | 0.052 | 2.0 | 0.134 | 5.0 | 0.377 | 4.5 | 0.274 | 3.1
StatiConF (ours) | 0.011 | 4.4 | 0.343 | 6.6 | 0.059 | 2.0 | 0.143 | 4.1 | 0.311 | 2.2 | 0.427 | 4.1
DynaConF (ours) | 0.010 | 2.6 | 0.338 | 6.4 | 0.058 | 2.1 | 0.140 | 3.9 | 0.308 | 2.2 | 0.296 | 3.8

##  5 Discussion

In this work, we addressed the problem of modeling and forecasting time series
with non-stationary conditional distributions. We proposed a new model,
DynaConF, that explicitly decouples the time-invariant conditional
distribution modeling and the time-variant non-stationarity modeling. We
designed specific architectures, developed two types of variational
posteriors, and employed Rao-Blackwellized particle filters to allow the model
to train efficiently on large multivariate time series and adapt online at
test time. Results on synthetic and real-world data show that our model can
learn and adapt to different types of changes in the conditional distribution
parameters and perform competitively or better than state-of-the-art time
series forecasting models.

Our model currently has the following limitations. (1) Since the variational
posterior model complexity scales in $O(T)$, it is challenging to train the
model on very long time series. (2) We used a simple observation distribution
family in this work for efficient inference, but there are cases where this
may impact performance. For future work, it would be interesting to develop
new variational posteriors and training algorithms that can scale better and
more flexible inference algorithms that can deal with more complex observation
distributions.

## References

  * Salinas et al. [2020] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski.  DeepAR: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ , 2020.
  * Lai et al. [2018] 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 on Research & Development in Information Retrieval_, 2018.
  * Salinas et al. [2019] David Salinas, Michael Bohlke-Schneider, Laurent Callot, Roberto Medico, and Jan Gasthaus.  High-dimensional multivariate forecasting with low-rank Gaussian Copula Processes.  In _Advances in Neural Information Processing Systems_ , 2019\.
  * Rasul et al. [2021a] Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs M. Bergmann, and Roland Vollgraf.  Multivariate probabilistic time series forecasting via conditioned normalizing flows.  In _International Conference on Learning Representations_ , 2021a.
  * Rangapuram et al. [2018] Syama Sundar Rangapuram, Matthias W. Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski.  Deep state space models for time series forecasting.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Brockwell and Davis [2009] Peter J. Brockwell and Richard A. Davis.  _Time Series: Theory and Methods_.  2009\.
  * Hamilton [1994] James Douglas Hamilton.  _Time Series Analysis_.  1994\.
  * Kwiatkowski et al. [1992] Denis Kwiatkowski, Peter CB Phillips, Peter Schmidt, and Yongcheol Shin.  Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root?  _Journal of Econometrics_ , 1992.
  * Dickey and Fuller [1979] David A. Dickey and Wayne A. Fuller.  Distribution of the estimators for autoregressive time series with a unit root.  _Journal of the American Statistical Association_ , 1979.
  * Kim et al. [2022] Taesung Kim, Jinhee Kim, Yunwon Tae, Cheonbok Park, Jang-Ho Choi, and Jaegul Choo.  Reversible instance normalization for accurate time-series forecasting against distribution shift.  In _International Conference on Learning Representations_ , 2022.
  * Box et al. [2015] George EP Box, Gwilym M. Jenkins, Gregory C. Reinsel, and Greta M. Ljung.  _Time Series Analysis: Forecasting and Control_.  2015\.
  * Engle [1982] Robert F. Engle.  Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation.  _Econometrica: Journal of the Econometric Society_ , 1982.
  * Bollerslev [1986] Tim Bollerslev.  Generalized autoregressive conditional heteroskedasticity.  _Journal of Econometrics_ , 1986.
  * Kalman [1960] Re E Kalman.  A new approach to linear filtering and prediction problems.  _Transactions of the ASME-Journal of Basic Engineering_ , 1960.
  * Hochreiter and Schmidhuber [1997] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural Computation_ , 1997.
  * Bai et al. [2018] Shaojie Bai, J. Zico Kolter, and Vladlen Koltun.  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.  _arXiv preprint arXiv:1803.01271_ , 2018.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in Neural Information Processing Systems_ , 2017.
  * Garnelo et al. [2018a] Marta Garnelo, Dan Rosenbaum, Christopher Maddison, Tiago Ramalho, David Saxton, Murray Shanahan, Yee Whye Teh, Danilo Rezende, and S. M. Ali Eslami.  Conditional neural processes.  In _International Conference on Machine Learning_ , 2018a.
  * Garnelo et al. [2018b] Marta Garnelo, Jonathan Schwarz, Dan Rosenbaum, Fabio Viola, Danilo J. Rezende, S. M. Eslami, and Yee Whye Teh.  Neural processes.  _arXiv preprint arXiv:1807.01622_ , 2018b.
  * Lim and Zohren [2020] Bryan Lim and Stefan Zohren.  Time series forecasting with deep learning: A survey.  _arXiv:2004.13408 [cs, stat]_ , 2020.
  * Holt [2004] Charles C. Holt.  Forecasting seasonals and trends by exponentially weighted moving averages.  _International Journal of Forecasting_ , 2004.
  * Cleveland et al. [1990] Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning.  STL: A seasonal-trend decomposition procedure based on loess.  _Journal of Official Statistics_ , 1990.
  * Smyl [2020] Slawek Smyl.  A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting.  _International Journal of Forecasting_ , 2020.
  * Bandara et al. [2020] Kasun Bandara, Christoph Bergmeir, and Slawek Smyl.  Forecasting across time series databases using recurrent neural networks on groups of similar series: A clustering approach.  _Expert Systems with Applications_ , 2020.
  * Oreshkin et al. [2020] Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.  _arXiv:1905.10437 [cs, stat]_ , 2020.
  * Wu et al. [2021] Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long.  Autoformer: Decomposition transformers with auto-correlation for long-term series forecasting.  In _Advances in Neural Information Processing Systems_ , 2021\.
  * Woo et al. [2022] Gerald Woo, Chenghao Liu, Doyen Sahoo, Akshat Kumar, and Steven Hoi.  ETSformer: Exponential Smoothing Transformers for Time-series Forecasting.  _arXiv preprint arXiv:2202.01381_ , 2022.
  * Corani et al. [2021] Giorgio Corani, Alessio Benavoli, and Marco Zaffalon.  Time series forecasting with Gaussian Processes needs priors.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , 2021.
  * Rasul et al. [2021b] Kashif Rasul, Calvin Seward, Ingmar Schuster, and Roland Vollgraf.  Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting.  In _Proceedings of the 38th International Conference on Machine Learning_ , 2021b.
  * Fraccaro et al. [2017] Marco Fraccaro, Simon Kamronn, Ulrich Paquet, and Ole Winther.  A disentangled recognition and nonlinear dynamics model for unsupervised learning.  In _Advances in Neural Information Processing Systems_ , 2017\.
  * de Bézenac et al. [2020] 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 H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , 2020.
  * Tang and Matteson [2021] Binh Tang and David S Matteson.  Probabilistic transformer for time series analysis.  In _Advances in Neural Information Processing Systems_ , 2021\.
  * Klushyn et al. [2021] Alexej Klushyn, Richard Kurle, Maximilian Soelch, Botond Cseke, and Patrick van der Smagt.  Latent matters: Learning deep state-space models.  In _Advances in Neural Information Processing Systems_ , 2021\.
  * Ansari et al. [2021] Abdul Fatir Ansari, Konstantinos Benidis, Richard Kurle, Ali Caner Turkmen, Harold Soh, Alexander J. Smola, Bernie Wang, and Tim Januschowski.  Deep explicit duration switching models for time series.  _Advances in Neural Information Processing Systems_ , 2021.
  * Kurle et al. [2020] Richard Kurle, Syama Sundar Rangapuram, Emmanuel de Bézenac, Stephan Günnemann, and Jan Gasthaus.  Deep Rao-Blackwellised particle filters for time series forecasting.  _Advances in Neural Information Processing Systems_ , 2020.
  * Nguyen and Quanz [2021] Nam Nguyen and Brian Quanz.  Temporal latent auto-encoder: A method for probabilistic multivariate time series forecasting.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 2021.
  * Le Guen and Thome [2020] Vincent Le Guen and Nicolas Thome.  Probabilistic time series forecasting with shape and temporal diversity.  _Advances in Neural Information Processing Systems_ , 2020.
  * Yu et al. [2021] Zhongjie Yu, Fabrizio G. Ventola, and Kristian Kersting.  Whittle networks: A deep likelihood model for time series.  In _International Conference on Machine Learning_ , 2021.
  * De Lange et al. [2021] Matthias De Lange, Rahaf Aljundi, Marc Masana, Sarah Parisot, Xu Jia, Ales Leonardis, Gregory Slabaugh, and Tinne Tuytelaars.  A continual learning survey: Defying forgetting in classification tasks.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2021.
  * Parisi et al. [2019] German I. Parisi, Ronald Kemker, Jose L. Part, Christopher Kanan, and Stefan Wermter.  Continual lifelong learning with neural networks: A review.  _Neural Networks_ , 2019.
  * Kirkpatrick et al. [2017] James Kirkpatrick, Razvan Pascanu, Neil Rabinowitz, Joel Veness, Guillaume Desjardins, Andrei A. Rusu, Kieran Milan, John Quan, Tiago Ramalho, Agnieszka Grabska-Barwinska, Demis Hassabis, Claudia Clopath, Dharshan Kumaran, and Raia Hadsell.  Overcoming catastrophic forgetting in neural networks.  _arXiv:1612.00796 [cs, stat]_ , 2017.
  * Gupta et al. [2022] Vibhor Gupta, Jyoti Narwariya, Pankaj Malhotra, Lovekesh Vig, and Gautam Shroff.  Continual learning for multivariate time series tasks with variable input dimensions.  _arXiv:2203.06852 [cs]_ , 2022.
  * Nguyen et al. [2018] Cuong V. Nguyen, Yingzhen Li, Thang D. Bui, and Richard E. Turner.  Variational continual learning.  In _International Conference on Learning Representations_ , 2018.
  * Kurle et al. [2019] Richard Kurle, Botond Cseke, Alexej Klushyn, Patrick van der Smagt, and Stephan Günnemann.  Continual learning with bayesian neural networks for non-stationary data.  In _International Conference on Learning Representations_ , 2019.
  * Kingma and Welling [2013] Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Kingma et al. [2017] Diederik P. Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling.  Improving variational inference with inverse autoregressive flow.  _arXiv:1606.04934 [cs, stat]_ , 2017.
  * Germain et al. [2015] Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle.  MADE: Masked autoencoder for distribution estimation.  In _Proceedings of the 32nd International Conference on Machine Learning_ , 2015.
  * Doucet et al. [2000] Arnaud Doucet, Nando de Freitas, Kevin Murphy, and Stuart Russell.  Rao-Blackwellised particle filtering for dynamic Bayesian networks.  In _Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence_ , UAI’00, 2000.
  * Alexandrov et al. [2020] Alexander Alexandrov, Konstantinos Benidis, Michael Bohlke-Schneider, Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Danielle C. Maddix, Syama Rangapuram, David Salinas, Jasper Schulz, Lorenzo Stella, Ali Caner Türkmen, and Yuyang Wang.  GluonTS: Probabilistic and Neural Time Series Modeling in Python.  _Journal of Machine Learning Research_ , 2020.
  * glu [2021] GluonTS, 2021.  URL https://github.com/awslabs/gluon-ts.
  * [51] Kashif Rasul.  PytorchTS.  URL https://github.com/zalandoresearch/pytorch-ts.
  * Matheson and Winkler [1976] James E. Matheson and Robert L. Winkler.  Scoring rules for continuous probability distributions.  _Management Science_ , 1976.
  * [53] Public time series datasets.  URL https://github.com/mbohlkeschneider/gluon-ts/tree/mv_release/datasets.
  * [54] Solar power dataset.  URL http://www.nrel.gov/grid/solar-power-data.html.
  * [55] Electricity dataset.  URL https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014.
  * [56] Traffic dataset.  URL http://pems.dot.ca.gov.
  * Smith and Topin [2017] Leslie N. Smith and Nicholay Topin.  Super-convergence: Very fast training of neural networks using large learning rates.  _arXiv preprint arXiv:1708.07120_ , 2017.
  * Kingma and Ba [2014] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.

##  Appendix A Normalizing-Flow-Based Variational Posterior

We develop an alternative variational posterior model based on Inverse
Autoregressive Flows (IAFs; [46]). In this model,
$\mathrm{q}(\bm{\chi}_{B:T})$ can be sampled jointly over $t=B:T$ in parallel
as follows. We again let $\mathrm{q}(\bm{\chi}_{B})=\mathrm{p}(\bm{\chi}_{B})$
and only use IAF to model $\mathrm{q}(\bm{\chi}_{B+1:T})$. Let
$T^{\prime}=T-B$. To sample the $i$-th dimension of $\bm{\chi}_{B+1:T}$, we
first sample a vector $\bm{z}_{0}$ of dimension $T^{\prime}$ from a standard
multivariate normal distribution

$\bm{z}_{0}\sim\mathcal{N}(\bm{0},\bm{I}).$ (14)

Then, this sample is transformed sequentially for $l=1,\ldots,L$ as

$\bm{z}_{l}=\bm{\mu}_{l}+\bm{\sigma}_{l}\odot\bm{z}_{l-1},$ (15)

where $\bm{\mu}_{l}$ and $\bm{\sigma}_{l}$ are the output from a neural
network, specifically MADE [47], with a soft-plus transformation on the latter
to make sure it is positive. The final output $\bm{z}_{L}$ is taken as the
sample of the $i$-th dimension of $\bm{\chi}_{B+1:T}$. The input of the MADE
at the $l$-th iteration consists of not only the previous output
$\bm{z}_{l-1}$ but also an additional learnable embedding for
$\bm{\chi}_{B+1:T,i}$, the $i$-th dimension of $\bm{\chi}_{B+1:T}$. We set
$L=3$ and use MADE of 2 hidden layers of size 1024 and learnable embeddings of
size 512 for all our experiments.

To compute $\mathrm{q}(\bm{\chi}_{B+1:T})$, we also need to compute the
determinants of the Jacobians
$|{\mathrm{d}\bm{z}_{l}}/{\mathrm{d}\bm{z}_{l-1}}|,l=1,\ldots,L$, in addition
to the density of the base distribution, which is the standard normal
$\mathcal{N}(\bm{0},\bm{I})$. Due to the structure of MADE and the definition
of IAF (Eq. 15), the Jacobian determinants can be easily computed from
$\bm{\sigma}_{l}$.

##  Appendix B Optimization Procedure

In contrast to existing time series models, we utilize a flexible variational
posterior with a large number of parameters in the order of $O(T)$ and a
structured prior model to account for conditional distribution changes over
time. However, jointly optimizing over these variational parameters and the
parameters in the conditional distribution model itself using stochastic
gradient descent can be prohibitively demanding on the computational
resources, especially GPU memory. Instead, we propose an alternative
optimization procedure to learn these parameters.

Specifically, instead of optimizing by stochastic gradient descent (SGD) over
all parameters in both the conditional distribution model and the prior and
variational posterior model, we propose to learn the model by alternating the
optimization of the parameters in the conditional distribution model and the
prior and variational posterior models. For the former, we condition on the
samples from the current posterior model while optimizing the conditional
distribution model on randomly sampled sub-sequences from the time series
using SGD. For the latter, we fix the conditional distribution model, sample
from the variational posterior model over the entire time series either
sequentially or in parallel, depending on the variational posterior model we
use, and compute the loss using the samples through the conditional
distribution, prior, and posterior models. Then we perform an update on the
prior and posterior model parameters using the gradient from the loss.

When the time series is high-dimensional and GPU memory becomes a constraint
during training, we randomly sample a subset of _observation_ dimensions for
each batch, since the loss decomposes over the observation dimensions in our
model.

##  Appendix C Hyperparameters

On synthetic data, we use the validation set to early stop and choose the best
model for both the baselines and StatiConF. For DynaConF, we keep training as
long as the loss is decreasing on the training set. We perform 50 updates per
epoch. We use 32 hidden units for our 2-layer MLP encoder. For the baselines,
we report their results with different hidden sizes, including the default
ones.

On real-world data, we use the “1cycle” learning rate scheduler [57] and
manually choose the number of epochs for StatiConF. After training StatiConF,
we reuse its encoders in DynaConF, so it only needs to learn the dynamic
model. For DynaConF, we train until the loss on the training set converges.

For our models, we use Adam [58] as the optimizer with the default learning
rate of $0.001$ unless the learning rate is controlled by the “1cycle”
learning rate scheduler. The dimension of the latent vector $\bm{z}_{t,i}$
(see Section 3.2) is set to $E=4$ across all the experiments.

Because of the diversity of the real-world datasets, we further apply
techniques to stablize training. Specifically, for all datasets, we use the
mean and standard deviation to shift and scale each dimension of the time
series. For Exchange, we use the mean and standard deviation of the recent
past data in a moving context window. For the other datasets, we simply use
the global mean and standard deviation of each dimension computed using the
whole training set. In all cases, for forecasting, the output from the model
is inversely scaled and shifted back for evaluation.

Different from the rest of the datasets, Taxi and Wikipedia consist of counts.
The time series in Taxi are relatively small counts, while the time series in
Wikipedia are very large and can have extreme values. To stabilize training,
we apply the dequantization technique used in [4] to Taxi, where we add a
small noise from the uniform distribution $\mathcal{U}(-0.5,0.5)$ to the
target time series $\bm{y}$. On Wikipedia, we use the quantiles (0.02 and
0.95) computed from the recent past data in a moving context window to
Winsorize extreme values.

##  Appendix D Experiment Setup Details

On synthetic data, we keep most baseline hyperparameters to their default
values but make the following changes to account for properties of our
synthetic data: (1) To reduce overfitting we remove any unnecessary input
features from the models and use only past observations with time lag $1$ as
input. (2) To allow the models to adapt to changes in the conditional
distribution we increase the context window size to $200$. This allows the
models to see enough observations generated with the latest ground-truth
distribution parameters, so the models have the necessary information to adapt
to the current distribution. For VAR(1) – Dynamic, we also tried extending it
to $500$, but it did not improve the performance. (3) DeepSSM allows modeling
of trend and seasonality, but since our synthetic data do not have those, we
explicitly remove those components from the model specification to avoid
overfitting; (4) DeepVAR allows modeling of different covariance structures in
the noise, such as diagonal, low-rank, and full-rank. Since our synthetic data
follow a diagonal covariance structure in the noise, we explicitly specify
that in DeepVAR.

We run all experiments for three different random seeds independently and
calculate and report the mean and standard deviation of each evaluation metric
for each model. On synthetic datasets, we use 1000 sample paths to empirically
estimate the predicted distributions for all models. On real-world datasets,
we use 100 sample paths.

We use three evaluation metrics: mean squared error (MSE), continuous ranked
probability score (CRPS)[52] and CRPSΣ. Assume that we observe $y$ at time $t$
but a probabilistic forecasting model predicts the distribution of $y$ to be
$F$. MSE is widely used for time series forecasting, and for a probabilistic
forecasting model, where the mean of the distribution is used for point
prediction, it is defined as

$\text{MSE}(F,y)=(\mathrm{E}_{z\sim F}[z]-y)^{2}$ (16)

for a single time point $t$ and averaged over all the time points in the test
set.

CRPS has been used for evaluating how close the predicted probability
distribution is to the ground-truth distribution and is defined as

$\text{CRPS}(F,y)=\int_{\mathbb{R}}(F(z)-\mathbb{I}[y\leq
z])^{2}\mathrm{d}{z},$ (17)

for a single time point $t$ and averaged over all the time points in the test
set, where $\mathbb{I}$ denotes the indicator function. Generally, $F(z)$ can
be approximated by the empirical distribution of the samples from the
predicted distribution. Both MSE and CRPS can be applied to multivariate time
series by computing the metric on each dimension and then averaging over all
the dimensions.

CRPSΣ is another metric that has been used in recent works [3, 4, 29] to
evaluate multivariate probabilistc forecasting results. To compute it, the sum
of the time series across all the dimensions is computed and then compared
with the sum of the prediction using CRPS. That is

$\text{CRPS}_{\Sigma}(F_{\Sigma},\sum_{i}y_{i})=\int_{\mathbb{R}}(F_{\Sigma}(z)-\mathbb{I}[\sum_{i}y_{i}\leq
z])^{2}\mathrm{d}{z},$ (18)

where $F_{\Sigma}$ is the distribution of the sum of the predicted values
$\sum_{i}z_{i},z_{i}\sim F_{i}(z)$ and $i$ denotes the dimension $i$ of the
time series. $F_{\Sigma}$ is usually approximated by the empirical
distribution of the samples summed across the dimensions.

##  Appendix E Additional Experiment Results

Table 4 and 5 show the MSE results of the baselines and our models on the
univariate and multivariate processes respectively. Table 6 shows the CRPSΣ
scores of our models on the real-world datasets, while Table 7 and 8 show the
full CRPS and MSE results of our models with means and standard deviations.
The results of the baselines on the real-world data are from [4, 29].

Table 4: Quantitative Evaluation (Synthetic Data). MSE results on univariate processes AR(1)-Flip/Sin/Dynamic. Values behind architectures indicate the number of hidden units. | Method | MSE
---|---|---
| AR(1)-Flip | AR(1)-Sin | AR(1)-Dynamic
| GroundTruth | 1.2$\scriptscriptstyle\pm\text{9.2e-4}$ | 1.6$\scriptscriptstyle\pm\text{3.0e-3}$ | 2.1$\scriptscriptstyle\pm\text{5.3e-3}$
| DeepAR – 10 | 1.3$\scriptscriptstyle\pm\text{1.7e-2}$ | 1.8$\scriptscriptstyle\pm\text{1.8e-2}$ | 3.2$\scriptscriptstyle\pm\text{7.2e-2}$
| DeepAR – 40 | 1.3$\scriptscriptstyle\pm\text{8.9e-3}$ | 1.8$\scriptscriptstyle\pm\text{1.3e-2}$ | 3.6$\scriptscriptstyle\pm\text{3.8e-1}$
| DeepAR – 160 | 1.3$\scriptscriptstyle\pm\text{2.8e-3}$ | 1.8$\scriptscriptstyle\pm\text{9.6e-3}$ | 3.5$\scriptscriptstyle\pm\text{3.3e-1}$
| DeepSSM | 1.3$\scriptscriptstyle\pm\text{2.2e-3}$ | 1.8$\scriptscriptstyle\pm\text{2.8e-3}$ | 3.3$\scriptscriptstyle\pm\text{3.0e-3}$
(ours) | StatiConF – MLP | 1.3$\scriptscriptstyle\pm\text{4.0e-3}$ | 1.9$\scriptscriptstyle\pm\text{7.7e-3}$ | 3.3$\scriptscriptstyle\pm\text{3.5e-2}$
StatiConF – PP | 1.3$\scriptscriptstyle\pm\text{4.7e-3}$ | 1.8$\scriptscriptstyle\pm\text{5.4e-3}$ | 3.3$\scriptscriptstyle\pm\text{6.2e-3}$
DynaConF – MLP | 1.3$\scriptscriptstyle\pm\text{4.0e-3}$ | 1.6$\scriptscriptstyle\pm\text{9.2e-2}$ | 2.6$\scriptscriptstyle\pm\text{5.2e-2}$
DynaConF – PP | 1.2$\scriptscriptstyle\pm\text{5.5e-3}$ | 1.6$\scriptscriptstyle\pm\text{2.7e-2}$ | 2.6$\scriptscriptstyle\pm\text{5.1e-2}$
Table 5: Quantitative Evaluation (Synthetic Data). MSE results on multivariate process VAR(1)-Dynamic. Values behind architectures indicate the number of hidden units. | Method | MSE
---|---|---
| GroundTruth | 2.8$\scriptscriptstyle\pm\text{7.6e-3}$
| DeepVAR – 10 | 8.4$\scriptscriptstyle\pm\text{5.6e-2}$
| DeepVAR – 40 | 8.4$\scriptscriptstyle\pm\text{2.7e-3}$
| DeepVAR – 160 | 8.3$\scriptscriptstyle\pm\text{2.2e-2}$
| TransformerMAF – 8 | 8.5$\scriptscriptstyle\pm\text{2.9e-2}$
| TransformerMAF – 32 | 8.5$\scriptscriptstyle\pm\text{3.7e-2}$
| TransformerMAF – 128 | 9.4$\scriptscriptstyle\pm\text{1.2e0}$
(ours) | StatiConF – MLP | 8.5$\scriptscriptstyle\pm\text{8.5e-2}$
StatiConF – PP | 8.5$\scriptscriptstyle\pm\text{2.3e-2}$
DynaConF – MLP | 7.9$\scriptscriptstyle\pm\text{5.7e-1}$
DynaConF – PP | 4.5$\scriptscriptstyle\pm\text{3.0e-1}$
Table 6: Quantitative Evaluation (Real-World Data). Full CRPSΣ results of StatiConF/DynaConF and 8 baselines with means and standard deviations on 6 publicly available datasets. Lower values are better. Method | Exchange | Solar | Electricity | Traffic | Taxi | Wikipedia
---|---|---|---|---|---|---
DeepVAR (I) | $0.008\scriptscriptstyle\pm\text{0.001}$ | $0.391\scriptscriptstyle\pm\text{0.017}$ | $0.025\scriptscriptstyle\pm\text{0.001}$ | $0.087\scriptscriptstyle\pm\text{0.041}$ | $0.506\scriptscriptstyle\pm\text{0.005}$ | $0.133\scriptscriptstyle\pm\text{0.002}$
DeepVAR (C) | $0.007\scriptscriptstyle\pm\text{0.000}$ | $0.319\scriptscriptstyle\pm\text{0.011}$ | $0.064\scriptscriptstyle\pm\text{0.008}$ | $0.103\scriptscriptstyle\pm\text{0.006}$ | $0.326\scriptscriptstyle\pm\text{0.007}$ | $0.241\scriptscriptstyle\pm\text{0.033}$
GP-Scaling | $0.009\scriptscriptstyle\pm\text{0.000}$ | $0.368\scriptscriptstyle\pm\text{0.012}$ | $0.022\scriptscriptstyle\pm\text{0.000}$ | $0.079\scriptscriptstyle\pm\text{0.000}$ | $0.183\scriptscriptstyle\pm\text{0.395}$ | $1.483\scriptscriptstyle\pm\text{1.034}$
GP-Copula | $0.007\scriptscriptstyle\pm\text{0.000}$ | $0.337\scriptscriptstyle\pm\text{0.024}$ | $0.024\scriptscriptstyle\pm\text{0.002}$ | $0.078\scriptscriptstyle\pm\text{0.002}$ | $0.208\scriptscriptstyle\pm\text{0.183}$ | $0.086\scriptscriptstyle\pm\text{0.004}$
LSTM-NVP | $0.006\scriptscriptstyle\pm\text{0.003}$ | $0.331\scriptscriptstyle\pm\text{0.020}$ | $0.024\scriptscriptstyle\pm\text{0.001}$ | $0.078\scriptscriptstyle\pm\text{0.001}$ | $0.175\scriptscriptstyle\pm\text{0.001}$ | $0.078\scriptscriptstyle\pm\text{0.001}$
LSTM-MAF | 0.005$\scriptscriptstyle\pm\text{0.003}$ | $0.315\scriptscriptstyle\pm\text{0.023}$ | $0.021\scriptscriptstyle\pm\text{0.000}$ | $0.069\scriptscriptstyle\pm\text{0.002}$ | $0.161\scriptscriptstyle\pm\text{0.002}$ | $0.067\scriptscriptstyle\pm\text{0.001}$
TransformerMAF | 0.005$\scriptscriptstyle\pm\text{0.003}$ | $0.301\scriptscriptstyle\pm\text{0.014}$ | $0.021\scriptscriptstyle\pm\text{0.000}$ | $0.056\scriptscriptstyle\pm\text{0.001}$ | $0.179\scriptscriptstyle\pm\text{0.002}$ | $0.063\scriptscriptstyle\pm\text{0.003}$
TimeGrad | $0.006\scriptscriptstyle\pm\text{0.001}$ | $0.287\scriptscriptstyle\pm\text{0.020}$ | $0.021\scriptscriptstyle\pm\text{0.001}$ | $0.044\scriptscriptstyle\pm\text{0.006}$ | 0.114$\scriptscriptstyle\pm\text{0.020}$ | 0.049$\scriptscriptstyle\pm\text{0.002}$
StatiConF | 0.006$\scriptscriptstyle\pm\text{0.001}$ | 0.255$\scriptscriptstyle\pm\text{0.045}$ | 0.020$\scriptscriptstyle\pm\text{0.004}$ | 0.033$\scriptscriptstyle\pm\text{0.002}$ | 0.149$\scriptscriptstyle\pm\text{0.013}$ | 0.069$\scriptscriptstyle\pm\text{0.009}$
DynaConF | 0.006$\scriptscriptstyle\pm\text{0.001}$ | 0.242$\scriptscriptstyle\pm\text{0.052}$ | 0.021$\scriptscriptstyle\pm\text{0.003}$ | 0.032$\scriptscriptstyle\pm\text{0.002}$ | 0.146$\scriptscriptstyle\pm\text{0.014}$ | 0.081$\scriptscriptstyle\pm\text{0.006}$
Table 7: Quantitative Evaluation (Real-World Data).. Full CRPS results with means and standard deviations on 6 publicly available datasets. Lower values are better. Method | Exchange | Solar | Electricity | Traffic | Taxi | Wikipedia
---|---|---|---|---|---|---
DeepVAR (I) | 0.013$\scriptscriptstyle\pm\text{0.000}$ | 0.434$\scriptscriptstyle\pm\text{0.012}$ | 1.059$\scriptscriptstyle\pm\text{0.001}$ | 0.168$\scriptscriptstyle\pm\text{0.037}$ | 0.586$\scriptscriptstyle\pm\text{0.004}$ | 0.379$\scriptscriptstyle\pm\text{0.004}$
DeepVAR (C) | 0.009$\scriptscriptstyle\pm\text{0.000}$ | 0.384$\scriptscriptstyle\pm\text{0.010}$ | 0.084$\scriptscriptstyle\pm\text{0.006}$ | 0.165$\scriptscriptstyle\pm\text{0.004}$ | 0.416$\scriptscriptstyle\pm\text{0.004}$ | 0.247$\scriptscriptstyle\pm\text{0.001}$
GP-Scaling | 0.017$\scriptscriptstyle\pm\text{0.000}$ | 0.415$\scriptscriptstyle\pm\text{0.009}$ | 0.053$\scriptscriptstyle\pm\text{0.000}$ | 0.140$\scriptscriptstyle\pm\text{0.002}$ | 0.346$\scriptscriptstyle\pm\text{0.348}$ | 1.549$\scriptscriptstyle\pm\text{1.017}$
GP-Copula | 0.008$\scriptscriptstyle\pm\text{0.000}$ | 0.371$\scriptscriptstyle\pm\text{0.022}$ | 0.056$\scriptscriptstyle\pm\text{0.002}$ | 0.133$\scriptscriptstyle\pm\text{0.001}$ | 0.360$\scriptscriptstyle\pm\text{0.201}$ | 0.236$\scriptscriptstyle\pm\text{0.000}$
LSTM-NVP | 0.010$\scriptscriptstyle\pm\text{0.001}$ | 0.365$\scriptscriptstyle\pm\text{0.020}$ | 0.059$\scriptscriptstyle\pm\text{0.001}$ | 0.172$\scriptscriptstyle\pm\text{0.001}$ | 0.327$\scriptscriptstyle\pm\text{0.001}$ | 0.333$\scriptscriptstyle\pm\text{0.001}$
LSTM-MAF | 0.012$\scriptscriptstyle\pm\text{0.003}$ | 0.378$\scriptscriptstyle\pm\text{0.032}$ | 0.051$\scriptscriptstyle\pm\text{0.000}$ | 0.124$\scriptscriptstyle\pm\text{0.002}$ | 0.314$\scriptscriptstyle\pm\text{0.003}$ | 0.282$\scriptscriptstyle\pm\text{0.002}$
TransformerMAF | 0.012$\scriptscriptstyle\pm\text{0.003}$ | 0.368$\scriptscriptstyle\pm\text{0.001}$ | 0.052$\scriptscriptstyle\pm\text{0.000}$ | 0.134$\scriptscriptstyle\pm\text{0.001}$ | 0.377$\scriptscriptstyle\pm\text{0.002}$ | 0.274$\scriptscriptstyle\pm\text{0.007}$
StatiConF | 0.011$\scriptscriptstyle\pm\text{0.001}$ | 0.343$\scriptscriptstyle\pm\text{0.031}$ | 0.059$\scriptscriptstyle\pm\text{0.003}$ | 0.143$\scriptscriptstyle\pm\text{0.001}$ | 0.311$\scriptscriptstyle\pm\text{0.006}$ | 0.427$\scriptscriptstyle\pm\text{0.001}$
DynaConF | 0.010$\scriptscriptstyle\pm\text{0.000}$ | 0.338$\scriptscriptstyle\pm\text{0.035}$ | 0.058$\scriptscriptstyle\pm\text{0.004}$ | 0.140$\scriptscriptstyle\pm\text{0.002}$ | 0.308$\scriptscriptstyle\pm\text{0.005}$ | 0.296$\scriptscriptstyle\pm\text{0.009}$
Table 8: Quantitative Evaluation (Real-World Data).. Full MSE results with means and standard deviations on 6 publicly available datasets. Lower values are better. Method | Exchange | Solar | Electricity | Traffic | Taxi | Wikipedia
---|---|---|---|---|---|---
[e-4] | [e+2] | [e+5] | [e-4] | [e+1] | [e+7]
DeepVAR (I) | 1.6 | 9.3 | 2.1 | 6.3 | 7.3 | 7.2
DeepVAR (C) | 1.9 | 29 | 55 | 15 | 5.1 | 3.8
GP-Scaling | 2.9 | 11 | 1.8 | 5.2 | 2.7 | 5.5
GP-Copula | 1.7 | 9.8 | 2.4 | 6.9 | 3.1 | 4.0
LSTM-NVP | 2.4 | 9.1 | 2.5 | 6.9 | 2.6 | 4.7
LSTM-MAF | 3.8 | 9.8 | 1.8 | 4.9 | 2.4 | 3.8
TransformerMAF | 3.4 | 9.3 | 2.0 | 5.0 | 4.5 | 3.1
StatiConF | 4.4$\scriptscriptstyle\pm\text{2.0e-4}$ | 6.6$\scriptscriptstyle\pm\text{1.2e2}$ | 2.0$\scriptscriptstyle\pm\text{2.8e4}$ | 4.1$\scriptscriptstyle\pm\text{3.1e-6}$ | 2.2$\scriptscriptstyle\pm\text{7.1e-1}$ | 4.1$\scriptscriptstyle\pm\text{4.6e5}$
DynaConF | 2.6$\scriptscriptstyle\pm\text{1.6e-5}$ | 6.4$\scriptscriptstyle\pm\text{1.3e2}$ | 2.1$\scriptscriptstyle\pm\text{5.1e4}$ | 3.9$\scriptscriptstyle\pm\text{4.0e-6}$ | 2.2$\scriptscriptstyle\pm\text{6.9e-1}$ | 3.8$\scriptscriptstyle\pm\text{1.6e5}$