# Simple Gaussian Process Simulation on R

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 <- 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()
`````` 