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=10, 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 10 Batches of 500-Realisations at time 1
Exact Estimate Bias Std.Error RMSE CI( 2.5 % , 97.5 % )
m 1.3248 1.3275 -0.00272 0.01621 0.04870 ( 1.29573 , 1.35927 )
S 1.3252 1.4101 -0.08497 0.06600 0.21546 ( 1.28077 , 1.53949 )
R> # Parallel MOnte Carlo for mod2
R> mcm.mod2 = MCM.sde(model=mod2,statistic=sde.fun1d,R=10, 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 10 Batches of 500-Realisations at time 1
Exact Estimate Bias Std.Error RMSE CI( 2.5 % , 97.5 % )
m 1.00000 0.99182 0.00818 0.01126 0.03477 ( 0.96975 , 1.01389 )
S 0.75505 0.69857 0.05648 0.03657 0.12339 ( 0.62689 , 0.77025 )
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=10,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 10 Batches of 500-Realisations at time 10
Exact Estimate Bias Std.Error RMSE CI( 2.5 % , 97.5 % )
m1 1.99991 1.99643 0.00348 0.00798 0.02420 ( 1.98079 , 2.01207 )
m2 18.00009 17.99736 0.00273 0.03196 0.09592 ( 17.93472 , 18.06 )
S1 0.25000 0.24152 0.00848 0.00342 0.01332 ( 0.23482 , 0.24822 )
S2 4.25005 4.33229 -0.08224 0.09641 0.30070 ( 4.14333 , 4.52125 )
C12 0.24998 0.22424 0.02574 0.01646 0.05568 ( 0.19198 , 0.2565 )
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=10,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 10 Batches of 500-Realisations at time 1
Estimate Std.Error CI( 2.5 % , 97.5 % )
mu1 -0.00283 0.00396 ( -0.01059 , 0.00493 )
mu2 0.00279 0.00560 ( -0.00819 , 0.01377 )
mu3 0.02416 0.01288 ( -0.00108 , 0.0494 )
mu4 0.11522 0.00233 ( 0.11065 , 0.11979 )
mu5 0.24563 0.00514 ( 0.23556 , 0.2557 )
mu6 0.00402 0.00489 ( -0.00556 , 0.0136 )
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:
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).type: type of SDEs to be used "ito" for It? 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
| 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"))
