The complete Gaussian kernel in the multi-factor Heston model: Option pricing and implied volatility applications

In this paper, we propose two new representation formulas for the conditional marginal probability density of the multi-factor Heston model. The two formulas express the marginal density as a convolution with suitable Gaussian kernels whose variances are related to the conditional moments of price returns. Via asymptotic expansion of the non-Gaussian function in the convolutions, we derive explicit formulas for European-style option prices and implied volatility. The European option prices can be expressed as Black–Scholes style terms plus corrections at higher orders in the volatilities of volatilities, given by the Black–Scholes Greeks. The explicit formula for the implied volatility clearly identiﬁes the effect of the higher moments of the price on the implied volatility surface. Further, we derive the relationship between the VIX index and the variances of the two Gaussian kernels. As a byproduct, we provide an explanation of the bias between the VIX and the volatility of total returns, which offers support to recently proposed methods to compute the variance risk premium. Via a series of numerical exercises, we analyse the accuracy of our pricing formula under different parameter settings for the one-and two-factor models applied to index options on the S&P500 and show that our approximation substantially reduces the computational time compared to that of alternative closed-form solution methods. In addition, we propose a simple approach to calibrate the parameters of the multi-factor Heston model based on the VIX index, and we apply the approach to the double and triple Heston models.


Motivation
The well-known Heston (1993) model provides a natural generalization of the Black and Scholes approach to option pricing by introducing stochastic dynamics for the volatility of returns.With its ability to reproduce several empirical features in the dynamics of asset prices, such as the leverage effect and the clustering of volatility, the Heston model has become one of the most widely used stochastic volatility models in the derivatives market.While the Heston model can generate smiles and smirks, it does not provide sufficient flexibility to capture the shape of the implied volatility surfaces, in particular, the largely independent fluctuations in slope and level over time.Another drawback of the He-Here, we propose two new closed analytical solutions based on two integral representation formulas for the conditional marginal density function, which we express as a convolution with suitable Gaussian kernels.These formulas provide a natural way to connect the margins of the multi-factor Heston model to the probability density function of the Black-Scholes model.Specifically, we extract two Gaussian kernels ( Theorems 2.1 and 2.2 ).The first kernel, G 0 , well known in the literature, has a variance 0 that is independent of the volatilities of volatilities (vols of vols) and is given by the integrated conditional mean of the point-in-time variance.The other kernel, G 2 , is hidden in the marginal probability density and has therefore not been explored in the literature.The variance 2 coincides with the conditional variance of continuously compounded returns.We show that G 2 is the "complete" kernel of the Heston model and its multi-factor generalization.Our representation formulas allow us to express the price of any derivative contract, not just vanilla contracts, as the price in a Black-Scholes world, with variance 0 or 2 , convolved with a suitable function that does not depend on the specific payoff of the contract.
The convolution formulas can be computed by solving the integrals numerically.While numerical integration methods are extremely powerful in terms of accuracy, they do not provide an explicit link between the structural properties of the model and the characteristics of the prices.In addition, while closed-form solutions are particularly useful for model calibration, common practice is to calibrate the implied volatility observed in the market, rather than the option prices, because implied volatility is a standardized measure of option prices that makes them comparable even when the underlying assets are not the same.Unfortunately, exact closed-form solutions for implied volatilities are not available in the Heston and multi-factor Heston framework.Therefore, easy-to-implement analytical approximations based on perturbation and asymptotic methods have become popular.Approximations not only help to accelerate the calibration to marketobserved quantities but also enhance the understanding of the analytical features of the model and the implied volatility surface.

Literature review
The earliest and best known asymptotic results are from Lewis (20 0 0) , who derived an asymptotic expansion for small values of the vols of vols.This result was followed by Lee (2001) , who obtain similar results assuming a slow mean reversion of volatility, and Fouque and Lorig (2011) , who assume fast mean-reverting volatility.Additionally, Antonelli and Scarlatti (2009) make an expansion around zero correlation.Friz, Gerhold, Gulisashvili, and Sturm (2011) derive an asymptotic expansion for the implied volatility of the Heston model for a large strike price.Forde and Jacquier (2009) obtain the small-time behavior of the implied volatility in the Heston model (with correlation), while Forde and Jacquier (2011) use large deviation techniques to obtain the small-time behavior of the implied volatility for general stochastic volatility models with zero correlation.Kristensen and Mele (2011) do not approximate the asset price directly but develop a power series expansion of the expected bias that would arise if the Black-Scholes model was used to price derivatives when the true market dynamics obeyed the Heston model.Drimus (2011) follows a similar approach using a different series expansion and shows how the convexity in volatility, measured by the Black-Scholes Volga, and the sensitivity of delta with respect to volatility, measured by the Black-Scholes Vanna, impact option prices in the Heston model.Fouque, Papanicolaou, Sircar, and Solna (2011) derive an asymptotic expansion for general multiscale stochastic volatility models using combined singular and regular perturbation theory.Bergomi and Guyon (2011) also consider multi-factor stochastic volatility models and derive an approximation for the volatility smile at the second order in the vols of vols.Their results coincide with those of Lewis (20 0 0) in the case of the Heston model.Lorig, Pagliarini, and Pascucci (2017) derive a family of asymptotic expansions for European-style option prices and implied volatilities for a general class of local stochastic volatility models.
Some authors have derived asymptotic expansions in a jumpdiffusion stochastic volatility setting (see Berestycki, Busca, & Florent, 2004;Medvedev & Scaillet, 2007 ).More recently, Jacquier and Lorig (2014) provide an explicit implied volatility approximation for any model with an analytically tractable characteristic function, which includes both affine stochastic volatility and exponential Lévy models.Nicolato and Sloth (2012) and Takahashi and Yamada (2012) develop asymptotic expansions around the Black-Scholes model for stochastic volatility models with jump diffusion.Pagliarini and Pascucci (2013) add jumps to a local-stochastic volatility model.Benhamou, Gobet, and Miri (2009) employ Malliavin calculus to develop an approximation formula under the one-factor Heston model with time-dependent parameters.Their option prices are given by a Black-Scholes term plus corrections related to the Greeks of the option.Nagashima, Chung, and Tanaka (2014) extend these results to the general multi-factor Heston model with timedependent parameters and find a similar expansion but with an extra term that captures the interaction between the different variance factors.Alós et al. (2012) use Malliavin calculus to study the short-term behavior of implied volatility for jump-diffusion models with stochastic volatility.Veng et al. (2019) derive an asymptotic expansion for put prices, extending the results of Benhamou et al. (2009) to the general multi-factor Heston model with timedependent parameters.
In line with this literature, we propose an asymptotic expansion of the conditional marginal density for small values of the vols of vols.The main difference with respect to the literature is that the expansion is done after extracting the Gaussian kernels, i.e., we expand only the function that is convolved with the Gaussian kernels.This approach yields particularly interesting results when the G 2 kernel is used, given that its dependence on the vols of vols is fully retained.With our approach, we naturally obtain option prices that can be expressed as the Black-Scholes price plus correction terms related to the Greeks of the options.Similarly, the implied volatility can be written as the square root of the integrated conditional variance plus corrections due to higherorder risks.These decompositions provide a clear understanding of how option prices and implied volatility respond to changes in the model parameters and underlying quantities, which is very important in practice for hedging purposes.

Main contribution
Our paper contributes to the existing literature in several respects.First, we provide two new exact formulas for the conditional marginal density of the multi-factor Heston model ( Theorems 2.1 and 2.2 ).As mentioned above, each formula expresses the marginal probability density as a convolution of a Gaussian kernel whose variance is related to the price-return process.While we do not compute these formulas numerically, this approach avoids some numerical challenges in computing the complex integrals involved in option pricing in the multi-factor Heston model.
Second, following the approach of Zhang, Shu, and M. (2010) ; Zhang, Zhen, Sun, and Zhao (2017) , we derive analytical formulas for the higher-order cumulants in the multi-factor Heston frame-

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ] work.These formulas show that the variance of the Gaussian kernel, G 2 , is given by the variance of the continuously compounded return.Third, we provide explicit approximation formulas for the marginal density function as a Gaussian kernel plus corrections using an asymptotic expansion for the non-Gaussian term in the convolution that defines the marginal probability.With this expansion, we obtain explicit formulas for European vanilla call and put options prices that can be expressed as the Black-Scholes price plus corrections at higher orders of the vols of vols (see Proposition 3.1 ).These formulas satisfy the put-call parity equation at any order of approximation.The formulas using the Gaussian kernel G 0 , while derived with a different approach, are equivalent to those of Bergomi and Guyon (2011) for the Heston model and are similar to those of Veng et al. (2019) for the multi-factor Heston model.However, these authors explicitly compute the first-order corrections in the slow and fast time scales, so they do not capture the effect of price skewness on the convexity of the implied volatility.The formulas using the Gaussian kernel G 2 are new and outperform those obtained using the Gaussian kernel G 0 for out-of-themoney options.
Finally, we derive an approximation for the implied volatility from a second-degree polynomial function of the forward moneyness.This formula allows the effect of price skewness on the asymmetry of the volatility smile to be clearly identified, in addition to the level, slope, and curvature of the implied volatility smirk.For the Heston model with zero drift, our formula is analogous to that of Bergomi and Guyon (2011) , and is in line with the approximation formula in a non-parametric setting proposed by Zhao, Zhang, and Chang (2013) .The key insight of our formula is explicit expressions for the level, slope and convexity in terms of the cumulative uncertainty of the asset price and the integrated volatility process.This is a new result that has not been reported in the literature.
Our work also contributes to the growing literature that explores the bias between the VIX index and the integrated conditional mean of the point-in-time variance, and it provides a possible interpretation of the variance risk premium.Defining the variance risk premium as the positive difference between the second cumulants in the physical and risk-neutral measures (in line with Zhao et al., 2013 ), we compute the premium explicitly in the multi-factor Heston model.In particular, we show that the risk-neutral second cumulant coincides with the square root of 2 /T , which is an implied volatility.Furthermore, we show that the VIX index can be associated not only with 0 (i.e., the integrated conditional variance) but also with 2 (i.e., the variance of the compounded return in the risk-neutral measure) opening the opportunity to calibrate the parameters of the multi-factor Heston model directly from the VIX.This result also provides support to recently proposed methods to compute the variance risk premium from model-free option-implied volatility measures.We also propose an explanation for the bias usually observed between the VIX index and the volatility of total returns.This topic is discussed further in Section 3.3 .
Finally, we provide a one-dimensional integral representation formula for European call and put options in the multi-factor Heston model following the approach in Recchioni and Sun (2016) .These formulas are used as an exact benchmark against which to test the accuracy of our option price approximations.
The rest of this paper is organized as follows.In Section 2 , we review the multi-factor Heston model, derive the main results of the paper, i.e., the two representation formulas for the conditional marginal density function, and introduce the two Gaussian kernels G 0 and G 2 .In Section 3 , we derive approximation formulas, in powers of the vols of vols, for option prices and implied volatility, and we provide an interpretation of the volatility smile.We then derive the relationship between the variances of the two kernels and the VIX index.In Section 4 , we present two simulation studies to assess the accuracy of our approximated formulas, and in Section 5 , we present empirical analyses to assess the effectiveness of our approach in terms of model calibration and forecasting option prices one day ahead.In Section 6 , we show empirically that the squared VIX is better approximated by 2 than 0 , and we use the VIX index to calibrate the parameters of the double and triple Heston models.Our results suggest that the dynamics of the third factor may be influenced by changes in macroeconomic conditions.Section 7 concludes.The proofs of the main results are given in Appendix A , while Appendices B and C report the derivation of the formulas for the option pricing with the expansion based on the Gaussian kernel G 0 ( Appendix B ) and the true marginal density ( Appendix B ). Supplementary material with detailed proofs and some additional results is available online.

Multi-factor Heston model treatment
In this section, we present the multi-factor Heston stochastic volatility model and the main theoretical results of the paper.The final goal is to derive an explicit, approximate expression for the price of European call and put options and for the implied volatility in the multi-factor Heston framework.The key results are two representation formulas for the conditional marginal density function (which is the starting point for the derivation of the option prices) associated with the multi-factor Heston model.The first representation formula shows that the conditional marginal density can be expressed as the convolution of a Gaussian kernel, that does not depend on the Heston vols of vols parameters, and a function that includes all the effects of the vols of vols.The second formula reveals the complete Gaussian kernel, i.e., the one that includes all the effects of the vols of vols and is able to fully capture the process dynamics.
The multi-factor Heston model ( Christoffersen et al., 2009 ) assumes the following stochastic volatility model: where x t denotes the log-price variable, v 1 ,t , . . ., v n,t is the corresponding variances, r(t ) is the instantaneous risk-free rate (assumed to be known in advance), χ j , v * j , and γ j are positive constants, and Z j,t , W j,t , j = 1 , 2 , . . ., n, are standard Wiener processes.All correlations among the Wiener processes are zero, except for E(dZ j,t , dW j,t ) = ρ j dt, where ρ j ∈ (−1 , 1) , j = 1 , 2 , . . ., n are constant correlation coefficients.Dividends are not included.The system of Eqs. ( 1) -( 2) is equipped with the following initial conditions: where S 0 and v j, 0 , j = 1 , 2 , . . ., n, are the initial spot price and variance respectively, which are assumed to be random variables concentrated at a point with probability one.
As specified in Heston (1993) , the quantities χ j are the speeds of mean reversion, v * j represents the long-term means, and γ j denotes the local variances (or volatilities of volatility) of each volatility process v j .These parameters are assumed to be positive, so the process is well defined.
Notably, if the Feller condition is enforced, i.e., 2 χ j v * j /γ 2 j > 1 , the variances v j,t are positive for any t > 0 with probability one
Furthermore, we use γ , v to denote the vectors containing the vols of vols, γ = (γ 1 , γ 2 , . . ., γ n ) , and the variances, v = (v 1 , v 2 , . . ., v n ) , respectively.The transition probability density function (pdf) associated with the stochastic differential system (1), ( 2) is denoted by where R denotes the set of real numbers, R n is the n -dimensional Euclidean vector space, and R n + the positive orthant.We also introduce the processes X t and Y t associated with the multi-factor Heston model (1) : where F t is the information set, i.e., the continuous σ -algebra generated by the point-in-time volatility processes, and E(v j,s |F t ) is the conditional mean of the point-in-time variance given by According to Zhang et al. (2017) with which is related to processes X t and Y t as follows In the following, we provide a representation formula for the conditional marginal density function, which enables extraction of the Gaussian kernel underlying the multi-factor Heston model ( Theorems 2.1 and 2.2 ).Specifically, we use G to denote the Gaussian kernel with variance (t , t ) , t < t , that is: We extract two Gaussian kernels 2 , identified in Theorems 2.1 and 2.2 , that we denote as the zero-order kernel G 0 and the secondorder kernel G 2 .As we show later, the terms "zero-order" and "second-order" reflect the fact that they contain, respectively, no powers of γ and all terms of second degree in γ .
Theorem 2.1 shows that the marginal probability density of the log-price variable can be written (see Eq. ( 18) ) as the convolution of the Gaussian kernel G 0 (independent of the vols of vols) and the function L γ , which accounts for the vols of vols effect.
Theorem 2.1.The marginal probability density of the log-price vari- where ı is an imaginary unit and E(v j,s | F t ) is the conditional mean (7) .Here, B j is given by where ζ j and ν j are the following quantities: Furthermore, M can be written as: where 0 (t , t ) is given by: where G 0 is the Gaussian kernel in (10) , computed for (t , t ) = 0 (t , t ) , and L γ is a function that accounts in full for the effects of the vols of vols: Proof.See Appendix A .
Building from the previous result, we derive Theorem 2.2 , which provides an alternative representation of the marginal density function expressed as the convolution of the Gaussian kernel G 2 and the function L * γ .
Theorem 2.2.The marginal probability density of the log-price vari- where L * γ is: Here, 2 is defined by

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ] and S 1 and S 2 are given by: while H j is given by: Furthermore, the following expansion holds: where S 1 and S 2 are given in Eqs. ( 21) and ( 22) , while S 2 c , S 3 c and S 3 d are: where ψ j is given by: Proof.See Appendix A .
We note that the expansion of H j in powers of γ j , as γ j → 0 + , with degree greater than two involves polynomial functions of k with degree greater or equal to three (see Appendix A ). Thus, the function L * γ contains only terms of order k or k n with n ≥ 3 .All terms of order k 2 are absorbed in G 2 .Therefore, we call G 2 the complete kernel of the multi-factor Heston model, and formula (18) is used to derive the option prices and implied volatility approximations in the next section of this paper 3 .
Te functions L γ in (17) and L * γ in (24) satisfy the following equation: 3 Higher-order Gaussian kernels have been discussed in Wand and Schucany (1990) .
where L γ and L * γ are the Fourier transforms of the functions L γ and L * γ with respect to the log-price, respectively.To provide intuition for the two Gaussian kernels, we derive the relationship between their variances, 0 and 2 , and the processes X t , Y t and R t t , defined by Eq. ( 5) , Eq. ( 6) and Eq. ( 8) .
Proposition 2.3.Let t < t and X t , Y t , and R t t be the processes in ( 5) , ( 6) and ( 8) .We have the following expressions for the conditional moments of X t and Y t : where 0 , S 1 , S 2 , S 2 c , and S 3 c are given in ( 16) , ( 21) , ( 22) , ( 25) and ( 26) 4 .Finally, in the multi-factor Heston model ( 1) , the conditional variance of the continuously compounded return R t t and the price skewness formula, as defined in Das and Sundaram (1999) , are: and where 2 is given in (20) .
Proof.The proof follows using the approach proposed in Zhang et al. (2017) .A detailed proof is given in the supplementary material online.
Notably, the proposed formulas hold for the expectation both in the risk-neutral and physical probability measures.This distinction is necessary in the discussion of the variance risk premium (in Section 3.3 ).
Interestingly, the variance of the Gaussian kernel G 0 is given by the second-order conditional moment of the process X t in (5) and is independent of the vols of vols, while the variance of the Gaussian kernel G 2 coincides with the conditional variance of the continuously compounded return R t t and, through its dependence on the vols of vols, γ j , fully captures the dynamics of the multi-factor Heston model.This makes G 2 the most natural kernel representation of the conditional marginal M. Furthermore, Eq. ( 33) shows that the function S 1 is responsible for the price skewness and the mixed moment between the cumulative uncertainty of the asset return and the uncertainty of the integrated variance process over the time interval [ t , t ] .In the next section, we show that the price skewness, Skewness DS , given in (33) , appears in the coefficient of the second-order term of the implied volatility in Eq. ( 63) and may cause the "volatility smile" convexity to change.
As a corollary of Theorem 2.2 , we provide expansions of the conditional marginal density M in powers of the vols-of-vols vector γ up to the third order.
Proof.The proof is based on the expansion in powers of the vols of vols of the function L γ , and it is available in the Supplementary material online.
We denote the approximations of the marginal density up to the third order as Proposition 2.5 below shows that the approximations of the marginal density in Eq. ( 38) satisfy the conditions that guarantee mass conservation, the martingale property (i.e., the asset price should be a martingale in the multi-factor Heston model) and the so-called symmetry condition.These conditions avoid normdefecting and martingale-defecting pdfs, as discussed in Lewis (20 0 0) Chapter 2.
Proposition 2.5.Let M 0 , M 1 , M 2 and M 3 be given in ( 38) .The following equations then hold which represent mass conservation ( 39) , the martingale property ( 40) and the symmetry condition ( 41) .These properties also hold for the marginal density M in (11) .
Proof.See the supplementary material online.
We conclude this section by emphasizing that Corollary 2.4 shows that the third-order expansion of the marginal density continues to involve only the Gaussian kernel G 2 , confirming that G 2 is the complete Gaussian kernel of the multi-factor Heston model.
In the next section, we use Corollary 2.4 to derive closed-form formulas for the option prices and implied volatility.

Option pricing
In this section, we derive explicit formulas for European vanilla call and put options starting from the representation of the multifactor Heston conditional marginal M provided in Theorem 2.2 and its approximations up to the third order in the vols of vols, given by Eq. ( 38) .The equivalent derivation starting from the representation of the conditional marginal M in terms of the Gaussian kernel G 0 , Eq. ( 15) is provided in Appendix B .
We use C(S 0 , T , E) and P (S 0 , T , E) to denote the price of European vanilla call and put options in the multi-factor Heston model, with spot price S 0 , maturity T , strike price E, and discount factor B (T ) , which is given by Specifically, C and P are defined as: and where v 0 is a vector of the variances at time t = 0 .
Furthermore, we use C BS (S 0 , T , E, T ) and P BS (S 0 , T , E, T ) to denote the classic Black-Scholes formulas for vanilla call and put options, where = (0 , T ) > 0 is the integrated variance over the time interval [0 , T ] , that is, and

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ] where N(x ) is given by and d 1 ( ) and d 2 ( ) are given by Proposition 3.1.Let C(S 0 , T , E) and P (S 0 , T , E) be the prices of European call and put options, respectively, with spot price S 0 , maturity T , strike price E and discount factor B (T ) , as given in Eqs. ( 43) -( 44) .
We have Here, 2 (0 , T ) is given by ( 20) , C BS and P BS denote the classic Black-Scholes formulas in ( 45) and ( 46) , and R 1 , R 2 and R 3 are corrections to the standard Black-Scholes formula due to the contribution of the first-, second-, and third-order correction terms of the expansion in powers of the vols of vols of the function L * γ (see, Eq. ( 19) ): and where S 1 , S 2 , S 2 c , S 3 c and S 3 d are given in ( 21) , ( 22) , ( 25) , ( 26) and ( 27) , respectively.The notation [ •] (•, •, •) in Eq. ( 53) and Eq. ( 54) means that the function in the square brackets is evaluated at the argument (•, •, •) , Note that for γ = 0 , G 2 coincides with G 0 , the correction terms R 1 , R 2 and R 3 become zero, and the option prices become the classic Black and Scholes prices for options with time-dependent but deterministic volatilities.
Dropping the arguments of 2 , S 1 , S 2 and S 2 c , Eqs. ( 52) , ( 53) can be rewritten as where m E is the log-moneyness associated with the forward price defined as Proof.See Appendix A .
We denote the approximated European vanilla option prices up to the the first-, second-and third-order approximations as Notably: (i) Proceeding further with the expansion of the function L * γ in the powers of the vols of vols, we can only add higherorder corrections to the option price approximations without affecting the zero-order contribution.The Black-Scholestype term is, in fact, determined by the Gaussian kernel G 2 , which is not affected by higher-order expansions in γ of L * γ .

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ] (ii) When the vols of vols go to zero, the option prices converge to the Black-Scholes-like term with volatility 0 (0 , T ) /T .
(iii) The Black-Scholes-type formulas for the European vanilla options (that is, the zero-order approximations) overprice at-the-money options.The first-order correction term R 1 , which affects the call and put options in the same way, can correct for this overpricing.In fact, R 1 is negative when the options are at-the-money (i.e., E/ (S 0 e T 0 r(s ) ds ) ≈ 1 ) and the correlations are negative.This finding indicates that in the case of the Heston model, where negative correlation values are usually observed, the prices of call and put options are smaller than those calculated using the standard Black-Scholes formulas for at-the-money options, thereby reducing the overpricing of the Black-Scholes formulas.(iv) The correction term, R 1 , shows why S 1 may be deemed responsible for the smile asymmetry.We observe zero price skewness (33) when ρ j = 0 , j = 1 , 2 , . . ., n .Indeed, zero correlation implies a null third-order correction to option pricing, indicating the crucial effect of non-null correlations.(v) The correction terms R m , m = 1 , 2 , 3 are the same for the call and put options.As a consequence, the pairs C m , P m , m = 1 , 2 , 3 , satisfy the put-call parity.In fact, The fact that the put-call parity holds is implied by the fact that the Fourier transforms of M 1 , M 2 and M 3 with respect to the log-price are equal to zero when the conjugate variable is equal to zero and to the imaginary unit (see Section 8 of the online supplementary material.)(vi) The correction terms R m , m = 1 , 2 , 3 , are linear in the Vega of the Black-Scholes formulas (see Eqs. ( 55) and ( 131) ).
Thus, small values of Vega imply small corrections.Note that for large values of γ j , the Vega goes to zero as e − 2 / 8 .Thus, for large values of γ j , the second-and third-order approximations of the option prices move toward the Black-Scholes-like term with volatility 2 (0 , T ) /T .A higherorder approximation is needed in this case to capture nonzero correction terms.In Section 4.1 , we numerically determine the range of values of γ that are coherent with expansion to the third order.(vi) Theorem 2.2 implies that any contract with maturity T and payoff P that allows for a closed or semi-closed form in the Black-Scholes framework can be written as a convolution of the Black-Scholes price with integrated variance 2 (0 , T ) and the function L * .Using the expansion in powers of the vols of vols of the Fourier transform of L * , which implies a representation of L * as a weighted sum of the derivatives of the Dirac delta function of the log-price, we obtain an expansion of the contract price given by the Black-Scholes price at zero order plus corrections at higher orders given by the Black-Scholes Greeks5 .Moreover, the current representation shows that the corrections to the Black-Scholes term are equal for the put and call options, at any order of approximation, implying that the put-call parity equation is satisfied at any order of approximation.

Implied volatility
The implied volatility in the multi-factor Heston model is defined as the quantity such that the following equality holds: We derive the first-and second-order approximations of as a function of the vols of vols (i.e., = ( γ )) by solving Proposition 3.2.The first-order, 1 ( γ ) , and second-order, 2 ( γ ) , approximations are given by Here, m E is the log-moneyness associated with the forward price (see Eq. ( 57) ), 0 and S 1 are defined in ( 16) and a 0 (T , γ ) , a 1 (T , γ ) and a 2 (T , γ ) are given by a 0 (T , γ ) = with S 2 and S 2 c given in ( 22) , ( 25) , respectively.Here, we have dropped the arguments (0 , T ) of the functions 0 , S 1 , S 2 and S 2 c .
Proof.See Appendix A .
The fact that the implied volatility expansion depends only on 0 and not on 2 is a consequence of the choice of γ = 0 as a base point of the Taylor expansion of the implied volatility.In fact, the same formula for the implied volatility can be derived using the second-order approximation to the call option price based on the Gaussian kernel G 0 .A suitable double expansion would allow a similar formula to be obtained for the implied volatility, with 0 replaced by 2 .This approach, however, is not reported in this paper, as it deserves further investigation.
Notably, the implied volatility resulting from the second-order approximation to the option price is a quadratic function of the forward moneyness and reduces to the approximation of Bergomi and Guyon (2011) in the case of the Heston model.The coefficients a 0 (T , γ ) and a 2 (T , γ ) are second-degree homogeneous functions of γ , while a 1 (T , γ ) is a homogeneous function of degree one.

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ] The expression of a 2 reveals that the convexity of the volatility smile depends on the function S 2 1 .This finding confirms that the quantity S 1 is responsible for the asymmetry in the smile since it allows for a concave smile.Bearing in mind that 0 , S 2 and S 2 c are non-negative for any time horizon and model parameters, Eq. ( 68) clearly shows the effect of price skewness on the volatility smile, i.e., large values of price skewness can destroy the U shape of implied volatility.Indeed, concave volatility smiles are allowed in mean-reverting underlying assets, where the option tenor is comparable to the characteristic reversion time of the asset 6 .
The second-order approximation of the implied volatility leads to the following approximation for the implied volatility skew: For null correlation coefficients S 1 = 0 , S 2 c = 0 , and IV skew = 0 , the second-order approximation, 2 , of the volatility surface is a strictly convex function with vertex at m E = 0 (i.e., when the option is at the money): Finally, a simple calculation proves that the implied volatility skew decays according to 1 / √ T as T → + ∞ since we have lim Therefore, we obtain lim i is the weight of the jth long-term variance mean.The limit for large maturity shows that, in the multi-factor Heston model, the interaction between the variances plays a crucial role in the implied volatility skew (see the squared term on the right-hand side of Eq. ( 71) ), as previously observed in Veng et al. (2019) .

The VIX index
The VIX volatility index, disseminated by the Chicago Board Options Exchange (CBOE), is built to provide a model-free, optionimplied, return volatility measure for the S&P 500 index.The 6 Some empirical evidence can be found at http://faculty.baruch.cuny.edu/jgatheral/Bachelier2008.pdf (see pages 53-56).
CBOE 7 computes the VIX from non-zero bid prices of European call and put options on the S&P 500 index using the formula where T is 30-day maturity, E i is the strike of the i th out-of-themoney option, F t ,t + τ = S t e rτ is the forward index quotation with constant interest rate, S t = e x t is the price at t, and E 0 is the first strike below the forward index level.The quantity Q (E i ) is the midpoint of the bid-ask spread of each option with strike E i .This definition is based on the following representation of the expected value of the future realized variance, given by Demeterfi, Derman, Kamal, and Zou (1999) : where P and C are put and call prices.Eq. ( 73) can then be rewritten (see Zhang et al., 2010 ) via a second-order Taylor approximation of the log function for E 0 ∼ F 0 ,T as: Jiang and Tian (20 05, 20 07) discuss the potential biases that can arise from approximating ( 74) with ( 72) such as (i) truncation errors (the minimum and maximum strikes are far from zero and infinity in practice); (ii) discretization errors (piecewise linear functions approximate the integrals in equation); (iii) expansion errors (the Taylor series expansion is truncated to the second order); and (iv) interpolation errors (linear interpolation of the maturities).In fact, a number of empirical studies indicate that the VIX overestimates the future volatility of the underlying assets.To improve the fit between the VIX index and the volatility of the underlying assets, Pacati, Renó, and Pompa (2018) proposed a new specification in the double Heston model that leads to a deterministic nonnegative shift, or displacement φ t , of the stochastic volatility level such that: Here, we take a different approach and show that the squared VIX can be associated with both 0 and 2 , which in our framework, are both candidates for the implied volatility.In fact, by taking the zero-order approximation in Eq. ( 60) , i.e., only the Black-Scholes term, it is trivial to derive 0 ( γ ) = 2 .Similarly, if we expand the option price formulas around the Gaussian kernel G 0 , we find 0 ( γ ) = 0 at the zero order.To derive the link between the VIX, 0 and 2 , we use different approximations of Eq. ( 74) .The starting point in both cases is to replace the following identities: Here, to maintain simple notation, f (S) denotes the price density.
Using ( 76) and (77) in Eq. ( 74) , we obtain: To derive the relationship between the VIX and 0 , we use the following expansion (see Zhang et al., 2010 ): In fact, neglecting the third-and higher-order terms in Eq. ( 79) and replacing the expansion of the second term of the right-hand side of Eq. ( 78) , we obtain Assuming a constant risk-free interest rate, we have Note that the recent paper by Huang, Schlag, Shaliastovich, and Thimme (2020) , which proposes a stochastic volatility model with stochastic vols of vols, reports a similar result.The squared VIX, in fact, is shown by the authors to be equal to the conditional mean of the integrated variance, which, as in our model, coincides with 0 .
To derive the relationship between the VIX and 2 , we first use the following Taylor expansion on the right-hand side of Eqs. ( 76) -( 77) , then the expansion of S T /E 0 around the conditional mean, denoted by Therefore, we first have where Var is the variance.Taking the variance of Eq. ( 82) , we obtain: Thus, when E 0 is very close to F 0 ,T , as in the CBOE computation, 2 (0 , T ) /T is a proxy for V IX 2 / 100 .This formula is in line with Theorem 1 from Chow, Jiang, and Li (2020) .Note that the variance is computed in the risk-neutral measure.
This finding suggests that the VIX index typically overestimates the conditional mean of the integrated spot volatility (given that normally 2 > 0 , as the correlation ρ j between prices and volatility is typically negative), and our analysis is in line with the results of Pacati et al. (2018) .In fact, our derivation provides an interpretation of the displacement parameter in Pacati et al. (2018) Our analysis, not only provides an explanation for the bias between the VIX and the conditional mean of the integrated spot volatility but also offers a theoretical foundation for recently proposed methods to compute the variance risk premium (VRP).Zhao et al. (2013) provide a natural definition of the variance risk premium, expressing it as the difference between the variances of the continuously compounded returns evaluated in the risk-neutral scenario and in the physical probability measures.In our framework, the variance risk premium for the multi-factor Heston model can be expressed as where Q and P denote, respectively, the risk-neutral and physical probability measures (here Q 2 is the same quantity denoted earlier as 2 ).In fact, as discussed above, Q 2 and P 2 , which are the variances of the Gaussian kernel underlying the price process in the risk-neutral and physical measures, coincide with the variances of the continuously compounded returns under the same two measures.Bon Bondarenko (2014) , Bollerslev, Tauchen, and Zhou (2009) , and Carr and Wu (2009) have proposed constructing the volatility risk premium based on the assumption that model-free option-implied volatility measures can provide a natural empirical analog to the market's risk-neutral expectation of the conditional total variation of returns.These authors had the correct intuition to use the VIX index as a proxy of the risk-neutral variance of returns, and by deriving this link explicitly, via 2 , we provide a theoretical justification for this approach.
Finally, we note that our results imply a nonlinear effect of the vols-of-vols risk on the VIX index, given that 2 is a quadratic function of the vols-of-vols parameters.Huang et al. (2020) also uncovered a nonlinear effect of the vols-of-vols risk on VIX options.

Accuracy of the option price approximations: simulation study
In this section, we study the accuracy of the approximation formulas derived in Sections 2 and 3 in reproducing European option prices and their performance in terms of computational time.As a benchmark, we compute the "true" European option prices by following the approach proposed in Recchioni and Sun (2016) (derivations are reported in Appendix C ).
In the following, we use the subscripts "H", "DH" and "TH" to denote option prices and their approximations in the Heston, double Heston and triple Heston frameworks.

Simulation study 1: Heston and double Heston on "reasonable" grid of parameters
We being this section by assessing the performance of the second-and third-order approximations, Eqs. ( 58) -( 59) , of the call and put option prices, C m,H , P m,H , with m = 2 , 3 , in the Heston framework.
The Heston exact formula is obtained by imposing n = 1 in Eq. ( 134) in Appendix C .Eqs. ( 132) and ( 133) in Appendix C are equal except for the values of q, which are valid over different intervals.In the following, we choose q = 1 .05 for a call option and q = −0 .05 for a put option8 .Eqs. ( 132) and ( 133) are defined via convergent integrals that can be computed accurately using a simple composite rectangular rule.
We evaluate the exact formulas C H and P H and the approximated formulas C m,H and P m,H for the points in the following set: These values of model parameters in grid M include those estimated by Christoffersen et al. (2009) in Section 4.2 (see, also the online supplementary material).Some descriptive statistics for the call and put option prices, computed with the exact formulas C H and P H , are shown in Table 1 .
Table 2 compares the exact option prices with their second-and third-order approximations.From left to right we report the vols of vols ( γ ), mean (mean C ), median (median C ), and standard deviation (std C ) of the relative call option errors, e C,m , and the mean (mean P ), median (median P ), and standard deviation (std P ) of the relative put option errors, e P,m , associated with the second-order approximation, m = 2 (in the top panel), and the third-order approximations, m = 3 (in the bottom panel).
The results in Table 2 show that while the quality of the approximations decreases as γ increases, the second-order approximation guarantees four correct significant digits up to a volatility of 50%.The average error of the second-order approximation is at most of the order of a percent for larger values of γ .The thirdorder approximation improves the estimation by less than one order of magnitude for the values of γ considered.Given that only marginal improvements are obtained with the third-order approximation, we focus on the second-order approximations when presenting numerical and empirical results in the following sections.
In the remainder of this section, we illustrate the computational advantages of using the second-order approximation formulas ( 58) and (59) in the Heston and double Heston models.To this end, we consider the same grid M as in Eq. ( 86) and compute 3125 call and put option prices for each value of γ , averaging over the other parameters on the grid, using the secondorder approximations C 2 ,H and P 2 ,H .In the case of the double Heston model, we have chosen This choice is made to limit the number of call and put options to be evaluated to 31250 for each value of γ , as in the Heston model.We then evaluate the number of points to be [m5G;December 23, 2020;19:58 ] Table 3 From left to right: average, min and max number of points, avg N p , min N p , and max N p , required by the rectangular rule to achieve the same accuracy as that of the second-order approximation formulas in the Heston (upper panel) and double  132) and (133) using 2 16 quadrature points.Table 3 shows that using formulas ( 58) and ( 59) allows considerable savings in computation time with respect to using the integral formulas for both the Heston model (top panel) and double Heston model (bottom panel).This computation time reduction is important because, for the same level of accuracy, the time required to evaluate option prices with the integral formulas in the double Heston model is, in the best case, approximately twice that needed for the Heston model.

Simulation study 2: Heston and double Heston with empirical parameters
In this subsection, we repeat the previous exercise using model parameters calibrated to real data.Specifically, we use the parameters estimated by Christoffersen et al. (2009) 9 for the Heston and double Heston models for the years 1990 to 2004.The spot variance of the Heston model is chosen to be 0.9, while the spot variances of the double Heston model are v 1 = 0 .13 and v 2 = 0 .75 .These choices are supported by the results of the empirical analysis discussed in Christoffersen et al. (2009) p. 1926. In fact, in Christoffersen et al. (2009) , the sum of the factor estimates v 1 , 0 and v 2 , 0 is 88% in the two-factor model, and the difference is approximately 62%, while it is 90% in the one-factor model.Finally, the risk-free interest rate is chosen to be 0.15.
As the first step, we compute the prices of 25 European vanilla call and put options with spot price S 0 = 100 , strike prices E = 80 + j/ 5 , j = 1 , 2 , . . ., 5 and time to maturity T = j/ 12 years, j = 1 , 2 , . . ., 5 using the exact integral formulas with 2 16 nodes.As in Section 4.1 , these values are denoted the "true values".Then, we compute the average relative errors, over the twenty-five options, for the put and call options second-order approximations in the Heston and double Heston frameworks.Finally, we determine the number of nodes necessary to achieve, with the integral formulas, the same average accuracy of the second-order approximation formulas and compare the computational times of the two methods.Tables 4 and 5 show the results of this experiment, respectively, for the Heston and double Heston frameworks.The columns in these tables are the same as those in Table 3 , with the only difference being that the time and accuracy are computed for a specific set of model parameters, average across the strike and time to maturity.
We observe that, on average, the relative error of the secondorder approximations is 0.02% for both put and call options in the Heston framework and 6.1% and 5.4%, respectively, in the double Heston model.These relative errors guarantee four correct significant digits for the Heston model and two correct significant digits for the double Heston model.The discrepancy in the accuracy between the two tables is due to the different magnitude of the vols positive probability unless, as remarked in Christoffersen et al. (2009) , the process satisfies a standard reflecting barrier at the origin.Interestingly, the Feller condition never holds in the case of process v 1 ,t .Violation of the Feller condition has also been noted in Pacati et al. (2018) .

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ] Table 4 From left to right: year, vol of vol γ 1 , number of points N p required by the rectangular rule to achieve the same accuracy as that of the second-order approximation formulas in the Heston model; Time H , Time 2 ,H time required to compute the fifty put and call options with the integral formulas and the BS-second-order approximations; Avg.rel.err P, • and Avg.rel.err.C, •: average relative errors on put and call options of the integral formulas with N p points and the second-order approximations.For each set of model parameters estimated by Christoffersen et al. (2009) (see Table 3 Panel A) over the years 1990-2004, we compute European put and call options with spot price S 0 = 100 , time to maturity T = j/ 12 years, and strike prices E j = 80 + 10( j − 1) , j = 1 , 2 , 3 , 4 , 5 .The risk-free interest rate is r = 0 .15 .The computation was conducted on an Intel CORE i7 (8th generation) processor.The true values are obtained with the integral formulas ( 132) and ( 133) using 2 16 points.From left to right: year, vols of vols γ 1 , γ 2 , number of points N p required by the rectangular rule to achieve the same accuracy as that of the second-order approximation formulas in the double Heston model; Time DH and Time 2 ,DH time required to compute the fifty put and call options with the integral formulas and the second-order approximations; Avg.rel.err P, • and Avg.rel.err.C, •: average relative errors on put and call options of the integral formulas with N p points and the second-order approximations.
For each set of model parameters estimated by Christoffersen et al. (2009) (see Table 3 , Panel B) over the years 1990-2004, we compute European put and call options with spot price S 0 = 100 , time to maturity T = j/ 12 years, and strike prices E j = 80 + 10( j − 1) , j = 1 , 2 , 3 , 4 , 5 .The risk-free interest rate is r = 0 .15 .The computation was conducted on an Intel CORE i7 (8th generation) processor.The true values are obtained with the integral formulas ( 132) and ( 133) using 2 16 points.In these years, we observe the largest relative errors for the double Heston model.In contrast to that of the double Heston model, the estimated vol of vol of the Heston model is always less than 80% and larger than 37% (see Tables 4 and 5 ).

Double Heston model
We conclude this section by comparing the relative errors of the call option prices given by the zero-order Black-Scholes-type term (with Gaussian kernel G 2 ) and by the first-, second-and third-order approximations as a function of vol of vol in the case of the Heston model and as a functions of the largest of the two vols of vols (i.e., γ 1 ) in the case of the double Heston model.In the left panel of Fig. 1 , the model parameters are the same as those used for Tables 4 and 5 .In this case, the prices for different levels of vols of vols are not perfectly comparable, as they also depend on the remaining model parameters.To isolate the effect of the vols of vols and expand the range of values considered for this parameter, in the central panel, we plot the approximated option prices when fixing the model parameters to the values estimated for the year 1990 (see Table 8 Year 1990), while the vol of vol γ 1 is chosen to be γ m 1 = e −3+ m/ 2 , m = 1 , 2 , . . ., 30 , and γ 2 = 0 .007 .The grid of strike prices and times to maturity, over which the average is taken, is the same for both panels.In the right panel, we show for the same parameters as the central panel, the correction terms R i , for i = 1 , 2 , 3 .The curves log-error vs log vol-of-vol (left and central panels) depicted in Fig. 1 show that the errors grow linearly with the vols of vols for values of γ up to approximately 200%.The figures also provide empirical evidence that while the third-order expansion slightly improves the approximation for volatilities up to approximately 200%, beyond this value, the second-and third-order approximations become indistinguishable from each other and con- verge to the Black and Scholes price, in line with the discussion in Section 3.1 .Interestingly, for larger values of the vols of vols, and in line with the theory of asymptotic series (as further explained in Section 7 of the supplementary material), the first-order approximation provides better estimates than higher-order approximations.This is a consequence of the fact that for large values of the vols of vols, the asymptotic expansion of L * γ may diverge.This is also signalled by the correction terms R 2 and R 3 becoming larger than R 1 .However, despite the non-convergence, the asymptotic expansion may still provide a satisfactory approximation when truncated to a finite number of terms.

Accuracy of the option price approximations: empirical calibration study
In this section, we assess the performance of the Heston second-order approximation formula ( 58) -( 59) to reproduce and to forecast traded European call and put option prices on the US S&P 500 index.In this exercise, the U.S. three-month government bond index is used as a proxy for the interest rate r.
The availability of an explicit and elementary formula for the implied volatility provides an advantage in terms of calibrating the model rather than estimating the parameters directly from the option prices.This is because it avoids biases caused by different magnitudes of option prices that are typically corrected by introducing appropriate weights in the optimization algorithm (i.e., the inverse of option Vegas, see Christoffersen et al., 2009 , or the bidask spread, see Date & Islyaev, 2015 ).Additionally, the simple link between implied volatility and model parameters allows for reliable estimates while accelerating the solution of the optimization problem.We note that, while formulas similar to ours for the implied volatility (i.e., Eq. ( 63) ) were derived by Bergomi and Guyon (2011) , their effectiveness for calibration purposes has not been tested in the literature.
Here, we provide empirical evidence that by using the secondorder approximations for the implied volatility 2 ,H , we can obtain "consistent" estimates of the Heston model parameters from both the call and put options.Typically, option prices are filtered to avoid inconsistency resulting from the simultaneous use of call and put option prices (see Pacati et al., 2018 ).We do not filter any observations that do not satisfy standard no-arbitrage conditions while investigating how this affects the model calibration.
Our dataset consists of 1200 European vanilla call and put options with four strike prices (i.e., n E = 4 ) and n T = 150 maturities.
Starting from the traded call option prices C o (S i , T i , E j ) with spot price S i , time to maturity T i and strike price E j , and using the U.S. three-month government bond yield as the risk-free interest rate, r, we compute the observed implied volatility, σ o C (S i , T i , E j ) , for i = 1 , 2 , . . ., n T , j = 1 , 2 , . . ., N E .This computation is performed using the Matlab function calcBSImpVol , which uses Li's rational function approximator for the initial estimate (see, Li, 2006;Li, 2008 ), followed by Householder's root finder of the third order to improve the convergence rate of the Newton-Raphson method.
For any time i = 1 , 2 , . . ., n T , we then estimate the Heston model parameters solve the optimization problem: min where 2 ,H is given in formula ( 63) with n = 1 and V is the following set of constraints: To solve problem (87) , we use a metric variable steepest descent algorithm (see, for example, Recchioni & Scoccia, 20 0 0 ; Fatone, Mariani, Recchioni, & Zirilli, 2013 ).This is an iterative algorithm that generates a sequence of points, k , k = 0 , 1 , . . ., belonging to the interior of the feasible region and moving opposite to the gradient vectors of the objective function computed in a suitable metric.
We then repeat the calibration procedure starting from the observed put prices P o (S i , T i , E j ) , where P o is the observed value of the put option, i = 1 , 2 , . . ., n T , and j = 1 , 2 , . . ., n E , and solve the problem min In this way, we obtain two optimal sets of model parameters, one starting from the call options, C , and the other starting from the put options, P .Some descriptive statistics for the estimated model parameters, initial variance, Feller ratio, objective function and observed implied volatility are given for the two sets in Table 6 .The values of the objective function compare favourably with those in Table 1 of Veng et al. (2019) .The two sets of parameters are almost identical, with the exception of the estimate of the long-term mean parameter.We argue that the difference in the v * parameter estimate from the call and put prices is due to market imperfections that lead to a spread between the implied volatility σ o of call and put options.In fact, the absolute value of the implied volatility spread, | derived from the call and put options is 0.04 on average, while the relative absolute spread (i.e., the ratio of the spread to implied volatility from the call) is 0.24.Interestingly, the absolute difference between the square root of the two long-term variance parameters is 0.05, and the ratio of this difference to the square root of the call variances is 0.29, thus mirroring the implied volatility spread.
To evaluate the model consistency, we compute the European call and put option prices using formulas C 2 ,H and P 2 ,H in Eqs. ( 58) -( 59) with both sets of estimated parameters.Fig. 2 shows

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ]   the observed and second-order (solid line and dotted line, respectively) call option prices.The approximations in Fig. 2 are obtained using the model parameters estimated by the observed implied volatility from call options (i.e., call set).The corresponding figures for the put prices, obtained using the model parameters estimated by the observed implied volatility from put options (i.e., put set) are available in the supplementary material.For each set, we compute the mean and standard deviation of the relative errors for the call options as: and we also compute the equivalent errors for the put options.The average relative errors E C, C and E P, P (i.e., when parameters are estimated starting from the corresponding option prices) are, respectively, 0.027 (i.e., 2.7%) and 0.031 (i.e., 3.1%).These errors are in line with those in Pacati et al. (2018) , where a double Heston model with jumps is used.By contrast, when using the model parameters of the put set to estimate the call prices, and vice versa, the relative errors E C, P and E P, C are, on average, 0.21 (i.e., 21%) for the call and 0.22 (i.e., 22%) for the put options.Thus, while the cross estimates produce a clear bias, the error is of the same order as the relative error in the implied volatility (i.e., 24%), suggesting that the bias is driven by market imperfections rather than an inconsistency with the methodology.We conclude this section by testing the potential of the calibrated parameters to forecast option prices one day ahead.Fig. 3 shows the one-day-ahead estimates for call (left panel) and put (right panel) option prices.Specifically, the option estimates at time t + 1 are calculated using the optimal parameter values at time t.The one-day-ahead estimated call prices are obtained using the model parameters C t , while the one-day-ahead estimated put prices are obtained from P t .The relative errors of the one-dayahead estimates are, on average, 4.67% for call options and 4.72% for put options.

Variance of the Gaussian kernels and the VIX index
In this section, we focus on the relationship between the VIX and the variances 0 and 2 in the Heston, double Heston, and triple Heston models.We use the VIX time series for the years 20 0 0, 20 01, 20 02 and 2003 provided by the CBOE10 and two time series for the realized variance, the median truncated realized variance and the 5-minute realized variance, both available from the Oxford-Man Institute11 .For the Heston model, we assume that the realized variance from the Oxford-Man Institute database plays the role of the spot variance; thus, v t = RV t .The use of the realized variance as a proxy for the short-term volatility factor is supported by the results illustrated in Corsi, Fusari, and La Vecchia (2013) .In the double Heston model, each factor variance is evaluated as a fraction of the total realized variance; thus, v 1 ,t = α 1 RV t and v 2 ,t = (1 − α 1 ) RV t .In the triple Heston model, the stochastic variances v j,t , j = 1 , 2 , 3 are chosen to be a fraction of the total ).Call price one-day-ahead estimates using the call set (left panel); Put price one-day-ahead estimates using the put set (on the right).The average relative errors of call and put options are 7.9% and 6.2%, respectively.

Gaussian kernels and the VIX index in the Heston and double Heston frameworks
We now empirically test the relationship between the VIX, 0 and 2 in the Heston and double Heston frameworks.For each fixed year, we use the model parameters in Christoffersen et al. (2009) , provided in Table 7 and Table 8 , to compute the kernel variance 0 and 2 , while the α 1 parameter is obtained by minimizing the squared residuals (sum of squared errors): across the four years (thus, α 1 is imposed to be the same for the four years).We find that the optimal value of α 1 is α 1 = 0 .15 when we use the median truncated realized variance and α 1 = 0 .06 in the case of 5-minute realized variance.We start by comparing, in Figs. 4 and 5 , the VIX time series (solid line) and ˜ 2 ,model (dotted line) in the Heston ( Fig. 4 ) and double Heston ( Fig. 5 ) models as a function of the day index for each year considered.The figures show that ˜ 2 ,DH (t) (see Eq. ( 90) ) more closely follows the VIX behavior for all years and both time series.This result is confirmed in Table 9 , which shows that the double Heston model outperforms the Heston model in terms of the sup-norm.The RMSE shown in Table 9 compares favourably with the results obtained by Corsi et al. (2013) (see Section 4.3, Table 4 ).
We then compare, in Fig. 6 , the fit between the VIX and ˜ 2 ,model versus ˜ 0 ,model .The figures clearly show, as discussed in Section 3.3 , that the VIX overestimates 0 ,model .
To provide further evidence of this point, we test for linear dependence between the VIX index and ˜ 0 ,model and ˜ 2 ,model with model = H and DH .This is done by regressing the daily VIX observations on the daily estimates of ˜ 0 ,model (t) and ˜ 2 ,model (t) when using, respectively, the median truncated realized variance (see, Table 10 , left panel) and the 5-minute realized variance (see Table 10 , right panel) as proxies of the spot variance process v t = n j=1 v j,t .The results of these zero-intercept regressions show that both ˜ 0 ,model (t) and ˜ 2 ,model (t) perform better than the naive linear model V IX t = β 1 RV t + noise .The coefficients are statistically significant at the 5% level.These results are in line with the findings of Huang et al. (2020) .These results also confirm our hypothesis that ˜ 2 ,model (t) in the Heston and double Heston models captures the VIX dynamics better than ˜ 0 ,model (t) .In fact, in both the Heston and double Heston models, the coefficient β 1 is, on average, closer to one when we regress on ˜ 2 ,model (t) rather than ˜ 0 ,model (t) .
We further investigate the quality of the VIX approximation by analysing the bias, i.e., E( ˜ m,model − V IX ) , for m = 0 , 2 , for the Heston and double Heston models.Table 11 shows that ˜ 0 ,model has a more pronounced bias than ˜ 2 ,model and that, particularly for the years characterized by large vols of vols, the use of ˜ 2 ,model substantially improves the fit.
To provide further intuition about the above results, we compare the accuracy of the double Heston call option pricing formulas when the expansion in the vols of vols is performed starting from the representation in Eq. ( 15) , after extracting the Gaussian kernel 0 , and from the representation in Eq. ( 18) , after extracting the Gaussian kernel 2 .We focus on out-of-the-money call options, which are the ones used to compute the VIX index.A better performance of the approximation formulas, written in terms of 2 , in pricing out-of-the-money options, would provide justification for the better performance of 2 itself in approximating the VIX.
As an illustration, we use the double Heston parameters estimated by Christoffersen et al. (2009) for European call options on the S&P500 in the year 2003 (see the last row of Table 8 ).The first factor of the double Heston model in 2003 is characterized by a high vols of vols (881%), a very small long-term mean (i.e., 0.33%) and a slow mean reverting speed (i.e., 0.1638), so the pa-  7 were used with the spot variance of the price log-return corresponding to the daily time series of the median truncated realized variance (left panels) and the 5-minute realized variance (right panels) from the Oxford-Man Institute.8 were used with the spot variance of the price log-return corresponding to the daily time series of the median truncated realized variance (left panels, α = 0 .15 ) and the 5-minute realized variance (right panels, α = 0 .06 ) from the Oxford-Man Institute.

Table 8
Estimated parameters for the double Heston model from Christoffersen et al. (2009) .

Double Heston model parameters
Year  (i.e., 0.4625).This suggests that the volatility remains closer to its long-term mean, with parameters satisfying the Feller condition.
The initial value of each factor variance is evaluated as a fraction of the total variance, so v 1 , 0 = α and v 2 , 0 = (1 − α) .The risk-free interest rate is chosen to be 0.15.We analyse the accuracy of the approximations as a function of α (which is the only parameter we estimated in the VIX exercise) to assess the sensitivity of the option price approximations to the choice of spot volatility.The comparison is presented in Fig. 7 , where the relative errors of the second-order call option approx-imations in the double Heston framework are shown in a logarithmic scale.The solid lines are the approximations obtained with the Gaussian kernel G 0 , while the dashed lines are those obtained with G 2 .In all panels, a logarithmic scale with logarithm-base 10 is used.The x -axis shows the values of the spot variance v 1 ,t = α.
Fig. 7 shows that the approximations obtained with the complete kernel G 2 are more accurate than those obtained with G 0 for the considered parameter values.Interestingly, some values of the spot variance v 1 , 0 (i.e., v 2 , 0 ) reduce the relative pricing errors.We note that as we move from v 1 , 0 = 0 to v 1 , 0 = 1 , we start closer to or further from the long-term means of the two factors.This movement has the effect of changing the relative contribution of each factor to the overall dynamics, which, as the process transitions [m5G;December 23, 2020;19:58 ]    from smooth mean reverting dynamics with a small vol of vol to dynamics with abrupt fluctuations, could make the option more difficult to price, explaining the observed changes in pricing errors.

ARTICLE IN PRESS
While we restrict the analysis to out-of-the-money options for a specific set of parameters, the comparison of the accuracy of the two approximations when expressed in terms of 0 or 2 deserves a full investigation, which will be the subject of future work.

Calibration of the double and triple Heston models from the VIX index
Given the encouraging results in the previous subsection, we explore the possibility of calibrating the parameters of the double and triple Heston models directly from the VIX daily data.
For this exercise, we use only the median truncated realized variance as a proxy of the spot variance v t .In this case, both the model parameters and α i are obtained by minimizing the squared residuals (SSE) separately for each year (thus, α i differs by year).
We use the Matlab lsqnonlin function to minimize the SSE.As a starting point, we the Heston and double Heston parameters from the previous section.We first calibrate the double Heston model.The results reported in Table 12 are similar to those of Christoffersen et al. (2009) for a large basket of options.The estimated parameters for the triple Heston model are reported in Table 13 , and the relative values of the coefficients α j , j = 1 , 2 , 3 are shown in Table 14 .14 , which is two orders of magnitude lower than the values in Table 11 .The coefficient β 1 is also closer to one.
Notably, the estimation performed using the VIX provides model parameter values for the first two factors that are similar to those obtained by Christoffersen et al. (2009) .In particular, the second factor, which is dominant, is consistently slowly mean reverting around its long-run mean, while the first factor has more volatile dynamics around its very small long-term mean across the whole period.By contrast, the temporal dynamics of the parameters characterizing the third factor appear to switch between these two types of behavior.Specifically, the vol of vol, γ 3 , and the speed of mean reversion, χ 3 , are lower in the years 20 01 and 20 02, while the long-term mean, v *

ARTICLE IN PRESS
JID: EOR [m5G;December 23, 2020;19:58 ] Table 10 Zero-intercept regression models with two proxies for spot variance.The model parameters of the Heston and double Heston models are taken from Table 3   Table 13 Triple Heston model parameters estimated from the VIX and the median truncated realized variance data by minimizing the SSE in Eq. ( 91) .

Conclusions
This paper introduces an approach to extract the Gaussian kernel behind the multi-factor Heston model, which allows a clear connection of the prices of European option contracts in the multifactor Heston framework to the corresponding prices in the Black-Scholes model.Our simple formulas illustrate how the option prices and implied volatility respond to changes in model parameters.A series of numerical exercises shows that our formulas are accurate, computationally efficient, and easy to calibrate.We numerically demonstrate that our approximations compare favourably with other pricing formulas available in the literature, such as those of Pacati et al. (2018) and Veng et al. (2019) , when we use them to calibrate the model parameters to the implied volatility and forecast the option prices one day ahead.

A.3. Proof of Corollary 2.4
We sketch the proof for the first three-order terms, G 2 , M 1 and M 2 .From Theorem 2.2 and using the expansion in formula (24) up to the second order, we have: e S 1 (t ,t )(ı k 3 + ı k )+ S 2 (t ,t )(k 4 −2 ık 3 −ık )+ S 2 c (t ,t )(k 4 −ı k 3 ) dk, (115 where S 1 is a linearly homogeneous function of the vols of vols, while S 2 and S 2 c are homogeneous functions of degree two.We compute the first three terms of the expansion in powers of the vols of vols of the function E ( γ ) = e S 1 (t ,t )(ı k 3 + ı k )+ S 2 (t ,t )(k 4 −2 ık 3 −ık )+ S 2 c (t ,t )(k 4 −ı k 3 ) . (116) The proof follows based on: and This concludes the proof.

A.4. Proof of Proposition 3.1
The proof follows by substituting M with its third-order approximation and integrating by parts.Details on why we obtain an explicit formula for the corrections terms are given in the online supplementary material.As mentioned above, the correction R m,MH , m = 1 , 2 , 3 for the call option is the same as the put correction since there are two changes of sign: one due to the payoff function and the other due to integration by parts over the interval ( −∞ , log E) rather than ( log E, + ∞ ).This concludes the proof.

ARTICLE IN PRESS
(121) To prove Eq. ( 63) , we proceed by computing the second-order derivatives, which are given by: An easy, but involved, computation illustrated in the online supplementary material shows that the addenda containing powers of (m E + 1 2 0 ) higher than two are cancelled by the addenda involving the Black-Scholes Vomma.In fact, we have: (123) Proceeding in a similar manner, we obtain mixed-order mixed derivatives: The thesis follows since we have S 2 = 1 .This concludes the proof.A more detailed proof is provided in the supplementary material section available online.

Appendix B. Formulas in terms of the Gaussian kernel G 0
In this section, we provide the second-order approximations of the option prices starting from the representation of the marginal density function given by formula (11) .
Corollary B.1.The following expansion of the conditional marginal M in powers of γ as γ → 0 holds: M(x, v , t, x , t ) = G 0 (x − x, t, t ) + M 1 , 0 (x, v , t, x , t ) + M 2 , 0 (x, v , t, x , t ) + o γ 2 , γ → 0 , (125) where M 1 , 0 and M 2 , 0 are given by M 1 , 0 (x, v , t, x , t ) (127) Here, S 1 is given by ( 21) , G 0 is the Gaussian kernel defined in ( 10) , and S 2 and S 2 c are given in ( 22) and ( 25 (129) Here, 0 (0 , T ) is given by ( 16) , C BS and P BS denote the classic Black-Scholes formulas, as in ( 45) and ( 46) , and R 1 , 0 and R 2 , 0 are the corrections to the standard Black-Scholes formula due to the contribution of the first-and second-order terms of the expansion in powers of the vols of vols of the marginal density function: (131) Proof.The proof is based on the proof of Corollary 2.4 and Proposition 3.2 for the expansion in powers the vols of vols considering that the Fourier transform of L γ is equal to the product of the Fourier transform of L * γ multiplied by e (k 2 −ık )(S 1 (t ,t ) −S 2 (t ,t )) .

Fig. 1 .
Fig. 1.Left and middle panels: relative errors of zero-, second-, and third-order approximations to the call options obtained in a log-log scale.Zero order: dash-dot line; first order: dashed line; second order: solid line; third order: dotted line.(Right) panel: correction terms R 1 (dashed line), R 2 (solid line), and R 3 (dotted line) as a function of the log of the vols of vols.

Fig. 3 .
Fig. 3. Observed option prices (solid line) and one-day-ahead estimates computed using second-order approximation (dotted line) for four different strike prices E 1 = 1900 , E 2 = 1975 , E 3 = 20 0 0 , and E 4 = 2025 and with expiry date T = December 19, 2015, versus time (September 1, 2014-March 30, 2015).Call price one-day-ahead estimates using the call set (left panel); Put price one-day-ahead estimates using the put set (on the right).The average relative errors of call and put options are 7.9% and 6.2%, respectively.

Fig. 4 .
Fig. 4.Each panel contains the VIX time series and the model implied volatility ˜ 2 ,H (i.e., Eq. (90) -Heston model) as a function of day.The model parameters in Table7were used with the spot variance of the price log-return corresponding to the daily time series of the median truncated realized variance (left panels) and the 5-minute realized variance (right panels) from the Oxford-Man Institute.

Fig. 5 .
Fig. 5.Each panel contains the VIX time series and the model implied volatility ˜ 2 ,DH (i.e., Eq. (90) -Double Heston model) as a function of day.The model parameters in Table8were used with the spot variance of the price log-return corresponding to the daily time series of the median truncated realized variance (left panels, α = 0 .15 ) and the 5-minute realized variance (right panels, α = 0 .06 ) from the Oxford-Man Institute.

Fig. 6 .
Fig. 6.Comparison of V IX, 0 /T and 2 /T (double Heston model) when the spot variance is computed with the RV median truncated realized variance.

Fig. 7 .
Fig. 7. Relative errors of the second-order approximations of call options in the double Heston framework obtained using the Gaussian kernels G 0 (solid-line) and G 2 (dashed line) as a function of the initial spot variance v 1 , 0 = α, v 2 , 0 = (1 − α) .The double Heston parameters used to compute the call option prices are those estimated in the year 2003 and shown in Table8.A logarithmic scale is used for the y -axis (logarithm base = 10).

Fig. 8
Fig. 8 presents the VIX time series (solid line) and ˜ 2 ,T H (dotted line) in the triple Heston model as a function of the day index for each year considered.The fit with ˜ 2 ,T H (t) appears to better capture the VIX behavior than the Heston and double Heston models presented in the previous subsection.The visual inspection is supported by the lower values of the average, minimum and maximum RMSE, which are 0.0187, 0.0126 and 0.0258, respectively, and of the bias, reported in the last column ofTable 14 , which is two orders of magnitude lower than the values in Table 11 .The coeffi-

Fig. 8 .
Fig. 8.Each panel contains the VIX time series and the model implied volatility ˜ 2 ,TH (i.e., Eq. (90) -Triple Heston model) as a function of day.The model parameters are in Table 13 .U.S. GDP growth, which registered 4.1% in 20 0 0, 1.1% in 2001, 1.7% in 2002 and 2.9% in 2003, i.e., the long-term volatility is higher in 2001 and 2002 when the GDP experiences a decline.The third factor may thus have macro-economic significance.This result is not surprising since other studies have shown a correlation between GDP growth and stock market returns (see, Ritter, 2005 and the references therein; Cournéde & Denk, 2015 ) 12 .

Table 1
Descriptive statistics for the exact call and put option prices evaluated on grid M of the Heston model.

Table 2
Descriptive statistics for the relative errors of second-and third-order option price approximations evaluated on grid M in the case of the Heston model.Second-order approximations in vols of vols ( C 2 ,H , P 2 ,H ) Heston (lower panel) models; the total time Time H (Time DH ) and Time 2 ,H (Time 2 ,DH ) required to compute the European options with integral formulas (132) and (133) with N p points and the second-order approximations C 2 ,H , P 2 ,H ( C 2 ,DH , P 2 ,DH ); Avg.rel.err.are the average relative errors of the put and call option approximations with integral formulas C H , P H ( P DH , C DH ) and second-order approximations P 2 ,H , C 2 ,H ( P 2 ,DH , C 2 ,DH).For a fixed vol of vol (or pair of vols of vols), the relative errors are computed by averaging over the remaining parameters in the set M for a total of 31250 option prices.The time is expressed in seconds.The computation was conducted on an Intel CORE i7 (8th generation) processor.The true values are obtained from integral formulas (132) and (133) using 2 16 quadrature points.
used in the quadrature rule to achieve the same level of accuracy when pricing the options with integral formulas Eq. (134) .In the top panel of Table3, we report, from left to right (and for γ = 0 .01,0 .05,0 .15,0 .25 , 0 .5 , 0 .8, 2 ), the average, mean and max number of points (truncated to the closest integer) required by the rectangular rule to achieve the same level of accuracy (i.e., relative error) in option prices as the second-order approximation formulas; the total times, Time H and Time 2 ,H , required to compute 31,250 European call and put options, respectively, with the integral formula and with formulas C 2 ,H and P 2 ,H ; and the average relative errors (i.e, Err.P H , Err.C H , Err.P 2 ,H and Err.C 2 ,H ) of the put and call options with the integral formulas P H and C H and with the second-order approximations P 2 ,H and C 2 ,H .The columns in the bottom panel of Table3are the same as those in the top panel, but the results correspond to the double Heston model.The computations were conducted on an Intel CORE i7 (8th generation) processor.The true values are obtained with the integral formulas (

Table 6
Descriptive statistics for estimated values of the model parameters and observed implied volatility σ o .

Table 9
Root mean square error (RMSE) obtained using ˜ 2 ,model to approximate the VIX index.

Table 11
Comparison of the bias in the estimates obtained with ˜ 0 ,H , ˜ 2 ,H , ˜ 0 ,DH and ˜ 2 ,DH using the two proxies for the spot variance.

Table 14
R 2 and Bias.