Confidence bands in case of fitting ARIMA in R?

Question

I want to look at the acf and pacf of my data, to identify the model for my mean equation, so I want to fit an ARMA for my mean equation and later on model the conditional variance by a ARCH/GARCH (I know I have to do jointly model estimation). In the first step I want to look at the ACF and PACF for identifying, I used the standard Acf of the forecast library, but I noticed, that the confidence bands in these plots are given for testing randomness and not for fitting a ARMA.

As wikipedia says:

Correlograms are also used in the model identification stage for fitting ARIMA models. In this case, a moving average model is assumed for the data and the following confidence bands should be generated:

$ \pm z_{1-\alpha/2}\sqrt{\frac{1}{N}\left(1+2\sum_{i=1}^{k} y_i^2\right)} $ How can I get these confidence bands, which increase as the lag increases?

Answer 1

Use acf from the stats package (or Acf from the forecast package) with ci.type="ma". Note that some people use the simpler approximation all the time - it's just to give an idea what models might be worth considering so accuracy isn't so important.

Bartlett's approximation (the one you quote from Wikipedia) is only relevant to examining the autocorrelation function: the confidence interval for a lag $q$ is given assuming, as a null hypothesis, a moving average process of order $q-1$; it's conditional upon the estimated autocorrelations of all previous lags. (So note that it's not especially relevant to deciding between, say, an ARMA(1,1) & an ARMA (1,2).)

You might suppose a similar formula for confidence intervals on the partial autocorrelations, mutatis mutandis; but you'd be wrong: if you assume an autoregressive process of order $p-1$, the standard errors on the partial autocorrelations are asymptotically $\frac{1}{\sqrt{n}}$, where $n$ is the number of observations. Quenouille (1949), "Approximate tests of correlation in time-series", JRSS B, 11, 1.

Answer 2

I've included an example below to show one method of how to calculate Bartlett's approximations and add them to a graph of the autocorrelation function.

In the example, I've done things by hand rather than rely on a particular package hence the code is longer than it perhaps need be. I don't claim the code to be efficient, but it does get the correct results.

The results in my example can be compared for accuracy to Figure C1.2. in Case 1 of Pankratz (1983) which uses the same dataset that I've used. Don't worry if you don't have the book because the content of Case 1 is available free to download.

Note that with slight modifications, the code below can be adapted to plot the standard errors on the pacf, which Scortchi has correctly pointed out should equal $n^{-1/2}$.

# Import data from the web
inventories <- scan("http://robjhyndman.com/tsdldata/books/pankratz.dat", skip=5, nlines=5, sep="")
# Calculate sample size and mean
n <- length(inventories)
mean.inventories <- sum(inventories)/n
# Express the data in deviations from the mean
z.bar <- rep(mean.inventories,n)
deviations <- inventories - z.bar
# Calculate the sum of squared deviations from the mean
squaredDeviations <- deviations^2
sumOfSquaredDeviations <-sum(squaredDeviations)
# Create empty vector to store autocorrelation coefficients
r <- c()
# Use a for loop to fill the vector with the coefficients
for (k in 1:n) {
  ends <- n - k
  starts <- 1 + k
  r[k] <- sum(deviations[1:(ends)]*deviations[(starts):(n)])/sumOfSquaredDeviations
}
# Create empty vector to store Bartlett's standard errors
bart.error <- c()
# Use a for loop to fill the vector with the standard errors
for (k in 1:n) {
  ends <- k-1
  bart.error[k] <- ((1 + sum((2*r[0:(ends)]^2)))^0.5)*(n^-0.5)
}
# Plot the autocorrelation function
plot(r[1:(n/4)], 
     type="h", 
     main="Autocorrelation Function", 
     xlab="Lag", 
     ylab="ACF",
     ylim=c(-1,1),
     las=1)
abline(h=0)
# Add Bartlett's standard errors to the plot
lines(2*bart.error[1:(n/4)])
lines(2*-bart.error[1:(n/4)])

After running the code you should see the following plot: