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> fpt1dItô 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.043660The following statistical measures (S3 method) for class fptsde1d() can be approximated for the first-passage-time \(\tau_{S(t)}\):
R> mean(fpt1d)[1] 0.19036R> moment(fpt1d , center = TRUE , order = 2) ## variance[1] 0.042175R> Median(fpt1d)[1] 0.11191R> Mode(fpt1d)[1] 0.058009R> quantile(fpt1d) 0% 25% 50% 75% 100%
0.0094354 0.0563074 0.1119142 0.2346988 0.9950344 R> kurtosis(fpt1d)[1] 6.1617R> skewness(fpt1d)[1] 1.8841R> cv(fpt1d)[1] 1.0793R> min(fpt1d)[1] 0.0094354R> max(fpt1d)[1] 0.99503R> moment(fpt1d , center= TRUE , order = 4)[1] 0.010982R> moment(fpt1d , center= FALSE , order = 4)[1] 0.03391The 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.00475Display 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> fpt1Itô 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.7198R> 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.84165By 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 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> fpt2Itô 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.2019R> 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.12710By 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()
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> fpt2dStratonovich 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.49838The 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.50334R> moment(fpt2d , center = TRUE , order = 2) ## variance[1] 0.000171725 0.000028003R> Median(fpt2d)[1] 0.13354 0.50320R> Mode(fpt2d)[1] 0.13396 0.50300R> 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.2046R> skewness(fpt2d)[1] 0.18393 0.17609R> cv(fpt2d)[1] 0.097906 0.010519R> min(fpt2d)[1] 0.089985 0.486055R> max(fpt2d)[1] 0.18469 0.52257R> moment(fpt2d , center= TRUE , order = 4)[1] 0.000000099928 0.000000002518R> moment(fpt2d , center= FALSE , order = 4)[1] 0.00034038 0.06422726The 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.00000Display 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]))
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> fpt3dItô 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.73533The 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.783810R> moment(fpt3d , center = TRUE , order = 2) ## variance[1] 0.0001904684 0.0000016756 0.0009550880R> Median(fpt3d)[1] 0.53100 0.02324 0.78545R> Mode(fpt3d)[1] 0.531049 0.023309 0.787171R> 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.2016R> skewness(fpt3d)[1] 0.22514 0.18703 -0.27444R> cv(fpt3d)[1] 0.025983 0.055621 0.039448R> min(fpt3d)[1] 0.495207 0.019306 0.680319R> max(fpt3d)[1] 0.577687 0.028735 0.869215R> moment(fpt3d , center= TRUE , order = 4)[1] 0.0000001117385477 0.0000000000083518 0.0000029263465403R> moment(fpt3d , center= FALSE , order = 4)[1] 0.08007775773 0.00000029944 0.38093388098The 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.00000Display 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")