Integration by parts formula for killed processes: A point of view from approximation theory

In this paper, we establish a probabilistic representation for two integration by parts formulas, one being of Bismut-Elworthy-Li's type, for the marginal law of a one-dimensional diffusion process killed at a given level. These formulas are established by combining a Markovian perturbation argument with a tailor-made Malliavin calculus for the underlying Markov chain structure involved in the probabilistic representation of the original marginal law. Among other applications, an unbiased Monte Carlo path simulation method for both integration by parts formula stems from the previous probabilistic representations.


Introduction
In this article, we consider the following one-dimensional stochastic differential equation (SDE in short) where the coefficients b, σ : R Ñ R are smooth and bounded functions and pW t q tě0 stands for a onedimensional Brownian motion on a given filtered probability space pΩ, F , pF t q tě0 , Pq. The aim of this paper is to provide a probabilistic representation for two integration by parts (IBP) formulas for the marginal law of the process X killed at a fixed given level L. To be more specific, for a starting point x ě L, let τ " inf tt ě 0 : X t ă Lu be the first hitting time of the level L by the onedimensional process X. For a given finite horizon T ą 0, we are interested in establishing probabilistic representations for IBP formulae related to the following quantities Erf 1 pX T q1 tτ ąT u s and B x Erf pX T q1 tτ ąT u s (2) where f is a real-valued smooth function defined on rL, 8q. Extensions to non-smooth functions or to the transition density of the killed process are also obtained.
In the recent past, IBP formulae have raised a lot of interest as these explicit formulae can be further analyzed to obtain properties of densities or for Monte Carlo simulation among other applications, see Nualart [27] or Malliavin and Thalmaier [25]. The former quantity in (2) is commonly considered in the literature on Malliavin calculus, while the latter quantity is referred in the literature to as the Bismut-Elworthy-Li (BEL for short) formula. The BEL formula is of interest for many practical applications such as the numerical computation of price sensitivities in finance for Delta hedging purpose. For a more detailed discussion on this topic, we refer the interested reader to Fournié and al. [13], Fournié and al. [12], Gobet et al. [18] [8] for a short sample.
In particular, in Section 2.6 of [25], the authors propose a continuous time version of the IBP formula for d-dimensional diffusion process X killed when it exits an open sub-domain D of R d . Denoting by τ the first exit time of X from D, the formula writes B x Erf pX T q1 tτ ąT u s " Erf pX T q1 tτ ąT u Hs where H has an explicit expression using stochastic integrals. Extensions have also been proposed by Arnaudon and Thalmaier [5], [6]. At this stage, it is important to observe that, from a numerical perspective, these formulae will inevitably involve a time discretization, thus introducing a bias, in order to devise a Monte Carlo simulation method, as it is already the case for the quantity Erf pX T q1 tτ ąT u s, see e.g. Gobet [17], Gobet and Menozzi [19]. As observed in [6], it may also require to compute the solution of a control problem which can be done off-line. One may thus claim as stated at the beginning of Section 2.6 in [25]: "Its Monte-Carlo implementation, which has not yet been done, seems to be relatively expensive in computing time".
Our approach is probabilistic and relies on a perturbation argument of Markov semigroups to derive a probabilistic representation for the marginal law of the killed process based on a simple Markov chain approximation scheme for which we develop an appropriate Malliavin calculus machinery.
The main novelty in comparison with the aforementioned previous works lies in the fact that an unbiased Monte Carlo simulation method directly stems from the integration by parts formulae derived here. One may thus devise an estimator which does not involve any bias but only a statistical error. To the best of our knowledge, this feature appears to be new. As a by product of our analysis, we propose a probabilistic representation for the two derivatives B x ppT, x, zq and B z ppT, x, zq where p0, 8qr L, 8q 2 Q pT, x, zq Þ Ñ ppT, x, zq stands for the transition density evaluated at terminal point z at time T of the process X starting from x and killed at level L. We also point out that devising a Monte Carlo estimator for IBP formulae without introducing a bias from the exact simulation methods of Jenkins [23] or Herrmann and Zucca [21] does not seem to be apparent. An extension of the exact simulation method introduced by Beskos and al. [9] to compute the two first derivatives with respect to the starting point x of the quantity Erf pX x T qs, X being a one-dimensional diffusion with constant diffusion coefficient, has been recently proposed by Tanré and Reutenauer [28]. However, it seems difficult to implement from this method a Monte Carlo estimator for the aforementioned derivatives of the transition density of the diffusion without introducing any bias.
The first step towards obtaining an IBP formula is to prove a probabilistic representation for the marginal law of the killed diffusion process, in the spirit of Bally and Kohatsu-Higa [7] which developed such a formula for multi-dimensional diffusion processes (without stopping) and some Lévy driven SDEs by means of a probabilistic perturbation argument for Markov semigroups. We also refer the reader to Labordère et al. [20] and Agarwal and Gobet [1] for some recent contributions in that direction for multi-dimensional diffusion processes.
Such representation involves a simple Markov chain structure evolving along a time grid given by the jump times of an independent renewal process. A similar representation was derived by the same authors in [14] by means of analytic arguments. However, the representation obtained here is different and more amenable for the implementation of Monte Carlo simulation methods or to establish IBP formulae.
Once such representation is established, we then want to prove suitable IBP formulae using the underlying Markov chain structure and eventually the noise provided by the jump times. For this purpose, we set up a tailor-made Malliavin calculus for this new approximation process and perform a careful propagation argument of spatial derivatives, backward in time for the first quantity in (2) and forward in time for the second quantity in (2).
These developments are not free of mathematical hurdles. In fact, the proposed methodology leads to the appearance of boundary terms which have to be treated carefully. The key idea that we develop to deal with this issue consists in using the noise provided by the jump times. Contrary to the IBP formulae developed here, we point out that in most cases the explicitness of the IBP formulae for diffusions killed at a boundary demands a number of simplifications and approximations. A technical argument commonly used consists in performing a localisation of the underlying process in order to ensure that it is not close to the boundary. This technique is successful but has some theoretical and practical limitations as shown in Delarue [11], [18] and Nakatsu [26]. This is one of the main reasons for the previously quoted statement in Section 2.6 of [25].
We finally emphasize that the variance of the Monte Carlo estimators associated to the IBP formulae established here tends to be large and even infinite. This feature is not new and appears to be reminiscent of the probabilistic representation originally obtained in [7]. Importance sampling or higher order methods have been proposed to circumvent this issue in the case of multi-dimensional diffusions, see Andersson and Kohatsu-Higa [3] and more recently Andersson et al. [4]. Following the ideas developed in [4], we show how to achieve finite variance for the Monte Carlo estimators obtained from the probabilistic representation formulas of the marginal law of the killed process and of both IBP formulas by employing an importance sampling scheme on 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 provide some basic definitions, assumptions and present a reflection principle based on a simple one step Markov chain that will play a central role in our probabilistic representations for the marginal law of the killed process and for our IBP formulae. In addition, we also construct the adequate Malliavin calculus machinery related to the underlying Markov chain upon which both IBP formulae are made. In Section 3, the probabilistic representation for the marginal law of the killed process, based on the Markov chain of Section 2, is established. The change from the process X to the Markov chain coming from the reflection principle of Section 2 simplifies our analysis as all the irregularities of the process appear as indicator functions. Section 4 is devoted to the main ingredients to obtain our first IBP formula, that is, the IBP formula for the quantity Erf 1 pX T q1 tτ ąT u s.
These ingredients are put to work in Section 5. Theorem 5 is the main result of this section. As a by product, we obtain a probabilistic representation for the first derivative of the transition density of the killed process with respect to its terminal point in Corollary 1. In Section 6, we establish the BEL formula for the law of the killed process. The main result of this section is Theorem 7 and, as a by product, we obtain a probabilistic representation for the first derivative of the transition density of the killed process with respect to the initial condition x in Corollary 2. Many of the proofs of Sections 3 and 4 are technical and postponed to the appendix in Section 10. In Section 7, we show how to achieve finite variance for our unbiased Monte Carlo estimators by an importance sampling technique that we briefly present. Some numerical results are presented in Section 8. Clearly, one needs to study numerical issues in more detail and these are left for later studies.
Notations. For a fixed given point z P R, the Dirac measure is denoted by δ z pdxq. Derivatives may be denoted by f 1 pxq in the one-dimensional case or by B i f pxq in the multi-dimensional case for the partial derivative with respect to the i-th variable appearing in the multivariate function f or also by B i x f pxq " f piq pxq where the latest is used in functions of one variable and the former is used mostly for multivariate functions.
We will often work with continuous or smooth functions defined on rL, 8q. In order to shorten notation, we will consider their extensions 1 on R. This is done just in order to shorten the length of equations and notation. Therefore all statements can be rewritten using the same assumptions but restricted to the domain rL, 8q.
The space of functions which are k-times continuously differentiable on a closed domain D is denoted by C k pDq. At the boundary points of D, all derivatives are considered as limits taken from the interior of D only. In the case that the derivatives are bounded, the space of corresponding functions is denoted by C k b pDq, k P N Y t0, 8u. In particular, note that functions in C 1 b may not be bounded but are at most linearly growing. Finally, C k p pDq denote the class of k-times differentiable functions with at most polynomial growth at infinity. The space of p-integrable random variables (r.v.'s) is denoted by L p with its extension L 8 for p " 8.
Given a measure space S, the space of real valued Borel measurable functions on S will be denoted by MpSq. We also introduce the simplex A n :" tt P p0, T s n ; 0 ă t 1 ă¨¨¨ă t n ď T u, n P N where T is fixed throughout the paper.
The transition density function at x of the standard Brownian motion at time t is denoted by gpt, xq " p2πtq´1 {2 expp´x 2 {p2tqq. Its associated Hermite polynomials of order i, i P N, are defined by H i pt, xq " pgpt, xqq´1B i x gpt, xq. In order to simplify lengthy equation, we may use the symbol E " to mean that two quantities are equal in expectation. Sometimes the same symbol maybe used for equality on conditional expectation. This will be clearly indicated at the point where it is used.
Generic constants are usually denoted by C and are independent of all variables unless otherwise explicitly stated. As usual they may change value from one line to the next.
As we will be using discrete Markov chains in this article, we will often have indexes whose range may be the set of integers. In order to shorten the length of statements, we will use the following notation for the most common set of indexes N n " t1, ..., nu orN n " t0, ..., nu with n PN " N Y t0u. In the case that n ď 0, then N n :" H.

Preliminaries
2.1. Assumptions and basic definitions. Throughout the article, we work on a probability space pΩ, F , Pq which is assumed to be rich enough to support all r.v's considered in what follows. In addition, we will work under the following assumptions on the coefficients: (i) The coefficients of the SDE (1) are smooth and bounded, in particular, b P C 8 b pRq and a P C 8 b pRq. (ii) The function σ is bounded and uniformly elliptic, that is, there exist a, a ą 0 such that for any x P R, a ď apxq " σ 2 pxq ď a. Therefore, without loss of generality, we will assume that σpxq ą 0.

2.2.
A reflection principle. Our probabilistic representation involves the following approximation pro-cessX where ρ is a r.v. distributed according to a Bernoulli(1/2) law, independent of W , namely Ppρ " 1q " Ppρ " 0q " 1{2. Under our assumption on the coefficients, the flow derivatives ofX s,x t exist. In particular, one has In the particular case that s " 0, we may use the simplified formX t "X 0,x t . At this point, we give a brief explanation about how the approximation processX appears in the forthcoming probabilistic representation. The proof of the following lemma is straightforward by using the reflection principle, see e.g. Karatzas and Shreve [24]. Lemma 1. Define the following approximation process: together with its associated exit timeτ :" inf t,Ȳ t " L ( . Then, for any bounded measurable function f , the following property is satisfied: 2.3. Basic Markov chain framework. We also consider a Poisson process with parameter λ ą 0, independent of the one-dimensional Brownian motion W with jump times T i , i P N and we set ζ i :" T i^T , i P N with the convention that ζ 0 " T 0 " 0. Define π to be the partition of r0, T s given by π :" t0 ": ζ 0 ă¨¨¨ă ζ NT ď T u. Associated with this set, we recall the definition of the simplex A n :" tt P p0, T s n ; 0 ă t 1 ă¨¨¨ă t n ď T u. For instance, on the set tN T " nu, n P N, we have pζ 1 , . . . , ζ n q P A n and ζ n`1 " T . In particular, for the set tN T " 0u (i.e. n " 0), we let π :" t0, T u and A 0 " H. In this sense, we will use throughout the rest of the paper the index n PN without any further mention of its range of values.
As it is the case in the previous observation, many proofs and definitions will be carried out conditioning on the set tN T " nu, n PN.
LetX :" pX i q iPN be the discrete time Markov chain starting at time 0 fromX 0 " x and evolving according toX where for simplicity we set σ i :" σpX i q, Z i`1 :" W ζi`1´Wζi " σ´1 i`X i`1´ρi`1Xi´p 1´ρ i`1 qp2L´X i qȃ nd tρ i ; i P Nu is an i.i.d. sequence of Bernoullip1{2q r.v.'s such that W, N and tρ i ; i P Nu are mutually independent. In what follows we use the notation h i " hpX i q, i PN n`1 for any function h : R Ñ R.
In particular, the reader may have noticed that we already used this notation in the above formula for h " σ. We also associate to the Markov chainX the following sets (7) D i,n :" tX i ě L, N T " nu, for i PN n`1 .
It is important to point out that given N T " n, the conditional distribution ofX n`1 givenX n is not the same as the conditional distribution ofX i`1 givenX i for i PN n´1 . This is due to the fact that the length of the last interval is T´ζ n , rather than the waiting time between two consecutive Poisson jumps. This remark applies to various definitions and results to be stated through the rest of the article.
We define the filtration G :" pG i q iPN where G i :" σpZ i , ζ i , ρ i q with the notation a i :" pa 1 , . . . , a i q for a " Z, ζ, ρ, i P N and G 0 defined as the trivial σ-field. We assume that the filtration G satisfies the usual conditions.
2.4. Simplified Malliavin Calculus for the underlying Markov chain. In this section we introduce the required material for our Malliavin calculus computations. Instead of using an infinite dimensional calculus as it is usually done in the literature, 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 Poisson process involved in the Markov chainX. In what follows, n PN unless stated otherwise.
We start by defining the following space of smooth r.v.'s.

Definition 1.
For i PN n , we define the set S i`1,n pXq as the subset of r.v.'s H P L 0 such that there exists a measurable function h : R 2ˆt 0, 1uˆA 2 Ñ R satisfying (2) For any r P t0, 1u and any ps, tq P A 2 , hp¨,¨, r, s, tq P C 8 p pR 2 q.
For a r.v. H P S i`1,n pXq, i PN n , we may sometimes abuse the notation and write that is the same symbol H may denote the r.v. or the function in the set S i`1,n pXq. One can easily define the flow derivatives for H P S i`1,n pXq as follows: We now define the derivative and integral operators for H P S i`1,n pXq, i PN n , as Note that due to the above definitions and Assumption pHq, we also have that I i`1 pHq, D i`1 H P S i`1,n pXq so that we can define iterations of the above operators, namely I ℓ`1 i`1 pHq " I i`1 pI ℓ i`1 pHqq and similarly D ℓ`1 i`1 pHq " D i`1 pD ℓ i`1 Hq, ℓ PN, with the convention I 0 i`1 pHq " D 0 i`1 pHq " H. Through this article, we will use the following notation for a certain type of conditional expectation that will appear frequently. For any X P L 1 and any i PN n , E i,n rXs :" ErX | G i , T n`1 , ρ n`1 , N T " ns.
With the above definitions, the following duality 2 is satisfied for any f P C 1 p pRq and any pi, ℓq PN nˆN : In order to obtain explicit norm estimates for r.v.'s in S i`1,n pXq it is useful to define for H P S i`1,n pXq, i PN n and p ě 1 }H} p p,i,n :" E i,n r|H| p s .
2 This duality is obtained using the Gaussian density ofX i`1 while in classical Malliavin calculus it is based on the density of the Wiener process. Therefore the derivative and integral, D i`1 and I i`1 , i PNn defined in the formula (11) are renormalizations of the usual duality principle in Malliavin calculus. In our case this notation simplifies greatly many equations.
Another useful formula that is used at several places later on is the following extraction formula for H 1 , H 2 P S i`1,n pXq : 1 Di,n |θ i | p ‰`1 ti"n`1u 1 Dn,n E n,n " 1 Dn`1,n |θ n`1 | p ‰ ď C.
As a consequence, for all p P r0, 2q, one has E "ˇˇś We importantly observe that the above probabilistic representation allows to implement an unbiased Monte-Carlo simulation method since one just has to simulate the Poisson process N , then the Markov chainX along the jump times of N and finally to compute the explicit product of weights. Contrary to our previous work [14], the Markov chainX is defined from the reflection principle introduced in Lemma 1 thus reducing further the number of variables in the problem. This probabilistic representation is thus more workable for the forthcoming IBP formulae.
Remark 3. We make several remarks before moving on: (i) The assumption f pLq " 0 can be avoided at the cost of longer formulae as one can obtain a probabilistic representation based on the same Markov chain for the probability Ppτ ě T q using the results in Section 3 of [14], but we do not pursue this goal here.
(ii) In the proof of the above result, one uses the following crucial property: For i PN n , the random variables 1 Di`1,n c i`1 1 , 1 Di`1,n BX i c i`1 2 P M i`1,n pX, 0q and 1 Di`1,n c i`1 2 P M i`1,n pX, 1{2q so that 1 Di`1,n I i`1 pc i`1 1 q, 1 Di`1,n I 2 i`1 pc i`1 2 q P M i`1,n pX,´1{2q and therefore 1 Di`1,nθi`1 P M i`1,n pX,´1{2q. (iii) The power p{2 appearing in the time degeneracy estimates in (19) is crucial in order to determine the integrability of the r.v. appearing on the right hand side of (16). This motivates the definition below.  20) in the case that i P N n and 1 Dn,n › › 1 Dn`1,n H › › p,n,n ď C in the case that i " n`1. At this stage, we find it useful to show graphically the dynamic structure of the Markov chain and the random weightsθ i , i P N NT`1 for the probabilistic representation (16). This will be important in order to understand the structure of the IBP formula.
In the above figure, one observes the evolution of the Markov chainX i , i PN n`1 , together with its associated weight sequence appearing in the probabilistic representation. We also note that at the end of the Markov chain evolution which always happens at time T , the test function f is evaluated atX n`1 . Furthermore, the right hand side of (16) is the product of all elements on top of the arrows on the first line of the above figure and the corresponding indicator functions of the sets D i,n defined by (7), with f pX n`1 q. The second line gives the time evolution followed by the jump times of the Poisson process N which coincide with the times at which transitions of the Markov chainX happen.
In other figures that appear later on, we will use the general symbol ‹ k , k P N n`1 defined in Figure  1 which stands for product ‹ k :" ś k i"1 1 Di,nθi . We remark here that as stated in (17), on the set tN T " nu, the definition of pX n`1 ,θ n`1 q differs from all other pX i ,θ i q, i P N n .

Ingredients for an IBP formula
In this section we give the main ingredients in order to establish an integration by parts formula for the marginal law taken at time T of the killed process.
In order to understand the main ingredients to be introduced in this section, one starts by supposing that the problem is to find an IBP formula for E " f 1 pX T q1 tτ ěT u ‰ . The BEL formula will be easier to handle and is tackled in Section 6. 4.1. Some Heuristics. The following heuristic arguments may help the reader to understand the strategy developed in the next sections to prove our IBP formulae.
Step 1: The first step was performed with the probabilistic representation established in Theorem 2 involving the Markov chainX evolving on a time grid governed by the set of jump times of the Poisson process N .
Step 2: As already explained in the introduction, the central idea is to reduce the infinite dimensional problem of finding an IBP formula for Erf 1 pX T q1 tτ ąT u s to a finite dimensional problem, namely finding an IBP formula for the quantity E " for which the dimension is random and given by the number of jumps of the Poisson process at time T .
At this stage, unfortunately, one cannot perform a standard integration by parts formula as in [27] on the whole time interval r0, T s for various reasons. For example, the Skorokhod integral of the product of weights ś NT`1 i"1 1 Di,N Tθ i will inevitably involve the Malliavin derivatives ofθ i and the indicator function 1 Di,N T , which in turn will raise integrability problems of the resulting Malliavin weight.
The key idea that we use consists in performing IBP formulae locally on each random intervals rζ i , ζ i`1 s, for i PN n on the set tN T " nu, that is, by using the noise of the Markov chain on this specific time interval and then by combining them in an suitable way. The case i " n is easy to consider because, by taking the conditional expectation E n,n r.s in the original probabilistic representation on the set tN T " nu, the IBP reduces to E n,n rf 1 pX NT`1 q1 Dn`1,nθn s " E n,n rf pX NT`1 q1 Dn`1,n I n`1 pθ n`1 qs by (11) and the fact that f pLq " 0. However, performing the IBP formula on a random interval rζ i , ζ i`1 s for 0 ď i ă n is more challenging.
In order to do it, our first ingredient consists in transferring the derivative operator appearing on the test function f backward in time from the last interval to the interval on which we perform the local IBP, say rζ i , ζ i`1 s. This operation will unfortunately generate a boundary term each time a transfer is performed, i.e. one for each time interval rζ j , ζ j`1 s, j " i,¨¨¨, n. As previously mentioned, the boundary terms will induce integrability problems. In order to circumvent this issue, we then introduce our second ingredient which consists in performing a boundary merging procedure of this term using the time randomness provided by the Poisson process. We refer the reader to Section 4.3 for a more detailed discussion of this issue. Then, an additional ingredient that we use is a time merging lemma about the reduction of jump times of the Poisson process. It is described in Lemma 13 and directly employed in the proof of our main result, namely Theorem 5.
As already explained above, the last ingredient consists in combining these various local IBP formulae in a suitable way. Roughly speaking, we consider a weighted sum of each integral operator, the weight being the length of the corresponding time interval. We remark that the last term in (21) is singular but well defined if the time variables are fixed which is the case here as there is a conditioning with respect to the Poisson process expressed with the conditional expectations E i,n in (21). Since the boundary weight i`1 deteriorates so that the resulting probabilistic representation will not be exploitable without the extra ingredient of the boundary merging lemmas. A further study is thus required before proving moment estimates. We will discuss this matter in detail in the next section.

4.3.
The boundary merging lemmas. The second ingredient that we need in order to establish our first IBP formula is the boundary merging lemma. Though the sequence of weights satisfies the time degeneracy estimate stated in Lemma 4, which is similar to the one satisfied by the sequence θ i , i P N NT ( in (19), the fact that it is multiplied by δ L pX i q increases its time degeneracy by a factor of pζ i´ζi´1 q´1 {2 .
The key idea in order to circumvent this issue consists in using the time randomness provided by the Poisson process in order to smooth this singularity. This will eventually allow us to retrieve a time degeneracy estimate similar to (19). We call this ingredient the boundary merging lemmas.
To be more specific, in order to circumvent the time degeneracy problem related to the boundary term on the right hand side of (21), one needs to consider two successive random time intervals, say rζ i´1 , ζ i s and rζ i , ζ i`1 s, and take expectations with respect to the intermediate time ζ i in order to regularize the term δ L pX i q that appears when one transfers the derivative using (21) on the time interval rζ i´1 , ζ i s. This operation will remove the time singularity induced by the Dirac delta distribution but will also reduce the number of jumps by one unit.
An important remark is that the boundary merging lemmas are not needed in the case that there are no jumps in the interval r0, T s, that is, on the set tN T " 0u, because of the condition f pLq " 0.
From the application of the boundary merging lemmas, a new Markov chain structure appears which we call the merged boundary (one-step) Markov chainX B whose transition on the time interval rζ j , ζ i`1 s, for j " i´1, i, is given bȳ with µpxq :"´σpLq{σpxq.
In the above one-step dynamics and in what follows, the index j " i´1 will be used to indicate the boundary merging of the underlying Markov chainX on the two consecutive time intervals rζ i´1 , ζ i s and rζ i , ζ i`1 s while the index j " i will be used in a second stage once the reduction of jumps (a.k.a. time merging) is performed in Section 5.1. As the proof of time degeneration estimates are similar and in order to simplify the presentation, we will deal with both cases at the same time in the next two Lemmas.
In the same way as it was done in Section 2.4, we define the associated space of smooth r.v.'s.

Definition 4.
For i PN n and j P ti´1, iu, we define the set S j,i`1,n pX B q as the subset of r.v.'s H P L 0 such that there exists a measurable functions h : (2) For any r P t0, 1u and any s P A i`1´j , one has hp¨,¨, r, sq P C 8 p pR 2 q. Similarly to the operators D and I, defined in Section 2.4, related to the Markov chainX, we introduce the operatorsD andĪ associated to the merged boundary processX B defined above. For any (smooth) r.v. H P S j,i`1,n pX B q and ℓ ě 1, we let j,i`1´p L`1´µpX j q˘`X j µpX j qqq P S j,i`1,n pX B q, it is clear that the r.v.'sĪ j,i`1 p1q andĪ 2 j,i`1 p1q also belong to S j,i`1,n pX B q so that they are explicit functions of the variablesX j ,X B j,i`1 , ρ i`1 , ζ j , ζ i and ζ i`1 . In the following result, we will also use new weights Ð Ý θ B˚e j,i`1 , for j " i´1, i, obtained by a merging procedure of the two weights for i " 0,¨¨¨, n´1. For the last interval, on tN T " nu, we have Ð Ý θ B˚e j,n`1 :" 4e λT a 1 pLq´bpLq aj , for j " n´1, n.
The cases j " i´1 for (27) and j " n´1 for Ð Ý θ B˚e n´1,n`1 are used in the following result where we perform the boundary merging procedure. The remaining cases j " i for Ð Ý θ B˚e i,i`1 and j " n for Ð Ý θ B˚e n,n`1 will be used once time merging is performed in Section 5. For this reason, some properties that we will employ later appear in the next lemma for these cases. The definition of the time degeneracy estimate in the sense of Definition 3 is naturally extended to this case if we define and replace (20) for H P S j,i`1,n pX B q by We are now in position to provide the second ingredient in order to establish our first IBP formula, namely the boundary merging lemma. Its proof is postponed to Section 10.3.2.
Lemma 5. Let f P C 0 p pRq and n PN. The following property is satisfied for any i P N n´1 Similarly, for the last time interval, one has Here, | Ð Ý θ B˚e j,n`1 | ď C a.s. for j " n´1, n. Moreover, with the above definitions Ð Ý θ B˚e j,i`1 P S j,i`1,n pX B q, j " i´1, i and Ð Ý θ B˚e j,i`1 satisfies the time degeneracy estimates in the sense of (28).
The above lemma can be illustrated using Figure 3 below. The boundary weight r.v. Ð Ý θ B i marked in blue is merged together with the weight Ð Ý θ e i`1 appearing in the next time interval, by taking expectations with respect to ζ i in (29). This operation leads to a new transition on the time interval rζ i´1 , ζ i`1 s for the underlying Markov chain. More precisely, the new r.v.X B i´1,i`1 , and the new weight, Ð Ý θ B˚e i´1,i`1 , both marked in red in the above figure, will replace the symbols in blue:  Figure 3, we see the result of the boundary merging procedure. Clearly the transition for the merged term has changed fromX i´1 toX B i´1,i`1 and then, for the next time interval rζ i`1 , ζ i`2 s, the transition of the underlying Markov chain remains unchanged and is given by but whose starting point isX B i´1,i`1 and its end point is Y i`2 . We remind the reader that we freely use function notation for smooth r.v.'s on the spaces S i`2,n pXq as in (8). Figure 3. Markov chain structure before merging of boundary terms on the left and after merging on the right.
So far, we have explained how to transfer the derivatives and how to deal with boundary terms. The last step consists in performing a local IBP formula on a fixed time interval, say rζ i , ζ i`1 s. This operation will involve the integral operator applied to corresponding weight, namely I i p1 Di,nθi q, and thus will inevitably increase the time singularity in the estimate as stated in (20). Moreover, by the extraction formula, the Malliavin derivative of 1 Di,n will also generate a boundary term that has to be carefully treated by a merging procedure. This is the purpose of the following second boundary merging Lemma that we now describe.
The proof of the following result is postponed to Section 10.3.3. As the resulting boundary merging weight is different, the notation for the boundary merged weight changes from B˚e to B f e. The merged boundary process does not change.
Lemma 6. Let f P C 0 p pRq and n PN. The following property is satisfied for any i P N n´1 Φpt, zq :" Similarly, for the last time interval, one has Moreover, for any i PN n , for any j " i´1, i, Ð Ý θ Bfe j,i`1 P S j,i`1,n pX B q and it satisfies the time degeneracy estimates in the sense of (28).
Remark 4. (i) We again emphasize the role of the two indexes j " i´1, i, for i P N n , in the definition of the boundary merged weights Ð Ý θ B˚e j,i`1 and Ð Ý θ Bfe j,i`1 of Lemmas 5 and 6. The index j " i´1 is used for the boundary merging operation while the index j " i will be used once the reduction of jumps operation (a.k.a. time merging) is performed in Section 5.1.
(ii) The reason why the termsΦg´1 appear in the definitions of the coefficientsd i`1 , contrary to the previous Lemma 5, is due to the time increments ζ i´ζi´1 and ζ n´ζn´1 on the left hand side of equalities (31) and (32).
(iii) An important technical remark which follows from the definitions of the merged weights Ð Ý θ B˚e j,i`1 and Ð Ý θ Bfe j,i`1 , j " i´1, i in Lemmas 5 and 6 is that they do not depend neither on ρ i or ρ i`1 .

A first IBP formula: Putting the ingredients to work
As explained at the beginning of Section 4, in order to obtain an IBP formula, one first has to take the conditional expectation w.r.t the Poisson process N inside E Once the jump times are fixed, one observes that there is a Markov chain structure to which one may apply the transfer of derivatives procedure described in Lemma 4. This leads to a tree structure of terms which are combined up to the time interval where we decide to stop the transfer of derivatives and to perform a local IBP formula using only the noise on this specific interval. Throughout this section we will often denote this interval by the general index k. For example, the tree structure is illustrated in Figure 4 where the IBP is performed on the time interval rζ k´1 , ζ k s for k " 2. The application of the IBP formula (11) on the time interval rζ k´1 , ζ k s will give after using the extraction formula (12) the new weight We thus see that a boundary term appears which is treated by the boundary merging procedure described in Lemma 6 and thus gives rise to the new weight associated to the symbol B f e.
The terms which contain boundary weights due to the successive application of the transfer of derivatives formula (21) are treated by the boundary merging procedure of Lemma 5. Figure 4 shows these terms before the boundary merging procedure in a simplified algebraic notation that will be introduced in this section.
This boundary merging operation will lead to the recovery of an integrable time degeneracy estimate and to the reduction of one jump unit (after an application of the so-called time merging procedure) in the underlying Poisson process. Finally, a local IBP formula is performed on the interval rζ k´1 , ζ k s where we have stopped the transfer of derivatives. As described above, after using the extraction formula, the new weight given by I k p1 D k,nθ k q will also generate a boundary term to which we apply the boundary merging procedure described by Lemma 6. As before, this will in turn reduce the number of jump times by one unit.
We importantly note that when one stops the (backward) transfer of derivatives procedure and decide to perform a local IBP formula on the time interval rζ k´1 , ζ k s, the weights 1 Dj,nθj for the preceding time intervals, namely rζ j´1 , ζ j s, for j " 1,¨¨¨, k´1, correspond to the original probabilistic representation in Theorem 2 and thus remain unchanged. A similar remark applies when a time merging of weights or a correction weight appears on one branch of the tree structure below.
The overall tree diagram before carrying out the boundary merging procedure can be schematically described as follows in the case of N T " 4 jumps of the Poisson process An important remark in order to understand the tree notation to be introduced in the next section is that the blue curved arrows will lead to the application of the boundary merging Lemma 5. Therefore, as the number of jumps will be reduced by the time merging procedure, those time merged weights will no longer be considered as weights associated to the set tN T " 4u but tN T " 3u. A similar remark applies when applying Lemma 6 to the boundary term generated by I k p1 D k,nθ k q which is denoted by the straight blue arrow in Figure 4.
We believe that the above explanations given before the proof of our main result are important for the reader because they are necessary to understand the following section where the various symbols for the above tree structure are introduced.
5.1. The tree structure in the IBP formula. In this section, we start by describing all the branches of the tree that are generated after performing transfer of derivatives and the two types of boundary merging previously described.
We denote by S n`1 , the set of all symbol sequences of length n`1 described by the following vectors of length n`1. More explicitly, when n " 0, we define S 1 :" tI 1 1 , B 1 1 u with I 1 1 :" pIq and B 1 1 :" pB f eq and for n PN, Therefore in S n`1 there are in total 2pn`1q`2n vectors. The above vectors give the description of a branch in the IBP formula and are defined as follows. On the one hand, we let C n`1 k :" p0, ..., 0, c, e, ..., eq and I n`1 k :" p0, ..., 0, I, e, ..., eq where the component c or I appears in the k-th coordinate for k " 2, ..., n`1 and we allow k " 1 only for the vectors I n`1 k . On the other hand, we denote by B n`1 k :" p0, ..., 0, Be , e, ..., eq, k " 2, . . . , n`1 and B n`1 k :" p0, ..., 0, B f e, e, ..., eq, k " 1, . . . , n`1, where the element B˚e or B f e appears in the k-th coordinate for k " 2, . . . , n and k " 1, ..., n respectively 3 .
As mentioned previously, the index k corresponds to the time interval on which we perform the local IBP formula. In other words, if the symbol 0 appears at the j-th coordinate of a vector, for 1 ď j ă k, this means that the weight corresponding to that time interval, namely rζ j´1 , ζ j s, remains the same as given in the probabilistic representation of Theorem 2, i.e. 1 Dj,nθj . The symbol c appearing at the k-th coordinate of the vector C n`1 k means that the new weight (associated to this vector and to the interval rζ k´1 , ζ k s) is 1 D k,n Ð Ý θ c k . The symbol I appearing at the k-th coordinate of I n`1 k stands for the IBP weight, i.e. 1 D k,n I k pθ k q. The symbol B˚e corresponds to the merging between exchange (denoted by e) and boundary (denoted by B) r.v.'s according to Lemma 5 while the symbol B f e corresponds to the same merging but occurring on the same interval as the one where we perform the IBP formula and is computed according to Lemma 6.
We remark here that in the last interval there is no boundary term because we always assume that the test function f vanishes at the boundary.
If we fix our attention on all the branches that are generated with the objective of carrying out the local IBP formula on the interval rζ k´1 , ζ k s, we find the two following subsets of S n`1 which correspond to the branches that finish with a correction weight Ð Ý θ c and the branches that finish with a time merging of two time intervals, before leaving unchanged the remaining weights on remaining intervals. For k " 1, . . . , n`1, we definē ( . For example, the following figure describes the situation in the case when originally there were n " 4 jump times. Observe that, on the one hand, due to the merging procedure, all sequences of arrows 3 In order to keep index notation short, so as not to have to always consider two cases when k takes as lowest value 1 or 2, we include C n`1 1 and B n`1 1 as symbols but any statement in this case should be taken as an empty statement or that the symbol corresponds to the zero (empty) element.  Figure 5. A tree in the case of N T " 4 jump times with an IBP on the interval rζ 1 , ζ 2 s after time merging. All branches that contain a merged terms are associated to the set tN T " 3u.

5.2.
The Markov chain and weights associated to the IBP tree branches. For the corresponding time partition π :" t0 " ζ 0 ă ... ă ζ n`1 " T u of the underlying Poisson process on the set tN T " nu, we now need to define the underlying Markov chainX s , for s P S n`1 or s P 9 S k n`1 , that will be used in the probabilistic representation of our first IBP formula.
It will be defined as the same process as the original Markov chainX except that on a certain time interval it may use the one-step transition of the merged boundary Markov chainX B given by (26). To be more specific, in the case that s " B n`1 k or B n`1 k , k " 1, ..., n`1, the Markov chainX s is defined bȳ k´1,k pX k´1 q for j " k, X j pk,X B k´1,k q for j " k`1, ..., n`1.
HereX B k´1,k pX k´1 q stands for the one-step transition defined in (26). Similarly, for j " k`1,¨¨¨, n`1, X j pk,X B k´1,k q is the flow notation for the scheme (6) taken at step j starting from the pointX B k´1,k at step k. For any other s P S n`1 of length n`1, the processX s corresponds to the original Markov chain dynamics, that is, we letX s "X. We also define the associated set D s i,n :" tX s i ě L, N T " nu.
We now introduce the weights corresponding to the Markov chainX s to be used in the IBP formula. For each s P S n`1 , we will define its associated weights p Ð Ý θ s 1 , ..., Ð Ý θ s n`1 q and its product 4 as We now proceed to define the weights for each element s " ps 1 , ..., s n`1 q P S n`1 as follows for i P N n`1 : As remarked at the beginning of Section 2.4, the weight Ð Ý θ s i is an explicit function of the underlying Markov chain and we will make use of such property in what follows. 4 Recall the standard convention ś H " 1.
In the case of time merging, we importantly refer the reader to Remark 4. Due to this, we know that the reduction of one jump does not affect the sequence ρ i , i P N n`1 as the boundary merged weight does not depend on it.
Theorem 5. Let f P C 1 b pRq such that f pLq " 0. Under assumption (H), the following IBP formula is satisfied: Moreover, the r.v. appearing inside the expectation of the right-hand side of the above equality belongs to L p pPq for p P r0, 2q.
Proof. We start from the Markov chain representation given by Theorem 2, namely where we notice that tN T " nu " tT n`1 ą T u X tT n ď T u. We remind the reader that E i,n rXs is the expectation of X conditional on tG i , T n`1 , ρ n`1 , N T " nu and therefore for the rest of the proof we will work on the set tN T " nu.
Step 1: IBP on the last interval. We start by proving the IBP for the last time interval. That is, by the tower property of conditional expectation and the integration by parts formula, (11) on the (deterministic) time interval rζ n , T s, noting that f pLq " 0, one has Erf 1 pX n`1 q1 Di,N T θ n`1 n ź i"1 1 Di,N Tθ i | T n`1 s " ErE n,n rD n`1 f pX n`1 q1 Dn`1,n θ n`1 s n ź i"1 1 Di,nθi | T n`1 s " Erf pX n`1 q1 Dn`1,n I n`1 pθ n`1 q n ź i"1 1 Di,nθi | T n`1 s.
Step 2: The transfer of derivatives. In this step, we will perform the transfer of derivatives from the last time interval rζ NT , T s to the time interval rζ k´1 , ζ k s. In order to carry this step, using the Markov property of the processX, we define for k P N n the functions: with the convention that F n`1 pX n`1 q " f pX n`1 q and ś H¨¨¨" 1. Note that the following recursive relation is satisfied for k P N n F k pX k q " E k,n rF k`1 pX k`1 q1 D k`1,n Ð Ý θ e k`1 s.
From the transfer of derivatives formula (25)
Next, we proceed using a backward induction argument by combining successive applications of the transfer of derivative formula of Lemma 4 with the tower property of conditional expectation. To be more specific, from Lemma 4, one has for k P N n 1 Di,nθi | T n`1 s for k " 1, . . . , n, with the convention ř H¨¨¨" 0.
Step 3: The local IBP on the interval rζ k´1 , ζ k s. To proceed, we first notice thatθ k P S k,n pXq, for k P N n`1 , is a smooth r.v.. Then from the tower property of conditional expectation (using E k´1,n r¨s), the integration by parts formula (11), the extraction formula (12) and (61) one obtains ErBX k F k pX k qθ k | G k´1 , T n`1 s " ErF k pX k q1 D k,n I k pθ k q | G k´1 , T n`1 s´ErF k pX k qδ L pX k qθ k | G k´1 , T n`1 s, Now, plugging the previous identity into (36) and using the recursive formula (34) yield for k P N n : 1 Di,nθi | T n`1 s (37)`E rf pX n`1 q1 Dn`1,n Ð Ý θ c n`1 n ź i"1 1 Di,nθi | T n`1 s.
Step 4: Combining all the local IBP formulae. Multiplying (33) by pT´ζ n q, (37) by pζ k´ζk´1 q and summing from k " 1 to n the resulting equalities, we obtain T Erf 1 pX n`1 q n`1 ź i"1 1 Di,nθi | T n`1 s " n`1 ÿ k"1 pζ k´ζk´1 qErBX n`1 f pX n`1 q n`1 ź i"1 1 Di,nθi | T n`1 s, "Erf pX n`1 q1 Dn`1,n pT´ζ n qI n`1 pθ n`1 q where we used the fact that T´ζ n`ř n k"1 pζ k´ζk´1 q " T´ζ 0 " T in the first equality. The final argument starts by multiplying the above identity by 1 tNT "nu and taking the expectation of both hand sides with respect to the Poisson process.
Step 5: The merging procedure. Now, we perform the boundary and time merging procedures for the fourth term appearing in the right-hand side of the above equality. This is done by first conditioning with respect to G j´1 , ζ j`1 , N T " n in the inside sum and then by applying Lemma 5, for j ą k and Lemma 6. More specifically, for any n ě 1, k P N n and k ă j ď n, We then apply Lemma 13 to the above identity. In order to do it, we make use of the following decomposition where G 1 and G 2 are the measurable functions defined for j ě k`1 by ı .

Now observe that the weights Ð Ý θ
B˚e j´1,j`1 , Ð Ý θ Bfe k´1,k`1 satisfy the time degeneracy estimates in the sense of (28) and so do the weights Ð Ý θ e i ,θ i , so that from the tower property of conditional expectation and Lemma 8 with p " 1, the following estimates hold with the convention ζ 0 " 0, ζ n`1 " T on tN T " nu.
Step 6: Integrability properties. The proof of the L p pPq-integrability, p P r0, 2q, follows from the time degeneracy estimates for each weight in the Lemmas 4, 5, 6 and 8 combined with a similar argument to the one employed at the end of the proof in Section 10 (see the discussion following Lemma 8). This also implies that the infinite sum over n converges absolutely and therefore after a re-ordering of the different terms one obtains the claimed formula. This concludes the proof. Remark 6. (i) The right hand side of the IBP formula may alternatively be written in a longer but maybe more appealing format as (ii) The restriction f pLq " 0, can be easily removed if one considers the test functionf pxq " f pxq´f pLq instead of f .
Note that, the previous theorem not only yields an IBP formula which is suitable for Monte Carlo simulation, it also provides a probabilistic representation for the derivative, with respect to the terminal point, of the transition density of the killed process at time T . To be more specific, under assumption (H), for any bounded measurable f defined on rL, 8q, one deduces Erf pX T q1 tτ ąT u s " ş 8 L f pzqppT, x, zq dz where p0, 8qˆrL, 8q 2 Q pT, x, zq Þ Ñ ppT, x, zq is the transition density of the killed process at time T starting from x at time 0. Moreover, from Theorem 5, by a standard approximation argument that we omit, one deduces that z Þ Ñ ppT, x, zq is differentiable on rL, 8q. 5 The next result provides a probabilistic representation of this derivative from which directly stems an unbiased Monte Carlo simulation method.

Corollary 1.
Under assumption (H), the transition density of the killed process at time T is differentiable with respect to its terminal point. Moreover, for all pT, x, zq P p0, 8qˆrL, 8q 2 the following probabilistic representation holds Here p s denotes the transition density of the Markov chainX s , that is, z Þ Ñ p s pT´ζ NT , x, zq is the density of the r.v.X s NT`1 conditional on X s NT " x ( . In particular, for s " B NT`1 NT`1 or B NT`1 NT`1 , the dynamics (26) readily gives p s pT´ζ NT ,X s NT , zq " gpapLqpζ NT`1´ζNT q, z´pLp1´µpX NT qq`X NT µpX NT qq sincē X s NT "X NT .

Bismut-Elworthy-Li type formula
In this section, we briefly derive the IBP formula for B x Erf pX T q1 tτ ěT u s commonly referred in the literature as the Bismut-Elworthy-Li formula. As we will see, obtaining this formula is simpler than the IBP formula derived in Theorem 5 as it does not involve boundary weights. Therefore the proof of this formula will not require any merging procedure contrary to the one of Theorem 5.
First, we give the transfer of derivatives lemma which is carried out forward in time in comparison with Lemma 4 where it is done backward in time. The proof is similar to the one of Lemma 4 given in Section 10.2 and therefore we omit it. Note that due to the change of time direction some changes of notation and sign occur among other changes in the formulae for the weights.
Theorem 7. Let f P C 1 b pRq such that f pLq " 0. Using the weights defined above, the following Bismut-Elworthy-Li formula is satisfied for any initial point x P rL, 8q: Moreover, the r.v. appearing inside the expectation in the right-hand side of the above equality belongs to L p pPq, for any p P r0, 2q.
In most of the arguments below, we will work on the set tN T " nu. In order to perform a forward induction argument through the Markov chain structure, we define for k P N n`1 the functions We let p F n`1 pX n`1 q :" f pX n`1 q and the following recursive relation is satisfied for k P N n Then, iterating the transfer of derivative formula in Lemma 7 for k P N n , we obtain 6 To further simplify the first term on the right-hand side of the above equation, we use the tower property of conditional expectation, the integration by parts formula (11) and the extraction formula (12) to obtain ErD k p F k pX k q1 D k,n Ý Ñ θ e k | G k´1 , T n`1 s " Er p F k pX k q1 D k,n I k p Ý Ñ θ e k q | G k´1 , T n`1 ś Er p F k pX k qδ L pX k q Ý Ñ θ e k | G k´1 , T n`1 s. In the case k " n, using the transfer of derivative formula of Lemma 7 on the last time interval and then performing the IBP formula (11), noting that f pLq " 0, we obtain the representation B x Erf pX n`1 q n`1 ź i"1 1 Di,nθi | T n`1 s " ErD n`1 f pX n`1 q n`1 ź where we remind the reader that we previously set Ý Ñ θ e n`1 " 0 in Lemma 7 for notational convenience. At this stage we emphasize that all the boundary terms in (40) vanish. In fact, as f pLq " 0, the boundary term in the last interval vanishes as stated in Lemma 7. For the other intervals, from (38), we claim that for j P N n , j | G j´1 , T n`1 s " 0. To see this, we note that forX j " L, one hasX j`1 " L`σpLqZ j`1 , which is independent of ρ j`1 , andθ j`1 given by (17) reduces toθ j`1 " 2p2ρ j`1´1 qλ´1p 1 2 papX j`1 q´apLqqI 2 j`1 p1q`pbpLq´a 1 pX j`1 qqI j`1 p1qà 2 pXj`1q 2´b 1 pX j`1 qq. Therefore the conclusion follows by conditioning with respect to X j " L ( in (41) and by noting that Er2ρ j`1´1 |G j´1 , T n`1 ,X j " Ls " 0.
From this property the identity (39) becomes Now, for each k P N n , one can multiply the above equality by the length of the interval on which the local IBP formula is performed, namely ζ k´ζk´1 and sum them over all k. For the last interval, we multiply (40) by T´ζ n . This gives Erf pX n`1 q n`1 ź i"j`1 From the above formula, the L p -moment estimate, p P r0, 2q, follows by similar arguments as described at the end of the proof of Theorem 2. Finally, one concludes by using the Lebesgue differentiation which yields T B x Erf pX n`1 q n`1 ź i"1 θ i |T n`1 s 1 tNT "nu ı and summing the previous formula over n.
Remark 8. (i) The right hand side of the IBP formula may alternatively be written as (ii) We note that the above formula does not involve any merging procedure. This is due to the fact that only the first derivative is being considered here. In fact, the key property (41) would not be satisfied if one considers second order derivatives. Therefore, a merging procedure similar to the one described in Section 4.3 would be necessary.
(iii) Similarly to the previous section, the above theorem yields a probabilistic representation for the derivative of the transition density of the killed process at time T from which stems an unbiased Monte Carlo simulation method. In contrast with Corollary 1, the derivative is taken with respect to the starting point. This can be seen by formally taking the Dirac mass at point z as a test function in Theorem 7. 7

Corollary 2.
Under assumption (H), the transition density of the killed process at time T is differentiable with respect to its starting point. Moreover, for all pT, x, zq P p0, 8qˆrL, 8q 2 the following probabilistic representation holds / -fi ffi fl .

Achieving finite variance by importance sampling
The previous probabilistic representations of Theorems 2, 5 and 7 as well as Corollaries 1 and 2 allow to devise an unbiased Monte Carlo simulation. However, in general, its use is hampered by the fact that the variance is infinite as suggested for instance by the moment estimate of Lemma 9. The main tool that we develop here in order to circumvent this issue consists in employing an importance sampling scheme on the jump times of the Poisson process N as originally proposed by Andersson and Kohatsu-Higa [3]. Since the arguments developed below follow similar lines of reasonings as those employed in [3], we will omit some technical details.
Let us first introduce a renewal process in the following sense: Definition 5. Let pT n q ně1 be a sequence of random variables such that pT n´Tn´1 q ně1 , with the convention T 0 " 0, are i.i.d. with density f and c.d.f: t Þ Ñ F ptq " ş t 8 f psq ds. Then, the renewal process J :" pJ t q tě0 with jump times pT n q ně1 is defined by J t :" ř ně1 1 tTnďtu . As previously done, we assume that J is independent of the Brownian motion W . It is readily seen that tJ t " nu " tT n ď t ă T n`1 u and by an induction argument that we omit, one may prove that the joint distribution of pT 1 ,¨¨¨, T n q is given by for any map Φ : ∆ n ptq Ñ R satisfying Er1 tJt"nu |ΦpT 1 ,¨¨¨, T n q|s ă 8. Usual choices that we will consider are the followings: Examples: (1) If the density function f is given by f ptq " λe´λ t 1 r0,8q ptq so that F ptq " 1´e´λ t , t ě 0, for some positive parameter λ, then J is a Poisson process with intensity λ. (2) If the density function f is given by f ptq " 1´ᾱ τ 1´α 1 t α 1 r0,τ s ptq, so that F ptq " pt{τ q 1´α , t P r0,τ s, for some parameters pα,τ q P p0, 1qˆp0, 8q, then J is a renewal process with r0,τs-valued Betap1ά , 1q jump times. (3) More generally, if the density function f is given by f ptq "τ 1´α´β Bpα,βq 1 t 1´α pτ´tq 1´β 1 r0,τ s ptq, so that F ptq " Bpt{τ , α, βq{Bpα, βq, r0, 1s Q x Þ Ñ Bpx, α, βq being the incomplete Beta function, for some parameters pα, β,τ q P p0, 1q 2ˆp 0, 8q, then J is a renewal process with r0,τ s-valued Betapα, βq jump times.
Having these definitions at hand, as done before, we define the partition π of r0, T s given by π :" t0 ": τ 0 ă¨¨¨ă τ JT ď T u with τ i :" T i^T . The probabilistic representation of Theorem 2 then becomes where the dynamics of the Markov chainX is given by (6) with innovations Z i`1 " W τi`1´Wτi , i " 0,¨¨¨, J T ( , with the set D i,n :" tX i ě L, J T " nu , i PN n`1 , and for i PN JT Moreover, the following estimates holds on the set tJ T " nu (47) ti"n`1u 1 Dn,n p1´F pT´τ n qq p E n,n where C is positive constant independent of n. Hence, using (47) and then (44), for all p ě 1, one gets As already mentioned before, the previous estimate is actually quite sharp and the series appearing in the right-hand side is finite for p P r0, 2q if J is a Poisson process. In order to achieve a finite variance, one has to select the law of the jump times suitably. Indeed, if for instance J is a renewal process with r0,τ s-valued Betap1´α, 1q jump times, α P p0, 1q,τ ą T , a simple computation shows that the above series is finite as soon as´p 2`α pp´1q ą´1, that is, pp 1 2´α q ă 1´α. In particular, taking α " 1{2, it is readily seen that the moment of order p of the random variable appearing inside the expectation in the right-hand side of (45) is finite for all p ě 1. Similarly, if J is a renewal process with r0,τs-valued Betap1´α, 1´βq jump times, pα, βq P p0, 1q 2 ,τ ą T , the integral appearing in the right-hand side of the above inequality is finite as soon as pp 1 2´α q ă 1´α and´pβ ă 1´β, the later condition being always satisfied. In particular, taking α " 1{2 and any β P p0, 1q, any moment of order p for p ě 1 is finite.
We now provide the probabilistic representation for the two IBP formulas using the above importance sampling technique. We thus redefine the weights appearing in the first IBP formula keeping in mind that the sequence pθ i q 1ďiďJT`1 is now given by (46). The new random variables p Ð Ý θ and Ð Ý θ B JT`1 :" 0, Ð Ý θ e JT`1 :" 2p1´F pT´τ JT qq´1 and Ð Ý θ c JT`1 :"´2p1´F pT´τ JT qq´1pσ 1 σq JT pTτ JT qI 2 JT`1 p1q. For the boundary merging weights of Lemmas 5 and 6, we first redefine the Markov chain X B with dynamics (26) by modifying the corresponding increments Z j,i`1 as done previously. We then set D B j,i`1,n :" tX B j,i`1 ě L, J T " nu, j " i´1, i. Finally, on the set tJ T " nu, pX j´L qd n`1 for j " n´1, n withd n`1 :" pΦg´1qpapLqpτ n`1´τj q, Z j,n`1 q. At this stage, it is important to remark that the above new weights satisfy a time degeneracy estimate, with a slight modification of Definition 3. The notation E i,n rXs now is used for the expectation of X conditional on G i , T n`1 , ρ n`1 , J T " n ( and one considers the corresponding norm }.} p,i,n . We now say that a weight H P S i,n satisfies the time degeneracy estimate if for all p ě 1 ,n ď Cf pτ i´τi´1 q´1pτ i´τi´1 q´1 2 in the case that i P N n and 1 Dn,n › › 1 Dn`1,n H › › p,n,n ď C in the case that i " n`1. In a completely analogous manner as done in Lemma 4, the weights Ð Ý θ a i for a P te, c, Bu satisfies the time degeneracy estimate (51). For the boundary merging weights, we replace (51) for H P S j,i`1,n pX B q by With the above new definitions and properties, we finally redefine the corresponding weights Ð Ý θ s with the related Markov chainX s for each s P S n`1 or s P 9 S k n`1 on the time partition π :" t0 ": τ 0 ă¨¨¨ă τ n`1 :" T u of the underlying renewal process on the set tJ T " nu. This is done in a completely analogous manner as presented in the subsection 5.2.
We can now restate Theorem 5 as follows. For any function f P C 1 b pRq satisfying f pLq " 0, Note also that the corresponding other formulation of Remark 6 also holds. Moreover, if J is a renewal process with r0,τs-valued Betap1´α, 1q jump times,τ ą T , with α satisfying pp 1 2´α q ă 1´α or if J is a renewal process with r0,τ s-valued Betap1´α, 1´βq jump times,τ ą T , with α and β such that pp 1 2´α q ă 1´α and β P p0, 1q, then the r.v. appearing inside the expectation of the right-hand side of the above equality belongs to L p pPq for any p ě 1. The proof follows similar lines of reasonings as those employed above in order to deal with the new probabilistic representation (45) and is thus omitted.
We proceed similarly for the BEL formula of Theorem 7. Namely, we redefine the weights Ý Ñ θ a , for a P te, c, Bu as follows Ý Ñ θ e i :"2pf pτ i´τi´1 qq´1` and also set Ý Ñ θ e n`1 :" 2p1´F pT´ζ n qq´1p1`p2ρ n`1´1 qσ 1 n Z n`1 q, Ý Ñ θ c n`1 :" 0. These new weights satisfy the time degeneracy estimate (51). Then, the following BEL formula is satisfied for any f P C 1 b pRq satisfying f pLq " 0 and any initial point x P rL, 8q: Moreover, the r.v. appearing inside the expectation in the right-hand side of the above equality belongs to L p pPq, for any p ě 1 in the case of r0,τs-valued Betap1´α, 1q jump times or r0,τ s-valued Betap1´α, 1´βq jump times under the condition pp 1 2´α q ă 1´α and β P p0, 1q. In particular, choosing α " 1{2, the L p pPq moment is finite for any p ě 1.

Numerical tests
In this section, we provide some numerical results for the unbiased Monte Carlo simulation method based on the probabilistic representation formula established in Theorem 2 for the marginal law of the killed process and the Bismut-Elworthy-Li (BEL for short) formula of Theorem 7. A similar numerical analysis could be done for the IBP formula established in Theorem 5 but we restrict to the two aforementioned case for sake of simplicity. For a one dimensional Brownian motion W , we thus consider the following one-dimensional SDE with dynamics σpX s qdW s , x P R and we choose the coefficients b, σ and the test function f as follows σpxq "σˆpsinpωxq`2q, bpxq "´x for some positive constantσ ą 0. With this particular choice, we first observe that assumption (H) is clearly satisfied and that a direct computation yields Lf pxq :" bpxqf 1 pxq`1 2 σ 2 pxqf 2 pxq " 0. The process pf pX t qq tě0 is thus a martingale and by Doob's stopping theorem one gets Erf pX τ^T qs " f px 0 q where τ " inf tt ě 0 : X t ď Lu. We now shift the function f by considering hpxq " f pxq´f pLq satisfying hpLq " 0 instead of f . It is readily seen that h satisfies Lhpxq " 0 so that ErhpX τ^T qs " ErhpX T q1 tτ ąT u s " hpx 0 q. Note that since h is not bounded but of polynomial growth, only Theorem 2 directly applies. However, the extension of Theorem 7 (and also of Theorem 5) to polynomially growing function can be performed by a standard approximation argument noting that all moments ofX NT`1 and X T are bounded. Having this extension in mind, by Theorem 7, one has We select the following parameters: T " 0.5, L " 0, c 0 " 0, c 1 " 1, c 3 " 1 and x 0 " 1 so that hpx 0 q " 2 and T h 1 px 0 q " 2. We use three different parameters sets forσ " ω " 0.1, 0.2, 0.3. We examine the performance of the proposed Monte Carlo estimator with respect to the previous sets of parameters when one uses the Exponential sampling (the distribution of the jump times is exponential with parameter λ) as it is written in Theorem 2 and Theorem 7 and when one uses an importance sampling technique with Beta distribution with parameters pγ,τ q for the jump times of the renewal process, as exposed in Section 7, which allows to achieve finite variance for our estimators. We note that though the variance is not finite in the case of Exponential time sampling, we include it here in order to compare its performance with the Beta sampling scheme. In both cases, we first select the optimal parameters which minimize the variance of the estimator using few samples, that is, the optimal λ in the case of Exponential sampling and the optimal pγ,τ q in the case of Beta sampling. We then use M " 4ˆ10 6 i.i.d. samples to estimate the considered quantities. The results are summarized in the two tables below. The first column of Table 1 and Table 2 provides the value of the two parametersσ " ω. The second column (resp. third column) of Table 1 provides the estimated value of the quantity ErhpX T q1 tτ ąT u s with its associated variance, L 1 pPq-error and 95%-confidence interval in the case of Exponential sampling (resp. Beta sampling). The second column (resp. third column) of Table 2 provides the estimated value of the quantity B x ErhpX T q1 tτ ąT u s with its associated variance, L 1 pPq-error and 95%-confidence interval in the case of Exponential sampling (resp. Beta sampling).
Most notably, we observe that in both tables the performance of our estimators quickly deteriorates asσ " ω increases. Actually, for large values ofσ, ω, say greater than 0.4, the variance becomes difficult to estimate from the simulations and the obtained estimates become unreliable. We also see that the  Table 2. Unbiased Monte Carlo estimation for the quantity B x ErhpX T q1 tτ ąT u s based on Theorem 7 by Exponential and Beta sampling with its associated variance and 95%confidence interval.
Beta time sampling method outperforms the Exponential time sampling for all values of the considered parameters especially for large values ofσ. This behavior was already observed in [3] and is reminiscent of unbiased simulation methods for multidimensional diffusion processes. It was thus expected here since there is no hope that our estimator will overcome this problem. To circumvent this issue, one may resort to more sophisticated method such as the second order approximation method developed by [4].

Some Conclusions
In the present work, we presented a probabilistic representation formula for the marginal law of a killed process based on a basic Markov chain which is obtained using the reflection principle. From this representation, we established two IBP formulae, one being of BEL's type, from which directly stem an unbiased Monte Carlo method. The main element used in this construction is a suitable tailor-made Malliavin calculus for the underlying Markov chain. For this reason, we do not need to use the full fledge power of Malliavin calculus by closing the derivative operator but just the concepts for simple discrete time Markov chain.
The methodology developed here seems to follow a general pattern that could be used to obtain IBP formulae for some other irregular functionals of the Wiener process for which boundary problems may appear such as the exit time, the local time, the running maximum or the occupation time of a multidimensional diffusion process. Although the problem investigated here focuses in the one dimensional case, we believe that the approach developed here also extends to some multi-dimensional cases for which the reflection principle is well understood and densities for basic approximation processes are known, see e.g. [2], [10] and the references therein.
On the other hand, it seems difficult at this moment to generalize the methods in [23] and [21] to the multi-dimensional case or even to obtain an amenable integration by parts formula based on a Markov chain using these formulations or a Lamperti like transform. We will discuss this extension in future works as well as their implementation for simulation purposes.  (1). Let P denote the semigroup operator associated with the killed process. That is, for a measurable and bounded function f , one defines P t f pxq " E " f pX t q1 tτ ątu ‰ . We remark that due to the indicator function, this semigroup is not conservative.
The heuristic argument in order to obtain the probabilistic representation is to use Itô's formula on an approximation process obtained from (1) by removing the drift and freezing the diffusion coefficient at the starting point. From Itô's formula, one obtains a one step expansion of the law of X around the law of the Markov chainX. Then, an IBP formula based on the Markov chainX has to be used to obtain the probabilistic representation of this first step expansion using one jump of the Poisson process N . Then, one just needs to iterate the first step expansion in order to obtain the full probabilistic representation.
In order to do this rigorously, one needs to use the regularity properties of the semigroup P which can be found in [15], Chapter VI and/or [16]. In particular, under assumption pHq on the coefficients, one obtains that if f is a smooth function such that lim xÓL f pxq " f pLq " 0 then P f P C 1,2 pp0, T sˆrL, 8qq and it satisfies B t P t f " LP t f on rL, 8q, t ą 0, where L " 1 2 aB 2 x`b B x is the infinitesimal generator of X. Moreover, one has sup 0ďtďT |B ℓ x P t f | 8 ď C, for ℓ " 1, 2, for some positive constant C :" CpT, a, bq. Under the condition that the test function f vanishes at L, we obtain that P satisfies its corresponding Dirichlet boundary condition P t f pLq " f pLq " 0 together with P 0 f pxq " f pxq for all x ě L.
The argument used to obtain the probabilistic representation starts by applying Itô's formula to pP T´t f pȲ t^τ qq tPr0,T s where the processȲ is defined in Lemma 1 on the interval r0, T s andτ is its associated exit time. We also recall the definition of the approximation process obtained from the reflection principle of Lemma 1, namelyX s,x t " ρx`p1´ρqp2L´xq`σpxqpW t´Ws q (with the shorten notation X t "X 0,x t ): u P u f pȲ s qˇˇu "T´s`1 2 apxqB 2 x P T´s f pȲ s q˙1 tτąsu ds E " P T f pxq`ż T 0ˆ1 2`a pxq´apȲ s q˘B 2 x P T´s f pȲ s q´bpȲ s qB x P T´s f pȲ s q˙1 tτ ąsu ds E " P T f pxq`2p2ρ´1q ż T 0ˆ1 2`a pxq´apX s q˘B 2 x P T´s f pX s q´bpX s qB x P T´s f pX s q˙1 tXsěLu ds.
We now rewrite the previous representation using the Markov chain pX i q 0ďiďNT`1 defined by (6) together with the Poisson process N . From the previous identity, we get P T f pxq E " f pX T q2p2ρ´1q 2p2ρ´1q ż T 0ˆ1 2`a pX s q´apxq˘B 2 x P T´s f pX s q`bpX s qB x P T´s f pX s q˙1 tXsěLu ds (54) E " f pX NT`1 qθ NT`1 1 tNT "0ù e λT 2λ´1p2ρ NT´1 q " 1 2 papX 1 q´apxqqB 2 x P T´ζ1 f pX 1 q`bpX 1 qB x P T´ζ1 f pX 1 q * 1 tX1ěLu 1 tNT "1u Next, we apply the IBP formula (11) with respect to the r.v.X 1 in the above expression. The formula is applied once to the terms associated with the drift coefficient b and two times with respect to the terms related to the diffusion coefficient a. In order to do that one first has to take the conditional expectation E 0,1 r.s in the second term of the above equality. As these IBPs involve the indicator function 1 tX1ěLu , the correct calculation has to be done using the theory in [22], Chapter V.9. 8 This yields P T f pxq E "f pX NT`1 qθ NT`1 1 tNT "0u`e λT 2λ´1p2ρ NT´1 qˆ1 2 I 1``a pX 1 q´apxq˘1 tX1ěLu˘Bx P T´ζ1 f pX 1 q`I 1`b pX 1 q1 tX1ěLu˘PT´ζ1 f pX 1 q˙1 tNT "1u E "f pX NT`1 qθ NT`1 1 tNT "0u`e λT 2λ´1p2ρ NT´1 qˆ1 2 I 1``a pX 1 q´apxq˘˘B x P T´ζ1 f pX 1 q`I 1`b pX 1 q˘P T´ζ1 f pX 1 q˙1 tX1ěLu 1 tNT "1u where we used the extraction formula (12) applied to the r.v. 1 tXsěLu and Lemma 11 (taking first the conditional expectation w.r.t ζ 1 ) for ℓ " 1, k " 0 for the last equality.
add them together. Therefore proving the following general statement will suffice: For any j " 0, 1, 2 and real-valued f P C 0 p pR 2 q, E " f pX i ,X i`1 q1 Di`1,n I j i`1 p1qδ L pX i qp2ρ i´1 qI i p1qˇˇG i´1 , ζ i`1 , N T " n ‰ (65) Step 2: Proof of (65). Case j " 0. In order to foster the understanding of the proof, let us first consider the case j " 0. The following arguments also gives rise to the definition of the merged boundary process.
First, observe that as p2ρ i´1 qσ i´1 Z i " L´X i´1 :" Y i´1 on the set X i " L ( , one has p2ρ i´1 qδ L pX i qI i p1q " δ L pX i q Y i´1 a i´1 pζ i´ζi´1 q . Next, we rewrite the right hand side of (65), namely, using the independence of W and N , we have We now use the following key property: for a real-valued bounded G i´1ˆB pR 2 q-measurable function h : ΩˆR 2 Ñ R, one has Erhpζ i , ζ i`1 q | G i´1 , ζ i`1 , N T " ns " Erhps`U, tq| G i´1 , ζ i`1 , N T " nsˇˇs "ζi´1, t"ζi`1 .
where U " Up0, t´sq is independent of G i´1 and N . More precisely, by using the tower property for conditional expectation and (66), we obtain E " f pX i ,X i`1 q1 Di`1,n δ L pX i qp2ρ i´1 qI i p1qˇˇG i´1 , ζ i`1 , ρ i , N T " n ‰ Y i´1 a i´1 s gpa i´1 s, Y i´1 qgpapLqppζ i`1´ζi´1 q´sq, z´Lq ds.
In order to compute the above time convolution and show that the above conditional expectation can be rewritten using the boundary processX B , we apply Lemma 10, in particular (??), noticing that Yi´1 ai´1s gpa i´1 s, Y i´1 q " B 2 gpa i´1 s, |Y i´1 |q, to obtain E " f pX i ,X i`1 q1 Di`1,n δ L pX i qp2ρ i´1 qI i p1qˇˇG i´1 , ζ i`1 , ρ i , N T " n ‰ " 1 a i´1 pζ i`1´ζi´1 q ż 8 L f pL, zqgpapLqpζ i`1´ζi´1 q, z´pLp1´µpX i´1 qq`X i´1 µpX i´1 qqq dz, This finishes the proof for the case j " 0.
Step 3: The general case: j " 1, 2. Using equality (62), one can similarly prove (65) in the case j " 1, 2. In fact, a useful formula to prove this general case directly from (62) is the following property for i P N n which was stated in (13).
Step 4: The general case follows by linearity. Note that in the general case of (29), we have that Ð Ý θ e i`1 is given by (22). This expression can be rewritten using the extraction formula (12) as linear combinations of terms of the type f j pX i ,X i`1 qI j i`1 p1q for some particular functions f j P C pR 2 q. Therefore the result in (65) applies.
One uses again (12) forĪ in order to rewrite the resulting expressions in a compact form. For this, note that we have conveniently definedD as the adjoint ofĪ for which the extraction formulae are satisfied.