Joint calibration to SPX and VIX options with signature-based models

Introduction

The joint calibration to SPXThe SPX is a theoretical index, in the sense that it is pegged to the value of the S&P 500 without being based on a held portfolio of stock shares. In the present paper we shall use SPX and S&P 500 interchangeably. and VIX options is a problem that has gained a lot of attention in quantitative finance since several years. It is sometimes still regarded as the holy grail of volatility modeling, even though significant progress has been made recently (see below in the literature overview).

One main reason for the increased interest is that the VIX index is no longer only used as indicator of volatility, but rather as important underlying for many derivatives. In fact, futures and options written on it are extensively used to hedge the volatility exposure of option portfolios, see e.g. Rhoads (2011). Suppose that an investor has a long position on the S&P 500 index. Although she might believe it has long-term prospects, it would be desirable to reduce her exposure to short-term volatility. Buying VIX derivatives with the belief that volatility is going to increase she might balance out these positions, while she was wrong, the losses to her VIX position could be mitigated by gains to the existing trade.

Historically speaking, in March 2004, futures on VIX, which trade on the forward 30-days realized volatility of the S&P 500, were launched on the Chicago Board Options Exchange (CBOE). Two years later, in February 2006, also trading in European options written on the VIX index started and has progressively increased since then. We address the reader to the official website of CBOEwww.cboe.com/tradable$\_$products/vix/ for details on how the VIX is computed and traded in practice, but just remark the following in view of expiration dates.

Note that for a given type of option just certain maturities are available. This in particular implies that on different days options with different times to maturity are traded. For example, options which expire exactly in one year are not available on a daily basis. Concerning our concrete applications, typically, VIX options expire on the Wednesday 30 days (or closest to 30 days) prior to the third Friday of the next calendar month. On the other hand a monthly SPX option typically expires on the third Friday of every month.

The steadily increasing liquidity on the market of VIX index products has marked the need to jointly calibrate to both, the prices of options on the underlying S&P 500 index and to prices of VIX derivatives.

From a data perspective the challenge, especially for short maturities, is to reconcile the large negative skew of SPX options’ implied volatilities with relatively lower implied volatilities arising from the VIX options. This has been discussed rigorously in Guyon (2020a), where the author shows a necessary condition to solve the joint calibration problem. As observed in the paper this translates, in the context of models with continuous trajectories, to a volatility process with high mean-reversion speed and large negative correlation to the S&P 500 index. In addition to that a high vol-of-vol is desirable in order to reproduce the aforementioned negative skew of at-the-money SPX options, but this can yield too high implied volatilities for the VIX options.

Inspired by Perez Arribas et al. (2020) and Cuchiero et al. (2022a), we consider here a new type of stochastic volatility model for the discounted price process $S=(S_{t})_{t\geq 0}$ with continuous trajectories. It is given by

Let us now highlight the implications of this modeling framework and the novelty of the present work.

The modeling framework can be seen as universal in the class of continuous non-rough stochastic volatility models, which is a consequence of the universal approximation result stated in Theorem 2.6.

It is not only universal in an approximate sense, but truly nests several classical models (see Remark 3.4) and for instance also the ‘quintic Ornstein-Uhlenbeck volatility model’, recently proposed by Abi Jaber et al. (2022b), which – with an additional input curve – is shown to fit SPX and VIX smiles well.

Up to our knowledge, it is the first signature-based model that is employed for pricing and calibration of VIX options as well as joint calibration, together with SPX options.

We illustrate that the joint calibration problem can be solved in this framework without jumps and rough volatility (compare also Guyon and Lekeufack (2022); Rømer (2022); Abi Jaber et al. (2022b)).

By using time-varying parameters we can go beyond short maturities both for SPX and VIX options (as classically tackled in the literature) and achieve a joint calibration also for longer maturities.

In order to achieve the highly accurate calibration results, illustrated in Section 5.4 and Section 7, we exploit the following mathematical and numerical properties.

Alternatively, a Fourier pricing approach for both VIX and SPX options can be used. Indeed, by building on the fact that the signature of $\widehat{Z}$ is an affine process (with values in the extended tensor algebra) as proved in Cuchiero et al. (2023), its Fourier-Laplace transform can be computed by solving an (extended tensor algebra valued) Riccati equation, which in turn can be used for Fourier-pricing as outlined in Section 6.1.

The remainder of the paper is organized as follows. Section 1.1 gives a review over the different contributions in the literature concerning the joint calibration problem. In Section 2 we introduce the signature in the context of continuous semimartingales, its main properties as well as notation used throughout the paper. Section 3 is dedicated to the introduction of our signature-based model and the connections to classical and also recent stochastic volatility models in the literature. Section 4 is then devoted to the discussion and proof of the matrix exponential formula for the (conditional) expected signature of polynomial processes. This result is at the core of Section 5, where we derive a tractable formula for the VIX, needed for pricing VIX options and VIX futures. Building on these formulas, our calibration results to VIX options are presented in Section 5.4.1. In Section 6, we then prove, similarly as for the VIX, a tractable expression for $S$. Additionally, we exploit in Section 6.1 the affine nature of the signature process (as proved in Cuchiero et al. (2023)), to obtain a Fourier pricing approach within our modeling choice for both VIX and SPX options. We finally present the numerical results of the joint calibration problem in Section 7, both in the case of constant parameters and with time-varying parameters, where the latter are introduced in Section 5.5 and Section 6.2.

The data used in Section 5.4.1 and Section 7.1 were purchased from OptionMetricshttps://optionmetrics.com/. Python codes used to produce the numerical results of the present work are available from the authors upon request.

This section is primarily dedicated to a literature review on the joint calibration problem and secondly, to a brief overview on signature methods in finance.

First attempts to solve the joint calibration problem appear in Gatheral (2008), with a double constant elasticity of variance model (CEV), which despite being rather flexible cannot fit accurately the implied volatilities of SPX and VIX options jointly. Later on, the belief that jumps in the SPX (or additionally also in the volatility) are necessary has lead to the following contributions, Sepp (2012); Papanicolaou and Sircar (2014); Baldeaux and Badran (2014); Kokholm and Stisen (2015); Pacati et al. (2018) and more recently Grzelak (2022) with a new perspective on randomization.

Continuous stochastic volatility models based on Markovian semimartingales have also been employed to solve the joint calibration problem. For instance, in Fouque and Saporito (2018) a Heston model with stochastic vol-of-vol has been calibrated, however only for maturities above 4 months where 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 non-standard transformation function. The well-working choice of two OU-processes illustrated there has been an inspiration for our concrete numerical implementations. We also point out that the (non-rough) model introduced in Abi Jaber et al. (2022a, b), where the volatility is described by a polynomial of order five in one single OU-process, falls (apart from the additional input curve) into this class of continuous Markovian models and is a particular instance of our framework. Let us also refer to the paper by Guyon and Mustapha (2022), where a neural SDE model has been successfully jointly calibrated. Within the class of continuous, however not necessarily Markovian models, Guyon and Lekeufack (2022) conduct an empirical and statistical analysis as well as a joint calibration for a family of models where the volatility depends on the paths of the asset. These models can be turned into Markovian ones by using exponential kernels instead of general ones.

Two further distinct and rather new lines of research are worth being mentioned as well: first, martingale optimal transport and second rough volatility.

The martingale optimal transport approach (see e.g. Guyon and Henry-Labordere (2013)) is used to calibrate discrete-time models as proposed in Guyon (2020b, 2021). These models are closely related to Schrödinger bridge problems, where the idea is to calibrate only the drift of the volatility while keeping the volatility of volatility unchanged, see e.g. Henry-Labordere (2019) or Guo et al. (2020, 2021) as well as the references therein regarding an optimal transport approach. Although the calibration within that setting is accurate, it is also computationally rather expensive and not amenable to calibrate to several maturities jointly. These computational challenges have been tackled recently in Guyon and Bourgey (2022) both in discrete and continuous time extending the contributions of Guyon (2020b, 2021).

In the area of rough volatility modeling, initiated by the seminal paper of Gatheral et al. (2018), the main idea is to replace the standard Brownian motion in the volatility process by a fractional Brownian motion. Even though the roughness of the trajectories found in Gatheral et al. (2018), can also be explained by the estimation procedure as discussed e.g. in Cont and Das (2023), the non-Markovianity given by the fractional Brownian motion with Hurst parameter $H<0.5$, manages to reproduce many stylized facts arising in financial data. Several classical models have been enhanced with rougher noise, but for simplicity we mention those employed in the SPX/VIX calibration. One example is the quadratic rough Heston model introduced in Gatheral et al. (2020), which was in turn calibrated in Rosenbaum and Zhang (2021) by relying on neural networks approaches, also exploited in e.g. Bayer et al. (2019); Cuchiero et al. (2020). For a large class of rough volatility models Jacquier et al. (2021) give new insights on the joint calibration problem by providing small-time formulas of the at-the-money implied volatility, skew and curvature for both European options on SPX and VIX options. In Rømer (2022) an exhaustive study of the flexibility of different rough volatility models to joint SPX/VIX calibration is carried out, including the rough Bergomi, the rough Heston and an extended rough Bergomi model. Some of these, for instance the rough Heston model, have an affine structure i.e., can be embedded in the class of affine Volterra processes, considered in Abi Jaber et al. (2019); Cuchiero and Teichmann (2019, 2020). In particular they allow for Fourier pricing after solving the related Riccati equations. This underlying structure is the building block of an extension with jumps investigated in Bondi et al. (2022a) and recently employed in the context of the joint calibration in Bondi et al. (2022b), where a rough Heston model with Hawkes-type jumps with Fourier pricing is employed.

Concerning our framework, signature-based methods provide a generic non-parametric way to extract characteristic features (linearly) and path-dependency from data, which is essential in (machine) learning and calibration tasks in finance. This explains why these techniques become more and more popular in mathematical finance, see e.g., Buehler et al. (2020); Kalsi et al. (2020); Perez Arribas et al. (2020); Lyons et al. (2020); Ni et al. (2021); Bayer et al. (2021); Min and Hu (2021); Akyildirim et al. (2022); Cuchiero et al. (2022a); Cuchiero and Möller (2023); Cuchiero et al. (2022c) and the references therein.

Signature: definition and properties

For a multi-index $I:=(i_{1},\ldots,i_{n})$ we set $|I|:=n$. We also consider the empty index $I:=\emptyset$ and set $|I|:=0$. If $n\geq 1$ or $n\geq 2$ we set $I^{\prime}:=(i_{1},\ldots,i_{n-1})$, and $I^{\prime\prime}:=(i_{1},\ldots,i_{n-2})$, respectively. We also use the notation

omitting the parameter $d$ whenever this does not introduce ambiguity. Observe that multi-indices can be identified with words, as it is done for instance in Lyons et al. (2020).

Observe in particular that ${\mathbf{b}}_{I}=\langle e_{I},\textbf{b}\rangle$.

A well-known and extremely useful property of the signature is that every polynomial function in the signature has a linear representation. For the precise statement we first need to introduce the following concept (see also Definition 2.4 in Lyons et al. (2020) or Section 2.2. in Bayer et al. (2021)).

For every two multi-indices $I$ and $J$ the shuffle product is defined recursively as

The proof of the paths version of the next result for can be found for instance in Ree (1958) or Lyons et al. (2007).

The result follows by induction using the chain rule for Stratonovich integrals. ∎

We recall now an important property of the signature. The result is known in the rough paths literature (see for instance Boedihardjo et al. (2016)), but can also be proved directly in the simpler situation of a continuous semimartingale that contains time as strictly monotone component.

In order to combine the value of the signature on different time intervals Chen’s identity going back to Chen (1957, 1977) turns out to be fundamental.

for each $0\leq s\leq u\leq t$. This can equivalently be written as

See Cuchiero et al. (2022a) for a direct proof using the definition of Stratonovich integrals. ∎

Let us recall also the universal approximation theorem of linear functions of the signature in the context of continuous semimartingales as stated in Cuchiero et al. (2022a). We refer to Cuchiero and Möller (2023), in the same context, for the situation involving not just the approximation of a functional up to a fixed final value $T$, but a uniform approximation on the whole time interval $[0,T]$.

is continuous for each multi-index $I$ and every $t\in[0,T]$.

See Theorem 2.12 and Remark 2.13 in Cuchiero et al. (2022a) for details concerning the proof and the choice of the metric $d_{{\mathcal{S}}^{(2)}}$. ∎

For the related concept of stochastic Taylor expansions and functional expansions we refer the reader to Section 2.3 in Cuchiero et al. (2022a) and to Dupire and Tissot-Daguette (2022), respectively.

Finally, we introduce the concept of polynomial diffusion process which will play a key role for the computation of conditional expected signatures. Here we denote by $\sqrt{{\,\cdot\,}}$ the matrix square root.

The model

As an alternative definition for the volatility process $(\sigma_{t}^{S})_{t\geq 0}$ one can set

for some fixed ${\varepsilon}>0$. In this case the value of the volatility process $\sigma^{S}$ at time $t$ does not depend on the whole trajectory of the primary process $X$, but just on its evolution from $t-{\varepsilon}$ to $t$. For an economically reasonable choice for ${\varepsilon}$ the lags used in Section 3.4 of Guyon and Lekeufack (2022) can be adapted to the current setting.

In the model given by (3.1) we describe the discounted prices and construct the VIX from them, in line with the definition of the CBOE for the computation of the VIX. However, contingent claims are often expressed in terms of undiscounted prices. If the dynamics of the discounted price process are given by (3.1), the undiscounted one fulfills

Possible choices for the primary process $X$ can be tractable processes such as an Ornstein-Uhlenbeck (OU) process or a Brownian motion. For the proposed analysis we indeed need that the corresponding conditional expected signature can be computed easily. Consider the following assumption.

$\widehat{X}=(t,X_{t})_{t\geq 0}$ is a polynomial diffusion process in the sense of Definition 2.7.

We illustrate here that several classical and also recently considered stochastic volatility models are nested within our modeling choice (3.2) under Assumption 3.3.

Expected signature of polynomial diffusion processes

Let $(Y_{t})_{t\geq 0}$ be a polynomial diffusion process in sense of Definition 2.7 whose dynamics are given by

We now explain how to employ the polynomial technology to compute the conditional expected signature of $(Y_{t})_{t\geq 0}$. Several representations of related quantities in particular for the Brownian case can be found in the literature, see for instance Fawcett (2003), Lyons and Victoir (2004), Lyons and Ni (2015), Boedihardjo et al. (2021), Rossi Ferrucci and Cass (2022). Our approach follows Cuchiero et al. (2023) and is based on the classical theory of polynomial processes (see Cuchiero et al. (2012) and Filipović and Larsson (2016)). Even though results for the corresponding infinite dimensional stochastic processes (see for instance Cuchiero and Svaluto-Ferro (2021); Cuchiero et al. (2021b)) are needed in the case of general signature SDEs considered in Cuchiero et al. (2023), the polynomial property of $(Y_{t})_{t\geq 0}$ here permits to stay in the finite dimensional setting.

Let $(Y_{t})_{t\geq 0}$ be the polynomial process given by (4.1) and $b$ and $a$ be the corresponding drift and diffusion coefficients. Then

Observe that the upper index on $Y_{0}^{k}$ and $Y_{0}^{h}$ refers to $Y$’s components and not to powers.

Let $\sigma_{j}(Y_{t})$ denote the $j$-th row of $\sigma(Y_{t})$. By definition of the signature, Stratonovich integral and by the shuffle property we can compute

the $d_{n}$-dimensional matrix representative of $L$.

where $e^{({\,\cdot\,})}$ denotes the matrix exponential.

For the present paper a crucial role is played by the polynomial process given by time, a $d$-dimensional OU process, and a Brownian motion. Specifically, we consider the process $\widehat{Z}_{t}:=(\widehat{X}_{t},B_{t})$ where $B$ is a Brownian motion and $\widehat{X}_{t}=(t,X_{t})$ with

for $a_{ij}(X_{t})=\sigma^{i}\sigma^{j}\rho_{ij}$, and $W$ being a $d$-dimensional Brownian motion. We denote by $\rho_{j(d+1)}$ the correlation between $X^{j}$ and $B$. Setting $\kappa^{d+1}:=0$ and $\sigma^{d+1}:=1$ we can see that $\widehat{Z}$ satisfies (4.1) in $d+2$ dimensions for

The corresponding ${\mathbf{b}}$ and ${\mathbf{a}}$ are given by ${\mathbf{b}}_{j}=e_{\emptyset}(1_{\{j=0\}}+\kappa^{j}(\theta^{j}-\widehat{Z}_{0}^{j})1_{\{j\neq 0\}})-e_{j}\kappa^{j}1_{\{j\neq 0\}}$ and ${\mathbf{a}}_{ij}=e_{\emptyset}\sigma^{i}\sigma^{j}\rho_{ij}1_{\{i,j\neq 0\}}$ and we thus get

An application of $L$ to the first basis elements yields the following results:

$L(e_{1})=e_{\emptyset}\kappa^{1}(\theta^{1}-X_{0}^{1})-e_{1}\kappa^{1}$;

$L(e_{I}\otimes e_{0})=e_{I}\shuffle{\mathbf{b}}_{0}+e_{I^{\prime}}\shuffle{\mathbf{a}}_{i_{|I|}0}=e_{I}$;

$L(e_{0}\otimes e_{1}\otimes e_{2})=e_{0}\otimes e_{1}\kappa^{2}(\theta^{2}-X_{0}^{2})-(e_{0}\otimes e_{1})\shuffle e_{2}\kappa^{2}+\frac{1}{2}e_{0}\sigma^{1}\sigma^{2}\rho_{12}$.

Letting $({\mathcal{F}}_{t})_{t\geq 0}$ be the filtration generated by $(\widehat{Z}_{t})_{t\geq 0}$ by Theorem 4.4 we can conclude that

Observe that given a subset $E\subseteq\{0,\ldots,d+1\}$, setting ${\mathcal{I}}_{E}:=\{I\colon i_{j}\in E\}$ it holds $L({\mathcal{I}}_{E})\subseteq{\mathcal{I}}_{E}$. This in particular implies that

for each $I\in{\mathcal{I}}_{E}$. Choosing $E=\{0,\ldots,d\}$, letting ${\mathscr{L}}_{E}:{\mathcal{I}}_{E}\to\{1,\ldots,(d+1)_{n}\}$ be a labelling injective function, and setting $G_{{\mathscr{L}}_{E}(I){\mathscr{L}}_{E}(J)}^{E}:=\eta_{IJ}$ we can see that (4.5) reduces to

To simplify the notation we often drop the $E$ from $G^{E}$ whenever this does not introduce any confusion.

VIX options with signatures

In this section we discuss the implication on pricing VIX options under the model (3.1)-(3.2) when Assumption 3.3 is in force. The implications of these on the log-price will be investigated in Section 6.

The CBOE Volatility Index, known by its ticker symbol VIX, is a popular measure of the market’s expected volatility of the S$\&$P 500 index, calculated and published by the Chicago Board Options Exchange (CBOE). The current VIX index value quotes the expected annualized change in the S&P 500 index over the following 30 days, based on options-based theory and current options-market data, more precisely

where $\Delta=30$ days and $S_{T}$ denotes the price process at time $T>0$. With the term VIX options we here usually refer to either put or calls written on VIX. In the present work we will take into account without loss of generality only call options.

Let $S=(S_{t})_{t\geq 0}$ be a price process described by

where $\sigma^{S}=(\sigma_{t}^{S})_{t\geq 0}$ denotes the volatility process, $B=(B_{t})_{t\geq 0}$ a one-dimensional Brownian motion.

Assume that $V=(\sigma^{S})^{2}$ satisfies

Part (i) follows directly from an application of Itô’s formula. Indeed, for any $t,T\geq 0$

where the integral with respect to the Brownian motion $B=(B_{t})_{t\geq 0}$ vanishes once we take the risk neutral conditional expectation, due to (5.2). This yields the alternative expression for ${\textrm{VIX}}_{T}^{2}$ and thus for ${\textrm{VIX}}_{T}$.

where for each $T>0$ the matrix $Q$ is given by

By Theorem 4.4 we can rewrite the matrix $Q$ as

Consider now the model described in Remark 3.1 and set for simplicity ${\varepsilon}\geq\Delta$. Then the results of Theorem 5.1(ii) still hold however with

Note that since the integration’s variable $t$ appears twice in this expression the time integral cannot be incorporated in the signature.

Observe that accounting for the scaling factor of 100, conventionally introduced by CBOE, the VIX index squared can equivalently be redefined (see e.g., Rosenbaum and Zhang (2021); Rømer (2022)) as

where $T,t>0$ and $\Delta=\frac{1}{12}$, i.e., approximately 30 days. Notice that since the expressions (5.3) and (5.6) differ only by a scaling factor, all the theoretical results of the present work hold true disregarding this scaling. For sake of simplicity we will always use (5.3). We address the reader to Chapter 11 in Gatheral (2011) for further details about the conventions of CBOE and its link with (5.1).

Approximation of the time integral: e.g., via the trapezoidal rule also applied for ${\textrm{VIX}}^{2}$ in Bourgey and De Marco (2021). Hence if we consider the shuffled coordinates ${\textup{\bf vec}}(e_{I}\shuffle e_{J})$ of the exponential matrix we can use the symmetry of the shuffle to reduce the number of integrals to be solved from $((d+1)_{2n})^{2}$ to $\frac{(d+1)_{n}((d+1)_{n}+1)}{2}\cdot(d+1)_{2n}$, instead of $(d_{2n})^{2}$. Observe that for our integral the error of such an approximation is given by

Approximation of the matrix exponential: we can avoid to approximate the integral by approximating the matrix exponential. Assuming that

this can for instance be done via its Taylor expansion:

Observe that (5.8) holds true whenever the spectral radius, i.e., the maximal eigenvalue in absolute value, of the matrix $G^{\top}\Delta$ is less than 1 (see for instance Theorem 1.5 in Quarteroni et al. (2010)). This requirement suggests that for numerical purposes the parameters of the primary underlying process have to be chosen accordingly.

An interesting example is given by the case where $X$ is a $d$-dimensional correlated Brownian motion, as considered for instance in Cuchiero et al. (2022a). In this case the process has no linear drift and the corresponding matrix $G$ is nilpotent, meaning that $G^{n}=0,$ for each $n$ big enough.

In general, this Taylor approach permits to avoid a numerical integration and produces an accurate approximation, allocating as few memory as possible.

where here $\lVert{\,\cdot\,}\lVert$ denotes the Euclidean norm. We stress the fact that the Cholesky decomposition can be carried out offline, and the computational benefit is immediate if several samples of the signature are considered.

In the following remark we discuss a possible dimension reduction technique from which one can benefit computationally. We follow the approach of Cuchiero et al. (2022b, 2021a), where by employing the Johnson-Lindenstrauss Lemma a random projection of the signature is considered, to which we refer as a randomized signature. A first way to use this tool is the following.

2 Options on VIX

We here briefly describe some important aspect of options on VIX. First of all, note that VIX options are written on VIX futures. The price process of a VIX future contract with maturity $T>0$, is given by

The price of a VIX option depends on the maturity of the corresponding VIX future.

When calibrating to VIX options, we stress that we additionally calibrate to VIX futures’ prices, see Section 5.4. This is important since future prices under the calibrated model are employed to compute its implied volatility surface. Including VIX futures in the calibration leads to a consistent model, both for VIX options and VIX futures, see e.g. Pacati et al. (2018); Guo et al. (2020); Guyon (2020a, 2021). Using market prices of the VIX futures to invert the implied volatility surface could lead to inconsistencies if one would like to price further derivatives with the calibrated model.

In this respect, let us here also comment on the computation of implied volatilities for VIX call options. Two equivalent approaches can be used to compute implied volatilities for a given maturity $T>0$:

use $F_{0}(T)$ in the so-called Black formula for options on futures (see Section 2.2 of Papanicolaou and Sircar (2014)).

Consider $e^{-rT}F_{0}(T)$, with $r>0$ the interest rate, in the Black-Scholes formula.

Recall that the Black formula coincides with Black-Scholes’ one when choosing the underlying in the Black-Scholes formula to be $e^{-r(T-t)}F_{t}(T)$. This implies that the two approaches are equivalent, i.e. lead to the same implied volatilities. In Section 5.4 we will only consider futures’ prices at time $t=0$ and maturity time $T>0$, hence for sake of simplicity we use the notation $F(T)$ instead of $F_{0}(T)$.

3 Variance reduction for pricing VIX options

We here discuss variance reduction techniques (see e.g. Glasserman (2004)) that can speed up the calibration in the subsequently applied Monte Carlo approach further. The key idea is to introduce a control variate, namely an easy to evaluate random variable $\Phi^{cv}$ such that given $T>0$ and $K>0$,

A well-working example of control variates used for pricing and calibrating neural SDE models can be found in Cuchiero et al. (2020); Gierjatowicz et al. (2020), where $\Phi^{cv}$ is constructed from hedging strategies.

In the following we describe two possible choices of control variates, which consist of polynomials on VIX futures. We stress the fact that these can be seen as linear functions of the signature of the primary process $\widehat{X}$, hence they belong to the class of sig-payoffs, see Lyons et al. (2020); Perez Arribas et al. (2020) and Section 4.2.2 in Cuchiero et al. (2022a).

The first example is to employ the VIX squared as main ingredient, see for instance Bourgey and De Marco (2021); Guerreiro and Guerra (2022) for a similar choice within a rough Bergomi model for pricing VIX options. This is particularly easy to treat in our set up, as for any given maturity $T>0$ we have

where the constant $c_{T,K}$ maximizing the variance reduction is given by:

Notice that also in this case both $Q$ and $Q^{cv}$ satisfy the condition for applying the Cholesky decomposition, leading to a faster evaluation of the control variate as discussed in Remark 5.5. Note that the Cholesky decomposition cannot be applied to $Q-Q^{cv}$, as this is in general an indefinite matrix.

As a second example we consider a generic polynomial in ${\textrm{VIX}}^{2}$ as control variate by defining

see for instance Section 4.1.1 in Glasserman (2004).

We stress the fact that for $m=1$ the two control variates introduced coincide.

4 Calibration to VIX options

In this section we focus on the calibration to VIX options only. Let $\mathcal{T}$ be a set of maturities and $\mathcal{K}$ a collection of strikes. Consider the model given by (3.1) and (3.2) and assume that Assumption 3.3 is in force.

Using Monte Carlo compute an approximation of option and futures’ prices with $N_{MC}>0$ samples, i.e.

Observe that an auxiliary randomization can be employed in every optimisation step as discussed in Remark 5.6. Moreover, if we want to use control variates to reduce the variance of the Monte Carlo estimator as described in the previous section, we would consider

Due to the variance reduction the number of samples needed is $N_{VR}\ll N_{MC}$ and $\Phi^{cv}$ is as in Section 5.3

The calibration to VIX call options and the corresponding futures on $\mathcal{T}$ and $\mathcal{K}$ consists in minimizing the functional

where $\mathcal{L}$ denotes a real-valued loss function, $F_{{\textrm{VIX}}}^{mkt}(T)$ the market’s futures’ prices and

the market’s option bid/ask prices $\pi_{\text{{{VIX}}}}^{mkt,b}(T,K),\pi_{\text{{{VIX}}}}^{mkt,a}(T,K)$, and bid/ask implied volatilities $\sigma_{\text{{{VIX}}}}^{mkt,b}(T,K),\sigma_{\text{{{VIX}}}}^{mkt,a}(T,K)$, respectively. We will specify the choice of the function $\mathcal{L}$ in Section 5.4.1 and Section 7.1.

Select $J^{\ast}\in(J_{i})_{i=1}^{m}$ such that

In the present section we report the results of the calibration to VIX options only. Here we consider call options written on the VIX on the trading day 02/06/2021, the same as in Guyon and Lekeufack (2022). We stress that for such recent dates the bid-ask spreads for VIX options are rather tight with respect to older dated options as considered for instance in Gatheral et al. (2020); Bondi et al. (2022b). The maturities are reported in the following table with the corresponding range of strikes (in percentage) with respect to the market’s futures prices.

We underline that the shortest maturity considered is 14 days. Regarding our modeling choice we fix $d=2$, $n=3$, which means to calibrate $40$ parameters. For $X$ we choose a $2$-dimensional Ornstein-Uhlenbeck processes, see Example 4.5, with the following (hyperparameter) configuration:

where the last column of $\rho$ are the correlations with the Brownian motion $B$ driving the price process $S$. The motivation of this parameters choice is to mimic a rough or strong mean-reverting model as suggested in Rogers (2023); Rømer (2022). We refer to Appendix A for numerical results where we use only a correlated $2$-dimensional Brownian motion as primary process, which yields significantly worse results. Before stating the loss function $\mathcal{L}$ that we employed in the calibration task, let us make the following remark.

where we recognize for the derivatives with respect to $\sigma$ and $\xi$, the Greeks Vega and Delta, respectively.

Motivated by Remark 5.8 we propose, for a fixed maturity and strike price, the following loss-function for $\beta\in\{0,1\}$

$\upsilon^{mkt}$ and $\delta^{mkt}$ denote the Vega and Delta of the option under the Black-Scholes model which depend on the maturity and on the strike price;

$F$ and $F^{mkt}$ denote futures with maturity $T$ such that the variables $\xi,\xi^{mkt}$ appearing in Remark 5.8 are $\xi=e^{-rT}F$ and $\xi^{mkt}=e^{-rT}F^{mkt}$, respectively;

We observe that by Remark 5.8 minimizing $\mathcal{L}^{0}$ is equivalent to minimizing an upper bound of the square of the right-hand side of (5.13) normalized by the bid-ask spread of the implied volatilities. Note that we slightly abused notation, since $\upsilon^{mkt}$ and $\delta^{mkt}$ of course depend on the strike and the maturity.

If our aim does not consist in calibrating to the mid-price or mid-implied-volatility precisely, but we merely want to be within the bid-ask spreads we can set $\beta=1$

For the next calibration result we minimize ${\mathcal{L}}^{1}$ as introduced above with $N_{MC}=80000$ Monte Carlo samples for the previous maturities and strikes.

We observe that the calibrated VIX smiles fall systematically in the bid-ask corridor for all the maturities considered. We report additionally in the next tables the relative absolute error between the market future prices and the calibrated ones for each maturity, i.e.,

5 The case of time-varying parameters

Therefore the variance process reads as follows,

Assume that for a set of maturities $\mathcal{T}^{{\textrm{VIX}}}$ it holds that $|T_{i}-T_{j}|\geq\Delta$ for all $i\neq j$.

Let $\mathcal{T}^{{\textrm{VIX}}}$ be a set of maturities on the VIX index and let $Q(T,\tau)$ be the matrix as defined in (5.5) (here for general $\tau>0$ instead of $\Delta$). Then, under (5.17) the VIX squared at time $T_{i}\in\mathcal{T}^{{\textrm{VIX}}}$ is given by

Note that, if $T_{i+1}-T_{i}>\Delta$ (which is in particular holds under Assumption 5.10) then,

By the definition of the VIX, it holds that

and hence the first statement follows by the definition of $Q$ in (5.5). ∎

Notice that also in the case of Proposition 5.11, Remark 5.5 applies.

SPX as a signature-based model

The goal of this section is to express the discounted price of the SPX, modeled via (3.1)-(3.2)

Recall that by (3.2) $\sigma^{S}$ is parametrized as follows

where $\widehat{X}_{t}=(t,X_{t})$ with $X$ a $d$-dimensional continuous semimartingale. Before addressing a more tractable expression for $S$, that allows to avoid (Euler) simulation schemes, we recall the following well-known integrability result.

In the following we suppose without loss of generality that $S_{0}=1$.

where for $\mathscr{L}:\{I:|I|\leq n\}\to\{1,\dots,(d+1)_{n}\},$

which is not finite in general, e.g. if $n=2$, $c=\sqrt{2}(n\text{!})$ then by Jensen’s inequality it follows

The key idea is to rewrite (3.1) as a type of signature-based model in sense of Cuchiero et al. (2022a) including $B=(B_{t})_{t\geq 0}$ as part of the primary process. This is possible since Itô integrals with respect to primary process’ components can be rewritten as linear functions of the signature of the primary process itself. To do so, we introduce Assumption 6.4 on $Z=(X,B)$ in order to describe the correlation structure between $B$ and $X$.

For all $i,j\in\{1,\ldots,d+1\}$ it holds

Let $S=(S_{t})_{t\geq 0}$ satisfy (3.1) with $S_{0}=1$, and $\sigma^{S}=(\sigma_{t}^{S})_{t\geq 0}$ satisfy (3.2). Suppose additionally that $Z=(X,B)$ satisfies Assumption 6.4. Then,

for a labeling function $\mathscr{L}:\{I:|I|\leq n\}\to\{1,\dots,(d+1)_{n}\}$.

Consider again the model described in Remark 3.1. Then the results of Proposition 6.5 still hold with

where $U^{0}_{t}$ is the upper-triangular matrix of the Cholesky decomposition of $Q^{0}(t)$.

Observe that a possible control variate for reducing the variance of the Monte Carlo estimator for pricing SPX options is the value at maturity of the log-price process. This means,

where ${\mathbf{b}}_{j}=\kappa^{j}\theta^{j}e_{\emptyset}-\kappa_{j}e_{j}$ and ${\mathbf{a}}_{jj}=(\sigma^{j})^{2}e_{\emptyset}$, using that (with a small abuse of notation) $\kappa^{0}\theta^{0}:=1$, $\kappa^{j}:=0$ and $\sigma^{0}:=0$. Using the notation of (2.2) consider then the Riccati operator ${\mathcal{R}}$ given by

where ${\bm{\psi}}$ is a solutionWe refer to Cuchiero et al. (2023) for the appropriate solution concept. of the extended tensor algebra valued Riccati equation

The representation of the Fourier-Laplace transform described above can then be used for Fourier pricing. We dedicate the remaining part of this section to illustrate how this can be done.

From Fourier analysis we know that for $K>0$ and $C<0$ it holds

Let us now consider the case of VIX options where Fourier pricing can be applied by computing the Fourier-Laplace transform of VIX squared, see also Sepp (2008); Papanicolaou and Sircar (2014); Cao et al. (2020); Bondi et al. (2022b) and references therein for a Fourier-based approach to pricing VIX options. Fix a labelling injective function ${\mathscr{L}}:\{I\colon|I|\leq n\}\to\{1,\ldots,(d+1)_{(2n+1)}\}$ as introduced before (2.1) and recall that by Theorem 5.1(ii) it holds

for each $y\geq 0$. This in particular implies that

Analogous formulations in terms of the error function are also possible, see for instance Cao et al. (2020); Bondi et al. (2022b). In the same spirit one can also obtain a representation of future prices. We here do not provide an implementation of this Fourier pricing approach but numerical experiments can be found in Cuchiero et al. (2023).

2 The case of time-varying parameters

Let $S=(S_{t})_{t\geq 0}$ satisfy (3.1) with $S_{0}=1$, and $(\sigma_{t}^{S})_{t\geq 0}$ satisfy (5.16) for a set of maturities ${\mathcal{T}}^{\textrm{VIX}}=\{T_{1},\dots,T_{N}\}$. Recall that in this case $V=(V_{t})_{t\geq 0}$ satisfy

Let $Z_{t}=(X_{t},B_{t})$ for all $t\geq 0$, where $X=(X_{t})_{t\geq 0}$ is a $d$-dimensional continuous semimartingale. Then, under Assumptions 6.4 we get the following recursion for the log-price process

for a labeling function $\mathscr{L}:\{I:|I|\leq n\}\to\{1,\dots,(d+1)_{2n+1}\}$.

and we will calculate each integral separately. We start with the first one.

Using similar arguments and Lemma 3.10 in Cuchiero et al. (2022a), the second integral yields

Joint calibration of SPX and VIX options

We here consider again the model introduced in (3.1)-(3.2) under Assumption 3.3. Note that we just work with call options, but the setup can easily be extended also to other liquid options on the market. Again we denote by ${\mathcal{T}}^{{\textrm{VIX}}}$ and ${\mathcal{T}}^{{\textrm{SPX}}}$ the maturities set for options written on VIX and SPX, respectively. Similarly we use the notation ${\mathcal{K}}^{{\textrm{VIX}}}$ and ${\mathcal{K}}^{{\textrm{SPX}}}$ for the corresponding strikes. The functional to be minimized in order to achieve a joint calibration of the SPX/VIX options reads as follows:

for a real-valued function $\mathcal{L}$. Observe that with a slight abuse of notation we denote this function as the one for $L_{{\textrm{VIX}}}$, but for SPX options we do not have to calibrate to futures, hence the last term of (5.14) vanishes.

By Proposition 6.5 the SPX call option payoff with maturity $T>0$ and a strike price $K>0$ reads in our model as follows

where $U_{T}$ denotes the upper-triangular matrix of the Cholesky decomposition of the symmetric positive semidefinite matrix $Q(T,\Delta)$.

Before presenting our numerical results, let us discuss two different ways of approaching the joint calibration problem that can be found in the recent literature.

The first approach consists in choosing for instance the first maturity of SPX and VIX to coincide (or differ up to two days, see Remark 1.1), i.e., $T_{1}^{{\textrm{SPX}}}=T_{1}^{{\textrm{VIX}}}$ and then for $j\geq 2$, $T_{j}^{{\textrm{SPX}}}=T_{j-1}^{{\textrm{VIX}}}+\Delta$, see for instance Guyon (2021); Guo et al. (2021); Guyon and Lekeufack (2022).

The second approach is to consider ${\mathcal{T}}^{{\textrm{SPX}}}={\mathcal{T}}^{{\textrm{VIX}}}$, i.e., to choose the same (or close together, see Remark 1.1) maturities both for SPX and VIX options. This perspective has been adopted for instance by Gatheral et al. (2018); Rosenbaum and Zhang (2021); Grzelak (2022); Bondi et al. (2022b).

Both approaches deal with the joint modeling of SPX and VIX options and in order to be consistent with both viewpoints taken in the literature, we show how our signature-based model solves the joint calibration within both settings. For this reason we split the rest of the section in two subsection and discuss them separately.

Here we consider call options for both indices on the trading day 02/06/2021, as in Guyon and Lekeufack (2022). Maturities are reported in the following tables with the corresponding range of strikes (in percentage) with respect to the spot and the market’s futures prices.

We stress that the shortest maturity considered is of 14 days for both SPX and VIX, then the second and third maturity of the SPX are 44 days and 58 days, respectively, and the second one for the VIX is 28 days. Moreover, we consider a high moneyness level (up to 220$\%$) for VIX options, usually rather difficult to fit. Regarding our modeling choice we fix $d=3$, $n=3$ and choose the primary process $X$ to be a three dimensional Ornstein-Uhlenbeck process (see Example 4.5) with parameters

We report also the relative absolute error between the market future prices and the calibrated ones as defined in (5.15):

Calibrate jointly $T_{1}^{{\textrm{SPX}}},T_{1}^{\textrm{VIX}}$ and $T_{2}^{\textrm{SPX}}$.

Use the parameters from the calibration of $T_{j}^{\textrm{SPX}}$ and $T_{j-1}^{\textrm{VIX}}$ to fit jointly the maturities $T_{j+1}^{\textrm{SPX}}$ and $T_{j}^{\textrm{VIX}}$ for $j=2,\dots,J$.

We consider $J=4$, where the last maturity for the SPX is 170 days, and the last maturity for the VIX is 77 days. For the first two maturities of the SPX and the first of the VIX we consider the same moneyness ranges as in Figure 4, hence we specify here only the ranges for the longer maturities:

We observe that for this choice of maturities Assumption 5.10 is satisfied. Hence the second expression for the time-varying VIX is used from Proposition (5.11). On the other hand in order to compute the price of the SPX options in the time-varying case we use the representation of the log-price provided in Proposition 6.10. In (7.1), we employ $\lambda=0.25$ for each calibration within the rolling procedure and we consider always as loss function ${\mathcal{L}}^{\beta}$ as introduced in (5.14) for $\beta=0$. It is worth mentioning that the initial parameter search discussed in Remark 5.7, has been employed for calibrating jointly $T_{1}^{{\textrm{SPX}}},T_{1}^{\textrm{VIX}}$ and $T_{2}^{\textrm{SPX}}$, whereas for the next slices we have considered the previously calibrated parameters as starting point of the optimization.

Finally we report the absolute relative error on the VIX futures’ prices:

1.2 Second approach

Let us now consider the second approach described at the beginning of Section 7.1. Specifically, we consider a unique set of maturities for both SPX and VIX on the trading day of 02/06/2021. For this study, we do not consider time-varying parameters. In the following table we report the moneyness ranges for SPX options in the second row and on the last row the ones for VIX options:

We consider $\lambda=0.5$ and as loss function $\mathcal{L}$ we employ (5.14) with $\beta=1$ for VIX options and the same (without futures) for SPX options.

We additionally report the relative absolute error of the calibrated VIX futures:

Appendix A Numerical results for the Brownian motion case

This appendix is dedicated to the calibration to VIX options only, similarly as in Section 5.4.1, however with the primary process $(X_{t})_{t\geq 0}$ being simply correlated Brownian motions (similarly as in Cuchiero et al. (2022a)) instead of OU-processes.

To be precise, we here model given by (3.1)-(3.2), where $(X_{t})_{t\geq 0}$ is $2$-dimensional Brownian motion. The correlation matrix of $Z=(X,B)$ is specified, as in Section 5.4.1, namely by

For the other parameters we consider a truncation’s level $n=3$, we sample $N_{MC}=80000$ trajectories for Monte Carlo pricing, and we minimize the loss function (5.14) with $\beta=1$ to fit the same data-set as in Section 5.4.1.

We observe that with this specification the model is neither able to calibrate to all future market prices (see Figure 9 below) nor to fit the market implied volatilites accurately. One can indeed see that the model implied volatilities often do not lie within the bid-ask spreads, in particular for high strikes and short maturities.

References