Year: 2021

Abstract:
Combinations of neural ODEs with recurrent neural networks (RNN), like GRUODE-Bayes or ODE-RNN are well suited to model irregularly observed time series. While those models outperform existing discrete-time approaches, no theoretical guarantees for their predictive capabilities are available. Assuming that the irregularly-sampled time series data originates from a continuous stochastic process, the L-optimal online prediction is the conditional expectation given the currently available information. We introduce the Neural Jump ODE (NJ-ODE) that provides a data-driven approach to learn, continuously in time, the conditional expectation of a stochastic process. Our approach models the conditional expectation between two observations with a neural ODE and jumps whenever a new observation is made. We define a novel training framework, which allows us to prove theoretical guarantees for the first time. In particular, we show that the output of our model converges to the L-optimal prediction. This can be interpreted as solution to a special filtering problem. We provide experiments showing that the theoretical results also hold empirically. Moreover, we experimentally show that our model outperforms the baselines in more complex learning tasks and give comparisons on real-world datasets.

heading: 1 INTRODUCTION
text: Stochastic processes are widely used in many fields to model time series that exhibit a random behaviour. In this work, we focus on processes that can be expressed as solutions of stochastic differential equations (SDE) of the form
dXt = µ(t,Xt)dt+ σ(t,Xt)dWt ,
with certain assumptions on the drift µ and the diffusion σ. With respect to the L2-norm, the best prediction of a future value of the process is provided by the conditional expectation given the current value. If the drift and diffusion are known or a good estimation is available, the conditional expectation can be approximated by a Monte Carlo (MC) simulation. However, since µ and σ are usually unknown, this approach strongly depends on the assumptions made on their parametric form. A more flexible approach is given by neural SDEs, where the drift µ and diffusion σ are modelled by neural networks (Tzen & Raginsky, 2019; Li et al., 2020; Jia & Benson, 2019). Nevertheless, modelling the diffusion can be avoided if one is only interested in forecasting the behaviour instead of sampling new paths.
An alternative widely used approach is to use Recurrent Neural Networks (RNN), where a neural network dynamically updates a latent variable with the observations of a discrete input time-series. RNNs are successfully applied to tasks for which time-series are regularly sampled, as for example speech or text recognition. However, often observations are irregularly observed in time. The standard approach of dividing the time-line into equally-sized intervals and imputing or aggregating observations might lead to a significant loss of information (Rubanova et al., 2019). Frameworks that overcome this issue are the GRU-ODE-Bayes (Brouwer et al., 2019) and the ODE-RNN (Rubanova et al., 2019), which combine a RNN with a neural ODE (Chen et al., 2018). In standard RNNs, the
hidden state is updated at each observation and constant in between. Conversely, in the GRU-ODEBayes and ODE-RNN framework, a neural ODE is trained to model the continuous evolution of the hidden state of the RNN between two observations. While GRU-ODE-Bayes and ODE-RNN both provide convincing empirical results, they lack thorough theoretical guarantees.
Contribution. In this paper, we introduce a mathematical framework to precisely describe the problem statement of online prediction and filtering of a stochastic process with temporal irregular observations. Based on this rigorous mathematical description, we introduce the Neural Jump ODE (NJ-ODE). The model architecture is very similar to the one of GRU-ODE-Bayes and ODE-RNN, however we introduce a novel training framework, which in contrast to them allows us to prove convergence guarantees for the first time. Moreover, we demonstrate empirically the capabilities of our model.
Precise problem formulation. We emphasize that a precise definition of all ingredients is needed, to be able to show theoretical convergence guarantees, which is the main purpose of this work. Since the objects of interest are stochastic processes, we use tools from probability theory and stochastic calculus. To make the paper more readable and comprehensible also for readers without background in these fields, the precise formulations and demonstrations of all claims are given in the appendix, while the main part of the paper focuses on giving well understandable heuristics.

heading: 2 PROBLEM STATEMENT
text: The problem we consider in this work, is the online forecasting of temporal data. We assume that we make observations of a Markovian stochastic process described by the stochastic differential equation (SDE)
dXt = µ(t,Xt)dt+ σ(t,Xt)dWt , (1)
at irregularly-sampled time points. Between those observation times, we want to predict the stochastic process, based only on the observations that we made previously in time, excluding the possibility to interpolate observations. Due to the Markov property, only the last observation is needed for an optimal prediction. Hence, after each observation we extrapolate the current observation into the future until we make the next observation. The time at which the next observation will be made is random and assumed to be independent of the stochastic process itself.
More precisely, we suppose to have a training set of N independent realisations of the RdX - dimensional stochastic process X defined in (1). Each realisation j is observed at nj random observation times t(j)1 , . . . , t (j) nj ∈ [0, T ] with values x (j) 1 , . . . , x (j) nj ∈ RdX . We assume that all coordinates of the vector x(j)i are observed. We are interested in forecasting how a new independent realization evolves in time, such that our predictions of X minimize the expected squared distance (L2-metric) to the true unknown path. The optimal prediction, i.e. the L2-minimizer, is the conditional expectation. Given that the value of the new realization at time t is xt, we are therefore interested in estimating the function
f(xt, t, s) := E[Xt+s|Xt = xt], s ≥ 0, (2)
which is the L2-optimal prediction until the next observation is made. To learn an approximation f̂ of f we make use of the N realisations of the training set. After training, f̂ is applied to the new realization. Hence, this can be interpreted as a special type of filtering problem. The following example illustrates the considered problem.
Example. A complicated to measure vital parameter of patients in a hospital is measured multiple times during the first 48 hours of their stay. For each patient, this happens at different times depending on the resources, hence the observation dates are irregular and exhibit some randomness. Patient 1 has n1 = 4 measurements at hours (t (1) 1 , t (1) 2 , t (1) 3 , t (1) 4 ) = (1, 14, 27, 34) where the values (x (1) 1 , x (1) 2 , x (1) 3 , x (1) 4 ) = (0.74, 0.65, 0.78, 0.81) are measured. Patient 2 only has n2 = 2 measurements at hours (t(2)1 , t (2) 2 ) = (3, 28) where the values (x (2) 1 , x (2) 2 ) = (0.56, 0.63) are measured. Similarly, the j-th patient has nj measurements at times (t (j) 1 , . . . , t (j) nj ) and has the measured values (x (j) 1 , . . . , x (j) nj ). Based on this data, we want to forecast the vital parameter of new patients coming to the hospital. In particular, for a patient with measured values x1 at time t1, we want to predict what
his values will likely be at any time t1 + s > t1. Importantly, we do not only focus on predicting the value at some t2 > t1, but we want to know the entire evolution of the value.

heading: 3 BACKGROUND
text: Recurrent Neural Network. The input to a RNN is a discrete time series of observations {x1, · · · , xn}. At each observation time ti+1, a neural network, the RNNCell, updates the latent variable h using the previous latent variable hi and the input xi+1 as
hi+1 := RNNCell(hi, xi+1).
Neural Ordinary Differential Equation. Neural ODEs (Chen et al., 2018) are a family of continuous-time models defining a latent variable ht := h(t) to be the solution to an ODE initial-value problem
ht := h0 + ∫ t t0 f(hs, s, θ)ds, t ≥ t0, (3)
where f(·, ·, θ) = fθ is a neural network with weights θ. Therefore, the latent variables can be updated continuously by solving this ODE (3). We can emphasize the dependence of ht on a numerical ODE solver by rewriting (3) as
ht := ODESolve(fθ, h0, (t0, t)) . (4)
ODE-RNN. ODE-RNN (Rubanova et al., 2019) is a mixture of a RNN and a neural ODE. In contrast to a standard RNN, we are not only interested in an output at the observation times ti, but also in between those times. In particular, we want to have an output stream that is generated continuously in time. This is achieved by using a neural ODE to model the latent dynamics between two observation times, i.e. for ti−1 < t < ti the latent variable is defined as in (3) and (4), with h0 and t0 replaced by hi−1 and ti−1. At the next observation time ti, the latent variable is updated by a RNN with the new observation xi. Fixing h0, the entire latent process can be computed by iteratively solving an ODE followed by applying a RNN. Rubanova et al. (2019) write this as{
h′i := ODESolve(fθ, hi−1, (ti−1, ti)) hi := RNNCell(h ′ i, xi) .
(5)
GRU-ODE-Bayes. The model architecture describing the latent variable in GRU-ODE-Bayes (Brouwer et al., 2019) is defined as a special case of the ODE-RNN architecture. In particular, a gated recurrent unit (GRU) is used for the RNN-cell and a continuous version of the GRU for the neural ODE fθ. Therefore, we focus on explaining the difference between our model architecture and the ODE-RNN architecture, in the following section.

heading: 4 PROPOSED METHOD – NEURAL JUMP ODE
text: Markovian paths. Our assumptions on the stochastic process X imply that it is a Markov process. In particular, the optimal prediction of a future state of X only depends on the current state rather than on the full history. Hence, the previous values do not provide any additional useful information for the prediction.
JumpNN instead of RNN. Using a neural ODE between two observation has the advantage that it allows to continuously model the hidden state between two observations. But since the underlying process is Markov, there is no need to use a RNN-cell to model the updates of the hidden state at each new observation. Instead, whenever a new observation is made, the new hidden state can solely be defined from this observation. Therefore, we replace the RNN-cell used in ODE-RNN by a standard neural network mapping the observation to the hidden state. We call this the jumpNN which can be interpreted as an encoder map. Compared to ODE-RNN, this architecture is easier to train.
Last observation and time increment as additional inputs for the neural ODE. The neural network fθ used in the neural ODE takes two arguments as inputs, the hidden state ht and the current time t. However, our theoretical problem analysis suggests, that instead of t the last observation time ti−1 and the time increment t− ti−1 should be used. Additionally the last observation xi−1 should also be part of the input.
NJ-ODE. Combining the ODE-RNN architecture (5) with the previous considerations, we introduce the modified architecture of Neural Jump ODE (NJ-ODE){
h′i := ODESolve(fθ, (hi−1, xi−1, ti−1, t− ti−1), (ti−1, ti)) hi := jumpNN(xi) .
(6)
An implementable version of this method is presented in the Algorithm 1. A neural ODE fθ transforms the hidden state between observations, and the hidden state jumps according to jumpNN when a new observation is available. The outputNN, a standard neural network, maps any hidden state ht to the output yt. To implement the continuous-in-time ODE evaluation, a discretization scheme is provided by the inner loop. In the training process, the weights of all three neural networks, jumpNN, the neural ODE fθ and outputNN are optimized.
Algorithm 1 The NJ-ODE. A small step size ∆t is fixed and we denote tn+1 := T . Input: Data points with timestamps {(xi, ti)}i=0...n, for i = 0 to n do hti = jumpNN(xi) . Update hidden state given next observation xi yti = outputNN(hti) . compute output s← ti while s+ ∆t ≤ ti+1 do hs+∆t = ODESolve(fθ, hs, xi, ti, s− ti, (s, s+ ∆t)) . get next hidden state ys+∆t = outputNN(hs+∆t) . compute output s← s+ ∆t
end while end for
Objective function. Our goal is to train the NJ-ODE model such that its output approximates the conditional expectation (2), which is the optimal prediction of the target process X with respect to the L2-norm. Therefore, we define a new objective function, with which we can prove convergence. Let yi− denote the output of the NJ-ODE at ti before the jump and yi the output at ti after the jump. Note that the outputs depend on parameters θ and the previously observed xi which are inputs to the model. Then the objective function is defined as
Φ̂N (θ) := 1
N N∑ j=1︸ ︷︷ ︸
paths
1
nj nj∑ i=1︸ ︷︷ ︸
dates
( |x(j)i − y
(j) i |︸ ︷︷ ︸
jump part at observations
+ |y(j)i − y (j) i− |︸ ︷︷ ︸
continuous part between two observations
)2 . (7)
We give an intuitive explanation for this definition. The “jump part” of the loss function forces the jumpNN to produce good updates based on new observations, while the other part forces the jump size to be small in (the empirical) L2-norm. Since the conditional expectation minimizes the jump size with respect to the L2-norm, this forces the neural ODE fθ to continuously transform the hidden state such that the output approximates the conditional expectation. Moreover, both parts of the loss function force the outputNN to reasonably transform the hidden state ht to the output yt.

heading: 5 MAIN RESULT – THEORETICAL CONVERGENCE GUARANTEE
text: In the following we informally state our main result. To formally state and prove the theorem, precise definitions of all ingredients are needed. This analysis is provided in the appendix where the following results is stated in Theorem E.2 and Theorem E.13. Theorem 5.1 (informal). We assume that for each number of paths N and for every size of the neural networks M , their weights are chosen optimally, as to minimize Φ̂N (θ). Then, if N and M tend to infinity, the output of NJ-ODE converges in mean (L1-convergence) to the conditional exception of the stochastic process X given the current information.
An intuitive explanation for this theorem was given with the definition of the objective function. In this result, the focus lies on the convergence analysis under the assumption that optimal weights are found. In the Appendix we discuss why this assumption is not restrictive.

heading: 6 EXPERIMENTS
text: For further details and results for all experiments see Appendix F.

heading: 6.1 TRAINING ON SYNTHETIC DATASETS
text: Evaluation metric. For synthetic datasets where an analytic formula for the conditional expectation exists, we can evaluate the distance of the model output to the target process (2). We use a sampling time grid with equidistant step size ∆t := TK , K ∈ N, on [0, T ]. On this grid, we compare for path j at time t, the true conditional expectation x̂(j)t with the predicted conditional expectation (the model output) y(j)t . For N2 test samples, the evaluation metric is defined as
eval(x̂, y) := 1
N2 N2∑ j=1 1 K + 1 K∑ i=0 ( x̂ (j) i∆t − y(j)i∆t )2 . (8)
Black-Scholes, Ornstein-Uhlenbeck and Heston. We test our algorithm on three scalar stochastic models, Black-Scholes, Ornstein-Uhlenbeck and Heston, with fixed parameters. For each model, we generate a dataset by sampling N = 20′000 paths on the time interval [0, 1] using the Euler-scheme with 100 time steps. Independently for each path, on average 10% of the grid points are randomly chosen as observation times. The NJ-ODE is trained on 80% of the data. On the remaining 20% the model is tested, by comparing the loss function (31) computed with the NJ-ODE to the loss function computed with the true conditional expectation (Figure 2). During training, the relative difference becomes very small, hence the true conditional expectation is nearly replicated.

heading: 6.2 FURTHER SYNTHETIC DATASETS
text: Heston model without Feller condition. We also train our model on a Heston dataset for which the Feller condition is not satisfied. As explained by Andersen (2007); Jean-François et al. (2015) this is a more delicate situation. We see that our model is very robust, since even in this critical case it learns to replicate the true conditional expectation process. In the Heston model, the variance of the stochastic process Xt is a stochastic process, vt. Here, we train our model to predict both processes at the same time. The training on this 2-dimensional dataset is successful as can be seen in Figure 3. The minimal evaluation metric after 200 epochs is 0.0983.
Dataset with changing regime. In this experiment we test how well our model can deal with stochastic processes, that undergo an (abrupt) change of regime at a certain point in time. Many real world time series might exhibit such a change of regime. Some examples are listed below.
• Longitudinal patient health recordings might experience changes depending on seasonal or longerterm influences, as for example due to the seasonal flue or currently the Covid-19 pandemic. • In many regions climate data has strong seasonal dependencies, that can lead to relatively abrupt changes as for example when the weather changes from dry to rain season.
• A stock market that suddenly changes from a bullish to a bearish market, for example due to a macro-economic event. An example for this would be the start of the Covid-19 crisis in the first quarter of 2020.
We test a change of regime by combining two synthetic datasets. On the first half of the time interval [0, 0.5] we use the Ornstein-Uhlenbeck and on the second half [0.5, 1] the Black-Scholes model. In Figure 4 we see that our model correctly learns the change of regime. The minimal evaluation metric after 200 epochs is 0.0463.
Dataset with explicit time dependence. Many real world datasets, have an explicit time dependence, i.e. the drift and diffusion of (9) explicitly depend on t. Examples are all datasets that have a certain periodicity, as for example weather data (seasonal and daily periodicity), intraday periodicity of stock prices (Andersen & Bollerslev, 1997) or prices of certain seasonal goods. We incorporate an explicit time dependence into the Black-Scholes dataset, by replacing the drift constant µ with the time dependent constant α2 (sin(βt) + 1), for α, β > 0. In Figure 5 we see that the model learns to adapt to the time-dependent coefficients. The minimal evaluation metric after 100 epochs is 0.0215 (β = 2π) and 0.02805 (β = 4π) respectively.

heading: 6.3 EMPIRICAL CONVERGENCE STUDY
text: We confirm the theoretical results of Theorem 5.1 by an empirical convergence study for growing numbers of training samples N1 and network size M , where the performance is measured by the evaluation metric (8). For each combination of the number of training samples N1 and the neural network size M , the NJ-ODE is trained 5 times on the Black-Scholes, Ornstein-Uhlenbeck and
Heston datasets. The means of the evaluation metric over the 5 trials with their standard deviations are shown in Figure 6 for Black-Scholes. Already 200 training samples lead to a small evaluation metric if the network size is big enough, suggesting that the Monte Carlo approximation of the loss function is good, already with a few samples.

heading: 6.4 COMPARISON TO GRU-ODE-BAYES ON SYNTHETIC DATA
text: On the Black-Scholes, Ornstein-Uhlenbeck and Heston datasets, we compare our model to GRUODE-Bayes (Brouwer et al., 2019), which is, to the best of our knowledge, the neural network based method with the most similar task to ours. Results of our comparison are shown in Table 6.4. On the Black-Scholes and Ornstein-Uhlenbeck dataset, our model performs similarly to GRU-ODE-Bayes. However, on the more complicated Heston dataset, the training of GRU-ODE-Bayes is unstable and does not converge. On the other hand, our model converges during training. Although the value of the evaluation metric is much higher than for Black-Scholes and Ornstein-Uhlebeck, the resulting model output is still meaningful, which is not the case for GRU-ODE-Bayes. Hence we conclude, that GRU-ODE-Bayes cannot be applied reasonably for the Heston dataset, while our method works as the theoretical results suggest.

heading: 6.5 REAL WORLD DATASETS WITH INCOMPLETE OBSERVATIONS
text: Self-imputation for incomplete observations. Until now we assumed that at each observation time all coordinates of the stochastic process X are observed. However, in many real world applications, the observations are incomplete, i.e. only for some of the coordinates an observation is available. To deal with such incomplete observations, we propose the following self-imputation method. Whenever
an observation xi ∈ RdX is made the mask mi ∈ {0, 1}dX tells whether the k-th coordinate of the vector xi is observed (mki = 1) or not (m k i = 0). Then we use the prediction yi− of our model to generate the imputed observation vector
x̃i := mt xi + (1dX −mi) yi−,
where is the element-wise multiplication (Hadamar product) and 1dX ∈ RdX is the one-vector. Instead of xi we use (x̃i,mi) as an input for the jump part jumpNN. The intuition behind this definition is the following. In the one dimensional case, if we do not make an observation, but input yt− as if it was an observation, we expect that this does not change the output ys for s ≥ t. From this point of view, yt− does not provide any additional information for the model. Similarly, we expect that imputing yt− for unobserved coordinates does not provide any information about this coordinate to the model. However, since the model might learn to transfer the information about an observed coordinate to an unobserved one, we extend the input to also include the information which coordinates were observed. For the ODE part we use yi, the prediction after processing the input (x̃i,mi), as input instead of xi. Here, the intuition is that if the model learns how to best use the incomplete observation, i.e. if the jump is good, then this is the best approximation for xi. Our objective function is adjusted by multiplying each term in the sum with the mask.
Climate forecast. We compare our model to GRU-ODE-Bayes on the USHCN daily dataset (Menne et al., 2016), using the same experimental setting as was used by Brouwer et al. (2019). We train a small (S) and a large (L) version of NJ-ODE with different total number of parameters. The validation set was used for early stopping after the first 100 epochs, where we trained for a total of 200 epochs. We see in Table 6.5 that our small version performs slightly worse than GRU-ODE-Bayes while our large version slightly outperforms it.
Physionet. We compare our model to the latent ODE on their extrapolation task on the PhysioNet Challenge 2012 dataset (Goldberger et al., 2000), using the same experimental setting as was used by Rubanova et al. (2019). The mean and standard deviation over 5 runs starting at different random initializations is reported for our model. We see in Table 6.5 that our model outperforms the latent ODE models although having only about a seventh of the trainable weights.

heading: 7 RELATED WORK
text: Stochastic filtering theory. Our main theorem is similar to the theoretical results of neural filtering (Lo, 2009), which is a neural network approach to stochastic filtering theory in discrete time. In stochastic filtering, potentially incomplete and noisy observations of X are available continuously in time. This observation process Y is usually described by the dynamics dYt = h(Xt)dt + dBt, where h is a measurable function and B is a Brownian motion. Stochastic filtering then estimates the conditional law of Xt given the noisy observations (Ys)0≤s≤t (Bain & Crisan, 2008, Def. 3.2) and therefore can provide conditional expectations. Comparably, we directly compute the conditional expectation given the last observation. Similar to our assumptions, in the neural filtering approach of GRU-ODE-Bayes (Brouwer et al., 2019) and (Ryder et al., 2018) observations are only available at irregular discrete time points. They approximate the conditional law of X given the last observation by a Gaussian distribution and learn its mean and variance parameters. In particular, the conditional expectation is then given by this mean parameter. In contrast, we do not make normality assumptions about the conditional distribution and we theoretically prove convergence to the true conditional expectation.
Neural ODEs with jumps. Except for GRU-ODE-Bayes (Brouwer et al., 2019) and ODE-RNN (Rubanova et al., 2019) another work studying a neural ODE with jumps is Neural Jump SDE (NJSDE) (Jia & Benson, 2019). Similar to the NJ-ODE framework (29), the latent process of NJSDE is described by a neural ODE with jumps at random times. This model is used to describe hybrid systems which evolve continuously in time but may also be interrupted by stochastic events. In contrast to that, we model the conditional expectation of a continuous stochastic process.

heading: 8 CONCLUSION
text: We presented the Neural Jump ODE, a data-driven framework for modelling the conditional expectation of a stochastic process given the previous observations. We introduced a rigorous mathematical description of our model and more generally for the class of neural ODE based models. Moreover, for the first time we provided theoretical guarantees for a model falling in this category. We evaluated our model empirically on six synthetic and two real world datasets. In comparison to the baselines GRU-ODE-Bayes and latent ODE, we achieved better results especially on complex datasets.

heading: ACKNOWLEDGEMENT
text: The authors thank Andrew Allan, Robert A. Crowell, Anastasis Kratsios and Pierre Ruyssen for helpful discussions, providing references and insights. Moreover, the authors thank the reviewers for their thoughtful feedback that contributed to significantly improve the paper. The authors gratefully acknowledge financial support coming from the Swiss National Science Foundation (SNF) under grant 179114.

heading: APPENDIX

heading: A SETUP

heading: A.1 STOCHASTIC PROCESS X
text: In this section we rigorously describe the process X and give the assumptions which are needed in order to derive the convergence results. Let dX , dW ∈ N and T > 0 be the fixed time horizon. Consider a filtered probability space (Ω,F ,F := {Ft}0≤t≤T ,P), on which an adapted dW -dimensional Brownian motion {Wt}t∈[0,T ] is defined. We define the stochastic process1 X := (Xt)t∈[0,T ] as the solution of the stochastic differential equation (SDE)
dXt = µ(t,Xt)dt+ σ(t,Xt)dWt , (9)
for all 0 ≤ t ≤ T , where X0 = x ∈ RdX is the starting point and the measurable functions µ : [0, T ] × RdX → RdX and σ : [0, T ] × RdX → RdX×dW are the drift and the diffusion respectively. We impose the following assumptions:
• X is continuous and square integrable, i.e. for P-a.e. ω ∈ Ω the map t 7→ Xt(ω) is continuous and E[X2t ] <∞ for every t ∈ [0, T ].
• µ and σ are both globally Lipschitz continuous in their second component, i.e. for ϕ ∈ {µ, σ} there exists a constant M̃ > 0 such that for all t ∈ [0, T ]
|ϕ(t, x)− ϕ(t, y)|2 ≤ M̃ |x− y|2 and |ϕ(t, x)|2 ≤ (1 + |x|2)M̃. (10)
In particular, their growth is at most linear in the second component.
• µ is bounded and continuous in its first component (t) uniformly in its second component (x), i.e. for every t ∈ [0, T ] and ε > 0 there exists a δ > 0 such that for all s ∈ [0, T ] with |t− s| < δ and all x ∈ RdX we have |µ(t, x)− µ(s, x)| < .
• σ is càdlàg (right-continuous with existing left-limit) in the first component and L2 integrable with respect to W , σ ∈ L2(W ), i.e.
E  dX∑ i=1 dW∑ j=1 ∫ T 0 sup x σi,j(x, t) 2d[W j ,W j ]t  = ∫ T 0 |sup x σ(x, t)|2F dt <∞, (11)
where |·|F denotes the Frobenius matrix norm. This is in particular implied if σ is bounded.

heading: A.2 RANDOM OBSERVATION DATES
text: In this section we describe how we model the observation dates. We want to treat irregularly observed time-series, i.e. without a fixed number of observations and not necessarily equally distributed. Therefore, we suppose that we have a possibly random number of n observations and that those observations are made at random times t1, . . . , tn. For simplicity we make the assumption that the observation dates are independent of the stochastic process itself. In future work, this will be generalized to observation dates correlated with X , that could be modelled by a point process on the same probability space. In contrast, here we hard-code the independence assumption by considering a second probability space (Ω̃, F̃ , P̃), on which the random observation times of the stochastic process X are defined. More precisely, we assume that:
• n : Ω̃→ N≥0 is a random variable with EP̃[n] <∞, the random number of observations,
• ti : Ω̃→ [0, T ] for 0 ≤ i ≤ n are sorted 2 random variables, the random observation times.
We denote the joint pushforward measure of n and {ti}0≤i≤n as P̃t := (n, t0, . . . , tn)#P̃. The random variable n can but does not have to be unbounded. If it is bounded, we define K :=
1A stochastic process is a collection of random variables Xt : Ω→ RdX , ω 7→ Xt(ω) for 0 ≤ t ≤ T . 2For all ω̃ ∈ Ω̃, 0 = t0 < t1(ω̃) < · · · < tn(ω̃)(ω̃) ≤ T .
max { k ∈ N | P̃(n ≥ k) > 0 } to be the maximal value of n, which otherwise is infinity. We use the
notation B([0, T ]) for the Borel σ-algebra of the set [0, T ]. Then we define for each 1 ≤ k ≤ K
λk : B([0, T ])→ [0, 1], B 7→ λk(B) := P̃(n≥k,(tk−)∈B)P̃(n≥k) ,
which is a probability measure on the time interval as shown in the lemma below. Here, (tk−) means the left-point of tk, for example, if t0 = 0, t1 = 1 we have t0 ∈ [0, 1), t1 /∈ [0, 1) but t0− /∈ [0, 1), t1− ∈ [0, 1). Lemma A.1. For 1 ≤ k < K + 1 the map λk defines a probability measure.
Proof. First we see that λk([0, T ]) = 1 since tk maps to [0, T ]. Furthermore, for any disjoint sets Bi ∈ B([0, T ]), i ∈ N, we have {n ≥ k, tk− ∈ ∪̇i≥1Bi} = ∪̇i≥1{n ≥ k, tk− ∈ Bi} and these sets are F̃-measurable, since they are defined through pre-images of random variables. Therefore, the additivity of P̃ implies that λk (∪̇i≥1Bi) = ∑ i≥1 λk (Bi).
Moreover, we define τ as the time of the last observation before a certain time t,
τ : [0, T ]× Ω̃→ [0, T ], (t, ω̃) 7→ τ(t, ω̃) := max{ti(ω̃)|0 ≤ i ≤ n(ω̃), ti(ω̃) ≤ t}.
Example A.2. We give two examples how the random observation dates could be defined.
• Let (Nt)t∈[0,T ] be a homogeneous Poisson point process with rate r > 0. Hence, N(0) = 0, t 7→ N(t) is constant except for jumps of size 1 at discrete random times and for any fixed t ∈ [0, T ], N(t) is Poisson distributed, i.e.
P̃[N(t) = k] = (rt) k k! e −rt for k ∈ N.
Then the number of observations is defined as n := N(T ) and satisfies EP̃[n] = rT . The observation dates t1, . . . , tn are defined as the discontinuity times of N , i.e. the times where N increases by 1. In particular, for 0 ≤ k ≤ n, tk := inf{t|N(t) = k} ≤ T . Therefore, for 0 ≤ a ≤ b ≤ T we can rewrite
λk((a, b)) = P̃(N(a)<k,N(b)≥k)
P̃(N(T )≥k) .
• We define n by one of the following options:
– as a constant in N>0. – as a Binomial random variable, n ∼ Binom(p, nmax), for p ∈ (0, 1), nmax ∈ N>0. – as a Geometric random variable, n ∼ Geom(p), for p ∈ (0, 1). – as a Poisson random variable, n ∼ Poi(r), for r > 0.
and t1, . . . , tn are defined by choosing n uniform random variables on [0, T ] and sorting them.
A.3 INFORMATION σ-ALGEBRA
In this section, we define a mathematical tool, the σ-algebra, that is essential to the description of the conditional expectation. This object describes which information is available at a certain time t. In the following, we leave away ω̃ ∈ Ω̃ whenever the meaning is clear. We define the filtration of the currently available information A := (At)t∈[0,T ] by
At := σ(Xti |ti ≤ t),
where ti are the observation times and σ(·) denotes the generated σ-algebra. By the definition of τ we have At = Aτ(t) for all t ∈ [0, T ]. (Ω × Ω̃,F ⊗ F̃ ,F ⊗ F̃ ,P × P̃) is the filtered product probability space which, intuitively speaking, combines the randomness of the stochastic process with the randomness of the observations. Here, F ⊗ F̃ consists of the tensor-product σ-algebras (F ⊗ F̃)t := Ft ⊗ F̃ for t ∈ [0, T ]. As explained in the remark below, At can be identified with a sub-σ-algebra of (F ⊗ F̃)t.
Remark A.3. The σ-fields At depend on ω̃ ∈ Ω̃ as well. If we look at the product probability space (Ω × Ω̃,F ⊗ F̃ ,F ⊗ F̃ ,P × P̃), where F ⊗ F̃ consists of the tensor-product σ-algebras (F ⊗ F̃)t := Ft ⊗ F̃ for t ∈ [0, T ], then
σ(Xti |ti ≤ t) := σ ({{ (ω, ω̃) ∈ Ω× Ω̃ ∣∣Xti(ω̃)(ω) ∈ A,n(ω̃) ≥ i, ti(ω̃) ≤ t} ∣∣∣A ∈ B(RdX ), i ∈ N})
is a well defined sub-σ-algebra of (F ⊗F̃)t. Furthermore, we can recover the ω̃-wise defined version of At by intersecting each set in it with Ω× {ω̃} and subsequently projecting the intersection on its first component. We use the notation Ãt := Ãt(ω̃) = σ(Xti(ω̃)|ti(ω̃) ≤ t) to distinguish this ω̃-wise definition from the definition as sub-σ-algebra of the product space given above. However, Lemma B.3 implies that for our considerations, both versions of this σ-algebra have the same effect, therefore we will simply write At for both versions, by abuse of notation.

heading: B OPTIMAL APPROXIMATION X̂ OF THE STOCHASTIC PROCESS X
text: We are interested in the “best” approximation (or prediction) X̂t of the process X that one can make at any time t ∈ [0, T ], given the currently available information At. For us “best” refers to the L2(Ω× Ω̃,F⊗ F̃ ,P× P̃)-minimizer, therefore, this approximation is given by the conditional expectation. Indeed, if we define ∆ := {(t, r) ∈ [0, T ]2|t+ r ≤ T}, and the function
µ̃ : ∆× RdX → RdX , ((t, r), ξ) 7→ E [µ(t+ r,Xt+r)|Xt = ξ] ,
this is proven in the following proposition.
Proposition B.1. The optimal (i.e. L2-norm minimizing) A-adapted process in L2(Ω×Ω̃,F⊗F̃ ,P× P̃) approximating (Xt)t∈[0,T ] is given by3 X̂ := (X̂t)t∈[0,T ] with X̂t := EP×P̃[Xt|At]. Moreover, this process is unique up to (P× P̃)-null-sets. In addition we have (ω̃-wise for ω̃ ∈ Ω̃) that
X̂t = Xτ(t) + ∫ t τ(t) µ̃ ( τ(t), s− τ(t), Xτ(t) ) ds, (12)
implying that X̂ is càdlàg.
The first part of the result follows from the elementary Proposition below (which is proven for example in (Durrett, 2010, Thm. 5.1.8) for R-valued random variables and can easily be extended to Rd-valued random variables when using the 2-norm). Proposition B.2. Given a probability space (Ω,F ,P) and a sub-σ-algebra A ⊂ F , the orthogonal projection of X ∈ L2(Ω,F ,P) on L2(Ω,A,P) is given by X̂ := E[X|A]. In particular, for every Z ∈ L2(Ω,A,P) with P(Z 6= X̂) > 0 we have
E[|X − Z|22] = E[|X − X̂|22] + E[|Z − X̂|22] > E[|X − X̂|22].
Proof Proposition B.1. In our case this means, that the optimal A-adapted process in L2(Ω× Ω̃,F⊗ F̃ ,P× P̃) approximating X is given by (X̂t)t∈[0,T ] with X̂t := EP×P̃[Xt|At]. Here, At is meant as a sub-σ-algebra of Ft × F̃ . This process is unique up to (P× P̃)-null-sets. Moreover, the following lemma shows that it coincides ω̃-wise with EP[Xt|Ãt](ω̃), where Ãt = Ãt(ω̃) is defined in Remark A.3.
Lemma B.3. For P̃-almost-every ω̃ ∈ Ω̃ we have EP×P̃[Xt|At](ω̃) = EP[Xt|Ãt](ω̃) P-almostsurely.
3While we give a pointwise definition, (Cohen & Elliott, 2015, Theorem 7.6.5) allows to define X̂ directly as the optional projection. By (Cohen & Elliott, 2015, Remark 7.2.2) this implies that the process X̂ is progressively measurable, in particular, jointly measurable in t and ω×ω̃. However, as we show below, even from the pointwise definition, it follows that X̂ is càdlàg, hence optional (Cohen & Elliott, 2015, Theorem 7.2.7).
Proof. Otherwise we have by Fubini’s theorem, Proposition B.2 and an argument similar to the one in Lemma E.4
EP×P̃ [∣∣Xt − EP×P̃[Xt|At]∣∣22] = EP̃ [EP [∣∣Xt − EP×P̃[Xt|At]∣∣22]]
> EP̃ [ EP [∣∣∣Xt − EP[Xt|Ãt]∣∣∣2
2
]] ,
which is a contradiction to Proposition B.2.
This proves the first part of Proposition B.1. The second part of this Proposition, i.e. (12), should be understood ω̃-wise, for ω̃ ∈ Ω̃. This is justified by Lemma B.3 and derived below. In particular, for the remainder of this section, all statements are meant ω̃-wise.
With the assumption that µ and σ are Lipschitz, (Protter, 2005, Thm. 7, Chap. V) implies that a unique solution of (9) exists. Furthermore, this solution is a Markov process as soon as the starting point x is fixed (Protter, 2005, Thm. 32, Chap. V). Hence, one can define a transition function P (compare (Protter, 2005, Chap. V.6)) such that for all s < t and φ bounded and measurable,
Ps,t(Xs, φ) := EP[φ(Xt)|σ(Xs)] = EP[φ(Xt)|Fs].
We have that Xτ(s) is Aτ(s)-measurable and therefore, since As = Aτ(s) ⊂ Fτ(s),
Pτ(s),t(Xτ(s), φ) = EP[φ(Xt)|As]. (13)
By our additional assumption on σ it follows from (Protter, 2005, Lem. before Thm. 28, Chap. IV) that
Mt := ∫ t 0 σ(s,Xs)dWs, 0 ≤ t ≤ T,
is a square integrable martingale, since the Browian motion W is square integrable. In particular, for 0 ≤ s ≤ t ≤ T we have EP[ ∫ t s σ(r,Xr)dWr|Fs] = EP[Mt −Ms|Fs] = 0. Moreover, the same is true when conditioning on As 4. Using the martingale property of M , we have (ω̃-wise) for every t ∈ [0, T ]
X̂t = EP[(Xt −Xτ(t)) +Xτ(t)|Aτ(t)]
= Xτ(t) + EP [∫ t τ(t) µ(r,Xr)dr ∣∣∣Aτ(t) ] + EP [∫ t τ(t) σ(r,Xr)dWr ∣∣∣Aτ(t) ]
= Xτ(t) + ∫ t τ(t) EP [ µ(r,Xr)|Aτ(t) ] dr,
(14)
where we used Fubini’s Theorem (for conditional expectations) in the last step. This is justified because EP[ ∫ T 0 |µ(r,Xr)|2dr] < ∞ follows from µ being bounded. Let us define ∆ := {(t, r) ∈ [0, T ]2|t+ r ≤ T} and the function
µ̃ : ∆× RdX → RdX , ((t, r), ξ) 7→ Pt,t+r(Xt, µ) ∣∣ Xt=ξ = EP [µ(t+ r,Xt+r)|Xt = ξ] ,
then we can use (13) to rewrite (14) as
X̂t = Xτ(t) + ∫ t τ(t) µ̃ ( τ(t), s− τ(t), Xτ(t) ) ds. (15)
This proves the second part of Proposition B.1.
Proposition B.4. The function µ̃ is (jointly) continuous. 4To see this, we choose a localizing sequence (τn)n∈N such that Mτn is bounded by n (works since M is continuous). Then the Markov property implies that E[Mτnt −Mτns |As] = E[M τn t −Mτns |Fs] = 0. Since Mτnt n→∞−−−−→Mt P-a.s. and since this sequence is dominated by the integrable random variable 1+supr≤t|Mr|21 (by Doob’s inequality and square integrability ofM ), dominated convergence implies that E[Mt−Ms|As] = 0.
Proof of Proposition B.4. For any fixed s ∈ [0, T ], x ∈ RdX , we define
ζs,·(x) : [0, T ]× Ω→ RdX , (t, ω) 7→ ζs,t(x)(ω) to be the solution of the SDE
ζs,t(ξ) = ξ + ∫ t s µ(r, ζs,r)dr + ∫ t s σ(r, ζs,r)dWr.
This solution exists and is unique by (Protter, 2005, Chap. V, Thm. 7), therefore we have that Xt = ζ0,t(x) P-almost surely. Furthermore, (Gubinelli, 2016, Thm. 4) implies that for s ≤ t we have Xt = ζs,t(ζ0,s(x)). Hence, for t = s+ r, we have the identity µ̃(s, r, ξ) = E[µ(t, ζs,t(ξ))].
Furthermore, by (Gubinelli, 2016, Thm. 8) we have for any ξ, ξ′ ∈ RdX and (s, r), (s′, r′) ∈ ∆ with t := s+ r, t′ := s′ + r′ ∈ [0, T ] that there exists some constant C such that
E [ |ζs,t(ξ)− ζs′,t′(ξ′)|22 ] ≤ C [ |ξ − ξ′|22 + (1 + |ξ|2 + |ξ′|2) 2 (|t− t′|+ |s− s′|) ] . (16)
Therefore, we have that
|µ̃(s, r, ξ)− µ̃(s′, r′, ξ′)|2 = |E [µ(t, ζs,t(ξ))]− E [µ(t′, ζs′,t′(ξ′))]|2 ≤ |E [µ(t, ζs,t(ξ))]− E [µ(t′, ζs,t(ξ))]|2
+ |E [µ(t′, ζs,t(ξ))]− E [µ(t′, ζs′,t′(ξ′))]|2 ≤ E [ |µ(t, ζs,t(ξ))− µ(t′, ζs,t(ξ))|22 ]1/2 + E [ |µ(t′, ζs,t(ξ))− µ(t′, ζs′,t′(ξ′))|22
]1/2 ≤ E [ |µ(t, ζs,t(ξ))− µ(t′, ζs,t(ξ))|22
]1/2 + M̃ E [ |ζs,t(ξ)− ζs′,t′(ξ′)|22 ]1/2 ,
(17)
where we used Jensen’s inequality in the second last and (10) in the last step. Hence, for (s′, r′, ξ′)→ (s, r, ξ) we have that the first term of (17) goes to zero due to continuity of µ in its first component uniformly in the second component. Moreover, the second term of (17) converges to zero by (16). Together, this proves continuity of µ̃.

heading: C RECALL: RNN, NEURAL ODE AND ODE-RNN - EQUIVALENT WAYS OF WRITING
text: We describe our model as the solution of the SDE (29), which is a compact way of writing. In the following, we first shortly recall recurrent neural networks (RNN) and the neural ODE and then recall the ODE-RNN. Furthermore, we give a step-by-step explanation, how the way ODE-RNN was formalized, can equivalently be written in terms of an SDE similar to (29). Finally we give the alternative way of writing our model and a short comparison to ODE-RNN.
Recurrent Neural Network. The input to a RNN is a discrete time series of observations {xt1 , · · · , xtn}. At each time ti+1, the latent variable h is updated using the previous latent variable hti and the input xti+1 as hti+1 := RNNCell(hti , xti+1), (18) where RNNCell is a neural network.
Neural Ordinary Differential Equation. Neural ODEs (Chen et al., 2018) are a family of continuous-time models defining a latent variable ht := h(t) to be the solution to an ODE initial-value problem (IVP):
ht := ht0 + ∫ t t0 f(hs, s, θ)ds, t ≥ t0, (19)
where f(·, ·, θ) = fθ is a neural network with weights θ. Therefore, the latent variables can be updated continuously by solving this ODE (19). We can emphasize the dependence of ht on a numerical ODE solver by rewriting (19) as
ht := ODESolve(fθ, ht0 , (t0, t)) . (20)
ODE-RNN. ODE-RNN (Rubanova et al., 2019) is a mixture of a RNN and a neural ODE. In contrast to a standard RNN, we are not only interested in an output at the observation times ti, but also in between those times. In particular, we want to have an output stream that is generated continuously in time. This is achieved by using a neural ODE to model the latent dynamics between two observation times, i.e. for ti < t < ti+1 the latent variable is defined as in (19) and (20), with ht0 and t0 replaced by hti and ti. At the next observation time ti+1, the latent variable is then updated by a RNN. Rubanova et al. (2019) write this as{
h′ti+1 := ODESolve(fθ, hti , (ti, ti+1)) hti+1 := RNNCell(h ′ ti+1 , xti+1) .
(21)
Therefore, fixing ht0 := h0, the entire latent process can be computed by iteratively solving an ODE followed by applying a RNN.
GRU-ODE-Bayes. The model architecture describing the latent variable in GRU-ODE-Bayes (Brouwer et al., 2019) is defined as the ODE-RNN but with the special choice of a gated recurrent unit (GRU) for the RNN-cell and the neural network fθ also being derived from a GRU.
ODE-RNN as càdlàg process. Thinking about the process h := (ht)t≥t0 defined in (21) as a (stochastic) process in time, it is defined to evolve continuously for ti ≤ t < ti+1 according to the ODE dynamics fθ and jumps at time ti+1 according to the RNN cell. In particular, it is defined to be a càdlàg5 process, for which hti+1− is the standard notation for the left limit, i.e. the last point before the jump at time ti+1. According to this notation we have hti+1− = h ′ ti+1 , hence, we can rewrite (21)
as { hti+1− := ODESolve(fθ, hti , (ti, ti+1)) hti+1 := RNNCell(hti+1− , xti+1) .
(22)

heading: C.1 RESIDUAL ODE-RNN AS A SPECIAL CASE OF CONTROLLED ODE
text: Residual ODE-RNN. We replace the standard RNN cell by a residual RNN cell (rRNN), as it was described e.g. in Yue et al. (2018). In particular, instead of applying the RNN cell such that hti = RNNCell(hti−, xi) we use a residual RNN cell to have hti = hti− + rRNNCell(hti−, xi). The residual RNN is as expressive as the standard RNN and was empirically shown to perform very similarly or even better than the standard framework (Yue et al., 2018). This way, the residual RNN cell models exactly the jump of the latent variable (i.e. the differences) that occurs at the time ti+1 when taking into account the next observation xti+1 . Therefore, we can rewrite the ODE-RNN (22) as {
hti+1− := ODESolve(fθ, hti , (ti, ti+1)) hti+1 := hti+1− + rRNNCell(hti+1− , xi) .
(23)
Controlled Ordinary Differential Equation. We briefly recall the definition of controlled ODEs as it was given in (Herrera et al., 2020, Section 4.1) and used in (Cuchiero et al., 2019; Herrera et al., 2020) to describe neural networks. We fix `, d ∈ N and define for 1 ≤ i ≤ d the vector fields
Vi : Θ× R≥0 × R` → R`, (θ, t, x) 7→ V θi (t, x),
which are càglàd6 in t and Lipschitz continuous in x. Furthermore, we define the scalar càdlàg control functions ui : R≥0 → R, t 7→ ui(t), which have finite variation and satisfy ui(0) = 0. Then we define the process Z := (Zt)t≥0 as the solution of the controlled ODE
dZt = d∑ i=1 V θi ( t, Zt− ) dui(t), Z0 = z, (24)
where z ∈ R` is some starting point. (24) is written in Itô’s differential notation for (stochastic) integrals. The solution of (24) exists and is unique under much more general assumptions than we made here (Protter, 2005, Chap. V, Thm. 7).
5i.e. right-continuous with existing left limits, also denoted as RCLL 6i.e. left continuous with existing right limits
Special Case: ODE-RNN. To write the ODE-RNN as controlled ODE, we set d = 2 and define V θ1 := fθ, V2 := rRNNCell and u1(t) := t. As pointed out by Herrera et al. (2020), one can take the ui to be semimartingales instead of deterministic functions. In line with this, we define u2 := u as the pure jump stochastic control process
u : Ω̃× [0, T ]→ R, (ω̃, t) 7→ ut(ω̃) := n(ω̃)∑ i=1 1[ti(ω̃),∞)(t). (25)
We note that u is an adapted process starting at 0 with finite variation on the product probability space (Ω× Ω̃,F ⊗ F̃ ,F⊗ F̃ ,P× P̃), since the total variation of u up to time T is n and EP×P̃[n] < ∞. The following result shows that the residual ODE-RNN can compactly be described as a controlled ODE. Proposition C.1. Using the vector fields and controls defined above, the latent variable process h = (ht)t≥t0 of the residual ODE-RNN can equivalent be written as the solution of the controlled ODE
dht = fθ(ht−, t−)dt+ rRNNCell(ht−, xt)du(t), ht0 = h0. (26)
Proof of Proposition C.1. By (Protter, 2005, Chap. II, Thm. 17), the stochastic integral is indistinguishable from the path-wise Lebesgue-Stieltjes integral if the integrator is of finite variation. Hence, we can assume that some ω̃ ∈ Ω̃ is fixed and that the following expressions are evaluated at this ω̃ whenever applicable. First note, that u is constant except at the ti where it increases by 1 (cf. Figure 7). In particular, the Lebesgue-Stieltjes integral of some càdlàg function g with respect to u is a sum, i.e. ∫ t
0 g(s−)dus = ∑ ti≤t g(ti−)∆uti , where ∆uti = uti − uti− = 1. Therefore, integrating (26)
from t0 to t we get
ht = ht0 + ∫ t t0 fθ(hs−, s−)ds+ ∫ t t0 rRNNCell(hs−, xs)du(s)
= ht0 + ∫ t t0 fθ(hs−, s−)ds+ ∑ ti≤t rRNNCell(hti−, xti).
In particular, since the first integral is continuous in t, we have for every 1 ≤ k ≤ n
htk = ht0 + ∫ tk− t0 fθ(hs−, s−)ds+ ∑ ti<tk rRNNCell(hti−, xti) + rRNNCell(htk−, xtk)
= htk− + rRNNCell(htk−, xtk),
(27)
and for tk < t < tk+1
ht = ht0 + ∫ tk t0 fθ(hs−, s−)ds+ ∑ ti≤tk rRNNCell(hti−, xti) + ∫ t tk fθ(hs−, s−)ds
= htk + ∫ t tk fθ(hs−, s−)ds. (28)
Together, (27) and (28) prove that (26) is equivalent to (23). We also emphasize that xt is used as input only for t = ti, 1 ≤ i ≤ n.

heading: D NEURAL JUMP ODE
text: We propose a model framework whose architecture can be interpreted as a simplification of the ODERNN (Rubanova et al., 2019) and the GRU-ODE-Bayes (Brouwer et al., 2019) model architecture.

heading: D.1 THE MODEL FRAMEWORK: NEURAL ODE BETWEEN JUMPS
text: We define X ⊂ RdX andH ⊂ RdH to be the observation and latent space for dX , dH ∈ N. Moreover, we define the following feed-forward neural networks with sigmoid activation functions:
• fθ1 : RdH × RdX × [0, T ]× [0, T ]→ RdH modelling the ODE dynamics, • ρθ2 : RdX → RdH modelling the jumps when new observations are made, and • gθ3 : RdH → RdY the readout map, mapping into the target space Y ⊂ RdY for dY ∈ N.
The trainable parameters of the neural networks are θ := (θ1, θ2, θ3) ∈ Θ := ∪M≥1ΘM . Here, for every M ∈ N, ΘM is defined to be the set of all parameters such that fθ1 and ρθ2 have hidden dimension M . We define the pure jump stochastic process (cf. Appendix C)
u : Ω̃× [0, T ]→ R, (ω̃, t) 7→ ut(ω̃) := n(ω̃)∑ i=1 1[ti(ω̃),∞)(t).
Then the Neural Jump ODE (NJ-ODE) is defined by the latent process H := (Ht)t∈[0,T ] and the output process Y := (Yt)t∈[0,T ] defined as solutions of the SDE system (cf. Appendix C)
H0 = ρθ2 (X0) , dHt = fθ1 ( Ht−, Xτ(t), τ(t), t− τ(t) ) dt+ (ρθ2 (Xt)−Ht−) dut,
Yt = gθ3(Ht).
(29)
Note that, only the values of X at the times ti, 0 ≤ i ≤ n are used as inputs. We will use the notation Hθt (X) and Y θ t (X) to emphasize the dependence of the latent process H and the output process Y on the model parameters θ and the input X .
Existence and uniqueness. We note that u is an adapted process starting at 0 with finite variation on the product probability space (Ω× Ω̃,F ⊗ F̃ ,F⊗ F̃ ,P× P̃), since the total variation of u up to time T is n and EP×P̃[n] <∞. Hence, since fθ1 and ρθ2 are Lipschitz continuous (as composition of Lipschitz continuous functions) a unique càdlàg solution Hθ exists, once an initial value is fixed (Protter, 2005, Thm. 7, Chap. V). Moreover, the resulting process (Y θt )t∈[0,T ] is also càdlàg and A-adapted.
Equivalent way of writing NJ-ODE. By applying the steps that were explained in Section C backwards to (29) (but without introducing the residual version of the RNN), it is easy to see that our model can equivalently be written (similar to (22)) as{
hti+1− := ODESolve(fθ1 , (ht, xti , ti, t− ti), (ti, ti+1)) hti+1 := ρθ2(xti+1) .
(30)
In particular, the main difference to (22) is that we use a modified input for fθ1 and that the neural network performing the jumps does not take hti+1− as an input, i.e. it is not a RNN.

heading: D.2 OBJECTIVE FUNCTION
text: We introduce a loss function, such that the output Y θ of the NJ-ODE can be trained to approximate X̂ , i.e. to model the best online prediction of the stochastic process X . We define the theoretical loss function and its Monte Carlo and Ergodic approximation.
Objective function. Let us define D to be the set of all RdX -valued A-adapted processes on the probability space (Ω× Ω̃,F ⊗ F̃ ,F⊗ F̃ ,P× P̃) . Then we define our objective functions
Ψ :D→ R, Z 7→ Ψ(Z) := EP×P̃
[ 1
n n∑ i=1 (|Xti − Zti |2 + |Zti − Zti−|2) 2
] , (31)
Φ :Θ→ R, θ 7→ Φ(θ) := Ψ(Y θ(X)), (32)
where Φ will be our (theoretical) loss function. Remark that from the definition of Y θ, it directly follows that it is an element of D, hence Φ is well-defined.
Monte Carlo approximation of the objective function. Let us assume, that we observe N ∈ N independent realisations of the path X at times (t̃(j)1 , · · · , t̃ (j) nj ), 1 ≤ j ≤ N , which are themselves independent realisations of the random vector (n, t1, · · · , tn). In particular, let us assume that X(j) ∼ X and (nj , t(j)1 , · · · , t (j) nj ) ∼ Pt are i.i.d. random processes (respectively variables) for 1 ≤ j ≤ N and that our training data is one realisation of them. Then
Φ̂N (θ) := 1
N N∑ j=1 1 nj nj∑ i=1 (∣∣∣∣X(j)t(j)i − Y θt(j)i (X(j)) ∣∣∣∣ 2 + ∣∣∣Y θ t (j) i (X(j))− Y θ t (j) i − (X(j)) ∣∣∣ 2 )2 , (33)
converges (P× P̃)-a.s. to Φ(θ) as N →∞, by the law of large numbers (cf. Theorem E.13). Ergodic approximation of the objective function. If we only observe one realization of the path X at times (t̃1, · · · , t̃N ), we can still approximate the objective function by assuming that µ and σ are time-independent and that the stochastic process X is ergodic in the following sense. We fix n = 1 and assume that the time increments ∆t̃j := t̃j − t̃j−1 are i.i.d. realizations of the probability distribution λ1. Furthermore, we consider each observation Xt̃j as one sample with initial condition Xt̃j−1 for which Y θ,j is the realization of Y θ. Then we approximate the objective function by
Φ̂N (θ) := 1
N N∑ j=1 (∣∣∣Xt̃j − Y θ∆t̃j (X)∣∣∣2 + ∣∣∣Y θ∆t̃j (X)− Y θ∆t̃j−(X)∣∣∣2)2 , (34) which is assumed to converge by the ergodicity assumption for N →∞ to
EP×P̃ [(∣∣Xt1 − Y θt1 ∣∣2 + ∣∣Y θt1 − Y θt1−∣∣2)2] . (35)
Instead of setting the random variable n = 1, one could similarly fix the time horizon T , and take for each sample all subsequent observations that lie in the time interval [tstart, tstart + T ], where tstart is the date of the first observation of this sample. The next sample would then start with the first observation after tstart + T .

heading: E THEORETICAL CONVERGENCE RESULTS
text: Our main results show, that the output of NJ-ODE converges to the conditional expectation when the size of the neural networks and the number of samples go to infinity. For each network size (and number of samples) we assume to have the weights minimizing the (Monte Carlo approximation of the) loss function. In particular, we do not consider the problem of finding those optimal weights, and therefore also do not analyse backpropagation through our model. In a similar setting, backpropagation was studied by Jia & Benson (2019).
For completeness we recall the definition of Lp-convergence. Definition E.1. Let 1 ≤ p <∞. Let (Ω,F ,P) be a probability space. Then a sequence of random variables (Xn)n∈N converges to a random variable X in Lp(Ω,F ,P) (or simply Lp), if
E[|Xn −X|p] n→∞−−−−→ 0.
We use the notation Xn
Lp−−→ X. L1-convergence is also denoted convergence in mean.

heading: E.1 CONVERGENCE WITH RESPECT TO THEORETICAL OBJECTIVE FUNCTION
text: Theorem E.2. Let θminM ∈ ΘminM := argminθ∈ΘM {Φ(θ)} for every M ∈ N. Then, for M →∞, the value of the loss function Φ (32) converges to the minimal value of Ψ (31) which is uniquely achieved by X̂ , i.e.
Φ(θminM ) M→∞−−−−→ min
Z∈D Ψ(Z) = Ψ(X̂).
Furthermore, for every 1 ≤ k ≤ K we have that Y θminM converges to X̂ as random variable in L1(Ω× [0, T ],P× λk). In particular, the limit process Y := limM→∞ Y θ min M equals X̂ (P× λk)almost surely as a random variable on Ω× [0, T ].
The idea of the proof is to split up the target jump process into its continuous-in-time parts and into the jumps. Both parts are continuous functions of their inputs and can therefore be approximated by neural networks (Hornik et al., 1989). More precisely, by Proposition B.1 and B.4 the continuous-intime part can be written as an integral over time of a function which is jointly continuous in its inputs. This jointly continuous function is approximated by a neural network, which is itself integrated over time – that is the neural ODE part of the NJ-ODE. The remainder of the proof shows L1-convergence of the output of NJ-ODE with optimal weights (with respect to the loss function) to the conditional expectation process.
For completeness we restate the straight forward generalization of the universal approximation result (Hornik et al., 1989, Thm. 2.4) to multidimensional output.
Theorem E.3 (Hornik). Let r, d ∈ N and σ be a sigmoid function. Let NN σr,d be the set of all 1-hidden layer neural networks mapping Rr to Rd. Then for every compact subset K ⊂ Rr, every > 0 and every f ∈ C(Rr,Rd) there exists a neural network g ∈ NN σr,d such that supx∈K |f(x)− g(x)|2 < .
The following Lemmas are used in the Proof of Theorem E.2.
Lemma E.4. Let 1 ≤ k ≤ K and let Z ∈ D be a process such that (P × λk)[X̂ 6= Z] > 0. Then there exists an ε > 0 such that B̃ := {t ∈ [0, T ] |EP[|Xt−Zt−|22] ≥ ε+EP[|Xt− X̂t−|22]} satisfies λk(B̃) > 0.
Proof. First remark that since X is continuous, we have Xt = Xt−. Let us define C := {(ω, t) ∈ Ω × [0, T ] | X̂t−(ω) 6= Zt−(ω)} and for each t ∈ [0, T ] let Ct := {ω ∈ Ω | (ω, t) ∈ C}. Then we have for B := {t ∈ [0, T ] |P(Ct) > 0} that λk(B) > 0, since otherwise by Fubini’s theorem
0 < (P× λk)[C] = ∫
[0,T ]
P(Ct)dλk(t) = 0,
which is a contradiction. Now Proposition B.2 yields that for each t ∈ B there exists some εt > 0 such that EP[|Xt− − Zt−|22] ≥ εt + EP[|Xt− − X̂t−|22]. This implies the claim. Indeed, assume no such ε > 0 exists, then we have for each n ∈ N that λk ( {t ∈ [0, T ] |EP[|Xt− − Zt−|22] ≥ 1n + EP[|Xt− − X̂t−| 2 2]} ) = 0. Therefore,
λk ( {t ∈ [0, T ] |EP[|Xt− − Zt−|22] > EP[|Xt− − X̂t−|22]} ) ≤ ∑ n∈N λk ( {t ∈ [0, T ] |EP[|Xt− − Zt−|22] ≥ 1n + EP[|Xt− − X̂t−| 2 2]} ) = 0,
which is a contradiction to λk(B) > 0.
Lemma E.5. For any A-adapted process Z it holds that
EP×P̃ [ 1 n n∑ i=1 |Xti − Zti−| 2 2 ] = EP×P̃ [ 1 n n∑ i=1 ∣∣∣Xti − X̂ti−∣∣∣2 2 ] + EP×P̃ [ 1 n n∑ i=1 ∣∣∣X̂ti− − Zti−∣∣∣2 2 ] .
Proof. First recall that by continuityXti = Xti−. Then the statement is a consequence of Proposition B.2, Lemma B.3 and Fubini’s theorem, which imply
EP×P̃ [ 1 n n∑ i=1 |Xti− − Zti−| 2 2 ] = EP̃ [ 1 n n∑ i=1 EP [ |Xti− − Zti−| 2 2 ]]
= EP̃ [ 1 n n∑ i=1 ( EP [∣∣∣Xti− − X̂ti−∣∣∣2 2 ] + EP [∣∣∣X̂ti− − Zti−∣∣∣2 2 ])]
= EP×P̃ [ 1 n n∑ i=1 ∣∣∣Xti− − X̂ti−∣∣∣2 2 ] + EP×P̃ [ 1 n n∑ i=1 ∣∣∣X̂ti− − Zti−∣∣∣2 2 ] .
Lemma E.6. Let 1 ≤ p < ∞. Let (Zn)n∈N be a sequence of random variables, and Z and Z̃ random variables defined on a common probability space such that Zn Lp−−→ Z and Zn Lp−−→ Z̃. Then

heading: Z = Z̃ almost surely.
text: Proof. With the triangle inequality it follows that E[|Z − Z̃|p]1/p = 0, which implies the claim.
Lemma E.7. The random variable ST := sup0≤t≤T |Xt|1 is square integrable and bounded in probability, i.e. for any ε > 0 exist some K > 0 such that P[ST > K] ≤ ε.
Proof. From the proof of Proposition B.1 we know that
Mt := ∫ t 0 σ(s,Xs)dWs, 0 ≤ t ≤ T,
is a square integrable martingale and that EP[sup0≤t≤T |Mt| j 1] ≤ cEP[|MT | j 1] < ∞, for j = 1, 2 and some constant c > 0 by Doob’s inequality. Moreover, µ is bounded, say by B, hence
|Xt|1 = ∣∣∣∣x+ ∫ t
0
µ(r,Xr)dr + ∫ t 0 σ(r,Xr)dWr ∣∣∣∣ 1
≤ |x|1 + ∫ t
0
|µ(r,Xr)|1 dr + |Mt|1
≤ |x|1 +B T + |Mt|1 .
Therefore, EP[ST ] ≤ |x|1 +B T + cEP[|MT |1] <∞ and similar EP[S2T ] <∞, which implies the claim. Indeed, if for a fixed ε > 0 no such K exists, then P[ST =∞] ≥ ε, which is a contradiction to integrability.
Proof of Theorem E.2. We start by showing that X̂ ∈ D is the unique minimizer of Ψ up to (P×λk)null-sets for any k ≤ K. First, we recall that for every ti we have X̂ti = Xti and by continuity of X that Xti = Xti−. Therefore,
Ψ(X̂) = EP×P̃
[ 1
n n∑ i=1 ∣∣∣Xti − X̂ti−∣∣∣2 2
]
= EP̃
[ 1
n n∑ i=1 EP [∣∣∣Xti − X̂ti−∣∣∣2 2 ]]
= min Z∈D EP×P̃
[ 1
n n∑ i=1 |Xti − Zti−| 2 2
]
≤ min Z∈D EP×P̃
[ 1
n n∑ i=1 (|Xti − Zti |2 + |Zti − Zti−|2) 2 ] = min
Z∈D Ψ(Z),
where we used Fubini’s theorem for the second line, Proposition B.2 for the third line and the triangle inequality for the fourth line. Hence, X̂ is a minimizer of Ψ. To see that it is unique (P× λk)-a.s., let Z ∈ D be a process such that (P × λk)[X̂ 6= Z] > 0. By Lemma E.4, this implies that there exists an ε > 0 such that B := {t ∈ [0, T ] |EP[|Xt− − Zt−|22] ≥ + EP[|Xt− − X̂t−|22]} satisfies λk(B) > 0.
Now recall that by definition of λk we have λk(B) = P̃(∪̇j≥k{n = j, tk− ∈ B})/P̃(n ≥ k) > 0. This implies that there exists j ∈ N≥k such that P̃(n = j, tk− ∈ B) > 0. Therefore,
EP̃
[ 1
n n∑ i=1 1{ti−∈B} ] ≥ EP̃ [ 1{n=j} 1 n n∑ i=1 1{ti−∈B} ] ≥ EP̃ [ 1{n=j} 1 j 1{tk−∈B}
] = 1j P̃(n = j, tk− ∈ B) > 0.
This inequality implies now that Z is not a minimizer of Ψ, because
Ψ(Z) = EP×P̃
[ 1
n n∑ i=1 (|Xti − Zti |2 + |Zti − Zti−|2) 2
]
≥ EP̃
[ 1
n n∑ i=1 EP [ |Xti − Zti−| 2 2
]]
= EP̃
[ 1
n n∑ i=1 ( 1{ti−∈B} + 1{ti−∈BC} ) EP [ |Xti − Zti−| 2 2
]]
≥ EP̃
[ 1
n n∑ i=1 ( ε1{ti−∈B} + EP [∣∣∣Xti − X̂ti−∣∣∣2 2 ])]
= εEP̃
[ 1
n n∑ i=1 1{ti−∈B} ] + min Z∈D Ψ(Z)
> min Z∈D Ψ(Z).
Next we show that (29) can approximate X̂ arbitrarily well. Since the dimension dH can be chosen freely, let us fix it to dH := dX . Furthermore, let us fix θ∗2 and θ ∗ 3 such that ρθ∗2 = gθ∗3 = id. From Theorem E.3 it follows that for any ε > 0 there exist M ∈ N and θ∗1 with (θ∗1 , θ∗2 , θ∗3) ∈ ΘM such that
sup (u,v,t,r)∈[−M,M ]dH×dX×∆ ∣∣fθ∗1 (u, v, t, r)− µ̃ (t, r, v)∣∣2 ≤ ε, (36) where we used that µ̃ is continuous by Proposition B.4. Since µ is bounded, also µ̃ is bounded, say by B − 1 > 0. On [−M,M ]dX we approximate µ̃ by the neural network fθ∗1 and outside we continuously extend fθ∗1 such that it is bounded by B. By abuse of notation we call this fθ∗1 again and use it as our neural network. Hence,
∣∣fθ∗1 − µ̃∣∣2 ≤ 2B. Using this, we can bound the distance between Y θ ∗ M
t and X̂ . In particular, if t ∈ {t1, · · · , tn}, we have∣∣∣Y θ∗Mt − X̂t∣∣∣ 2 = ∣∣(Ht− + (ρθ∗M (Xt)−Ht−))−Xt∣∣2 = 0,
and if t not in {t1, · · · , tn}, then (12), (36) and the previous bound yield for ST := sup0≤t≤T |Xt|1∣∣∣Y θ∗Mt − X̂t∣∣∣ 1 ≤ ∣∣∣Y θ∗Mτ(t) − X̂τ(t)∣∣∣ 1
+ ∫ t τ(t) ∣∣fθ∗1 (Hs−, Xτ(s), τ(s), s− τ(s))− µ̃ (τ(s), s− τ(s), Xτ(s))∣∣1 ds ≤ ε T1{ST≤M} + 2BT1{ST>M} ≤ ε T + 2BT1{ST>M}.
Here we used that ST ≤ M implies that (36) can be used for all τ(t) ≤ s ≤ t. Moreover, by equivalence of the 1- and 2-norm, there exists a constant c > 0 such that∣∣∣Y θ∗Mt − X̂t∣∣∣
2 ≤ c ε T + 2cBT1{ST>M}.
With Lemma E.7 we know that EP[1{ST>M}] = P[ST > M ] =: M M→∞−−−−→ 0. Now we can show the convergence of Φ(θminM ) using these two bounds and Xti = X̂ti . Indeed,
min Z∈D
Ψ(Z) ≤ Φ(θminM ) ≤ Φ(θ∗M )
= EP×P̃
[ 1
n n∑ i=1 (∣∣∣Xti − Y θ∗Mti ∣∣∣ 2 + ∣∣∣Y θ∗Mti − Y θ∗Mti− ∣∣∣ 2
)2]
≤ EP×P̃
[ 1
n n∑ i=1 (∣∣∣X̂ti − Y θ∗Mti ∣∣∣ 2 + ∣∣∣Y θ∗Mti − X̂ti∣∣∣ 2 + ∣∣∣X̂ti − X̂ti−∣∣∣ 2 + ∣∣∣X̂ti− − Y θ∗Mti− ∣∣∣ 2
)2]
≤ EP×P̃
[ 1
n n∑ i=1 (∣∣∣Xti − X̂ti−∣∣∣ 2 + c(ε T + 2BT1{ST>M}) )2]
≤ EP̃
[ 1
n n∑ i=1 EP [(∣∣∣Xti − X̂ti−∣∣∣ 2 + c(ε T + 2BT1{ST>M}) )2]]
≤ EP̃  1 n n∑ i=1 ( EP [∣∣∣Xti − X̂ti−∣∣∣2 2 ]1/2 + EP [ (c ε T + 2cBT1{ST>M}) 2 ]1/2)2 ,
where we used the triangle-inequality for the L2-norm in the last step. We can bound EP [ (c ε T + 2cBT1{ST>M}) 2 ] ≤ 3(c ε T )2 + 3(2cBT )2 M =: c2M ,
which is a constant converge to 0 as M →∞. Using that for a ∈ R, we have a ≤ a2 + 1, we get
min Z∈D
Ψ(Z) ≤ Φ(θminM ) ≤ Φ(θ∗M )
≤ EP̃  1 n n∑ i=1 ( EP [∣∣∣Xti − X̂ti−∣∣∣2 2 ]1/2 + c(M) )2 ≤ EP̃ [ 1
n n∑ i=1 ( EP [∣∣∣Xti − X̂ti−∣∣∣2 2 ] + 2cMEP [∣∣∣Xti − X̂ti−∣∣∣2 2 ]1/2 + c2M )]
≤ (1 + 2cM )EP×P̃
[ 1
n n∑ i=1 ∣∣∣Xti − X̂ti−∣∣∣2 2 ] + 2cM + c 2 M
≤ (1 + 2cM ) min Z∈D Ψ(Z) + 2cM + c 2 M M→∞−−−−→ min Z∈D Ψ(Z),
where we used Ψ(X̂) = minZ∈D Ψ(Z) and that cM converges to 0.
In the last step we show that the limits limM→∞ Y θ min M and limM→∞ Y θ ∗ M exist as limits in the Banach space L := L1(Ω× [0, T ],F⊗ B([0, T ]),P× λk), for every k ≤ K, and that they are both equal to X̂ . Let us fix k ≤ K. First we note that for every B ∈ B([0, T ]) we have
Eλk [1B ] = λk(B) = P̃(n≥k,tk−∈B) P̃(n≥k) = EP̃[1{n≥k}1{tk−∈B}]
P̃(n ≥ k) .
Using “measure theoretic induction” (Durrett, 2010, Case 1-4 of Proof of Thm. 1.6.9) this yields for c := (P(n ≥ k))−1 and a B([0, T ])-measurable function Z : [0, T ]→ R, t 7→ Zt := Z(t) that
Eλk [Z] = cEP̃[1{n≥k}Ztk−]. (37)
Moreover, the triangle inequality and Lemma E.5 yield
Φ(θ∗M )−Ψ(X̂) ≥ EP×P̃
[ 1
n n∑ i=1 ∣∣∣Xti − Y θ∗Mti− ∣∣∣2 2 ] −Ψ(X̂)
= EP×P̃
[ 1
n n∑ i=1 ∣∣∣X̂ti− − Y θ∗Mti− ∣∣∣2 2 ] .
(38)
For any RdX -valued Z ∈ L the Hölder inequality, together with the fact that n ≥ 1, yields
EP×P̃ [|Z|2] = EP×P̃ [√ n√ n |Z|2 ] ≤ EP×P̃ [n] 1/2 EP×P̃ [ 1 n |Z|22 ]1/2 . (39)
Together, this implies that limM→∞ Y θ ∗ M = X̂ as a L-limit. Indeed, with c̃ := EP×P̃ [n] 1/2 <∞ we have
EP×λk [∣∣∣X̂ − Y θ∗M ∣∣∣
2
] = cEP×P̃ [ 1{n≥k} ∣∣∣X̂tk− − Y θ∗Mtk−∣∣∣2] ≤ c c̃EP×P̃ [ 1{n≥k} 1
n ∣∣∣X̂tk− − Y θ∗Mtk−∣∣∣22 ]1/2
≤ c c̃EP×P̃ [ 1{n≥k} 1
n n∑ i=1 ∣∣∣X̂ti− − Y θ∗Mti− ∣∣∣2 2
]1/2
≤ c c̃EP×P̃
[ 1
n n∑ i=1 ∣∣∣X̂ti− − Y θ∗Mti− ∣∣∣2 2 ]1/2 ≤ c c̃ ( Φ(θ∗M )−Ψ(X̂)
)1/2 M→∞−−−−→ 0, where we used first (37) and (39) followed by two simply upper bounds and (38) in the last step. The same argument can be applied to show that limM→∞ Y θ min M = X̂ as a L-limit. In particular, this proves that the limit Y := limM→∞ Y θ min M exists as L-limit and by Lemma E.6 it equals X̂ (P× λk)-almost surely, for any k ≤ K.
Remark E.8. This result can be extended to any other neural network architecture for which a universal approximation theorem equivalent to (Hornik et al., 1989, Thm. 2.4) exists. Moreover, the stochastic process X defined in (9) can be chosen more general, in particular, the diffusion part∫ σdW can be replaced by any martingale, as long as the resulting process still is a Markov process and µ̃ stays continuous.
Remark E.9. If we used the modified loss function Ψ̃ which is identical to Ψ except that we drop the factor 1n , everything would work similarly and we could show L
2-convergence instead of L1convergence of Y θ min M to X̂ . However, we remark that there might exists Z ∈ D, such that Ψ(Z) <∞
while Ψ̃(Z) = ∞. In particular, if some moment of n does not exist, such a process can be constructed. Remark E.10. It follows directly from Theorem E.2 that limM→∞ Y θ min M = X̂ as random variables on Ω× [0, T ] except on sets which are null sets with respect to every product measure P× λk for 1 ≤ k ≤ K. Remark E.11. The result of Theorem E.2 does not imply that X̂ and limM→∞ Y θ min M are modifications or indistinguishable. For example, if B ⊂ [0, T ] is a subset such that no left-point of the observation times (tk−) lies in B with probability greater 0, i.e. λk(B) = 0 for 1 ≤ k ≤ K, then Theorem E.2 does not tell us how close (limM→∞ Y θ min M )t is to X̂t for t ∈ B. In particular, it does not tell us whether they are equal P-almost surely. Furthermore, such a set B always exists, since there has to exists t ∈ [0, T ] such that B := {t} has measure 0 for all k.
In the following corollary we show, that Theorem E.2 can be extended to show convergence to the conditional expectation of ϕ(X), for some function ϕ ∈ C2,b(RdX ,R), i.e. a function that is twice continuously differentiable with bounded derivatives.
Corollary E.12. Let ϕ ∈ C2,b(RdX ,R), then the statement of Theorem E.2 holds equivalently, when replacing X in the loss functions Ψ and Φ by Γ = ϕ(X) and X̂ by the conditional expectation Γ̂, where Γ̂t := EP×P̃[ϕ(Xt)|At].
Corollary E.12 combined with the monotone convergence theorem for conditional expectation theoretically enables us to make statements about the conditional law and conditional moments of X under some a priori integrability assumptions.
Proof of Corollary E.12. We first remark that Proposition B.2 and Lemma E.4, E.5, B.3 hold similarly for the conditional expectation Γ̂. Hence, the same argument as in the proof of Theorem E.2 implies that Γ̂ is the unique A-adapted minimizer of Ψ. For simplicity of the notation we assume that dX = dY = 1, i.e. that the process X and the Brownian motion W are 1-dimensional. However, the following works as well in the general case, where the correlations of the Brownian motion components have to be taken into account. By Itô’s Formula (Protter, 2005, Chap. II, Thm. 32), Γ = ϕ(X) is the solution of the SDE
dϕ(X)t = ϕ ′(Xt)dXt + 1 2ϕ ′′(Xt)d[X,X]t
= ϕ′(Xt)µ(t,Xt)dt+ ϕ ′(Xt)σ(t,Xt)dWt + 1 2ϕ ′′(Xt)σ(t,Xt) 2dt
= α(t,Xt)dt+ β(t,Xt)dWt,
for α(t,Xt) := ϕ′(Xt)µ(t,Xt) + 12ϕ ′′(Xt)σ(t,Xt) 2 and β(t,Xt) := ϕ′(Xt)σ(t,Xt). Defining α̃ similar to µ̃ as
α̃ : ∆× RdX → RdX , (t, r, ξ) 7→ Pt,t+r(Xt, α) ∣∣ Xt=ξ = EP [α(t+ r,Xt+r)|Xt = ξ] ,
one can use the boundedness of of ϕ′ and ϕ′′ to show that it is continuous and that Γ̂t = E[ϕ(X)t|At] = Γτ(t) + ∫ t τ(t) α̃ ( τ(t), s− τ(t), Xτ(t) ) ds.
In particular, Proposition B.1 and B.4 hold equivalently for Γ. Similar to (36), the neural network parameters can be chosen such that
sup (u,v,t,r)∈[−M,M ]dH×dX×∆ ∣∣fθ∗1 (u, v, t, r)− α̃ (t, r, v)∣∣2 ≤ ε, (40) which then implies the statement of the Corollary similar as in the proof of Theorem E.2.

heading: E.2 CONVERGENCE OF THE MONTE CARLO APPROXIMATION
text: In the following, we assume the size of the neural network M is fixed and we study the convergence with respect to the number of samples N . Moreover, we show that both types of convergence can be combined. To do so, we define Θ̃M := {θ ∈ ΘM | |θ|2 ≤M}, which is a compact subspace of ΘM . It is straight forward to see, that ΘM in Theorem E.2 can be replaced by Θ̃M . Indeed, if the needed neural network weights for an ε-approximation have too large norm, then one can increase M until it is sufficiently big. The following convergence analysis is based on (Lapeyre & Lelong, 2019, Chapter 4.3).
Theorem E.13. Let θminM,N ∈ ΘminM,N := arg infθ∈Θ̃M {Φ̂N (θ)} for every M,N ∈ N. Then, for every M ∈ N, (P× P̃)-a.s.
Φ̂N N→∞−−−−→ Φ uniformly on Θ̃M .

heading: Moreover, for every M ∈ N, (P× P̃)-a.s.
text: Φ(θminM,N ) N→∞−−−−→ Φ(θminM ) and Φ̂N (θminM,N ) N→∞−−−−→ Φ(θminM ).
In particular, one can define an increasing sequence (NM )M∈N in N such that for every 1 ≤ k ≤ K we have that Y θ min M,NM converges to X̂ for M →∞ as random variable in L1(Ω× [0, T ],P×λk). In particular, the limit process Y := limM→∞ Y θminM,NM equals X̂ (P× λk)-almost surely as a random variable on Ω× [0, T ].
The following Monte Carlo convergence analysis is based on (Lapeyre & Lelong, 2019, Section 4.3). In comparison to them, we do not need the additional assumptions that were essential in (Lapeyre & Lelong, 2019, Section 4.3), i.e. that all minimizing neural network weights generate the same neural network output. This assumption is not needed, because we do not aim to show that Y θ min M,N converges to Y θ min M .
We define the separable Banach space S := {x = (xi)∈N ∈ `1(RdX ) | ‖x‖`1 <∞} with the norm ‖x‖`1 := ∑ i∈N|xi|2.

heading: E.2.1 CONVERGENCE OF OPTIMIZATION PROBLEMS
text: Consider a sequence of real valued functions (fn)n defined on a compact set K ⊂ Rd. Define, vn = infx∈K fn(x) and let xn be a sequence of minimizers fn(xn) = infx∈K fn(x).
From (Rubinstein & Shapiro, 1993, Theorem A1 and discussion thereafter) we have the following Lemma.
Lemma E.14. Assume that the sequence (fn)n converges uniformly on K to a continuous function f . Let v∗ = infx∈K f(x) and S∗ = {x ∈ K : f(x) = v∗}. Then vn → v∗ and d(xn,S∗)→ 0 a.s.
The following lemma is a consequence of (Ledoux & Talagrand, 1991, Corollary 7.10) and (Rubinstein & Shapiro, 1993, Lemma A1).
Lemma E.15. Let (ξi)i≥1 be a sequence of i.i.d random variables with values in S and h : Rd×S → R be a measurable function. Assume that a.s., the function θ ∈ Rd 7→ h(θ, ξ1) is continuous and for all C > 0, E(sup|θ|2≤C |h(θ, ξ1)|) < +∞. Then, a.s. θ ∈ R d 7→ 1N ∑N i=1 h(θ, ξi) converges locally uniformly to the continuous function θ ∈ Rd 7→ E(h(θ, ξ1)),
lim N→∞ sup |θ|2≤C ∣∣∣∣∣ 1N N∑ i=1 h(θ, ξi)− E(h(θ, ξ1)) ∣∣∣∣∣ = 0 a.s.

heading: E.2.2 STRONG LAW OF LARGE NUMBERS
text: Let us define F (x, y, z) := |x− y|2 + |y − z|2
and ξj := (X j tj1 , . . . , Xj tj nj , 0, . . . ), where Xj tji are random variables describing the realizations of the training data, as defined in Section D.2. By this definition we have nj := nj(ξj) := maxi∈N{ξj,i 6= 0} P-almost-surely and we know that ξj are i.i.d. random variables taking values in S . Furthermore, let us write Y θt (ξ) to make the dependence of Y on the input and the weight θ explicit. Then we define
h(θ, ξj) := 1
nj nj∑ i=1 F ( Xtji , Y θ tji (ξj), Y θ tji− (ξj) )2 .
Lemma E.16. The following properties are satisfied.
(P1) There exists κ > 0 such that for all S ∈ S and θ ∈ Θ̃M we have |Y θt (S)|2 ≤ κ (1+|Xτ(t)|2) for all t ∈ [0, T ].
(P2) Almost-surely the random function θ ∈ Θ̃M 7→ Y θt is uniformly continuous for every t ∈ [0, T ].
(P3) We have EP×P̃ [ 1 n ∑n i=1 |Xti |22 ] <∞ and EP×P̃ [ 1 n ∑n i=1 |Xti−|22 ] <∞.
Proof. By definition of the neural networks with sigmoid activation functions (in particular having bounded outputs), all neural network outputs are bounded in terms of the norm of the network weights, which is assumed to be bounded, not depending on the norm of the input. Since after a jump at τ(t), Y has the value Xτ(t), we can find κ depending on T , such that the claimed bound is satisfied for all t, proving (P1).
Since the activation functions are continuous, also the neural networks are continuous with respect to their weights θ, which implies that also θ ∈ Θ̃M 7→ Y θt is continuous. Since Θ̃M is compact, this automatically yields uniform continuity and therefore finishes the proof of (P2). (P3) follows directly from the stronger result in Lemma E.7.
Proof of Theorem E.13. We apply Lemma E.15 to the sequence of i.i.d random function h(θ, ξj). From (P1) of Lemma E.16 we have that
F (Xtji , Ytji , Ytji− )2 = (∣∣∣Xtji − Y θtji ∣∣∣2 + ∣∣∣Y θtji − Y θtji−∣∣∣2)2 ≤ 4 ( |Xtji | 2 2 + |Y θtji | 2 2 + |Y θtji | 2 2 + |Y θtji−| 2 2
) ≤ 4 ( 3|Xtji | 2 2 + κ(1 + |Xtji−1 | 2 2) ) .
Hence, we obtain that
h(θ, ξj) = 1
nj nj∑ i=1 F (Xtji , Ytji , Ytji− )2 ≤ 12 + 4κ nj nj∑ i=1 ∣∣∣Xtji ∣∣∣22 + 4κ+ |x|22, implying that
EP×P̃ [ sup θ∈Θ̃M h(θ, ξj) ] ≤ (12 + 4κ)EP×P̃ [ 1 n n∑ i=1 |Xti |22 ] + 4κ+ |x|22 <∞, (41)
using (P3) of Lemma E.16. By (P2) of Lemma E.16, the function θ 7→ h(θ) is continuous. Therefore, we can apply Lemma E.15, yielding that almost-surely for N →∞ the function
θ 7→ 1 N N∑ j=1 h(θ, ξj) = Φ̂N (θ) (42)
converges uniformly on Θ̃M to
θ 7→ EP×P̃[h(θ, ξ1)] = Φ(θ). (43)
We deduce from Lemma E.14 that d(θminM,N ,Θ min M ) → 0 a.s. when N → ∞. Then there exists a sequence (θ̂minM,N )N∈N in Θ min M such that |θminM,N − θ̂minM,N |2 → 0 a.s. for N → ∞. The uniform continuity of the random functions θ 7→ Y θt on Θ̃M implies that |Y θminM,N t − Y θ̂minM,N t |2 → 0 a.s. when N →∞ for all t ∈ [0, T ]. By continuity of F this yields |h(θminM,N , ξ1)− h(θ̂minM,N , ξ1)|2 → 0 a.s. as N →∞. With (41) we can apply dominated convergence which yields
lim N→∞
EP×P̃ [ |h(θminM,N , ξ1)− h(θ̂minM,N , ξ1)|2 ] = 0.
Since for every integrable random variableZ we have 0 ≤ |E[Z]|2 ≤ E[|Z|2] and since θ̂minM,N ∈ ΘminM we can deduce
lim N→∞ Φ(θminM,N ) = lim N→∞
EP×P̃ [ h(θminM,N , ξ1) ] = lim N→∞ EP×P̃ [ h(θ̂minM,N , ξ1) ] = Φ(θminM ). (44)
Now by triangle inequality,
|Φ̂N (θminM,N )− Φ(θminM )| ≤ |Φ̂N (θminM,N )− Φ(θminM,N )|+ |Φ(θminM,N )− Φ(θminM )|. (45)
(42), (43) and (44) imply that both terms on the right hand side converge to 0 when N →∞, which finishes the proof of the first part of the Theorem.
We define N0 := 0 and for every M ∈ N NM := min { N ∈ N | N > NM−1, |Φ(θminM,N )− Φ(θminM )| ≤ 1M + |Φ(θ min M )−Ψ(X̂)| } ,
which is possibly due to (44). Then Theorem E.2 implies that
|Φ(θminM,NM )−Ψ(X̂)| ≤ 1 M + 2|Φ(θ min M )−Ψ(X̂)|
M→∞−−−−→ 0.
Therefore, we can apply the same arguments as in the proof of Theorem E.2 (starting from (38)) to show that
EP×λk [∣∣∣X̂ − Y θminM,NM ∣∣∣
2
] ≤ c c̃ ( Φ(θminM,NM )−Ψ(X̂) )1/2 M→∞−−−−→ 0, for every 1 ≤ k ≤ K.
Corollary E.17. In the setting of Theorem E.13, we also have that (P× P̃)-a.s.
Φ(θminM,NM ) M→∞−−−−→ Ψ(X̂) and Φ̂ÑM (θ min M,ÑM ) M→∞−−−−→ Ψ(X̂),
where (ÑM )M∈N is another increasing sequence in N.
Proof. The first convergence result was already shown in the proof of Theorem E.13 and the second one can be shown similarly, when defining ÑM by Ñ0 := 0 and for every M ∈ N
ÑM := min { N ∈ N | N > ÑM−1, |Φ̂N (θminM,N )− Φ(θminM )| ≤ 1M + |Φ(θ min M )−Ψ(X̂)| } ,
which is possibly due to (45).

heading: E.3 DISCUSSION ABOUT OPTIMAL WEIGHTS
text: In Theorem E.2 and E.13, the focus lies on the convergence analysis under the assumption that optimal weights are found. Below we discuss, why this assumption is not restrictive in theory.
Global versus local optima. The assumption that the optimal weights are found, is typical for a convergence analysis of a neural network based algorithm, since the objective function is highly complex and non-convex with respect to the weights. In particular, it is well known that the standard choice of (stochastic) gradient descent optimization methods do in general only find local and not global minima. Since the difference between any local minimum and the global minimum can not generally be bounded, it is unrealistic to hope for a theoretical proof of convergence with respect to such optimisation schemes. On the other hand, global optimization methods as for example simulated annealing provably convergence (in probability) to a global optimum (Locatelli, 2000; Lecchini-Visintini et al., 2008). Hence, combining those with our result, convergence in probability of our model output to the conditional expectation can be established, without the assumption that the optimal weights are found. However, these global optimization schemes come at the cost of much slower training compared to (stochastic) gradient descent methods when applied in practice. Moreover, several works have focused on showing that most local optima of neural networks are nearly global, see for example (Feizi et al., 2017) and the related work therein. Hence, using (stochastic) gradient descent optimization methods likely yield nearly globally optimal weights much more efficiently. In our case, this is also supported by our empirical convergence studies in Section 6.3.

heading: F EXPERIMENTAL DETAILS
text: All implementations were done using PyTorch. The code is available at https://github.com/ HerreraKrachTeichmann/NJODE.
F.1 IMPLEMENTATION DETAILS
Dataset. For each of the SDE models (Black-Scholes, Ornstein-Uhlenbeck, Heston) a dataset was generated by samplingN = 20′000 paths of the SDE using the Euler-scheme. We used an equidistant time grid of mesh 0.01 between time 0 and T = 1. Independently for each path, observation times were sampled from Pt, by using each of the grid points with probability 0.1 as an observation time. In particular, n ∼ Bin(100, 0.1), t0 = 0 and the observation times {ti}1≤i≤n were chosen uniformly on the time grid. Hence, 10% of the grid points were used on average. This way, n and ti are defined as a discretized version of those given in Example A.2 where n is binomially distributed and ti chosen uniformly on [0, T ]. For each of these datasets, the samples were used in a 80%/20% split for training and testing. The SDEs of the dataset models and the chosen parameters are described below.
• Black-Scholes: – SDE: dXt = µXtdt+ σXtdWt, where W is a 1-dimensional Brownian motion – conditional expectation: E(Xt+s|Xt) = Xteµs
– used parameters: µ = 2, σ = 0.3, X0 = 1 • Ornstein-Uhlenbeck:
– SDE: dXt = −k(Xt −m)dt+ σdWt, where W is a 1-dimensional Brownian motion – conditional expectation: E(Xt+s|Xt) = Xte−ks +m ( 1− e−ks ) – used parameters: k = 2,m = 4, σ = 0.3, X0 = 1
• Heston: – SDE: for W and Z 1-dimensional Brownian motions
dXt = µXtdt+ √ vtXtdWt
dvt = −k(vt −m)dt+ σ √ vtdZt
– conditional expectation7: E(Xt+s|Xt) = Xteµs
– used parameters: µ = 2, σ = 0.3, X0 = 1 k = 2, m = 4, v0 = 4, ρ = Corr(W,Z) = 0.5
Architecture. In our experiments we choose the dimension of the latent variable to be dH = 10. For fθ1 , g̃θ3 and ρ̃θ2 we use 2-hidden-layer feed-forward neural networks, with 50 nodes in each hidden layer and tanh activation functions. Then the neural networks gθ3 and ρθ2 are defined as residual versions of g̃θ3 and ρ̃θ2 , by adding a residual shortcut between the input and the output of the neural networks. Dropout was applied after each non-linearity with a rate of 0.1. To scale the possibly unbounded inputs of the neural networks to a bounded hypercube, we applied tanh component-wise to x and h in every neural network. This was done, because neural networks sometimes become unstable when their inputs become large. To solve the ODE of the neural ODE part, the simple Euler-method was used.
Training. The neural networks were trained using the Adam optimizer (Kingma & Ba, 2014) with a learning rate of 0.001 and weight decay 0.0005 for 200 epochs using a batch size of 200. A random initialization was used and no hyper-parameter optimization was needed.
Further training results. Further training results on test samples are shown in Figure 8, 9, 10.

heading: F.2 EXPERIMENTS ON OTHER DATASETS
text: We test our framework on additional synthetic datasets.
7see (Rujivan & Zhu, 2012, Equation 2.9)

heading: F.2.1 HESTON MODEL WITHOUT FELLER CONDITION
text: If the Feller condition 2 km ≥ σ2
is satisfied in the Heston model, it is known that the variance process vt is always strictly bigger than 0 and that the Euler-scheme works well to sample from it. However, if the Feller condition is not satisfied, the variance process can touch 0, where the process is reflected deterministically. Then the Euler-scheme can not longer be used to sample from the model. Moreover, this situation is generally considered as more delicate. Close to 0 the distribution of vt behaves differently than further away of 0 (Andersen, 2007, Section 3). As explained by Andersen (2007), there are differently well performing sampling schemes for a Heston model where the Feller condition is not satisfied. We use the simplest to implement, which is a slight extension of the Euler scheme, where values of vt below 0 are replaced by 0 (Andersen, 2007, Section 2.3). Although there is empirical evidence, that the resulting sampling distribution of vt close to 0 is not correctly replicating the true distribution (Jean-François et al., 2015, Figure 2,3), this method already produces sufficiently good sample paths.
Dataset. Heston model, sampled as in Section F.1 with the extension described above. Used parameters: µ = 2, σ = 3, X0 = 1 k = 2, m = 1, v0 = 0.5, ρ = Corr(W,Z) = 0.5. Hence the Feller condition and also the weaker condition 4 km ≥ σ2 discussed in (Jean-François et al., 2015, Section 3.2), are both not satisfied. We produce two datasets, a 1-dimensional one similar to before, where only X is stored and a 2-dimensional one, where both X and v are stored, hence also v is a target for prediction. Note that v has the same conditional expectation as the Ornstein-Uhlenbeck SDE. In the 2-dimensional dataset, X and v are always observed at the same time.
Architecture & Training. Same as in Section F.1, but with batch size 100.
Results. The model learns to replicate the true conditional expectation process, which is analytic, hence not effected by the sampling scheme. In particular, we see that our model is very robust, since even in the delicate case where the Feller condition is not satisfied and a sampling scheme is used, that does not perfectly replicate the true distribution, our model still works well. Due to the very similar results, in Figure 11 we only show plots on test samples of the 2-dimensional dataset, where X and v are predicted.

heading: F.2.2 DATASET WITH CHANGING REGIME
text: Dataset. To be able to evaluate the performance of our model, we test a change of regime by combining two synthetic datasets. On the first half of the time interval [0, 0.5] we use the OrnsteinUhlenbeck and on the second half [0.5, 1] the Black-Scholes model. The Black-Scholes process takes as starting point, the last point of the Ornstein-Uhlenbeck process. We use the same hyperparameters for the dataset generation as in Section F.1, except that we set the parameter m = 10 in the Ornstein-Uhlenbeck model to make the two parts act on similar scales.
Architecture & Training. Same as in Section F.1, but with batch size 100. Moreover, we used 100 neurons in each hidden layer, to account for the more complicated setting, where also a time dependence has to be learnt.
Results. In the plots on test samples of the dataset shown in Figure 12 we see that our model correctly learns the change of regime.

heading: F.2.3 DATASET WITH EXPLICIT TIME DEPENDENCE
text: Dataset. We use the Black-Scholes datasets of Section F.1, where we replace the constant µ by the time dependent constant α2 (sin(βt) + 1), for α, β > 0. The conditional expectation changes accordingly. For the data generation we use the same hyper-parameters as in Section F.1, except that we use instead of µ the parameter α = 2 and β ∈ {2π, 4π}. Architecture & Training. Same as in Section F.1, but with batch size 100. Moreover, we used 400 neurons in each hidden layer, to account for the more complicated setting, where an explicit time dependence has to be learnt.
Results. We show plots on test samples of the datasets in Figure 5. We see that the model learns to adapt to the time-dependent coefficients.

heading: F.3 DETAILS ON CONVERGENCE STUDY
text: Evaluation metric. For the sampling time grid with equidistant step size ∆t := Tν , ν ∈ N, on [0, T ] and the true and predicted conditional expectation for path j ∈ N, X̂j and Y j respectively, we define the evaluation metric as
eval(X̂, Y ) := 1
N2 N2∑ j=1 1 ν + 1 ν∑ i=0 ( X̂ji∆t − Y j i∆t )2 , (46)
where N2 is the number of test samples, and j accordingly iterates over the paths in the test set.
Increasingly big training sets. We use the following procedure to create increasingly big training sets, while keeping the exactly same test set for evaluation. Out of the initial 20′000 paths, we take N2 := 4′000 paths, which are fixed as the testing set. Out of the remaining 16′000 paths, we randomly choose N1 training paths for N1 ∈ {200, 400, 800, 1′600, 3′200, 6′400, 12′800}. Increasing neural network sizes. The increasingly big neural networks are defined as follows. For all involved networks, we use the feed-forward 2-layer architecture with tanh activations (cf. Appedix F.1), where each hidden layer has the same size M for M ∈ {10, 20, 40, 80, 160, 320}. Results on Black-Scholes dataset. In accordance with the theoretical results in Theorem E.2 and E.13, we see that the evaluation metric decreases whenN1 andM increase (Figure 13). It is important to notice, that already a quite small number of samples can lead to a good approximation of the conditional expectation, if the network is big enough. In particular, the Monte Carlo approximation of the theoretical loss function is good already with a few samples. On the other hand, even with a large number of samples, the evaluation metric does not become so small, if the network size is not big enough. From a practitioners point of view, this is good news, since increasing the network size is often much easier than collecting more training data.
Results on Ornstein-Uhlenbeck and Heston dataset. We get very similar results on the Ornstein Uhlenbeck (Figure 14) and Heston dataset (Figure 15). Similar as for Black Scholes, also for Ornstein-Uhlenbeck the training size N1 is not as important as the network size M . For all network sizes, increasing the training size further than 1600 hardly changes the performance, while increasing M is crucial to get better performance. In contrast to this, for the Heston dataset we see that a large
number of training samples is more important to get a smaller convergence metric. This reflects the fact, that the Heston dataset is more complex and more difficult to learn.

heading: F.4 DETAILS ON COMPARISON TO GRU-ODE-BAYES
text: GRU-ODE-Bayes (Brouwer et al., 2019). To the best of our knowledge, this is the neural network based method with the most similar task to ours. In particular, this continuous time model is trained to learn the unknown temporal parameters of a normal distribution, best describing the conditional distribution of X given the previous observations. This distribution is given by the Fokker-Planck equation. Brouwer et al. (2019) outlined, that their model can exactly represent the Fokker-Planck dynamics of the Ornstein-Uhlenbeck process, since the corresponding distribution is Gaussian.
Implementation of GRU-ODE-Bayes. We use the code of the official implementation of (Brouwer et al., 2019)8 and slightly adjust it for our purpose. In particular, we do not use incomplete observations, hence the input mask used for this task has always only 1-entries. Furthermore, we slightly changed the scheme how the time steps are taken, to be the same as in our implementation, so that comparisons can be made. Besides these minor changes, the original model is used and
8https://github.com/edebrouwer/gru_ode_bayes
trained on all 3 datasets (Black-Scholes, Ornstein-Uhlenbeck and Heston). We tried out all combinations of the following parameters and always chose the best performing one for comparison to our model: impute ∈ {True,False}, logvar ∈ {True,False}, mixing ∈ {0.0001, 0.5} and hiddensize ∈ {50, 100}, whereby phidden and prephidden were chosen to be equal to hiddensize. The first parameter choices were the ones used in the official implementation. The model was always trained using the Euler-method for solving ODEs and with dropout = 0.1, to be comparable to our implementation. Furthermore, we always used the full GRU-ODE-Bayes implementation, since it should be the more powerful one. For the comparison, we only use the estimated mean of the normal distribution (the estimated variance is not used), which is precisely the estimated conditional expectation of the model.
Implementation of NJ-ODE. Same as described in Section F.1. For each dataset only one model was trained, since we already saw the convergence properties before and wanted to have a qualitative comparison to GRU-ODE-Bayes. In particular, we did not try to optimize our hyper-parameters for best performance (e.g. by choosing different number of layers or neurons or different activation functions) and used about 10K trainable parameters compared to the best performing GRU-ODEBayes models of our study which used 112K trainable parameters.
Datasets. Same as described in Section F.1. The train and test sets were fixed to be the same for all trained models.
Training. All models were trained using the Adam optimizer (Kingma & Ba, 2014) with a learning rate of 0.001 and weight decay 0.0005 for 100 epochs using a batch size of 20. A random initialization was used.

heading: F.4.1 FURTHER RESULTS AND DISCUSSION OF THE COMPARISON
text: Evolution of losses and evaluation metric during training. In Figure 16 and 17 we see, that already after a few epochs the NJ-ODE model finds close to optimal weights for the given network size on the Black-Scholes and Orstein-Uhlenbeck dataset and oscillates around this optimum. The Heston dataset is considerably more difficult and the model slowly converges to close-to-optimal weights for the given network size. In comparison to this, it takes the GRU-ODE-Bayes model longer to converge to its close-to-optimal weights on the Black-Scholes and Orstein-Uhlenbeck dataset. Moreover, the model does not converge to close-to-optimal weights on the Heston dataset, but rather oscillates between bad weights. Due to some very large outliers, this is not directly visible in the plot, but can be deduced from Table 6.4.
Comparison of predicted paths. For each dataset we show 5 paths that were predicted with NJODE and with GRU-ODE-Bayes, first at the optimal epoch, i.e. where the test loss was minimal during training (Figures 18, 19, 20), and then at the last epoch (Figures 21, 22, 23). For the sake of comparison, in each row the performance on the same test sample is shown. The results are very similar for NJ-ODE and GRU-ODE-Bayes on the Black-Scholes and Ornstein-Uhlenbeck dataset, with sometimes the one and sometimes the other having a slightly better prediction. On the Heston
0 20 40 60 80 100
0.1
0.2
tra in
_lo ss
0 20 40 60 80 100 0.15
0.16 0.17 ev al _lo ss
0 20 40 60 80 100 epoch
0.01
0.02
ev al
_m et
ric 0 20 40 60 80 100
0.0075
0.0100
tra in
_lo ss
0 20 40 60 80 100
0.008
0.009
ev al
_lo ss
0 20 40 60 80 100 epoch
0.0005
0.0010
0.0015
ev al
_m et
ric
0 20 40 60 80 100 0
500
tra in
_lo ss
0 20 40 60 80 100 47
48
ev al
_lo ss
0 20 40 60 80 100 epoch
2
3
ev al
_m et
ric
Figure 16: NJ-ODE on Black-Scholes, Ornstein-Uhlenbeck and Heston (from left to right). Blue (1st row): training loss, orange (2nd row): evaluation loss, green (3rd row): evaluation metric.
0 20 40 60 80 100
50
0
50
tra in
_lo ss
0 20 40 60 80 100
7500
5000
2500
ev al
_lo ss
0 20 40 60 80 100 epoch
0.00
0.01
ev al
_m et
ric
0 20 40 60 80 100 0
5000
10000
tra in
_lo ss
0 20 40 60 80 100
0
200000
ev al
_lo ss
0 20 40 60 80 100 epoch
0.00
0.01
0.02
ev al
_m et
ric 0 20 40 60 80 100
0
50000
tra in
_lo ss
0 20 40 60 80 100 0
1
ev al
_lo ss
1e7
0 20 40 60 80 100 epoch
0
500
ev al
_m et
ric
Figure 17: GRU-ODE-Bayes’ best performing models on Black-Scholes, Ornstein-Uhlenbeck and Heston (from left to right). Blue (1st row): training loss, orange (2nd row): evaluation loss, green (3rd row): evaluation metric.
dataset, the results of NJ-ODE are good predictions, being correct whenever a new observation is made and not too far away even after longer periods without observations. On the other hand, GRU-ODE-Bayes is not even correct at times of new observations and learns an incorrect behaviour in between observations (e.g. making kinks where there should not be any). For the last epoch, we see that this malfunction is amplified.

heading: F.4.2 FURTHER EXPERIMENTS
text: Since the results of GRU-ODE-Bayes were unexpectedly bad on the Heston dataset, we performed additional tests. In particular, we retrained the best combinations of GRU-ODE-Bayes with the larger batch size 100. For smaller versions of the network (M), i.e. with hiddensize = prephidden = phidden = 50 this stabilized the training, such that the models converged. Still, even at the best epoch, the models suffered from the same difficulties as explained in Section F.4.1. Increasing the hiddensize by factor 2 to 100 (L), made the training unstable again, and the models did not converge any more. This gives more empirical evidence, that GRU-ODE-Bayes can not be reliably trained on more complex datasets, where the target conditional distributions differs to much from a normal distribution. In contrast to this, retraining NJ-ODE with batch size 100, once for the smaller version described in Section F (S), once for a larger version with 100 instead of 50 neurons in all hidden layers (M) and once with 200 neurons (L), yielded the expected results of better performance with larger networks. In particular, there is no instability for larger networks. In Table F.4.2 we show results of our model and of the best GRU-ODE-model for the given size.

heading: F.5 REAL WORLD DATASETS WITH INCOMPLETE OBSERVATIONS

heading: F.5.1 LOSS FUNCTION FOR INCOMPLETE OBSERVATIONS
text: In accordance with the self-imputation scheme, the loss function is adjusted to only use the nonimputed coordinates, i.e. if m is the random process in {0, 1}dX describing which coordinates are observed we have
Ψ(Z) := EP×P̃
[ 1
n n∑ i=1 (|mti (Xti − Zti)|2 + |mti (Zti − Zti−)|2) 2
] .

heading: F.5.2 CLIMATE FORECASTING DETAILS
text: Dataset. We use the publicly available United State Historical Climatology Network (USHCN) daily dataset (Menne et al., 2016) together with all pre-processing steps as they were provided by Brouwer et al. (2019). In particular, there are 5 sporadically observed (i.e. incomplete observations) climate variables (daily temperatures, precipitation, and snow) measured at 1′114 stations scattered over the United States during an observation window of 4 years (between 1996 and 2000) where each station has an average of 346 observations over those 4 years. For 5 folds, the data is split into train (70%), validation (20%) and test (10%) sets. The task is to predict the next 3 measurements after the first 3 years of observation. The mean squared error between the prediction and the correct values is computed on the validation and test set.
Baselines. We use the results reported in (Brouwer et al., 2019, Table 1) as baselines for our comparison and perform the exact same 5-fold cross validation using the same folds with the same train, validation and test sets. For completeness, we give all results of the table together with our results in Table 6.5. We only show the mean squared error (MSE), since our model does not provide the negative log-likelihood.
Implementation of NJ-ODE. We once use the architecture described in Section F.1 (S) and once use the exact same architecture, but with hidden size dH = 50 and 400 instead of 50 nodes in each hidden layer (L). To deal with the incomplete observations, the self-imputation scheme described in Section 6.5 is used. In particular, the data is self-imputed and passed together with the observation mask as input to the network ρ̃θ2 . Moreover, the network fθ1 uses the NJ-ODE output at the last observation time instead of the last observation as input.
Training. Both versions of NJ-ODE were trained using the Adam optimizer (Kingma & Ba, 2014) with a learning rate of 0.001 and weight decay 0.0005 for 100 epochs using a batch size of 100. A random initialization was used. The performance on the validation set was used for early stopping after the first 100 epochs, i.e. early stopping was possible at any epoch between {101, . . . , 200}. Results. The results of our model and all models reported in (Brouwer et al., 2019, Table 1) are shown in Table 6.5.

heading: F.5.3 PHYSIONET PREDICTION DETAILS
text: Dataset. We use the publicly available PhysioNet Challenge 2012 dataset (Goldberger et al., 2000) together with all pre-processing steps as they were described and provided by Rubanova et al. (2019). In particular, there are 41 features of 8000 patients that are observed irregularly over a 48 hour time period. The observations are put on a time grid with step size 0.016 hours leading to 3000 grid points. While in their paper Rubanova et al. (2019) say that they use 2880 grid points (i.e. minute wise) in their implementation they used 3000. Moreover, in contrast to what was written in the paper, the 4 constant features were not excluded in their implementation, hence we also keep them. Furthermore, we also rescale the time-grid to [0, 1] and normalize each feature as Rubanova et al. (2019) did it for training the model. The dataset is split with the same fixed seed into 80% training and 20% test set. In particular, no cross validation but only multiple runs with new random initializations are performed to be exactly comparable to the results reported by Rubanova et al. (2019). On the test set, the observation paths are split in half. The first half (first 1500 time steps) is used as input for the model, from which the second half (second 1500 time steps) should be forecast. The masked mean squared error (MSE) with a certain balancing procedure between the predictions and the true values is reported. For comparability we use the exact same implementation of the MSE as Rubanova et al. (2019).
Baselines. We compare the performance of our model to latent ODE on the extrapolation task (as described above) on physionet. As baseline for our comparison we use the results reported in (Rubanova et al., 2019, Table 5). For completeness, we show all extrapolation results of the table together with our results in Table 6.5. We shortly outline the different approach of latent ODE compared to our model for the given extrapolation task. Latent ODE also splits the training samples similar to the test samples in half, using the first half as input and the second half as target. It is trained as an encoder-decoder, encoding the observations in the first half and reconstructing (decoding) the second half. This falls in the standard supervised learning framework. In particular, this approach cannot straight forward be extended for online forecasting. Moreover, this approach might learn certain path dependencies. On the other hand, our model is trained as always, online forecasting after each observation until the next observation is made. Instead of splitting the training samples we use the entire path as input for our unsupervised training framework. Our model is based on the assumption that paths are Markov, therefore it cannot learn path dependencies, i.e. dependencies on more than just the last observation. However, by training the model also on the second half of the training samples, it learns the underlying behaviour there, which should be helpful for the extrapolation task.
Implementation of NJ-ODE. We use the architecture described in Section F.1, but with hidden size dH = 41. To deal with the incomplete observations, the self-imputation scheme described in Section 6.5 is used. In particular, the data is self-imputed and passed together with the observation mask as input to the network ρ̃θ2 . Moreover, the network fθ1 uses the NJ-ODE output at the last observation time instead of the last observation as input.
Training. The NJ-ODE was trained using the Adam optimizer (Kingma & Ba, 2014) with a learning rate of 0.001 and weight decay 0.0005 for 175 epochs using a batch size of 50. 5 runs with random initialization were performed over which the mean and standard deviations were calculated. In particular, these runs always used the same training and test set, specified by the same random seed as in thee implementation of Rubanova et al. (2019).
Results. The results of our model and all models reported in (Rubanova et al., 2019, Table 5, extrapolation) are shown in Table 6.5. We report the minimal MSE on the test set during the 175 epochs, since it was not differently specified in (Rubanova et al., 2019). However, if the MSE of the epoch is used where the training loss is minimal, the results are nearly the same with 1.986± 0.058 (×10−3). Moreover, we trained a larger model for 120 epochs, where 200 nodes were used instead of 50 (187′323 parameters in total), leading to slightly better results of 1.934± 0.007 (×10−3) (at minimal MSE) and 1.982± 0.027 (×10−3) (at minimal training loss).

Review summaries:
This paper introduces Neural Jump Ordinary Differential Equations as a method for learning models of continuous-time stochastic processes sampled at random time epochs. Specifically, the paper studies the problem of estimating the marginal conditional expectation (i.e., the L2 optimal approximation conditional on the available information) by estimating an auxiliary stochastic differential equation, parameterized by  neural networks, that approximates the conditional expectation of the process of interest at each point in time. The neural networks are trained by using a “randomized” mean squared-loss objective. The main theoretical results in the paper include asymptotic consistency of the optimal objective value in the limit of a large neural network, as well as consistency of a Monte Carlo sample average estimator of the value. The paper also establishes the L2 convergence of the estimated auxiliary solution to the marginal conditional expectation.

The authors propose a method for learning the conditional expectation of stochastic process in an online fashion. The paper bears a considerable theoretical treatment, derived from the stochastic filtering literature, which is present both in the main body of the paper and the appendix. Besides the model, the paper also aims to provide a theoretical justification of the convergence of their method.