MCM.sde: Parallel Monte-Carlo Methods for SDE's
In Sim.DiffProc: Simulation of Diffusion Processes

Description Usage Arguments Details Value Note Author(s) References See Also Examples

Generate R Monte-Carlo (version parallel) replicates of a statistic applied to SDE's (1,2 and 3 dim) for the two cases Ito and Stratonovich interpretations.

MCM.sde(model, ...)

## Default S3 method:
MCM.sde(model, statistic, R = 100, time, exact = NULL, 
         names = NULL, level = 0.95, parallel = c("no", "multicore", "snow"),
         ncpus = getOption("ncpus", 1L), cl = NULL, ...)

## S3 method for class 'MCM.sde'
plot(x,index = 1,type=c("all","hist","qqplot","boxplot","CI"), ...)

`model`	an object from class `snssde1d`, `snssde2d` and `snssde3d`.
`statistic`	a function which when applied to model returns a vector containing the statistic(s) of interest.
`R`	the number of Monte-Carlo replicates. Usually this will be a single positive integer "R > 1".
`time`	the time when estimating the statistic(s) of interesttime between `t0` and `T`. The default `time = T`.
`exact`	a named list giving the exact statistic(s) if it exists otherwise `exact = NULL`.
`names`	named the statistic(s) of interest. The default `names=c("mu1","mu2",...)`.
`level`	the confidence level(s) of the required interval(s).
`parallel`	the type of parallel operation to be used (if any). The default `parallel = "no"`.
`ncpus`	integer: number of processes to be used in parallel operation: typically one would chose this to the number of available CPUs.
`cl`	an optional parallel or snow cluster for use if `parallel = "snow"`.
`x`	an object inheriting from class `"MCM.sde"`.
`index`	the index of the variable of interest within the output of `"MCM.sde"`.
`type`	the type of plot of the Monte-Carlo estimation of the variable of interest. The default `type = "all"`.
`...`	potentially further arguments for (non-default) methods.

We have here developed Monte-Carlo methods whose essence is the use of repeated experiments to evaluate a statistic(s) of interest in SDE's. For example estimation of moments as: mean, variance, covariance (and other as median, mode, quantile,...). With the standard error and the confidence interval for these estimators.

An overview of this package, see browseVignettes('Sim.DiffProc') for more informations.

The returned value is an object of class "MCM.sde", containing the following components:

`mod`	The SDE's used (class: `snssde1d`, `snssde2d` and `snssde3d`).
`dim`	Dimension of the model.
`call`	The original call to `"MCM.sde"`.
`Fn`	The function statistic as passed to `"MCM.sde"`.
`ech`	A matrix with `sum(R)` column each of which is a Monte-Carlo replicate of the result of calling statistic.
`time`	The time when estimating the statistic(s) of interest.
`name`	named of statistic(s) of interest.
`MC`	Table contains simulation results of statistic(s) of interest: Estimate, Bias (if exact available), Std.Error and Confidence interval.

When parallel = "multicore" is used are not available on Windows, parallel = "snow" is primarily intended to be used on multi-core Windows machine where parallel = "multicore" is not available. For more details see Q.E.McCallum and S.Weston (2011).

A.C. Guidoum, K. Boukhetala.

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. doi:10.18637/jss.v096.i02

Paul Glasserman (2003). Monte Carlo Methods in Financial Engineering. Springer-Verlag New York.

Jun S. Liu (2004). Monte Carlo Strategies in Scientific Computing. Springer-Verlag New York.

Christian Robert and George Casella (2010). Introducing Monte Carlo Methods with R. Springer-Verlag New York.

Nick T. Thomopoulos (2013). Essentials of Monte Carlo Simulation: Statistical Methods for Building Simulation Models. Springer-Verlag New York.

Q. Ethan McCallum and Stephen Weston (2011). Parallel R. O'Reilly Media, Inc.

MEM.sde moment equations methods for SDE's.

## Example 1 : (1 dim)
## dX(t) = 3*(1-X(t)) dt + 0.5 * dW(t), X(0)=5, t in [0,10]
## set the model 1d
f <- expression(3*(1-x));g <- expression(0.5)
mod1d <-  snssde1d(drift=f,diffusion=g,x0=5,T=10,M=50)

## function of the statistic(s) of interest.
sde.fun1d <- function(data, i){
  d <- data[i, ]
  return(c(mean(d),Mode(d),var(d)))
}

mc.sde1d = MCM.sde(model=mod1d,statistic=sde.fun1d,R=100,exact=list(Me=1,Mo=1,Va=0.5^2/6),
                 names=c("Me(10)","Mo(10)","Va(10)"))
mc.sde1d
plot(mc.sde1d,index=1)
plot(mc.sde1d,index=2)
plot(mc.sde1d,index=3)

## Example 2 : with Parallel computing
## Not run: 
mod1d <-  snssde1d(drift=f,diffusion=g,x0=5,T=10,M=1000)
## On Windows or Unix
mc.sde1d = MCM.sde(model=mod1d,statistic=sde.fun1d,R=1000,exact=list(Me=1,Mo=1,Va=0.5^2/6),
           names=c("Me(10)","Mo(10)","Va(10)"),parallel="snow",ncpus=parallel::detectCores())
mc.sde1d
## On Unix only
mc.sde1d = MCM.sde(model=mod1d,statistic=sde.fun1d,R=1000,exact=list(Me=1,Mo=1,Va=0.5^2/6),
           names=c("Me(10)","Mo(10)","Va(10)"),parallel="multicore",ncpus=parallel::detectCores())
mc.sde1d

## End(Not run)

## Example 3: (2 dim)
## dX(t) = 1/mu*(theta-X(t)) dt + sqrt(sigma) * dW1(t),
## dY(t) = X(t) dt + 0 * dW2(t)
## Not run: 
## Set the model 2d
mu=0.75;sigma=0.1;theta=2
x0=0;y0=0;init=c(x=0,y=0)
f <- expression(1/mu*(theta-x), x)  
g <- expression(sqrt(sigma),0)
OUI <- snssde2d(drift=f,diffusion=g,M=1000,Dt=0.01,x0=init)

## function of the statistic(s) of interest.
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)))
}
## Monte-Carlo at time = 5
mc.sde2d_a = MCM.sde(model=OUI,statistic=sde.fun2d,R=100,time=5,
                  parallel="snow",ncpus=parallel::detectCores())
mc.sde2d_a
## Monte-Carlo at time = 10
mc.sde2d_b = MCM.sde(model=OUI,statistic=sde.fun2d,R=100,time=10,
                  parallel="snow",ncpus=parallel::detectCores())
mc.sde2d_b

## Compared with exact values at time 5 and 10
E_x <- function(t) theta+(x0-theta)*exp(-t/mu)
V_x <- function(t) 0.5*sigma*mu *(1-exp(-2*(t/mu)))
E_y <- function(t) y0+theta*t+(x0-theta)*mu*(1-exp(-t/mu))
V_y <- function(t) sigma*mu^3*((t/mu)-2*(1-exp(-t/mu))+0.5*(1-exp(-2*(t/mu))))
cov_xy <- function(t) 0.5*sigma*mu^2 *(1-2*exp(-t/mu)+exp(-2*(t/mu)))

## at time=5
mc.sde2d_a = MCM.sde(model=OUI,statistic=sde.fun2d,R=100,time=5,
    exact=list(m1=E_x(5),m2=E_y(5),S1=V_x(5),S2=V_y(5),C12=cov_xy(5)),
	parallel="snow",ncpus=parallel::detectCores())
mc.sde2d_a
plot(mc.sde2d_a,index=1)
plot(mc.sde2d_a,index=2)
## at time=10
mc.sde2d_b = MCM.sde(model=OUI,statistic=sde.fun2d,R=100,time=10,
    exact=list(m1=E_x(10),m2=E_y(10),S1=V_x(10),S2=V_y(10),C12=cov_xy(10)),
	parallel="snow",ncpus=parallel::detectCores())
mc.sde2d_b
plot(mc.sde2d_b,index=1)
plot(mc.sde2d_b,index=2)

## End(Not run)

## Example 4: (3 dim)
## dX(t) = sigma*(Y(t)-X(t)) dt + 0.1 * dW1(t)
## dY(t) = (rho*X(t)-Y(t)-X(t)*Z(t)) dt + 0.1 * dW2(t)
## dZ(t) = (X(t)*Y(t)-bet*Z(t)) dt + 0.1 * dW3(t)
## W1(t), W2(t) and W3(t) are three correlated Brownian motions with Sigma
## Not run: 
## Set the model 3d
sigma=10;rho=28; bet=8/3
f <- expression(sigma*(y-x),rho*x-y-x*z,x*y-bet*z)
g <- expression(0.1,0.1,0.1)
# correlation matrix
Sigma <-matrix(c(1,0.3,0.5,0.3,1,0.2,0.5,0.2,1),nrow=3,ncol=3) 
mod3d <- snssde3d(x0=rep(0,3),drift=f,diffusion=g,M=1000,Dt=0.01,corr=Sigma)

## function of the statistic(s) of interest.
sde.fun3d <- function(data, i){
  d <- data[i,]
  return(c(mean(d$x),mean(d$y),mean(d$z)))
}
## Monte-Carlo at time = 10
mc.sde3d = MCM.sde(mod3d,statistic=sde.fun3d,R=100,parallel="snow",ncpus=parallel::detectCores())
mc.sde3d

## End(Not run)