R Time Series Issue: Autocorrelation Function

In two earlier blog posts, I talked about using R for time series analysis and in particular I pointed out that the arima() function does some slightly strange things when calculating residuals.

In this post, I'd like to draw attention to the acf() function which plots the autocorrelation function for a particular time series. In actual fact, the default setting is that acf() will plot the autocorrelation function, but the plotting can be suppressed by specifying the argument plot=FALSE. If the plotting is suppressed then R will instead report the values of autocorrelation coefficients for up to 30 lags. If you want to analyse more than 30 lags, say 50, then you need to specify max.lags=50.

The standard formula for calculating autocorrelation coefficients is \[r_{k} = \frac{\sum_{t=1}^{n-k}(z_{t} - \bar{z})(z_{t+k} - \bar{z})}{\sum_{t=1}^{n} (z_{t} - \bar{z})^{2}}\] where $r_{k}$ is the autocorrelation coefficient at lag $k$, $z_{t}$ is a time-series, and $\bar{z}$ is the mean of the time-series.

The formula above is not the forumla for calculating autocorrelation coefficients because there are other ways of calculating $r_{k}$. Pankratz (1983) (p.35) refers me to Jenkins and Watts (1968) for a discussion and evaluation of the alternatives.

A problem is that the exact equation used by acf() itself is given in Venables and Ripley's Modern Applied Statistics with S and since I do not have access to this book, how do I know which formula is being computed when I call acf() in R? Reviewing the source code in RStudio by pressing F2 does not readily reveal the formula used by the function. Not knowing the exact details behind the acf() function is kind of like working with a black-box!

To bring the point further, if I were to present a piece of research that involved estimating autocorrelation coefficients using the acf() function and then someone came along and quizzed me about the formulas used to carry out the research, how could I respond? Er... I just used a function and I have no clue how it was programmed, sorry about that!? That wouldn't be a nice situation to find oneself in.

The basic point I'm trying to make is that it's good practice to be aware of the formulas being used in the analysis that you're carrying out. Having at least a fair idea about how the formulas could be used "by hand" is a welcomed bonus.

I'm not alone in trying to decipher what acf() actually does. Professor of Statistics Galit Shmuéli asked a question on Cross Validated also with the intention of finding out exactly what it does.

To get around the problem I decided to write my own function in R that will calculate autocorrelation coefficients for a series of data using the aforementioned standard formula. I've tested the function in order to confirm that it works correctly by replicating various results for estimated acf's reported in Pankratz (1983) - as well as some other books and papers.

Here is an R script which contains the function that I wrote to calculate autocorrelation coefficients. I've written the code in such a way that the code and the actual formula look similar to each other. To test the function, I've used two datasets from Box and Jenkins (1976). After you run the script, you can compare the results to those found in the book if you have a copy of it.

I've also devised a basic plotting scheme that overcomes a weird plotting convention associated with acf(). Professor Stoffer at the Department of Statistics at the University of Pittsburgh has already shone light on this issue on his website (see Issue #6).