# R code for Ornstein-Uhlenbeck process

Can any one help me with some R code to run Ornstein-Uhlenbeck process?

The code of Euler Maruyama simulation method is pretty simple (nu is long run mean, lambda is mean reversion speed):

ornstein_uhlenbeck <- function(T,n,nu,lambda,sigma,x0){
dw  <- rnorm(n, 0, sqrt(T/n))
dt  <- T/n
x <- c(x0)
for (i in 2:(n+1)) {
x[i]  <-  x[i-1] + lambda*(nu-x[i-1])*dt + sigma*dw[i-1]
}
return(x);
}


Take a look at the sde package; specifically the dcOU and dsOU functions. You may also find some examples on the R-SIG-Finance mailing list, which would be in the results of a search on www.rseek.org.

You can also use the Sim.DiffProc package.

Have a look at this document:
Sim.DiffProc: A Package for Simulation of Diffusion Processes in R

See esp. chapter 2.1.2

There is even a Graphical User Interface (GUI) available for some functions:
http://cran.r-project.org/web/packages/Sim.DiffProcGUI/index.html

See chapter 4 in the above document for details.

• Ive used the SimDiffProc library to do the same but I feel the simulations are wrong. The Euler simulation process gives better results. – Mahesh Feb 23 '13 at 0:39
• @Mahesh: What makes you feel, that the simulations are wrong? – vonjd Dec 4 '17 at 11:03

The Euler method is simple but it gives an approximate distribution. The method implemented below gives an exact distribution of $X_{t_i}$ and exact conditional distributions $(X_{t_j} \mid X_{t_i})$.

rOU <- function(npaths, T, nsteps, x0, theta1, theta2, theta3){
dt <- T/nsteps
r <- theta1/theta2
s <- theta3*sqrt(-expm1(-2*theta2*dt)/2/theta2)
e <- exp(-theta2*dt)
out <- rbind(x0, matrix(NA_real_, nsteps, npaths))
for(i in 2:(nsteps+1)){
out[i,] <- rnorm(npaths, r+e*(out[i-1,]-r), s)
}
out
}


Let's check the covariance $Cov(X_{t_1}, X_{t_2})$:

 theta1 <- 1; theta2 <- 2; theta3 <- 3
> nsteps <- 10
> sims <- rOU(npaths=500000, T=1, nsteps=nsteps, x0=0,
+             theta1=theta1, theta2=theta2, theta3=theta3)
> # check covariance
> t1 <- 1/2; t2 <- 1
> cov(sims[1+nsteps*t1,], sims[1+nsteps*t2,]) # estimated
 0.713272
> theta3^2/2/theta2 *
+   (exp(-theta2*(t2-t1)) - exp(-theta2*(t2+t1))) #exact
 0.7157078