--- title: "Constructs and Analysis of Bridges Stochastic Differential Equations" 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{Constructs and Analysis of Bridges Stochastic Differential Equations} %\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) ``` # Bridges SDE's Consider now a $d$-dimensional stochastic process $X_{t}$ defined on a probability space $(\Omega, \mathfrak{F},\mathbb{P})$. We say that the bridge associated to $X_{t}$ conditioned to the event $\{X_{T}= a\}$ is the process: $$ \{\tilde{X}_{t}, t_{0} \leq t \leq T \}=\{X_{t}, t_{0} \leq t \leq T: X_{T}= a \} $$ where $T$ is a deterministic fixed time and $a \in \mathbb{R}^d$ is fixed too. # The `bridgesdekd()` functions The (S3) generic function `bridgesdekd()` (where `k=1,2,3`) for simulation of 1,2 and 3-dim bridge stochastic differential equations,Ito or Stratonovich type, with different methods. The main arguments consist: - The `drift` and `diffusion` coefficients as R expressions that depend on the state variable `x` (`y` and `z`) and time variable `t`. - `corr` the correlation structure of standard Wiener process ; must be a real symmetric positive-definite square matrix (2D and 3D, default: `corr=NULL`). - The number of simulation steps `N`. - The number of the solution trajectories to be simulated by `M` (default: `M=1`). - Initial value `x0` at initial time `t0`. - Terminal value `y` final time `T` - The integration step size `Dt` (default: `Dt=(T-t0)/N`). - The choice of process types by the argument `type="ito"` for Ito or `type="str"` for Stratonovich (by default `type="ito"`). - The numerical method to be used by `method` (default `method="euler"`). By Monte-Carlo simulations, the following statistical measures (`S3 method`) for class `bridgesdekd()` (where `k=1,2,3`) can be approximated for the process at any time $t \in [t_{0},T]$ (default: `at=(T-t0)/2`): * The expected value $\text{E}(X_{t})$ at time $t$, using the command `mean`. * The variance $\text{Var}(X_{t})$ at time $t$, using the command `moment` with `order=2` and `center=TRUE`. * The median $\text{Med}(X_{t})$ at time $t$, using the command `Median`. * The mode $\text{Mod}(X_{t})$ at time $t$, using the command `Mode`. * The quartile of $X_{t}$ at time $t$, using the command `quantile`. * The maximum and minimum of $X_{t}$ at time $t$, using the command `min` and `max`. * The skewness and the kurtosis of $X_{t}$ at time $t$, using the command `skewness` and `kurtosis`. * The coefficient of variation (relative variability) of $X_{t}$ at time $t$, using the command `cv`. * The central moments up to order $p$ of $X_{t}$ at time $t$, using the command `moment`. * The result summaries of the results of Monte-Carlo simulation at time $t$, using the command `summary`. We can just make use of the `rsdekd()` function (where `k=1,2,3`) to build our random number for class `bridgesdekd()` (where `k=1,2,3`) at any time $t \in [t_{0},T]$. the main arguments consist: - `object` an object inheriting from class `bridgesdekd()` (where `k=1,2,3`). - `at` time between $s=t0$ and $t=T$. The function `dsde()` (where `k=1,2,3`) approximate transition density for class `bridgesdekd()` (where `k=1,2,3`), the main arguments consist: - `object` an object inheriting from class `bridgesdekd()` (where `k=1,2,3`). - `at` time between $s=t0$ and $t=T$. - `pdf` probability density function `Joint` or `Marginal`. The following we explain how to use this functions. # `bridgesde1d()` Assume that we want to describe the following bridge sde in Ito form: \begin{equation}\label{eq0166} dX_t = \frac{1-X_t}{1-t} dt + X_t dW_{t},\quad X_{t_{0}}=3 \quad\text{and}\quad X_{T}=1 \end{equation} We simulate a flow of $1000$ trajectories, with integration step size $\Delta t = 0.001$, and $x_0 = 3$ at time $t_0 = 0$, $y = 1$ at terminal time $T=1$. ```{r} set.seed(1234, kind = "L'Ecuyer-CMRG") f <- expression((1-x)/(1-t)) g <- expression(x) mod <- bridgesde1d(drift=f,diffusion=g,x0=3,y=1,M=1000) mod ``` In Figure 1, we present the flow of trajectories, the mean path (red lines) of solution of $X_{t}|X_{0}=3,X_{T}=1$: ```{r 01,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} plot(mod,ylab=expression(X[t])) lines(time(mod),apply(mod$X,1,mean),col=2,lwd=2) legend("topleft","mean path",inset = .01,col=2,lwd=2,cex=0.8,bty="n") ``` Hence we can just make use of the `rsde1d()` function to build our random number generator for $X_{t}|X_{0}=3,X_{T}=1$ at time $t=0.55$: ```{r} x <- rsde1d(object = mod, at = 0.55) head(x, n = 3) ``` The function `dsde1d()` can be used to show the kernel density estimation for $X_{t}|X_{0}=3,X_{T}=1$ at time $t=0.55$ (`hist=TRUE` based on `truehist()` function in [MASS](https://cran.r-project.org/package=MASS) package): ```{r 04,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} dens <- dsde1d(mod, at = 0.55) plot(dens,hist=TRUE) ## histgramme plot(dens,add=TRUE) ## kernel density ``` ```{r 33, echo=FALSE, fig.cap='Bridge sde 1D. Histgramme and kernel density estimation for $X_{t}|X_{0}=3,X_{T}=1$', fig.env='figure*',fig.width=7,fig.height=7} knitr::include_graphics(c("Figures/fig03.png","Figures/fig1008.png")) ``` [Return to bridgesde1d()](#bridgesde1d) # `bridgesde2d()` Assume that we want to describe the following $2$-dimensional bridge SDE's in Stratonovich form: \begin{equation}\label{eq:09} \begin{cases} dX_t = -(1+Y_{t}) X_{t} dt + 0.2 (1-Y_{t})\circ dB_{1,t},\quad X_{t_{0}}=1 \quad\text{and}\quad X_{T}=1\\ dY_t = -(1+X_{t}) Y_{t} dt + 0.1 (1-X_{t}) \circ dB_{2,t},\quad Y_{t_{0}}=-0.5 \quad\text{and}\quad Y_{T}=0.5 \end{cases} \end{equation} with $(B_{1,t},B_{2,t})$ are two correlated standard Wiener process: $$ \Sigma= \begin{pmatrix} 1 & 0.3\\ 0.3 & 1 \end{pmatrix} $$ We simulate a flow of $1000$ trajectories, with integration step size $\Delta t = 0.01$: ```{r} set.seed(1234, kind = "L'Ecuyer-CMRG") fx <- expression(-(1+y)*x , -(1+x)*y) gx <- expression(0.2*(1-y),0.1*(1-x)) Sigma <-matrix(c(1,0.3,0.3,1),nrow=2,ncol=2) mod2 <- bridgesde2d(drift=fx,diffusion=gx,x0=c(1,-0.5),y=c(1,0.5),Dt=0.01,M=1000,type="str",corr=Sigma) mod2 ``` In Figure 2, we present the flow of trajectories of $X_{t}|X_{0}=1,X_{T}=1$ and $Y_{t}|Y_{0}=-0.5,Y_{T}=0.5$: ```{r 06,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} plot(mod2,col=c('#FF00004B','#0000FF82')) ``` ```{r 333, echo=FALSE, fig.cap=' Bridge sde 2D ', fig.env='figure*',fig.width=5,fig.height=5} knitr::include_graphics("Figures/fig04.png") ``` Hence we can just make use of the `rsde2d()` function to build our random number generator for the couple $X_{t},Y_{t}|X_{0}=1,Y_{0}=-0.5,X_{T}=1,Y_{T}=0.5$ at time $t=5$: ```{r} x2 <- rsde2d(object = mod2, at = 5) head(x2, n = 3) ``` The marginal density of $X_{t}|X_{0}=1,X_{T}=1$ and $Y_{t}|Y_{0}=-0.5,Y_{T}=0.5$ at time $t=5$ are reported using `dsde2d()` function. A `contour` plot of joint density obtained from a realization of the couple $X_{t},Y_{t}|X_{0}=1,Y_{0}=-0.5,X_{T}=1,Y_{T}=0.5$ at time $t=5$, see e.g. Figure 3: ```{r 09,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} ## Marginal denM <- dsde2d(mod2,pdf="M",at = 5) plot(denM, main="Marginal Density") ## Joint denJ <- dsde2d(mod2, pdf="J", n=100,at = 5) plot(denJ,display="contour",main="Bivariate Transition Density at time t=5") ``` ```{r 103, echo=FALSE, fig.cap='The marginal and joint density of $X_{t}|X_{0}=1,X_{T}=1$ and $Y_{t}|Y_{0}=-0.5,Y_{T}=0.5$ at time $t=5$', fig.env='figure*',fig.width=7,fig.height=7} knitr::include_graphics(c("Figures/fig1009.png","Figures/fig1010.png")) ``` A $3$D plot of the transition density at $t=5$ obtained with: ```{r 11,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} plot(denJ,main="Bivariate Transition Density at time t=5") ``` ```{r 33303, echo=FALSE, fig.cap='$3$D plot of the transition density of $X_{t}|X_{0}=1,X_{T}=1$ and $Y_{t}|Y_{0}=-0.5,Y_{T}=0.5$ at time $t=5$ ', fig.env='figure*',fig.width=5,fig.height=5} knitr::include_graphics("Figures/fig1011.png") ``` We approximate the bivariate transition density over the set transition horizons $t\in [1,9]$ with $\Delta t = 0.005$ using the code: ```{r ,eval=FALSE, include=TRUE} for (i in seq(1,9,by=0.005)){ plot(dsde2d(mod2, at = i,n=100),display="contour",main=paste0('Transition Density \n t = ',i)) } ``` [Return to bridgesde2d()](#bridgesde2d) # `bridgesde3d()` Assume that we want to describe the following bridges SDE's (3D) in Ito form: \begin{equation} \begin{cases} dX_t = -4 (1+X_{t}) Y_{t} dt + 0.2 dW_{1,t},\quad X_{t_{0}}=0 \quad\text{and}\quad X_{T}=0\\ dY_t = 4 (1-Y_{t}) X_{t} dt + 0.2 dW_{2,t},\quad Y_{t_{0}}=-1 \quad\text{and}\quad Y_{T}=-2\\ dZ_t = 4 (1-Z_{t}) Y_{t} dt + 0.2 dW_{3,t},\quad Z_{t_{0}}=0.5 \quad\text{and}\quad Z_{T}=0.5 \end{cases} \end{equation} We simulate a flow of $1000$ trajectories, with integration step size $\Delta t = 0.001$. ```{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) mod3 <- bridgesde3d(x0=c(0,-1,0.5),y=c(0,-2,0.5),drift=fx,diffusion=gx,M=1000) mod3 ``` For plotting (back in time) using the command `plot`, and `plot3D` in space the results of the simulation are shown in Figure 5: ```{r 12,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} plot(mod3) ## in time plot3D(mod3,display = "persp",main="3D Bridge SDE's") ## in space ``` ```{r 3333, echo=FALSE, fig.cap=' Bridge sde 3D ', fig.env='figure*',fig.width=7,fig.height=7} knitr::include_graphics(c("Figures/fig05.png","Figures/fig06.png")) ``` Hence we can just make use of the `rsde3d()` function to build our random number generator for the triplet $X_{t},Y_{t},Z_{t}|X_{0}=0,Y_{0}=-1,Z_{0}=0.5,X_{T}=0,Y_{T}=-2,Z_{T}=0.5$ at time $t=0.75$: ```{r} x3 <- rsde3d(object = mod3, at = 0.75) head(x3, n = 3) ``` the density of $X_{t}|X_{0}=0,X_{T}=0$, $Y_{t}|Y_{0}=-1,Y_{T}=-2$ and $Z_{t}|Z_{0}=0.5,Z_{T}=0.5$ process at time $t=0.75$ are reported using `dsde3d()` function. For an approximate joint density for triplet $X_{t},Y_{t},Z_{t}|X_{0}=0,Y_{0}=-1,Z_{0}=0.5,X_{T}=0,Y_{T}=-2,Z_{T}=0.5$ at time $t=0.75$ (for more details, see package [sm](https://cran.r-project.org/package=sm) or [ks](https://cran.r-project.org/package=ks).) ```{r 15,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE} ## Marginal denM <- dsde3d(mod3,pdf="M",at =0.75) plot(denM, main="Marginal Density") ## Joint denJ <- dsde3d(mod3,pdf="J",at=0.75) plot(denJ,display="rgl") ``` [Return to bridgesde3d()](#bridgesde3d) # 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. Bladt, M. and Sorensen, M. (2007). Simple simulation of diffusion bridges with application to likelihood inference for diffusions. Working Paper, University of Copenhagen. 2. 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.