# Parallel Monte-Carlo and Moment Equations for SDEs

#### A.C. Guidoum1 and K. Boukhetala2

#### 2020-05-04

# 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> 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=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.3227 -0.00205   0.01438 0.06270 ( 1.29455 , 1.35091 )
    S 1.3252   1.3354  0.01024   0.07996 0.34869 ( 1.17868 , 1.49212 )

    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.00000   1.7869 0.78688   0.01597 0.78996 ( 1.75558 , 1.81818 )
    S 0.75505   2.6798 1.92479   0.15066 2.03374 ( 2.38455 , 2.97513 )

## 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=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  1.99647 -0.00344   0.00483 0.02135
    m2  18.00009 18.00058  0.00049   0.02262 0.09858
    S1   0.25000  0.24924 -0.00076   0.00258 0.01129
    S2   4.25005  4.31008  0.06003   0.06261 0.27943
    C12  0.24998  0.24965 -0.00033   0.01129 0.04921
           CI( 2.5 % , 97.5 % )
    m1      ( 1.987 , 2.00594 )
    m2  ( 17.95625 , 18.04491 )
    S1     ( 0.24418 , 0.2543 )
    S2    ( 4.18737 , 4.43279 )
    C12   ( 0.22752 , 0.27178 )

## 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> 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.07146   0.00319 ( -0.07771 , -0.06521 )
    mu2 -0.05744   0.00512  ( -0.06748 , -0.0474 )
    mu3 -0.02513   0.02182   ( -0.0679 , 0.01764 )

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

# 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 Ito and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc. Accept Submission to Journal of Statistical Software.

2. Guidoum AC, Boukhetala K (2020). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.6, 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)↩