R Under development (unstable) (2024-03-03 r86036 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Sim.DiffProc)
Package 'Sim.DiffProc', version 4.9
browseVignettes('Sim.DiffProc') for more informations.
> 
> ## 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
Stratonovich Sde 2D:
	| dX(t) = Y(t) * dt + 0 o dW1(t)
	| dY(t) = (4 * (1 - X(t)^2) * Y(t) - X(t)) * dt + 0.2 o dW2(t)
Method:
	| Euler scheme with order 0.5
Summary:
	| Size of process	| N  = 10001.
	| Number of simulation	| M  = 1.
	| Initial values	| (x0,y0) = (0,0).
	| Time of process	| t in [0,100].
	| Discretization	| Dt = 0.01.
> plot(mod2d2,pos=2)
> dev.new()
dev.new(): using pdf(file="Rplots1.pdf")
> plot(mod2d2,union = FALSE)
dev.new(): using pdf(file="Rplots2.pdf")
> dev.new()
dev.new(): using pdf(file="Rplots3.pdf")
> plot2d(mod2d2,type="n") ## in plane (O,X,Y)
> points2d(mod2d2,col=rgb(0,100,0,50,maxColorValue=255), pch=16)
> 
> 
> ## dX(t) = 4*(-1-X(t))*Y(t) dt + 0.2 dW1(t)
> ## dY(t) = 4*(1-Y(t))*X(t) dt + 0.2 dW2(t)
> set.seed(1234)
> 
> fx <- expression(4*(-1-x)*y , 4*(1-y)*x )
> gx <- expression(0.25*y,0.2*x)
> 
> Sigma <- matrix(c(0.9,0.3,0.3,0.4),nrow=2,ncol=2) # correlation matrix
> mod2d1 <- snssde2d(drift=fx,diffusion=gx,corr=Sigma,x0=c(x0=1,y0=-1),M=500)
> mod2d1
Itô Sde 2D:
	| dX(t) = 4 * (-1 - X(t)) * Y(t) * dt + 0.25 * Y(t) * dB1(t)
	| dY(t) = 4 * (1 - Y(t)) * X(t) * dt + 0.2 * X(t) * dB2(t)
	| Correlation structure:                    
             0.9 0.3
             0.3 0.4
Method:
	| Euler scheme with order 0.5
Summary:
	| Size of process	| N  = 1001.
	| Number of simulation	| M  = 500.
	| Initial values	| (x0,y0) = (1,-1).
	| Time of process	| t in [0,1].
	| Discretization	| Dt = 0.001.
> summary(mod2d1)

	Monte-Carlo Statistics for (X(t),Y(t)) at time t = 1
                          X         Y
Mean             -0.6926907 0.4460955
Variance          0.0064178 0.0203865
Median           -0.6902043 0.4470989
Mode             -0.6865062 0.4770332
First quartile   -0.7446477 0.3533673
Third quartile   -0.6367109 0.5307240
Minimum          -1.0023128 0.0465707
Maximum          -0.4513685 1.2314521
Skewness         -0.2418736 0.5798726
Kurtosis          3.4039294 5.8010997
Coef-variation   -0.1156520 0.3200692
3th-order moment -0.0001244 0.0016879
4th-order moment  0.0001402 0.0024110
5th-order moment -0.0000090 0.0011300
6th-order moment  0.0000058 0.0009784
> ##
> dev.new()
dev.new(): using pdf(file="Rplots4.pdf")
> plot(mod2d1,type="n")
> mx <- apply(mod2d1$X,1,mean)
> my <- apply(mod2d1$Y,1,mean)
> lines(time(mod2d1),mx,col=1)
> lines(time(mod2d1),my,col=2)
> legend("topright",c(expression(E(X[t])),expression(E(Y[t]))),lty=1,inset = .01,col=c(1,2),cex=0.95)
> ##
> dev.new()
dev.new(): using pdf(file="Rplots5.pdf")
> plot2d(mod2d1) ## in plane (O,X,Y)
> lines(my~mx,col=2)
> 
> ## dX(t) = 4*(-1-X(t))*Y(t) dt + 0.2 dW1(t)
> ## dY(t) = 4*(1-Y(t))*X(t) dt + 0.2 dW2(t)
> set.seed(1234)
> 
> fx <- expression(4*(-1-x)*y , 4*(1-y)*x )
> gx <- expression(0.25*y,0.2*x)
> 
> mod2d1 <- snssde2d(drift=fx,diffusion=gx,x0=c(x0=1,y0=-1),M=500)
> mod2d1
Itô Sde 2D:
	| dX(t) = 4 * (-1 - X(t)) * Y(t) * dt + 0.25 * Y(t) * dW1(t)
	| dY(t) = 4 * (1 - Y(t)) * X(t) * dt + 0.2 * X(t) * dW2(t)
Method:
	| Euler scheme with order 0.5
Summary:
	| Size of process	| N  = 1001.
	| Number of simulation	| M  = 500.
	| Initial values	| (x0,y0) = (1,-1).
	| Time of process	| t in [0,1].
	| Discretization	| Dt = 0.001.
> summary(mod2d1)

	Monte-Carlo Statistics for (X(t),Y(t)) at time t = 1
                          X          Y
Mean             -0.6952580  0.4472993
Variance          0.0057844  0.0207342
Median           -0.6964787  0.4532228
Mode             -0.6898841  0.4635814
First quartile   -0.7468943  0.3592755
Third quartile   -0.6451885  0.5386053
Minimum          -0.9951340 -0.0873122
Maximum          -0.4436166  1.0762313
Skewness         -0.0257725 -0.0227105
Kurtosis          3.4648585  4.3630757
Coef-variation   -0.1093915  0.3219185
3th-order moment -0.0000113 -0.0000678
4th-order moment  0.0001159  0.0018757
5th-order moment -0.0000024  0.0001410
6th-order moment  0.0000046  0.0003928
> ##
> dev.new()
dev.new(): using pdf(file="Rplots6.pdf")
> plot(mod2d1,type="n")
> mx <- apply(mod2d1$X,1,mean)
> my <- apply(mod2d1$Y,1,mean)
> lines(time(mod2d1),mx,col=1)
> lines(time(mod2d1),my,col=2)
> legend("topright",c(expression(E(X[t])),expression(E(Y[t]))),lty=1,inset = .01,col=c(1,2),cex=0.95)
> ##
> dev.new()
dev.new(): using pdf(file="Rplots7.pdf")
> plot2d(mod2d1) ## in plane (O,X,Y)
> lines(my~mx,col=2)
> 
> ## Now W1(t) and W2(t) are correlated.
> 
> set.seed(1234)
> Sigma <- matrix(c(0.9,0.3,0.3,0.4),nrow=2,ncol=2) # correlation matrix
> mod2d1 <- snssde2d(drift=fx,diffusion=gx,corr=Sigma,x0=c(x0=1,y0=-1),M=1000)
> mod2d1
Itô Sde 2D:
	| dX(t) = 4 * (-1 - X(t)) * Y(t) * dt + 0.25 * Y(t) * dB1(t)
	| dY(t) = 4 * (1 - Y(t)) * X(t) * dt + 0.2 * X(t) * dB2(t)
	| Correlation structure:                    
             0.9 0.3
             0.3 0.4
Method:
	| Euler scheme with order 0.5
Summary:
	| Size of process	| N  = 1001.
	| Number of simulation	| M  = 1000.
	| Initial values	| (x0,y0) = (1,-1).
	| Time of process	| t in [0,1].
	| Discretization	| Dt = 0.001.
> summary(mod2d1)

	Monte-Carlo Statistics for (X(t),Y(t)) at time t = 1
                          X         Y
Mean             -0.6932348 0.4455943
Variance          0.0065102 0.0172920
Median           -0.6896905 0.4399882
Mode             -0.6896284 0.4258462
First quartile   -0.7473370 0.3627950
Third quartile   -0.6404650 0.5321686
Minimum          -1.0387092 0.0584500
Maximum          -0.4364747 1.1318649
Skewness         -0.2131556 0.1203688
Kurtosis          3.2057513 3.6055238
Coef-variation   -0.1163899 0.2951089
3th-order moment -0.0001120 0.0002737
4th-order moment  0.0001359 0.0010781
5th-order moment -0.0000081 0.0001448
6th-order moment  0.0000054 0.0001734
> ##
> dev.new()
dev.new(): using pdf(file="Rplots8.pdf")
> plot(mod2d1,type="n")
> mx <- apply(mod2d1$X,1,mean)
> my <- apply(mod2d1$Y,1,mean)
> lines(time(mod2d1),mx,col=1)
> lines(time(mod2d1),my,col=2)
> legend("topright",c(expression(E(X[t])),expression(E(Y[t]))),lty=1,inset = .01,col=c(1,2),cex=0.95)
> ##
> dev.new()
dev.new(): using pdf(file="Rplots9.pdf")
> plot2d(mod2d1) ## in plane (O,X,Y)
> lines(my~mx,col=2)
> 
> 
> ## dX(t) = Y(t)* dt  + 0.2 o dW3(t)         
> ## dY(t) = (4*( 1-X(t)^2 )* Y(t) - X(t))* dt + 0.2 o dW2(t)
> ## dZ(t) = (4*( 1-X(t)^2 )* Z(t) - X(t))* dt + 0.2 o dW3(t)
> ## W1(t), W2(t) and W3(t) are three correlated Brownian motions with Sigma
> 
> fx <- expression( y , (4*( 1-x^2 )* y - x), (4*( 1-x^2 )* z - x))
> gx <- expression( 0.2 , 0.2, 0.2)
> # correlation matrix
> Sigma <-matrix(c(1,0.3,0.5,0.3,1,0.2,0.5,0.2,1),nrow=3,ncol=3) 
> 
> mod3d2 <- snssde3d(drift=fx,diffusion=gx,N=10000,T=100,type="str",corr=Sigma)
> mod3d2
Stratonovich Sde 3D:
	| dX(t) = Y(t) * dt + 0.2 o dB1(t)
	| dY(t) = (4 * (1 - X(t)^2) * Y(t) - X(t)) * dt + 0.2 o dB2(t)
	| dZ(t) = (4 * (1 - X(t)^2) * Z(t) - X(t)) * dt + 0.2 o dB3(t)
	| Correlation structure:                        
             1.0 0.3 0.5
             0.3 1.0 0.2
             0.5 0.2 1.0
Method:
	| Euler scheme with order 0.5
Summary:
	| Size of process	| N  = 10001.
	| Number of simulation	| M  = 1.
	| Initial values	| (x0,y0,z0) = (0,0,0).
	| Time of process	| t in [0,100].
	| Discretization	| Dt = 0.01.
> 
> proc.time()
   user  system elapsed 
   2.89    0.71    3.59