The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
MCM.sde() functionR> 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)).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".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.mod1Itô 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.mod2Stratonovich 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 )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.mod2dItô 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 )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.mod3dStratonovich 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 )MEM.sde() functionR> 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.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.mod1Itô 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.3252R> mem.mod2 = MEM.sde(drift=fx,diffusion=gx,type="str",solve = TRUE,parms = c(theta=0.75), init = start, time = t)
R> mem.mod2Stratonovich 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.3257R> 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)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.mod2dItô 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> 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.mod3dItô 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.2500snssdekd() &
dsdekd() & rsdekd()- Monte-Carlo
Simulation and Analysis of Stochastic Differential Equations.bridgesdekd() &
dsdekd() & rsdekd() - Constructs and
Analysis of Bridges Stochastic Differential Equations.fptsdekd() &
dfptsdekd() - Monte-Carlo Simulation and Kernel Density
Estimation of First passage time.MCM.sde() &
MEM.sde() - Parallel Monte-Carlo and Moment Equations for
SDEs.TEX.sde() - Converting
Sim.DiffProc Objects to LaTeX.fitsde() - Parametric Estimation
of 1-D Stochastic Differential Equation.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)↩︎
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.