Linear regression with autocorrelated noise

Consider two time series and such that:

where is noise:

By iteration of Equation 2, we see that has gaussian unconditional distribution:

so that individual observations of are distributed according to a perfectly specified linear model.

This does not mean that, given observational data , we are allowed to make standard linear model assumptions to perform valid inference on the parameters and of Equation 1 and Equation 3. Since the noise terms are not independent draws from a single distribution, but are rather autocorrelated, the usual OLS variance estimate under linear model assumptions will be biased, as we show below ¹.

It is fairly easy to work out the consequences of autocorrelation. Suppose, more generally, that the error term is a stationary time series with unconditional mean and unconditional variance . The OLS estimate of is²:

which is unbiased since . The estimate of the noise variance , on the other hand:

where as usual, and we have used the fact that (since each has the same unconditional variance ). Hence the OLS estimate is biased if .

Similarly, the variance-covariance matrix of the OLS estimator is:

whereas its OLS estimate is:

which is biased for .

Even though the variance estimators are themselves biased, the biases could still vanish in the asymptotic limit. This is the case for , as we can see by rewriting:

where we have used the projector properties of to recast the trace in terms of a symmetric operator. In principle, nothing prevents the operator above to have eigenvalues, which would make the estimator asymptotically biased³. In realistic cases, one expects the correlations to decay exponentially with ⁴ , in which case the trace is bounded to be of , and as .

For things are not so favorable. It is enough to consider a special case of a plain intercept term: . In this case, we find with some manipulations:

Since , we see that:

which amounts to say that is asymptotically biased⁵.

Illustration

The (foldable) block below defines helpers to simulate the results of linear regression on data generated according to . These are the same functions used in my previous post on misspecification and sandwich estimators – slightly adapted to the current case.

library(dplyr)
library(ggplot2)

rxy_fun <- function(rx, f, reps) {
    res <- function(n) {
        x <- rx(n)  # X has marginal distribution 'rx'
        y <- f(x) + reps(x)  # Y has conditional mean 'f(x)' and noise 'reps(x)'
        return(tibble(x = x, y = y))  
    }
    return(structure(res, class = "rxy"))
}

plot.rxy <- function(x, N = 1000, seed = 840) {
    set.seed(seed)
    
    ggplot(data = x(N), aes(x = x, y = y)) +
        geom_point(alpha = 0.3) + 
        geom_smooth(method = "lm", se = FALSE)
}

lmsim <- function(rxy, N = 100, vcov = stats::vcov, B = 1e3, seed = 840) 
{ 
    set.seed(seed)
    
    res <- list(
        coef = matrix(nrow = B, ncol = 2), 
        vcov = vector("list", B),
        sigma2 = numeric(B)
        )
    colnames(res$coef) <- c("(Intercept)", "x")
    class(res) <- "lmsim"
                                
    for (b in 1:B) {
        .fit <- lm(y ~ ., data = rxy(N))
        res$coef[b, ] <- coef(.fit)  # Store intercept and slope in B x 2 matrix
        res$vcov[[b]] <- vcov(.fit)  # Store vcov estimates in length B list.
        res$sigma2[[b]] <- sigma(.fit) ^ 2
    }
    
    return(res)
}

print.lmsim <- function(x) 
{
    cat("Simulation results:\n\n")
    cat("* Model-trusting noise variance:\n ")
    print( mean(x$sigma2) )
    cat("* Model-trusting vcov of coefficient estimates:\n")
    print( avg_est_vcov <- Reduce("+", x$vcov) / length(x$vcov) )
    cat("\n* Simulation-based vcov of coefficient estimates:\n")
    print( emp_vcov <- cov(x$coef))
    cat("\n* Ratio (Model-trusting / Simulation):\n")
    print( avg_est_vcov / emp_vcov )
    return(invisible(x))
}

We simulate linear regression on data generated according to:

where . The noise is , and results in the unconditional variance of the corresponding linear model , twice the conditional variance .

rxy_01 <- rxy_fun(
    rx = \(n) 1 + arima.sim(list(order = c(1,0,0), ar = 0.4), n = n),
    f = \(x) 1 + x,
    reps = \(x) arima.sim(
        list(order = c(1,0,0), ar = 1/sqrt(2)), 
        n = length(x) 
        )
)

plot(rxy_01)

From the simulation below, we see that with serial observations, is relatively close to , but the grossly underestimates all entries (as can be seen from the last line of the output of lmsim() below).

lmsim(rxy_01, N = 100)

Simulation results:

* Model-trusting noise variance:
 [1] 1.870606
* Model-trusting vcov of coefficient estimates:
            (Intercept)           x
(Intercept)  0.03583159 -0.01663739
x           -0.01663739  0.01659708

* Simulation-based vcov of coefficient estimates:
            (Intercept)           x
(Intercept)  0.15486131 -0.02665435
x           -0.02665435  0.02978162

* Ratio (Model-trusting / Simulation):
            (Intercept)         x
(Intercept)   0.2313786 0.6241905
x             0.6241905 0.5572928

To correctly estimate , we could try using the “autocorrelation-consistent” sandwich estimator sandwich::vcovHAC() ⁶. It turns out that, even with a relatively simple example like the present one, the sample size required for the HAC estimator’s bias to die out is unreasonably large (see below). With such large samples, one can probably obtain much better results by leaving out some data for model building, performing inference on the remaining data with a proper time-series model.

lmsim(rxy_01, vcov = sandwich::vcovHAC, N = 100)

Simulation results:

* Model-trusting noise variance:
 [1] 1.870606
* Model-trusting vcov of coefficient estimates:
            (Intercept)           x
(Intercept)  0.08787795 -0.02323242
x           -0.02323242  0.02339146

* Simulation-based vcov of coefficient estimates:
            (Intercept)           x
(Intercept)  0.15486131 -0.02665435
x           -0.02665435  0.02978162

* Ratio (Model-trusting / Simulation):
            (Intercept)         x
(Intercept)   0.5674623 0.8716182
x             0.8716182 0.7854329

lmsim(rxy_01, vcov = sandwich::vcovHAC, N = 500)

Simulation results:

* Model-trusting noise variance:
 [1] 1.974131
* Model-trusting vcov of coefficient estimates:
             (Intercept)            x
(Intercept)  0.023032270 -0.005723968
x           -0.005723968  0.005684149

* Simulation-based vcov of coefficient estimates:
             (Intercept)            x
(Intercept)  0.029600757 -0.005993161
x           -0.005993161  0.006152216

* Ratio (Model-trusting / Simulation):
            (Intercept)         x
(Intercept)   0.7780973 0.9550834
x             0.9550834 0.9239189

lmsim(rxy_01, vcov = sandwich::vcovHAC, N = 1000)

Simulation results:

* Model-trusting noise variance:
 [1] 1.98033
* Model-trusting vcov of coefficient estimates:
             (Intercept)            x
(Intercept)  0.011878079 -0.002771484
x           -0.002771484  0.002849089

* Simulation-based vcov of coefficient estimates:
             (Intercept)            x
(Intercept)  0.015085291 -0.002844716
x           -0.002844716  0.002855522

* Ratio (Model-trusting / Simulation):
            (Intercept)         x
(Intercept)   0.7873948 0.9742566
x             0.9742566 0.9977471

@online{gherardi2023,
  author = {Gherardi, Valerio},
  title = {Linear Regression with Autocorrelated Noise},
  date = {2023-05-25},
  url = {https://vgherard.github.io/posts/2023-05-20-linear-regression-with-autocorrelated-noise/linear-regression-with-autocorrelated-noise.html},
  langid = {en}
}