snssde.R
In Sim.DiffProc: Simulation of Diffusion Processes

## ----setup, echo = F, message = F, results = 'hide',screenshot.force=FALSE----
library(Sim.DiffProc)
library(knitr)
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,
        width = 70)

## -------------------------------------------------------------------
set.seed(1234, kind = "L'Ecuyer-CMRG")
theta = 0.5
f <- expression( (0.5*theta^2*x) )
g <- expression( theta*x )
mod1 <- snssde1d(drift=f,diffusion=g,x0=10,M=1000,type="ito") # Using Ito
mod2 <- snssde1d(drift=f,diffusion=g,x0=10,M=1000,type="str") # Using Stratonovich 
mod1
mod2

## -------------------------------------------------------------------
summary(mod1, at = 1)
summary(mod2, at = 1)

## -------------------------------------------------------------------
x1 <- rsde1d(object = mod1, at = 1)  # X(t=1) | X(0)=x0 (Ito SDE)
x2 <- rsde1d(object = mod2, at = 1)  # X(t=1) | X(0)=x0 (Stratonovich SDE)
head(data.frame(x1,x2),n=5)

## ----01,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE-----
#  mu1 = log(10); sigma1= sqrt(theta^2)  # log mean and log variance for mod1
#  mu2 = log(10)-0.5*theta^2 ; sigma2 = sqrt(theta^2) # log mean and log variance for mod2
#  AppdensI <- dsde1d(mod1, at = 1)
#  AppdensS <- dsde1d(mod2, at = 1)
#  plot(AppdensI , dens = function(x) dlnorm(x,meanlog=mu1,sdlog = sigma1))
#  plot(AppdensS , dens = function(x) dlnorm(x,meanlog=mu2,sdlog = sigma2))

## ----001, echo=FALSE, fig.cap='Approximate transitional density for $X_{t}|X_{0}$ at time $t-s=1$ with log-normal curves. mod1: Ito and mod2: Stratonovich  ', fig.env='figure*',fig.width=10,fig.height=10----
knitr::include_graphics(c("Figures/fig007.png","Figures/fig008.png"))

## ----02,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE-----
#  ## Ito
#  plot(mod1,ylab=expression(X^mod1))
#  lines(time(mod1),apply(mod1$X,1,mean),col=2,lwd=2)
#  lines(time(mod1),apply(mod1$X,1,bconfint,level=0.95)[1,],col=4,lwd=2)
#  lines(time(mod1),apply(mod1$X,1,bconfint,level=0.95)[2,],col=4,lwd=2)
#  legend("topleft",c("mean path",paste("bound of", 95,"% confidence")),inset = .01,col=c(2,4),lwd=2,cex=0.8)
#  ## Stratonovich
#  plot(mod2,ylab=expression(X^mod2))
#  lines(time(mod2),apply(mod2$X,1,mean),col=2,lwd=2)
#  lines(time(mod2),apply(mod2$X,1,bconfint,level=0.95)[1,],col=4,lwd=2)
#  lines(time(mod2),apply(mod2$X,1,bconfint,level=0.95)[2,],col=4,lwd=2)
#  legend("topleft",c("mean path",paste("bound of",95,"% confidence")),col=c(2,4),inset =.01,lwd=2,cex=0.8)

## ----100, echo=FALSE, fig.cap='mod1: Ito and mod2: Stratonovich  ', fig.env='figure*',fig.width=7,fig.height=7----
knitr::include_graphics(c("Figures/fig07.png","Figures/fig08.png"))

## -------------------------------------------------------------------
set.seed(1234, kind = "L'Ecuyer-CMRG")
x0=5;y0=0
mu=3;sigma=0.5
fx <- expression(-(x/mu),x)  
gx <- expression(sqrt(sigma),0)
mod2d <- snssde2d(drift=fx,diffusion=gx,Dt=0.01,M=1000,x0=c(x0,y0),method="smilstein")
mod2d

## ----eval=FALSE, include=TRUE---------------------------------------
#  summary(mod2d, at = 10)

## ----03,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE-----
#  ## in time
#  plot(mod2d)
#  ## in plane (O,X,Y)
#  plot2d(mod2d,type="n")
#  points2d(mod2d,col=rgb(0,100,0,50,maxColorValue=255), pch=16)

## ----102, echo=FALSE, fig.cap=' Ornstein-Uhlenbeck process and its integral ', fig.env='figure*',fig.width=10,fig.height=10----
knitr::include_graphics(c("Figures/fig09.png","Figures/fig009.png"))

## -------------------------------------------------------------------
out <- rsde2d(object = mod2d, at = 10) 
head(out,n=3)

## ----04,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE-----
#  ## the marginal density
#  denM <- dsde2d(mod2d,pdf="M",at =10)
#  plot(denM, main="Marginal Density")
#  ## the Joint density
#  denJ <- dsde2d(mod2d, pdf="J", n=100,at =10)
#  plot(denJ,display="contour",main="Bivariate Transition Density at time t=10")

## ----1002, echo=FALSE, fig.cap='Marginal and Joint density at time t=10 ', fig.env='figure*',fig.width=10,fig.height=10----
knitr::include_graphics(c("Figures/fig1001.png","Figures/fig1002.png"))

## ----07,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE-----
#  plot(denJ,display="persp",main="Bivariate Transition Density at time t=10")

## ----10002, echo=FALSE, fig.cap='Marginal and Joint density at time t=10 ', fig.env='figure*',fig.width=10,fig.height=10----
knitr::include_graphics(c("Figures/fig1003.png"))

## ----eval=FALSE, include=TRUE---------------------------------------
#  for (i in seq(1,10,by=0.005)){
#  plot(dsde2d(mod2d, at = i,n=100),display="contour",main=paste0('Transition Density \n t = ',i))
#  }

## -------------------------------------------------------------------
set.seed(1234, kind = "L'Ecuyer-CMRG")
mu = 4; sigma=0.1
fx <- expression( y ,  (mu*( 1-x^2 )* y - x)) 
gx <- expression( 0 ,2*sigma)
mod2d <- snssde2d(drift=fx,diffusion=gx,N=10000,Dt=0.01,type="str",method="rk1")

## ----9,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE------
#  plot(mod2d,ylim=c(-8,8))   ## back in time
#  plot2d(mod2d)              ## in plane (O,X,Y)

## ----100002, echo=FALSE, fig.cap='The stochastic Van-der-Pol equation', fig.env='figure*',fig.width=10,fig.height=10----
knitr::include_graphics(c("Figures/fig1004.png","Figures/fig1005.png"))

## -------------------------------------------------------------------
set.seed(1234, kind = "L'Ecuyer-CMRG")
mu = 1.2; sigma=0.1;nu=2;theta=0.5
fx <- expression( mu*x ,nu*(theta-y)) 
gx <- expression( x*sqrt(y) ,sigma*sqrt(y))
Sigma <- matrix(c(1,0.3,0.3,2),nrow=2,ncol=2) # correlation matrix
HM <- snssde2d(drift=fx,diffusion=gx,Dt=0.001,x0=c(100,1),corr=Sigma,M=1000)
HM

## -------------------------------------------------------------------
out <- rsde2d(object = HM, at = 1) 
head(out,n=3)

## ----eval=FALSE, include=TRUE---------------------------------------
#  denJ <- dsde2d(HM,pdf="J",at =1,lims=c(-100,900,0.4,0.75))
#  plot(denJ,display="contour",main="Bivariate Transition Density at time t=10")
#  plot(denJ,display="persp",main="Bivariate Transition Density at time t=10")

## -------------------------------------------------------------------
set.seed(1234, kind = "L'Ecuyer-CMRG")
fx <- expression(4*(-1-x)*y , 4*(1-y)*x , 4*(1-z)*y) 
gx <- rep(expression(0.2),3)
mod3d <- snssde3d(x0=c(x=2,y=-2,z=-2),drift=fx,diffusion=gx,M=1000)
mod3d

## ----eval=FALSE, include=TRUE---------------------------------------
#  s = 1
#  mean(mod3d, at = s)
#  moment(mod3d, at = s , center = TRUE , order = 2) ## variance
#  Median(mod3d, at = s)
#  Mode(mod3d, at = s)
#  quantile(mod3d , at = s)
#  kurtosis(mod3d , at = s)
#  skewness(mod3d , at = s)
#  cv(mod3d , at = s )
#  min(mod3d , at = s)
#  max(mod3d , at = s)
#  moment(mod3d, at = s , center= TRUE , order = 4)
#  moment(mod3d, at = s , center= FALSE , order = 4)

## ----eval=FALSE, include=TRUE---------------------------------------
#  summary(mod3d, at = t)

## ----10,fig.env='figure*', fig.cap=' ',eval=FALSE, include=TRUE-----
#  plot(mod3d,union = TRUE)         ## back in time
#  plot3D(mod3d,display="persp")    ## in space (O,X,Y,Z)

## ----103, echo=FALSE, fig.cap=' Flow of $1000$ trajectories of $(X_t ,Y_t ,Z_t)$ ', fig.env='figure*',fig.width=7,fig.height=7----
knitr::include_graphics(c("Figures/fig10.png","Figures/fig11.png"))

## -------------------------------------------------------------------
out <- rsde3d(object = mod3d, at = 1) 
head(out,n=3)

## ----11,fig.env='figure*', fig.cap=' Marginal density of $X_t$, $Y_t$ and $Z_t$ at time $t=1$ ',fig.width=3.5,fig.height=3.5----
den <- dsde3d(mod3d,pdf="M",at =1)
plot(den, main="Marginal Density") 

## ----111,fig.env='figure*', fig.cap='  ',eval=FALSE, include=TRUE,fig.width=5,fig.height=5----
#  denJ <- dsde3d(mod3d,pdf="J")
#  plot(denJ,display="rgl")

## -------------------------------------------------------------------
set.seed(1234, kind = "L'Ecuyer-CMRG")
K = 4; s = 1; sigma = 0.2
fx <- expression( (-K*x/sqrt(x^2+y^2+z^2)) , (-K*y/sqrt(x^2+y^2+z^2)) , (-K*z/sqrt(x^2+y^2+z^2)) ) 
gx <- rep(expression(sigma),3)
mod3d <- snssde3d(drift=fx,diffusion=gx,N=10000,x0=c(x=1,y=1,z=1))

## ----12,fig.env='figure*',fig.width=3.5,fig.height=3.5, fig.cap=' Attractive model for 3D diffusion processes '----
plot3D(mod3d,display="persp",col="blue")

## -------------------------------------------------------------------
set.seed(1234, kind = "L'Ecuyer-CMRG")
fx <- expression(y,0,0) 
gx <- expression(z,1,1)
Sigma <-matrix(c(1,0.2,0.5,0.2,1,-0.7,0.5,-0.7,1),nrow=3,ncol=3)
modtra <- snssde3d(drift=fx,diffusion=gx,M=1000,corr=Sigma)
modtra

## ----14, fig.cap='  ', fig.env='figure*',eval=FALSE, include=TRUE----
#  X <- rsde3d(modtra,at=1)$x
#  MASS::truehist(X,xlab = expression(X[t==1]));box()
#  lines(density(X),col="red",lwd=2)
#  legend("topleft",c("Distribution histogram","Kernel Density"),inset =.01,pch=c(15,NA),lty=c(NA,1),col=c("cyan","red"), lwd=2,cex=0.8)
#  ## Cov-Matrix
#  color.palette=colorRampPalette(c('white','green','blue','red'))
#  filled.contour(time(modtra), time(modtra), cov(t(modtra$X)), color.palette=color.palette,plot.title = title(main = expression(paste("Covariance empirique:",cov(X[s],X[t]))),xlab = "time", ylab = "time"),key.title = title(main = ""))

## ----1000002, echo=FALSE, fig.cap='The histogram and kernel density of $X_t$ at time $t=1$. Emprical variance-covariance matrix', fig.env='figure*',fig.width=10,fig.height=10----
knitr::include_graphics(c("Figures/fig1007.png","Figures/fig1006.png"))