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

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.