Матожидание ошибок при симуляции МНК

Я больше года мучился, не понимал, почему аддитивность/мультипликативность нормальных ошибок играет такую большую роль, и даже разъяснение преподавателя не помогало. Почти год назад я написал этот пост: http://kostyrka.ru/blog/archives/1246 — и, оказывается, впустую! Сегодня на эконометрике ко мне пришло видение этой проблемы в свете того, что матожидание от функции не равно функции от матожидания. Конкретно это имеет значение, если вас просят оценить уравнение, в котором ошибка мультипликативна и, следовательно, при логарифмировании её матожидание должно дать ноль. Что делает плохой студент? Генерирует мультипликативную ошибку со средним 1 из первого попавшегося распределения!

Предположим, вы хотите научиться оценивать уравнение
\[Q = a K^\alpha L^\beta \cdot \varepsilon,\]
где \(\mathbb{E}(\varepsilon\mid K, L)=1\). Уже и компьютер включили, и gretl запустили. Какое распределение \(\varepsilon\) надо взять?

По прочтении условия очень хочется взять распределение ошибок \(\chi^2_1\), так как \(\mathbb{E}(\chi^2_1)=1\). Ещё очень хочется взять \(\mathcal{N}(1;\sigma^2)\). Некоторые кулибины берут \(\mathcal{U}[0{,}5;1{,}5]\) или \(\mathcal{U}[0;2]\). Так почему же нельзя брать ни то, ни другое, ни третье?

Ответ: так как по МНК уравнение будет оцениваться в логарифмах, то важно, чтобы не \(\mathbb{E}(\varepsilon)\) было единицей, а \(\mathbb{E}(\ln \varepsilon)\) было нулём. Матожидание не всепроникающее, поэтому матожидание логарифма не равно логарифму матожидания. Если бы выполнялась такая глупость, что \(\mathbb{E}\bigl(f(X)\bigr) = f \bigl( \mathbb{E}(X) \bigr)\), то было бы и \(\mathbb{E}(\ln X) = \ln \mathbb{E}(X)\), и \(\mathbb{E}(X^2) = \bigl(\mathbb{E}(X)\bigr)^2\), дисперсия стала бы тождественным нулём, и наступил бы конец света.

 

Хи-квадрат после логарифмирования начнёт себя вести очень плохо. Настолько плохо, что формулу для его матожидания я запишу в попустительской нотации:
\[ \varepsilon \sim \chi^2_1 \quad \Rightarrow \quad \mathbb{E}(\ln \varepsilon) = \mathbb{E}\bigl(\ln \mathcal{N}^2(0;1)\bigr) =\mathbb{E}\bigl(2\ln |\mathcal{N}(0;1)|\bigr) = 2\mathbb{E}\bigl(\ln |\mathcal{N}(0;1)|\bigr) \]

Во-первых, отрицательные значения перейдут вправо (распределение аргумента станет положительным). Во-вторых, из оставшихся логарифмов ошибок более 2/3 будут отрицательными (вероятность, что величина из нормального распределения по модулю будет меньше единицы, равна 68,27 %), причём некоторые из них будут убийственно отрицательными: логарифм близкой к нулю величины уходит глубоко под землю. Если уж идти до конца, то
\[ \varepsilon \sim \chi^2_1 \quad \Rightarrow \quad \mathbb{E}(\ln \varepsilon) = 2\mathbb{E}\bigl(\ln |\mathcal{N}(0;1)|\bigr) = \sqrt{\frac{2}{\pi}} \int\limits_{-\infty}^{+\infty} \ln x \cdot e^{ -\frac{x^2}{2}} \,\mathrm{d}x = -\gamma -\ln 2 \approx -1{,}270\,36, \]
где \(\gamma\) — постоянная Эйлера—Маскерони.

В общем случае матожидание логарифма хи-квадрата с \(k\) степенями свободы выражается через полигамма-функцию (даже дигамма-функцию):
\[ \mathbb{E}(\chi^2_k) = \ln 2 + \psi^{(0)}\left(\frac{k}{2} \right), \quad \psi^{(0)}(x) \equiv \frac{\Gamma'(x)}{\Gamma(x)} \]
Biased log of chi square
Если хотите, чтобы матожидание логарифма было равно нулю, то надо брать хи-квадрат с \(1{,}866\,025\) степенями свободы.

 

Нормальное распределение может выдать абсолютно любые значения случайной величины, так как функция плотности определена на \(\mathbb{R}\). Поэтому если ошибка мультипликативна и нормальна, то, вообще говоря, пропадают все значения, где ошибка получилась меньше нуля. Кроме того, если ошибка распределена как \(\mathcal{N}(1;\sigma^2)\), а наивные пользователи думают, что её логарифм в среднем даст ноль, то спешу их разочаровать следующей картиной:
Biased Log of Normal Distribution
Здесь avgleps — это \(\overline{\ln \varepsilon_i} = \mathbb{E}(\ln\varepsilon_i \mid \sigma)\), где \(\varepsilon_i \sim \mathcal{N}(1;\sigma^2)\). Любой желающий может в этом убедиться (если скрипт виснет, надо уменьшить размер выборки до, скажем, nulldata 5000):
nulldata 100000
scalar step=0.01
#
loop for (se=0.001; se<10; se+=step) --progressive series eps = randgen(N,1,se) series leps = ln(eps) genr avgleps = mean(leps) store biasedln.gdt se avgleps endloop # open biasedln.gdt gnuplot avgleps se --with-lines --output=display --suppress-fitted

На самом деле безусловное матожидание логарифма нормального распределения — это следующее аналитическое выражение:
\[ \varepsilon \sim \mathcal{N}(1;\sigma^2) \quad \Rightarrow \quad \mathbb{E}\bigl(\ln \varepsilon \bigr) = \frac{1}{2} \left(-{}_1\mathrm{F}_1\left(0;\frac{1}{2};-\frac{1}{2 \sigma ^2}\right)+\ln \left(\frac{\sigma ^2}{2}\right)-\gamma \right) + \frac12 i \pi \mathrm{erfc}\left(\frac{1}{\sqrt{2} \sigma}\right), \]
где \({}_1\mathrm{F}_1\) — вырожденная гипергеометрическая функция Куммера (Kummer confluent hypergeometric function).
Biased expectation
Чтобы избавиться от мнимой единицы, некоторые товарищи могут взять условное матожидание (\(x>0\)) или выкинуть все отрицательные значения (как на графике в gretl’е), то и тогда у них матожидание будет равно ещё более жуткому выражению (раскрываю свои карты) — и притом обе эти функции будут давать абсолютно одинаковые графики:
\[ \text{Expectation}[\log (x)\unicode{f3d3}x>0,x\approx \text{NormalDistribution}[1,\sigma ]] \]
\[ \frac{e^{-\frac{1}{2 \sigma ^2}} \left(-e^{\frac{1}{2 \sigma ^2}} \sigma \text{Hypergeometric1F1}^{(1,0,0)}\left(0,\frac{1}{2},-\frac{1}{2 \sigma ^2}\right)+\sqrt{\frac{2}{\pi }} \text{Hypergeometric1F1}^{(1,0,0)}\left(1,\frac{3}{2},\frac{1}{2 \sigma ^2}\right)+\gamma e^{\frac{1}{2 \sigma ^2}} \sigma \text{erfc}\left(\frac{1}{\sqrt{2} \sigma }\right)-e^{\frac{1}{2 \sigma ^2}} \sigma \log (2) \text{erfc}\left(\frac{1}{\sqrt{2} \sigma }\right)-2 e^{\frac{1}{2 \sigma ^2}} \sigma \text{erfc}\left(\frac{1}{\sqrt{2} \sigma }\right) \log (\sigma )-2 \gamma e^{\frac{1}{2 \sigma ^2}} \sigma +4 e^{\frac{1}{2 \sigma ^2}} \sigma \log (\sigma )\right)}{2 \sigma \left(\text{erf}\left(\frac{1}{\sqrt{2} \sigma }\right)+1\right)} \]

Вам это надо? Не надо. Симуляция уже показала вам, насколько велико смещение. Если \(\sigma\not\approx 1{,}092\,340\) (приблизительное решение этой страшной аналитической вещи относительно \(\sigma\)), то наличествует систематическая ошибка, и коэффициенты модели неверны, причём по вышеприведённым графикам можно оценить смещение. Даже если смещение равно 0,1, но наблюдений в выборке 100 000, то...

Совсем плохой случай: если в исходной модели ошибки мультипликативны и нормальны с нулевым матожиданием и дисперсией \(\sigma^2\), то при переходе к логарифмам их матожидание станет равным \(\frac{1}{2} \left(-\gamma + \ln \frac{\sigma^2}{2} \right)\), смещение принимает какое угодно значение, смещается оценка свободного члена, и ни при какой огромной выборке остальные истинные коэффициенты получить не удастся!

 

Если кто-то возьмёт равномерное распределение ошибок с матожиданием 1, то ясно, что для получения осмысленных ошибок нижняя граница интервала распределения может изменяться от 0 до 1, а верхняя — от 2 до 1 (зеркально). Посмотрим на матожидание ошибки при различных параметрах генерирования:
Biased Log of Uniform Distribution
nulldata 100000
scalar step=0.001
#
loop for (bnd=0.001; bnd<1; bnd+=step) --progressive series eps = randgen(u,bnd,2-bnd) series leps = ln(eps) genr avgleps = mean(leps) store biasedln.gdt bnd avgleps endloop # open biasedln.gdt gnuplot avgleps bnd --with-lines --output=display --suppress-fitted

Уже лучше, но всё равно плохо. Если ошибка имеет равномерное распределение от \(b\) до \((2-b)\), то её среднее равно одному, однако среднее её логарифма нулю никак не равно. А вот чему оно равно:
\[\varepsilon \sim \mathcal{U}[b;2-b] \quad \Rightarrow \quad \mathbb{E}\bigl(\ln \varepsilon\bigr) = \frac{-2 b+b \ln (2-b)+b \ln b-2 \ln (2-b)+2}{2 (b-1)} \]
Интересная функция. Её предел в нуле равен \(\ln 2 - 1 \approx -0{,}307\), а при стремлении к единице её значение стремится к нулю. Однако это значит, что для минимизации смещения требуется брать значения границ, равные \([0{,}95;1{,}05]\) или ещё у́же, а до такого додумается далеко не каждый: все будут бояться, что при таком малом разбросе ошибок получатся гигантские (\(\geqslant50\)) t-статистики, а сама регрессия потеряет всякий смысл, исчезнет случайность, а «наблюдаемые» значения будут лежать почти в одной гиперплоскости. При этом все забывают, что только лишь очень близкие к единице \(\varepsilon\) позволяют задействовать отношение эквивалентности \(\ln (1 + \delta) \sim \delta\) при \(\delta\approx 0\).
Biased log of uniform distribution

 

Мораль 1. При генерировании искусственных данных помните, что в модели с мультипликативными ошибками очень желательна нормальность их логарифмов с центром массы в нуле. Если матожидание ошибок, преобразованных в аддитивные, не равно нулю, коэффициенты МНК будут в порядке, но оценка константы будет смещённой и несостоятельной, что критично при оценке технологического фактора в модели Кобба—Дугласа.

Мораль 2. Решаете задачку, гоняете циферки, а вам нужны ошибки, которые не сместят коэффициентов? Берите \(\varepsilon = \exp\bigl(\mathcal{N}(0;1)\bigr)\).

 

P.S. Мне кажется, что при симуляции Монте-Карло имеет смысл поиграться с функциональной формой ошибок и добиться того, чтобы проверка выполнения условия \(\sum(y_i^* - {\boldsymbol{x}^*_i}'\boldsymbol\beta)=0\) начиналась на шаг ранее, когда исходная форма \(y_i = f(\boldsymbol{x}_i, \boldsymbol\beta, \varepsilon_i)\) преобразуется в \(y_i^* = f^{-1}(y_i) = f^{-1}\bigl(f(\boldsymbol{x}_i, \boldsymbol\beta, \varepsilon_i)\bigr) = {\boldsymbol{x}^*_i}'\boldsymbol\beta + \varepsilon_i^*\), и чтобы особое внимание уделялось условию \(\mathbb{E}(\varepsilon_i^*)=0\). Конечно, это не так важно, но всё-таки есть такое эмпирическое правило, что наличие значимой константы говорит о пропущенных переменных; поэтому если в истинной модели есть константа, но систематическая ошибка её сильно снижает, то может возникнуть ложное ощущение того, что в модели учтено достаточное количество влияющих факторов, в то время как друг на друга накладываются большая необъяснённая остаточная вариация и систематическая ошибка.

About Andreï Kostyrka

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

6 Responses to Матожидание ошибок при симуляции МНК

  1. Егор says:

    Кажется, нарушение не самое страшное. МНК-оценки $\alpha$ и $\beta$ по крайней мере останутся несмещенными и состоятельными. В оценку мультипликативной константы подмешается матожидание экспоненты ошибки логарифмической модели, но так ли это важно экономисту?

    • А то, что нарушается условие первого порядка \(\sum \varepsilon_i=0\), это ничего? МНК не приходит в точку минимума суммы квадратов ошибок, коэффициенты в жизни не будут несмещёнными и состоятельными, причём даже асимптотически!
      Ах да, я забыл, это ж макроэкономисты любят строить межстрановые модели с мультипликативными факторами по 19 точкам и при наличии ошибок в измерении иксов, простите.

  2. Вольный дрыномашец says:

    Во-первых, условия первого порядка должно записывать в моментах, а не в суммах, то есть $E(\epsilon_i)=0$.
    Во-вторых, условие первого порядка меняется после перехода к логарифмам, и значит то условие, которое указали вы, выполняться не обязано.
    В-третьих, когда в модели есть константа, ошибка прекрасно нормализуется, какое бы матожидание она не имела. Если её матожидание 1, вычтите 1 и добавьте эту единицу к константе. Кстати говоря, когда константа вовсе не включена в модель, условие первого порядка, на которое вы ссылаетесь, вообще излишне (просто по построению). Это помогает понять, почему величина матожидания никак не влияет на состоятельность. А для асимптотики вроде асимптотических тестов сгодятся обычные оценки Уайта, нечего тут мудрить.
    В-четвертых, бравые макроэкономисты уже как пару десятилетий гоняют панельные регрессии с инструментами, надежность которых , в том числе в плане оценки производственной функции в разных странах, не вызывает каких-то (фундаментальных) сомнений.
    В-пятых, ошибак в иксах нет только в сказках.

  3. Егор says:

    Условие первого порядка выполнено априори, мы ведь получаем из него оценки, утверждая для МНК, что соответствующая производная равна нулю.

    Пример такой: посмотрим модель данных $Y = A L^\alpha \e$, $\e>0$. Ей соответствует модель в логарифмах $\ln Y = \ln A + \alpha \ln L + \ln \e$. МНК-оценка $\hat{\alpha} = \frac{\sum (\ln L_i – \bar{\ln L})\ln Y}{\sum (\ln L_i – \bar{\ln L})^2}$, подставляем туда $\ln Y$, получаем:
    $\hat{\alpha} = \frac{\sum (\ln L_i – \bar{\ln L})( \ln A + \alpha \ln L_i + \ln \e_i)}{\sum (\ln L_i – \bar{\ln L})^2} = \alpha + $\frac{\sum (\ln L_i – \bar{\ln L})\ln \e}{\sum (\ln L_i – \bar{\ln L})^2}$. Берем УМО по L, получаем, что вне зависимости от $\E(\ln \e|L)$ выполнено:
    $\E(\hat{\alpha}|L) = \alpha + \E(\ln \e|L) \frac{\sum (\ln L_i – \bar{\ln L})}{\sum (\ln L_i – \bar{\ln L})^2} = \alpha$. Какое бы ни было матожидание ошибки (кроме бесконечного), оценка наклона остается несмещенной. Состоятельность аналогично.

    Стоит подумать над тем, как располагая оценкой $\hat{\ln Y}$ из такой модели построить состоятельный прогноз $\hat{\Y}$ (как раз из того, что $\E(e^{\hat{\ln Y}}|L) < \E(Y|L)$).

    • Егор Дмитриевич,

      Прошу Вашего прощения. Я был глубоко неправ, полагая, что коэффициенты получатся смещёнными при нулевом матожидании эпсилон. Скоро напишу эррату. Пока что кратко покаюсь:

      — Коэффициенты будут состоятельными. Я только что проверил аналитически и симуляционно.

      — Единственное, что нельзя точно оценить в этой модели, — это мультипликатор технологии, \(a\). Всё остальное отлично получается, особенно когда число наблюдений достаточно велико (хотя бы тысяча… не как у Баумоля — 72 точки и ожесточённая полемика с де Лонгом, который нашёл у него ошибки в регрессорах).

      Забудьте последние два абзаца поста, их здесь никогда не было.

  4. Борис says:

    Я бы уточнил в строчке “Здесь avgleps…”, что математическое ожидание и среднее — условные, при положительных эпсилон. И заодно указал бы, какая строчка кода неявно выполняет отбор :)

Leave a Reply

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