Using White's Robust Co-variance Matrix vs Weighted Least Squares to correct for heteroscedasticity

Question

I've been trying to figure this out for a bit. How does using White's robust co-variance matrix in OLS vs weighted least squares affect mean response confidence intervals?

I've experimented with both of these and the mean responses are very similar, but my confidence intervals for the mean responses vary quite a bit.

Here is some sample code. I am using python statsmodels library. My features consist of 36 qualitative inputs and two quantitative inputs with 516 observations. I estimated the standard deviation function by regressing the residuals against the fitted values using OLS, then used those as the weights (as laid out by Kutner et al.) I use the get_prediction method of my 'reg' model instance below to get the mean and confidence interval.

Using White's Robust Covariance:

import statsmodels.api as sm
mod = sm.OLS(y, X)
reg = mod.fit(cov_type='HC0')
reg.get_prediction(exog=data).summary_frame()

Using Weighted Least Squares:

mod = sm.WLS(y, X, weights=weights)
reg = mod.fit()
reg.get_prediction(exog=data).summary_frame()

Answer 1

Under the usual Gauss-Markov Model where $\boldsymbol{Y}=\boldsymbol{X\beta}+\boldsymbol{\epsilon}$, it is assummed $\boldsymbol{\epsilon}\sim N\left(0,\boldsymbol{I}\sigma^{2}\right)$, where $\boldsymbol{I}$ is an $n\times n$ identity matrix, where $n$ is the number of observations in your dataset. This implies $Var(\boldsymbol{Y})=Var\left(\boldsymbol{X\beta}+\boldsymbol{\epsilon}\right)=Var(\boldsymbol{\epsilon})=\boldsymbol{I}\sigma^{2}.$ Then, under this model:

\begin{eqnarray*} Var\left(\hat{\boldsymbol{\beta}}\right) & = & Var\left[\boldsymbol{\left(X^{\prime}X\right)^{-}X^{\prime}Y}\right]\\ & = & Var\left[\boldsymbol{AY}\right]\\ & = & \boldsymbol{A}Var(\boldsymbol{Y})\boldsymbol{A^{\prime}}\\ & = & \boldsymbol{A}\left(\boldsymbol{I}\sigma^{2}\right)\boldsymbol{A}^{\prime}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(*)\\ & = & \left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]\left(\boldsymbol{I}\sigma^{2}\right)\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]^{\prime}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\left[\boldsymbol{X}^{\prime}\boldsymbol{X}\right]\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\\ & = & \sigma^{2}\boldsymbol{I}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\\ & = & \sigma^{2}\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-} \end{eqnarray*}

And normally we don't know $\sigma^{2}$ so it is estimated from the residuals to give us: $Var\left(\boldsymbol{\hat{\beta}}\right)=\hat{\sigma}^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]$. Most statistical software will calculate the variance of the estimated coefficients by using this formula. But what happens now, if we use a model where we still have $\boldsymbol{Y}=\boldsymbol{X\beta}+\boldsymbol{\epsilon}$, but this time it is assummed $\boldsymbol{\epsilon}\sim N\left(0,\boldsymbol{\Sigma}\right)$? Well, in this case, the variance formula changes in (*) above to be:

\begin{eqnarray*} Var\left(\hat{\boldsymbol{\beta}}\right) & = & \boldsymbol{A\Sigma}\boldsymbol{A}^{\prime}\\ & = & \left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]\boldsymbol{\Sigma}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]^{\prime}\\ & = & \left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{\Sigma X}\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-} \end{eqnarray*}

So, now, you see that the formula has changed a bit. Since we no longer can assume homoscedastic variances, we have the covariance matrix, $\boldsymbol{\Sigma}$ sandwiched into the middle of the formula for the variance of the estimated coefficients (hence the nickname the 'Huber-White "Sandwich" Estimator'). Again, we typically do not know this, so we usually estimate $\boldsymbol{\Sigma}$ with $\hat{\boldsymbol{\Sigma}}$ from the residuals from an initial fitting by Ordinary Least Squares Regression and then once $\hat{\boldsymbol{\Sigma}}$ is estimated, the $Var\left(\hat{\boldsymbol{\beta}}\right)$ can be estimated. You may obtain the specifics of how each program goes about estimating $\hat{\boldsymbol{\Sigma}}$ exactly by referring to their documentation.

Now, under weighted least squares, the estimate of $\boldsymbol{\beta}$ is now weighted and takes on different form. It is of the form:

\begin{eqnarray*} \hat{\boldsymbol{\beta}}_{weighted} & = & \left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{WY} \end{eqnarray*}

where $\boldsymbol{W}$ is a diagonal matrix of weights (usually formed by taking the inverse of the fitted value from some variance function). The variance of this estimate is different. It is of the form:

\begin{eqnarray*} Var\left(\hat{\boldsymbol{\beta}}_{weighted}\right) & = & Var\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{WY}\right]\\ & = & Var\left[\boldsymbol{AY}\right]\\ & = & \boldsymbol{AY}\boldsymbol{A}^{\prime}\\ & = & \boldsymbol{A}\sigma^{2}\boldsymbol{W}\boldsymbol{^{-}A}^{\prime}\\ & = & \sigma^{2}\boldsymbol{A}\boldsymbol{W}^{-}\boldsymbol{A}^{\prime}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]\boldsymbol{W}^{-}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]^{\prime}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]\boldsymbol{W}^{-}\boldsymbol{W}^{\prime}\boldsymbol{X}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]\boldsymbol{X}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\right]\left[\boldsymbol{X}^{\prime}\boldsymbol{W}\boldsymbol{X}\right]\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\boldsymbol{I}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-} \end{eqnarray*}

and $\sigma^{2}$ is can be estimated here by a number of different methods, including use of replicates or near replicates, or by modelling the variance function. See Chapter 11 of Applied Linear Statistical Models by Kutner, Nachtsheim, Neter, and Li, 5th ed. for additional details.