This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Functional quantization of rough volatility and applications to volatility derivatives

Ofelia Bonesini Ofelia Bonesini, Department of Mathematics, University of Padova [email protected] Giorgia Callegaro Giorgia Callegaro, Department of Mathematics, University of Padova, Via Trieste 63, 35121 Padova, Italy. [email protected]  and  Antoine Jacquier Antoine Jacquier, Department of Mathematics, Imperial College London, London SW7 1NE, UK. [email protected]
Abstract.

We develop a product functional quantization of rough volatility. Since the optimal quantizers can be computed offline, this new technique, built on the insightful works by Luschgy and Pagès [31, 32, 35], becomes a strong competitor in the new arena of numerical tools for rough volatility. We concentrate our numerical analysis on the pricing of options on the VIX and realized variance in the rough Bergomi model [4] and compare our results to other benchmarks recently suggested.

Keywords: Riemann-Liouville process, Volterra process, functional quantization, series expansion, rough volatility, VIX options.

1. Introduction

Gatheral, Jaisson and Rosenbaum [18] recently introduced a new framework for financial modelling. To be precise — according to the reference website https://sites.google.com/site/roughvol/home — almost twenty-four hundred days have passed since instantaneous volatility was shown to have a rough nature, in the sense that its sample paths are α\alpha-Hölder-continuous with α<12\alpha<\frac{1}{2}. Many studies, both empirical [8, 15, 16] and theoretical [14, 3], have confirmed this, showing that these so-called rough volatility models are a more accurate fit to the implied volatility surface and to estimate historical volatility time series.

On equity markets, the quality of a model is usually measured by its ability to calibrate not only to the SPX implied volatility but also VIX Futures and the VIX implied volatility. The market standard models had so far been Markovian, in particular the double mean-reverting process [19, 24], Bergomi’s model [9] and, to some extent, jump models [10, 29]. However, they each suffer from several drawbacks, which the new generation of rough volatility models seems to overcome. For VIX Futures pricing, the rough version of Bergomi’s model was thoroughly investigated in [26], showing accurate results. Nothing comes for free though and the new challenges set by rough volatility models lie on the numerical side, as new tools are needed to develop fast and accurate numerical techniques. Since classical simulation tools for fractional Brownian motions are too slow for realistic purposes, new schemes have been proposed to speed it up, among which the Monte Carlo hybrid scheme [8, 33], a tree formulation [22], quasi Monte-Carlo methods [7] and Markovian approximations [1, 11].

We suggest here a new approach, based on product functional quantization [35]. Quantization was originally conceived as a discretization technique to approximate a continuous signal by a discrete one [38], later developed at Bell Laboratory in the 1950s for signal transmission [20]. It was however only in the 1990s that its power to compute (conditional) expectations of functionals of random variables [21] was fully understood. Given an d\mathbb{R}^{d}-valued random vector on some probability space, optimal vector quantization investigates how to select an d\mathbb{R}^{d}-valued random vector X^\widehat{X}, supported on at most NN elements, that best approximates XX according to a given criterion (such as the LrL^{r}-distance, r1r\geq 1). Functional quantization is the infinite-dimensional version, approximating a stochastic process with a random vector taking a finite number of values in the space of trajectories for the original process. It has been investigated precisely [31, 35] in the case of Brownian diffusions, in particular for financial applications [36]. However, optimal functional quantizers are in general hard to compute numerically and instead product functional quantizers provide a rate-optimal (so, in principle, sub-optimal) alternative often admitting closed-form expressions [32, 36].

In Section 2 we briefly review important properties of Gaussian Volterra processes, displaying a series expansion representation, and paying special attention to the Riemann-Liouville case in Section 2.2. This expansion yields, in Section 3, a product functional quantization of the processes, that shows an L2L^{2}-error of order log(N)H\log(N)^{-H}, with NN the number of paths and HH a regularity index. We then show, in Section 3.1, that these functional quantizers, although sub-optimal, are stationary. We specialise our setup to the generalized rough Bergomi model in Section 4 and show how product functional quantization applies to the pricing of VIX Futures and VIX options, proving in particular precise rates of convergence. Finally, Section 5 provides a numerical confirmation of the quality of our approximations for VIX Futures and Call Options on the VIX in the rough Bergomi model, benchmarked against other existing schemes. In this Section, we also discuss how product functional quantization of the Riemann-Liouville process itself can be exploited to price options on realized variance.

Notations.

We set \mathbb{N} as the set of strictly positive natural numbers. We denote by 𝒞[0,1]\mathcal{C}[0,1] the space of real-valued continuous functions over [0,1][0,1] and by L2[0,1]L^{2}[0,1] the Hilbert space of real-valued square integrable functions on [0,1][0,1], with inner product f,gL2[0,1]:=01f(t)g(t)𝑑t\langle f,g\rangle_{L^{2}[0,1]}:=\int_{0}^{1}f(t)g(t)dt, inducing the norm fL2[0,1]:=(01|f(t)|2𝑑t)1/2\|f\|_{L^{2}[0,1]}:=(\int_{0}^{1}|f(t)|^{2}dt)^{1/2}, for each f,gL2[0,1]f,g\in L^{2}[0,1]. L2()L^{2}(\mathbb{P}) denotes the space of square integrable (with respect to \mathbb{P}) random variables.

2. Gaussian Volterra processes on +\mathbb{R}_{+}

For clarity, we restrict ourselves to the time interval [0,1][0,1]. Let {Wt}t[0,1]\{W_{t}\}_{t\in[0,1]} be a standard Brownian motion on a filtered probability space (Ω,,{t}t[0,1],)(\Omega,{\mathcal{F}},\{{\mathcal{F}}_{t}\}_{t\in[0,1]},\mathbb{P}), with {t}t[0,1]\{{\mathcal{F}}_{t}\}_{t\in[0,1]} its natural filtration. On this probability space we introduce the Volterra process

Zt:=0tK(ts)𝑑Ws,t[0,1],Z_{t}:=\int_{0}^{t}K(t-s)dW_{s},\qquad t\in[0,1], (1)

and we consider the following assumptions for the kernel KK:

Assumption 2.1.

There exist α(12,12){0}\alpha\in\left(-\frac{1}{2},\frac{1}{2}\right)\setminus\{0\} and L:(0,1](0,)L:(0,1]\to(0,\infty) continuously differentiable, slowly varying at 0, that is, for any t>0t>0, limx0L(tx)L(x)=1\lim_{x\downarrow 0}\frac{L(tx)}{L(x)}=1, and bounded away from 0 function with |L(x)|C(1+x1)|L^{\prime}(x)|\leq C(1+x^{-1}), for x(0,1]x\in(0,1], for some C>0C>0, such that

K(x)=xαL(x),x(0,1].K(x)=x^{\alpha}L(x),\qquad x\in(0,1].

This implies in particular that KL2[0,1]K\in L^{2}[0,1], so that the stochastic integral (1) is well defined. The Gamma kernel, with K(u)=eβuuαK(u)=e^{-\beta u}u^{\alpha}, for β>0\beta>0 and α(12,12){0}\alpha\in(-\frac{1}{2},\frac{1}{2})\setminus\{0\}, is a classical example satisfying Assumption 2.1. Straightforward computations show that the covariance function of ZZ reads

RZ(s,t)=0tsK(tu)K(su)𝑑u,s,t[0,1].R_{Z}(s,t)=\int_{0}^{t\wedge s}K(t-u)K(s-u)du,\quad s,t\in[0,1]. (2)

Under Assumption 2.1, ZZ is a Gaussian process admitting a version which is ε\varepsilon-Hölder continuous for any ε<12+α=H\varepsilon<\frac{1}{2}+\alpha=H and hence also admits a continuous version [8, Proposition 2.11].

2.1. Series expansion

We introduce a series expansion representation for the centered Gaussian process ZZ in (1), which will be key to develop its functional quantization. Inspired by [32], introduce the stochastic process

Yt:=n1𝒦[ψn](t)ξn,t[0,1],Y_{t}{:=}\sum_{n\geq 1}\mathcal{K}[\psi_{n}](t)\xi_{n},\qquad t\in[0,1], (3)

where {ξn}n1\{\xi_{n}\}_{n\geq 1} is a sequence of i.i.d. standard Gaussian random variables, {ψn}n1\{\psi_{n}\}_{n\geq 1} denotes the orthonormal basis of L2[0,1]L^{2}[0,1]:

ψn(t)=2cos(tλn), with λn=4(2n1)2π2,\psi_{n}(t)=\sqrt{2}\cos\left(\frac{t}{\sqrt{\lambda_{n}}}\right),\quad\text{ with }\lambda_{n}=\frac{4}{(2n-1)^{2}\pi^{2}}, (4)

and the operator 𝒦:L2[0,1]𝒞[0,1]\mathcal{K}:{L}^{2}[0,1]\to\mathcal{C}[0,1] is defined for fL2[0,1]f\in{L}^{2}[0,1] as

𝒦[f](t):=0tK(ts)f(s)𝑑s,for all t[0,1].\mathcal{K}[f](t):=\int_{0}^{t}K(t-s)f(s)ds,\qquad\text{for all }t\in[0,1]. (5)
Remark 2.2.

The stochastic process YY in (3) is defined as a weighted sum of independent centered Gaussian variables, so for every t[0,1]t\in[0,1] the random variable YtY_{t} is a centered Gaussian random variable and the whole process YY is Gaussian with zero mean.

We set the following assumptions on the functions {𝒦[ψn]}n\{\mathcal{K}[\psi_{n}]\}_{n\in\mathbb{N}}:

Assumption 2.3.

There exists H(0,12)H\in(0,\frac{1}{2}) such that

  • (A)

    there is a constant C1>0C_{1}>0 for which, for any n1n\geq 1, 𝒦[ψn]\mathcal{K}[\psi_{n}] is (H+12)(H+\frac{1}{2})-Hölder continuous, with

    sups,t[0,1],st|𝒦[ψn](t)𝒦[ψn](s)||ts|H+12C1n;\sup_{s,t\in[0,1],s\neq t}\frac{|\mathcal{K}[\psi_{n}](t)-\mathcal{K}[\psi_{n}](s)|}{|t-s|^{H+\frac{1}{2}}}\leq C_{1}n;
  • (B)

    there exists a constant C2>0C_{2}>0 such that

    supt[0,1]|𝒦[ψn](t)|C2n(H+12), for all n1.\sup_{t\in[0,1]}|\mathcal{K}[\psi_{n}](t)|\leq C_{2}n^{-(H+\frac{1}{2})},\quad\text{ for all }n\geq 1.

Notice that under these assumptions, the series (3) converges both almost surely and in L2()L^{2}(\mathbb{P}) for each t[0,1]t\in[0,1] by Khintchine-Kolmogorov Convergence Theorem [12, Theorem 1, Section 5.1].

It is natural to wonder whether Assumption 2.1 implies Assumption 2.3 given the basis functions (4). This is far from trivial in our general setup and we provide examples and justifications later on for models of interest. Similar considerations with slightly different conditions can be found in [32]. We now focus on the variance-covariance structure of the Gaussian process YY.

Lemma 2.4.

For any s,t[0,1]s,t\in[0,1], the covariance function of YY is given by

RY(s,t):=𝔼[YsYt]=0tsK(tu)K(su)𝑑u.R_{Y}(s,t):=\mathbb{E}[Y_{s}Y_{t}]=\int_{0}^{t\wedge s}K(t-u)K(s-u)du.
Proof.

Exploiting the definition of YY in (3), the definition of 𝒦\mathcal{K} in (5) and the fact that the random variable ξn\xi_{n}’s are i.i.d. standard Normal, we obtain

RY(s,t)\displaystyle R_{Y}(s,t) =\displaystyle= 𝔼[YsYt]=𝔼[(n1𝒦[ψn](s)ξn)(m1𝒦[ψm](t)ξm)]=n1𝒦[ψn](s)𝒦[ψn](t)\displaystyle\mathbb{E}[Y_{s}Y_{t}]=\mathbb{E}\Bigg{[}\Big{(}\sum_{n\geq 1}\mathcal{K}[\psi_{n}](s)\xi_{n}\Big{)}\Big{(}\sum_{m\geq 1}\mathcal{K}[\psi_{m}](t)\xi_{m}\Big{)}\Bigg{]}=\sum_{n\geq 1}\mathcal{K}[\psi_{n}](s)\mathcal{K}[\psi_{n}](t)
=\displaystyle= n1(01K(su)𝟙[0,s](u)ψn(u)𝑑u01K(tr)𝟙[0,t](r)ψn(r)𝑑r)\displaystyle\sum_{n\geq 1}\Big{(}\int_{0}^{1}K(s-u)\mathbb{1}_{[0,s]}(u)\psi_{n}(u)du\int_{0}^{1}K(t-r)\mathbb{1}_{[0,t]}(r)\psi_{n}(r)dr\Big{)}
=\displaystyle= n1K(s)𝟙[0,s](),ψnL2[0,1]K(t)𝟙[0,t](),ψnL2[0,1]\displaystyle\sum_{n\geq 1}\langle K(s-\cdot)\mathbb{1}_{[0,s]}(\cdot),\psi_{n}\rangle_{L^{2}[0,1]}\cdot\langle K(t-\cdot)\mathbb{1}_{[0,t]}(\cdot),\psi_{n}\rangle_{L^{2}[0,1]}
=\displaystyle= n1K(t)𝟙[0,t](),K(s)𝟙[0,s](),ψnL2[0,1]ψnL2[0,1]\displaystyle\sum_{n\geq 1}\Big{\langle}K(t-\cdot)\mathbb{1}_{[0,t]}(\cdot),\langle K(s-\cdot)\mathbb{1}_{[0,s]}(\cdot),\psi_{n}\rangle_{L^{2}[0,1]}\psi_{n}\Big{\rangle}_{L^{2}[0,1]}
=\displaystyle= K(t)𝟙[0,t](),n1K(s)𝟙[0,s](),ψnL2[0,1]ψnL2[0,1]\displaystyle\Big{\langle}K(t-\cdot)\mathbb{1}_{[0,t]}(\cdot),\sum_{n\geq 1}\langle K(s-\cdot)\mathbb{1}_{[0,s]}(\cdot),\psi_{n}\rangle_{L^{2}[0,1]}\psi_{n}\Big{\rangle}_{L^{2}[0,1]}
=\displaystyle= K(t)𝟙[0,t](),K(s)𝟙[0,s]()L2[0,1]\displaystyle\langle K(t-\cdot)\mathbb{1}_{[0,t]}(\cdot),K(s-\cdot)\mathbb{1}_{[0,s]}(\cdot)\rangle_{L^{2}[0,1]}
=\displaystyle= 01K(su)𝟙[0,s](u)K(tu)𝟙[0,t](u)𝑑u=0tsK(tu)K(su)𝑑u.\displaystyle\int_{0}^{1}K(s-u)\mathbb{1}_{[0,s]}(u)K(t-u)\mathbb{1}_{[0,t]}(u)du=\int_{0}^{t\wedge s}K(t-u)K(s-u)du.

Remark 2.5.

Notice that the centered Gaussian stochastic process YY admits a continuous version, too. Indeed, we have shown that YY has the same mean and covariance function as ZZ and, consequently, that the increments of the two processes share the same distribution. Thus, [8, Proposition 2.11] applies to YY as well, yielding that the process admits a continuous version. This last key property of YY can be alternatively proved directly as done in Appendix A.2.

Lemma 2.4 implies that 𝔼[YsYt]=𝔼[ZsZt],\mathbb{E}[Y_{s}Y_{t}]=\mathbb{E}[Z_{s}Z_{t}], for all s,t[0,1]s,t\in[0,1]. Both ZZ and YY are continuous, centered, Gaussian with the same covariance structure, so from now on we will work with YY, using

Z=n1𝒦[ψn]ξn,-a.s.Z=\sum_{n\geq 1}\mathcal{K}[\psi_{n}]\xi_{n},\quad\mathbb{P}\text{-a.s.} (6)

2.2. The Riemann - Liouville case

For K(u)=uH12K(u)=u^{H-\frac{1}{2}}, with H(0,12)H\in(0,\frac{1}{2}), the process (1) takes the form

ZtH:=0t(ts)H12dWs,t[0,1],Z^{H}_{t}:=\int_{0}^{t}(t-s)^{H-\frac{1}{2}}dW_{s},\qquad t\in[0,1], (7)

where we add the superscript HH to emphasise its importance. It is called a Riemann-Liouville process (henceforth RL) (also known as Type II fractional Brownian motion or Lévy fractional Brownian motion), as it is obtained by applying the Riemann-Liouville fractional operator to the standard Brownian motion, and is an example of a Volterra process. This process enjoys properties similar to those of the fractional Brownian motion (fBM), in particular being HH-self-similar and centered Gaussian. However, contrary to the fractional Brownian motion, its increments are not stationary. For a more detailed comparison between the fBM and ZHZ^{H} we refer to [37, Theorem 5.1]. In the RL case, the covariance function RZH(,)R_{Z^{H}}(\cdot,\cdot) is available [25, Proposition 2.1] explicitly as

RZH(s,t)=1H+12(st)H+12(st)H12F12(1,12H;2H+1;stst),s,t[0,1],R_{Z^{H}}(s,t)=\frac{1}{H+\frac{1}{2}}(s\land t)^{H+\frac{1}{2}}(s\lor t)^{H-\frac{1}{2}}\ {}_{2}F_{1}\left(1,\frac{1}{2}-H;2H+1;\frac{s\land t}{s\vee t}\right),\quad s,t\in[0,1],

where F12(a,b;c;z){}_{2}F_{1}(a,b;c;z) denotes the Gauss hypergeometric function  [34, Chapter 5, Section 9]. More generally, [34, Chapter 5, Section 11], the generalized Hypergeometric functions Fqp(z){}_{p}F_{q}(z) are defined as

Fqp(z)=Fqp(a1,a2,,ap;c1,c2,,cq;z):=k=0(a1)k(a2)k(ap)k(c1)k(c2)k(cq)kzk!,{}_{p}F_{q}(z)={}_{p}F_{q}(a_{1},a_{2},\dots,a_{p};c_{1},c_{2},\dots,c_{q};z):=\sum_{k=0}^{\infty}\frac{{(a_{1})}_{k}{(a_{2})}_{k}\cdots{(a_{p})}_{k}}{{(c_{1})}_{k}{(c_{2})}_{k}\cdots{(c_{q})}_{k}}\frac{z}{k!}, (8)

with the Pochammer’s notation (a)0:=1{(a)}_{0}:=1 and (a)k:=a(a+1)(a+2)(a+k1){(a)}_{k}:=a(a+1)(a+2)\cdots(a+k-1), for k1k\geq 1, where none of the ckc_{k} are negative integers or zero. For pqp\leq q the series (8) converges for all zz and when p=q+1p=q+1 convergence holds for |z|<1|z|<1 and the function is defined outside this disk by analytic continuation. Finally, when p>q+1p>q+1 the series diverges for nonzero zz unless one of the aka_{k}’s is zero or a negative integer.

Regarding the series representation (3), we have, for t[0,1]t\in[0,1] and n1n\geq 1,

𝒦H[ψn](t):\displaystyle\mathcal{K}_{H}[\psi_{n}](t): =20t(ts)H12cos(sλn)𝑑s\displaystyle=\sqrt{2}\int_{0}^{t}(t-s)^{H-\frac{1}{2}}\cos\Big{(}\frac{s}{\sqrt{\lambda_{n}}}\Big{)}ds (9)
=221+2HtH+12F21(1;34+H2,54+H2;t24λn).\displaystyle=\frac{2\sqrt{2}}{1+2H}\ t^{H+\frac{1}{2}}\ {}_{1}F_{2}\left(1;\frac{3}{4}+\frac{H}{2},\frac{5}{4}+\frac{H}{2};-\frac{t^{2}}{4\lambda_{n}}\right). (10)

Assumption 2.3 holds in the RL case here using [32, Lemma 4] (identifying 𝒦H[ψn]\mathcal{K}_{H}[\psi_{n}] to fnf_{n} from [32, Equation (3.7)]). Assumption 2.3 (B) implies that, for all t[0,1]t\in[0,1],

n1𝒦H[ψn](t)2n1(supt[0,1]|𝒦H[ψn](t)|)2C22n11n1+2H<,\sum_{n\geq 1}\mathcal{K}_{H}[\psi_{n}](t)^{2}\leq\sum_{n\geq 1}\left(\sup_{t\in[0,1]}|\mathcal{K}_{H}[\psi_{n}](t)|\right)^{2}\leq C_{2}^{2}\sum_{n\geq 1}\frac{1}{n^{1+2H}}<\infty,

and therefore the series (3) converges both almost surely and in L2()L^{2}(\mathbb{P}) for each t[0,1]t\in[0,1] by Khintchine-Kolmogorov Convergence Theorem [12, Theorem 1, Section 5.1].

Remark 2.6.

The expansion (3) is in general not a Karhunen-Loève decomposition [36, Section 4.1.1]. In the RL case, it can be numerically checked that the basis {𝒦H[ψn]}n\{\mathcal{K}_{H}[\psi_{n}]\}_{n\in\mathbb{N}} is not orthogonal in L2[0,1]L^{2}[0,1] and does not correspond to eigenvectors for the covariance operator of the Riemann-Liouville process. In his PhD Thesis [13], Corlay exploited a numerical method to obtain approximations of the first terms in the K-L expansion of processes for which an explicit form is not available.

3. Functional quantization and error estimation

Optimal (quadratic) vector quantization was conceived to approximate a square integrable random vector X:(Ω,,)dX:(\Omega,\mathcal{F},\mathbb{P})\rightarrow\mathbb{R}^{d} by another one X^\widehat{X}, taking at most a finite number NN of values, on a grid ΓN:={x1N,x2N,,xNN}\Gamma^{N}:=\{x_{1}^{N},x_{2}^{N},\dots,x_{N}^{N}\}, with xiNd,i=1,,Nx_{i}^{N}\in\mathbb{R}^{d},i=1,\dots,N. The quantization of XX is defined as X^:=ProjΓN(X)\widehat{X}:=\mathrm{Proj}_{\Gamma^{N}}(X), where ProjΓN:dΓN\mathrm{Proj}_{\Gamma^{N}}:\mathbb{R}^{d}\rightarrow\Gamma^{N} denotes the nearest neighbour projection. Of course the choice of the NN-quantizer ΓN\Gamma^{N} is based on a given optimality criterion: in most cases ΓN\Gamma^{N} minimizes the distance 𝔼[|XX^|2]1/2\mathbb{E}[|X-\widehat{X}|^{2}]^{1/2}. We recall basic results for one-dimensional standard Gaussian, which shall be needed later, and refer to [21] for a comprehensive introduction to quantization.

Definition 3.1.

Let ξ\xi be a one-dimensional standard Gaussian on a probability space (Ω,,)(\Omega,{\mathcal{F}},\mathbb{P}). For each nn\in\mathbb{N}, we define the optimal quadratic nn-quantization of ξ\xi as the random variable ξ^n:=ProjΓn(ξ)=i=1nxin1Ci(Γn)(ξ)\widehat{\xi}^{n}:=\mathrm{Proj}_{\Gamma^{n}}(\xi)=\sum_{i=1}^{n}x_{i}^{n}1_{C_{i}(\Gamma^{n})}(\xi), where Γn={x1n,,xnn}\Gamma^{n}=\{x_{1}^{n},\dots,x_{n}^{n}\} is the unique optimal quadratic nn-quantizer of ξ\xi, namely the unique solution to the minimization problem

minΓn,Card(Γn)=n𝔼[|ξProjΓn(ξ)|2],\min_{\Gamma^{n}\subset\mathbb{R},\mathrm{Card}(\Gamma^{n})=n}\mathbb{E}[|\xi-\mathrm{Proj}_{\Gamma^{n}}(\xi)|^{2}],

and {Ci(Γn)}i{1,,n}\{C_{i}(\Gamma^{n})\}_{i\in\{1,\dots,n\}} is a Voronoi partition of \mathbb{R}, that is a Borel partition of \mathbb{R} that satisfies

Ci(Γn){y:|yxin|=min1jn|yxjn|}C¯i(Γn),C_{i}(\Gamma^{n})\subset\left\{y\in\mathbb{R}:|y-x_{i}^{n}|=\min_{1\leq j\leq n}|y-x_{j}^{n}|\right\}\subset\overline{C}_{i}(\Gamma^{n}),

where the right-hand side denotes the closure of the set in \mathbb{R}.

The unique optimal quadratic nn-quantizer Γn={x1n,,xnn}\Gamma^{n}=\{x_{1}^{n},\dots,x_{n}^{n}\} and the corresponding quadratic error are available online, at http://www.quantize.maths-fi.com/gaussian_database for n{1,,5999}n\in\{1,\dots,5999\}.

Given a stochastic process, viewed as a random vector taking values in its trajectories space, such as L2[0,1]L^{2}[0,1], functional quantization does the analogue to vector quantization in an infinite-dimensional setting, approximating the process with a finite number of trajectories. In this section, we focus on product functional quantization of the centered Gaussian process ZZ from (1) of order NN (see [35, Section 7.4] for a general introduction to product functional quantization). Recall that we are working with the continuous version of ZZ in the series (6). For any m,Nm,N\in\mathbb{N}, we introduce the following set, which will be of key importance all throughout the paper:

𝒟mN:={𝐝m:i=1md(i)N}.{\mathcal{D}}_{m}^{N}:=\left\{\boldsymbol{\mathrm{d}}\in\mathbb{N}^{m}:\prod_{i=1}^{m}d(i)\leq N\right\}. (11)
Definition 3.2.

A product functional quantization of ZZ of order NN is defined as

Z^t𝐝:=n=1m𝒦[ψn](t)ξ^nd(n),t[0,1],\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t}:=\sum_{n=1}^{m}\mathcal{K}[\psi_{n}](t)\widehat{\xi}_{n}^{d(n)},\qquad t\in[0,1], (12)

where 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N}, for some mm\in\mathbb{N}, and for every n{1,,m},n\in\{1,\dots,m\}, ξ^nd(n)\widehat{\xi}_{n}^{d(n)} is the (unique) optimal quadratic quantization of the standard Gaussian random variable ξn\xi_{n} of order d(n)d(n), according to Definition 3.1.

Remark 3.3.

The condition i=1md(i)N\prod_{i=1}^{m}d(i)\leq N in Equation (11) motivates the wording ‘product’ functional quantization. Clearly, the optimality of the quantizer also depends on the choice of mm and 𝐝\boldsymbol{\mathrm{d}}, for which we refer to Proposition 3.6 and Section  5.1.

Before proceeding, we need to make precise the explicit form for the product functional quantizer of the stochastic process ZZ:

Definition 3.4.

The product functional 𝐝\boldsymbol{\mathrm{d}}-quantizer of ZZ is defined as

χi¯𝐝(t):=n=1m𝒦[ψn](t)xind(n),t[0,1],i¯=(i1,,im),\chi_{\underline{i}}^{\boldsymbol{\mathrm{d}}}(t):=\sum_{n=1}^{m}\mathcal{K}[\psi_{n}](t)\ x_{i_{n}}^{d(n)},\qquad t\in[0,1],\quad\underline{i}=(i_{1},\dots,i_{m}), (13)

for 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N} and 1ind(n)1\leq i_{n}\leq d(n) for each n=1,,m.n=1,\dots,m.

Remark 3.5.

Intuitively, the quantizer is chosen as a Cartesian product of grids of the one-dimensional standard Gaussian random variables. So, we also immediately find the probability associated to every trajectory χi¯𝐝\chi_{\underline{i}}^{\boldsymbol{\mathrm{d}}}: for every i¯=(i1,,im)n=1m{1,,d(n)}\underline{i}=(i_{1},\dots,i_{m})\in\prod_{n=1}^{m}\{1,\dots,d(n)\},

(Z^𝐝=χi¯𝐝)=n=1m(ξnCin(Γd(n))),\mathbb{P}(\widehat{Z}^{\boldsymbol{\mathrm{d}}}=\chi_{\underline{i}}^{\boldsymbol{\mathrm{d}}})=\prod_{n=1}^{m}\mathbb{P}(\xi_{n}\in C_{i_{n}}(\Gamma^{d(n)})), (14)

where Cj(Γd(n))C_{j}(\Gamma^{d(n)}) is the jj-th Voronoi cell relative to the d(n)d(n)-quantizer Γd(n)\Gamma^{d(n)} in Definition 3.1.

The following, proved in Appendix A.1, deals with the quantization error estimation and its minimization and provides hints to choose (m,𝐝)(m,\boldsymbol{\mathrm{d}}). A similar result on the error can be obtained applying [32, Theorem 2] to the first example provided in the reference. For completeness we preferred to prove the result in an autonomous way in order to further characterize the explicit expression of the rate optimal parameters. Indeed, we then compare these rate optimal parameters with the (numerically computed) optimal ones in Section 5.1. The symbol \lfloor\cdot\rfloor denotes the lower integer part.

Proposition 3.6.

Under Assumption 2.3, for any N1N\geq 1, there exist m(N)m^{*}(N)\in\mathbb{N} and C>0C>0 such that

𝔼[Z^𝐝NZL2[0,1]2]12Clog(N)H,\mathbb{E}\left[\left\|\widehat{Z}^{\boldsymbol{\mathrm{d}}^{*}_{N}}-Z\right\|^{2}_{L^{2}[0,1]}\right]^{\frac{1}{2}}\leq C\log(N)^{-H},

where 𝐝N𝒟m(N)N\boldsymbol{\mathrm{d}}^{*}_{N}\in{\mathcal{D}}_{m^{*}(N)}^{N} and with, for each n=1,,m(N)n=1,\dots,m^{*}(N),

dN(n)=N1m(N)n(H+12)(m(N)!)2H+12m(N).d^{*}_{N}(n)=\Big{\lfloor}N^{\frac{1}{m^{*}(N)}}n^{-(H+\frac{1}{2})}\left(m^{*}(N)!\right)^{\frac{2H+1}{2m^{*}(N)}}\Big{\rfloor}. (15)

Furthermore m(N)=𝒪(log(N))m^{*}(N)=\mathcal{O}(\log(N)).

Remark 3.7.

In the RL case, the trajectories of Z^H,𝐝\widehat{Z}^{H,\boldsymbol{\mathrm{d}}} are easily computable and they are used in the numerical implementations to approximate the process ZHZ^{H}. In practice, the parameters mm and 𝐝=(d(1),,d(m))\boldsymbol{\mathrm{d}}=(d(1),\dots,d(m)) are chosen as explained in Section 5.1.

3.1. Stationarity

We now show that the quantizers we are using are stationary. The use of stationary quantizers is motivated by the fact that their expectation provides a lower bound for the expectation of convex functionals of the process (Remark 3.9) and they yield a lower (weak) error in cubature formulae [35, page 26]. We first recall the definition of stationarity for the quadratic quantizer of a random vector [35, Definition 1].

Definition 3.8.

Let XX be an d\mathbb{R}^{d}-valued random vector on (Ω,,)(\Omega,{\mathcal{F}},\mathbb{P}). A quantizer Γ\Gamma for XX is stationary if the nearest neighbour projection X^Γ=ProjΓ(X)\widehat{X}^{\Gamma}=\mathrm{Proj}_{\Gamma}(X) satisfies

𝔼[X|X^Γ]=X^Γ.\mathbb{E}\left[X|\widehat{X}^{\Gamma}\right]=\widehat{X}^{\Gamma}. (16)
Remark 3.9.

Taking expectation on both sides of (16) yields 𝔼[X]=𝔼[𝔼[X|X^Γ]]=𝔼[X^Γ].\mathbb{E}[X]=\mathbb{E}[\mathbb{E}[X|\widehat{X}^{\Gamma}]]=\mathbb{E}[\widehat{X}^{\Gamma}]. Furthermore, for any convex function f:df:\mathbb{R}^{d}\to\mathbb{R}, the identity above, the conditional Jensen’s inequality and the tower property yield

𝔼[f(X^Γ)]=𝔼[f(𝔼[X|X^Γ])]𝔼[𝔼[f(X)|X^Γ]]=𝔼[f(X)].\mathbb{E}[f(\widehat{X}^{\Gamma})]=\mathbb{E}[f(\mathbb{E}[X|\widehat{X}^{\Gamma}])]\leq\mathbb{E}[\mathbb{E}[f(X)|\widehat{X}^{\Gamma}]]=\mathbb{E}[f(X)].

While an optimal quadratic quantizer of order NN of a random vector is always stationary [35, Proposition 1(c)], the converse is not true in general. We now present the corresponding definition for a stochastic process.

Definition 3.10.

Let {Xt}t[T1,T2]\{X_{t}\}_{t\in[T_{1},T_{2}]} be a stochastic process on (Ω,,{t}t[T1,T2],)(\Omega,{\mathcal{F}},\{{\mathcal{F}}_{t}\}_{t\in[T_{1},T_{2}]},\mathbb{P}). We say that an NN-quantizer ΛN:={λ1N,,λNN}L2[T1,T2]\Lambda^{N}:=\{\lambda_{1}^{N},\cdots,\lambda_{N}^{N}\}\subset L^{2}[T_{1},T_{2}], inducing the quantization X^=X^ΛN\widehat{X}=\widehat{X}^{\Lambda^{N}}, is stationary if 𝔼[Xt|X^t]=X^t\mathbb{E}[X_{t}|\widehat{X}_{t}]=\widehat{X}_{t}, for all t[T1,T2]t\in[T_{1},T_{2}].

Remark 3.11.

To ease the notation, we omit the grid ΛN\Lambda^{N} in X^ΛN\widehat{X}^{\Lambda^{N}}, while the dependence on the dimension NN remains via the superscript 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N} (recall (12)).

As was stated in Section 2.1, we are working with the continuous version of the Gaussian Volterra process ZZ given by the series expansion (6). This will ease the proof of stationarity below (for a similar result in the case of the Brownian motion [35, Proposition 2]).

Proposition 3.12.

The product functional quantizers inducing Z^𝐝\widehat{Z}^{\boldsymbol{\mathrm{d}}} in (12) are stationary.

Proof.

For any t[0,1]t\in[0,1], by linearity, we have the following chain of equalities:

𝔼[Zt|{ξ^nd(n)}1nm]=𝔼[k1𝒦[ψk](t)ξk|{ξ^nd(n)}1nm]=k1𝒦[ψk](t)𝔼[ξk|{ξ^nd(n)}1nm].\mathbb{E}\left[Z_{t}|\{\widehat{\xi}_{n}^{d(n)}\}_{1\leq n\leq m}\right]=\mathbb{E}\left[\sum_{k\geq 1}\mathcal{K}[\psi_{k}](t)\xi_{k}\Big{|}\{\widehat{\xi}_{n}^{d(n)}\}_{1\leq n\leq m}\right]=\sum_{k\geq 1}\mathcal{K}[\psi_{k}](t)\mathbb{E}\left[\xi_{k}\Big{|}\{\widehat{\xi}_{n}^{d(n)}\}_{1\leq n\leq m}\right].

Since the 𝒩(0,1)\mathcal{N}(0,1)-Gaussian ξn\xi_{n}’s are i.i.d., by definition of optimal quadratic quantizers (hence stationary), we have 𝔼[ξk|ξ^id(i)]=δikξ^id(i)\mathbb{E}[\xi_{k}|\widehat{\xi}_{i}^{d(i)}]=\delta_{ik}\widehat{\xi}_{i}^{d(i)}, for all i,k{1,,m}i,k\in\{1,\dots,m\}, and therefore

𝔼[ξk|{ξ^nd(n)}1nm]=𝔼[ξk|ξ^kd(k)]=ξ^kd(k), for all k{1,,m}.\mathbb{E}\left[\xi_{k}\Big{|}\{\widehat{\xi}_{n}^{d(n)}\}_{1\leq n\leq m}\right]=\mathbb{E}\left[\xi_{k}\Big{|}\widehat{\xi}_{k}^{d(k)}\right]=\widehat{\xi}_{k}^{d(k)},\text{ for all }k\in\{1,\dots,m\}.

Thus, we obtain

𝔼[Zt|{ξ^nd(n)}1nm]=k1𝒦[ψk](t)ξ^kd(k)=Z^t𝐝.\mathbb{E}\left[Z_{t}\Big{|}\{\widehat{\xi}_{n}^{d(n)}\}_{1\leq n\leq m}\right]=\sum_{k\geq 1}\mathcal{K}[\psi_{k}](t)\widehat{\xi}_{k}^{d(k)}=\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t}.

Finally, exploiting the tower property and the fact that the σ\sigma-algebra generated by Z^t𝐝\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t} is included in the σ\sigma-algebra generated by {ξ^nd(n)}n{1,,m}\{\widehat{\xi}_{n}^{d(n)}\}_{n\in\{1,\dots,m\}} by Definition 3.2, we obtain

𝔼[Zt|Z^t𝐝]\displaystyle\mathbb{E}\left[Z_{t}\Big{|}\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t}\right] =𝔼[𝔼[Zt|{ξ^nd(n)}n{1,,m}]|Z^t𝐝]=𝔼[Z^t𝐝|Z^t𝐝]=Z^t𝐝,\displaystyle=\mathbb{E}\left[\mathbb{E}\left[Z_{t}\Big{|}\{\widehat{\xi}_{n}^{d(n)}\}_{n\in\{1,\dots,m\}}\right]\Big{|}\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t}\right]=\mathbb{E}\left[\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t}\Big{|}\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t}\right]=\widehat{Z}^{\boldsymbol{\mathrm{d}}}_{t},

which concludes the proof. ∎

4. Application to VIX derivatives in rough Bergomi

We now specialize the setup above to the case of rough volatility models. These models are extensions of classical stochastic volatility models, introduced to better reproduce the market implied volatility surface. The volatility process is stochastic and driven by a rough process, by which we mean a process whose trajectories are HH-Hölder continuous with H(0,12)H\in(0,\frac{1}{2}). The empirical study [18] was the first to suggest such a rough behaviour for the volatility, and ignited tremendous interest in the topic. The website https://sites.google.com/site/roughvol/home contains an exhaustive and up-to-date review of the literature on rough volatility. Unlike continuous Markovian stochastic volatility models, which are not able to fully describe the steep implied volatility skew of short-maturity options in equity markets, rough volatility models have shown accurate fit for this crucial feature. Within rough volatility, the rough Bergomi model [4] is one of the simplest, yet decisive frameworks to harness the power of the roughness for pricing purposes. We show how to adapt our functional quantization setup to this case.

4.1. The generalized Bergomi model

We work here with a slightly generalised version of the rough Bergomi model, defined as

{Xt=120t𝒱s𝑑s+0t𝒱s𝑑Bs,X0=0,𝒱t=v0(t)exp{γZtγ220tK(ts)2𝑑s},𝒱0>0,\left\{\begin{array}[]{rcll}X_{t}&=&\displaystyle-\frac{1}{2}\int_{0}^{t}{{\mathcal{V}_{s}}}ds+\int_{0}^{t}\sqrt{{{\mathcal{V}_{s}}}}dB_{s},&X_{0}=0,\\ \mathcal{V}_{t}&=&\displaystyle v_{0}(t)\exp\left\{\gamma Z_{t}-\frac{\gamma^{2}}{2}\int_{0}^{t}K(t-s)^{2}ds\right\},&\mathcal{V}_{0}>0,\end{array}\right. (17)

where XX is the log-stock price, 𝒱\mathcal{V} the instantaneous variance process driven by the Gaussian Volterra process ZZ in (1), γ>0\gamma>0 and BB is a Brownian motion defined as B:=ρW+1ρ2WB:=\rho W+\sqrt{1-\rho^{2}}W^{\perp} for some correlation ρ[1,1]\rho\in[-1,1] and W,WW,W^{\perp} orthogonal Brownian motions. The filtered probability space is therefore taken as t=tWtW{\mathcal{F}}_{t}={\mathcal{F}}_{t}^{W}\lor{\mathcal{F}}_{t}^{W^{\perp}}, t0t\geq 0. This is a non-Markovian generalization of Bergomi’s second generation stochastic volatility model [9], letting the variance be driven by a Gaussian Volterra process instead of a standard Brownian motion. Here, vT(t)v_{T}(t) denotes the forward variance for a remaining maturity tt, observed at time TT. In particular, v0v_{0} is the initial forward variance curve, assumed to be 0{\mathcal{F}}_{0}-measurable. Indeed, given market prices of variance swaps σT2(t)\sigma_{T}^{2}(t) at time TT with remaining maturity tt, the forward variance curve can be recovered as vT(t)=ddt(tσT2(t))v_{T}(t)=\frac{d}{dt}\left(t\sigma_{T}^{2}(t)\right), for all t0t\geq 0, and the process {vs(ts)}0st\{v_{s}(t-s)\}_{0\leq s\leq t} is a martingale for all fixed t>0t>0.

Remark 4.1.

With K(u)=uH12K(u)=u^{H-\frac{1}{2}}, γ=2νCH\gamma=2\nu C_{H}, for ν>0\nu>0, and CH:=2HΓ(3/2H)Γ(H+1/2)Γ(22H)C_{H}:=\sqrt{\frac{2H\Gamma(3/2-H)}{\Gamma(H+1/2)\Gamma(2-2H)}}, we recover the standard rough Bergomi model [4].

4.2. VIX Futures in the generalized Bergomi

We consider the pricing of VIX Futures (www.cboe.com/tradable_products/vix/) in the rough Bergomi model. They are highly liquid Futures on the Chicago Board Options Exchange Volatility Index, introduced on March 26, 2004, to allow for trading in the underlying VIX. Each VIX Future represents the expected implied volatility for the 30 days following the expiration date of the Futures contract itself. The continuous version of the VIX at time TT is determined by the continuous-time monitoring formula

VIXT2:\displaystyle\mathrm{VIX}_{T}^{2}: =𝔼T[1ΔTT+ΔdXs,Xs]=1ΔTT+Δ𝔼[𝒱s|T]𝑑s\displaystyle=\mathbb{E}_{T}\left[\frac{1}{\Delta}\int_{T}^{T+\Delta}d\langle X_{s},X_{s}\rangle\right]=\frac{1}{\Delta}\int_{T}^{T+\Delta}\mathbb{E}[\mathcal{V}_{s}|{\mathcal{F}}_{T}]ds (18)
=1ΔTT+Δ𝔼T[v0(s)eγZsγ220sK(su)2𝑑u]𝑑s\displaystyle=\frac{1}{\Delta}\int_{T}^{T+\Delta}\mathbb{E}_{T}\left[v_{0}(s)e^{\gamma Z_{s}-\frac{\gamma^{2}}{2}\int_{0}^{s}K(s-u)^{2}du}\right]ds (19)
=1ΔTT+Δv0(s)eγ0TK(su)𝑑Wuγ220sK(su)2𝑑u𝔼T[eγTsK(su)𝑑Wu]𝑑s\displaystyle=\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(s)e^{\gamma\int_{0}^{T}K(s-u)dW_{u}-\frac{\gamma^{2}}{2}\int_{0}^{s}K(s-u)^{2}du}\mathbb{E}_{T}\left[e^{\gamma\int_{T}^{s}K(s-u)dW_{u}}\right]ds (20)
=1ΔTT+Δv0(s)eγ0TK(su)𝑑Wuγ220sK(su)2𝑑ueγ22TsK(su)2𝑑u𝑑s,\displaystyle=\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(s)e^{\gamma\int_{0}^{T}K(s-u)dW_{u}-\frac{\gamma^{2}}{2}\int_{0}^{s}K(s-u)^{2}du}e^{\frac{\gamma^{2}}{2}\int_{T}^{s}K(s-u)^{2}du}ds, (21)

similarly to [26], where Δ\Delta is equal to 3030 days, and we write 𝔼T[]:=𝔼[|T]\mathbb{E}_{T}[\cdot]:=\mathbb{E}[\cdot|{\mathcal{F}}_{T}] (dropping the subscript when T=0T=0). Thus, the price of a VIX Future with maturity TT is given by

𝒫T:=𝔼[VIXT]=𝔼[(1ΔTT+Δv0(t)eγZtT,Δ+γ22(0tTK(s)2𝑑s0tK(s)2𝑑s)𝑑t)12],\mathcal{P}_{T}:=\mathbb{E}\left[\mathrm{VIX}_{T}\right]=\mathbb{E}\left[\left(\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(t)e^{\gamma Z^{T,\Delta}_{t}+\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)}dt\right)^{\frac{1}{2}}\right], (22)

where the process (ZtT,Δ)t[T,T+Δ](Z^{T,\Delta}_{t})_{t\in[T,T+\Delta]} is given by

ZtT,Δ=0TK(ts)𝑑Ws,t[T,T+Δ].Z^{T,\Delta}_{t}=\int_{0}^{T}K(t-s)dW_{s},\qquad t\in[T,T+\Delta]. (23)

To develop a functional quantization setup for VIX Futures, we need to quantize the process ZT,ΔZ^{T,\Delta}, which is close, yet slightly different, from the Gaussian Volterra process ZZ in (1).

4.3. Properties of ZTZ^{T}

To retrieve the same setting as above, we normalize the time interval to [0,1][0,1], that is T+Δ=1T+\Delta=1. Then, for TT fixed, we define the process ZT:=ZT,1TZ^{T}:=Z^{T,1-T} as

ZtT:=0TK(ts)𝑑Ws,t[T,1],Z^{T}_{t}:=\int_{0}^{T}K(t-s)dW_{s},\qquad t\in[T,1], (24)

which is well defined by the square integrability of KK. By definition, the process ZTZ^{T} is centered Gaussian and Itô isometry gives its covariance function as

RZT(t,s)=0TK(tu)K(su)𝑑u,t,s[T,1].R_{Z^{T}}(t,s)=\int_{0}^{T}K(t-u)K(s-u)du,\qquad t,s\in[T,1].

Proceeding as previously, we introduce a Gaussian process with same mean and covariance as those of ZTZ^{T}, represented as a series expansion involving standard Gaussian random variables; from which product functional quantization follows. It is easy to see that the process ZTZ^{T} has continuous trajectories. Indeed, (ZtTZsT)2𝔼[|ZtZs|2|TW](Z_{t}^{T}-Z_{s}^{T})^{2}\leq\mathbb{E}[|Z_{t}-Z_{s}|^{2}|\mathcal{F}^{W}_{T}], by conditional Jensen’s inequality since ZtT=𝔼[Zt|TW]Z^{T}_{t}=\mathbb{E}[Z_{t}|\mathcal{F}^{W}_{T}]. Then, applying tower property, for any Ts<t1T\leq s<t\leq 1,

𝔼[|ZtTZsT|2]𝔼[|ZtZs|2],\displaystyle\mathbb{E}\left[\left|Z^{T}_{t}-Z^{T}_{s}\right|^{2}\right]\leq\mathbb{E}\left[|Z_{t}-Z_{s}|^{2}\right],

and therefore the H-Hölder regularity of ZZ (Section 2) implies that of ZTZ^{T}.

4.3.1. Series expansion

Let {ξn}n1\{\xi_{n}\}_{n\geq 1} be an i.i.d. sequence of standard Gaussian and {ψn}n1\{\psi_{n}\}_{n\geq 1} the orthonormal basis of L2[0,1]L^{2}[0,1] from (4). Denote by 𝒦T()\mathcal{K}^{T}(\cdot) the operator from L2[0,1]L^{2}[0,1] to 𝒞[T,1]\mathcal{C}[{T},1] that associates to each fL2[0,1]f\in L^{2}[0,1],

𝒦T[f](t):=0TK(ts)f(s)𝑑s,t[T,1].\mathcal{K}^{T}[f](t):=\int_{0}^{T}K(t-s)f(s)ds,\quad t\in[T,1]. (25)

We define the process YTY^{T} as (recall the analogous  (3)):

YtT:=n1𝒦T[ψn](t)ξn,t[T,1].Y^{T}_{t}:=\sum_{n\geq 1}\mathcal{K}^{T}[\psi_{n}](t)\xi_{n},\qquad t\in[T,1]. (26)

The lemma below follows from the corresponding results in Remark 2.2 and Lemma 2.4:

Lemma 4.2.

The process YTY^{T} is centered, Gaussian and with covariance function

RYT(s,t):=𝔼[YsTYtT]=0TK(tu)K(su)𝑑u,for all s,t[T,1].R_{Y^{T}}(s,t):=\mathbb{E}\left[Y^{T}_{s}Y^{T}_{t}\right]=\int_{0}^{T}K(t-u)K(s-u)du,\qquad\text{for all }s,t\in[T,1].

To complete the analysis of ZTZ^{T}, we require an analogue version of Assumption 2.3.

Assumption 4.3.

Assumption 2.3 holds for the sequence (𝒦T[ψn])n1(\mathcal{K}^{T}[\psi_{n}])_{n\geq 1} on [T,1][T,1] with the constants C1C_{1} and C2C_{2} depending on TT.

4.4. The truncated RL case

We again pay special attention to the RL case, for which the operator (25) reads, for each nn\in\mathbb{N},

𝒦HT[ψn](t):=0T(ts)H12ψn(s)𝑑s,for all t[T,1],\mathcal{K}_{H}^{T}[\psi_{n}](t):=\int_{0}^{T}(t-s)^{H-\frac{1}{2}}\psi_{n}(s)ds,\qquad\text{for all }t\in[T,1],

and satisfies the following, proved in Appendix A.4:

Lemma 4.4.

The functions {𝒦HT[ψn]}n1\{\mathcal{K}_{H}^{T}[\psi_{n}]\}_{n\geq 1} satisfy Assumption 4.3.

A key role in this proof is played by an intermediate lemma, proved in Appendix A.3, which provides a convenient representation for the integral 0T(tu)H12e𝚒πu𝑑u\int_{0}^{T}(t-u)^{H-\frac{1}{2}}\mathrm{e}^{\mathtt{i}\pi u}du, tT0t\geq T\geq 0, in terms of the generalised Hypergeometric function F21(){}_{1}{F}_{2}(\cdot).

Lemma 4.5.

For any tT0t\geq T\geq 0, the representation

0T(tu)H12e𝚒πu𝑑u=e𝚒πt[(ζ12(t,h1)ζ12((tT),h1))𝚒π(ζ32(t,h2)ζ32((tT),h2))]\int_{0}^{T}(t-u)^{H-\frac{1}{2}}\mathrm{e}^{\mathtt{i}\pi u}du={\mathrm{e}^{\mathtt{i}\pi t}}\left[\Big{(}\zeta_{\frac{1}{2}}(t,h_{1})-\zeta_{\frac{1}{2}}((t-T),h_{1})\Big{)}-\mathtt{i}\pi\Big{(}\zeta_{\frac{3}{2}}(t,h_{2})-\zeta_{\frac{3}{2}}((t-T),h_{2})\Big{)}\right]

holds, where h1:=12(H+12)h_{1}:=\frac{1}{2}(H+\frac{1}{2}) and h2=12+h1h_{2}=\frac{1}{2}+h_{1}, χ(z):=14π2z2\chi(z):=-\frac{1}{4}\pi^{2}z^{2} and

ζk(z,h):=z2h2hF21(h;k,1+h;χ(z)),for k{12,32}.\zeta_{k}(z,h):=\displaystyle\frac{z^{2h}}{2h}{}_{1}{F}_{2}\left(h;k,1+h;\chi(z)\right),\qquad\text{for }k\in\left\{\frac{1}{2},\frac{3}{2}\right\}. (27)
Remark 4.6.

The representation in Lemma 4.5 can be exploited to obtain an explicit formula for 𝒦HT[ψn](t)\mathcal{K}_{H}^{T}[\psi_{n}](t), t[T,1]t\in[T,1] and nn\in\mathbb{N}:

𝒦HT[ψn](t)=2mH+120mT(mtu)H12cos(πu)𝑑u=2mH+12{0mT(mtu)H12e𝚒πu𝑑u}\displaystyle\mathcal{K}_{H}^{T}[\psi_{n}](t)=\frac{\sqrt{2}}{m^{H+\frac{1}{2}}}\int_{0}^{mT}(mt-u)^{H-\frac{1}{2}}\cos(\pi u)du=\frac{\sqrt{2}}{m^{H+\frac{1}{2}}}\Re\left\{\int_{0}^{mT}(mt-u)^{H-\frac{1}{2}}\mathrm{e}^{\mathtt{i}\pi u}du\right\}
=2mH+12{e𝚒πmt[(ζ12(mt,h1)ζ12(m(tT),h1))𝚒π(ζ32(mt,h2)ζ32(m(tT),h2))]}\displaystyle=\frac{\sqrt{2}}{m^{H+\frac{1}{2}}}\Re\Bigg{\{}{\mathrm{e}^{\mathtt{i}\pi mt}}\bigg{[}\Big{(}\zeta_{\frac{1}{2}}(mt,h_{1})-\zeta_{\frac{1}{2}}(m(t-T),h_{1})\Big{)}-\mathtt{i}\pi\Big{(}\zeta_{\frac{3}{2}}(mt,h_{2})-\zeta_{\frac{3}{2}}(m(t-T),h_{2})\Big{)}\bigg{]}\Bigg{\}}
=2mH+12{cos(mtπ)(ζ12(mt,h1)ζ12(m(tT),h1))+πsin(mtπ)(ζ32(mt,h2)ζ32(m(tT),h2))},\displaystyle=\frac{\sqrt{2}}{m^{H+\frac{1}{2}}}\bigg{\{}\cos(mt\pi)\Big{(}\zeta_{\frac{1}{2}}(mt,h_{1})-\zeta_{\frac{1}{2}}(m(t-T),h_{1})\Big{)}+\pi\sin(mt\pi)\Big{(}\zeta_{\frac{3}{2}}(mt,h_{2})-\zeta_{\frac{3}{2}}(m(t-T),h_{2})\Big{)}\bigg{\}},

with m:=n12m:=n-\frac{1}{2} and ζ12()\zeta_{\frac{1}{2}}(\cdot), ζ32()\zeta_{\frac{3}{2}}(\cdot) in (27). We shall exploit this in our numerical simulations.

4.5. VIX Derivatives Pricing

We can now introduce the quantization for the process ZT,ΔZ^{T,\Delta}, similarly to Definition 3.2, recalling the definition of the set 𝒟mN{\mathcal{D}}_{m}^{N} in (11):

Definition 4.7.

A product functional quantization for ZT,ΔZ^{T,\Delta} of order NN is defined as

Z^tT,Δ,𝐝:=n=1m𝒦T,Δ[ψnT,Δ](t)ξ^nd(n),t[T,T+Δ],\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}}_{t}:=\sum_{n=1}^{m}\mathcal{K}^{T,\Delta}[\psi_{n}^{T,\Delta}](t)\widehat{\xi}_{n}^{d(n)},\qquad t\in[T,T+\Delta], (28)

where 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N}, for some mm\in\mathbb{N}, and for every n{1,,m},n\in\{1,\dots,m\}, ξ^nd(n)\widehat{\xi}_{n}^{d(n)} is the (unique) optimal quadratic quantization of the Gaussian variable ξn\xi_{n} of order d(n)d(n).

The sequence {ψnT,Δ}n\{\psi_{n}^{T,\Delta}\}_{n\in\mathbb{N}} denotes the orthonormal basis of L2[0,T+Δ]L^{2}[0,T+\Delta] given by

ψnT,Δ(t)=2T+Δcos(tλn(T+Δ)), with λn=4(2n1)2π2,\psi_{n}^{T,\Delta}(t)=\sqrt{\frac{{2}}{T+\Delta}}\cos\left(\frac{t}{\sqrt{\lambda_{n}}(T+\Delta)}\right),\quad\text{ with }\lambda_{n}=\frac{4}{(2n-1)^{2}\pi^{2}}, (29)

and the operator 𝒦T,Δ:L2[0,T+Δ]𝒞[T,T+Δ]\mathcal{K}^{T,\Delta}:{L}^{2}[0,T+\Delta]\to\mathcal{C}[T,T+\Delta] is defined for fL2[0,T+Δ]f\in{L}^{2}[0,T+\Delta] as

𝒦T,Δ[f](t):=0TK(ts)f(s)𝑑s,t[T,T+Δ].\mathcal{K}^{T,\Delta}[f](t):=\int_{0}^{T}K(t-s)f(s)ds,\qquad t\in[T,T+\Delta].

Adapting the proof of Proposition 3.12 it is possible to prove that these quantizers are stationary, too.

Remark 4.8.

The dependence on Δ\Delta is due to the fact that the coefficients in the series expansion depend on the time interval [T,T+Δ][T,T+\Delta].

In the RL case for each nn\in\mathbb{N}, we can write, using Remark 4.6, for any t[T,T+Δ]t\in[T,T+\Delta]:

𝒦HT,Δ[ψnT,Δ](t)\displaystyle\mathcal{K}_{H}^{T,\Delta}[\psi_{n}^{T,\Delta}](t) =2T+Δ0T(ts)H12cos(sλn(T+Δ))𝑑s,\displaystyle=\sqrt{\frac{{2}}{T+\Delta}}\int_{0}^{T}(t-s)^{H-\frac{1}{2}}\cos\left(\frac{s}{\sqrt{\lambda_{n}}(T+\Delta)}\right)ds,
=2(T+Δ)H(n1/2)H+120(n1/2)T+ΔT((n1/2)T+Δtu)H12cos(πu)𝑑u\displaystyle=\frac{\sqrt{2}(T+\Delta)^{H}}{(n-1/2)^{H+\frac{1}{2}}}\int_{0}^{\frac{(n-1/2)}{T+\Delta}T}\Big{(}\frac{(n-1/2)}{T+\Delta}t-u\Big{)}^{H-\frac{1}{2}}\cos(\pi u)du
=2(T+Δ)H(n12)H+12{cos((n12)T+Δtπ)(ζ12((n12)T+Δt,h1)ζ12((n12)T+Δ(tT),h1))\displaystyle=\frac{\sqrt{2}(T+\Delta)^{H}}{(n-\frac{1}{2})^{H+\frac{1}{2}}}\bigg{\{}\cos\Big{(}\frac{(n-\frac{1}{2})}{T+\Delta}t\pi\Big{)}\Big{(}\zeta_{\frac{1}{2}}\Big{(}\frac{(n-\frac{1}{2})}{T+\Delta}t,h_{1}\Big{)}-\zeta_{\frac{1}{2}}\Big{(}\frac{(n-\frac{1}{2})}{T+\Delta}(t-T),h_{1}\Big{)}\Big{)}
+πsin((n12)T+Δtπ)(ζ32((n12)T+Δt,h2)ζ32((n12)T+Δ(tT),h2))}.\displaystyle\qquad+\pi\sin\Big{(}\frac{(n-\frac{1}{2})}{T+\Delta}t\pi\Big{)}\Big{(}\zeta_{\frac{3}{2}}\Big{(}\frac{(n-\frac{1}{2})}{T+\Delta}t,h_{2}\Big{)}-\zeta_{\frac{3}{2}}\Big{(}\frac{(n-\frac{1}{2})}{T+\Delta}(t-T),h_{2}\Big{)}\Big{)}\bigg{\}}.

We thus exploit Z^T,Δ,𝐝\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}} to obtain an estimation of VIXT\mathrm{VIX}_{T} and of VIX Futures through the following

VIX^T𝐝:=(1ΔTT+Δv0(t)exp{γZ^tT,Δ,𝐝+γ22(0tTK(s)2𝑑s0tK(s)2𝑑s)}𝑑t)12,\displaystyle\qquad\widehat{\mathrm{VIX}}_{T}^{\boldsymbol{\mathrm{d}}}:=\left(\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(t)\exp\left\{\gamma\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}}_{t}+\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)\right\}dt\right)^{\frac{1}{2}}, (30)
𝒫^T𝐝:=𝔼[(1ΔTT+Δv0(t)exp{γZ^tT,Δ,𝐝+γ22(0tTK(s)2𝑑s0tK(s)2𝑑s)}𝑑t)12].\displaystyle\qquad\widehat{\mathcal{P}}_{T}^{\boldsymbol{\mathrm{d}}}:=\mathbb{E}\left[\left(\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(t)\exp\left\{\gamma\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}}_{t}+\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)\right\}dt\right)^{\frac{1}{2}}\right]. (31)
Remark 4.9.

The expectation above reduces to the following deterministic summation, making its computation immediate:

𝒫^T𝐝\displaystyle\widehat{\mathcal{P}}_{T}^{\boldsymbol{\mathrm{d}}} =\displaystyle= 𝔼[(1ΔTT+Δv0(t)eγn=1m𝒦T,Δ[ψnT,Δ](t)ξ^nd(n)+γ22(0tTK(s)2𝑑s0tK(s)2𝑑s)𝑑t)12]\displaystyle\mathbb{E}\left[\left(\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(t)e^{\gamma\sum_{n=1}^{m}\mathcal{K}^{T,\Delta}[\psi_{n}^{T,\Delta}](t)\widehat{\xi}_{n}^{d(n)}+\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)}dt\right)^{\frac{1}{2}}\right]
=\displaystyle= i¯Id(1ΔTT+Δv0(t)eγn=1m𝒦T,Δ[ψnT,Δ](t)xind(n)+γ22(0tTK(s)2𝑑s0tK(s)2𝑑s)𝑑t)12\displaystyle\sum_{\underline{i}\in I^{d}}\left(\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(t)e^{\gamma\sum_{n=1}^{m}\mathcal{K}^{T,\Delta}[\psi_{n}^{T,\Delta}](t)x_{i_{n}}^{d(n)}+\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)}dt\right)^{\frac{1}{2}}
n=1m(ξnCin(Γd(n))),\displaystyle\cdot\prod_{n=1}^{m}\mathbb{P}(\xi_{n}\in C_{i_{n}}(\Gamma^{d(n)})),

where ξ^nd(n)\widehat{\xi}_{n}^{d(n)} is the (unique) optimal quadratic quantization of ξn\xi_{n} of order d(n)d(n), Cj(Γd(n))C_{j}(\Gamma^{d(n)}) is the jj-th Voronoi cell relative to the d(n)d(n)-quantizer (Definition 3.1), with j=1,,d(n)j=1,\dots,d(n) and i¯=(i1,,im)j=1m{1,,d(j)}\underline{i}=(i_{1},\dots,i_{m})\in\prod_{j=1}^{m}\{1,\dots,d(j)\}. In the numerical illustrations displayed in Section 5, we exploited Simpson rule to evaluate these integrals. In particular, we used simps function from scipy.integrate with 300300 points.

4.6. Quantization error of VIX Derivatives

The following L2L^{2}-error estimate is a consequence of Assumption 4.3 (B) and its proof is omitted since it is analogous to that of Proposition 3.6:

Proposition 4.10.

Under Assumption 4.3, for any N1N\geq 1, there exist mT(N)m_{T}^{*}(N)\in\mathbb{N}, C>0C>0 such that

𝔼[Z^T,Δ,𝐝T,NZT,ΔL2[T,T+Δ]2]12Clog(N)H,\mathbb{E}\left[\left\|\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}^{*}_{T,N}}-Z^{T,\Delta}\right\|^{2}_{L^{2}[T,T+\Delta]}\right]^{\frac{1}{2}}\leq C\log(N)^{-H},

for 𝐝T,N𝒟mT(N)N\boldsymbol{\mathrm{d}}^{*}_{T,N}\in{\mathcal{D}}_{m_{T}^{*}(N)}^{N} and with, for each n=1,,mT(N)n=1,\dots,m_{T}^{*}(N),

dT,N(n)=N1mT(N)n(H+12)(mT(N)!)2H+12mT(N).d^{*}_{T,N}(n)=\Big{\lfloor}N^{\frac{1}{m_{T}^{*}(N)}}n^{-(H+\frac{1}{2})}\left(m_{T}^{*}(N)!\right)^{\frac{2H+1}{2m_{T}^{*}(N)}}\Big{\rfloor}.

Furthermore mT(N)=𝒪(log(N))m_{T}^{*}(N)=\mathcal{O}(\log(N)).

As a consequence, we have the following error quantification for European options on the VIX:

Theorem 4.11.

Let F:F:\mathbb{R}\to\mathbb{R} be a globally Lipschitz-continuous function and 𝐝m\boldsymbol{\mathrm{d}}\in\mathbb{N}^{m} for some mm\in\mathbb{N}. There exists 𝔎>0\mathfrak{K}>0 such that

|𝔼[F(VIXT)]𝔼[F(VIX^T𝐝)]|𝔎𝔼[ZT,ΔZ^T,Δ,𝐝L2([T,T+Δ])2]12.\left|\mathbb{E}\left[F\left(\mathrm{VIX}_{T}\right)\right]-\mathbb{E}\left[F\left(\widehat{\mathrm{VIX}}_{T}^{\boldsymbol{\mathrm{d}}}\right)\right]\right|\leq\mathfrak{K}\ \mathbb{E}\left[\left\|Z^{T,\Delta}-\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}}\right\|_{L^{2}([T,T+\Delta])}^{2}\right]^{\frac{1}{2}}. (32)

Furthermore, for any N1N\geq 1, there exist mT(N)m_{T}^{*}(N)\in\mathbb{N} and >0\mathfrak{C}>0 such that, with 𝐝T,N𝒟mT(N)N\boldsymbol{\mathrm{d}}^{*}_{T,N}\in{\mathcal{D}}_{m_{T}^{*}(N)}^{N},

|𝔼[F(VIXT)]𝔼[F(VIX^T𝐝T,N)]|log(N)H.\left|\mathbb{E}\left[F\left(\mathrm{VIX}_{T}\right)\right]-\mathbb{E}\left[F\left(\widehat{\mathrm{VIX}}_{T}^{\boldsymbol{\mathrm{d}}^{*}_{T,N}}\right)\right]\right|\leq\mathfrak{C}\log(N)^{-H}. (33)

The upper bound in (33) is an immediate consequence of  (32) and Proposition 4.10. The proof of (32) is much more involved and is postponed to Appendix A.5.

Remark 4.12.
  • When F(x)=1F(x)=1, we obtain the price of VIX Futures and the quantization error

    |𝒫T𝒫^T𝐝|𝔎𝔼[ZT,ΔZ^T,Δ,𝐝L2([T,T+Δ])2]12,\left|\mathcal{P}_{T}-\widehat{\mathcal{P}}_{T}^{\boldsymbol{\mathrm{d}}}\right|\leq\mathfrak{K}\ \mathbb{E}\left[\left\|Z^{T,\Delta}-\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}}\right\|_{L^{2}([T,T+\Delta])}^{2}\right]^{\frac{1}{2}},

    and, for any N1N\geq 1, Theorem 4.11 yields the existence of mT(N)m_{T}^{*}(N)\in\mathbb{N}, >0\mathfrak{C}>0 such that

    |𝒫T𝒫^T𝐝T,N|log(N)H.\left|\mathcal{P}_{T}-\widehat{\mathcal{P}}_{T}^{\boldsymbol{\mathrm{d}}^{*}_{T,N}}\right|\leq\mathfrak{C}\log(N)^{-H}.
  • Since the functions F(x):=(xK)+F(x):=(x-K)_{+} and F(x):=(Kx)+F(x):=(K-x)_{+} are globally Lipschitz continuous, the same bounds apply for European Call and Put options on the VIX.

5. Numerical results for the RL case

We now test the quality of the quantization on the pricing of VIX Futures in the standard rough Bergomi model, considering the RL kernel in Remark 4.1.

5.1. Practical considerations for mm and 𝕕\mathbb{d}

Proposition 3.6 provides, for any fixed N,N\in\mathbb{N}, some indications on m(N)m^{*}(N) and 𝐝N𝒟mN\boldsymbol{\mathrm{d}}^{*}_{N}\in\mathcal{D}_{m}^{N} (see (11)), for which the rate of convergence of the quantization error is log(N)H\log(N)^{-H}. We present now a numerical algorithm to compute the optimal parameters. For a given number of trajectories NN\in\mathbb{N}, the problem is equivalent to finding mm\in\mathbb{N} and 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N} such that 𝔼[ZHZ^H,𝐝L2[0,1]2]\mathbb{E}[\|Z^{H}-\widehat{Z}^{H,\boldsymbol{\mathrm{d}}}\|_{L^{2}[0,1]}^{2}] is minimal. Starting from (A.1) and adding and subtracting the quantity n=1m(01𝒦H[ψn](t)2𝑑t)\sum_{n=1}^{m}(\int_{0}^{1}\mathcal{K}_{H}[\psi_{n}](t)^{2}dt), we obtain

𝔼[ZHZ^H,𝐝L2[0,1]2]\displaystyle\mathbb{E}\left[\left\|Z^{H}-\widehat{Z}^{H,\boldsymbol{\mathrm{d}}}\right\|_{L^{2}[0,1]}^{2}\right] =n=1m(01𝒦H[ψn](t)2𝑑t)[𝜺d(n)(ξn)]2+km+101𝒦H[ψk](t)2𝑑t\displaystyle=\sum_{n=1}^{m}\left(\int_{0}^{1}\mathcal{K}_{H}[\psi_{n}](t)^{2}dt\right)[\boldsymbol{\varepsilon}^{d(n)}(\xi_{n})]^{2}+\sum_{k\geq m+1}\int_{0}^{1}\mathcal{K}_{H}[\psi_{k}](t)^{2}dt
=n=1m(01𝒦H[ψn](t)2𝑑t){[𝜺d(n)(ξn)]21}+k101𝒦H[ψk](t)2𝑑t,ix\displaystyle=\sum_{n=1}^{m}\left(\int_{0}^{1}\mathcal{K}_{H}[\psi_{n}](t)^{2}dt\right)\left\{\left[\boldsymbol{\varepsilon}^{d(n)}(\xi_{n})\right]^{2}-1\right\}+\sum_{k\geq 1}\int_{0}^{1}\mathcal{K}_{H}[\psi_{k}](t)^{2}dt,ix (34)

where 𝜺d(n)(ξn)\boldsymbol{\varepsilon}^{d(n)}(\xi_{n}) denotes the optimal quadratic quantization error for the quadratic quantizer of order d(n)d(n) of the standard Gaussian random variable ξn\xi_{n} (see Appendix A.1 for more details). Notice that the last term on the right-hand side of (5.1) does not depend on mm, nor on 𝐝\boldsymbol{\mathrm{d}}. We therefore simply look for mm and 𝐝\boldsymbol{\mathrm{d}} that minimize

A(m,𝐝):=n=1m(01𝒦H[ψn]2(t)𝑑t)([𝜺d(n)(ξn)]21).A(m,\boldsymbol{\mathrm{d}}):=\sum_{n=1}^{m}\left(\int_{0}^{1}\mathcal{K}_{H}[\psi_{n}]^{2}(t)dt\right)\left([\boldsymbol{\varepsilon}^{d(n)}(\xi_{n})]^{2}-1\right).

This can be easily implemented: the functions 𝒦H[ψn]\mathcal{K}_{H}[\psi_{n}] can be obtained numerically from the Hypergeometric function and the quadratic errors 𝜺d(n)(ξn)\boldsymbol{\varepsilon}^{d(n)}(\xi_{n}) are available at www.quantize.maths-fi.com/gaussian_database, for d(n){1,,5999}d(n)\in\{1,\dots,5999\}. The algorithm therefore reads as follows

  1. (i)

    fix mm;

  2. (ii)

    minimize A(m,𝐝)A(m,\boldsymbol{\mathrm{d}}) over 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in\mathcal{D}_{m}^{N} and call it A~(m)\widetilde{A}(m);

  3. (iii)

    minimize A~(m)\widetilde{A}(m) over mm\in\mathbb{N}.

The results of the algorithm for some reference values of NN\in\mathbb{N} are available in Table 1, where N¯traj:=i=1m¯(N)d¯N(i)\overline{N}_{traj}:=\prod_{i=1}^{\overline{m}(N)}\overline{d}_{N}(i) represents the number of trajectories actually computed in the optimal case. In Table 2, we compute the rate optimal parameters derived in Proposition 3.6: the column ‘Relative error’ contains the normalized difference between the L2L^{2}-quantization error made with the optimal choice of m¯(N)\overline{m}(N) and 𝐝¯N\overline{\boldsymbol{\mathrm{d}}}_{N} in Table 1 and the L2L^{2}-quantization error made with m(N)m^{*}(N) and 𝐝N{\boldsymbol{\mathrm{d}}}^{*}_{N} of the corresponding line of the table, namely |ZHZ^H,𝐝¯NL2[0,1]ZHZ^H,𝐝NL2[0,1]|ZHZ^H,𝐝¯NL2[0,1]\frac{|\|Z^{H}-\widehat{Z}^{H,\overline{\boldsymbol{\mathrm{d}}}_{N}}\|_{L^{2}[0,1]}-\|Z^{H}-\widehat{Z}^{H,\boldsymbol{\mathrm{d}}^{*}_{N}}\|_{L^{2}[0,1]}|}{\|Z^{H}-\widehat{Z}^{H,\overline{\boldsymbol{\mathrm{d}}}_{N}}\|_{L^{2}[0,1]}}. In the column Ntraj:=i=1m(N)dN(i){N}^{*}_{traj}:=\prod_{i=1}^{{m}^{*}(N)}{d}^{*}_{N}(i) we display the number of trajectories actually computed in the rate-optimal case. The optimal quadratic vector quantization of a standard Gaussian of order 11 is the random variable identically equal to zero and so when d(i)=1d(i)=1 the corresponding term is uninfluential in the representation.

Table 1. Optimal parameters.
NN m¯(N)\overline{m}(N) 𝐝¯N\overline{\boldsymbol{\mathrm{d}}}_{N} N¯traj\overline{N}_{traj}
1010 22 55 - 22 1010
10210^{2} 44 88 - 33 - 22 - 22 9696
10310^{3} 66 1010 - 44 - 33 - 22 - 22 - 22 960960
10410^{4} 88 1010 - 55 - 44 - 33 - 22 - 22 - 22 - 22 96009600
10510^{5} 1010 1414 - 66 - 44 - 33 - 33 - 22 - 22 - 22 - 22 - 22 9676896768
10610^{6} 1212 1414 - 66 - 55 - 44 - 33 - 33 - 22 - 22 - 22 - 22 - 22 - 22 967680967680
Table 2. Rate-optimal parameters.
NN m(N)=log(N)m^{*}(N)=\lfloor\log(N)\rfloor Relative error 𝐝N\boldsymbol{\mathrm{d}}^{*}_{N} NtrajN_{traj}^{*}
1010 22 2.75%2.75\% 33 - 22 66
10210^{2} 44 1.30%1.30\% 55 - 33 - 22 - 22 6060
10310^{3} 66 1.09%1.09\% 66 - 44 - 33 - 22 - 22 - 22 576576
10410^{4} 99 3.08%3.08\% 66 - 44 - 33 - 22 - 22 - 22 - 22 - 11 - 11 11521152
10510^{5} 1111 3.65%3.65\% 77 - 44 - 33 - 33 - 22 - 22 - 22 - 22 - 11 - 11 - 11 40324032
10610^{6} 1313 2.80%2.80\% 88 - 55 - 44 - 33 - 33 - 22 - 22 - 22 - 22 - 22 - 11 - 11 - 11 4608046080
NN m(N)=log(N)m^{*}(N)=\lfloor\log(N)\rfloor - 1 Relative error 𝐝N\boldsymbol{\mathrm{d}}^{*}_{N} NtrajN_{traj}^{*}
1010 11 2.78%2.78\% 1010 1010
10210^{2} 33 1.13%1.13\% 66 - 44 - 33 7272
10310^{3} 55 1.22%1.22\% 77 - 44 - 33 - 33 - 22 504504
10410^{4} 88 1.35%1.35\% 77 - 44 - 33 - 33 - 22 - 22 - 22 - 22 40324032
10510^{5} 1010 2.29%2.29\% 77 - 55 - 44 - 33 - 22 - 22 - 22 - 22 - 22 - 11 1344013440
10610^{6} 1212 2.25%2.25\% 88 - 55 - 44 - 33 - 33 - 22 - 22 - 22 - 22 - 22 - 22 - 11 9216092160
NN m(N)=log(N)m^{*}(N)=\lfloor\log(N)\rfloor - 2 Relative error 𝐝N\boldsymbol{\mathrm{d}}^{*}_{N} NtrajN_{traj}^{*}
10210^{2} 22 2.53%2.53\% 1212 - 88 9696
10310^{3} 44 1.44%1.44\% 99 - 55 - 44 - 33 540540
10410^{4} 77 1.46%1.46\% 77 - 55 - 44 - 33 - 22 - 22 - 22 33603360
10510^{5} 99 1.57%1.57\% 88 - 55 - 44 - 33 - 33 - 22 - 22 - 22 - 22 2304023040
10610^{6} 1111 1.48%1.48\% 99 - 66 - 44 - 33 - 33 - 33 - 22 - 22 - 22 - 22 - 22 186624186624

5.2. The functional quantizers

The computations in Section 2 and 3 for the RL process, respectively the ones in Section 4.3 and 4.4 for ZH,TZ^{H,T}, provide a way to obtain the functional quantizers of the processes.

5.2.1. Quantizers of the RL process

For the RL process, Definition 3.4 shows that its quantizer is a weighted Cartesian product of grids of the one-dimensional standard Gaussian random variables. The time-dependent weights 𝒦H[ψn]()\mathcal{K}_{H}[\psi_{n}](\cdot) are computed using (9), and for a fixed number of trajectories NN, suitable m¯(N)\overline{m}(N) and 𝐝¯N𝒟m¯(N)N\overline{\boldsymbol{\mathrm{d}}}_{N}\in{\mathcal{D}}_{\overline{m}(N)}^{N} are chosen according to the algorithm in Section 5.1. Not surprisingly, Figures 1 show that as the paths of the process get smoother (HH increases) the trajectories become less fluctuating and shrink around zero. For H=0.5H=0.5, where the RL process reduces to the standard Brownian motion, we recover the well-known quantizer from [35, Figures 7-8]. This is consistent as in that case 𝒦H[ψn](t)=λn2sin(tλn)\mathcal{K}_{H}[\psi_{n}](t)=\sqrt{\lambda_{n}}\sqrt{2}\sin\left(\frac{t}{\sqrt{\lambda_{n}}}\right), and so YHY^{H} is the Karhuenen-Loève expansion for the Brownian motion [35, Section 7.1].

Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Figure 1. Product functional quantizations of the RL process with N-quantizers, for H{0.1,0.25,0.5}H\in\{0.1,0.25,0.5\}, for N=10N=10 and N=100N=100.

5.2.2. Quantizers of ZH,TZ^{H,T}

A quantizer for ZH,TZ^{H,T} is defined analogously to that of ZHZ^{H} using Definition 3.4. The weights 𝒦HT[ψn]()\mathcal{K}_{H}^{T}[\psi_{n}](\cdot) in the summation are available in closed form, as shown in Remark 4.6. It is therefore possible to compute the NN-product functional quantizer, for any NN\in\mathbb{N}, as Figure  2 displays.

Refer to caption
Refer to caption
Figure 2. Product functional quantization of ZH,TZ^{H,T} via N-quantizers, with H=0.1H=0.1, T=0.7T=0.7, for N{10,100}N\in\{10,100\}.

5.3. Pricing and comparison with Monte Carlo

In this section we show and comment some plots related to the estimation of prices of derivatives on the VIX and realized variance. We set the values H=0.1H=0.1 and ν=1.18778\nu=1.18778 for the parameters and investigate three different initial forward variance curves v0()v_{0}(\cdot), as in [26]:

  • Scenario 1. v0(t)=0.2342v_{0}(t)=0.234^{2};

  • Scenario 2. v0(t)=0.2342(1+t)2v_{0}(t)=0.234^{2}(1+t)^{2};

  • Scenario 3. v0(t)=0.23421+tv_{0}(t)=0.234^{2}\sqrt{1+t}.

The choice of such ν\nu is a consequence of the choice η=1.9\eta=1.9, consistently with [8], and of the relationship ν=η2H2CH\nu=\eta\frac{\sqrt{2H}}{2C_{H}}. In all these cases, v0v_{0} is an increasing function of time, whose value at zero is close to the square of the reference value of 0.250.25.

5.3.1. VIX Futures Pricing

One of the most recent and effective way to compute the price of VIX Futures is a Monte-Carlo-simulation method based on Cholesky decomposition, for which we refer to [26, Section 3.3.2]. It can be considered as a good approximation of the true price when the number MM of computed paths is large. In fact, in [26] the authors tested three simulation-based methods (Hybrid scheme + forward Euler, Truncated Cholesky, SVD decomposition) and ‘all three methods seem to approximate the prices similarly well’. We thus consider the truncated Cholesky approach as a benchmark and take M=106M=10^{6} trajectories and 300300 equidistant point for the time grid.

Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Refer to caption
Figure 3. VIX Futures prices (left) and relative error (right) computed with quantization and with Monte-Carlo as a function of the maturity TT, for different numbers of trajectories, for each forward variance curve scenario.

In Figure 3, we plot the VIX Futures prices as a function of the maturity TT, where TT ranges in {1,2,3,6,9,12}\{1,2,3,6,9,12\} months (consistently with actual quotations) on the left, and the corresponding relative error w.r.t. the Monte Carlo benchmark on the right. It is clear that the quantization approximates the benchmark from below and that the accuracy increases with the number of trajectories.
We highlight that the quantization scheme for VIX Futures can be sped up considerably by storing ahead the quantized trajectories for ZH,T,ΔZ^{H,T,\Delta}, so that we only need to compute the integrations and summations in Remark 4.9, which are extremely fast.

Table 3. Grid organization times (in seconds) as a function of the maturity (rows, in months) and of the number of trajectories (columns).
Grid organisation time
10210^{2} 10310^{3} 10410^{4} 10510^{5} 10610^{6}
1 0.474 0.491 0.99 4.113 37.183
2 0.476 0.487 0.752 4.294 39.134
3 0.617 0.536 0.826 4.197 37.744
6 0.474 0.475 0.787 4.432 37.847
9 0.459 0.6 0.858 3.73 41.988
12 0.498 0.647 1.016 3.995 38.045

Furthermore, the grid organization time itself is not that significant. In Table 3 we display the grid organization times (in seconds) as a function of the maturity (rows) expressed in months and of the number of trajectories (columns). From this table one might deduce that the time needed for the organization of the grids is suitable to be performed once per day (say every morning) as it should be for actual pricing purposes. It is interesting to note that the estimations obtained with quantization (which is an exact method) are consistent in that they mimick the trend of benchmark prices over time even for very small values of NN. However, as a consequence of the variance in the estimations, the Monte Carlo prices are almost useless for small values of MM. Moreover, improving the estimations with Monte Carlo requires to increase the number of points in the time grid with clear impact on computational time, while this is not the case with quantization since the trajectories in the quantizers are smooth. Indeed, the trajectories in the quantizers are not only smooth but also almost constant over time, hence reducing the number of time steps to get the desired level of accuracy. Notice that here we may refer also to the issue of complexity related to discretization: a quadrature formula over nn points has a cost 𝒪(n)\mathcal{O}(n), while the simulation with a Cholesky method over the same grid has cost 𝒪(n2)\mathcal{O}(n^{2}). Finally, our quantization method does not require much RAM. Indeed, all the simulations performed with quantization can be easily run on a personal laptop111The personal computer used to run the quantization codes has the following technical specifications: RAM: 8.00 GB, SSD memory: 512 GB, Processor: AMD Ryzen 7 4700U with Radeon Graphics 2.00 GHz., while this is not the case for the Monte Carlo scheme proposed here222The computer used to run the Monte Carlo codes is a virtual machine (OpenStack/Nova/KVM/Qemu, www.openstack.org) with the following technical specifications: RAM: 32.00 GB, CPU: 8 virtual cores, Hypervisor CPU: Intel(R) Xeon(R) CPU E5-2650 v3 @ 2.30GHz, RAM 128GB, Storage: CEPH cluster (www.ceph.com).. For the sake of completeness, we also recall that combining Monte Carlo pricing of VIX futures/options with an efficient control variate speeds up the computations significantly [23].

Refer to caption
Refer to caption
Figure 4. Log-log (natural logarithm) plot of the empirical absolute error with the theoretically predicted one for Scenario 1, with T{1,12}T\in\{1,12\} months.

In Figure 4, we show some plots comparing the behaviour of the empirical error with the theoretically predicted one. We have decided to display only a couple of maturities for the first scenario since the other plots are very similar. The figures display in a clear way that the order of convergence of the empirical error should be bigger than the theoretically predicted one: in particular, we expect it to be 𝒪(log(N)1)\mathcal{O}(\log(N)^{-1}).

5.3.2. VIX Options Pricing

To complete the discussion on VIX Options pricing, we present in Figure 5 the approximation of the prices of ATM Call Options on the VIX obtained via quantization as a function of the maturity TT and for different numbers of trajectories against the same price computed via Monte Carlo simulations with M=106M=10^{6} trajectories and 300300 equidistant point for the time grid, as a benchmarch. Each plot represents a different scenario for the initial forward variance curve. For all scenarios, as the number NN of trajectories goes to infinity, the prices in Figure 5 are clearly converging, and the limiting curve is increasing in the maturity, as it should be.

all Refer to caption

Refer to caption
Refer to caption
Figure 5. Prices of ATM Call Options on the VIX via quantization.

5.3.3. Pricing of Continuously Monitored Options on Realized Variance

Product functional quantization of the process (ZtH)t[0,T](Z_{t}^{H})_{t\in[0,T]} can be exploited for (meaningful) pricing purposes, too. We first price variance swaps, whose price is given by the following expression

𝔖T:=𝔼[1T0T𝒱tdt|0].\mathfrak{S}_{T}:=\mathbb{E}\left[\frac{1}{T}\int_{0}^{T}\mathcal{V}_{t}\text{d}t\bigg{|}\mathcal{F}_{0}\right].

Let us recall that, in the rough Bergomi model,

𝒱t=v0(t)exp(2νCHZtHν2CH2Ht2H),\mathcal{V}_{t}=v_{0}(t)\exp\Big{(}2\nu C_{H}Z_{t}^{H}-\frac{\nu^{2}C_{H}^{2}}{H}t^{2H}\Big{)},

where CH=2HΓ(3/2H)Γ(H+1/2)Γ(22H)C_{H}=\sqrt{\frac{2H\Gamma(3/2-H)}{\Gamma(H+1/2)\Gamma(2-2H)}}, ν>0\nu>0 is an endogenous constant and v0(t)v_{0}(t) being the initial forward variance curve. Thus, exploiting the fact that, for any fixed t[0,T]t\in[0,T], ZtHZ^{H}_{t} is distributed according to a centred Gaussian random variable with variance 0t(ts)2H1𝑑s=t2H2H\int_{0}^{t}(t-s)^{2H-1}ds=\frac{t^{2H}}{2H}, the quantity 𝔖T\mathfrak{S}_{T} can be explicitly computed:

𝔖T=1T0Tv0(t)dt.\mathfrak{S}_{T}=\frac{1}{T}\int_{0}^{T}v_{0}(t)\text{d}t.

This is particularly handy and provides us a simple benchmark. The price 𝔖T\mathfrak{S}_{T} is, then, approximated via quantization through

𝔖^T𝕕=i¯Id(1T0Tv0(t)exp(2νCHn=1m𝒦H[ψn](t)xind(n)ν2CH2Ht2H)𝑑t)n=1m(ξnCin(Γd(n))).\displaystyle\widehat{\mathfrak{S}}_{T}^{\mathbb{d}}=\sum_{\underline{i}\in I^{d}}\left(\frac{1}{T}\int_{0}^{T}v_{0}(t)\exp\left(2\nu C_{H}\sum_{n=1}^{m}\mathcal{K}_{H}[\psi_{n}](t)x_{i_{n}}^{d(n)}-\frac{\nu^{2}C_{H}^{2}}{H}t^{2H}\right)dt\right)\prod_{n=1}^{m}\mathbb{P}(\xi_{n}\in C_{i_{n}}(\Gamma^{d(n)})). (35)

Numerical results are presented in Figure 6. On the left-hand side we display a table with the approximations (depending on NN, the number of trajectories) of the price of a swap on the realized variance in Scenario 1, for T=1T=1, and the true value v0=0.2342v_{0}=0.234^{2}. On the right-hand side a log-log (natural logarithm) plot of the error against the function clog(N)Hc\log(N)^{-H}, with cc being a suitable positive constant. For variance swaps the error is not performing very well. It is indeed very close to the upper bound clog(N)Hc\log(N)^{-H} that we have computed theoretically. One possible theoretical motivation for this behaviour lies in the difference between strong and weak error rates. Weak error and strong error do not necessarily share the same order of convergence, being the weak error faster in general. See [5, 6, 17] for recent developments on the topic in the rough volatility framework. For pricing purposes, we are interested in weak error rates. Indeed, the pricing error should in principle have the following form 𝔼[f(ZH)]𝔼[f(Z^H)]\mathbb{E}[f(Z^{H})]-\mathbb{E}[f(\widehat{Z}^{H})], where Z^H\widehat{Z}^{H} is the process that we are using to approximate the original ZHZ^{H} and ff is a functional that comes from the payoff function and that we can interpret as a test function. Thus, the functional ff has a smoothing effect. On the other hand, the upper bound for the quantization error we have computed is a strong error rate. This theoretical discrepancy motivates the findings in Figure 4 when pricing VIX Futures and other options on the VIX: the empirical error seems to converge with order 𝒪(log(N)1)\mathcal{O}(\log(N)^{-1}), while the predicted order is 𝒪(log(N)H)\mathcal{O}(\log(N)^{-H}). The different empirical rates that are seen in Figure 4 for VIX futures (roughly 𝒪(log(N)1)\mathcal{O}(\log(N)^{-1}))) and in Figure 6 for variance swaps (much closer to 𝒪(log(N)H)\mathcal{O}(\log(N)^{-H})) could be also related to the different degree of pathwise regularity of the processes ZZ and ZTZ^{T} . While tZt=0tK(ts)𝑑Wst\rightarrow Z_{t}=\int_{0}^{t}K(t-s)dWs is a.s. (Hϵ)(H-\epsilon)-Hölder, for fixed TT, the trajectories tZtT=0TK(ts)𝑑Wst\rightarrow Z_{t}^{T}=\int_{0}^{T}K(t-s)dWs of ZTZ^{T} are much smoother when t(T,T+Δ)t\in(T,T+\Delta) and tt is bounded away from TT. When pricing VIX derivatives, we are quantizing almost everywhere a smooth Gaussian process (hence error rate of order log(N)1)\log(N)^{-1}), while when pricing derivatives on realized variance, we are applying quantization to a rough Gaussian process (hence error rate of order 𝒪(log(N)H)\mathcal{O}(\log(N)^{-H})), resulting in a deteriorated accuracy for the prices of realized volatility derivatives such as the variance swaps in Figure 6.

Furthermore, it can be easily shown that, for any 𝕕𝒟mN\mathbb{d}\in\mathcal{D}_{m}^{N} and for any m,Nm,N\in\mathbb{N}, with m<Nm<N, 𝔖^T𝕕\widehat{\mathfrak{S}}_{T}^{\mathbb{d}} always provides a lower bound for the true price 𝔖T{\mathfrak{S}}_{T}. Indeed, since the quantizers Z^H,𝕕\widehat{Z}^{H,\mathbb{d}} of the process ZHZ^{H} are stationary (cfr. Proposition 3.12), an application of Remark 3.9 to the convex function f(x)=exp(2νCHx)f(x)=\exp(2\nu C_{H}x) together with the positivity of v0(t)exp(ν2CH2t2HH)v_{0}(t)\exp(-\frac{\nu^{2}C_{H}^{2}t^{2H}}{H}), for any t[0,T]t\in[0,T], yields

𝔖^T𝕕\displaystyle\widehat{\mathfrak{S}}_{T}^{\mathbb{d}} =𝔼[1T0Tv0(t)exp(ν2CH2t2HH)exp(2νCHZ^TH,𝕕)dt|0]\displaystyle=\mathbb{E}\left[\frac{1}{T}\int_{0}^{T}v_{0}(t)\exp\left(-\frac{\nu^{2}C_{H}^{2}t^{2H}}{H}\right)\exp\left(2\nu C_{H}\widehat{Z}_{T}^{H,\mathbb{d}}\right)\text{d}t\bigg{|}\mathcal{F}_{0}\right]
=1T0Tv0(t)exp(ν2CH2t2HH)𝔼0[exp(2νCHZ^TH,𝕕)]dt\displaystyle=\frac{1}{T}\int_{0}^{T}v_{0}(t)\exp\left(-\frac{\nu^{2}C_{H}^{2}t^{2H}}{H}\right)\mathbb{E}_{0}\left[\exp\left(2\nu C_{H}\widehat{Z}_{T}^{H,\mathbb{d}}\right)\right]\text{d}t
1T0Tv0(t)exp(ν2CH2t2HH)𝔼0[exp(2νCHZTH)]dt=𝔖T.\displaystyle\leq\frac{1}{T}\int_{0}^{T}v_{0}(t)\exp\left(-\frac{\nu^{2}C_{H}^{2}t^{2H}}{H}\right)\mathbb{E}_{0}\left[\exp\left(2\nu C_{H}{Z}_{T}^{H}\right)\right]\text{d}t=\mathfrak{S}_{T}.
True price 0.05480.0548
Quantization, N=102N=10^{2} 0.02300.0230
Quantization, N=103N=10^{3} 0.02460.0246
Quantization, N=104N=10^{4} 0.02570.0257
Quantization, N=105N=10^{5} 0.02660.0266
Quantization, N=106N=10^{6} 0.02730.0273
Refer to caption
Figure 6. Prices and errors for variance swaps.

To complete this section, we plot in Figure 7 approximated prices of European Call Options on the realized variance via quantization with N{102,103,104,105,106}N\in\{10^{2},10^{3},10^{4},10^{5},10^{6}\} trajectories and via Monte Carlo with M=106M=10^{6} trajectories, as a benchmark. In order to take advantage of the trajectories obtained, we compute the price of a realized variance Call option with strike KK and maturity T=1T=1 as

𝒞(K,T)=𝔼[(1T0T𝒱tdtK)+|0],\mathcal{C}(K,T)=\mathbb{E}\left[\left(\frac{1}{T}\int_{0}^{T}\mathcal{V}_{t}\text{d}t-K\right)_{+}\bigg{|}\mathcal{F}_{0}\right],

and we approximate it via quantization through

𝒞^𝕕(K,T)=i¯Id(1T0Tv0(t)exp(2νCHn=1m𝒦H[ψn](t)xind(n)ν2CH2Ht2H)𝑑tK)+n=1m(ξnCin(Γd(n))).\widehat{\mathcal{C}}^{\mathbb{d}}(K,T)=\sum_{\underline{i}\in I^{d}}\Bigg{(}\frac{1}{T}\int_{0}^{T}v_{0}(t)\exp\left(2\nu C_{H}\sum_{n=1}^{m}\mathcal{K}_{H}[\psi_{n}](t)x_{i_{n}}^{d(n)}-\frac{\nu^{2}C_{H}^{2}}{H}t^{2H}\right)dt-K\Bigg{)}_{+}\prod_{n=1}^{m}\mathbb{P}(\xi_{n}\in C_{i_{n}}(\Gamma^{d(n)})).

The three plots in Figure 7 display the behaviour of the price of a European Call on the realized variance as a function of the strike price KK (close to the ATM value) for the three scenarios considered before.

Refer to caption
Refer to caption
Refer to caption
Figure 7. Prices of European Call Option on realized variance computed via Monte Carlo with M=106M=10^{6} trajectories and via quantization with N{102,103,104,105,106}N\in\{10^{2},10^{3},10^{4},10^{5},10^{6}\} trajectories, as a function of KK.

5.3.4. Quantization and MC comparison

In order to make a fair comparison between quantization and Monte Carlo simulations, we present a figure to display, for each methodology, the computational work needed for a given error tolerance for the pricing of VIX Futures. The plots in Figure 8 should be read as follows. First, for any M,N{102,103,104,105,106}M,N\in\{10^{2},10^{3},10^{4},10^{5},10^{6}\}, we have computed the corresponding pricing errors: εMC(M):=|PriceMC(M)RefPrice|\varepsilon^{MC}(M):=|\textrm{Price}^{MC}(M)-\textrm{RefPrice}| and εQ(N):=|PriceQ(N)RefPrice|\varepsilon^{Q}(N):=|\textrm{Price}^{Q}(N)-\textrm{RefPrice}| where PriceMC(M)\textrm{Price}^{MC}(M) is the Monte Carlo price obtained via truncated Cholesky with MM trajectories, PriceQ(N)\textrm{Price}^{Q}(N) is the price computed via quantization with NN trajectories and RefPrice comes from the lowerbound in Equation (3.4)(3.4) in [26] and the associated computational time in seconds tMC(M)t^{MC}(M) and tQ(N)t^{Q}(N), respectively for Monte Carlo simulation and quantization. Then, each point in the plot is associated either to a value of MM in case of Monte Carlo (the circles in Figure 8), or NN in case of quantization (the triangles in Figure 8), and its xx-coordinate provides the absolute value of the associated pricing error, while its yy-coordinate represents the associated computational cost in seconds.

These plots lead to the following observations:

  • For quantization, which is an exact method, the error is strictly monotone in the number of trajectories.

  • When a small number of trajectories is considered, quantization provides a lower error with respect to Monte Carlo, at a comparable cost.

  • For large numbers of trajectories Monte Carlo overcomes quantization both in terms of accuracy and of computational time.

To conclude, quantization can always be run with an arbitrary number of trajectories and furthermore for N{102,103,104}N\in\{10^{2},10^{3},10^{4}\} it leads to a lower error with respect to Monte Carlo, at a comparable computational cost, as it is visible from Figure 8. This makes quantization particularly suitable to be used when dealing with standard machines, i.e., laptops with a RAM memory smaller or equal to 1616GB.

Refer to caption
Refer to caption
Figure 8. Computational costs for quantization vs Monte Carlo for Scenario 11, with T=1T=1 month (left-hand side) and T=12T=12 months (right-hand side). The number of trajectories, MM for Monte Carlo and NN for quantization, corresponding to a specific dot is displayed above it.

6. Conclusion

In this paper we provide, on the theoretical side, a precise and detailed result on the convergence of product functional quantizers of Gaussian Volterra processes, showing that the L2L^{2}-error is of order log(N)H\log(N)^{-H}, with NN the number of trajectories and HH the regularity index.

Furthermore, we explicitly characterize the rate optimal parameters, mNm^{*}_{N} and 𝕕N\mathbb{d}^{*}_{N}, and we compare them with the corresponding optimal parameters, m¯N\overline{m}_{N} and 𝕕¯N\overline{\mathbb{d}}_{N}, computed numerically.

In the rough Bergomi model, we apply product functional quantization to the pricing of VIX options, with precise rates of convergence, and of options on realized variance, comparing those – whenever possible – to standard Monte Carlo methods.

The thorough numerical analysis carried out in the paper shows that unfortunately, despite the conceptual promise of functional quantization, while the results on the VIX are very promising, other types of path-dependent options seem to require machine resources way beyond the current requirements of standard Monte Carlo schemes, as shown precisely in the case of variance swaps. While product functional quantization is an exact method, the analysis provided here does not however promise a bright future in the context of rough volatility. It may nevertheless be of practical interest when machine resources are limited and indeed the results for VIX Futures pricing are strongly encouraging in this respect. Functional quantization for rough volatility can however be salvaged when used as a control variate tool to reduce the variance in classical Monte Carlo simulations.

Appendix A Proofs

A.1. Proof of Proposition 3.6

Consider a fixed N1N\geq 1 and (m,𝐝)(m,\boldsymbol{\mathrm{d}}) for 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N}. We have

𝔼[ZZ^𝐝L2[0,1]2]\displaystyle\mathbb{E}\left[\left\|Z-\widehat{Z}^{\boldsymbol{\mathrm{d}}}\right\|_{L^{2}[0,1]}^{2}\right] =𝔼[n1𝒦[ψn]()ξnn=1m𝒦[ψn]()ξ^nd(n)L2[0,1]2]\displaystyle=\mathbb{E}\left[\left\|\sum_{n\geq 1}\mathcal{K}[\psi_{n}](\cdot)\xi_{n}-\sum_{n=1}^{m}\mathcal{K}[\psi_{n}](\cdot)\widehat{\xi}_{n}^{d(n)}\right\|_{L^{2}[0,1]}^{2}\right]
=𝔼[n=1m𝒦[ψn]()(ξnξ^nd(n))+km+1𝒦[ψk]()ξkL2[0,1]2]\displaystyle=\mathbb{E}\left[\left\|\sum_{n=1}^{m}\mathcal{K}[\psi_{n}](\cdot)(\xi_{n}-\widehat{\xi}_{n}^{d(n)})+\sum_{k\geq m+1}\mathcal{K}[\psi_{k}](\cdot)\xi_{k}\right\|_{L^{2}[0,1]}^{2}\right]
=𝔼[01|n=1m𝒦[ψn](t)(ξnξ^nd(n))+km+1𝒦[ψk](t)ξk|2𝑑t]\displaystyle=\mathbb{E}\left[\int_{0}^{1}\left|\sum_{n=1}^{m}\mathcal{K}[\psi_{n}](t)(\xi_{n}-\widehat{\xi}_{n}^{d(n)})+\sum_{k\geq m+1}\mathcal{K}[\psi_{k}](t)\xi_{k}\right|^{2}dt\right]
=01(n=1m𝒦[ψn]2(t)𝔼[|ξnξ^nd(n)|2]+km+1𝒦[ψk]2(t))𝑑t\displaystyle=\int_{0}^{1}\left(\sum_{n=1}^{m}\mathcal{K}[\psi_{n}]^{2}(t)\mathbb{E}\left[|\xi_{n}-\widehat{\xi}_{n}^{d(n)}|^{2}\right]+\sum_{k\geq m+1}\mathcal{K}[\psi_{k}]^{2}(t)\right)dt
=01(n=1m𝒦[ψn]2(t)𝜺d(n)(ξn)2+km+1𝒦[ψk]2(t))𝑑t,\displaystyle=\int_{0}^{1}\left(\sum_{n=1}^{m}\mathcal{K}[\psi_{n}]^{2}(t)\boldsymbol{\varepsilon}^{d(n)}(\xi_{n})^{2}+\sum_{k\geq m+1}\mathcal{K}[\psi_{k}]^{2}(t)\right)dt, (36)

using Fubini’s Theorem and the fact that {ξn}n1\{\xi_{n}\}_{n\geq 1} is a sequence of i.i.d. Gaussian and where 𝜺d(n)(ξn):=inf(α1,,αd(n))d(n)𝔼[min1id(n)|ξnαi|2]\boldsymbol{\varepsilon}^{d(n)}(\xi_{n}):=\inf_{(\alpha_{1},\dots,\alpha_{d(n)})\in\mathbb{R}^{d(n)}}\sqrt{\mathbb{E}[\min_{1\leq i\leq d(n)}|\xi_{n}-\alpha_{i}|^{2}]}. The Extended Pierce Lemma [35, Theorem 1(b)] ensures that 𝜺d(n)(ξn)Ld(n)\boldsymbol{\varepsilon}^{d(n)}(\xi_{n})\leq\frac{L}{d(n)} for a suitable positive constant LL. Exploiting this error bound and the property (B) for 𝒦[ψn]\mathcal{K}[\psi_{n}] in Assumption 2.3, we obtain

𝔼[ZZ^𝐝L2[0,1]2]\displaystyle\mathbb{E}\left[\|Z-\widehat{Z}^{\boldsymbol{\mathrm{d}}}\|_{L^{2}[0,1]}^{2}\right] =n=1m(01𝒦[ψn]2(t)𝑑t)𝜺d(n)(ξn)2+km+101𝒦[ψk]2(t)𝑑t\displaystyle=\sum_{n=1}^{m}\left(\int_{0}^{1}\mathcal{K}[\psi_{n}]^{2}(t)dt\right)\boldsymbol{\varepsilon}^{d(n)}(\xi_{n})^{2}+\sum_{k\geq m+1}\int_{0}^{1}\mathcal{K}[\psi_{k}]^{2}(t)dt (37)
C22{n=1mn(2H+1)𝜺d(n)(ξn)2+km+1k(2H+1)}\displaystyle\leq C_{2}^{2}\left\{\sum_{n=1}^{m}n^{-(2H+1)}\boldsymbol{\varepsilon}^{d(n)}(\xi_{n})^{2}+\sum_{k\geq m+1}k^{-(2H+1)}\right\} (38)
C22{n=1mn(2H+1)L2d(n)2+km+1k(2H+1)}\displaystyle\leq C_{2}^{2}\left\{\sum_{n=1}^{m}n^{-(2H+1)}\frac{L^{2}}{d(n)^{2}}+\sum_{k\geq m+1}k^{-(2H+1)}\right\} (39)
C~(n=1m1n2H+1d(n)2+km+1k(2H+1)),\displaystyle\leq\widetilde{C}\left(\sum_{n=1}^{m}\frac{1}{n^{2H+1}d(n)^{2}}+\sum_{k\geq m+1}k^{-(2H+1)}\right), (40)

with C~=max{L2C22,C22}\widetilde{C}=\max\{L^{2}C_{2}^{2},C_{2}^{2}\}. Inspired by  [31, Section 4.1], we now look for an “optimal” choice of mm\in\mathbb{N} and 𝐝𝒟mN\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N}. This reduces the error in approximating ZZ with a product quantization of the form in (12). Define the optimal product functional quantization  Z^N,\widehat{Z}^{N,\star} of order NN as the Z^𝐝\widehat{Z}^{\boldsymbol{\mathrm{d}}} which realizes the minimal error:

𝔼[ZZ^N,L2[0,1]2]=min{𝔼[ZZ^𝐝L2[0,1]2],m,𝐝𝒟mN}.\mathbb{E}\left[\left\|Z-\widehat{Z}^{N,\star}\right\|_{L^{2}[0,1]}^{2}\right]=\min\left\{\mathbb{E}\left[\left\|Z-\widehat{Z}^{\boldsymbol{\mathrm{d}}}\right\|_{L^{2}[0,1]}^{2}\right],m\in\mathbb{N},\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N}\right\}.

From (37) we deduce

𝔼[ZZ^N,L2[0,1]2]C~infm{km+11k2H+1+inf{n=1m1n2H+1d(n)2,𝐝𝒟mN}}.\mathbb{E}\left[\left\|Z-\widehat{Z}^{N,\star}\right\|_{L^{2}[0,1]}^{2}\right]\leq\widetilde{C}\inf_{m\in\mathbb{N}}\left\{\sum_{k\geq m+1}\frac{1}{k^{2H+1}}+\inf\left\{\sum_{n=1}^{m}\frac{1}{n^{2H+1}d(n)^{2}},\boldsymbol{\mathrm{d}}\in{\mathcal{D}}_{m}^{N}\right\}\right\}. (41)

For any fixed mm\in\mathbb{N} we associate to the internal minimization problem the one we get by relaxing the hypothesis that d(n)d(n)\in\mathbb{N}:

:=inf{n=1m1n2H+1z(n)2,{z(n)}n=1,,m(0,):n=1mz(n)N}.\mathfrak{I}:=\inf\Big{\{}\sum_{n=1}^{m}\frac{1}{n^{2H+1}z(n)^{2}},\{z(n)\}_{n=1,\dots,m}\in(0,\infty):\prod_{n=1}^{m}z(n)\leq N\Big{\}}. (42)

For this infimum, we derive a simple solution exploiting the arithmetic-geometric inequality using Lemma B.2. Setting z~(n):=γN,mn(H+12)\widetilde{z}(n):=\gamma_{N,m}n^{-(H+\frac{1}{2})}, with γN,m:=N1m(j=1mj(2H+1))12m,\gamma_{N,m}:=N^{\frac{1}{m}}\Big{(}\prod_{j=1}^{m}j^{-(2H+1)}\Big{)}^{-\frac{1}{2m}}, n=1,,mn=1,\dots,m, we get

=n=1m1n2H+1z~(n)2=N2mm(n=1mn(2H+1))1m,\mathfrak{I}=\sum_{n=1}^{m}\frac{1}{n^{2H+1}\widetilde{z}(n)^{2}}=N^{-\frac{2}{m}}m\Big{(}\prod_{n=1}^{m}n^{-(2H+1)}\Big{)}^{\frac{1}{m}},

and notice that the sequence {z~(n)}\{\widetilde{z}(n)\} is decreasing. Since ultimately the vector 𝐝\boldsymbol{\mathrm{d}} consists of integers, we use d~(n)=z~(n)\widetilde{d}(n)=\lfloor\widetilde{z}(n)\rfloor, n=1,,mn=1,\dots,m. In fact, this choice guarantees that

n=1md~(n)=n=1mz~(n)n=1mz~(n)=N.\prod_{n=1}^{m}\widetilde{d}(n)=\prod_{n=1}^{m}\lfloor\widetilde{z}(n)\rfloor\leq\prod_{n=1}^{m}\widetilde{z}(n)=N.

Furthermore, setting d~(j)=z~(j)\widetilde{d}(j)=\lfloor\widetilde{z}(j)\rfloor for each j{1,,m}j\in\{1,\dots,m\}, we obtain

d~(j)+1(j(2H+1))12=jH+12(z~(j)+1)jH+12z~(j)=jH+12N1mjH+12{n=1m1n2H+1}12m=N1m{n=1m1n2H+1}12m.\frac{\widetilde{d}(j)+1}{(j^{-(2H+1)})^{\frac{1}{2}}}=j^{H+\frac{1}{2}}(\lfloor\widetilde{z}(j)\rfloor+1)\geq j^{H+\frac{1}{2}}\widetilde{z}(j)=\frac{j^{H+\frac{1}{2}}N^{\frac{1}{m}}}{j^{H+\frac{1}{2}}}\left\{\prod_{n=1}^{m}\frac{1}{n^{2H+1}}\right\}^{-\frac{1}{2m}}=N^{\frac{1}{m}}\left\{\prod_{n=1}^{m}\frac{1}{n^{2H+1}}\right\}^{-\frac{1}{2m}}.

Ordering the terms, we have (d~(j)+1)2N2m(n=1mn(2H+1))1mj(2H+1)(\widetilde{d}(j)+1)^{2}N^{-\frac{2}{m}}\Big{(}\prod_{n=1}^{m}n^{-(2H+1)}\Big{)}^{\frac{1}{m}}\geq j^{-(2H+1)}, for each j{1,,m}j\in\{1,\dots,m\}. From this we deduce the following inequality (notice that the left-hand side term is defined only if d~(1),,d~(m)>0\widetilde{d}(1),\dots,\widetilde{d}(m)>0):

j=1mj(2H+1)d~(j)2\displaystyle\sum_{j=1}^{m}j^{-(2H+1)}\widetilde{d}(j)^{-2} j=1m(d~(j)+1d~(j))2N2m(n=1mn(2H+1))1m\displaystyle\leq\sum_{j=1}^{m}\Big{(}\frac{\widetilde{d}(j)+1}{\widetilde{d}(j)}\Big{)}^{2}N^{-\frac{2}{m}}\Big{(}\prod_{n=1}^{m}n^{-(2H+1)}\Big{)}^{\frac{1}{m}} (43)
=N2m(n=1mn(2H+1))1mj=1m(d~(j)+1d~(j))2\displaystyle=N^{-\frac{2}{m}}\Big{(}\prod_{n=1}^{m}n^{-(2H+1)}\Big{)}^{\frac{1}{m}}\sum_{j=1}^{m}\Big{(}\frac{\widetilde{d}(j)+1}{\widetilde{d}(j)}\Big{)}^{2} (44)
4mN2m(n=1mn(2H+1))1m.\displaystyle\leq 4mN^{-\frac{2}{m}}\Big{(}\prod_{n=1}^{m}n^{-(2H+1)}\Big{)}^{\frac{1}{m}}. (45)

Hence, we are able to make a first error estimation, placing in the internal minimization of the right-hand side of (41) the result of inequality in (43).

𝔼[ZZ^N,L2[0,1]2]\displaystyle\mathbb{E}\left[\left\|Z-\widehat{Z}^{N,\star}\right\|_{L^{2}[0,1]}^{2}\right] C~inf{km+11k2H+1+4mN2m(n=1mn(2H+1))1m,mI(N)}\displaystyle\leq\widetilde{C}\inf\left\{\sum_{k\geq m+1}\frac{1}{k^{2H+1}}+4mN^{-\frac{2}{m}}\left(\prod_{n=1}^{m}n^{-(2H+1)}\right)^{\frac{1}{m}},m\in I(N)\right\} (46)
Cinf{km+11k2H+1+mN2m(n=1mn(2H+1))1m,mI(N)},\displaystyle\leq C^{\prime}\inf\left\{\sum_{k\geq m+1}\frac{1}{k^{2H+1}}+mN^{-\frac{2}{m}}\left(\prod_{n=1}^{m}n^{-(2H+1)}\right)^{\frac{1}{m}},m\in I(N)\right\}, (47)

where C=4C~C^{\prime}=4\widetilde{C} and the set

I(N):={m:N2mm(2H+1)(n=1mn(2H+1))1m1},I(N):=\{m\in\mathbb{N}:N^{\frac{2}{m}}m^{-(2H+1)}\Big{(}\prod_{n=1}^{m}n^{-(2H+1)}\Big{)}^{-\frac{1}{m}}\geq 1\}, (48)

which represents all mm’s such that all d~(1),,d~(m)\widetilde{d}(1),\dots,\widetilde{d}(m) are positive integers. This is to avoid the case where i=1md~(i)N\prod_{i=1}^{m}\widetilde{d}(i)\leq N holds only because one of the factors is zero. In fact, for all n{1,,m}n\in\{1,\dots,m\}, d~(n)=z~(n)\widetilde{d}(n)=\lfloor\widetilde{z}(n)\rfloor is a positive integer if and only if z~(n)1\widetilde{z}(n)\geq 1. Thanks to the monotonicity of {z(n)}n=1,,m\{z(n)\}_{n=1,\dots,m}, we only need to check that

z~(m)=N1mm(H+12)(n=1mn(2H+1))12m1.\widetilde{z}(m)=N^{\frac{1}{m}}m^{-(H+\frac{1}{2})}\Big{(}\prod_{n=1}^{m}n^{-(2H+1)}\Big{)}^{-\frac{1}{2m}}\geq 1.

First, let us show that I(N)I(N), defined in (48) for each N1N\geq 1, is a non-empty finite set with maximum given by m(N)m^{*}(N) of order log(N)\log(N). We can rewrite it as I(N)={m1:amlog(N)}I(N)=\{m\geq 1:a_{m}\leq\log(N)\}, where

an=12log(j=1nn2H+1j2H+1).a_{n}=\frac{1}{2}\log\left(\prod_{j=1}^{n}\frac{n^{2H+1}}{j^{2H+1}}\right). (49)

We can now verify that the sequence ana_{n} is increasing in nn\in\mathbb{N}:

anan+1\displaystyle a_{n}\leq a_{n+1}
\displaystyle\Longleftrightarrow j=1nlog(j(2H+1))nlog(n(2H+1))j=1n+1log(j(2H+1))(n+1)log((n+1)(2H+1))\displaystyle\sum_{j=1}^{n}\log\left(j^{-(2H+1)}\right)-n\log\left(n^{-(2H+1)}\right)\leq\sum_{j=1}^{n+1}\log\left(j^{-(2H+1)}\right)-(n+1)\log\left((n+1)^{-(2H+1)}\right)
\displaystyle\Longleftrightarrow nlog(n(2H+1))log((n+1)(2H+1))(n+1)log((n+1)(2H+1))\displaystyle-n\log\left(n^{-(2H+1)}\right)\leq\log\left((n+1)^{-(2H+1)}\right)-(n+1)\log\left((n+1)^{-(2H+1)}\right)
\displaystyle\Longleftrightarrow log(n(2H+1))log((n+1)(2H+1)),\displaystyle\log\left(n^{-(2H+1)}\right)\geq\log\left((n+1)^{-(2H+1)}\right),

which is obviously true. Furthermore the sequence (an)n(a_{n})_{n} diverges to infinity since

j=1nn(2H+1)j(2H+1)=n(2H+1)nj=1n1j(2H+1)n(2H+1)nj=2n1j(2H+1)n(2H+1)n1n(2H+1)(n1)n(2H+1).\prod_{j=1}^{n}\frac{n^{(2H+1)}}{j^{(2H+1)}}={n^{(2H+1)n}}\prod_{j=1}^{n}\frac{1}{j^{(2H+1)}}\geq{n^{(2H+1)n}}\prod_{j=2}^{n}\frac{1}{j^{(2H+1)}}\geq{n^{(2H+1)n}}\frac{1}{n^{(2H+1)(n-1)}}\geq{n^{(2H+1)}}.

and H(0,12)H\in(0,\frac{1}{2}). We immediately deduce that I(N)I(N) is finite and, since {1}I(N)\{1\}\subset I(N), it is also non-empty. Hence I(N)={1,,m(N)}I(N)=\{1,\dots,m^{*}(N)\}. Moreover, for all N1N\geq 1, am(N)log(N)<am(N)+1a_{m^{*}(N)}\leq\log(N)<a_{m^{*}(N)+1}, which implies that m(N)=𝒪(log(N))m^{*}(N)=\mathcal{O}(\log(N)).
Now, the error estimation in (46) can be further simplified exploiting the fact that, for each mI(N)m\in I(N),

mN2m(n=1mn(2H+1))1m=mm(2H+1)(m(2H+1)N2m(n=1mn(2H+1))1m)1m2H.mN^{-\frac{2}{m}}\left(\prod_{n=1}^{m}n^{-(2H+1)}\right)^{\frac{1}{m}}=mm^{-(2H+1)}\left(m^{-(2H+1)}N^{\frac{2}{m}}\left(\prod_{n=1}^{m}n^{-(2H+1)}\right)^{-\frac{1}{m}}\right)^{-1}\leq m^{-2H}.

The last inequality is a consequence of the fact that (n=1mn(2H+1))1m1\left(\prod_{n=1}^{m}n^{-(2H+1)}\right)^{-\frac{1}{m}}\geq 1 by definition. Hence,

𝔼[ZZ^N,L2[0,1]2]Cinf{km+11k2H+1+m2H,mI(N)},\mathbb{E}\Big{[}\|Z-\widehat{Z}^{N,\star}\|_{L^{2}[0,1]}^{2}\Big{]}\leq{C^{\prime}}\inf\Bigg{\{}\sum_{k\geq m+1}\frac{1}{k^{2H+1}}+m^{-2H},m\in I(N)\Bigg{\}}, (50)

for some suitable constant C>0C^{\prime}>0.

Consider now the sequence {bn}n\{b_{n}\}_{n\in\mathbb{N}}, given by bn=kn+11k2H+1+n2Hb_{n}=\sum_{k\geq n+1}\frac{1}{k^{2H+1}}+n^{-2H}. For n1n\geq 1,

bn+1bn=kn+21k2H+1+1(n+1)2H[kn+11k2H+1+1n2H]=1(n+1)2H+1(n+1)2H+11n2H0,b_{n+1}-b_{n}=\sum_{k\geq n+2}\frac{1}{k^{2H+1}}+\frac{1}{(n+1)^{2H}}-\left[\sum_{k\geq n+1}\frac{1}{k^{2H+1}}+\frac{1}{n^{2H}}\right]=-\frac{1}{(n+1)^{2H}}+\frac{1}{(n+1)^{2H+1}}-\frac{1}{n^{2H}}\leq 0,

so that the sequence is decreasing and the infimum in (50) is attained at m=m(N)m=m^{*}(N). Therefore,

𝔼[ZZ^N,L2[0,1]2]\displaystyle\mathbb{E}\left[\|Z-\widehat{Z}^{N,\star}\|_{L^{2}[0,1]}^{2}\right] \displaystyle\leq Cinf{km+11k2H+1+m2H,mI(N)}\displaystyle{C^{\prime}}\inf\left\{\sum_{k\geq m+1}\frac{1}{k^{2H+1}}+m^{-2H},m\in I(N)\right\}
=\displaystyle= C(km(N)+11k2H+1+m(N)2H)C(m(N)2H1+1+m(N)2H)\displaystyle C^{\prime}\left(\sum_{k\geq m^{*}(N)+1}\frac{1}{k^{2H+1}}+m^{*}(N)^{-2H}\right)\leq C^{\prime}\left({m^{*}(N)^{-2H-1+1}}+m^{*}(N)^{-2H}\right)
=\displaystyle= 2Cm(N)2HClog(N)2H.\displaystyle 2C^{\prime}m^{*}(N)^{-2H}\leq C\log(N)^{-2H}.

A.2. Proof of Remark 2.5

This can be proved specializing the computations done in [32, page 656]. Consider an arbitrary index n1.n\geq 1. For all t,s[0,1]t,s\in[0,1], exploiting Assumption 2.3, we have that, for any ρ[0,1]\rho\in[0,1],

|𝒦[ψn](t)𝒦[ψn](s)|\displaystyle\left|\mathcal{K}[\psi_{n}](t)-\mathcal{K}[\psi_{n}](s)\right| =\displaystyle= |𝒦[ψn](t)𝒦[ψn](s)|ρ|𝒦[ψn](t)𝒦[ψn](s)|1ρ\displaystyle\big{|}\mathcal{K}[\psi_{n}](t)-\mathcal{K}[\psi_{n}](s)\big{|}^{\rho}\big{|}\mathcal{K}[\psi_{n}](t)-\mathcal{K}[\psi_{n}](s)\big{|}^{1-\rho}
\displaystyle\leq (supu,v[0,1],uv|𝒦[ψn](u)𝒦[ψn](v)||uv|H+12|ts|H+12)ρ(2supt[0,1]𝒦[ψn](t))1ρ\displaystyle\left(\sup_{u,v\in[0,1],u\neq v}\frac{|\mathcal{K}[\psi_{n}](u)-\mathcal{K}[\psi_{n}](v)|}{|u-v|^{H+\frac{1}{2}}}|t-s|^{H+\frac{1}{2}}\right)^{\rho}\left(2\sup_{t\in[0,1]}\mathcal{K}[\psi_{n}](t)\right)^{1-\rho}
\displaystyle\leq (C1n)ρ(2C2n(H+12))1ρ|ts|ρ(H+12)=Cρnρ(H+32)(H+12)|ts|ρ(H+12),\displaystyle(C_{1}n)^{\rho}(2C_{2}n^{-(H+\frac{1}{2})})^{1-\rho}|t-s|^{\rho(H+\frac{1}{2})}=C_{\rho}n^{\rho(H+\frac{3}{2})-(H+\frac{1}{2})}|t-s|^{\rho(H+\frac{1}{2})},

where Cρ:=C1ρ(2C2)1ρ<.C_{\rho}:=C_{1}^{\rho}(2C_{2})^{1-\rho}<\infty. Therefore

[𝒦[ψn]]ρ(H+12)=supts[0,1]|𝒦[ψn](t)𝒦[ψn](s)||ts|ρ(H+12)Cρnρ(H+32)(H+12).\left[\mathcal{K}[\psi_{n}]\right]_{\rho(H+\frac{1}{2})}=\sup_{t\neq s\in[0,1]}\frac{\left|\mathcal{K}[\psi_{n}](t)-\mathcal{K}[\psi_{n}](s)\right|}{|t-s|^{\rho(H+\frac{1}{2})}}\leq C_{\rho}n^{\rho(H+\frac{3}{2})-(H+\frac{1}{2})}. (51)

Notice that ρ(H+32)(H+12)<12\rho(H+\frac{3}{2})-(H+\frac{1}{2})<-\frac{1}{2} when ρ[0,HH+3/2]\rho\in[0,\frac{H}{H+3/2}] so that (51) implies

n=1[𝒦[ψn]]ρ(H+12)2Cρ2n=1n2ρ(H+32)2(H+12)Cρ2n=1n(1+ε)=K<.\sum_{n=1}^{\infty}\left[\mathcal{K}[\psi_{n}]\right]_{\rho(H+\frac{1}{2})}^{2}\leq C_{\rho}^{2}\sum_{n=1}^{\infty}n^{2\rho(H+\frac{3}{2})-2(H+\frac{1}{2})}\leq C_{\rho}^{2}\sum_{n=1}^{\infty}n^{-(1+\varepsilon)}=K<\infty. (52)

In particular,

𝔼[|YtYs|2]=n=1|𝒦[ψn](t)𝒦[ψn](s)|2n=1[𝒦[ψn]]ρ(H+12)2|ts|2ρ(H+12)K|ts|2ρ(H+12).\mathbb{E}\left[|Y_{t}-Y_{s}|^{2}\right]=\sum_{n=1}^{\infty}\left|\mathcal{K}[\psi_{n}](t)-\mathcal{K}[\psi_{n}](s)\right|^{2}\leq\sum_{n=1}^{\infty}\left[\mathcal{K}[\psi_{n}]\right]_{\rho(H+\frac{1}{2})}^{2}|t-s|^{2\rho(H+\frac{1}{2})}\leq K|t-s|^{2\rho(H+\frac{1}{2})}.

As noticed in Remark 2.2 the process YY is centered Gaussian. Hence, for each t,s[0,1]t,s\in[0,1] so is YtYsY_{t}-Y_{s}. Proposition B.1 therefore implies that, for any rr\in\mathbb{N},

𝔼[|YtYs|2r]=𝔼[|YtYs|2]r(2r1)!!K|ts|2rρ(H+12),\mathbb{E}\left[|Y_{t}-Y_{s}|^{2r}\right]=\mathbb{E}\left[|Y_{t}-Y_{s}|^{2}\right]^{r}(2r-1)!!\leq K^{\prime}|t-s|^{2r\rho(H+\frac{1}{2})}, (53)

where K=Kr(2r1)!!K^{\prime}=K^{r}(2r-1)!!, yielding existence of a continuous version of YY since choosing rr\in\mathbb{N} such that 2rρ(H+12)>12r\rho(H+\frac{1}{2})>1, Kolmogorov continuity theorem [27, Theorem 3.23] applies directly.

A.3. Proof of Lemma 4.5

Let H+:=H+12H_{+}:=H+\frac{1}{2}. Using [28, Corollary 1, Equation (12)] (with ψ=b2+b1a>1/2\psi=b_{2}+b_{1}-a>1/2), the identity

F21(a,b1,b2,r)=Γ(b1)Γ(b2)Γ(a)π01G2,22,0([b1,b2],[a,12],u)cos(2ru)duu,{}_{1}{F}_{2}(a,b_{1},b_{2},-r)=\frac{\Gamma(b_{1})\Gamma(b_{2})}{\Gamma(a)\sqrt{\pi}}\int_{0}^{1}G_{2,2}^{2,0}\left([b_{1},b_{2}],\left[a,\frac{1}{2}\right],u\right)\cos\left(2\sqrt{ru}\right)\frac{du}{u},

holds for all r>0r>0, where GG denotes the Meijer-G function, generally defined through the so-called Mellin-Barnes type integral [30, Equation (1), Section 5.2]) as

Gp,qm,n([a1,,ap],[b1,,bq],z)=12πiLj=1mΓ(bjs)j=1nΓ(1aj+s)j=m+1qΓ(1bj+s)j=n+1pΓ(ajs)zs𝑑s.G_{p,q}^{\,m,n}\!\left([a_{1},\dots,a_{p}],[b_{1},\dots,b_{q}],z\right)={\frac{1}{2\pi i}}\int_{L}{\frac{\prod_{j=1}^{m}\Gamma(b_{j}-s)\prod_{j=1}^{n}\Gamma(1-a_{j}+s)}{\prod_{j=m+1}^{q}\Gamma(1-b_{j}+s)\prod_{j=n+1}^{p}\Gamma(a_{j}-s)}}\,z^{s}\,ds. (54)

This representation holds if z0z\neq 0, 0mq0\leq m\leq q and 0np0\leq n\leq p, for integers m,n,p,qm,n,p,q, and akbj1,2,3,a_{k}-b_{j}\neq 1,2,3,\dots, for k=1,2,,nk=1,2,\dots,n and j=1,2,,mj=1,2,\dots,m. The last constraint is set to prevent any pole of any Γ(bjs),j=1,2,,m,\Gamma(b_{j}-s),j=1,2,\dots,m, from coinciding with any pole of any Γ(1ak+s),k=1,2,,n\Gamma(1-a_{k}+s),k=1,2,\dots,n. With a>0a>0, b2=1+ab_{2}=1+a and b1=12b_{1}=\frac{1}{2}, since G2,22,0([12,a+1],[a,12],u)=uaG_{2,2}^{2,0}\left(\left[\frac{1}{2},a+1\right],\left[a,\frac{1}{2}\right],u\right)=u^{a}, we can therefore write

01ua1cos(2ru)𝑑u=1aF21(a;12,a+1;r).\int_{0}^{1}u^{a-1}\cos\left(2\sqrt{ru}\right)du=\frac{1}{a}{}_{1}{F}_{2}\left(a;\frac{1}{2},a+1;-r\right). (55)

Similarly, using integration by parts and properties of generalised Hypergeometric functions,

01ua1sin(2ru)𝑑u\displaystyle\int_{0}^{1}u^{a-1}\sin\left(2\sqrt{ru}\right)du =sin(2r)ara01ua12cos(2ru)𝑑u\displaystyle=\frac{\sin(2\sqrt{r})}{a}-\frac{\sqrt{r}}{a}\int_{0}^{1}u^{a-\frac{1}{2}}\cos(2\sqrt{ru})du (56)
=sin(2r)ara(a+12)F21(a+12;12,a+32;r)\displaystyle=\frac{\sin(2\sqrt{r})}{a}-\frac{\sqrt{r}}{a(a+\frac{1}{2})}{}_{1}{F}_{2}\left(a+\frac{1}{2};\frac{1}{2},a+\frac{3}{2};-r\right) (57)
=2ra+12F21(a+12;32,a+32;r),\displaystyle=\frac{2\sqrt{r}}{a+\frac{1}{2}}{}_{1}{F}_{2}\left(a+\frac{1}{2};\frac{3}{2},a+\frac{3}{2};-r\right), (58)

where the last step follows from the definition of generalized sine function sin(z)=zF10(32,14z2)\sin(z)=z\,{}_{0}F_{1}(\frac{3}{2},-\frac{1}{4}z^{2}). Indeed, exploiting (8), we have

sin(2r)a\displaystyle\frac{\sin(2\sqrt{r})}{a} \displaystyle- ra(a+12)F21(a+12;12,a+32;r)\displaystyle\frac{\sqrt{r}}{a(a+\frac{1}{2})}{}_{1}{F}_{2}\left(a+\frac{1}{2};\frac{1}{2},a+\frac{3}{2};-r\right)
=\displaystyle= 2raF10(32,r)ra(a+12)F21(a+12;12,a+32;r)\displaystyle\frac{2\sqrt{r}}{a}{}_{0}F_{1}\left(\frac{3}{2},-r\right)-\frac{\sqrt{r}}{a(a+\frac{1}{2})}{}_{1}{F}_{2}\left(a+\frac{1}{2};\frac{1}{2},a+\frac{3}{2};-r\right)
=\displaystyle= 2ra(a+12)[(a+12)F10(32;r)12F21(a+12;12,a+32;r)]\displaystyle\frac{2\sqrt{r}}{a\left(a+\frac{1}{2}\right)}\left[\left(a+\frac{1}{2}\right){}_{0}F_{1}\left(\frac{3}{2};-r\right)-\frac{1}{2}{}_{1}{F}_{2}\left(a+\frac{1}{2};\frac{1}{2},a+\frac{3}{2};-r\right)\right]
=\displaystyle= 2ra(a+12)[(a+12)k=0(r)kk!(3/2)k12k=0(a+1/2)kk!(1/2)k(a+3/2)k(r)k]\displaystyle\frac{2\sqrt{r}}{a\left(a+\frac{1}{2}\right)}\left[\left(a+\frac{1}{2}\right)\sum_{k=0}^{\infty}\frac{(-r)^{k}}{k!(3/2)_{k}}-\frac{1}{2}\sum_{k=0}^{\infty}\frac{(a+1/2)_{k}}{k!(1/2)_{k}(a+3/2)_{k}}(-r)^{k}\right]
=\displaystyle= 2ra(a+12)k=01k![(a+1/2)(3/2)k1/2(a+1/2)k(1/2)k(a+3/2)k](r)k\displaystyle\frac{2\sqrt{r}}{a\left(a+\frac{1}{2}\right)}\sum_{k=0}^{\infty}\frac{1}{k!}\left[\frac{(a+1/2)}{(3/2)_{k}}-\frac{1/2(a+1/2)_{k}}{(1/2)_{k}(a+3/2)_{k}}\right](-r)^{k}
=\displaystyle= 2ra(a+12)k=01k![a(a+1/2)k(3/2)k(a+3/2)k](r)k\displaystyle\frac{2\sqrt{r}}{a\left(a+\frac{1}{2}\right)}\sum_{k=0}^{\infty}\frac{1}{k!}\left[\frac{a(a+1/2)_{k}}{(3/2)_{k}(a+3/2)_{k}}\right](-r)^{k}
=\displaystyle= 2r(a+12)k=01k!(a+1/2)k(3/2)k(a+3/2)k(r)k=2r(a+12)F21(a+12;32,a+32;r).\displaystyle\frac{2\sqrt{r}}{\left(a+\frac{1}{2}\right)}\sum_{k=0}^{\infty}\frac{1}{k!}\frac{(a+1/2)_{k}}{(3/2)_{k}(a+3/2)_{k}}(-r)^{k}=\frac{2\sqrt{r}}{\left(a+\frac{1}{2}\right)}{}_{1}F_{2}\left(a+\frac{1}{2};\frac{3}{2},a+\frac{3}{2};-r\right).

Letting α:=H12\alpha:=H-\frac{1}{2}, τ:=tT\tau:=t-T, and mapping v:=tuv:=t-u, w:=vtw:=\frac{v}{t} and y:=w2y:=w^{2}, we write

0T(tu)αe𝚒πu𝑑u\displaystyle\int_{0}^{T}(t-u)^{\alpha}\mathrm{e}^{\mathtt{i}\pi u}du =e𝚒πt(tT)tvαe𝚒πv𝑑v=e𝚒πt[0tvαe𝚒πv𝑑v0τvαe𝚒πv𝑑v]\displaystyle=\mathrm{e}^{\mathtt{i}\pi t}\int_{(t-T)}^{t}v^{\alpha}\mathrm{e}^{-\mathtt{i}\pi v}dv=\mathrm{e}^{\mathtt{i}\pi t}\left[\int_{0}^{t}v^{\alpha}\mathrm{e}^{-\mathtt{i}\pi v}dv-\int_{0}^{\tau}v^{\alpha}\mathrm{e}^{-\mathtt{i}\pi v}dv\right]
=e𝚒πt[t1+α01wαe𝚒πwt𝑑wτ1+α01wαe𝚒πwτ𝑑w]\displaystyle=\mathrm{e}^{\mathtt{i}\pi t}\left[t^{1+\alpha}\int_{0}^{1}w^{\alpha}\mathrm{e}^{-\mathtt{i}\pi wt}dw-\tau^{1+\alpha}\int_{0}^{1}w^{\alpha}\mathrm{e}^{-\mathtt{i}\pi w\tau}dw\right]
=e𝚒πt2[t1+α01yα12e𝚒πty𝑑yτ1+α01yα12e𝚒πyτy𝑑y]\displaystyle=\frac{\mathrm{e}^{\mathtt{i}\pi t}}{2}\left[t^{1+\alpha}\int_{0}^{1}y^{\frac{\alpha-1}{2}}\mathrm{e}^{-\mathtt{i}\pi t\sqrt{y}}dy-\tau^{1+\alpha}\int_{0}^{1}y^{\frac{\alpha-1}{2}}\mathrm{e}^{-\mathtt{i}\pi y\tau\sqrt{y}}dy\right]
=e𝚒πt2[I(t)I(τ)],\displaystyle=\frac{\mathrm{e}^{\mathtt{i}\pi t}}{2}\left[I(t)-I(\tau)\right], (59)

where I(z):=z1+α01vα12e𝚒πzv𝑑vI(z):=z^{1+\alpha}\int_{0}^{1}v^{\frac{\alpha-1}{2}}\mathrm{e}^{-\mathtt{i}\pi z\sqrt{v}}dv.
We therefore write, for z{t,τ}z\in\{t,\tau\}, using (55)-(56), πz=2r\pi z=2\sqrt{r}, and identifying a1=α12a-1=\frac{\alpha-1}{2},

I(z)\displaystyle I(z) =z1+α01vα12e𝚒πzv𝑑v=z1+α01vα12cos(πzv)𝑑v𝚒z1+α01vα12sin(πzv)𝑑v\displaystyle=z^{1+\alpha}\int_{0}^{1}v^{\frac{\alpha-1}{2}}\mathrm{e}^{-\mathtt{i}\pi z\sqrt{v}}dv=z^{1+\alpha}\int_{0}^{1}v^{\frac{\alpha-1}{2}}\cos(\pi z\sqrt{v})dv-\mathtt{i}z^{1+\alpha}\int_{0}^{1}v^{\frac{\alpha-1}{2}}\sin(\pi z\sqrt{v})dv
=2z1+αH+F21(H+2;12,1+H+2;r)𝚒zH+4r1+H+F21(12+H+2;32,32+H+2;r)\displaystyle=\frac{2z^{1+\alpha}}{H_{+}}{}_{1}{F}_{2}\left(\frac{H_{+}}{2};\frac{1}{2},1+\frac{H_{+}}{2};-r\right)-\mathtt{i}z^{H_{+}}\frac{4\sqrt{r}}{1+H_{+}}{}_{1}{F}_{2}\left(\frac{1}{2}+\frac{H_{+}}{2};\frac{3}{2},\frac{3}{2}+\frac{H_{+}}{2};-r\right)
=zH+h1F21(h1;12,1+h1;π2z24)𝚒πz1+H+h2F21(h2;32,1+h2;π2z24),\displaystyle=\frac{z^{H_{+}}}{h_{1}}{}_{1}{F}_{2}\left(h_{1};\frac{1}{2},1+h_{1};-\frac{\pi^{2}z^{2}}{4}\right)-\mathtt{i}\frac{\pi z^{1+H_{+}}}{h_{2}}{}_{1}{F}_{2}\left(h_{2};\frac{3}{2},1+h_{2};-\frac{\pi^{2}z^{2}}{4}\right),

since α=H12=H+1\alpha=H-\frac{1}{2}=H_{+}-1, h1=H+2h_{1}=\frac{H_{+}}{2} and h2=12+h1h_{2}=\frac{1}{2}+h_{1}. Plugging these into (A.3), we obtain

0T(tu)αe𝚒πu𝑑u\displaystyle\int_{0}^{T}(t-u)^{\alpha}\mathrm{e}^{\mathtt{i}\pi u}du =e𝚒πt2[I(t)I(τ)]\displaystyle=\frac{\mathrm{e}^{\mathtt{i}\pi t}}{2}\left[I(t)-I(\tau)\right]
=e𝚒πt2[zH+h1F21(h1;12,1+h1;π2z24)𝚒πz1+H+h2F21(h2;32,1+h2;π2z24)]z=t\displaystyle=\frac{\mathrm{e}^{\mathtt{i}\pi t}}{2}\Bigg{[}\frac{z^{H_{+}}}{h_{1}}{}_{1}{F}_{2}\left(h_{1};\frac{1}{2},1+h_{1};-\frac{\pi^{2}z^{2}}{4}\right)-\mathtt{i}\frac{\pi z^{1+H_{+}}}{h_{2}}{}_{1}{F}_{2}\left(h_{2};\frac{3}{2},1+h_{2};-\frac{\pi^{2}z^{2}}{4}\right)\Bigg{]}_{z=t}
e𝚒πt2[zH+h1F21(h1;12,1+h1;π2z24)𝚒πz1+H+h2F21(h2;32,1+h2;π2z24)]z=τ\displaystyle\quad-\frac{\mathrm{e}^{\mathtt{i}\pi t}}{2}\left[\frac{z^{H_{+}}}{h_{1}}{}_{1}{F}_{2}\left(h_{1};\frac{1}{2},1+h_{1};-\frac{\pi^{2}z^{2}}{4}\right)-\mathtt{i}\frac{\pi z^{1+H_{+}}}{h_{2}}{}_{1}{F}_{2}\left(h_{2};\frac{3}{2},1+h_{2};-\frac{\pi^{2}z^{2}}{4}\right)\right]_{z=\tau}
=e𝚒πt2h1[(t)H+F21(h1;12,1+h1;π2t24)(τ)H+F21(h1;12,1+h1;π2τ24)]\displaystyle=\frac{\mathrm{e}^{\mathtt{i}\pi t}}{2h_{1}}\left[(t)^{H_{+}}{}_{1}{F}_{2}\left(h_{1};\frac{1}{2},1+h_{1};-\frac{\pi^{2}t^{2}}{4}\right)-(\tau)^{H_{+}}{}_{1}{F}_{2}\left(h_{1};\frac{1}{2},1+h_{1};-\frac{\pi^{2}\tau^{2}}{4}\right)\right]
𝚒πe𝚒πt2h2[(t)1+H+F21(h2;32,1+h2;π2t24)(τ)1+H+F21(h2;32,1+h2;π2τ24)]\displaystyle\quad-\mathtt{i}\frac{\pi\mathrm{e}^{\mathtt{i}\pi t}}{2h_{2}}\left[(t)^{1+H_{+}}{}_{1}{F}_{2}\left(h_{2};\frac{3}{2},1+h_{2};-\frac{\pi^{2}t^{2}}{4}\right)-(\tau)^{1+H_{+}}{}_{1}{F}_{2}\left(h_{2};\frac{3}{2},1+h_{2};-\frac{\pi^{2}\tau^{2}}{4}\right)\right]
=e𝚒πt[ζ12(t,h1)ζ12(τ,h1)𝚒π(ζ32(t,h2)ζ32(τ,h2))],\displaystyle={\mathrm{e}^{\mathtt{i}\pi t}}\left[\zeta_{\frac{1}{2}}(t,h_{1})-\zeta_{\frac{1}{2}}(\tau,h_{1})-\mathtt{i}\pi\left(\zeta_{\frac{3}{2}}(t,h_{2})-\zeta_{\frac{3}{2}}(\tau,h_{2})\right)\right],

where χ(z):=14π2z2\chi(z):=-\frac{1}{4}\pi^{2}z^{2} and ζ12\zeta_{\frac{1}{2}} and ζ32\zeta_{\frac{3}{2}} as defined in the lemma.

A.4. Proof of Lemma 4.4

We first prove (A). For each nn\in\mathbb{N} and all t[T,1]t\in[T,1], recall that

𝒦HT[ψn](t)=20T(tu)H12cos(uλn)𝑑u=2tTtvH12cos(tvλn)𝑑v,\mathcal{K}_{H}^{T}[\psi_{n}](t)=\sqrt{2}\int_{0}^{T}(t-u)^{H-\frac{1}{2}}\cos\left(\frac{u}{\sqrt{\lambda_{n}}}\right)du=\sqrt{2}\int_{t-T}^{t}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv,

with the change of variables v=tuv=t-u. Assume Ts<t1{T}\leq s<t\leq 1. Two situations are possible:

  • If 0sT<tTs<t10\leq s-T<{t-T}\leq s<t\leq 1, we have

    |𝒦HT[ψn](t)𝒦HT[ψn](s)|\displaystyle\left|\mathcal{K}_{H}^{T}[\psi_{n}](t)-\mathcal{K}_{H}^{T}[\psi_{n}](s)\right| =\displaystyle= 2|tTtvH12cos(tvλn)𝑑vsTsvH12cos(svλn)𝑑v|\displaystyle\sqrt{2}\left|\int_{t-T}^{t}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv-\int_{s-T}^{s}v^{H-\frac{1}{2}}\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)dv\right|
    \displaystyle\leq 2(|tTsvH12(cos(tvλn)cos(svλn))dv|\displaystyle\sqrt{2}\Bigg{(}\left|\int_{t-T}^{s}v^{H-\frac{1}{2}}\left(\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)-\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)\right)dv\right|
    +|stvH12cos(tvλn)dv|+|sTtTvH12cos(svλn)dv|)\displaystyle+\left|\int_{s}^{t}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv\right|+\left|\int_{s-T}^{t-T}v^{H-\frac{1}{2}}\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)dv\right|\Bigg{)}
    \displaystyle\leq 2(tTsvH12|cos(tvλn)cos(svλn)|dv\displaystyle\sqrt{2}\Bigg{(}\int_{t-T}^{s}v^{H-\frac{1}{2}}\left|\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)-\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)\right|dv
    +stvH12dv+sTtTvH12dv)\displaystyle+\int_{s}^{t}v^{H-\frac{1}{2}}dv+\int_{s-T}^{t-T}v^{H-\frac{1}{2}}dv\Bigg{)}
    \displaystyle\leq 2(tTsvH12|tsλn|𝑑v+K|ts|H+12+K|ts|H+12)\displaystyle\sqrt{2}\Bigg{(}\int_{t-T}^{s}v^{H-\frac{1}{2}}\left|\frac{t-s}{\sqrt{\lambda_{n}}}\right|dv+K|t-s|^{H+\frac{1}{2}}+K|t-s|^{H+\frac{1}{2}}\Bigg{)}
    \displaystyle\leq 2(|ts|λntTsvH12𝑑v+2K|ts|H+12)\displaystyle\sqrt{2}\Bigg{(}\frac{|t-s|}{\sqrt{\lambda_{n}}}\int_{t-T}^{s}v^{H-\frac{1}{2}}dv+2K|t-s|^{H+\frac{1}{2}}\Bigg{)}
    \displaystyle\leq 2(|ts|λn()H12L1[0,1]+2K|ts|H+12)C~1T|ts|H+12,\displaystyle\sqrt{2}\Bigg{(}\frac{|t-s|}{\sqrt{\lambda_{n}}}\|(\cdot)^{H-\frac{1}{2}}\|_{L^{1}[0,1]}+2K|t-s|^{H+\frac{1}{2}}\Bigg{)}\leq\widetilde{C}_{1}^{T}|t-s|^{H+\frac{1}{2}},

    with C~1T=max{22K,2λn()H12L1[0,1]}=max{22K,2(2n1)π2()H12L1[0,1]}\widetilde{C}_{1}^{T}=\max\left\{2\sqrt{2}K,\sqrt{\frac{{2}}{\lambda_{n}}}\|(\cdot)^{H-\frac{1}{2}}\|_{L^{1}[0,1]}\right\}=\max\left\{2\sqrt{2}K,\frac{\sqrt{2}(2n-1)\pi}{2}\|(\cdot)^{H-\frac{1}{2}}\|_{L^{1}[0,1]}\right\}, since cos()\cos(\cdot) is Lipschitz on any compact and 0vH12𝑑v\int_{0}^{\cdot}v^{H-\frac{1}{2}}dv is (H+12)\left(H+\frac{1}{2}\right)-Hölder continuous.

  • If 0sTstTt10\leq s-T\leq s\leq t-T\leq t\leq 1,

    |𝒦HT[ψn](t)𝒦HT[ψn](s)|\displaystyle\left|\mathcal{K}_{H}^{T}[\psi_{n}](t)-\mathcal{K}_{H}^{T}[\psi_{n}](s)\right| =\displaystyle= 2|tTtvH12cos(tvλn)𝑑vsTsvH12cos(svλn)𝑑v|\displaystyle\sqrt{2}\left|\int_{t-T}^{t}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv-\int_{s-T}^{s}v^{H-\frac{1}{2}}\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)dv\right|
    =\displaystyle= 2|tTtvH12cos(tvλn)𝑑vsTsvH12cos(svλn)𝑑v\displaystyle\sqrt{2}\Bigg{|}\int_{t-T}^{t}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv-\int_{s-T}^{s}v^{H-\frac{1}{2}}\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)dv
    +stTvH12cos(tvλn)𝑑vstTvH12cos(tvλn)𝑑v\displaystyle+\int_{s}^{t-T}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv-\int_{s}^{t-T}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv
    +stTvH12cos(svλn)dvstTvH12cos(svλn)dv|\displaystyle+\int_{s}^{t-T}v^{H-\frac{1}{2}}\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)dv-\int_{s}^{t-T}v^{H-\frac{1}{2}}\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)dv\Bigg{|}
    \displaystyle\leq 2(|stTvH12(cos(tvλn)cos(svλn))dv|\displaystyle\sqrt{2}\Bigg{(}\left|\int_{s}^{t-T}v^{H-\frac{1}{2}}\left(\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)-\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)\right)dv\right|
    +|stvH12cos(tvλn)dv|+|sTtTvH12cos(svλn)dv|)\displaystyle+\left|\int_{s}^{t}v^{H-\frac{1}{2}}\cos\left(\frac{t-v}{\sqrt{\lambda_{n}}}\right)dv\right|+\left|\int_{s-T}^{t-T}v^{H-\frac{1}{2}}\cos\left(\frac{s-v}{\sqrt{\lambda_{n}}}\right)dv\right|\Bigg{)}
    \displaystyle\leq C~1T|ts|H+12,\displaystyle\dots\leq\widetilde{C}_{1}^{T}|t-s|^{H+\frac{1}{2}},

    where the dots correspond to the same computations as in the previous case and leads to the same estimation with the same constant C~1T\widetilde{C}_{1}^{T}.

This proves (A).

To prove (B), recall that, for T[0,1]T\in[0,1] and nn\in\mathbb{N}, the function 𝒦HT[ψn]:[T,1]\mathcal{K}^{T}_{H}[\psi_{n}]:[T,1]\to\mathbb{R} reads

𝒦HT[ψn](t)\displaystyle\mathcal{K}^{T}_{H}[\psi_{n}](t) =20T(ts)H12cos((n12)πs)𝑑s\displaystyle=\sqrt{2}\int_{0}^{T}(t-s)^{H-\frac{1}{2}}\cos\left(\left(n-\frac{1}{2}\right)\pi s\right)ds
=2mH+120mT(mtu)H12cos(πu)du=:Φm(t).\displaystyle=\frac{\sqrt{2}}{m^{H+\frac{1}{2}}}\int_{0}^{mT}(mt-u)^{H-\frac{1}{2}}\cos\left(\pi u\right)du=:\Phi_{m}(t). (60)

with the change of variable u=(n12)s=:msu=(n-\frac{1}{2})s=:ms. Denote from now on ~:={m=n12,n}\widetilde{\mathbb{N}}:=\{m=n-\frac{1}{2},n\in\mathbb{N}\}. From (A.4), we deduce, for each m~m\in\widetilde{\mathbb{N}} and t[T,1]t\in[T,1],

mH+12Φm(t)=20mT(mtu)H12cos(πu)du=:2ϕm(t).m^{H+\frac{1}{2}}\Phi_{m}(t)=\sqrt{2}\int_{0}^{mT}(mt-u)^{H-\frac{1}{2}}\cos\left(\pi u\right)du=:\sqrt{2}\phi_{m}(t). (61)

To end the proof of (B), it therefore suffices to show that (ϕm(t))m~,t[T,1](\phi_{m}(t))_{m\in\widetilde{\mathbb{N}},t\in[T,1]} is uniformly bounded since, in that case we have

𝒦HT[ψn]\displaystyle\|\mathcal{K}_{H}^{T}[\psi_{n}]\|_{\infty} =\displaystyle= supt[T,1]|𝒦HT[ψn](t)|=supt[T,1]|Φn12(t)|=2(n12)H+12supt[T,1]|ϕn12(t)|\displaystyle\sup_{t\in[T,1]}|\mathcal{K}_{H}^{T}[\psi_{n}](t)|=\sup_{t\in[T,1]}|\Phi_{n-\frac{1}{2}}(t)|=\frac{\sqrt{2}}{(n-\frac{1}{2})^{H+\frac{1}{2}}}\sup_{t\in[T,1]}|\phi_{n-\frac{1}{2}}(t)|
\displaystyle\leq 2(n12)H+12supt[T,1],m~|ϕm(t)|2(n12)H+12CC2Tn(H+12),\displaystyle\frac{\sqrt{2}}{(n-\frac{1}{2})^{H+\frac{1}{2}}}\sup_{t\in[T,1],m\in\widetilde{\mathbb{N}}}|\phi_{m}(t)|\leq\frac{\sqrt{2}}{(n-\frac{1}{2})^{H+\frac{1}{2}}}C\leq C_{2}^{T}n^{-(H+\frac{1}{2})},

for some C2T>0C_{2}^{T}>0, proving (B). The following guarantees the uniform boundedness of ϕx\phi_{x} in (61).

Proposition A.1.

For any T[0,1]T\in[0,1], there exists C>0C>0 such that |ϕx(t)|C|\phi_{x}(t)|\leq C for all x0x\geq 0, t[T,1]t\in[T,1].

Proof.

For x>0x>0, we write

ϕx(t)=0xT(xtu)H12cos(πu)𝑑u={0xT(xtu)H12e𝚒πu𝑑u}.\phi_{x}(t)=\int_{0}^{xT}(xt-u)^{H-\frac{1}{2}}\cos\left(\pi u\right)du=\Re\left\{\int_{0}^{xT}(xt-u)^{H-\frac{1}{2}}\mathrm{e}^{\mathtt{i}\pi u}du\right\}.

Using the representation in Lemma 4.5, we are thus left to prove that the maps ζ12(,h1)\zeta_{\frac{1}{2}}(\cdot,h_{1}) and ζ32(,h2)\zeta_{\frac{3}{2}}(\cdot,h_{2}), defined in (27), are bounded on [0,)[0,\infty) by, say L12L_{\frac{1}{2}} and L32L_{\frac{3}{2}}. Indeed, in this case,

supx>0,t[T,1]|ϕx(t)|\displaystyle\sup_{x>0,t\in[T,1]}|\phi_{x}(t)| =\displaystyle= supx>0,t[T,1]|0xT(xtu)H12e𝚒πu𝑑u|\displaystyle\sup_{x>0,t\in[T,1]}\left|\int_{0}^{xT}(xt-u)^{H-\frac{1}{2}}\mathrm{e}^{\mathtt{i}\pi u}du\right|
\displaystyle\leq supx>0,t[T,1]|e𝚒πxt2[(ζ12(xt,h1)ζ12(x(tT),h1))𝚒π(ζ32(xt,h2)ζ32(x(tT),h2))]|\displaystyle\sup_{x>0,t\in[T,1]}\Bigg{|}\frac{\mathrm{e}^{\mathtt{i}\pi xt}}{2}\Bigg{[}\Big{(}\zeta_{\frac{1}{2}}(xt,h_{1})-\zeta_{\frac{1}{2}}(x(t-T),h_{1})\Big{)}-\mathtt{i}\pi\Big{(}\zeta_{\frac{3}{2}}(xt,h_{2})-\zeta_{\frac{3}{2}}(x(t-T),h_{2})\Big{)}\Bigg{]}\Bigg{|}
\displaystyle\leq 12supy,z[0,)|(ζ12(y,h1)ζ12(z,h1))𝚒π(ζ32(y,h2)ζ32(z,h2))|\displaystyle\frac{1}{2}\sup_{y,z\in[0,\infty)}\Bigg{|}\Big{(}\zeta_{\frac{1}{2}}(y,h_{1})-\zeta_{\frac{1}{2}}(z,h_{1})\Big{)}-\mathtt{i}\pi\Big{(}\zeta_{\frac{3}{2}}(y,h_{2})-\zeta_{\frac{3}{2}}(z,h_{2})\Big{)}\Bigg{|}
\displaystyle\leq π{supy[0,)|ζ12(y,h1)|+supy[0,)|ζ32(y,h2)|}L12+L32=C<+.\displaystyle\pi\left\{\sup_{y\in[0,\infty)}\left|\zeta_{\frac{1}{2}}(y,h_{1})\right|+\sup_{y\in[0,\infty)}\left|\zeta_{\frac{3}{2}}(y,h_{2})\right|\right\}\leq L_{\frac{1}{2}}+L_{\frac{3}{2}}=C<+\infty.

The maps ζ12(,h1)\zeta_{\frac{1}{2}}(\cdot,h_{1}) and ζ32(,h2)\zeta_{\frac{3}{2}}(\cdot,h_{2}) are both clearly continuous. Moreover, as zz tends to infinity ζk(z,h)\zeta_{k}(z,h) converges to a constant ckc_{k}, for (k,h)({12,32},{h1,h2})(k,h)\in(\{\frac{1}{2},\frac{3}{2}\},\{h_{1},h_{2}\}). The identities

F21(h;12,1+h;x)h=01cos(2xu)u1h𝑑uandF21(h;32,1+h;x)h=12x01sin(2xu)u3/2h𝑑u\frac{{}_{1}F_{2}\left(h;\frac{1}{2},1+h;-x\right)}{h}=\int_{0}^{1}\frac{\cos(2\sqrt{xu})}{u^{1-h}}du\quad\text{and}\quad\frac{{}_{1}F_{2}\left(h;\frac{3}{2},1+h;-x\right)}{h}=\frac{1}{2\sqrt{x}}\int_{0}^{1}\frac{\sin(2\sqrt{xu})}{u^{3/2-h}}du

hold (this can be checked with Wolfram Mathematica for example) and therefore,

ζ12(z,h1)\displaystyle\zeta_{\frac{1}{2}}(z,h_{1}) =\displaystyle= z2h12h1F21(h1;12,1+h1;π2z24)=z2h1201uh11cos(πzu)𝑑u\displaystyle\frac{z^{2h_{1}}}{2h_{1}}{}_{1}F_{2}\left(h_{1};\frac{1}{2},1+h_{1};-\frac{\pi^{2}z^{2}}{4}\right)=\frac{z^{2h_{1}}}{2}\int_{0}^{1}u^{h_{1}-1}\cos(\pi z\sqrt{u})du
=\displaystyle= z2h120πzx2(h11)(πz)2(h11)cos(x)2xπ2z2𝑑x=1π2h10πzx2h11cos(x)𝑑x,\displaystyle\frac{z^{2h_{1}}}{2}\int_{0}^{\pi z}\frac{x^{2(h_{1}-1)}}{(\pi z)^{2(h_{1}-1)}}\cos(x)\frac{2x}{\pi^{2}z^{2}}dx=\frac{1}{\pi^{2h_{1}}}\int_{0}^{\pi z}x^{2h_{1}-1}\cos(x)dx,

where, in the second line, we used the change of variables x=πzux=\pi z\sqrt{u}. In particular, as zz tends to infinity, this converges to π2h10+x2h11cos(x)dx=cos(πh1)π2h1Γ(2h1)=:c1/20.440433\pi^{-2h_{1}}\int_{0}^{+\infty}x^{2h_{1}-1}\cos(x)dx=\frac{\cos(\pi h_{1})}{\pi^{2h_{1}}}\Gamma(2h_{1})=:c_{1/2}\approx 0.440433.
Analogously, for k=32k=\frac{3}{2},

ζ32(z,h2)\displaystyle\zeta_{\frac{3}{2}}(z,h_{2}) =\displaystyle= z2h22h2F21(h2;32,1+h2;π2z24)=z2h22πz01uh23/2sin(πzu)𝑑u\displaystyle\frac{z^{2h_{2}}}{2h_{2}}{}_{1}F_{2}\left(h_{2};\frac{3}{2},1+h_{2};-\frac{\pi^{2}z^{2}}{4}\right)=\frac{z^{2h_{2}}}{2\pi z}\int_{0}^{1}u^{h_{2}-3/2}\sin(\pi z\sqrt{u})du
=\displaystyle= z2h212π0πzx2h23)(πz)2h23)sin(x)2xπ2z2𝑑x=1π2h20πzx2(h21)sin(x)𝑑x,\displaystyle\frac{z^{2h_{2}-1}}{2\pi}\int_{0}^{\pi z}\frac{x^{2h_{2}-3)}}{(\pi z)^{2h_{2}-3)}}\sin(x)\frac{2x}{\pi^{2}z^{2}}dx=\frac{1}{\pi^{2h_{2}}}\int_{0}^{\pi z}x^{2(h_{2}-1)}\sin(x)dx,

with the same change of variables as before. This converges to π2h20+x2h22sin(x)dx=cos(πh2)π2h2Γ(2h21)=:c3/20.193\pi^{-2h_{2}}\int_{0}^{+\infty}x^{2h_{2}-2}\sin(x)dx=\frac{-\cos(\pi h_{2})}{\pi^{2h_{2}}}\Gamma(2h_{2}-1)=:c_{3/2}\approx 0.193 as zz tends to infinity. For k>0k>0, ζk(z,h)=z2h(1+𝒪(z2))\zeta_{k}(z,h)=z^{2h}(1+\mathcal{O}(z^{2})) at zero. Since H(0,12)H\in(0,\frac{1}{2}), the two functions are continuous and bounded and the proposition follows. ∎

A.5. Proof of Theorem 4.11

We only provide the proof of (32) since, as already noticed, that of (33) follows immediately. Suppose that F:F:\mathbb{R}\to\mathbb{R} is Lipschitz continuous with constant MM. By Definitions (18) and (30), we have

|𝔼[F(VIXT)]𝔼[F(VIX^T𝐝)]|\displaystyle\Big{|}\mathbb{E}\left[F\left(\mathrm{VIX}_{T}\right)\right]-\mathbb{E}\left[F\left(\widehat{\mathrm{VIX}}_{T}^{\boldsymbol{\mathrm{d}}}\right)\right]\Big{|}
=\displaystyle= |𝔼[F(|1ΔTT+Δv0(t)exp{γZtT,Δ+γ22(0tTK(s)2ds0tK(s)2ds)}dt|12)]\displaystyle\Bigg{|}\mathbb{E}\left[F\left(\left|\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(t)\exp\left\{\gamma Z^{T,\Delta}_{t}+\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)\right\}dt\right|^{\frac{1}{2}}\right)\right]
𝔼[F(|1ΔTT+Δv0(t)exp{γZ^tT,Δ,𝐝+γ22(0tTK(s)2ds0tK(s)2ds)}dt|12)]|.\displaystyle-\mathbb{E}\left[F\left(\left|\frac{1}{\Delta}\int_{T}^{T+\Delta}v_{0}(t)\exp\left\{\gamma\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}}_{t}+\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)\right\}dt\right|^{\frac{1}{2}}\right)\right]\Bigg{|}.

For clarity, let Z:=ZT,ΔZ:=Z^{T,\Delta}, Z^:=Z^T,Δ,𝐝\widehat{Z}:=\widehat{Z}^{T,\Delta,\boldsymbol{\mathrm{d}}}, :=TT+Δh(t)eγZt𝑑t\mathfrak{H}:=\int_{T}^{T+\Delta}h(t)e^{\gamma Z_{t}}dt and ^:=TT+Δh(t)eγZ^t𝑑t\widehat{\mathfrak{H}}:=\int_{T}^{T+\Delta}h(t)e^{\gamma\widehat{Z}_{t}}dt, with

h(t):=v0(t)Δexp{γ22(0tTK(s)2𝑑s0tK(s)2𝑑s)},for t[T,T+Δ].h(t):=\frac{v_{0}(t)}{\Delta}\exp\left\{\frac{\gamma^{2}}{2}\left(\int_{0}^{t-T}K(s)^{2}ds-\int_{0}^{t}K(s)^{2}ds\right)\right\},\qquad\text{for }t\in[T,T+\Delta].

We can therefore write, using the Lipschitz property of FF (with constant MM) and Lemma B.3,

|𝔼[F(VIXT)]𝔼[F(VIX^T𝐝)]|\displaystyle\left|\mathbb{E}\left[F\left(\mathrm{VIX}_{T}\right)\right]-\mathbb{E}\left[F\left(\widehat{\mathrm{VIX}}_{T}^{\boldsymbol{\mathrm{d}}}\right)\right]\right| =|𝔼[F(12)]𝔼[F(^12)]|𝔼[|F(12)F(^12)|]\displaystyle=\left|\mathbb{E}\left[F\left(\mathfrak{H}^{\frac{1}{2}}\right)\right]-\mathbb{E}\left[F\left(\widehat{\mathfrak{H}}^{\frac{1}{2}}\right)\right]\right|\leq\mathbb{E}\left[\left|F\left(\mathfrak{H}^{\frac{1}{2}}\right)-F\left(\widehat{\mathfrak{H}}^{\frac{1}{2}}\right)\right|\right]
M𝔼[|12^12|]M𝔼[(1+1^)|^|]\displaystyle\leq M\mathbb{E}\left[\left|\mathfrak{H}^{\frac{1}{2}}-\widehat{\mathfrak{H}}^{\frac{1}{2}}\right|\right]\leq M\mathbb{E}\Bigg{[}\left(\frac{1}{\mathfrak{H}}+\frac{1}{\widehat{\mathfrak{H}}}\right)\left|\mathfrak{H}-\widehat{\mathfrak{H}}\right|\Bigg{]}
=:M𝔼[A|^|]M𝔼[ATT+Δh(t)|eγZteγZ^t|dt]\displaystyle=:M\mathbb{E}\left[A\left|\mathfrak{H}-\widehat{\mathfrak{H}}\right|\right]\leq M\mathbb{E}\Bigg{[}A\int_{T}^{T+\Delta}h(t)\left|e^{\gamma Z_{t}}-e^{\gamma\widehat{Z}_{t}}\right|dt\Bigg{]}
M𝔼[ATT+Δh(t)γ(eγZt+eγZ^t)|ZtZ^t|𝑑t].\displaystyle\leq M\mathbb{E}\Bigg{[}A\int_{T}^{T+\Delta}h(t)\gamma\left(e^{\gamma Z_{t}}+e^{\gamma\widehat{Z}_{t}}\right)\left|Z_{t}-\widehat{Z}_{t}\right|dt\Bigg{]}.

Now, an application of Hölder’s inequality yields

|𝔼[F(VIXT)]𝔼[F(VIX^T𝐝)]|\displaystyle\left|\mathbb{E}\left[F\left(\mathrm{VIX}_{T}\right)\right]-\mathbb{E}\left[F\left(\widehat{\mathrm{VIX}}_{T}^{\boldsymbol{\mathrm{d}}}\right)\right]\right| M𝔼[γA|TT+Δh(t)2(eγZt+eγZ^t)2𝑑t|12|TT+Δ|ZtZ^t|2𝑑t|12]\displaystyle\leq M\mathbb{E}\left[\gamma A\left|\int_{T}^{T+\Delta}h(t)^{2}\left(e^{\gamma Z_{t}}+e^{\gamma\widehat{Z}_{t}}\right)^{2}dt\right|^{\frac{1}{2}}\left|\int_{T}^{T+\Delta}\left|Z_{t}-\widehat{Z}_{t}\right|^{2}dt\right|^{\frac{1}{2}}\right] (62)
M𝔼[(γA)2TT+Δh(t)2(eγZt+eγZ^t)2𝑑t]12𝔼[TT+Δ|ZtZ^t|2𝑑t]12\displaystyle\leq M\mathbb{E}\left[(\gamma A)^{2}\int_{T}^{T+\Delta}h(t)^{2}\left(e^{\gamma Z_{t}}+e^{\gamma\widehat{Z}_{t}}\right)^{2}dt\right]^{\frac{1}{2}}\mathbb{E}\left[\int_{T}^{T+\Delta}\left|Z_{t}-\widehat{Z}_{t}\right|^{2}dt\right]^{\frac{1}{2}} (63)
=𝔎𝔼[TT+Δ|ZtZ^t|2𝑑t]12,\displaystyle=\mathfrak{K}\ \mathbb{E}\left[\int_{T}^{T+\Delta}\left|Z_{t}-\widehat{Z}_{t}\right|^{2}dt\right]^{\frac{1}{2}}, (64)

where 𝔎:=M𝔼[γ2A2TT+Δh(t)2(eγZt+eγZ^t)2𝑑t]12\mathfrak{K}:=M\mathbb{E}[\gamma^{2}A^{2}\int_{T}^{T+\Delta}h(t)^{2}(e^{\gamma Z_{t}}+e^{\gamma\widehat{Z}_{t}})^{2}dt]^{\frac{1}{2}}. It remains to show that 𝔎\mathfrak{K} is a strictly positive finite constant. This follows from the fact that {Zt}t[T,T+Δ]\{Z_{t}\}_{t\in[T,T+\Delta]} does not explode in finite time (and so does not its quantization Z^\widehat{Z} either). The identity (a+b)22(a2+b2)(a+b)^{2}\leq 2(a^{2}+b^{2}) and Hölder’s inequality imply

𝔎2\displaystyle\mathfrak{K}^{2} \displaystyle\leq 4M2γ2𝔼[(1+1^)TT+Δh(t)2(e2γZt+e2γZ^t)𝑑t]\displaystyle 4M^{2}\gamma^{2}\mathbb{E}\left[\left(\frac{1}{\mathfrak{H}}+\frac{1}{\widehat{\mathfrak{H}}}\right)\int_{T}^{T+\Delta}h(t)^{2}\left(e^{2\gamma Z_{t}}+e^{2\gamma\widehat{Z}_{t}}\right)dt\right]
\displaystyle\leq 4M2γ2𝔼[|1+1^|2]12𝔼[|TT+Δh(t)2(e2γZt+e2γZ^t)𝑑t|2]12\displaystyle 4M^{2}\gamma^{2}\mathbb{E}\left[\left|\frac{1}{\mathfrak{H}}+\frac{1}{\widehat{\mathfrak{H}}}\right|^{2}\right]^{\frac{1}{2}}\mathbb{E}\left[\left|\int_{T}^{T+\Delta}h(t)^{2}\left(e^{2\gamma Z_{t}}+e^{2\gamma\widehat{Z}_{t}}\right)dt\right|^{2}\right]^{\frac{1}{2}}
\displaystyle\leq 16M2γ2𝔼[12+1^2]12𝔼[|TT+Δh(t)2e2γZt𝑑t|2+|TT+Δh(t)2e2γZ^t𝑑t|2]12\displaystyle 16M^{2}\gamma^{2}\mathbb{E}\left[\frac{1}{\mathfrak{H}^{2}}+\frac{1}{\widehat{\mathfrak{H}}^{2}}\right]^{\frac{1}{2}}\mathbb{E}\left[\left|\int_{T}^{T+\Delta}h(t)^{2}e^{2\gamma Z_{t}}dt\right|^{2}+\left|\int_{T}^{T+\Delta}h(t)^{2}e^{2\gamma\widehat{Z}_{t}}dt\right|^{2}\right]^{\frac{1}{2}}
=\displaystyle= :16M2γ2(A1+A2)12(B1+B2)12.\displaystyle:16M^{2}\gamma^{2}(A_{1}+A_{2})^{\frac{1}{2}}(B_{1}+B_{2})^{\frac{1}{2}}.

We only need to show that A1,A2,B1A_{1},A_{2},B_{1} and B2B_{2} are finite. Since hh is a positive continuous function on the compact interval [T,T+Δ][T,T+\Delta], we have

\displaystyle\mathfrak{H} TT+Δinfs[T,T+Δ](h(s)eγZs)dtΔinfs[T,T+Δ]h(s)eγZs\displaystyle\geq\int_{T}^{T+\Delta}\inf_{s\in[T,T+\Delta]}\left(h(s)e^{\gamma Z_{s}}\right)dt\geq\Delta\inf_{s\in[T,T+\Delta]}h(s)e^{\gamma Z_{s}} (65)
Δinft[T,T+Δ]h(t)infs[T,T+Δ]eγZsΔh~exp{γinfs[T,T+Δ]Zs},\displaystyle\geq\Delta\inf_{t\in[T,T+\Delta]}h(t)\inf_{s\in[T,T+\Delta]}e^{\gamma Z_{s}}\geq\Delta\widetilde{h}\exp\left\{\gamma\inf_{s\in[T,T+\Delta]}Z_{s}\right\}, (66)

with h~:=inft[T,T+Δ]h(t)>0\widetilde{h}:=\inf_{t\in[T,T+\Delta]}h(t)>0. The inequality (65) implies

A1\displaystyle A_{1} =\displaystyle= 𝔼[2]𝔼[exp{2γinfs[T,T+Δ]Zs}]Δ2h~2=𝔼[exp{2γsups[T,T+Δ](Zs)}]Δ2h~2\displaystyle\mathbb{E}\left[\mathfrak{H}^{-2}\right]\leq\frac{\mathbb{E}\left[\exp\left\{-2\gamma\inf_{s\in[T,T+\Delta]}Z_{s}\right\}\right]}{\Delta^{2}\widetilde{h}^{2}}=\frac{\mathbb{E}\left[\exp\left\{2\gamma\sup_{s\in[T,T+\Delta]}(-Z_{s})\right\}\right]}{\Delta^{2}\widetilde{h}^{2}}
=\displaystyle= 1Δ2h~2𝔼[exp{2γsups[T,T+Δ]Zs}],\displaystyle\frac{1}{\Delta^{2}\widetilde{h}^{2}}\mathbb{E}\left[\exp\left\{2\gamma\sup_{s\in[T,T+\Delta]}Z_{s}\right\}\right],

since Z-Z and ZZ have the same law. The process Z=(Zt)t[T,T+Δ]Z=(Z_{t})_{t\in[T,T+\Delta]} is a continuous centered Gaussian process defined on a compact set. Thus, by Theorem 1.5.4 in [2], it is almost surely bounded there. Furthermore, exploiting Lemma B.4 and Borel-TIS inequality [2, Theorem 2.1.1], we have

𝔼[e2γsups[T,T+Δ]Zs]=:𝔼[e2γZ]=0+(e2γZ>u)du=0+(Z>log(u)2γ)du\displaystyle\mathbb{E}\left[e^{2\gamma\sup_{s\in[T,T+\Delta]}Z_{s}}\right]=:\mathbb{E}\left[e^{2\gamma\|Z\|}\right]=\int_{0}^{+\infty}\mathbb{P}\left(e^{2\gamma\|Z\|}>u\right)du=\int_{0}^{+\infty}\mathbb{P}\left(\|Z\|>\frac{\log(u)}{2\gamma}\right)du
=0e2γ𝔼[Z]𝑑u+e2γ𝔼[V]+(Z>log(u)2γ)𝑑u=e2γ𝔼[Z]+e2γ𝔼[V]+e12(12γlog(u)𝔼[Z]σT)2𝑑u\displaystyle=\int_{0}^{e^{2\gamma\mathbb{E}[\|Z\|]}}du+\int_{e^{2\gamma\mathbb{E}[\|V\|]}}^{+\infty}\mathbb{P}\left(\|Z\|>\frac{\log(u)}{2\gamma}\right)du=e^{2\gamma\mathbb{E}[\|Z\|]}+\int_{e^{2\gamma\mathbb{E}[\|V\|]}}^{+\infty}e^{-\frac{1}{2}\left(\frac{\frac{1}{2\gamma}\log(u)-\mathbb{E}[\|Z\|]}{\sigma_{T}}\right)^{2}}du
e2γ𝔼[Z]+0+e12(12γlog(u)𝔼[Z]σT)2𝑑u,\displaystyle\leq e^{2\gamma\mathbb{E}[\|Z\|]}+\int_{0}^{+\infty}e^{-\frac{1}{2}\left(\frac{\frac{1}{2\gamma}\log(u)-\mathbb{E}[\|Z\|]}{\sigma_{T}}\right)^{2}}du, (67)

with Z:=sups[T,T+Δ]Zs\|Z\|:=\sup_{s\in[T,T+\Delta]}Z_{s} and σT2:=supt[T,T+Δ]𝔼[Zt2]\sigma_{T}^{2}:=\sup_{t\in[T,T+\Delta]}\mathbb{E}[Z_{t}^{2}]. The change of variable log(u)2γ=v\frac{\log(u)}{2\gamma}=v in the last term in (A.5) yields

0+e12(12γlog(u)𝔼[Z]σT)2𝑑u=2γe12(v𝔼[Z]σT)2e2γv𝑑v=2π2γ𝔼[e2γY],\int_{0}^{+\infty}e^{-\frac{1}{2}\left(\frac{\frac{1}{2\gamma}\log(u)-\mathbb{E}[\|Z\|]}{\sigma_{T}}\right)^{2}}du=2\gamma\int_{\mathbb{R}}e^{-\frac{1}{2}\left(\frac{v-\mathbb{E}[\|Z\|]}{\sigma_{T}}\right)^{2}}e^{2\gamma v}dv=\sqrt{2\pi}2\gamma\mathbb{E}[e^{2\gamma Y}],

since Y𝒩(𝔼[Z],σT)Y\sim\mathcal{N}(\mathbb{E}[\|Z\|],\sigma_{T}), and hence A1A_{1} is finite. Now, notice that, in analogy to the last line of the proof of Proposition 3.12, for any t[T,T+Δ]t\in[T,T+\Delta], we have

𝔼[Zt|(Z^s)s[T,T+Δ]]\displaystyle\mathbb{E}\left[Z_{t}\Big{|}(\widehat{Z}_{s})_{s\in[T,T+\Delta]}\right] =𝔼[𝔼[Zt|{ξ^nd(n)}n=1,,m]|(Z^s)s[T,T+Δ]]=𝔼[Z^t|(Z^s)s[T,T+Δ]]=Z^t,\displaystyle=\mathbb{E}\left[\mathbb{E}\left[Z_{t}\Big{|}\{\widehat{\xi}_{n}^{d(n)}\}_{n=1,\dots,m}\right]\Big{|}(\widehat{Z}_{s})_{s\in[T,T+\Delta]}\right]=\mathbb{E}\left[\widehat{Z}_{t}\Big{|}(\widehat{Z}_{s})_{s\in[T,T+\Delta]}\right]=\widehat{Z}_{t}, (68)

since the sigma-algebra generated by (Z^s)s[T,T+Δ](\widehat{Z}_{s})_{s\in[T,T+\Delta]} is included in the sigma-algebra generated by {ξ^nd(n)}n=1,,m\{\widehat{\xi}_{n}^{d(n)}\}_{n=1,\dots,m}. Now, exploiting, in sequence, (68), the conditional version of supt[T1,T2]𝔼[ft]𝔼[supt[T1,T2]ft]\sup_{t\in[T_{1},T_{2}]}\mathbb{E}[f_{t}]\leq\mathbb{E}[\sup_{t\in[T_{1},T_{2}]}f_{t}], conditional Jensen’s inequality together with the convexity of xeγxx\mapsto e^{\gamma x}, for γ>0\gamma>0 and the tower property, we obtain

𝔼[exp{γsupt[T,T+Δ]Z^t}]\displaystyle\mathbb{E}\left[\exp\left\{\gamma\sup_{t\in[T,T+\Delta]}\widehat{Z}_{t}\right\}\right] =𝔼[exp{γsupt[T,T+Δ]𝔼[Zt|(Z^s)s[T,T+Δ]]}]\displaystyle=\mathbb{E}\left[\exp\left\{\gamma\sup_{t\in[T,T+\Delta]}\mathbb{E}\left[Z_{t}\Big{|}(\widehat{Z}_{s})_{s\in[T,T+\Delta]}\right]\right\}\right]
𝔼[exp{γ𝔼[supt[T,T+Δ]Zt|(Z^s)s[T,T+Δ]]}]\displaystyle\leq\mathbb{E}\left[\exp\left\{\gamma\mathbb{E}\left[\sup_{t\in[T,T+\Delta]}Z_{t}\Big{|}(\widehat{Z}_{s})_{s\in[T,T+\Delta]}\right]\right\}\right]
𝔼[𝔼[exp{γsupt[T,T+Δ]Zt}|(Z^s)s[T,T+Δ]]]\displaystyle\leq\mathbb{E}\left[\mathbb{E}\left[\exp\left\{\gamma\sup_{t\in[T,T+\Delta]}Z_{t}\right\}\Big{|}(\widehat{Z}_{s})_{s\in[T,T+\Delta]}\right]\right]
=𝔼[exp{γsupt[T,T+Δ]Zt}].\displaystyle=\mathbb{E}\left[\exp\left\{\gamma\sup_{t\in[T,T+\Delta]}Z_{t}\right\}\right]. (69)

Thus, we have

A2=𝔼[^2]1Δ2h~2𝔼[exp{γsupt[T,T+Δ]Z^t}]1Δ2h~2𝔼[exp{γsupt[T,T+Δ]Zt}],A_{2}=\mathbb{E}\left[\widehat{\mathfrak{H}}^{-2}\right]\leq\frac{1}{\Delta^{2}\widetilde{h}^{2}}\mathbb{E}\left[\exp\left\{\gamma\sup_{t\in[T,T+\Delta]}\widehat{Z}_{t}\right\}\right]\leq\frac{1}{\Delta^{2}\widetilde{h}^{2}}\mathbb{E}\left[\exp\left\{\gamma\sup_{t\in[T,T+\Delta]}Z_{t}\right\}\right],

which is finite because of the proof of the finiteness of A1A_{1}, above.

Exploiting Fubini’s theorem we rewrite B1B_{1} as

B1=𝔼[(TT+Δh(t)2e2γZt𝑑t)2]=TT+ΔTT+Δh(t)2h(s)2𝔼[e2γ(Zt+Zs)]𝑑t𝑑s.B_{1}=\mathbb{E}\left[\left(\int_{T}^{T+\Delta}h(t)^{2}e^{2\gamma Z_{t}}dt\right)^{2}\right]=\int_{T}^{T+\Delta}\int_{T}^{T+\Delta}h(t)^{2}h(s)^{2}\mathbb{E}\left[e^{2\gamma(Z_{t}+Z_{s})}\right]dtds. (70)

Since (Zt)t[T,T+Δ](Z_{t})_{t\in[T,T+\Delta]} is centered Gaussian with covariance 𝔼[ZtZs]=0TK(tu)K(su)𝑑u\mathbb{E}[Z_{t}Z_{s}]=\int_{0}^{T}K(t-u)K(s-u)du, then (Zt+Zs)𝒩(0,g(t,s))(Z_{t}+Z_{s})\sim\mathcal{N}(0,g(t,s)), with g(t,s):=𝔼[(Zt+Zs)2]=0T(K(tu)+K(su))2𝑑ug(t,s):=\mathbb{E}[(Z_{t}+Z_{s})^{2}]=\int_{0}^{T}(K(t-u)+K(s-u))^{2}du and therefore

B1=TT+ΔTT+Δh(t)2h(s)2e2γ2g(t,s)𝑑t𝑑sB_{1}=\int_{T}^{T+\Delta}\int_{T}^{T+\Delta}h(t)^{2}h(s)^{2}e^{2\gamma^{2}g(t,s)}dtds (71)

is finite since both hh and gg are continuous on compact intervals. Finally, for B2B_{2} we have

B2\displaystyle B_{2} =\displaystyle= 𝔼[(TT+Δh(t)2e2γZ^t𝑑t)2]=TT+ΔTT+Δh(t)2h(s)2𝔼[e2γ(Z^t+Z^s)]𝑑t𝑑s\displaystyle\mathbb{E}\left[\left(\int_{T}^{T+\Delta}h(t)^{2}e^{2\gamma\widehat{Z}_{t}}dt\right)^{2}\right]=\int_{T}^{T+\Delta}\int_{T}^{T+\Delta}h(t)^{2}h(s)^{2}\mathbb{E}\left[e^{2\gamma(\widehat{Z}_{t}+\widehat{Z}_{s})}\right]dtds
\displaystyle\leq TT+ΔTT+Δh(t)2h(s)2𝔼[e2γ(Zt+Zs)]𝑑t𝑑s=B1,\displaystyle\int_{T}^{T+\Delta}\int_{T}^{T+\Delta}h(t)^{2}h(s)^{2}\mathbb{E}\left[e^{2\gamma({Z}_{t}+{Z}_{s})}\right]dtds=B_{1},

where we have used the fact that for all t,s[T,T+Δ]t,s\in[T,T+\Delta], (Z^t+Z^s)(\widehat{Z}_{t}+\widehat{Z}_{s}) is a stationary quantizer for (Zt+Zs)(Z_{t}+Z_{s}) and so 𝔼[e2γ(Z^t+Z^s)]𝔼[e2γ(Zt+Zs)]\mathbb{E}[e^{2\gamma(\widehat{Z}_{t}+\widehat{Z}_{s})}]\leq\mathbb{E}[e^{2\gamma({Z}_{t}+{Z}_{s})}] since f(x)=e2γxf(x)=e^{2\gamma x} is a convex function (see Remark 3.9 in Section 3.1). Therefore B2B_{2} is finite and the proof follows.

Appendix B Some useful results

We recall some important results used throughout the text. Straightforward proofs are omitted.

Proposition B.1.

For a Gaussian random variable Z𝒩(μ,σ)Z\sim\mathcal{N}(\mu,\sigma),

𝔼[|Zμ|p]={(p1)!!σp,if p is even,0,if p is odd.\mathbb{E}\left[|Z-\mu|^{p}\right]=\left\{\begin{array}[]{ll}(p-1)!!\sigma^{p},&\text{if }p\text{ is even},\\ 0,&\text{if }p\text{ is odd.}\end{array}\right.

We recall [39, Problem 8.5], correcting a small error, used in the proof of Proposition 3.6:

Lemma B.2.

Let m,Nm,N\in\mathbb{N} and p1,,pmp_{1},\dots,p_{m} positive real numbers. Then

inf{n=1mpnxn2:x1,,xm(0,),n=1mxnN}=mN2m(j=1mpj)1m,\inf\left\{\sum_{n=1}^{m}\frac{p_{n}}{x_{n}^{2}}:\quad x_{1},\dots,x_{m}\in(0,\infty),\quad\prod_{n=1}^{m}x_{n}\leq N\right\}=mN^{-\frac{2}{m}}\left(\prod_{j=1}^{m}p_{j}\right)^{\frac{1}{m}},

where the infimum is attained for xn=N1mpn12(j=1mpj)12mx_{n}=N^{\frac{1}{m}}p_{n}^{\frac{1}{2}}\left(\prod_{j=1}^{m}p_{j}\right)^{-\frac{1}{2m}}, for all n{1,,m}n\in\{1,\dots,m\}.

Proof.

The general arithmetic-geometric inequalities imply

1mn=1mpnxn2(n=1mpnxn2)1m=(n=1mpn)1m(n=1m1xn2)1m(n=1mpn)1mN2m,\frac{1}{m}\sum_{n=1}^{m}\frac{p_{n}}{x_{n}^{2}}\geq\left(\prod_{n=1}^{m}\frac{p_{n}}{x_{n}^{2}}\right)^{\frac{1}{m}}=\left(\prod_{n=1}^{m}p_{n}\right)^{\frac{1}{m}}\left(\prod_{n=1}^{m}\frac{1}{x_{n}^{2}}\right)^{\frac{1}{m}}\geq\left(\prod_{n=1}^{m}p_{n}\right)^{\frac{1}{m}}N^{-\frac{2}{m}},

since n=1mxnN\prod_{n=1}^{m}x_{n}\geq N by assumption. The right-hand side does not depend on x1,,xmx_{1},\dots,x_{m}, so

inf{n=1mpnxn2:x1,,xm(0,),n=1mxnN}m(n=1mpn)1mN2m.\inf\left\{\sum_{n=1}^{m}\frac{p_{n}}{x_{n}^{2}}:\quad x_{1},\dots,x_{m}\in(0,\infty),\quad\prod_{n=1}^{m}x_{n}\leq N\right\}\geq m\left(\prod_{n=1}^{m}p_{n}\right)^{\frac{1}{m}}N^{-\frac{2}{m}}.

Choosing x~n=N1mpn12(j=1mpj)12m\widetilde{x}_{n}=N^{\frac{1}{m}}p_{n}^{\frac{1}{2}}\left(\prod_{j=1}^{m}p_{j}\right)^{-\frac{1}{2m}}, for all n{1,,m}n\in\{1,\dots,m\}, we obtain

m(n=1mpnN2)1m=n=1mpnx~n2inf{n=1mpnxn2:x1,,xm(0,),n=1mxnN}m(n=1mpnN2)1m,\displaystyle m\left(\prod_{n=1}^{m}\frac{p_{n}}{N^{2}}\right)^{\frac{1}{m}}=\sum_{n=1}^{m}\frac{p_{n}}{\widetilde{x}_{n}^{2}}\geq\inf\left\{\sum_{n=1}^{m}\frac{p_{n}}{x_{n}^{2}}:x_{1},\dots,x_{m}\in(0,\infty),\prod_{n=1}^{m}x_{n}\leq N\right\}\geq m\left(\prod_{n=1}^{m}\frac{p_{n}}{N^{2}}\right)^{\frac{1}{m}},

which concludes the proof. ∎

Lemma B.3.

The following hold:

  1. (i)

    For any x,y>0x,y>0, |xy|(1x+1y)|xy||\sqrt{x}-\sqrt{y}|\leq\left(\frac{1}{\sqrt{x}}+\frac{1}{\sqrt{y}}\right)|x-y|.

  2. (ii)

    Set C>0C>0. For any x,yx,y\in\mathbb{R}, |eCxeCy|C(eCx+eCy)|xy||e^{Cx}-e^{Cy}|\leq C\left(e^{Cx}+e^{Cy}\right)|x-y|.

Lemma B.4.

For a positive random variable XX on (Ω,,)(\Omega,{\mathcal{F}},\mathbb{P}), 𝔼[X]=0+(X>u)𝑑u\mathbb{E}[X]=\int_{0}^{+\infty}\mathbb{P}(X>u)du.

Acknowledgements

The authors would like to thank Andrea Pallavicini and Emanuela Rosazza Gianin for fruitful discussions. The second author was supported by the Grant BIRD190200 “Term Structure Dynamics in Interest Rate and Energy Markets: Modelling and Numerics”.

References

  • [1] Abi Jaber, E. and El Euch, O. (2019): Multifactor approximation of rough volatility models, SIAM Journal on Financial Mathematics, 10(2), pp. 309-349.
  • [2] Adler, R.J. and Taylor, J.E. (2007): Random Fields and Geometry, Springer Monographs in Mathematics, New York, Springer-Verlag.
  • [3] Alòs, E.; León, J. A. and Vives J. (2007): On the short-time behavior of the implied volatility for jump-diffusion models with stochastic volatility, Finance and Stochastics, 11(4), pp. 571-589.
  • [4] Bayer, C.; Friz, P.K. and Gatheral, J. (2016): Pricing under rough volatility, Quantitative Finance, 16(6), pp. 887-904.
  • [5] Bayer, C.; Fukasawa, M. and Nakahara, S. (2022): On the weak convergence rate in the discretization of rough volatility models, ArXiV preprint, https://arxiv.org/abs/2203.02943.
  • [6] Bayer, C.; Hall, E.J. and Tempone, R. (2021): Weak error rates for option pricing under linear rough volatility, arXiv preprint, https://arxiv.org/abs/2009.01219.
  • [7] Bayer, C.; Hammouda, C.B. and Tempone, R. (2020): Hierarchical adaptive sparse grids and quasi Monte Carlo for option pricing under the rough Bergomi model, Quantitative Finance, 20(9), pp. 1457-1473.
  • [8] Bennedsen, M.; Lunde, A. and Pakkanen, M.S. (2017): Hybrid scheme for Brownian semistationary processes, Finance and Stochastics, 21, pp. 931-965.
  • [9] Bergomi; L. (2005); Smile dynamics II, Risk, pp. 67-73.
  • [10] Carr, P.P. and Madan, D.B. (2014): Joint modeling of VIX and SPX options at a single and common maturity with risk management applications, IIE Transactions, 46(11), pp. 1125-1131.
  • [11] Chen, W.; Langrené, N.; Loeper, G. and Zhu Q. (2021): Markovian approximation of the rough Bergomi model for Monte Carlo option pricing, Mathematics, 9(5), pp. 528.
  • [12] Chow, Y.S. and Teichner E. (1997): Probability Theory, Springer Texts in Statistics, New York, Springer-Verlag.
  • [13] Corlay, S. (2011): Quelques aspects de la quantification optimale, et applications en finance (in English, with French summary), PhD Thesis, Université Pierre et Marie Curie.
  • [14] Fukasawa, M. (2011): Asymptotic analysis for stochastic volatility: martingale expansion, Finance and Stochastics, 15(4), pp. 635-654.
  • [15] Fukasawa, M.; Takabatake, T. and Westphal, R. (2021): Is volatility rough?, Mathematical Finance, to appear.
  • [16] Fukasawa, M. (2021): Volatility has to be rough, Quantitative Finance, 21, pp. 1-8.
  • [17] Gassiat, P. (2022): Weak error rates of numerical schemes for rough volatility, arXiv preprint, https://arxiv.org/abs/2203.09298.
  • [18] Gatheral, J.; Jaisson, T. and Rosenbaum, M. (2018): Volatility is rough, Quantitative Finance, 18(6), pp. 933-949.
  • [19] Gatheral, J. (2008): Consistent modelling of SPX and VIX options, Presentation, Bachelier Congress, London.
  • [20] Gersho, A. and Gray, R.M. (1992): Vector Quantization and signal compression, New York, Kluwer Academic Publishers.
  • [21] Graf, S. and Luschgy., H. (2007): Foundations of quantization for probability distributions, Lecture Notes in Mathematics, 1730, Berlin Heidelberg, Springer.
  • [22] Horvath, B.; Jacquier, A. and Muguruza A. (2019): Functional central limit theorems for rough volatility, arxiv.org/abs/1711.03078.
  • [23] Horvath, B.; Jacquier, A. and Tankov P. (2020): Volatility options in rough volatility models, SIAM Journal on Financial Mathematics, 11(2).
  • [24] Huh, J.; Jeon, J. and Kim, J.H. (2018): A scaled version of the double-mean-reverting model for VIX derivatives, Mathematics and Financial Economics, 12(4), pp. 495-515.
  • [25] Jacquier, A.; Pakkanen, M. S. and Stone, H. (2018): Pathwise large deviations for the rough Bergomi model, Journal of Applied Probability, 55(4), pp. 1078-1092.
  • [26] Jacquier, A.; Martini, C. and Muguruza, A. (2018): On VIX Futures in the rough Bergomi model, Quantitative Finance, 18(1), pp. 45-61.
  • [27] Kallenberg, O. (2002): Foundations of Modern Probability, 2nd edition, Probability and Its Applications, New York, Springer-Verlag.
  • [28] Karp, D.B. (2015): Representations and inequalities for generalized Hypergeometric functions, J. Math. Sci. (N.Y.), 207, pp. 885-897.
  • [29] Kokholm, T. and Stisen, M. (2015): Joint pricing of VIX and SPX options with stochastic volatility and jump models, Journal of Risk Finance, 16(1), pp. 27-48.
  • [30] Luke, Y.L. (1969): The special functions and their approximations, Volume 1, Academic Press, New York and London.
  • [31] Luschgy, H. and Pagès, G. (2002): Functional quantization of Gaussian processes, Journal of Functional Analysis, 196(2), pp. 486-531.
  • [32] Luschgy, H. and Pagès, G. (2007): High-resolution product quantization for Gaussian processes under sup-norm distortion, Bernoulli, 13(3), pp. 653-671.
  • [33] McCrickerd, R. and Pakkanen, M.S. (2018): Turbocharging Monte Carlo pricing for the rough Bergomi model, Quantitative Finance, 18(11), pp. 1877-1886.
  • [34] Olver, F.W.J. (1997): Asymptotics and special functions, 2nd Edition, A.K. Peters / CRC Press.
  • [35] Pagès, G. (2007): Quadratic optimal functional quantization of stochastic processes and numerical applications, Monte Carlo and Quasi-Monte Carlo Methods 2006, Springer, Berlin Heidelberg, pp. 101-142.
  • [36] Pagès, G. and Printems, J. (2005): Functional quantization for numerics with an application to option pricing, Monte Carlo Methods and Applications, 11(4), pp. 407-446.
  • [37] Picard, J.(2011): Representation formulae for the fractional Brownian motion, Séminaire de Probabilités XLIII. Lecture Notes in Mathematics, 2006, Springer-Verlag, Berlin Heidelberg, pp. 3-70.
  • [38] Sheppard, W.F. (1897): On the calculation of the most probable values of frequency-constants, for data arranged according to equidistant division of a scale, Proc. Lond. Math. Soc. (3), 1(1), pp. 353-380.
  • [39] Steele, J. M. (2004): The Cauchy-Schwarz Master-Class, Cambridge University Press.