Generating correlation matrix for AR(1) model

Assume that we are in the time series data setting, where we have data at equally-spaced times which we denote by random variables . The AR(1) model, commonly used in econometrics, assumes that the correlation between and is , where is some parameter that usually has to be estimated.

If we were writing out the full correlation matrix for consecutive data points , it would look something like this:

(Side note: This is an example of a correlation matrix which has Toeplitz structure.)

Given , how can we generate this matrix quickly in R? The function below is my (current) best attempt:

ar1_cor <- function(n, rho) {
exponent <- abs(matrix(1:n - 1, nrow = n, ncol = n, byrow = TRUE) - 
    (1:n - 1))
rho^exponent
}

In the function above, n is the number of rows in the desired correlation matrix (which is the same as the number of columns), and rho is the parameter. The function makes use of the fact that when subtracting a vector from a matrix, R automatically recycles the vector to have the same number of elements as the matrix, and it does so in a column-wise fashion.

Here is an example of how the function can be used:

ar1_cor(4, 0.9)
#       [,1] [,2] [,3]  [,4]
# [1,] 1.000 0.90 0.81 0.729
# [2,] 0.900 1.00 0.90 0.810
# [3,] 0.810 0.90 1.00 0.900
# [4,] 0.729 0.81 0.90 1.000

Such a function might be useful when trying to generate data that has such a correlation structure. For example, it could be passed as the Sigma parameter for MASS::mvrnorm(), which generates samples from a multivariate normal distribution.

Can you think of other ways to generate this matrix?