Why would bootstrap OLS standard errors differ from ML estimate?

Question

Let's say I have a regression dataset (paired x and y) such that the response variable (y) has an unknown distribution (but definitely not Gaussian) and is large enough such that the central limit theorem holds. If I do OLS with bootstrapping and observe that the bootstrap standard error of the coefficient is far from the ML estimate of the standard error, does that indicate that Gauss-Markov assumptions are violated? If so, let's then assume that homoskedasticity is violated and thus the standard error estimate is not BLUE. If I use an approach like adjusting the OLS with robust standard errors, should you expect both standard error estimates to agree?

Answer 1

You might just as well ask, "Why is the sample correlation between two independent variables always non-zero?" and the answer is the same: samples never conserve the properties of the probability models that generate them unless they are, somehow, perfectly balanced.

Consider this example where I generate data that have errors generated straight from empirical normal quantiles. That is that the ECDF of the residual very closely matches the normal distribution.

library(boot)
set.seed(123)
x <- seq(-3, 3, 0.1)
e <- qnorm(1:99/100)
data <- expand.grid('x'=x,'e'=e)
data$y <- data$x + data$e
f <- lm(y ~ x, data=data)
confint(f)
betahat <- function(data,i) cov(data[i, 'x'],data[i, 'y'])/var(data[i, 'x'])
b <- boot(data, betahat, R=1000, stype='i')
boot.ci(b, type='norm')

Gives:

> confint(f)
                  2.5 %     97.5 %
(Intercept) -0.02422861 0.02422861
x            0.98623908 1.01376092
> boot.ci(b, type='norm')
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = b, type = "norm")

Intervals : 
Level      Normal        
95%   ( 0.9861,  1.0140 )  
Calculations and Intervals on Original Scale

Which agrees out to nearly 4 decimal places. The asymptotic result is that these are the same, but in finite samples, you'll never have perfectly normal residuals.

It's true that Gauss-Markov does not require that errors are normally distributed. It is an asymptotic result. But the question bears consideration: if the bootstrapped error estimate is very different from the model-based error estimate, should you go with the bootstrapped error estimate because the probability model is misspecified? It can go both ways. Change the above to a simulation to consider the 80% coverage of the CIs.

library(boot)
set.seed(123)
`%has%` <- function(ci, true) true > ci[1] & true < ci[2]    

res <- replicate(1000, {
  x <- seq(-3, 3, 1)
  e <- rnorm(7)
  data <- data.frame(x=x, e=e)
  data$y <- data$x + data$e
  f <- lm(y ~ x, data=data)
  
  betahat <- function(data,i) cov(data[i, 'x'],data[i, 'y'])/var(data[i, 'x'])
  b <- boot(data, betahat, R=1000, stype='i')
  
  c(
    olscover = confint(f, level = 0.8, 'x') %has% 1,
    bscover = boot.ci(b, type='norm', conf=0.8)$normal[-1] %has% 1
  )

})

Gives:

> rowMeans(res)
olscover  bscover 
   0.795    0.782 

In other words, the bootstrap fails to produce CIs that cover at the nominal 80% rate. The difference is not substantial, and the simulation isn't exactly fast, but you can play with it to understand the risks. The bootstrap is just slower to converge, so in small sample like this N=7 regression, when the OLS assumptions hold, the OLS is the better answer. On the other hand, you don't know when those assumptions are true, and if you have decent power, the bootstrap is a robust safeguard against model misspecifications.