# Constructs and Analysis of Bridges Stochastic Differential Equations" In Sim.DiffProc: Simulation of Diffusion Processes

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: $$\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$$ 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$. 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$: 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$:

x <- rsde1d(object = mod, at = 0.55)


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 package):

dens <- dsde1d(mod, at = 0.55)
plot(dens,hist=TRUE) ## histgramme


r|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: \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} 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$: plot(mod2,col=c('#FF00004B','#0000FF82'))  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$: 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: ## 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|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
plot(denJ,main="Bivariate Transition Density at time t=5")


r|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 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() # bridgesde3d() Assume that we want to describe the following bridges SDE's (3D) in Ito form: $$\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}$$ We simulate a flow of$1000$trajectories, with integration step size$\Delta t = 0.001$. 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: plot(mod3) ## in time plot3D(mod3,display = "persp",main="3D Bridge SDE's") ## in space  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$: 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 or ks.)

## 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")


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

## Try the Sim.DiffProc package in your browser

Any scripts or data that you put into this service are public.

Sim.DiffProc documentation built on May 29, 2024, 8:09 a.m.