The moving-eigenvalue method: hitting time for Itô processes and moving boundaries

We present simple solutions of first-passage and first-exit time problems for general moving boundaries and general Itô processes in one dimension, including diffusion processes with convection. The approach uses eigenfunction expansion, despite the boundary time-variability that, until now, has been an obstacle for spectral methods. The eigenfunction expansion enables the analytical reduction of the problem to a set of equivalent ordinary differential equations, which can be input directly to readily available solvers. The method is thus suitable as a basis for efficient numerical computation. We illustrate the technique by application to Wiener and Ornstein–Uhlenbeck processes for a variety of moving boundaries, including cases for which exact results are known.


Introduction
Consider a stochastic process X t starting at t = 0 in the moving (time-variable) range [a(t), b(t)]. What is the probability S(t) that the process will remain in the interval at least until time t (figure 1)? When the boundaries a and b are finite, we call this the first-exit time problem. When b is infinite, we have the closely-related first-passage time problem. We refer to both of these problems as hitting-time problems.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Finding S(t), known as the survivor function, survival probability, or persistence, is a common problem in many fields, such as earth quake prediction in physics [1], reaction times in chemistry [2], neuron spiking in biology [3], structure monitoring in engineering [4], and risk assessment in economics [5], but problems of this type are notoriously difficult [6][7][8][9][10].
Here, we consider any process X t in one spatial dimension defined by the Itô stochastic differential equation (SDE) where W t is a standard Wiener process (Brownian motion), i.e., a continuous process having W 0 = 0 and independent Gaussian increments W u+t − W t ∈ N(0, u) for all t > 0, u 0; the drift μ(x, t) is continuously differentiable with respect to x, and the diffusion σ(x, t) > 0 is twice continuously differentiable with respect to x. In the common case where μ and σ are independent of t, we call the process time-homogeneous. We allow the boundaries a and b to be any regulated functions, meaning functions that can be written as uniform limits of step functions [11]. This wide class includes all differentiable functions with finite numbers of finite discontinuities.

Solutions
We use the term solution to mean the combination of two stages: first, a mathematically rigorous reduction of the original problem to an equivalent but numerically more tractable form, and then the numerical processing of this reduced form. This second stage can be purely numeric or partly symbolic, but it will always involve varying degrees of approximation. The quality of the first stage result is determined by the applicability of an efficient second stage, and how well it lends itself to finding abstract features of the solution, such as derivatives and asymptotic properties. The quality of the second stage is determined by its computational efficiency and accuracy.
Our contribution focuses on the first stage, where we show how hitting-time problems for moving boundaries can be reduced to sets of ordinary differential equations (ODEs) without resorting to approximations. These ODEs can then be solved easily and efficiently by a diversity of available solver programs. Although the solution of the ODEs is not the primary concern of this paper, we illustrate it by several computed examples, and provide some practical advice for applications.

Moving boundaries
Moving boundaries in models are obviously important because boundaries are not always constant in applications, but there is another, more profound, reason. A variable transformation can sometimes convert a complicated stochastic process into a simpler process, such as a standard Ornstein-Uhlenbeck process, but a common side-effect is for a constant boundary in the original problem to be changed into a moving boundary in the transformed problem.
Example 1 (Doob's transformation). Doob's transformation [12, section 5.9] converts the standard Wiener process W u with boundary b(u) by u = e 2t and X t e −t W u to a standard Ornstein-Uhlenbeck process satisfying the Langevin equation with the boundary Here, a boundary b that was originally constant is transformed into an exponential boundary, and a square root boundary b(u) = √ u is transformed into a constant boundary.
Example 2 (Lamperti's transformation). The noise term in (1) is multiplicative, or statedependent; a modelling example would be a machine which becomes noisier as it runs faster. Under mild conditions [13], Lamperti's transformation can convert multiplicative noise to additive by the substitution The new equation becomes Example 3 (SDE with space-independent coefficients). The process defined by the SDE can be transformed to a Wiener process W τ with variance 2 by the substitution with the side effect that the constant boundary b(t) = b 0 changes to the time-variable boundary [14].
. This changes the constant boundary b(t) = b 0 to the timevariable boundary U(t, b 0 ).

Scope and related work
Because of the huge number of publications in the field of hitting-time problems (for reviews, see, e.g. [3,6,[16][17][18][19][20][21]), we propose a simple taxonomy to introduce the relevant work and to indicate the scope of this paper.
Furthermore, we introduce minor divisions according to: • Process type: in discrete time, the process X t is most often a random walk, but birth-death processes are also common, especially in queuing theory, where some of the earliest hitting-time problems appeared (e.g., [32][33][34][35][36][37][38]). In continuous time, the process is usually a Wiener or an Ornstein-Uhlenbeck process, but many other types have been studied, e.g., diffusion processes [22,25,[39][40][41][42], Itô processes [43], piece-wise constant diffusion [29], time-dependent diffusion [44,45], constant diffusion with drift that solves Burger's equation [46], continuous Gauss-Markov processes [47], stationary Markov processes [48], Lévy processes [12,21,43,49], integrated Wiener processes [8], integrated Lévy processes [21], stationary Gaussian processes [50], Feller processes [16], Wiener processes with sine modulation [51], Markov additive diffusion [52], and relaxation processes, which do not terminate upon reaching a given threshold for the first time [53,54]. Discontinuous processes such as Lévy flights may jump over the border-known as a first arrival-before the first passage, but time-homogeneous Itô processes are always continuous. In this paper, we consider the class of Itô processes, which includes diffusion processes, such as Wiener and Ornstein-Uhlenbeck processes. Although we consider only continuous-time processes here, these can be seen as limiting cases of discrete-time processes, and Wiener and Ornstein-Uhlenbeck processes can be seen as limiting cases of random walks and birth-death processes, respectively. The relationship between discrete and continuous problems is close and has been addressed in detail by, for example, Borovkov [55] and Whitt [56]. • Boundary type: studies of moving-boundary problems are usually limited to specific types, such as linear [7], piece-wise linear [57], cubic [58], piece-wise cubic [8], square root [59], and power boundaries [6,17], but many authors consider larger classes of boundaries, typically smooth functions [16, 38, 40-43, 46, 47, 60-62]. Here, we allow the comprehensive class of regulated boundaries [11], which includes piece-wise continuously differentiable boundaries with finite discontinuities as a special case. • Initial conditions: in the literature, a deterministic initial value x 0 = X 0 of the process is almost always given as the initial condition, but there are a few authors who allow a given initial distribution π(x) for X 0 , e.g. [14,17,20,38,47,61,63]. Of course, a solution S(t, x 0 ) for arbitrary x 0 provides a Green's function that can be integrated to provide a solution for π(x) using S(t) = b(0) a(0) π(x)S(t, x) dx, but this has some subtle disadvantages: when we want to compute a numerical result, we may have faster convergence if we use the initial distribution directly in the computation rather than an additional integration. This becomes particularly pronounced for methods that perform piece-wise (spline) approximations of boundaries. Here, we assume that π(x) = p(x, 0), the distribution of X 0 , is given as the initial condition.
Only a few methods are capable of finding the survivor function for generic moving boundaries and Itô processes. These methods belong mainly to three groups: Monte Carlo, integral equations, and methods based on the Fokker-Planck partial differential equation (PDE).
The simplest method is Monte Carlo simulation, where the Itô equation (1) is used directly to generate sample paths [15]. Unfortunately, this method is slow, requires very short time steps for the accurate detection of boundary crossings, and is impractical for large t.
Methods in the second group convert the problem to an integral equation [41,42,47], but the kernel of this equation contains a weak singularity, causing complications. It is possible to formulate an integral equation including a function chosen to eliminate the singularity [42], but finding such a function is difficult.
The third group first converts the Itô equation into a Fokker-Planck equation [68]. This can be solved exactly only in rare cases [6], and it normally requires PDE-solving techniques, of which single-domain spectral methods [69] currently represent the state of the art. However, these methods are difficult to apply in the non-rectangular domains generated by moving boundaries, and one typically has to resort to domain decomposition approaches, such as finite-element methods [70].

Outline
Our method can certainly be classified as a spectral method because it uses eigenfunction expansion, but whereas this is often considered to be synonymous with the separation of variables, we show in this paper that eigenfunction expansion is more widely applicable and admits even inseparable PDEs.
In short, we look at the Fokker-Planck equation for fixed t as a Sturm-Liouville system. We then observe that moving boundary conditions can in fact be expressed neatly in terms of the time-variable eigenvalues and eigenfunctions of the Sturm-Liouville operator, enabling a transformation of the Fokker-Planck equation into a set of ODEs notwithstanding the moving boundaries. Importantly, the boundary condition representation in the transformed equations is concise and enables an elegant convergence proof. For example, in the time-homogeneous first-passage case, these equations can be expressed as where ω is a vector representing the probability distribution of the process at time t, a is the boundary (which is a function of t), Λ is the diagonal matrix of Sturm-Liouville eigenvalues λ k (which are functions of a), and C is the skew-symmetric boundary sensitivity matrix defined by After recollecting a few important facts from Sturm-Liouville theory in section 2, we solve the first-exit time problem for time-homogeneous SDEs (Itô diffusions) in section 3, and then the first-passage problem in section 4. In section 5, we generalize this to time-variable SDEs (Itô processes). We give examples by applying the theory to Wiener and Ornstein-Uhlenbeck processes in sections 6 and 7, respectively. Section 8 discusses the results, including some limitations and practical aspects, and concludes the paper.

Preliminaries: hitting time for constant boundaries
In this section we introduce our notation while recapitulating some useful results from Sturm-Liouville theory.

The Fokker-Planck equation
For brevity, we write the partial derivative ∂y/∂x of a function y(x, t) with respect to the space variable x as y .

Sturm-Liouville boundary value problems
A differential operator L of the form is known as a Sturm-Liouville operator, and the differential equation as a Sturm-Liouville differential equation. This equation, together with the Dirichlet boundary conditions is an example of a Sturm-Liouville boundary value problem, for which there is a wealth of theoretical results, such as the following well-known lemma: Lemma 1 (Green's formula). Given two solutions ϕ μ , ϕ ν of the Sturm-Liouville equation (5), corresponding to eigenvalues λ = μ and λ = ν, respectively, we have for μ = ν where C denotes an arbitrary constant. For μ = ν, Proof. Multiplying the Sturm-Liouville equation for ϕ μ with ϕ ν and vice versa gives Subtracting the equations and integrating, Rewriting and letting μ → ν in (7), we obtain (8).

Example 5 (Fokker-Planck equation in Sturm-Liouville form). After multiplying the Fokker-Planck equation (3) by the integrating factor
we can write the RHS of (3) in Sturm-Liouville form (5), The operator L for the boundary conditions (6) is a self-adjoint operator of the Hilbert space We define the norm of a function f ∈ L 2 w (I) to be f w f | f , where we sometimes drop the weight index w if there is no risk of confusion. We denote an eigenfunction corresponding to the eigenvalue λ k as ϕ λ k , or just as ϕ k ; a Greek index denotes an eigenvalue, whereas a Roman letter indicates the order of the eigenvalue. An eigenfunction ϕ k is normalized to Example 6 (Norm of eigenfunction over semi-infinite interval). Suppose that an eigenfunction ϕ ν satisfies the Dirichlet condition ϕ ν (a) = 0 and ϕ ν (x), ϕ ν (x) → 0 when x → ∞. By (8), the squared norm of ϕ ν is Differentiating the equation and substituting into (11), Let y 1 and y 2 be linearly independent solutions of a Sturm-Liouville equation (5), and define the eigenfunctions by The eigenvalues λ can then be solved from the equation ϕ λ (b) = 0. The derivative ϕ λ (a) = y 1 (a)y 2 (a) − y 2 (a)y 1 (a) = W{y 1 , y 2 }(a) (12) equals the Wronskian of y 1 and y 2 at x = a. It is easily verified that (PW) vanishes for the Wronskian W of a Sturm-Liouville equation, implying that P(a)ϕ λ (a) is independent of a.

Liouville's transformation
When P, w > 0 and P, P , w, w are locally absolutely continuous on I = [a, b], it is possible to simplify (4) and (5) by Liouville's transformation [71]. Define The transformed equation becomes where In the expression for Q, the derivatives are taken with respect to the original variable x.

Example 7 (Liouville transformation of the Fokker-Planck equation). Define
Then, for the Fokker-Planck equation (3), and Here, Example 8 (Liouville transformation for Ornstein-Uhlenbeck process). For the standard Ornstein-Uhlenbeck process defined by the Langevin equation (2),

Eigenfunction expansion for finite intervals
A classical result in Sturm-Liouville theory is that for finite intervals I = (a, b), the eigenvalues are real, simple, and countable; and the set of eigenfunctions of the Sturm-Liouville system (5) and (6) forms a complete, orthogonal basis for the weighted Hilbert space L 2 w (I), provided that 1/P, q, w ∈ L 1 (I) [72, section 4.6]. It follows from the Sturm separation theorem [72, section 2.6] that the zeros of an eigenfunction are simple, because they interleave the zeros of a linearly independent solution to the same differential equation.
The eigenvalues can be arranged in a sequence −∞ < λ 1 < λ 2 < . . ., and the eigenfunction ϕ k associated with the eigenvalue λ k has exactly k − 1 zeros in the open interval I. Any function f ∈ L 2 w (I) can be series expanded in terms of the eigenfunctions ϕ k , meaning that given we have or convergence in the norm of L 2 w (I). The solutions of (5) are continuously differentiable with respect to x and λ when P > 0 and P , q, w are continuous [73].
If P, q, and w are continuous at a and b, then λ k is differentiable there and [72, section 4.4] When P, w > 0, and k → ∞ [72, section 4.3], With certain smoothness constraints on f, the series expansion (15) converges absolutely and locally uniformly; see the next section for details.

Eigenfunction expansion for semi-infinite intervals
When the interval is semi-infinite, i.e., b = ∞, the spectrum is not necessarily discrete. However, Weyl [74] showed that if the Liouville transformation is applicable, Q is bounded from below, and Q(x) → ∞ when x → ∞, then the spectrum is still discrete and we can series expand any function f ∈ L 2 w (I) in the same way as in (15). A slightly more complex condition, the Molchanov criterion [75], is both necessary and sufficient: given that Q is bounded from below, the spectrum is discrete and bounded from below if and only if, for some h, Example 9 (Spectrum for Wiener processes). For the standard Wiener process, Q = 0, showing that the spectrum is not purely discrete for any semi-infinite interval.
Example 10 (Spectrum for Ornstein-Uhlenbeck processes). For the Ornstein-Uhlenbeck process defined by (2), we already found in example 8 above that Q(x) = x 2 /4 − 1/2, which is bounded from below. The spectrum is discrete and bounded from below, because Q(x) → ∞ when x → ∞.
Weyl [74, chapter II, theorem 4; chapter III, theorem 7] showed that the series expansion converges absolutely and locally uniformly if f and Lf are continuous and in L 2 w (I). For the case P = w = 1, Titchmarsh [76, chapter II, section 2.7] presented the same convergence result, but required instead that q be continuous, f be the integral of an absolutely continuous function, L f ∈ L 2 w (I), and for all non-real λ. By appropriately padding f, P, q, and w from x = b to x = ∞, the result can be applied to eigenfunction expansion in the finite interval case. Some alternative conditions for discreteness are given by Müller-Pfeiffer [77].

Solving parabolic PDEs for rectangular domains
The well-known separation of variables method for solving PDEs such as (9) for a constant boundary or rectangular domain in the time-homogeneous case expresses the solution in terms of an orthonormal basis {ψ k } of eigenfunctions of the corresponding Sturm-Liouville system, here and then assumes that there is a solution in the form We shall use a similar approach, and for regularity conditions and questions of existence and uniqueness, we refer to Evans [78] and Lieberman [79]. Using (19), the PDE (9) is reduced to a set of equations, one for each eigenvalue, Projecting the equation on ψ k , i.e., taking the inner product with ψ k , we find The initial values ω k (0) are the Fourier coefficients of p at t = 0, For short time intervals Δt, we have

The survivor function
Our final target is the survivor function S, which expresses the probability that the process has not yet reached the boundary, where Ψ k are the integrals of the normalized eigenfunctions, provided that the integration and summation order can be switched, that is, when the series is absolutely convergent and its integral is finite.

First-exit time for time-homogeneous processes
We , and this is the reason why we call the proposed method 'the movingeigenvalue method'.

One constant and one step function boundary
We first assume that the upper boundary b is constant, while the lower boundary a(t) < b is a single-step function. We denote the kth eigenvalue by λ k (a, b), corresponding to the eigenfunction ϕ k . In this section and the next, we abbreviate the partial derivative ∂λ k (a, b)/∂a as λ k . Starting with (19), suppose the boundary a takes a step Δa > 0 up at t ( figure 3). This step causes a change in the eigenfunction basis and the inner product. Before the step, the kth eigenfunction is ϕ k = ϕ λ k (a,b) , whereas after the step, the eigenfunction changes to ϕ + k = ϕ λ k (a+Δa,b) . The inner product (10) changes to The step effectively 'truncates' the bottom part (x < a + Δa) of the density p(x, t) to and We can simplify this expression using lemma 1 and the fact that By choosing k = m and letting Δa → 0 we have Similarly,  (24),

Substituting (25)-(27) into
If, instead, Δa < 0, the boundary takes a step down (figure 4). Denote by ϕ k the zeropadded extension of ϕ k to [a + Δa, b]. We now have Equations (28), (29), and (21), which holds for any differentiable λ(t), allow us to solve the first-exit time problem for step function boundaries; where the boundary is constant, we use (21), whereas when the boundary has a step, we use (28) or (29).

One constant and one smooth boundary
Suppose that the boundary can be sandwiched between an outer step function and an inner step function ( figure 5). The first-passage probability for an inner boundary will at any time be larger than or equal to the probability for an outer boundary. Therefore, for a sandwiched boundary, the first-passage probability will be bounded above and below by the hitting times of the two step functions, so any function that can be sandwiched between step functions converging uniformly to a common limit allows us to compute the first-exit time by taking the limit for the sequence of step function boundaries, if such a limit exists. The existence is guaranteed for the class of regulated functions [11], which includes differentiable functions as special cases.
For small positive Δa in (28), whereas for small negative Δa in (29), We have for both negative and positive Δa, We note that c km = −c mk when k = m. The first factor disambiguates the choice of polarity of the eigenfunctions. Assuming that ϕ k is chosen so that ϕ k (a) > 0, the Sturm separation theorem implies that its zeros are simple (section 2.4) and that the first factor in the expression for c km reduces to one in the limit as Δa approaches zero. It follows for the special case k = m that c kk = 1 + O(Δa 2 ).
Combining (22) and (23), or We can summarize these equations in a vector formula. Let a(t) be a differentiable boundary and let Δt approach zero. In the limit, where ω = (ω 1 , ω 2 , . . .) T and Λ = diag({λ m }). We call the matrix C a the boundary sensitivity matrix corresponding to the boundary a, and define it as the skew-symmetric matrix whose off-diagonal components are

Convergence
Do the solutions ω n of the system truncated to n dimensions converge to a limit ω when n approaches infinity? We know that |ω(0)| must be finite if p(x, 0) ∈ L 2 w . By Peano's theorem, the continuity of C a guarantees that ω n exists in a neighbourhood of zero. For |ω n |, we have By Grönwall's lemma, |ω n (t)| |ω n (0)|e −λ 1 t for positive t. By making n large enough, |ω n (0) − ω n+m (0)| can be made arbitrarily small. Because λ 1 is bounded from below, Cauchy's convergence criterion implies that ω n and its derivative dω n /dt converge locally uniformly to limit functions ω and dω/dt, respectively.

Two smooth boundaries
Holding boundary a constant and instead taking a step Δb < 0 (down) at the upper boundary b, we have ϕ k (a) = ϕ k (b) = ϕ + m (a) = ϕ + m (b + Δb) = 0. By an argument analogous to that in the previous two sections, we have In the same way as in section 3.1, we find that the formulas for the upper boundary that correspond to (28) and (29) for the lower boundary are Here, λ refers to the partial derivative of λ with respect to boundary b. For small negative Δb in (33), whereas for small positive Δb in (34), so we have for both negative and positive Δb, Assuming that ϕ k is chosen so that ϕ k (a) > 0, the first factor in the expression for c km reduces to (−1) k−m when Δb approaches zero. Consider now a step in a immediately followed by a step in b, which is in turn followed by a period of constant boundary. Starting from the coefficient ω k (t), the step Δa in a leads to ω + k (t). The immediate subsequent change Δb in b then leads to ω ++ k (t). Finally, the boundaries remain constant for Δt, leading to ω k (t + Δt). In vector form, where we have defined C b as in (32), mutatis mutandis. Summing these, while observing that Dividing by Δt, and taking the limit when Δt approaches zero, The convergence is shown in the same way as in section 3.3.
We are now ready to summarize the results from sections 3.1-3.4 above into the following theorem: Theorem 1 (First-exit time for time-homogeneous processes). Let X t be a stochastic process defined by the Itô time-homogeneous stochastic differential equation Then the probability that X t remains within D at time t is where ω (ω 1 , ω 2 , . . .) T , ω k (0) = π(x)|ψ k (x, 0) , Λ diag({λ m }), k, m = 1, 2, . . .and Remarks. Here we have required, for convenience, that the boundary is continuously differentiable, but it is easy to see how the theorem extends to any regulated function boundary. When there are steps, equations (28), (29) and (33), (34) can be used to find the change in ω. Note that these equations are exact and do not involve any approximations. Border discontinuities translate directly to survivor function discontinuities. Equation (37) is handy for manipulation. For example, if we are interested in the larget or asymptotic properties of S(t), we can introduce the substitution s = log t which adds an exponential factor e s in front of Λ, Another example is when the boundaries converge to fixed values, implying that da/dt, db/dt → 0 when t → ∞ and that the equation approaches the well-known fixed-boundary equation dω dt = − Λω.

First-passage time for time-homogeneous processes
A solution for the first-passage problem can be obtained by making b infinite in theorem 1: Corollary 1 (First − passage time for time − homogeneous processes). Suppose that the conditions in theorem 1 hold, but that b is a natural boundary, i.e., b = ∞ and D = {(x, t)|t 0, a(t) x}. The same conclusion (35) can then be drawn, except that Remarks. The first-passage problem is more difficult than the first-exit problem because the corresponding Sturm-Liouville system is now singular, due to the unboundedness of the spatial interval (b = ∞). Therefore, the spectrum of the Sturm-Liouville operator is not necessarily discrete, and the solution cannot always be expressed as a series expansion in terms of eigenfunctions. However, the problem can sometimes be transformed into one that has a Sturm-Liouville operator with a discrete spectrum, as will be illustrated in the first-passage problem for Wiener processes below.

Hitting time for time-inhomogeneous processes
Time-homogeneous Itô processes are always continuous and Markovian, but this does not hold when the drift and diffusion depend on time, μ = μ(x, t) and σ = σ(x, t). In this section, we investigate how far we can push the moving-eigenvalue method in this more difficult case. As a consequence of the time-inhomogeneity, the eigenfunctions for constant boundaries also depend on time, ψ k = ψ k (x, t), and the eigenfunction expansion of p (19) becomes even for rectangular domains. Differentiating this with respect to t produces a new term compared to (20), Projecting this on ψ k , we obtain or, in vector form, where F is the matrix of Fourier coefficients of the eigenfunction time derivatives, Using the same procedure as in section 3, we find that the equation for moving boundaries, corresponding to (37), is Writing d dt proves convergence in the same way as in section 3.3 when ω T (F T + F + 2Λ)ω/|ω| 2 is locally bounded from below. The symmetric matrix F + F T can be expressed as

Wiener processes
We are now ready to demonstrate the method by application to the two most common Itô processes: Wiener processes in this section, and Ornstein-Uhlenbeck processes in the next section.

First-exit time
The SDE for the standard Wiener process is corresponding to the Fokker-Planck equation The eigenvalues of the Sturm-Liouville system −p /2 = λp where p(a) = p(b) = 0, the normalized eigenfunctions, their integrals, and the boundary sensitivity matrices are  curves become indistinguishable to the naked eye. Computation of the theoretical curve used the ten first eigenmodes.
We now compare computed solutions with the known explicit solution S 0 (t) = 1/(t + 1) for the initial probability p(x, 0) = (1 − x 2 )/2 exp (1 − x 2 )/2 and the square root boundaries ± √ t + 1. By comparing the exact survivor function with the solutions S n (t) for the system (37) truncated to the first n dimensions, we can see how the truncation error depends on n. In addition, using (41), we can study the asymptotic behaviour, which is infeasible with Monte Carlo methods. The result is shown in figure 7. The left column shows the agreement between S 0 and S n for n = 1, 2, 4, 8, at both linear and logarithmic scale. To illustrate the differences better, the right column shows the logarithms of the relative errors ΔS n /S 0 = (S 0 − S n )/S 0 . The results for n = 1 and n = 2 are identical because the problem is symmetric.

First-passage time
As example 9 showed, the spectrum corresponding to the first-passage problem for the Wiener process is not discrete, but we can use Doob's transformation (example 1) to convert the problem into a first-passage time problem for an Ornstein-Uhlenbeck process.

First-exit time
The Langevin equation (2) defining the Ornstein-Uhlenbeck process corresponds to the Fokker-Planck equation The right-hand side can be written in Sturm-Liouville form (5) as The Liouville transformation Y = p exp(x 2 /4) results in (14), which has the general solution where U and V are the Weber parabolic cylinder functions [80, chapter 19]. The general solution to (44) is thus and for the boundary conditions p(a) = p(b) = 0, we can express the eigenfunctions by and the kth eigenvalue λ k is the kth solution of the equation Y λ (b) = 0. Dividing ϕ k by the norm ϕ k (17), produces the normalized eigenfunctions where we used (12) in the last step and the fact that Y k (a) = W{U, V} = 2/π [80, section 19.4]. By observing that exp(−x 2 /4) = U (−1/2, x) is a solution of (14) corresponding to the eigenvalue λ 0 = 0, and applying lemma 1, we can compute the integrals Ψ k of the normalized eigenfunctions, By theorem 1, for k = m and ρ = a, b, the elements of the boundary sensitivity matrices are Examples of a few first-exit survivor functions for the standard Ornstein-Uhlenbeck process are shown in figure 8, using the same boundaries and other parameters as the examples in section 6.1.  The easiest way to obtain the normalized eigenfunctions and their integrals is to let b approach infinity in (46) and (47), leading to

First-passage time
The boundary sensitivity matrix is Examples of a few first-passage survivor functions for the standard Ornstein-Uhlenbeck process are shown in figure 9. Here, there is no upper boundary, whereas the lower boundaries In figure 10, we compare computed survivor functions S n with the known-exact survivor function given the exponential boundary a(t) = −e −t and the initial probability Here, the survivor function decays very rapidly, but the ODE solver is able to keep up until t ≈ 30 when the survivor function approaches 10 −20 .

Discussion and conclusions
The main contribution of this paper is the exact reduction of hitting-time problems to systems of ODEs expressed by a simple vector equation. This is in spite of the boundary time-variability, which hitherto has prevented a spectral approach. We have demonstrated the simplicity and effectiveness of the method by application and numerical solution of sample moving-boundary hitting-time problems involving Wiener and Ornstein-Uhlenbeck processes. The differential equation (37) can be solved quickly and easily by conventional numerical methods with about the same efficiency as the integral in (48). In addition, the form (37) is adequate for describing the properties of the solution when, for example, we have a set of firstpassage time measurements, a parametric model to fit, and want to compute the Cramér-Rao lower bound for the parameter estimates [82]. Approximations. Although theorem 1 is the primary result of the paper, one major application is as a basis for computation, so we shall briefly discuss some numerical aspects.
If the initial distribution p(x, 0) contains only the fundamental eigenmode, i.e., if ω k = 0 for k 2, such as after a long period of nearly constant boundary, and the boundary changes slowly, truncation to a two-or even one-dimensional system appears to suffice in many practical situations. When truncating to a single dimension, the survivor function can be approximated as where we recognize λ 1 as the hazard function. It is interesting to note that Bass and Erickson [40] found an asymptotic version of this formula for the first-exit time problem with symmetric increasing boundaries and diffusion processes using a completely different approach. During periods when the non-fundamental eigenmodes are significant, such as when the initial probability distribution has sharp peaks or discontinuities, or when the boundaries change steeply, we temporarily need to include more eigenmodes. The higher modes will decay quickly; we can see from (18) that the eigenvalues λ k for the first-exit problem tend to increase as k 2 . Possible alternatives to using high-order modes are interpolation, a Hopf-Cole transformation [10], or techniques for convergence acceleration such as deferred extrapolation to the limit [83].
It is common in applications that the initial distribution consists of a spike at a determinate position x 0 , or, formally, a Dirac distribution δ(x − x 0 ), which is not in L 2 w . However, we can approximate this arbitrarily well by a narrow pulse, as we have done for the upper solid curves in figures 6-9. Implementation. Solving the ODE (37) numerically is straightforward. Most of the effort lies in computing the eigenvalue functions λ k (a, b), which depend only on a and b, and are otherwise independent of time in the time-homogeneous case. A useful approach is first to tabulate λ k , using, for example, Prüfer transformations [84], and then to interpolate. Eigenfunctions are only needed for conversion of solutions between the original space and the transformed space, and typically do not expect very small computational error. For debugging, it is helpful to solve a model problem additionally by Monte Carlo simulation and the corresponding Fokker-Planck equation by finite element methods. Limitations. The most serious limitation of the proposed method is for the first-passage problem when the spectrum is not discrete. A related difficulty occurs when the initial distribution is tail heavy and is not in L 2 w . Sometimes a transformation can avoid the problem, such as Doob's transformation in the case of Wiener processes (example 1), but finding an appropriate transformation is an open problem, to the best of our knowledge. An alternative approach could be to approximate a continuous spectrum by a discrete one [85]. Generalizations. The proposed approach can be described as a Galerkin spectral method that remains single domain despite its non-rectangularity. The moving boundary just leads to some additional terms in the ODEs in the spectral representation. Our results suggest that the moving-eigenvalue method proposed here is applicable to other partial differential equations of parabolic type with non-rectangular domains, such as the Schrödinger equation, and with some modification, also with Neumann or mixed boundary conditions.