# snssde1d()

Assume that we want to describe the following SDE:

Itô form3:

$$$\label{eq:05} dX_{t} = \frac{1}{2}\theta^{2} X_{t} dt + \theta X_{t} dW_{t},\qquad X_{0}=x_{0} > 0$$$ Stratonovich form: $$$\label{eq:06} dX_{t} = \frac{1}{2}\theta^{2} X_{t} dt +\theta X_{t} \circ dW_{t},\qquad X_{0}=x_{0} > 0$$$

In the above $$f(t,x)=\frac{1}{2}\theta^{2} x$$ and $$g(t,x)= \theta x$$ ($$\theta > 0$$), $$W_{t}$$ is a standard Wiener process. To simulate this models using snssde1d() function we need to specify:

• The drift and diffusion coefficients as R expressions that depend on the state variable x and time variable t.
• The number of simulation steps N=1000 (by default: N=1000).
• The number of the solution trajectories to be simulated by M=1000 (by default: M=1).
• The initial conditions t0=0, x0=10 and end time T=1 (by default: t0=0, x0=0 and T=1).
• The integration step size Dt=0.001 (by default: Dt=(T-t0)/N).
• The choice of process types by the argument type="ito" for Ito or type="str" for Stratonovich (by default type="ito").
• The numerical method to be used by method (by default method="euler").
R> theta = 0.5
R> f <- expression( (0.5*theta^2*x) )
R> g <- expression( theta*x )
R> mod1 <- snssde1d(drift=f,diffusion=g,x0=10,M=1000,type="ito") # Using Itô
R> mod2 <- snssde1d(drift=f,diffusion=g,x0=10,M=1000,type="str") # Using Stratonovich
R> mod1
Itô Sde 1D:
| dX(t) = (0.5 * theta^2 * X(t)) * dt + theta * X(t) * dW(t)
Method:
| Euler scheme with order 0.5
Summary:
| Size of process   | N  = 1001.
| Number of simulation  | M  = 1000.
| Initial value     | x0 = 10.
| Time of process   | t in [0,1].
| Discretization    | Dt = 0.001.
R> mod2
Stratonovich Sde 1D:
| dX(t) = (0.5 * theta^2 * X(t)) * dt + theta * X(t) o dW(t)
Method:
| Euler scheme with order 0.5
Summary:
| Size of process   | N  = 1001.
| Number of simulation  | M  = 1000.
| Initial value     | x0 = 10.
| Time of process   | t in [0,1].
| Discretization    | Dt = 0.001.

Using Monte-Carlo simulations, the following statistical measures (S3 method) for class snssde1d() can be approximated for the $$X_{t}$$ process at any time $$t$$:

• The expected value $$\text{E}(X_{t})$$ at time $$t$$, using the command mean.
• The variance $$\text{Var}(X_{t})$$ at time $$t$$, using the command moment with order=2 and center=TRUE.
• The median $$\text{Med}(X_{t})$$ at time $$t$$, using the command Median.
• The mode $$\text{Mod}(X_{t})$$ at time $$t$$, using the command Mode.
• The quartile of $$X_{t}$$ at time $$t$$, using the command quantile.
• The maximum and minimum of $$X_{t}$$ at time $$t$$, using the command min and max.
• The skewness and the kurtosis of $$X_{t}$$ at time $$t$$, using the command skewness and kurtosis.
• The coefficient of variation (relative variability) of $$X_{t}$$ at time $$t$$, using the command cv.
• The central moments up to order $$p$$ of $$X_{t}$$ at time $$t$$, using the command moment.
• The empirical $$\alpha \%$$ confidence interval of expected value $$\text{E}(X_{t})$$ at time $$t$$ (from the $$2.5th$$ to the $$97.5th$$ percentile), using the command bconfint.
• The result summaries of the results of Monte-Carlo simulation at time $$t$$, using the command summary.

The following statistical measures (S3 method) for class snssde1d() can be approximated for the $$X_{t}$$ process at any time $$t$$, for example at=1:

R> s = 1
R> mean(mod1, at = s)
[1] 11.466
R> moment(mod1, at = s , center = TRUE , order = 2) ## variance
[1] 33.47
R> Median(mod1, at = s)
[1] 10.125
R> Mode(mod1, at =s)
[1] 7.6499
R> quantile(mod1 , at = s)
     0%     25%     50%     75%    100%
2.6288  7.2415 10.1250 14.2546 43.1338 
R> kurtosis(mod1 , at = s)
[1] 6.4091
R> skewness(mod1 , at = s)
[1] 1.4907
R> cv(mod1 , at = s )
[1] 0.5048
R> min(mod1 , at = s)
[1] 2.6288
R> max(mod1 , at = s)
[1] 43.134
R> moment(mod1, at = s , center= TRUE , order = 4)
[1] 7194.3
R> moment(mod1, at = s , center= FALSE , order = 4)
[1] 64144

The summary of the results of mod1 and mod2 at time $$t=1$$ of class snssde1d() is given by:

R> summary(mod1, at = 1)

Monte-Carlo Statistics for X(t) at time t = 1

Mean                  11.4664
Variance              33.5040
Median                10.1250
Mode                   7.6499
First quartile         7.2415
Third quartile        14.2546
Minimum                2.6288
Maximum               43.1338
Skewness               1.4907
Kurtosis               6.4091
Coef-variation         0.5048
3th-order moment     289.0886
4th-order moment    7194.3134
5th-order moment  155589.8634
6th-order moment 4047545.0993
R> summary(mod2, at = 1)

Monte-Carlo Statistics for X(t) at time t = 1

Mean                  10.46259
Variance              32.17703
Median                 8.96700
Mode                   7.50721
First quartile         6.39204
Third quartile        13.19293
Minimum                1.52840
Maximum               43.27138
Skewness               1.45170
Kurtosis               5.88001
Coef-variation         0.54217
3th-order moment     264.96957
4th-order moment    6087.93213
5th-order moment  120792.75449
6th-order moment 2960112.24586

Hence we can just make use of the rsde1d() function to build our random number generator for the conditional density of the $$X_{t}|X_{0}$$ ($$X_{t}^{\text{mod1}}| X_{0}$$ and $$X_{t}^{\text{mod2}}|X_{0}$$) at time $$t = 1$$.

R> x1 <- rsde1d(object = mod1, at = 1)  # X(t=1) | X(0)=x0 (Itô SDE)
R> x2 <- rsde1d(object = mod2, at = 1)  # X(t=1) | X(0)=x0 (Stratonovich SDE)
R> head(x1,n=10)
 [1]  4.8225  6.7801  7.0047 16.4111  5.9737  5.6141 19.1933  4.5259
[9] 16.5058 10.2209
R> head(x2,n=10)
 [1] 12.6038 21.4670 12.9827  6.6244 13.5564  8.1077 11.9912 14.5045
[9]  8.2545  4.9719
R> summary(data.frame(x1,x2))
       x1              x2
Min.   : 2.63   Min.   : 1.53
1st Qu.: 7.24   1st Qu.: 6.39
Median :10.12   Median : 8.97
Mean   :11.47   Mean   :10.46
3rd Qu.:14.26   3rd Qu.:13.19
Max.   :43.13   Max.   :43.27  

The function dsde1d() can be used to show the Approximate transitional density for $$X_{t}|X_{0}$$ at time $$t-s=1$$ with log-normal curves:

R> mu1 = log(10); sigma1= sqrt(theta^2)  # log mean and log variance for mod1
R> mu2 = log(10)-0.5*theta^2 ; sigma2 = sqrt(theta^2) # log mean and log variance for mod2
R> AppdensI <- dsde1d(mod1, at = 1)
R> AppdensS <- dsde1d(mod2, at = 1)
R> plot(AppdensI , dens = function(x) dlnorm(x,meanlog=mu1,sdlog = sigma1))
R> plot(AppdensS , dens = function(x) dlnorm(x,meanlog=mu2,sdlog = sigma2))

In Figure 2, we present the flow of trajectories, the mean path (red lines) of solution of and , with their empirical $$95\%$$ confidence bands, that is to say from the $$2.5th$$ to the $$97.5th$$ percentile for each observation at time $$t$$ (blue lines):

R> ## Itô
R> plot(mod1,ylab=expression(X^mod1))
R> lines(time(mod1),apply(mod1$X,1,mean),col=2,lwd=2) R> lines(time(mod1),apply(mod1$X,1,bconfint,level=0.95)[1,],col=4,lwd=2)
R> lines(time(mod1),apply(mod1$X,1,bconfint,level=0.95)[2,],col=4,lwd=2) R> legend("topleft",c("mean path",paste("bound of", 95,"% confidence")),inset = .01,col=c(2,4),lwd=2,cex=0.8) R> ## Stratonovich R> plot(mod2,ylab=expression(X^mod2)) R> lines(time(mod2),apply(mod2$X,1,mean),col=2,lwd=2)
R> lines(time(mod2),apply(mod2$X,1,bconfint,level=0.95)[1,],col=4,lwd=2) R> lines(time(mod2),apply(mod2$X,1,bconfint,level=0.95)[2,],col=4,lwd=2)
R> legend("topleft",c("mean path",paste("bound of",95,"% confidence")),col=c(2,4),inset =.01,lwd=2,cex=0.8)

# snssde2d()

The following $$2$$-dimensional SDE’s with a vector of drift and a diagonal matrix of diffusion coefficients:

Itô form: $$$\label{eq:09} \begin{cases} dX_t = f_{x}(t,X_{t},Y_{t}) dt + g_{x}(t,X_{t},Y_{t}) dW_{1,t}\\ dY_t = f_{y}(t,X_{t},Y_{t}) dt + g_{y}(t,X_{t},Y_{t}) dW_{2,t} \end{cases}$$$ Stratonovich form: $$$\label{eq:10} \begin{cases} dX_t = f_{x}(t,X_{t},Y_{t}) dt + g_{x}(t,X_{t},Y_{t}) \circ dW_{1,t}\\ dY_t = f_{y}(t,X_{t},Y_{t}) dt + g_{y}(t,X_{t},Y_{t}) \circ dW_{2,t} \end{cases}$$$

$$W_{1,t}$$ and $$W_{2,t}$$ is a two independent standard Wiener process. To simulate $$2d$$ models using snssde2d() function we need to specify:

• The drift (2d) and diffusion (2d) coefficients as R expressions that depend on the state variable x, y and time variable t.
• The number of simulation steps N (default: N=1000).
• The number of the solution trajectories to be simulated by M (default: M=1).
• The initial conditions t0, x0 and end time T (default: t0=0, x0=c(0,0) and T=1).
• The integration step size Dt (default: Dt=(T-t0)/N).
• The choice of process types by the argument type="ito" for Ito or type="str" for Stratonovich (default type="ito").
• The numerical method to be used by method (default method="euler").

## Ornstein-Uhlenbeck process and its integral

The Ornstein-Uhlenbeck (OU) process has a long history in physics. Introduced in essence by Langevin in his famous 1908 paper on Brownian motion, the process received a more thorough mathematical examination several decades later by Uhlenbeck and Ornstein (1930). The OU process is understood here to be the univariate continuous Markov process $$X_t$$. In mathematical terms, the equation is written as an Ito equation: $$$\label{eq016} dX_t = -\frac{1}{\mu} X_t dt + \sqrt{\sigma} dW_t,\quad X_{0}=x_{0}$$$ In these equations, $$\mu$$ and $$\sigma$$ are positive constants called, respectively, the relaxation time and the diffusion constant. The time integral of the OU process $$X_t$$ (or indeed of any process $$X_t$$) is defined to be the process $$Y_t$$ that satisfies: $$$\label{eq017} Y_{t} = Y_{0}+\int X_{t} dt \Leftrightarrow dY_t = X_{t} dt ,\quad Y_{0}=y_{0}$$$ $$Y_t$$ is not itself a Markov process; however, $$X_t$$ and $$Y_t$$ together comprise a bivariate continuous Markov process. We wish to find the solutions $$X_t$$ and $$Y_t$$ to the coupled time-evolution equations: $$$\label{eq018} \begin{cases} dX_t = -\frac{1}{\mu} X_t dt + \sqrt{\sigma} dW_t\\ dY_t = X_{t} dt \end{cases}$$$

We simulate a flow of $$1000$$ trajectories of $$(X_{t},Y_{t})$$, with integration step size $$\Delta t = 0.01$$, and using second Milstein method.

R> x0=5;y0=0
R> mu=3;sigma=0.5
R> fx <- expression(-(x/mu),x)
R> gx <- expression(sqrt(sigma),0)
R> mod2d <- snssde2d(drift=fx,diffusion=gx,Dt=0.01,M=1000,x0=c(x0,y0),method="smilstein")
R> mod2d
Itô Sde 2D:
| dX(t) = -(X(t)/mu) * dt + sqrt(sigma) * dW1(t)
| dY(t) = X(t) * dt + 0 * dW2(t)
Method:
| Second-order Milstein scheme
Summary:
| Size of process   | N  = 1001.
| Number of simulation  | M  = 1000.
| Initial values    | (x0,y0) = (5,0).
| Time of process   | t in [0,10].
| Discretization    | Dt = 0.01.

The following statistical measures (S3 method) for class snssde2d() can be approximated for the $$(X_{t},Y_{t})$$ process at any time $$t$$, for example at=5:

R> s = 5
R> mean(mod2d, at = s)
[1]  1.0069 12.3221
R> moment(mod2d, at = s , center = TRUE , order = 2) ## variance
[1] 0.74904 6.87244
R> Median(mod2d, at = s)
[1]  1.0139 12.2275
R> Mode(mod2d, at = s)
[1]  0.96332 11.76403
R> quantile(mod2d , at = s)
$x 0% 25% 50% 75% 100% -1.26391 0.41043 1.01394 1.60721 3.78758$y
0%     25%     50%     75%    100%
3.0444 10.6755 12.2275 14.0689 20.0437 
R> kurtosis(mod2d , at = s)
[1] 2.8872 3.0034
R> skewness(mod2d , at = s)
[1]  0.037126 -0.093784
R> cv(mod2d , at = s )
[1] 0.85998 0.21286
R> min(mod2d , at = s)
[1] -1.2639  3.0444
R> max(mod2d , at = s)
[1]  3.7876 20.0437
R> moment(mod2d, at = s , center= TRUE , order = 4)
[1]   1.6232 142.1364
R> moment(mod2d, at = s , center= FALSE , order = 4)
[1]     7.3044 29372.9369

The summary of the results of mod2d at time $$t=5$$ of class snssde2d() is given by:

R> summary(mod2d, at = s)

Monte-Carlo Statistics for (X(t),Y(t)) at time t = 5
X          Y
Mean              1.00689   12.32207
Variance          0.74979    6.87932
Median            1.01394   12.22751
Mode              0.96332   11.76403
First quartile    0.41043   10.67551
Third quartile    1.60721   14.06895
Minimum          -1.26391    3.04445
Maximum           3.78758   20.04373
Skewness          0.03713   -0.09378
Kurtosis          2.88723    3.00340
Coef-variation    0.85998    0.21286
3th-order moment  0.02410   -1.69217
4th-order moment  1.62316  142.13642
5th-order moment  0.32989 -140.26206
6th-order moment  5.57797 4744.19317

For plotting (back in time) using the command plot, the results of the simulation are shown in Figure 3.

R> plot(mod2d)

Take note of the well known result, which can be derived from either this equations. That for any $$t > 0$$ the OU process $$X_t$$ and its integral $$Y_t$$ will be the normal distribution with mean and variance given by: $\begin{cases} \text{E}(X_{t}) =x_{0} e^{-t/\mu} &\text{and}\quad\text{Var}(X_{t})=\frac{\sigma \mu}{2} \left (1-e^{-2t/\mu}\right )\\ \text{E}(Y_{t}) = y_{0}+x_{0}\mu \left (1-e^{-t/\mu}\right ) &\text{and}\quad\text{Var}(Y_{t})=\sigma\mu^{3}\left (\frac{t}{\mu}-2\left (1-e^{-t/\mu}\right )+\frac{1}{2}\left (1-e^{-2t/\mu}\right )\right ) \end{cases}$

Hence we can just make use of the rsde2d() function to build our random number for $$(X_{t},Y_{t})$$ at time $$t = 10$$.

R> out <- rsde2d(object = mod2d, at = 10)
R> head(out,n=10)
           x       y
1  -0.688860 17.6594
2   0.503306 11.4270
3   0.845208 18.3864
4  -0.386227 23.0378
5   0.066081  9.3175
6  -0.018649 19.9152
7  -1.566948 11.4247
8   1.277140 20.4721
9   0.039498  8.4714
10 -1.127192  6.2833
R> summary(out)
       x                y
Min.   :-3.471   Min.   :-1.2
1st Qu.:-0.439   1st Qu.:11.3
Median : 0.173   Median :14.7
Mean   : 0.159   Mean   :14.8
3rd Qu.: 0.762   3rd Qu.:18.3
Max.   : 3.125   Max.   :29.7  
R> cov(out)
        x       y
x 0.79238  2.2424
y 2.24243 24.9690

Figure 4, show simulation results for moments of system :

R> mx <- apply(mod2d$X,1,mean) R> my <- apply(mod2d$Y,1,mean)
R> Sx <- apply(mod2d$X,1,var) R> Sy <- apply(mod2d$Y,1,var)
R> Cxy <- sapply(1:1001,function(i) cov(mod2d$X[i,],mod2d$Y[i,]))
R> out_b <- data.frame(mx,my,Sx,Sy,Cxy)
R> matplot(time(mod2d), out_b, type = "l", xlab = "time", ylab = "",col=2:6,lwd=2,lty=2:6,las=1)
R> legend("topleft",c(expression(hat(E)(X[t]),hat(E)(Y[t]),hat(Var)(X[t]),hat(Var)(Y[t]),hat(Cov)(X[t],Y[t]))),inset = .05,col=2:6,lty=2:6,lwd=2,cex=0.9)

For each SDE type and for each numerical scheme, the density of $$X_t$$ and $$Y_t$$ at time $$t=10$$ are reported using dsde2d() function, see e.g. Figure 5: the marginal density of $$X_t$$ and $$Y_t$$ at time $$t=10$$.

R> denM <- dsde2d(mod2d,pdf="M",at =10)
R> denM

Marginal density of X(t-t0)|X(t0)=5 at time t = 10

Data: x (1000 obs.);    Bandwidth 'bw' = 0.2012

x                f(x)
Min.   :-4.0747   Min.   :0.00002
1st Qu.:-2.1238   1st Qu.:0.00245
Median :-0.1730   Median :0.04302
Mean   :-0.1730   Mean   :0.12802
3rd Qu.: 1.7778   3rd Qu.:0.26022
Max.   : 3.7287   Max.   :0.42696

Marginal density of Y(t-t0)|Y(t0)=0 at time t = 10

Data: y (1000 obs.);    Bandwidth 'bw' = 1.13

y               f(y)
Min.   :-4.591   Min.   :0.000004
1st Qu.: 4.827   1st Qu.:0.001830
Median :14.245   Median :0.013030
Mean   :14.245   Mean   :0.026520
3rd Qu.:23.663   3rd Qu.:0.054167
Max.   :33.080   Max.   :0.076298  
R> plot(denM, main="Marginal Density")

Created using dsde2d() plotted in (x, y)-space with dim = 2. A contour and image plot of density obtained from a realization of system $$(X_{t},Y_{t})$$ at time t=10.

R> denJ <- dsde2d(mod2d, pdf="J", n=100,at =10)
R> denJ

Joint density of (X(t-t0),Y(t-t0)|X(t0)=5,Y(t0)=0) at time t = 10

Data: (x,y) (2 x 1000 obs.);

x                 y               f(x,y)
Min.   :-3.4710   Min.   :-1.2016   Min.   :0.000000
1st Qu.:-1.8220   1st Qu.: 6.5217   1st Qu.:0.000015
Median :-0.1730   Median :14.2450   Median :0.000758
Mean   :-0.1730   Mean   :14.2450   Mean   :0.004797
3rd Qu.: 1.4760   3rd Qu.:21.9683   3rd Qu.:0.005098
Max.   : 3.1249   Max.   :29.6916   Max.   :0.037473  
R> plot(denJ,display="contour",main="Bivariate Transition Density at time t=10")
R> plot(denJ,display="image",main="Bivariate Transition Density at time t=10")

A $$3$$D plot of the transition density at $$t=10$$ obtained with:

R> plot(denJ,main="Bivariate Transition Density at time t=10")

We approximate the bivariate transition density over the set transition horizons $$t\in [1,10]$$ by $$\Delta t = 0.005$$ using the code:

R> 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))
+ }

## The stochastic Van-der-Pol equation

The Van der Pol (1922) equation is an ordinary differential equation that can be derived from the Rayleigh differential equation by differentiating and setting $$\dot{x}=y$$, see Naess and Hegstad (1994); Leung (1995) and for more complex dynamics in Van-der-Pol equation see Jing et al. (2006). It is an equation describing self-sustaining oscillations in which energy is fed into small oscillations and removed from large oscillations. This equation arises in the study of circuits containing vacuum tubes and is given by: $$$\label{eq:12} \ddot{X}-\mu (1-X^{2}) \dot{X} + X = 0$$$ where $$x$$ is the position coordinate (which is a function of the time $$t$$), and $$\mu$$ is a scalar parameter indicating the nonlinearity and the strength of the damping, to simulate the deterministic equation see Grayling (2014) for more details. Consider stochastic perturbations of the Van-der-Pol equation, and random excitation force of such systems by White noise $$\xi_{t}$$, with delta-type correlation function $$\text{E}(\xi_{t}\xi_{t+h})=2\sigma \delta (h)$$ $$$\label{eq:13} \ddot{X}-\mu (1-X^{2}) \dot{X} + X = \xi_{t},$$$ where $$\mu > 0$$ . It’s solution cannot be obtained in terms of elementary functions, even in the phase plane. The White noise $$\xi_{t}$$ is formally derivative of the Wiener process $$W_{t}$$. The representation of a system of two first order equations follows the same idea as in the deterministic case by letting $$\dot{x}=y$$, from physical equation we get the above system: $$$\label{eq:14} \begin{cases} \dot{X} = Y \\ \dot{Y} = \mu \left(1-X^{2}\right) Y - X + \xi_{t} \end{cases}$$$ The system can mathematically explain by a Stratonovitch equations: $$$\label{eq:15} \begin{cases} dX_{t} = Y_{t} dt \\ dY_{t} = \left(\mu (1-X^{2}_{t}) Y_{t} - X_{t}\right) dt + 2 \sigma \circ dW_{2,t} \end{cases}$$$

Implemente in R as follows, with integration step size $$\Delta t = 0.01$$ and using stochastic Runge-Kutta methods 1-stage.

R> mu = 4; sigma=0.1
R> fx <- expression( y ,  (mu*( 1-x^2 )* y - x))
R> gx <- expression( 0 ,2*sigma)
R> mod2d <- snssde2d(drift=fx,diffusion=gx,N=10000,Dt=0.01,type="str",method="rk1")
R> mod2d
Stratonovich Sde 2D:
| dX(t) = Y(t) * dt + 0 o dW1(t)
| dY(t) = (mu * (1 - X(t)^2) * Y(t) - X(t)) * dt + 2 * sigma o dW2(t)
Method:
| Runge-Kutta method with order 1
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.

For plotting (back in time) using the command plot, and plot2d in plane the results of the simulation are shown in Figure 8.

R> plot2d(mod2d) ## in plane (O,X,Y)
R> plot(mod2d)   ## back in time

# snssde3d()

The following $$3$$-dimensional SDE’s with a vector of drift and a diagonal matrix of diffusion coefficients:

Itô form: $$$\label{eq17} \begin{cases} dX_t = f_{x}(t,X_{t},Y_{t},Z_{t}) dt + g_{x}(t,X_{t},Y_{t},Z_{t}) dW_{1,t}\\ dY_t = f_{y}(t,X_{t},Y_{t},Z_{t}) dt + g_{y}(t,X_{t},Y_{t},Z_{t}) dW_{2,t}\\ dZ_t = f_{z}(t,X_{t},Y_{t},Z_{t}) dt + g_{z}(t,X_{t},Y_{t},Z_{t}) dW_{3,t} \end{cases}$$$ Stratonovich form: $$$\label{eq18} \begin{cases} dX_t = f_{x}(t,X_{t},Y_{t},Z_{t}) dt + g_{x}(t,X_{t},Y_{t},Z_{t}) \circ dW_{1,t}\\ dY_t = f_{y}(t,X_{t},Y_{t},Z_{t}) dt + g_{y}(t,X_{t},Y_{t},Z_{t}) \circ dW_{2,t}\\ dZ_t = f_{z}(t,X_{t},Y_{t},Z_{t}) dt + g_{z}(t,X_{t},Y_{t},Z_{t}) \circ dW_{3,t} \end{cases}$$$

$$W_{1,t}$$, $$W_{2,t}$$ and $$W_{3,t}$$ is a 3 independent standard Wiener process. To simulate this system using snssde3d() function we need to specify:

• The drift (3d) and diffusion (3d) coefficients as R expressions that depend on the state variables x, y , z and time variable t.
• The number of simulation steps N (default: N=1000).
• The number of the solution trajectories to be simulated by M (default: M=1).
• The initial conditions t0, x0 and end time T (default: t0=0, x0=c(0,0,0) and T=1).
• The integration step size Dt (default: Dt=(T-t0)/N).
• The choice of process types by the argument type="ito" for Ito or type="str" for Stratonovich (default type="ito").
• The numerical method to be used by method (default method="euler").

## Basic example

Assume that we want to describe the following SDE’s (3D) in Itô form: $$$\label{eq0166} \begin{cases} dX_t = 4 (-1-X_{t}) Y_{t} dt + 0.2 dW_{1,t}\\ dY_t = 4 (1-Y_{t}) X_{t} dt + 0.2 dW_{2,t}\\ dZ_t = 4 (1-Z_{t}) Y_{t} dt + 0.2 dW_{3,t} \end{cases}$$$

We simulate a flow of $$10000$$ trajectories, with integration step size $$\Delta t = 0.001$$.

R> fx <- expression(4*(-1-x)*y , 4*(1-y)*x , 4*(1-z)*y)
R> gx <- rep(expression(0.2),3)
R> mod3d <- snssde3d(x0=c(x=2,y=-2,z=-2),drift=fx,diffusion=gx,M=1000)
R> mod3d
Itô Sde 3D:
| 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)
| dZ(t) = 4 * (1 - Z(t)) * Y(t) * dt + 0.2 * dW3(t)
Method:
| Euler scheme with order 0.5
Summary:
| Size of process   | N  = 1001.
| Number of simulation  | M  = 1000.
| Initial values    | (x0,y0,z0) = (2,-2,-2).
| Time of process   | t in [0,1].
| Discretization    | Dt = 0.001.

The following statistical measures (S3 method) for class snssde3d() can be approximated for the $$(X_{t},Y_{t},Z_{t})$$ process at any time $$t$$, for example at=1:

R> s = 1
R> mean(mod3d, at = s)
[1] -0.79844  0.88219  0.79220
R> moment(mod3d, at = s , center = TRUE , order = 2) ## variance
[1] 0.0095972 0.1060181 0.0101603
R> Median(mod3d, at = s)
[1] -0.80094  0.84685  0.80026
R> Mode(mod3d, at = s)
[1] -0.79349  0.78589  0.81356
R> quantile(mod3d , at = s)
$x 0% 25% 50% 75% 100% -1.11138 -0.86893 -0.80094 -0.73366 -0.45868$y
0%       25%       50%       75%      100%
0.0074104 0.6667170 0.8468463 1.0820634 2.1778814

$z 0% 25% 50% 75% 100% 0.41228 0.73257 0.80026 0.86194 1.10379  R> kurtosis(mod3d , at = s) [1] 3.0392 3.3565 3.2356 R> skewness(mod3d , at = s) [1] 0.21678 0.47794 -0.40900 R> cv(mod3d , at = s ) [1] -0.12276 0.36927 0.12730 R> min(mod3d , at = s) [1] -1.1113848 0.0074104 0.4122806 R> max(mod3d , at = s) [1] -0.45868 2.17788 1.10379 R> moment(mod3d, at = s , center= TRUE , order = 4) [1] 0.00028049 0.03780189 0.00033469 R> moment(mod3d, at = s , center= FALSE , order = 4) [1] 0.44275 1.19684 0.43112 The summary of the results of mod3d at time $$t=1$$ of class snssde3d() is given by: R> summary(mod3d, at = s)  Monte-Carlo Statistics for (X(t),Y(t),Z(t)) at time t = 1 X Y Z Mean -0.79844 0.88219 0.79220 Variance 0.00961 0.10612 0.01017 Median -0.80094 0.84685 0.80026 Mode -0.79349 0.78589 0.81356 First quartile -0.86893 0.66672 0.73257 Third quartile -0.73366 1.08206 0.86194 Minimum -1.11138 0.00741 0.41228 Maximum -0.45868 2.17788 1.10379 Skewness 0.21678 0.47794 -0.40900 Kurtosis 3.03922 3.35648 3.23564 Coef-variation -0.12276 0.36927 0.12730 3th-order moment 0.00020 0.01652 -0.00042 4th-order moment 0.00028 0.03780 0.00033 5th-order moment 0.00002 0.01752 -0.00004 6th-order moment 0.00001 0.02559 0.00002 For plotting (back in time) using the command plot, and plot3D in space the results of the simulation are shown in Figure 9. R> plot(mod3d,union = TRUE) ## back in time R> plot3D(mod3d,display="persp") ## in space (O,X,Y,Z) Hence we can just make use of the rsde3d() function to build our random number for $$(X_{t},Y_{t},Z_{t})$$ at time $$t = 1$$. R> out <- rsde3d(object = mod3d, at = s) R> head(out,n=10)  x y z 1 -0.85699 1.00533 0.89567 2 -0.96721 0.93228 0.84998 3 -0.64995 0.68665 0.75584 4 -0.78889 0.80664 0.86423 5 -0.80354 1.51718 0.97352 6 -0.82615 1.29612 0.95568 7 -0.67545 0.58208 0.87521 8 -0.88578 1.18595 0.87038 9 -0.88253 1.34115 0.94522 10 -0.72538 1.25413 0.86337 R> summary(out)  x y z Min. :-1.111 Min. :0.00741 Min. :0.412 1st Qu.:-0.869 1st Qu.:0.66672 1st Qu.:0.733 Median :-0.801 Median :0.84685 Median :0.800 Mean :-0.798 Mean :0.88219 Mean :0.792 3rd Qu.:-0.734 3rd Qu.:1.08206 3rd Qu.:0.862 Max. :-0.459 Max. :2.17788 Max. :1.104  R> cov(out)  x y z x 0.0096068 -0.017303 -0.0041835 y -0.0173035 0.106124 0.0198063 z -0.0041835 0.019806 0.0101705 For each SDE type and for each numerical scheme, the marginal density of $$X_t$$, $$Y_t$$ and $$Z_t$$ at time $$t=1$$ are reported using dsde3d() function, see e.g. Figure 10. R> den <- dsde3d(mod3d,pdf="M",at =1) R> den  Marginal density of X(t-t0)|X(t0)=2 at time t = 1 Data: x (1000 obs.); Bandwidth 'bw' = 0.02216 x f(x) Min. :-1.17786 Min. :0.0002 1st Qu.:-0.98145 1st Qu.:0.0677 Median :-0.78503 Median :0.5272 Mean :-0.78503 Mean :1.2716 3rd Qu.:-0.58862 3rd Qu.:2.5228 Max. :-0.39221 Max. :4.0482 Marginal density of Y(t-t0)|Y(t0)=-2 at time t = 1 Data: y (1000 obs.); Bandwidth 'bw' = 0.07007 y f(y) Min. :-0.20281 Min. :0.00007 1st Qu.: 0.44492 1st Qu.:0.01690 Median : 1.09265 Median :0.17064 Mean : 1.09265 Mean :0.38559 3rd Qu.: 1.74037 3rd Qu.:0.70762 Max. : 2.38810 Max. :1.33052 Marginal density of Z(t-t0)|Z(t0)=-2 at time t = 1 Data: z (1000 obs.); Bandwidth 'bw' = 0.02182 z f(z) Min. :0.34681 Min. :0.0002 1st Qu.:0.55242 1st Qu.:0.0387 Median :0.75803 Median :0.5519 Mean :0.75803 Mean :1.2147 3rd Qu.:0.96365 3rd Qu.:2.1898 Max. :1.16926 Max. :4.1812  R> plot(den, main="Marginal Density")  For an approximate joint transition density for $$(X_t,Y_t,Z_t)$$ (for more details, see package sm or ks.) R> denJ <- dsde3d(mod3d,pdf="J") R> plot(denJ,display="rgl") Return to snssde3d() ## Attractive model for 3D diffusion processes If we assume that $$U_w( x , y , z , t )$$, $$V_w( x , y , z , t )$$ and $$S_w( x , y , z , t )$$ are neglected and the dispersion coefficient $$D( x , y , z )$$ is constant. A system becomes (see Boukhetala,1996): $\begin{eqnarray}\label{eq19} % \nonumber to remove numbering (before each equation) \begin{cases} dX_t = \left(\frac{-K X_{t}}{\sqrt{X^{2}_{t} + Y^{2}_{t} + Z^{2}_{t}}}\right) dt + \sigma dW_{1,t} \nonumber\\ dY_t = \left(\frac{-K Y_{t}}{\sqrt{X^{2}_{t} + Y^{2}_{t} + Z^{2}_{t}}}\right) dt + \sigma dW_{2,t} \\ dZ_t = \left(\frac{-K Z_{t}}{\sqrt{X^{2}_{t} + Y^{2}_{t} + Z^{2}_{t}}}\right) dt + \sigma dW_{3,t} \nonumber \end{cases} \end{eqnarray}$ with initial conditions $$(X_{0},Y_{0},Z_{0})=(1,1,1)$$, by specifying the drift and diffusion coefficients of three processes $$X_{t}$$, $$Y_{t}$$ and $$Z_{t}$$ as R expressions which depends on the three state variables (x,y,z) and time variable t, with integration step size Dt=0.0001. R> K = 4; s = 1; sigma = 0.2 R> 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)) ) R> gx <- rep(expression(sigma),3) R> mod3d <- snssde3d(drift=fx,diffusion=gx,N=10000,x0=c(x=1,y=1,z=1)) R> mod3d Itô Sde 3D: | dX(t) = (-K * X(t)/sqrt(X(t)^2 + Y(t)^2 + Z(t)^2)) * dt + sigma * dW1(t) | dY(t) = (-K * Y(t)/sqrt(X(t)^2 + Y(t)^2 + Z(t)^2)) * dt + sigma * dW2(t) | dZ(t) = (-K * Z(t)/sqrt(X(t)^2 + Y(t)^2 + Z(t)^2)) * dt + sigma * dW3(t) Method: | Euler scheme with order 0.5 Summary: | Size of process | N = 10001. | Number of simulation | M = 1. | Initial values | (x0,y0,z0) = (1,1,1). | Time of process | t in [0,1]. | Discretization | Dt = 0.0001. The results of simulation are shown: R> plot3D(mod3d,display="persp",col="blue") Return to snssde3d() ## Transformation of an SDE one-dimensional Next is an example of one-dimensional SDE driven by three independent Brownian motions ($$W_{1,t}$$,$$W_{2,t}$$,$$W_{3,t}$$), as follows: $$$\label{eq20} dX_{t} = \mu W_{1,t} dt + \sigma W_{2,t} \circ dW_{3,t}$$$ To simulate the solution of the process $$X_t$$, we make a transformation to a system of three equations as follows: $\begin{eqnarray}\label{eq21} \begin{cases} % \nonumber to remove numbering (before each equation) dX_t = \mu Y_{t} dt + \sigma Z_{t} \circ dW_{3,t} \nonumber\\ dY_t = dW_{1,t} \\ dZ_t = dW_{2,t} \nonumber \end{cases} \end{eqnarray}$ run by calling the function snssde3d() to produce a simulation of the solution, with $$\mu = 1$$ and $$\sigma = 1$$. R> fx <- expression(y,0,0) R> gx <- expression(z,1,1) R> modtra <- snssde3d(drift=fx,diffusion=gx,M=1000,type="str") R> modtra Stratonovich Sde 3D: | dX(t) = Y(t) * dt + Z(t) o dW1(t) | dY(t) = 0 * dt + 1 o dW2(t) | dZ(t) = 0 * dt + 1 o dW3(t) Method: | Euler scheme with order 0.5 Summary: | Size of process | N = 1001. | Number of simulation | M = 1000. | Initial values | (x0,y0,z0) = (0,0,0). | Time of process | t in [0,1]. | Discretization | Dt = 0.001. R> summary(modtra)  Monte-Carlo Statistics for (X(t),Y(t),Z(t)) at time t = 1 X Y Z Mean -0.02621 0.01128 -0.00189 Variance 0.75292 0.90784 0.94256 Median 0.01526 -0.03007 0.01229 Mode 0.20482 -0.34317 0.06299 First quartile -0.53526 -0.58976 -0.66554 Third quartile 0.48897 0.64642 0.69749 Minimum -4.40578 -3.13029 -3.29677 Maximum 3.15060 3.03938 3.04791 Skewness -0.46812 0.00297 0.01273 Kurtosis 5.01693 3.01359 2.96985 Coef-variation -33.10636 84.50443 -513.94224 3th-order moment -0.30583 0.00257 0.01165 4th-order moment 2.84403 2.48370 2.63848 5th-order moment -4.78544 -0.16245 -0.13013 6th-order moment 26.27889 10.91901 12.27323 the following code produces the result in Figure 12. R> plot(modtra$X,plot.type="single",ylab=expression(X[t]))
R> lines(time(modtra),apply(modtra$X,1,mean),col=2,lwd=2) R> legend("topleft",c("mean path"),col=2,lwd=2,cex=0.8) The histogram and kernel density of $$X_t$$ at time $$t=1$$ are reported using dsde3d() function, see e.g. Figure 13. R> den <- dsde3d(modtra,pdf="Marginal",at=1) R> den$resx

Call:
density.default(x = x, na.rm = TRUE)

Data: x (1000 obs.);    Bandwidth 'bw' = 0.1728

x                y
Min.   :-4.924   Min.   :0.00003
1st Qu.:-2.776   1st Qu.:0.00357
Median :-0.628   Median :0.01677
Mean   :-0.628   Mean   :0.11626
3rd Qu.: 1.521   3rd Qu.:0.17825
Max.   : 3.669   Max.   :0.51356  
R> MASS::truehist(den$ech$x,xlab = expression(X[t==1]));box()
R> lines(den$resx,col="red",lwd=2) R> 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) Figure 14 and 15, show approximation results for $$m_{1}(t)= \text{E}(X_{t})$$, $$S_{1}(t)=\text{V}(X_{t})$$ and $$C(s,t)=\text{Cov}(X_{s},X_{t})$$: R> m1 <- apply(modtra$X,1,mean) ## m1(t)
R> S1  <- apply(modtra$X,1,var) ## s1(t) R> out_a <- data.frame(m1,S1) R> matplot(time(modtra), out_a, type = "l", xlab = "time", ylab = "", col=2:3,lwd=2,lty=2:3,las=1) R> legend("topleft",c(expression(m[1](t),S[1](t))),inset = .09,col=2:3,lty=2:3,lwd=2,cex=0.9) R> color.palette=colorRampPalette(c('white','green','blue','red')) R> 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 = ""))

# References

1. Boukhetala K (1996). Modelling and Simulation of a Dispersion Pollutant with Attractive Centre, volume 3, pp. 245-252. Computer Methods and Water Resources, Computational Mechanics Publications, Boston, USA.

2. Guidoum AC, Boukhetala K (2018). Performing Parallel Monte Carlo and Moment Equations Methods for Itô and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc. Preprint submitted to Journal of Statistical Software.

3. Guidoum AC, Boukhetala K (2018). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.2, URL https://cran.r-project.org/package=Sim.DiffProc.

1. Department of Probabilities & Statistics, Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, E-mail (acguidoum@usthb.dz)

2. 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)

3. The equivalently of $$X_{t}^{\text{mod1}}$$ the following Stratonovich SDE: $$dX_{t} = \theta X_{t} \circ dW_{t}$$.