Year: 2021

Abstract:
We propose a novel probabilistic framework for modeling stochastic dynamics with the rigorous use of stochastic optimal control theory. The proposed model called the neural Markov controlled stochastic differential equation (CSDE) overcomes the fundamental and structural limitations of conventional dynamical models by introducing the following two components: (1) Markov dynamic programming to efficiently train the proposed CSDE and (2) multi-conditional forward-backward losses to provide information for accurate inference and to assure theoretical optimality. We demonstrate that our dynamical model efficiently generates a complex time series in the data space without extra networks while showing comparable performance against existing model-based methods on several datasets.

heading: 1 INTRODUCTION
text: Recently, there has been interest in using continuous dynamical systems to approximate complex time series. Neural ODE Chen et al. (2018), which opened the way for continuous representation of neural networks, have been widely investigated and thoroughly analyzed by Massaroli et al. (2020). As the stochastic generalization of ODE, Neural SDEs Li et al. (2020) have been proposed by regarding intrinsic stochasticity in data representations (e.g., stock market data). Since the conventional Neural ODE/SDEs only utilize the initial information of trajectories when propagating dynamics, modelling complex time-series with naive Neural ODE/SDEs has been regarded as inefficient and undesirable choices, as pointed out by Kidger et al. (2020).
To address these problems, Rubanova et al. (2019) presented an auto-regressive model to generalize recurrent neural networks (RNNs) to have continuous hidden dynamics with neural ODE. Furthermore, Chen et al. (2018) proposed an encoder-decoder structure with Neural ODE in the latent space to reconstruct/predict complex data representation. Although the aforementioned approaches produce remarkable results, they focus on suggesting additional probabilistic structures rather than improving the learnability of the Neural ODE model itself. Compared to aforementioned approaches, we focus on solving the fundamental issues of Neural ODE/SDEs. First, we raise two important questions.
Q1) How can we construct an efficient network architecture for Neural ODE/SDE models that do not require additional recurrent networks to model complex time series?
Q2) How can we train Neural ODE/SDEs that can utilize richer information of observed sequences to accurately generate complex time series?
As SDEs can be posed as stochastic generalizations of ODEs, we focus on a stochastic framework and adopt the stochastic optimal control theory as our primary analysis tool for the rigorous and systematic analysis of the aforementioned problems. Keeping this in mind, the contributions of our paper are to answer the above two questions. A1) Novel probabilistic framework for stochastic dynamics. We propose a novel neural controlled stochastic differential equation (CSDE) to model the complex stochastic time series, where multiple control agents are defined to construct local dynamics in their own private temporal states. With this property, the proposed CSDE incorporates Markov dynamic programming, enables our model to directly infer the complex trajectory on data space rather than the latent space without any extra network (e.g., encoders/decoders), and shows remarkable efficiency compared to existing methods.
A2) Novel conditional losses. We introduce a novel Markov forward conditional (MFcond) loss to utilize multi-conditioned dynamics instead of the conventional dynamics determined by partial initial conditions. The proposed MFcond loss makes our method to model the complex information of time-sequential data. To impose regularization and to ensure the optimality of control agents, we also suggest a novel Markov backward conditional (MBcond) loss.

heading: 2 RELATED WORK
text: ODE As a Latent Probabilistic Model. Rubanova et al. (2019) suggested an ODE-RNN by combining RNN with the latent dynamics induced by the Neural ODE. To deal with irregular time-stamps, exponential-decaying of the hidden states was also discussed by Che et al. (2018). De Brouwer et al. (2019) assumed that the observations are sampled from the stochastic dynamics induced from SDEs and introduced GRU-ODE to approximate the observed stochastic time series.
SDE As a Latent Probabilistic Model. Liu et al. (2021) incorporated Neural SDEs with recurrent models as a primary probabilistic dynamical model to generate stochastic continuous-time latent variables. While this SDE model could describe the stochastic dynamics on the latent space with recurrent structures (e.g., RNN encoder/decoder), it required a whole sequence of historical observations as inputs to the model. Unfortunately, this type of formulation leads to non-Markov types of SDEs, which makes it difficult to analyze the probabilistic characteristics of the dynamics. Unlike this model, we focus on the Markov SDEs while maintaining identical objectives.
Neural CDE and RDE. Kidger et al. (2020) proposed a data-driven neural controlled differential equation called Neural CDE to incorporate a rough-path analysis theory and model complex time series. Morrill et al. (2021) extended the rough-path theory with a Neural RDE to deal with the continuous time series over long time.
Generative SDE Models. Recently, Kidger et al. (2021) suggested SDE-based generative adversarial networks (GANs). Park et al. (2021) utilized the temporal conditional Wasserstein distance to construct GANs for time-series generation.
Please refer to Appendix A.1 for additional discussion on related works.

heading: 3 MARKOV NEURAL CONTROLLED SDE
text: In Section 3.1, we introduce a novel SDE model that considers temporally private agents. In Section 3.2, we propose the Markov-DP-TP framework to efficiently solve the stochastic optimal control problem with the proposed neural SDE model. Finally, we suggest novel Markov conditional forward and backward losses in Section 3.3 and 3.4, respectively. In the Appendix, we provided the detailed technical definitions.

heading: 3.1 CONTROLLED STOCHASTIC DIFFERENTIAL EQUATIONS
text: The basic object of our interest is a controlled Ft-adapted process Xαt with multiple control agents α = {α1, · · · , αM} ∈ A where A denotes the set of admissible control agents. In particular, the stochastic process Xαt is defined as a solution to the following CSDE:
dXαt = M∑ i=1 wi(t)b i ( t,Xαt , α i ) dt+ M∑ i=1 wi(t)σ i ( t,Xαt , α i ) dWt, (1)
where b and σ : [0, T ]×Rd×A→ Rd are the drift and diffusion functions, respectively. Each control agent αi : [0, T ]× Rd × Rm, αi = αi(t,Xt; θi),∀1 ≤ i ≤ M is defined as a Markov closed-loop feedback control, which is parameterized by the neural network θi. While every agent is defined as a closed-loop feedback-type Carmona (2016b), the solution to the CSDE above, Xαt , is the Markov process, which means that process Xαt is propagated using the information of the current state.
Let T = {tk}1≤k≤N be a set of ordered times1 such that 0 = t1 < · · · < tk < tl < · · · < tN = T . The set of functions {wi(t)}1≤i≤M is defined as an indicator function on the intervals, wi(t) = 1tk≤t≤tl with predetermined starting/ending points tk, tl in T. We call this function temporal privacy (TP) because it represents each agent’s attention on different sub-intervals. Overall, in (1), the stochastic process Xαt is propagated by summing M -number of individual agent’s weighted attentions {∑M wib i(·, ·, αi), ∑M wiσ i(·, ·, αi) } . To understand the behavior of the proposed
CSDE more deeply, we consider the following detailed example:
1The time interval dt ≈ ∆t = |tk − tl| for any k, l can be set regularly/irregularly in our method.
Role of Temporal Privacy. We define wr(s) = 1t≤s≤u, t, u ∈ T with r ≤ M . Then, Xαu in (1) given Xt at an interval [t, u] can be equivalently rewritten in the integration form:
Xα=[α 1,··· ,αM ]
u = X α t + ∫ u t M∑ i wi(s)b i(s,Xαs , α i)ds+ ∫ u t M∑ i wi(s)σ i(s,Xαs , α i)dWs
= Xαt + ∫ wr(s) br(s,Xαs , α r)ds+ ∫ wr(s) σr(s,Xαt , α r)dWs = X αr u .
(2)
In (2), the activated control agent to evaluate the stochastic process Xαu for the interval [t, u] is only αr (i.e., Xαu = X αr
u ) owing to the definition of the weighting function w(·)(t). This means that the remaining control agents {αj}j 6=r are not used for the evaluation of the stochastic process in the sub-interval [t, u]. While each agent αi is activated at its own private sub-interval, this leads our method to adopt dynamic programming (DP) to train Neural CSDEs in the form of (1). In this paper, we aim to solve the optimal control problem via DP with multiple agents, where each agent specializes in solving a particular sub-problem in its private interval.

heading: 3.2 MARKOV DYNAMIC PROGRAMMING PRINCIPLES
text: The dynamic programming principle is one of the fundamental philosophies for dealing with stochastic optimal control problems. Its basic idea is to consider a family of sub-problems with different initial times/states and establish the relation among the sub-problems to systemically solve them. Using the mathematical property of the proposed CSDE with TP, we present an efficient learning strategy to solve stochastic optimal control problems via Markov dynamic programming (Markov-DP).
In this paper, we aim to solve the stochastic optimal control problem by training control agents α = [α1, · · · , αM ] and minimizing the cost functional J(t,Xαt ) : [0, T ]× Rd → R+:
J(t,Xαt ) = E [∫ T t l(s,Xαs )ds+ Ψ(X α T ) ∣∣∣∣∣Ft ] = E [∫ u t l(s,Xαs )ds+ J(u,X α u ) ∣∣∣∣Xαt ] , (3) where l : [0, T ]×Rd → R+ is the running cost (e.g., L2 loss) that computes the discrepancy between the propagated process Xαt and the observed data point yt at each time t, Ψ(X α T ) : Rd → R+ is the terminal cost that estimates the discrepancy between the terminal state and the data yT . To evaluate the cost functional J(t,Xαt ) at time t with control agents α, the running cost is integrated over the time interval [t, T ] conditioned on filtration Ft. Note that the expectation conditioned on Ft in (3) can be substituted to the expectation conditioned on Xαt in light of the Markov property presented in Section A.2, and the cost functional at time t only depends on the current state of the process Xαt .
Markov-DP with Temporal Privacy. By combining the tower property of the conditional expectations with the dynamic programming principle and Itô’s formula (Oksendal (1992)), one can show that a minimization problem can be recursively decomposed into sub-problems owing to the property of TP in our proposed CSDE:
V (t,Xαt ) , inf α J(t,Xαt ) = inf α E [∫ u
t
l(s,Xαs )ds+ J(u,X α u ) ∣∣∣∣Xαt ]︸ ︷︷ ︸ (A)
= inf αr
E [∫ u
t
l(s,Xαs )ds ∣∣∣∣Xαt ]︸ ︷︷ ︸ (B) + inf α(−r) E[J(u,Xαu )|Xαt ]︸ ︷︷ ︸ (B’) ,
(4)
where V is an optimal cost functional (i.e., value function), αr denotes the r-th control agent, and α(−r) = [α1, · · · , 0, · · · , αM ] indicates the set of remaining agents (the r-th component is zero). In (4), the minimization problem (A) over α is divided into two sub-problems using the dynamic programming principle, which are (B) and (B’). Because the minimization problem (B) is only dependent on the control agent αr parameterized by the neural network θr, we compute the gradient descent of θr to solve the sub-problem (B):
θrk+1 = θ r k −
∂
∂θr E [∫ {s:wr(s)=1} l ( s,X αr(·,·,θrk) s ) ds ∣∣∣∣∣Xαt ] , (5)
where wr(s) = 1t≤s≤u is the TP function at an interval [t, u] and k is the index for the learning iterations. In (5), the r-th control agent αr minimizes the cost functional using the gradient descent scheme at its own temporal sub-interval. As the remaining sub-problem (B’) over agents α(−r) can also be recursively decomposed into smaller sub-problems using the dynamic programming principle, the original problem (A) is solved separately with M -number of control agents α = {α1, · · · , αM} with the M -number of gradient descent schemes. This indicates that we can obtain the set of agents α? = {αi(·, ·; θi?)} by collecting individual optimal agents with sub-problems. In this paper, we combine the Markov-DP with M gradient descent schemes in (5) and CSDE with TP in (1) and introduce a novel Markov-DP-TP framework. In the numerical experiments in Section 4.4, we show that the proposed Markov-DP-TP framework remarkably increases the model efficiency compared to conventional non-DP naive approaches, which makes our method directly model the complex time series in the data space. However, despite the improvements with our novel Markov-DP-TP framework, there exist remaining practical/theoretical issues that should be addressed to solve the optimal control problem with complex datasets.
1) Conditional Dependency. The main practical issue in implementing the Markov-DP-TP framework is that explicit conditional states are not given, e.g., Xαt in (5). As different initial/terminal conditions of SDE lead to totally different behaviors of induced dynamics, well-designed conditional information is a crucial factor in training the Neural CSDE for specific applications. In Section 3.3, we introduce the Markov Forward conditional (MFcond) loss to train the Neural CSDE with well-posed conditional information that ensures accurate network predictions.
2) Theoretical Optimality. In the optimal control theory, there are well-known partial differential equations called Hamiltonian-Jacobi-Bellman (HJB) equations, which assure the theoretical optimality of control agents. If the control agents can solve the HJB equation, the proposed method attains the optimal state Vt(Xαt ) = infα Jt(X α t ) = Jt(X α? t ). However, the optimal agents α? of the proposed CSDE with gradient descent are not generally equivalent to the solution to the HJB equation. In Section 3.4, we propose the Markov Backward conditional (MBcond) loss to assure the optimality of control agents and to provide information in backward dynamics for regularization.

heading: 3.3 MARKOV FORWARD CONDITION
text: In this section, we first raise the important question: Why is the well-posed conditional estimation in cost functional important to accurately train Neural SDE (CSDE) models? To elucidate the importance of this question, we consider the following minimization problem with the cost functional with naive partial information:
inf α L(α) = inf α Ey0 [∫ T 0 l(s,Xαs )ds+ Ψ(X α T ) ∣∣∣∣∣X0 = y0 ] , (6)
where y(·) = {yt}t∈[0,T ] denotes a set of observed data, and y0 is the initial data at time t = 0. In (6), the conditional expectation is taken to the single initial state X0 = y0, and the control agents minimize the accumulated losses using this partial information. As pointed out by Kidger et al. (2020), this naive cost functional causes a problem when dealing with high-dimensional complex datasets. This is because the Neural CSDE should disentangle the inherent latent information of complex high-dimensional data to generate accurate results, but the control agents are trained with only the restrictive and partial information of the observed data (i.e., initial condition X0 = y0). To solve this problem, we introduce a novel loss function called the MFcond loss that can fully exploit the information of the given observed data y(·), while keeping the Markov structure of Xαt : Definition 1. (MFcond loss) We define the prediction operator T αs,t as follows, for s < t,
T αs,t := 1 |I(s, t)| ∑
m∈I(s,t)
[ Xαtm + ∫ t tm M∑ i=1 wib i(u,Xαu , α)du+ ∫ t tm M∑ i=1 wiσ i(u,Xαu , α)dW (m) u ∣∣∣∣Xαtm = ytm ] ,
(7) where I(s, t) := {m : s ≤ tm < t}, |I(s, t)| is the cardinality of I(s, t), and {W (m)u }m∈I(s,t) denotes the Wiener processes with respect to time u. Let us define a random stopping time τs such that τs := inft{t : l(t, T αs,t) > } for the pre-determined threshold 2. Then, we can define the MFcond loss with the stopping time τ(·), as follows:
2Please refer to Appendix (A.7) for detailed information
Lf ( α,y(·) ) = Ey(·) [∫ T t l ( τs, T αs,τs ) χ(s)ds+ Ψ(XαT ) ] , (8)
where χ(s) is an indicator function that produces values at the observed time (i.e., χ(s) = 1 if ys is observed at s; otherwise, χ(s) = 0). This function is used to consider the irregularly sampled data points. In (8), naive running cost l of (6) is replaced with l ◦ T αs,τs , in which the MFcond loss recursively accumulates the expected future losses l ◦ T αs,τs conditioned on multiple observations. At each time s, stopping time τs decides the future time to stop the CSDE propagation by determining if the accumulated losses are larger than the predetermined threshold or not. While the proposed loss requires a set of multiple conditions on the Markov process Xαt to train control agents, information is utilized to generate time-sequential data, and complex dynamics can be expressed. A conceptual illustration of the proposed MFcond loss is shown in Figure 1-(a).
The main idea of our MFcond loss in (8) is to minimize the differences between the future estimations Xα,su for any given s ≤ u. In other words, the proposed CSDE is trained to generate an identical future estimation of Xαu given any past initial conditions X α (·) = y(·), i.e., (X α,s u ≈ Xα,tu ,∀s ≤ t ≤ u) to estimate network inference with multiple conditions in the test time. This idea is used to introduce a novel inference procedure to overcome the raised issues on the partial information.
Network Inference. Let {ytm} be the observed data sequences until the current time t in the test dataset. Our objective is to predict the future points {ŷtk}, (tm ≤ t < tk). Our model generates the stochastic estimation X̂tk to approximate ŷtk at a future time given multiple initial conditions ŷtm :
ŷtk ≈ X̂ α tk = T α tm,tk =
1 |I| ∑
s∈I(tm,tk)
[ Xαts + M∑ i=1 ∫ tk ts wib i(t,Xαt , α i)dt+ M∑ i=1 ∫ tk ts wiσ i(t,Xαt , α i)dW (s) t ∣∣∣∣Xαts = ŷts ] (9)
Network Inference
In (9), each control agent makes decisions on its specialized temporal state and collaborates to generate a stochastic conditional estimation X̂αtk and approximate ŷtk . As our MFcond loss induces identical estimations Xα,ŷtmtk for any tm, X̂ α tk
utilizes multiple conditions {ŷtm} and fully exploits the past information to predict/estimate future values. A conceptual illustration of the network
inference is shown in Figure 1-(b). While the proposed inference mechanism utilizes enlarged information3 compared to a single initial condition, it can model the complex time-sequential data.
If the control agents are trained with the naive cost functional, the terminal states Xα,su (conditioned on initial state Xs = ys) and Xα,tu (conditioned on initial state Xt = yt) are largely different, which causes problems when we generate complex time-sequential data during the test time, whereas our inference mechanism introduced in (9) utilizes averaged multi-decisions Xαtk given different initial conditions. Thus, the MFcond loss is essential for utilizing the proposed inference procedure.
Unlike the dynamical auto-regressive probabilistic models (e.g., ODE-RNNs) that encode whole (or partial) data sequences, as shown in (1), the proposed Markovian CSDE model only uses the current observation to propagate stochastic dynamics. An additional inference mechanism coordinates the multi-conditioned trajectories to utilize information and produces complex time series.

heading: 3.4 MARKOV BACKWARD CONDITION
text: In the previous section, we suggested the Markov forward conditional loss that exploits the entire information of time-series data to generate accurate results. Aside from its empirical benefits to some applications, no theoretical/empirical optimality of (4) is assured by minimizing the MFcond loss in general. To tackle this problem, in this section, we further introduce the additional stochastic dynamics relating optimality of proposed CSDE-TP.
Let us define the auxiliary process Zt = V (t,Xα?t ) with a value function V , where α? denotes the optimal control agents. Subsequently, we consider the following forward-backward stochastic differential equations (FBSDEs):
(Xα?t , Zt) =  dXα?t = ∑M i=1 wib i(t,Xα?t , α i ?)dt+ ∑M i=1 wiσ i(t,Xα?t , α i ?)dWt dZt = −l(s,Xα?t )dt+ ∑M i=1∇V (t,Xt)wiσi(t,X α? t , α i ?)dWt
ZT = Ψ(X α? T )
(10)
FBSDEs
The first SDE (i.e., Xα?t ) called the forward SDE has an identical form of (1) and propagates stochastic evaluation in the forward direction with optimal control agents. The second SDE (i.e., Zt) called backward SDE recursively subtracts the running cost from the terminal state Ψ(Xα?T ) in the backward direction using forward estimations Xα?t and cancels the effect of martingales in the diffusion term. We utilize the property of backward dynamics Zt to train the control agents for the following reasons.
1) Backward Multi-conditions. Like the MFcond loss with multi-conditions in the forward direction, we want to provide additional information to backward dynamics to train the control agents.
2) Approximated Solution of HJBE. The auxiliary process Zt gives the theoretical optimality for control agents related to the HJB equation based on the results developed in Yong & Zhou (1999); Pardoux & Tang (1999), where the process Zt = V (t, ·) admits a solution of the HJB equation in (11) and induces an optimal solution for the minimization problem infα J in (4).
∂V (t, x)
∂t +
1 2 Tr[σTσ(t, x, α?)∇2V (t, x)] +∇V (t, x)T b(t, x, α?) + l(t, x) = 0, (11)
where V (T, x) = Ψ(x). In (11), we want to approximate Zt using control agents for optimality. However, the process Zt requires optimal control agents α? that cannot be obtained during the training time. To overcome this problem, we approximate the auxiliary process Zt with Zαt parameterized by neural control agents α(·, ·, θ), which is defined as the modified version of Zt. In particular, Zt can be expressed in the following integral form:
Zαt = Ψ(X α T )− ∫ t T M∑ i wi(s)l(s,X α s )ds+ ∫ t T M∑ i wi(s)σ i(s,Xαs , α i)∇J(s,Xαs )dWTs , (12)
where J is the cost functional defined in (3), and∇J denotes the gradient of the cost functional with respect to its spatial axis. Using the proposed process Zαt , we introduce a novel loss function called the MBcond loss to satisfy the two objectives discussed above.
3Please refer to detailed explanation in Appendix A.3.
Algorithm 1 Neural Markov CSDE-TP Require: γ = 0.95,
for k = 1 to K (i.e., the total number of training iterations) do 1) Simulate forward controlled SDE with Markov control agents 1-1) dXαkt = ∑M i=1 wib i(t,X αk t , α i k)dt+ ∑M i=1 wiσ i(t,X αk t , α i k)dWt
1-2) Evaluate each decision of control agents αik = αik(t,X αk t ; θ i k) 1-3) Compute the MFcond loss for M -control agents {Lf (αik(·, ·, θik))} with stopping time τ(·) 1-4) Update threshold for random stopping time, k+1 ← 12 max l ( t, T αks,t (ys) ) 2) Simulate backward controlled SDE 2-1) dZαkt = − ∑M i wil(s,X αk t )dt+ ∑M i=1∇J(t,X αk t )wiσ
idWt 2-2) Evaluate the MBcond loss for M -control agents, {Lb(αik(·, ·, θik))}1≤i≤M
3) Update control agents with Markov-DP 3-1) θik+1 = θik − γ∇θiLf (αi(·, ·, θik))− (1− γ)∇θiLb(αi(·, ·, θik))
end for
Definition 2. (MBcond loss) Let us define the auxiliary process Zαt as the solution to (12). Then, the MBcond loss can be defined as follows:
Lb(α) = Ey(·),t∈[0,T ] [ |Zαt |2 ∣∣∣Xt = yt]. (13) Theoretically, if we optimize the MBcond loss (13) according to the proposed backward dynamics Zαt , the PDE reformulation of backward dynamics, called Non-linear Feynmann-Kac, have the identical solution4 to HJB equation in (11). Thus, our method can attain the optimal solution of original problem posed in section 3.2.
Intuitively saying, one can show that the MBcond loss is equivalent to the reformulation of the minimization problems in (4) using Itô’s formula. Thus, solving the minimization problem infα Lb induces an identical effect to solve the original problem infα J . The only difference is that we utilize multiple conditions to provide conditional information on the backward dynamics Zαt for the regularization of control agents trained with forward conditional dynamics and to impose constraints on control agents, which induces an approximated solution to the HJB equation.

heading: 3.5 OBJECTIVE FUNCTION
text: In this section, we describe the overall training procedure, which incorporates all the proposed components (i.e., Markov-DP with CSDE-TP, MFcond loss, and MBcond loss) as follows:
inf α L(α)︸ ︷︷ ︸
MFBcond
= inf α=[α1,··· ,αM ] γLf (α)︸ ︷︷ ︸ MFcond + (1− γ)Lb(α)︸ ︷︷ ︸ MBcond
CSDE-TP ≈ M∑ i inf αi γLf ([αi, α(−i)]) + [ (1− γ)Lb([αi, α(−i)]) ] ,
(14)
where Lf and Lb are defined in (8) and (13), respectively, and γ is a balancing hyperparameter. In (14), the control agents α = [α1, · · · , αM ] are trained with a convex combination of MFcond and MBcond losses. Using the property of CSDE-TP with Markov-DP, the original problem is approximated with the collection of M sub-problems, and each control agent is separately trained with M gradient descent schemes. Algorithm 1 describes the detailed procedure of our method.

heading: 4 EXPERIMENTS
text: Network structure of control agents. The neural network structure for each agent control consists of 2-layers of fully-connected layers, where each module has 128 latent dimensions. For the activation units, we used the specialized module LipSwish, Chen et al. (2019); Kidger et al. (2021), to stabilize the FBSDEs during training. Please refer to Appendix A.6 for detailed information on the network architecture. Datasets. For the evaluations, we used PhysioNet, Speech Commands, Beijing Air-Quality, and S&P500 Stock Market datasets. Refer to Appendix A.5 for data statistics and prepossessing procedures.
4Please refer to Appendix A.4 for the discussion on theoretical optimality induced by the MBcond loss.

heading: 4.1 TIME-SEQUENTIAL DATA RECONSTRUCTION
text: In this experiment, we compared our model against baseline dynamic models: [Latent ODE, Chen et al. (2018)], [Latent SDE, Li et al. (2020)], [ODE-RNN, Rubanova et al. (2019)], [GRU-D, Che et al. (2018)], [mTAND, Shukla & Marlin (2021)], and [ODE2VAE, Çağatay Yıldız et al. (2019)]. We used open-source codes provided by the authors for comparison. For the Latent ODE (SDE) methods, RNN and ODE-RNN were used for the encoder structures, where the decoder structures were identically set to ODE (SDE). Table 1 shows the performance of all baseline methods compared to the proposed CSDE-TP for the reconstruction tasks. As evaluation metrics, we used mean squared errors (MSE) and negative log-likelihood (NLL) with open-source code in Rubanova et al. (2019). As shown in Table 1, the proposed method consistently outperformed the baseline methods by a large margin. In this experiment, we observed that latent dynamics-based methods (e.g., Latent ODE/SDE with RNN and ODE-RNN encoders) on models attained similar performances. We set the latent dimensions of each control agent to 128 for both the reconstruction and prediction experiments. In the experiments on both datasets, the Mckean-Vlasov (MV) type of the SDE model slightly improved the performance, where it subtracted the mean (i.e., mean-shifting) of the control agent outputs to normalize/reduce the intrinsic volatility in the inferred process X̂αtk .

heading: 4.2 TIME-SEQUENTIAL DATA PREDICTION
text: In this experiment, all the methods were trained with the first 36 time-sequences (e.g., 80%) on the training dataset, and evaluated to predict the remaining 12 time sequences (e.g., 20%) on the test dataset. We followed a procedure similar to that suggested by Rubanova et al. (2019) to predict the future time series. Table 2 compares the performance in terms of MSE and NLL, in which the proposed method considerably outperforms the baseline methods. In the experiment on the Air quality dataset, we set a smaller latent dimension of 64 as the feature dimension compared to those of the other datasets. Interestingly, the CSDE-TP-MV model exhibited worse performance than the vanilla CSDE-TP model when the data dimension was very small. This result indicates that the mean-shifting in a low-dimensional data space reduces the expressiveness of the CSDEs.

heading: 4.3 UNCERTAINTY ESTIMATION ON STOCK MARKET DATASET
text: When high volatility is observed over the temporal/spatial axes, conventional evaluation metrics such as MSEs hardly capture the stochastic property of the time-series variations. Thus, to capture the stochasticity, we evaluated the distance between the distributions of the test data and the inferred/generated data using the maximum mean discrepancy (MMD). We followed the protocol
suggested by Li et al. (2017) to evaluate the MMD distance, where we used two Gaussian RBF kernels with bandwidths of [5.0, 10.0]. Using this evaluation metric, we experimented on reconstruction tasks using the S&P-500 Stock Market dataset. Table 3 shows that the proposed CSDE-TP outperforms baselines and effectively recovers the distributional information of stock prices with the stochastic property of the SDE models and the proposed optimization framework. Interestingly, the latent SDE model attains better performance compared to the Latent ODE, as it utilizes an additional Wiener process to model the data uncertainty. The performance improvement of the Latent SDE vanishes when we remove the diffusion term (σ = 0) of the latent SDE.

heading: 4.4 EMPIRICAL STUDY
text: Efficiency of the Markov-DP-TP framework. To show the empirical advantages of our CSDE-TP model with Markov-DP learning schemes, we evaluated our CSDE-TP according to a different number of control agents on the prediction task using the Air Quality dataset. Figure 2-(a) shows the training MSEs for several variants of the proposed model in the first 20 epochs, where CSDE-TPShallow1, -Shallow2, and -Deep (i.e., black, blue, and red lines) denote the proposed models with a different number of control agents, i.e., M = 2, 8, and 48, respectively. The standard CSDE model (i.e., the black dashed line) utilized a single agent,M = 1. For all models, the total number of training parameters was equivalently set to ≈ 40K, and the number of parameters was normalized. As shown in Figure 2-(a), despite using the same number of parameters, employing multiple agents clearly outperforms the standard CSDE in terms of the learning curve. From this fact, we can conclude that the Markov-DP-TP significantly increased the network efficiency compared to the standard CSDE, which indicates that our Markov-DP framework is crucial for training controlled dynamics models. Efficiency of the MFcond loss. In this experiment, we show the empirical advantages of the multiconditioned CSDE in (8) against the naive partial-conditioned CSDE in (6). Similar to previous experiments, the results were obtained for the prediction task with the Air Quality dataset. Figure 2-(b) shows the model confidence in testing MSEs for the first 50 epochs, where shaded areas indicate the confidence regions (i.e., ± std). The proposed MFcond loss exhibits considerable performance improvement (.08 .87) compared to the conventional native cost functional and reduces the variances in loss landscapes with stable learning. With the theoretical discussion in Appendix A.3, we conclude that the proposed CSDE actively exploits the information of the complex time series with multiple conditions to accurately generate complex time-series.

heading: 5 CONCLUSION
text: In this paper, we introduce a novel Markov-type CSDE with the TP function that records the individual attention of each control agent at sub-intervals along the temporal axis. Using the properties of the CSDE and TP, we suggest Markov DP to efficiently train the control agents by decomposing the original problem into smaller sub-problems. To overcome the practical/theoretical issues, we propose two novel losses, namely, MFcond and MBcond losses. The MFcond loss captures the future time to estimate the running costs, while multiple conditions are actively provided to forward dynamics. The MBcond loss assures the theoretical optimality of the control agents and imposes regularization by providing additional information to backward dynamics. Experimental results demonstrate the efficiency of the proposed method for various tasks using real datasets.

heading: A APPENDIX

heading: A.1 DETAILED COMPARISON TO EXISTING METHODS
text: In this section, we investigate the relation between our method and existing methods.
Reverse SDE vs. Backward SDE. Song et al. (2020) suggested a novel SDE called reverse SDE, which shares semantically similar idea with BSDE: both reverse/backward SDEs enhance the forward SDE by providing additional information to drift/diffusion functions in forward dynamics.
The mathematical motivation of the reverse SDE in Anderson (1982) is to pose the SDEs with Wiener processes Wt, Ŵt with respect to these minimal increasing/decreasing sigma algebras At, Ât and define the relation between them:
dŴt = 1
pt(Xt) ∇[pt(Xt)σ(t,Xt)]dt+ dWt, (15)
where Xt is a solution to the forward SDE and pt is the probability density of Xt. Using the relation in (15), the reverse SDE transforms the prior distribution (e.g., Gaussian noise distribution) back into the data distribution (e.g., 2D images) by gradually removing the noises and reconstruct the original data with the well-designed score function (i.e., ∇pt(x)) in backward dynamics. In contrast to the reverse SDE, the role of backward SDE in this paper is to consider the probabilistic reformulation to access the cost functional to provide the additional information in backward dynamics.
Stacked ODE vs. CSDE-TP. Massaroli et al. (2020) suggested the stacked Neural ODE that shares similar idea with the proposed CSDE-TP, where temporally piece-wise neural nets are considered to model the complex dynamics. However, the stacked ODE faces the aforementioned problem on partial conditional information when generating complex data as their models only take initial values to propagate dynamics. As opposed to their models, the proposed model is trained with multiple observations and directly generates time series in data space without any latent embedding network. Furthermore, we generalize their optimal control problems to the stochastic version and propose the Markov-DP-TP framework that can systemically solve the problem.
DDPMs vs. CSDE-TP. Tashiro et al. (2021) suggested denoising denoising diffusion models (DDPMs) that are conditioned on the set of the observed data, where the generated sequential data is assumed to be gradually transformed from an initial state in the forward direction and the backward process is parameterized by the neural network and trained to minimize specific ELBOs.
Specifically, the transition probability pθ(Xt−1|Xt) in the backward process is defined as a parameterized Gaussian distribution:
pθ(Xt|Xt−1) = N (Xt−1;µθ(t,Xt), σθ(t,Xt)), (16) where mean and covariance (µθ, σθ) are parameterized by the neural network θ. Similar to the proposed CSDE, their parameterized functions are closed-loop type processes and the whole probabilistic sequential model pθ is posed as the Markov chain: pθ(X0:T ) = pθ(XT ) ∏T t=1 pθ(Xt−1|Xt).
In contrast to the DDPM, the probability transition in the proposed CSDE is defined as the continuous generalization called controlled Fokker-Planck equation (CFPE):
∂ ∂t pαθ (x, t|y, s) = −∇f(x, t, αθ)pt(x) + 1 2 Tr [ ∇2σσT (x, t, αθ) · pt(x) ] , (17)
where t > s ∈ [0, T ) and x, y ∈ Rd, pt ∼ Xαt is the probability distribution ofXαt . The CFPE in (17) one-to-one corresponds to the CSDE (i.e., Xαt ). Compared to the discrete-time Gaussian transition model, this conditional probability can express complex continuous-time probability transitions while maintaining the Markov structure.

heading: A.2 NOTATIONS AND BACKGROUND
text: We first define the basic definitions of probabilistic objects: Definition 3. A filtration {Ft} is an increasing sequence of {Ft} of σ-algebra such that F0 · · · ,⊂ Ft ⊂ F . The triplet (Ω,Ft,P) is called a filtered probability space. Definition 4. The filtration generated by Wiener process Wt is defined as FWt = σ{W0, · · · ,Wt}. In this case, Wt is naturally FWt -adapted by construction.
Definition 5. The stochastic process {Xt} is called {Ft}-adapted if Xt is Ft measurable for every 0 ≤ t ≤ T .
Throughout this paper, we work on the filtered probability space (Ω, {Ft}t∈[0,T ],P) with the ddimensional Ft-Wiener process Wt and natural filtration FWt . We assume that αi for all 1 ≤ i ≤M is admissible Markov control, (i.e., αi is Ft-adapted and αi ∈ Ai, Xαt has a unique solution). Definition 6. (Markov Process) Let Xt be a Ft-adapted stochastic process. Then, Xt is the Markov process if the following equality holds:
E[Xt|Fs] = E[Xt|Xs], ∀s ≤ t. (18) Definition 7. (Controlled Stochastic Differential Equation)
Xαt = Xs + ∫ t s b (u,Xαu , α) du+ ∫ t s σ (u,Xu, α) dWu, for 0 ≤ s ≤ t ≤ T. (19)
The solution to the above CSDE is denoted as Xα,st . If the initial states are specified (i.e., starting point Xs = x), we denote the solution as X α,s,x t . By the definition of Markovian control agents, in all cases, the solution to the proposed CSDE in (1) is a Markov process.
Mathematical Assumptions. In this paper, we assume that functions b, σ are uniformly Lipschitz continuous along its spatial axis and bounded ‖b(t, 0; ·)‖ , ‖σ(t, 0; ·)‖ at the entire interval [0, T ]. We assume that each functions bi(·, x, ·), σi(·, x, ·),Ψ(x), l(·, x) are twice differentiable for all 1 ≤ i ≤ M , (i.e., , bi, σi,Ψ, l ∈ C2(Rn). and both drift and diffusion functions are uniformly Lipschitz on its spatial axis, i.e., bi, ∂xbi, ∂2xb i, σi, ∂xσ i, ∂2xσ
i ∈ Lip), and the trainable parameters of the control agents αi, θi are lying in the compact subset C of its ambient space (i.e., θi ∈ C ⊂ Rm).

heading: A.3 ENLARGED INFORMATION BY COLLECTION OF OBSERVED DATA
text: In the proposed inference procedure, we define a novel operator T in (9) to consider the multiconditioned dynamics with the Markov-type SDE model. Although this operator plays a central role in the paper, its mathematical properties have not been carefully dealt and investigated thoroughly. In this section, we discuss the relation between this operator and the enlarged information that is obtained by collecting past observations. In addition, we generalize the inference mechanism in (9) to a mathematically rigorous form and discuss the effect of the proposed operator T by showing some probability inequality.
Suppose that we have two observed conditional states {Xtm}, {Xtn} until the current time t, (tn, tm < t < tk) and the objective is to predict/generate the future value ytk using this information. We consider the deterministic time tk by replacing random stopping time τtm to simplify the discussion. First, we define the two-parameter stochastic process Y to model the proposed operator T in an alternative way:
Ttm,tk = Y (tm, tn)(w) , 1
2
( Xα,tm1,tk (w) +X α,tn 2,tk (w) ) , (20)
where w ∈ Ω takes a value in the probability space. The stochastic process Y is the (Ftm ∨ Ftn)valued random variable for any fixed tm, tn < t by the definition, whereFtm∨Ftn , Σ(FM1∪FM2) is the composited smallest sigma algebra with two filtrations. In the definition, we assume that processes Xα,tm1,tk , X α,tn 2,tk
are derived from two independent Wiener processes Wt and Ŵt. Then, we can define the two-parameter martingale Zakai (1981); Khoshnevisan (2003) in the following form:
M(tm, tn)(w) = E [l(tk, Y (tm, tn))|Ftm ∨ Ftn ] . (21) By the definition of M, it can be easily shown that M is the reformulation of the MFcond loss for some fixed number of past observations. Note that M is truly a martingale because conditional estimations are summed in the definition of T . The control agents are trained to minimize the M given the information induced by past observations (i.e., composited filtration { ∨ tm<t
FMtm }), which indicates that the proposed inference procedure can infer the future value X̂αtk according to the enlarged information { ∨ tm<t
FMtm }. By the fact that M is a martingale with respect to the composited filtration, we obtain the following result using Doob’s maximal inequality:
1− 1 2η ( E [∥∥Xα,t1,tk − ytk∥∥]+ E [∥∥Xα,t2,tk − ytk∥∥]) ≤ P [sup tn<t sup tm<t E[l ◦ T |Ftm ∨ Ftn ] ≤ η ] ,
(22)
where the inequality shows that errors between the future value ytk and the generated samples X α,t tk at time tk are bounded by the maximal perturbation probability. As the control agents are trained to minimize the MFcond loss (i.e.,M) in the right-hand side of inequality, it renders the probabilistic bound of L2 errors at future time tk.

heading: A.4 DETAILED DISCUSSIONS ON THE MBCOND LOSS
text: In this section, we investigate the detailed theoretical structure of the MBcond loss and its fundamental rationale for the optimality of control agents. For this, we rephrase the cost functional in the general form:
J(t, x) = E [∫ T t l(s,Xαs ) + Ψ(X α T ) ∣∣∣∣∣Xt = x ] . (23)
The classical non-linear Feynman-Kac theorem in Yong & Zhou (1999) states that given the cost functional J with the control agents α, one can obtain the second-order parabolic partial differential equation from (23):
∂J ∂t + 〈∇J, b(t, x, α)〉+ 1 2 Tr [ σσT (t, x, α)∇2J ] + l(s, x) = 0, (24)
where 〈·, ·〉 denotes the inner product and the boundary condition is given as J(T, x) = Ψ(x). Subsequently, by applying Itô’s formula to (23), we obtain the following probabilistic formulation:
Ψ(XαT ) = J(s,Xs) + ∫ T t [ ∂J ∂t (s,Xs) + 1 2 Tr[σσT (s,Xs, α)∇2J ] + 〈b(s,Xs, α),∇J〉 ] dt
+ ∫ T t 〈σT (s,Xs, α)∇J(s,Xs), dWt〉 = − ∫ T t l(s,Xs)ds+ ∫ T s 〈σT (s,Xs, α)∇J(s,Xs), dWs〉.
(25)
By rearranging each term above, the backward stochastic differential equation is induced directly.
Zαt = J(t,X α t ) = Ψ(X α T ) + ∫ T t l(s,Xs)ds− ∫ T t 〈σT (t,Xt, α)∇J(s,Xs), dWs〉. (26)
Note that, in the main paper, we use the inverse sign convention, ∫ T t (·) = − ∫ t T
(·), to emphasize the backward direction. By using the formulations, the MBcond loss in (13) can be rewritten in full description as follows:
Lb(α,Zαt ) = ∫ [0,T ] Ey(·) [∣∣∣∣Ψ(XαT ) + ∫ T t l(s,Xαs )ds− ∫ T t 〈σT (t,Xαt , α)∇J(s,Xαs ), dWs〉 ∣∣∣∣2 ∣∣∣∣∣Xαt = yt ] dt, (27) where y(·) = (y1, · · · , yT ) ∼ p(y1, · · · , yT ) denotes the set of observed data. The regularization effect comes from the expectation evaluation of the third term in (26). Specifically, one can obtain the following equality by using the Itô’s isometry:
E ∣∣∣∣∣ ∫ T t 〈σT (t,Xt, α)∇J(s,Xs), dWs〉 ∣∣∣∣∣ 2 |Xt  = E[∫ T t ∥∥σT (t,Xt, α)∇J(s,Xs)∥∥2 dt|Xt] . (28) Because the MBcond loss is posed to minimize this additional martingale term in (28) in backward dynamics according to the forward dynamics Xαt , it reduces the over-confidence of generated timeseries. By the relation (t, x)→ J(t, x)→ Zαt → Lb(t, x) for any (t, x) ∈ [0, T ]× R+, the update rule for the MBcond loss can be expressed as follows:
θrk+1 = θ r k −
∂
∂θr
[ Lb ( s,X αr(·,·,θrk) s ) ds ] , (29)
where this formulation is similar to (5) and shows that gradient descent with respect to θr for the MBcond loss can be explicitly defined.
Admissible Control Set A. In previous discussions, we show the relation between J with BSDE dynamics Zt and the well-defined gradient descent. The next step is to define the proper control set A to relate the gradient descent with optimality.
Let us define the Hilbert space L2 , {ϕ(t, x; θ̄);Rn-valued Ft-progressively measurable ∀θ̄ ∈ C} with the norm ‖ϕ‖2L2 = E [∫ T 0 |ϕ(t, x; θ)|2dt ] < ∞. We assume that control agents αr is Lαr
Lipschitz on the parameter variable, i.e., ∥∥∥αr(·, ·; θrk,1)− αr(·, ·; θrk,2)∥∥∥L2 ≤ Lαr ∥∥∥θrk,1 − θrk,2∥∥∥ for any θrk,1 6= θrk,2 ∈ Rm and any 1 ≤ k ≤ K. In all case, we assume that any θrk lies in the compact subset C of Rm. Each functions bi(·, x, ·), σi(·, x, ·),Ψ(x), l(·, x) are twice differentiable for all 1 ≤ i ≤ M , (i.e., , bi, σi,Ψ, l ∈ C2(Rn). and both drift and diffusion functions are uniformly Lipschitz on its spatial axis, i.e., bi, ∂xbi, ∂2xb i, σi, ∂xσ i, ∂2xσ
i ∈ Lip). As we defined Ψ and l as usual Euclidean distance, regularity/uniform Lipschitzness for these functions are trivial.
For the fixed parameter θ, we define the r-th control agent as θr → αr(·, ·, θr) , αr(θ) ∈ L2. Truly, the image space of αr(θ) is the closed subspace of the Hilbert space L2 due to the Lipschitzness with compactness of θ.
Let θr(k) : N → Rm be the trajectory for the training parameters of the r-th control agent at learning iteration k. Without loss of generality, Jr[αr] = J(t,X (αr,α(−r)) t ). We define the Euclidean closed metric balls {Bδkr }k∈N centered at θ(k) with the radius δ k r < ∞ such that Bδkr = {ϑ ∈ Rm; ‖ϑ− θr(k)‖ ≤ δkr , θr(k) is local minimum of Jr[θ(k)]}. Let us consider the sub-sequence {θ(k̄)}k̄∈N̄ ⊆ {θr(k)}k∈N , which induces the strictly-decreasing cost functional {Jr[θ(k̄)]}k̄∈N̄ with the ordered index set N̄ . Then, the admissible control set A is defined as follows:
α , [α1, · · · , αr · · · , αM ] ∈ A = M⊗ r=1  K⋂ K̄ K̄⋃ j̄=1 αr(·, ·, Bδj̄ ); j̄ ≤ K̄ ∈ N̄  ⊂ M⊗ r=1 L2, (30)
where K̄ is the maximal element in N̄ and the constant K ∈ N indicates the last iteration index of training defined in Algorithm 1. Intuitively, the control set A can be understood as a collection of local minimum obtained by M gradient descent schemes during training.
V , J [α?] = J [α(θ(K))] = inf α∈A J [α(θ)]. (31)
where V ∈ C1,2([0, T ],Rd). By the definition of metric balls {Bδk} and strictly-decreasing properties, the infimum in (31) is attained when θ(K) = θ and the control agent α(θ(K)) is optimal in this control set.
Relation to Stochastic Maximum Principle (SMP). We consider an arbitrary control in the convex set K ∈ A with β ∈ K and the optimal control α(θ(K)). Let DJ |β = dd J(θ(K) + (β − α))| =0 be the Gâteaux derivative (this can be defined while the control set is a vector sub-space, A ⊂ L2). By the result of Pontryagin maximum principle, Theorem (4.12) in Carmona (2016a), one can obtain the inequality as follows:∥∥∥∥ ∂∂αH ∥∥∥∥ · (∨Mr Lαr ) ‖θ(K)− θ‖ ≥ ∥∥∥∥ ∂∂αH ∥∥∥∥ · ‖[α(t,Xt, θ(K))− β(t,Xt, θ)]‖ ≥ DJ |β ≥ 0 (32)
for t ∈ [0, T ] almost surely, where we defineH , H(t,Xt, Yt, Zt, αt) for the Hamiltonian systemH with adjoint variables Yt, Zt and define the arbitrary control β = α(·, ·, θ) ∈ A with some θ. The first inequality is satisfied due to the definition of Lipschitzian control agents. The optimality condition indicates converging upper-bound of DJ |β to 0. In our method, the optimality condition of the proposed learning framework is bounded by the Euclidean distance between θ(K) and θ in parameter space. Thus, the proposed framework poses a fundamentally different approach to interact with optimality conditions in SMP. As we define that the θ(K) is a local minimum of J with the inequality ‖θ(K)− θ‖ < δKθ , the gradient descent scheme that induces the tight radius {δkθ}k∈N+ assures optimality by the relation 0 ≤ DJ |β ≈ δKθ . Relation to the HJB equation. We consider the infinitesimal generator Lt of the non-homogeneous controlled Markov process Xt as Lαt f = 〈∇f, b(t, x, α)〉+ 12Tr [ σσT (t, x, α)∇2f ] . We show the important relation between the proposed MBcond loss with the HJB equation as follows:
∂J ∂t (t, x) + Lα(θ(K))t J(t, x) + l(t, x)︸ ︷︷ ︸
Non-linear Feynman-Kac, MBcond loss
= 0 = ∂V
∂t (t, x) + inf α∈A [Lαt V (t, x) + l(t, x)]︸ ︷︷ ︸
HJB equation, exact solution
(33)
Equivalence Relation
In the left-hand side of (33), the PDE formula is directly consequence of the non-linear Feynman-Kac theorem that we derive in (24). The distinct point is that control agents are obtained by the gradient descent of the MBcond loss with BSDE (i.e., Zt). Note that, as shown in (31), θ(K) is actually an optimal control. This means that, without heavy calculations to solve the PDEs, the gradient descent algorithm also assures optimality of control agents in the proposed control set A.
In contrast, the HJB equation in the right-hand side states that the optimal control agent can be obtained by solving the second-order parabolic formula and the infimum is taken by considering algebraic properties of candidates for the exact solution. If the solution to HJBE exists in the control set A, the PDE in the left hand side of (33) approximates the solution to the HJB. Overall, we argue that the MFBcond loss can provide a novel deep learning-based paradigm to adopt/solve the conventional stochastic optimal control problem in a feasible way (i.e., well-defined loss functions with the gradient descent scheme).

heading: A.5 DATA PREPOSSESSING
text: PhysioNet dataset, Silva et al. (2012), contains overall 8000 multivariate time series obtained for the first 48 hours of a different patient’s admission to intensive care unit (ICU). Each patient has a set of 35 various clinical features. We normalized features of all patients in the dataset to have zero mean and unit variance. We used a half of time-series as the training dataset and the remaining parts as the test dataset.
Speech Commands dataset, Warden (2018), consists of one-second audio records of various spoken words such as “Yes”, “No”, “Up”, and “Down”. Since there were nearly 100,000 record samples, we sub-sampled the dataset due to the dimensionality of training instances on two conflicting classes (i.e., “Right” and “Left”). Overall, 6950 time-series records were selected, where 80% were used as training dataset and the remaining parts as the test dataset. We pre-processed these time series by computing Mel-frequency cepstrum coefficients from the audio signal, so that each time series was spaced with 65 and 54 channels. Then, we normalized each channel of all signals in the dataset to have zero mean and unit variance.
Beijing Air-Quality dataset, Zhang et al. (2017), consists of multi-year recordings of air quality data across different locations in Beijing. Each sample contains 6-dimensional time series features of PM2.5, PM10, SO2, NO2, CO, and O3, which are recorded per hour. We segmented data to have the length of 48 and normalized each feature of all data in the dataset to have zero mean and unit variance.
S&P-500 Stock Market dataset consists of stock market data with 6-dimensional feature vectors (i.e., [High, Low, Open, Close, Volume, Adj Close]). For the complete data acquisitions, we excluded enterprises with incomplete recordings during sampling duration, thus total 381 enterprises are selected. The time-series are sampled every 30-min with T = 48 temporal states. Similar to Speech commands dataset, we used first 80% of temporal states to train the model and the remaining parts are used for prediction task.

heading: A.6 EXPERIMENTS DETAILS
text: Different SDE candidates for CSDE. Owing to the abstract form of the proposed CSDE in (1), various types of drift and diffusion functions (i.e., b and σ) can be selected according to different applications. In Table 4, we enumerate candidate functions. In the experiments, we adopted two models: Vanilla and Mckean-Vlasov (MV) SDEs.
Hyperparameters. For the running and terminal costs (l and Ψ, respectively), we used the l2 distance, i.e., l(s, x) = ‖x− ys‖22 and Ψ(x) = ‖x− yT ‖
2. In all experiments, γ is set to 0.95. To estimate the gradient of the MBcond loss, we estimated numerical gradients with the auto-grad library in Pytorch (Paszke et al. (2019)).
Network Architecture for Neural Control Agents. Each control agent αi(t,Xt; θi) has an identical neural network architecture, which consists of linear layers and non-linear units. Figure 3 shows the detailed network architecture. Each agent takes concatenation of temporal/spatial tensors (t,Xt) as its input, where the temporal tensor t is transformed into new form t′ by the time inhomogeneous embedding layer. We followed the setting suggested in Park et al. (2021) for this embedding. After time embedding, concatenated tensor (t′, Xt) is fed into two Linear layers with non-
linearity units (i.e., LipSwish in Chen et al. (2019); Kidger et al. (2021)). Finally, the transformed tensors are split into the control terms for drift and diffusion functions. The diffusion functions are defined as non-degenerate types, where σi(t,Xt, αi) = Diag(zt) and zt is the output of the last linear layer. The latent dimension of each Linear layer was set to 128 in all experiments except for the prediction task with the Air Quality dataset (= 64). Thus, a total number of training parameters for single control agent αi is ≈ 11K. Simulation of CSDE and Temporal Privacy Function. Let T = t ∈ {tk}1≤k≤N for the pre-fixed time interval ∆t. We apply the Euler-Maruyama scheme to approximately simulate the proposed CSDE:
Xαt+∆t = X α t + M∑ i=1 wi(t)bi(t,Xαt , α i(t,Xαt ; θ i))∆t + M∑ i=1 wi(t)σi(t,Xαt , α i(t,Xαt ; θ i))Z, (34)
where Z , Z(0, √
∆tId) is a d-dimensional Gaussian random variable with zero mean and covariance√ ∆tId.
Analysis of Instability at Contact Points. In every time stamps t, drift and diffusion functions are controlled over neural control agents αi where we assume that t− and t+ are adjacent points of contact point t with infinitesimally small duration. The process At indicates drift integral term and σαs denotes the diffusion term in our forward CSDE dynamics. As shown in the inequality, the Markovian property is still preserved, and the magnitude of jumps are controlled by Lipschitzness of drift/diffusion functions. E [∥∥Xt− −Xt+∥∥2 |Ft−] ≤ E[‖At‖2 |Xt− ] + E [∫ t t− ‖σαs ‖ 2 ds|Xt− ] + E [∫ t+ t ∥∥σβs ∥∥2 ds|Xt−] . (35)
In a probabilistic point of view, the set of contact points may be regarded as measure-zero, and the probabilistic evaluation is not changed.
Figure 4 shows a particular example, where 14 temporal states (i.e., |T| = 14) with 4 control agents are considered. In the figure, the black line indicates the trajectory of time series, blue dots denote the observed data points, and shaded grey dots denote the missing data points. Each control agent takes 5 data points, where 2 temporal states are shared to other agents. In the experiments, the total number of temporal privacy functions are maximally set to M = |T|/2, where each control agent shares 2 points for smooth transitions of stochastic dynamics.

heading: A.7 ADDITIONAL EMPIRICAL STUDY
text: Effect of hyper-parameter γ. In Figure 5-(a), the effect of hyper-parameter γ is shown. Similar to Fig 2-(b), the results were obtained for the prediction task with the Air Quality dataset. Each red, black, and blue line indicates the test MSE for different γ ∈ [0.0, 0.95, 1.0] over 50 epochs. If the MFcond is deactivated during the training time i.e., γ = 0.0, only MBcond loss is utilized to train the proposed CSDE-TP, and the model produces poor results. As our inference procedure requires the model to train with multiple conditions, the obtained result seems obvious. If the MBcond loss is deactivated during the training time i.e., γ = 1.0, multi-conditioned information in backward dynamics Zαt are canceled, and the performance is decreased significantly i.e., 1.277→ 2.003. This clearly shows that MBcond loss boost the performance.
Effect of random stopping time. In Figure 5-(b), the effect of strategy to select is shown. If we select threshold as an uniform random variable ∼ U [s, T ] which is independent to Xt, then the network quickly falls into instability as shown in the red line of Figure 5-(b). This shows that the well-designed strategy for selecting threshold is crucial factor to stabilize the network learning landscape. Contrary to the random sampling strategy, our method defined in Algorithm 1 select half value of maximal MFcond loss in last learning steps as the threshold for random stopping time (i.e., 12 max lk−1 → k). As the threshold is always bounded above the maximal loss in the last steps, random stopping time at iteration k is decided in the time set:
τks ∈ {t : l(t, T αk s,t ) >
1 2 max l(t, T αk−1s,t )}, (36)
where τks denotes the stopping time at learning iteration k. If the network trains the MFcond loss so that lk , l(t, T αks,t )→ 0 as training proceeds k →∞, then it is clear that the stopping time vanishes τk→∞s ∈ ∅. Thus, the strategy in (36) is well-defined.

heading: A.8 DETAILED EXPLANATIONS OF MARKOV DYNAMIC PROGRAMMING WITH TEMPORAL PRIVACY
text: For the clear explanation of proposed Markov-DP-TP, let us consider the detailed example. we decompose sub-problem (B′) in (4) into another smaller sub-problems:
inf α(−r) E[J(u,Xαu )]︸ ︷︷ ︸ (B′) = inf β E
[∫ u′ u l(s,Xαs )ds ] ︸ ︷︷ ︸
(C)
+ inf β(−r′) E [J(u′, Xαu′)]︸ ︷︷ ︸ (C′) , (37)
where we set α(−r) = β, wr(s) = 1t≤s≤u. In this case, the problem (B) on interval [u, T ] is now decomposed into smaller sub-problems (C), (C′) on two intervals [u, u′] and [u′, T ]. Similarly to u in (4), another auxiliary time index u′ is considered here for additional problem (C). The corresponding new temporal privacy function wr′(s) = 1u≤s≤u′ is defined on the interval [u, u′].
By repeating temporal decomposition of original problem (A) M times, one can find the following hierarchical relations:
• P1). original problem, T = [t, · · · , T ], α, no temporal privacy • P2). Two sub-problems (B) + (B′) in (4),
Time set, T = [t, · · · , u, · · · , T ], control agents, α = [αr, α(−r)], temporal privacy functions = {wr}
• P3). Three sub-problems (B) + (C) + (C ′), Time set, T = [t, · · · , u, · · · , u′, · · · , T ], control agents, α = [αr, β, β(−r
′)], temporal privacy functions = {wr, wr′}
• P4). M sub-problems, (A) + (B) + (C) + . . . , Time set, T = [t, · · · , T−tM , · · · , r ∗ T−t M , · · · , T ],
α = [α1, α2, · · · , αr, · · · , αM ], temporal privacy functions = {w1, w2, · · · , wr, · · · , wM}
The role of u in (3) and (4) is replaced to u and u′ in (P3), and replaced to (r ∗ T−tM ) in (P4) in the table if the time interval is assumed to be regularly sampled. Similarly, the role of r in (3) and (4) is replaced to r and r′ in (P3).

heading: A.9 TOY EXAMPLE ON SYNTHETIC DATA
text: In this section, we conduct the reconstruction experiment on synthetic data to show the different behaviors and demonstrate the advantages of the proposed CSDE compared to previous methods.
Stochastic Trigonometric Data. In this experiment, we define the 100-dimensional stochastic process with composition of trigonometric functions (i.e., sin, cos) as follows:
Yt =
[ 1
2 sin(5πt+ Z1t) + 0.25 cos
( 13
5 πt+ Z2t
) + Z3 ] ∈ R100, (38)
where we assume t ∈ [0, 1.0] and the total number of temporal states are set to 48 (i.e., |T| = 48). In the definition of synthetic process Yt, both the period and amplitude are randomized with mean-zero Gaussian random variables (i.e., Z1 ∼ N (0, 1.0), Z2 ∼ N (0, 2.0), Z3 ∼ N (0, 12Id)). With the effect of Gaussian random variables, the process contains high volatility in both the spatial/temporal axes. We compare our method to the auto-regressive ODE-RNN Rubanova et al. (2019) model using the open-source code implemented by the authors. To observe the fundamental difference between ODE-RNN and CSDE-TP, we stop the training procedure when the estimated MSEs of both models attained the threshold (≤ .07). In Figures 6 and 7, the first axis of trigonometric data are visualized. The results of each model are indicated as the blue lines (i.e., Xt) and the synthetic trigonometric data is indicated as the red lines (i.e., Yt). The 95%-confidence regions (i.e., CR-95) of both the test and predicted time-series are shown as red and blue shaded regions, respectively.
ODE-RNN. Figure 6 shows the results of the ODE-RNN model. Although the ODE-RNN model attains relatively similar MSEs compared to the proposed model, there are two main issues in their model to be discussed.
1) It hardly captures the vertical perturbation of test data induced by Z3 and the obtained result produces a small variance at every temporal states.
2) It hardly captures the horizontal perturbation of test data induced by Z1, Z2, and the obtained result produces the temporally unmatched trajectories.
These phenomenons occurred due to the deterministic property of the ODE-RNN model, where the dynamical transition in their model is posed as the ODEs that cannot express the stochastic variation.
CSDE-TP. Figure 7 shows the result of the proposed CSDE-TP model and shows the advantages of adopting the SDE in modelling stochastic dynamics. Compared to the results of the ODE-RNN, the proposed method accurately captures both the vertical/horizontal perturbations and recover the 95% confidence region. It is clear that our CSDE-TP delicately expresses the complex volatility of stochastic trajectories.
Discussions. As aforementioned in Section 4.3, experimental results on synthetic stochastic data show that the MSE is not the best metric to train/evaluate the time-series models if there exists the
high volatility in the dataset. In this case, distributional metrics such as MMD and Wasserstein distance can be good substitutes for training/evaluating stochastic data.

heading: A.10 FUTURE WORK
text: We plan to extend the proposed CSDE model to a general controlled Markov Itô-Lévy jump diffusion model (Øksendal & Sulem (2007)) to delicately express the complex time-sequential data. For example, the proposed CSDE can be generalized to the Markov Itô-Lévy jump diffusion of the following form:
dXαt = b(t,X α t , α(θ))dt+ σ(t,X α t , α(θ))dWt +
∫ Γ(t, Z)N(dt, dZ), (39)
where N(t, z) = ∑
0<s≤t Xz∈U (ηs − ηs−), and Poisson random measure ηt. As the previous work in Jia & Benson (2019) show the effectiveness of the jump process in modelling complex discontinuous dynamics, we believe this generalization will produce the comparable results and broaden our understanding in modelling dynamical systems for time-series data.

Review summaries:
This paper presents a framework for learning a stochastic dynamics from observed trajectories. As opposed to conventional strategies based on recurrent neural networks, the procedure learns a model of the stochastic dynamics in real space. Secondly, in contrast to recent approaches based on neural ODEs / SDEs, the training is temporally localized through windowing functions that the authors refer to as temporal privacy. The main technical innovations involve the development of new loss functions that make learning with temporal privacy tractable. 

The authors propose a novel method to model stochastic dynamic of complex time series. They build connections between neural SDE and stochastic optimal control theory. By using the proposed MFcond and MBcond losses, they can train control agents to learn dynamics of time series more accurately than existing methods. 

This paper proposes to model a time series directly with a neural controlled stochastic differential equation with multiple control agents. It introduces a concept of temporal privacy, which defines the attention of each control agent to be a certain time interval. Then the authors introduce Markov dynamic programming to efficiently minimize a loss function defined in terms of trajectory of the stochastic differential equations. In order to make the proposed temporal privacy work, the authors proposed two loss function: 1) MFcond minimizes the incoherence between the future estimates starting at any intermediate time stamp, 2) MBcond is to ensure the theoretical optimality of the control agents. 

This paper presents a new approach to train stochastic dynamical systems using a set of tools from stochastic optimal control theory. The methodology is built upon _controlled stochastic differential equations_ with multiple control agents that modulate the dynamical evolution. Each agent is deliberately chosen to be "active" in a certain time interval, which leads to so-called "temporally private" dynamics that allow for dynamic programming based optimization. Since the temporal segmentation of optimization objective requires arbitrary intermediate state $X_s$ with $0<s<T$ as initial values, the Markov forward condition (MFcond), which computes some sort of a mean value for future states, is proposed. To further stabilize the training and compensate for the methodology's lacking in theoretical optimality guarantees, the authors augment the original loss with the Markov backward conditional (MBcond) loss. The method is shown to outperform competing continuous-time methods on standard benchmarks.