R: Autocorrelation in a matrix

Question

i have do to a monte carlo approach for AR(1) time series. I have to generate 10,000 time series of length 100 and afterwards i have to get the first step autocorrelation rho_1 for every time series. My problem is that i just get NA values for the autocorrelation and the calculation takes way to much time. I have no problem with computing the AR(1) time series. Thank you for your help :)

gen_ar <- function(a,b,length,start)
{
 z<-rep(0,length)                                 
 e<-rnorm(n=length,sd=1)  
 z[1]<-start
  for (i in 2:length)
  {
    z[i]<-a+b*z[i-1]+e[i]
  }
 z
}

mc <- matrix(c(rep(0,10000000)),nrow=10000)
 for (i in 1:10000)
 {
  mc[i,] <- gen_ar(0.99,1,100,0)
 }

 ac <- matrix(c(rep(0,10000)),nrow=1)
  for (i in 1:10000){
   for (j in 1:99){
    ac[i] <- cor(mc[i,j],mc[i,j+1])
   }
  }

Answer 1

Statistics aside, I think this achieves your goals, and I don't get NA's. I changed the way it was done b/c you said it was going slow.

mc <- matrix(rep(NA,1E5), nrow=100)
for(i in seq_len(100)){
    mc[,i] <- arima.sim(model=list(ar=0.99), n=100, sd=1) + 1
}

myAR <- function(x){
    cor(x[-1], x[-length(x)])
}

answer <- apply(mc, 2, myAR)

I skipped the last set of nested for loops and replaced them with apply(). It seems easier to read, and is likely faster. Also, to use apply(), I created a function called myAR, which carries out the same calculation that cor() did in your for() loops.

Now, there are a couple of statistical adjustments that I made. Primarily, these were in the simulation step.

First, your simulated AR(1) process has a coefficient that is equal to 1, which seems odd to me (this would not be stationary, and arima.sim() won't even let you simulate this type of process).

Moreover, your "a" parameter adds 1 to the time series at each time step. In other words, your time series is monotonically increasing from 1 to 100 because the coefficient is equal to 1. This too would make your time series nonstationary, and with such a strong positive slope the cor() function would likely return 1 as the estimated correlation, regardless of the value of the simulated AR coefficient. I assume that you wanted the long-term mean to hover near 1, so the 1 is simply added to the entire time series after it is simulated, not iteratively at each time step.

Assuming that you did want to generate a nonstationary time series by adding some constant (a) at each time step, you could do the following:

myInnov <- function(N=100, a=1, SD=1) {a + rnorm(n=N, sd=SD)}
mc2 <- matrix(rep(NA,1E7), nrow=100)
for(i in seq_len(1E5)){
    mc2[,i] <- arima.sim(model=list(ar=0.99), n=100, innov=myInnov(a=1, N=100, SD=1)) + 1
}

I hope that this helps.