2.1. Step (i)

We are interested in the boundary crossing probability of a one-dimensional diffusion process Y, whose dynamics are governed by the following stochastic differential equation:

(1) \begin{equation} \begin{cases} \mathrm{d} Y(t) = \mu_Y(t,Y(t))\,\mathrm{d} t + \sigma_Y(t,Y(t))\,\mathrm{d} W(t), & t\in (0,1],\\ Y(0) = y_0, \end{cases}\end{equation}

where $\{W_t\}_{t\geq 0}$ is a standard Brownian motion process defined on a stochastic basis $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},\mathbb{P})$ and $y_0$ is constant. We assume that $\mu_Y\colon[0,1] \times \mathbb{R} \to \mathbb{R}$ and $\sigma_Y \colon [0,1] \times \mathbb{R} \to (0,\infty)$ satisfy the following conditions sufficient for the uniqueness and existence of a strong solution to (1) [Reference Ethier and Kurtz18, pp. 297–298]:

Condition 1. The functions $\mu_Y$ and $\sigma_Y$ are continuous in both variables and such that, for some $K < \infty$ , $x \mu_Y(t,x) \leq K(1 +x^2)$ and $\sigma_Y^2(t,x) \leq K(1 + x^2)$ for all $(t,x)\in [0,1]\times \mathbb{R}$ and, for any bounded open set $U \subset \mathbb{R}$ , there exists a $K_U < \infty$ such that

\begin{equation*} \lvert\mu_Y(t,x) - \mu_Y(t,y)\rvert + \lvert\sigma_Y^2(t,x) - \sigma_Y^2(t,y)\rvert \leq K_U\lvert x - y\rvert \end{equation*}

for all $t\in [0,1]$ , $x,y\in U$ .

For the unit-diffusion transformation ${\psi}_t(y) \,:\!=\, \int_0^y\sigma_Y(t,u)^{-1}\,\mathrm{d} u$ to be well defined, we will also assume the following.

By Itô’s lemma, the transformed process $X =\{X(t) \,:\!=\, {\psi}_t(Y(t))\}_{t\in[0,1]}$ is a diffusion process with a unit diffusion coefficient (see, e.g., [Reference Kloeden and Platen35, Section 4.7] or [Reference Rogers54]),

(2) \begin{equation} \begin{cases} \mathrm{d} X(t) = {\mu_X}(t,X(t))\,\mathrm{d} t + \mathrm{d} W(t),&t\in (0,1], \\ X(0) = x_0 \,:\!=\, {\psi}_0(y_0), \end{cases}\end{equation}

where ${\mu_X}(t,x) = (\partial_t {\psi}_t + {\mu_Y}/{\sigma_Y} - \tfrac12 \partial_x\sigma_Y)\circ {\psi}_t^{-1}(x)$ and ${\psi}_t^{-1}(z)$ is the inverse of $z = {\psi}_t(y)$ in y.

Denote by ${C}=C([0,1])$ the space of continuous functions $x \colon [0,1]\to\mathbb{R}$ equipped with the uniform norm $\|{x}\|_{\infty} \,:\!=\, \sup_{t\in[0,1]}\lvert x(t)\rvert$ . For a fixed $x_0 \in \mathbb{R}$ , consider the class

\begin{equation*} \mathcal{G}\,:\!=\, \{(\,f^-,f^{+})\colon f^{\pm} \in {C},\,f^-(0) < x_0 < f^+(0),\, \inf_{0 \leq t \leq 1}(\,f^+(t) - f^-(t)) > 0\}\end{equation*}

of pairs of functions from C and introduce the notation

\begin{equation*} S(\,f^-,f^{+})\,:\!=\, \{x \in {C}\colon f^-(t) < x(t) < f^+(t),\,t\in [0,1] \},\qquad (\,f^-,f^{+}) \in \mathcal{G}.\end{equation*}

The problem we deal with in this paper is how to compute the probability of the form $\mathbb{P}(Y \in S(g_0^-,g_0^{+}))$ for $(g_0^-,g_0^{+}) \in \mathcal{G}$ . Clearly, the desired probability coincides with $\mathbb{P}(X \in G)$ , $G\,:\!=\, S(g^-,g^{+})$ , where $g^{\pm}(t) \,:\!=\, {\psi}_t(g_0^{\pm}(t))$ , $t\in[0,1]$ . It is also clear that $(g^{-},g^{+}) \in \mathcal{G}$ due to Condition 2. Henceforth, we work exclusively with the process X and the boundaries $g^{\pm}$ .

2.2. Step (ii)

In the context of curvilinear boundary crossing probabilities, [Reference Fu and Wu19] approximated the Wiener process by discrete Markov chains with transition probabilities computed by first taking the values of the Wiener process transition densities on a lattice and then normalising these to obtain a probability distribution on that lattice (there is also a small adjustment of that distribution to perfectly match the first two moments of the original transition probabilities and the ‘discretised’ ones).

In the general diffusion case, due to the absence of closed-form expressions for the transition density of X, we use transition probabilities of the weak Taylor approximations of X to construct the approximating discrete process. As the second-order expansion is to be used, we further require that the following condition is met.

Condition 3. For any fixed $x \in \mathbb{R}$ , ${\mu_X}(\,\cdot\,,x) \in C^1([0,1])$ , and for any fixed $t \in [0,1]$ , ${\mu_X}(t,\,\cdot\,) \in C^2(\mathbb{R})$ . Moreover, for any $r>0$ , there exists a $K_r <\infty$ such that

\begin{equation*} \lvert {\mu_X}(t,x)\rvert + \lvert \partial_t {\mu_X}(t,x)\rvert + \lvert \partial_{x} {\mu_X}(t,x) \rvert + \lvert \partial_{xx} {\mu_X}(t,x) \rvert \leq K_r , \qquad t\in[0,1], \,\lvert x \rvert \leq r. \end{equation*}

Next, for any $n \geq 1$ , let $t_{n,k}\,:\!=\, k/n$ , $k=0,1,\ldots, n$ , be the uniform partition of [0,1] of rank $\Delta_{n}\,:\!=\,1/n$ . Introduce discrete scheme drift ${\beta}_{n,k}$ and diffusion ${\alpha}_{n,k}^{1/2}$ coefficients by setting, for $k=1,2,\ldots, n$ ,

(3) \begin{align} {\beta}_{n,k}(x) & \,:\!=\, \big({\mu_X} + \tfrac12\Delta_{n}(\partial_t{\mu_X} + {\mu_X}\partial_x{\mu_X} + \tfrac12\partial_{xx}{\mu_X})\big)(t_{n,k-1},x), \end{align}

(4) \begin{align} {\alpha}_{n,k}^{1/2}(x) & \,:\!=\, 1 + \tfrac12\Delta_{n}\partial_x {\mu_X} (t_{n,k-1},x). \end{align}

For each fixed $n \geq 1$ , the nth second-order weak Taylor approximation of the diffusion (2) is defined as the discrete-time process $\zeta_n \,:\!=\,\{\zeta_{n,k}\}_{k=0}^n$ , specified by $\zeta_{n,0} = x_0$ and

\begin{equation*} \zeta_{n,k}= \zeta_{n,k-1} + {\beta}_{n,k}(\zeta_{n,k-1})\Delta_{n} + \alpha_{n,k}^{1/2}(\zeta_{n,k-1}) \Delta_{n}^{1/2} Z_{n,k}, \qquad k=1,\ldots,n, \end{equation*}

where $\{Z_{n,k}\}_{k=1}^n$ is a triangular array of independent standard normal random variables. For more detail on weak Taylor approximations to solutions of stochastic differential equations, see, e.g., [Reference Kloeden and Platen35, Chapter 14]. Clearly, the conditional distributions of the increments $\zeta_{n,k}- \zeta_{n,k-1}$ given $\zeta_{n,k-1}$ are Gaussian, and so the transition probabilities of the discrete-time process $\zeta_{n}$ can be easily obtained.

2.3. Step (iii)

Next, we further approximate the discrete-time continuous state space process $\zeta_n$ with a discrete-time discrete-state-space Markov chain $\xi_n\,:\!=\,\{\xi_{n,k}\}_{k=0}^n$ , whose transition probabilities are based on the normalised values of the transition density of $\zeta_n$ .

To improve the convergence rate for our approximations to $\mathbb{P}(X \in G)$ , we construct our Markov chains $\xi_n$ choosing, generally speaking, different state spaces for each time step. Namely, the state space $E_{n,k}$ for $\xi_{n,k}$ , $k=0,1,\ldots,n$ , is chosen to be a lattice such that $g^{\pm}(t_{n,k}) \in E_{n,k}$ . This modification improves upon the Markov chain approximation suggested in [Reference Fu and Wu19], and is widely used for accelerating the convergence rate of numerical schemes in barrier option pricing [Reference Cheuk and Vorst13]. More precisely, the spaces $E_{n,k}$ are constructed as follows. Let $g_{n,k}^{\pm} \,:\!=\, g^{\pm}(t_{n,k})$ , $k=1,\ldots,n$ , and, for fixed $\delta\in (0,\tfrac12]$ and $\gamma>0$ , set

\begin{equation*} {w}_{n,k} \,:\!=\, \begin{cases} \dfrac{(g_{n,k}^+ -g_{n,k}^{-})/\Delta_{n}^{1/2 +\delta} }{\lfloor{\gamma(g_{n,k}^+ -g_{n,k}^{-})/\Delta_{n}^{1/2 +\delta}}\rfloor}, & 1 \leq k<n, \\ \dfrac{(g^+(1) -g^-(1))/\Delta_{n} }{\lfloor{\gamma(g^+(1) -g^-(1))/\Delta_{n} }\rfloor}, & k=n, \end{cases}\end{equation*}

assuming that n is large enough that the integer parts in all the denominators are non-zero. We set the time-dependent space lattice step sizes to be

(5) \begin{equation} h_{n,k} \,:\!=\, \begin{cases} {w}_{n,k}\Delta_{n}^{1/2 +\delta}, & 1 \leq k < n, \\ {w}_{n,n}\Delta_{n}, & k = n. \end{cases}\end{equation}

The state space for $\xi_{n,k}$ is the $h_{n,k}$ -spaced lattice

(6) \begin{equation} E_{n,k}\,:\!=\, \{ g_{n,k}^+ - j h_{n,k} \colon j \in \mathbb{Z}\}, \qquad k = 1,\ldots,n.\end{equation}

We also put $E_{n,0}\,:\!=\,\{x_0\}$ and $E_{n}\,:\!=\, E_{n,0}\times E_{n,1} \times \cdots \times E_{n,n}$ . Note that

(7) \begin{equation} \max_{1\leq k\leq n}\lvert {w}_{n,k} - \gamma^{-1}\rvert \to 0\quad \mbox{ as $n\to\infty$} ,\end{equation}

and $\gamma^{-1}\leq {w}_{n,k}\leq 2\gamma^{-1}$ for all $1\leq k\leq n$ (assuming, as above, that n is large enough).

Further, for $k=1,\ldots,n,$ recalling (3) and (4), we let

(8) \begin{equation} {\mu}_{n,k}(x) \,:\!=\, {\beta}_{n,k}(x)\Delta_{n}, \qquad \sigma^{2}_{n,k}(x)\,:\!=\, {\alpha}_{n,k}(x)\Delta_{n},\end{equation}

and define $\xi_{n,k}$ , $k=1,\ldots,n$ , $n\geq 1$ , as a triangular array of random variables, where each row forms a Markov chain with one-step transition probabilities given by

\begin{equation*} p_{n,k}(x,y) \,:\!=\, \mathbb{P}(\xi_{n,k} = y \mid \xi_{n,k-1}=x) = \varphi(y \mid x + {\mu}_{n,k}(x),\, \sigma^{2}_{n,k}(x))\frac{h_{n,k}}{{C_{n,k}}(x)}\end{equation*}

for $(x,y) \in E_{n,k-1}\times E_{n,k}$ , where $\varphi(x\mid\mu,\sigma^2)\,:\!=\, (2\pi \sigma^2)^{-1/2}\mathrm{e}^{-(x -\mu)^2/(2\sigma^2) }$ , $x\in\mathbb R$ , denotes the Gaussian density with mean $\mu$ and variance $\sigma^2$ , and

(9) \begin{equation} {C_{n,k}}(x) \,:\!=\, \sum_{y \in E_{n,k}}\varphi(y \mid x + {\mu}_{n,k}(x),\sigma^{2}_{n,k}(x))h_{n,k}, \qquad x \in \mathbb{R}.\end{equation}

2.4. Step (iv)

As an initial approximation for $\mathbb{P}(X \in G )$ we could take $\mathbb{P}(g_{n,k}^- < X(t_{n,k}) < g_{n,k}^+, k = 1,\ldots,n)$ , which can in turn be approximated using the Chapman–Kolmogorov equations by

\begin{equation*} \mathbb{P}(g_{n,k}^- < \xi_{n,k} < g_{n,k}^+, \, k = 1,\ldots,n) = \sum_{\boldsymbol{{x}}\in E_n^G} \prod_{k=1}^n p_{n,k}(x_{k-1},x_k),\end{equation*}

where $\boldsymbol{{x}} = (x_0,x_1,\ldots,x_n)$ and $E_{n}^G\,:\!=\, E_{n,0}\times E_{n,1}^G \times \cdots \times E_{n,n}^G$ ,

\begin{equation*} E_{n,k}^G\,:\!=\, \{x \in E_{n,k} \colon g^-_{n,k} < x < g^+_{n,k}\},\qquad k = 1,\ldots,n.\end{equation*}

Without loss of generality, we can assume that n is large enough that none of $E_{n,k}^G$ is empty.

An approximation of this kind was used in [Reference Fu and Wu19] to approximate the boundary crossing probability of the Brownian motion. Such discrete-time approximations of boundary crossing probabilities are well known to converge slowly to the respective probability in the ‘continuously monitored’ case since they fail to account for the probability of boundary crossing by the continuous-time process at an epoch inside a time interval between consecutive points on the time grid used.

In order to correct for this, so-called ‘continuity corrections’ have been studied in the context of sequential analysis [Reference Siegmund and Yuh58] and, more recently, in the context of barrier option pricing [Reference Broadie, Glasserman and Kou9]. These types of corrections have also been used in [Reference Pötzelberger and Wang50] to correct for discrete-time monitoring bias in Monte Carlo estimates of the boundary crossing probabilities of the Brownian motion. Without such a correction, using the classical result from [Reference Nagaev43, Reference Nagaev44], the convergence rate of the approximation from [Reference Fu and Wu19] was shown to be $O(n^{-1/2})$ . However, as our numerical experiments demonstrate, using our more accurate approximations of the transition probabilities in conjunction with the Brownian bridge correction greatly improves it from $O(n^{-1/2})$ to $O(n^{-2})$ .

In the case of standard Brownian motion, the correction consists of simply multiplying the one-step transition probabilities by the non-crossing probability of a suitably pinned Brownian bridge. Due to the local Brownian nature of a diffusion process, it was shown in [Reference Baldi and Caramellino3] that the leading-order term of the diffusion bridge crossing probability is given by an expression close to that of the Brownian bridge. Thus, to account for the possibility of the process X crossing the boundary inside time intervals $[t_{n,k-1},t_{n,k}]$ , we define a process ${\widetilde{X}_{n}} \,:\!=\, \{{\widetilde{X}_{n}}(t)\}_{t\in[0,1]}$ which interpolates between the subsequent nodes $(t_{n,k}, \xi_{n,k})$ with a collection of independent Brownian bridges:

(10) \begin{equation} {\widetilde{X}_{n}}(t) \,:\!=\, B_{n,k}^{ \xi_{n,k-1}, \xi_{n,k} }(t), \qquad t \in [t_{n,k-1},t_{n,k}], \ k=1,\ldots,n,\end{equation}

where $B_{n,k}^{x,y}(t) \,:\!=\, B_{n,k}^{\circ}(t) + x + n(t-t_{n,k-1})(y-x)$ , $x,y\in\mathbb R$ , and $B_{n,k}^{\circ}(t)$ , $t\in [t_{n,k-1},t_{n,k}]$ , are independent Brownian motions ‘pinned’ at the time–space points $(t_{n,k-1},0)$ and $(t_{n,k},0)$ , these Brownian bridges being independent of $\xi_n$ . Analogous to [Reference Novikov, Frishling and Kordzakhia45, Theorem 1], using the Chapman–Kolmogorov equations, the non-crossing probability of the boundaries $g^{\pm}$ for ${\widetilde{X}_{n}}$ can be expressed as

(11) \begin{equation} \mathbb{P}({\widetilde{X}_{n}} \in G)= \mathbb{E}\prod_{k=1}^{n}\left(1-{\pi}_{n,k}(g^-,g^+\mid \xi_{n,k-1},\xi_{n,k})\right),\end{equation}

where

\[ {\pi}_{n,k}(g^-,g^+\mid x,y) \,:\!=\, \mathbb{P} \bigg(\sup_{t\in [t_{n,k-1},t_{n,k}]} (B_{n,k}^{x,y}(t) -g^+(t))(B_{n,k}^{x,y}(t) -g^-(t))\ge 0 \bigg)\]

is the probability that the trajectory of a Brownian motion process pinned at the points $(t_{n,k-1},x)$ and $(t_{n,k},y)$ will be outside of the ‘corridor’ between the boundaries $g^{\pm}(t)$ at some point during the time interval $[t_{n,k-1},t_{n,k}]$ .

2.5. Step (v)

The above representation can be equivalently written as a matrix product:

(12) \begin{equation} \mathbb{P}( {\widetilde{X}_{n}} \in G ) = \textbf{T}_{n,1}\textbf{T}_{n,2}\cdots \textbf{T}_{n,n} \textbf{1}^\top,\end{equation}

where $\textbf{1}=(1,\ldots,1)$ is a row vector of length $\lvert E_{n,n}^G \rvert$ , and the sub-stochastic matrices $\textbf{T}_{n,k}$ of dimensions $\lvert E_{n,k-1}^G \rvert \times \lvert E_{n,k}^G \rvert$ have entries equal to the respective taboo transition probabilities

\begin{equation*} (1-{\pi}_{n,k}(g^-,g^+ \mid x,y))p_{n,k}(x,y), \qquad (x,y) \in E_{n,k-1}^G \times E_{n,k}^G.\end{equation*}

Unfortunately, closed-form expressions for curvilinear boundary Brownian bridge crossing probabilities ${\pi}_{n,k}(g^-,g^+\mid x,y)$ are known in a few special cases only, so we approximate the original boundaries $g^{\pm}$ with piecewise linear functions $f_n^{\pm}$ , which linearly interpolate between the subsequent nodes $(t_{n,k},g_{n,k}^{\pm})$ for $k=0,\ldots,n$ . In the special case of a one-sided boundary (when $g^-=-\infty$ ), the expression for the linear boundary crossing probability of the Brownian bridge is well known [Reference Borodin and Salminen7, p. 63]:

\begin{equation*} {\pi}_{n,k}({-}\infty ,f_n^+ \mid x,y) = \exp\bigg\{\frac{-2}{\Delta_{n}}(g_{n,k-1}^+ - x)(g_{n,k}^+ - y)\bigg\}, \qquad x<g_{n,k-1}^+,y<g_{n,k}^+.\end{equation*}

In the case when both the upper and lower boundaries $f_n^{\pm}$ are linear, if the time interval $\Delta_{n}$ is sufficiently small, we can approximate the Brownian bridge crossing probability with the sum of one-sided crossing probabilities:

\begin{equation*} {\pi}_{n,k}(\,f_n^-,f_n^+\mid x,y) = {\pi}_{n,k}(\,f_{n}^-,\infty \mid x,y) + {\pi}_{n,k}({-}\infty,f_{n}^+ \mid x,y) - \vartheta(x,y,\Delta_{n}),\end{equation*}

where the positive error term $\vartheta$ admits the obvious upper bound

\begin{align*} \vartheta(x,y,{\Delta_n}) & \leq \pi_{n,k}(\,f_n^-,\infty\mid x,y)\max_{t\in[t_{n,k-1},t_{n,k}]}P_{t,f_n^-(t);t_{n,k},y} \biggl(\sup_{s\in [t,t_{n,k}]}(W(s)-f_n^+(s))\geq 0\biggr) \\ & \quad + \pi_{n,k}({-}\infty,f_n^+\mid x,y)\max_{t\in[t_{n,k-1},t_{n,k}]}P_{t,f_n^+(t);t_{n,k},y} \biggl(\inf_{s\in [t,t_{n,k}]}(W(s)-f_n^-(s))\leq 0\biggr),\end{align*}

where $P_{s,a;t,b}(\cdot) \,:\!=\, \mathbb{P}(\,\cdot\mid W(s)=a,W(t)=b)$ . An infinite series expression for ${\pi}_{n,k}(\,f_n^-,f_n^+\mid x,y)$ can be found, e.g., in [Reference Anderson2, Reference Hall25].

We can further apply our method to approximate probabilities of the form $\mathbb{P}(X \in G, X(1) \in [a,b])$ for some $[a,b] \subseteq [g^-(1), g^+(1)]$ . We first replace the final space grid $E_{n,n}^G$ with

\begin{equation*} E_{n,n}^{[a,b]} \,:\!=\, \bigg\{a \leq x \leq b \colon x = b - j\frac{b-a}{\lfloor{\gamma (b-a)/{\Delta_n}}\rfloor}, j\in \mathbb{Z}\bigg\}.\end{equation*}

Then, instead of (12) we have $\mathbb{P}({\widetilde{X}_{n}} \in G,{\widetilde{X}_{n}}(1) \in [a,b]) = \textbf{T}_{n,1}\textbf{T}_{n,2} \cdots \textbf{T}_{n,n}\textbf{1}_{[a,b]}^{\top}$ , where $\textbf{T}_{n,n}$ is now of dimension $\lvert E_{n,n-1}^G\rvert \times \lvert E_{n,n}^{[a,b]}\rvert$ and $\textbf{1}_{[a,b]} = (1,\ldots,1)$ is a row vector of length $\lvert E_{n,n}^{[a,b]}\rvert$ .

2.6. Remarks

We conclude this section with several remarks aimed at clarifying the choices we made in our construction and the connections of our method with existing approaches in the literature.

Remark 1. Our main results actually hold in a more general setting where, instead of using uniform partitions, we employ general ones $0 = t_{n,0} < t_{n,1} <\cdots < t_{n,n} = 1$ subject to the condition that, for some constants ${\eta_2} \geq {\eta_1} > 0$ and all $n\ge 1$ ,

(13) \begin{equation} \frac{{\eta_1}}{n} \leq t_{n,k} - t_{n,k-1} \leq \frac{{\eta_2}}{n}, \qquad k=1,2,\ldots,n. \end{equation}

All the proofs extend to this case in a straightforward way.

Using a more general sequence of partitions satisfying (13) helps to improve convergence rates by choosing a higher frequency of partition nodes on time intervals where the second derivatives of the boundaries $g^{\pm}$ are large.

Remark 4. If we used $\delta = 0$ in (5), there would be no convergence of the sequence of processes $\{\xi_n\}_{n \geq 1}$ to the desired diffusion limit. This is because there would be no convergence of the moments of the increments (cf. Lemma 3). Note also that, instead of using the power function (5), we could take $h_{n,k}\,:\!=\, w_{n,k}({\Delta_n})^{1/2}\psi({\Delta_n})$ for some $\psi(x) \to 0$ as $x \downarrow 0$ , with $w_{n,k} \,:\!=\, v_{n,k}/\lfloor{v_{n,k}}\rfloor$ , where

\begin{equation*} v_{n,k} \,:\!=\, \frac{g_{n,k}^+ -g_{n,k}^-}{(\Delta_{n} )^{1/2}\psi(\Delta_{n} )}, \end{equation*}

(in our case, $\psi(x) = x^{\delta}/\gamma$ ). We chose the power function for simplicity’s sake.

5.1. The Wiener process with one-sided boundary

Using the method of images, [Reference Daniels14] obtained a closed-form expression for the crossing probability of the boundary

\begin{equation*} g_\mathrm{D}(t) \,:\!=\, \tfrac12 - t \ln \big(\tfrac{1}{4}(1 + \sqrt{1 + 8\mathrm{e}^{-1/t}})\big), \qquad t>0,\end{equation*}

for the standard Wiener process $W\,:\!=\,\{W(t) \colon t\geq 0\}$ .

In order to make the state space $E_{n}^G$ finite in the case of the one-sided boundary where $g^-(t) = -\infty$ , we insert an absorbing lower boundary at a low enough fixed level $L<x_0$ and replace $E_{n,k}^G$ with

\begin{align*} E_{n,k}^{G,L} & \,:\!=\, \{x \in E_{n,k} \colon L < x < g_{n,k}^+ \}, \qquad k=1,\ldots,n, \\ E_{n}^{G,L} & \,:\!=\, E_{n,0} \times E_{n,1}^{G,L} \times \cdots \times E_{n,n}^{G,L}.\end{align*}

We approximate $\mathbb{P}(W(t) < g_\mathrm{D}(t),\, t\in [0,1])$ with $\big(\prod_{k=1}^n\textbf{T}_{n,k}^L\big)\textbf{1}^\top$ , where $\textbf{T}_{n,k}^L$ are substochastic matrices of dimensions $\big(\big\lvert E_{n,k-1}^{G,L}\big\rvert + 1\big) \times \big(\big\lvert E_{n,k}^{G,L}\big\rvert + 1\big)$ with entries equal to the respective transition probabilities

\begin{equation*} \begin{cases} q_{n,k}(x,y), & (x,y) \in E_{n,k-1}^{G,L} \times E_{n,k}^{G,L}, \\ \sum_{z\in E_{n,k}\cap({-}\infty,L]}q_{n,k}(x,z), & x \in E_{n,k-1}^{G,L},\, y = L,\\ 1, & x = L,\,y = L,\\ 0, & \text{otherwise}, \end{cases}\end{equation*}

where we put $f_n^-(t) = -\infty$ in the definition of $q_{n,k}$ in (26). This approximation assumes that the lower auxiliary boundary L is sufficiently far away from the initial point $x_0$ and the upper boundary, such that after a sample path crosses the lower boundary it is highly unlikely that it will cross the upper boundary in the remaining time. In our example, we took $L=-3$ . The probability of the Wiener process first hitting this level and then crossing $g_\mathrm{D}$ prior to time $t=1$ is less than $1.26\times 10^{-6}$ . Further, we chose $x_0 =0$ , $\delta = 0$ , and $\gamma = 2$ . To guarantee convergence of the scheme, $\delta$ must be strictly positive; however, for the values of n we are interested in and the larger value of $\gamma$ compared to the one in the example from Fig. 1, the error is negligible.

Figure 2. Approximation of the boundary $g_\mathrm{D}$ non-crossing probabilities for the Wiener process. The exact Daniels boundary crossing probability in this case is $0.479\,749\,35\ldots$ The left pane shows a log–log plot of the absolute approximation error as a function of n. The right pane shows the time evolution of the taboo transition density using the Markov chain approximation ${\widetilde{X}_{n}}$ with $n = 20$ .

The left pane in Fig. 2 shows a log–log plot of the absolute approximation error as a function of n. The crosses $\times$ show the absolute error of the Markov chain approximation without the Brownian bridge correction, while the bullets $\bullet$ show the error when using the Brownian bridge correction. The upper and lower dashed lines correspond to $C_1n^{-1/2}$ and $C_2n^{-2}$ respectively, where $C_1$ and $C_2$ are fitted constants. We see that the convergence rate $\lvert \mathbb{P}(W\in G) - \mathbb{P}({\widetilde{X}_{n}} \in G_n)\rvert$ is of the order $O(n^{-2})$ , which is the same as the boundary approximation order of the error $\lvert \mathbb{P}(W \in G) - \mathbb{P}(W \in G_n)\rvert = O(n^{-2})$ proved in [Reference Borovkov and Novikov8] in the case of twice continuously differentiable boundaries. It appears that, due to the high accuracy of approximation of the increments’ moments (Lemma 3) and the Brownian bridge correction applied to the transition probabilities, we achieve a much faster convergence rate compared to the convergence rate $O(n^{-1/2})$ achieved in [Reference Fu and Wu19].

The right pane in Fig. 2 shows the time evolution of the taboo transition density using the Markov chain approximation ${\widetilde{X}_{n}}$ with $n = 20$ . The positions of the nodes on the surface correspond to the points from the respective $E_{n,k}^{G,-3}$ . Note from (5) that the spacing between the nodes at the final time step $t = 1$ is finer compared to earlier time steps. This is crucial for the observed improved convergence rate.

5.2. The Ornstein–Uhlenbeck process

Let X be the Ornstein–Uhlenbeck (OU) process satisfying the stochastic differential equation

\begin{equation*} \begin{cases} \mathrm{d} X(t) = -X(t)\,\mathrm{d} t + \mathrm{d} W(t), \qquad t \in (0,1], \\ X(0) = 0. \end{cases}\end{equation*}

The usual approach for computing the boundary crossing probability of the OU process is to express the process in terms of a time-changed Brownian motion. This is achieved by using the time substitution $\theta(t) \,:\!=\, (\mathrm{e}^{2 t}-1)/2$ , so that we can write $X(t) = \mathrm{e}^{-t}W(\theta(t))$ .

To illustrate the effectiveness of our approximation, we consider the following two-sided boundary for which explicit boundary crossing probabilities are available for the OU process:

\begin{equation*} g_{\psi}^{\pm}(t) \,:\!=\, \mathrm{e}^{-t}\psi_{\pm}(\theta(t)), \quad \mbox{where } \psi_{\pm}(t) = \pm \tfrac12 t \cosh^{-1}(\mathrm{e}^{4/t}), \qquad t>0.\end{equation*}

Letting $t \,:\!=\, \theta(s)$ , we obtain

\begin{align*} P_{\psi}(T) & \,:\!=\, \mathbb{P}(\psi_{-}(t) < W(t) < \psi_{+}(t), \, 0 \leq t \leq T) \\ & = \mathbb{P}(\mathrm{e}^{-s}\psi_{-}(\theta(s)) < \mathrm{e}^{-s}W(\theta(s)) < \mathrm{e}^{-s}\psi_{\pm}(\theta(s)), \, 0 \leq s \leq \theta^{-1}(T)) \\ & = \mathbb{P}(g_{\psi}^{-}(s) < X(s) < g_{\psi}^{+}(s), \, s \in [0,\theta^{-1}(T)]),\end{align*}

where we set $T \,:\!=\, \theta(1)$ so that $s\in[0,1]$ . A closed-form expression for $P_{\psi}(T)$ can be found [Reference Lerche38, p. 28]. The exact boundary crossing probability in this case is $0.750\,502\,88\ldots$

Using (3) and (4), the approximate drift and diffusion coefficients of the weak second-order Itô–Taylor expansion for the OU process are given by ${\beta}_{n,k}(x) = - x + \tfrac12\Delta_{n}x$ , $\alpha^{1/2}_{n,k}(x) = 1 - \tfrac12\Delta_{n}$ . From the numerical results below, it appears that it is sufficient to use the weak second-order Itô–Taylor expansion of transition densities instead of the true transition density of the OU process for our Markov chain approximation to maintain a $O(n^{-2})$ convergence rate of the boundary crossing probabilities.

In the log–log plot in the left pane of Fig. 3, the crosses $\times$ show the absolute error of the Markov chain approximation without the Brownian bridge correction, while the bullets $\bullet$ show the error with the Brownian bridge correction. The markers $\bigtriangleup$ and $\circ$ show the absolute error of the Markov chain approximation using the exact transition density of the OU process and the Euler–Maruyama approximation instead of the transition density from the Itô–Taylor expansion, respectively. From this plot, we empirically observe that the convergence rate $\lvert \mathbb{P}(X\in G) - \mathbb{P}({\widetilde{X}_{n}} \in G_n)\rvert$ is of the order $O(n^{-2})$ , which is the same as the boundary approximation order of error $\lvert \mathbb{P}(X \in G) - \mathbb{P}(X \in G_n)\rvert = O(n^{-2})$ proved in [Reference Downes and Borovkov15] in the case of twice continuously differentiable boundaries. It appears that, at this level of accuracy, we might ignore the higher-order terms in the diffusion bridge crossing probability derived in [Reference Baldi and Caramellino3].

Figure 3. Approximation of non-crossing probabilities of the OU process with boundaries $g_{\psi}^{\pm}$ . The left pane shows a log–log plot of the absolute approximation error as a function of n. The upper and lower dashed lines correspond to $C_1n^{-1/2}$ , $C_2n^{-1}$ , and $C_3n^{-2}$ respectively, where $C_1$ , $C_2$ , and $C_3$ are fitted constants. The right pane depicts the time evolution of the taboo transition density using the Markov chain approximation ${\widetilde{X}_{n}}$ with $n = 20$ .

5.3. The Bessel process

To test our method on a non–Gaussian diffusion process, we chose the Bessel process, due to the availability of a closed-form expression for its transition density. Note also that the Cox–Ingersoll–Ross (CIR) process popular in mathematical finance is a suitably time–space-transformed Bessel process, and hence boundary crossing probabilities for the CIR process can be immediately obtained from those for the Bessel case.

For $\nu \geq 0$ , the Bessel process of order $\nu$ can be defined as the (strong) solution of the stochastic differential equation [Reference Borodin and Salminen7, p. 66]

\begin{equation*} \mathrm{d} X(t) = \frac{2\nu + 1}{2X(t)}\,\mathrm{d} t + \mathrm{d} W(t), \qquad t >0; \quad X(0) = x_0 >0.\end{equation*}

The approximate drift and diffusion coefficients of the weak second-order Itô–Taylor expansion for the Bessel process are given by

\begin{equation*} \beta_{n,k}(x) = \frac{2\nu + 1}{2 x} + \frac{1-4\nu^2}{8x^3}\Delta_n, \qquad \alpha^{1/2}_{n,k} (x) = 1 - \frac{2 \nu +1}{4x^2}\Delta_n.\end{equation*}

For numerical illustration purposes, we tested our algorithm in the case when $\nu = \frac12$ , $x_0 = 1$ , $g^{-}(t) = 0.7$ and $g^+(t) = +\infty$ , for which a closed-form expression is available for benchmarking [Reference Hamana and Matsumoto26, Theorem 2.2]). The results of the numerical calculations are shown in Fig. 4. The left pane presents a log–log plot of the absolute approximation error. The $\times$ markers show the Markov chain approximation error without the Brownian bridge correction. The markers $\circ$ , $\bullet$ , and $\Delta$ show the error of the Markov chain approximation with the Brownian bridge correction when using the Euler–Maruyama approximation, the Itô–Taylor approximation, and the exact transition density respectively. We can see that the observed convergence rate is the same, $O(n^{-2})$ , whether one uses the Itô–Taylor approximation or the exact transition densities, but the rate drops to $O(n^{-1})$ when using the Euler–Maruyama approximation. The right pane shows the time evolution of the taboo transition density using the Markov chain approximation ${\widetilde{X}_{n}}$ with $n = 20$ . An artificial absorbing boundary has been inserted at $x = 7$ to ensure that probability mass is not lost.

Figure 4. Approximation of the boundary crossing probabilities for the Bessel process with initial condition $X(0) = 1$ , starting above the boundary $g^-(t) = 0.7$ . The exact boundary crossing probability in this case is $0.534\,924\ldots$ The left pane shows a log–log plot of the absolute approximation error as a function of n. The dashed lines correspond to $C_1n^{-1/2}$ , $C_2n^{-1}$ , $C_3n^{-2}$ , and $C_4n^{-2}$ from top to bottom, where $C_1$ , $C_2$ , $C_3$ , and $C_4$ are fitted constants. The right pane depicts the time evolution of the taboo transition density using the Markov chain approximation ${\widetilde{X}_{n}}$ with $n = 20$ .

5.4. Comparison with other methods

In this subsection we will comment on the relative performance of our approach compared to the previously proposed ones in the special basic case of one-sided boundary g for the Brownian motion process.

Based on the results of our numerical experiments, one of the most efficient methods for numerical computation of the boundary crossing probability of the Brownian motion is the method of integral equations [Reference Peskir49]. The method proposed in [Reference Loader and Deely41, Reference Park and Schuurmann47] is based on numerically solving (using ‘midpoint’ quadratures) a Volterra integral equation of the first kind with the kernel $K (t,s): = \overline{\Phi}((g(t)-g(s))/\sqrt{t-s})$ , $0<s<t$ , where $\overline{\Phi}$ is the standard normal ‘tail’. Unfortunately, the convergence rate for that method was discussed in neither [Reference Park and Schuurmann47] nor [Reference Loader and Deely41]. If the above kernel K had no end-point singularity, it would follow from the proof of [Reference Linz40, Theorem 9.1] that the convergence rate is $O(n^{-2})$ , n being the number of steps on the uniform time grid used, since that rate is essentially determined by the quadrature error. In view of the results of [Reference Lyness and Ninham42, Section 7], the proof of the above-mentioned theorem from [Reference Linz40] indicates that, in the case of our kernel K, the convergence rate will be $O(n^{-3/2})$ . And indeed, our numerical experiments show that the convergence rate for that method is of that order and, moreover, that the error may have the form $C_0 n^{-3/2} + C_0 n^{- 2} + C_0 n^{-5/2} + \cdots$ . Hence, in our comparison, when comparing the results of different methods after applying Richardson’s extrapolation [Reference Sidi57, p. 27], we assumed that the error dependence on the small parameter $n^{-1}$ for the above-mentioned integral equation method is of that form.

In [Reference Ricciardi, Sacerdote and Sato52] the first-passage-time density of the boundary g was shown to satisfy a Volterra integral equation of the second type. Its solution was approximated using the trapezoidal rule in [Reference Buonocore, Nobile and Ricciardi11]; however, the convergence rates where not discussed in that paper. In [Reference Sacerdote and Tomassetti55], the authors obtained a series expansion for the solution to the integral equation using the method of successive approximations, with each term in the series containing an integral. Error analysis was done on the truncation of the series expansion. However, the integrals in the series expansion were obtained numerically using quadratures, and there was no mention of the convergence rate for the quadratures.

Note that both Volterra integral equations mentioned above are special cases of the spectrum of integral equations given in [Reference Peskir49, Theorem 6.1]. Our numerical experiments showed that the approximation from [Reference Park and Schuurmann47] performs better, and hence we chose it to be the benchmark for the comparison in this subsection.

In Fig. 5 we show a comparison between the Markov chain approximation and the integral equation method proposed in [Reference Loader and Deely41], which is a version of the approach from [Reference Park and Schuurmann47] (labeled ‘IE’). The label ‘MC (No BB)’ refers to the Markov chain approximation results obtained without using the Brownian bridge correction, whereas just ‘MC’ corresponds to the Markov chain method with that correction.

Figure 5. Log–log plot of the absolute approximation errors versus the computational time for our Markov chain approximation and the integral equation method applied to the standard Brownian motion crossing Daniels’ boundary before time $T=1$ . Labels ending with R $(p_1,p_2,\ldots,p_k)$ indicate the use of Richardson’s extrapolation applied repeatedly to the respective sequence A(h) with an error of the form $a_1 h^{p_1} + a_2 h^{p_2} + \cdots + a_kh^{p_k} + o(h^{p_k})$ , h being the appropriate small parameter in the scheme used.

We see from Fig. 5 that applying the Brownian bridge correction dramatically improves the efficiency of our scheme. Pre-Richardson’s extrapolation, the computational times for the integral equation method are generally lower, which is unsurprising since our method discretises both space and time variables, while the integral equation method only discretises time. However, our scheme becomes competitive once we apply Richardson’s extrapolation. Furthermore, and most importantly, our method works for general diffusion processes whereas the method of integral equations requires an explicit expression for the transition density.

Another important advantage of our approach is that it can easily be modified to calculate expressions that involve the space variable, because our scheme essentially approximates taboo transition probabilities.