Using White's Robust Co-variance Matrix vs Weighted Least Squares to correct for heteroscedasticity 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()

 A: 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. 
