Decomposition formula for rough Volterra stochastic volatility models

The research presented in this article provides an alternative option pricing approach for a class of rough fractional stochastic volatility models. These models are increasingly popular between academics and practitioners due to their surprising consistency with financial markets. However, they bring several challenges alongside. Most noticeably, even simple non-linear financial derivatives as vanilla European options are typically priced by means of Monte-Carlo (MC) simulations which are more computationally demanding than similar MC schemes for standard stochastic volatility models. In this paper, we provide a proof of the prediction law for general Gaussian Volterra processes. The prediction law is then utilized to obtain an adapted projection of the future squared volatility -- a cornerstone of the proposed pricing approximation. Firstly, a decomposition formula for European option prices under general Volterra volatility models is introduced. Then we focus on particular models with rough fractional volatility and we derive an explicit semi-closed approximation formula. Numerical properties of the approximation for a popular model -- the rBergomi model -- are studied and we propose a hybrid calibration scheme which combines the approximation formula alongside MC simulations. This scheme can significantly speed up the calibration to financial markets as illustrated on a set of AAPL options.


Introduction
It is well known that the main issue of the Black-Scholes model lies in its assumptions about volatility of the modelled asset. Opposed to the model assumptions, the realized volatility time series tend to cluster, depend on the spot asset level and certainly they do not take a constant value within a reasonable time-frame (see e.g. Cont (2001)).
To deal with the aforementioned inconsistencies, stochastic volatility (SV) models have been proposed originally by Hull and White (1987) and later e.g. by Heston (1993). These models do not only assume that the asset price follows a specific stochastic process, but also that the instantaneous volatility of asset returns is of random nature as well. Especially, the latter approach by Heston became popular in the eyes of both practitioners and academics. Several modifications of this model have been proposed over the last 20 years: models with jump-diffusion dynamics (Bates 1996;Duffie, Pan, and Singleton 2000), with time-dependent parameters (Mikhailov and Nögel 2003;Elices 2008;Benhamou, Gobet, and Miri 2010), with fractionally scaled volatility (Comte and Renault 1998;Alòs, León, and Vives 2007;El Euch and Rosenbaum 2019) and models with several aspects combined Baustian, Mrázek, Pospíšil, and Sobotka 2017).
The original pricing approach of Heston (1993) was several times revisited, e.g by Lewis (2000), Attari (2004) and Albrecher, Mayer, Schoutens, and Tistaert (2007) with focus on semi-closed form Fourier transform solutions by Kahl and Jäckel (2006), Alfonsi (2010) with respect to Monte-Carlo simulation techniques and last but not least by Alòs (2012) who introduced an analytical approach to option pricing approximation. This approach improves on the techniques introduced by Hull and White (1987) and shows how an adaptive projection of future volatility can be used to price European options under Heston (1993) model. Many other papers generalized this idea, see e.g. Alòs, de Santiago, and Vives (2015); Merino and Vives (2015), or recently also a paper by Merino, Pospíšil, Sobotka, and Vives (2018). In this article, we revisit the latter approach and we come up with the approximation technique for SV models with volatility process driven by the fractional Brownian motion. This includes exponential rough fractional volatility models introduced by Rosenbaum (2014, 2018).
Although many SV models have been proposed since the original Hull and White (1987) model, it seems that none of them can be considered as the universal best market practice approach. Several models might perform well for calibration to complex volatility surfaces, but can suffer from over-fitting or they might not be robust in the sense described by Pospíšil, Sobotka, and Ziegler (2018). Also a model with a good fit to implied volatility surface might not be in-line with the observed time-series properties.
Pioneers of the fractional SV models - Comte and Renault (1998), see also Comte, Coutin, and Renault (2012) -assumed the so-called Hurst parameter 1 ranged within H ∈ (1/2, 1) which implies that the spot variance evolution is represented by a persistent process, i.e. it would have a long-memory property. In Alòs, León, and Vives (2007), a mean reverting fractional stochastic volatility model with H ∈ (0, 1) was presented. Gatheral, Jaisson, and Rosenbaum (2018) and Bayer, Friz, and Gatheral (2016) came up with a more detailed analysis of rough fractional volatility models which should be consistent with market option prices (Bayer, Friz, and Gatheral 2016), with realized volatility time series and also they could provide superior volatility prediction results to several other models (Bennedsen, Lunde, and Pakkanen 2017). An approach considering a two factor fractional volatility model, combining a rough term (H < 1 2 ) and a persistent one (H > 1 2 ), was presented in Funahashi and Kijima (2017). Recently, also an approximation for target-volatility options under log-normal fractional SABR model was studied by Alòs, Chatterjee, Tudor, and Wang (2019) who use the Malliavin caluculus techniques to derive the decomposition formula.
A typical problem of rough fractional models lies in their tractability -only cumbersome simulation techniques seem to be available for vanilla European options at the time of writing this article. This serves as a motivation to develop a pricing approximation and we introduce an approach based on the works of Alòs (2012) and Merino and Vives (2015).
The structure of the paper is the following. Section 2 is devoted to preliminaries. In Section 3 we present a generic decomposition formula of the vanilla European option fair value. In Setion 4 we present Volterra volatility models. In particular in Section 4.1 we prove the prediction law for general Gaussian Volterra processes, in Section 4.2 we obtain the decomposition formula for the exponential Volterra volatility model including the error bound for the approximation formula and in Section 4.3 we provide particular results for exponential fracional volatility models (including the standard Brownian motion case). In Section 5, we provide a numerical comparison of approximated prices against Monte Carlo simulations. All obtained results are concluded in Section 6.

Preliminaries and notation
Let S = (S t , t ∈ [0, T ]) be a strictly positive asset price process under a market chosen risk neutral probability measure P that follows the model: where S 0 is the current price, r ≥ 0 is the interest rate, W t andW t are independent standard Wiener processes defined on a probability space (Ω, F, P) and ρ ∈ (−1, 1). In the following, we will denote by F W and FW the filtrations generated by W andW respectively. Moreover, we define F := F W ∨ FW . The volatility process σ t is a square-integrable process assumed to be adapted to the filtration generated by W and its trajectories are assumed to be a.s. càdlàg and stricly positive a.e.. Note that ρ is the correlation between the price and the volatility processes. Without any loss of generality, it will be convenient in the following sections to make the change of variable X t = log S t , t ∈ [0, T ], and write Recall that Z := ρW + 1 − ρ 2W is a standard Wiener process.
The following notation will be used throughout the paper: • We will denote by BS(t, x, y) the price of a plain vanilla European call option under the classical Black-Scholes model with constant volatility y, current log stock price x, time to maturity τ = T − t, strike price K and interest rate r. In this case, where Φ(·) denotes the cumulative distribution function of the standard normal law and • In our setting, the call option price is given by where E t is the conditional expectation respect to the σ−algebra F t .
• Recall that from the Feynman-Kac formula for the model (2), the operator satisfies L y BS(t, x, y) = 0.

A generic decomposition formula
In this section, we provide an insight on a generic decomposition formula based on the work of Alòs (2012), Merino and Vives (2015) and Merino, Pospíšil, Sobotka, and Vives (2018). In particular, we recover the results for a generic stochastic volatility model presented in Merino, Pospíšil, Sobotka, and Vives (2018).
It is well known that if the stochastic volatility process is independent of the price process, the pricing formula of a vanilla European call option is given by whereσ 2 t is the so called average future variance that is defined bȳ Naturally,σ t is called the average future volatility. We consider the adapted projection of the future variance and the average future variance as to obtain a decomposition of V t in terms of v t . This idea switches an anticipative problem related to the anticipative processσ t into a non-anticipative one related to the adapted process v t . Taking into account M defined in (4), we can write In this paper, we will utilize the following lemma which is proved in Alòs (2012), p. 406.
where a 2 t is defined by (5). Remark 3.2. It is easy to see that the previous Lemma holds for put options and for several other non-path dependent options as well (e.g. Gap options). Now we use a generic decomposition formula proved in Merino, Pospíšil, Sobotka, and Vives (2018).
Theorem 3.3 (Generic decomposition formula). Let B t be a continuous semimartingale with respect to the filtration F W , let A(t, x, y) be a C 1,2,2 ([0, T ] × [0, ∞) × [0, ∞)) function and let v 2 t and M t be defined as above. Then we are able to formulate the expectation of e −rT A(T, X T , v 2 T )B T in the following way: Proof. See the proof of Theorem 3.1 in Merino, Pospíšil, Sobotka, and Vives (2018).
Using the previous decomposition formula, we find that Corollary 3.4 (BS decompostion formula). Under assumptions of the Theorem 3.3, we can obtain a decomposition of European option price V t as: Proof. Using Theorem 3.3 with A(t, X t , v 2 t ) = BS(t, X t , v t ) and B ≡ 1, the proof follows in a straightforward way.
The terms (I) and (II) are not easy to evaluate. Therefore, it becomes important to find simpler approximations to (I) and (II) and estimate the error terms. In order to find these approximations, we are going to apply Theorem 3.3 to find a decomposition formula for the terms (I) and (II). Using a decomposition of the term (I) can be found, and using After that process, we can approximate the price of a call option by where t denotes error terms. Terms of t under a general setting for σ t are provided in Appendix A. We note that the error term will depend on the assumed volatility dynamics.

General Volterra volatility model
In this section, we apply the generic decomposition formula to model (2) with general Volterra volatility process defined as where g : [0, +∞)×R → [0, +∞) is a deterministic function such that σ t belongs to L 1 (Ω×[0, +∞)) and Y = (Y t , t ≥ 0) is the Gaussian Volterra process where K(t, s) is a kernel such that for all t > 0 and denote the autocovariance function of process Y t and be the variance function (i.e. the second moment). Extending the Theorem 3.1 in Sottinen and Viitasaari (2017) enables us to rephrase the adapted projection of the future squared volatility.
Theorem 4.1 (Prediction law for Gaussian Volterra processes). Let (Y t , t ≥ 0) be the Gaussian Volterra process (7) satisfying assumptions (A1) and (A2). Then, the conditional process and deterministic covariance function In the upcoming sections, we will denoter(u|t) :=r (u, u|t). Under the general volatility process (6), we have In the upcoming lemma, we express the conditional expectation of the future squared volatility in terms of the mean functionm t (u).
Let us denote: Lemma 4.2 (Auxiliary terms in the decomposition formula for the general volatility model). Let Proof. Let 0 ≤ t ≤ u and Theorem 4.1 implies that is a Gaussian density function with stochastic meanm t (u) and deterministic variancer(u|t). To calculate the quadratic variation, we note that where Since we retrieve the following expression by Itô's formula, and consequently, Due to: More generally, we have and consequently, Let F = F (t, m) be a C 1,2 -function of time t and "spot" m =m t (u) of the prediction martingale. Because of the filtrations are the same, we have, in general, and applying the Itô formula, we obtain (12). Moreover,
Lemma 4.3 (Auxiliary terms in the decomposition formula for the exponential Volterra volatility model). Let σ t be as in (19) and 0 ≤ t ≤ u. Then Proof. Letφ t (u, z) be given by (16). Then It is now easy to calculate the partial derivative ∂ 2 F . We get Changing variables v = z−mt(u) √r (u|t) and dz = r(u|t) dv, we obtain Remaining formulas follow accordingly.
Lemma 4.5. Let σ t be as in (19) and 0 ≤ t ≤ u. Then, we can re-write F (t,m t (u)) as Moreover, we also have the following equalities Proof. The calculations to obtain these statements are straightforward.
Proposition 4.6 (Terms in the approximation formula for the exponential Volterra volatility model). Let σ t be as in (19) and 0 ≤ t ≤ u. Then and In particular, and Proof. We have that Similarly, we have that For the exponential Volterra volatility model we can determine an upper error bound for the price approximation in the following way.
Theorem 4.7 (Upper error bound for the exponential Volterra volatility model). Let X t be a log-price process (2) with σ t being the exponential Volterra volatility process (19). Then we can express the call option fair value V t using the processes R t , U t from Proposition 4.6. In particular, where t are error terms of order O ξ 2 + ρξ 2 .
Proof. Note that using (24) we have that Applying the Jensen's inequality to (27), we can see that Then, it is easy to find that is a Gaussian process. Therefore 1 a 2 t has finite moments of all orders. Using the error terms specified in the Appendix A and Lemma 3.1, we find the following decompositions for each term Assuming that all previous conditional expectations are finite, the error order is O ξ 2 + ρξ 2 .
Further, we express differentials with respect to the n th -power of the exponential Volterra volatility process when Y t is a semimartingale.
Lemma 4.8. Let σ t be as in (19) and Y t a semimartingale. Let n ≥ 1, we have that Proof. The formula is an immediate consequence of the Itô formula.
Lemma 4.9. Let σ t be as in (19) and Y t is a semimartingale. We can calculate dU t and dR t . In order to simplify the notation, we define the two following functions We define the following auxiliary function Then, it is easier to see that the covariations are the following Proof.
Now, we can re-write U t as and R t as We have that Using Lemma 4.8, we obtain We have that Using Lemma 4.8, we obtain Remark 4.10. Note that the Theorem 4.7 requires that the conditional expectations of (32) and (33) be finite. In the case where σ t is as in (19) and Y t is a semimartingale, it easy to see that these conditions are met using Lemma 4.9.

Exponential fractional volatility model
Let us now focus on a very important example of Gaussian Volterra processes, the fractional Brownian motion (fBm), which is a process with a Hurst parameter H ∈ (0, 1), with covariance function and, in particular, with variance One of the first applications of fractional Brownian motion is credited to Hurst (1951) who modelled the long term storage capacity of reservoirs along the Nile river. However, the idea of this concept goes back to Kolmogorov (1940), who studied Wiener spirals and some other curves in Hilbert spaces. Later, Lévy (1953) used the Riemann-Liouville fractional integral to define the process asB where H may be any positive number. This type of integral turned out to be ill-suited to applications of fractional Brownian motion because of its over-emphasis on the origin for many applications. In their highly cited work, Mandelbrot and Van Ness (1968) introduced the Weyl's representation of the fractional Brownian motion: where and W t is the standard Wiener process. Nowadays, the most widely used representation of fBm is the one by Molchan and Golosov (1969) where .
To understand the connection between Molchan-Golosov and Mandelbrot-Van Ness representations of fBm we refer readers to the paper by Jost (2008). Despite of the above mentioned arguments, Alòs, Mazet, and Nualart (2000) proposed to consider a processB t = t 0 (t − s) H−1/2 dW s instead of B H t in fractional stochastic calculus, since Z t has absolutely continuous trajectories. Since B H t is not a semimartingale, the processB t = Γ(H + 1/2)B H t − Z t is also not a semimartingale. Later on, Thao (2006) introduced the so called approximate fractional Brownian motion process aŝ and showed that for every ε > 0 the processB ε t is a semimartingale and it converges toB t in L 2 (Ω) when ε tends to zero. This convergence is uniform with respect to t ∈ [0, T ] (Thao 2006, Theorem 2.1).
Let us now consider the exponential fractional volatility process where (B H t , t ≥ 0) is one of the above mentioned representations of fBm. We are especially interested in the "rough" models, i.e. when H < 1/2. In this case, we call the model (1) with volatility process (40) the rough fractional stochastic volatility model (αRFSV). Note that if α = 1, we get the rBergomi model (Bayer, Friz, and Gatheral 2016;Gatheral, Jaisson, and Rosenbaum 2018), if α = 0, we get the original exponential fractional volatility model. Values of α between 0 and 1 give us a new degree of freedom that can be viewed as a weight between these two models.
Example 4.11 (Volatility driven by the approximate fractional Brownian motion). Let us consider model (2) with volatility process Then Note that if ε = 0, we get exactly the variance r(t) = t 2H , that it is the variance of the standard fractional Brownian motion. Further we havê and thus Example 4.12 (Volatility driven by the standard Wiener process). If in the previous Example 4.11 we take H = 1/2 and ε = 0, we get model (2) with exponential Wiener volatility process is the standard Wiener process, i.e. whereK(t, s) = 1 {s≤t} . In this case, we have that It is easy to see that dM t = 2ξσ 2 t dW t φ(t, T , α), and thus and For a model without exponential drift (α = 0) these formulas simplify to For matter of convenience, we define the functions and We can re-write U t and R t as and It is easy to find the dU t and dR t , Remark 4.13. We can do a Taylor expansion of U 0 and R 0 to understand better their dependencies, doing that we obtain Theorem 4.14 (Decomposition formula for exponential Wiener volatility model). Let X t be the log-price process (2) with σ t being the exponential Wiener volatility process defined in (41).
Assuming without any loss of generality that the options starts at time 0, then we can express the call option fair value V 0 using the processes U 0 , R 0 from (44) and (46) respectively. In particular, where denotes error terms and for α ≥ 0, | | is at most of the order Cξ( √ T +ρξ 2 )T 3/2 Π(α, T, ξ, ρ). The exact bound is given in Appendix B.
Proof. The detailed proof is given in Appendix B, where we also examine the order of magnitude by the first Taylor term of the integrals.
Remark 4.15. It is worth to mention that the order of the error bound from Theorem 4.14 is better than the general estimate from Theorem 4.7, where the time dependency is not considered. To get finer estimates also for the exponential fractional model (case H = 1/2), a proof similar to the one in Appendix B would have to be performed with more complicated but still tractable calculations.
Example 4.16 (Volatility driven by the standard fractional Brownian motion). Let us consider a model with volatility process where B H t is the standard fractional Brownian motion as defined in (38), i.e. with the Molchan-Golosov kernel (39). Then, the formulas for U 0 and R 0 are given in Proposition 4.6 with particular kernel (39), autocovariance function (36) andr(t, s|u) as in Theorem 4.1. In this case, we do not give the formulas for U 0 and R 0 after substituting the Molchan-Golosov kernel, since these formulas are too long. However, it is worth to mention that the formulas are explicit and numerical evaluation requires only the computation of some multiple Gaussian integrals.

Numerical comparison of approximation formula
In this section, we focus on numerical aspects of the introduced approximation formula. We detail on its numerical implementation and a comparison with the Monte Carlo (MC) simulation framework introduced by Bennedsen, Lunde, and Pakkanen (2017) will be provided.
In the second part of this section, we also introduce two interesting outcome analysis for rBergomi model. In particular, we show how the model can be efficiently calibrated using the approximation formula to short maturity smiles. We remark that classical SV models (e.g. Heston model) might fail to fit the short term smiles, unless they exploit high volatility of volatility levels for which they would be typically inconsistent with the long term skew of the volatility surface.
In what follows, we will inspect the approximation quality for rBergomi model and time to maturity / volatility of volatility ξ scenarios. Based on the nature of error terms (see Appendix A) those two factors should play prominent role when it comes to approximation quality.

On implementation of the approximation formula
We note that for the models studied in this paper, we have obtained either a semi-closed form or analytical formula for standard Wiener case (H=0.5). Moreover, for the class of exponential fractional models -represented by the αRFSV model -we only need to numerically evaluate multiple integrals in R 0 and U 0 .
In our case, this was done using a trapezoid quadrature routine -not necessarily the most efficient approach, but easy to implement. We used a discretisation 2 of integrands such that the numerical error doesn't affect the results in a significant way. I.e. to be lower than standard MC errors when comparing to simulated prices or lower than the expected approximation error.
For benchmarking we use a first-order hybrid MC scheme introduced by Bennedsen, Lunde, and Pakkanen (2017) alongside 50000 MC sample paths. Similarly to the implementation of the approximation formula, we remark that this scheme could be also improved as described in McCrickerd and Pakkanen (2018). 5.2 Sensitivity analysis for rBergomi (α = 1) approximation w.r.t. increasing ξ and time to maturity τ In this section, we illustrate the approximation quality for European call options under various model regimes / data set properties as described in Table 1. We use option expiries up to 1Y -we are expecting a loss of approximation quality, based on the nature of the approximation formula.
Since we utilized a first order approximation arguments with respect to volatility of volatility, we are also expecting more pronounced differences between MC simulations and the formula for large values of ξ.
2 Typically we used from 1000 up to 27000 points for 3D integrals. In Figure 1, we illustrate the approximation quality of the rBergomi approximation for low ξ values. We can observe an expected behaviour: a very good match upto 3M expiry and almost linear deterioration of the quality with increasing τ . Also the approximation formula provides a similar scale of errors across the tested moneyness.
For different moneyness regimes and 1M time to maturity, we obtained the following discrepancies between the MC trials and the introduced formula, measured in the relative option fair value (FV) 3 :
Although the introduced approximation is typically not suitable for calibrations to the whole volatility surface -due to the deterioration of approximation quality when increasing time to maturity -we will illustrate how it can significantly speed-up MC calibration to the provided forward at-the-money (ATMF) backbone.

Short-tenor calibration and a hybrid calibration to ATMF backbone
Unlike previous analyses, which were based on artificial data / model parameter values, we inspect an application fo the formula on the calibration to real option market data. In particular, we utilize four data sets of AAPL options which were analysed in detail by Pospíšil, Sobotka, and Ziegler (2018). Descriptive statistics of the data sets are provided in Table 2. The following calibration test trials will be considered.
• Calibration to short maturity smiles: This should illustrate how well the model can fit short maturity smiles using the introduced approximation formula without exploiting too high volatility of volatility values (ξ). For each data set we selected the shortest maturity slice with more than one traded option. The values were not interpolated by any model, i.e. we calibrated to discrete close mid-prices of traded options. We also confirm, that both MC simulation and the formula reprice the smile with the final calibrated parameters without significant differences.
• Hybrid calibration to the ATMF backbone: In the second trial, we calibrate to each data set ATMF backbone. We note that because we have only a discrete set of traded options we might not have for each maturity an option with strike equal to the corresponding forward. Hence, we take an option with the closest strike to the forward value for each expiry. We use the proposed approximation formula only for τ < 0.2, for longer time to maturities we price by MC simulations.
In both cases, the calibration routine was formulated as a standard least-square optimization problem. I.e. to obtain calibrated parameters, we numerically evaluated where N is the total number of contracts for the calibration, Mid i is the mid price of the i th option and rBergomi i (Θ) represents the corresponding model price based on parameter set Θ. The model price is either obtained by the approximation formula or by means of MC simulations otherwise. The optimization is performed using Matlab's local search trust region optimizer which also needs an initial guess to start with.  All following results will be quoted in relative FV: e.g. rBergomi i (Θ)/S 0 and also differences between market and the calibrated model will be denoted using this measure. For the calibration to the whole surface of European options, errors in FV below 0.5% are typically considered to be acceptable, whereas anything exceeding 1% difference is considered as a significant model inconsistency.
Firstly, we display the results for the short-maturity calibration. In Figure 2, we illustrate that even with a not well suited initial guess for Data # 4, we can obtain satisfactory results (i.e. errors significantly below 0.5% mark). In fact, for all tested data sets, obtained values of the calibrated parameters were not very sensitive to the initial guess (only a number of iterations differed). This is a desired feature which typically is not present under classical SV models, see e.g. Mrázek, Pospíšil, and Sobotka (2016). For two other sets we have obtained qualitatively similar fitting errors, however for the data set # 3 we have retrieved 4 errors (out of 22) with absolute FV difference greater than 0.5% and two even greater than 1%, see Figure 3. We conclude that this was partly caused due to high ξ values compared to other three calibration trials and also slightly longer maturity -the data set #3 includes only one option at the shortest maturity, hence the second shortest was used. Still we can conclude that the obtained errors (also verified by using MC simulations) are overall acceptable, although might not be optimal. We have observed that at least short maturity smiles (< 1M) can be efficiently calibrated using the proposed approximation formula, while for expiries greater than 1M, one would need to stay within a low volatility of volatility regime, otherwise the discrepancies would lead to a nonoptimal solution when recomputed using a more precise (but much more costly) MC simulations. To illustrate efficiency of the approximation will show how the approximation can speed-up ATMFbackbone calibration.
We will now inspect the hybrid calibration where we switch between the approximation and a MC pricer based on properties of options being priced. In particular, we focus on data set from 15 th May (Data # 4) and for τ < 0.2 we will use the approximation formula and MC simulations otherwise. For the calibrated parameters, we will also measure the time spent computing FVs by each pricer. For completeness, we remark that both implementations are not perfect. I.e. for MC simulations under the rBergomi model, the "turbocharging" improvements were introduced by McCrickerd and Pakkanen (2018) recently. On the other hand, numerical integrations in the approximation could be performed by some adaptive quadrature and could be vectorized to improve the computation speed.  In Figure 4, we illustrate calibration fit to the ATMF backbone of the option price surface. We conclude that we have retrieved similar errors for both the prices computed using the proposed approximation and the longer maturity option prices quantified by MC simulations. The final fit of the calibrated model (recomputed by MC simulations) is very good, especially considering that the studied model has only 4 parameters. Moreover, only a fraction of the time spent by MC pricer was needed to computed all FV using the approximation formula. In particular, 98.43% of the pricing time 4 we were computing MC simulation estimates of FVs. We also note that 7/9 of total evaluations were computed by the approximation formula.

Conclusion
In previous sections we studied an approximation approach for the pricing of European options under rough stochastic volatility dynamics. Our approach is based on the option price decomposition results obtained by Alòs (2012) for the standard Heston SV model and on its recent generalisation to other SV models by Merino, Pospíšil, Sobotka, and Vives (2018). The main contribution of our research is to derive pricing formulas suitable for various practical applications by applying the general decomposition on a class of Volterra volatility models. Our main focus is laid on the rough volatility models firstly introduced by Gatheral, Jaisson, and Rosenbaum (2018) and Bayer, Friz, and Gatheral (2016).
In particular, a prediction law for Gaussian Volterra processes was proved and an adapted projection of future volatility was obtained in Section 4 for the class of exponential Volterra volatility models, where the volatility process can be expressed as a positive L 1 function of the time and the Volterra process Y t , i.e. σ t = g(t, Y t ). We focused on one particular example of Gaussian Volterra processes, namely on a (rough) fractional Brownian motion B H t . "Roughness" of sample paths is determined by the Hurst parameter, H ∈ (0, 1/2], where for H = 0.5 we would recover a standard Wiener process.
The pricing formula for European options, which is numerically tractable, is then derived under a newly introduced αRFSV model, i.e. for where α ∈ [0, 1] and r(t) is the corresponding autocovariance function defined in Section 4. This newly introduced model seems to have interesting special cases: • for α = 0 it reverts to a simple exponential RFSV model • and for α = 1 to rBergomi model respectively.
Both special cases of the αRFSV model are studied and praised for their surprising consistency with various financial markets in (Bayer, Friz, and Gatheral 2016) and several other more recent papers.
In Section 5, the newly obtained approximation formula was compared to the MC pricing approach introduced by Bennedsen, Lunde, and Pakkanen (2017) and McCrickerd and Pakkanen (2018). This enabled us to numerically verify the obtained solution, to quantify its approximation errors under various settings and, last but not least, to comment on suitability of the rBergomi model for calibration tasks to real market data based on AAPL stock options 5 . The following conclusions were drawn: • The approximation error is well behaved for short maturities (typically for less than 1M) and the error increases with time to maturity and ξ parameter.
• For medium-term expiries we are able to obtain well approximated prices only under low ξ regimes.
• Although, the approximation under a rough Volterra process involves several numerical integration procedures, it is much faster than MC simulation approach implemented in the same environment 6 . For the standard Wiener case (and in particular for the original Bergomi model), we can have an analytical pricing formula, but also more efficient and well developed MC simulation schemes.
• Considering that the modelling approach studied has only few parameters, we were able to fit the sample market data surprisingly well. For calibration to short maturity smiles, we can use just the approximation formula -this was verified by recomputing calibration errors using MC simulations. For calibration to the whole surface, one can utilize a newly introduced hybrid scheme which consist of a combination of approximation and simulation techniques. The idea is quite simple -to use approximation formula for low maturities or for low ξ values and MC simulations for the remaining computations. Suitability of this scheme was judged by a simple calibration to an ATMF-like backbone. We retrieved a well calibrated model, while saving a significant computational time compared to the calibration based on MC simulations only.
However, we remark that implementation of the approximation could be further improvedfor simplicity we used a simple trapezoidal quadrature to numerically evaluate integrals appearing in U t and R t expressions. Another possible improvement could focus on postulating variance dynamics (in our setting redefining function g) to have a variance process which would be consistent with observed variance swap prices. This would typically require a more complex structure than a simple exponential drift term under the rBergomi model and, possibly, a calibration to vanilla options without less liquid wings of the volatility surface. We conclude that once another volatility function g is postulated, it should be only a matter of algebraic operations to obtain the corresponding approximation formula.

Funding
The work of Jan Pospíšil was partially supported by the Czech Science Foundation (GAČR) grant no. GA18-16680S "Rough models of fractional stochastic volatility". The work of Josep Vives was partially supported by Spanish grant MEC MTM 2016-76420-P.

A Decomposition formulas in the general model
In this appendix, we give the error terms for a decomposition of the general model.
The term (I) can be decomposed as The term (II) can be decomposed as B Upper-bound for decomposition formula for the exponential Wiener volatility model In this appendix we obtain the upper-bound for the decomposition formula for the exponential Wiener volatility model in Example 4.12.
The above upper-bound error is difficult to interpret. In order to do this, we derive a Taylor expansion for one of the terms. Then, the following error behaviour is retrieved: C ρξ 3 σ t T 3/2 6(α − 2) 2 32 √ 2ρ −(α − 2)ξ 2 + 21 √ T .