Operator Deep Smoothing for Implied Volatility

## Abstract

Abstract We devise a novel method for nowcasting implied volatility based on neural operators.
Better known as implied volatility smoothing in the financial industry, nowcasting of implied volatility means constructing a smooth surface that is consistent with the prices presently observed on a given option market.
Option price data arises highly dynamically in ever-changing spatial configurations, which poses a major limitation to foundational machine learning approaches using classical neural networks.
While large models in language and image processing deliver breakthrough results on vast corpora of raw data, in financial engineering the generalization from big historical datasets has been hindered by the need for considerable data pre-processing.
In particular, implied volatility smoothing has remained an instance-by-instance, hands-on process both for neural network-based and traditional parametric strategies.
Our general operator deep smoothing approach, instead, directly maps observed data to smoothed surfaces.
We adapt the graph neural operator architecture to do so with high accuracy on ten years of raw intraday S&P 500 options data, using a single model instance.
The trained operator adheres to critical no-arbitrage constraints and is robust with respect to subsampling of inputs (occurring in practice in the context of outlier removal).
We provide extensive historical benchmarks and showcase the generalization capability of our approach in a comparison with classical neural networks and SVI, an industry standard parametrization for implied volatility.
The operator deep smoothing approach thus opens up the use of neural networks on large historical datasets in financial engineering.

## 1 Introduction

Options trading experienced phenomenal growth in recent years.
In its 2023 trading volume report (Cboe Global Markets, Inc.,, ),
the CBOE announced the fourth consecutive year of record-breaking volumes on its options exchanges, citing a record-breaking number of transactions for European options on the S&P 500 index.
European options are financial derivative contracts that give their holder the right, but not the obligation, to either buy or sell an underlying asset at a predetermined price (the *strike*) at a predetermined time (the *expiry*).
An option specifying the right to buy (respectively to sell) is called a *Call* (respectively *Put*) option.
Options are traded on a wide range of underlyings, including stocks, indices, currencies and commodities, and can be used to hedge against or speculate on the price movements of the underlying asset.

A key concept in options trading is the so-called *implied volatility*,
which transforms the nominal price of an option into a conceptually and numerically convenient metric.
The *implied volatility surface* is the collection of implied volatilities as observed at a specific point in time, visualized in three-dimensional space as a surface over the (strike, expiry)-domain.
It provides an intuitive representation of the current state of the options market and is crucial for hedging and risk management.
The extraction of a smooth surface from quoted option prices is called *implied volatility smoothing* and allows to infer (or *nowcast*) theoretical option prices for interpolated strike values and expiry times.
It remains one of the key challenges in options trading.

Conventionally, implied volatility smoothing relies on parametric surfaces whose parameters are optimized based on the distance to observed prices while adhering to *absence-of-arbitrage* conditions,
which ensure the consistency of prices extrapolated from the smoothed surface.
The development of such ad-hoc models for implied volatility traces back to *SVI* (Gatheral,, ),
which models implied volatility slice-wise for each maturity and successfully captures its key features on Equity indices.
A continuous interpolation scheme for SVI slices was provided in Gatheral and Jacquier, (), yielding a full surface.
Nowadays, sophisticated market makers employ custom parametrizations, which can be considered proprietary trading secrets and reduce SVI to a benchmark role.

Regardless of the particular parametric surface model used, the conventional smoothing approach boils down to the continued execution of a numerical optimization routine:
A smoothed surface expires as soon as quotes are updated (whenever markets move), necessitating the re-calibration of parameters.
Success and duration of this routine is sensitive to initial conditions, search heuristics, and termination criteria, which exposes practitioners to considerable process uncertainties during trading hours (or *online*).
In response, we introduce a novel *operator deep smoothing* approach, replacing the instance-by-instance optimization with a single evaluation of a neural network.
This greatly simplifies online calibration, at the upfront cost of training the network *offline* from historical data (in the spirit of Hernandez, (); Horvath et al., (); Liu et al., ()).
Our unique use of *neural operators* (Kovachki et al.,, ) is fundamentally directed by the observation that the raw inputs for volatility smoothing (the collections of observed volatilities) over time vary in size and spatial arrangement:
Options expire, new maturities and strikes become available, and the coordinates of existing options evolve continuously in the domain of the implied volatility surface (Figure 1(b)).
This setting excludes classical neural networks – which required fixed-size inputs – from direct application.
Neural operators, instead, conceptualize observed data as point-wise discretizations of latent functions in implicit infinite-dimensional function spaces and are well suited for the task.

#### Contributions

We introduce operator deep smoothing, a general approach for discretization-invariant data interpolation based on neural operators, and apply it to implied volatility smoothing.
Our technique transcends traditional parametric smoothing and directly maps observed volatilities to smoothed surfaces.
Comparable neural network based approaches are limited to certain option markets (e.g. FX markets, where options by default spread out on fixed rectilinear grids, as relied upon by Bergeron et al., () for its VAE approach) or require data pre-processing (as in Cont and Vuletić, (), which achieves fixed rectilinear grids by linear interpolation of market values, setting aside questions related to no-arbitrage constraints).
Instead, our technique novelly adapts the graph neural operator (GNO) architecture (Anandkumar et al.,, ) to consistently smooth input data of any size and spatial arrangement.
While neural operators have successfully been used in Physics to numerically solve partial differential equations, our application is the first in financial engineering and highlights the values of their discretization-invariance properties, so far rather under-explored.
We employ our method on ten years of intraday S&P 500 options data, smoothing more than 60 million volatility datapoints using a single model instance with around 100 thousand trained parameters.
We report errors that substantially improve on the SVI industry benchmark and are highly competitive with Ackerer et al., (), which performs smoothing by training one classical neural network per volatility surface.
We proceed to successfully demonstrate the generalization capabilities of our model for end-of-day options data of the S&P 500 as well as three further major US indices.
No data from these three indices has been used for training.

We explore the technical implications of our method in Section 5.
Here, we discuss the broader impact of our contributions:

- •
Operator Deep Smoothing for Implied Volatility – Our method massively simplifies volatility smoothing, an area where effective methods mean competitive advantage and frequently remain trade secrets.
Therefore, we believe that our approach lowers the entry-barrier for effective volatility smoothing even among industry professionals.
Practitioners and researchers that are not directly involved in options trading frequently employ rudimentary methods (SVI or linear/spline interpolation).
Here, our trained operator, served as a hands-off tool, could provide “cheap” and accurate surfaces for use in downstream tasks.
Ultimately, our method may be useful for all participants of option markets.
This includes the general public, whose trading in such markets has been increasing substantially (Doherty et al.,, ) and which benefits from broadly accessible investment tools.
- •
Neural Operators for Discretization-Invariant Interpolation – Our operator deep smoothing approach constitutes the first application of neural operators for interpolation/extrapolation tasks and paves the way to future research on the versatility of the discretization-invariance of neural operators in industrial applications characterized by dynamic and spatially irregular data.
At least in financial engineering, surfaces similar to the implied volatility surface (and higher-dimensional equivalents, as for example the *volatility cube*) are ubiquitous. We expect our technique to be transferable and to streamline and robustify engineers’ and researchers’ algorithms and data pipelines.

#### Literature Review

The aforementioned SVI was developed for internal use at Merrill Lynch in 1999 and later advocated in Gatheral, ().
Its extension to surface-based SSVI in Gatheral and Jacquier, () has been eagerly adopted by practitioners, which have since contributed to its robust calibration and generalizations (Corbetta et al.,, ; Hendriks and Martini,, ; Guo et al.,, ).
It was augmented in Ackerer et al., () by a multiplicative neural network corrector, based on guided network training by means of no-arbitrage soft constraints from Zheng, ().
The absence-of-arbitrage conditions for implied volatility surfaces – providing safeguards for option pricing – were formulated in Roper, (),
and we provide an equivalent formulation, based on Fukasawa, (); Lucic, (), for practical purposes.
In Chataigner et al., () static arbitrage constraints were used to perform option calibration (with an additional regularization technique),
which can be considered to be instance-by-instance smoothing of nominal price data. In
Bergeron et al., () a classical VAE (variational autoencoder) was applied to implied volatility smoothing on FX markets, where strikes of quoted options are tied to a fixed grid of *deltas*.(^1^11Delta is the derivative of the price with respect to the underlying asset and is standard in FX strike quoting.)
This specificity of FX markets allows the use of a conventional feedforward neural network based decoder.
Recent option calibration approaches based on neural networks have been proposed in Baschetti et al., (); Hernandez, (); Horvath et al., (); Van Mieghem et al., ().

A comprehensive account on neural operators is given in Kovachki et al., (),
unifying previous research on different neural operator architectures and techniques (Anandkumar et al.,, ; Li et al.,, ).
Subsequent developments investigating the expressivity of these architectures as well as their generalizations include Hao et al., (); Huang et al., (); Lanthaler et al., (); Li et al., (); Lingsch et al., (); Tran et al., ().

#### Outline

We review financial concepts
and the challenges of implied volatility smoothing in Section 2.
In Section 3, we provide a review of neural operators (Section 3.1)
and introduce our operator deep smoothing approach
for general interpolation tasks (Section 3.2).
In Section 4, we perform experiments for implied volatility smoothing of S&P 500 options data.
Finally, Section 5
gathers limitations as well as outlooks regarding the use of neural operators for interpolation purposes.

#### Code

We make code for the paper available at the location [https://github.com/rwicl/operator-deep-smoothing-for-implied-volatility](https://github.com/rwicl/operator-deep-smoothing-for-implied-volatility).
In particular, the code repository contains a general *PyTorch* (Paszke et al.,, ) implementation of the graph neural operator architecture for operator deep smoothing.

## 2 Background: Implied Volatility

We consider a market of European options written on an underlying asset, which we observe at a given instant $T_{0}$, and denote the time-$T$ forward price of the underlying asset by $F_{T_{0},T}$.

#### European Call Options

The option market consists of a finite collection of *European Call options*,(^2^22In practice, market participants trade both Call and Put options, which are mathematically equivalent through the well-known *Put-Call parity*. The latter thus allows to speak in terms of Call options only.)
each identified by its expiry $T\in(T_{0},\infty)$ and its strike $K\in(0,\infty)$,
and we write $C(T,K)$ for its (undiscounted) price.
In practice, these are traded for fixed expiries $T_{1},\ldots,T_{m}$;
for each $T_{i}$, only a finite range of strikes $K_{1}^{i},\ldots,K_{n_{i}}^{i}$ is available,
typically widening as the expiry increases (Figure 1).

#### Black-Scholes

The *Black-Scholes model* is the simplest diffusive asset model and captures the volatility of the underlying asset with a single parameter $v\in(0,\infty)$.
Its popularity stems from the closed-form expression it admits for the price of a European Call option with *time-to-expiry* $\tau=T-T_{0}$ and *log-moneyness* $k=\log(K/F_{T_{0},T})$:

$$ $\mathrm{BS}(\tau,k,v)=\Phi\left(d_{1}\left(\tau,k,v\right)\right)-e^{k}\,\Phi \left(d_{2}\left(\tau,k,v\right)\right)$ $$

(in units of the time-$T$ forward of the underlying).
Here, $\Phi$ denotes the cumulative distribution function of the standard Normal distribution, while

$$ $d_{1}(\tau,k,v)=\frac{-k}{v\sqrt{\tau}}+\frac{1}{2}v\sqrt{\tau},\quad d_{2}( \tau,k,v)=\frac{-k}{v\sqrt{\tau}}-\frac{1}{2}v\sqrt{\tau}.$ (1) $$

#### Implied Volatility

While not able (any longer) to fit market data, the mathematical tractability of the Black-Scholes model gave rise to the concept of implied volatility:
Given a Call option with price $C(T,K)$, its implied volatility $v(\tau,k)$ is defined by $C(T,K)=F_{T_{0},T}\operatorname{BS}(\tau,k,v(\tau,k))$.
By using time-to-expiry/log-moneyness coordinates, implied volatility provides a universal way to consistently compare the relative expensiveness of options of across different strikes, expiries, underlyings, and interest rate environments.
Its characteristic shape helps traders make intuitive sense of the instantaneous state of option markets relative to a (flat) Black-Scholes model baseline.

#### Implied Volatility Smoothing

This refers to fitting a smooth surface $\hat{v}\colon(0,\infty)\times\mathbb{R}\to(0,\infty)$ to a collection $\mathbf{v}=\{v(\tau_{l},k_{l})\}_{l=1}^{p}$ of observed implied volatilities.
Naive strategies such as cubic interpolation are ill-fated:
Surfaces generated by such interpolation rules will in general correspond to Call option prices that are exploitable by so-called *arbitrage*, namely cost-less trading strategies generating a guaranteed profit.
In option markets, an arbitrage is called *static* when set up solely from fixed positions in options and a dynamically (but a finite number of times) readjusted position in the underlying.
Beyond simple interpolations, practitioners have devised ad-hoc parametrizations for implied volatility, in particular the aforementioned SVI, which are not expected to perfectly match all reference prices.
Instead, the model parameters are optimized with respect to an objective function that measures market price discrepancy and includes penalization terms ruling out static arbitrage.
These penalization terms are commonly formulated on the basis of the following theorem, which summarizes the shape constraints of the implied volatility surface (Gatheral and Jacquier,, ; Lucic,, ; Roper,, ).

###### Theorem 2.1 (Volatility Validation) .

Let $\hat{v}\colon(0,\infty)\times\mathbb{R}\to(0,\infty)$ be continuous and satisfying

- (i)
Strike arbitrage: For every $\tau>0$, the slice $\hat{v}_{\tau}=\hat{v}(\tau,\makebox[4.30554pt]{{$\cdot$}})$ is of class $C^{2}$ with
$\operatorname{But}(\tau,\makebox[4.30554pt]{{$\cdot$}},\hat{v}_{\tau},\partial
_{k}\hat{v}_{\tau},\partial^{2}_{k}\hat{v}_{\tau})\geq 0$
(2)
and $\limsup_{k\uparrow\infty}\frac{\hat{v}^{2}_{\tau}(k)}{k}<\frac{2}{\tau}$,
where
$\operatorname{But}(\tau,k,v_{0},v_{1},v_{2})=\left(1+d_{1}(\tau,k,v_{0})v_{1}
\sqrt{\tau}\right)\left(1+d_{2}(\tau,k,v_{0})v_{1}\sqrt{\tau}\right)+v_{0}v_{2
}\tau.$
- (ii)
Calendar arbitrage: For each $k\in\mathbb{R}$, $\hat{v}(\cdot,k)\sqrt{\cdot}$ is non-decreasing and vanishes at the origin.

Then, $(T,K)\mapsto\operatorname{BS}(\tau,k,\hat{v}(\tau,k))$ defines a Call price surface that is free from static arbitrage.

Condition 2.1(ii) is equivalent to prices increasing in maturity (uncertain increases as time passes), while Condition 2.1(i) arises when computing the *implied probability density* $f_{\tau}$ of the underlying:

$$ $f_{\tau}(\makebox[4.30554pt]{{$\cdot$}})=\frac{\varphi(-d_{2}(\tau,\makebox[4. 30554pt]{{$\cdot$}},v))}{v}\operatorname{But}(\tau,\makebox[4.30554pt]{{$\cdot $}},v,\partial_{k}v,\partial_{k}^{2}v).$ (3) $$

Since a density needs to be non-negative, (3) explains why Condition 2.1(i) above is required.

Figure: (a) Relative trading volume per quoted interval, averaged over S&P 500 dataset 2012-2021. Left: Rectangular domain w.r.t. time-to-expiry/log-moneyness (99.9% of trading volume with maturities ¡1 year). Right: Rectangular domain w.r.t. transformed coordinates (96.6% of trading volume while more evenly populated); boundary delineated in Left.
Refer to caption: extracted/5912674/resources/trade_volume.png

###### Theorem 2.1 (Volatility Validation) .

## 3 Neural Operators for Discretization-Invariant Smoothing

### 3.1 Background: Neural Operators

We provide full details about notations, terms, and additional context in Appendix A.

#### Philosophy

The development of *neural operators* is based on the philosophy that observed data $\mathbf{a}=\{a_{l}\}_{l=1}^{p}$ arises as the evaluation of a latent function $a\colon D\to\mathbb{R}^{c_{\text{in}}}$, defined on some domain $D\subseteq\mathbb{R}^{d}$, at a discretization $\pi=\{x_{l}\}_{l=1}^{p}$ of $D$.
That is, $\mathbf{a}=a|_{\pi}$, or

$$ $a_{l}=a(x_{l}),\quad l=1,\dots,p.$ (4) $$

An input-output relationship of data $\mathbf{a}\mapsto\mathbf{u}$ is then ”really” described by an operator $F\colon\mathcal{A}\to\mathcal{U}$ between function spaces $\mathcal{A}$ and $\mathcal{U}$.
Neural operators are abstract neural network architectures $F^{\theta}\colon\mathcal{A}\to\mathcal{U}$, with implementations that integrate (4) consistently across the variable discretization $\pi$.

#### Technicalities

We review the core concepts of neural operators from Kovachki et al., ().
Let $D$ be a bounded domain in $\mathbb{R}^{d}$
and $\mathcal{A}$ and $\mathcal{U}$ be Banach spaces of functions mapping from $D$ to $\mathbb{R}^{c_{\text{in}}}$ and $\mathbb{R}^{c_{\text{out}}}$, respectively.
Neural operators are finitely parametrized mappings $F^{\theta}\colon\mathcal{A}\to\mathcal{U}$ with *universality* for continuous target operators and with *discretization-invariant* implementations $\tilde{F}^{\theta}$.
In the space $C(\mathcal{A},\mathcal{U})$ topologized by uniform convergence on compacts, the architecture $F^{\theta}$ is called *universal* if $\{F^{\theta}\}_{\theta\in\Theta}$ is dense in $C(\mathcal{A},\mathcal{U})$, with $\Theta$ the parameter set.
An implementation of $F^{\theta}$ is an algorithm $\tilde{F}^{\theta}$ which accepts observed data $\mathbf{a}=a|_{\pi}$ and outputs a function $u\in\mathcal{U}$ and is such that
$\tilde{F}^{\theta}_{\pi}(\makebox[4.30554pt]{{$\cdot$}})=\tilde{F}^{\theta}(
\makebox[4.30554pt]{{$\cdot$}}|_{\pi})\in C(\mathcal{A},\mathcal{U})$.
Now, $\tilde{F}^{\theta}$ is called *discretization-invariant* if $\lim_{n\uparrow\infty}\tilde{F}^{\theta}_{\pi^{(n)}}=F^{\theta}$ in $C(\mathcal{A},\mathcal{U})$,
given a *discrete refinement*(^3^33A *discrete refinement* of $D$ is a nested sequence $(\pi^{(n)})_{n\in\mathbb{N}}$ of discretizations of $D$ for which for every $\varepsilon>0$, there exists $N\in\mathbb{N}$ such that $\{\mathbb{B}_{\mathbb{R}^{d}}(x,\varepsilon)\colon x\in\pi^{(N)}\}$ covers $D$.) of $D$.

Let $K$ be a set of input functions, compact in $\mathcal{A}$, and let $\varepsilon>0$.
In combination, universality and discretization-invariance allow to posit the existence of parameters $\theta$, such that for all $a\in K$,

$$ $\|\tilde{F}^{\theta}(a|_{\pi})-F(a)\|_{\mathcal{U}}\leq\varepsilon,$ (5) $$

irrespective of the particular discretization $\pi$ given that the data $\mathbf{a}=a|_{\pi}$ is scattered sufficiently densely across $D$.
The training of neural operators is analogous to the classical finite-dimensional setting.
It happens in the context of an implicit *training distribution* $\mu$ on the input space $\mathcal{A}$ and aims at minimizing the generalization error

$$ $R_{\mu}\colon\theta\mapsto\mathbb{E}_{a\sim\mu}\|\tilde{F}^{\theta}(a|_{\pi})- F(a)\|_{\mathcal{U}},$ (6) $$

through the use of gradient descent methods applied to empirical estimates of (6).
These estimates are constructed from a training dataset $\mathcal{D}=\{(\mathbf{a}^{(i)},\mathbf{u}^{(i)})\}_{i=1}^{n}$ of *features* $\mathbf{a}^{(i)}=a^{(i)}|_{\pi^{(i)}}$ and *labels* $\mathbf{u}^{(i)}=F(a^{(i)})|_{\pi^{(i)}}$ on the basis of (mini) batching heuristics and are frequently transformed or augmented by additional terms through the use of apposite loss functions.

### 3.2 Operator Deep Smoothing

Figure: (a) Smoothed surface with input volatility quotes (blue). $\delta_{\text{abs}}=0.004$, $\delta_{\text{spr}}(\hat{v},\mathbf{v})=1.061$.
Refer to caption: extracted/5912674/resources/op_deep_smoothing_surface.png

Let $\mathbf{v}=\{v(x_{l})\}_{l=1}^{p}$ be the collection of observed data,
for example implied volatilities as in Section 2.
This notation silently adopts the neural operator philosophy,
connecting the data point $v(x_{l})$ with coordinates $x_{l}$, hinting at a latent function $v\colon D\to\mathbb{R}$ giving rise to the observed values.
The *smoothing* or *interpolation* task consists of constructing an appropriately regular candidate $\hat{v}$ for $v$ from the data $\mathbf{v}=v|_{\pi}$, with $\pi=\{x_{1},\dots,x_{p}\}$.
The operator deep smoothing approach uses a neural operator $\tilde{F}^{\theta}$, trained using historical data, to generate $\hat{v}:=\tilde{F}^{\theta}(\mathbf{v})$.
It fundamentally leverages the discretization-invariance to produce consistent results even if the ”sensors” $x_{l}$ of the latent function $v$ continuously change in availability and/or location in the domain $D$.
This is the situation for volatility smoothing, where the sensors $x_{l}=(\tau_{l},k_{l})$ are the (time-to-expiry, log-moneyness) coordinates of the quoted options,
which – as noted in Section 1 and illustrated in Figure 1(b) –
move with the market,
thus motivating the operator deep smoothing approach.
We now describe the methodology in more detail. For a first glance at the results we refer to Figure 2.

#### Methodology

Translating the smoothing task to an operator learning problem, we find that the target operator is the (continuous) identity operator $F=F_{\text{id}}\colon\mathcal{A}\ni v\mapsto v\in\mathcal{U}$,
with equal input and output spaces $\mathcal{A}=\mathcal{U}=C(\overline{D})$.
The training dataset is the collection $\mathcal{D}=\{\mathbf{v}^{(i)}\}$ of historical data (labels and features coincide),
and we train a suitable neural operator architecture $\tilde{F}^{\theta}$ such that

$$ $|\tilde{F}^{\theta}(\mathbf{v})(x)-v(x)|\leq\varepsilon,\qquad\text{for all }x \in D,$ (7) $$

for some given error tolerance $\varepsilon$ for the interpolation task,
while $v$ is the latent function of which we observe $\mathbf{v}=v|_{\pi}$.
Instead of minimizing empirical estimates of the $\infty$-norm $\|\tilde{F}^{\theta}(\mathbf{v})-v\|_{\mathcal{U}}$, however, we suggest a *fitting loss* based on the *root mean square relative error*,

$$ $\mathcal{L}_{\text{fit}}(\theta,\mathbf{v}):=\sqrt{\frac{1}{|\pi|}\sum_{x\in \pi}\bigg{(}\frac{|\tilde{F}^{\theta}(\mathbf{v})(x)-v(x)|}{1\vee|v(x)|}\bigg{ )}^{2}},$ (8) $$

for its smoothness and invariance to the scale of the data.(^4^44Empirical estimates of $\infty$-norms and $L^{2}$-norms are equivalent loss functions on finite-dimensional spaces.)
Additional engineering techniques, such as sub-sampling of inputs during training, are explored in our practical investigation in Section 4.

#### Practical Constraints

Depending on the application,
the smoothing task may be subject to constraints.
For volatility smoothing, the smoothed surface $\hat{v}^{\theta}=\tilde{F}^{\theta}(\mathbf{v})$ must be free of static arbitrage.
This is effectively enforced by augmenting the loss function with additional penalization terms, moving away from a pure operator learning problem.(^5^55It is not the goal of volatility smoothing to learn the identity operator throughout the entire input space.)
This does not only promote the relevant properties in the neural operator output but can also help define it when faced with sparsity of data in the domain $D$
(in this context see also Li et al., ()).
From Theorem 2.1,
these penalization terms naturally motivate our choices of $\mathcal{L}_{\text{but}}$ and $\mathcal{L}_{\text{cal}}$.
The strike arbitrage constraint 2.1(i) is handled via

$$ $\mathcal{L}_{\text{but}}(\theta;\mathbf{v})=\left\|\left(\operatorname{But}( \makebox[4.30554pt]{{$\cdot$}},\hat{v}^{\theta},\partial_{k}\hat{v}^{\theta}, \partial_{k}^{2}\hat{v}^{\theta})-\varepsilon\right)^{-}\right\|_{1},$ (9) $$

where we ignore the asymptotic condition since our experiments are focused on the bounded domain $D$
(Figure 1(a)).
The inclusion of $\varepsilon$ promotes strictly positive implied densities (we will use $\varepsilon=10^{-3}$), while we choose the $1$-norm to induce sparsity in the constraint violation.
The calendar arbitrage constraint 2.1(ii) can be tackled analogously with

$$ $\mathcal{L}_{\text{cal}}(\theta;\mathbf{v})=\left\|\left(\partial_{\tau}\left[ (\tau,k)\mapsto v^{\theta}(\tau,k)\sqrt{\tau}\right]-\varepsilon\right)^{-} \right\|_{1},$ (10) $$

where again we ignore the asymptotic condition since $D$ is bounded away from zero time-to-expiry.

#### Interpolating Graph Neural Operator

Various neural operator architectures exist, mostly arising from the *kernel integral transform* framework of Kovachki et al., ().
Most prominently, these include Fourier neural operators (FNO), delivering state-of-the-art results on fixed grid data, as well as graph neural operators (GNO), able to handle arbitrary mesh geometries (both reviewed in Appendix A.1).
While highly effective with documented universality, these neural operators are not directly applicable for interpolation tasks as their layers include a pointwise-applied linear transformation, which limits the output to the set of the input data locations.
Dropping this *local* transformation results in an architecture proved to retain universality Kovachki et al., () and that – at least for its implementation as a GNO – allows to interpolate functions.
On the other hand, Lanthaler et al., () proves universality for the architecture combining the local linear transformation with a simple averaging operation, suggesting the fundamental importance of the collaboration of local and non-local components for the expressivity of neural operators.
This was noted in Kovachki et al., (), for whom retaining the local components can be “beneficial in practice”, and confirmed in our experiments, where a purely non-local architecture led to substantially reduced performance.

We therefore propose a new architecture for operator deep smoothing leveraging GNOs’ unique ability to handle irregular mesh geometries.
We use a purely non-local first layer (dropping the pointwise linear transformation), and use it to produce hidden state at all required output locations, enabling subsequent layers to retain their local transformations.
Since GNOs do not theoretically guarantee a smooth output,
we augment the training with additional regularization terms such as
$\mathcal{L}_{\text{reg}}(\theta;\mathbf{v})=\|\Delta\hat{v}^{\theta}\|_{2}$,
with $\Delta$ the Laplace operator,
and provide a full description in Appendix B.

## 4 Experiments

We detail our practical investigation of the operator deep smoothing approach for implied volatility.

### 4.1 Model Training

#### Dataset and Splits

We perform our numerical experiments using 20-minute intervals of CBOE S&P 500 Index Option data from 2012 to 2021.
The dataset amounts to a collection of 49089 implied volatility surfaces and just above 60 million individual volatility datapoints (after domain truncation).
We refer the reader to Appendix C.1 for full
details on the preparation of the dataset.
We allocate the first nine years of data (2012 to 2020) to training, keeping 750 randomly drawn surfaces for validation purposes,
and use the final year of the dataset (2021) for testing.
This yields a training dataset $\mathcal{D}_{\text{train}}$ containing $n_{\text{train}}=43442$ surfaces, a validation dataset $\mathcal{D}_{\text{val}}$ containing $n_{\text{val}}=750$ surfaces and a test dataset $\mathcal{D}_{\text{test}}$ with $n_{\text{test}}=4897$ surfaces.

#### Data Transformation

Motivated by Figure 1(a), we transform time-to-expiry and log-moneyness via $\rho=\sqrt{\tau}$ and $z=k/\rho$.
Intuitively, this transformation converts the“natural” scaling of implied volatility by the square root of time-to-expiry to a scaling of the input domain.
From here on, we consider the domain in these coordinates, setting $D=(\rho_{\text{min}},\rho_{\text{max}})\times(z_{\text{min}},z_{\text{max}})=(
0.01,1)\times(-1.5,0.5)$.
In $(\tau,k)$-coordinates, $D$ becomes a cone-shaped region, that, on average, contains 96.6% of traded options (with time-to-expiry below one year) and, with respect to $(\rho,z)$-coordinates, is more evenly populated, improving the numerics.

#### Model Configuration

We remind that we rely on the interpolating graph neural operator introduced in Section 3.2 and described in detail in Appendix B.
The model hyperparameters (giving rise to 102529 trainable parameters in total) were identified by manual experimentation and are detailed in Appendix C.2.

#### Loss Function

We implement a *Vega*-weighted (see Appendix C.3) version of the fitting loss $\mathcal{L}_{\text{fit}}$ from (8).
We compute $\mathcal{L}_{\text{but}}$ directly as (9)
(in the transformed coordinates) on a synthetic grid using finite differences.
For $\mathcal{L}_{\text{cal}}$, we implement a multiplicative version that is invariant to the level of implied volatility.
We provide a precise description in Appendix C.3.

#### Model Training

We train the GNO over $500$ epochs on the training dataset using the AdamW optimizer with learning rate $\lambda=10^{-4}$ and weight decay rate $\beta=10^{-5}$,
and use a pseudo batch size of 64 by accumulating gradients.
We randomly sub-sample the inputs $\mathbf{v}$
and randomize the grids on which we compute the arbitrage losses.
The training is performed in around 200 hours using a NVIDIA Quadro RTX 6000 GPU.
The validation loss is reported in Appendix C.3 in Table 2.
Code and trained model weights are provided in the supplementary material.

### 4.2 Results

Figure: Figure 3: Benchmark metrics (surface-averages, between 25%- and 75%-quantiles) over training period (computed from $\mathcal{D}_{\text{val}}$) and testing period, as well as surface-averaged spread. Resampled, monthly for training period and weekly for testing period.
Refer to caption: extracted/5912674/resources/finetune/benchmark_quantiles.png

#### Benchmark Metrics

Let $\mathbf{v}=\{v(x)\}_{x\in\pi}$ be the collection of observed implied volatilties and $\hat{v}$ the smoothed surface as produced by a given method.
We measure absolute relative error:

$$ $\delta_{\text{abs}}(\hat{v}(x),v(x))=\frac{|\hat{v}(x)-v(x)|}{v(x)}.$ $$

The average of $\delta_{\text{abs}}(\hat{v},\mathbf{v})$ is known as the *mean absolute percentage error* (or *MAPE*).
As in Corbetta et al., (), we moreover realize the importance of analyzing the smoothing algorithm in terms of nominal price error relative to the size of the bid-ask spread $s(x)=\operatorname{BS}_{\pm}(x,v_{\text{ask}}(x))-\operatorname{BS}_{\pm}(x,v_
{\text{bid}}(x))$.(^6^66$v_{\text{bid}}(x)$ and $v_{\text{ask}}(x)$ are the implied volatilities corresponding to Bid and Ask option prices while $\operatorname{BS}_{\pm}$ is the Black-Scholes formula for Call (resp. Put) options for positive (resp. negative) log-moneyness values: $\operatorname{BS}_{\pm}(\tau,k,v)=\begin{cases}\Phi(d_{1}(\tau,k,v))-e^{k}\Phi (d_{2}(\tau,k,v)),&k>0\\ e^{k}\Phi(-d_{2}(\tau,k,v))-\Phi(-d_{1}(\tau,k,v)),&k\leq 0\end{cases}.$)
We define

$$ $\delta_{\text{spr}}(\hat{v}(x),v(x))=\frac{2}{s(x)}\big{|}\operatorname{BS}_{ \pm}(x,\hat{v}(x))-\operatorname{BS}_{\pm}(x,v(x))\big{|}.$ $$

Since we use simple mid reference prices, $\delta_{\text{spr}}(\hat{v}(x),v(x))\leq 1$ indicates that the prediction $\hat{v}(x)$ for the option $x$ lies within the bid-ask spread.

#### Evaluation and Model Finetuning

During production use, the GNO would be retrained regularly using the most recent available data.
We emulate this procedure during our evaluation of the benchmark metrics over the test dataset $\mathcal{D}_{\text{test}}$ (containing the year 2021):
Following the evaluation of the first month’s test data, the GNO is trained for 10 epochs on the test data just reviewed, with each mini-batch augmented by an equal amount of data from the training dataset $\mathcal{D}_{\text{train}}$.
We repeat this, progressively incorporating an additional month of data, until the entire dataset $\mathcal{D}_{\text{test}}$ is assessed.
This finetuning-evaluation procedure takes circa 1.8 GPU hours per month.

#### Analysis

It is apparent from Figure 3 that our operator deep smoothing approach substantially improves on SVI’s smoothing capabilities, with respect to both $\delta_{\text{abs}}$ and $\delta_{\text{spr}}$.(^7^77We produce the SVI benchmark as described in Section E.)
Our approach, with monthly finetuning, smooths the volatility surface with a MAPE of around 0.5%, while SVI fluctuates between 1% and 2%.
The various figures in Appendix D illustrate the qualitative improvements of our method versus SVI.
Moreover, our approach appears highly competitive with Ackerer et al., (), which performs instance-by-instance volatility smoothing using classical neural networks and reports a MAPE of around 1% for synthetically generated data.(^8^88Compare Table 1 of Ackerer et al., (). We note that Ackerer et al., () does not perform a similar restriction of the domain of the volatility surface.)
We reenact its backtesting for the period Jan-Apr 2018 using our method and summarize the results in Table 5 in Appendix C.5.(^9^99We emphasize the following aspects. First, the backtest involves dropping half of input datapoints for each surface, and our trained operator continues to perform accurate smoothing, a strong indicator of the robustness of our approach with respect to subsampling of inputs afforded by neural operator’s discretization invariance. Second, while Ackerer et al., () requires to train 61 neural networks to perform the backtest once, our operator approach enables us to run 25 repetitions in around two minutes on a consumer grade laptop CPU, which is the average time that it takes Ackerer et al., () to train one network.)

An important fact to consider when analyzing Figure 3 is the historical tightening of the bid-ask spread (displayed in the third row), driven by increases in competition on the S&P 500 option market.
This explains why $\delta_{\text{spr}}$ is very small early in the training dataset, both for operator deep smoothing and SVI, while $\delta_{\text{abs}}$ is large:
Wide spreads make $\delta_{\text{spr}}$ more lenient an error metric but are accompanied by noisier prices, necessitating greater need for correction by the smoothing algorithm, in turn captured by $\delta_{\text{abs}}$.(^10^1010This is visually discernible from a comparison of Figure 8 and Figure 10.)
This argumentation extends to spikes in the spread, which indicate periods of market stress, and then allows to explain spikes in $\delta_{\text{abs}}$.

Complementary to Figure 3, Figure 4 resolves the error metrics as well as the terms controlling the absence of arbitrage spatially, averaged over time.
$\delta_{\text{abs}}$ tends to be larger on the Call side (positive log-moneyness), in accordance with Call option’s noisier prices (Call options experience less trading than Put options).
Moreover, we discern that, on average, the smoothed surfaces are completely free of arbitrage (indicated by non-negativity).

Figure: Figure 4: Average spatial distribution of benchmark metrics and arbitrage terms over $\mathcal{D}_{\text{test}}$.
Refer to caption: extracted/5912674/resources/finetune/spatial_metrics_test.png

#### Generalization

To test the generalization capabilities of our approach, we procure end-of-day options data for the S&P 500 (SPX), the NASDAQ-100 (NDX), the Dow Jones Industrial Average (DJX), and the Russell 2000 (RUT) for the month of January 2021 from the *OptionMetrics Ivy DB US* database, accessed by us through the *Wharton Research Data Services* (WRDS).
We evaluate the trained operator on the data for the month of January 2021 (right after the training period on the S&P 500 intraday data), and report the average error metrics $\delta_{\text{abs}}$ and $\delta_{\text{spr}}$ as well as the average arbitrage losses $\mathcal{L}_{\text{cal}}$ and $\mathcal{L}_{\text{but}}$ in Table 1.
Firstly, our method maintains its performance on end-of-day S&P 500 data, validating the soundness of our approach: While end-of-day data is slightly different from intraday data, our method still yields small error metrics and arbitrage-free prices.
Secondly, the method generalizes well to other indices.
We want to stress the fact that our operator has solely been trained on intraday S&P 500 data.
Its accurate and virtually arbitrage-free output on end-of-day data of other indices is a strong indicator of the robustness of our approach.
We provide further example plots for these datasets in Appendix D.3.

**Table 1: Average error metrics and arbitrage losses for end-of-day options data for US indices in January 2021. The GNO has been trained solely on intraday S&P 500 data from before 2021.**
|  | SPX | NDX | DJX | RUT |
| --- | --- | --- | --- | --- |
| $\delta_{\text{abs}}$ | 0.00272 | 0.01057 | 0.01629 | 0.00885 |
| $\delta_{\text{spr}}$ | 0.64423 | 1.53306 | 0.20736 | 1.03183 |
| $\mathcal{L}_{\text{but}}$ | 6.41e-05 | 2.16e-04 | 1.61e-03 | 4.48e-06 |
| $\mathcal{L}_{\text{cal}}$ | 0.00000 | 0.00000 | 3.13e-08 | 0.00000 |

## 5 Discussion

#### Summary

We provide a novel method for implied volatility smoothing, resulting from an application of our general operator deep smoothing approach for discretization-invariant data interpolation.
The approach leverages a graph neural operator to directly map given data – consistently across size and spatial arrangement – to smoothed surfaces, transcending classical parametric smoothing techniques.
In the example of volatility smoothing, benefits include a massively simplified online calibration process.

#### Learning from Large Datasets

By moving the application of neural networks from the instance-by-instance level (Ackerer et al.,, ) to the “operator level”, we leverage the information contained in the entire training dataset for the smoothing of every single surface.
In other words, our method “unlocks” large historical options datasets for volatility smoothing.
We argue that our substantial outperformance against Ackerer et al., () in the “Extrapolation-Test”-setting of the benchmark detailed in Table 5 of Appendix C.5 owes to this circumstance.

#### Subsampling of Inputs

The discretization-invariance of the GNO entails that our method is robust with respect to subsampling of inputs.
In practice, subsampling of inputs occurs in the context of outlier removal.
In the example of volatility smoothing, certain quotes may be determined spurious.
Simply removing anomalous datapoints from the input is compatible with our method (moreover, we leverage this fact during operator training to improve generalization, compare Appendix C.3).

#### Compression

Figure 3 makes the compression qualities of the operator deep smoothing approach apparent:
We compute the entire historical timeseries using a single GNO instance, with around 100 thousand parameters.
Evaluating the SVI benchmark, on the other hand, requires 61454 model instances (one per slice), or a total of 307270 parameters.
A comparison with Ackerer et al., (), which for each smoothed surface trains a new neural network of around 5085 parameters,(^11^1111Computed as the sum of $120=3\times 40$ parameters for the input layer, three times $1640=41\times 40$ parameters for the hidden layers, $41$ parameters for the output layers, plus $4$ additional parameters of the SSVI prior and a scaling parameter.) is striking:
Smoothing of the CBOE dataset 2012–2021 at its 20-minute interval frequency would require more than 200 million parameters (more with rising frequency).
At the same time, we expect our GNO to perform accurate smoothing over the entire training period and beyond (with regular finetuning), and our model instance remains fixed, even when moving to higher-frequency data.

#### Limitations and Perspectives

Compared to ad-hoc volatility parametrizations like SVI, the operator deep smoothing approach loses interpretability of parameters, which for some practitioners may be a stringent requirement.
This disadvantage is generally shared by neural network based engineering solutions.
Moreover, in some situations dimensionality reduction (even without interpretability of parameters) may be a desirable additional feature that is not directly achieved by our operator deep smoothing approach.
Combining the VAE method (Bergeron et al.,, ) with our operator deep smoothing approach could lead to further promising potential applications of neural operators.
Huang et al., () introduces *neural mappings*, which generalize neural operators to mixed infinite-/finite-dimensionality for input or output spaces.
This motivates a discretization-invariant GNO-based encoder, fit to handle raw incoming market data, and a classical decoder to extend the operator deep smoothing approach to a VAE-like architecture.

## Acknowledgements

Access to data was made possible through the Clearify Project funding by an Imperial College SME engagement grant jointly with Zeliade Systems, as well as the MSc in Mathematics and Finance, Imperial College London.
AJ gratefully acknowledges financial support from the EPSRC grants EP/W032643/1 and EP/T032146/1.
Computational resources and support were provided by the Imperial College Research Computing Service (Harvey,, ).
For the purpose of open access, the author(s) has applied a Creative Commons Attribution (CC BY) licence (where permitted by UKRI, ‘Open Government Licence’ or ‘Creative Commons Attribution No-derivatives (CC BY-ND) licence’ may be stated instead) to any Author Accepted Manuscript version arising.

## Appendix A Neural Operators

We give a review of the kernel integral transform neural operator framework of Kovachki et al., (), and expand in more detail on its graph neural operator.

#### Notation and Terms

Let $\mathcal{A}$ and $\mathcal{U}$ be input and output space of an operator learning problem, as introduced in Section 3.1.
Then, $\mathcal{A}$ and $\mathcal{U}$ are Banach spaces of functions $D\to\mathbb{R}^{c_{\text{in}}}$ and $D\to\mathbb{R}^{c_{\text{out}}}$, respectively, where $D$ is a (bounded) domain in $\mathbb{R}^{d}$.
The mathematical analysis of neural operators in Kovachki et al., () summarized hereafter, as well as the definition of universality and discretization-invariance in Section 3.1, make use of the following terms and notations.

In the context of Kovachki et al., (), a domain is a bounded and connected open set that is topologically regular (in the sense that it is the interior of its closure).
A domain $D$ is *Lipschitz* if its boundary locally is the graph of a Lipschitz continuous function defined on an open ball of $\mathbb{R}^{d-1}$.
An *open ball* – for any metric space $\mathcal{X}=(\mathcal{X},d)$ – is the set

$$ $\mathbb{B}_{\mathcal{X}}(x,\varepsilon)=\{y\in\mathcal{X}\colon d(y,x)< \varepsilon\}.$ $$

A *discrete refinement* of $\mathcal{X}$ is a nested sequence $(\pi_{n})$ of discretizations of $\mathcal{X}$ (finite subsets of $\mathcal{X}$), such that for every $\varepsilon>0$ there is $N\in\mathbb{N}$ such that $\{\mathbb{B}(x,\varepsilon)\colon x\in\pi_{N}\}$ covers $\mathcal{X}$.

We consider the space $C(\mathcal{A},\mathcal{U})$ of continuous operators between $\mathcal{A}$ and $\mathcal{U}$.
$C(\mathcal{A},\mathcal{U})$ is topologized by uniform convergence on compact sets.
With respect to this topology, a sequence $(F_{n})_{n\in\mathbb{N}}$ in $C(\mathcal{A},\mathcal{U})$ converges with limit $F\in C(\mathcal{A},\mathcal{U})$, if, for every $\varepsilon>0$ and every compact set $K$ in $\mathcal{A}$, it holds

$$ $\lim_{n\to\infty}\|F_{n}-F\|_{\infty,K}=0.$ $$

Here,

$$ $\|H\|_{\infty,K}=\sup_{a\in K}\|H(a)\|_{\mathcal{U}},\quad H\in C(\mathcal{A}, \mathcal{U}).$ $$

It is well known that this topology on $C(\mathcal{A},\mathcal{U})$ is induced by the metric

$$ $\rho(F,G)=\sum_{n=0}^{\infty}\frac{\|G-F\|_{\infty,\overline{\mathbb{B}}_{C( \mathcal{A},\mathcal{U})}(0,n)}}{1\vee\|G-F\|_{\infty,\overline{\mathbb{B}}_{C (\mathcal{A},\mathcal{U})}(0,n)}},\quad F,G,\in C(\mathcal{A},\mathcal{U}).$ $$

Therefore, the notion of *density* in $C(\mathcal{A},\mathcal{U})$, as used to define universality of neural operators in Section 3.1, is well defined.

### A.1 Kernel Integral Neural Operators

#### Kernel Integral Transform Neural Operators and Universality

A kernel integral transform neural operator consists of the sequential application of:

- 1)
A *lifting layer*
$L_{\mathcal{P}}\colon[a\colon D\to\mathbb{R}^{c_{\text{in}}}]\mapsto[h_{0}
\colon D\to\mathbb{R}^{c_{0}},\;h_{0}(x)=\mathcal{P}(a(x))],$
given by the pointwise application of a function $\mathcal{P}\colon\mathbb{R}^{c_{\text{in}}}\to\mathbb{R}^{c_{0}}$.
- 2)
The forward propagation through $J$ neural operator layers $L_{0},\dots,L_{J-1}$:
$[h_{0}\colon D\to\mathbb{R}^{c_{0}}]\xmapsto{L_{0}}[h_{1}\colon D\to\mathbb{R}
^{c_{1}}]\xmapsto{L_{1}}\dots\xmapsto{L_{J-1}}[h_{J-1}\colon D\to\mathbb{R}^{c
_{J}}];$
each layer $L_{j}$ operates as
$h_{j+1}(y)=\sigma_{j}\left(W_{j}h_{j}(y)+\int_{D}\kappa_{j}(y,x)h_{j}(x)dx+b_{
j}(y)\right),\quad y\in D,$
(11)
where
•
$W_{j}\in\mathbb{R}^{c_{j+1}\times c_{j}}$ is a weight matrix applied pointwise,
•
$\kappa_{j}\in C(D\times D,\mathbb{R}^{c_{j+1}\times c_{j}})$ is a kernel function parametrizing the integral transform and subject to integrability conditions,
•
The bias term $b_{j}$ is itself a function from $D$ to $\mathbb{R}^{c_{j+1}}$,
•
$\sigma_{j}$ is a classical neural network activation function.
- 3)
A *projection layer*
$L_{\mathcal{Q}}\colon[h_{J}\colon D\to\mathbb{R}^{c_{J}}]\mapsto[u\colon D\to
\mathbb{R}^{c_{\text{out}}},\;u(x)=\mathcal{Q}(h_{J}(x))],$
given by the pointwise application of a function $\mathcal{Q}\colon\mathbb{R}^{c_{J}}\to\mathbb{R}^{c_{\text{out}}}$.

In practice, all components (lifting, kernel functions, projection) are implemented as classical feedforward neural networks (FFNs).
This neural operator architecture is universal in the following sense.

###### Theorem A.1 (Universal Approximation; Theorem 11 of Kovachki et al., ( 2023 ) ) .

Let $D$ be a (bounded) Lipschitz domain.
Assume:

- •
$\mathcal{A}=W^{k_{1},p_{1}}(D)$ for $k\in\mathbb{N}_{\geq 0}$ and $1\leq p_{1}<\infty$, or $\mathcal{A}=C(\overline{D})$.
- •
$\mathcal{U}=W^{k_{2},p_{2}}(D)$ for $k\in\mathbb{N}_{\geq 0}$ and $1\leq p_{2}<\infty$, or $\mathcal{U}=C(\overline{D})$.

Then, a subset of kernel integral neural operators, with kernel functions and bias functions taken from a suitable set of FNNs, is dense in $C(\mathcal{A},\mathcal{U})$.

###### Theorem A.1 (Universal Approximation; Theorem 11 of Kovachki et al., ( 2023 ) ) .

#### Discretization-Invariant Implementations

Consider a neural operator $F^{\theta}$ and let $\pi=\{x_{l}\}_{l=1}^{p}$ be a discretization of $D$.
To make sense of a basic discretization-invariant implementation for $F^{\theta}$, associate with $\pi$ a partition $(D_{1},\dots,D_{p})$ of $D$ for which $\lambda_{d}(D_{l})>0$ and $x_{l}\in D_{l}$ for $l=1,\dots,p$. Here $\lambda_{d}$ denotes the Lebesgue measure on $\mathbb{R}^{d}$.
Consider the following implementation of $F^{\theta}$ (written in terms of a single constituent layer $L=(W,\kappa,b,\sigma)$):

$$ $\tilde{L}(h|_{\pi})(y)=\sigma\left(Wh(y)+\sum_{l=1}^{p}\kappa(y,x_{l})h(x_{l}) \lambda_{d}(D_{l})+b(y)\right),\quad y\in D.$ (12) $$

Kovachki et al., () establishes the following.

###### Theorem A.2 (Discretization Invariance; Theorem 8 of Kovachki et al., ( 2023 ) ) .

Let $F^{\theta}\colon\mathcal{A}\to\mathcal{U}$ be a kernel integral neural operator, where $\mathcal{A}$ and $\mathcal{U}$ both continuously embed into $C(\overline{D})$.
Then, the implementation of $F^{\theta}$ based on (12) is discretization-invariant as defined in Section 3.1.

(12) suggests the straightforward (quasi) Monte-Carlo inspired implementation

$$ $\tilde{L}(h|_{\pi})(y)=\sigma\left(Wh(y)+\frac{\lambda_{d}(D)}{|\pi|}\sum_{x \in\pi}\kappa(y,x)h(x)+b(y)\right),\quad y\in D.$ (13) $$

Most effectively, $\pi$ is a low-discrepancy sequence in $D$.

###### Theorem A.2 (Discretization Invariance; Theorem 8 of Kovachki et al., ( 2023 ) ) .

### A.2 Graph Neural Operators

The curse of dimensionality makes the direct implementation (13) prohibitively expensive in practice.
Instead, Anandkumar et al., () introduces graph neural operators (or, GNOs, for short) which replace the kernel integral operation at the heart of the framework by a sum approximation and organizes the constituent terms using a directed graph structure:
The discretization $\pi=\{x_{1},\dots,x_{m}\}$ of the input data $\mathbf{h}=h|_{\pi}$ is enriched with a directed graph structure $G_{\mathbf{h}}=(V,E)$, allowing the following implementation $\tilde{F}^{\theta}$ of $F^{\theta}$:

$$ $\tilde{L}(h|_{\pi})(y)=\sigma\left(Wh(y)+\frac{1}{|\mathcal{N}_{\text{in}}(y)| }\sum_{x\in\mathcal{N}_{\text{in}}(y)}\kappa(y,x)h(x)+b(y)\right),\quad y\in V.$ (14) $$

Here, $\mathcal{N}_{\text{in}}(y)$ is the set of so-called in-neighbors of $y$ in the graph $G_{\mathbf{h}}$: $x\in\mathcal{N}_{\text{in}}(y)$ iff $(x,y)\in E$.
It is clear that, to compute output at $y$, the point $y$ must be included as a node into the graph $G_{\mathbf{h}}$.
On the other hand, it is necessary to drop the local linear transformation with $W$ if $y$ is not part of the input data locations $\pi$ (compare our discussion in 3).
It is important to reconcile the input and output locations of the layers when creating the graph structure to enable an efficient implementation using *message passing* algorithms.

The graph structure can be meticulously adjusted to implement various complexity-reducing techniques like Nyström approximation or integration domain truncation that effectively aim at a systematic reduction of the size of $\mathcal{N}_{\text{in}}(y)$; the naive implementation (12) is recovered for the case of a complete directed graph (with self-loops) for which $\mathcal{N}_{\text{in}}(y)=\pi$.
Note that choosing $\mathcal{N}_{\text{in}}(y)$ as a strict subset of $\pi$ breaks the guaranteed smoothness in the GNO output.

## Appendix B Interpolation Graph Neural Operator

We detail our modifications of the GNO architecture.
A generic Pytorch implementation is available at [https://github.com/rwicl/operator-deep-smoothing-for-implied-volatility](https://github.com/rwicl/operator-deep-smoothing-for-implied-volatility).

Let $\mathbf{v}=\{v(x)\}_{x\in\pi}$ be the given data.

#### Graph Construction

Arguably part of the model architecture is the graph construction:
Compile the set $\pi_{\text{out}}=\{y_{l}\}_{l=0}^{q}$ of points $y\in D$ at which to compute the smoothed surface $\hat{v}(y)$.
During operator training, this will be the set of input locations (to compute the fitting loss) as well as any additional locations needed to compute auxiliary loss terms (the arbitrage losses $\mathcal{L}_{\text{but}}(\theta,\mathbf{v}))$ and $\mathcal{L}_{\text{cal}}(\theta,\mathbf{v})$ in the case of volatility smoothing).
For each $y\in\pi_{\text{out}}$, we construct the set of in-neighbors from the set of input data locations:

$$ $\mathcal{N}_{\text{in}}(y)\subseteq\pi_{\text{in}}.$ (15) $$

In other words, we employ a Nyström approximation with nodes limited to the input data locations, which allows us to employ kernel functions with input skip connections.
We set $G_{\mathbf{v}}=(\pi_{\text{out}},E)$, where

$$ $E=\bigcup_{y\in\pi_{\text{out}}}\{(x,y)\colon x\in\mathcal{N}_{\text{in}}(y)\}.$ $$

#### Forward Propagation

Given $G_{\mathbf{v}}$, we perform the first step of the forward propagation as follows:

$$ $\left\{\begin{aligned} \tilde{h}_{0}(x)&=\mathcal{P}_{0}(v(x)),&x\in\pi_{\text {in}}\\ h_{1}(y)&=(\sigma_{0}\circ\mathcal{Q}_{0})\left(\mathcal{K}(\tilde{\mathbf{h}} _{0};\mathbf{v})(y)+b_{0}\right),&y\in\pi_{\text{out}}.\end{aligned}\right.$ $$

For the subsequent layers $j=1,\dots,J-1$, we then proceed using the classical scheme:

$$ $\left\{\begin{aligned} \tilde{h}_{j}(y)&=\mathcal{P}_{j}(h_{j}(y)),\\ h_{j+1}(y)&=(\sigma_{j}\circ\mathcal{Q}_{j})\left(W_{j}\tilde{h}_{j}(y)+ \mathcal{K}_{j}(\tilde{\mathbf{h}}_{j};\mathbf{v})(y)+b_{j}\right)\end{aligned }\right.,\quad y\in\pi_{\text{out}}.$ $$

In the above:

- •
$\mathcal{P}_{j}\colon\mathbb{R}^{c_{j}}\to\mathbb{R}^{\tilde{c}_{j}}$ and $\mathcal{Q}_{j}\colon\mathbb{R}^{\tilde{c}_{j+1}}\to\mathbb{R}^{c_{j+1}}$ are layer-individual lifting and projection, implemented as FNNs.
- •
$W_{j}\in\mathbb{R}^{\tilde{c}_{j+1}\times\tilde{c}_{j}}$ is a weight matrix (not present for $j=0$), while $b_{j}\in\mathbb{R}^{\tilde{c}_{j+1}}$ is a constant bias term.
- •
$\mathcal{K}_{j}$ is the sum approximation of the kernel integral with kernel weight function $\kappa^{W}_{j}\colon D^{2}\times\mathbb{R}^{\tilde{c}_{j}}\times\mathbb{R}^{c_
{0}}\to\mathbb{R}^{\tilde{c}_{j+1}\times\tilde{c}_{j}}$ and kernel bias function $\kappa^{b}_{j}\colon D^{2}\times\mathbb{R}^{\tilde{c}_{j}}\times\mathbb{R}^{c_
{0}}\to\mathbb{R}^{\tilde{c}_{j+1}}$ (both with state and input skip connections):
$\mathcal{K}_{j}(\tilde{\mathbf{h}}_{j};\mathbf{v})(y)=\frac{1}{|\mathcal{N}_{
\text{in}}(y)|}\sum_{x\in\mathcal{N}_{\text{in}}(y)}\kappa^{W}_{j}(y,x,\tilde{
h}_{j}(x);v(x))\tilde{h}_{j}(x)+\kappa^{b}_{j}(y,x,\tilde{h}_{j}(x);v(x)).$
Both $\kappa^{W}_{j}$ and $\kappa^{b}_{j}$ are implemented as FNNs.

Note that omitting the local linear transformation in the first layer allows to extract the fist hidden state $h_{1}(y)$ for all $y\in\pi_{\text{out}}$ from the lifted input $\tilde{\mathbf{h}}_{0}$, which is defined solely for the input locations $x\in\pi_{\text{in}}$.
Providing each layer with its own lifting and projection allows to separate the hidden channel size $c_{0},\dots,c_{J}$ from the the dimensions $\tilde{c}_{0},\dots,\tilde{c}_{J}$ of the space in which the integral transformation is performed.
Moreover, the individual lifting and projection help re-parametrize the state before performing the integral transform (inspired by the succesful Transformer architecture), which allows to keep the size of the kernel weight matrix low.

## Appendix C Supplementary Information: Data, Model, Training, Evaluation

This section contains additional information regarding our empirical study of the operator deep smoothing method for implied volatility smoothing.

### C.1 Data

#### Data Source

Our numerical experiments are based on the ”Option Quotes” dataset product available for purchase from the CBOE.
Our version of the dataset contains relevant data features for S&P 500 Index options for the years 2012 through 2021 and is summarized on a 20-minute interval basis.

#### Data Preparation

We compute *simple mids* for options and underlying by averaging bid and ask quotes and use this aggregate as a reference price for all our subsequent computations.
We calculate discount factors and forward prices from Put-Call parity using the industry standard technique based on linear regression.
We compute time-to-expiry in units of one year as well as log-moneyness as defined in Section 2.
We extract implied volatilities using the *py-vollib-vectorized* project available in the *Python Package Index* at the location [https://pypi.org/project/py-vollib-vectorized/](https://pypi.org/project/py-vollib-vectorized/).
*py-vollib-vectorized* implements a vectorized version of Jäckel, ()’s *Let’s-be-rational* state-of-the-art method for computing implied volatility.
We discard all implied volatilities of in-the-money options, or, in other words, we compose our implied volatility surface from Put options for non-positive log-moneyness values and Call options for positive log-moneyness values.

### C.2 Model

We configure the modified GNO architecture introduced in Appendix B as follows:

- •
We impose the following restriction on the in-neighborhood sets $\mathcal{N}_{\text{in}}(y)$:
$\mathcal{N}_{\text{in}}(y)\subseteq\{(\rho_{l},z_{l})\in\pi\colon|\rho_{l}-
\rho|\leq\overline{\rho}\}.$
(16)
This choice stems from the consideration that implied volatility smoothing requires limited global informational exchange along the time-to-expiry axis.
We then generate $\mathcal{N}_{\text{in}}(y)$ of certain size $K$ using a high-discrepancy subsampling heuristic detailed in the supplemented code material.
Specifically, we use the values $\overline{\rho}=0.3$ and $K=50$.
- •
We employ three hidden layers and a channel size of 16: $J=4$, and $c_{1},c_{2},c_{3}=16$ ($c_{0}$ and $c_{J}$ are determined as $1$ by the scalar dimension of volatility data). We use GELU-activations for the hidden layers and a Softplus-activation for the output layer (to ensure the positivity of the smoothed surfaces): $\sigma_{0},\dots,\sigma_{J-1}=\operatorname{GELU}$, and $\sigma_{J}=\operatorname{Softplus}$.
- •
We retain $\mathcal{P}_{0},\dots,\mathcal{P}_{J-1}$ and $\mathcal{Q}_{J}$ as single-hidden layer FNNs with 64 hidden nodes and GELU-activations for the hidden layers. The remaining lifting and projections remain unutilized. In particular, $\tilde{c}_{0},\dots,\tilde{c}_{J}=16$.
- •
We implement the kernel weight and bias functions as two-hidden layer FNNs with 64 hidden nodes and GELU-activations for the hidden layers.

This configuration has been converged upon using manual experimentation and amounts to a total number of 102529 trainable parameters.

### C.3 Training

#### Loss Function

To ease notation we write $\hat{v}^{\theta}=\tilde{F}^{\theta}(\mathbf{v})$.
We implement a Vega-weighted version of the fitting loss (8):

$$ $\mathcal{L}_{\text{fit}}(\theta;\mathbf{v})=\left(\frac{1}{|\pi_{\mathbf{v}}|} \sum_{x\in\pi_{\mathbf{v}}}w_{\mathcal{V}}(x;\mathbf{v})\left|(\hat{v}^{\theta }(x)-v(x))/v(x)\right|^{2}\right)^{1/2}.$ $$

Here,

$$ $w_{\mathcal{V}}(x;\mathbf{v})=\frac{\mathcal{V}(x,v(x))}{\frac{1}{|\pi|}\sum_{ x\in\pi}\mathcal{V}(x,v(x))}\vee 1,$ $$

where $\mathcal{V}(x,v(x))$ is the *Black-Scholes Vega*, the sensitivity of the Black-Scholes option price with respect to its volatility parameter:

$$ $\mathcal{V}(\tau,x,v)=\partial_{v}\mathrm{BS}(\tau,k,v)=\varphi(d_{1}(\tau,k,v ))\sqrt{\tau}.$ (17) $$

defined in Section 1 (adjusted for transformed coordinates).

For the implementation of the no-arbitrage penalization terms $\mathcal{L}_{\text{but}}$ and $\mathcal{L}_{\text{cal}}$, we first generate discretizations $\pi_{\rho}=\{\rho_{1},\dots,\rho_{m}\}$ and $\pi_{z}=\{z_{1},\dots,z_{n}\}$ of $[\rho_{\text{min}},\rho_{\text{max}}]$ and $[z_{\text{min}},z_{\text{max}}]$.
We resolve the derivative terms $\partial_{z}v^{\theta}$ and $\partial_{z}^{2}v^{\theta}$ on the synthetic rectilinear grid $\pi=\pi_{\rho}\times\pi_{z}$ using (central) finite differences.
Then, we translate $\mathcal{L}_{\text{but}}$ directly from (9) as

$$ $\mathcal{L}_{\text{but}}(\theta;\mathbf{v},\pi)=\frac{1}{|\pi|}\sum_{x\in\pi} \left(\operatorname{But}(x;\tilde{v}^{\theta}(x),\Delta_{z,\pi}\tilde{v}^{ \theta}(x),\Delta_{z,\pi}^{2}\tilde{v}^{\theta}(x))-\varepsilon\right)^{-},$ $$

where $\operatorname{But}$ is made consistent with the transformed coordinates and we used obvious notation for the finite differences.
We use $\varepsilon=10^{-3}$.
On the other hand, we enforce the monotonicity constraint of Theorem 2.1 using

$$ $\mathcal{L}_{\text{cal}}(\theta;\mathbf{v},\pi_{\rho},\pi_{z})=\frac{1}{mn} \sum_{i=1}^{m}\sum_{j=1}^{n}\left(\frac{\tilde{v}^{\theta}(\rho_{i+1},z_{j})}{ \tilde{v}^{\theta}(\rho_{i},(\rho_{i+1}z_{j})/\rho_{i})}-\frac{\rho_{i}}{\rho_ {i+1}}-\varepsilon\right)^{-}.$ $$

Compared to a derivative based implementation, this implementation is independent of the scale and – in our empirical experiments – has provided an improved signal.
Since the Nyström approximation employed by the graph neural operator (as well as the choice (16)) break the guaranteed smoothness of the operator output, we additionally introduce $\|\partial_{\rho}^{2}\hat{v}^{\theta}\|_{2}$ and $\|\partial^{2}_{z}\hat{v}^{\theta}\|_{2}$ as regularization terms:

$$ $\mathcal{L}_{\text{reg-$\rho$}}(\theta;\mathbf{v},\pi)=\sqrt{\frac{1}{|\pi|} \sum_{x\in\pi}|\Delta_{\rho,\pi}^{2}\tilde{v}^{\theta}(x)|^{2}},\quad\mathcal{ L}_{\text{reg-$z$}}(\theta;\mathbf{v},\pi)=\sqrt{\frac{1}{|\pi|}\sum_{x\in\pi} |\Delta_{z,\pi}^{2}\tilde{v}^{\theta}(x)|^{2}}.$ $$

We compose the final loss function as a weighted sum of all terms introduced:

$$ $\mathcal{L}(\theta;\mathbf{v},\pi_{\rho},\pi_{z})=\sum\left\{\begin{aligned} & \lambda_{\text{fit}}\mathcal{L}_{\text{fit}}(\theta;\mathbf{v}),\\ &\lambda_{\text{but}}\mathcal{L}_{\text{but}}(\theta;\mathbf{v},\pi_{\rho} \times\pi_{z}),\\ &\lambda_{\text{cal}}\mathcal{L}_{\text{cal}}(\theta;\mathbf{v},\pi_{\rho},\pi _{z}),\\ &\lambda_{\text{reg-$\rho$}}\mathcal{L}_{\text{reg-$\rho$}}(\theta;\mathbf{v}, \pi_{\rho}\times\pi_{z}),\\ &\lambda_{\text{reg-$z$}}\mathcal{L}_{\text{reg-$z$}}(\theta;\mathbf{v},\pi_{ \rho}\times\pi_{z}).\end{aligned}\right.$ $$

The specific weights are

| $\lambda_{\text{fit}}$ | $\lambda_{\text{cal}}$ | $\lambda_{\text{but}}$ | $\lambda_{\text{reg-$\rho$}}$ | $\lambda_{\text{reg-$z$}}$ |
| --- | --- | --- | --- | --- |
| $1$ | $10$ | $10$ | $0.01$ | $0.01$ |

The particular weighting of the individual terms has been retrieved by manual experimentation, led by the findings of Ackerer et al., ().

#### Validation Loss

Table 2 displays descriptive statistics of the validation losses.

**Table 2: Validation Loss**
|  | mean | std | min | 25% | 50% | 75% | max |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $\mathcal{L}$ | 0.062 | 0.032 | 0.029 | 0.053 | 0.059 | 0.064 | 0.506 |
| $\mathcal{L}_{\text{fit}}$ | 0.022 | 0.016 | 0.009 | 0.016 | 0.020 | 0.024 | 0.232 |
| $\mathcal{L}_{\text{but}}$ | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
| $\mathcal{L}_{\text{cal}}$ | 0.000 | 0.001 | 0.000 | 0.000 | 0.000 | 0.000 | 0.008 |
| $\mathcal{L}_{\text{reg-$\rho$}}$ | 0.644 | 0.137 | 0.196 | 0.538 | 0.641 | 0.747 | 0.996 |
| $\mathcal{L}_{\text{reg-$z$}}$ | 0.143 | 0.254 | 0.048 | 0.083 | 0.101 | 0.134 | 5.154 |

### C.4 Evaluation

Here we provide additional results to supplement our performance evaluation of operator deep smoothing (OpDS) for implied volatility.
Table 3 and Table 4 display descriptive statistics of OpDS versus SVI.
OpDS* refers to benchmarking without monthly finetuning.
Moreover, we include Figures 5 and 7, which are versions of Figures 3 and  4 with the only difference being that Figures 5 and 7 were created using OpDS* instead of OpDS. Finally, Figure 6 shows the average spatial distribution of benchmark metrics and arbitrage term over the training dataset, which complements the same averages on the test dataset shown in Figures 3.

**Table 3: Descriptive statistics $\delta_{\text{abs}}$ over $\mathcal{D}_{\text{val}}$/$\mathcal{D}_{\text{test}}$.**
|  | mean | std | 1% | 25% | 50% | 75% | 99% |
| --- | --- | --- | --- | --- | --- | --- | --- |
| OpDS | 0.009/0.005 | 0.007/0.001 | 0.003/0.003 | 0.006/0.004 | 0.008/0.005 | 0.010/0.005 | 0.021/0.007 |
| OpDS* | 0.009/0.007 | 0.007/0.001 | 0.003/0.003 | 0.006/0.07 | 0.008/0.007 | 0.010/0.008 | 0.021/0.012 |
| SVI | 0.021/0.015 | 0.006/0.002 | 0.007/0.010 | 0.016/0.013 | 0.020/0.014 | 0.025/0.016 | 0.034/0.020 |

**Table 4: Descriptive statistics $\delta_{\text{spr}}$ over $\mathcal{D}_{\text{val}}$/$\mathcal{D}_{\text{test}}$.**
|  | mean | std | 1% | 25% | 50% | 75% | 99% |
| --- | --- | --- | --- | --- | --- | --- | --- |
| OpDS | 0.479/1.265 | 0.662/0.347 | 0.193/0.609 | 0.274/1.025 | 0.330/1.240 | 0.550/1.451 | 1.526/2.453 |
| OpDS* | 0.479/1.866 | 0.662/0.574 | 0.193/0.731 | 0.274/1.457 | 0.330/1.826 | 0.550/2.233 | 1.526/3.445 |
| SVI | 1.124/3.382 | 0.877/0.826 | 0.301/1.464 | 0.492/2.827 | 0.715/3.320 | 1.646/3.914 | 3.710/5.247 |

Figure: Figure 5: OpDS*: Benchmark metrics (surface-averages, between 25%- and 75%-quantiles) over training period (computed from validation dataset, resampled monthly) and testing period (resampled every two days).
Refer to caption: extracted/5912674/resources/raw/benchmark_quantiles.png

Figure: Figure 6: OpDS: Average spatial distribution of benchmark metrics and arbitrage term over train dataset.
Refer to caption: extracted/5912674/resources/raw/spatial_metrics_train.png

Figure: Figure 7: OpDS*: Average spatial distribution of benchmark metrics and arbitrage terms over test dataset (non-finetuned model).
Refer to caption: extracted/5912674/resources/raw/spatial_metrics_test.png

### C.5 Comparison to Classical Neural Networks

We reproduce the experiment underlying Table 1 of Ackerer et al., () using our operator deep smoothing approach.
Given observed option quotes, it involves dropping 50% of datapoints, and then measuring the MAPE of the smoothed surface at retained datapoints (“Train”) as well as dropped datapoints (”Test”).
“Interpolate” and “Extrapolate” are different settings dictating how exactly the datapoints which to drop are selected (for details refer to Ackerer et al., ()).
The experiment is performed on end-of-day S&P 500 data in the period from January to April 2018 and averaged percentiles are reported.
Table 5 reproduces the relevant row from Table 1 of Ackerer et al., () (“DS”) as well as our results averaged over 25 repetitions (“OpDS”).

**Table 5: Backtesting results of Operator Deep Smoothing vs. Deep Smoothing (“DS”; taken from [2]); quantiles in %, Jan-Apr 2018 end-of-day SPX data.**
|  | Interpolation | Extrapolation |  |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| MAPE | Train | Test | Train | Test |  |  |  |  |  |  |  |  |  |
|  | $\lambda$ | $q_{05}$ | $q_{50}$ | $q_{95}$ | $q_{05}$ | $q_{50}$ | $q_{95}$ | $q_{05}$ | $q_{50}$ | $q_{95}$ | $q_{05}$ | $q_{50}$ | $q_{95}$ |
| OpDS | 10 | 0.5 | 0.7 | 1.0 | 0.5 | 0.7 | 1.1 | 0.5 | 0.7 | 1.0 | 0.7 | 0.9 | 1.3 |
| DS | 10 | 0.5 | 0.7 | 1.2 | 0.5 | 0.8 | 1.2 | 0.4 | 0.6 | 0.9 | 1.2 | 1.7 | 2.4 |

## Appendix D Additional Plots

### D.1 Example Plots

We plot the results of operator deep smoothing (OpDS) vs. SVI on example inputs.
To aid the visual clarity of our plots, we display only every third market datapoint.

Figure: Figure 8: Smoothing of quotes $\mathbf{v}\in\mathcal{D}_{\text{val}}$ from 20.07.2012 at 10:50:00.
Refer to caption: extracted/5912674/resources/svi_vs_opds_2012.png

Figure: Figure 9: Smoothing of quotes $\mathbf{v}\in\mathcal{D}_{\text{val}}$ from 21.10.2016 at 13:10:00.
Refer to caption: extracted/5912674/resources/svi_vs_opds_2016.png

Figure: Figure 10: Smoothing of quotes $\mathbf{v}\in\mathcal{D}_{\text{train}}$ from 04.01.2021 at 10:50:00.
Refer to caption: extracted/5912674/resources/svi_vs_opds_2021.png

### D.2 Monthly Breakdown of Spatial Distributions of Benchmark Metrics on Test Dataset

Figure: Figure 11: Average spatial distribution of $\delta_{\text{abs}}$ on $\mathcal{D}_{\text{test}}$ for OpDS with monthly finetuning, per month. Blank cells indicate that no data was available for the particular region in the respective month.
Refer to caption: extracted/5912674/resources/finetune/monthly_mape.png

Figure: Figure 12: Average spatial distribution of $\delta_{\text{fit}}$ on $\mathcal{D}_{\text{test}}$ for OpDS with monthly finetuning, per month. Blank cells indicate that no data was available for the particular region in the respective month.
Refer to caption: extracted/5912674/resources/finetune/monthly_spread.png

### D.3 Example Plots: Option Metrics End-Of-Day US Index Options Data

Figure: Figure 13: Smoothing of SPX end-of-day data from 07.01.2021. Every third datapoint displayed.
Refer to caption: extracted/5912674/resources/wrds/spx_2021-01-07.png

Figure: Figure 14: Smoothing of NDX end-of-day data from 07.01.2021. Every second datapoint displayed.
Refer to caption: extracted/5912674/resources/wrds/ndx_2021-01-07.png

Figure: Figure 15: Smoothing of DJX end-of-day data from 07.01.2021. Every datapoint displayed.
Refer to caption: extracted/5912674/resources/wrds/djx_2021-01-07.png

Figure: Figure 16: Smoothing of RUT end-of-day data from 07.01.2021. Every datapoint displayed.
Refer to caption: extracted/5912674/resources/wrds/rut_2021-01-07.png

## Appendix E SVI

SVI, originally devised by Gatheral (Gatheral,, ), stands for *Stochastic Volatility Inspired* and is a low-dimensional parametrization for implied volatility slices (namely for each maturity).
It captures the key features of implied volatility of Equity indices,
and has become the industry-wide benchmark for implied volatility smoothing on such markets.
Its ”raw” variant parametrizes the ”slice” of implied volatility at a given time-to-expiry $\tau$ as follows:

$$ $\hat{v}_{\tau}(k)=\sqrt{\frac{a+b\left(\rho(k-m)+\sqrt{(k-m)^{2}+\sigma^{2}} \right)}{\tau}},\quad\text{for all }k\in\mathbb{R},$ $$

where $a,b,k,m,\sigma$ are real parameter values.

#### Calibration

While stylistically accurate, SVI does not easily guarantee absence of static arbitrage opportunities,
and several authors have investigated this issue (Gatheral and Jacquier,, ; Martini and Mingone,, ; Mingone,, ; Martini and Mingone,, ).
To produce our SVI benchmark we therefore rely on the constrained SLSQP optimizer provided by the *SciPy* scientific computing package for Python,
with the mean square error objective, a positivity constraint and the constraint (2) (computed in closed form), and the following parameter bounds:

$$ $a\in\mathbb{R},\quad b\in[0,1],\quad\rho\in[-1,1],\quad m\in[-1.5,0.5],\quad \sigma\in[10^{-8},2].$ $$

We ignore the calendar arbitrage condition.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 6
presentation: 3
questions: See Weaknesses.
soundness: 3
strengths: - **Novelty:** This paper proposes a novel application of GNO architectures to implied volatility smoothing, a task traditionally reliant on parametric models like SVI. By leveraging the discretization-invariance of neural operators, the method effectively addresses the challenges of dynamic and irregular financial data, marking a significant step forward in financial engineering.

- **Practical Significance:** The operator deep smoothing approach eliminates the need for instance-by-instance recalibration, streamlining the online calibration process. This drastically reduces computational overhead, making the method highly practical for real-time applications in trading, risk management, and other financial operations.

- **Robust Validation:** The method is extensively validated on a decade of intraday S&P 500 options data, demonstrating superior accuracy compared to the SVI benchmark and other machine learning techniques. Its strong generalization performance on unseen datasets, including end-of-day options for other indices, further highlights its robustness and adaptability.

 - **Reproducibility:** The authors provide open-source code, model weights, and detailed implementation details, including architecture, loss functions, and training setups. These resources ensure the experiments are easy to replicate and extend.
summary: This paper introduces a novel approach to smoothing implied volatility using neural operators, focusing on Graph Neural Operator (GNO) architectures. Unlike traditional parametric models or classical neural networks, the proposed operator deep smoothing method directly maps observed implied volatilities to smooth, arbitrage-free surfaces, eliminating the need for instance-by-instance recalibration. By harnessing the unique capability of neural operators to handle data of varying sizes and spatial arrangements, the model effectively addresses the challenges posed by the dynamic and irregular nature of real-world financial data. Validated on a decade of intraday S&P 500 options data, the method outperforms the SVI benchmark and other machine learning techniques while generalizing effectively to unseen datasets, including end-of-day options for other indices, without retraining.
weaknesses: - **Limited Benchmark Comparisons:** The paper benchmarks its approach against SVI [1] and Ackerer et al. [2] but does not include comparisons with other key methods, such as SSVI [3] and VAE-based approaches [4]. Incorporating these would provide a more comprehensive evaluation. Additionally, using synthetic data, as in [2], could further strengthen the experimental validation.

 - **Insufficient Analysis of Computational Efficiency:** While the paper highlights the elimination of instance-by-instance recalibration, it lacks a detailed discussion on computational complexity, including runtime and memory requirements. These metrics are crucial for assessing the method’s practicality in high-frequency financial contexts.

 - **Abstract Treatment of Discretization-Invariance:** The explanation of discretization-invariance is largely theoretical. Providing concrete financial scenarios, such as handling abrupt market shifts, would better illustrate its practical significance and strengthen the narrative.

- **Clarity of Experimental Figures**: Figures 3 and 4 have unclear legends and insufficient annotations, which reduce their interpretability. Improving their clarity and labeling would make the results more accessible and impactful.

[1] A parsimonious arbitrage-free implied volatility parameterization with application to the valuation of volatility derivatives.

[2] Deep Smoothing of the Implied Volatility Surface.

[3] Arbitrage-free SVI volatility surfaces.

[4] Variational Autoencoders: A Hands-Off Approach to Volatility.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
questions: please refer to the weaknesses part
soundness: 3
strengths: 1. This paper is well-written and easy to follow.

2. This proposed method is resilient to input subsampling, aligns with no-arbitrage conditions, and eliminates the need for extensive data pre-processing.
summary: This paper proposes a novel approach to nowcast implied volatility using neural operators, advancing the technique of implied volatility smoothing by constructing consistent, smooth surfaces that reflect current option market prices. Unlike traditional machine learning models, which struggle with dynamically changing spatial configurations in option price data, this work leverages a graph neural operator architecture to achieve robust and accurate predictions. Extensive benchmarks showcase that this work achieves superior generalization and accuracy over conventional neural networks and the SVI parametrization.
weaknesses: 1. The contribution is limited. It is focused on applying the Graph Neural Operator (GNO) architecture to the specific task of nowcasting implied volatility, without modifications to the operator itself.

2. It's important to discuss various neural operators and clarify why the Graph Neural Operator (GNO) was chosen for this task.

---

confidence: 2
contribution: 3
initial_rating: 5
presentation: 3
questions: See Weaknesses section.
soundness: 3
strengths: The proposed method appears to be an interesting approach to model implied volatility surfaces for option pricing. The method through composite loss functions that include constraints such as no-arbitrage allows for direct training of a neural network architecture that directly provides smooth surfaces that directly satisfy these constraints. While I cannot comment directly how difficult/hands-on other alternative approaches are in this domain, the proposed model appears to have the flexibility to introduce other loss terms or vary the graph structure to incorporate other structural properties directly into the model. The model leads to lower errors compared to the SVI method.
summary: **Disclaimer**: I have no background on financial engineering so I cannot comment on the originality/significance aspect of this work. I might have missed essential points and only base my review on technical soundness and empirical validation

This paper introduces a method to nowcast implied volatility based on neural operators. This is achieved by a graph neural operator with MLPs as the lifting/projection layers and the kernel functions. The model uses a composite loss function that includes the fitting loss, no-arbitrage constraints, enforcing monotonicity and regularization terms to ensure smoothness. The method improves over the baseline SVI model.
weaknesses: The paper gives a good detail on the exact instantiation and implementation of the proposed method. I understand that there might be little prior work to base this paper on. However, it is unclear how the exact choices have been made and what would be their impact if changed. Note that the choices here might be purely made by domain knowledge, which I cannot judge and therefore might miss the specific considerations made here. However, I think the reader would benefit from more detail on the following: 

**In-Neighborhood sets**: The choice of in-heighborhood sets is arguably key model choice in this work. However, this work lacks detail on how the choice has been made and the subsampling heuristic has been devised. The authors argue that the choice is made because volatility smoothing requires limited global information exchange. I would then expect that the error of this method would go up if larger neighborhood sets are used. This paper would benefit from an additional analysis on changing the $\bar{p}$ and $K$ parameters. 

**Other model choices**: Again, there is little information about all the other model choices and/or instantiations. In particular, the choice of the kernel functions as MLP (as opposed to other choices) would be interesting for the reader. I'm also wondering whether the MLP can output negative values? I understand that there is a softplus operation in the architecture and ensures positive values but again the choice of a MLP would benefit from further explanation/ablation. 

**Loss function weighting coefficients:** I would kindly ask the authors to include more detail how the weighting of the loss function has been performed. As it is stated right now, there is no information on how to select/adjust the weightings on a new pr.oblem. It would help the dissemination the method if the authors would share their methodology.

**Ablations**: The claims in the paper could be strengthened by providing ablations In particular, different (wider) choices of the graph neighborhood would justify the choice of the authors. Additionally, omitting the no-arbitrage term and any of the other auxiliary terms would give more insight on the effect of the individual chosen components. 

I would reconsider my score if additional justification on the specific choices are included in the paper and the additional ablations on in-neighborhood and loss function weighting are performed.

Meta-review: The authors proposed a neural operator approach for implied volatility surface smoothing, and the trained model can directly map observed data to implied volatility surfaces, eliminating the need for calibration every time. Some reviewers pointed out missing baseline methods, such as the VAE approach, which was a flaw, but not very crucial. I recommended an acceptance as all reviewers left positive feedback in the end.

Decision: Accept (Poster)