Fast exact joint S&P 500/VIX smile calibration in discrete and continuous time

An arbitrage-free model for exotic options that captures smiles and futures is presented

CLICK HERE TO DOWNLOAD THE PDF

Florian Bourgey and Julien Guyon introduce a novel discrete-time-continuous-time exact calibration method. They build an S&P 500/VIX jointly calibrated discrete-time model, which they later extend to continuous time by martingale interpolation. The benefit of this technique is that both steps can be made much faster than the known methods that calibrate a continuous-time model directly

In previous work, Guyon (2020, 2024) showed how to build a nonparametric discrete-time arbitrage-free model that perfectly matches market data on Standard & Poor’s 500 index value (SPX) futures, SPX options, Chicago Board Options Exchange Volatility Index (VIX) futures and VIX options. The probability distribution is built by minimising the relative entropy with respect to a reference probability measure, and this Schrödinger problem is solved numerically using an extended Sinkhorn algorithm. This provided the first exact solution to this joint calibration problem, a difficult problem (especially for short maturities) that had eluded quants for many years. Jointly calibrating to SPX and VIX futures and options is important to prevent arbitrage and ensure accurate pricing of liquid hedging instruments; calibrating to VIX derivatives means incorporating market information on SPX forward volatilities. Figure 1 showcases the additional information contributed by the VIX by quantifying how model-free bounds for SPX path-dependent payoffs tighten when the prices of VIX futures and VIX options are included. The extra information is also seen from table A, where we compare the prices of several options in models that are all calibrated to SPX smiles but not all calibrated to VIX futures and smiles. Table A shows that to avoid mispricing some payoffs (particularly forward-starting payoffs, which are sensitive to forward volatilities) it is important that the model also fits VIX futures and options, even when the payoff depends only on SPX prices.

The aim of this paper is twofold:

• to speed up the construction of the discrete-time model by turning to Newton-type methods for solving the Schrödinger system (in the next section, we show numerically that a mixed Newton–Sinkhorn method and an implied Newton method converge much faster than the Sinkhorn algorithm); and

• to quickly build a continuous-time extension of the model, allowing the pricing of options depending on SPX values at any date t, while ensuring calibration to the market smiles of SPX and VIX.

The first step draws inspiration from De March (2018), who explored the entropic approximation for discrete-time multi-dimensional martingale optimal transport, excluding VIX, using Newton’s method. Our continuous-time construction in the second step may appear similar to the Bass local volatility of Conze & Henry-Labordère (2022) but it is fundamentally different. First, unlike the construction by Conze & Henry-Labordère (2022), our purely forward Markov functional construction does not require the solution of a fixed-point problem; it is thus much faster. This is because we first build an arbitrage-free multi-marginal discrete-time model consistent with market data. Second, our continuous-time interpolated model fits not only SPX option prices but also VIX market data.

By following both steps, we thus quickly (in less than a minute) build a continuous-time model that is, by construction, exactly calibrated to SPX and VIX smiles and futures. By contrast, other known continuous-time exact solutions to the joint calibration problem (Guo et al 2022; Guyon 2022a) are more involved and demand significantly more computation time. Approximate parametric continuous-time solutions, including those based on rough or rough-like path-dependent volatility models, classical stochastic volatility models or signature-based models, are also costly in terms of computation time; we refer the reader to the extended online version of this article (Bourgey & Guyon 2022) for references. The main benefit of our novel discrete-time-continuous-time calibration method is that both steps are much faster than the known methods that directly calibrate a continuous-time model.

A natural practical application of our continuous-time model is the pricing and hedging of structured products by exotics desks (see table A). With our model, the pricing and hedging of structured products on the SPX indeed take into account all the information given by SPX smiles (the risk-neutral distributions of future SPX values) as well as all the information brought by VIX futures and VIX smiles (the risk-neutral distributions of some future SPX implied volatilities). Once the model has been calibrated (this takes less than a minute per VIX expiry), it is straightforward to implement and use, as it is a Markov functional model that involves simulating only one Brownian motion, along with the VIX at VIX future expiries. The model can also be used for computing reserves and other valuation adjustments and for assessing model risk.

Fast exact joint calibration in discrete time

Setting and notation

Let $T_1>0$ denote a VIX future maturity and set $T_2=T_1+\tau$, where $\tau =30$ days. For simplicity, assume zero interest rates, repos and dividends. We take as given the full market smiles of the SPX index $S$ at $T_1$ and $T_2$ (ie, the full continuum of SPX call prices $C_i(K)$ of maturity $T_i$ for any $i\in \{1,2\}$ and all strikes $K\ge 0$) as well as the full market smile of the VIX index $V$ at $T_1$, (ie, the full continuum of VIX call prices $C_V(K)$ for all strikes $K\ge 0$). For $i\in \{1,2\}$, we use the shorthand notation $S_i:=S_{T_i}$. We use the term forward-starting log contract (FSLC) to denote the financial derivative that pays $-(2/\tau )\ln (S_2/S_1)$ at $T_2$. From the VIX definition (substituting the strip of out-of-the-money options with the log contract for simplicity), the price of the FSLC at $T_1$ is $V^2$.

Figure 1: Model-free bounds of forward-starting calls {\left(S_{T_2}/S_{T_1}-K\right)}_{+}, with and without VIX options data. LB denotes the lower bound. UB denotes the upper bound. LV denotes the local volatility model. LV + Bergomi 2F denotes the two-factor Bergomi stochastic local volatility model, and {\mu }^{*} denotes our discrete minimum-entropy jointly calibrated model (Guyon 2020). We refer the reader to Bourgey & Guyon (2022) and Henry-Labordère (2017) for details on how to compute model-free bounds.

Table A: Prices of various options in the Dupire local volatility model, local two-factor Bergomi model and our continuous-time model. [We choose the same parameters as Guyon (2022b, table 4) for the local two-factor Bergomi, jointly calibrated to the term-structures of SPX ATM skew and VIX${ }^{\text{2}}$ implied volatility: $k_{\text{1}}=\text{21.91}$, $k_{\text{2}}=\text{1.04}$, ${\rho }_{XY}=\text{1}$, ${\rho }_{SX}=-\text{1}$, ${\rho }_{SY}=-\text{1}$, ${\theta }_{\text{1}}=\text{0.77}$, $\omega =\text{6.64}$. We define $M_{t,T}:={\max }_{t\le u\le T}{S}_u$, $A_{t,T}:=(\text{1}/(T-t)){\int }_t^T{S}_{u}du$]

			Our continuous-
Price	LV	LV + Bergomi 2F	time model
${(M_{\text{0},T_{\text{2}}}-S_{\text{0}})}_{+}$	$\text{0}\text{76.81}\pm \text{0.30}$	$\text{0}\text{73.13}\pm \text{0.28}$	$\text{0}\text{75.37}\pm \text{0.36}$
${(M_{T_{\text{1}},T_{\text{2}}}-S_{\text{0}})}_{+}$	$\text{0}\text{70.53}\pm \text{0.32}$	$\text{0}\text{67.20}\pm \text{0.30}$	$\text{0}\text{68.92}\pm \text{0.38}$
${(M_{T_{\text{1}},T_{\text{2}}}-S_{T_{\text{1}}})}_{+}$	$\text{0}\text{61.21}\pm \text{0.28}$	$\text{0}\text{56.16}\pm \text{0.23}$	$\text{0}\text{62.99}\pm \text{0.34}$
$\text{100}\times \frac{M_{T_{\text{1}},T_{\text{2}}}}{S_{T_{\text{1}}}}$	$\text{101.896}\pm \text{0.011}$	$\text{101.707}\pm \text{0.09}\text{0}$	$\text{101.941}\pm \text{0.012}$
${(A_{T_{\text{1}},T_{\text{2}}}-S_{\text{0}})}_{+}$	$\text{0}\text{32.02}\pm \text{0.22}$	$\text{0}\text{31.96}\pm \text{0.22}$	$\text{0}\text{33.22}\pm \text{0.26}$
${(A_{T_{\text{1}},T_{\text{2}}}-S_{T_{\text{1}}})}_{+}$	$\text{0}\text{19.51}\pm \text{0.16}$	$\text{0}\text{19.05}\pm \text{0.15}$	$\text{0}\text{23.14}\pm \text{0.17}$
${\left(\frac{A_{T_{\text{1}},T_{\text{2}}}}{S_{T_{\text{1}}}}-\text{1}\right)}_{+}$	$(\text{7.05}\pm \text{0.06})$	$(\text{6.87}\pm \text{0.05})$	$(\text{8.27}\pm \text{0.06})$
	$\times {\text{10}}^{-\text{3}}$	$\times {\text{10}}^{-\text{3}}$	$\times {\text{10}}^{-\text{3}}$

For each maturity $T_i$, $i\in \{1,2\}$, the absence of static SPX arbitrage is equivalent to the existence of a risk-neutral measure ${\mu }_i:={\partial }^2{C}_i/\partial K^2$, $i\in \{1,2\}$, such that the price of any vanilla option $u_i(\cdot )$ written on $S_i$ is the expectation ${𝔼}^i[u_i(S_i)]:={𝔼}^{{\mu }_i}[u_i(S_i)]$ of the payoff under ${\mu }_i$. Similarly, by the absence of static VIX arbitrage, there exists a risk-neutral measure ${\mu }_V:={\partial }^2{C}_V/\partial K^2$ such that the price of any vanilla option $u_V(\cdot )$ written on $V$ is the expectation ${𝔼}^V[u_V(V)]:={𝔼}^{{\mu }_V}[u_V(V)]$ of the payoff under ${\mu }_V$.

In the absence of dynamic SPX arbitrage (or calendar arbitrage), ${\mu }_1$ and ${\mu }_2$ are of convex order (ie, ${𝔼}^1[f(S_1)]\le {𝔼}^2[f(S_2)]$ for any convex function $f:{ℝ}_{ >0}\to ℝ$), even if we allow trading in the FSLC at $T_1$. By the absence of arbitrage, the price of $S_i$ at time $0$ is the initial SPX spot value $S_0>0$ (ie, ${𝔼}^i[S_i]=S_0$). Furthermore, ${𝔼}^V[V]=F_V\ge 0$, where $F_V$ is the value at time $0$ of the VIX future maturing at $T_1$. Finally, in order for the log contracts and the VIX squared to have finite prices, in the remainder of this section we assume the following.

Assumption 1

For any $i\in \{1,2\}$, the given marginals ${\mu }_1$, ${\mu }_V$ and ${\mu }_2$ satisfy $E^i[S_i]=S_0$,  and $E^V[V]=F_V$, .

Let $𝒳:={ℝ}_{ >0}\times {ℝ}_{ \ge 0}\times {ℝ}_{ >0}$ and define the strictly convex function $L:x\in {ℝ}_{ >0}\mapsto -(2/\tau )\ln x$. For a probability distribution $\rho$ on $ℝ$, we denote the associated cumulative distribution by $F_{\rho }$ (ie, $F_{\rho }(x)=\rho ((-\infty ,x])$ for every $x\in ℝ$). Let $𝒫({ℝ}_{ >0}^2)$ (respectively, $𝒫(𝒳)$) denote the set of all probability measures on ${ℝ}_{ >0}^2$ (respectively, $𝒳$).

Let ${𝒰}^V$ be the set of all measurable functions $u_1,u_2:{ℝ}_{ >0}\to ℝ$, $u_V:{ℝ}_{ \ge 0}\to ℝ$, ${\Delta }_S,{\Delta }_L:{ℝ}_{ >0}\times {ℝ}_{ \ge 0}\to ℝ$ satisfying $u_i\in L^1({\mu }_i)$ for $i\in \{1,V,2\}$, and ${\Delta }_S$, ${\Delta }_L$ bounded. We use the shorthand notation:

	${\Delta }_S^{(S)}(s_1,v,s_2)$	$:={\Delta }_S(s_1,v)(s_2-s_1)$
	${\Delta }_L^{(L)}(s_1,v,s_2)$	$:={\Delta }_L(s_1,v)(L(s_2/s_1)-v^2)$

to denote the profit-and-losses (P&Ls) from delta-hedging at time $T_1$ in the SPX and the log-contract, respectively. Finally, let ${ℳ}_c({\mu }_1,{\mu }_V,{\mu }_2)$ denote the set of all VIX-constrained martingale probability measures:

	${ℳ}_c({\mu }_1,{\mu }_V,{\mu }_2):=\{\mu \in 𝒫(𝒳):S_1\sim 𝜇{\mu }_1,V\sim 𝜇{\mu }_V,S_2\sim 𝜇{\mu }_2,$
	$\mathrm{\hspace{1em}\hspace{1em}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{𝔼}^{\mu }[S_2\mid S_1,V]=S_1,{𝔼}^{\mu }[L(S_2/S_1)|S_1,V]=V^2\}$

Entropy minimisation

Solving the joint calibration problem is equivalent to building a probability measure $\mu \in {ℳ}_c({\mu }_1,{\mu }_V,{\mu }_2)$. In the absence of joint SPX/VIX arbitrage, there may exist an infinite number of models $\mu$ within the convex set ${ℳ}_c({\mu }_1,{\mu }_V,{\mu }_2)$ of jointly calibrated models. To build a specific, meaningful jointly calibrated model, in the spirit of Avellaneda et al (1997), Guyon (2020, 2024) suggests a minimum entropy approach: choose a reasonable (though not jointly calibrated) model $\mu ¯$ on $𝒳$, possibly one derived from a model already in use at the financial institution, and build the probability measure $\mu \in {ℳ}_c({\mu }_1,{\mu }_V,{\mu }_2)$ that is closest to $\mu ¯$ in the entropic sense (ie, $\mu$ has minimum relative entropy with respect to $\mu ¯$):

(M)

From Guyon (2024), if the minimum entropy problem is finite, there exists a unique minimiser ${\mu }^{*}\in {ℳ}_c({\mu }_1,{\mu }_V,{\mu }_2)$ of the form:

	$$\frac{d{\mu }^{*}}{d\mu ¯}=e_u(S_1,V,S_2)$$
	$$e_u(s_1,v,s_2):=e^{u_1(s_1)+u_V(v)+u_2(s_2)+{\Delta }_S^{(S)}(s_1,v,s_2)+{\Delta }_L^{(L)}(s_1,v,s_2)}$$

where $u:=(u_1,u_V,u_2,{\Delta }_S,{\Delta }_L)$, if they exist, are maximisers (called Schrödinger potentials) of the dual problem:

(P)

Additionally, in any case, $D_{\mu }=P_{\mu }$. Note that (P) is an unconstrained concave maximisation problem. Both problems are dual to each other; following a terminology proposed by Dupire, (M) is a measure problem while (P) is a portfolio problem.

As, in practice, only a finite number of SPX and VIX vanilla options are available for trading, we consider vanilla payoffs $u_1$, $u_V$ and $u_2$ that are linear combinations of finitely many call options, along with one position in the bond, one position in $S_1$ and one position in the VIX futures. Therefore, we consider a market data set $𝒦$ composed of call options on $S_1$, $V$ and $S_2$, which we denote by ${(C_K^1)}_{K\in {𝒦}_1}$, ${(C_K^V)}_{K\in {𝒦}_V}$ and ${(C_K^2)}_{K\in {𝒦}_2}$ with respective strikes ${𝒦}_1$, ${𝒦}_V$ and ${𝒦}_2$, and we build a model of the form:

	$$\frac{d{\mu }_{𝒦,\theta }}{d\mu ¯}=e_{\theta }(S_1,V,S_2)$$

where $\theta$ is an element of the set $\Theta$ of all $\theta :=(c,{\Delta }_S^0,{\Delta }_V^0,a^1,a^V,a^2,{\Delta }_S,{\Delta }_L)$ such that $c,{\Delta }_S^0,{\Delta }_V^0\in ℝ$, $a^1\in {ℝ}^{{𝒦}_1}$, $a^V\in {ℝ}^{{𝒦}_V}$, $a^2\in {ℝ}^{{𝒦}_2}$ and ${\Delta }_S,{\Delta }_L:{ℝ}_{ >0}\times {ℝ}_{ \ge 0}\to ℝ$ are bounded measurable functions of $(s_1,v)$. The measure ${\mu }_{𝒦,\theta }$ is then a consistent, arbitrage-free model that is jointly calibrated to the market prices of SPX/VIX futures and options if and only if $\theta$ solves the so-called $𝒦$-Schrödinger system:

(1)

The first equation states that ${\mu }_{𝒦,\theta }$ is a probability measure, while the others state that it belongs to the set ${ℳ}_{c,𝒦}({\mu }_1,{\mu }_V,{\mu }_2)$ of probability measures $\mu$ satisfying:

	$${𝔼}^{\mu }[S_1]=S_0,{𝔼}^{\mu }[V]=F_V,\forall K\in {𝒦}_1,{𝔼}^{\mu }[{(S_1-K)}_{+}]=C_K^1$$
	$$\forall K\in {𝒦}_V,{𝔼}^{\mu }[{(V-K)}_{+}]=C_K^V,\forall K\in {𝒦}_2,{𝔼}^{\mu }[{(S_2-K)}_{+}]=C_K^2$$
	$${𝔼}^{\mu }[S_2\mid S_1,V]=S_1,{𝔼}^{\mu }[L(S_2/S_1)\mid S_1,V]=V^2$$

Remark 1

The value of the prior measure $\mu ¯$ is chosen by the modeller. Examples include a lognormal prior:

$$\mu ¯(d{s}_1,dv,d{s}_2)=\nu (d{s}_1,dv)T(s_1,v,d{s}_2)$$

where $\nu ={\mu }_1\otimes {\mu }_V$ and $T(s_1,v,d{s}_2)$ is the distribution of:

$$s_1\exp (v\sqrt{\tau }G-\frac{1}{2}{v}^2\tau )$$

with $G\sim 𝒩(0,1)$;¹¹ 1 Assuming that $S_2$ is lognormally conditioned on $S_1$ and that $V$ is financially natural, which may not be the best choice in practice. Indeed, in this case, ${𝔼}^{\mu }[e^{\delta S_2}\mid S_1,V]=+\infty$ for $\delta >0$, so the expectations in (1) may not be well defined. In practice, we avoid those integrability issues by working with a finite-support approximation of $\mu ¯$ stemming from the Gaussian quadrature approximation of the expectations (see remark 2). or the independent prior (product measure) $\mu ¯={\mu }_1\otimes {\mu }_V\otimes {\mu }_2$.

Remark 2

To evaluate the expectations arising in (1), we use a Gauss–Legendre quadrature when integrating with respect to $s_1$ (respectively, $v$) with grid ${𝒢}_1:=\{s_1^{(1)}\le \dots \le s_1^{(n_1)}\}$ (respectively, ${𝒢}_V:=\{v^{(1)}\le \dots \le v^{(n_V)}\}$). For $s_2$, in the case of the lognormal prior, we use a Gauss–Hermite quadrature with knots $\{z^{(1)},\dots ,z^{(n_2)}\}.$

Solving the Schrödinger system (1)

The Sinkhorn algorithm

The classical method for solving Schrödinger systems is the Sinkhorn algorithm, an iterative method that sequentially solves the individual equations in the Schrödinger system (here, (1)), which converge to the optimiser ${\theta }^{*}$ of the whole system. This algorithm has recently gained popularity in machine learning, where it is used to quickly compute Wasserstein distances and more generally to solve optimal transport problems, via a small entropic penalty. It has also been applied to martingale optimal transport problems (De March 2018) and, in particular, in quantitative finance to quickly build arbitrage-free smiles (De March & Henry-Labordère 2019). In Guyon (2024) the Sinkhorn algorithm was extended to accommodate the martingality and consistency constraints in (1), and shown to converge toward a jointly calibrated model. However, the convergence was somewhat slow (see the ‘Comparison of the different algorithms’ section below). In the following sections, we present two faster alternatives for solving (1) numerically; both rely on solving the portfolio problem (P).

The Newton–Sinkhorn algorithm

Observe that if we define the concave function:

	$J_{\mu ¯,𝒦}(\theta )$	$:=c+{\Delta }_S^0{S}_0+{\Delta }_V^0{F}_V$
		$\mathrm{\hspace{1em}\hspace{1em}}+\sum _{i\in \{1,V,2\}}\sum _{K\in {𝒦}_i}{a}_K^i{C}_K^i-{𝔼}^{\mu }[e_{\theta }(S_1,V,S_2)]+1$

solving the $𝒦$-Schrödinger system is equivalent to cancelling the gradient of $J_{\mu ¯,𝒦}$. Hence, to solve the system and build ${\mu }^{*}$, one can directly solve the portfolio problem:

$$P_{\mu ¯,𝒦}:=\underset{\theta \in \Theta }{sup}{J}_{\mu ¯,𝒦}(\theta )$$

(2)

which is the ‘finitely many payoffs’ version of (P). To this end, we suggest the following Newton–Sinkhorn algorithm. Each iteration involves a Newton step followed by a Sinkhorn step.

Newton step. Starting from an initial guess ${\theta }^{(0)}$, we first solve for every iteration $n\in \mathbb{N}$ and $(s_1,v)\in {𝒢}_1\times {𝒢}_V$ the portfolio problem (2):

${\theta }^{-\Delta ,(n+1)}$

$=\underset{{\theta }^{-\Delta }\in {\Theta }^{-\Delta }}{\arg \max }{J}_{\mu ¯,𝒦}({\theta }^{-\Delta },{\Delta }_S^{(n)}(s_1,v),{\Delta }_L^{(n)}(s_1,v))$

(3)

where:

$${\theta }^{-\Delta }:=(c,{\Delta }_S^0,{\Delta }_V^0,a^1,a^V,a^2)$$

Since the Hessian of $J_{\mu ¯,𝒦}({\theta }^{-\Delta },{\Delta }_S^{(n)}(s_1,v),{\Delta }_L^{(n)}(s_1,v))$ is known in closed form, this step is extremely fast. We solve (3) using the function scipy.optimize.minimize (method="trust-exact") from the scipy library.

Sinkhorn step. Then, for all $(s_1,v)\in {𝒢}_1\times {𝒢}_V$, we jointly solve for ${\Delta }_S^{(n+1)}(s_1,v)$ and ${\Delta }_L^{(n+1)}(s_1,v)$ the two-dimensional nonlinear system:

	$f_{s_1,v}({\Delta }_S^{(n+1)}(s_1,v),{\Delta }_L^{(n+1)}(s_1,v),a^{2,(n+1)})$	$=0,$		(4)
	$g_{s_1,v}({\Delta }_S^{(n+1)}(s_1,v),{\Delta }_L^{(n+1)}(s_1,v),a^{2,(n+1)})$	$=0,$		(5)

where $a^{2,(n+1)}$ is the optimal vector $a^2$ from the previous step, (3), and for all $(x,y)\in {ℝ}^2$:

	$f_{s_1,v}(x,y,a^2)$
	$:={\int }_{ℝ}(s_2-s_1)e^{{\sum }_{K\in {𝒦}_2}{a}_K^2{(s_2-K)}_{+}+x(s_2-s_1)+y(L(\frac{s_2}{s_1})-v^2)}$
	$\mathrm{\hspace{1em}\hspace{1em}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \mu ¯(s_1,v,d{s}_2)$
	$g_{s_1,v}(x,y,a^2)$
	$:={\int }_{ℝ}\left(L\left(\frac{s_2}{s_1}\right)-v^2\right)e^{{\sum }_{K\in {𝒦}_2}{a}_K^2{(s_2-K)}_{+}+x(s_2-s_1)+y(L(\frac{s_2}{s_1})-v^2)}$
	$\mathrm{\hspace{1em}\hspace{1em}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \mu ¯(s_1,v,d{s}_2)$

This is the Sinkhorn step, where the last two equations of (1) are both solved.

We use the Levenberg–Marquardt algorithm to solve (4) and (5) via scipy.optimize.root(method="lm"). A full Newton algorithm with parameterised deltas has also been considered in the extended online version of this article (Bourgey & Guyon 2022).

The implied Newton algorithm

Inspired by De March (2018), we observe that:

	${\theta }^{*}$	$=\underset{\theta \in \Theta }{\arg \max }{J}_{\mu ¯,𝒦}(\theta )$
		$=\underset{{\theta }^{-\Delta }\in {\Theta }^{-\Delta }}{\arg \max }{J}_{~}({\theta }^{-\Delta },{\Delta }_S^{*}(\cdot ,\cdot ,a^2),{\Delta }_L^{*}(\cdot ,\cdot ,a^2))$

where $({\Delta }_S^{*}(\cdot ,\cdot ,a^2),{\Delta }_L^{*}(\cdot ,\cdot ,a^2))$ solves the two-dimensional nonlinear system:

	$f_{s_1,v}({\Delta }_S^{*}(\cdot ,\cdot ),{\Delta }_L^{*}(\cdot ,\cdot ),a^2)$	$=0,$
	$g_{s_1,v}({\Delta }_S^{*}(\cdot ,\cdot ),{\Delta }_L^{*}(\cdot ,\cdot ),a^2)$	$=0.$

That is, for each ${\theta }^{-\Delta }$, we first optimise over ${\Delta }_S(\cdot ,\cdot )$ and ${\Delta }_L(\cdot ,\cdot )$, and then we optimise over ${\theta }^{-\Delta }$. Note that the inner optimisation depends on ${\theta }^{-\Delta }$ only through $a^2$.

Figure 2: Comparison of the performance of the Newton–Sinkhorn, implied Newton (with a 10-iteration pure Sinkhorn warm-start) and Sinkhorn algorithms. We plot the common logarithm of the calibration error as a function of the computational time as of April 20, 2020, with T_{\text{1}}=\text{30} days.

As with $J_{\mu ¯,𝒦}$, the gradient and Hessian of $J_{~}$ are known in closed form. In fact, $J_{~}$ has the same gradient and Hessian as $J_{\mu ¯,𝒦}$, except for the terms involving differentiation with respect to $a^2$, whose exact computation is detailed in Bourgey & Guyon (2022).

Comparison of the different algorithms

All the numerical tests were performed on a MacBook Pro laptop with a 2.6 GHz six-core Intel Core i7 processor and 32 GB memory using Python 3.9.6.

To compare the various algorithms, we plot the logarithm (in base 10) of the calibration error for the SPX smiles at $T_1,T_2$ and the VIX smile at $T_1$ as a function of the computational time as of April 20, 2020 ($T_1=30$ days); see figure 2. Other calibration dates in 2018 and 2022 are considered by Bourgey & Guyon (2022). The calibration error is computed as the sum of the absolute relative errors of the three futures, the mean of the three smiles and the total mass of ${\mu }^{*}$, which must be a probability measure.

We choose the lognormal prior for the reference measure $\mu ¯$ (see remark 1) and we take $n_1=n_v=45$ knots for the integration with respect to $s_1$ and $v$, and $n_2=25$ knots for the integration with respect to $s_2$. We set $s_1^{(1)}=F_{{\mu }_1}^{-1}(q)$, $v^{(1)}=F_{{\mu }_V}^{-1}(q)$ and $s_1^{(n_1)}=F_{{\mu }_1}^{-1}(1-q)$, $v^{(n_V)}=F_{{\mu }_V}^{-1}(1-q)$ with $q={10}^{-3}$ for the lowest and highest values of the quadrature grids of $s_1$ and $v$.

We compare the calibration speed of the Sinkhorn (S), Newton–Sinkhorn (NS) and implied Newton (IN) algorithms. For IN, as suggested by De March (2018, section 7.1), we initialise the algorithm with 10 iterations of a (pure) S algorithm (warm-start procedure). A warm-start initialisation of NS yielded no improvement.

Figure 3: Futures and smiles of S_{\text{1}}, V, S_{\text{2}} after 35 seconds in the discrete-time model calibrated with the implied Newton algorithm versus market smiles as of April 20, 2020, with T_{\text{1}}=\text{30} days. (a) Smile of SPX, T_{\text{1}}=\text{30} days. (b) Smile of VIX, T_{\text{1}}=\text{30} days. (c) Smile of SPX, T_{\text{2}}=\text{60} days.

Figure 4: SPX smiles at T_{\text{1}} and T_{\text{2}} and VIX smile at T_{\text{1}} (market, discrete-time model computed using the implied Newton algorithm and continuous-time extension), as of August 1, 2018. (a) Smile of SPX, T_{\text{1}}=\text{21} days. (b) Smile of VIX, T_{\text{1}}=\text{21} days. (c) Smile of SPX, T_{\text{2}}=\text{51} days.

In figure 2, we observe that IN is the fastest, and S the slowest. In figure 3, we plot the futures and smiles for $S_1$, $V$ and $S_2$ obtained with the IN algorithm after 60 seconds as well as the market smiles: the fits are perfect. The martingality and VIX constraints (reported in Bourgey & Guyon (2022)) are also perfectly satisfied. Similar plots are obtained for the other two calibration dates by Bourgey & Guyon (2022).

We also tested the three methods when the reference measure is the product measure (independent prior). The numerical results are reported in Bourgey & Guyon (2022). We observe that the three methods seem to be less stable with the independent prior and typically require a higher number of nodes; we chose $n_1=n_V=n_2=45$. With the independent prior, IN needs more steps to converge.

Continuous-time extension

We now extend the discrete-time model ${\mu }^{*}$ (and the probability space) to build a continuous-time model for ${(S_t)}_{t\in [0,T_2]}$. The model will also include the VIX at $T_1$, $V$. That is, we build a probability $\mathbb{P}$ on $C([0,T_2],ℝ)\times {ℝ}_{+}$ representing the distribution of $({(S_t)}_{t\in [0,T_2]},V)$. Here, $V$ plays the role of a discrete-time stochastic volatility, representing the stochastic volatility anticipated at $T_1$ for the $[T_1,T_2]$ period, but our model involves no continuous-time stochastic volatility process. Similarly to that of Conze & Henry-Labordère (2022), our model is computationally efficient, as it only requires simulating one Brownian motion (and $V$).

The key advantage of our construction compared with the Bass local volatility of Conze & Henry-Labordère (2022) is that it starts directly from a joint discrete distribution of $(S_1,S_2)$, making the continuous-time interpolation a purely forward construction with no need to solve a fixed-point problem. Moreover, it includes a stochastic volatility component, $V$, for the calibration to VIX futures and VIX smiles in addition to SPX smiles. As a result, in contrast with Conze & Henry-Labordère (2022) and in line with the path-dependency observed in financial markets (Guyon & Lekeufack 2023), our model is path-dependent: the SPX dynamics after $T_1$ depend on both $S_1$ and $V$.

Step 1: simulation of {\left(S_t\right)}_{t\in [0,T_1]}

We want ${(S_t)}_{t\in [0,T_1]}$ to be a $\mathbb{P}$-martingale and $S_1$ to have distribution ${\mu }_1$ under $\mathbb{P}$. To achieve this, one possible choice is to use a Markov functional model $S_t=u(t,W_t)$, $t\in [0,T_1]$, where $W$ is a $\mathbb{P}$-Brownian motion and $u$ satisfies the heat equation:

$${\partial }_{t}u+\frac{1}{2}{\partial }_x^{2}u=0$$

Further, define the filtration:

In such a case, ${(S_t)}_{t\in [0,T_1]}$ is an $(({ℱ}_t),\mathbb{P})$-martingale, and the terminal condition $u(T_1,\cdot )=g$ is determined via quantiles so that $u(T_1,W_{T_1})$ has distribution ${\mu }_1$, ie, for all $x\in ℝ$, we set:

$$g(x):=F_{{\mu }_1}^{-1}\left(F_{𝒩(0,1)}\left(\frac{x}{\sqrt{T_1}}\right)\right).$$

The solution $u$ to the heat equation is explicit and given by:

	$u(t,x)$	$=𝔼[g(W_{T_1})\mid W_t=x]$
		$=𝔼[g(x+W_{T_1}-W_t)]=(g\ast K_{T_1-t})(x)$		(6)

where $\ast$ is the convolution operator and where, for $t>0$, $K_t(x)=e^{-x^2/2t}/\sqrt{2\pi t}$ is the heat kernel.

Step 2: simulation of V given {\left(S_t\right)}_{t\in [0,T_1]}

At this stage, we have simulated ${(S_t)}_{t\in [0,T_1]}$ such that ${(S_t)}_{t\in [0,T_1]}$ is an $(({ℱ}_t),\mathbb{P})$-martingale and $S_1$ has distribution ${\mu }_1$ under $\mathbb{P}$. Now, we simulate $V$ given $\sigma ({(W_s)}_{s\in [0,T_1]})$ under $\mathbb{P}$ as follows: the distribution of $V$ given $\sigma ({(W_s)}_{s\in [0,T_1]})$ under $\mathbb{P}$ is assumed to depend only on $S_1$ and is taken equal to the distribution of $V$ given $S_1$ under ${\mu }^{*}$. Since $S_1$ has distribution ${\mu }_1$ under both ${\mu }^{*}$ and $\mathbb{P}$, this means that the distribution of $(S_1,V)$ is the same under ${\mu }^{*}$ and $\mathbb{P}$; in particular, $V$ has the distribution ${\mu }_V$ under $\mathbb{P}$.

Step 3: simulation of ${(S_t)}_{t\in [T_1,T_2]}$ given

$ℱ$${ }_{{𝑻}_{\mathrm{𝟏}}}$.  In this last step, we build dynamics for ${(S_t)}_{t\in [T_1,T_2]}$ conditional on ${ℱ}_{T_1}$ such that ${(S_t)}_{t\in [T_1,T_2]}$ is an $(({ℱ}_t),\mathbb{P})$-martingale starting from $S_1$, and $S_2$ has distribution ${\mu }_2$. We use once again a Markov functional construction. Given $S_1=s$ and $V=v$, we consider:

$$g_{s,v}(x):=F_{{\mu }_{2|s,v}}^{-1}\left(F_{𝒩(0,1)}\left(\frac{x}{\sqrt{\tau }}\right)\right)$$

for all $x\in ℝ$ and where ${\mu }_{2|s,v}$ is the distribution of $S_2$ given $S_1=s$ and $V=v$ under ${\mu }^{*}$. Then, given ${ℱ}_{T_1}$, we define for $t\in (T_1,T_2]$:

$$S_t:=u_{S_1,V}(t,W_t-W_{T_1})$$

where for every $s,v>0$:

	$u_{s,v}(t,x)$	$=𝔼[g_{s,v}(W_{T_2}-W_{T_1})\mid W_t-W_{T_1}=x]$
		$=𝔼[g_{s,v}(x+W_{T_2}-W_t)]=(g_{s,v}*K_{T_2-t})(x)$		(7)

It is easy to check that ${(S_t)}_{t\in [T_1,T_2]}$ is an $(({ℱ}_t),\mathbb{P})$-martingale starting from $S_1$. Moreover, the distribution of $S_2$ given $(S_1,V)$ is the same under ${\mu }^{*}$ and $\mathbb{P}$. As the distribution of $(S_1,V)$ is the same under ${\mu }^{*}$ and $\mathbb{P}$, we conclude that the distribution of $(S_1,V,S_2)$ is the same under ${\mu }^{*}$ and $\mathbb{P}$. In particular $S_1$, $V$ and $S_2$ have distributions ${\mu }_1$, ${\mu }_V$ and ${\mu }_2$ under $\mathbb{P}$. We have thus built a model $\mathbb{P}$ on $({(S_t)}_{t\in [0,T_2]},V)$ such that (a) $S_1$, $V$ and $S_2$ have distributions ${\mu }_1$, ${\mu }_V$ and ${\mu }_2$ under $\mathbb{P}$; (b) ${(S_t)}_{t\in [0,T_2]}$ is an $(({ℱ}_t),\mathbb{P})$-martingale; and (c) $V$ is the VIX at $T_1$, since by construction:

	${𝔼}^{\mathbb{P}}\left[L\left(\frac{S_2}{S_1}\right)\right|{ℱ}_{T_1}]$	$={𝔼}^{\mathbb{P}}\left[L\left(\frac{S_2}{S_1}\right)\right|S_1,V]$
		$={𝔼}^{{\mu }^{*}}\left[L\left(\frac{S_2}{S_1}\right)\right|S_1,V]=V^2$

Remark 3 (Extension to several VIX maturities)

Note that this approach can easily be iterated on intervals $[T_i,T_{i+1}]$. For example, after step 3, we will have generated $S_2\sim {\mu }_2$. Then, we just need to generate $V_{T_2}\sim {\mu }_{V_{T_2}}$ and repeat the same procedure. Here, we disregard the Wednesday/Friday issue and make the approximation that the VIX future maturities are exactly 30 days apart.

The numerical implementation of the continuous-time extension is easy; for details we refer the reader to Bourgey & Guyon (2022). The resulting smiles (market, discrete and continuous time) for the SPX at $T_1$ and $T_2$ and that of the VIX at $T_1$ are displayed in figure 4; they were computed by Monte Carlo simulation with ${10}^5$ paths.

Pricing

Our continuous-time model can be used to price path-dependent options on the SPX with the guarantee that the model exactly matches the SPX smiles at $T_1$ and $T_2$ and the VIX future and VIX smile at $T_1$, thus taking into account information about the forward volatility at $T_1$ that is not included in the SPX smiles. As a pricing exercise, we compare it with models commonly used by practitioners: the Dupire local volatility model (Dupire 1994) and the local stochastic version of the two-factor Bergomi model (Bergomi 2005); those models are calibrated to the full SPX implied volatility surface but not to VIX smiles. We consider spot-starting, forward-starting and mixed versions of lookback and Asian options. The prices are reported in table A along with their 95% confidence interval. We used ${10}^5$ Monte Carlo paths and a trapezoidal rule to approximate the time integral for the Asian option. Note that for forward-starting options the price in our jointly calibrated model is always larger than in the other two models: ignoring the VIX information leads to underpricing these forward-starting payoffs.

Conclusion

In this article, we have:

• improved model-free bounds on SPX options by incorporating VIX options data;

• built the minimum-entropy jointly calibrated discrete-time model {\mu }^{*} very fast using an implied Newton method;

• seamlessly extended this discrete-time model to continuous time in a purely forward fashion, using Markov functionals.

Thus, we have established a swift process for creating an arbitrage-free continuous-time model for SPX that accurately calibrates to SPX smiles, VIX futures and VIX smiles. Such a model can be used for pricing and hedging exotic options, computing reserves or valuation adjustments, and assessing model risk. Our main methodological contribution is that we first build a jointly calibrated discrete-time model $(S_{T_i})$ (where the $T_i$ are the calibrated maturities) that is later extended to continuous time using an arbitrage-free martingale time interpolation. Since the discrete-time model can be exactly calibrated much faster than continuous-time models, and since extremely fast extrapolations exist, this novel approach seems to be a promising new avenue for calibrating models.

Florian Bourgey is a quantitative researcher at Bloomberg in New York, while Julien Guyon is a professor of applied mathematics at École des Ponts ParisTech in Paris. They appreciate the valuable feedback from the anonymous referees.

Email: [e-mail], [e-mail].

References

• Avellaneda M, C Friedman, R Holmes and D Samperi, 1997 
  Calibrating volatility surfaces via relative-entropy minimization 
  Applied Mathematical Finance 4(1), pages 37–64 
• Bergomi L, 2005 
  Smile dynamics II 
  Risk October, pages 67–73 
• Bourgey F and J Guyon, 2022 
  Fast exact joint S&P 500/VIX smile calibration in discrete and continuous time 
  Preprint, SSRN 4315084 
• Conze A and P Henry-Labordère, 2022 
  A new fast local volatility model 
  Risk April, http://www.risk.net/7944906 
• De March H, 2018 
  Entropic approximation for multi-dimensional martingale optimal transport 
  Preprint, arXiv:1812.11104 
• De March H and P Henry-Labordère, 2019 
  Building arbitrage-free implied volatility: Sinkhorn’s algorithm and variants 
  Preprint, SSRN 3326486 
• Dupire B, 1994 
  Pricing with a smile 
  Risk January, pages 18–20 
• Guo I, G Loeper, J Obłój and S Wang, 2022 
  Joint modeling and calibration of SPX and VIX by optimal transport 
  SIAM Journal on Financial Mathematics 13(1), pages 1–31 
• Guyon J, 2020 
  The joint S&P 500/VIX smile calibration puzzle solved 
  Risk April, http://www.risk.net/7518926 
• Guyon J, 2022a 
  Dispersion-constrained martingale Schrödinger bridges: joint entropic calibration of stochastic volatility models to S&P 500 and VIX smiles 
  Preprint, SSRN 3853237 
• Guyon J, 2022b 
  The VIX future in Bergomi models: fast approximation formulas and joint calibration with S&P 500 skew 
  SIAM Journal on Financial Mathematics 13(4), pages 1,418–1,485 
• Guyon J, 2024 
  Dispersion-constrained martingale Schrödinger problems and the exact joint S&P 500/VIX smile calibration puzzle 
  Finance and Stochastics 28(1), pages 27–79 
• Guyon J and J Lekeufack, 2023 
  Volatility is (mostly) path-dependent 
  Quantitative Finance 23(9), pages 1,221–1,258 
• P Henry-Labordère, 2017 
  Model-Free Hedging: A Martingale Optimal Transport Viewpoint 
  Chapman and Hall/CRC