We have seen examples of regression where the basic assumption of uncorrelated residuals is compromised. Finding the autocorrelation of the residuals using Durbin–Watson is one way to diagnose the correlation. Here, we perform a step-by-step estimation.

Step 1: Plot the data

plot(Nif_data$Year, Nif_data$Index, xlab = "Year", ylab = "Index")

Step 2: Develop a regression model

fit <- lm(Index ~ Year, data=Nif_data)
summary(fit)

Call:
lm(formula = Index ~ Year, data = Nif_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3410.3  -544.5   -96.5   507.6  5603.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -438.128     27.801  -15.76   <2e-16 ***
Year         566.726      2.243  252.65   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1020 on 5352 degrees of freedom
Multiple R-squared:  0.9226,	Adjusted R-squared:  0.9226 
F-statistic: 6.383e+04 on 1 and 5352 DF,  p-value: < 2.2e-16

Step 3: Estimate Residuals

Nif_data$resid <- resid(fit)

Step 4: Durbin–Watson (D-W) Statistics

D-W statistic is the sum of differences between successive residuals squared divided by the sum of residuals squared.

D-W Statistics = sum (ei - ei-1)2 / sum(ei2)

sum(diff(Nif_data$resid)^2) / sum(Nif_data$resid^2)

0.006301032

R can do better – using the ‘durbinWatsonTest’ function from the library ‘car’.

library(car)
fit <- lm(Index ~ Year, data=Nif_data)
durbinWatsonTest(fit)

lag Autocorrelation D-W Statistic p-value
   1       0.9936623   0.006301032       0
 Alternative hypothesis: rho != 0

We continue from the Aus beer production data set. In the last post, we manually built scatter plots at a few lags (lag = 0, 1, 2, and 4) and introduced the ‘acf’ function.

There is another R function, ‘gg_lag’, which can be an even better illustration.

new_production %>% gg_lag(Beer)

For example, the lag = 1 plot (top left). The yellow cluster represents the plot of Q4 on the Y-axis against Q3 of the same year on the X-axis. On the other hand, the purple cluster represents the plot of Q1 on the Y-axis against Q4 of the previous year on the X-axis.

Similarly, for the lag = 2 plot (top middle), The yellow represents the plot of Q4 against Q2 (two quarters earlier) of the same year. The purple represents the plot of Q1 against Q3 of the previous year.

The lag = 4 plot shows a strong positive correlation (everyone on the diagonal line). This is not surprising; it plots Q4 of one year against Q4 of the previous year, Q1 of one year with Q1 of the previous year, etc.

As we have seen, autocorrelation is the correlation of a variable in a time series with itself but at a time lag. For example, how are the variables at time ts correlated to the t-1s? We will use the aus_production dataset from the R library ‘fpp3’ to illustrate the concept of autocorrelation.

library(fpp3)
str(aus_production)

 $ Quarter    : qtr [1:74] 1992 Q1, 1992 Q2, 1992 Q3, 1992 Q4, 1993 Q1, 1993 Q2, 1993 Q3, 1993 Q4, 1994 Q1, 1994 Q2, 1994 Q3, 1994 Q4, 1995 Q1, 1995 Q2, 1995 Q3,...
 $ Beer       : num [1:74] 443 410 420 532 433 421 410 512 449 381 ...
 $ Tobacco    : num [1:74] 5777 5853 6416 5825 5724 ...
 $ Bricks     : num [1:74] 383 404 446 420 394 462 475 443 421 475 ...
 $ Cement     : num [1:74] 1289 1501 1539 1568 1450 ...
 $ Electricity: num [1:74] 38332 39774 42246 38498 39460 ...
 $ Gas        : num [1:74] 117 151 175 129 116 149 163 138 127 159 ...

We will use the beer production data.

Let’s plot the production data with itself without any lag in time.

plot(B_data, B_data, xlab = "Beer Production at i", ylab = "Beer Production at i")

There is no surprise here; the data is in perfect correlation. In the next step, we will give a lag of one time interval, i.e., a plot of (2,1), (3,2), (4,3), etc. The easiest way to achieve this is to remove the first element of the vector and plot against the vector with the last element removed.

plot(B_data[-1], B_data[-n], xlab = "Beer Production at i - 1", ylab = "Beer Production at i")

You clearly see a lack of correlation. What about a plot with a lag of 2 and 4?

There is a negative correlation compared with the time series 2 quarters ago.

There is a good correlation compared with the time series 4 quarters ago.

The whole process is established using the autocorrelation function (ACF).

autocorrelation <- acf(B_data, lag.max=10, plot=FALSE)
plot(autocorrelation,
     main="Autocorrelation",
     xlab="Lag Parameter",
     ylab="ACF")

ACF at large parameter 1 indicates how successive values of beer production relate to each other, 2 indicates how production two periods apart relate to each other, etc.