Junchao Chen
email: juchen@lpsm.paris

Noufel Frikha
email: frikha@lpsm.paris

Houzhi Li
email: houzhi.li@lpsm.paris

Probabilistic representation of integration by parts formulae for some stochastic volatility models with unbounded drift

In this paper, we establish a probabilistic representation as well as some integration by parts formulae for the marginal law at a given time maturity of some stochastic volatility model with unbounded drift. Relying on a perturbation technique for Markov semigroups, our formulae are based on a simple Markov chain evolving on a random time grid for which we develop a tailor-made Malliavin calculus. Among other applications, an unbiased Monte Carlo path simulation method stems from our formulas so that it can be used in order to numerically compute with optimal complexity option prices as well as their sensitivities with respect to the initial values or Greeks in finance, namely the Delta and Vega, for a large class of non-smooth European payoff. Numerical results are proposed to illustrate the efficiency of the method.

Introduction

In this work, we consider a two dimensional stochastic volatility model given by the solution of the following stochastic differential equation (SDE for short) with dynamics (1.1)

   S t = s 0 + t 0 rS s ds + t 0 σ S (Y s )S s dW s , Y t = y 0 + t 0 b Y (Y s ) ds + t 0 σ Y (Y s ) dB s , d B, W s = ρ ds
where the coefficients b Y , σ S , σ Y : R → R are smooth functions, W and B are one-dimensional standard Brownian motions with correlation factor ρ ∈ (-1, 1) both being defined on some probability space (Ω, F, P) .

The aim of this article is to prove a probabilistic representation formula for two integration by parts (IBP) formulae for the marginal law of the process (S, Y ) at a given time maturity T . To be more specific, for a given starting point (s 0 , y 0 ) ∈ (0, ∞) × R and a given finite time horizon T > 0, we establish two Bismut-Elworthy-Li (BEL) type formulae for the two following quantities (1.2) ∂ s0 E [h(S T , Y T )] and

∂ y0 E [h(S T , Y T )]
where h is a real-valued possibly non-smooth payoff function defined on [0, ∞) × R. Such IBP formulae have attracted a lot of interest during the last decades both from a theoretical and a practical point of views as they can be further analyzed to derive properties related to the transition density of the underlying process or to develop Monte Carlo simulation algorithm among other practical applications, see e.g. Nualart [START_REF] Nualart | The Malliavin calculus and related topics[END_REF], Malliavin and Thalmaier [START_REF] Malliavin | Stochastic Calculus of Variations in Mathematical Finance[END_REF] and the references therein. They are also of major interest for computing sensitivities, also referred as to Greeks in finance, of arbitrage price of financial derivatives which is the keystone for hedging purpose, i.e. for protecting the value of a portfolio against some possible changes in sources of risk. The two quantities appearing in (1.2) corresponds respectively to the Delta and Vega of the European option with payoff h(S T , Y T ). For a more detailed discussion on this topic, we refer the interested reader to Fournié and al. [START_REF] Fournié | Applications of Malliavin calculus to Monte Carlo methods in finance[END_REF], [START_REF] Fournié | Applications of Malliavin calculus to Monte-Carlo methods in finance[END_REF] for IBP formulae related to European, Asian options and conditional expectations, Gobet and al. [START_REF] Gobet | Computation of Greeks for barrier and look-back options using Malliavin calculus[END_REF], [START_REF] Bernis | Monte Carlo evaluation of Greeks for multidimensional barrier and lookback options[END_REF] for IBP formulae related to some barrier or lookback options. Let us importantly point out that, from a numerical point of view, the aforementioned IBP formulae will inevitably involve a time discretization procedure of the underlying process and Malliavin weights, thus introducing two sources of error given by a bias and a statistical error, as it is already the case for the computation of the price E[h(S T , Y T )].

Relying on a perturbation argument for the Markov semigroup generated by the couple (X, Y ), we first establish a probabilistic representation formula for the marginal law (S T , Y T ) for a fixed prescribed maturity T > 0 based on a simple Markov chain evolving along a random time grid given by the jump times of an independent renewal process. Such probabilistic representation formula was first derived in Bally and Kohatsu-Higa [START_REF] Bally | A probabilistic interpretation of the parametrix method[END_REF] for the marginal law of a multi-dimensional diffusion process and of some Lévy driven SDEs with bounded drift, diffusion and jump coefficients. Still in the case of bounded coefficients, it was then further investigated in Labordère and al. [START_REF] Henry-Labordère | Unbiased simulation of stochastic differential equations[END_REF], Agarwal and Gobet [START_REF] Agarwal | Finite variance unbiased estimation of stochastic differential equations[END_REF] for multi-dimensional diffusion processes and in Frikha and al. [START_REF] Frikha | Integration by parts formula for killed processes: a point of view from approximation theory[END_REF] for one-dimensional killed processes. The major advantage of the aforementioned probabilistic formulae lies in the fact that an unbiased Monte Carlo simulation method directly stems from it. Thus, it may be used to numerically compute an option price with optimal complexity since its computation will be only affected by the statistical error. However, let us emphasize that in general the variance of the Monte Carlo estimator tends to be large or even infinite. In order to circumvent this issue, an importance sampling scheme based on the law of the jump times of the underlying renewal process has been proposed in Anderson and Kohatsu-Higa [START_REF] Andersson | Unbiased simulation of stochastic differential equations using parametrix expansions[END_REF] in the multi-dimensional diffusion framework and in [START_REF] Frikha | Integration by parts formula for killed processes: a point of view from approximation theory[END_REF] for one-dimensional killed processes.

The main novelty of our approach in comparison with the aforementioned works is that we allow the drift coefficient b Y to be possibly unbounded as it is the case in most stochastic volatility models (Stein-Stein, Heston, ...). Such boundedness condition on the drift coefficient has appeared persistently in the previous contributions and is actually essential since basically it allows to remove the drift in the choice of the approximation process in order to derive the probabilistic representation formula. The key ingredient that we here develop in order to remove this restriction consists in choosing adequatly the approximation process around which the original perturbation argument of the Markov semigroup (X, Y ) is done by taking into account the transport of the initial condition by the deterministic ordinary differential equation (ODE) having unbounded coefficient 1 . The approximation process, or equivalently the underlying Markov chain on which the probabilistic representation is based, is then obtained from the original dynamics (1.1) by freezing the coefficients b Y , σ S and σ Y along the flow of this ODE. We stress that the previous choice is here crucial since it provides the adequate approximation process on which some good controls can be established. To the best of our knowledge, this feature appears to be new in this context.

Having this probabilistic representation formula at hand together with the tailor-made Malliavin calculus machinery for this well-chosen underlying Markov chain, in the spirit of the BEL formula established in [START_REF] Frikha | Integration by parts formula for killed processes: a point of view from approximation theory[END_REF] for killed diffusion processes with bounded drift coefficient, we rely on a propagation of the spatial derivatives forward in time then perform local IBP formulas on each time interval of the random time grid and finally merge them in a suitable manner in order to establish the two BEL formulae for the two quantities (1.2). Following the ideas developed in [START_REF] Andersson | Unbiased simulation of stochastic differential equations using parametrix expansions[END_REF], we achieve finite variance for the Monte Carlo estimators obtained from the probabilistic representation formulas of the couple (S T , Y T ) and of both IBP formulae by selecting adequatly the law of the jump times of the renewal process. We finally provide some numerical tests illustrating our previous analysis.

The article is organized as follows. In Section 2, we introduce our assumptions on the coefficients, present the approximation process that will be the main building block for our perturbation argument as well as the Markov chain that will play a central role in our probabilistic representation for the marginal law of the process (X, Y ) and for our IBP formulae. In addition, we construct the taillor-made Malliavin calculus machinery related to the underlying Markov chain upon which both IBP formulae are made. In Section 3, relying on the Markov chain introduced in Section 2, we establish in Theorem 3.1 the probabilistic representation formula for the coupled (S T , Y T ). In Section 4, we establish the BEL formulae for the two quantities appearing in (1.2). The main result of this section is Theorem 4.1. Some numerical results are presented in Section 5. The proofs of Theorem 3.1 and of some other technical but important results are postponed to the appendix of Section A.

Notations: For a fixed time T and positive integer n, we will use the following notation for time and space variables

s n = (s 1 , • • • , s n ), x n = (x 1 , • • • , x n ), the differentials ds n = ds 1 • • • ds n , dx n = dx 1 • • • dx n and also introduce the simplex ∆ n (T ) := {s n ∈ [0, T ] n : 0 ≤ s 1 < • • • s n ≤ T }.
In order to deal with time-degeneracy estimates, we will often use the following space-time inequality:

(1.3) ∀p, q > 0, ∀x ∈ R, |x| p e -q|x| 2 ≤ (p/(2qe)) p/2 .

For two positive real numbers α and β, we define the Mittag-Leffler function z → E α,β (z) = ∞ k=0 z k /Γ(αk + β). For a positive integer d, we denote by C ∞ p (R d ) the space of real-valued functions which are infinitely differentiable on R d with derivatives of any order having polynomial growth.

Preliminaries: assumptions, definition of the underlying Markov chain and related Malliavin calculus

2.1. Assumptions. Throughout the article, we work on a probability space (Ω, F, P) which is assumed to be rich enough to support all random variables that we will consider in what follows. We will work under the following assumptions on the coefficients:

(AR) The coefficients σ S and σ Y are bounded and smooth, in particular σ S and σ Y belong to C ∞ b (R). The drift coefficient b Y belongs to C ∞ (R) and admits bounded derivatives of any order greater than or equal to one. In particular, the drift coefficient b Y may be unbounded.

(ND) There exists κ ≥ 1 such that for all x ∈ R, κ -1 ≤ a S (x) ≤ κ, κ -1 ≤ a Y (x) ≤ κ where a S = σ 2 S and a Y = σ 2 Y . Therefore, without loss of generality, we will assume that both σ S and σ Y are positive function.

Apply Itô's Lemma to X t = ln(S t ). We get (2.1)

     X t = x 0 + t 0 r -1 2 a S (Y s ) ds + t 0 σ S (Y s ) dW s , Y t = y 0 + t 0 b Y (Y s ) ds + t 0 σ Y (Y s ) dB s , d B, W s = ρ ds,
with x 0 = ln(s 0 ). Without loss of generality, we will thus work with the Markov semigroup associated to the process (X, Y ), namely

P t h(x 0 , y 0 ) = E[h(X t , Y t )].
2.2. Choice of the approximation process. As already mentioned in the introduction, our strategy here is based on a probabilistic representation of the marginal law, in the spirit of the unbiased simulation method introduced for diffusion processes by Bally and Kohatsu-Higa [START_REF] Bally | A probabilistic interpretation of the parametrix method[END_REF], see also Labordère and al. [START_REF] Henry-Labordère | Unbiased simulation of stochastic differential equations[END_REF], and investigated from a numerical perspective by Andersson and Kohatsu-Higa [START_REF] Andersson | Unbiased simulation of stochastic differential equations using parametrix expansions[END_REF]. We also mention the recent contribution of one the author with Kohatsu-Higa and Li [START_REF] Frikha | Integration by parts formula for killed processes: a point of view from approximation theory[END_REF] for IBP formulae for the marginal law of one-dimensional killed diffusion processes.

However, at this stage, it is important to point out that our choice of approximation process significantly differs from the four aforementioned references. Indeed, in the previous contributions, the drift is assumed to be bounded and basically plays no role so that one usually removes it in the dynamics of the approximation process. In order to handle the unbounded drift term b Y appearing in the dynamics of the volatility process, one has to take into account the transport of the initial condition by the ODE obtained by removing the noise in the dynamics of Y . To be more specific, we denote by (m t (s, y)) t∈[s,T ] , 0 ≤ s ≤ T , the unique solution to the ODE ṁt = b Y (m t ) with initial condition m s = y. Observe that by time-homogeneity of the coefficient b Y , one has m t (s, y) = m t-s (0, y). We will simplify the notation when s = 0 and write m t (y 0 ) for m t (0, y 0 ). When there is no ambiguity, we will often omit the dependence with respect to the initial point y 0 and we only write m t for m t (y 0 ). We now introduce the approximation process ( X, Ȳ ) defined by

(2.2)    Xx0 t = x 0 + t 0 (r -1 2 a S (m s )) ds + t 0 σ S (m s ) dW s , Ȳ y0 t = y 0 + t 0 b Y (m s ) ds + t 0 σ Y (m s ) dB s , d B, W s = ρ ds.
We will make intensive use of the explicit form of the Markov semigroup ( Pt ) t∈[0,T ] defined for any bounded measurable map h : R 2 → R by Pt h(x 0 , y

0 ) = E[h( Xx0 t , Ȳ y0 t )]. Lemma 2.1. Let (x 0 , y 0 ) ∈ R 2 , ρ ∈ (-1, 1
) and t ∈ (0, ∞). Then, for any bounded and measurable map h :

R 2 → R, it holds Pt h(x 0 , y 0 ) = R 2 h(x, y) p(t, x 0 , y 0 , x, y) dxdy (2.3) with p(t, x 0 , y 0 , x, y) = 1 2πσ S,t σ Y,t 1 -ρ 2 t exp - 1 2 (x -x 0 -(rt -1 2 a S,t )) 2 a S,t (1 -ρ 2 t ) - 1 2 (y -m t ) 2 a Y,t (1 -ρ 2 t ) × exp ρ t (x -x 0 -(rt -1 2 a S,t ))(y -m t ) σ S,t σ Y,t (1 -ρ 2
t ) where we introduced the notations

a S,t = a S,t (y 0 ) := σ 2 S,t := t 0 a S (m s (y 0 )) ds, a Y,t = a Y,t (y 0 ) := σ 2 Y,t := t 0 a Y (m s (y 0 )) ds, σ S,Y,t = σ S,Y,t (y 0 ) := t 0 (σ S σ Y )(m s (y 0 )) ds, ρ t := ρσ S,Y,t /(σ S,t σ Y,t ).
Moreover, for any t ∈ (0, T ], there exists some positive constant C := C(T, ρ, a, r, κ) such that

(2.4) ∀t ∈ (0, T ], p(t, x 0 , y 0 , x, y) ≤ C q4κ (t, x 0 , y 0 , x, y)
where, for a positive parameter c, we introduced the density function

(2.5) (x, y) → qc (t, x 0 , y 0 , x, y) := 1 2πct exp - (x -x 0 ) 2 2ct - (y -m t ) 2 2ct .
Proof. We write

( Xx0 t , Ȳ y0 t ) = x 0 + rt - 1 2 a S,t + t 0 σ S (m s ) dW s , m t + t 0 σ Y (m s ) ρdW s + 1 -ρ 2 d W s
where W is a one-dimensional standard Brownian motion independent of W . We thus deduce that ( Xx0

t , Ȳ y0 t ) ∼ N (µ(t, x 0 , y 0 ), Σ t ) with µ(t, x 0 , y 0 ) = x 0 + rt - 1 2 a S,t , m t and Σ t = a S,t ρσ S,Y,t ρσ S,Y,t a Y,t .
The expression of the transition density then readily follows. Now, from (ND), it is readily seen that a S,t , a Y,t ≤ κt so that using the inequalities (a

-b) 2 ≥ 1 2 a 2 -b 2 and ρ 2 t ≤ ρ 2 , it follows p(t, x 0 , y 0 , x, y) = 1 2πσ S,t σ Y,t 1 -ρ 2 t exp - 1 2 (x -x 0 -(rt -1 2 a S,t )) 2 a S,t (1 -ρ 2 t ) - 1 2 (y -m t ) 2 a Y,t (1 -ρ 2 t ) × exp ρ t 1 -ρ 2 t (x -x 0 -(rt -1 2 a S,t ))(y -m t ) σ S,t σ Y,t ≤ C 1 2π(2κ)t exp -(4κ) -1 (x -x 0 ) 2 2t -(4κ) -1 (y -m t ) 2 2t =: C q4κ (t, x 0 , y 0 , x, y)
for some positive constants C := C(T, λ, ρ, a, r, κ).

We will also use the notation ( Xs,x t , Ȳ s,y t ) t≥s for the approximation process starting from (x, y) at time s and with coefficients frozen along the deterministic flow {m t (s, y) = m t-s (y), t ≥ s}. Note that the corresponding Markov semigroup satisfies Ps,t h(x, y) := E h( Xs,y t , Ȳ s,y t ) = E h( Xy t-s , Ȳ y t-s ) = Pt-s h(x, y).

2.3.

Markov chain on random time grid. The first tool that we will employ is a renewal process N that we now introduce.

Definition 2.1. Let τ := (τ n ) n≥1 be a sequence of random variables such that (τ n -τ n-1 ) n≥1 , with the convention τ 0 = 0, are i.i.d. with positive density function f and cumulant distribution function t → F (t) = t -∞ f (s) ds. Then, the renewal process N := (N t ) t≥0 with jump times τ is defined by N t := n≥1 1 {τn≤t} .

It is readily seen that, for any t > 0, {N t = n} = {τ n ≤ t < τ n+1 } and by an induction argument that we omit, one may prove that the joint distribution of (τ 1 , • • • , τ n ) is given by

P(τ 1 ∈ ds 1 , • • • , τ n ∈ ds n ) = n-1 j=0 f (s i+1 -s i )1 {0<s1<•••<sn} which in turn implies E[1 {Nt=n} Φ(τ 1 , • • • , τ n )] = E[1 {τn≤t<τn+1} Φ(τ 1 , • • • , τ n )] = ∞ t ∆n(t) Φ(s 1 , • • • , s n ) n j=0 f (s j+1 -s j ) ds n+1
with the convention s 0 = 0. Hence, by Fubini's theorem, it holds

E[1 {Nt=n} Φ(τ 1 , • • • , τ n )] = ∆n(t) Φ(s 1 , • • • , s n )(1 -F (t -s n )) n-1 j=0 f (s j+1 -s j ) ds n (2.6) for any measurable map Φ : ∆ n (t) → R satisfying E[1 {Nt=n} |Φ(τ 1 , • • • , τ n )|] < ∞.
Usual choices that we will consider are the followings.

Example 2.2.

(1) If the density function f is given by f (t) = λe -λt 1 [0,∞) (t) for some positive parameter λ, then N is a Poisson process with intensity λ.

(2) If the density function f is given by f

(t) = 1-α τ 1-α 1 t α 1 [0,τ ] (t)
for some parameters (α, τ ) ∈ (0, 1) × (0, ∞), then N is a renewal process with [0, τ ]-valued Beta(1 -α, 1) jump times.

(3) More generally, if the density function f is given by f

(t) = τ 1-α-β B(α,β) 1 t 1-α (τ -t) 1-β 1 [0,τ ] (t) for some parameters (α, β, τ ) ∈ (0, 1) 2 × (0, ∞), then N is a renewal process with [0, τ ]-valued Beta(α, β) jump times.
Given a sequence Z = (Z 1 n , Z 2 n ) n≥1 of i.i.d. random vector with law N (0, I 2 ) and a renewal process N independent of Z with jump times (τ i ) i≥0 , we set ζ i = τ i ∧ T , with the convention ζ 0 = 0, and we consider the two-dimensional Markov chain ( X, Ȳ ) with ( X0 , Ȳ0 ) = (x 0 , y 0 ) at time 0 (evolving on the random time grid (ζ i ) i≥0 ) and with dynamics for any 0

≤ i ≤ N T (2.7)    Xi+1 = Xi + r(ζ i+1 -ζ i ) -1 2 a S,i + σ S,i Z 1 i+1 , Ȳi+1 = m i + σ Y,i ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 ,
where we introduced the notations

a S,i := σ 2 S,i := a S,ζi+1-ζi ( Ȳi ) = ζi+1-ζi 0 a S (m s ( Ȳi )) ds, a Y,i := σ 2 Y,i := a Y,ζi+1-ζi ( Ȳi ) = ζi+1-ζi 0 a Y (m s ( Ȳi )) ds, σ S,Y,i := ζi+1-ζi 0 (σ S σ Y )(m s ( Ȳi )) ds, ρ i := ρ ζi+1-ζi ( Ȳi ) = ρ σ S,Y,i σ S,i σ Y,i , m i := m ζi+1-ζi ( Ȳi ).
We will denote by σ S,i the first derivative of y → σ S,i (y) taken at Ȳi and proceed similarly for the quantities σ Y,i , σ S,Y,i , ρ i and m i . We define the filtration G = (G i ) i≥0 where G i = σ(Z 1 j , Z 2 j , 0 ≤ j ≤ i), for i ≥ 1 and G 0 stands for the trivial σ-field. We assume that the filtration G satisfies the usual conditions. For an integer n, we will use the notations

ζ n = (ζ 0 , • • • , ζ n ) and τ n = (τ 0 , • • • , τ n ).
2.4. Tailor-made Malliavin calculus for the Markov chain ( X, Ȳ ). In this section we introduce a tailormade Malliavin calculus for the underlying Markov chain ( X, Ȳ ) defined by (2.7) which will be employed in order to establish our IBP formulae. Instead of using an infinite dimensional calculus as it is usually done in the literature, see e.g. Nualart [START_REF] Nualart | The Malliavin calculus and related topics[END_REF], the approach developed below is based on a finite dimensional calculus for which the dimension is given by the number of jumps of the underlying renewal process involved in the Markov chain ( X, Ȳ ). Definition 2.3. Let n ∈ N. For any i ∈ {0, • • • , n}, we define the set S i,n ( X, Ȳ ), as the space of random variables H such that

• H = h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 ), on the set {N T = n}, where we recall ζ n+1 := (ζ 0 , • • • , ζ n+1 ) = (0, ζ 1 , • • • , ζ n , T ). • For any s n+1 ∈ ∆ n+1 (T ), the map h(., ., ., ., s n+1 ) ∈ C ∞ p (R 4 ).
For a r.v. H ∈ S i,n ( X, Ȳ ), we will often abuse the notations and write

H ≡ H( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 )
that is the same symbol H may denote the r.v. or the function in the set S i,n ( X, Ȳ ). One can easily define the flow derivatives for H ∈ S i,n ( X, Ȳ ) as follows

∂ Xi+1 H = ∂ 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 ), ∂ Ȳi+1 H = ∂ 4 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 ), ∂ Xi H = ∂ 1 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 ) + ∂ 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 )∂ Xi Xi+1 , ∂ Ȳi H = ∂ 2 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 ) + ∂ 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 )∂ Ȳi Xi+1 + ∂ 4 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 )∂ Ȳi Ȳi+1 ,
and from the dynamics (2.7)

∂ Xi Xi+1 = 1, ∂ Ȳi Ȳi+1 = m i + σ Y,i ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 + σ Y,i ρ i 1 -ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 , (2.8) ∂ Ȳi Xi+1 = - 1 2 a S,i + σ S,i Z 1 i+1 = - 1 2 a S,i + σ S,i σ S,i Xi+1 -Xi -r(ζ i+1 -ζ i ) - 1 2 a S,i . (2.9)
We now define the integral and derivative operators for H ∈ S i,n ( X, Ȳ ), as

I (1) i+1 (H) = H Z 1 i+1 σ S,i (1 -ρ 2 i ) - ρ i 1 -ρ 2 i ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 σ S,i -D (1)
i+1 H, (2.10)

I (2) i+1 (H) = H ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 σ Y,i (1 -ρ 2 i ) - ρ i 1 -ρ 2 i Z 1 i+1 σ Y,i -D (2)
i+1 H, (2.11)

D (1) i+1 H = ∂ Xi+1 H, (2.12) D (2) i+1 H = ∂ Ȳi+1 H. (2.13)
Note that due to the above definitions and assumption (H), it is readily checked that

I (1) i+1 (H), I (2) i+1 (H), D (1)
i+1 H and D

(2) i+1 H are elements of S i,n ( X, Ȳ ) so that we can define iterations of the above operators. Namely, by induction, for a multi-index α = (α 1 , • • • , α p ) of length p with α i ∈ {1, 2} and α p+1 ∈ {1, 2}, we define

I (α,αp+1) i+1 (H) = I (αp+1) i+1 (I (α) i+1 (H)), D (α,αp+1) i+1 H = D (αp+1) i+1 (D (α) i+1 H) with the intuitive notation (α, α p+1 ) = (α 1 , • • • , α p+1 ).
Throughout the article, we will use the following notation for a certain type of conditional expectation that will be frequently employed. For any X ∈ L 1 (P) and any

i ∈ {0, • • • , n}, E i,n [X] = E[X|G i , τ n+1 , N T = n]
where we recall that we employ the notation τ n+1 = (τ 0 , • • • , τ n+1 ). Having the above definitions and notations at hand, the following duality formula is satisfied: for any non-empty multi-index α of length p, with α i ∈ {1, 2}, for any i ∈ {1, • • • , p}, p being a positive integer, it holds

(2.14) E i,n D (α) i+1 f ( Xi+1 , Ȳi+1 )H = E i,n f ( Xi+1 , Ȳi+1 )I (α)
i+1 (H) .

In order to obtain explicit norm estimates for random variables in S i,n ( X, Ȳ ), it is useful to define for

H ∈ S i,n ( X, Ȳ ), i ∈ {0, • • • , n} and p ≥ 1 H p p,i,n = E i,n [|H| p
]. We will also employ a chain rule formula for the integral operators defined above.

Lemma 2.2. Let H ≡ H( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 ) ∈ S i,n ( X, Ȳ ), for some i ∈ {0, • • • , n}. The following chain rule formulae hold for any (α 1 , α 2 ) ∈ {1, 2} 2 ∂ Xi I (α1) i+1 (H) = I (α1) i+1 (∂ Xi H), ∂ Xi I (α1,α2) i+1 (H) = I (α1,α2) i+1 (∂ Xi H). (2.15)
Moreover, one has

∂ Ȳi I (1) i+1 (H) = I (1) i+1 (∂ Ȳi H) - σ S,i σ S,i I (1)
i+1 (H) -

ρ i 1 -ρ 2 i σ Y,i σ S,i I (2)
i+1 (H), (2.16)

∂ Ȳi I (2) i+1 (H) = I (2) i+1 (∂ Ȳi H) - σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i I (2)
i+1 (H), (2.17)

∂ Ȳi I (1,1) i+1 (H) = I (1,1) i+1 (∂ Ȳi H) -2 σ S,i σ S,i I (1,1) i+1 (H) - ρ i 1 -ρ 2 i σ Y,i σ S,i I (1,2) i+1 (H) + I (2,1)
i+1 (H) , (2.18)

∂ Ȳi I (2,2) i+1 (H) = I (2,2) i+1 (∂ Ȳi H) -2 σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i I (2,2)
i+1 (H), (2.19)

∂ Ȳi I (1,2) i+1 (H) = I (1,2) i+1 (∂ Ȳi H) - σ S,i σ S,i + σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i I (1,2) i+1 (H) - ρ i 1 -ρ 2 i σ Y,i σ S,i I (2,2)
i+1 (H). (2.20)

Proof. Observe that from the very definitions (2.10) and (2.11), one directly gets

∂ Xi I (1) i+1 (1) = ∂ Xi I (2)
i+1 (1) = 0 while, also by direct computation, we obtain

∂ Ȳi I (1) i+1 (1) = - σ S,i σ S,i I (1)
i+1 (1) -

ρ i 1 -ρ 2 i σ Y,i σ S,i I (2)
i+1 (1),

∂ Ȳi I (2) i+1 (1) = - σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i I (2)
i+1 [START_REF] Agarwal | Finite variance unbiased estimation of stochastic differential equations[END_REF].

We thus deduce

∂ Xi I (α1) i+1 (H) = ∂ Xi HI (α1) i+1 (1) + H∂ Xi I (α1) i+1 (1) -∂ Xi D (α1) i+1 H = ∂ Xi HI (α1) i+1 (1) -D (α1) i+1 (∂ Xi H) = I (α1) i+1 (∂ Xi H)
where we used the fact

D (α1) i+1 ∂ Xi H = ∂ Xi D (α1)
i+1 H which easily follows by direct computation. As a consequence, it is readily seen

∂ Xi I (α1,α2) i+1 (H) = ∂ Xi I (α2) i+1 (I (α1) i+1 (H)) = I (α2) i+1 (∂ Xi I (α1) i+1 (H)) = I (α2) i+1 (I (α1) i+1 (∂ Xi H)) = I (α1,α2) i+1 (∂ Xi H).
This concludes the proof of (2.15). The chain rule formulae (2.16), (2.17), (2.18)

i+1 ∂ Ȳi H = ∂ Ȳi D (α1) i+1 H. Indeed, by standard computations, it holds ∂ Ȳi D (1) i+1 H = ∂ Ȳi ∂ Xi+1 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n ) = ∂ 2 2,3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n ) + ∂ 2 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Ȳi Xi+1 + ∂ 2 4,3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Ȳi Ȳi+1 , D (1)
i+1 ∂ Ȳi H = ∂ Xi+1 ∂ Ȳi h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n ) = ∂ Xi+1 (∂ 2 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n ) + ∂ 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Ȳi Xi+1 + ∂ 4 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Ȳi Ȳi+1 ) = ∂ 2 3,2 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n ) + ∂ 2 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Ȳi Xi+1 + ∂ 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Xi+1 ∂ Ȳi Xi+1 + ∂ 2 3,4 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Ȳi Ȳi+1 + ∂ 4 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Xi+1 ∂ Ȳi Ȳi+1 = ∂ Ȳi D (1) i+1 H + ∂ 3 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Xi+1 ∂ Ȳi Xi+1 + ∂ 4 h( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n )∂ Xi+1 ∂ Ȳi Ȳi+1 = ∂ Ȳi D (1) i+1 H + D (1) i+1 H∂ Xi+1 ∂ Ȳi Xi+1 + D (2) i+1 H∂ Xi+1 ∂ Ȳi Ȳi+1 = ∂ Ȳi D (1) i+1 H + σ S,i σ S,i D (1)
i+1 H + ρ i 1 -ρ 2 i σ Y,i σ S,i D (2)
i+1 H

where we used the two identities

∂ Xi+1 ∂ Ȳi Xi+1 = σ S,i σ S,i and ∂ Xi+1 ∂ Ȳi Ȳi+1 = ρ i 1-ρ 2 i σ Y,i
σ S,i which readily stems from (2.8), (2.9) and the dynamics (2.7).

From (2.10) and the previous identity, we thus obtain

∂ Ȳi I (1) i+1 (H) = ∂ Ȳi I (1) i+1 (1)H + I (1) i+1 (1)∂ Ȳi H -∂ Ȳi D (1) i+1 H = - σ S,i σ S,i I (1)
i+1 (1)H -

ρ i 1 -ρ 2 i σ Y,i σ S,i I (2)
i+1 (1)H + I

(1)

i+1 (1)∂ Ȳi H -D (1)
i+1 ∂ Ȳi H + σ S,i σ S,i D (1)
i+1 H

+ ρ i 1 -ρ 2 i σ Y,i σ S,i D (2) i+1 H = - σ S,i σ S,i I (1)
i+1 (1)H -D

(1)

i+1 H + I (1) i+1 (1)∂ Ȳi H -D (1) i+1 ∂ Ȳi H - ρ i 1 -ρ 2 i σ Y,i σ S,i I (2) i+1 (1)H -D (2) i+1 H = I (1) i+1 (∂ Ȳi H) - σ S,i σ S,i I (1)
i+1 (H) -

ρ i 1 -ρ 2 i σ Y,i σ S,i I (2)
i+1 (H).

Similarly, after some algebraic manipulations using (2.7) and (2.8), we get

∂ Ȳi+1 ∂ Ȳi Ȳi+1 = σ Y,i σ Y,i - ρ i ρi 1-ρ 2 i so that D (2) i+1 ∂ Ȳi H = ∂ Ȳi D (2) i+1 H + D (2) i+1 H∂ Ȳi+1 ∂ Ȳi Ȳi+1 = ∂ Ȳi D (2) i+1 H + σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i D (2)
i+1 H so that, omitting some technical details, we get

∂ Ȳi I (2) i+1 (H) = I (2) i+1 (∂ Ȳi H) - σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i I (2) i+1 (1)H + σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i D (2) i+1 H = I (2) i+1 (∂ Ȳi H) - σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i I (2)
i+1 (H).

The identities (2.18), (2.19) and (2.20) eventually follows from (2.16) and (2.17) using some simple algebraic computations.

We conclude this section by introducing the following space of random variables which satisfy some time regularity estimates.

Definition 2.4. Let ∈ Z and n ∈ N. For any i ∈ {0, • • • , n}, we define the space M i,n ( X, Ȳ , /2) as the set of finite random variables H ∈ S i,n ( X, Ȳ ) satisfying the following time regularity estimate: for any p ≥ 1, for any c > 0, there exists some positive constants C := C(T ), c , T → C(T ) being non-decreasing and c being independent of T , such that for any (x i , y i , x i+1 , y i+1 ,

s n+1 ) ∈ R 4 × ∆ n+1 (T ), |H(x i , y i , x i+1 , y i+1 , s n+1 )| p qc (s i+1 -s i , x i , y i , x i+1 , y i+1 ) (2.21) ≤ C(s i+1 -s i ) p /2 qc (s i+1 -s i , x i , y i , x i+1 , y i+1 )
where the density function R 2 (x i+1 , y i+1 ) → qc (s i+1 -s i , x i , y i , x i+1 , y i+1 ) is defined in Lemma 2.1.

We again remark that since the space M i,n ( X, Ȳ , /2) is a subset of S i,n ( X, Ȳ ), when we say that a random variable M i,n ( X, Ȳ , /2) this statement is always understood on the set {N T = n}.

Before proceeding, let us provide a simple example of some random variables that belong to the aforementioned space. From (2.10) and the dynamics (2.7) of the Markov chain ( X, Ȳ ), it holds

I (1) i+1 (1) = Xi+1 -Xi -(r(ζ i+1 -ζ i ) -1 2 a S,i ) a S,i (1 -ρ 2 i ) - ρ i 1 -ρ 2 i Ȳi+1 -m i σ S,i σ Y,i ,
I (1,1) i+1 (1) = (I (1) i+1 (1)) 2 -D 1 i+1 (I (1)
i+1 (1)) = (I

(1) i+1 (1)) 2 - 1 a S,i (1 -ρ 2 i )
, so that, I

i+1 (1) and I

(1,1)

i+1 (1) belong to S i,n ( X, Ȳ ). Moreover, under (ND), for any p ≥ 1, it holds

I (1) i+1 (1)(x i , y i , x i+1 , y i+1 , s n+1 ) p ≤ C 1 + |x i+1 -x i | p (s i+1 -s i ) p + |y i+1 -m i (y i )| p (s i+1 -s i ) p
and similarly,

I (1,1) i+1 (1)(x i , y i , x i+1 , y i+1 , s n+1 ) p ≤ C 1 + 1 (s i+1 -s i ) p + |x i+1 -x i | 2p (s i+1 -s i ) 2p + |y i+1 -m i (y i )| 2p (s i+1 -s i ) 2p .
Hence, from the space-time inequality (1.3), for any c > 0 and any c > c, it holds

I (1) i+1 (1)(x i , y i , x i+1 , y i+1 , s n+1 ) p qc (s i+1 -s i , x i , y i , x i+1 , y i+1 ) ≤ C(s i+1 -s i ) -p/2 qc (s i+1 -s i , x i , y i , x i+1 , y i+1 )
and

I (1,1) i+1 (1)(x i , y i , x i+1 , y i+1 , s n+1 ) p qc (s i+1 -s i , x i , y i , x i+1 , y i+1 ) ≤ C(s i+1 -s i ) -p qc (s i+1 -s i , x i , y i , x i+1 , y i+1 )
for some positive constant C := C(T ), T → C(T ) being non-decreasing. We thus conclude that I

(1)

i+1 (1) ∈ M i,n ( X, Ȳ , -1/2) and I (1,1) i+1 (1) ∈ M i,n ( X, Ȳ , -1) for any i ∈ {0, • • • , n}.
A straightforward generalization of the above example is the following property that will be frequently used in the sequel. We omit its proof.

Lemma 2.3. Let n ∈ N and i ∈ {0, • • • , n}. Assume that H ∈ M i,n ( X, Ȳ , /2) and D (α1) i H ∈ M i,n ( X, Ȳ , /2),
for some ( , ) ∈ Z 2 and some α 1 ∈ {1, 2}. Then, it holds

I (α1) i (H) ∈ M i,n ( X, Ȳ , (( -1) ∧ )/2). Additionally, if D (α2) i H ∈ M i,n ( X, Ȳ , /2) and D (α1,α2) i H ∈ M i,n ( X, Ȳ , /2), for some ∈ Z and α 2 ∈ {1, 2}, then it holds I (α1,α2) i (H) ∈ M i,n ( X, Ȳ , (( -2) ∧ ( -1) ∧ )/2). Finally, if H 1 ∈ M i,n ( X, Ȳ , 1 /2) and H 2 ∈ M i,n ( X, Ȳ , 2 /2) for some ( 1 , 2 ) ∈ Z 2 then H 1 H 2 ∈ M i,n ( X, Ȳ , ( 1 + 2 )/2) and (ζ i+1 -ζ i )I (α1) i (H 1 ) ∈ M i,n ( X, Ȳ , ( 1 + 1)/2).
Finally, we importantly emphasize that if H ∈ M i,n ( X, Ȳ , /2) for some n ∈ N, i ∈ {0, • • • , n} and ∈ Z, then, its conditional L p (P)-moment is finite and also satisfies a time regularity estimate. More precisely, for any p ≥ 1, it holds

(2.22) H p,i,n ≤ C(ζ i+1 -ζ i ) /2
for some positive constant C := C(T ), T → C(T ) being non-decreasing. Indeed, using the fact that the sequence Z is independent of N as well as the upper-estimate (2.4) of Lemma 2.1 and finally (2.21), one directly gets

H p p,i,n = E i,n |H( Xi , Ȳi , Xi+1 , Ȳi+1 , ζ n+1 )| p Xi , Ȳi , τ n+1 , N T = n = R 2 |H( Xi , Ȳi , x i+1 , y i+1 , ζ n+1 )| p p(ζ i+1 -ζ i , Xi , Ȳi , x i+1 , y i+1 ) dx i+1 dy i+1 ≤ C R 2 |H( Xi , Ȳi , x i+1 , y i+1 , ζ n+1 )| p q4κ (ζ i+1 -ζ i , Xi , Ȳi , x i+1 , y i+1 ) dx i+1 dy i+1 ≤ C(ζ i+1 -ζ i ) p /2
so that (2.22) directly follows. The previous conditional L p (P)-moment estimate will be used at several places in the sequel.

Probabilistic representation for the couple (S T , Y T ).

In this section, we establish a probabilistic representation for the marginal law (S T , Y T ), or equivalently, for the law of (X T , Y T ) which is based on the Markov chain ( X, Ȳ ) introduced in the previous section. For γ > 0, we denote by B γ (R 2 ) the set of Borel measurable map h : R 2 → R satisfying the following exponential growth assumption at infinity, namely, for some positive constant C, for any (x,

y) ∈ R 2 , (3.1) |h(x, y)| ≤ C exp(γ(|x| 2 + |y| 2 )).
Theorem 3.1. Let T > 0. Under assumptions (AR) and (ND), the law of the couple (X T , Y T ) given by the unique solution to the SDE (2.1) at time T starting from (x 0 = ln(s 0 ), y 0 ) at time 0 satisfies the following probabilistic representation: there exists a positive constant c := c(T, b Y , κ) such that for any 0 ≤ γ < c -1 and any h ∈ B γ (R 2 ), it holds

(3.2) E[h(X T , Y T )] = E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i
where the random variables θ i ∈ S i-1,n ( X, Ȳ ) are defined by

θ i = (f (ζ i -ζ i-1 )) -1 I (1,1) i (c i S ) -I (1)
i

(c i S ) + I (2,2) i (c i Y ) + I (2) i (b i Y ) + I (1,2) i (c i Y,S ) , 1 ≤ i ≤ N T , (3.3) θ N T +1 = (1 -F (T -ζ N T )) -1 , (3.4) with c i S := 1 2 a S ( Ȳi ) -a S (m i-1 ) , c i Y := 1 2 a Y ( Ȳi ) -a Y (m i-1 ) , b i Y := b Y ( Ȳi ) -b Y (m i-1 ), c i Y,S := ρ((σ S σ Y )( Ȳi ) -(σ S σ Y )(m i-1 )).
Assume furthermore that N is a renewal process with Beta(α, 1) jump times. For any p ≥ 1 satisfying p( 1 2 -α) ≤ 1 -α, for any γ such that 0 ≤ pγ < c -1 and any h ∈ B γ (R 2 ), the random variable appearing inside the expectation in the right-hand side of (3.2) admits a finite L p (P)-moment. In particular, if α = 1/2 then for any p ≥ 1, for any h ∈ B γ (R 2 ) with 0 ≤ pγ < c -1 , the L p (P)-moment is finite.

The proof of Theorem 3.1 is postponed to Appendix A.1.

Integration by parts formulae

In this section, we establish two IBP formulae for the law of the couple (S T , Y T ). More precisely, we are interested in providing a Bismut-Elworthy-Li formula for the two quantities

∂ s0 E[h(S T , Y T )], ∂ y0 E[h(S T , Y T )].
Our strategy is divided into two steps as follows:

Step 1: The first step was performed with the probabilistic representation established in Theorem 3.1 for the couple (X T , Y T ) involving the two-dimensional Markov chain ( X, Ȳ ) evolving on a time grid governed by the jump times of the renewal process N . Introducing h(x, y) = f (e x , y) and assuming that f is of polynomial growth at infinity, it is sufficient to consider the two quantities

∂ s0 E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i , ∂ y0 E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i for h ∈ B γ (R 2 )
for some γ > 0 recalling that x 0 = ln(s 0 ).

Step 2: At this stage, one might be tempted to perform a standard IBP formula as presented in Nualart [START_REF] Nualart | The Malliavin calculus and related topics[END_REF] on the whole time interval [0, T ]. However, such a strategy is likely to fail. The main reason is that the Skorokhod integral of the product of weights N T i=1 θ i will inevitably involve the Malliavin derivative of θ i which will in turn raise some integrability issues of the resulting Malliavin weight. The key idea that we use in order to circumvent this issue consists in performing local IBP formulae on each of the random intervals

[ζ i , ζ i+1 ], i = 0, • • • , N T ,
that is, by using the noise of the Markov chain on this specific time interval and then by combining all these local IBP formulae in a suitable way.

To implement successfully our strategy, two main ingredients are needed. Our first ingredient consists in transferring the partial derivatives ∂ s0 and ∂ y0 on the expectation forward in time from the first time interval [0, ζ 1 ] to the interval on which we perform the local IBP formula, say [ζ i , ζ i+1 ]. Our second ingredient consists in combining these various local IBP formulae in an adequate manner. Roughly speaking, we will consider a weighted sum of each IBP formula, the weight being precisely the length of the corresponding time interval.

The transfer of derivative formula

. Lemma 4.1. Let h ∈ C 1 p (R 2 ) and n ∈ N. The maps R 2 (x, y) → E i,n h( Xi+1 , Ȳi+1 )θ i+1 |( Xi , Ȳi ) = (x, y) , i ∈ {0, • • • , n}, belong to C 1 p (R 2 ) a.s.
Moreover, the following transfer of derivative formulae hold

∂ s0 E 0,n h( X1 , Ȳ1 )θ 1 = E 0,n ∂ X1 h( X1 , Ȳ1 ) θ 1 s 0 (4.1) while for 1 ≤ i ≤ n, ∂ Xi E i,n h( Xi+1 , Ȳi+1 )θ i+1 = E i,n ∂ Xi+1 h( Xi+1 , Ȳi+1 )θ i+1 . (4.2)
Similarly, the following transfer of derivative formulae hold: for any 0 ≤ i ≤ n -1

∂ Ȳi E i,n h( Xi+1 , Ȳi+1 )θ i+1 = E i,n ∂ Ȳi+1 h( Xi+1 , Ȳi+1 ) - → θ e,Y i+1 + E i,n ∂ Xi+1 h( Xi+1 , Ȳi+1 ) - → θ e,X i+1 + E i,n h( Xi+1 , Ȳi+1 ) - → θ c i+1 (4.3) with - → θ e,Y i+1 = (f (ζ i+1 -ζ i )) -1 I (1,1) i+1 (d i+1 S ) + I (2,2) i+1 (d i+1 Y ) + I (1) i+1 (e Y,i+1 S ) + I (2) i+1 (e Y,i+1 Y ) + I (1,2) i+1 (d i+1 Y,S ) , - → θ e,X i+1 = (f (ζ i+1 -ζ i )) -1 I (1) i+1 (e X,i+1 S ), - → θ c i+1 = I (1) i+1 ∂ Ȳi Xi+1 θ i+1 - - → θ e,X i+1 + ∂ Ȳi θ i+1 + I (2) i+1 m i θ i+1 - - → θ e,Y i+1 + σ Y,i ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 + σ Y,i ρ i 1 -ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 θ i+1 , d i+1 S = m i c i+1 S , d i+1 Y = m i c i+1 Y , d i+1 Y,S = m i c i+1 Y,S , e Y,i+1 S = -m i c i+1 S + ∂ Ȳi c i+1 Y,S , e Y,i+1 Y = m i b i+1 Y + ∂ Ȳi c i+1 Y , e X,i+1 S = ∂ Ȳi c i+1 S .
For i = n, one also has

∂ Ȳn E n,n h( Xn+1 , Ȳn+1 )θ n+1 = E n,n ∂ Ȳn+1 h( Xn+1 , Ȳn+1 ) - → θ e,Y n+1 + E n,n ∂ Xn+1 h( Xn+1 , Ȳn+1 ) - → θ e,X n+1 + E n,n h( Xn+1 , Ȳn+1 ) - → θ c n+1 (4.4)
with

- → θ e,Y n+1 = (1 -F (T -ζ n )) -1 m n + σ Y,n ρ n Z 1 n+1 + 1 -ρ 2 n Z 2 n+1 + σ Y,n ρ n 1 -ρ 2 n 1 -ρ 2 n Z 1 n+1 -ρ n Z 2 n+1 , - → θ e,X n+1 = (1 -F (T -ζ n )) -1 - 1 2 a S,n + σ S,n Z 1 n+1 ,
and we set -→ θ c n+1 = 0 for notational convenience. Finally, the weight sequences (

- → θ e,Y i ) 1≤i≤n+1 , ( - → θ e,X i ) 1≤i≤n+1 and ( - → θ c i ) 1≤i≤n+1 defined above satisfy f (ζ i -ζ i-1 ) - → θ e,Y i , f (ζ i -ζ i-1 ) - → θ c i ∈ M i-1,n ( X, Ȳ , -1/2), f (ζ i -ζ i-1 ) - → θ e,X i ∈ M i-1,n ( X, Ȳ , 0), i ∈ {1, • • • , n} and (1 -F (T -ζ n )) - → θ e,Y n+1 ∈ M n,n ( X, Ȳ , 0), (1 -F (T -ζ n )) - → θ e,X n+1 , ∈ M n,n ( X, Ȳ , 1/2).
The proof of Lemma 4.1 is postponed to Appendix A.2. The transfer of derivative procedure starts on the first time interval [0, ζ 1 ] according to formulae (4.1) and (4.3) (for i = 0). It expresses the fact that the flow derivatives ∂ s0 and ∂ y0 of the conditional expectations on the left-hand side of the equations are transferred to derivative operators ∂ X1 and ∂ Ȳ1 on the test function h appearing on the right-hand side. Remark that the first derivatives of h have been written ubiquitously as ∂ Xi+1 h( Xi+1 , Ȳi+1 ) and ∂ Ȳi+1 h( Xi+1 , Ȳi+1 ).

Then, by the Markov property satisfied by the process ( X, Ȳ ), the function h appearing inside the (conditional) expectations on the right-hand side of (4.1) and (4.3) (for i = 0) will be given by the conditional expectation appearing on the left-hand side of the same equations but for i = 1. The transfer of derivative formulae for the following time intervals are obtained by induction using (4.2) and (4.3) up to the last time interval. Doing so, we obtain various transfer of derivative formulae by transferring successively the derivative operators through all intervals forward in time.

4.2. The integration by parts formulae. We first define the weights that will be used in our IBP formulae. For an integer n, on the set {N T = n}, for any k ∈ {1, • • • , n + 1} and any j ∈ {1, • • • , k}, we define

- → θ I (1),n+1 k := n+1 i=k+1 θ i × I (1) k (θ k ) × k-1 i=1 θ i , - → θ C n+1 j := n+1 i=j+1 θ i × - → θ c j × j-1 i=1 - → θ e,Y i , - → θ I (2),n+1 k := n+1 i=k+1 θ i × I (2) k ( - → θ e,Y k ) × k-1 i=1 - → θ e,Y i , - → θ I (1),n+1 k j := n+1 i=k+1 θ i × I (1) k (θ k ) × k-1 i=j+1 θ i × - → θ e,X j × j-1 i=1 - → θ e,Y i , j = 1, • • • , k -1, - → θ I (1),n+1 k k := n+1 i=k+1 θ i × I (1) k ( - → θ e,X k ) × k-1 i=1 - → θ e,Y i .
With the above definitions at hand, we are now able to state our IBP formulae.

Theorem 4.1. Let T > 0. Under assumptions (AR) and (ND), the law of the couple (X T , Y T ), given by the unique solution to the SDE (2.1) at time T starting from (x 0 = ln(s 0 ), y 0 ) at time 0, satisfies the following Bismut-Elworthy-Li IBP formulae: there exists some positive constant c := c(T, b Y , κ) such that for any 0 ≤ γ < c -1 and any h ∈ B γ (R 2 ), for any (s 0 , y 0 ) ∈ (0, ∞) × R, it holds

s 0 T ∂ s0 E h(X T , Y T ) = E h( XN T +1 , ȲN T +1 ) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (1),N T +1 k (4.5)
and

T ∂ y0 E h(X T , Y T ) = E   h( XN T +1 , ȲN T +1 ) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (2),N T +1 k + k j=1 - → θ C N T +1 j + - → θ I (1),N T +1 k j   . (4.6)
Moreover, if N is a renewal process with Beta(α, 1) jump times, then, for any p ≥ 1 satisfying p( 1 2 -α) ≤ 1 -α, for any γ such that 0 ≤ pγ < c -1 and any h ∈ B γ (R 2 ), the random variables appearing inside the expectation in the right-hand side of (4.5) and (4.6) admit a finite L p (P)-moment. In particular, if α = 1/2 then for any p ≥ 1, for any h ∈ B γ (R 2 ) with 0 ≤ pγ < c -1 , the L p (P)-moment is finite.

Proof. We only prove the IBP formula (4.6). The proof of (4.5) follows by completely analogous (and actually more simple) arguments and is thus omitted.

Step 1: proof of the IBP formula (4.6)

for h ∈ C 1 b (R 2 ). Let h ∈ C 1 b (R 2 )
. From Theorem 3.1 and Fubini's theorem, we write

E[h(X T , Y T )] = n≥0 E E h( Xn+1 , Ȳn+1 ) n+1 i=1 θ i |τ n+1 1 {N T =n} (4.7)
where we used the fact that {N T = n} = {τ n+1 > T } ∩ {τ n ≤ T }. In most of the arguments below, we will work on the set {N T = n}. In order to perform our induction argument forward in time through the Markov chain structure, we define for k ∈ {0, • • • , n + 1} the functions

H k ( Xk , Ȳk ) := E k,n h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i = E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i | Xk , Ȳk , τ n+1 , N T = n .
We also let H n+1 ( Xn+1 , Ȳn+1 ) := h( Xn+1 , Ȳn+1 ). Note that we omit the dependence w.r.t the sequence τ n+1 in the definition of the (random) maps (H k ) 0≤k≤n+1 . From the above definition and using (ND), (AR), it follows that the map H k belongs to C 1 p (R 2 ) a.s. for any 0 ≤ k ≤ n + 1. Moreover, from the tower property of conditional expectation the following relation is satisfied for any

k ∈ {0, • • • , n} H k ( Xk , Ȳk ) = E k,n [H k+1 ( Xk+1 , Ȳk+1 )θ k+1 ]. (4.8)
Now, using first the Lebesgue differentiation theorem and then iterating the transfer of derivative formula (4.3) in Lemma 4.1, we obtain2 for any k

∈ {1, • • • , n}, ∂ y0 E h( Xn+1 , Ȳn+1 ) n+1 i=1 θ i τ n+1 = ∂ y0 E E 0,n H 1 ( X1 , Ȳ1 )θ 1 τ n+1 = E ∂ y0 E 0,n H 1 ( X1 , Ȳ1 )θ 1 τ n+1 = E D (2) k H k ( Xk , Ȳk ) k i=1 - → θ e,Y i τ n+1 + k j=1 E H j ( Xj , Ȳj ) - → θ c j j-1 i=1 - → θ e,Y i τ n+1 + k j=1 E D (1) j H j ( Xj , Ȳj ) - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 . (4.9)
To further simplify the first term appearing on the right-hand side of (4.9), we use the tower property of conditional expectation (w.r.t E k-1,n [.]) and the integration by parts formula (2.14). For any k ∈ {1, • • • , n}, we obtain

E D (2) k H k ( Xk , Ȳk ) - → θ e,Y k G k-1 , τ n+1 = E H k ( Xk , Ȳk )I (2) k ( - → θ e,Y k ) G k-1 , τ n+1 .
We also simplify the third term appearing on the right-hand side of (4.9), by using the transfer of derivatives formula (4.2) up to the time interval

[ζ k-1 , ζ k ]. For any j ∈ {1, • • • , k}, it holds E D (1) j H j ( Xj , Ȳj ) - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 = E D (1) k H k ( Xk , Ȳk ) k i=j+1 θ i - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 , so that, if j ∈ {1, • • • , k -1}, taking conditional expectation (using again E k-1,n [.]
) and then performing an IBP formula on the last time interval

[ζ k-1 , ζ k ] yield E D (1) k H k ( Xk , Ȳk ) k i=j+1 θ i - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 = E H k ( Xk , Ȳk )I (1) k (θ k ) k-1 i=j+1 θ i - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 , while if j = k, we obtain E D (1) k H k ( Xk , Ȳk ) k i=j+1 θ i - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 = E H k ( Xk , Ȳk )I (1) k ( - → θ e,X k ) k-1 i=1 - → θ e,Y i τ n+1 .
Coming back to (4.9) and using the definition of the maps (H k ) 0≤k≤n+1 , we thus deduce

∂ y0 E h( Xn+1 , Ȳn+1 ) n+1 i=1 θ i τ n+1 = E H k ( Xk , Ȳk )I (2) k ( - → θ e,Y k ) × k-1 i=1 - → θ e,Y i τ n+1 + k j=1 E H j ( Xj , Ȳj ) - → θ c j × j-1 i=1 - → θ e,Y i τ n+1 + k-1 j=1 E H k ( Xk , Ȳk )I (1) k (θ k ) × k-1 i=j+1 θ i × - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 + E H k ( Xk , Ȳk )I (1) k ( - → θ e,X k ) k-1 i=1 - → θ e,Y i τ n+1 = E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (2) k ( - → θ e,Y k ) × k-1 i=1 - → θ e,Y i τ n+1 + k j=1 E h( Xn+1 , Ȳn+1 ) n+1 i=j+1 θ i × - → θ c j × j-1 i=1 - → θ e,Y i τ n+1 + k-1 j=1 E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (1) k (θ k ) × k-1 i=j+1 θ i - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 + E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (1) k ( - → θ e,X k ) k-1 i=1 - → θ e,Y i τ n+1 . (4.10)
In the case k = n + 1, using the transfer of derivative formula (4.4) of Lemma 4.1 on the last time interval and then performing the IBP formula (2.14), we obtain the representation

∂ y0 E h( Xn+1 , Ȳn+1 ) n+1 i=1 θ i τ n+1 = E D (2) n+1 h( Xn+1 , Ȳn+1 ) n+1 i=1 - → θ e,Y i τ n+1 + n+1 j=1 E H j ( Xj , Ȳj ) - → θ c j j-1 i=1 - → θ e,Y i τ n+1 + n+1 j=1 E D (1) j H j ( Xj , Ȳj ) - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 = E h( Xn+1 , Ȳn+1 )I (2) n+1 ( - → θ e,Y n+1 ) n i=1 - → θ e,Y i τ n+1 + n+1 j=1 E H j ( Xj , Ȳj ) - → θ c j j-1 i=1 - → θ e,Y i τ n+1 + n+1 j=1 E D (1) n+1 h( Xn+1 , Ȳn+1 ) n+1 i=j+1 θ i - → θ e,X j j-1 i=1 - → θ e,Y i τ n+1 = E h( Xn+1 , Ȳn+1 )I (2) n+1 ( - → θ e,Y n+1 ) × n i=1 - → θ e,Y i τ n+1 + n+1 j=1 E h( Xn+1 , Ȳn+1 ) n+1 i=j+1 θ i × - → θ c j × j-1 i=1 - → θ e,Y i τ n+1 + n j=1 E h( Xn+1 , Ȳn+1 )I (1) n+1 (θ n+1 ) × n i=j+1 θ i × - → θ e,X j × j-1 i=1 - → θ e,Y i τ n+1 + E h( Xn+1 , Ȳn+1 )I (1) n+1 ( - → θ e,X n+1 ) × n i=1 - → θ e,Y i τ n+1 (4.11)
where, for the last term appearing in the right-hand side of the above identities, we employed the transfer of derivative formula (4.2) up to the last time interval and then performed an IBP formula. Now, the key point in order to establish the IBP formula (4.6) is to combine in a suitable way the identities (4.10) and (4.11). For each k ∈ {0, • • • , n}, we multiply the above formulae by the length of the interval on which the local IBP formula is performed, namely we multiply by ζ k -ζ k-1 both sides of (4.10), k = 1, • • • , n -1, and we multiply by T -ζ n both sides of (4.11). We then sum them over all k. Recalling that

n+1 k=1 ζ k -ζ k-1 = T -ζ 0 = T , we deduce T ∂ y0 E h( Xn+1 , Ȳn+1 ) n+1 i=1 θ i τ n+1 = n+1 k=1 (ζ k -ζ k-1 )E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (2) k ( - → θ e,Y k ) × k-1 i=1 - → θ e,Y i τ n+1 + n+1 k=1 (ζ k -ζ k-1 ) k j=1 E h( Xn+1 , Ȳn+1 ) n+1 i=j+1 θ i × - → θ c j × j-1 i=1 - → θ e,Y i τ n+1 + n+1 k=1 (ζ k -ζ k-1 ) k-1 j=1 E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (1) k (θ k ) × k-1 i=j+1 θ i × - → θ e,X j × j-1 i=1 - → θ e,Y i τ n+1 + E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (1) k ( - → θ e,X k ) × n i=1 - → θ e,Y i τ n+1 = E h( Xn+1 , Ȳn+1 ) n+1 k=1 (ζ k -ζ k-1 ) - → θ I (2),n+1 k + k j=1 - → θ C n+1 j + - → θ I (1),n+1 k j τ n+1 .
We now provide a sharp upper-estimate for the above quantity. From Lemma B.2 and Lemma 4.1, it follows that

f (ζ i -ζ i-1 )θ i , f (ζ i -ζ i-1 ) - → θ e,Y i , f (ζ i -ζ i-1 ) - → θ c i ∈ M i-1,n ( X, Ȳ , -1/2) and f (ζ i -ζ i-1 ) - → θ e,X i ∈ M i-1,n ( X, Ȳ , 0) for any i ∈ {1, • • • , n}.
Moreover, from the very definition of the weights θ i , -→ θ e,X i and -→ θ e,Y i , after some simple but cumbersome computations that we omit, one has f

(ζ i -ζ i-1 )D (1) i (θ i ), f (ζ i -ζ i-1 )D (2) i ( - → θ e,Y i ) ∈ M i-1,n ( X, Ȳ , -1) and f (ζ i -ζ i-1 )D (1) i ( - → θ e,X i ) ∈ M i-1,n ( X, Ȳ , -1/2) so that from Lemma 2.3 we conclude f (ζ i -ζ i-1 )(ζ i -ζ i-1 )I (1) i (θ i ) ∈ M i-1,n ( X, Ȳ , 0), f (ζ i -ζ i-1 )(ζ i -ζ i-1 )I (2) i ( - → θ e,Y i ) ∈ M i-1,n ( X, Ȳ , 0) and f (ζ i -ζ i-1 )(ζ i -ζ i-1 )I (1) i ( -
→ θ e,X i ) ∈ M i-1,n ( X, Ȳ , 1/2).
Hence, from the boundedness of h, the tower property of conditional expectation and (2.22), it holds

(ζ k -ζ k-1 )E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (2) k ( - → θ e,Y k ) × k-1 i=1 - → θ e,Y i τ n+1 ≤ C n+1 (1 -F (T -ζ n )) -1 n i=k+1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1 2 (f (ζ k -ζ k-1 )) -1 k-1 i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1 2
so that using the identity (2.6)

n≥0 E n+1 k=1 (ζ k -ζ k-1 )E h( Xn+1 , Ȳn+1 ) n+1 i=k+1 θ i × I (2) k ( - → θ e,Y k ) × k-1 i=1 - → θ e,Y i τ n+1 1 {N T =n} ≤ n≥0 C n+1 n+1 k=1 E (1 -F (T -ζ n )) -1 (f (ζ k -ζ k-1 )) -1 n i=1,i =k (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1 2 1 {N T =n} ≤ n≥0 C n+1 n+1 k=1 ∆n(T ) n i=1,i =k (s i -s i-1 ) -1/2 ds n ≤ n≥0 (n + 1)C n+1 T (n+1)/2 Γ n (1/2) Γ(1 + n/2) < ∞.
From similar arguments that we omit, it follows

(ζ k -ζ k-1 ) k j=1 E h( Xn+1 , Ȳn+1 ) - → θ C n+1 j + - → θ I (1),n+1 k j τ n+1 ≤ C n+1 (ζ k -ζ k-1 ) k j=1 (1 -F (T -ζ n )) -1 n i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1/2 [1 + 1 {i=k} (ζ i -ζ i-1 ) -1/2 ].
so that using again the identity (2.6)

n≥0 E n+1 k=1 (ζ k -ζ k-1 ) k j=1 E h( Xn+1 , Ȳn+1 ) - → θ C n+1 j + - → θ I (1),n+1 k j τ n+1 1 {N T =n} ≤ n≥0 C n+1 n+1 k=1 E (ζ k -ζ k-1 ) k j=1 (1 -F (T -ζ n )) -1 n i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1/2 [1 + 1 {i=k} (ζ i -ζ i-1 ) -1/2 ]1 {N T =n} ≤ n≥0 C n+1 (n + 1)(n + 2)T (n+1)/2 Γ n (1/2) Γ(1 + n/2) < ∞.
The preceding estimates combined with (4.7) and the Lebesgue dominated convergence theorem allows to conclude that

y 0 → E[h(X T , Y T )] is continuously differentiable with T ∂ y0 E[h(X T , Y T )] = T ∂ y0 E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i = n≥0 E T ∂ y0 E h( Xn+1 , Ȳn+1 ) n+1 i=1 θ i τ n+1 1 {N T =n} = n≥0 E E h( Xn+1 , Ȳn+1 ) n+1 k=1 (ζ k -ζ k-1 ) - → θ I (2),n+1 k + k j=1 - → θ C n+1 j + - → θ I (1),n+1 k j τ n+1 1 {N T =n} = E h( XN T +1 , ȲN T +1 ) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (2),N T +1 k + k j=1 - → θ C n+1 j + - → θ I (1),N T +1 k j
where we used Fubini's theorem for the last equality. This completes the proof of the IBP formula (4.6) for

h ∈ C 1 b (R 2 ).
Step 2: Extension to h ∈ B γ (R 2 ) for some positive γ.

We now extend the two IBP formulae that we have established in the previous step to the case of a test function h ∈ B γ (R 2 ) for some sufficiently small γ > 0. Let us note that under assumption (H), from Kusuoka and Stroock [START_REF] Kusuoka | Applications of the Malliavin calculus[END_REF], Corollary (3.25) and the upper-estimate (3.27) therein, the process (X t , Y t ) t≥0 admits a smooth transition density (t, x 0 , y 0 , x, y) → p(t, x 0 , y 0 , x, y)

∈ C ∞ ((0, ∞) × R 2 × R 2 ) and for any h ∈ C 1 b (R 2 ), it holds ∂ α s0 ∂ β y0 E[h(X T , Y T )] = R 2 h(x, y) ∂ α s0 ∂ β y0 p(T, x 0 , y 0 , x, y) dxdy
for any T > 0 and any integers α and β.

We then proceed as in step 2 of the proof of Theorem 3.1. Namely, we prove that

T ∂ y0 E[h(X T , Y T )] = E h( XN T +1 , ȲN T +1 ) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (2),N T +1 k + k j=1 - → θ C n+1 j + - → θ I (1),N T +1 k j = R 2 h(x, y) E p(T -ζ N T , XN T , ȲN T , x, y) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (2),N T +1 k + k j=1 - → θ C n+1 j + - → θ I (1),N T +1 k j dxdy (4.12) for any h ∈ C 1 b (R 2 ). Indeed, since f (ζ i -ζ i-1 )θ i ∈ M i-1,n ( X, Ȳ , -1/2) and f (ζ k -ζ k-1 )I (2) k ( - → θ e,Y k ) ∈ M k-1,n ( X, Ȳ , -1), for some c := c(T, b Y , κ) > 4κ, it holds E p(T -ζ n , Xn , Ȳn , x, y) n+1 k=1 (ζ k -ζ k-1 ) - → θ I (2),n+1 k τ n+1 ≤ C n+1 (R 2 ) n q4κ (T -ζ n , x n , y n , x, y) n+1 k=1 (ζ k -ζ k-1 )(1 -F (T -ζ n )) -1 n i=k+1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1/2 × (f (ζ k -ζ k-1 )) -1 (ζ k -ζ k-1 ) -1 k-1 i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1/2 n i=1 q4κ (ζ i -ζ i-1 , x i-1 , y i-1 , x i , y i ) dx n dy n ≤ C n+1 qc (T, x 0 , y 0 , x, y) n+1 k=1 (1 -F (T -ζ n )) -1 n i=1 (f (ζ i -ζ i-1 )) -1 n i=1,i =k (ζ i -ζ i-1 ) -1/2 (4.13)
where, for the first inequality we used the upper-estimate (2.4) and for the last inequality we used Lemma B.3.

From similar arguments, one gets

E p(T -ζ n , Xn , Ȳn , x, y) n+1 k=1 (ζ k -ζ k-1 ) k j=1 - → θ C n+1 j + - → θ I (1),N T +1 k j τ n+1 ≤ C n+1 qc (T, x 0 , y 0 , x, y) n+1 k=1 (ζ k -ζ k-1 ) k j=1 (1 -F (T -ζ n )) -1 n i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1/2 [1 + 1 {i=k} (ζ i -ζ i-1 ) -1/2 ]. (4.14)
Now, from the upper-bounds (4.13) and (4.14) as well as the identity (2.6), we conclude

n≥0 E E p(T -ζ n , Xn , Ȳn , x, y) n+1 k=1 (ζ k -ζ k-1 ) - → θ I (2),n+1 k + k j=1 - → θ C n+1 j + - → θ I (1),n+1 k j τ n+1 1 {N T =n} ≤ qc (T, x 0 , y 0 , x, y) n≥0 C n+1 E n+1 k=1 (1 -F (T -ζ n )) -1 n i=1 (f (ζ i -ζ i-1 )) -1 n i=1,i =k (ζ i -ζ i-1 ) -1/2 + n+1 k=1 (ζ k -ζ k-1 ) k j=1 (1 -F (T -ζ n )) -1 n i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1/2 [1 + 1 {i=k} (ζ i -ζ i-1 ) -1/2 ]
≤ qc (T, x 0 , y 0 , x, y)

n≥0 C n+1 [(n + 1) + (n + 1)(n + 2)/2]T (n+1)/2 Γ n (1/2) Γ(1 + n/2)
= CT 1/2 qc (T, x 0 , y 0 , x, y). (4.15)

From the preceding inequality and Fubini's theorem, we thus get

E p(T -ζ N T , XN T , ȲN T , x, y) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (2),N T +1 k + k j=1 - → θ C N T +1 j + - → θ I (1),N T +1 k j ≤ CT 1/2 qc (T, x 0 , y 0 , x, y) (4.16)
for some positive constant C := C(T ) such that T → C(T ) is non-decreasing. Applying again Fubini's theorem allows to complete the proof of (4.12). Hence,

T ∂ y0 E[h(X T , Y T )] = R 2 h(x, y) ∂ y0 p(T, x 0 , y 0 , x, y) dxdy = R 2 h(x, y) E p(T -ζ N T , XN T , ȲN T , x, y) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (2),N T +1 k + k j=1 - → θ C N T +1 j + - → θ I (1),N T +1 k j dxdy for any h ∈ C 1 b (R 2 ).
A monotone class argument allows to conclude that the preceding identity is still valid for any bounded and measurable map h defined over R 2 and a standard approximation argument allows to extend it to h ∈ B γ (R 2 ) for any 0 ≤ γ < c -1 , c being the positive constant appearing in (4.16). We eventually conclude from the preceding identity, (4.15) combined with Fubini's theorem that

T ∂ y0 E[h(X T , Y T )] = E h( XN T +1 , ȲN T +1 ) N T +1 k=1 (ζ k -ζ k-1 ) - → θ I (2),N T +1 k + k j=1 - → θ C N T +1 j + - → θ I (1),N T +1 k j for any h ∈ B γ (R 2 ) such that 0 ≤ γ < c -1 .
Step 3: L p (P)-moments for a renewal process with Beta jump times.

From the above formula, the proof of the L p (P)-moment estimate when N is a renewal process with Beta jump times follows by similar arguments as those employed at step 3 of the proof of Theorem 3.1. We omit the remaining technical details.

Numerical Results

In this section, we provide some numerical results for the unbiased Monte Carlo algorithm that stems from the probabilistic representation formula established in Theorem 3.1 and the Bismut-Elworthy-Li formulae of Theorem 4.1 for the couple (S T , Y T ) that allows to compute the Delta and the Vega related to the option price of the vanilla option with payoff h(S T ). We here consider the unique strong solution associated to the SDE (1.1) for three different models corresponding to three different diffusion coefficient function σ S and two different options, namely Call and digital Call options with maturity T and strike K, with payoff functions h(x, y) = (exp(x) -K) + and h(x, y) = 1 {exp(x)≥K} respectively. For these three models, the drift function of the volatility process is defined by b Y (x) = λ Y (µ -x) and we fix the parameters as follows: T = 0.5, r = 0.03, K = 1.5,

x 0 = ln(s 0 ) = 0.4, Y 0 = 0.2, σ Y (.) ≡ σ Y = 0.2, λ Y = 0.5, µ = 0.
3 and ρ = 0.6. We also consider two type of renewal process N : a Poisson process with intensity parameter λ = 0.5 and a renewal process with Beta(1 -α, 1) jump times with parameters α = 0.1 and τ = 2.

5.1. Black-Scholes Model. We first consider the simple (toy) example corresponding to the Black-Scholes dynamics

dS t = rS t dt + σ S S t dW t , dY t = b Y (Y t ) dt + σ Y (Y t )dB t , d B, W t = ρdt, ρ ∈ (-1, 1
). with constant diffusion coefficient function σ S (.) ≡ σ S > 0. The law of (S T , Y T ) can be computed explicitly so that analytical formulas are available for the price, Delta and Vega. Note that the discount factor e -rT has been added in our probabilistic representation formula for comparison purposes. In this example, we importantly remark that the dynamics of the Euler scheme writes (5.1)

   Xi+1 = Xi + r -1 2 a S (ζ i+1 -ζ i ) + σ S ζ i+1 -ζ i Z 1 i+1 , Ȳi+1 = m i + σ Y ζ i+1 -ζ i ρZ 1 i+1 + 1 -ρ 2 Z 2 i+1 ,
with m i = m ζi+1-ζi ( Ȳi ) = µ+( Ȳi -µ)e -λ(ζi+1-ζi) . Also, the weights (θ i ) 1≤i≤N T +1 in the probabilistic representation (3.2) of Theorem 3.1 greatly simplifies, namely

θ i = (f (ζ i -ζ i-1 )) -1 I (2) i (b i Y ), 1 ≤ i ≤ N T , and θ N T +1 = (1 -F (T -ζ N T )) -1
. We perform M 1 = 10 7 Monte Carlo path simulations to approximate the price as well as the two Greeks and compare them with the corresponding values obtained using the standard Monte Carlo method combined with an Euler-Maruyama approximation scheme for the dynamics (1.1) with M 2 = 160000 Monte Carlo simulations paths and mesh size δ = T /n where n = 200. The Delta and Vega are obtained using the Monte Carlo finite difference approach combined with the Euler-Maruyama discretization scheme, that is, denoting by E n M2 (s 0 , y 0 ) the Monte Carlo estimator with Euler-Maruyama scheme, we compute (E n M2 (s 0 +ε, y 0 )-E n M2 (s 0 , y 0 ))/ε and (E n M2 (s 0 , y 0 +ε)-E n M2 (s 0 , y 0 ))/ε respectively with ε = 10 -2 . The numerical results for the three different quantities are summarized in Table 1, Table 2, Table 3 respectively. The first column provides the value of the parameter σ S . The second column stands for the value of the price, Delta or Vega obtained by the corresponding Black-Scholes formula. The third and fourth columns correspond to the value obtained by the Euler-Maruyama discretization scheme together with its 95% confidence interval. The fifth and sixth (resp. seventh and eighth ) columns provide the estimated value with its 95% confidence interval by our method in the case of Exponential sampling (resp. Beta sampling). We observe a good behaviour of the unbiased estimators for all three quantities and for all the values of the parameter σ S .

aS,i = ζ i+1 -ζ i 0 σ1 µ + ( Ȳi -µ)e -λ Y s + σ2 2 ds, = (σ1µ + σ2) 2 (ζi+1 -ζi) + σ 2 1 ( Ȳi -µ) 2 1 -e -2λ Y (ζ i+1 -ζ i ) 2λY + 2σ1(σ1µ + σ2)( Ȳi -µ) 1 -e -λ Y (ζ i+1 -ζ i ) λY , a S,i = σ 2 1 ( Ȳi -µ) 1 -e -2λ Y (ζ i+1 -ζ i ) λY + 2σ1(σ1µ + σ2) 1 -e -λ Y (ζ i+1 -ζ i ) λY , ρi = ρ ζ i+1 -ζ i 0 α µ + ( Ȳi -µ)e -λs + C ds σS,i √ ζi+1 -ζi = ρ α( Ȳi -µ)(1 -e -λ(ζ i+1 -ζ i ) )/λ + (σ1µ + σ2)(ζi+1 -ζi) σS,i √ ζi+1 -ζi , ρ i = ρ σS,i σ1(1 -e -λ Y (ζ i+1 -ζ i ) )/λY -σ S,i σ1( Ȳi -µ)(1 -e -λ Y (ζ i+1 -ζ i ) )/λY + (σ1µ + σ2)(ζi+1 -ζi) aS,i √ ζi+1 -ζi .
The parameters for the unbiased Monte Carlo method and the Monte Carlo method combined with an Euler-Maruyama approximation scheme are chosen as in the first example. The numerical results related to the price, Delta and Vega are provided in Table 4, Table 5, Table 6 respectively for the Call option and in Table 7, Table 8, Table 9 for the digital Call option. In spite of the fact that the main assumptions are not satisfied, we again observe a good performance of the unbiased estimators for all three quantities and for all the values of the parameters σ 1 , σ 2 , except for the computation of the Vega of a Call option for large values of σ 1 and σ 2 .

σ

5.3.

A model with a periodic diffusion coefficient function. In our last example, the volatility of spot price takes the following form σ S (x) = σ 1 cos(x) + σ 2 where σ 1 and σ 2 are two positive constants such that σ 2 -σ 1 > 0 in order to ensure that (ND) is satisfied. Here, the coefficients appearing in the dynamics (2.7) write

aS,i = ζ i+1 -ζ i 0 σ1 cos µ + ( Ȳi -µ)e -λ Y s + σ2 2 ds, a S,i = -2α ζ i+1 -ζ i 0 e -λ Y s sin µ + ( Ȳi -µ)e -λ Y s σ1 cos µ + ( Ȳi -µ)e -λ Y s + σ2 ds, ρi = ρ ζ i+1 -ζ i 0 σ1 cos µ + ( Ȳi -µ)e -λ Y s + σ2 ds σS,i √ ζi+1 -ζi , ρ i = -ρ σ1σS,i ζ i+1 -ζ i 0 e -λ Y s sin µ + ( Ȳi -µ)e -λ Y s ds + σ S,i ζ i+1 -ζ i 0 σ1 cos µ + ( Ȳi -µ)e -λ Y s + σ2 ds aS,i √ ζi+1 -ζi
and no analytical expressions are available. However, a simple numerical integration method can be employed for the computation of the above integrals. We here use Simpson's 3/8 rule which for a real-valued C 4 ([0, T ]) function g writes as follows

∀t ∈ [0, T ], t 0 g(s)ds ≈ t 8 g(0) + 3g t 3 + 3g 2t 3 + g(t)
with an error given by g (4) (t )T 5 /6480 for some t ∈ [0, T ]. The parameters of the unbiased Monte Carlo method and the Monte Carlo Euler-Maruyama scheme remain unchanged. The numerical results related to the price, Delta and Vega are provided in Table 10, Table 11

Lf (x, y) = (r - 1 2 a S (y))∂ x f (x, y) + 1 2 a S (y)∂ 2 x f (x, y) + b Y (y)∂ y f (x, y) + 1 2 a Y (y)∂ 2 y f (x, y) + ρ(σ S σ Y )(y)∂ 2 x,y f (x, y), Lt f (x, y) = (r - 1 2 a S (m t (y 0 )))∂ x f (x, y) + 1 2 a S (m t (y 0 ))∂ 2 x f (x, y) + b Y (m t (y 0 ))∂ y f (x, y) + 1 2 a Y (m t (y 0 ))∂ 2 y f (x, y) + ρ(σ S σ Y )(m t (y 0 ))∂ 2
x,y f (x, y)

for any f ∈ C 2 b (R 2 ).
Step 1: Probabilistic representation for a bounded and continuous map h

We establish a first order expansion of the Markov semigroup (P t ) t≥0 around ( Pt ) t≥0 . We apply Itô's rule to the map [0, t] × R 2 (s, x, y)

→ P t-s h(x, y) ∈ C 1,2 b ([0, t] × R 2 ) for h ∈ C ∞ b (R 2 ), observing that ∂ s P t-s h(x, y) = -LP t-s h(x, y). We obtain h( Xt , Ȳt ) = P t h(x 0 , y 0 ) + t 0 ∂ s P t-s h( Xs , Ȳs ) + Ls P t-s h( Xs , Ȳs ) ds + M t = P t h(x 0 , y 0 ) + t 0 ( Ls -L)P t-s h( Xs , Ȳs ) ds + M t
where M := (M t ) t≥0 is a square integrable martingale. We then take expectation in the previous expression, make use of Fubini's theorem and finally let t ↑ T by dominated convergence theorem so that

P T h(x 0 , y 0 ) = PT h(x 0 , y 0 ) + T 0 E[(L -Ls )P T -s h( Xs , Ȳs )] ds = E[h( Xx0 T , Ȳ y0 T )] + T 0 E 1 2 (a S ( Ȳ y0 s ) -a S (m s (y 0 ))[∂ 2 x P T -s h( Xx0 s , Ȳ y0 s ) -∂ x P T -s h( Xx0 s , Ȳ y0 s ))] ds + T 0 E 1 2 (a Y ( Ȳ y0 s ) -a Y (m s (y 0 )))∂ 2 y P T -s h( Xx0 s , Ȳ y0 s ) + (b Y ( Ȳ y0 s ) -b Y (m s (y 0 )))∂ y P T -s h( Xx0 s , Ȳ y0 s ) ds + T 0 E ρ((σ S σ Y )( Ȳ x0 s ) -(σ S σ Y )(m s (y 0 )))∂ 2 x,y P T -s h( Xx0 s , Ȳ y0 s ) ds. (A.1)
We now rewrite the previous first order expansion using the Markov chain ( Xi , Ȳi ) 0≤i≤N T +1 and the renewal process N . From the previous identity, the definition of θ N T +1 in (3.4) and the identity (2.6), we directly obtain

P T h(x 0 , y 0 ) = E[h( XN T +1 , ȲN T +1 )θ N T +1 1 {N T =0} ] (A.2) + E ((1 -F (T -ζ 1 ))f (ζ 1 )) -1 1 {N T =1} 1 2 (a S ( Ȳ1 ) -a S (m 0 ))D (1,1) 1 P T -ζ1 h( X1 , Ȳ1 ) - 1 2 (a S ( Ȳ1 ) -a S (m 0 ))D (1)
1 P T -ζ1 h( X1 , Ȳ1 ) + 1 2 (a Y ( Ȳ1 ) -a Y (m 0 ))D (2,2) 1 P T -ζ1 h( X1 , Ȳ1 ) +(b Y ( Ȳ1 ) -b Y (m 0 ))D (2)
1 P T -ζ1 h( X1 , Ȳ1 ) + ρ((σ S σ Y )( Ȳ1 ) -(σ S σ Y )(m 0 ))D (1,2) 1 P T -ζ1 h( X1 , Ȳ1 ) = E[h( XN T +1 , ȲN T +1 )θ N T +1 1 {N T =0} ] + E ((1 -F (T -ζ 1 ))f (ζ 1 )) -1 1 {N T =1} c 1 S D (1,1) 1 P
T -ζ1 h( X1 , Ȳ1 ) -c 1 S D (1)
1 P T -ζ1 h( X1 , Ȳ1 ) + c 1 Y D (2,2) 1 P T -ζ1 h( X1 , Ȳ1 ) + b 1 Y D (2) 1 P T -ζ1 h( X1 , Ȳ1 ) + c 1 S,Y D (1,2) 1
P T -ζ1 h( X1 , Ȳ1 ) .

(A.3)
Next, we apply the IBP formula (2.14) with respect to the random vector ( X1 , Ȳ1 ) in the above expression. In order to do that rigorously, one first has to take the conditional expectation E 0,1 [.] in the second term of the above equality. We thus obtain

E 0,1 c 1 S D (1,1) 1 P
T -ζ1 h( X1 , Ȳ1 ) -c 1 S D (1)
1 P T -ζ1 h( X1 , Ȳ1 ) + c 1 Y D (2,2) 1 P T -ζ1 h( X1 , Ȳ1 ) + b 1 Y D (2) 1 P T -ζ1 h( X1 , Ȳ1 ) + c 1 S,Y D (1,2) 1 P
T -ζ1 h( X1 , Ȳ1 ) = E 0,1 I (1,1) 1 (c 1 S ) -I (1) 1 (c 1 S ) + I (2,2) 1 (c 1
Y ) + I (2) 1 (b 1
Y ) + I (1,2) 1 (c 1
S,Y ) P T -ζ1 h( X1 , Ȳ1 ) .
From Lemma B.2 and the estimate (2.22), we get

(A.4) E 0,1 I (1,1) 1 (c 1 S ) -I (1) 1 (c 1 S ) + I (2,2) 1 (c 1
Y ) + I (2) 1 (b 1
Y ) + I (1,2) 1 (c 1
S,Y ) P T -ζ1 h( X1 , Ȳ1 ) ≤ C T |h| ∞ ζ -1/2 1
for some positive constant C T such that T → C T is non-decreasing. The previous estimate yields an integrable time singularity. Indeed, from the previous estimate and (2.6), one directly gets

E ((1 -F (T -ζ 1 ))f (ζ 1 )) -1 1 {N T =1} E 0,1 I (1,1) 1
(c 1 S ) -I

(1)

1 (c 1 S ) + I (2,2) 1 (c 1 Y ) + I (2) 1 (b 1 Y ) + I (1,2) 1 (c 1 S,Y ) P T -ζ1 h( X1 , Ȳ1 ) ≤ CE[((1 -F (T -ζ 1 ))f (ζ 1 )) -1 ζ -1/2 1 1 {N T =1} ] ≤ C T 0 s -1/2 1 ds 1 < ∞.
Coming back to (A.3), we thus derive

P T h(x 0 , y 0 ) = E[h( XN T +1 , ȲN T +1 )θ N T +1 1 {N T =0} ] + E ((1 -F (T -ζ 1 ))f (ζ 1 )) -1 1 {N T =1} × I (1,1) 1 (c 1 S ) -I (1)
1 (c 1 S ) + I (2,2) 1 (c 1 Y ) + I (2) 1 (b 1 Y ) + I (1,2) 1 (c 1 S,Y ) P T -ζ1 h( X1 , Ȳ1 ) = E[h( XN T +1 , ȲN T +1 )θ N T +1 1 {N T =0} ] + E P T -ζ1 h( X1 , Ȳ1 )θ 2 θ 1 1 {N T =1} . (A.5)
Our aim now is to iterate the above first order expansion. We prove by induction the following formula: for any positive integer n, one has

P T h(x 0 , y 0 ) = n-1 j=0 E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i 1 {N T =j} + E P T -ζn h( Xn , Ȳn ) n+1 i=1 θ i 1 {N T =n} . (A.6)
The case n = 1 corresponds to (A.5). We thus assume that (A.6) holds at step n. We expand the last term appearing in the right-hand side of the previous equality using again (A.1), by then applying Lemma B.1 and by finally performing IBPs as before.

To be more specific, using the notations introduced in Subsection 2.2, from (A.1) and a change of variable, for any ζ ∈ [0, T ], one has

P T -ζ h(x, y) = E[h( Xζ,x T , Ȳ ζ,x T )] + T ζ E 1 2 (a S ( Ȳ ζ,y s ) -a S (m s-ζ (y))[∂ 2 x P T -s h( Xζ,x s , Ȳ ζ,y s ) -∂ x P T -s h( Xζ,x s , Ȳ ζ,y s ))] ds + T ζ E 1 2 (a Y ( Ȳ ζ,y s ) -a Y (m s-ζ (y)))∂ 2 y P T -s h( Xζ,x s , Ȳ ζ,y s ) + (b Y ( Ȳ ζ,y s ) -b Y (m s-ζ (y)))∂ y P T -s h( Xζ,x s , Ȳ ζ,y s ) ds + T ζ E ρ((σ S σ Y )( Ȳ ζ,x s ) -(σ S σ Y )(m s-ζ (y)))∂ 2 x,y P T -s h( Xζ,x s , Ȳ ζ,y s ) ds.
We take ζ = ζ n , (x, y) = ( XN T , ȲN T ) in the previous equality, then multiply it by n+1 i=1 θ i 1 {N T =n} and finally take expectation. We obtain

E P T -ζn h( Xn , Ȳn ) n+1 i=1 θ i 1 {N T =n} = E h( Xζn, Xn T , Ȳ ζn, Ȳn T ) n+1 i=1 θ i 1 {N T =n} + E n+1 i=1 θ i 1 {N T =n} T ζn 1 2 (a S ( Ȳ ζn, Ȳn s ) -a S (m s-ζn ( Ȳn ))[∂ 2 x P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) -∂ x P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s
))]ds

+ E n+1 i=1 θ i 1 {N T =n} T ζn 1 2 (a Y ( Ȳ ζn, Ȳn s ) -a Y (m s-ζn ( Ȳn )))∂ 2 y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s )ds +
E n+1 i=1 θ i 1 {N T =n} T ζn (b Y ( Ȳ ζn, Ȳn s ) -b Y (m s-ζn ( Ȳn )))∂ y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s )ds +
E n+1 i=1 θ i 1 {N T =n} T ζn ρ((σ S σ Y )( Ȳ ζn, Ȳn s ) -(σ S σ Y )(m s-ζn ( Ȳn )))∂ 2 x,y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) ds . (A.7)
Now, from the very definition of the Markov chain ( Xi , Ȳi ) 0≤i≤N T +1 and of the weight sequence (θ i ) 1≤i≤N T +1 of Theorem 3.1, the first term of the above equality can be written as

(A.8) E h( Xζn, Xn T , Ȳ ζn, Ȳn T ) n+1 i=1 θ i 1 {N T =n} = E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i 1 {N T =n} .
We now look at the second, third, fourth and fifth terms. Let us deal with the third and fourth terms. The others are treated in a similar manner and we will omit some technical details. We first take its conditional expectation w.r.t

{ζ 1 = t 1 , • • • , ζ n = t n , N T = n} and introduce the measurable function G(t 1 , • • • , t n , s, T ) := E n+1 i=1 θ i T ζn 1 2 (a Y ( Ȳ ζn, Ȳn s ) -a Y (m s-ζn ( Ȳn )))∂ 2 y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) + (b Y ( Ȳ ζn, Ȳn s ) -b Y (m s-ζn ( Ȳn )))∂ y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) ds ζ 1 = t 1 , • • • , ζ n = t n , N T = n which satisfies |G(t 1 , • • • , t n , s, T )| ≤ CE n+1 i=1 |θ i | 1 + |W s -W ζn | + | W s -W ζn | ζ 1 = t 1 , • • • , ζ n = t n , N T = n ≤ CE n+1 i=1 |θ i ||ζ 1 = t 1 , • • • , ζ n = t n , N T = n
where we used the boundedness of a Y , the Lipschitz regularity of b Y , the inequalities sup 0≤t≤T |∂ y P t h| ∞ ≤ C for = 1, 2 and, for the last inequality the fact that, conditionally on {ζ

1 = t 1 , • • • , ζ n = t n , N T = n}, the random variables (W s -W ζn , W s -W ζn ) are independent of the sigma-field G n . Recall now that P(N T = n, ζ 1 ∈ dt 1 , • • • , ζ n ∈ dt n ) = (1 -F (T -t n )) n-1 j=0 f (t j+1 -t j )dt 1 , • • •
, dt n on the set ∆ n (T ) so that from Lemma B.2 and the estimate (2.22), we obtain

E T ζn |G(ζ 1 , • • • , ζ n , s, T )|1 {N T =n} ds ≤ C T tn ∆n(T ) n i=1 (t i -t i-1 ) -1/2 dt 1 • • • dt n dt n+1 < ∞.
Hence, by Lemma B.1, it holds

E 1 {N T =n} n+1 i=1 θ i T ζn 1 2 (a Y ( Ȳ ζn, Ȳn s ) -a Y (m s-ζn ( Ȳn )))∂ 2 y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) + (b Y ( Ȳ ζn, Ȳn s ) -b Y (m s-ζn ( Ȳn )))∂ y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) ds = E 1 {N T =n} T ζn G(ζ 1 , • • • , ζ n , s, T ) ds = E n i=1 θ i (1 -F (T -ζ n+1 )) -1 (f (ζ n+1 -ζ n )) -1 1 2 (a Y ( Ȳn+1 ) -a Y (m n ))D (2,2) n+1 P T -ζn+1 h( Xn+1 , Ȳn+1 ) + (b Y ( Ȳn+1 ) -b Y (m n ))D (2)
n+1 P T -s h( Xn+1 , Ȳn+1 ) 1 {N T =n+1} .
Finally, we take the conditional expectation E n,n+1 [.] inside the above expectation and then employ the IBP formula (2.14), two times w.r.t. the diffusion coefficient and one time w.r.t the drift coefficient as done before. We obtain

E n i=1 θ i (1 -F (T -ζ n+1 )) -1 (f (ζ n+1 -ζ n )) -1 1 2 (a Y ( Ȳn+1 ) -a Y (m n ))D (2,2)
n+1 P T -ζn+1 h( Xn+1 , Ȳn+1 ) + (b Y ( Ȳn+1 ) -b Y (m n ))D (2)
n+1 P T -ζn+1 h( Xn+1 , Ȳn+1 ) 1 {N T =n+1} = E n i=1 θ i (1 -F (T -ζ n+1 )) -1 (f (ζ n+1 -ζ n )) -1 [I (2,2) n+1 (c n+1 Y
) + I

n+1 (b n+1 Y )]P T -ζn+1 h( Xn+1 , Ȳn+1 )1 {N T =n+1} . (2)
In a completely analogous manner, we derive

E n+1 i=1 θ i 1 {N T =n} T ζn 1 2 (a S ( Ȳ ζn, Ȳn s ) -a S (m s-ζn ( Ȳn ))[∂ 2 x P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) -∂ x P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ))]ds = E n i=1 θ i (1 -F (T -ζ n+1 )) -1 (f (ζ n+1 -ζ n )) -1 [I (1,1) n+1 (c n+1 S
) -I

n+1 (c n+1 S )]P T -ζn+1 h( Xn+1 , Ȳn+1 )1 {N T =n+1} and E n+1 i=1 θ i 1 {N T =n} T ζn ρ((σ S σ Y )( Ȳ ζn, Ȳn s ) -(σ S σ Y )(m s-ζn ( Ȳn )))∂ 2 x,y P T -s h( Xζn, Xn s , Ȳ ζn, Ȳn s ) ds = E n i=1 θ i (1 -F (T -ζ n+1 )) -1 (f (ζ n+1 -ζ n )) -1 I (1,2) n+1 (c n+1 Y,S )P T -ζn+1 h( Xn+1 , Ȳn+1 )1 {N T =n+1} . (1)
Summing the three previous identities, we obtain that the sum of the second, third, fourth and fifth term in the right-hand side of (A.7) is equal to

E n i=1 θ i (1 -F (T -ζ n+1 )) -1 (f (ζ n+1 -ζ n )) -1 × I (1,1) i+1 (a n+1 S ) -I (1) n+1 (a n+1 S ) + I (2,2) n+1 (a n+1 Y ) + I (2) n+1 (b n+1 Y ) + I (1,2) n+1 (a n+1 Y,S ) P T -ζn+1 h( Xn+1 , Ȳn+1 )1 {N T =n+1} = E n+2 i=1 θ i P T -ζn+1 h( Xn+1 , Ȳn+1 )1 {N T =n+1}
where we used the very definitions (3.3) and (3.4) of the weights (θ i ) 1≤i≤N T +1 on the set {N T = n + 1}. This concludes the proof of (A.6) at step n + 1.

To conclude it remains to prove the absolute convergence of the first sum and the convergence to zero of the last term in (A.6). These two results follow directly from the boundedness of h and the general estimates on the product of weights established in Lemma B.2.

Indeed, from Lemma B.2, the estimate (2.22), the tower property of conditional expectation and the identity (2.6), we obtain

E h( XN T +1 , ȲN T +1 ) N T +1 i=1 |θ i | 1 {N T =j} ≤ C j |h| ∞ E (1 -F (T -ζ j )) -1 j i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1 2 1 {N T =j} = C j |h| ∞ ∆j (T ) j i=1 (s i -s i-1 ) -1 2 ds j = C j |h| ∞ T j 2 Γ j (1/2) Γ(1 + j/2) which in turn yields n-1 j=0 E h( XN T +1 , ȲN T +1 ) N T +1 i=1 |θ i | 1 {N T =j} ≤ |h| ∞ n≥0 (CT 1/2 ) n Γ(1 + n/2) = |h| ∞ E 1/2,1 (CT 1/2 )
so that the series converge absolutely. Similarly,

E P T -ζn h( Xn , Ȳn ) n+1 i=1 θ i 1 {N T =n} ≤ C n |h| ∞ T n 2 Γ n (1/2) Γ(1 + n/2)
so that the remainder indeed vanishes as n goes to infinity. We thus conclude (A.9)

P T h(x 0 , y 0 ) = n≥0 E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i 1 {N T =n} = E h( XN T +1 , ȲN T +1 ) N T +1 i=1 θ i for any h ∈ C 2 b (R 2
). We eventually extend the above representation formula to any bounded and continuous function h using a standard approximation argument. The remaining technical details are omitted.

Step 2: Extension to measurable map h satisfying the growth assumption (3.1) We first extend the previous result to any bounded and measurable h. This follows from a monotone class argument that we now detail.

Let us first consider h ∈ C b (R 2 ). From Fubini's theorem, it holds

E h( Xn+1 , Ȳn+1 ) n+1 i=1 θ i N T = n, ζ n = R 2 h(x, y) E p(T -ζ n , Xn , Ȳn , x, y) n+1 i=1 θ i ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n ) N T = n,
n+1 i=1 θ i ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n ) N T = n, ζ n ≤ (R 2 ) n p(T -ζ n , x n , y n , x, y) n i=1 |θ i (x i-1 , y i-1 , x i , y i , ζ n )|p(ζ i -ζ i-1 , x i-1 , y i-1 , x i , y i ) dx n dy n ≤ (C T ) n (R 2 ) n q4κ (T -ζ n , x n , y n , x, y) n i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1 2 q4κ (ζ i -ζ i-1 , x i-1 , y i-1 , x i , y i ) dx n dy n ≤ (C T ) n n i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1 2 qc (T, x 0 , y 0 , x, y)
for some c := c(T, b Y , κ) ≥ 4κ. Hence, from (A.9) and again Fubini's theorem, justified by the previous estimate and the fact that

E (C T ) N T N T i=1 (f (ζ i -ζ i-1 )) -1 (ζ i -ζ i-1 ) -1 2
< ∞, one has

P T h(x 0 , y 0 ) = R 2 h(x, y)E p(T -ζ N T , XN T , ȲN T , x, y) N T +1 i=1 θ i dxdy (A.10)
for any h ∈ C b (R 2 ). Moreover, from the previous computations, the following upper-bound holds

E p(T -ζ N T , XN T , ȲN T , x, y) N T +1 i=1 θ i = n≥0 ∆n(T ) E p(T -s n , Xn , Ȳn , x, y) n+1 i=1 θ i N T = n, ζ n = (s 1 , • • • , s n ) (1 -F (T -s n )) n i=1 f (s i -s i-1 ) ds n ≤ n≥0 ∆n(T ) (C T ) n n i=1 (s i -s i-1 ) -1
2 ds n qc (T, x 0 , y 0 , x, y)

= E 1/2,1 (CT 1/2 )q c (T, x 0 , y 0 , x, y).

(A.11)

It now follows from (A.10) and a monotone class argument that the probabilistic representation formula (3.2) is valid for any real-valued bounded and measurable map h defined over R 2 . The extension to any measurable map h satisfying the growth assumption: |h(x, y)| ≤ C exp(γ(|x| 2 + |y| 2 )) for any 0 ≤ γ < (2κ) -1 , follows from the integral representation (A.10), the upper-bound (A.11) combined with a standard approximation argument. Remaining technical details are omitted.

Step 3: Finite L p (P)-moment for the probabilistic representation If N is a renewal process with Beta(α, 1) jump times then f (s i -

s i-1 ) = 1-α τ 1-α 1 (si-si-1) α 1 [0,τ ] (s i -s i-1 ) and 1 -F (T -s n ) = 1 -T -sn τ 1-α ≥ 1 -( T τ ) 1-α ,
n+1 i=1 |θ i ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n )| p N T = n, ζ n ≤ (R 2 ) n p(T -ζ n , x n , y n , x, y) n i=1 |θ i (x i-1 , y i-1 , x i , y i , ζ n )| p p(ζ i -ζ i-1 , x i-1 , y i-1 , x i , y i ) dx n dy n ≤ C n (R 2 ) n q4κ (T -ζ n , x n , y n , x, y) n i=1 (f (ζ i -ζ i-1 )) -p (ζ i -ζ i-1 ) -p 2 q4κ (ζ i -ζ i-1 , x i-1 , y i-1 , x i , y i ) dx n dy n ≤ C n n i=1 (ζ i -ζ i-1 ) αp (ζ i -ζ i-1 ) -p 2 qc (T, x 0 , y 0 , x, y)
for some c := c(T, b Y , κ) ≥ 4κ. Now, using the fact that

E C n N T i=1 (ζ i -ζ i-1 ) αp (ζ i -ζ i-1 ) -p 2 < ∞ as soon as p( 1 2 -α) < 1 -α and that h ∈ B γ (R 2 )
, from the previous computation, we obtain

E[|h( XN T +1 , ȲN T +1 )| p N T +1 i=1 |θ i | p ] = E E |h( Xn+1 , Ȳn+1 )| p n+1 i=1 |θ i | p N T , ζ n ≤ E C n N T i=1 (ζ i -ζ i-1 ) αp (ζ i -ζ i-1 ) -p 2 R 2
e γp(|x| 2 +|y| 2 ) qc (T, x 0 , y 0 , x, y) dxdy.

To conclude the proof, it suffices to note that the above space integral is finite as soon as 0 ≤ γp < c -1 .

A.2. Proof of Lemma 4.1. Since h ∈ C 1 p (R 2 ) and E i,n | Xi+1 | q +| Ȳi+1 | q < ∞ a.s., for any q ≥ 1, under (AR), one may differentiate under the (conditional) expectation and deduce that (x, y)

→ E i,n h( Xi+1 , Ȳi+1 )θ i+1 |( Xi , Ȳi ) = (x, y) ∈ C 1 p (R 2 ) for any i ∈ {0, • • • , n}.
The rest of the proof is divided into three parts.

Step 1: proofs of (4.1) and (4.

2)

The transfer of derivatives formulae (4.1) and (4.2) are easily obtained by differentiating under expectation (which is allowed by the polynomial growth at infinity of h) noting from the definition of the Markov chain X that

∂ s0 X1 = ∂ s0 ln(s 0 ) = 1 s0 and ∂ Xi h( Xi+1 , Ȳi+1 ) = ∂ Xi+1 h( Xi+1 , Ȳi+1 )∂ Xi Xi+1 = ∂ Xi+1 h( Xi+1 , Ȳi+1
). Observe as well that from (2.15), the fact that

∂ Xi c i+1 S = ∂ Xi c i+1 Y = ∂ Xi b i+1 Y = ∂ Xi c i+1 Y,S = ∂ Xi I (1) i+1 (1) = ∂ Xi I (2)
i+1 (1) = 0 and the very definition of the random variables (θ i ) 1≤i≤n+1 , one has ∂ Xi θ i+1 = 0. This gives the identities (4.1) and (4.2).

Step 2: proofs of (4.3) and (4.4)

The proofs of (4.3) and (4.4) are more involved. Let us prove (4.3). We proceed by considering the difference between the term appearing on the left-hand side and the first two terms appearing on the right-hand side of (4.3). On the one hand, using the IBP formula (2.14) and (2.8), we get

∂ Ȳi E i,n h( Xi+1 , Ȳi+1 )θ i+1 = E i,n ∂ Xi+1 h( Xi+1 , Ȳi+1 )∂ Ȳi Xi+1 θ i+1 + E i,n ∂ Ȳi+1 h( Xi+1 , Ȳi+1 )∂ Ȳi Ȳi+1 θ i+1 + E i,n h( Xi+1 , Ȳi+1 )∂ Ȳi θ i+1 = E i,n h( Xi+1 , Ȳi+1 ) I (1) i+1 (∂ Ȳi Xi+1 θ i+1 ) + I (2) i+1 (m i θ i+1 ) + E i,n h( Xi+1 , Ȳi+1 )I (2) i+1 σ Y,i ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 θ i+1 + E i,n h( Xi+1 , Ȳi+1 )I (2) i+1 σ Y,i ρ i 1 -ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 θ i+1 + E i,n h( Xi+1 , Ȳi+1 )∂ Ȳi θ i+1 .
On the other hand, again from the IBP formula (2.14), we obtain

E i,n ∂ Xi+1 h( Xi+1 , Ȳi+1 ) - → θ e,X i+1 + E i,n ∂ Ȳi+1 h( Xi+1 , Ȳi+1 ) - → θ e,Y i+1 + E i,n h( Xi+1 , Ȳi+1 ) - → θ c i+1 = E i,n h( Xi+1 , Ȳi+1 ) I (1) i+1 ( - → θ e,X i+1 ) + I (2) i+1 ( - → θ e,Y i+1 ) + - → θ c i+1 .
Combining the two previous identities, we see that the difference

∂ Ȳi E i,n h( Xi+1 , Ȳi+1 )θ i+1 -E i,n ∂ Xi+1 h( Xi+1 , Ȳi+1 ) - → θ e,X i+1 +E i,n ∂ Ȳi+1 h( Xi+1 , Ȳi+1 ) - → θ e,Y i+1 +E i,n h( Xi+1 , Ȳi+1 ) - → θ c
i+1 can be written as

E h( Xi+1 , Ȳi+1 )I (2) i+1 m i θ i+1 - - → θ e,Y i+1 + E h( Xi+1 , Ȳi+1 )∂ Ȳi θ i+1 -E i,n h( Xi+1 , Ȳi+1 )I (1) i+1 ( - → θ e,X i+1
)

+ E i,n h( Xi+1 , Ȳi+1 )I (1) i+1 (∂ Ȳi Xi+1 θ i+1 ) + E h( Xi+1 , Ȳi+1 )I (2) i+1 σ Y,i ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 θ i+1 + E h( Xi+1 , Ȳi+1 )I (2) i+1 σ Y,i ρ i 1 -ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 θ i+1 (A.12) -E h( Xi+1 , Ȳi+1 ) - → θ c i+1 .
Before proceeding, we provide the explicit expression for the quantity ∂ Ȳi θ i+1 . Using the chain rule formula of Lemma 2.2, after some standard but cumbersome computations, we obtain

∂ Ȳi θ i+1 = (f (ζ i+1 -ζ i )) -1 I (1,1) i+1 (∂ Ȳi c i+1 S ) -I (1) i+1 (∂ Ȳi c i+1 S ) + I (2,2) i+1 (∂ Ȳi c i+1 Y ) + I (2) i+1 (∂ Ȳi b i+1 Y ) + I (1,2) i+1 (∂ Ȳi c i+1 Y,S ) - σ S,i σ S,i 2I (1,1) i+1 (c i+1 S ) -I (1) i+1 (c i+1 S ) + I (1,2) i+1 (c i+1 Y,S ) + σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i 2I (2,2) i+1 (c i+1 Y ) + I (2) i+1 (b i+1 Y ) + I (1,2) i+1 (c i+1 Y,S ) + ρ i 1 -ρ 2 i σ Y,i σ S,i I (1,2) i+1 (c i+1 S ) + I (2,1) i+1 (c i+1 S ) -I (2)
i+1 (c i+1 S ) .

Also, after some simple algebraic simplifications using the definitions of -→ θ e,Y i+1 and -→ θ e,X i+1 in (4.3), one obtains

I (2) i+1 m i θ i+1 - - → θ e,Y i+1 = -(f (ζ i+1 -ζ i )) -1 I (2,2) i+1 ∂ Ȳi c i+1 Y + I (1,2) i+1 ∂ Ȳi c i+1 Y,S and
I (1) i+1 ( - → θ e,X i+1 ) = (f (ζ i+1 -ζ i )) -1 I (1,1) i+1 ∂ Ȳi c i+1 S .
Combining the three previous identities and gathering similar terms, we obtain

I ( 2
)
i+1 m i θ i+1 - - → θ e,Y i+1 + ∂ Ȳi θ i+1 -I (1) i+1 ( - → θ e,X i+1 ) = (f (ζ i+1 -ζ i )) -1 -I (1) i+1 (∂ Ȳi c i+1 S ) + I (2) i+1 (∂ Ȳi b i+1 Y ) - σ S,i σ S,i 2I
(1,1)

i+1 (c i+1 S ) -I (1) i+1 (c i+1 S ) + I (1,2) i+1 (c i+1 Y,S ) + σ Y,i σ Y,i - ρ i ρ i 1 -ρ 2 i 2I (2,2) i+1 (c i+1 Y ) + I (2) i+1 (b i+1 Y ) + I (1,2) i+1 (c i+1 Y,S ) (A.13) + ρ i 1 -ρ 2 i σ Y,i σ S,i I (1,2) i+1 (c i+1 S ) + I (2,1) i+1 (c i+1 S ) -I (2)
i+1 (c i+1 S ) .

The previous identity will be used in the next step of the proof. Coming back to (A.12) and using the definition of the weight -→ θ c i+1 allows to conclude the proof of the identity (4.3).

Step

i+1 d i+1 S , D (1,1)
i+1 d i+1 S , D (2)
i+1 d i+1 Y , D (2,2)
i+1 d i+1 Y , D (1)
i+1 d i+1 Y , D (2)
i+1 d i+1 Y , D (1,2)
i+1 d i+1 Y , e Y,i+1 S , D (1)
i+1 e Y,i+1 S , e Y,i+1 Y , D (2) i+1 e Y,i+1 Y , D (1)
i+1 e Y,i+1 S ∈ M i,n ( X, Ȳ , 0) as well as Lemma 2.3, we conclude

f (ζ i+1 -ζ i ) - → θ e,Y i+1 ∈ M i,n ( X, Ȳ , -1/2), i ∈ {0, • • • , n -1} . Note also that e X,i+1 S = ∂ Ȳi c i+1 S = 1 2 (a S ( Ȳi+1 )∂ Ȳi Ȳi+1 -a S (m i )m i ) = 1 2 (a S ( Ȳi+1 ) -a S (m i ))∂ Ȳi Ȳi+1 + 1 2 a S (m i )(∂ Ȳi Ȳi+1 -m i )
so that, using on the one hand the Lipschitz regularity of a S and on the other hand (2.8), from similar arguments as those used in the proof of Lemma B.2, we conclude that

1 2 a S ( Ȳi+1 ) -a S (m i ))∂ Ȳi Ȳi+1 , 1 2 a S (m i )(∂ Ȳi Ȳi+1 -m i ) ∈ M i,n ( X, Ȳ , 1/2)
which in turn implies that e X,i+1 S ∈ M i,n ( X, Ȳ , 1/2). Moreover, standard computations that we omit show that D

(1) i+1 e X,i+1 S ∈ M i,n ( X, Ȳ , 0) so that by Lemma 2.3 we deduce

f (ζ i+1 -ζ i ) - → θ e,X i+1 ∈ M i,n ( X, Ȳ , 0). We now prove that f (ζ i+1 -ζ i ) - → θ c i+1 ∈ M i,n ( X, Ȳ , -1/2) for any i ∈ {0, • • • , n -1}. We use the decomposition f (ζ i+1 -ζ i ) - → θ c i+1 = f (ζ i+1 -ζ i ) I (2) i+1 m i θ i+1 - - → θ e,Y i+1 + ∂ Ȳi θ i+1 -I (1) i+1 ( - → θ e,X i+1 ) + I (1) i+1 ∂ Ȳi Xi+1 f (ζ i+1 -ζ i )θ i+1 + I (2) i+1 σ Y,i ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 + σ Y,i ρ i 1 -ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 f (ζ i+1 -ζ i )θ i+1 . We first prove that f (ζ i+1 -ζ i ) I (2) i+1 m i θ i+1 - - → θ e,Y i+1 + ∂ Ȳi θ i+1 -I (1)
i+1 ( -→ θ e,X i+1 ) ∈ M i,n ( X, Ȳ , -1/2). We investigate each term appearing on the right-hand side of (A.13).

In particular, we first use the fact that

c i+1 S , c i+1 Y , b i+1 Y , c i+1 Y,S ∈ M i,n ( X, Ȳ , 1/2), ∂ Ȳi c i+1 S , ∂ Ȳi b i+1 Y ∈ M i,n ( X, Ȳ , 0)
and the fact that when one applies the differential operators D to these elements the resulting random variables belong to M i,n ( X, Ȳ , 0) for any (α 1 , α 2 ) ∈ {1, 2}

2 . From Lemma 2.3, we thus conclude that the elements I

(1)

i+1 (∂ Ȳi c i+1 S ), I (1,1) i+1 (c i+1 S ), I (1,2) i+1 (c i+1 S ), I (2,1)
i+1 (c i+1 S ), I

i+1 (∂ Ȳi b i+1 Y ), I (2,2) i+1 (c i+1 Y ), I (1,2) i+1 (c i+1 Y,S ) belong to M i,n ( X, Ȳ , -1/2) and that I (1) i+1 (c i+1 S ), I (2)
i+1 (b i+1 Y ) belong to M i,n ( X, Ȳ , 0). Moreover, using (ND), one gets that there exists C > 0 such that for any

i ∈ {0, • • • , n -1}, |σ S,i /σ S,i |+|σ Y,i /σ Y,i |+|σ Y,i /σ S,i |+|ρ i /(1-ρ 2 i )|+|ρ i ρ i /(1- ρ 2 i )| ≤ C. We thus conclude that f (ζ i+1 -ζ i ) I (2) i+1 m i θ i+1 - - → θ e,Y i+1 + ∂ Ȳi θ i+1 -I (1) i+1 ( - → θ e,X i+1 ) ∈ M i,n ( X, Ȳ , -1/2). It thus suffices to prove I (1) i+1 ∂ Ȳi Xi+1 f (ζ i+1 -ζ i )θ i+1 , I (2) i+1 σ Y,i (ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 )f (ζ i+1 -ζ i )θ i+1 and I (2) i+1 σ Y,i ρ i √ 1-ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 f (ζ i+1 -ζ i )θ i+1 belong to M i,n ( X, Ȳ , -1/2).
In order to do this, we remark that

∂ Ȳi Xi+1 = - 1 2 a S,i + σ S,i Z 1 i+1 ∈ M i,n ( X, Ȳ , 1/2), D (1)
i+1 (∂ Ȳi Xi+1 ) = σ S,i σ S,i ∈ M i,n ( X, Ȳ , 0), σ Y,i (ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 ) ∈ M i,n ( X, Ȳ , 1/2), D (2)
i+1 σ Y,i (ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 = σ Y,i σ Y,i ∈ M i,n ( X, Ȳ , 0), σ Y,i ρ i 1 -ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 ∈ M i,n ( X, Ȳ , 1/2), D (2)
i+1 σ Y,i ρ i 1 -ρ 2 i 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 = - ρ i ρ i 1 -ρ 2 i ∈ M i,n ( X, Ȳ , 0) and from Lemma B.2, f (ζ i+1 -ζ i )θ i+1 ∈ M i,n ( X, Ȳ , -1/2). From Lemma 2.3, it follows that f (ζ i+1 -ζ i )θ i+1 ∂ Ȳi Xi+1 ∈ M i,n ( X, Ȳ , 0) and f (ζ i+1 -ζ i )θ i+1 σ Y (m i )m i (ρ(W ζi+1 -W ζi ) + 1 -ρ 2 ( Wζi+1 -Wζi )) ∈ M i,n ( X, Ȳ ,
0). Now following similar computations as those employed in the proof of Lemma B.2 and omitting some technical details we obtain D

(α) i+1 (f (ζ i+1 -ζ i )θ i+1 ) ∈ M i,n (
X, Ȳ , -1) so that from the chain rule formula and Lemma 2.3, the random variables

D (α) i+1 (f (ζ i+1 -ζ i )θ i+1 ∂ Ȳi Xi+1 ), D (α) i+1 (σ Y,i (ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 )f (ζ i+1 -ζ i )θ i+1 ) and D (α) i+1 (σ Y,i ρ i √ 1-ρ 2 i ( 1 -ρ 2 i Z 1 i+1 - ρ i Z 2 i+1 )f (ζ i+1 -ζ i )θ i+1 ) belong to M i,n ( X, Ȳ , -1/2). From Lemma 2.3, we thus conclude that I (1) i+1 ∂ Ȳi Xi+1 f (ζ i+1 - ζ i )θ i+1 , I (2) i+1 σ Y,i (ρ i Z 1 i+1 + 1 -ρ 2 i Z 2 i+1 )f (ζ i+1 -ζ i )θ i+1 and I (2) i+1 σ Y,i ρ i √ 1-ρ 2 i ( 1 -ρ 2 i Z 1 i+1 -ρ i Z 2 i+1 )f (ζ i+1 - ζ i )θ i+1 belong to M i,n ( X, Ȳ , -1/2).
From the preceding arguments, we eventually deduce that f

(ζ i+1 -ζ i ) - → θ c i+1 ∈ M i,n ( X, Ȳ , -1/2) for any i ∈ {0, • • • , n -1}.
Finally, from the very definition of the weights on the last time interval -→ θ e,Y n+1 and -→ θ e,X n+1 one directly gets that

(1 -F (T -ζ n )) - → θ e,Y n+1 = m n + σ Y,n ρ n Z 1 n+1 + 1 -ρ 2 n Z 2 n+1 + σ Y,n ρ n 1 -ρ 2 n 1 -ρ 2 n Z 1 n+1 -ρ n Z 2 n+1 belongs to M n,n ( X, Ȳ , 0) and that (1 -F (T -ζ n )) - → θ e,X n+1 = - 1 2 a S,n + σ S,n Z 1 n+1 belongs to M n,n ( X, Ȳ , 1/2
E T ζn G(ζ 1 , . . . , ζ n , s, T )1 {N T =n} ds = E G(ζ 1 , . . . , ζ n , ζ n+1 , T )(1 -F (T -ζ n+1 )) -1 (1 -F (T -ζ n ))(f (ζ n+1 -ζ n )) -1 1 {N T =n+1} .
Proof. The proof follows by rewriting the above expectations using (2.6). We rewrite the expectation on the right-hand side in integral form. By Fubini's theorem, we obtain

E G(ζ 1 , . . . , ζ n , ζ n+1 , T )(1 -F (T -ζ n+1 )) -1 (1 -F (T -ζ n ))(f (ζ n+1 -ζ n )) -1 1 {N T =n+1} = ∆n+1(T ) G(s 1 , • • • , s n+1 , T )(1 -F (T -s n+1 )) -1 (1 -F (T -s n ))(f (s n+1 -s n )) -1 × (1 -F (T -s n+1 )) n j=0 f (s j+1 -s j ) ds n+1 = ∆n(T ) T sn G(s 1 , • • • , s n+1 , T ) ds n+1 (1 -F (T -s n )) n-1 j=0 f (s j+1 -s j ) ds n .
This completes the proof.

Lemma B.2. Let n ∈ N. On the set {N T = n}, the sequence of weights (θ i ) 1≤i≤n+1 defined by (3.3) and (3.4) satisfy:

∀i ∈ {1, • • • , n} , f (ζ i -ζ i-1 )θ i ∈ M i-1,n ( X, Ȳ , -1/2), (1 -F (T -ζ n ))θ n+1 ∈ M n,n ( X, Ȳ , 0). (B.1)
Proof. We investigate each term appearing in the definition of θ i ∈ S i-1,n ( X, Ȳ ) and seek to apply Lemma 2.3. From the Lipschitz property of a S and the space-time inequality (1.3), for any c > 0 and any c > c, the map

(x i-1 , y i-1 , x i , y i , s n+1 ) → c i S (x i-1 , y i-1 , x i , y i , s n+1 ) satisfies |c i S (x i-1 , y i-1 , x i , y i , s n+1 )| p qc (s i -s i-1 , x i-1 , y i-1 , x i , y i ) ≤ C|y i -m i-1 (y i-1 )| p qc (s i -s i-1 , x i-1 , y i-1 , x i , y i ) ≤ C(s i -s i-1 ) p/2 qc (s i -s i-1 , x i-1 , y i-1 , x i , y i ) so that, the random variables c i S ∈ M i-1,n ( X, Ȳ , 1/2), for any i ∈ {1, • • • , n + 1}.
Moreover, from the boundedness of the first and second derivatives of a S , it follows

D (1) i c i S = ∂ Ȳi c i S ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n+1 ) = ∂ 4 c i S ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n+1 ) = 1 2 a S ( Ȳi ) ∈ M i-1,n ( X, Ȳ , 0) and D
(1,1) i c i S = ∂ 2 Ȳi c i S ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n+1 ) = ∂ 2 4 c i S ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n+1 ) = 1 2 a S ( Ȳi ) ∈ M i-1,n ( X, Ȳ , 0).
From Lemma 2.3, we thus conclude

I (1) i (c i S ) ∈ M i,n ( X, Ȳ , 0) and I (1,1) i (c i S ) ∈ M i,n ( X, Ȳ , -1/2), i ∈ {1, • • • , n} .
In a completely analogous manner, omitting some technical details, we derive

I (2) i (b i Y ) ∈ M i,n ( X, Ȳ , 0), and
I (1,2) i (c i Y,S ), I (2,2) i (c i Y ) ∈ M i,n ( X, Ȳ , -1/2).
Hence, we obtain f

(ζ i -ζ i-1 )θ i ∈ M i-1,n ( X, Ȳ , -1/2), for any i ∈ {1, • • • , n}. We finally observe that (1 -F (T - ζ n ))θ n+1 = 1 ∈ M n,n ( X, Ȳ , 0). The proof is now complete.
Lemma B.3. Let T > 0 and n a positive integer. For any s n = (s 1 , • • • , s n ) ∈ ∆ n (T ), for any (x, y) ∈ R 2 , for any positive constant c there exist two positive constants C and c > c such that the transition density (t, x, y) → qc (t, x 0 , y 0 , x, y) defined by (2.5) satisfies the following semigroup property: (R 2 ) n qc (T -s n , x n , y n , x, y)×q c (s n -s n-1 , x n-1 , y n-1 , x n , y n )ו • •×q c (s 1 , x 0 , y 0 , x 1 , y 1 ) dx n dy n ≤ C n qc (T, x 0 , y 0 , x, y).

Proof. The dx 1 • • • dx n integrals are treated using the standard semigroup property of Gaussian kernels so that from the very definition of qc , it directly follows for some positive constant C ≥ 1 and then the semigroup property satisfied by Gaussian kernels. This completes the proof.

Appendix C. Some useful formulas

We here provide some useful formulas in order to device the unbiased Monte Carlo algorithms based on Theorem 3.1 and Theorem 4.1. Their proofs follow from standard computations as those used in subsection 2.4 and are omitted.

The following formulae are required in order to compute the weights (θ i ) 1≤i≤N T appearing in the identity (3.2). Note that in our examples since a Y (.) is constant, one has c i Y (.) ≡ 0 for i ∈ {1, • • • , N T }. Hence, for i ∈ {1, • • • , N T }, one has:

I (1) i (c i S ) = c i S I (1)
i (1),

I (1,1) i (c i S ) = c i S ((I (1)
i (1)) 2 -D

i I

i (1)), I

i (1)I

i (1) -I

i (1)D

(2)

i (c i Y,S ) -c i Y,S D (2)
i I

i (1).

The following formulae are needed in order to compute the weights for the Delta appearing in the identity (4.5), for i ∈ {1, • • • , N T } one has:

D (1) i (I (1,1) i (c i S )) = 2c i S I (1)
i (1)D

(1) i I

i (1), D

i (1), D

i (1), D

i (1)I

(2)

i (1) + c i Y,S D (1)
i I

i (1)I

i (1) -D

i I

i (1)D

i (c i Y,S ). The following formulae are required for the computation of the weights for the Vega appearing in the identity (4.6), for i ∈ {1, • • • , N T } it holds:

I (1,1) i+1 (d i+1 S ) = m i I (1,1)
i+1 (c i+1 S ), I

i (c i+1 Y,S ), I

i (c i+1 S ), D

i (I

(1,1) i (c i S )) = D (2)
i (c i S )(I

i (1)) 2 + 2c i S I

(1)

i (1)D

(2) i I

i (1) -D

i (c i S )D

(1) i I

i (1), D

(2)

i (I (1) i (c i S )) = D (2)
i (c i S )I

i (1) + c i S D

(2) i I

i (1), D

(2)

i (I (2) i (b i Y )) = D (2) i b i Y I (2)
i (1) + b i Y D

(2) i I

(2)

i (1) -b Y ( Ȳi), D (2)
i (I

(1,2) i (c i Y,S )) = D (2)
i (c i Y,S )I

i (1)I

(2)

i (1) + c i Y,S D (2) i I
(2)

i (1)I

i (1) + c i Y,S I

i (1)D

(2) i I

i (1) -I

(1)

i (1)D (2,2) i (c i Y,S ) -2D (2)
i (c i Y,S )D

i (1), D

i-1 (I

S )) = D (1,1) i (c i
i-1 (c i S )(I

i (1)) 2 + 2c i S I

i (1)D

i-1 I

i (1) -D

(2)

i-1 (c i S )D (1) i I (1)
i (1) -c i S D

i-1 (D

i (1)), D

i (1) + c i S D

(2)

i-1 I

(1)

i (1), D

i-1 (I

(2) i (b i Y )) = b i Y D (2)
i-1 I

(2)

i (1) + I

i (1)D

(2)

i-1 I

i (1) + I

i (1)D

(2)

i-1 I

(2)

i (1)) -I

i-1 (I

i (1))D

(2)

i (c i Y,S ) -c i Y,S D (2)
i-1 D

(2) i I

(1)

i (1) -D (2) i-1 (c i Y,S )D (2)
i I

i (1).

D

(2)

i θi =(f (ζi -ζi-1)) -1 D (2)
i (I

(1,1) i (c i S )) -D (2)
i (I

i (c i S )) + D

(2)

i (I (2)
i (b i Y )) + D

(2)

i (I (1,2) i (c i Y,S )) , D (2)
i-1 θi =(f (ζi -ζi-1)) -1 D

(2)

i-1 (I

(1,1) i (c i S )) -D (2)
i-1 (I

i (c i S )) + D

i-1 (I

i (b i Y )) + D

i θi + (f (ζi -ζi-1)) -1 D

(2)

i (D (2)
i-1 (c i Y,S ))I

i (1) -I

i (e X,i S ) = (f (ζi -ζi-1)) -1 (I

i (1)D

(2) i-1 (c i S ) -D

(1) i D

(2) i-1 (c i S )), D

-→ θ e,X i = (f (ζi -ζi-1)) -1 e X,i S D

(1) i I

(1)

i (1), I

i (θiD

i-1 Xi) = (θiI

i (1) -D

i θi)D

(2)

i-1 Xi -D

i (D

i-1 Xi)θi,

I (2) i ( 1 -ρ 2 i-1 Z 1 i -ρi-1Z 2 i )θi = 1 -ρ 2 i-1 Z 1 i -ρi-1Z 2 i (I (2)
i (1)θi -D

(2)

i θi) + ρi-1θi σY,i-1 1 -ρ 2 i-1 , -→ θ c i =I (2)
i (m i-1 θi --→ θ e,Y i ) + σY,i-1 ρ i-1

1 -ρ 2 i-1

I (2) i ( 1 -ρ 2 i-1 Z 1 i -ρi-1Z 2 i )θi + I (1)
i (θiD

i-1 Xi) -I

i ( -→ θ e,X i

) + D (2)
i-1 θi,

-→ θ e,Y N T +1 =θN T +1 m N T + σY,N T ρ N T 1 -ρ 2 N T ( 1 -ρ 2 N T Z 1 N T +1 -ρN T Z 2 N T +1 ) , -→ θ e,X N T +1 =θN T +1 - 1 2 a S,N T + σ S,N T Z 1 N T +1 , -→ θ c N T +1 =0, I (2)
N T +1 ( -→ θ e,Y N T +1 ) = -→ θ e,Y N T +1 I (2)
N T +1 (1) -D (2)

N T +1 (1) + θN T +1 ρ N T ρN T 1 -ρ 2 N T , I (2)
N T +1 ( -→ θ e,X N T +1 ) = -→ θ e,X N T +1 I (2)
N T +1 (1) -D

(2)

N T +1 -→ θ e,X N T +1 = -→ θ e,X N T +1 I (2)
N T +1 (1) -θN T +1 σ S,N T σS,N T .

Table 9 .

 9 [-0.00499082, 0.00595295] -0.000755278 [-0.00588091, 0.00437035] 0.1 0.15 -0.0200101 [-0.0308873, -0.00913292] -0.0249364 [-0.0346885, -0.0151842] -0.0286496 [-0.0358769, -0.0214223] 0.2 0.25 -0.0246278 [-0.0366948, -0.0125608] -0.032211 [-0.0436935, -0.0207285] -0.0311689 [-0.0392428, -0.023095] 0.3 0.4 -0.0354025 [-0.04987, -0.0209349] -0.0422004 [-0.0518535, -0.0325472] -0.0413018 [-0.0489346, -0.0336691] 0.4 0.5 -0.0492556 [-0.0663201, -0.0321911] -0.0512594 [-0.0638074, -0.0387114] -0.0517876 [-0.0597881, -0.0437871] Comparison between the unbiased Monte Carlo estimation for the Vega of a digital Call option in the Stein-Stein type model for different values of the parameters σ 1 and σ 2 .

  00773549 [-0.00782648, -0.00764449] -0.00805159 [-0.00985368, -0.0062495] -0.00846248 [-0.0101504, -0.00677453] 0.2 0.25 -0.0156691 [-0.0158849, -0.0154532] -0.0161045 [-0.0194751, -0.0127339] -0.0137565 [-0.0169305, -0.0105825] 0.3 0.4 -0.0235822 [-0.0240098, -0.0231547] -0.0177797 [-0.0236385, -0.0119209] -0.0232616 [-0.0288379, -0.0176852] 0.4 0.5 -0.030774 [-0.0314484, -0.0300996] -0.031729 [-0.0405267, -0.0229313] -0.0327252 [-0.0402293, -0.0252211]

  similarly to step 2, by Fubini's theorem, we getE |h( Xn+1 , Ȳn+1 )| p n+1 i=1 |θ i | p N T = n, ζ n = R 2 |h(x, y)| p E p(T -ζ n , Xn , Ȳn , x, y) n+1 i=1 |θ i ( Xi-1 , Ȳi-1 , Xi , Ȳi , ζ n )| p N T = n, ζ n dxdy.The above formula is justified by Lemma B.2 and Lemma B.3 which yield E p(T -ζ n , Xn , Ȳn , x, y)

(R 2 2 2c≤

 22 ) n qc (T -s n , x n , y n , x, y) × qc (s n -s n-1 , x n-1 , y n-1 , x n , y n ) × • • • × qc (s 1 , x 0 , y 0 , x 1 , y 1 ) dx n dy n -s n ) e -(y-m T -sn (yn ))We now perform the change of variablesy 1 = m s1 (z 1 ), y 2 = m s2 (z 2 ), • • • , y n = m sn (z n ).Observe that since b Y admits a bounded first derivative the determinants of the Jacobians J s1 (z1 ) := ∂ x m s1 (z 1 ), • • • , J T -sn (z n ) = ∂ x m T -sn (z n ) are (locally) uniformly bounded for any (s 1 , • • • , s n ) ∈ ∆ n (T ).Remark also that from the semigroup property m si+1-si (m si (z i )) = m si+1 (z i ), for 1 ≤ i ≤ n with the convention s n+1 = T . Hence, for some positive constants C and c > c that may change from line to line, we get1 √ 2πcT e -(x-x 0 ) 2 2cT R n 1 2πc(T -s n ) e -(y-m T -sn (yn )) 2 2c(T -sn ) -s n ) e -(y-m T (zn )) 2 2c(T -sn ) × • • • × 1 √2πcs 1 e -(ms 1 (z 1 )-ms 1 (y 0 )C n qc (T x 0 , y 0 , x 1 , y 1 )where we first used the bi-Lipschitz property of the flow (s, x) → m s (x) which yields ∀t ∈ [0, T ], C -1 |x -z| 2 ≤ |m t (x) -m t (z)| 2 ≤ C|x -z| 2

D ( 1 )N

 1 i θi = (f (ζi -ζi-1)) -1 D T +1 (θN T +1) = θN T +1I (1) N T +1 (1) -D (1) N T +1 θN T +1 = θN T +1I (1)N T +1 (1).

i ( 1 )

 1 (b Y ( Ȳi)∂ Ȳi-1 Ȳi -b Y (mi-1)m i-1 ) -b Y ( Ȳi)∂ Ȳi-1 Ȳi, D

=

  m i-1 θi + (f (ζi -ζi-1)) -1 (I

i

  (m i-1 θi --→ θ e,Y i ) = -(f (ζi -ζi-1)) -1 I (1,2) i (D (2) i-1 (c i Y,S )) = -(f (ζi -ζi-1)) -1 D (2) i-1 (c i Y,S )I(1)i (1)I

  ζi -ζi-1)) -1 I

Table 1 .

 1 Comparison between the unbiased Monte Carlo estimation and the Monte Carlo Euler-Maruyama scheme for the price of a Call option in the Black-Scholes model for different values of σ S .

	σ S	B-S formula	Price	Euler Scheme 95% CI	Exponential sampling Price 95% CI	Price	Beta sampling 95% CI
	0.25 0.111804 0.111467 [0.110699, 0.112235] 0.112285 [0.111781, 0.112789] 0.112159 [0.111734, 0.112584]
	0.3 0.132621 0.13293 [0.132337, 0.133524] 0.133054 [0.132394, 0.133713] 0.132954 [0.132482, 0.133425]
	0.4 0.174152 0.17392 [0.173103, 0.174737] 0.175346 [0.174557, 0.176135] 0.174584 [0.173912, 0.175255]
	0.6 0.256572 0.258063 [0.256727, 0.259399] 0.255934 [0.254592, 0.257277] 0.256514 [0.255419, 0.257608]
	σ S	B-S formula	Delta	Euler Scheme 95% CI	Exponential sampling Delta 95% CI	Delta	Beta sampling 95% CI
	0.25 0.556589 0.555686 [0.553178, 0.558194] 0.554613 [0.551336, 0.557891] 0.556314 [0.553141, 0.559488]
	0.3 0.560018 0.561099 [0.559455, 0.562742] 0.557398 [0.553517, 0.56128] 0.557561 [0.554848, 0.560274]
	0.4 0.569512 0.570293 [0.568533, 0.572053] 0.569098 [0.565706, 0.572489] 0.56731 [0.564279, 0.570341]
	0.6 0.592743 0.594988 [0.592957, 0.59702] 0.586245 [0.582428, 0.590062] 0.588015 [0.584572, 0.591457]

Table 2 .

 2 Comparison between the unbiased Monte Carlo estimation and the Monte Carlo Euler-Maruyama scheme for the Delta of a Call option in the Black-Scholes model for different values of σ S .

	σ S	B-S formula	Vega	Exponential sampling 95% CI	Vega	Beta sampling 95% CI
	0.25	0	0.000745386 [-0.00102979, 0.00252057] -0.000438032 [-0.00211468, 0.00123862]
	0.3	0	-0.0013932 [-0.0036299, 0.000843502] -0.000491083 [-0.00249688, 0.00151471]
	0.4	0	0.00331309 [0.000258292, 0.00636788] -0.00117019 [-0.00393975, 0.00159938]
	0.6	0	-0.00286877 [-0.00777679, 0.00203925] -0.0027807 [-0.00718374, 0.00162235]

Table 3 .

 3 Comparison between the unbiased Monte Carlo estimation for the Vega of a Call option in the Black-Scholes model for different values of σ S .5.2.A Stein-Stein type model. In this second example, we consider a Stein-Stein type model where the diffusion coefficient function for the spot price is an affine function, namely σ S (x) = σ 1 x+σ 2 where σ 1 and σ 2 are two positive constants. Note carefully that σ S is not uniformly elliptic and bounded so that (AR) and (ND) are not satisfied. However, we heuristically choose σ 1 and σ 2 so that σ S (Y t ) is bounded and strictly positive with high probability. Also, analytical expressions for the coefficients are available, namely

Table 4 .

 4 Comparison between the unbiased Monte Carlo estimation for the price of a Call option in the Stein-Stein type model for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Delta	Euler Scheme 95% CI	Exponential sampling Delta 95% CI	Delta	Beta sampling 95% CI
	0.1 0.15 0.545257 [0.542601, 0.547914] 0.542838 [0.534227, 0.551448] 0.540304 [0.535519, 0.545088]
	0.2 0.25 0.548642 [0.54574, 0.551544] 0.541611 [0.533812, 0.549409] 0.535165 [0.51798, 0.552351]
	0.3 0.4 0.566919 [0.563629, 0.570208] 0.555688 [0.548827, 0.56255] 0.558808 [0.553545, 0.56407]
	0.4 0.5 0.579445 [0.575861, 0.583028] 0.569003 [0.559972, 0.578034] 0.568666 [0.561328, 0.576004]

Table 5 .

 5 Comparison between the unbiased Monte Carlo estimation for the Delta of a Call option in the Stein-Stein type model for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Vega	Euler Scheme 95% CI	Exponential sampling Vega 95% CI	Vega	Beta sampling 95% CI
	0.1 0.15 0.0370801 [0.0367679, 0.0373923] 0.0340152 [0.0317984, 0.036232] 0.0350342 [0.0333134, 0.036755]
	0.2 0.25 0.0738723 [0.0731769, 0.0745676] 0.0704958 [0.0662385, 0.074753] 0.0652897 [0.0612527, 0.0693267]
	0.3 0.4 0.11114 [0.109907, 0.112373] 0.0899367 [0.0830599, 0.0968136] 0.10303 [0.0912086, 0.114851]
	0.4 0.5 0.14385 [0.142055, 0.145645] 0.122496 [0.109106, 0.135885] 0.133235 [0.125356, 0.141114]

Table 6 .

 6 Comparison between the unbiased Monte Carlo estimation for the Vega of a Call option in the Stein-Stein type model for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Price	Euler Scheme 95% CI	Exponential sampling Price 95% CI	Price	Beta sampling 95% CI
	0 0.3 0.469652 [0.467241, 0.472063] 0.469251 [0.468668, 0.469834] 0.468994 [0.468439, 0.469548]
	0.1 0.15 0.491121 [0.488708, 0.493535] 0.489251 [0.488158, 0.490344] 0.489883 [0.489044, 0.490723]
	0.2 0.25 0.459518 [0.45711, 0.461926] 0.458577 [0.457324, 0.459829] 0.458555 [0.457574, 0.459535]
	0.3 0.4 0.430451 [0.428057, 0.432845] 0.428559 [0.42734, 0.429779] 0.428941 [0.428074, 0.429807
	0.4 0.5 0.408908 [0.406529, 0.411286] 0.40788 [0.404907, 0.410852] 0.409511 [0.408526, 0.410496]

Table 7 .

 7 Comparison between the unbiased Monte Carlo estimation for the price of a digital Call option in the Stein-Stein type model for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Delta	Euler Scheme 95% CI	Exponential sampling Delta 95% CI	Delta	Beta sampling 95% CI
	0 0.3 1.22307 [1.19252, 1.25363] 1.24579 [1.24326, 1.24833] 1.24408 [1.24165, 1.24651]
	0.1 0.15 2.17706 [2.13691, 2.21721] 2.17577 [2.16695, 2.18459] 2.18049 [2.17398, 2.18701]
	0.2 0.25 1.29839 [1.26695, 1.32984] 1.26832 [1.26267, 1.27397] 1.26854 [1.26428, 1.27279]
	0.3 0.4 0.776519 [0.752002, 0.801036] 0.792788 [0.789598, 0.795978] 0.793139 [0.790702, 0.795577]
	0.4 0.5 0.606688 [0.58496, 0.628416] 0.618031 [0.610424, 0.625638] 0.621753 [0.619061, 0.624446]

Table 8 .

 8 Comparison between the unbiased Monte Carlo estimation for the Delta of a digital Call option in the Stein-Stein type model for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Vega	Euler Scheme 95% CI	Vega	Exponential sampling 95% CI	Vega	Beta sampling 95% CI
	0 0.

Table 10 .

 10 Comparison between the unbiased Monte Carlo estimation for the price of a Call option in the model with σ S (x) = σ 1 cos(x) + σ 2 for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Delta	Euler Scheme 95% CI	Exponential sampling Delta 95% CI	Delta	Beta sampling 95% CI
	0.1 0.15 0.556917 [0.554118, 0.559717] 0.560077 [0.556733, 0.563422] 0.555364 [0.552636, 0.558092]
	0.2 0.25 0.577937 [0.574727, 0.581148] 0.579704 [0.575622, 0.583785] 0.577287 [0.574331, 0.580243]
	0.3 0.4 0.605788 [0.601947, 0.60963] 0.604575 [0.602771, 0.606379] 0.601188 [0.599354, 0.603021]
	0.4 0.5 0.62865 [0.624204, 0.633096] 0.623698 [0.618519, 0.628878] 0.626259 [0.622246, 0.630271]

Table 11 .

 11 Comparison between the unbiased Monte Carlo estimation for the Delta of a Call option in the model with σ S (x) = σ 1 cos(x) + σ 2 for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Vega	Euler Scheme 95% CI	Vega	Exponential sampling 95% CI	Vega	Beta sampling 95% CI
	0.1 0.15 -0.

Table 12 .

 12 Comparison between the unbiased Monte Carlo estimation for the Vega of a Call option in the model with σ S (x) = σ 1 cos(x) + σ 2 for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Price	Euler Scheme 95% CI	Exponential sampling Price 95% CI	Price	Beta sampling 95% CI
	0 0.3 0.470206 [0.467795, 0.472617] 0.468756 [0.468174, 0.469338] 0.46883 [0.468273, 0.469387]
	0.1 0.15 0.481972 [0.479559, 0.484385] 0.481373 [0.480778, 0.481968] 0.481499 [0.480937, 0.482061]
	0.2 0.25 0.446163 [0.443761, 0.448566] 0.445241 [0.444661, 0.445821] 0.445228 [0.444679, 0.445778]
	0.3 0.4 0.407307 [0.40493, 0.409684] 0.407842 [0.407286, 0.408398] 0.407422 [0.40689, 0.407954]
	0.4 0.5 0.379459 [0.37711, 0.381808] 0.380161 [0.379604, 0.380717] 0.379669 [0.379146, 0.380193]

Table 13 .

 13 Comparison between the unbiased Monte Carlo estimation for the price of a digital Call option in the model with σ S (x) = σ 1 cos(x) + σ 2 for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Delta	Euler Scheme 95% CI	Exponential sampling Delta 95% CI	Delta	Beta sampling 95% CI
	0 0.3 1.24454 [1.21372, 1.27535] 1.24282 [1.24031, 1.24533] 1.24579 [1.24335, 1.24822]
	0.1 0.15 1.5064 [1.47264, 1.54016] 1.51062 [1.50751, 1.51373] 1.51253 [1.50958, 1.51549]
	0.2 0.25 0.846474 [0.820904, 0.872044] 0.835956 [0.834251, 0.83766] 0.834907 [0.833284, 0.836531]
	0.3 0.4 0.512589 [0.492588, 0.53259] 0.528498 [0.527433, 0.529563] 0.527844 [0.526825, 0.528863]
	0.4 0.5 0.399506 [0.381818, 0.417194] 0.40319 [0.402367, 0.404012] 0.402428 [0.401642, 0.403214]

Table 14 .

 14 Comparison between the unbiased Monte Carlo estimation for the Delta of a digitalCall option in the model with σ S (x) = σ 1 cos(x) + σ 2 for different values of the parameters σ 1 and σ 2 .

	σ 1 σ 2	Vega	Euler Scheme 95% CI	Vega	Exponential sampling 95% CI	Vega	Beta sampling 95% CI
	0 0.3	0	[0, 0]	0.00122669 [-0.00426466, 0.00671804] -0.00228675 [-0.00742253, 0.00284903]
	0.1 0.15 0.00769619 [-0.000285705, 0.0156781] 0.00900717 [0.00339831, 0.014616] 0.00730275 [0.0021288, 0.0124767]
	0.2 0.25 0.0138531 [0.00480268, 0.0229036] 0.0146584 [0.00906018, 0.0202566] 0.0131819 [0.00819735, 0.0181665]
	0.3 0.4 0.0107747 [0.00279286, 0.0187565] 0.00807909 [0.00292412, 0.0132341] 0.0116897 [0.00689717, 0.0164823]
	0.4 0.5 0.0153924 [0.00585238, 0.0249324] 0.0152859 [0.0102097, 0.0203621] 0.0164134 [0.0117414, 0.0210853]

Table 15 .

 15 Comparison between the unbiased Monte Carlo estimation for the Vega of a digitalCall option in the model with σ S (x) = σ 1 cos(x) + σ 2 for different values of the parameters σ 1 and σ 2 .

  ζ n dxdy which can be justified as follows. From Lemma B.2, Lemma B.3 and the upper-bound estimate (2.4), it holds E p(T -ζ n , Xn , Ȳn , x, y)

  Emergence of jumps in the renewal process N . The first result is used in the proof of the probabilistic representation in Theorem 3.1 and is used to express that time integrals add jumps to the renewal process N . In what follows, N is a renewal process in the sense of Definition 2.1.Lemma B.1. Let n ∈ N and G : {(t 1 , . . . , t n+2 ) : 0 < t 1 < • • • < t n+1 < t n+2 := T } → Rbe a measurable function such that E T ζn |G(ζ 1 , . . . , ζ n , s, T )|1 {N T =n} ds < ∞. Then, it holds

). The proof is now complete. Appendix B. Some technical results B.1.

This dynamical system is obtained by removing the noise, that is, by setting σ Y ≡ 0, from the dynamics of Y in (1.1).

As before, we use the convention∅ • • • = 0, ∅ • • • = 1.

respectively for the Call option and in Table

13, Table14, Table15for the digital Call option. Here again, the unbiased estimators perform very well for all range of values of the parameters.

Appendix A. Proof of Theorem 3.1 and Lemma 4.1 A.1. Proof of Theorem 3.1. The proof is divided into three steps. In the first part, we establish the probabilistic representation for a bounded and continuous function h. We then provide the extension to measurable maps satisfying the growth condition 3.1. We eventually conclude by establishing the L p -moments when the jump times are distributed according to the Beta law.

Denote by L and ( Lt ) t≥0 the infinitesimal generators of (P t ) t≥0 and ( Pt ) t≥0 respectively given by