1 Introduction

The generalized Ornstein–Uhlenbeck (OU) model describes a stochastic process with at least one equilibrium point. It provides a framework for a wide range of physical, biological and social systems, wherein stabilization is viewed in terms of a potential minimum, characterized by a negative Lyapunov exponent, and high-frequency fluctuations are interpreted in terms of specific noise forcing. For example, an OU process is used to study neuronal activity [1] and the time-evolution of trait values towards their evolutionary optima [2]. In a clinical setting the health of the hepatic dynamic equilibrium is fit to an OU process using maximum likelihood estimation [3]. Stochastic volatility, crucial for deducing stock returns or option pricing, is treated in terms of an OU process [4], as are the noise spectra of climate observations [5,6,7].

A particular stochastic model is commonly studied in terms of the survival probability, which is associated with the probability of one or more events occurring, or for the system to reach a defined threshold for the first time starting from a given initial position. Due to the generality of the question it addresses, survival analysis has been widely used in science and engineering. Examples include Feshbach resonances and the quantum Zeno effect (for example [8,9,10]), engineering reliability analysis [11], financial risk management [12], and event history analysis in sociology [13]. Moreover, in the specific case of an OU process survival analyses from neuroscience [14] and epidemiology [15, 16] to quantitative finance [17,18,19] and extreme value statistics of correlated random variables [20] demonstrate the ubiquity of the approach. Our survival analysis is broadly relevant to all systems that can be described by an OU process. For example, it can be shown [21] that the survival probability of Brownian motion with an absorbing boundary that moves in time t as \(\propto \sqrt{t}\) can be recast as an OU process with a fixed absorbing boundary using Lamperti’s (or Doob’s) transformation.

Let \(\{X(t),\, t \ge 0\}\) be an OU process starting at x. The time it takes for the state of a system to encounter a threshold \({X(t)=\beta }\) for the first time is variously called the first hitting time or first passage time, whose probability distribution function \({\zeta (t,x)}\), is defined as

$$\begin{aligned} {\zeta (t,x)=\frac{\partial }{\partial t}\text {Prob}\{T \le t\}}, \end{aligned}$$
(1)

where

$$\begin{aligned} {T\equiv \text {inf}\{t\,:\,X(t)>\beta ;X(0)=x\}.} \end{aligned}$$
(2)

The time integral of the probability distribution of the first passage time is the survival probability (see for example [22,23,24,25] for reviews).

Despite the broad applicability of the survival probability, a simple accurate analytical expression has been lacking, which is the primary motivation of our work. The exact mathematical form is constructed using the Laplace transform and its inverse, which results in a series expansion of special functions [26]. However, its complexity confounds practical implementation. In particular, when the initial data are close to the boundary of the confining potential [27], or when the boundary itself is near the equilibrium point [28], one must retain a considerable number of terms in the expansion.

Fig. 1
figure 1

Schematic of the problem under study; what is the probability that a particle reaches the edge of the potential well? We divide the potential into an “outer” region I and a boundary layer II, or the “inner” region, of the potential, in each of which the asymptotically dominant solutions are determined and then matched

Full size image

Our treatment of survival probability is shown in the schematic potential of Fig. 1, which contains a particle governed by a one-dimensional Ornstein–Uhlenbeck process through an overdamped autonomous Langevin equation given by

$$\begin{aligned} {{\dot{X}}(t)=-aX(t)+\sqrt{2 b}\xi (t)}, \end{aligned}$$
(3)

where a and b are positive constants and \(\xi (t)\) is Gaussian white noise with a zero mean and \(\langle \xi (t)\xi (s)\rangle =\delta (t-s)\).

The survival probability commonly characterizes the anomalous or abnormal behavior of a system. This motivates our consideration of a threshold state (\({X=\beta }\)) far from the equilibrium state (\({X(t)=0}\)), so that it is rare that the system will reach the threshold.

The next section is devoted to our overall approach. In Sect. 2.1 we summarize our method. In Sect. 2.2 we briefly review previous results on the survival probability of an OU process that will be used in Sect. 2.3 to derive a simple analytical expression that reproduces the survival probability for large values of the domain boundaries and when the initial data are close to the boundary. In Sect. 2.4 we compare the analytical and numerical results.

2 Using Matched Asymptotic Methods to Determine the Survival Probability

2.1 Summary of Analytical Method

Our approach is as follows. As shown in Fig. 1 we divide the domain into two regions [29]: a broad O(1) region (I) containing the minimum of the potential, \({X=0}\), and a narrow \(O\left( {1}/{\beta }\right) \) boundary layer near \({X=\beta }\). We solve the limiting differential equations in these regions, from which we develop a uniform composite solution for the probability density of the survival probability using asymptotic matching, which is a highly accurate general analytical method [30], recently used to obtain the solution to the related problem of stochastic resonance [31].

2.2 Previous Results on the Survival Probability of an OU Process

The canonical approach of finding the survival probability of this OU process is to determine the probability distribution for the first hitting time in Laplace space

$$\begin{aligned} {{\tilde{\zeta }}(\lambda ,x)=\int _{0}^\infty e^{-\lambda t}\zeta (t,x) dt,} \end{aligned}$$
(4)

which leads to the following analytical expression [32,33,34],

$$\begin{aligned} {\tilde{\zeta }}(\lambda ,x)= \, e^{\frac{a(x^2-\beta ^2)}{4b}}\frac{D_{-\lambda / a}(-x \sqrt{a/b})}{D_{-\lambda / a}(-\beta \sqrt{a/b})}, \end{aligned}$$
(5)

where \(D_{\lambda }(z)\) is the parabolic cylinder function, \(\beta \) is the boundary position and x is the initial position.

This equation can be written as a spectral decomposition [35], wherein the countable number of eigenvalues correspond to the zeros of the denominator of Eq. (5). The final expression for the probability distribution of the first hitting time can be obtained by inverting each term of this spectral decomposition, which can be written as a weighted sum of exponential functions as

$$\begin{aligned} \begin{aligned} \zeta (t,x)=\sum _{p=0}^\infty e^{\lambda _p(\beta ) t}\frac{\varPhi (\lambda _p(\beta ) / a,x \sqrt{a/b})}{\varPhi '(\lambda _p(\beta ) / a,\beta \sqrt{a/b})}, \end{aligned}\end{aligned}$$
(6)

where \(\lambda _p(\beta )\) is the solution of

$$\begin{aligned} \varPhi (\lambda _p(\beta ) / a,\beta \sqrt{a/b})=0, \end{aligned}$$
(7)

and \(\varPhi \) and \(\varPhi '\) are defined in [35] in terms of Kummer and digamma functions.

For large values of t and \(\beta -x\), Eq. (6) simplifies considerably. In fact, because all the eigenvalues, \(\lambda _p(\beta )\), are negative, only that with the smallest magnitude will contribute significantly as \(t \rightarrow \infty \). Moreover, for large values of \(\beta -x\), Eq. (6) becomes

$$\begin{aligned} \begin{aligned} \lim _{\beta -x\rightarrow \infty }\zeta (t,x)&=\frac{1}{t_1(\beta ;0)} \exp \left[ -\frac{t}{t_1(\beta ;0)}\right] \qquad \text {with}\\ t_1(\beta ;0)&= \frac{1}{a\beta } \sqrt{\frac{2 \pi b}{a}}\exp \left[ \frac{a\beta ^2}{2b}\right] , \end{aligned}\end{aligned}$$
(8)

where \(t_1(\beta ;0)\) is the mean first passage time of the process from \(X(t=0)=0\) to the boundary \(\beta \) [36].

Importantly, however, the asymptotic expression in Eq. (8) is not valid when \(x\approx \beta \). Indeed, as we show here, when the process starts in the neighborhood of \(x=\beta \) there is a non-trivial leakage of probability. This leakage is not taken into account by Eq. (8) and transpires very rapidly, on a time scale of order \(1/\beta ^2\). Therefore, taking this approach requires that one calculate an enormous number of eigenvalues in Eq. (6), which is computationally inefficient. Our approach avoids this problem.

2.3 Detailed Development of the Approach

We derive an asymptotic expression for the survival probability, which is the time integral of the probability distribution of the first passage time, and in the large \(\beta \) limit it is trivial to evaluate when \(x\approx \beta \).

The probability density, \(\rho (y,t;x,s)\), of the OU process in Eq. (3) is described by the Kolmogorov forward (KFE) and backward (KBE) equations, the former of which is

$$\begin{aligned} \partial _t \rho (y,t;x,s)&= a~\partial _y[y(t)\rho (y,t;x,s)]\nonumber \\&\qquad +\,b~\partial _{yy}\rho (y,t;x,s)\,\equiv {\mathcal {L}}_y\rho (y,t;x,s), \end{aligned}$$
(9)

where the operator \({\mathcal {L}}_y\) is the generator of the OU process viz., \([{\mathcal {L}}_y f](y,t)=a~\partial _y [y f(y,t)]+b~\partial _{yy}f(y,t)\). The KBE follows by replacing \({\mathcal {L}}_y\) with its adjoint, \({\mathcal {L}}_x^*\), defined as \([{\mathcal {L}}_x^* f](x,s)=-ax~\partial _x f(x,s)+b~\partial _{xx}f(x,s)\).

The KFE gives the evolution of the probability density of the process when the initial position and time, (x, s), are known, while the KBE treats the evolution when the final position and time, (y, t), are known. Both equations have initial condition

$$\begin{aligned} \rho (y,u;x,u)=\delta (y-x), \end{aligned}$$
(10)

and boundary conditions

$$\begin{aligned} \rho (\pm \infty ,t;x,s)=0 \qquad \text {and} \qquad \rho (y,t;\pm \infty ,s)=0 \end{aligned}$$
(11)

for the KFE and the KBE respectively.

The survival probability within the interval \((\alpha ,\beta )\) is defined in terms of the KBE density, \(\rho _{KBE}\), as

$$\begin{aligned} S(t;x,s)=\int _{\alpha }^{\beta }\rho _{KBE}(y,t;x,s)dy, \end{aligned}$$
(12)

which satisfies

$$\begin{aligned} \begin{aligned} -\partial _s S(t;x,s)=-a~x(s)\partial _x S(t;x,s)+b~\partial _{xx} S(t;x,s), \end{aligned}\end{aligned}$$
(13)

with initial condition

$$\begin{aligned} S(t=s;x,s)=\varTheta (x-\alpha )\varTheta (\beta -x), \end{aligned}$$
(14)

and boundary conditions

$$\begin{aligned} S(t;x=\alpha ,s)=S(t;x=\beta ,s)=0 \;\;\; \forall \;s \le t, \end{aligned}$$
(15)

where \(\varTheta (\cdot )\) is the Heaviside theta function.

Let \(P_i^\alpha (P_i^\beta )\) be the probability of hitting \(x=\alpha (\beta )\) for the first time after i time steps and let \(S^\alpha (S^\beta )\) be the survival probability with an absorbing boundary at \(x=\alpha (\beta )\). Hence, if we discretize the stochastic process, the survival probability after n time steps is

$$\begin{aligned} \begin{aligned} S_n&=\prod _{i=1}^n(1-P_i^\alpha -P_i^\beta )=\prod _{i=1}^n[(1-P_i^\alpha )(1-P_i^\beta ) -P_i^\alpha P_i^\beta ]\\&\simeq \prod _{i=1}^n(1-P_i^\alpha )\prod _{i=1}^n(1-P_i^\beta )=S_n^\alpha S_n^\beta . \end{aligned}\end{aligned}$$
(16)

In the second line of Eq. (16), we have neglected the term \(P_i^\alpha P_i^\beta \) for \(|\alpha |\gg 1\) and \(|\beta | \gg 1\), allowing us to write the survival probability in the interval \((\alpha , \beta )\) as the product of the two survival probabilities in the two intervals \((-\infty , \beta )\) and \((\alpha , \infty )\) with \(\alpha <\beta \). From this point on we will only consider the survival probability in the interval \((-\infty , \beta )\), and note that the derivation for the interval \((\alpha , \infty )\) follows in straightforward analogy.

We now rewrite Eq. (13) in a rescaled form,

$$\begin{aligned} \begin{aligned} \partial _t S(x,t)=-x(t)\partial _x S(x,t)+\partial _{xx} S(x,t), \end{aligned}\end{aligned}$$
(17)

with the new variables,

$$\begin{aligned} x \rightarrow x\sigma , \;\;\; \beta \rightarrow \beta \sigma ~~\text {and}~~ s \rightarrow -\frac{t}{a}, \end{aligned}$$
(18)

in which the spatial coordinate is expressed in terms of the standard deviation \(\sigma =\sqrt{b/a}\) of a stationary OU process obtained from the solution of Eq. (9) in the limit \(t \rightarrow \infty \). The new initial and boundary conditions are

$$\begin{aligned} {\left\{ \begin{array}{ll} S(x,t=0)=\varTheta (\beta -x), \\ S(x=\beta ,t)=0,\hbox {and}\\ S(x=-\infty ,t)=1. \end{array}\right. } \end{aligned}$$
(19)

We solve the limiting differential equations within the two regions, from which we construct an approximate uniform solution by asymptotic matching. We denote \(S_{\text {out}}(x,t)\) and \(S_{\text {in}}(x,t)\) as the solutions in region I and II respectively.

Region I. The outer solution is obtained by imposing \(\beta -x \gg 1\). In this limit the probability distribution of the first hitting time in Eq. (8) is valid, and its time integral gives the survival probability as

$$\begin{aligned} S_{\text {out}}(t){\simeq }\,\text {e}^{-\left( \frac{\beta }{\sqrt{2 \pi }} \text {e}^{-\frac{\beta ^2}{2} } \right) t} , \end{aligned}$$
(20)

where the rescaled \(\beta \) and t of Eq. (18) have been used.

Region II. In the boundary layer, or the inner region, we have \(x \sim \beta \), where the approximation of Region I is no longer valid. We let \(\epsilon \equiv 1/\beta \ll 1\) and introduce the following stretched coordinates,

$$\begin{aligned} \eta =\frac{x-\frac{1}{\epsilon }}{\epsilon } ~~\text {and}~~ \theta =\frac{t}{\epsilon ^2}, \end{aligned}$$
(21)

that we use to rewrite Eq. (17) as

$$\begin{aligned} \frac{1}{\epsilon ^2}\partial _\theta S_{\text {in}}(\eta ,\theta )=&-\left[ \frac{1}{\epsilon ^2}+\eta (\theta )\right] \partial _\eta S_{\text {in}}(\eta ,\theta ) \nonumber \\&+\frac{1}{\epsilon ^2}\partial _{\eta \eta } S_{\text {in}}(\eta ,\theta ), \end{aligned}$$
(22)

which at leading-order becomes

$$\begin{aligned} \partial _\theta S_{\text {in}}(\eta ,\theta )=-\partial _\eta S_{\text {in}}(\eta ,\theta )+\partial _{\eta \eta } S_{\text {in}}(\eta ,\theta ). \end{aligned}$$
(23)

Clearly, Eq. (23) is a diffusion equation for \(S_{\text {in}}(\eta ,\theta )\) along the characteristics

$$\begin{aligned} \frac{d\eta }{d\theta }=1 ~~\text {and}~~ \frac{d\rho }{d\theta }=1, \end{aligned}$$
(24)

which we can then write as

$$\begin{aligned} \partial _\rho S_{\text {in}}(\mu ,\rho )=\partial _{\mu \mu } S_{\text {in}}(\mu ,\rho ), \end{aligned}$$
(25)

wherein \(\mu \equiv \eta -\theta \) and \(\rho \equiv \theta \), so that the boundary condition becomes \(S_{\text {in}}(\mu =-\rho ,\rho )=0\).

We solve Eq. (25) by first finding its Green’s function, \(G(\rho ,\mu ;\nu )\), which satisfies

$$\begin{aligned} \partial _\rho G(\mu ,\rho ;\nu )-\partial _{\mu \mu } G(\mu ,\rho ;\nu )=\delta (\mu -\nu )\delta (\rho ). \end{aligned}$$
(26)

This Green’s function is associated with the probability density \(\rho _{KBE}\), satisfying the KBE, and hence the survival probability through Eq. (12). Hence, this density satisfies the following conditions;

$$\begin{aligned} {\left\{ \begin{array}{ll} \rho _{KBE}(\mu ,\rho =0;\nu )=\delta (\mu -\nu ), \\ \rho _{KBE}(\mu =0,\rho ;\nu )=\rho _{KBE}(\mu =-\infty ,\rho ;\nu )=0. \end{array}\right. } \end{aligned}$$
(27)

The Green’s function is (e.g., [37])

$$\begin{aligned} \begin{aligned} G(\rho ,\mu ;\nu )=&\frac{1}{\sqrt{4\pi \,\rho }} \bigg (\exp \left[ -\frac{(\nu -\mu )^2}{4\,\rho }\right] \\&-\exp \left[ -\frac{(\nu +\mu )^2}{4\,\rho }-\nu \right] \bigg ), \end{aligned}\end{aligned}$$
(28)

which, as noted above, now allows us to write the solution of Eq. (25) as

$$\begin{aligned} \begin{aligned} S_{\text {in}}(\mu ,\rho )&=\int _{-\infty }^0d\nu \int _{-\infty }^\infty d\phi \,G(\mu ,\rho ;\phi )\delta (\phi -\nu )\\ {}&=\int _{-\infty }^0 d\nu \,G(\mu ,\rho ;\nu ). \end{aligned}\end{aligned}$$
(29)

Upon integration and reversion to the original variable t and to the stretched coordinate \(\eta \), we find

$$\begin{aligned} \begin{aligned} S_{\text {in}}(\eta , t){\simeq }&\frac{K}{2}\text {e}^{\eta }\bigg (-\text {erfc}\bigg [\frac{1}{2}\sqrt{\frac{1}{t}}(-t-\eta )\bigg ]\\ {}&+\text {e}^{- \eta }\text {erfc}\bigg [\frac{1}{2}\sqrt{\frac{1}{t}}(-t+\eta )\bigg ]\bigg )\,{\equiv \, K \,S_{\text {in}}'(\eta , t)}. \end{aligned}\end{aligned}$$
(30)

We determine the constant K by requiring the outer limit of the inner solution to equal the outer solution;

$$\begin{aligned} \begin{aligned}&K=\lim _{\eta \rightarrow -\infty }S_{\text {in}}(\eta ,t)=S_{\text {out}}(t)\,{\simeq }\,\text {e}^{-\left( \frac{\beta }{\sqrt{2 \pi }} \text {e}^{-\frac{\beta ^2}{2} } \right) t} . \end{aligned} \end{aligned}$$
(31)

Therefore, the uniformly valid approximate composite analytical solution for the survival probability is

$$\begin{aligned} \begin{aligned} S(x,t)&= S_{\text {in}}(x,t)+S_{\text {out}}(t)-K\\&{\simeq } \frac{1}{2}\text {e}^{-\left( \frac{\beta }{\sqrt{2 \pi }}\text {e}^{-\frac{\beta ^2}{2} } \right) t} \text {e}^{\beta (x-\beta )} \\&\quad \times \bigg ({-~}\text {erfc}\bigg [\frac{1}{2}\sqrt{\frac{1}{t}}(-t-\beta (x-\beta ))\bigg )\\&\quad +\,\text {e}^{-\beta (x-\beta )}\text {erfc}\bigg [\frac{1}{2}\sqrt{\frac{1}{t}}(-t+\beta (x-\beta ))\bigg ]\bigg )\\&\equiv S{'}_{\text {in}}(x,t)S_{\text {out}}(t), \end{aligned} \end{aligned}$$
(32)

written in terms of the original spatial coordinate x rather that the stretched coordinate \(\eta \).

Fig. 2
figure 2

a–c Plot of Eq. (32), S(x, t), versus \(x \in [0,\beta ]\) for different values of t (solid lines) compared to the numerical solution of Eq. (17) (crosses). We use three different values of the boundary position \(\beta \); a 3.5, b 1.96, and c 1.645, corresponding to the probability of finding the system below \(x=\beta \) as \(t \rightarrow \infty \) and without absorbing boundaries of a \(99.95\%\), b \(99\%\), and c \(90\%\) respectively. d The numerical solution of Eq. (17) with two boundaries (red circles), \(S_{\alpha ,\beta }(x,t)\), \(S_{\alpha }(x,t) S_{\beta }(x,t)\) (blue crosses), and the approximate analytical solution in Eq. (33)\(, S'_{\alpha }(x,t) S'_{\beta }(x,t)\). Note that \(S_{\alpha ,\beta }(x,t)\) overlaps exactly with \(S_{\alpha }(x,t)S_{\beta }(x,t)\). Here \(\alpha =\beta =3.5\) and \(x=0\)

Full size image

2.4 Comparing Analytical and Numerical Results

We compare Eq. (32) with the numerical solution of Eq. (17) for different values of the parameter \(\beta \) in Fig. 2, and find excellent agreement even when \(\beta \) is not asymptotically large. Indeed, the smallest value of the distance to the boundary is \(\beta =1.645 \sigma \), corresponding to a \(90\%\) probability that an unbounded stationary OU process remains below the boundary.

When \(t\gg 1/\beta ^2\) the left hand side of Eq. (23) is negligible and Eq. (32) becomes \(S_{\text {out}}(t)[1-e^{\beta (x-\beta )}]\) and depends on time solely through \(S_{\text {out}}(t)\), which is the prefactor of \(S_{\text {in}}(x,t)\). Thus, depending on x the rate of the decrease in the survival probability is controlled by \(S_{\text {in}}(x,t)\), decaying more rapidly near the boundary for early times.

The accuracy of the asymptotic solutions over a wide range of the \(\beta \) facilitates simple and wide ranging applications. For example, determining the input parameters of a leaky integrate-and-fire (LIF) neural model (see e.g., [38] and refs therein) based on experimentally observable interspike intervals [28]. Of particular contemporary relevance is deducing the “critical community size” in disease epidemiology [27], or the population threshold below which infections do not persist.

Finally, we appeal to Eq. (16) to form the asymptotic expression for the survival probability of an Ornstein–Uhlenbeck process with two absorbing boundaries as

$$\begin{aligned} S_{\alpha ,\beta }(x,t) = S{'}^{\alpha }_{\text {in}}(x,t) S^{\alpha }_{\text {out}}(t)S{'}^{\beta }_{\text {in}}(x,t) S^{\beta }_{\text {out}}(t)\; \forall t>s. \end{aligned}$$
(33)

We note that only one of the boundary regions will contribute significantly viz.,

$$\begin{aligned} \begin{aligned}&S{'}_{\text {in}}^{\alpha ,\beta }(x,t) =S{'}^{\alpha }_{\text {in}}(x,t)S{'}^{\beta }_{\text {in}}(x,t) = S{'}^{\gamma }_{\text {in}}(x,t), \\&\gamma = \text {min}_{|h_\alpha |, |h_\beta |}(\alpha ,\beta ), \end{aligned}\end{aligned}$$
(34)

where \(h_{\alpha (\beta )}=x-\alpha (\beta )\).

In Fig. 2d we show that the numerical solution with two boundaries matches that obtained by considering the two boundaries separately and with the analytical solution in Eq. (33).

3 Conclusion

We have obtained an asymptotic analytical solution for the survival probability of an Ornstein–Uhlenbeck process for large potentials. We divide the potential into two layers near the boundaries and a broad region between, which contains the origin, where we use the solution of Ricciardi and Sato [35]. The uniformly continuous solution is obtained by matching the two approximate solutions in the boundary layers with that in the broad region between. The solution agrees extremely well with both the numerical solution and with the more restricted asymptotic expression for the survival probability known in literature. Importantly, our analysis remains valid even when the initial position of the stochastic process is close to one of the boundaries, and furthermore it takes into account the non-negligible leakage of the probability early in the time evolution.

Despite the analysis using the assumption of asymptotically large boundaries, we showed that it agrees well with the numerical solution, even when the boundary positions are the same order of magnitude as the standard deviation of the stationary probability distribution function of the OU process. Indeed, we demonstrate consistency when there is a \(90 \%\) probability that the system without absorbing boundaries is found at a position less than \(\beta \); when reaching the boundary can no longer be considered a rare event. Therefore, our method can be easily extended to more general (less restrictive) settings. Finally, our compact analytical solution provides a computationally trivial framework for survival analysis of use across the broad spectrum of stochastic systems where the Ornstein–Uhlenbeck process arises.