[2019-05-22] Всё написанное ниже является пережитком эпохи. От читателя требуется понимать публикацию с точностью до наоборот.
В эконометрике очень полезно, если ошибки имеют нормальное распределение. В таком случае эконометристы сразу могут получить надёжные оценки для ошибок коэффициентов и интервалы для коэффициентов. При этом студентов обычно коварно успокаивают, мол, «если у вас очень много наблюдений, то тогда нормальность появится сама». Хо-хо, чёрта с два. Если у меня величина имеет \(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
А вот и экспонента: 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
Так, это хорошо. Теперь представим, что после долгих лет упорного дата-майнинга (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
На диаграмме разброса покажем только первые десять тысяч точек, так как оно всё рябит.
Выполнены почти все условия теоремы Гаусса—Маркова ((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
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
Очень похоже на распределение Стьюдента, но я не могу знать наверняка, так как и пик острее и выше, и «хвосты» тяжелее.
Для полного счастья запишем, как будут выглядеть истинные ошибки в модели с функцией Кобба—Дугласа и прибавленной ошибкой:
\[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))
Вообще говоря, ничто не может быть нагляднее, чем график. Особенно трёхмерный. Конечно, я не могу построить график с миллионом точек. Даже сохранить данные не могу (программа зависла, пришлось к данному моменту всё заново делать... чёрт). Поэтому построим график с уравнениями конечных регрессий по выборке размером в две тысячи наблюдений.
Для «правильной» регрессии l_Y график выглядит так:
Видно, что «пух» есть, что сильный, что на миллионе наблюдений он вообще превращается в персидский ковёр, но... Жизнь есть жизнь.
Для регрессии с неправильно специфицированными ошибками график выглядит так:
«Пуха» меньше, но и цена этого выше.
Вывод. Допущение о «нормальности» этого мира — абсурд. Несмотря на «нормальность» распределения ростов, на стьюдентовость распределения масс и на то, что уже на двадцати пяти монетках можно получить нормальное распределение, глупо, смешно и опасно предполагать, что парочка безобидных преобразований просто увеличит стандартные ошибки. Даже если это утверждал сам Андрей Александрович Мамонтов, посмеиваясь над скептически настроенными студентами. Правильная спецификация регрессоров даже асимптотически не может исцелить от ненормальности ошибок. 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
Постсубскриптум. Вы только что прослушали «Фантазию» для двух фортепиано А. Н. Скрябина не сравнение двух моделей, а сравнение двух спецификаций, так как в обоих случаях регрессоры были одними и теми же.
