Monte-Carlo Simulation and Analysis of Stochastic Differential Equations

A.C. Guidoum1 and K. Boukhetala2

2017-09-16

snssde1d()

Assume that we want to describe the following SDE:

Itô form3:

\[\begin{equation}\label{eq:05} dX_{t} = \frac{1}{2}\theta^{2} X_{t} dt + \theta X_{t} dW_{t},\qquad X_{0}=x_{0} > 0 \end{equation}\] Stratonovich form: \[\begin{equation}\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 \end{equation}\]

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:

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=10000,type="ito") # Using Itô
R> mod2 <- snssde1d(drift=f,diffusion=g,x0=10,M=10000,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  = 10000.
    | 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  = 10000.
    | 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 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.242
R> moment(mod1, at = s , center = TRUE , order = 2) ## variance
[1] 35.939
R> Median(mod1, at = s)
[1] 9.9465
R> Mode(mod1, at =s)
[1] 7.3233
R> quantile(mod1 , at = s)
     0%     25%     50%     75%    100% 
 1.6649  7.0876  9.9465 13.7595 60.5808 
R> kurtosis(mod1 , at = s)
[1] 8.3946
R> skewness(mod1 , at = s)
[1] 1.7425
R> cv(mod1 , at = s )
[1] 0.53331
R> min(mod1 , at = s)
[1] 1.6649
R> max(mod1 , at = s)
[1] 60.581
R> moment(mod1, at = s , center= TRUE , order = 4)
[1] 10845
R> moment(mod1, at = s , center= FALSE , order = 4)
[1] 70949

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.24150
Variance               35.94280
Median                  9.94653
Mode                    7.32326
First quartile          7.08756
Third quartile         13.75953
Minimum                 1.66495
Maximum                60.58077
Skewness                1.74254
Kurtosis                8.39459
Coef-variation          0.53331
3th-order moment      375.49254
4th-order moment    10844.84504
5th-order moment   311518.31516
6th-order moment 11081160.20971
R> summary(mod2, at = 1)

    Monte-Carlo Statistics for X(t) at time t = 1
                              
Mean                   9.93652
Variance              28.00642
Median                 8.78940
Mode                   7.32584
First quartile         6.26875
Third quartile        12.25361
Minimum                1.56960
Maximum               59.93883
Skewness               1.76422
Kurtosis               8.77355
Coef-variation         0.53259
3th-order moment     261.48030
4th-order moment    6881.61371
5th-order moment  181855.81506
6th-order moment 5991430.46438

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]  6.6897 13.4024  6.1000 11.6483 11.3254  8.6207 18.5029 16.6387
 [9] 12.8814  7.4597
R> head(x2,n=10)
 [1] 13.0135  2.9416  6.7040 16.4737  9.7548 10.2588 11.7907 12.8285
 [9]  7.7060  4.9572
R> summary(data.frame(x1,x2))
       x1              x2       
 Min.   : 1.66   Min.   : 1.57  
 1st Qu.: 7.09   1st Qu.: 6.27  
 Median : 9.95   Median : 8.79  
 Mean   :11.24   Mean   : 9.94  
 3rd Qu.:13.76   3rd Qu.:12.25  
 Max.   :60.58   Max.   :59.94  

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)
mod1: Itô and mod2: Stratonovich

mod1: Itô and mod2: Stratonovich

Return to snssde1d()

snssde2d()

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

Itô form: \[\begin{equation}\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} \end{equation}\] Stratonovich form: \[\begin{equation}\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} \end{equation}\]

\(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:

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: \[\begin{equation}\label{eq016} dX_t = -\frac{1}{\mu} X_t dt + \sqrt{\sigma} dW_t,\quad X_{0}=x_{0} \end{equation}\] 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: \[\begin{equation}\label{eq017} Y_{t} = Y_{0}+\int X_{t} dt \Leftrightarrow dY_t = X_{t} dt ,\quad Y_{0}=y_{0} \end{equation}\] \(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: \[\begin{equation}\label{eq018} \begin{cases} dX_t = -\frac{1}{\mu} X_t dt + \sqrt{\sigma} dW_t\\ dY_t = X_{t} dt \end{cases} \end{equation}\]

We simulate a flow of \(10000\) 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=10000,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  = 10000.
    | 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]  0.95114 12.20124
R> moment(mod2d, at = s , center = TRUE , order = 2) ## variance
[1] 0.73635 7.21815
R> Median(mod2d, at = s)
[1]  0.95203 12.20370
R> Mode(mod2d, at = s)
[1]  1.1064 12.2104
R> quantile(mod2d , at = s)
$x
      0%      25%      50%      75%     100% 
-2.41344  0.37763  0.95203  1.52497  4.52125 

$y
     0%     25%     50%     75%    100% 
 2.3806 10.3706 12.2037 13.9978 22.5233 
R> kurtosis(mod2d , at = s)
[1] 3.070 2.997
R> skewness(mod2d , at = s)
[1]  0.0030583 -0.0014066
R> cv(mod2d , at = s )
[1] 0.90223 0.22021
R> min(mod2d , at = s)
[1] -2.4134  2.3806
R> max(mod2d , at = s)
[1]  4.5212 22.5233
R> moment(mod2d, at = s , center= TRUE , order = 4)
[1]   1.6649 156.1783
R> moment(mod2d, at = s , center= FALSE , order = 4)
[1]     6.4876 28764.6272

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              0.95114   12.20124
Variance          0.73642    7.21887
Median            0.95203   12.20370
Mode              1.10642   12.21043
First quartile    0.37763   10.37056
Third quartile    1.52497   13.99781
Minimum          -2.41344    2.38061
Maximum           4.52125   22.52335
Skewness          0.00306   -0.00141
Kurtosis          3.06997    2.99697
Coef-variation    0.90223    0.22021
3th-order moment  0.00193   -0.02728
4th-order moment  1.66491  156.17834
5th-order moment  0.01862   17.01957
6th-order moment  6.45525 5602.95438

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

R> plot(mod2d)
Ornstein-Uhlenbeck process and its integral

Ornstein-Uhlenbeck process and its integral

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   1.697779 14.5391
2   0.028943 18.2030
3   0.062434 18.4875
4  -1.093699 11.8186
5   0.377475 14.1720
6  -1.061190  4.5399
7   0.527083 12.9464
8   1.701165 20.9626
9   1.279216 18.5345
10  1.134661 18.8035
R> summary(out)
       x                y        
 Min.   :-2.913   Min.   :-5.39  
 1st Qu.:-0.382   1st Qu.:11.14  
 Median : 0.196   Median :14.60  
 Mean   : 0.192   Mean   :14.54  
 3rd Qu.: 0.772   3rd Qu.:17.95  
 Max.   : 3.652   Max.   :35.31  
R> cov(out)
        x       y
x 0.75085  2.0889
y 2.08894 25.9761

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 (10000 obs.);   Bandwidth 'bw' = 0.1229

       x                f(x)        
 Min.   :-3.2814   Min.   :0.00000  
 1st Qu.:-1.4559   1st Qu.:0.00363  
 Median : 0.3696   Median :0.05137  
 Mean   : 0.3696   Mean   :0.13682  
 3rd Qu.: 2.1951   3rd Qu.:0.25752  
 Max.   : 4.0206   Max.   :0.45864  

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

Data: y (10000 obs.);   Bandwidth 'bw' = 0.7247

       y               f(y)         
 Min.   :-7.563   Min.   :0.000001  
 1st Qu.: 3.698   1st Qu.:0.000214  
 Median :14.959   Median :0.007709  
 Mean   :14.959   Mean   :0.022178  
 3rd Qu.:26.220   3rd Qu.:0.042179  
 Max.   :37.482   Max.   :0.080483  
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 10000 obs.);

       x                 y              f(x,y)        
 Min.   :-2.9126   Min.   :-5.389   Min.   :0.000000  
 1st Qu.:-1.2715   1st Qu.: 4.785   1st Qu.:0.000001  
 Median : 0.3696   Median :14.959   Median :0.000156  
 Mean   : 0.3696   Mean   :14.959   Mean   :0.003668  
 3rd Qu.: 2.0107   3rd Qu.:25.133   3rd Qu.:0.002791  
 Max.   : 3.6519   Max.   :35.308   Max.   :0.041450  
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))
+ }

Return to snssde2d()

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: \[\begin{equation}\label{eq:12} \ddot{X}-\mu (1-X^{2}) \dot{X} + X = 0 \end{equation}\] 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)\) \[\begin{equation}\label{eq:13} \ddot{X}-\mu (1-X^{2}) \dot{X} + X = \xi_{t}, \end{equation}\] 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: \[\begin{equation}\label{eq:14} \begin{cases} \dot{X} = Y \\ \dot{Y} = \mu \left(1-X^{2}\right) Y - X + \xi_{t} \end{cases} \end{equation}\] The system can mathematically explain by a Stratonovitch equations: \[\begin{equation}\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} \end{equation}\]

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

Return to snssde2d()

snssde3d()

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

Itô form: \[\begin{equation}\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} \end{equation}\] Stratonovich form: \[\begin{equation}\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} \end{equation}\]

\(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:

Basic example

Assume that we want to describe the following SDE’s (3D) in Itô form: \[\begin{equation}\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} \end{equation}\]

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=10000)
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  = 10000.
    | 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.79598  0.87648  0.79364
R> moment(mod3d, at = s , center = TRUE , order = 2) ## variance
[1] 0.0096599 0.1032059 0.0103867
R> Median(mod3d, at = s)
[1] -0.80283  0.85439  0.80044
R> Mode(mod3d, at = s)
[1] -0.81892  0.77165  0.80357
R> quantile(mod3d , at = s)
$x
      0%      25%      50%      75%     100% 
-1.09932 -0.86489 -0.80283 -0.73250 -0.33409 

$y
       0%       25%       50%       75%      100% 
-0.053244  0.651313  0.854387  1.079257  2.315140 

$z
     0%     25%     50%     75%    100% 
0.29484 0.73033 0.80044 0.86553 1.10380 
R> kurtosis(mod3d , at = s)
[1] 3.2415 3.2116 3.4146
R> skewness(mod3d , at = s)
[1]  0.35784  0.37135 -0.44778
R> cv(mod3d , at = s )
[1] -0.12348  0.36655  0.12842
R> min(mod3d , at = s)
[1] -1.099317 -0.053244  0.294841
R> max(mod3d , at = s)
[1] -0.33409  2.31514  1.10380
R> moment(mod3d, at = s , center= TRUE , order = 4)
[1] 0.00030253 0.03421494 0.00036845
R> moment(mod3d, at = s , center= FALSE , order = 4)
[1] 0.43737 1.14324 0.43484

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.79598  0.87648  0.79364
Variance          0.00966  0.10322  0.01039
Median           -0.80283  0.85439  0.80044
Mode             -0.81892  0.77165  0.80357
First quartile   -0.86489  0.65131  0.73033
Third quartile   -0.73250  1.07926  0.86553
Minimum          -1.09932 -0.05324  0.29484
Maximum          -0.33409  2.31514  1.10380
Skewness          0.35784  0.37135 -0.44778
Kurtosis          3.24149  3.21159  3.41460
Coef-variation   -0.12348  0.36655  0.12842
3th-order moment  0.00034  0.01231 -0.00047
4th-order moment  0.00030  0.03421  0.00037
5th-order moment  0.00004  0.01270 -0.00005
6th-order moment  0.00002  0.02106  0.00003

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)
3D SDEs

3D SDEs

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.84120 0.96564 0.78370
2  -0.87580 1.16739 0.85986
3  -0.73101 1.04363 0.76730
4  -0.92454 1.28851 0.93054
5  -0.65252 0.21122 0.69777
6  -0.98157 1.95078 0.84848
7  -0.86326 1.18198 0.84539
8  -0.75253 1.03912 0.80682
9  -0.79263 0.59630 0.78578
10 -0.94422 1.24745 0.88755
R> summary(out)
       x                y                 z        
 Min.   :-1.099   Min.   :-0.0532   Min.   :0.295  
 1st Qu.:-0.865   1st Qu.: 0.6513   1st Qu.:0.730  
 Median :-0.803   Median : 0.8544   Median :0.800  
 Mean   :-0.796   Mean   : 0.8765   Mean   :0.794  
 3rd Qu.:-0.733   3rd Qu.: 1.0793   3rd Qu.:0.866  
 Max.   :-0.334   Max.   : 2.3151   Max.   :1.104  
R> cov(out)
           x         y          z
x  0.0096608 -0.018038 -0.0042498
y -0.0180380  0.103216  0.0193423
z -0.0042498  0.019342  0.0103877

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 (10000 obs.);   Bandwidth 'bw' = 0.01402

       x                 f(x)       
 Min.   :-1.14138   Min.   :0.0000  
 1st Qu.:-0.92904   1st Qu.:0.0190  
 Median :-0.71670   Median :0.3941  
 Mean   :-0.71670   Mean   :1.1762  
 3rd Qu.:-0.50436   3rd Qu.:2.2485  
 Max.   :-0.29203   Max.   :4.2014  

 Marginal density of Y(t-t0)|Y(t0)=-2 at time t = 1

Data: y (10000 obs.);   Bandwidth 'bw' = 0.04555

       y                 f(y)        
 Min.   :-0.18991   Min.   :0.00001  
 1st Qu.: 0.47052   1st Qu.:0.01225  
 Median : 1.13095   Median :0.14253  
 Mean   : 1.13095   Mean   :0.37817  
 3rd Qu.: 1.79137   3rd Qu.:0.72765  
 Max.   : 2.45180   Max.   :1.26758  

 Marginal density of Z(t-t0)|Z(t0)=-2 at time t = 1

Data: z (10000 obs.);   Bandwidth 'bw' = 0.01439

       z                f(z)       
 Min.   :0.25167   Min.   :0.0000  
 1st Qu.:0.47549   1st Qu.:0.0242  
 Median :0.69932   Median :0.3290  
 Mean   :0.69932   Mean   :1.1158  
 3rd Qu.:0.92315   3rd Qu.:2.0677  
 Max.   :1.14698   Max.   :4.0055  
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: \[\begin{equation}\label{eq20} dX_{t} = \mu W_{1,t} dt + \sigma W_{2,t} \circ dW_{3,t} \end{equation}\] 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=10000,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  = 10000.
    | 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.00009   -0.00421   0.00286
Variance              0.86357    1.01953   1.03221
Median                0.00062   -0.00279   0.00230
Mode                 -0.03670    0.07252  -0.00533
First quartile       -0.56728   -0.69611  -0.68297
Third quartile        0.57921    0.67839   0.67776
Minimum              -5.32251   -3.33941  -3.73809
Maximum               5.18314    3.52985   3.91080
Skewness             -0.06089   -0.00759   0.01373
Kurtosis              4.37318    2.90888   3.02486
Coef-variation   -10658.65427 -239.69710 354.91565
3th-order moment     -0.04886   -0.00781   0.01440
4th-order moment      3.26128    3.02364   3.22286
5th-order moment     -1.03009   -0.11215   0.14934
6th-order moment     29.64147   14.25180  16.66042

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)
Simulation of X_t

Simulation of \(X_t\)

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 (10000 obs.);   Bandwidth 'bw' = 0.122

       x                 y          
 Min.   :-5.6886   Min.   :0.00000  
 1st Qu.:-2.8792   1st Qu.:0.00062  
 Median :-0.0697   Median :0.00855  
 Mean   :-0.0697   Mean   :0.08890  
 3rd Qu.: 2.7398   3rd Qu.:0.10898  
 Max.   : 5.5493   Max.   :0.47644  
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 = ""))

Return to snssde3d()

Further reading

  1. snssdekd() & dsdekd() & rsdekd()- Monte-Carlo Simulation and Analysis of Stochastic Differential Equations.
  2. bridgesdekd() & dsdekd() & rsdekd() - Constructs and Analysis of Bridges Stochastic Differential Equations.
  3. fptsdekd() & dfptsdekd() - Monte-Carlo Simulation and Kernel Density Estimation of First passage time.
  4. MCM.sde() & MEM.sde() - Parallel Monte-Carlo and Moment Equations for SDEs.
  5. fitsde() - Parametric Estimation of 1-D Stochastic Differential Equation.

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 (2017). 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 (2017). Sim.DiffProc: Simulation of Diffusion Processes. R package version 3.8, 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}\).