Exact simulation of the first-passage time of diffusions

TL;DR

This paper introduces an exact rejection sampling algorithm for simulating the first-passage time of general one-dimensional diffusion processes efficiently, avoiding full path simulation and leveraging Girsanov's transformation.

Contribution

The authors develop a novel rejection sampling method for exact simulation of first-passage times in diffusions, extending beyond Brownian cases.

Findings

Algorithm achieves exact simulation without full path construction.

Efficiency is supported by theoretical analysis and numerical tests.

Method applies to a wide class of one-dimensional diffusions.

Abstract

Since diffusion processes arise in so many different fields, efficient tech-nics for the simulation of sample paths, like discretization schemes, represent crucial tools in applied probability. Such methods permit to obtain approximations of the first-passage times as a by-product. For efficiency reasons, it is particularly challenging to simulate directly this hitting time by avoiding to construct the whole paths. In the Brownian case, the distribution of the first-passage time is explicitly known and can be easily used for simulation purposes. The authors introduce a new rejection sampling algorithm which permits to perform an exact simulation of the first-passage time for general one-dimensional diffusion processes. The efficiency of the method, which is essentially based on Girsanov's transformation , is described through theoretical results and numerical examples.

Equations

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

$${\tau }_L:=\inf \{t\ge 0:X_t=L\}.$$

$$\displaystyle\begin{array}[]{l}\mbox{Step 1. Simulate}\ T_{L}\ \mbox{the first hitting time of the Brownian motion}\\ \mbox{Step 2. If the event }A\ \mbox{is true, then}\ Y=T_{L}\ \mbox{otherwise go to Step 1.}\end{array}$$

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

$$Y_t=\eta (X_t):={\int }_x^{X_t}\frac{1}{\sigma (u)}\mathrm{d}u$$

$$b_0(x)=\frac{b\circ {\eta }^{-1}(x)}{\sigma \circ {\eta }^{-1}(x)}-\frac{1}{2}{\sigma }^{\prime }\circ {\eta }^{-1}(x).$$

$$d{X}_t=b(X_t)dt+d{B}_t,t\ge 0,X_0=x.$$

$$\beta (x)={\int }_0^{x}b(y)\mathrm{d}y,p(x)={\int }_0^x{e}^{-\beta (y)}\mathrm{d}y\text{mbox}and\gamma :=\frac{b^2+b^{\prime }}{2}.$$

$$\mathbb{E}_{\mathbb{P}}[\psi(\tau_{L})1_{\{\tau_{L}<\infty\}}]=\mathbb{E}_{\mathbb{Q}}[\psi(\tau_{L})\eta(\tau_{L})]\exp\Big{\{}\beta(L)-\beta(x)\Big{\}},$$

$$\eta(t):=\mathbb{E}\Big{[}\exp-\int_{0}^{t}\gamma(L-R_{s})\mathrm{d}s\Big{|}R_{t}=L-x\Big{]}.$$

$${\mathbb{E}}_{\mathbb{P}}[\psi ({\tau }_L)1_{\{{\tau }_L<\infty \}}]$$

$$\displaystyle=\mathbb{E}_{\mathbb{Q}}\Big{[}\psi(\tau_{L})\exp\left(\beta(B_{\tau_{L}})-\beta(x)-\int_{0}^{\tau_{L}}\gamma(B_{s})\mathrm{d}s\right)\Big{]}.$$

$${\mathbb{E}}_{\mathbb{P}}[\psi ({\tau }_L)1_{\{{\tau }_L<\infty \}}]$$

$$\displaystyle=e^{\beta(L)-\beta(x)}\mathbb{E}_{\mathbb{Q}}\Big{[}\psi(\tau_{L})\eta(\tau_{L})\Big{]},$$

$$\eta(t):=\mathbb{E}_{\mathbb{Q}}\Big{[}\exp-\int_{0}^{t}\gamma(B_{s})\mathrm{d}s\Big{|}\tau_{L}=t\Big{]}.$$

$$v(x)={\int }_0^x{\int }_0^y\frac{2p^{\prime }(y)}{p^{\prime }(z)}\mathrm{d}z\mathrm{d}y,$$

$$D_{R,T}:=\Big{\{}(t,v)\in[0,T]\times\mathbb{R}_{+}:\ v\leq\gamma(L-R_{t})\Big{\}}.$$

$$f_Y(t)=\frac{1}{\Xi }\eta (t)f_T(t),$$

$$\Xi :={\int }_0^{\infty }\eta (t)f_T(t)\mathrm{d}t.$$

$$\Xi=\exp-\Big{\{}\beta(L)-\beta(x)\Big{\}}.$$

$$\mathbb{E}[\psi (Y)]$$

$$\displaystyle=\sum_{i=1}^{\infty}\mathbb{E}\Big{[}\psi(T^{(1)})1_{\overline{E}_{1}}\Big{]}\mathbb{P}(E_{1})^{i-1}=\mathbb{E}\Big{[}\psi(T^{(1)})1_{\overline{E}_{1}}\Big{]}(1-\mathbb{P}(E_{1}))^{-1}.$$

$$\mathbb{E}[\psi (Y)]$$

$$\mathbb{P}(N(D_{R,T})=0|T)$$

$$\displaystyle=\mathbb{E}\Big{[}\exp-\int_{0}^{T}\gamma(L-R_{t})\mathrm{d}t\Big{|}R_{T}=L-x,T\Big{]}=\eta(T).$$

$$\mathbb{E}[\psi (Y)]=\frac{1}{\mathbb{P}(N(D_{R,T})=0)}\mathbb{E}[\psi (T)\eta (T)].$$

$$0\le \gamma (x)\le \kappa ,\text{mbox}forallx\in ]-\infty ,L].$$

$$\beta \leftarrow \frac{T-{\mathcal{E}}_1}{T-{\mathcal{E}}_0}\beta +\sqrt{\frac{(T-{\mathcal{E}}_1)({\mathcal{E}}_1-{\mathcal{E}}_0)}{T-{\mathcal{E}}_0}}G.$$

$$\mbox{if}\quad\kappa V\leq\gamma\Big{(}L-\|\mathcal{E}_{1}(L-x)(1,0,0)/T+\beta\|\Big{)}\ \mbox{then}\ \mathcal{W}=1\ \mbox{else}\ \mathcal{W}=0,$$

$$\beta =\overleftarrow{U_2}+\frac{\overrightarrow{U_2}-\overleftarrow{U_2}}{\overrightarrow{U_1}-\overleftarrow{U_1}}(U-\overleftarrow{U_1})+\sqrt{\frac{(\overrightarrow{U_1}-U)(U-\overleftarrow{U_1})}{\overrightarrow{U_1}-\overleftarrow{U_1}}}G.$$

$$\mbox{if}\ \mathcal{E}\leq\gamma\Big{(}L-\|U(L-x)(1,0,0)/T+\beta\|\Big{)}\ \mbox{then}\ \mathcal{W}=1\ \mbox{else}\ \mathcal{W}=0.$$

Full text

Key words and phrases: first-passage time, Brownian motion, diffusion processes, Girsanov’s transformation, exact simulation, randomized algorithm.

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

Introduction

First-passage times of diffusion processes appear in several fields, i.e. economics [21], mathematical finance [23, 28], neuroscience [9, 34], physics [1], psychology [29] and reliability theory [30]. In neuroscience particular diffusion processes model the evolution of the membrane potential and the first-passage times of the considered process throught a suitable boundary describes the neuronal firing times (cf. [9, 34]). Likewise, in metrology or in quality control, diffusion processes are used to describe the error of some instrument which is required to be bounded. The first-passage time through a suitable boundary reproduces the first time the error becomes out of control [40]. In this setting, rapid detection of anomalies corresponds to recognize the optimal stopping time of a diffusion process [41]. The importance of practical applications dealing with the first passage time motivated many mathematical studies on the subject.

To introduce the corresponding mathematical problem, let us consider the stochastic process $(X_{t},\ t\geq 0)$ the solution of the SDE:

[TABLE]

where $(B_{t},\ t\geq 0)$ stands for the standard one-dimensional Brownian motion. We denote by $\tau_{L}$ the first time the diffusion reaches the level $L$:

[TABLE]

The study of the hitting times $\tau_{L}$ and of its approximations for general diffusion processes is an active area of research.

Explicit expressions or even asymptotic solutions to boundary crossing problems are only available in rare situations which moreover do not correspond to interesting instances for applications. Consequently studies on this subject have been directed to the development of alternative methods like numerical or simulative ones.

As far as numerical methods are concerned, they are efficient but need specific algorithms for each function of the hitting time $\tau_{L}$. Most of the studies focus on the probability density function $p$ defined by $p(t)dt=\mathbb{P}(\tau_{L}\in dt)$ which satisfies Volterra-type integral equations. It suffices therefore to approximate certain integral. Examples of studies in this direction are Durbin [11, 12], Ferebee [13], Ricciardi et al. [33]; Giorno et al. [14], Sacerdote and Tomassetti [36]. The numerical approach proposed in [8] seems to be particularly efficient. This type of approach is also used in [3, 35] to deal with first passage times of a two-dimensional correlated diffusion process. For particular families of diffusion which can be constructed using simple transformations of a Brownian motion, Pötzelberger and Wang [37, 32, 38] proposed a different approach based on the explicit Brownian crossing probabilities.

The second important course of action uses simulations to determine the hitting times $\tau_{L}$. These methods often require a strong computational effort and are affected by a non negligible approximation error. Since stopping times for time discrete Markov processes can be quite easily simulated, a natural way to deal with diffusion processes is to introduce a discretization scheme for the corresponding stochastic differential equation. However, only an upper bound of the stopping time can be determined through this approximation. It is therefore important to improve the algorithm. In [7] and [19] a shift of the boundary is proposed to improve the approximation, while in [15, 16, 17] the proposed method compute on each small interval of the time discretization, the probability for the process to cross the boundary in this time window. Such a method needs precise asymptotics of hitting probabilities for pinned diffusions. Gobet [18] exteded the study to a multidimensional diffusion and described the error for both the discrete and the continuous Euler schemes. Precise asymptotics for general pinned diffusions are pointed out by Baldi and Caramellino [2] permitting to deal with a larger class of diffusions.

Let us just notice that besides the simulated approximation methods there also exist simulated exact methods dealing with the sample of diffusion paths on a fixed time interval, first introduced by Beskos & Roberts [6]. Several modifications of this algorithm have been proposed [4, 5, 25]. However, the two main approaches are not exclusive: recently Herrmann and Tanré [20] proposed a new simulation method which has no link with the two approaches just described. It is based on an iterative algorithm and leads to an efficient approximation of Brownian hitting times for curved boundaries, and, by extension, can deal with diffusion processes expressed as functionals of a Brownian motion.

This paper is a contribution in the framework of the simulated exact methods. Here we present an efficient method for simulating directly the first-passage time $\tau_{L}$ without building the whole paths of the diffusion process and without any approximation error. The method can be summarized as follows: we need an algorithm producing a random variable $Y$ such that $Y$ and $\tau_{L}$ are identically distributed. In particular, this approach permits to reduce the approximation error of $\mathbb{E}[\psi(\tau_{L})]$, where $\psi:\mathbb{R}\to\mathbb{R}$ is any measurable and bounded function, to the only Monte-Carlo one. Our approach rests on two important facts: it is easy to simulate the first-passage time of the standard Brownian motion on one hand and Girsanov’s transformation allows to link the Brownian paths to any diffusion paths on the other hand. Combining these features, we highlight the following acceptance-rejection algorithm (A0):

[TABLE]

The first challenge is to find an appropriate event $A$. Subsequently we present alternatives in order to reduce the averaged number of steps of this acceptance/rejection algorithm. Our main idea originates in two different mathematical studies; let us now focus our attention on them.

• •

Ichiba & Kardaras [22] proposed an efficient estimation of the first-passage density using on one hand Girsanov’s transformation and on the second hand a Monte Carlo procedure in order to estimate the Radon-Nikodym derivative which involves a Bessel bridge of dimension $\delta=3$. This interesting approach permits to deal with a large family of diffusion processes but introduces a Monte-Carlo approximation error.

• •

The second pillar is the procedure of exact simulation introduced by Beskos & Roberts [6]. Using a judicious rejection sampling algorithm, they simulate in a surprisingly simple way a diffusion process on any finite time interval. They propose also to use their algorithm in order to simulate extremes and hitting times by building a skeleton of the whole trajectory.

The algorithm proposed here, based on these two pillars, outperform [6] when the objective is the simulation of $\tau_{L}$ instead of the whole diffusion paths.

The material is organized as follows: after reminding the famous Girsanov transformation, we present in details the rejection sampling algorithm in Section 1 and prove that the outcome $Y$ has the expected distribution. In Section 2, we focus our attention on the efficiency of the procedure: the number of steps is described and upper-bounds are computed. Unfortunately our procedure can not be applied to any diffusion processes since it requests different technical assumptions. We discuss in Section 3 and Section 4 two modifications of the main algorithm in order to weaken these assumptions. Finally, for completeness, numerical examples permit to appreciate and illustrate the efficiency of the method (Section 5).

1 Rejection sampling method for first-passage time of diffusions

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. We denote by $\tau_{L}$ the first time the diffusion reaches the level $L$ as defined in (0.1). It suffices to assume the existence of a unique weak solution to the equation, see for instance [26] 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 Lamperti transform, we observe that this restriction is not sharp at all. Indeed if $X$ is solution to (1.1) then by fixing $x\in\mathbb{R}$, we define

[TABLE]

which satisfies therefore the equation: $dY_{t}=b_{0}(Y_{t})\,dt+dB_{t}$ with

[TABLE]

That is why, from now on, $X$ stands for the unique solution of

[TABLE]

and the level to reach satisfies $L>x$. Let us note that, for the particular Brownian case: $b(x)\equiv 0$, the first-passage time is extremely simple to simulate since $\tau_{L}$ and $(L-x)^{2}/G^{2}$ have the same distribution where $G$ is a standard gaussian random variable. The aim is to use this simple form when considering the general diffusion case and the important tool for such a strategy is Girsanov’s formula. We first recall it before introducing the main rejection algorithm.

1.1 Girsanov’s transformation

Let us first introduce $3$ important expressions related to the drift term $b$ :

[TABLE]

Assumption 1.1**.**

The drift term $b\in\mathcal{C}^{1}(]-\infty,L])$.

Grisanov’s transformation permits to link the diffusion process $X$ solution of (1.2) to the Brownian motion (a slight modification of Proposition 2.1 in [22]).

Proposition 1.1**.**

Under Assumption 1.1, for any bounded measurable function $\psi:\mathbb{R}\to\mathbb{R}$, we obtain

[TABLE]

where $\mathbb{P}$ (resp. $\mathbb{Q}$) corresponds to the distribution of the diffusion $X$ (resp. the Brownian motion $B$) and

[TABLE]

Here $(R_{t},\,t\geq 0)$ stands for a $3$-dimensional Bessel process with $R_{0}=0$.

Proof.

Let us denote $\zeta$ the explosion time of the diffusion $X$ defined by (1.2). The Feller test for explosion (see for instance Proposition 5.22 in [26]) claims that Assumption 1.1 leads to $\tau_{L}<\zeta$, $\mathbb{P}$-almost surely. By applying successively Girsanov’s transformation (Proposition 1.7.5.4 in [24]) and Itô’s formula, we obtain

[TABLE]

Under $\mathbb{Q}$, $(B_{t})$ corresponds to a one-dimensional Brownian motion with starting point $x$ and $\tau_{L}=\inf\{t\geq 0:\ B_{t}=L\}<\infty$ a.s. (recurrence of the Brownian paths). Hence $B_{\tau_{L}}=L$ and consequently

[TABLE]

where

[TABLE]

In order to conclude it suffices to note that, given $\{\tau_{L}=t\}$, $R_{s}:=L-B_{t-s}$ is a Bessel bridge [39] for $s\in[0,t]$. ∎

1.2 Description of the algorithm

The rejection sampling algorithm for the simulation of $\tau_{L}$ is based on Proposition 1.1. We need to introduce the following

Assumption 1.2**.**

The function $\gamma$ defined by (1.3) is non-negative and the first-passage time satisfies $\tau_{L}<\infty$ a.s.

These main conditions permit the function $\eta$ defined by (1.5) to belong to the interval $[0,1]$ and therefore to represent a probability of rejection in the procedure.

Remark 1.2**.**

Let us just note that $\tau_{L}<\infty$ a.s. if the drift term satisfies one of the additional conditions: either $\displaystyle\lim_{x\to-\infty}p(x)=-\infty$ (see (1.3) for the definition of $p$) or $\displaystyle\lim_{x\to-\infty}v(x)<\infty$ where

[TABLE]

according to the Feller test for explosion (Proposition 5.5.32 in [26]).

Let us now focus our attention to a first theoretical description of the acceptance-rejection algorithm.

Algorithm (A0).

Let us fix the level $L>x$.

Step 1: Simulate a non-negative random variable $T$ with p.d.f. $f_{T}$.

Step 2: Simulate a $3$-dimensional Bessel process $(R_{t})$ on the time interval $[0,T]$ with endpoint $R_{T}=L-x$. We define by $D_{R,T}$ the stochastic domain:*

[TABLE]

*This domain depends on both random elements $(R_{t})$ and $T$.

Step 3: Simulate a Poisson point process $N$ on the state space $[0,T]\times\mathbb{R}_{+}$, independent of the Bessel process, whose intensity measure is the Lebesgue one.

Step 4: If $N(D_{R,T})=0$ then set $Y=T$ otherwise go to Step 1.

Outcome: the random variable $Y$.*

We shall discuss later on, how this algorithm can be applied in practice.

Theorem 1.3**.**

If $\gamma$ is a non-negative function then the p.d.f. $f_{Y}$ of the outcome variable $Y$ satisfies

[TABLE]

where $\eta$ is given by (1.5) and $\Xi$ stands for the normalization coefficient

[TABLE]

In particular, under Assumption 1.1 & 1.2, if $T/(L-x)^{2}$ has the same distribution as $1/G^{2}$ where $G$ is a standard Gaussian r.v., then $Y$ and $\tau_{L}$, defined by (0.1) & (1.2), are identically distributed.

Let us observe that in the particular case $T\sim\frac{(L-x)^{2}}{G^{2}}$, (1.4) leads to the following identify:

[TABLE]

Proof.

The arguments presented here are quite classical for a rejection sampling method. We introduce $\mathcal{I}$ the number of iterations in Algorithm (A0). For each iteration $\mathcal{I}=i$, we denote by $T^{(i)}$, $R^{(i)}_{t}$ and $N^{(i)}$ the generated random variables or processes. Of course each iteration corresponds to independent variables. Let $\psi$ be a measurable bounded function. We get

[TABLE]

where E_{i}:=\Big{\{}N^{(i)}(D_{R^{(i)},T^{(i)}})\neq 0\Big{\}} and $\overline{E}_{i}=\Omega\setminus E_{i}$. Hence, using the independence property, we obtain

[TABLE]

From now on, we delete the index for notational simplicity. Therefore

[TABLE]

The distribution of the Poisson point process leads to:

[TABLE]

We deduce

[TABLE]

For $\psi\equiv 1$, we obtain $\mathbb{P}(N(D_{R,T})=0)=\mathbb{E}[\eta(T)]$ which implies the statement of the proposition. ∎

Algorithm (A0) seems quite difficult to achieve in practice since we do not know how to simulate a Poisson point process on the state space $[0,T]\times\mathbb{R}_{+}$ on one hand and since the domain $D_{R,T}$ depends on the whole trajectory of $(R_{t})$ on the other hand. In order to overcome the first difficulty, we reduce the domain to a bounded one introducing the following assumption:

Assumption 1.3**.**

there exists a constant $\kappa>0$ such that

[TABLE]

It suffices therefore to simulate the Poisson point process on the reduced space $[0,T]\times[0,\kappa]$. The second difficulty is not insurmontable. In fact we don’t actually need to precisely describe the entirely domain $D_{R,T}$: we only need to know if the points of the Poisson process belong to it. This can fortunately be obtained in a quite simple way, that’s why we propose two modifications of the theoretical Algorithm (A0) – Algorithm (A1) and (A2) – which are easy to implement. These modifications essentially concern the way used to simulate the Poisson point process.

Algorithm (A1).

Let us fix the level $L>x$.

Step 1: Simulate a non-negative random variable $T$ with p.d.f. $f_{T}$.

Initialization: $\beta=(0,0,0)$, $\mathcal{W}=0$, $\mathcal{E}_{0}=0$ and $\mathcal{E}_{1}$ an exponentially distributed random variable of average $1/\kappa$.

Step 2: While $\mathcal{E}_{1}\leq T$ and $\mathcal{W}=0$,

Step 2.1: Simulate three independent random variables: a $3$-dimensional gaussian random variable $G$, an exponentially distributed variable $e$ with average $1/\kappa$ and a standard uniform random variable $V$.

Step 2.2: Compute*

[TABLE]

Change the value of $\mathcal{W}$, $\mathcal{E}_{0}$ and $\mathcal{E}_{1}$ in the following way:

[TABLE]

$\mathcal{E}_{0}\leftarrow\mathcal{E}_{1}$ and $\mathcal{E}_{1}\leftarrow\mathcal{E}_{1}+e$.

Step 3: If $\mathcal{W}=0$ then set $Y=T$ otherwise go to Step 1.

Outcome: the r.v. $Y$.**

In order to present the second modified algorithm, we introduce different notations. Let us consider a finite set $\mathcal{S}\subset\mathbb{R}\times\mathbb{R}^{3}$ with the property: if both $(t,x)$ and $(t,y)$ belong to $\mathcal{S}$ then $x=y$. Let us moreover define $\underline{\mathcal{S}}=\inf\{t:\ (t,x)\in\mathcal{S}\}$ and $\overline{\mathcal{S}}=\sup\{t:\ (t,x)\in\mathcal{S}\}$. For any $t\in[\underline{\mathcal{S}},\overline{\mathcal{S}}]$ we introduce $\overleftarrow{t}:=(\overleftarrow{t_{1}},\overleftarrow{t_{2}})$ corresponding to the couple $(s,x)\in\mathcal{S}$ such that $\overleftarrow{t_{1}}=\sup\{s:(s,x)\in\mathcal{S},\ s\leq t\}$. Similarly we define $\overrightarrow{t}:=(\overrightarrow{t_{1}},\overrightarrow{t_{2}})$ the couple $(s,x)\in\mathcal{S}$ such that $\overrightarrow{t_{1}}=\inf\{s:(s,x)\in\mathcal{S},\ s>t\}$.

Algorithm (A2).

Let us fix the level $L>x$.

Step 1: Initialization: $\mathcal{E}=0$ and $\mathcal{W}=0$. Simulate a non-negative random variable $T$ with p.d.f. $f_{T}$. Set $\mathcal{S}=\{(0,(0,0,0)),(T,(0,0,0))\}$.

Step 2: While $\mathcal{E}\leq\kappa$ and $\mathcal{W}=0$,

Step 2.1: Simulate two independent random variables: an exponentially distributed variable $e$ with average $1/T$ and a uniformly distributed variable $U$ on the interval $[0,T]$. We change the value of $\mathcal{E}$: $\mathcal{E}\leftarrow\mathcal{E}+e$.

Step 2.2: If $\mathcal{E}\leq\kappa$ then we simulate a $3$-dimensional gaussian random variable $G$ with covariance ${\rm Id}_{3}$ and we define:*

[TABLE]

We increase the set $\mathcal{S}$ in the following way: $\mathcal{S}\leftarrow\mathcal{S}\cup\{(U,\beta)\}$ and achieve the test:

[TABLE]

Step 3: If $\mathcal{W}=0$ then $Y=T$ otherwise return to Step 1.

Outcome: $Y$.**

Proposition 1.4**.**

Under Assumption 1.1–1.3, if $T/(L-x)^{2}$ has the same distribution as $1/G^{2}$ where $G$ is a standard Gaussian r.v., then $Y$, the outcome of Algorithm (A1) (respectively Algorithm (A2)) and $\tau_{L}$, defined by (0.1) & (1.2), are identically distributed.

Proof.

Algorithms (A1) and (A2) are just two different modifications of Algorithm (A0).

Step 1. Let us first present a quite different expression for $\eta(t)$ defined by (1.5) (these arguments were already developed in [22]). Let us define

[TABLE]

for any non-negative functional $F$. Since the Bessel process $(R_{t},\,t\geq 0)$ has the same distribution as the norm of the 3-dimensional Brownian motion $(B_{t},\ t\geq 0)$, we obtain:

[TABLE]

where

[TABLE]

The Brownian path can be decomposed as follows:

[TABLE]

Here $(\beta_{s},0\leq s\leq 1)$ stands for a standard $3$-dimensional Brownian bridge with $\beta_{0}=\beta_{1}=0$. We obtain therefore a new expression for the functional $G$:

[TABLE]

Since the Brownian bridge is rotationally invariant, we observe $G(x)=G(\|x\|e_{1})$ where $e_{1}=(1,0,0)$. We deduce

[TABLE]

and in particular:

[TABLE]

Step 2. Using the new expression of the weight $\eta(t)$, we propose the same kind of rejection algorithm as (A1). It suffices therefore to consider the domain:

[TABLE]

and the associated event $N(D_{\beta,T})=0$ where $N$ is a Poisson point process on the state space $[0,T]\times[0,\kappa]$. As usual, there exist three methods of simulation for the Poisson process.

One method consists in first sampling the number of points in the rectangle using the Poisson distribution and secondly in choosing their position uniformly on the rectangle. An algorithm based on such a procedure would be time consuming and requires a huge memory space. It is not reasonable in practice. 2. 2.

The second method consists in the simulation of the point process in the infinite domain $[0,T]\times\mathbb{R}_{+}$: we simulate a sequence of independent couples $(e_{k},U_{k})_{k\geq 1}$, $e_{k}$ and $U_{k}$ being independent random variables. $e_{k}$ is exponentially distributed with average $1/T$ and $U_{k}$ is uniformly distributed on the interval $[0,T]$ for any $k\geq 1$. Let us define:

[TABLE]

The point process therefore corresponds to the couples $(e_{1},U_{1}),\ldots,(e_{1}+e_{2}+\ldots+e_{\eta-1},U_{\eta-1})$ belonging to the domain $[0,T]\times[0,\kappa]$. Let us note that $\eta=1$ corresponds to the case where no point of the Poisson process belongs to the domain. Algorithm (A2) is based on this simulation. 3. 3.

A third method consists in using similar arguments as those just described but we interchange the space and the time variables. That means that a cumulated sum of exponentially distributed variables with average $1/\kappa$ until it overcomes the level $T$ are used to describe the time variable and uniform variables on the interval $[0,\kappa]$ for the space variable: we obtain Algorithm (A1).

Once the Poisson process is realized on the domain $[0,T]\times[0,\kappa]$, it suffices to test each point in order to observe if it belongs or not to the domain $D_{\beta,T}$ defined by (1.11). For this test we use the classical Brownian bridge simulation as described in the fourth step of Algorithm (A1) or Step 2.2 in Algorithm (A2). ∎

2 Efficiency of the algorithms

Let us denote by $\mathcal{I}$ the number of iterations observed in order to simulate the hitting time $\tau_{L}$ defined by (0.1) & (1.2) and $\mathcal{N}_{1},\ldots,\mathcal{N}_{\mathcal{I}}$ the numbers of random points (Poisson process) used for each iteration. Of course the efficiency is directly linked to

[TABLE]

the total number of random variables. It is therefore quite challenging to modify the algorithms in order to reduce it. Focusing on the efficiency, we first describe the number $\mathcal{I}$, secondly we propose some modification in order to reduce it and finally present upper-bound for the total number of random variables $\mathcal{N}_{\Sigma}$.

2.1 Number of iterations

Let us first note that the distribution of the number of iterations $\mathcal{I}$ (number of random variates $T$ produced in order to obtain the outcome $Y$) does not depend on the choice of the algorithm (A0), (A1) or (A2). The variable $\mathcal{I}$ is exponentially large as soon as $L-x$ becomes large:

Proposition 2.1**.**

Under Assumptions 1.1–1.3, the following upper-bound holds:

[TABLE]

Proof.

Let us just note that $\mathcal{I}$ is geometrically distributed with

[TABLE]

Since $\tau_{L}<\infty$ a.s., (1.4) leads to $\mathbb{E}[\mathcal{I}]=\mathbb{E}[\eta(\tau_{L})]^{-1}=\exp\{\beta(L)-\beta(x)\}.$ Moreover, due to the condition $\gamma\leq\kappa$, we obtain

[TABLE]

where $G$ is a standard gaussian r.v. ∎

In the simulation framework, such an exponential large number of iterations seems crippling in practice for large values of $L$ and stimulates one to find interesting modifications of Algorithms (A1)–(A2).

Space splitting

In order to reduce the number of iterations, we could replace the simulation of $\tau_{L}$ for a diffusion starting in $x$ (denoted here by $\tau(x\to L)$) by the simulation of $k$ independent first-passage times $\tau(x+(i-1)(L-x)/k\to x+i(L-x)/k)$ for $i=1,\ldots,k$. It suffices then to take the sum of all these quantities. The global number of iterations for this new algorithm, called Algorithm ${\rm(A1)_{split}}$ or ${\rm(A2)_{split}}$, becomes:

[TABLE]

The previous inequality is related to the particular choice $k=(\lfloor(L-x)\sqrt{2\kappa}\rfloor+1)$. Such a procedure permits to replace an exponentially large number of iterations into a linear one, as $L$ becomes large.

Shifting the function $\gamma$

An other way to reduce $\mathcal{I}$ is to decrease the parameter $\kappa$. Such a trick is possible in the following situation.

Assumption 2.1**.**

There exist two constants $\kappa>\gamma_{0}>0$ such that

[TABLE]

The modification is based on the following observation: (1.4) can be replaced by

[TABLE]

$\eta_{0}(t)$ being defined in a similar way as $\eta(t)$ in (1.5) just by replacing $\gamma(t)$ by the shifted function $\gamma(t)-\gamma_{0}$. In other words, the upper-bound parameter $\kappa$ can be replaced by the shifted value $\kappa-\gamma_{0}$ and consequently the number of iterations is reduced.

Proposition 2.2**.**

Let Algorithm $(A1)_{\rm shift}$, respectively $(A2)_{\rm shift}$, defined in a similar way as Algorithm $(A1)$, resp. $(A2)$, just by replacing the function $\gamma(\cdot)$ by the shifted function $\gamma(\cdot)-\gamma_{0}$. Under Assumption 1.1, 1.2 and 2.1, if $T$ is an inverse gaussian random variable ${\rm IG}(\frac{L-x}{\sqrt{2\gamma_{0}}},(L-x)^{2})$ then the outcome of the modified algorithm $Y$ and $\tau_{L}$ defined by (0.1) & (1.2) are identically distributed. Moreover this modification reduces the number of iterations:

[TABLE]

Proof.

Let us consider a non-negative bounded function $f$ and introduce

[TABLE]

where $\tau_{L}$ is the Brownian first-passage time, that is $\tau_{L}$ and $(L-x)^{2}/G^{2}$ are identically distributed (here $G$ stands for a standard gaussian r.v.). Hence

[TABLE]

We deduce that $E_{f}/E_{1}=\mathbb{E}[f(T)]$ where $T\sim{\rm IG}(\frac{L-x}{\sqrt{2\gamma_{0}}},(L-x)^{2})$. Let us denote by $f_{T}$ the p.d.f. of the variable $T$. By Theorem 1.3, we get

[TABLE]

Let us recall that $\mathbb{E}_{x}[\eta(\tau_{L})]=\exp-\{\beta(L)-\beta(x)\}$. Then (2.3) permits to point out the required distribution. Let us just notice that using the inverse gaussian variable increases the conditional probability of acceptance:

[TABLE]

Consequently the probability of acceptance/rejection at each step changes: $\mathcal{I}_{\rm mod}$ is geometrically distributed and

[TABLE]

∎

Remark 2.3**.**

Using a many-to-one transformation, Michael, Schucany and Haas introduced a simple generator of inverse gaussian distributions (see, for instance, [10] p.149). We can apply the following procedure to simulate $T\sim{\rm IG}(\mu,\lambda)$: let $N$ be a standard gaussian r.v. and $U$ uniformly distributed on $[0,1]$, independent of $N$.

Set $X=\mu+\frac{\mu^{2}N^{2}}{2\lambda}-\frac{\mu}{2\lambda}\,\sqrt{4\mu\lambda N^{2}+\mu^{2}N^{4}}$. 2. 2.

If $U\leq\frac{\mu}{\mu+X}$ then $T=X$ else $T=\frac{\mu^{2}}{X}$.

2.2 Efficiency of Algorithm (A1)

In this section, we are looking for an upper-bound of the number of random variates $\mathcal{N}_{\Sigma}$ used in order to simulate one variate $\tau_{L}$. First we shall focus on the first iteration in Algorithm (A1) and secondly adapt the procedure to the whole number of variates.

Proposition 2.4**.**

Let $\gamma$ satisfy Assumptions 1.1–1.3. We assume moreover the existence of $C_{\gamma}>0$ and $r<1$ such that

[TABLE]

Then there exist two constants $M_{\gamma,1}>0$ and $M_{\gamma,2}>0$ such that the number of random points in the first iteration of Algorithm (A1), satisfies

[TABLE]

Proof.

Let us denote by $\mathbb{P}_{c}$ (resp. $\mathbb{E}_{c}$) the conditional probability measure (resp. expectation) given both $\tau_{L}=T$ and the Bessel bridge trajectory $(R_{t},0\leq t\leq T)$. We observe that

[TABLE]

where $(U_{k})_{k}$ is a sequence of uniformly distributed r.v. on $[0,\kappa]$ and $(\epsilon_{k})_{k}$ is a sequence of exponentially distributed r.v. with average $\frac{1}{\kappa}$, both sequences being independent. We denote $\mathcal{E}_{k}=\epsilon_{1}+\ldots+\epsilon_{k}$. Hence

[TABLE]

By similar computations, for $k\geq 1$,

[TABLE]

where $(U^{(1)},\ldots,U^{(k)})$ stands for a set of reordered uniform random variables and $\mathcal{E}_{k}$ is gamma distributed $\gamma(k,\kappa)$. Hence

[TABLE]

Consequently, using a change of variable, we obtain

[TABLE]

By integration by parts, (2.2) becomes $\mathbb{E}_{c}[\mathcal{N}_{1}]=H_{T}+\kappa I_{T}$ where

[TABLE]

Let us first study the expression $I_{T}$. Using the scaling property of the Bessel bridge, we have

[TABLE]

where $(\hat{R}_{w}^{T},\,0\leq w\leq 1)$ is a $3$-dimensional Bessel bridge with condition $\hat{R}_{1}^{T}=(L-x)/\sqrt{T}$. Let us observe that

[TABLE]

Here $(\bar{R})$ is a standard $3$-dimensional Bessel bridge. This stochastic ordering result can be proven using the following arguments: first the Bessel bridge $\hat{R}^{T}$ has the same distribution as the Euclidian norm of $\frac{(L-x)w}{\sqrt{T}}\,e_{1}+\beta_{w}$ where $\beta_{w}$ is a standard $3$-dimensional Brownian bridge, see (1.9) and (1.10). Secondly the triangle inequality leads to

[TABLE]

In order to prove the stochastic ordering (2.8) it suffices then to note that $\|\beta_{w}\|$ and $\bar{R}_{w}$ are identically distributed.

Let us now introduce some parameter $\alpha\in]0,1[$ and define

[TABLE]

Then

[TABLE]

By (2.8) we obtain:

[TABLE]

Using the agreement formula (see, for instance, Theorem 2 in [31]), we obtain

[TABLE]

Here $C_{3}=\sqrt{2}/\Gamma(3/2)$ and $\bar{\tau}=\tau+\hat{\tau}$ where $\tau$ is the first hitting time of the level $1$ for a $3$-dimensional Bessel process and $\hat{\tau}$ an independent copy of $\tau$. Combining classical inequalities and the explicit expression of the Laplace transform corresponding to the distribution of $\tau$ (see [27]) leads to

[TABLE]

for any $\lambda>0$. $I_{\nu}$ stands for the Bessel function of the first kind. In particular $I_{1/2}(x)=\sqrt{\frac{2}{\pi x}}\,\sinh x$. The particular choice $\lambda=T^{2\alpha}/2$ leads to

[TABLE]

There exists a constant $C_{\alpha}>0$ such that

[TABLE]

Let us go back to the upper-bound (2.2). By definition, $\gamma_{T}\geq C_{\gamma}|x-T^{1/2+\alpha}|^{-r}$ as soon as $T\geq(x+1)^{2/(1+2\alpha)}$ and therefore

[TABLE]

Let us recall that $T$ is the Brownian first-passage time through the level $L$ starting in $x$: $T$ and $\frac{(L-x)^{2}}{G^{2}}$ are identically distributed, $G$ being a standard gaussian variable. Taking now the expectation with respect to $T$, we get

[TABLE]

Since $r<1$, it is possible to choose $\alpha\in]0,1[$ small enough such that $(1+2\alpha)r<1$ and therefore $\mathbb{E}[G^{-(1+2\alpha)r}]<\infty$. Let us choose for instance $\alpha=\frac{1-r}{4r}$ that means $(1+2\alpha)r=\frac{1+r}{2}$ ∎

Proposition 2.5**.**

Let us consider Algorithm $(A1)_{\rm shift}$ with input variable $T\sim{\rm IG}(\frac{L-x}{\sqrt{2\gamma_{0}}},(L-x)^{2})$ as in Proposition 2.2. Then the number of random points associated to the first iteration satisfies:

[TABLE]

Proof.

The arguments are similar to those developed in the proof of Proposition 2.4. The formula (2.7) leads to $\mathbb{E}_{c}[\mathcal{N}_{1}]\leq 1+\kappa T$. Since $T$ is an inverse gaussian r.v. we obtain the announced upper-bound:

[TABLE]

∎

The main description of the algorithm efficiency encompasses both the information about the first iteration and the description of the number of iterations.

Theorem 2.6**.**

Let us denote by $\mathcal{N}_{\Sigma}$ the total number of random points used in the simulation of the diffusion first-passage time.

Under condition (2.5), the number of points of Algorithm (A1), in the context of Proposition 1.4, satisfies:

[TABLE]

$M_{\gamma,1}$* and $M_{\gamma,2}$ being the constants introduced in Proposition 2.4.* 2. 2.

The number of points of the modified algorithm $(A1)_{\rm shift}$, in the context of Proposition 2.2, satisfies:

[TABLE]

Let us notice in particular that the average number of random points is finite which actually plays an important role for numerical purposes.

Proof.

We shall just prove the first item, the second one is based on similar arguments. Let us recall that $\mathcal{I}$ is the number of iterations in Algorithm (A1) and $A$ stands for the event Acceptance in the acceptance-rejection method. Therefore

[TABLE]

Let us denote by $p=e^{-(\beta(L)-\beta(x))}$ and observe that $\mathbb{P}(A)=p$. Moreover we get the following upper bound

[TABLE]

We deduce that (2.2) becomes

[TABLE]

Proposition 2.4 permits to conclude the proof. ∎

2.3 Efficiency of Algorithm (A2)

In this section, $\mathcal{N}_{1}$ still represents the number of random points used for the first iteration of the algorithm.

Theorem 2.7**.**

Let $\gamma$ satisfies Assumptions 1.1–1.3 and the particular condition (2.5). We consider Algorithm (A2) with $T$ corresponding to the Brownian first-passage time to the level $L$, with $x<L$. Then there exist two constants $M_{\gamma,1}>0$ and $M_{\gamma,2}>0$ such that the total number of random points $\mathcal{N}_{\sigma}$ in Algorithm (A2) satisfies

[TABLE]

The statement of Theorem 2.7 is similar to Theorem 2.6 but the arguments of proof are slightly different.

Proof.

We recall that $\mathbb{P}_{c}$ stands for the conditional probability given the first-passage time $T$ and the Bessel bridge $(R_{t})$ ending with $R_{T}=L-x$. We observe that

[TABLE]

where $U_{1}$ is uniformly distributed on the interval $[0,T]$ and $e_{1}$ is an exponentially distributed random variable independent of $U_{1}$ and with average $1/T$. Hence (2.12) becomes

[TABLE]

Let us generalize these computations. We introduce $(e_{n},\ n\geq 0)$ a sequence of exponentially distributed random variables with average $1/T$, $\mathcal{E}_{n}:=e_{1}+\ldots+e_{n}$ and $(U_{n},\ n\geq 0)$ a sequence of uniform distributed independent r.v. on the interval $[0,T]$, all r.v. being independent. Moreover let us denote by $F_{c}(t)=\gamma(L-R_{t})$ for any $t\in[0,T]$ (deterministic function under the conditional distribution $\mathbb{P}_{c}$) and $G_{c}(y)=\frac{1}{T}\int_{0}^{T}1_{\{y\geq F_{c}(t)\}}\,dt$. For $k>1$, we get

[TABLE]

Let us note that, under $\mathbb{P}_{c}$, we have on one hand $\mathcal{E}_{k}\sim\gamma(k,\frac{1}{T})$ and on the other hand:

[TABLE]

where $(V^{(1)},V^{(2)},\dots,V^{(k-1)})$ are ordered variables corresponding to the sequence of uniform random on $[0,1]$: $(V_{1},V_{2},\ldots,V_{k-1}).$

[TABLE]

Here $\mathcal{G}_{c}$ stands for $\mathcal{G}_{c}(y)=\int_{0}^{1}G_{c}(ty)\,\mathrm{d}t$.

[TABLE]

We deduce

[TABLE]

Taking the expectation with respect to both $T$ and the Bessel bridge $(R_{t})$, we obtain

[TABLE]

Let us now consider some upper-bound. First we note that $0\leq G_{c}(x)\leq 1$ for any $x\in\mathbb{R}$. Let us now consider $\mathcal{G}_{c}(x)$ using a change of variable:

[TABLE]

where $(\hat{R}_{w}^{T},\,0\leq w\leq 1)$ is a $3$-dimensional Bessel bridge with condition $\hat{R}_{1}^{T}=(L-x)/\sqrt{T}$, as considered in the proof of Proposition 2.4. For $x>0$,

[TABLE]

By (2.15) and the previous inequalities, we obtain:

[TABLE]

Let us just note that the variable $\Pi_{T}$ depends on $x$. Using the obvious upper-bound $Te^{-T\gamma_{T}}\leq\int_{0}^{1}Te^{-T\gamma_{T}u}\,\mathrm{d}u$ and the definition of $\gamma_{T}$ in (2.9), we observe that all arguments developed in the proof of Proposition 2.4 and in Theorem 2.6, from (2.9) on, can be adapted to (2.16). ∎

3 When the diffusion hits the level $L$ before a given fixed time $t_{0}$

Both Algorithm (A1) and (A2) concern first-passage times for the diffusion $(X_{t},\,t\geq 0)$ defined in (1.2) under some particular conditions. The boundedness of $\gamma:=\frac{b^{2}+b^{\prime}}{2}$ on $]-\infty,L]$ is indisputably the most restrictive condition (Assumption 1.3). Unfortunately the result presented so far cannot be generalized to the particular situation when $\gamma$ takes some negative values. Nevertheless we can adapt the algorithms by considering $\tau_{L}|\tau_{L}\leq t_{0}$ instead of $\tau_{L}$, defined by (0.1), for any fixed time $t_{0}$.

Let us consider a modification of the theoretical Algorithm $(A_{0})$. The idea is to focus our attention to random times which are a.s. smaller than some given value $t_{0}$.

Assumption 3.1**.**

There exists two constants $\kappa>0>-m$ such that

[TABLE]

Algorithm (A3).

Let us fix the level $L>x$, $x$ being the initial value of the diffusion.

Step 1: Simulate a non-negative random variable $T$ with p.d.f. $f_{T}$ of support $[0,t_{0}]$.

Step 2: Simulate a $3$-dimensional Bessel process $(R_{t})$ on the time interval $[0,T]$ with endpoint $R_{T}=L-x$. We define by $D_{R,T}^{m}$ the stochastic domain satisfying:*

[TABLE]

Step 3: Simulate a Poisson point process $N$ on the state space $[0,T]\times\mathbb{R}_{+}$, independent of the Bessel process, whose intensity measure is the Lebesgue one.

Step 4: If $N(D_{R,T}^{m})=0$ then set $Y=T$ otherwise go to Step 1.

Outcome: the random variable $Y$.**

It is obvious that the support of the random variable $Y$ is contained in $[0,t_{0}]$ by construction.

Proposition 3.1**.**

Let $\gamma$ a function satisfying Assumptions 1.1, 1.2 & 3.1 and $G$ a standard gaussian r.v. If $T$ has the same distribution as the conditional distribution of $(L-x)^{2}/G^{2}$ given $(L-x)^{2}/G^{2}\leq t_{0}$, then the random variable $Y$ obtained by Algorithm (A3) has the same distribution as $\tau_{L}$ given $\tau_{L}\leq t_{0}$.

Proof.

Step 1. Let us first describe the dependence of the distribution of the outcome $Y$ with respect to the distribution of the input $T$. Let us also note that, due to the condition on the function $\gamma$, $\gamma(L-R_{t})+\frac{mt_{0}}{T}\geq 0$ for any $t$ and $T\leq t_{0}$. We deduce that

[TABLE]

Following the arguments developed in Theorem 1.3, we obtain, for any non-negative or bounded function $\psi$,

[TABLE]

where $\eta$ is defined in (1.5). In other words, the p.d.f. of the outcome $Y$ is given by

[TABLE]

where $\Xi$ is the normalization coefficient. The result is quite similar to the statement of Theorem 1.3. The differences are the following: on one hand, the function $\gamma$ can take negative values, on the other hand the support of the distribution of $T$ has to be compact.

Step 2. Let us now link the conditional distributions of $Y$ and $\tau_{L}$. Girsanov’s transformation, see Proposition 1.1, implies

[TABLE]

where $\mathbb{P}$ (resp. $\mathbb{Q}$) corresponds to the distribution of the diffusion $X$ (resp. the Brownian motion $B$). The particular case $\psi\equiv 1$ leads to:

[TABLE]

Combining the ratio between (3.2) and (3.3) on one side and the equation (3) on the other side leads to

[TABLE]

where $T$ has the same law as $\tau_{L}|\tau_{L}\leq t_{0}$ under the probability measure $\mathbb{Q}$ (Brownian hitting time). In order to conclude, it suffices to note that, under $\mathbb{Q}$, the random variable $\tau_{L}$ and $(L-x)^{2}/G^{2}$ are identically distributed. ∎

Of course, for simulation purposes, we don’t use directly Algorithm (A3) since it is impossible to simulate the Poisson point process $N$ on the whole unbounded state space $[0,T]\times\mathbb{R}_{+}$. So we apply exactly the same procedures developed in Section 1 which transform Algorithm (A0) into (A1) or (A2). Such procedures modify Algorithm (A3) and permit to easily execute numerical simulations (Section 5.2).

4 Unbounded drift terms

The aim of this section is to handle with drift terms $b$ which do not satisfy the condition $\gamma(x)\leq\kappa$ for any $x\leq L$. From now on, we shall just assume that

[TABLE]

The algorithms (A1) and (A2) cannot be used in such a situation since the domain associated to the Poisson point process is unbounded. One way to overcome this difficulty is to replace the drift term of the diffusion by a modified version: for a given large parameter $\rho$,

[TABLE]

In other words, it suffices to keep the function $b$ on the interval $[-\rho,L]$ and to extend this function in a regular way to $]-\infty,L]$ in order to get a bounded function. Here $\lim_{x\to-\infty}b_{\rho}(x)=b(-\rho)$, the boundedness being verified.

Obviously $b_{\rho}$ is a $\mathcal{C}^{1}$-continuous function and both Algorithm (A1) and (A2) can be applied to the diffusion associated to this new drift term. The diffusion is denoted $X^{\rho}$.

Proposition 4.1**.**

Let the function $\gamma$ satisfy Assumption 1.1. Moreover we assume that $\lim_{x\to-\infty}p(x)=-\infty$ where $p$ is defined in (1.3). Then the distribution of $\tau_{L}^{\rho}$ converges towards the distribution of $\tau_{L}$ as $\rho\to\infty$. Moreover the Kolmorogov distance satisfies:

[TABLE]

Here $F_{Y}(x)$ stands for the cumulative distribution function of the r.v. $Y$.

In practice, in order to approximate $\tau_{L}$ we simulate $\tau_{L}^{\rho}$ using one of the algorithms presented before, with the parameter $\rho$ large enough. Let us notice that the condition $\lim_{x\to-\infty}p(x)=-\infty$ appearing in the statement of Proposition 4.1 implies Assumption 1.2 as mentioned in Remark 1.2.

Proof.

The arguments are quite straightforward: let us consider $x\in\mathbb{R}_{+}$, then

[TABLE]

Since $X^{\rho}$ and $X$ have the same distribution before exiting from the interval $[-\rho,L]$ (same drift term on this space subset), we obtain

[TABLE]

∎

The efficiency of this algorithm is already described in Section 2 and numerical results are presented in Section 5.3.

5 Examples and numerics

In this section we study the performances of the different algorithms proposed through the paper in order to simulate $\tau_{L}$ the first-passage time of the diffusion starting in $x$ with $x<L$. Let us just give a quick review in the following table describing which algorithm should be used when the function $\gamma$ satisfies particular bounds on the interval $]-\infty,L]$.

Condition on $\gamma$ r.v. simulated Algorithm and statement

$0\leq\gamma(x)\leq\kappa$ $\tau_{L}$ (A1) or (A2) in Proposition 1.4

$0<\gamma_{0}\leq\gamma(x)\leq\kappa$ $\tau_{L}$ ${\rm(A1)_{shift}}$ or ${\rm(A2)_{shift}}$ in Prop. 2.2

$-m\leq\gamma(x)\leq\kappa$ $\tau_{L}$ given $\tau_{L}\leq t_{0}$ (A3) in Prop. 3.1

$0\leq\gamma(x)$ $\tau_{L}^{\rho}$ (approx.) ${\rm(A1)^{\rho}}$ or ${\rm(A2)^{\rho}}$ in Prop. 4.1

Let us recall that the efficiency of all these algorithms can be increased using the space splitting argument proposed in Section 2.1.

In order to illustrate results stem from simulations, we consider the first hitting time of a Brownian motion with drift $\mu$. In this case, a closed form of the p.d.f. for $\tau_{L}$ is available (see, for instance [26] p. 197):

[TABLE]

Since the level $L$ is larger than the starting position $x$, Assumption 1.2 will be satisfied if and only if $\mu>0$.

In the opposite figure, there is a comparison between the theoretical probability density function of $\tau_{L}$ and its normalized histogram ($10\,000$ simulations) obtained by Algorithm (A1).

Here $b(x)\equiv 1$, $x=0$ and level $L=2$.

5.1 Algorithms (A1) and (A2)

We shall compare Algorithm (A1) and (A2), both of them simulating exactly the first-passage time $\tau_{L}$ for a diffusion process satisfying Assumption 1.1–1.3.

Example 1. We consider the following stochastic differential equation:

[TABLE]

We first observe that $\gamma(x)=(b^{2}(x)+b^{\prime}(x))/2=((2+\sin(x))^{2}+\cos(x))/2$ satisfies $0\leq\gamma\leq 5$. We use Algorithm (A1) in order to produce an histogram of the hitting time distribution (Figure 1).

The total number of random variables used in this algorithm is of course random and the associated histogram (Figure 2) illustrates its distribution.

Let us just note that Algorithm (A1) is quite time consuming since simulating $10\,000$ hitting times (corresponding to $L=2$ and $X_{0}=0$) requires about $300$ cpu. This time increases very fast as $L$ becomes large: the figure (right) presents the number of iterations versus $L$. In order to increase the algorithm efficiency and consequently reduce the computation time, one way is to shift the function $\gamma$ as described in Proposition 2.2. The function $\gamma$ introduced in (5.1) satisfies $\gamma(x)\geq\gamma_{0}:=1/4$. The comparison of the efficiency of Algorithm $(A1)$ (solid curve) and Algorithm ${\rm(A1)_{shift}}$ (dashed curve) is presented in Figure 2.

The example (5.1) also permits to illustrate the second modification of the algorithm developed in Section 2.1: the space splitting method. Let us recall that we can replace the simulation of $\tau_{L}$ for a diffusion starting in $x$ (denoted here by $\tau(x\to L)$) by the simulation of $k$ independent first-passage times $\tau(x+(i-1)(L-x)/k\to x+i(L-x)/k)$ for $i=1,\ldots,k$. It suffices then to take the sum of all these quantities. Figure 3 (left) represents the number of random variables for the particular case: $L=2$, $X_{0}=0$ and $k=20$. We can notice that such a procedure is slightly more efficient for this particular case: the average number of random variables (estimation with $10\,000$ simulations) namely goes from about $1791$ down to $102$. It permits also to deal with large level values (Figure 3 – right) which is a huge task with Algorithm (A1). For each level $L$ there is an optimal splitting (optimal value $k$) illustrates by Figure 4.

Up to this point, all simulations use Algorithm (A1). Of course, it is also possible to deal with Algorithm (A2) in such a situation. We obtain the same hitting time distribution and can therefore compare the number of iterations and the number of random variables used. Let us notice that both Algorithms (A1) and (A2) are based on the simulation of a Poisson point process defined on the rectangle $[0,T]\times[0,\kappa]$: for one algorithm, the simulation starts with the point with smallest abscissa and then goes from one point to the next one by increasing the abscissa, for the other one, the order of simulation depends on the ordinate. It is possible of course to have a one-to-one correspondance between these two simulation technics: each point $(t,x)\in[0,T]\times[0,\kappa]$ of the Poisson process is then transformed as follows

[TABLE]

Consequently the same sequence of random points can be used for both Algorithms and it becomes easy to compare the total number of random points. Let us introduce $\mathcal{N}_{1}$ the total number of random variables used by Algorithm (A1) and respectively $\mathcal{N}_{2}$ by Algorithm (A2), then we define the difference

[TABLE]

In the diffusion case (5.1), the estimated mean $\overline{\Delta}$ is equal to $-196.0$ and its standard deviation equals $1853.2$ when observing $10\,000$ simulations with $x=0$ and $L=2$. We claim therefore that Algorithm (A2) is more efficient than Algorithm (A1) for Example 1 since $\overline{\Delta}$ is significantly negative (p-value: $1.9E-26$). However this difference is small in comparison to the mean number of random variables: $\overline{\Delta}/\overline{\mathcal{N}}_{1}\approx-0.11$. In other situations, this difference can be much stronger or positive. It essentially depends on the drift term $b$ of the diffusion. Let us consider a second situation.

Example 2. Let us consider the following SDE:

[TABLE]

associated to the first-passage time for $L=1$. Let us note that $\tau_{1}<\infty$ a.s.

We observe that $\gamma(x)=\frac{1}{2}\left((1+\arctan(1-x))^{2}-\frac{1}{1+(1-x)^{2}}\right)$ satisfies $0\leq\gamma\leq\frac{1}{2}(1+\pi/2)^{2}$. Both Algorithms (A1) and (A2) can be used in order to simulate the first-passage time $\tau_{L}$ and the same kind of study is achieved: here, the estimated difference $\overline{\Delta}=-169.7$, its standard deviation $235.2$ and finally the reduction rate $\overline{\Delta}/\overline{\mathcal{N}}_{1}\approx-0.67$ for $x=0$, $L=1$ and $10\,000$ simulations. The choice between $(A1)$ and $(A2)$, in such a situation, becomes an important question (the hitting time distribution is illustrated by Figure 5).

5.2 Simulation of $\tau_{L}$ given $\tau_{L}\leq t_{0}.$

In this section, we illustrate Algorithm (A3) presented in Section 3 which permits to overcome the condition $\gamma\geq 0$. Indeed this condition can be weaken (Assumption 3.1) provided the study only concerns the distribution of the first-passage time $\tau_{L}$ conditioned on the event $\tau_{L}\leq t_{0}$. Here $t_{0}$ is any fixed time.

We therefore introduce the following SDE:

[TABLE]

We focus our attention to the case $X_{0}=0$, $L=1$ and different values for $t_{0}$. Let us observe that $\gamma(x)=(\arctan(x)^{2}-1/(1+x^{2}))/2$ satisfies $-1/2\leq\gamma(x)\leq\pi^{2}/8$ and that the first-passage time is almost surely finite (Assumption 1.2). As described in the statement of Proposition 3.1, the simulation of $\tau_{L}$ given $\tau_{L}\leq t_{0}$ is based on a two-steps acceptance/rejection method:

• •

In order to simulate $G^{2}$ given $G^{2}\geq(L-x)^{2}/t_{0}$ (here $G$ stands for a standard gaussian variable), we use shifted exponentially distributed random variables with an optimal (classical) rejection rule.

• •

Using the acceptance/rejection procedure introduced in Algorithm (A3) we transform the distribution of $(L-x)^{2}/G^{2}$ given $(L-x)^{2}/G^{2}\leq t_{0}$ into $\tau_{L}$ given $\tau_{L}\leq t_{0}$.

In this illustration, we choose for the Poisson point process appearing in Algorithm (A3) the same kind of simulation as Algorithm (A1). Of course, it is also possible to adapt the procedure in order to follow the technics of Algorithm (A2) and finally to compare the efficiency of both Algorithms (we let this task to the reader).

5.3 Modified Algorithm ${\rm(A1)}^{\rho}$

In this section we consider diffusion processes $(X_{t},\ t\geq 0)$ satisfying (1.1) with unbounded drift terms $b$. We apply the modified algorithm described in Section 4 in order to approximate the first-passage time through the level $L$. Let us consider as example an Ornstein-Uhlenbeck process with $b(x)=-\alpha x+\beta$ where $\alpha$ and $\beta$ are positive constants. We compute easily

[TABLE]

and observe that in general $\gamma$ is neither a bounded function on the whole interval $]-\infty,L]$ nor a positive one. The following choice of parameters $\alpha=0.3$, $\beta=1$ with starting position $X_{0}=0$ and boundary $L=1$ ensures that $\gamma$ is a positive function but $b$ remains unbounded, that’s why we replace the original drift term by its modified version (4.2):

[TABLE]

The modified $\gamma$ satisfies $\gamma_{\rho}(x)=\gamma(x)$ for $x\in[-\rho,L]$ and

[TABLE]

The function $\gamma_{\rho}$ is now positive on the whole interval $]-\infty,L]$ and admits the following upper-bound:

[TABLE]

To obtain such an expression, it suffices to upper-bound separately both main terms of (5.7). We can therefore apply Algorithm (A1) or Algorithm (A2) in order to simulate the approximated first-passage time $\tau_{L}^{\rho}$.

Figure 7 presents a comparison between the histograms of the hitting time distribution ($10\,000$ independent simulations and $\rho=5$) with the approximation of the hitting time density obtained via the numerical algorithm proposed in [8]. We can notice that the histogram and the numerical pdf perfectly fit. Moreover, the choice of different larger values for $\rho$ does not significantly affect the results.