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

A.C. Guidoum1 and K. Boukhetala2

2017-09-16

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=10000,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 9709 among 10000.
R> head(fpt1d$fpt, n = 10)
 [1] 0.055751 0.335977 0.429426 0.270945 0.065990 0.073951 0.058403
 [8] 0.677860 0.220476 0.132349

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

R> mean(fpt1d)
[1] 0.19478
R> moment(fpt1d , center = TRUE , order = 2) ## variance
[1] 0.041274
R> Median(fpt1d)
[1] 0.1151
R> Mode(fpt1d)
[1] 0.049867
R> quantile(fpt1d)
       0%       25%       50%       75%      100% 
0.0056766 0.0545646 0.1150954 0.2599597 0.9996590 
R> kurtosis(fpt1d)
[1] 5.6836
R> skewness(fpt1d)
[1] 1.745
R> cv(fpt1d)
[1] 1.0431
R> min(fpt1d)
[1] 0.0056766
R> max(fpt1d)
[1] 0.99966
R> moment(fpt1d , center= TRUE , order = 4)
[1] 0.0096844
R> moment(fpt1d , center= FALSE , order = 4)
[1] 0.031922

The result summaries of the first-passage-time \(\tau_{S(t)}\):

R> summary(fpt1d)

Monte-Carlo Statistics of F.P.T:
|T(S(t),X(t)) = inf{t >=  0 : X(t) >=  2 * (1 - sinh(0.5 * t)) }
                        
Mean             0.19478
Variance         0.04128
Median           0.11510
Mode             0.04987
First quartile   0.05456
Third quartile   0.25996
Minimum          0.00568
Maximum          0.99966
Skewness         1.74498
Kurtosis         5.68361
Coef-variation   1.04306
3th-order moment 0.01463
4th-order moment 0.00968
5th-order moment 0.00594
6th-order moment 0.00396

Display the exact first-passage-time \(\tau_{S(t)}\), see Figure 1:

R> plot(time(mod1d),mod1d$X[,1],type="l",lty=3,ylab="X(t)",xlab="time",axes=F)
R> curve(2*(1-sinh(0.5*x)),add=TRUE,col=2)
R> points(fpt1d$fpt[1],2*(1-sinh(0.5*fpt1d$fpt[1])),pch=19,col=4,cex=0.5)
R> lines(c(fpt1d$fpt[1],fpt1d$fpt[1]),c(0,2*(1-sinh(0.5*fpt1d$fpt[1]))),lty=2,col=4)
R> axis(1, fpt1d$fpt[1], bquote(tau[S(t)]==.(fpt1d$fpt[1])),col=4,col.ticks=4)
R> legend('topleft',col=c(1,2,4),lty=c(1,1,NA),pch=c(NA,NA,19),legend=c(expression(X[t]),expression(S(t)),expression(tau[S(t)])),cex=0.8,bty = 'n')
R> box()

The kernel density approximation of ‘fpt1d’, using dfptsde1d() function (hist=TRUE based on truehist() function in MASS package), see e.g. Figure 2.

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=10000)
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> fpt1
Itô Sde 1D:
    | dX(t) = 0.48 * X(t) * dt + 0.07 * X(t) * dW(t)
    | t in [0,10].
Boundary:
    | S(t) = 7 + 3.2 * t + 1.4 * t * sin(1.75 * t)
F.P.T:
    | T(S(t),X(t)) = inf{t >=  0 : X(t) >=  7 + 3.2 * t + 1.4 * t * sin(1.75 * t) }
    | Crossing realized 10000 among 10000.
R> head(fpt1$fpt, n = 10)
 [1] 5.8433 5.9046 8.5778 5.8437 6.3338 8.2663 6.0841 6.0792 5.9392
[10] 5.9133
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.50656
Variance          0.90314
Median            6.11345
Mode              6.02315
First quartile    5.94489
Third quartile    6.37952
Minimum           5.33516
Maximum           9.00606
Skewness          1.48294
Kurtosis          3.48450
Coef-variation    0.14606
3th-order moment  1.27278
4th-order moment  2.84217
5th-order moment  5.50503
6th-order moment 11.29587

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=10000)
R> ## Set the boundary S(t)
R> St  <- expression( 12 )
R> ## Generate the fpt
R> fpt2 <- fptsde1d(mod2, boundary = St)
R> fpt2
Itô Sde 1D:
    | dX(t) = theta1 * 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)
    | t in [1,4].
Boundary:
    | S(t) = 12
F.P.T:
    | T(S(t),X(t)) = inf{t >=  1 : X(t) >=  12 }
    | Crossing realized 9239 among 10000.
R> head(fpt2$fpt, n = 10)
 [1] 2.6047 1.3890 2.0845 2.3062 2.1464 1.2869 2.7639 1.6551 2.3065
[10] 1.9044
R> summary(fpt2)

Monte-Carlo Statistics of F.P.T:
|T(S(t),X(t)) = inf{t >=  1 : X(t) >=  12 }
                        
Mean             2.17508
Variance         0.48925
Median           2.07380
Mode             1.42680
First quartile   1.53478
Third quartile   2.66124
Minimum          1.09920
Maximum          3.99968
Skewness         0.59228
Kurtosis         2.40614
Coef-variation   0.32158
3th-order moment 0.20269
4th-order moment 0.57596
5th-order moment 0.55317
6th-order moment 1.02681

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

The following \(2\)-dimensional SDE’s with a vector of drift and a diagonal matrix of diffusion coefficients:

\[\begin{equation}\label{eq:09} \begin{cases} dX_t = f_{x}(t,X_{t},Y_{t}) dt + g_{x}(t,X_{t},Y_{t}) dW_{1,t}\\ dY_t = f_{y}(t,X_{t},Y_{t}) dt + g_{y}(t,X_{t},Y_{t}) dW_{2,t} \end{cases} \end{equation}\]

\(W_{1,t}\) and \(W_{2,t}\) is a two independent standard Wiener process. First passage time (2D) \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})})\) is defined as:

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

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=10000,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> fpt2d
Stratonovich Sde 2D:
    | dX(t) = 5 * (-1 - Y(t)) * X(t) * dt + 0.5 * Y(t) o dW1(t)
    | dY(t) = 5 * (-1 - X(t)) * Y(t) * dt + 0.5 * X(t) o dW2(t)
    | t in [0,1].
Boundary:
    | S(t) = sin(2 * pi * t)
F.P.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) }
    | Crossing realized 10000 among 10000.
R> head(fpt2d$fpt, n = 10)
         x       y
1  0.14163 0.50239
2  0.15048 0.49761
3  0.14318 0.50680
4  0.13822 0.49820
5  0.11952 0.50771
6  0.13421 0.50393
7  0.14127 0.50701
8  0.12712 0.49730
9  0.13485 0.50973
10 0.14422 0.50425

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)
[1] 0.13380 0.50332
R> moment(fpt2d , center = TRUE , order = 2) ## variance
[1] 0.000171600 0.000027741
R> Median(fpt2d)
[1] 0.13350 0.50333
R> Mode(fpt2d)
[1] 0.13444 0.50385
R> quantile(fpt2d)
$x
      0%      25%      50%      75%     100% 
0.078654 0.124879 0.133500 0.141967 0.203225 

$y
     0%     25%     50%     75%    100% 
0.48107 0.49981 0.50333 0.50674 0.52376 
R> kurtosis(fpt2d)
[1] 3.3049 3.2157
R> skewness(fpt2d)
[1] 0.266992 0.047234
R> cv(fpt2d)
[1] 0.097910 0.010465
R> min(fpt2d)
[1] 0.078654 0.481071
R> max(fpt2d)
[1] 0.20322 0.52376
R> moment(fpt2d , center= TRUE , order = 4)
[1] 0.0000000973370 0.0000000024752
R> moment(fpt2d , center= FALSE , order = 4)
[1] 0.00033934 0.06421941

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.13380 0.50332
Variance         0.00017 0.00003
Median           0.13350 0.50333
Mode             0.13444 0.50385
First quartile   0.12488 0.49981
Third quartile   0.14197 0.50674
Minimum          0.07865 0.48107
Maximum          0.20322 0.52376
Skewness         0.26699 0.04723
Kurtosis         3.30487 3.21570
Coef-variation   0.09791 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

Display the exact first-passage-time \(\tau_{S(t)}\), see Figure 5:

R> plot(ts.union(mod2d$X[,1],mod2d$Y[,1]),col=1:2,lty=3,plot.type="single",type="l",ylab= "",xlab="time",axes=F)
R> curve(sin(2*pi*x),add=TRUE,col=3)
R> points(fpt2d$fpt$x[1],sin(2*pi*fpt2d$fpt$x[1]),pch=19,col=4,cex=0.5)
R> lines(c(fpt2d$fpt$x[1],fpt2d$fpt$x[1]),c(sin(2*pi*fpt2d$fpt$x[1]),-10),lty=2,col=4)
R> axis(1, fpt2d$fpt$x[1], bquote(tau[X[S(t)]]==.(fpt2d$fpt$x[1])),col=4,col.ticks=4)
R> points(fpt2d$fpt$y[1],sin(2*pi*fpt2d$fpt$y[1]),pch=19,col=5,cex=0.5)
R> lines(c(fpt2d$fpt$y[1],fpt2d$fpt$y[1]),c(sin(2*pi*fpt2d$fpt$y[1]),-10),lty=2,col=5)
R> axis(1, fpt2d$fpt$y[1], bquote(tau[Y[S(t)]]==.(fpt2d$fpt$y[1])),col=5,col.ticks=5)
R> legend('topright',col=1:5,lty=c(1,1,1,NA,NA),pch=c(NA,NA,NA,19,19),legend=c(expression(X[t]),expression(Y[t]),expression(S(t)),expression(tau[X[S(t)]]),expression(tau[Y[S(t)]])),cex=0.8,inset = .01)
R> box()

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

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

The following \(3\)-dimensional SDE’s with a vector of drift and a diagonal matrix of diffusion coefficients:

\[\begin{equation}\label{eq17} \begin{cases} dX_t = f_{x}(t,X_{t},Y_{t},Z_{t}) dt + g_{x}(t,X_{t},Y_{t},Z_{t}) dW_{1,t}\\ dY_t = f_{y}(t,X_{t},Y_{t},Z_{t}) dt + g_{y}(t,X_{t},Y_{t},Z_{t}) dW_{2,t}\\ dZ_t = f_{z}(t,X_{t},Y_{t},Z_{t}) dt + g_{z}(t,X_{t},Y_{t},Z_{t}) dW_{3,t} \end{cases} \end{equation}\]

\(W_{1,t}\), \(W_{2,t}\) and \(W_{3,t}\) is a 3 independent standard Wiener process. First passage time (3D) \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})},\tau_{(S(t),Z_{t})})\) is defined as:

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

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 dW_{1,t}\\ dY_t = 4 (1-Y_{t}) X_{t} dt + 0.2 dW_{2,t}\\ dZ_t = 4 (1-Z_{t}) Y_{t} dt + 0.2 dW_{3,t} \end{cases} \end{equation}\]

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> mod3d <- snssde3d(drift=fx,diffusion=gx,x0=c(x=2,y=-2,z=0),M=10000)

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> fpt3d
Itô Sde 3D:
    | dX(t) = 4 * (-1 - X(t)) * Y(t) * dt + 0.2 * dW1(t)
    | dY(t) = 4 * (1 - Y(t)) * X(t) * dt + 0.2 * dW2(t)
    | dZ(t) = 4 * (1 - Z(t)) * Y(t) * dt + 0.2 * dW3(t)
    | t in [0,1].
Boundary:
    | S(t) = -1.5 + 3 * t
F.P.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 }
    | Crossing realized 10000 among 10000.
R> head(fpt3d$fpt, n = 10)
         x        y       z
1  0.52621 0.023815 0.81590
2  0.54118 0.026534 0.77393
3  0.51448 0.022113 0.81521
4  0.54429 0.022144 0.80633
5  0.52074 0.023054 0.74946
6  0.53391 0.023566 0.81329
7  0.55821 0.021651 0.75795
8  0.52328 0.024877 0.74534
9  0.53361 0.022729 0.80295
10 0.53325 0.022675 0.75404

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)
[1] 0.53174 0.02324 0.78391
R> moment(fpt3d , center = TRUE , order = 2) ## variance
[1] 0.0001848793 0.0000016501 0.0009496877
R> Median(fpt3d)
[1] 0.531274 0.023217 0.784788
R> Mode(fpt3d)
[1] 0.529649 0.023378 0.790696
R> quantile(fpt3d)
$x
     0%     25%     50%     75%    100% 
0.48742 0.52249 0.53127 0.54078 0.59382 

$y
      0%      25%      50%      75%     100% 
0.018707 0.022367 0.023217 0.024085 0.028541 

$z
     0%     25%     50%     75%    100% 
0.66037 0.76422 0.78479 0.80516 0.90264 
R> kurtosis(fpt3d)
[1] 3.0297 3.0221 3.1072
R> skewness(fpt3d)
[1]  0.15353  0.12206 -0.17140
R> cv(fpt3d)
[1] 0.025572 0.055277 0.039314
R> min(fpt3d)
[1] 0.487421 0.018707 0.660366
R> max(fpt3d)
[1] 0.593821 0.028541 0.902641
R> moment(fpt3d , center= TRUE , order = 4)
[1] 0.0000001035772487 0.0000000000082302 0.0000028029868596
R> moment(fpt3d , center= FALSE , order = 4)
[1] 0.08026053818 0.00000029706 0.38111380516

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.53174 0.02324  0.78391
Variance         0.00018 0.00000  0.00095
Median           0.53127 0.02322  0.78479
Mode             0.52965 0.02338  0.79070
First quartile   0.52249 0.02237  0.76422
Third quartile   0.54078 0.02408  0.80516
Minimum          0.48742 0.01871  0.66037
Maximum          0.59382 0.02854  0.90264
Skewness         0.15353 0.12206 -0.17140
Kurtosis         3.02971 3.02211  3.10722
Coef-variation   0.02557 0.05528  0.03931
3th-order moment 0.00000 0.00000 -0.00001
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

Display the exact first-passage-time \(\tau_{S(t)}\), see Figure 9:

R> plot(ts.union(mod3d$X[,1],mod3d$Y[,1],mod3d$Z[,1]),col=1:3,lty=3,plot.type="single",type="l",ylab="",xlab="time",axes=F)
R> curve(-1.5+3*x,add=TRUE,col=4)
R> points(fpt3d$fpt$x[1],-1.5+3*fpt3d$fpt$x[1],pch=19,col=5,cex=0.5)
R> lines(c(fpt3d$fpt$x[1],fpt3d$fpt$x[1]),c(-1.5+3*fpt3d$fpt$x[1],-10),lty=2,col=5)
R> axis(1, fpt3d$fpt$x[1], bquote(tau[X[S(t)]]==.(fpt3d$fpt$x[1])),col=5,col.ticks=5)
R> points(fpt3d$fpt$y[1],-1.5+3*fpt3d$fpt$y[1],pch=19,col=6,cex=0.5)
R> lines(c(fpt3d$fpt$y[1],fpt3d$fpt$y[1]),c(-1.5+3*fpt3d$fpt$y[1],-10),lty=2,col=6)
R> axis(1, fpt3d$fpt$y[1], bquote(tau[Y[S(t)]]==.(fpt3d$fpt$y[1])),col=6,col.ticks=6)
R> points(fpt3d$fpt$z[1],-1.5+3*fpt3d$fpt$z[1],pch=19,col=7,cex=0.5)
R> lines(c(fpt3d$fpt$z[1],fpt3d$fpt$z[1]),c(-1.5+3*fpt3d$fpt$z[1],-10),lty=2,col=7)
R> axis(1, fpt3d$fpt$z[1], bquote(tau[Z[S(t)]]==.(fpt3d$fpt$z[1])),col=7,col.ticks=7)
R> legend('topright',col=1:7,lty=c(1,1,1,1,NA,NA,NA),pch=c(NA,NA,NA,NA,19,19,19),legend=c(expression(X[t]),expression(Y[t]),expression(Z[t]),expression(S(t)),expression(tau[X[S(t)]]),expression(tau[Y[S(t)]]),expression(tau[Z[S(t)]])),cex=0.8,inset = .01)
R> box()

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, see e.g. Figure 10.

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

Marginal density for the F.P.T of X(t)
    | T(S,X) = inf{t >= 0 : X(t) <= -1.5 + 3 * t}

Data: out[, "x"] (10000 obs.);  Bandwidth 'bw' = 0.0019396

       x                f(x)        
 Min.   :0.48160   Min.   : 0.0002  
 1st Qu.:0.51111   1st Qu.: 0.1347  
 Median :0.54062   Median : 3.0044  
 Mean   :0.54062   Mean   : 8.4636  
 3rd Qu.:0.57013   3rd Qu.:16.1477  
 Max.   :0.59964   Max.   :29.7329  

Marginal density for the F.P.T of Y(t)
    | T(S,Y) = inf{t >= 0 : Y(t) >= -1.5 + 3 * t}

Data: out[, "y"] (10000 obs.);  Bandwidth 'bw' = 0.00018285

       y                 f(y)        
 Min.   :0.018158   Min.   :  0.003  
 1st Qu.:0.020891   1st Qu.:  1.928  
 Median :0.023624   Median : 34.167  
 Mean   :0.023624   Mean   : 91.390  
 3rd Qu.:0.026357   3rd Qu.:174.679  
 Max.   :0.029090   Max.   :309.870  

Marginal density for the F.P.T of Z(t)
    | T(S,Z) = inf{t >= 0 : Z(t) <= -1.5 + 3 * t}

Data: out[, "z"] (10000 obs.);  Bandwidth 'bw' = 0.0043582

       z                 f(z)        
 Min.   :0.018158   Min.   :  0.003  
 1st Qu.:0.020891   1st Qu.:  1.928  
 Median :0.023624   Median : 34.167  
 Mean   :0.023624   Mean   : 91.390  
 3rd Qu.:0.026357   3rd Qu.:174.679  
 Max.   :0.029090   Max.   :309.870  
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. 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 (2017). Performing Parallel Monte Carlo and Moment Equations Methods for Itô and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc. Preprint submitted to Journal of Statistical Software.

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

  6. 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.

  7. 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.

  8. 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)