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)$, $1\leq 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?
 A: 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} = \exp\Big(-\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.
A: Calculate the decomposed covariance matrix by incomplete Cholesky decomposition or any other matrix decomposition technique. Decomposed matrix should be TxM, where M is only a fraction of T.
http://en.wikipedia.org/wiki/Incomplete_Cholesky_factorization
