Why does heteroskedasticity not affect $R^2$ and why does it make estimated regression more statistically significant?

Question

I am studying what the consequences of heteroskedasticity are. And i found that assuming that the model is linear in the parameters (i.e $Y=X\beta+\epsilon$), is identifiable, has no perfect collinearity (i.e $X$ is full column rank), is exogenous (i.e. $E[\epsilon|X]=0$) and has no serial correlations (i.e. $Var[\epsilon|X]=\Sigma$, where $\Sigma$ is diagonal) we have :

1. $E[\hat{\beta}] = \beta$, because $\hat{\beta} = \beta + (X^TX)^{-1}X^T\epsilon$ and $E[\epsilon|X]=0$
   

2. OLS is no longer BLUE because the homoskedasticity assumption is absent (and hence Gauss-Markov theorem cannot be applied).
   

3. Estimated variances of regression coefficients will be biased (because $E\bigg[\hat{Var}[\hat{\beta}_i|X]\bigg|X\bigg] = (X^TX)^{-1}_{ii}E\bigg[MS_{res}\bigg|X\bigg]\neq (X^TX)^{-1}_{ii}\Sigma_{ii}=Var[\hat{\beta}_i|X]$, where $MS_{res} = \dfrac{e^Te}{n-p}$, where $e = Y-\hat{Y}=$ the residuals)
   

4. The corresponding t-statistics of the estimated regression coefficients will appear to be more significant than they really are (page 4). This implies that $MS_{res}<\Sigma_{ii}$. I do not understand why this should be the case.
   

5. Since $R^2=1-\dfrac{\epsilon^T(I-H)\epsilon}{Y^T(I-J_n)Y}$ is based on overall sums of squares, it is unaffected by heteroskedasticity (page 5). Again I do not understand why this should be the case. Specifically, in what sense is $R^2$ unaffected? Like if I somehow made a different linear model to deal with the heteroskedasticity, then would this new linear model have the same $R^2$? If so how?
   

Could you please verify points 1,2 and 3, and give non-heuristic explanations for points 4 and 5?

Edit1:

@ChristophHanck,
I dont understand why $R^2$ would converge at all under heteroskedasticity. I mean, under homoskedasticity, assuming $e_i\stackrel{iid}{\sim}(\text{mean}=0,\text{variance}= \sigma^2)$ and $n>>p$, we have $\dfrac{e^Te}{n-p}\approx\dfrac{e^Te}{n}=\dfrac{e_1^2+\cdots+e_n^2}{n}\rightarrow E[\epsilon_i^2]=\sigma^2$, by weak law of large numbers. But under heteroskedasticity, assuming $e_i\stackrel{independent}{\sim}(\text{mean}=0,\text{variance}= \sigma^2_i)$ and $n>>p$, why should $\dfrac{e_1^2+\cdots+e_n^2}{n}$ be expected to converge at all?
And in terms of this answer to your cited question, why should $\dfrac{\sum_{i=1}^N{m_{ii}\sigma^2_i}}{N - K}$ be expected to converge at all?
also, under heterskedasticity, the $y_i$'s are not iid. So, why should $\sum(y_i-\bar{y})^2$ converge at all?

Answer 1

As to 5, $R^2$ will tend to $$ 1-\frac{Var(u)}{Var(y)} $$ under homoskedasticity. Under heteroskedasticity (no detailed derivation here), the numerator will be related to some function of the average variance (by suitably modifying the derivation here: How do I use the standard regression assumptions to prove that $\hat{\sigma}^2$ is an unbiased estimator of $\sigma^2$?)

In a linear model, we have, under suitable assumptions, $$ Var(y)=\beta^2Var(x)+Var(u), $$ so that $$ R^2\to_p1-\frac{Var(u)}{\beta^2Var(x)+Var(u)} $$ So provided the average variance under heteroskedasticity is similar to the variance under homosekdasticity, we may expect similar results for $R^2$.

Illustration:

n <- 10000

u.hom <- rnorm(n)
u.het <- rnorm(n, sd=runif(n, 0.5, 1.5)) #  i.e., on average, the s.d. will also be one
x <- rnorm(n)

y.hom <- 2*x + u.hom
y.het <- 2*x + u.het

> summary(lm(y.hom~x))

Call:
lm(formula = y.hom ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5527 -0.6754 -0.0023  0.6854  3.8040 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.008325   0.010040  -0.829    0.407    
x            2.009258   0.009825 204.494   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.004 on 9998 degrees of freedom
Multiple R-squared:  0.807, Adjusted R-squared:  0.807 
F-statistic: 4.182e+04 on 1 and 9998 DF,  p-value: < 2.2e-16


> summary(lm(y.het~x))

Call:
lm(formula = y.het ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0801 -0.6242 -0.0082  0.6446  4.9125 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.003788   0.010510    0.36    0.719    
x           2.011002   0.010286  195.51   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.051 on 9998 degrees of freedom
Multiple R-squared:  0.7927,    Adjusted R-squared:  0.7927 
F-statistic: 3.823e+04 on 1 and 9998 DF,  p-value: < 2.2e-16

Answer 2

As to 4, that is indeed, while empirically often the case, not necessarily true.

Consider as a tractable example the expressions for the standard and robust variance estimator (the correct one under heteroskedasticity) for a regression on a constant and a dummy investigated in these answers: Eicker-Huber-White Robust Variance Estimator and How to prove equality of standard errors for two-sample t-test and linear regression?

In particular, for $n_1$ "treated" observations for which the dummy equals one and the rest of $n$ being the other $n_0$ observations, we have for the robust variance estimator of the slope coefficient that $$ \widehat{\mathbb{V}}_{\text{EHW}}=\frac{n_1^2s_0^2+n_0^2s_1^2}{n_0^2n_1^2}, $$ with $s_t^2$, $t=0,1$, the sum of squared residuals belonging to the two groups of observations, and $s^2=s_0^2+s_1^2$.

Also (here, we work with the ML estimator for the error variance for simplicity which divides the sum of squared residuals with $n$ instead of $n-2$ to ease the algebra; that does not matter for large $n$; also note the definition of $s^2$ in the second link is different from the one employed here, just scale by $n$), the standard OLS slope variance estimator is $$ \widehat{\mathbb{V}}_{\text{standard}}=\frac{s_0^2+s_1^2}{n_0n_1}, $$

We have $$ \widehat{\mathbb{V}}_{\text{EHW}}<\widehat{\mathbb{V}}_{\text{standard}} $$ when (unless $n_0= n_1$, in which case we always have $\widehat{\mathbb{V}}_{\text{EHW}}=\widehat{\mathbb{V}}_{\text{standard}}$) $$ s_1^2>\frac{n_1}{n_0}s_0^2 $$ or $$ \frac{s_1^2}{n_1}>\frac{s_0^2}{n_0}, $$ i.e., when the variance of the "treated" observations is higher than that of the untreated.

Take, e.g., $n_0=n/4$ and hence $n_1=3n/4$. Then, $$\widehat{\mathbb{V}}_{\text{EHW}}<\widehat{\mathbb{V}}_{\text{standard}}$$ will result when $$ s_1^2>3s_0^2, $$ i.e., when the variance associated with the dummy's observations equal to 1 is sufficiently larger that the others, unless there are relatively too many such observations (i.e., when $n_1/n_0$ is large). Illustration:

library(sandwich)
library(lmtest)

n <- 2000
n0 <- 1/4*n
n1 <- 3/4*n

y <- c(rnorm(n0, sd=1), rnorm(n1, sd=1.05))  # i.e. more variance when dummy=1
x1 <- c(rep(0, n0), rep(1,n1)) # a dummy regressor

reg <- lm(y~x1)

> summary(reg)

Call:
lm(formula = y ~ x1)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3243 -0.6771 -0.0006  0.6766  3.0950 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03402    0.04677  -0.727    0.467
x1           0.02879    0.05400   0.533    0.594

Residual standard error: 1.046 on 1998 degrees of freedom
Multiple R-squared:  0.0001422, Adjusted R-squared:  -0.0003582 
F-statistic: 0.2842 on 1 and 1998 DF,  p-value: 0.594


> robust.se <- coeftest(reg, vcov = vcovHC(reg, "HC0")) 

> robust.se

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.034021   0.045556 -0.7468   0.4553
x1           0.028788   0.053065  0.5425   0.5875


> residuals <- resid(reg)

> 3*sum(residuals[1:n0]^2) # n_1/n_0*s_0^2
[1] 1556.493

> sum(residuals[(n0+1):n]^2) # s_1^2
[1] 1666.235

Answer 3

Interesting question! @christoph hanck has treated a lot of stuff here. In your edit, you ask a question: why does the average of $e_i^2$, with iid $e_i\sim(0,\sigma_i^2)$ converge? Here I provide an answer to this question. First we have \begin{align} \frac 1n \sum_i e_i^2 &= \frac 1n \sum_i (e_i^2-\sigma_i^2) + \frac 1n \sum_i \sigma_i^2 \end{align} Let's look at the first term. Given the indepence assumed on $e_i$, $z_i=e_i^2-\sigma_i^2$ is also indepndent and has a zero mean. Therefore as long as you have a moment of $2+\delta$ for $e_i$, the result (that $\frac 1n \sum_i e_i^2$ converges) would follow. For simplicity, however, let's assume finite fourth moments for $e_i$ such that we could invoke elementary results here; for some real $a>0$, we have \begin{align} \Pr\left(|\bar z|>a\right)&\leq \frac{var(\bar z)}{a^2}\\ &= \frac{1}{n^2a^2}\sum_ivar(z_i)\\ &\leq \frac{1}{n^2a^2}nc\\ &=\frac c{a^2}n^{-1}\end{align} for some constant c. This in turn implies that $\frac 1n \sum_i (e_i^2-\sigma_i^2) \overset{p}{\to}0$. Therefore, $$\frac 1n \sum_i e_i^2\overset{p}{\to} \lim_{n\to \infty}\frac 1n \sum_i \sigma_i^2.$$

Note, however, that in a regression framework and when the focus is the OLS estimators, one focuses often on the standard errors of the OLS estimator and that would be slightly different.