# Probabilistic Learning of Multivariate Time Series with Temporal
Irregularity

Yijun LI
<EMAIL_ADDRESS>& Cheuk Hang LEUNG
<EMAIL_ADDRESS>

School of Data Science
City University of Hong Kong
& Qi WU
<EMAIL_ADDRESS>The corresponding author.

###### Abstract

Multivariate sequential data collected in practice often exhibit temporal
irregularities, including nonuniform time intervals and component
misalignment. However, if uneven spacing and asynchrony are endogenous
characteristics of the data rather than a result of insufficient observation,
the information content of these irregularities plays a defining role in
characterizing the multivariate dependence structure. Existing approaches for
probabilistic forecasting either overlook the resulting statistical
heterogeneities, are susceptible to imputation biases, or impose parametric
assumptions on the data distribution. This paper proposes an end-to-end
solution that overcomes these limitations by allowing the observation arrival
times to play the central role of model construction, which is at the core of
temporal irregularities. To acknowledge temporal irregularities, we first
enable unique hidden states for components so that the arrival times can
dictate when, how, and which hidden states to update. We then develop a
conditional flow representation to non-parametrically represent the data
distribution, which is typically non-Gaussian, and supervise this
representation by carefully factorizing the log-likelihood objective to select
conditional information that facilitates capturing time variation and path
dependency. The broad applicability and superiority of the proposed solution
are confirmed by comparing it with existing approaches through ablation
studies and testing on real-world datasets.

Robotics Autonomous systems Financial markets Multivariate time series
Asynchronous data Irregular sampling Non-uniform time interval Serial
dependence Multivariate dependence Neural ODEs Recurrent neural networks Gated
Recurrent Unit Long short-term memory Normalizing flow models Conditional
distribution Joint distribution Time-varying dynamics Likelihood estimation
Probabilistic forecasting

##  1 Introduction

Irregularly sampled time series are sequential recordings with unevenly spaced
intervals between successive observation points. In the multivariate context,
these time series can exhibit asynchrony because observations of individual
component series may not align due to unevenly spaced intervals. These
temporal irregularities arise whenever "missing" values occur in applications
[kremer2010detecting, spadon2021pay, xu2013modeling]. They could occur due to
insufficient observation. For instance, in applications such as autonomous
systems, tracking the position and speed of a group of moving objects often
encounter missing values when certain objects travel behind obstacles or when
the return signals attenuate differently. In this case, the observation
process is constrained, hindering regular monitoring of continuously varying
phenomena. They could also occur if observable events themselves arrive at
inherently random times or with varying frequencies. For example, a liquid
stock can trade on a regular trading day at a millisecond time scale. However,
options contracts written on the same underlying stock may only trade once
every few minutes. In this case, temporal irregularity is an intrinsic feature
of the data-generating process. It is the liquidity disparity that results in
the differences in time rather than a lack of observation.

Temporal irregularities, namely non-uniform time intervals, varying number of
observations, and lack of alignment, pose challenges to the probabilistic
forecasting of multivariate time series (MTS) data [shukla2020survey].
Probabilistic forecasting [benidis2018deep] involves specifying the target
variables’ joint distribution, learning its representation from data, and
forecasting its time evolution. It is a fundamental task that includes the
prediction problem of the mean, median, or specific quantile level as special
cases. It also provides confidence intervals crucial for mission-critical
downstream decision-making tasks [ban2018machine]. However, the presence of
temporal irregularities makes it challenging to model the distribution of MTS
data, a defining characteristic of which is the existence of dependence
structure [zhang2021memory].

At the component level, the dependence structure of MTS data is serial,
meaning the marginal distributions of a given component variable conditional
on different time points are not independent. In addition, different
components can have different degrees and different natures of serial
dependencies. Some exhibit short-ranged but nonlinear serial dependencies,
while others are linear but long-ranged [zhang2021memory]. At the group level,
the dependence structure refers to quantities such as the copula function of
the joint distribution [zhiri2022unified]. Specific components may depend more
strongly on each other than on others. Moreover, both the component-level
serial dependencies and the group-level inter-dependencies can vary over time;
for instance, the autocorrelations of individual components and the copula
function of the group data can be path dependent.

In other words, the distributional properties of MTS data can exhibit
heterogeneities and vary over time, regardless of whether they are irregularly
sampled. However, in cases where temporal irregularities are inherent to the
data generation process, these irregularities reflect the underlying cause of
these properties. Therefore, attributing them to missing values would be
inappropriate. An ideal framework for probabilistic learning should
acknowledge the intrinsic nature of uneven spacing and misalignments in the
data and represent the associated data distribution accordingly. Before
introducing our approach, we first review existing approaches to handling
temporal irregularity and joint distribution modeling.

###  1.1 On Handling Temporal Irregularity

There are three primary approaches to handling irregular sampling. The first
approach involves converting an irregularly sampled time series into one with
evenly spaced time intervals before making predictions. The discretization
method selects a larger uniformly spaced time interval, with each interval
containing several observations, and computes the mean of these observations
to represent the interval’s value. However, this method loses local
information due to averaging. In contrast, the imputation method interpolates
the missing values of the lower-frequency components instead of averaging the
higher-frequency components. It keeps the local information of the higher-
frequency components intact and uses models such as the Gaussian Process
regression model [williams2006gaussian], the Recurrent Neural Networks (RNNs)
[cao2018brits], and the Generative Adversarial Networks (GANs)
[luo2018multivariate] to impute the missing values of the lower-frequency
components.

The next approach proposes using end-to-end models to avoid the "interpolate
first and predict later" idea. This approach modifies classical recurrent
architectures to encode the information embedded in irregular temporal
patterns. For example, Che et al. [che2018recurrent] added an exponential
decay mechanism in the hidden state. Neil et al. [neil2016phased] extended the
Long short-term memory (LSTM) unit by adding a new time gate controlled by a
parametrized oscillation with a frequency range. Additionally, Mozer et al.
[mozer2017discrete] incorporated multiple time scales of the hidden state and
made a context-dependent selection of time scales for information storage and
retrieval. However, harnessing the uneven spacing property in the data and its
asynchronous property is indirect using this approach and is less coherent
compared to the Neural Ordinary Differential Equations (Neural ODEs) approach.

Chen et al. [chen2018neural] introduced the Neural ODE framework by extending
discrete neural networks into continuous-time networks, which makes it a
natural candidate for handling data with arbitrarily small time intervals. For
instance, De Brouwer et al. [de2019gru] embedded the Neural ODEs in the
classical Gated Recurrent Unit (GRU) cell and derived the dynamics of the
hidden state; Rubanova et al. [rubanova2019latent] proposed the Latent ODE and
were the first to embed Neural ODEs in a Variational Autoencoder to address
the problem of irregularly sampled time series. Unlike the classical GRU cell,
which keeps the hidden state constant in the absence of observations, the
continuous-time GRU cell learns to evolve the hidden state using Neural ODEs.
Building upon similar ideas, Lechner et al. [lechner2020learning] transformed
the standard LSTM into a continuous version to address the issues of gradient
vanishing and exploding. In section 2.1 and 2.2, we will provide more
background on Neural ODEs and how to utilize them to model irregularly sampled
data.

###  1.2 On Modeling Joint Data Distribution

For point estimation tasks, vanilla recurrent architectures, including RNN,
GRU, and LSTM, can capture different aspects of the aforementioned properties
of the MTS data. However, they are not directly applicable to distribution
prediction due to the deterministic nature of the transition functions of
their hidden states, which do not account for modeling uncertainties
[chung2015recurrent].

One class of models modifies the output function of neural networks to model
the joint distribution or quantile function. For example, the models in
[salinas2020deepar, de2019gru, salinas2019high] assume the data-generating
process follows parametric distribution, such as the multivariate Gaussian
(for continuous variables) and multivariate Negative Binomial (for discrete
variables). To make the distribution time-varying, they condition the
distribution parameters on hidden states and forecast them dynamically.
Alternatively, researchers use quantile regression methods to fit the quantile
function of the joint distribution. They use the quantile loss [wen2017multi,
yan2019cross] or the continuous ranked probability score (CRPS)
[laio2007verification] as the objective function to train the model and
predict multiple quantile points simultaneously conditional on the hidden
states.

Recently, unsupervised deep learning models have been utilized to learn the
joint distribution of data, including combining Variational Autoencoders
(VAEs) with RNNs [chung2015recurrent] and integrating Normalizing Flows
[rasul2020multivariate] into RNNs. Among these, two notable models are Real
NVP (Real-valued Non-Volume Preserving) and MAF (Masked Autoregressive Flow).
Real NVP defines a series of masked affine transformations, each of which
splits the input data into two parts, applying a nonlinear transformation to
one part while leaving the other unchanged. These invertible transformations
facilitate both forward and inverse passes during training and generation. On
the other hand, MAF is a normalizing flow model that employs autoregressive
models to parameterize the inverse transformations. It transforms the input
data through a series of invertible masked autoregressive layers, with each
layer predicting the parameters of the conditional distribution for each
dimension based on the previously transformed dimensions. By conditioning on
the previous dimensions, MAF captures complex dependencies in the data
distribution.

The flow models are flexible in capturing intricate and evolving dependence
structures and impose minimal assumptions about the functional form of the
joint data distribution. These characteristics make them attractive for
dealing with complex data, although they do not specifically address the
structural aspects of irregular sampling. In section 2.3, we shall further
detail the background of representing data distribution using the normalizing
flow approach.

###  1.3 Our Approach and Contributions

Discussions in 1.1 and 1.2 unveil dislocations and disparities among ideas of
handling temporal irregularities in the data and ideas to model its joint
distribution. This paper bridges this gap by introducing a deep learning
solution called the Recurrent Flow Network (RFN). It can seamlessly integrate
the treatment of temporal irregularities with the learning of joint data
distribution. Its novelties are as follows.

(i) The proposed RFN framework formulates a two-layer representation that
distinguishes marginal learning of component dynamics from multivariate
learning of joint data distribution. It is a versatile methodology that can be
trained end-to-end and accommodate synchronous and asynchronous data
structures. It also broadly applies to underlying recurrent architectures,
including RNNs, GRUs, and LSTMs. Once the joint data distribution is learned,
it is ready for sampling despite the distribution being non-parametrically
represented via neural networks.

(ii) The marginal learning block addresses the component-level heterogeneity
in the presence of unevenly spaced intervals. It achieves this by assigning
unique sets of hidden states to individual variables, overcoming the drawback
of disregarding statistical heterogeneities when hidden states are shared
(e.g., in models like GRU-D [che2018recurrent], ODERNN [rubanova2019latent],
GRU-ODE-Bayes [de2019gru], and ODELSTM [lechner2020learning]), and by allowing
the arrival times of component observations to determine how hidden states are
updated, effectively avoiding data imputation.

(iii) The multivariate learning block then resolves the struggle faced by
existing distribution models [salinas2020deepar, de2019gru, salinas2019high]
in achieving a non-parametric representation of non-Gaussian data
distribution, simultaneously with a flexible choice of information to generate
time variation. The conditional CNF (Continuous Normalizing Flow)
representation we develop enables one to choose what information to use to
drive the time variation of the base distribution and the flow map without
compromising any non-parametric capacity to represent the non-Gaussian data
distribution.

(iv) Building upon (ii) and (iii), we strategically condition the log-
likelihood objective on the observation times and along the component
direction. This conditioning structure enables the optimizer o fully
acknowledge and account for both the uneven spacing aspect and the asynchrony
aspect of temporal irregularity in the MTS data. By conditioning on the
observation times, the RFN ensures that the model incorporates the specific
time points at which the data is observed. Additionally, conditioning along
the component direction acknowledges the asynchrony between different
variables, which arises exclusively in the multivariate context.

We validate the novelties mentioned above through three ablation studies using
simulations and demonstrate the overall performance of the RFNs on three real-
world datasets. The ablation studies simulate sample paths of a multivariate
correlated Geometric Brownian Motion process to investigate the effects of
sharing hidden states, imputing unobserved variables, and making different
choices of conditional information. Meanwhile, the experimental datasets
include the physical activities of the human body from the MuJoCo module
[rubanova2019latent], the climate records of weather from the USHCN dataset
[menne2010long], and a three-year-long transaction record from HKEX comprising
minute-level transactions for the stock with the largest market capitalization
and its options. We compare four baseline models between their performance in
vanilla forms and the performances utilizing the RFN specification. The
vanilla versions of these four baselines are GRU-D, ODERNN, GRU-ODE-Bayes, and
ODELSTM. They are all designed to target multivariate irregularly sampled time
series. The results show that the RFN framework has broad applicability and
outperforms existing approaches significantly.

The rest of this paper is organized as follows. Section 2 provides the
background knowledge that the RFNs build upon. Section 3 presents the data
structure and necessary notations. Section 4 introduces model specifications
for the marginal learning block. Sections 5 and 6 present methodologies of the
multivariate learning block for the synchronous case and the asynchronous
case. Section 7 discusses the ablation studies and the experiments. Finally,
Section 8 concludes the paper. Additional formulations, proofs, extended
numerics, and algorithms are provided in the Appendix A, B, C, and D
respectively.

##  2 Backgrounds

Understanding the working mechanism of the proposed RFN model requires
knowledge of Neural ODEs, their applications to model irregularly sampled
data, and the flow representation of distributions. This section summarizes
these subjects to make the paper self-contained. Throughout, we denote random
variables as follows: $X$ for scalar, $\mathsf{X}$ for vector, and
$\mathbf{X}$ for matrix, and their corresponding sample values as $x$,
$\mathsf{x}$, and $\mathbf{x}$, respectively.

###  2.1 Neural Ordinary Differential Equations

Neural ODEs were developed as the continuous limit of the ResNet model. The
ResNet model solves the degradation problem of neural network [he2016deep,
srivastava2015highway, simonyan2015very, szegedy2015going] where researchers
noticed that, as the network layers go deeper, the training loss begins to
increase steadily once the network depth crosses a certain threshold. Consider
a $L$-layer network, with $\mathsf{x}_{0}$ being the input and
$\mathsf{x}_{l}$ being the output of each layer $l\in\\{1,\cdots,L\\}$.
Instead of learning the mapping from $\mathsf{x}_{0}$ to $\mathsf{x}_{L}$
directly, ResNet learns the difference between the input and output of each
layer,

$\mathsf{x}_{l}=\mathsf{x}_{l-1}+\mathsf{f}_{\theta_{l}}(\mathsf{x}_{l-1}).$
(1)

The shortcut connections ensure that when adding one more layer, the worst
scenario is that the learned residual is zero, resulting in the corresponding
layer being an identity map. As a result, the training error will not increase
when the model becomes deeper.

Chen et al. [chen2018neural] proposed that taking the limit of the number of
layers to infinity shall turn discrete layers into continuous layers. The
resulting continuous limit of the recursive equation (1) is an ODE

$\frac{d\mathsf{x}(l)}{dl}=\mathsf{f}_{\theta}(\mathsf{x}(l),l).$ (2)

Solving (2) with initial condition $\mathsf{x}(0)$ is equivalent to the
forward pass of ResNet. One can use numerical methods such as the Euler and
the Runge–Kutta methods to solve (2). Similar to the role of backpropagation
in discrete neural networks, the adjoint sensitivity method
[pontryagin1987mathematical] computes gradients by solving an augmented ODE
backward and updates parameters in continuous neural networks. The advantage
of the adjoint sensitivity method is that it does not require storing
intermediate quantities from the forward pass, enabling the training of Neural
ODEs with constant memory.

###  2.2 Irregularly Sampled Time Series Modeling

Building on the insights from the Neural ODEs, researchers applied the
continuous extension idea to model irregularly sampled MTS data using other
recurrent architectures, such as ODERNN [rubanova2019latent], GRU-ODE-Bayes
[de2019gru], and ODELSTM [lechner2020learning]. For instance, De Brouwer et
al. [de2019gru] derived the continuous-time dynamics of the hidden states for
the vanilla GRU cell and used ODE solvers to obtain its solution. The vanilla
GRU has the following updating formulas

$\displaystyle\mathsf{r}_{t}$
$\displaystyle=\sigma(\mathbf{w}_{r}\mathsf{x}_{t}+\mathbf{u}_{r}\mathsf{h}_{t-1}+\mathsf{b}_{r}),$
$\displaystyle\mathsf{z}_{t}$
$\displaystyle=\sigma(\mathbf{w}_{z}\mathsf{x}_{t}+\mathbf{u}_{z}\mathsf{h}_{t-1}+\mathsf{b}_{z}),$
$\displaystyle\tilde{\mathsf{h}}_{t}$
$\displaystyle=\tanh(\mathbf{w}_{h}\mathsf{x}_{t}+\mathbf{u}_{h}(\mathsf{r}_{t}\odot\mathsf{h}_{t-1})+\mathsf{b}_{h}),$
$\displaystyle\mathsf{h}_{t}$
$\displaystyle=(1-\mathsf{z}_{t})\odot\tilde{\mathsf{h}}_{t}+\mathsf{z}_{t}\odot\mathsf{h}_{t-1},$
(3)

where $\mathsf{x}_{t}\in\mathbb{R}^{D}$ is the time-$t$ MTS data;
$\mathsf{r}_{t}$, $\mathsf{z}_{t}$, $\tilde{\mathsf{h}}_{t}$,
$\mathsf{h}_{t}\in\mathbb{R}^{H}$ are vectors denoting the reset gates, the
update gates, the candidate update gates, and the hidden states;
$\mathbf{w}_{\\{r,z,h\\}}\in\mathbb{R}^{H\times D}$,
$\mathbf{u}_{\\{r,z,h\\}}\in\mathbb{R}^{H\times H}$, and
$\mathsf{b}_{\\{r,z,h\\}}\in\mathbb{R}^{H}$ are the matrix parameters to be
trained; $\sigma$ and $\tanh$ are the sigmoid and the hyperbolic tangent
activation function; the operator $\odot$ denotes element-wise multiplication.

Subtracting $\mathsf{h}_{t-1}$ on both sides of (3) leads to

$\displaystyle\mathsf{h}_{t}-\mathsf{h}_{t-1}$
$\displaystyle=\left(1-\mathsf{z}_{t}\right)\odot\left(\tilde{\mathsf{h}}_{t}-\mathsf{h}_{t-1}\right).$

Then by taking the limit as the time difference between $t-1$ and $t$ tends to
zero, one arrives at the continuous-time dynamics of the hidden states

$\displaystyle\frac{d\mathsf{h}(t)}{dt}$
$\displaystyle=(1-\mathsf{z}(t))\odot(\tilde{\mathsf{h}}(t)-\mathsf{h}(t)).$
(4)

One can use (4) to evolve the hidden states during periods without
observations and (3) to update the hidden state at observation time points.
The continuous-time dynamics thus seamlessly interpolates hidden states
between any two consecutive times of observations, accommodating unevenly
spaced intervals. Figure 1 illustrates this process.

Figure 1: The evolution and update of hidden state in GRU-ODE-Bayes. Solid
dots represent observed data points of $X^{1}_{t},X^{2}_{t},X^{3}_{t}$. GRU-
ODE evolves the hidden state in the absence of observations for all variables,
while GRU-Bayes updates the hidden state when at least one variable has an
observed value.

###  2.3 Unconditional Normalizing Flow

For tasks related to probabilistic forecasts, one needs a representation of
the data distribution. Let $p(\mathsf{x})$, $\mathsf{x}\in\mathbb{R}^{D}$ be
the probability density of the data-generating distribution and
$p(\mathsf{z})$, $\mathsf{z}$ $\in\mathbb{R}^{D}$ be the probability density
of the base distribution which is typically set as the standard normal, i.e.,
$\mathsf{Z}\sim N(\mathsf{0},\mathbb{I}_{D})$. The idea of the normalizing
flow model is to find a differentiable bijective function $\mathsf{f}$ which
can map samples from $\mathsf{Z}$ to $\mathsf{X}$ [dinh2014nice,
rezende2015variational]

$\displaystyle\mathsf{f}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D};\quad\mathsf{f}(\mathsf{z})=\mathsf{x}.$

In the discrete formulation, $\mathsf{f}$ is typically specified as a sequence
of neural networks,
$\mathsf{f}=\mathsf{f}_{1}\circ\cdots\circ\mathsf{f}_{M-1}\circ\mathsf{f}_{M}.$
Given $\mathsf{f}$ and $p(\mathsf{z})$, one obtains $p(\mathsf{x})$ and $\log
p(\mathsf{x})$ as

$\displaystyle
p(\mathsf{x})=p(\mathsf{f}^{-1}(\mathsf{x}))\left|\operatorname{det}[\partial_{\mathsf{z}}\mathsf{f}(\mathsf{z})]\right|^{-1},$
$\displaystyle\log p(\mathsf{x})=\log
p(\mathsf{f}^{-1}(\mathsf{x}))-\log\left|\operatorname{det}[\partial_{\mathsf{z}}\mathsf{f}(\mathsf{z})]\right|,$

via the change of variable formula, where
$\operatorname{det}[\partial_{\mathsf{z}}\mathsf{f}(\mathsf{z})]$ is the
determinant of the Jacobian matrix of $\mathsf{f}$ whose $(i,j)$ entries are
$\partial_{z_{j}}x_{i};$ $i,j=1,\cdots,D$.

However, designing the architectures of $\mathsf{f}_{1},\cdots,\mathsf{f}_{M}$
is challenging because they need to satisfy three conditions: being bijective,
differentiable, and facilitating the computation of the determinant of the
Jacobian of the function $\mathsf{f}$. Several models have thus been proposed,
such as Real NVP [dinh2016density], MAF [papamakarios2017masked], and GLOW
[kingma2018glow]. See [kobyzev2020normalizing] for a complete introduction to
the flow models.

The continuous normalizing flow model (CNF) [grathwohl2018ffjord] offers a
solution to this challenge by extending the composition of discrete maps into
a continuous map, whose differential form reads

$\displaystyle\frac{\partial\mathsf{z}(s)}{\partial
s}=\mathsf{f}(\mathsf{z}(s),s;\theta),s\in[s_{0},s_{1}],$ (5a)
$\displaystyle\text{where}\quad\mathsf{z}(s)|_{s=s_{0}}=\mathsf{z},\,\,\mathsf{z}(s)|_{s=s_{1}}=\mathsf{x}.$
Unlike the physical time $t$, $s$ is called the flow time of the dynamics (5).
At the initial flow time $s_{0}$, the value of the flow $\mathsf{z}(s_{0})$ is
set as $\mathsf{z}$, which samples from the base distribution $p(\mathsf{z})$
of the base random variable $\mathsf{Z}$. At the terminal flow time $s_{1}$,
the value of the flow $\mathsf{z}(s_{1})$ is set to equal $\mathsf{x}$, which
is the observed sample from the distribution of the true data-generating
distribution $p(\mathsf{x})$.

The discrete formulation requires careful design of the weight matrices of
$\mathsf{f}_{j},$ $1\leq j\leq M$ to be triangular to facilitate computing the
Jacobian’s determinant easily. The continuous formulation, however, replaces
the computation of the Jacobian’s determinant with relatively cheap trace
operations, thanks to the Instantaneous Change of Variables theorem
[chen2018neural]. The log-density of the continuous flow therefore follows

$\displaystyle\frac{\partial\log p(\mathsf{z}(s))}{\partial
s}=-\operatorname{Tr}\left[\partial_{\mathsf{z}(s)}\mathsf{f}\right],s\in[s_{0},s_{1}],$
(5b) $\displaystyle\text{where}\,\quad
p(\mathsf{z}(s))|_{s=s_{0}}=p(\mathsf{z}),\;p(\mathsf{z}(s))|_{s=s_{1}}=p(\mathsf{x}).$

Solving equations (LABEL:eqt:ffjord1) and (LABEL:eqt:ffjord2) together, we
have

$\begin{bmatrix}\mathsf{x}\\\ \log
p(\mathsf{x})\end{bmatrix}=\begin{bmatrix}\mathsf{z}\\\ \log
p(\mathsf{z})\end{bmatrix}+\int_{s_{0}}^{s_{1}}\begin{bmatrix}\mathsf{f}(\mathsf{z}(s),s;\theta)\\\
-\operatorname{Tr}\left[\partial_{\mathsf{z}(s)}\mathsf{f}\right]\end{bmatrix}ds.$
(5)

##  3 Data Structure

With enough background knowledge, we are ready to describe the data structure
pertaining to the RFNs. Consider an MTS dataset containing $N$ instances. Each
instance is a $D$-dimensional sample path of the true data-generating process.
All instances span the same $[0,T]$ period.

To account for the presence of temporal irregularity in a given instance
$i\in\\{1,\cdots,N\\}$, we first collect all time points at which at least one
component series has an observation and define this collection as the time
vector of observations,

$\mathsf{t}^{i}:=[t_{1}^{i},\cdots,t_{K_{i}}^{i}],\quad 0\leq
t_{1}^{i}\leq\cdots\leq t_{K_{i}}^{i}\leq T.$

For example, consider a simple dataset containing two univariate sample paths.
The first sample path has four observations at physical time $1,3,5,6$ in
units of seconds, and the second has three observations at $2,4,6$. Then
$K_{1}=4$ and $K_{2}=3$ with

$\displaystyle\mathsf{t}^{1}=[t^{1}_{1},\cdots,t^{1}_{K_{1}}]=[t^{1}_{1},t^{1}_{2},t^{1}_{3},t^{1}_{4}]=[1,3,5,6],$
$\displaystyle\mathsf{t}^{2}=[t^{2}_{1},\cdots,t^{2}_{K_{2}}]=[t^{2}_{1},t^{2}_{2},t^{2}_{3}]=[2,4,6].$

At a particular observation time $t\in\mathsf{t}^{i}$, we use
$\mathsf{x}_{i,t}\in\mathbb{R}^{D\times 1}$ to denote the time-$t$ sample
values of the random vector $\mathsf{X}_{t}$ in the $i^{\textrm{th}}$
instance, within which we use $x_{i,t}^{d}\in\mathbb{R}$ to denote the
$d^{\textrm{th}}$ component, which is the sample value of the
$d^{\textrm{th}}$ random scalar $X^{d}_{t}$,

$\mathsf{x}_{i,t}:=[x_{i,t}^{1},\cdots,x_{i,t}^{d},\cdots,x_{i,t}^{D}]^{\prime},\,t\in\mathsf{t}^{i};$

we then aggregate $\mathsf{x}_{i,t}$ from all observation times
$t=t^{i}_{1},t^{i}_{2},\cdots,t^{i}_{K_{i}}$ of the instance $i$ to form the
instance-level data matrix $\mathbf{x}_{i}\in\mathbb{R}^{D\times K_{i}}$,

$\mathbf{x}_{i}:=[\mathsf{x}_{i,t_{1}^{i}};\mathsf{x}_{i,t_{2}^{i}};\cdots;\mathsf{x}_{i,t_{K_{i}}^{i}}]^{\prime},\,i\in[1,\cdots,N];$

finally, the entire MTS dataset is the collection of all $N$ instances
$\\{\\{\mathbf{x}_{1}\\};\\{\mathbf{x}_{2}\\};\cdots;\\{\mathbf{x}_{N}\\}\\}$.

An instance $\mathbf{x}_{i}$ can be synchronous or asynchronous depending on
whether $x^{d}_{i,t}$ is observed or not for any combination of the variable
dimension $d\in\\{1,\cdots,D\\}$ and the observation time
$t\in\mathsf{t}^{i}$.

######  Definition 1.

Synchronous multivariate time series (Syn-MTS):
An instance $\mathbf{x}_{i}$ where all of its $D$ component series have
observations at each and every time points $t\in\mathsf{t}^{i}$ (see Figure
2(a)).

######  Definition 2.

Asynchronous multivariate time series (Asyn-MTS):
An instance $\mathbf{x}_{i}$ where at least one of its $D$ component series
does not have all observations at all time point $t\in\mathsf{t}^{i}$ (see
Figure 2(b)).

The following mask matrix $\mathbf{m}_{i}$ makes the distinction between Asyn-
MTS and Syn-MTS precise. For each instance $i$, given the data matrix
$\mathbf{x}_{i}$, its mask matrix is defined as

$\mathbf{m}_{i}:=[\mathsf{m}_{i,t_{1}^{i}};\mathsf{m}_{i,t_{2}^{i}};\cdots;\mathsf{m}_{i,t_{K_{i}}^{i}}]^{\prime},$

where each vector $\mathsf{m}_{i,t}$ denotes whether or not the constituent
component variables are observed at time $t\in\mathsf{t}^{i}$,

$\displaystyle\mathsf{m}_{i,t}:$
$\displaystyle=[m^{1}_{i,t},\cdots,m^{d}_{i,t},\cdots,m^{D}_{i,t}]^{\prime},\text{where}$
$\displaystyle m^{d}_{i,t}$ $\displaystyle=\begin{cases}1,\quad&\text{if
}x^{d}_{i,t}\text{ is observed};\\\ 0,\quad&\text{if }x^{d}_{i,t}\text{ is
unobserved}.\end{cases}$

In the sequel, we shall drop the instance script $i$ to lighten notations,
e.g., $t^{i}_{k}$, $\mathsf{x}_{i,t}$ and $\mathsf{m}_{i,t}$ shall become
$t_{k}$, $\mathsf{x}_{t}$ and $\mathsf{m}_{t}$, whenever the context is clear.

(a) Syn-MTS

(b) Asyn-MTS

Figure 2: (a) and (b) are sample path examples of Syn-MTS and Asyn-MTS where
observed data points are marked as solid circle dots. While the time intervals
between consecutive observation times are unevenly spaced in both cases,
component observations of the Syn-MTS sample path are always aligned. In
contrast, in the Asyn-MTS case, no observation time has complete observations.
This demonstrates that uneven spacing originates at the univariate level,
while asynchrony arises exclusively in the multivariate context.

In the upcoming sections, we will see how this collection of data observation
times permeates the marginal representation of individual variables, the
formulation of cross-variable joint distribution, and the factorization of the
log-likelihood, playing a pivital role in coherently modeling the temporal
irregularities at all stages.

##  4 Marginal Learning

The marginal learning block aims to acknowledge heterogeneities and temporal
irregularities in the component data by addressing two drawbacks of the hidden
state-sharing feature in existing models. These drawbacks are present in
vanilla recurrent architectures such as RNN, GRU, and LSTM, as well as models
that build upon these vanilla recurrent architectures and target irregularly
sampled time series such as those discussed in Section 1.1 and 2.2, i.e.,
ODERNN [rubanova2019latent], GRU-ODE-Bayes [de2019gru], and ODELSTM
[lechner2020learning].

Firstly, the statistical heterogeneities in the marginal dynamics may not be
adequately captured when the hidden states are shared among component
variables. The optimization process aims at minimizing the overall prediction
error for all the variables. The shared hidden neurons are not optimized
towards capturing the individual statistical features both statically, such as
the maximum, the minimum, and the mean of each variable
[harutyunyan2019multitask], and dynamically such as the range of serial
dependencies [belletti2018factorized].

Consider a multivariate financial series where one component variable is the
stock price series, and the rest are the price series of its options at
different strikes and maturities. The one-day price change of the APPLE stock,
for example, is usually no more than a few percentages; however, the daily
price swings of its options can be as wide as several times the base prices if
strikes are far out of the money.

Secondly, when temporal irregularities are present, imputing unobserved
variables introduces unnecessary biases if the hidden states are shared
[zhang2021life]. We illustrate this phenomenon in Figure 3. At time $t_{1}$ in
the left figure, only $x^{2}_{t_{1}}$ has an observation while $x^{1}_{t_{1}}$
and $x^{3}_{t_{1}}$ are not observed. If we generate imputed values
$\tilde{x}^{1}_{t_{1}}$ and $\tilde{x}^{3}_{t_{1}}$ and stack them together
with the observed $x^{2}_{t_{1}}$, input data are biased by definition. The
hidden state learned at $t_{1}$ will inherit these biases and carry them into
the future due to its shared nature, even if all three components can be
observed later such as $t_{2}$. See the figure on the right.

Figure 3: The left sub-figure draws a sample path of Asyn-MTS. A subset of
variables have no observations at time $t_{1},t_{3},t_{4}$; The right sub-
figure draws the corresponding shared hidden states. The learned hidden states
are biased because the input data contain imputed, hence biased, values at
time points $t_{1},t_{3},t_{4}$.

Our solution is to assign unique sets of hidden states to component variables
and let the arrival times of observations determine when, how, and which set
of hidden states to update. This solution makes full use of the heterogeneous
and asynchronous characteristics of the data but does not require imputing
constituent variables in the absence of observations to complete the data. We
use Figure 5 and the necessary formulation to demonstrate its working
mechanism.

Figure 5(a) draws a synchronous sample path, and Figure 5(b) draws an
asynchronous sample path. The sample paths are unevenly spaced in both cases,
i.e., the lengths of time intervals $(t_{1},t_{2})$, $(t_{2},t_{3})$, and
$(t_{3},t_{4})$ are all different. However, data points in the synchronous
case Figure 5(a) are aligned at all times $t_{1},t_{2},t_{3},t_{4}$. In
contrast, no observation time has complete observations in the asynchronous
case Figure 5(b). The following discussion is based on the asynchronous case.
The synchronous case becomes a special case once the asynchronous case is
clear.

Let the hidden states of the three component variables be
$\mathsf{h}^{1}_{t}$, $\mathsf{h}^{2}_{t}$, and $\mathsf{h}^{3}_{t}$. Recall
that the vector $\mathsf{t}:=[t_{1},\cdots,t_{K}]$ is the collection of all
time points at which at least one component series has an observation. We
specify equations governing the evolution of the hidden states in two regimes.

The Continuous Updating regime corresponds to continuous time intervals from
$t_{1+}$ to $t_{2-}$, $t_{2+}$ to $t_{3-}$, $t_{3+}$ to $t_{4-}$ in Figure
5(b) when no component variables have observations. In this regime, we use
equation (6) to evolve the marginal hidden states for all variables
continuously. Here, $t_{-}$ denotes the time point right before $t$ and
$t_{+}$ means the time point immediately after $t$.

The Discrete Updating regime corresponds to the set of discrete time points
$t=t_{1},t_{2},t_{3},t_{4}$ at which at least one component variable is
observed. In this case, we only update the marginal hidden states of those
observed and keep the marginal hidden states of the unobserved variables
unchanged per equation (7). For example, at $t_{1}$, only $x^{2}_{t_{1}}$ and
$x^{3}_{t_{1}}$, but not $x^{1}_{t_{1}}$, are observed. We therefore update
$h^{2}_{t_{1+}}=\mathsf{GRU}^{2}(h^{2}_{t_{1-}},x^{2}_{t_{1}})$ and
$h^{3}_{t_{1+}}=\mathsf{GRU}^{3}(h^{3}_{t_{1-}},x^{3}_{t_{1}})$, and leave
$h^{1}_{t_{1}}$ unchanged, meaning $h^{1}_{t_{1+}}=h^{1}_{t_{1-}}$.

In general, the continuous updating regulates the continuous no-observation
periods $t\in[0,T]\backslash\mathsf{t}$ through Neural ODE (equation (6)),
whereas the discrete updating governs the set of discrete times
$t\in\mathsf{t}$ at which observations arrive via GRU (equation (7)).
Together, they cover the entire evolution of the hidden states through time in
$[0,T]$.

The complete formulation is as follows. For any $t\in[0,T]$, let
$\mathsf{h}^{d}(t)\in\mathbb{R}^{S\times 1}$ be the $d^{\text{th}}$ set of
hidden states that is responsible for learning the $d$-th component variable,
$S$ be the dimension of each set of hidden states and
$\mathbf{h}(t)\in\mathbb{R}^{S\times D}$ be the hidden states matrix of all
components,

$\displaystyle\mathbf{h}(t)=[\mathsf{h}^{1}(t);\cdots;\mathsf{h}^{d}(t);\cdots;\mathsf{h}^{D}(t)],\,t\in[0,T].$

Correspondingly, the reset gates $\mathbf{r}^{\\{c,u\\}}(t)$, update gates
$\mathbf{z}^{\\{c,u\\}}(t)$, and candidate update gates
$\tilde{\mathbf{h}}^{\\{c,u\\}}(t)$ in both updating schemes, i.e. highlighted
with the superscripts $c$ and $u$, are all $\mathbb{R}^{S\times D}$ matrices.

Continuous Updating:

$\displaystyle\forall t\in[0,T]\backslash[t_{1},\cdots,t_{K}],$
$\displaystyle\frac{d\mathbf{h}(t)}{dt}=(1-\mathbf{z}^{c}(t))\odot(\mathbf{\tilde{h}}^{c}(t)-\mathbf{h}(t)),\textit{where}$
(6) $\displaystyle\quad\left\\{\begin{aligned}
\mathbf{r}^{c}(t)&=\sigma\left(\mathcal{W}^{c}_{r}\circledast\mathsf{x}(t)+\mathcal{U}^{c}_{r}\circledast\mathbf{h}(t)+\mathbf{b}^{c}_{r}\right),\\\
\mathbf{z}^{c}(t)&=\sigma\left(\mathcal{W}^{c}_{z}\circledast\mathsf{x}(t)+\mathcal{U}^{c}_{z}\circledast\mathbf{h}(t)+\mathbf{b}^{c}_{z}\right),\\\
\mathbf{\tilde{h}}^{c}(t)&=\tanh\left(\mathcal{W}^{c}_{h}\circledast\mathsf{x}(t)+\mathcal{U}^{c}_{h}\circledast(\mathbf{r}^{c}(t)\odot\mathbf{h}(t))+\mathbf{b}^{c}_{h}\right).\end{aligned}\right.$

Discrete Updating:

$\displaystyle\forall t\in[t_{1},\cdots,t_{K}],$
$\displaystyle\mathbf{h}(t_{+}):=\mathbf{GRU}\left(\mathbf{h}(t_{-}),\mathsf{x}_{t}\right)\odot\mathbf{m}_{t}+\mathbf{h}(t_{-})\odot(1-\mathbf{m}_{t}),$
(7) $\displaystyle\textit{where}\,\,\mathbf{m}_{t}\in\mathbb{R}^{S\times
D}\,\,\textit{and}\,\,\mathbf{GRU}\left(\mathbf{h}(t_{-}),\mathsf{x}_{t}\right)\in\mathbb{R}^{S\times
D}\,\textit{are}$ $\displaystyle\,\left\\{\begin{aligned}
&\mathbf{m}_{t}:=[\mathsf{m}^{1}_{t};\cdots;\mathsf{m}^{D}_{t}]^{\prime},\\\
&\mathbf{GRU}\left(\mathbf{h}(t_{-}),\mathsf{x}_{t}\right):=\Big{[}\mathsf{GRU}^{1}\left(\mathsf{h}^{1}(t_{-}),x_{t}^{1}\right);\cdots;\\\
&\mathsf{GRU}^{d}\left(\mathsf{h}^{d}(t_{-}),x_{t}^{d}\right);\cdots;\mathsf{GRU}^{D}\left(\mathsf{h}^{D}(t_{-}),x_{t}^{D}\right)\Big{]},\end{aligned}\right.$
and the update formulas for GRU are $\displaystyle\,\left\\{\begin{aligned}
&\mathbf{r}^{u}(t_{-})=\sigma\left(\mathcal{W}^{u}_{r}\circledast\mathsf{x}_{t}+\mathcal{U}^{u}_{r}\circledast\mathbf{h}(t_{-})+\mathbf{b}^{u}_{r}\right),\\\
&\mathbf{z}^{u}(t_{-})=\sigma\left(\mathcal{W}^{u}_{z}\circledast\mathsf{x}_{t}+\mathcal{U}^{u}_{z}\circledast\mathbf{h}(t_{-})+\mathbf{b}^{u}_{z}\right),\\\
&\mathbf{\tilde{h}}^{u}(t_{-})=\tanh\left(\mathcal{W}^{u}_{h}\circledast\mathsf{x}_{t}+\mathcal{U}^{u}_{h}\circledast(\mathbf{r}^{u}(t_{-})\odot\mathbf{h}(t_{-}))+\mathbf{b}^{u}_{h}\right),\\\
&\mathbf{h}(t_{+})=(1-\mathbf{z}^{u}(t_{-}))\odot\mathbf{\tilde{h}}^{u}(t_{-})+\mathbf{z}^{u}(t_{-})\odot\mathbf{h}(t_{-}).\end{aligned}\right.$

Figure 4: Schematic visualization of the two $\circledast$-related
calculations in (8).

During the period from right after the $k^{\text{th}}$ observation time to
just before the next one, i.e., $[t_{k_{+}},t_{{k+1}_{-}}]$, all component-
wise hidden states $\mathsf{h}^{d}(t)$ evolve in the same fashion according to
equation (6). However, when the physical time $t$ passes through any of the
observation times $t_{k}$, we check additionally which component variables
have observations and only update those using equation (7) via the standard
GRUs,

$\displaystyle\forall t\in[t_{1},\cdots,t_{K}],\,\forall d=1,2,\cdots,D$
$\displaystyle\mathsf{h}^{d}(t_{+})=\begin{cases}\begin{aligned}
&\mathsf{GRU}^{d}\left(\mathsf{h}^{d}(t_{-}),x_{t}^{d}\right),&\text{if}\
m_{t}^{d}=1;\\\ &\mathsf{h}^{d}\left(t_{-}\right),&\text{if}\
m_{t}^{d}=0.\end{aligned}\end{cases}$

The set of parameters to be trained are

$\zeta^{\\{c,u\\}}=\left\\{\mathcal{W}^{\\{c,u\\}}_{\\{r,z,h\\}},\mathcal{U}^{\\{c,u\\}}_{\\{r,z,h\\}},\mathbf{b}^{\\{c,u\\}}_{\\{r,z,h\\}}\right\\},$

where $\mathbf{b}^{\\{c,u\\}}_{\\{r,z,h\\}}$ are matrices in
$\mathbb{R}^{S\times D}$. The parameters
$\mathcal{W}^{\\{c,u\\}}_{\\{r,z,h\\}}$ and
$\mathcal{U}^{\\{c,u\\}}_{\\{r,z,h\\}}$ are tensors in $\mathbb{R}^{D\times
S\times 1}$ and $\mathbb{R}^{D\times S\times S}$ respectively

$\displaystyle\mathcal{W}$
$\displaystyle=\\{\mathbf{w}^{1};\cdots;\mathbf{w}^{D}\\},\,\mathbf{w}^{d}\in\mathbb{R}^{S\times
1}$ $\displaystyle\mathcal{U}$
$\displaystyle=\\{\mathbf{u}^{1};\cdots;\mathbf{u}^{D}\\},\,\mathbf{u}^{d}\in\mathbb{R}^{S\times
S}.$

The operator $\odot$ is the usual element-wise multiplication. The operator
$\circledast$ is defined as follows, resulting in $\mathbb{R}^{S\times D}$
matrices,

$\displaystyle\left.\begin{aligned}
&\mathcal{W}\circledast\mathsf{x}(t):=[\mathbf{w}^{1}x^{1}(t);\cdots;\mathbf{w}^{D}x^{D}(t)],\\\
&\mathcal{U}\circledast\mathbf{h}(t):=[\mathbf{u}^{1}\mathsf{h}^{1}(t);\cdots;\mathbf{u}^{D}\mathsf{h}^{D}(t)],\end{aligned}\right.$
(8)

whose columns $\mathbf{w}^{d}x^{d}(t)$ and $\mathbf{u}^{d}\mathsf{h}^{d}(t)$
are in $\mathbb{R}^{S\times 1}$ for all $d$ $\in\\{1,\cdots,D\\}$.

Figure 4 visualizes the calculations in (8) schematically. Although each
variable has distinct GRU parameters, there is no need to train them
individually and separately, as it would be time-consuming and impractical for
high-dimensional data. The tensor operations in (8), combined with the matrix
formulations in (6-7), facilitate simultaneous training of all parameters
together.

The formulations we provide so far use GRU and Neural ODE. However, the idea
of using unique sets of hidden states and letting the arrival times of
observations determine their updating is broadly applicable. We demonstrate
its generality by providing formulations when the underlying is replaced by
GRU-D, ODERNN, and ODELSTM in Appendix A.

(a) Synchronous case

(b) Asynchronous case

Figure 5: The framework of RFNs for (a) Syn-MTS and (b) Asyn-MTS. Both (a) and
(b) have three component variables $X^{1}_{t},X^{2}_{t},X^{3}_{t}$. The solid
dots are observations. Different colors refer to different variables. In the
marginal block, each variable has its own set of hidden states
$\mathsf{h}^{1}(t),\mathsf{h}^{2}(t),\mathsf{h}^{3}(t)$. The hidden states of
a particular variable will be updated only when it has an observation. In the
multivariate block, the time-$t_{k}$ base distribution parameters,
$\mu_{t_{k}}$ and $\Sigma_{t_{k}}$, are functions of the hidden states
$\mathbf{h}_{{t_{k}}_{-}}=[\mathsf{h}^{1}_{{t_{k}}_{-}};\mathsf{h}^{2}_{{t_{k}}_{-}};\mathsf{h}^{3}_{{t_{k}}_{-}}]$
for the Syn-MTS case and additionally the component variables
$\mathsf{X}_{t_{k}}^{-d}:=[x^{1}_{t_{k}},\cdots,x^{d-1}_{t_{k}}]$ for the
Asyn-MTS case.

##  5 Multivariate Synchronous Learning

Once the marginal learning block is set up, we seek to maximize the log-
likelihood of all observations,

$\displaystyle\max_{\Phi}\log
p(\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{K}};\Phi),$

where $p$ is the joint density to be estimated from the observed data at all
observation times $[t_{1},\cdots,t_{K}]$ and $\Phi$ is the parameter set.

###  5.1 Factorizing Log-likelihood Objective

If the multivariate sequence is synchronous, observing one variable implies
the observation of all other component variables, despite potentially
unevenly-spaced time intervals. Nevertheless, the synchronous property ensures
alignments, enabling the use of the chain rule of probability,

$\displaystyle\quad\,\,p(\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{K}};\Phi)$
$\displaystyle=p(\mathsf{x}_{t_{1}};\phi_{1})\cdotp
p(\mathsf{x}_{t_{2}}|\mathsf{x}_{t_{1}};\phi_{2})\cdotp\,\,\cdots\,\,\cdotp
p(\mathsf{x}_{t_{K}}|\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{K-1}};\phi_{K})$
$\displaystyle=p(\mathsf{x}_{t_{1}};\phi_{1})\cdotp\prod_{k=2}^{K}\,p(\mathsf{x}_{t_{k}}|\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{k-1}};\phi_{k}),$
(9)

which can be used to factorize the maximization of the joint log-likelihood at
all observation times into optimizing conditional log-likelihoods
individually,

$\displaystyle\quad\,\max_{\Phi}\log
p(\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{K}};\Phi)$
$\displaystyle=\max_{\phi_{1}}\log
p(\mathsf{x}_{t_{1}};\phi_{1})+\sum_{k=2}^{K}\,\max_{\phi_{k}}\log
p(\mathsf{x}_{t_{k}}|\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{k-1}};\phi_{k}),$

where $\phi_{k},k=1,\cdots,K$ is the parameter set of the conditional joint
density at time $t_{k}$.

Note that the hidden states contain the exogenous information extracted from
data. In particular, the amount of exogenous information carried by the
historical observations up to $t_{k-1}$ does not increase until the next
observation arrives at $t_{k}$. Therefore, conditioning on the sample path
$\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{k-1}}$ can be faithfully
approximated by conditioning on the hidden state at time $t_{k-}$,

$\displaystyle
p\left(\mathsf{x}_{t_{k}}|\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{k-1}}\right)\simeq
p\left(\mathsf{x}_{t_{k}}|\mathbf{h}_{t_{k-}}\right),$

under mild assumptions, except that the hidden states just before $t_{1}$ do
not carry any exogenous data information as it is evolved from the initial
hidden state $\mathbf{h}_{t_{0}}$ by the Neural ODE. Thus, the time-$t_{1}$
conditional density $p(\mathsf{x}_{t_{1}}|\mathbf{h}_{t_{1-}};\phi_{1})$
degenerates to its unconditional counterpart $p(\mathsf{x}_{t_{1}};\phi_{1})$.
We can therefore maximize the hidden-states-conditioned log-likelihood $\log
p(\mathsf{x}_{t_{k}}|\mathbf{h}_{t_{k-}};\phi_{k})$ instead,

$\displaystyle\quad\,\max_{\Phi}\log
p(\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{K}};\Phi)$
$\displaystyle=\sum_{k=1}^{K}\max_{\phi_{k}}\log
p\left(\mathsf{x}_{t_{k}}|\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{k-1}};\phi_{k}\right)$
$\displaystyle\simeq\sum_{k=1}^{K}\max_{\phi_{k}}\log
p(\mathsf{x}_{t_{k}}|\mathbf{h}_{t_{k-}};\phi_{k}).$ (10)

###  5.2 Conditional Flow Representation

Optimizing (10) requires a conditional representation of the data
distribution, wherein the log-densities can be tailored to suit one’s specific
conditional choices. For this purpose, we develop the following conditional
formulation of the continuous normalizing flow.

Given two random variables $\mathsf{X}$ and $\mathsf{Y}$, it represents the
conditional log-density $\log p(\mathsf{x}|\mathsf{y})$ of $\mathsf{X}$
through the non-parametric map $\mathsf{f}$ under the control of $\mathsf{Y}$.
We shall call this representation the Conditional CNF.

######  Lemma 1.

Let $\tilde{\mathsf{z}}(s)=[\mathsf{z}(s),\mathsf{y}(s)]^{T}$ be a finite
continuous random variable, and the probability density function of
$\tilde{\mathsf{z}}(s)$ is
$p(\tilde{\mathsf{z}}(s))=p(\mathsf{z}(s),\mathsf{y}(s))$ which depends on
flow time $s$, where $s_{0}\leq s\leq s_{1}$. Given the governing dynamics of
$\tilde{\mathsf{z}}(s)$

$\displaystyle\frac{\partial\tilde{\mathsf{z}}(s)}{\partial
s}=\left[\begin{array}[]{c}\frac{\partial\mathsf{z}(s)}{\partial s}\\\
\frac{\partial\mathsf{y}(s)}{\partial
s}\end{array}\right]=\left[\begin{array}[]{c}\mathsf{f}(\mathsf{z}(s),s,\mathsf{y}(s);\theta)\\\
0\end{array}\right],\text{ where }s_{0}\leq s\leq s_{1},$
$\displaystyle\mathsf{z}(s_{1})=\mathsf{x},\quad\mathsf{y}(s_{1})=\mathsf{y},$

and $\mathsf{f}$ is Lipschitz continuous in $\mathsf{z}$ and continuous in $s$
for any $\mathsf{y}$. We have

$\displaystyle\log p(\mathsf{x},\mathsf{y})=\log
p(\mathsf{z}(s_{0}),\mathsf{y})-\int_{s_{0}}^{s_{1}}\operatorname{Tr}[\partial_{\mathsf{z}(s)}\mathsf{f}(\mathsf{z}(s),s,\mathsf{y};\theta)]ds.$
(11)

Proof: See the Appendix B.

######  Proposition 5.1.

Let the assumptions in Lemma 1 hold. The conditional log-density
$p(\mathsf{x}|\mathsf{y})$ is given by

$\displaystyle\log p(\mathsf{x}|\mathsf{y})=\log
p(\mathsf{z}(s_{0})|\mathsf{y})-\int_{s_{0}}^{s_{1}}\operatorname{Tr}[\partial_{\mathsf{z}(s)}\mathsf{f}(\mathsf{z}(s),s,\mathsf{y};\theta)]ds.$
(12)

Proof: Subtracting $\log p(\mathsf{y})$ on both sides of equation (11) in
Lemma 1 and using the fact that $\log p(\mathsf{x},\mathsf{y})-\log
p(\mathsf{y})=\log\frac{p(\mathsf{x},\mathsf{y})}{p(\mathsf{y})}=\log
p(\mathsf{x}|\mathsf{y})$ gives (12).

Proposition 5.1 overcomes the limitations of existing distribution models,
which are either non-parametric111By non-parametric, we meant the specific
functional form of a distribution is not explicitly known. With this
understanding, the data distribution represented by a neural network
parameterized flow map is considered non-parametric despite neural nets being
technically parameterized functions that are nonlinear. but static, such as
models for images, or time-varying but parametric, such as models in
[salinas2020deepar, de2019gru, salinas2019high] in which data distributions
are assumed Gaussian path-dependent parameters. By extending the unconditional
CNF formulation (see Section 2.3) to the conditional formulation, (12) allows
the data distribution to be both non-parametric via the flow map $\mathsf{f}$
and time-varying, e.g., by choosing the conditional information $\mathsf{y}$
properly.

###  5.3 Time-varying Specification

Our choice of the conditional information $\mathsf{y}$ is the hidden states.
We shall let the hidden states determine the parameter values of the Gaussian
base distribution. Because hidden states contain information of historical
data and vary at different observation times, this choice of conditioning
leads to path-dependent parameters of the base distribution. Consequently,
upon transformation via the flow map, the data distribution also becomes path
dependent.

Specifically, we use the standard Multilayer Perceptron (MLP) to map the
time-${t_{k-}}$ hidden states $\mathbf{h}_{t_{k-}}$ learned from the marginal
learning block to the mean vector $\mu_{t_{k}}$ and covariance matrix
$\Sigma_{t_{k}}$ of the Gaussian base

$\displaystyle\mathsf{Z}_{t_{k}}|\mathbf{h}_{t_{k-}}\sim
N(\mu_{t_{k}},\Sigma_{t_{k}}),\,\text{where}$ (13)
$\displaystyle\\{\mu_{t_{k}},\Sigma_{t_{k}}\\}=\textbf{MLP}(\mathbf{h}_{t_{k-}};\psi),$

resulting in the following conditional base density

$\displaystyle
p(\mathsf{z}_{t_{k}}|\mathbf{h}_{t_{k-}};\mu_{t_{k}},\Sigma_{t_{k}})=\frac{\exp\Big{[}-\frac{1}{2}\big{(}\mathsf{z}_{t_{k}}-\mu_{t_{k}}\big{)}^{\prime}\Sigma^{-1}_{t_{k}}\big{(}\mathsf{z}_{t_{k}}-\mu_{t_{k}}\big{)}\Big{]}}{(2\pi)^{D/2}\text{det}(\Sigma_{t_{k}})^{1/2}}.$

This is in contrast with the unconditional flow models where the base
distribution parameters are constants as it is a Standard Normal
$\mathsf{Z}\sim N(\mathsf{0},\mathbb{I}_{D})$.

For example, at $t_{1}$ in Figure 5(a), the marginal hidden states
$\\{\mathsf{h}^{1}_{t_{1-}},\mathsf{h}^{2}_{t_{1-}},\mathsf{h}^{3}_{t_{1-}}\\}$
determine the $t_{1}$ Gaussian parameters $\\{\mu_{t_{1}},\Sigma_{t_{1}}\\}$
of $\mathsf{Z}_{t_{1}}$. The flow map $\mathsf{f}$ then transform samples from
the Gaussian base density $p(\mathsf{z}_{t_{1}}|\mathbf{h}_{t_{1-}})$ non-
parametrically and concurrently into $t_{1}$ samples of the data density
$p(\mathsf{x}_{t_{1}}|\mathbf{h}_{t_{1-}})$.

It is worth noting that previous studies [salinas2020deepar, de2019gru,
salinas2019high] also use hidden states to learn Gaussian parameters. However,
those Gaussian parameters are associated with the data distribution. Our
representation shifts the Gaussian assumption from the data distribution to
the base distribution, and it is the base distribution parameters that depend
on hidden states. This distinction is crucial because the resulting data
distribution can now be non-parametric, non-Gaussian, and time-varying
simultaneously.

Next, we use the following gated mechanism to incorporate the hidden state
dependence into the flow map

$\displaystyle\mathsf{f}(\mathsf{z}(s),s,\mathsf{h}_{t_{k-}};\theta)=(\mathbf{w}_{z}\mathsf{z}(s)+\mathbf{w}_{h}\mathsf{h}_{t_{k-}}+\mathsf{b}_{z})\sigma(w_{s}s+b_{s}),$
(14)

where $\mathsf{h}_{t_{k-}}\in\mathbb{R}^{H\times 1}$ is reshaped from the
matrix $\mathbf{h}_{t_{k-}}\in\mathbb{R}^{S\times D}$ and $\sigma$ is the
sigmoid activation function. The set of trainable parameters $\theta$ includes
$\mathbf{w}_{h}\in\mathbb{R}^{D\times H}$, $\mathsf{b}_{z}\in\mathbb{R}^{D}$,
and $\\{w_{s},b_{s}\\}\in\mathbb{R}$.

The first term
$(\mathbf{w}_{z}\mathsf{z}(s)+\mathbf{w}_{h}\mathsf{h}_{t_{k-}}+\mathsf{b}_{z})$
transforms the Gaussian samples from the base distribution to samples of the
data-generating distribution subject to the influence of conditional
information $\mathsf{h}_{t_{k-}}$. The second term $\sigma(w_{s}s+b_{s})$ is
the gate that modulates the magnitude multiplicatively of the first term
through the flow time $s$. Together, they determine the instantaneous change
of $\mathsf{z}(s)$

$\displaystyle\frac{\partial\mathsf{z}(s)}{\partial
s}=\mathsf{f}(\mathsf{z}(s),s,\mathsf{h}_{t_{k-}};\theta),s\in[s_{0},s_{1}],$
$\displaystyle\text{where}\quad\mathsf{z}(s)|_{s=s_{0}}=\mathsf{z}_{t_{k}},\,\,\mathsf{z}(s)|_{s=s_{1}}=\mathsf{x}_{t_{k}},$

where the initial value of the flow $\mathsf{z}(s)|_{s=s_{0}}$ is set as the
Gaussian base sample $\mathsf{z}_{t_{k}}$, sampled from
$\mathsf{Z}_{t_{k}}|\mathbf{h}_{t_{k-}}$ per (LABEL:Eq13), and the terminal
flow value $\mathsf{z}(s)|_{s=s_{1}}$ is set to equal the time-$t_{k}$
observed data $\mathsf{x}_{t_{k}}$. Its solution, in integral form, maps
samples from the base distribution to samples of the data distribution
concurrently,

$\displaystyle\mathsf{x}_{t_{k}}=\mathsf{z}_{t_{k}}+\int_{s_{0}}^{s_{1}}\mathsf{f}(\mathsf{z}(s),s,\mathsf{h}_{t_{k-}};\theta)ds,$
(15)

Note that the mapping between $\mathsf{z}_{t_{k}}$ and $\mathsf{x}_{t_{k}}$ in
(15) varies as the physical observation time $t_{k}$ varies.

We summarize the time-varying specification of the conditional CNF
representation for Syn-MTS data as follows. For
$t_{k}\in[t_{1},\cdots,t_{K}]$,

$\displaystyle\quad\,\begin{bmatrix}\mathsf{x}_{t_{k}}\\\ \log
p(\mathsf{x}_{t_{k}}|\mathbf{h}_{t_{k-}};\phi_{t_{k}})\end{bmatrix}=\begin{bmatrix}\mathsf{z}_{t_{k}}\\\
\log
p(\mathsf{z}_{t_{k}}|\mathbf{h}_{t_{k-}};\mu_{t_{k}},\Sigma_{t_{k}})\end{bmatrix}+\int_{s_{0}}^{s_{1}}\begin{bmatrix}\mathsf{f}\big{(}\mathsf{z}(s),s,\mathsf{h}_{t_{k-}};\theta\big{)}\\\
-\operatorname{Tr}[\partial_{\mathsf{z}(s)}\mathsf{f}]\end{bmatrix}ds.$ (16)

It has two essential components, the time-varying base distribution driven by
hidden states where $\psi$ in (LABEL:Eq13) is the parameter to be trained, and
the time-varying flow map indexed by observation arrival times with $\theta$
in (14) being the set of trainable parameters.

###  5.4 Training and Sampling

To train the FRN model, we maximize the conditional log-likelihoods (10)
across all sample instances at all observation times. We set the dataset level
log-likelihood $\mathcal{L_{\text{Syn-MTS}}}$ as the sample average of the $N$
sample instances $\mathcal{L}^{i}_{\text{Syn-MTS}}$, $i=1,\cdots,N$, weighted
by the number of observation times $K_{i}$ of each instance,

$\displaystyle\quad\,\mathcal{L_{\text{Syn-MTS}}}$
$\displaystyle=\frac{1}{N}\underset{i=1}{\overset{N}{\sum}}\mathcal{L}^{i}_{\text{Syn-
MTS}}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{K_{i}}\underset{k=1}{\overset{K_{i}}{\sum}}\log
p(\mathsf{x}_{i,t_{k}^{i}}|\mathbf{h}_{i,t_{k-}^{i}};\phi_{i,t_{k}^{i}})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{K_{i}}\underset{k=1}{\overset{K_{i}}{\sum}}\begin{pmatrix}\log
p(\mathsf{z}_{i,t_{k}^{i}}|\mathbf{h}_{i,t_{k-}^{i}};\mu_{i,t_{k}^{i}},\Sigma_{i,t_{k}^{i}})+\qquad\\\
\quad\int_{s_{1}}^{s_{0}}\operatorname{Tr}\big{[}\partial_{\mathsf{z}(s)}\mathsf{f}(\mathsf{z}(s),s,\mathbf{h}_{i,t_{k-}^{i}};\theta)\big{]}ds\end{pmatrix}.$
(17)

Once the model is trained, we can forecast the joint data distribution at time
point $t_{k}$ given the observations
$\\{\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{k-1}}\\}$ of the $i^{\text{th}}$
instance. The confidence interval of predictions can be obtained by sampling
from the learned joint distribution.

For sampling, we first obtain the hidden states $\mathbf{h}_{t_{k-}}$ for the
$i^{\mathrm{th}}$ instance based on the observations
$\\{\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{k-1}}\\}$ per (6) and (7), and
use $\mathbf{h}_{t_{k-}}$ to predict the base distribution parameters
$\\{\mu_{t_{k}},\Sigma_{t_{k}}\\}$ per (LABEL:Eq13). We then sample a given
number of points $\mathsf{z}_{t_{k}}$ from the base distribution with the
predicted parameters. Finally, we transform these base sample points using the
Conditional CNF map (15) to obtain concurrent samples that follow the data
distribution.

##  6 Multivariate Asynchronous Learning

If the multivariate sequence is asynchronous, it implies the possibility that
not all components of the vector $\mathsf{x}_{t_{k}}$ are observable at any
given time $t_{k}\in\mathsf{t}$ (as defined in Definition 2). Note that the
asynchrony aspect of temporal irregularities arises exclusively in the
multivariate context, whereas the uneven spacing aspect originates essentially
at the univariate level.

###  6.1 Additional Factorization of Log-likelihoods

To address the additional complication associated with asynchrony, we apply
the chain rule of probability again to (9) along the component dimension to
further factorize each conditional log-densities,

$\displaystyle\log p(\mathsf{x}_{t_{k}}|\mathbf{h}_{t_{k-}};\phi_{k})$
$\displaystyle=\log\left(\prod_{d=1}^{D}\,p\big{(}x_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},x_{t_{k}}^{1},\cdots,x_{t_{k}}^{d-1};\phi_{t_{k}}^{d}\big{)}\right)$
$\displaystyle=\sum_{d=1}^{D}\log
p\big{(}x_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},x_{t_{k}}^{1},\cdots,x_{t_{k}}^{d-1};\phi_{t_{k}}^{d}\big{)},$

and define $\mathsf{x}^{-d}_{t_{k}}$ as follows to lighten the notation

$\mathsf{x}^{-d}_{t_{k}}:=[x^{1}_{t_{k}},\cdots,x^{d-1}_{t_{k}}]^{\prime}\in\mathbb{R}^{(d-1)\times
1}$

The maximization of the joint log-likelihood at all observation times then
becomes

$\displaystyle\max_{\Phi}\log
p(\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{K}};\Phi)$
$\displaystyle=\sum_{k=1}^{K}\max_{\\{\phi_{k}^{1},\cdots,\phi_{k}^{D}\\}}\bigg{(}\log\prod_{d=1}^{D}p\big{(}x_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},\mathsf{x}^{-d}_{t_{k}};\phi_{t_{k}}^{d}\big{)}\bigg{)}$
$\displaystyle=\sum_{k=1}^{K}\max_{\\{\phi_{k}^{1},\cdots,\phi_{k}^{D}\\}}\bigg{(}\sum_{d=1}^{D}\log
p\big{(}x_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},x_{t_{k}}^{1},\cdots,x_{t_{k}}^{d-1};\phi_{t_{k}}^{d}\big{)}\bigg{)}$
$\displaystyle=\sum_{k=1}^{K}\sum_{d=1}^{D}\max_{\phi_{k}^{d}}\log
p\big{(}x_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},\mathsf{x}^{-d}_{t_{k}};\phi_{t_{k}}^{d}\big{)}.$
(18)

Factorization along the component dimension enables optimizing only the
conditional distribution of variables with observations at a given observation
point $t_{k}$, rather than optimizing for all $D$ variables.

Note that different from the Syn-MTS case (equation (10)), the conditional
information in (18) contains not only the historical information embedded in
the hidden states $\mathbf{h}_{t_{k-}}$, but also the concurrent information
from the cross-variables $\mathsf{x}^{-d}_{t_{k}}$ at $t_{k}$. Below we show
how to embed the concurrent information into the base distribution and the
flow map.

###  6.2 Conditional Information Specification

Firstly, for the base distribution, unlike the Syn-MTS case that uses the
standard fully-connected MLP to map hidden states $\mathbf{h}_{t_{k-}}$ to the
distribution parameters $\mu_{t_{k}}$ and $\Sigma_{t_{k}}$, we shall design a
masked-MLP (m-MLP) that keeps only the first $d$$-$$1$ variables
$[x^{1}_{t_{k}},\cdots,x^{d-1}_{t_{k}}]$ by muting the rest
$[x^{d+1}_{t_{k}},\cdots,x^{D}_{t_{k}}]$,

$\displaystyle\\{\mu^{d}_{t_{k}},\Sigma^{d}_{t_{k}}\\}=\textbf{m-MLP}(\mathbf{h}_{t_{k-}},\mathsf{x}_{t_{k}}^{-d};\psi_{m}),$
(19)

where $\mu_{t_{k}}^{d}\in\mathbb{R}^{1}$ and
$\Sigma_{t_{k}}^{d}\in\mathbb{R}^{1\times 1}$. Figure 6 draws the architecture
of the masked-MLP, where the inputs are the hidden states and all variables,
and the outputs are the parameters of the base distribution with different
conditioned information.

Secondly, for the flow model, we reshape the hidden state matrix
$\mathbf{h}_{t_{k-}}\in\mathbb{R}^{S\times D}$ into vector
$\mathsf{h}_{t_{k-}}\in\mathbb{R}^{H\times 1}$, and encode it into flow model,
similar to treatment in the Syn-MTS case. In addition, in Asyn-MTS, we also
need to fuse concurrent variables $\mathsf{x}^{-d}_{t_{k}}$ into the flow map,

$\displaystyle\quad\,\,\mathsf{f}(z^{d}(s),s,\mathsf{h}_{t_{k-}},\mathsf{x}^{-d}_{t_{k}};\theta)=(w_{z}z^{d}(s)+\mathsf{w}_{h}\mathsf{h}_{t_{k-}}+\mathsf{w}_{x}\mathsf{x}^{-d}_{t_{k}}+b_{z})\sigma(w_{s}s+b_{s}),$
(20)

where $\\{w_{z},b_{z},w_{s},b_{s}\\}\in\mathbb{R}$,
$\mathsf{w}_{x}\in\mathbb{R}^{1\times(d-1)}$,
$\mathsf{w}_{h}\in\mathbb{R}^{1\times H}$ are trainable parameters.

Figure 6: The architecture of the masked-MLP (m-MLP). The Bottom two squares
are weight matrices for the first layer and second layer, in which the shaded
squares represent the trainable parameters and the unshaded squares represent
the parameters set to zero in the training process.

Note that the observation time $t_{k}$ is defined as a time whenever
constituent variables have observations. It does not guarantee every
constituent variable has full observations. In this case, we use the marginal
predictions

$\displaystyle\tilde{\mathsf{x}}_{t_{k}}=\mathsf{x}_{t_{k}}\odot\mathsf{m}_{t_{k}}+\hat{\mathsf{x}}_{t_{k}}\odot(1-\mathsf{m}_{t_{k}})$
$\displaystyle\text{where}\quad\hat{\mathsf{x}}_{t_{k}}=g(\mathbf{h}_{t_{k-}};\eta),$

to address the issue that certain variables in the vector
$\mathsf{x}^{-d}_{t_{k}}$ are still "missing" even if we only optimize the
conditional distribution of the observed variables at this time. Here
$g:\mathbb{R}^{S\times D}\rightarrow\mathbb{R}^{D}$ is a function that maps
the hidden states to values in data series, with $\eta$ being the trainable
parameter. Equation (18) therefore becomes

$\displaystyle\quad\,\max_{\Phi}\log
p(\mathsf{x}_{t_{1}},\cdots,\mathsf{x}_{t_{K}};\Phi)=\sum_{k=1}^{K}\sum_{d=1}^{D}\max_{\phi_{t}^{d}}\log
p(x_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},\tilde{\mathsf{x}}_{t_{k}}^{-d};\phi_{t_{k}}^{d}),$
(21)

where
$\tilde{\mathsf{x}}_{t_{k}}^{-d}:=[\tilde{x}_{t_{k}}^{1},\cdots,\tilde{x}_{t_{k}}^{d-1}]^{\prime}$.

We summarize the conditional CNF representation for Asyn-MTS data as

$\displaystyle\quad\,\begin{bmatrix}x^{d}_{t_{k}}\\\ \log
p(x_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},\tilde{\mathsf{x}}_{t_{k}}^{-d};\phi_{t_{k}}^{d})\end{bmatrix}=\begin{bmatrix}z^{d}_{t_{k}}\\\
\log
p\big{(}z_{t_{k}}^{d}|\mathbf{h}_{t_{k-}},\tilde{\mathsf{x}}_{t_{k}}^{-d};\mu^{d}_{t_{k}},\Sigma^{d}_{t_{k}}\big{)}\end{bmatrix}+\int_{s_{0}}^{s_{1}}\begin{bmatrix}\mathsf{f}(z^{d}(s),s,\mathsf{h}_{t_{k-}},\tilde{\mathsf{x}}^{-d}_{t_{k}};\theta)\\\
-\operatorname{Tr}\big{[}\partial_{z^{d}(s)}\mathsf{f}\big{]}\end{bmatrix}ds$
(22)

where the scalar sample $z^{d}_{t_{k}}$ is a realization of

$Z^{d}_{t_{k}}|\mathbf{h}_{t_{k-}},\tilde{\mathsf{x}}_{t}^{-d}\;\sim\;N(\mu^{d}_{t_{k}},\Sigma^{d}_{t_{k}}),$

whose parameters $\psi_{m}$ are learned through (19).

###  6.3 Training

There are two tasks in training the Asyn-MTS model, maximizing the log-
likelihood of observations and imputing the missing values accurately.
Therefore, different from the objective function (equation (18)) of the Syn-
MTS model, we add the Mean Square Error (MSE) term to the objective function
of the Asyn-MTS model,

$\displaystyle\mathcal{L}^{i}_{\text{Asyn-
MTS}}=\sum_{k=1}^{K_{i}}\sum_{d=1}^{D}\frac{1}{K_{i}D}\bigg{[}(x_{i,t_{k}^{i}}^{d}-\hat{x}_{i,t_{k}^{i}}^{d})^{2}\odot
m_{i,t_{k}^{i}}^{d}+\log
p\big{(}z_{i,t_{k}^{i}}^{d}|\mathbf{h}_{i,t_{k-}^{i}},\tilde{\mathsf{x}}^{-d}_{i,t_{k}^{i}};\mu_{i,t_{k}^{i}}^{d},\Sigma_{i,t_{k}^{i}}^{d}\big{)}+\int_{s_{1}}^{s_{0}}\operatorname{Tr}[\partial_{z^{d}(s)}\mathsf{f}]ds\bigg{]}.$
(23)

Here, we only compute the MSE associated with the variables that have
observations since we have no ground truth of the unobserved one. The dataset
level log-likelihood is then the sum of $\mathcal{L}^{i}_{\text{Asyn-MTS}}$
for all sample instances.

###  6.4 Sampling

The sampling process at time $t_{k}$ of the Asyn-MTS model differs from that
of the Syn-MTS model and is similar to the Gibbs sampling. We first predict
the parameters of base distribution of variable $X_{t_{k}}^{1}$ conditional on
the hidden state $\mathbf{h}_{t_{k-}}$ from marginal learning block, i.e.,
$\\{\mu_{t_{k}}^{1},\Sigma_{t_{k}}^{1}\\}$
$=\textbf{m-MLP}(\mathbf{h}_{t_{k-}};\psi_{m})$. Then, we sample a given
number of points $z_{t_{k}}^{1}$ from this base distribution, and we transform
these points using the flow model. Next, we predict the parameters of base
distribution of variable $X_{t_{k}}^{2}$ conditional on the hidden state
$\mathbf{h}_{t_{k-}}$ and the generated samples of $X_{t_{k}}^{1}$, i.e.,
$\\{\mu_{t_{k}}^{2},\Sigma_{t_{k}}^{2}\\}$
$=\textbf{m-MLP}(\mathbf{h}_{t_{k-}},x^{1}_{t_{k}};\psi_{m})$. Then, we obtain
samples $z_{t_{k}}^{2}$ from the second base distribution and transform these
points using the flow model. The process continues until the last variable is
sampled.

We verify the components proposed above by further ablation studies in section
7.3. The algorithms for the training and sampling procedures of the Syn-MTS
model and the Asyn-MTS model are provided in Appendix D.

##  7 Experiments

We conduct three ablation studies to assess the efficacy of the model
specification, followed by three experiments using real-world datasets to
evaluate the overall performance of the RFNs. In all studies and experiments,
we use 70% of the sample instances for training, 15% for validation, and 15%
for testing.

All studies and experiments are run on Dell Precision 7920 Workstations with
Intel(R) Xeon(R) Gold 6256 CPU at 3.60GHz, and three sets of NVIDIA Quadro
GV100 GPUs. All models are implemented in Python 3.8. The versions of the main
packages of our code are: Pytorch 1.8.1+cu102, torchdiffeq: 0.2.2, Sklearn:
0.23.2, Numpy: 1.19.2, Pandas: 1.1.3, Matplotlib: 3.3.2.

###  7.1 Baseline Models

We set up four baseline models, GRU-D, ODERNN, GRUODE, and ODELSTM, with and
without an RFN specification. In their vanilla forms, these baselines all aim
to model multivariate irregularly sampled time series but assume observations
follow a multivariate normal distribution and predict parameters of the data
distribution based on the hidden state. To compare the vanilla form of a
baseline model, e.g., GRU-D, with its RFN counterpart, e.g., termed as RFN-
GRU-D, we use the corresponding RFN specifications laid out in Appendix A as
the marginal learning block, for instance, (A.1) for GRU-D, (A.2) for ODERNN,
(A.3) for ODELSTM, and use the formulations in Section 5 and 6 for the
multivariate learning. For GRUODE, its RFN specification has been discussed in
length in the main body of the paper (see (6) and (7)). Hyperparameters, such
as learning rate and batch size, are tuned using a random search for all
experiments. The ’dopri5’ method is employed as the numerical solver for the
ODE-based models. The parameters of the best model obtained during training
are saved for testing purposes.

###  7.2 Evaluation Metric

We use the Continuous Ranked Probability Score (CRPS) [matheson1976scoring] as
the evaluation metric to measure the proximity between the learned
distribution and the empirical distribution of data. The CRPS metric is the
suitable scoring function for the evaluation since it attains the minimum when
the learned distribution of data aligns with its empirical distribution. The
CRPS is defined as

$\displaystyle\operatorname{CRPS}(\hat{F},x)=\int_{\mathbb{R}}(\hat{F}(z)-\mathbb{I}\\{x\leq
z\\})^{2}\mathrm{\leavevmode\nobreak\ d}z,$

where $\hat{F}$ is the estimated CDF of random variable $X$ and $x$ is the
realized observations of $X$. In our dataset, the observations are irregularly
sampled, and some values are missing, and thus we only compute $\mathrm{CRPS}$
when the variables can be observed.

Since $\mathrm{CRPS}$ targets to evaluate the estimated distribution of the
univariate random variable, [salinas2019high] suggests using
$\mathrm{CRPS}_{\text{sum }}$ as a new proper multivariate-scoring rule to
evaluate the estimated distribution of multivariate random variables,

$\displaystyle\mathrm{CRPS}_{\text{sum
}}=\mathbb{E}_{t}\biggl{[}\operatorname{CRPS}\biggl{(}\hat{F}_{\text{sum
}}(t),\sum_{d=1}^{D}x_{t}^{d}\biggl{)}\biggl{]},$

where $\hat{F}_{\text{sum }}$ is computed by adding the estimated CDF of all
variables together. In practice, the empirical CDF is always used to represent
$\hat{F}$ and $\hat{F}_{\text{sum }}$. In our experiments, we take 100 samples
from trained models to estimate the empirical CDF. For Asyn-MTS datasets, some
variables are not observable at an observation time; we thus only add the
estimated CDF of observed variables as $\hat{F}_{\text{sum }}$.

###  7.3 Ablation Studies via Simulation

(a) Data Generating Process. We set the Data Generation Process (DGP)
$\mathsf{X_{t}}=[X^{1}_{t},\cdots,X^{5}_{t}]$ as the following continuous-time
correlated Geometric Brownian Motion (GBM) [glasserman2004monte] to generate
sample instances,

$\displaystyle
d\mathsf{X}_{t}=\mathsf{\mu}\odot\mathsf{X}_{t}dt+\text{diag}(\mathsf{\sigma}\odot\mathsf{X}_{t})d\mathsf{W}_{t},\,t\geq
0,$ (24)

where $\mathsf{\mu}=[\mu^{1},\cdots,\mu^{5}]^{\prime}$ and
$\sigma=[\sigma^{1},\cdots,\sigma^{5}]^{\prime}$ are the drift term and
volatility term of variables, and
$\mathsf{W}_{t}=[W^{1}_{t},\cdots,W^{5}_{t}]^{\prime}$ is the multivariate
Brownian Motion process such that $\langle
dW_{t}^{i},dW^{j}_{t}\rangle=\rho_{ij}(t)dt$ where

$[\rho_{ij}(t)]_{1\leq i,j\leq 5}:=\begin{aligned}
\sin\left(\frac{\pi}{2}t\right)\begin{bmatrix}1&\rho_{1}&0&0&0\\\
\rho_{1}&1&0&0&0\\\ 0&0&1&\rho_{2}&\rho_{2}\\\ 0&0&\rho_{2}&1&\rho_{2}\\\
0&0&\rho_{2}&\rho_{2}&1\\\ \end{bmatrix}.\end{aligned}$

We devise such a correlated GBM as the DGP for several reasons. Firstly, the
GBM is a continuous-time stochastic process where the instantaneous changes of
state variables are explicitly governed by stochastic differential equations
(SDEs). It is a fundamental model in describing diffusive phenomena in physics
and random fluctuations of asset prices in financial markets, for example, in
the Black-Scholes model [black1973pricing]. Secondly, the correlated GBM
admits log-normal distribution for each variable, distinguishing it from the
normal distribution. The flexibility in setting the drift, volatility, and
correlation terms allows one to generate diverse statistical properties of the
sample paths. Thirdly, the correlation among the five variables is not set to
static but time-varying, gradually increasing as time progresses from $0$ to
$1$ due to the sinusoidal function $\sin(\frac{\pi}{2}t)$. This guarantees
that the joint distribution of the five variables is constantly changing.
Figure 7 showcases three representative sample paths from the dataset.
Additionally, Figure 7 illustrates the sample correlation matrix at observed
time points 0.3, 0.6, and 0.9, validating the dynamic and gradual increase in
correlation.

Figure 7:  (a) Three representative sample paths. (b) The sample correlation
matrix at observed time points 0.3, 0.6, 0.9.

Figure 8: The sample correlation matrix recovered from the learned Syn-MTS
model (a) and the learned Asyn-MTS model (b) at 0.3, 0.6, and 0.9.

(a) Syn-MTS

(b) Asyn-MTS

(c) Zoom-in of Syn-MTS

(d) Zoom-in of Asyn-MTS

Figure 9: Predicted quantile intervals (shaded area) from 20% to 80% of Syn-
MTS data (a,c) and Asyn-MTS data (b,d) using the RFN-GRUODE specification.
Solid dots represent observations. Broader quantile ranges are available in
Appendix C.

(b) Simulation Settings. We first simulate 1,000 sample paths using (24) with
uniform time intervals. The Brownian parameters, i.e., the drift terms and the
diffusion terms, of the 1,000 sample paths are set as

$\displaystyle\mu^{1}=\mu^{2}\sim\text{Unif}(-0.2,-0.05),$
$\displaystyle\mu^{3}=\mu^{4}=\mu^{5}\sim\text{Unif}(0.05,0.2),$
$\displaystyle\sigma^{1}=\sigma^{2}\sim\text{Unif}(0.15,0.3),$
$\displaystyle\sigma^{3}=\sigma^{4}=\sigma^{5}\sim\text{Unif}(0.15,0.3),$

where $\text{Unif}(\cdot,\cdot)$ stands for a uniform distribution.

We create the Syn-MTS and Asyn-MTS datasets through random sampling from the
simulated dataset with complete observation. In the Syn-MTS dataset, we
randomly choose half of the time points and sample the corresponding
observations $\mathsf{X}_{t}$ at the chosen time points. To generate the Asyn-
MTS dataset, we randomly eliminate half of the observations for each variable.
The elimination is independent between any two variables. Therefore, some
variables have observations at an observed time point, and some variables do
not.

Table 1: Geometric Brownian Motions: Baseline models in vanilla form (column
2,3) vs. using RFN specification (column 5,6). Numbers underneath
$\mathrm{CRPS},\mathrm{CRPS}_{\text{sum }}$ are displayed in the form of mean
$\pm$ std, averaged over five experiments. The last two columns
$\mathrm{CRPS}$$\uparrow$ and $\mathrm{CRPS}_{\text{sum }}$$\uparrow$ report
percentage increases of performance using RFN specification.

model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ | model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ |  $\mathrm{CRPS}$$\uparrow$ |  $\mathrm{CRPS}_{\text{sum }}$$\uparrow$
---|---|---|---|---|---|---|---
Syn-MTS
GRUODE | 0.2428 $\pm$ 0.0135 | 0.9423 $\pm$ 0.0831 | RFN-GRUODE | 0.1947 $\pm$ 0.0048 | 0.6116 $\pm$ 0.0255 | 19.83% | 35.01%
ODELSTM | 0.2556 $\pm$ 0.0101 | 1.0219 $\pm$ 0.0456 | RFN-ODELSTM | 0.1914 $\pm$ 0.0044 | 0.5946 $\pm$ 0.0062 | 25.12% | 41.81%
GRU-D | 0.2227 $\pm$ 0.0148 | 0.8253 $\pm$ 0.0948 | RFN-GRU-D | 0.1903 $\pm$ 0.0029 | 0.5922 $\pm$ 0.0059 | 14.54% | 28.24%
ODERNN | 0.2543 $\pm$ 0.0230 | 1.0025 $\pm$ 0.1439 | RFN-ODERNN | 0.1940 $\pm$ 0.0012 | 0.6077 $\pm$ 0.0074 | 23.72% | 39.38%
Asyn-MTS
GRUODE | 0.2475 $\pm$ 0.0061 | 0.4982 $\pm$ 0.0172 | RFN-GRUODE | 0.1884 $\pm$ 0.0020 | 0.3410 $\pm$ 0.0039 | 23.90% | 31.54%
ODELSTM | 0.3164 $\pm$ 0.0639 | 0.6185 $\pm$ 0.1106 | RFN-ODELSTM | 0.1927 $\pm$ 0.0028 | 0.3472 $\pm$ 0.0021 | 39.09% | 43.86%
GRU-D | 0.2344 $\pm$ 0.0049 | 0.4670 $\pm$ 0.0105 | RFN-GRU-D | 0.1902 $\pm$ 0.0016 | 0.3506 $\pm$ 0.0054 | 18.87% | 24.92%
ODERNN | 0.2315 $\pm$ 0.0057 | 0.4597 $\pm$ 0.0105 | RFN-ODERNN | 0.1909 $\pm$ 0.0021 | 0.3446 $\pm$ 0.0045 | 17.56% | 25.04%

(c) Data and Model Validation. The experiment results in Table 1 demonstrate
that RFNs (e.g., RFN-GRUODE, RFN-ODELSTM, RFN-GRU-D, RFN-ODERNN) outperform
the corresponding baselines with smaller CRPS and $\text{CRPS}_{\text{sum}}$.
This indicates the superior ability of the RFNs to capture the joint
distribution of multivariate irregularly sampled time series at each observed
time point and recover its dynamics. Figure 8 and 8 illustrate the correlation
matrices sampled from the estimated Syn-MTS and Asyn-MTS models at observed
time points 0.3, 0.6, and 0.9, closely resembling the ground truth in Figure
7. The correlation coefficients from the recovered matrices gradually
transition from weak to strong as time progresses, aligning with the ground
truth in Figure 7. This dynamic capture of the joint distribution validates
the effectiveness of the RFN specification. Lastly, we employ the RFNs for
interval estimation of correlated GBM processes. Figure 9(a) and Figure 9(b)
display prediction intervals for Syn-MTS and Asyn-MTS, respectively, while
Figure 9(c) and Figure 9(d) provide a closer view of the predicted Syn-MTS and
Asyn-MTS.

(d) Effects of Sharing Hidden States. We first examine the proposition that
training each component variable independently in the marginal learning block
can mitigate the biases introduced by sharing hidden states, therefore better
capturing component heterogeneities.

The comparison is between training an RFN model using one shared set of hidden
states in the marginal learning block versus assigning unique sets of hidden
states to variables. Results from Figure 11 and 11 indicate that training
variables independently yield consistent and noticeable positive improvements
in the Syn-MTS model. Since the data is synchronous, component variables
differ only in their statistical features, statically and dynamically. This
shows the importance of preserving heterogeneities in the marginal
distributions and their dynamics. Meanwhile, the improvements become much more
significant when data is asynchronous, implying that the biases caused by and
spread from sharing hidden states among variables (see Figure 3 for an
illustration) are almost unbearable.

(a) $\mathrm{CRPS}$ (b) $\mathrm{CRPS}_{\text{sum }}$

Figure 10: Syn-MTS: Training variables independently (left bars in lighter
green) vs. Sharing hidden states (right bars in darker green). The smaller the
values of $\mathrm{CRPS}$ and $\mathrm{CRPS}_{\text{sum }}$ are, the better.

(a) $\mathrm{CRPS}$ (b) $\mathrm{CRPS}_{\text{sum }}$

Figure 11: Asyn-MTS: Training variables independently (left bars in lighter
green) vs. Sharing hidden states (right bars in darker green).

(e) Effects of Imputing Unobserved Variables. When we specify the multivariate
asynchronous learning block in Section 6, we posit that there are two ways to
handle concurrent "missing" values across variables at a given observation
time.

We can either impute the unobserved variables, which completes the
asynchronous dataset into a synchronous one and train them using the Syn-MTS
model, or utilize our Asyn-MTS model, where we factorize the conditional log-
likelihood along the component dimension one more time at any observation
time, and only compute the resulting factorized conditional log-likelihoods
for variables that have observations. In other words, the objective function
(18) avoids computing any loss terms for unobserved variables. The unobserved
variables only serve as concurrent conditional information for computing the
conditional likelihoods of observed variables. Figure 13 indicates clearly
that our approach (the left bars in lighter green), without completing the
data, is indeed robustly superior than the "imputing before training"
approach, where the unobserved variables are imputed by their marginal block
forecasts, despite it being more intuitive.

(a) $\mathrm{CRPS}$ (b) $\mathrm{CRPS}_{\text{sum }}$

Figure 12: Asynchronous data: Imputing unobserved variables (right bars in
darker green) vs Not imputing them (left bars in lighter green).

(a) $\mathrm{CRPS}$ (b) $\mathrm{CRPS}_{\text{sum }}$

Figure 13: Asyn-MTS model: Using concurrent variables (left bars in lighter
green) vs. Without using them (right bars in darker green).

(f) Effects of Conditional Choices. Third, in the Asyn-MTS model, we model the
distribution of observed variables conditional on the historical information
via hidden states $\mathbf{h}_{t_{k}}$ and the concurrent information
$\tilde{\mathsf{x}}_{t_{k}}^{-d}$, i.e., $\log
p(\mathsf{x}_{t_{k}}^{d}|\mathbf{h}_{t_{k}},\tilde{\mathsf{x}}_{t_{k}}^{-d})$
in equation (21). We conduct an ablation study that computes the distribution
only conditional on the past values and drops the concurrent values, i.e.,
$\log p(\mathsf{x}_{t_{k}}^{d}|\mathbf{h}_{t_{k}})$. Figure 13 shows that
improvements from conditional on both past and concurrent values are robust.
The improvement is modest because in this study the hidden state already
captures the dependence structure of all variables from the previous time
step. In other words, if the dependence structure remains relatively stable at
the current time step, the added value of incorporating the new incoming data
is incremental.

###  7.4 Dataset of Physical Activities (MuJoCo)

In applications such as robotics, multivariate time series are commonly used
to capture the position, speed, and trajectory of objects, exhibiting a high
degree of correlation among variables. It is crucial to understand the
dependence structure within these time series. However, irregular and
asynchronous frequency of measurements frequently arises due to limitations in
measurement devices. This experiment tests the performance of using the RFN
specification to predict object positions against baseline models in their
vanilla forms.

The MuJoCo physics dataset is first introduced in [rubanova2019latent] to
verify that the ODE-based models can learn an approximation of Newtonian
physics. The dataset is created by the “Hopper” model from the Deepmind
Control Suite [tassa2018deepmind], and it contains 14 variables in total,
where the first 7 variables and the last 7 variables control the position and
the velocities of the “Hopper”, respectively.

By following the generation and preprocessing of [rubanova2019latent], we
randomly sample the initial position and velocities of “Hopper” such that the
hopper rotates in the air and falls on the ground. We generate 5,000
regularly-sampled time series instances with 150 observed time points for each
instance, and randomly sample 50% values to generate the Syn-MTS and Asyn-MTS
datasets.

Table 2 shows that the RFNs perform robustly better than the baselines in
their vanilla form. As a validation, we use the trained model to forecast the
future joint distribution at t=1.40 based on a given instance’s past ten
observations (ranging from t=0.09 to t=1.17). Subsequently, we generate 10
samples from the estimated distribution at t=1.40 and compare them with the
ground truth at t=1.40. Figure 14 displays the 10 samples at t=1.40 from both
the Asyn-MTS model and the Syn-MTS model. They all resemble the ground truth
closely, affirming the RFN’s proficiency in predicting the joint distribution.

(a) t=0.09 (b) t=0.21 (c) t=0.36 (d) t=0.47

(e) t=0.67 (f) t=0.70 (g) t=0.87 (h) t=0.88

(i) t=1.03 (j) t=1.17 (k) t=1.40

Figure 14: (a)-(j) represent data observed at historical times with unevenly
spaced intervals. The objective is to estimate the data distribution at
t=1.40. Fig. (k) displays the observed ground truth at t=1.40. We train the
Syn-MTS and Asyn-MTS models to forecast the desired distribution at t=1.40.
The first and last five figure in the second row displays five samples
generated from the trained Syn-MTS model (Asyn-MTS model) at t=1.40.  Table 2:
Physical Activities of Human Body (MuJoCo). (Same reporting formats as in
Table 1)

model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ | model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ |  $\mathrm{CRPS}$$\uparrow$ |  $\mathrm{CRPS}_{\text{sum }}$$\uparrow$
---|---|---|---|---|---|---|---
Syn-MTS
GRUODE | 0.2312 $\pm$ 0.0056 | 1.0877 $\pm$ 0.0776 | RFN-GRUODE | 0.2336 $\pm$ 0.0032 | 1.0570 $\pm$ 0.0392 | -1.04% | 2.82%
ODELSTM | 0.2599 $\pm$ 0.0021 | 1.3501 $\pm$ 0.0417 | RFN-ODELSTM | 0.2593 $\pm$ 0.0064 | 1.2163 $\pm$ 0.0563 | 0.24% | 9.91%
GRU-D | 0.2465 $\pm$ 0.0058 | 1.1992 $\pm$ 0.0484 | RFN-GRU-D | 0.2229 $\pm$ 0.0033 | 1.0140 $\pm$ 0.0409 | 9.58% | 15.44%
ODERNN | 0.2258 $\pm$ 0.0077 | 1.0571 $\pm$ 0.0578 | RFN-ODERNN | 0.2204 $\pm$ 0.0059 | 0.9929 $\pm$ 0.0194 | 2.37% | 6.07%
Asyn-MTS
GRUODE | 0.2950 $\pm$ 0.0066 | 1.1519 $\pm$ 0.0212 | RFN-GRUODE | 0.2183 $\pm$ 0.0025 | 0.8382 $\pm$ 0.0099 | 25.98% | 27.24%
ODELSTM | 0.3259 $\pm$ 0.0088 | 1.2946 $\pm$ 0.0479 | RFN-ODELSTM | 0.2403 $\pm$ 0.0065 | 0.8911 $\pm$ 0.0279 | 26.26% | 31.17%
GRU-D | 0.2584 $\pm$ 0.0064 | 1.0170 $\pm$ 0.0352 | RFN-GRU-D | 0.2112 $\pm$ 0.0027 | 0.8155 $\pm$ 0.0019 | 18.27% | 19.81%
ODERNN | 0.2955 $\pm$ 0.0042 | 1.1476 $\pm$ 0.0179 | RFN-ODERNN | 0.2403 $\pm$ 0.0050 | 0.9060 $\pm$ 0.0192 | 18.67% | 21.05%

###  7.5 Dataset of Climate Records (USHCN)

Variables in natural phenomena, such as climate data, also display strong
serial and cross correlations. For instance, temperature values from one
season provide valuable information about temperature patterns in the
following season. Additionally, precipitation, such as rain and snow, can
cause temperature drops and affect humidity levels. Missing observations are
also common in climate data due to inclement weather or equipment
malfunctions. In this experiment, we evaluate the performance of the RFN
specification on the USHCN dataset.

The United States Historical Climatology Network (USHCN) dataset
[menne2010long] consists of daily measurements from 1,218 centers across the
country. It includes 5 variables: precipitation, snowfall, snow depth, maximum
temperature, and minimum temperature. Following the preprocessing approach of
[de2019gru], we select the training data from the first quarter of the last
four years (1996-2000). This yields 4,494 instances, with each instance’s time
period normalized from 0 to 12.5 at intervals of 0.1. The observations for
each variable have uneven spacing, making the dataset asynchronous.
Consequently, we exclusively employ the Asyn-MTS model for training this
dataset.

The results presented in Table 3 demonstrate the superior performance of the
RFN specification over its vanilla counterpart for all baselines when applied
to the climate dataset. These findings reinforce that climate data does not
conform to a multivariate Gaussian distribution and highlight the importance
of capturing the dependence structure among variables for accurate weather
forecasting.

Table 3: Climate Records of Weather (USHCN). (Same reporting formats as in
Table 1)

model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ | model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ |  $\mathrm{CRPS}$$\uparrow$ |  $\mathrm{CRPS}_{\text{sum }}$$\uparrow$
---|---|---|---|---|---|---|---
GRUODE | 0.3129 $\pm$ 0.0024 | 0.5914 $\pm$ 0.0031 | RFN-GRUODE | 0.2605 $\pm$ 0.0153 | 0.5385 $\pm$ 0.0229 | 16.75% | 8.94%
ODELSTM | 0.3113 $\pm$ 0.0049 | 0.5909 $\pm$ 0.0075 | RFN-ODELSTM | 0.2723 $\pm$ 0.0222 | 0.5494 $\pm$ 0.0264 | 12.54% | 7.02%
GRU-D | 0.3048 $\pm$ 0.0030 | 0.5829 $\pm$ 0.0066 | RFN-GRU-D | 0.2636 $\pm$ 0.0133 | 0.5382 $\pm$ 0.0179 | 13.51% | 7.66%
ODERNN | 0.3124 $\pm$ 0.0022 | 0.5917 $\pm$ 0.0039 | RFN-ODERNN | 0.2544 $\pm$ 0.0095 | 0.5276 $\pm$ 0.0139 | 18.57% | 10.84%

Figure 15: Transaction prices and times of the Tencent stock and its options
with moneyness levels 0.96 (in the money), 1.00 (at the money), and 1.10 (out
the money) from 01/Dec/2014 to 31/Dec/2017. Values are "missing" means there
are no transaction records during the 3-minute intervals.

###  7.6 Dataset of Financial Transactions (Stock Options)

Stock options are financial derivatives used for speculation and hedging
purposes, allowing investors to buy or sell company stock at a predetermined
price (the strike) within a specified time period (the maturity). Models in
financial mathematics assume that the underlying stock follows a specific
stochastic differential equation, such as the Heston model [heston1993closed]
or the SABR model [hagan2002managing], but require a concurrent set of option
prices at different strikes to calibrate the parameters. However, options with
different strikes have varying liquidities in reality, with ATM options
trading more frequently than ITM and OTM ones222At-the-money (ATM) options are
options whose strikes are equal to the underlying stock price, whereas in-the-
money (ITM) and out-of-the-money (OTM) refer to cases where strike prices are
below or above the stock price.. This liquidity disparity is illustrated in
Figure 15 using the transaction records of Tencent, the stock with the largest
market capitalization in the HKEX, and its options with moneyness (defined as
the ratio of a strike over by the stock price) of 0.96, 1.00, and 1.10. Figure
15 clearly shows that options with different moneyness do not align, and they
are all unevenly spaced. In particular, the ATM option has a higher trading
frequency than the ITM and OTM options. It is also evident that the price
movements of options are strongly correlated with that of the underlying
stock.

We use the Asyn-MTS model to simultaneously forecast the joint distribution of
the ATM, OTM, and ITM option prices and the underlying stock price, each of
which is considered a unique variable. The dataset is obtained from HKEX,
comprising 3-minute level trading records for the stock and its options from
01/Dec/2014 to 31/Dec/2017. We divide the transaction data for each day into
one instance, resulting in a total of 720 instances, with each instance having
a time length of 110 (where each time step represents 3 minutes, and the
trading time for one day is 5.5 hours). We construct 16 variables, where one
variable represents the underlying stock price and the other 15 variables
represent options with moneyness ranging from 0.96 to 1.10.

Figure 16: Predicted ranges of quantile (shaded areas)of the Tencent stock
price and its options prices. Solid dots represent observed price points.

Table 4 presents the experiment results, showing that RFNs achieve lower
$\mathrm{CRPS}$ and $\mathrm{CRPS}_{\text{sum}}$ compared to the baselines.
The predicted quantile intervals for one-day transaction records are
visualized in Figure 16. Although options with moneyness 0.96 and 1.10 have no
trading, we can also predict their trends by leveraging the cross-sectional
information from other options. These results suggest that the RFN
specification applies well to option price prediction by capturing the
dependence structure between irregularly transacted prices of the underlying
stock and its options.

Table 4: Transaction Records of the Tencent Stock and its Options (HKEX).
(Same reporting formats as in Table 1)

model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ | model | $\mathrm{CRPS}$ | $\mathrm{CRPS}_{\text{sum }}$ |  $\mathrm{CRPS}$$\uparrow$ |  $\mathrm{CRPS}_{\text{sum }}$$\uparrow$
---|---|---|---|---|---|---|---
GRUODE | 0.1960 $\pm$ 0.0103 | 0.5645 $\pm$ 0.0435 | RFN-GRUODE | 0.1645 $\pm$ 0.0035 | 0.4221 $\pm$ 0.0305 | 16.05% | 25.22%
ODELSTM | 0.1947 $\pm$ 0.0116 | 0.5494 $\pm$ 0.0463 | RFN-ODELSTM | 0.1654 $\pm$ 0.0074 | 0.4138 $\pm$ 0.0239 | 15.06% | 24.69%
GRU-D | 0.2044 $\pm$ 0.0081 | 0.5677 $\pm$ 0.0268 | RFN-GRU-D | 0.1675 $\pm$ 0.0101 | 0.4199 $\pm$ 0.0344 | 18.07% | 26.02%
ODERNN | 0.1850 $\pm$ 0.0161 | 0.5208 $\pm$ 0.0720 | RFN-ODERNN | 0.1633 $\pm$ 0.0041 | 0.4070 $\pm$ 0.0179 | 11.73% | 21.85%

##  8 Conclusion

Forecasting the joint distribution of multivariate time series with temporal
irregularities and heterogeneous dependence structure is challenging. These
irregularities can arise due to insufficient observation or inherent
randomness in event occurrences. The perspective we take is that we allow the
arrival times of observations, which are the core of the temporal
irregularities, to play a central role in the model building process.

The proposed RFN framework consists of a marginal learning block and a
multivariate learning block. The marginal learning block captures individual
serial dependencies by assigning unique hidden states to component variables
and updating them based on arrival times of events. This approach acknowledges
the intrinsic nature of temporal irregularities in the data and avoids using
imputation methods. The multivariate learning block uses a conditional
formulation of the Continuous Normalizing Flow model to incorporate the time-
varying hidden states learned from the marginal learning block, as well as the
concurrent information, into learning the non-parametric distribution. The
log-likelihood objective of RFNs fully acknowledges the irregular spacing and
asynchronous pattern in the data.

In the comprehensive evaluations of RFNs using real-world datasets from
application domains such as finance, physics, and climate, RFNs consistently
outperform models that assume Gaussian distributions for modeling dependence
structures in irregularly sampled time series across various state-of-the-art
specifications, e.g., GRU-D, GRU-ODE-Bayes, ODELSTM, and ODERNN. These
findings highlight the robustness and effectiveness of RFNs.

##  Appendix A Marginal learning blocks

Throughout (A.1),(A.2),(A.3), the operator $\odot$ is the element-wise
multiplication and the operator $\circledast$ is defined as in (8).

###  A.1 GRU-D

GRU-D [che2018recurrent] modifies the classical GRU by adding the trainable
exponential decays,

$\gamma_{t}=\exp\\{-\max(\mathsf{0},\mathbf{w}_{\gamma}\delta_{t}+\mathsf{b}_{\gamma})\\},$

where $\delta_{t}$ is the time intervals between two observations, and
$\mathbf{w}_{\gamma}$ and $\mathsf{b}_{\gamma}$ are trainable parameters.
First, it imputes the missing values with the weighted average between the
last observations and the empirical means,

$\hat{x}_{t}^{d}=m_{t}^{d}x_{t}^{d}+\left(1-m_{t}^{d}\right)\left(\gamma_{x_{t}}^{d}x_{t^{\prime}}^{d}+\left(1-\gamma_{x_{t}}^{d}\right)\tilde{x}^{d}\right),$

where $m_{t}^{d}$ is the mask of $d^{\text{th}}$ variable at time $t$,
$\gamma_{x_{t}}^{d}$ is the trainable decay, $x_{t^{\prime}}^{d}$ is the
previous observation ($t^{\prime}<t$), and $\tilde{x}^{d}$ is its empirical
mean. Second, it employs hidden state decay to further capture the missing
patterns in the hidden state,

$\hat{\mathbf{h}}_{t-1}=\gamma_{h_{t}}\odot\mathbf{h}_{t-1}.$

Applying the marginal learning idea to GRU-D, the updating equations are as
follows

$\displaystyle\mathbf{r}_{t}$
$\displaystyle=\sigma\left(\mathcal{W}_{r}\circledast\mathsf{\hat{x}}_{t}+\mathcal{U}_{r}\circledast\mathbf{\hat{h}}_{t-1}+\mathcal{V}_{r}\circledast\mathbf{m}_{t}+\mathbf{b}_{r}\right),$
$\displaystyle\mathbf{z}_{t}$
$\displaystyle=\sigma\left(\mathcal{W}_{z}\circledast\mathsf{\hat{x}}_{t}+\mathcal{U}_{z}\circledast\mathbf{\hat{h}}_{t-1}+\mathcal{V}_{z}\circledast\mathbf{m}_{t}+\mathbf{b}_{z}\right),$
$\displaystyle\mathbf{\tilde{h}}_{t}$
$\displaystyle=\tanh\left(\mathcal{W}_{h}\circledast\mathsf{\hat{x}}_{t}+\mathcal{U}_{h}\circledast(\mathbf{r}_{t}\odot\mathbf{\hat{h}}_{t-1})+\mathcal{V}_{h}\circledast\mathbf{m}_{t}+\mathbf{b}_{h}\right),$
$\displaystyle\mathbf{h}_{t}$
$\displaystyle=\mathbf{z}_{t}\odot\mathbf{h}_{t-1}+(1-\mathbf{z}_{t})\odot\mathbf{\tilde{h}}_{t},$

where
$\mathcal{W}_{\\{r,z,h\\}},\mathcal{U}_{\\{r,z,h\\}},\mathcal{V}_{\\{r,z,h\\}},\mathbf{b}_{\\{r,z,h\\}}$
are parameters to be trained,

$\displaystyle\mathcal{W}_{\\{r,z,h\\}}$
$\displaystyle=\\{\mathbf{w}_{\\{r,z,h\\}}^{1};\cdots;\mathbf{w}_{\\{r,z,h\\}}^{D}\\}\in\mathbb{R}^{D\times
S\times 1},$ $\displaystyle\mathcal{U}_{\\{r,z,h\\}}$
$\displaystyle=\\{\mathbf{u}_{\\{r,z,h\\}}^{1};\cdots;\mathbf{u}_{\\{r,z,h\\}}^{D}\\}\in\mathbb{R}^{D\times
S\times S},$ $\displaystyle\mathcal{V}_{\\{r,z,h\\}}$
$\displaystyle=\\{\mathbf{v}_{\\{r,z,h\\}}^{1};\cdots;\mathbf{v}_{\\{r,z,h\\}}^{D}\\}\in\mathbb{R}^{D\times
S\times 1}.$

###  A.2 ODERNN

ODERNN [rubanova2019latent] proposes handling irregularly sampled time series
by modeling the hidden state evolution using a Neural ODE in the absence of
observations. Compared to GRU-ODE-Bayes [de2019gru], such a Neural ODE does
not admit explicit solutions; instead, it is modeled using a standard MLP
network. When there are observations, a classical GRU cell will update the
hidden state.

When implementing ODERNN to learn the marginal block, we assign each variable
with unique training parameters, and thus the evolution formula of ODERNN in
the marginal learning block is

$\frac{d\mathsf{h}^{d}(t)}{dt}=\mathsf{f}^{d}_{\theta}(\mathsf{h}^{d}(t),t).$

Let $\mathbf{h}_{t_{-}}=\mathbf{h}(t_{-})$ and
$\mathbf{h}_{t_{+}}=\mathbf{h}(t_{+})$, the updating equations of ODERNN for
the marginal learning are

$\displaystyle\mathbf{r}_{t_{-}}$
$\displaystyle=\sigma\left(\mathcal{W}_{r}\circledast\mathsf{x}_{t}+\mathcal{U}_{r}\circledast\mathbf{h}_{t_{-}}+\mathbf{b}_{r}\right),$
$\displaystyle\mathbf{z}_{t_{-}}$
$\displaystyle=\sigma\left(\mathcal{W}_{z}\circledast\mathsf{x}_{t}+\mathcal{U}_{z}\circledast\mathbf{h}_{t_{-}}+\mathbf{b}_{z}\right),$
$\displaystyle\mathbf{\tilde{h}}_{t_{-}}$
$\displaystyle=\tanh\left(\mathcal{W}_{h}\circledast\mathsf{x}_{t}+\mathcal{U}_{h}\circledast(\mathbf{r}_{t}\odot\mathbf{h}_{t_{-}})+\mathbf{b}_{h}\right),$
$\displaystyle\mathbf{h}_{t_{+}}$
$\displaystyle=\mathbf{z}_{t}\odot\mathbf{h}_{t_{-}}+(1-\mathbf{z}_{t})\odot\mathbf{\tilde{h}}_{t-},$

where
$\mathcal{W}_{\\{r,z,h\\}},\mathcal{U}_{\\{r,z,h\\}},\mathbf{b}_{\\{r,z,h\\}}$
are the trainable parameters,

$\displaystyle\mathcal{W}_{\\{r,z,h\\}}$
$\displaystyle=\\{\mathbf{w}_{\\{r,z,h\\}}^{1};\cdots;\mathbf{w}_{\\{r,z,h\\}}^{D}\\}\in\mathbb{R}^{D\times
S\times 1},$ $\displaystyle\mathcal{U}_{\\{r,z,h\\}}$
$\displaystyle=\\{\mathbf{u}_{\\{r,z,h\\}}^{1};\cdots;\mathbf{u}_{\\{r,z,h\\}}^{D}\\}\in\mathbb{R}^{D\times
S\times S}.$

###  A.3 ODELSTM

ODELSTM [lechner2020learning] is designed to learn the long-term dependencies
in irregularly sampled time series. In the absence of observations, the hidden
states are evolved by a Neural ODE, and when there is an observation, the
hidden state is updated by the classical LSTM. Unlike the ODERNN and GRU-ODE-
Bayes, the cell state of LSTM will contain long-term memory. To cast the
ODELSTM in the RFN framework, we assign each variable a unique set of training
parameters. It has the following evolution equation

$\frac{d\mathsf{h}^{d}(t)}{dt}=\mathsf{f}^{d}_{\theta}(\mathsf{h}^{d}(t),t).$

Let $\mathbf{h}_{t_{-}}=\mathbf{h}(t_{-})$ and
$\mathbf{h}_{t_{+}}=\mathbf{h}(t_{+})$. The updating equations of ODELSTM for
the marginal learning are

$\displaystyle\mathbf{f}_{t_{-}}$
$\displaystyle=\sigma\left(\mathcal{W}_{f}\circledast\mathsf{x}_{t}+\mathcal{U}_{f}\circledast\mathbf{h}_{t_{-}}+\mathbf{b}_{f}\right),$
$\displaystyle\mathbf{i}_{t_{-}}$
$\displaystyle=\sigma\left(\mathcal{W}_{i}\circledast\mathsf{x}_{t}+\mathcal{U}_{i}\circledast\mathbf{h}_{t_{-}}+\mathbf{b}_{i}\right),$
$\displaystyle\mathbf{o}_{t_{-}}$
$\displaystyle=\sigma\left(\mathcal{W}_{o}\circledast\mathsf{x}_{t}+\mathcal{U}_{o}\circledast\mathbf{h}_{t_{-}}+\mathbf{b}_{o}\right),$
$\displaystyle\mathbf{\tilde{c}}_{t_{-}}$
$\displaystyle=\tanh\left(\mathcal{W}_{c}\circledast\mathsf{x}_{t}+\mathcal{U}_{c}\circledast\mathbf{h}_{t_{-}}+\mathbf{b}_{c}\right),$
$\displaystyle\mathbf{c}_{t_{+}}$
$\displaystyle=\mathbf{f}_{t_{-}}\odot\mathbf{c}_{(t-1)_{+}}+\mathbf{i}_{t_{-}}\odot\mathbf{\tilde{c}}_{t_{-}}$
$\displaystyle\mathbf{h}_{t_{+}}$
$\displaystyle=\mathbf{o}_{t_{-}}\odot\tanh(\mathbf{c}_{t_{+}}),$

where $\mathcal{W}_{(f,i,o,c)},\mathcal{U}_{(f,i,o,c)},\mathbf{b}_{(f,i,o,c)}$
are the parameters to be trained,

$\displaystyle\mathcal{W}_{(f,i,o,c)}$
$\displaystyle=\\{\mathbf{w}_{(f,i,o,c)}^{1};\cdots;\mathbf{w}_{(f,i,o,c)}^{D}\\}\in\mathbb{R}^{D\times
S\times 1},$ $\displaystyle\mathcal{U}_{(f,i,o,c)}$
$\displaystyle=\\{\mathbf{u}_{(f,i,o,c)}^{1};\cdots;\mathbf{u}_{(f,i,o,c)}^{D}\\}\in\mathbb{R}^{D\times
S\times S}.$

##  Appendix B Proof of Lemma 5.1

###### Proof.

Consider the random variables $\mathsf{X}\in\mathbb{R}^{D}$ and
$\mathsf{Y}\in\mathbb{R}^{C}$. Let
$\tilde{\mathsf{Z}}(s)=[\mathsf{Z}(s),\mathsf{Y}(s)]^{T}$ be a finite
continuous random variable. Let the realizations of $\tilde{\mathsf{Z}}(s)$ be
$\tilde{\mathsf{z}}(s)=[\mathsf{z}(s),\mathsf{y}(s)]^{T}$, and the probability
density of $\tilde{\mathsf{Z}}(s)$ be
$p(\tilde{\mathsf{z}}(s))=p(\mathsf{z}(s),\mathsf{y}(s))$ which depends on the
flow time $s$, where $s_{0}\leq s\leq s_{1}$. We assume $\mathsf{z}(s)$
evolves continuously in the real space, starting from a sample point of a pre-
defined distribution at $s=s_{0}$ and ending at a sample point $\mathsf{x}$ of
$\mathsf{X}$ at $s=s_{1}$. According to (LABEL:eqt:ffjord1), the governing
equation of $\tilde{\mathsf{z}}(s)$ can be written as

$\displaystyle\frac{\partial\tilde{\mathsf{z}}(s)}{\partial
s}=\left[\begin{array}[]{c}\frac{\partial\mathsf{z}(s)}{\partial s}\\\
\frac{\partial\mathsf{y}(s)}{\partial
s}\end{array}\right]=\left[\begin{array}[]{c}\mathsf{f}(\mathsf{z}(s),s,\mathsf{y}(s);\theta)\\\
0\end{array}\right],\text{ where }s_{0}\leq s\leq s_{1},$

$\mathsf{y}(s)=\mathsf{y}$ and $\mathsf{z}(s_{1})=\mathsf{x}$. From the
differential equation $\frac{\partial\mathsf{y}(s)}{\partial s}=0$, we have
$\mathsf{y}(s)=\mathsf{y}$. Hence, we are going to find the density
$p\left(\tilde{\mathsf{z}}_{s}\right)=p\left(\mathsf{z}_{s},\mathsf{y}\right)$
such that

$\frac{\partial\mathsf{z}(s)}{\partial
s}=\mathsf{f}(\mathsf{z}(s),s,\mathsf{y};\theta).$

The proof requires the change of variables in the probability density theorem.
We restate the result below.

######  Result.

Suppose that $G(\cdot)$ is a bijective function and differentiable. Given the
variables $(\mathsf{Z},\mathsf{Y})$ and the corresponding density function
$p(\mathsf{z},\mathsf{y})$, the density function of
$(G(\mathsf{Z}),\mathsf{Y})$ is

$\displaystyle
p(G(\mathsf{z}),\mathsf{y})=p(\mathsf{z},\mathsf{y})\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}}G(\mathsf{z})&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}^{-1}$ (25)
$\displaystyle\Rightarrow$ $\displaystyle\log p(G(\mathsf{z}),\mathsf{y})=\log
p(\mathsf{z},\mathsf{y})-\log\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}}G(\mathsf{z})&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}.$

Now, we consider $\frac{\partial\log p(\tilde{\mathsf{z}}(s))}{\partial
s}=\frac{\partial\log p(\mathsf{z}(s),\mathsf{y})}{\partial s}$. We write
$\tilde{\mathsf{z}}(s+\epsilon)=[\mathsf{z}(s+\epsilon),\mathsf{Y}]^{{}^{\prime}}=[T_{\epsilon}(\mathsf{z}(s)),\mathsf{Y}]^{{}^{\prime}}$.
From the first principle of derivatives, we have

$\displaystyle\partial_{s}\log p(\tilde{\mathsf{z}}(s))$
$\displaystyle=\underset{\epsilon\rightarrow 0^{+}}{\lim}\frac{\log
p(\tilde{\mathsf{z}}(s+\epsilon))-\log p(\tilde{\mathsf{z}}(s))}{\epsilon}$
$\displaystyle=\underset{\epsilon\rightarrow 0^{+}}{\lim}\frac{\log
p(T_{\epsilon}(\mathsf{z}(s)),\mathsf{y})-\log
p(\mathsf{z}(s),\mathsf{y})}{\epsilon}.$

We simplify the quantity $\underset{\epsilon\rightarrow 0^{+}}{\lim}\frac{\log
p(T_{\epsilon}(\mathsf{z}(s)),\mathsf{y})-\log
p(\mathsf{z}(s),\mathsf{y})}{\epsilon}$. Applying equation
(LABEL:eqt:log_density_difference), we have

$\displaystyle\partial_{s}\log p(\tilde{\mathsf{z}}(s))$
$\displaystyle=\underset{\epsilon\rightarrow
0^{+}}{\lim}\frac{-\log\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{\mathsf{z}}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}}{\epsilon}$
$\displaystyle\overset{\text{L'H\^{o}pital}}{=}-\underset{\epsilon\rightarrow
0^{+}}{\lim}\partial_{\epsilon}\log\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times
D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}=-\underset{\epsilon\rightarrow
0^{+}}{\lim}\frac{\partial_{\epsilon}\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times
D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}}{\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}}$
$\displaystyle=\underbrace{\frac{\underbrace{-\underset{\epsilon\rightarrow
0^{+}}{\lim}\partial_{\epsilon}\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times
D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}}_{\text{bounded}}}{\underset{\epsilon\rightarrow
0^{+}}{\lim}\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times
D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}}}_{1}=-\underset{\epsilon\rightarrow
0^{+}}{\lim}\partial_{\epsilon}\biggl{|}\det\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times D}&\mathbb{I}_{C}\end{bmatrix}\biggl{|}.$

Applying the Jacobi’s formula, we have

$\displaystyle\partial_{s}\log p(\tilde{\mathsf{z}}(s))$
$\displaystyle=-\underset{\epsilon\rightarrow
0^{+}}{\lim}\operatorname{Tr}\biggl{(}\operatorname{adj}\bigg{(}\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times
D}&\mathbb{I}_{C}\end{bmatrix}\bigg{)}\times\partial_{\epsilon}\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times D}&\mathbb{I}_{C}\end{bmatrix}\biggl{)}$
$\displaystyle=-\operatorname{Tr}\biggl{(}\underbrace{\underset{\epsilon\rightarrow
0^{+}}{\lim}\operatorname{adj}\bigg{(}\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times
D}&\mathbb{I}_{C}\end{bmatrix}\bigg{)}}_{\mathbb{I}_{C}}\times\underset{\epsilon\rightarrow
0^{+}}{\lim}\frac{\partial}{\partial\epsilon}\begin{bmatrix}&\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))&\mathsf{0}_{D\times
C}\\\ &\mathsf{0}_{C\times D}&\mathbb{I}_{C}\end{bmatrix}\biggl{)}$
$\displaystyle=-\operatorname{Tr}\biggl{(}\underset{\epsilon\rightarrow
0^{+}}{\lim}\bigg{(}\partial_{\epsilon}\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))\bigg{)}\biggl{)}.$

Applying Taylor series expansion on $T_{\epsilon}(\mathsf{z}(s))$ and taking
the limit, we have

$\displaystyle\partial_{s}\log p(\tilde{\mathsf{z}}(s))$
$\displaystyle=-\operatorname{Tr}\biggl{(}\underset{\epsilon\rightarrow
0^{+}}{\lim}\bigg{(}\partial_{\epsilon}\partial_{\mathsf{z}(s)}T_{\epsilon}(\mathsf{z}(s))\bigg{)}\biggl{)}$
$\displaystyle=-\operatorname{Tr}\biggl{(}\underset{\epsilon\rightarrow
0^{+}}{\lim}\bigg{(}\partial_{\epsilon}\partial_{\mathsf{z}(s)}(\mathsf{z}(s)+\partial_{s}\mathsf{z}(s)\epsilon+O(\epsilon^{2}))\bigg{)}\biggl{)}$
$\displaystyle=-\operatorname{Tr}\biggl{(}\underset{\epsilon\rightarrow
0^{+}}{\lim}\biggl{(}\partial_{\epsilon}(\mathbb{I}+\partial_{\mathsf{z}(s)}f(\mathsf{z}(s),s,\mathsf{y};\theta)\epsilon+O(\epsilon^{2}))\bigg{)}\biggl{)}$
$\displaystyle=-\operatorname{Tr}\biggl{(}\partial_{\mathsf{z}(s)}f(\mathsf{z}(s),s,\mathsf{y};\theta)\biggl{)}.$

As such, we have

$\displaystyle\int_{s_{0}}^{s_{1}}\partial_{s}\log
p(\tilde{\mathsf{z}}(s))ds=\int_{s_{0}}^{s_{1}}-\operatorname{Tr}\biggl{(}\partial_{\mathsf{z}(s)}\mathsf{f}(\mathsf{z}(s),s,\mathsf{y};\theta)\biggl{)}ds$
$\displaystyle\Rightarrow$ $\displaystyle\log
p(\tilde{\mathsf{z}}(s_{1}))=\log
p(\tilde{\mathsf{z}}(s_{0}))+\int_{s_{1}}^{s_{0}}\operatorname{Tr}\biggl{(}\partial_{\mathsf{z}(s)}\mathsf{f}(\mathsf{z}(s),s,\mathsf{y};\theta)\biggl{)}ds$
$\displaystyle\Rightarrow$ $\displaystyle\log
p(\mathsf{z}(s_{1}),\mathsf{y})=\log
p(\mathsf{z}(s_{0}),\mathsf{y})+\int_{s_{1}}^{s_{0}}\operatorname{Tr}\biggl{(}\partial_{\mathsf{z}(s)}\mathsf{f}(\mathsf{z}(s),s,\mathsf{y};\theta)\biggl{)}ds.$

∎

##  Appendix C More Experiment Results

In section 7.3 (c) Data and Model Validation, the prediction intervals of
Figure 9 are restricted to quantiles ranging from 0.2 to 0.8. Here, we present
results with a broader quantile range in Figure 17.

(a) Syn-MTS with $0.1\leq q\leq 0.9$ (b) Asyn-MTS with $0.1\leq q\leq 0.9$ (c)
Syn-MTS with $0.2\leq q\leq 0.8$ (d) Asyn-MTS with $0.2\leq q\leq 0.8$

(e) Syn-MTS with $0.3\leq q\leq 0.7$ (f) Asyn-MTS with $0.3\leq q\leq 0.7$ (g)
Syn-MTS with $0.4\leq q\leq 0.6$ (h) Asyn-MTS with $0.4\leq q\leq 0.6$

Figure 17: Predicted ranges of quantile (shaed areas) of the simulated Syn-MTS
and Asyn-MTS at various probability intervals $q$.

##  Appendix D Algorithms of Training and Sampling of Syn-MTS and Asyn-MTS

The training process should proceed sequentially to compute the hidden state
and the joint distribution of all variables at each observed time point.

For example, in a Syn-MTS instance with two variables observed at time points
$t_{1}$, $t_{3}$, $t_{5}$, we first compute the hidden state
$\mathbf{h}_{t_{1-}}$ and feed the observation $\mathsf{x}_{t_{1}}$ into the
flow model $\mathsf{f}_{\theta}$ to obtain the joint distribution at $t_{1}$.
Then, we evolve the hidden states to $t_{3}$ and feed the observation
$\mathsf{x}_{t_{3}}$ into the flow model $\mathsf{f}_{\theta}$ to obtain the
joint distribution at $t_{3}$, and so on so force. However, this process is
time-consuming since we must call the flow model $t_{K}$ times.

To reduce the complexity, we first evolve the hidden states to the final time
point and store all the necessary parameters at each time point. For example,
we compute and store the base distribution parameters at each time point,
e.g.,
$\\{\mu_{t_{1}},\Sigma_{t_{1}}\\}=\textbf{MLP}(\mathbf{h}_{t_{1-}};\psi)$,
$\\{\mu_{t_{3}},\Sigma_{t_{3}}\\}=\textbf{MLP}(\mathbf{h}_{t_{3-}};\psi)$, and
$\\{\mu_{t_{5}},\Sigma_{t_{5}}\\}=\textbf{MLP}(\mathbf{h}_{t_{5-}};\psi)$. We
then concatenate all the observations and the stored parameters together,
i.e.,
$\mathsf{x}_{t}=\\{\mathsf{x}_{t_{1}},\mathsf{x}_{t_{3}},\mathsf{x}_{t_{5}}\\},\mu_{t}=\\{\mu_{t_{1}},\mu_{t_{3}},\mu_{t_{5}}\\},\Sigma_{t}=\\{\Sigma_{t_{1}},\Sigma_{t_{3}},\Sigma_{t_{5}}\\}$.
Finally, we feed $\mathsf{x}_{t}$ to the flow model $\mathsf{f}_{\theta}$ to
transform the values at all time points simultaneously. This approach requires
only one call to the flow model.

Algorithm 1 Algorithm for training and sampling process of Syn-MTS (example of
RFN-GRUODE)

1:Training (Syn-MTS):

2:Input: Observations: $\\{\mathbf{x}_{i}\\}_{i=1}^{N}$; Masks:
$\\{\mathbf{m}_{i}\\}_{i=1}^{N}$; Observed time points:
$\\{\mathsf{t}_{i}=[t_{1},\cdots,t_{K_{i}}]\\}_{i=1}^{N}$; Flow time interval
$[s_{1},s_{0}]$

3:Initialize: $\text{time}=0$, $\mathbf{h}_{0}$, and all trainable parameters
$\zeta^{c},\zeta^{u},\psi,\theta$

4:for $i=1$ to $N$ do

5:for $k=1$ to $K_{i}$ do

6:$\mathbf{h}(t_{k-})=\text{Continuous
Updating}(\mathbf{h}(t_{(k-1)+}),\text{time},t_{k};\zeta^{c})$ % marginal
hidden state evolves to $t_{k}$

7:$\mu_{t_{k}},\Sigma_{t_{k}}=\textbf{MLP}(\mathbf{h}_{t_{k-}};\psi)$ %
predict and store the parameters of base distribution at $t_{k}$

8:$\text{time}=t_{k}$

9:$\mathbf{h}(t_{k+})=\text{Discrete
Updating}(\mathbf{h}(t_{k-}),\mathsf{x}_{t_{k}},\mathsf{m}_{t_{k}};\zeta^{u})$
% marginal memory updates at $t_{k}$

10:end for

11:% concatenate observations and parameters at all time points and tranform
them at the same time by the flow model

12:Set
$t\in\\{t_{1}^{i},\cdots,t_{K_{i}}^{i}:\mathsf{m}_{t_{k}}=1,k\in\\{1,\cdots,K_{i}\\}\\}$
% all the observed time points

13:$\mathsf{z}_{t}=\mathsf{x}_{t}+\int_{s_{1}}^{s_{0}}\mathsf{f}(\mathsf{z}(s),s,\mathsf{h}_{t-};\theta)ds$
% transform observations from data distribution to base distribution (from
$s_{1}$ to $s_{0}$)

14:$\log p(\mathsf{z}_{t}|\mathbf{h}_{t-})=\log
p(\mathsf{z}_{t};\mu_{t},\Sigma_{t})$ % compute the log-likelihood of
transformed observations in base distribution

15:$\log p(\mathsf{x}_{t}|\mathbf{h}_{t-})=\log
p(\mathsf{z}_{t}|\mathbf{h}_{t-})+\int_{s_{1}}^{s_{0}}\operatorname{Tr}[\partial_{\mathsf{z}(s)}\mathsf{f}]ds$
% compute the log-likelihood of observed data points

16:$\mathcal{L}^{i}_{\text{Syn-MTS}}=-\log
p(\mathsf{x}_{t}|\mathbf{h}_{t_{-}})$ % compute the loss of sample $i$ by
equation (17)

17:end for

18:$\mathcal{L}_{\text{Syn-
MTS}}=\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}^{i}_{\text{Syn-MTS}}$ % compute the
total loss by averaging the loss of all the samples

19:$\zeta^{c},\zeta^{u},\psi,\theta\leftarrow\arg\min_{\zeta^{c},\zeta^{u},\psi,\theta}\mathcal{L}_{\text{Syn-
MTS}}$ % optimize the training parameters via stochastic gradient descent
algorithm

20:

21:Sampling (Syn-MTS):

22:Input: Observations $\mathbf{x}_{i}$; Masks $\mathbf{m}_{i}$; Observed time
points $\mathsf{t}_{i}=[t_{1},\cdots,t_{K_{i}}]$; Flow time interval
$[s_{0},s_{1}]$; Trained model $\mathsf{f}$

23:Initialize: $\text{time}=0$, $\mathbf{h}_{0}$

24:for $k=1$ to $K_{i}$ do

25:$\mathbf{h}(t_{k-})=\text{Continuous
Updating}(\mathbf{h}(t_{(k-1)+}),\text{time},t_{k})$ % marginal hidden state
evolves to $t_{k}$

26:$\mu_{t_{k}},\Sigma_{t_{k}}=\textbf{MLP}(\mathbf{h}_{t_{k-}};\psi)$ %
predict the parameters of base distribution at $t_{k}$

27:$\mathsf{z}_{t_{k}}\sim N(\mu_{t_{k}},\Sigma_{t_{k}})$ % sampling samples
in predicted base distribution

28:$\mathsf{x}_{t_{k}}=\mathsf{z}_{t_{k}}+\int_{s_{0}}^{s_{1}}\mathsf{f}(\mathsf{z}(s),s,\mathsf{h}_{t_{k-}};\theta)ds$
% transform the samples from base distribution to data distribution (from
$s_{0}$ to $s_{1}$)

29:$\text{time}=t_{k}$

30:$\mathbf{h}(t_{k+})=\text{Discrete
Updating}(\mathbf{h}(t_{k-}),\mathsf{x}_{t_{k}},\mathsf{m}_{t_{k}})$ %
marginal memory updates at $t_{k}$

31:end for

32:return $\mathsf{x}_{t}$

Algorithm 2 Algorithm for training and sampling process of Asyn-MTS (example
of RFN-GRUODE)

1:Training (Asyn-MTS):

2:Input: Observations: $\\{\mathbf{x}_{i}\\}_{i=1}^{N}$; Masks:
$\\{\mathbf{m}_{i}\\}_{i=1}^{N}$; Observed time points:
$\\{\mathsf{t}_{i}=[t_{1},\cdots,t_{K_{i}}]\\}_{i=1}^{N}$; Flow time interval
$[s_{1},s_{0}]$

3:Initialize: $\text{time}=0$, $\mathbf{h}_{0}$, and all trainable parameters
$\zeta^{c},\zeta^{u},\psi_{m},\theta,\eta$

4:for $i=1$ to $N$ do

5:for $k=1$ to $K_{i}$ do

6:$\mathbf{h}(t_{k-})=\text{Continuous
Updating}(\mathbf{h}(t_{(k-1)+}),\text{time},t_{k};\zeta^{c})$ % marginal
hidden state evolves to $t_{k}$

7:$\mu_{t_{k}}^{d},\Sigma_{t_{k}}^{d}=\textbf{m-MLP}(\mathbf{h}_{t_{k-}},\mathsf{x}_{t_{k}}^{-d};\psi_{m})$
% predict the parameters of base distribution at $t_{k}$, $\mu$ and $\Sigma$
are masked-MLPs

8:$\hat{\mathsf{x}}_{t_{k}}=g(\mathbf{h}_{t_{k-}};\eta)$ % predict the values
of variables at $t_{k}$

9:$\tilde{\mathsf{x}}_{t_{k}}=\mathsf{x}_{t_{k}}\odot\mathsf{m}_{t_{k}}+\hat{\mathsf{x}}_{t_{k}}\odot(1-\mathsf{m}_{t_{k}})$
% impute the variables that are unobserved at $t_{k}$

10:$\text{time}=t_{k}$

11:$\mathbf{h}(t_{k+})=\text{Discrete
Updating}(\mathbf{h}(t_{k-}),\mathsf{x}_{t_{k}},\mathsf{m}_{t_{k}};\zeta^{u})$
% marginal memory updates at $t_{k}$

12:end for

13:% concatenate observations and parameters at all time points and tranform
them at the same time by the flow model

14:for $d=1$ to $D$ do

15:Set
$t\in\\{t_{1}^{i},\cdots,t_{K_{i}}^{i}:m^{d}_{t_{k}}=1,k\in\\{1,\cdots,K_{i}\\}\\}$
% all the observed time points for variable $d$

16:$z^{d}_{t}=x^{d}_{t}+\int_{s_{1}}^{s_{0}}\mathsf{f}(z^{d}(s),s,\mathsf{h}_{t_{k-}},\tilde{\mathsf{x}}^{-d}_{t};\theta)ds$
% transform data distribution of observed variable $d$ to base distribution

17:$\log p(z_{t}^{d}|\mathbf{h}_{t_{-}},\tilde{\mathsf{x}}_{t}^{-d})=\log
p(z_{t}^{d};\mu_{t}^{d},\Sigma_{t}^{d})$ % compute the log-likelihood of
observations of variable $d$ in base distribution

18:$\log p(x_{t}^{d}|\mathbf{h}_{t_{-}},\tilde{\mathsf{x}}_{t}^{-d})=\log
p(z_{t}^{d}|\mathbf{h}_{t_{-}},\tilde{\mathsf{x}}_{t}^{-d})+\int_{s_{1}}^{s_{0}}\operatorname{Tr}[\partial_{z^{d}(s)}\mathsf{f}]ds$
% compute the log-likelihood of observed data points

19:$\mathcal{L}_{\text{Asyn-MTS}}^{d,i}=-\log
p(x_{t}^{d}|\mathbf{h}_{t_{-}},\tilde{\mathsf{x}}_{t}^{-d})$ % only compute
the negative log-likelihood loss for the observed variables

20:$\mathcal{L}_{\text{MSE}}^{d,i}=(x_{t}^{d}-\hat{x}_{t}^{d})^{2}$ % compute
the MSE loss for imputation

21:end for

22:$\mathcal{L}^{i}_{\text{Asyn-MTS}}=\sum_{d=1}^{D}(\mathcal{L}_{\text{Asyn-
MTS}}^{d,i}+\mathcal{L}_{\text{MSE}}^{d,i})$ % compute the loss of sample $i$
by equation (22)

23:end for

24:$\mathcal{L}_{\text{Asyn-
MTS}}=\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}^{i}_{\text{Asyn-MTS}}$ % compute
the total loss by averaging the loss of all the samples

25:$\zeta^{c},\zeta^{u},\psi_{m},\theta,\eta\leftarrow\arg\min_{\zeta^{c},\zeta^{u},\psi_{m},\theta,\eta}\mathcal{L}_{\text{Asyn-
MTS}}$ % optimize the training parameters via stochastic gradient descent
algorithm

26:

27:Sampling (Asyn-MTS):

28:Input: Observations: $\mathbf{x}_{i}$; Masks: $\mathbf{m}_{i}$; Observed
time points: $\mathsf{t}_{i}=[t_{1},\cdots,t_{K_{i}}]$; Flow time interval
$[s_{0},s_{1}]$; Trained model $\mathsf{f}$

29:Initialize: $\text{time}=0$, $\mathbf{h}_{0}$

30:for $k=1$ to $K_{i}$ do

31:$\mathbf{h}(t_{k-})=\text{Continuous
Updating}(\mathbf{h}(t_{(k-1)+}),\text{time},t_{k};\zeta^{c})$ % marginal
hidden state evolves to $t_{k}$

32:for $d=1$ to $D$ do

33:
$\mu_{t_{k}}^{d},\Sigma_{t_{k}}^{d}=\textbf{m-MLP}(\mathbf{h}_{t_{k-}},\mathsf{x}^{-d}_{t_{k}})$
% predict the parameters of base distribution for variable $d$ at $t_{k}$

34: $z_{t_{k}}^{d}\sim N(\mu_{t_{k}}^{d},\Sigma_{t_{k}}^{d})$ % sampling
samples in predicted base distribution

35:
$x_{t_{k}}^{d}=z_{t_{k}}^{d}+\int_{s_{0}}^{s_{1}}\mathsf{f}(z^{d}(s),s,\mathsf{h}_{t_{k-}},\mathsf{x}_{t_{k}}^{-d};\theta)ds$
% transform the samples from base distributions to data distributions

36:end for

37:$\text{time}=t_{k}$

38:$\mathbf{h}(t_{k+})=\text{Discrete
Updating}(\mathbf{h}(t_{k-}),\mathsf{x}_{t_{k}},\mathsf{m}_{t_{k}};\zeta^{u})$
% marginal memory updates at $t_{k}$

39:end for

40:return $\mathsf{x}_{t}$