I'm working on a time series problem where the spacing between observations is usually 12 or 24 hours, but this is not guaranteed. I'd really like to estimate the auto-correlation function, and I've coded up a solution in R (shown at the bottom of this question). Basically, I'm looping through all the observations in the dataset and for each observation I'm looking for other observations 1,2,3,... days in the past. If that observation exists, then I use this pair in my computation of the acf, otherwise I don't.

To verify my function works, I've been comparing it with R's base acf() on a simulated dataset where the spacing is equal. The agreement is close but not perfect (for n=10000 there's hardly any difference in the estimated acf's, but at n=100 I'm seeing differences as large as 0.05).

My question is: how is the acf computed, exactly? I realize it's a correlation between lagged observations, but what estimates do we use? For example, I'm using:

$$S_{X,Y} = \frac{\sum X_{i+t} X_i - (\sum X_{i+t}) (\sum X_i)/n}{n-1}$$ $$S_{X,X} = \frac{\sum X_i^2 - (\sum X_i)^2/n}{n-1}$$ $$\widehat{AR(1)} = S_{X,Y}/S_{X,X}$$

where $X_i$ is the current observation, $X_{i+t}$ is the observation that is exactly 1 day (or 2,3,...) in the past, and all summations are over pairs that have exactly that difference. So, I'm guessing that one of my $n$'s should be $n-1$ or vice versa, as asymptotically my estimator agrees with R's, but I can't figure it out. Any suggestions?

acfUnequal = function( data, data.col=5, maxLags=10 )
  lags = 1:maxLags
  data = data[!is.na(data[,data.col]),] #Can't use for acf anyways, and causes problems in calcs
  reqCols = c("Date", "Hour")
  test = reqCols %in% colnames(data)
  if( any(!test) )
    stop(paste0("Missing the following columns: ", paste(reqCols[!test],collapse=", ")))
    stop("data.col is larger than ncol(data)")

  diffT = as.numeric(diff(data$Date)) + diff(data$Hour)/2400
  Exy = rep(0,length(lags)) #E(XY), computed by adding up all X*Y then dividing by count
  Exx = rep(0,length(lags)) #E(X^2), computed by adding up all X^2 then dividing by count
  Eyy = rep(0,length(lags)) #E(Y^2), computed by adding up all Y^2 then dividing by count
  Ex = rep(0,length(lags)) #E(X), computed by adding up all X then dividing by count
  cnt = rep(0,length(lags)) #Count, used to compute E(XY), E(X), E(Y)
  for(i in 2:nrow(data))
    diffCurr = 0
    j = i
    for(lag in lags)
      while(diffCurr<lag-0.05 & j>=2)
        j = j-1
        diffCurr = diffCurr + diffT[j]
      if(!diffCurr>lag+.05) #time lag within 0.05 days detected, use data to compute acf
        Exy[lag] = Exy[lag] + data[i,data.col]*data[j,data.col]
        Exx[lag] = Exx[lag] + data[i,data.col]*data[i,data.col]
        Eyy[lag] = Eyy[lag] + data[j,data.col]*data[j,data.col]
        Ex[lag] = Ex[lag] + data[i,data.col]
        Ey[lag] = Ey[lag] + data[j,data.col]
        cnt[lag] = cnt[lag] + 1
  sxy = (Exy - Ex*Ey/cnt)/(cnt-1)
  sxx = (Exx - Ex^2/cnt)/(cnt-1)
  acf = sxy/sxx

Ok, I figured out the problem. In my formula, I'm essentially estimating $E(X)$ twice: once by computing $\sum X_i$ and once by computing $\sum X_{i+t}$. In R, $E(X)$ is estimated by using all the observations instead of these two subsets, and the divisor is the number of observations used to estimate $E(X)$. Also, $S_{XX}$ is estimated using all the observations as well.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.