1. Introduction
The interest in understanding the first passage time (FPT) can be traced back to the early 20th century [Reference Bachelier4], [Reference Schrödinger46]. Known also as the first hitting time, the FPT is defined as a random time at which a stochastic process visits a predefined state. The phenomenon of uncertainty in time is often observed in natural or social science. Therefore, for over a century, the FPT has been actively studied in economics, physics, biology, etc. [Reference Dassios and Lim14], [Reference Novikov, Frishling and Kordzakhia37], [Reference Redner43], [Reference Ricciardi, Crescenzo, Giorno and Nobile45].
Depending on various types of underlying processes and hitting boundaries, the FPT itself involves a large cluster of different research topics. We refer to [Reference Arbib3], [Reference Blake and Lindsey7], [Reference Novikov38] and [Reference Siegert51] for a non-conclusive review. Among these topics, especially in the area of mathematical finance and insurance, the single-side constant-barrier crossing problem is one of the most commonly studied; see, e.g., [Reference Baldi, Caramellino and Iovino5] and [Reference Dassios and Zhang17]. A general approach to solving such problems starts with finding the Laplace transform (LT) of the FPT density (FPTD). The LT usually comes from a unique solution to a second-order non-homogeneous ordinary differential equation (ODE) with Dirichlet-type boundary values [Reference Doob18], [Reference Ikeda and Watanabe26]. For many familiar diffusion processes, the LTs have been solved and are listed in [Reference Borodin and Salminen9]. However, these LTs usually are expressed in terms of special functions, and only a few of them have explicit inverse transforms. Therefore, many efforts have been made on the numerical inverse side. We refer to [Reference Abate and Whitt1] for more details. Alternatively, using the spectral theorem on linear operators [Reference Itô and McKean27], [Reference Kent30], [Reference Kent31], one can simplify the original LT. Under certain circumstances, closed-form FPTDs can be obtained through series representations [Reference Alili, Patie and Pedersen2], [Reference Linetsky34]. But one may find that the spectral decomposition approach has convergence issues for small t. In the present paper, our object is to apply perturbation theory and solve explicit FPTD approximations for general single-side level crossing problems.
Consider a filtered probability space 
$\left(\Omega,\mathcal{F},\left\{\mathcal{F}_t\right\}_{t\ge 0},\mathbb{P}\right)$
which satisfies the usual conditions and is generated by a standard Brownian motion. Let I be an open interval on 
$\mathbb{R}$
, and let 
$h(\cdot)$
be a real-valued continuous function defined on I. Our underlying process is from a class of stochastic differential equations (SDEs) which have at least weak solutions and are strong Markov:
\begin{equation}X_t=\epsilon h(X_t)dt + dW_t,\ X_0=x \in I.\end{equation}
In our setting, 
$\epsilon\in{\mathbb R}$
, and it should properly define 
${\left\{{X_t}\right\}}_{t\ge 0}$
on its domain. For convenience of deduction, we set the volatility to be constant. If a time-homogeneous diffusion coefficient 
$\sigma(x)$
is given, one may refer to [Reference Revuz and Yor44, Theorem 1.6] to retrieve an SDE in (1) by using a time-changed Brownian motion. Also, given a hitting level 
$a\in \mathbb{R}$
, we specify two types of boundaries on I, namely boundaries for upper and lower regions:
\begin{equation*}\partial I^u_a \coloneqq \left\{a,+\infty\right\},\ \partial I^l_a \coloneqq \left\{-\infty, a\right\}.\end{equation*}
As shorthand, we use 
$\partial I_a$
to represent a single-side boundary without labelling the direction. Suppressing x and a, we define the FPT of 
${\left\{{X_t}\right\}}_{t\ge 0}$
from x to a as
\begin{equation*}\tau \coloneqq \inf\left\{t> 0\,{:}\,X_t\in \partial I_a\right\}.\end{equation*}
Note that the Brownian filtration 
${\left\{{\mathcal{F}_t}\right\}}_{t\ge 0}$
is continuous. Therefore, according to [Reference Peskir and Shiryaev42], 
$\tau$
is well defined (i.e. regular at 
$\partial I_a$
). In addition, for 
$x\in I$
, it is guaranteed that 
$\mathbb{P}_x\left(\tau>0\right)=1$
.
For those FPTs which are almost surely (a.s.) finite, i.e. 
$\mathbb{P}_x\left(\tau<+\infty\right)=1$
, we are interested in finding their explicit distributions. Clearly, when 
$h(x)\equiv 0$
(a standard Brownian motion), the distribution of 
$\tau$
is given by an inverse Gaussian (or inverse gamma, equivalently) [Reference Borodin and Salminen9]. However, for most nontrivial drifts, there is no closed-form solution. An example is 
$h(x)=x$
, which corresponds to the Ornstein–Uhlenbeck (OU) process. In this case, the explicit density is only available when restricting to 
$a=0$
[Reference Göing-Jaeschke and Yor20].
In this paper, we apply the perturbation technique [Reference Holmes24] to solve Dirichlet-type boundary value problems (BVPs). By inverting the perturbed LTs from the frequency domain, where those LTs usually have much simpler forms, we then are able to derive closed-form asymptotic densities in the time domain. The main contribution of this paper is to provide a unified recursive framework for solving the single-side barrier hitting problem. By using the Green’s function representation and potential theory [Reference Peskir and Shiryaev42], we prove the convergence of the perturbation series and demonstrate the error of a truncated series, respectively. As illustrations, we show the perturbed FPTDs of the OU, the Bessel, the exponential-Shiryaev [Reference Dassios and Li13], and the hypergeometric diffusion processes [Reference Borodin8] in this paper, alongside their numerical comparisons with inverse Laplace transform algorithms [Reference Abate and Whitt1] and Monte Carlo simulations.
The rest of the paper is organised as follows: Section 2 introduces our main results; Section 3 demonstrates applications on various diffusion processes; in Section 4 we provide numerical comparisons; Section 5 concludes. Apart from where it is mentioned in the main body, all proofs can be found in [Reference Li33] or the appendix of this paper.
2. Main results
2.1. Perturbed Dirichlet problem
We follow our previous set-up. Further, let 
$C^2$
denote 
$C^2(I)$
. For any 
$f\in C^2$
, we also assume that the infinitesimal generator 
$\mathcal{A}\,f(x)$
of 
${\left\{{X_t}\right\}}_{t\ge 0}$
exists for all 
$x\in I$
:
\begin{equation*}\mathcal{A}\,f(x)=\epsilon h(x)f^{\prime}(x) +\mathcal{G}f(x),\end{equation*}
where 
$\mathcal{G}\cdot$
is the infinitesimal generator for a standard Brownian motion,
\begin{equation}\mathcal{G}f(x)=\frac{1}{2}f^{\prime\prime}(x).\end{equation}
Consider 
$\beta\in\mathbb{C}^{+}$
(i.e. 
$\beta \in \mathbb{C}$
with 
$Real(\beta)> 0$
), and define
\begin{equation}f(x,\beta) \coloneqq \mathbb{E}_x\left[e^{-\beta \tau}V\left(X_\tau\right)\right],\end{equation}
where 
$V\left(\cdot\right)$
is a finite function. The first step of our work is to find a proper BVP which is satisfied by 
$f(x,\beta)$
. To see this, note that 
${\left\{{X_t}\right\}}_{t\ge 0}$
is continuous over all stopping times; on the other hand, by our assumption 
${\left\{{X_t}\right\}}_{t\ge 0}$
is a strong Markov process. According to potential theory [Reference Peskir and Shiryaev42], 
$f(x,\beta)$
is the unique solution to the following Dirichlet problem:
\begin{equation}\mathcal{A}\,f(x)=\beta f(x),\ x\in I.\end{equation}
Moreover, the corresponding boundary conditions are given by
\begin{equation}f(\partial I_a)=\underline{V}.\end{equation}
In the notation above, 
$\underline{V} \coloneqq \left[V(a),V(\pm \infty)\right]^T$
is a vector of the boundary values depending on the direction of crossing. Refer to (3); by setting 
$V(a)=1$
and 
$V(\pm\infty)=0$
, we immediately find that the solution to the BVP (4) and (5) is the LT for the density of 
$\tau$
:
\begin{equation*}f(x,\beta)=\mathbb{E}_x\left[e^{-\beta \tau}\right].\end{equation*}
In the second step, we apply perturbations on 
$\epsilon$
and find perturbed BVPs accordingly. The perturbation approach is a common technique in solving asymptotics for complex systems. It has been successfully applied in quantum physics and mathematical finance [Reference Dassios and Wu16], [Reference Fouque, Papanicolaou, Sircar and Sølna19], [Reference Schrödinger47]. Traditionally, it is required that the perturbation parameter should be small. However, we will show later that this is not necessary in our case.
For brevity, we ignore the function arguments in following. By default all operations are with respect to x. Consider a sequence of 
$C^2$
functions 
$\left\{f_i\right\}_{i\ge 0}$
such that f can be expressed as
\begin{equation}f=\sum_{i=0}^\infty \epsilon^i f_i.\end{equation}
Substituting (6) into (4) yields
\begin{equation}\sum_{i=0}^{\infty}\epsilon^i\left(\epsilon h f_i^{\prime} +\mathcal{G}f_i\right)=\sum_{i=0}^\infty \epsilon^i \beta f_i,\ \forall x\in I.\end{equation}
Rearranging terms in (7), we further get
\begin{equation*}\mathcal{G}f_0-\beta f_0+\sum_{i=1}^{\infty}\epsilon^i\left(h f_{i-1}^{\prime}+\mathcal{G}f_i-\beta f_i\right)=0,\ \forall x\in I.\end{equation*}
Note that by extracting the 0th order and assigning proper boundary conditions we can obtain the BVP for the standard Brownian motion (where the LT inverse is already known). Higher-order BVPs can be solved via a recursive system which accumulates information from 
$f_0$
and the drift function h.
Denote the BVP with 
$i=0$
by the term o(1); by assigning the same boundary conditions as in the initial problem, we have
\begin{align*}o(1)\,:\,\mathcal{G}f_0 & =\beta f_0,\ x\in I,\\[3pt]
f_0(\partial I_a) & =[1,0]^T.\end{align*}
For 
$i\ge 1$
, we use the notation 
$o(\epsilon^i)$
and define
\begin{align*}o\left(\epsilon^i\right)\,:\, \mathcal{G}f_i & =\beta f_i-h\cdot f_{i-1}^{\prime},\ x\in I,\\[3pt]
f_i(\partial I_a) & =[0,0]^T.\end{align*}
In practice, it is not realistic to have infinite-order solutions. We are more interested in knowing if the series (6) converges, or, given a truncation order 
$N\in \mathbb{N}$
, what exactly the error terms are in the remaining higher-order ODEs. These two questions are answered, respectively, by the following two subsections.
2.2. Convergence of the perturbation series
Proposition 2.1. (Sufficient conditions for the convergence of the perturbation.) Let 
$\beta\in\mathbb{C}^{+}$
, and let h be the real-valued continuous drift function defined on I. If both h and 
$h^{\prime}$
are bounded on I and 
$\beta$
is large enough, then the perturbation series (6) is convergent.
Proof. Without loss of generality, we consider the domain 
$I=(0,+\infty)$
(this requirement is only for the simplicity of the presentation; depending on the hitting direction the domain could be chosen as 
$(a,+\infty)$
or 
$(\!-\infty, a)$
).
Define 
$G_\beta$
to be the Green’s operator of 
$\beta-\mathcal{G}$
(cf. Equation (2)) in I such that
\begin{equation*}G_\beta u(x) \coloneqq \int_0^{+\infty}G_\beta(x,y)u(y)dy,\end{equation*}
where 
$G_\beta(x,y)$
is the Green’s function of the linear operator 
$\beta-\mathcal{G}$
. By the definition of the Green’s function, we can show that 
$G_\beta(x,y)$
is of the following explicit form:
\begin{equation*}G_\beta(x,y)=\begin{cases}\sqrt{\frac{2}{\beta}}\sinh\big(\sqrt{2\beta}y\big)e^{-\sqrt{2\beta}x},\ y\le x,\\[9pt]
\sqrt{\frac{2}{\beta}}\sinh\big(\sqrt{2\beta}x\big)e^{-\sqrt{2\beta}y},\ y> x.\end{cases}\end{equation*}
We further introduce the multiplication operator 
$Hu(x) \coloneqq h(x)u(x)$
and the derivative operator 
$D u(x) \coloneqq u^{\prime}(x)$
. By considering the solution of the 
$o(\epsilon^i)$
-ODE and using the Green’s function, one can check that
\begin{equation}f_i=(G_\beta H D)^i f_0.\end{equation}
Substituting (8) into (6) further yields
\begin{equation}f=\sum_{i=0}^\infty(\epsilon G_\beta H D)^i f_0.\end{equation}
As a sufficient condition for (9) to converge, we only need to show that the norm of the operator 
$\epsilon G_\beta H D$
defined on a proper function space is less than one. To see this, consider the Banach space 
$L^\infty(I)$
and 
$u\in L^\infty(I)$
. By definition,
\begin{align}|| \epsilon G_\beta H D u||_\infty=& ||\epsilon\sqrt{\frac{2}{\beta}}\bigg(\int_0^x \sinh\big(\sqrt{2\beta}y\big)e^{-\sqrt{2\beta}x}h(y)u^{\prime}(y)dy\end{align}
\begin{align}\hspace*{42pt} &+ \int_x^\infty \sinh\big(\sqrt{2\beta}x\big)e^{-\sqrt{2\beta}y} h(y)u^{\prime}(y) dy\bigg)||_\infty. \hspace*{-25pt}\end{align}
Using integration by parts, the integral in the right-hand side of (10) becomes
\begin{align*} \hspace*{-45pt}\int_0^x \sinh\big(\sqrt{2\beta}y\big)e^{-\sqrt{2\beta}x}h(y)u^{\prime}(y)dy\hspace*{45pt}\end{align*}
\begin{align}= e^{-\sqrt{2\beta}x}\int_0^x (u(x)-u(z))\sqrt{2\beta}\cosh(\sqrt{2\beta}z)h(z)dz\end{align}
\begin{align} +e^{-\sqrt{2\beta}x}\int_0^x (u(x)-u(z))\sinh\big(\sqrt{2\beta}z\big)h^{\prime}(z)dz.\end{align}
Similarly, the integral in (11) becomes
\begin{align*}\hspace*{-45pt} \int_x^\infty \sinh\big(\sqrt{2\beta}x\big)e^{-\sqrt{2\beta}y} h(y)u^{\prime}(y) dy \hspace*{45pt}\end{align*}
\begin{align}= \sinh\big(\sqrt{2\beta}x\big)\int_x^\infty (u(z)-u(x))\sqrt{2\beta}e^{-\sqrt{2\beta}z} h(z) dz\end{align}
\begin{align} -\sinh\big(\sqrt{2\beta}x\big)\int_x^\infty (u(z)-u(x))e^{-\sqrt{2\beta}y} h^{\prime}(y) dy.\end{align}
We observe that (12) is bounded by
\begin{equation}\frac{1-e^{-2\sqrt{2\beta}x}}{2}(||u||_\infty-u(x))||h||_\infty,\end{equation}
(13) is bounded by
\begin{equation}\frac{1+e^{-2\sqrt{2\beta}x}-2e^{-\sqrt{2\beta}x}}{2\sqrt{2\beta}}(||u||_\infty-u(x))||h^{\prime}||_\infty,\end{equation}
(14) is bounded by
\begin{equation}\frac{1-e^{-2\sqrt{2\beta}x}}{2}(||u||_\infty-u(x))||h||_\infty,\end{equation}
and (15) is bounded by
\begin{align}\frac{1-e^{-2\sqrt{2\beta}x}}{2\sqrt{2\beta}}(||u||_\infty-u(x))||h^{\prime}||_\infty.\end{align}
Now (16) and (18) together are bounded by
\begin{equation*}\left(1-e^{-2\sqrt{2\beta}x}\right)(||u||_\infty-u(x))||h||_\infty\le ||u||_\infty||h||_\infty,\end{equation*}
and (17) and (19) together are bounded by
\begin{equation*}\frac{2-2e^{-\sqrt{2\beta}x}}{2\sqrt{2\beta}}(||u||_\infty-u(x))||h^{\prime}||_\infty\le\frac{1}{\sqrt{2\beta}} ||u||_\infty||h^{\prime}||_\infty.\end{equation*}
Hence,
\begin{equation*}|| \epsilon G_\beta H D u||_\infty\le \epsilon \left(\sqrt{\frac{2}{\beta}}||h||_\infty+\frac{1}{\beta}||h^{\prime}||_\infty\right)||u||_\infty,\end{equation*}
and as long as h and 
$h^{\prime}$
are bounded on I and 
$\beta$
is large enough, the norm of the operator is guaranteed to be less than one. This concludes our proof.
Remark 2.1. Using the operator representation, Equation (8) produces the solution to the 
$o(\epsilon^i)$
-ODE. Alternatively, the solution can also be written recursively as
\begin{equation*}f_i=f_0\left[\int_{0}^x 2e^{2\sqrt{2\beta} y} \left(\int_{0}^y h(z)k_i(z,\beta)e^{-2\sqrt{2\beta} z}dz + C_1\right)dy +C_2\right],\end{equation*}
where
\begin{equation*}k_i(z,\beta) \coloneqq \sqrt{2\beta} f_0^{-1}(z,\beta)f_{i-1}(z,\beta)-\partial_z\left(f_0^{-1}(z,\beta)f_{i-1}(z,\beta)\right).\end{equation*}
$C_1$
and 
$C_2$
are constants subject to 
$f_i(\partial I,\beta)=[0,0]^T$
. We refer the reader to [Reference Li33, Lemma 3.4.2, p. 38] for more details of the recursive representation.
2.3 Error function
In this subsection, we provide an error function for the perturbed density function with a truncation order 
$N\in \mathbb{N}$
.
Denote the Nth-order truncated LT by
\begin{equation*}f^{(N)} \coloneqq \sum_{i=0}^N \epsilon^i f_i.\end{equation*}
We further assume that the inverse LTs of f, 
$f^{(N)}$
, and 
$\partial_x f_N(x,\beta)$
exist, and these are denoted by 
$p_\tau(t)$
, 
$p^{(N)}_\tau(t)$
, and 
$\eta_N(x,t)$
, respectively. Define the difference (absolute error) between the original and the perturbed FPTDs by
\begin{equation}q_\tau(t) \coloneqq p_\tau(t)-p_\tau^{(N)}(t).\end{equation}
Then we have the following result.
Proposition 2.2. (Probabilistic representation of the truncation error.) For all 
$t\in (0,+\infty)$
and all 
$\beta\in\mathbb{C}^{+}$
, if
\begin{equation}\int_0^{+\infty}e^{-\beta t}\mathbb{E}_x\left[\int_0^{\tau\wedge t}\left|h(X_u)\eta_N\left(X_u,t-u\right)\right|du\right]dt <+\infty,\end{equation}
then
\begin{equation}q_\tau(t)=\epsilon^{N+1}\mathbb{E}_x\left[\int_0^{\tau\wedge t}h(X_u)\eta_N\left(X_u,t-u\right)du\right].\end{equation}
Proof. Please refer to the proof of Proposition 3.3.1, pp. 35–37 of [Reference Li33].
Remark 2.2. The error function (22) relies on the 
$L^1$
-condition given by (21). The applications we present later show that the 
$L^1$
-condition is easy to justify in practice. On the other hand, by referring to the proof in [Reference Li33], we can also derive an operator representation of (22). But the corresponding Green’s function is for the operator 
$\beta - \mathcal{A}$
, given which it is very difficult to find a unified sufficient condition satisfying (21).
3. Applications
In this section, we use five diffusion processes to illustrate the applications of our perturbation framework. In the examples below, some of the closed-form density functions may already have been found, but because of the complexity of the FPT problem, those densities may either be valid only for special cases (e.g. the OU FPTD in the special case 
$\theta = 0$
[Reference Göing-Jaeschke and Yor20]) or suffer computational efficiency issues (e.g. Bessel FPTDs in [Reference Hamana and Matsumoto21]). For these well-studied diffusion processes, in this section and the following numerical illustration section, we will show that the perturbation can provide accurate density asymptotics while maintaining a much faster computing speed. On the other hand, considering the fact that there are still many diffusion processes whose closed-form FPTDs have not been found yet, another more important purpose of this section is to demonstrate that the perturbation framework can be applied to a wide branch of diffusion processes.
In order to keep the paper concise, in the contents below we provide only the perturbed density functions, with comments wherever necessary. More detailed discussions and proofs of the results can be found in [Reference Li33, Chapters 4--6] or the appendix of this paper.
3.1. Ornstein–Uhlenbeck process
The OU process was first introduced to describe the velocity of a particle that follows a Brownian motion movement [Reference Uhlenbeck and Ornstein52]. Later the process appeared widely in neural science [Reference Lánský, Sacerdote and Tomassetti32], [Reference Wan and Tuckwell54] and mathematical finance [Reference Heston23], [Reference Jarrow and Turnbull28], [Reference Schwartz48], [Reference Vasicek53]. The FPT of the OU process has been extensively studied (see [Reference Patie40] for a brief review), and the LT of the FPT can be found in [Reference Borodin and Salminen9, Eq. 2.0.1, p. 542]. In our setting, the drift function of the OU process is given by
\begin{equation}h(x)=\theta -x,\ \theta\in\mathbb{R},\ x\in \mathbb{R}.\end{equation}
The OU process has a unique strong solution, and it is recurrent in 
$\mathbb{R}$
; see [Reference Uhlenbeck and Ornstein52] and [Reference Wang and Uhlenbeck55]. Our framework therefore can be applied. For 
$a=0$
and 
$x>0$
, we consider the hitting-from-above problem, i.e. 
$I=(0,+\infty)$
. (Note that the generalised case 
$a<x$
and the hitting-from-below problem (
$I=(\!-\infty, a)$
) can be retrieved by using an affine transformation and a reflecting process, respectively.)
Proposition 3.1. (Nth-order perturbed FPTD of the OU process.) Let 
$N\in {\mathbb N}$
. The Nth-order perturbed downward FPTD of the OU process is given by
\begin{equation*}{p}_{\tau}^{(N)}(t)=\frac{e^{-\frac{x^2}{4t}}}{\sqrt{2\pi}}\sum_{n=0}^{2N-1}h_nt^{\frac{n}{2}-1}D_{-n+1}\left(\frac{x}{\sqrt{t}}\right),\end{equation*}
where 
$D_\cdot(\cdot)$
is the parabolic cylinder function,
\begin{equation*}h_n \coloneqq \sum_{\substack{i,j,k;\\2i-j-k=n}}\epsilon^i c_k^{(i,j)}\theta^kx^j,\end{equation*}
and 
$\{c_k^{(i,j)}\}$
is a triple-indexed real sequence
\begin{equation*}\{c_k^{(i,j)}\,{:}\, 0\le i\le N,\ 0\le j \le 2i,\ 0\le k \le (2i-j)\wedge i\}\end{equation*}
whose recursion representation is given in Appendix A.
Proof. Please refer to [Reference Li33, Proposition 4.1.3, p. 46].
Remark 3.1. Recall that in Proposition 2.1, a sufficient condition for the perturbation to converge is that h should be bounded in I. However, referring to (23), apparently this is not the case. In fact, according to [Reference Li33, Proposition 4.1.4, p. 47], the right-tail asymptotic of the perturbed FPTD is 
$t^{N-\frac{3}{2}}$
. Therefore, when 
$N\ge 2$
the perturbed density diverges when 
$t\uparrow +\infty$
.
Remark 3.2. On the other hand, when 
$N=1$
, the right-tail asymptotic converges to 0 at the rate of 
$t^{-\frac{1}{2}}$
. In this case, even though h is not bounded, by using Proposition 2.2, we have shown in [Reference Li33, Proposition 4.2.2, p. 52] that the perturbed density converges in 
$t\in[0,T]$
for any 
$T>0$
at the rate of 
$O(\epsilon^2)$
. Moreover, the error representation in Equation (22) is valid and the error function 
$\eta_1(x,t)$
is given by [Reference Li33, Lemma 4.1.5, p. 49].
Returning to Proposition 3.1, in view of the close relation between the D-functions and the Hermite polynomials [Reference Lozier35, Sections 12.7.1–12.7.2], given 
$N=1$
, we can further write 
$p_\tau^{(N)}(t)$
explicitly in the following way:
\begin{equation}{p}_{\tau}^{(1)}(t)=\Big(1+\epsilon \Big(\frac{1}{2}x^2 - \theta x+\frac{1}{2} t\Big)\Big)\frac{x}{\sqrt{2\pi t^3}}e^{-\frac{x^2}{2t}}.\end{equation}
3.2. Bessel process with 
$0<n<2$
The Bessel process was introduced in [Reference McKean36] as the norm of an n-dimensional Brownian motion. Denoted by BES(n) (sometimes by 
$BES(\nu)$
with 
$\nu = \frac{n-2}{2}$
), the basic properties of the process have been discussed in [Reference Revuz and Yor44, Chapter XI]. In mathematical finance, the family of Bessel processes is closely related to models of short rates and stochastic volatilities [Reference Carr and Linetsky10], [Reference Chen and Scott11], [Reference Cox, Ingersoll and Ross12], [Reference Dassios, Qu and Lim15], [Reference Heston23]. Like that of the OU process, the FPT of the Bessel process has been extensively studied. One can find the LTs for 
$n\ge 2$
in [Reference Borodin and Salminen9]; in [Reference Hamana and Matsumoto21], the LTs of 
$0<n<2$
(and other general n) and some explicit FPTDs are given.
We consider the class of Bessel processes with orders 
$n=1+\epsilon$
where 
$-1<\epsilon<1$
. For 
$BES(1+\epsilon)$
, the h-function is specified by
\begin{equation}h(x)=\frac{1}{2x},\ x\in (0,+\infty).\end{equation}
Let 
$I=(a,+\infty)$
with 
$0<a<x$
. The initial LT is given by the ratio of modified Bessel functions of the second kind; see [Reference Kent29] (and also [Reference Hamana and Matsumoto21, Eq. (2.5)]). In fact, according to [Reference Hamana and Matsumoto21, Eq. (2), Theorem 2.2], the explicit FPTD of the downward hitting is given by the following (note that 
$\nu<0$
in our case):
\begin{align}p_\tau(t)=&\left(1- \Big(\frac{a}{x}\Big)^{\frac{\epsilon-1}{2}}\int_0^{+\infty}\frac{L_{\frac{1-\epsilon}{2},\frac{x}{a}}(y)}{y}e^{-\frac{y(x-a)}{a}}dy\right)\frac{x-a}{\sqrt{2\pi t^3}}e^{-\frac{(x-a)^2}{2t}}\nonumber\\[6pt]
& + \Big(\frac{a}{x}\Big)^{\frac{\epsilon-1}{2}}\int_0^t \frac{x-a}{\sqrt{2\pi s^3}}e^{-\frac{(x-a)^2}{2s}}\int_0^{+\infty}L_{\frac{1-\epsilon}{2},\frac{x}{a}}(y)e^{-\frac{y(x-a)\sqrt{t}}{a\sqrt{s}}}\frac{(x-a)}{2a\sqrt{s t}}dy dt,\end{align}
where
\begin{equation*}L_{\mu,c}(y)=\frac{\cos(\pi \mu)\left(I_{\mu}(cy)K_\mu(y)-I_\mu(y)K_\mu(cy)\right)}{K^2_\mu(y)+\pi^2I_\mu^2(y)+2\pi\sin(\pi\mu)K_\mu(y)I_\mu(y)},\end{equation*}
and 
$I_\cdot(\cdot)$
, 
$K_\cdot(\cdot)$
are the modified Bessel functions.
From (26) above we can see that the density involves a convolution of Bessel functionals. Therefore, the computing speed might be limited in practice (assuming there is no Monte Carlo (M.C.) simulation/importance sampling or parallel computation used).
As an alternative, we provide the perturbed FPTD in the next proposition. Note that, in the case 
$I=(a,+\infty)$
, the drift function in (25) and its derivative are bounded on I. According to Proposition 2.1, the perturbation series converges. But for simplicity of calculation, the result in the following proposition is only calculated up to the first order. (Actually, our numerical tests show that the first-order perturbation is accurate enough in providing density asymptotics.) Higher-order results are left to future work.
Proposition 3.2. (First-order perturbed FPTD of 
$BES(1+\epsilon)$
.) The first-order perturbed downward FPTD of 
$BES(1+\epsilon)$
is given by
\begin{align}
p_\tau^{(1)}(t) & = \left(1+ \frac{\epsilon}{2}\ln\Big( \frac{a}{x}\Big)\right)\frac{x-a}{\sqrt{2\pi t^3}}e^{-\frac{(x-a)^2}{2t}}\nonumber\\[8pt]
&\quad +\frac{\epsilon(x-a)}{\sqrt{2\pi t}}\int_{(x-a)^2}^\infty \frac{1}{2t(\sqrt{y}-x+3a)(\sqrt{y}+x+a)}e^{-\frac{y}{2t}}dy.\end{align}
Proof. Please refer to the proof of Proposition 5.1.3, pp. 72–73 of [Reference Li33].
Remark 3.3. The perturbed FPTD (27) maintains the inverse gamma part as in the original density function (26). From the proposition above, we can clearly see that the first-order truncation is indeed a simplification of the higher orders, which were originally represented using the Bessel functionals. In terms of the truncation error, we have shown in [Reference Li33, Proposition 5.2.2, p. 77] that 
$p_\tau^{(1)}$
converges at the rate of 
$O(\epsilon^2)$
, and the error function 
$\eta_1(x,t)$
can be found therein.
3.3. Exponential-Shiryaev process
In this subsection, we consider a newly introduced three-parameter diffusion process:
\begin{equation}dX_t=\epsilon \big(e^{-2\alpha X_t}-c\big)dt+dW_t,\ X_0=x,\end{equation}
with 
$\epsilon\ge 0$
, 
$\alpha\ge 0$
, 
$0\le c\le 1$
, and 
$x\in{\mathbb R}$
. This process was discovered by the authors in the course of this research. The motivation for studying such a process is due to its special drift function properties. Note that the exponential functions are closed under differentiation and integration; therefore, in applying the perturbation mechanism, we expected the perturbed density to have a neat mathematical form. We have conducted a fundamental analysis of the process and realised that the sample path of 
${\left\{{X_t}\right\}}_{t\ge 0}$
can be used to describe the log-price of economic bubbles (cf. [Reference Dassios and Li13] and [Reference Li33, Chapter 6, p. 87]).
The Shiryaev process [Reference Shiryaev49], [Reference Shiryaev50] refers to the following SDE, which was derived by A. N. Shiryaev in the context of sequential analysis:
\begin{equation*}dZ_t = (1+\mu Z_t)dt + \sigma dW_t,\ \mu\in{\mathbb R},\ \sigma>0.\end{equation*}
Later, in the paper [Reference Peskir41], G. Peskir named this process the Shiryaev process and analysed its transition density. The transition density of the Shiryaev process is linked to the Hartman–Watson distribution [Reference Hartman and Watson22], which, moreover, is involved in the Asian-option pricing problem. The reason for calling the SDE (28) the exponential-Shiryaev process is that the exponential transform of 
${\left\{{X_t}\right\}}_{t\ge 0}$
,
\begin{equation*}Y_t \coloneqq e^{2\alpha X_t},\end{equation*}
is a scaled version of 
${\left\{{Z_t}\right\}}_{t\ge 0}$
; i.e.
\begin{equation}Y_t=2\alpha \epsilon Z_t\end{equation}
with 
$\mu=\frac{\alpha}{\epsilon}-c$
and 
$\sigma = \frac{1}{\epsilon}$
.
Our present paper focuses on illustrating the perturbed FPTD of the exponential-Shiryaev process (SDE (28)). (In fact, the FPTD of the Shiryaev process can then be found based on (29) and using the FPTD of 
${\left\{{X_t}\right\}}_{t\ge 0}$
.) We have shown in [Reference Li33, Chapter 6, p. 87] that 
${\left\{{X_t}\right\}}_{t\ge 0}$
has a unique strong solution, it is strong Markov, and its FPT is finite a.s. A unique solution of the BVP therefore exists. For a hitting level 
$a\in{\mathbb R}$
, the original LTs are given by
\begin{equation}f(x,\beta)=\begin{cases}e^{\lambda(x)-\lambda(a)}M\left(m,n,\psi(x)\right)/M\left(m,n,\psi(a)\right),\ x>a,\\[5pt]
e^{\lambda(x)-\lambda(a)} U(m,n,\psi(x))/U(m,n,{\hat\psi(a)}),\ x<a,\end{cases}\end{equation}
where 
$M(\cdot, \cdot,\cdot)$
and 
$U(\cdot, \cdot,\cdot)$
are the Kummer’s functions and
\begin{equation*}\begin{cases}m \coloneqq (\sqrt { {c}^{2}{\epsilon}^{2}+2 \beta}- \epsilon c)/(2\alpha),\\[4pt]
n \coloneqq (\sqrt {{c}^{2}{\epsilon}^{2}+2\beta}+ \alpha)/\alpha,\\[4pt]
\psi(x) \coloneqq \epsilon {{ e}^{-2 \alpha x}}/{\alpha},\\[4pt]
\lambda(x) \coloneqq -2\alpha m x.\end{cases}\end{equation*}
(See [Reference Li33, Section 6.4.1, pp. 105–107].)
For a quick illustration of the perturbation, we consider the hitting-from-above problem. Without loss of generality, we let 
$a=0$
(the case of general 
$a\in{\mathbb R}$
can be solved by using an affine transformation) and 
$I=(0,+\infty)$
. According to Proposition 2.1, the drift function
\begin{equation*}h(x)=e^{-2\alpha x}-c,\ x\in I,\end{equation*}
and its derivative are bounded, and therefore the perturbation series converges. Similarly as in the previous subsection, in the proposition below we provide the first-order truncated result. In the special case 
$c=0$
, a recursive solution of the higher-order ODEs and the inverse transforms can be found in [Reference Li33, Appendix 6.A, pp. 124–127].
Proposition 3.3. (First-order perturbed FPTD of the exponential-Shiryaev process.) The first-order perturbed downward FPTD of 
${\left\{{X_t}\right\}}_{t\ge 0 }$
is given by
\begin{align}
p_\tau^{(1)}(t) & = \bigg(1+\epsilon \bigg(cx+ \frac{\left(1-e^{-2\alpha x}\right)(\alpha t-x)}{2\alpha x}\bigg)\bigg)\frac{x}{\sqrt{2\pi t^3}}e^{-\frac{x^2}{2t}}\nonumber\\[5pt]
&\quad-\epsilon \frac{\alpha}{4}(1-e^{-2\alpha x}){{e}^{\alpha x ( \frac{\alpha t}{2x}+1 ) }}{\textit{Erfc}}\bigg(\frac{x}{\sqrt{2t}}+\alpha \sqrt {\frac{t}{2}}\bigg),\end{align}
where 
$\text{Erfc}(\cdot)$
is the complementary error function.
Proof. Please refer to the proof of Proposition 6.4.1, pp. 108–109 in [Reference Li33].
3.4. Hyperbolic Ornstein–Uhlenbeck process
In the following two subsections, we consider two special processes from the hypergeometric diffusion class [Reference Borodin8]. The name ‘hypergeometric diffusion’ comes from the fact that the LTs of their FPTs to a constant level are connected to the hypergeometric functions. The hypergeometric diffusion is defined as
\begin{equation}dZ_t = \Big(\Big(\nu+\frac{1}{2}\Big)c\coth(cZ_t)-\frac{\rho}{c}\tanh(cZ_t)\Big)dt + dW_t,\ Z_0>0,\end{equation}
with 
$\nu\ge -\frac{1}{2},\ \rho\ge 0,\ c\ge 0$
(see [Reference Borodin8, Eq. (4.1)]).
Note that the hypergeometric diffusion can be regarded as a generalisation of other familiar stochastic processes. For example, if we let 
$c\downarrow 0$
with 
$\nu \neq -\frac{1}{2}$
and 
$\rho\neq 0$
, we then have the radial OU process. If we keep 
$c\downarrow 0$
and 
$\nu\neq-\frac{1}{2}$
but let 
$\rho = 0$
, we get the Bessel process. Similarly, if we let 
$c\downarrow 0$
, 
$\rho\neq0$
, but 
$\nu = -\frac{1}{2}$
, we get the OU process.
In this subsection, we focus on the case of 
$c> 0$
and 
$\nu=-\frac{1}{2}$
. Let 
$\rho = \epsilon>0$
and further consider 
$\theta\in{\mathbb R}$
; we define the hyperbolic OU process as follows:
\begin{equation}dX_t=\epsilon\Big(\theta - \frac{\tanh(cX_t)}{c}\Big)dt + dW_t,\ X_0\in{\mathbb R}.\end{equation}
The drift function of 
${\left\{{X_t}\right\}}_{t\ge 0}$
is correspondingly given by
\begin{equation}h(x)=\theta - \frac{\tanh(cx)}{c},\ x\in{\mathbb R}.\end{equation}
Again, for simplicity of presentation, we demonstrate only the hitting-from-above case. The hitting-from-below case can be retrieved by using the symmetry of 
$\tanh(\cdot)$
. Let 
$a>0$
be the hitting level, and let 
$I=(a,+\infty)$
. (The case 
$a<0$
can be solved as well; however, it involves a different series representation in the perturbed density. For simplicity we skip this case, but the result can be provided upon request.) After tedious calculations, we have the following LT of the true FPTD (solution to the original BVP):
\begin{equation}f(x,\beta)=\frac{L(x,\beta)}{L(a,\beta)},\end{equation}
where
\begin{align*}L(x,\beta) \coloneqq &\frac{ \left(\tanh \left( cx \right) -1 \right) ^{\frac{1}{2}(A_1+A_2+\epsilon(1+\theta c)/c^2)}}{\left( \tanh \left( cx \right) +1 \right) ^{{\frac {\epsilon\,\theta\,c+{c}^{2}+\epsilon}{2{c}^{2}}}}}\\
& \big( {\mbox{$_2\textit{F}_1$}(A_1,B_1;\,C_1;\,1)}{\mbox{$_2\textit{F}_1$}(A_2,B_2;\,C_2;\,Z(x))} \left( 2\,\tanh \left( cx \right) +2 \right) ^{\frac{1}{2}C_2}{4}^{C_1-1}\\
&-{\mbox{$_2\textit{F}_1$}(A_2,B_2;\,C_2;\,1)}{\mbox{$_2\textit{F}_1$}(A_1,B_1;\,C_1;\,Z(x))} \left( 2\,\tanh \left( cx \right) +2 \right) ^{\frac{1}{2}C_1} \big),\end{align*}
with
\begin{equation*}\begin{cases}A_1 \coloneqq {\frac {\sqrt {\left( c\theta-1 \right) ^{2}{\epsilon}^{2}+2\,\beta\,{c}^{2}}-2\,\epsilon+\sqrt {\left( c\theta+1 \right) ^{2}{\epsilon}^{2}+2\,\beta\,{c}^{2}}}{2{c}^{2}}},\\[4pt]
A_2 \coloneqq {\frac {\sqrt {\left( c\theta-1 \right) ^{2}{\epsilon}^{2}+2\,\beta\,{c}^{2}}-2\,\epsilon-\sqrt {\left( c\theta+1 \right) ^{2}{\epsilon}^{2}+2\,\beta\,{c}^{2}}}{2{c}^{2}}},\\[4pt]
B_1 \coloneqq 1+\frac{2\epsilon}{c^2}+A_1,\\[4pt]
B_2 \coloneqq 1+\frac{2\epsilon}{c^2}+A_2,\\[4pt]
C_1 \coloneqq {1+\frac {\sqrt {\left( c\theta+1 \right) ^{2}{\epsilon}^{2}+2\,\beta\,{c}^{2}}}{{c}^{2}}},\\[4pt]
C_2 \coloneqq 2-C_1,\\Z(x) \coloneqq \frac{1}{2}\tanh \left( cx \right) +\frac{1}{2},\end{cases}\end{equation*}
and with 
$ {}_{2}\textrm{F}_{1}(\cdot,\cdot;\cdot;\cdot.)$
denoting the Gauss hypergeometric function.
The LT in (35) is different from the one shown in [Reference Borodin8, Section 6]. The latter, in fact, only gives the LT of the reflecting hyperbolic OU process in the special case 
$\theta = 0$
. From the equations above, we can see that the LT of the hyperbolic OU FPTD is already complicated enough. It can be imagined that finding the explicit inverse would be even more difficult or even impossible. Therefore, we provide the first-order perturbed FPTD below as an asymptotic to the original density (note that the drift function in (34) and its derivative are bounded in 
${\mathbb R}$
, and according to Proposition 2.1 the perturbed series converges).
Proposition 3.4. (First-order perturbed FPTD of the hyperbolic OU process.) The first-order perturbed downward FPTD of the hyperbolic OU process is given by
\begin{equation}p^{(1)}_\tau(t)=\Big(1+\epsilon(x-a)\Big(\frac{1}{c}-\theta\Big)\Big)\frac{x-a}{\sqrt{2\pi t^3}}e^{-\frac{(x-a)^2}{2t}}-\frac{\epsilon}{c}T(t),\end{equation}
where
\begin{align*}T(t) &\coloneqq (A(x)-A(a))\sqrt{\frac{2}{\pi}}t^{-\frac{5}{2}}((x-a)^2-t)e^{-\frac{(x-a)^2}{2t}}\\[3pt]
&\quad +\sum_{n=1}^\infty \Bigg[(B(x,n)-B(a,n))\bigg(\frac{1}{\sqrt{2\pi}}\Big(\frac{(x-a)t^{-\frac{3}{2}}}{2}-cnt^{-\frac{1}{2}}\Big)e^{-\frac{(x-a)^2}{2t}}\\[3pt]
&\quad +\frac{c^2n^2}{2}e^{cn(x-a)+\frac{c^2n^2t}{2}}{Erfc}\left(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\right)\bigg)\\[3pt]
&\quad -(C(x,n)-C(a,n))\bigg(\sqrt{\frac{2}{\pi}}t^{-\frac{3}{2}}\left(\frac{(x-a)^2}{t}-1-cn(x-a)+c^2n^2t\right)e^{-\frac{(x-a)^2}{2t}}\\[3pt]
&\quad -c^3n^3e^{cn(x-a)+\frac{c^2n^2t}{2}}{Erfc}\left(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\right)\bigg)\Bigg],\end{align*}
given that
\begin{equation*}\begin{cases}A(x) \coloneqq {\frac {2{{\rm e}^{cx}}\arctan ( {{\rm e}^{cx}} ) }{{c}^{2}}}-{\frac {\ln ( ( {{\rm e}^{cx}} ) ^{2}+1 ) }{{c}^{2}}}-\frac{ \pi e^{cx}}{c^2}+\frac{2x}{c},\\[3pt]
B(x,n) \coloneqq (\!-1)^{n} \frac{e^{-2ncx} }{cn(2n+1)},\\[3pt]
C(x,n) \coloneqq (\!-1)^{n} \frac{e^{-2ncx} }{c^2n(2n+1)}.\end{cases}\end{equation*}
Proof. Please refer to Appendix B.
Remark 3.4. Referring to the proof in Appendix B, we see that the series representation of the perturbed density is from the expansion of the 
$\arctan$
function (Equations (B4) and (B5)). In fact, the series converges very fast, especially if we notice that
\begin{align*}e^{cn(x-a)+\frac{c^2n^2t}{2}}\text{Erfc}\left(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\right)\sim &e^{cn(x-a)+\frac{c^2n^2t}{2}}\frac{e^{-\Big(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\Big)^2}}{\Big(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\Big)\sqrt{\pi}}\\
=&\frac{e^{-\frac{(x-a)^2}{2t}}}{\Big(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\Big)\sqrt{\pi}},\ n\uparrow +\infty\end{align*}
(see [Reference Olver39, 7.1.23, p. 298]).
3.5 Hyperbolic Bessel process
Recall (32). By letting 
$\nu=-\frac{1}{2}+\epsilon$
with 
$\epsilon>0$
we have the hyperbolic Bessel process:
\begin{equation}dX_t=\epsilon c\coth(cX_t)dt + dW_t,\ X_0>0.\end{equation}
From [Reference Borodin and Salminen9, No. 33, p. 70], it can be seen that the hyperbolic Bessel process is used in the path decomposition of the drifted Brownian motion (with a nonnegative drift) conditioning on not hitting 0. More detailed discussions and transition densities of the hyperbolic Bessel process can be found in [Reference Borodin8, Section 5].
In this subsection, we consider the hitting-from-above case. Let 
$a>0$
and 
$I=(a,+\infty)$
. The h-function of 
${\left\{{X_t}\right\}}_{t\ge 0}$
is given by
\begin{equation}h(x)=c\coth(cx),\ x>0,\end{equation}
which, together with its derivative, is bounded in I. Again, by Proposition 2.1 we know that the perturbation series converges.
The original LT of the downward FPTD is given by
\begin{equation}f(x,\beta)=\frac{\cosh^{2\alpha}(ca){\mbox{$_2\text{F}_1$}(\alpha,\alpha+\frac{1}{2};\,2\alpha-\epsilon+1;\,1-\tanh^2(cx))}}{\cosh^{2\alpha}(cx){\mbox{$_2\text{F}_1$}(\alpha,\alpha+\frac{1}{2};\,2\alpha-\epsilon+1;\,1-\tanh^2(ca))}},\end{equation}
where
\begin{equation*}\alpha = \frac{\sqrt{\epsilon^2+\frac{2\beta}{c^2}}+\epsilon}{2}.\end{equation*}
(See [Reference Borodin8, Eqs. (1.3), (4.3), (4.4)].) As with the hyperbolic OU FPTD, in view of the complex structure of the LT, it seems to be not an easy task to find the explicit inverse of (39). In the next proposition, as an alternative, we provide the first-order perturbed FPTD.
Proposition 3.5. (First-order perturbed FPTD of the hyperbolic Bessel process.) The first-order perturbed downward FPTD of the hyperbolic Bessel process is given by
\begin{equation}p^{(1)}_\tau(t)=\left(1-\epsilon c(x-a)\right)\frac{x-a}{\sqrt{2\pi t^3}}e^{-\frac{(x-a)^2}{2t}}+\epsilon \tilde{T}(t),\end{equation}
where
\begin{align*}\tilde{T}(t) \coloneqq\, &(\tilde{A}(x)-\tilde{A}(a))\sqrt{\frac{2}{\pi}}t^{-\frac{5}{2}}\left((x-a)^2-t\right)e^{-\frac{(x-a)^2}{2t}}\\[5pt]
&+\sum_{n=1}^\infty \Bigg[(\tilde{B}(x,n)-\tilde{B}(a,n))\bigg(\frac{1}{\sqrt{2\pi}}\bigg(\frac{(x-a)t^{-\frac{3}{2}}}{2}-cnt^{-\frac{1}{2}}\bigg)e^{-\frac{(x-a)^2}{2t}}\\[5pt]
&+\frac{c^2n^2}{2}e^{cn(x-a)+\frac{c^2n^2t}{2}}\text{Erfc}\left(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\right)\bigg)\\[5pt]
&-(\tilde{C}(x,n)-\tilde{C}(a,n))\bigg( \sqrt{\frac{2}{\pi}}t^{-\frac{3}{2}}\left(\frac{(x-a)^2}{t}-1-cn(x-a)+c^2n^2t\right)e^{-\frac{(x-a)^2}{2t}}\\[5pt]
&-c^3n^3e^{cn(x-a)+\frac{c^2n^2t}{2}}\text{Erfc}\left(\frac{(x-a)}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\right)\bigg)\Bigg],\end{align*}
and
\begin{equation*}\begin{cases}\tilde{A}(x) \coloneqq -\frac{1-e^{cx}}{c}\ln(1-e^{-cx})-\frac{1+e^{cx}}{c}\ln(1+e^{-cx}),\\[5pt]
\tilde{B}(x,n) \coloneqq \frac{e^{-2ncx} }{n(2n+1)},\\[5pt]
\tilde{C}(x,n) \coloneqq \frac{e^{-2ncx} }{nc(2n+1)}.\end{cases}\end{equation*}
Proof. Please refer to Appendix B.
Remark 3.5. The two T-functions in Propositions 3.4 and 3.5 are very similar (the convergence in Proposition 3.5 is guaranteed by Remark 3.4), apart from the fact that the series functions 
$A,\ B,\ C$
and 
$\tilde A,\ \tilde B,\ \tilde C$
are different. In fact, the series functions are representations of integrals involving the 
$\ln$
(cf. Equations (B12) and (B13)) and 
$\arctan$
functions, respectively, while the inverse LTs in the two problems are the same. Instead of using the series representation, one can also refer to the proof in Appendix B to derive an integral form. We leave this to the reader.
Remark 3.6. Propositions 3.4 and 3.5 also provide building blocks for the perturbed FPTD of the hyperbolic radial OU process (Hyp-ROU). If we define 
$h_1$
and 
$h_2$
to be the drift functions (cf. Equations (34) and (37)) of the hyperbolic OU (Hyp-OU) and the hyperbolic Bessel (Hyp-Bes) processes, respectively, then the h-function of the Hyp-ROU is given by 
$h=h_1+h_2$
. Refer to the recursion formula in Remark 2.1; the first-order perturbed LT of the Hyp-ROU is
\begin{equation*}f_1^{Hyp-ROU}(x,\beta)=f_1^{Hyp-OU}(x,\beta)+f_1^{Hyp-Bes}(x,\beta),\end{equation*}
where 
$f_1^{Hyp-OU},\ f_1^{Hyp-Bes}$
are the first-order perturbed LTs of the Hyp-OU and the Hyp-Bes (cf. Equations (B9) and (B14)). Then, after a reorganisation of the proofs in Appendix B, the first-order perturbed FPTD of the Hyp-ROU is given by
\begin{equation*}p_\tau^{(1),\ Hyp-OU}=\Big(1+\epsilon(x-a)\Big(\frac{1}{c}-\theta-c\Big)\Big)\frac{x-a}{\sqrt{2\pi t^3}}e^{-\frac{(x-a)^2}{2t}}+\epsilon\Big( \tilde{T}(t)-\frac{T(t)}{c}\Big).\end{equation*}
A similar result can be obtained for the radial OU process (ROU; cf. Sections 3.1 and 3.2):
\begin{align*}{p}_{\tau}^{(1),\ ROU}(t)& = \Big(1+\epsilon \Big(\frac{1}{2}(x-a)^2 - \theta (x-a)+\frac{1}{2} t+\frac{1}{2}\ln\Big( \frac{a}{x}\Big)\Big)\Big)\frac{x-a}{\sqrt{2\pi t^3}}e^{-\frac{(x-a)^2}{2t}}\\
& \quad +\frac{\epsilon(x-a)}{\sqrt{2\pi t}}\int_{(x-a)^2}^\infty \frac{1}{2t(\sqrt{y}-x+3a)(\sqrt{y}+x+a)}e^{-\frac{y}{2t}}dy.\end{align*}
4. Numerical examples
In this section, we demonstrate six numerical examples (Figure 1) comparing the first-order perturbed FPTDs (cf. Propositions 3.1–3.5) to other well-studied methods. The perturbation parameter in all examples is chosen to be 
$\epsilon=0.1$
, and other parameters of different diffusions are listed in captions under Figures 1a–1f. Computing times of different methods are reported in Table 1, and deviation statistics for each method compared to the true values are reported in Table 2.
First-order perturbed FPTDs and comparisons with 
$\epsilon=0.1$
.
Among these six figures, two (Figures 1a and 1b) illustrate the FPTDs of the OU process. In Figure 1a, we consider the special case 
$\theta=a$
, and the explicit density formula is given in [Reference Göing-Jaeschke and Yor20]; the comparison is made among the perturbation/inverse LT (ILT; cf. the Talbot approach in [Reference Abate and Whitt1]), explicit density, spectral decomposition, and three-dimensional Brownian bridge simulation (BB-simulation) approaches (cf. [Reference Alili, Patie and Pedersen2] and [Reference Ichiba and Kardaras25], respectively). Figure 1b presents a general scenario where 
$\theta\neq a$
. Note that in this case no closed-form FPTD has yet been found.
In Figure 1c, the Bessel process with a non-integer order 
$n=1.1$
is studied. Unlike in other examples, the starting point and the hitting level are chosen as 
$x=0.7$
and 
$a=0.1$
. The reason for not selecting large levels (e.g. 
$x=2$
and 
$a=1.2$
) is that the Bessel process behaves very much like the Brownian motion when 
$X_t>>0$
. In this example, we include the explicit density function (cf. Equation (26)) as the benchmark for the perturbation and the Talbot inverse algorithm.
Figures 1d–1f are similar in their presentation formats. We consider respectively the FPTDs of the exponential-Shiryaev, the hyperbolic-OU, and the hyperbolic-Bessel processes. Apart from providing their perturbed and ILT FPTDs, we also include the brute-force Monte Carlo simulation (M.C.) and the standard Brownian motion (BM) FPTDs. (The ILT in Figure 1d is from the Talbot algorithm, while computations in Figures 1d and 1f are from the Gaver–Stehfest algorithm [Reference Abate and Whitt1]. The change in these two hypergeometric diffusions is due to the restriction of Python in computing hypergeometric functions with complex numbers.) Note that the black dashed curves in the last three graphs are the same, since the BM FPTD is only determined by the starting point and the hitting level.
The results from Figure 1 show that the first-order perturbation works well for all diffusion processes discussed above. Also, Figures 1d–1f indicate that, even for small 
$\epsilon$
, one cannot use the BM FPTD to obtain the desired approximation to the FPTDs of structurally different diffusions (although we do not provide the graphs for the OU and the Bessel processes, this is also the case there); the M.C. approach cannot give accurate estimates either. On the other hand, by checking the computing times in Table 1 (note that the times for Spectral* and BB-Sim* are scaling-up estimates from the 10 density points in Figure 1b), one may realise that the perturbation possesses the fastest computing speed in almost all scenarios (except that it is slower than the ILT in the Bessel FPTD).
In Table 2, we further compare each method’s deviation from the true values (explicit solutions wherever applicable). All statistics are calculated in an absolute manner, but note that for the spectral decomposition and the BB-simulation, there are only 10 density points (as opposed to 100 density points for the other methods). The table shows that the ILT produces more accurate results than our perturbation method in both processes. The average deviation shows that the errors of the perturbation are at the magnitude of 10e-4, although, at the total variation level, the perturbation generates about 3–5% error. However, as the purpose of the current paper is to produce the density curve instead of the probability curve, we would not consider the total deviation here to be a major benchmark. (In principle, the probability curve can be retrieved using the perturbation on different BVPs. We will leave this to future work.) On the other hand, as can be found in [Reference Alili, Patie and Pedersen2], the spectral decomposition has a convergence issue when t is small, and this explains the large deviations observed in Table 2. The BB-simulation produces the smallest error among the four methods, but by referring to Table 1 we see that this is at the cost of an extremely slow computation speed.
Computing times (seconds) of density curves in Figure 1.
Statistics of absolute deviations in Figure 1.
Another potential question is about the error behaviour on different model parameters. In [Reference Li33, Chapters 4–6], we have demonstrated perturbation error functions for different diffusion processes. But unfortunately, a uniform error bound has not been found yet. In general, by Proposition 2.2, we know for sure that the smaller the 
$\epsilon$
, the smaller the error. Apart from tuning 
$\epsilon$
, based on observations we have also found that smaller t and smaller ratios 
$x/a$
generate smaller errors. Taking the benchmark OU process (benchmark cases where 
$\theta = a=1.2$
) as an example, if we choose larger 
$\epsilon$
, smaller ratio 
$x/a$
, and smaller t (
$\epsilon=1.2,\ x=1.3,\ t\le 1$
), the average deviation reported in Table 2 becomes 7.00e-3; and if we choose smaller 
$\epsilon$
, larger ratio 
$x/a$
, and larger t (
$\epsilon=0.3,\ x=9,\ t\le 5$
), the average deviation increases to 3.26e-2. Apparently there is a complicated relationship among 
$\epsilon$
, x, 
$\theta$
, and a. But to understand it better, more serious study needs to be conducted. We leave this to future work.
As a last remark in this section, we highlight a few findings regarding higher-order OU perturbations. Recall that in Remark 3.1, we mentioned that the higher-order OU perturbations diverge on the right tail with 
$t^{N-\frac{3}{2}}$
. But this does not necessarily mean that the higher-order densities diverge for large t (at least not large enough for the divergence to appear); in fact, depending on the choices of 
$\epsilon$
and other parameters, higher-order densities may also produce more accurate results. Using the parameters 
$x=2,\ \theta=a=1.2$
again, in Table 3 we report average deviations and large t deviations (
$t=10$
years) for perturbed density curves compared to the true values, given 
$\epsilon\in{\left\{{0.1,0.3}\right\}}$
. The time range is chosen to be 
$t\in[0,10]$
. Results show that for 
$\epsilon=0.1$
, the perturbation converges up to 
$t=10$
, but for 
$\epsilon=0.3$
it slowly diverges.
OU benchmark deviations with different perturbation orders.
5. Conclusion
In this paper, we provide a systematic approach for solving the closed-form FPTD asymptotics of diffusion processes. We show a sufficient condition for the convergence of the perturbation series and derive a probabilistic representation for the error. The closed-form solution resulting from the perturbation not only increases computational efficiency, but also provides analytical tractability in understanding the FPTDs at extreme times. Using the framework we demonstrate valid approximations to FPTDs of various diffusion processes. Potential applications of this paper could be found in survival analysis, financial mathematics, and many other fields. Further work could be done in exploring the perturbed FPTDs of other diffusion processes, analysing error behaviours with different model parameters, and finding FPTDs via simulations using the probabilistic representation in Proposition 2.2.
Appendix A. Recursive structure in the OU perturbation
Result A.1. (Decomposition structure I.) For 
$i=1$
and 
$i=2$
, 
$\big\{c_k^{(i,j)}\big\}$
is explicitly given by the following:
\begin{align*}i=1\,{:}\,&\begin{cases}&c_0^{(1,2)}=\frac{1}{2}\\[7pt]
&c_0^{(1,1)}=\frac{1}{2};\ c_1^{(1,1)}=-1\\[7pt]
\end{cases}\\[10pt]
i=2\,{:}\,&\begin{cases}&c_0^{(2,4)}=\frac{1}{8}\\[7pt]
&c_0^{(2,3)}=\frac{1}{12};\ c_1^{(2,3)}=-\frac{1}{2}\\[7pt]
&c_0^{(2,2)}=-\frac{1}{8};\ c_1^{(2,2)}=0;\ c_2^{(2,2)}=\frac{1}{2}\\[7pt]
&c_0^{(2,1)}=-\frac{1}{8};\ c_1^{(2,1)}=\frac{1}{2};\ c_2^{(2,1)}=-\frac{1}{2}\end{cases}\end{align*}
Result A.2. (Decomposition structure II.) For 
$i\ge 3$
, 
$\left\{c_k^{(i,j)} \right\}$
is recursively determined by the following:
\begin{equation*}\begin{cases}j=2i\,{:}\,&\ c_0^{(i,2i)}=\frac{c_0^{(i-1,2i-2)}}{2i}\\
\\
j=2i-1\,{:}\,&\ \begin{cases}c_0^{(i,2i-1)}=\frac{1}{(2i-1)}c_0^{(i-1,2i-3)}-\frac{2i-3}{2(2i-1)}c_0^{(i-1,2i-2)}\\[5pt]
c_1^{(i,2i-1)}=\frac{1}{(2i-1)}\left(c_1^{(i-1,2i-3)}-c_0^{(i-1,2i-2)}\right)\\
\end{cases}\\\\
2<j<2i-1\,{:}\,&\ \begin{cases}k=(2i-j)\wedge i\,{:}\,\ \begin{cases}j>i\,{:}\, c_{2i-j}^{(i,j)}=\frac{1}{j}\left(c_{2i-j}^{(i-1,j-2)}-c_{2i-j-1}^{(i-1,j-1)}\right)\\[9pt]
j=i\,{:}\ c_{i}^{(i,i)}=-\frac{1}{i}c_{i-1}^{(i-1,i-1)}\\[5pt]
j<i\,{:}\ c_{i}^{(i,j)}=\frac{1}{2}(j+1)c_i^{(i,j+1)}\\[5pt]
\qquad -\frac{1}{j}c_{i-1}^{(i-1,j-1)}+c_{i-1}^{(i-1,j)}\end{cases}\\\\
0<k<(2i-j)\wedge i\,{:}\ c_{k}^{(i,j)}=\frac{1}{2}(j+1)c_{k}^{(i,j+1)}+\frac{1}{j}c_{k}^{(i-1,j-2)}\\[8pt]
\qquad\qquad\qquad -\frac{j-1}{j}c_{k}^{(i-1,j-1)}-\frac{1}{j}c_{k-1}^{(i-1,j-1)}+c_{k-1}^{(i-1,j)}\\\\
k=0\,{:}\ c_{0}^{(i,j)}=\frac{1}{2}(j+1)c_0^{(i,j+1)}+\frac{1}{j}c_0^{(i-1,j-2)}-\frac{j-1}{j}c_0^{(i-1,j-1)}\\
\end{cases}\\\\
j=2\,{:}&\ \begin{cases}k=i\,{:}&\ c_i^{(i,2)}=\frac{3}{2}c_i^{(i,3)}-\frac{1}{2}c_{i-1}^{(i-1,1)}+c_{i-1}^{(i-1,2)}\\[5pt]
0<k<i\,{:}&\ c_k^{(i,2)}=\frac{3}{2}c_k^{(i,3)}-\frac{1}{2}c_k^{(i-1,1)}-\frac{1}{2}c_{k-1}^{(i-1,1)}+c_{k-1}^{(i-1,2)}\\[5pt]
k=0\,{:}&\ c_0^{(i,2)}=\frac{3}{2}c_0^{(i,3)}-\frac{1}{2}c_0^{(i-1,1)}\end{cases}\\\\
j=1\,{:}&\ \begin{cases}k=i\,{:}&\ c_i^{(i,1)}=c_i^{(i,2)}+c_{i-1}^{(i-1,1)}\\[5pt]
0<k<i\,{:}&\ c_k^{(i,1)}=c_k^{(i,2)}+c_{k-1}^{(i-1,1)}\\[5pt]
k=0\,{:}&\ c_0^{(i,1)}=c_0^{(i,2)}\end{cases}\\
\end{cases}\end{equation*}
Appendix B. Proofs of Propositions 3.4 and 3.5
Proof. Part I (hyperbolic OU perturbed LT). We first consider a simplified version of (34). Let 
$h(x)=\tanh(cx)$
and define 
$\gamma \coloneqq \sqrt{2\beta}$
. According to the formula in Remark 2.1, the first-order solution is
\begin{equation}\tilde{f}_1= f_0\left[\int_{a}^x2e^{2\gamma y}\left(\gamma\int_{a}^y \tanh(cz)e^{-2\gamma z}dz+C_1\right)dy+C_2\right].\end{equation}
Rewrite
\begin{equation}\gamma\int_{a}^y \tanh(cz)e^{-2\gamma z}dz=\gamma \int_{a}^y \tanh(cz)e^{-cz}e^{cz-2\gamma z}dz,\end{equation}
and note that
\begin{equation*}\int \tanh(cz)e^{-cz}dz = -{\frac {\sinh \left( cz \right) }{c}}+2\,{\frac {\arctan \left( {{\rm e}^{cz}} \right) }{c}}+{\frac {\cosh \left( cz \right) }{c}}.\end{equation*}
Then integrating by parts in (B2) yields
\begin{align}&\gamma \int_{a}^y \tanh(cz)e^{-cz}e^{cz-2\gamma z}dz=-\gamma(c-2\gamma)\int_{a}^y 2\,{\frac {\arctan \left( {{\rm e}^{cz}} \right) }{c}}e^{cz-2\gamma z} dz\nonumber\\[2pt]
&+\gamma\left( -{\frac {\sinh \left( cz \right) }{c}}+2\,{\frac {\arctan \left( {{\rm e}^{cz}} \right) }{c}}+{\frac {\cosh \left( cz \right) }{c}}\right)e^{cz-2\gamma z}\bigg|_{z=a}^y\nonumber\\[2pt]
&-\gamma(c-2\gamma){\frac { \left( -\cosh \left( cz \right) +\sinh \left( cz \right) \right) {{\rm e}^{cz-2\,\gamma\,z}}}{2\gamma\,c}}\bigg|_{z=a}^y.
\end{align}
The first term on the right-hand side needs special care. Note that, for 
$m>1$
,
\begin{equation}\arctan(m)=\frac{\pi}{2}-\sum_{n=0}^\infty (\!-1)^n \frac{1}{(2n+1)m^{2n+1}}.\end{equation}
Substitute (B4) into (B3). After tedious calculations, we get
\begin{equation}\gamma\int_{a}^y \tanh(cz)e^{-2\gamma z}dz = I_1(y)-I_1(a),\end{equation}
where
\begin{align*}I_1(y)&=\gamma\left( -{\frac {\sinh \left( cy \right) }{c}}+2\,{\frac {\arctan \left( {{\rm e}^{cy}} \right) }{c}}+{\frac {\cosh \left( cy \right) }{c}}\right)e^{cy-2\gamma y}\nonumber\\[5pt]
&\quad -\gamma(c-2\gamma)\bigg[{\frac { \left( -\cosh \left( cy \right) +\sinh \left( cy \right) \right) {{\rm e}^{cy-2\,\gamma\,y}}}{2\gamma\,c}}\nonumber\\[5pt]
&\quad +\frac{2}{c^2} \left(\frac{\pi e^{cy(\!-\frac{2\gamma }{c}+1)}}{2(\!-\frac{2\gamma }{c}+1)}-\sum_{n=0}^\infty (\!-1)^n \frac{ {e^{-(2n+\frac{2\gamma }{c})cy} }}{(2n+1)(\!-2n-\frac{2\gamma }{c})}\right)\bigg].\end{align*}
Now substitute (B5) into (B1) and let 
$C_1=I_1(a)$
; then
\begin{align*}&\int_{a}^x2e^{2\gamma y}\left(\gamma\int_{a}^y \tanh(cz)e^{-2\gamma z}dz+C_1\right)dy=\int_{a}^x2e^{2\gamma y}\gamma I_1(y)dy\\[5pt]
&\quad =2\gamma \int_{a}^x\left( -{\frac {\sinh \left( cy \right) }{c}}+2\,{\frac {\arctan \left( {{\rm e}^{cy}} \right) }{c}}+{\frac {\cosh \left( cy \right) }{c}}\right)e^{cy}\\[5pt]
&\qquad -(c-2\gamma)\bigg[{\frac { \left( -\cosh \left( cy \right) +\sinh \left( cy \right) \right) {{\rm e}^{cy}}}{2\gamma\,c}}\\[5pt]
&\qquad +\frac{2}{c^2} \bigg(\frac{\pi e^{cy}}{2(\!-\frac{2\gamma }{c}+1)}-\sum_{n=0}^\infty (\!-1)^n \frac{ {e^{-2ncy} }}{(2n+1)(\!-2n-\frac{2\gamma }{c})}\bigg)\bigg]dy.\end{align*}
The integral of each of the terms on the right-hand side can easily be calculated. Let 
$C_2=0$
; summarising the results of the calculations, we get the first-order solution as
\begin{equation}\tilde{f}_1(x,\beta)=f_0(x,\beta)({g}_1(x,\beta)-{g}_1(a,\beta)),\end{equation}
where
\begin{align*}g_1(x,\beta)&=2\sqrt{2\beta} \left(2\,{\frac {{{\rm e}^{cx}}\arctan \left( {{\rm e}^{cx}} \right) }{{c}^{2}}}-{\frac {\ln \left( \left( {{\rm e}^{cx}} \right) ^{2}+1 \right) }{{c}^{2}}}-\frac{\pi e^{cx}}{c^2} \right)+x\\[5pt]
&\quad -\frac{4\sqrt{2\beta}(c-2\sqrt{2\beta})}{c^2} \left(\frac{cx}{2\sqrt{2\beta} }-\sum_{n=1}^\infty (\!-1)^{n+1} \frac{e^{-2ncx} }{c2n(2n+1)\Big(\!-2n-\frac{2\sqrt{2\beta} }{c}\Big)}\right).\end{align*}
Note that (B6) is not yet the LT to be inverted. It corresponds to the drift function 
$h(x)=\tanh(cx)$
, while what we want is 
$h(x)=\theta - \frac{\tanh(cx)}{c}$
. But we are almost there. By the linearity of the perturbed ODE solution (cf. Remark 2.1), the first-order solution of 
$h(x)=- \frac{\tanh(cx)}{c}$
is given by
\begin{equation}-f_0(x,\beta)\frac{{g}_1(x,\beta)-{g}_1(a,\beta)}{c}.\end{equation}
On the other hand, by computing the first-order solution of a drifted Brownian motion with 
$h(x)=\theta$
, we have
\begin{equation}\hat{f}_1(x,\beta)=-f_0(x,\beta)\theta(x-a).\end{equation}
Using the linearity property again, and combining (B7) and (B8), we obtain the final first-order solution of the hyperbolic OU process in (33) as
\begin{equation}{f}_1(x,\beta)=f_0(x,\beta)\Big(\!-\frac{{g}_1(x,\beta)-{g}_1(a,\beta)}{c}-\theta(x-a)\Big).\end{equation}
Part II (hyperbolic Bessel perturbed LT). Similarly as in Part I, we first consider the simplified drift function 
$h(x)=\coth(cx)$
and then generalise the result to 
$h(x)=c\coth(cx)$
. Using the equation in Remark 2.1 again, for the drift function 
$h(x)=\coth(cx)$
we have
\begin{equation}\tilde{f}_1=f_0\left[\int_{a}^x2e^{2\gamma y}\left(\gamma\int_{a}^y \coth(cz)e^{-2\gamma z}dz+C_1\right)dy+C_2\right].\end{equation}
Repeating the same trick as in (B2), and applying integration by parts, we then get
\begin{align}
\gamma \int_{a}^y \coth(cz)e^{-2\gamma z}dz & = \gamma\left(\frac {e^{-cz}+\ln(1-e^{-cz})-\ln(1+e^{-cz}) }{c}\right)e^{cz-2\gamma z}\bigg|_{z=a}^y\nonumber\\[8pt]
&\quad -\gamma(c-2\gamma)\int_{a}^y \left(\frac {e^{-cz}+\ln(1-e^{-cz})-\ln(1+e^{-cz}) }{c}\right)e^{cz-2\gamma z} dz.\end{align}
For the second term on the right-hand side, we consider the expansions
\begin{equation}\ln(1-e^{-cz})=-\sum_{n=1}^\infty \frac{e^{-ncz}}{n},\ \ln(1+e^{-cz})=\sum_{n=1}^\infty (\!-1)^{n+1}\frac{e^{-ncz}}{n}.\end{equation}
Substituting the expansions into (B11) and repeating the same calculations as in Part I, we find
\begin{equation}\gamma \int_{a}^y \coth(cz)e^{-2\gamma z}dz=I_2(y)-I_2(a)\end{equation}
with
\begin{align*}I_2(y) &=\gamma\left(\frac {e^{-cy}+\ln(1-e^{-cy})-\ln(1+e^{-cy}) }{c}\right)e^{cy-2\gamma y}\\&\quad -\frac{\gamma(c-2\gamma)}{c}\left(\!-\frac{e^{-2\gamma y}}{2\gamma}+2\sum_{n=0}^\infty \frac{e^{-(2nc+2\gamma)y}}{(2n+1)(2nc+2\gamma)}\right).\end{align*}
The remaining calculations follow the same routines as in Part I: substituting (B13) into (B10), letting 
$C_1=I_2(a)$
, and continuing the outer-integral calculations; setting 
$C_2=0$
and using again the linearity of the perturbed LT. In the end, for the drift function 
$h(x)=c\coth(cx)$
, the first-order solution is given by
\begin{equation}f_1(x,\beta) = f_0(x,\beta)(g_1(x,\beta) - g_1(a,\beta)),\end{equation}
where 
$g_1(x,\beta)$
is redefined as
\begin{align*}{g}_1(x,\beta)&=2\gamma\bigg(\!-x-{\frac {1-{e^{c x}}}{c}\ln \left( 1- { e^{-c x}} \right) }-{\frac {1+{ e^{c x}}}{c}}\ln \left( 1+ { e^{-c x}} \right)\\[5pt]
&\quad+(c-2\gamma)\bigg(\!-\frac{x}{2\gamma }+\sum_{n=1}^\infty \frac{ e^{-2ncx}}{nc(2n+1)(2nc+2\gamma)}\bigg)\bigg).\end{align*}
Part III (explicit inverse of perturbed LTs). Recall that in (6), the first-order perturbed LT is
\begin{equation}f^{(1)}=f_0+\epsilon f_1,\end{equation}
where, depending on the underlying process, 
$f_1$
is given by (B9) or (B14).
Referring to the 
$g_1$
-functions in (B9) and (B14), one can easily find that the LT parameter 
$\beta$
involved in the equation above belongs to the following structures (recall that 
$f_0=(x-a)/\sqrt{2\pi t^3}e^{-(x-a)^2/2/t}$
):
\begin{equation*}e^{-\sqrt{\beta}\sqrt{\alpha}},\ e^{-\sqrt{\beta}\sqrt{\alpha}}\sqrt{\beta},\ e^{-\sqrt{\beta}\sqrt{\alpha}}\frac{\sqrt{\beta}}{\sqrt{\beta}+\xi},\ e^{-\sqrt{\beta}\sqrt{\alpha}}\frac{\beta}{\sqrt{\beta}+\xi},\end{equation*}
where 
$\alpha=2(x-a)^2$
and 
$\xi=\frac{cn}{\sqrt{2}}$
are constant variables in the inverse transform. By [Reference Bateman6, Section 5.6, Eqs. (1), (5), (15), and (13)], we have
\begin{equation}{\mathcal{L}^{-1}}{\left\{{e^{-\sqrt{\beta}\sqrt{\alpha}}}\right\}}(t)=\frac{|x-a|}{\sqrt{2\pi }}t^{-\frac{3}{2}}e^{-\frac{(x-a)^2}{2t}},\end{equation}
\begin{equation}{\mathcal{L}^{-1}}{\left\{{e^{-\sqrt{\beta}\sqrt{\alpha}}\sqrt{\beta}}\right\}}(t)=\frac{t^{-\frac{5}{2}}}{2\sqrt{\pi}}\left((x-a)^2-t\right)e^{-\frac{(x-a)^2}{2t}},\end{equation}
\begin{align}
{\mathcal{L}^{-1}}{\left\{{e^{-\sqrt{\beta}\sqrt{\alpha}}\frac{\sqrt{\beta}}{\sqrt{\beta}+\xi}}\right\}}(t)&= \frac{1}{\sqrt{2\pi}}\left(\frac{|x-a|t^{-\frac{3}{2}}}{2}-cnt^{-\frac{1}{2}}\right)e^{-\frac{(x-a)^2}{2t}}\nonumber\\
& \quad +\frac{c^2n^2}{2}e^{cn|x-a|+\frac{c^2n^2t}{2}}\text{Erfc}\left(\frac{|x-a|}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\right),\end{align}
and
\begin{align}
{\mathcal{L}^{-1}}{\left\{{e^{-\sqrt{\beta}\sqrt{\alpha}}\frac{\beta}{\sqrt{\beta}+\xi}}\right\}}(t)&= \frac{t^{-\frac{3}{2}}}{2{\sqrt{\pi}}}\left(\frac{(x-a)^2}{t}-1-cn|x-a|+c^2n^2t\right)e^{-\frac{(x-a)^2}{2t}}\nonumber\\
&\quad-\frac{c^3n^3}{2\sqrt{2}}e^{cn|x-a|+\frac{c^2n^2t}{2}}\text{Erfc}\left(\frac{|x-a|}{\sqrt{2t}}+\frac{cn}{\sqrt{2}}\sqrt{t}\right).\end{align}
Finally, substituting (B16)–(B19) into the inverse of (B15) with the corresponding 
$f_1$
-functions from either (B9) or (B14) and summarising the results, one finds the perturbed density functions in Propositions 3.4 and 3.5. We leave the calculations to the reader and conclude the proof here.
Acknowledgements
The authors thank the referees and the associate editor of the Applied Probability journals for providing insightful comments and helpful suggestions on the first two versions of this paper.
