Я больше года мучился, не понимал, почему аддитивность/мультипликативность нормальных ошибок играет такую большую роль, и даже разъяснение преподавателя не помогало. Почти год назад я написал этот пост: 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)} \]
Если хотите, чтобы матожидание логарифма было равно нулю, то надо брать хи-квадрат с \(1{,}866\,025\) степенями свободы.
Нормальное распределение может выдать абсолютно любые значения случайной величины, так как функция плотности определена на \(\mathbb{R}\). Поэтому если ошибка мультипликативна и нормальна, то, вообще говоря, пропадают все значения, где ошибка получилась меньше нуля. Кроме того, если ошибка распределена как \(\mathcal{N}(1;\sigma^2)\), а наивные пользователи думают, что её логарифм в среднем даст ноль, то спешу их разочаровать следующей картиной:
Здесь 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).
Чтобы избавиться от мнимой единицы, некоторые товарищи могут взять условное матожидание (\(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 (зеркально). Посмотрим на матожидание ошибки при различных параметрах генерирования:
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\).
Мораль 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\). Конечно, это не так важно, но всё-таки есть такое эмпирическое правило, что наличие значимой константы говорит о пропущенных переменных; поэтому если в истинной модели есть константа, но систематическая ошибка её сильно снижает, то может возникнуть ложное ощущение того, что в модели учтено достаточное количество влияющих факторов, в то время как друг на друга накладываются большая необъяснённая остаточная вариация и систематическая ошибка.
