Speeding up the Euler scheme for killed diffusions

Let $X$ be a linear diffusion taking values in $(\ell,r)$ and consider the standard Euler scheme to compute an approximation to $\mathbb{E}[g(X_T)\mathbf{1}_{[T<\zeta]}]$ for a given function $g$ and a deterministic $T$, where $\zeta=\inf\{t\geq 0: X_t \notin (\ell,r)\}$. It is well-known since \cite{GobetKilled} that the presence of killing introduces a loss of accuracy and reduces the weak convergence rate to $1/\sqrt{N}$ with $N$ being the number of discretisatons. We introduce a drift-implicit Euler method to bring the convergence rate back to $1/N$, i.e. the optimal rate in the absence of killing, using the theory of recurrent transformations developed in \cite{rectr}. Although the current setup assumes a one-dimensional setting, multidimensional extension is within reach as soon as a systematic treatment of recurrent transformations is available in higher dimensions.


Introduction
Let X be a diffusion on some filtered probability space taking values in ( , r) and solving where B is a Brownian motion, and ζ := inf{t ≥ 0 : X t / ∈ ( , r)} is the first exit time from the interval ( , r).The process is killed at ζ and sent to a cemetery state.
Let's assume that at least one of the boundaries are accessible and ζ is finite a.s. and consider E[g(X T )1 [T <ζ] ] for a given function g and a deterministic T .Letting g take the value 0 at the cemetery state, one can rewrite this expression in terms of the killed diffusion as E[g(X T )1 [T <ζ] ].Such computations appear very naturally in many applied problems of science, engineering, and finance.For instance, in Mathematical Finance theory, such an expectation corresponds to the price of a barrier option with payoff g and maturity T written on a stock whose price process is given by X.The barrier feature renders the option worthless if the stock price hits one of the accessible boundaries before the maturity of the option.
A closed form expression for E[g(X T )1 [T <ζ] ] is rarely available even in this one-dimensional setting.Thus, one needs to resort to an approximation scheme for an answer.Arguably the easiest approach is to run a standard Euler-Maruyama scheme on the SDE (1.1) by setting where X0 = x, t 0 = 0, N > 0 is an integer, t n = nT N for n = 1, . . .N , and compute E[g( XT )1 [T <τ ] ], where τ is the first time that the discrete-time process ( Xtn ) N n=0 hits any of the barriers.Under standard regularity conditions on the diffusion process and g, such a scheme indeed converges as N → ∞.However, it converges at a rate much slower than a standard Euler-Maruyama scheme applied to a diffusion process that is not killed at accessible boundaries.Indeed it was shown by Gobet [16] that under standard hypothesis the above scheme for the killed diffusion converges weakly at rate N −1/2 as opposed to N −1 , which is the rate of weak convergence for the Euler-Maruyama scheme in the absence of killing (see, e.g., Talay and Tubaro [37] or Mikulevičius and Platen [29]).This rate is optimal since it is reached when X is a Brownian motion and g is an indicator function of a set strictly contained in ( , r) (see Siegmund and Yuh [36]).
C ¸etin [6] conjectured that using a recurrent transformation would bring the convergence rate back to N −1 .A recurrent transformation at heart is a change of measure that keeps the Markovian structure intact while transforming the process into a recurrent one.In particular X never touches the boundaries of ( , r) under the new measure Q. [6] shows that Q is locally absolutely continuous with respect to the original measure P, and X follows for some function h and a Q-Brownian motion W .That the above claim was a conjecture and not following immediately from the standard results on Euler-Maruyama schemes is that h h is explosive near boundaries and is not Lipschitz, which is in fact needed for X not to touch the previously accessible boundaries after the measure change.This can create significant difficulties with approximation and may even lead to divergence (see, e.g., the potential issues that may arise with non-Lipschitz drivers and methods on how to resolve them in [22] and [23]).
In this paper we prove this conjecture with a slight "twist."Note that if one applies the Euler-Maruyama scheme naively to (1.2), one obtains, as usual, a Brownian motion with drift whose parameters change at times of discretisation.This process will hit finite boundaries with positive probability, and therefore will exit the state space of X with positive probability.One way to overcome this is to impose an ad hoc reflection on the boundaries.However, this will introduce a local time term in computations requiring additional estimates on its convergence rate to 0. Moreover, it is far from obvious that reflection is the optimal resolution of problems arising from the discretised process exiting the domain.
We instead study a drift-implicit method that keeps the state space intact after discretisation.To see this, suppose that b ≡ 0, which can be obtained by changing the scale if necessary, and consider the backward Euler-Maruyama scheme where h becomes a concave function.Note that different than what one would expect from a backward scheme (see, e.g.Mao and Szpruch [28], Alfonsi [2], Alfonsi [3], and Neunkirch and Szpruch [30] to name a few) the σ 2 -term in the drift of (1.2) is still evaluated at X tn .This stems from the fact that (1.2) with b ≡ 0 should be viewed as a time-changed version of where the time change is given by t 0 σ 2 (Y s )ds.We make an extensive use of this correspondence in our proofs.
Our main result is Theorem 4.1 which proves that the rate of weak convergence of the above backward Euler-Maruyama scheme is N −1 under standard assumptions on the diffusion process.Moreover, there is no single h function that achieves this rate.We show that any nonnegative concave h vanishing at accessible boundaries can be used to obtain this convergence rate as long as it satisfies some mild growth conditions.Such functions are easy to construct and we study in Section 5 the construction of some particular h-functions to compute approximate prices for barrier options in a Black-Scholes framework.Our numerical results are very promising and error terms very rapidly converge to 0 even with a small number of iterations.Moreover, in the case of a particular local volatility model with double barriers, our method yields smaller error terms than the so-called Brownian bridge method when the number of discretisations is reasonably large.
We are not the first to consider implicit schemes for studying diffusions with infinite lifetime and taking values in a strict subset of R. Alfonsi in [2,3] and Neunkirch and Szpruch [30] consider such scalar processes whose SDE representation is given by and f satisfying the conditions of a Feller test ensuring that Y takes values in ( , r) (see also [11] in the special case of Cox-Ingersoll-Ross (CIR) process).[3] and [30] show that a the drift implicit Euler scheme for Y converge strongly with rate N −1 if f satisfies certain integrability conditions including However, this condition cannot be satisfied by h that paves the way for recurrent transformation rendering X recurrent and following (1.2).Indeed, if the dynamics of X are given by (1.2) where h is a function satisfying the condition of Theorem 3.2 in [6], b ≡ 0, and σ ≡ 1, then where C is an adapted, continuous and increasing process.Note that h is concave and Moreover, h never vanishes at the boundary points where h does.Thus, (1.5) implies that the local martingale in the above Q-Doob-Meyer decomposition of 1 h(X) is a true martingale, which in turn will yield 1  h exp(−A) is a true martingale, where h(Xt) .But this would imply that P ∼ Q (when restricted to F t , with (F t ) representing the underlying filtration) in view of the absolute continuity relationship manifested in Theorem 3.2 in [6].This is not possible since for an arbitrary t > 0 Q(ζ < t) = 0 while P(ζ < t) > 0.
The estimates obtained by the authors in [3] and [30] rely on the Burkholder-Davis-Gundy (BDG) inequality which requires the corresponding local martingale be a true martingale.As 1 h(X) is a strict local submartingale under Q, one needs to develop new techniques to arrive at the needed estimate for convergence theorem.
This brings to the fore another novelty of our paper.Given the impossibility of the use of BDG inequality we use potential theoretic methods that yield the boundedness of inverse moments of h(X) under Q, which is crucial for obtaining the weak convergence result in our paper (or a strong convergence type results considered by Alfonsi, Neunkirch and Szpruch).We use the theory of Kato class potentials to show the boundedness of required moments.Kato potentials are one of the fundamental objects in the study of Schrödinger operators (see, e.g.[1], [10], [8], [9]).We show in Theorem 2.1 that the additive functional dA t = − 1 2 h (Xt) h(Xt) belongs to a particular Kato class defined in [8], which in turn yields the boundedness of the inverse moment of 1 h ( X tn ) (uniformly in N ) in conjunction with a comparison argument via Lemma 3.2.The potential theory also helps us to prove uniform bounds on the moments of integral functionals of h −2−p ( X t ) (see Theorem 3.1 for an exact description).
Our methodology offers hope to study the convergence rates for CIR processes that do not satisfy (1.5).We also show in this paper that if one considers the 3-dimensional Bessel process, the implicit scheme in (1.3) converges weakly at rate N −1 .Clearly, (1.5) is violated since the reciprocal of a 3-dimensional Bessel process is a prime example of a strict local martingale.This process satisfies the conditions of Theorem 4.1 and one obtains the optimal convergence rate for sufficiently smooth g.We leave the investigation of the convergence rate for conservative diffusions on (0, ∞) satisfying (1.4) in the absence of condition (1.5) to a future study.
Although our analysis assumes a one-dimensional framework, a close look into our technical analysis reveals that our convergence result does not depend heavily on this assumption apart from the comparison argument used in Lemma 3.2.In particular it is relatively clear how to obtain a version of Theorem 2.3 in the multidimensional case using well-known potential theoretic arguments on Kato classes.However, our main obstacle in not being able to extend our results to a multidimensional setting is the absence of a systematic study of recurrent transformations in higher dimensions.Also note that Lemma 3.2 is only used to obtain estimates on h( X), which is always a one-dimensional object with X referring to the continuous Euler scheme.Such a study and its applications to Euler methods for killed diffusions will be the subject of future research.
The outline of the paper is as follows.Section 2 fixes the setting, gives a brief summary of results for recurrent transformations needed for this paper together with novel inverse moment estimates, and introduces the backward Euler-Maruyama scheme that is tailored for our purposes.Section 3 obtains the moment estimates that will be needed for the weak convergence analysis.Theoretical predictions are confirmed via numerical studies in Section 5 and Section 6 concludes.

Preliminaries
Let X be a regular diffusion on ( , r), where −∞ ≤ < r ≤ ∞.We assume that infinite boundaries are inaccessible and if any of the boundaries are reached in finite time, the process is killed and sent to the cemetery state ∆.This is the only instance when the process can be 'killed', we do not allow killing inside ( , r).The set of points that can be reached in finite time starting from the interior of ( , r) and entrance boundaries will be denoted by I.That is, I is the union of ( , r) with the regular, exit and entrance boundaries.The law induced on C(R + , I), the space of I-valued continuous functions on [0, ∞), by X with X 0 = x will be denoted by P x as usual, while ζ will correspond to its lifetime, i.e. ζ := inf{t > 0 : X t / ∈ ( , r)}.We also introduce the set I ∆ := I ∪ {∆} and extend any I-valued Borel measurable function f to I ∆ by setting f (∆) = 0 unless stated otherwise.The filtration (F 0 t ) t≥0 will correspond to the natural filtration of X, Ft will be the universal completion of F 0 t , and F t = Ft+ so that (F t ) t≥0 is a right continuous filtration.We will also set F := t≥0 F t .We refer the reader to [5] for a summary of results and references on one-dimensional diffusions.The definitive treatment of such diffusions is, of course, contained in [24].
Since we are only concerned with the diffusion process until it is killed, we can assume without any loss of generality that X is on natural scale.The extra regularity conditions imposed in the following assumption are standard in the theory of Euler discretisations for SDEs.Assumption 2.1.X is a regular one-dimensional diffusion on ( , r) such that where σ : ( , r) → (0, ∞) is continuously differentiable with a bounded derivative, B is a standard Brownian motion, and ζ = inf{t > 0 : X t ∈ {l, r}}.Moreover, σ( +) (resp.σ(r−)) exist and is finite if (resp.r ) is finite.
Note that the speed measure m associated with X is given by m(dx) = 2σ −2 (x)dx on the Borel subsets of ( , r).
Since we are interested in diffusions with killing, the following assumption is needed to ensure that we are not dealing with a vacuous problem.Assumption 2.2.P x (ζ < ∞) > 0 for each x ∈ ( , r).
Let I 0 be the set of points in I that can be reached from its interior in finite time.Note that under Assumptions 2.1 and 2.2 there are only two cases to consider: Case 1: Both and r are accessible, which in turn implies and r are finite and Case 2: Only one of and r is accessible, which can be assumed to be without any loss of generality.
In particular, As is always finite as a result of the above convention, the following will also be assumed for convenience: As a transient diffusion on ( , r), X has a finite potential density, u : ( , r) 2 → R + , with respect to its speed measure (see Paragraph 11 in Section II.1 of [5]).That is, for any nonnegative and measurable f vanishing at accessible boundaries The potential density is symmetric and is explicitly known in terms of the scale function and the speed measure of X.This leads to the following specification of the potential density: The following is a direct consequence of Theorem 3.2 in [6].The reader is referred to [6] for all unexplained terminology.Theorem 2.1.Suppose that Assumptions 2.1-2.3 are in force.Let f : (0, r) → (0, ∞) be a continuous function such that (0,r) f (y)m(dy) < ∞ and (0,r) yf (y)m(dy) < ∞, and define h(x) := ( ,r) u(x, y)f (y)m(dy).
Then, the following hold: (1) (h, M ) is a recurrent transform of X, where ds .
(2) There exists a probability measure Q h,x on F that is locally absolutely continuous with respect to P x such that and W is an Q h,x -Brownian motion.(3) If S is a stopping time such that Q h,x (S < ∞) = 1, then for any F ∈ F S the following identity holds: where E h,x is the expectation operator with respect to the probability measure Q h,x .In particular, Note that h constructed above is a concave function that is twice continuously differentiable and satisfies on ( , r) 3) The class of concave functions h such that h = U f where f is a continuous function satisfying the conditions of Theorem 2.1 will be denoted by H 0 .
We shall also consider the following h-transformation when r = ∞: Theorem 2.2.Suppose that Assumptions 2.1-2.3 are in force and r = ∞.Let h(x) := x.Then, the following hold: (1) There exists a probability measure Q h,x on F that is locally absolutely continuous with respect to P x such that and W is a Q h,x -Brownian motion.(2) If S is a stopping time such that Q h,x (S < ∞) = 1, then for any F ∈ F S the following identity holds: , where E h,x is the expectation operator with respect to the probability measure Q h,x .In particular, The above result is well-known and the reader is referred to Theorem 6.2 in [14] for a proof in a much more general setting.Note that the h-transform of Theorem 2.2 does not produce a recurrent diffusion.Indeed, Q h,x (lim t→∞ X t = ∞) = 1 since the corresponding scale function is finite at ∞.
For ease of later reference define the set H to be the union of H 0 and the set that contains only the identity function when r = ∞.If r is finite, set H = H 0 .Lemma 2.1.Let h ∈ H.
(1) For any given z > 0 consider the function H defined by H is strictly increasing and H((0, r)) = R.
(2) h is increasing if r = ∞.However, h is bounded.In particular, for h ∈ H 0 , we have Proof.
Next, suppose h ∈ H 0 .Then, the dominated convergence theorem implies that h(0) = 0 as well as h(r) = 0 if r < ∞ since the potential density vanishes at finite endpoints.Moreover, as h is strictly concave and never vanishes in the interior of the state space, h (0) > 0. Thus, This proves the desired range for H when r = ∞.Indeed, in this case h is increasing, which in turn yields x ≥ 1.
If r < ∞, similar considerations imply h (r) < 0, and therefore This completes proof of the first assertion.
which is nonnegative and finite by the assumption on f .In particular, Thus, This yields the desired boundedness and the boundary levels for the derivatives.(3) First suppose r < ∞.Since h does not vanish at the boundaries, u(y, y)/h(y) is bounded. Moreover, This proves the claim when r < ∞.Now suppose that r = ∞.Thus, u(y, y) = y and The first integral on the right hand side is finite since f is m-integrable and y/h(y) is bounded on [0, 1] as h (0) > 0. The second integral is also finite since h(∞) > 0 and yf (y)m(dy) < ∞ by assumption.
Let g : I 0 → R be a continuous function vanishing at accessible boundaries and Then for h ∈ H and a deterministic T > 0 we have In order to approximate the expectation on the right side of (2.7) we shall use a backward Euler-Maruyama (BEM) scheme: Let N > 1 be an integer and define t n := n N T for n = 0, . . ., N .Set X0 = X 0 and proceed inductively by setting,Then proceed inductively by setting Note that in view of Lemma 2.1 the mapping x → x − z h h (x) is one-to-one and onto for any given z > 0. Thus, the above scheme is well-defined since σ(x) > 0 for all x ∈ (0, r).
As we shall see in Section 3 the following type of diffusion processes on (0, r) will play a crucial role: where ζ(Y ) denotes the first hitting time of 0 or r.Note that c = 0 corresponds to the recurrent transform defined above.
Theorem 2.3.Suppose that Assumptions 2.1-2.3 are in force, h ∈ H, and Y is a process defined by (2.9) with Y 0 = X 0 .Assume further that c ≤ 0 if r = ∞, and c = 0 if h(x) = x for all x.Then the following statements are valid: (2) For any stopping time S that is bounded Q h,X 0 -a.s.there exists a constant K that does not depend on X 0 such that .
(1) First observe that a scale function and speed measure for Y can be chosen as where d ∈ (0, r).Since s y (0) = −∞, 0 is an inaccessible boundary for Y .By the same token, r is also an inaccessible boundary when s y (r) = ∞, which will be valid when r < ∞ or c ≤ 0. (2) Define Z by and note that Z is a nonnegative Q h,X 0 -local martingale by a straightforward application of Ito's formula.By Theorem 62.19 in [35] there exists a probability measure P such that where β is a P -Brownian motion, and whenever S is a stopping time that is finite Q h,X 0 -a.s., one has where x − denotes the negative part of x and we drop the dependency on Y for ζ to ease the exposition.Suppose that S < R, Q h,X 0 , a.s.where R is a deterministic constant, and note that Let W c,y denote the law of the process Ỹ starting at y, where d Ỹt = dβ t + cdt and gets killed at hitting 0 or r.Thus, where C is the positive continuous additive functional of Ỹ with , , and d m is the associated speed measure of Y .Since a scale function and a speed measure of Ỹ can be chosen as (2.11) by another application of (2.6) due to the bounds obtained via (2.10) and (2.11), and the assumption on the choice of c when r = ∞.As and the right side converges to 0 by the dominated convergence theorem due to (2.12), we establish that µ C ∈ K 1 ( Ỹ ) (see Definition 2.2 in Chen [8]) by Proposition 2.4 in [8].Therefore, by Proposition 2.3 in [9] we arrive at sup y∈(0,r) for some constants d 1 and d 2 .This proves the claim.(3) Since the semigroup is self-dual with respect to the speed measure, for any nonnegative measurable f we have where In particular, when f (y) = q(ε, y, y * ) for some ε > 0, where q is the transition density of Y with respect to its speed measure, we obtain where L y * is the diffusion local time with respect to the speed measure.Letting ε → 0 we arrive at h 2+p (Ys) ds is lower semi-continuous.Note that the finiteness of the integral on the right hand side follows from the fact that |h (y)| ≥ α for some α > 0 on D.
Observe that where where v 1 is the 1-potential density of Y .Since v 1 is jointly continuous (see Paragraphs 10-11 in Chapter II of [5]), the claimed semi-continuity follows.Since for some K, the claim follows from the arbitrariness of y * .

Moment estimates for the continuous BEM scheme
In this section we will obtain some moment estimates, including inverse ones, that will be necessary to establish the weak rate of convergence.We start with the following consequence of Ito's formula.
Proof.The decomposition (3.1) follows from Ito's formula and straightforward calculations regarding the derivatives of the inverse function.
To prove the second assertion first observe that where we drop the dependency on t n and z to ease the exposition.
Observe that and that the claim follows immediately if h(x) = x since the term in the parenthesis in (3.2) becomes nonnegative.Thus, it remains to show the assertion when h ∈ H 0 .First consider the case r = ∞, and let u := h h and note that lim x→∞ u (x) = 0 by Lemma 2.1.
where the second equality is an application if L'Hospital's rule.Thus, for some c < 0 where x * := inf{x : h (x) = 0} > 0 by Lemma 2.1.An alternative representation for µ is given by Thus, we will be done if is bounded from below on (0, x * 2 ).Indeed, as h is bounded away from 0 on this interval, the hypothesis on h implies When r < ∞, we have in particular that σ is bounded.Moreover, for some constant K, which renders r bounded.Observing that the remaining terms in (3.4) has the correct sign completes the proof.
The next result is a key comparison result that relates the inverse moments of the BEM scheme to those of the process (2.9) and thereby provide estimates that are valid uniformly in N .Lemma 3.2.Suppose that h satisfies the conditions of Lemma 3.1, σ is bounded, r = ∞, and consider the BEM scheme defined by (2.8) for h ∈ H. Then for any non-decreasing and measurable function φ that does not change sign, we have where Y is the process defined by (2.9) with c = c 1 , c 1 is as in Lemma 3.1, and A is a continuous time-change defined by A 0 = 0 and Proof.Consider the process Y defined by Y t = X A −1 t .Dambis, Dubins-Schwarz Theorem (cf.Theorem V.1.6 in [32]) yields where µ t ≥ c 1 and β is a standard Brownian motion adapted to the filtration (F A −1 t ) t≥0 .Then the comparison theorem for stochastic differential equations (cf.Theorem 2.10 in [7]) show that where This completes the proof.The main moment estimates are collected in the following theorem.Theorem 3.1.Suppose that h satisfies the conditions of Lemma 3.1, σ is bounded, and consider the BEM scheme defined by (2.8).Then for any T > 0 and p ∈ [0, 1), the following statements are valid: (1) For each m ∈ N (2) For each n where m ≥ 0 is an integer and Proof of the above theorem is lengthy and is delegated to the Appendix.We end this section with the following lemma that will be useful in our PDE approach to weak convergence rate in the following section.Lemma 3.3.Suppose that h satisfies the conditions of Lemma 3.1, σ is bounded, and consider the BEM scheme defined by (2.8).Then for any T > 0 the following statements are valid: (1) Let p ∈ [0, 1) and m ≥ 0 be an integer.For each n with K being a constant independent of n.
(3) Suppose f and h satisfy the conditions of the previous part and b ∈ C 2 b ((0, r), R).Then for each n and t ∈ [t n , t n+1 ] , for some constant K independent of n. Proof.
(1) It follows directly from the definition of µ and the hypothesis on h that for some K. Also note that if m ≥ 1 and r = ∞, there exists a K such that x m h −(2+p) ≤ Kh m−(2+p) for x ∈ [0, 1].An analogous bound can be obtained near r when r is finite.Thus, (2) Let µ s := µ(t n , X tn ; s, X s ), u := h h , and η s := H x (t n , X tn ; s, X s ).Then Ito's formula yields where M is a local martingale with M tn = 0 since η tn = 1, and Observe that the hypothesis on h implies that for some other constant K that does not depend on s.
Thus, combined with the assumption on f we arrive at for some constant K.This in particular implies M is a martingale since we can deduce from the estimates (3.6) and (3.7) that the set {f ] is a stopping time} is uniformly integrable as soon as we once again recall that |h /h| < Kh −p for some p < 1.
Hence, the claim holds in view of (3.10).
(3) Applying Ito's formula and repeating the similar estimates yields the claim.
Observe that by a Girsanov transformation we can rewrite the above expression in terms of a diffusion process satisfying the conditions in earlier sections.Indeed, defining Q on G via where and c ∈ (0, r).Thus, we may assume µ ≡ 0 and consider where X is a process satisfying Assumption 2.1, b is bounded, ε < σ < K σ and g is sufficiently regular.
for some K h and p ∈ (0, 1), g ∈ C 6 b ((0, r), R) is a bounded function with g (k) (0) = 0 (and Then Moreover, v and v t are uniformly bounded and there exists a constant K such that , where Note that u(t, 0) = 0 for t ≤ T .Moreover, it follows from Theorem 5.2 in [27] that u is the unique solution of and that sup t≤T,x∈(0,r) Also note that since we have (4.9)In particular, u t (•, 0) = 0, which in turn implies u xx (•, 0) = 0. Analogous boundary conditions also holds at r if r is finite.
Let w := u t and note that w solves (4.7) with the boundary condition w(t, 0) = 0 and w(T, •) = − 1 2 σ 2 g − σ 2 gb.Using the stochastic representation in (4.9) and analogous arguments we again arrive at w t vanishing at finite boundaries .
Using the PDE for u it is straightforward to establish that v solves (4.5) and is bounded.Moreover, as v x = hux−uh h 2 , using integration by parts we arrive at Since h (0) < ∞ and u and u xx vanish at 0 (and are jointly continuous near t = T ), there exists a neighbourhood of 0 in which |h |(y) ≤ Kh 1−p (y) ≤ K 2 y, |u(•, y)| + |u xx (•, y)| < Ky (due to Lipschitz continuity), and h(y) > cy.Thus, whenever x belongs to this neighbourhood, we have Thus, v x /h 1−p is bounded near 0. Analogous considerations when r < ∞ shows that the ratio is bounded over (0, r).
Next observe that v t is bounded since u t vanishes at finite boundaries and u tx is bounded.In particular, v t h p remain bounded near finite boundaries (uniformly in t).Multiplying (4.5) by h p and using the fact that v x /h 1−p is bounded demonstrate that sup t≤T,x∈(0,r) Finally, since v t = w h , repeating the above arguments and using the fact that w xx vanish at finite boundaries and is Lipschitz continuous in view of (4.8), we deduce v tx /h 1−p is bounded.Similar arguments (due to the boundedness of w tx = u ttx in view of (4.8) also lead to sup t≤T,x∈(0,r) In view of the above proposition, and for the convenience of the reader, we collect all the assumptions needed in Assumption 4.1 below to prove our convergence result.
Assumption 4.1.The functions σ, b, h and g satisfy the following regularity conditions.
(5) The function v defined by (4.4) belongs to C 1,4 ((0, r), R), satisfies (4.5) as well as the growth conditions for some constant K and integer m ≥ 0.
Remark 4.1.The first condition on the derivatives of h is not restrictive for practical purposes.Indeed, if a given h ∈ H ∩ C 4 ((0, r), (0, ∞)) does not satisfy this condition, one can always linearise this concave function near the boundaries at which h vanishes to obtain a new concave function satisfying the stated condition.
Theorem 4.1.Consider the BEM scheme defined by (2.8) as well as the associated error for some constant K independent of N under Assumption 4.1.
Proof.Let π 0 (s) = 1, with the convention that π k = π k (T N −1 ), and observe that Moreover, in view of (4.5) (in fact dividing both sides of the equality by σ 2 ) we have where M is a local martingale and First note that M is martingale due to (3.6) by the hypothesis on v and that h is bounded.Moreover, Lemma 3.3 shows (for a generic constant K that may change from line to line albeit remaining bounded uniformly in N ) such that x (t n , X tn ; s, X s ) where we have used the boundedness of π n+1 /π n several times and the last two lines follow from the interchange of the order of integration on the third and the fourth lines.This proves the assertion in view of Theorem 3.1 and, in particular (3.6) and (3.8), since (π n )s are non-negative and uniformly bounded, and H x ≥ 1.

Numerical analysis
This section is dedicated to the numerical experiments illustrating the above technical analysis.As we shall see, one does not really need to satisfy all the conditions assumed in Theorem 4.1 in order to achieve the advertised convergence rate in practice.The experiments below will compare the our methodology developed in this paper to standard numerical approaches for pricing of barrier options.
We shall consider the classical Black-Scholes model in the first part.As barrier option values are quite sensitive to the market skew/smile of volatility, the time-homogeneous hyperbolic local volatility model will also be studied and the corresponding results will be reported in the second part.
5.1.Black-Scholes model for barrier options.For expository purposes1 , let's assume that the asset price follows with volatility σ > 0 under risk neutral probability P. The value of the barrier option with payoff g is given by with ζ := inf{t > 0 : S t / ∈ (e , e r )}, where e represents the down barrier and e r the up barrier for −∞ ≤ < r ≤ ∞.
For the volatility to be bounded away from 0 (cf.point (2) in Assumption 4.1), we perform a change of variable X t = ln(S t ).Equations (5.1) and (5.2) then become with ζ = inf{t > 0 : X t / ∈ ( , r)} and g(x) still denotes the payoff function in x variable by an abuse of notation.
To remove the drift in (5.3), we follow the Girsanov transformation described at beginning of Section 4. We thus obtain and g(x) = g(x)e − 1 2 x .
We shall perform a path transformation method described in earlier section that either produces a recurrent process (see Theorem 2.1) or generates a transient process with infinite lifetime (see Theorem 2.2).

Specification of the recurrent transformation.
• Double barrier case with l and r finite.
We shall pick to construct the function h via (2.3).This in particular yields h (3) is bounded in ( , r), which in turn implies the boundedness of h (2)  h by means of L'Hopital's rule.In particular, (1) in Assumption 4.1 is satisfied.Double integration from (2.3) gives where a and b are given by • Single barrier case with l finite and r = +∞ We shall choose with h (x) = e −x and h (x) = −e −x .Note that with this choice of h (1) in Assumption 4.1 is only partially satisfied as |h (x)| h(x) is unbounded for x around l.
We will apply the implicit scheme (2.8) so that the price (5.4) is approximated by Remark 5.1.In the Black-Scholes model with the change of variable X t = ln(S t ), the H function is identical at each time step and needs to be computed once.In the implementation, we introduce a dense grid covering the interval (l, r), calculate the values of H on these points and H −1 is computed by piecewise constant approximation.

The transient transformation.
In the single barrier case of a down-and-out option that will constitute a part of our experiments we can also consider transformation via h(x) = x − l when l is finite and r = +∞, as in Theorem 2.2.Under Q h,x , the the process X defined in (5.5) follows One advantage of this transformation is that the inverse of the function H appearing in the implicit scheme (2.8) can be computed analytically and is given by

Down and out put option.
For a down-and-out put barrier option, the payoff is given by max(K − S T , 0)1 ζ>T where ζ := inf{t > 0 : S t / ∈ (b, +∞)}, 0 < b(= e 0, K is the option strike and T the maturity.In Black-Scholes model, standard barrier option prices are given analytically and are provided for completeness (see, e.g, p.153 of [18]).It uses a common set of factors: where N is the cumulative distribution function of a standard Normal, and H = {b, B}.For a down-and-out put barrier option with S 0 > H = b and K > H the price is given analytically by: price As mentioned at the beginning of this section, to put our methodology in perspective we have also implemented two other approaches to the numerical pricing of the barrier option: • Standard Euler without hitting probability: It consists of discretizing the SDE (5.3) according to the Euler scheme ).This numerical scheme for barrier option pricing had been studied in [16], where it was shown to have a convergence rate of O( 1 √ N ).This loss of accuracy is mainly due to the fact that it is possible for X to cross the barriers l or r at some time t between grid points t i and t i+1 and never be below the barrier at any of the dates t i for i = 1, .., N .

• Standard Euler with hitting probability:
Although this is still based on the Euler scheme simulations (5.7), it applies a further correction to remove the barrier crossing biases via the conditional no-hitting probability pi using the Brownian bridge technique (see e.g, p.169 of [16]).More precisely, the pi are defined and can be computed analytically as pi := P(∀t where the process ( X t ) 0≤t≤T is the continuous Euler scheme which interpolates ( X t i ) 0≤i≤N in the following way: It then corrects the payoff g(X T )1 ζ>T by considering instead g( X t N ) As shown in [17], this bias correction brings the convergence rate back to of order N −1 , which is the rate of weak convergence for the Euler-Maruyama scheme in the absence of killing.Moreover, in this specific Black-Scholes implementation, the simulation is exact, i.e no discretisation error occurs due to constant σ.
We shall next summarise the experiments details and comparison results.

5.2.1.
Set of parameters.The numerical experiments are conducted using the following values for the parameters: S 0 = 1, T = 1 year, l = log(b = 0.8), r = +∞ and σ = 20%.For thoroughness, we have considered in-the-money (K = 1.2), at-the-money (K = 1) and out-the-money (K = 0.9) options.To reduce statistical noise, the simulations are run with 1 million Monte Carlo paths.The benchmark price is calculated analytically with formula (5.7).
As our final results do not show any significant dependency on the moneyness of the option, we shall only report the results for at-the-money (ATM) options.In particular the discrepancy between benchmark prices and the numerical value for ATM down-and-out put options is shown in Figure 1.We have not observed any stability issues with any of our h-transformation schemes.As discussed earlier, the standard Euler with hitting probability method has no discretisation error.The discrepancy is therefore essentially the statistical noise.
Our numerical results show the rapid convergence of the numerical approximation of prices given by the recurrent and transient transforms via the implicit scheme and demonstrate clearly its effectiveness over the standard Euler scheme without hitting probability correction.This confirms the findings of our theoretical analysis even without satisfying all the conditions of Theorem 4.1.
Moreover, the prices given by the recurrent and transient transforms are quite comparable as predicted by the theoretical analysis.Figures 2 and 3 show the log-log plot of the discrepancy associated to the recurrent and transient transforms respectively for ATM down-and-out put option, respectively.The respective numerical rates of convergence observed are 0.95 and 0.9.5.3.Down and up out double barrier call option.For a down-and-up-out barrier call option, the payoff is given by max(S T − K, 0)1 ζ>T where ζ := inf{t > 0 : S t / ∈ (b, B), 0 < b(= e ) < B(= e r ) < ∞, K is the option strike and T the maturity.In Black-Scholes model, with b < S 0 < B, the price can be computed using Ikeda and Kunitomo formula (see Theorem 3.2 in [26]):  where Note that the option price is expressed as an infinite series invoving weighted normal distribution functions.However, numerical studies in [26] show the convergence of the formula is rapid and it is suggested that it suffices to calculate the leading two or three terms for most cases.Here, we use the Excel spreadsheet provided in [18] which computes each series above with n from −5 to 5.
For the standard Euler with hitting probability correction, the no-hitting probability pi is also given as an infinite series in [16] The practical studies on this formula again suggest it is perfectly sufficient for numerical purposes to calculate the leading two or three terms in most cases.To be conservative, in our experiments, the pi are estimated using n from −5 to 5.

Set of parameters.
The numerical experiments are conducted using the following values for the parameters: S 0 = 1, T = 1 year, b = 0.85, B = 1.25 and σ = 20%.For thoroughness, we consider in-the-money (K = 0.9), at-the-money (K = 1) and out-the-money (K = 1.05) options.To reduce statistical noise, the simulations are run with 1 million Monte Carlo paths.The benchmark price is calculated with formula (5.7) with truncation by keeping terms from n = −5 to n = 5.
As in the previous study, no significant difference is observed by changing the moneyness of the option.Thus, we we will again report the results pertaining to the ATM options.The discrepancies between benchmark prices and numerical methods for the ATM down-and-up-out call options are shown in Figure 4. We have not observed any stability issues with the recurrent transform method.As discussed, the standard Euler with hitting probability correction has only truncation error in the 2 up to a typographical error.computation of no hitting probability (5.8), which we believe to be negligible.The discrepancy can then be attributed essentially to the statistical noise.As in the previous experiment our numerical results align with the theoretical predictions.In particular Figure 5 shows the log-log plot of the discrepancy associated to the recurrent transform method for ATM down-and-up-out call options with numerical rate of convergence of 0.81.5.4.Time-homogeneous hyperbolic local volatility model.Since the advent of the Black-Scholes option pricing formula, the study of implied volatility has become a central preoccupation for both academics and practitioners.It is well known, actual option prices rarely conform to the predictions of explicit formulas because the idealized assumptions required for it to hold don't apply in the real world.Consequently, implied volatility (the volatility input to the Black-Scholes formula that generates the market European Call or Put price) in general depends on the strike K and the maturity of the option T .The collection of all such implied volatilities is known as the volatility surface.For example, the effect that implied volatility σ im (T, K) is a decreasing function of strike is called skew and is usually observed in equity derivatives market.This means that the underlying asset price process cannot be explained using the Black-Scholes model, for which the implied volatility does not depend on the strike.This motivates the researchers to find a convenient model for the underlying asset to evaluate contingent claim prices.Local volatility models, either parametric or non-parametric, (see e.g [13,12,34]), arguably capture the surface of implied volatilities more precisely than other approaches such as stochastic volatility models (see e.g [31,33]).Needless to say, the volatility surface has a significant impact on barrier option valuation.Indeed, the barrier hitting probability depends strongly on the dynamics of the volatility of the spot pricess (see, e.g., [15]).
For our analysis, we consider the time homogeneous hyperbolic local volatility model (HLV), which is widely used in quantitative finance community to capture the market skew.It corresponds to a parametric local volatility-type model in which the dynamics of the underlying under the risk neutral measure P is given by Here ν > 0 is the level of volatility, β ∈ (0, 1] is the skew parameter.First introduced in [25] it behaves similarly to the Constant Elasticity of Variance (CEV) model and has been used for numerical experiments in, e.g., [20,21,19].A practical advantage of this model is that zero is not an attainable boundary, which in turn avoids some numerical instabilities present in the CEV model when the underlying asset price is close to zero (see e.g.[4]).It corresponds to the Black-Scholes model for β = 1 and exhibits a skew for the implied volatility surface when β = 1. Figure 6 illustrates the impact of the parameter β on the skew of the volatility surface.We observe that the skew increases significantly with decreasing value of β.For example with ν = 0.3, β = 0.2, the difference in volatility between strikes at 50% and at 100% is about 15%.5.4.1.Down and up out double barrier call option.In this implementation we shall set h(x) = (x − l)(r − x) and the associated BEM scheme will be then solved using bisection method with Octave vectorization for faster code execution.Consequently the price is approximated by For comparison, we compute also the numerical price given by the standard Euler scheme with hitting probability.The scheme is given by equation (5.7) and the no hitting probability formula by (5.8), where σ is computed using the parametric local volatility function (5.9).Experiment details and comparison results are described below.For thoroughness, we consider in-the-money (K = 0.9), at-the-money (K = 1) and out-the-money (K = 1.05) options.The benchmark prices for each numerical method are computed by the method itself with very dense time grid and high number of Monte Carlo paths.
In this case we observed some differences regarding the moneyness of the option in our numerical results.More precisely, the method performed relative poorly for the ATM option.For this reason we report below the results in all three cases and provide an explanation for the seemingly poor performance for the ATM option.
The discrepancies between benchmark prices and numerical methods for ITM, ATM and OTM double barrier call options are shown respectively in Figures 7, 8 and 9.We have not observed any stability issues with the recurrent transform scheme.Interestingly, our recurrent transformation has a much smaller error than the explicit Euler method with hitting probability correction when the number of discretisations is reasonably large.More importantly, this outperformance is still valid even if the number of Monte Carlo simulations for the explicit Euler method is increased five times.Having said that, one should still treat such a conclusion with caution as our benchmark price and hitting probabilities are calculated by applying a truncation and, thus, is subject to error.Nevertheless, the outperformance is still promising as our truncation is no coarser than the common industry practice.Figures 10,11 and 12 show, respectively, the log-log plot of the discrepancy associated to the recurrent transform method for ITM, ATM and OTM double barrier call options.The numerical rate of convergence are respectively 0.91, 0.63 and 1, using 2 × 10 5 Monte Carlo simulations.Although the rate of convergence for the ATM option is far from the theoretical rate of 1, a closer look at Figure 8 reveals a clue.Note that the error of approximation converges very rapidly to zero after a few iterations and further discretisations do not significanly alter the already very small error term.This indicates that the observed error in this case can be mostly attributed to the statistical noise and the simple regression to obtain the convergence rate does not work well.
When we run the same experiment for the Euler scheme with hitting probability correction with 2 × 10 5 Monte Carlo simulations, we observe a similar drop in the performance and the convergence rates are found to be 0.50, 0.59 and 0.61, respectively.However, the convergence rates for the latter scheme increases to 0.83, 0.83 and 0.77, respectively, when the number of simulations are increased five-fold.

Conclusion
We have introduced a novel backward Euler-Maruyama method to increase the weak convergence rate of approximations in the presence of killing.The numerical experiments confirm our theoretical prediction that the convergence rate is of order 1/N , where N is the number of discretisations.Moreover, the numerical studies suggest that one does not need a large N to obtain a sufficiently  close approximations as all numerical studies indicate errors terms diminishing very rapidly with a small number of iterations.The numerical experiments also suggested our method outperforming the Brownian bridge method in certain cases although such a statement does not currently have any theoretical backing.However, we believe that the method developed in this paper will perform better when applied to a higher order Euler-scheme such as the Milstein scheme.Such investigations are left for future research.
Moreover, a close look into our technical analysis reveals that our convergence result does not depend heavily on the one-dimensional nature of the problem.In particular it is relatively clear how to obtain a version of Theorem 2.3 in the multidimensional case using well-established potential theoretic arguments.However, our main obstacle in not being able to immediately obtain a multidimensional version of Theorem 4.1 is the absence of a systematic study of recurrent transformations in higher dimensions.Such a study and its applications to the Euler methods for killed diffusions will be the subject of future research.
Appendix A. Proof of Theorem 3.1 Proof will be divided into several steps considering first the case of r = ∞ and making use of the comparison Lemma 3.2.In what follows K denotes a generic constant independent of N .
(1) First suppose r = ∞.Since 1/h is decreasing, Lemma 3.2 and Theorem 2.3 imply Moreover, Lemma 3.2 also yields where Y is a process that shares the same law with the process in Theorem 2.3 with c = c 1 and the last inequality follows from Theorem 2.3.Similarly, by considering instead the time change , where Y is a process such that Y tn = X tn and with β being a standard Brownian motion.Consequently, Theorem 2.3 yields ) Now consider the case r < ∞ and set x 1 := inf{x ≥ 0 : h (x) = 0} and x 2 := inf{x ≥ x 1 : h (x) < 0}.Then, there exist functions h 1 and h 2 such that h = h 1 h 2 , h 1 (resp.h 2 ) is non-decreasing (resp.non-increasing) and constant on (x Applying Ito formula to ((x 2 − Y 1 t ) + ) 2 and ((x 2 − X t ) + ) 2 , the comparison theorem employed in Lemma 3.2 shows that P h,X 0 ( Y 1 t ∧ x 2 ≤ X t ∧ x 2 , t ≤ T ) = 1.An analogous argument also shows that As h 1 is non-decreasing, h 2 is non-increasing and h h 1 (resp.h h 2 ) is constant on (0, x 2 ) (resp.(x 1 , r)), the above comparisons imply that Analogous considerations also yield ess sup τ ∈Tn E h,X 0 1 h ( X τ ) F tn < ∞.
(3) We shall now show the boundedness of the moments.Note that there is nothing to show when r < ∞.So, let's assume that r = ∞.Recall that 3 Although h2 does not quite satisfy the condition therein, we obtain the result that we need by a change of scale and considering instead the function h defined by h(x) = h2((r − x) + ).Note that h is still continuously differentiable.
Thus, E h,X 0 ( X t ) ≤ E h,X 0 ( X tn ) + K(t − t n ) for some K due to the boundedness of σ and h as well as the uniform bound on the inverse moment of h( X t ).This shows that sup t≤T,M E h,X 0 ( X t ) ≤ X 0 + KT.where Z is a local martingale.

Now
Next observe that for m ≥ 1 as h(0) = 0, h (0) > 0 and h /h ≤ 1 x .The last identity follows from the fact that Observe that (τ k ) k≥1 , where τ k := inf{t ≥ t n : X t ≥ k} is a localising sequence for Z.Therefore, a standard localisation argument, (A.1) and (A.2) together imply for t ∈ (t n , t n+1 ] E h,X 0 ( X m+1 t ) ≤ E h,X 0 ( X m+1 tn ) + (t − t n )KE(m − 1), in view of the Fatou's lemma for some constant K, which in turn yields Finally, note that this in particular implies that Z is a true martingale.Thus, for τ ∈ T n and m ≥ 2 Taking conditional expectations show E h,X 0 X m τ |F tn ≤ X m tn + KE h,X 0 t n+1 tn X m−1 t dt F tn , (A.3) yielding (3.7).To establish (3.8) we need the following lemma.
Lemma A.1.Suppose that h satisfies the conditions of Lemma 3.1, σ is bounded, and consider the BEM scheme defined by (2.8).For any p ∈ [0, 1), any n and t n ≤ s ≤ t < t n+1 we have E h,X 0 h −p ( X t )|F s ≤ h −p ( X s ) exp(K(t − s)), for some constant K > 0 that is independent of n.
Proof.Let µ t := µ(t n , X tn ; t, X t ).A straightforward application of Ito's formula yields Next note that τ k := inf{t ≥ t n : X t < 1/k} is a localising seqeunce for M .Thus, using the optional stopping theorem and Fatou's lemma and monotone convergence we arrive at E h,X 0 h −p ( X t )|F s ≤ h −p ( X s ) + KE h,X 0 t s h −p ( X u )du F s , t n ≤ s ≤ t ≤ t n+1 , for some constant K in view of the boundedness of σ.We deduce the claim by Gronwall's lemma.Now we return to the proof of the estimate (3.8).Observe that the hypothesis on h implies 1 − exp((s − t n )σ 2 ( X tn ) h 2h ( X tn )) ≤ K T N 1 h p ( X tn ) for some K > 0. Without loss of generality let's also suppose that h ≤ 1.Thus, where the second line follows from Lemma A.1 and that H x ≥ 1.
Next suppose m ≥ 1.Note that the calculations similar to the ones leading to (A.3) imply that E h,X 0 X m t |F tn ≤ X m tn + KE h,X 0 t n+1 tn X m−1 s ds F tn , Thus, the elementary inequality x m−1 ≤ 1 + x m and Gronwall's lemma yield Therefore, where the last line follows from (3.10).
Combining above estimates, we arrive at the claimed result via (3.6).
1 , r) (resp.(0, x 2 )).Let's define the processes Y i , where Y i 0 = X 0 and ) + c i dt, t ∈ (t n , t n+1 ]. ( X tn ) H 2 x (t n , X tn ; t, X t ) h h ( X t ) + µ(t n , X tn ; t, X t ) dt n , X tn ; t, X t ) dt, X tn ) H 2x (t n , X tn ; t, X t )ph −p ( X t ) 2µ t h ( X t ) + h ( X t ) 2 ( X t ) dtwhere M is a local martingale and α 1 and α 2 are positive constants depending on the bounds on h and h since µ t > c 1 (resp.µ t < c 2 ) whenever h ( X t ) > 0 (resp.h ( X t ) < 0) by Lemma 3.1 and h never vanishes at the same time as h.Thus, there exists a constant K that depends only on h, p and c 1 and c 2 such that dh −p ( X t ) ≤ dM t + K σ 2 ( X tn )h −p ( X t ) h n , X tn ; t, X t ) dt since −α 1 x + (1 − p)α 2 x 2 is bounded from below.