version 4.4

Author: ac-guidoum

demo/Sim.DiffProc.R

@@ -1,4 +1,4 @@
-## Wed Oct 17 12:16:40 2018
+## Fri Sep 15 23:08:50 2017
 options(prompt="R> ",scipen=20,digits=5,scientific = 5,width = 70,warning=FALSE, message=FALSE)
 ############################################################################
 #                               Demo 1                                     # 

tests/bridgesde.R

@@ -0,0 +1,21 @@
+library(Sim.DiffProc)
+
+
+### 2-dim Ito bridge sde
+set.seed(1234)
+
+fx <- expression(4*(-1-x) , x)
+gx <- expression(0.2 , 0)
+res <- bridgesde2d(drift=fx,diffusion=gx,Dt=0.005,M=2000)
+res
+summary(res) ## Monte-Carlo statistics at time T/2=2.5
+summary(res,at=1) ## Monte-Carlo statistics at time 1
+summary(res,at=4) ## Monte-Carlo statistics at time 4
+##
+plot(res,type="n")
+lines(time(res),apply(res$X,1,mean),col=3,lwd=2)
+lines(time(res),apply(res$Y,1,mean),col=4,lwd=2)
+legend("topright",c(expression(E(X[t])),expression(E(Y[t]))),lty=1,inset = .7,col=c(3,4))
+##
+plot2d(res)
+

tests/fitsde.R

@@ -0,0 +1,29 @@
+library(Sim.DiffProc)
+
+
+## Application to real data
+## CKLS modele vs CIR modele 
+## CKLS (mod1):  dX(t) = (theta1+theta2* X(t))* dt + theta3 * X(t)^theta4 * dW(t)
+## CIR  (mod2):  dX(t) = (theta1+theta2* X(t))* dt + theta3 * sqrt(X(t))  * dW(t)
+set.seed(1234)
+
+data(Irates)
+rates <- Irates[,"r1"]
+rates <- window(rates, start=1964.471, end=1989.333)
+
+fx1 <- expression(theta[1]+theta[2]*x)
+gx1 <- expression(theta[3]*x^theta[4])
+gx2 <- expression(theta[3]*sqrt(x))
+
+fitmod1 <- fitsde(rates,drift=fx1,diffusion=gx1,pmle="euler",start = list(theta1=1,theta2=1,
+                  theta3=1,theta4=1),optim.method = "L-BFGS-B")
+fitmod2 <- fitsde(rates,drift=fx1,diffusion=gx2,pmle="euler",start = list(theta1=1,theta2=1,
+                  theta3=1),optim.method = "L-BFGS-B")	
+summary(fitmod1)
+summary(fitmod2)
+coef(fitmod1)
+coef(fitmod2)
+confint(fitmod1,parm=c('theta2','theta3'))
+confint(fitmod2,parm=c('theta2','theta3'))
+AIC(fitmod1)
+AIC(fitmod2)

tests/fptsde.R

@@ -0,0 +1,35 @@
+library(Sim.DiffProc)
+
+
+## Example 1: 
+
+# SDE's 2d
+fx <- expression(5*(-1-y)*x , 5*(-1-x)*y)
+gx <- expression(0.5 , 0.5)
+mod2d <- snssde2d(drift=fx,diffusion=gx,x0=c(2,-2),M=2000)
+
+# boundary
+
+St <- expression(-1+5*t)
+
+# random fpt
+
+out <- fptsde2d(mod2d,boundary=St)
+out
+summary(out)
+
+# Marginal density 
+
+denM <- dfptsde2d(out,pdf="M")
+denM
+plot(denM)
+
+# Joint density
+
+denJ <- dfptsde2d(out,pdf="J",n=200,lims=c(0.28,0.4,0.04,0.13))
+denJ
+plot(denJ)
+plot(denJ,display="image")
+plot(denJ,display="image",drawpoints=TRUE,cex=0.5,pch=19,col.pt='green')
+plot(denJ,display="contour")
+plot(denJ,display="contour",color.palette=colorRampPalette(c('white','green','blue','red')))

tests/rsde.R

@@ -0,0 +1,25 @@
+library(Sim.DiffProc)
+
+## 2-dim SDE
+set.seed(1234)
+
+# SDE's 2d
+fx <- expression(3*(2-y),2*x)
+gx <- expression(1,y)
+mod2d <- snssde2d(drift=fx,diffusion=gx,x0=c(1,2),M=5000)
+
+# random 
+r2d <- rsde2d(mod2d,at=0.5)
+summary(r2d)
+
+# Marginal density 
+
+denM <- dsde2d(mod2d,pdf="M", at=0.5)
+denM
+plot(denM)
+
+# Joint density
+denJ <- dsde2d(mod2d,pdf="J", at= 0.5)
+denJ
+plot(denJ)
+plot(denJ,display="contour")

tests/snssde.R

@@ -0,0 +1,18 @@
+library(Sim.DiffProc)
+
+## sde 2-dim
+set.seed(1234)
+
+fx <- expression( y , (4*( 1-x^2 )* y - x))
+gx <- expression( 0 , 0.2)
+
+mod2d2 <- snssde2d(drift=fx,diffusion=gx,type="str",T=100,N=10000)
+mod2d2
+plot(mod2d2,pos=2)
+dev.new()
+plot(mod2d2,union = FALSE)
+dev.new()
+plot2d(mod2d2,type="n") ## in plane (O,X,Y)
+points2d(mod2d2,col=rgb(0,100,0,50,maxColorValue=255), pch=16)
+
+