Exact simulation of first exit times for one-dimensional diffusion processes

TL;DR

This paper introduces an exact simulation algorithm for first exit times of one-dimensional diffusion processes, avoiding errors from discretization by using advanced stochastic techniques.

Contribution

A novel acceptance-rejection algorithm that precisely simulates first exit times without approximation errors for one-dimensional diffusions.

Findings

Algorithm accurately simulates first exit times

The method is efficient and theoretically justified

Numerical examples demonstrate practical effectiveness

Abstract

The simulation of exit times for diffusion processes is a challenging task since it concerns many applications in different fields like mathematical finance, neuroscience, reliability... The usual procedure is to use discretiza-tion schemes which unfortunately introduce some error in the target distribution. Our aim is to present a new algorithm which simulates exactly the exit time for one-dimensional diffusions. This acceptance-rejection algorithm requires to simulate exactly the exit time of the Brownian motion on one side and the Brownian position at a given time, constrained not to have exit before, on the other side. Crucial tools in this study are the Girsanov transformation, the convergent series method for the simulation of random variables and the classical rejection sampling. The efficiency of the method is described through theoretical results and numerical examples.

Equations

$$d{X}_t=\mu (X_t)dt+\sigma (X_t)d{B}_t,X_0=x\in (a,b),$$

$${\tau }_{a,b}(X):=\inf \{t>0:X_t\notin [a,b]\}.$$

$$d{X}_t=\mu (X_t)dt+d{B}_t,X_0=0,$$

$$\displaystyle\mathbb{E}[F(B_{t},\,t\leq\tau_{a,b}(B))|B_{0}=y]=\mathbb{E}\Big{[}F\Big{(}\frac{b-a}{2}\ B_{t}+\frac{a+b}{2},\,t\leq\tau_{-1,1}(B)\Big{)}\Big{|}B_{0}=x\Big{]},$$

$$\Big{(}\tau_{a,b}(B^{y}),B_{\tau_{a,b}(B^{y})}^{y}\Big{)}\stackrel{{\scriptstyle(d)}}{{=}}\Big{(}\frac{(b-a)^{2}}{4}\,\tau_{-1,1}(B^{x}),\frac{b-a}{2}B^{x}_{\tau_{-1,1}(B^{x})}+\frac{a+b}{2}\Big{)}.$$

$${\tau }^x=\inf \{t>0:B_t^x\notin I\}\text{mbox}withI=[-1,1].$$

$$\mathbb{P}(B_{\tau }^x=1)=\frac{\mathfrak{s}(-1)-\mathfrak{s}(x)}{\mathfrak{s}(-1)-\mathfrak{s}(1)}=\frac{x+1}{2},\forall x\in [-1,1].$$

$$\mathbb{E}[\psi (B_t^x),{\tau }^x>t]={\int }_{-1}^1\psi (y)p(t,x,y)\mathrm{d}y,x\in [-1,1].$$

$$\left\{\displaystyle\begin{array}[]{l}\frac{\partial p}{\partial t}(t,x,y)=\frac{1}{2}\,\frac{\partial^{2}p}{\partial y^{2}}(t,x,y)\\[5.0pt] p(t,x,y)\to\delta_{x}(y)\quad\mbox{as}\ t\to 0,\\[5.0pt] p(t,x,y)=0\quad\mbox{for}\ y\in\{-1,1\},\ t>0.\end{array}\right.$$

$$p(t,x,y)=\frac{1}{\sqrt{t}}\sum_{n=-\infty}^{+\infty}\Big{\{}\phi((x-y-4n)/\sqrt{t})-\phi((x+y-2-4n)/\sqrt{t})\Big{\}}$$

$$p(t,x,y)=\sum_{n\geq 1}\exp\Big{(}-\frac{n^{2}\pi^{2}t}{8}\Big{)}\sin\Big{(}\frac{n\pi}{2}\ (x+1)\Big{)}\sin\Big{(}\frac{n\pi}{2}\ (y+1)\Big{)}.$$

$$f(y)=\sum _{n\ge 0}{f}_n(y).$$

$$f(y)\le \kappa h(y),y\in \mathbb{R}.$$

$$|R_n(y)|\le r_n,\forall y\in \mathbb{R}\text{mbox}where{R}_n(y):=\sum _{k\ge n+1}{f}_k(y).$$

$$p(t,x,y)$$

$$a_{-n}(t,x,y)<0\text{mbox}and{a}_n(t,x,y)>0,\forall n\ge 1,\forall (x,y)\in [-1,1{]}^2.$$

$$(x+y-2-4n)-(x-y-4n)=(x+y-2+4n)-(x-y+4n)=2y-2.$$

$$(x+y-2-4n)<(x-y-4n)<0<(x+y-2+4n)<(x-y+4n),$$

$$n_0=\lfloor \sqrt{t}/4\rfloor +2.$$

$$R_n(y):=\sum _{k\ge n+1}(a_k(t,x,y)+a_{-k}(t,x,y)),n\ge 0.$$

$$\frac{1}{\sqrt{t}}\sum _{k\ge n+1}\varphi ((\alpha +4k)/\sqrt{t})\le \frac{1}{\sqrt{t}}{\int }_n^{\infty }\varphi ((\alpha +4z)/\sqrt{t})\mathrm{d}z,$$

$$\frac{1}{\sqrt{t}}\sum _{k\ge n+1}\varphi ((\alpha -4k)/\sqrt{t})\le \frac{1}{\sqrt{t}}{\int }_n^{\infty }\varphi ((\alpha -4z)/\sqrt{t})\mathrm{d}z,$$

$$|R_n(y)|$$

$$\le \frac{1}{\sqrt{t}}{\int }_n^{\infty }\varphi ((x-y-4z)/\sqrt{t})\mathrm{d}z+\frac{1}{\sqrt{t}}{\int }_n^{\infty }\varphi ((x+y-2+4z)/\sqrt{t})\mathrm{d}z.$$

$$|R_n(y)|$$

$$\displaystyle\leq r_{n}(t):=\frac{1}{4}\Big{(}1-{\rm erf}\Big{(}\frac{4n-2}{\sqrt{2t}}\Big{)}\Big{)}.$$

$$p(t,x,y)\le \kappa (t,x)h_{t,x}(y),\forall y\in [-1,1],$$

$$h_{t,x}(y)=\frac{1}{\sqrt{t}}\varphi ((x-y)/\sqrt{t}),$$

$$\frac{p(t,x,y)}{h_{t,x}(y)}\le 1+\sum _{n\ge 1}\frac{a_{-n}(t,x,y)}{h_{t,x}(y)}+\frac{a_n(t,x,y)}{h_{t,x}(y)}.$$

$$\frac{p(t,x,y)}{h_{t,x}(y)}$$

$$\frac{p(t,x,y)}{h_{t,x}(y)}\le 1+\sum _{n=1}^{n_0}\frac{\varphi ((x-y-4n)/\sqrt{t})}{\varphi ((x-y)/\sqrt{t})}.$$

Full text

Key words and phrases: Exit time, Brownian motion, diffusion processes, Girsanov’s transformation, rejection sampling, exact simulation, randomized algorithm, conditioned Brownian motion.

2010 AMS subject classifications: primary 65C05; secondary: 65N75, 60G40.

Introduction

First exit time distributions for general stochastic processes are of prime importance in many contexts. In mathematical finance, they permit to appreciate the risk of default for a given path-dependent option; in neuroscience, they characterize for instance the interspike time distribution… Since diffusion processes (solutions of stochastic differential equations) form an important family of stochastic processes, we aim to describe quite precisely their exit times. Unfortunately a simple and explicit expression of their distribution is not generally available which leads us to consider numerical approximations. Our task is to point out an algorithm which permits to simulate the exit time of the diffusion. This objective was already concerned by several studies introducing a discretization scheme for the corresponding stochastic differential equation. Most of them are based on improvements of the classical Euler scheme (see for instance [7], [12], [11]) which essentially consist in reducing the error stem from the approximation procedure. Let us also note another point of view which consists in approximating the probability density function of the exit time and therefore to approximate the solution of an integral equation [19].

Our approach is completely different since we emphasize an exact simulation procedure: the distribution of the random outcome of the algorithm is identical to the distribution of the first exit time of the diffusion process. For such simulation methods based on an acceptance-rejection sampling, the challenge is to describe and reduce if possible the time consumption of the simulation. Beskos & Roberts, in their founding work [5], already introduced the exact simulation for diffusion paths on some fixed time interval. Meanwhile several modifications of this algorithm have been proposed [3, 4, 15]. The basic idea of the rejection sampling is to sequentially observe independent random objects generated according to a proposal distribution until a condition is satisfied. That means that each object is accepted or rejected according to a certain probability weight.

For diffusion processes on a fixed interval $[0,T]$, the proposed object pointed out by Beskos & Roberts is the Brownian paths (easy to simulate) and the rejection weight can be computed using the Girsanov transformation which essentially requires the simulation of Brownian bridges. In a previous work [13], the authors proposed a similar procedure in order to exactly simulate the first passage time of a one-dimensional diffusion through a given threshold: the proposal object is then the Brownian first passage time (inverse gaussian random variable) and the rejection weight requires the simulation of Brownian bridges conditioned to stay under a given threshold (Bessel processes). In order to adapt such a procedure to exit times, it suffices to generate Brownian exit times and to use the rejection weight suggested by the Girsanov transformation. Unfortunately this random weight requires the simulation of a Brownian bridge conditioned to stay in a given interval, which corresponds to a SDE with singular coefficients: there is no way to exactly simulate such a stochastic paths. Such an intuitive generalization leads therefore to a deadlock.

That’s why we aim to present a quite different rejection sampling based on a similar concept (Girsanov’s transformation) and avoiding the simulation of the whole conditioned Brownian paths. Let us consider the stochastic process $(X_{t},\ t\geq 0)$, solution of the SDE:

[TABLE]

where $(B_{t},\ t\geq 0)$ stands for the standard one-dimensional Brownian motion. We denote by $\tau_{a,b}$ the first time the diffusion exits from the interval $[a,b]$:

[TABLE]

The aim of the study is to propose an efficient algorithm in order to simulate the first exit time $\tau_{a,b}(X)$. Under suitable conditions, the Lamperti transform permits to reduce the area of investigation to the constant diffusion case: $\sigma(x)\equiv 1$ providing to change the interval $[a,b]$. That’s why we shall only focus our attention on the diffusion process:

[TABLE]

for $t\leq\tau_{a,b}(X)$ with $a<0<b$. The exact algorithm presented here (Section 3) is essentially based on the Girsanov transformation which permits to relate the diffusion (0.2) to a standard one-dimensional Brownian motion. Let us roughly describe the crucial tools needed to execute such an algorithm. The proposal random variable is the Brownian exit time of the interval $[a,b]$, denoted $\mathcal{T}_{\rm prop}$. The weight used in order to accept or reject this proposal depends both on the conditional distribution of the Brownian paths $(B_{t})$ given that $\tau_{a,b}(B)=\mathcal{T}_{\rm prop}$ and on a Poisson process defined on the time axis $\mathbb{R}_{+}$ which is independent of the Brownian motion $(B_{t})$ and whose events occur at time $T_{1},\ldots,T_{n},\ldots$ Then conditionally to $\tau_{a,b}(B)=\mathcal{T}_{\rm prop}$, only the values of $B_{T_{k}}$, for any $T_{k}\leq\mathcal{T}_{\rm prop}$, and $B_{\mathcal{T}_{\rm prop}}$ (dots in Figure 1) are involved in the specific conditions associated to the acceptance of $\mathcal{T}_{\rm prop}$.

The exact simulation of $\tau_{a,b}(X)$ therefore requires to describe:

The conditional distribution of $B_{t}$ for a given $t>0$, given that the first exit of the interval $\tau_{a,b}(B)$ is larger than $t$ (Section 1). Using the classical convergent series method, we propose an algorithm denoted by CONDITIONAL_DISTR (x,[a,b],t) for the simulation of the constrained Brownian motion $B_{t}$ with initial condition $x$. 2. 2.

The distribution of the Brownian exit time of the interval $[a,b]$ as the Brownian trajectory starts in $x$. The corresponding algorithm, presented in Section 2, is denoted BROWNIAN_EXIT_ASYMM  (x,[a,b]) .

The two main algorithms (DET) and (GDET) presented in Section 3 use these two previous algorithms as basic components. They permit the exact simulation of diffusion exit times under weak conditions on the coefficient $\mu$. Efficiency results are pointed out (Theorem 3.2 and Proposition 3.4) and numerical illustrations complete the study in Section 3.4

1 Brownian motion constrained to stay in an interval $[a,b]$

Let us first consider the Brownian motion in the space interval $[a,b]$. We need to describe the distribution of $B_{t}$ for a given $t>0$, given that the first exit of the interval $\tau_{a,b}(B)$ is larger than $t$. We fix the starting position of the Brownian motion $B_{0}=y\in(a,b)$. Due to both the translation invariance and the scaling property of the Brownian motion, we get

[TABLE]

where $x=\frac{2y-a-b}{b-a}$. So if we denote $B^{x}$ the Brownian motion starting in $x$, we deduce easily the following distribution identity

[TABLE]

In other words, it suffices to study precisely the first Brownian exit problem from the normalized interval $[-1,1]$ with the initial condition $B_{0}=x\in[-1,1]$. For notational convenience, we denote by

[TABLE]

Let us first emphasize a classical result on the exit position: since $(B_{t},\ t\geq 0)$ is a martingale, the identity function $\mathfrak{s}(x)=x$ corresponds to the scale function of the Brownian motion. The optimal stopping theorem permits therefore to describe the probability to exit from the interval on one particular side:

[TABLE]

1.1 Distribution of $B_{t}$ given $\tau>t$

Let us now focus our attention to the distribution of the constrained Brownian motion. We introduce the killed Brownian motion: as soon as the process hits the boundary of the interval, it jumps to a cemetery point $x^{*}\notin[-1,1]$. Indeed let us denote by $p(t,x,dy):=\mathbb{P}(B^{x}_{t}\in dy,\tau>t)$ the transition probability of the Brownian motion joint with $\tau>t$. It can be written as $p(t,x,y)=q(t,x,y)\mathbb{P}(\tau^{x}>t)$ where $q(t,x,y)$ is the transition probability of the Brownian motion conditioned to stay in the interval $[-1,1]$ till time $t$. For any non negative or bounded measurable function $\psi$, we get

[TABLE]

It is well known (see, for instance [14] Section 4.11 or [8] Section 5.7) that $(t,y)\mapsto p(t,x,y)$ satisfies an Initial-Boundary value problem associated to the heat equation:

[TABLE]

Milstein and Tretyakov [18] recall that the solution of (1.4) takes two different expressions. The first one is obtained by the method of images and the second one is based on a spectral decomposition of the heat equation. We introduce the standard gaussian pdf and cdf: $\phi(x)=\frac{1}{\sqrt{2\pi}}\,e^{-x^{2}/2}$ and $\Phi(x)=\int_{-\infty}^{x}\phi(y)\mathrm{d}y$. Then

[TABLE]

which is a convenient formula for small values of $t$. For large times $t$, we would prefer the second formula:

[TABLE]

Let us fix $(t,x)$: since we get an explicit expression of $p(t,x,y)$, we can point out exact simulation algorithms based on the acceptance/rejection method even if $p(t,x,y)$ is not a probability density function. Indeed it suffices to divide by $\mathbb{P}(\tau^{x}>t)$ in order to obtain such a density (cf. Section 2.1). The method applied here is a particular acceptance/rejection method presented in [10] as the convergent series method.

1.2 Inequalities related to the series expansion

Let us consider $f$ a non negative function satisfying a series expansion:

[TABLE]

We assume that $I(f):=\int_{\mathbb{R}}f(y)\,dy<\infty$. The convergent series method permits to simulate a random variable with the corresponding probability distribution function $f/I(f)$. It requires two important features:

• •

the existence of both a pdf $h$ and a constant $\kappa>0$ such that

[TABLE]

• •

the existence of a positive sequence $(r_{n})_{n\geq 0}$ which converges toward [math] and satisfies the following reminder upper-bound

[TABLE]

CONVERGENT SERIES METHOD

Initialization. $\mathcal{N}_{c}=0$.

Step 1. Generate a random variable $Y$ with density $h$.

Step 2. Generate a random variable $U$ uniformly distributed on $[0,1]$ and define $W=\kappa Uh(Y)$.

Step 3. Initialization. $\mathcal{N}_{l}=0$, $S=0$, ${\rm Test}=0$.

Step 4. While $({\rm Test}=0)$ do:

• •

$S\leftarrow S+f_{\mathcal{N}_{l}}(Y)$

• •

${\rm Test}=1_{\{|S-W|>r_{\mathcal{N}_{l}}\}}$

• •

$\mathcal{N}_{l}\leftarrow\mathcal{N}_{l}+1$ and $\mathcal{N}_{c}\leftarrow\mathcal{N}_{c}+1$.

Step 5. If $W\leq S$ then $X=Y$ otherwise go to Step 1.

Outcome: the random variable $X$ with density $f/I(f)$ and the global number of terms of the series expansions $\mathcal{N}_{c}$ used.

In order to simulate $q(t,x,y)$, i.e. the conditional distribution of $B_{t}$ given $t<\tau$, we need to choose the density $h$, to point out some constant $\kappa$ and the sequence $(r_{n})_{n\geq 0}$ for each explicit expression of $f$ given by (1.5) and (1.6). Let us just note that it is challenging to find the smallest constant $\kappa$ as possible since the number of iterations of Step 1 in order to get an outcome is geometrically distributed with average $\kappa/I(f)$.

Some comments on the series expansion (1.5)

We note that (1.5) is not an alternating series and can be rewritten as

[TABLE]

where $a_{n}(t,x,y)=\frac{1}{\sqrt{t}}(\phi((x-y-4n)/\sqrt{t})-\phi((x+y-2-4n)/\sqrt{t}))$. We first observe that

[TABLE]

Moreover both increments of the function $\phi$ are computed on intervals whose length does not depend on the variables $x$ and $n$:

[TABLE]

For $n$ large enough, the increments are decreasing and therefore $a_{n}(t,x,y)+a_{-n}(t,x,y)$ becomes negative. Indeed by considering the four terms:

[TABLE]

it is straightforward to see that the smallest one in absolute value is $|x+y-2+4n|$. Since the change of convexity of the curve $x\mapsto\phi(x)$ appears for $x=1$, the sum $a_{n}(t,x,y)+a_{-n}(t,x,y)$ is negative as soon as $x+y-2+4n\geq\sqrt{t},$ which is satisfied in particular if $n\geq n_{0}$ with

[TABLE]

Since the terms of the series are negative for $n\geq n_{0}$, the series (1.5) is obviously not an alternating series.

Bound of the series reminder in (1.5)

Let us fix $(t,x)\in\mathbb{R}_{+}^{*}\times[-1,1]$. We introduce the reminder of the series (1.5):

[TABLE]

In order to apply the convergent series method, we need to bound this reminder. Since $\phi$ is a decreasing function, we obtain the following bound:

[TABLE]

as soon as $\alpha+4n\geq 0.$ Moreover, by symmetry,

[TABLE]

as soon as $\alpha-4n\leq 0.$ If $\alpha=x-y$ or $\alpha=x+y-2$, then the inequalities are satisfied for any $n\geq 1$. Finally, for $n\geq 1$, we get

[TABLE]

Let us now observe that we can obtain a bound which is uniform with respect to both the variable $x$ and $y$.

[TABLE]

Let us note that this upper-bound is efficient for small values of $t$.

Proposal distribution for the rejection method using (1.5)

The aim is to use the acceptance/rejection algorithm in order to simulate a random variable with the target density $q(t,x,y)=p(t,x,y)/\mathbb{P}(\tau^{x}>t)$ where $p(t,x,y)$ is given by (1.5). That’s why we are looking for a proposal pdf $h$ and a constant $\kappa(t,x)$ independent of $y$ satisfying

[TABLE]

as explained in (1.7). We suggest here the following choice of the proposal distribution:

[TABLE]

that means that we choose a gaussian distribution centered in $x$ with variance $t$: it corresponds to the distribution of $B_{t}$ without any conditioning. The rejection method shall permit to go from this initial distribution to the conditional distribution with respect to the event $\{\tau>t\}$. The algorithm is therefore particularly efficient if the considered event is satisfied with a large probability, this is namely the case for small values of $t$.

Let us now focus our attention to the upper-bound

[TABLE]

Since the sum $a_{-n}(t,x,y)+a_{n}(t,x,y)$ is negative as soon as $n\geq n_{0}$ (see the definition of $n_{0}$ in (1.10)), we obtain

[TABLE]

Moreover $a_{n}(t,x,y)$ is negative, and therefore

[TABLE]

In fact all the ratios considered just above have the same nice property: they are smaller than $1$. It suffices to notice that $(x-y)\in]-2,2[$ implies $\phi((x-y)/\sqrt{t})\geq\phi(\alpha/\sqrt{t})$ for any $\alpha\notin]-2,2[$. Hence, for $\alpha=(x-y-4n)$ with $n\geq 1$, the ratios in the r.h.s. of (1.11) are smaller than 1. Using the definition of $n_{0}$ in (1.10), (1.11) becomes

[TABLE]

Let us note that this bound is uniform with respect to the variable $x$ (starting value of the Brownian paths).

Some comments on the series expansion (1.6)

Since $q(t,x,y)$ can be characterized by two different series expansions, it is useful to understand which kind of acceptance/rejection algorithm each series produces. Let us now focus on the second series (1.6):

[TABLE]

where $\beta_{n}(t)=e^{-n^{2}\pi^{2}t/8}$ and a_{n}(x,y)=\sin\Big{(}\frac{n\pi}{2}\ (x+1)\Big{)}\sin\Big{(}\frac{n\pi}{2}\ (y+1)\Big{)}. Let us first notice that, for fixed $t>0$, $(\beta_{n}(t))_{n\geq 1}$ is a decreasing sequence of positive real numbers which converges towards [math]. Since $a_{n}$ is defined as the product of two sine functions, one idea to bound the reminder of the series is to use the Abel transform and to prove that $\sum_{n=1}^{N}a_{n}(x,y)$ is bounded as $N$ tends to infinity. It is quite easy to obtain some bound when $x\neq y$ but for the particular case $x=y$, the sequence of partial sums is increasing (sum of squares) and tends to infinity. In conclusion, we cannot apply technics based on the Abel transform for the series (1.6). We shall therefore present an other approach.

Bound of the series reminder in (1.6)

Let us define the reminder of the series (1.6):

[TABLE]

Let us find a bound of this reminder. Since $|a_{k}(x,y)|\leq 1$ for any $k\in\mathbb{N}$ and any $(x,y)\in]-1,1[^{2}$, we obtain

[TABLE]

Moreover, $u\mapsto e^{-u^{2}\pi^{2}t/8}$ is a decreasing function on $\mathbb{R}_{+}$, and therefore we get the following bound

[TABLE]

Note that this bound is sharp for large value of $t$.

Proposal distribution for the rejection method using (1.6)

The aim is to use an acceptance /rejection method in order to simulate a random variable with the target density $p(t,x,y)/\mathbb{P}(\tau>t)$ based on the equation (1.6), as already done for (1.5). We are looking for a probability distribution function $h(y)$ and a constant $\kappa(t,x)>0$ independent of $y$ such that:

[TABLE]

as explained in (1.7). Our particular choice for the function $h$ is

[TABLE]

which corresponds to the invariant probability measure of the Brownian motion conditioned to stay in the interval $[-1,1]$. Let us now determine a constant $\kappa(t,x)$. We have

[TABLE]

where

[TABLE]

Using the formula $a^{n}-b^{n}=(a-b)(a^{n-1}+a^{n-2}b+\ldots+b^{n-1})$ applied to the sine function, we observe that $|s_{n}(x)|\leq n$ for all $x\in[-1,1]$. Hence

[TABLE]

Moreover the function $x\mapsto x^{2}\,e^{-x^{2}\pi^{2}t/8}$ is decreasing as soon as $x\geq\frac{2\sqrt{2}}{\pi\sqrt{t}}$, we obtain the following upper bound:

[TABLE]

where $n_{0}=\left\lfloor\frac{2\sqrt{2}}{\pi\sqrt{t}}\right\rfloor+1$.

The integration by parts formula and a change of variables lead to

[TABLE]

since $\sqrt{\frac{2}{\pi t}}\leq\sqrt{\pi}n_{0}/2\leq n_{0}$. We just note that

[TABLE]

as $t$ becomes large. Moreover

[TABLE]

Finally we can choose the constant $\kappa(t,x):=\frac{4}{\pi}(C_{1}(t)+C_{2}(t))\ \sin(\frac{\pi}{2}\,(x+1)).$ Let us note that $\kappa(t,x)\sim\frac{4}{\pi}\sin(\frac{\pi}{2}\,(x+1))e^{-\pi^{2}t/8}$ as $t$ tends to infinity. Since the averaged number of iterations in the acceptance/rejection method corresponds to $\kappa(t,x)/\mathbb{P}_{x}(\tau>t)$ and since both $\kappa(t,x)$ and $\mathbb{P}_{x}(\tau>t)$ are of the same order in the large time limit, the efficiency of the algorithm using (1.6) still remains strong when the time variable enlarges. Indeed let us decompose the following probability:

[TABLE]

where

[TABLE]

Using the definition (1.14), we obtain

[TABLE]

The inequality $|s_{n}(x)|\leq n$ for any $x\in[-1,1]$ and $n\geq 1$ leads to

[TABLE]

We observe that each term of the series converges in a monotonous way towards [math] as $t$ tends to $\infty$ which implies that the ratio tends towards $1$ by the Lebesgue theorem. We therefore deduce that for any $\theta>1$, there exists $t_{\theta}>0$ such that the average number of iterations in the acceptance/rejection algorithm is smaller than $\theta$ as soon as $t\geq t_{\theta}$.

1.3 Algorithm and numerics

Let us now describe the algorithms used in the following in order to simulate $q(t,x,dy)$, the conditional distribution of the Brownian motion at time $t$ given $\tau>t$. We know that $p(t,x,dy)$ admits two different series expansion presented in (1.5) and (1.6). As introduced in Section 1.2, the classical convergent series method requires both a proposal distribution $h$, satisfying the inequality (1.7) that is $p(t,x,y)\leq\kappa(t,x)h_{t,x}(y)$ (also necessary for classical acceptance/rejection methods), and a precise description of the series reminder characterized by the sequence $(r_{n}(t))_{n}$ and $(\hat{r}_{n}(t))_{n}$, see (1.8).

We just recall the results obtained in the previous section:

It is straightforward that the convergent series method using (1.5) is convenient for small values of $t$ while (1.6) is rather convenient for large $t$. That’s why we choose a threshold $t_{\rm c}>0$ (threshold for the conditional distribution) such that (1.5) is used for $t\leq t_{\rm c}$ and (1.6) otherwise. For practical purposes, we fix $t_{\rm c}=0.7$, this choice is motivated by the curves in Fig. 2 and will be held for all numerical illustrations presented in this study. Applying the convergence series algorithm either for small times or large times permits to obtain the simulations presented in Fig. 3. We observe that, even if $x\neq 0$, the condition $t<\tau$ leads to a distribution which looks like symmetric as $t$ becomes large and which converges towards $h_{t,x}(y)$ the invariant measure of the diffusion conditioned to stay in the interval $[-1,1]$.

Efficiency of the algorithms

Let us just recall classical results concerning the efficiency of the convergent series method introduced in Section 1.2. We introduce $\mathcal{N}_{c}$ the random number of computations of terms $f_{n}$ used in order to simulate just one random variable $X$. In fact $\mathcal{N}_{c}$ also corresponds to the number of random variable $S$ generated before the algorithm halts. We also define $\mathcal{N}_{l}$ the local counter depending on $Y$ which represents the number of terms $f_{n}$ used till the decision of acceptance or rejection of the variable $Y$ can be taken. Theorem IV.5.2 in Devroye [10] emphasizes the following upper-bound:

[TABLE]

where $R_{n}(y)$ is the reminder of the series expansion defined in the general framework (1.8). Let us also notice that the starting idea behind this bound is quite classical: it suffices to use the classical expansion: $\mathbb{E}[\mathcal{N}_{l}|Y]\leq\sum_{n=0}^{\infty}\mathbb{P}(\mathcal{N}_{l}>n|Y)$. In the study developed in the previous section, we obtained precise bounds for the reminder terms only for $n\geq 1$. So we shall modify the bound isolating the term $n=0$ - instead of (1.17) - which plays a crucial role for the description of the algorithm efficiency:

[TABLE]

Since the uniform bound $|R_{n}(y)|\leq r_{n}$ holds for $n\geq 1$, and since the average number of iterations is $\kappa/I(f)$, Wald’s inequality immediately implies that the total number $\mathcal{N}_{c}$ satisfies for general positive integrable functions $f$:

[TABLE]

Let us describe the consequences of this general statement to our particular algorithms.

For small values of the time variable $t$.

For small $t$, it is useful to use the series expansion (1.5) and the associated reminder bounds presented in Table 1.

Proposition 1.1**.**

For $t>0$ and $x\in[-1,1]$, we observe the following bound for the number of computations in the convergent series algorithm associated to the expansion (1.5):

[TABLE]

where $\vartheta_{2}$ stands for the Jacobi theta function given by

[TABLE]

see [20] for instance. In particular the r.h.s of the previous inequality tends to $3$ as $t$ tends to [math].

Proof.

The arguments are based on the upper-bound (1.18). Let us recall that

[TABLE]

Using the bound

[TABLE]

we obtain

[TABLE]

In order to conclude, it suffices to replace $\kappa$ by (1.12) and to notice that $I(f)=\mathbb{P}_{x}(\tau>t)$. ∎

For large values of the time variable $t$.

For large values of $t$ it is more convenient to use the series expansion (1.6). A similar approach to Proposition 1.1 easily leads to the upper bound:

Proposition 1.2**.**

Let $t\geq 0$ and $x\in[-1,1]$. The number of computations $\mathcal{N}_{c}$ for the convergent series algorithm associated to the expansion (1.6) satisfies

[TABLE]

where the constant $\kappa(t,x)$ associated to (1.6) is described in Section 1.3 and $\vartheta_{3}$ stands for the Jacobi theta function given by

[TABLE]

Since the proofs of Proposition 1.2 and Proposition 1.1 are similar, we let the details of the proof to the reader.

2 Distribution of the Brownian first exit time

In this section, we focus our attention to the distribution of the first Brownian exit $(\tau,B^{x}_{\tau})$, where $B^{x}$ stands for the one-dimensional Brownian motion starting in $x$ and $\tau$ the exit time of the interval $[-1,1]$ as defined in (1.2).

2.1 Distribution of the first exit time $\tau$

Since $p(t,x,y)$ satisfies two different series expansions (1.5) and (1.6), we can also obtain two series for the cumulative distribution of the exit time $\tau$:

[TABLE]

Therefore (1.5) becomes

[TABLE]

We deduce easily the expression of the pdf:

[TABLE]

Of course the particular case $x=0$ ensures simplifications: the symmetry of the function $\phi$ implies $\phi(-(2n+1)/\sqrt{t})=\phi(-(2\times(-n-1)+1)/\sqrt{t})$ and therefore

[TABLE]

Such an expression for the pdf of $\tau$ is of prime interest for simulation purposes. Let us note that the reminder of the series is small for small values of $t$. Indeed the sequence $(R_{1}(2n+1,t))_{n\geq 1}$ is a decreasing sequence under the assumption $t\leq 9$. For large values $t$, it is more convenient to consider (1.6) which leads to

[TABLE]

(see also (3.3) in Milstein and Tretyakov [18]). We deduce the expression

[TABLE]

Here also the case $x=0$ plays a crucial role due to its simple expression:

[TABLE]

This expression is of prime interest as soon as the variable $t$ is large. Indeed the series becomes alternating with a decreasing sequence $(R_{2}(2n+1))_{n\geq 1}$ as soon as $t\geq 4/(9\pi^{2})$.

2.2 Algorithms and numerics

The aim is to describe an algorithm which permits to simulate both the Brownian exit time $\tau$ and the exit position $B^{x}_{\tau}$ of the interval $[-1,1]$. In the previous section, we have pointed out two different series (2.1) and (2.3) corresponding to the pdf of $\tau$. We could therefore apply the classical convergent series method for simulation purposes. But here, we prefer to introduce another method based on an iterative procedure, the advantage of our approach is to deal with alternative series rather than general convergent series and therefore the reminder of the series is easier to bound.

First of all, we shall focus our attention to the case $x=0$. In this case the interval is symmetric and the expression of the pdf of $\tau$ is simplified, see (2.2) and (2.4). Moreover these series are alternating for suitable conditions on the time variable $t$. So we can apply the alternating series method for the simulation of the exit time $\tau$, the exit location being just uniformly distributed in $\{-1,1\}$. The alternating series method is an acceptance/rejection method. We need therefore a proposal distribution and an acceptance procedure. For the proposal distribution, let us first fix $t_{\rm e}>0$ (threshold for the exit distribution). Let us consider a standard gaussian random variable $G$. We define a new variable $Y$ as follows: if $1/G^{2}\leq t_{\rm e}$ then $Y=1/G^{2}$ otherwise $Y=t_{\rm e}-\frac{8}{\pi^{2}}\log(U)$ where $U$ is uniformly distributed on $[0,1]$ and independent of $G$. The density function of $Y$ which corresponds to the proposal distribution of the algorithm satisfies:

[TABLE]

where $\kappa^{-1}=\pi\ {\rm erf}(\sqrt{1/(2t_{\rm e})})\ e^{\pi^{2}t_{\rm e}/8}/2$. Using (2.2) and (2.4), we deduce the following link between the proposal distribution and the target distribution:

[TABLE]

Let us present now the acceptance/rejection method applied to this particular situation which permits to simulate the exit time $\tau$ in the symmetric case.

BROWNIAN EXIT TIME FOR SYMM. INTERVALS (with parameter $t_{\rm e}$)

BROWNIAN_EXIT_SYMMETRIC

First initialization: $\mathcal{N}_{s}=0$.

Step 0: Second initialization. *$n=0$, ${\rm Test}=0$, $L_{0}=0$, $U_{0}=1$.

Step 1. Generate a random variable $Y$ with pdf $\hat{h}$ given by (2.5). If $Y\leq t_{\rm e}$ then set $i=1$ and $C=1$ else set $i=2$ and $C=\kappa$.

Step 2. Generate a random variable $V$ uniformly distributed on $[0,1]$.

Step 3. While $(V<CU_{n})\ \&\ ({\rm Test}=0)$ do:*

• •

$n\leftarrow n+1$* and $\mathcal{N}_{s}\leftarrow\mathcal{N}_{s}+1$,*

• •

$L_{n}=U_{n-1}-\frac{R_{i}(4n-1,Y)}{R_{i}(1,Y)}$, $U_{n}=L_{n}+\frac{R_{i}(4n+1,Y)}{R_{i}(1,Y)}$

• •

${\rm Test}=1_{\{V\leq CL_{n}\}}$**

*Step 4. If ${\rm Test}=1$ then $Z=Y$ otherwise go to Step 0.

Outcome: the random variable $Z$ with density $p_{\tau}$ and the number of incrementations needed $\mathcal{N}_{s}$.

---

The outcome variable distribution $p_{\tau}$ corresponds to the target one as soon as the sequences $(R_{i}(2n+1,t)/R_{i}(1,t))_{n\geq 1}$ appearing in (2.8) are decreasing. This is an easy adaptation of the classical alternating series method (proof left to the reader). Such a property leads to the condition:

[TABLE]

Let us just note that the probability of acceptance $A=\{{\rm Accept}\ Y\}$ in the acceptance/rejection algorithm satisfies:

[TABLE]

We deduce that the number of random variables $Y$ simulated in order to obtain $Z$ is geometrically distributed with average $2$ and does not depend on the choice of $t_{\rm e}$. Nevertheless the parameter $t_{\rm e}$ has an influence on the efficiency of the algorithm as illustrated in Fig. 5. This figure represents the averaged value of the number of iterations of Step 3 needed by the algorithm in order to simulate one r.v. with the p.d.f. $p_{\tau}$. This figure suggest to choose a parameter $t_{\rm e}$ of the order of $1/2$. For this parameter, we obtain:

[TABLE]

Proposition 2.1**.**

Let us note that $\mathcal{N}_{s}$ is the random number of iterations of Step 3 used in the previous algorithm called BROWNIAN_EXIT_SYMMETRIC in order to simulate just one random variable $Z$ with density $p_{\tau}$. Then

[TABLE]

where $x\mapsto{\rm erfc}(x)$ is the complementary error function and $t_{\rm e}$ is the parameter appearing in the algorithm and satisfying (2.9).

In particular, for $t_{\rm e}=1/2$, the value of the upper bound is approximatively $1.027$ which emphasizes that the inequality is quite sharp in comparison with the estimated number of iterations in Figure 5 (solid line/stars). Moreover the bound is reasonably small which permits to confirm the efficiency of the proposed algorithm.

Proof.

The arguments are quite similar to those developed in Theorem 5.1 in Devroye [10] for the efficiency of the alternating series method. Let us denote by $\mathcal{N}_{s}^{\rm loop}$ the number of steps of type 3 used in order to go from Step 1 to Step 4, in other words: the number of increments $n\leftarrow n+1$ during one loop. By Wald’s equation and since the number of Steps 0 in this algorithm is geometrically distributed with parameter $1/2$ (see (2.2)), we get

[TABLE]

Let us first compute the probability of the following event $\{\mathcal{N}_{s}^{\rm loop}>0\}$. This event corresponds to $\{V<C\}$. In the following, we shall just recall that the constant $C$ so as $L_{1}$, $U_{1}$… depend on $Y$, that’s why we use from now on the notation $C=C_{Y}$. We obtain:

[TABLE]

since $0\leq C_{Y}\leq 1$ a.s. Let us now consider the event $\{\mathcal{N}_{s}^{\rm loop}>1\}$. This event corresponds to $\{CL_{1}<V<CU_{1}\}$ and therefore

[TABLE]

The definition of $L_{1}$ and $U_{1}$ leads to

[TABLE]

where $i=1$ for $Y\leq t_{\rm e}$ and $i=2$ otherwise. The same arguments lead to

[TABLE]

We deduce

[TABLE]

Using the distribution of the random variable $Y$ with density $\hat{h}$ in (2.5), we obtain

[TABLE]

Combining the definition of $R_{1}$ in (2.2) and the change of variable $w=\frac{y}{(4k+1)^{2}}$ leads to

[TABLE]

By Fubini’s theorem, we obtain

[TABLE]

By the change of variable $z=1/w$, we get

[TABLE]

Moreover, using the change of variable $w=r^{-2}/2$, we have

[TABLE]

Combining (2.13) and (2.14) we obtain the following upper-bound:

[TABLE]

Let us note that $\mathcal{A}_{1}(t_{\rm e})$ becomes small as $t_{\rm e}$ becomes small. Let us now focus our attention to $\mathcal{A}_{2}$ defined in (2.2).

[TABLE]

Combining (2.15) and (2.2) permits to obtain the announced bound. ∎

Let us now come to the exit simulation in the asymmetric case. Since we know how to handle with symmetric intervals, we shall use this first algorithm BROWNIAN_EXIT_SYMMETRIC in an iteration procedure in order to solve the asymmetric case. The main idea consists in the following :

• •

starting in $x$, we consider the largest interval centered in $x$ of the type $[x-\delta,x+\delta]$ and included in $[-1,1]$.

• •

We simulate the exit time of this interval denoted by $T_{1}$ and the exit position will be uniformly distributed in $\{x-\delta,x+\delta\}$.

• •

If $B^{x}_{T_{1}}\in\{-1,1\}$ then we set $T=T_{1}$ and $X=B^{x}_{T_{1}}$ else we start a new simulation for the exit time and position of a Brownian motion $\tilde{B}^{x}$ starting in $x=B_{T_{1}}$ from the largest symmetric interval centered in $x$ and included in $[-1,1]$. This exit time is denoted by $T_{2}$.

• •

As above, if $\tilde{B}^{x}_{T_{2}}\in\{-1,1\}$ then we set $X=\tilde{B}^{x}_{T_{2}}$ and $T=T_{1}+T_{2}$ else we start again with a new initial position in the interval $[-1,1]$.

BROWNIAN EXIT TIME FOR ASYMMETRIC INTERVALS $[a,b]$

BROWNIAN_EXIT_ASYMM

Input: $x$ (initial value of the Brownian paths) and $[a,b]$.

**Step 0: Initialization.

*$X=x$, $T=0$, ${\rm test}=0$, $L=a$, $U=b$ and $\mathcal{N}_{\rm as}=0$.

While (test=0) **do:

** Step 1. Set $D=\min(X-L,U-X)$. Generate a random variable $\tau$ (and the number of iterations $\mathcal{N}_{s}$) using the algorithm BROWNIAN_EXIT_SYMMETRIC and define $S=D^{2}\tau$. Set $T\leftarrow T+S$ and $\mathcal{N}_{\rm as}\leftarrow\mathcal{N}_{\rm as}+\mathcal{N}_{\rm s}$.

Step 2. Generate a random variable $V$ uniformly distributed on $\{-D,D\}$. Set $X\leftarrow X+V$

Step 3. If $X\in\{a,b\}$ then ${\rm test}=1$ else set $L\leftarrow X-L$ and $U\leftarrow U-X$.

**End While

** Outcome: the exit time $T$ and the exit location $X$ of the interval $[a,b]$ and the total number of iterations $\mathcal{N}_{\rm as}$.

---

The number of iterations is stochastically upper-bounded by a geometrically distribution with parameter $1/2$ since we only deal with symmetric intervals. Moreover this random number is independent of the generation cost of any Brownian symmetric exit time and position. Let us finally note that $(T,X)$ and $(\tau,B^{x}_{\tau})$ are obviously identically distributed.

This algorithm can be illustrated by Fig. 6 for the standard Brownian exit time of the asymmetric interval $[-1.5,2]$. Of course the algorithm is not restricted to the standard Brownian case, it takes into account any initial position $x$ belonging to the interval $[a,b]$.

Corollary 2.2**.**

The algorithm BROWNIAN_EXIT_ASYMM involves a random number of calls to BROWNIAN_EXIT_SYMMETRIC whose efficiency is characterized by the number of iterations $\mathcal{N}_{s}$. That’s why the total number of iterations $\mathcal{N}_{\rm as}$ is directly linked to the efficiency of BROWNIAN_EXIT_ASYMM. We have

[TABLE]

where $x\mapsto{\rm erfc}(x)$ is the complementary error function and $t_{\rm e}$ is the parameter appearing in BROWNIAN_EXIT_SYMMETRIC.

The statement is a direct consequence of the geometrical distributed upper-bound of the number of calls to the symmetric case algorithm on one hand (Wald’s identity therefore leads to $\mathbb{E}[\mathcal{N}_{\rm as}]\leq 2\mathbb{E}[\mathcal{N}_{s}]$) and of Proposition 2.1 on the other hand.

3 First exit time for one-dimensional diffusion

This section is concerned with the exit problem for a one-dimensional diffusion. Let us first consider $(X_{t},\ t\geq 0)$ the solution of the following stochastic differential equation:

[TABLE]

where $(B_{t},\ t\geq 0)$ stands for the standard one-dimensional Brownian motion. As already introduced in (0.1), we denote by $\tau_{a,b}(X)$ the first exit time of the interval $[a,b]$. It suffices to assume the existence of a unique weak solution to the equation, see for instance [16] for the corresponding conditions. In order to simplify the presentation of all algorithms, we restrict our study to the constant diffusion case: $\sigma(x)\equiv 1$. Using the classical Lamperti transform, we observe that this restriction is not sharp at all. That is why, from now on, $X$ stands for the unique solution of

[TABLE]

Let us note that, for the particular Brownian case: $\mu(x)\equiv 0$, the first exit time has already been presented in Section 2. The aim is to use the results developed in the previous sections when considering the general diffusion case and the important tool for such a strategy is Girsanov’s formula.

3.1 First Exit Time Algorithm (DET)

Let us now consider the exact simulation algorithm which permits to handle with the diffusion exit problem.

DIFFUSION EXIT TIME (DET)

Parameter: $\gamma_{0}$, input functions $\gamma(\cdot)$ and $\beta(\cdot)$

First initialization: $\mathcal{N}_{\rm tot}=0$.

Step 0: Initialization. *$Z=x$, $T=0$, ${\rm test}=0$. Here $x$ stands for the initial value of the diffusion.

While (test=0) **do:

** Step 1. Generate an expon. distr. random variable $E$ with parameter $\gamma_{0}$ and $U$ and $V$ two random variables uniformly distributed on $[0,1]$, the variables $E$, $U$ and $V$ being independent.

Step 2. Simulate the Brownian exit time and location*

[TABLE]

*and set $\mathcal{N}_{\rm tot}\leftarrow\mathcal{N}_{\rm tot}+\mathcal{N}_{\rm as}$.

Step 3. If $S<E$ then*

• •

if* $U\leq\beta(Y)$ then set ${\rm test}=1$, $Z\leftarrow Y$ and *$T\leftarrow T+S$

*                  else go to Step 0 end if,*

else**

• •

simulate the Brownian location at time $E$ given that the exit time is larger than $E$: $(Y_{c},\mathcal{N}_{c})=\mbox{\scriptsize CONDITIONAL\_DISTR\,(Z,[a,b],E)}$, and set $\mathcal{N}_{\rm tot}\leftarrow\mathcal{N}_{\rm tot}+\mathcal{N}_{c}$.

• •

if $\gamma_{0}V\leq\gamma(Y_{c})$ then go to Step 0 else set $Z\leftarrow Y_{c}$ and $T\leftarrow T+E$.

**end if.

*End While

Outcome: the exit location $Z$, the exit time $T$ of the interval $[a,b]$ and the efficiency index $\mathcal{N}_{\rm tot}$.

---

Let $\rho\geq 0$ such that $\gamma(x):=\frac{\mu^{2}(x)+\mu^{\prime}(x)}{2}+\rho$ is a non negative function on the interval $[a,b]$. We introduce $\gamma_{+}=\sup_{x\in[a,b]}\gamma(x)$ and define $\Delta=\int_{a}^{b}\mu(y)dy$ and

[TABLE]

Let us just notice that $0\leq\beta(a)\leq 1$ and $0\leq\beta(b)\leq 1$.

Theorem 3.1**.**

The outcome $(Z,T)$ of the Algorithm (DET) with parameter $\gamma_{0}=\gamma_{+}$ and input functions $\gamma(\cdot)$ and $\beta(\cdot)$ satisfies: for any non-negative measurable functions $f$ and $g$,

[TABLE]

where $X$ stands for the diffusion defined by (3.2). In particular, if $\mu^{\prime}+\mu^{2}$ is a non-negative function on the interval $[a,b]$ (we set $\rho=0$), then the outcome $(Z,T)$ has the same distribution as $(X_{\tau_{a,b}(X)},\tau_{a,b}(X))$.

Proof.

The proof of Theorem 3.1 is organized as follows: first we introduce a stochastic theoretical model which is based on the one-dimensional Brownian motion and on an independent Poisson process. This model is closely related to the exact simulation introduced by Beskos and Roberts in [5]. Secondly we consider in details the outcome of Algorithm (DET) and point out the strong link between the outcome and random variables derived from the theoretical model of the first step. Finally we shall use the Girsanov transformation in order to conclude the proof.

Step 1. A theoretical model. Let us consider for any non-negative functions $f$ and $g$:

[TABLE]

where $\tau_{B}$ means $\tau_{a,b}(B)$ for notational simplicity and $B$ represents a Brownian motion starting in $x\in[a,b]$ and $\gamma$ a non negative continuous function. Let us now introduce a Poisson point process $N$ on the space $\mathbb{R}_{+}\times[0,\gamma_{+}]$ with Lebesgue intensity measure and independent of the Brownian paths. If we reorganize the points of the Poisson process with respect to the abscissa, we obtain a sequence of points defined by $(\xi_{n},U_{n})_{n\geq 1}$ where $(U_{n})$ are independent uniformly distributed r.v. on $[0,\gamma_{+}]$ and $\xi_{n}:=\sum_{j=1}^{n}e_{j}$ with $(e_{j})_{j\geq 1}$ a sequence of independent exponentially distributed r.v. with parameter $\gamma_{+}$. Both sequences $(U_{n})_{n}$ and $(e_{n})_{n}$ are independent. Therefore, for any subset $A$ of $\mathbb{R}_{+}\times[0,\gamma_{+}]$, we get

[TABLE]

where $\lambda(\cdot)$ is the Lebesgue measure. In particular, if $A$ is defined as follows:

[TABLE]

the independence of both the Brownian motion and the Poisson process implies

[TABLE]

where

[TABLE]

Let us first observe that $\{\tau_{B}\leq\xi_{1}\}\subset\{N(A)=0\}$ which implies that

[TABLE]

For $n\geq 1$, due to the Markov property of the Brownian path and since the Poisson process is independent of the Brownian motion, we have

[TABLE]

Since $\xi_{1}$ is exponentially distributed and independent of the Brownian motion, we get

[TABLE]

where $\Gamma(y):=\mathbb{P}(U_{1}>\gamma(y))$ and $p(t,x)=\mathbb{P}_{x}(\tau_{B}>t)$. The probability measure $\mathcal{M}(x,dt,dy)$ represents the distribution of the couple $(T,Y)$ where $T$ is an exponentially distributed r.v. with parameter $\gamma_{+}$ and $Y$ represents the distribution of $B_{T}$ the value of the Brownian motion at time $T$ starting in $x$ and given $\{\tau_{B}>T\}$, $T$ and $(B_{t})$ being independent. The recurrence relations (3.1) and (3.1) are satisfied for any $n\geq 1$.

Step 2. Relation between Algorithm (DET) and the Brownian-Poisson theoretical model presented in Step 1.

If we denote by $\mathcal{N}_{0}$ the number of Step 0 used during the procedure which leads to the computation of the outcome $(Z,T)$ and $\mathcal{N}_{1}$ the number of exponentially distributed random variables generated (Step 1), then we obviously obtain

[TABLE]

using (3.1).

Let us now consider that $\mathcal{N}_{1}>1$. On the event $\{\mathcal{N}_{0}=1\}\cap\{\mathcal{N}_{1}>1\}$, we observe that

• •

the first exponentially distr. r.v. $E$ satisfies $E<S$ where $S$ is given by Step 2 and corresponds to the Brownian exit time of the interval $[a,b]$.

• •

$Y_{c}$ given that $E<S$ satisfies $V>\gamma(Y_{c})$, where $V$ is uniformly distributed and $Y_{c}$ becomes the new starting point for the next use of the exit problem.

In other words, if we denote by $\mathcal{I}_{n}:=\mathbb{E}_{x}[f(T)g(Z)1_{\{\mathcal{N}_{0}=1,\ \mathcal{N}_{1}=n+1\}}]$ then (3.6) leads to $\mathcal{I}_{0}(x,f,g)=I_{0}(x,f,g\cdot\beta)$. Moreover by (3.1)

[TABLE]

and using the same arguments, we prove easily that $\mathcal{I}_{n}$ satisfies the recurrence relations (3.1) and (3.1). We deduce that

[TABLE]

Therefore

[TABLE]

Since the Algorithm (DET) is an acceptance/rejection algorithm, (3.10) leads to

[TABLE]

By (3.5)

[TABLE]

Step 3. The Girsanov transformation. In this last part of the proof, the aim is to link the distribution of $(T,Z)$ described in (3.11) to the distribution of $(X_{\tau_{a,b}(X)},\tau_{a,b}(X))$ where $X$ is the diffusion defined by (3.2). Using the definition of the function $\beta$ and Itô’s formula, we obtain that

[TABLE]

is an exponential martingale. Moreover the stopped martingale $(M_{t\wedge\tau_{B}})_{t\geq 0}$ is bounded. The stopping theorem therefore implies $\mathbb{E}_{x}[M_{\tau_{B}}]=1$ which means that

[TABLE]

The expression (3.11) and Girsanov’s transformation permits to obtain

[TABLE]

where $X$ stands for the diffusion defined by (3.2). In particular, if $\rho=0$ that is $\mu^{2}+\mu^{\prime}$ is a non-negative function, we get

[TABLE]

∎

3.2 Efficiency of the algorithm

Let us now focus our attention to the analysis of Algorithm (DET). Of course since the algorithm permits to simulate exactly the random variables desired, there is no error terms to deal with, it suffices therefore to describe the time needed by the algorithm. We introduced in Algorithm (DET) the random variable $\mathcal{N}_{\rm tot}$ which permits to have a precise idea of the efficiency. Let us just add the information concerning the starting position of the Brownian motion $B_{0}=x$ with the following notation $\mathcal{N}_{\rm tot}=\mathcal{N}_{\rm tot}^{x}$. Let us also note that the efficiency of Algorithm (DET) in particular depends on two different parameters: $t_{\rm e}$ which appears in the use of the algorithm BROWNIAN_EXIT_SYM  and therefore also in BROWNIAN_EXIT_ASYMM  and $t_{\rm c}\in\mathbb{R}_{+}$ which appears in CONDITIONAL_DISTR (if $t\leq t_{\rm c}$ we consider the algorithm associated to the small values of $t$ and for $t>t_{\rm c}$ the algorithm associated to the large values).

Using informations concerning the cost of each part of the algorithm, namely Proposition 1.1, Proposition 1.2 and Proposition 2.1, we obtain a bound for $\mathbb{E}[\mathcal{N}_{\rm tot}^{x}]$.

Theorem 3.2**.**

The random variable $\mathcal{N}_{\rm tot}^{x}$ which is one of the outcomes of Algorithm (DET) and represents its cost satisfies the following bound: there exists a constant $C(t_{\rm c},t_{\rm e})>0$ independent of both the interval $[a,b]$ and the starting position $x$, such that

[TABLE]

We recall that $\tau_{a,b}$ stands for the first exit time of the interval $[a,b]$. The parameters $\rho$ and $\gamma_{+}$ are defined in the introduction of Theorem 3.1.

Each term appearing in the denominator, that is $\mathbb{E}_{x}[e^{-\rho\tau_{a,b}(X)}]$ on one hand and $\mathbb{E}_{(a+b)/2}[e^{-\gamma_{+}\tau_{a,b}(B)}]$ on the other hand, tends to [math] when the interval size $b-a$ tends to infinity. It is therefore important to choose $\rho$ and $\gamma_{+}$ as small as possible in order to obtain a sharper bound. The exit problem of the Brownian motion can be precisely described using classical results on Laplace transforms and differential equations, see for instance [9] or [6]. We obtain

[TABLE]

Let us also note that $\mathbb{E}_{x}[e^{-\rho\tau_{a,b}(X)}]$ can be linked to the two linearly independent solutions of the differential equation (see, for instance [9])

[TABLE]

which permits to describe the asymptotic behaviour as $b-a$ tends to infinity. For the link between the Laplace transform and the speed measure, see for instance [2].

Proof.

The counter $\mathcal{N}_{\rm tot}$ introduced in the algorithm can be decomposed as follows:

[TABLE]

where $\mathcal{N}_{\rm tot}^{x,k}$ represents the number of counter increases observed in-between the $k$-th and $(k+1)$-th passage through the item Step 0. We recall that we defined $\mathcal{N}_{0}$ in the proof of Theorem 3.1: it corresponds to the number of Step 0 necessary to obtain the desired outcome. Since Algorithm (DET) is an acceptance-rejection algorithm, the random variable $\mathcal{N}_{0}$ is geometrically distributed. Let us also note that $\mathcal{N}_{\rm tot}^{x,k}=0$ a.s. on the event $\{\mathcal{N}_{0}<k\}$ and conditionally to $\{\mathcal{N}_{0}\geq k\}$, $\mathcal{N}_{\rm tot}^{x,k}$ has the same distribution as $\mathcal{N}_{\rm tot}^{x,1}$. Hence

[TABLE]

The last equality is related to the third step in the proof of Theorem 3.1. Let us now describe $\mathbb{E}[\mathcal{N}_{\rm tot}^{x,1}]$. Let $(S,Y,\mathcal{N}_{\rm as})$ stands for the result of the first use of the function BROWNIAN_EXIT_ASYMM  and $(Y_{c},\mathcal{N}_{c})$ of the first use of CONDITIONAL_DISTR , we can therefore distinguish two different cases.

• •

if $S<E$ then $\mathcal{N}^{x,1}_{\rm tot}=\mathcal{N}_{\rm as}$.

• •

if $S>E$ then $\mathcal{N}^{x,1}_{\rm tot}=\mathcal{N}_{\rm as}+\mathcal{N}_{c}+\widehat{\mathcal{N}}^{Y_{c},1}_{\rm tot}1_{\{\gamma_{0}V>\gamma(Y_{c})\}}$ where $\widehat{\mathcal{N}}^{x,1}_{\rm tot}$ is an independent copy of $\mathcal{N}^{x,1}_{\rm tot}$. Such a property is essentially based on the Markov property of the Brownian paths.

We deduce

[TABLE]

• •

Let us note that $\mathbb{P}(S>E)=\mathbb{P}_{x}(\tau_{a,b}(B)>\xi)$ where $\tau_{a,b}(B)$ is the exit time of the Brownian motion from the interval $[a,b]$ and $\xi$ is exponentially distributed with parameter $\gamma_{+}$. Using the scaling property (1.1), we have

[TABLE]

where $\tau$ is the first Brownian exit time of the normalized interval $[-1,1]$ and $y=\frac{2x-a-b}{b-a}$. Since $y\mapsto\mathbb{P}_{y}(\tau>4\xi(b-a)^{-2})$ is a concave function whose derivative vanishes for $y=0$ (see the expression (2.3)), we get

[TABLE]

• •

Since $\mathcal{N}_{\rm as}$ is the cost of the function $\mbox{\scriptsize BROWNIAN\_EXIT\_ASYMM\,}(x,[a,b])$ and using the scaling property (1.1), we obtain that it is equal to the cost of $\mbox{\scriptsize BROWNIAN\_EXIT\_ASYMM\,}(y,[-1,+1])$ which satisfies due to Corollary 2.2: $\mathbb{E}[\mathcal{N}_{\rm as}]\leq C_{0}(t_{\rm e})$ with

[TABLE]

where $x\mapsto{\rm erfc}(x)$ is the complementary error function and $t_{\rm e}$ is the parameter appearing in BROWNIAN_EXIT_SYMMETRIC. This parameter $t_{\rm e}$ satisfying (2.9) can be chosen in order to minimize this average.

• •

Moreover we need some information on $\mathbb{E}_{x}[\mathcal{N}_{c}|E=t]$ where $\mathcal{N}_{c}$ is the cost of the function $\mbox{\scriptsize CONDITIONAL\_DISTR\,}(x,[a,b],t)$. Due to the scaling property (1.1), we know that the cost of this function is identical as the cost of $\mbox{\scriptsize CONDITIONAL\_DISTR\,}(y,[-1,1],4t(b-a)^{-2})$. We choose a parameter $t_{\rm c}>0$ such that $4t(b-a)^{-2}\leq t_{\rm c}$ corresponds to the algorithm for small time values and $4t(b-a)^{-2}>t_{\rm c}$ corresponds to large times. Combining Proposition 1.1 and Proposition 1.2 permits to have the following bound: for any $x\in[a,b]$,

[TABLE]

Let us introduce $C_{1}(t_{\rm c}):=\sup_{t\geq 0}\Theta(t)<\infty$, then

[TABLE]

Using these three items and (3.2), we obtain

[TABLE]

The identity (3.2) permits to conclude the proof.

∎

3.3 Modifications of Algorithm (DET)

Under the drift condition $\mu^{\prime}+\mu^{2}\geq 0$, Algorithm (DET) permits to simulate in an exact way the exit time of the interval $[a,b]$ for diffusion processes. Let us now improve this algorithm in order to deal with any one dimensional diffusion. This generalization first requires a modified algorithm with outcome $(X_{\tau_{a,b}(X)\wedge\kappa},\tau_{a,b}(X)\wedge\kappa)$ for any $\kappa>0$. The simulation of $(X_{\tau_{a,b}(X)},\tau_{a,b}(X))$ can then be obtained by iteration due to the diffusion Markov property.

Let us present Algorithm ($\kappa$-DET).

MODIFIED DIFFUSION EXIT TIME ($\kappa$-DET)

Parameters: $\rho$ and $\gamma_{0}$, input functions $\gamma(\cdot)$ and $\beta_{m}(\cdot)$

First initialization: $\mathcal{N}_{\rm tot}=0$.

Step 0: Initialization. *$K=\kappa$, $Z=x$, $T=0$, ${\rm test}=0$. Here $x$ stands for the initial value of the diffusion and $\kappa$ the time upper-bound.

While (test=0) **do:

** Step 1 Generate an expon. distr. random variable $E$ with parameter $\gamma_{0}$ and $U$, $V$ and $W$ three random variables uniformly distributed on $[0,1]$, the variables $E$, $U$, $V$ and $W$ being independent.

Step 2. Simulate the Brownian exit time and location*

[TABLE]

*and set $\mathcal{N}_{\rm tot}\leftarrow\mathcal{N}_{\rm tot}+\mathcal{N}_{\rm as}$.

Step 3. If $S=\min(K,E,S)$ then*

• •

if* $U\leq\beta_{m}(Y)$ and $W\leq e^{-\rho\,(K-S)}$ *

*                  then set ${\rm test}=1$, $Z\leftarrow Y$ and *$T\leftarrow T+K$

*                  else go to Step 0 end if,*

elseif $K=\min(K,E,S)$**

• •

simulate the conditional Brownian location at time $K$:

$(Y_{c},\mathcal{N}_{c})=\mbox{\scriptsize CONDITIONAL\_DISTR\,(Z,[a,b],K)}$* and set $\mathcal{N}_{\rm tot}\leftarrow\mathcal{N}_{\rm tot}+\mathcal{N}_{c}$.*

• •

if* $U\leq\beta_{m}(Y_{c})$ then set ${\rm test}=1$, $Z\leftarrow Y$ and *$T\leftarrow T+S$

*                  else go to Step 0 end if,*

elseif $E=\min(K,E,S)$**

• •

simulate the conditional Brownian location at time $E$:

$(Y_{c},\mathcal{N}_{c})=\mbox{\scriptsize CONDITIONAL\_DISTR\,(Z,[a,b],E)}$* and set $\mathcal{N}_{\rm tot}\leftarrow\mathcal{N}_{\rm tot}+\mathcal{N}_{c}$.*

• •

if $V\leq\gamma(Y_{c})$ then go to Step 0 else set $Z\leftarrow Y_{c}$, $T\leftarrow T+E$ and $K\leftarrow K-E$.

**end if.

*End While

Outcome: the stopping time $T$ corresponding to the minimum between the diffusion exit time of the interval $[a,b]$ and the constant time $\kappa$, the location $Z$ of the stopped diffusion at time $T$ and the efficiency index $\mathcal{N}_{\rm tot}$.***

Let $\rho\geq 0$ such that $\gamma(x):=\frac{\mu^{2}(x)+\mu^{\prime}(x)}{2}+\rho$ is a non negative function on the interval $[a,b]$. We recall that $\gamma_{+}$ and $\beta$ are defined in (3.3). We introduce

[TABLE]

Let us just notice that $0\leq\beta_{m}(x)\leq 1$ for any $x\in[a,b]$.

Theorem 3.3**.**

The outcome $(Z,T)$ of the Algorithm ($\kappa$-DET) with parameter $\rho$, $\gamma_{0}=\gamma_{+}$ and input functions $\gamma(\cdot)$ and $\beta_{m}(\cdot)$ has the same distribution as $(X_{\tau_{a,b}(X)\wedge\kappa},\tau_{a,b}(X)\wedge\kappa)$, where $X$ stands for the diffusion defined by (3.2).

Proposition 3.4**.**

The random variable $\mathcal{N}_{\rm tot}^{x}$ which is one of the outcomes of Algorithm ($\kappa$-DET) satisfies the following bound: there exists a constant $C(t_{\rm c},t_{\rm e})>0$ independent of both the interval $[a,b]$ and the starting position $x$, such that

[TABLE]

We recall that $\tau_{a,b}$ stands for the first exit time of the interval $[a,b]$. The parameters $\rho$ and $\gamma_{+}$ are defined in the introduction of Theorem 3.1.

The proof of Proposition 3.4 is just a slight modification of the proof of Theorem 3.2. The details are left to the reader (see Appendix).

Proof of Theorem 3.3.

The key arguments are similar to those used in the proof of Theorem 3.1 and the structure of the proof is the same.

*Step 1. The Poisson-Brownian model. * Let us define the following expression depending on the Brownian paths $(B_{t})_{t\geq 0}$, on the first exit time of the interval $[a,b]$ denoted by $\tau_{B}$ and on the independent Poisson process $N$:

[TABLE]

with

[TABLE]

Let us introduce the following series expansion:

[TABLE]

the definition of $I_{n}(x,f,g,\kappa)$ is similar to the definition appearing in the proof of Theorem 3.1, it suffices to replace in the definition $\tau_{B}$ by $\tau_{B}\wedge\kappa$. Therefore

[TABLE]

where $(\xi_{1},U_{1})$ stands for the coordinates of the Poisson process point with the smallest abscissa. Moreover we obtain the following step by step property:

[TABLE]

Step 2. Relation between Algorithm ($\kappa$-DET) and the Brownian-Poisson theoretical model presented in Step 1.

If we denote by $\mathcal{N}_{0}$ the number of Step 0 used during the procedure which leads to the computation of the outcome $(Z,T)$ and $\mathcal{N}_{1}$ the number of exponentially distributed random variables generated (Step 1), then we obviously obtain for $\mathcal{I}_{n}(x,f,g,\kappa):=\mathbb{E}_{x}[f(T)g(Z)1_{\{\mathcal{N}_{0}=1,\ \mathcal{N}_{1}=n+1\}}]$,

[TABLE]

Using similar arguments, we can prove that $\mathcal{I}_{n}$ satisfies the same step by step relation as(3.16). So by identification, we get $\mathcal{I}_{n}(x,f,g,\kappa)=I_{n}(x,f\times e^{-\rho(\kappa-\cdot)},g\times\beta_{m},\kappa)$ for all $n\geq 0$ and therefore

[TABLE]

Since Algorithm ($\kappa$-DET) is an acceptance rejection algorithm, we have

[TABLE]

The link between $\beta_{m}$ and $\beta$ leads to

[TABLE]

where $M_{t}$ is the exponential martingale defined in (3.1). It is actually important to note that since the time interval is bounded (upper-bounded by $\kappa$) we can use the time reversal expression $e^{-\rho(\kappa-\cdot)}$ and the ratio in the previous equation verifies a cancellation of the terms $e^{-\rho\kappa}$.

Finally it suffices therefore to use the Girsanov transformation in order to obtain the announced result: $\mathbb{E}[f(T)g(Z)]=\mathbb{E}_{x}[f(\tau_{X}\wedge\kappa)g(X_{\tau_{X}\wedge\kappa})].$ ∎

Let us now focus our attention on the exact simulation of $(X_{\tau_{a,b}},\tau_{a,b})$ for any diffusion process with regular drift term and constant diffusion coefficient. We have already seen that Algorithm (DET) permits to reach such objective however it is restricted to drift terms satisfying $\mu^{\prime}+\mu^{2}\geq 0$ on the interval $[a,b]$. If such a condition is not verified, we propose the following procedure: first we choose some time parameter $\kappa>0$, then we apply:

GENERAL DIFFUSION EXIT TIME (GDET)

parameters: $\kappa$, $a$ and $b$.

Initialization. *$Z=x$, $T=0$, $\mathcal{N}_{\rm it}=0$

While $Z\notin\{a,b\}$ do:*

• •

Simulate $(S,Y)=\kappa\mbox{{\scriptsize-DET}}(Z,\kappa)$

• •

$T\leftarrow T+S$, $Z\leftarrow Y$ and $\mathcal{N}_{\rm it}\leftarrow\mathcal{N}_{\rm it}+1$.

**End While

Outcome: the couple of random variables $(Z,T)$ and the number of iterations $\mathcal{N}_{\rm it}$.***

Due to the Markov property, it is obvious that the outcome of such an algorithm and $(X_{\tau_{a,b}(X)},\tau_{a,b}(X))$ are identically distributed. Moreover the number of iterations of GDET has the same distribution as $\lfloor\frac{\tau_{a,b}(X)}{\kappa}\rfloor+1$. Hence

[TABLE]

It is evident that the number of iterations decreases as $\kappa$ becomes large but we should be careful for a clever choice of $\kappa$ since the number of rejections in Algorithm ($\kappa$-DET) grows exponentially fast when $\kappa$ enlarges, see Proposition 3.4. A reasonable choice is therefore $\kappa\approx\rho^{-1}$.

3.4 Examples and numerics

The aim of this section is to emphasize the efficiency of Algorithm (DET) and Algorithm ($\kappa$-DET) through the analysis of two examples. The first situation concerns a diffusion whose drift term $b$ satisfies the condition $\mu^{\prime}+\mu^{2}\geq 0$ and consequently only requires the basic Algorithm (DET). The second situation concerns the Ornstein-Uhlenbeck process which plays an essential role in several applications namely in neuroscience. For both examples, we set the parameters appearing in the algorithms: $t_{\rm c}=0.7$ and $t_{\rm e}=0.5$.

3.4.1 Example with the Algorithm (DET)

We first consider a stochastic differential equation which was already presented in [13] for the simulation of the first passage time. Here the objective is clearly different since we focus our attention to the exit time and exit position of the diffusion and the algorithm is different too. Let us also note that a similar diffusion process was also introduced in [5].

We consider the following stochastic differential equation:

[TABLE]

We first observe that $\gamma(x)=(\mu^{2}(x)+\mu^{\prime}(x))/2=((2+\sin(x))^{2}+\cos(x))/2$ satisfies $0\leq\gamma\leq 5$. We deduce that we can apply Theorem 3.1 with the particular choice $\rho=0$: the outcome $(Z,T)$ of Algorithm (DET) has therefore the same distribution as $(X_{\tau_{a,b}(X)},\tau_{a,b}(X))$. The algorithm permits to obtain the histograms of respectively the exit time (Fig. 7 left) and the counter $\mathcal{N}_{\rm tot}$ illustrating the efficiency of the algorithm DET (right). These histograms use a sample of size $100\,000$ and concerns the interval $[a,b]=[-0.5,0.5]$. The average value of the counter is $8.5$ and its estimated standard deviation is $8.88$. Let us just compare the approximated results obtained by Algorithm (DET) for the interval $[a,b]=[-0.5,0.5]$ with a classical Euler method with step size $0.0001$ and $100\,000$ samples.

Algo. $\mathbb{E}[\tau_{a,b}]$ $\sigma(\tau_{a,b})$ $\mathbb{E}[\tau_{a,b}1_{\{X_{\tau_{a,b}}=a\}}]$ $\mathbb{P}(X_{\tau_{a,b}}=a)$

Euler method 0.18262 0.13796 0.18446 0.12530

(DET) 0.17927 0.13667 0.18018 0.12685

When considering larger intervals like for instance $[-1,2]$ (Fig. 8), the counter becomes large (average: 1205) and the algorithm DET is rather time consuming (C++ progaming: CPU 1,77 sec for $100\,000$ samples of the exit from the interval $[-0.5,0.5]$ and CPU 230,8 sec for the interval $[-1,2]$). Figure 9 is an illustration of the high level of rejection for the algorithm (DET) as the interval size increases.

3.4.2 The Ornstein-Uhlenbeck case.

Let us now consider the Ornstein-Uhlenbeck process given by the following SDE

[TABLE]

We shall determine the distribution of the exit time and location from the interval $[-a,a]$ with $a>0$. The drift term is defined by $\mu(x)=-\mu_{0}x$, we can therefore find $\rho$ such that

[TABLE]

It suffices to choose $\rho=\frac{\mu_{0}}{2}$. Such a choice implies

[TABLE]

Since the coefficient $\rho$ is strictly positive (for positive $\mu_{0}>0$), we cannot use the (DET)-algorithm in order to simulate the exit time. We shall therefore use the ($\kappa$-DET) algorithm in order to simulate exactly the couple $(X_{\tau_{-a,a}(X)\wedge\kappa},\tau_{-a,a}(X)\wedge\kappa)$ for any given constant time $\kappa>0$. Fig. 10 represents the distribution of $X_{\tau_{-a,a}(X)\wedge\kappa}$ and Fig. 11 illustrates the time distribution and describes the counter values $\mathcal{N}_{\rm tot}$.

In order to simulate exactly the exit time for the Ornstein-Uhlenbeck process, we use the (GDET)-algorithm. The histogram in Figure 12 emphasizes the distribution of the exit time for the particular case: $\mu_{0}=2$ and $a=1$ and the efficiency of the algorithm for such a simulation. We can easily observe that the exit time increases as the parameter $\mu_{0}$ increases, see Fig. 13.

Acknowledgment: We thank Antoine Lejay for helpful discussions about the Brownian exit times and their simulation. We particularly appreciated his working paper [17].

Appendix

Proof of Proposition 3.4.

The counter $\mathcal{N}_{\rm tot}$ depends on the parameter $\kappa$. So we set $\mathcal{N}_{\rm tot}=\mathcal{N}_{\kappa,{\rm tot}}$. This number can be decomposed as follows:

[TABLE]

where $\mathcal{N}_{\kappa,{\rm tot}}^{x,k}$ represents the number of counter increases observed in-between the $k$-th and $(k+1)$-th passage through the item Step 0. Since Algorithm ($\kappa$-DET) is an acceptance-rejection algorithm, the random variable $\mathcal{N}_{0}$ is geometrically distributed. Let us also note that $\mathcal{N}_{\kappa,{\rm tot}}^{x,k}=0$ a.s. on the event $\{\mathcal{N}_{0}<k\}$ and conditionally to $\{\mathcal{N}_{0}\geq k\}$, $\mathcal{N}_{\kappa,{\rm tot}}^{x,k}$ has the same distribution as $\mathcal{N}_{\kappa,{\rm tot}}^{x,1}$. Hence

[TABLE]

since $\mathbb{P}(\mathcal{N}_{0}=1)=I(x,e^{-\rho(\kappa-\cdot)},\beta_{m},\kappa)=e^{-\rho\kappa}\beta_{m}(x)\mathbb{E}_{x}[M_{\tau_{B}\wedge\kappa}]=e^{-\rho\kappa}\beta_{m}(x)$, the martingale $(M_{t})$ being defined in (3.1).

Let us now describe $\mathbb{E}[\mathcal{N}_{\kappa,{\rm tot}}^{x,1}]$. Let $(S,Y,\mathcal{N}_{\rm as})$ stands for the result of the first use of the function BROWNIAN_EXIT_ASYMM  and $(Y_{c}^{t},\mathcal{N}_{c}^{t})$ of the first use of CONDITIONAL_DISTR  at time $t$, we can therefore distinguish three different cases.

• •

if $S<E\wedge\kappa$ then $\mathcal{N}^{x,1}_{\kappa,{\rm tot}}=\mathcal{N}_{\rm as}$.

• •

if $\kappa<E\wedge S$ then $\mathcal{N}^{x,1}_{\kappa,{\rm tot}}=\mathcal{N}_{\rm as}+\mathcal{N}_{c}^{\kappa}$.

• •

if $E<S\wedge\kappa$ then $\mathcal{N}^{x,1}_{\kappa,{\rm tot}}=\mathcal{N}_{\rm as}+\mathcal{N}_{c}^{E}+\widehat{\mathcal{N}}^{Y_{c}^{E},1}_{\kappa-E,{\rm tot}}1_{\{\gamma_{0}V>\gamma(Y_{c}^{E})\}}$ where $\widehat{\mathcal{N}}^{x,1}_{\kappa,{\rm tot}}$ is an independent copy of $\mathcal{N}^{x,1}_{\kappa,{\rm tot}}$.

We deduce

[TABLE]

Using the same arguments as those developed in the proof of Theorem 3.2, we get the existence of two constants $C_{0}(t_{\rm e})>0$ and $C_{1}(t_{\rm c})>0$ such that $\mathbb{E}[\mathcal{N}_{\rm as}]\leq C_{0}(t_{\rm e})$ and $\mathbb{E}[\mathcal{N}_{c}^{t}1_{\{S>t\}}]\leq C_{1}(t_{\rm c})$ for all $t\geq 0$. We then deduce

[TABLE]

This result can be generalized for any $K\leq\kappa$ and therefore we obtain

[TABLE]

Let us note that $\mathbb{P}_{(a+b)/2}(S\wedge\kappa\leq E)=\mathbb{P}_{(a+b)/2}(\tau_{a,b}(B)\wedge\kappa\leq\xi)$ where $\tau_{a,b}(B)$ is the exit time of the Brownian motion from the interval $[a,b]$ and $\xi$ is exponentially distributed with parameter $\gamma_{+}$. Using the scaling property (1.1), we have

[TABLE]

Combining (3.21) (3.22) and (3.19), we obtain

[TABLE]

∎