Under the usual 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.

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.

Under the usual 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.

Under the usual 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.