Susanne Fuchs
email: fuchs@zas.gwz-berlin.de

Ralf Winkler

Pascal Perrier

Do Speakers' Vocal Tract Geometries Shape their Articulatory Vowel Space?

This study investigates the relation between parameters describing differences between speaker-specific vocal tract geometries and articulatory distances between the corner vowels based on MRI data of 9

French speakers. Results provide evidence that speaker with a longer pharynx produce larger displacements between low back and high front vowels. Preliminary modeling results are also presented with the aim to study the relation between motor commands, articulation and acoustics.

Introduction

Speaker-specific articulatory behaviour seems to be an inherent characteristic of speech production. We are particularly interested in the question where inter-speaker variability stems from and draw our attention to the potential inter-speaker differences in vocal tract morphology. One of the difficulties in such an undertaking is to disentangle between the potential factors making one speaker distinguishable from another. Particularly, with regard to articulatory control and its variability, it is quite challenging to tease apart whether two speakers differ, since this is the intended goal (sociolinguistic and dialectal factors may play a role) or whether the differ due to differences in their vocal tract morphology.

The approach we favour to deal with such an issue is to gather experimental data for a large set of speakers and compare them with simulations carried out with realistic physical models of their vocal tract. Modelling has the advantage that different factors underlying the produced output can be controlled separately. However, it has to be considered carefully since it is based on hypotheses, estimations, theories how reality works, but it is not a copy of the reality.

Hypothesis and assumptions

The hypothesis and assumption we will make, are based on previous studies, which are briefly outlined below:

Honda and colleagues [START_REF] Honda | Human palate and related structures: their articulatory consequences[END_REF] analysed vocal tract geometries on the basis of x-ray data for 10 Japanese and 10 Caucasian American English speakers. They found a reciprocal relationship between pharyngeal distance (horizontal structure) and lower facial height (vertical structure). Speakers with longer horizontal structures tend to have smaller vertical structures and vice versa. Honda et al. also found some dependencies of vowel production on vocal tract morphology, but not on consonant production.

As a starting point we hypothesise that the shape of the vocal tract partly determines the organization of the articulatory vowel space. More specifically, speakers with a long pharynx (large vertical dimension) and a short palate (short horizontal dimension) use a larger articulatory distance between high and low vowels and a smaller articulatory distance between front and back vowels than speakers with a short pharynx. This hypothesis is based on the assumption that speakers primarily aim towards auditory goals (although the representation of speech is in general multi-modal, see [START_REF] Perrier | Control and representations in speech production[END_REF]), i.e. they adapt their articulation to the respective vocal tract boundary in order to reach the intended goal.

Methods

Nine French speakers were recorded by means of Magnetic Resonance Imaging (MRI) with a Philips Gyroscan T10-NT Powertrack 1000 scanner generating a static longitudinal magnetic field of 1.0 Tesla. An anterior neck coil was used. The repetition time was 1660 ms and the echo delay time was 9 ms. The image matrix was composed of 256x256 pixels, each with a spatial resolution of 1mm in the y-direction and 1.4mm in the x-direction. Data were originally collected for 10 isolated vowels /i e ɛ a y ø oe u o ɔ/ to study inter-speaker acoustic and articulatory variability [START_REF] Apostol | Étude et simulation des caractéristiques individuelles des locuteurs par modélisation du processus de production de la parole[END_REF]. For each vowel three 18 slice series of 3.6mm thick parallel sections with a 0.4mm distance between slices were collected. The acquisition lasted 48s for each sound. The air way was segmented from the surrounding tissues in a manual procedure. We used the itk-SNAP (version 1.4) software for segmentation. So far the corner vowels /i,a,u/ and schwa have been analysed in the 3D space.

Quantifying differences in vocal tract shapes

Quantifying the differences between individual vocal tract geometries turned out to be more complicated than we expected. The literature is full of parameters describing vocal tract growing curves from newborns to adults. However, it is sparse with respect to a parameterisation of inter-speaker differences during adulthood, going beyond the general statement that males have a longer pharynx than females.

We decided to run a Principal Component Analysis (hereafter PCA) on the dataset corresponding to the mid-sagittal outer vocal tract contour. Only the data for the neutral vowel schwa were included in order to avoid effects known for vocal tract length variation in different vowel productions. Using absolute coordinates in a PCA requires a joint spatial reference system of the speakers vocal tracts. The most consistent feature for all speakers' vocal tracts was the posterior pharyngeal wall, which could be approximated by a straight line. The origin of the new reference coordinate system was the lowest identifiable point of the pharyngeal wall. A second point was selected which occurred as high as possible along the pharyngeal wall. These two points were necessary for the translation and rotation of the 2D vocal tract contours so that the posterior pharyngeal wall always corresponded to the y-axis.

Hereafter a PCA was performed. The first factor explained 64 percent of the variance in the data and accounted for differences in pharyngeal length (see Fig. 1, left). The second factor explained 29 percent of the variance in the data. It is less clear how to explain it, but we interpret it with respect to differences in palatal length and height (Fig. 1, right).

Quantifying the articulatory vowel space

To quantify the articulatory vowel space between the corner vowels, we adapted and modified a method proposed by [START_REF] Whalen | Predicting midsagittal pharynx shape from tongue position during vowel production[END_REF] who defined flesh points (hereafter coils) on tongue surface configurations of MRI data (see Fig. 2).

The first coil (coil 1) was placed exactly at the tongue tip and all of the following 5 coils (coil 2-6) were placed in equal distances of 1cm (large dots in Fig. 3). Coils 4 to 6 correspond approximately to the tongue dorsum and tongue back sensors in EMA recordings. They are the most relevant in vowel production. An example of the coils' placement for the three corner vowels is given in Fig. 2.

Next, we drew triangles (corresponding to the vowel space of the corner vowels) for each coil respectively (an example is given in Fig. 3 for coil 3, black triangle in bold). Euclidian distances were calculated for the /ai/, /au/ and /ui/ distances of each coil. In order to normalise for the global differences in vocal tract size while taking into account the relation between high versus low and front versus back vowel distances, we additionally calculated the ratio (hereafter ratio) between the /au/ and /ui/ vowel distances. This ratio was in general >1, since distances between high and low vowels were larger than between front and back vowels. The higher the ratio value the larger the high-low distance relative to the front-back distance.

Results: Relation between geometrical vocal tract parameters and articulatory vowel space

The following table lists only significant results for the correlations between vocal tract parameters (factor f1 and f2) and the articulatory measures (distances: ai, au, ui; ratio). The number of samples for all correlation is always 9 (= 9 speakers).

Factor 1 showed a significant correlation with the /ai/-distance, i.e. the longer the pharynx, the larger the Euclidian distance for high versus low vowels. A similar finding could not be found for the /au/ distance, which may on the one hand be explained with respect to the potential motor equivalence strategies occurring in /u/ (lip protrusion, tongue retraction or laryngeal lowering may equally well contribute to the low second formant values of /u/) but on the other hand it may also be related to the overall shape of the velo-pharyngeal part (the bend from the hard palate to the pharynx). Factor 2 shows a correlation with the /au/-distance, i.e. the shorter the palate, the longer the distance between the low back and the back high vowels. The correlation of the factors with the ratio of the articulatory vowel distances shows in general a good agreement. Significant correlations of f1 with the ratio are found for the more anterior coils at the tongue (2, 3). Factor 2 correlates with the ratio for coils 4 to 6, the ones which are important for vowel production.

Building speaker specific biomechanical tongue models

In a next step we implemented the individual vocal tract configurations and the individual tongue surface contours for the neutral vowel (corresponding to the rest position in the model) in the 2D biomechanical tongue model described in [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF] to build individual models. The implementation involved the following procedures:

1. A rough division of the individual outer vocal tract contour into different anatomical sections in the midsagittal plane (e.g. pharynx, velum, palate etc.) 2. A definition of the beginning and end of the tongue surface contour for the neutral vowel con-figuration and an adaptation of the associated finite element mesh 3. A calculation of the λ commands at rest (according to [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF] recruitment thresholds λ; values equal the length of the muscle fibers at rest, no force generation.) Forces are generated as soon as the tongue shape changes slightly. Four examples are shown in Figure 4. These individual models will allow us to discover the nonlinear relations between motor commands, kinematic output, and acoustics, and to account for the potential effects emerging from the supine position of the subject during the MRI recording. From the 3D MRI volume data during schwa production we also computed the acoustics (first four formants) in order to check the actual realization of the vowel corresponding to the rest position in the model.

So far we run simulations for the model cb (short pharynx) and gg (long pharynx). To do so, we used the λ commands for the rest position and the corner vowels /a, i, u/ provided in [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF] as a reference and calculated the differences between λ at rest and λ at one of the corner vowels. The differences in λ rather than absolute λ commands were applied to the cb and gg model. Starting from the rest position 50 simulations were run for each of the corner vowels. λ commands were randomly chosen among all the possible combinations in the specified range of the target command.

Figure 5 displays the results of our simulations showing /a/ variations around the target position. Tongue variability for /a/ looks similar for both speakers, but for cb the pharynx gets nearly closed due to the downward and backward movement of the tongue whereas for gg it gets more constricted. This is in agreement with our hypothesis that speakers with a short pharynx and a long palate have to control their vertical movements precisely since their vocal tract constraints them to move the tongue in the vertical direction. However, there is probably a tradeoff between tongue and jaw movements. With an open jaw the tongue can move downwards without a pharyngeal constriction. To what extent more or less vigorous jaw movement may be related to the individual vocal tract geometry seems to be an interesting question for future studies.

Figure 1 .

 1 Figure 1. Mean outer vocal tract contour (red line) varied solely by the full range of values for factor 1 (left) and solely by the full range of values for factor 2 (right).

Figure 2 .

 2 Figure 2. Schematic view of flesh point markers (coils) defined on the tongue contour in the mid-sagittal plane (schwa, speaker fs).

Figure 3 .

 3 Figure 3. Example of coil placement: green/highest tongue contour = /i/ configuration, blue/middle tongue contour = /u/ configuration, red/lowest tongue contour = /a/ configuration, in black mid-outer vocal tract contour

Figure 4 .

 4 Figure 4. Individual 2D biomechanical tongue models based on the MRI data (gg, cs = long pharynx, cb, av = short pharynx)

Figure 5 .

 5 Figure 5. Simulations of /a/ variability using a similar set of λ variation: left for cb and right for gg

Table 1 .

 1 Significant correlations between factor scores (f1 and f2) and articulatory distances

	Significant	Coil	R	P-value
	correlation
	f1 -ai	5	-0.716	0.027
		6	-0.719	0.026
	f2 -au	6	-0.714	0.028
	f1 -ratio	2	-0.703	0.032
		3	-0.684	0.041
	f2 -ratio	4	-0.837	0.003
		5	-0.846	0.002
		6	-0.803	0.007

Acknowledgements

We like to thank Lian Apostol for providing the MRI data and Mark Tiede for the idea running a PCA on the vocal tract data in order to parameterize the differences. This work was primarily supported by the Christian Benoit Association by means of an award dedicated to the first author. It is dedicated to Christian Benoit and Dieter Fuchs.


Eric Benhamou
email: eric.benhamou@pricingpartners.com

Emmanuel Gobet
email: emmanuel.gobet@imag.fr

Mohammed Miri
email: mohammed.miri@pricingpartners.com

TIME DEPENDENT HESTON MODEL

Keywords: asymptotic expansion, Malliavin calculus, small volatility of volatility, time dependent Heston model AMS subject classifications. 60J75, 60HXX

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction.

Stochastic volatility modeling has emerged in the late nineties as a way to manage the smile. In this work, we focus on the Heston model which is a lognormal model where the square of volatility follows a CIR 1 process. The call (and put) price has a closed formula in this model thanks to a Fourier inversion of the characteristic function (see Heston [START_REF] Heston | A closed-form solutions for options with stochastic volatility with applications to bond and currency options[END_REF], Lewis [START_REF] Lewis | Option valuation under stochastic volatility with Mathematica code[END_REF] and Lipton [START_REF] Lipton | The vol smile problem[END_REF]). When the parameters are piecewise constant, one can still derive a recursive closed formula using a PDE method (see Mikhailov and Nogel [START_REF] Mikhailov | Hestons stochastic volatility model implementation, calibration and some extensions[END_REF]) or a Markov argument in combination with affine models (see Elices [START_REF] Elices | Models with time-dependent parameters using transform methods: application to Heston's model[END_REF]), but formula evaluation becomes increasingly time consuming. However, for general time dependent parameters there is no analytical formula and one usually has to perform Monte Carlo simulations. This explains the interest of recent works for designing more efficient Monte Carlo simulations: see Broadie and Kaya [START_REF] Broadie | Exact simulation of stochastic volatility and other affine jump diffusion processes[END_REF] for an exact simulation and bias-free scheme based on Fourier integral inversion; see Andersen [START_REF] Andersen | Efficient simulation of the Heston stochastic volatility model[END_REF] based on a Gaussian moment matching method and a user friendly algorithm; see Smith [START_REF] Smith | An almost exact simulation method for the Heston model[END_REF] relying on an almost exact scheme; see Alfonsi [START_REF] Alfonsi | High order discretization schemes for the CIR process: application to Affine Term Structure and Heston models[END_REF] using higher order schemes and a recursive method for the CIR process. For numerical partial differential equations, we refer the reader to Kluge's doctoral dissertation [START_REF] Kluge | Pricing derivatives in stochastic volatility models using the Finite Difference method[END_REF].

Comparison with the literature. A more recent trend in the quantitative literature has been the use of the so called approximation method to derive analytical formulae. This has led to an impressive number of papers, with many original ideas. For instance, Al òs et al. [START_REF] Al Òs | On the short-time behaviour of the implied volatility for jump-diffusion models with stochastic volatility[END_REF] have been studying the short time behavior of implied volatility for stochastic volatility using an extension of It ô's formula. Another trend has focused on analytical techniques to derive the asymptotic expansion of the implied volatility near expiry (see for instance Berestycki et al. [START_REF] Berestycki | Computing the implied volatility in stochastic volatility models[END_REF], [START_REF] Henry Labordère | A general asymptotic implied volatility for stochastic volatility models[END_REF], Hagan et al. [START_REF] Hagan | Managing smile risk[END_REF], Lewis [START_REF] Lewis | Geometries and smile asymptotics for a class of stochastic volatility models[END_REF], Osajima [START_REF] Osajima | The asymptotic expansion formula of implied volatility for dynamic SABR model and FX hybrid model[END_REF] or Forde [START_REF] Forde | The small-time behavior of diffusion and time-changed diffusion processes on the line[END_REF]). But in these works the implied volatility near expiry does not have a closed formula because the related geodesic distance is not explicit. It can, however, be approximated by a series expansion [START_REF] Lewis | Geometries and smile asymptotics for a class of stochastic volatility models[END_REF]. The drawback to these methods is their inability to handle non-homogeneous (that is to say time dependent) parameters. For long maturities, another approach has been the asymptotic expansion w.r.t. the mean reversion parameter of the volatility as shown in [START_REF] Fouque | Derivatives in financial Markets with stochastic volatility[END_REF]. In the case of zero correlation, averaging techniques as exposed in [START_REF] Piterbarg | Time to smile[END_REF] and [START_REF] Piterbarg | Stochastic volatility model with time-dependent skew[END_REF] can be used. Antonelli and Scarlatti take another view in [START_REF] Antonelli | Pricing options under stochastic volatility: a power series approach[END_REF] and have suggested price expansion w.r.t. correlation. For all of these techniques, the domain of availability of the expansion is restricted to either short or long maturities, to zero correlation, or to homogeneous parameters. In our work, we aim to give an analytical formula which covers both short and long maturities, that also handles time inhomogeneous parameters as well as non-null correlations. As a difference with several previously quoted papers, our purpose consists also of justifying mathematically our approximation.

The results closest to ours are probably those based on an expansion w.r.t. the volatility of volatility by Lewis [Lew00]: it is based on formal analytical arguments and is restricted to constant parameters. Our formula can be viewed as an extension of Lewis' formula in order to address a time dependent Heston model, using a direct probabilistic approach. In addition, we prove an error estimate which shows that our approximation formula for call/put is of order 2 w.r.t. the volatility of volatility. The advantage of this current approximation is that the evaluation is about 100 to 1000 times quicker than a Fourier based method (see our numerical tests). tial equation (SDE):

dX t = √ v t dW t - v t 2 dt, X 0 = x 0 , (1.1) dv t = κ(θ t -v t )dt + ξ t √ v t dB t , v 0 , (1.2) d W, B t = ρ t dt,
where (B t , W t ) 0≤t≤T is a two-dimensional correlated Brownian motion on a filtered probability space (Ω, F , (F t ) 0≤t≤T , P) with the usual assumptions on filtration (F t ) 0≤t≤T .

In our setting, (X t ) t is the log of the forward price and (v t ) t is the square of the volatility which follows a CIR process with an initial value v 0 > 0, a positive mean reversion κ, a positive long-term level (θ t ) t , a positive volatility of volatility (ξ t ) t and a correlation (ρ t ) t . These time dependent parameters are assumed to be measurable and bounded on [0, T].

To develop our approximation method, we will examine the following perturbed process w.r.t. ǫ ∈ [0, 1]:

dX ǫ t = v ǫ t dW t - v ǫ t 2 dt, X ǫ 0 = x 0 , dv ǫ t = κ(θ t -v ǫ t )dt + ǫξ t v ǫ t dB t , v ǫ 0 = v 0 , (1.3)
so that our perturbed process coincides with the initial one for ǫ = 1 : X 1 t = X t , v 1 t = v t . For the existence of the solution v ǫ , we refer to Chapter IX in [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF] (moreover, the process is non-negative for kθ t ≥ 0, see also the proof of Lemma 4.2). Our main purpose is to give an accurate analytic approximation, in a certain sense, of the expected payoff of a put option : g(ǫ) = e -T 0 r t dt E[(Ke T 0 (r t -q t )dt+X ǫ

T ) + ] (1.4)
where r (resp. q) is the risk-free rate (resp. the dividend yield), T is the maturity and ǫ = 1. Extensions to call options and other payoffs are discussed later.

Outline of the paper. In Section 2, we explain the methodology of the small volatility of volatility expansion. An approximation formula is then derived in Theorem 2.3 and its accuracy stated in Theorem 2.4. This section ends by explicitly expressing the formula's coefficients for general time dependent parameters (constant, smooth and piecewise constant). Our expansion allows us to give equivalent constant parameters for the time dependent Heston model (see Subsection 2.6). As a second corollary, the options calibration for Heston's model using only one maturity becomes an ill-posed problem; we give numerical results to confirm this situation. In section 3, we provide numerical tests to benchmark our formula with the closed formula in the case of constant and piecewise constant parameters. In Section 4, we prove the accuracy of the approximation stated in Theorem 2.4: this section is the technical core of the paper. In Section 5, we establish lemmas used to make the calculation of the correction terms explicit (those derived in Theorem 2.3). In Section 6, we conclude this work and give a few extensions. In the appendix, we recall details about the closed formula (of Heston [START_REF] Heston | A closed-form solutions for options with stochastic volatility with applications to bond and currency options[END_REF] and Lewis [START_REF] Lewis | Option valuation under stochastic volatility with Mathematica code[END_REF]) in the case of constant (and piecewise constant) parameters.

Smart Taylor expansion.

Notations. N Extremes of deterministic functions.

For a càdlàg function l : [0, T] R, we denote l In f = inf t∈[0,T] l t and l Sup = sup t∈[0,T] l t .

N 2.2. Differentiation.

(i) For a smooth function x → l(x), we denote by l (i) (x) its i-th derivative.

(ii) Given a fixed time t and for a function ǫ f ǫ t , we denote (if it has a meaning) the i th derivative at ǫ = 0 by f i,t =

∂ i f ǫ t ∂ǫ i | ǫ=0 .
2.2. Definitions. In order to make the approximation explicit, we introduce the following family of operators indexed by maturity T.

D 2.1. Integral Operator. We define the integral operator ω (.,.)

.,T as follows: (i) For any real number k and any integrable function l, we set

ω (k,l) t,T = T t e ku l u du, ∀t ∈ [0, T]. (ii) For any real numbers (k 1 , • • • , k n ) and for any integrable functions (l 1 , • • • , l n ), the n-times iteration is given by ω (k 1 ,l 1 ),••• ,(k n ,l n ) t,T = ω (k 1 ,l 1 ω (k 2 ,l 2 ),••• ,(kn,ln) .,T ) t,T , ∀t ∈ [0, T].
(iii) When the functions (l 1 , • • • , l n ) are equal to the unity constant function 1, we simply write

ωk 1 ,••• ,k n t,T = ω (k 1 ,1),••• ,(k n ,1) t,T , ∀t ∈ [0, T].

About the CIR process. Assumptions.

In order to bound the approximation errors, we need a positivity assumption for the CIR process.

Assumption (P). The parameters of the CIR process (1.2) verify the following conditions:

ξ In f > 0, ( 2κθ ξ 2 ) In f ≥ 1.
This assumption is crucial to ensure the positivity of the process on [0, T], which is stated in detail in Lemma 4.2 (remember that v 0 > 0). We have

P(∀t ∈ [0, T] : v t > 0) = 1.
When the functions θ and ξ are constant, Assumption (P) coincides with the usual Feller test condition 2κθ ξ 2 ≥ 1 (see [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]). Note that the above assumption ensures that the positivity property also holds for the perturbed CIR process (1.3): for any ǫ ∈ [0, 1], we have

P(∀t ∈ [0, T] : v ǫ t > 0) = 1 (see Lemma 4.2).
We also need a uniform bound of the correlation in order to preserve the non degeneracy of the SDE (1.1) conditionally on (B t ) 0≤t≤T . Assumption (R). The correlation is bounded away from -1 and +1:

|ρ| Sup < 1.

Taylor Development.

In this paragraph, we present the main steps leading to our results. Complete proofs are given later.

If (F B t ) t denotes the filtration generated by the Brownian motion B, the distribution of X ǫ T conditionally to F B T is a Gaussian distribution with mean

x 0 + T 0 ρ t v ǫ t dB t - 1 2 T 0 v ǫ t dt and variance T 0 (1 -ρ 2 t )v ǫ t dt (ǫ ∈ [0, 1]
). Therefore, the function (1.4) can be expressed as follows:

g(ǫ) = E[P BS (x 0 + T 0 ρ t v ǫ t dB t - T 0 ρ 2 t 2 v ǫ t dt, T 0 (1 -ρ 2 t )v ǫ t dt)], (2.1)
where the function (x, y) P BS (x, y) is the put function price in a Black-Scholes model with spot e x , strike K, total variance y, risk-free rate r eq = T 0 r(t)dt T

, dividend yield

q eq = T 0 q(t)dt T
and maturity T. For the sake of completeness, we recall that P BS (x, y) has the following explicit expression Ke -r eq T N 1 √ y log( Ke -r eq T e x e -q eq T ) +

1 2 √ ye x e -q eq T N 1 √ y log( Ke -r eq T e x e -q eq T ) -

1 2 √ y .
In the following, we expand P BS (., .) with respect to its two arguments. For this, we note that P BS is a smooth function (for y > 0). In addition, there is a simple relation between its partial derivatives:

∂P BS ∂y (x, y) = 1 2 ( ∂ 2 P BS ∂x 2 (x, y) - ∂P BS ∂x (x, y)), ∀x ∈ R, ∀y > 0, (2.2)
which can be proved easily by a standard calculation left to the reader.

Under assumption (P), for any t, v ǫ t is C 2 w.r.t ǫ at ǫ = 0 (differentiation in L p -sense). This result will be shown later. In addition, v ǫ does not vanish (for any ǫ ∈ [0, 1]).

Hence, by putting v

ǫ i,t = ∂ i v ǫ t ∂ǫ i , we get dv ǫ 1,t = -κv ǫ 1,t dt + ξ t v ǫ t dB t + ǫξ t v ǫ 1,t 2 v ǫ t dB t , v ǫ 1,0 = 0, dv ǫ 2,t = -κv ǫ 2,t dt + ξ t v ǫ 1,t v ǫ t dB t + ǫξ t v ǫ 2,t 2 v ǫ t dB t -ǫξ t [v ǫ 1,t ] 2 4[v ǫ t ] 3/2 dB t , v ǫ 2,0 = 0.
From the definitions v i,t ≡

∂ i v ǫ t ∂ǫ i | ǫ=0 , we easily deduce v 0,t = e -κt (v 0 + t 0 κe κs θ s ds), v 1,t = e -κt t 0 e κs ξ s √ v 0,s dB s , (2.3) v 2,t = e -κt t 0 e κs ξ s v 1,s (v 0,s ) 1 2 dB s .
(2.4)

Note that v 0,t coincides also with the expected variance E(v t ) because of the linearity of the drift coefficient of (v t ) t . Now, to expand g(ǫ), we use the Taylor formula twice, first applied to ǫ v ǫ t and v ǫ t at ǫ = 1 using derivatives computed at ǫ = 0:

v 1 t = v 0,t + v 1,t + v 2,t 2 + • • • , v 1 t = √ v 0,t + v 1,t 2(v 0,t ) 1 2 + v 2,t 4(v 0,t ) 1 2 - v 2 1,t 8(v 0,t ) 3 2 + • • • ,
secondly for the smooth function P BS at the second order w.r.t. the first and second variable around (x 0

+ T 0 ρ t √ v 0,t dB t - T 0 ρ 2 t 2 v 0,t dt, T 0 (1 -ρ 2 t )v 0,t dt).
For convenience, we simply write

PBS =P BS (x 0 + T 0 ρ t √ v 0,t dB t - T 0 ρ 2 t 2 v 0,t dt, T 0 (1 -ρ 2 t )v 0,t dt), (2.5) ∂ i+j PBS ∂x i y j = ∂ i+j P BS ∂x i y j (x 0 + T 0 ρ t √ v 0,t dB t - T 0 ρ 2 t 2 v 0,t dt, T 0 (1 -ρ 2 t )v 0,t dt).
Then, one gets

g(1) =E[ PBS ] (2.6) +E[ ∂ PBS ∂x ( T 0 ρ t ( v 1,t 2(v 0,t ) 1 2 + v 2,t 4(v 0,t ) 1 2 - v 2 1,t 8(v 0,t ) 3 2 )dB t - T 0 ρ 2 t 2 (v 1,t + v 2,t 2 )dt)] (2.7) +E[ ∂ PBS ∂y T 0 (1 -ρ 2 t )(v 1,t + v 2,t 2 )dt] (2.8) + 1 2 E[ ∂ 2 PBS ∂x 2 ( T 0 ρ t v 1,t 2(v 0,t ) 1 2 dB t - T 0 ρ 2 t 2 v 1,t dt) 2 ]
(2.9)

+ 1 2 E[ ∂ 2 PBS ∂y 2 ( T 0 (1 -ρ 2 t )v 1,t dt) 2 ] (2.10) +E[ ∂ 2 PBS ∂xy ( T 0 (1 -ρ 2 t )v 1,t dt)( T 0 ρ t v 1,t 2(v 0,t ) 1 2 dB t - T 0 ρ 2 t 2 v 1,t dt)] (2.11) +E (2.12)
where E is the error in our Taylor expansion. In fact, we notice that:

E[ PBS ] = E[E[e -T 0 r t dt (K -e x 0 + T 0 (r t -q t - v 0,t 2 )dt+ T 0 √ v 0,t (ρ t dB t + 1-ρ 2 t dB ⊥ t ) ) + |F B T ]] = P BS (x 0 , T 0 v 0,t dt),
where B ⊥ is a Brownian motion independent on F B T . Furthermore, the relation (2.2) remains the same for PBS and this enables us to simplify the expansion above. This gives:

P 2.2. The approximation (2.12) is equivalent to

g(1) = P BS (x 0 , T 0 v 0,t dt) + E[ ∂ PBS ∂y T 0 (v 1,t + v 2,t )dt] + 1 2 E[ ∂ 2 PBS ∂y 2 ( T 0 v 1,t dt) 2 ] + E.
The details of the proof are given in Subsection 5.2. At first sight, the above formula looks like a Taylor formula of P BS w.r.t. the cumulated variance. In fact, it is different, note that the coefficient of v 2,t is not 1/2 but 1. We do not have any direct interpretation of this formula. The next step consists of making explicit the correction terms as a combination of Greeks of the BS price.

T 2.3. Under assumptions (P) and (R), the put2 price is approximated by

e -T 0 r t dt E[(K -e T 0 (r t -q t )dt+X 1 T ) + ] =P BS (x 0 , var T ) + 2 i=1 a i,T ∂ i+1 P BS ∂x i y (x 0 , var T ) + 1 i=0 b 2i,T ∂ 2i+2 P BS ∂x 2i y 2 (x 0 , var T ) + E, (2.13)
where

var T = T 0 v 0,t dt, a 1,T = ω (κ,ρξv 0,. ),(-κ,1) 0,T , a 2,T = ω (κ,ρξv 0,. ),(0,ρξ),(-κ,1) 0,T , b 0,T = ω (2κ,ξ 2 v 0,. ),(-κ,1),(-κ,1) 0,T , b 2,T = a 2 1,T 2 .
The proof is postponed to Subsection 5.3. Finally, we give an estimate regarding the error E arising in the above theorem. T 2.4. Under assumptions (P) and (R), the error in the approximation (2.13) is estimated as follows:

E = O [ξ Sup √ T] 3 √ T .
In view of Theorem 2.4, we may refer to the formula (2.13) as a second order approximation formula w.r.t. the volatility of volatility.

Computation of coefficients.

Constant parameters. The case of constant parameters (θ, ξ, ρ) gives us the coefficients a and b explicitly. Indeed in this case, the operator ω is a simple iterated integration of exponential functions. Using Mathematica, we derive the following explicit expressions.

P 2.5. Explicit computations. For constant parameters, one has:

var T =m 0 v 0 + m 1 θ, a 1,T =ρξ(p 0 v 0 + p 1 θ), a 2,T =(ρξ) 2 (q 0 v 0 + q 1 θ), b 0,T =ξ 2 (r 0 v 0 + r 1 θ).
where

m 0 = e -κT -1 + e κT κ , m 1 = T - e -κT -1 + e κT κ , p 0 = e -κT -κT + e κT -1 κ 2 , p 1 = e -κT κT + e κT (κT -2) + 2 κ 2 , q 0 = e -κT -κT(κT + 2) + 2e κT -2 2κ 3
, q 1 = e -κT 2e κT (κT -3) + κT(κT + 4) + 6 2κ 3 ,

r 0 = e -2κT -4e κT κT + 2e 2κT -2 4κ 3 , r 1 = e -2κT 4e κT (κT + 1) + e 2κT (2κT -5) + 1 4κ 3 .
R 2.1. In the case of constant parameters (θ, ξ, ρ), we retrieve the usual Heston model. In this particular case, our expansion coincides exactly with Lewis' volatility of volatility series expansion (see Equation (3.4), page 84 in [START_REF] Lewis | Option valuation under stochastic volatility with Mathematica code[END_REF] for Lewis' expansion formula and page 93 in [START_REF] Lewis | Option valuation under stochastic volatility with Mathematica code[END_REF] for the explicit calculation of the coefficients J (i) with ϕ = 1 2 ). Using his notation, we have a 1,T = J (1) , a 2,T = J (4) and b 0,T = J (3) .

Smooth parameters. In this case, we may use a Gauss-Legendre quadrature formula for the computation of the terms a and b.

Piecewise constant parameters. The computation of the variance var T is straightforward. Thus, it remains to provide explicit expressions of a and b as a function of the piecewise constant data. Let

T 0 = 0 ≤ T 1 ≤ • • • ≤ T n = T such that θ, ρ, ξ are constant on each interval ]T i , T i+1
[ and are equal respectively to θ T i+1 , ρ T i+1 , ξ T i+1 . Before giving the recursive relation, we need to introduce the following functions: ω1,t = ω (κ,ρξv 0,. ) 0,t , ω2,t = ω (2κ,ξ 2 v 0,. ) 0,t , α t = ω (k,ρξv 0,. ),(0,ρξ) 0,t , β t = ω (2k,ξ 2 v 0,. ),(-κ,1) 0,t . P 2.6. Recursive calculations. For piecewise constant coefficients, one has:

a 1,T i+1 = a 1,T i + ω-κ T i ,T i+1 ω1,T i + ρ T i+1 ξ T i+1 f 1 κ,v 0,T i (θ T i+1 , T i , T i+1 ), a 2,T i+1 = a 2,T i + ω-κ T i ,T i+1 α T i + ρ T i+1 ξ T i+1 ω0,-κ T i ,T i+1 ω1,T i + (ρ T i+1 ξ T i+1 ) 2 f 2 κ,v 0,T i (θ T i+1 , T i , T i+1 ), b 0,T i+1 = b 0,T i + ω-κ T i ,T i+1 β T i + ω-κ,-κ T i ,T i+1 ω2,T i + ξ 2 T i+1 f 0 κ,v 0,T i (θ T i+1 , T i , T i+1 ), α T i+1 = α T i + ρ T i+1 ξ T i+1 (T i+1 -T i ) ω1,T i + ρ 2 T i+1 ξ 2 T i+1 g 1 κ,v 0,T i (θ T i+1 , T i , T i+1 ), β T i+1 = β T i + ω-κ T i ,T i+1 ω2,T i + ξ 2 T i+1 g 2 κ,v 0,T i (θ T i+1 , T i , T i+1 ), ω1,T i+1 = ω1,T i + ρ T i+1 ξ T i+1 h 1 κ,v 0,T i (θ T i+1 , T i , T i+1 ), ω2,T i+1 = ω2,T i + ξ 2 T i+1 h 2 κ,v 0,T i (θ T i+1 , T i , T i+1 ), v 0,T i+1 = e -κ(T i+1 -T i ) (v 0,T i -θ T i+1 ) + θ T i+1 ,
where f 0 κ,v 0 (θ, t, T) = e -2κT (e 2κt (θ-2v 0 )+e 2κT ((-2κt+2κT-5)θ+2v 0 )+4e κ(t+T) ((-κt+κT+1)θ+κ(t-T)v 0 )) 4κ 3

, f 1 κ,v 0 (θ, t, T) = e -κT (e κT ((-κt+κT-2)θ+v 0 )-e κt ((κt-κT-2)θ-κtv 0 +κTv 0 +v 0 ))

κ 2 , f 2 κ,v 0 (θ, t, T) = e -κ(t+3T) (2e κ(t+3T) ((κ(T-t)-3)θ+v 0 )+e 2κ(t+T) ((κ(κ(t-T)-4)(t-T)+6)θ-(κ(κ(t-T)-2)(t-T)+2)v 0 )) 2κ 3 , g 1 κ,v 0 (θ, t, T) = 2e κT θ+e κt (κ 2 (t-T) 2 v 0 -(κ(κ(t-T)-2)(t-T)+2)θ) 2κ 2 , g 2 κ,v 0 (θ, t, T) = e -κT (e 2κT θ-e 2κt (θ-2v 0 )+2e κ(t+T) (κ(t-T)(θ-v 0 )-v 0 )) 2κ 2 , h 1 κ,v 0 (θ,
a 1,T i+1 = T i 0 e κt ρ t ξ t v 0,t ω (-κ,1) t,T i+1 dt + T i+1 T i e κt ρ t ξ t v 0,t ω (-κ,1) t,T i+1 dt =a 1,T i + T i 0 e κt ρ t ξ t v 0,t ω (-κ,1) T i ,T i+1 dt + T i+1 T i e κt ρ t ξ t v 0,t ω (-κ,1) t,T i+1 dt =a 1,T i + ω (-κ,1) T i ,T i+1 T i 0 e κt ρ t ξ t v 0,t dt + T i+1 T i e κt ρ t ξ t v 0,t ω (-κ,1) t,T i+1 dt =a 1,T i + ω-κ T i ,T i+1 ω1,T i + ρ T i+1 ξ T i+1 f 1 κ,v 0,T i (θ i+1 , T i , T i+1 ),
where the functions f 1 κ,v 0 and ω-κ are calculated analytically using Mathematica. The other terms are calculated analogously.

Corollaries of the approximation formula (2.13).

Averaging Heston's model parameters. We derive a first corollary of the approximation formula in terms of equivalent Heston models. As explained in [START_REF] Piterbarg | Time to smile[END_REF], this averaging principle may facilitate efficient calibration. Namely, we search for equivalent constant parameters κ, θ, ξ, ρ for the Heston model3

d Xt = √ vt dW t - vt 2 dt, X0 = x 0 , d vt = κ( θt -vt )dt + ξ √ vt dB t , v0 = v 0 , d W, B t = ρdt,
that equalize the price of call/put options maturing at T in the time dependent model (equality up to the approximation error E). The following rules give the equivalent parameters as a function of the variance var T and the coefficients a 1,T , a 2,T , b 0,T that are computed in the time dependent model. Results are expressed using

a = a 2,T m 1 m 1 q 0 -m 0 q 1 , b = - a 1,T m 1 m 1 p 0 -m 0 p 1 , c = var T ( p 1 m 1 p 0 -m 0 p 1 - q 1 m 1 q 0 -m 0 q 1 ),
where m 0 , m 1 , p 0 , p 1 , q 0 , q 1 , r 0 and r 1 are given in Proposition 2.5.

Averaging rule in the case of zero correlation. If ρ t ≡ 0, the equivalent constant parameters (for maturity T) are

κ = κ, θ = var T -m 0 v 0 m 1 , ξ = b 0,T r 0 v 0 + r 1 θ , ρ = 0.
Proof. Two sets of prices coincide at maturity T if they have the same approximation formula (2.13). In this case a 1,T = a 2,T = b 2,T = 0, thus the approximation formula depends only on two quantities var T and b 0,T . It is quite clear that there is not a single choice of parameters to fit these two quantities. A simple solution results from the choice of κ = κ and ρ = 0: then, using Proposition 2.5, we obtain the announced parameters θ and ξ.

R 2.2. In this case of zero correlation and θ = v 0 = θ, we exactly retrieve Piterbarg's results for the averaged volatility of volatility ξ ( see [START_REF] Piterbarg | Time to smile[END_REF]).

Averaging rule in the case of non zero correlation. We follow the same arguments as before. Now the approximation formula also depends on the four quantities var T , a 1,T , a 2,T and b 2,T . Thus, equalizing call/put prices at maturity T is equivalent to equalizing these four quantities in both models, by adjusting κ, θ, ξ and ρ. Unfortunately, we have not found a closed expression for these equivalent parameters. An alternative and simpler way of proceeding consists of modifying the unobserved initial value v0 of the variance process while keeping κ = κ. For non-vanishing correlation (ρ t ) t , it leads to two possibilities

v0 =b (b ± √ b 2 -4ac) 2a - p 1 var T m 1 p 0 -m 0 p 1 , θ = var T -m 0 v0 m 1 , ξ = b 0,T r 0 v0 + r 1 θ , ρ = - 2a ξ(b ± √ b 2 -4ac) .
In practice, only one solution gives realistic parameters. However, this rule is heuristic since there is a priori no guarantee that these averaged parameters satisfy the assumption (P), which is the basis for the arguments correctness.

Proof. Using Proposition 2.5, one has to solve the following system of equations

var T =m 0 v0 + m 1 θ, a 1,T = ρ ξ(p 0 v0 + p 1 θ),
a 2,T =( ρ ξ) 2 (q 0 v0 + q 1 θ), b 0,T = ξ2 (r 0 v0 + r 1 θ).

The first equation gives θ =

var T -m 0 v0 m 1

. Replacing this identity in a 1,T and a 2,T gives

v0 =( a 1,T ( ρ ξ) - p 1 var T m 1 ) m 1 p 0 m 1 -p 1 m 0 , v0 =( a 2,T ( ρ ξ) 2 - q 1 var T m 1 ) m 1 q 0 m 1 -q 1 m 0 .
It readily leads to a quadratic equation ax 2 + bx + c = 0 with x = 1 ρ ξ . By solving this equation, we easily complete the proof of the result.

Collinearity effect in the Heston model. Another corollary of the approximation formula (2.13) is that we can obtain the same vanilla prices at time T with different sets of parameters. For instance, take on the one hand v 0 = θ = 4%, κ 1 = 2 and ξ 1 = 30% (model M 1 ) and on the other hand v 0 = θ = 4%, κ 2 = 3 and ξ 2 = 38.042% (model M 2 ), both models having zero correlation. The resulting error between implied volatilities within the two models are presented in Table 2.1: they are so small that prices can be considered as equal. Actually, this kind of example is easy to create even with non-null correlation: as before, in view of the approximation formula (2.13), it is sufficient to equalize the four quantities var T , a 1,T , a 2,T and b 2,T .

As a consequence, calibrating a Heston model using options with a single maturity is an ill-posed problem, which is not a surprising fact.

Numerical accuracy of the approximation.

We give numerical results of the performance of our method. In what follows, the spot S 0 , the risk-free rate r and the dividend yield q are set respectively to 100, 0% and 0%. The initial value of the variance process is set to v 0 = 4% (initial volatility equal to 20%). Then we study the numerical accuracy w.r.t. K, T, κ, θ, ξ and ρ by testing different values for these parameters. In order to present more interesting results for various relevant maturities and strikes, we allow the range of strikes to vary over the maturities. The strike values evolve approximately as S 0 exp(c √ θT) for some real numbers c and θ = 6%. The extreme values of c are chosen to be equal to ±2.57, which represents the 1%-99% quantile of the standard normal distribution. This corresponds to very out-of-themoney options or very deep-in-the-money options. The set of pairs (maturity, strike) chosen for the tests is given in Table 3.1.

Constant parameters. In Table 3.2, we report the numerical results when θ = 6%, κ = 3, ξ = 30% and ρ = 0%, giving the errors of implied Black-Scholes volatilities between our approximation formula (see Equation (2.13)) and the price calculated using the closed formula (see appendix), for the maturities and strikes of Table 3.1. The table should be read as follows: for example, for one year maturity and strike equal to 170, the implied volatility is equal to 24.14% using the closed formula and 24.20% with the approximation formula, giving an error of -6.33 bps. In Table 3.2, we observe that the errors do not exceed 7 bps for a large range of strikes and maturities. We notice that the errors are surprisingly higher for short maturities. At first sight, it is counterintuitive as one would expect our perturbation method to work better for short maturities and worse for long maturities, since the difference between our proxy model (BS with volatility (v 0,t ) t ) and the original one is increasing w.r.t. time. In fact, this intuition is true for prices but not for implied volatilities. When we T 3.2 Implied Black-Scholes volatilities of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 3, ξ = 30% and ρ = 0%.

3M

23 compare the price errors (in Price bp 4 ) for the same data, we observe in Table 3.3 that the error terms are not any bigger for short maturities but vary slightly over time with two observed effects. The error term first increases over time as the error between the proxy and the original model increases over time, as forecasted. But for long maturities, presumably because the volatility converges to its stationary regime, errors decrease. When we convert these prices to implied Black-Scholes volatilities, these error terms are dramatically amplified for short maturities due to very small vega. Finally, note that for fixed maturity, price errors are quite uniform w.r.t. strike K. Impact of the correlation. Analogous results for correlations equal to -20%, 20% and -50% are reported in Tables 3.4-3.5, 3.6-3.7 and 3.8-3.9. We notice that the errors are smaller for a correlation close to zero and become larger when the absolute value of the correlation increases. However, for realistic correlation values (-50% for instance), the accuracy for the usual maturities and strikes remains excellent (error smaller than 20 bps), except for very extreme strikes.

Impact of the volatility of volatility. In view of Theorem 2.4, the smaller the volatility of volatility, the more accurate the approximation. In the following numerical tests, we increase ξ, while maintaining Assumption (P). Thus, the new Heston's parameters are κ = 10, ξ = 1 and ρ = -50%, the other parameters remaining unchanged. The comparative results on implied volatilties and prices are presented in Table 3.10 and 3.11. As expected, the approximation is less accurate than for ξ = 30%, but still accurate enough to be efficiently used for fast calibration. The results for prices are 4 Error price bp= more satisfactory than for implied volatilities. Once again, for short maturities, the errors in implied volatilities may be quite significant, except for options not-far-fromthe-money. Impact of the assumption (P). The assumption (P) is a technical assumption that we use to establish error estimates for the approximation formula (2.13). In the test that follows, we relax this assumption by taking new parameters θ = 3%, κ = 2, ξ = 40% and ρ = 0% for which the ratio 2κθ/ξ 2 = 0.75 < 1. Results are reported in Tables 3.12 and 3.13. We observe that the approximation formula still works (errors are smaller than 20 bps) but it is less accurate (compare with Table 3.2 for which the ratio 2κθ/ξ 2 is equal to 4 > 1). An extension of the validity of our formula by relaxing Assumption (P) is presumably relevant. This investigation is left for further research.

Piecewise constant parameters. Heston's constant parameters have been set to: v 0 = 4%, κ = 3. In addition, the piecewise constant functions θ, ξ and ρ are equal respectively at each interval of the form ] i 4 , i+1 4 [ to 4% + i × 0.05%, 30% + i × 0.5% and -20% + i × 0.35%. In the same Tables 3.15 and 3.16, we report values using three different formulas. For a given maturity, the first row is obtained using the closed formula with piecewise constant parameters (see appendix), the second row uses our approximation formula (2.13) and the third row uses the closed formula with constant parameters computed by averaging (see Section 2.6). In order to give complete information on our tests, we also report in Table 3.14 the values used for the averaging parameters (following Section 2.6).

Of course, the quickest approach is the use of the approximation formula (2.13). As before, its accuracy is very good, except for very extreme strikes. It is quite interesting is to observe that the averaging rules that we propose are extremely T 3.4 Implied Black-Scholes volatilities of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 3, ξ = 30% and ρ = -20%.

3M

24 accurate.

Computational time. Regarding the computational time, the approximation formula (2.13) yields essentially the same computational cost as the Black-Scholes formula, while the closed formula requires an additional space integration involving many exponential and trigonometric functions for which evaluation costs are higher. For instance, using a 2, 6 GHz Pentium PC, the computations of the 64 numerical values in Table 3.2 (3.4, 3.6 or 3.8) take 4.71 ms using the approximation formula and 301ms using the closed formula. For the example with time dependent coefficients (reported in Table 3.15), the computational time for the 64 prices is about 40.2 ms using the approximation formula and 2574 ms using the closed formula. Roughly speaking, the use of the approximation formula enables us to speed up the valuation (and thus the calibration) by a factor 100 to 600.

Proof of Theorem 2.4.

The proof is divided into several steps. In Subsection 4.1 we give the upper bounds for derivatives of the put function P BS , in Subsection 4.2 the conditions for positivity of the squared volatility process v, in Subsection 4.3 the upper bounds for the negative moments of the integrated squared volatility T 0 v t dt, in Subsection 4.4 the upper bounds for derivatives of functionals of the squared volatility process v. Finally, in Subsection 4.5, we complete the proof of Theorem 2.4 using the previous Subsections.

Notations. In order to alleviate the proofs, we introduce some notations specific to this section. Differentiation. For every process Z ǫ , we write (if these derivatives have a meaning):

(i) Z i,t = ∂ i Z ǫ t ∂ǫ i | ǫ=0 , T 3.5
Put prices of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 3, ξ = 30% and ρ = -20%. (ii) the i th Taylor residual by R Z ǫ i,t = Z ǫ t -i j=0 ǫ j j! Z j,t . Generic constants. We keep the same notation C for all non-negative constants (i) depending on universal constants, on a number p ≥ 1 arising in L p estimates, on θ In f , v 0 and K;

3M

(ii) depending in a non decreasing way on κ,

1 1-|ρ| 2 Sup , θ Sup , ξ Sup , ξ Sup ξ In f

and T.

We write A = O(B) when |A| ≤ CB for a generic constant.

Miscellaneous.

(i) We write σ ǫ t = v ǫ t for the volatility for the perturbed process. (ii) if (Z) t∈[0,T] is a càdlàg process, we denote by Z * its running extremum:

Z * t = sup s≤t |Z s |, ∀t ∈ [0, T].
(iii) The L p norm of a random variable is denoted, as usual, by

Z p = E[|Z| p ] 1/p .

Upper bounds for put derivatives.

L 4.1. For every (i, j) ∈ N 2 , there exists a polynomial P with positive coefficients such that:

sup x∈R ∂ i+j P BS ∂x i y j (x, y) ≤ P( √ y) y (2 j+i-1) + 2 .
Proof. Note that it is enough to prove the estimates for j = 0, owing to the relation (2.2). We now take j = 0. For i = 0, the inequality holds because P BS is bounded. Thus T 3.6 Implied Black-Scholes volatilities of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 3, ξ = 30% and ρ = 20%.

3M

21 consider i ≥ 1. Then by differentiating the payoff, one gets:

∂ i P BS ∂x i (x, y) = ∂ i x E[e -T 0 r t dt (K -e x+ T 0 (r t -q t )dt- y 2 + y T W T ) + ] = -∂ i-1 x E[½ (e x+ T 0 (r t -q t )dt- y 2 + y T W T ≤K) e x- T 0 q t dt- y 2 + y T W T ] = -∂ i-1 x E[Ψ(x + G)]
where Ψ is a bounded function (by K) and G is a Gaussian variable with zero mean and variance equal to y. For such a function, we write

E[Ψ(x + G)] = R Ψ(z) e -(z-x) 2 /(2y) √ 2πy dz
and from this, it follows by a direct computation that

∂ i-1 x E[Ψ(x + G)] ≤ C y i-1 2
for any x and y. We have proved the estimate for j = 0 and i ≥ 1.

Positivity of the squared volatility process v.

For a complete review related to time homogeneous CIR processes, we refer the reader to [START_REF] Öing-Jaeschke | A survey and some generalizations of Bessel processes[END_REF]. For time dependent CIR process, see [START_REF] Maghsoodi | Solution of the extended CIR term structure and bond option valuation[END_REF] where the existence and representation using squared Bessel processes are provided.

To prove the positivity of the process v, we show that it can be bounded from below by a suitable time homogeneous CIR process, time scale being the only difference (see definition 5.1.2. in [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF]). The arguments are quite standard, but since we T 3.7 Put prices of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 3, ξ = 30% and ρ = 20%.

= v 0 + A t 0 (κ(θ s -v s )ds + ξ s √ v s dB s ) = v 0 + t 0 ( κ ξ 2 A s (θ A s -ṽs )ds + √ ṽs d Bs ).
Now we apply a comparison result for SDEs twice (see Proposition 5.2.18 in [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]).

1. First, one gets ṽt ≥ n t , where (n s ) s is the (unique) solution of

n t = 0 + t 0 - κ ξ 2 A s n s ds + √ n s d Bs , because v 0 ≥ 0 and κ ξ 2 As (θ A s -x) ≥ -κ ξ 2

As

x, for all x ∈ R and s ∈ [0, A -1 T ]. Of course n t = 0, thus ṽt is non-negative.

2. Secondly, using the non-negativity of ṽ, we only need to compare drift coefficients for the non-negative variable x. Under (P), since

κ ξ 2 A s (θ A s -x) ≥ 1 2 - κ ξ 2 In f x ∀x ≥ 0, ∀s ∈ [0, A -1 T ],
we obtain ṽt ≥ y t for t ∈ [0, A -1 T ] a.s. Moreover, the positivity of y (and consequently that of v) is standard: indeed, y is a 2-dimensional squared Bessel process with a time/space scale change (see [START_REF] Öing-Jaeschke | A survey and some generalizations of Bessel processes[END_REF], or the proof of Lemma 4.3 below).

Upper bound for negative moments of the integrated squared volatility process

T 0 v t dt. L 4.3. Assume (P). Then for every p > 0, one has:

sup 0≤ǫ≤1 E[( T 0 v ǫ t dt) -p ] ≤ C T p .
T 3.9 Put prices of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 3, ξ = 30% and ρ = -50%. Before proving the result, we mention that analogous estimates appear in [START_REF] Bossy | An efficient discretization scheme for one dimensional SDEs with a diffusion coefficient function of the form |x| a[END_REF] (Lemmas A.1 and A.2): some exponential moments are stated under stronger conditions than those in assumption (P). In addition, the uniformity of the estimates w.r.t. ξ (or equivalently w.r.t. ǫ) is not emphasized. In our study, it is crucial to get uniform estimates w.r.t. ǫ.

3M

Proof. Fix p ≥ 1 2 (for 0 < p < 1 2 , we derive the result from the case p = 1 2 using the H ölder inequality). The proof is divided into two steps. We first prove the estimates in the case of constant coefficients κ, θ, ξ with κθ = 1 2 , ǫ = 1 and ξ = 1. Then, using the time change of Lemma 4.2, we derive the result for (v ǫ t ) t . The critical point is to get estimates that are uniform w.r.t. ǫ.

Step 1. Take θ t ≡ θ, ξ t ≡ 1, κθ = 1 2 , ǫ = 1 and consider

dy t = ( 1 2 -κy t )dt + √ y t dB t , y 0 = v 0 ,
for a standard Brownian motion B. We represent y as a time space transformed squared Bessel process (see [START_REF] Öing-Jaeschke | A survey and some generalizations of Bessel processes[END_REF])

y t = e -κt z (e κt -1) 4κ
where z is a 2-dimensional squared Bessel process. Therefore, using a change of variable and the explicit expression of Laplace transform for the integral of z (see [START_REF] Borodin | Handbook of Brownian motion-facts and formulae[END_REF] p.377), one obtains for any u ≥ 0

E[exp(-u T 0 y t dt)] ≤ E[exp(-4ue -2κT (e κT -1) 4κ 0 z s ds)] ≤ cosh( √ 2u(1 -e -κT ) 2κ ) -1 exp(- √ 2ue -κT v 0 tanh( √ 2u(1 -e -κT ) 2κ
)).

Combining this with the identity x -p = 1 Γ(p) ∞ 0 u p-1 e -ux du for x = T 0 y t dt, one gets:

E[( T 0 y t dt) -p ] ≤ 1 Γ(p) ∞ 0 u p-1 cosh( √ 2u(1 -e -κT ) 2κ ) -1 exp(- √ 2ue -κT v 0 tanh( √ 2u(1 -e -κT ) 2κ
))du.

Define the parameter λ 2 = (e κT -1) 2κv 0 and the new variable n

= √ 2u(1-e -κT ) 2κ = v 0 e -κT λ 2 √ 2u. It readily follows that E[( T 0 y t dt) -p ] ≤ C( e κT λ 2 ) 2p ∞ 0 n 2p-1 cosh(n) -1 exp(- tanh(n)n λ 2 )dn,
where C is a constant depending only on v 0 and p. We upper bound the above integral differently according to the value of λ. (ii) If λ ≤ 1, split the integral into two parts, n ≤ arctanh(λ) and n ≥ arctanh(λ). For the first part, simply use n ≥ tanh(n) for any n. For the second part, use tanh(n) ≥ λ and cosh(n) -1 ≤ 1. This gives

(i) If λ ≥ 1, then E[( T 0 y t dt) -p ] ≤ C( e κT λ 2 ) 2p ∞ 0 n 2p-1 cosh(n) -1 dn ≤ Ce 2pκT . ( 4
E[( T 0 y t dt) -p ] ≤C ( e κT λ 2 ) 2p arctanh(λ) 0 n 2p-1 cosh(n) -1 exp(- tanh 2 (n) λ 2 )dn + ( e κT λ 2 ) 2p ∞ arctanh(λ) n 2p-1 exp(- n λ )dn := C[T 1 + T 2 ]. (4.3)
We upper bound the two terms separately. 1. First term T 1 . Using the change of variable m = tanh(n)

λ , one has:

T 1 ≤ e 2pκT λ -4p+1 1 0 arctanh(λm) 2p-1 cosh(arctanh(λm)) exp(-m 2 )dm.
Because of λ ≤ 1, we have the following inequalities for m ∈ [0, 1[: arctanh(λm) ≤ λarctanh(m), cosh(arctanh(λm)) ≤ cosh(arctanh(m)).

Using 2p -1 ≥ 0, it readily follows that

T 1 ≤ ( e 2κT λ 2 ) p 1 0 arctanh(m) 2p-1 cosh(arctanh(m)) exp(-m 2 )dm. (4.4)
2. Second term T 2 . Clearly, we have

T 2 ≤ ( e κT λ 2 ) 2p ∞ 0 n 2p-1 exp(- n λ )dn = ( e 2κT λ 2 ) p ∞ 0 v 2p-1 e -v dv
E[( T 0 v ǫ t dt) -p ] ≤ (ǫξ Sup ) 2p E( A -1 ǫ,T 0 y ǫ s ds) -p ≤ C(ǫξ Sup ) 2p e 2p κ (ǫξ In f ) 2 A -1 ǫ,T (1 + 1 [A -1 ǫ,T ] p ) ≤ Ce 2pκ ξ 2 Sup ξ 2 In f T (ξ 2p Sup + ξ 2p Sup ξ 2p In f 1 T p )
where we have used

ǫ 2 ξ 2 In f T ≤ A -1 ǫ,T ≤ ǫ 2 ξ 2 Sup T.
Note that the upper bound does not depend on ǫ ∈]0, 1]. For ǫ = 0, the upper bound in Lemma 4.3 is also true because (v 0 t ) t is deterministic and

max(v 0 , θ Sup ) ≥ v 0 t ≥ min(v 0 , θ In f ) > 0. (4.8)
4.4. Upper bound for residuals of the Taylor development of g(ǫ) defined in (1.4). Throughout the following paragraph, we assume that (P) is in force. We define the variables:

P ǫ T = T 0 ρ t (σ ǫ t -σ 0,t )dB t - T 0 ρ 2 t 2 (v ǫ t -v 0,t )dt, Q ǫ T = T 0 (1 -ρ 2 t )(v ǫ t -v 0,t )dt. Notice that (x 0 + T 0 ρ t √ v 0,t dB t - T 0 ρ 2 t 2 v 0,t dt, T 0 (1-ρ 2 t )v 0,t dt)+(P 1 T , Q 1 T ) = (x 0 + T 0 ρ t v 1 t dB t - T 0 ρ 2 t 2 v 1 t dt, T 0 (1 -ρ 2 t )v 1 t dt).
The main result of this subsection is the following proposition, the statement of which uses the notation introduced at the beginning of Section 4. P 4.4. One has the following estimates for every p ≥ 1

P 1 T p ≤C(ξ Sup √ T) √ T, R P 1 2,T p ≤C(ξ Sup √ T) 3 √ T, R (P 1 ) 2 2,T p ≤C(ξ Sup √ T) 3 T, Q 1 T p ≤C(ξ Sup √ T)T, R Q 1 2,T p ≤C(ξ Sup √ T) 3 T, R (Q 1 ) 2 2,T p ≤C(ξ Sup √ T) 3 T 2 , R P 1 Q 1 2,T p ≤C(ξ Sup √ T) 3 T 3 2 .
To estimate the derivatives and the residuals for the variables P ǫ T and Q ǫ T , we need first to prove the existence of the derivatives and the residuals of the volatility process σ ǫ t = v ǫ t and its square v ǫ . Finally we prove Proposition 4.4. 4.4.1. Upper bounds for derivatives of σ ǫ and v ǫ . Under assumption (P), the volatility process σ ǫ t is governed by the SDE:

dσ ǫ t = (( κθ t 2 - ǫ 2 ξ 2 t 8 ) 1 σ ǫ t - κ 2 σ ǫ t )dt + ǫξ t 2 dB t , σ ǫ 0 = √ v 0 , (4.9)
where we have used Ito's Lemma and positivity of v ǫ t (see Lemma 4.2). In order to estimate R σ ǫ 0,t , we are going to prove that it verifies a linear equation (Lemma 4.5) from which we deduce an a priori upper bound (Proposition 4.6). We iterate the same analysis for the residuals R σ ǫ 1,t (Proposition 4.7) and R σ ǫ 2,t (Proposition 4.8). Analogously, we give upper bounds for the residuals of v ǫ t (Proposition 4.9). L 4.5. Under (P), the process (R σ ǫ 0,t = σ ǫ t -σ 0 t ) 0≤t≤T is given by

R σ ǫ 0,t = U ǫ t t 0 (U ǫ s ) -1 (- ǫ 2 ξ 2 s 8σ 0,s ds + ǫξ s 2 dB s ),
where

dU ǫ t = -α ǫ t U ǫ t dt, U ǫ t = 1, α ǫ t = ( κθ t 2 - ǫ 2 ξ 2 t 8 ) 1 σ ǫ t σ 0,t + κ 2 .
Proof. From the definition (σ 0,t ) t = (σ 0 t ) t and the equation (4.9), one obtains the SDE

dσ 0,t = ( κθ t 2σ 0,t - κ 2 σ 0,t )dt, σ 0,0 = √ v 0 .
Substitute this equation in (4.9) to obtain

dR σ ǫ 0,t = -α ǫ t R σ ǫ 0,t dt - ǫ 2 ξ 2 t 8σ 0,t dt + ǫξ t 2 dB t , R σ ǫ 0,0 = 0. (4.10)
Note that R σ ǫ 0,. is the solution of a linear SDE. Hence, it can be explicitly represented using the process U ǫ (see Th. 52 in [START_REF] Protter | Stochastic integration and differential equations[END_REF]):

R σ ǫ 0,t = U ǫ t t 0 (U ǫ s ) -1 (- ǫ 2 ξ 2 s 8σ 0,s ds + ǫξ s 2 dB s ).
P 4.6. Under (P), for every p ≥ 1 one has

(R σ ǫ 0,. ) * t p ≤ Cǫξ Sup √ t.
In particular, the application ǫ σ ǫ t is continuous5 at ǫ = 0 in L p . Proof. At first sight, the proof seems to be straightforward from Lemma 4.5. But actually, the difficulty lies in the fact that one can not uniformly in ǫ upper bound U ǫ t in L p (because of the term with 1/σ ǫ t in α ǫ t ). Using Lemma 4.5 and Ito's formula for the product (U ǫ t ) -1 (

t 0 ǫξ s 2 dB s ), one has R σ ǫ 0,t = U ǫ t t 0 (U ǫ s ) -1 (- ǫ 2 ξ 2 s 8σ 0,s ds) + t 0 ǫξ s 2 dB s -U ǫ t t 0 ( s 0 ǫξ u 2 dB u )d(U ǫ s ) -1 .
Under (P), one has α ǫ t ≥ κ/2 > 0, which implies that t → U ǫ t is decreasing and t → (U ǫ t ) -1 is increasing. Thus, 0 ≤ U ǫ t (U

ǫ 2 ξ 2 s 8σ 0,s + ǫσ 1,s (( α ǫ s σ 0,s - κ 2σ 0,s )R σ ǫ 0,s + ǫ 2 ξ 2 s 8σ 0,s ))ds.
Moreover, for every p ≥ 1, one has

(R σ ǫ 1,. ) * t p ≤ C(ǫξ Sup √ t) 2 .
In particular, the application ǫ σ ǫ t is C 1 at ǫ = 0 in L p sense with the first derivative at ǫ = 0 equal to σ 1,t (justifying a posteriori the definition R σ ǫ 1,. ). Proof. From Equations (4.10) and (4.12), it readily follows that

dR σ ǫ 1,t = -α ǫ t R σ ǫ 1,t dt -ǫσ 1,t (α ǫ t - κθ t 2(σ 0,t ) 2 - κ 2 )dt - ǫ 2 ξ 2 t 8σ 0,t dt, R σ ǫ 1,0 = 0.
Because of the identity

-(α ǫ t - κθ t 2(σ 0,t ) 2 - κ 2 ) = (( α ǫ t σ 0,t - κ 2σ 0,t )R σ ǫ 0,t + ǫ 2 ξ 2 t 8(σ 0,t ) 2 ),
one deduces the equality

R σ ǫ 1,t = U ǫ t t 0 (U ǫ s ) -1 (- ǫ 2 ξ 2 s 8σ 0,s + ǫσ 1,s (( α ǫ s σ 0,s - κ 2σ 0,s )R σ ǫ 0,s + ǫ 2 ξ 2 s 8(σ 0,s ) 2 ))ds.
Then

|R σ ǫ 1,t | ≤ t 0 U ǫ t (U ǫ s ) -1 ( ǫ 2 ξ 2 s 8σ 0,s + ǫ|σ 1,s |(( α ǫ s σ 0,s + κ 2σ 0,s )|R σ ǫ 0,s | + ǫ 2 ξ 2 s 8(σ 0,s ) 2 ))ds ≤ t 0 U ǫ t (U ǫ s ) -1 ( ǫ 2 ξ 2 s 8σ 0,s + ǫ|σ 1,s |( κ 2σ 0,s |R σ ǫ 0,s | + ǫ 2 ξ 2 s 8(σ 0,s ) 2 ))ds + ǫ t 0 U ǫ t (U ǫ s ) -1 α ǫ s σ 0,s |σ 1,s R σ ǫ 0,s |ds ≤ t 0 ( ǫ 2 ξ 2 s 8σ 0,s + ǫ|σ 1,s |( κ 2σ 0,s |R σ ǫ 0,s | + ǫ 2 ξ 2 s 8(σ 0,s ) 2 ))ds + ǫ( σ 1,. R σ ǫ 0,. σ 0,. ) * t ,
where we have used

U ǫ t (U ǫ s ) -1 ≤ 1 for every s ∈ [0, t] and U ǫ t t 0 α ǫ s (U ǫ s ) -1 ds = 1 -U ǫ t ≤
R σ ǫ 2,t = U ǫ t t 0 (U ǫ s ) -1 [ǫ 2 (( α ǫ s σ 0,s - κ 2σ 0,s )R σ ǫ 0,s + ǫ 2 ξ 2 s 8(σ 0,s ) 2 )( σ 2,s 2 - (σ 1,s ) 2 σ 0,s ) + ǫ(( α ǫ s σ 0,s - κ 2σ 0,s )R σ ǫ 1,s + ǫ 2 ξ 2 s 8(σ 0,s ) 2 )σ 1,s ]ds.
Moreover, for every p ≥ 1, one has

(R σ ǫ 2,. ) * t p ≤ C(ǫξ Sup √ t) 3 .
In particular, the application ǫ σ ǫ t is C 2 at ǫ = 0 in L p sense with the second derivative at ǫ = 0 equal to σ 2,t .

Proof. The equality is easy to check. The estimate is proved in the same way as in the proof of Proposition 4.7, we therefore skip the details.

C 4.9. The application ǫ v ǫ t is C 2 at ǫ = 0 in L p sense. The residuals for the squared volatility satisfy the following inequalities: for every p ≥ 1, one has

(R v ǫ 0,. ) * t p ≤ Cǫξ Sup √ t, (R v ǫ 1,. ) * t p ≤ C(ǫξ Sup √ t) 2 , (R v ǫ 2,. ) * t p ≤ C(ǫξ Sup √ t) 3 .
Proof. Note that v ǫ t = (σ ǫ t ) 2 = (σ 0,t + R σ ǫ 0,t ) 2 = v 0,t + 2σ 0,t R σ ǫ 0,t + (R σ ǫ 0,t ) 2 . Thus, we have R v ǫ 0,t = 2σ 0,t R σ ǫ 0,t +(R σ ǫ 0,t ) 2 , which leads to the required estimate using σ 0,t ≤ max( √ v 0 , θ Sup ) and Proposition 4.6. The other estimates are proved analogously using Propositions 4.7 and 4.8 and Inequalities (4.13) and (4.15).

Proof of Proposition 4.4. We can write

P 1 T = T 0 ρ t R σ 1 0,t dB t - T 0 ρ 2 t 2 R v 1 0,t dt, R P 1 2,T = T 0 ρ t R σ 1 2,t dB t - T 0 ρ 2 t 2 R v 1 2,t dt.
Then, using Propositions 4.6, 4.8 and Corollary 4.9, we prove the two first estimates of Proposition 4.4. The others inequalities are proved in the same way.

4.5. Proof of Theorem 2.4. For convenience, we introduce the following notation for λ ∈ [0, 1]:

PBS (λ) =P BS x 0 + T 0 ρ t ((1 -λ) √ v 0,t + λ v 1 t )dB t - T 0 ρ 2 t 2 ((1 -λ)v 0,t + λv 1 t )dt , T 0 (1 -ρ 2 t )((1 -λ)v 0,t + λv 1 t )dt , ∂ i+j PBS ∂x i y j (λ) = ∂ i+j P BS ∂x i y j x 0 + T 0 ρ t ((1 -λ) √ v 0,t + λ v 1 t )dB t - T 0 ρ 2 t 2 ((1 -λ)v 0,t + λv 1 t )dt , T 0 (1 -ρ 2 t )((1 -λ)v 0,t + λv 1 t )dt .
Notice that PBS (see (2.5)) is a particular case of PBS for λ = 0:

PBS = PBS (0), ∂ i+j PBS ∂x i y j = ∂ i+j PBS ∂x i y j (0).
Now, we represent the error E in (2.12) using the previous notations. A second order Taylor expansion leads to

g(1) = E( PBS (1)) = E( PBS (0) + ∂ λ PBS (0) + 1 2 ∂ 2 λ PBS (0) + 1 0 dλ (1 -λ) 2 2 ∂ 3 λ PBS (λ)).
The first term E( PBS (0)) is equal to (2.6). Approximations of the three above derivatives contribute to the error E. 1. We have

E(∂ λ PBS (0)) = E( ∂ PBS ∂x P 1 T + ∂ PBS ∂y Q 1 T )
. These two terms are equal to (2.7) and (2.8) plus an error equal to

E( ∂ PBS ∂x R P 1 2,T + ∂ PBS ∂y R Q 1 2,T ).
2. Regarding the second derivatives, we have

E( 1 2 ∂ 2 λ PBS (0)) = E( 1 2 ∂ 2 PBS ∂x 2 (P 1 T ) 2 + 1 2 ∂ 2 PBS ∂y 2 (Q 1 T ) 2 + ∂ 2 PBS ∂xy P 1 T Q 1 T )
. These terms are equal to (2.9), (2.10) and (2.11), plus an error equal to

E( 1 2 ∂ 2 PBS ∂x 2 R (P 1 ) 2 2,T + 1 2 ∂ 2 PBS ∂y 2 R (Q 1 ) 2 2,T + ∂ 2 PBS ∂xy R P 1 Q 1 2,T ).
3. The last term with ∂ 3 λ PBS is neglected and thus is considered as an error.

To sum up, we have shown that

E = 1 i=0 E[ ∂ 1 PBS ∂x i y 1-i (0)R (P 1 ) i (Q 1 ) 1-i 2,T ] + 2 i=0 C i 2 2 E[ ∂ 2 PBS ∂x i y 2-i (0)R (P 1 ) i (Q 1 ) 2-i 2,T ] + 1 0 (1 -λ) 2 2 3 i=0 C i 3 E[ ∂ 3 PBS ∂x i y 3-i (λ)(P 1 T ) i (Q 1 T ) 3-i ]dλ.
Using Lemma 4.1 and Assumption (R), one obtains for all λ ∈ [0, 1]

∂ i+j PBS ∂x i y j (λ) 2 ≤ C ( T 0 ((1 -λ)v 0,t + λv 1 t )dt) -(2 j+i-1) + 2 4 ≤ C((1 -λ) ( T 0 v 0,t dt) -(2 j+i-1) + 2 4 + λ ( T 0 v 1 t dt) -(2 j+i-1) + 2 4 )
where we have applied a convexity argument. Finally, apply Lemma 4.3 with ǫ = 0 and ǫ = 1 to conclude that

∂ i+j PBS ∂x i y j (λ) 2 ≤ C ( √ T) (2 j+i-1) + , uniformly w.r.t. λ ∈ [0, 1]. Combining this with Proposition 4.4 yields that |E| ≤ C        1 i=0 (ξ Sup √ T) 3 T 1-i/2 ( √ T) 1-i + 2 i=0 (ξ Sup √ T) 3 T 2-i/2 ( √ T) 3-i + 3 i=0 (ξ Sup √ T) 3 T 3-i/2 ( √ T) 5-i        ≤ C(ξ Sup √ T) 3 √ T.
Theorem 2.4 is proved.

5. Proof of Proposition 2.2 and Theorem 2.3.

Preliminary results.

In this section, we bring together the results (and their proofs) which allow us to derive the explicit terms in the formula (2.13). In the following, α t (resp. β t ) is a square integrable and predictable process (resp. deterministic) and l is a smooth function with derivatives having, at most, exponential growth. For the next Malliavin calculus computations, we freely use standard notations from [START_REF] Nualart | Malliavin calculus and related topics[END_REF].

L 5.1.

(Lemma 1.2.1 in [Nua06]) Let G ∈ D 1,∞ (Ω). One has E[G t 0 α s dB s ] = E[ t 0 α s D B s (G)ds],
where

D B (G) = (D B s (G)) s≥0 is the first Malliavin derivative of G w.r.t. B. Taking G = l( T 0 ρ t √ v 0,t dB t )
gives the following result. L 5.2. One has:

E[( T 0 α t dB t )l( T 0 ρ t √ v 0,t dB t )] = E[( T 0 ρ t √ v 0,t α t dt)l (1) ( T 0 ρ t √ v 0,t dB t )].
L 5.3. For any deterministic integrable function f and any continuous semimartingale Z vanishing at t=0, one has:

T 0 f (t)Z t dt = T 0 ω (0, f ) t,T dZ t .
Proof. This is an application of the It ô formula to the product ω (0, f ) t,T Z t . L 5.4. One has:

E[l( T 0 ρ t √ v 0,t dB t ) T 0 β t v 1,t dt] = ω (κ,ρξv 0,. ),(-κ,β) 0,T E[l (1) ( T 0 ρ t √ v 0,t dB t )], E[l( T 0 ρ t √ v 0,t dB t ) T 0 β t v 2 1,t dt] = ω (2κ,ξ 2 v 0,. ),(-2κ,β) 0,T E[l( T 0 ρ t √ v 0,t dB t )] + 2ω
(κ,ρξv 0,. ),(κ,ρξv 0,. ),(-2κ,β) 0,T

E[l (2) ( T 0 ρ t √ v 0,t dB t )], E[l( T 0 ρ t √ v 0,t dB t ) T 0 β t v 2,t dt] = ω (κ,ρξv 0,. ),(0,ρξ),(-κ,β) 0,T E[l (2) ( T 0 ρ t √ v 0,t dB t )].
Proof. Using Lemmas 5.2 ( f (t) = e -κt β t , Z t = t 0 e κs ξ s √ v 0,s dB s ) and 5.3, one has:

E[l( T 0 ρ t √ v 0,t dB t ) T 0 β t v 1,t dt] = E[l( T 0 ρ t √ v 0,t dB t ) T 0 e -κt β t t 0 e κs ξ s √ v 0,s dB s dt] = E[l( T 0 ρ t √ v 0,t dB t ) T 0 ω (-κ,β) t,T e κt ξ t √ v 0,t dB t ] = E[l (1) ( T 0 ρ t √ v 0,t dB t )] T 0 ω (-κ,β) t,T
e κt ρ t ξ t v 0,t dt, which gives the first equality. The second and the third are proved in the same way.

L 5.5. One has

E[ ∂ i+j PBS ∂x i y j ] = ∂ i+j P BS ∂x i y j (x 0 , T 0 v 0,t dt).
Proof. One has

E[ ∂ i PBS ∂x i ] = ∂ i x=x 0 E[P BS (x 0 + T 0 ρ t √ v 0,t dB t - T 0 ρ 2 t 2 v 0,t dt, T 0 (1 -ρ 2 t )v 0,t dt)] = ∂ i P BS ∂x i (x 0 , T 0 v 0,t dt).
Since PBS verifies the following relation

∂ PBS ∂y = 1 2 ( ∂ 2 PBS ∂x 2 - ∂ PBS ∂x ), (5.1)
we immediately obtain the result.

Proof of Proposition 2.2. One has

E[ ∂ PBS ∂x ( T 0 ρ t ( v 1,t 2(v 0,t ) 1 2 + v 2,t 4(v 0,t ) 1 2 - v 2 1,t 8(v 0,t ) 3 2 )dB t - T 0 ρ 2 t 2 (v 1,t + v 2,t 2 )dt)] = E[ 1 2 ( ∂ 2 PBS ∂x 2 - ∂ PBS ∂x ) T 0 ρ 2 t (v 1,t + v 2,t 2 )dt] -E[ ∂ 2 PBS ∂x 2 T 0 ρ 2 t v 2 1,t 8v 0,t dt] = E[ ∂ PBS ∂y T 0 ρ 2 t (v 1,t + v 2,t 2 )dt] -E[ ∂ 2 PBS ∂x 2 T 0 ρ 2 t v 2 1,t 8v 0,t dt],
where we have used Lemma 5.2 at the first equality and identity (5.1) at the second one. Plugging this relation into the approximation (2.12) and summing the second and third line, one has

g(1) =E[ PBS ] + E[ ∂ PBS ∂y T 0 (v 1,t + v 2,t 2 )dt] -E[ ∂ 2 PBS ∂x 2 T 0 ρ 2 t v 2 1,t 8v 0,t dt] + 1 2 E[ ∂ 2 PBS ∂x 2 ( T 0 ρ t v 1,t 2(v 0,t ) 1 2 dB t - T 0 ρ 2 t 2 v 1,t dt) 2 ] + 1 2 E[ ∂ 2 PBS ∂y 2 ( T 0 (1 -ρ 2 t )v 1,t dt) 2 ] +E[ ∂ 2 PBS ∂xy ( T 0 (1 -ρ 2 t )v 1,t dt)( T 0 ρ t v 1,t 2(v 0,t ) 1 2 dB t - T 0 ρ 2 t 2 v 1,t dt)] + E.
(5.2)

In addition, one has

-E[ ∂ 2 PBS ∂x 2 T 0 ρ 2 t v 2 1,t 8v 0,t dt] + 1 2 E[ ∂ 2 PBS ∂x 2 ( T 0 ρ t v 1,t 2(v 0,t ) 1 2 dB t - T 0 ρ 2 t 2 v 1,t dt) 2 ] = E[ ∂ 2 PBS ∂x 2 T 0 ( t 0 ρ s v 1,s 2(v 0,s ) 1 2 dB s - t 0 ρ 2 s 2 v 1,s ds)(ρ t v 1,t 2(v 0,t ) 1 2 dB t - ρ 2 t 2 v 1,t dt)] = E[ 1 2 ( ∂ 3 PBS ∂x 3 - ∂ 2 PBS ∂x 2 ) T 0 ( t 0 ρ s v 1,s 2(v 0,s ) 1 2 dB s - t 0 ρ 2 s 2 v 1,s ds)ρ 2 t v 1,t dt] = E[ ∂ 2 PBS ∂xy T 0 ( t 0 ρ s v 1,s 2(v 0,s ) 1 2 dB s - t 0 ρ 2 s 2 v 1,s ds)ρ 2 t v 1,t dt],
where we have used Ito's Lemma for the square at the first equality, Lemma 5.2 at the second and Identity (5.1) at the third one. Substituting this relation in the approximation (5.2) and summing the second and fourth line, one gets

g(1) =E[ PBS ] + E[ ∂ PBS ∂y T 0 (v 1,t + v 2,t 2 )dt] +E[ ∂ 2 PBS ∂xy ( T 0 ( t 0 ρ s v 1,s 2(v 0,s ) 1 2 dB s - t 0 ρ 2 s 2 v 1,s ds)ρ 2 t v 1,t dt + ( T 0 (1 -ρ 2 t )v 1,t dt)( T 0 ρ t v 1,t 2(v 0,t ) 1 2 dB t - T 0 ρ 2 t 2 v 1,t dt))] + 1 2 E[ ∂ 2 PBS ∂y 2 ( T 0 (1 -ρ 2 t )v 1,t dt) 2 ] + E.
(5.3)

We now study the second term of (5.3). In the computations below, we use Ito's Lemma for the second equality, Lemma 5.2 and Identity (5.1) for the third equality and Lemma 5.1 (G = ∂ 2 PBS ∂xy v 1,t ) for the fourth one; it gives , a 2,T = ω (κ,ρξv 0,. ),(0,ρξ),(-κ,1) 0,T .

A =E[ ∂ 2 PBS ∂xy ( T 0 ( t 0 ρ s v 1,s 2(v 0,s ) 1 2 dB s - t 0 ρ 2 s 2 v 1,s ds)ρ 2 t v 1,t dt + ( T 0 (1 -ρ 2 t )v 1,t dt)( T 0 ρ t v 1,t 2(v 0,t ) 1 2 dB t - T 0 ρ 2 t 2 v 1,t dt))] =E[ ∂ 2 PBS ∂xy ( T 0 ( t 0 ρ s v 1,s 2(v 0,s ) 1 2 dB s - t 0 ρ 2 s 2 v 1,s ds)(ρ 2 t + 1 -ρ 2 t )
Actually, the result is an immediate application of Lemma 5.4 and Lemma 5.5.

Step 2: We show the equality (κ,ρξv 0,. ),(-κ,1),(κ,ρξv 0,. ),(-κ,1) 0,T + 2ω

(κ,ρξv 0,. ),(κ,ρξv 0,. ),(-κ,1),(-κ,1) 0,T Putting f (t) = ρ t ξ t v 0,t e kt and g(t) = e -kt in the previous equality readily gives b 2,T = a 2 1,T 2 , which finishes the proof.

Conclusion.

We have established an approximation pricing formula for call/put options in the time dependent Heston models. We prove that the error is of order 3 w.r.t. the volatility of volatility and 2 w.r.t. the maturity. In practice, taking the Fourier method as a benchmark, the accuracy is excellent for a large range of strikes and maturities. In addition, the computational time is about 100 to 1000 times smaller than using an efficient Fourier method.

Following the arguments in [START_REF] Benhamou | Smart expansion and fast calibration for jump diffusion[END_REF], our formula extends immediately to other payoffs depending on S T (note that the identities (2.2) and (5.1) are valid for any payoff of this type). As explained in [START_REF] Benhamou | Smart expansion and fast calibration for jump diffusion[END_REF], the smoother the payoff, the higher the error order w.r.t. T; the less smooth the payoff, the lower the error order w.r.t. T. For digital options, the error order w.r.t. T becomes 3/2 instead of 2.

Extensions to exotic options and to the third order expansion formula w.r.t. the volatility of volatility are left for further research.

Appendix: closed formulas in Heston model.

There are few closed representations for the call/put prices written on the asset S t = e t 0 (r s -q s )ds e X t in the Heston model (defined in (1.1) and (1.2)). We focus on the Heston formula [START_REF] Heston | A closed-form solutions for options with stochastic volatility with applications to bond and currency options[END_REF] and on the Lewis formula [START_REF] Lewis | Option valuation under stochastic volatility with Mathematica code[END_REF]. Both of them rely on the knowledge of the characteristic function of the log-asset price (X t ) t and on Fourier transform-based approaches.

(i) In [START_REF] Heston | A closed-form solutions for options with stochastic volatility with applications to bond and currency options[END_REF], Heston obtains a representation in a Black-Scholes form:

Call Heston (t, S t , v t ; T, K) = S t e -T t q s ds P 1 -Ke -T t r s ds P 2 , where both probabilities P 1 and P 2 are equal to a one-dimensional integral of characteristic functions.

(ii) In [START_REF] Lewis | Option valuation under stochastic volatility with Mathematica code[END_REF], Lewis takes advantage of the generalized Fourier transform, by using an integration along a straight line in the complex plane parallel to the real axis. It is important to detect the strip where the integration is safe. Lewis suggests the use of complex numbers z such that Im(z) = 1 2 . His formula writes Call Heston (t, S t , v t ; T, K) = S t e -T t q s ds -Ke -T t r s ds

2π i 2 +∞ i 2 -∞
e -izX φ T (-z) dz z 2iz where X = log S t e -T t qsds Ke -T t rsds and φ T (z) = E(e z(X T -X t ) |F t ). Then, the above integral is evaluated by numerical integration. Using PDE arguments in combination with affine models, we can obtain an explicit formula for φ T (z) in the case of constant Heston parameters. In addition, it can be computed without discontinuities in z, following the arguments in [START_REF] Jackel | Not-So-Complex Logarithms in the Heston Model[END_REF]. For piecewise constant parameters, the characterictic function φ T (z) can be computed recursively using nested Riccatti equations with constant coefficients: we refer to the work by Mikhailov and Nogel [START_REF] Mikhailov | Hestons stochastic volatility model implementation, calibration and some extensions[END_REF].

In our numerical tests, we prefer the Lewis formula which gives better numerical results, in particular for very small or very large strikes, compared to the Heston formula.

  t, T) =

	e κT θ+e κt ((κt-κT-1)θ+κ(T-t)v 0 ) κ (e κt -e κT )(e κt (θ-2v 0 )-e κT θ) , 2κ t (T) = -e tu +e Tu κ,v 0 (θ, t, T) = h 2 and ωu u , ω0,u t (T) = e Tu (-tu+Tu-1)+e tu , u 2 Proof. According to Theorem 2.3, one has : , ωu,u t (T) =	(e tu -e Tu ) 2u 2	2	.

  Error in implied Black-Scholes volatilities (in bp) between the closed formulas (see appendix) of the two models M 1 and M 2 expressed as relative strikes. Maturity is equal to one year.

	strikes K	80%	90%	100%	110%	120%
	model M 1	20.12%	19.64%	19.50%	19.62%	19.92%
	model M 2	20.11%	19.65%	19.51%	19.62%	19.92%
	errors (bp)	0.69	-0.35	-0.81	-0.42	0.34

T 2.1

  Set of maturities and strikes used for the numerical tests.

	T/K
	3M	70	80	90	100	110	120	125	130
	6M	60	70	80	100	110	130	140	150
	1Y	50	60	80	100	120	150	170	180
	2Y	40	50	70	100	130	180	210	240
	3Y	30	40	60	100	140	200	250	290
	5Y	20	30	60	100	150	250	320	400
	7Y	10	30	50	100	170	300	410	520
	10Y	10	20	50	100	190	370	550	730

T 3.1

  Put prices of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 3, ξ = 30% and ρ = 0%.

					T 3.3
	3M	30.00	20.08	10.87	4.22	1.14	0.24	0.10	0.04
		30.00	20.08	10.87	4.22	1.14	0.24	0.10	0.04
		0.03	-0.11	0.06	0.08	0.09	-0.15	-0.14	-0.07
	6M	40.01	30.07	20.52	6.20	2.72	0.40	0.14	0.05
		40.01	30.08	20.52	6.19	2.71	0.40	0.14	0.05
		0.05	-0.16	-0.18	0.26	0.26	-0.34	-0.29	-0.08
	1Y	50.01	40.11	21.84	9.12	3.08	0.51	0.15	0.09
		50.01	40.11	21.84	9.11	3.07	0.52	0.16	0.09
		0.04	-0.21	0.06	0.44	0.23	-0.51	-0.29	-0.12
	2Y	60.03	50.20	32.08	13.26	4.71	0.79	0.28	0.11
		60.03	50.20	32.08	13.26	4.71	0.79	0.29	0.11
		0.00	-0.18	-0.03	0.38	0.17	-0.43	-0.29	-0.06
	3Y	70.02	60.15	41.70	16.39	5.73	1.21	0.36	0.15
		70.02	60.15	41.70	16.39	5.73	1.21	0.37	0.15
		0.01	-0.09	-0.08	0.27	0.11	-0.31	-0.22	-0.07
	5Y	80.01	70.15	43.80	21.26	8.50	1.61	0.58	0.21
		80.01	70.15	43.80	21.26	8.50	1.61	0.58	0.21
		0.01	-0.04	0.01	0.15	0.08	-0.19	-0.15	-0.04
	7Y	90.00	70.42	53.15	25.14	9.32	1.97	0.66	0.26
		90.00	70.42	53.15	25.14	9.32	1.97	0.67	0.26
		0.00	-0.04	-0.01	0.09	0.04	-0.14	-0.10	-0.03
	10Y	90.01	80.23	55.22	29.92	11.49	2.62	0.84	0.33
		90.01	80.23	55.22	29.92	11.49	2.62	0.84	0.33
		0.00	-0.02	0.01	0.06	0.03	-0.09	-0.07	-0.02
		Price Approximation-True Price Spot	× 10000

  Put prices of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 6%, κ = 10, ξ = 1 and ρ = -50%.

					T 3.11
	3M	30.05	20.30	11.28	4.35	0.95	0.13	0.04	0.01
		30.04	20.35	11.31	4.33	0.87	0.08	0.05	0.05
		0.99	-4.95	-2.80	2.41	7.37	4.62	-0.31	-3.51
	6M	40.06	30.28	20.96	6.40	2.54	0.21	0.05	0.01
		40.05	30.32	21.02	6.36	2.47	0.15	0.06	0.06
		0.92	-3.90	-5.83	3.18	6.51	5.23	-1.28	-4.72
	1Y	50.08	40.31	22.37	9.29	2.71	0.24	0.04	0.02
		50.07	40.33	22.40	9.26	2.65	0.21	0.07	0.06
		0.41	-2.60	-2.58	2.39	5.95	3.19	-2.52	-3.89
	2Y	60.10	50.39	32.54	13.33	4.18	0.41	0.09	0.02
		60.10	50.40	32.56	13.31	4.14	0.38	0.10	0.05
		-0.01	-1.58	-1.82	1.35	3.67	2.26	-1.17	-2.89
	3Y	70.06	60.28	42.09	16.38	5.09	0.71	0.13	0.03
		70.05	60.28	42.11	16.37	5.06	0.69	0.14	0.06
		0.17	-0.73	-1.46	0.94	2.74	2.19	-0.86	-2.22
	5Y	80.04	70.25	44.15	21.15	7.76	1.02	0.26	0.06
		80.03	70.25	44.16	21.15	7.74	1.01	0.27	0.08
		0.11	-0.36	-0.57	0.61	1.72	1.49	-0.38	-1.68
	7Y	90.00	70.56	53.46	24.96	8.45	1.31	0.32	0.09
		90.00	70.57	53.46	24.96	8.44	1.30	0.32	0.10
		0.06	-0.44	-0.47	0.47	1.42	1.16	-0.46	-1.42
	10Y	90.02	80.31	55.47	29.67	10.51	1.85	0.44	0.13
		90.02	80.31	55.48	29.67	10.50	1.84	0.45	0.14
		0.03	-0.19	-0.20	0.36	1.09	0.95	-0.42	-1.23
									.2)

  Scholes volatilities of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 3%, κ = 2, ξ = 40% and ρ = 0%. Put prices of the closed formula, of the approximation formula and related errors (in bp), expressed as a function of maturities in fractions of years and relative strikes. Parameters: θ = 3%, κ = 2, ξ = 40% and ρ = 0%.

								T 3.12 T 3.13
	3M Implied Black-3M 30.00 23.27% 30.00	20.06 21.25% 20.07	10.67 19.59% 10.67	3.76 18.86% 3.76	0.88 19.47% 0.88	0.18 20.64% 0.18	0.08 21.25% 0.08	0.03 21.85% 0.03
		22.35% 0.12	21.48% -0.46	19.56% 0.24	18.85% 0.32	19.43% 0.45	20.83% -0.88	21.48% -0.58	21.94% -0.13
	6M	92.35 40.01	-22.93 30.05	2.23 20.32		1.62 5.13	3.51 1.91	-18.90 0.25	-22.93 0.09	-9.79 0.04
	6M	24.10% 40.00	22.05% 30.05	20.22% 20.35	18.21% 5.11	18.68% 1.89	20.75% 0.27	21.78% 0.10	22.72% 0.03
		22.52% 0.29	22.26% -0.37	20.50% -2.15	18.14% 1.99	18.59% 2.06	21.16% -2.72	22.10% -1.01	22.51% 0.29
	1Y	158.79 50.01	-20.74 40.06	-28.41 20.98		7.08 7.01	9.08 1.73	-40.69 0.25	-32.03 0.08	20.96 0.05
	1Y	23.96% 50.00	22.01% 40.06	18.89% 21.00	17.60% 6.95	18.51% 1.73	20.84% 0.29	22.23% 0.08	22.85% 0.04
		22.20% 0.55	22.14% -0.27	18.99% -1.80	17.45% 5.92	18.48% 0.70	21.42% -4.43	22.20% 0.07	22.30% 0.98
	2Y	175.41 60.01	-12.81 50.07	-10.17 30.89		14.90 9.70	2.60 2.28	-57.41 0.28	2.27 0.10	54.89 0.04
	2Y	22.72% 60.01	21.05% 50.07	18.61% 30.93	17.24% 9.62	18.04% 2.29	20.26% 0.33	21.42% 0.10	22.42% 0.03
		21.40% 0.63	21.20% -0.39	18.83% -4.17	17.10% 8.03	18.06% -0.51	20.72% -4.44	21.32% 0.40	21.42% 1.61
	3Y	132.35 70.01	-14.49 60.04	-22.09 40.60		14.34 11.82	-1.35 2.56	-46.40 0.36	9.96 0.10	100.04 0.04
	3Y	22.44% 70.00	20.84% 60.04	18.66% 40.64	17.16% 11.75	17.88% 2.58	19.60% 0.41	20.84% 0.09	21.67% 0.03
		20.74% 0.36	20.67% 0.28	18.93% -4.17	17.06% 6.94	17.91% -1.42	19.96% -4.48	20.67% 0.71	20.79% 1.54
	5Y	170.16 80.00	16.92 70.03	-27.04 41.45		10.16 15.21	-3.17 3.75	-36.03 0.38	16.92 0.11	87.99 0.04
	5Y	21.56% 80.00	20.09% 70.02	17.86% 41.47	17.16% 15.16	17.61% 3.75	19.08% 0.40	19.94% 0.11	20.75% 0.02
		20.03% 0.15	19.88% 0.26	17.92% -1.76	17.10% 4.61	17.62% -0.34	19.28% -2.79	19.83% 0.58	20.03% 1.24
	7Y	153.81 90.00	20.49 70.07	-5.89 51.04		5.27 17.97	-0.54 3.78	-19.86 0.39	11.43 0.10	72.25 0.03
	7Y	21.93% 90.00	19.01% 70.08	17.88% 51.06	17.17% 17.94	17.60% 3.80	18.76% 0.41	19.54% 0.09	20.16% 0.02
		19.51% 0.01	19.09% -0.24	17.95% -1.63	17.14% 3.25	17.62% -1.07	18.88% -1.94	19.39% 0.70	19.58% 1.00
	10Y	241.42 90.00	-7.47 80.03	-6.53 51.95		3.16 21.43	-1.53 4.55	-12.47 0.46	14.41 0.10	58.41 0.03
	10Y	20.21% 90.00	18.92% 80.03	17.58% 51.96	17.20% 21.41	17.53% 4.56	18.42% 0.48	19.09% 0.09	19.61% 0.02
		19.24% 0.04	18.88% 0.06	17.60% -0.64	17.18% 2.14	17.54% -0.80	18.49% -1.44	18.97% 0.61	19.16% 0.76
		96.63	4.46	-1.61		1.76	-0.93	-7.64	11.97	44.80
								T 3.14
				Equivalent averaged parameters.
	Combining (4.3), (4.4) and (4.5), we obtain E[( T v0 θ		ξ	ρ
		3M 6M 1Y	4 % 3.97 % 3.28 %					(e κT -1) 2κv 0 ≥ T 30 % 4.04 % 4 % 30.12 % 2v 0 , which gives -20 % -19.93 % 4.38 % 30.89 % -19.72 %
		2Y	4.64 %					4.02 %		31.12 %	-18.95 %
		3Y	56.24 %				4.04 %		32.10 %	-18.20 %
		5Y 7Y 10Y	E[( 28.58 % 84.92 % 14.54 %	0	T	y t dt) -p ] ≤ C 2.68 % 0.59 % 4.57 %	33.63 % e 2pκT T p , 35.41 % 39.98 %	-16.52 % -14.80 % -12.32 %	(4.6)
	available when λ ≤ 1. ( A -1 ǫ,t 0 y ǫ s ds)/(ǫξ Sup ) 2 and in view of (4.7), it follows that
	To sum up (4.2) and (4.6), we have proved that
			E[(	0	T	y t dt) -p ] ≤ Ce 2pκT (1 +	1 T p ),	(4.7)
	. t ≥ y ǫ Step 2. Take ǫ ∈]0, 1]. We apply Lemma 4.2 to v ǫ , in order to write v ǫ for a constant C depending only on p and v 0 . A -1 ǫ,t t = A ǫ,t 0 (ǫξ s ) 2 ds and dy ǫ t = ( 1 2 -κ (ǫξ In f ) 2 y ǫ t )dt + y ǫ t d Bǫ t , y ǫ 0 = y 0 . Thus, we get T 0 v ǫ (4.5) where t dt ≥

T 0 y t dt) -p ] ≤ C( e 2κT

λ 2 ) p . In view of the inequality (e x -1 ≥ x, x ≥ 0), we have λ 2 =

  T 3.16 Put prices of the closed formula, of the approximation formula and of the averaging formula, expressed as a function of maturities in fractions of years and relative strikes. Piecewise constant parameters.

	3M	30.00	20.07	10.78	3.93	0.87	0.13	0.05	0.02
		30.00	20.08	10.78	3.93	0.87	0.13	0.05	0.02
		30.00	20.07	10.78	3.93	0.87	0.13	0.05	0.02
	6M	40.01	30.06	20.41	5.53	2.06	0.18	0.05	0.01
		40.00	30.06	20.42	5.53	2.05	0.18	0.05	0.02
		40.01	30.06	20.41	5.53	2.06	0.18	0.05	0.01
	1Y	50.01	40.07	21.33	7.85	1.97	0.17	0.03	0.02
		50.01	40.07	21.35	7.84	1.95	0.18	0.04	0.02
		50.01	40.07	21.33	7.85	1.97	0.17	0.03	0.02
	2Y	60.02	50.11	31.38	11.23	2.92	0.24	0.06	0.01
		60.01	50.11	31.39	11.23	2.90	0.25	0.07	0.02
		60.02	50.11	31.38	11.23	2.92	0.24	0.06	0.01
	3Y	70.01	60.07	41.07	13.92	3.55	0.41	0.08	0.02
		70.01	60.07	41.08	13.92	3.54	0.42	0.09	0.03
		70.01	60.07	41.07	13.92	3.55	0.41	0.08	0.02
	5Y	80.01	70.07	42.64	18.37	5.74	0.61	0.15	0.04
		80.01	70.07	42.64	18.36	5.72	0.61	0.16	0.04
		80.01	70.07	42.64	18.37	5.74	0.61	0.15	0.04
	7Y	90.00	70.24	52.22	22.15	6.46	0.86	0.21	0.06
		90.00	70.24	52.22	22.15	6.45	0.86	0.21	0.07
		90.00	70.24	52.22	22.15	6.46	0.86	0.21	0.06
	10Y	90.01	80.14	54.13	27.17	8.71	1.42	0.35	0.11
		90.01	80.14	54.13	27.16	8.70	1.42	0.36	0.12
		90.01	80.14	54.13	27.17	8.71	1.42	0.35	0.11

  ǫ s ) -1 ≤ 1 for s ∈ [0, t]. Consequently, we deduce

	|R σ ǫ 0,t | ≤	0	t	ǫ 2 ξ 2 s 8σ 0,s	ds + (	0	.		ǫξ s 2	dB s ) * t + (	0	.	ǫξ s 2	dB s ) * t (1 -U ǫ t )
	≤	0	t	ǫ 2 ξ 2 s 8σ 0,s	ds + (	0	.	ǫξ s dB s ) * t .	(4.11)
	Now we easily complete the proof by observing that σ 0,s ≥ min( θ In f , ( . 0 ξ s dB s ) * √ t. t p ≤ Cξ Sup We define	√ v 0 ) and
							σ 1,t = U 0 t		0	t	(U 0 s ) -1 ξ s 2	dB s .
	Therefore, (σ 1,t ) 0≤t≤T solves the following SDE:
	dσ 1,t = -(	κθ t 2(σ 0,t ) 2 +	κ 2	)σ 1,t dt +	ξ t 2	dB t , σ 1,0 = 0,	(4.12)
	and for every p ≥ 1						(σ 1,. ) * t p ≤ Cξ Sup	√ t.	(4.13)
	P 4.7. Under (P), the process (R σ ǫ 1,t = σ ǫ t -σ 0 t -ǫσ 1,t ) 0≤t≤T fulfills the equality:
	R σ ǫ 1,t = U ǫ t	0	t	(U ǫ s ) -1 (-

  From Equation (2.3), one has D B s v 1,t = e -kt e ks ξ s √ v 0,s . Hence it is deterministic. Thus, using Identity (5.1) and Lemma 5.2 for the first equality and Equation (2.4) for the second equality, one has:

	Now, plug this last equality into (5.3) and use the identity
		0	T	((	0	t	ρ 2 s v 1,s ds)v 1,t dt + (	0	t	(1 -ρ 2 s )v 1,s ds)ρ 2 t v 1,t dt) +	1 2	(	0	T	(1 -ρ 2 t )v 1,t dt) 2 =
		0	T	((	0	t	ρ 2 s v 1,s ds)v 1,t dt + (	0	t	(1 -ρ 2 s )v 1,s ds)(ρ 2 t + 1 -ρ 2 t )v 1,t dt) =
		0	T	((	0	t	(ρ 2 s + 1 -ρ 2 s ))v 1,s ds)v 1,t dt =	1 2	(	0	T	v 1,t dt) 2 ;
	it immediately gives the result.
	5.3. Proof of Theorem 2.3. Proof. Step 1: We show the equality
										E[		∂ PBS ∂y	0	T	(v 1,t + v 2,t )dt] =	2 i=1	a i,T	∂ i+1 P BS (x 0 , ∂x i y 0 v 0,t dt) T	,
	where
															(κ,ρξv 0,. ),(-κ,1) a 1,T = ω 0,T
																	v 1,t dt
						+			0	T	(			0	t	(1 -ρ 2 s )v 1,s ds)(ρ t	v 1,t 2(v 0,t )	1 2	dB t -	ρ 2 t 2	v 1,t dt))]
					=	0	T	E[	∂ 2 ∂xy PBS	v 1,t (	0	t	ρ s	v 1,s 2(v 0,s )	1 2	dB s -	0	t	ρ 2 s 2	v 1,s ds)]dt
						+ E[	∂ 2 ∂y 2 PBS	0	T	(	0	t	(1 -ρ 2 s )v 1,s ds)ρ 2 t v 1,t dt]
					=	0	T	E[	∂ 2 ∂xy PBS	(v 1,t (-	0	t	ρ 2 s 2	v 1,s ds) +	0	t	ρ s	2	v 1,s √ v 0,s	D B s v 1,t ds)
						+	∂ 3 ∂x 2 y PBS	v 1,t	0	t	ρ 2 s 2	v 1,s ds]dt + E[	∂ 2 ∂y 2 PBS	0	T	(	0	t	(1 -ρ 2 s )v 1,s ds)ρ 2 t v 1,t dt].
	A =E[	∂ 2 ∂y 2 PBS			0	T	((	0	t	ρ 2 s v 1,s ds)v 1,t dt + (	0	t	(1 -ρ 2 s )v 1,s ds)ρ 2 t v 1,t dt)]
	+ E[	∂ PBS ∂y				0	T	(		0	t	v 1,s 2v 0,s	e -kt e ks ξ s	√ v 0,s dB s )dt]
	=E[	∂ 2 ∂y 2 PBS			0	T	((	0	t	ρ 2 s v 1,s ds)v 1,t dt + (	0	t	(1 -ρ 2 s )v 1,s ds)ρ 2 t v 1,t dt)] + E[	∂ PBS ∂y	0	T	v 2,t 2	dt].

  -κs ds)(e κt v 2 1,t dt + ξ t √ v 0,t e κt ( where we have used Lemma 5.3 ( f (t) = e -κt , Z t = ( t 0 v 1,s ds)(e κt v 1,t )) for the second equality and Lemmas 5.2 and 5.3 ( f (t) = ( T t e -κs ds)ρ t ξ t v 0,t e κt , Z t = t 0 v 1,s ds) for the last one. An application of the first and second equality in Lemma 5.4 gives the announced result. Actually, it remains to show that b 2,T = a 2 1,T /2. Indeed, consider two càdlàg functions f and g : [0, T] R. Then t 4 dt 4 dt 2 dt 3 dt 1 .

			(	T 0 f t (	T t g s ds)dt) 2 2	=	T 0	T 0 f t 1 (	T t 1	g t 3 dt 3 ) f t 2 ( 2	T t 2	g t 4 dt 4 )dt 2 dt 1
														T	T	T	T
										=		0	f t 1 (	t 1	t 1	g t 3 f t 2 (	t 2	g t 4 dt 4 )dt 3 dt 2 )dt 1
														T	T	T	T
										=		0	f t 1 (	t 1	f t 2	t 2	t 2	g t 3 g t 4 dt 3 dt 4 dt 2
														T	T	T
														+	t 1	g t 3	t 3	f t 2	t 2	g t 4 dt 4 dt 2 dt 3 )dt 1
														T	T	T	T
										=2			0	f t 1	t 1	f t 2	t 2	g t 3	t 3	g t 4 dt 3 dt 4 dt 2 dt 1
														T	T	T	T
														+	0	f t 1	t 1	g t 3	t 3	f t 2	t 2
														=	a 2 1,T 2	.
	Indeed, one has
	1 2	E[	∂ 2 ∂y 2 ( PBS	0	T	v 1,t dt) 2 ] = E[	∂ 2 ∂y 2 PBS	0	T	(	0	t	v 1,s ds)v 1,t dt]
				= E[	∂ 2 ∂y 2 PBS	0	T	(		t	T	0 e t	v 1,s ds)dB t )]
				= E[	∂ 2 ∂y 2 PBS	0	T	(		t	T	e -κs ds)e κt v 2 1,t dt] + E[	∂ 3 ∂xy 2 PBS	0	T	ω (κ,ρξv 0,. ),(-κ,1) t,T	v 1,t dt],

g

The approximation formula for the call price is obtained using the call/put parity relation: in (2.13), it consists of replacing on the l.h.s. the put payoff by the call one, and on the r.h.s., the put price function P BS by the similar call price function, while coefficients remain the same.

In this approach, we leave the initial value v0 equal to v 0 . Indeed, it is not natural to modify its value since it is not a parameter, but rather an unobserved factor.

Note that from the upper bound (4.11) in the proof, we easily obtain that the continuity also holds a.s., and not only in L p . Since only the latter is needed in what follows, we do not go into detail.