Lecture 34: Diagnosing Autocorrelated Errors

161.251 Regression Modelling

Diagnosing Autocorrelated Errors

In Lectures 32/33 we used trend, seasonal terms and weather predictors to explain the mean response through time. A fitted curve can follow the data quite closely and still leave several consecutive observations above it, followed by several below it.

If today’s response is higher than the model predicts, does that tell us anything about tomorrow’s deviation? If it does, treating those observations as independent gives us the wrong calculation of precision.

By the end of this lecture, students should be able to:

1. Define lagged autocorrelation and distinguish errors from residuals.
2. Use a variance calculation to explain why dependence changes precision.
3. Explain and use residual plots, runs tests, Durbin-Watson, ACF plots, and Ljung-Box tests.
4. Recognise when dependence remains after a reasonable mean model has been fitted.

What Is Autocorrelation?

Write our regression model as

\[Y_t = \mu_t + \varepsilon_t, \qquad \mu_t = \mathbf{x}_t^\top\boldsymbol{\beta}.\]

The predictors in \(\mathbf{x}_t\) describe the mean. Autocorrelation describes a relationship between errors at different times, after accounting for that mean.

We cannot observe \(\varepsilon_t\), because we do not know the true mean. We work with the fitted residuals

\[e_t = y_t - \widehat{\mu}_t.\]

A high \(R^2\) describes how closely the fitted values track the responses. It does not tell us whether the remaining deviations are related through time.

What Does a Lag Mean?

Lag \(k\) means a separation of \(k\) observation intervals. For quarterly data, lag 1 is one quarter and lag 4 is one year. The observations must be equally spaced for this interpretation.

Suppose the errors have constant variance \(\sigma^2\) and their covariance depends on the time gap, rather than the calendar date. Then

\[\rho_k = \operatorname{Corr}(\varepsilon_t,\varepsilon_{t-k}) = \frac{\operatorname{Cov}(\varepsilon_t,\varepsilon_{t-k})}{\sigma^2}.\]

• \(\rho_1 > 0\): nearby errors tend to have the same sign, producing persistent runs.
• \(\rho_1 < 0\): nearby errors tend to alternate in sign.
• Independent errors have \(\rho_k = 0\) for every \(k > 0\). Zero correlation alone does not establish independence.

A Classical Example: Consumer Expenditure and Money Stock

This example studies consumer expenditure (ConsExp) and the stock of money (MoneyStk) in billions of US dollars.

Download MoneyStock.csv

MoneyStock <- read.csv("../data/MoneyStock.csv", header = TRUE)
MoneyStock |>
    head() |>
    kable()

Year	Quarter	ConsExp	MoneyStk	Year.Q
1952	1	214.6	159.3	1952.00
1952	2	217.7	161.2	1952.25
1952	3	219.6	162.8	1952.50
1952	4	227.2	164.6	1952.75
1953	1	230.9	165.9	1953.00
1953	2	233.3	167.9	1953.25

pairs(MoneyStock[3:5])

Fit the Regression Model

For quarter \(t\), we fit

\[\mathrm{ConsExp}_t = \beta_0 + \beta_1\mathrm{MoneyStk}_t + \varepsilon_t.\]

MoneyStock_lm1 <- lm(ConsExp ~ MoneyStk, data = MoneyStock)
summary(MoneyStock_lm1)

Call:
lm(formula = ConsExp ~ MoneyStk, data = MoneyStock)

Residuals:
   Min     1Q Median     3Q    Max 
-7.176 -3.396  1.396  2.928  6.361 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -154.7192    19.8500  -7.794 3.54e-07 ***
MoneyStk       2.3004     0.1146  20.080 8.99e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.983 on 18 degrees of freedom
Multiple R-squared:  0.9573,    Adjusted R-squared:  0.9549 
F-statistic: 403.2 on 1 and 18 DF,  p-value: 8.988e-14

MoneyStk and ConsExp have a strong linear association: \(R^2\) is about 0.96. The reported slope also looks precise. That precision calculation assumes that knowing the error in one quarter tells us nothing about another quarter’s error.

Residual Plot Over Time

Before running a test, follow the residuals from left to right. Look for stretches above or below zero, alternation, and any remaining trend or seasonal pattern.

plot(residuals(MoneyStock_lm1), type = "b", ylab = "Residuals", xlab = "Quarter")
abline(h = 0, lty = 2)

The residuals wander slowly above and below zero. Adjacent quarters often have deviations in the same direction. Some runs occur by chance even with independent errors; we need to judge whether this much persistence is plausible.

Why Dependence Changes Precision

Start with a simpler calculation. Suppose two observations measure the same mean \(\mu\), each with error variance \(\sigma^2\) and correlation \(\rho\). Their average has variance

\[\begin{aligned} \operatorname{Var}\left(\frac{Y_1+Y_2}{2}\right) &= \frac{1}{4}\left[\operatorname{Var}(Y_1)+\operatorname{Var}(Y_2) +2\operatorname{Cov}(Y_1,Y_2)\right] \\ &= \frac{\sigma^2}{2}(1+\rho). \end{aligned}\]

With independent errors, \(\rho=0\) and averaging halves the variance. With \(\rho=0.8\), the variance is \(0.9\sigma^2\): the second observation adds little precision because the two errors tend to move together.

The standard error is then \(\sqrt{1.8}\approx1.34\) times the value calculated under independence. Counting two observations does not give us two independent pieces of information.

What Changes in Regression?

Remember the matrix form \(\mathbf{Y}=X\boldsymbol{\beta}+\boldsymbol{\varepsilon}\). Treating \(X\) as fixed, let \(\Sigma=\operatorname{Var}(\boldsymbol{\varepsilon}\mid X)\). The OLS estimate still satisfies

\[\widehat{\boldsymbol{\beta}}=\boldsymbol{\beta}+(X^\top X)^{-1}X^\top\boldsymbol{\varepsilon}.\]

Its variance is

\[\operatorname{Var}(\widehat{\boldsymbol{\beta}}\mid X) =(X^\top X)^{-1}X^\top\Sigma X(X^\top X)^{-1}.\]

Independent errors with equal variance give \(\Sigma=\sigma^2 I\), reducing this to \(\sigma^2(X^\top X)^{-1}\). Autocorrelation introduces non-zero off-diagonal covariances into \(\Sigma\).

If \(E(\boldsymbol{\varepsilon}\mid X)=\mathbf{0}\), OLS remains unbiased. The usual standard errors can still be wrong, so the reported confidence intervals and tests become unreliable. For regression coefficients, the direction of that error depends on both the predictors and the dependence pattern.

Diagnostic Tools

We have seen a pattern that could make our precision calculation unreliable. Each diagnostic summarises a different feature of that pattern.

Diagnostic	Question about the ordered residuals
Runs test	Do residuals stay on the same side of their median unusually long, or switch sides unusually often?
Durbin-Watson	How different is each residual from the one immediately before it?
ACF	How strong is the correlation at each time gap?
Ljung-Box	Is there evidence of correlation across a chosen group of lags?

Keep the time order. Rearranging observations leaves the fitted regression unchanged but changes these diagnostics. Also check the mean model: a missing trend or seasonal term can produce a residual pattern of its own.

Runs Test

A run is an uninterrupted sequence on the same side of a threshold. For example, \(+++\;--\;++\;-\) contains four runs, although it contains eight observations.

The default RunsTest() divides numerical residuals at their median. Write \(n_+\) for the number above the median and \(n_-\) for the number at or below it. Under random ordering, the expected number of runs is

\[E(R)=1+\frac{2n_+n_-}{n_++n_-}.\]

Too few runs indicate clustering on the same side; too many indicate excessive alternation. The test compares the observed count with its distribution under random ordering.

Runs in the MoneyStock Residuals

There are ten residuals on each side of the median, so under random ordering we expect

\[E(R)=1+\frac{2(10)(10)}{20}=11\text{ runs}.\]

RunsTest(residuals(MoneyStock_lm1))

    Runs Test for Randomness

data:  residuals(MoneyStock_lm1)
runs = 5, m = 10, n = 10, p-value = 0.008985
alternative hypothesis: true number of runs is not equal the expected number
sample estimates:
median(x) 
 1.395704 

We observe only five runs, with \(P\approx0.009\). This supports the impression of long stretches on the same side. The default test is two-sided, so it can detect either unusually few or unusually many runs. It does not measure the size of the residuals or identify the process producing them.

Durbin-Watson Test

Runs use only which side each residual falls on. Durbin-Watson uses the sizes of the changes between neighbouring residuals:

\[d=\frac{\sum_{t=2}^{n}(e_t-e_{t-1})^2}{\sum_{t=1}^{n}e_t^2}.\]

The denominator measures the overall residual variation. The numerator measures how much residuals change from one observation to the next.

• Positive autocorrelation gives small neighbouring differences, so \(d\) tends towards 0.
• With no lag-1 autocorrelation, \(d\) is typically near 2.
• Negative autocorrelation gives larger neighbouring differences, so \(d\) tends towards 4.

These values guide interpretation; they are not significance cut-offs.

Durbin-Watson for MoneyStock

We test \(H_0:\rho_1=0\) against \(H_A:\rho_1>0\). The one-sided alternative targets the positive autocorrelation suggested by the time plot.

dwtest(MoneyStock_lm1, alternative = "greater")

    Durbin-Watson test

data:  MoneyStock_lm1
DW = 0.32821, p-value = 2.303e-08
alternative hypothesis: true autocorrelation is greater than 0

The statistic is \(d\approx0.33\), with a very small P-value. The residuals change much less between adjacent quarters than we would expect from their overall variation.

In dwtest(), alternative = "greater" refers to correlation greater than zero. Evidence for that alternative comes from a small value of \(d\). The test’s reference distribution accounts for the regression predictors; we do not simply declare every value below 2 significant.

Measuring Correlation at Each Lag

Durbin-Watson focuses on neighbouring observations. The autocorrelation function (ACF) lets us examine longer gaps as well. The sample autocorrelation at lag \(k\) is

\[r_k=\frac{\sum_{t=k+1}^{n}(e_t-\bar e)(e_{t-k}-\bar e)} {\sum_{t=1}^{n}(e_t-\bar e)^2}.\]

The numerator pairs residuals \(k\) observations apart. Products are positive when both deviations from \(\bar e\) have the same sign, and negative when they have opposite signs. The denominator puts the result on a correlation scale. Lag 0 is always 1.

For OLS with an intercept, \(\bar e=0\) apart from numerical rounding. Expanding the squared differences in Durbin-Watson then gives

\[d=2(1-r_1)-\frac{e_1^2+e_n^2}{\sum_{t=1}^{n}e_t^2} \approx2(1-r_1).\]

The approximation drops the two endpoint terms. It can be noticeably different in a short series, so use it to understand the direction of \(d\), rather than to replace the calculation.

ACF Plot

acf(residuals(MoneyStock_lm1), main = "ACF of MoneyStock residuals")

The lag-1 value is about 0.75, and positive correlation persists across several quarters. This agrees with the slow movement in the time plot. We should also check whether a missing time trend could explain part of it.

Reading the ACF Reference Lines

Under an uncorrelated white-noise reference model, the default approximate 95% reference limits at an individual non-zero lag are

\[\pm\frac{1.96}{\sqrt n}.\]

Here \(n=20\), giving limits of about \(\pm0.44\). The lag-1 correlation is well outside them. A sequence of large correlations deserves more attention than one isolated bar just beyond a line.

These are approximate, pointwise limits, not a simultaneous 95% band for the whole plot. With many lags, some bars can cross by chance. Fitting the regression also affects the behaviour of its residuals, so treat the lines as a visual reference. The acf() plotting documentation describes the white-noise assumption behind them.

Ljung-Box Test

If several ACF bars look large, we can combine them. For a chosen maximum lag \(m\), the hypotheses are

\[H_0:\rho_1=\cdots=\rho_m=0, \qquad H_A:\text{at least one of these correlations is non-zero}.\]

The Ljung-Box statistic is

\[Q=n(n+2)\sum_{k=1}^{m}\frac{r_k^2}{n-k}.\]

Squaring prevents positive and negative correlations from cancelling. Larger correlations give a larger \(Q\) and more evidence against the null hypothesis. Choose \(m\) using the time scale and the available sample size, before inspecting a collection of P-values.

Ljung-Box for MoneyStock

Here we use \(m=10\), combining correlations over gaps from one to ten quarters.

Box.test(residuals(MoneyStock_lm1), lag = 10, type = "Ljung-Box")

    Box-Ljung test

data:  residuals(MoneyStock_lm1)
X-squared = 64.397, df = 10, p-value = 5.286e-10

We obtain \(Q\approx64.4\) and a very small P-value. With the default fitdf = 0, Box.test() compares this with an approximate \(\chi^2_{10}\) distribution.

There are only 20 observations here, and the residuals come from an estimated regression, so treat that P-value as an approximate diagnostic summary. The time plot, runs test and Durbin-Watson already give substantial evidence of persistence. This result does not establish that every lag is correlated or tell us which change to the model will address it.

Work Through the Diagnostics

1. If \(d=0.6\), what lag-1 correlation does \(d\approx2(1-r_1)\) suggest?

. . .

Rearranging gives \(r_1\approx1-d/2=0.7\), suggesting strong positive correlation between neighbours. This is an approximation, not a P-value.

. . .

2. If \(r_1=0\) but \(r_4\) is large, which diagnostic displays that pattern, and which can test a group of lags containing it?

. . .

The ACF displays the lag-4 correlation. Ljung-Box with \(m\geq4\) includes it in a joint test; a lag-1 diagnostic alone can miss it.

. . .

3. Why would shuffling complete rows of MoneyStock leave the coefficients and \(R^2\) unchanged but alter the residual diagnostics?

. . .

OLS minimises the same sum of squares after reordering. Runs, neighbouring differences and lagged products depend on which observations follow one another, so the diagnostics need the actual time order.

A Modern Comparison: Seoul Bike Residuals

The same diagnostic logic applies to the richer daily bike model from Lectures 32/33, which combined two Fourier pairs with the time trend and weather predictors.

Download seoul_bike_daily.csv

seoul_daily <- read_csv("../data/seoul_bike_daily.csv", show_col_types = FALSE) |>
    mutate(date = as.Date(date)) |>
    arrange(date) |>
    filter(functioning_day_all == "yes") |>
    mutate(S1 = sin(2 * pi * day_index/365), C1 = cos(2 * pi * day_index/365), S2 = sin(4 *
        pi * day_index/365), C2 = cos(4 * pi * day_index/365))

seoul_mean <- lm(rented_bike_count ~ day_index + S1 + C1 + S2 + C2 + mean_temperature +
    total_rainfall + total_solar_radiation, data = seoul_daily)

seoul_diagnostics <- seoul_daily |>
    mutate(residual = residuals(seoul_mean))

We fitted the mean model to all 352 functioning days. The dates remain visible here, including the gaps left by 13 outage days.

Seoul: One Uninterrupted Period

For the daily diagnostics, select the longest uninterrupted period: 1 December 2017 through 10 April 2018. Otherwise, an ordinary ACF would treat consecutive retained rows as though they were always one day apart.

Refit the same mean formula on these 131 consecutive days. All four diagnostics below use residuals from this refitted model, so the Durbin-Watson test and the other summaries refer to the same observations and fit.

seoul_run_data <- seoul_daily |>
    filter(date >= as.Date("2017-12-01"), date <= as.Date("2018-04-10"))

seoul_run_model <- update(seoul_mean, data = seoul_run_data)
seoul_run_residuals <- residuals(seoul_run_model)

Seoul: Runs Above and Below the Median

Divide the residuals at their median, 126.6 bikes. Each uninterrupted colour block is one run.

There are 65 residuals above the median and 66 at or below it. Compare the actual run count with \(E(R)=1+2n_+n_-/n\).

seoul_runs <- RunsTest(seoul_run_residuals)

Observed runs	Expected runs	Test statistic z	P-value
55	66.5	-1.93	0.0537

The count is 55 runs; \(z\) standardises its departure from the expected count, with a continuity correction. The two-sided P-value is just above 0.05. This test uses only the median-side sequence, so it can miss dependence visible in the residual magnitudes.

Seoul: Durbin-Watson

Use the sizes of the changes between consecutive daily residuals:

\[d=\frac{\sum_{t=2}^{n}(e_t-e_{t-1})^2}{\sum_{t=1}^{n}e_t^2} =\frac{609135793}{453162261} \approx 1.34.\]

seoul_dw <- dwtest(seoul_run_model, alternative = "greater")
seoul_dw

    Durbin-Watson test

data:  seoul_run_model
DW = 1.3442, p-value = 7.693e-06
alternative hypothesis: true autocorrelation is greater than 0

The statistic is below 2, with strong evidence against \(H_0:\rho_1=0\) in favour of positive autocorrelation. Unlike the runs test, Durbin-Watson uses how far neighbouring residuals move, as well as their direction.

Seoul: ACF of the Daily Residuals

acf(seoul_run_residuals, main = "ACF of Seoul residuals: refitted on consecutive days")

The lag-1 correlation is \(r_1\approx 0.327\). The approximate pointwise reference limits are

\[\pm\frac{1.96}{\sqrt{131}}\approx\pm 0.171.\]

Neighbouring daily residuals remain positively correlated after refitting the mean. Read the other lags from the plot as well; the reference lines are approximate and do not form a simultaneous band.

Seoul: Ljung-Box Across Two Weeks

Combine the correlations at lags 1 through 14, rather than testing only neighbouring days:

\[Q=n(n+2)\sum_{k=1}^{14}\frac{r_k^2}{n-k}.\]

seoul_ljung_box <- Box.test(seoul_run_residuals, lag = 14, type = "Ljung-Box")

Statistic Q	Reference df	P-value
45.19	14	3.79e-05

The small P-value gives evidence of remaining autocorrelation across these lags. It does not imply that every lag is correlated. This is an approximate diagnostic on fitted regression residuals, using the default fitdf = 0.

The runs test gave weaker evidence because it discarded residual magnitudes. Taken together, the time plot, Durbin-Watson, ACF and Ljung-Box suggest that refitting the mean has not removed the dependence. These results concern the 131-day fit, rather than establishing adequacy over the entire gappy year.

Diagnostic Summary

The time plot comes first, then the ACF. The runs test, Durbin-Watson and Ljung-Box each summarise a different feature of the ordered residuals, so a collection of small P-values supports what we can already see in the plots. None of them identifies the error process that produced the pattern.

Before changing the error model, check the mean: a missing trend, seasonal term or outage can leave a residual pattern of its own. For both MoneyStock and Seoul we fitted reasonable mean models and the dependence remained, so we treat it as part of the errors and model it directly.

Looking Ahead

Every model we have fitted this semester assumed independent errors. Flexible curves, factors and weights changed the mean or the variance; the correlation between errors never entered the model. After this lecture we know how to detect that correlation, and in MoneyStock and Seoul it is clearly there.

Next lecture we model it. An AR(1) error model, fitted by generalised least squares (GLS), keeps the mean formula and adds a correlation structure to the errors instead. We will fit the tourism model with gls() and corAR1(), compare the standard errors with the ordinary fit, and check whether the normalised residuals look closer to white noise.