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


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 & Poors 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 T1>0 denote a VIX future maturity and set T2=T1+τ, where τ=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 T1 and T2 (ie, the full continuum of SPX call prices Ci(K) of maturity Ti for any i{1,2} and all strikes K0) as well as the full market smile of the VIX index V at T1, (ie, the full continuum of VIX call prices CV(K) for all strikes K0). For i{1,2}, we use the shorthand notation Si:=STi. We use the term forward-starting log contract (FSLC) to denote the financial derivative that pays -(2/τ)ln(S2/S1) at T2. 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 T1 is V2.

Figure 1: Model-free bounds of forward-starting calls (ST2/ST1-K)+, 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 μ* 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.

Risk 0224 Guyon technical fig 1

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 2 implied volatility: k1=21.91k2=1.04ρXY=1ρSX=-1ρSY=-1θ1=0.77ω=6.64. We define Mt,T:=maxtuTSuAt,T:=(1/(T-t))tTSudu]

   Our continuous-
PriceLVLV + Bergomi 2Ftime model
  ×10-3 ×10-3 ×10-3

For each maturity Tii{1,2}, the absence of static SPX arbitrage is equivalent to the existence of a risk-neutral measure μi:=2Ci/K2i{1,2}, such that the price of any vanilla option ui() written on Si is the expectation 𝔼i[ui(Si)]:=𝔼μi[ui(Si)] of the payoff under μi. Similarly, by the absence of static VIX arbitrage, there exists a risk-neutral measure μV:=2CV/K2 such that the price of any vanilla option uV() written on V is the expectation 𝔼V[uV(V)]:=𝔼μV[uV(V)] of the payoff under μV.

In the absence of dynamic SPX arbitrage (or calendar arbitrage), μ1 and μ2 are of convex order (ie, 𝔼1[f(S1)]𝔼2[f(S2)] for any convex function f: >0), even if we allow trading in the FSLC at T1. By the absence of arbitrage, the price of Si at time 0 is the initial SPX spot value S0>0 (ie, 𝔼i[Si]=S0). Furthermore, 𝔼V[V]=FV0, where FV is the value at time 0 of the VIX future maturing at T1. 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{1,2}, the given marginals μ1μV and μ2 satisfy Ei[Si]=S0Ei[|lnSi|]< and EV[V]=FVEV[V2]<.

Let 𝒳:= >0× 0× >0 and define the strictly convex function L:x >0-(2/τ)lnx. For a probability distribution ρ on , we denote the associated cumulative distribution by Fρ (ie, Fρ(x)=ρ((-,x]) for every x). Let 𝒫( >02) (respectively, 𝒫(𝒳)) denote the set of all probability measures on  >02 (respectively, 𝒳).

Let 𝒰V be the set of all measurable functions u1,u2: >0uV: 0ΔS,ΔL: >0× 0 satisfying uiL1(μi) for i{1,V,2}, and ΔSΔL bounded. We use the shorthand notation:

 ΔS(S)(s1,v,s2) :=ΔS(s1,v)(s2-s1) 
 ΔL(L)(s1,v,s2) :=ΔL(s1,v)(L(s2/s1)-v2) 

to denote the profit-and-losses (P&Ls) from delta-hedging at time T1 in the SPX and the log-contract, respectively. Finally, let c(μ1,μV,μ2) denote the set of all VIX-constrained martingale probability measures:


Entropy minimisation

Solving the joint calibration problem is equivalent to building a probability measure μc(μ1,μV,μ2). In the absence of joint SPX/VIX arbitrage, there may exist an infinite number of models μ within the convex set c(μ1,μV,μ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 μ¯ on 𝒳, possibly one derived from a model already in use at the financial institution, and build the probability measure μc(μ1,μV,μ2) that is closest to μ¯ in the entropic sense (ie, μ has minimum relative entropy with respect to μ¯):

 Dμ¯ :=infμc(μ1,μV,μ2)H(μμ¯)H(μμ¯) :={𝔼μ[ln(dμdμ¯)]=𝔼μ¯[dμdμ¯ln(dμdμ¯)]if μμ¯+otherwise} (M)

From Guyon (2024), if the minimum entropy problem is finite, there exists a unique minimiser μ*c(μ1,μV,μ2) of the form:


where u:=(u1,uV,u2,ΔS,ΔL), if they exist, are maximisers (called Schrödinger potentials) of the dual problem:

  Pμ¯:=supu𝒰VJμ¯(u) Jμ¯(u) :=𝔼1[u1(S1)]+𝔼V[uV(V)] +𝔼2[u2(S2)]-𝔼μ¯[eu(S1,V,S2)]+1} (P)

Additionally, in any case, Dμ¯=Pμ¯. 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 u1uV and u2 that are linear combinations of finitely many call options, along with one position in the bond, one position in S1 and one position in the VIX futures. Therefore, we consider a market data set 𝒦 composed of call options on S1V and S2, which we denote by (CK1)K𝒦1(CKV)K𝒦V and (CK2)K𝒦2 with respective strikes 𝒦1𝒦V and 𝒦2, and we build a model of the form:

 eθ(s1,v,s2):=exp(c+ΔS0s1+ΔV0v+K𝒦1aK1(s1-K)+ +K𝒦VaKV(v-K)++K𝒦2aK2(s2-K)+ +ΔS(S)(s1,v,s2)+ΔL(L)(s1,v,s2)) 

where θ is an element of the set Θ of all θ:=(c,ΔS0,ΔV0,a1,aV,a2,ΔS,ΔL) such that c,ΔS0,ΔV0a1𝒦1aV𝒦Va2𝒦2 and ΔS,ΔL: >0× 0 are bounded measurable functions of (s1,v). The measure μ𝒦,θ 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 θ solves the so-called 𝒦-Schrödinger system:

  𝔼μ¯[dμ𝒦,θdμ¯]=1,𝔼μ¯[S1dμ𝒦,θdμ¯]=S0,𝔼μ¯[Vdμ𝒦,θdμ¯]=FV 𝔼μ¯[(S1-K)+dμ𝒦,θdμ¯]=CK1,K𝒦1 𝔼μ¯[(V-K)+dμ𝒦,θdμ¯]=CKV,K𝒦V 𝔼μ¯[(S2-K)+dμ𝒦,θdμ¯]=CK2,K𝒦2 𝔼μ¯[(S2-S1)dμ𝒦,θdμ¯|S1=s1,V=v]=0 𝔼μ¯[(L(S2S1)-V2)dμ𝒦,θdμ¯|S1=s1,V=v]=0 s1>0,v0} (1)

The first equation states that μ𝒦,θ is a probability measure, while the others state that it belongs to the set c,𝒦(μ1,μV,μ2) of probability measures μ satisfying:

Remark 1

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


where ν=μ1μV and T(s1,v,ds2) is the distribution of:


with G𝒩(0,1);11 1 Assuming that S2 is lognormally conditioned on S1 and that V is financially natural, which may not be the best choice in practice. Indeed, in this case, 𝔼μ¯[eδS2S1,V]=+ for δ>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 μ¯ stemming from the Gaussian quadrature approximation of the expectations (see remark 2). or the independent prior (product measure) μ¯=μ1μVμ2.

Remark 2

To evaluate the expectations arising in (1), we use a Gauss–Legendre quadrature when integrating with respect to s1 (respectively, v) with grid 𝒢1:={s1(1)s1(n1)} (respectively, 𝒢V:={v(1)v(nV)}). For s2, in the case of the lognormal prior, we use a Gauss–Hermite quadrature with knots {z(1),,z(n2)}.

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 θ* 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μ¯,𝒦(θ) :=c+ΔS0S0+ΔV0FV 

solving the 𝒦-Schrödinger system is equivalent to cancelling the gradient of Jμ¯,𝒦. Hence, to solve the system and build μ*, one can directly solve the portfolio problem:

 Pμ¯,𝒦:=supθΘJμ¯,𝒦(θ) (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 θ(0), we first solve for every iteration n and (s1,v)𝒢1×𝒢V the portfolio problem (2):

 θ-Δ,(n+1) =argmaxθ-ΔΘ-ΔJμ¯,𝒦(θ-Δ,ΔS(n)(s1,v),ΔL(n)(s1,v)) (3)



Since the Hessian of Jμ¯,𝒦(θ-Δ,ΔS(n)(s1,v),ΔL(n)(s1,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 (s1,v)𝒢1×𝒢V, we jointly solve for ΔS(n+1)(s1,v) and ΔL(n+1)(s1,v) the two-dimensional nonlinear system:

 fs1,v(ΔS(n+1)(s1,v),ΔL(n+1)(s1,v),a2,(n+1)) =0, (4)
 gs1,v(ΔS(n+1)(s1,v),ΔL(n+1)(s1,v),a2,(n+1)) =0, (5)

where a2,(n+1) is the optimal vector a2 from the previous step, (3), and for all (x,y)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:

 θ* =argmaxθΘJμ¯,𝒦(θ) 

where (ΔS*(,,a2),ΔL*(,,a2)) solves the two-dimensional nonlinear system:

 fs1,v(ΔS*(,),ΔL*(,),a2) =0, 
 gs1,v(ΔS*(,),ΔL*(,),a2) =0. 

That is, for each θ-Δ, we first optimise over ΔS(,) and ΔL(,), and then we optimise over θ-Δ. Note that the inner optimisation depends on θ-Δ only through a2.

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 T1=30 days.

Risk 0224 Guyon technical fig 2

As with Jμ¯,𝒦, the gradient and Hessian of J~μ¯,𝒦 are known in closed form. In fact, J~μ¯,𝒦 has the same gradient and Hessian as Jμ¯,𝒦, except for the terms involving differentiation with respect to a2, 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 T1,T2 and the VIX smile at T1 as a function of the computational time as of April 20, 2020 (T1=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 μ*, which must be a probability measure.

We choose the lognormal prior for the reference measure μ¯ (see remark 1) and we take n1=nv=45 knots for the integration with respect to s1 and v, and n2=25 knots for the integration with respect to s2. We set s1(1)=Fμ1-1(q)v(1)=FμV-1(q) and s1(n1)=Fμ1-1(1-q)v(nV)=FμV-1(1-q) with q=10-3 for the lowest and highest values of the quadrature grids of s1 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 S1VS2 after 35 seconds in the discrete-time model calibrated with the implied Newton algorithm versus market smiles as of April 20, 2020, with T1=30 days. (a) Smile of SPXT1=30 days. (b) Smile of VIXT1=30 days. (c) Smile of SPXT2=60 days.

Risk 0224 Guyon technical fig 3

Figure 4: SPX smiles at T1 and T2 and VIX smile at T1 (market, discrete-time model computed using the implied Newton algorithm and continuous-time extension), as of August 1, 2018. (a) Smile of SPXT1=21 days. (b) Smile of VIXT1=21 days. (c) Smile of SPXT2=51 days.

Risk 0224 Guyon technical fig 4

In figure 2, we observe that IN is the fastest, and S the slowest. In figure 3, we plot the futures and smiles for S1V and S2 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 n1=nV=n2=45. With the independent prior, IN needs more steps to converge.

Continuous-time extension

We now extend the discrete-time model μ* (and the probability space) to build a continuous-time model for (St)t[0,T2]. The model will also include the VIX at T1V. That is, we build a probability  on C([0,T2],)×+ representing the distribution of ((St)t[0,T2],V). Here, V plays the role of a discrete-time stochastic volatility, representing the stochastic volatility anticipated at T1 for the [T1,T2] 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 (S1,S2), 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 T1 depend on both S1 and V.

Step 1: simulation of (St)t[0,T1]

We want (St)t[0,T1] to be a -martingale and S1 to have distribution μ1 under . To achieve this, one possible choice is to use a Markov functional model St=u(t,Wt)t[0,T1], where W is a -Brownian motion and u satisfies the heat equation:


Further, define the filtration:

 t:={σ((Ws)s[0,t])if t[0,T1)σ((Ws)s[0,t],V)if t[T1,T2]. 

In such a case, (St)t[0,T1] is an ((t),)-martingale, and the terminal condition u(T1,)=g is determined via quantiles so that u(T1,WT1) has distribution μ1, ie, for all x, we set:


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

   =𝔼[g(x+WT1-Wt)]=(gKT1-t)(x) (6)

where  is the convolution operator and where, for t>0Kt(x)=e-x2/2t/2πt is the heat kernel.

Step 2: simulation of V given (St)t[0,T1]

At this stage, we have simulated (St)t[0,T1] such that (St)t[0,T1] is an ((t),)-martingale and S1 has distribution μ1 under . Now, we simulate V given σ((Ws)s[0,T1]) under  as follows: the distribution of V given σ((Ws)s[0,T1]) under  is assumed to depend only on S1 and is taken equal to the distribution of V given S1 under μ*. Since S1 has distribution μ1 under both μ* and , this means that the distribution of (S1,V) is the same under μ* and ; in particular, V has the distribution μV under .

Step 3: simulation of (St)t[T1,T2] given

 𝑻𝟏.  In this last step, we build dynamics for (St)t[T1,T2] conditional on T1 such that (St)t[T1,T2] is an ((t),)-martingale starting from S1, and S2 has distribution μ2. We use once again a Markov functional construction. Given S1=s and V=v, we consider:


for all x and where μ2|s,v is the distribution of S2 given S1=s and V=v under μ*. Then, given T1, we define for t(T1,T2]:


where for every s,v>0:

   =𝔼[gs,v(x+WT2-Wt)]=(gs,v*KT2-t)(x) (7)

It is easy to check that (St)t[T1,T2] is an ((t),)-martingale starting from S1. Moreover, the distribution of S2 given (S1,V) is the same under μ* and . As the distribution of (S1,V) is the same under μ* and , we conclude that the distribution of (S1,V,S2) is the same under μ* and . In particular S1V and S2 have distributions μ1μV and μ2 under . We have thus built a model  on ((St)t[0,T2],V) such that (a) S1V and S2 have distributions μ1μV and μ2 under ; (b) (St)t[0,T2] is an ((t),)-martingale; and (c) V is the VIX at T1, since by construction:

Remark 3 (Extension to several VIX maturities)

Note that this approach can easily be iterated on intervals [Ti,Ti+1]. For example, after step 3, we will have generated S2μ2. Then, we just need to generate VT2μVT2 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 T1 and T2 and that of the VIX at T1 are displayed in figure 4; they were computed by Monte Carlo simulation with 105 paths.


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 T1 and T2 and the VIX future and VIX smile at T1, thus taking into account information about the forward volatility at T1 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 105 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.


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 μ* 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 (STi) (where the Ti 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.



  • 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, 
  • 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, 
  • 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 

Only users who have a paid subscription or are part of a corporate subscription are able to print or copy content.

To access these options, along with all other subscription benefits, please contact or view our subscription options here:

You are currently unable to copy this content. Please contact to find out more.

Most read articles loading...

You need to sign in to use this feature. If you don’t have a account, please register for a trial.

Sign in
You are currently on corporate access.

To use this feature you will need an individual account. If you have one already please sign in.

Sign in.

Alternatively you can request an individual account here