R Under development (unstable) (2025-12-12 r89163 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 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)
> 
> ## 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=200)
> 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  = 200.
	| 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.6959631  0.4524414
Variance          0.0069413  0.0189029
Median           -0.6921780  0.4527953
Mode             -0.6888319  0.4511555
First quartile   -0.7468119  0.3705037
Third quartile   -0.6320845  0.5207363
Minimum          -1.0183596 -0.0644030
Maximum          -0.4960188  0.9405639
Skewness         -0.4688270  0.0486878
Kurtosis          3.4847591  4.3331869
Coef-variation   -0.1197112  0.3038798
3th-order moment -0.0002711  0.0001265
4th-order moment  0.0001679  0.0015483
5th-order moment -0.0000223 -0.0000193
6th-order moment  0.0000087  0.0002385
> ##
> 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=200)
> 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  = 200.
	| 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.6948387  0.4452270
Variance          0.0060344  0.0189548
Median           -0.6925074  0.4510093
Mode             -0.6980912  0.4692483
First quartile   -0.7410837  0.3649011
Third quartile   -0.6388846  0.5233289
Minimum          -0.9516900 -0.0185297
Maximum          -0.4883747  0.8885734
Skewness         -0.2595185  0.0786553
Kurtosis          3.1911652  4.0353213
Coef-variation   -0.1117979  0.3092273
3th-order moment -0.0001217  0.0002053
4th-order moment  0.0001162  0.0014498
5th-order moment -0.0000069  0.0000289
6th-order moment  0.0000040  0.0001849
> ##
> 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 
   3.48    1.12    4.59