# Estimating auto-correlation with unequally spaced data

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=", ")))
if(ncol(data)<data.col)
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
return(acf)
}


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.