Conceptio › Archive › arXiv CS

arXiv CSopen access

Variational Inference for Lévy Process-Driven SDEs via Neural Tilting

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Variational Inference for Lévy Process-Driven SDEs via Neural Tilting

Yaman Kındap1 , Manfred Opper2 , Benjamin Dupuis3 , Umut Şimşekli3 , Tolga Birdal1

arXiv:2605.10934v1 [cs.LG] 11 May 2026

1 Imperial College London, UK 2 Technical University of Berlin, Germany INRIA, CNRS, Département d’Informatique de l’Ecole Normale Supérieure / PSL, France

Abstract Modelling extreme events and heavy-tailed phenomena is central to building reliable predictive systems in domains such as finance, climate science, and safetycritical AI. While Lévy processes provide a natural mathematical framework for capturing jumps and heavy tails, Bayesian inference for Lévy-driven stochastic differential equations (SDEs) remains intractable with existing methods: Monte Carlo approaches are rigorous but lack scalability, whereas neural variational inference methods are efficient but rely on Gaussian assumptions that fail to capture discontinuities. We address this tension by introducing a neural exponential tilting framework for variational inference in Lévy-driven SDEs. Our approach constructs a flexible variational family by exponentially reweighting the Lévy measure using neural networks. This parametrization preserves the jump structure of the underlying process while remaining computationally tractable. To enable efficient inference, we develop a quadratic neural parametrization that yields closed-form normalization of the tilted measure, a conditional Gaussian representation for stable processes that facilitates simulation, and symmetry-aware Monte Carlo estimators for scalable optimization. Empirically, we demonstrate that the method accurately captures jump dynamics and yields reliable posterior inference in regimes where Gaussian-based variational approaches fail, on both synthetic and real-world datasets.

Introduction

Real-world stochastic systems, from financial markets to climate dynamics and safety-critical AI, generate observations exhibit sudden discontinuities, asymmetric shocks, and extreme events whose probability decays as a power law x−α rather than exponentially. This heavy-tailed structure is a defining property of many consequential phenomena, not an anomalous departure from Gaussianity. Models that approximate it with Gaussian noise are therefore systematically miscalibrated at the extremes, where the cost of error is highest. The natural mathematical language for this full class of stochastic behaviour is the theory of infinitely divisible distributions and associated Lévy processes, which form the complete family of processes with independent and stationary increments [19, 27]. Providing tractable approximate Bayesian inference over latent trajectories in Lévy processdriven stochastic differential equations (SDEs) would therefore yield a principled, uncertainty-aware framework for modelling temporal dynamics across any domain where the data-generating process departs from Gaussianity. Neural SDE and ODE approaches have transformed continuous-time sequence modelling by providing end-to-end differentiable, scalable frameworks competitive on standard benchmarks [9, 55, 28, 39, 12], with recent extensions such as Neural Jump SDEs [25] and Neural MJD [18] explicitly incorporating discontinuities and achieving strong results on time series with abrupt changes. Deep time-series forecasters, including DeepAR, DLinear, and N-HiTS [49, 60, 8], achieve state-of-the-art performance on standard benchmarks through flexible, data-driven parametrizations without any differential Preprint.

Figure 1: Posterior sample paths for two representative realisations of the OU (left pair) and doublewell (right pair) systems, comparing the Gaussian SDE and our tilted-stable model. equation structure. Yet all of these methods share a fundamental limitation. They model noise, including jump magnitudes, with light-tailed distributions and fit their parameters by maximum likelihood. This commits each model to a parametric noise structure calibrated to average behaviour, yielding predictions that are systematically overconfident at extreme quantiles, the failure mode that heavy-tailed phenomena directly expose, as shown in Fig. 1. This leaves a critical gap: no tractable framework exists for posterior inference over the latent dynamics of a Lévy-driven SDE, leaving jump structure as a fixed parametric assumption rather than a latent quantity to be inferred from data. Existing variational inference (VI) methods for stochastic processes obtain tractable path-space objectives through Brownian change-of-measure formulas, but these cannot produce posterior inference over the Lévy measure. Continuous variational processes are singular with respect to jump processes, and Brownian drift corrections in jump-diffusions leave the jump law unchanged. We close this gap by learning a neural exponential tilt of the prior Lévy measure, preserving the end-to-end differentiability and scalability of neural methods while yielding a tractable variational posterior with explicit jump and drift contributions to the ELBO. Our contributions are: • A path-space VI framework for Lévy-driven SDEs that derives the optimal Markov variational posterior as an exponential tilt of the prior Lévy measure. • Neural Tilting: a quadratic neural parametrisation of the tilting function that preserves flexible stateand time-dependent posterior jump behaviour while yielding closed-form normalising constants. • Efficient simulation and ELBO optimisation algorithms that exploit the conditionally Gaussian representation of stable processes and the symmetry of their Lévy measures. • Empirical validation on synthetic and real-world datasets showing that the framework recovers heavy-tailed jump dynamics and improves tail calibration where Gaussian baselines fail. Our source code will be made publicly available under: circle-group.github.io/research/NeuralTilting.

Related work

Simulation and data-driven methods for Lévy systems. Numerical methods for Lévy-driven SDEs build on time-discretisation schemes [42, 17], with variance reduction via multilevel Monte Carlo [15, 24], adaptive grids [31], and exact simulation [41]. Continuous-time formulations [7, 30] and conditionally Gaussian jump representations [21, 29] improve efficiency; the latter technique is extended in our work from the prior Lévy measure to the exponentially tilted posterior. Data-driven approaches recover Lévy dynamics from observations via stable distribution properties [34] and Koopman operator methods [35]. Despite exactness or arbitrarily small discretisation error, all of these methods lack end-to-end differentiability, precluding integration with modern ML pipelines. Neural differential equation models. Neural ODEs [9] and their stochastic extensions provide scalable, end-to-end differentiable continuous-time sequence models, with scalable adjoint gradients [33] and expressive latent architectures [28]. Extensions incorporate fractional white noise [54] and rough-path-driven jump dynamics [22], and training efficiency has been substantially improved via stochastic optimal control [12]. All of these methods are restricted to continuous-path dynamics, precluding heavy-tailed jump behaviour. 2

Jump-diffusion and Markov jump process models. A complementary line introduces jump structure into neural dynamical models: [25] extends neural ODEs with event-driven latent discontinuities, [18] combines a neural Itô diffusion with a compound Poisson process, and [61] reformulates temporal point processes as neural jump-diffusion SDEs. All restrict jump magnitudes to light-tailed distributions, with no capacity to represent heavy-tailed Lévy geometry. A parallel line targets discrete-state Markov jump processes (MJPs) via variational mean-field approximations [40], neural ODE rate functions [50], expectation propagation [16], and zero-shot foundation models [4]; these address discrete state spaces rather than continuous heavy-tailed dynamics. Deep time-series forecasters. Recurrent networks with parametric likelihoods [49], hierarchical interpolation [8], linear models [60], and normalising flows conditioned on sequential context [44] achieve strong benchmark performance. Score-based diffusion models [43, 53, 37, 38] and interpretable diffusion architectures [58] provide richer generative distributions. Distribution-free calibration via conformal prediction [52, 56, 59] provides coverage guarantees independently of the noise model. All of these methods are restricted to light-tailed parametric noise distributions, with no capacity to represent heavy-tailed jump geometry. Variational inference (VI) for stochastic processes. VI for stochastic processes typically relies on the standard Girsanov theorem, which provides a tractable Radon–Nikodym derivative between path measures by changing the drift under a Brownian SDE, yielding a quadratic ELBO in closed form for Brownian motion [39] and its fractional counterpart [11], with black-box variants using neural approximate posteriors [48, 55]. For Lévy-driven SDEs, Fokker–Planck-based variational methods have been used to estimate drift functions under specified α-stable noise [13], but their objective is drift recovery rather than posterior inference over latent paths or the Lévy measure. The resulting limitation is structural: existing stochastic process VI methods provide no variational mechanism for learning posterior heavy-tailed jump behaviour.

Background

Definition 3.1 (Lévy process [1]). Let (Ω, F, P) be a probability space with filtration {Ft }t≥0 . An adapted stochastic process {Lt : t ≥ 0} with values in Rd is a Lévy process if it satisfies: (i) L0 = 0 a.s., (ii) independent increments, (iii) stationary increments, and (iv) stochastic continuity. The stationarity and independence of increments suggest that Lévy processes may be interpreted as a family of random walks in continuous time [5, 27, 10]. Remark 3.2 (Connection to infinitely divisible distributions). At any fixed time, the class of Lévy processes corresponds one-to-one with infinitely divisible distributions (IDDs), which are probability measures µ on Rd that can be written as the n-fold convolution µ = µn ∗ · · · ∗ µn for any positive integer n. This correspondence provides access to a rich family including Gaussian (smooth diffusion), Poisson (discrete events), and stable processes (self-similar, heavy-tailed). Brownian motion is the only Lévy process with continuous sample paths; all other settings yield jump-diffusion processes [27]. Definition 3.3 (Characteristic function [26]). A Lévy process {Lt }t≥0 in R without a Brownian component is defined through its characteristic function as Z iuLt iuy E e = exp t e − 1 − iy 1|y|<1 ν(dy) , (1) R0

R where ν is a Lévy measure on R0 := R \ {0} satisfying R0 min(1, y 2 ) ν(dy) < ∞. This representation shows that the distribution of a Lévy process is uniquely determined by its Lévy measure ν, so learning ν is sufficient to characterise the process. R Remark 3.4 (Truncated Lévy measures). We focus on infinite-activity processes with R ν(dy) = ∞. In practice, we work with a truncated measure restricted to Rτ := R \ (−τ, τ ) for a threshold τ > 0. This removes the infinite activity from small jumps while preserving the large jumps that characterise rare discontinuities. The contribution from small jumps below τ can be well-approximated by Gaussian noise, providing a justification for this choice in practical settings [3, 20]. Definition 3.5 (Symmetric α-stable Lévy process). The most prominent example of a heavy-tailed Lévy process is the symmetric α-stable process {Lα t } with Lévy measure ν(dy) = |y|−1−α dy,

α ∈ (0, 2).

(2)

The parameter α controls tail heaviness; stable processes have infinite variance for all α and a finite mean only for α ∈ (1, 2). The truncated version is ντ (dy) = |y|−1−α 1|y|≥τ dy, which is finite since 3

R∞

y −1−α dy = α−1 τ −α < ∞. Although ντ is a finite measure, the p-th moment of the jump size distribution remains infinite for p ≥ α, reflecting preservation of the power-law tail under truncation. τ

Variational Inference for Lévy processes (via Neural Tilting)

We consider stochastic processes X0:T := {Xt : 0 ≤ t ≤ T } with Xt ∈ Rd and fixed initial state X0 = x0 , on Skorokhod spaces (the function space of right-continuous paths with left limits, appropriate for processes with jump discontinuities) driven by general Lévy processes: dXt = ftθ (Xt ) dt + dLt + σ(Xt ) dBt ,

(3)

where ftθ (·) is a neural drift function with parameters θ, σ(·) ∈ Rd×d is a diffusion matrix with D(x) = σ(x)σ(x)⊤ , {Lt } is a pure-jump Lévy process, and {Bt } is a standard Brownian motion. d For d > 1 we assume the driving Lévy noise has independent components, so the Lévy measure on R Nd factorises as ν(dy) = i=1 ν (i) (dyi );the fully correlated case is a direction for future work. Given noisy observations Yti ∼ N Xti , σε2 at times {ti }ni=1 , our goal is to approximate the posterior distribution p(X0:T | Y0:T ) and learn the parameters θ. We focus on the purely non-Gaussian case σ(·) = 0 in our experiments in order to isolate the contribution of the jump structure, though the framework extends to joint jump-diffusion dynamics.

Since the log-likelihood log p(Y0:T | θ) associated with noisy observations generated from Eq. (3) is intractable, we propose a variational approach based on the evidence lower bound (ELBO): ( " n # ) X θ arg sup EQ log pθ (Yti | Xti ) − DKL Q ∥ P , (4) Q

i=1

where P θ (X0:T ) is the prior path measure induced by Eq. (3) and the supremum is over trial path measures Q, equaling the log-likelihood log pθ (Y0:T | θ), attained when Q = pθ (X0:T | Y0:T ). 4.1 Variational optimisation over Markov path measures We derive the variational family by optimising over Markov trial path measures, following the Doob h-transform structure of conditioned Markov processes [45]. Theorem 4.1 (Optimal variational posterior family). Let P θ be the prior path measure of the Lévy– diffusion SDE in Eq. (3), with generator LP t in Eq. (17). Optimising the ELBO in Eq. (4) over Markov trial path measures yields the variational posterior SDE as dXt = ftθ (Xt ) + D∇ϕt (Xt ) dt + σ(Xt ) dBt + dLϕt , (5) where ϕt is the variational potential function, {Lϕt } is the tilted Lévy process1 with jump measure ν̃ such that ν̃(dy, t, Xt ) = eϕt (Xt +y)−ϕt (Xt ) ντ (dy). (6) and D = σσ ⊤ is the diffusion matrix. For smooth test functions G, the posterior generator is ⊤ θ 2 1 LQ t G(x) = ft (x) + D∇ϕt (x) ∇G(x) + 2 tr D ∇ G(x) Z + G(x + y) − G(x) eϕt (x+y)−ϕt (x) ντ (dy).

(7)

Compared to Eq. (17), the optimal variational family retains the same diffusion coefficient but acquires a Brownian drift correction D∇ϕt and a tilted Lévy measure. Setting ϕt ≡ 0 in Eq. (7) recovers the prior generator LP t of Eq. (17) (see Sec. A). For pure-jump priors (σ = 0), both the Brownian term and the drift correction D∇ϕt vanish, recovering the form used in Sec. 5. The result holds for any Lévy measure ντ . Proof. The detailed, step-by-step proof is provided in Sec. A. Corollary 4.2 (KL divergence between Q and P θ ). For the variational posterior process in Theorem 4.1, the path-space KL divergence is "Z Z # Z T 1 T θ ⊤ DKL Q ∥ P = EQ f (y, t, Xt ) ντ (dy) dt + ∇ϕt (Xt ) D ∇ϕt (Xt ) dt , (8) 2 0 0 Rd 1 {Lϕ } is not a Lévy process since its jump measure depends on both time and the current state X . t t

where f (y, t, Xt ) = Ht (Xt , y) ln Ht (Xt , y) − Ht (Xt , y) + 1, Ht (x, y) = eϕt (x+y)−ϕt (x) . Under the factorised Lévy measure of Eq. (3), the jump integral decomposes into d independent scalar integrals. For pure-jump priors (σ = 0, hence D = 0), the Brownian term vanishes. The framework extends naturally to priors with state-dependent jump sizes γ(Xt− )y; see Sec. E for details and a discussion of the simulation adaptation this requires. 4.2

Quadratic parametrization for tractability

Practical implementation requires addressing the computational challenges of (a) sampling from the tilted process Lϕt and (b) evaluating the intractable expectations in Eq. (8). We address both through a quadratic parametrization that retains neural flexibility while enabling analytical simplifications. For clarity, we specialise to d = 1; the multivariate extension under the independence assumption follows from Rem. 4.3. Consider the quadratic tilting function: 2

ϕt (x′ ) = At x′ + Bt x′ ,

(9)

where At and Bt are time-dependent scalar coefficients parametrized by neural networks. For a jump from state x to x + y, the tilting factor becomes: Ht (x, y) = eϕt (x+y)−ϕt (x) = exp At (2xy + y 2 ) + Bt y . (10) The quadratic structure concentrates the time dependence into the scalar coefficients At and Bt , evaluated once per time step, while the dependence on state x and jump size y takes a fixed quadratic form. This makes Ht (x, y) a Gaussian-type function of y for any fixed t and x, enabling analytical computation of the normalising constants required for simulation. Remark 4.3 (Multivariate extension). For a d-dimensional process, the quadratic parametrization extends under the assumption of independent state dimensions. The tilting function decomposes as Pd Pd (i) (i) (i) ϕt (x′ ) = i=1 ϕt (x′i ) = i=1 At (x′i )2 + Bt x′i , where each dimension has its own scalar (i) (i) ′ coefficients At and Bt . This is equivalent to the matrix form ϕt (x′ ) = x′⊤ At x′ + b⊤ t x with (1) (d) (1) (d) ⊤ diagonal At = diag(At , . . . , At ) and bt = (Bt , . . . , Bt ) . The tilting factor then factorises Qd (i) (i) across dimensions as Ht (x, y) = i=1 exp At (2xi yi + yi2 ) + Bt yi , so the KL integral and the normalising constants of Sec. 4.4 remain tractable dimension-by-dimension. Adaptive temporal encoding. The coefficients At and Bt are outputs of two MLPs whose shared input is a learned embedding of time. For a query time t, the embedding is e(t) =

N X i=1

wi (t) vi ,

e−s | t−τi | wi (t) = PN , −s | t−τj | j=1 e

(11)

de where {τi }N are learnable embedding vectors, i=1 ⊂ [0, T ] are learnable reference times, {vi } ⊂ R and s > 0 is a learnable sharpness (stored in log-space to ensure positivity). This is a Nadaraya– Watson estimator with a Laplacian kernel: the weights wi (t) decay exponentially with temporal distance from each reference point. Crucially, both the reference locations and the embeddings are optimised jointly with the rest of the model, so the τi migrate during training toward times at which the ELBO is sensitive to the tilting, in practice toward the locations of jumps, without any explicit supervision of jump timing. The coefficients are then At = −(amin + softplus(fA (e(t)))) < 0 and Bt = fB (e(t)), where fA and fB are MLPs and amin > 0 is a small fixed positive constant that enforces a strict lower bound on |At |; the reparametrisation guarantees strict negativity for all t, as required for integrability of the tilted measure (Rem. 4.5).

Adaptive tempering and finite moment properties. For symmetric α-stable priors, the constraint At < 0, enforced by the parametrisation, ensures the tilted Lévy measure is integrable, a necessary condition for the variational SDE to be well-posed. The resulting exponentially truncated power-law tail behaviour is analogous to tempered stable processes [47], though here the tempering adapts with time and state rather than being fixed, ensuring all moments of the tilted process are finite. This lighter-tailed character is a consequence of the quadratic parametrisation: for instance, ϕt for which ϕt (Xt + y) − ϕt (Xt ) remains bounded in y yields a posterior with the same tail index α as the prior. 5

Algorithm 1 KL approximation for tilted Algorithm 2 Tilted Lévy SDE simulation −1 symmetric Lévy measures Input: Initial state Xt0 ; time grid {tj }N j=0 ; neural network (m)

parameters

Input: Process paths {Xtj }; neural outputs Atj , Btj 1: for each time step tj and each path m do 2: Evaluate Atj , Btj from the neural networks. Draw K samples {yk }K 3: k=1 via yk = τ (1 − uk )−1/α , uk ∼ Uniform(0, 1). (m) 4: Compute fk+ = f (yk , tj , Xtj ) and

1: for each time step tj do 2: Evaluate neural outputs Atj , Btj .

(m)

fk− = f (−yk , tj , Xtj ). Approximate: PK + − k=1 (fk + fk ). 6: end for

(m) I(tj , Xtj )

Draw jump y (i) ∼ N µy (r(i) ), σy2 (r(i) ) . 8: end for PNj (i) 9: Update: Xtj+1 = Xtj + ftθj (Xtj ) ∆t + i=1 y . 10: end for 7:

≈

C K

4.3

Draw K samples {yk } via yk = τ (1 − uk )−1/α , uk ∼ Uniform(0, 1). 4: Compute Λj via Eq. (15); draw Nj ∼ Poisson(Λj ). 5: for i = 1, . . . , Nj do 6: Propose r∗ ∼ r−1−α ; accept as r(i) with probability C(r∗ , t, Xt )/M (t, Xt ). 3:

Approximation of the KL divergence (m)

With M Monte Carlo (MC) paths {Xtj } drawn from the variational posterior (Sec. 4.4), the ELBO Eq. (4) is approximated as " n # N −1 M X 1 X X (m) (m) ∆tj Î tj , Xtj log pθ Yti | Xti − , (12) L(θ, ϕ) ≈ M m=1 i=1 j=0 R where Î(tj , Xtj ) is a Monte Carlo estimate of I(tj , Xtj ) := R f (y, t, Xt ) ντ (dy) at time tj . The remainder of this subsection develops Î; path simulation is detailed in Sec. 4.4. We exploit the measure to reduce the integration domain to I(t, Xt ) = R ∞ symmetry of the truncated stable −1−α f (y, t, X ) + f (−y, t, X ) y dy. We express this as an expectation over jump sizes t t τ y ∼ p(y) = ατ α y −1−α (supported on [τ, ∞)), with normalisation constant C = α−1 τ −α . Samples from p(y) are drawn via inverse-CDF: y = τ (1 − u)−1/α where u ∼ Uniform(0, 1). Alg. 1 summarizes the resulting MC procedure. The total computational cost is O(M N K), where M is the number of MC paths, N the number of time steps, and K the number of jump samples. 4.4

Forward simulation of tilted Lévy processes

The variational SDE Eq. (5) requires forward simulation of the tilted process Lϕt with state- and time-dependent jump measure ν̃. We exploit the conditionally Gaussian structure of stable processes. Theorem 4.4 (Conditionally Gaussian representation of tilted R ∞ stable jumps). The tilted jump measure ν̃(y, t, Xt ) admits a disintegration ν̃(y, t, Xt ) = 0 σ̃(r; y, t, Xt ) π̃τ (r, t, Xt ) dr, where σ̃(r; ·, t, Xt ) is a probability kernel and π̃τ is a truncated non-negative mixing measure: 2 eϕt (Xt +y)−ϕt (Xt ) N y, r2 σG C(r, t, Xt ) r−1−α σ̃(r; y, t, Xt ) = , π̃τ (r, t, Xt ) = 1r≥τ , (13) C(r, t, Xt ) K(α, t, Xt ) with normalising constants C(r, t, Xt ) and K(α, t, Xt ), and σG > 0 the Gaussian scale of the variance-mean mixture representation of the stable process. Proof. The result extends the disintegration of α-stable jump measures established by [46, 32] to include the tilting factor eϕt (Xt +y)−ϕt (Xt ) . See Sec. B for a detailed proof. Remark 4.5 (Analytical tractability under quadratic parametrization). Under Eq. (9), the normalising constant C(r, t, Xt ) of the conditional kernel is available in closed form. Setting K1 = 2At Xt + Bt and K2 = At − 2r21σ2 (K2 < 0 is enforced by At < 0, as guaranteed by the parametrisation of G K2 Sec. 4.2), C(r, t, Xt ) = √ 1 2 2 exp − 4K12 . Thus, the conditional kernel σ̃(r; ·, t, Xt ) reduces −2K2 r σG to N µy (r), σy2 (r) with µy (r) = −K1 /(2K2 ) and σy2 (r) = −1/(2K2 ).

Table 1: Synthetic experiment results (σε = 0.10), reporting means over 50 realisations per α. Held-out CRPS, mean absolute parameter recovery errors, and jump CRPS at thresholds (p97.5 , p99 ) for the OU and double-well (DW) systems; lower is better. Best result per metric in bold. Ornstein–Uhlenbeck CRPS ↓

|θ̂−θ | ↓

|µ̂−µ | ↓

p97.5 ↓

p99 ↓

CRPS ↓

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

2.88 0.89 4.63 2.27 2.97 5.43 2.15 2.10 2.46 2.41 0.41 0.18 1.46 0.76 1.06 0.20 0.85 0.60 1.35 1.05 1.63 0.63 4.32 1.98 3.01 6.11 1.12 0.56 1.23 0.67 0.34 0.14 1.81 0.79 1.41 0.33 0.58 0.37 0.85 0.62 0.48 0.37 5.25 2.43 2.57 4.54 2.35 2.36 3.01 3.20 0.37 0.16 1.95 1.05 1.72 0.33 0.55 0.35 0.77 0.54 0.72 0.61 6.01 2.63 2.66 4.84 0.82 0.79 1.25 1.26 0.36 0.15 1.72 1.24 1.62 0.28 0.60 0.38 0.85 0.60 0.62 0.45 6.11 1.96 2.62 3.51 0.86 0.71 1.14 0.98 0.34 0.14 1.15 1.39 1.07 0.30 0.47 0.29 0.61 0.41 0.41 0.23 4.43 2.04 2.57 3.56 0.44 0.33 0.52 0.38 0.36 0.14 1.42 1.16 1.21 0.27 0.47 0.27 0.59 0.38 0.42 0.28 5.79 2.38 2.29 4.41 0.52 0.43 0.58 0.48 0.37 0.15 1.08 1.26 0.76 0.25 0.47 0.28 0.56 0.38 0.36 0.19 5.23 2.16 1.72 2.20 0.91 0.28 0.98 0.34 0.37 0.15 1.11 1.19 1.03 0.30 0.47 0.26 0.55 0.34 0.39 0.22 4.99 1.62 2.04 3.36 0.51 0.31 0.54 0.33 0.38 0.16 2.01 1.14 1.49 0.39 0.46 0.26 0.52 0.32

∗

Double well |θ̂1 −θ1∗ | ↓

∗

|θ̂2 −θ2∗ | ↓

p97.5 ↓

p99 ↓

All α 0.88 0.43 5.19 2.16 2.50 4.21 1.08 0.87 1.30 1.12 0.37 0.15 1.52 1.11 1.26 0.29 0.55 0.34 0.74 0.52

Theorem 4.6 (Exact rejection sampler for the tilted mixing measure). Under the quadratic parametrisation of Sec. 4.2, the tilted mixing measure π̃τ (r, t, Xt ) ∝ C(r, t, Xt ) r−1−α of Theorem 4.4 admits an exact rejection sampler with proposal q(r) ∝ r−1−α , r ≥ τ . The normalising constant C(r, t, Xt ) of Rem. 4.5 satisfies the state-dependent envelope C(r, t, Xt ) ≤ M (t, Xt ) = exp K12 /4|At | , (14) so each proposal r∗ is accepted with probability C(r∗ , t, Xt )/M (t, Xt ). The sampler is exact and requires no discretisation, gradient computation, or step-size tuning. Proof. See Sec. C. Jump intensity estimation. For time interval (tj , tj+1 ) of length ∆t, the total jump intensity is approximated as K ∆t C X Λj ≈ Ht (Xtj , yk ) + Ht (Xtj , −yk ) , (15) K k=1

reusing the samples {yk } from the KL computation. The number of jumps in the interval is then Nj ∼ Poisson(Λj ). Alg. 2 assembles the complete forward simulation procedure. The total computational cost of Alg. 2 is O(M N (K + N̄j )) per iteration, where N̄j is the expected number of jumps per time interval; the rejection sampler adds a constant expected factor per jump equal to the reciprocal of the mean acceptance probability. In training, the K intensity samples are shared with Alg. 1, reducing the incremental simulation cost to O(M N N̄j ).

Experiments

We validate our framework on two settings that both exhibit pronounced heavy-tailed dynamics: (i) synthetic data with ground-truth parameters and (ii) challenging real-world forecasting task. Datasets. In the synthetic setting, data are generated from stable-process-driven SDEs with known parameters and trajectories, enabling direct assessment of posterior quality and parameter recovery. We consider two drift families: a linear Ornstein–Uhlenbeck (OU) system (parameters θ, µ) and a nonlinear double-well potential (parameters θ1 , θ2 ).2 For each stability index α ∈ {1.1, 1.2, . . . , 1.9} we generate 50 independent realisations with randomly drawn drift parameters. In the financial setting, we work with hourly log-prices of ten technology stocks (NVDA, GOOGL, MSFT, AAPL, AMZN, META, TSLA, AMD, NFLX, INTC) spanning June 2024 to April 2026. Models are trained on rolling 30-day windows and evaluated on the subsequent 2-day forecast horizon, yielding between 302 and 314 non-overlapping evaluation periods depending on ticker availability. We report results on both univariate (per-ticker) and multivariate (d = 10) forecasting tasks. Baselines. For the synthetic experiments we compare against a Gaussian SDE with an identical drift parametrisation, isolating the effect of the noise model. For the financial experiments we additionally compare against DeepAR [49], N-HiTS [8], DLinear [60], Neural Jump SDEs [25], Neural MJD [18], and a Gaussian SDE. 2 OU drift: f θ (x) = θ(µ − x) with θ > 0. Double-well drift: f θ (x) = θ x − θ x3 with θ , θ > 0. 1 2 1 2

Figure 2: Reliability diagrams. Empirical coverage vs. nominal level; the diagonal denotes perfect calibration. TS tracks the diagonal closely across all settings, while all probabilistic baselines severely under-cover. Left: NVDA univariate (314 evaluation windows). Centre: GOOGL univariate (304 windows). Right: Multivariate (d=10, 302 windows). Metrics. Our primary metric is the Continuous Ranked Probability Score (CRPS) [36], which jointly rewards sharpness and calibration and applies to any model that generates sample paths. For synthetic experiments we additionally report mean absolute parameter recovery error |θ̂ − θ∗ |. For financial experiments we supplement CRPS with MSE and MAE to enable comparison with deterministic baselines, and report a jump CRPS computed on the subset of price increments exceeding the p-th percentile for p ∈ {90, 95, 97.5, 99} to specifically assess tail performance. Uncertainty calibration is assessed via reliability diagrams. Full implementation details, including wall-clock training times (Tab. 3), are provided in Sec. G. 5.1

Synthetic experiments

Linear drift (Ornstein–Uhlenbeck) Tab. 1 reports held-out CRPS, parameter recovery errors at σε = 0.10 and jump CRPS; additional details and results at σε = 0.05 are in Sec. H. Our model achieves lower CRPS than the Gaussian baseline at every α value, with the improvement most pronounced for small α, reaching 3.2× at α = 1.1, and narrowing as α → 2 as expected. For the mean-reversion rate θ, our model substantially reduces recovery error (2.16 vs. 5.19 averaged over all α), confirming that correctly modelling heavy-tailed noise prevents the drift from absorbing spurious jump contributions; the global mean µ is harder to recover due to partial trade-offs with the tilting function when observations are sparse. Non-linear drift (double well). The double-well system shows consistent CRPS improvement across all α (0.15 vs. 0.37, a 2.4× ratio). The most striking result is for the cubic parameter θ2 , whose recovery error falls from 1.26 to 0.29, a 4.3× improvement, confirming that without a correct heavy-tailed noise model the drift absorbs rare large-amplitude excursions as apparent curvature. Recovery of the linear parameter θ1 is also improved overall (1.11 vs. 1.52), with the tail advantage confirmed at every threshold in Tab. 4. Results at σε = 0.05 in Sec. H show the same pattern. 5.2

Financial forecasting

We evaluate on NVDA hourly log-prices across 314 rolling 30-day training windows with a 2day forecast horizon. Each window proceeds in two phases. In the training phase, the ELBO is maximised over θ and ϕ, serving as a system identification procedure, with α fixed by grid search over {1.1, . . . , 1.9}. Forecasts are then generated by the learned prior SDE (with learned θ and ϕ = 0), initialised from posterior samples at the window boundary. The Gaussian SDE follows the same protocol, isolating the noise model as the sole difference. Tab. 2 reports mean CRPS alongside jump CRPS at four tail thresholds. Results for GOOGL follow the same qualitative pattern, confirming that the findings generalise across stocks; full results are provided in Sec. I. Overall CRPS. Our model achieves the lowest mean CRPS (0.432), a 21% improvement over the next-best probabilistic baseline (Neural MJD, 0.545). The Gaussian SDE, which shares the same drift architecture and differs only in the noise model, scores 54% worse (0.663), isolating the benefit of heavy-tailed noise. The deterministic baselines N-HiTS and DLinear produce point forecasts so their CRPS reduces to MAE; they rank fifth and sixth. Tail performance. Tab. 2 shows that our advantage is concentrated in the tails. The gap over the next-best baseline (DeepAR) grows in absolute terms as the threshold rises, reaching 0.277 at p99 compared to 0.186 at p90 , consistent with the synthetic finding that a correctly specified heavy-tailed noise model specifically improves fit to rare large-amplitude moves. 8

Table 2: Financial forecasting results; lower is better for all metrics. Best in bold. † Deterministic model; CRPS equals MAE. Jump CRPS computed on price increments exceeding the indicated percentile threshold. Left: Univariate (NVDA), 314 evaluation windows. Right: Multivariate (d=10), 302 windows; NJ-SDE excluded (diverged in all runs). Jump CRPS ↓

Jump CRPS ↓ Model

CRPS ↓

p90

p95

p97.5

p99

NJ-SDE 0.849±1.125 1.385 1.615 1.769 1.973 DLinear† 0.802±0.708 1.357 1.600 1.637 2.086 N-HiTS† 0.746±0.635 1.406 1.594 1.690 2.108 Gaussian SDE 0.663±0.984 1.113 1.233 1.436 1.929 DeepAR 0.556±0.509 1.000 1.145 1.334 1.749 Neural MJD 0.545±0.530 1.059 1.300 1.495 2.021 TS (ours) 0.432±0.427 0.814 0.970 1.133 1.472

Model

CRPS ↓

Energy ↓

p90

p95

p97.5

p99

N-HiTS† 0.690±0.297 2.775 1.031 1.132 1.290 1.453 Gaussian SDE 0.556±0.241 2.237 0.766 0.827 0.926 1.093 DeepAR 0.549±0.302 2.222 0.874 0.952 1.095 1.261 DLinear† 0.526±0.260 2.137 0.850 0.952 1.135 1.288 Neural MJD 0.504±0.260 2.024 0.799 0.899 1.061 1.224 TS (ours) 0.499±0.247 2.020 0.725 0.788 0.886 1.015

Calibration. Our model produces well-calibrated predictive intervals: empirical coverage closely tracks the nominal level across the full range of prediction intervals, as illustrated by the reliability diagram in Fig. 2. All probabilistic baselines severely under-cover at every nominal level, indicating overconfident predictive distributions. 5.3 Multivariate financial forecasting We evaluate on all ten stocks jointly (d = 10) across 302 rolling evaluation windows, using the energy score as the primary multivariate scoring rule alongside per-dimension CRPS. NJ-SDE diverged in every multivariate run and is excluded. Tab. 2 reports the results. Results. Our model achieves state-of-the-art average performance, matching Neural MJD on CRPS (0.499 vs 0.504) and energy score (2.020 vs 2.024), with the distinction emerging in tail performance. Tail performance. Tab. 2 shows our model leads at every jump CRPS threshold, and the absolute gap over the next-best model (Gaussian SDE) widens from 0.041 at p90 to 0.078 at p99 , consistent with the pattern observed in both the synthetic and univariate financial experiments. Neural MJD, despite being near-tied on overall CRPS, drops to third on jump CRPS at every threshold, suggesting its mixture approximation under-represents the most extreme joint moves. Calibration. Our model remains well-calibrated in the multivariate setting (50% PI → 45.8% empirical coverage; 90% PI → 82.7%). The Gaussian SDE is also reasonably calibrated (54.1%, 86.4%), while Neural MJD and DeepAR severely under-cover (9.6% and 14.6% at the 50% level respectively), indicating that their probabilistic forecasts are overconfident. Reliability diagrams for both GOOGL and the multivariate setting are shown in Fig. 2.

Conclusion

We introduced a variational inference framework for Lévy-driven SDEs in which the approximate posterior is expressed as a tilted Lévy process, with the prior jump measure reweighted by a learned exponential factor and the Brownian component acquiring the standard score correction. Although our experiments focus on stable processes, the change-of-measure construction applies to the broader class of Lévy processes, opening a route to scalable posterior inference in models with general non-Gaussian jump structure. The quadratic neural parametrisation is the computational device that makes this tractable, yielding closed-form normalising constants and an exact rejection sampler for the tilted jump law. Across all experimental settings the framework delivers 2–4× CRPS improvements over Gaussian baselines on synthetic systems, a 21% improvement over the next-best probabilistic model on financial data, and the strongest tail performance at every jump-CRPS threshold in both univariate and multivariate forecasting. Parameter recovery confirms that without a correctly specified heavytailed noise model, the drift absorbs rare large-amplitude excursions as spurious curvature, distorting inferred dynamics even when the overall fit appears reasonable. Limitations and future work. The quadratic parametrisation tempers posterior tails and may underrepresent asymmetric or multi-modal jump laws; richer flow-based tilting functions are a natural next step. Our derivation uses a time-discretised change-of-measure argument; integrability conditions for a fully formal treatment are discussed in Rem. A.1. The multivariate formulation currently assumes independent stable dimensions; incorporating Lévy copulas would allow cross-dimensional tail dependence to be modelled explicitly. The stability index α is fixed by grid search; joint estimation within the ELBO remains future work. 9

Acknowledgments and Disclosure of Funding T. B. was supported by a UKRI Future Leaders Fellowship (MR/Y018818/1). The authors acknowledge support from the UK AI Research Resource (AIRR Isambard AI) through grant 0251-45840945-1 - TopoFound. U.Ş. is partially supported by the French government under the management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). B.D. and U.Ş. are partially supported by the European Research Council Starting Grant DYNASTY – 101039676.

References [1] David Applebaum. Lévy processes and stochastic calculus. Cambridge university press, 2009. [2] Ran Aroussi. yfinance: Download market data from Yahoo! Finance’s API, 2019. [3] Søren Asmussen and Jan Rosiński. Approximations of small jumps of Lévy processes with a view towards simulation. Journal of Applied Probability, 38(2):482–493, 2001. [4] David Berghaus, Kostadin Cvejoski, Patrick Seifner, César Ojeda, and Ramsés J. Sánchez. Foundation inference models for Markov jump processes. In Advances in Neural Information Processing Systems, volume 37, 2024. [5] J. Bertoin. Lévy Processes. Cambridge Tracts in Mathematics. Cambridge University Press, 1996. [6] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Yash Katariya, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye WandermanMilne, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. [7] Peter J Brockwell. Lévy-driven CARMA processes. Annals of the Institute of Statistical Mathematics, 53(1):113–124, 2001. [8] Cristian Challu, Kin G. Olivares, Boris N. Oreshkin, Federico Garza, Max MergenthalerCanseco, and Artur Dubrawski. N-HiTS: Neural hierarchical interpolation for time series forecasting. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, 2023. [9] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018. [10] R Cont and P Tankov. Financial Modelling with Jump Processes. Chapman & Hall/CRC, 2003. [11] Rembert Daems, Manfred Opper, Guillaume Crevecoeur, and Tolga Birdal. Variational inference for sdes driven by fractional noise. In The Twelfth International Conference on Learning Representations (ICLR 2024), 2024. [12] Rembert Daems, Manfred Opper, Guillaume Crevecoeur, and Tolga Birdal. Efficient Training of Neural SDEs Using Stochastic Optimal Control. In European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN), 2025. [13] Min Dai, Jinqiao Duan, Jianyu Hu, and Xiangjun Wang. Variational inference of the drift function for stochastic differential equations driven by Lévy processes. Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(6):061103, 2022. [14] DeepMind, Igor Babuschkin, Kate Baumli, Alison Bell, Surya Bhupatiraju, Jake Bruce, Peter Buchlovsky, David Budden, Trevor Cai, Aidan Clark, Ivo Danihelka, Antoine Dedieu, Claudio Fantacci, Jonathan Godwin, Chris Jones, Ross Hemsley, Tom Hennigan, Matteo Hessel, Shaobo Hou, Steven Kapturowski, Thomas Keck, Iurii Kemaev, Michael King, Markus Kunesch, Lena Martens, Hamza Merzic, Vladimir Mikulik, Tamara Norman, George Papamakarios, John Quan, Roman Ring, Francisco Ruiz, Alvaro Sanchez, Laurent Sartran, Rosalia Schneider, Eren Sezener, Stephen Spencer, Srivatsan Srinivasan, Miloš Stanojević, Wojciech Stokowiec, Luyu Wang, Guangyao Zhou, and Fabio Viola. The DeepMind JAX Ecosystem, 2020. 10

[15] Steffen Dereich and Felix Heidenreich. A multilevel Monte Carlo algorithm for Lévy-driven stochastic differential equations. Stochastic Processes and their Applications, 121(7):1565– 1587, 2011. [16] Yannick Eich, Bastian Alt, and Heinz Koeppl. Entropic matching for expectation propagation of Markov jump processes. In Proceedings of the 28th International Conference on Artificial Intelligence and Statistics (AISTATS), 2025. [17] Nicolas Fournier. Simulation and approximation of Lévy-driven stochastic differential equations. ESAIM: Probability and Statistics, 15:233–248, 2011. [18] Yuanpei Gao, Qi Yan, Yan Leng, and Renjie Liao. Neural Non-Stationary Merton Jump Diffusion for Time Series Prediction. In Advances in Neural Information Processing Systems, volume 38, 2025. [19] Boris Vladimirovich Gnedenko and Andrey Nikolaevich Kolmogorov. Limit distributions for sums of independent random variables, volume 2420. Addison-wesley, 1968. [20] Simon Godsill, Ioannis Kontoyiannis, and Marcos Tapia Costa. Generalised shot-noise representations of stochastic systems driven by non-Gaussian Lévy processes. Advances in Applied Probability, 56(4):1215–1250, 2024. [21] Simon Godsill, Marina Riabiz, and Ioannis Kontoyiannis. The Lévy State Space Model. In 2019 53rd Asilomar Conference on Signals, Systems, and Computers, pages 487–494, 2019. [22] Christian Holberg and Cristopher Salvi. Exact gradients for stochastic spiking neural networks driven by rough signals. In Advances in Neural Information Processing Systems, volume 37, 2024. [23] Jean Jacod and Albert Shiryaev. Limit theorems for stochastic processes, volume 288. Springer Science & Business Media, 2013. [24] Ajay Jasra, Kody JH Law, and Prince Peprah Osei. Multilevel particle filters for Lévy-driven stochastic differential equations. Statistics and Computing, 29(4):775–789, 2019. [25] Junteng Jia and Austin R Benson. Neural jump stochastic differential equations. Advances in Neural Information Processing Systems, 32, 2019. [26] O. Kallenberg. Foundations of Modern Probability. Springer-Verlag, 2nd edition, 2002. [27] Sato Ken-Iti. Lévy Processes and Infinitely Divisible Distributions, volume 68. Cambridge University Press, 1999. [28] Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for irregular time series. Advances in neural information processing systems, 33:6696– 6707, 2020. [29] Yaman Kındap and Simon Godsill. Generalised hyperbolic state-space models for inference in dynamic systems. IEEE Open Journal of Signal Processing, 5:132–139, 2023. [30] David Kleinhans and Rudolf Friedrich. Continuous-time random walks: Simulation of continuous trajectories. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 76(6):061102, 2007. [31] Arturo Kohatsu-Higa and Peter Tankov. Jump-adapted discretization schemes for Lévy-driven SDEs. Stochastic Processes and their Applications, 120(11):2258–2285, 2010. [32] Tatjana Lemke, Marina Riabiz, and Simon J Godsill. Fully bayesian inference for α-stable distributions using a poisson series representation. Digital Signal Processing, 47:96–115, 2015. [33] Xuechen Li, Ting-Kam Leonard Wong, Ricky T. Q. Chen, and David Duvenaud. Scalable gradients for stochastic differential equations. In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS), 2020. 11

[34] Yang Li and Jinqiao Duan. Extracting stochastic dynamical systems with α-stable Lévy noise from data. Journal of Statistical Mechanics: Theory and Experiment, 2022(2):023405, 2022. [35] Yubin Lu and Jinqiao Duan. Discovering transition phenomena from data of stochastic dynamical systems with Lévy noise. Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(9):093110, 2020. [36] James E Matheson and Robert L Winkler. Scoring rules for continuous probability distributions. Management science, 22(10):1087–1096, 1976. [37] Gabriel Nobis, Maximilian Springenberg, Marco Aversa, Michael Detzel, Rembert Daems, Roderick Murray-Smith, Shinichi Nakajima, Sebastian Lapuschkin, Stefano Ermon, Tolga Birdal, et al. Generative fractional diffusion models. Advances in neural information processing systems, 37:25469–25509, 2024. [38] Gabriel Nobis, Maximilian Springenberg, Arina Belova, Rembert Daems, Christoph Knochenhauer, Manfred Opper, Tolga Birdal, and Wojciech Samek. Fractional diffusion bridge models. In Advances in neural information processing systems, 2025. [39] Manfred Opper. Variational inference for stochastic differential equations. Annalen der Physik, 531(3):1800233, 2019. [40] Manfred Opper, Andreas Ruttor, and Guido Sanguinetti. Approximate inference in continuous time Gaussian-jump processes. In Advances in Neural Information Processing Systems, volume 23, 2010. [41] Murray Pollock, Adam M Johansen, and Gareth O Roberts. On the exact and ε-strong simulation of (jump) diffusions. Bernoulli, 22(2):794–856, 2016. [42] Philip Protter and Denis Talay. The Euler scheme for Lévy driven stochastic differential equations. The Annals of Probability, 25(1):393–423, 1997. [43] Kashif Rasul, Calvin Seward, Ingmar Schuster, and Roland Vollgraf. Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting. In Proceedings of the 38th International Conference on Machine Learning (ICML), 2021. [44] Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann, and Roland Vollgraf. Multivariate probabilistic time series forecasting via conditioned normalizing flows. In The Ninth International Conference on Learning Representations (ICLR), 2021. [45] L Chris G Rogers and David Williams. Diffusions, Markov processes, and martingales, volume 2. Cambridge university press, 2000. [46] Jan Rosiński. Series representations of Lévy processes from the perspective of point processes. In Lévy processes: theory and applications, pages 401–415. Birkhäuser Boston Boston, MA, 2001. [47] Jan Rosiński. Tempering stable processes. Stochastic processes and their applications, 117(6):677–707, 2007. [48] Thomas Ryder, Andrew Golightly, A. Stephen McGough, and Dennis Prangle. Black-box variational inference for stochastic differential equations. In Proceedings of the 35th International Conference on Machine Learning (ICML), 2018. [49] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. DeepAR: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 36(3):1181–1191, 2020. [50] Patrick Seifner and Ramsés J. Sánchez. Neural Markov jump processes. In Proceedings of the 40th International Conference on Machine Learning (ICML), 2023. [51] Dario Shariatian, Umut Simsekli, and Alain Oliviero Durmus. Denoising Lévy probabilistic models. In The Thirteenth International Conference on Learning Representations, 2025. 12

[52] Kamilė Stankevičiūtė, Ahmed M. Alaa, and Mihaela van der Schaar. Conformal time-series forecasting. In Advances in Neural Information Processing Systems, volume 34, 2021. [53] Yusuke Tashiro, Jiaming Song, Yang Song, and Stefano Ermon. CSDI: Conditional score-based diffusion models for probabilistic time series imputation. In Advances in Neural Information Processing Systems, volume 34, 2021. [54] Anh Tong, Thanh Nguyen-Tang, Toan Tran, and Jaesik Choi. Learning fractional white noises in neural stochastic differential equations. In Advances in Neural Information Processing Systems, volume 35, 2022. [55] Belinda Tzen and Maxim Raginsky. Neural stochastic differential equations: Deep latent Gaussian models in the diffusion limit. arXiv preprint arXiv:1905.09883, 2019. [56] Chen Xu and Yao Xie. Conformal prediction interval for dynamic time-series. In Proceedings of the 38th International Conference on Machine Learning (ICML), 2021. [57] Eun Bi Yoon, Keehun Park, Sungwoong Kim, and Sungbin Lim. Score-based generative models with Lévy processes. Advances in Neural Information Processing Systems, 36:40694–40707, 2023. [58] Xinyu Yuan and Yan Qiao. Diffusion-TS: Interpretable diffusion for general time series generation. In The Twelfth International Conference on Learning Representations (ICLR), 2024. [59] Margaux Zaffran, Olivier Féron, Yannig Goude, Julie Josse, and Aymeric Dieuleveut. Adaptive conformal predictions for time series. In Proceedings of the 39th International Conference on Machine Learning (ICML), 2022. [60] Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu. Are transformers effective for time series forecasting? In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, 2023. [61] Shuai Zhang, Chuan Zhou, Yang Liu, Peng Zhang, Xixun Lin, and Zhi-Ming Ma. Neural jump-diffusion temporal point processes. In Proceedings of the 41st International Conference on Machine Learning (ICML), 2024.

Appendix A

Variational derivation of the optimal Markov posterior

Both Theorem 4.1 and its Corollary derive from a single key object: the optimal local transition ratio, identified by constrained variational optimisation over Markov path measures. We first present a proof sketch covering both results, then give complete step-by-step derivations. Proof sketch. We begin with a time-discretised decomposition of the KL divergence: X K−1 DKL Q ∥ P θ =

k+1 k qtk (x) DKL (x) dx,

k=0 t

k+1 k where DKL (x) =

(x′ |x)

qtk+1 ,tk (x′ |x) ln ptk+1 ,tk (x′ |x) dx′ . Incorporating the Chapman–Kolmogorov k+1

constraints via Lagrange multipliers ϕt (x), the constrained variational optimisation yields the optimal transition ratio ′ qt+∆t,t (x′ |x) eϕt (x ) . = (16) pt+∆t,t (x′ |x) Ep eϕt (Xt+∆t ) | x Theorem 4.1 (posterior generator and tilted Lévy measure). The infinitesimal generator of Eq. (3) for sufficiently smooth test functions G : Rd → R is Z P θ ⊤ Lt G(x) = ft (x) ∇G(x) + G(x + y) − G(x) ντ (dy) + 12 tr D(x) ∇2 G(x) . (17) Rd

This generator characterisation exists as Lévy processes are Markov, so their dynamics are completely specified by their infinitesimal behaviour. When σ = 0, the diffusion term vanishes and the generator reduces to a drift plus a pure integral operator. Using the optimal transition ratio to rewrite Eq [G(Xt+∆t ) | x] as a ratio of prior expectations, expanding numerator and denominator via the prior generator expansion, and equating the O(∆t) coefficient with LQ t G(x) yields the conjugate relation −ϕt (x) P ϕt (·) ϕt (·) LQ Lt e G(·) (x) − G(x) e−ϕt (x) LP (x). t G(x) = e t e Substituting the prior generator Eq. (17) into this relation and expanding the jump and diffusion parts separately yields Eq. (7); the tilted Lévy measure Eq. (6) is identified directly from the jump integrand. Corollary (KL formula). Expanding the per-step KL using the same ratio and taking ∆t → 0 gives n −ϕt (x) P ϕt (·) o t+∆t,t ϕt DKL (x) = ∆t e−ϕt (x) LP e ϕ (x) − 1 + ϕ (x) e Lt e (x) + o(∆t). t t t R Substituting the generator Eq. (17), the jump integral contributes Rd f (y, t, x) ντ (dy) per unit time BM after simplification. For the diffusion part, the Brownian component of LP . A let L 2 denote 1 BM ϕ ϕ 1 ⊤ direct calculation usingL (e ) = e 2 tr(D∇ ϕ) + 2 ∇ϕ D∇ϕ and LBM (eϕ ϕ) = eϕ (1 + ϕ) LBM ϕ + ∇ϕ⊤ D∇ϕ yields, after cancellation of the (1 + ϕ) terms, the residual contribution 1 ⊤ 2 ∇ϕt (x) D ∇ϕt (x) per unit time. Summing over time steps, converting to path integrals under Q, and taking ∆t → 0 yields Eq. (8). We now give complete step-by-step derivations. Step 1 is common to both results; Step 2 proves Theorem 4.1; Step 3 proves the Corollary. Throughout, ϕ ≡ ϕt (x) and all generator actions are evaluated at a fixed x unless otherwise noted. Step 1: Optimal transition ratio. We minimise the KL divergence subject to the constraint that Q corresponds to a valid Markov process. The Chapman–Kolmogorov equations require Z ′ qt+∆t (x ) = qt+∆t,t (x′ |x) qt (x) dx. 14

Introducing Lagrange multipliers ϕt (x′ ) for each constraint and taking the functional derivative of the augmented objective with respect to qt+∆t,t (x′ |x) yields the first-order condition qt+∆t,t (x′ |x) = ϕt (x′ ) − ln Zt (x), pt+∆t,t (x′ |x) where Zt (x) = Ep eϕt (Xt+∆t ) | x is the normalising constant imposed by the constraint. Exponentiating gives the optimal transition ratio Eq. (16). ln

Step 2: Proof of Theorem 4.1. We derive the posterior generator from the optimal transition ratio. For any test function G, the posterior generator is defined by Eq [G(Xt+∆t ) | x] = G(x) + ∆t LQ t G(x) + o(∆t). Using the optimal transition ratio Eq. (16), we can also write this expectation as a ratio of prior expectations: Ep eϕt (Xt+∆t ) G(Xt+∆t ) | x Eq [G(Xt+∆t ) | x] = . Ep eϕt (Xt+∆t ) | x Applying the prior generator expansion Ep [G(Xt+∆t ) | x] = G(x) + ∆t LP t G(x) + o(∆t) to numerator and denominator: ϕt (·) eϕ G(x) + ∆t LP G(·) (x) + o(∆t) t e . Eq [G(Xt+∆t ) | x] = ϕt (·) (x) + o(∆t) eϕ + ∆t LP t e Factoring out eϕ and performing a first-order Taylor expansion of the denominator, (1 + ϕt −1 ϕt ∆t e−ϕ LP = 1 − ∆t e−ϕ LP t e (x)) t e (x) + o(∆t), then expanding and retaining O(∆t) terms gives h i ϕt −ϕ P ϕt Eq [G(Xt+∆t ) | x] = G(x) + ∆t e−ϕ LP Lt e (x) + o(∆t). t e G (x) − G(x) e Equating the O(∆t) coefficient with LQ t G(x) yields the conjugate relation −ϕ P ϕt (·) ϕt (·) LQ Lt e G(·) (x) − G(x) e−ϕ LP (x). t G(x) = e t e

(18)

P P θ We now substitute the prior generator LP t = LJ + LBM + ft · ∇ into Eq. (18).

Jump part. Setting H ≡ Ht (x, y) = eϕt (x+y)−ϕ : Z Z −ϕ P ϕ −ϕ P ϕ e LJ (e G) − G e LJ e = HG(x + y) − G(x) ντ (dy) − G (H − 1)ντ (dy) R ZR = H G(x + y) − G(x) ντ (dy) ZR = G(x + y) − G(x) eϕt (x+y)−ϕt (x) ντ (dy). R

The tilted Lévy measure is identified directly from this integrand as ν̃(dy, t, x) eϕt (x+y)−ϕt (x) ντ (dy), establishing Eq. (6).

Diffusion part. Since ∂i ∂j (eϕ G) = eϕ [∂i ϕ∂j ϕG + ∂i ∂j ϕG + ∂i ϕ∂j G + ∂j ϕ∂i G + ∂i ∂j G]: X ϕ −ϕ P e−ϕ LP LBM (eϕ ) = 12 Dij ∂i ϕ ∂j G + ∂j ϕ ∂i G + ∂i ∂j G . BM (e G) − G e ij

With D symmetric, ij Dij ∂j ϕ ∂i G (D∇ϕ)⊤ ∇G + 21 tr(D ∇2 G).

(D∇ϕ)⊤ ∇G, so the diffusion contribution is

Drift part. The drift term passes through the conjugate relation unchanged, since e−ϕ (ftθ · ∇)(eϕ G) − G e−ϕ (ftθ · ∇)eϕ = ftθ · ∇G. Adding all three contributions confirms Eq. (7). Step 3: Proof of the Corollary. With the optimal transition ratio in hand, the per-step KL is t+∆t,t DKL (x) = Eq [ϕt (Xt+∆t ) | x] − ln Zt (x).

Applying the generator expansion Ep [G(Xt+∆t ) | x] = G(x) + ∆t LP t G(x) + o(∆t) to both terms: ϕ Ep e t ϕt (Xt+∆t ) | x Eq [ϕt (Xt+∆t ) | x] = Ep eϕt (Xt+∆t ) | x h i ϕt (·) ϕt (·) = ϕ + ∆t e−ϕ LP ϕt (·))(x) − ϕ e−ϕ LP (x) + o(∆t), t (e t e ϕt (·) ln Zt (x) = ϕ + ∆t e−ϕ LP (x) + o(∆t). t e

Subtracting gives o n t+∆t,t ϕt −ϕ P ϕt DKL (x) = ∆t e−ϕ LP Lt e (x) + o(∆t). t (e ϕt )(x) − (1 + ϕ) e

(19)

R Jump contribution. Let LP J G(x) = R (G(x + y) − G(x)) ντ (dy) be the jump part of the generator. We evaluate each term in Eq. (19) for LP J: Z h i ϕt e−ϕ LP (e ϕ )(x) = eϕt (x+y)−ϕ ϕt (x + y) − ϕ ντ (dy), t J R Z h i −ϕ P ϕt (1 + ϕ) e LJ e (x) = (1 + ϕ) eϕt (x+y)−ϕ − 1 ντ (dy). R

Setting H ≡ Ht (x, y) = eϕt (x+y)−ϕ and subtracting: ϕt −ϕ P ϕt LJ e (x) e−ϕ LP J (e ϕt )(x) − (1 + ϕ) e Z h i = Hϕt (x + y) − ϕ − (1 + ϕ)(H − 1) ντ (dy) ZR h i = Hϕt (x + y) − Hϕ − H + 1 ντ (dy) ZR h i = H(ϕt (x + y) − ϕ) − H + 1 ντ (dy) ZR h Z i = H ln H − H + 1 ντ (dy) = f (y, t, x) ντ (dy), R

where the penultimate equality uses ln H = ϕt (x + y) − ϕ. 1 2 Diffusion contribution. Let LP BM G(x) = 2 tr(D ∇ G(x)) be the Brownian part. We compute each factor by applying the product rule twice. ϕ ϕ ϕ For LP BM (e ): since ∂i ∂j e = e (∂i ϕ∂j ϕ + ∂i ∂j ϕ), ϕ ⊤ 2 1 1 e−ϕ LP BM (e ) = 2 ∇ϕ D ∇ϕ + 2 tr(D ∇ ϕ). ϕ ϕ ϕ ϕ ϕ For LP BM (e ϕ): since ∂j (e ϕ) = e ∂j ϕ(1 + ϕ), applying ∂i gives ∂i ∂j (e ϕ) = e [∂i ϕ∂j ϕ(2 + ϕ) + ∂i ∂j ϕ(1 + ϕ)], so ϕ ⊤ 2 1 1 e−ϕ LP BM (e ϕ) = (2 + ϕ) 2 ∇ϕ D ∇ϕ + (1 + ϕ) 2 tr(D ∇ ϕ).

Subtracting (1 + ϕ) times the first from the second, the tr(D ∇2 ϕ) terms cancel exactly: ϕ −ϕ P e−ϕ LP LBM (eϕ ) = 12 (2 + ϕ − (1 + ϕ)) ∇ϕ⊤D ∇ϕ = 12 ∇ϕ⊤D ∇ϕ. BM (e ϕ) − (1 + ϕ) e

Combining and taking the limit. Adding the jump and diffusion contributions, the full per-step KL Eq. (19) becomes Z t+∆t,t DKL (x) = ∆t f (y, t, x) ντ (dy) + 12 ∇ϕt (x)⊤D ∇ϕt (x) + o(∆t). R

Summing over time steps and converting the sum to a path integral under Q via EQ [ · ] in the limit ∆t → 0 yields Eq. (8). 16

k qtk (x)∆t →

Remark A.1 (Formal status of the continuous-time limit). The derivation passes from the discretetime KL decomposition to the continuous-time path integral Eq. (8) via ∆t → 0. A fully rigorous treatment requires two things: (i) establishing absolute continuity of Q with respect to P θ on path space, so that dQ/dP θ is well-defined; and (ii) verifying convergence of the discrete KL sum to the path integral. For (i), the Radon-Nikodym derivative takes the form of a Doléans-Dade stochastic exponential of a local martingale, which must be shown to be a true martingale rather than merely a local one. For Brownian SDEs, the Girsanov theorem accomplishes this via the Novikov condition " !# Z T ⊤ 1 EP exp 2 ∇ϕt (Xt ) D ∇ϕt (Xt ) dt < ∞; 0

the exponent matches the integrand of the quadratic KL term in Eq. (8). For pure-jump processes, the analogous requirement is a Novikov-type condition on the stochastic exponential now driven by the tilting factor Ht rather than ∇ϕt ; the complete theory is developed in [23]. The At < 0 constraint ensures the tilted measure has finite total mass and all finite moments, consistent with such integrability requirements, but a complete formal proof is left to future work.

Conditionally Gaussian representations of Lévy measures

It can be shown that a disintegration of a Lévy measure ν(dy) can be formed such that [46] Z ∞ ν(dy) = σ(r; dy) π(r) dr, 0

where σ(·; ·) is a probability kernel for fixed r and π(r) is the Lévy measure of a subordinator process. Hence ν(dy) is modelled as a mixture measure with mixing measure π(r). This implies that a sequence of jumps {yi } can be represented as conditionally independent random variables given {ri }. Symmetric stable processes. For the symmetric α-stable case the mixing measure is π(r) = r−1−α (corresponding to a non-negative stable subordinator). Choosing the probability kernel as a zero-mean 2 Gaussian N 0, r2 σG [32], the disintegration for a single jump y is Z ∞ −1−α 2 ν(y) dy ∝ N y, r2 σG r dr dy. 0

In practice we use a truncated mixing measure πτ (r) = ατ α r−1−α (supported on [τ, ∞)), yielding a proper probability density. Tilted symmetric stable processes. The Lévy measure of the tilted process is ν̃(y, t, Xt ) = eϕt (Xt +y)−ϕt (Xt ) ν(dy). Introducing the conditionally Gaussian disintegration and the quadratic parametrization of Eq. (9), the normalizing constant of the conditional kernel is Z 2 C(r, t, Xt ) = eϕt (Xt +y)−ϕt (Xt ) N y, r2 σG dy. R

With K1 = 2At Xt + Bt and K2 = At − 2r21σ2 , completing the square gives the closed form in G the normalizing constant for the truncated tilted mixing measure is K(α, t, Xt ) = R ∞ main text. The C(r, t, Xt ) r−1−α dr, and the resulting probability density is τ Z ∞ ν̃τ (y, t, Xt ) = σ̃(r; y, t, Xt ) π̃τ (r, t, Xt ) dr, τ R with σ̃ and π̃τ as defined in Theorem 4.4. One can verify R ν̃τ (dy, t, Xt ) = 1 by exchanging the order of integration and using the fact that σ̃ is a probability kernel.

Exact rejection sampler for the tilted mixing measure

Proof of Theorem 4.6. We derive a uniform envelope for C(r, t, Xt ) under the quadratic parametrisation Eq. (9). From Rem. 4.5, with K1 = 2At Xt + Bt and K2 = At − 2r21σ2 , G 1 K12 C(r, t, Xt ) = p exp − . 2 4K2 −2K2 r2 σG 17

2 2 Introduce u = 2|At |r2 σG ≥ 0. Since At < 0, we have −2K2 r2 σG = 1 + u and −K12 /(4K2 ) = 2 K1 u / (4|At |(1 + u)), giving 2 K1 u 1 exp · . C(r, t, Xt ) = √ 4|At | 1 + u 1+u √ For all u ≥ 0, both 1/ 1 + u ≤ 1 and u/(1 + u) ≤ 1, so 2 K1 C(r, t, Xt ) ≤ exp =: M (t, Xt ). 4|At |

The bound M (t, Xt ) is independent of r and hence a valid envelope. Proposing r∗ ∼ q(r) ∝ r−1−α 1r≥τ via the inverse-CDF transform and accepting with probability C(r∗ , t, Xt )/M (t, Xt ) yields an exact draw from π̃τ (r, t, Xt ) ∝ C(r, t, Xt ) r−1−α by standard rejection-sampling correctness.

Analytical integration of the KL integrand

To compute the KL loss component at each time step, the following integral must be evaluated: Z ∞ f (y, t, Xt ) + f (−y, t, Xt ) y −1−α dy, (20) I(t, Xt ) = τ

where f (y, t, Xt ) = Ht (Xt , y) ln Ht (Xt , y) − Ht (Xt , y) + 1 and Ht (x, y) = exp(At (2xy + y 2 ) + Bt y) under the quadratic parametrization. Setting K1 = 2At Xt + Bt , the integrand decomposes into four components: Z ∞ 2 I1+ (t, Xt ) = K1 eK1 y+At y y −α dy, Z τ∞ 2 + I2 (t, Xt ) = At eK1 y+At y y 1−α dy, Z ∞τ 2 I3+ (t, Xt ) = eK1 y+At y y −1−α dy, τ Z ∞ τ −α I4+ (t, Xt ) = y −1−α dy = . α τ Each of I1± , I2± , I3± can be expressed as a convergent series using the substitution η = At y 2 (with P∞ K n At < 0) and upper incomplete gamma functions. Expanding eK1 y = n=0 n!1 y n yields I1+ (t, Xt ) =

∞ n−α+1 K1 X K1n (−At )− 2 Γ n−α+1 , −At τ 2 , 2 2 n=0 n!

and analogously for the remaining terms. The total integral at each step is I(t, Xt ) = I + (t, Xt ) + I − (t, Xt ), where I ± denotes the contributions from +y and −y respectively, truncated to a finite number of series terms in practice. Comparison with Monte Carlo estimation. This series representation provides an exact alternative to Alg. 1. In practice, however, for the small truncation thresholds τ usedpin our experiments, the series estimator exhibits larger variance than direct Monte Carlo. When |K1 |/ |At | is non-negligible, early terms in the series can be large in magnitude before the factorial denominator dominates, causing numerical instability that persists even after series truncation. Alg. 1 is therefore used throughout this paper; the series representation is retained here as a reference and may be preferable when τ is large relative to the jump scale.

Extension to state-dependent jump sizes

The derivation in Sec. 4.1 assumes that jump sizes enter the prior SDE additively and independently of the current state. Here we show that the variational framework extends to the more general prior Z dXt = ftθ (Xt ) dt + σ(Xt ) dBt + γ(Xt− ) y Ñ (dy, dt), (21) R

where γ : R → R>0 is a smooth state-dependent scaling function and Ñ is the compensated Poisson random measure with intensity ντ (dy) dt. The infinitesimal generator of Eq. (21) is Z P θ 2 1 Lt G(x) = ft (x) ∇G(x) + 2 tr(D ∇ G(x)) + G(x + γ(x)y) − G(x) ντ (dy). (22) R

Proposition E.1 (Tilted variational family under state-dependent jump sizes). For a prior with generator Eq. (22), the KL divergence between the prior and the optimal Markov posterior is "Z Z # Z T 1 T ⊤ θ ∇ϕt (Xt ) D ∇ϕt (Xt ) dt , (23) DKL Q ∥ P = EQ fγ (y, t, Xt ) ντ (dy) dt + 2 0 0 R where fγ (y, t, x) = Hγ ln Hγ − Hγ + 1 with Hγ ≡ Htγ (x, y) = eϕt (x+γ(x)y)−ϕt (x) . The posterior generator is ⊤ θ 2 1 LQ t G(x) = ft (x) + D∇ϕt (x) ∇G(x) + 2 tr(D ∇ G(x)) Z (24) + G(x + γ(x)y) − G(x) eϕt (x+γ(x)y)−ϕt (x) ντ (dy), R

and the variational SDE is dXt = ftθ (Xt ) + D∇ϕt (Xt ) dt + σ(Xt ) dBt +

γ(Xt− ) y Ñ ϕ (dy, dt),

(25)

where Ñ ϕ is the compensated eϕt (Xt− +γ(Xt− )y)−ϕt (Xt− ) ντ (dy) dt.

Poisson

random

measure

with

tilted

intensity

Proof. The proof follows Steps 1 and 3 of Sec. A verbatim, with the sole substitution G(x + y) → G(x + γ(x)y) throughout the jump generator. The Brownian contribution is Runchanged. For the jump contribution, the same cancellation as the jump calculation in Step 3 gives R fγ (y, t, x) ντ (dy) with Hγ = eϕt (x+γ(x)y)−ϕt (x) in place of H. Simulation. The conditionally Gaussian disintegration of Theorem 4.4 adapts to the state-dependent case: the tilted kernel becomes a Gaussian in y with mean and variance scaled by γ(x)−1 and γ(x)−2 respectively. The quadratic parametrization and the rejection sampler of Sec. 4.4 remain tractable under this scaling with straightforward modifications. All experiments in this paper use γ(·) ≡ 1, which recovers the formulation of Sec. 4; adapting the simulation to γ(·) ̸= 1 is left to future work.

Connection to score-based generative modelling

The variational family has a precise structural connection to score-based generative modelling. For Gaussian SDEs the drift correction D∇ϕt equals the score of the marginal density, enabling learned reverse processes. The analogous formula for Lévy systems replaces the jump measure by [pt (x + y)/pt (x)] ν(dy); at optimum, the tilted measure derived here reduces to exactly this ratio, identifying neural tilting as the structural backbone of Lévy diffusion generative models. Two recent works have constructed Lévy generative models. [57] derive an exact reverse-time SDE for an isotropic α-stable forward process, in which the learnt correction enters only the reverse drift via a fractional score function of the current state, whilst the driving noise in the reverse process remains a fixed isotropic α-stable process. [51] take a discrete-time route, replacing Gaussian increments in DDPM with α-stable ones and using a normal variance-mixture representation of α-stable random variables to recover tractable backward kernels; here too the background jump law is fixed and the generative capability resides in the denoising function. A generative model built on the variational family of Theorem 4.1 would differ structurally: the reverse process would carry a state-conditioned tilted jump measure eϕt (x+y)−ϕt (x) ν(dy), so the jump law itself adapts to the current state rather than remaining a fixed isotropic noise corrected by a learnt drift. Remark F.1 (Joint parametrization in the jump-diffusion setting). In the joint setting (σ ̸= 0), the Brownian drift correction D∇ϕt (Xt ) and the jump tilting factor eϕt (Xt +y)−ϕt (Xt ) are both determined by the same potential ϕt . This coupling is a structural consequence of Theorem 4.1: independently parametrizing the two channels breaks the shared ϕt structure and forfeits the optimality 19

guarantee. Under the quadratic parametrization of Sec. 4.2, both the gradient ∇ϕt and the increment ϕt (x + y) − ϕt (x) are linear functions of the same neural outputs (At , Bt ), so joint estimation of the Brownian and jump corrections is automatic. Developing richer parametrizations that balance the two channels with greater flexibility while preserving the Markov optimality structure is a natural direction for future work.

Implementation Details

All tilted stable (TS) and Gaussian SDE experiments were implemented in JAX [6], and parameter updates were carried out using the Optax optimisation library [14]. Financial data were sourced via the yfinance API [2]. G.1

Wall-clock training time

To give a practical sense of the computational cost, we report representative wall-clock training times measured on a single NVIDIA GH200 Grace–Hopper GPU. For the synthetic experiments, times are reported per dataset realisation; for the financial experiments, times are reported per rolling train-and-forecast window. The TS and Gaussian SDE models use the same optimisation protocol, so their runtimes are directly comparable and isolate the additional cost of heavy-tailed posterior inference. The remaining baselines are included as context, though they solve different optimisation problems and belong to different model classes. The additional cost of TS relative to the Gaussian SDE is consistent across settings and reflects the overhead of sampling from the tilted Lévy measure and estimating the jump contribution to the ELBO. In the synthetic setting, TS is approximately 5× slower than the Gaussian SDE while remaining highly stable across both drift families. In finance, the univariate TS model requires approximately 84 minutes per rolling window, increasing to approximately 247 minutes in the d = 10 multivariate setting. This scaling is consistent with the added cost of simulating and optimising heavytailed jump structure in higher dimension.

Table 3: Representative single-GPU wall-clock training times. Synthetic times are reported per run (3000 optimisation steps). Financial times are reported per rolling window. The multivariate NJ-SDE time is estimated from the 90 of 302 windows that completed. Setting

Model

Time

Synthetic (OU / double well) Synthetic (OU / double well)

TS (ours) Gaussian SDE

∼ 122 min/run ∼ 25–26 min/run

Financial, univariate Financial, univariate Financial, univariate Financial, univariate Financial, univariate Financial, univariate Financial, univariate

TS (ours) Gaussian SDE DeepAR NJ-SDE Neural MJD N-HiTS DLinear

∼ 84 min/window ∼ 20 min/window ∼ 1.2 min/window ∼ 2.4 min/window ∼ 0.3 min/window ∼ 0.3 min/window < 1 s/window

Financial, multivariate (d = 10) Financial, multivariate (d = 10) Financial, multivariate (d = 10) Financial, multivariate (d = 10) Financial, multivariate (d = 10) Financial, multivariate (d = 10) Financial, multivariate (d = 10)

TS (ours) Gaussian SDE DeepAR NJ-SDE Neural MJD N-HiTS DLinear

∼ 247 min/window ∼ 34 min/window ∼ 3.3 min/window ∼ 11 min/window ∼ 0.3 min/window ∼ 0.7 min/window ∼ 0.1 min/window

Although TS is substantially slower than forecasting-oriented baselines such as DLinear, N-HiTS, and Neural MJD, these models do not perform posterior inference over a Lévy-driven latent path measure. We therefore view the runtime comparison with the Gaussian SDE as the most informative like-for-like measure of the computational price of modelling heavy-tailed jump dynamics. G.2

Optimisation details

Both TS and Gaussian SDE models use the same neural drift parametrisation ftθ : a one-hiddenlayer MLP of width 32. We found that increasing the drift width consistently degraded predictive performance, suggesting that a small drift network acts as an effective regulariser by preventing the drift from absorbing variability that should instead be attributed to the latent noise process. Tilted Stable architecture. The TS variational posterior uses the quadratic parametrisation from Sec. 4.2, with the coefficients At and Bt generated by MLPs of width 256 and 5 hidden layers. The temporal encoder uses 100 learnable reference times and embedding dimension 64. The lower-bound constant in the reparametrisation of At is fixed at amin = 0.001. The truncation threshold is fixed at τ = 0.01 across all experiments. 20

Remark G.1 (Role and sensitivity of τ ). Setting τ > 0 replaces the infinite-activity component of the prior with a finite compound Poisson process, which is required for the variational framework of Sec. 4.1 to be well-posed without compensators. The heavy-tailed character of the prior is not affected: as established in Dfn. 3.5, the power-law tail and the tail index α are preserved for any τ > 0. The contribution of the discarded small jumps (|y| < τ ) is well-approximated by Gaussian noise [3], as discussed in Rem. 3.4. Consequently, results are insensitive to the precise value of τ provided it is small relative to the typical jump scale; τ = 0.01 satisfies this condition in all experimental settings considered here. Discretisation and Monte Carlo estimation. In the synthetic experiments, both TS and Gaussian SDE models are trained for 3000 optimisation iterations. In the financial experiments, the univariate TS model is also trained for 3000 iterations per rolling window, while the multivariate TS model uses 10000 iterations because convergence is substantially slower in the d = 10 setting. The Gaussian SDE is trained for 3000 iterations in all experiments, which we found sufficient for reliable convergence. All remaining baselines are likewise trained to convergence, which in practice was achieved within 3000 iterations. In each iteration, the ELBO (or Gaussian analogue) is estimated using 500 independent simulated latent paths, and the Euler discretisation uses 1000 latent time steps between the start and end of the observed training window. For TS, the KL jump integral is approximated using 1000 Monte Carlo samples per time step as described in Sec. 4.3. The Gaussian SDE uses the same Euler discretisation and numbers of simulated paths. Regularisation. We apply ℓ2 regularisation to all trainable parameters for every model. The regularisation scale is tuned manually as a model-specific hyperparameter for all methods considered in the experiments. Optimiser and gradient stabilisation. The TS model is trained with RMSProp and no learning-rate decay. The learning rate is set to 10−4 in the synthetic and univariate financial experiments, and to 10−3 in the multivariate financial experiment. We found the absence of decay important in practice: with decay, the combination of heavy-tailed Monte Carlo noise and already conservative secondmoment normalisation often led to updates that became too small late in training. Compared with Adam at similar learning rates, RMSProp produced materially more stable optimisation, whereas Adam frequently led to pronounced oscillations. A plausible explanation is that under heavy-tailed Monte Carlo noise, first-moment momentum amplifies rare but very large gradient excursions instead of damping them, while RMSProp’s second-moment normalisation is more conservative. The Gaussian SDE and all remaining benchmarks are trained with Adam using learning rate 10−4 and exponential decay factor 0.95. In addition, TS uses a custom layerwise gradient rescaling step before the RMSProp update. For each parameter tensor g, we compute its robust p scale via the empirical 0.95-quantile of |g|, and compare it with the root-mean-square norm ∥g∥2 / |g|. The tensor is then rescaled by ! p ∥g∥2 / |g| g ← g / max 1, . q0.95 (|g|) + ε This transformation preserves the relative magnitudes of coordinates within a layer, unlike coordinatewise clipping, but suppresses updates when the overall layerwise norm is dominated by a small number of extreme entries. Empirically, this was important for stable optimisation of the TS objective, whose heavy-tailed jump samples can induce occasional gradient spikes. We interpret the method as a robust compromise between no clipping and hard clipping: it retains directional information while reducing the influence of rare, disproportionately large gradients.

Additional results on synthetic data

Tabs. 4 to 6 report the full synthetic jump-CRPS results at σε = 0.10 and the full synthetic results at σε = 0.05. The qualitative findings mirror those at σε = 0.10: our model achieves lower CRPS than the Gaussian baseline at every α in both systems, substantially reduces θ and θ2 recovery error, and leads on jump CRPS at all thresholds. The µ (OU) and θ1 (double well) recovery results follow the same pattern as at σε = 0.10: the Gaussian baseline performs better on these parameters at most α values, reflecting the trade-off between the tilting function and the global drift parameters discussed 21

Table 4: Jump CRPS (σε = 0.10) at increasing percentile thresholds for each α; lower is better. Best result per metric in bold. Ornstein–Uhlenbeck p90

p95

Double well

p97.5

p99

p90

p95

p97.5

p99

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

1.40 0.90 1.10 0.53 0.55 0.39 0.43 0.84 0.45

1.15 0.38 0.97 0.39 0.39 0.26 0.31 0.24 0.24

1.70 1.00 1.63 0.63 0.66 0.41 0.47 0.85 0.48

1.52 0.46 1.54 0.53 0.50 0.29 0.38 0.26 0.27

2.15 1.12 2.35 0.82 0.86 0.44 0.52 0.91 0.51

2.10 0.56 2.36 0.79 0.71 0.33 0.43 0.28 0.31

2.46 1.23 3.01 1.25 1.14 0.52 0.58 0.98 0.54

2.41 0.67 3.20 1.26 0.98 0.38 0.48 0.34 0.33

0.54 0.42 0.43 0.43 0.39 0.40 0.40 0.41 0.41

0.31 0.22 0.23 0.23 0.20 0.20 0.21 0.20 0.21

0.66 0.48 0.48 0.49 0.42 0.43 0.43 0.43 0.43

0.42 0.28 0.27 0.29 0.23 0.23 0.23 0.23 0.23

0.85 0.58 0.55 0.60 0.47 0.47 0.47 0.47 0.46

0.60 0.37 0.35 0.38 0.29 0.27 0.28 0.26 0.26

1.35 0.85 0.77 0.85 0.61 0.59 0.56 0.55 0.52

1.05 0.62 0.54 0.60 0.41 0.38 0.38 0.34 0.32

All α 0.73 0.48 0.87 0.64 1.08 0.87 1.30 1.12 0.43 0.22 0.47 0.27 0.55 0.34 0.74 0.52

Table 5: Synthetic experiment results (σε = 0.05), reporting means over 50 realisations per α. Held-out CRPS and mean absolute parameter recovery errors for the OU and double-well (DW) systems; lower is better. Best result per metric in bold. Ornstein–Uhlenbeck ∗

CRPS ↓

|θ̂−θ | ↓

Double well ∗

|µ̂−µ | ↓

|θ̂1 −θ1∗ | ↓

CRPS ↓

|θ̂2 −θ2∗ | ↓

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

0.99 0.81 0.55 0.42 0.39 0.35 0.36 0.80 0.41

0.63 0.30 0.38 0.25 0.23 0.21 0.22 0.22 0.20

5.13 7.59 9.01 7.48 7.18 6.16 9.61 8.57 6.98

2.17 2.55 2.42 1.91 2.64 2.41 2.65 2.51 1.68

3.99 3.16 2.67 2.60 2.83 3.54 2.20 2.05 2.02

7.07 3.80 5.14 5.12 4.26 3.74 2.57 3.22 2.56

0.41 0.36 0.34 0.36 0.33 0.34 0.36 0.40 0.39

0.15 0.13 0.12 0.12 0.11 0.12 0.13 0.15 0.14

1.10 0.92 1.73 1.42 1.16 0.98 0.82 1.02 1.72

0.83 1.10 1.39 1.62 1.88 1.82 2.25 1.26 1.32

1.93 2.17 3.35 3.26 3.42 3.08 2.91 2.21 2.37

0.24 0.29 0.33 0.31 0.44 0.30 0.44 0.29 0.51

All α 0.56 0.29 7.52 2.33 2.78 4.16 0.37 0.13 1.21 1.50 2.75 0.35

in Sec. 5.1. With lower observation noise the jump CRPS advantage for the double-well system is now fully resolved (no missing entries). Fig. 3 shows posterior sample paths for nine additional realisations, further illustrating the consistent advantage of the tilted-stable model across seeds and both dynamical systems.

Table 6: Jump CRPS (σε = 0.05) at increasing percentile thresholds for each α; lower is better. Best result per metric in bold. Ornstein–Uhlenbeck p90

p95

Double well

p97.5

p99

p90

p95

p97.5

p99

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

1.40 0.90 1.10 0.53 0.55 0.39 0.43 0.84 0.45

1.15 0.38 0.97 0.39 0.39 0.26 0.31 0.24 0.24

1.70 1.00 1.63 0.63 0.66 0.41 0.47 0.85 0.48

1.52 0.46 1.54 0.53 0.50 0.29 0.38 0.26 0.27

2.15 1.12 2.35 0.82 0.86 0.44 0.52 0.91 0.51

2.10 0.56 2.36 0.79 0.71 0.33 0.43 0.28 0.31

2.46 1.23 3.01 1.25 1.14 0.52 0.58 0.98 0.54

2.41 0.67 3.20 1.26 0.98 0.38 0.48 0.34 0.33

0.55 0.46 0.41 0.43 0.39 0.37 0.41 0.43 0.42

0.27 0.22 0.19 0.18 0.17 0.16 0.18 0.18 0.18

0.67 0.55 0.47 0.49 0.44 0.40 0.44 0.46 0.45

0.38 0.30 0.24 0.24 0.21 0.19 0.21 0.21 0.21

0.89 0.72 0.57 0.60 0.52 0.46 0.50 0.51 0.49

0.57 0.44 0.34 0.33 0.29 0.24 0.27 0.26 0.26

1.41 1.12 0.82 0.85 0.72 0.58 0.64 0.62 0.60

1.03 0.80 0.57 0.55 0.47 0.36 0.40 0.38 0.38

All α 0.73 0.48 0.87 0.64 1.08 0.87 1.30 1.12 0.43 0.19 0.49 0.24 0.58 0.33 0.82 0.55

Additional financial results

Tab. 7 reports results for GOOGL, using the same evaluation protocol as NVDA (304 rolling windows, 2-day forecast horizon). N-HiTS was not evaluated on GOOGL. The ranking is consistent with NVDA: our model achieves the lowest CRPS and leads on jump CRPS at every threshold. Reliability diagrams for GOOGL and the multivariate (d=10) setting are shown in Fig. 2. Table 7: GOOGL financial forecasting results averaged over 304 evaluation windows; lower is better for all metrics. Best in bold. † Deterministic model; CRPS equals MAE. Jump CRPS is computed on price increments exceeding the indicated percentile threshold. N-HiTS was not evaluated on GOOGL. Jump CRPS ↓ Model NJ-SDE DLinear† Gaussian SDE Neural MJD DeepAR TS (ours)

CRPS ↓

p90

p95

p97.5

p99

1.199±1.258 0.745±0.609 0.610±0.670 0.542±0.454 0.536±0.512 0.435±0.419

1.609 1.257 1.166 0.999 1.002 0.843

1.805 1.498 1.369 1.184 1.223 1.053

2.312 1.919 1.690 1.455 1.529 1.430

2.570 2.240 1.827 1.881 1.823 1.709

Figure 3: Posterior sample paths for nine additional synthetic realisations, comparing the Gaussian SDE and the tilted-stable model on the OU (left pair) and double-well (right pair) systems. Each row corresponds to an independent realisation; the layout within each panel follows Fig. 1.

Record · ID 175238 · SHA-256 dbcb4414c04cb4ac
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.