Código R para el proceso Ornstein-Uhlenbeck

Question

¿Puede alguien ayudarme con algún código R para ejecutar el proceso Ornstein-Uhlenbeck?

Answer 1

El código de Método de simulación Euler Maruyama es bastante simple ( nu es la media de largo plazo, lambda es la velocidad de reversión media):

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);
}

Answer 2

Echa un vistazo a la sde en concreto, el paquete dcOU y dsOU funciones. También puede encontrar algunos ejemplos en la lista de correo de R-SIG-Finance, que estaría en los resultados de una búsqueda en www.rseek.org .

Answer 3

También puede utilizar el Sim.DiffProc paquete.

Echa un vistazo a este documento:
Sim.DiffProc: Un paquete para la simulación de procesos de difusión en R

Véase especialmente el capítulo 2.1.2

Incluso hay una interfaz gráfica de usuario (GUI) disponible para algunas funciones:
http://cran.r-project.org/web/packages/Sim.DiffProcGUI/index.html

Véase el capítulo 4 del documento anterior para más detalles.

Answer 4

El método de Euler es sencillo pero da una distribución aproximada. El método implementado a continuación da una distribución exacta de $X_{t_i}$ y las distribuciones condicionales exactas $(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
}

Comprobemos la covarianza $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
[1] 0.713272
> theta3^2/2/theta2 * 
+   (exp(-theta2*(t2-t1)) - exp(-theta2*(t2+t1))) #exact
[1] 0.7157078