Implementing Neural SDEs for Data-Driven Dynamics of the Bitcoin Option Surface

1. Introduction

1.1. Motivation

Accurately modeling the dynamics of options is a core function of modern quantitative finance, underpinning both valuation and risk management in numerous asset classes. Traditional approaches to option pricing, such as Black-Scholes (Black and Scholes 1973), Heston (Heston 1993) and SABR (Hagan et al. 2002) have been extensively developed and tested in the context of mature financial markets such as equities or foreign exchange. As with any modeling approach, these frameworks have simplifying assumptions, in particular in the form of parametric assumptions for the drift and diffusion processes.

Cryptocurrency markets, and by extension their derivatives markets, represent a markedly different environment than established financial markets. It is because of these altered dynamics that the assumptions in traditional options pricing models may be even more restrictive when applied to cryptocurrency as opposed to traditional asset classes. Ahmed et al. (2024) provide a review paper which highlights that high volatility seems to be a structural feature in cryptocurrency markets. They also suggest that machine learning methods, in particular deep learning models, may better capture non-linear dynamics and handle high-frequency data.

Furthermore, there has been little work done in the space of cryptocurrency option pricing within the literature. While all options modelling techniques are domain agnostic (i.e., the market itself doesn’t matter, the mathematics of the models do not change), as far as we are aware there has not been any work done on the modelling of a multi-strike, multi-maturity options book (i.e., the so-called volatility surface) for cryptocurrencies.

1.2. Goals

The broad goal of this project is to use raw options quotes, obtained from the Deribit exchange, and create a data-driven model for the joint dynamics of the options observable in the market. This can be broken down into three points.

• Data driven method to model dynamics of the observable call options surface. We use non-parametric machine learning methods to extract the model dynamics from the observable call options on Deribit. Rather than assuming a parametric form for the drift and diffusion of the underlying, the model seeks to learn the evolution of option prices directly from historical time-series data across strikes and maturities. By focusing on the evolution of the entire surface rather than individual options, the model captures the co-movements and structural relationships inherent in a live options market, enabling a more realistic depiction of the dynamics observed in high-frequency cryptocurrency data.
• Ensure the model does not admit arbitrage. Absence of arbitrage is a basic consistency requirement for financial models. We want to ensure our model accounts for this by ensuring that outputs admitting arbitrage are either eliminated entirely, or at the least minimised. To ensure theoretical and practical validity, the model is explicitly designed to satisfy static, and minimise any violation of dynamic, no-arbitrage conditions. Static arbitrage constraints impose shape restrictions on the option surface—such as monotonicity in strike and convexity, and calendar spread consistency—ensuring that prices at a single point in time form a feasible and economically consistent surface. This is the more fundamental requirement in the sense that violations of static arbitrage constraints lead to immediate, model–free arbitrage opportunities. Dynamic arbitrage constraints, in contrast, govern the temporal evolution of option prices, ensuring that the model’s drift and diffusion dynamics are consistent with an equivalent martingale measure. What each of these constraints mean definitionally will be discussed in more depth in Section 2.2 and Section 2.3 below. The static conditions are incorporated directly into the learning process through constrained optimisation and diffusion shrinkage methods derived from the Arbitrage-Free Neural SDE framework of Cohen et al. (2023a).
• Use of the model for risk management and hedging. As application examples of the model fitted to the market in the present paper, we build on the work of Cohen et al. (see Cohen et al. (2022) and Cohen et al. (2023b)) to hedge options portfolios and provide value at risk (VaR) estimates via market simulation.

1.3. Significance and Contributions

This work contributes to the emerging literature on data-driven financial modelling by applying the arbitrage-free Neural SDE framework to the cryptocurrency derivatives market. It provides the first empirical application of an arbitrage-consistent neural SDE to intraday, multi-strike, multi-maturity crypto options data. Beyond methodological adaptation, the model demonstrates practical relevance by producing risk and hedging metrics that are derived from the same underlying stochastic dynamics, bridging the gap between theoretical consistency and applied market utility.

1.4. Overview of Dataset Used

The dataset has been obtained from the Deribit exchange, which is a cryptocurrency options exchange, where at each timestep we see every available option on Bitcoin from the exchange. This provides effectively a timeseries of options prices, though we need to perform some pre-processing to create a static grid of strikes and maturities for the time-series. This will be explored more in depth in Section 3. An important feature to note of the raw data is that it is high-frequency, meaning that each consecutive snapshot is provided at sub-second level, see Figure 1. There are $118,193$ unique timestamps in one dataset which spans 24 hours.

The date for the market data used in the initial experiment is May 22nd 2025. For context, the price of BTCUSD on that day is plotted in Figure 2.

We also perform an extended analysis on a dataset following the same schema, i.e., over a week’s worth of data as opposed to a day. This extended analysis covers the dates 7th to 14th (inclusive) of September 2025. The distribution of timestep lengths is seen in Figure 3 and the movement of underlying asset is seen in Figure 4. This dataset has 842,940 unique timesteps.

2. Background

2.1. Market Models

The present paper follows the approach of so-called market models, which model the joint dynamics of all liquidly tradable assets simultaneously and impose hard and/or soft constraints on the modelling process to ensure absence of arbitrage. This idea is similar to the interest rate modelling framework of Heath-Jarrow-Morton (HJM) (Heath et al. 1992), though models for options are inherently more complex due to the constrained state-space of options prices created through no-arbitrage relationships.

The literature for market models contains works that have been done in special cases. Schönbucher (1999) and Lyons (1995) have proposed a model for the implied volatility of single options. For example, denote by $S_t$ the price of the underlying asset at time t, $C_t$ the price of a call option, K the strike price and T the maturity. If there is only one option to model, one has the static bounds ${(S_t-K)}^{+}\le C_t\le S_t$ and $X=C_T={(S_T-K)}^{+}$, where the latter corresponds to the terminal pay-off condition. The authors of these papers do this by reparameterising the option price by its implied volatility. By contrast, modelling the joint dynamics of multiple options across strikes and maturities introduces a much richer set of static arbitrage constraints that prices must satisfy (Carr and Madan 2005,Davis and Hobson 2007,Gatheral 2011). To our knowledge, only two substantial contributions address this setting: the arbitrage-consistent approach of Wissel (2008) and the data-driven neural SDE framework of Cohen et al. (2023a). The latter introduces the concept of a “liquid lattice”—a structured, yet arbitrage-consistent, state space for option prices—allowing dynamic modelling of an entire option surface rather than a single contract. This addresses the need for having a data-driven model for options that can take in a time-series of data, so the present paper builds on the methods provided by Cohen et al. (Cohen et al. 2020,Cohen et al. 2022,Cohen et al. 2023b).

2.2. Static Arbitrage

Static arbitrage refers to the ability to make a risk free profit by exploiting the prices of options at a specific point in time (Carr and Madan 2005,Davis and Hobson 2007,Gatheral 2011). Static arbitrage strategies are immediate and model–free; thus there is a strong requirement that option price surfaces should not allow static arbitrage. Effectively, these static arbitrage constraints are a set of shape constraints on the price surface. Specifically, for call options the conditions are:

1.

Monotonicity in maturity. Whenever two call prices share the same moneyness m but have different maturities ${\tau }_i<{\tau }_j$, the call price must satisfy

$$C({\tau }_j,m)\ge C({\tau }_i,m).$$

2.

Monotonicity in moneyness. At any fixed maturity $\tau$, when two moneyness values satisfy $m_k<m_{\ell }$, the corresponding call prices must satisfy

$$C(\tau ,m_k)\ge C(\tau ,m_{\ell }),$$

These conditions are derived from the absence of ’calender spread’ arbitrage formalised by Carr and Madan (2005) and further reviewed by Gatheral (2011).

3.

Convexity in moneyness. Also at a fixed $\tau$, for any three consecutive moneyness points $m_{k-1}<m_k<m_{k+1}$, convexity requires (see e.g. Breeden and Litzenberger 1978 and Roper and Rutkowski 2009)

$$C(\tau ,m_{k-1})-2C(\tau ,m_k)+C(\tau ,m_{k+1})\ge 0.$$

4.

Nonnegativity. For every node $({\tau }_i,m_i)$, the call price must be nonnegative:

$$C({\tau }_i,m_i)\ge 0.$$

As in Cohen et al. (2020), we can collect all the linear inequalities arising from conditions (1)–(3) into a single matrix inequality

$$A\mathbf{c}\ge \mathbf{0},$$

giving a set of R linear constraints in ${\mathbb{R}}^n$. In particular, each row of A corresponds to exactly one of the following three types of finite-difference expressions:

$$\begin{array}{cc}& \left(\mathrm{a}\right)\mathrm{If}({\tau }_j,m)\mathrm{and}({\tau }_i,m)\mathrm{are}\mathrm{adjacent}\mathrm{in}\tau ,C({\tau }_j,m)-C({\tau }_i,m)\ge 0,\hfill \\ & \left(\mathrm{b}\right)\mathrm{If}(\tau ,m_k)\mathrm{and}(\tau ,m_{\ell })\mathrm{are}\mathrm{adjacent}\mathrm{in}m,C(\tau ,m_k)-C(\tau ,m_{\ell })\ge 0,\hfill \\ & \left(\mathrm{c}\right)\mathrm{If}(\tau ,m_{k-1}),(\tau ,m_k),(\tau ,m_{k+1})\mathrm{are}\mathrm{consecutive}\mathrm{in}m,\hfill \\ & C(\tau ,m_{k-1})-2C(\tau ,m_k)+C(\tau ,m_{k+1})\ge 0.\hfill \end{array}$$

Meanwhile, the nonnegativity requirement is captured by the coordinate-wise condition $\mathbf{c}\ge \mathbf{0}$. We use this parameterisation in the data preprocessing step, for detecting and removing arbitrage in our historical option prices. Similar arbitrage-cleaning procedures are applied in Cohen et al. (2020) and Wissel (2008) to enforce shape constraints before model calibration.

2.3. Dynamic Arbitrage

While static arbitrage concerns the shape of the option price surface at a single point in time, dynamic arbitrage refers to risk–free profits that can be generated by trading options (and the underlying asset) over time as their prices evolve. In other words, even if an option surface satisfies all static no–arbitrage conditions at each date, its temporal evolution may still violate arbitrage–free dynamics if there is an inconsistent drift or diffusion structure across maturities and strikes (Derman and Kani 1994,Dupire et al. 1994,Fukasawa 2011).

Formally, a family of call price processes ${\left\{C_t(\tau ,m)\right\}}_{t\ge 0}$ is ensured to be dynamically arbitrage–free if there exists a risk–neutral probability measure $\mathbb{Q}$ under which the discounted price processes of the underlying and of all replicable portfolios of options are martingales. This implies that:

$$e^{-{\int }_0^t{r}_{s}ds}{S}_t\mathrm{and}{e}^{-{\int }_0^t{r}_{s}ds}{C}_t(\tau ,m)\mathrm{are}\mathbb{Q}-\mathrm{martingales}.$$

2.4. Neural Networks in Options Pricing & Modelling

The use of neural networks in options pricing was first pursued by Hutchinson et al. (1994), who trained feedforward networks to learn the pricing function of S&P 500 index options. The resurgence of neural networks in the 2010s brought new techniques for option pricing. Major advances include Neural SDEs, generative models, and the use of physics informed neural nets (PINNs).

2.4.1. Dynamic Hedging via Deep Learning

The concept of using deep learning for hedging, or ’deep hedging’ was introduced by Buehler et al. (2019). Here, the hedging problem is treated as a reinforcement learning task (Sutton et al. 1998), where an agent (the hedging strategy) learns to optimise a reward or risk objective over time. Buehler et al. (2019) used a synthetic dataset, with transaction costs, and found that the learned hedging strategy outperforms the traditional Black-Scholes delta hedging in terms of risk-adjusted returns. Similar work by Kolm and Ritter (2019) applied reinforcement learning (RL) to dynamic replication and hedging, training a hedging agent that learns how to optimally trade to manage option risk. These RL-based approaches mark a transition from using neural networks merely as function approximators to using them as decision-makers in a sequential setting.

2.4.2. Generative Models

Another use of neural networks is generative models, which seek to capture the joint distribution of asset prices and implied volatility surfaces, enabling realistic market simulation and data augmentation. Generative Adversarial Networks (GANs) (Goodfellow et al. 2020), and other related generative models, have been used to learn distributions of option prices across strikes, maturities and time, with the aim to produce synthetic data that is the same statistically as real market data. Work by Wiese et al. (2019), building on deep hedging, proposed a GAN-based simulator for equity options. The approach of Wiese et al. (2019) involved a generator network producing a sequence of market states (underlying price and entire implied volatility surface over time) which a discriminator network tries to distinguish from real historical market data. By training these networks against each other, the generator learns to create highly realistic option price dynamics. Similarly, the model introduced as ’VolGAN’ by Vuletić and Cont (2024) has a GAN-based architecture trained on historical sequences of implied volatility surfaces and underlying prices, capable of simulating realistic joint dynamics of the volatility surface and the underlying while respecting static arbitrage constraints. Cont and Vuletić (2023) also introduce another, tractable, method capable of simulating arbitrage free volatility surfaces using more conventional stochastic and factor dynamics with constraints and calibration steps.

2.4.3. Neural SDEs

Noteable examples of neural SDE martingale models include Cuchiero et al. (2020), using GANs to match risk-neutral probability distributions, and Gierjatowicz et al. (2022), who train the neural SDE in a more traditional way. The thesis by Phytides (2023) demonstrates that neural SDEs, when properly constrained with no-arbitrage rules, can accurately replicate a benchmark process (like Black-Scholes and constant elasticity of variance (CEV)). Another approach to note is the work by DeLise (2021), who treats neural SDEs as universal Itô process approximators. Using standard martingale theory and Girsanov’s Theorem, these trained neural SDEs reliably reproduce the Black-Scholes option price (i.e. he trains the neural SDE on just the underlying geometric Brownian motion dynamics and has it reproduce Black-Scholes option prices). In contrast, Lei Fan (2024) use real world historical S&P 500 and S&P 100 index options to train neural SDE models for pricing and hedging.

Cohen et al. (2023a) developed an arbitrage-free neural SDE market model for European options, incorporating explicit constraints to enforce static arbitrage-free relationships among options. In their framework, a set of latent factors (such as implied volatility levels or slopes) evolve in time along with the underlying asset price, and these jointly drive the entire surface of option prices. The drift and volatility of these factors are given by neural networks, which are trained on time-series data of options and underlying prices. To ensure the learned model does not violate no-arbitrage conditions, the neural network outputs are passed through transformations that enforce the linear inequalities reflecting static arbitrage constraints. This constrained neural network approach guarantees that every price surface produced by the model is free of calendar arbitrage and butterfly arbitrage by design. The result is a nonparametric yet arbitrage-consistent market model that can be calibrated to the market’s entire option book. Such a model can serve as a market simulator, generating realistic paths of future implied volatility surfaces for risk management use, or be used to price illiquid exotic derivatives by Monte Carlo simulation in the learned dynamics.

2.4.4. Physics Informed Neural Networks (PINNs) for Options Pricing

PINNs offer a powerful way to reduce the “black-box” nature of deep learning by embedding known physical or theoretical laws directly into the learning process. In PINNs, the solution to a (partial) differential equation (PDE) is approximated by a neural network ${\mu }_{\theta }(x,t)$, and the governing PDE is enforced by penalising any deviation from its expected behavior.

There have been several attempts at applying PINNs to the option pricing problem. The thesis by Tanios (2021) demonstrates the feasibility of applying PINNs to high-dimensional option pricing problem by solving both the forward and inverse formulations of the Black-Scholes and Heston partial differential equations (PDEs). In his work, PINNs are used not only to accurately approximate the solution to these PDEs, but also to compute sensitivities (Greeks) via automatic differentiation, thereby overcoming the limitations of traditional numerical methods in high-dimensional settings. Similarly, recent work by Dhiman and Hu (2023) extends the PINN framework to both European and American options, showing that embedding financial constraints—such as no-arbitrage conditions—directly into the neural network’s loss function can lead to improved accuracy and computational efficiency. While Tanios uses synthetically generated data, Dhiman and Hu (2023) include the application of PINNs to real market data.

2.5. Choice of Modelling Approach

Static constraints define a feasible state space for option prices, while dynamic constraints govern admissible drifts/diffusions over time. The “liquid lattice”, i.e., linear-inequality parameterisation, enables tractable enforcement of static constraints during learning and cleaning. For our use case, these constraints are not optional regularisation but necessary domain knowledge that stabilises learning and ensures economically meaningful outputs.

Deep hedging (i.e., RL) targets hedging strategy optimisation, not market option surface dynamics. Generative models (e.g., GANs) excel at sample realism, but can struggle with enforcing the absence of arbitrage and with the interpretability of calibrated factor dynamics. PINNs embed the PDE structure, but align best with parametric PDE models. Thus, we will focus on Neural SDEs, which provide a stochastic, non-parametric, time-series view with continuous-time semantics. Combined with explicit no-arbitrage conditions, they can model the joint evolution of the entire option surface while remaining economically consistent. Building on the approach of Cohen et al. (2023a), we will construct a model with dynamics which are learnt directly from the data, for the temporal evolution of the options surface.

3. Data Preprocessing

This section details the pre-processing steps required to take raw options quotes, obtained from Deribit, and clean and transform them so that they are ready to be fed into the neural SDE model. This involves the creation of a static lattice of options at high frequency, arbitrage repair of these option price surfaces and a transformation of the surfaces into a lower dimensional representation.

3.1. Static Grid for Options

Over time, the time-to-maturity/moneyness coordinates of liquid options change stochastically. To facilitate training our model, it is hence required to create a static lattice of time-to-maturity/moneyness coordinates, to which market options are then matched. Let $\tau =T-t$ denote time to maturity. m is forward log moneyness, i.e., $m=ln\left({\displaystyle \frac{K}{F_t\left(T\right)}}\right)$, with K being the strike and $F_t$ being the forward price.

It is important to note that the entire option surface is never observable at once; it is only ever the currently liquid options for which we see quotes within the order book. Because of this, for each time to maturity $\tau$ there is a different set of log-moneynesses m. In particular, options further from expiry have a larger range of traded moneyness values in comparison to options traded closer to maturity. This represents what Cohen et al. (2023a) call the ’liquid lattice’ of options. When constructing our static grid, it is important to capture this aspect, as we are interested in modelling the movement of tradeable options that have sufficient liquidity. Conversely, we are not interested in modelling a very deep out-of-the-money option that has a one day expiry, because we are very unlikely, given most market conditions, to observe this option ever being actively traded.

3.1.1. Constructing a Static Grid at High-Frequency

Prior “market model” work in the literature has also parameterised the price of options in this fashion. Wissel (2008) constructs a square lattice (i.e. each maturity has the same distribution of moneyness), while Cohen et al. (2023a) use end-of-day closing data for STOXX 50 index and DAX options, which naturally provide a distribution of moneyness for each maturity and hence the shape of the liquid lattice is chosen heuristically.

However, at high frequency/intraday frequency each maturity has effectively a continuous distribution of traded moneynesses. Figure 5 highlights the observed distribution of available log-moneyness values for each maturity cluster in the Deribit data for a 24-hour period on May 22nd 2025.

This leads to a method that captures the distribution of m for each $\tau$, or the shape of the “liquid lattice”. We are primarily interested in ensuring that this liquid lattice reflects the liquidity profile within our observed dataset. We construct a lattice over the joint $(\tau ,m)$ space as follows. First, let ${\left\{{\overline{\tau }}_i\right\}}_{i=1}^p$ be the sorted unique maturities observed in the dataset. For each ${\overline{\tau }}_i$, compute empirical moneyness percentiles $m_{\mathrm{lo},i}=P_1(m\mid \tau ={\overline{\tau }}_i)$ and $m_{\mathrm{hi},i}=P_{99}(m\mid \tau ={\overline{\tau }}_i)$, and form an $N_m$-point linear grid

$$m_{ij}=m_{\mathrm{lo},i}+(j-1){\displaystyle \frac{m_{\mathrm{hi},i}-m_{\mathrm{lo},i}}{N_m-1}},j=1,\dots ,N_m.$$

The lattice nodes are the Cartesian set $\mathcal{N}=\left\{({\overline{\tau }}_i,m_{ij})\right\}$. We then fit a 1-nearest-neighbor map in the 2D space (Euclidean metric) and assign each quote k with features $({\tau }_k,m_k)$ directly to its nearest node index ${\ell }_k=arg{min}_{({\overline{\tau }}_i,m_{ij})\in \mathcal{N}}\parallel ({\tau }_k,m_k)-({\overline{\tau }}_i,m_{ij})\parallel$. For each calendar time t and node in $\mathcal{N}$, if multiple quotes map to the same node, we keep the one with the highest liquidity (measured by trading volume). Finally, we pivot to a time-by-node panel and fill any missing entries by linear interpolation over t, followed by forward/backward fill as needed.

This method is implemented by Algorithm 1, noting that within the 24 hours of data, we take all unique maturities rather than performing clustering on $\tau$. Let $\mathcal{Q}={\left\{q_k\right\}}_{k=1}^M$ be quote observations, where each

$$q_k=(t_k,T_k,K_k,B_k,A_k,S_k,V_k,{\mathrm{type}}_k),$$

with timestamp $t_k$, expiry $T_k$, strike $K_k$, best bid $B_k$, best ask $A_k$, underlying spot $S_k$, liquidity proxy $V_k$ (USD traded volume), and option type. Fix reference date $t_{\mathrm{ref}}$, rates $(r,q)$, and lattice hyperparameters $n_{\tau }\in \mathbb{N}$, $n_m\in \mathbb{N}$, and $K_{max}\in \mathbb{N}$ (default 50).

For each retained quote k, define

$${\tau }_k:={\displaystyle \frac{T_k-t_{\mathrm{ref}}}{365.25}},F_k:=S_k{e}^{(r-q){\tau }_k},m_k:=log\left({\displaystyle \frac{K_k}{F_k}}\right),$$

$${\mathrm{mid}}_k:={\displaystyle \frac{B_k+A_k}{2}},c_k:={\mathrm{mid}}_k{e}^{-(r-q){\tau }_k}(=C_k/F_k).$$

Retain calls only (${\mathrm{type}}_k=\mathrm{C}$) with $V_k>0$.

Algorithm 1: Static Lattice Construction

Require:  Cleaned quotes ${\left\{(t_k,{\tau }_k,m_k,c_k,V_k)\right\}}_{k=1}^M$, $n_{\tau }$, $n_m$, $K_{max}$, fill method (linear) Ensure:  Reduced node set ${\mathcal{L}}_{\mathrm{sub}}$, dense panel $C(t,\ell )$ on ${\mathcal{L}}_{\mathrm{sub}}$ 1: $\mathcal{U}\leftarrow \mathrm{sort}\left(\mathrm{unique}{\left\{{\tau }_k\right\}}_{k=1}^M\right)$ 2: Fit 1D K-means with $n_{\tau }$ clusters on $\mathcal{U}$; obtain centroids ${\left\{{\tilde{\tau }}_i\right\}}_{i=1}^{n_{\tau }}$ 3: Sort centroids increasingly: ${\overline{\tau }}_1<\dots <{\overline{\tau }}_{n_{\tau }}$ and relabel quote cluster indices accordingly 4: for$i=1,\dots ,n_{\tau }$do 5:     ${\mathcal{M}}_i\leftarrow \{m_k:\mathrm{quote}k\mathrm{assigned}\mathrm{to}\mathrm{cluster}i\}$ 6:     if ${\mathcal{M}}_i=⌀$ then 7:         $({\underline{m}}_i,{\overline{m}}_i)\leftarrow \left(P_1\left(\left\{m_k\right\}\right),P_{99}\left(\left\{m_k\right\}\right)\right)$ 8:     else 9:         $({\underline{m}}_i,{\overline{m}}_i)\leftarrow \left(P_1\left({\mathcal{M}}_i\right),P_{99}\left({\mathcal{M}}_i\right)\right)$ 10:     end if 11:     for $j=1,\dots ,n_m$ do 12:         $m_{ij}\leftarrow {\underline{m}}_i+{\displaystyle \frac{j-1}{n_m-1}}({\overline{m}}_i-{\underline{m}}_i)$ 13:         Create node $\ell =(i,j)$ with coordinates $x_{\ell }:=({\overline{\tau }}_i,m_{ij})$ 14:     end for 15: end for 16: $\mathcal{L}\leftarrow \left\{x_{\ell }\right\}$, $N:=\left|\mathcal{L}\right|=n_{\tau }{n}_m$ 17: Fit 1-nearest-neighbor model on $\mathcal{L}$ in ${\mathbb{R}}^2$ (Euclidean distance) 18: for  $k=1,\dots ,M$do 19:     ${\ell }_k\leftarrow arg{min}_{\ell \in \{1,\dots ,N\}}{\parallel ({\tau }_k,m_k)-x_{\ell }\parallel }_2$ 20: end for 21: for each pair $(t,\ell )$ observed among $\left\{(t_k,{\ell }_k)\right\}$ do 22:     ${\mathcal{I}}_{t,\ell }\leftarrow \{k:t_k=t,{\ell }_k=\ell \}$ 23:     $k^{★}(t,\ell )\leftarrow arg{max}_{k\in {\mathcal{I}}_{t,\ell }}{V}_k$ 24:     $C^{\mathrm{sparse}}(t,\ell )\leftarrow c_{k^{★}(t,\ell )}$ 25: end for 26: Set all unfilled entries of $C^{\mathrm{sparse}}$ to NaN 27: Node coverage count: $n_{\ell }\leftarrow \#\{t:C^{\mathrm{sparse}}(t,\ell )\mathrm{is}\mathrm{observed}\}$ 28: Keep ${\mathcal{L}}_K$ = indices of the $K_{max}$ largest $n_{\ell }$ 29: Restrict $C^{\mathrm{sparse}}$ to columns in ${\mathcal{L}}_K$; drop columns entirely NaN 30: Let remaining column indices be ${\mathcal{L}}_{\mathrm{sub}}$ and corresponding nodes be ${\left\{x_{\ell }\right\}}_{\ell \in {\mathcal{L}}_{\mathrm{sub}}}$ 31: for each $\ell \in {\mathcal{L}}_{\mathrm{sub}}$ do 32:     Interpolate $C^{\mathrm{sparse}}(·,\ell )$ linearly in time 33:     Apply forward-fill then backward-fill at time boundaries 34: end for 35: Denote resulting dense panel by $C(t,\ell )$ 36: return ${\mathcal{L}}_{\mathrm{sub}}$, $C(t,\ell )$

If train/test splitting is used, build $(\mathcal{L},\mathrm{NN})$ on the training quotes only, apply the snapping-and-panel step separately to train and test, then keep the intersection of node indices:

$${\mathcal{L}}_{\cap }={\mathcal{L}}_{\mathrm{sub}}^{\mathrm{train}}\cap {\mathcal{L}}_{\mathrm{sub}}^{\mathrm{test}},$$

and reindex both panels to ${\mathcal{L}}_{\cap }$ with the same linear + forward/backward fill.

The resulting liquid lattice from this process is displayed in Figure 6, which is the lattice from 24 hours of data from the 22nd of May 2025. Here we can see that, while the lattice is irregular and asymmetric, it does contain the general shape of an options lattice we would expect, with options closer to maturity having a smaller distribution support of traded moneyness values. Additionally, one can speculate that running this process over a longer time-period, say one month, the lattice would look more symmetric as one would observe more traded combinations of $(\tau ,m)$.

3.2. Arbitrage Repair of Historical Option Prices

When training an options model, we want to ensure that the input data is arbitrage free. Thus cleaning input data of arbitrage violations is a logical next step in the data processing. We may encounter seeming arbitrage violations in the input data for two reasons:

1.

Using the mid-price to construct our model. Within our data, and in the microstructure of markets, we do not see a single price for any asset, but rather two quotes: the bid, which is the highest price a buyer is willing to pay, and the ask, which is the lowest price a seller is willing to accept. The information contained within the spread of these relates to the liquidity of the option’s price and is captured within the lattice creation, via volume weighting, and the arbitrage cleaning process. The existence of an arbitrage-free surface between the spread is discussed in Appendix A.

2.

Data artefacts of our lattice creation method. The lattice of $(\tau ,m)$ pairs that are created using Algorithm 1 are not guaranteed to satisfy our static arbitrage constraints. Even if the mid-price captured a static arbitrage-free surface, we have interpolated these prices to common nodes which could potentially break convexity and/or monotonicity conditions.

Thus it is important to note that what may seem like arbitrage violations in the data are unlikely to be opportunities to make actual arbitrage profits in the market.

From Algorithm 1, we have a timeseries of option surfaces. We can use methods from Cohen et al. (2020) to “repair” these constructed surfaces of static arbitrage violations. Note that Cohen et al. (2020) applied this method to end-of-day closing data, while in our context it is applied to intraday, high frequency cryptocurrency options data.

Recall the representation of static arbitrage constraints reproduced from Cohen et al. (2020) in Section 2.2 above:

$$Ac\ge 0$$

where A represents a matrix encoding each of the options shape constraints and c is a stacked vector of call option prices, across moneyness/time-to-maturity pairs. Let $A\in {\mathbb{R}}^{m\times N}$ and $b\in {\mathbb{R}}^m$ encode the shape constraints (following Cohen et al. (2020)):

$$\begin{array}{cc}\hfill \left(\mathrm{i}\right)\mathrm{Monotone}\mathrm{in}\mathrm{maturity}\text{:}& c({\tau }_{i+1},m)-c({\tau }_i,m)\ge 0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \left(\mathrm{ii}\right)\mathrm{Monotone}\mathrm{in}\mathrm{strike}\text{:}& c(\tau ,m_j)-c(\tau ,m_{j+1})\ge 0,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \left(\mathrm{iii}\right)\mathrm{Convex}\mathrm{in}\mathrm{strike}\text{:}& {\displaystyle \frac{c(\tau ,m_j)-c(\tau ,m_{j-1})}{x_j-x_{j-1}}}\le {\displaystyle \frac{c(\tau ,m_{j+1})-c(\tau ,m_j)}{x_{j+1}-x_j}},x_j:=e^{m_j},\hfill \\ \hfill \left(\mathrm{iv}\right)\mathrm{Bounds}\mathrm{and}\mathrm{intrinsic}& \hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \mathrm{lower}\mathrm{envelope}\text{:}& 0\le c(\tau ,m)\le 1,c(\tau ,m)\ge {\left(1-e^m\right)}_{+}.\hfill \end{array}$$

(4)

On a non-uniform $x=e^m$ grid, () is implemented as a linear inequality

$${\displaystyle \frac{1}{x_j-x_{j-1}}}{c}_{j-1}-\left({\displaystyle \frac{1}{x_j-x_{j-1}}}+{\displaystyle \frac{1}{x_{j+1}-x_j}}\right)c_j+{\displaystyle \frac{1}{x_{j+1}-x_j}}{c}_{j+1}\ge 0.$$

When adjacent $\tau$ rows have misaligned m grids, (1) evaluates $c({\tau }_{i+1},m)$ by linear interpolation on the ${\tau }_{i+1}$ row.

To speed up repeated projections, redundant rows of $Ac\ge b$ are eliminated by LP tests, a method which comes from Caron et al. (1989): For each row i, solve

$$\underset{c}{max}{A}_{i}c\mathrm{s}.\mathrm{t}.A_{-i}c\ge b_{-i}.$$

where $A_{-i}$ and $b_{-b}$ refer to the matrix A with its i-th row removed and the vector b with its i-th element removed. We are trying to test if the i-th constraint is redundant, so if the optimum $\le b_i+\mathrm{tol}$, row i is redundant and removed.

Subsequently, each time slice is repaired by the nearest feasible point

$$\begin{array}{c}\hfill {\widehat{c}}_t=arg\underset{c\in {\mathbb{R}}^N}{min}\frac{1}{2}{\parallel c-c_t^{\mathrm{raw}}\parallel }_2^2\mathrm{s}.\mathrm{t}.Ac\ge b,0\le c\le 1.\end{array}$$

(5)

This strictly enforces all static constraints while minimally perturbing the raw panel. When bid/ask data is available, we perform a spread-aware L¹ projection, biasing adjustments to remain within the spread by introducing non-negative up/down slack ${ϵ}^{+},{ϵ}^{-}$ and weights $w^{+},w^{-}$ proportional to inverse half-spreads:

$$\begin{array}{cc}\hfill \underset{c,{ϵ}^{\pm }\ge 0}{min}& \sum _{n=1}^N\left(w_n^{+}{ϵ}_n^{+}+w_n^{-}{ϵ}_n^{-}\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& c=c_t^{\mathrm{mid}}+{ϵ}^{+}-{ϵ}^{-},Ac\ge b,0\le c\le 1,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & w_n^{+}={\displaystyle \frac{1}{max\{c_{t,n}^{\mathrm{ask}}-c_{t,n}^{\mathrm{mid}},\delta \}}},w_n^{-}={\displaystyle \frac{1}{max\{c_{t,n}^{\mathrm{mid}}-c_{t,n}^{\mathrm{bid}},\delta \}}}.\hfill \end{array}$$

(8)

The spread-aware L¹ adjustment matches the method of Cohen et al. (2020) exactly.

After these steps, the repaired surfaces ${\left\{{\widehat{c}}_t\right\}}_t$ are static-arbitrage free and feed all downstream steps. Because each time slice solves a convex program with warm starts, the method scales to high-frequency data while remaining faithful to the observed bid-ask data. Figure 7 highlights an example of the adjustments made to create an arbitrage-free surface. We can see that the deep out-of-the-money options are the ones that have the largest corrections. This can make sense, as these options are typically the ones with the largest bid-ask spread and also the ones with the most unreliable (and potentially stale) quotes.

3.3. Linear Factor Representation

Now that we have a series of arbitrage free options surfaces, we notice that it is high-dimensional and that there are likely redundancies in information contained within the high-dimensional representation, in the sense that we can expect prices on the surface to be correlated in some way. As such it is useful and efficient to decompose the surface into a low-dimensional representation.

We follow the method presented in Cohen et al. (2023a), where it is labelled Algorithm 1. This algorithm is designed to decompose the surface into latent factors that are designed to minimise arbitrage in the reconstruction.

Suppose that the surface can be described by d latent factors ${\xi }_{1t},\dots ,{\xi }_{dt}$, collected into the factor vector ${\Xi }_t={({\xi }_{1t},\dots ,{\xi }_{dt})}^{\top }\in {\mathbb{R}}^d$. We link these factors to the vector of call option prices via a linear decomposition:

$$\begin{array}{c}\hfill \tilde{c}=G_0(\tau ,m)+\sum _{i=1}^d{G}_i(\tau ,m){\xi }_{it}\end{array}$$

(9)

Here $G_0(\tau ,m)\in {\mathbb{R}}^N$ is a baseline price vector (the temporal mean of the arbitrage-repaired surfaces) and each $G_i(\tau ,m)\in {\mathbb{R}}^N$, $i=1,\dots ,d$, is a vector of linear loading coefficients indexed over the N lattice nodes $(\tau ,m)$, mapping the scalar factor ${\xi }_{it}$ to price perturbations across the surface. Stacking the loading vectors column-wise gives the $N\times d$ loading matrix $M=[G_1\mid G_2\mid \dots \mid G_d]$, so that (9) can be written compactly as $\tilde{c}=G_0+M{\Xi }_t$.

The loading vectors $G_i$ and associated factor scores ${\xi }_{it}$ are obtained via a modified principal components analysis (PCA) procedure, following Algorithm 1 of Cohen et al. (2023a). The decomposition is performed sequentially, with each successive factor chosen to optimise a different objective:

1.

Minimal Dynamic Arbitrage — The first loading vector $G_1$ (denoted $G_{\mathrm{dyn}}$) is chosen to align with the leading principal component of the empirical z-factor, i.e., the deterministic drift operator applied to the reconstructed surface (see Appendix B). This ensures the dominant direction of surface evolution is representable by the factor model, thereby minimising violation of the HJM-type drift restriction.

2.

Statistical Accuracy — The second loading vector $G_2$ (denoted $G_{\mathrm{stat}}$) is the leading principal component of the residual price variation after projecting out $G_1$. This captures the direction of maximum remaining variance, ensuring accurate reconstruction.

3.

Minimal Static Arbitrage — The third loading vector $G_3$ (denoted $G_{\mathrm{sa}}$) is chosen from the residual subspace to minimise the static arbitrage violations (measured via the constraint matrix A) in the reconstructed surfaces. A hinge-penalty objective penalises violations of $A\tilde{c}\ge 0$ during the optimisation of this direction.

The loading matrix is therefore partitioned as $M=[G_{\mathrm{dyn}}\mid G_{\mathrm{stat}}\mid G_{\mathrm{sa}}]$ under a three-factor model ($d=3$).

Figure 8 illustrates the three loading vectors and the baseline $G_0$ under a three-factor model. While the stochastic loading vectors are difficult to interpret individually, the baseline $G_0$ exhibits the standard call surface shape one would expect. Similar to the approach by Cohen et al. (2023a), all loading vectors $G_i$ ($i\ge 1$) are hard coded to zero at maturity $\tau =0$ to ensure that the terminal payoff condition is satisfied.

4. Neural SDE Modelling

This section builds on the method of Cohen et al. (2023a), applying it to high-frequency cryptocurrency options data to obtain a data-driven, arbitrage-free model taking in the timeseries of available cryptocurrency option market data.

4.1. Stochastic Dynamics of Market Factors

We first posit an SDE for the evolution of our market factors under the statistical (a.k.a. “objective”, “physical” or “real–world”) probability measure, commonly denoted as the $\mathbb{P}$ measure. Specifically, the dynamics of the d-dimensional factor vector ${\Xi }_t={({\xi }_{1t},\dots ,{\xi }_{dt})}^{\top }$ introduced in Section 3.3 are assumed to be governed by the SDE,

$$d{\Xi }_t={\mu }_{\theta }({\Xi }_t,t)dt+{\sigma }_{\varphi }({\Xi }_t,t)d{W}_t$$

(10)

where ${\mu }_{\theta }$ is the drift vector and ${\sigma }_{\varphi }$ is the diagonal diffusion matrix and $W_t$ is a standard Brownian motion. Each of ${\mu }_{\theta }$, ${\sigma }_{\varphi }$ and $W_t$ are d-dimensional, where d refers to the number of factors chosen from our factor decoding, as described in Section 3.3. Our diffusion matrix is assumed diagonal, meaning we assume that the noise terms between our $\Xi$ factors are uncorrelated, and with the 24 hours of data this assumption seems to hold relatively well when analysing the residual factor plots, see Section 4.4.2.

In our data-driven approach, we represent the drift and diffusion components of this SDE with neural networks with weights $\theta$ and $\varphi$ respectively. The specific architecture of the neural networks is detailed in Appendix E.

We can formulate the loss function for training the neural networks in this model by considering the Euler-Maryuama discretisation for a stochastic differential equation. Under this expansion, we have a defined transition. Let $\Delta t_k:=t_{k+1}-t_k$ and $\Delta W_k:=W_{t_{k+1}}-W_{t_k}\sim \mathcal{N}(0,\Delta t_k{I}_d)$. The Euler–Maruyama step for the Itô SDE is

$$\begin{array}{cc}\hfill {\Xi }_{k+1}& ={\Xi }_k+{\mu }_{\theta }({\Xi }_k,t_k)\Delta t_k+{\sigma }_{\varphi }({\Xi }_k,t_k)\Delta W_k\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & ={\Xi }_k+{\mu }_{\theta }({\Xi }_k,t_k)\Delta t_k+{\sigma }_{\varphi }({\Xi }_k,t_k)\sqrt{\Delta t_k}{\epsilon }_k,{\epsilon }_k\sim \mathcal{N}(0,I_d).\hfill \end{array}$$

(12)

Conditionally on ${\Xi }_k$, the transition density is Gaussian:

$$\begin{array}{c}{\Xi }_{k+1}\mid {\Xi }_k\sim \mathcal{N}\left(m_k,\Delta t_k{\Gamma }_k\right),\hfill \\ \hfill m_k:={\Xi }_k+{\mu }_{\theta }({\Xi }_k,t_k)\Delta t_k,{\Gamma }_k:={\sigma }_{\varphi }({\Xi }_k,t_k){\sigma }_{\varphi }{({\Xi }_k,t_k)}^{\top }.\end{array}$$

(13)

We can then formulate an expectation maximisation (EM) loss function. Define the observed increment $\Delta {\Xi }_k:={\Xi }_{k+1}-{\Xi }_k$. The one-step negative log-likelihood (NLL) is

$$\begin{array}{c}{\ell }_k(\theta ,\varphi )=\frac{1}{2}[{\left(\Delta {\Xi }_k-{\mu }_{\theta }({\Xi }_k,t_k)\Delta t_k\right)}^{\top }{\left(\Delta t_k{\Gamma }_k\right)}^{-1}\left(\Delta {\Xi }_k-{\mu }_{\theta }({\Xi }_k,t_k)\Delta t_k\right)\hfill \\ \hfill +logdet\left(\Delta t_k{\Gamma }_k\right)+dlog\left(2\pi \right)],\end{array}$$

(14)

so that for a path ${\left\{{\Xi }_k\right\}}_{k=0}^K$, the total loss function value is

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{NLL}}(\theta ,\varphi )=\sum _{k=0}^{K-1}{\ell }_k(\theta ,\varphi ).\end{array}$$

(15)

This NLL function now formulates our test loss, where minimising the NLL is equivalent to maximising the likelihood of seeing the next step given the current data. This follows the classical EM framework from Dempster et al. (1977) and is the neural SDE learning strategy from Dridi et al. (2021). This is structurally analogous to the maximisation/inference loop approach from Greff et al. (2017).

This training scheme does make the assumption that the transition/discretisation of the process is Gaussian and at ultra-high-frequency (sub-second level intervals) this assumption does not hold. When using 24 hours of high-frequency cryptocurrency options data, the prices were aggregated, by taking the median price, to minutely increments in order for goodness of fit tests to look reasonable. The training results at high frequency and associated discussion can be found in Appendix C.

One can also consider using different expansions or quasi-likelihoods as the loss term for the SDE. Aït-Sahalia (2002) develops a closed-form approximation for the likelihood of discretely observed diffusion processes, addressing the fact that most continuous-time SDEs lack an analytically tractable transition density. The approach constructs a sequence of explicit Hermite polynomial expansions that converge uniformly to the true but unknown transition density. By successively correcting the leading Gaussian term through higher-order terms in the sampling interval $\Delta t$, this method yields a quasi-maximum likelihood estimator that remains asymptotically equivalent to the true MLE. Conceptually, it extends the Euler–Maruyama Gaussian transition by incorporating higher-order corrections that depend on the local drift, diffusion, and their spatial derivatives, thereby capturing the nonlinearity and skewness inherent in finite-time transitions. We explored using this approximation as a higher-order likelihood correction to the Gaussian transition in (13), attempting to capture the non-Gaussianity of the high frequency intervals, and the training results in this case are presented and discussed in Appendix D.

4.2. Arbitrage Free Factor Representation

In order to constrain the outputs of our model to be free of static arbitrage, we must first consider what these constraints mean in our factor space. As our constraints are linear inequalities and our factor representation of options prices is also linear, these constraints have a representation in the factor space. Recall the linear constraints in Section 2.2, $Ac\ge 0$ and the linear factor decomposition, $\tilde{c}=G_0(\tau ,m)+{\sum }_{i=1}^d{G}_i(\tau ,m){\xi }_{it}$. Substituting this decomposition into the constraint system gives

$$A\tilde{c}=A{G}_0+\sum _{i=1}^{d}A{G}_i{\xi }_{it}\ge 0.$$

Hence, in the factor space the arbitrage-free region is characterised by all factor vectors ${\Xi }_t={({\xi }_{1t},\dots ,{\xi }_{dt})}^{\top }$ satisfying

$$A_G{\Xi }_t\ge b_G,A_G:=[A{G}_1,A{G}_2,\dots ,A{G}_d],b_G:=-A{G}_0.$$

This defines a convex polytope, or region of space, in the factor space, within which any combination of factors yields option prices that respect all static arbitrage constraints. During training or simulation, we therefore enforce that the decoded factors ${\Xi }_t$ remain inside this feasible region

$${\mathcal{P}}_{\mathrm{arb}}=\left\{{\Xi }_t|A_G{\Xi }_t\ge b_G\right\},$$

ensuring that the reconstructed surface is statically arbitrage-free.

Removing redundant linear constraints as described in Section 3.2, we can plot our static arbitrage free convex polytope and the associated $\xi$ factors to visually inspect if our factor decoding process in Section 3.3 has been done correctly. Figure 9 presents scores in a two factor model and the corresponding static no-arbitrage polygon. Inspecting this image we can see that the vast majority of our latent $\xi$ factors lie within the acceptable polygon. Those few scores, which lie either on the boundary or just outside the polygon, are there due to the hinge-penalty-based optimisation done to reduce static arbitrage in option price reconstruction in the factor decoding process given in Section 3.3.

4.3. Constrained Networks and No-Static-Arbitrage Guarantees

Given the specified region of space in which the $\Xi$ factors can evolve without option prices allowing static arbitrage, we can enforce constraints on the Neural SDE to ensure that the process remains within the required polytope. Cohen et al. (2023a) provide the method to achieve this using Friedman and Pinsky conditions for the non-attainability of a boundary for an SDE (Friedman and Pinsky 1973).

Definition: [Friedman–Pinsky Non-Attainability Conditions]

Consider an Itô diffusion in ${\mathbb{R}}^d$,

$$d{\Xi }_t=\mu \left({\Xi }_t\right)dt+\sigma \left({\Xi }_t\right)d{W}_t,$$

and let $\mathcal{P}=\{\Xi \in {\mathbb{R}}^d:a_i^{\top }\Xi \ge b_i,i=1,\dots ,n_f\}$ be a convex polytope defined by $n_f$ linear constraints (i.e. $n_f$ faces). For each face i with boundary function

$$q_i(\Xi ):=a_i^{\top }\Xi -b_i=0,$$

define the normal drift and diffusion components near the boundary as

$${\mu }_i(\Xi ):=a_i^{\top }\mu (\Xi ),{\sigma }_i^2(\Xi ):=a_i^{\top }\Gamma (\Xi )a_i,\Gamma (\Xi ):=\sigma (\Xi )\sigma {(\Xi )}^{\top }.$$

Then, a sufficient condition for the boundary $\partial \mathcal{P}$ to be non-attainable (i.e. the process remains in $\mathcal{P}$ almost surely if started inside) is that, for each i,

$$\begin{array}{cc}\hfill & \left(\mathrm{i}\right){\mu }_i(\Xi )\ge 0\mathrm{whenever}{q}_i(\Xi )=0,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & \left(\mathrm{ii}\right){\sigma }_i^2(\Xi )=0\mathrm{whenever}{q}_i(\Xi )=0.\hfill \end{array}$$

(17)

That is, the drift of the SDE points inward at the boundary and the diffusion vanishes in the direction normal to it.

Cohen et al. (2023a) enforce of these conditions via a smooth Lipschitz function. This function will shrink the diffusion as it gets closer to the boundary and slowly add drift in the opposite direction. Choosing such a function to be smooth and Lipschitz is important so as not to break backpropagation for training the neural networks. Recall the linear price map $\tilde{c}=G_0+M{\Xi }_t$ from (9), where ${\Xi }_t\in {\mathbb{R}}^d$ is the factor vector and $M=[G_{\mathrm{dyn}}\mid G_{\mathrm{stat}}\mid G_{\mathrm{sa}}]\in {\mathbb{R}}^{N\times d}$ is the loading matrix. Given the lattice constraints $A\tilde{c}\ge b$, we map them into factor space as

$$H{\Xi }_t\ge h,H:=AM,h:=b-A{G}_0,$$

and augmented with a data-driven box $L\le {\Xi }_t\le U$. The box bounds $L,U\in {\mathbb{R}}^d$ are obtained component-wise from the training factor scores as $L_k=Q_{\alpha }\left({\left\{{\xi }_{k,t}\right\}}_t\right)$ and $U_k=Q_{1-\alpha }\left({\left\{{\xi }_{k,t}\right\}}_t\right)$ for a small quantile level $\alpha$ (defaulting to componentwise min/max when quantile trimming is not applied). These box constraints are converted into additional half-space rows appended to $H{\Xi }_t\ge h$, so that the diffusion shrinkage and inward drift correction treat empirical range boundaries identically to the static-arbitrage faces. Let $n_f$ denote the total number of half-space constraints (arbitrage faces plus box faces) after this augmentation. Rows of H are row-normalised when building the shrink geometry, so the signed distance of ${\Xi }_t$ to face i is

$${\rho }_i\left({\Xi }_t\right):={\displaystyle \frac{H_i{\Xi }_t-h_i}{\parallel H_i\parallel }}=H_i{\Xi }_t-h_i(\sin \mathrm{ce}\parallel H_i\parallel =1).$$

Following, Cohen et al. (2023a) the Lipschitz shrink profile is given by

$$\varphi \left(x\right)={\displaystyle \frac{x_{+}}{1+x_{+}}},x_{+}:=max\{x,0\},$$

Near the boundary ${\rho }_i\left({\Xi }_i\right)\downarrow 0$, we have $\varphi \left({\rho }_i\right)\downarrow 0$; far inside the polytope $\varphi \left({\rho }_i\right)\uparrow 1$.

At a given state ${\Xi }_t$ we select the $k=min\{d,n_f\}$ most relevant faces (smallest ${\rho }_i$), build an orthonormal frame $Q\left({\Xi }_t\right)\in {\mathbb{R}}^{d\times d}$ whose first k columns point along the corresponding unit normals, and compute shrink coefficients $\epsilon \left({\Xi }_t\right)\in {[0,1]}^k$ with components ${\epsilon }_j\left({\Xi }_t\right)=\varphi \left({\rho }_{i_j}\left({\Xi }_t\right)\right)$. This yields two equivalent shrink operators (chosen by mode):

- Matrix shrink:

$$P\left({\Xi }_t\right)=Q\left({\Xi }_t\right)\left[\begin{array}{cc}\mathrm{diag}\left(\epsilon \left({\Xi }_t\right)\right)& 0\\ 0& I_{d-k}\end{array}\right]Q{\left({\Xi }_t\right)}^{\top }\in {\mathbb{R}}^{d\times d},$$

so the effective diffusion is

$${\Gamma }_{\mathrm{eff}}\left({\Xi }_t\right)=P\left({\Xi }_t\right)\Gamma \left({\Xi }_t\right)P{\left({\Xi }_t\right)}^{\top }.$$

- Diagonal shrink (default, diagonal noise): compute elementwise multipliers $s\left({\Xi }_t\right)\in {[0,1]}^d$ so that $\mathrm{diag}\left(s\left({\Xi }_t\right)\right)\sigma \left({\Xi }_t\right)$ best matches $P\left({\Xi }_t\right)\mathrm{diag}\left(\sigma \left({\Xi }_t\right)\right)$ in a covariance sense. The effective diffusion used inside the model is

$${\sigma }_{\mathrm{eff}}\left({\Xi }_t\right)=s\left({\Xi }_t\right)\odot \sigma \left({\Xi }_t\right),$$

In both cases, along any active face normal $a_i$ one has $a_i^{\top }{\Gamma }_{\mathrm{eff}}\left({\Xi }_t\right)a_i\to 0$ as ${\rho }_i\left({\Xi }_t\right)\downarrow 0$, satisfying (). The same shrink can also be applied at the loss level: rather than modifying the diffusion $\sigma$ directly, the constraint geometry is absorbed into a state-dependent preconditioner $\Omega \left({\Xi }_t\right)\in {\mathbb{R}}^{d\times d}$ (derived from $P\left({\Xi }_t\right)$) that transforms the observed increments and likelihood. Specifically, the tuple $(\Omega ,det\Omega ,{\Omega }^{-\top }\Delta \Xi )$ is passed to the negative log-likelihood (14), so that the Gaussian transition is evaluated in the preconditioned coordinate system where the boundary constraints are automatically respected.

To enforce (16) via an inward drift correction, we precompute, along the observed factor path $X_t$, per–time inward directions and normalised boundary distances using

$$({\mathrm{dirs}}_t,{\epsilon }_t^{\mu })=(W,b,X_{·};{\rho }_{★},{\epsilon }_{★}^{\mu }),$$

where $(W,b)$ are the row–normalised halfspaces from the polytope output. At training time, for a batch indexed by t, the raw drift $\mu (·)$ is corrected by

$${\mu }_{\mathrm{corr}}\left(X_t\right)=\mu \left(X_t\right)+\Delta {\mu }_t,\Delta {\mu }_t=\sum _k{w}_{t,k}{\mathrm{dirs}}_{t,k},w_{t,k}={\displaystyle \frac{{\lambda }_{\mathrm{bc}}}{{\epsilon }_{t,k}^{\mu }+{\epsilon }_0}}\mathbf{1}\{{\epsilon }_{t,k}^{\mu }\le {\tau }_{\mathrm{cut}}\},$$

with user-selectable parameters $({\lambda }_{\mathrm{bc}},{\epsilon }_0,{\tau }_{\mathrm{cut}})$ and an optional norm cap $\parallel \Delta {\mu }_t\parallel \le \mathrm{cap}$.

Together, the smooth diffusion shrink guarantees $a_i^{\top }{\Gamma }_{\mathrm{eff}}\left({\Xi }_t\right)a_i\to 0$ as the boundary is approached, and the inward drift correction ensures $a_i^{\top }{\mu }_{\mathrm{corr}}\left({\Xi }_t\right)\ge 0$ there; hence the Friedman–Pinsky non-attainability conditions are satisfied for the learned SDE on the static-arbitrage polytope.

4.4. Training Results

This section will discuss and review training results and goodness of fit tests of this Neural SDE model on our Deribit cryptocurrency options data.

4.4.1. Loss Curve

Illustrated in Figure 10, each loss curve, with and without the drift/diffusion constraints, suggest convergence. However, the loss with the constraints is far lower than the loss without them. Additionally, the gap between train and test loss is smaller. This is interpreted as suggesting that the constraints are indeed working as expected, reducing the state space of our $\Xi$ factor values.

Here, the train/test split is 80/20. Meaning, the first 80% of the 24-hour period is used as training and the remaining 20% is the test, or out-of-sample, region.

4.4.2. Residuals Plots

Following the estimation of the Neural SDE, we assess the adequacy of the fitted drift and diffusion functions by examining the one-step standardised residuals. Under the Euler–Maruyama discretisation of the stochastic differential equation

$$d{Y}_t={\mu }_{\theta }\left(Y_t\right)dt+{\sigma }_{\varphi }\left(Y_t\right)d{W}_t,$$

the discrete-time transition between observations can be approximated as

$$Y_{t+\Delta t}\approx Y_t+{\mu }_{\theta }\left(Y_t\right)\Delta t+{\sigma }_{\varphi }\left(Y_t\right)\sqrt{\Delta t}{\epsilon }_t,{\epsilon }_t\sim \mathcal{N}(0,I_d),$$

where ${\epsilon }_t$ are the model innovations implied by the fitted dynamics. Given observed data ${\left\{Y_t\right\}}_{t=1}^T$, we therefore compute the empirical residuals as

$${\widehat{\epsilon }}_t={\displaystyle \frac{Y_{t+\Delta t}-Y_t-{\mu }_{\theta }\left(Y_t\right)\Delta t}{{\sigma }_{\varphi }\left(Y_t\right)\sqrt{\Delta t}}}.$$

If the estimated SDE correctly captures the conditional mean and variance of the data, then ${\widehat{\epsilon }}_t$ should be approximately independent and standard normal across time and dimensions.

These residuals thus serve as a key diagnostic measure: quantile–quantile (QQ) plots of ${\widehat{\epsilon }}_t$ against a standard normal distribution test for deviations from Gaussianity (such as skewness or excess kurtosis), while pairwise scatter plots of ${\widehat{\epsilon }}_t^{\left(i\right)}$ and ${\widehat{\epsilon }}_t^{\left(j\right)}$ assess the cross-sectional independence across latent factors.

From the QQ Plots in Figure 11, we can see that the assumption of Gaussian transition increments from our Euler-Maryuama training step is approximately correct, though not perfect. There is still some overdispersion at the tails and underdispersion at centrality. This suggests that deviations from our models predictions are still possible, which make sense considering the cryptocurrency market’s higher volatility.

In addition to marginal normality checks, we also compute the pairwise correlations between residual components to evaluate the cross-sectional independence of the model’s innovations. For each pair of factors $(i,j)$, the sample correlation coefficient is given by

$${\rho }_{ij}={\displaystyle \frac{\mathrm{Cov}({\widehat{\epsilon }}^{\left(i\right)},{\widehat{\epsilon }}^{\left(j\right)})}{\sqrt{\mathrm{Var}\left({\widehat{\epsilon }}^{\left(i\right)}\right)\mathrm{Var}\left({\widehat{\epsilon }}^{\left(j\right)}\right)}}}.$$

Under a correctly specified diffusion, the residuals should be uncorrelated across factors, reflecting that each dimension of the learned SDE captures an independent source of risk or market variation.

Analysing the residuals as plotted in Figure 12, we can see that there is some slight correlation between the ’dynamic arbitrage’ factor and the ’statistical accuracy’ factor. In a sense, this is not unexpected, as in testing we notice that the vast majority of variance is explained with a 2-factor model. When training with two factors (one statistical-accuracy factor and one static-arbitrage factor), we observe an $R^2$ value of 99.5–99.8%, indicating that almost all variation in the observed option prices is captured by the model’s reconstructed surfaces. In this context, $R^2$ measures the proportion of total price variance explained by the factor representation, reflecting both the completeness and goodness-of-fit of the decomposition. In a way this is unsurprising, as other martingale stochastic volatility models are also two-factor models (see e.g. Heston (1993); Hagan et al. (2002)). Aside from that, the correlations are low and residuals are center aligned.

5. Applications

This section considers two fundamental downstream applications: hedging and risk management, in the form of Value at Risk (VaR). These represent the most immediate and practically meaningful uses of a market pricing model. Both tasks provide direct validation of the model’s ability to capture market dynamics: hedging evaluates how well the learned dynamics can replicate real-world sensitivities and risk exposures, while VaR assesses the model’s capacity to simulate realistic loss distributions under stochastic evolution. However, this is not to say that this is an exhaustive list of the model’s capabilities as it provides a complete stochastic representation of how the implied volatility and option prices evolve through time. In principle, this framework can be extended to a wide range of quantitative finance tasks — including scenario generation, risk-neutral valuation, stress testing, derivative portfolio simulation, and even data-driven model calibration or discovery of latent market factors.

5.1. Market Simulation and VaR

(Cohen et al. 2023b) detail the methodology on how to create VaR estimates for options portfolios using Neural SDE market models. In essence, the Neural SDE is a generative model, meaning that it can simulate the learned $\Xi$ market factors. Via a factor decoding process (i.e. moving from the low dimensional factor space back to a higher dimensional options surface), simulated prices for various option portfolios can be obtained. As we now have a distribution of potential returns for the options on our defined liquid lattice, we can naturally compute VaR estimates (and hence other risk metrics such as expected short-fall) using this distribution.

5.1.1. Market Simulation Methodology

Let ${\Xi }_t\in {\mathbb{R}}^d$ denote the vector of latent market factors, evolved according to the trained Neural–SDE:

$$d{\Xi }_t={\mu }_{\theta }\left({\Xi }_t\right)dt+{\sigma }_{\varphi }\left({\Xi }_t\right)d{W}_t,{\sigma }_{\varphi }\left({\Xi }_t\right):=diag\left(\sigma \left({\Xi }_t\right)\right)\ge 0,$$

where ${\mu }_{\theta }$ and ${\sigma }_{\varphi }$ represent the learned drift and diffusion networks, and $W_t$ is a standard Brownian motion. During simulation, we again ensure that the factor process ${\Xi }_t$ remains within the static–arbitrage polytope defined by the half-space constraints:

$${\mathcal{P}}_{\mathrm{arb}}=\{{\Xi }_t\mid A_G{\Xi }_t\ge b_G\},$$

For each face i of the polytope, define the signed distance

$${\rho }_i\left({\Xi }_t\right)=A_{G,i}{\Xi }_t-b_{G,i},i=1,\dots ,m,$$

which measures proximity to the ith constraint boundary.

As before, we apply a smooth Lipschitz function ${\phi }_{\sigma }:{\mathbb{R}}_{+}\to [0,1]$ to modulate the effective diffusion near the boundary:

$${\epsilon }_i\left({\Xi }_t\right)={\phi }_{\sigma }\left({\rho }_i\left({\Xi }_t\right)\right)={\displaystyle \frac{{\rho }_i\left({\Xi }_t\right)}{1+{\rho }_i\left({\Xi }_t\right)}}.$$

Let $Q\left({\Xi }_t\right)$ be an orthonormal basis spanning the k most active faces (smallest ${\rho }_i$), and define

$$P\left({\Xi }_t\right)=Q\left({\Xi }_t\right)\left[\begin{array}{cc}diag\left(\epsilon \left({\Xi }_t\right)\right)& 0\\ 0& I_{d-k}\end{array}\right]Q{\left({\Xi }_t\right)}^{\top }.$$

The effective diffusion is then given by

$${\sigma }_{\mathrm{eff}}\left({\Xi }_t\right)=diag\left(s\left({\Xi }_t\right)\odot \sigma \left({\Xi }_t\right)\right),$$

where $s\left({\Xi }_t\right)\in {[0,1]}^d$ are elementwise shrink coefficients matching $P\left({\Xi }_t\right)\sigma \left({\Xi }_t\right)$ in covariance. By construction, the normal component $A_{G,i}^{\top }{\Gamma }_{\mathrm{eff}}\left({\Xi }_t\right)A_{G,i}\to 0$ as ${\rho }_i\left({\Xi }_t\right)\downarrow 0$, satisfying the diffusion non-attainability condition.

Again following our previous procedure, a smooth correction is added to ensure the drift points inward at the boundary. Define a Lipschitz function ${\phi }_{\mu }\left(\rho \right)={\displaystyle \frac{1}{1+\rho }}$ and inward directions $d_i$ associated with each active face normal $A_{G,i}$. The corrected drift is:

$${\mu }_{\mathrm{corr}}\left({\Xi }_t\right)={\mu }_{\theta }\left({\Xi }_t\right)+\sum _{i\in \mathcal{I}\left({\Xi }_t\right)}{\omega }_i\left({\Xi }_t\right)d_i,{\omega }_i\left({\Xi }_t\right)={\displaystyle \frac{{\left[-A_{G,i}^{\top }{\mu }_{\theta }\left({\Xi }_t\right)-{\phi }_{\mu }\left({\rho }_i\left({\Xi }_t\right)\right)\right]}_{+}}{\langle d_i,A_{G,i}\rangle }},$$

ensuring $A_{G,i}^{\top }{\mu }_{\mathrm{corr}}\left({\Xi }_t\right)\ge 0$ along all active faces.

Given a discretisation step $\Delta t$, each Euler-Maryuama step evolves as:

$${\Xi }_{t+\Delta t}={\Xi }_t+{\mu }_{\mathrm{corr}}\left({\Xi }_t\right)\Delta t+{\sigma }_{\mathrm{eff}}\left({\Xi }_t\right)\sqrt{\Delta t}{\epsilon }_t,{\epsilon }_t\sim \mathcal{N}(0,I_d).$$

Running B Monte Carlo paths of ${\Xi }_t$ generates a sample distribution of factor evolutions ${\left\{{\Xi }_t^{\left(b\right)}\right\}}_{b=1}^B$.

5.1.2. VaR of Option Portfolios

Table 1 presents VaR results for various option portfolios, for both the Neural SDE simulation and a standard historical simulation approach.

The results in Table 1 show that across all portfolios, the Neural–SDE simulated Value-at-Risk (VaR) estimates are consistently higher in magnitude than those obtained via historical simulation. This indicates that the Neural–SDE captures fatter tails and greater short-term variability in the factor dynamics, reflecting more realistic risk in high-frequency markets such as cryptocurrency options. The largest VaR differences appear in short-maturity, short-volatility positions (e.g., short strangles and short condors), where non-linear payoffs are most sensitive to volatility clustering and jump-like behaviour. Conversely, for longer-dated or more balanced portfolios (e.g., butterflies or calendar spreads), the Neural–SDE and historical methods yield more comparable results, suggesting that the learned diffusion smooths risk over longer horizons.

Plotting a few examples of loss distributions in Figure 13 and comparing against a distribution from historical simulation, we can see that the Neural SDE does indeed provide a different distribution. We interpret these as reflecting data-driven and arbitrage-consistent dynamics of the underlying surface. Noting that traditional historical simulation approaches assume that options, or market factors, move individually, however we know that options move as a surface with constraints to their state, as discussed in Section 2.2. The historical simulation approach taken here in comparison, is a simulation of the market factors, that is to say that we look at the historical distribution of returns we have (in the test period) and simulate returns by sampling from this distribution. Option prices/portfolios are simulated by a factor decoding process.

5.2. Hedging of Option’s Portfolios

Cohen et al. (2022) highlight how can we use the learned model to hedge a portfolio of options. While their work presents both sensitivity based hedging and minimum-variance (MV) hedging, in our case the applications to cryptocurrency options will be focused on the MV hedge as it directly uses the learned Neural SDE model, whereas the sensitivity based hedge only uses the factor representation of option prices.

5.2.1. Minimum Variance Hedging

Let ${\Xi }_t\in {\mathbb{R}}^d$ denote the latent market factors driving the option surface, and recall the linear decoder from Section 3.3,

$$\tilde{c}\left({\Xi }_t\right)=G_0+\sum _{i=1}^d{G}_i{\xi }_{it}=G_0+M{\Xi }_t,M\in {\mathbb{R}}^{N\times d},$$

(18)

where $M=[G_1\mid \dots \mid G_d]$ is the loading matrix defined in Section 3.3. For a portfolio defined by node weights $w\in {\mathbb{R}}^N$, the portfolio’s factor exposure vector is given by

$$s=M^{\top }w\in {\mathbb{R}}^d,$$

(19)

To first order, the one-step change in portfolio value is given by

$$\Delta V_t\approx s^{\top }\Delta {\Xi }_t,\Delta {\Xi }_t:={\Xi }_{t+1}-{\Xi }_t.$$

(20)

From the learned Neural–SDE model, the drift and diffusion functions ${\mu }_{\theta }\left({\Xi }_t\right)$ and ${\sigma }_{\varphi }\left({\Xi }_t\right)$ define the instantaneous covariance of the latent factors as

$${\Sigma }_t^{\mathrm{base}}={\sigma }_{\varphi }\left({\Xi }_t\right){\sigma }_{\varphi }{\left({\Xi }_t\right)}^{\top }.$$

(21)

To ensure that this covariance structure respects the static arbitrage geometry defined by the polytope constraints $A_G{\Xi }_t\ge b_G$, we apply a state-dependent linear map ${\Omega }_t$ that shrinks the diffusion along the active boundary directions, as detailed in Section 17. The resulting adjusted covariance is

$${\Sigma }_t={\Omega }_t{\Sigma }_t^{\mathrm{base}}{\Omega }_t^{\top }.$$

(22)

This procedure enforces that factor dynamics remain within the feasible arbitrage-free region during simulation and hedging.

Now, consider a set of k hedging instruments indexed by $\mathcal{H}=\{j_1,\dots ,j_k\}$. Let $B\in {\mathbb{R}}^{d\times k}$ denote their factor exposures,

$$B:=M_{:,\mathcal{H}}^{\top }.$$

(23)

The static minimum-variance (MV) hedge solves

$$u^{★}=arg\underset{u\in {\mathbb{R}}^k}{min}Var\left({(s-Bu)}^{\top }\Delta {\Xi }_t\right)={\left(B^{\top }\overline{\Sigma }B\right)}^{-1}{B}^{\top }\overline{\Sigma }s,$$

(24)

where $\overline{\Sigma }={\displaystyle \frac{1}{T-1}}{\sum }_t{\Sigma }_t$ is the time-averaged model-implied factor covariance. A dynamic version of this hedge recomputes $u_t^{★}$ at each time step according to

$$u_t^{★}={\left(B^{\top }{\Sigma }_{t}B\right)}^{-1}{B}^{\top }{\Sigma }_{t}s,$$

(25)

yielding time-varying hedge weights that evolve consistently with the latent factor covariance implied by the Neural–SDE. In practice, the ${\Sigma }_t$ matrices are computed per step using the endogenous shrinkage ${\Omega }_t$ and the trained model parameters $(\theta ,\varphi )$.

5.2.2. Hedging an Example Portfolio (A 30-Day ATM Option)

We present hedging results for a portfolio consisting an at-the-money option with a 30-day expiry. Figure 14 showcases the root mean squared error of the Neural SDE MV hedge compared with standard Black-Scholes delta and delta-vega hedging.

While we cannot say that the MV hedging procedure from the learnt Neural SDE model outperforms the delta or delta-vega hedge, we can state that it is comparable to these traditional hedging methods. These hedging results on cryptocurrency are similar to the results from Cohen et al. (2022), who present hedging results from traditional equities and forex markets.

6. Conclusions

This paper provides a synthesis of several methods in the literature and presented a pipeline that takes in raw options quotes and produces a non-parametric model for the entire options surface. The presented Neural SDE model adheres to static arbitrage constraints and the applications to risk management and hedging have been showcased. One can conclude that combining machine learning with known financial laws, e.g., arbitrage constraints, can provide non-parameteric alternatives to traditional models.

Our results demonstrate the feasibility of applying the arbitrage-free Neural SDE market model of Cohen et al. (2023a) to an entirely new asset class—high-frequency cryptocurrency options—where the microstructure and volatility regimes differ substantially from traditional markets. In particular, this extends the model’s applicability to ultra-high-frequency data, developing preprocessing methods for constructing static liquid lattices and repairing arbitrage in sparse intraday panels. The elements brought together here form a complete pipeline from raw data to an operational, data-driven, arbitrage-free simulator of the Bitcoin options market.

Appendix B.1. Drift of the Option Surface Under the Physical Measure

From the factor representation of normalised call prices:

$${\tilde{c}}_t(\tau ,m)=G_0(\tau ,m)+\sum _{i=1}^d{\xi }_{i,t}{G}_i(\tau ,m),$$

(A9)

where $\tau =T-t$ and $m=log(K/F_t)$, applying Itô’s lemma yields the differential of the call surface evaluated along the moving coordinates $({\tau }_t,m_t)$:

$$d{\tilde{c}}_t=\sum _i{G}_i({\tau }_t,m_t)d{\xi }_{i,t}+{\displaystyle \frac{\partial {\tilde{c}}_t}{\partial \tau }}d{\tau }_t+{\displaystyle \frac{\partial {\tilde{c}}_t}{\partial m}}d{m}_t+{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2{\tilde{c}}_t}{\partial m^2}}d{\langle m\rangle }_t.$$

(A10)

Substituting $d{\tau }_t=-dt$, $d{m}_t=(-{\alpha }_t+\frac{1}{2}{\gamma }_t^2)dt-{\gamma }_{t}d{W}_{0,t}$, and $d{\langle m\rangle }_t={\gamma }_t^{2}dt$ gives the drift of surface evolution:

$${\displaystyle \frac{d{\tilde{c}}_t}{dt}}=-{\displaystyle \frac{\partial {\tilde{c}}_t}{\partial \tau }}-{\displaystyle \frac{1}{2}}{\gamma }_t^2{\displaystyle \frac{\partial {\tilde{c}}_t}{\partial m}}+{\displaystyle \frac{1}{2}}{\gamma }_t^2{\displaystyle \frac{{\partial }^2{\tilde{c}}_t}{\partial m^2}}+\sum _i{\mu }_{i,t}{G}_i({\tau }_t,m_t).$$

(A11)

Appendix B.2. The Infinitesimal Drift Operator and the z-Factor

We define the deterministic differential operator

$$Z:=-{\displaystyle \frac{\partial }{\partial \tau }}-{\displaystyle \frac{1}{2}}{\gamma }_t^2{\displaystyle \frac{\partial }{\partial m}}+{\displaystyle \frac{1}{2}}{\gamma }_t^2{\displaystyle \frac{{\partial }^2}{\partial m^2}},$$

(A12)

and apply it to the reconstructed option surface to obtain the z-factor:

$$z_t=Z\left[{\tilde{c}}_t\right]({\tau }_t,m_t).$$

(A13)

This represents the infinitesimal deterministic drift of the option surface induced by the decay of time to maturity and the diffusion of moneyness. Intuitively, $z_t$ describes how the option surface should evolve in $(\tau ,m)$-space in order to remain dynamically arbitrage-free.

Appendix B.3. Dynamic No-Arbitrage Under the Basis Expansion

Under the basis expansion of Eq. (A11), the HJM drift restriction becomes:

$$G^{\top }{\mu }_t-z_t=G^{\top }{\sigma }_t{\psi }_t,$$

(A14)

where ${\mu }_t$ and ${\sigma }_t$ are respectively the drift and diffusion coefficients of the factor process $d{\xi }_t={\mu }_{t}dt+{\sigma }_{t}d{W}_t$, and ${\psi }_t$ is the vector of market prices of risk. Eq. (A14) states that, up to a change of measure, the observed drift of the factors $G^{\top }{\mu }_t$ should be explainable by the deterministic surface drift $z_t$ and the model’s diffusion structure $G^{\top }{\sigma }_t$.

Appendix B.4. Principal Component Analysis of the z-Fact or

We compute $z_t$ empirically for each observation of the reconstructed surface and perform Principal Component Analysis (PCA) over its temporal evolution. The PCA decomposes the variability of $z_t$ into orthogonal modes corresponding to consistent directions of deterministic surface motion (e.g., level, skew, curvature shifts). If the market is dynamically arbitrage-free, these principal directions should align closely with the factor basis $\left\{G_i\right\}$.

Residual components of $z_t$ that are orthogonal to $\mathrm{span}\left\{G_i\right\}$ indicate systematic drifts of the surface that are not representable by the model’s factors, implying violations of the HJM drift restriction. By quantifying the variance explained by the first d principal components of $z_t$, we obtain a diagnostic of the extent to which the Neural-SDE factor model captures the dynamically consistent evolution of the option surface.

Appendix B.5. Minimising Dynamic Arbitrage

By training the model and selecting factor structures such that the leading principal components of $z_t$ align with the learned factors $G_i$, we effectively minimise the unexplained component of the HJM drift restriction. This ensures that the model-implied surface evolution under the physical measure $\mathbb{P}$ is as close as possible to one that admits an equivalent martingale measure $\mathbb{Q}$, that is, the dynamic arbitrage in the model is minimised.

Appendix B.6. Empirical Computation of the ’z’ Factor

While the work by Cohen et al. (2023a) uses splines to calculate the derivatives in A12, we have access to high-frequency data so we use finite difference methods. The ${\gamma }_t$ parameter is computed by fitting an initial, ’Stage 0’, neural SDE to just the underlying process.

Figure A1. Surface derivatives at t=0, 22nd of May 2025 Deribit call options.

Figure A1. Surface derivatives at t=0, 22nd of May 2025 Deribit call options.

Our estimated surface derivatives are jagged/irregular due to the sparsity of our lattice. However, they do approximately resemble the surface derivatives presented by Cohen et al. (2023a) in shape.