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

A.C. Guidoum1 and K. Boukhetala2

2018-11-28

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 975 among 1000.
R> head(fpt1d$fpt, n = 10)
 [1] 0.096261 0.464695 0.509866 0.027795 0.044845 0.036563 0.100767
 [8] 0.020404 0.079737 0.043660

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.19036
R> moment(fpt1d , center = TRUE , order = 2) ## variance
[1] 0.042175
R> Median(fpt1d)
[1] 0.11191
R> Mode(fpt1d)
[1] 0.058009
R> quantile(fpt1d)
       0%       25%       50%       75%      100% 
0.0094354 0.0563074 0.1119142 0.2346988 0.9950344 
R> kurtosis(fpt1d)
[1] 6.1617
R> skewness(fpt1d)
[1] 1.8841
R> cv(fpt1d)
[1] 1.0793
R> min(fpt1d)
[1] 0.0094354
R> max(fpt1d)
[1] 0.99503
R> moment(fpt1d , center= TRUE , order = 4)
[1] 0.010982
R> moment(fpt1d , center= FALSE , order = 4)
[1] 0.03391

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.19036
Variance         0.04222
Median           0.11191
Mode             0.05801
First quartile   0.05631
Third quartile   0.23470
Minimum          0.00944
Maximum          0.99503
Skewness         1.88408
Kurtosis         6.16170
Coef-variation   1.07935
3th-order moment 0.01634
4th-order moment 0.01098
5th-order moment 0.00698
6th-order moment 0.00475

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=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> 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 1000 among 1000.
R> head(fpt1$fpt, n = 10)
 [1] 8.7914 5.8383 6.1322 6.6407 5.9285 8.5290 6.1353 6.1978 8.2015
[10] 5.7198
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.54698
Variance          0.91449
Median            6.14013
Mode              6.07867
First quartile    5.96915
Third quartile    6.43835
Minimum           5.62479
Maximum           8.98987
Skewness          1.43627
Kurtosis          3.34637
Coef-variation    0.14607
3th-order moment  1.25604
4th-order moment  2.79853
5th-order moment  5.34044
6th-order moment 10.84165

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> 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 925 among 1000.
R> head(fpt2$fpt, n = 10)
 [1] 2.3961 1.3939 1.6272 2.6868 2.4874 1.4654 2.7769 1.3799 1.9608
[10] 2.2019
R> summary(fpt2)

Monte-Carlo Statistics of F.P.T:
|T(S(t),X(t)) = inf{t >=  1 : X(t) >=  12 }
                        
Mean             2.13599
Variance         0.47712
Median           2.03809
Mode             1.44323
First quartile   1.52003
Third quartile   2.57763
Minimum          1.10178
Maximum          3.98158
Skewness         0.68436
Kurtosis         2.59034
Coef-variation   0.32338
3th-order moment 0.22554
4th-order moment 0.58969
5th-order moment 0.62583
6th-order moment 1.12710

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=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> 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 1000 among 1000.
R> head(fpt2d$fpt, n = 10)
         x       y
1  0.14772 0.50296
2  0.12583 0.49905
3  0.15479 0.49246
4  0.14647 0.50153
5  0.15205 0.50438
6  0.15565 0.49350
7  0.12404 0.50258
8  0.14260 0.50512
9  0.12047 0.50759
10 0.14081 0.49838

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.13391 0.50334
R> moment(fpt2d , center = TRUE , order = 2) ## variance
[1] 0.000171725 0.000028003
R> Median(fpt2d)
[1] 0.13354 0.50320
R> Mode(fpt2d)
[1] 0.13396 0.50300
R> quantile(fpt2d)
$x
      0%      25%      50%      75%     100% 
0.089985 0.125564 0.133543 0.142065 0.184686 

$y
     0%     25%     50%     75%    100% 
0.48606 0.49984 0.50320 0.50673 0.52257 
R> kurtosis(fpt2d)
[1] 3.3818 3.2046
R> skewness(fpt2d)
[1] 0.18393 0.17609
R> cv(fpt2d)
[1] 0.097906 0.010519
R> min(fpt2d)
[1] 0.089985 0.486055
R> max(fpt2d)
[1] 0.18469 0.52257
R> moment(fpt2d , center= TRUE , order = 4)
[1] 0.000000099928 0.000000002518
R> moment(fpt2d , center= FALSE , order = 4)
[1] 0.00034038 0.06422726

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.13391 0.50334
Variance         0.00017 0.00003
Median           0.13354 0.50320
Mode             0.13396 0.50300
First quartile   0.12556 0.49984
Third quartile   0.14206 0.50673
Minimum          0.08999 0.48606
Maximum          0.18469 0.52257
Skewness         0.18393 0.17609
Kurtosis         3.38183 3.20459
Coef-variation   0.09791 0.01052
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=1000)

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 1000 among 1000.
R> head(fpt3d$fpt, n = 10)
         x        y       z
1  0.51515 0.022702 0.77990
2  0.53911 0.024360 0.85046
3  0.52964 0.021566 0.74697
4  0.51553 0.023404 0.77999
5  0.52309 0.020915 0.76343
6  0.49688 0.022795 0.80233
7  0.53436 0.024997 0.79127
8  0.51361 0.024625 0.82525
9  0.50504 0.022368 0.75422
10 0.53545 0.025108 0.73533

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.531420 0.023284 0.783810
R> moment(fpt3d , center = TRUE , order = 2) ## variance
[1] 0.0001904684 0.0000016756 0.0009550880
R> Median(fpt3d)
[1] 0.53100 0.02324 0.78545
R> Mode(fpt3d)
[1] 0.531049 0.023309 0.787171
R> quantile(fpt3d)
$x
     0%     25%     50%     75%    100% 
0.49521 0.52228 0.53100 0.53988 0.57769 

$y
      0%      25%      50%      75%     100% 
0.019306 0.022376 0.023240 0.024121 0.028735 

$z
     0%     25%     50%     75%    100% 
0.68032 0.76400 0.78545 0.80508 0.86922 
R> kurtosis(fpt3d)
[1] 3.0739 2.9687 3.2016
R> skewness(fpt3d)
[1]  0.22514  0.18703 -0.27444
R> cv(fpt3d)
[1] 0.025983 0.055621 0.039448
R> min(fpt3d)
[1] 0.495207 0.019306 0.680319
R> max(fpt3d)
[1] 0.577687 0.028735 0.869215
R> moment(fpt3d , center= TRUE , order = 4)
[1] 0.0000001117385477 0.0000000000083518 0.0000029263465403
R> moment(fpt3d , center= FALSE , order = 4)
[1] 0.08007775773 0.00000029944 0.38093388098

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.53142 0.02328  0.78381
Variance         0.00019 0.00000  0.00096
Median           0.53100 0.02324  0.78545
Mode             0.53105 0.02331  0.78717
First quartile   0.52228 0.02238  0.76400
Third quartile   0.53988 0.02412  0.80508
Minimum          0.49521 0.01931  0.68032
Maximum          0.57769 0.02874  0.86922
Skewness         0.22514 0.18703 -0.27444
Kurtosis         3.07389 2.96872  3.20162
Coef-variation   0.02598 0.05562  0.03945
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"] (1000 obs.);   Bandwidth 'bw' = 0.0029707

       x                f(x)        
 Min.   :0.48629   Min.   : 0.0016  
 1st Qu.:0.51137   1st Qu.: 0.9135  
 Median :0.53645   Median : 6.2259  
 Mean   :0.53645   Mean   : 9.9600  
 3rd Qu.:0.56152   3rd Qu.:18.3844  
 Max.   :0.58660   Max.   :29.8335  

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"] (1000 obs.);   Bandwidth 'bw' = 0.00029278

       y                 f(y)        
 Min.   :0.018428   Min.   :  0.015  
 1st Qu.:0.021224   1st Qu.:  1.322  
 Median :0.024021   Median : 37.723  
 Mean   :0.024021   Mean   : 89.312  
 3rd Qu.:0.026817   3rd Qu.:167.293  
 Max.   :0.029614   Max.   :291.043  

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"] (1000 obs.);   Bandwidth 'bw' = 0.0069309

       z                 f(z)        
 Min.   :0.018428   Min.   :  0.015  
 1st Qu.:0.021224   1st Qu.:  1.322  
 Median :0.024021   Median : 37.723  
 Mean   :0.024021   Mean   : 89.312  
 3rd Qu.:0.026817   3rd Qu.:167.293  
 Max.   :0.029614   Max.   :291.043  
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 (2018). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.3, 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)