Constructs and Analysis of Bridges Stochastic Differential Equations

A.C. Guidoum1 and K. Boukhetala2

2020-03-27

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,It? 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 It? 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> f <- expression((1-x)/(1-t))
R> g <- expression(x)
R> mod <- bridgesde1d(drift=f,diffusion=g,x0=3,y=1,M=1000,method="milstein")
R> mod
Itô Bridge Sde 1D:
    | dX(t) = (1 - X(t))/(1 - t) * dt + X(t) * dW(t)
Method:
    | First-order Milstein scheme
Summary:
    | Size of process   | N = 1001.
    | Crossing realized | C = 974 among 1000.
    | Initial value     | x0 = 3.
    | Ending value      | y = 1.
    | Time of process   | t in [0,1].
    | Discretization    | Dt = 0.001.
R> summary(mod) ## default: summary at time = (T-t0)/2

Monte-Carlo Statistics for X(t) at time t = 0.5
    | Crossing realized 974 among 1000
                          
Mean               1.98205
Variance           1.48837
Median             1.67933
Mode               1.28981
First quartile     1.15515
Third quartile     2.48642
Minimum            0.38182
Maximum           10.27865
Skewness           1.94209
Kurtosis           8.81024
Coef-variation     0.61552
3th-order moment   3.52644
4th-order moment  19.51694
5th-order moment 107.39652
6th-order moment 688.14117

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")
 Bridge sde 1D

Bridge sde 1D

Figure 2, show approximation results for \(m(t)=\text{E}(X_{t}|X_{0}=3,X_{T}=1)\) and \(S(t)=\text{V}(X_{t}|X_{0}=3,X_{T}=1)\):

R> m  <- apply(mod$X,1,mean) 
R> S  <- apply(mod$X,1,var)
R> out <- data.frame(m,S)
R> matplot(time(mod), out, type = "l", xlab = "time", ylab = "", col=2:3,lwd=2,lty=2:3,las=1)
R> legend("topright",c(expression(m(t),S(t))),col=2:3,lty=2:3,lwd=2,bty="n")

The following statistical measures (S3 method) for class bridgesde1d() can be approximated for the \(X_{t}|X_{0}=3,X_{T}=1\) process at any time \(t\), for example at=0.55:

R> s = 0.55
R> mean(mod, at = s)
[1] 1.9035
R> moment(mod, at = s , center = TRUE , order = 2) ## variance
[1] 1.5899
R> Median(mod, at = s)
[1] 1.5775
R> Mode(mod, at = s)
[1] 1.1897
R> quantile(mod , at = s)
      0%      25%      50%      75%     100% 
 0.32299  1.07949  1.57748  2.35522 12.37977 
R> kurtosis(mod , at = s)
[1] 14.852
R> skewness(mod , at = s)
[1] 2.6276
R> cv(mod , at = s )
[1] 0.66275
R> min(mod , at = s)
[1] 0.32299
R> max(mod , at = s)
[1] 12.38
R> moment(mod, at = s , center= TRUE , order = 4)
[1] 37.619
R> moment(mod, at = s , center= FALSE , order = 4)
[1] 125.48

The result summaries of the \(X_{t}|X_{0}=3,X_{T}=1\) process at time \(t=0.55\):

R> summary(mod, at = 0.55)

Monte-Carlo Statistics for X(t) at time t = 0.55
    | Crossing realized 974 among 1000
                           
Mean                1.90351
Variance            1.59150
Median              1.57748
Mode                1.18973
First quartile      1.07949
Third quartile      2.35522
Minimum             0.32299
Maximum            12.37977
Skewness            2.62755
Kurtosis           14.85244
Coef-variation      0.66275
3th-order moment    5.27547
4th-order moment   37.61937
5th-order moment  292.25260
6th-order moment 2544.14187

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 = s) 
R> head(x, n = 10)
 [1] 2.85901 0.93959 1.48894 0.47225 2.11850 0.72946 0.71184 1.64045
 [9] 1.20608 2.10977
R> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.323   1.079   1.577   1.903   2.355  12.380 

Display the random number generator for \(X_{t}|X_{0}=3,X_{T}=1\), see Figure 3:

R> plot(time(mod),mod$X[,1],type="l",ylab="X(t)",xlab="time",axes=F,lty=3)
R> points(s,x[1],pch=19,col=2,cex=0.5)
R> lines(c(s,s),c(0,x[1]),lty=2,col=2)
R> lines(c(0,s),c(x[1],x[1]),lty=2,col=2)
R> axis(1, s, bquote(at==.(s)),col=2,col.ticks=2)
R> axis(2, x[1], bquote(X[t==.(s)]),col=2,col.ticks=2)
R> legend('topright',col=2,pch=19,legend=bquote(X[t==.(s)]==.(x[1])),bty = 'n')
R> box()

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 = s)
R> dens

 Density of X(t-t0)|X(t0) = 3, X(T) = 1 at time t = 0.55

Data: x (974 obs.); Bandwidth 'bw' = 0.2164

       x                f(x)        
 Min.   :-0.3261   Min.   :0.00000  
 1st Qu.: 3.0126   1st Qu.:0.00095  
 Median : 6.3514   Median :0.00547  
 Mean   : 6.3514   Mean   :0.07480  
 3rd Qu.: 9.6901   3rd Qu.:0.05867  
 Max.   :13.0289   Max.   :0.50888  
R> plot(dens,hist=TRUE) ## histgramme
R> plot(dens,add=TRUE)  ## kernel density

Approximate the transitional densitie of \(X_{t}|X_{0}=3,X_{T}=1\) at \(t-s = \{0.25,0.75\}\):

R> plot(dsde1d(mod,at=0.75))
R> plot(dsde1d(mod,at=0.25),add=TRUE)
R> legend('topright',col=c('#0000FF4B','#FF00004B'),pch=15,legend=c("t-s=0.25","t-s=0.75"),bty = 'n')
 Transitional densitie at time $t-s = 0.25,0.75$

Transitional densitie at time \(t-s = 0.25,0.75\)

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 dW_{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 dW_{2,t},\quad Y_{t_{0}}=-0.5 \quad\text{and}\quad Y_{T}=0.5 \end{cases} \end{equation}\]

We simulate a flow of \(1000\) trajectories, with integration step size \(\Delta t = 0.01\), and using Runge-Kutta method order 1:

R> fx <- expression(-(1+y)*x , -(1+x)*y)
R> gx <- expression(0.2*(1-y),0.1*(1-x))
R> mod2 <- bridgesde2d(drift=fx,diffusion=gx,x0=c(1,-0.5),y=c(1,0.5),Dt=0.01,M=1000,type="str",method="rk1")
R> mod2
Stratonovich Bridge Sde 2D:
    | dX(t) = -(1 + Y(t)) * X(t) * dt + 0.2 * (1 - Y(t)) o dW1(t)
    | dY(t) = -(1 + X(t)) * Y(t) * dt + 0.1 * (1 - X(t)) o dW2(t)
Method:
    | Runge-Kutta method with order 1
Summary:
    | Size of process   | N = 1001.
    | Crossing realized | C = 996 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.
R> summary(mod2) ## default: summary at time = (T-t0)/2

Monte-Carlo Statistics for (X(t),Y(t)) at time t = 5
    | Crossing realized 996 among 1000
                        X        Y
Mean              0.00687  0.00304
Variance          0.02097  0.00521
Median            0.00925  0.00291
Mode              0.01899  0.00049
First quartile   -0.08950 -0.04446
Third quartile    0.10620  0.05169
Minimum          -0.48106 -0.26840
Maximum           0.45477  0.27567
Skewness         -0.01526 -0.05365
Kurtosis          2.99882  3.31363
Coef-variation   21.08027 23.77779
3th-order moment -0.00005 -0.00002
4th-order moment  0.00132  0.00009
5th-order moment -0.00001  0.00000
6th-order moment  0.00014  0.00000

In Figure 6, 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

Bridge sde 2D

Figure 7, show approximation results for \(m_{1}(t)=\text{E}(X_{t}|X_{0}=1,X_{T}=1)\), \(m_{2}(t)=\text{E}(Y_{t}|Y_{0}=-0.5,Y_{T}=0.5)\),and \(S_{1}(t)=\text{V}(X_{t}|X_{0}=1,X_{T}=1)\), \(S_{2}(t)=\text{V}(Y_{t}|Y_{0}=-0.5,Y_{T}=0.5)\), and \(C_{12}(t)=\text{COV}(X_{t},Y_{t}|X_{0}=1,Y_{0}=-0.5,X_{T}=1,Y_{T}=0.5)\):

R> m1  <- apply(mod2$X,1,mean) 
R> m2  <- apply(mod2$Y,1,mean) 
R> S1  <- apply(mod2$X,1,var)
R> S2  <- apply(mod2$Y,1,var)
R> C12 <- sapply(1:dim(mod2$X)[1],function(i) cov(mod2$X[i,],mod2$Y[i,]))
R> out2 <- data.frame(m1,m2,S1,S2,C12)
R> matplot(time(mod2), out2, type = "l", xlab = "time", ylab = "", col=2:6,lwd=2,lty=2:6,las=1)
R> legend("top",c(expression(m[1](t),m[2](t),S[1](t),S[2](t),C[12](t))),col=2:6,lty=2:6,lwd=2,bty="n")

The following statistical measures (S3 method) for class bridgesde2d() can be approximated for the \(X_{t}|X_{0}=1,X_{T}=1\) and \(Y_{t}|Y_{0}=-0.5,Y_{T}=0.5\) process at any time \(t\), for example at=6.75:

R> s = 6.75
R> mean(mod2, at = s)
[1] 0.0287775 0.0091243
R> moment(mod2, at = s , center = TRUE , order = 2) ## variance
[1] 0.0193018 0.0045472
R> Median(mod2, at = s)
[1] 0.0276037 0.0094694
R> Mode(mod2, at = s)
[1] 0.0073237 0.0145955
R> quantile(mod2 , at = s)
$x
       0%       25%       50%       75%      100% 
-0.483710 -0.067725  0.027604  0.122543  0.473686 

$y
        0%        25%        50%        75%       100% 
-0.1956590 -0.0347267  0.0094694  0.0554088  0.2198906 
R> kurtosis(mod2 , at = s)
[1] 3.0214 3.0475
R> skewness(mod2 , at = s)
[1] 0.029355 0.018179
R> cv(mod2 , at = s )
[1] 4.8302 7.3942
R> min(mod2 , at = s)
[1] -0.48371 -0.19566
R> max(mod2 , at = s)
[1] 0.47369 0.21989
R> moment(mod2 , at = s , center= TRUE , order = 4)
[1] 0.001127918 0.000063139
R> moment(mod2 , at = s , center= FALSE , order = 4)
[1] 0.001233587 0.000065621

The result summaries of the \(X_{t}|X_{0}=1,X_{T}=1\) and \(Y_{t}|Y_{0}=-0.5,Y_{T}=0.5\) process at time \(t=6.75\):

R> summary(mod2, at = 6.75)

Monte-Carlo Statistics for (X(t),Y(t)) at time t = 6.75
    | Crossing realized 996 among 1000
                        X        Y
Mean              0.02878  0.00912
Variance          0.01932  0.00455
Median            0.02760  0.00947
Mode              0.00732  0.01460
First quartile   -0.06773 -0.03473
Third quartile    0.12254  0.05541
Minimum          -0.48371 -0.19566
Maximum           0.47369  0.21989
Skewness          0.02935  0.01818
Kurtosis          3.02141  3.04746
Coef-variation    4.83019  7.39417
3th-order moment  0.00008  0.00001
4th-order moment  0.00113  0.00006
5th-order moment  0.00001  0.00000
6th-order moment  0.00011  0.00000

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=6.75\):

R> x2 <- rsde2d(object = mod2, at = s) 
R> head(x2, n = 10)
           x         y
1   0.194766 -0.080807
2  -0.115444 -0.067434
3   0.061045 -0.114263
4  -0.211655 -0.073105
5  -0.179774  0.032211
6   0.042803 -0.073156
7   0.102801 -0.087403
8   0.177113 -0.065050
9   0.448498  0.021871
10 -0.084539 -0.042989
R> summary(x2)
       x                 y           
 Min.   :-0.4837   Min.   :-0.19566  
 1st Qu.:-0.0677   1st Qu.:-0.03473  
 Median : 0.0276   Median : 0.00947  
 Mean   : 0.0288   Mean   : 0.00912  
 3rd Qu.: 0.1225   3rd Qu.: 0.05541  
 Max.   : 0.4737   Max.   : 0.21989  

Display the random number generator for the couple \(X_{t},Y_{t}|X_{0}=1,Y_{0}=-0.5,X_{T}=1,Y_{T}=0.5\), see Figure 8:

R> plot(ts.union(mod2$X[,1],mod2$Y[,1]),col=1:2,lty=3,plot.type="single",type="l",ylab= "",xlab="time",axes=F)
R> points(s,x2$x[1],pch=19,col=3,cex=0.8)
R> points(s,x2$y[1],pch=19,col=4,cex=0.8)
R> lines(c(s,s),c(-10,x2$x[1]),lty=2,col=6)
R> lines(c(0,s),c(x2$x[1],x2$x[1]),lty=2,col=3)
R> lines(c(0,s),c(x2$y[1],x2$y[1]),lty=2,col=4)
R> axis(1, s, bquote(at==.(s)),col=6,col.ticks=6)
R> axis(2, x2$x[1], bquote(X[t==.(s)]),col=3,col.ticks=3)
R> axis(2, x2$y[1], bquote(Y[t==.(s)]),col=4,col.ticks=4)
R> legend('topright',legend=bquote(c(X[t==.(s)]==.(x2$x[1]),Y[t==.(s)]==.(x2$y[1]))),bty = 'n')
R> box()

For each SDE type and for each numerical scheme, the density of \(X_{t}|X_{0}=1,X_{T}=1\) and \(Y_{t}|Y_{0}=-0.5,Y_{T}=0.5\) at time \(t=6.75\) are reported using dsde2d() function, see e.g. Figure 9:

R> denM <- dsde2d(mod2,pdf="M",at =s)
R> denM

 Marginal density of X(t-t0)|X(t0) = 1, X(T) = 1 at time t = 6.75

Data: x (996 obs.); Bandwidth 'bw' = 0.03145

       x                 f(x)        
 Min.   :-0.57806   Min.   :0.00014  
 1st Qu.:-0.29153   1st Qu.:0.02891  
 Median :-0.00501   Median :0.32700  
 Mean   :-0.00501   Mean   :0.87168  
 3rd Qu.: 0.28151   3rd Qu.:1.70604  
 Max.   : 0.56803   Max.   :2.77730  

 Marginal density of Y(t-t0)|Y(t0) = -0.5, Y(T) = 0.5 at time t = 6.75

Data: y (996 obs.); Bandwidth 'bw' = 0.01522

       y                  f(y)       
 Min.   :-0.241315   Min.   :0.0004  
 1st Qu.:-0.114600   1st Qu.:0.1648  
 Median : 0.012116   Median :1.0318  
 Mean   : 0.012116   Mean   :1.9710  
 3rd Qu.: 0.138831   3rd Qu.:3.6994  
 Max.   : 0.265547   Max.   :5.9032  
R> plot(denM, main="Marginal Density")

Created using dsde2d() plotted in (x, y)-space with dim = 2. A contour and image plot of 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=6.75.

R> denJ <- dsde2d(mod2, pdf="J", n=100,at =s)
R> denJ

 Joint density of (X(t-t0),Y(t-t0)|X(t0)=1,Y(t0)=-0.5,X(T)=1,Y(T)=0.5) at time t = 6.75

Data: (x,y) (2 x 996 obs.);

       x                  y                 f(x,y)       
 Min.   :-0.48371   Min.   :-0.195659   Min.   : 0.0000  
 1st Qu.:-0.24436   1st Qu.:-0.091772   1st Qu.: 0.1040  
 Median :-0.00501   Median : 0.012116   Median : 0.6052  
 Mean   :-0.00501   Mean   : 0.012116   Mean   : 2.4535  
 3rd Qu.: 0.23434   3rd Qu.: 0.116003   3rd Qu.: 3.2169  
 Max.   : 0.47369   Max.   : 0.219891   Max.   :15.5664  
R> plot(denJ,display="contour",main="Bivariate Transition Density at time t=6.755")
R> plot(denJ,display="image",main="Bivariate Transition Density at time t=6.755")

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

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

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 It? 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> 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 = 1000 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.
R> summary(mod3) ## default: summary at time = (T-t0)/2

Monte-Carlo Statistics for (X(t),Y(t),Z(t)) at time t = 0.5
    | Crossing realized 1000 among 1000
                       X        Y        Z
Mean             0.69964  0.50782  0.10111
Variance         0.00888  0.00727  0.01506
Median           0.69431  0.50637  0.11076
Mode             0.66826  0.48946  0.13381
First quartile   0.63724  0.44969  0.02482
Third quartile   0.76375  0.56653  0.18193
Minimum          0.41309  0.20705 -0.36741
Maximum          0.99111  0.73060  0.45619
Skewness         0.12460 -0.03772 -0.32465
Kurtosis         3.01795  2.84608  3.41476
Coef-variation   0.13469  0.16790  1.21361
3th-order moment 0.00010 -0.00002 -0.00060
4th-order moment 0.00024  0.00015  0.00077
5th-order moment 0.00001  0.00000 -0.00010
6th-order moment 0.00001  0.00000  0.00007

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

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

  Bridge sde 3D

Bridge sde 3D

Figure 13, show approximation results for \(m_{1}(t)=\text{E}(X_{t}|X_{0}=0,X_{T}=0)\), \(m_{2}(t)=\text{E}(Y_{t}|Y_{0}=-1,Y_{T}=-2)\), \(m_{3}(t)=\text{E}(Z_{t}|Z_{0}=0.5,Z_{T}=0.5)\) and \(S_{1}(t)=\text{V}(X_{t}|X_{0}=0,X_{T}=0)\), \(S_{2}(t)=\text{V}(Y_{t}|Y_{0}=-1,Y_{T}=-2)\), \(S_{3}(t)=\text{V}(Z_{t}|Z_{0}=0.5,Z_{T}=0.5)\),

R> m1  <- apply(mod3$X,1,mean) 
R> m2  <- apply(mod3$Y,1,mean) 
R> m3  <- apply(mod3$Z,1,mean) 
R> S1  <- apply(mod3$X,1,var)
R> S2  <- apply(mod3$Y,1,var)
R> S3  <- apply(mod3$Z,1,var)
R> out3 <- data.frame(m1,m2,m3,S1,S2,S3)
R> matplot(time(mod3), out3, type = "l", xlab = "time", ylab = "", col=2:7,lwd=2,lty=2:7,las=1)
R> legend("bottom",c(expression(m[1](t),m[2](t),m[3](t),S[1](t),S[2](t),S[3](t))),col=2:7,lty=2:7,lwd=2,bty="n")

The following statistical measures (S3 method) for class bridgesde3d() can be approximated for the \(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 any time \(t\), for example at=0.75:

R> s = 0.75
R> mean(mod3, at = s)
[1]  1.99472  0.12727 -0.49999
R> moment(mod3, at = s , center = TRUE , order = 2) ## variance
[1] 0.0112535 0.0046449 0.0299859
R> Median(mod3, at = s)
[1]  1.99507  0.12781 -0.49998
R> Mode(mod3, at = s)
[1]  1.98686  0.12874 -0.48712
R> quantile(mod3 , at = s)
$x
    0%    25%    50%    75%   100% 
1.6845 1.9229 1.9951 2.0680 2.4141 

$y
       0%       25%       50%       75%      100% 
-0.154991  0.083464  0.127811  0.172711  0.334215 

$z
      0%      25%      50%      75%     100% 
-1.03031 -0.61367 -0.49998 -0.37916  0.14141 
R> kurtosis(mod3 , at = s)
[1] 2.9072 3.2388 2.9359
R> skewness(mod3 , at = s)
[1]  0.0184680 -0.1088275  0.0026912
R> cv(mod3 , at = s )
[1]  0.053208  0.535778 -0.346506
R> min(mod3 , at = s)
[1]  1.68447 -0.15499 -1.03031
R> max(mod3 , at = s)
[1] 2.41405 0.33422 0.14141
R> moment(mod3 , at = s , center= TRUE , order = 4)
[1] 0.000368914 0.000070017 0.002645087
R> moment(mod3 , at = s , center= FALSE , order = 4)
[1] 16.1007584  0.0007662  0.1100921

The result summaries of the \(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\):

R> summary(mod3, at = 0.75)

Monte-Carlo Statistics for (X(t),Y(t),Z(t)) at time t = 0.75
    | Crossing realized 1000 among 1000
                       X        Y        Z
Mean             1.99472  0.12727 -0.49999
Variance         0.01126  0.00465  0.03002
Median           1.99507  0.12781 -0.49998
Mode             1.98686  0.12874 -0.48712
First quartile   1.92294  0.08346 -0.61367
Third quartile   2.06803  0.17271 -0.37916
Minimum          1.68447 -0.15499 -1.03031
Maximum          2.41405  0.33422  0.14141
Skewness         0.01847 -0.10883  0.00269
Kurtosis         2.90722  3.23885  2.93586
Coef-variation   0.05321  0.53578 -0.34651
3th-order moment 0.00002 -0.00003  0.00001
4th-order moment 0.00037  0.00007  0.00265
5th-order moment 0.00001  0.00000  0.00006
6th-order moment 0.00002  0.00000  0.00039

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 = s) 
R> head(x3, n = 10)
        x        y        z
1  2.1601 0.123274 -0.38506
2  1.8655 0.219702 -0.42943
3  2.0614 0.175151 -0.53321
4  2.0773 0.239680 -0.62917
5  2.0383 0.131557 -0.50760
6  1.8922 0.014355 -0.38828
7  1.9735 0.158596 -0.46338
8  1.8926 0.100864 -0.49191
9  1.9648 0.124344 -0.46666
10 2.0839 0.165809 -0.94800
R> summary(x3)
       x              y                 z         
 Min.   :1.68   Min.   :-0.1550   Min.   :-1.030  
 1st Qu.:1.92   1st Qu.: 0.0835   1st Qu.:-0.614  
 Median :2.00   Median : 0.1278   Median :-0.500  
 Mean   :1.99   Mean   : 0.1273   Mean   :-0.500  
 3rd Qu.:2.07   3rd Qu.: 0.1727   3rd Qu.:-0.379  
 Max.   :2.41   Max.   : 0.3342   Max.   : 0.141  

Display the random number generator 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\): , see Figure 14:

R> plot(ts.union(mod3$X[,1],mod3$Y[,1],mod3$Z[,1]),col=1:3,lty=3,plot.type="single",type="l",ylab= "",xlab="time",axes=F)
R> points(s,x3$x[1],pch=19,col=4,cex=0.8)
R> points(s,x3$y[1],pch=19,col=5,cex=0.8)
R> points(s,x3$z[1],pch=19,col=6,cex=0.8)
R> lines(c(s,s),c(-10,x3$x[1]),lty=2,col=7)
R> lines(c(0,s),c(x3$x[1],x3$x[1]),lty=2,col=4)
R> lines(c(0,s),c(x3$y[1],x3$y[1]),lty=2,col=5)
R> lines(c(0,s),c(x3$z[1],x3$z[1]),lty=2,col=6)
R> axis(1, s, bquote(at==.(s)),col=7,col.ticks=7)
R> axis(2, x3$x[1], bquote(X[t==.(s)]),col=4,col.ticks=4)
R> axis(2, x3$y[1], bquote(Y[t==.(s)]),col=5,col.ticks=5)
R> axis(2, x3$z[1], bquote(Z[t==.(s)]),col=6,col.ticks=6)
R> legend("bottomleft",legend=bquote(c(X[t==.(s)]==.(x3$x[1]),Y[t==.(s)]==.(x3$y[1]),Z[t==.(s)]==.(x3$z[1]))),bty = 'n',cex=0.75)
R> box()

For each SDE type and for each numerical scheme, 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, see e.g. Figure 15:

R> denM <- dsde3d(mod3,pdf="M",at =s)
R> denM

 Marginal density of X(t-t0)|X(t0) = 0, X(T) = 0 at time t = 0.75

Data: x (1000 obs.);    Bandwidth 'bw' = 0.02399

       x               f(x)       
 Min.   :1.6125   Min.   :0.0002  
 1st Qu.:1.8309   1st Qu.:0.0225  
 Median :2.0493   Median :0.5061  
 Mean   :2.0493   Mean   :1.1436  
 3rd Qu.:2.2677   3rd Qu.:2.1594  
 Max.   :2.4860   Max.   :3.6310  

 Marginal density of Y(t-t0)|Y(t0) = -1, Y(T) = -2 at time t = 0.75

Data: y (1000 obs.);    Bandwidth 'bw' = 0.01506

       y                 f(y)       
 Min.   :-0.20016   Min.   :0.0003  
 1st Qu.:-0.05527   1st Qu.:0.0541  
 Median : 0.08961   Median :0.6216  
 Mean   : 0.08961   Mean   :1.7238  
 3rd Qu.: 0.23450   3rd Qu.:3.2131  
 Max.   : 0.37939   Max.   :6.0600  

 Marginal density of Z(t-t0)|Z(t0) = 0.5, Z(T) = 0.5 at time t = 0.75

Data: z (1000 obs.);    Bandwidth 'bw' = 0.03917

       z                 f(z)        
 Min.   :-1.14781   Min.   :0.00012  
 1st Qu.:-0.79613   1st Qu.:0.02155  
 Median :-0.44445   Median :0.34207  
 Mean   :-0.44445   Mean   :0.71017  
 3rd Qu.:-0.09277   3rd Qu.:1.43614  
 Max.   : 0.25891   Max.   :2.16290  
R> plot(denM, main="Marginal Density")

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> 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.5, 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 ()