Heteroskedastyczność — diagnoza, konsekwencje i leczenie

Średniozaawansowany
Streszczenie

Heteroskedastyczność od podstaw: definicja i źródła, dowód, że estymator MNK pozostaje nieobciążony lecz traci efektywność i ma błędne błędy standardowe (wariancja kanapkowa), test Breuscha-Pagana i White'a, odporne błędy HC0/HC3 z dowodem zgodności, uogólniona i ważona MNK (GLS/WLS) jako BLUE — z wykresami w skali szarości, szczegółowo policzonymi przykładami i kodem R oraz Python.

Założenie sferyczności składnika losowego — stała wariancja i brak korelacji — jest jednym z filarów klasycznego modelu MNK. Heteroskedastyczność narusza jego pierwszą część: wariancja błędu przestaje być stała. Rozważymy, co dokładnie z tego wynika dla estymatora i wnioskowania, jak naruszenie wykryć oraz jak je skorygować — z pełnymi wyprowadzeniami i rachunkiem na konkretnych danych.

Definicja i źródła

Definicja
Homo- i heteroskedastyczność

Model spełnia założenie homoskedastyczności, gdy warunkowa wariancja błędu jest stała:

$$ \mathrm{Var}(\varepsilon_i \mid \mathbf{x}_i) = \sigma^2 \quad\text{dla wszystkich } i. $$

Heteroskedastyczność to sytuacja, w której wariancja zależy od obserwacji (zwykle od zmiennych objaśniających):

$$ \mathrm{Var}(\varepsilon_i \mid \mathbf{x}_i) = \sigma_i^2,\qquad \sigma_i^2 \neq \text{const}. $$

W zapisie macierzowym $\mathrm{Var}(\boldsymbol\varepsilon\mid\mathbf{X})=\boldsymbol\Omega=\operatorname{diag}(\sigma_1^2,\dots,\sigma_n^2)$ zamiast $\sigma^2\mathbf{I}$.

Typowym przykładem jest regresja dochodu na wykształcenie: wśród osób słabo wykształconych dochody są zbliżone (mała wariancja), wśród wysoko wykształconych rozpiętość jest ogromna — od nisko opłacanego absolwenta po przedsiębiorcę (duża wariancja). Heteroskedastyczność pojawia się najczęściej w danych przekrojowych o dużym zróżnicowaniu jednostek (firmy, kraje), przy zmiennych ilorazowych (dochód, PKB, cena) oraz w danych panelowych.

Reszty względem dopasowanych — homoskedastyczność a lejek heteroskedastyczności
Najważniejsze narzędzie diagnostyczne: wykres reszt względem wartości dopasowanych. Po lewej reszty mają stały rozrzut wokół zera (homoskedastyczność). Po prawej rozrzut rośnie wraz z $\hat y$ — charakterystyczny lejek, sygnał heteroskedastyczności.

Konsekwencje dla estymatora MNK

Kluczowe jest rozróżnienie dwóch własności estymatora: nieobciążoności (która przeżywa heteroskedastyczność) oraz efektywności i poprawności błędów standardowych (które nie).

Twierdzenie
Skutki heteroskedastyczności

Przy heteroskedastyczności estymator MNK $\hat{\boldsymbol\beta}$ pozostaje nieobciążony i zgodny, lecz:

  1. traci efektywność — nie jest już najlepszym liniowym nieobciążonym estymatorem (BLUE);
  2. jego prawdziwa wariancja to macierz kanapkowa $(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\Omega\,\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}$, a nie $\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}$;
  3. klasyczny estymator błędów standardowych jest więc niezgodny, co psuje statystyki $t$, $F$ i przedziały ufności.
Dowód
Nieobciążoność przeżywa, lecz wariancja jest kanapkowa
  1. Nieobciążoność nie korzysta z homoskedastyczności. Podstawiając $\mathbf{y}=\mathbf{X}\boldsymbol\beta+\boldsymbol\varepsilon$, mamy $\hat{\boldsymbol\beta}-\boldsymbol\beta=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\varepsilon$. Stąd $$ \mathbb{E}[\hat{\boldsymbol\beta}\mid\mathbf{X}]=\boldsymbol\beta+(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\underbrace{\mathbb{E}[\boldsymbol\varepsilon\mid\mathbf{X}]}_{=\,\mathbf{0}}=\boldsymbol\beta. $$ Wykorzystano jedynie egzogeniczność (A3), nie zaś sferyczność (A4) — nieobciążoność jest więc nienaruszona.
  2. Prawdziwa wariancja. Z tej samej tożsamości, podstawiając $\mathrm{Var}(\boldsymbol\varepsilon\mid\mathbf{X})=\boldsymbol\Omega$: $$ \mathrm{Var}(\hat{\boldsymbol\beta}\mid\mathbf{X})=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\Omega\,\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}. $$
  3. Kiedy upraszcza się do wzoru klasycznego. Tylko przy $\boldsymbol\Omega=\sigma^2\mathbf{I}$ środkowy czynnik daje $\sigma^2\mathbf{X}^\top\mathbf{X}$ i całość redukuje się do $\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}$. Przy heteroskedastyczności człony się nie skracają, więc estymator $s^2(\mathbf{X}^\top\mathbf{X})^{-1}$ szacuje **niewłaściwą** macierz — może być zarówno zaniżony, jak i zawyżony.
  4. Utrata efektywności. Dowód [twierdzenia Gaussa-Markowa](/ekonometria/zalozenia-kmnk/) wymaga $\boldsymbol\Omega=\sigma^2\mathbf{I}$; bez tego MNK przestaje mieć najmniejszą wariancję w klasie liniowych nieobciążonych estymatorów. Optymalny estymator podaje GLS (poniżej).
∎

Wykrywanie heteroskedastyczności

Punktem wyjścia jest wykres reszt; formalne potwierdzenie dają testy oparte na regresji pomocniczej kwadratów reszt.

Twierdzenie
Test Breuscha-Pagana i White'a

Breusch-Pagan. Estymujemy regresję pomocniczą kwadratów reszt na zmienne objaśniające:

$$ \hat\varepsilon_i^2=\delta_0+\delta_1 x_{i1}+\dots+\delta_p x_{ip}+u_i. $$

Statystyka mnożnika Lagrange’a $\mathrm{LM}=n\cdot R^2_{\text{aux}}\sim\chi^2_p$ przy $H_0$ (homoskedastyczność, $\delta_1=\dots=\delta_p=0$).

White. Ta sama konstrukcja, lecz regresja pomocnicza zawiera dodatkowo kwadraty i iloczyny zmiennych — wychwytuje również heteroskedastyczność nieliniową, kosztem mocy w małej próbie.

Regresja kwadratów reszt na x — test Breuscha-Pagana
Test Breuscha-Pagana na zbiorze z lejkiem: kwadraty reszt $\hat\varepsilon_i^2$ wyraźnie rosną ze zmienną $x$, a dopasowana prosta pomocnicza ma dodatnie, istotne nachylenie. Im lepiej regresja pomocnicza tłumaczy $\hat\varepsilon^2$ (większe $R^2_{\text{aux}}$), tym większe $\mathrm{LM}=nR^2_{\text{aux}}$.
Przykład
Breusch-Pagan: dane jednorodne

Dla regresji MNK $n=5$ reszty wynoszą $\hat\varepsilon=(-0{,}8;\ 0{,}6;\ 1{,}0;\ -0{,}6;\ -0{,}2)$, więc $\hat\varepsilon^2=(0{,}64;\ 0{,}36;\ 1{,}00;\ 0{,}36;\ 0{,}04)$. Regresja pomocnicza $\hat\varepsilon^2$ na $x$ daje $R^2_{\text{aux}}=0{,}278$, zatem

$$ \mathrm{LM}=n\cdot R^2_{\text{aux}}=5\cdot 0{,}278=1{,}39,\qquad \chi^2_{1;\,0{,}05}=3{,}84. $$

Ponieważ $1{,}39<3{,}84$, nie odrzucamy homoskedastyczności — dane nie wykazują heteroskedastyczności.

Przykład
Breusch-Pagan: dane z lejkiem

Rozważmy zbiór $x=(1,\dots,8)$, $y=(3,5,6,11,9,16,12,22)$. Regresja $y$ na $x$ daje $\hat y=0+\tfrac{7}{3}x$, reszty o rosnącym rozrzucie i kwadraty $\hat\varepsilon^2=(0{,}44;\ 0{,}11;\ 1{,}00;\ 2{,}78;\ 7{,}11;\ 4{,}00;\ 18{,}78;\ 11{,}11)$. Regresja pomocnicza $\hat\varepsilon^2$ na $x$ ma nachylenie $2{,}16$ i $R^2_{\text{aux}}=0{,}664$, więc

$$ \mathrm{LM}=8\cdot 0{,}664=5{,}31,\qquad \chi^2_{1;\,0{,}05}=3{,}84. $$

Ponieważ $5{,}31>3{,}84$ ($p=0{,}021$), odrzucamy homoskedastyczność — heteroskedastyczność jest istotna.

Pomocny jest też wykres skala–położenie ($\sqrt{|\hat\varepsilon|}$ względem $\hat y$): rosnący trend potwierdza, że rozrzut zależy od poziomu prognozy.

Wykres skala-położenie z rosnącym trendem
Wykres skala–położenie dla zbioru z lejkiem: pierwiastek z wartości bezwzględnej reszty rośnie z $\hat y$. Linia trendu o dodatnim nachyleniu jest graficznym odpowiednikiem dodatniego współczynnika w regresji pomocniczej testu BP.

Strukturę wariancji rosnącej z wieloma regresorami można zobrazować jako powierzchnię nad płaszczyzną zmiennych objaśniających.

Powierzchnia wariancji błędu rosnąca ze zmiennymi objaśniającymi
Wariancja błędu $\sigma_i^2$ jako rosnąca powierzchnia nad płaszczyzną dwóch regresorów. Przy homoskedastyczności powierzchnia byłaby płaska na wysokości $\sigma^2$; tutaj rośnie ku narożnikom dużych wartości $x_1,x_2$ — to właśnie kształt, który testy BP i White’a starają się wykryć.

Leczenie I: odporne błędy standardowe

Najprostsze i najpowszechniejsze podejście nie zmienia estymatora $\hat{\boldsymbol\beta}$, lecz zastępuje błędną macierz wariancji jej odpornym estymatorem, zgodnym nawet przy nieznanej heteroskedastyczności.

Twierdzenie
Estymator White'a (HC)

Macierz kanapkową szacuje się, zastępując nieznane $\sigma_i^2$ przez $\hat\varepsilon_i^2$:

$$ \widehat{\mathrm{Var}}_{\text{HC}}(\hat{\boldsymbol\beta})=(\mathbf{X}^\top\mathbf{X})^{-1}\!\left(\sum_{i} w_i\,\hat\varepsilon_i^2\,\mathbf{x}_i\mathbf{x}_i^\top\right)\!(\mathbf{X}^\top\mathbf{X})^{-1}, $$

gdzie $w_i=1$ daje HC0 (oryginał White’a), a $w_i=(1-h_{ii})^{-2}$ — HC3, zalecany w małych próbach ($h_{ii}$ to elementy diagonalne macierzy kapeluszowej).

Dowód
Dlaczego White działa, mimo że σ²ᵢ są nieznane
  1. Cel. Wariancja kanapkowa wymaga środkowego czynnika $\mathbf{X}^\top\boldsymbol\Omega\mathbf{X}=\sum_i\sigma_i^2\,\mathbf{x}_i\mathbf{x}_i^\top$. Pojedynczych $\sigma_i^2$ nie da się estymować (jedna obserwacja na parametr).
  2. Kluczowe spostrzeżenie. Nie potrzebujemy poszczególnych $\sigma_i^2$, lecz całej **sumy** $\sum_i\sigma_i^2\mathbf{x}_i\mathbf{x}_i^\top$ — obiektu o stałym wymiarze $k\times k$. Ponieważ $\mathbb{E}[\hat\varepsilon_i^2\mid\mathbf{x}_i]\approx\sigma_i^2$, średnia $\tfrac1n\sum_i\hat\varepsilon_i^2\mathbf{x}_i\mathbf{x}_i^\top$ jest zgodnym estymatorem $\tfrac1n\sum_i\sigma_i^2\mathbf{x}_i\mathbf{x}_i^\top$ (prawo wielkich liczb dla $k\times k$ wyrazów).
  3. Korekta HC3. Reszty MNK są przeciętnie „za małe" o czynnik $(1-h_{ii})$ (rzut zmniejsza ich wariancję: $\mathbb{E}[\hat\varepsilon_i^2]=(1-h_{ii})\sigma_i^2$ przy homoskedastyczności). Wagi $(1-h_{ii})^{-2}$ kompensują to obciążenie, dając lepsze własności w małej próbie.
  4. Wniosek. Wstawienie zgodnego estymatora środka w miejsce $\mathbf{X}^\top\boldsymbol\Omega\mathbf{X}$ daje zgodny estymator całej macierzy wariancji — bez żadnego założenia o kształcie heteroskedastyczności.
∎
Przykład
Błędy klasyczne a HC3

Dla zbioru z lejkiem ($\hat\beta_1=\tfrac73=2{,}333$) porównajmy błąd standardowy nachylenia:

Estymator wariancji$\operatorname{SE}(\hat\beta_1)$$t=\hat\beta_1/\operatorname{SE}$
klasyczny $s^2(\mathbf{X}^\top\mathbf{X})^{-1}$$0{,}424$$5{,}50$
White HC0$0{,}394$$5{,}93$
HC3 (zalecany)$0{,}612$$3{,}81$

Klasyczny błąd standardowy zaniża niepewność, zawyżając statystykę $t$ z $3{,}81$ (HC3) do $5{,}50$. Wniosek o istotności pozostaje, lecz przy słabszym efekcie klasyczne SE mogłyby dać fałszywą istotność. HC3 jest tu wyraźnie ostrożniejszy niż HC0 — stąd zalecenie dla małych prób.

Leczenie II: uogólniona i ważona MNK

Gdy znamy (lub potrafimy oszacować) strukturę wariancji, możemy odzyskać efektywność, przekształcając model tak, by błąd znów był sferyczny.

Twierdzenie
GLS Aitkena i WLS

Dla $\mathrm{Var}(\boldsymbol\varepsilon\mid\mathbf{X})=\boldsymbol\Omega$ estymator uogólnionej MNK (GLS)

$$ \hat{\boldsymbol\beta}_{\text{GLS}}=(\mathbf{X}^\top\boldsymbol\Omega^{-1}\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\Omega^{-1}\mathbf{y} $$

jest BLUE. Przy heteroskedastyczności $\boldsymbol\Omega$ jest diagonalna, więc GLS sprowadza się do ważonej MNK (WLS) z wagami $w_i=1/\sigma_i^2$ — minimalizacji $\sum_i w_i(y_i-\mathbf{x}_i^\top\boldsymbol\beta)^2$.

Dowód
GLS jest BLUE — przez sprowadzenie do modelu sferycznego
  1. Przekształcenie. Niech $\mathbf{P}=\boldsymbol\Omega^{-1/2}$ (istnieje, bo $\boldsymbol\Omega$ jest dodatnio określona). Mnożymy model lewostronnie: $\mathbf{P}\mathbf{y}=\mathbf{P}\mathbf{X}\boldsymbol\beta+\mathbf{P}\boldsymbol\varepsilon$.
  2. Błąd staje się sferyczny. $\mathrm{Var}(\mathbf{P}\boldsymbol\varepsilon)=\mathbf{P}\,\boldsymbol\Omega\,\mathbf{P}^\top=\boldsymbol\Omega^{-1/2}\boldsymbol\Omega\,\boldsymbol\Omega^{-1/2}=\mathbf{I}$. Model przekształcony spełnia więc założenia klasyczne.
  3. MNK na przekształconym modelu jest BLUE (twierdzenie Gaussa-Markowa). Po podstawieniu: $$ \hat{\boldsymbol\beta}=\big((\mathbf{P}\mathbf{X})^\top\mathbf{P}\mathbf{X}\big)^{-1}(\mathbf{P}\mathbf{X})^\top\mathbf{P}\mathbf{y}=(\mathbf{X}^\top\boldsymbol\Omega^{-1}\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\Omega^{-1}\mathbf{y}. $$
  4. Przypadek diagonalny. Dla $\boldsymbol\Omega=\operatorname{diag}(\sigma_i^2)$ mamy $\boldsymbol\Omega^{-1}=\operatorname{diag}(1/\sigma_i^2)$, więc formuła to dokładnie ważona MNK z wagami $w_i=1/\sigma_i^2$ — obserwacjom o większej wariancji nadajemy mniejszą wagę.
∎
Przykład
WLS na zbiorze z lejkiem

Przyjmijmy, że wariancja rośnie liniowo, $\sigma_i^2\propto x_i$, więc wagi $w_i=1/x_i$. Rozwiązując ważone równania normalne $\mathbf{X}^\top\mathbf{W}\mathbf{X}\hat{\boldsymbol\beta}=\mathbf{X}^\top\mathbf{W}\mathbf{y}$ z $\mathbf{W}=\operatorname{diag}(1/x_i)$ dla danych z przykładu 2, otrzymujemy

$$ \hat{\boldsymbol\beta}_{\text{WLS}}=(0{,}547;\ 2{,}212), $$

wobec MNK $\hat{\boldsymbol\beta}=(0;\ 2{,}333)$. WLS ufa bardziej obserwacjom o małym $x$ (mniejsza wariancja), nieco zmieniając oba współczynniki; jest przy tym efektywniejszy od MNK, jeśli założenie o wagach jest trafne. W praktyce $\sigma_i^2$ estymuje się z reszt (FGLS), gdyż rzadko są znane.

Kiedy co stosować

SytuacjaZalecenie
Brak pewności co do problemudomyślnie HC (HC3) — bezpieczne nawet bez heteroskedastyczności
Potwierdzona heteroskedastyczność, nieznana strukturaHC3
Znana lub dobrze oszacowana struktura $\sigma_i^2$WLS / FGLS (efektywniejszy)
Dane prawoskośne, $y>0$transformacja logarytmiczna $\ln y$
Dane panelowebłędy klasterowane po jednostkach

Odporne błędy nie naprawiają obciążenia.

Odporne błędy standardowe korygują wyłącznie wnioskowanie (statystyki $t$, $F$, przedziały). Nie usuwają obciążenia z pominiętej zmiennej ani błędnej formy funkcyjnej — te wymagają zmiany specyfikacji modelu, a nie korekty macierzy wariancji. Heteroskedastyczność to problem efektywności i błędów standardowych, nie poprawności punktowych oszacowań.

Podsumowanie

Definicja
Heteroskedastyczność w pigułce
  • Definicja: $\mathrm{Var}(\varepsilon_i\mid\mathbf{x}_i)=\sigma_i^2$ zmienne — naruszenie sferyczności.
  • Skutki: MNK nadal nieobciążony i zgodny, lecz nie-BLUE; klasyczne SE niezgodne (prawdziwa wariancja jest kanapkowa).
  • Diagnoza: wykres reszt (lejek), wykres skala–położenie, testy Breuscha-Pagana ($\mathrm{LM}=nR^2_{\text{aux}}$) i White’a.
  • Leczenie: odporne błędy HC0/HC3 (bez założeń o strukturze) albo WLS/GLS (efektywne, gdy struktura znana).
  • Granica: korekta dotyczy wnioskowania, nie usuwa obciążenia ze złej specyfikacji.

Dalej: Autokorelacja · Założenia KMNK · Testy hipotez · Szeregi czasowe

Literatura uzupełniająca
  • Podręcznik: J. M. Wooldridge, Introductory Econometrics: A Modern Approach, rozdz. 8
  • Podręcznik: W. H. Greene, Econometric Analysis, rozdz. 9 (heteroskedastyczność, GLS)
  • H. White, „A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity", Econometrica 48(4), 1980, s. 817–838
  • T. S. Breusch, A. R. Pagan, „A Simple Test for Heteroscedasticity and Random Coefficient Variation", Econometrica 47(5), 1979, s. 1287–1294
  • Pakiety R: sandwich, lmtest
Oprogramowanie

R:

library(sandwich); library(lmtest)

## Zbiór z lejkiem (przyklad z tekstu)
x <- 1:8; y <- c(3, 5, 6, 11, 9, 16, 12, 22)
m <- lm(y ~ x)

## Test Breuscha-Pagana: LM = n * R^2 z regresji e^2 ~ x
bptest(m)                                   # LM = 5,31; df = 1; p = 0,021

## Test White'a (z kwadratami i interakcjami)
bptest(m, varformula = ~ x + I(x^2))

## Odporne bledy standardowe
coeftest(m, vcov = vcovHC(m, type = "HC0"))  # White
coeftest(m, vcov = vcovHC(m, type = "HC3"))  # zalecane w malych probach

## WLS przy zalozeniu sigma_i^2 ~ x_i (wagi 1/x)
m_wls <- lm(y ~ x, weights = 1 / x)          # beta = (0,547; 2,212)

Python:

import numpy as np, statsmodels.api as sm
import statsmodels.stats.diagnostic as diag

x = np.arange(1, 9); y = np.array([3, 5, 6, 11, 9, 16, 12, 22])
X = sm.add_constant(x)
m = sm.OLS(y, X).fit()

# Breusch-Pagan
bp = diag.het_breuschpagan(m.resid, X)
print(f"LM = {bp[0]:.3f}, p = {bp[1]:.4f}")   # LM = 5,31; p = 0,021

# Odporne bledy standardowe
print(sm.OLS(y, X).fit(cov_type='HC0').bse)
print(sm.OLS(y, X).fit(cov_type='HC3').bse)   # wieksze niz klasyczne

# WLS (wagi 1/x)
print(sm.WLS(y, X, weights=1.0 / x).fit().params)