# simulate autoregressive data that is also multivariate normal

I am trying to simulate data that is correlated to varying degrees. However, the data itself will have a degree of autocorrelation as well. I can get the first part of the problem with mvrnorm function using this code:

library(MASS)

set.seed(1531488844)

# Generate random data values

r <- 0.95
xrow = 10
ySigma <- matrix(
c(1  , r  , r^2, r^3, r^4, r^5, r^6, r^7, r^8, r^9,
r  , 1  , r  , r^2, r^3, r^4, r^5, r^6, r^7, r^8,
r^2, r  , 1  , r  , r^2, r^3, r^4, r^5, r^6, r^7,
r^3, r^2, r  , 1  , r  , r^2, r^3, r^4, r^5, r^6,
r^4, r^3, r^2, r  , 1  , r  , r^2, r^3, r^4, r^5,
r^5, r^4, r^3, r^2, r  , 1  , r  , r^2, r^3, r^4,
r^6, r^5, r^4, r^3, r^2, r  , 1  , r  , r^2, r^3,
r^7, r^6, r^5, r^4, r^3, r^2, r  , 1  , r  , r^2,
r^8, r^7, r^6, r^5, r^4, r^3, r^2, r  , 1  , r  ,
r^9, r^8, r^7, r^6, r^5, r^4, r^3, r^2, r  , 1  ),
nrow = xrow
)

rawvars <- mvrnorm(n=10000, mu=rep(0,xrow), Sigma=ySigma)

pairs(rawvars[1:100,])


BUT the problem is that each individual column is white noise and looks like this:

What I would like to have is each column of data follow an autoregressive process so it looks like this:

I know that mvrnorm can't create autoregressive data - but I was thinking there may be a way to accomplish what I am after by combining the two functions.