Bridges SDE’s

Consider now a d-dimensional stochastic process X_t defined on a probability space (Ω, 𝔉, ℙ). We say that the bridge associated to X_t conditioned to the event {X_T = a} is the process: {X̃_t, t₀ ≤ t ≤ T} = {X_t, t₀ ≤ t ≤ T : X_T = a} where T is a deterministic fixed time and a ∈ ℝ^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 ∈ [t₀, T] (default: at=(T-t0)/2):

• The expected value E(X_t) at time t, using the command mean.
• The variance Var(X_t) at time t, using the command moment with order=2 and center=TRUE.
• The median Med(X_t) at time t, using the command Median.
• The mode 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 ∈ [t₀, 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: We simulate a flow of 1000 trajectories, with integration step size Δt = 0.001, and x₀ = 3 at time t₀ = 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₀ = 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₀ = 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₀ = 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

Return to bridgesde1d()

bridgesde2d()

Assume that we want to describe the following 2-dimensional bridge SDE’s in Stratonovich form:

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 Δ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₀ = 1, X_T = 1 and Y_t|Y₀ = −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₀ = 1, Y₀ = −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₀ = 1, X_T = 1 and Y_t|Y₀ = −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₀ = 1, Y₀ = −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")

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

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

3D plot of the transition density of X_t|X₀ = 1, X_T = 1 and Y_t|Y₀ = −0.5, Y_T = 0.5 at time t = 5

We approximate the bivariate transition density over the set transition horizons t ∈ [1, 9] with Δ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:

We simulate a flow of 1000 trajectories, with integration step size Δ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 

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, Y₀ = −1, Z₀ = 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, X_T = 0, Y_t|Y₀ = −1, Y_T = −2 and Z_t|Z₀ = 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, Y₀ = −1, Z₀ = 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 (2024). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.9, URL https://cran.r-project.org/package=Sim.DiffProc.