Почему так важна нормальность ошибок в регрессии

В эконометрике очень полезно, если ошибки имеют нормальное распределение. В таком случае эконометристы сразу могут получить надёжные оценки для ошибок коэффициентов и интервалы для коэффициентов. При этом студентов обычно коварно успокаивают, мол, «если у вас очень много наблюдений, то тогда нормальность появится сама». Хо-хо, чёрта с два. Если у меня величина имеет \(t\)-распределение или, того хуже, равномерное, то нормальности у меня автоматически никакой не появится. Тогда гегемоны-семинаристы говорят: «Ой, ну ладно вам придираться, был доверительный интервал шириной 2,8 — стал шириной 2,9».

Сегодня я покажу, что отсутствие нормальности приводит к куда более страшным последствиям.

Рассмотрим функцию Кобба—Дугласа, и пусть, как всегда, отдача от труда равна \(\frac13\), от капитала — \(\frac23\), производственная функция описывается формулой
\[Y = A \cdot L^\alpha \cdot K^\beta \cdot f(\varepsilon). \]

Очень важно определить, что такое \(f(\varepsilon)\). Если мы хотим получить линейную модель вида
\[\ln Y =\ln A + \alpha \ln L + \beta \ln K +\varepsilon, \]
то при потенцировании должно быть
\[Y = A \cdot L^\alpha \cdot K^\beta \cdot e^{\varepsilon}.\]

Если \(\varepsilon\sim \mathcal{N}(0; 1)\), то как выглядит распределение \(e^{\varepsilon}\)? Давайте посмотрим.

В gretl сгенерируем Add — Random Variable — Normal со средним 0 и стандартным отклонением 1 (команда randgen(N,0,1)). Пускай для надёжности это будет десять тысяч, сто тысяч, миллион случайных величин! Нам же всего лишь график посмотреть, правда?

Summary statistics, using the observations 1 - 1000000
for the variable 'epsilon' (1000000 valid observations)

  Mean                     -0.0013720
  Median                  -0.00026809
  Minimum                     -4.9740
  Maximum                      4.7117
  Standard deviation           1.0007
  C.V.                         729.32
  Skewness                 -0.0015528
  Ex. kurtosis              0.0075161
  5% percentile               -1.6493
  95% percentile               1.6438
  Interquartile range          1.3480
  Missing obs.                      0

Test for normality of epsilon:
 Doornik-Hansen test = 2.76293, with p-value 0.251211
 Shapiro-Wilk W = 0.999998, with p-value 0.999982
 Lilliefors test = 0.000685939, with p-value ~= 0.3
 Jarque-Bera test = 2.75571, with p-value 0.252118

Normal distribution

Normal distribution 2

А вот и экспонента: Add — Define new variable: expepsilon=exp(epsilon).

Summary statistics, using the observations 1 - 1000000
for the variable 'expepsilon' (1000000 valid observations)

  Mean                         1.6477
  Median                      0.99973
  Minimum                   0.0069156
  Maximum                      111.24
  Standard deviation           2.1661
  C.V.                         1.3146
  Skewness                     6.1000
  Ex. kurtosis                 91.825
  5% percentile               0.19218
  95% percentile               5.1750
  Interquartile range          1.4510
  Missing obs.                      0

Exponent of normal distribution

Exponent of normal distribution 2

Так, это хорошо. Теперь представим, что после долгих лет упорного дата-майнинга (data mining) мы достали данные о миллионе рабочих и миллионе рабочих мест, оборудованных капиталом. Теперь сгенерируем рабочую силу и капитал для функции Кобба—Дугласа. Пусть они все будут равномерны, труд — от 3 до 15 (randgen(u,3,15)), капитал — от 5 до 20 (randgen(u,5,20)). Ничего, регрессия только выиграет от такой независимости. Add — Random variable — Uniform. Почему такие значения? Во-первых, все логарифмы больше единицы. Во-вторых, чем шире разброс регрессоров, тем лучше. Корреляцию необходимо тоже проверить статистически.

Да, и не забываем менять seed! Tools — Seed for random numbers. А ещё лучше записывать, какие числа используем.

Summary statistics, using the observations 1 - 1000000
for the variable 'L' (1000000 valid observations)
  Mean                         8.9972
  Median                       8.9914
  Minimum                      3.0000
  Maximum                      15.000
  Standard deviation           3.4639
  C.V.                        0.38499
  Skewness                  0.0028195
  Ex. kurtosis                -1.2002
  5% percentile                3.6035
  95% percentile               14.402
  Interquartile range          5.9986
  Missing obs.                      0

Summary statistics, using the observations 1 - 1000000
for the variable 'K' (1000000 valid observations)
  Mean                         12.500
  Median                       12.501
  Minimum                      5.0000
  Maximum                      20.000
  Standard deviation           4.3319
  C.V.                        0.34655
  Skewness                4.5311e-006
  Ex. kurtosis                -1.2022
  5% percentile                5.7482
  95% percentile               19.247
  Interquartile range          7.5089
  Missing obs.                      0

corr(L, K) = -0.00050070
Under the null hypothesis of no correlation:
 t(999998) = -0.500703, with two-tailed p-value 0.6166

Uniform factor distribution 1

Uniform factor distribution 2

На диаграмме разброса покажем только первые десять тысяч точек, так как оно всё рябит.

Factor distribution

Выполнены почти все условия теоремы Гаусса—Маркова ((1) отсутствие корреляции ошибок благодаря качественному генератору случайных величин (статистика Дарбина—Уотсона), (2) детерминированность и независимость \(X_i\), (3) отсутствие систематической ошибки \((\mathbb{E}(\varepsilon)=0)\), (4) гомоскедастичность ошибок (в силу, опять же, генератора)) — осталось только заняться спецификацией.

Всё отлично. Теперь придумываем модель. Пусть
\[Y = 2{,}71828 \cdot L^{0{,}333} \cdot K^{0{,}667} \cdot e^{\varepsilon},\]
\[Z = 2{,}71828 \cdot L^{0{,}333} \cdot K^{0{,}667} + 5\varepsilon.\]

В данном случае мы немного чувствуем себя богами, и это оттого, что нам известны истинные \(\varepsilon\). Что же, генерируем желаемые переменные в Add — Define new variable:

Y=2.71828*L^0.333*K^0.667*expepsilon
Z=abs(2.71828*L^0.333*K^0.667+5*epsilon)

Почему я применил abs() в последней формуле? Потому что кое-где эти дурацкие нормальные ошибки уводят наш чудесный выпуск в минус. На миллион наблюдений и при моих параметрах модели таких случаев едва ли наберётся сто, поэтому можно успешно их увести в положительную область. Всё равно они по модулю маленькие — действительность отражают неплохо.

Смотреть там особо не на что, пока мы не получим результаты регрессионного анализа. Выделяем все четыре. Нажимаем Variables — Add — Logs of selected variables. Теперь идём в Model — Ordinary Least Squares и задаём параметры. Зависим l_Y, независимы l_K и l_L.

Model 1: OLS, using observations 1-1000000
Dependent variable: l_Y
             coefficient   std. error   t-ratio   p-value
  -------------------------------------------------------
  const       1.00204      0.00811912    123.4    0.0000  ***
  l_L         0.333701     0.00229338    145.5    0.0000  ***
  l_K         0.665010     0.00261911    253.9    0.0000  ***

Mean dependent var   3.340627   S.D. dependent var   1.042599
Sum squared resid     1001331   S.E. of regression   1.000667
R-squared            0.078822   Adjusted R-squared   0.078820
F(2, 999997)         42783.15   P-value(F)           0.000000
Log-likelihood       -1419603   Akaike criterion      2839213
Schwarz criterion     2839248   Hannan-Quinn          2839223
Log-likelihood for Y = -4.76023e+006

Всё отлично. Правда, R² низкий. Но это ещё ничего не значит. Как видите, этот тест отлично показал, что R² — это вообще не индикатор. Это регрессия начинает извиняться, что в каждом наблюдении у неё есть ошибка, а ошибок много (очень много), а они ещё все скачут экспоненциально (вернее, дают большую ошибку в оригинальный \(Y\) и оставляют этот гадкий след даже после логарифмирования). Коэффициенты, конечно же, точные. Дополнительно исчисленный Дарбин—Уотсон (в таблице не показан) равен 1,998 792. Теперь вторая модель.

Model 2: OLS, using observations 1-1000000
Dependent variable: l_Z
             coefficient   std. error    t-ratio   p-value
  --------------------------------------------------------
  const       0.860155     0.00174908     491.8    0.0000  ***
  l_L         0.350244     0.000494058    708.9    0.0000  ***
  l_K         0.700475     0.000564230   1241      0.0000  ***

Mean dependent var   3.320815   S.D. dependent var   0.375988
Sum squared resid    46470.85   S.E. of regression   0.215571
R-squared            0.671274   Adjusted R-squared   0.671273
F(2, 999997)          1021021   P-value(F)           0.000000
Log-likelihood       115526.5   Akaike criterion    -231047.0
Schwarz criterion   -231011.5   Hannan-Quinn        -231037.2
Log-likelihood for Z = -3.20529e+006

Что?! Что я вижу?! На миллионе наблюдений в допущении экспоненты нормальных ошибок мы получили асимптотически смещённые оценки! Позор! Позор! Позор этому миру. Можно уже начинать сжигать книги.

Проведём тест на сумму коэффициентов. Я понимаю, что на миллионе в первой регрессии и так уже коэффициенты близки к истинным, поэтому не так хорошо видно, что́ я на самом деле пытаюсь понять, ибо \(0{,}333\,7+0{,}665\,0\approx\ldots1\) — а насколько именно приближённо? Вот на эту единицу мы и проверим. Tests — Linear restrictions. Вводим мантру:

b[l_L]+b[l_K]=1

Ответ для регрессии \(Y\):

Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 999997) = 0.137098, with p-value = 0.711183

Restricted estimates:
             coefficient   std. error   t-ratio   p-value
  -------------------------------------------------------
  const       0.999065     0.00116651    856.5    0.0000  ***
  l_L         0.334261     0.00172455    193.8    0.0000  ***
  l_K         0.665739     0.00172455    386.0    0.0000  ***
  Standard error of the regression = 1.00067

Смотрим p-value: 0,711 — конечно, при миллионе-то наблюдений! Отлично, тест пройден.

А если ввести линейные ограничения в модели с кривыми ошибками (регрессия \(Z\))?

Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 999997) = 4569, with p-value = 0

Restricted estimates:
             coefficient   std. error    t-ratio   p-value
  --------------------------------------------------------
  const       0.977157     0.000251871   3880      0.0000  ***
  l_L         0.328230     0.000372363    881.5    0.0000  ***
  l_K         0.671770     0.000372363   1804      0.0000  ***
  Standard error of the regression = 0.216063

Удивительно, но наши смещённые коэффициенты чудесным образом «пофиксились» самостоятельно! Другое дело, что p-value равно нулю, что мы улетели в «хвост» распределения, гипотеза о равенстве суммы единице отвергается (ошибка первого рода), поэтому данный тест показывает нам жирную дулю, что ненормальные ошибки — бич регрессий.

Давайте сохраним остатки (Save — Residuals) и посмотрим на их распределение. Заодно и потестируем. Я назвал их eY и eZ соответственно.

Test for normality of eY:
 Doornik-Hansen test = 2.76177, with p-value 0.251356
 Shapiro-Wilk W = 0.999998, with p-value 0.999989
 Lilliefors test = 0.000695252, with p-value ~= 0.28
 Jarque-Bera test = 2.75455, with p-value 0.252266

Good residuals 1

Good residuals 2

Test for normality of eZ:
 Doornik-Hansen test = 994212, with p-value 0
 Shapiro-Wilk W = 0.913314, with p-value 1.50967e-149
 Lilliefors test = 0.0623541, with p-value ~= 0
 Jarque-Bera test = 2.59783e+007, with p-value 0

Bad residuals 1

Bad residuals 2

Очень похоже на распределение Стьюдента, но я не могу знать наверняка, так как и пик острее и выше, и «хвосты» тяжелее.

Для полного счастья запишем, как будут выглядеть истинные ошибки в модели с функцией Кобба—Дугласа и прибавленной ошибкой:
\[Y = AL^{\alpha}K^{\beta} + \varepsilon = AL^{\alpha}K^{\beta}\left(1 + \frac{\varepsilon}{AL^{\alpha}K^{\beta}} \right)\]

Так как последняя скобка довольно близка к единице (уж слишком на многое делят несчастный эпсилон), то действует выражение эквивалентности \(x \stackrel{x\rightarrow0}{\sim} \ln (1+x)\). Но я есмь человек, любящий точность, а оттого посмотрю на оба распределения. Как видно ниже, они почти сливаются.

error1=epsilon/(2.71828*L^0.333*K^0.667)
error2=ln(1+epsilon/(2.71828*L^0.333*K^0.667))

error1
error2
Bad error

Вообще говоря, ничто не может быть нагляднее, чем график. Особенно трёхмерный. Конечно, я не могу построить график с миллионом точек. Даже сохранить данные не могу (программа зависла, пришлось к данному моменту всё заново делать... чёрт). Поэтому построим график с уравнениями конечных регрессий по выборке размером в две тысячи наблюдений.

Для «правильной» регрессии l_Y график выглядит так:
Volumetric normal distribution 1

Volumetric normal distribution 2

Volumetric normal distribution 3

Видно, что «пух» есть, что сильный, что на миллионе наблюдений он вообще превращается в персидский ковёр, но... Жизнь есть жизнь.

Для регрессии с неправильно специфицированными ошибками график выглядит так:

Volumetric not normal distribution

Volumetric not normal distribution 2

Volumetric not normal distribution 3

«Пуха» меньше, но и цена этого выше.

Вывод. Допущение о «нормальности» этого мира — абсурд. Несмотря на «нормальность» распределения ростов, на стьюдентовость распределения масс и на то, что уже на двадцати пяти монетках можно получить нормальное распределение, глупо, смешно и опасно предполагать, что парочка безобидных преобразований просто увеличит стандартные ошибки. Даже если это утверждал сам Андрей Александрович Мамонтов, посмеиваясь над скептически настроенными студентами. Правильная спецификация регрессоров даже асимптотически не может исцелить от ненормальности ошибок. Memento mori.

Постскриптум. «Иногда малая выборка лучше большой» — это верно или неверно? Неверно, но доказательство очень трудное. Тогда я скажу по-другому. Можно ли подпортить картину торжества модели с правильной спецификацией? Можно. Построим при помощи МНК те же две регрессии, но только на выборке из 800 наблюдений. Так как выборка — случайная величина, то и бояться её нечего.

Model 1: OLS, using observations 1-800
Dependent variable: l_Y
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       0.755572     0.294980      2.561    0.0106    **
  l_L         0.309048     0.0845139     3.657    0.0003    ***
  l_K         0.796715     0.0943801     8.442    1.47e-016 ***

Mean dependent var   3.375052   S.D. dependent var   1.057215
Sum squared resid    806.8346   S.E. of regression   1.006151
R-squared            0.096536   Adjusted R-squared   0.094269
F(2, 797)            42.58025   P-value(F)           2.69e-18
Log-likelihood      -1138.554   Akaike criterion     2283.107
Schwarz criterion    2297.161   Hannan-Quinn         2288.506
Log-likelihood for Y = -3838.6

Model 2: OLS, using observations 1-800
Dependent variable: l_Z
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       0.786741     0.0624576     12.60    2.61e-033 ***
  l_L         0.348992     0.0178945     19.50    1.48e-069 ***
  l_K         0.732322     0.0199836     36.65    4.14e-173 ***

Mean dependent var   3.333156   S.D. dependent var   0.379279
Sum squared resid    36.17173   S.E. of regression   0.213037
R-squared            0.685294   Adjusted R-squared   0.684505
F(2, 797)            867.7622   P-value(F)           8.2e-201
Log-likelihood       103.3827   Akaike criterion    -200.7654
Schwarz criterion   -186.7115   Hannan-Quinn        -195.3665
Log-likelihood for Z = -2563.14

Какая регрессия кажется лучше? Коэффициенты одинаково кривоваты и почти равны (и это при 800 наблюдениях! Вы там думали, что склеить пять наборов по 24 точки — это асимптотика, но нет, взгляните горькой правде в глаза), но у второй куда больше дьявольский R², да и ошибки коэффициентов поменьше... Однако отсутствие асимптотики не мешает второй модели по-прежнему запарывать тест линейного ограничения:

(Регрессия Y)
Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 797) = 0.702842, with p-value = 0.402082

Restricted estimates:
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       1.00047      0.0410358    24.38     1.49e-098 ***
  l_L         0.262031     0.0632130     4.145    3.76e-05  ***
  l_K         0.737969     0.0632130    11.67     3.49e-029 ***
  Standard error of the regression = 1.00596

(Регрессия Z)
Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 797) = 9.2669, with p-value = 0.00241004

Restricted estimates:
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       0.975023     0.00873522   111.6     0.0000    ***
  l_L         0.312843     0.0134561     23.25    1.05e-091 ***
  l_K         0.687157     0.0134561     51.07    1.13e-253 ***
  Standard error of the regression = 0.214138

Постсубскриптум. Вы только что прослушали «Фантазию» для двух фортепиано А. Н. Скрябина не сравнение двух моделей, а сравнение двух спецификаций, так как в обоих случаях регрессоры были одними и теми же.

About Andreï Kostyrka

Науколюб, грамматический нацист, антитеист. Пишу стихотворения, сочиняю музыку, верстаю книги, занимаюсь эконометрикой и настраиваю фортепиано.
Bookmark the permalink.

Leave a Reply

Your email address will not be published. Required fields are marked *