Constructs and Analysis of Bridges Stochastic Differential Equations

A.C. Guidoum1 and K. Boukhetala2

2020-11-07

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:

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

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:

The function dsde() (where k=1,2,3) approximate transition density for class bridgesdekd() (where k=1,2,3), the main arguments consist:

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")
R> f <- expression((1-x)/(1-t))
R> g <- expression(x)
R> mod <- bridgesde1d(drift=f,diffusion=g,x0=3,y=1,M=1000)
R> mod
Itô Bridge Sde 1D:
    | dX(t) = (1 - X(t))/(1 - t) * dt + X(t) * dW(t)
Method:
    | Euler scheme with order 0.5
Summary:
    | Size of process   | N = 1001.
    | Crossing realized | C = 978 among 1000.
    | Initial value     | x0 = 3.
    | Ending value      | y = 1.
    | Time of process   | t in [0,1].
    | Discretization    | Dt = 0.001.

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> plot(mod,ylab=expression(X[t]))
R> lines(time(mod),apply(mod$X,1,mean),col=2,lwd=2)
R> 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) 
R> head(x, n = 3)
[1] 0.72282 0.96118 0.94990

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

R> dens <- dsde1d(mod, at = 0.55)
R> plot(dens,hist=TRUE) ## histgramme
R> plot(dens,add=TRUE)  ## kernel density

Bridge sde 1D. Histgramme and kernel density estimation for X_{t}|X_{0}=3,X_{T}=1Bridge sde 1D. Histgramme and kernel density estimation for X_{t}|X_{0}=3,X_{T}=1

Return to 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")
R> fx <- expression(-(1+y)*x , -(1+x)*y)
R> gx <- expression(0.2*(1-y),0.1*(1-x))
R> Sigma <-matrix(c(1,0.3,0.3,1),nrow=2,ncol=2)
R> 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)
R> mod2
Stratonovich Bridge Sde 2D:
    | dX(t) = -(1 + Y(t)) * X(t) * dt + 0.2 * (1 - Y(t)) o dB1(t)
    | dY(t) = -(1 + X(t)) * Y(t) * dt + 0.1 * (1 - X(t)) o dB2(t)
    | Correlation structure:                    
             1.0 0.3
             0.3 1.0
Method:
    | Euler scheme with order 0.5
Summary:
    | Size of process   | N = 1001.
    | Crossing realized | C = 1000 among 1000.
    | Initial values    | x0 = (1,-0.5).
    | Ending values     | y = (1,0.5).
    | Time of process   | t in [0,10].
    | Discretization    | Dt = 0.01.

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> plot(mod2,col=c('#FF00004B','#0000FF82'))

Bridge sde 2D

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) 
R> head(x2, n = 3)
         x          y
1 0.051742  0.0095811
2 0.135792  0.0436799
3 0.021494 -0.0348084

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> ## Marginal 
R> denM <- dsde2d(mod2,pdf="M",at = 5)
R> plot(denM, main="Marginal Density")
R> ## Joint
R> denJ <- dsde2d(mod2, pdf="J", n=100,at = 5)
R> plot(denJ,display="contour",main="Bivariate Transition Density at time t=5")

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=5The 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

A \(3\)D plot of the transition density at \(t=5\) obtained with:

R> plot(denJ,main="Bivariate Transition Density at time t=5")

\(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\)

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{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")
R> fx <- expression(-4*(1+x)*y, 4*(1-y)*x, 4*(1-z)*y)
R> gx <- rep(expression(0.2),3)
R> mod3 <- bridgesde3d(x0=c(0,-1,0.5),y=c(0,-2,0.5),drift=fx,diffusion=gx,M=1000)
R> mod3
Itô Bridge 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)
Method:
    | Euler scheme with order 0.5
Summary:
    | Size of process   | N = 1001.
    | Crossing realized | C = 998 among 1000.
    | Initial values    | x0 = (0,-1,0.5).
    | Ending values     | y  = (0,-2,0.5).
    | Time of process   | t in [0,1].
    | Discretization    | Dt = 0.001.

For plotting (back in time) using the command plot, and plot3D in space the results of the simulation are shown in Figure 5:

R> plot(mod3) ## in time
R> plot3D(mod3,display = "persp",main="3D Bridge SDE's") ## in space 

Bridge sde 3DBridge sde 3D

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) 
R> head(x3, n = 3)
       x        y        z
1 2.0636 0.074044 -0.30984
2 1.9445 0.111705 -0.26487
3 1.9288 0.054322 -0.50509

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

R> ## Marginal
R> denM <- dsde3d(mod3,pdf="M",at =0.75)
R> plot(denM, main="Marginal Density")
R> ## Joint
R> denJ <- dsde3d(mod3,pdf="J",at=0.75)
R> plot(denJ,display="rgl")

Return to bridgesde3d()

Further reading

  1. snssdekd() & dsdekd() & rsdekd()- Monte-Carlo Simulation and Analysis of Stochastic Differential Equations.
  2. bridgesdekd() & dsdekd() & rsdekd() - Constructs and Analysis of Bridges Stochastic Differential Equations.
  3. fptsdekd() & dfptsdekd() - Monte-Carlo Simulation and Kernel Density Estimation of First passage time.
  4. MCM.sde() & MEM.sde() - Parallel Monte-Carlo and Moment Equations for SDEs.
  5. TEX.sde() - Converting Sim.DiffProc Objects to LaTeX.
  6. fitsde() - Parametric Estimation of 1-D Stochastic Differential Equation.

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 (2020). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.8, URL https://cran.r-project.org/package=Sim.DiffProc.


  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 ()↩︎