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: In these equations, μ and σ 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: 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:

We simulate a flow of 1000 trajectories of (X_t, Y_t), with integration step size Δt = 0.01, and using second Milstein method.

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
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 summary of the results of mod2d at time t = 10 of class snssde2d() is given by:

R> summary(mod2d, at = 10)

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

R> ## in time
R> plot(mod2d)
R> ## in plane (O,X,Y)
R> plot2d(mod2d,type="n") 
R> points2d(mod2d,col=rgb(0,100,0,50,maxColorValue=255), pch=16)

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=3)

        x       y
1 0.32480 12.2224
2 0.91370  9.3092
3 0.49834 14.4406

The density of X_t and Y_t at time t = 10 are reported using dsde2d() function, see e.g. Figure 4: the marginal density of X_t and Y_t at time t = 10. For 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, see:

R> ## the marginal density
R> denM <- dsde2d(mod2d,pdf="M",at =10)
R> plot(denM, main="Marginal Density")
R> ## the Joint density
R> denJ <- dsde2d(mod2d, pdf="J", n=100,at =10)
R> plot(denJ,display="contour",main="Bivariate Transition Density at time t=10")

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

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

Marginal and Joint density at time t=10

We approximate the bivariate transition density over the set transition horizons t ∈ [1, 10] by Δ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 ẋ = 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: where x is the position coordinate (which is a function of the time t), and μ 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 ξ_t, with delta-type correlation function E(ξ_tξ_t + h) = 2σδ(h) where μ > 0 . It’s solution cannot be obtained in terms of elementary functions, even in the phase plane. The White noise ξ_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 ẋ = y, from physical equation we get the above system: The system can mathematically explain by a Stratonovitch equations:

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

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
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")

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

R> plot(mod2d,ylim=c(-8,8))   ## back in time
R> plot2d(mod2d)              ## in plane (O,X,Y)

Return to snssde2d()

The Heston Model

Consider a system of stochastic differential equations:

Conditions to ensure positiveness of the volatility process are that 2νθ > σ², and the two Brownian motions (B_1, t, B_2, t) are correlated. Σ to describe the correlation structure, for example: $$ \Sigma= \begin{pmatrix} 1 & 0.3 \\ 0.3 & 2 \end{pmatrix} $$

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
R> mu = 1.2; sigma=0.1;nu=2;theta=0.5
R> fx <- expression( mu*x ,nu*(theta-y)) 
R> gx <- expression( x*sqrt(y) ,sigma*sqrt(y))
R> Sigma <- matrix(c(1,0.3,0.3,2),nrow=2,ncol=2) # correlation matrix
R> HM <- snssde2d(drift=fx,diffusion=gx,Dt=0.001,x0=c(100,1),corr=Sigma,M=1000)
R> HM

Itô Sde 2D:
    | dX(t) = mu * X(t) * dt + X(t) * sqrt(Y(t)) * dB1(t)
    | dY(t) = nu * (theta - Y(t)) * dt + sigma * sqrt(Y(t)) * dB2(t)
    | Correlation structure:                    
             1.0 0.3
             0.3 2.0
Method:
    | Euler scheme with order 0.5
Summary:
    | Size of process   | N  = 1001.
    | Number of simulation  | M  = 1000.
    | Initial values    | (x0,y0) = (100,1).
    | Time of process   | t in [0,1].
    | Discretization    | Dt = 0.001.

Hence we can just make use of the rsde2d() function to build our random number for (X_t, Y_t) at time t = 1.

R> out <- rsde2d(object = HM, at = 1) 
R> head(out,n=3)

       x       y
1 181.42 0.56526
2 137.85 0.52612
3 243.57 0.62700

The density of X_t and Y_t at time t = 1 are reported using dsde2d() function. See:

R> denJ <- dsde2d(HM,pdf="J",at =1,lims=c(-100,900,0.4,0.75))
R> plot(denJ,display="contour",main="Bivariate Transition Density at time t=10")
R> plot(denJ,display="persp",main="Bivariate Transition Density at time t=10")

Return to snssde2d()

Basic example

Assume that we want to describe the following SDE’s (3D) in Ito form: with (W_1, t, W_2, t, W_3, t) are three indpendant standard Wiener process.

We simulate a flow of 1000 trajectories, with integration step size Δt = 0.001.

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
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)
R> moment(mod3d, at = s , center = TRUE , order = 2) ## variance
R> Median(mod3d, at = s)
R> Mode(mod3d, at = s)
R> quantile(mod3d , at = s)
R> kurtosis(mod3d , at = s)
R> skewness(mod3d , at = s)
R> cv(mod3d , at = s )
R> min(mod3d , at = s)
R> max(mod3d , at = s)
R> moment(mod3d, at = s , center= TRUE , order = 4)
R> moment(mod3d, at = s , center= FALSE , order = 4)

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

R> summary(mod3d, at = t)

For plotting (back in time) using the command plot, and plot3D in space the results of the simulation are shown in Figure 7.

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 = 1) 
R> head(out,n=3)

         x       y       z
1 -0.73929 0.65851 0.64650
2 -0.63795 0.41731 0.79531
3 -0.80249 1.16501 0.83721

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 8.

R> den <- dsde3d(mod3d,pdf="M",at =1)
R> plot(den, main="Marginal Density") 

Marginal density of X_t, Y_t and Z_t at time t = 1

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): with initial conditions (X₀, Y₀, Z₀) = (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> set.seed(1234, kind = "L'Ecuyer-CMRG")
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))

The results of simulation (3D) are shown in Figure 9:

R> plot3D(mod3d,display="persp",col="blue")

Attractive model for 3D diffusion processes

Return to snssde3d()

Transformation of an SDE one-dimensional

Next is an example of one-dimensional SDE driven by three correlated Wiener process (B_1, t,B_2, t,B_3, t), as follows: with: $$ \Sigma= \begin{pmatrix} 1 & 0.2 &0.5\\ 0.2 & 1 & -0.7 \\ 0.5 &-0.7&1 \end{pmatrix} $$ To simulate the solution of the process X_t, we make a transformation to a system of three equations as follows: run by calling the function snssde3d() to produce a simulation of the solution, with μ = 1 and σ = 1.

R> set.seed(1234, kind = "L'Ecuyer-CMRG")
R> fx <- expression(y,0,0) 
R> gx <- expression(z,1,1)
R> Sigma <-matrix(c(1,0.2,0.5,0.2,1,-0.7,0.5,-0.7,1),nrow=3,ncol=3)
R> modtra <- snssde3d(drift=fx,diffusion=gx,M=1000,corr=Sigma)
R> modtra

Itô Sde 3D:
    | dX(t) = Y(t) * dt + Z(t) * dB1(t)
    | dY(t) = 0 * dt + 1 * dB2(t)
    | dZ(t) = 0 * dt + 1 * dB3(t)
    | Correlation structure:                          
             1.0  0.2  0.5
             0.2  1.0 -0.7
             0.5 -0.7  1.0
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.

The histogram and kernel density of X_t at time t = 1 are reported using rsde3d() function, and we calculate emprical variance-covariance matrix C(s, t) = Cov(X_s, X_t), see e.g. Figure 10.

R> X <- rsde3d(modtra,at=1)$x
R> MASS::truehist(X,xlab = expression(X[t==1]));box()
R> lines(density(X),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)
R> ## Cov-Matrix
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()