CODEWITHSUNDEEP
rstochasticstochastic-process
Ticktock
Ticktock

Reputation: 3

Simple Gaussian Process Simulation on R

How to simulate a Gaussian process X(t), t = 1, . . . , 200, with mean value function m(t) = 0 and covariance function r(h) = Cov(t, t + h) = exp(-|h|). I Know that this process is sometimes referred to as the Ornstein-Uhlenbeck process but how to plot the simulated process.

Thanks in anticipation

Upvotes: 0

Views: 1243

Answers (2)

Abdessabour Mtk
Abdessabour Mtk

Reputation: 3888

Following the Wikipedia definition an Ornstein–Uhlenbeck process is defined by the following stochastic differential equation: Where is a wiener process, one of it's properties is that it has Gaussian increments i.e

the aforementioned equation can be discretized in the following fashion: Where And due to the Gaussian increments property of the Wiener process we have that . That means that the values of the increments can be generated using sqrt(dt)*rnorm(1)

I coded the following function in R that takes the time vector, the mean of the process, the standard deviation and the value of theta.

simulate <- function(t, mean=0, std=1, x0=mean, theta=1, number.of.points=length(t)){
  # calculate time differences
  dt <- diff(t)
  X <- vector("numeric", length=number.of.points)
  X[1] <- x0
  for(i in 1:(number.of.points-1)){
    X[i+1] <- X[i] + theta * (mean-X[i])*dt[i] + std * sqrt(dt[i])* rnorm(1)
  }
  data.frame(x=t, y=X)
}
simulate(t=1:200) %>%  ggplot(aes(x,y)) + geom_line()

Another implementation Using purrr

simulate <- function(t, mean=0, sd=1, theta=1, number.of.points=length(t)){
  stopifnot(!missing(t) | !missing(number.of.points))
  if(missing(t)){
    t <- 1:number.of.points
  } 

  unlist(purrr::accumulate2(vector("numeric", length=number.of.points-1), diff(t), function(x, o, y) {
    x + theta*(mean - x)* y + sqrt(y)*rnorm(1)
  }, .init=x0), use.names=F) -> X
  
  data.frame(x=t, y=X)
}
simulate(number.of.points=200) %>%  ggplot(aes(x,y)) + geom_line()

Upvotes: 1

tester
tester

Reputation: 1692

Using the function from here: https://quant.stackexchange.com/questions/1260/r-code-for-ornstein-uhlenbeck-process

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

test <- ornstein_uhlenbeck(200, 200, 0, 0.8, 1, 0)
plot(x = seq_along(test), y = test, type = 'l')

(Note that it gives you only an approximate distribution as noted in one of the answers to the question in the link.)

I assumed T = 200, n = 200, nu = 0 (as you mentioned), mean reversion parameter of 0.8, sigma of 1 and the process starts at 0.

Upvotes: 0

Related Questions