Table Of ContentON VIX FUTURES IN THE ROUGH BERGOMI MODEL
ANTOINEJACQUIER,CLAUDEMARTINI,ANDAITORMUGURUZA
Abstract. The rough Bergomi model introduced by Bayer, Friz and Gatheral [3] has been outperforming
7 conventionalMarkovianstochasticvolatilitymodelsbyreproducingimpliedvolatilitysmilesinaveryrealistic
1
manner,inparticularforshortmaturities. WeinvestigateherethedynamicsoftheVIXandtheforwardvariance
0
curvegeneratedbythismodel,anddevelopefficientpricingalgorithmsforVIXfuturesandoptions. Wefurther
2
analysethevalidityoftheroughBergomimodeltojointlydescribetheVIXandtheSPX,andpresentajoint
n
a calibrationalgorithmbasedonthehybridschemebyBennedsen,LundeandPakkanen[4].
J
6
1. Introduction
1
] Volatility, though not directly observed nor traded, is a fundamental object on financial markets, and has
R
beenthecentreofattentionofdecadesoftheoreticalandpracticalresearch,bothtoestimateitandtouseitfor
P
. tradingpurposes. Theformergoalhasusuallybeencarriedoutunderthehistoricalmeasure(P)whilethelatter,
n
i through the introduction of volatility derivatives (VIX and related family), has been evolving under the pricing
f
- measure Q. Most models used for pricing purposes (Heston [19], SABR [17], Bergomi [5]) are constructed
q
[ under Q and are of Markovian nature (making pricing, and hence calibration, easier). Recently, Gatheral,
1 Jaisson and Rosenbaum [14] broke this routine and introduced a fractional Brownian motion as driving factor
v
of the volatility process. This approach (Rough Fractional Stochastic Volatility, RFSV for short) opens the
0
6 door to revisiting classical pricing and calibration conundrums. They, together with the subsequent paper by
2
Bayer,FrizandGatheral,(seealso[1,12])inparticularshowedthatthesemodelswereabletocapturetheextra
4
0 steepness of the implied volatility smile in Equity markets for short maturities, which continuous Markovian
.
1 stochastic volatility models fail to describe. The icing on the cake is the (at last!!) reconciliation between the
0
7 twomeasuresPandQwithinagivenmodel,showingremarkableresultsbothforestimationandforprediction.
1
One of the key issues in Equity markets is, not only to fit the (SPX) implied volatility smile, but to do so
:
v
jointlywithacalibrationoftheVIX(Futuresandideallyoptions). Gatheral’s[15]doublemeanrevertingprocess
i
X
is the leading (Markovian) continuous model in this direction, while models with jumps have been proposed
r
a abundantly by Carr and Madan [9] and Kokholm and Stisen [20]. This issue was briefly tackled by Bayer, Friz
and Gatheral [3] for a particular rough model (rough Bergomi), and we aim here at providing a deeper analysis
of VIX dynamics under this rough model and at implementing pricing schemes for VIX Futures and options.
Our main contribution is a precise link between the forward variance curve (ξ (·)) and the initial forward
T T≥0
variancecurveξ (·)intheroughBergomimodel. Thisinturn,allowsusnotonlytoprovidesimulationmethods
0
Date:January17,2017.
2010 Mathematics Subject Classification. 91G20,91G99,91G60,91B25.
Key words and phrases. Impliedvolatility,fractionalBrownianmotion,roughBergomi,VIXFutures,VIXsmile.
TheauthorswouldliketothankChristianBayer,JimGatheral,MikkoPakkanenandMathieuRosenbaumforusefuldiscussions.
AJacknowledgesfinancialsupportfromtheEPSRCFirstGrantEP/M008436/1. Thenumericalimplementationshavebeencarried
outonthecollaborativeplatformZanadu(www.zanadu.io).
1
2 ANTOINEJACQUIER,CLAUDEMARTINI,ANDAITORMUGURUZA
for the VIX, but also to refine the log-normal approximation of [3] for VIX Futures, matching exactly the first
twomoments. Finally,wedevelopanefficientalgorithmforVIXFuturescalibration,uponwhichwebuildajoint
calibration method with the SPX. As opposed to the Cholesky approach in [3], we adapt the hybrid-scheme by
Bennedsen,LundeandPakkanen[4]withbettercomplexityO(nlogn). AssumingtheuniversalityoftheHurst
parameter H across VIX and SPX allows us to compute efficiently prices recursively with complexity O(n). In
passing, we also investigate the joint consistency of VIX and SPX in the market. The organisation of the paper
follows accordingly: we first introduce the rough Bergomi model and its main properties (Section 2), before
presenting its pricing power for VIX Futures (Section 3), and finally develop the joint calibration algorithm in
Section 4.
2. Rough volatility and the rough Bergomi model
Comte and Renault[7] were the first to propose a stochastic volatility model in which the instantaneous
volatility is driven by a fractional Brownian motion, with a Hurst index restricted to be greater than 1/2.
RecentlyGatheral,JaissonandRosenbaum[14]presentedanewapproachwithaHurstindexsmallerthan1/2,
producingextremelygoodfitstoobservedvolatilitydataunderthephysicalmeasureP. Thesemodelsformthe
so-calledRoughFractionalStochasticVolatility(RFSV)familythatisunderstoodasanaturalextensionofthe
classicalvolatilitymodelsdrivenbystandardBrownianmotion. OurworkfocusesonthepricingmeasureQand
we assume through this paper that the model presented by Gatheral, Jaisson and Rosenbaum [14] under P is a
reasonable model. Finally, and most importantly, we follow the recent paper by Bayer, Friz and Gatheral [3],
in order to extend the RFSV model to pricing schemes under the measure Q. More precisely, Bayer, Friz and
Gatheral [3] proposed the following model for the log stock price process X :=log(S):
1 (cid:112)
dX =− V dt+ V dW , X =0
(2.1) t 2 t t t 0
V =ξ (t)E(2νC V ), V >0,
t 0 H t 0
(cid:113)
with ν,ξ (·) > 0, E(·) is the Dol´eans-Dade [10] stochastic exponential and C := 2HΓ(2−H+) , where, for
0 H Γ(H+)Γ(2−2H)
notational convenience (throughout the paper), we use the symbols H :=H±1. All the processes are defined
± 2
onagivenfilteredprobabilityspace(Ω,F,(F ) ,Q)supportingthetwostandardBrownianmotionsW andZ
t t≥0
(see below). The initial forward variance curve is observed at inception, and we therefore assume without loss
of generality that it is F -measurable. The process V, defined as
0
(cid:90) t
(2.2) V := (t−u)H−dZ ,
t u
0
is a centred Gaussian process with covariance structure
(cid:90) 1(cid:18)t (cid:19)H−
E(V V )=s2H −u (1−u)H−du, for any s,t∈[0,1].
t s s
0
We shall also introduce, for any 0≤T ≤t the notations
(cid:90) t (cid:90) T
(2.3) V := (t−u)H−dZ and VT := (t−u)H−dZ .
t,T u t u
T 0
NoteinparticularthatVT =V . ThetwostandardBrownianmotionsW andZ arecorrelatedwithcorrelation
T T
parameter ρ∈(−1,1). Here, ξ (t) denotes the forward variance observed at time T for a maturity equal to t.
T
ON VIX FUTURES IN THE ROUGH BERGOMI MODEL 3
More precisely, if σ2(t) denotes the fair strike of a variance swap observed at time T and maturing at t, then
T
1 (cid:90) t
σ2(t)= ξ (u)du,
T t−T T
T
or equivalently ξ (t) = d (cid:0)(t−T)σ2(t)(cid:1). For any fixed t > 0, the process (ξ (t)) , is a martingale, i.e.
T dt T s s≤t
E[ξ (t)|F ]=ξ (t), for all u≤s≤t. Furthermore, VT is a centred Gaussian process with variance
s u u t
t2H −(t−T)2H
(2.4) V(VT)= , for t≥T,
t 2H
and covariance structure
(2.5) E(cid:0)VTVT(cid:1)=(cid:90) T [(t−u)(s−u)]H−du= (s−t)H− (cid:26)tH+F(cid:18) −t (cid:19)−(t−T)H+F(cid:18)T −t(cid:19)(cid:27),
t s H s−t s−t
0 +
for any t<s, where we introduce the function F:R →R as
−
(2.6) F(u):= F (−H ,H ,1+H ,u),
2 1 − + +
and F is the hypergeometric function [2, Chapter 15]. Finally, the quadratic variation of V is given by
2 1
t2H
(2.7) [V] = , for t≥0.
t 2H
2.1. Hybrid simulation scheme. Bayer, Friz and Gatheral [3] present a Cholesky method to simulate the
rough Bergomi model. Although exact, this method is very slow and other approaches need to be considered
for calibration purposes. Recently, Bennedsen, Lunde and Pakkanen [4] presented a new simulation scheme for
Brownian semistationary (BSS) processes. This method, as opposed to Cholesky, is an approximate method.
However, in [4] the authors show that the method yields remarkable results in the case of the rough Bergomi
model. In addition, their approach leads to a natural simulation of both the Volterra process V and the stock
price S and yields a computational complexity of order O(nlogn).
Definition 2.1. LetW beastandardBrownianmotiononagivenfilteredprobabilityspace(Ω,F,(F ) ,P).
t t≥0
(cid:82)t
AtruncatedBrowniansemistationary(BSS)processisdefinedasB(t)= g(t−s)σ(s)dW ,fort≥0,whereσ
0 s
is (F ) -predictable with locally bounded trajectories and finite second moments, and g : (0,∞) → [0,∞) is
t t≥0
Borel measurable and square integrable. We shall call it a BSS(α,W) process if furthermore
(i) there exists α∈(cid:0)−1,1(cid:1)\{0} such that g(x)=xαL (x) for all x∈(0,1], where L ∈C1((0,1]→[0,∞)),
2 2 g g
is slowly varying1 at the origin and bounded away from zero. Moreover, there exists a constant C > 0
such that |L(cid:48)(x)|≤C(1+x−1) for all x∈(0,1];
g
(ii) the function g is differentiable on (0,∞).
Under this assumption, the hybrid scheme, proposed in [4] and recalled in Appendix A, provides an efficient
way to simulate BSS processes. It applies in particular to the rough Bergomi model:
Proposition 2.2. The Volterra process V in (2.2) is a truncated BSS(H ,Z) process.
−
Proof. From(2.2),g(x)≡xH− andσ(·)≡1asinDefinition2.1,sothatV isaBSS process. SinceH− ∈(−21,0),
then L ≡1, and V satisfies Definition 2.1(i); Definition 2.1(i) trivially holds, and so does the corollary. (cid:3)
g
1AmeasurablefunctionL:(0,1]→[0,∞)isslowlyvarying[6]at0ifforanyt>0,limL(tx)/L(x)=1.
x↓0
4 ANTOINEJACQUIER,CLAUDEMARTINI,ANDAITORMUGURUZA
The corollary implies that we can apply the hybrid scheme to V. In particular, for κ = 1 the matrix form
representation of the scheme reads (recall that n :=(cid:98)nT(cid:99))
T
Z 0 ··· ··· 0
V(cid:0)n1(cid:1) Z0,1 Z ... ... 0 (cid:0)1b∗1(cid:1)H−
V(cid:0)..n2(cid:1)= Z12,,11 Z01 ... ... ... (cid:0)nn1b1∗2(cid:1)H− ,
. .
V(cid:0)nnT(cid:1) ... ... ... ... ... (cid:0)1b∗ .. (cid:1)H−
ZnT−1,1 ZnT−2 ··· Z1 Z0 n nT−1
where the coefficients {b∗} are defined in (A.1). This matrix multiplication is, by brute force, of order O(n2),
i
however using discrete convolution we may use FFT to reduce it to O(nlogn) as suggested in [4] .
3. Rough Bergomi and VIX
Bayer,FrizandGatheral[3]brieflydiscussthelackofconsistencyoftheroughBergomimodelwithobserved
VIX options data, leading to an incorrect term structure of the VIX. In this section, we investigate in detail
the dynamics of the VIX, and propose a log-normal approximation. Additionally, we investigate the viability
of the model in terms of VIX Futures and options, and compare it to the approximation in [3].
3.1. VIX Futures in the rough Bergomi model. From now on, we fix a given maturity T ≥0, and define
the VIX at time T via the continuous-time monitoring formula
(cid:32) (cid:12) (cid:33)
1 (cid:90) T+∆ (cid:12)
VIX2 :=E d(cid:104)X ,X (cid:105)ds(cid:12)F ,
T ∆ s s (cid:12) T
T (cid:12)
where ∆ is equal to 30 days. The risk-neutral formula for the VIX future V with maturity T is then given by
T
(cid:115) (cid:12) (cid:115) (cid:12)
(3.1) VT :=E(VIXT|F0)=E ∆1 (cid:90) T+∆E(d(cid:104)Xs,Xs(cid:105)|FT)ds(cid:12)(cid:12)(cid:12)(cid:12)F0=E ∆1 (cid:90) T+∆ξT(s)ds(cid:12)(cid:12)(cid:12)(cid:12)F0.
T (cid:12) T (cid:12)
Note that, when T >0, ξ (s) is a market input which is not F -measurable, and is hence difficult to interpret
T 0
only knowing F . We shall make repeated use of the following random variable defined for any t≥T, by
0
(3.2) η (t):=exp(cid:0)2νC VT(cid:1).
T H t
Proposition 3.1. The VIX dynamics are given by
VIX =(cid:40)1 (cid:90) T+∆ξ (t)η (t)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19)dt(cid:41)1/2.
T ∆ 0 T H
T
Proof. Using Fubini’s theorem and the instantaneous variance representation in (2.1), we can write
1 (cid:90) T+∆ 1 (cid:90) T+∆ (cid:104) (cid:105)
VIX2 = E(V |F )ds= E ξ (t)E(2νC V )|F dt
T ∆ s T ∆ 0 H t T
T T
= ∆1 (cid:90) T+∆E(cid:20)ξ0(t)ηT(t)exp(cid:18)2νCHVt,T − ν2CHH2t2H(cid:19)(cid:12)(cid:12)(cid:12)(cid:12)FT(cid:21)dt,
T
with V defined in (2.3). Since η (t)∈F and ξ (t)∈F , this expression simplifies to
t,T T T 0 0
VIX2T = ∆1 (cid:90) T+∆ξ0(t)ηT(t)E(cid:20)exp(cid:18)2νCHVt,T − ν2CHH2t2H(cid:19)(cid:12)(cid:12)(cid:12)(cid:12)FT(cid:21)dt.
T
ON VIX FUTURES IN THE ROUGH BERGOMI MODEL 5
The proposition follows since V is centred Gaussian, independent of F , with variance given in (2.4), and
t,T T
E(cid:0)e2νCHVt,T(cid:12)(cid:12)FT(cid:1)=E(cid:0)e2νCHVt,T(cid:1)=exp(cid:16)ν2HCH2 (t−T)2H(cid:17). (cid:3)
The main challenge for simulation is η (t). However, since the latter is independent of ξ (·), robustness of
T 0
simulationschemesfortheVIXwillnotbeaffectedbythequalitativepropertiesoftheinitialvariancecurveξ .
0
Proposition 3.2. The forward variance curve ξ in the rough Bergomi model admits the representation
T
ξ (t)=ξ (t)η (t)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19), for any t≥T.
T 0 T H
Proof. Since E(V |F )=ξ (t) by (3.1), the proposition follows from Proposition 3.1 and the equality
t T T
E(V |F )=ξ (t)η (t)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19).
t T 0 T H
(cid:3)
Bayer, Friz and Gatheral [3] did not derive such a representation for ξ , and their approach for pricing VIX
T
derivativesreliesonanapproximationwhichavoidsthecomputationsdevelopedinthissection. Proposition3.2
allows for a better understanding of the process ξ , and for an innovative approach to price VIX derivatives.
T
3.2. Upper and lower bounds for VIX Futures.
Theorem 3.3. The following bounds hold for VIX Futures:
(3.3) 1 (cid:90) T+∆(cid:112)ξ (t)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19)dt≤V ≤(cid:40)1 (cid:90) T+∆ξ (s)ds(cid:41)1/2.
∆ 0 4H T ∆ 0
T T
Proof. The conditional Jensen’s inequality gives
VT =E(VIXT|F0)=E(cid:115)∆1 (cid:90)TT+∆ξT(s)ds(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)F0≤(cid:118)(cid:117)(cid:117)(cid:116)E(cid:32)∆1 (cid:90)TT+∆ξT(s)ds(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)F0(cid:33).
Furthermore,sinceξ isF -adapted,Fubini’stheoremalongwiththemartingalepropertyofξ yieldtheupper
0 0 T
(cid:113)
bound V = E(VIX |F ) ≤ ∆−1(cid:82)T+∆ξ (s)ds. To obtain a lower bound we use the representation in
T T 0 T 0
Proposition 3.2, and Cauchy-Schwarz’s inequality, and Fubini’s theorem, so that
(cid:115) (cid:12)
VT =E(VIXT|F0)=E ∆1 (cid:90) T+∆ξ0(t)ηT(t)exp(cid:18)ν2HCH2 [(t−T)2H −t2H](cid:19)dt(cid:12)(cid:12)(cid:12)(cid:12)F0
T (cid:12)
≥E(cid:34) 1 (cid:90) T+∆(cid:112)ξ (t)η (t)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19)dt(cid:12)(cid:12)(cid:12)F (cid:35)
∆ 0 T 2H (cid:12) 0
T (cid:12)
= 1 (cid:90) T+∆(cid:112)ξ (t)E(cid:16)(cid:112)η (t)(cid:17)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19)dt
∆ 0 T 2H
T
= 1 (cid:90) T+∆(cid:112)ξ (t)exp(cid:18)ν2CH2 (cid:2)t2H −(t−T)2H(cid:3)(cid:19)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19)dt
∆ 0 4H 2H
T
= 1 (cid:90) T+∆(cid:112)ξ (t)exp(cid:18)ν2CH2 (cid:2)(t−T)2H −t2H(cid:3)(cid:19)dt,
∆ 0 4H
T
since ηT(t)1/2 is log-normal (Proposition 3.1), so that E((cid:112)ηT(t))=exp(cid:16)ν24CHH2 (cid:2)t2H −(t−T)2H(cid:3)(cid:17). (cid:3)
6 ANTOINEJACQUIER,CLAUDEMARTINI,ANDAITORMUGURUZA
We perform a numerical experiment to check the tightness of the bounds obtained in Proposition 3.3). For
this analysis we consider three qualitative scenarios for the initial forward variance curve:
√
(3.4) Scenario 1: ξ (t)=0.2342; Scenario 2: ξ (t)=0.2342(1+t)2; Scenario 3: ξ (t)=0.2342 1+t.
0 0 0
Figures 1 suggest that the lower bound given in Proposition 3.3) is surprisingly tight for very different shapes
VIX futures scenario 1 Absolute difference
0.24 0.0020
Absolute difference
0.23
0.0015
VIX futures price000...222012 MLUopowpneeterr -bbCooauurlnnodd Absolute difference0.0010
0.0005
0.19
0.108.0 0.5 1.0 1.5 2.0 0.00000.0 0.5 1.0 1.5 2.0
Maturity Maturity
VIX futures scenario 2 Absolute difference
0.8 0.0025
Monte-Carlo Absolute difference
0.7 LUopwpeerr bboouunndd 0.0020
VIX futures price000...465 Absolute difference00..00001105
0.3 0.0005
0.02.0 0.5 1.0 1.5 2.0 0.00000.0 0.5 1.0 1.5 2.0
Maturity Maturity
VIX futures scenario 3 Absolute difference
0.32 0.0012
Monte-Carlo Absolute difference
0.30 LUopwpeerr bboouunndd 0.0010
VIX futures price00..2268 Absolute difference000...000000000468
0.24 0.0002
0.202.0 0.5 1.0 1.5 2.0 0.00000.0 0.5 1.0 1.5 2.0
Maturity Maturity
Figure 1. Bounds vs. Monte Carlo (Truncated Cholesky) in all three scenarios.
ofξ . Thiscanbeexplainedwiththefollowingargument: considerasimplifiedanddeterministicversionofthe
0
VIX futures price in Proposition 3.1), denoted by
(cid:115)1 (cid:90) T+∆ (cid:114) ∆ ∆2
φ(T):= f(t)dt= f(T)+ f(cid:48)(T)+ f(cid:48)(cid:48)(T)+O(∆3),
∆ 2 6
T
for some strictly positive (deterministic) function f ∈C2(R). We further introduce
1 (cid:90) T+∆(cid:112) (cid:112) ∆ f(cid:48)(T)
ψ(T):= f(t)dt= f(T)+ +O(∆2),
(cid:112)
∆ 4 f(T)
T
which is the corresponding lower bound by Cauchy-Schwarz’s (or Jensen’s) inequality, so that
(cid:18)f(cid:48)(cid:48)(T) f(cid:48)(T)2 (cid:19)
φ(T)2−ψ(T)2 =∆2 − +O(∆3).
6 16f(T)
ON VIX FUTURES IN THE ROUGH BERGOMI MODEL 7
Hence, we observe that for small ∆, as it is the case in VIX futures, the lower bound is very close from the
original value which explains (at least for the deterministic case) the behaviour observed in Figure 1.
3.3. NumericalimplementationofVIXprocess. Inthissection,weinvestigatedifferentsimulationschemes
for the VIX in the rough Bergomi model.
3.3.1. HybridschemeandforwardEulerapproach. Inordertosimulatetheprocess(VT) itisimportant
t t∈[T,T+∆]
tonoticefrom(2.3)thatthekernelhasasingularityonlyforVT,hencewemayovercomethisbysimulatingthe
T
process using the hybrid scheme from Section 2.1. Then, we may easily extract Z and simulate (VT)
t t∈(T,T+∆]
using the forward Euler scheme with complexity O(n). A forward Euler scheme is chosen in this case due
to the fact that the kernel in (2.3) no longer has a singularity for t ∈ (T,T +∆] and this method is faster
than the hybrid scheme. Once (VT) is simulated, numerical integration routines may be used to
t t∈[T,T+∆]
simulate the VIX process using the expression in Proposition 3.1. It must be pointed out that this approach
is computationally expensive and memory consuming since it involves to simulate the Volterra process using
the hybrid scheme with complexity of O(nlogn) and additionally, each VT by forward Euler. The simulation
t
algorithm can be summarised as follows:
Algorithm 3.4 (VIX simulation in the rough Bergomi model). Fix a grid T ={t } and κ≥1.
i i=0,...,nT
(1) Simulate the Volterra process (V ) using the hybrid scheme in Appendix A, yielding a sample of
t t∈[0,T]
the random variable VT =V ;
T T
(2) extract the path of the Brownian motion Z driving the Volterra process.
Zti =Zti−1 +nH−(V(ti)−V(ti−1)), for i=1,...,κ,
Z =Z +Z , for i>κ;
ti ti−1 i−1
(3) fix a grid T = {τ } on [T,T + ∆] and approximate the continuous-time process VT by the
j j=0,...,N
discrete-time version V(cid:101)T defined via the following forward Euler scheme:
V(cid:101)T :=VT and V(cid:101)T :=(cid:88)nT Zti −Zti−1 , for j =1,...,N;
τ0 T τj (τj −ti−1)−H−
i=1
(4) compute the VIX process via numerical integration, for example using a composite trapezoidal rule:
1/2
VIX ≈1 N(cid:88)−1Q2T,τj +Q2T,τj+1(τ −τ ) ,
T ∆ 2 j j−1
j=0
where Q2T,τj :=ξ0(τj)exp(cid:16)2vCHV(cid:101)τTj(cid:17)exp(cid:16)ν2HCH2 (cid:0)(τj −T)2H −τj2H(cid:1)(cid:17).
Remark 3.5. Step 4 may obviously be replaced by any available numerical integration routine, but one must
then carefully choose the partition in Step 3.
3.3.2. Truncated Cholesky approach. Alternatively, one could use the more expensive, yet exact, Cholesky de-
composition to simulate VT on [T,T +∆] since its covariance structure is known from (2.5). However, com-
putational complexity aside, with the same grid T as in Algorithm 3.4, numerical experiments suggest that
the determinant of the covariance matrix is equal to zero when using more than n = 8 discretisation points.
T
Hence, although valid in theory, the Cholesky approach is not feasible numerically. This fact implies that there
exists strong linear dependence. In fact, for any ε > 0, the strict inequality corr(VT,VT ) < corr(VT,VT )
t1 t1+ε t2 t2+ε
holds for all T <t <t , as well as the following:
1 2
8 ANTOINEJACQUIER,CLAUDEMARTINI,ANDAITORMUGURUZA
Proposition 3.6. The limit limcorr(VT,VT )=1 holds for any t∈[T,T +∆].
t t+ε
ε↓0
Proof. This follows readily from the continuity in L(cid:32) 2 of the map t(cid:55)→VT:
t
(cid:90) T
E[(VT −VT)2]= [(t+ε−u)H− −(t−u)H−]2du
t+ε t
0
(cid:20)(T +ε)1+2H+ −ε1+2H+ +T1+2H+ 2T1+H+εH+ (cid:18) T(cid:19)(cid:21)
= − F −H ,1+H ,2+H ,− .
1+2H 1+H 2 1 + + + ε
+ +
Applying the identity [2, page 564]
(cid:18) (cid:19) (cid:18) (cid:19)
Γ(c)Γ(b−a) 1 Γ(c)Γ(a−b) 1 1
F (a,b,c;z)= (−z)−a F a,a−c+1,a−b+1; + F b,b−c+1,b−a+1;
2 1 Γ(b)Γ(c−a) 2 1 z Γ(a)Γ(c−b)(−z)b2 1 z
to the case (a,b,c)=(−H ,1+H ,2+H ) and using properties of the Gamma function along with the fact
+ + +
that F (a,0,c,z)≡1, we obtain
2 1
(cid:18) T(cid:19) 1+H (cid:18)T(cid:19)H+ (cid:18) −ε(cid:19)
F −H ,1+H ,2+H ,− = + F −H ,2+2H ,−2H , .
2 1 + + + ε 1+2H ε 2 1 + + + T
+
Finally, the series representation of F [2, Chapter 15.1.1] implies that F (cid:0)−H ,2+2H ,−2H ,−ε(cid:1) con-
2 1 2 1 + + + T
verges to 1 as ε tends to zero, and hence that E[(VT −VT)2] tends to zero. (cid:3)
t+ε t
InlightofProposition3.6,wemodelexactlythedependencestructureonthefirst8gridpointst ,...,t ,then
1 8
truncatetheCholeskydecompositionandcomputethecorrelationsρ :=corr(VT ,VT ),fori=0,...,n−9
i t8+i t8+i+1
to approximate the process by adequately rescaling and correlating each pair of subsequent grid points. In
contrasttothehybrid+forwardEulerscheme,thecomputationalcomplexityismuchlower,sincetheCholesky
method is truncated with only 8 components. The VIX simulation algorithm therefore reads as follows:
Algorithm 3.7 (VIX simulation (truncated Cholesky)). Fix a grid T={τ } on [T,T +∆],
j j=0,...,N
(i) compute the covariance matrix of (VT) using (2.5);
τj i=j,...,8
(ii) generate {VT} by correlating and rescaling using (2.5):
τj j=1,...,N
(cid:113) ρ(VT ,VT)VT (cid:113)
VτTj = V(VτTj) (cid:113)τj−V1(VτTj )τj−1 + 1−ρ(VτTj−1,VτTj)2N(0,1), for j =9,...,N;
τj−1
(iii) compute the VIX via numerical integration as in Algorithm 3.4(4).
3.3.3. Numerical experiment. We compute the price of VIX Futures using the simulation algorithm introduced
in the previous section. We set the same parameters as in [3] and [4]:
√
C 2H
(3.5) ξ =0.2352, H =0.07, ν =1.9 H ≈1.2287, κ=2.
0 2
We perform 105 simulations for the hybrid scheme + forward Euler (HS+FE) method, while 106 simulations
are used for the truncated Cholesky. Figure 2 suggests that both methods agree qualitatively and converge
to a similar output. In particular, the truncated Cholesky approach seems to suffer from larger oscillations as
maturity increases, even with 106 simulations. Nevertheless, Figure 2 indicates that the Monte-Carlo variance
increases in T for both schemes, which is confirmed in Figure 3, where the error also increases with maturity.
On the other hand, Figure 4 shows that the HS+FE methods is slower than the truncated Cholesky method,
which is consistent with the computational complexities discussed in the previous section. In particular, the
computational time of the Cholesky method is almost constant when using parallel computing. Figure 4 also
ON VIX FUTURES IN THE ROUGH BERGOMI MODEL 9
suggests that large simulations are needed to obtain precise prices. In light of this analysis, both methods
seem to approximate the required output in a decent manner. Even if the truncated Cholesky approach gives
a considerably fast output for each maturity, it is not considered for calibration, since its computational time
grows linearly in the number of maturities, making the algorithm too slow for reasonable calibration. Instead,
we will use the truncated Cholesky approach as a benchmark for the upcoming approximations.
VIX futures simulation methods
0.230
Hybrid+Forward Euler
0.225
Truncated Cholesky
0.220
VIX futures price0000....222201010055
0.195
0.190
0.1805.0 0.5 1.0 1.5 2.0
Time
Figure 2. VIX Futures using HS+FE and truncated Cholesky
Monte-Carlo standard deviations as a function of maturity
0.00075
0.00070
0.00065
Standard deviation0000....00000000000056450055 HTryubnrcida+teFdo rCwhaorlde sEkuyler
0.00040
0.00035
0.000300.0 0.5 1.0 1.5 2.0
Time
Figure 3. Monte-Carlo standard deviations, with 105 simulations.
Monte-Carlo standard deviations (10log-scale) Computational time
0.025 250
Hybrid+forward Euler Hybrid+forward Euler
Truncated Cholesky Truncated Cholesky
0.020 200
Standard deviation00..001105 Seconds110500
0.005 50
0.0020.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 20.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0
Simulations Simulations
Figure 4. Monte-Carlo standard deviations and computational times for a fixed maturity
(T =2 years) for both methods using efficient parallel computing.
10 ANTOINEJACQUIER,CLAUDEMARTINI,ANDAITORMUGURUZA
3.4. The VIX process and log-normal approximations. We now investigate approximate methods to
price VIX futures and options. We define the F -measurable random variable ∆VIX2 =(cid:82)T+∆ξ (t)dt. In [3],
T T T T
Bayer, Friz and Gatheral assumed that log(∆VIX2) follows a Gaussian distribution, and computed directly its
T
first and second moments, which hence fully characterise the distribution of ∆VIX2. However, since the sum
T
of log-normal random variables is not (in general) log-normal, the fact that ξ is log-normal does not imply
T
that ∆VIX2 is. In a different context (geometric Brownian motion and Asian options), Dufresne [11] proved
T
that, under certain conditions, an integral of log-normal variables asymptotically converges to a log-normal.
This approximation has been widely used for many applications [11], and Dufresne’s result motivates Bayer-
Friz-Gatheral’sassumption. Weprovidehereexactformulaeforthemeanandvarianceofthisdistribution,and
compare them numerically to those by Bayer-Friz-Gatheral.
Proposition 3.8. The following holds:
σ2 :=V(log(∆VIX2))=−2logE(∆VIX2)+logE(cid:0)(∆VIX2)2(cid:1),
T T T
σ2
µ:=E(log(∆VIX2))=logE(∆VIX2)− .
T T 2
Furthermore, E(∆VIX2)=(cid:82)T+∆ξ (t)dt and
T T 0
(3.6) E(cid:2)(∆VIX2)2(cid:3)=(cid:90) ξ (u)ξ (t)exp(cid:26)ν2CH2 (cid:2)(u−T)2H +(t−T)2H −u2H −t2H(cid:3)(cid:27)eΘu,tdudt.
T 0 0 H
[T,T+∆]2
where Θ is equal to zero if u=t and otherwise equal to Θ , where
u,t u∨t,u∧t
(cid:26)u2H −(u−T)2H +t2H −(t−T)2H (u−t)H− (cid:20) (cid:18) −t (cid:19) (cid:18)T −t(cid:19)(cid:21)(cid:27)
Θ :=2ν2C2 +2 tH+F −(t−T)H+F .
u,t H 2H H u−t u−t
+
Remark 3.9. Since all the integrals in the proposition are computed over compact intervals, they are finite as
long as ξ is well behaved on [T,T +∆]. Assuming this is indeed the case is not restrictive in practice as ξ
0 0
represents the initial forward variance curve; note in particular that continuity of ξ is sufficient.
0
Remark 3.10. The reason why Θ is defined that way is for numerical purposes. Indeed, when u < t,
u,t
Θ is not well defined since F(x) only makes sense when x ≤ 1 (details on the radius of convergence of
u,t
hypergeometric functions can be found in [2][Page 556]), even though the integral representation (3.6) is still
welldefined. However,mostnumericalpackagesimplementhypergeometricfunctionsviaseriesexpansions. The
trick from Θ to Θ allows us to bypass this issue.
u,t u,t
Proof of Proposition 3.8. The expectation follows directly from the tower property and Fubini’s theorem. For
the second moment, we use the decomposition
(cid:32) (cid:33)
(cid:90) T+∆(cid:90) T+∆ (cid:90) T+∆(cid:90) T+∆
E(cid:0)(∆VIX2)2(cid:1)=E ξ (u)ξ (t)dtdu = E(ξ (u)ξ (t))dtdu
T T T T T
T T T T
where in the last step we use that ξ (s) is F -measurable in order to apply Fubini. Using the representation
T T
obtained in Proposition 3.2, we get
(3.7) E(ξ (u)ξ (t))=ξ (u)ξ (t)exp(cid:26)ν2CH2 (cid:2)(u−T)2H −u2H +(t−T)2H −t2H(cid:3)(cid:27)E(cid:0)eϑu,t(cid:1),
T T 0 0 H
where the random variable
(cid:90) T
ϑ :=2νC (cid:2)(u−s)H− +(t−s)H−(cid:3)dZ
u,t H s
0
