Sampling from an AR(1) process using normal samples

Question

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?

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
library(ts.extend)

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

#Plot the series
plot(SERIES)