# Solved – Simulating a Gaussian (Ornstein Uhlenbeck) process with an exponentially decaying covariance function

I'm trying to generate many draws (i.e., realizations) of
a Gaussian process \$e_i(t)\$, \$1leq t leq T\$ with mean 0 and covariance
function \$gamma(s,t)=exp(-|t-s|)\$.

Is there an efficient way to do this that wouldn't involve computing the square
root of a \$T times T\$ covariance matrix? Alternatively can
anyone recommend an `R` package to do this?

Contents

Yes. There is a very efficient (linear time) algorithm, and the intuition for it comes directly from the uniformly-sampled case.

Suppose we have a partition of \$[0,T]\$ such that \$0=t_0 < t_1 < t_2 < cdots < t_n = T\$.

Uniformly sampled case

In this case we have \$t_i = i Delta\$ where \$Delta = T/n\$. Let \$X_i := X(t_i)\$ denote the value of the discretely sampled process at time \$t_i\$.

It is easy to see that the \$X_i\$ form an AR(1) process with correlation \$rho = exp(-Delta)\$. Hence, we can generate a sample path \${X_t}\$ for the partition as follows \$\$ X_{i+1} = rho X_i + sqrt{1-rho^2} Z_{i+1} >, \$\$ where \$Z_i\$ are iid \$mathcal N(0,1)\$ and \$X_0 = Z_0\$.

General case

We might then imagine that it could be possible to do this for a general partition. In particular, let \$Delta_i = t_{i+1} – t_i\$ and \$rho_i = exp(-Delta_i)\$. We have that \$\$ gamma(t_i,t_{i+1}) = rho_i >, \$\$ and so we might guess that \$\$ X_{i+1} = rho_i X_i + sqrt{1-rho_i^2} Z_{i+1} >. \$\$

Indeed, \$mathbb E X_{i+1} X_i = rho_i\$ and so we at least have the correlation with the neighboring term correct.

The result now follows by telescoping via the tower property of conditional expectation. Namely, \$\$ newcommand{e}{mathbb E} e X_i X_{i-ell} = e( e(X_i X_{i-ell} mid X_{i-1} )) = rho_{i-1} mathbb E X_{i-1} X_{i-ell} = cdots = prod_{k=1}^ell rho_{i-k} >, \$\$ and the product telescopes in the following way \$\$ prod_{k=1}^ell rho_{i-k} = expBig(-sum_{k=1}^ell Delta_{i-k}Big) = exp(t_{i-ell} – t_i) = gamma(t_{i-ell},t_i) >. \$\$

This proves the result. Hence the process can be generated on an arbitrary partition from a sequence of iid \$mathcal N(0,1)\$ random variables in \$O(n)\$ time where \$n\$ is the size of the partition.

NB: This is an exact sampling technique in that it provides a sampled version of the desired process with the exactly correct finite-dimensional distributions. This is in contrast to Euler (and other) discretization schemes for more general SDEs, which incur a bias due to the approximation via discretization.

Rate this post