Linear regression with autocorrelated noise
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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 1.
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 is2:
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 biased3. In realistic cases, one expects the correlations to decay exponentially with 4 , 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 biased5.
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.
Code
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() 6. 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
Footnotes
For the linear model assumptions to hold, the
pairs should come from independent realizations of the same time series, which is of course not the type of data we are usually presented with.↩︎As usual we stack observations vertically in the
and matrices.↩︎For an extreme case, suppose that
(no covariate except for an intercept term), and let the noise term be , where and are independent -scores. One can easily see that, in this setting, and .↩︎For instance, for the
noise of Eq. Equation 2, we have .↩︎The difference
decays as , which is of the same order of the estimation target . Not sure I’m using standard terminology here.↩︎Disclaimer: I haven’t read any theory about the HAC estimator, so I may be misusing it here, but I would have expected it to work relatively well on such an “easy” example. For illustrations on how to use sandwich estimators for first- and second-order linear model misspecification, you can read this post of mine.↩︎
Reuse
Citation
@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}
}
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.