Simulating a Gaussian 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?

-

2
It's a stationary process (looks close to a simple version of an OU process). Is it uniformly sampled? – cardinal May 24 '12 at 13:53
The R package `mvtnorm` has `rmvnorm(n, mean, sigma)` where `sigma` is the covariance matrix; you'd have to construct the covariance matrix for your sampled / selected $t$s yourself, though. – jbowman May 24 '12 at 14:06

2
@jb Presumably $T$ is huge, otherwise the OP wouldn't be asking to avoid the matrix decomposition (which is implicit in `rmvnorm`). – whuber♦ May 24 '12 at 14:31

2 Answers

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.

-

Just a few more remarks. (1) To get a good idea of what the continuous time process looks like, $n$ and $T$ must be chosen so that $\Delta$ is small, say less than $0.1$. (2) The  inverse covariance (precision) matrix for the timeseries vector is tri-diagonal, as is its Cholesky root. – Yves May 24 '12 at 17:55
@Yves: Thanks for your comments. To be clear, the procedure I've outlined gives an exact realization of the continuous-time process sampled on the corresponding partition; in particular, there is no discretization error like there is in typical Euler-scheme approximation to more general SDEs. The inverse Cholesky, as shown by the construction in the answer has nonzero terms only on the diagonal and lower off-diagonal, so it's a little simpler than tridiagonal. – cardinal May 24 '12 at 18:13

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

-

2
Can you give an explicit form of the Cholesky decomposition here? I think that the answer by cardinal achieves just that, if you think about it, by expressing $X_i$ as a function of the history. – StasK May 25 '12 at 2:54
The algorithm is a little too long to summarize. You can find an excellent description here: Kernel ICA, page 20. Note that this algorithm is incomplete, meaning it doesn't calculate the entire decomposition but rather an approximation (hence it is much faster). I have published code for this algorithm in the KMBOX toolbox, you can download it here: km_kernel_icd‌​. – Steven Jun 6 '12 at 20:45
default
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 29, "mathjax_display_tex": 5, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0}, "config": {"markdown_headings": false, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "english_confidence": 0.927247941493988, "perplexity_flag": "head"}