--- title: "Monte-Carlo Simulation and Kernel Density Estimation of First passage time" author: - A.C. Guidoum^[Department of Mathematics and Computer Science, Faculty of Sciences and Technology, University of Tamanghasset, Algeria, E-mail (acguidoum@univ-tam.dz)] and K. Boukhetala^[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)] date: "`r Sys.Date()`" output: knitr:::html_vignette: toc: yes vignette: > %\VignetteIndexEntry{Monte-Carlo Simulation and Kernel Density Estimation of First passage time} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, echo = F, message = F, results = 'hide',screenshot.force=FALSE} library(Sim.DiffProc) library(knitr) knitr::opts_chunk$set(comment="", prompt=TRUE, fig.show='hold',warning=FALSE, message=FALSE) options(prompt="R> ",scipen=16,digits=5,warning=FALSE, message=FALSE, width = 70) ``` # 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: - `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} set.seed(1234, kind = "L'Ecuyer-CMRG") f <- expression( (1-0.5*x) ) g <- expression( 1 ) 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)) ) fpt1d <- fptsde1d(mod1d, boundary = St) fpt1d head(fpt1d$fpt, n = 5) ``` The following statistical measures (`S3 method`) for class `fptsde1d()` can be approximated for the first-passage-time $\tau_{S(t)}$: ```{r,eval=FALSE, include=TRUE} mean(fpt1d) moment(fpt1d , center = TRUE , order = 2) ## variance Median(fpt1d) Mode(fpt1d) quantile(fpt1d) kurtosis(fpt1d) skewness(fpt1d) cv(fpt1d) min(fpt1d) max(fpt1d) moment(fpt1d , center= TRUE , order = 4) moment(fpt1d , center= FALSE , order = 4) ``` The kernel density approximation of 'fpt1d', using `dfptsde1d()` function (`hist=TRUE` based on `truehist()` function in [MASS](https://cran.r-project.org/package=MASS) package) ```{r 2,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} plot(dfptsde1d(fpt1d),hist=TRUE,nbins="FD") ## histogramm plot(dfptsde1d(fpt1d)) ## kernel density ``` Since [fptdApprox](https://cran.r-project.org/package=fptdApprox) and [DiffusionRgqd](https://cran.r-project.org/package=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](https://cran.r-project.org/package=fptdApprox) to approximate the first passage time density: ```{r,eval=FALSE, include=TRUE} require(fptdApprox) x <- character(4) x[1] <- "m * x" x[2] <- "(sigma^2) * x^2" x[3] <- "dnorm((log(x) - (log(y) + (m - sigma^2/2) * (t- s)))/(sigma * sqrt(t - s)),0,1)/(sigma * sqrt(t - s) * x)" x[4] <- "plnorm(x,log(y) + (m - sigma^2/2) * (t - s),sigma * sqrt(t - s))" Lognormal <- diffproc(x) 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](https://cran.r-project.org/package=Sim.DiffProc) package: ```{r} ## Set the model X(t) f <- expression( 0.48*x ) g <- expression( 0.07*x ) mod1 <- snssde1d(drift=f,diffusion=g,x0=1,T=10,M=1000) ## Set the boundary S(t) St <- expression( 7 + 3.2 * t + 1.4 * t * sin(1.75 * t) ) ## Generate the fpt fpt1 <- fptsde1d(mod1, boundary = St) head(fpt1$fpt, n = 5) summary(fpt1) ``` By plotting the approximations: ```{r 3,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} 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) plot(dfptsde1d(fpt1,bw="bcv"),add=TRUE) legend('topright', lty = c(1, NA), col = c(1,'#BBCCEE'),pch=c(NA,15),legend = c('Approx.fpt.density()', 'fptsde1d()'), lwd = 2, bty = 'n') ``` ```{r 33, echo=FALSE, fig.cap=' `fptsde1d()` vs `Approx.fpt.density()` ', fig.env='figure*',fig.width=7,fig.height=7} knitr::include_graphics("Figures/fig01.png") ``` ### `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](https://cran.r-project.org/package=DiffusionRgqd) to approximate the first passage time density: ```{r,eval=FALSE, include=TRUE} require(DiffusionRgqd) G1 <- function(t) { theta[1] * (10+0.2 * sin(2 * pi * t) + 0.3 * prod(sqrt(t), 1+cos(3 * pi * t))) } G2 <- function(t){-theta[1]} Q2 <- function(t){0.1} res2 = GQD.TIpassage(8, 12, 1, 4, 1 / 100, theta = c(0.5)) ``` Using `fptsde1d()` and `dfptsde1d()` functions in the [Sim.DiffProc](https://cran.r-project.org/package=Sim.DiffProc) package: ```{r} ## Set the model X(t) theta1=0.5 f <- expression( theta1*x*(10+0.2*sin(2*pi*t)+0.3*sqrt(t)*(1+cos(3*pi*t))-x) ) g <- expression( sqrt(0.1)*x ) mod2 <- snssde1d(drift=f,diffusion=g,x0=8,t0=1,T=4,M=1000) ## Set the boundary S(t) St <- expression( 12 ) ## Generate the fpt fpt2 <- fptsde1d(mod2, boundary = St) head(fpt2$fpt, n = 5) summary(fpt2) ``` By plotting the approximations (`hist=TRUE` based on `truehist()` function in [MASS](https://cran.r-project.org/package=MASS) package): ```{r 4,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} plot(dfptsde1d(fpt2),hist=TRUE,nbins = "Scott",main = 'Approximation First-Passage-Time Density', ylab = 'Density', xlab = expression(tau[S(t)]), cex.main = 0.95) lines(res2$density ~ res2$time, type = 'l',lwd=2) legend('topright', lty = c(1, NA), col = c(1,'#FF00004B'),pch=c(NA,15),legend = c('GQD.TIpassage()', 'fptsde1d()'), lwd = 2, bty = 'n') ``` ```{r 44, echo=FALSE, fig.cap='`fptsde1d()` vs `GQD.TIpassage()` ', fig.env='figure*',fig.width=7,fig.height=7} knitr::include_graphics("Figures/fig02.png") ``` ## FPT for 2-Dim SDE's Assume that we want to describe the following Stratonovich SDE's (2D): \begin{equation}\label{eq016} \begin{cases} dX_t = 5 (-1-Y_{t}) X_{t} dt + 0.5 Y_{t} \circ dW_{1,t}\\ dY_t = 5 (-1-X_{t}) Y_{t} dt + 0.5 X_{t} \circ dW_{2,t} \end{cases} \end{equation} and \[ S(t)=\sin(2\pi t) \] Set the system $(X_t , Y_t)$: ```{r} set.seed(1234, kind = "L'Ecuyer-CMRG") fx <- expression(5*(-1-y)*x , 5*(-1-x)*y) gx <- expression(0.5*y,0.5*x) 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)) fpt2d <- fptsde2d(mod2d, boundary = St) head(fpt2d$fpt, n = 5) ``` 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,eval=FALSE, include=TRUE} mean(fpt2d) moment(fpt2d , center = TRUE , order = 2) ## variance Median(fpt2d) Mode(fpt2d) quantile(fpt2d) kurtosis(fpt2d) skewness(fpt2d) cv(fpt2d) min(fpt2d) max(fpt2d) moment(fpt2d , center= TRUE , order = 4) moment(fpt2d , center= FALSE , order = 4) ``` The result summaries of the couple \((\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})})\): ```{r} summary(fpt2d) ``` The marginal density of \((\tau_{(S(t),X_{t})}\) and \(\tau_{(S(t),Y_{t})})\) are reported using `dfptsde2d()` function. ```{r 6,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} denM <- dfptsde2d(fpt2d, pdf = 'M') 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 7,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} denJ <- dfptsde2d(fpt2d, pdf = 'J',n=100) plot(denJ,display="contour",main="Bivariate Density of F.P.T",xlab=expression(tau[x]),ylab=expression(tau[y])) 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 8, echo=TRUE, fig.cap=' ', fig.env='figure*', message=FALSE, warning=FALSE,eval=FALSE, include=TRUE} plot(denJ,display="persp",main="Bivariate Density of F.P.T",xlab=expression(tau[x]),ylab=expression(tau[y])) ``` [Return to fptsde2d()](#FPT for 2-Dim SDE) ## FPT for 3-Dim SDE's Assume that we want to describe the following SDE's (3D): \begin{equation}\label{eq0166} \begin{cases} dX_t = 4 (-1-X_{t}) Y_{t} dt + 0.2 dB_{1,t}\\ dY_t = 4 (1-Y_{t}) X_{t} dt + 0.2 dB_{2,t}\\ dZ_t = 4 (1-Z_{t}) Y_{t} dt + 0.2 dB_{3,t} \end{cases} \end{equation} with $(B_{1,t},B_{2,t},B_{3,t})$ are three correlated standard Wiener process: $$ \Sigma= \begin{pmatrix} 1 & 0.3 &-0.5\\ 0.3 & 1 & 0.2 \\ -0.5 &0.2&1 \end{pmatrix} $$ and $$ S(t)=-1.5+3t $$ Set the system $(X_t , Y_t , Z_t)$: ```{r} set.seed(1234, kind = "L'Ecuyer-CMRG") fx <- expression(4*(-1-x)*y , 4*(1-y)*x , 4*(1-z)*y) gx <- rep(expression(0.2),3) Sigma <-matrix(c(1,0.3,-0.5,0.3,1,0.2,-0.5,0.2,1),nrow=3,ncol=3) mod3d <- snssde3d(drift=fx,diffusion=gx,x0=c(x=2,y=-2,z=0),M=1000,corr=Sigma) ``` Generate the triplet $(\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})},\tau_{(S(t),Z_{t})})$, with `fptsde3d()` function:: ```{r} St <- expression(-1.5+3*t) fpt3d <- fptsde3d(mod3d, boundary = St) head(fpt3d$fpt, n = 5) ``` 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,eval=FALSE, include=TRUE} mean(fpt3d) moment(fpt3d , center = TRUE , order = 2) ## variance Median(fpt3d) Mode(fpt3d) quantile(fpt3d) kurtosis(fpt3d) skewness(fpt3d) cv(fpt3d) min(fpt3d) max(fpt3d) moment(fpt3d , center= TRUE , order = 4) moment(fpt3d , center= FALSE , order = 4) ``` The result summaries of the triplet $(\tau_{(S(t),X_{t})},\tau_{(S(t),Y_{t})},\tau_{(S(t),Z_{t})})$: ```{r} summary(fpt3d) ``` 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. ```{r 10,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} denM <- dfptsde3d(fpt3d, pdf = "M") 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](https://cran.r-project.org/package=sm) or [ks](https://cran.r-project.org/package=ks).) ```{r 111,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} denJ <- dfptsde3d(fpt3d,pdf="J") plot(denJ,display="rgl") ``` [Return to fptsde3d()](#fptsde3d) # Further reading 1. [`snssdekd()` & `dsdekd()` & `rsdekd()`- Monte-Carlo Simulation and Analysis of Stochastic Differential Equations](snssde.html). 2. [`bridgesdekd()` & `dsdekd()` & `rsdekd()` - Constructs and Analysis of Bridges Stochastic Differential Equations](bridgesde.html). 3. [`fptsdekd()` & `dfptsdekd()` - Monte-Carlo Simulation and Kernel Density Estimation of First passage time](fptsde.html). 4. [`MCM.sde()` & `MEM.sde()` - Parallel Monte-Carlo and Moment Equations for SDEs](mcmsde.html). 5. [`TEX.sde()` - Converting Sim.DiffProc Objects to LaTeX](sdetotex.html). 6. [`fitsde()` - Parametric Estimation of 1-D Stochastic Differential Equation](fitsde.html). # References 1. Boukhetala K (1996). Modelling and Simulation of a Dispersion Pollutant with Attractive Centre, volume 3, pp. 245-252. Computer Methods and Water Resources, Computational Mechanics Publications, Boston, USA. 2. Boukhetala K (1998). Estimation of the first passage time distribution for a simulated diffusion process. Maghreb Mathematical Review, 7, pp. 1-25. 3. Boukhetala K (1998). Kernel density of the exit time in a simulated diffusion. The Annals of The Engineer Maghrebian, 12, pp. 587-589. 4. Guidoum AC, Boukhetala K (2024). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.9, URL https://cran.r-project.org/package=Sim.DiffProc. 5. Pienaar EAD, Varughese MM (2016). DiffusionRgqd: An R Package for Performing Inference and Analysis on Time-Inhomogeneous Quadratic Diffusion Processes. R package version 0.1.3, URL https://CRAN.R-project.org/package=DiffusionRgqd. 6. Roman, R.P., Serrano, J. J., Torres, F. (2008). First-passage-time location function: Application to determine first-passage-time densities in diffusion processes. Computational Statistics and Data Analysis. 52, 4132-4146. 7. Roman, R.P., Serrano, J. J., Torres, F. (2012). An R package for an efficient approximation of first-passage-time densities for diffusion processes based on the FPTL function. Applied Mathematics and Computation, 218, 8408-8428.