First-passage Time Estimation of Di ff usion Processes through Time-Varying Boundaries with an Application in Finance

In this paper, we develop a Monte Carlo based algorithm for estimating the FPT (first passage time) density of the solution of a one-dimensional time-homogeneous SDE (stochastic differential equation) through a time-dependent frontier. We consider Brownian bridges as well as local Daniels curve approximations to obtain tractable estimations of the FPT probability between successive points of a simulated path of the process. Under mild assumptions, a (unique) Daniels curve local approximation can easily be obtained by explicitly solving a non-linear system of equations.


Introduction
Let X be a time homogeneous one-dimensional diffusion process which is the unique (strong) solution of the following stochastic differential equation: that is the functions µ and σ satisfy regularity conditions as described for example in Karatzas and Shreve (1991).
If S is a time dependent boundary, we are interested in estimating either the pdf or cdf of the first-passage time of the diffusion process through this boundary that is we will study the following random variable: In general, there is no explicit expression for the first-passage time density of a diffusion process through a time-varying boundary.To this date, few specific cases provide closed-form formulas for some classes of time-varying boundaries for example when the process is Gaussian (Durbin (1985), Durbin and Williams (1992), DiNardo et al. (2001)) or a Bessel process (Deaconu and Herrmann (2013)).Thus, we mainly rely on simulation techniques to estimate this density for the more general case.Besides standard Monte Carlo estimation, only a few papers suggest alternatives or improvements such as considering upcrossing probabilities in between simulation points (Giraudo and Sacerdote (1999), Giraudo, Sacerdote, and Zucca (2001)).
The main goal of this work is to construct a computationally efficient algorithm that will provide reliable FPT density and probability estimates.The paper is organized as follows.In section 2, we review existing Monte Carlo based techniques.
In section 3, we establish the mathematical foundations leading to a novel algorithm.Section 4 is devoted to various examples enabling us to evaluate the algorithm's performance.Section 5 consists of an empirical comparison study of Monte Carlo based methods.Finally, section 6 illustrates the pertinence of FPT estimation in the context of portfolio management.

Crude Monte Carlo
This is the simplest and best-known approach based on the law of large numbers.After fixing a time interval, basically we divide the latter into smaller ones, simulate a path of the process along those time points and, if it occurs, note the subinterval where the first upcrossing occurs.Generally, the midpoint of this subinterval forms the estimated first-passage time of this simulated path.We repeat the process a large number of time in order to construct a pdf or cdf estimate of this stopping time.
Usually, we need sizeable amounts of simulations and an extremely fine partition of time intervals in order to have a suitable estimation of first-passage time probabilities.Another drawback of the crude Monte Carlo approach is that it tends to overestimates the true value of the first-passage time since an upcrossing may occur earlier in between simulated points of a complete path as illustrated in figure 1.

Monte Carlo Approach with Upcrossings
Instead of continuously repeating the whole Monte Carlo procedure with an even finer interval partition to obtain better estimates, let us see how one could improve on the initial estimates without discarding the original simulated paths.
An astute idea that have been put forward is to ideally obtain the probability law of an upcrossing between simulated points, thus if p k is the probability of an upcrossing in the time interval [t k , t k+1 ], then one would simply generate a value U k taken from a uniform random variable on [0, 1] and assert that there is an upcrossing if U k p k .
Unfortunately the FPT of a diffusion bridge will more than often not be available for more general diffusion processes, thus we need to propose an adequate estimation of this probability.
For each subinterval, one could simulate paths of tied-down processes to obtain FPT estimates for each subintervals along simulated paths of the process.This approach was used for example in Giraurdo and Sacerdote (1999) and Giraurdo, Sacerdote and Zucca (2001) where Kloeden-Platen numerical schemes were used for path simulations to construct approximation of diffusion bridges.
Another alternative, as first proposed in Strittmatter (1987) is to consider that for a small enough interval, the diffusion part of the process should remain fairly constant and then look at a Brownian bridge approximation of the diffusion bridge and exploit known results on the FPT of Brownian bridges.
For the Brownian bridge to constitute an efficient approximation, the diffusion part σ(X t ) of the SDE (stochastic differential equation) should be constant.Therefore one should apply a Lamperti transform on both the original process and the frontier as described for example in Iacus (2008).Indeed , define and apply it on the original process and the time-varying boundary.Assuming that F is one-to-one, then the original problem is equivalent to finding the FPT density of where the new boundary is given by S * (t) = F (S (t)) and, by Itô's formula, the diffusion process Y follows the dynamics ) Another advantage of opting for the Lamperti transformed FPT problem in this form is that one can readily make use of results from Downes and Borovkov (2008) on bounds on the FPT pdf.
Following the Brownian bridge approximation, we need to compute the upcrossing probabilities through the new boundary S * .Dating back to Durbin's (1985) paper, we know that the FPT density of continuous Gaussian processes through a timevarying boundary must satisfy a Volterra type integral equation of the second kind.Oftentimes, we cannot explicitly solve this equation and therefore we must rely on some numerical approximation.For example, Durbin and Williams (1992) proposed an expansion series where each term involve computing multiple integrals while Di Nardo et al. (2001) proposed a more tractable iterative scheme to numerically solve the equation.

Local Daniels Curve Approximation Approach
While still favoring Brownian bridge approximations of the diffusion bridges after a Lamperti transform as described previously, instead of attempting to approximate the FPT of the brownian bridge through the original boundary using a costly computational method, we could instead introduce carefully chosen successive local curve approximations of the given frontier.
A first attempt would be to consider piecewise linear approximations of the time-varying boundary since explicit formulas are readily available in this case.Although computationally fast, since it does not take into account local curvatures, this technique may sometimes underestimate (locally convex) or overestimate (locally concave) the Brownian bridge FPT upcrossing probability and therefore should be used for a fine partition of the time interval.
In our approach, we propose to focus on local Daniels curve approximations of the time-varying boundary.Indeed, the more flexible three-parameter Daniels curve of the form also produces explicit tractable formulas of the first passage time probability, thus one would typically get a better approximation of the true probability p k .
First we must show, under mild assumptions, that a (unique) Daniels curve approximation can easily be obtained by simply taking the endpoints of the segment and the value at midpoint (or another point of our choosing).Thus the following proposition will provide useful.
Proposition 1.Let [0, ∆t] be a time interval, consider the points (0, a), (∆t/2, b),(∆t, c) and set A = e then there is a unique Daniels curve (4) passing through the three points with parameters Proof.The set of points generate the following non-linear system of equations: obviously the first equation gives α = 2a > 0, while simple algebraic manipulations on the last two equations lead us to solve the following linear system which can be rewritten in the form ) > 0 then there exists a unique solution given by: this would constitute the solution to the original system provided that β > 0 and lim t→∆t β 2 + 4γe −α 2 /t > 0.
Notice first that 1 A then γ 0 and clearly lim t→∆t β 2 + 4γe −α 2 /t > 0 is satisfied.So if we assume now that B C > 1 A then γ < 0 , thus we need to verify that β 2 + 4γ A 2 > 0 which is the case since The final step is to make sure that it solves the original system.Substituting back in system ( 7), (where only positive square roots are involved), we see that is the case only if which is verified through (5).
Next, by applying Theorem 3.4 of Di Nardo et al. (2001) to the special case of Brownian bridges, the following proposition provides us with an explicit expression for the FPT upcrossing probabilities.
Proposition 2. Consider a Brownian bridge W 0,∆t define on an time interval [0, ∆t] and S a Daniels curve defined by ( 4) where α, β > 0, γ ∈ R and Finally, proposition 1 and 2 suggest that once a path of the diffusion process is produced, for each pair of simulated points on a subinterval of length ∆t, one simply applies a linear translation of the line segment between these points onto the interval [0, ∆t] while translating the corresponding boundary segment.
This naturally leads us to the following algorithm:  Note that step 9 includes extreme cases where the middle point of the frontier in a subinterval may not be reached by a Daniels curve, thus we use the closest curve possible.

Numerical Examples
We will focus our examples on diffusion processes which paths can be simulated exactly.Therefore with known results on FPT density and bounds, it will allow us to better visualize the approximation error due essentially to the algorithm.
Example 1 Consider the following Ornstein-Uhlenbeck process and time varying boundary This diffusion process is a Gauss-Markov process and according to Di Nardo et al. ( 2001) the chosen boundary allows us to obtain an explicit FPT density given by where φ 0 is the probability density function of the Ornstein-Uhlenbeck process starting at X (0) = 1.6.
Figure 2 compares the true FPT density with the empirical density histogram obtained through our algorithm using a time step discretization of 0.01 of the time interval [0, 1] and 10 000 simulated paths.Furthermore, the algorithm gives us a FPT probability estimate of 0.9619 over the whole interval compared to the true value of 0.9608 representing a relative error of about 0.12%.Example 2 Consider the following squared Bessel process and time-varying boundary By applying the Lamperti transform to both the process and boundary we obtain respectively This transformed diffusion process is the Bessel process and, according to Deaconu and Herrmann (2013) the chosen boundary allows us to obtain an explicit FPT density given by.
Figure 3 compares the true FPT density with the empirical density histogram obtained through our algorithm using a time step discretization of 0.01 of the time interval [0, 1] and 10000 simulated paths.Furthermore, the algorithm gives us a FPT probability estimate of 0.7432 over the whole interval compared to the true value of 0.7500 representing a relative error of about 0.91%.) In this case, it is worth noting that since the true density does not satisfy the initial condition S (0) > X (0), therefore to increase the algorithm's accuracy near the origin one should envisage a finer partition around this point.
Example 3 Consider the following geometric Brownian process and linear boundary By applying the Lamperti transform to both the process and boundary we obtain respectively As in example 1, this transformed diffusion process is a Gauss-Markov process and, although the new frontier does not allow an explicit FPT density, using the deterministic algorithm in Di Nardo et al. ( 2001) with a 0.01 time step discretization of the time interval [0, 1], we can obtain a reliable approximation.
Figure 4 compares the Di Nardo FPT density approximation with the empirical density histogram obtained through our algorithm using the same time step discretization with 10 000 simulated paths.
By applying the Lamperti transform to both the process and boundary we obtain respectively As opposed to the preceding examples, this diffusion process has neither an explicit expression nor probability law however using Beskos and Roberts' ( 2005) exact algorithm we can simulate exact sample paths.Although an explicit FPT density is not available, using results of Downes and Borovkov (2008) we can, in this case, obtain the following lower and upper bounds: φ W (S (t)) Figure 5 compares the FPT bounds with the empirical density histogram obtained through our algorithm using a 0.01 time step discretization of the time interval [0, 1] starting initially with 15000 simulations and obtaining 11768 valid paths through the exact algorithm.

Empirical Efficiency Comparison of Monte Carlo Based Techniques
As of now, essentially three different approaches exist in estimating the FPT probabilities of a general time-homogeneous one-dimensional diffusion process through a time-varying boundary, all of them are Monte Carlo based and involve estimating upcrossing probabilities of a diffusion bridge in between simulation points.
While Giraurdo, Sacerdote and Zucca (2001) suggest Monte Carlo subsimulations of these diffusion bridges, Giraurdo and Sacerdote (1999) and the present authors propose first a Brownian bridge approximation then followed by different deterministic numerical estimations of the upcrossing probabilites.
We will focus our comparison study on the total FPT probability of the Ornstein-Uhlenbeck process of our first numerical example.This choice is mainly motivated by the fact that both the diffusion process and associated diffusion bridges can be simulated exactly and furthermore the FPT pdf of the diffusion process has a known closed-form expression.
Since here we are only interested in comparing the total FPT probability on the time interval [0,1], once the original simulated paths are drawn only paths for which no barrier crossings were detected will require upcrossing probability estimations.

Comparison with Simulated Diffusion Bridges Approach
Basically, once standard Monte Carlo path simulations are drawn, we fix the number m of time steps discretization of all subintervals as well as the number M of simulated diffusion bridge paths in each subinterval.In Giraurdo, Sacerdote and Zucca's (2001) paper they devise an elaborate scheme for a general diffusion bridge through Kloeden-Platen approximation schemes.In our case an Ornstein-Uhlenbeck bridge can directly be simulated (see Bladt and Sorensen (2014)), thus we will use the later since we will gain in accuracy as well as computing time.Using N = 10000 simulated paths with time steps of length ∆ = 0.01, we obtained the following results: We notice that while roughly both methods give us the same relative errors for the FPT estimate, our method produces a reasonable estimate at least four times faster.Although the Ornstein-Uhlenbeck bridge simulations would give us eventually the true upcrossing probabilities with larger values of M and m it seems there would be a heavy price to pay in regards to computing time.

Comparison with FPT Fdf Approximation of Brownian Bridges Approach
Again we start with the usual Monte Carlo path simulations but, following the Brownian bridge approach, we need a FPT probabilty approximation of an upcrossing in between simulation points.Giraurdo and Sacerdote (1999) proposed an iterative scheme to produce m estimated points of the true FPT pdf in each subinterval, it is based on approximating the Volterra integral equation associated with the FPT pdf.As a standard for comparison, since our method makes use of an estimation based on the values of only three points taken from the boundary in each subinterval, we will do the same for their method.Therefore three-point Simpson's rules will estimate the upcrossing probabilities.Using each time N = 10000 simulated paths, we obtained the following results: Here we achieved a relative error of less than 0.2% as fast as about one second of computing time while the three-point FPT pdf approximation needed as much as 80 seconds to achieve the same level of accuracy.Of course one could increase the number of FPT pdf estimation points in order to achieve better accuracy but the immediate drawback seems to be generating sensibly greater computing time.

Application in Portfolio Management
Consider a risky strategy on an horizon [0, T ], the investor may encounter a specific instant t when the amount of wealth x(t) be sufficient enough so that he may, at this point, safely reinvest all of his money in a simple bank account with constant interest rate r and the resulting terminal wealth x(T ) will attain his financial goal z.So we consider the following and we naturally want to compute the probability P(τ z ≤ T ) of such an event.If x 0 > 0 is his initial wealth then we will assume that z > x 0 e rT so that the investor cannot achieve his financial goal by simply placing his initial investment in a bank account.
In order to investigate this goal-achieving problem, we must first define a mathematical setting for the dynamics of the financial market.We will consider here the celebrated Black-Scholes model that we next describe.The first asset is a bank account whose price at time t, S 0 (t), is the solution to the following ODE (ordinary differential equation): The next asset consist of a stock whose price S 1 (t) at time t are the solutions to the following SDE (stochastic differential equation): where {W(t), t≥0} is a standard Brownian motion and b > r.
Let u(t), 0 ≤ t ≤ T be a financial strategy (or portfolio) where u(t) is the amount placed in the stock.If we assume that all strategies u(t) are self-financed (no outside injection of funds to the investors) and with no transaction costs then the wealth dynamic at time t is given by the following SDE: dx (t) = (rx (t) + σθu (t)) dt + u (t) σdW (t) where θ = b−r σ .Finally, among all the possible strategies, we will focus on the one generated by a family of stochastic control problems defined by min V AR (x (T )) subject to E (x (T )) = z These are known as mean-variance problems and are considered the cornerstone of modern portfolio management.
In the case where the portfolio is unconstrained, the optimal wealth process is given by: Clearly we conclude in our example that, in the case of mean-variance strategies under a no-bankruptcy restriction, the FPT probability is a decreasing function of the targeted wealth.

Conclusion
In this paper, we proposed an innovative Monte Carlo based algorithm for estimating the first-passage time density of a one-dimensional diffusion process through a deterministic time-dependent boundary.After considering a Lamperti transform, the main idea relied on a piecewise Daniels curve approximation of the boundary as well as Brownian bridge approximations of diffusion bridges.Through specific examples we illustrated the efficiency of our method by comparison with either known theoretical results or an existing deterministic method.We followed with an empirical comparison of accuracy and speed of our method with other Monte Carlo based methods.We concluded with a practical application pertaining to a portfolio management problem.Finally, an interesting question would be to see if our methodology could be extended to other types of processes with continuous paths such as superdiffusion and subdiffusion processes.

Figure 1 .
Figure 1.Undetected prior upcrossing through a basic Monte Carlo path simulation

[
Step 1:] Apply the Lamperti transform (2) to the original diffusion process (1) and frontier S to obtain the new process (3) and boundary S * = F (S (t)) [Step 2:] Select a time interval [T l , T u ] and construct a partition T l = t 0 < t 1 < . . .< t n = T u [Step 3:] Initialize FPT vector counter to τ := {0, . . ., 0} [Step 4:] Initialize path counter to k := 1 WHILE k is less than N the number of desired paths DO the following: [Step 5:] Simulate a path of the process {Y (t 1 ) , . . ., Y (t n )} [Step 6:] Initialize subinterval counter to i := 1 WHILE i is less than n the number of desired subintervals DO the following: [Step 7:] IF Y (t i ) S * (t i ) THEN set i th FPT vector component to τ i := τ i + 1 and path counter to k := k + 1, GO TO Step 5 [Step 8:] Set ∆

∆[
Step 10:] Set probability upcrossing to p := c 1 e − α 2 2∆ + c 2 e − 2α 2 Step 11:] Generate a value U taken from a uniform random variable [Step 12:] IF U ≤ p THEN set i th FPT vector component to τ i := τ i + 1 and path counter to k := k + 1, GO TO Step 5, ELSE set i := i + 1, GO TO Step 7

Figure 4 .
Figure 4. Geometric Brownian Motion FPT pdf through the boundary 1.0 + 2.0t Example 4 Consider the modified Cox-Ingersoll-Ross process and linear boundary

xFigure 6 .
Figure 6.Comparison of goal-achieving FPT probabilities of mean-variance strategies