In acguidoum/Sim.DiffProc: Simulation of Diffusion Processes

library(Sim.DiffProc)
library(knitr)
library(deSolve)
knitr::opts_chunk$set(comment="", prompt=TRUE, fig.show='hold',warning=FALSE, message=FALSE)
options(prompt="R> ",scipen=16,digits=5,warning=FALSE, message=FALSE,mc.cores=2)

The `MCM.sde()` function

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

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

set.seed(1234, kind = "L'Ecuyer-CMRG")
theta = 0.75; x0 = 1
fx <- expression( 0.5*theta^2*x )
gx <- expression( theta*x )
mod1 <- snssde1d(drift=fx,diffusion=gx,x0=x0,M=500,type="ito")
mod2 <- snssde1d(drift=fx,diffusion=gx,x0=x0,M=500,type="str")
## True values of means and variance for mod1 and mod2
E.mod1 <- function(t) x0 * exp(0.5 * theta^2 * t)
V.mod1 <- function(t) x0^2 * exp(theta^2 * t) * (exp(theta^2 * t) - 1)
E.mod2 <- function(t) x0 * exp(theta^2 * t)
V.mod2 <- function(t) x0^2 * exp(2 * theta^2 * t) * (exp(theta^2 * t) - 1)
## function of the statistic(s) of interest.
sde.fun1d <- function(data, i){
     d <- data[i, ]
     return(c(mean(d),var(d)))
}
# Parallel MOnte Carlo for mod1
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)
mcm.mod1
# Parallel MOnte Carlo for mod2
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)
mcm.mod2

Two-dimensional SDEs

set.seed(1234, kind = "L'Ecuyer-CMRG")
mu=1;sigma=0.5;theta=2
x0=0;y0=0;init=c(x0,y0)
f <- expression(1/mu*(theta-x), x)  
g <- expression(sqrt(sigma),0)
OUI <- snssde2d(drift=f,diffusion=g,M=500,Dt=0.015,x0=c(x=0,y=0))
## true values of first and second moment at time 10
Ex <- function(t) theta+(x0-theta)*exp(-t/mu)
Vx <- function(t) 0.5*sigma*mu *(1-exp(-2*(t/mu)))
Ey <- function(t) y0+theta*t+(x0-theta)*mu*(1-exp(-t/mu))
Vy <- function(t) sigma*mu^3*((t/mu)-2*(1-exp(-t/mu))+0.5*(1-exp(-2*(t/mu))))
covxy <- function(t) 0.5*sigma*mu^2 *(1-2*exp(-t/mu)+exp(-2*(t/mu)))
tvalue = list(m1=Ex(10),m2=Ey(10),S1=Vx(10),S2=Vy(10),C12=covxy(10))
## 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)))
}
## Parallel Monte-Carlo of 'OUI' at time 10
mcm.mod2d = MCM.sde(OUI,statistic=sde.fun2d,time=10,R=20,exact=tvalue,parallel="snow",ncpus=2)
mcm.mod2d

Three-dimensional SDEs

set.seed(1234, kind = "L'Ecuyer-CMRG")
mu=0.5;sigma=0.25
fx <- expression(mu*y,0,0) 
gx <- expression(sigma*z,1,1)
Sigma <-matrix(c(1,0.3,-0.5,0.3,1,0.2,-0.5,0.2,1),nrow=3,ncol=3)
modtra <- snssde3d(drift=fx,diffusion=gx,M=500,type="str",corr=Sigma)
## function of the statistic(s) of interest.
sde.fun3d <- function(data, i){
  d <- data[i,]
  return(c(mean(d$x),median(d$x),Mode(d$x)))
}
## Monte-Carlo at time = 10
mcm.mod3d = MCM.sde(modtra,statistic=sde.fun3d,R=10,parallel="snow",ncpus=2)
mcm.mod3d

The `MEM.sde()` function

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

fx <- expression( 0.5*theta^2*x )
gx <- expression( theta*x )
start = c(m=1,S=0)
t = seq(0,1,by=0.001)
mem.mod1 = MEM.sde(drift=fx,diffusion=gx,type="ito",solve = TRUE,parms = c(theta=0.75), init = start, time = t)
mem.mod1
mem.mod2 = MEM.sde(drift=fx,diffusion=gx,type="str",solve = TRUE,parms = c(theta=0.75), init = start, time = t)
mem.mod2

plot(mem.mod1$sol.ode, mem.mod2$sol.ode,ylab = c("m(t)"),select="m", xlab = "Time",main="",col = 2:3,lty=1)
legend("topleft",c(expression(m[mod1](t),m[mod2](t))),inset = .05, col=2:3,lty=1)
plot(mem.mod1$sol.ode, mem.mod2$sol.ode,ylab = c("S(t)"),select="S", xlab = "Time",main="",col = 2:3,lty=1)
legend("topleft",c(expression(S[mod1](t),S[mod2](t))),inset = .05, col=2:3,lty=1)

Two-dimensional SDEs

fx <- expression(1/mu*(theta-x), x)  
gx <- expression(sqrt(sigma),0)
start = c(m1=0,m2=0,S1=0,S2=0,C12=0)
t = seq(0,10,by=0.001)
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)
mem.mod2d

matplot.0D(mem.mod2d$sol.ode,main="")

Three-dimensional SDEs

fx <- expression(mu*y,0,0) 
gx <- expression(sigma*z,1,1)
RHO <- expression(0.75,0.5,-0.25)
start = c(m1=5,m2=0,m3=0,S1=0,S2=0,S3=0,C12=0,C13=0,C23=0)
t = seq(0,1,by=0.001)
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)
mem.mod3d

matplot.0D(mem.mod3d$sol.ode,main="",select=c("m1","m2","m3"))
matplot.0D(mem.mod3d$sol.ode,main="",select=c("S1","S2","S3"))
matplot.0D(mem.mod3d$sol.ode,main="",select=c("C12","C13","C23"))

References

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

acguidoum/Sim.DiffProc documentation built on March 8, 2024, 8:50 p.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

Tweet to @rdrrHQ

GitHub issue tracker

ian@mutexlabs.com