A simple method to calculate first-passage time densities with arbitrary initial conditions

Numerous applications all the way from biology and physics to economics depend on the density of first crossings over a boundary. Motivated by the lack of general purpose analytical tools for computing first-passage time densities (FPTDs) for complex problems, we propose a new simple method based on the independent interval approximation (IIA). We generalise previous formulations of the IIA to include arbitrary initial conditions as well as to deal with discrete time and non-smooth continuous time processes. We derive a closed form expression for the FPTD in z and Laplace-transform space to a boundary in one dimension. Two classes of problems are analysed in detail: discrete time symmetric random walks (Markovian) and continuous time Gaussian stationary processes (Markovian and non-Markovian). Our results are in good agreement with Langevin dynamics simulations.


Introduction
When the electric potential between the interior and exterior of a neurone exceed a certain threshold, the neurone fires. After firing, the interior potential is abruptly reset to its rest value and the process starts over. How often it starts over depends on external stimuli (e.g. light and touch) and firing frequencies of neighbouring neurones. To better understand neurone firing, and ultimately how neurones work, researchers in the field [1,2] use stochastic models to calculate how long it takes for the interior potential to pass the firing threshold for the first time.
Neurone dynamics is not the only case where first-passage problems arise. Such problems frequently occur in physics, chemistry, biology, ecology and economics [3,4] and is one of the reasons why first-passage problems are so heavily studied. But despite enormous interest there are surprisingly few cases where we know the probability density of first-passage times analytically.
Most analytical results are for Markov processes that mainly comes from two approaches. In the first approach, the so-called method of images, one solves the Fokker-Planck equation with absorbing boundaries [5,6]. Even though conceptually simple, it is limited to symmetric problems such as when the absorbing boundary is at the bottom of a symmetric potential well. The second approach is renewal theory [5,7]. It works for non-symmetric problems but often leads to expressions in Laplace-space that cannot be inverted analytically. Even though useful, both these approaches are in practice limited to simple problems. In fact, neither of them can provide the first-passage time density (FPTD) for a Brownian particle in a harmonic potential for a general boundary and starting point. Thus, in order to address more complex first-passage problems we need better analytical methods.
Another class of useful methods have been developed to solve persistence problems. In persistence problems one wishes to know the probability S(t) that a stochastic variable remains below or above a boundary from the start up to some time point t. The FPTD ( ) r t is simply related to the persistence according to To calculate the persistence, methods that involve the probability for all trajectories with an even number of boundary crossings has been used [8][9][10][11]. But apart from a few special cases, these crossing probabilities cannot be calculated exactly and one needs approximations. One approximation scheme that gained interest is the independent interval approximation (IIA) [12][13][14], which assumes that the length of time intervals between consecutive boundary crossings are independent. However, in its present formulation the IIA assumes that the processes has a well defined continuous velocity which means that it cannot deal with nonsmooth processes, such as discrete time processes or Brownian motion. Furthermore, in previous studies using the IIA, the initial condition was drawn from the equilibrium distribution. To apply IIA to a wider class of systems, these shortcoming must be remedied.
In this paper we generalise the IIA to discrete time series and non-smooth processes with arbitrary initial conditions. Starting with the discrete time case, we find a simple expression for the probability density of firstpassage times to a boundary in z-transform space. We then generalise our equations to the continuous time case and obtain a similar expression but now in Laplace transform space. The expressions are based on return probability densities to the boundary and the probability that the stochastic variable is above the boundary at some time. To show the applicability of our results we study the discrete time symmetric random walk and nonsmooth continuous time Gaussian stationary processes (GSP), both Markov and non-Markov systems.

Methods
In this section we outline the IIA framework and derive an expression for the FPTD for continuous and discrete processes x(t) in one dimension. We denote the FPTD by ( | ) is the starting point and x=B is the location of the absorbing boundary. In discrete time we let n be the number of time steps and = D t n t , where Dt is the time increment. We develop the IIA framework for discrete processes and then show how to generalise it to continuous processes.
The IIA equations herein relates three core quantities: the FPTD ( | ) , and the return probability density that x(n) returns to B after a B-crossing either from above, ( ) y + n , or from below, ( ) yn , after n steps. The quantities ( ) w > n and ( ) y  n are inputs to our framework which one needs to calculate on case by case basis. The probability ( ) w > n is in general simple to calculate. We find it by integrating the probability density function ( | ) P x n x , 0 of x(n): The averaged return probability densities ( ) y  n on the other hand are more complicated and needs to be discussed further (averaging procedure detailed below).

IIA for discrete time processes
To better understand ( ) y  n , consider a discrete time process that pass through B repeatedly (see figure 1, top).
The number of steps that x(n) remains below B is denoted by ¼ T T , , 1 3 , and above B by ¼ T T , , 2 4 . The time to the first-passage, T 1 , is special because it depends on x 0 . The density of T 1 is simply ( | ) r T x B 1 0 . After the first B-crossing at T 1 , ( ) x T 1 ends up at some distance  D 0 1 above B, rarely precisely on B (i.e. ) D = 0 1 . To calculate the distribution of T 2 , we must consider the trajectory from + D B 1 back across B. We denote the distribution of T 2 by ( ) y D + T , 2 1 where we assume that the length of T 2 is independent on T 1 . This is the core assumption of the IIA and is true for Markov processes. At T 2 , the process crossed B from above and is below B by D 2 . To find the number of steps until the next crossing, T 3 , we must consider the trajectory from -D In appendix B we show simulation results for the overshoot distribution for the discrete Brownian walk which is well described by ( ) ( ) l p D = D 2 erfc 2 . Working with the averaged return probability densities ( ) y  n instead of ( ) y D  n, implies that we ignore fluctuations in áDñ and approximate the original process x(n) by a clipped process. The dynamics of the clipped process is: when x(n) crosses B, draw Δ from ( ) l D  , make a jump to  D B , and continue (see figure 1, bottom). The clipped process is obviously different from the true x(n) but simpler to handle analytically. But the difference is small. We show in appendix B for the discrete time Brownian walk that ( ) w > n for the clipped process is practically indistinguishable from the true one. Below we formulate the IIA equations based on the clipped process. We may calculate ( ) in terms of the number of B-crossings. Our derivation below is the discrete time version of the continuous time case in [8,13]. Note, however, that in [8,13]  To calculate ( ) p N 1 , assume that the first up-cross occurred at < n N 1 and that there is no down-cross between n 1 and N. Since n 1 can be anywhere from 0 to N, this gives , assume that the first up-cross occurred at < n N 1 , the first down-cross occurred at n 2 between n 1 and N, that the second up-cross happened at n 3 between n 2 and N, and no down-cross between n 3 and N. This gives Summing over all odd number of B-crossings we get ( ) w > N (see equation (3)). To solve equation (3) for , > z 1) which turns the convolutions in equation (7) into products. After summing the resulting geometric series we obtain This equation relates the FPTD to the probability ( ) w > z of being above the boundary, and to the return probability densities ( ) y + z and ( ) yz . This constitutes one of our main results in this paper.

IIA for continuous time processes
In order to obtain the continuous time version of equation (8) we proceeds as follows. When B is reached from below, the trajectory makes a jump to , where ò is a small constant. As   0 we approach the continuous time case. The overshoot distributions ( ) l D  for this cases is a Dirac delta function, which leads to In appendix B we show explicitly how ( ) w > t for the clipped Ornstein-Uhlenbeck process (OUP) converge to the real one as   0. To derive the IIA equations in the continuous time case, we proceed in the same way as for the discrete time case using p k (t), but with sums in equation (7) changed to integrals ( ) (see equation (3)) we obtain a similar geometric series as before that leads to After thermally averaging over the initial position x 0 , this result is equal to the one obtained in [8,13] where also expressions for the return densities can be found. However, in contrast to the results here, the derivations in [8,13] assume smooth processes.

Results
In this section we apply our main results (equations (8) and (10)) to (A) discrete time symmetric random walks and (B) non-smooth continuous time GSPs, a Markovian case and a non-Markovian case. We also show that our equations lead to Kramers escape (appendix C), and that they are consistent with the method of images (appendix D) in continuous time. To test the validity of our results we compare them to Langevin dynamics simulations (see appendix A for simulation details).

Discrete time symmetric random walks
A prominent example of a discrete time process is the Markovian symmetric random walk that evolves via . The jump length η is an independent and identically distributed random variable drawn from a symmetric distribution ( ) f h of mean zero. This distribution could for example be Gaussian, Lévy, exponential or uniform.
The FPTD for the symmetric random walk is however not known for general B and x 0 except in terms of a double Laplace transform, the so-called Pollaczeck-Spitzer formula, that no one thus far have been able to invert analytically [9]. The one exact result that exists is the universal Sparre-Andersen theorem [15] which applies to all symmetric and continuous jump distributions ( ) f h . It says that the persistence to stay above (or below) the Here we use Q(n) and the IIA formalism to put forward a simple summation formula for the FPTD for general B and x 0 . However, we first need to find the return probability densities ( ) y  n .

Analytical predictions
In this problem the process x(n) behaves in the same way on both sides of the boundary, which means that the return probability densities on either side of B are equal, ( ) ( ) y y =  n n. We approximate ( ) y n with the discrete derivative of Q(n) (see equation (5) 1. Here ( ) Q n is the unit step function (discrete Heaviside step function) that takes care of the initial condition ( which after inversion leads to . Indeed, expanding equation (13) for large n we get the random walk result: (14) is a generalisation, albeit approximative, of the Sparre-Andersen theorem to general boundary and initial conditions.

Simulations and numerical results
In figure 2 we compare equation (14) to simulations when ( ) f h is Gaussian. Overall we find good correspondence, especially as -B x 0 increases. But as it decreases, we start to see deviations for small times, e.g for -= B x 1 0 . The reason is that the overshooting start to play a role and the derivative of the persistence is no longer a good approximation to ( ) y n . From simulations we find that the average overshooting length is 0.63 that is comparable to -= B x 1 0 .

Analytical predictions
A zero mean GSP is completely characterised by its correlator ( ) If f(t) decays exponentially for large t it follows that the probability that x(t) does not change sign during a time interval of length t is asymptotically equal to [8,9,14] where the rate r is the so-called persistence exponent that in general is non-trivial to calculate. Since the process x (t) is symmetric around B=0, we have ( ) ( )  y y + -t t (for Markov processes they are identical). Denoting them by ( ) y t and using that ( ) , see equation (5) Note that for non-Markovian systems, ( ) y  t will in principle remember all B-crossings and therefore depend on x 0 . However, since we work in the asymptotic limit (see equation (15)) we neglect this contribution.
With ( ) y  t at hand, we may calculate the FPTD for a general GSP from equation (10). First we use that ( ) in Laplace space. Second we put ( ) y  s in g(s) (see equation (10)) so that . After inversion we find the expression ( | ) that is valid for any GSP with correlator f(t) decaying exponentially (Markovian or non-Markovian, smooth or non-smooth). For Gaussian processes ( ) w > t is given by are the mean and variance, respectively. To determine the crossing rates  r , we proceed as follows. First we get + r from the normalisation condition (10), which together with equation (19) leads to . To close the system we need one of the crossing rates  r . However, they cannot be calculated from our IIA equations directly and must be acquired from either experiments, simulations or other analytical means. In this paper, we will assume that B and x 0 are not too close to each other which means that we may estimater as the reciprocal of the mean-first passage time This can be understood as follows. When B is far away from the potential minimum x=0, up-crossing events are rare. When they are rare, the distribution of times between up-crossing events is approximately the same as the FPTD. Because the equilibration time of the process is much shorter than τ, the FPDT is asymptotically equal 1 for long times. In the subsequent section we apply equations (18)- (21) to (i) the Markovian OUP and (ii) a non-Markovian case of two coupled Ornstein-Uhlenbeck systems. But before we move on to specific examples, we clarify some of our method's limitations.
Note that equation (10) is exact within the IIA but in general we do not know ( ) y  t . If we assume that ( ) y  t decays exponentially we arrive at equation (18). This assumption is only asymptotically true for a general GSP. This means that equation ( . However, this would introduce new parameters as well as conditions for when to truncate the sum. Using the asymptotic behaviour of ( ) y  t makes our approach free of cut-off parameters.  (18), 'o' are results from Langevin dynamics simulations (averaged over 10 6 realisations), and 'x' is the approximation in [23] (see equation (E.1)). The insets show the behaviour at short times. 3 The equilibrium time is of order one whereas τ is about ten times larger already when » B 2, where ( ) t~B exp 2

Markov case: the OUP
The OUP is the only Markovian GSP [7] (i.e up to a trivial scaling in time and position (see equation (A.2))). The FPTD for the OUP is difficult to calculate explicitly [18]. The one exception is the symmetric case when the absorbing boundary is at the bottom of the harmonic well which can be solved with the method of images [19]. However, this does not work for the general problem. Instead, several efforts focused on the renewal equation in Laplace space. But because the renewal equation cannot generally be inverted analytically [20,21] numerical inversion [22] and series expansion around the poles [23] have been used. The last example [23] currently holds the best analytical approximation. But even though in principle exact, none of their expressions are on closed form and must be evaluated numerically. To work with their expressions one must specify at least one cut-off parameter (sometimes two) which in practise must be done by trail and error. We compare our results to [23] using their so-called integral representation (see appendix E for details).
To use equation (18), we need the mean, variance, and mean-first passage time ( -r 1 ) to B. For the OUP they are given by Equation (23) follows from the backward Fokker-Planck equation [24]. Below we use these relations in equations (18)- (20) to compare with Langevin dynamics simulations.

Non-Markov case: two coupled Ornstein-Uhlenbeck systems
To see how our IIA equations perform for a non-Markovian process we study the simple coupled system , Gaussian white noise. We are particularly interested in the variable x(t) that we want to represent a non-smooth GSP. This means that we have to choose the constant parameters a i and b i ( = i 1, 2, 3) with care. A non-smooth GSP is characterised by its small time behaviour of the correlator f(t). Indeed, if with the constant ¹ a 0 then the process is non-smooth [14]. Moreover, since the OUP is the only Markovian GSP all others must be non-Markovian. This follows from Doob's theorem which says that a GSP is Markovian only if its correlator f(t) is a single exponential [9]. With the following choices 4 . This is not a single exponential and since ( ) 1 5 2 we conclude that x(t) is a non-smooth and a non-Markovian GSP. We also find the mean and variance to be As a consistency check we see that ( , which is what a zero mean GSP requires. Although, we point out that there are other sets of parameters that lead to the same correlator as above. r t x B 0 from our IIA framework we need the crossing rates r ± . However, it is non-trivial to calculate the mean-first passage time for non-Markovian systems analytically (see e.g. [25]). But because the main purpose of this example is to investigate how well the IIA equations work for non-Markovian systems we will for simplicity extractr from simulations as the reciprocal of the mean first-passage time. With this the other crossing rate + r follows from equation (20).

Simulations and numerical results
To validate our method we compare ( | ) r t x For the OUP (figure 3) our analytical FPTD decays exponentially in the long time limit. This agrees with Kramers escape rate and our IIA equations capture this regime well for all values of -B x 0 . For short times there is a discrepancy between simulations and our IIA equations when -B x 0 gets smaller. This is because a considerable amount of first-passage events occur for short times before ( ) y  t has attained its exponential form, as discussed above. Also, the reason why our IIA equations systematically underestimates the simulations for short times is the following. Assume for simplicity that = x 0 0 and that B is infinitesimally above x 0 . In this case, the first-passage dynamics is practically indistinguishable from Brownian motion for short times. For Markov processes we can show that ( ) y  t is a Dirac delta function (appendix D) which means that there is an infinite number of boundary crossings once the boundary is crossed. This is very different from ( )  y  - t e r t that says that the average time between two boundary crossings is  r 1 .
In figure 3 we also included one of the best analytical approximations for the FPTD [23] (see appendix E). Their formula approximates the short time dynamics better than our method while for long times both approaches match well with each other.
It is straightforward to generalise our method to two boundaries, B we need four return probability densities which our method is unable to handle. But for the symmetric case, | | = ¢ B B , ( ) y  t are enough to describe the crossing in and out of the region   -B x B. To calculater we use a generalisation of equation (21) to two boundaries (see e.g. [24]), and to get + r we use equation (20), where we replace ( ) . Figure 4 shows that our IIA equations match simulations for two boundaries with similar accuracy as for one boundary.
In figure 5 we compare our analytical non-Markovian FPTD to Langevin dynamics simulations where we keep = y 0 0 in all cases. This explains why we see better correspondence with simulations when = x 1 0 than for =x 1 0 (see equation (27)). For simplicity we extractedr from simulations. Overall, the IIA follow the simulated data in a similar way as for the simple OUP. That is, as the value of B increases we get better correspondence to simulations, see figure 6. The reason our method works so well for non-Markovian dynamics can be traced back to the value of B and the initial condition. For large values of B, which we mainly consider, the process spends on average long times below the boundary such that up-crossing events becomes Poisson distributed [8]. This deletes the memory of the considered non-Markovian model and therefore it becomes well described by our IIA formula.  In figures 3-5 we kept B fixed and changed initial conditions. We set B=3 where our approximation hold well. However, as B gets closer to zero, we expect to see deviations. To investigate this we keep initial conditions fixed and change B, see figure 6. In both the Markovian and the non-Markovian case we see that theory and simulations do not correspond well as B gets smaller. This is expected because (i) ( ) y  t are not decaying exponentially and (ii) -r is now not well estimated by the mean-first passage time. Also, the non-Markovian case is clearly more sensitive to small values of B than the OUP.

Summary and outlook
There are plenty of examples where one wants to know the probability density of first-passage times to a boundary. To find this density, one often use the method of images or renewal theory. These approaches are however in practise limited to simple cases. To find better methods, we improved the so-called IIA, developed for persistence problems, that is limited to smooth stochastic processes where trajectories have well defined velocities and thermally averaged initial conditions. This excludes for example Brownian motion. We generalised the IIA to discrete time and continuous time processes with non-smooth trajectories that start from fixed initial conditions. In our derivation we replace the original trajectory by a clipped process. After this replacement, we can deal with the IIA equations exactly in discrete time since we can analytically handle the overshooting length Δ from each B-crossing. Our procedure is further extended to continuous time, where a regularisation scheme is required. From our IIA formalism we derive a simple expression for the FPTD to a general boundary and initial condition in one dimension [see equations (8) and (10)]. This expression relies on that we know the functional form of the return probability densities. But once it is known our approach is parameter free. To show the validity of our expression we applied it to the discrete time Brownian walk and in continuous time to the OUP, the only Markovian GSP as well as a non-Markovian GSP. In discrete time we use the Sparre-Andersen theorem which yields the return probability densities for all symmetric random walks. For the GSPs we use that the return probability densities decays exponentially for long times [8,9,14], , where we estimater (up-crossing rate) as the reciprocal of the mean firstpassage time which is known for the OUP [24]. Then, from the normalisation condition of the FPTD we determine the last parameter, + r (down-crossing rate). All cases match well with Langevin dynamics simulations. We also show that in continuous time our IIA equations (i) reproduces Kramers expression for the escape of a Brownian particle out of a harmonic potential, and (ii) that it becomes equivalent to the method of images for Markovian process with symmetry around x=B, e.g. when the boundary is at the bottom of a potential well.
Our IIA equations are new and we anticipate that they will have a wide applicability to previously intractable first-passage and escape problems. We further hope that this study will spur further activities into the development of general purpose analytical frameworks for dealing with first-passage time problems.

Appendix B. Supplementary figures
In figure 1 we illustrated how we approximate the real process x(n) by the so-called clipped process. Here in figure 7 we show explicitly how the clipped process differs from the real one for the continuous time OUP and the discrete time Brownian walk. To the left in figure 7 we show the continuous time case. Clearly we approach the analytical curve for ( ) w > t (see equation (19)) as we let   0 (see equation (9)) in our simulations. To the

Appendix C. Kramers escape
For long times our method is consistent with Kramers escape theory. To see this we use the final value theorem [27] and find that equation (10) in the limit  s 0 (long times) becomes after Laplace inversion Using equation (20) to eliminate + r gives When B is large enough we get that ( ) w -¥ » > 1 1, Kramers expression [17,28]. As an example, for the harmonic potential we find that ( ) w ¥ »

Appendix D. Method of images
The method of images can successfully give the FPTD for symmetric problems. Notably, it only works for Markov processes. Here we prove that our IIA formalism is consistent with this method. To show this, we first use the renewal approach to find the return probabilities. Consider a Markovian particle that diffuses between two boundaries separated by a small distance ò (see figure 8). We denote an 'up-cross' by crossing B from below and a 'down-cross' by crossing B from above. When x=B is crossed the particle jumps immediately to   B , as explained in section 2. When symmetry around x=B holds the distribution of times between an up-cross and a down-cross are both equal to the first-passage to a point that is a distance  away Using this in the renewal equation we find If we now let   0, up-and down-crossings occur to the same boundary and therefore ( ) With ( ) y t at hand we may derive the method of images formula. First we set ( ) ( ) y y =  s s in equation (10) and then take   0, which after inversion leads to  where in the last step we use the symmetry around x=B which gives the method of images formula.

Appendix E. Alili's formula
To compare our result for the OUP to one of the best known approximations we have implemented one of the formulas from [23]. In our notation it reads   figure 3.