Calculating Autocorrelation Function from scratch in R

Question

I want to understand how to calculate Autocorrelation Function from scratch in R. How can I use cor(x=y, y=lag(x=y, k=2)) to get ACF when y is a ts object?

I've already tried all choices for use argument [use = "complete.obs" # c("everything", "all.obs", "complete.obs", "na.or.complete", "pairwise.complete.obs")].

set.seed(1)
y <-ts(data = rnorm(20), start = c(2010, 1), frequency = 4)
y

# Qtr1        Qtr2        Qtr3        Qtr4
# 2010  0.91897737  0.78213630  0.07456498 -1.98935170
# 2011  0.61982575 -0.05612874 -0.15579551 -1.47075238
# 2012 -0.47815006  0.41794156  1.35867955 -0.10278773
# 2013  0.38767161 -0.05380504 -1.37705956 -0.41499456
# 2014 -0.39428995 -0.05931340  1.10002537  0.76317575

acf(x=y, plot=FALSE)
# Autocorrelations of series ‘y’, by lag
# 
# 0.00   0.25   0.50   0.75   1.00   1.25   1.50   1.75   2.00   2.25   2.50   2.75   3.00   3.25 
# 1.000 -0.122 -0.185 -0.049  0.147 -0.283 -0.255  0.212  0.097 -0.120 -0.181  0.286 -0.063  0.094 

cor(
      x   = y
    , y   = lag(x=y, k=2)
    , use = "complete.obs" # c("everything", "all.obs", "complete.obs", "na.or.complete", "pairwise.complete.obs")
    )

# [1] 1

Answer 1

It looks like you have two problems.

First problem: cor(x, lag(x, k = k)) will always be 1 because the two series are not temporally aligned. You need to use ts.union or cbind.ts before using cor

X <- ts.union(yt = y, yt2 = lag(x = y, k =  2))
head(X)
##            yt      yt2
## [1,]       NA -0.62645
## [2,]       NA  0.18364
## [3,] -0.62645 -0.83563
## [4,]  0.18364  1.59528
## [5,] -0.83563  0.32951
## [6,]  1.59528 -0.82047
tail(X)
##              yt      yt2
## [17,]  1.124931 -0.01619
## [18,] -0.044934  0.94384
## [19,] -0.016190  0.82122
## [20,]  0.943836  0.59390
## [21,]  0.821221       NA
## [22,]  0.593901       NA

The problem is that without time index, the raw data are the same (with a temporal shift). You can test yourself

x1 <- as.vector(y)
x2 <- as.vector(lag(y, k = 2))
all.equal(x1, x2)
## [1] TRUE

Therefore, if you to compute the correlation coefficient between a a time series and its lag, you can use X (created using ts.union)

cor(X[, 1], X[, 2], use = "complete.obs")
## [1] -0.19018

Still the results is different from acf(y, plot = FALSE)$acf[3]

acf(y, plot = FALSE)$acf[3]
## [1] -0.18521

Which bring us to the second reason why you can't use cor to compute acf:

The mathematical definition of the Acf suppose a least a second-order stationarity (mean and variance equal for each lag):

But if you use the standard implementation of cor, for each series, a different value for mean and variance will be computed (the value for the lag will be different from the original series).

c0 <- var(y)
m <- mean(y)
n <- length(y)
ct <- sum((X[, 1] - m) * (X[, 2] - m), na.rm = TRUE) / (n - 1)
(rt <- ct / c0)
## [1] -0.18521

all.equal(rt, acf(y, plot = FALSE)$acf[3])
## [1] TRUE

Answer 2

Got the correct answer following the link provided by @khashaa. Thanks

acf(x=y, plot=FALSE)
# Autocorrelations of series ‘y’, by lag
# 
# 0.00   0.25   0.50   0.75   1.00   1.25   1.50   1.75   2.00   2.25   2.50   2.75   3.00   3.25 
# 1.000 -0.122 -0.185 -0.049  0.147 -0.283 -0.255  0.212  0.097 -0.120 -0.181  0.286 -0.063  0.094 

sum((y-mean(y))*(lag(x=y, k=0)-mean(y)))/sum((y-mean(y))^2)
[1] 1
sum((y-mean(y))*(lag(x=y, k=1)-mean(y)))/sum((y-mean(y))^2)
[1] -0.1222859
sum((y-mean(y))*(lag(x=y, k=2)-mean(y)))/sum((y-mean(y))^2)
[1] -0.1852114
sum((y-mean(y))*(lag(x=y, k=3)-mean(y)))/sum((y-mean(y))^2)
[1] -0.04940401