Parallel Monte-Carlo and Moment Equations for SDEs

A.C. Guidoum1 and K. Boukhetala2

2017-09-16

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> 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> E.mod2 <- function(t) x0
R> V.mod1 <- function(t) x0^2*exp(theta^2*t)*(exp(theta^2*t)-1)
R> V.mod2 <- function(t) x0^2 *(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=100, 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].

MCM Based on 100 Batches of 500-Realisations at time 1

   Exact Estimate    Bias Std.Error    RMSE  CI( 2.5 % , 97.5 % )
m 1.3248   1.3210 0.00380   0.00467 0.04643 ( 1.31183 , 1.33013 )
S 1.3252   1.2792 0.04601   0.03231 0.32147 ( 1.21582 , 1.34248 )
R> # Parallel MOnte Carlo for mod2
R> mcm.mod2 = MCM.sde(model=mod2,statistic=sde.fun1d,R=100, 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].

MCM Based on 100 Batches of 500-Realisations at time 1

    Exact Estimate     Bias Std.Error    RMSE  CI( 2.5 % , 97.5 % )
m 1.00000  1.00288 -0.00288   0.00379 0.03774 ( 0.99545 , 1.01031 )
S 0.75505  0.77451 -0.01946   0.02064 0.20540 ( 0.73406 , 0.81496 )
R> # plot(s) of Monte Carlo outputs of mod1
R> plot(mcm.mod1,index = 1)  # mean
R> plot(mcm.mod1,index = 2)  # variance
 MC output of mean and variance of `mod1` MC output of mean and variance of `mod1`

MC output of mean and variance of mod1

R> # plot(s) of Monte Carlo outputs of mod2
R> plot(mcm.mod2,index = 1)  # mean
R> plot(mcm.mod2,index = 2)  # variance
 MC output of mean and variance of `mod2` MC output of mean and variance of `mod2`

MC output of mean and variance of mod2

Two-dimensional SDEs

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=100,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].

MCM Based on 100 Batches of 500-Realisations at time 10

       Exact Estimate    Bias Std.Error    RMSE    CI( 2.5 % , 97.5 % )
m1   1.99991  1.99579 0.00412   0.00179 0.01785    ( 1.99228 , 1.9993 )
m2  18.00009 17.99768 0.00241   0.00818 0.08142 ( 17.98165 , 18.01371 )
S1   0.25000  0.24999 0.00001   0.00159 0.01585   ( 0.24687 , 0.25311 )
S2   4.25005  4.22234 0.02771   0.02676 0.26623   ( 4.16989 , 4.27479 )
C12  0.24998  0.24573 0.00425   0.00490 0.04876   ( 0.23613 , 0.25533 )

Three-dimensional SDEs

R> mu=0.5;sigma=0.25
R> fx <- expression(mu*y,0,0) 
R> gx <- expression(sigma*z,1,1)
R> modtra <- snssde3d(drift=fx,diffusion=gx,M=500,type="str")
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),var(d$x),cov(d$x,d$y),cov(d$x,d$z)))
+ }
R> ## Monte-Carlo at time = 10
R> mcm.mod3d = MCM.sde(modtra,statistic=sde.fun3d,R=100,parallel="snow",ncpus=2)
R> mcm.mod3d
Stratonovich Sde 3D:
 | dX(t) = mu * Y(t) * dt + sigma * Z(t) o dW1(t)
 | dY(t) = 0 * dt + 1 o dW2(t)
 | dZ(t) = 0 * dt + 1 o dW3(t)
 | t in [0,1].

MCM Based on 100 Batches of 500-Realisations at time 1

    Estimate Std.Error    RMSE   CI( 2.5 % , 97.5 % )
mu1  0.00265   0.00151 0.01503 ( -0.00031 , 0.00561 )
mu2  0.00276   0.00195 0.01943 ( -0.00106 , 0.00658 )
mu3  0.00601   0.00599 0.05962 ( -0.00573 , 0.01775 )
mu4  0.11472   0.00080 0.00796  ( 0.11315 , 0.11629 )
mu5  0.24996   0.00211 0.02099   ( 0.24582 , 0.2541 )
mu6 -0.00317   0.00186 0.01846 ( -0.00682 , 0.00048 )

The MEM.sde() function

R> MEM.sde(drift, diffusion, 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
 | dS(t) = 0.5625 * S(t) + 0.5625 * m(t)^2

Approximation of moment at time 1
 | m(1) = 1
 | S(1) = 0.75505
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> 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,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) * dW1(t)
 | dY(t) = 0 * dt + 1 * dW2(t)
 | dZ(t) = 0 * dt + 1 * dW3(t)
 | t in [0,1].

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

Approximation of moment at time 1                                                      
   | m1(1)  =  5 | S1(1)  =  0.11458 | C12(1)  =  0.25
   | m2(1)  =  0 | S2(1)  =  1.00000 | C13(1)  =  0.25
   | m3(1)  =  0 | S3(1)  =  1.00000 | C23(1)  =  1.00
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. fitsde() - Parametric Estimation of 1-D Stochastic Differential Equation.

References

  1. Guidoum AC, Boukhetala K (2017). Performing Parallel Monte Carlo and Moment Equations Methods for Itô and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc. Preprint submitted to Journal of Statistical Software.

  2. Guidoum AC, Boukhetala K (2017). Sim.DiffProc: Simulation of Diffusion Processes. R package version 3.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 (acguidoum@usthb.dz)

  2. Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, E-mail (kboukhetala@usthb.dz)