Monte-Carlo Simulation and Kernel Density Estimation of First passage time

A.C. Guidoum1 and K. Boukhetala2

2020-05-04

The fptsdekd() functions

A new algorithm based on the Monte Carlo technique to generate the random variable FPT of a time homogeneous diffusion process (1, 2 and 3D) through a time-dependent boundary, order to estimate her probability density function.

Let \(X_t\) be a diffusion process which is the unique solution of the following stochastic differential equation:

\[\begin{equation}\label{eds01} dX_t = \mu(t,X_t) dt + \sigma(t,X_t) dW_t,\quad X_{t_{0}}=x_{0} \end{equation}\]

if \(S(t)\) is a time-dependent boundary, we are interested in generating the first passage time (FPT) of the diffusion process through this boundary that is we will study the following random variable:

\[ \tau_{S(t)}= \left\{ \begin{array}{ll} inf \left\{t: X_{t} \geq S(t)|X_{t_{0}}=x_{0} \right\} & \hbox{if} \quad x_{0} \leq S(t_{0}) \\ inf \left\{t: X_{t} \leq S(t)|X_{t_{0}}=x_{0} \right\} & \hbox{if} \quad x_{0} \geq S(t_{0}) \end{array} \right. \]

The main arguments to ‘random’ fptsdekd() (where k=1,2,3) consist:

The following statistical measures (S3 method) for class fptsdekd() can be approximated for F.P.T \(\tau_{S(t)}\):

The main arguments to ‘density’ dfptsdekd() (where k=1,2,3) consist:

Examples

FPT for 1-Dim SDE

Consider the following SDE and linear boundary:

\[\begin{align*} dX_{t}= & (1-0.5 X_{t}) dt + dW_{t},~x_{0} =1.7.\\ S(t)= & 2(1-sinh(0.5t)) \end{align*}\]

Generating the first passage time (FPT) of this model through this boundary: \[ \tau_{S(t)}= \inf \left\{t: X_{t} \geq S(t) |X_{t_{0}}=x_{0} \right\} ~~ \text{if} \quad x_{0} \leq S(t_{0}) \]

Set the model \(X_t\):

R> f <- expression( (1-0.5*x) )
R> g <- expression( 1 )
R> mod1d <- snssde1d(drift=f,diffusion=g,x0=1.7,M=1000,method="taylor")

Generate the first-passage-time \(\tau_{S(t)}\), with fptsde1d() function ( based on density() function in [base] package):

R> St  <- expression(2*(1-sinh(0.5*t)) )
R> fpt1d <- fptsde1d(mod1d, boundary = St)
R> fpt1d
Itô Sde 1D:
    | dX(t) = (1 - 0.5 * X(t)) * dt + 1 * dW(t)
    | t in [0,1].
Boundary:
    | S(t) = 2 * (1 - sinh(0.5 * t))
F.P.T:
    | T(S(t),X(t)) = inf{t >=  0 : X(t) >=  2 * (1 - sinh(0.5 * t)) }
    | Crossing realized 967 among 1000.
R> head(fpt1d$fpt, n = 5)
[1] 0.015654 0.034897 0.121027 0.027400 0.188879

The following statistical measures (S3 method) for class fptsde1d() can be approximated for the first-passage-time \(\tau_{S(t)}\):

R> mean(fpt1d)
R> moment(fpt1d , center = TRUE , order = 2) ## variance
R> Median(fpt1d)
R> Mode(fpt1d)
R> quantile(fpt1d)
R> kurtosis(fpt1d)
R> skewness(fpt1d)
R> cv(fpt1d)
R> min(fpt1d)
R> max(fpt1d)
R> moment(fpt1d , center= TRUE , order = 4)
R> moment(fpt1d , center= FALSE , order = 4)

The kernel density approximation of ‘fpt1d’, using dfptsde1d() function (hist=TRUE based on truehist() function in MASS package)

R> plot(dfptsde1d(fpt1d),hist=TRUE,nbins="FD")  ## histogramm
R> plot(dfptsde1d(fpt1d))              ## kernel density

Since fptdApprox and DiffusionRgqd packages can very effectively handle first passage time problems for diffusions with analytically tractable transitional densities we use it to compare some of the results from the Sim.DiffProc package.

fptsde1d() vs Approx.fpt.density()

Consider for example a diffusion process with SDE:

\[\begin{align*} dX_{t}= & 0.48 X_{t} dt + 0.07 X_{t} dW_{t},~x_{0} =1.\\ S(t)= & 7 + 3.2 t + 1.4 t \sin(1.75 t) \end{align*}\]

The resulting object is then used by the Approx.fpt.density() function in package fptdApprox to approximate the first passage time density:

R> require(fptdApprox)
R> x <- character(4)
R> x[1] <- "m * x"
R> x[2] <- "(sigma^2) * x^2"
R> x[3] <- "dnorm((log(x) - (log(y) + (m - sigma^2/2) * (t- s)))/(sigma * sqrt(t - s)),0,1)/(sigma * sqrt(t - s) * x)"
R> x[4] <- "plnorm(x,log(y) + (m - sigma^2/2) * (t - s),sigma * sqrt(t - s))"
R> Lognormal <- diffproc(x)
R> res1 <- Approx.fpt.density(Lognormal, 0, 10, 1, "7 + 3.2 * t + 1.4 * t * sin(1.75 * t)",list(m = 0.48,sigma = 0.07))

Using fptsde1d() and dfptsde1d() functions in the Sim.DiffProc package:

R> ## Set the model X(t)
R> f <- expression( 0.48*x )
R> g <- expression( 0.07*x )
R> mod1 <- snssde1d(drift=f,diffusion=g,x0=1,T=10,M=1000)
R> ## Set the boundary S(t)
R> St  <- expression( 7 + 3.2 * t + 1.4 * t * sin(1.75 * t) )
R> ## Generate the fpt
R> fpt1 <- fptsde1d(mod1, boundary = St)
R> head(fpt1$fpt, n = 5)
[1] 8.5774 8.3330 6.0229 6.2040 8.4028
R> summary(fpt1)

Monte-Carlo Statistics of F.P.T:
|T(S(t),X(t)) = inf{t >=  0 : X(t) >=  7 + 3.2 * t + 1.4 * t * sin(1.75 * t) }
                         
Mean              6.51486
Variance          0.92541
Median            6.11030
Mode              6.01900
First quartile    5.94570
Third quartile    6.36463
Minimum           5.43196
Maximum           8.93411
Skewness          1.44150
Kurtosis          3.33182
Coef-variation    0.14766
3th-order moment  1.28326
4th-order moment  2.85329
5th-order moment  5.45631
6th-order moment 11.07543

By plotting the approximations:

R> plot(res1$y ~ res1$x, type = 'l',main = 'Approximation First-Passage-Time Density', ylab = 'Density', xlab = expression(tau[S(t)]),cex.main = 0.95,lwd=2)
R> plot(dfptsde1d(fpt1,bw="bcv"),add=TRUE)
R> legend('topright', lty = c(1, NA), col = c(1,'#BBCCEE'),pch=c(NA,15),legend = c('Approx.fpt.density()', 'fptsde1d()'), lwd = 2, bty = 'n')
`fptsde1d()` vs `Approx.fpt.density()`

fptsde1d() vs Approx.fpt.density()

fptsde1d() vs GQD.TIpassage()

Consider for example a diffusion process with SDE:

\[\begin{align*} dX_{t}= & \theta_{1}X_{t}(10+0.2\sin(2\pi t)+0.3\sqrt(t)(1+\cos(3\pi t))-X_{t}) ) dt + \sqrt(0.1) X_{t} dW_{t},~x_{0} =8.\\ S(t)= & 12 \end{align*}\]

The resulting object is then used by the GQD.TIpassage() function in package DiffusionRgqd to approximate the first passage time density:

R> require(DiffusionRgqd)
R> G1 <- function(t)
+      {
+  theta[1] * (10+0.2 * sin(2 * pi * t) + 0.3 * prod(sqrt(t),
+  1+cos(3 * pi * t)))
+  }
R> G2 <- function(t){-theta[1]}
R> Q2 <- function(t){0.1}
R> res2 = GQD.TIpassage(8, 12, 1, 4, 1 / 100, theta = c(0.5))

Using fptsde1d() and dfptsde1d() functions in the Sim.DiffProc package:

R> ## Set the model X(t)
R> theta1=0.5
R> f <- expression( theta1*x*(10+0.2*sin(2*pi*t)+0.3*sqrt(t)*(1+cos(3*pi*t))-x) )
R> g <- expression( sqrt(0.1)*x )
R> mod2 <- snssde1d(drift=f,diffusion=g,x0=8,t0=1,T=4,M=1000)
R> ## Set the boundary S(t)
R> St  <- expression( 12 )
R> ## Generate the fpt
R> fpt2 <- fptsde1d(mod2, boundary = St)
R> head(fpt2$fpt, n = 5)
[1] 2.3610 3.4198 1.5236 3.2258 1.3588
R> summary(fpt2)

Monte-Carlo Statistics of F.P.T:
|T(S(t),X(t)) = inf{t >=  1 : X(t) >=  12 }
                        
Mean             2.18183
Variance         0.49703
Median           2.07252
Mode             1.48321
First quartile   1.55813
Third quartile   2.67935
Minimum          1.15626
Maximum          3.99508
Skewness         0.60216
Kurtosis         2.33319
Coef-variation   0.32313
3th-order moment 0.21100
4th-order moment 0.57639
5th-order moment 0.54952
6th-order moment 1.00384

By plotting the approximations (hist=TRUE based on truehist() function in MASS package):

R> plot(dfptsde1d(fpt2),hist=TRUE,nbins = "Scott",main = 'Approximation First-Passage-Time Density', ylab = 'Density', xlab = expression(tau[S(t)]), cex.main = 0.95)
R> lines(res2$density ~ res2$time, type = 'l',lwd=2)
R> legend('topright', lty = c(1, NA), col = c(1,'#FF00004B'),pch=c(NA,15),legend = c('GQD.TIpassage()', 'fptsde1d()'), lwd = 2, bty = 'n')
`fptsde1d()` vs `GQD.TIpassage()`

fptsde1d() vs GQD.TIpassage()

FPT for 2-Dim SDE’s

Assume that we want to describe the following Stratonovich SDE’s (2D):

\[\begin{equation}\label{eq016} \begin{cases} dX_t = 5 (-1-Y_{t}) X_{t} dt + 0.5 Y_{t} \circ dW_{1,t}\\ dY_t = 5 (-1-X_{t}) Y_{t} dt + 0.5 X_{t} \circ dW_{2,t} \end{cases} \end{equation}\]

and \[ S(t)=\sin(2\pi t) \]

Set the system \((X_t , Y_t)\):

R> fx <- expression(5*(-1-y)*x , 5*(-1-x)*y)
R> gx <- expression(0.5*y,0.5*x)
R> mod2d <- snssde2d(drift=fx,diffusion=gx,x0=c(x=1,y=-1),M=1000,type="str")

Generate the couple \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})})\), with fptsde2d() function::

R> St <- expression(sin(2*pi*t))
R> fpt2d <- fptsde2d(mod2d, boundary = St)
R> head(fpt2d$fpt, n = 5)
        x       y
1 0.13129 0.49686
2 0.13770 0.50351
3 0.12611 0.50800
4 0.13025 0.49508
5 0.12735 0.50360

The following statistical measures (S3 method) for class fptsde2d() can be approximated for the couple \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})})\):

R> mean(fpt2d)
R> moment(fpt2d , center = TRUE , order = 2) ## variance
R> Median(fpt2d)
R> Mode(fpt2d)
R> quantile(fpt2d)
R> kurtosis(fpt2d)
R> skewness(fpt2d)
R> cv(fpt2d)
R> min(fpt2d)
R> max(fpt2d)
R> moment(fpt2d , center= TRUE , order = 4)
R> moment(fpt2d , center= FALSE , order = 4)

The result summaries of the couple \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})})\):

R> summary(fpt2d)

Monte-Carlo Statistics for the F.P.T of (X(t),Y(t))
    | T(S(t),X(t)) = inf{t >=  0 : X(t) <=  sin(2 * pi * t) }
    |    And
    | T(S(t),Y(t)) = inf{t >=  0 : Y(t) >=  sin(2 * pi * t) }
                  T(S,X)  T(S,Y)
Mean             0.13314 0.50327
Variance         0.00016 0.00003
Median           0.13267 0.50317
Mode             0.12910 0.50334
First quartile   0.12443 0.49978
Third quartile   0.14144 0.50644
Minimum          0.09399 0.48584
Maximum          0.19188 0.52176
Skewness         0.21425 0.04935
Kurtosis         3.40307 3.11558
Coef-variation   0.09431 0.00998
3th-order moment 0.00000 0.00000
4th-order moment 0.00000 0.00000
5th-order moment 0.00000 0.00000
6th-order moment 0.00000 0.00000

The marginal density of \((\tau_{(S(t),X_{t})}\) and \(\tau_{(S(t),Y_{t})})\) are reported using dfptsde2d() function.

R> denM <- dfptsde2d(fpt2d, pdf = 'M')
R> plot(denM)

A contour and image plot of density obtained from a realization of system \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})})\).

R> denJ <- dfptsde2d(fpt2d, pdf = 'J',n=100)
R> plot(denJ,display="contour",main="Bivariate Density of F.P.T",xlab=expression(tau[x]),ylab=expression(tau[y]))
R> plot(denJ,display="image",main="Bivariate Density of F.P.T",xlab=expression(tau[x]),ylab=expression(tau[y]))

A \(3\)D plot of the Joint density with:

R> plot(denJ,display="persp",main="Bivariate Density of F.P.T",xlab=expression(tau[x]),ylab=expression(tau[y]))

Return to fptsde2d()

FPT for 3-Dim SDE’s

Assume that we want to describe the following SDE’s (3D): \[\begin{equation}\label{eq0166} \begin{cases} dX_t = 4 (-1-X_{t}) Y_{t} dt + 0.2 dB_{1,t}\\ dY_t = 4 (1-Y_{t}) X_{t} dt + 0.2 dB_{2,t}\\ dZ_t = 4 (1-Z_{t}) Y_{t} dt + 0.2 dB_{3,t} \end{cases} \end{equation}\]

with \((B_{1,t},B_{2,t},B_{3,t})\) are three correlated standard Wiener process: \[ \Sigma= \begin{pmatrix} 1 & 0.3 &-0.5\\ 0.3 & 1 & 0.2 \\ -0.5 &0.2&1 \end{pmatrix} \] and \[ S(t)=-1.5+3t \]

Set the system \((X_t , Y_t , Z_t)\):

R> fx <- expression(4*(-1-x)*y , 4*(1-y)*x , 4*(1-z)*y) 
R> gx <- rep(expression(0.2),3)
R> Sigma <-matrix(c(1,0.3,-0.5,0.3,1,0.2,-0.5,0.2,1),nrow=3,ncol=3)
R> mod3d <- snssde3d(drift=fx,diffusion=gx,x0=c(x=2,y=-2,z=0),M=1000,corr=Sigma)

Generate the triplet \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})},\tau_{(S(t),Z_{t})})\), with fptsde3d() function::

R> St <- expression(-1.5+3*t)
R> fpt3d <- fptsde3d(mod3d, boundary = St)
R> head(fpt3d$fpt, n = 5)
        x        y       z
1 0.52386 0.022183 0.80782
2 0.52280 0.025187 0.85021
3 0.51602 0.023298 0.80947
4 0.54177 0.023816 0.78090
5 0.53506 0.024629 0.77339

The following statistical measures (S3 method) for class fptsde3d() can be approximated for the triplet \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})},\tau_{(S(t),Z_{t})})\):

R> mean(fpt3d)
R> moment(fpt3d , center = TRUE , order = 2) ## variance
R> Median(fpt3d)
R> Mode(fpt3d)
R> quantile(fpt3d)
R> kurtosis(fpt3d)
R> skewness(fpt3d)
R> cv(fpt3d)
R> min(fpt3d)
R> max(fpt3d)
R> moment(fpt3d , center= TRUE , order = 4)
R> moment(fpt3d , center= FALSE , order = 4)

The result summaries of the triplet \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})},\tau_{(S(t),Z_{t})})\):

R> summary(fpt3d)

Monte-Carlo Statistics for the F.P.T of (X(t),Y(t),Z(t))
    | T(S(t),X(t)) = inf{t >=  0 : X(t) <=  -1.5 + 3 * t }
    |    And
    | T(S(t),Y(t)) = inf{t >=  0 : Y(t) >=  -1.5 + 3 * t }
    |    And
    | T(S(t),Z(t)) = inf{t >=  0 : Z(t) <=  -1.5 + 3 * t }
                  T(S,X)  T(S,Y)   T(S,Z)
Mean             0.53109 0.02321  0.78215
Variance         0.00014 0.00000  0.00100
Median           0.53083 0.02316  0.78309
Mode             0.53106 0.02300  0.78850
First quartile   0.52263 0.02229  0.76112
Third quartile   0.53910 0.02411  0.80285
Minimum          0.49262 0.01888  0.67295
Maximum          0.56859 0.02869  0.86988
Skewness         0.12203 0.08678 -0.15598
Kurtosis         2.89278 3.08794  3.08544
Coef-variation   0.02233 0.05763  0.04036
3th-order moment 0.00000 0.00000  0.00000
4th-order moment 0.00000 0.00000  0.00000
5th-order moment 0.00000 0.00000  0.00000
6th-order moment 0.00000 0.00000  0.00000

The marginal density of \(\tau_{(S(t),X_{t})}\) ,\(\tau_{(S(t),Y_{t})}\) and \(\tau_{(S(t),Z_{t})})\) are reported using dfptsde3d() function.

R> denM <- dfptsde3d(fpt3d, pdf = "M")
R> plot(denM)

For an approximate joint density for \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})},\tau_{(S(t),Z_{t})})\) (for more details, see package sm or ks.)

R> denJ <- dfptsde3d(fpt3d,pdf="J")
R> plot(denJ,display="rgl")

Return to fptsde3d()

Further reading

  1. snssdekd() & dsdekd() & rsdekd()- Monte-Carlo Simulation and Analysis of Stochastic Differential Equations.
  2. bridgesdekd() & dsdekd() & rsdekd() - Constructs and Analysis of Bridges Stochastic Differential Equations.
  3. fptsdekd() & dfptsdekd() - Monte-Carlo Simulation and Kernel Density Estimation of First passage time.
  4. MCM.sde() & MEM.sde() - Parallel Monte-Carlo and Moment Equations for SDEs.
  5. TEX.sde() - Converting Sim.DiffProc Objects to LaTeX.
  6. fitsde() - Parametric Estimation of 1-D Stochastic Differential Equation.

References

  1. Boukhetala K (1996). Modelling and Simulation of a Dispersion Pollutant with Attractive Centre, volume 3, pp. 245-252. Computer Methods and Water Resources, Computational Mechanics Publications, Boston, USA.

  2. Boukhetala K (1998). Estimation of the first passage time distribution for a simulated diffusion process. Maghreb Mathematical Review, 7, pp. 1-25.

  3. Boukhetala K (1998). Kernel density of the exit time in a simulated diffusion. The Annals of The Engineer Maghrebian, 12, pp. 587-589.

  4. Guidoum AC, Boukhetala K (2020). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.6, URL https://cran.r-project.org/package=Sim.DiffProc.

  5. Pienaar EAD, Varughese MM (2016). DiffusionRgqd: An R Package for Performing Inference and Analysis on Time-Inhomogeneous Quadratic Diffusion Processes. R package version 0.1.3, URL https://CRAN.R-project.org/package=DiffusionRgqd.

  6. Roman, R.P., Serrano, J. J., Torres, F. (2008). First-passage-time location function: Application to determine first-passage-time densities in diffusion processes. Computational Statistics and Data Analysis. 52, 4132-4146.

  7. Roman, R.P., Serrano, J. J., Torres, F. (2012). An R package for an efficient approximation of first-passage-time densities for diffusion processes based on the FPTL function. Applied Mathematics and Computation, 218, 8408-8428.

  1. Department of Probabilities & Statistics, Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, E-mail (acguidoum@usthb.dz)↩

  2. Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, E-mail (kboukhetala@usthb.dz)↩