Parallel Monte-Carlo and Moment Equations for SDEs

A.C. Guidoum1 and K. Boukhetala2

2020-11-07

The MCM.sde() function

R> MCM.sde(model, statistic, R = 1000, time, exact = NULL, names = NULL,level = 0.95, 
+         parallel = c("no", "multicore", "snow"),ncpus = getOption("ncpus", 1L), cl = NULL, ...)

The main arguments of MCM.sde() function in Sim.DiffProc package consist:

R> plot(x,index = 1,type=c("all","hist","qqplot","boxplot","CI"), ...)

This takes a MCM.sde() object and produces plots for the R replicates of the interesting quantity.

One-dimensional SDE

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
R> theta = 0.75; x0 = 1
R> fx <- expression( 0.5*theta^2*x )
R> gx <- expression( theta*x )
R> mod1 <- snssde1d(drift=fx,diffusion=gx,x0=x0,M=500,type="ito")
R> mod2 <- snssde1d(drift=fx,diffusion=gx,x0=x0,M=500,type="str")
R> ## True values of means and variance for mod1 and mod2
R> E.mod1 <- function(t) x0 * exp(0.5 * theta^2 * t)
R> V.mod1 <- function(t) x0^2 * exp(theta^2 * t) * (exp(theta^2 * t) - 1)
R> E.mod2 <- function(t) x0 * exp(theta^2 * t)
R> V.mod2 <- function(t) x0^2 * exp(2 * theta^2 * t) * (exp(theta^2 * t) - 1)
R> ## function of the statistic(s) of interest.
R> sde.fun1d <- function(data, i){
+      d <- data[i, ]
+      return(c(mean(d),var(d)))
+ }
R> # Parallel MOnte Carlo for mod1
R> mcm.mod1 = MCM.sde(model=mod1,statistic=sde.fun1d,R=20, exact=list(m=E.mod1(1),S=V.mod1(1)),parallel="snow",ncpus=2)
R> mcm.mod1
Itô Sde 1D:
 | dX(t) = 0.5 * theta^2 * X(t) * dt + theta * X(t) * dW(t)
 | t in [0,1] with mesh equal to 0.001

PMCM Based on 20 Batches with 500-Realisations at time 1:

   Exact Estimate    Bias Std.Error    RMSE  CI( 2.5 % , 97.5 % )
m 1.3248   1.3505 0.02571   0.01070 0.05327 ( 1.32952 , 1.37146 )
S 1.3252   1.3326 0.00742   0.03637 0.15872  ( 1.2613 , 1.40386 )
R> # Parallel MOnte Carlo for mod2
R> mcm.mod2 = MCM.sde(model=mod2,statistic=sde.fun1d,R=20, exact=list(m=E.mod2(1),S=V.mod2(1)),parallel="snow",ncpus=2)
R> mcm.mod2
Stratonovich Sde 1D:
 | dX(t) = 0.5 * theta^2 * X(t) * dt + theta * X(t) o dW(t)
 | t in [0,1] with mesh equal to 0.001

PMCM Based on 20 Batches with 500-Realisations at time 1:

   Exact Estimate    Bias Std.Error    RMSE  CI( 2.5 % , 97.5 % )
m 1.7550   1.7889 0.03383   0.01418 0.07045 ( 1.76109 , 1.81667 )
S 2.3257   2.3365 0.01081   0.06376 0.27812 ( 2.21157 , 2.46151 )

Two-dimensional SDEs

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
R> mu=1;sigma=0.5;theta=2
R> x0=0;y0=0;init=c(x0,y0)
R> f <- expression(1/mu*(theta-x), x)  
R> g <- expression(sqrt(sigma),0)
R> OUI <- snssde2d(drift=f,diffusion=g,M=500,Dt=0.015,x0=c(x=0,y=0))
R> ## true values of first and second moment at time 10
R> Ex <- function(t) theta+(x0-theta)*exp(-t/mu)
R> Vx <- function(t) 0.5*sigma*mu *(1-exp(-2*(t/mu)))
R> Ey <- function(t) y0+theta*t+(x0-theta)*mu*(1-exp(-t/mu))
R> Vy <- function(t) sigma*mu^3*((t/mu)-2*(1-exp(-t/mu))+0.5*(1-exp(-2*(t/mu))))
R> covxy <- function(t) 0.5*sigma*mu^2 *(1-2*exp(-t/mu)+exp(-2*(t/mu)))
R> tvalue = list(m1=Ex(10),m2=Ey(10),S1=Vx(10),S2=Vy(10),C12=covxy(10))
R> ## function of the statistic(s) of interest.
R> sde.fun2d <- function(data, i){
+   d <- data[i,]
+   return(c(mean(d$x),mean(d$y),var(d$x),var(d$y),cov(d$x,d$y)))
+ }
R> ## Parallel Monte-Carlo of 'OUI' at time 10
R> mcm.mod2d = MCM.sde(OUI,statistic=sde.fun2d,time=10,R=20,exact=tvalue,parallel="snow",ncpus=2)
R> mcm.mod2d
Itô Sde 2D:
 | dX(t) = 1/mu * (theta - X(t)) * dt + sqrt(sigma) * dW1(t)
 | dY(t) = X(t) * dt + 0 * dW2(t)
 | t in [0,15] with mesh equal to 0.015

PMCM Based on 20 Batches with 500-Realisations at time 10:

       Exact Estimate    Bias Std.Error    RMSE
m1   1.99991  2.00256 0.00265   0.00475 0.02087
m2  18.00009 18.04024 0.04015   0.01598 0.08038
S1   0.25000  0.25229 0.00229   0.00397 0.01746
S2   4.25005  4.29577 0.04572   0.05856 0.25934
C12  0.24998  0.25877 0.00879   0.01266 0.05588
       CI( 2.5 % , 97.5 % )
m1    ( 1.99325 , 2.01187 )
m2  ( 18.00892 , 18.07156 )
S1    ( 0.24451 , 0.26007 )
S2    ( 4.18099 , 4.41055 )
C12   ( 0.23396 , 0.28358 )

Three-dimensional SDEs

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
R> mu=0.5;sigma=0.25
R> fx <- expression(mu*y,0,0) 
R> gx <- expression(sigma*z,1,1)
R> Sigma <-matrix(c(1,0.3,-0.5,0.3,1,0.2,-0.5,0.2,1),nrow=3,ncol=3)
R> modtra <- snssde3d(drift=fx,diffusion=gx,M=500,type="str",corr=Sigma)
R> ## function of the statistic(s) of interest.
R> sde.fun3d <- function(data, i){
+   d <- data[i,]
+   return(c(mean(d$x),median(d$x),Mode(d$x)))
+ }
R> ## Monte-Carlo at time = 10
R> mcm.mod3d = MCM.sde(modtra,statistic=sde.fun3d,R=10,parallel="snow",ncpus=2)
R> mcm.mod3d
Stratonovich Sde 3D:
 | dX(t) = mu * Y(t) * dt + sigma * Z(t) o dB1(t)
 | dY(t) = 0 * dt + 1 o dB2(t)
 | dZ(t) = 0 * dt + 1 o dB3(t)
 | t in [0,1] with mesh equal to 0.001
 | Correlation structure:                    
        1.0 0.3 -0.5
        0.3 1.0  0.2
       -0.5 0.2  1.0

PMCM Based on 10 Batches with 500-Realisations at time 1:

    Estimate Std.Error    CI( 2.5 % , 97.5 % )
mu1 -0.06544   0.00325 ( -0.07181 , -0.05907 )
mu2 -0.05929   0.00555 ( -0.07017 , -0.04841 )
mu3 -0.04464   0.01661  ( -0.0772 , -0.01208 )

The MEM.sde() function

R> MEM.sde(drift, diffusion, corr = NULL, type = c("ito", "str"), solve = FALSE, parms = NULL, init = NULL, time = NULL, ...)

The main arguments of MEM.sde() function in Sim.DiffProc package consist:

One-dimensional SDE

R> fx <- expression( 0.5*theta^2*x )
R> gx <- expression( theta*x )
R> start = c(m=1,S=0)
R> t = seq(0,1,by=0.001)
R> mem.mod1 = MEM.sde(drift=fx,diffusion=gx,type="ito",solve = TRUE,parms = c(theta=0.75), init = start, time = t)
R> mem.mod1
Itô Sde 1D:
 | dX(t) = 0.5 * 0.75^2 * X(t) * dt + 0.75 * X(t) * dW(t)
 | t in [0,1].

Moment equations: 
 | dm(t) = 0.28125 * m(t)
 | dS(t) = 0.5625 * m(t)^2 + 1.125 * S(t)

Approximation of moment at time 1
 | m(1) = 1.3248
 | S(1) = 1.3252
R> mem.mod2 = MEM.sde(drift=fx,diffusion=gx,type="str",solve = TRUE,parms = c(theta=0.75), init = start, time = t)
R> mem.mod2
Stratonovich Sde 1D:
 | dX(t) = 0.5 * 0.75^2 * X(t) * dt + 0.75 * X(t) o dW(t)
 | t in [0,1].

Moment equations: 
 | dm(t) = 0.5625 * m(t)
 | dS(t) = 0.5625 * m(t)^2 + 1.6875 * S(t)

Approximation of moment at time 1
 | m(1) = 1.755
 | S(1) = 2.3257
R> plot(mem.mod1$sol.ode, mem.mod2$sol.ode,ylab = c("m(t)"),select="m", xlab = "Time",main="",col = 2:3,lty=1)
R> legend("topleft",c(expression(m[mod1](t),m[mod2](t))),inset = .05, col=2:3,lty=1)
R> plot(mem.mod1$sol.ode, mem.mod2$sol.ode,ylab = c("S(t)"),select="S", xlab = "Time",main="",col = 2:3,lty=1)
R> legend("topleft",c(expression(S[mod1](t),S[mod2](t))),inset = .05, col=2:3,lty=1)

Two-dimensional SDEs

R> fx <- expression(1/mu*(theta-x), x)  
R> gx <- expression(sqrt(sigma),0)
R> start = c(m1=0,m2=0,S1=0,S2=0,C12=0)
R> t = seq(0,10,by=0.001)
R> mem.mod2d = MEM.sde(drift=fx,diffusion=gx,type="ito",solve = TRUE,parms = c(mu=1,sigma=0.5,theta=2), init = start, time = t)
R> mem.mod2d
Itô Sde 2D:
 | dX(t) = 1/1 * (2 - X(t)) * dt + sqrt(0.5) * dW1(t)
 | dY(t) = X(t) * dt + 0 * dW2(t)
 | t in [0,10].

Moment equations: 
 | dm1(t)  = 2 - m1(t)
 | dm2(t)  = m1(t)
 | dS1(t)  = 0.5 - 2 * S1(t)
 | dS2(t)  = 2 * C12(t)
 | dC12(t) = S1(t) - C12(t)

Approximation of moment at time 10                                                              
  | m1(10)  =   1.9999 | S1(10)  =  0.25 | C12(10)  =  0.24998
  | m2(10)  =  18.0001 | S2(10)  =  4.25                      
R> matplot.0D(mem.mod2d$sol.ode,main="")

Three-dimensional SDEs

R> fx <- expression(mu*y,0,0) 
R> gx <- expression(sigma*z,1,1)
R> RHO <- expression(0.75,0.5,-0.25)
R> start = c(m1=5,m2=0,m3=0,S1=0,S2=0,S3=0,C12=0,C13=0,C23=0)
R> t = seq(0,1,by=0.001)
R> mem.mod3d = MEM.sde(drift=fx,diffusion=gx,corr=RHO,type="ito",solve = TRUE,parms = c(mu=0.5,sigma=0.25), init = start, time = t)
R> mem.mod3d
Itô Sde 3D:
 | dX(t) = 0.5 * Y(t) * dt + 0.25 * Z(t) * dB1(t)
 | dY(t) = 0 * dt + 1 * dB2(t)
 | dZ(t) = 0 * dt + 1 * dB3(t)
 | t in [0,1].
 | Correlation structure: E(dB1dB2) = 0.75 * dt
                        : E(dB1dB3) = 0.5 * dt
                        : E(dB2dB3) = -0.25 * dt

Moment equations: 
 | dm1(t)  = 0.5 * m2(t)
 | dm2(t)  = 0
 | dm3(t)  = 0
 | dS1(t)  = 0.0625 * S3(t) + 0.0625 * m3(t)^2 + C12(t)
 | dS2(t)  = 1
 | dS3(t)  = 1
 | dC12(t) = 0.1875 * m3(t) + 0.5 * S2(t)
 | dC13(t) = 0.125 * m3(t) + 0.5 * C23(t)
 | dC23(t) = -0.25

Approximation of moment at time 1                                                         
   | m1(1)  =  5 | S1(1)  =  0.11458 | C12(1)  =   0.2500
   | m2(1)  =  0 | S2(1)  =  1.00000 | C13(1)  =  -0.0625
   | m3(1)  =  0 | S3(1)  =  1.00000 | C23(1)  =  -0.2500
R> matplot.0D(mem.mod3d$sol.ode,main="",select=c("m1","m2","m3"))
R> matplot.0D(mem.mod3d$sol.ode,main="",select=c("S1","S2","S3"))
R> matplot.0D(mem.mod3d$sol.ode,main="",select=c("C12","C13","C23"))

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. Guidoum AC, Boukhetala K (2020). “Performing Parallel Monte Carlo and Moment Equations Methods for Itô and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc”. Journal of Statistical Software, 96(2), 1–82. https://doi.org/10.18637/jss.v096.i02

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