Monte-Carlo Simulation and Analysis of Stochastic Differential Equations

A.C. Guidoum1 and K. Boukhetala2

2020-03-27

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=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 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.105
R> moment(mod1, at = s , center = TRUE , order = 2) ## variance
[1] 32.835
R> Median(mod1, at = s)
[1] 10.002
R> Mode(mod1, at =s)
[1] 7.8095
R> quantile(mod1 , at = s)
     0%     25%     50%     75%    100% 
 1.9254  7.1387 10.0024 13.3241 51.1648 
R> kurtosis(mod1 , at = s)
[1] 8.2632
R> skewness(mod1 , at = s)
[1] 1.7425
R> cv(mod1 , at = s )
[1] 0.51627
R> min(mod1 , at = s)
[1] 1.9254
R> max(mod1 , at = s)
[1] 51.165
R> moment(mod1, at = s , center= TRUE , order = 4)
[1] 8926.5
R> moment(mod1, at = s , center= FALSE , order = 4)
[1] 63011

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.10464
Variance              32.86746
Median                10.00239
Mode                   7.80952
First quartile         7.13866
Third quartile        13.32414
Minimum                1.92537
Maximum               51.16476
Skewness               1.74255
Kurtosis               8.26317
Coef-variation         0.51627
3th-order moment     328.34764
4th-order moment    8926.45431
5th-order moment  233172.08753
6th-order moment 7316857.36310
R> summary(mod2, at = 1)

    Monte-Carlo Statistics for X(t) at time t = 1
                             
Mean                  9.85097
Variance             24.27731
Median                8.86305
Mode                  7.43570
First quartile        6.30123
Third quartile       12.35974
Minimum               2.20022
Maximum              31.71660
Skewness              1.26052
Kurtosis              5.09477
Coef-variation        0.50017
3th-order moment    150.78260
4th-order moment   3002.79303
5th-order moment  44410.08970
6th-order moment 815890.65508

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.1029  5.6148 11.1593  6.9677 18.0892  4.6807  9.5430  7.5471
 [9]  7.1909  6.5665
R> head(x2,n=10)
 [1]  3.8854  8.4079  3.7702 10.0380  6.4972 13.6009  7.8279 11.8803
 [9] 11.3957 14.8763
R> summary(data.frame(x1,x2))
       x1              x2       
 Min.   : 1.92   Min.   : 2.20  
 1st Qu.: 7.14   1st Qu.: 6.30  
 Median :10.00   Median : 8.86  
 Mean   :11.11   Mean   : 9.85  
 3rd Qu.:13.32   3rd Qu.:12.36  
 Max.   :51.16   Max.   :31.72  

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 \(1000\) trajectories of \((X_{t},Y_{t})\), with integration step size \(\Delta t = 0.01\), and using second Milstein method.

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:

[1]  0.95195 12.16559
[1] 0.72754 7.18564
[1]  0.92433 12.14727
[1]  0.8266 12.0993
$x
      0%      25%      50%      75%     100% 
-2.15377  0.35614  0.92433  1.53354  3.52599 

$y
     0%     25%     50%     75%    100% 
 3.4447 10.2635 12.1473 13.9465 21.3874 
[1] 2.9428 2.9313
[1] 0.13690 0.14656
[1] 0.89647 0.22045
[1] -2.1538  3.4447
[1]  3.526 21.387
[1]   1.5608 151.6568
[1]     6.6618 28574.6041

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


    Monte-Carlo Statistics for (X(t),Y(t)) at time t = 5
                        X          Y
Mean              0.95195   12.16559
Variance          0.72827    7.19283
Median            0.92433   12.14727
Mode              0.82660   12.09929
First quartile    0.35614   10.26354
Third quartile    1.53354   13.94653
Minimum          -2.15377    3.44466
Maximum           3.52599   21.38740
Skewness          0.13690    0.14656
Kurtosis          2.94280    2.93131
Coef-variation    0.89647    0.22045
3th-order moment  0.08508    2.82727
4th-order moment  1.56080  151.65678
5th-order moment  0.23980  198.90502
6th-order moment  5.50616 5308.99475

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

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

           x       y
1   0.923572 14.2278
2   1.072368 11.5379
3   0.119491 20.1443
4   0.138435 13.8444
5  -1.209009  4.8711
6  -1.541834 11.1854
7   1.110304 13.9646
8  -0.786921 15.8038
9  -0.545338 13.0437
10 -0.053273 12.9980
       x                y        
 Min.   :-2.316   Min.   :-1.42  
 1st Qu.:-0.370   1st Qu.:11.17  
 Median : 0.165   Median :14.54  
 Mean   : 0.185   Mean   :14.62  
 3rd Qu.: 0.767   3rd Qu.:17.96  
 Max.   : 3.222   Max.   :33.00  
        x       y
x 0.72687  1.9337
y 1.93372 25.7796

Figure 4, show simulation results for moments of system :

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


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

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

       x                f(x)        
 Min.   :-2.8915   Min.   :0.00002  
 1st Qu.:-1.2192   1st Qu.:0.00924  
 Median : 0.4531   Median :0.07572  
 Mean   : 0.4531   Mean   :0.14935  
 3rd Qu.: 2.1254   3rd Qu.:0.28944  
 Max.   : 3.7977   Max.   :0.48332  

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

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

       y               f(y)         
 Min.   :-4.856   Min.   :0.000004  
 1st Qu.: 5.468   1st Qu.:0.001104  
 Median :15.791   Median :0.010757  
 Mean   :15.791   Mean   :0.024193  
 3rd Qu.:26.115   3rd Qu.:0.047181  
 Max.   :36.438   Max.   :0.074897  

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.


 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.   :-2.3158   Min.   :-1.420   Min.   :0.000000  
 1st Qu.:-0.9313   1st Qu.: 7.186   1st Qu.:0.000078  
 Median : 0.4531   Median :15.791   Median :0.000852  
 Mean   : 0.4531   Mean   :15.791   Mean   :0.005122  
 3rd Qu.: 1.8375   3rd Qu.:24.396   3rd Qu.:0.005889  
 Max.   : 3.2220   Max.   :33.002   Max.   :0.040300  

A \(3\)D plot of the transition density at \(t=10\) obtained with:

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

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.

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.

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

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:

[1] -0.79493  0.86777  0.79348
[1] 0.0093895 0.1039466 0.0104875
[1] -0.80094  0.85686  0.79924
[1] -0.80855  0.86337  0.81022
$x
      0%      25%      50%      75%     100% 
-1.06097 -0.86137 -0.80094 -0.73601 -0.33834 

$y
       0%       25%       50%       75%      100% 
-0.030089  0.648630  0.856856  1.076947  2.175117 

$z
     0%     25%     50%     75%    100% 
0.37315 0.72625 0.79924 0.86569 1.04056 
[1] 3.4267 3.5087 3.2125
[1]  0.41420  0.40985 -0.37923
[1] -0.12196  0.37172  0.12913
[1] -1.060969 -0.030089  0.373149
[1] -0.33834  2.17512  1.04056
[1] 0.00030271 0.03798713 0.00035404
[1] 0.43402 1.12242 0.43508

The summary of the results of mod3d at time \(t=1\) of class snssde3d() is given by:


  Monte-Carlo Statistics for (X(t),Y(t),Z(t)) at time t = 1
                        X        Y        Z
Mean             -0.79493  0.86777  0.79348
Variance          0.00940  0.10405  0.01050
Median           -0.80094  0.85686  0.79924
Mode             -0.80855  0.86337  0.81022
First quartile   -0.86137  0.64863  0.72625
Third quartile   -0.73601  1.07695  0.86569
Minimum          -1.06097 -0.03009  0.37315
Maximum          -0.33834  2.17512  1.04056
Skewness          0.41420  0.40985 -0.37923
Kurtosis          3.42670  3.50871  3.21249
Coef-variation   -0.12196  0.37172  0.12913
3th-order moment  0.00038  0.01376 -0.00041
4th-order moment  0.00030  0.03799  0.00035
5th-order moment  0.00004  0.01708 -0.00005
6th-order moment  0.00002  0.02735  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.

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

          x       y       z
1  -0.90979 0.98650 0.76635
2  -0.80140 0.49268 0.77932
3  -0.76156 0.77070 0.93226
4  -0.65886 0.89811 0.93844
5  -0.88491 1.01199 0.85718
6  -0.63544 0.66637 0.75138
7  -0.90259 0.94879 0.88718
8  -0.62886 0.34688 0.62979
9  -0.94196 0.89948 0.80507
10 -0.58044 0.51197 0.60664
       x                y                 z        
 Min.   :-1.061   Min.   :-0.0301   Min.   :0.373  
 1st Qu.:-0.861   1st Qu.: 0.6486   1st Qu.:0.726  
 Median :-0.801   Median : 0.8569   Median :0.799  
 Mean   :-0.795   Mean   : 0.8678   Mean   :0.793  
 3rd Qu.:-0.736   3rd Qu.: 1.0770   3rd Qu.:0.866  
 Max.   :-0.338   Max.   : 2.1751   Max.   :1.041  
           x         y          z
x  0.0093989 -0.017872 -0.0043331
y -0.0178716  0.104051  0.0203124
z -0.0043331  0.020312  0.0104980

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.


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

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

       x                 f(x)       
 Min.   :-1.12442   Min.   :0.0002  
 1st Qu.:-0.91204   1st Qu.:0.0241  
 Median :-0.69966   Median :0.3668  
 Mean   :-0.69966   Mean   :1.1760  
 3rd Qu.:-0.48728   3rd Qu.:2.1560  
 Max.   :-0.27490   Max.   :4.2248  

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

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

       y                 f(y)        
 Min.   :-0.24687   Min.   :0.00007  
 1st Qu.: 0.41282   1st Qu.:0.01523  
 Median : 1.07251   Median :0.16975  
 Mean   : 1.07251   Mean   :0.37859  
 3rd Qu.: 1.73221   3rd Qu.:0.71164  
 Max.   : 2.39190   Max.   :1.26484  

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

Data: z (1000 obs.);    Bandwidth 'bw' = 0.02316

       z                f(z)       
 Min.   :0.30366   Min.   :0.0003  
 1st Qu.:0.50526   1st Qu.:0.0343  
 Median :0.70686   Median :0.6145  
 Mean   :0.70686   Mean   :1.2389  
 3rd Qu.:0.90845   3rd Qu.:2.4851  
 Max.   :1.11005   Max.   :3.7382  

For an approximate joint transition density for \((X_t,Y_t,Z_t)\) (for more details, see package sm or ks.)

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.

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:

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

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.

  Monte-Carlo Statistics for (X(t),Y(t),Z(t)) at time t = 1
                         X         Y        Z
Mean               0.00210   0.00711  0.02497
Variance           0.75810   1.04206  0.98254
Median            -0.00023  -0.02394  0.02131
Mode              -0.02666  -0.21236 -0.01505
First quartile    -0.51844  -0.72702 -0.63515
Third quartile     0.49242   0.75232  0.71023
Minimum           -4.00577  -3.02569 -3.20162
Maximum            3.53355   3.69549  2.78420
Skewness          -0.05380   0.09120  0.00749
Kurtosis           4.32863   2.78478  2.92680
Coef-variation   414.81602 143.58877 39.69596
3th-order moment  -0.03551   0.09701  0.00729
4th-order moment   2.48771   3.02396  2.82549
5th-order moment  -0.67727   0.83774 -0.31149
6th-order moment  16.81061  14.67331 12.82519

the following code produces the result in Figure 12.

 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.


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

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

       x                y          
 Min.   :-4.517   Min.   :0.00003  
 1st Qu.:-2.377   1st Qu.:0.00299  
 Median :-0.236   Median :0.02468  
 Mean   :-0.236   Mean   :0.11667  
 3rd Qu.: 1.905   3rd Qu.:0.17103  
 Max.   : 4.045   Max.   :0.51682  

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})\):

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. TEX.sde() - Converting Sim.DiffProc Objects to LaTeX.
  6. 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 (2020). Performing Parallel Monte Carlo and Moment Equations Methods for Ito and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc. Accept Submission to Journal of Statistical Software.

  3. Guidoum AC, Boukhetala K (2020). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.5, 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 ()

  2. Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, E-mail ()

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