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?