This is probably a very straight forward question but I want to verify how I should sample from an AR(1) process in R using just the rnorm() function in R (or any similar function in another language). Say I want to sample 100 samples from the following AR(1) model:

$$ X_t = \rho X_{t-1} + \epsilon_t, \qquad \epsilon_t \sim N(0, \sigma^2_\epsilon). $$

Then the joint distribution of the vector $(X_1,\dots,X_{100})$ is multivariate Gaussian with zero mean and covariance matrix $\Sigma$, where

$$ BB^T :=\Sigma = \frac{\sigma^2_\epsilon}{1-\rho^2} \begin{bmatrix} 1 & \rho & \rho^2 & \rho^3 & \rho^4&\cdots\\ \rho & 1 &\rho & \rho^2 & \rho^3 &\cdots\\ \rho^2 &\rho & 1 &\rho & \rho^2&\cdots\\ \rho^3 &\rho^2 &\rho & 1 & \rho&\cdots\\ \vdots &\vdots &\vdots& \vdots &\ddots& \vdots\\ \end{bmatrix} $$ So I should be able to sample in R using the following:

X = matrix(rnorm(100, 0, sigmaeps), 1, 100) %*% B

Is this approach correct?


1 Answer 1


Your general approach is fine, but obviously you will need to find the matrix $B$ from your stated form for $\Sigma$. (Be careful here --- unless you use the symmetric square root as $B$ you will need to amend your code to pre-multiply by $B$ instead of post-multiplying.)

In any case, if you want to save yourself some time you can use the rGARMA function in the ts.extend package where this stuff is already programmed. This function will produce random vectors from any stationary Gaussian ARMA model, and it can also generate from the conditional distribution if some of the time-series values are given. Here is an example:

#Load the package

#Produce n = 16 time-series vectors from AR(1) process with length m = 30
m        <- 30
AR       <- 0.7
MEAN     <- 4
SERIES   <- rGARMA(n = 16, m = m, ar = AR, mean = MEAN, errorvar = ERRORVAR)

#Plot the series

enter image description here

  • $\begingroup$ do you mind explaining a bit more about pre and post multiplying with B, I was planning to take the Cholesky decomposition: $\Sigma = RR^T$ where R is the upper triangular $\endgroup$ Dec 29, 2020 at 2:45
  • $\begingroup$ The variance rules for linear transforms are $\mathbb{V}(B \epsilon) = B \mathbb{V}(\epsilon) B^\text{T}$ and $\mathbb{V}(\epsilon B) = B^\text{T} \mathbb{V}(\epsilon) B$, so the variance you will get with your formula (using post-multiplication) is $B^\text{T} B \neq B B^\text{T}$. Unless I am mistaken, the Cholesky decomposition doesn't allow you to swap the order of the parts. $\endgroup$
    – Ben
    Dec 29, 2020 at 4:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.