Monte-Carlo Simulation and Kernel Density Estimation of First passage time (original) (raw)

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() (wherek=1,2,3) consist:

• object an object inheriting from classsnssde1d, snssde2d andsnssde3d.
• boundary an expression of a constant or time-dependent boundary \(S(t)\).

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

• The expected value \(\text{E}(\tau_{S(t)})\), using the commandmean.
• The variance \(\text{Var}(\tau_{S(t)})\), using the command moment with order=2 andcenter=TRUE.
• The median \(\text{Med}(\tau_{S(t)})\), using the command Median.
• The mode \(\text{Mod}(\tau_{S(t)})\), using the command Mode.
• The quartile of \(\tau_{S(t)}\), using the command quantile.
• The maximum and minimum of \(\tau_{S(t)}\), using the commandmin and max.
• The skewness and the kurtosis of \(\tau_{S(t)}\), using the commandskewness and kurtosis.
• The coefficient of variation (relative variability) of \(\tau_{S(t)}\), using the commandcv.
• The central moments up to order \(p\) of \(\tau_{S(t)}\), using the commandmoment.
• The result summaries of the results of Monte-Carlo simulation, using the command summary.

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

• object an object inheriting from classfptsdekd() (where k=1,2,3).
• pdf probability density function Joint orMarginal.

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> set.seed(1234, kind = "L'Ecuyer-CMRG")
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 966 among 1000.

R> head(fpt1d$fpt, n = 5)

[1] 0.038850 0.042709 0.113998 0.066493 0.244749

The following statistical measures (S3 method) for classfptsde1d() 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’, usingdfptsde1d() function (hist=TRUE based ontruehist() function in MASS package)

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

Since fptdApprox andDiffusionRgqdpackages 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 theApprox.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.DiffProcpackage:

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.4491 5.8324 6.2194 6.0628 6.0337

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.45436
Variance          0.82529
Median            6.09538
Mode              6.00222
First quartile    5.93615
Third quartile    6.33197
Minimum           5.47910
Maximum           8.92931
Skewness          1.62490
Kurtosis          3.97197
Coef-variation    0.14075
3th-order moment  1.21824
4th-order moment  2.70530
5th-order moment  5.36133
6th-order moment 11.11820

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 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 theGQD.TIpassage() function in package DiffusionRgqdto 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.DiffProcpackage:

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] 1.9278 1.3564 2.7549 1.8239 1.5315

Monte-Carlo Statistics of F.P.T:
|T(S(t),X(t)) = inf{t >=  1 : X(t) >=  12 }
                        
Mean             2.20795
Variance         0.50426
Median           2.10384
Mode             1.45117
First quartile   1.55048
Third quartile   2.68149
Minimum          1.09881
Maximum          3.99568
Skewness         0.53741
Kurtosis         2.31695
Coef-variation   0.32162
3th-order moment 0.19243
4th-order moment 0.58915
5th-order moment 0.52001
6th-order moment 0.99788

By plotting the approximations (hist=TRUE based ontruehist() 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()

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> set.seed(1234, kind = "L'Ecuyer-CMRG")
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.13622 0.50310
2 0.14040 0.50428
3 0.13260 0.50229
4 0.14323 0.50805
5 0.13881 0.51709

The following statistical measures (S3 method) for classfptsde2d() 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})})\):

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.13470  0.50335
Variance         0.00019  0.00003
Median           0.13418  0.50319
Mode             0.13300  0.50272
First quartile   0.12511  0.49986
Third quartile   0.14343  0.50675
Minimum          0.09431  0.48489
Maximum          0.22080  0.52167
Skewness         0.38774 -0.00054
Kurtosis         4.33349  3.16471
Coef-variation   0.10167  0.01046
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 usingdfptsde2d() 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> set.seed(1234, kind = "L'Ecuyer-CMRG")
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.51850 0.024239 0.75841
2 0.52495 0.024041 0.80606
3 0.53565 0.021115 0.75086
4 0.53324 0.023411 0.74892
5 0.54566 0.023167 0.76537

The following statistical measures (S3 method) for classfptsde3d() 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})})\):

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.53198 0.02323  0.78404
Variance         0.00015 0.00000  0.00093
Median           0.53138 0.02321  0.78402
Mode             0.53097 0.02353  0.78322
First quartile   0.52328 0.02231  0.76216
Third quartile   0.53975 0.02407  0.80498
Minimum          0.49196 0.01912  0.68248
Maximum          0.56943 0.02792  0.88831
Skewness         0.17921 0.16558 -0.01077
Kurtosis         2.96629 2.99563  2.99648
Coef-variation   0.02298 0.05735  0.03897
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 usingdfptsde3d() 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()

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 (2024). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.9, 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.