First Passage Time Problems

Examples

Scalar FPT problem

Consider a time-homogeneous diffusion process with dynamics governed by the SDE: \[ dX_t = 0.5X_t(1-X_t^2)dt+0.5dB_t. \] Suppose we are interested in the first exit time of the process \(X_t\) from the strip \([-0.5,2]\). we can then calculate the survival probability surface for all initial values on the strip \([-0.5,2]\) by solving a backward PDE (Tuckwell and Wan 1984). In R this can be achieved using the method of lines by way of the MOL.passage() function:

library("DiffusionRimp")

mu  <- function(X){(0.5*X*(1-X^2))}
sig <- function(X) {0.5}
res <- MOL.passage(Xs = 0.5,t = 10, limits=c(-0.5, 2), N = 51, delt = 1/250)
  
persp(res$Xt, res$time, res$surface, col = 'white', xlab = 'X_s',
      ylab = 'Time (t)', zlab = 'Survival Prob.', border = NA, shade = 0.5,
      theta = 145, phi = 30, ticktype = 'detailed')

plot(res$density~res$time, col = '#222299', type='l', main = 'First passage time density')

#-------------------------------------------------------------------------------

By supplying MOL.passage() with an initial value, the function will calculate a numerical derivative of the survival probability surface in order to calculate the first passage time density function.

Bivariate FPT problem

In similar fashion, we can calculate the survival probability surface for a bivariate diffusion escaping from an enclosed region. Consider a bivariate diffusion \[ \begin{aligned} dX_t &=[0.5X_t(1 -X_t^2)+Y_t]dt + dB_t^{(1)}\\ dY_t &=[0.5Y_t(1 -Y_t^2)-X_t]dt + dB_t^{(2)}\\ \end{aligned} \] on the square region \([-2, 2]\times [-2, 2]\). Then we can calculate the evolution of the survival probability surface in R using:

#===============================================================================
# Bi-cubic diffusion with concentration in the even quadrants.
#===============================================================================
 
 # Define drift and diffusion terms.
 mu1   <- function(X,Y){0.5*(1 - X^2)*X - Y}
 mu2   <- function(X,Y){0.5*(1 - Y^2)*Y - X}
 sig11 <- function(X,Y){1}
 sig22 <- function(X,Y){1}

 # Parameters of the problem.
 Xs    <- 0.5      # Starting X-coordinate
 Ys    <- 0.5      # Starting Y-coordinate
 t     <- 10       # Final horizon time
 Xbar  <- c(-2, 2) # Limits in X-dim
 Ybar  <- c(-2, 2) # Limits in Y-dim
 Nodes <- 51       # How many nodes x nodes (incl. ends)
 delt  <- 1/250    # Time step size

 res <- BiMOL.passage(Xs, Ys, t, c(Xbar, Ybar), Nodes, delt)
 
 time.sequence <- c(0.5, 2, 4, 6, 8, 10)
 k = 0
 for(i in time.sequence)
 {
  k = k+1
  persp(res$Xt, res$Yt, res$surface[,,i/delt], col='white',
        xlab = 'State (X_0)', ylab = 'State (Y_0)', 
        zlab = 'Survival Probability', border = NA, 
        shade = 0.5, theta = 145 + 180*k/6)
 }

  plot(res$density~res$time, col = '#222299', type = 'l',
       main = 'First passage time density')

#-------------------------------------------------------------------------------

First escape from a non-square region

Consider a bivariate Ornstein-Uhlenbeck process:

\[ \begin{aligned} dX_t &=[0.5X_t(1 -X_t^2)+Y_t]dt + dB_t^{(1)}\\ dY_t &=[0.5Y_t(1 -Y_t^2)-X_t]dt + dB_t^{(2)}.\\ \end{aligned} \]

We can approximate the survival probability surface for the process \((X_t,Y_t)\) escaping from a unit circle centred at \((0.5,0.5)\) by passing an indicator function to the BiMOL.passage() function. Here, the indicator function assigns points of the discretization lattice values of 0 or 1 depending on whether a given point on the lattice is within the region of interest or not. In R:

# Define drift and diffusion terms:
mu1   <- function(X,Y){-1*X}
mu2   <- function(X,Y){-1*Y}
sig11 <- function(X,Y){0.5}
sig22 <- function(X,Y){0.5}

Nodes <- 101
delt  <- 1/500

# Define a region in the x-y plane:
region <- function(x,y){sqrt((x-0.5)^2+(y-0.5)^2)<1}
res    <- BiMOL.passage(0.5, 0.5, 10, c(-1,2,-1,2), Nodes, delt, Phi = region)

# Draw a nice plot of the evolution of the survival prob.:
library("colorspace")
colpal=function(n){rev(sequential_hcl(n,power=0.8,l=c(40,100)))}
time.sequence <- c(0.1, 0.25, 0.5, 1)
k = 0
for(i in time.sequence)
{
  k = k+1
  filled.contour(res$Xt, res$Yt, res$surface[,,i/delt], color.palette = colpal,
                main = paste0('Survival probability \n (t = ', i,')'),
                xlab = expression(X[0]), ylab = expression(Y[0]),
                plot.axis = {
                th =seq(0,1,1/100)
                lines(0.5+sin(2*pi*th),0.5+cos(2*pi*th))
                abline(h = 0, v = 0, lty = 'dashed')
                })

}

We can compare this to simulated first passage times as follows… NOTE: The simulation algorithm used here will typically over-estimate first passage times. Consequently, the frequency distribution will be slightly off/biased (See Giraudo and Sacerdote (1999)). See for example the effect of changing the stepsize.

sim  = function(N = 50000, plt = FALSE)
{
   set.seed(200707881)
   delt = 1/2000
   times = rep(0, N)
   X = rep(0.5, N)
   Y = rep(0.5, N)
   k = 1
   d = 0
   while(N >1)
   {
     X = X + mu1(X,Y)*delt + sig11(X,Y)*rnorm(N,sd = sqrt(delt))
     Y = Y + mu2(X,Y)*delt + sig22(X,Y)*rnorm(N,sd = sqrt(delt))
     d = d + delt
     wh = (sqrt((X-0.5)^2+(Y-0.5)^2)>=1)
     if(plt)
     {
       plot(X,Y,pch = 20,xlim = c(-1,2), ylim =c(-1,2))
                     th =seq(0,1,1/100)
                lines(0.5+sin(2*pi*th),0.5+cos(2*pi*th))
     }
     if(any(wh))
     {
        X = X[which(!wh)]
        Y = Y[which(!wh)]
        m = N - length(X)
        N = length(X)
        times[k:(k+m-1)] =  d
        k = k+m
     }
   }
   return(times)
}

sm = sim()

hist(sm, freq = FALSE, border = 'white', col= 'gray75',
     ylim = range(res$density), xlim = c(0, 6))
lines(res$density~res$time, col = '#222299', lwd = 2)

To get a visualization of the problem, you can run a small scale simulation with a plot using the command sim(50, TRUE).

References

Giraudo, Maria Teresa, and Laura Sacerdote. 1999. “An Improved Technique for the Simulation of First Passage Times for Diffusion Processes.” Communications in Statistics-Simulation and Computation 28 (4). Taylor & Francis: 1135–63.

Tuckwell, Henry C, and Frederic YM Wan. 1984. “First-Passage Time of Markov Process to Moving Barriers.” Journal of Applied Probability. JSTOR, 695–709.

Solving First Passage Time Problems Using the Method of Lines

Etienne A.D. Pienaar

February, 2016