Regresja liniowa
Kompletny rozdział podręcznikowy: od intuicji ekonomicznej przez równania normalne, postać macierzową, twierdzenie Gaussa-Markowa z dowodem, wnioskowanie statystyczne, diagnostykę modelu, predykcję a przyczynowość, po trzy studia przypadków z kodem R.
Wymagana znajomość: algebra liniowa (podstawy), statystyka opisowa, rachunek różniczkowy jednej zmiennej.
Rozdział prowadzi od najprostszego modelu z jednym regresorem do postaci macierzowej, twierdzenia Gaussa-Markowa, wnioskowania statystycznego, diagnostyki i trzech studiów przypadków z kodem R. Regresja liniowa jest jednocześnie metodą opisu danych, procedurą estymacji parametrów modelu statystycznego i geometrycznym rzutem wektora na podprzestrzeń kolumnową — każde z tych ujęć pojawia się w kolejnych rozdziałach.
Motywacja i pierwszy model
Losowość danych empirycznych
Rozważamy przekrojowe dane o gospodarstwach domowych: dla każdego obserwujemy miesięczny dochód i miesięczną konsumpcję. Zależność między tymi wielkościami jest dodatnia, lecz niedeterministyczna — gospodarstwa o tym samym dochodzie różnią się konsumpcją wskutek odmiennej struktury wydatków, stopy oszczędzania, liczby domowników i wielu innych czynników nieobserwowanych.
Dane empiryczne nie tworzą funkcji — tworzą chmurę punktów. Regresja liniowa wyznacza prostą opisującą przeciętny kierunek tej zależności.
Rozważmy mały zbiór danych. Dochód i konsumpcja są wyrażone w tysiącach złotych miesięcznie.
Tabela 1.1. Przykładowe dane o dochodzie i konsumpcji.
| Gospodarstwo | Dochód $x_i$ | Konsumpcja $y_i$ | Uwagi ekonomiczne |
|---|---|---|---|
| 1 | 3,0 | 2,7 | niski dochód, małe oszczędności |
| 2 | 3,5 | 3,0 | wydatki podstawowe dominują |
| 3 | 4,0 | 3,2 | część dochodu oszczędzana |
| 4 | 4,5 | 3,8 | wyższe wydatki bieżące |
| 5 | 5,0 | 4,0 | przeciętne gospodarstwo |
| 6 | 5,5 | 4,2 | stabilna konsumpcja |
| 7 | 6,0 | 4,7 | wzrost wydatków |
| 8 | 6,5 | 5,0 | wyższy standard życia |
| 9 | 7,0 | 5,1 | większa stopa oszczędzania |
| 10 | 7,5 | 5,6 | konsumpcja rośnie mniej niż dochód |
Zależność jest dodatnia, ale nie idealna. Regresja liniowa nie mówi, że każde gospodarstwo o dochodzie 6 tysięcy złotych wyda dokładnie pewną kwotę. Mówi, jaka jest przeciętna wartość konsumpcji przy danym dochodzie, zgodnie z przyjętym liniowym opisem.
Pierwszy model
Najprostszy model regresji liniowej zapisujemy jako:
$$Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i, \qquad i = 1, \ldots, n.$$W przykładzie konsumpcji:
$$\mathit{Consumption}_i = \beta_0 + \beta_1 \mathit{Income}_i + \varepsilon_i.$$Każdy element tego zapisu ma znaczenie.
- $Y_i$ to zmienna objaśniana, czyli wielkość, którą chcemy opisać lub przewidzieć.
- $X_i$ to zmienna objaśniająca, czyli wielkość, za pomocą której opisujemy $Y_i$.
- $\beta_0$ to wyraz wolny.
- $\beta_1$ to nachylenie prostej regresji.
- $\varepsilon_i$ to składnik losowy, czyli wszystko to, czego nie uwzględniliśmy w modelu.
Składnik losowy nie oznacza błędu technicznego badacza — skupia on pominięte czynniki, błędy pomiaru i losowe zaburzenia, których model nie wyjaśnia.
Interpretacja parametru nachylenia
Parametr $\beta_1$ mierzy zmianę wartości oczekiwanej $Y$ związaną ze wzrostem $X$ o jedną jednostkę. W modelu konsumpcji mówi, o ile tysięcy złotych wzrasta przeciętna konsumpcja, gdy dochód wzrasta o jeden tysiąc złotych.
Jeżeli:
$$\widehat{\mathit{Consumption}} = 1{,}1 + 0{,}62 \cdot \mathit{Income},$$to interpretacja współczynnika $0{,}62$ brzmi:
Wzrost dochodu o 1 tysiąc złotych wiąże się przeciętnie ze wzrostem konsumpcji o 0,62 tysiąca złotych, czyli o 620 złotych.
Słowo „wiąże się" jest tu celowe. Sama regresja nie dowodzi jeszcze przyczynowości. Do wniosku przyczynowego potrzebujemy dodatkowej argumentacji: jakie zmienne zostały pominięte, czy dochód jest egzogeniczny, czy nie występuje odwrócona zależność, czy pomiar jest poprawny.
Ilustracja graficzna
Estymator metody najmniejszych kwadratów
Kryterium najmniejszych kwadratów
Mamy dane $(x_1, y_1), \ldots, (x_n, y_n)$ i chcemy wybrać prostą:
$$\hat{y}_i = b_0 + b_1 x_i.$$Dla każdej obserwacji powstaje reszta:
$$e_i = y_i - \hat{y}_i = y_i - b_0 - b_1 x_i.$$Regresja metodą najmniejszych kwadratów wybiera takie $b_0$ i $b_1$, aby minimalizować sumę kwadratów reszt:
$$S(b_0, b_1) = \sum_{i=1}^{n}(y_i - b_0 - b_1 x_i)^2.$$Kwadraty mają trzy zalety. Po pierwsze, dodatnie i ujemne odchylenia nie znoszą się. Po drugie, duże błędy są karane silniej niż małe. Po trzecie, funkcja celu jest gładka i wygodna matematycznie.
Równania normalne
Różniczkujemy funkcję $S$ względem $b_0$:
$$\frac{\partial S}{\partial b_0} = \sum_{i=1}^{n} 2(y_i - b_0 - b_1 x_i)(-1).$$Warunek pierwszego rzędu wymaga:
$$\sum_{i=1}^{n}(y_i - b_0 - b_1 x_i) = 0.$$Po rozdzieleniu sum:
$$\sum_{i=1}^{n} y_i - nb_0 - b_1 \sum_{i=1}^{n} x_i = 0.$$Stąd:
$$\sum_{i=1}^{n} y_i = nb_0 + b_1 \sum_{i=1}^{n} x_i.$$Dzielimy przez $n$:
$$\bar{y} = b_0 + b_1\bar{x}.$$Otrzymujemy:
$$b_0 = \bar{y} - b_1\bar{x}.$$To oznacza, że prosta regresji przechodzi przez punkt średnich $(\bar{x}, \bar{y})$.
Teraz różniczkujemy względem $b_1$:
$$\frac{\partial S}{\partial b_1} = \sum_{i=1}^{n} 2(y_i - b_0 - b_1 x_i)(-x_i).$$Warunek pierwszego rzędu:
$$\sum_{i=1}^{n} x_i(y_i - b_0 - b_1 x_i) = 0.$$Po rozdzieleniu:
$$\sum_{i=1}^{n} x_i y_i = b_0 \sum_{i=1}^{n} x_i + b_1 \sum_{i=1}^{n} x_i^2.$$Podstawiając $b_0 = \bar{y} - b_1\bar{x}$ i porządkując wyrazy, otrzymujemy:
$$\hat{\beta}_1 = b_1 = \frac{\displaystyle\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\displaystyle\sum_{i=1}^{n}(x_i - \bar{x})^2}.$$A wyraz wolny:
$$\hat{\beta}_0 = b_0 = \bar{y} - \hat{\beta}_1\bar{x}.$$Nachylenie regresji jest kowariancją $x$ i $y$ podzieloną przez wariancję $x$: gdy obie zmienne odchylają się od swoich średnich w tym samym kierunku, licznik jest dodatni i prosta rośnie; gdy w przeciwnych — licznik jest ujemny.
Przykład
Rozważmy pięć obserwacji:
$$x = (1, 2, 3, 4, 5), \qquad y = (2, 4, 5, 4, 5).$$Obliczamy średnie: $\bar{x} = 3$, $\bar{y} = 4$.
Odchylenia: $x_i - \bar{x} = (-2, -1, 0, 1, 2)$.
Licznik:
$$\sum(x_i - \bar{x})(y_i - \bar{y}) = (-2)(-2) + (-1)(0) + (0)(1) + (1)(0) + (2)(1) = 6.$$Mianownik:
$$\sum(x_i - \bar{x})^2 = 4 + 1 + 0 + 1 + 4 = 10.$$Nachylenie:
$$\hat{\beta}_1 = \frac{6}{10} = 0{,}6.$$Wyraz wolny:
$$\hat{\beta}_0 = 4 - 0{,}6 \cdot 3 = 2{,}2.$$Linia regresji: $\hat{y}_i = 2{,}2 + 0{,}6 \cdot x_i$.
Własności algebraiczne estymatora MNK
Reszty $e_i = y_i - \hat{y}_i$ mają kilka własności wynikających wyłącznie z algebry minimalizacji — bez żadnych założeń statystycznych.
- Suma reszt wynosi zero: $\sum_{i=1}^{n} e_i = 0$.
- Reszty są nieskorelowane ze zmienną $X$: $\sum_{i=1}^{n} x_i e_i = 0$.
- Prosta przechodzi przez punkt $(\bar{x}, \bar{y})$.
Należy odróżnić te własności od założenia egzogeniczności. Reszty $e_i$ z estymowanej regresji są z definicji nieskorelowane z $X$ w próbie — wynika to wyłącznie z algebry minimalizacji. Założenie egzogeniczności dotyczy nieobserwowalnego składnika losowego $\varepsilon_i$, a nie obliczonych reszt.
Model statystyczny i założenia klasyczne
Model probabilistyczny a dopasowanie algebraiczne
Estymator MNK wyznacza prostą dla każdego zbioru danych, niezależnie od tego, jak dane zostały wygenerowane. Aby mówić o własnościach estymatora — nieobciążoności, wariancji, przedziałach ufności — potrzebujemy modelu statystycznego: zbioru założeń o procesie generującym dane.
Założenia klasycznego modelu regresji liniowej
Model regresji liniowej w postaci macierzowej $\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ nazywamy klasycznym, jeżeli spełnione są następujące warunki:
- model jest liniowy względem parametrów;
- macierz $\mathbf{X}$ ma pełny rząd kolumnowy;
- $E(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \mathbf{0}$;
- $\mathrm{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \sigma^2 \mathbf{I}$;
- do dokładnego wnioskowania w małych próbach często dodaje się normalność: $\boldsymbol{\varepsilon} \mid \mathbf{X} \sim N(\mathbf{0}, \sigma^2 \mathbf{I})$.
Liniowość względem parametrów
Model $Y_i = \beta_0 + \beta_1 \ln(X_i) + \varepsilon_i$ jest liniowy względem parametrów, bo $\ln(X_i)$ można traktować jako nową zmienną objaśniającą. Model $Y_i = \frac{1}{\beta_0 + \beta_1 X_i} + \varepsilon_i$ nie jest liniowy względem parametrów.
Pełny rząd macierzy $X$
Pełny rząd oznacza brak doskonałej współliniowości. Żadna kolumna macierzy $X$ nie może być dokładną kombinacją liniową pozostałych kolumn. Jeżeli w modelu umieścimy stałą, zmienną $\mathit{Male}$ i zmienną $\mathit{Female}$, gdzie $\mathit{Male} + \mathit{Female} = 1$ dla każdej osoby, to mamy doskonałą współliniowość. Jedna z tych zmiennych musi zostać pominięta.
Egzogeniczność
Warunek $E(\varepsilon_i \mid X_i) = 0$ mówi, że składnik losowy nie jest systematycznie powiązany z regresorem. Jeżeli badamy wpływ edukacji na płace, ale pomijamy zdolności. Jeżeli zdolności wpływają na edukację i płace, to $E(\varepsilon_i \mid \mathit{Education}_i)$ nie jest równe zero.
Homoskedastyczność i brak autokorelacji
Warunek $\mathrm{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \sigma^2 \mathbf{I}$ zawiera dwa elementy: jednakową wariancję $\sigma^2$ dla każdej obserwacji (homoskedastyczność) i brak korelacji między składnikami losowymi (brak autokorelacji). Naruszenie pierwszego warunku nazywamy heteroskedastycznością, naruszenie drugiego — autokorelacją.
Postać macierzowa regresji
Zapis modelu wielorakiego
W regresji wielorakiej mamy więcej niż jedną zmienną objaśniającą:
$$Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_k X_{ik} + \varepsilon_i.$$W zapisie macierzowym:
$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},$$gdzie:
$$\mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, \qquad \mathbf{X} = \begin{bmatrix} 1 & x_{11} & \cdots & x_{1k} \\ 1 & x_{21} & \cdots & x_{2k} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & \cdots & x_{nk} \end{bmatrix}, \qquad \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_k \end{bmatrix}.$$Estymator MNK w postaci macierzowej
Minimalizujemy:
$$S(\mathbf{b}) = (\mathbf{y} - \mathbf{X}\mathbf{b})'(\mathbf{y} - \mathbf{X}\mathbf{b}).$$Rozwijamy:
$$S(\mathbf{b}) = \mathbf{y}'\mathbf{y} - 2\mathbf{b}'\mathbf{X}'\mathbf{y} + \mathbf{b}'\mathbf{X}'\mathbf{X}\mathbf{b}.$$Gradient funkcji $S$ względem $\mathbf{b}$ przyjmuje postać:
$$\frac{\partial S}{\partial \mathbf{b}} = -2\mathbf{X}'\mathbf{y} + 2\mathbf{X}'\mathbf{X}\mathbf{b}.$$Warunek pierwszego rzędu:
$$\mathbf{X}'\mathbf{X}\hat{\boldsymbol{\beta}} = \mathbf{X}'\mathbf{y}.$$Jeżeli $\mathbf{X}'\mathbf{X}$ jest odwracalna:
$$\hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}.$$Macierz projekcji
Wartości dopasowane:
$$\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}.$$Macierz:
$$\mathbf{P} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$$nazywa się macierzą projekcji albo macierzą kapeluszową. Nazwa pochodzi stąd, że przekształca $\mathbf{y}$ w $\hat{\mathbf{y}}$:
$$\hat{\mathbf{y}} = \mathbf{P}\mathbf{y}.$$Macierz reszt:
$$\mathbf{M} = \mathbf{I} - \mathbf{P}.$$Reszty:
$$\hat{\boldsymbol{\varepsilon}} = \mathbf{y} - \hat{\mathbf{y}} = \mathbf{M}\mathbf{y}.$$Interpretacja geometryczna
Kolumny macierzy $X$ rozpinają podprzestrzeń w $\mathbb{R}^n$. Wektor obserwacji $\mathbf{y}$ zwykle nie leży dokładnie w tej podprzestrzeni. MNK znajduje punkt w tej podprzestrzeni, który jest najbliżej $\mathbf{y}$ w sensie odległości euklidesowej. Tym punktem jest $\hat{\mathbf{y}}$.
Własności estymatora MNK
Rozkład estymatora wokół prawdziwego parametru
Podstawiamy model $\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ do estymatora:
$$\hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'(\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}).$$Po rozwinięciu:
$$\hat{\boldsymbol{\beta}} = \boldsymbol{\beta} + (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\boldsymbol{\varepsilon}.$$Wzór ten ma kluczowe znaczenie. Pokazuje, że błąd estymacji wynosi:
$$\hat{\boldsymbol{\beta}} - \boldsymbol{\beta} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\boldsymbol{\varepsilon}.$$Jeżeli składnik losowy nie jest systematycznie powiązany z $X$, przeciętnie błąd estymacji znika.
Nieobciążoność
Warunkowo względem $X$:
$$E(\hat{\boldsymbol{\beta}} \mid \mathbf{X}) = \boldsymbol{\beta} + (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'E(\boldsymbol{\varepsilon} \mid \mathbf{X}).$$Jeżeli:
$$E(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \mathbf{0},$$to:
$$E(\hat{\boldsymbol{\beta}} \mid \mathbf{X}) = \boldsymbol{\beta}.$$Estymator MNK jest wtedy nieobciążony.
Wariancja estymatora
Zakładamy:
$$\mathrm{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \sigma^2 \mathbf{I}.$$Wtedy:
$$\mathrm{Var}(\hat{\boldsymbol{\beta}} \mid \mathbf{X}) = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}.$$W regresji prostej wariancja nachylenia wynosi:
$$\mathrm{Var}(\hat{\beta}_1 \mid X) = \frac{\sigma^2}{\displaystyle\sum_{i=1}^{n}(x_i - \bar{x})^2}.$$Wynik ten ma prostą interpretację. Estymator jest dokładniejszy, gdy wariancja błędu $\sigma^2$ jest mniejsza oraz gdy zmienna $X$ ma większą zmienność w próbie.
Zgodność
Nieobciążoność dotyczy średniej po wielu hipotetycznych próbach o tej samej liczebności. Zgodność dotyczy zachowania estymatora, gdy liczba obserwacji rośnie. W uproszczeniu:
$$\hat{\boldsymbol{\beta}} \xrightarrow{p} \boldsymbol{\beta}.$$Aby estymator MNK był zgodny, potrzebujemy między innymi, aby:
$$\frac{1}{n}\mathbf{X}'\boldsymbol{\varepsilon} \xrightarrow{p} 0$$oraz aby:
$$\frac{1}{n}\mathbf{X}'\mathbf{X}$$miała granicę odwracalną.
Twierdzenie Gaussa-Markowa
Dowód. Niech inny liniowy estymator ma postać:
$$\tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{y}.$$Estymator MNK można zapisać:
$$\hat{\boldsymbol{\beta}} = \mathbf{A}\mathbf{y}, \qquad \mathbf{A} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'.$$Zapiszmy:
$$\mathbf{C} = \mathbf{A} + \mathbf{D}.$$Nieobciążoność $\tilde{\boldsymbol{\beta}}$ wymaga:
$$E(\tilde{\boldsymbol{\beta}} \mid X) = \mathbf{C}\mathbf{X}\boldsymbol{\beta} = \boldsymbol{\beta}$$dla każdego $\boldsymbol{\beta}$. Stąd:
$$\mathbf{C}\mathbf{X} = \mathbf{I}.$$Ponieważ $\mathbf{A}\mathbf{X} = \mathbf{I}$, otrzymujemy:
$$\mathbf{D}\mathbf{X} = \mathbf{0}.$$Wariancja alternatywnego estymatora:
$$\mathrm{Var}(\tilde{\boldsymbol{\beta}} \mid X) = \sigma^2 \mathbf{C}\mathbf{C}'.$$Podstawiamy $\mathbf{C} = \mathbf{A} + \mathbf{D}$:
$$\mathbf{C}\mathbf{C}' = \mathbf{A}\mathbf{A}' + \mathbf{A}\mathbf{D}' + \mathbf{D}\mathbf{A}' + \mathbf{D}\mathbf{D}'.$$Ponieważ $\mathbf{D}\mathbf{X} = \mathbf{0}$, mamy $\mathbf{D}\mathbf{A}' = \mathbf{D}\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1} = \mathbf{0}$. Analogicznie $\mathbf{A}\mathbf{D}' = \mathbf{0}$. Zatem:
$$\mathrm{Var}(\tilde{\boldsymbol{\beta}} \mid X) = \mathrm{Var}(\hat{\boldsymbol{\beta}} \mid X) + \sigma^2 \mathbf{D}\mathbf{D}'.$$Macierz $\mathbf{D}\mathbf{D}'$ jest nieujemnie określona. Dlatego wariancja dowolnego innego liniowego nieobciążonego estymatora jest nie mniejsza niż wariancja MNK. $\square$
Twierdzenie Gaussa-Markowa dotyczy wyłącznie klasy estymatorów liniowych i nieobciążonych przy klasycznych założeniach. Przy kryterium minimalizacji błędu średniokwadratowego procedury obciążone — takie jak ridge regression — mogą w pewnych konfiguracjach uzyskiwać niższy MSE.
Dopasowanie modelu i współczynnik determinacji
Rozkład zmienności
W modelu ze stałą całkowitą zmienność $Y$ mierzymy przez:
$$TSS = \sum_{i=1}^{n}(y_i - \bar{y})^2.$$Zmienność wyjaśniona:
$$ESS = \sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2.$$Zmienność niewyjaśniona:
$$RSS = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2 = \sum_{i=1}^{n} e_i^2.$$Dla MNK ze stałą zachodzi:
$$TSS = ESS + RSS.$$
Współczynnik $R^2$
Współczynnik determinacji definiujemy jako:
$$R^2 = \frac{ESS}{TSS} = 1 - \frac{RSS}{TSS}.$$$R^2$ mówi, jaka część zmienności zmiennej objaśnianej została odtworzona przez liniowe dopasowanie modelu.
Wysokie $R^2$ nie jest warunkiem wystarczającym do uznania modelu za poprawny przyczynowo — regresje na zmiennych trendowych dają wysokie $R^2$ bez żadnej zależności ekonomicznej. Niskie $R^2$ nie wyklucza z kolei istotności szacowanego parametru: w badaniach mikroekonomicznych wiele wartościowych modeli wyjaśnia poniżej 30% wariancji zmiennej zależnej.
Skorygowane $R^2$
Dodanie kolejnej zmiennej do modelu nigdy nie zwiększa $RSS$, dlatego zwykłe $R^2$ nie maleje po dodaniu regresora. Skorygowane $R^2$ karze za liczbę parametrów:
$$\bar{R}^2 = 1 - \frac{RSS/(n-k-1)}{TSS/(n-1)}.$$Jest użyteczne do porównywania modeli o różnej liczbie zmiennych, ale nie powinno zastępować rozumowania ekonomicznego.
Wnioskowanie statystyczne
Estymacja wariancji składnika losowego
Ponieważ $\sigma^2$ nie jest znane, estymujemy je za pomocą reszt:
$$\hat{\sigma}^2 = \frac{RSS}{n - k - 1}.$$Dzielimy przez $n-k-1$, a nie przez $n$, ponieważ estymujemy $k+1$ parametrów. Liczba $n - k - 1$ to liczba stopni swobody reszt.
Błędy standardowe
Estymowana macierz wariancji MNK:
$$\widehat{\mathrm{Var}}(\hat{\boldsymbol{\beta}}) = \hat{\sigma}^2(\mathbf{X}'\mathbf{X})^{-1}.$$Błąd standardowy $\hat{\beta}_j$ to pierwiastek z $j$-tego elementu diagonalnego tej macierzy:
$$se(\hat{\beta}_j) = \sqrt{\widehat{\mathrm{Var}}(\hat{\beta}_j)}.$$Test $t$
Aby sprawdzić hipotezę:
$$H_0: \beta_j = \beta_{j,0},$$liczymy statystykę:
$$t = \frac{\hat{\beta}_j - \beta_{j,0}}{se(\hat{\beta}_j)}.$$W klasycznym modelu z normalnością błędów statystyka ma rozkład $t$-Studenta z $n - k - 1$ stopniami swobody.
Przedział ufności
Przedział ufności dla $\beta_j$ na poziomie $1 - \alpha$:
$$\hat{\beta}_j \pm t_{1-\alpha/2,\, n-k-1} \cdot se(\hat{\beta}_j).$$Interpretacja nie brzmi: „prawdopodobieństwo, że parametr leży w tym konkretnym przedziale, wynosi 95%". Parametr jest stały. Poprawna interpretacja częstościowa brzmi: gdybyśmy wielokrotnie losowali próbę i budowali przedziały tą samą metodą, około 95% takich przedziałów pokrywałoby prawdziwy parametr.
Test $F$
Test $F$ służy do testowania wielu ograniczeń jednocześnie, na przykład:
$$H_0: \beta_2 = \beta_3 = \beta_4 = 0.$$Porównujemy model nieograniczony z modelem ograniczonym. Statystyka:
$$F = \frac{(RSS_R - RSS_U)/q}{RSS_U/(n-k-1)},$$gdzie $q$ to liczba ograniczeń, $RSS_R$ to suma kwadratów reszt w modelu ograniczonym, a $RSS_U$ w modelu nieograniczonym.
Regresja wieloraka i interpretacja ceteris paribus
Znaczenie dodatkowych zmiennych objaśniających
Regresja prosta często daje mylące wyniki, ponieważ pomija inne ważne czynniki. Przykład: płaca zależy nie tylko od edukacji, ale także od doświadczenia, branży, miasta, płci, umiejętności, znajomości języków, stanowiska i wielu innych cech.
Model wieloraki:
$$\mathit{Wage}_i = \beta_0 + \beta_1 \mathit{Education}_i + \beta_2 \mathit{Experience}_i + \beta_3 \mathit{Female}_i + \varepsilon_i.$$Współczynnik $\beta_1$ interpretujemy jako różnicę w przeciętnej płacy związaną z jednym dodatkowym rokiem edukacji przy ustalonym doświadczeniu i płci.
Znaczenie „przy pozostałych zmiennych stałych"
Warunek ceteris paribus nie oznacza, że w rzeczywistości obserwujemy dwie identyczne osoby różniące się tylko edukacją. Oznacza, że model statystycznie porównuje osoby o takim samym poziomie zmiennych kontrolnych uwzględnionych w modelu. Jeśli ważna cecha nie została uwzględniona, interpretacja ceteris paribus może być fałszywa.
Twierdzenie Frisch-Waugh-Lovella
Regresja wieloraka mierzy wpływ tej części $X_1$, która nie jest liniowo wyjaśniona przez pozostałe zmienne kontrolne — stąd mowa o efekcie „po oczyszczeniu" z innych regresorów.
Zmienne jakościowe, interakcje i postaci funkcyjne
Zmienne zero-jedynkowe
Zmienna zero-jedynkowa przyjmuje wartości 0 albo 1. Przykład:
$$\mathit{Female}_i = \begin{cases} 1, & \text{jeżeli osoba jest kobietą,} \\ 0, & \text{w przeciwnym razie.} \end{cases}$$Model:
$$\mathit{Wage}_i = \beta_0 + \beta_1 \mathit{Education}_i + \beta_2 \mathit{Female}_i + \varepsilon_i.$$Współczynnik $\beta_2$ oznacza różnicę w przeciętnej płacy między grupą $\mathit{Female} = 1$ a grupą bazową $\mathit{Female} = 0$, przy tym samym poziomie edukacji.
Pułapka zmiennych fikcyjnych
Jeżeli mamy trzy regiony: A, B i C, nie możemy w modelu ze stałą umieścić trzech zmiennych zero-jedynkowych dla wszystkich regionów. Musimy pominąć jedną kategorię bazową. Wtedy współczynniki przy pozostałych kategoriach interpretujemy względem tej bazy.
Interakcje
Interakcja pozwala, aby wpływ jednej zmiennej zależał od poziomu drugiej. Model:
$$\mathit{Wage}_i = \beta_0 + \beta_1 \mathit{Education}_i + \beta_2 \mathit{Female}_i + \beta_3(\mathit{Education}_i \cdot \mathit{Female}_i) + \varepsilon_i.$$Dla grupy bazowej $\mathit{Female} = 0$ wpływ edukacji wynosi $\beta_1$. Dla $\mathit{Female} = 1$ wpływ edukacji wynosi $\beta_1 + \beta_3$.
Logarytmy
Modele logarytmiczne są w ekonomii bardzo częste. Jeśli:
$$\ln(Y_i) = \beta_0 + \beta_1 X_i + \varepsilon_i,$$to wzrost $X$ o jednostkę wiąże się w przybliżeniu ze zmianą $Y$ o $100\beta_1$ procent. Jeśli:
$$\ln(Y_i) = \beta_0 + \beta_1 \ln(X_i) + \varepsilon_i,$$to $\beta_1$ jest elastycznością: wzrost $X$ o 1% wiąże się ze zmianą $Y$ o $\beta_1$ procent.
Składniki kwadratowe
Model:
$$Y_i = \beta_0 + \beta_1 X_i + \beta_2 X_i^2 + \varepsilon_i$$pozwala na zależność nieliniową względem $X$, ale nadal jest liniowy względem parametrów. Efekt krańcowy $X$ wynosi:
$$\frac{\partial E(Y \mid X)}{\partial X} = \beta_1 + 2\beta_2 X.$$
Diagnostyka modelu
Rola diagnostyki modelu
Oszacowanie regresji jest początkiem analizy, a nie jej końcem. Model może mieć istotne współczynniki i wysokie $R^2$, a jednocześnie być źle wyspecyfikowany. Diagnostyka służy do sprawdzenia, czy wyniki są wiarygodne.
Wykres reszt względem wartości dopasowanych
Na osi poziomej umieszczamy $\hat{y}_i$, a na osi pionowej reszty $e_i$. Dobry model nie powinien pokazywać systematycznego wzorca. Jeżeli widzimy krzywiznę, model liniowy może być zbyt prosty. Jeżeli rozrzut reszt rośnie wraz z dopasowaniem, mamy podejrzenie heteroskedastyczności.
Heteroskedastyczność
Heteroskedastyczność oznacza, że wariancja błędu zależy od poziomu zmiennych objaśniających:
$$\mathrm{Var}(\varepsilon_i \mid X_i) \neq \sigma^2.$$Klasyczne błędy standardowe mogą być wtedy błędne. Rozwiązaniem często są odporne błędy standardowe White’a albo estymacja ważona.
Autokorelacja
W danych czasowych błędy mogą być skorelowane w czasie:
$$\mathrm{Cov}(\varepsilon_t, \varepsilon_{t-1}) \neq 0.$$Wtedy klasyczne błędy standardowe nie są wiarygodne. Stosuje się między innymi błędy standardowe Neweya-Westa albo modele dynamiczne.
Współliniowość
Wysoka współliniowość oznacza silną zależność między regresorami. Nie powoduje obciążenia MNK, ale zwiększa wariancję estymatorów. Współczynniki mogą być niestabilne i mieć duże błędy standardowe.
Miarą diagnostyczną jest VIF:
$$\mathit{VIF}_j = \frac{1}{1 - R_j^2},$$gdzie $R_j^2$ pochodzi z regresji $X_j$ na pozostałe regresory.
Obserwacje wpływowe
Niektóre obserwacje mogą silnie wpływać na wynik regresji. Analizuje się między innymi leverage $h_{ii}$, reszty studentyzowane i odległość Cooka:
$$D_i = \frac{e_i^2}{(k+1)\hat{\sigma}^2} \cdot \frac{h_{ii}}{(1-h_{ii})^2}.$$Obserwacji wpływowych nie wolno automatycznie usuwać. Najpierw trzeba sprawdzić, czy są błędami pomiaru, czy ważnymi przypadkami ekonomicznymi.
Predykcja a interpretacja przyczynowa
Dwa różne cele
Regresję można stosować do predykcji albo do wnioskowania przyczynowego. Te cele są różne.
W predykcji pytamy:
Jak dobrze model przewiduje $Y$ dla nowych obserwacji?
W interpretacji przyczynowej pytamy:
Co stałoby się z $Y$, gdybyśmy zmienili $X$, trzymając inne czynniki pod kontrolą?
Model predykcyjny może działać dobrze, choć nie ma interpretacji przyczynowej. Model przyczynowy może mieć umiarkowane $R^2$, ale odpowiadać na ważne pytanie ekonomiczne.
Pominięte zmienne
Jeżeli prawdziwy model to:
$$Y_i = \beta_0 + \beta_1 X_i + \gamma Z_i + u_i,$$a estymujemy:
$$Y_i = \alpha_0 + \alpha_1 X_i + v_i,$$to obciążenie pominiętej zmiennej w regresji prostej wynosi:
$$\mathrm{plim}(\hat{\alpha}_1) = \beta_1 + \gamma \frac{\mathrm{Cov}(X, Z)}{\mathrm{Var}(X)}.$$Jeżeli $Z$ wpływa na $Y$ i jest skorelowane z $X$, współczynnik przy $X$ przejmuje część wpływu $Z$.
Odwrócona przyczynowość
Jeżeli badamy wpływ liczby policjantów na przestępczość, regresja może pokazać dodatnią korelację. Nie oznacza to, że policjanci powodują przestępczość. Może być odwrotnie: miasta z większą przestępczością zatrudniają więcej policjantów.
Błąd pomiaru
Jeżeli zmienna $X$ jest mierzona z błędem, klasyczny błąd pomiaru zwykle osłabia współczynnik, przesuwając go w kierunku zera. To kolejny powód, dla którego regresja wymaga wiedzy o danych, nie tylko komendy komputerowej.
Kompletny projekt w R: dane symulowane
Uzasadnienie podejścia symulacyjnego
Symulacja pozwala zbadać zachowanie MNK w kontrolowanych warunkach, ponieważ prawdziwy proces generujący dane jest znany badaczowi. W przypadku danych empirycznych wartość prawdziwego parametru pozostaje nieznana.
packages <- c("tidyverse", "modelsummary", "sandwich", "lmtest", "car", "broom")
installed <- rownames(installed.packages())
for (p in packages) {
if (!(p %in% installed)) install.packages(p)
}
library(tidyverse)
library(modelsummary)
library(sandwich)
library(lmtest)
library(car)
library(broom)
set.seed(123)
n <- 500
income <- rnorm(n, mean = 6, sd = 1.5)
family_size <- sample(1:5, n, replace = TRUE)
city <- rbinom(n, size = 1, prob = 0.55)
# Prawdziwy model: konsumpcja zależy od dochodu, liczby osób i miasta.
epsilon <- rnorm(n, mean = 0, sd = 0.6)
consumption <- 0.8 + 0.55 * income + 0.35 * family_size + 0.25 * city + epsilon
dane <- tibble(consumption, income, family_size, city)
model_simple <- lm(consumption ~ income, data = dane)
model_multi <- lm(consumption ~ income + family_size + city, data = dane)
summary(model_simple)
summary(model_multi)
modelsummary(list(
"Regresja prosta" = model_simple,
"Regresja wieloraka" = model_multi
))
# Odporne błędy standardowe
coeftest(model_multi, vcov = vcovHC(model_multi, type = "HC1"))
# Diagnostyka
par(mfrow = c(2, 2))
plot(model_multi)
par(mfrow = c(1, 1))
# VIF
vif(model_multi)
# Tidy output
tidy(model_multi, conf.int = TRUE)
glance(model_multi)
augment(model_multi) |> head()
Wykresy w R
# Wykres rozrzutu i prosta regresji
ggplot(dane, aes(x = income, y = consumption)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "Dochód a konsumpcja",
subtitle = "Prosta MNK pokazuje przeciętny kierunek zależności",
x = "Dochód miesięczny w tys. zł",
y = "Konsumpcja miesięczna w tys. zł"
) +
theme_minimal()
# Wykres diagnostyczny reszt
aug <- augment(model_multi)
ggplot(aug, aes(x = .fitted, y = .resid)) +
geom_point(alpha = 0.6) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(
title = "Reszty względem wartości dopasowanych",
x = "Wartości dopasowane",
y = "Reszty"
) +
theme_minimal()
Studium przypadku I: keynesowska funkcja konsumpcji
Model ekonomiczny
Keynesowska funkcja konsumpcji opisuje zależność między dochodem rozporządzalnym a konsumpcją:
$$C_i = \beta_0 + \beta_1 Y_i + \varepsilon_i.$$Parametr $\beta_1$ interpretujemy jako krańcową skłonność do konsumpcji. Jeżeli $\beta_1 = 0{,}7$, to wzrost dochodu o 1 zł wiąże się przeciętnie ze wzrostem konsumpcji o 70 groszy.
Dane przykładowe
| Rok | Dochód | Konsumpcja | Oszczędności |
|---|---|---|---|
| 1 | 100 | 92 | 8 |
| 2 | 110 | 99 | 11 |
| 3 | 120 | 107 | 13 |
| 4 | 130 | 113 | 17 |
| 5 | 140 | 121 | 19 |
| 6 | 150 | 128 | 22 |
| 7 | 160 | 134 | 26 |
| 8 | 170 | 143 | 27 |
| 9 | 180 | 149 | 31 |
| 10 | 190 | 157 | 33 |
Kod R
cons <- tibble(
year = 1:10,
income = seq(100, 190, by = 10),
consumption = c(92, 99, 107, 113, 121, 128, 134, 143, 149, 157)
)
m_cons <- lm(consumption ~ income, data = cons)
summary(m_cons)
ggplot(cons, aes(income, consumption)) +
geom_point(size = 2) +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = "Funkcja konsumpcji",
x = "Dochód",
y = "Konsumpcja"
) +
theme_minimal()
Studium przypadku II: płace, edukacja i doświadczenie
Pytanie badawcze
Chcemy opisać zależność płac od edukacji i doświadczenia. Naturalny model logarytmiczny:
$$\ln(\mathit{Wage}_i) = \beta_0 + \beta_1 \mathit{Education}_i + \beta_2 \mathit{Experience}_i + \beta_3 \mathit{Experience}_i^2 + \varepsilon_i.$$Składnik kwadratowy pozwala na malejące przyrosty doświadczenia.
Interpretacja
Jeżeli $\beta_1 = 0{,}08$, to dodatkowy rok edukacji wiąże się w przybliżeniu z płacą wyższą o 8%, przy tym samym doświadczeniu.
Efekt krańcowy doświadczenia:
$$\frac{\partial E(\ln \mathit{Wage} \mid X)}{\partial \mathit{Experience}} = \beta_2 + 2\beta_3 \mathit{Experience}.$$Jeżeli $\beta_2 > 0$ i $\beta_3 < 0$, doświadczenie zwiększa płace, ale coraz wolniej.
set.seed(321)
n <- 1000
education <- sample(10:20, n, replace = TRUE)
experience <- pmax(0, round(rnorm(n, 12, 6)))
female <- rbinom(n, 1, 0.5)
u <- rnorm(n, 0, 0.25)
log_wage <- 1.8 + 0.075 * education + 0.045 * experience -
0.0012 * experience^2 - 0.10 * female + u
wage <- exp(log_wage)
wages <- tibble(wage, log_wage, education, experience, female)
m_wage <- lm(log(wage) ~ education + experience + I(experience^2) + female,
data = wages)
summary(m_wage)
modelsummary(m_wage)
Studium przypadku III: ceny mieszkań
Regresja hedoniczna
Regresja hedoniczna opisuje cenę dobra jako funkcję jego cech. Dla mieszkań:
$$\ln(\mathit{Price}_i) = \beta_0 + \beta_1 \mathit{Area}_i + \beta_2 \mathit{Rooms}_i + \beta_3 \mathit{Distance}_i + \beta_4 \mathit{New}_i + \varepsilon_i.$$Cena zależy od powierzchni, liczby pokoi, odległości od centrum, standardu budynku, piętra, dzielnicy i wielu innych cech.
Logarytmiczna transformacja ceny
Ceny mieszkań są dodatnie i zwykle prawoskośne. Logarytm stabilizuje wariancję i pozwala interpretować współczynniki jako przybliżone zmiany procentowe.
set.seed(2026)
n <- 800
area <- runif(n, 25, 120)
rooms <- pmax(1, round(area / 25 + rnorm(n, 0, 0.7)))
distance <- runif(n, 1, 25)
new_building <- rbinom(n, 1, 0.35)
noise <- rnorm(n, 0, 0.12)
log_price <- 11.2 + 0.012 * area + 0.035 * rooms -
0.018 * distance + 0.09 * new_building + noise
price <- exp(log_price)
flats <- tibble(price, area, rooms, distance, new_building)
m_flat <- lm(log(price) ~ area + rooms + distance + new_building,
data = flats)
summary(m_flat)
coeftest(m_flat, vcov = vcovHC(m_flat, type = "HC1"))
Pełna procedura badawcza
Procedura badawcza
Rzetelna analiza regresji nie zaczyna się od wywołania funkcji lm() w R. Zaczyna się od sformułowania pytania badawczego. Procedura:
- Zdefiniuj zmienną objaśnianą.
- Zdefiniuj główną zmienną objaśniającą.
- Wyjaśnij, dlaczego oczekujesz określonego znaku współczynnika.
- Wskaż zmienne kontrolne.
- Ustal, czy interesuje Cię predykcja, opis, czy interpretacja przyczynowa.
- Oszacuj model bazowy.
- Sprawdź diagnostykę reszt i obserwacji wpływowych.
- Sprawdź odporność wyniku na alternatywne specyfikacje.
- Zinterpretuj współczynniki w jednostkach ekonomicznych.
- Sformułuj ograniczenia badania.
Raportowanie wyników
W pracy akademickiej samo przedstawienie wydruku funkcji summary() jest niewystarczające. Wyniki należy opisać i zinterpretować. Dobra tabela regresji powinna zawierać:
- estymowane współczynniki;
- błędy standardowe;
- liczbę obserwacji;
- $R^2$ i skorygowane $R^2$;
- informację o typie błędów standardowych;
- informację o zmiennych kontrolnych;
- jasny opis zmiennej zależnej.
modelsummary(
list(
"Model 1" = model_simple,
"Model 2" = model_multi
),
statistic = "({std.error})",
stars = TRUE,
gof_map = c("nobs", "r.squared", "adj.r.squared")
)
Najczęstsze błędy i błędne intuicje
Błąd: współczynnik zawsze oznacza przyczynę
Regresja mierzy zależność warunkową. Przyczynowość wymaga dodatkowych założeń. Jeżeli badamy wpływ edukacji na płace i pomijamy zdolności, współczynnik przy edukacji może być obciążony.
Błąd: wysoki poziom istotności oznacza ważność ekonomiczną
W dużych próbach bardzo małe efekty mogą być statystycznie istotne. Ważność ekonomiczna wymaga interpretacji wielkości efektu.
Błąd: brak istotności oznacza brak efektu
Brak istotności może wynikać z małej próby, dużej wariancji, współliniowości lub błędnego pomiaru. Nie jest automatycznym dowodem braku zależności.
Błąd: usuwanie obserwacji odstających bez uzasadnienia
Obserwacje odstające mogą być błędami, ale mogą też być najważniejszymi obserwacjami w danych. Decyzja o ich usunięciu musi być udokumentowana.
Błąd: mechaniczne dodawanie zmiennych
Więcej zmiennych nie zawsze oznacza lepszy model. Zmienna kontrolna może być mediatorem, koliderem albo źródłem dodatkowego błędu interpretacyjnego.
Zadania z odpowiedziami
Zadanie 1: regresja prosta
Dane:
$$x = (2, 4, 6, 8), \qquad y = (3, 5, 7, 10).$$Oblicz $\hat{\beta}_0$ i $\hat{\beta}_1$.
Odpowiedź
Średnie: $\bar{x} = 5$, $\bar{y} = 6{,}25$.
Licznik:
$$\sum(x_i - \bar{x})(y_i - \bar{y}) = (-3)(-3{,}25) + (-1)(-1{,}25) + (1)(0{,}75) + (3)(3{,}75) = 23.$$Mianownik:
$$\sum(x_i - \bar{x})^2 = 9 + 1 + 1 + 9 = 20.$$Nachylenie:
$$\hat{\beta}_1 = 23/20 = 1{,}15.$$Wyraz wolny:
$$\hat{\beta}_0 = 6{,}25 - 1{,}15 \cdot 5 = 0{,}5.$$Zadanie 2: interpretacja logarytmu
W modelu:
$$\ln(\mathit{Wage}_i) = 1{,}5 + 0{,}09 \cdot \mathit{Education}_i + \varepsilon_i$$zinterpretuj współczynnik przy edukacji.
Odpowiedź
Dodatkowy rok edukacji wiąże się w przybliżeniu z płacą wyższą o 9%. Dokładna zmiana procentowa wynosi:
$$100(e^{0{,}09} - 1) \approx 9{,}42\%.$$Zadanie 3: $R^2$
Jeżeli $TSS = 500$ i $RSS = 125$, oblicz $R^2$.
Odpowiedź
$$R^2 = 1 - \frac{RSS}{TSS} = 1 - \frac{125}{500} = 0{,}75.$$Model wyjaśnia 75% zmienności $Y$ w próbie.
Zadanie 4: test $t$
Estymator wynosi $\hat{\beta} = 0{,}8$, błąd standardowy $0{,}2$. Testujemy $H_0: \beta = 0$. Oblicz statystykę $t$.
Odpowiedź
$$t = \frac{0{,}8 - 0}{0{,}2} = 4.$$Zadanie 5: pominięta zmienna
Prawdziwy model zawiera zdolności, ale badacz estymuje płace tylko na edukacji. Wyjaśnij kierunek obciążenia, jeśli zdolności dodatnio wpływają na płace i są dodatnio skorelowane z edukacją.
Odpowiedź
Obciążenie jest dodatnie. Współczynnik przy edukacji przejmuje część dodatniego wpływu zdolności, ponieważ osoby lepiej wykształcone mają przeciętnie wyższe zdolności, a zdolności zwiększają płace.
Słownik pojęć
Regresja liniowa — Model opisujący zmienną objaśnianą jako liniową funkcję parametrów i zmiennych objaśniających.
MNK — Metoda najmniejszych kwadratów, czyli procedura minimalizacji sumy kwadratów reszt.
Reszta — Różnica między obserwacją a wartością dopasowaną.
Składnik losowy — Nieobserwowana część modelu obejmująca pominięte czynniki i zaburzenia losowe.
Egzogeniczność — Warunek, że składnik losowy ma wartość oczekiwaną zero warunkowo względem regresorów.
Homoskedastyczność — Stała wariancja składnika losowego.
Autokorelacja — Korelacja składników losowych między obserwacjami, szczególnie częsta w szeregach czasowych.
Współliniowość — Silna zależność między zmiennymi objaśniającymi.
$R^2$ — Udział zmienności zmiennej objaśnianej wyjaśniony przez model.
BLUE — Najlepszy liniowy nieobciążony estymator.
Bibliografia
- Aitken, A. C. (1935). On Least Squares and Linear Combinations of Observations. Proceedings of the Royal Society of Edinburgh, 55, 42–48.
- Angrist, J. D., Pischke, J.-S. (2009). Mostly Harmless Econometrics. Princeton University Press.
- Angrist, J. D., Pischke, J.-S. (2015). Mastering Metrics. Princeton University Press.
- Davidson, R., MacKinnon, J. G. (2004). Econometric Theory and Methods. Oxford University Press.
- Galton, F. (1886). Regression Towards Mediocrity in Hereditary Stature. Journal of the Anthropological Institute, 15, 246–263.
- Gauss, C. F. (1809). Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium. Hamburg.
- Goldberger, A. S. (1991). A Course in Econometrics. Harvard University Press.
- Greene, W. H. (2018). Econometric Analysis. Pearson.
- Hayashi, F. (2000). Econometrics. Princeton University Press.
- Hansen, B. E. (2022). A Modern Gauss-Markov Theorem. Econometrica, 90(3), 1283–1294.
- Legendre, A.-M. (1805). Nouvelles méthodes pour la détermination des orbites des comètes. Paris: F. Didot.
- Stock, J. H., Watson, M. W. (2020). Introduction to Econometrics. Pearson.
- Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data. MIT Press.
- Wooldridge, J. M. (2020). Introductory Econometrics: A Modern Approach. Cengage.