Exit problem for Ornstein-Uhlenbeck processes: a random walk approach

TL;DR

This paper extends a random walk-based algorithm to efficiently approximate the exit time of Ornstein-Uhlenbeck processes, generalizing a method previously used for Brownian motion.

Contribution

It introduces a generalized Walk on Moving Spheres algorithm for Ornstein-Uhlenbeck processes and analyzes its efficiency.

Findings

The algorithm effectively approximates exit times for Ornstein-Uhlenbeck processes.

The method demonstrates computational efficiency and accuracy.

Extension of Brownian motion techniques to Ornstein-Uhlenbeck processes.

Abstract

In order to approximate the exit time of a one-dimensional diffusion process, we propose an algorithm based on a random walk. Such an algorithm so-called Walk on Moving Spheres was already introduced in the Brownian context. The aim is therefore to generalize this numerical approach to the Ornstein-Uhlenbeck process and to describe the efficiency of the method.

Equations

$$d{X}_t=b(t,X_t)dt+\sigma (t,X_t)d{W}_t,t\ge 0,$$

$$\mathcal{T}=\inf \{t\ge 0:X_t\notin [a,b]\}.$$

$${\psi }_{\pm }(t)=\pm \sqrt{t\log \left(\frac{d^2}{t}\right)},\text{for }t\in [0,d^2],$$

$$p(t)=\frac{1}{d\sqrt{2\pi }}\sqrt{\frac{1}{t}\log \left(\frac{d^2}{t}\right)},t\ge 0.$$

$$|{\psi }_{\pm }(t)|⩽\frac{d}{\sqrt{e}},\forall t\in [0,d^2].$$

$$d{X}_t=-\theta (X_t-\mu )dt+\sigma d{W}_t,t\ge 0,$$

$$d{X}_t=b(t,X_t)dt+\sigma (t,X_t)d{W}_t,t\ge 0.$$

$$|\sigma (s,x)-\sigma (s,y){|}^2⩽\rho (|x-y|),\forall x,y\in ]0,+\infty [,\forall t\in {\mathbb{R}}^{+}.$$

$$|b(s,x)-b(s,y)|⩽K_t|x-y|,\forall x,y\in H,s⩽t$$

$$X_t=X_0{e}^{-\theta t}+\mu (1-e^{-\theta t})+\sigma e^{-\theta t}{\int }_0^t{e}^{\theta s}d{W}_s,t⩾0.$$

$$X_t=X_0{e}^{-\theta t}+\mu (1-e^{-\theta t})+\frac{\sigma e^{-\theta t}}{\sqrt{2\theta }}{V}_{e^{2\theta t}-1}.$$

$${\psi }_{\mathrm{OU}}^{\pm }(t,x)=e^{-\theta t}\left(\frac{\sigma }{\sqrt{2\theta }}{\psi }_{\pm }(e^{2\theta t}-1)+x\right),$$

$${\tau }_{\mathrm{OU}}=\frac{\log (\tau +1)}{2\theta }a.s.$$

$${\tau }_{\mathrm{OU}}$$

$$=\inf \left\{t>0:x{e}^{-\theta t}+\frac{\sigma }{\sqrt{2\theta }}{e}^{-\theta t}{V}_{e^{2\theta t}-1}\notin [{\psi }_{\mathrm{OU}}^{-}(t,x),{\psi }_{\mathrm{OU}}^{+}(t,x)]\right\}$$

$$=\inf \left\{t>0:\frac{\sigma }{\sqrt{2\theta }}{e}^{-\theta t}{V}_{e^{2\theta t}-1}\notin [{\psi }_{\mathrm{OU}}^{-}(t,x)-x{e}^{-\theta t},{\psi }_{\mathrm{OU}}^{+}(t,x)-x{e}^{-\theta t}]\right\}$$

$$=\inf \left\{\frac{\log (u+1)}{2\theta }>0:V_u\notin [{\psi }_{-}(u),{\psi }_{+}(u)]\right\}=\frac{\log (\tau +1)}{2\theta }.$$

$$d=\left\{\begin{array}[]{c @{\text{ if }} l}\sqrt{2\theta e}\min\left(\frac{(b_{\gamma,x}-x)}{\sigma},\frac{2(x-a_{\gamma,x})}{\sqrt{\sigma^{2}+4\theta e\,x(x-a_{\gamma,x})}+\sigma}\right)&x\geqslant 0\\ \sqrt{2\theta e}\min\left(\frac{(x-a_{\gamma,x})}{\sigma},\frac{2(b_{\gamma,x}-x)}{\sqrt{\sigma^{2}-4\theta e\,x(b_{\gamma,x}-x)}+\sigma}\right)&x\leqslant 0\end{array}\right.$$

$$-\frac{\sigma d}{\sqrt{2\theta e}}+\frac{x}{\sqrt{1+d^2}}⩽{\psi }_{\mathrm{OU}}^{-}(t,x)⩽{\psi }_{\mathrm{OU}}^{+}(t,x)⩽\frac{\sigma d}{\sqrt{2\theta e}}+x.$$

$${\psi }_{\mathrm{OU}}^{-}(t,x)⩾-\frac{\sigma d}{\sqrt{2\theta e}}+\frac{x}{\sqrt{1+d^2}}⩾-\frac{\sigma d}{\sqrt{2\theta e}}+x(1-\frac{d^2}{2}).$$

$$P(d)=x\frac{d^2}{2}+\frac{\sigma }{\sqrt{2\theta e}}d+(a_{\gamma ,x}-x).$$

$$d_l=\frac{1}{x}\sqrt{\frac{{\sigma }^2}{2\theta e}+2x(x-a_{\gamma ,x})}-\frac{\sigma }{x\sqrt{2\theta e}}.$$

$$d_u=(b_{\gamma ,x}-x)\frac{\sqrt{2\theta e}}{\sigma }.$$

$$-\frac{\sigma d}{\sqrt{2\theta e}}⩽{\psi }_{\mathrm{OU}}^{-}(t,0)⩽{\psi }_{\mathrm{OU}}^{+}(t,0)⩽\frac{\sigma d}{\sqrt{2\theta e}}.$$

$$\mathbb{E}[N_{ϵ}]⩽\delta |\log (ϵ)|,\forall ϵ⩽{ϵ}_0.$$

$${\mathbb{E}}_x[N]⩽U(x),\forall x\in K.$$

$$\frac{{\sigma }^2}{2}{u}^{\prime \prime }-\theta x{u}^{\prime }=\frac{-1}{(x-a{)}^2(x-b{)}^2},\text{for }x\in ]a,b[.$$

$$u^{\prime }(x)=C(x)e^{\frac{2\theta }{{\sigma }^2}x}.$$

$$C(x)=-\frac{2}{{\sigma }^2}{\int }_0^x\frac{e^{-\frac{2\theta }{{\sigma }^2}s}}{(s-a{)}^2(s-b{)}^2}ds.$$

$$u(x)=-\frac{2}{{\sigma }^2}{\int }_0^x{e}^{\frac{2\theta }{{\sigma }^2}u}{\int }_0^u\frac{e^{-\frac{2\theta }{{\sigma }^2}s}}{(s-a{)}^2(s-b{)}^2}dsdu,\text{for }x\in ]a,b[.$$

$$T_n=\sum _{k=1}^n{\tau }_k^{\mathrm{OU}}$$

$$\mathbb{E}[Z_{n+1}-Z_n|{\mathcal{F}}_{T_n}]$$

Full text

Key words and phrases: Exit time, Ornstein-Uhlenbeck processes, generalized spheroids, WOMS algorithm.

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

1 Introduction

Simulating the first exit time for a diffusion from a given domain is primordial since these times appear in many domains. In mathematical finance, for instance, studying barrier options requires to estimate if the underlying stock price stays in a given interval. In the simple Black-Scholes model, the distribution of the first exit time is therefore well-known. In more complex models corresponding to general diffusion processes, such an explicit expression is not available and requires the use of numerical approximations.

Several methods have been introduced in order to approximate first exit times. The classical and most common approximation method is the Euler–Maruyama scheme based on a time discretization procedure. The exit time of the diffusion process is in that case replaced by the exit time of the scheme. The approximation is quite precise but requires to restrict the study on a given fixed time interval on one hand and to describe precisely the probability for the diffusion to exit inbetween two consecutive nodes of the time grid on the other hand.

In this study, we aim to introduce a random walk in order to approximate the diffusion exit time from a given interval. Let us introduce $(X_{t},\,t\geq 0)$ the unique solution of a stochastic differential equation:

[TABLE]

where $(W_{t},\,t\geq 0)$ stands for a one-dimensional Brownian motion. Let us also fix some interval $I=[a,b]$ which strictly contains the starting position $X_{0}=x$. We denote by $\mathcal{T}$ the diffusion first exit time:

[TABLE]

Our approach consists in constructing a random walk $(T_{n},X_{n})_{n\geq 0}$ on $\mathbb{R}_{+}\times\mathbb{R}$ which corresponds to a skeleton of the Brownian paths. In other words, the sequence $(T_{n},X_{n})$ belongs to the graph of the trajectory. Moreover we construct the walk in such a way that $(T_{n},X_{n})$ converges as time elapses towards the exit time and location $(\mathcal{T},X_{\mathcal{T}})$. It suffices therefore to introduce a stopping procedure in the algorithm to achieve the approximation scheme. Of course, such an approach is interesting provided that $(T_{n},X_{n})$ is easy to simulate numerically. For the particular Brownian case, the distribution of the exit time from an interval has a quite complicated expression which is difficult to use for simulation purposes (see, for instance [14]) whereas the exit distribution from particular time-dependent domains, for instance the spheroids also called heat balls, can be precisely determined. These time-dependent domains are characterized by their boundaries:

[TABLE]

where the parameter $d>0$ corresponds to the size of the spheroid. The first time the Brownian motion path $(t,W_{t})$ exits from the domain $\{(t,x):\ |x|\leq\psi_{+}(t)\}$, denoted by $\tau$, is well-known. Its probability density function [7] is given by

[TABLE]

It is therefore easy to generate such an exit time since $\tau$ and $d^{2}U^{2}e^{-N^{2}}$are identically distributed. Here $U$ and $N$ are independent random variables, $U$ is uniformly distributed on $[0,1]$ and $N$ is a standard gaussian random variable. Let us notice that the boundaries of the spheroids satisfy the following bound:

[TABLE]

This remark permits to explain the general idea of the algorithm. First we consider $(T_{0},X_{0})$ the starting time and position of the Brownian paths, that is $(0,x)$. Then we choose the largest parameter $d$ possible such that the spheroid starting in $(T_{0},X_{0})$ is included in the domain $\mathbb{R}_{+}\times[a,b]$. We observe the first exit time of this spheroid and its corresponding exit location, this couple is denoted by $(T_{1},X_{1})$. Due to the translation invariance of the Brownian motion, we can construct an iterative procedure, just considering $(T_{1},X_{1})$ like a starting time and position for the Brownian motion. So we consider a new spheroid included in the interval and $(T_{2},X_{2})$ shall correspond to the exit of this second spheroid and so on. Step by step we construct a random walk on spheroids also called WOMS algorithm (Walk On Moving Spheres) which converges towards the exit time and position $(\mathcal{T},W_{\mathcal{T}})$. This sequence is stopped as soon as the position $X_{n}$ is close enough to the boundary of the considered interval. The idea of this algorithm lies in the definition of spherical processes and the walk on spheres introduced by Müller [9] and used in the sequel by Motoo [8] and Sabelfeld [12] [13]. It permits also in some more technical advanced way to simulate the first passage time for Bessel processes [4].

In this study, we focus our attention on a particular family of diffusions which is strongly related to the Brownian motion: the Ornstein-Uhlenbeck processes. The idea is to use this link to adapt the Brownian algorithm in an appropriate way. This link implies changes on the time-dependent domains for which the exit problem can be expressed in a simpler way. We present the random walk algorithm (WOMS) for the Ornstein-Uhlenbeck process, describe the approximation error depending on the stopping procedure and emphasize the efficiency of the method. We describe the mean number of generalized spheroids necessary to obtain the approximated exit time.

2 The Ornstein-Ulhenbeck processes

Let us first recall the definition of the Ornstein-Uhlenbeck process and present different essential properties which permit to link this diffusion to a standard Brownian motion.

Let $\theta\in\mathbb{R}^{+}$, $\sigma\in\mathbb{R}^{+}$, $\mu\in\mathbb{R}$. The Ornstein-Uhlenbeck process (O.U.) starting in $x_{0}$ with parameters $\theta$, $\mu$, and $\sigma$ is the unique solution of the following stochastic differential equation (SDE):

[TABLE]

where $W$ stands for a standard one-dimensional Brownian motion. Existence and uniqueness for equation (2.1) can be easily deduced from a general statement concerning SDE, see for instance Revuz, Yor, Chap. IX [11]. Let us just recall this result.

Proposition 2.1**.**

Consider the following stochastic differential equation

[TABLE]

If there exists a Borel function $\rho:]0,+\infty[\rightarrow]0,+\infty[$ satisfying $\int_{0}^{+\infty}\frac{dx}{\rho(x)}=+\infty$ and such that

[TABLE]

and if, for each compact set $H$ and each $t\geqslant 0$, there exists a constant $K_{t}>0$ such that

[TABLE]

then pathwise uniqueness holds for equation (2.2).

Since obviously the drift and diffusion coefficients of the O.U. process satisfy the hypotheses of Proposition 2.1, pathwise uniqueness holds for (2.1). Let us now present an explicit expression of this solution. The Ornstein-Uhlenbeck process can be written as a stochastic integral with respect to the Brownian motion:

[TABLE]

Levy’s theorem permits to replace the stochastic integral by a time-changed Brownian motion. We obtain therefore another expression for the process which is more handy to manipulate.

Since $\theta>0$, there exists a standard Brownian motion $(V_{t})_{t\geq 0}$ such that

[TABLE]

This simplified expression is a crucial tool for the construction of the algorithm in the exit problem framework as it clearly appears in the forthcoming statements.

Remark 2.2**.**

In following computations, we put $\mu=0$. This restriction is only motivated by notational simplification and the study can easily be extended to the general case.

Let us now describe how such a strong relation between the Brownian motion and the Ornstein-Uhlenbeck process permits to emphasize a time-dependent domain of $\mathbb{R}$ whose exit time can be easily and exactly simulated.

3 Exit time of generalized spheroids

Let us consider the spheroids defined by the boundaries $\psi^{\pm}(t)$ in (1.1). We recall that the Brownian exit problem of a such a spheroid is completely explicit, so that the simulation of the exit time $\tau$ is rather simple. Due to the symmetry property of the spheroid, the conditional probability distribution of the exit location $W_{\tau}$ given $\tau$ is equal to $\frac{1}{2}\,\delta_{\psi^{+}(\tau)}+\frac{1}{2}\,\delta_{\psi^{-}(\tau)}$. For the Ornstein-Uhlenbeck process, we can obtain some similar information due to the strong relation with the Brownian motion.

Let us introduce two new boundaries defined by:

[TABLE]

where $\theta$ and $\sigma$ correspond to the parameters of the O.U-process $(X_{t})_{t\geq 0}$ in (2.1). We call generalized spheroid the domain defined by these boundaries.

We introduce the exit time $\tau_{\rm OU}=\inf\{t>0:\ X_{t}\notin[\psi_{\rm OU}^{-}(t,x),\psi_{\rm OU}^{+}(t,x)]\}$.

Proposition 3.1**.**

Let $\tau=\inf\{t>0:\ V_{t}\notin[\psi_{-}(t),\psi_{+}(t)]\}$ the first time the Brownian motion $(V_{t})_{t\geq 0}$ defined in (2.4) exits from the spheroid. Then the exit time $\tau_{OU}$ satisfies:

[TABLE]

Proof.

Using both the definition of $\tau_{\rm OU}$ and the expression of $X_{t}$ with respect to the Brownian motion $V_{t}$, we obtain

[TABLE]

∎

This statement is a crucial tool for simulation purposes. It permits first to simulate a Brownian exit time from a spheroid, then to use Proposition 3.1 to obtain the O.U. exit time from the generalized spheroid. Let us notice that the shape of the generalized spheroid depends on the O.U. starting position. Therefore, if we define a WOMS, the shape of the spheroids will change at each step of the algorithm. In the Brownian motion context, the spheroids are symmetric and their extremas can be computed easily. This important advantage permits to compute easily the maximal size of the spheroids included in the interval $[a,b]$ and is not fulfilled in the O.U. case. It is therefore an harder work to determine the optimal size of the generalized spheroid. This can be achieved by finding an upper-bound for the upper boundary and a lower-bound for the lower boundary. As a consequence, we determine a parameter characterizing the generalized spheroid which guaranties that it remains fully contained in the interval $[a,b]$. Since the bounds are quite rough, the boundaries of the generalized spheroid are unfortunately not tangent to the interval bounds. The algorithm shall be therefore a little slowed down.

Proposition 3.2**.**

Let $\gamma>0$, and $x\in[a,b]$ the starting point of the spheroid, that is $\psi^{\pm}_{OU}(0,x)=x$. Let us set $a_{\gamma,x}=a+\gamma(x-a)$ and $b_{\gamma,x}=b-\gamma(b-x)$. We define

[TABLE]

For such a choice of parameter, the generalized spheroid is fully contained in the interval $[a_{\gamma,x},b_{\gamma,x}]$.

In the following statements we denote by $d_{x}$ the parameter associated to the spheroid with initial point $x$.

Proof.

Let us first consider the case: $x>0$. Combining the upper bound of the function $\psi^{+}$ presented in (1.3) and the definition of $\psi_{\rm OU}$, we obtain

[TABLE]

We keep the upper bound found previously and focus on the lower bound:

[TABLE]

The determination of a convenient choice for the parameter $d>0$ requires to find the positive solution of the equation $P(d)=0$ where

[TABLE]

Consequently we obtain

[TABLE]

The identification with the upper bound gives us

[TABLE]

Hence setting $d=\min(d_{u},d_{l})$ permits the generalized spheroid to belong to the interval $[a_{\gamma,x},b_{\gamma,x}]$.

The case $x<0$ uses similar arguments since we observe a symmetry with respect to the origin between the generalized spheroid starting in $x$ and the one starting in $-x$. We use the results previously computed for $|x|$ and $[-b_{\gamma,x},-a_{\gamma,x}]$ which leads to the statement. The case $x=0$ is simple to handle with, since the previous boundaries (3.3) become

[TABLE]

It suffices to set $d=\displaystyle{\frac{\sqrt{2\theta e}}{\sigma}}\min(|a_{\gamma,0}|,b_{\gamma,0})$, which corresponds to the limit case as $x$ tends to [math] in both results previously established.

∎

4 WOMS for the Ornstein-Uhlenbeck processes

Let us now present the approximation procedure of the Ornstein-Uhlenbeck exit time from a given interval $[a,b]$. This algorithm is based on a walk on generalized spheroids (WOMS) described in the previous section.

ALGORITHM (O.U. WOMS)

Initialization: Let: $X_{0}=x_{0}$, $\mathcal{T}_{\epsilon}=0$

From step $n$ to step $n+1$:

While $X_{n}\leqslant b-\epsilon$ and $X_{n}\geqslant a+\epsilon$ do

$\bullet$

Generate the Brownian exit time from the spheroid with parameter $d_{X_{n}}$ defined in (3.2). We denote this stopping time by $\tau_{n+1}$.

$\bullet$

We set $\tau^{OU}_{n+1}=\frac{\log(\tau_{n+1}+1)}{2\theta}$.

$\bullet$

Generate a Bernoulli distributed r.v. $\mathcal{B}\sim\mathcal{B}(\frac{1}{2})$, if $\mathcal{B}=1$ then set $X_{n+1}=\psi_{OU}^{-}(\tau^{OU}_{n+1},X_{n})$ otherwise set $X_{n+1}=\psi_{OU}^{+}(\tau^{OU}_{n+1},X_{n})$.

$\bullet$

$\mathcal{T}_{\epsilon}\leftarrow\mathcal{T}_{\epsilon}+\tau_{n+1}^{OU}$.

Outcome: $\mathcal{T}_{\epsilon}$ the approximated O.U.-exit time from the interval $[a,b]$.

The CPU efficiency of such an algorithm shall be compared to the efficiency of classical approaches in the exit time approximation framework. Let us consider a particular situation: the exit time from the interval $[3,5]$ for the Ornstein-Uhlenbeck process starting in $4$ with $\theta=5$ and $\sigma=7$. We use an improved Euler method based on the correction by means of the sharp large deviations estimate of the exit probability. Such a method takes into account the probability for the diffusion path to exit inbetween two neighboring gridpoints (see the procedure described in [1]). The simulation of $100\,000$ samples with the step size $10^{-4}$ requires 64,7 seconds for this improved Euler method whereas the WOMS algorithm presented in this paper requires about 2,19 seconds for the corresponding choice $\epsilon=10^{-2}$ (here $\gamma=10^{-6}$).

Even if the study presented here concerns the exit time of some given interval $[a,b]$ denoted by $\tau_{[a,b]}$, let us just mention the possible link with first passage times (FPT). Intuitively for negative $a$ with large value $|a|$, the exit time of the interval can be approximated by the first passage time of the level $b$ denoted by $\tau_{b}$ i.e. $\lim_{|a|\to\infty}\mathbb{P}(\tau_{[a,b]}=\tau_{b})=1$. Several approaches permit to describe quite precisely the probability distribution of the Ornstein-Uhlenbeck FPT. In Figure 3, we illustrate that the distributions of both the exit time (histogram) and the first passage time (p.d.f.) present a thight fit. The histogram corresponds to the exit time obtained for an OU process starting in $-3$ with coefficients $\theta=1$ and $\sigma=1$ and observed on the interval with bounds $a=-10$ and $b=-1$. The curve corresponds to a numerical approximation of the first passage time density presented by Buonocore, Nobile and Ricciardi in [3]. An other approximation procedure for the FPT simulation is proposed by Herrmann and Zucca in [5]: it consists in simulating exactly the PFT of a slightly modified diffusion process. This modified diffusion has the following property: its drift term is bounded and coincides with the Ornstein-Uhlenbeck drift on the interval $[a,b]$ with $|a|$ large. Numerical comparisons permit to observe that the simulation of the exit time with the WOMS algorithm is highly more efficent than the method proposed in [5]: the simulation of a sample of size $100\,000$ takes a total of 3,7 seconds of CPU time with the first method and 197,2 seconds with the second one. Here the OU-process starts in $-3$ with coefficients $\theta=1$ and $\sigma=1$ and is observed on the interval $[-10,-1]$.

Let us now describe the WOMS algorithm for the Ornstein-Uhlenbeck process and especially emphasize its efficiency through theoretical results. We study how the strong relation between our process and the Brownian motion affects the statements obtained in the Brownian motion case. Let us just recall that the efficiency of the walk on spheres in the particular Brownian case is quite strong: the averaged number of steps is of the order $|\log(\epsilon)|$ (see for instance [2], for an overview of the convergence rate). In the Ornstein-Uhlenbeck case, we reach a similar efficiency result.

Average number of steps

Theorem 4.1**.**

Let $N_{\epsilon}$ be the random number of steps observed in the algorithm. Then there exist a constant $\delta>0$ and $\epsilon_{0}>0$ such that

[TABLE]

The statement is similar to the Brownian motion case, and the proofs are based on similar arguments. To prove this statement, we introduce a result coming from the potential theory and using Markov chains.

Let us consider a Markov chain $(X_{n})_{n\in\mathbb{N}}$ defined on a state space $I$ decomposed into two distinct subsets $K$ and $\partial K$, $\partial K$ being the so-called frontier. Let us define $N=\inf\{n\in\mathbb{N},X_{n}\in\partial K\}$ the hitting time of $\partial K$. We assume that $N$ is a.s. finite, then the following statement holds:

Proposition 4.2**.**

If there exists a function $U$ s.t. the sequence $(U(X_{n\wedge N}))_{n\in\mathbb{N}}$ is non negative and if the sequence $(U(X_{n\wedge N})+{n\wedge N})_{n\in\mathbb{N}}$ represents a super-martingale adapted to the natural filtration of the considered Markov chain $(X_{n})$, then

[TABLE]

The proof of this classical upper-bound is left to the reader, it is essentially based on the optimal stopping theorem and on the monotone convergence theorem (see, for instance,[10], p139).

Proof of Theorem 4.1.

Step 1. Let us first introduce a function $u$ which plays an important role in the construction of a super-martingale linked to the random walk.

We consider the following differential equation:

[TABLE]

This second order differential equation can be solved in a classical way. Let us first solve the related homogeneous equation: we obtain

[TABLE]

The method of variation of parameters leads to

[TABLE]

Integrating $u^{\prime}$ one more time implies an explicit expression of one particular solution (4.2).

[TABLE]

Step 2. We consider now the sequence $(T_{n},X_{n})_{n\in\mathbb{N}}$ of cumulative exit times, i.e.

[TABLE]

and exit location given by the WOMS algorithm for the Ornstein-Uhlenbeck process.

Let us introduce $Z_{n}=u(X_{n})+cn$ where $c$ is a positive constant (which shall be determined in the following calculus) and $u$ is the function detailed in Step 1 of the proof. We shall prove that this process is a super-martingale with respect to the filtration $(\mathcal{F}_{T_{n}})_{n\in\mathbb{N}}$ induced by $(\mathcal{F}_{t})$, the natural filtration of the Brownian motion $(V_{t})_{t\geqslant 0}$ enlightened in (2.4).

By Itô’s formula we obtain

[TABLE]

where $(M_{n})_{n\in\mathbb{N}}=\left(\int^{T_{n}}_{0}\sigma u^{\prime}(X_{s})dW_{s}\right)_{n\in\mathbb{N}}$ is a martingale and $D_{[a,b]}(x)=(x-a)(b-x)$ for $x\in[a,b]$. Remark now that

[TABLE]

where $\tilde{X}_{s}:=X_{T_{n}+s}$ has the same distribution as the Ornstein-Uhlenbeck starting in $X_{n}$. We now upper bound this term: we consider in a first time that $X_{n}$ is positive. By Proposition 3.2 we are then allowed to compute the corresponding coefficient $d_{X_{n}}$ which we denote by $d_{n}>0$ for notation simplicity. Let us fix some parameter $\Delta\in]0,1[$.

First case: $d_{n}\leqslant\Delta$, that is satisfied either if

[TABLE]

or

[TABLE]

with $b_{\gamma}=b_{\gamma,X_{n}}$ and $a_{\gamma}=a_{\gamma,X_{n}}$.

We first consider that $X_{n}$ is close enough to $b_{\gamma}$. Using (3.3), we have for any $t\in\text{Supp}\,(\psi_{OU}^{\pm})=\left[0,\frac{\log(1+d_{n}^{2})}{2\theta}\right]$:

[TABLE]

Since $d_{n}\leq\Delta<1$, we have $d_{n}^{2}\leq d_{n}$. Moreover $X_{n}\leq b_{\gamma}$ so that

[TABLE]

The last upper-bound uses the definition of $d_{n}$ presented in Proposition 3.2 Hence we have

[TABLE]

We then write, using the fact that $\tau^{\rm OU}_{n+1}$ is independent of $\mathcal{F}_{T_{n}}$,

[TABLE]

where $\tau_{n}$ denotes the exit time for Brownian motion from the spheroid of parameter $d_{n}$. If $\tau$ denotes the Brownian exit time of the generalized spheroid of normalized size ($d=1$), then the scaling property of Brownian motion implies that $\tau_{n}$ and $d^{2}_{n}\tau$ are identically distributed. Hence, noticing that $\tau\leqslant 1$ and recalling that $d_{n}^{2}\leqslant 1$, we obtain

[TABLE]

In the considered case, we know that

[TABLE]

which implies

[TABLE]

In the other case ($X_{n}$ close to $a$) the arguments already used just above lead to a similar upper-bound. We observe for any $t\in\left[0,\frac{\log(1+d_{n}^{2})}{2\theta}\right]$:

[TABLE]

This upper bound leads to the same result as (4.10) just replacing $\beta$ by another positive constant $\tilde{\beta}$. Combining both inequalities, for $d_{n}$ smaller than $\Delta$, we get

[TABLE]

Second case: $d_{n}>\Delta$

In this case, we use the upper-bound:

[TABLE]

We deduce

[TABLE]

Both inequalities (4.11) and (4.13) suggest the existence of a constant $\tilde{c}>0$ such that $\Xi(X_{n})\leqslant-\tilde{c}$.

Finally, using the symmetry property of the considered spheroid, the case $x$ negative is treated similarly, leading to a positive constant $c$ such that

[TABLE]

In conclusion, the stochastic process $Z_{n}=u(X_{n})+cn$ is a super-martingale due to the combination of (4.5) and (4.14).

Step 3. In order to apply the optimal stopping theorem described in Proposition 2.6., we need on one hand that $(U(X_{n})+cn)_{n\geqslant 0}$ is a super-martingale but also on the other hand that $(U(X_{n}))_{n\geqslant 0}$ is a non negative sequence. For the first property we could choose $U=u+\kappa$, $u$ being the function introduced in (4.3) and $\kappa$ a constant. For the second property we need to have a non negative sequence, so we have to choose in a suitable way the constant $\kappa$. Let us note that the function $u$ satisfies $u(0)=0$ and is a concave function. So in order to obtain a positive function on the interval $[a_{\gamma,x},b_{\gamma,x}]$ it suffices to choose $\kappa=-\min(u(b_{\gamma,x}),u(a_{\gamma,x}))$.

Consequently we need to study the behavior of $u$ at the frontiers of $[a_{\gamma,x},b_{\gamma,x}]$ that is for $x=b-\epsilon$ and $x=a+\epsilon$. Putting $b_{\gamma}:=b_{\gamma,b-\epsilon}$, we obtain

[TABLE]

Using Taylor’s expansion

[TABLE]

where $\xi$ belongs to $[0,b_{\gamma}-s]$. Moreover

[TABLE]

where $c_{i},i\in\{1,2,3,4\}$ are positive constants and $c_{2}=c_{4}=\displaystyle{\frac{1}{(b-a)^{2}}}$.

[TABLE]

where $I_{i,j,k}=\displaystyle{\int^{b_{\gamma}}_{0}\frac{(b_{\gamma}-s)^{k}}{(b-s)^{i}(s-a)^{j}}ds}$. We can notice that

[TABLE]

and

[TABLE]

Let us bound the last integral, using once again the partial fraction decomposition

[TABLE]

As in the previous computations, it is possible to take an equivalent as $\epsilon$ tends to zero, that is there exists $\delta>0$ such that

[TABLE]

Combining (4.16) and (4.17), and taking an equivalent when $\epsilon$ tends to [math] leads to state that

[TABLE]

A similar computation on $u(a_{\gamma,a+\epsilon})$ gives us some $\tilde{D}^{\prime}$ and $\tilde{\delta}$. Setting $\kappa=\max(\tilde{D}+\delta,\tilde{D}^{\prime}+\tilde{\delta})$ and $U(x)=u(x)-\kappa$ permits to obtain the positivity of the sequence $(U(X_{n}))_{n\geq 1}$.

Step 4. The statement of the theorem is a direct consequence of the optimal stopping theorem Proposition 2.6. If $N_{\epsilon}$ is almost surely finite, then

[TABLE]

In order to finish the proof, it remains to justify that $N_{\epsilon}$ is almost surely finite. Since $b_{\gamma,x}-x\geq(1-\gamma)\epsilon$ and $x-a_{\gamma,x}\geq(1-\gamma)\epsilon$ for any $x\in[a+\epsilon,b-\epsilon]$, we deduce that there exists a strictly positive lower-bound $d_{\epsilon}$ such that $d_{X_{n}}\geq d_{\epsilon}$ for any $n$. Introducing $(s_{n})$ a sequence of independent and identically distributed random variables corresponding to Brownian exits of a unit spheroid, we deduce that $T_{n}$ is stochastically lower-bounded by

[TABLE]

Moreover $S_{n}$ tends to infinity almost surely as $n\to\infty$. By Lemma 4.3 and by construction, $T_{n}$ is stochastically inbetween $S_{n}$ and $\mathcal{T}$ (an almost surely finite random variable) for any $n\leq N_{\epsilon}$. The stopping rule $N_{\epsilon}$ is therefore almost surely finite. ∎

Lemma 4.3**.**

The sequence of cumulative times $(T_{n})_{n\geqslant 1}$ appearing in the algorithm and defined by (4.4) are stochastically smaller than $\mathcal{T}$ the first exit time of the Ornstein-Uhlenbeck process.

Proof.

We need to emphasize the link between the markov chain induced by the algorithm, denoted $((T_{n},X_{n}))_{n\in\mathbb{N}}$ with $(T_{0},X_{0})=(0,0)$, and a path of the Ornstein-Ulhenbeck process.

At the starting point of the Ornstein-Uhlenbeck trajectory, we introduce a spheroid of maximum size contained in the interval $[a,b]\times\mathbb{R}_{+}$. The intersection of this spheroid and the path corresponds to the point $(t_{1},z_{1})$. Then this construction leads us to state that $(t_{1},z_{1})$ has the same distribution as $(T_{1},X_{1})$. Hence, from $(t_{1},z_{1})$ we can construct a maximum size spheroid and consider the intersection $(t_{2},z_{2})$ between the trajectory after $t_{1}$ and this second spheroid. Once again we get from the construction that $(t_{2},z_{2})$ and $(T_{2},X_{2})$ are identically distributed. We can therefore step by step build a sequence $((t_{n},z_{n}))_{n\in\mathbb{N}}$ of intersections between the considered trajectory and the spheroids. We obtain that the skeleton of the trajectory $(t_{n},z_{n})_{n\in\mathbb{N}}$ and the sequence $(T_{n},X_{n})_{n\in\mathbb{N}}$ are identically distributed. By construction, we also note that $t_{n}\leqslant\mathcal{T}$ for all $n\in\mathbb{N}$, which implies the announced result. ∎

Bounds for the Exit-Time distribution

Let us now precise the rate of convergence for the algorithm based on the random walk. We should describe how far the outcome of the algorithm and the diffusion exit time are. We recall that the outcome depends on the parameter $\epsilon$.

Theorem 4.4**.**

We consider $0<\gamma<2$ and $\delta=\epsilon^{\gamma}$. We denote by $F$ the cumulative distribution function of the exit time from the interval $[a,b]$ and $F_{\epsilon}$ the distribution function of the algorithm outcome. Then for any $\rho>1$, there exists $\epsilon_{0}>0$ such that

[TABLE]

for all $t\in\mathbb{R}$ and $\epsilon\leq\epsilon_{0}$.

In other words, the precision of the approximation pointed out in Theorem 4.4 is characterized by the following error bound:

[TABLE]

Figure 5 presents the dependence of this bound with respect to $\epsilon$ and $\theta$, all other parameters being fixed.

Such a statement is directly related to properties of the Ornstein-Uhlenbeck process and its strong link with the Brownian motion.

Proof.

As in Lemma 4.3, we build step by step a sequence $((t_{n},z_{n}))_{n\in\mathbb{N}}$ of intersections between the path of the Ornstein-Uhlenbeck process and the spheroids in such a way that the sequences $((t_{n},z_{n}))_{n\geq 0}$ and $((T_{n},X_{n}))_{n\geq 0}$ are identically distributed.

If we introduce $N_{\epsilon}$ the stopping time appearing in the stopping procedure of the algorithm and $\tilde{N}_{\epsilon}=\inf\{n\in\mathbb{N},z_{n}\notin[a+\epsilon,b-\epsilon]\}$, the identity in law of those random variables yields. By construction, $t_{n}\leqslant\mathcal{T}$ for all $n\in\mathbb{N}$, where $\mathcal{T}$ stands for the diffusion first exit time from the interval $[a,b]$. This inequality remains true when $t_{n}$ is replaced by the random stopping time $t_{\tilde{N}_{\epsilon}}$.

Hence

[TABLE]

We focus on the first term of this upper bound. Using the strong Markov property, we obtain

[TABLE]

For any $y\in[a,a+\epsilon]\cup[b-\epsilon,b]$ we write

[TABLE]

We first consider the case $y\in[b-\epsilon,b]$, the previous inequality becomes

[TABLE]

In order to handle with the first term in the right hand side, we introduce $s$ the scale function of the O.-U.-process:

[TABLE]

It has been shown in Karatzas, 5.5 [6] that

[TABLE]

Since $y\in[b-\epsilon,b]$ and since the integrated function is non negative and increasing we obtain

[TABLE]

We now focus on the second term in the r.h.s. of (4.24): $\mathbb{P}_{y}(\mathcal{T}_{b}>\delta)\leqslant\mathbb{P}_{b-\epsilon}(\mathcal{T}_{b}>\delta)$ for all $y\in[b-\epsilon,b]$. We denote by $\tilde{X}$ the Ornstein-Uhlenbeck process starting in $b-\epsilon$. We obtain

[TABLE]

Let us assume that $b>0$. In this case, the following asymptotic property holds:

[TABLE]

Using the particular form of $\delta=\epsilon^{\gamma}$, we obtain

[TABLE]

A similar bound can be obtained for b negative and also for $y\in[a,a+\epsilon]$.

Finally combining this result with (4.22), (4.23) and (4.27) leads to the announced statement. ∎

Remark 4.5**.**

Let us note that all the results presented so far, that is the efficiency of the algorithm and the convergence rate, concern the family of Ornstein-Uhlenbeck processes with parameter $\mu=0$ in (2.1). It is straightforward to extend the statements to the general case: it suffices to replace the interval $[a,b]$ by a time-dependent interval $[a-\mu(1-e^{-\theta t}),b-\mu(1-e^{-\theta t})]$.