# Deep Probabilistic Time Series Forecasting using Augmented Recurrent Input
for Dynamic Systems

Haitao Liu<EMAIL_ADDRESS>Changjun Liu<EMAIL_ADDRESS>Xiaomo Jiang<EMAIL_ADDRESS>Xudong Chen<EMAIL_ADDRESS>Shuhua Yang<EMAIL_ADDRESS>Xiaofang Wang<EMAIL_ADDRESS>School of
Energy and Power Engineering, Dalian University of Technology, China, 116024
Digital Twin Laboratory for Industrial Equipment, Dalian University of
Technology, China, 116024.  Shenyang Blower Works Group Corporation, China,
110869.

###### Abstract

The demand of probabilistic time series forecasting has been recently raised
in various dynamic system scenarios, for example, system identification and
prognostic and health management of machines. To this end, we combine the
advances in both deep generative models and state space model (SSM) to come up
with a novel, data-driven deep probabilistic sequence model. Specifically, we
follow the popular encoder-decoder generative structure to build the recurrent
neural networks (RNN) assisted variational sequence model on an augmented
recurrent input space, which could induce rich stochastic sequence dependency.
Besides, in order to alleviate the inconsistency issue of the posterior
between training and predicting as well as improving the mining of dynamic
patterns, we (i) propose using a lagged hybrid output as input for the
posterior at next time step, which brings training and predicting into
alignment; and (ii) further devise a generalized auto-regressive strategy that
encodes all the historical dependencies for the posterior. Thereafter, we
first investigate the methodological characteristics of the proposed deep
probabilistic sequence model on toy cases, and then comprehensively
demonstrate the superiority of our model against existing deep probabilistic
SSM models through extensive numerical experiments on eight system
identification benchmarks from various dynamic systems. Finally, we apply our
sequence model to a real-world centrifugal compressor forecasting problem, and
again verify its outstanding performance by quantifying the time series
predictive distribution.

###### keywords:

State space model , Dynamic system , Recurrent neural networks , Variational
inference , Variational autoencoder , Compressor

††journal: Mechanical Systems and Signal Processing

##  1 Introduction

Time series forecasting is a long standing problem in widespread decision-
making scenarios, for example, dynamic system identification [1, 2],
prognostic and health management (PHM) of machines [3, 4, 5], and business
demand forecasting [6]. It thus has gained extensive attention from academic
and industrial community over decades.

In order to learn the underlying dynamic and temporal pattern from historical
time series data and perform desirable forecasting in the future, classic time
series methodologies, for example, exponential smoothing [7], and auto-
regressive integrated moving average (ARIMA) model [8, 9], have been studied
for a long history. They usually incorporate user’s prior knowledge by
decomposing the time series structure into trend, seasonality and so on.
Consequently, these methods have high interpretability while suffering from
poor capability of learning complex and long time dependency. Moreover, state
space model (SSM) [7, 10, 11], which is similar to hidden Markov models [12]
and could include the above classic methods as special cases [9], offers a
more principled way to learn dynamic patterns from data.

For a forecasting model, an important feature is that it could not only
provide the prediction mean but also quantify the associated uncertainty,
which is essential for risk management, for instance, the estimation of
remaining useful life [13]. To this end, an alternative way is interpreting
the SSM paradigm from Bayesian view, for example, the well-known Kalman filter
[14] and the dynamic linear models [15]. Besides, it is popular to combine SSM
with Gaussian process (GP) [16], which is theoretically equivalent to Kalman
filter when using specific kernels. The Gaussian process state space model
(GPSSM) [17, 18] places GP prior over both the transition and measurement
functions of SSM, thus being capable of exporting probabilistic forecasting.
The prominent weakness of GPSSM however is the poor scalability to tackle a
large amount of time series data, due to the cubic time complexity [19].
Though we could resort to the sparse and distributed approximations [20, 21,
22, 23] to improve the scalability of GPSSM, it on the other hand may limit
the model capability due to the approximations.

Alternatively, the deep learning community has followed the SSM framework to
extensively exploit the recurrent neural networks (RNNs), for example, gated
recurrent unit (GRU) [24], long short term memory (LSTM) [25] and their
variants [26, 27], in order to perform time series point prediction. Due to
their high capability of extracting high-dimensional features for learning
complicated patterns within time series data, the RNNs have reported
remarkable and success stories in various domains [28, 29]. In order to enable
probabilistic forecasting using RNN, it is suggested to combine RNN with
generative models. For example, we could combine the auto-regressive strategy
with RNN and train the model through maximum likelihood [30, 31].
Alternatively, we could adopt RNN to learn the time-varying coefficients of
linear/nonlinear stochastic SSM [32, 33], thus enjoying both probabilistic
forecasting and good interpretability. Another line of work builds the
generative GP on the top of RNN as measurement function to produce stochastic
outputs [34, 35]. Finally, combining RNN with variational autoencoder (VAE)
[36] has been exploited recently to build connections between SSM and VAE [37,
38, 39, 40, 31, 41]. These models encode the sequence characteristics and
dependency on the manifold embedded in a low-dimensional latent space, and are
usually trained through variational inference (VI) by deriving the evidence
lower bound (ELBO) as objective, thus resulting in the stochastic or
variational RNN.

Along the line of deep SSM which gains benefits from two communities, we
devise a novel, data-driven deep variational sequence model in an augmented
recurrent input space, named VRNNaug, to perform probabilistic time series
forecasting. The main contributions of this paper are three-fold.

  * 1.

We propose a VAE-type sequence model on an augmented input space, the inputs
of which consist of the original input signal and the additional latent input
which encodes all the sequence characteristics. The augmented input space
could induce rich stochastic sequence dependency.

  * 2.

We alleviate the inconsistency issue of the posterior between training and
predicting as well as improving the modeling of dynamic patterns in the VI
framework by (i) feeding the lagged hybrid output as input for the posterior
at next time step, which brings training and predicting into alignment; and
(ii) presenting a generalized auto-regressive strategy which encodes all the
historical dependencies for the posterior. This is found to greatly improve
the quality of future forecasting.

  * 3.

We perform extensive numerical experiments to verify the superiority of our
model on eight system identification benchmarks from various dynamic systems,
and thereafter apply the proposed sequence model on a real-world centrifugal
compressor case to illustrate the desirable time series forecasting.

The remainder of this paper is organized as follows. Section 2 first defines
the time series forecasting problem, and then section 3 presents the proposed
deep variational sequence model based on augmented recurrent input space.
Thereafter, section 4 conducts a comprehensive comparison study on eight
system identification benchmarks as well as a centrifugal compressor case to
verify the superiority of proposed sequence model. Finally, concluding remarks
of this paper are provided in section 5. Note that the acronyms and notations
used in this paper are summarized in Appendices A and B.

##  2 Preliminaries

It is assumed that the sequence data
$\mathcal{D}=\\{\mathbf{u}_{1:T}\in\mathbb{R}^{T\times
d_{\mathbf{u}}},\mathbf{y}_{1:T}\in\mathbb{R}^{T\times d_{\mathbf{y}}}\\}$ at
$T$ discrete and evenly distributed time points are sampled from an unknown
dynamic system. We attempt to model it in order to learn the underlying
dynamic and temporal patterns. Thereafter, we use the fitted sequence model
$\mathcal{M}$ to forecast the predictive distribution
$p(\mathbf{y}_{T+1:T+F}|\mathcal{D},\mathbf{u}_{T+1:T+F})$ of future window
$\mathbf{y}_{T+1:T+F}$ given the history data $\mathcal{D}$, the forecast
horizon $F$, and the further input signals (also known as covariates or
exogenous variables) $\mathbf{u}_{T+1:T+F}$.

Specifically, we assume that the target dynamic system could be described by
the state space model as

$\displaystyle\mathbf{h}_{t}=$ $\displaystyle
f(\mathbf{h}_{t-1},\mathbf{u}_{t}),$ (1) $\displaystyle\mathbf{y}_{t}=$
$\displaystyle g(\mathbf{h}_{t}),$ (2)

where the subscript $t$ $(1\leq t\leq T)$, which is a non-negative integer,
indexes time; the variable $\mathbf{h}_{t}\in\mathbb{R}^{d_{\mathbf{h}}}$ is
the system hidden state to be inferred at time $t$ and $\mathbf{h}_{0}$
represents the initial hidden state; the input signal (action)
$\mathbf{u}_{t}$ is known at all time points;
$\mathbf{y}_{t}\in\mathbb{R}^{d_{\mathbf{y}}}$ is the measurement
(observation) at time $t$; the transition function
$f:\mathbb{R}^{d_{\mathbf{h}}}\mapsto\mathbb{R}^{d_{\mathbf{h}}}$ evolves the
hidden state $\mathbf{h}_{t-1}$ to the next state $\mathbf{h}_{t}$;111When the
transition $f$ is a linear time-invariant mapping, we arrive at the linear
SSM. and the measurement function (also known as emission function)
$g:\mathbb{R}^{d_{\mathbf{h}}}\mapsto\mathbb{R}^{d_{\mathbf{y}}}$ maps the
hidden state to the associated output.

It is observed that when we are describing the functions $f$ and $g$ using
deep mappings, e.g., neural networks, it results in the so-called deep SSM
which takes the advances from two fields: the high flexibility from neural
networks and the high expressivity via hidden states. When future combining
SSM with generative neural networks, for example, VAE, we arrive at the
subclass of deep SSM, which is capable of performing probabilistic time series
forecasting and will be investigated below in our work. Note that the possible
perturbations, for example, modeling error and measurement error, have been
absorbed in (1) and (2).

##  3 Variational sequence model using augmented stochastic inputs

In order to incorporate stochasticity into the sequence model, we propose to
use additional random variable to augment the input space as
$\mathbf{u}_{t}\in\mathbb{R}^{d_{\mathbf{u}}}\mapsto[\mathbf{u}_{t},\mathbf{z}_{t}]\in\mathbb{R}^{d_{\mathbf{u}}+d_{\mathbf{z}}}$.
Consequently, rich stochastic sequence representation could be induced in the
augmented input space for improving both the model capability and the
quantification of uncertainty. The proposed model as well as the inference and
implementation details will be elaborated as below.

###  3.1 Model definition

Under the framework of latent variable model, the joint prior with additional
latent variables $\mathbf{z}_{1:T}$ is assumed to factorize over time as

$\displaystyle\begin{aligned}
p(\mathbf{y}_{1:T},\mathbf{z}_{1:T}|\mathbf{u}_{1:T})=&\prod_{t=1}^{T}p(\mathbf{y}_{t}|\mathbf{z}_{t},\mathbf{u}_{t})p(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t}).\end{aligned}$
(3)

In (3), the latent sequence $\mathbf{z}$ starts from the initial input state
$\mathbf{z}_{0}$ which could take for example the zero initialization for cold
starting. It is observed that the Gaussian state transition
$p(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t})$ $(1\leq t\leq T)$ (i)
evolves over time for encoding the sequence dependency; and (ii) it is
particularly conditioned on all the historical input signals
$\mathbf{u}_{1:t}$ until time $t$ and all the historical latent variables
$\mathbf{z}_{0:t-1}$ until time $t-1$, not only the information limited at
adjacent time points $t-1$ and $t$. The consideration of all the historical
information in $p(\mathbf{z}_{t}|.)$ is expected to help learn long sequence
patterns. Besides, the Gaussian likelihood
$p(\mathbf{y}_{t}|\mathbf{z}_{t},\mathbf{u}_{t})$ in (3) plays as the role of
decoder for mapping the encoded sequence inputs to the related output at time
$t$.222Note that the likelihood is not limited to Gaussian, it can adopt other
task-dependent distributions. This joint prior can be estimated by Monte Carlo
sampling conducted recursively over time.

For model inference, we employ the type-II maximum likelihood to tune model
parameters by maximizing the marginal likelihood

$\displaystyle p(\mathbf{y}_{1:T}|\mathbf{u}_{1:T})=\int
p(\mathbf{y}_{1:T},\mathbf{z}_{1:T}|\mathbf{u}_{1:T})d\mathbf{z}_{1:T}.$ (4)

When performing prediction, we require the posterior
$p(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})$, the distribution of
which however is unknown. We hence resort to variational inference by
introducing Gaussian variational posterior as an approximation. Specifically,
similar to the prior
$p(\mathbf{z}_{1:T}|\mathbf{u}_{1:T})=\prod_{t=1}^{T}p(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t})$
in (3), the variational posterior could be factorized over time as

$\displaystyle\begin{aligned}
q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})=&\prod_{t=1}^{T}q(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\mathbf{y}_{1:t}),\end{aligned}$
(5)

which describes the generation of augmented latent states given observations
and input signals up to current time, and each state transition
$q(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\mathbf{y}_{1:t})$ takes
the Gaussian form.

Again, it is worth noting that the additional encoding of all the previous
variables in (3) and (5) enables the learning of long term dependency. The
encoding method using generalized auto-regressive strategy will be detailed in
section 3.3, and its superiority over the usage of only adjacent variables
will be verified in the ablation study in section 4.3.

###  3.2 Derivation of ELBO

To simultaneously infer the interested marginal likelihood
$p(\mathbf{y}_{1:T}|\mathbf{u}_{1:T})$ and the variational posterior
$q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})$, we resort to
variational inference. To this end, we minimize the Kullback-Leibler (KL)
divergence between the variational posterior
$q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})$ and the exact
posterior $p(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})$, which is
formulated as

$\displaystyle\begin{aligned}
&\mathrm{KL}[q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})||p(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})]\\\
=&\mathbb{E}_{q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})}\left[\log\frac{q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})}{p(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})}\right]\\\
=&\log
p(\mathbf{y}_{1:T}|\mathbf{u}_{1:T})+\mathbb{E}_{q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})}\left[\log\frac{q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})}{p(\mathbf{z}_{1:T},\mathbf{y}_{1:T}|\mathbf{u}_{1:T})}\right].\end{aligned}$
(6)

Since the KL term (6) is always non-negative, we derive the following ELBO for
the log marginal likelihood $\log p(\mathbf{y}_{1:T}|\mathbf{u}_{1:T})$ as

$\displaystyle\begin{aligned} \log
p(\mathbf{y}_{1:T}|\mathbf{u}_{1:T})\geq\mathcal{L}=&\mathbb{E}_{q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})}\left[\log\frac{p(\mathbf{y}_{1:T},\mathbf{z}_{1:T}|\mathbf{u}_{1:T})}{q(\mathbf{z}_{1:T}|\mathbf{u}_{1:T},\mathbf{y}_{1:T})}\right].\end{aligned}$
(7)

By inserting the joint prior (3) and the variational posterior (5) into
$\mathcal{L}$, we arrive at the factorization over time as

$\displaystyle\begin{aligned}
\mathcal{L}=&\sum_{t=1}^{T}\mathbb{E}_{q(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\mathbf{y}_{1:t})}[\log
p(\mathbf{y}_{t}|\mathbf{z}_{t},\mathbf{u}_{t})]\\\
&-\mathrm{KL}[q(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\mathbf{y}_{1:t})||p(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t})],\end{aligned}$
(8)

where the latent prior usually takes the unit Gaussian
$p(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t})=\mathcal{N}(\mathbf{z}_{t}|\mathbf{0},\mathbf{I})$
for simplicity; and for the likelihoods and variational posteriors, all of
them take the Gaussian form amortized through neural networks, which will be
elaborated in section 3.3. Besides, given the Gaussian posterior
$q(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\mathbf{y}_{1:t})=\mathcal{N}(\mathbf{z}_{t}|\bm{\mu}_{\mathbf{z}_{t}},\bm{\Sigma}_{\mathbf{z}_{t}})$
with the diagonalized variance
$\bm{\Sigma}_{\mathbf{z}_{t}}=\mathrm{diag}[\bm{\nu}_{\mathbf{z}_{t}}]$, the
KL terms in the right-hand side of (8) have analytical expressions as

$\displaystyle\begin{aligned}
&\mathrm{KL}[q(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\mathbf{y}_{1:t})||p(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t})]\\\
=&\frac{1}{2}\left[\mathrm{Tr}(\bm{\Sigma}_{\mathbf{z}_{t}})+\bm{\mu}_{\mathbf{z}_{t}}^{\mathsf{T}}\bm{\mu}_{\mathbf{z}_{t}}-d_{\mathbf{z}}-\log|\bm{\Sigma}_{\mathbf{z}_{t}}|\right].\end{aligned}$
(9)

Finally, it is found that minimizing the KL divergence in (6) is equivalent to
maximizing the ELBO in (8), which (i) pushes the bound $\mathcal{L}$ towards
the marginal likelihood and (ii) infers the interested variational posterior
$q(\mathbf{z}_{t}|.)$.

The estimation of ELBO $\mathcal{L}$, especially the first likelihood
expectation in the right-hand side of (8), could be conducted through the re-
parameterization trick [36] due to the latent Gaussian random variables,
therefore resulting in a differentiable objective for optimization. For
instance, instead of directly sampling from the variational Gaussian posterior
$q(\mathbf{z}_{t}|.)$, we sample from another random variable
$\bm{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ without trainable
parameters, and then use it to obtain
$\hat{\mathbf{z}}_{t}=\bm{\mu}_{\mathbf{z}_{t}}+\sqrt{\bm{\nu}_{\mathbf{z}_{t}}}\odot\bm{\epsilon}$
where the symbol $\odot$ represents element-wise product.

Besides, directly evaluating the recurrent ELBO $\mathcal{L}$ unrolled over a
long time series requires high memory consumption due to the many historical
states. More efficiently, we employ the shingling technique proposed in [42]
to convert the long time series into several short time series chunks.
Specifically, given a long time series $\mathcal{D}$ with length of $T$, we
start from an arbitrary time index $t$ and obtain a short time series with
sequence length $W$ ($W<T$). Consequently, we could obtain $J=T-W+1$ short
time series chunks $\mathcal{D}_{c}=\\{\mathcal{D}_{c}^{j}\\}_{j=1}^{J}$, each
of which has the sequence length of $W$. This shingling data preprocessing
allows modeling many independent time series with varying length
simultaneously. Besides, the above time series split could naturally support
an efficient and unbiased estimation of $\mathcal{L}$ on a subset
$\mathcal{B}\subseteq\mathcal{D}_{c}$ of the short time series chunks with
$|\mathcal{B}|\ll J$, thus allowing using mini-batch optimizer, e.g., Adam
[43]. That is, an unbiased estimation of $\mathcal{L}$ on the subset
$\mathcal{B}$ can be expressed as

$\displaystyle\tilde{\mathcal{L}}=\frac{J}{|\mathcal{B}|}\sum_{\mathcal{D}_{c}^{j}\in\mathcal{B}}\mathcal{L}^{j},$
(10)

where $\mathcal{L}^{j}$ indicates the ELBO evaluated on the short time series
$\mathcal{D}_{c}^{j}$ via (8).

###  3.3 Amortized VRNNaug

After deriving the ELBO (8) for the proposed VRNNaug, the next is to elaborate
the practical implementations for constructing the variational posterior
$q(\mathbf{z}_{t}|.)$ and the likelihood $p(\mathbf{y}_{t}|.)$. The detailed
implementations of the proposed VRNNaug sequence model are depicted in Fig. 1.
It is observed that the whole model follows the encoder-decoder structure,
wherein the encoder embeds the temporal patterns into the stochastic, time-
aware manifold represented by $\mathbf{z}_{t}$, while the decoder generates
outputs from the time-aware manifold. Furthermore, we employ amortized
variational inference, i.e., amortizing the parameters over neural networks,
to build and train our proposed VRNNaug model efficiently. Particularly, (i) a
hybrid strategy is utilized to address the inconsistency of
$q(\mathbf{z}_{t}|.)$ between training and predicting, and (ii) a generalized
auto-regressive strategy is proposed in building $q(\mathbf{z}_{t}|.)$ in
order to summarize the long historical patterns, which will be elaborated in
what follows.

Figure 1: The inference and generative networks of the proposed VRNNaug model
for probabilistic time series forecasting. Note that the gray nodes represent
input signals and observations. Particularly, the transient gray node
$\dot{\mathbf{y}}_{t-1}$ indicates the mixture of observation
$\mathbf{y}_{t-1}$ and sample $\hat{\mathbf{y}}_{t-1}$. Finally, it is found
that in comparison to the inference network, the major difference in the
generative network is that it is using the sample $\hat{\mathbf{y}}_{t-1}$
since we are inaccessible to the true observation when predicting.

Inconsistency issue. Before proceeding to the amortized modeling, we first
discuss the inconsistency issue occurred in the proposed model. It is found
that the stochastic sequence model defined in (3) and (5) however may raise
inconsistency between training and predicting when we resort to the
variational posterior
$q(\mathbf{z}_{t}|\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\mathbf{y}_{1:t})$. That
is, different from the training phase, when we are going to predict at the
future time $T+1$, we have to use the prior
$p(\mathbf{z}_{T+1}|\mathbf{z}_{0:T},\mathbf{u}_{1:T+1})$ instead of the
variational posterior
$q(\mathbf{z}_{T+1}|\mathbf{z}_{0:T},\mathbf{u}_{1:T+1},\mathbf{y}_{1:T+1})$,
since there is no access to the further observation $\mathbf{y}_{T+1}$. The
gap raised by this inconsistency might deteriorate the quality of forecasting.

Inspired by [30], in order to alleviate the inconsistency issue, we could use
the lagged output $\mathbf{y}_{t-1}$ instead of $\mathbf{y}_{t}$ in the
conditions of the variational posterior $q(\mathbf{z}_{t}|.)$ with $t\leq T$
during the training phase. Besides, we sample $\hat{\mathbf{y}}_{t-1}\sim
p(\mathbf{y}_{t-1}|.)$ from the previous probabilistic output when predicting
at time $t$ with $T+1\leq t\leq T+F$. This auto-regressive strategy extracts
the known previous output or sample and feeds it as input for
$q(\mathbf{z}_{t}|.)$ at next time point, which is available for both training
and predicting. However, this strategy introduces another inconsistency: we
use the ground truth $\mathbf{y}_{t-1}$ during training, but the sample
$\hat{\mathbf{y}}_{t-1}$ during predicting. This discrepancy has been found to
raise error accumulation, thus deteriorating the RNN performance on natural
language processing (NLP) tasks [44].

Therefore, inspired by [45], we decide to use a simple trick to alleviate the
discrepancy by taking the following lagged hybrid

$\displaystyle\dot{\mathbf{y}}_{t-1}=\frac{1}{2}(\mathbf{y}_{t-1}+\hat{\mathbf{y}}_{t-1}),\quad
1\leq t\leq T,$ (11)

and feeding it as input at next time point for the posterior
$q(\mathbf{z}_{t}|.)$ during training, see Fig. 1(a). The ground truth
$\mathbf{y}_{t-1}$ in (11) plays as the teaching force to accelerate model
convergence, while the usage of predictive sample $\hat{\mathbf{y}}_{t-1}$
plays as the regularizer to reduce the gap between training and predicting.
Besides, akin to $\mathbf{z}_{0}$, the initials $\mathbf{y}_{0}$ and
$\hat{\mathbf{y}}_{0}$ in (11) could simply take zero initialization for cold
starting, thus resulting in $\dot{\mathbf{y}}_{0}=\mathbf{0}$. Note that when
the learned output distributions are close to the true distributions of
observations in the later stage of training, we have
$\dot{\mathbf{y}}_{t-1}\approx\mathbf{y}_{t-1}$.

NN-assisted modeling. To perform amortized variational inference, we attempt
to parameterize the Gaussian variational posterior $q(\mathbf{z}_{t}|.)$ and
the Gaussian likelihood $p(\mathbf{y}_{t}|.)$ in ELBO (8) through neural
networks. Firstly, the mean $\bm{\mu}_{\mathbf{z}_{t}}$ and variance
$\bm{\nu}_{\mathbf{z}_{t}}$ of $q(\mathbf{z}_{t}|.)$ depend on the historical
data $\mathbf{z}_{0:t-1}$, $\mathbf{u}_{1:t}$ and $\dot{\mathbf{y}}_{0:t-1}$,
and are parameterized through a multi-layer perceptron (MLP) as

$\displaystyle\begin{aligned}
\left[\bm{\mu}_{\mathbf{z}_{t}},\log\bm{\nu}_{\mathbf{z}_{t}}\right]=&f^{z}(\mathbf{z}_{0:t-1},\mathbf{u}_{1:t},\dot{\mathbf{y}}_{0:t-1})=\mathrm{MLP}^{z}(\bar{\mathbf{z}}_{t-1},\bar{\mathbf{u}}_{t},\bar{\mathbf{y}}_{t-1})\end{aligned}$
(12)

In (12), the recurrent transformations
$\bar{\mathbf{z}}_{t-1}\in\mathbb{R}^{d_{\mathbf{z}}}$,
$\bar{\mathbf{u}}_{t}\in\mathbb{R}^{d_{\mathbf{u}}}$ and
$\bar{\mathbf{y}}_{t-1}\in\mathbb{R}^{d_{\mathbf{y}}}$ respectively employ the
generalized auto-regressive strategy to evolve over all the historical states
$\mathbf{z}_{0:t-1}$, $\mathbf{u}_{1:t}$ and $\dot{\mathbf{y}}_{0:t-1}$ as

$\displaystyle\bar{\mathbf{z}}_{t-1}=$
$\displaystyle\phi^{z}(\mathbf{z}_{0:t-1}),$ (13)
$\displaystyle\bar{\mathbf{u}}_{t}=$
$\displaystyle\phi^{u}(\mathbf{u}_{1:t}),$ (14)
$\displaystyle\bar{\mathbf{y}}_{t-1}=$
$\displaystyle\phi^{y}(\dot{\mathbf{y}}_{0:t-1}),$ (15)

in order to extract long temporal patterns. Specifically, the recurrent
mappings $\phi^{z}(.)$, $\phi^{u}(.)$ and $\phi^{y}(.)$ shown in Fig. 1 are
achieved by building upon the outputs of for example RNN models as

$\displaystyle\mathbf{h}^{z}_{t-1}$
$\displaystyle=\mathrm{RNN}^{z}(\mathbf{h}^{z}_{t-2},\mathbf{z}_{t-1}),\quad\bar{\mathbf{z}}_{t-1}=\mathrm{MLP}^{\bar{z}}(\mathbf{h}^{z}_{t-1}),$
(16) $\displaystyle\mathbf{h}^{u}_{t}$
$\displaystyle=\mathrm{RNN}^{u}(\mathbf{h}^{u}_{t-1},\mathbf{u}_{t}),\quad\bar{\mathbf{u}}_{t}=\mathrm{MLP}^{\bar{u}}(\mathbf{h}^{u}_{t}),$
(17) $\displaystyle\mathbf{h}^{y}_{t-1}$
$\displaystyle=\mathrm{RNN}^{y}(\mathbf{h}^{y}_{t-2},\dot{\mathbf{y}}_{t-1}),\quad\bar{\mathbf{y}}_{t-1}=\mathrm{MLP}^{\bar{y}}(\mathbf{h}^{y}_{t-1}).$
(18)

The above RNN cells, for instance, GRU [24] and LSTM [25], naturally encode
the sequence characteristics. It is notable that the recurrent mappings
$\phi^{z}(.)$, $\phi^{u}(.)$ and $\phi^{y}(.)$ are not limited to RNN models.
The recently developed sequence extractors, for example, temporal
convolutional networks (TCN) [46] and transformer [47], are also available
here.

Besides, it is worth noting that the recurrent transformation
$\bar{\mathbf{y}}_{t-1}$ in (15) and the hidden state $\mathbf{h}^{y}_{t-1}$
in (18) solely hold for training, see the inference network in Fig. 1(a). When
predicting, we have
$\bar{\mathbf{y}}_{t-1}=\phi^{y}(\hat{\mathbf{y}}_{0:t-1})$ that accepts the
samples from previous output distributions, with the hidden state now
formulated as

$\displaystyle\mathbf{h}^{y}_{t-1}=\mathrm{RNN}^{y}(\mathbf{h}^{y}_{t-2},\hat{\mathbf{y}}_{t-1}),$
(19)

see the generative network in Fig. 1(b).

Secondly, as for the Gaussian likelihood (decoder) in the ELBO (8), it writes
as

$\displaystyle
p(\mathbf{y}_{t}|\mathbf{z}_{t},\mathbf{u}_{t})=\mathcal{N}(\mathbf{y}_{t}|\bm{\mu}_{\mathbf{y}_{t}},\bm{\Sigma}_{\mathbf{y}_{t}}),$
(20)

where the diagonalized variance
$\bm{\Sigma}_{\mathbf{y}_{t}}=\mathrm{diag}[\bm{\nu}_{\mathbf{y}_{t}}]$. Since
$\mathbf{z}_{t}$ is mapped from the recurrent inputs $\bar{\mathbf{z}}_{t-1}$,
$\bar{\mathbf{u}}_{t}$ and $\bar{\mathbf{y}}_{t-1}$ via $\phi^{z}(.)$,
$\phi^{u}(.)$ and $\phi^{y}(.)$, we follow the architecture of Dense-Net [48]
by extracting and concatenating all the previous inputs in Fig. 1 in order to
parameterize the Gaussian parameters through MLP as

$\displaystyle[\bm{\mu}_{\mathbf{y}_{t}},\log\bm{\nu}_{\mathbf{y}_{t}}]=\mathrm{MLP}^{y}(\mathbf{z}_{t},\mathbf{u}_{t},\bar{\mathbf{z}}_{t-1},\bar{\mathbf{u}}_{t},\bar{\mathbf{y}}_{t-1}).$
(21)

This sort of dense connection between encoder and decoder eases model training
because of the block-cross feature sharing. Finally, we again highlight that
the likelihood is not limited to Gaussian, it can adopt other task-dependent
distributions. For instance, we could choose the negative-binomial likelihood
for positive count data, and the Bernoulli likelihood for binary data [30].

###  3.4 Prediction

As shown in Fig. 1(b), given the trained deep sequence model $\mathcal{M}$, in
order to perform recursive prediction over the future time $t\in[T+1,T+F]$
with forecast horizon $F$, we first recursively sample $\hat{\mathbf{z}}_{t}$
from the state transition $q(\mathbf{z}_{t}|.)$. Note that the recurrent
transformation $\bar{\mathbf{y}}_{t-1}$ in (12) now is unrolled over the
previous output samples $\hat{\mathbf{y}}_{0:t-1}$ via (19) since the true
output observations are inaccessible when predicting. Thereafter, we pass the
latent sample through the decoder $p(\mathbf{y}_{t}|.)$ to output the
prediction sample $\hat{\mathbf{y}}_{t}$ at current time point. It is notable
that although the model employs Gaussian variables, the final predictive
distribution is not Gaussian. Hence, we could use Monte Carlo sampling to
output $K$ samples for empirically quantifying the implicit predictive
distribution.

Finally, pseudo implementations of the proposed VRNNaug model for training and
predicting are provided in Algorithms 1 and 2, respectively. Note that the
training process of VRNNaug in Algorithm 1 is conducted on the whole time
series with length $T$. One could employ the aforementioned shingling
technique to make data preprocessing and then use the mini-batch ELBO (10) to
train the model more efficiently. Besides, for line 4 in the prediction
process of Algorithm 2, when $t=T+1$, we sample from previous state and output
at $t=T$, which requires additional memory consumption; alternatively, one
could start from cold by simply adopting the zero initialization. Finally, it
is observed that the major difference of the for-loop in Algorithms 1 and 2 is
the estimation of the recurrent transformation $\bar{\mathbf{y}}_{t-1}$. When
training, we use the mixture of observation $\mathbf{y}_{t-1}$ and prediction
$\hat{\mathbf{y}}_{t-1}$ in (11) and (18); contrarily, when predicting,
without the access to observations, we only use the prediction
$\hat{\mathbf{y}}_{t-1}$ in (19).

Algorithm 1 Training of the proposed VRNNaug

1: Inputs: Training time series data
$\mathcal{D}=\\{\mathbf{u}_{1:T}\in\mathbb{R}^{T\times
d_{\mathbf{u}}},\mathbf{y}_{1:T}\in\mathbb{R}^{T\times d_{\mathbf{y}}}\\}$

2: Output: Sequence model $\mathcal{M}$ for future forecasting

3: repeat

4: for $t=1:T$ do

5: Perform recurrent transformations
$\bar{\mathbf{z}}_{t-1},\bar{\mathbf{u}}_{t},\bar{\mathbf{y}}_{t-1}$ via
(16)-(18)

6: Estimate mean $\bm{\mu}_{\mathbf{z}_{t}}$ and variance
$\bm{\nu}_{\mathbf{z}_{t}}$ of $q(\mathbf{z}_{t}|.)$ via (12)

7: Obtain mean $\bm{\mu}_{\mathbf{y}_{t}}$ and variance
$\bm{\nu}_{\mathbf{y}_{t}}$ of $p(\mathbf{y}_{t}|.)$ via (21)

8: end for

9: Calculate ELBO via (8)

10: Update model parameters via stochastic optimizer

11: until Optimizer has been terminated

Algorithm 2 Prediction of the proposed VRNNaug

1: Inputs: Sequence model $\mathcal{M}$, forecast horizon $F$, and input
signals $\\{\mathbf{u}_{t}\\}_{t=T+1}^{T+F}$

2: Output: $K$ prediction samples

3: for $t=T+1:T+F$ do

4: Perform recurrent transformations
$\bar{\mathbf{z}}_{t-1},\bar{\mathbf{u}}_{t},\bar{\mathbf{y}}_{t-1}$ via (16),
(17) and (19)

5: Estimate mean $\bm{\mu}_{\mathbf{z}_{t}}$ and variance
$\bm{\nu}_{\mathbf{z}_{t}}$ of $q(\mathbf{z}_{t}|.)$ via (12)

6: Obtain mean $\bm{\mu}_{\mathbf{y}_{t}}$ and variance
$\bm{\nu}_{\mathbf{y}_{t}}$ of $p(\mathbf{y}_{t}|.)$ via (21)

7: Sample from $p(\mathbf{y}_{t}|.)$ to produce $K$ prediction samples at time
$t$

8: end for

###  3.5 Discussions

So far, we have described the proposed deep sequence model with known input
signal $\mathbf{u}_{t}$ at any time point, which includes additional covariate
that may affect the output. Note that our model is mainly designed for free
forecasting, i.e., multi-step ahead prediction. If one wants to conduct one-
step ahead prediction, we could simply feed the ground truth
$\mathbf{y}_{t-1}$ into both the inference and generative networks in Fig. 1.
For time series without additional control signals, one could simply use the
time features, for example, hour-of-day and day-of-week, as the additional
$\mathbf{u}_{t}$ for modeling.

Finally, the relation and difference of our model to existing literature
should be highlighted. The most closely related works come from [37, 40, 31,
41, 30]. In comparison to the similar VAE-type probabilistic sequence models,
the proposed VRNNaug brings several improvements: (i) it adopts hybrid output
in (11) during training for greatly alleviating the inconsistency issue; and
(ii) it employs the generalized auto-regressive strategy in (13)-(15) to form
the recurrent inputs for better capturing long temporal dependency. The
benefits brought by these improvements will be verified by the comparative
study in section 4.2 and the ablation study in section 4.3.

##  4 Numerical experiments

This section first showcases the methodological characteristics of proposed
VRNNaug model on two toy cases, followed by the extensive comparison against
existing RNN assisted stochastic sequence models on eight system
identification benchmarks and a real-world centrifugal compressor sensor data.

We implement our model using PyTorch on a Linux workstation with TITAN RTX
GPU. The detailed experimental configurations are provided in Appendix C. As
for the probabilistic forecasting provided by our model, we quantify the
quality of predictive distribution by the $p$50- and $p$90-quantile losses
[32, 30, 49]. Generally, given the $i$-th time series output $\mathbf{y}^{i}$
($1\leq i\leq d_{\mathbf{y}}$) and its corresponding $\rho$-quantile
prediction $\hat{\mathbf{y}}^{i}_{\rho}$, the $\rho$-quantile loss is
expressed as

$\displaystyle\mathrm{QL}_{\rho}(\mathbf{y}^{i},\hat{\mathbf{y}}^{i}_{\rho})=2\times\frac{\sum_{t}P_{\rho}(y^{i}_{t},\hat{y}^{i}_{\rho,t})}{\sum_{t}|y^{i}_{t}|},\quad
P_{\rho}=\left\\{\begin{matrix}&\rho(y-\hat{y}_{\rho}),\quad&y>\hat{y}_{\rho},\\\
&(1-\rho)(\hat{y}_{\rho}-y),\quad&\mathrm{otherwise},\\\ \end{matrix}\right.$
(22)

where the parameter $\rho\in[0,1]$. It is found that the $p50$ criterion using
$\rho=0.5$ quantifies the relative error of prediction since
$\mathrm{QL}_{0.5}(\mathbf{y}^{i},\hat{\mathbf{y}}^{i}_{\rho})=\sum_{t}|y^{i}_{t}-\hat{y}^{i}_{t}|/\sum_{t}|y^{i}_{t}|$
wherein $\hat{y}^{i}_{t}$ now is the medium of prediction samples at time $t$.
As for the $p90$ criterion, it could represent the coverage of predictive
distribution since if the observation falls below the 90%-quantile prediction,
we have a small $P_{\rho}=0.1(\hat{y}_{\rho}-y)$; otherwise, a large
$P_{\rho}=0.9(y-\hat{y}_{\rho})$ should be used. Moreover, given desirable
predictions, the predictive distribution with uncertainty (i.e., a large
variance) might not be favored by the $p90$ criterion because of the large
difference $|\hat{y}_{\rho}-y|$ in $P_{\rho}$. Overall, lower $p50$ and $p90$
values indicate better quality of predictive distribution. Hence, in the
following experiments, we use $p50=\mathrm{QL}_{0.5}(.,.)$ to quantify the
accuracy of prediction medium, and $p90=\mathrm{QL}_{0.9}(.,.)$ to represent
the coverage of predictive distribution.

###  4.1 Toy cases

We below study the proposed VRNNaug model on two toy cases (the linear
Gaussian system and the heteroscedastic motorcycle case) with different
features.

Linear Gaussian system. We first consider a one-dimensional linear Gaussian
system [41] expressed as

$\displaystyle\mathbf{h}_{t+1}=$ $\displaystyle\left[\begin{matrix}0.7&0.8\\\
0&0.1\\\ \end{matrix}\right]\mathbf{h}_{t}+\left[\begin{matrix}-1\\\ 0.1\\\
\end{matrix}\right]u_{t}+\bm{\epsilon}_{h},$ (23) $\displaystyle y_{t}=$
$\displaystyle\left[\begin{matrix}1&0\\\
\end{matrix}\right]\mathbf{h}_{t}+\epsilon,$ (24)

where $\mathbf{h}\in\mathbb{R}^{2}$ is a two-dimensional hidden state,
$\bm{\epsilon}_{h}\sim\mathcal{N}(\mathbf{0},0.5\mathbf{I})$ is the process
noise, and $\epsilon\sim\mathcal{N}(0,1)$ is the measurement noise. We use
this dynamic model to sequentially generate 2000 training points, 2000
validation points and 5000 test points. Note that when generating training and
validation signals, an excitation input signal $u_{t}\in[-2.5,2.5]$ polluted
with uniform random noise is employed for the linear system. Differently, the
input signal for the test data writes as

$\displaystyle
u_{t}=\sin\left(\frac{2t\pi}{10}\right)+\sin\left(\frac{2t\pi}{25}\right).$

Figure 2: The predictive distribution of VRNNaug on the linear Gaussian
system. The red curve represents 50-quantile prediction, while the associated
shaded region indicates the interval bounded by 95-quantile and 5-quantile
predictions. The blue curve is the true observation with the shaded region
representing 95% confidence interval.

The predictive distributions of the proposed VRNNaug together with the true
observations are depicted in Fig. 2. It is observed that for this linear
system, the flexible and nonlinear VRNNaug model successfully captures the
underlying linear dynamics. This linear case has also been used in [41] to
test another deep SSM model, named STORN [37], and the gray-box model SSEST
[50] with two latent states. Though successfully capturing the dynamics, they
conservatively estimate the uncertainty of prediction. Our model however is
found to well estimate the uncertainty in Fig. 2.

Heteroscedastic motorcycle data. We next study a more complicated motorcycle
time series dataset [51] composed of 133 points with time input and
acceleration output. The motorcycle data takes from an experiment on the
efficacy of crash helmets, and yields multiple acceleration observations at
some time points. In comparison to the previous linear Gaussian system with
homoscedastic noise, this nonlinear dataset contains more challenging time-
varying noise.

Figure 3: The predictive distribution of VRNNaug on the heteroscedastic
motorcycle dataset. The red curve represents 50-quantile prediction, while the
associated shaded region indicates the interval bounded by 95-quantile and
5-quantile predictions. The circles represent 133 data points.

We seek to verify the capability of VRNNaug for capturing nonlinear dynamic
patterns with heteroscedastic noise on this case. We here use all the 133 data
points for training, validation and testing, and extract the time feature as
input signal for our VRNNaug model. As depicted in Fig. 3, the proposed
VRNNaug captures the non-linear dynamics of underlying acceleration curve, for
example, the constant behavior before $t=15$ ms and the step behavior around
$t=15$ ms. More interestingly, it well describes the time-varying noise, for
example, the pretty small perturbation before $t=15$ ms and the large and
time-varying noise thereafter.

###  4.2 System identification benchmarks

This section verifies the performance of proposed deep sequence model on eight
popular system identification benchmarks from dynamics systems like hydraulic
actuator, electric motors, hair dryer, furnace, tank, cascaded tanks system
[52], F-16 aircraft ground vibration test [53], and wiener-hammerstein process
noise system [54].333The first five datasets are available at
https://homes.esat.kuleuven.be/~smc/daisy/daisydata.html and the last three
datasets are collected at http://nonlinearbenchmark.org/index.html. Table 1
summarizes the sequence length and dimensions of the considered time series
datasets. It is found that all these datasets except f16gvt have a single
input signal, and the particular tank dataset has two outputs. The goal of
sequence modeling is to infer the underlying high-dimensional latent state and
the relevant nonlinear dynamics from data, and use them for future
forecasting.

Table 1: Characteristics of the system identification benchmarks. dataset | $T$ | $d_{\mathbf{u}}$ | $d_{\mathbf{y}}$
---|---|---|---
actuator | 1024 | 1 | 1
cascaded tank (ctank) | 2028 | 1 | 1
drive | 500 | 1 | 1
dryer | 1000 | 1 | 1
f16gvt | 32768 | 2 | 1
gas furnace (furnace) | 296 | 1 | 1
tank | 2500 | 1 | 2
wiener hammerstein (hammerstein) | 32768 | 1 | 1

The performance of VRNNaug is compared against state-of-the-art RNN assisted
deep stochastic SSM models, including (i) the linear SSM model with parameters
learned through neural networks (Deep-LSSM) [32]; (ii) the VAE-RNN model with
the encoder built upon RNN [40]; (iii) the stochastic RNN (STORN) following
the encoder-decoder structure [37]; (iv) the deep auto-regressive model (Deep-
AR) [30]; and finally, (v) the deep SSM for probabilistic forecasting, donated
as DSSMF [31].444The DSSMF in [31] additionally adopts an automatic relevance
determination network to identify the importance of dimensions of the input
signal $\mathbf{u}_{t}$, which is not implemented in our comparison. Tables 2
and 3 report the $p50$ and $p90$ results of various deep stochastic SSM models
over ten runs, respectively. Note that the best and second-best results on the
eight benchmarks are marked as gray and light gray, respectively.

Table 2: The $p$50 results of various deep stochastic SSM models on the eight system identification benchmarks. The best and second-best results on each dataset are marked as gray and light gray, respectively. Datasets | Deep-LSSM | VAE-RNN | STORN | Deep-AR | DSSMF | VRNNaug
---|---|---|---|---|---|---
actuator | 0.3866±0.0225 | 0.3529±0.0263 | 0.3503±0.0104 | 0.3839±0.0407 | 0.3737±0.0536 | 0.3241±0.0265
ctank | 0.2777±0.0469 | 0.1558±0.0174 | 0.1065±0.0092 | 0.0965±0.0082 | 0.1028±0.0333 | 0.0715±0.0192
drive | 0.3157±0.0359 | 0.4925±0.0057 | 0.3792±0.0599 | 0.5083±0.0072 | 0.5056±0.0119 | 0.2098±0.0186
dryer | 0.0201±0.0011 | 0.0203±0.0013 | 0.0162±0.0010 | 0.0179±0.0016 | 0.0153±0.0025 | 0.0158±0.0017
f16gvt | 0.7108±0.1586 | 0.0723±0.0072 | 0.0735±0.0044 | 0.2973±0.0111 | 0.3021±0.0124 | 0.1540±0.0140
furnace | 0.0238±0.0021 | 0.0436±0.0147 | 0.0225±0.0006 | 0.0231±0.0005 | 0.0242±0.0009 | 0.0226±0.0008
tank ($y_{1}$) | 0.1005±0.0203 | 0.1885±0.0155 | 0.0785±0.0212 | 0.0539±0.0169 | 0.0529±0.0072 | 0.0591±0.0154
tank ($y_{2}$) | 0.0597±0.0229 | 0.0615±0.0104 | 0.0559±0.0083 | 0.0414±0.0085 | 0.0437±0.0042 | 0.0429±0.0063
hammerstein | NA | 0.9400±0.1816 | 0.7364±0.0526 | 0.4676±0.0582 | 0.5088±0.0972 | 0.2261±0.0777

Table 3: The $p$90 results of various deep stochastic SSM models on the eight system identification benchmarks. The best and second-best results on each dataset are marked as gray and light gray, respectively. Datasets | Deep-LSSM | VAE-RNN | STORN | Deep-AR | DSSMF | VRNNaug
---|---|---|---|---|---|---
actuator | 0.2059±0.0226 | 0.2302±0.0280 | 0.1894±0.0124 | 0.1965±0.0214 | 0.1925±0.0172 | 0.1776±0.0158
ctank | 0.2209±0.0826 | 0.0979±0.0111 | 0.0656±0.0093 | 0.0758±0.0179 | 0.0611±0.0192 | 0.0358±0.0138
drive | 0.1610±0.0164 | 0.1966±0.0039 | 0.1737±0.0276 | 0.2173±0.0101 | 0.2245±0.0118 | 0.1029±0.0058
dryer | 0.0079±0.0010 | 0.0093±0.0015 | 0.0060±0.0009 | 0.0082±0.0012 | 0.0062±0.0015 | 0.0067±0.0012
f16gvt | 0.3180±0.0664 | 0.0342±0.0026 | 0.0353±0.0018 | 0.1310±0.0053 | 0.1356±0.0059 | 0.0733±0.0065
furnace | 0.0210±0.0009 | 0.0176±0.0013 | 0.0221±0.0007 | 0.0229±0.0012 | 0.0222±0.0010 | 0.0197±0.0009
tank ($y_{1}$) | 0.0386±0.0089 | 0.1219±0.0127 | 0.0605±0.0182 | 0.0592±0.0281 | 0.0566±0.0115 | 0.0318±0.0159
tank ($y_{2}$) | 0.0287±0.0160 | 0.0286±0.0084 | 0.0344±0.0109 | 0.0463±0.0136 | 0.0504±0.0072 | 0.0285±0.0091
hammerstein | NA | 0.3992±0.0603 | 0.3689±0.0288 | 0.2024±0.0313 | 0.2122±0.0339 | 0.1443±0.0212

It is first observed that due to the linear SSM framework, the Deep-LSSM is
hard to learn the nonlinear dynamic behaviors from sequence data, even though
it attempts to learn the local, time-varying coefficients. For instance,
though the Deep-LSSM forecasts the rough profile on the ctank dataset, its
prediction however leaves far away from the ground truth in comparison to the
others, which is indicated by the poor $p50$ results in Table 2. Furthermore,
it fails on the large and complicated hammerstein dataset. Therefore, in
comparison to other competitors, the Deep-LSSM showcases the worst performance
in general in terms of both $p50$ and $p90$. Despite the linear model
structure, the undesirable performance of Deep-LSSM may be attributed to the
type-I maximum likelihood framework wherein the latent variables are not
integrated out.

As for the VAE-type deep SSM models, the more complicated STORN is found to
outperform the simple VAE-RNN on most cases in terms of both $p50$ and $p90$.
This is attributed to the description of the recurrency in previous outputs
$\mathbf{d}_{t}=\mathrm{RNN}(\mathbf{y}_{1:t})$ when determining the posterior
$q(\mathbf{z}_{t}|.)$. But as has been discussed, the usage of current output
$\mathbf{y}_{t}$ at time $t$ raises the inconsistency issue: when forecasting
time series through generative networks, we have no access to current output
and then have to sample $\hat{\mathbf{z}}_{t}$ from prior. This thus makes
STORN’s $p50$ results in Table 3 slightly worse than that of Deep-AR and DSSMF
which use the auto-regressive strategy to alleviate the issue. Besides, the
sampling from prior when generating time series data makes VAE-RNN and STORN
tend to provide relatively large uncertainty estimation, which will be
illustrated in Fig. 4 as well as the compressor case in next section.

As for the last three deep SSM models, all of them employ sort of auto-
regressive strategy to avoid the inconsistency issue. In comparison to Deep-AR
and DSSMF, the proposed VRNNaug adopts (i) the hybrid output to further
alleviate the gap between training and predicting; and (ii) the generalized
auto-regressive strategy for extracting long term historical patterns. As a
result, it performs the best on these system identification benchmarks in
terms of both $p50$ and $p90$.555The VRNNaug forecasting in future time period
on the eight system identification benchmarks is respectively illustrated in
Appendix D. Besides, it is found that akin to Deep-LSSM, the Deep-AR also
trains the model in the simple type-I maximum-likelihood framework, thus
making it perform slightly worse than the competitor DSSMF derived through
variational inference, especially in terms of the $p90$ criterion.

Finally, to fully represent the quality of prediction distribution of the
above probabilistic sequence models, Fig. 4 adopts the empirical coverage
probability (ECP) criterion used in [33] to compare their performance on the
actuator and drive datasets. It is found that the proposed VRNNaug achieves
good empirical coverage on the two datasets. As for VAE-RNN, its ECP is higher
than the theoretical coverage on the actuator dataset, which implies
conservative uncertainty that has also been identified by the worst $p90$
value in Table 3. As for the challenging drive dataset, since the predictions
of all the models except VRNNaug are usually far away from the
observations,666For example, the VAE-RNN fails on the drive dataset by
providing almost the constant predictions. indicated by the poor $p50$ results
in Table 2, their ECP results in Fig. 4(b) are not good. The exception of good
90% ECP of VAE-RNN in Fig. 4(b) is solely attributed to the large uncertainty.

Figure 4: Coverage plots of predictive distributions of various probabilistic
sequence models on the (a) actuator and (b) drive datasets.

###  4.3 Ablation study

With the good performance reported in Tables 2 and 3, we here perform ablation
study of the proposed VRNNaug model on the drive dataset. The ablation study
includes two additional VRNNaug variants in order to reveal individual effects
of components: (i) the VRNNaug-v1 which resembles Deep-AR by using only the
auto-regressive strategy via $\mathbf{y}_{t-1}$; and (ii) the VRNNaug-v2 which
uses the generalized auto-regressive strategy in (16) and (17), but simply
uses $\mathbf{y}_{t-1}$ instead of the hybrid $\dot{\mathbf{y}}_{t-1}$ in (18)
to determine $\bar{\mathbf{y}}_{t-1}$. Differently, the complete VRNNaug model
adopts both the previous hybrid outputs and the generalized auto-regressive
strategy to alleviate the inconsistency issue and improve the model
capability.

Figure 5: Ablation study of the proposed VRNNaug model on the drive dataset.

Fig. 5 depicts the comparative results of three VRNNaug models in terms of
both the $p50$ and $p90$ criteria on the drive dataset, and has the following
findings.

  * 1.

By solely using the auto-regressive strategy akin to Deep-AR and DSSMF, the
VRNNaug-v1 is incapable of capturing the underlying nonlinear dynamical
behavior on the challenging drive dataset. Contrarily, the generalized auto-
regressive strategy summarizes all the long term historical patterns of not
only $\mathbf{y}$ but also $\mathbf{z}$ and $\mathbf{u}$. Consequently, both
VRNNaug-v2 and VRNNaug significantly outperform the basic VRNNaug-v1;

  * 2.

The hybrid output $\dot{\mathbf{y}}_{t-1}$ used in training helps further fill
the gap between inference and generation, since now we always consider the
prediction samples in the two phases. Hence, it brings considerable
improvements for VRNNaug over VRNNaug-v2.

###  4.4 Axis displacement forecasting of centrifugal compressor

This section verifies the performance of the proposed VRNNaug on an
11-variable time series data collected from a real-world centrifugal
compressor with the rotating speed of 5556 RPM [55]. As illustrated in Fig. 6,
this time series data was recorded every one hour from March to October, 2013,
by placing sensors at the inlet, the inlet guide vanes, the exhaust, the
bearing, and the axis of centrifugal compressor, thus resulting in 3710 raw
data points for each variable.

Figure 6: Illustration of sensors layout for collecting the compressor time
series data.

As for this real-world time series task, we perform data preprocessing before
applying deep stochastic SSM models on it for forecasting. These 11-variable
time series data is collected manually from sensors, thus may containing human
errors, outliers, noise and missing values. Hence, we first perform outlier
analysis to identify and remove the outliers in each raw time series variable
via the box-plot method, thus resulting in 3524 data points. Secondly, in
order to eliminate the effect of magnitudes within multivariate time series,
we normalize all the variables into the same range $[0,1]$. Thirdly, we
conduct discrete wavelet packet transform (DWPT) analysis by the Daubechies
wavelet of eight-order [56] as well as the Bayesian thresholding technique in
order to denoise our time series data.

By processing and denoising the raw data, we finally obtain the compressor
time series containing 3524 data points. Following the experimental
configurations elaborated in Appendix C, our aim is to predict the axis
displacement at the last 10% time points, given the remaining ten variables as
covariates.

Table 4: The $p50$ and $p90$ results of various deep stochastic SSM models on the real-world compressor dataset. Note that the best and second-best results on each dataset are marked as gray and light gray, respectively. Method | $p50$ | $p90$
---|---|---
Deep-LSSM | 0.0597±0.0119 | 0.0428±0.0106
VAE-RNN | 0.0563±0.0093 | 0.0259±0.0040
STORN | 0.0540±0.0067 | 0.0270±0.0055
Deep-AR | 0.0674±0.0294 | 0.0531±0.0371
DSSMF | 0.0600±0.0188 | 0.0581±0.0351
VRNNaug | 0.0488±0.0050 | 0.0255±0.0068

Table 4 reports the forecasting results of VRNNaug against that of other deep
stochastic SSM models on the real-world compressor dataset. The superiority of
VRNNaug over other competitors has again been observed in terms of both the
$p50$ and $p90$ criteria for forecasting the future axis displacement.
Moreover, Fig. 7 illustrates the forecasting of these deep stochastic SSM
models for the future axis displacement. It is observed that the proposed
VRNNaug provides more accurate predictions. Besides, the uncertainty of Deep-
AR and DSSMF is relatively small, thus indicated by the poor $p90$ results in
Table 4. Contrarily, the uncertainty of VAE-RNN and STORN is relatively large,
which may be attributed to the sampling from prior $p(\mathbf{z}_{t}|.)$
rather than the informative posterior $q(\mathbf{z}_{t}|.)$ in the generative
network.777Note that in comparison to the small uncertainty of Deep-AR and
DSSMF, this conservative uncertainty is favored by the $p90$ criterion since
it might cover most observations on the compressor case. The proposed VRNNaug
achieves a trade-off between the two extreme cases, thus providing desirable
estimation of uncertainty, which is revealed by the best $p90$ results.

Figure 7: The future axis displacement of centrifugal compressor forecasted by
various deep SSM models. The red curve represents 50-quantile prediction,
while the associated shaded region indicates the interval bounded by
95-quantile and 5-quantile predictions.

##  5 Conclusions

In this paper, we present a variational RNN model on the augmented recurrent
space for probabilistic time series forecasting. When using variational
inference to derive the ELBO as objective for model training, we figure out
the issue of inconsistency that would deteriorate the model performance.
Hence, we propose (i) feeding the hybrid output as input at next time step,
and (ii) presenting the generalized auto-regressive strategy in order to
better extract temporal patterns. Through extensive experiments on various
time series tasks against existing deep stochastic SSM models, we have
showcased the superiority of proposed VRNNaug model by quantifying predictive
distribution. Future work would attempt to extend the proposed model to
challenging scenarios in dynamic systems with for example missing outputs,
irregular time points, and knowledge transfer across multivariate time series.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China
(52005074), and the Fundamental Research Funds for the Central Universities
(DUT19RC(3)070). Besides, it was partially supported by the Research and
Innovation in Science and Technology Major Project of Liaoning Province
(2019JH1-10100024), and the MIIT Marine Welfare Project (Z135060009002).

##  Appendix A Acronyms

$\begin{array}[]{ll}\mbox{ARIMA}&\mbox{Auto-Regressive Integrated Moving
Average}\\\ \mbox{CFP}&\mbox{Empirical Coverage Probability}\\\ \mbox{Deep-
AR}&\mbox{Deep Auto-regressive
Model~{}\cite[cite]{[\@@bibref{Number}{salinas2020deepar}{}{}]}}\\\
\mbox{Deep-LSSM}&\mbox{Deep Linear State Space
Model~{}\cite[cite]{[\@@bibref{Number}{rangapuram2018deep}{}{}]}}\\\
\mbox{DSSMF}&\mbox{Deep State Space Model for Probabilistic
Forecasting~{}\cite[cite]{[\@@bibref{Number}{li2021learning}{}{}]}}\\\
\mbox{DWPT}&\mbox{Discrete Wavelet Packet Transform}\\\
\mbox{ELBO}&\mbox{Evidence Lower Bound}\\\ \mbox{GRU}&\mbox{Gated Recurrent
Unit}\\\ \mbox{KL}&\mbox{Kullback-Leibler}\\\ \mbox{LSTM}&\mbox{Long Short
Term Memory}\\\ \mbox{MLP}&\mbox{Multi-layer Perceptron}\\\
\mbox{NLP}&\mbox{Natural Language Processing}\\\ \mbox{PHM}&\mbox{Prognostic
and Health Management}\\\ \mbox{RNN}&\mbox{Recurrent Neural Networks}\\\
\mbox{SSM}&\mbox{State Space Model}\\\ \mbox{STORN}&\mbox{Stochastic
RNN~{}\cite[cite]{[\@@bibref{Number}{bayer2014learning}{}{}]}}\\\
\mbox{VAE}&\mbox{Variational
Autoencoder~{}\cite[cite]{[\@@bibref{Number}{fraccaro2018deep}{}{}]}}\\\
\mbox{VAE-RNN}&\mbox{Variational Autoencoder built upon RNN}\\\
\mbox{VI}&\mbox{Variational Inference}\\\ \mbox{VRNNaug}&\mbox{Variational RNN
via Augmented Recurrent Input Space}\end{array}$

##  Appendix B Notations

$\begin{array}[]{ll}\mathcal{B}&\mbox{Subset of short time series chunks}\\\
\mathcal{D}&\mbox{Time series dataset}\\\ \mathcal{D}_{c}&\mbox{Converted
short time series chunks}\\\ \mathcal{D}_{c}^{j}&\mbox{$j$-th short time
series chunk}\\\ d_{\mathbf{u}}&\mbox{Dimensionality of input signal}\\\
{d_{\mathbf{z}}}&\mbox{Dimensionality of latent state}\\\
d_{\mathbf{h}}&\mbox{Dimensionality of hidden state}\\\
d_{\mathbf{y}}&\mbox{Dimensionality of output}\\\ F&\mbox{Forecast horizon}\\\
f&\mbox{Transition function}\\\ g&\mbox{Measurement function}\\\
\mathbf{h}_{t}&\mbox{Hidden state of RNN or SSM at time $t$}\\\ J&\mbox{Number
of short time series chunks}\\\ K&\mbox{Number of Monte Carlo samples}\\\
\mathcal{L}&\mbox{Evidence lower bound calculated on $\mathcal{D}$}\\\
\tilde{\mathcal{L}}&\mbox{Unbiased estimation of evidence lower bound on
$\mathcal{D}_{c}$}\\\ \mathcal{L}^{j}&\mbox{Evidence lower bound calculated on
$\mathcal{D}_{c}^{j}$}\\\ \mathcal{M}&\mbox{Sequence model}\\\
\mathcal{N}(.,.)&\mbox{Gaussian distribution}\\\ p50&\mbox{50\%-quantile loss
of future forecasting}\\\ p90&\mbox{90\%-quantile loss of future
forecasting}\\\ \mathrm{QL}_{\rho}&\mbox{$\rho$-quantile loss}\\\
T&\mbox{Sequence length of training time series}\\\ \mathbf{u}_{t}&\mbox{Input
signal at time $t$}\\\ \bar{\mathbf{z}}_{t}&\mbox{Recurrent latent state until
time $t$}\\\ W&\mbox{Chunk size}\\\ \mathbf{y}_{t}&\mbox{Output observation at
time $t$}\\\ \bar{\mathbf{y}}_{t}&\mbox{Recurrent output until time $t$}\\\
\hat{\mathbf{y}}_{t}&\mbox{Output sample at time $t$}\\\
\dot{\mathbf{y}}_{t}&\mbox{Hybrid output at time $t$}\\\
\mathbf{z}_{t}&\mbox{Latent state at time $t$}\\\
\hat{\mathbf{z}}_{t}&\mbox{Latent sample at time $t$}\\\
\bar{\mathbf{z}}_{t}&\mbox{Recurrent latent state until time $t$}\\\
\bm{\epsilon}&\mbox{Random variable following unit Gaussian distribution}\\\
\phi^{u}(.),\phi^{z}(.),\phi^{y}(.)&\mbox{Recurrent mappings}\\\
\bm{\mu}_{\mathbf{z}_{t}},\bm{\nu}_{\mathbf{z}_{t}},\bm{\Sigma}_{\mathbf{z}_{t}}&\mbox{Mean,
variance and diagonalized variance of $q(\mathbf{z}_{t}|.)$}\\\
\bm{\mu}_{\mathbf{y}_{t}},\bm{\nu}_{\mathbf{y}_{t}},\bm{\Sigma}_{\mathbf{y}_{t}}&\mbox{Mean,
variance and diagonalized variance of $p(\mathbf{y}_{t}|.)$}\\\
\odot&\mbox{Element-wise product operator}\end{array}$

##  Appendix C Experimental details

Toy cases. For the two toy cases, we adopt the same settings for the following
system identification benchmarks except that (i) the chunk size is set as
$W=133$ and the latent dimensionality takes $d_{\mathbf{z}}=20$ on the
motorcycle dataset, and (ii) the Adam optimizer is ran over 200 epochs.

System identification benchmarks. The experimental configurations on the eight
system identification datasets ($\mathtt{actuator}$, $\mathtt{ctank}$,
$\mathtt{drive}$, $\mathtt{dryer}$, $\mathtt{f16gvt}$, $\mathtt{furnace}$,
$\mathtt{tank}$, and $\mathtt{hammerstein}$) are detailed as below.

Firstly, we perform standardization on the sequence inputs and outputs of time
series data to have zero mean and unit variance. We then split the whole data
by choosing the first 50% as training set, the middle 20% as validation set,
and the final 30% as testing set. Besides, we use the shingling strategy [42]
to further split the training and validation datasets with the chunk size
(i.e., sequence length) $W=64$, thus resulting into many short time series
sets.

As for model configuration of the proposed VRNNaug, we build three GRU models
in (16)-(18) using a single hidden layer with 100 units. The weight parameters
of these GRU models are initialized by the orthogonal initialization [57], and
the bias parameters are initialized with zeros. As for the MLPs in VRNNaug, we
adopt the fully-connected (FC) neural network with three hidden layers, the
ReLU activation, and skip connection [58] if possible. The number of units for
the hidden layers is set as $\mathrm{max}[d_{\mathbf{x}},50]$ wherein
$d_{\mathbf{x}}$ is the number of input features. Finally, for the VAE
structure in Fig. 1, we set the latent dimensionality of $\mathbf{z}$ as
$d_{\mathbf{z}}=10$.

As for model training, the Adam optimizer is employed with the batch size of
$|\mathcal{B}|=128$ and the maximum number of epochs as 100. Particularly, we
have a scheduled learning strategy. That is, we start using a learning rate of
$1\times 10^{-3}$ within the first ten epochs, and then check it every ten
epochs: if the validation loss cannot decrease within the last ten epochs, we
use a decreasing factor 0.5 to adjust the learning rate. We stop the training
when (i) it reaches the maximum number of epochs; or (ii) the adjusted
learning rate is below $1\times 10^{-6}$.

As for model forecasting, we start from cold and recursively sample $K=100$
points from the predictive distributions of deep stochastic SSM models at the
future time period $[T+1,T+F]$, and then use the prediction samples to
evaluate the related $p50$ and $p90$ quantile criteria.

It is finally notable that we repeat the experiments on each time series task
ten times by starting with different random seeds in order to comprehensively
evaluate the performance of deep sequence models.

Centrifugal compressor sensor data. After obtaining the processed compressor
sequence data with $d_{\mathbf{u}}=10$ and $d_{\mathbf{y}}=1$, we adopt the
similar experimental settings on the above system identification benchmarks
for model training and forecasting. The differences are that we here split the
whole data by choosing the first 80% as training set, the middle 10% as
validation set, and the final 10% as testing set; and we use the shingling
strategy to split the training and validation datasets with the chunk size
$W=128$.

##  Appendix D Illustration of VRNNaug on benchmarks

Fig. 8 depicts the predictive distributions of the proposed VRNNaug on eight
system identification benchmarks, respectively. Note that for the benchmarks
with many test time points (e.g., $F\geq 500$), we only illustrate the
predictions at the first 500 time points. It is observed that the predictions
of VRNNaug well describe the sequence behaviors on most benchmarks.

Figure 8: Illustration of the VRNNaug forecasting on eight system
identification benchmarks. The red curve represents 50-quantile prediction,
while the associated shaded region indicates the interval bounded by
95-quantile and 5-quantile predictions.

## References

## References

  * [1] K. J. Åström, P. Eykhoff, System identification—A survey, Automatica 7 (2) (1971) 123–162.
  * [2] M. Verhaegen, V. Verdult, Filtering and system identification: A least squares approach, Cambridge university press, 2007.
  * [3] R. Zhao, R. Yan, Z. Chen, K. Mao, P. Wang, R. X. Gao, Deep learning and its applications to machine health monitoring, Mechanical Systems and Signal Processing 115 (2019) 213–237.
  * [4] W. Yu, I. Y. Kim, C. Mechefske, Remaining useful life estimation using a bidirectional recurrent neural network based autoencoder scheme, Mechanical Systems and Signal Processing 129 (2019) 764–780.
  * [5] Y. Lei, B. Yang, X. Jiang, F. Jia, N. Li, A. K. Nandi, Applications of machine learning to machine fault diagnosis: A review and roadmap, Mechanical Systems and Signal Processing 138 (2020) 106587.
  * [6] A. Tadjer, A. Hong, R. B. Bratvold, Machine learning based decline curve analysis for short-term oil production forecast, Energy Exploration & Exploitation (2021) 01445987211011784.
  * [7] R. Hyndman, A. B. Koehler, J. K. Ord, R. D. Snyder, Forecasting with exponential smoothing: the state space approach, Springer Science & Business Media, 2008.
  * [8] G. E. Box, G. M. Jenkins, Some recent advances in forecasting and control, Journal of the Royal Statistical Society: Series C (Applied Statistics) 17 (2) (1968) 91–109.
  * [9] J. Durbin, S. J. Koopman, Time series analysis by state space methods, Oxford university press, 2012.
  * [10] T. B. Schön, A. Wills, B. Ninness, System identification of nonlinear state-space models, Automatica 47 (1) (2011) 39–49.
  * [11] H. Hanachi, J. Liu, A. Banerjee, Y. Chen, Sequential state estimation of nonlinear/non-gaussian systems with stochastic input for turbine degradation estimation, Mechanical Systems and Signal Processing 72 (2016) 32–45.
  * [12] L. Rabiner, B. Juang, An introduction to hidden Markov models, IEEE ASSP Magazine 3 (1) (1986) 4–16.
  * [13] Y. Liao, L. Zhang, C. Liu, Uncertainty prediction of remaining useful life using long short-term memory network based on bootstrap method, in: IEEE International Conference on Prognostics and Health Management, IEEE, 2018, pp. 1–8.
  * [14] R. E. Kalman, A new approach to linear filtering and prediction problems, Journal of Basic Engineering 82 (1) (1960) 35–45.
  * [15] M. West, J. Harrison, Bayesian forecasting and dynamic models, Springer Science & Business Media, 2006.
  * [16] C. K. Williams, C. E. Rasmussen, Gaussian processes for machine learning, MIT press Cambridge, MA, 2006.
  * [17] R. Frigola, Y. Chen, C. Rasmussen, Variational Gaussian process state-space models, in: Advances in Neural Information Processing Systems, Vol. 27, 2014, pp. 3680–3688.
  * [18] A. D. Ialongo, M. van der Wilk, J. Hensman, C. E. Rasmussen, Overcoming mean-field approximations in recurrent Gaussian process models., in: International Conference on Machine Learning, 2019, pp. 2931–2940.
  * [19] H. Liu, Y.-S. Ong, X. Shen, J. Cai, When gaussian process meets big data: A review of scalable GPs, IEEE Transactions on Neural Networks and Learning Systems 31 (11) (2020) 4405–4423.
  * [20] E. Snelson, Z. Ghahramani, Sparse gaussian processes using pseudo-inputs, in: Advances in Neural Information Processing Systems, MIT Press, 2006, pp. 1257–1264.
  * [21] M. Titsias, Variational learning of inducing variables in sparse Gaussian processes, in: Artificial Intelligence and Statistics, 2009, pp. 567–574.
  * [22] J. Hensman, N. Fusi, N. D. Lawrence, Gaussian processes for big data, in: Uncertainty in Artificial Intelligence, Citeseer, 2013, p. 282.
  * [23] H. Liu, J. Cai, Y. Wang, Y.-S. Ong, Generalized robust Bayesian committee machine for large-scale Gaussian process regression, in: International Conference on Machine Learning, 2018, pp. 1–10.
  * [24] K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, Y. Bengio, Learning phrase representations using RNN encoder–decoder for statistical machine translation, in: Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing, 2014, pp. 1724–1734.
  * [25] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural Computation 9 (8) (1997) 1735–1780.
  * [26] I. Sutskever, O. Vinyals, Q. V. Le, Sequence to sequence learning with neural networks, arXiv preprint arXiv:1409.3215.
  * [27] J. Gehring, M. Auli, D. Grangier, D. Yarats, Y. N. Dauphin, Convolutional sequence to sequence learning, in: International Conference on Machine Learning, PMLR, 2017, pp. 1243–1252.
  * [28] Z. C. Lipton, J. Berkowitz, C. Elkan, A critical review of recurrent neural networks for sequence learning, arXiv preprint arXiv:1506.00019.
  * [29] W. Yu, I. Y. Kim, C. Mechefske, Analysis of different rnn autoencoder variants for time series classification and machine prognostics, Mechanical Systems and Signal Processing 149 (2021) 107322.
  * [30] D. Salinas, V. Flunkert, J. Gasthaus, T. Januschowski, DeepAR: Probabilistic forecasting with autoregressive recurrent networks, International Journal of Forecasting 36 (3) (2020) 1181–1191.
  * [31] L. Li, J. Yan, X. Yang, Y. Jin, Learning interpretable deep state space model for probabilistic time series forecasting, arXiv preprint arXiv:2102.00397.
  * [32] S. S. Rangapuram, M. W. Seeger, J. Gasthaus, L. Stella, Y. Wang, T. Januschowski, Deep state space models for time series forecasting, Advances in Neural Information Processing Systems 31 (2018) 7785–7794.
  * [33] A. K. Yanchenko, S. Mukherjee, Stanza: A nonlinear state space model for probabilistic inference in non-stationary time series, arXiv preprint arXiv:2006.06553.
  * [34] M. Al-Shedivat, A. G. Wilson, Y. Saatchi, Z. Hu, E. P. Xing, Learning scalable deep kernels with recurrent structure, The Journal of Machine Learning Research 18 (1) (2017) 2850–2886.
  * [35] Q. She, A. Wu, Neural dynamics discovery via Gaussian process recurrent neural networks, in: Uncertainty in Artificial Intelligence, PMLR, 2020, pp. 454–464.
  * [36] D. P. Kingma, M. Welling, Auto-encoding variational Bayes, arXiv preprint arXiv:1312.6114.
  * [37] J. Bayer, C. Osendorfer, Learning stochastic recurrent networks, arXiv preprint arXiv:1411.7610.
  * [38] R. G. Krishnan, U. Shalit, D. Sontag, Deep kalman filters, arXiv preprint arXiv:1511.05121.
  * [39] R. G. Krishnan, U. Shalit, D. Sontag, Structured inference networks for nonlinear state space models, in: Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, 2017, pp. 2101–2109.
  * [40] M. Fraccaro, Deep latent variable models for sequential data, Ph.D. thesis, Department of Applied Mathematics and Computer Science, Technical University of Denmark (2018).
  * [41] D. Gedon, N. Wahlström, T. B. Schön, L. Ljung, Deep state space models for nonlinear system identification., arxiv:eess.SY.
  * [42] J. Leskovec, A. Rajaraman, J. D. Ullman, Mining of massive data sets, Cambridge university press, 2014.
  * [43] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980.
  * [44] S. Bengio, O. Vinyals, N. Jaitly, N. Shazeer, Scheduled sampling for sequence prediction with recurrent neural networks, in: Advances in Neural Information Processing Systems, 2015, pp. 1171–1179.
  * [45] W. Zhang, Y. Feng, F. Meng, D. You, Q. Liu, Bridging the gap between training and inference for neural machine translation, in: Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, 2019, pp. 4334–4343.
  * [46] S. Bai, J. Z. Kolter, V. Koltun, An empirical evaluation of generic convolutional and recurrent networks for sequence modeling, arXiv preprint arXiv:1803.01271.
  * [47] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, I. Polosukhin, Attention is all you need, in: Advances in Neural Information Processing Systems, 2017, pp. 6000–6010.
  * [48] G. Huang, Z. Liu, L. Van Der Maaten, K. Q. Weinberger, Densely connected convolutional networks, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 4700–4708.
  * [49] Y. Chen, Y. Kang, Y. Chen, Z. Wang, Probabilistic forecasting with temporal convolutional neural network, Neurocomputing 399 (2020) 491–501.
  * [50] L. Lennart, System identification: Theory for the user, PTR Prentice Hall, Upper Saddle River, NJ 28.
  * [51] B. W. Silverman, Some aspects of the spline smoothing approach to non-parametric regression curve fitting, Journal of the Royal Statistical Society: Series B (Methodological) 47 (1) (1985) 1–21.
  * [52] M. Schoukens, P. Mattson, T. Wigren, J.-P. Noel, Cascaded tanks benchmark combining soft and hard nonlinearities, in: Workshop on nonlinear system identification benchmarks, 2016, pp. 20–23.
  * [53] J.-P. Noël, M. Schoukens, F-16 aircraft benchmark based on ground vibration test data, in: 2017 Workshop on Nonlinear System Identification Benchmarks, 2017, pp. 19–23.
  * [54] M. Schoukens, J.-P. Noel, Wiener-hammerstein benchmark with process noise, in: Workshop on nonlinear system identification benchmarks, 2016, pp. 15–19.
  * [55] X. Shengli, J. Xiaomo, H. Jinzhi, Y. Shuhua, W. Xiaofang, Bayesian wavelet pca methodology for turbomachinery damage diagnosis under uncertainty, Mechanical Systems and Signal Processing 80 (Dec) (2016) 1–18.
  * [56] I. Daubechies, Ten lectures on wavelets, SIAM, 1992.
  * [57] A. M. Saxe, J. L. McClelland, S. Ganguli, Exact solutions to the nonlinear dynamics of learning in deep linear neural networks, arXiv preprint arXiv:1312.6120.
  * [58] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, 2016, pp. 770–778.