Neural SDDEs: Stochastic Delay Differential Equations

• Neural SDDEs are continuous-time stochastic models where neural networks parameterize drift and diffusion while incorporating delayed, history-dependent dynamics.
• They split uncertainty into aleatoric and epistemic components, providing immediate uncertainty estimates without needing repeated sampling.
• Applications include irregular astronomical time series and multivariate forecasting, with robust numerical methods ensuring stability and improved convergence.

Neural Stochastic Delay Differential Equations (Neural SDDEs) are continuous-time stochastic dynamical models in which drift and diffusion terms are parameterized by neural networks while the dynamics depend explicitly on delayed states or selected history segments. In the formulations studied in "Delay-SDE-net: A deep learning approach for time series modelling with memory and uncertainty estimates" (Eggen et al., 2023) and "Modeling Irregular Astronomical Time Series with Neural Stochastic Delay Differential Equations" (Oh et al., 24 Aug 2025), Neural SDDEs combine neural parameterization of dynamics, stochastic evolution to represent uncertainty or noise, and explicit delay dependence to capture memory effects and long-range temporal interactions. Their state is naturally history-dependent rather than Markovian in the current point alone, so both analysis and computation are organized around delayed trajectories, segment processes, and numerical schemes that preserve access to past states.

1. Formal structure and relation to other neural differential equations

A central SDDE formulation with multiple delays is given by the segment process

$X_t(u)\coloneqq X(t+u), \qquad u\in[-\tau,0],$

together with a projection onto selected delay points

$$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$

The corresponding delayed stochastic system is written as

$X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$

where $\Pi_1$ and $\Pi_2$ may select different sets of delays for drift and diffusion (Eggen et al., 2023). This makes the model non-Markovian in the current state alone and places the relevant state on a path space.

The latent-state astronomical formulation uses a fixed delay $\tau>0$, an initial history function $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$, a drift $\gamma$ depending on current and delayed latent states, and a diffusion $\sigma$ parameterized by neural networks with weights $\theta_\gamma$ and $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$0. The paper emphasizes that $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$1 is not a literal physical delay extracted from the data, but a model “lookback window” that enables the network to learn delayed interactions. It also states the functional-analytic interpretation that the SDDE state is not just a point $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$2, but a history segment $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$3 in the function space $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$4, so $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$5 is a $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$6-valued Markov process (Oh et al., 24 Aug 2025).

The relation to nearby model classes is explicit.

Model class	Defining dependence	Property emphasized
Neural ODE	Current state only	Deterministic and memoryless
Neural SDE	Current state plus diffusion	Stochastic but Markovian in the state
Neural DDE	Current and delayed states	Delayed and deterministic
Neural SDDE	Current and delayed states plus diffusion	Stochasticity plus delay

In this taxonomy, Neural SDDEs are positioned as the combination of the two missing pieces: stochasticity plus delay (Oh et al., 24 Aug 2025). A plausible implication is that they are most natural when the signal is simultaneously noisy, irregularly sampled, and history-dependent.

2. Neural parameterization, latent augmentation, and uncertainty decomposition

The most explicit neural construction in the time-series setting is the Delay-SDE-net, which uses three neural components: $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$7 for drift, $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$8 for aleatoric diffusion, and $$\Pi(\eta)\coloneqq (\eta(u_1),\ldots,\eta(u_p))\in\mathbb{R}^{dp}.$$9 for epistemic diffusion, with

$X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$0

The continuous-time Neural SDDE is

$X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$1

In the theoretical analysis, $X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$2, $X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$3, and $X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$4 are taken as two-layer neural networks, while the empirical studies do not restrict them to two layers. A two-layer network is written as

$X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$5

and the paper uses the Barron-space framework to control approximation error for such networks (Eggen et al., 2023).

A distinctive architectural choice is the split of uncertainty into aleatoric and epistemic components. Aleatoric uncertainty captures intrinsic randomness in the data-generating process and is trained from prediction residuals after the drift has been fitted. Epistemic uncertainty is modeled as uncertainty due to limited knowledge or OOD behavior, and $X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$6 is trained as a classifier separating in-distribution and out-of-distribution data. The final uncertainty prediction is therefore $X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$7, and the model can report both components instantly. The paper emphasizes that this is a single deterministic method for uncertainty quantification: once the three networks are trained, prediction and uncertainty are obtained immediately, unlike Bayesian or ensemble methods that require repeated sampling or multiple forward passes (Eggen et al., 2023).

The astronomical Neural SDDE augments the latent dynamics with a controlled-path representation inspired by Neural CDEs. It defines an extended state $X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$8, where $X(t)= \begin{cases} \eta(0)+\displaystyle\int_0^t f\big(s,\Pi_1(X_s)\big)\,ds +\int_0^t g\big(s,\Pi_2(X_s)\big)\,dW(s), & t\ge 0,\[1ex] \eta(t), & -\tau\le t<0, \end{cases}$9 is a continuous interpolation of the raw irregular observations, and $\Pi_1$0 combines the current time, latent state, and observed control path. The hidden state at $\Pi_1$1 is initialized by $\Pi_1$2, and the authors note that in sparse or noisy settings, a simple static initialization $\Pi_1$3 on $\Pi_1$4 can be preferable for stability and ease of implementation (Oh et al., 24 Aug 2025).

Training procedures differ accordingly. Delay-SDE-net uses a sequential and deliberately decoupled procedure: first train the drift $\Pi_1$5, then the aleatoric diffusion $\Pi_1$6, and finally the epistemic diffusion $\Pi_1$7, using the soft Brownian offset method to generate OOD samples (Eggen et al., 2023). The astronomical model uses supervised cross-entropy for regular classification and a maximum-softmax-probability outlier score for novelty detection, with $\Pi_1$8 given as an example threshold (Oh et al., 24 Aug 2025).

3. Discretization, approximation theory, and solver design

Delay-SDE-net adopts a simplified data-generating SDDE in which the diffusion is state-independent in the continuous formulation and depends only on the initial path. The paper states that this allows the Milstein correction term to vanish, so in the discrete approximation the scheme reduces to Euler–Maruyama. Its theoretical analysis decomposes the total error into approximation error from replacing the true SDDE by the continuous-time neural SDDE and discretization error from replacing the continuous-time neural SDDE by its numerical scheme. Under the assumptions that $\Pi_1$9, the initial path $\Pi_2$0 is of bounded variation and $\Pi_2$1-Hölder continuous, and the neural approximants $\Pi_2$2, bounded with bounded derivatives, the paper proves

$\Pi_2$3

and

$\Pi_2$4

In numerical convergence experiments with synthetic SDDE data, a log-log regression yields a convergence rate

$\Pi_2$5

which the authors state is better than the classical Euler–Maruyama rate $\Pi_2$6, partly due to the state-independent diffusion structure in their setup (Eggen et al., 2023).

For general delayed stochastic systems, the numerical-analysis literature establishes several solver backbones that are directly about classical SDDEs rather than neural networks. For neutral SDDEs driven by Brownian motions, the Euler–Maruyama convergence rate is one half in strong $\Pi_2$7-sup norm under polynomial-growth assumptions on the neutral, drift, and diffusion terms (Ji et al., 2015). For SDDEs with a non-global Lipschitz drift term and a non-constant delay, an adaptive Euler–Maruyama method with a left-neighbor approximation of the delayed state is strongly convergent, and under stronger assumptions its rate is

$\Pi_2$8

(Liu et al., 2024). For nonlinear SDDEs under generalized Khasminskii-type conditions, truncated Euler–Maruyama schemes yield strong convergence and $\Pi_2$9-order rates, while also preserving exponential stability properties in mean square and almost surely (Guo et al., 2017, Song et al., 2020).

A complementary line of work treats rougher coefficient classes. For Carathéodory-type drift coefficients that are only measurable in time and Lipschitz in the current state but merely Hölder continuous in the delay variable, a randomized Euler scheme attains the blockwise rate

$\tau>0$0

and the paper interprets the factor $\tau>0$1 as degradation caused by error propagation through successive delay blocks (Przybyłowicz et al., 2023). On the long-time side, strong convergence on the infinite horizon can be obtained with constants independent of time; for example, the truncated Euler–Maruyama analysis proves

$\tau>0$2

under polynomial-growth and contractivity assumptions, and segment-process convergence is established for backward Euler–Maruyama (Wang et al., 20 May 2025). This suggests that Neural SDDE training and rollout inherit a substantial solver-design question: the choice between explicit, adaptive, truncated, randomized, or implicit schemes is analytically tied to delay structure, regularity, and long-time stability rather than being a purely implementation-level detail.

4. Non-Markovianity, segment processes, and density evolution

A defining analytical fact is that SDDEs are generally not Markovian in the current state alone. The density-evolution paper states that Fokker–Planck equations do not exist for SDDEs because the solutions are not Markov processes and have no corresponding infinitesimal generators. For the discrete-delay SDDE

$\tau>0$3

the core construction lifts the delayed system to a coupled Markovian system on a larger state space by segmenting time into delay-length blocks, so that the lifted transition density $\tau>0$4 satisfies a classical Fokker–Planck equation in the enlarged phase space (Zheng et al., 2016).

The resulting density formulas are recursive rather than local in $\tau>0$5. The one-time density $\tau>0$6 of the SDDE solution is represented through the transitional densities $\tau>0$7 of the lifted system, with different expressions on the first delay interval, at exact multiples of the delay, and on intermediate intervals. The paper’s governing equation for the SDDE density is therefore not a single local PDE in $\tau>0$8, but a recursive representation built from Markovian transition densities of the lifted system (Zheng et al., 2016).

This perspective is directly compatible with Neural SDDE formulations that treat the history segment as the state variable. The astronomical paper states that $\tau>0$9 is a $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$0-valued Markov process, while the density-evolution paper argues that, once a Neural SDDE is specified, its density evolution can in principle be studied through the lifted Fokker–Planck system with neural-network drift and diffusion terms. The stated possibilities include density-based training objectives, uncertainty quantification, model analysis at the density level, and numerical inference by solving associated Fokker–Planck PDEs block by block when smoothness and ellipticity assumptions are available (Oh et al., 24 Aug 2025, Zheng et al., 2016). A plausible implication is that probabilistic analysis of Neural SDDEs is naturally more path-space oriented than the corresponding analysis of Neural SDEs.

5. Stability transfer, long-time behavior, and numerical equivalence

For delayed stochastic systems, stability is not only a property of the continuous-time model; it can be transferred to and from the numerical method under strong convergence and moment-boundedness assumptions. A general theorem for autonomous SDDEs proves that $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$1th moment exponential stability of the SDDE is equivalent to $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$2th moment exponential stability of its numerical method, provided the numerical method has bounded finite-time $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$3th moments and strong convergence on finite time intervals. The paper treats the truncated Euler–Maruyama method as an example and emphasizes that the result does not depend on the detailed structure of the numerical method beyond those generic hypotheses (Bao et al., 2019).

This stability-equivalence principle extends to more specialized uncertainty models. In the $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$4-framework, practical exponential stability in $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$5-th moment is shown to be equivalent among a $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$6-SDDE, an auxiliary nondelay $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$7-SDE, and their Euler–Maruyama approximations, provided the delay $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$8 or the step size $\phi\in \mathcal{C}([-\tau,0),\mathbb{R}^{d_z})$9 is sufficiently small (Lu, 2024). The paper’s conclusion is that careful simulations can be used to examine the practical exponential stability of the underlying delayed system under some reasonable assumptions.

Long-time approximation introduces the segment process as the correct state variable. The infinite-horizon convergence paper proves that one can obtain bounds of the form

$\gamma$0

with $\gamma$1 independent of time, and then use these results to approximate invariant measures of SDDEs through the numerical segment process (Wang et al., 20 May 2025). This suggests that for Neural SDDEs, long-run empirical distributions and stability properties are meaningful only after specifying both the continuous delayed model and a solver whose segment-process dynamics satisfy the needed moment and contractivity estimates.

6. Empirical regimes, benchmark behavior, and open analytical constraints

The two principal application regimes represented in the literature are multivariate time-series forecasting with uncertainty decomposition and irregular astronomical sequence modeling. Delay-SDE-net is evaluated against a VAR(4) model and the original SDE-net on simulated nonlinear two-dimensional SDDE data with memory and heteroskedasticity. The paper reports that Delay-SDE-net is consistently best for drift prediction RMSE, aleatoric variance RMSE, and epistemic OOD detection via ROCAUC, while noting that VAR outperforms SDE-net on the low-dimensional simulated problem because it captures lagged dependence directly, and that SDE-net performs worse because it lacks explicit memory (Eggen et al., 2023).

On real data, Delay-SDE-net is tested on daily zonal-mean stratospheric U wind and temperature at $\gamma$2 hPa and $\gamma$3N, trained on 1979–2008 and tested on 2009–2018. For 1–7 day ahead prediction, the paper states that Delay-SDE-net gives the best RMSE for wind prediction and the best uncertainty calibration, while SDE-net performs worst and VAR is competitive but generally inferior. It further highlights that $\gamma$4 captures seasonal variance well, $\gamma$5 rises during unusual events like sudden stratospheric warmings, and the combined uncertainty gives reasonable 95% confidence intervals (Eggen et al., 2023).

The astronomical study evaluates Neural SDDEs on the LSST dataset derived from the PLAsTiCC challenge, with 4925 instances, 6 input dimensions corresponding to the ugrizy filters, 36 sequences, and 14 classes, using a 70:15:15 train/validation/test ratio. Across four scenarios—regular fully labeled classification, regular training with irregular testing at 50% missing rate, fully irregular training and testing with 50% missing rate, and novelty detection with partial labels—Neural SDDE is reported as best across all four scenarios. In regular classification it reaches 0.706 accuracy and 0.684 weighted F1; in the irregular testing scenario it achieves 0.644 accuracy and 0.631 weighted F1; in the fully irregular scenario it obtains 0.656 accuracy and 0.626 weighted F1; and for novelty detection it achieves AUROC 0.707 and AUPRC 0.932, outperforming Neural LSDE, Neural SDE, Neural DDE, SDE-Delay-Net, and a broad set of recurrent and continuous-time baselines (Oh et al., 24 Aug 2025).

The same paper reports that delayed models outperform corresponding no-delay versions overall, particularly under irregular sampling and missing-label settings, whereas under fully labeled and regularly sampled conditions performance is relatively insensitive to $\gamma$6. It also identifies several limitations: the model relies on a fixed delay parameter $\gamma$7, training is computationally heavier than simpler discrete-time baselines because it requires continuous-time integration, storage of past segments, and reverse-time gradient computation through both drift and diffusion, and the paper does not present a detailed ablation isolating each neural subnetwork architecture independently (Oh et al., 24 Aug 2025). More generally, the classical numerical-analysis papers make clear that a generic neural parameterization does not automatically satisfy the growth, Lipschitz, monotonicity, or moment conditions used in SDDE convergence theorems. For neutral and jump systems in particular, rigorous transfer to neural parameterizations would require network-induced growth control, Lipschitz or monotonicity constraints, moment stability of the learned system, approximation-error decomposition, and, for jump models, integrability of neural jump maps against the Lévy measure (Ji et al., 2015). This suggests that the mathematically decisive unresolved issue for Neural SDDEs is not only expressivity, but the alignment of neural architectures with the regularity structure demanded by delay-stochastic analysis.