1.4pt

# Pricing and calibration in the 4-factor path-dependent volatility model

Guido Gazzani University of Verona, Department of Economics, Via Cantarane 24,
37129 Verona, Italy<EMAIL_ADDRESS>Julien Guyon École des Ponts
ParisTech, Centre d’Enseignement et de Recherche en Mathématiques et Calcul
Scientifique (CERMICS), Bâtiment Coriolis, 6-8 avenue Blaise Pascal, 77455
Marne la Vallée Cedex 2, France<EMAIL_ADDRESS>
The authors acknowledge financial support from the BNP Paribas Chair _Futures
of Quantitative Finance_. The first author acknowledges financial support by
the “PHC AMADEUS” program (project number: 47561RJ), funded by the French
Ministry for Europe and Foreign Affairs, the French Ministry for Higher
Education and Research, and the Austrian Ministry for Higher Education. The
present work was initiated when the first author was affiliated to École des
Ponts ParisTech (CERMICS).

###### Abstract

We consider the path-dependent volatility (PDV) model of Guyon and Lekeufack
(2023), where the instantaneous volatility is a linear combination of a
weighted sum of past returns and the square root of a weighted sum of past
squared returns. We discuss the influence of an additional parameter that
unlocks enough volatility on the upside to reproduce the implied volatility
smiles of S&P 500 and VIX options. This PDV model, motivated by empirical
studies, comes with computational challenges, especially in relation to VIX
options pricing and calibration. We propose an accurate neural network
approximation of the VIX which leverages on the Markovianity of the 4-factor
version of the model. The VIX is learned as a function of the Markovian
factors and the model parameters. We use this approximation to tackle the
joint calibration of S&P 500 and VIX options.

Keywords: path-dependent volatility, calibration of financial models, neural
networks, S&P 500/VIX joint calibration
MSC (2020) Classification: 91B70, 91G20, 91G30, 91G60, 65C20.

##  1 Introduction

###  1.1 Path-dependent volatility models

In this article, we are interested in calibrating the path-dependent
volatility (PDV) model of Guyon and Lekeufack (2023) to S&P 500 (SPX) and VIX
futures and options. PDV models aim to provide an explicit and accurate
description of the purely endogenous part of the joint dynamics of an asset
price and its volatility (“spot-vol dynamics”): the volatility is described as
a deterministic, explicit function of past asset returns. This feedback effect
from past returns to volatility, which in turn impacts future returns, creates
very rich, intricate, nonlinear dynamics of asset prices.

PDV models have recently received increased attention for the following
reasons:

  * •

Empirical data confirms the path-dependent nature of volatility, see for
instance Zumbach (2009, 2010); Chicheportiche and Bouchaud (2014); Blanc et
al. (2017); Guyon and Lekeufack (2023); Andrès et al. (2023), as well as the
ARCH-GARCH literature.

  * •

PDV models can naturally reproduce all the following important stylized facts
about volatility: leverage effect, volatility clustering, heavy distribution
tails, roughness at the daily scale, weak and strong Zumbach effects, jump-
like behavior with very fast large volatility spikes followed by slower
decays.

Like local volatility models, PDV models are complete models, but they can
generate much more realistic joint spot-vol dynamics, as they capture the
path-dependent nature of volatility. Of course, volatility is not purely path-
dependent: some unanticipated, exogenous news have an impact on the volatility
of financial markets. Stochastic volatility (SV) models start from the
exogenous part, postulate a dynamics for the instantaneous volatility which
typically depends on unobservable factors, such as Ornstein-Uhlenbeck factors,
and can only generate some implicit, complicated path-dependency by
correlating the Brownian motions that drive the dynamics of the asset price
with those that drive the dynamics of the SV. Of course, one can build a fully
PDV model from an SV model by fully correlating those Brownian motions; see
for instance the quadratic rough Heston model by Gatheral et al. (2020) and
the path-dependent two-factor Bergomi model in Guyon (2022a) where it is shown
that the joint SPX/VIX smile calibration of the two-factor Bergomi model
degenerates the model into its path-dependent version. However, PDV is not
about fully correlating the Brownian motions of an underlying SV model. The
PDV paradigm is different. It suggests another, more natural approach to
volatility modeling:

  1. (i)

First, model the (large) purely endogenous part of volatility _explicitly_ as
well as possible, using PDV models that depend on _observable_ factors, such
as past asset returns.

  2. (ii)

Then, model the (smaller) exogenous part, for example by analyzing empirical
residuals of the (pure) PDV model.

We call the resulting model a _path-dependent stochastic volatility_ (PDSV)
model. The class of explicit PDV models includes for instance: many
contributions in the ARCH-GARCH literature, e.g., Sentana (1995), the complete
model of Hobson and Rogers (1998) also discussed by Foschi and Pascucci
(2005), the diagonal QARCH model of Chicheportiche and Bouchaud (2014), the
ZHawkes process of Blanc et al. (2017), the PDV model of Guyon and Lekeufack
(2023), as well as some version of the EWMA Heston of Parent (2023). In
particular, the PDV model of Guyon and Lekeufack (2023) was shown to have a
better predictive power of both implied volatility and future realized
volatility, across equity indexes and train/test periods, than the best
competing models available in the literature.

In the continuous-time version of the PDV model of Guyon and Lekeufack (2023),
the instantaneous volatility $\sigma_{t}$ takes a very simple form: it is an
affine combination of a weighted sum $R_{1}$ of past returns, and the square
root of a weighted sum $R_{2}$ of past squared returns:

$\sigma_{t}=\beta_{0}+\beta_{1}R_{1,t}+\beta_{2}\sqrt{R_{2,t}}.$

This particular, simple specification of a PDV came as the result of a
thorough empirical study where more complicated features and functional forms
(in particular, neural networks, including LSTM networks) were tested.

One distinctive feature of this PDV model is that it is homogeneous in
volatility, in particular capturing the dependency on past realized volatility
through the term $\sqrt{R_{2}}$. The other PDV models that have been studied
in the literature are either homogeneous in variance, or mix terms homogeneous
to a volatility with terms homogeneous to a variance, and/or lack the
important $R_{2}$ feature.

This distinctive feature may explain the higher predictive power of the model.
However, it also means that the model comes with computational challenges:
unlike other models which are chosen because of their mathematical
tractability, such as affine or polynomial models, it does not seem to have a
nice mathematical structure allowing for fast pricing. In particular, in this
PDV model, the forward instantaneous variances and the VIX are not directly
accessible, since the instantaneous variance $\sigma_{t}^{2}$ contains
$R_{1}\sqrt{R_{2}}$ and $\sqrt{R_{2}}$ terms whose conditional expectations
cannot be exactly computed.

The main objective of this article is to address those computational
challenges, with a view towards calibrating the model to SPX and VIX smiles.
We propose an accurate neural network approximation of the VIX which leverages
on the Markovianity of the 4-factor version of the model. This approximation
is then used to tackle the joint calibration of SPX and VIX options.

###  1.2 The joint calibration problem

Jointly calibrating to SPX and VIX futures and options is important to prevent
arbitrage and ensure accurate pricing of liquid hedging instruments. Indeed,
while calibrating to SPX options means incorporating market information on SPX
spot volatilities, calibrating to VIX derivatives means incorporating market
information on SPX _forward_ volatilities, an information that is not
contained in SPX option prices. For instance, Bourgey and Guyon (2024) show
the additional information contributed by VIX futures and VIX smiles, by
quantifying how model-free bounds for SPX path-dependent payoffs tighten when
the prices of VIX futures and VIX options are included. They also show that,
to avoid mispricing some payoffs (particularly forward-starting payoffs, which
are most sensitive to forward volatilities), it is important that the model
not only fits SPX options but also fits VIX options, even when the payoff
depends only on SPX prices.

In the rest of this section we review some of the literature on the joint
calibration problem. We refer the reader to Di Nunno et al. (2023) for a
recent survey on volatility modeling and advances on the joint calibration
problem. We split the literature in two main streams of research: exact
(nonparametric) and approximate (parametric) joint fits.

The first exact joint fit was obtained in Guyon (2020b), where a jointly
calibrated nonparametric discrete-time model was built by minimum entropy,
i.e., by solving a Schrödinger problem; see also Guyon (2024). Efficient
algorithms for this problem (in particular, faster than the Sinkhorn
algorithm) have recently been presented in Bourgey and Guyon (2024), who also
extended the model to continuous time by martingale interpolation. A direct
continuous-time Schrödinger bridge approach, inspired by Henry-Labordère
(2019), is presented in Guyon (2022b). A similar direct nonlinear optimal
transport approach was suggested in Guo et al. (2022). The direct continuous-
time approaches are much more computationally demanding than the novel
discrete-time-continuous-time procedure of Bourgey and Guyon (2024). Note that
an exact joint calibration ensures the absence of joint SPX/VIX arbitrages.

The other line of research concerns approximate fits of parametric models. The
PDV model of Guyon and Lekeufack (2023) falls into this category. The benefit
of low-parametric models is their interpretability. However, they cannot fit
market smiles exactly, since the low-parametric structure induces rigidities.
The main objective of a modeler is therefore to design parametric models that
are able to accurately fit market SPX and VIX smiles. As explained in Guyon
(2020a), the main challenge is to fit at the same time large short-term SPX
at-the-money (ATM) skews with relatively low VIX implied volatilities. In a
first attempt, Gatheral (2008) used a double constant elasticity of variance
model, which, despite its quite large number of parameters, cannot jointly fit
the implied volatilities of SPX and VIX options accurately enough for trading
purposes. Later on, several authors have included jumps in the dynamics of the
asset price and/or its volatility, e.g., the forward variance model of Cont
and Kokholm (2013) described as the exponential of an affine process with Lévy
jumps which allows for Fourier pricing, a regime-switching enhancement of the
classical Heston model by Papanicolaou and Sircar (2014), the 3/2 model with
jumps in the asset price of Baldeaux and Badran (2014), in the volatility
(Kokholm and Stisen (2015)), or with co-jumps and idiosyncratic jumps in the
volatility (Pacati et al. (2018)).

Continuous-paths models have also been employed to try to solve the joint
calibration problem. For instance, in Fouque and Saporito (2018), a Heston
model with stochastic vol-of-vol is calibrated, but only for maturities above
4 months when VIX options are less liquid. More recently, Rømer (2022)
considered a model where the volatility is driven by two Ornstein-Uhlenbeck
(OU) processes using a hyperbolic transformation function. Unfortunately, the
study of Rømer (2022) does not address the calibration to VIX futures, which
is crucial as VIX futures are the tradable assets on which VIX options are
written. Guyon and Lekeufack (2023) show that the 4-factor PDV (4FPDV) model,
a low-parametric Markovian model with only 10 parameters, is able to jointly
fit SPX smiles, VIX futures, and VIX smiles with good accuracy. Some
parametric models use a much large number of parameters, such as signature-
based models or neural SDE models. The use of two OU processes inspired the
first contribution with signature-based models, see Cuchiero et al. (2023b,
a). These models address the joint calibration problem by training only a
subset of model parameters, namely the linear read out map, as is usually done
in reservoir computing. The class of Gaussian polynomial processes introduced
in Abi Jaber et al. (2022), which can be seen as a special case of signature-
based models, lends itself to a quantization pricing technique for SPX and VIX
options. Using neural SDEs, Guyon and Mustapha (2023) show that the joint
calibration problem can be very accurately solved by a one-factor stochastic
local volatility model, provided enough flexibility is allowed in the
dynamics.

Rough volatility models, e.g., the quadratic rough Heston model (Gatheral et
al. (2020); Rosenbaum and Zhang (2022)) have also been used to tackle the
joint calibration problem with moderate success, due to their very
parsimonious, low-parametric nature. Rømer (2022) indeed shows that the rough
Bergomi model, the rough Heston model, and an extended rough Bergomi model are
all outperformed in the joint calibration task by a hyperbolic 2-OU model; see
also Abi Jaber and Li (2024) where different Bergomi-type models are compared
in this regard have been compared in terms of their calibration performances
on SPX options only. The rough Heston model belongs to the class of affine
Volterra processes considered in Abi Jaber et al. (2019); Cuchiero and
Teichmann (2019), and thus allows for Fourier pricing after solving the
corresponding Riccati equations. This underlying structure is the building
block of an extension with jumps (see Bondi et al. (2024a)) that is employed
in the context of the joint calibration by Bondi et al. (2024b), where a rough
Heston model with Hawkes-type jumps is shown, using Fourier pricing, to solve
the joint calibration problem for short maturities with good accuracy.

###  1.3 Main contributions

The main contributions of this article are the following.

  * •

The joint calibration procedure in Guyon and Lekeufack (2023) was manual, and
the optimization depended strongly on an initial guess. In this work we
present an automatic, more robust procedure.

  * •

Moreover, the calibration in Guyon and Lekeufack (2023) was very slow, since
the VIX was computed by nested Monte Carlo. In order to build a faster
calibration procedure, we provide a novel deep learning approach to simulating
the VIX under a given parametric Markov model, in our case the 4FPDV model, to
efficiently price VIX futures and VIX options.

  * •

With the addition of a single parameter we generate enough flexibility in the
4FPDV model to solve the joint SPX/VIX calibration problem for short
maturities, leveraging on the neural approximation of the VIX.

  * •

We also show that the 4FPDV model calibrates the SPX implied volatility
surface very well, over a large range of maturities (up to at least 1 year).

  * •

Finally, we provide an analysis of the stability over time of the calibrated
parameters.

The remainder of this article is structured as follows. Section 2 is a
reminder about the 4FPDV model. The novel deep learning estimation of the VIX
is presented in Section 3. Calibration results are presented in Section 4 for
the calibration to the SPX implied volatility surface only and in Section 5
for the joint SPX/VIX calibration. Finally, the stability analysis is carried
out in Section 6. The market data used in our numerical experiments is
provided by BNP Paribas via the Chair _Futures of Quantitative Finance._

##  2 The model

Throughout this article, $(r_{t})_{t\geq 0}$ and $(q_{t})_{t\geq 0}$ denote
the interest rate and the dividend yield curves, respectively, and are assumed
deterministic. The dividend yield includes the repo rate. We recall the
4-factor path-dependent volatility (4FPDV) model of Guyon and Lekeufack
(2023). The price $(S_{t})_{t\geq 0}$ of the asset is modeled under the unique
risk-neutral probability measure ${\mathbb{Q}}$ by

${\mathrm{d}}S_{t}=(r_{t}-q_{t})S_{t}\,{\mathrm{d}}t+S_{t}\sigma_{t}\,{\mathrm{d}}W_{t},$
(2.1)

where $S_{0}>0$, $\sigma=(\sigma_{t})_{t\geq 0}$ is the instantaneous
volatility process to be specified, and $W=(W_{t})_{t\geq 0}$ is a standard
one-dimensional ${\mathbb{Q}}$-Brownian motion. The volatility is path-
dependent and takes the following simple explicit form as a function of past
returns:

$\displaystyle\sigma_{t}$ $\displaystyle:=\sigma(R_{1,t},R_{2,t}),$
$\displaystyle\sigma(R_{1},R_{2})$
$\displaystyle:=\beta_{0}+\beta_{1}R_{1}+\beta_{2}\sqrt{R_{2}}+\beta_{1,2}R_{1}^{2}{\bf
1}_{\\{R_{1}>0\\}},$ (2.2) $\displaystyle R_{1,t}$
$\displaystyle:=(1-\theta_{1})R_{1,0,t}+\theta_{1}R_{1,1,t},$ $\displaystyle
R_{2,t}$ $\displaystyle:=(1-\theta_{2})R_{2,0,t}+\theta_{2}R_{2,1,t},$

where $R_{1,0,t}$, $R_{1,1,t}$, $R_{2,0,t}$, $R_{2,1,t}$ are four path-
dependent factors defined by111 By $(\frac{\mathrm{d}S_{u}}{S_{u}})^{2}$ we
mean $\sigma_{u}^{2}\,{\mathrm{d}}u$.

$R_{n,p,t}=\int_{-\infty}^{t}\lambda_{n,p}e^{-\lambda_{n,p}(t-u)}\left(\frac{\mathrm{d}S_{u}}{S_{u}}\right)^{n},\qquad
n\in\\{1,2\\},\quad p\in\\{0,1\\}.$

The two factors $R_{1,p,t}$ are exponentially weighted moving averages of past
returns, measuring the recent trend in the asset price, while the two factors
$R_{2,p,t}$ are exponentially weighted moving averages of past squared
returns, measuring the recent volatility in the asset price, with exponential
weights $\lambda_{n,p}>0$. Note that

$R_{n,t}=\int_{-\infty}^{t}K_{n}(t-u)\left(\frac{\mathrm{d}S_{u}}{S_{u}}\right)^{n},\qquad
n\in\\{1,2\\},$

is a weighted average of past returns (when $n=1$) or past squared returns
(when $n=2$), where the convolution kernels

$K_{n}(t):=(1-\theta_{n})\lambda_{n,0}e^{-\lambda_{n,0}t}+\theta_{n}\lambda_{n,1}e^{-\lambda_{n,1}t}$
(2.3)

are convex combinations of two exponential kernels. This parametrization of
the kernels allows us to

  1. (i)

mix short memory (large $\lambda_{n,0}$) and long memory (smaller
$\lambda_{n,1}$)—in particular, such kernels are well approximated by power
laws over a large range of maturities;

  2. (ii)

build a Markovian model, since

$\displaystyle\mathrm{d}R_{1,p,t}$
$\displaystyle=-\lambda_{1,p}R_{1,p,t}\,{\mathrm{d}}t+\lambda_{1,p}\sigma(R_{1,t},R_{2,t})\,{\mathrm{d}}W_{t},$
$\displaystyle\mathrm{d}R_{2,p,t}$
$\displaystyle=\lambda_{2,p}(\sigma(R_{1,t},R_{2,t})^{2}-R_{2,p,t})\,{\mathrm{d}}t.$

The parameters have the following interpretation:

  * •

$\lambda_{1,0}$ captures the dependence of $R_{1}$ on recent returns; the
larger $\lambda_{1,0}$, the more weight is given to recent returns;

  * •

$\lambda_{1,1}<\lambda_{1,0}$ captures the dependence of $R_{1}$ on older
returns; the smaller $\lambda_{1,1}$, the more weight is given to old returns;

  * •

$\lambda_{2,0}$ captures the dependence of $R_{2}$ on recent squared returns;
the larger $\lambda_{2,0}$, the more weight is given to recent squared
returns;

  * •

$\lambda_{2,1}<\lambda_{2,0}$ captures the dependence of $R_{2}$ on older
squared returns; the smaller $\lambda_{2,1}$, the more weight is given to old
squared returns;

  * •

$\theta_{1}\in[0,1]$ (resp., $\theta_{2}\in[0,1]$) mixes the dependence on
recent and older returns (resp., squared returns) to form the summary variable
$R_{1}$ (resp., $R_{2}$);

  * •

$\beta_{0}>0$ is the baseline instantaneous volatility;

  * •

$\beta_{1}<0$ and $\beta_{2}>0$ are the sensitivities of the volatility to the
trend $R_{1}$ and the historical volatility $\sqrt{R_{2}}$, respectively;

  * •

$\beta_{1,2}\geq 0$ is an additional parameter which produces volatility when
the asset price trends up—see below.

#### The role of $\beta_{1,2}$.

The original 4FPDV model, which aimed to explain a single global level of
implied volatility (e.g., the VIX) using past daily returns, has
$\beta_{1,2}=0$. As explained in Guyon and Lekeufack (2023), allowing
$\beta_{1,2}$ to vary did not significantly increase the predictive power of
the model. However, the continuous-time limit of the original 4FPDV model,
that is, Model (2) with $\beta_{1,2}=0$, generates out-the-money (OTM) call
implied volatilities that are too small and do not increase with strike for
large strikes. In order for the continuous-time model to fit not only one
single global level of implied volatility but full smiles of implied
volatility, Guyon and Lekeufack (2023) added the term
$\beta_{1,2}R_{1}^{2}{\bf 1}_{\\{R_{1}>0\\}}$, with $\beta_{1,2}\geq 0$.
Indeed, compared to the case where $\beta_{1,2}=0$, this term generates added
volatility when the asset price experiences a positive trend. This allows the
asset price to take large values with higher probability, and this results in
increased risk-neutral density for large asset price values, therefore
increased implied volatility for large strikes.

Table 1:

$\lambda_{1,0}=62.11$ | $\lambda_{1,1}=32.25$ | $\theta_{1}=0.23$ | $\lambda_{2,0}=9.57$ | $\lambda_{2,1}=3.51$ | $\theta_{2}=0.99$
---|---|---|---|---|---

$\beta_{0}=0.026$ | $\beta_{1}=-0.138$ | $\beta_{2}=0.69$
---|---|---

Figure 1: Impact of $\beta_{1,2}$ on model SPX implied volatilities at
maturity 16 days (top left), the corresponding model SPX density (top right),
and model VIX futures and VIX implied volatilities (bottom). Model parameters
are reported in Table 1, except $\beta_{1,2}$ which varies in
$\\{0,0.05,0.1,0.15\\}$. Prices of call options on $S$ are computed by Monte
Carlo with variance reduction (see Section 4) using $N_{\text{MC}}=10^{5}$
trajectories and a discretization step $\Delta t=\frac{1}{2520}$. VIX call
option prices are computed using $3\cdot 10^{3}$ nested Monte Carlo paths.
Initial values of the factors: $R_{1,0,0}=0.2988$, $R_{1,1,0}=0.2397$,
$R_{2,0,0}=0.016$, $R_{2,1,0}=0.02$.

This is illustrated in Figure 1 (top), where we plot the model SPX smile and
risk-neutral distribution at maturity 16 days. We use the model parameters
from Table 8 in Guyon and Lekeufack (2023), calibrated to market smiles,
except that we vary $\beta_{1,2}\in\\{0,0.05,0.1,0.15\\}$; the parameters are
reported in Table 1. Figure 1 (top left) shows that the extra term
$\beta_{1,2}R_{1}^{2}{\bf 1}_{\\{R_{1}>0\\}}$ allows us to reproduce the
hockey-stick shape that is typical of market SPX smiles; Figure 1 (top right)
illustrates the increased risk-neutral density on the large SPX side. Figure 1
(bottom) shows the model VIX smile with maturity 14 days. The increased VIX
future value, denoted by a dashed line, results from the additional
volatility.

The dynamics of the instantaneous volatility in model (2) is given by the
following lemma. For ease of notation, we denote
$R_{t}:=(R_{1,0,t},R_{1,1,t},R_{2,0,t},R_{2,1,t})$ and introduce the following
quantities:

$\bar{\lambda}_{n}:=(1-\theta_{n})\lambda_{n,0}+\theta_{n}\lambda_{n,1},\quad\bar{R}_{n,t}:=\frac{(1-\theta_{n})\lambda_{n,0}R_{n,0,t}+\theta_{n}\lambda_{n,1}R_{n,1,t}}{\bar{\lambda}_{n}},\quad
n\in\\{1,2\\}.$

######  Lemma 2.1.

Let $\sigma=(\sigma_{t})_{t\geq 0}$ satisfy (2). Then

$\mathrm{d}\sigma_{t}=\mu(R_{t})\,{\mathrm{d}}t+\nu(R_{t})\,{\mathrm{d}}W_{t},$

where

$\displaystyle\mu(R_{t})$ $\displaystyle=-(\beta_{1}+2\beta_{1,2}R_{1,t}{\bf
1}_{\\{R_{1,t}>0\\}})\bar{\lambda}_{1}\bar{R}_{1,t}+\beta_{2}\bar{\lambda}_{2}\frac{(\sigma_{t}^{2}-\bar{R}_{2,t})}{2\sqrt{R_{2,t}}}+\beta_{1,2}\bar{\lambda}_{1}^{2}\sigma_{t}^{2},$
$\displaystyle\nu(R_{t})$ $\displaystyle=(\beta_{1}+2\beta_{1,2}R_{1,t}{\bf
1}_{\\{R_{1,t}>0\\}})\bar{\lambda}_{1}\sigma_{t}.$

###### Proof.

The proof is a straightforward application of Itô’s formula. ∎

When $\beta_{1,2}=0$, the instantaneous (lognormal) volatility of the
instantaneous volatility is constant (equal to
$|\beta_{1}|\bar{\lambda}_{1}$), like in Bergomi models. When $\beta_{1,2}>0$,
it is equal to $|(\beta_{1}+2\beta_{1,2}R_{1,t}{\bf
1}_{\\{R_{1,t}>0\\}})|\bar{\lambda}_{1}$ and thus depends on the trend
$R_{1,t}$. It vanishes when $R_{1,t}=-\frac{\beta_{1}}{2\beta_{1,2}}>0$, and
can get large when the asset price trends very positively.

##  3 Learning the VIX with neural networks

Our goal in this article is to build a computer procedure to jointly calibrate
the 4FPDV model to SPX smiles, VIX futures, and VIX smiles. To this end, we
first need to be able to fastly compute the VIX in the model at any date in
the future. In view of calibrating the model to VIX futures and VIX smiles, we
will not only learn VIXT as a function of the Markovian factors
$(R_{1,0,T},R_{1,1,T},R_{2,0,T},R_{2,1,T})$, but also as a function of the
model parameters.

###  3.1 A brief reminder on the VIX

The CBOE Volatility Index, also known as the VIX, is a popular measure of the
market’s expected volatility of the SPX index over the next 30 days. It is
calculated and published by the Chicago Board Options Exchange (CBOE). The VIX
is meant to represent the 30-day implied volatility of the log-contract on the
SPX index. We thus take as stylized definition of the VIX

${\textrm{VIX}}_{T}^{2}:=\text{Price}_{T}\left[-\frac{2}{\Delta}\log\left(\frac{S_{T+\Delta}}{F_{T}^{T+\Delta}}\right)\right],$
(3.1)

where $\Delta=30$ days and $F_{t}^{T+\Delta}$ denotes the price at $t$ of the
SPX future with maturity $T+\Delta$. Within a pricing model, we therefore
define the VIX by

${\textrm{VIX}}_{T}^{2}:=\mathbb{E}\left[-\frac{2}{\Delta}\log\left(\frac{S_{T+\Delta}}{F_{T}^{T+\Delta}}\right)\middle|\mathcal{F}_{T}\right],$
(3.2)

where ${\mathbb{E}}[\cdot]$ denotes the expectation under the pricing measure
${\mathbb{Q}}$ and $(\mathcal{F}_{t})$ denotes the filtration representing the
flow of available information in the market. Note that the nonnegativity of
the r.h.s. of (3.1) results from absence of arbitrage, while that of the
r.h.s. of (3.2) results from the martingale property of
$(F_{t}^{T+\Delta})_{0\leq t\leq T+\Delta}$, the concavity of the logarithm,
and Jensen’s inequality.

In order to compute the VIX in model (2), we thus need to compute the
conditional expectation (3.2), where $(\mathcal{F}_{t})$ is the filtration
generated by the Brownian motion $W$. From Neuberger (1994); Dupire (1994),
since the SPX model price has continuous paths,

${\textrm{VIX}}_{T}^{2}=\mathbb{E}\left[\frac{1}{\Delta}\int_{T}^{T+\Delta}\sigma_{t}^{2}\,{\mathrm{d}}t\middle|\mathcal{F}_{T}\right].$
(3.3)

This equality, which immediately follows from (3.2) by applying the Itô
formula to $\log(F_{t}^{T+\Delta})$ between $T$ and $T+\Delta$, explains the
choice (3.1) to represent the expected volatility of the SPX index over the
next 30 days. We observe that when $\sigma_{t}$ satisfies (2), the expression
(3.3) is not known in closed form, i.e., one cannot explicitly compute

${\textrm{VIX}}_{T}^{2}=\frac{1}{\Delta}\mathbb{E}\left[\int_{T}^{T+\Delta}\big{(}\beta_{0}+\beta_{1}R_{1,t}+\beta_{2}\sqrt{R_{2,t}}+\beta_{1,2}R_{1,t}^{2}{\bf
1}_{\\{R_{1,t}>0\\}}\big{)}^{2}\,{\mathrm{d}}t\middle|\mathcal{F}_{T}\right].$
(3.4)

Nevertheless, due to the Markovianity of
$R_{t}=(R_{1,0,t},R_{1,1,t},R_{2,0,t},R_{2,1,t})$ and the time-homogeneity of
its dynamics, the VIX squared at time $T>0$ is a measurable function
$f(\Theta,R_{T})$, where $\Theta$ denotes the collection of the model
parameters,

$\Theta=(\lambda_{1,0},\lambda_{1,1},\theta_{1},\lambda_{2,0},\lambda_{2,1},\theta_{2},\beta_{0},\beta_{1},\beta_{2},\beta_{1,2})\in\mathbb{R}^{10}.$

Our key idea is to learn the function $f$ by parameterizing it as a feed-
forward neural network.

###  3.2 The learning procedure

The steps of the learning procedure are outlined as follows:

  1. (i)

For each of the 10 model parameters $\Theta_{i}$, sample $N>0$ independent
configurations $\Theta_{i}(\omega_{j})$ accounting for the natural constraint
$\Theta_{i}\in[a_{i},b_{i}]$:

$\forall\ i\in\\{1,\dots,10\\},\ \forall\
j\in\\{1,\dots,N\\},\qquad\Theta_{i}(\omega_{j})\sim{\mathcal{U}}[a_{i},b_{i}],$

where ${\mathcal{U}}[a,b]$ denotes the uniform distribution over the interval
$[a,b]$.

  2. (ii)

For each parameter combination $j\in\\{1,\dots,N\\}$, simulate one random
realization of the vector
$R_{t}(\omega_{j}):=(R_{1,0,t}^{j},R_{1,1,t}^{j},R_{2,0,t}^{j},R_{2,1,t}^{j})$
in the 4FPDV model (2) for a grid of future dates
$t\in\\{t_{1}(\omega_{j}),\dots,t_{M}(\omega_{j})\\}$ to be specified later.

  3. (iii)

Also store the corresponding $\widehat{\textrm{VIX}}_{t}(\omega_{j})$ for
$t\in\\{t_{1}(\omega_{j}),\dots,t_{M}(\omega_{j})\\}$ and
$j\in\\{1,\dots,N\\}$, where $\widehat{{\textrm{VIX}}}$ denotes the square
root of the random variable (3.4) estimated using nested Monte Carlo.

  4. (iv)

Randomly split the total number of parameter configurations into
$N_{1}+N_{2}=N$. The first $N_{1}$ samples consist of the training set, while
the next $N_{2}$ samples form the validation set.

  5. (v)

For a limited number of epochs, train a neural network
${\mathcal{N}}{\mathcal{N}}:(\Theta,R)\in{\mathbb{R}}^{14}\mapsto\textrm{VIX}\in{\mathbb{R}}^{+}$
on $N_{1}$ while tuning the architecture based on the validation set $N_{2}$
using a cross-validation approach, e.g., Random Search, Bayesian Optimization,
or Hyperband; see Li et al. (2017). We use as loss function the root mean
squared error (RMSE), i.e., the square root of

$\frac{1}{NM}\sum_{j=1}^{N}\sum_{k=1}^{M}\left(\widehat{\textrm{VIX}}_{t_{k}(\omega_{j})}(\omega_{j})-\mathcal{N}\mathcal{N}(\Theta(\omega_{j}),R_{t_{k}(\omega_{j})}(\omega_{j}))\right)^{2}.$

The inputs of the network are standardized to help the training process.

######  Remark 3.1.

We did not consider the loss function

$\frac{1}{NM}\sum_{j=1}^{N}\sum_{k=1}^{M}\left(\frac{1}{\Delta}\int_{t_{k}(\omega_{j})}^{t_{k}(\omega_{j})+\Delta}\sigma_{t}^{2}(\omega_{j})\,{\mathrm{d}}t-\mathcal{N}\mathcal{N}(\Theta(\omega_{j}),R_{t_{k}(\omega_{j})}(\omega_{j}))\right)^{2}$

(which does not require simulating nested paths) to learn the VIX squared due
to the very large conditional variance of the integrated instantaneous
variance. Using averages $\widehat{\textrm{VIX}}_{t}(\omega_{j})$ improves the
learning, as analyzed in Alfonsi et al. (2023).

Note that this methodology naturally extends to any parametric Markovian
volatility model. The above routine is first used to select the optimal
architecture for the neural network. The training further continues to achieve
higher accuracy. In the next sections, we report our numerical results.

###  3.3 Training and testing of the network

We consider $N_{1}=6.8\cdot 10^{5}$, $N_{2}=1.2\cdot 10^{5}$ configurations of
the model parameters to train and validate the network, respectively. Natural
constraints yielding enough flexibility for the volatility smiles generated by
the 4FPDV model (2) are given by

$\displaystyle\beta_{0}$
$\displaystyle\in[0,0.2],\quad\beta_{1}\in[-0.25,0),\quad\beta_{2}\in[0,1),\quad\beta_{1,2}\in[0,0.3],$
$\displaystyle\lambda_{n,p}$
$\displaystyle\in[1,100],\quad\theta_{n}\in[0,1],\qquad n\in\\{1,2\\},\quad
p\in\\{0,1\\}.$

We require that $\lambda_{n,0}>\lambda_{n,1}$ so that $\lambda_{n,0}$ encodes
the short memory and $\lambda_{n,1}$ the long memory. Notice that possible
combinations of the above may give rise to extremely large values of the
instantaneous volatility of the instantaneous volatility, hence in the
simulation of the training set we remove the combinations for which
$|\beta_{1}|\bar{\lambda}_{1}>10$. We observe the four factors
$R_{t}(\omega_{j})$ and the VIX, namely
$\widehat{\textrm{VIX}}_{t}(\omega_{j})$, at $M=200$ evenly spaced time points
$t_{1}(\omega_{j})<\cdots<t_{M}(\omega_{j})$, with
$t_{1}(\omega_{j})=\max\\{1/\lambda_{1,0}(\omega_{j}),1/\lambda_{2,0}(\omega_{j})\\}$
and $t_{M}(\omega_{j})=1$. This choice of $t_{1}(\omega_{j})$ ensures that the
four factors are observed close to or in their stationary states. We compute
$\widehat{{\textrm{VIX}}}_{t}$ by nested Monte Carlo: $10^{4}$ nested paths
are used to simulate the conditional expectation in (3.4), where the integral
is approximated with Riemann sums. Both processes $\sigma$ and $R$ are sampled
with discretization step $\Delta t=\frac{1}{2520}$.

To assess the goodness of the training, we consider two tests. First, we build
the histogram of the average absolute error over 10,000 different Monte Carlo
samples between the optimal neural network approximation
${\mathcal{N}}{\mathcal{N}}^{\star}$ and $\widehat{{\textrm{VIX}}}$ for 1,000
random configurations of the model parameters. In Figure 2, we observe that
for 99% of the tested configurations, the mean absolute error is smaller than
0.55 (around half a volatility point).

Figure 2: Histogram of the average absolute error between
${\mathcal{N}}{\mathcal{N}}^{\star}$ and $\widehat{{\textrm{VIX}}}$ for
$10^{3}$ configurations of parameters of the model. The average is taken over
$10^{4}$ Monte Carlo samples.

Some of the 1,000 configurations of model parameters tested in Figure 2
generate unrealistic VIX smiles and may produce the largest average absolute
errors. Therefore, in our second test, we compare the VIX implied volatility
smile computed using the neural approximation
${\mathcal{N}}{\mathcal{N}}^{\star}$ of the VIX with the one computed using
$\widehat{{\textrm{VIX}}}$ for a few sets of market-calibrated parameters.
Figure 3 shows an excellent agreement between the two VIX smiles and the two
VIX futures for one set of realistic parameters.

Figure 3: Model parameters: $\lambda_{1,0}=29$, $\lambda_{1,1}=20$,
$\theta_{1}=0.69$, $\lambda_{2,0}=81$, $\lambda_{2,1}=66$, $\theta_{2}=0.25$,
$\beta_{0}=0.11$, $\beta_{1}=-0.057$, $\beta_{2}=0.1$, $\beta_{1,2}=0.256$.
Top left: VIX smiles and VIX future (dashed) with maturity 16 days. In blue:
the VIX is computed via the pre-trained neural network
${\mathcal{N}}{\mathcal{N}}^{\star}$. In red: the VIX is computed via nested
Monte Carlo, $\widehat{{\textrm{VIX}}}$. Top right: A total of $10^{6}$
predictions of ${\mathcal{N}}{\mathcal{N}}^{\star}$ vs.
$\widehat{{\textrm{VIX}}}$. Bottom: histogram of the absolute error between
${\mathcal{N}}{\mathcal{N}}^{\star}$ and $\widehat{{\textrm{VIX}}}$.

In Figure 3 we additionally plot the predictions of the neural network against
the nested Monte Carlo benchmark as well as the histogram of the absolute
error between the two. For this particular set of realistic parameters, the
average absolute error is 0.2, which is in the bulk of the distribution of
average absolute errors plotted in Figure 2. For 99% of the $10^{6}$ Monte
Carlo samples, the (pathwise) absolute error between
${\mathcal{N}}{\mathcal{N}}^{\star}$ and $\widehat{{\textrm{VIX}}}$ is smaller
than 0.65.

In Appendix A, we report similar results for parameters jointly calibrated to
the SPX and VIX smiles on June 2, 2021, and June 3, 2021. This informs us that
for market-calibrated parameters the optimal neural network accurately learns
the VIX computed with nested Monte Carlo and can be used to price VIX futures
and VIX options. In any case, when we calibrate model parameters to VIX market
data using the ${\mathcal{N}}{\mathcal{N}}^{\star}$ approximation of the VIX,
we systematically check that the VIX smile computed using
${\mathcal{N}}{\mathcal{N}}^{\star}$ and the one computed using the nested
estimator $\widehat{{\textrm{VIX}}}$ match each other closely (see Figures 9
and 10).

###  3.4 The architecture of the neural network

N∈N Figure 4: Example of architecture for a feed-forward neural network.

Here, we report the architecture of the network obtained after the
hyperparameter tuning via Bayesian optimization. We asked the KerasTuner to
select from 1 to 5 dense hidden layers for a feed-forward neural network, with
a number of nodes ranging from 64 to 512 with 32 as step size, and two
possible activation functions, namely

$\text{ReLU}(x):=\max\\{x,0\\},\quad\text{tanh}(x):=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}},\qquad\forall
x\in{\mathbb{R}}.$

We report the optimal values in Table 2. Note that the optimal architecture
uses the full 5 hidden layers, with varying numbers of nodes per layer, and
that it mixes both activation functions. Moreover, we tuned the learning rate
of the Adam optimizer between $10^{-5}$ and $10^{-2}$ with log-sampling. The
optimal learning rate is $4.2\cdot 10^{-5}$.

Table 2:  Architecture of our neural network Layer type | Activation | No. Nodes
---|---|---
Input | ReLU | 14
Hidden 1 | tanh | 448
Hidden 2 | tanh | 64
Hidden 3 | ReLU | 224
Hidden 4 | tanh | 416
Hidden 5 | ReLU | 128
Output | linear | 1

###  3.5 Pricing and calibration of VIX futures and VIX options

Pricing. Suppose we are given a set ${\mathcal{T}}^{{\textrm{VIX}}}$ of VIX
option maturities. For each $T\in{\mathcal{T}}^{{\textrm{VIX}}}$, we consider
a collection ${\mathcal{K}}_{T}^{{\textrm{VIX}}}$ of strike prices. Using the
classical Monte Carlo method and the neural approximation
${\mathcal{N}}{\mathcal{N}}^{\star}$ of the VIX under model (2) with
parameters $\Theta$, the price of the VIX future and the price of a VIX call
are respectively approximated by

$\displaystyle F_{\Theta}^{{\textrm{VIX}}}(T)$
$\displaystyle=\frac{1}{N_{MC}}\sum_{j=1}^{N_{MC}}{\mathcal{N}}{\mathcal{N}}^{\star}(\Theta,R_{T}(\omega_{j})),$
$\displaystyle C_{\Theta}^{{\textrm{VIX}}}(T,K)$
$\displaystyle=e^{-\int_{0}^{T}r_{s}\,{\mathrm{d}}s}\frac{1}{N_{MC}}\sum_{j=1}^{N_{MC}}({\mathcal{N}}{\mathcal{N}}^{\star}(\Theta,R_{T}(\omega_{j}))-K)^{+},$

for all $T\in{\mathcal{T}}^{{\textrm{VIX}}}$ and for all
$K\in{\mathcal{K}}_{T}^{{\textrm{VIX}}}$.

Calibration. We now turn to the problem of calibrating the model parameters to
VIX futures and VIX options. Notice that since the actual tradable underlying
of VIX options are the VIX futures with the same expiry as the option, it is
crucial that the VIX futures be very well calibrated, as already argued in,
e.g., Pacati et al. (2018); Guyon (2020a); Cuchiero et al. (2023a). To be
consistent, the model VIX implied volatilities are then of course computed by
inverting the Black formula using the corresponding model VIX future.

Since we have learned the VIX not only as a function of the four Markovian
factors but also as a function of the 10 model parameters, $\Theta$, we can
use the loss function

$\displaystyle L_{{\textrm{VIX}}}(\Theta)$
$\displaystyle=\omega_{F}\frac{1}{\\#{\mathcal{T}}^{{\textrm{VIX}}}}\sum_{T\in{\mathcal{T}}^{{\textrm{VIX}}}}\ell\left(F_{\Theta}^{{\textrm{VIX}}}(T),F^{{\textrm{VIX}}}(T)\right)$
(3.5)
$\displaystyle\quad+\omega_{{\textrm{VIX}}}\frac{1}{\\#{\mathcal{T}}^{{\textrm{VIX}}}}\sum_{T\in{\mathcal{T}}^{{\textrm{VIX}}}}\frac{1}{\\#{\mathcal{K}}_{T}^{{\textrm{VIX}}}}\sum_{K\in{\mathcal{K}}_{T}^{{\textrm{VIX}}}}\gamma_{T,K}^{{\textrm{VIX}}}\ell\left(C_{\Theta}^{{\textrm{VIX}}}(T,K),C^{{\textrm{VIX}}}(T,K)\right)$

to calibrate the model to VIX futures and VIX options, where

  * •

$\omega_{F},\ \omega_{{\textrm{VIX}}}>0$ are hyperparameters weights;

  * •

$\\#$ denotes the cardinality of a set;

  * •

$(\gamma_{T,K}^{{\textrm{VIX}}})_{T,K}$ are normalized vega weights; given
${\mathcal{K}}(T)$ the strikes available for maturity $T>0$,

$\gamma_{T,K}:=\frac{{\mathcal{V}}_{T,K}}{\sum_{K^{\prime}\in{\mathcal{K}}(T)}{\mathcal{V}}_{T,K^{\prime}}},$

where ${\mathcal{V}}_{T,K}$ denotes the Black vega of the call option with
maturity $T$ and strike $K$, computed using the corresponding market mid
implied volatility.

  * •

$F^{{\textrm{VIX}}},C^{{\textrm{VIX}}}$ denote the market futures and the
market call option prices, respectively;

  * •

$\ell:\mathbb{R}^{+}\times\mathbb{R}^{+}\to\mathbb{R}^{+}$ is a score
function; we pick

$\ell(x,y)=\left(\frac{x}{y}-1\right)^{2}.$ (3.6)

Note that we calibrate the model to VIX option prices, not VIX implied
volatilities, as model VIX implied volatilities depend on model VIX futures;
this prevents a mismatch between model and market VIX futures to cascade into
a mismatch between model and market VIX implied volatilities. The weights
$(\gamma_{T,K}^{{\textrm{VIX}}})_{T,K}$ are often taken to be the inverse of
the market bid-ask spreads, which gives more weight to those options where the
bid-ask spread is tighter (see, e.g., Cont and Ben Hamida (2005) for a
discussion on this topic). In Figure 5 we compare the normalized vega weights
with the inverse of the market bid-ask spreads as of June 2, 2021, for VIX
options with 33 days to expiration. The vega weights are a smoothed version of
the inverse of the bid-ask spread.

Figure 5: In blue: normalized vega weights; in red: normalized bid-ask spread
weights. Data as of June 2, 2021.

######  Remark 3.2.

An important remark is that we use neural networks to directly learn the VIX
at future dates, pathwise, as a function of the model parameters (and the
Markovian factors)—not the prices at time 0 of VIX futures and VIX options.
Our use of neural networks thus differs from the so-called “deep pricing” and
“deep calibration” approaches. Let $\mathcal{T}$ and $\mathcal{K}$ be
collections of maturities and strikes, respectively. Using the notation
$\phi_{T,K}$ to represent the price of a European vanilla option or its
implied volatility with maturity $T\in\mathcal{T}$ and strike
$K\in\mathcal{K}$, the methodologies can be summarized in three network
parametrizations:

  * •

Global deep pricing:
${\mathcal{N}}{\mathcal{N}}:\Theta\mapsto(\phi_{T,K}^{\text{model}})_{T\in\mathcal{T},K\in\mathcal{K}}$
as in Horvath et al. (2021); Rosenbaum and Zhang (2022); Rømer (2022), for
fixed $\mathcal{T}$ and $\mathcal{K}$.

  * •

Pointwise deep pricing:
${\mathcal{N}}{\mathcal{N}}:(\Theta,T,K)\mapsto\phi_{T,K}^{\text{model}}$ as
in Bayer and Stemper (2018); Baschetti et al. (2024).

  * •

Deep calibration:
${\mathcal{N}}{\mathcal{N}}:(\phi_{T,K}^{\text{model}})_{T\in\mathcal{T},K\in\mathcal{K}}\mapsto\Theta$
as for instance in Hernandez (2016); Gambara and Teichmann (2020) to directly
solve the calibration problem.

##  4 Calibration to SPX options

###  4.1 Pricing and calibration of SPX options

Suppose we are given a set ${\mathcal{T}}^{{\textrm{SPX}}}$ of SPX option
maturities. Then, for each $T\in{\mathcal{T}}^{{\textrm{SPX}}}$, we consider a
collection ${\mathcal{K}}_{T}^{{\textrm{SPX}}}$ of strike prices. Under model
(2.1) with parameters $\Theta$, the price of an SPX call or put option is
approximated by

$\displaystyle C_{\Theta}^{{\textrm{SPX}}}(T,K)$
$\displaystyle=e^{-\int_{0}^{T}r_{s}\,{\mathrm{d}}s}\frac{1}{N_{MC}}\sum_{j=1}^{N_{MC}}(S_{T}(\omega_{j})-K)^{+},$
(4.1) $\displaystyle P_{\Theta}^{{\textrm{SPX}}}(T,K)$
$\displaystyle=e^{-\int_{0}^{T}r_{s}\,{\mathrm{d}}s}\frac{1}{N_{MC}}\sum_{j=1}^{N_{MC}}(K-S_{T}(\omega_{j}))^{+},$
(4.2)

for all $T\in{\mathcal{T}}^{{\textrm{SPX}}}$ and for all
$K\in{\mathcal{K}}_{T}^{{\textrm{SPX}}}$ respectively. Denote by
$\sigma_{\text{IV},\Theta}^{{\textrm{SPX}}}(T,K)$ the model implied volatility
under the parameters $\Theta$ and for a fixed maturity and strike price,
$T,K>0$ computed for OTM options. Since by construction the model and market
SPX futures curve agree we can consider a loss function which aims to measure
the discrepancy between either option prices (see, e.g., Cont and Ben Hamida
(2005); Papanicolaou and Sircar (2014); Cuchiero et al. (2023a)) or implied
volatilities (as, e.g., in Pacati et al. (2018); Abi Jaber et al. (2022);
Bondi et al. (2024b); Guyon and Mustapha (2023)). In the following we consider
the latter approach and introduce a loss function $L_{{\textrm{SPX}}}$ which
reads as follows

$L_{{\textrm{SPX}}}(\Theta)=\omega_{{\textrm{SPX}}}\frac{1}{\\#{\mathcal{T}}^{{\textrm{SPX}}}}\sum_{T\in{\mathcal{T}}^{{\textrm{SPX}}}}\frac{1}{\\#{\mathcal{K}}_{T}^{{\textrm{SPX}}}}\sum_{K\in{\mathcal{K}}_{T}^{{\textrm{SPX}}}}\ell\left(\sigma_{\text{IV},\Theta}^{{\textrm{SPX}}}(T,K),\sigma_{\text{IV}}^{{\textrm{SPX}}}(T,K)\right),$
(4.3)

where $\omega_{{\textrm{SPX}}}>0$ is a hyperparameter weight and $\ell$ is the
score function (3.6).

###  4.2 Calibration to the SPX surface

In Figures 6 and 7 we report the results of two calibration exercises: as of
June 3, 2021, we calibrate to monthly SPX options with maturities between 15
and 351 days; and as of October 25, 2023, to monthly SPX options with
maturities from 23 days to 296 days. At each iteration of the optimizer we
simulate $N_{MC}=2\cdot 10^{5}$ trajectories with discretization time step
$\Delta t=\frac{1}{504}$. The calibrated parameters of the path-dependent
volatility model are reported in Tables 3 and 4, respectively. With only 10
parameters, the model is able to accurately fit smiles for both short and long
maturities. In particular, the two-exponential kernels allow us to accurately
capture the market term-structure of the ATM skew (see Figure 8). We observe
that two different mean reversion speeds need to be coupled, see the values of
the $\lambda$s and $\theta$s in Table 3.

Figure 6:  Calibration as of June 3, 2021. Left: in blue dots, the calibrated
implied volatility smiles under the PDV model; in red stars, the market mid
implied volatility smiles. Right: the interpolated relative error between the
calibrated implied volatility smiles and the market mid ones displayed on the
left. We report the initial values of the factors given the calibrated
parameters in the Table 3. Recall that these are uniquely determined by the
calibrated $\lambda$s and the observed past daily returns for June 3, 2021
with a cut-off at 1000 days: $R_{1,0}=0.0894$, $R_{1,1}=-0.1602$,
$R_{2,0}=0.0031$, $R_{2,1}=0.0476$.  Table 3:  Calibrated parameters of the
4FPDV model. Calibration to the SPX surface as of June 3, 2021.

$\lambda_{1,0}=34.39$ | $\lambda_{1,1}=13.26$ | $\theta_{1}=0.501$ | $\lambda_{2,0}=95.63$ | $\lambda_{2,1}=1.428$ | $\theta_{2}=0.448$
---|---|---|---|---|---

$\beta_{0}=0.0493$ | $\beta_{1}=-0.1999$ | $\beta_{2}=0.5479$ | $\beta_{1,2}=0.2285$
---|---|---|---

Figure 7:  Calibration as of October 25, 2023. See caption of Figure 6 for
color reference. The initial values of the factors with a cut-off at 1000 days
are given by: $R_{1,0}=-1.2054$, $R_{1,1}=-0.2428$, $R_{2,0}=0.0178$,
$R_{2,1}=0.0166$. Table 4:  Calibrated parameters of the 4FPDV model.
Calibration to the SPX surface as of October 25, 2023.

$\lambda_{1,0}=53.03$ | $\lambda_{1,1}=6.031$ | $\theta_{1}=0.685$ | $\lambda_{2,0}=12.03$ | $\lambda_{2,1}=8.325$ | $\theta_{2}=0.2876$
---|---|---|---|---|---

$\beta_{0}=0.0381$ | $\beta_{1}=-0.1483$ | $\beta_{2}=0.7097$ | $\beta_{1,2}=0.1671$
---|---|---|---

(a)

(b)

Figure 8: Comparison of the market and model term-structures of ATM skew.

The calibration to the SPX surface up to the one-year maturity takes
approximately 30 minutes using Cuda 11.4 with PyTorch on a GPU-NVIDIA-SMI. The
optimization was carried out, starting from a randomized offline initial
configuration, using the derivative-free optimizer Py-BOBYQA (see Cartis et
al. (2019, 2022)) relying on a trust-region method which proved particularly
efficient with respect to other optimizers. We enforced natural bounds on the
parameters during optimization and then scaled model parameters to help the
learning.

##  5 Joint calibration of SPX and VIX options

We now report the results of the joint calibration of SPX and VIX monthly
options as of June 2, 2021 (Figure 9 and Table 5) and June 3, 2021 (Figure 10
and Table 6). As loss function we employ the sum of the two loss functions
introduced in (3.5) and (4.3). Recall that SPX monthly options expire the
third Friday on the expiration month, while the monthly VIX futures and VIX
options mature 30 days before the third Friday of the following month.

Figure 9:  Joint calibration as of June 2, 2021. Comparison of calibrated and
market SPX/VIX implied volatility smiles. Top left and bottom: SPX smiles. Top
right: the market VIX future is displayed with dotted red lines while the
model VIX future is indicated by the dotted blue line; in black (resp., dotted
black) the VIX smile (resp., VIX future) computed with nested Monte Carlo. The
initial values of the factors given the calibrated parameters in Table 5 are
$R_{1,0}=0.2689$, $R_{1,1}=0.2375$, $R_{2,0}=0.0249$, $R_{2,1}=0.02491$. Table
5:  Calibrated parameters of the 4FPDV model. Joint calibration to the SPX and
VIX smiles as of June 2, 2021.

$\lambda_{1,0}=44.42$ | $\lambda_{1,1}=33.19$ | $\theta_{1}=0.398$ | $\lambda_{2,0}=4.311$ | $\lambda_{2,1}=3.254$ | $\theta_{2}=0.72$
---|---|---|---|---|---

$\beta_{0}=0.0254$ | $\beta_{1}=-0.1602$ | $\beta_{2}=0.6922$ | $\beta_{1,2}=0.1639$
---|---|---|---

Figure 10:  Joint calibration as of June 3, 2021. See Figure 9. The initial
values of the factors given the calibrated parameters in Table 6 are
$R_{1,0}=0.0669$, $R_{1,1}=0.0916$, $R_{2,0}=0.02197$, $R_{2,1}=0.02725$.
Table 6:  Calibrated parameters of the 4FPDV model. Joint calibration to the
SPX and VIX smiles as of June 3, 2021.

$\lambda_{1,0}=42.78$ | $\lambda_{1,1}=31.51$ | $\theta_{1}=0.389$ | $\lambda_{2,0}=3.694$ | $\lambda_{2,1}=3.693$ | $\theta_{2}=0.698$
---|---|---|---|---|---

$\beta_{0}=0.0264$ | $\beta_{1}=-0.1665$ | $\beta_{2}=0.6829$ | $\beta_{1,2}=0.1628$
---|---|---|---

In Figures 9 and 10 we observe that we are able to calibrate jointly the first
two monthly smiles of the SPX (top left and bottom) and the first monthly VIX
future and VIX smile (top right) with a very good accuracy. Note that for both
calibration dates, $\lambda_{1,0}\approx\lambda_{1,1}$ and
$\lambda_{2,0}\approx\lambda_{2,1}$, which means that one-exponential kernels
seem enough for this fit. This is due to the fact that we only calibrate to
short-dated options here. Recall that, by contrast, in Section 4.2, we have
seen that two-exponential kernels are needed to jointly fit short-dated and
long-dated SPX implied volatilities. For each iteration of the optimizer we
simulated $N_{MC}=2\cdot 10^{5}$ trajectories with discretization step $\Delta
t=\frac{1}{504}$. The hyperparameters of the loss functions for SPX and VIX
options are fixed to be
$(\omega_{{\textrm{SPX}}},\omega_{{\textrm{VIX}}},\omega_{F})=(10,5,20)$. Like
for the calibration to the SPX surface, we use the Py-BOBYQA optimizer and
randomize the initial guess within the natural bound of the parameters. The
joint calibration takes around 8 minutes to complete.

In Figures 9 and 10, we do not only report the VIX future and VIX implied
volatilities in the calibrated model when we use our neural network
approximation $\mathcal{N}\mathcal{N}^{\star}$ of the VIX. We also report
those quantities when we use nested Monte Carlo paths to estimate the VIX.
Figures 9 and 10 show that for the jointly calibrated parameters our neural
network approximation of the VIX is accurate enough for trading purposes.

##  6 On the stability of the calibrated parameters

In this section we address the stability of the calibrated parameters for the
4FPDV model. Stability of the calibrated parameters is desirable as, usually,
models are periodically recalibrated for hedging purposes; parameter stability
prevents oscillating hedge ratios and higher transaction costs. A similar
analysis has been carried out in Cuchiero et al. (2023a).

We split our analysis in two parts. In the first part (Section 6.1), we
calibrate OTM SPX monthly options every day of the fourth trading week of
October 2023. In the second part (Section 6.2), we compare the parameters
found for the two joint calibration exercises in Section 5, i.e., for the two
consecutive days June 2, 2021 and June 3, 2021. In both cases we display the
corresponding kernels of the 4FPDV model as defined in (2.3).

###  6.1 Calibration to SPX options only

MAE
---
October 23, 2023 $2.23\cdot 10^{-3}$
October 24, 2023 $2.26\cdot 10^{-3}$
October 25, 2023 $2.12\cdot 10^{-3}$
October 26, 2023 $2.63\cdot 10^{-3}$
October 27, 2023 $2.32\cdot 10^{-3}$
Table 7: Mean absolute error of the SPX implied volatility fit. The range of
moneyness considered for each maturity $T>0$, is given by
$[K_{\min},K_{\max}]$ where $K_{\min}=1-0.4\sqrt{T}$ and
$K_{\max}=1+0.25\sqrt{T}$.

In this section, we calibrate OTM SPX monthly options every day of the fourth
trading week of October 2023, that is, October 23 to October 27, 2023. We
calibrate to monthly maturities up to one year. Figure 7 shows the result of
the calibration as of October 25, 2023. Table 7, in which we report the mean
absolute error (MAE) of the implied volatility fit, shows that the fits as of
October 23, 24, 26, and 27, not reported here, achieve the same degree of
accuracy.

The values of the calibrated parameters on the five consecutive trading days
are reported in Figure 11. The $\beta$ parameters
$(\beta_{0},\beta_{1},\beta_{2},\beta_{1,2})$ (Figure 11(a)) are remarkably
stable over time. The kernel parameters
$(\lambda_{1,0},\lambda_{1,1},\theta_{1},\lambda_{2,0},\lambda_{2,1},\theta_{2})$
(Figure 11(b)) appear less stable, in particular the $K_{2}$ parameters.
However, note that different values of
$(\lambda_{1,0},\lambda_{1,1},\theta_{1})$ (resp.,
$(\lambda_{2,0},\lambda_{2,1},\theta_{2})$) produce similar two-exponential
kernels. This is illustrated in Figure 12, and is mostly observed for $K_{1}$.

(a) Sensitivities $(\beta_{0},\beta_{1},\beta_{2},\beta_{1,2})$ of the
volatility for five consecutive calibration dates.

(b) Parameters of the kernels $K_{1}$ and $K_{2}$ for five consecutive
calibration dates. L.h.s.: parameters
$(\lambda_{1,0},\lambda_{1,1},\theta_{1})$ of the kernel $K_{1}$. R.h.s.:
parameters $(\lambda_{2,0},\lambda_{2,1},\theta_{2})$ of the kernel $K_{2}$.
For each plot, the $y$-axis on the l.h.s. refers to the natural scale of the
$\lambda_{n,p}$s while the $y$-axis on the r.h.s. refers to the values of the
$\theta_{p}$s.

Figure 11:  Figure 12: Comparison of the kernels on five consecutive trading
days. L.h.s.: $K_{1}$. R.h.s.: $K_{2}$

###  6.2 Joint SPX/VIX calibration

In Figure 13, we plot the jointly calibrated parameters of Tables 5 (June 2,
2021) and 6 (June 3, 2021). We observe that the calibrated parameters are very
stable from one calibration date to the next, including the kernel parameters.
In particular Figure 14 illustrates that both kernels $K_{1}$ and $K_{2}$ do
not change drastically from one calibration date to the next.

(a) Sensitivities $(\beta_{0},\beta_{1},\beta_{2},\beta_{1,2})$ of the
volatility for two consecutive calibration dates.

(b) Parameters of the kernels $K_{1}$ and $K_{2}$ for two consecutive
calibration dates. L.h.s.: parameters
$(\lambda_{1,0},\lambda_{1,1},\theta_{1})$ of the kernel $K_{1}$. R.h.s.:
parameters $(\lambda_{2,0},\lambda_{2,1},\theta_{2})$ of the kernel $K_{2}$.
For each plot, the $y$-axis on the l.h.s. refers to the natural scale of the
$\lambda_{n,p}$s while the $y$-axis on the r.h.s. refers to the values of the
$\theta_{p}$s.

Figure 13:  Figure 14: Comparison of the calibrated exponential kernels
jointly to SPX and VIX. The two vertical dotted lines denote the smallest
maturity considered for calibration (13 days) and the largest one (44 days).

##  Appendix A Appendix: Additional tests of the neural network

In this appendix, we report the error between the neural network predictions
${\mathcal{N}}{\mathcal{N}}^{\star}$ and the nested Monte Carlo estimator
$\widehat{{\textrm{VIX}}}$ for market-calibrated parameters, see Figures 15
and 16.

(a)

(b)

Figure 15:  Model parameters are jointly calibrated on June 2, 2021. Left:
comparison between $N_{MC}=3\cdot 10^{5}$ realizations of the nested Monte
Carlo estimator $\widehat{{\textrm{VIX}}}$ and the neural approximation
${\mathcal{N}}{\mathcal{N}}^{\star}$. Here the $\widehat{{\textrm{VIX}}}$ is
computed with $\Delta t=\frac{1}{2520}$, $3\cdot 10^{5}$ outer paths and
$10^{4}$ nested paths. Right: histogram of the absolute error between the
neural network predictions ${\mathcal{N}}{\mathcal{N}}^{\star}$ and
$\widehat{{\textrm{VIX}}}$. The average absolute error is 0.202 $\%$.

(a)

(b)

Figure 16: Same as Figure 15, with model parameters jointly calibrated on June
3, 2021. The average absolute error between the neural network predictions
${\mathcal{N}}{\mathcal{N}}^{\star}$ and $\widehat{{\textrm{VIX}}}$ is 0.185
$\%$.

## References

  * Abi Jaber and Li (2024) E. Abi Jaber and S. Li.  Volatility models in practice: Rough, Path-dependent or Markovian?  _Preprint arXiv:2401.03345_ , 2024.
  * Abi Jaber et al. (2019) E. Abi Jaber, M. Larsson, and S. Pulido.  Affine Volterra processes.  _Annals of Applied Probability_ , 29(5):3155–3200, 2019.
  * Abi Jaber et al. (2022) E. Abi Jaber, C. Illand, and S. Li.  Joint SPX-VIX calibration with Gaussian polynomial volatility models: deep pricing with quantization hints.  _Preprint arXiv:2212.08297_ , 2022.
  * Alfonsi et al. (2023) A. Alfonsi, B. Lapeyre, and J. Lelong.  How many inner simulations to compute conditional expectations with least-square Monte Carlo?  _Methodology and Computing in Applied Probability_ , 25(3):71, 2023.
  * Andrès et al. (2023) H. Andrès, A. Boumezoued, and B. Jourdain.  Implied volatility (also) is path-dependent.  _Preprint arXiv:2312.15950_ , 2023.
  * Baldeaux and Badran (2014) J. Baldeaux and A. Badran.  Consistent modelling of VIX and equity derivatives using a 3/2 plus jumps model.  _Applied Mathematical Finance_ , 21(4):299–312, 2014.
  * Baschetti et al. (2024) F. Baschetti, G. Bormetti, and P. Rossi.  Deep calibration with random grids.  _Quantitative Finance_ , 1(1):1–23, 2024.
  * Bayer and Stemper (2018) C. Bayer and B. Stemper.  Deep calibration of rough stochastic volatility models.  _Preprint arXiv:1810.03399_ , 2018.
  * Blanc et al. (2017) P. Blanc, J. Donier, and J.-P. Bouchaud.  Quadratic Hawkes processes for financial prices.  _Quantitative Finance_ , 17(2):171–188, 2017\.
  * Bondi et al. (2024a) A. Bondi, G. Livieri, and S. Pulido.  Affine volterra processes with jumps.  _Stochastic Processes and their Applications_ , 168:104264, 2024a.
  * Bondi et al. (2024b) A. Bondi, S. Pulido, and S. Scotti.  The rough Hawkes Heston stochastic volatility model.  _Mathematical Finance_ , pages 1–45, 2024b.
  * Bourgey and Guyon (2024) F. Bourgey and J. Guyon.  Fast exact joint S&P 500/VIX smile calibration in discrete and continuous time.  _Risk_ , February, 2024.
  * Cartis et al. (2019) C. Cartis, J. Fiala, B. Marteau, and L. Roberts.  Improving the flexibility and robustness of model-based derivative-free optimization solvers.  _ACM Transactions on Mathematical Software (TOMS)_ , 45(3):1–41, 2019.
  * Cartis et al. (2022) C. Cartis, L. Roberts, and O. Sheridan-Methven.  Escaping local minima with local derivative-free methods: a numerical investigation.  _Optimization_ , 71(8):2343–2373, 2022.
  * Chicheportiche and Bouchaud (2014) R. Chicheportiche and J. Bouchaud.  The fine-structure of volatility feedback I: Multi-scale self-reflexivity.  _Physica A: Statistical Mechanics and its Applications_ , 410:174–195, 2014.
  * Cont and Ben Hamida (2005) R. Cont and S. Ben Hamida.  Recovering volatility from option prices by evolutionary optimization.  _Journal of Computational Finance_ , 8(4):43–76, 2005.
  * Cont and Kokholm (2013) R. Cont and T. Kokholm.  A consistent pricing model for index options and volatility derivatives.  _Mathematical Finance_ , 23(2):248–274, 2013\.
  * Cuchiero and Teichmann (2019) C. Cuchiero and J. Teichmann.  Markovian lifts of positive semidefinite affine Volterra type processes.  _Decisions in Economics and Finance_ , 42(2):407–448, 2019.
  * Cuchiero et al. (2023a) C. Cuchiero, G. Gazzani, J. Möller, and S. Svaluto-Ferro.  Joint calibration to SPX and VIX options with signature-based model.  _Preprint arXiv:2301.13235_ , 2023a.
  * Cuchiero et al. (2023b) C. Cuchiero, G. Gazzani, and S. Svaluto-Ferro.  Signature-based models: Theory and calibration.  _SIAM Journal on Financial Mathematics_ , 14(3):910–957, 2023b.
  * Di Nunno et al. (2023) G. Di Nunno, K. Kubilius, Y. Mishura, and A. Yurchenko-Tytarenko.  From constant to rough: A survey of continuous volatility modeling.  _Mathematics_ , 11(19):4201, 2023.
  * Dupire (1994) B. Dupire.  Pricing with a smile.  _Risk_ , 7(1):18–20, 1994.
  * Foschi and Pascucci (2005) P. Foschi and A. Pascucci.  Calibration of the Hobson and Rogers model: empirical tests.  _Computational Statiscs & Data Analysis_, 53:2219–2235, 2005.
  * Fouque and Saporito (2018) J.-P. Fouque and Y. Saporito.  Heston stochastic vol-of-vol model for joint calibration of VIX and S&P 500 options.  _Quantitative Finance_ , 18(6):1003–1016, 2018\.
  * Gambara and Teichmann (2020) M. Gambara and J. Teichmann.  Consistent recalibration models and deep calibration.  _Preprint arXiv:2006.09455_ , 2020.
  * Gatheral (2008) J. Gatheral.  Consistent modeling of SPX and VIX options.  _Bachelier Congress_ , 2008.
  * Gatheral et al. (2020) J. Gatheral, P. Jusselin, and M. Rosenbaum.  The quadratic rough Heston model and the joint S&P 500/VIX smile calibration problem.  _Risk_ , May, 2020.
  * Guo et al. (2022) I. Guo, G. Loeper, J. Obloj, and S. Wang.  Optimal transport for model calibration.  _Risk_ , January, 2022.
  * Guyon (2020a) J. Guyon.  Inversion of convex ordering in the VIX market.  _Quantitative Finance_ , 20(10):1597–1623, 2020a.
  * Guyon (2020b) J. Guyon.  The joint S&P 500/VIX smile calibration puzzle solved.  _Risk_ , April, 2020b.
  * Guyon (2022a) J. Guyon.  The VIX future in Bergomi models: Fast approximation formulas and joint calibration with S&P 500 skew.  _SIAM Journal on Financial Mathematics_ , 13(4):1418–1485, 2022a.
  * Guyon (2022b) J. Guyon.  Dispersion-constrained martingale Schrödinger bridges: Joint entropic calibration of stochastic volatility models to S&S 500 and VIX smiles.  _Available at SSRN 4165057_ , 2022b.
  * Guyon (2024) J. Guyon.  Dispersion-constrained martingale Schrödinger problems and the exact joint S&P 500/VIX smile calibration puzzle.  _Finance and Stochastics_ , 28(1):27–79, 2024\.
  * Guyon and Lekeufack (2023) J. Guyon and J. Lekeufack.  Volatility is (mostly) path-dependent.  _Quantitative Finance_ , 23(9):1221–1258, 2023\.
  * Guyon and Mustapha (2023) J. Guyon and S. Mustapha.  Neural joint S&P 500/VIX smile calibration.  _Risk_ , December, 2023.
  * Henry-Labordère (2019) P. Henry-Labordère.  From (martingale) Schrödinger bridges to a new class of stochastic volatility models.  _Preprint arXiv:1904.04554_ , 2019.
  * Hernandez (2016) A. Hernandez.  Model calibration with neural networks.  _Available at SSRN 2812140_ , 2016.
  * Hobson and Rogers (1998) D. Hobson and L. Rogers.  Complete models with stochastic volatility.  _Mathematical Finance_ , 8(1):27–48, 1998.
  * Horvath et al. (2021) B. Horvath, A. Muguruza, and M. Tomas.  Deep learning volatility: a deep neural network perspective on pricing and calibration in (rough) volatility models.  _Quantitative Finance_ , 21(1):11–27, 2021.
  * Kokholm and Stisen (2015) T. Kokholm and M. Stisen.  Joint pricing of VIX and SPX options with stochastic volatility and jump models.  _The Journal of Risk Finance_ , 16(1):27–48, 2015\.
  * Li et al. (2017) L. Li, K. Jamieson, G. De Salvo, A. Rostamizadeh, and A. Talwalkar.  Hyperband: A novel bandit-based approach to hyperparameter optimization.  _The journal of machine learning research_ , 18(1):6765–6816, 2017.
  * Neuberger (1994) A. Neuberger.  The log contract.  _Journal of portfolio management_ , 20:74–74, 1994.
  * Pacati et al. (2018) C. Pacati, G. Pompa, and R. Renò.  Smiling twice: the Heston++ model.  _Journal of Banking & Finance_, 96:185–206, 2018.
  * Papanicolaou and Sircar (2014) A. Papanicolaou and R. Sircar.  A regime-switching Heston model for VIX and S&P 500 implied volatilities.  _Quantitative Finance_ , 14(10):1811–1827, 2014\.
  * Parent (2023) L. Parent.  The EWMA Heston model.  _Quantitative Finance_ , 23(1):71–93, 2023.
  * Rømer (2022) S. E. Rømer.  Empirical analysis of rough and classical stochastic volatility models to the SPX and VIX markets.  _Quantitative Finance_ , 22(10):1805–1838, 2022\.
  * Rosenbaum and Zhang (2022) M. Rosenbaum and J. Zhang.  Deep calibration of the quadratic rough Heston model.  _Risk_ , September, 2022.
  * Sentana (1995) E. Sentana.  Quadratic ARCH models.  _The Review of Economic Studies_ , 62(4):639–661, 1995.
  * Zumbach (2009) G. Zumbach.  Time reversal invariance in finance.  _Quantitative Finance_ , 9(5):505–515, 2009.
  * Zumbach (2010) G. Zumbach.  Volatility conditional on price trends.  _Quantitative Finance_ , 10(4):431–442, 2010\.