# Linear model estimation in the presence of heteroscedasticity

Assuming a sample of random variables where the error terms for each random variable ($y_{i}$) are given by $\epsilon_{1}, \dots, \epsilon_{n} \sim N(0, \sigma^{2})$, a linear model is developed such that:

$$y_{i} = \beta_{0} + \epsilon_{i}$$

then, using OLS, finding,

$$\frac{d}{d\hat{\beta}_{0}} \sum_{i}^{n}(y_{i} - \hat{\beta}_{0})^{2} = 0$$

and solving for $\hat{\beta}_{0}$ yields, $\hat{\beta}_{0} = \frac{1}{n}\sum_{i}^{n} y_{i} = \bar{y}$

However, now, let's assume that there are two random variables where the error terms for the corresponding random variable ($y_{i}$) are given by $\epsilon_{1} \sim N(0, \sigma^{2})$ and $\epsilon_{2} \sim N(0, n\sigma^{2})$.

Once again, the a linear model is developed such that:

$$y_{i} = \beta_{0} + \epsilon_{i}$$

Can the estimator for $\beta_{0}$ given above be used for this model? If not, in light of the heteroscedasticity involved in the error terms, how is an estimator for $\beta_{0}$ developed?

I've seen linear models that, given a data point, use indicator functions to determine which terms in the model are non-zero. Perhaps this applies here?

If you know which distribution your samples are coming from, you can use weighted least squares. In particular,

The Gauss–Markov theorem shows that, when this is so, $\hat{\boldsymbol{\beta}}$ is a best linear unbiased estimator (BLUE). If, however, the measurements are uncorrelated but have different uncertainties, a modified approach might be adopted. Alexander Aitken showed that when a weighted sum of squared residuals is minimized, $\hat{\boldsymbol{\beta}}$ is BLUE if each weight is equal to the reciprocal of the variance of the measurement.

so give weights $w_i=\frac{1}{{\sigma_i}^2}$.

I've usually heard weighted least squares for heteroskedasticity motivated as user1448319 quoted -- as a solution for predictable sampling variance.

But I'm more accustomed to simply running OLS and adjusting the standard errors by calculating the variance-covariance matrix as

$$(X'X)^{-1}X'\epsilon \epsilon' X (X'X)^{-1}$$

...which gives consistent standard errors in the presence of arbitrary forms of heteroskedasticity. The point estimates will be identical, because the vcov is calculated from the fitted OLS regression.h

Alternatively, you can use the wild bootstrap (or the wild cluster bootstrap for repeated-observation or panel data).

A fairly readably explaination is here

I'll expand:

First off, you will want to notice that the expected value of the error term is still sufficient to achieve an unbiased estimate. That means - in expectation - your estimate is still correct.
So the answer to your question is yes, you can still use your estimator to estimate the parameter.

What is not correct anymore, however, is the estimate of the standard error and with that the estimate of the variance of the estimator $\beta_0$.
Intuitively: If there are two different error term variances in your data. Any single "estimator" you will try to find will therefore lie between those true values. It will be wrong. The trouble is that you need this estimate of $V(\hat{\beta})$ to calculate, for example, the interval estimations or do hypothesis testing. And in general to know how accurate your model is. Without a correct estimate, your estimator may still be correct in expectation, but there is no way you may know how far off it is in reality.

There are few ways around this though. In your case if you can segment the data into the two different sets of data you can use the weighted least squares option. This assumes that you can assign weights to all observations which correct for the error, which is assumed to be different for all observations but still uncorrelated with each other (so no correlation between the errors of observations - no autocorrelation).
Mathematically it says the covariance matrix is a product of your $\sigma^2$ and a diagonal matrix so that $V[\epsilon] = \sigma^2\Omega$. In this case this matrix would include either the term $1$ or $n$ depending on the $i$.
In this case all you need to do is to estimate your model as you would otherwise, but use a different formula . This is easiest to do in matrix notation. In essence, you transform your model with a matrix which look like this $$W = \begin{pmatrix} w_1 & 0 & 0 \\ 0 & w_2 & 0 \\ 0 & 0 & w_i \end{pmatrix} = \Omega^{-1}$$ where the $w_i$ are the weighting factors which transform the standard error of each observation into a homoscedastic error. For the ones which are have the error of the first type, the $w_i$ is equal to $1/\sigma^2$. For the other ones you use $w_i=1/(n\sigma^2)$.
Next, you plug this matrix in the WLS estimator $$\hat{\beta} = (X'WX)^{-1}X'Wy$$ Your estimator is now distributed as such: $\widehat \beta_{WLS} ~\sim N(\beta , \sigma^2(X'W^{-1}X)^{-1})$ and this is homoscedastic. As such, this is then the BLUE estimator - all you need.
But if you are just going to point estimate your $\beta$, you will notice that this may even lead to the same estimate as your original estimator without using $W$ (or not). What you need this transformation for, as hinted above, are the variance estimation.
To dig a bit deeper, WLS is actually just a special case of Generalized Least Squares. With GLS you can also get rid of autocorrelation by using the off-diagonal terms in your matrix $W$. What GLS does is to transform the model, multiplying it with a matrix $P$ which is chosen so that $P'P = W$ (don't worry about how to get $P$ exactly). As it turns out if you multiply your model by $P$, you can then use OLS as usual on THAT data and you will have the unbiased estimate you originally wanted. Anyway this allows for a quick and easy estimate of our variance of beta: $$V[\hat{\beta}] = \sigma^2(X'WX)^{-1}$$ Since you assume to know all these values, you can calculate the standard errors with this.

In reality, your $\Omega$ may be unknown. This is when you can use Feasible generalized least squares. It works almost the same only that you estimate your missing values first. You do an OLS regression and just use the residuals as estimators for your Matrix.

But really, in that case it would be wise to use error terms which are already consistent. You can use the vcov command in R for this. The so called "White HC" error terms are consistent under heteroscedasticity. If your sample isn't too small and if you don't know the exact structure of the heteroscedasticity, this is a good bet.