Heteroskedastyczność — diagnoza i leczenie

Definicja heteroskedastyczności

W klasycznym modelu regresji zakładamy, że reszty mają stałą wariancję (homoskedastyczność):

$$Var(\varepsilon_i \mid X_i) = \sigma^2 \quad \text{dla wszystkich } i$$

Heteroskedastyczność to naruszenie tego założenia — wariancja błędów zmienia się wraz ze zmiennymi objaśniającymi:

$$Var(\varepsilon_i \mid X_i) = \sigma_i^2 \neq \text{stała}$$

Intuicja

Wyobraź sobie regresję dochodu na wykształcenie:

Osoby z małym wykształceniem zarabiają podobnie (mała wariancja)
Osoby z wysokim wykształceniem — od lekarza do artysty — mają ogromne rozpiętości dochodów (duża wariancja)

To klasyczny przykład heteroskedastyczności rosnącej.

Najczęstsze źródła

Dane przekrojowe z dużym zróżnicowaniem (firmy różnej wielkości, kraje)
Zmienne ilorazowe (PKB, dochód, cena)
Modele liniowe dla zmiennych procentowych
Dane panelowe — heteroskedastyczność między jednostkami

Konsekwencje dla regresji MNK

Dobra wiadomość: Estymatory MNK są nadal nieobciążone i zgodne.

Zła wiadomość: Estymatory MNK nie są efektywne (nie są BLUE) — błędy standardowe są błędnie oszacowane.

Konkretnie:

Standardowe $SE(\hat{\beta})$ są za małe lub za duże
Statystyki $t$ i $F$ mają błędny rozkład
Przedziały ufności i testy hipotez są nieprawidłowe

Główne zagrożenie: Możemy stwierdzić istotność tam, gdzie jej nie ma (lub odwrotnie).

Wizualna diagnoza

Pierwsza i najważniejsza: wykres reszt vs. dopasowane wartości.

plot(model, which = 1)  # reszty vs. dopasowane
plot(model, which = 3)  # skala-lokacja (sqrt(|reszty|) vs. dopasowane)

Homoskedastyczność: reszty rozłożone symetrycznie, bez wzorca, wokół zera.

Heteroskedastyczność: rozrzut reszt rośnie (lub maleje) wraz z $\hat{y}$ — charakterystyczny “lejek”.

Testy statystyczne

Test Breusha-Pagana

Regresja kwadratów reszt na zmienne objaśniające:

$$\hat{\varepsilon}_i^2 = \delta_0 + \delta_1 x_{i1} + \cdots + \delta_k x_{ik} + u_i$$

Statystyka: $LM = n \cdot R^2_{\hat{\varepsilon}^2} \sim \chi^2_k$

$H_0$: homoskedastyczność (wszystkie $\delta_j = 0$)

Test White’a

Rozbudowana wersja BP — testuje też nieliniowe związki:

$$\hat{\varepsilon}_i^2 = \delta_0 + \sum_j \delta_j x_{ij} + \sum_j \delta_{jj} x_{ij}^2 + \sum_{jBardziej ogólny, ale traci moc przy małych próbach.

Interpretacja

$p < 0.05$ → odrzucamy homoskedastyczność → heteroskedastyczność jest problem
Ale testy mają słabą moc przy małych $n$ — patrz przede wszystkim na wykres!

Rozwiązania

1. Odporne błędy standardowe (HC Standard Errors)

Najpopularniejsze rozwiązanie — nie zmienia estymatorów, tylko koryguje ich błędy standardowe.

$$\widehat{Var}_{HC3}(\hat{\beta}) = (\mathbf{X}^T\mathbf{X})^{-1}\left[\sum_i \frac{\hat{\varepsilon}_i^2}{(1-h_{ii})^2} x_i x_i^T\right](\mathbf{X}^T\mathbf{X})^{-1}$$

gdzie $h_{ii}$ — elementy diagonalne macierzy hat ($H = X(X^TX)^{-1}X^T$).

Typy:

HC0 — oryginalne White’a
HC1 — korygowane o stopnie swobody
HC3 — zalecane w małych próbach (koryguje o leverage)

# R
coeftest(model, vcov = vcovHC(model, type = "HC3"))

# Python
model.fit(cov_type='HC3')

Zaleta: Zawsze bezpieczne stosować — nawet jeśli nie ma heteroskedastyczności, wyniki będą podobne do klasycznych SE.

2. WLS — Ważona Metoda Najmniejszych Kwadratów

Gdy znamy (lub dobrze estymujemy) strukturę wariancji $\sigma_i^2 \propto h(x_i)$:

Minimalizujemy ważoną sumę kwadratów reszt:

$$\sum_{i=1}^n \frac{1}{\sigma_i^2} (y_i - \mathbf{x}_i^T \boldsymbol{\beta})^2$$

Odpowiednik: mnożymy każdą obserwację przez $1/\sigma_i$.

Problem: Rzadko znamy dokładnie $\sigma_i^2$. Musimy je estymować — to FGLS (Feasible GLS).

# Estymacja wag z pomocniczej regresji
model_aux <- lm(log(residuals(model)^2) ~ log(x1) + log(x2))
sigma_hat <- exp(fitted(model_aux))
model_wls <- lm(y ~ x1 + x2, weights = 1/sigma_hat^2)

3. Transformacja logarytmiczna

Często wystarczy zlogarytmować zmienną zależną:

$$\ln(y_i) = \beta_0 + \beta_1 x_{i1} + \cdots + \varepsilon_i$$

Logarytm “ściska” wartości — rozkłady prawoskośne stają się bardziej symetryczne i homoskedastyczne.

Stosuj gdy: $y > 0$, dane są prawoskośne, wahania procentowe są ważniejsze od bezwzględnych.

Heteroskedastyczność w danych panelowych

W danych panelowych ($N$ jednostek, $T$ okresów) heteroskedastyczność ma specyficzną strukturę:

Między jednostkami: firmy różnej wielkości mają różne wariancje błędów
Wewnątrz jednostki: wariancja zmienia się w czasie

library(plm)
model_panel <- plm(wage ~ educ + exper, data = panel_data,
                   index = c("id", "year"), model = "within")
                   
# Odporne błędy na heteroskedastyczność i autokorelację
coeftest(model_panel, vcov = vcovHC(model_panel, type = "HC3",
                                     cluster = "group"))

Clustered SE — odporne błędy z klasteringiem po jednostkach: korygują i za heteroskedastyczność, i za autokorelację wewnątrz jednostki.

Podsumowanie: kiedy co stosować

Sytuacja	Zalecenie
Nie wiesz czy jest problem	Zawsze używaj HC3
Potwierdzono heteroskedastyczność	HC3 lub WLS
Znasz strukturę $\sigma_i^2$	WLS (efektywniejszy)
Dane prawoskośne ($y > 0$)	Transformacja log
Dane panelowe	Clustered SE

Złota zasada: Jeśli próba jest duża, HC3 chroni cię bez kosztów. Stosuj domyślnie.

Następnie: Szeregi czasowe — analiza i prognozowanie