Neural ODEs: Does Continuous Time Beat a Delta-Time Feature?

Neural Ordinary Differential Equations (Chen et al., NeurIPS 2018) parameterize the derivative of a hidden state with a neural network and let an ODE solver compute the output. The model becomes continuous-depth and, more usefully for markets, continuous-time: the hidden state can be evaluated at any $t$, not only at the grid points your data happens to sit on.

That property has an obvious pitch for financial data, which arrives at random intervals. But the pitch hides a much weaker claim than it appears to make. A plain GRU handed an extra input feature — $\log(\Delta t)$ since the previous observation — also "knows" how much time passed. The real question is not whether continuous-time models can ingest irregular data. It is whether learned continuous dynamics between observations extract anything a delta-time feature does not already give you, at ten to a hundred times the inference cost.

This article sets up that experiment, gives the full toolchain to run it, and reports what the comparison would have to show to justify the machinery. The bar-construction side of the argument — that trade timing itself carries signal — is established with measured evidence across 17 bar types in beyond time bars; here the only consequence that matters is a modeling one: the solver evaluates the hidden state at arbitrary $t$, so no interpolation or padding is needed.

The Claim Under Test

Three nested hypotheses, each strictly stronger than the last:

1. H1 — An ODE-RNN fit on a raw irregular BTC trade stream beats a GRU fit on the same stream with no time information at all. This is nearly certain and nearly worthless: the ODE-RNN simply has information the baseline lacks.
2. H2 — The ODE-RNN beats the same GRU once the GRU is given $\log(\Delta t)$ as an input feature. This is the honest test. It is the claim the Neural ODE literature is usually read as making, and it is the one almost never demonstrated on financial data.
3. H3 — The H2 advantage survives the latency budget, i.e. it holds when the ODE-RNN is deployed with a fixed-step solver fast enough for live inference rather than adaptive dopri5.

Experiment Protocol

Same tick stream, same target, same tuning budget for both arms. No resampling to 1m bars — resampling destroys the exact structure under test.

A negative result on H2 is fully publishable and arguably the more valuable outcome — the same standard applied in the honest negative result and deflated Sharpe under multiple testing.

Secondary Experiment: CNF vs Student-t

If the ODE-RNN comparison is too expensive, the cheaper empirical anchor is distributional: fit a continuous normalizing flow on real BTC daily returns and compare its fitted density against a Student-t fit and against the empirical tail, reporting tail-quantile errors at the 1% and 5% levels. That connects directly to the distributional work already measured in GARCH volatility forecasting and asymmetric GARCH.

The rest of this article is the background and the toolchain needed to run the above.

Background: From ResNets to Continuous Dynamics

A residual network computes $h_{t+1} = h_t + f(h_t, \theta_t)$. Shrink the step size and increase the number of layers and you approach a continuous limit:

$\frac{dh(t)}{dt} = f(h(t), t, \theta)$

Instead of $T$ discrete layers with separate parameters, a single network $f$ specifies the instantaneous rate of change. The output at time $T$ solves an initial value problem:

$h(T) = h(0) + \int_0^T f(h(t), t, \theta) \, dt$

A black-box solver (Euler, Runge-Kutta, Dormand-Prince) computes the integral numerically, adaptively choosing step sizes from local error estimates.

The Adjoint Method

Backpropagating through every solver step stores all intermediate states, costing memory proportional to step count. Chen et al. instead solve a backward ODE with $O(1)$ memory. Define the adjoint state $a(t) = -\partial L / \partial h(t)$, which satisfies

$\frac{da(t)}{dt} = -a(t)^T \frac{\partial f(h(t), t, \theta)}{\partial h}$

and accumulate the parameter gradient along the backward pass:

$\frac{dL}{d\theta} = -\int_T^0 a(t)^T \frac{\partial f(h(t), t, \theta)}{\partial \theta} \, dt$

The backward pass solves an augmented system computing $h(t)$, $a(t)$ and $dL/d\theta$ simultaneously, running the solver in reverse from $T$ to $0$.

The trade-off: reconstructing $h(t)$ backwards accumulates numerical error, badly so for stiff or chaotic dynamics. Checkpointing is the middle ground — store $h(t)$ at a few intermediate times, recompute between them. Price processes are continuous semimartingales and moderately smooth, so the adjoint generally holds up; stiffness around microstructure events may force adaptive solvers or hybrids.

ODE-RNN: Continuous Hidden State Between Observations

The ODE-RNN is the architecture the main experiment turns on. Between observations the hidden state evolves under the ODE:

$h(t) = \text{ODESolve}(f_\theta, h(t_i), t_i, t)$

When observation $x_{i+1}$ arrives at $t_{i+1}$, a discrete update fires:

$h(t_{i+1}^+) = \text{RNNCell}(h(t_{i+1}^-), x_{i+1})$

The superscripts denote state just before and just after the observation. Between observations, learned continuous dynamics; at observations, new information.

The mechanism the H2 ablation probes is this: a long gap means the hidden state has evolved a long way under $f_\theta$ — decaying toward a baseline, or diverging — and the shape of that evolution is learned rather than supplied as a scalar. Whether that expressiveness pays for itself on real trade data is exactly what is unmeasured.

Natural fits beyond ticks: multi-asset portfolios where each asset has its own observation schedule and the latent state of correlated assets keeps evolving while only one is observed; and event-driven signals (news, earnings, macro releases) arriving at irregular times.

Neural SDEs: Adding Stochastic Components

Neural ODEs are deterministic and cannot represent the noise in a price path. The classical treatment — geometric Brownian motion, risk-neutral drift, and the ways constant-volatility assumptions fail against the smile — is covered in Black-Scholes options pricing. Neural SDEs replace the parametric form with a learned diffusion term:

$dh(t) = f(h(t), t, \theta) \, dt + g(h(t), t, \phi) \, dW(t)$

$f$ is the drift, $g$ the diffusion, $W(t)$ a Wiener process, both $f$ and $g$ neural networks.

Architecture: Drift Net and Diffusion Net

• Drift net $f_\theta$: expected trajectory — trend, mean reversion, momentum. Trained to minimize prediction error.
• Diffusion net $g_\phi$: noise magnitude, learned as a function of state. High in volatile regimes, low in calm ones.

The split mirrors classical quant finance: drift is the risk-neutral (or P-measure) dynamics, diffusion is the volatility surface.

Training Neural SDEs

The stochastic integral $\int g\,dW$ is not classically differentiable. Three approaches:

1. Pathwise gradients: reparameterization trick — sample Brownian paths, differentiate through the solver treating noise as fixed input.
2. Score matching: estimate $\nabla_h \log p(h)$ and train via denoising objectives — the same machinery used as a standalone generative model in diffusion models for crypto prediction, here demoted to a training option for the SDE.
3. Finite-dimensional distribution matching: match marginals at observation times rather than full path measures.

Pathwise is the practical default. torchsde implements Euler-Maruyama, Milstein and stochastic Runge-Kutta with autodiff support.

Learning the Volatility Surface

Classical models (Heston, SABR) impose parametric forms on the diffusion coefficient. A Neural SDE learns $g_\phi(S, t)$ as an arbitrary function:

$dS(t) = \mu_\theta(S(t), t) \, dt + \sigma_\phi(S(t), t) \, S(t) \, dW(t)$

This is a universal approximator for diffusion processes — any Ito process to arbitrary accuracy given enough capacity.

Latent ODEs for Missing and Asynchronous Data

The Latent ODE (Rubanova, Chen, Duvenaud, 2019) pairs a Neural ODE with a VAE: financial series are noisy observations of a smooth underlying process, learned in a low-dimensional latent space.

1. Recognition network: an ODE-RNN running backward through observations to produce $q(z_0 | x_{1:N})$.
2. Latent dynamics: $z(t) = z(t_0) + \int_{t_0}^{t} f_\theta(z(s), s)\, ds$.
3. Decoder: $\hat{x}(t) = \text{Decoder}(z(t))$, evaluable at any requested time.

The loss is the standard ELBO:

$\mathcal{L} = \mathbb{E}_{q(z_0)} \left[ \sum_{i=1}^{N} \log p(x_i | z(t_i)) \right] - \text{KL}(q(z_0 | x_{1:N}) \| p(z_0))$

Portfolio state estimation: from sparse asynchronous observations across assets, infer a continuous latent state capturing joint dynamics — essentially a nonlinear, learned version of the Kalman filter.

Missing data imputation: halts and weekend gaps get a smoothly interpolated latent trajectory; the decoder produces plausible paths through the gap with uncertainty from the VAE posterior.

Multi-frequency fusion: daily closes, intraday VWAP and tick data in one model without a shared time grid.

Continuous Normalizing Flows for Return Distributions

CNFs use a Neural ODE to transform a simple base distribution into a complex target. The transformation is $\frac{dz(t)}{dt} = f(z(t), t, \theta)$, and the log-density follows the instantaneous change of variables:

$\frac{\partial \log p(z(t))}{\partial t} = -\text{tr}\left(\frac{\partial f}{\partial z(t)}\right)$

Integrate state and log-density forward from $z(0) \sim p_0$ to get $z(T)$ and $\log p(z(T))$. The Jacobian trace is estimated with Hutchinson's stochastic estimator for scalability.

Crypto returns are leptokurtic and negatively skewed — measured on real data in GARCH volatility forecasting and asymmetric GARCH and the leverage effect — and a CNF learns that shape instead of assuming it. Conditioning the flow on state variables lets the shape change with regime. Because the density is exact rather than approximate, it can feed the tail metrics computed in Monte Carlo and bootstrap backtests and the CVaR construction in the HRP/CVaR portfolio pipeline; joint tail structure is handled separately in copula models for joint risk.

Conditional Density Estimation

The useful framing is a direct contrast between three published approaches to the same problem. TFT gives you a fixed set of quantiles from a quantile output layer. Conformal prediction gives you calibrated intervals with a coverage guarantee. A conditional CNF gives you an exact, differentiable density:

$\frac{dz(t)}{dt} = f(z(t), t, \text{features}_t, \theta)$

Differentiability is the distinguishing property: the density can sit inside a downstream objective and be backpropagated through, which neither fixed quantiles nor conformal intervals support. Whether the exact density is worth its cost against TFT quantiles on the same target is untested here.

Comparison with Discrete Alternatives

The LSTM-versus-attention half of this table is argued in detail, with benchmarks, in the TFT article — see its "TFT vs LSTM vs Vanilla Transformer" and "When LSTM Still Wins" sections. The columns that are new here are the two on the right, and the rows that decide the H2/H3 questions are the last two.

Python Implementation with torchdiffeq

torchdiffeq provides ODE solvers with adjoint backpropagation.

Basic Neural ODE for Price Dynamics

import torch
import torch.nn as nn
from torchdiffeq import odeint_adjoint as odeint

class PriceDynamics(nn.Module):
    """Neural network defining dh/dt = f(h, t)."""

    def __init__(self, hidden_dim: int = 64):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(hidden_dim, 128),
            nn.Tanh(),
            nn.Linear(128, 128),
            nn.Tanh(),
            nn.Linear(128, hidden_dim),
        )

    def forward(self, t, h):
        return self.net(h)


class NeuralODEPredictor(nn.Module):
    """
    Encode observed features -> latent state,
    evolve via Neural ODE,
    decode to price prediction.
    """

    def __init__(self, input_dim: int, hidden_dim: int = 64):
        super().__init__()
        self.encoder = nn.Linear(input_dim, hidden_dim)
        self.dynamics = PriceDynamics(hidden_dim)
        self.decoder = nn.Linear(hidden_dim, 1)

    def forward(self, x0, eval_times):
        """
        x0:         (batch, input_dim) features at t=0
        eval_times: (T,) times at which to evaluate the ODE
        Returns:    (T, batch, 1) predictions
        """
        h0 = self.encoder(x0)                          # (batch, hidden_dim)
        h_traj = odeint(self.dynamics, h0, eval_times,
                        method='dopri5', rtol=1e-5, atol=1e-7)
        return self.decoder(h_traj)

ODE-RNN for Irregular Tick Data

This is the experimental arm. The baseline it must beat is nn.GRUCell on the same stream with log(t_next - t_prev) concatenated to x.

class ODERNNCell(nn.Module):
    """Single step: ODE-evolve, then RNN-update."""

    def __init__(self, input_dim: int, hidden_dim: int = 64):
        super().__init__()
        self.dynamics = PriceDynamics(hidden_dim)
        self.gru_cell = nn.GRUCell(input_dim, hidden_dim)

    def forward(self, h, x, t_prev, t_next):
        times = torch.tensor([t_prev, t_next], dtype=torch.float32)
        h_evolved = odeint(self.dynamics, h, times,
                           method='dopri5')[-1]  # state at t_next
        h_updated = self.gru_cell(x, h_evolved)
        return h_updated


class ODERNN(nn.Module):
    """Process irregularly-sampled sequence."""

    def __init__(self, input_dim: int, hidden_dim: int = 64):
        super().__init__()
        self.cell = ODERNNCell(input_dim, hidden_dim)
        self.decoder = nn.Linear(hidden_dim, 1)
        self.hidden_dim = hidden_dim

    def forward(self, observations, times):
        """
        observations: list of (batch, input_dim) tensors
        times:        list of floats, observation timestamps
        """
        batch_size = observations[0].shape[0]
        h = torch.zeros(batch_size, self.hidden_dim)

        outputs = []
        for i in range(len(observations)):
            t_prev = 0.0 if i == 0 else times[i - 1]
            h = self.cell(h, observations[i], t_prev, times[i])
            outputs.append(self.decoder(h))

        return torch.stack(outputs)  # (seq_len, batch, 1)

Training Loop

def train_neural_ode(model, train_loader, epochs=100, lr=1e-3):
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(
        optimizer, T_max=epochs
    )

    for epoch in range(epochs):
        epoch_loss = 0.0
        for batch in train_loader:
            features, times, targets = batch
            optimizer.zero_grad()

            predictions = model(features, times)
            loss = torch.nn.functional.mse_loss(predictions, targets)

            loss.backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
            optimizer.step()

            epoch_loss += loss.item()

        scheduler.step()

        if (epoch + 1) % 10 == 0:
            avg_loss = epoch_loss / len(train_loader)
            print(f"Epoch {epoch+1}/{epochs}, Loss: {avg_loss:.6f}")

Continuous Normalizing Flow for Return Distributions

from torchdiffeq import odeint

class CNFDynamics(nn.Module):
    """Dynamics for continuous normalizing flow."""

    def __init__(self, dim: int = 1, hidden_dim: int = 64):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(dim + 1, hidden_dim),  # +1 for time
            nn.Tanh(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, dim),
        )
        self.dim = dim

    def forward(self, t, state):
        z = state[..., :self.dim]
        t_expand = t.expand(z.shape[0], 1)
        zt = torch.cat([z, t_expand], dim=-1)
        dz = self.net(zt)

        e = torch.randn_like(z)
        e_dz = torch.autograd.grad(
            dz, z, e, create_graph=True
        )[0]
        trace_jac = (e_dz * e).sum(dim=-1, keepdim=True)

        return torch.cat([dz, -trace_jac], dim=-1)


class ReturnDistributionCNF(nn.Module):
    """Model return distributions with continuous normalizing flows."""

    def __init__(self, dim: int = 1):
        super().__init__()
        self.dynamics = CNFDynamics(dim)
        self.dim = dim

    def log_prob(self, x):
        """Compute log probability of observed returns."""
        log_p0 = torch.zeros(x.shape[0], 1)
        state0 = torch.cat([x, log_p0], dim=-1)
        state0.requires_grad_(True)

        times = torch.tensor([1.0, 0.0])  # backward
        state_T = odeint(self.dynamics, state0, times,
                         method='dopri5')[-1]

        z_T = state_T[..., :self.dim]
        delta_log_p = state_T[..., self.dim:]

        log_p_base = -0.5 * (z_T ** 2 + torch.log(
            torch.tensor(2 * torch.pi)
        )).sum(dim=-1, keepdim=True)

        return log_p_base + delta_log_p

    def sample(self, n_samples: int):
        """Generate samples from learned distribution."""
        z0 = torch.randn(n_samples, self.dim)
        log_p0 = torch.zeros(n_samples, 1)
        state0 = torch.cat([z0, log_p0], dim=-1)

        times = torch.tensor([0.0, 1.0])  # forward
        state_T = odeint(self.dynamics, state0, times,
                         method='dopri5')[-1]

        return state_T[..., :self.dim]

Practical Notes

These are the things that will bite you while running the experiment above.

Solver Choice and Speed

dopri5 gives accuracy guarantees but variable compute cost, which is why H3 is a separate hypothesis from H2: an edge that only exists under an adaptive solver may not survive a live latency budget. Fixed-step solvers (Euler, RK4) give predictable latency at the cost of accuracy. Practical compromise: train with dopri5, deploy with a fixed-step solver calibrated to match the adaptive solver's output on representative data — and measure the gap rather than assuming it is small.

odeint(f, h0, t, method='dopri5', rtol=1e-6, atol=1e-8)

odeint(f, h0, t, method='euler', options={'step_size': 0.1})

For stiff problems near discontinuous jumps, method='implicit_adams' or method='scipy_solver'.

Numerical Stability

If $f_\theta$ outputs large values the state diverges. Mitigations:

• Spectral normalization on layers of $f_\theta$ to control the Lipschitz constant.
• Gradient clipping during training (shown in the training loop above).
• Time normalization: scale timestamps to $[0, 1]$.
• Regularization on dynamics norm: add $\lambda \|f_\theta(h, t)\|^2$ to the loss.

Integration Interval

For data spanning months, do not integrate from $t = 0$ to $t = 10{,}000$ minutes:

t_normalized = (timestamps - timestamps[0]) / (timestamps[-1] - timestamps[0])

This keeps the solver in a numerically friendly regime and makes the learned dynamics scale-invariant. It also matters for the comparison: an unnormalized ODE-RNN can lose to the GRU baseline for purely numerical reasons, which would be a measurement artifact rather than a finding.

Handling Multiple Time Scales

Markets have dynamics at microsecond, second, minute and daily scales at once, and a single Neural ODE may struggle to hold all of them. Options: stack multiple ODE blocks at different time scales; expand the hidden dimension for both fast and slow capacity; or run separate Neural ODEs per frequency band and fuse outputs. (This is a model-capacity question, distinct from the backtest-fidelity question in adaptive resolution drill-down.)

Open Research Directions

Neural Jump SDEs: add a learned jump component for sudden dislocations (flash crashes, earnings surprises):

$dh = f \, dt + g \, dW + \int_{\mathbb{R}} \gamma(h, z) \, \tilde{N}(dt, dz)$

where $\tilde{N}$ is a compensated Poisson random measure and $\gamma$ a learned jump kernel.

Neural Controlled Differential Equations (Neural CDEs): replace the Wiener process with a general driving signal, letting observed data streams drive the dynamics. Natural for order flow, where the stream of trades and quotes drives latent market state.

Differentiable Market Simulation: use a Neural SDE as the generative model inside a differentiable simulator and train strategies end-to-end by backpropagating through it via the adjoint. Strategy and market model co-evolve.

Physics-Informed Neural ODEs: PINNs embed a differential-equation residual in the loss — the technique itself is introduced in the Navier-Stokes article. The financial application is to impose no-arbitrage, put-call parity and martingale conditions as penalty terms or hard constraints on the learned dynamics.

Conclusion

The continuous-time framing is structurally correct: markets are continuous processes observed at discrete, irregular times, and models should respect that. The toolchain is mature — torchdiffeq, torchsde, plain PyTorch for the rest — and the known costs are solver speed and numerical stability over long integration intervals.

What is not established is the part that decides whether any of this belongs in a production stack. Structural correctness is not evidence of predictive advantage, and the specific advantage claimed for continuous dynamics over a $\log(\Delta t)$ feature has not been measured on real trade data here. Until the H2 table above is filled in, treat everything below the fold as a well-specified hypothesis with a working implementation attached, not as a result.

References

• Chen, R.T.Q., Rubanova, Y., Bettencourt, J., Duvenaud, D. (2018). Neural Ordinary Differential Equations. NeurIPS 2018. arXiv:1806.07366
• Rubanova, Y., Chen, R.T.Q., Duvenaud, D. (2019). Latent ODEs for Irregularly-Sampled Time Series. NeurIPS 2019. arXiv:1907.03907
• Jia, J., Benson, A.R. (2019). Neural Jump Stochastic Differential Equations. NeurIPS 2019.
• Kidger, P., Morrill, J., Foster, J., Lyons, T. (2020). Neural Controlled Differential Equations for Irregular Time Series. NeurIPS 2020.
• Hasan, A., Pereira, J.M., Farsiu, S., Carin, L. (2021). Neural Network Stochastic Differential Equation Models with Applications to Financial Data Forecasting. arXiv:2111.13164
• torchdiffeq: github.com/rtqichen/torchdiffeq
• UvA Deep Learning Tutorials — Neural ODEs: uvadlc-notebooks.readthedocs.io