Vincent Lemaire

Thibaut Montes

Gilles Pagès

Stationary Heston model: Calibration and Pricing of exotics using Product Recursive Quantization

A major drawback of the Standard Heston model is that its implied volatility surface does not produce a steep enough smile when looking at short maturities. For that reason, we introduce the Stationary Heston model where we replace the deterministic initial condition of the volatility by its invariant measure and show, based on calibrated parameters, that this model produce a steeper smile for short maturities than the Standard Heston model. We also present numerical solution based on Product Recursive Quantization for the evaluation of exotic options (Bermudan and Barrier options).

Introduction

Originally introduced by Heston in [START_REF] Steven L Heston | A closed-form solution for options with stochastic volatility with applications to bond and currency options[END_REF], the Heston model is a stochastic volatility model used in Quantitative Finance to model the joint dynamics of a stock and its volatility, denoted pS pxq t q tě0 and pv x t q tě0 , respectively, where v x 0 " x is the initial condition of the volatility. Historically, the initial condition of the volatility x is considered as deterministic and is calibrated in the market like the other parameters of the model. This model received an important attention among practitioners for two reasons: first, it is a stochastic volatility model, hence it introduces smile in the implied volatility surface as observed in the market, which is not the case of models with constant volatility and second, in its original form, we have access to a semi closed-form formula for the characteristic function which allows us to price European options (Call & Put) almost instantaneously using the Fast Fourier approach (Carr & Madan in [START_REF] Carr | Option valuation using the fast fourier transform[END_REF]). Yet, a complaint often heard about the Heston model is that it fails to fit the implied volatility surface for short maturities because the model cannot produce a steep-enough smile for those maturities (see [START_REF] Gatheral | The volatility surface: a practitioner's guide[END_REF]).

Noticing that the volatility process is ergodic with a unique invariant distribution ν " Γpα, βq where the parameters α and β depend on the volatility diffusion parameters, it has been first proposed by Pagès & Panloup in [START_REF] Pagès | Approximation of the distribution of a stationary markov process with application to option pricing[END_REF] to directly consider that the process evolves under its stationary regime in place of starting it at time 0 from a deterministic value. We denote by pS pνq t q tě0 and pv ν t q tě0 the couple asset-volatility in the Stationary Heston model. Replacing the initial condition of the volatility by the stationary measure does not modify the long-term behavior of the implied volatility surface but does inject more randomness into the volatility for short maturities. This tends to produce a steeper smile for short maturities, which is the kind of behavior we are looking for. Later, the short-time and long-time behavior of the implied volatility generated by such model has been studied by Jacquier & Shi in [START_REF] Jacquier | The randomized heston model[END_REF]. Another extension of the Heston model have been suggested and extensivel analyzed in order to reproduce the slope of the skew for short-term expiring options: the Rough Heston model where the volatility satisfies a Voltera equation driven by a "rough" Liouville process with H-Hölder paths, H " 0.1 (see [START_REF] Jaisson | Rough fractional diffusions as scaling limits of nearly unstable heavy tailed hawkes processes[END_REF][START_REF] Hamza Guennoun | Asymptotic behavior of the fractional heston model[END_REF][START_REF] Gatheral | Volatility is rough[END_REF][START_REF] Callegaro | Fast hybrid schemes for fractional riccati equations (rough is not so tough[END_REF][START_REF] Gatheral | Rational approximation of the rough heston solution[END_REF] for details on the model and numerical solutions).

An other extension of the Heston model have been suggested in order to be able to reproduce the slope of the skew for short-term expiring options: the Rough Heston model (see [START_REF] Jaisson | Rough fractional diffusions as scaling limits of nearly unstable heavy tailed hawkes processes[END_REF][START_REF] Hamza Guennoun | Asymptotic behavior of the fractional heston model[END_REF][START_REF] Gatheral | Volatility is rough[END_REF][START_REF] Callegaro | Fast hybrid schemes for fractional riccati equations (rough is not so tough[END_REF][START_REF] Gatheral | Rational approximation of the rough heston solution[END_REF] for details on the model and numerical solutions).

In the beginning of the paper, we briefly recall the well-known methodology used for the pricing of European option in the Standard Heston model. Based on that, we express the price I 0 of a European option on the asset S pνq T as

I 0 " E " e ´rT ϕpS pνq T q ‰ " E " f pv ν 0 q ‰ (0.1)
where f pvq is the price of the European option in the Standard Heston model for a given set of parameters. The last expectation can be computed efficiently using quadrature formulas either based on optimal quantization of the Gamma distribution or on Laguerre polynomials.

Once we are able to price European options, we can think of calibrating our model to market data. Indeed the parameters of the model are calibrated using the implied volatility surface observed in the market. However, the calibration of the Standard Heston model is highly depending on the initial guess we choose in the minimization problem. This is due to an overparametrization of the model (see [START_REF] Gauthier | Fitting the smile, smart parameters for sabr and heston[END_REF]). Hence, when we consider the Heston model in its stationary regime, there is one parameter less to calibrate as the initial value of the volatility is no longer deterministic. The stationary model tends to be more robust when it comes to calibration.

In the second part of paper, we deal with the pricing of Exotic options such as Bermudan and Barrier options. We propose a method based on hybrid product recursive quantization. The "hybrid" term comes from the fact that we use two different types of schemes for the discretization of the volatility and the asset (Milstein and Euler-Maruyama). Recursive quantization was first introduced by Pagès & Sagna in [START_REF] Pagès | Recursive marginal quantization of the euler scheme of a diffusion process[END_REF]. It is a Markovian quantization (see [START_REF] Pagès | Optimal Quantization Methods and Applications to Numerical Problems in Finance[END_REF]) drastically improved by the introduction of fast deterministic optimization procedure of the quantization grids and the transition weights. This optimization allows them to drastically reduce the time complexity by an order of magnitude and build such trees in a few seconds. Originally devised for Euler-Maruyama scheme of one dimensional Brownian diffusion, it has been extended to one-dimensional higher-order schemes by [START_REF] Thomas A Mcwalter | Recursive marginal quantization of higher-order schemes[END_REF] and to medium dimensions using product quantization (see [FSP18, RMKP17, CFG18, CFG17, PS18]). Then, once the quantization tree is built, we proceed by a backward induction using the Backward Dynamic Programming Principle for the price of Bermudan options and using the methodology detailed in [START_REF] Sagna | Pricing of barrier options by marginal functional quantization[END_REF][START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF] based on the conditional law of the Brownian Bridge for the price of Barrier options.

The paper is organized as follows. First, in Section 1, we recall the definition of the Standard Heston model and the interesting features of the volatility diffusion which bring us to define the Stationary Heston model. In Section 2, we give a fast solution for the pricing of European options in the Stationary Heston model when there exists methods for the Standard model. Finally, once we are able to price European options, we can define the optimization problem of calibration on implied volatility surface. We perform the calibration of both models and compare their induced smile for short maturities options. Once this model has been calibrated, in Section 3, we propose a numerical method based on hybrid product recursive quantization for the pricing of exotic financial products: Bermudan and Barrier options. For this method, we give an estimate of the L 2 -error introduced by the approximation.

The Heston Model

The Standard Heston model is a two-dimensional diffusion process pS pxq t , v is the dynamic of the risky asset,

• v x t is the dynamic of the volatility process,

• S pxq 0 " s 0 ě 0 is the initial value of the process, • r P R denotes the interest rate,

• q P R is the dividend rate,

• ρ P r´1, 1s is the correlation between the asset and the volatility,

• pW, Ă W q is a two-dimensional standard Brownian motion,

• θ ě 0 the long run average price variance,

• κ ě 0 the rate at which v x t reverts to θ,

• ξ ě 0 is the volatility of the volatility,

• v x 0 " x ě 0 is the deterministic initial condition of the volatility. This model is widely used by practitioner for various reasons. One is that it leads to semiclosed forms for vanilla options based on a fast Fourier transform. The other is that it represents well the observed mid and long-term market behavior of the implied volatility surface observed on the market. However, it fails producing or even fitting to the smile observed for short-term maturities.

Remark 1.1 (The volatility). One can notice that the volatility process is autonomous thence we are facing a one dimensional problem. Moreover, the volatility process is following a Cox-Ingersoll-Ross (CIR) diffusion also known as the square root diffusion. Existence and uniqueness of a strong solution to this stochastic differential equation has been first shown in [START_REF] Ikeda | Stochastic differential equations and diffusion processes[END_REF], if x ě 0. Moreover, it has been shown, see [START_REF] Lamberton | Introduction to stochastic calculus applied to finance[END_REF], that if the Feller condition holds, namely ξ 2 ď 2κθ, for every x ą 0, then the unique solution pv x t q tě0 satisfies @t ě 0,

Ppτ x 0 " `8q " 1 (1.2)
where τ x 0 is the first hitting time defined by

τ x 0 " inftt ě 0 | v x t " 0u where inf H " `8. (1.3)
Moreover, the CIR diffusion admits, as a Markov process, a unique stationary regime, characterized by its invariant distribution ν " Γpα, βq (1.4) where α " θβ and β " 2κ{ξ 2 .

(1.5) Based on the above remarks, the idea is to precisely consider the volatility process under its stationary regime, i.e., replacing the deterministic initial condition from the Standard Heston model by a ν-distributed random variable independent of pW, Ă W q. We will refer to this model as the Stationary Heston model. Our first aim is to inject more randomness for short maturities (t small) into the volatility but also to reduce the number of free parameters to stabilize and robustify the calibration of the Heston model which is commonly known to be overparametrized (see e.g. [START_REF] Gauthier | Fitting the smile, smart parameters for sabr and heston[END_REF]).

This model was first introduced by [START_REF] Pagès | Approximation of the distribution of a stationary markov process with application to option pricing[END_REF] (see also [START_REF] Ikeda | Stochastic differential equations and diffusion processes[END_REF], p. 221). More recently, [START_REF] Jacquier | The randomized heston model[END_REF] studied its small-time and large-time behaviors of the implied volatility. The dynamic of the asset price pS pνq t q tě0 and its stochastic volatility pv ν t q tě0 in the Stationary Heston model are given by

$ ' ' & ' ' % dS pνq t S pνq t " pr ´qqdt `av ν t `ρd Ă W t `a1 ´ρ2 dW t dv ν t " κpθ ´vν t qdt `ξa v ν t d Ă W t (1.6)
where v ν 0 " Lpνq " Γpα, βq with β " 2κ{ξ 2 , α " θβ. S pνq 0 , r and q are the same parameters as those defined in (1.1) and the parameters ρ, θ, κ, θ and ξ can be described as in the Standard Heston model.

Pricing of European Options and Calibration

In this section, we first calibrate both Stationary and Standard Heston models and then compare their short-term behaviors of their resulting implied volatility surfaces. For that purpose we relied on a dataset of options price on the Euro Stoxx 50 observed the 26th of September 2019 (see Figure 1). This is why, as a preliminary step we briefly recall the well-known methodology for the evaluation of European Call and Put in the Standard Heston model. Based on that, we outline how to price these options in the Stationary Heston model. Then, we describe the methodology employed for the calibration of both models: the Stationary Heston model (1.6) and the Standard Heston model (1.1) and then we discuss the obtained parameters and compare their short-term behaviors.

European options

The price of the European option with payoff ϕ on the asset S pνq T , under the Stationary Heston model, exercisable at time T is given by

I 0 " E " e ´rT ϕpS pνq T q ‰ . (2.1)
After preconditioning by v ν 0 , we have

I 0 " E " E " e ´rT ϕpS pνq T q | σpv ν 0 q ‰ ı " E " f pv ν 0 q ‰ (2.2)
where f pvq is the price of the European option in the Standard Heston model with deterministic initial conditions for the set of parameters λpvq " ps 0 , r, q, θ, κ, ξ, ρ, vq.

Example 2.1 (Call). If ϕ is the payoff of a Call option then f is simply the price given by Fourier transform in the Standard Heston model of the European Call Option. The price at time 0, for a spot price s 0 , of an European Call Cpλpvq, K, T q with expiry T and strike K under the Standard Heston model with parameters λpvq " ps 0 , r, q, θ, κ, ξ, ρ, vq is

Cpλpvq, K, T q " E " e ´rT pS pvq T ´Kq `‰ " e ´rT ´E " S pvq T 1 S pvq T ěK ‰ ´K E " 1 S pvq T ěK ‰ "
s 0 e ´qT P 1 `λpvq, K, T ˘´K e ´rT P 2 `λpvq, K, T ˘(2.3)

with P 1 `λpvq, K, T ˘and P 2 `λpvq, K, T ˘given by

P 1 `λpvq, K, T ˘" 1 2 `1 π ż `8 0 Re ˆe´iu logpKq iu ψ `λpvq, u ´i, T s0
e pr´qqT ˙du P 2 `λpvq, K, T ˘" 1 2 `1 π ż `8 0 Re ˆe´iu logpKq iu ψ `λpvq, u, T ˘˙du (2.4)
where i is the imaginary unit s.t. i 2 " ´1, ψ `λpvq, u, T ˘is the characteristic function of the logarithm of the stock price process at time T . Several representations of the characteristic function exist, we choose to use the one proposed by [SST04, Gat11, AMST07], which is numerically more stable. It reads

ψ `λpvq, u, T ˘" E " e iu logpS pvq T q | S pvq 0 , x ‰ " e iuplogps 0 q`pr´qqT q ˆeθκξ ´2`p κ´ρξui´dqT ´2 logpp1´g e ´dt q{p1´gqq
ȇv 2 ξ ´2pκ´ρξui´dqp1´e ´dt q{p1´g e ´dt q

(2.5) with d " a pρξui ´κq 2 ´ξ2 p´ui ´u2 q and g " pκ ´ρξui ´dq{pκ ´ρξui `dq.

(2.6) Hence, in (2.2), f pvq can be replaced by C `λpvq, K, T ˘, which yields

I 0 " E " e ´rT pS pνq T ´Kq `‰ " E " C `λpv ν 0 q, K, T ˘ı. (2.7)
Now, we come to the pricing of European options in the Stationary Heston model, using the expression of the density of v ν 0 " Γpα, βq, (2.2) reads

I 0 " E " f pv ν 0 q ‰ " ż `8 0 f pvq β α Γpαq v α´1 e ´βv dv.
(2.8)

Now, several approaches exists in order to approximate this integral on the positive real line.

• Quantization based quadrature formulas. One could use a quantization-based cubature formula with an optimal quantizer of v ν 0 with the methodology detailed in Appendix D. Given that optimal quantizer of size N , p v N 0 , we approximate I 0 by p

I N 0 p I N 0 " E " f pp v N 0 q ‰ " N ÿ i"1 f pv N 0,i q P `p v N 0 " v N 0,i ˘.
(2.9)

Remarks 2.2. In one dimension, the minimization problem, that consists in building an optimal quantizer, is invariant by linear transformation. Hence applying a linear transformation to an optimal quantizer preserves its optimality. For example, if we consider an optimal quantization p X N of a standard normal distribution N p0, 1q then µ `σ p X N is an optimal quantizer of a normal distribution N pµ, σ 2 q and the associated probabilities of each Voronoï centroid stay the same.

In our case, noticing that if we consider a Gamma random variable X " Γpα, 1q then the rescaling of X by 1{β yields X{β " Γpα, βq. Hence, for building the optimal quantizer p v N 0 of v ν 0 , we can build an optimal quantizer of X " Γpα, 1q and then rescale it by 1{β, yielding p v N 0 " p X N {β. Our numerical tests showed that it is numerically more stable to use this approach.

In order to build the optimal quantizer, we use Lloyd's method detailed in Appendix D to X " Γpα, 1q with the cumulative distribution function F X pxq " PpX ď xq and the partial first moment K X pxq " ErX 1 Xďx s given by @x ą 0,

F X pxq " 1 Γpαq γpα, xq, K X pxq " αF X pxq ´xα e ´x Γpαq ,
otherwise, F X pxq " 0, K X pxq " 0,

(2.10)

where γpα, xq " ş x 0 t α´1 e ´t dt is the lower gamma function. And the associated probabilities of the optimal quantizer p v N 0 are given by (D.10)

P `p v N 0 " v N 0,i ˘" P `p X N " x N i ˘" F X `xN i`1{2 ˘´F X `xN i´1{2 ˘(2.11)
where @i P 2, N , x N i´1{2 "

x N i´1 `xN i 2
and x N 1{2 " 0 and x N N `1{2 " `8.

• Quadrature formula from Laguerre polynomials. One could also use an algorithm based on fixed point quadratures for the numerical integration. Indeed, noticing that the density we are integrating against is a gamma density which is exactly the Laguerre weighting function (up to a rescaling). Then, I 0 rewrites

I 0 " ż `8 0 f pvq β α Γpαq v α´1 e ´βv dv " β α Γpαq ż `8 0 f pvqωpvqdv (2.12)
where ωpvq " v α´1 e ´βv is the Laguerre weighting function. Then, for a fixed integer n ě 11 , I 0 is approximated by

r I n 0 " β α Γpαq n ÿ i"1 ω i f pv i q (2.13)
where the ω i 's are the Laguerre weights and the v i 's are the associated Laguerre nodes.

Calibration

Now that we are able to compute the price of European options, we define the problem of minimization we wish to optimize in order to calibrate our models parameters. Let P SH be the set of parameters of the Stationary Heston model that needs to be calibrated, defined by (S 0 " 3541, r " ´0.0032 and q " 0.00225) The expiries T are given in days and the strikes K in percentage of the spot.

P SH " pθ, κ,
The others parameters are directly inferred from the market: we get S 0 " 3541, r " ´0.0032 and q " 0.00225. In our case, we calibrate to option prices all having the same maturity. The problem can be formulated as follows: we search for the set of parameters φ ‹ P P that minimizes the relative error between the implied volatility observed on the market and the implied volatility produced by the model for the given set of parameters, such that P " P SH for the Stationary Heston model and P " P H for the Standard Heston model. There is no need to calibrate the parameters s 0 , r and q since they are directly observable in the market.

Being interested in the short-term behaviors of the models, it is natural to calibrate both models based on options prices at a small expiry. Once the optimization procedures have been performed, we compare their performances for small expiries. For that, we calibrate using only the data on the volatility surface in Figure 1 with expiry 50 days (T " 50{365) and then we compare both models to the market implied volatility at expiry 22 days which is the smallest available in the data set.

Remark 2.3. The calibration is performed in C++ on a laptop with a 2,4 GHz 8-Core Intel Core i9 CPU using the randomized version of the simplex algorithm of [START_REF] John | A simplex method for function minimization[END_REF] proposed in the C++ library GSL. This algorithm is a derivative-free optimization method. It uses only the value of the function at each evaluation point. The computation time for calibrating the Standard Heston model is around 20s and a bit more than a minute for the Stationary model. However, these computation times need to be considered carefully because the calibration time highly depends on the initial condition we choose for the minimizer and on the implementation of the Call pricer in the Standard Heston model.

Optimization without penalization

We want to find the set of parameter φ ‹ that minimizes the relative error between the volatilities observed in the market and the ones generated by the model, hence leading to the following minimization problem

min φPP ÿ K ˆσMarket iv pK, T q ´σModel iv pφ, K, T q σ M arket iv pK, T q ˙2 (2.16)
where T is the expiry of the chosen options chosen a priori and K are their strikes. σ M arket iv pK, T q is the Mark-to-Market implied volatility taken from the observed implied volatility surface and the implied volatility σ M odel iv pφ, K, T q is the Black-Scholes volatility σ that matches the European Call price in this model to the price given by the Standard or Stationary Heston model with the set of parameters φ.

In all the following figures, the strike K is given in percentage of the spot S 0 . It is clear in Figure 2 (right) that both models fit really well to the market data and more precisely, the Stationary model succeeds to calibrate with the same precision as the Standard one with one less parameter. Moreover, one notices that even for 22 days maturity options, the Standard Heston model tends to over-estimate the implied volatility and fails to produce the right smile whereas the Stationary Heston model is closer to the market observations. Now, we extrapolate the implied volatility surfaces, given by the two models, for even smaller maturities (7 and 14 days) in order to analyze the behavior of each model for short-term expiries. Now, we investigate how these models behave for longer maturities. Do they succeed in preserving the general shape of the market volatility surface or are they only correctly fitting the maturity on which we calibrated them?

Figure 5 represents the relative error between the implied volatility given by the market and the one given by the models calibrated models at 50 days. Clearly, one notices that the Standard Heston model only fits at this expiry. Indeed, when looking at the expiry 22 days or for long-term maturities, the relative error explodes. The term-structure of the implied volatility surface of the market is not preserved when using the Standard Heston model. However, the Stationary Heston model does fit well at both short and long term expiries. The Stationary model produces a steep smile for very short maturities and flattens correctly to the appropriate mean for long expiries. Table 1: Parameters obtained for both models after calibration without penalization for options with maturity 50 days (S 0 " 3541, r " ´0.0032 and q " 0.00225).

However, looking closely at the parameters obtained after calibration (which are summarized in Table 1), one notices that both sets of calibrated parameters are far from satisfying the Feller condition. And we have to keep in mind that the calibration procedure is performed in order to price path-dependent or American style derivatives using Monte-Carlo simulation or alternative numerical methods, as developed in the next Section. Hence, the Feller condition has to be satisfied, this is the reason why we add a constraint to the minimization problem in order to penalize the sets of parameters not satisfying the condition.

Optimization with penalization using the Feller condition

The minimization problem becomes

min φPP ÿ K ˆσMarket iv pK, T q ´σModel iv pφ, K, T q σ M arket iv pK, T q ˙2 `λ maxpξ 2 ´2κθ, 0q (2.17)
where λ is the penalization factor to be adjusted during the procedure. The obtained parameters after calibration are summarized in Table 2. The Feller condition is still not fulfilled for both models but it is not far from being satisfied. We choose λ " 0.01 which seems to be right the compromise in order to avoid underfitting the model because of the constraint.

φ ‹ ρ v 0 θ κ ξ
Standard Heston ´0.83 0.0045 0.17023 2.19 1.04

Stationary Heston ´0.99 0.02691 19.28 1.15

Table 2: Parameters obtained for both models after calibration with penalization (λ " 0.01) for options with maturity 50 days (S 0 " 3541, r " ´0.0032 and q " 0.00225).

Figure 6 displays the resulting implied volatility curves at 50 days and 22 days for both calibrated models and observed in the market with calibration at 50 days. Adding a penalization term deteriorates the calibration results compared to the non-penalized case (see Figure 2 (right)) but the results are still acceptable. Now, again, we extrapolate the implied volatility of both models for very short term maturities in Figure 7. The Stationary Heston model produces the desired smile, however the Standard Heston model fails to produce prices sensibly different than 0 for strikes higher than 105 with this set of parameters, this is why there is no values in implied volatility curves. Figure 8 represents, as in the non-penalized case, the relative error between the implied volatility given by the market and the one given by the models calibrated models at 50 days using a penalization. The Standard Heston model completely fails to preserve the term-structure while being calibrated at 50 days. In comparison, the Stationary Heston behaves much better and the relative error does not explodes for long-term expiries, meaning that the long run average price variance is well caught.

Toward the pricing of Exotic Options

In this Section, we evaluate first Bermudan options and then Barrier options under the Stationary Heston model. For both products, the pricing rely on a Backward Dynamic Programming Principle. The numerical solution we propose is based on a two-dimensional product recursive quantization scheme. We extend the methodology previously developed by [START_REF] Fiorin | Product markovian quantization of a diffusion process with applications to finance[END_REF][START_REF] Callegaro | American quantized calibration in stochastic volatility[END_REF][START_REF] Callegaro | Pricing via recursive quantization in stochastic volatility models[END_REF], where they considered an Euler-Maruyama scheme for both components. In this paper, we consider a hybrid scheme made up with an Euler-Maruyama scheme for the log-stock price dynamics and a Milstein scheme for the (boosted) volatility process. Finally, we apply the backward algorithm that corresponds to the financial product we are dealing with (the Quantized Backward Dynamic Programming Principle for Bermudan Options, see [START_REF] Bally | A quantization algorithm for solving multidimensional discrete-time optimal stopping problems[END_REF][START_REF] Bally | A quantization tree method for pricing and hedging multi-dimensional american options[END_REF][START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF] and the algorithm by [START_REF] Sagna | Pricing of barrier options by marginal functional quantization[END_REF][START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF] for Barrier Options based on the conditional law of the Brownian motion).

Discretization scheme of a stochastic volatility model

We first present the time discretization schemes we use for the asset-volatility couple pS pνq t , v ν t q tPr0,T s . For the volatility, we choose a Milstein on a boosted version of the process in order to preserve the positivity of the volatility and we select an Euler-Maruyama scheme for the log of the asset.

The boosted volatility. Based on the discussion in Appendix A, we will work with the following boosted volatility process: Y t " e κt v ν t , t P r0, T s for some κ ą 0, whose diffusion is given by dY t " e κt κθdt `ξ e κt{2 a Y t d Ă W t .

(3.1)

The Milstein discretization scheme of Y t is given by Other schemes could have been used, see [START_REF] Alfonsi | On the discretization schemes for the CIR (and bessel squared) processes[END_REF] for an extensive review of the existing schemes for the discretization of the CIR model, but in our case we needed one allowing us to use the fast recursive quantization, i.e., where we can express explicitly and easily the cumulative distribution function and the first partial moment of the scheme, which is the case of the Milstein scheme (we give more details in SubSection 3.2).

s Y t k`1 " M r b,r σ `tk , s Y t k , Z 2 k`1 ˘(3.2) with t k " T k n and
Hence, as our time-discretized scheme is well defined because its positivity is ensured if the Feller condition is satisfied, we can start to think of the time-discretization of our process pS pνq t k q kP 0,n .

The log-asset. For the asset, the standard approach is to consider the process which is the logarithm of the asset X t " logpS t q. Applying Itô's formula, the dynamics of X t is given by dX t " ´r ´q ´vt 2 ¯dt `?v t dW t .

(3.6) Now, using a standard Euler-Maruyama scheme for the discretization of X t , we have

# s X t k`1 " E b,σ `tk , s X t k , s Y t k , Z 1 k`1 s Y t k`1 " M r b,r σ `tk , s Y t k , Z 2 k`1 ˘(3.7)
where

Z 1 k`1 " N p0, 1q, Z 2 k`1 " N p0, 1q, CorrpZ 1 k`1 , Z 2 k`1 q " ρ and E b,σ pt, x, y, zq " x `bpt, x, yqh `σpt, x, yq ? h z (3.8)
with bpt, x, yq " r ´q ´e´κt y 2 and σpt, x, yq " e ´κt{2 ? y.

(3.9)

Hybrid Product Recursive Quantization

In this part, we describe the methodology used for the construction of the product recursive quantization tree of the couple log asset-boosted volatility in the Heston model. In Figure 9, as an example, we synthesise the main idea behind the recursive quantization of a diffusion v t which has been time-discretized with F 0 pt, x, zq. We start at time t 0 " 0 with a quantizer p v 0 taking values in the grid Γ t 0 " tv 0 1 , . . . , v 0 10 u of size 10, where each point is represented by a black bullet (') with probability p 0 i " Ppp v 0 " v 0 i q is represented by a bar. In the Stationary Heston model, p v 0 is an optimal quantization of the Gamma distribution given by (1.4) and (1.5). Then, starting from this grid, we simulate the process from time t 0 to time t 1 " 5 days with our chosen time-discretization scheme F 0 pt, x, zq, yielding r v 1 " F 0 pt 0 , p v 0 , Z 1 q, where Z 1 is a standardized Gaussian random variable. Each trajectory starts from point v 0 i with probability p 0 i . And finally we project the obtained distribution at time t 1 onto a grid Γ t 1 " tv 1 1 , . . . , v 1 10 u of cardinality 10, represented by black triangles (IJ) such that p v 1 is an optimal quantizer of the discretized and simulated process starting from quantizer p v 0 at time t 0 " 0.

Remark 3.1. In practice, for low dimensions, we do not simulate trajectories. We use the information on the law of r v 1 conditionally of starting from p v 0 . The knowledge of the distribution allows us to use deterministic algorithms during the construction of the optimal quantizer of r v 1 that are a lot faster than algorithms based on simulation.

v1 v 1 i p 1 i = ℙ( v1 = v 1 i ) v 0 i p 0 i = ℙ( v0 = v 0 i )
Figure 9: Example of recursive quantization of the volatility process in the Heston model for one time-step.

In our case, we consider the following stochastic volatility system

# dX t " bpt, X t , Y t qdt `σpt, X t , Y t qdW t dY t " r bpt, Y t qdt `r σpt, Y t qd Ă W t (3.10)
where W t and Ă W t are two correlated Brownian motions with correlation ρ P r´1, 1s, b and σ are defined in (3.9) and r b and r σ are defined in (3.3). Our aim is to build a quantization tree of the couple pX t , Y t q at given dates t k , k " 0, . . . , n based on a recursive product quantization scheme. The product recursive quantization of such diffusion system has already been studied by [START_REF] Callegaro | Pricing via recursive quantization in stochastic volatility models[END_REF] and [START_REF] Rudd | Fast quantization of stochastic volatility models[END_REF] in the case case where both processes are discretized using an Euler-Maruyama scheme.

One can notice that building the quantization tree p p Y k q kP 0,n approximating pY t q tPr0,T s is a one dimensional problem as the diffusion of Y t is autonomous. Hence, based on our choice of discretization scheme, we will apply the fast recursive quantization (detailed above in Figure 9) that was introduced in [PS15] for one dimensional diffusion discretized by an Euler-Maruyama discretization scheme and then extended to higher order schemes, still in one dimension, by [START_REF] Thomas A Mcwalter | Recursive marginal quantization of higher-order schemes[END_REF]. The minor difference with existing literature is that, in our problem, the initial condition y 0 is not deterministic.

Then, using the quantization tree of p p Y k q kP 0,n we will be able to build the tree p p X k q kP 0,n following ideas developed in [FSP18, RMKP17, CFG18, CFG17]. Indeed, once the quantization tree of the volatility is built, we are in a one-dimensional setting and we are able to use fast deterministic algorithms.

Quantizing the volatility (a one-dimensional case)

Let pY t q tPr0,T s be a stochastic process in R and solution to the stochastic differential equation

dY t " r bpt, Y t qdt `r σpt, Y t qd Ă W t (3.11)
where Y 0 has the same law than the stationary measure ν: LpY 0 q " ν. In order to approximate our diffusion process, we choose a Milstein scheme for the time discretization, as defined in 3.4 and we build recursively the Markovian quantization tree p p Y t k q kP 0,n where p

Y t k`1 is the Voronoï quantization of r Y t k`1 defined by r Y t k`1 " M r b,r σ `tk , p Y t k , Z 2 k`1 ˘, p Y t k`1 " Proj Γ Y N 2,k`1 `r Y t k`1 ˘(3.12)
and the projection operator

Proj Γ Y N 2,k`1 p¨q is defined in (D.2), Γ Y N 2,k`1 " y k`1 1 , . . . , y k N 2,k`1
( is the grid of the optimal quantizer of r Y t k`1 and Z 2 k`1 " N p0, 1q. In order alleviate the notations, we will denote r

Y k and p Y k in place of r Y t k and p Y t k .
The first step consists in building p Y 0 , an optimal quantizer of size N 2,0 of Y 0 . Noticing that Y 0 " v ν 0 , we use the optimal quantizer we built for the pricing of European options. Then, we build recursively p p Y k q k"1,...,n , where the N 2,k -tuple are defined by y k 1:N 2,k " `yk 1 , . . . , y k N 2,k ˘, by solving iteratively the minimization problem defined in the Appendix D in (D.6), with the help of Lloyd's method I. Replacing X by r Y k`1 in (D.6) yields

y k`1 j " E " r Y k`1 1 Y k`1 P C j `ΓY N 2,k`1 ˘ı P ´r Y k`1 P C j `ΓY N 2,k`1 ˘" E " M r b,r σ `tk , p Y k , Z 2 k`1 ˘1M r b,r σ `tk , p Y k ,Z 2 k`1 ˘P C j `ΓY N 2,k`1 ˘ı P ´Mr b,r σ `tk , p Y k , Z 2 k`1 ˘P C j `ΓY N 2,k`1
˘¯.

(3.13) Now, preconditioning by p Y k in the numerator and the denominator and using

p k i " P `p Y k " y k i ˘,
we have

y k`1 j " E " E " M r b,r σ `tk , p Y k , Z 2 k`1 ˘1M r b,r σ `tk , p Y k ,Z 2 k`1 ˘P C j `ΓY N 2,k`1 ˘| p Y k ı  E " P ´Mr b,r σ `tk , p Y k , Z 2 k`1 ˘P C j `ΓY N 2,k`1 ˘| p Y k ¯ " N 2,k ÿ i"1 E " M r b,r σ `tk , y k i , Z 2 k`1 ˘1M r b,r σ `tk ,y k i ,Z 2 k`1 ˘P C j `ΓY N 2,k`1 ˘ı p k i N 2,k ÿ i"1 P ´Mr b,r σ `tk , y k i , Z 2 k`1 ˘P C j `ΓY N 2,k`1 ˘¯p k i " N 2,k ÿ i"1 ´Kk i `yk`1 j`1{2 ˘´K k i `yk`1 j´1{2 ˘¯p k i N 2,k ÿ i"1 ´F k i `yk`1 j`1{2 ˘´F k i `yk`1 j´1{2 ˘¯p k i (3.14) where C j `ΓY N 2,k`1 ˘" `yk`1 j´1{2 , y k`1 j`1{2
‰ is defined in (D.1). F k i and K k i are the cumulative distribution function and the first partial moment function of

U k i " µ k i `κk i pZ 1 k`1
`λk i q 2 respectively with

κ k j " pr σr σ 1 x qpt k , y k j qh 2 , λ k j " 1 ? hr σ 1 x pt k , y k j q
, and µ k j " y k j ´σpt k , y k j q 2r σ 1

x pt k , y k j q `hˆr bpt k , y k j q ´pr σr σ 1 x qpt k , y k j q 2 ˙.

(3.15)

The functions F k i and K k i can explicitly be determined in terms of the density and the cumulative distribution function of the normal distribution.

Lemma 3.2. Let U " µ `κpZ `λq 2 , with µ, κ, λ P R, λ ě 0, κ ą 0 and Z " N p0, 1q then the cumulative distribution function F X and the first partial moment K U of U are given by

F U pxq " `FZ px `q ´FZ px ´q˘1 xąµ K U pxq " ˆFU pxq `µ `κpλ 2 `1q ˘`κ ? 2π ´x´e ´x2 2 ´x`e ´x2 2 ¯˙1 xąµ (3.16) where x `" b x´µ κ ´λ, x
´" ´b x´µ κ ´λ and F Z is the cumulative distribution function of Z.

Finally, we can apply the Lloyd algorithm defined in Appendix D.9 with F X and K X defined by

F X pxq " N 2,k ÿ i"1 p k i F k i pxq and K X pxq " N 2,k ÿ i"1 p k i K k i pxq.
(3.17)

In order to be able to build recursively the tree quantization p p Y k q k"0,...,n , we need to have access to the weights p k i " P `p Y k " y k i ˘, which can be themselves computed recursively, as well as the conditional probabilities

p k ij " P `p Y k`1 " y k`1 j | p Y k " y k i ˘.
Lemma 3.3. The conditional probabilities p k ij are given by

p k ij "F k i `yk`1 j`1{2 ˘´F k i `yk`1 j´1{2 ˘. (3.18)
And the probabilities p k`1 j are given by

p k`1 j " N 2,k ÿ i"1 p k i p k ij . (3.19) Proof. The p k ij " P `p Y k`1 " y k`1 j | p Y k " y k i " P ´Mr b,r σ `tk , p Y k , Z 2 k`1 ˘P C j `ΓY N 2,k`1 ˘| p Y k " y k i " P ´Mr b,r σ `tk , y k i , Z 2 k`1 ˘P C j `ΓY N 2,k`1 ˘"F k i `yk`1 j`1{2 ˘´F k i `yk`1 j´1{2 ȃnd p k`1 j " P `p Y k`1 " y k`1 j ˘" N 2,k ÿ i"1 P `p Y k`1 " y k`1 j | p Y k " y k i ˘P `p Y k " y k i " N 2,k ÿ i"1 p k i p k ij .
As an illustration, we display in Figure 10 the rescaled grids obtained after recursive quantization of the boosted-volatility, where p v k " e ´κt k p Y k and p p Y k q k"1,...,n are the quantizers built using the fast recursive quantization approach.

Quantizing the asset (a one-dimensional case again)

Now, using the fact that pY t q t has already been quantized and the Euler-Maruyama scheme of pX t q t , as defined (3.8), we define the Markov quantized scheme

r X t k`1 " E b,σ `tk , p X t k , p Y t k , Z 1 k`1 ˘, p X t k`1 " Proj Γ X N 1,k`1 `r X t k`1 ˘(3.20)
where the projection operator where the N 1,k -tuple are defined by x k

Proj Γ X N 1,k`1 p¨q is defined in (D.2), Γ X N 1,k`1 is the optimal N 1,k`1 - quantizer of r X t k`1 and Z 1 k`1 " N p0,
1:N 1,k " `xk 1 , . . . , x k N 1,k ˘. Replacing X by r X k in (D.6) yield x k`1 j 1 " E " E b,σ `tk , p X t k , p Y t k , Z 1 k`1 ˘1E b,σ `tk , p Xt k , p Yt k ,Z 1 k`1 ˘P C j 1 `ΓX N 1,k`1 ˘ı P ´Eb,σ `tk , p X t k , p Y t k , Z 1 k`1 ˘P C j 1 `ΓX N 1,k`1 ˘" N 1,k ÿ i 1 "1 N 2,k ÿ i 2 "1 E " E b,σ `tk , x k i 1 , y k i 2 , Z 1 k`1 ˘1E b,σ `tk ,x k i 1 ,y k i 2 ,Z 1 k`1 ˘P C j 1 `ΓX N 1,k`1 ˘ı p k pi 1 ,i 2 q N 1,k ÿ i 1 "1 N 2,k ÿ i 2 "1 P ´Eb,σ `tk , x k i 1 , y k i 2 , Z 1 k`1 ˘P C j 1 `ΓX N 1,k`1 ˘¯p k pi 1 ,i 2 q " N 1,k ÿ i 1 "1 N 2,k ÿ i 2 "1 ´Kk pi 1 ,i 2 q `xk`1 j 1 `1{2 ˘´K k pi 1 ,i 2 q `xk`1 j 1 ´1{2 ˘¯p k pi 1 ,i 2 q N 1,k ÿ i 1 "1 N 2,k ÿ i 2 "1 ´F k pi 1 ,i 2 q `xk`1 j 1 `1{2 ˘´F k pi 1 ,i 2 q `xk`1 j 1 ´1{2 ˘¯p k pi 1 ,i 2 q (3.21)
where p k pi 1 ,i 2 q " P `p X k " x k i 1 , p Y k " y k i 2 ˘and F k pi 1 ,i 2 q and K k pi 1 ,i 2 q are the cumulative distribution function and the first partial moment function of the normal distribution

µ k pi 1 ,i 2 q `Z1 k`1 σ k pi 1 ,i 2 q
and they are defined by

F k pi 1 ,i 2 q pxq " F Z ˆx ´µk pi 1 ,i 2 q σ k pi 1 ,i 2 q Kk pi 1 ,i 2 q pxq " µ k pi 1 ,i 2 q F Z ˆx ´µk pi 1 ,i 2 q σ k pi 1 ,i 2 q ˙`σ k pi 1 ,i 2 q K Z ˆx ´µk pi 1 ,i 2 q σ k pi 1 ,i 2 q ˙(3.22) with µ k pi 1 ,i 2 q " x k i 1 `bpt k , x k i 1 , y k i 2 qh and σ k pi 1 ,i 2 q " σpt k , x k i 1 , y k i 2 q ? h (3.23)
and F Z and K Z are the cumulative distribution function and the first partial moment of the standard normal distribution. Finally, we apply the Lloyd method defined in Appendix (D.9) with F X and K X defined by

F X pxq " N 1,k ÿ i 1 "1 N 2,k ÿ i 2 "1 p k pi 1 ,i 2 q F k pi 1 ,i 2 q pxq and K X pxq " N 1,k ÿ i 1 "1 N 2,k ÿ i 2 "1 p k pi 1 ,i 2 q K k pi 1 ,i 2 q pxq. (3.24)
The sensitive part concerns the computation of the joint probabilities p k pi 1 ,i 2 q . Indeed, they are needed at each step in order to be able to design recursively the quantization tree.

Lemma 3.4. The joint probabilities p k pi 1 ,i 2 q are given by the following forward induction

p k`1 pj 1 ,j 2 q " N 1,k ÿ i"1 N 2,k ÿ j"1 p k pi 1 ,i 2 q P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 ˘(3.25)
where the joint conditional probabilities

P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 ˘are
given by the formulas below, depending on the correlation

• if CorrpZ 1 k`1 , Z 2 k`1 q " ρ " 0 P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 ˘" p k i 2 j 2 " N `xk i 1 ,i 2 ,j 1 ,`˘´N `xk i 1 ,i 2 ,j 1 ,´˘ı , (3.26
) where p k i 2 j 2 is defined in (3.18) and

x k i 1 ,i 2 ,j 1 ,´" x k`1 j 1 ´1{2 ´µk pi 1 ,i 2 q σ k pi 1 ,i 2 q , x k i 1 ,i 2 ,j 1 ,`" x k`1 j 1 `1{2 ´µk pi 1 ,i 2 q σ k pi 1 ,i 2 q , (3.27)
with µ k pi 1 ,i 2 q and σ k pi 1 ,i 2 q defined in (3.23).

• if CorrpZ 1 k`1 , Z 2 k`1 q " ρ ‰ 0 P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 " P ´Z1 k`1 P `xk i 1 ,i 2 ,j 1 ,´, x k i 1 ,i 2 ,j 1 ,`‰ , Z 2 k`1 P ´by k i 2 ,j 2 ,´´λ k i 2 , b y k i 2 ,j 2 ,`´λ k i 2 ıP ´Z1 k`1 P `xk i 1 ,i 2 ,j 1 ,´, x k i 1 ,i 2 ,j 1 ,`‰ , Z 2 k`1 P " ´by k i 2 ,j 2 ,`´λ k i 2 , ´by k i 2 ,j 2 ,´´λ k i 2 ¯( 3.28)
where

y k i 2 ,j 2 ,´" 0 _ y k`1 j 2 ´1{2 ´µk i 2 κ k i 2 , y k i 2 ,j 2 ,`" 0 _ y k`1 j 2 `1{2 ´µk i 2 κ k i 2 , (3.29)
with µ k i 2 , κ k i 2 and λ k i 2 defined in (3.15).

Remark 3.5. The probability in the right hand side of (3.28) can be computed using the cumulative distribution function of a correlated bivariate normal distribution2 . Indeed, let

F ρ px 1 , x 2 q " PpX 1 ď x 1 , X 2 ď x 2 q
the cumulative distribution function of the correlated centered Gaussian vector pX 1 , X 2 q with unit variance and correlation ρ, we have P `X1 P ra, bs, X 2 P rc, ds ˘" F ρ pb, dq ´Fρ pb, cq ´Fρ pa, dq `Fρ pa, cq (3.30) with a, c ě ´8 and b, d ď `8.

Proof.

p k`1 pj 1 ,j 2 q " P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 " N 1,k ÿ i"1 N 2,k ÿ j"1 P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 ˘P `p X k " x k i 1 , p Y k " y k i 2 " N 1,k ÿ i"1 N 2,k ÿ j"1 p k pi 1 ,i 2 q P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 ˘. • if CorrpZ 1 k`1 , Z 2 k`1 q " ρ " 0 P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 " p k i 2 j 2 P `p X k`1 " x k`1 j 1 | p X k " x k i 1 , p Y k " y k i 2 " p k i 2 j 2 P ´s X k`1 P `xk`1 j 1 ´1{2 , x k`1 j 1 `1{2 ‰ | p X k " x k i 1 , p Y k " y k i 2 " p k i 2 j 2 P ´Eb,σ `tk , x k i 1 , y k i 2 , Z 1 k`1 ˘P `xk`1 j 1 ´1{2 , x k`1 j 1 `1{2 ‰ " p k i 2 j 2 " N `xk i 1 ,i 2 ,j 1 ,`˘´N `xk i 1 ,i 2 ,j 1 ,´˘ı , • if CorrpZ 1 k`1 , Z 2 k`1 q " ρ ‰ 0 P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2 " P ´Eb,σ `tk , x k i 1 , y k i 2 , Z 1 k`1 ˘P `xk`1 j 1 ´1{2 , x k`1 j 1 `1{2 ‰ , M r b,r σ `tk , y k i 2 , Z 2 k`1 ˘P `yk`1 j 2 ´1{2 , y k`1 j 2 `1{2 ‰ " P ´µk pi 1 ,i 2 q `σk pi 1 ,i 2 q Z 1 k`1 P `xk`1 j 1 ´1{2 , x k`1 j 1 `1{2 ‰ , µ k i 2 `κk i 2 pZ 2 k`1 `λk i 2 q 2 P `yk`1 j 2 ´1{2 , y k`1 j 2 `1{2
‰ "

P ´Z1 k`1 P `xk i 1 ,i 2 ,j 1 ,´, x k i 1 ,i 2 ,j 1 ,`‰ , pZ 2 k`1 `λk i 2 q 2 P `yk i 2 ,j 2 ,´, y k i 2 ,j 2 ,`‰ " P ´Z1 k`1 P `xk i 1 ,i 2 ,j 1 ,´, x k i 1 ,i 2 ,j 1 ,`‰ , Z 2 k`1 P ´by k i 2 ,j 2 ,´´λ k i 2 , b y k i 2 ,j 2 ,`´λ k i 2 ıP ´Z1 k`1 P `xk i 1 ,i 2 ,j 1 ,´, x k i 1 ,i 2 ,j 1 ,`‰ , Z 2 k`1 P " ´by k i 2 ,j 2 ,`´λ k i 2 , ´by k i 2 ,j 2 ,´´λ k i 2 ¯¯.
Remark 3.6. Another possibility for the quantization of the Stationary Heston model could be to use optimal quantizers for the volatility at each date t k in place of using recursive quantization. Indeed, the volatility pv t q t being stationary and the fact that we required the volatility to start at time 0 from the invariant measure, we could use the grid of the optimal quantization p v 0 of size N of the stationary measure with its associated weights for every dates, hence setting p v k " p v 0 . We need as well the transitions from time t k to t k`1 defined by

P `p v k`1 " v k`1 j 2 | p v k " v k i 2 ˘. (3.31)
These probabilities can be computed using the conditional law of the CIR process described in [START_REF] John C Cox | A theory of the term structure of interest rates[END_REF][START_REF] Andersen | Efficient simulation of the heston stochastic volatility model[END_REF], which is a non-central chi-square distribution. Then, we would build the recursive quantizer of the log-asset at date p X k`1 with the standard methodology of recursive quantization using the already built quantizers of the volatility p v k and the log-asset p

X k at time t k , i.e. r X k`1 " E b,σ `tk , p X k , p v k , Z 1 k`1 ˘and p X k`1 " Proj Γ X N 1,k`1 `r X k`1 ˘(3.32)
where, this time, the Euler scheme is not defined in function of the boosted-volatility but directly in function of the volatility and is given by

E b,σ `t, x, v, z ˘" x `h´r ´q ´v 2 ¯`? v ? hz. (3.33)
However, the difficulties with this approach come from the computation of the couple transitions

P `p X k`1 " x k`1 j 1 , p v k`1 " v k`1 j 2 | p X k " x k i 1 , p v k " v k i 2 ˘. (3.34)
Indeed, these probability weights would not be as straightforward to compute as the methodology we adopt in this paper, namely using time-discretization schemes for both components. Our approach allows us to express the conditional probability of the couple as the probability that a correlated bi-variate Gaussian vector lies in a rectangle domain and this can be easily be computed numerically.

About the L 2 -error

In this part, we study the L 2 -error induced by the product recursive quantization approximation

p U k " p p X k , p Y k q of s U k " p s X k , s Y k q, the
time-discretized processes defined in (3.2) and (3.7) by

s U k " F k´1 p s U k´1 , Z k q (3.35)
where Z k " pZ 1 k , Z 2 k q is a standardized correlated Gaussian vector and the hybrid discretization scheme F k pu, Zq is given by

F k pu, Zq " ˜Eb,σ `tk , x, y, Z 1 k`1 Mr b,r σ `tk , y, Z 2 k`1 ˘¸. (3.36)
We recall the definition of the product recursive quantizer p

U k " p p X k , p Y k q. Its first component p X k is the projection of r X k onto Γ X N 1,k and the second component p Y k is the projection of r Y k onto Γ Y N 2,k , i.e., p X k`1 " Proj Γ X N 1,k`1 `r X k`1 ˘and p Y k`1 " Proj Γ Y N 2,k`1 `r Y k`1 ˘(3.37)
where r X k and r Y k are defined in (3.12) and (3.20), respectively. Moreover, if we consider the couple r U k " p r X k , r Y k q, using the above notations we have

r U k " F k´1 p p U k´1 , Z k q. (3.38)
It has been shown in [START_REF] Fiorin | Product markovian quantization of a diffusion process with applications to finance[END_REF][START_REF] Pagès | Weak and strong error analysis of recursive quantization: a general approach with an application to jump diffusions[END_REF] that if, for all k " 0, . . . , n ´1, the schemes F k pu, zq are Lipschitz in u, then there exists constants j " 1, . . . , n, C j ă `8 such that

} p U k ´s U k } 2 ď k ÿ j"1 C j `N1,j ˆN2,j ˘´1{2 (3.39)
where p U k and s U k are the processes defined in (3.37) and (3.38). The proof of this result is based on the extension of Pierce's lemma to the case of product quantization (see Lemma 2.3 in [START_REF] Pagès | Weak and strong error analysis of recursive quantization: a general approach with an application to jump diffusions[END_REF]).

In our case, the diffusion of the boosted volatility in the CIR model does not have Lipschitz drift and volatility components, hence the above result from [START_REF] Fiorin | Product markovian quantization of a diffusion process with applications to finance[END_REF][START_REF] Pagès | Weak and strong error analysis of recursive quantization: a general approach with an application to jump diffusions[END_REF] does not apply in our context. Even if we can hope to obtain similar results by applying the same kind of arguments, the results we obtain have to considered carefully. Indeed, when we take the limit in n Ñ `8, the number of time-step, the error upper-bound term goes to infinity. However, in practice, we consider h " kT {n fixed and then study the behavior of p U k in function of N 1,j and N 2,j for j ě k. The proof of the following proposition is given in Appendix C. Proposition 3.7. Let b, σ, r b and r σ, defined by (3.3) and (3.9), the coefficients of the log-asset and the boosted-volatility of the Heston model. Let, for every k " 0, . . . , n, p U k the hybrid recursive product quantizer at level N 1,k ˆN2,k of s U k . Then, for every k " 0, . . . , n (3.42)

} p U k ´s U k } 2 ď k ÿ j"0 r A j,k `N1,j ˆN2,j ˘´1{2 `Bk ? h (3.40) where r A j,k " 2 p´2 2p C 2 p A j,k ˆ2p p 2 ´1qj β j p } p U 0 } p 2 `αp 1 ´2p p 2 ´1qj β j p 1 ´2 p 2 ´1β p ˙1{p (
where ř H " 0 by convention and C T phq " Op1q.

Backward algorithm for Bermudan and Barrier options

Bermudan Options A Bermudan option is a financial derivative product that gives the right to its owner to buy or sell (or to enter to, in the case of a swap) an underlying product with a given payoff ψ t p¨, ¨q at predefined exercise dates tt 0 , ¨¨¨, t n u. Its price, at time t 0 " 0, is given by

sup τ Ptt 0 ,¨¨¨,tnu E " e ´rτ ψ τ pX τ , Y τ q | F t 0 ı
where X t and Y t are solutions to the system defined in (3.10).

In this part, we follow the numerical solution first introduced by [START_REF] Bally | A quantization tree method for pricing and hedging multi-dimensional american options[END_REF][START_REF] Bally | A quantization algorithm for solving multidimensional discrete-time optimal stopping problems[END_REF]. They proposed to solve discrete-time optimal stopping problems using a quantization tree of the risk factors X t and Y t .

Let F X,Y " pFq 0ďkďn the natural filtration of X and Y . Hence, we can define recursively the sequence of random variable

L p -integrable pV k q 0ďkďn # V n " e ´rtn ψ n pX n , Y n q, V k " max `e´rt k ψ k pX k , Y k q, ErV k`1 | F k s ˘, 0 ď k ď n ´1 (3.43) called Backward Dynamic Programming Principle. Then V 0 " sup Ere ´rτ ψ τ pX τ , Y τ q | F 0 s, τ P Θ 0,n (
with Θ 0,n the set of all stopping times taking values in tt 0 , ¨¨¨, t n u. The sequence pV k q 0ďkďn is also known as the Snell envelope of the obstacle process `e´rt k ψ k pX k , Y k q ˘0ďkďn . In the end, ErV 0 s is the quantity we are interested in. Indeed, ErV 0 s is the price of the Bermudan option whose payoff is ψ k and is exercisable at dates tt 1 , ¨¨¨, t n u.

Following what was defined in (3.43), in order to compute ErV 0 s, we will need to use the previously defined quantizer of X k and Y k : p X k and p Y k . Hence, for a given global budget N " N 1,0 N 2,0 `¨¨¨`N 1,n N 2,n , the total number of nodes of the tree by the couple p p X k , p Y k q 0ďkďn , we can approximate the Backward Dynamic Programming Principle (3.43) by the following sequence involving the couple p p

X k , p Y k q 0ďkďn # p V n " e ´rtn ψ n p p X n , p Y n q, p V k " max `e´rt k ψ k p p X k , p Y k q, Er p V k`1 | p p X k , p Y k qs ˘, k " 0, . . . , n ´1.
(3.44)

Remark 3.8. A direct consequence of choosing recursive Markovian Quantization to spatially discretize the problem is that the sequence p p X k , p Y k q 0ďkďn is Markovian. Hence p p V k q 0ďkďn defined in (3.44) obeying a Backward Dynamic Programming Principle is the Snell envelope of

`e´rt k ψ k p p X k , p Y k q ˘0ďkďn
. This is the main difference with the first approach of [START_REF] Bally | A quantization tree method for pricing and hedging multi-dimensional american options[END_REF][START_REF] Bally | A quantization algorithm for solving multidimensional discrete-time optimal stopping problems[END_REF],

where in there case they only had a pseudo-Snell envelope of

`e´rt k ψ k p p X k , p Y k q ˘0ďkďn .
Using the discrete feature of the quantizers, (3.44) can be rewritten

$ ' ' ' ' ' & ' ' ' ' ' % p v n px n i 1 , y n i 2 q " e ´rtn ψ n px n i 1 , y n i 2 q, i 1 " 1, . . . , N 1,n i 2 " 1, . . . , N 2,n p v k px k i 1 , y k i 2 q " max ´e´rt k ψ k px k i 1 , y k i 2 q, N 1,k`1 ÿ j 1 "1 N 2,k`1 ÿ j 2 "1 π k pi 1 ,i 2 q,pj 1 ,j 2 q p v k`1 px k`1 j 1 , y k`1 j 2 q ¯, k " 0, . . . , n ´1 i 1 " 1, . . . , N 1,k i 2 " 1, . . . , N 2,k (3.45) where π k pi 1 ,i 2 q,pj 1 ,j 2 q " P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2
˘is the conditional probability weight given in (3.28). Finally, the approximation of the price of the Bermudan option is given by

E " p v 0 px 0 , p Y 0 q ‰ " N 2,0 ÿ i"1 p i p v 0 px 0 , y 0 i q (3.46)
with p i " P `p Y 0 " y 0 i ˘given by (2.11).

Barrier Options A Barrier option is a path-dependent financial product whose payoff at maturity date T depends on the value of the process X T at time T and its maximum or minimum over the period r0, T s. More precisely, we are interested by options with the following types of payoff h h " f pX T q 1 tsup tPr0,T s XtPIu or h " f pX T q 1 tinf tPr0,T s XtPIu (3.47)

where I is an unbounded interval of R, usually of the forme p´8, Ls or rL, `8q (L is the barrier) and f can be any vanilla payoff function (Call, Put, Spread, Butterfly, ...).

In this part, we follow the methodology initiated in [START_REF] Sagna | Pricing of barrier options by marginal functional quantization[END_REF] in the case of functional quantization. This work is based on the Brownian bridge method applied to the Euler-Maruyama scheme as described e.g. in [START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF]. We generalize it to stochastic volatility models and product Markovian recursive quantization. X t being discretized by an Euler-Maruyama scheme, yielding s X k with k " 0, . . . , n, we can determine the law of max tPr0,T s s X t and min tPr0,T s s X t given the values s

X k " x k , s Y k " y k , k " 0, . . . , n L ´max tPr0,T s s X t | s X k " x k , s Y k " y k , k " 0, . . . , n ¯" L ´max k"0,...,n´1 `Gk px k ,y k q,x k`1 ˘´1 pU k q ¯(3.48) and L ´min tPr0,T s s X t | s X k " x k , s Y k " y k , k " 0, . . . , n ¯" L ´max k"0,...,n´1 `F k px k ,y k q,x k`1 ˘´1 pU k q ¯(3.49)
where pU k q k"0,...,n´1 are i.i.d uniformly distributed random variables over the unit interval and pG k px,yq,z q ´1 and pF k px,yq,z q ´1 are the inverse of the conditional distribution functions G k px,yq,z and F k px,yq,z defined by

G k px,yq,z puq " ´1 ´e´2n px´uqpz´uq T σ 2 pt k ,x,yq ¯1tuěmaxpx,zqu (3.50)
and

F k px,yq,z puq " 1 ´´1 ´e´2n px´uqpz´uq T σ 2 pt k ,x,yq ¯1tuďminpx,zqu . (3.51)
Now, using the resulting representation formula for E f p s X T , max tPr0,T s s X t q (see e.g. [Sag10, Pag18]), we have a new representation formula for the price of up-and-out options s P U O and down-and-out options s

P DO s P U O " e ´rT E " f p s X T q 1 sup tPr0,T s s XtďL ‰ " e ´rT E " f p s X T q n´1 ź k"0 G k pX k ,Y k q, s X k`1 pLq  (3.52)
and

s P DO " e ´rT E " f p s X T q 1 inf tPr0,T s s XtěL ‰ " e ´rT E " f p s X T q n´1 ź k"0 ´1 ´F k p s X k , s Y k q, s X k`1 pLq ¯ (3.53)
where L is the barrier.

# p V n " e ´rT f p p X n q, p V k " E " g k p p X k , p Y k , p X k`1 q p V k`1 | p p X k , p Y k q ‰ , 0 ď k ď n ´1 (3.54) that can be rewritten $ ' ' ' ' ' & ' ' ' ' ' % p v n px n i 1 , y n i 2 q " e ´rT f px n i q, i " 1, . . . , N 1,n j " 1, . . . , N 2,n p v k px k i 1 , y n i 2 q " N 1,k`1 ÿ j 1 "1 N 2,k`1 ÿ j 2 "1 π k pi 1 ,i 2 q,pj 1 ,j 2 q p v k`1 px k`1 j 1 , y k`1 j 2 qg k px k i 1 , y k i 2 , x k`1 j 1 q, k " 0, . . . , n ´1 i " 1, . . . , N 1,k j " 1, . . . , N 2,k (3.55) with π k pi 1 ,i 2 q,pj 1 ,j 2 q " P `p X k`1 " x k`1 j 1 , p Y k`1 " y k`1 j 2 | p X k " x k i 1 , p Y k " y k i 2
˘the conditional probabilities given in (3.28) and g k px, y, zq is either equal to G k px,yq,z pLq or 1 ´F k px,yq,z pLq depending on the option type. Finally, the approximation of the price of the barrier option is given by

Er p V 0 s " E " p v 0 px 0 , p Y 0 q ‰ " N 2,0 ÿ i"1 p i p v 0 px 0 , y 0 i q (3.56)
with p i " P `p Y 0 " y 0 i ˘given by (2.11).

Numerical illustrations

In this part, we deal with numerical experiments in the Stationary Heston model. We will apply the methodology based on hybrid product recursive quantization to the pricing of European, Bermudan and Barrier options. For the model parameters, we consider the parameters given in Table 2 obtained after the penalized calibration procedure and instead of considering the market value for S 0 , we take S 0 " 100 in order to get prices of an order we are used to. For the size of the quantization grids, we consider grids of constant size for all time-steps: for all k " 0, . . . , n, we take N 1,k " N 1 and N 2,k " N 2 where n is the number of time steps. During the numerical tests, we vary the tuple values pn, N 1 , N 2 q.

All the numerical tests have been carried out in C++ on a laptop with a 2,4 GHz 8-Core Intel Core i9 CPU. The computations of the transition probabilities are parallelized on the CPU.

European options First, we compare, in Table 3, the price of European options with maturity t n " T " 0.5 (6 months) computed using the quantization tree to the benchmark price computed using the methodology based on the quadrature formula (the quadrature formula with Laguerre polynomials) explained in Section 2. In place of using the backward algorithm (3.44) (without the function max) for computing the expectation at the expiry date, we use the weights p k pi 1 ,i 2 q defined in (3.25) and built by forward induction, in order to compute

E " e ´rtn ψ n p p X n , p Y n q ‰ " e ´rtn N 1,n ÿ i 1 "1 N 2,n ÿ i 2 "1 ψ n px n i 1 , y n i 2 q. (3.57)
We give, in parenthesis, the relative error induced by the quantization-based approximation. We compare the behavior of the pricers with different size of grids and numbers of discretization steps. We notice that the main part of the error is explained by the size of the time-step n. Table 4: Comparison between European options prices, with maturity T " 0.5 (6 months), given by quantization and the benchmark, in function of the strike K and of the size n where we set pN 1 , N 2 q " p50, 10q.

pN 1 , N 2 q K Benchmark p20, 5q p50
Bermudan options Then, in Figure 11, we display the prices of monthly exercisable Bermudan options with maturity T " 0.5 (6 months) for Call and Put of strikes K " 100. The prices are computed by quantization and we compare the behavior of the pricer for different choices of time-step n and sizes of the asset grids N 1 where we set N 2 " 10. Again, we notice that the choice of n has a high impact on the price given by quantization compared to the choice of the grid size. Barrier options Finally, in Figure 12, we display the prices of an up-and-out Barrier option with strike K " 100, maturity T " 0.5 (6 months), barrier L " 115 and N 2 " 10 computed with quantization. Again, we can notice the impact of n on the approximated price.

Appendices

A Discretization scheme for the volatility preserving the positivity

We recall the dynamics of the volatility

dv t " κpθ ´vt qdt `ξ? v t d Ă W t
with κ ą 0, θ ą 0 and ξ ą 0. In this section, we discuss the choice of the discretization scheme under the Feller condition, which ensures the positivity of the process.

Euler-Maruyama scheme. Discretizing the volatility using an Euler-Maruyama scheme

s v t k`1 " s v t k `κpθ ´s v t k qh `ξa s v t k ? h Z 2 k`1 with t k " kh, h " T {n and Z 2 k`1 " p Ă W t k`1 ´Ă W t k q{ ?
h may look natural. However, such a scheme clearly does not preserve positivity of the process even if the Feller condition is fulfilled since

P `s v t 1 ă 0 ˘" P ˆZ ă ´v0 ´κpθ ´v0 qh ξ ? v 0 ? h ˙ą 0
with Z " N p0, 1q. This suggests to introduce the Milstein scheme which is quite tractable in one dimension in absence of Lévy areas.

Milstein scheme. The Milstein scheme of the stochastic volatility is given by

s v t k`1 " M b,σ `tk , s v t k`1 , Z 2 k`1
where

(see (3.4)) M b,σ pt, x, zq " x ´σpxq 2σ 1 x pxq `hˆb pt, xq ´pσσ 1 x qpxq 2 ˙`pσσ 1 x qpxqh 2 ˆz `1 ? hσ 1 x pxq ˙2.
with bpxq " κpθ ´xq, σpxq " ξ ?

x and σ 1 x pxq " ξ 2 ?

x . Consequently, under the Feller condition, the positivity of M b,σ pt, x, zq is ensured if

x ě σpxq 2σ 1 x pxq ě 0, bpt, xq ě pσσ 1 x qpxq 2 ě 0.
In The last inequality is satisfied thanks to the condition ξ 2 2κθ ď 1 ensuring the positivity of the scheme.

B L p -linear growth of the hybrid scheme

The aim of this section is to show the L p -linear growth of the scheme F k pu, zq with u " px, yq defined by

F k pu, Zq " ˜Eb,σ `tk , x, y, Z 1 k`1 Mr b,r σ `tk , y, Z 2 k`1 ˘¸. (B.1)
where the schemes E b,σ and M r b,r σ are defined in (3.8) and (3.4), respectively. The results on the L p -linear growth of the schemes are essentially based on the key Lemma 2.1 proved in [START_REF] Pagès | Weak and strong error analysis of recursive quantization: a general approach with an application to jump diffusions[END_REF]

Z 2 k`1 pZ 2 k`1 q 2 ´1'
.

First, we show the linear growth of apuq C Proof of the L 2 -error estimation of Proposition 3.7

|apuq| " ´ˇˇx `h`r ´e´κt k y 2 ˘ˇˇ2 `ˇy `eκt k κθh ˇˇ2 ¯1{2 " ´|x| 2 `|y| 2 `h2 ´r2 `e´2κt k 4 |y| 2 ¯`e 2κt k κ 2 θ 2 h 2 ¯1{2 ď ´|u| 2 ´1 `h2 e ´2κt k 4 ¯`h 2 `r2 `e2κt k κ 2 θ 2 ˘¯1{2 ď |u| ´1 `h2 e ´2κt k 4 ¯1{2 `h`r 2 `e2κt k κ 2 θ 2 ˘1{2 ď |u| ´1 `h h 2 ¯`h `r2 `e2κT κ 2 θ 2 ˘1{2 ď |u|p1 `Lhq `Lh where L " max ´1 2 , `r2 `e2κT κ 2 θ 2 ˘1{2 ¯.
We have, for every k " 0, . . . , n ´1

p U k`1 ´s U k`1 " p U k`1 ´r U k`1 `r U k`1 ´s U k`1 " p U k`1 ´r U k`1 `Fk p p U k , Z k`1 q ´Fk p s U k , Z k`1 q (C.1)
by the very definition of r U k`1 and s U k`1 . Hence,

} p U k`1 ´s U k`1 } 2 ď } p U k`1 ´r U k`1 } 2 `} r U k`1 ´s U k`1 } 2 ď } p U k`1 ´r U k`1 } 2 `}F k p p U k , Z k`1 q ´Fk p s U k , Z k`1 q} 2 . (C.2)
Using the definition of Milstein scheme of the boosted -volatility models M r b,r σ in (3.5), the 1 2 -Hölder property of ?

x, for every y, y 1 P R `one has ˇˇM r b,r σ `t, y, z ˘´M r b,r σ `t, y ´?hq. Wo choose λ of this order because we wish to divide equally the impact of h and get ? h on each side. Hence, we have (C.9) and ř H " 0 by convention. Now, we follow the lines of the proof developed in [START_REF] Pagès | Weak and strong error analysis of recursive quantization: a general approach with an application to jump diffusions[END_REF], we apply the revisited Pierce's lemma for product quantization (Lemma 2.3 in [START_REF] Pagès | Weak and strong error analysis of recursive quantization: a general approach with an application to jump diffusions[END_REF]) with r " 2 and let p ą r " 2, which yields

› › F k p p U k , Z k`1 q ´Fk p s U k , Z k`1 q › › 2 ď › › E b,σ `tk , p X k , p Y k , Z 1 k`1 ˘´E b,σ `tk , s X k , s Y k , Z 1 k`1 ˘› › 2 `› › M r b,r σ `tk , p Y k , Z 2 k`1 ˘´M r b,r σ `tk , s Y k , Z 2 k`1 ˘› › 2 ď } p X k ´s X k } 2 `´1 `e´κt k 2 h ¯} p Y k ´s Y k } 2 `› › › ? h `ξ e κt k {2 `e´κt k {2 ˘b| p Y k ´s Y k | › › › 2 ď } p X k ´s X k } 2 `´1 `e´κt k 2 h ¯} p Y k ´s Y k } 2 `› › › a 2hpξ 2 e κt k `e´κt k q b | p Y k ´s Y k | › › › 2 . (C.
a 2hpξ 2 e κt k `e´κt k q b | p Y k ´s Y k | ď 1 2 ˆ2hpξ 2 e κt k `e´κt k q λ `| p Y k ´s Y k |λ ˙. (C.6) Then, › › F k p p U k , Z k`1 q ´Fk p s U k , Z k`1 q › › 2 ď } p X k ´s X k } 2 `´1 `e´κt k 2 h ¯} p Y k ´s Y k } 2 `› › › 1 2 ˆ2hpξ 2 e κt k `e´κt k q λ `| p Y k ´s Y k |λ ˙› › › 2 ď } p X k ´s X k } 2 `´1 `e´κt k 2 h `λ 2 ¯} p Y k ´s Y k } 2 `
} p U k`1 ´s U k`1 } 2 ď 2 p´2 2p C p k`1 ÿ j"0 A j,k`1 } r U j } p `N1,j ˆN2,j ˘´1{2 `Bk`1 ? h (C.10)
where C p " 2C 1,p and C 1,p is the constant appearing in Pierce lemma (see the second item in Theorem D.8 and [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF] for further details) and we used that } r U j } p ě σ p p r U j q " inf aPR 2 } r U j ´a} p . Moreover, noting that the hybrid discretization scheme F k has an L p -linear growth, (see Appendix B), i.e. @k " 0, . . . , n ´1, @u P R Furthermore, E " | p U j | p ‰ can be upper-bounded using Jensen's inequality and the stationary property satisfied by p X j and p Y j independently. Indeed, they are one-dimensional quadratic optimal quantizers of r X j and r Y j , respectively, hence they are stationary in the sense of Proposition D.6.

} p U j } p p ď 2 p 2 ´1´E " | p X j | p ‰ `E " | p Y j | p ‰ ď 2 p 2 ´1ˆE " ˇˇE " r X j | p X j ‰ˇˇp ı `E " ˇˇE " r Y j | p Y j ‰ˇˇp ı ď 2 p 2 ´1´E " | r X j | p ‰ `E " | r Y j | p ‰ " 2 p 2 ´1} r U j } p p ď 2 p 2 ´1} r U j } p 2 .
(C.13) (a) Sharp rate [START_REF] Zador | Asymptotic quantization error of continuous signals and the quantization dimension[END_REF][START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]. Let X P L p`δ R pPq for some δ ą 0. Let P X pdξq " ϕpξq λpdξq `νpdξq, where ν K λ i.e., is singular with respect to the Lebesgue measure λ on R. Then, there is a constant r J p,1 P p0, `8q such that where, for r P p0, `8q, σ r pXq " min aPR }X ´a} r ă `8 is the L r -pseudo-standard deviation. Now, we will be interested by the construction of such quadratic optimal quantizer. We differentiate Q 2,N , whose gradient is given by where K X p¨q and F X p¨q are the first partial moment and the cumulative distribution respectively, function of X, i.e. K X pxq " E " X 1 Xďx ‰ and F X pxq " P `X ď x ˘. (D.7)

∇ Q 2,
Hence, one can notices that the optimal quantizer that cancel the gradient defined in (D.6), hence is an optimal quantizer, is a stationary quantizer in the following sense

E " p x N | X ‰ " p X N . (D.8)
The last equality in (D.6) was the starting point to the development of the first method devoted to the numerical computation of optimal quantizers: the Lloyd's method I. This method was first devised in 1957 by S.P. Lloyd and published later [START_REF] Lloyd | Least squares quantization in pcm[END_REF]. Starting from a sorted Ntuple x r0s and with the knowledge of the first partial moment K X and the cumulative distribution function F X of X, the algorithm, which is essentially a deterministic fixed point method, is defined as follows In the seminal paper of [START_REF] John | Exponential rate of convergence for lloyd's method i[END_REF], it has been shown that `xrns ˘ně1 converges exponentially fast toward x, the optimal quantizer, when the density ϕ of X is log-concave and not piecewise affine. Numerical optimizations can be made in order to increase the rate of convergence to the optimal quantizer such as fixed point search acceleration, for example the Anderson acceleration (see [START_REF] Donald G Anderson | Iterative procedures for nonlinear integral equations[END_REF] for the original paper and [START_REF] Homer | Anderson acceleration for fixed-point iterations[END_REF] for details on the procedure).

Of course, other algorithms exist, such as the Newton Raphson zero search procedure or its variant the Levenberg-Marquardt algorithm which are deterministic procedures as well if the density, the first partial moment and the cumulative distribution function of X are known. Additionally, we can cite stochastic procedures such as the CLVQ procedure (Competitive Learning Vector Quantization) which is a zero search stochastic gradient and the randomized version of the Lloyd's method I. For more details, the reader can refer to [START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF][START_REF] Pagès | Pointwise convergence of the lloyd i algorithm in higher dimension[END_REF].

Once the algorithm (D.9) has been converging, we have at hand the quadratic optimal quantizer p X N of X and its associated probabilities given by

P `p X N " x n i ˘" F X `xN i`1{2 ˘´F X `xN i´1{2 ˘, i "

Figure 1 :

 1 Figure 1: Implied volatility surface of the Euro Stoxx 50 as of the 26th of September 2019.(S 0 " 3541, r " ´0.0032 and q " 0.00225) The expiries T are given in days and the strikes K in percentage of the spot.

Figure 2 :

 2 Figure 2: Implied volatilities for 22 (left) and 50 (right) days expiry options after calibration at 50 days without penalization.

Figure 3 :

 3 Figure 3: Implied volatilities for 7 (left) and 14 (right) days expiry options after calibration at 50 days without penalization.

Figure 4 :

 4 Figure 4: Term-structure of the volatility in function of T and K of both models (left: Standard Heston and right: Stationary Heston) after calibration at 50 days without penalization.

Figure 5 :

 5 Figure 5: pK, T q ÝÑ |σ M arket iv

Figure 6 :

 6 Figure 6: Implied volatilities for 22 (left) and 50 (right) days expiry options after calibration at 50 days with penalization.

Figure 7 :

 7 Figure 7: Implied volatilities for 7 (left) and 14 (right) days expiry options after calibration at 50 days with penalization.

Figure 8 :

 8 Figure 8: pK, T q ÝÑ |σ M arket iv

Figure 11 :

 11 Figure 11: Prices of Bermudan options in the stationary Heston model given by product hybrid recursive quantization with fixed value N 2 " 10.

Figure 12 :

 12 Figure12: Prices of Barrier options with strike K " 100 in the stationary Heston model given by product hybrid recursive quantization with fixed value N 2 " 10.

Γ

  N ĂR,|Γ N |ďN }X ´p X N } p " Non asymptotic upper-bound[START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF][START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF]. Let δ ą 0. There exists a real constant C 1,p P p0, `8q such that, for every R-valued random variable X, @N ě 1, minΓ N ĂR,|Γ N |ďN }X ´p X N } p ď C 1,p σ δ`p pXqN ´1 (D.4)

  1q. Again, in order to simplify the notations, r X t k and p X t k are denoted in what follows by r X k and p X k . Note that we are still in an one-dimensional case, hence we can apply the same methodology as developed in Appendix D and build recursively the quantization `p X k ˘k"0,...,n as detailed above,

						0.30
						0.10 0.15 0.20 0.25	Probab ility p k i
						0.05
						0.00
						50 60
						40
						10 20 30 Time
	0.00	0.05	Quan tizer vk 0.10 0.15	0.20	0.25	0

Figure 10: Rescaled Recursive quantization of the boosted-volatility process with its associated weights from t " 0 to t " 60 days with a time step of 5 days with grids of size N " 10. The recursive quantization methodology is applied to p Y k and then we display the rescaled volatility p v k " e ´κt k p Y k .

Table 3 :

 3 Comparison between European options prices, with maturity T " 0.5 (6 months), given by quantization and the benchmark, in function of the strike K and pN 1 , N 2 q where we set n " 180.

	, 10q	p100, 10q	p150, 10q

  our case, if the first condition holds true since In order to solve this problem, we consider the following boosted volatility process Y Milstein scheme for the boosted volatility. Let Y t " e κt v t , t P r0, T s for some κ ą 0, which satisfies, owing to Itô's formuladY t " e κt κθdt `ξ e κt{2 a Y t d Ă W t .Remark A.1. The process pY t q tPr0,T s will have a higher variance but, having in mind a quantized scheme, this has no real impact (by contrast with a Monte Carlo simulation).

	σpxq 2σ 1 x pxq pσσ 1 x qpxq 2 " " ξ s Y t k`1 " M r b,r σ `tk , s ξ ? x 2 ξ 2 ? x ? x ξ 2 ? 2 Y t k , Z 2 " x x " k`1 ξ 2 4 ȗsing the notation defined in (3.4) where drift and volatility terms of the boosted process, now the second one fails. Indeed time-dependents, are given by r bpt, xq " e κt κθ, r σpt, xq " ξ ? x e κt{2 and r σ 1 x pt, xq " ξ e κt{2 2 ? x . Under the Feller condition, the positivity of the scheme is ensured, since can be bigger than bpt, xq. Now, if we look at the Milstein discretization scheme of Y t r σpt, xq 2r σ 1 x pt, xq " x and pr σr σ 1 x qpt, xq 2 " ξ 2 e

t " e κt v t , t P r0, T s. (A.1) κt 4 ď r bpt, xq " e κt κθ.

  in R d that we recall below. Lemma B.1. (a) Let u P R d and Apuq be a d ˆq-matrix and let apuq P R d . Let p P r2, 3q.For any centered random vector ζ P L p R d pΩ, A, Pq, one has for every h P p0, `8qNow, we will apply Lemma B.1 to F k pu, zq defined in (B.1) further on in order to show its L p -linear growth. Let apuq P R 2 and let Apuq be a 2 ˆ3-matrix defined by

			apuq "	˜x	`h`r ´e´κt k y 2 `eκt k κθh	y ¸, Apuq "	y ˜e´κt k {2 ? 0	0 y e κt k {2 ? ?	h ξ 2 e κt k 4 0	and
											¨Z1
											k`1
											ζ "
	E	"	|apuq `?hApuqζ| p ‰	ď ˆ1	`pp ´1qpp ´2q 2	h ˙|apuq| p `h`1 `p	`h p 2	´1˘} Apuq} p E	"	|ζ| p ‰
											(B.2)
	where }Apuq} " `TrpApuqA ‹ puqq ˘1{2 .
	(b) In particular, if |apuq| ď |u|p1 `Lhq `Lh and }Apuq} p ď 2 p´1 Υ p p1 `|u| p q, then
						E	"	|apuq `?hApuqζ| p ‰	ď `eκph L `Kp ˘h ``e κph `Kp h ˘|u| p ,	(B.3)
	where
			κ p "	pp ´1qpp ´2q 2	`2pL and K p " 2 p´1 Υ p `1 `p	`h p 2	´1˘E "	|ζ| p ‰	.	(B.4)

  Then, we study }Apuq} p }Apuq} p " ´e´κt k |y| `|y| e κt k `h ξ 4 e 2κt k 16 Z k`1 q| p ‰ ď α p `βp |u| p with α p " `eκph L `Kp ˘h and β p " e κph `Kp h (B.5) where K p and κ p are defined in the Lemma B.1.

											¯p{2
			" ´|y|pe ´κt k `eκt k q	16 `h ξ 4 e 2κt k	¯p{2
			ď 2	p 2	´1´|	y|	p 2 pe ´κt k `eκt k q	p 2 `h p 2	4 p ξ 2p e pκt k	ď
			2	p 2	´1´| y| p `1 2	pe ´κt k `eκt k q	p 2 `h p 2	4 p ξ 2p e pκt k	ď
			2	p 2	´1 p1 `eκT q 2	p 2	´|y| p `1	`h p 2	2 2p´1 ξ 2p e pκT	p1 `eκT q 1	2 p	ď
			2 p´1 Υ p `1 `|u| p where
	Υ p " p1`e κT q 2 L p -linear growth	p 2	`h p 2 ξ 2p e pκT 2 2p . Hence, by Lemma B.1, the discretization scheme F k has an
					E	" |F k pu,

  1 , z ˘ˇ" ˇˇˇ´z ξ e κt{2 ? Maruyama scheme of the log-asset E b,σ defined in (3.8) we have, for any x, x 1 , y, y 1 P R

	h ˇˇ`| z|ξ e κt{2 ? `?y ¯2 ´´z h|z|ξ e κt{2 `|y ´y1 | ξ e κt{2 ? h 2 ? ˇˇ? y ´ay 1 2 a |y ´y1 | h `?y `ay 1	¯2ˇˇˇď ď `ay 1 (C.3)
	and using the definition of the Euler-Èb,σ
	`t, x, y, z ˘´E b,σ `t, x 1 , y 1 , z ˘ˇď |x ´x1 |	`e´κt 2	h|y ´y1 | `e´κt{2 ?	h|z| a	|y ´y1 |.	(C.4)
	Now, when we replace x, y, x 1 , y 1 by p X k , p Y k , s X k , s Y k in the last expression, we get an upper-
	bound for the last term of (C.2)

  pξ 2 e κt k `e´κt k q T phq " p1 `ξ2 e κT qp1 ´?hq ´1 " Op1q. Finally, (C.2) is upper-bounded by} p U k`1 ´s U k`1 } 2 ď } p U k`1 ´r U k`1 } 2 `?2 ´1

	where	ď ď	? ?	h λ U k ´s U k } 2 `CT phq ¯} p ? h ¯k´j`1 `?hC T phq k ÿ j"0 2 k´j 2 ´1 ? 2 `?h h λ ´1 U k ´s U k } 2 `CT ¯} p ¯} p `h 2 `λ 2 `?h 2 ď k`1 ÿ j"0 } p U j ´r U j } 2 2 k´j`1 2 2 ď k`1 ÿ j"0 A j,k`1 } p U j ´r U j } 2 `Bk`1 h j"0 U `?h A j,k " 2 2 ´1 2 ´1 k´j 2 e ? h 2 pk´jq and B k " C T phq k´1 ÿ 2 k´1´j 2 e ? h 2 pk´1´jq	`?h 2	¯k´j (C.8)

k ´s U k } 2 `CT phq ? h (C.7)

where C

  2 , E " |F k pu, Z k`1 q| p ‰ ď α p `βp |x| p , (C.11)where the coefficients α p and β p are defined in (B.5). Hence, for all j " 0, . . . , n ´1, we have

	} r U j`1 } p p " E	"	E	"	|F j p p U j , Z j`1 q| p | p U j	‰‰	ď α p `βp } p U j } p p .	(C.12)

  N pxq "

		ˆE " px N i ´Xq 1 X P `xN i´1{2 ,x N i`1{2	‰	˙i"1,...,N ı	.	(D.5)
	Moreover, if x is solution to the distortion minimization problem then it satisfies
			"	ı
	∇ Q 2,N pxq " 0	ðñ x N i "	E P ´X P `xN X 1 X P `xN i´1{2 ,x N i`1{2 i´1{2 , x N i`1{2 ‰ ‰ ¯,		i " 1, . . . , N	(D.6)
		ðñ x N i "	K X	`xN i`1{2 ˘´K X `xN i`1{2 ˘´F X `xN `xN i´1{2 i´1{2 ˘, FX i " 1, . . . , N

  1, . . . , n. (D.10)

In practice, we choose n "

This number of points allows us to reach a high precision while keeping the computation time under control.

C++ implementation of the upper right tail of a bivariate normal distribution can be found in John Burkardt's website https://people.sc.fsu.edu/~jburkardt/cpp_src/toms462/toms462.html.

Acknowledgment

The authors wish to thank Guillaume Aubert for fruitful discussion on the Heston model and Jean-Michel Fayolle for his advice on the calibration of the models. PhD thesis of Thibaut Montes is funded by a CIFRE grand from The Independent Calculation Agent (The ICA) and French ANRT.

Let X be a R-valued random variable with distribution P X defined on a probability space pΩ, A, Pq such that X P L 2 R pΩ, A, Pq. Definition D.1. Let Γ N " tx N 1 , . . . , x N N u Ă R be a subset of size N , called N -quantizer. A Borel partition pC i pΓ N qq iPt1,...,N u of R is a Voronoï partition of R induced by the N -quantizer Γ N if, for every i P t1, . . . , N u,

The Borel sets C i pΓ N q are called Voronoï cells of the partition induced by Γ N .

Remark D.2. Any such N -quantizer is in correspondence with the N -tuple x " px N 1 , . . . , x N N q P pRq N as well as with all N -tuples obtained by a permutation of the components of x. This is why we will sometimes replace Γ N by x.

If the quantizers are in non-decreasing order:

then the Voronoï cells are given by

where @i PP t2, . . . , N u, x N i´1{2 "

and its associated probabilities, also called weights, are given by

Definition D.4. The quadratic distortion function at level N induced by an N -tuple x " px N 1 , . . . , x N N q is given by

Of course, the above result can be extended to the L p case by considering the L p -mean quantization error in place of the quadratic one.

We briefly recall some classical theoretical results, see [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF][START_REF] Pagès | Numerical Probability: An Introduction with Applications to Finance[END_REF] for further details. The first one treats of existence of optimal quantizers. Theorem D.5. (Existence of optimal N -quantizers) Let X P L 2 R pPq and N P N ‹ .

(a) The quadratic distortion function Q 2,N at level N attains a minimum at a N -tuple x ‹ " px N 1 , . . . , N N q and Γ ‹ N " x N i , i P t1, . . . , N u ( is a quadratic optimal quantizer at level N .

(b) If the support of the distribution P X of X has at least N elements, then x ‹ " px N 1 , . . . , x N N q has pairwise distinct components, P X `Ci pΓ ‹ N q ˘ą 0, i P t1, . . . , N u. Furthermore, the sequence N Þ Ñ inf xPpRq N Q 2,N pxq converges to 0 and is decreasing as long as it is positive.

A really interesting and useful property concerning quadratic optimal quantizers is the stationary property, this property is deeply connected to the addressed problem after for the optimization of the quadratic optimal quantizers in (D.6).

Proposition D.6. (Stationarity) Assume that the support of P X has at least N elements. Any L 2 -optimal N -quantizer Γ N P pRq N is stationary in the following sense: for every Voronoï quantization p

Moreover P `X P Ť i"1,...,N BC i pΓ N q ˘" 0, so all optimal quantization induced by Γ N a.s. coincide.

The uniqueness of an optimal N -quantizer, due to Kieffer [START_REF] John | Exponential rate of convergence for lloyd's method i[END_REF], was shown in dimension one under some assumptions on the density of X.

Theorem D.7. (Uniqueness of optimal N -quantizers see [START_REF] John | Exponential rate of convergence for lloyd's method i[END_REF]) If P X pdξq " ϕpξqdξ with log ϕ concave, then for every N ě 1, there is exactly one stationary N -quantizer (up to the permutations of the N -tuple). This unique stationary quantizer is a global (local) minimum of the distortion function, i.e. @N ě 1, arg min

In what follows, we will drop the star notation (‹) when speaking of optimal quantizers, x ‹ and Γ ‹ N will be replaced by x and Γ N . The next result elucidates the asymptotic behavior of the distortion. We saw in Theorem D.5 that the infimum of the quadratic distortion converges to 0 as N goes to infinity. The next theorem, known as Zador's Theorem, establishes the sharp rate of convergence of the L p -mean quantization error.

Theorem D.8. (Zador's Theorem) Let p P p0, `8q.


M Usman

M Farooq

K Hanna

Existence of SARS-CoV-2 in wastewater: Implications for its environmental 1 transmission in developing communities 2

Australia, 2 Italy, [START_REF] Rosa | First detection of SARS-CoV-2 in 103 untreated wastewaters in Italy[END_REF] USA, 4 and France. [START_REF] Wurtzer | Time course 109 quantitative detection of SARS-CoV-2 in Parisian wastewaters correlates with 110 COVID-19 confirmed cases[END_REF] The first report of SARS-CoV-2 detection in 32 sewage came from the Netherlands suggesting its technical feasibility. 1 A study in the 33 Paris area of France reported the existence of SARS-CoV-2 in all raw wastewater 34 samples, and in 6 out of 8 treated wastewater samples. [START_REF] Wurtzer | Time course 109 quantitative detection of SARS-CoV-2 in Parisian wastewaters correlates with 110 COVID-19 confirmed cases[END_REF] Similarly, 6 out of 12 untreated 35 wastewater samples in Italy were positive for SARS-CoV-2. [START_REF] Rosa | First detection of SARS-CoV-2 in 103 untreated wastewaters in Italy[END_REF] A study in the USA 36 reported higher levels of SARS-CoV-2 in wastewater than would have been expected 37 from clinically-confirmed cases. 4 They estimated that roughly 5% of all faecal samples 38 were positive for SARS-CoV-2 in the testing period whereas the number of clinically-39 confirmed cases was 0.026%. In Australia, the number of infected individuals, 40 estimated from untreated wastewater, were in reasonable agreement with clinical 41 observations. 2 Analysing wastewater for SARS-CoV-2 has therefore been proposed 42 as a warning tool and complimentary approach to track the prevalence of infection. aerosols, created by a defective wastewater plumbing system, were identified as a 63 potential transmission route within a housing block in Hong Kong. [START_REF] Hung | The SARS epidemic in Hong Kong: what lessons have we 115 learned?[END_REF] Finally, the use of 64 wastewater as irrigation water has the potential to impact on soil quality and potential 65 uptake in crops or the contamination of groundwater resources by the virus if present.

66

In the developing world, an additional transmission route may be via faecal-oral 67 contact. Although hypothetical at this stage, the faecal-oral transmission has been 68 reported for many viral diseases such as Ebola Haemorrhagic fever, Hepatitis A,

69

Hepatitis E. [START_REF] Heller | COVID-19 faecal-oral transmission: Are 117 we asking the right questions?[END_REF] It becomes a crucial consideration when over half of the global population

  of coronavirus disease 2019 (COVID-19) is a public health 21 emergency of international concern. The disease, caused by severe acute respiratory 22 syndrome coronavirus 2 (SARS-CoV-2), emerged in China in December 2019 and was 23 declared a pandemic by the World Health Organization (WHO) on March 11 th , 2020 24 due to its world-wide spread. This disease has now been reported in over 213 25 counties/territories, with more than 5.5 million confirmed cases and over 0.35 million 26 deaths globally.

27SARS-

  CoV-2 is transmitted person-to-person via the inhalation of aerosol/droplets or 28 through contact with contaminated surfaces. Since the detection of SARS-CoV-2 in 29 faecal samples, it has also become clear that coronavirus is also present in human 30 A c c e p t e d m a n u s c r i p t

1- 5 43

 5 According to Reuters news (April 16, 2020), Australia's Government is already 44 planning to monitor its sewage for the presence of coronavirus . Wastewater 45 surveillance for SARS-CoV-2 could provide an unbiased opportunity to track its 46 epidemiology in countries having limited resources for clinical diagnosis.47 The presence of SARS-CoV-2 in wastewater may also have consequences for public 48 health in developing countries with poor water and sewage infrastructure, inadequate 49 institutional and technical disinfection capabilities, and lack of financing. According to 50 2017's World Water Development Report by United Nations, 80% of wastewater 51 worldwide (>95% in some developing countries) is released to the environment without 52 adequate treatment. Wastewater-related exposure to SARS-CoV-2 remains a 53 significant possibility in such vulnerable communities. For example, it is a common 54 practice to dump wastewater without adequate treatment into surface water bodies 55 (e.g. rivers, canals etc.) which downstream are being used for drinking water or 56 irrigation. Cases of accidental contamination of drinking water with raw sewage is 57 commonplace in developing countries, and even reported in the developed world. 6 58 A c c e p t e d m a n u s c r i p t wastewater could potentially aerosolize the virus if present. Although there is currently 60 no evidence of SARS-CoV-2 transmission via exposure to aerosolized wastewater, 61 this transmission route was identified during the SARS outbreak in 2003. Wastewater

  62

  70

( 4 . 9 71

 49 2 billion) lacks access to safely-managed sanitation. This highlights the need to consider the potential risks to developing communities from 72 environmental transmissions differently than those who are served by adequate 73 disinfection. We, therefore, call governments to integrate safety of wastewater 74 management, drinking water supply and recreational water environments in their fight 75 against the virus. Investments into wastewater infrastructure would be essential to 76 contribute in controlling outbreaks of waterborne diseases. Effective wastewater 77 disinfection to treat and minimize onward environmental transmission of infected 78 wastewater is also essential. A final disinfection step should be imposed urgently if 79 existing wastewater treatment plants are not optimized to remove viruses. There is 80 also an urgent need to provide clear information about the efficiency of the current 81 water disinfection treatments to stop the viral spread. The health of sanitation workers 82 should also be protected by following best safety practices.83In summary, a significant knowledge gap exists regarding transmission, persistence, 84 and fate of SARS-CoV-2 in wastewater and the environment requiring greater

  s c r i p t

surfaces and its persistence, mobility and fate in soils and the aquatic environment,

and potentially the food chain. 89 90