# 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:

$$$\label{eds01} dX_t = \mu(t,X_t) dt + \sigma(t,X_t) dW_t,\quad X_{t_{0}}=x_{0}$$$

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:

• object an object inheriting from class snssde1d, snssde2d and snssde3d.
• boundary an expression of a constant or time-dependent boundary $$S(t)$$.

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

• The expected value $$\text{E}(\tau_{S(t)})$$, using the command mean.
• The variance $$\text{Var}(\tau_{S(t)})$$, using the command moment with order=2 and center=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 command min and max.
• The skewness and the kurtosis of $$\tau_{S(t)}$$, using the command skewness and kurtosis.
• The coefficient of variation (relative variability) of $$\tau_{S(t)}$$, using the command cv.
• The central moments up to order $$p$$ of $$\tau_{S(t)}$$, using the command moment.
• The result summaries of the results of Monte-Carlo simulation, using the command summary.

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

• object an object inheriting from class fptsdekd() (where k=1,2,3).
• pdf probability density function Joint or Marginal.

# 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 966 among 1000.
R> head(fpt1d$fpt, n = 10)  [1] 0.1339457 0.0403355 0.0185658 0.1145573 0.0091086 0.1376444 [7] 0.6836460 0.2273584 0.0287167 0.0967543 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.20261 R> moment(fpt1d , center = TRUE , order = 2) ## variance [1] 0.046399 R> Median(fpt1d) [1] 0.11883 R> Mode(fpt1d) [1] 0.056178 R> quantile(fpt1d)  0% 25% 50% 75% 100% 0.0069061 0.0539587 0.1188284 0.2503842 0.9993877  R> kurtosis(fpt1d) [1] 5.5196 R> skewness(fpt1d) [1] 1.7513 R> cv(fpt1d) [1] 1.0637 R> min(fpt1d) [1] 0.0069061 R> max(fpt1d) [1] 0.99939 R> moment(fpt1d , center= TRUE , order = 4) [1] 0.011908 R> moment(fpt1d , center= FALSE , order = 4) [1] 0.039228 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.20261 Variance 0.04645 Median 0.11883 Mode 0.05618 First quartile 0.05396 Third quartile 0.25038 Minimum 0.00691 Maximum 0.99939 Skewness 1.75134 Kurtosis 5.51959 Coef-variation 1.06372 3th-order moment 0.01753 4th-order moment 0.01191 5th-order moment 0.00749 6th-order moment 0.00508 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> 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)]==.(fpt1dfpt[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(fpt1fpt, n = 10)
 [1] 5.9592 6.2240 6.0662 6.1222 6.1849 6.0208 5.9970 6.0740 6.1962
[10] 6.3749
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.51837
Variance          0.90712
Median            6.11896
Mode              6.03207
First quartile    5.95884
Third quartile    6.40436
Minimum           5.40261
Maximum           8.95967
Skewness          1.45679
Kurtosis          3.40877
Coef-variation    0.14611
3th-order moment  1.25861
4th-order moment  2.80494
5th-order moment  5.37535
6th-order moment 10.95455

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> 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 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 923 among 1000.
R> head(fpt2$fpt, n = 10)  [1] 2.6795 1.4390 1.8486 1.4483 1.4806 1.3327 2.5660 1.6606 3.2797 [10] 2.1307 R> summary(fpt2)  Monte-Carlo Statistics of F.P.T: |T(S(t),X(t)) = inf{t >= 1 : X(t) >= 12 } Mean 2.18245 Variance 0.51251 Median 2.06456 Mode 1.46359 First quartile 1.54973 Third quartile 2.65918 Minimum 1.15462 Maximum 3.99958 Skewness 0.66462 Kurtosis 2.52209 Coef-variation 0.32802 3th-order moment 0.24385 4th-order moment 0.66246 5th-order moment 0.70278 6th-order moment 1.30438 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') ## FPT for 2-Dim SDE’s The following $$2$$-dimensional SDE’s with a vector of drift and a diagonal matrix of diffusion coefficients: $$$\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}$$$ $$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): $$$\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}$$$ 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.12680 0.50495
2  0.13276 0.50512
3  0.14800 0.50806
4  0.12442 0.50139
5  0.13105 0.50873
6  0.14841 0.50785
7  0.14134 0.49099
8  0.14555 0.51156
9  0.13435 0.51354
10 0.11551 0.49764

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.13377 0.50309
R> moment(fpt2d , center = TRUE , order = 2) ## variance
[1] 0.000185360 0.000028352
R> Median(fpt2d)
[1] 0.13266 0.50319
R> Mode(fpt2d)
[1] 0.13091 0.50198
R> quantile(fpt2d)
$x 0% 25% 50% 75% 100% 0.094199 0.124116 0.132663 0.142240 0.186445$y
0%     25%     50%     75%    100%
0.48659 0.49997 0.50319 0.50659 0.52154 
R> kurtosis(fpt2d)
[1] 3.3985 3.2624
R> skewness(fpt2d)
[1]  0.378860 -0.084235
R> cv(fpt2d)
[1] 0.101825 0.010589
R> min(fpt2d)
[1] 0.094199 0.486586
R> max(fpt2d)
[1] 0.18645 0.52154
R> moment(fpt2d , center= TRUE , order = 4)
[1] 0.0000001170023 0.0000000026277
R> moment(fpt2d , center= FALSE , order = 4)
[1] 0.00034078 0.06410149

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.13377  0.50309
Variance         0.00019  0.00003
Median           0.13266  0.50319
Mode             0.13091  0.50198
First quartile   0.12412  0.49997
Third quartile   0.14224  0.50659
Minimum          0.09420  0.48659
Maximum          0.18645  0.52154
Skewness         0.37886 -0.08423
Kurtosis         3.39853  3.26240
Coef-variation   0.10182  0.01059
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> 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]))

## FPT for 3-Dim SDE’s

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

$$$\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}$$$ $$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): $$$\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}$$$ 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.55604 0.024212 0.74398 2 0.57666 0.023595 0.76097 3 0.52223 0.024165 0.76671 4 0.51569 0.025876 0.81817 5 0.54366 0.020797 0.72875 6 0.52437 0.025312 0.78193 7 0.54578 0.024291 0.85570 8 0.55933 0.023881 0.74460 9 0.53802 0.021959 0.81993 10 0.52354 0.026541 0.79178 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.531421 0.023273 0.783244 R> moment(fpt3d , center = TRUE , order = 2) ## variance [1] 0.0001797741 0.0000016317 0.0010413796 R> Median(fpt3d) [1] 0.531635 0.023249 0.784739 R> Mode(fpt3d) [1] 0.532304 0.023071 0.786603 R> quantile(fpt3d) $x
0%     25%     50%     75%    100%
0.48934 0.52283 0.53163 0.53961 0.57670

$y 0% 25% 50% 75% 100% 0.018989 0.022437 0.023249 0.024127 0.027370$z
0%     25%     50%     75%    100%
0.66295 0.76184 0.78474 0.80541 0.87300 
R> kurtosis(fpt3d)
[1] 3.3214 2.9400 3.0973
R> skewness(fpt3d)
[1]  0.078856  0.055384 -0.234933
R> cv(fpt3d)
[1] 0.025243 0.054913 0.041222
R> min(fpt3d)
[1] 0.489340 0.018989 0.662954
R> max(fpt3d)
[1] 0.57670 0.02737 0.87300
R> moment(fpt3d , center= TRUE , order = 4)
[1] 0.000000107559478 0.000000000007843 0.000003365691761
R> moment(fpt3d , center= FALSE , order = 4)
[1] 0.0800593358 0.0000002987 0.3801585251

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.02327  0.78324
Variance         0.00018 0.00000  0.00104
Median           0.53163 0.02325  0.78474
Mode             0.53230 0.02307  0.78660
First quartile   0.52283 0.02244  0.76184
Third quartile   0.53961 0.02413  0.80541
Minimum          0.48934 0.01899  0.66295
Maximum          0.57670 0.02737  0.87300
Skewness         0.07886 0.05538 -0.23493
Kurtosis         3.32143 2.93997  3.09733
Coef-variation   0.02524 0.05491  0.04122
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.0028298

x                f(x)
Min.   :0.48085   Min.   : 0.002
1st Qu.:0.50693   1st Qu.: 0.564
Median :0.53302   Median : 4.417
Mean   :0.53302   Mean   : 9.575
3rd Qu.:0.55910   3rd Qu.:17.384
Max.   :0.58519   Max.   :34.315

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

y                 f(y)
Min.   :0.018133   Min.   :  0.02
1st Qu.:0.020656   1st Qu.:  3.27
Median :0.023179   Median : 50.70
Mean   :0.023179   Mean   : 98.99
3rd Qu.:0.025702   3rd Qu.:177.25
Max.   :0.028225   Max.   :317.12

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

z                 f(z)
Min.   :0.018133   Min.   :  0.02
1st Qu.:0.020656   1st Qu.:  3.27
Median :0.023179   Median : 50.70
Mean   :0.023179   Mean   : 98.99
3rd Qu.:0.025702   3rd Qu.:177.25
Max.   :0.028225   Max.   :317.12  
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")

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

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