Sylvain Corlay

Joachim Lebovits

Jacques Lévy Véhel

Multifractional stochastic volatility models

Keywords: Hull & White model, functional quantization, vector quantization, Karhunen-Loève, Gaussian process, fractional Brownian motion, multifractional Brownian motion, white noise theory, S-transform, Wick-Itô integral, stochastic differential equations

.

Introduction

Volatility in financial markets is both of crucial importance and hard to model in an accurate way. It has been long known that a constant volatility as in the Black & Scholes model (see [START_REF] Black | The pricing of options and corporate liabilities[END_REF][START_REF] Merton | Theory of rational option pricing[END_REF]) is not consistent with empirical findings, such as the smile effect (i.e. the fact that volatility depends on both strike and maturity of an option). More basically, there is no reason to expect that instantaneous volatility should be constant. Since the late 80's, several models allowing for a varying volatility have appeared. The most popular ones include ARCH models and their generalizations (see [START_REF] Engle | Co-integration and error correction: Representation, estimation, and testing[END_REF] as well as [36, chap. 20 & 21] and references therein) in discrete time, and stochastic volatility models [START_REF] Hull | Pricing interest-rate derivative securities[END_REF][START_REF] Heston | A closed-form solution for options with stochastic volatility, with application to bond and currency options[END_REF] and the local volatility model [START_REF] Dupire | Pricing with a smile[END_REF] in continuous time. The local volatility model, in particular, is the only Markov diffusion process allowing one to exactly calibrate the marginals of the risk-neutral probability and thus to to reproduce observed implied volatility smiles. However, this model does not take into account another well documented fact: while stocks do not typically exhibit correlations, volatility does display long-range correlations (see, e.g. [START_REF] Andersen | The distribution of realized stock return volatility[END_REF]). Stochastic volatility models, in contrast, are able to incorporate this feature, provided an adequate driving noise is used. In [START_REF] Comte | Long memory in continuous-time stochastic volatility models[END_REF][START_REF] Comte | Affine fractional stochastic volatility models[END_REF], this is achieved through fractional integration. More precisely, the model considered in [START_REF] Comte | Long memory in continuous-time stochastic volatility models[END_REF] for the dynamic of the price of a risky asset reads as follows:

ß dS t = µ(t, S t )dt + S t σ t dW t , d ln(σ t ) = θ (µln(σ t )) dt + γdB H t , σ 0 > 0,

where W is a Brownian motion and B H t is an independent fractional Brownian motion (fBm) with parameter H under the historical probability. Section 2 recalls some basic facts about fBm. For now, it suffices to remind that its increments display long range dependence when H > 1/2, and that its pointwise Hölder regularity is almost surely equal to H at all times. Except when H = 1/2, fBm is not a semi-martingale, and thus neither will be ln(σ t ). However, this does not raise any problem: while price processes must be semimartingales due to absence of arbitrage constraints (see, e.g. [START_REF] Delbaen | The mathematics of arbitrage[END_REF]), such a requirement need not be imposed on instantaneous volatility. As explained in [START_REF] Comte | Long memory in continuous-time stochastic volatility models[END_REF], prices arising from Model (1) are indeed semi-martingales.

A prominent feature of this model is that it is consistent with the slow decay in the correlations of volatility observed in practice. It also accounts for two features related to the measured smile effect: the volatility process is less persistent in the short term than a standard diffusion, while it is more persistent in the long run ([19, p. 3]). We verify this fact in Section 6 by solving (1) numerically using recently developed functional quantization-based cubature methods.

By the very nature of this model, the evolution in time of the smile is governed by the single parameter H. This does not permit enough flexibility to fit volatility surfaces. In addition, (1) implies that the volatility has constant regularity equal to H. We provide in Section 1 empirical evidence that this is typically not the case by analysing records of S&P 500 and VIX volatility index. In addition, local regularity estimated on these data are often smaller than 1/2, thus ruling out the desirable feature of long range dependence if one insists on using fBm as a model. Although rather succinct, the numerical experiments of Section 1 indicate that a stochastic process with varying local regularity would provide a better fit to volatility.

The main idea of this work is to replace the fBm appearing in (1) with a more general process called multifractional Brownian motion (mBm). This is an extension of fBm where the parameter H is replaced by a smooth function h. By doing so, one obtains a process that has, at each time t, pointwise Hölder regularity equal to h(t), and, no matter the value of h in (0, 1), always display increments long-range dependence as long as h is not constant. Thus, replacing fBm by mBm in (1) enables the possibility of fitting the nonstationary local regularity of volatility as measured on data, as well as of maintaining long-range dependence properties. In addition, as we will show from numerical experiments in Section 6, the model (written in a risk-neutral setting):

   dF t = F t σ t dW t , d ln(σ t ) = θ (µln(σ t )) dt + γ h d ⋄ B h t + γ σ dW σ t , σ 0 > 0, d W, W σ t = ρdt, [START_REF] Alòs | Stochastic calculus with respect to Gaussian processes[END_REF] where B h t is an mBm, yields shapes of the smile at maturity T that are governed by a weighted average of the values of the function h up to time T : thus, by adequately choosing h, one may mimic a given implied volatility surface more faithfully than with the Hull & White model driven by fBm (1).

In order to give a rigorous meaning to the model above, a stochastic integral with respect to mBm must be defined. Multifractional and fractional Brownian motion are not semimartingales, thus classical Itô theory does not apply to them. At the time [START_REF] Comte | Long memory in continuous-time stochastic volatility models[END_REF] was written, no theory for integration with respect to fBm was available yet. Various approaches have been developed since, based mainly on Malliavin calculus [START_REF] Decreusefond | Stochastic analysis of the fractional Brownian motion[END_REF][START_REF] Alòs | Stochastic calculus with respect to Gaussian processes[END_REF], pathwise integrals [START_REF] Zähle | Integration with respect to fractal functions and stochastic calculus[END_REF], rough paths [START_REF] Coutin | An introduction to (stochastic) calculus with respect to fractional Brownian motion[END_REF], and white noise theory [START_REF] Elliott | A general fractional white noise theory and applications to finance[END_REF][START_REF] Biagini | An introduction to white-noise theory and Malliavin calculus for fractional Brownian motion[END_REF][START_REF] Bender | An S-transform approach to integration with respect to a fractional Brownian motion[END_REF]. As far as mBm is concerned, we note first that the Skohorod integral developed for instance in [START_REF] Alòs | Stochastic calculus with respect to Gaussian processes[END_REF] does not seem to be easily adaptable to mBm, as it would require writing B h as a Wiener integral over a finite interval. This appears to be an elusive task (see [START_REF] Herbin | Stochastic integration with respect to multifractional Brownian motion via tangent fractional Brownian motion[END_REF]Section 5] for details and [START_REF] Boufoussi | Local time and Tanaka formula for a Volterra-type multifractional Gaussian process[END_REF] for a related approach). As for the pathwise approach of [START_REF] Zähle | Integration with respect to fractal functions and stochastic calculus[END_REF], it extends immediately to the case of mBm, as it is only relies on the regularity properties of the couple integrator -integrand (see [START_REF] Zähle | Integration with respect to fractal functions and stochastic calculus[END_REF] and [START_REF] Herbin | Stochastic integration with respect to multifractional Brownian motion via tangent fractional Brownian motion[END_REF] for more details). However, solving stochastic differential equations (S.D.E.) typically requires H > 1/2 or h(t) > 1/2, which precludes its use in our context. The last approach, namely the white noise, lends itself naturally to an extension to mBm. In addition, this integral offers several advantages over both the pathwise and Skohorod ones already in the case of fBm: it allows to deal with any H ∈ (0, 1); contrarily to the pathwise integral, it is centred, and probabilistic quantities such as variances and expectations are easily computed; furthermore, the white noise integral is a generalization of the Skohorod one in the case of fBm, in the sense that if the latter one exists, then so does the former one, and both coincide [START_REF] Bender | An S-transform approach to integration with respect to a fractional Brownian motion[END_REF]. The white noise-based, or Wick-Itô, integral with respect to fBm was developed in [START_REF] Elliott | A general fractional white noise theory and applications to finance[END_REF][START_REF] Bender | An Itô formula for generalized functionals of a fractional Brownian motion with arbitrary Hurst parameter[END_REF][START_REF] Biagini | An introduction to white-noise theory and Malliavin calculus for fractional Brownian motion[END_REF], and applied to option pricing in a fractional Black and Scholes model in [START_REF] Elliott | A general fractional white noise theory and applications to finance[END_REF] (such a fractional Black and Scholes model raises both financial and economical issues, see [START_REF] Bender | Arbitrage with fractional Brownian motion? Theory Stoch[END_REF][START_REF] Björk | A note on Wick products and the fractional Black-Scholes model[END_REF]). It was extended to mBm in [START_REF] Lebovits | White noise-based stochastic calculus with respect to multifractional Brownian motion[END_REF]. This is the theory we will use in order to study precisely our stochastic volatility models.

The reader interested in the links between the Wick-Itô and various others integrals with respect to fBm may consult [START_REF] Bender | An S-transform approach to integration with respect to a fractional Brownian motion[END_REF]Theorem 6.2], [START_REF] Nualart | A white noise approach to fractional Brownian motion[END_REF]Proposition 8] or [START_REF] Coutin | An introduction to (stochastic) calculus with respect to fractional Brownian motion[END_REF]Corollary 8]. About stochastic differential equations driven by fBm in the Wick-Itô sense, see [START_REF] Holden | Stochastic partial differential equations[END_REF][START_REF] Kuo | White noise distribution theory[END_REF].

While we focus here on the multifractional stochastic volatility model (2) (we also briefly consider a multifractional SABR model with β = 1 in Section 5.3), we would like to mention that mBm is useful in a variety of applications in finance (see [1] for a partial list of articles dealing with mBm in this field).

In order to assess the relevance of our model, we compute numerically in Section 6 the smiles at different maturities. Since the solution cannot be written in an explicit form, we need to resort to approximations. In our case, this is made possible by recent advances in the theory of functional quantization of Gaussian processes.

Functional quantization of Gaussian processes has become an active field of research in recent years since the seminal article [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF]. As far as applications are concerned, cubature methods [START_REF] Pagès | Functional quantization for numerics with an application to option pricing[END_REF][START_REF] Corlay | Some aspects of optimal quantization and applications to finance[END_REF] and variance reduction methods [START_REF] Corlay | Functional quantization-based stratified sampling methods[END_REF][START_REF] Lejay | A variance reduction technique using a quantized Brownian motion as a control variate[END_REF] based on functional quantization have been proposed. However, as the numerical use of functional quantizers requires the evaluation of the Karhunen-Loève eigenfunctions, this method was restricted to processes for which a closed-form expression for this expansion is known, such as Brownian motion. In [START_REF] Corlay | The Nyström method for functional quantization with an application to the fractional Brownian motion[END_REF], a numerical method was proposed to perform numerical quadratic functional quantization of more general Gaussian processes, which will be applied here to multifractional Brownian motion.

We show that we can handle a fast and accurate forward start option pricing in this model thanks to a functional quantization-based cubature method similar to the one proposed in [START_REF] Pagès | Functional quantization for numerics with an application to option pricing[END_REF] and in [START_REF] Corlay | Some aspects of optimal quantization and applications to finance[END_REF]. This allows us to study the dependency of the smile dynamics on the functional parameter of the considered mBm.

The remaining of this paper is organized as follows. Section 1 proposes a short empirical study of the volatility of the S&P 500, which shows that its regularity evolves in time. We recall in Section 2 basic facts about fBm and mBm. In Section 3, we explain how to perform functional quantization of mBm and investigate the rate of decay of the corresponding quantization error. Quantization-based cubature is also addressed in this section. In Section 4, we provide some background on the white noise-based stochastic integral with respect to mBm. It also shows how to solve simple S.D.E. in this frame and presents general remarks on the quantization of solutions of S.D.E. A detailed treatment of the multifractional Hull & White and SABR models is proposed in Section 5. Numerical experiments, displaying the evolution of prices as a function of strike, as well as conclusions are gathered in Section 6.

A short empirical study of volatility

The aim of this section is to provide empirical evidence that volatility in financial markets is irregular, and that its local regularity, as measured by the pointwise Hölder exponent, evolves in time. It is by no means intended to present a complete statistical study of volatility, which is outside the scope of this work.

Local regularity and its estimation

Let us first recall the definition of the pointwise Hölder exponent of a real stochastic process X whose trajectories are continuous and nowhere differentiable. This is the stochastic process α X defined for every t as

α X (t) = sup ß β, lim sup h→0 |X(t + h) -X(t)| |h| β = 0 ™ .
When there is no risk of confusion, we shall write α in place of α X . If one wants to put this notion to use in practice, the first problem is of that of estimating α from discrete data. There is a wealth of works dealing with this issue. In this study, we shall use the so-called increment ratio statistic developed in [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional processes[END_REF]. The reader interested in other estimation methods in this area may consult the very complete list of reference in this article. Let us just mention here [START_REF] Peng | Inférence statistique pour des processus multifractionnaires cachés dans un cadre de modèles à volatilité stochastique[END_REF], which deals with statistical issues in the frame of a multifractional stochastic volatility model.

Experiments with S&P 500 data

The first problem we face is that we wish to estimate the local regularity of the volatility of stock prices. Volatility is not directly observed on the market. We have tested two classical strategies to obtain volatility signals. In the first one, studied in this subsection, historical volatility is estimated directly from high frequency records of prices. The second approach uses quoted vanilla option prices to obtain an estimation of the integrated local volatility, and is dealt with in the next subsection.

Our raw data in this section are minute quotes of the S&P 500, recorded from February 2, 2012 to July 23, 2012, which amounts to a total of 47748 samples. To estimate the historical volatility, we use essentially the same procedure as in [START_REF] Andersen | The distribution of realized stock return volatility[END_REF] (see also references therein). More precisely, we first compute the returns by taking logarithms of differences. We then process the data to remove the high frequency market microstructure noise using a low-pass filter. Samples are then grouped into blocks corresponding to a time period of four hours. The volatility attached to a block is then estimated as the standard deviation of the filtered samples contained in this block. See Figure 1 for a graph of the original S&P 500 series along with the associated estimated volatility. Note that our estimated historical volatility bears some resemblance with the ones displayed e.g. in [START_REF] Andersen | The distribution of realized stock return volatility[END_REF]. In both cases, the volatility appears to be highly irregular. We then estimate the local regularity of the volatility using the increment ratio statistic of [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional processes[END_REF], to obtain the results displayed on Figure 2. As we shall use this regularity as an input for our model below in Section 6, we need an analytical expression for it. We have thus regressed the raw regularity using a simple sine function, also shown on Figure 2. We denote this regressed function ◊ h V olSP . This is the function that will be used in our numerical experiments below.

Local regularity estimated in this way on the volatility of the S&P 500 is clearly not constant in time. It seems to oscillate with a period of roughly six weeks, and ranges approximately between 0.2 and 0.8.

Experiments with VIX data

Instead of estimating the volatility from raw data as above, another possibility is to use the VIX index, which is a popular measure of the implied volatility of the S&P 500 index options. Let us recall some basic facts about this index. Consider a filtration F = (F t ) t∈[0,T ] . We assume that the forward price of maturity T follows the dynamics dF t = σ t dW t on [0, T ], where W is a standard Brownian motion adapted to F under the riskneutral measure and where (σ t ) t∈[0,T ] is a stochastic process adapted to F . Applying Itô's lemma, one gets

E [ln(F T /F 0 )] = - 1 2 T 0 E[σ 2 s ]ds.
More generally, at any time t, the conditional expectation E î t+T t σ 2 s ds F t ó can be deduced from the logarithmic profile, which can be represented as a combination of Call and Put payoffs of the same maturity

ln(F/F 0 ) = F -F 0 F 0 + F0 0 (K -F ) + K 2 dK - ∞ F0 (F -K) + K 2 dK.
In practice, one can approximate this logarithm profile by a discrete combination of available Call and Put payoff. This is what was done by the Chicago Board Options Exchange to design the VIX volatility index, which is an approximation of the one-month variance swap rate of the S&P 500 index

t+τ t E[σ 2 s |F t ]
ds, where τ is equal to one month (see [START_REF]The CBOE volatility[END_REF] for the precise definition of the VIX index). A convex combination of option prices of the preceding and following listed maturities are used to approximate the desired one-month maturity option prices.

If one assumes that the VIX index is actually equal to the quantity

t+τ t E[σ 2 s |F t ]
ds where σ 2 s is the instantaneous variance of the S&P 500, then it seems intuitively plausible that its pointwise regularity at time t is 1 + min(α vol (t), α vol (t+ τ )), since integration typically amounts to adding one to the local regularity (a theoretical justification of this fact is still an open problem). Thus, if estimating α V IX (t) on VIX data yields a curve ranging in (1, 2) and varying in time, then such an empirical result would again support the use of mBm rather than fBm to model volatility.

Figure 3 displays our data, the VIX minute quotes from February 2, 2012 to June 20, 2012 (392865 data points). As above, we remove the market microstructure noise by low-pass filtering these data. In order to allow for a fair comparison with the previous experiment, we then subsample the series in order to obtain a four hours time period record. Finally, we estimate the pointwise regularity using again the increment ratio statistic, also shown on Figure 3.

As in the previous subsection, the estimated pointwise regularity for the volatility seems to clearly vary in time. It ranges between 0.2 and 0.55. Although an oscillatory behaviour is less pronounced than in the previous experiment, it is still apparent. The graphs in figure 2 and 3 indicate that volatility is indeed irregular, with varying pointwise Hölder exponent. As a consequence, it cannot be adequately modelled by fBm, since this process has constant regularity, as we recall in the next section. In contrast, mBm is able to capture this and other properties of volatility, such as long range dependence, as we explain below.

Background on multifractional Brownian motion

Fractional Brownian motion (fBm) [START_REF] Kolmogorov | Wienersche spiralen und einige andere interessante kurven im hilbertschen raume[END_REF][START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF] is a centred Gaussian process with features that make it a useful model in various applications such as financial and teletraffic modelling, image analysis and synthesis, geophysics and more. These features include self-similarity, long-range dependence and the ability to match any prescribed constant local regularity. Fractional Brownian motion depends on a parameter, usually denoted by H and called the Hurst exponent, that belongs to (0, 1). Its covariance function R H reads:

R H (t, s) := γ H 2 Ä |t| 2H + |s| 2H -|t -s| 2H ä ,
where γ H is a positive constant. A normalized fBm is one for which γ H = 1. Obviously, when H = 1 2 , fBm reduces to standard Brownian motion. While fBm is a useful model, the fact that most of its properties are governed by the single number H restricts its application in some situations. In particular, its Hölder exponent remains the same all along its trajectory. Thus, for instance, long-range dependent fBm, which require H > 1 2 , must have smoother paths than Brownian motion. Multifractional Brownian motion [START_REF] Peltier | Multifractional Brownian motion: definition and preliminary results[END_REF][START_REF] Benassi | Gaussian processes and pseudo-differential elliptic operators[END_REF] was introduced to overcome these limitations. The basic idea is to replace the real H by a function t → h(t) ranging in (0, 1).

The construction of mBm is best understood through the introduction of a fractional Brownian field. Fix a positive real T . A fractional Brownian field on [0, T ]×(0, 1) is a Gaussian field, denoted (B(t, H)) (t,H)∈[0,T ]×(0,1) , such that for every H in (0, 1) the process (B H t ) t∈[0,T ] , where B H t := B(t, H), is a fractional Brownian motion with Hurst parameter H.

For a deterministic continuous function h : [0, T ] → (0, 1), we call multifractional Brownian motion with functional parameter h the Gaussian process B h := (B h t ) t∈[0,T ] defined by B h t := B(t, h(t)). We say that h is the regularity function of the mBm. The fractional field (B(t, H)) (t,H)∈[0,T ]×(0,1) is termed normalized when, for all H in (0, 1), (B H t ) t∈[0,T ] is a normalized fBm. In this case we will also say that B h is normalized. In order for mBm to posses interesting properties, we need some regularity of B(t, H) with respect to H. More precisely, we will always assume that B(t, H) satisfies the following condition:

∀T ∈ R * + , ∀[c, d] ⊂ (0, 1), ∃(Λ, δ) ∈ (R * + ) 2 such that E[(B(t, H) -B(t, H ′ )) 2 ] ≤ Λ |H -H ′ | δ for every (t, H, H ′ ) in [0, T ] × [c, d] 2 . (H)
Under this assumption, and if the functional parameter h is continuous, then the associated mBm has a continuous modification. The class of mBm is rather large, since there is some freedom in choosing the correlations between the fBms composing the fractional field B. For definiteness, we will mostly consider in this work the so-called "wellbalanced" multifractional Brownian motion. Essentially the same analysis could be conducted with other versions. More precisely, a well-balanced mBm is obtained from the field B(t, H) := 1 cH R e itu -1 |u| H+1/2 W (du) where W denotes a complex-valued Gaussian measure (cf. [START_REF] Stoev | How rich is the class of multifractional Brownian motions? Stochastic Process[END_REF] for more details). We show in Proposition 3.5 that assumption (H) is satisfied by the well-balanced fractional Brownian field (in fact, it is verified by all mBms considered so far in the literature).

The proof of the following proposition can be found in [START_REF] Ayache | The covariance of the multifractional brownian motion and applications to long-range dependence[END_REF]:

Proposition 2.1 (Covariance function of well-balanced mBm). The covariance function R h of well-balanced mBm is given by

R h (t, s) = c 2 ht,s c h(t) c h(s) Å 1 2 |t| 2ht,s + |s| 2ht,s -|t -s| 2ht,s ã , (3)
where h t,s := h(t)+h(s)

2 and c x := Å 2π Γ(2x+1) sin(πx) ã 1 2 .
The other main properties of mBm are as follows: the pointwise Hölder exponent at any point t of B (h) is almost surely equal to h(t) ∧ β h (t), where β h (t) is the pointwise Hölder exponent of h at t [START_REF] Herbin | Stochastic 2-microlocal analysis[END_REF]. For a smooth h, one thus may control the local regularity of the paths by the value of h. In addition, the increments of mBm display long range dependence for all non-constant h(t) [START_REF] Ayache | The covariance of the multifractional brownian motion and applications to long-range dependence[END_REF]. Finally, when h is C 1 , mBm is tangent to fBm with exponent h(u) in the neighbourhood of any u in the following sense [START_REF] Falconer | Multifractional, multistable, and other processes with prescribed local form[END_REF]:

® B h u+rt -B h u r h(u) ; t ∈ [a, b] ´law ----→ r→0 + {B h(u) t ; t ∈ [a, b]}.
These properties show that mBm is a more versatile model that fBm: in particular, it is able to mimic in a more faithful way local properties of financial records, Internet traffic and natural landscapes [START_REF] Bianchi | Pathwise identification of the memory function of multifractional Brownian motion with application to finance[END_REF][START_REF] Li | Towards describing multi-fractality of traffic using local Hurst function[END_REF][START_REF] Echelard | Terrain modeling with multifractional brownian motion and self-regulating processes[END_REF] by matching their local regularity. This is important e.g. for purposes of detection or real-time control. The price to pay is of course that one has to deal with the added complexity brought by having a functional parameter instead of a single number. In general, the increments of multifractional Brownian motion are neither independent nor stationary. Since an mBm B h is an fBm of Hurst index H when h is constant and equal to H, there is no risk of confusion by denoting B H the fractional Brownian motion with Hurst index H.

We end this paragraph by noting that the use of a fractional Brownian field permits further generalizations: for instance, a multifractional process with random exponent is defined in [START_REF] Ayache | Multifractional processes with random exponent[END_REF], and a self-regulating process is considered in [START_REF] Echelard | Terrain modeling with multifractional brownian motion and self-regulating processes[END_REF].

Functional quantization of multifractional Brownian motion

In this section, we first present some background on functional quantization and the numerical methods used to obtain the quadratic optimal product quantization of mBm. We provide some numerical results in the special case of the well-balanced multifractional Brownian motion for certain examples of the functional parameter h. Then, in Section 3.2, we investigate the rate of decay of the quantization error for mBm. Section 3.3 presents the functional quantization-based cubature formulas that we use to devise a numerical scheme for the computation of option prices in the proposed multifractional stochastic volatility model.

Computation of the quantization

The quantization of a random variable X valued in a reflexive separable Banach space (E, | • |) consists in its approximation by a random variable Y that is measurable with respect to X and that takes finitely many values in E. The resulting error of the discretization is usually measured by the L p norm of the difference |X -Y |. If we settle on a fixed maximum cardinal N for Y (Ω), the minimization of the error reduces to the following optimization problem:

min ¶ |X -Y | p , Y : Ω → E measurable with respect to X, card(Y (Ω)) ≤ N © . ( 4
)
As Y is supposed to be measurable with respect to X, there exists a Borel map Proj :

E → E valued in a finite subset Γ of E such that Y = Proj(X). The finite subset Γ is called the codebook. Hence if Γ = {γ 1 , • • • , γ N }, there exists a Borel partition C = {C 1 , • • • , C N } of E such that Proj = N i=1 γ i 1 Ci .
In other words, Proj performs the process of mapping the continuous set X(Ω) to the finite set Γ. Let Proj Γ denote a nearest neighbour projection on Γ. Clearly,

|X -Proj Γ (X)| ≤ |X -Proj(X)| so that |X -Proj Γ (X)| p ≤ |X -Proj(X)| p .
Hence, in order to minimize the quantization error, it is optimal to use a nearest neighbour projection on the codebook Γ. A solution of ( 4) is called an L p -optimal quantizer of X. An elementary property of an L 2 -optimal quantizer is stationarity:

E[X|Y ] = Y . We denote by E N,p (X, | • |) the minimal L p quantization error for the random variable X and the norm | • |: E N,p (X, | • |) = min ¶ |X -Y | p , Y measurable with respect to X and |Y (Ω)| ≤ N © We now assume that X is a bi-measurable stochastic process on [0, T ] verifying T 0 E |X t | 2 dt
< ∞, that we see as a random variable valued in the Hilbert space H = L 2 ([0, T ]). Suppose that its covariance function Γ X is continuous. In [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF], it is shown that, in the centred Gaussian case, linear subspaces U of H spanned by N -stationary quantizers correspond to principal components of X, in other words, are spanned by eigenvectors of the covariance operator of X. Thus, the quadratic optimal quantization of Gaussian processes consists in using its Karhunen-Loève decomposition e X n , λ X n n≥1 . To perform optimal quantization, the Karhunen-Loève expansion is first truncated at a fixed order m and then the R m -valued Gaussian vector constituted of the m first coordinates of the process on its Karhunen-Loève decomposition is quantized. To reach optimal quantization, we have to determine both the optimal rank of truncation d X (N ) (the quantization dimension) and the optimal d X (N )-dimensional Gaussian quantizer corresponding to the first coordinates,

d X (N ) j=1 N 0, λ X j .
We have the following representation of the quadratic

distortion E N (X) := E N,2 X, | • | L 2 ([0,T ]) : E N (X) 2 = j≥m+1 λ X j + E N m j=1 N 0, λ X j 2 .
From a numerical viewpoint, we are thus confronted on the one hand with the finite-dimensional quantization of the Gaussian distribution

m j=1 N 0, λ X j
and on the other hand with the numerical evaluation of the first Karhunen-Loève eigenfunctions e X n 1≤n≤d X (N ) . Various numerical algorithms have been developed to deal with the first point. Let us mention Lloyd's algorithm and the Competitive Learning Vector Quantization (CLVQ). A review of these methods is available in [START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF]. As far as the evaluation of the first Karhunen-Loève eigenfunctions is concerned, closed-expressions are available for standard Brownian motion, standard Brownian bridge and Ornstein-Uhlenbeck process. Other examples of explicit Karhunen-Loève expansions may be found in [START_REF] Deheuvels | A Karhunen-Loeve decomposition of a Gaussian process generated by independent pairs of exponential random variables[END_REF] and [START_REF] Pycke | Explicit Karhunen-Loève expansions related to the Green function of the Laplacian[END_REF]. In the general case, the so-called Nyström method for approximating the solution of the associated integral equation may be used. It reads

T 0 Γ X (•, s)e X k (s)ds = λ X k e X k , k ≥ 1, (5)
where both the eigenvalues and the eigenvectors have to be determined. The Nyström method relies on the use of a quadrature scheme to approximate the integral, so that it turns into a matrix eigensystem.

When dealing with the midpoint quadrature rule, and for sufficiently regular kernels Γ X , the error admits an asymptotic expansion in the form of the sum of even powers of the step size, for both the eigenvalues and the eigenfunctions. We take advantage of this asymptotic expansion by using Richardson-Romberg extrapolation methods. This method has been benchmarked against the Karhunen-Loève eigensystems of standard Brownian motion, Brownian bridge and Ornstein-Uhlenbeck process in [START_REF] Corlay | The Nyström method for functional quantization with an application to the fractional Brownian motion[END_REF].

Instead of using an optimal quantization for the distribution

d X (N ) j=1 N 0, λ X j ,
another possibility is to use a product quantization, that is to use the Cartesian product of the optimal quadratic quantizers of the standard one-dimensional Gaussian distributions N 0, λ X j 1≤j≤d X (N ) . In the case of independent marginals, this yields a stationary quantizer, i.e. a quantizer Y of X which satisfies E[X|Y ] = Y . This property, shared with optimal quantizers, results in a convergence rate of a higher order for the quantization-based cubature scheme, as explained in [START_REF] Pagès | Functional quantization for numerics with an application to option pricing[END_REF]. An advantage of this approach is that one-dimensional Gaussian quantization is a fast procedure.

In [START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF], deterministic optimization methods (e.g. Newton-Raphson) are shown to converge rapidly to the unique optimal quantizer of the one-dimensional Gaussian distribution. A sharply optimized database of quantizers of standard univariate and multivariate Gaussian distributions is available for download on the web site [59], which is devoted to optimal quantization. One still has to determine the quantization level for each dimension to obtain optimal product quantization. In this case, the minimization of the distortion becomes:

Ä E prod N (X) ä 2 := min    d j=1 E 2 Nj N 0, λ X j + j≥d+1 λ X j , N 1 × • • • × N d ≤ N, d ≥ 1    . (6)
A solution of ( 6) is called an optimal K-L product quantizer. This problem can be solved by the "blind optimization procedure", which consists in computing the criterion for every possible decomposition

N 1 ו • •× N d with N 1 ≥ • • • ≥ N d .
The result of this procedure can be stored for future use. Optimal decompositions for a wide range of values of N for both Brownian motion and Brownian bridge are available on the web site [59]. Another fact on quadratic functional product quantization is that it is shown to be rate-optimal under certain assumptions on the K-L eigenvalues (see Theorem 3.1).

Quadratic product quantizers of fBms and well-balanced mBms for different H and h are displayed on Figures 4 and5. A fixed product decomposition is used for simplicity.

These graphs reflect, to some extent, the features of the quantized process, in particular its correlation and regularity properties:

In the case of fBm (Figure 4), when H increases, the rate of decay of the Karhunen-Loève eigenvalues also increases (and so does the pathwise Hölder regularity of the underlying process), so that even though we do not change the quantization dimension in this example, the contribution of higher-order Karhunen-Loève eigenvalues decreases. In Figure 4, one can see that the curves of the functional quantizer localize around the quantization of the first Karhunen-Loève coordinate when H = 0.75, while this is not the case when H = 0.25.

In addition, the distribution of the curves on the plane is related to the fact that the almost sure Hausdorf dimension of the paths of fBm is 2-H: for small H, we expect the set of curves to be more space-filling than for large H, a feature that can be indeed be verified on the figure. In addition, the long-term correlation of fBm for H > 1/2, which results in paths typically having strong trends, translates here into curves which are roughly monotonous. Conversely, the negative correlations which characterizes the case H < 1/2 is reflected in the more oscillatory behaviour of the curves in the left pane.

The case of mBm (Figure 5) makes even clearer the relation between the properties of the process and the optimal quantizer. In the right pane of Figure 5, the function h decreases linearly from 0.9 to 0.1. One can see that, for small t, both the distribution of the curves and their trend look like the ones of fBm with large H. As t increases, the curves become more space-filling and oscillatory, in agreement with the fact that, for t close to 1, the corresponding mBm has larger local Hausdorff dimension. Similar remarks hold for the case where h is an increasing function (right pane of Figure 5).

Another way of interpreting these figures is to recall that mBm is tangent, at each t, to fBm with exponent H = h(t). The behaviour of the curves on Figure 5 is a translation of this fact in the quantization domain.

Finally, note that the shape of the convex envelopes in each of the four figures roughly matches the time evolution of the variances of the corresponding processes, i.e. t 0.25 , t 0.75 , t 0.1+0.8t and t 0.9-0.8t .

Rate of decay of the quantization error for mBm

The rate of decay of the quadratic functional quantization error was first investigated in [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF]. More precise results were then established in [START_REF] Luschgy | Sharp asymptotics of the functional quantization problem for Gaussian processes[END_REF]. These results rely on assumptions on the asymptotic behaviour of the Karhunen-Loève eigenvalues of the considered process. In Subsection 3.2.1, we recall the main result involving the rate of decay of these eigenvalues, leading to sharp rates of convergence for the quantization of fBm.

Unfortunately, such asymptotics for the Karhunen-Loève eigenvalues are not known at this time in the case of mBm. However, since the regularity of mBm is known, one may use another, less precise, type of results: these yield an upper estimate on the rate of decay of the quantization error [START_REF] Luschgy | Functional quantization rate and mean regularity of processes with an application to Lévy processes[END_REF]. This is explained in Subsection 3.2.2.

In the following, for two positive sequences (x n ) n∈N and (y n ) n∈N , we write

x n ∼ n→∞ y n if lim n→∞ xn yn = 1. The symbol x n n→∞ y n means that lim n→∞ xn yn ≤ 1. Finally, x n ≍ n→∞ y n means that x n = O(y n ) and y n = O(x n ) as n → ∞.

Sharp rates based on asymptotics of Karhunen-Loève eigenvalues

Recall the following well-known definition:

Definition 1 (Regularly varying function at infinity). A measurable function φ : (s, ∞) → (0, ∞), (s > 0) is regularly varying at infinity with index b ∈ R if for every t > 0, lim x→∞ φ(tx) φ(x) = t b .
Let X be a bi-measurable centred Gaussian process on [0, T ] with a continuous covariance function Γ X and such that

T 0 E[X 2 s ]ds < ∞.
Denote by e X n , λ X n n≥1 its Karhunen-Loève eigensystem.

Theorem 3.1 (Quadratic quantization error asymptotics [START_REF] Luschgy | Sharp asymptotics of the functional quantization problem for Gaussian processes[END_REF]). Assume that λ

X n ∼ φ(n) as n → ∞, where φ : (s, ∞) → (0, ∞) is a decreasing regularly varying function of index -b < -1 and s > 0. Set ψ(x) := 1 xφ(x) . Then E N (X) ∼ Ç Å b 2 ã b-1 b b -1 å 1/2 ψ(log(N )) -1/2 as N → ∞.
Moreover, the optimal product quantization dimension m X (N ) verifies m X (N ) ∼ 2 b log(N ) as N → ∞, and the optimal product quantization error E prod N (X) of level N satisfies

E prod N (X) Ç Å b 2 ã b-1 b b -1 + C(1) å 1/2 ψ(log(N )) -1/2 as N → ∞,
where C(1) is a universal positive constant.

Though the optimal product quantization is not asymptotically optimal, it still provides a rate-optimal sequence of quantizers. In the case where b = 1, a similar result is true, with the additional property that the optimal product quantization does yield an asymptotically optimal quadratic quantization error.

The case of fractional Brownian motion

In [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF][START_REF] Bronski | Asymptotics of Karhunen-Loeve eigenvalues and tight constants for probability distributions of passive scalar transport[END_REF], it is shown that the Karhunen-Loève eigenvalues of fBm on [0, T ] verify

λ B H n ∼ ν H n 2H+1 as n → ∞,
where ν H is a positive constant. Thus, fBm satisfies the hypotheses of Theorem 3.1 and

E N B H ∼ K H log(N ) H as N → ∞ for some K H > 0, and E prod N B H ≍ 1 log(N ) H as N → ∞.

Mean regularity and domination of the functional quantization rate

We recall the definition of regular variation at 0: Definition 2 (Regularly varying function at zero). A measurable function φ : (0, s) → (0, ∞), (s > 0) is regularly varying at 0 with index b ≥ 0 if for every t > 0, lim

x→0 φ(tx) φ(x) = t b .
Definition 3 (The φ-Lipschitz assumption). Let X be a bi-measurable process on [0, T ]. We say that X satisfies the φ-Lipschitz assumption for ρ > 0, which we denote by (L φ,ρ ), if there is a non-decreasing function φ : R + → [0, ∞], continuous at 0 with φ(0) = 0, such that

(L φ,ρ ) ≡    ∀(s, t) ∈ [0, T ] 2 , E [|X t -X s | ρ ] ≤ (φ(|t -s|)) ρ , if ρ ≥ 1 ∀t ∈ [0, T ], ∀h ∈ (0, T ], E ñ sup t≤s≤(t+h)∧T |X s -X t | ρ ô ≤ (φ(h)) ρ if 0 < ρ < 1. Remark 1. The φ-Lipschitz assumption implies that E î |X| ρ L ρ ([0,T ]) ó < ∞ so that P-almost surely, t → X t lies in L ρ ([0, T ]).
Theorem 3.2 (Mean regularity and quantization rate). Let X be a bi-measurable process on [0, T ] such that X t ∈ L ρ for every t ∈ [0, T ], ρ > 0. Assume that X satisfies (L φ,ρ ) where φ is regularly varying at 0 with index b. Then

∀(r, p) ∈ [0, ρ] 2 , E N,r X, | • | L p ([0,T ]) ≤ C r,p ß φ(1/ log(N )), if b > 0, ψ(1/ log(N )), if b = 0, with ψ(x) := x 0 (φ(ξ)) (r∧1) ξ dξ 1 r∧1
, assuming in addition that

x 0 (φ(ξ)) (r∧1) ξ dξ < ∞ if b = 0. In particular, if φ(u) = cu b , b > 0, then E N,r X, | • | L p ([0,T ]) = O(log(N ) -b ).

The case of multifractional Brownian motion

Recall that a function h is said to be β-Hölder continuous (with β > 0) if there exists η in R * + such that, forall (s, t) in [0, T ] 2 , |h(s)h(t)| ≤ η |s -t| β . Theorem 3.3 (L 2 -mean regularity of multifractional Brownian motion). Let B h be an mBm with functional parameter h satisfying assumption (H). Assume that h is β-Hölder continuous, then there exists a positive constant M such that

∀(s, t) ∈ [0, T ] 2 , E B h t -B h s 2 ≤ M |t -s| 2 inf u∈[0,T ] h(u)∧βδ , ( 7
)
where δ is given in assumption (H).

Proof: We may assume that the fractional field (B(t, H)

) (t,H)∈[0,T ]×[c,d] is normalized. For (t, s) in [0, T ] 2 : E B h t -B h s 2 ≤ 2 E î (B(t, h(t)) -B(s, h(t))) 2 ó + 2 E î (B(s, h(t)) -B(s, h(s))) 2 ó ≤ 2 Ä |t -s| 2h(t) + Λ |h(t) -h(s)| δ ä ≤ 2 Ä |t -s| 2H1 (1 + T 2(H2-H1) ) + Λ η β |t -s| βδ ä ≤ 2 (1 + T 2(H2-H1) ) (1 + Λ η β ) |t -s| 2H1 + |t -s| βδ ≤ M |t -s| 2H1∧βδ ,
where

H 1 := inf u∈[0,T ] h(u), H 2 := sup u∈[0,T ] h(u) and M := 2(1 + T 2(H2-H1) ) (1 + T 2H1∨βδ-2H1∧βδ ) (1 + Λ η β ).
Corollary 3.4 (Upper bound on the quantization error for multifractional Brownian motion). With the same notations and assumptions as in Theorem 3.3:

E N,r B h , | • | L p ([0,T ]) = O Ä log(N ) -(H1∧ βδ 2 )
ä , for every (r, p) in (R + ) 2 .

Proof: Since B h is a Gaussian process, Theorem 3.3 shows that B h fulfils the φ-Lipschitz assumption for every integer ρ of the form ρ := 2p where p is a positive integer and for the continuous function φ ρ defined on R + by φ ρ (0) := 0 and φ ρ (x) := (κ ρ/2 ) 1/2ρ √ M x H1∧ βδ 2 . We have denoted, for n in N, κ n the number such that E Y 2n = κ n E Y 2 n for the centred Gaussian random variable Y . It is clear that φ ρ is regularly varying with index H 1 ∧ βδ 2 , which is positive. The result then follows from Theorem 3.2. Remark 2. Corollary 3.4 extends to every process

V h := (V h t ) t∈[0,T ] , of the form V h t := Z(t, h(t)) where Z := (Z(t, H)) (t,H)∈R×[H1,H2] is a Gaussian field such that one can find (Λ, γ, δ) in (R * + ) 3 with ∀(s, t, H, H ′ ) ∈ [0, T ] 2 × [H 1 , H 2 ] 2 , E î Z(t, H) -Z(s, H ′ ) 2 ó ≤ Λ (|t -s| γ + |H -H ′ | δ ),
In this case, for every (r, p)

in (R * + ) 2 , we get E N,r V h , | • | L p ([0,T ]) = O log(N ) -( γ 2 ∧ βδ 2 ) .
Proposition 3.5. The well-balanced mBm

B h t := 1 c h(t) R e itξ -1 |u| h(t)+1/2 W (dξ)
satisfies assumption (H).

Proof: One computes:

I H,H ′ t := E î (B(t, H) -B(t, H ′ )) 2 ó = R e itξ -1 c H |ξ| H+1/2 - e itξ -1 c H ′ |ξ| H ′ +1/2 2 du = R e itξ -1 ξ 2 1 cH |ξ| 1/2-H -1 c H ′ |ξ| 1/2-H ′ 2 dξ. (8)
For every ξ in R * , the map

f ξ : [c, d] → R + , defined by f ξ (H) := 1 cH |ξ| 1/2-H is C 1 since H → c H is C 1 on (0, 1).
Thus there exists a positive real D such that

∀(ξ, H) ∈ R * × [c, d], |f ′ ξ (H)| ≤ D |ξ| 1/2-H (1 + | ln(|ξ|)|) ≤ D Ä |ξ| 1/2-c + |ξ| 1/2-d ä (1 + | ln(|ξ|)|).
Thanks to the mean-value theorem, (8) yields

I H,H ′ t ≤ D 2 |H -H ′ | 2 R |e itξ -1| 2 |ξ| 2 Ä |ξ| 1/2-c + |ξ| 1/2-d ä 2 (1 + | ln(|ξ|)|) 2 dξ ≤ D 2 |H -H ′ | 2 Ç 2 3 |ξ|>1 (1 + | ln(|ξ|)|) 2 |ξ| 1+2c dξ + (2t) 2 |ξ|≤1 |ξ| 1-2d (1 + | ln(|ξ|)|) 2 dξ å ≤ (2 3 + T 2 ) D 2 Ç |ξ|>1 (1 + | ln(|ξ|)|) 2 |ξ| 1+2c dξ + |ξ|≤1 |ξ| 1-2d (1 + | ln(|ξ|)|) 2 dξ å |H -H ′ | 2 .
Since the two integrals in the last line are finite, (H) is verified with δ = 2.

Remark 3 (Quantization error and small ball probability). In the case of Gaussian measures, upper bounds and lower bounds of the quantization error can be related to lower and upper bounds for small ball probabilities respectively [START_REF] Dereich | On the link between small ball probabilities and the quantization problem for Gaussian measures on Banach spaces[END_REF], and a converse relationship was obtained in [START_REF] Graf | Functional quantization and small ball probabilities for Gaussian processes[END_REF]. As a consequence, the knowledge of logarithmic small ball asymptotics gives asymptotics for the quantization error. Let us mention that tight asymptotics of the L ∞ small ball probability were obtained for a special case of multifractional Brownian motion in [START_REF] Ayache | The small ball behavior of a non-stationary increments process: the multifractional Brownian motion[END_REF].

Quantization-based cubature

In this section, we first recall the error bounds on quantization-based cubature formulas. We then handle the case of exponentials of continuous centred Gaussian processes.

Basic formula and related inequalities in the case of Lipschitz continuous functionals

The idea of quantization-based cubature methods is to approximate the distribution of the random variable X by the distribution of a quantizer Y of X. As Y is a discrete random variable, we have

P Y = N i=1 p i δ yi . Therefore, if F : E → R is a Borel functional, E[F (Y )] = N i=1 p i F (y i ). (9)
Hence, the weighted discrete distribution (y i , p i ) 1≤i≤N of Y allows us to compute the sum [START_REF] Bender | An Itô formula for generalized functionals of a fractional Brownian motion with arbitrary Hurst parameter[END_REF]. We review some error bounds that can be derived when approximating E[F (X))] by [START_REF] Bender | An Itô formula for generalized functionals of a fractional Brownian motion with arbitrary Hurst parameter[END_REF]. See [START_REF] Pagès | Functional quantization for numerics with an application to option pricing[END_REF] for more details.

1.

If X ∈ L 2 , Y a quantizer of X of size N and F is Lipschitz continuous, then |E[F (X)] -E[F (Y )]| ≤ [F ] Lip X -Y 2 . ( 10
)
where

[F ] Lip is the Lipschitz constant of F . In particular, if (Y N ) N ≥1 is a sequence of quantizers such that lim N →∞ X -Y N 2 = 0, then the distribution N i=1 p N i δ x N i
of Y N weakly converges to the distribution

P X of X as N → ∞. 2. If Y is a stationary quantizer of X, i.e. Y = E[X|Y ],
and F is differentiable with a Lipschitz continuous derivative DF , then

|E[F (X)] -E[F (Y )]| ≤ [DF ] Lip X -Y 2 2 , (11)
where [DF ] Lip is the Lipschitz constant of DF . If F is twice differentiable and D 2 F is bounded, then we can replace [DF ] Lip by 1 2 D 2 F ∞ .

3.

If F is a semi-continuous convex functional and Y is a stationary quantizer of X,

E[F (Y )] ≤ E[F (X)]. (12)
This is a simple consequence of Jensen's inequality. Indeed,

E[F (Y )] Stationarity = E[F (E[X|Y ])] Jensen ≤ E[E[F (X)|Y ])] = E[F (X)].
Remark 4. In the infinite-dimensional case, convexity does not imply continuity. In infinite-dimensional Banach spaces, a semi-continuity hypothesis is necessary for Jensen's inequality which is the reason why we had to make this additional hypothesis on F . See [START_REF] Zapała | Jensen's inequality for conditional expectations in Banach spaces[END_REF] for more details.

The case of exponentials of continuous centred Gaussian processes

Let (X s ) s∈[0,T ] be a continuous centred Gaussian process on [0, T ]. Then the covariance function of X is also continuous. In addition, Fernique's theorem entails that E î T 0 X 2 s ds ó is finite. We view X as a random variable valued in the separable Banach space (C([0, T ], R), • ∞ ). Let " X be a stationary quantizer of X.

By the mean-value theorem, for all (x, y) ∈ R 2 , |e xe y | ≤ e |x|+|y| |x -y|. Consequently, for p ≥ 1, using Hölder's inequality:

E T 0 e Xs -e Xs p ds 1/p ≤ E ï T 0 e p|Xs|+p Xs X s -" X s p ds ò 1/p ≤ E ï T 0 e p p|Xs| e p p Xs ds ò 1 2p p E ï T 0 X s -" X s pq ds ò 1 2p q ,
where (p, q) ∈ (1, ∞) 2 are conjugate exponents. For ǫ > 0, we choose (p, q) such that pq = p + ǫ. This gives q = 1 + ǫ/p and p = 1 + p/ǫ.

By

We shall apply [START_REF] Bianchi | Pathwise identification of the memory function of multifractional Brownian motion with application to finance[END_REF] with p = 2ǫ in Section 5.2: this will allow us to control the L 2-ǫ quantization error of the exponential of a centred continuous Gaussian process X by the L 2 quantization of X.

Richardson-Romberg extrapolation

The accuracy of quantization-based cubature formulas can be dramatically improved by the use of Richardson-Romberg extrapoltion methods, with respect to the quantization error, of with respect to the quantization level. Here is a brief presentation of these methods.

With respect to the quantization error

In the general setting of a non-uniform random variable X, a quadratic optimal N -quantizer Y N of X and a C 1 functional with Lipschitz continuous derivative, Equation [START_REF] Bender | Arbitrage with fractional Brownian motion? Theory Stoch[END_REF] does not provide a true asymptotic expansion which would allow one to use a Richardson-Romberg expansion, but it suggests the use a higher-order Taylor expansion of F (X) -F (Y N ) to get one.

It follows from Taylor's formula that there exists a vector

ζ ∈ [X, Y N ] such that E[F (X)] = E [F (Y N )] + E [ DF (Y N ), X -Y N ] =E[DF (YN ).E[X-YN |YN ]]=0 by stationarity. + 1 2 E D 2 F (Y N )(X -Y N ) ⊗2 + 1 6 E ζ(X -Y N ) ⊗3 + o E |X -Y N | 3 = E [F (Y N )] + 1 2 E D 2 F (Y N )(X -Y N ) ⊗2 + O Ä E î |X -Y N | 3 óä . (14)
In [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF], it is proved that the asymptotics of the L s quantization error induced by a sequence of L r -optimal quantizers remains rate-optimal in the case of probability distributions on R d , with s < r + d for a class of distributions including the Gaussian distribution. This leads to E |X -

Y N | 3 = O E |X -Y N | 2 3 2
. This holds e.g. for Brownian motion.

Unfortunately, no sharp equivalence between X -Y N

E[F (Y l )]E 2 k -E[F (Y k )]E 2 l E 2 k -E 2 l . (15)
Although this kind of Richardson-Romberg extrapolation has not received a full theoretical justification yet, it does dramatically increase the efficiency of quantization-based cubature formulas.

With respect to the quantization level

When the value of E 2 k is not known, one may rely on an asymptotic expansion with respect to the quantization level.

Remark 5 (Romberg extrapolation with respect to the quantization level). In Section 3.2, we have seen that under some assumptions on the eigenvalues of the convergence operator, the rate of convergence of optimal quantizers and K-L optimal product quantizers is (ln(N ) -α ) for some α ∈ (0, 1). Replacing the distortion E N by its asymptotics 1 ln(N ) α as N → ∞ in Equation ( 15) suggests to approximate E[F (X)] by the quantity

E[F (Y l )](ln l) 2α -E[F (Y k )](ln k) 2α (ln l) 2α -(ln k) 2α . (16)

Multidimensional Richardson-Romberg extrapolation

Let X 1 and X 2 be two independent random variables. We wish to estimate the expectation E F X 1 , X 2 for some regular functional F . In that view, one may use a cubature based on a product quantization Ä " X 1 , " X 2 ä of (X 1 , X 2 ), and perform a multidimensional Richardson-Romberg extrapolation. This amounts to performing two Richardson-Romberg extrapolations as described already, one related to the quantization error of X 1 between quantization levels N 1 and M 1 , and one related to the quantization error of X 2 between quantization levels N 2 and M 2 . This leads to approximating E[F (X 1 , X 2 )] by the quantity

E 2 M1 E 2 M2 F N1,N2 -E 2 N1 E 2 M2 F M1,N2 -E 2 M1 E 2 N2 F N1,M2 + E 2 N1 E 2 N2 F M1,M2 E 2 M1 -E 2 N1 E 2 M2 -E 2 N2 , (17)
where F p,q denotes the estimated expectation obtained with the quantization-based cubature and quantization levels of p and q for X 1 and X 2 respectively. In other words, F p,q is defined by

F p,q := E F " X 1 p , " X 2 q
where " X 1 p , " X 2 q are quantizers of levels p and q for X 1 and X 2 respectively. In Equation ( 17), E Mi and E Ni denote the quadratic quantization error of level M i and N i for X i .

Stochastic calculus with respect to mBm

From now on and until the end of the work, we fix our mBm to be the well-balanced multifractional Brownian motion defined in Section 2. In addition, we will always assume that h is a C 1 function with derivative bounded on R.

Some backgrounds on white noise theory

This section provides the necessary background on white noise theory to define a stochastic integral and to handle S.D.E.s driven by mBm.

Define the probability space as Ω := S ′ (R) and let F := B(S ′ (R)) be the σ-algebra of Borel sets. There exists a probability measure µ such that, for every f in L 2 (R), the map •, f : Ω → R defined by

•, f (ω) = ω, f (where ω, f is by definition ω(f ), i.e the action of the distribution ω on the function f ) is a centred Gaussian random variable with variance equal to f 2 L 2 (R) under µ. For every n in N, denote e n (x) := (-1) n π -1/4 (2 n n!) -1/2 e x 2 /2 d n dx n (e - As is customary, we denote (L 2 ) the space L 2 (Ω, G, µ) where G is the σ-field generated by ( •, f ) f ∈L 2 (R) . For every random variable Φ of (L 2 ) there exists, according to the Wiener-Itô theorem, a unique sequence

(f n ) n∈N of functions f n in L 2 (R n ) such that Φ can be decomposed as Φ = +∞ n=0 I n (f n ), where L 2 (R n ) denotes
the set of all symmetric functions f in L 2 (R n ) and I n (f ) denotes the nth multiple Wiener-Itô integral of f with the convention that I 0 (f 0 ) = f 0 for constants f 0 . Moreover we have the equality

E[Φ 2 ] = +∞ n=0 n! f n 2 L 2 (R n )
where E denotes the expectation with respect to µ. For any Φ :=

+∞ n=0 I n (f n ) satisfying the condition +∞ n=0 n! |A ⊗n f n | 2 0 < +∞, define the element Γ(A)(Φ) of (L 2 ) by Γ(A)(Φ) := +∞ n=0 I n (A ⊗n f n )
, where A ⊗n denotes the nth tensor power of the operator A (see [START_REF] Janson | Gaussian Hilbert Spaces[END_REF]Appendix E] for more details about tensor products of operators). The operator Γ(A) is densely defined on (L 2 ). It is invertible and its inverse Γ(A)

-1 is bounded. Let us denote ϕ 2 0 := ϕ 2 (L 2 ) for ϕ in (L 2 ) and let Dom(Γ(A) n ) be the domain of the nth iteration of Γ(A). Define the family of norms ( • p ) p∈Z by:

Φ p := Γ(A) p Φ 0 = Γ(A) p Φ (L 2 ) , ∀p ∈ Z, ∀Φ ∈ (L 2 ) ∩ Dom(Γ(A) p ).
For any p in N, let (S p ) := {Φ ∈ (L 2 ) : Γ(A) p Φ exists and belongs to (L 2 )} and define (S -p ) as the completion of the space (L 2 ) with respect to the norm • -p . As in [START_REF] Kuo | White noise distribution theory[END_REF], we let (S) denote the projective limit of the sequence ((S p )) p∈N and (S)

* the inductive limit of the sequence ((S -p )) p∈N . The space (S) is called the space of stochastic test functions and (S) * the space of Hida distributions. One can show that, for any p in N, the dual space (S p ) * of S p is (S -p ). Thus we will write (S -p ), in the sequel, to denote the space (S p ) * . Note also that (S) * is the dual space of (S). We will note •, • the duality bracket between (S) * and (S). If Φ belongs to (L 2 ) then we have the equality Φ, ϕ = Φ, ϕ (L 2 ) = E[Φ ϕ]. Since we have defined a topology given by a family of norms on the space (S) * it is possible to define a derivative and an integral in (S) * . A function Φ : R → (S)

* is called a stochastic distribution process, or an (S) * -process, or a Hida process.

The Hida process Φ is said to be differentiable at

t 0 if lim r→0 r -1 (Φ(t 0 + r) -Φ(t 0 )) exists in (S) * .
Moreover we may also define an integral of an Hida process:

Theorem 4.1 (Integral in (S) * ). Assume that Φ : R → (S) * is weakly in L 1 (R, dt), i.e. assume that for all ϕ in (S), the mapping u → Φ(u), ϕ from R to R belongs to L 1 (R, dt). Then, there exists a unique element in (S) * , denoted by R Φ(u)du, such that

≠ ≠ R Φ(u)du, ϕ ∑ ∑ = R Φ(u), ϕ du for all ϕ in (S). ( 18
)
One says that Φ is (S) * -integrable on R or integrable on R in the Pettis sense.

For every f in L 2 (R), define the Wick exponential of •, f , noted : e •,f :, as the (L 2 ) random variable equal to e •,f -1 2 |f | 2 0 . The S-transform of an element Φ of (S * ), noted S(Φ), is defined to be the function from S (R) to R given by S(Φ)(η) := Φ, : e •,η for every η in S (R). Finally for every (Φ, Ψ) ∈ (S) * × (S) * , there exists a unique element of (S) * , called the Wick product of Φ and Ψ and noted Φ ⋄ Ψ, such that S(Φ ⋄ Ψ)(η) = S(Φ)(η) S(Ψ)(η); for every η in S (R).

The map S : Φ → S(Φ), from (S) * to (S) * , is injective. Furthermore, let Φ : R → (S) * be a fixed (S) * process. If Φ is (S) * -integrable over R then for all η in S (R),

S( R Φ(u) du)(η) = R S(Φ(u))(η) du. If Φ is (S) * -differentiable over R then for all η in S (R), S[ dΦ dt (t)](η) = d dt [SΦ(t)](η) .
For any Φ in (S) * and k in N * , let Φ ⋄k denote the element

k times Φ ⋄ • • • ⋄ Φ of (S) * .
One can generalize the definition of exp ⋄ to the case where Φ belongs to (S) * . Indeed, for any Φ in (S) * such that the sum

+∞ k=0 Φ ⋄k k! converges in (S) * , define the element exp ⋄ Φ of (S) * by setting exp ⋄ Φ := +∞ k=0 Φ ⋄k k! . It is called Wick exponential of Φ.
For f in L 2 (R) and Φ := •, f , it is easy to verify that exp ⋄ Φ exists and coincides with : e •,f : defined at the beginning of this section.

Fractional and multifractional white noise

Operators M H and ∂MH ∂H . Let H belong to (0, 1). Following [START_REF] Elliott | A general fractional white noise theory and applications to finance[END_REF], the operator M H is defined in the Fourier domain by

◊ M H (u)(y) := √ 2π cH |y| 1/2-H u(y), ∀y ∈ R * .
This operator is well defined on the homogeneous Sobolev space of order

1/2 -H noted L 2 H (R) and defined by L 2 H (R) := {u ∈ S ′ (R) : u = T f ; f ∈ L 1 loc (R) and u H < +∞}. The norm • H derives from the inner product •, • H defined on L 2 H (R) by: u, v H := 1 c 2 H R |ξ| 1-2H " u (ξ)" v (ξ)dξ where c H is defined right after Definition 2.1.
The definition of the operator ∂MH ∂H is quite similar [START_REF] Lebovits | White noise-based stochastic calculus with respect to multifractional Brownian motion[END_REF]. Precisely, define, for H in (0, 1), the space Γ H (R) := {u ∈ S ′ (R) : u = T f ; f ∈ L 1 loc (R) and u δH (R) < +∞}, where the norm

• δH (R) derives from the inner product •, • δH defined on Γ H (R) by u, v δH := 1 c 2 H R (β H + ln |ξ|) 2 |ξ| 1-2H " u (ξ) dξ. The operator ∂MH ∂H , from Γ H (R), •, • δH (R) to L 2 (R), •, • L 2 (R)
, is defined in the Fourier domain by:

ÿ ∂MH ∂H (u)(y) := -(β H + ln |y|) √ 2π cH |y| 1/2-H u(y), ∀y ∈ R * .

Fractional and multifractional white noise

For any measurable function h : R → (0, 1), it is easily seen that the process B h := B h t t∈R defined by

∀(ω, t) ∈ Ω × R, B h t := +∞ k=0 Ç t 0 M h(t) (e k )(s)ds å
•, e k is an mBm. Assuming that h is differentiable, we define the (S) * -valued function W h := (W h t ) t∈R by

W h t := +∞ k=0 ï d dt Å t 0 M h(t) (e k )(s) ds ãò •, e k . ( 19
)
The following theorem states that, for all real t, the right-hand side of ( 19) does indeed belong to (S) * and is exactly the (S) * -derivative of B h at point t.

Theorem-Definition 4.1 ([47, Theorem-definition 5.1]). Let h : R → (0, 1) be a C 1 deterministic function such that its derivative function h ′ is bounded. The process W h defined by ( 19) is an (S) * -process which verifies the following equality (in (S) * ):

W h t = +∞ k=0 M h(t) (e k )(t) •, e k + h ′ (t) +∞ k=0 Å t 0 ∂MH ∂H (e k )(s) H=h(t) ds ã •, e k . ( 20
)
Moreover the process B h is (S) * -differentiable on R and verifies dB h dt (t) = W h t in (S) * . When the function h is constant, identically equal to H, we will denote W H := W H t t∈R and call the process W h fractional white noise. This process was defined and studied in [START_REF] Elliott | A general fractional white noise theory and applications to finance[END_REF][START_REF] Biagini | An introduction to white-noise theory and Malliavin calculus for fractional Brownian motion[END_REF].

Stochastic integral with respect to mBm

Using the tools presented above, we can now recall the definition of the Wick-Itô integral with respect to mBm from [START_REF] Lebovits | White noise-based stochastic calculus with respect to multifractional Brownian motion[END_REF]. Theorem 4.3 in this section will be instrumental to solve the S.D.E. encountered later Definition 4 (Multifractional Wick-Itô integral). Let B h be a normalized multifractional Brownian motion and Y : R → (S) * be a process such that the process t → Y t ⋄ W h t is (S) * -integrable on R. The process Y is said to be d ⋄ B h -integrable on R or integrable on R with respect to mBm B h . Moreover, the integral on R of Y with respect to B h is defined by

R Y s d ⋄ B h s := R Y s ⋄ W h s ds. ( 21
)
For an interval

I of R, I Y s d ⋄ B h s := R 1 I (s) Y s d ⋄ B h s .
When the function h is constant over R, equal to H, the multifractional Wick-Itô integral coincides with the fractional Wick-Itô integral defined in [START_REF] Elliott | A general fractional white noise theory and applications to finance[END_REF], [START_REF] Biagini | An introduction to white-noise theory and Malliavin calculus for fractional Brownian motion[END_REF], [START_REF] Bender | An Itô formula for generalized functionals of a fractional Brownian motion with arbitrary Hurst parameter[END_REF] and [START_REF] Bender | An S-transform approach to integration with respect to a fractional Brownian motion[END_REF]. In particular, when Y is adapted and when the function h is identically equal to 1/2, (21) is nothing but the classical Itô integral with respect to Brownian motion.

The multifractional Wick-Itô integral verifies the following properties: Proposition 4.2. Let B h be an mBm and I be an interval of R.

• For all (a, b) in R 2 such that a < b, b a 1 d ⋄ B h u = B h b -B h a almost surely.
• Let X : I → (S) * be a d ⋄ B h -integrable process over I.

If I X s d ⋄ B h s belongs to (L 2 ), then E[ I X s d ⋄ B h s ] = 0.

Multifractional Wick-Itô integral of deterministic elements

In order to solve differential equations driven by an mBm that will be encountered below, it is necessary to know the exact nature of multifractional Wick-Itô integrals of deterministic elements.

For H in (0, 1) and f in S (R), define the function g f : R × (0, 1) → R by g f (t, H) := t 0 M H (f )(x)dx where M H is the operator defined in at the beginning of Section 4.1.1. It has been shown that g f belongs to C ∞ (R × (0, 1), R) (cf. [START_REF] Lebovits | White noise-based stochastic calculus with respect to multifractional Brownian motion[END_REF]Lemma 5.5]). The main result on the multifractional Wick-Itô integral of deterministic elements is the following: Theorem 4.3. ([47, Theorem 5.25]) Let h : R → (0, 1) be a C 1 deterministic function and let f : R → R be a measurable deterministic function which belongs to L 1 loc (R). Let Z := (Z t ) t∈R be the process defined by

Z t := t 0 f (s) d ⋄ B h s .
Then Z is an (S) * -process which verifies the following equality in (S)

* t 0 f (s) d ⋄ B h s = +∞ k=0 Ç t 0 f (s) d ds [g e k (s, h(s))] ds å •, e k . ( 22
)
Moreover Z is a (centred) Gaussian process if and only if

+∞ k=0 Ä t 0 f (s) d ds [g e k (s, h(s))] ds ä 2 < +∞, for all t in R.
In this case we have, for every t in R,

Z t = t 0 f (s) d ⋄ B h s L ∼ N 0, +∞ k=0 Ç t 0 f (s) d ds [g e k (s, h(s))] ds å 2 . ( 23
)
In particular, Z is a Gaussian process when f belongs to C 1 (R; R).

Deriving the quantity E Z 2 t in the previous theorem might be complicated using Equation [START_REF] Corlay | Functional quantization-based stratified sampling methods[END_REF]. However, when f is a C 1 function, thanks to the Itô formula with respect to mBm [47, Theorem 6.9], we obtain the following integration-by-parts formula

t 0 f (s) d ⋄ B h s (L 2 ) = f (t) B h t - t 0 f ′ (s) B h s ds, (24)
which leads to

E Z 2 t = f (t) 2 t 2h(t) + t 0 t 0 f ′ (s) f ′ (u) R h (s, u) ds du -2f (t) t 0 f ′ (s) R h (t, s) ds. (25)
Remark 6. The integration-by-parts formula (24) allows to identify almost surely

t 0 f (s) d ⋄ B h s with the quantity I f t (B h ) where the map I f t : C 0 ([0, t]; R) → R is defined by I f t : g → Ç f (t)g(t) - t 0 f ′ (s)g(s) ds å . (26)

Stochastic differential equations

We solve in this subsection the two stochastic differential equations that define the multifractional stochastic volatility models presented in Section 5.

Mixed multifractional Brownian S.D.E.

Let us consider the following mixed multifractional stochastic differential equation, where γ 1 and γ 2 are positive constants and B t is a Brownian motion:

® dX t = X t γ 1 d ⋄ B t + γ 2 d ⋄ B h t , X 0 = x 0 ∈ R. (27)
Of course ( 27) is a shorthand notation for the equation

X t = x 0 + γ 1 t 0 X s d ⋄ B s + γ 1 t 0 X s d ⋄ B h s , X 0 = x 0 ∈ R,
where the previous equality holds in (S) * . A solution of this equation will be called geometric mixed multifractional Brownian motion. Rewriting the previous equation in terms of derivatives in (S) * , we get:

dX t dt = X t ⋄ Ä γ 1 W 1/2 t + γ 2 W h t ä , x 0 ∈ R. (28)
Theorem 4.4 (Geometric mixed multifractional Brownian motion). The (S) * -process (X t ) t∈[0,T ] defined by

X t := x 0 exp ⋄ γ 1 B t + γ 2 B h t , (29)
is the unique solution of (28) in (S) * .

Proof: Applying the S-transform to both sides of Equation ( 28) and denoting by y η the map t → S(X t )(η), for every η in S (R), we get:

y ′ η (t) = y η (t) Å γ 1 M 1/2 (η)(t) + γ 2 d dt [g η (t, h(t))] ã , y η (0) = x 0 .
This equation admits a unique solution which verifies y η (t) = x 0 exp {γ

1 t 0 M 1/2 (η)(u)du + γ 2 t 0 d du [g η (u, h (u)
)]du}. Using (i) and (ii) of [47, Theorem 5.12] we hence get, for every η in S (R),

y η (t) = x 0 exp {γ 1 S(B t )(η) + γ 2 S(B h t )(η)} = S x 0 exp ⋄ {γ 1 B t + γ 2 B h t } (η).
The injectivity of the S-transform allows us to conclude that [START_REF] Janson | Gaussian Hilbert Spaces[END_REF]Equality (3.16)], one sees that X is an (L 2 )-valued process that may be represented as:

X t = x 0 exp ⋄ γ 1 B t + γ 2 B h t for every t in [0, T ]. Remark 7. (i) Using
X t = x 0 exp ¶ γ 1 B t + γ 2 B h t -1 2 Ä γ 2 1 t + γ 2 2 t 2h(t)
ä© .

(ii) The theorem is also a consequence of [41, Theorem 3.1.2].

Mixed multifractional Ornstein-Uhlenbeck S.D.E.

Let us now consider the following mixed stochastic differential equation:

® dU t = θ(µ -U t )dt + (α 1 d ⋄ B t + α 2 d ⋄ B h t ) U 0 = u 0 ∈ R, (30)
where (B t ) t∈R and B h t t∈R are independent, θ ≥ 0 and µ, α 1 , α 2 belong to R. A solution of this equation will be called a mixed multifractional Ornstein-Uhlenbeck process. Theorem 4.5 (Mixed multifractional Ornstein-Uhlenbeck process). The L 2 (Ω)-valued process (U t ) t∈R defined by

U t := u 0 e -θt + µ 1 -e -θt + α 1 t 0 e θ(s-t) d ⋄ B s + α 2 t 0 e θ(s-t) d ⋄ B h s , (31)
is the unique solution of the stochastic differential equation [START_REF] Echelard | Terrain modeling with multifractional brownian motion and self-regulating processes[END_REF].

Proof: The proof that the process U defined by [START_REF] Elliott | A general fractional white noise theory and applications to finance[END_REF] is the unique solution of ( 30) is very similar to the one of Theorem 4.4. Indeed, setting y η (t) := S(U t )(η) for every (t, η) in R × R and applying the S-transform to both sides of (30) we get, for every η in S (R), the ordinary differential equation

y ′ η (t) = θ(µ -y η (t)) + α 1 M 1/2 (η)(t) + α 2 d dt [g η (t, h(t))], y η (0) = u 0 . (32)
Its unique solution is

y η (t) = u 0 e -θt + e -θt t 0 e θs θµ + α 1 M 1/2 (η)(s) + α 2 d ds [g η (s, h(s))] ds, y η (0) = u 0 .
Again, one concludes using the injectivity of the S-transform.

Quantization of solutions of S.D.E. driven by mBm

Quantizing a Gaussian process X often yields as well a quantization of the solutions of stochastic differential equations driven by X: indeed, in many cases, these solutions may be expressed as functionals of X. A quantizer of the solution can then be obtained by simply plugging the quantizer of X into the functional.

In the one-dimensional setting, under rather general conditions on the diffusion coefficients and if X is a continuous semimartingale, this functional is easily determined using the Lamperti transform (see [START_REF] Luschgy | Functional quantization of a class of Brownian diffusions: a constructive approach[END_REF]). In this case, the corresponding quantizer of the stochastic differential equation is obtained by plugging the Gaussian quantizer in the S.D.E. written in the Stratonovich sense, leading to a finite set of ordinary differential equations. This leads to a simple and general constructive method to build a functional quantization of the solution of an S.D.E.

Unfortunately, no such result is available in the case of an S.D.E. driven by multifractional Brownian motion (or even by fractional Brownian motion). However, in some situations, and in particular when an explicit solution is known, one may sometimes still use the procedure just described: if the functional giving the solution is regular enough, quantization-based cubatures can then be used. This is for instance the case of geometric mixed multifractional Brownian motion defined in Section 4.3.1, which is a simple functional of a Brownian motion and a multifractional Brownian motion (see Remark 7 and section 3.3.2). We describe two other favourable situations in the next subsections.

The case of a Wiener integral

An easy case is the one of a Wiener integral

t 0 f (s)d ⋄ B h
s where f is a C 1 deterministic function. The integration-by-parts formula for mBm [START_REF] Coutin | An introduction to (stochastic) calculus with respect to fractional Brownian motion[END_REF]

reads t 0 f (s)d ⋄ B h s a.s. = f (t)B h t - t 0 f ′ (s)B h s ds. Thus, for p ≥ 1, the stochastic process t → t 0 f (s)d ⋄ B h
s , seen as a random variable valued in L p (0, T ), is the image of B h by the map

J f : L p ([0, T ]) → L p ([0, T ]) g → f (•)g(•) - • 0 f ′ (s)g(s)ds.
In other words we have (dt-almost everywhere) J f (g)(t) = I f t (g) where I f t defined in Remark 6.

Proposition 4.6 (L p -regularity of the Wiener map). For every p ≥ 1, the map J f is Lipschitz continuous on L p ([0, T ]).

Proof: It is straightforward that for (g 1 , g 2 ) ∈ L p ([0, T ]) 2 J f (g 1 ) -J f (g 2 ) p ≤ f (g 1 -g 2 ) p + • 0 f ′ (s)(g 1 (s) -g 2 (s))ds p ≤ ( f ∞ + f ′ ∞ T ) g 1 -g 2 p .
In Appendix A, we prove that if h is C 1 , the Karhunen-Loève eigenfunctions of a well-balanced mBm B h have bounded variations, and thus stationary quantizers of B h have bounded variations as well (because they lie on a subspace of L 2 ([0, T ]) spanned by a finite number of Karhunen-Loève eigenfunctions, as already mentioned). In this setting, another integration by parts gives

I f t Ä " B h ä = t 0 f (s)d " B h
s where d " B h s (ω) stands for the signed measure associated with the function of bounded variation s → " B h s (ω).

The case of certain simple diffusions

Another easy case is the one of an S.D.E. of the form

Y t = y 0 + t 0 β(s, Y s )ds + X t , (33)
where β(s, y) is assumed to be Lipschitz continuous in y uniformly in s. This setting is addressed in [52, p. 20-21] (actually, in [START_REF] Luschgy | Functional quantization of a class of Brownian diffusions: a constructive approach[END_REF], the Lamperti transform is used to reduce a general Brownian diffusion to this case), where the authors consider the associated integral equation

y(t) = y 0 + t 0 β(s, y(s))ds + g(t), (34)
where g ∈ L p ([0, T ]) is fixed. The existence and uniqueness in L p ([0, T ]) of a solution for the integral equation [START_REF] Graf | Functional quantization and small ball probabilities for Gaussian processes[END_REF] follows from the same approach used for ordinary differential equations. Then the solution of the associated S.D.E. [START_REF] Falconer | Multifractional, multistable, and other processes with prescribed local form[END_REF] simply reads U t = Ψ β p (X) t , where Ψ β p : L p ([0, T ]) → L p ([0, T ]) is the functional that maps g ∈ L p ([0, ]) to the unique solution in L p ([0, T ]) of Equation [START_REF] Graf | Functional quantization and small ball probabilities for Gaussian processes[END_REF]. In [START_REF] Luschgy | Functional quantization of a class of Brownian diffusions: a constructive approach[END_REF], the map Ψ β p is showed to be Lipschitz continuous in L p ([0, T ]). More precisely, one has

c([β] Lip , T ) g 1 -g 2 p p ≤ Ψ β p (g 1 ) -Ψ β p (g 2 ) p p ≤ C([β] Lip , T ) g 1 -g 2 p p , with c([β] Lip , T ) = 1 2 p-1 (1-[β] p Lip T p ) and C([β] Lip , T ) = e 2 p-1 [β] Lip T p-1 .
Mixed multifractional Ornstein-Uhlenbeck process, defined in Section 4.3.2, is of the form [START_REF] Graf | Functional quantization and small ball probabilities for Gaussian processes[END_REF], with β(s, u) = θ(µu) and X = α 1 B + α 2 B h .

Multifractional stochastic volatility models

Subsection 5.1 is devoted to a short recall on implied forward start volatility. In Subsections 5.2 and 5.3 we propose two multifractional stochastic volatility models, the multifractional Hull & White and the multifractional SABR volatility models. The former generalizes the fractional long-memory stochastic volatility model presented in [START_REF] Comte | Long memory in continuous-time stochastic volatility models[END_REF]Paragraph 2] and the latter is an extension of the celebrated SABR stochastic volatility model [START_REF] Hagan | Managing smile risk[END_REF] in the case where β is equal to 1. We propose a numerical scheme based on the functional quantization of mBm for the computation of the price of forward start options in the two cases. For any ε > 0 we provide an upper bound for the L 2-ε quantization error of the instantaneous volatility process σ in the two models.

The implied forward start volatility

Vanilla option prices are typically converted in terms of Black Scholes implied volatility by practitioners, because this quantity can be easily interpreted as the crude option price. As we devised a numerical scheme for forward start options, we give here the associated notion of "forward implied volatility", with a special attention because the term "forward volatility" is used for different notions in the literature.

The forward start option price.

Let W be a standard Brownian motion on [0, T ] and τ ∈ (0, T ). Let us consider the stochastic differential equation dS t = S t σ t dW t (with (σ t ) t∈[0,T ] a deterministic process) whose solution is a geometric Brownian motion S t = S 0 exp Ä t 0 σ s dW s -1 2 t 0 σ 2 s ds ä . The forward start Call option price FSPrimeBS(σ, τ, T, K) is given by

FSPrimeBS(σ, τ, T, K) = E ñ Å S T S τ -K ã + ô = N (d 1 ) -KN (d 2 ),
where

d 1 := σ √ T -τ 2 + ln(K) σ √ T -τ , d 2 := d 1 -σ √ T -τ 2 and σ 2 := 1 T -τ T τ σ 2 s ds.
In other words, we have

FSPrimeBS(σ, τ, T, K) = PrimeBS(1, σ, T -τ, K), (35)
where (S 0 , Vol, Mat, Strike) → PrimeBS(S 0 , Vol, Mat, Strike) is the closed-form expression for the vanilla Call option price in the Black & Scholes model.

The implied forward start volatility.

In the Black & Scholes model, where the asset price follows a geometric Brownian motion with a constant volatility, the forward start Call (or Put) option price is an increasing function of the volatility (if the strike is not zero). Conversely, for a given forward start Call (or Put) option price, the Black & Scholes implied volatility is the unique value of the volatility for which the Black & Scholes formula recovers the price; in other words, the implied forward start volatility associated with a given forward F 0 , a forward start date τ , a maturity T > τ , a strike K, and an option price P is defined by P = FSPrimeBS (ImpliedFSVolBS (τ, T, K, P ) , τ, T, K) .

(

) 36
Using Equation [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF], this yields

ImpliedFSVolBS (τ, T, K, P ) = ImpliedVolBS (1, T -τ, K, P ) , (37)
where ImpliedVolBS (Fwd, maturity, Strike, Price) is the Black & Scholes implied volatility a certain forward, maturity, strike and option price.

Multifractional Hull & White stochastic volatility model

We assume that, under the risk-neutral measure, the forward price of a risky asset is the solution of the S.D.E.

ß dF t = F t σ t dW t , d ln(σ t ) = θ (µ -ln(σ t )) dt + γ h d ⋄ B h t + γ σ dW σ t , σ 0 > 0, (38)
where θ ≥ 0 and where W and W σ are two standard Brownian motions and B h is a well-balanced multifractional Brownian motion independent of W and W σ with functional parameter h assumed to be continuously differentiable. We assume that W is decomposed into ρdW σ t + 1ρ 2 dW F t , where W F is a Brownian motion independent of W σ . Hence, [START_REF] Herbin | Stochastic integration with respect to multifractional Brownian motion via tangent fractional Brownian motion[END_REF] writes

® dF t = F t σ t Ä ρdW σ t + 1 -ρ 2 dW F t ä d ln(σ t ) = θ (µ -ln(σ t )) dt + γ h d ⋄ B h t + γ σ dW σ t , σ 0 > 0. (39)
We denote respectively by F σ , F F and F h the natural filtrations of W σ , W F and B h . We define the filtration F σ,h by F σ,h t = σ F σ t , F h t and F F,σ,h by F F,σ,h t = σ F F t , F σ t , F h t . The unique solution of [START_REF] Herbin | Stochastic integration with respect to multifractional Brownian motion via tangent fractional Brownian motion[END_REF] reads

® F t = F 0 exp Ä t 0 σ s dW s -1 2 t 0 σ 2 s ds ä σ s = exp ln(σ 0 )e -θs + µ 1 -e -θs + γ σ s 0 e θ(u-s) dW σ u + γ h s 0 e θ(u-s) d ⋄ B h u . (40)
In other words, ln(σ t ) is a mixed multifractional Ornstein-Uhlenbeck process. Note that, although the volatility process is not a semimartingale, the process (F t ) t∈[0,T ] remains a (positive) F F,σ,h -local martingale, and thus a super-martingale. The same proof as in [START_REF] Jourdain | Loss of martingality in asset price models with lognormal stochastic volatility[END_REF] shows that, if ρ = 0, this local martingale is indeed a martingale. Numerical experiments seem to indicate that this property still holds for ρ < 0, a fact that remains to be proved. We now consider the problem of pricing a forward start call option (the put case is handled similarly). The payoff of this option writes Ä The following decomposition holds:

F t = F 0 exp Ç ρ t 0 σ s dW σ s - ρ 2 2 t 0 σ 2 s ds å measurable with respect to F σ,h t exp Ç 1 -ρ 2 t 0 σ s dW F s - 1 -ρ 2 2 t 0 σ 2 s ds å . Conditioning by F σ,h T yields E Ä FT Fτ -K ä + = E E Ä FT Fτ -K ä + F σ,h T = E E Ä F τ,T exp Ä 1 -ρ 2 t 0 σ t dW F t -1-ρ 2 2 t 0 σ 2 s ds ä -K ä + F σ,h T = E PrimeBS F τ,T , (1 -ρ 2 ) 1 T -τ T τ σ 2 s ds =:I σ τ,T 1 2 , T -τ, K , (41)
where s ds as explicit functionals of the paths of W σ and B h in L 2 ([0, T ]). Recall that σ is the exponential of a mixed multifractional Ornstein-Uhlenbeck process:

F τ,T := exp Ä ρ T τ σ s dW σ s -ρ 2
σ t = exp Ä ln(σ 0 )e -θt + µ 1 -e -θt + γ σ e -θt I e θ• t (W σ ) + γ h e -θt I e θ• t (B h ) ä . (42)
This yields an explicit functional form for T τ σ 2 s ds as a function of the paths of W σ and B h . Denote (p h j ) 1≤j≤N1 and (χ h j ) 1≤j≤N1 the weights and the paths of the quantizer " B h of B h , and (p σ j ) 1≤j≤N2 and (χ σ j ) 1≤j≤N2 the weights and the paths of the quantizer W σ of W σ . Conditionally on B h = χ h i , one has I σ τ,T = I σ i τ,T , where

I σ i τ,T := T τ σ i s dW σ s and σ i t = exp ln(σ 0 )e -θt + µ 1 -e -θt + γ σ t 0 e θ(s-t) dW σ s + γ h e -θt I g t χ h i .
Appendix B shows that χ h i has bounded variations. This entails that σ i is a semimartingale. Define σ i , W σ τ,T := σ i , W σ Tσ i , W σ τ , where •, • denotes the semimartingale bracket and let us denote by T τ σ i s • dW σ s the Stratonovich integral of σ i . Then, I σ i τ,T reads

I σ i τ,T = T τ σ i s • dW σ s - 1 2 σ i , W σ τ,T .
Itô's formula yields

T τ σ i t dW σ t = σ i T -σ τ γ σ - 1 γ σ T τ σ i t θ µ -ln σ i t dt - γ h γ σ T τ σ i t dχ h i (t) = T τ σ i t • dW σ t - γ σ 2 T τ σ i t dt = 1 2 σ i ,W σ τ,T . Moreover, T τ " σ i t d W σ t = " σ i T -σ τ γ σ - 1 γ σ T τ " σ i t θ Ä µ -ln Ä " σ i t ää dt - γ h γ σ T τ " σ i t dχ h i (t).
This shows that

:= T τ " σ i s d W σ s -γσ 2
T τ " σ i s ds. The cubature formula is then fully explicit and one finally obtains the following approximation:

E ñ Å F T F τ -K ã + ô ≈ N1 i=1 N2 j=1 p h i p σ j PrimeBS Ñ F i,j τ,T , Ç 1 -ρ 2 1 T -τ T τ σ i,j (s) 2 ds å 1 2 , T -τ, K é , where F i,j τ,T = exp Ç ρ T τ σ i,j (s)dχ σ j (s) -ργ σ 1 2 T τ σ i,j (s)ds - ρ 2 2 T τ σ i,j (s) 2 ds å , and
σ i,j (t) := exp Ä ln(σ 0 )e -θt + µ 1 -e -θt + γ σ e -θt I e θ• t (χ σ j ) + γ h e -θt I e θ• t (χ h j ) ä .
ln(σ i,j ) 1≤i≤N1,1≤j≤N2 and (p h i p σ j ) 1≤i≤N1,1≤j≤N2 are the paths and weights of a stationary quantizer of the mixed multifractional Ornstein-Uhlenbeck process ln(σ). The results stated in Section 4.4.2 allow us to control its quadratic quantization error with the quantization error of W σ and " B h . We then apply (13) to get an upper bound for the L 2-ǫ quantization error of the process σ on [0, T ], for any ǫ > 0 .

Multifractional SABR model

We consider the case where, under the risk-neutral measure, the forward price of a risky asset is the solution of the S.D.E.

ß dF t = F t σ t dW t , dσ t = σ t γ h d ⋄ B h t + γ σ dW σ t , σ 0 > 0, (43)
where W and W σ are two standard Brownian motions and B h is a well-balanced multifractional Brownian motion independent of W and W σ with functional C 1 parameter h. We assume that W is decomposed into ρdW σ t + 1ρ 2 dW F t , where W F is a Brownian motion independent of W σ . We use the same notations as in the previous section for F σ , F F , F h , F σ,h and F F,σ,h . Hence, [START_REF] Janson | Gaussian Hilbert Spaces[END_REF] writes

® dF t = F t σ t Ä ρdW σ t + 1 -ρ 2 dW F t ä dσ t = σ t γ h d ⋄ B h t + γ σ dW σ t , σ 0 > 0. (44)
This is an extension of the SABR model, in when the β parameter is equal to 1. This model can be handled in the same way as the multifractional Hull & White model. The solution of the stochastic differential equation verified by σ, established in Theorem 4.4, is

σ t = σ 0 exp ⋄ γ σ W σ t + γ h B h t = σ 0 exp Å γ σ W σ t + γ h B h t - 1 2 Ä γ 2 σ t + γ 2 h t 2h(t) ä ã . ( 45
)
Reasoning as in the case of the Hull & White model presented in Section 5.2, it can be shown that F is an F F,σ,h -martingale for ρ = 0. In addition, the same numerical procedures as above may be used.

6 Numerical experiments

Variance reduction method for the quantization-based cubature

Here, we present a simple kind of control variate method for the quantization-based cubature method that we use, which slightly improves the numerical accuracy of the method. Numerical experiments carried out in [START_REF] Corlay | Some aspects of optimal quantization and applications to finance[END_REF] showed that, in the case of vanilla options, computing the implied volatility using the estimated forward instead of the theoretical forward in the Black & Scholes formula improves the accuracy. The counterpart of this method in the frame of forward start options is to replace the "1" appearing in Formula (37) by the quantity

I N1,N2 := 1≤i≤N1, 1≤j≤N2 p h i p σ j F i,j τ,T (46)
This also holds when using Richardson-Romberg extrapolation: in this case, one uses the extrapolated value of I N1,N2 instead of 1 in Formula (37). These methods were used to generate the numerical results presented below.

Numerical results

We present results on the multifractional Hull & White model. We have computed the price as a function of strike for different maturities: 1, 2.5, 5 and 10 years. Driving noises were chosen in the class of fBms and mBms. More precisely, we display results of our experiments with:

1. An fBm with H = 0.2.

2. An fBm with H = 0.5.

3. An fBm with H = 0.75.

4. An fBm with H = 0.9. 9. And, finally, an mBm with h = ◊ h V olSP , which corresponds to the regularity estimated on the S&P 500 trace that was analysed in Section 1.

An mBm with h

= h 1 = 0.
The four functions h 1 to h 4 are plotted on Figure 6. The values of the other parameters are γ h = 0.3, γ σ = ρ = 0 (except for the experiments displayed on Figure 9), θ = 0.3, µ = ln(0.2), σ 0 = 0.2 and F 0 = 100. The results displayed on Figures 7 and8 provide an experimental justification to the claims made in the introduction. Indeed, one sees that, for the short maturity T = 1 year, in the fractional Hull & White model (i.e. with h constant), the smiles are more pronounced for small H and decrease as H increase, while the reverse is true for all maturities larger than one year (Figure 7). Thus, stronger correlations in the driving noise do translate in this model into a slower decrease of the smile as maturities increase, as noted in [START_REF] Comte | Affine fractional stochastic volatility models[END_REF]. However, with such an fBm-based model, an H larger than 1/2 is needed to ensure long-range dependence and thus a more realistic evolution of the smile as compared to the Brownian case. As mentioned above, this is not compatible with empirical graphs of the volatility which show a very irregular behaviour, and would require a small H, a fact which was confirmed in Section 1 through an estimation of the local regularity. In addition, the local regularity of the volatility evolves in time, calling for a varying H, i.e. an mBm.

Another aspect is that a fixed H, as in a modeling with fBm, does not allow to control independently the shape of the smiles at different maturities. This is possible with mBm, where the smile at maturity T depends on a weighted average of the values of h up to time T , as can be inferred from equalities [START_REF] Coutin | An introduction to (stochastic) calculus with respect to fractional Brownian motion[END_REF] and [START_REF] Heston | A closed-form solution for options with stochastic volatility, with application to bond and currency options[END_REF]. This is apparent on Figure 8. We have compared fBms and mBms at various maturities T , where H and h are chosen such that h(t) = H, or, for the bottom right plot, h 1 (t) = h 4 (t). One sees that the shape of the smile depends on a weighted average of past values of h. For instance, in the bottom left plot, the values of h before T = 2.5 are in average smaller than 0.9, resulting in a flatter smile. The fact that a weighted average must be considered is apparent on the bottom right plot: indeed, the smile is more pronounced for h 1 , although the average from 0 to 5 this function is smaller than the one of h 4 . In contrast, the values in the immediate past of t = 5 are larger for h 1 than for h 4 , as may be checked on Figure 6. An adequate choice of h may thus allow one to better approximate a whole implied volatility surface. This topic will be addressed in a future work.

Figure 9 displays:

• an example with ◊ h V olSP , that regresses the regularity function of the volatility of S&P 500 estimated in Section 1,

• an example with ρ = 0 and h = h 2 .

As one can see, smiles computed with a regularity function obtained from market data indeed display all the features observed empirically as detailed above. We believe this provides a further justification to the relevance of our model. Top left: fBm with H = 0.9 and mBm with function h 1 at T = 5 (h 1 (5) = 0.9). Top right: fBm with H = 0.2 and mBm with function h 2 at T = 5 (h 2 (5) = 0.2). Bottom left: fBm with H = 0.9 and mBm with function h 3 at T = 2.5 (h 3 (2.5) = 0.9). Bottom right: mBm with function h 1 and mBm with function h 4 at T = 5 (h 1 (5) = h 4 (5) = 0.9).

We show that F i has bounded variations for every i in {1, 2, 3}. The cases of F 1 , F 2 and F 3 are similar, and we only treat here F 1 . Let (t i ) 0≤i≤N be a sequence of elements of [0, T ] such that 0 = t 0 < t 1 < x 2H2-1x 2H1-1 dx.

Define the family of maps (g α ) α∈R * + from R + to R + , by g α (x) := α x if x > 0 and g α (x) := 1 if x = 0. Let 49) and [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF] we finally obtain

N i=1 |F 1 (t i ) -F 1 (t i-1 )| ≤ 2K 7 Å 1 + 1 2H 1 ã T + T 2H1 + T 2H2 < +∞,
which ends the proof.

Figure 1 :

 1 Figure 1: S&P 500 minute data (left) and estimated volatility for time periods of four hours (right).

Figure 2 :

 2 Figure 2: Estimated regularity of the volatility of the S&P 500 minute data (blue) and its regression (green).

Figure 3 :

 3 Figure 3: VIX minute data (left) and estimated regularity for time periods of four hours (right).

Figure 4 :

 4 Figure 4: Quadratic 5 × 2 × 2-product quantizer of fBm on [0, 1] with H = 0.25 (left) and H = 0.75 (right).

Figure 5 :

 5 Figure 5: Quadratic 5 × 2 × 2-product quantizer of mBm on [0, 1] with h(t) := 0.1 + 0.8t (left) and h(t) := 0.9 -0.8t (right).

.

  Define the map φ : C([0, T ], R) → C([0, T ], R) by φ(f ) := T 0 e 2p p|f (s)| ds. It is easily shown that φ is convex and continuous on (C([0, T ], R), • ∞ ). Hence, Inequality (12) yields E

2 2 and E D 2 F

 22 (Y N )(X -Y N ) ⊗2 has been established yet. Still, Equation (14) suggests to use a Richardson-Romberg extrapolation with respect to the quantization error E 2 N := X -Y N 2 . The two-steps extrapolation between N = k and N = l leads to approximate E[F (X)] by the quantity

x 2 )(2k + 2 ) 2p f, e k 2 L 2 (

 2222 the nth Hermite function. Let (| • | p ) p∈Z be the family of norms defined by |f | 2 p := +∞ k=0 R) , for all (p, f ) in Z × L 2 (R). The operator A defined on S (R) by A := -d 2 dx 2 + x 2 + 1 admits the sequence (e n ) n∈N as eigenfunctions and the sequence (2n + 2) n∈N as eigenvalues.

FT Fτ -K ä +FT

 ä for some fixed maturity τ ∈ [0, T ]. We need to compute the risk-neutral expectation E Ä

2 T τ σ 2

 22 s ds ä and PrimeBS is the closed-form expression for the price of a Call option in the Black & Scholes model, detailed in Appendix 5.1. The aim is to estimate the expectation (41) by a quantization-based cubature associated with the functional quantization of B h and W σ . We thus need to write the terms F τ,T and T τ σ 2

3 Figure 6 :

 36 Figure 6: left: functions h 1 and h 4 ; middle: function h 2 ; right: function h 3 .

Figure 7 :

 7 Figure 7: Comparisons of vanilla option volatility smiles for fBm with H = 0.2, H = 0.5, H = 0.7 and H = 0.9 at different maturities. Top left: T = 1. Top right: T = 2.5. Bottom left: T = 5. Bottom right: T = 10.

  0

Figure 8 :

 8 Figure8: Comparisons of vanilla option volatility smiles for various fBm and mBm at several maturities. Top left: fBm with H = 0.9 and mBm with function h 1 at T = 5 (h 1 (5) = 0.9). Top right: fBm with H = 0.2 and mBm with function h 2 at T = 5 (h 2 (5) = 0.2). Bottom left: fBm with H = 0.9 and mBm with function h 3 at T = 2.5 (h 3 (2.5) = 0.9). Bottom right: mBm with function h 1 and mBm with function h 4 at T = 5 (h 1 (5) = h 4 (5) = 0.9).

Figure 9 :

 9 Figure 9: Vanilla option volatility smiles in the multifractional Hull & White model, with γ h = 0, γ σ = 0.3, θ = 0.3, µ = ln(0.2), σ 0 = 0.2 and F 0 = 100, and h = ◊ h V olSP for maturities T = 1, T = 2.5 and T = 5 (left), and with γ h = 0.3, γ σ = 0.3 ρ = -0.5, θ = 0.3, µ = ln(0.2), σ 0 = 0.2 and F 0 = 100, and h = h 2 for maturities T = 1, T = 2.5 and T = 5 (right).

  K 5 := sup α∈[0,T ]| ln(α)| (e 2H1 ln(α) + e 2H2 ln(α) ). The mean-value theorem applied to g α yields∀s ∈ [0, T ], |D i (s)| ≤ 2 -1 K 5 |2h ti,s -2h ti-1,s | ≤ K 5 sup u∈[0,T ] |h ′ (u)| |t it i-1 | =: K 6 |t it i-1 |.We hence have shown that∀i ∈ {1; • • • ; N }, G i ≤ =:K7 (1 + T ) (1 + 2H 2 ) (1 + K 4 ) (1 + K 6 ) Ç |t it i-1 | + ti ti-1 x 2H2-1x 2H1-1 dx å .

  • • • < t N = T . For any i in {1, • • • , N } we get, |F 1 (t i ) -F 1 (t i-1 )| ≤ ds ≤ T (1 + K 2 ) |t it i-1 | e 2H1T + e 2H2T =: K 3 |t it i-1 |.

						=:K1
			sup s∈[0,T ]	e k (s) c h(s)		0	T	c 2 ht i ,s c h(ti)	t	2ht i ,s i	-	c 2 ht i-1 ,s c h(ti-1)	t	2ht i-1 ,s i-1	ds
		≤ K 1	Å T 0	c 2 ht i ,s c h(ti)	t	2ht i ,s i	-t	2ht i-1 ,s i-1	ds	+	0	T	c 2 ht i ,s c h(ti)	-	c 2 ht i-1 ,s c h(ti-1)	t	2ht i-1 ,s i-1	ds	ã	. (48)
														=:Gi			=:Li
	Since the map (s, t) →	c 2 h t,s
				c 2 ht i ,s c h(ti)	-	c 2 ht i-1 ,s c h(ti-1)		≤ sup s∈[0,T ]
	ï	ò																	c 2 h t,s c h(t) . Setting [H 1 , H 2 ] :=
	inf u∈[0,T ] h(u), sup u∈[0,T ] h(u)		, one gets:
	L T 0		t	2ht i-1 ,s i-1										(49)
	Besides, G i ≤	sup (t,s)∈[0,T ] 2	c 2 h t,s c h(t)		T 0	t i	2ht i ,s -t	2ht i-1 ,s i-1	ds =: K 4	T 0	t	2ht i ,s i	-t	2ht i-1 ,s i-1	ds.
	Now, writing
				t i	2ht i ,s -t	2ht i-1 ,s i-1	= t i	2ht i ,s -t	2ht i ,s i-1	-t	2ht i ,s i-1 -t	2ht i-1 ,s i-1	,
																			:=Ci(s)	:=Di(s)
	we easily get that
																				ti
				∀s ∈ [0, T ], |C i (s)| ≤ 2H 2	ti-1

c h(t) is C 1 as soon as h is C 1 , the mean-value theorem yields |f ′ s (t)| |t it i-1 | =: K 2 |t it i-1 |,

where f ′ s (t) denotes, for every s in [0, T ], the derivative, at point t, of the map t → i ≤ K 2 |t it i-1 |

Acknowledgments

The authors are thankful to F. Abergel, B. Dupire, G. Pagès and M. Yor for many helpful remarks and comments.