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:

model: an object from classes snssde1d(), snssde2d() and snssde3d().
statistic: a function which when applied to the model (SDEs) returns a vector containing the statistic(s) of interest. Any further arguments can be passed to statistic(s) through the ... argument.
R: number of Monte Carlo replicates (R batches), this will be a single positive integer.
time: fixed time at which the estimate of the statistic(s).
exact: a named list giving the exact statistic(s), if it exists the bias calculation will be performed.
names: named the statistic(s) of interest. Default names=c("mu1","mu2",...).
level: confidence level of the required interval(s).
parallel: the type of parallel operation to be used. "multicore" does not work on Microsoft Windows operating systems, but on Unix is allowed and uses parallel operations. Default parallel="no".
ncpus: an integer value specifying the number of cores to be used in the parallelized procedure. Default is 1 core of the machine.
cl: an optional parallel cluster for use if parallel = "snow". Default cl = makePSOCKcluster(rep("localhost", ncpus)).

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.

x: an object from the class MCM.sde().
index: the index of the variable of interest within the output of class MCM.sde().
type: type of plots. Default type="all".

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:

drift: an R vector of expressions which contains the drift specification (1D, 2D and 3D).
diffusion: an R vector of expressions which contains the diffusion specification (1D, 2D and 3D).
corr: the correlation coefficient ‘|corr|<=1’ of \(W_{1}(t)\) and \(W_{2}(t)\) (2D) must be an expression length equal 1. And for 3D \((W_{1}(t),W_{2}(t),W_{3}(t))\) an expressions length equal 3.
type: type of SDEs to be used "ito" for Ito form and "str" for Stratonovich form. The default type="ito".
solve: if solve=FALSE only the symbolic computational of system will be made. And if solve=TRUE a numerical approximation of the obtained system will be performed.
parms: parameters passed to drift and diffusion.
init: initial (state) values of system.
time: time sequence (vector) for which output is sought; the first value of time must be the initial time.
...: arguments to be passed to ode() function available in deSolve package, if solve=TRUE.

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

Department of Mathematics and Computer Science, Faculty of Sciences and Technology, University of Tamanghasset, Algeria, E-mail (acguidoum@univ-tam.dz)↩︎
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)↩︎

Parallel Monte-Carlo and Moment Equations for SDEs

A.C. Guidoum¹ and K. Boukhetala²

2024-03-05

The `MCM.sde()` function

One-dimensional SDE

Two-dimensional SDEs

Three-dimensional SDEs

The `MEM.sde()` function

One-dimensional SDE

Two-dimensional SDEs

Three-dimensional SDEs

Further reading

References

Parallel Monte-Carlo and Moment Equations for SDEs

A.C. Guidoum1 and K. Boukhetala2

2024-03-05

The MCM.sde() function

One-dimensional SDE

Two-dimensional SDEs

Three-dimensional SDEs

The MEM.sde() function

One-dimensional SDE

Two-dimensional SDEs

Three-dimensional SDEs

Further reading

References

A.C. Guidoum¹ and K. Boukhetala²

The `MCM.sde()` function

The `MEM.sde()` function