Быстрый и грубый ответ на вопрос: сколько людей в России заболеет коронавирусной инфекцией COVID-19?

Оговорка. Данная публикация является выражением личного мнения автора, которое может совпадать или не совпадать с мнением экспертов ВОЗ, представителей государственных организаций, медицинских работников или других вовлечённых учреждений. Для получения последних обновлений о ситуации и профессиональных советов перейдите на официальный веб-сайт Всемирной организации здравоохранения .

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

В России часто слышится вопрос:

— Когда наша страна пройдёт пик?

Это сложный вопрос, поскольку под «пиком» обычно подразумевается «день с наибольшим количеством новых случаев», и на основе только лишь информации о пике нельзя принять никакие ответственные решения. Более общий вопрос обычно звучит так: «Как будет выглядеть динамика через 2 месяца?» При этом многие статьи — даже вроде бы научные! — говорят о какой-то эфемерной «логистической кривой». Тем не менее, у этой кривой одна проблема: вылезать труднее, чем влезать, поэтому динамика до и после пика может отличаться от зеркальной. Скажите, эти два объекта похожи?

Простая логистическая кривая и заболеваемость COVID 19 в Италии

Вот почему асимметрию надо моделировать в явном виде. Мы применим подход, который был использован в работе Fernández & Steel (1998), где авторы обобщают симметричные функции плотности распределения до асимметричных с одним параметром: левая часть растягивается, а правая — сжимается на один и тот же коэффициент, а затем функция нормализуется так, чтобы совпадали первые два момента исходного и нового распределений. Ради лёгкости интерпретации мы ничего не будем нормализовать: мы не моделируем плотность — мы просто берём логистическую кривую и смещаем, масштабируем и скашиваем её, чтобы подогнать под реальные данные. Кроме того, логистическая кривая имеет экспоненциально затухающие хвосты, что не слишком реалистично из-за того, что люди впадают в сладостное забвение, расслабляются, думают, что если в день 10–15 новых случаев, то всё уже закончилось, и начинают слишком близко общаться, и это-то и оборачивается катастрофой. Моделирование тяжести хвоста (скажем, с помощью t-распределения Стьюдента) может быть рассмотрено, однако следует опасаться малого числа наблюдений и вытекающую из этого ненадёжность оценки этого параметра. Формула для обобщённой логистической кривой, описывающей количество новых случаев в день, имеет следующий вид:
\[
g(x; \mu, \sigma, s, h) \mathrel{\stackrel{\text{def}}{=}}
\begin{cases}
h \cdot f \left( \frac{x - \mu}{\sigma \cdot s} \right), & \frac{x - \mu}{\sigma} > 0, \\
h \cdot f \left( \frac{s(x - \mu)}{\sigma} \right), & \frac{x - \mu}{\sigma} \le 0.
\end{cases}
\]
Здесь \(\mu\) — это мода (точка максимума функции, её пик), \(\sigma\) — горизонтальное растяжение, \(s\) — коэффициент скошенности (асимметрии), принимающий значения из диапазона \((0, +\infty)\) и показывающий относительную тяжесть правого хвоста (\(s > 1\) соответствует растяжению правого хвоста ещё правее), а \(h\) отражает вертикальное масштабирование (количество случаев меняется от страны к стране, поэтому абсолютная шкала должна тоже органически проистекать из данных). Помните, что f может быть любой симметричной функцией — мы используем функции, похожие на функции плотности, из соображений теории... и удобства. Нам нужен быстрый и простой ответ!

Обобщённая логистическая кривая для моделирования динамики коронавируса в России

Для начала скачаем мартовские и апрельские данные с Википедии со страниц вида ‘2020 coronavirus pandemic in X’, где X — название страны. Рассмотрим такие страны, как Испания, Германия, Италия, Люксембург, Швеция, Чехия, Турция и Россия. Собранные данные приводятся в таблице ниже; оговоримся, что это информация на 1 мая 2020 г. и что она может поменяться позже из-за обновления источников или изменения методологии.

Date spa ger ita lux swe cze tur rus
2020-03-01 85 117 1,694 0 7 3 0 0
2020-03-02 121 150 2,036 0 7 3 0 3
2020-03-03 166 188 2,502 0 11 5 0 3
2020-03-04 228 240 3,089 0 16 5 0 3
2020-03-05 282 349 3,858 2 17 8 0 4
2020-03-06 401 534 4,636 3 22 19 0 10
2020-03-07 525 684 5,883 4 29 26 0 14
2020-03-08 674 847 7,375 5 38 32 0 17
2020-03-09 1,231 1,112 9,172 5 54 38 0 20
2020-03-10 1,695 1,460 10,149 7 108 63 0 20
2020-03-11 2,277 1,884 12,462 7 181 94 1 28
2020-03-12 3,146 2,369 15,113 26 238 116 1 34
2020-03-13 5,232 3,062 17,660 34 266 141 5 45
2020-03-14 6,391 3,795 21,157 51 292 189 6 59
2020-03-15 7,988 4,838 24,747 77 339 298 18 63
2020-03-16 9,942 6,012 27,980 81 384 383 47 93
2020-03-17 11,826 7,156 31,506 140 425 450 98 114
2020-03-18 14,769 8,198 35,713 203 511 560 191 147
2020-03-19 18,077 10,999 41,035 335 603 765 359 199
2020-03-20 21,571 13,957 47,021 484 685 889 670 253
2020-03-21 25,496 16,662 53,578 670 767 1,047 947 306
2020-03-22 29,909 18,610 59,138 798 845 1,161 1,236 367
2020-03-23 35,480 22,672 63,927 875 927 1,287 1,529 438
2020-03-24 42,058 27,436 69,176 1,099 1,027 1,472 1,872 495
2020-03-25 50,105 31,554 74,386 1,333 1,125 1,763 2,433 658
2020-03-26 57,786 36,508 80,539 1,453 1,228 2,022 3,629 840
2020-03-27 65,719 42,288 86,498 1,605 1,317 2,395 5,698 1,036
2020-03-28 73,232 48,582 92,472 1,831 1,392 2,657 7,402 1,264
2020-03-29 80,110 52,547 97,689 1,950 1,450 2,817 9,217 1,534
2020-03-30 87,956 57,298 101,739 1,988 1,566 3,001 10,827 1,836
2020-03-31 95,923 61,913 105,792 2,178 1,683 3,308 13,531 2,337
2020-04-01 104,118 67,366 110,574 2,319 1,828 3,589 15,679 2,777
2020-04-02 112,065 73,522 115,242 2,487 1,982 3,858 18,135 3,548
2020-04-03 119,199 79,696 119,827 2,612 2,179 4,190 20,921 4,149
2020-04-04 126,168 85,778 124,632 2,729 2,253 4,472 23,934 4,731
2020-04-05 131,646 91,714 128,948 2,804 2,359 4,587 27,069 5,389
2020-04-06 136,675 95,391 132,547 2,843 2,566 4,822 30,217 6,343
2020-04-07 141,942 99,225 135,586 2,970 2,763 5,017 34,109 7,497
2020-04-08 148,220 103,228 139,422 3,034 2,890 5,312 38,226 8,672
2020-04-09 153,222 108,202 143,626 3,115 3,036 5,569 42,282 10,131
2020-04-10 158,273 113,525 147,577 3,223 3,100 5,732 47,029 11,917
2020-04-11 163,027 117,658 152,271 3,270 3,175 5,902 52,167 13,584
2020-04-12 166,831 120,479 156,363 3,281 3,254 5,991 56,956 15,770
2020-04-13 170,099 123,016 159,516 3,292 3,369 6,059 61,049 18,328
2020-04-14 174,060 125,098 162,488 3,307 3,539 6,141 65,111 21,102
2020-04-15 180,659 127,584 165,155 3,373 3,665 6,301 69,392 24,490
2020-04-16 184,948 130,450 168,941 3,444 3,785 6,433 74,193 27,938
2020-04-17 190,839 133,830 172,434 3,480 3,909 6,549 78,546 32,008
2020-04-18 194,416 137,439 175,925 3,537 3,964 6,654 82,329 36,793
2020-04-19 198,674 139,897 178,972 3,550 4,027 6,746 86,306 42,853
2020-04-20 200,210 141,672 181,228 3,558 4,173 6,900 90,980 47,121
2020-04-21 204,178 143,457 183,957 3,618 4,284 7,033 95,591 52,763
2020-04-22 208,389 145,694 187,327 3,654 4,406 7,132 98,674 57,999
2020-04-23 213,024 148,046 189,973 3,665 4,515 7,187 101,790 62,773
2020-04-24 219,764 150,383 192,994 3,695 4,610 7,273 104,912 68,622
2020-04-25 223,759 152,438 195,351 3,711 4,686 7,352 107,773 74,588
2020-04-26 226,629 154,175 197,675 3,723 4,756 7,404 110,130 80,949
2020-04-27 229,422 155,193 199,414 3,729 4,884 7,445 112,261 87,147
2020-04-28 232,128 156,337 201,505 3,741 4,980 7,504 114,653 93,558
2020-04-29 236,899 157,641 203,591 3,769 5,040 7,579 117,589 99,399
2020-04-30 239,369 159,119 205,463 3,784 5,051 7,682 120,204 106,498

Теперь мы должны выбрать метод оценивания. В принципе, для каждой страны у нас есть пары наблюдений \((1, Y_1), (2, Y_2), \ldots, (61, Y_{61})\), и предполагается, что функциональная форма известна с точностью до 4 (логистическая) или 5 (Стьюдент) параметров. Самым очевидным способом получения оценки является минимизация некоторой штрафной статистики, основанной на функции от расхождений между реальными данными и предсказаниями модели; самый простой случай — это сумма (или среднее) квадратов остатков. Однако на практике часто используются более робастные штрафные функции, а иногда и более робастные статистики, поэтому мы также попробуем для получения оценки параметров проминимизировать по оным сумму модулей (абсолютных отклонений, как в нелинейной квантильной регрессии), а также медиану квадратов ошибок (как в регрессии наименьшей медианы квадратов авторства Rousseeuw). Мы будем минимизировать эти 3 критерия. Обратите внимание, что последний из них не является гладким по параметрам, поэтому может существовать целый континуум значений параметров, дающий одно и то же значение критерия (даже глобальный минимум).
\[
C_1 = \frac{1}{n} \sum_{t=1}^T (Y_t - g(t; \theta))^2, \quad C_2 = \frac{1}{n} \sum_{t=1}^T |Y_t - g(t; \theta)|, \quad C_3 = \mathrm{median}[(Y_t - g(t; \theta))^2]
\]
В-третьих, поскольку мы моделируем параметрические кривые с четырьмя параметрами, нужно прикинуть приблизительные оценки, чтобы близко угадать динамику. Все вышеописанные критерии нелинейны по параметрам, поэтому мы вручную подберём реалистичные значения для каждой из стран, а затем размахнёмся от них вверх-вниз, сделаем просторный гиперкуб и попробуем с помощью стохастического оптимизатора найти в нём точку-кандидата для оптимума. После этого мы доведём его в точный максимум с помощью оптимизатора Нельдера—Мида (с добавлением разумных ограничений).

Запрограммируем общую кривую с помощью нашей собственной функции, которая будет принимать на вход независимую переменную времени (количество дней, прошедших с 1 марта), тип распределения (логистическое или Стьюдента, но при наличии времени и желания можно прописать любые) и параметры этого распределения (вектор длины 4 или 5: пик, горизонтальный масштаб, скошенность, вертикальный масштаб, а для распределения Стьюдента — количество степеней свободы).

Отметим, что это не функция плотности и что не нужно думать, что для неё будут выполняться свойства логистического или t-распределения, так как она имеет другой масштаб, и вообще, это просто кривая, которая соответствует количеству людей, то есть случаев заболевания. Так выглядит динамика для всех стран на одном графике.

Новые случаи заражения коронавирусов как доля от максимума за всю историю

Теперь выберем хорошие стартовые значения для оптимизатора. Для численной стабильности (особенно в градиентных методах) желательно, чтобы у всех параметров был один и тот же порядок величины (можно всё перемасштабировать), и параметр вертикального масштаба, который измеряет людей, гуляет где-то в тысячах или десятках тысяч, поэтому поделим его на 10 000. Для каждой страны мы прикинули на глаз параметры \((\tilde \mu, \tilde\sigma, \tilde s, \tilde h\ [, \widetilde{df}])\), а для стохастического оптимизатора будем использовать гиперкуб, пропорциональный нашим прикидкам: от \((\tilde\mu/2, \tilde\sigma/3, \tilde s/3, \tilde h / 4\ [, 0.1])\) до \((1.5\tilde\mu, 3\tilde\sigma, 3\tilde s, 3\tilde h\ [, 50])\).

Начальные параметры для оценивания динамики коронавируса

Для оптимизации мы используем пакет hydroPSO, так как в нём всё великолепно параллелизуется и гибко настраивается. Пространства высокой размерности становятся более разреженными быстрее, чем кажется, поэтому для плотного из покрытия мы используем начальную популяцию точек размером 400 или 600. Сходимость довольно быстрая для логистической модели и более медленная для модели Стьюдента, поэтому мы используем 50–100 итераций для первой и 75–125 итераций для второй (это всё же грубая предварительная оценка первого шага). Одна итерация занимает 0,25–0,5 секунды, поэтому оценка всех моделей заняла 8 ⋅ (50 + 75 + 100 + 75 + 100 + 125) ⋅ 0,35 ≈ 1500 секунд, или примерно полчаса, — такова цена снижения риска скатывания в локальный оптимум.

Вы думаете, что всё прямо так сразу можно оценить?.. Не тут-то было!

Нам мешает одно яблочко... то есть облачко на горизонте: пока что в России не было пика! Это создаёт ворох проблем: строго говоря, параметр асимметрии и параметр масштаба вместе не идентифицированы, поскольку мы наблюдаем только левую часть распределения: можно умножить параметр асимметрии на 2 и параметр масштабирования на 2 и получить абсолютно такой же левый хвост. В этом случае опасно надеяться, что в самом правом конце данных содержится пик. Известный факт, что оценивание параметров усечённого нормального распределения, когда среднее лежит за пределами диапазона усечения, очень ненадёжно; в нашей ситуации пик также находится за пределами диапазона, что усугубляет положение. Это означает, что получения оценок для России параметр асимметрии должен быть зафиксирован (как параметр других стран), поэтому будет оцениваться только горизонтальный масштаб. Например, если бы у России была динамика затухания, как в Люксембурге, то каким было бы общее число случаев? Иными словами, мы можем делать сценарные прогнозы: турецкий сценарий, итальянский сценарий, среднеевропейский сценарий.

Максимизация второго шага из возможного решения считается за мгновение, и оценки не сильно меняются. Однако результаты оказались более нестабильными для Швеции, одна модель дала странные результаты для Турции, а в случае в Россией вообще всё поехало — потому что мы ещё не применили трюк с фиксированием параметра.

Оценённая логистическая кривая количества случаев заражения коронавирусом

Видим, что среди оценённых моделей нет явного лидера, так как и логистическая, так и стьюдентовская подгонка дают примерно одинаковые результаты. Воспользуемся принципом бритвы Оккама и во избежание перепараметризации ограничим наш анализ только логистическим распределеним и забудем про оценку тяжести хвостов.

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

Вот и настал момент, которого все ждали с нетерпением: сколько же случаев заболевания коронавирусной инфекцией ожидается в России в 2020 году? Когда настанет тот день, когда будет ожидаться менее одного нового заболевшего? Ещё лучше, когда ждать менее чем одного случая в течение двух недель (то есть 1/14 случаев в день), и каково будет количество случаев до этой даты — когда закончится пандемия? Для этого мы можем вызвать вычислить очень простую вещь: идти вправо во времени до тех пор, пока наша функция прироста случаев с оценёнными параметрами не примет значение 1 или 1/14. Общее ожидаемое количество инфицированных будет кумулятивной суммой вплоть до этой даты.

Помните, мы рассматриваем только логистические модели для прогнозирования количества случаев в каждой стране до дня, когда будет менее 1 случая или менее 1/14 случаев в день. Для каждой страны мы выбираем в качестве лучшей модели (из трёх логистических с разными оценивателями) ту, которая точнее всего предсказывает общее количество случаев по состоянию на последнюю дату (то есть на конец выборки, 30 апреля).

Случаев к 30 апреля Сред. наименьш. квадр. Сред. наименьш. модулей Медиан. наименьш. квадратов
Испания 239,369 237,304 242,555 203,910
Германия 159,119 158,465 164,505 183,371
Италия 205,463 203,246 204,339 212,175
Люксембург 3,784 3,749 4,108 5,147
Швеция 5,051 5,076 5,046 5,160
Чехия 7,682 7,646 7,710 8,033
Турция 120,204 121,944 116,104 102,600
Итого к 30 апреля Дней с 1 марта до 1 чел./д. Дней с 1 марта до 1 чел./14 д. Случаев до 1 чел./д. Случаев до 1 чел./14 д.
Испания 239,369 189 231 282,401 282,415
Германия 159,119 142 172 171,369 171,379
Италия 205,463 201 249 241,977 241,993
Люксембург 3,784 83 109 3,838 3,847
Швеция 5,051 148 196 6,752 6,768
Чехия 7,682 95 120 8,029 8,037
Турция 120,204 140 166 145,922 145,930
Пик (μ) Гориз. масшт. (σ) Скошенность (s) Верт. масшт. (h)
Испания 25.383 7.859 2.026 2.854
Германия 27.533 7.050 1.612 2.180
Италия 20.519 9.619 1.881 2.103
Люксембург 20.131 3.999 2.410 0.068
Швеция 32.858 12.033 1.519 0.052
Чехия 27.818 6.764 1.400 0.113
Турция 41.888 7.754 1.277 1.829

Но постойте! Забыли про Россию!

У нас 7 стран, значит, 7 разных сценариев. Эти сценарии в терминах параметров описываются днём ​​пика (\(\mu\)), горизонтальным размахом (\(\sigma\)), асимметрией (\(s\)), вертикальным масштабом (\(h\)), а для России мы не можем знать степень асимметрии. Сравните эти два графика.

Асимметричная логистическая кривая, дающая различные хвосты

Для одного и того же левого хвоста мы получаем разные правые хвосты: в примере выше красный в 4 раза тяжелее синего. Это очень сильно сказывается на прогнозах: синяя кривая соответствует пандемии, которая затухает до 1 нового случая в день за 66 дней и подкашивает за всё время 37 469 человек, а красная кривая — пандемии, которая тянется 232 дня и охватывает 127 451 случай! Это означает, что начальные значения для России должны быть изменены в соответствии со всеми семью сценариями для других стран.

После переоценки логистической модели для России с асимметрией, зафиксированной на значении 1 (минимизация среднего модуля отклонений дала немного более точное общее число случаев с относительной погрешностью чуть менее 1%), мы можем умножить и \(\sigma\), и \(s\) на оценку параметра \(s\) для других стран, рассчитать два красных дня календаря (менее 1 случая за 1 день или за 14 дней) и получить общее количество случаев.

Оценки для России: \(\hat\mu = 56.1\), \(\hat \sigma / s = 6.33\), \(\hat h = 2.49\). Это означает, что, несмотря на отсутствие чёткого пика на 61-й день, он таки имел быть на 56-й день, непосредственно перед окончанием сбора данных.

Дней до 1 чел./д. Дней до 1 чел./14 д. Случаев до 1 чел./д. Случаев до 1 чел./14 д.
Сценарий Испании 320 388 402,422 402,445
Сценарий Германии 223 267 283,730 283,745
Сценарий Италии 284 343 357,887 357,907
Сценарий Люксембурга 429 526 536,540 536,573
Сценарий Швеции 205 243 260,712 260,724
Сценарий Чехии 182 215 233,328 233,339
Сценарий Турции 161 188 207,443 207,452
Средний сценарий (s = 1.73) 249 299 315,371 315,387

Прогноз нового количеств случаев заражения коронавирусом в России

Прогноз общего числа случаев заражения COVID 19 в России

Прежде чем перейти к заключению, обязательно оговорим несколько «но»:

  1. Оценка смещения кривой (\(\hat\mu\)) по левой ея части очень ненадёжна, поэтому вокруг этой оценки существует огромный доверительный интервал. Насколько огромный? Любопытствующий читатель может даже рассчитать стандартные ошибки, но о нормальности оценки не может быть и речи, поскольку хорошо известно, что поведение оценок на границе области допустимых параметров далеко от асимптотически нормального (как, например, метод максимального правдоподобия не работает для оценки экстремальных порядковых статистик).
  2. Сам пик, который мы определили как максимальное число случаев, соответствует пику сглаженной кривой, даже если отдельные дневные наблюдения могут скакать выше или ниже — предполагается, помимо всего прочего, что в среднем ошибка ноль. Если существует систематическое занижение отчётности из-за недостаточного количества доступных тестов, а реальные цифры выше, то эта упрощённая модель не может дать ничего полезного исходя из имеющихся данных, и требуется использовать другие, более сложные методы.
  3. Множество бессимптомных случаев, которые мало где регистрируются, составляет подводную часть айсберга, поэтому, если этот невидимый статистике айсберг полностью поднять над уровнем моря, то вершина сместится вправо.
  4. Дневные цифры сильно зависят от тактики тестирования, хорошей погоды (это не «экономика солнечных пятен» Джевонса, а корреляция с количеством людей, собирающихся на улице), выступлений президента и иных неучтённых переменных.

Все эти оговорки свидетельствуют о том, что полученная нами оценка, скорее всего, недооценивает положение реального пика, но для решения проблем, изложенных выше, нужно обзавестись более качественными пространственными данными, собрать больше информации о плотности населения, географическом распределении — настоящая работа для настоящих эпидемиологов.

Выводы — если Россия будет повторять динамику некоторых европейских стран и если забыть вышеизложенные оговорки — таковы:

  1. Пик, скорее всего, уже наступил — он пришёлся на самый конец апреля 2020 года;
  2. Общее количество заражённых COVID-19 россиян составит от 200 до 500 тысяч;
  3. При среднем европейском сценарии и средней асимметрии динамики (при которой уменьшение количества случаев в день после пика примерно в 3 раза медленнее, чем нарастание до пика) общее число случаев в России за всё время составит примерно 315 тысяч;
  4. Ожидается, что в начале ноября 2020 г. число новых случаев в день станет меньше 1, а к концу декабря будет меньше 1 случая в течение 2 недель.

Только время покажет, насколько точны эти примитивные оценки, поэтому рекомендуется повторить оценивание на новых данных (возможно, на большем количестве стран) и перестроить прогнозы.

Скопируйте приведённую выше таблицу в R в объект ds и используйте приведённый ниже код для репликации результатов из данной статьи.

library(stargazer)
library(hydroPSO)
 
cs <- c("spa", "ger", "ita", "lux", "swe", "cze", "tur", "rus")
csn <- c("Spain", "Germany", "Italy", "Luxembourg", "Sweden", "Czech Republic", "Turkey", "Russia")
 
mycurve <- function(x, distrib, pars) {
  pars[4] <- pars[4] * 1e4
  xstd <- (x - pars[1]) / pars[2]
  f <- switch(distrib,
              logistic = ifelse(xstd > 0, dlogis(xstd/pars[3]), dlogis(xstd*pars[3])) * pars[4],
              student = ifelse(xstd > 0, dt(xstd/pars[3], df = pars[5]), dt(xstd*pars[3], df = pars[5])) * pars[4]
  )
  return(f)
}
 
png("plot1-logistic-curve-coronavirus-italy.png", 640, 400, type = "cairo", pointsize = 16)
plot(diff(ds$ita), type = "l", xlab = "Days since the 1st of March", ylab = "New cases per day", bty = "n", lwd = 2, main = "Logistic curve for Italy")
lines(1:60, mycurve(1:60, "logistic", c(25, 6, 1, 2.45)), col = "red", lty = 2, lwd = 2)
text(50, 1500, "WTF?!", col = "red", cex = 2)
dev.off()
 
mycols <- c("#000000", "#008800", "#AA0000")
 
png("plot2-density-logistic-curve-for-corona-virus-flattening.png", 720, 360, pointsize = 16)
par(mar = c(2, 2, 3, 1), mfrow = c(1, 2))
curve(dlogis(x), -6, 6, n = 1001, main = "Skewed logistic density", bty = "n", lwd = 2, ylim = c(0, dlogis(0)), ylab = "", xlab = "")
curve(mycurve(x, "logistic", c(0, 1, 1.5, 1e-4)), -6, 6, n = 1001, add = TRUE, col = mycols[2], lwd = 2, lty = 2)
curve(mycurve(x, "logistic", c(0, 1, 0.5, 1e-4)), -6, 6, n = 1001, add = TRUE, col = mycols[3], lwd = 2, lty = 2)
legend("topright", c("log.", "log., skew = 1.5", "log., skew = 0.5"), col = mycols, pch = 16, bty = "n")
 
curve(dt(x, 5), -6, 6, n = 1001, main = "Skewed Student’s density", bty = "n", lwd = 2, ylim = c(0, dt(0, 5)))
curve(mycurve(x, "student", c(0, 1, 1.5, 1e-4, 5)), -6, 6, n = 1001, add = TRUE, col = mycols[2], lwd = 2, lty = 2)
curve(mycurve(x, "student", c(0, 1, 0.5, 1e-4, 5)), -6, 6, n = 1001, add = TRUE, col = mycols[3], lwd = 2, lty = 2)
legend("topright", c("t(5)", "t(5, skew = 1.5)", "t(5, skew = 0.5)"), col = mycols, pch = 16, bty = "n")
dev.off()
 
ds[, (length(cs)+2):(length(cs)*2+1)] <- sapply(ds[, 2:(length(cs)+1)], function(x) diff(c(NA, x)))
colnames(ds)[(length(cs)+2):(length(cs)*2+1)] <- paste0("d", cs)
ds <- ds[2:nrow(ds), ]
 
mycols <- rainbow(length(cs), end = 0.8, v = 0.8)
png("plot3-new-covid-19-cases-as-a-fraction-of-maximum.png", 640, 400, pointsize = 16)
plot(NULL, NULL, xlim = c(1, nrow(ds)), ylim = c(0, 1), xlab = "Days since March, 1", ylab = "Cases / Cases at peak", bty = "n", main = "Scaled dynamics of new cases per day")
for (i in (length(cs)+2):(length(cs)*2+1)) lines(ds[, i] / max(ds[, i]), type = "l", lwd = 2, col = mycols[i - length(cs) - 1])
legend("topleft", cs, pch = 16, col = mycols, bty = "n")
dev.off()
 
start.pars <- list(list(c(27, 8, 1.8, 3), c(26, 11, 2, 2, 5)), # Spain
                   list(c(26, 6, 2, 2.3), c(27, 9, 2, 1.5, 5)), # Germany
                   list(c(20, 9, 2.1, 2.2), c(20, 12, 2.1, 1.5, 5)), # Italy
                   list(c(25, 5, 1.5, 0.08), c(25, 7, 1.5, 0.055, 5)), # Luxembourg
                   list(c(35, 8, 1.2, 0.07), c(35, 12, 1.2, 0.045, 5)), # Sweden
                   list(c(30, 7, 1.3, 0.11), c(30, 11, 1.3, 0.07, 5)), # Czech,
                   list(c(43, 7, 1.3, 1.8), c(43, 11, 1.3, 1.2, 5)), # Turkey
                   list(c(60, 7, 1, 2.8), c(58, 10, 1, 1.8, 5)) # Russia
                   )
start.pars <- lapply(start.pars, function(x) {names(x) <- c("logistic", "student"); x})
 
leg <- function(i) {
  c(parse(text = paste0('mu == ', start.pars[[i]][[1]][1], '*" | "*', start.pars[[i]][[2]][1])),
    parse(text = paste0('sigma == ', start.pars[[i]][[1]][2], '*" | "*', start.pars[[i]][[2]][2])),
    parse(text = paste0('s == ', start.pars[[i]][[1]][3], '*" | "*', start.pars[[i]][[2]][3])),
    parse(text = paste0('h == ', start.pars[[i]][[1]][4], '*" | "*', start.pars[[i]][[2]][4])),
    parse(text = paste0('df(t) == ', start.pars[[i]][[2]][5])))
}
 
png("plot4-initial-values-for-search.png", 720, 1200, pointsize = 16)
par(mfrow = c(4, 2))
for (i in 1:length(cs)) {
  plot(ds[, paste0("d", cs[i])], type = "l", lwd = 2, xlab = "Days since March, 1", ylab = "Registered cases", bty = "n", main = paste0(csn[i], " initial values"), ylim = c(0, max(ds[, paste0("d", cs[i])])))
  lines(mycurve(1:nrow(ds), "logistic", start.pars[[i]][[1]]), col = mycols[i], lwd = 2, lty = 2)
  lines(mycurve(1:nrow(ds), "student", start.pars[[i]][[2]]), col = mycols[i], lwd = 2, lty = 3)
  legend("topleft", legend = leg(i), bty = "n")
}
dev.off()
 
criterion <- function(pars, observed, distrib, fun = function(x) mean(x^2)) {
  theor <- mycurve(x = 1:length(observed), distrib = distrib, pars = pars)
  return(fun(theor - observed))
}
 
rls <- function(x) sqrt(mean(x^2)) # Root mean least square; root is here for the purposes of comparison
nla <- function(x) mean(abs(x))
nms <- function(x) sqrt(median(x^2)) # Root median of squares
opt.res <- vector("list", length(cs))
for (i in 1:length(cs)) {
  sl <- start.pars[[i]][[1]]
  st <- start.pars[[i]][[2]]
  ll <- start.pars[[i]][[1]] / c(2, 3, 3, 4) # Lower bound for logistic
  ul <- start.pars[[i]][[1]] * c(1.5, 3, 3, 3) # Upper bound for logistic
  lt <- c(start.pars[[i]][[2]][1:4] / c(2, 3, 3, 4), 0.1) # Lower bound for Student
  ut <- c(start.pars[[i]][[2]][1:4] * c(1.5, 3, 3, 3), 50) # Upper bound for Student
  contl1 <- list(npart = 400, REPORT = 1, parallel = "parallel", par.nnodes = 6, maxit = 50, plot = TRUE)
  contl2 <- list(npart = 400, REPORT = 1, parallel = "parallel", par.nnodes = 6, maxit = 75, plot = TRUE)
  contl3 <- list(npart = 400, REPORT = 1, parallel = "parallel", par.nnodes = 6, maxit = 100, plot = TRUE)
  contt1 <- list(npart = 600, REPORT = 1, parallel = "parallel", par.nnodes = 6, maxit = 75, plot = TRUE)
  contt2 <- list(npart = 600, REPORT = 1, parallel = "parallel", par.nnodes = 6, maxit = 100, plot = TRUE)
  contt3 <- list(npart = 600, REPORT = 1, parallel = "parallel", par.nnodes = 6, maxit = 125, plot = TRUE)
  set.seed(1); cat(csn[i], "logistic: least squares\n")
  opt.log.nls <- hydroPSO(par = sl, fn = function(x) criterion(x, ds[, paste0("d", cs[i])], "logistic", rls), lower = ll, upper = ul, control = contl1)
  set.seed(1); cat(csn[i], "logistic: least absolute deviations\n")
  opt.log.nla <- hydroPSO(par = sl, fn = function(x) criterion(x, ds[, paste0("d", cs[i])], "logistic", nla), lower = ll, upper = ul, control = contl2)
  set.seed(1); cat(csn[i], "logistic: median of least squares\n")
  opt.log.nms <- hydroPSO(par = sl, fn = function(x) criterion(x, ds[, paste0("d", cs[i])], "logistic", nms), lower = ll, upper = ul, control = contl3)
  set.seed(1); cat(csn[i], "student: least squares\n")
  opt.t.nls <- hydroPSO(par = st, fn = function(x) criterion(x, ds[, paste0("d", cs[i])], "student", rls), lower = lt, upper = ut, control = contt1)
  set.seed(1); cat(csn[i], "student: least absolute deviations\n")
  opt.t.nla <- hydroPSO(par = st, fn = function(x) criterion(x, ds[, paste0("d", cs[i])], "student", nla), lower = lt, upper = ut, control = contt2)
  set.seed(1); cat(csn[i], "student: median of least squares\n")
  opt.t.nms <- hydroPSO(par = st, fn = function(x) criterion(x, ds[, paste0("d", cs[i])], "student", nms), lower = lt, upper = ut, control = contt3)
 
  opt.res[[i]] <- list(opt.log.nls, opt.log.nla, opt.log.nms, opt.t.nls, opt.t.nla, opt.t.nms)
  save(opt.res, file = "opt.RData", compress = "xz")
}
 
ui4 <- diag(4)
ci4 <- c(10, 1, 0.01, 1e-5)
ui5 <- diag(5)
ci5 <- c(10, 1, 0.01, 1e-5, 0.1)
 
opt.refined <- vector("list", length(cs))
for (i in 1:length(cs)) {
  cat(csn[i], "\n")
  sl1 <- opt.res[[i]][[1]]$par # Replace with start.pars[[i]]$logistic if you did not do the stochastic optimisation
  sl2 <- opt.res[[i]][[2]]$par # Replace with start.pars[[i]]$logistic if you did not do the stochastic optimisation
  sl3 <- opt.res[[i]][[3]]$par # Replace with start.pars[[i]]$logistic if you did not do the stochastic optimisation
  st1 <- opt.res[[i]][[4]]$par # Replace with start.pars[[i]]$student if you did not do the stochastic optimisation
  st2 <- opt.res[[i]][[5]]$par # Replace with start.pars[[i]]$student if you did not do the stochastic optimisation
  st3 <- opt.res[[i]][[6]]$par # Replace with start.pars[[i]]$student if you did not do the stochastic optimisation
 
  opt.log.nls <- constrOptim(sl1, function(x) criterion(x, ds[, paste0("d", cs[i])], "logistic", rls), grad = NULL, ui = ui4, ci = ci4, outer.iterations = 300)
  opt.log.nla <- constrOptim(sl2, function(x) criterion(x, ds[, paste0("d", cs[i])], "logistic", nla), grad = NULL, ui = ui4, ci = ci4, outer.iterations = 300)
  opt.log.nms <- constrOptim(sl3, function(x) criterion(x, ds[, paste0("d", cs[i])], "logistic", nms), grad = NULL, ui = ui4, ci = ci4, outer.iterations = 300)
  opt.t.nls <- constrOptim(st1, function(x) criterion(x, ds[, paste0("d", cs[i])], "student", rls), grad = NULL, ui = ui5, ci = ci5, outer.iterations = 300)
  opt.t.nla <- constrOptim(st2, function(x) criterion(x, ds[, paste0("d", cs[i])], "student", nla), grad = NULL, ui = ui5, ci = ci5, outer.iterations = 300)
  opt.t.nms <- constrOptim(st2, function(x) criterion(x, ds[, paste0("d", cs[i])], "student", nms), grad = NULL, ui = ui5, ci = ci5, outer.iterations = 300)
 
  opt.refined[[i]] <- list(opt.log.nls, opt.log.nla, opt.log.nms, opt.t.nls, opt.t.nla, opt.t.nms)
}
save(opt.res, opt.refined, file = "opt.RData", compress = "xz")
 
 
xmax <- 120
xinf <- 1000
 
mycols <- c("#CC0000", "#0000FFCC")
png("plot5-fitted-logistic-curve-total-infected-coronavirus-cases.png", 720, 1200, type = "cairo", pointsize = 16)
par(mfrow = c(4, 2))
for (i in 1:length(cs)) {
  plot(ds[, paste0("d", cs[i])], xlim = c(1, xmax), type = "l", lwd = 3, xlab = "Days since March, 1", ylab = "Cases", bty = "n", main = paste0("Fitted curves for ", csn[i])) 
  for (j in 1:3) {
    lines(mycurve(1:xmax, distrib = "logistic", pars = opt.refined[[i]][[j]]$par), lwd = 2, lty = j, col = mycols[1])
    lines(mycurve(1:xmax, distrib = "student", pars = opt.refined[[i]][[j + 3]]$par), lwd = 2, lty = j, col = mycols[2])
  }
  legend("topleft", c("Logis.", "Stud."), pch = 16, col = mycols, bty = "n")
  legend("topright", c("Mean(least squares)", "Mean(least abs. dev.)", "Median(least abs. sq.)"), lty = 1:3, lwd = 2, bty = "n")
}
dev.off()
 
res.tab <- vector("list", length(cs))
for (i in 1:length(cs)) {
  cat(csn[i])
  a <- cbind(c(start.pars[[i]]$logistic, NA), c(opt.refined[[i]][[1]]$par, NA), c(opt.refined[[i]][[2]]$par, NA), c(opt.refined[[i]][[3]]$par, NA),
             start.pars[[i]]$student, opt.refined[[i]][[4]]$par, opt.refined[[i]][[5]]$par, opt.refined[[i]][[6]]$par)
  colnames(a) <- c("Log. init. val.", "Log. LS", "Log. LAD", "Log. LMS", "Stud. init. val.", "Stud. LS", "Stud. LAD", "Stud. LMS")
  row.names(a) <- c("Peak", "Horiz. scale", "Skew", "Vert. scale", "Stud. DoF")
  res.tab[[i]] <- a
  # stargazer(a, out = paste0(csn[i], ".html"), type = "html", summary = FALSE, digits = 2)
}
 
find.day <- function(cases, distrib, pars) {
  uniroot(function(x) mycurve(x, distrib, pars) - cases, lower = pars[1], upper = pars[1] + 180, extendInt = "downX", maxiter = 200, tol = 1e-8)$root
}
 
get.cases <- function(xmax, distrib, pars) {
  round(sum(mycurve(1:xmax, distrib, pars)))
}
 
days.to.one <- days.to.tw <- cases.one <- cases.tw <- cases.insample <- matrix(NA, nrow = length(cs), ncol = 6)
best.mod <- best.log <- numeric(length(cs))
for (i in 1:8) {
  days.to.one[i, ] <- c(sapply(1:3, function(j) find.day(1, "logistic", opt.refined[[i]][[j]]$par)),
                        sapply(4:6, function(j) find.day(1, "student", opt.refined[[i]][[j]]$par)))
  days.to.tw[i, ] <- c(sapply(1:3, function(j) find.day(1/14, "logistic", opt.refined[[i]][[j]]$par)),
                        sapply(4:6, function(j) find.day(1/14, "student", opt.refined[[i]][[j]]$par)))
  cases.one[i, ] <- c(sapply(1:3, function(j) get.cases(ceiling(days.to.one[i, j]), "logistic", opt.refined[[i]][[j]]$par)),
                      sapply(4:6, function(j) get.cases(ceiling(days.to.one[i, j]), "student", opt.refined[[i]][[j]]$par)))
  cases.tw[i, ] <- c(sapply(1:3, function(j) get.cases(ceiling(days.to.tw[i, j]), "logistic", opt.refined[[i]][[j]]$par)),
                      sapply(4:6, function(j) get.cases(ceiling(days.to.tw[i, j]), "student", opt.refined[[i]][[j]]$par)))
  cases.insample[i, ] <- c(sapply(1:3, function(j) get.cases(nrow(ds), "logistic", opt.refined[[i]][[j]]$par)),
                           sapply(4:6, function(j) get.cases(nrow(ds), "student", opt.refined[[i]][[j]]$par)))
  best.mod[i] <- which.min(abs(ds[nrow(ds), cs[i]] - cases.insample[i, ]))
  best.log[i] <- which.min(abs(ds[nrow(ds), cs[i]] - cases.insample[i, 1:3]))
  cat(csn[i], "has the best in-sample fit from model", best.mod[i], "and among logistic models,", best.log[i], "\n")
}
colnames(days.to.one) <- colnames(days.to.tw) <- c("Log. LS", "Log. LAD", "Log. LMS", "Stud. LS", "Stud. LAD", "Stud. LMS")
row.names(days.to.one) <- row.names(days.to.tw) <- csn
 
mod.inds <- c("Mean least squares", "Mean abs. dev.", "Median least sq.")
cases.ins <- cbind(as.numeric(ds[nrow(ds), 2:(length(cs))]), cases.insample[1:7, 1:3])
colnames(cases.ins) <- c("Real cases by April, 30", paste0(mod.inds, " forecast"))
row.names(cases.ins) <- csn[1:7]
stargazer(cases.ins, type = "html", out = "insample.html", summary = FALSE)
best.log
 
best.pars <- lapply(1:7, function(i) opt.refined[[i]][[best.log[i]]]$par)
best.pars <- do.call(rbind, best.pars)
row.names(best.pars) <- csn[1:7]
colnames(best.pars) <- c("Peak (mu)", "Horiz. scale (sigma)", "Skewness (s)", "Vert. scale (h)")
stargazer(best.pars, type = "html", out = "bestpars.html", summary = FALSE)
 
days.to <- ceiling(cbind(sapply(1:7, function(i) days.to.one[i, best.log[i]]), sapply(1:7, function(i) days.to.tw[i, best.log[i]])))
cases <- ceiling(cbind(sapply(1:7, function(i) cases.one[i, best.log[i]]), sapply(1:7, function(i) cases.tw[i, best.log[i]])))
dcases <- cbind(as.numeric(ds[nrow(ds), 2:length(cs)]), days.to, cases)
colnames(dcases) <- c("Cases by April, 30", "Days till 1 new case/day", "Days till 1 new case/14 days", "Cases till 1 new case/day", "Cases till 1 new case/14 days")
row.names(dcases) <- csn[1:7]
stargazer(dcases, type = "html", out = "days-to-one.html", summary = FALSE)
 
png("plot6-density-logistic-curve-yielding-different-tails.png", 600, 400, pointsize = 16)
par(mar = c(2, 2, 3, 1))
curve(mycurve(x, "logistic", c(10, 3, 2, 1)), 0, 10, n = 1001, col = "blue", lwd = 2, bty = "n", xlab = "n", ylab = "n", main = "Identification problems without the right tail (logistic)", xlim = c(0, 40))
curve(mycurve(x, "logistic", c(10, 6, 4, 1)), 0, 10, n = 1001, add = TRUE, col = "red", lwd = 2)
curve(mycurve(x, "logistic", c(10, 3, 2, 1)), 10, 40, n = 1001, add = TRUE, col = "blue", lwd = 2, lty = 2)
curve(mycurve(x, "logistic", c(10, 6, 4, 1)), 10, 40, n = 1001, add = TRUE, col = "red", lwd = 2, lty = 2)
legend("right", c(expression(sigma == 3 ~ ", " ~ s == 2), expression(sigma == 6 ~ ", " ~ s == 4)), bty = "n", col = c("blue", "red"), lwd = 2)
legend("topleft", c(expression(mu == 10), expression(h == 1)), bty = "n")
dev.off()
 
d.good <- find.day(1, "logistic", c(10, 3, 2, 1))
d.bad <- find.day(1, "logistic", c(10, 6, 4, 1))
get.cases(ceiling(d.good), "logistic", c(10, 3, 2, 1))
get.cases(ceiling(d.bad), "logistic", c(10, 6, 4, 1))
 
criterion.log.rest <- function(pars, skew, observed, fun = function(x) mean(x^2)) {
  theor <- mycurve(x = 1:length(observed), distrib = "logistic", pars = c(pars[1:2], skew, pars[3]))
  return(fun(theor - observed))
}
set.seed(1); cat("Russia logistic: least squares\n")
opt.log.nls <- hydroPSO(par = start.pars[[8]]$logistic, fn = function(x) criterion.log.rest(x, 1, ds[, "drus"], rls), lower = start.pars[[8]]$logistic[c(1, 2, 4)] / 4, upper = start.pars[[8]]$logistic[c(1, 2, 4)] * 4, control = contt3)
set.seed(1); cat("Russia logistic: least absolute deviations\n")
opt.log.nla <- hydroPSO(par = start.pars[[8]]$logistic[c(1, 2, 4)], fn = function(x) criterion.log.rest(x, 1, ds[, "drus"], nla), lower = start.pars[[8]]$logistic[c(1, 2, 4)] / 4, upper = start.pars[[8]]$logistic[c(1, 2, 4)] * 4, control = contt3)
 
rus.nls <- c(opt.log.nls$par[1:2], 1, opt.log.nls$par[3])
rus.nla <- c(opt.log.nla$par[1:2], 1, opt.log.nla$par[3])
rus.nls.cases <- get.cases(nrow(ds), "logistic", rus.nls)
rus.nla.cases <- get.cases(nrow(ds), "logistic", rus.nla)
c(ds[nrow(ds), "rus"], ds[nrow(ds), "rus"] - rus.nls.cases, ds[nrow(ds), "rus"] - rus.nla.cases)
 
# Least absolute deviations is better
 
russia <- matrix(NA, nrow = 8, ncol = 4)
mycols <- rainbow(length(cs), end = 0.8, v = 0.8)
for (i in 1:7) {
  these.pars <- c(rus.nla[1], rus.nla[2] * best.pars[i, 3], best.pars[i, 3], rus.nla[4])
  to1 <- ceiling(find.day(1, "logistic", these.pars))
  to14 <- ceiling(find.day(1/14, "logistic", these.pars))
  c1 <- get.cases(to1, "logistic", these.pars)
  c14 <- get.cases(to14, "logistic", these.pars)
  russia[i, ] <- c(to1, to14, c1, c14)
}
these.pars <- c(rus.nla[1], rus.nla[2] * mean(best.pars[, 3]), mean(best.pars[, 3]), rus.nla[4])
to1 <- ceiling(find.day(1, "logistic", these.pars))
to14 <- ceiling(find.day(1/14, "logistic", these.pars))
c1 <- get.cases(to1, "logistic", these.pars)
c14 <- get.cases(to14, "logistic", these.pars)
russia[8, ] <- c(to1, to14, c1, c14)
 
 
png("plot7-russia-coronavirus-covid-19-cases-forecast-daily.png", 640, 400, pointsize = 14)
plot(ds[, "drus"], xlim = c(1, 305), ylim = c(0, max(ds[, "drus"])), type = "l", lwd = 3, bty = "n", main = "Conditional forecast COVID-19 pandemic in Russia for 2020", xlab = "", ylab = "Daily cases", xaxt = "n")
axis(1, seq(1, 305, length.out = 11), c(month.name[3:12], ""), las = 2)
for (i in 1:7) {
  these.pars <- c(rus.nla[1], rus.nla[2] * best.pars[i, 3], best.pars[i, 3], rus.nla[4])
  lines(nrow(ds):305, mycurve(nrow(ds):305, "logistic", these.pars), lwd = 2, col = mycols[i])
}
these.pars <- c(rus.nla[1], rus.nla[2] * mean(best.pars[, 3]), mean(best.pars[, 3]), rus.nla[4])
lines(nrow(ds):305, mycurve(nrow(ds):305, "logistic", these.pars), lwd = 3, lty = 2, col = "black")
legend("topright", paste0(c(csn[1:7], "Average"), " scenario"), col = c(mycols[1:7], "black"), lwd = 2, lty = c(rep(1, 7), 2), bty = "n")
dev.off()
 
png("plot8-russia-coronavirus-covid-19-cases-forecast-cumulative.png", 640, 400, pointsize = 14)
plot(cumsum(ds[, "drus"]), xlim = c(1, 305), ylim = c(0, max(russia)), type = "l", lwd = 3, bty = "n", main = "Forecast dynamics of COVID-19 pandemic in Russia", xlab = "", ylab = "Total cases, thousands", xaxt = "n", yaxt = "n")
axis(1, seq(1, 305, length.out = 11), c(month.name[3:12], ""), las = 2)
axis(2, seq(0, 6e5, 1e5), as.character(seq(0, 6e2, 1e2)), las = 2)
for (i in 1:7) {
  these.pars <- c(rus.nla[1], rus.nla[2] * best.pars[i, 3], best.pars[i, 3], rus.nla[4])
  lines(nrow(ds):305, cumsum(mycurve(nrow(ds):305, "logistic", these.pars)) + ds[nrow(ds)-1, "rus"], lwd = 2, col = mycols[i])
}
these.pars <- c(rus.nla[1], rus.nla[2] * mean(best.pars[, 3]), mean(best.pars[, 3]), rus.nla[4])
lines(nrow(ds):305, cumsum(mycurve(nrow(ds):305, "logistic", these.pars)) + ds[nrow(ds)-1, "rus"], lwd = 3, lty = 2, col = "black")
legend("topleft", paste0(c(csn[1:7], "Average"), " scenario"), col = c(mycols[1:7], "black"), lwd = 2, lty = c(rep(1, 7), 2), bty = "n")
dev.off()
 
 
mean(best.pars[, 3])
 
colnames(russia) <- colnames(dcases)[2:5]
row.names(russia) <- c(paste0(csn[1:7], " scenario"), paste0("Average scenario (s = ", round(mean(best.pars[, 3]), 2), ")"))
stargazer(russia, type = "html", out = "russia-final.html", summary = FALSE)

References:

Fernández, C., Steel, M.F.J., 1998. On bayesian modeling of fat tails and skewness. Journal of the American Statistical Association 93, 359–371.

Как сделать коэффициенты со звёздочками уровня значимости в LaTeX’е

In this quick entry, it will be shown how it is possible to format regression output tables in LaTeX and automatically print significance codes (in the form of stars) without any symbols placed a priori.

First of all, one should obtain estimates and standard errors from any statistical package (preferably R, but as long as homicide is considered a crime punishable by incarceration or death, and using SPSS is not, one can use whatever they want—although the author wishes it were the other way about).

Second, one should use regular expressions or any other dark wizardry in order to convert such table into a list of LaTeX control sequences taking two argument each, possible delimited by an ampersand. Here is an example how this can be done by taking raw Excel output (tab-delimited) and running any regex engine on it.

Before:

market	0.0552	0.0055	0.00988	0.00149
ret	-1.974	0.0329	-0.767	0.00993
sls	-1.31	0.36	-0.342	0.128
age	-0.0109	0.0246	-0.00381	0.00663
roa	-0.107	0.107	-0.0733	0.0281
Intercept	-0.586	0.138	-0.165	0.036

Latex table before regular expressions

Regular expression:

s/([A-Za-z]+)\t([-0-9.]+)\t([0-9.]+)\t([-0-9.]+)\t([0-9.]+)/\1 & \\coef{\2}{\3} & \\coef{\4}{\5} \\\\/g

After:

market & \coef{0.0552}{0.0055} & \coef{0.00988}{0.00149} \\
ret & \coef{-1.974}{0.0329} & \coef{-0.767}{0.00993} \\
sls & \coef{-1.31}{0.36} & \coef{-0.342}{0.128} \\
age & \coef{-0.0109}{0.0246} & \coef{-0.00381}{0.00663} \\
roa & \coef{-0.107}{0.107} & \coef{-0.0733}{0.0281} \\
Intercept & \coef{-0.586}{0.138} & \coef{-0.165}{0.036} \\

Latex table after regular expressions

Do not be surprised that now there 3 columns, not 5. Suppose one would like to have the standard error placed beneath the coefficient. Then what?

This is the standard preamble.

\documentclass[11pt, a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{microtype}
\usepackage{booktabs}
\usepackage[british]{babel}
\usepackage{fp}
\newlength\coeflength
\setlength\coeflength{2cm}

Now, it is time to write a command that would typeset a fixed-width box (because all entries in a column have the same width). The basic building block will have 3 arguments: width (optional), estimate (numeric), and standard error (numeric).

\newcommand{\simplestuff}[3][\coeflength]{\parbox[t]{#1}{\hfil$#2$}\parbox[t]{#1}{\hfil$(#3)$}}

Now, the stars. fp package provides facilities to work with fixed-point arithmetic, which is more than enough in our case. Three significance levels—10%, 5%, and 1%—and normal distribution quantiles will be used because this was initially written for applied people in accounting. For asymptotically normal estimators, the critical values can be obtained via the R command qnorm(), e. g. qnorm(1-0.001/2), which should give 3.29 as the critical value for a two-sided t test of significance at 0.1% level.

\newcommand{\coef}[3][\coeflength]{%
\FPset\critsss{2.575829}% 0.995 quantile of the standard normal distribution
\FPset\critss{1.959964}% 0.975 quantile of the standard normal distribution
\FPset\crits{1.644854}% 0.95 quantile of the standard normal distribution
\FPset\coefest{#2}% The first argument is the estimate
\FPset\coefse{#3}% The second argument is the standard error
\FPdiv\zstat{\coefest}{\coefse}% This is the observed t-statistic
\FPabs\zabs{\zstat}% This is its absolute value
\FPiflt\zabs{\crits}\def\kostyrkastar{\relax}\else% Minimising the risk of overwriting an existing command
\FPiflt\zabs{\critss}\def\kostyrkastar{*}\else%
\FPiflt\zabs{\critsss}\def\kostyrkastar{**}\else\def\kostyrkastar{***}%
\fi%
\makebox[#1]{\null\hfill$#2$\kern.1em\rlap{${}^{\kostyrkastar}$}}\makebox[#1]{\null\hfill$(#3)$}%
}

Since hard def’s are used for brevity, in order to minimise the chance that the star-drawing command overwrites any existing one, one should consider using their family name or the name of their pet as a command prefix—above is an example of precautiousness rather than hubris.

LaTeX output table coefficient significance stars based on the t statistic

\setlength{\coeflength}{1.5cm}
 
Look how stars appear when the absolute value of the $t$ statistic increases.
 
$\beta_0$ \coef{0.00}{1.0} ($t=0$) \par
$\beta_1$ \coef{1.10}{2.2} ($t=0.5$) \par
$\beta_2$ \coef{2.00}{2.0} ($t=1$) \par
$\beta_3$ \coef{-3.00}{2.0} ($t=-1.5$) \par
$\beta_4$ \coef{5.10}{3.0} ($t=1.7$) \par
$\beta_5$ \coef{-3.80}{2.0} ($t = -1.9$) \par
$\beta_6$ \coef{4.92}{2.5} ($t = 1.968$) \par
$\beta_8$ \coef{19.60}{7.0} ($t = 2.8$) \par
$\beta_9$ \coef{-800}{8.0} ($t = -100$)

LaTeX output table coefficient significance stars in two columns

\begin{table}[htbp]
	\setlength{\coeflength}{2.2cm}
	\setlength{\tabcolsep}{3pt}
	\begin{tabular}{lll}
		\toprule
		& \null\hfill NCSKEW\hfill\null &\null\hfill DUVOL\hfill\null \\
		\midrule
		Market value & \coef{0.0552}{0.0055} & \coef{0.00988}{0.00149} \\
		Return & \coef{-1.974}{0.0329} & \coef{-0.767}{0.00993} \\
		Sales  & \coef{-1.310}{0.36} & \coef{-0.342}{0.128} \\
		Age & \coef{-0.0109}{0.0246} & \coef{-0.00381}{0.00663} \\
		ROA & \coef{-0.107}{0.107} & \coef{-0.0733}{0.0281} \\
		(Constant) & \coef{-0.586}{0.138} & \coef{-0.165}{0.036} \\
		Sector dummies &\null\hfill Yes\hfill\null &\null\hfill Yes\hfill\null \\
		Year dummies &\null\hfill Yes\hfill\null &\null\hfill Yes\hfill\null \\
		\midrule
		Observations &\null\hfill 14\,689\hfill\null &\null\hfill 14\,694\hfill\null \\
		\bottomrule
		\multicolumn{3}{l}{\footnotesize Significance codes: * $p<0.10$, ** $p<0.05$, *** $p<0.01$.} \\[-.5ex]
		\multicolumn{3}{l}{\footnotesize Cluster-robust (by country) standard errors in brackets.} \\
	\end{tabular}
\caption{OLS estimates}
\end{table}

If one is working in particle physics and needs to use sigma levels, one can be more creative with their notation. For example, one might need more markers or vertcial placement. This can be achieved in an instant.

\newcommand{\coefnuclear}[3][\coeflength]{%
	\FPset\coefest{#2}%
	\FPset\coefse{#3}%
	\FPdiv\zstat{\coefest}{\coefse}%
	\FPabs\zabs{\zstat}%
	\FPiflt\zabs{0.5}\def\mystar{\relax}\else\FPiflt\zabs{1}\def\mystar{.}\else%
	\FPiflt\zabs{1.5}\def\mystar{*}\else\FPiflt\zabs{2}\def\mystar{*.}\else%
	\FPiflt\zabs{2.5}\def\mystar{**}\else\FPiflt\zabs{3}\def\mystar{**.}\else%
	\FPiflt\zabs{3.5}\def\mystar{***}\else\FPiflt\zabs{4}\def\mystar{***.}\else%
	\FPiflt\zabs{4.5}\def\mystar{****}\else\FPiflt\zabs{5}\def\mystar{****.}\else
    \def\mystar{\text{Eureka!}}\fi%
	\parbox[t]{#1}{\hfil$#2$\kern.1em\rlap{${}^{\mystar}$}\\%
	\null\hfil\scriptsize$(#3)$}
}

LaTeX output table coefficient significance stars and standard errors underneath

\begin{table}[htbp]
	\setlength{\coeflength}{2cm}
	\setlength{\tabcolsep}{3pt}
	\begin{tabular}{lllll}
		\toprule
		& \null\hfill LHC\hfill\null & \null\hfill Sein\hfill\null & \null\hfill Pectus\hfill\null & \null\hfill Corps-nichons \hfill\null\\
		\midrule
		Experiment 1 & \coefnuclear{4.9908}{1.3099} & \coefnuclear{7.4037}{1.2957} & \coefnuclear{5.6249}{1.3124} & \coefnuclear{1.2275}{1.2893} \\[2ex]
		Experiment 2 & \coefnuclear{4.3236}{1.2972} & \coefnuclear{10.0386}{1.3176} & \coefnuclear{3.7301}{1.3056} & \coefnuclear{6.7275}{1.2844} \\[2ex]
		Experiment 3 & \coefnuclear{6.1690}{1.2955} & \coefnuclear{7.9955}{1.2917} & \coefnuclear{1.8410}{1.2883} & \coefnuclear{8.7442}{1.3116} \\[2ex]
		Experiment 4 & \coefnuclear{0.6228}{1.3083} & \coefnuclear{3.0899}{1.2977} & \coefnuclear{9.3049}{1.3027} & \coefnuclear{2.4693}{1.2962} \\
		\bottomrule
		\multicolumn{3}{l}{\footnotesize Significance codes: ${}^* = 1\sigma$, ${}^. = 0.5\sigma$.} \\
	\end{tabular}
	\caption{Trying to uncover Higgs' bosom}
\end{table}

One can change the placement of \hfill and \null to change the alignment. It is, however, a very delicate matter, and good alignment differs from table to table since sometimes, the decimal separator should serve as the pivot. A good (but rather restrictive in real world) idea would be using the same number of digits for all table entries.

Another good example would be using bootstrapped percentile confidence intervals and checking whether zero is inside it or not. One can supply multiple intervals. Making that kind of output is extremely easy if one uses boot package. However, examples of such style in real scientific journals are virtually unseen—because either their typesetters have not read this entry (yet) and are blissfully unaware of this method, or this sort of acrobatics is not deemed æsthetically pleasing or practical.

\newcommand{\coefbootci}[4][\coeflength]{%
	\FPset\coefest{#2}%
	\FPset\coeflb{#3}%
	\FPset\coefub{#4}%
	\parbox[t]{#1}{$#2$\kern0.5em$[#3\hfill%
		\FPiflt\coeflb{0}\FPifgt\coefub{0}\relax\else**\fi%%
		\hfill#4]$}%
}

\coefbootci[4cm]{0.034}{-0.01}{0.09} \par
\coefbootci[4cm]{3.511}{2.048}{5.194}

It is possible to make any formatting aspect adjustable. The last example shows how to make coefficients significant at 5% bold.

\newcommand{\coefbold}[2]{%
	\FPset\coefest{#1}%
	\FPset\coefse{#2}%
	\FPdiv\zstat{\coefest}{\coefse}%
	\FPabs\zabs{\zstat}%
	\FPiflt\zabs{2}\ensuremath{#1}\else\ensuremath{\mathbf{#1}}\fi%
}

LaTeX conditional formatting based on t statistics

Virgin insignificant coefficient with $t=1.5$: \coefbold{4.20}{2.8} \par
Chad significant coefficient with $t=2.5$: \coefbold{6.66}{2.664}

P.S. In no way is this entry published in order to promote hacking p-values or abusing the data and/or methodology in order to get a sweet succulent p-value that is less than 0.05. Researchers are encouraged to avoid such filthy terms as ‘a convincing trend towards significance’ and other ones, and honestly provide the results as is, using asymptotic refinement via bootstrap, or empirical-likelihood-based with proper calibration, or, in small samples with limited-value variables, exact values for non-parametric small-sample statistics.

Stochastic portfolio optimisation and out-of-sample testing using R

Author: Andreï Victorovitch Kostyrka, University of Luxembourg.

Last updated: 2018-05-31.

One might ask a question, “What is it like to try to pick an optimal portfolio and invest in stocks?” Well, here I have prepared some R code that does the following thing:

  • Downloads stock data from Yahoo Finance from 2014 for a selection of stocks;
  • Puts all stocks together across the maximum set of available dates and nicely imputes missing values using Kalman filter with structural modelling;
  • Tries all possible combinations or portfolios (as many as possible, using Halton pseudo-random low-discrepancy sequences, doing many iterations and saving the best result from every iteration to save memory);
  • Creates amazing colourful plots for subsequent analysis;
  • Returns statistics for equally weighted portfolios, max-return under equally weighted risk, min-risk under equally weighted return, and maximum Sharpe ratio (for simplicity, the ratio of returns to risk, assuming there is no risk-free asset);
  • Tests these 4 portfolios out of sample for 1 year (“if one started trading one year ago, what would one have obtained”).

Two sets of stocks were chosen: most lucrative U.S. companies and most lucrative Russian companies. This is what one gets if one chooses 13 top U.S. companies:

library(BatchGetSymbols)
library(zoo)
library(imputeTS) # For imputation
library(randtoolbox) # For Sobel sampling
library(parallel)
 
# http://fortune.com/2017/06/07/fortune-500-companies-profit-apple-berkshire-hathaway/
# https://www.swiftpassportservices.com/blog/profitable-companies-russia/
 
syms <- c("GOOG", "AAPL", "JPM", "BRK-A", "WFC", "MCD", "BAC", "MSFT", "AMZN", "JNJ", "C", "MO", "GE")
symnames <- c("Google", "Apple", "J.P. Morgan", "Berkshire Hathaway", "Wells Fargo", "McDonald's",
              "Bank of America", "Microsoft", "Amazon", "Johnson and Johnson", "Citigroup", "Altria group", "General Electric")
cat("Calculating optimal portfolio for US stocks\n")
 
# syms <- c("SGTZY", "ROSN.ME", "OGZPY", "LUKOY", "SBRCY", "TRNFP.ME", "NILSY", "ALRS.ME", "YNDX")
# symnames <- c("SurgutNefteGas", "RosNeft", "GazProm", "LukOil", "SberBank", "TransNeft",
#               "Norilsk Nickel", "Alrosa", "Yandex")
# cat("Calculating optimal portfolio for Russian stocks\n")
 
first.date <- as.Date("2014-01-01")
last.date <- as.Date("2018-08-01")
end.ins <- as.Date("2017-08-01")
 
dat <- BatchGetSymbols(tickers = syms, first.date = first.date, last.date = last.date)
 
allt <- dat$df.tickers
allt <- allt[!is.na(allt$price.adjusted), ]
dates <- sort(unique(allt$ref.date))
prices <- data.frame(dates)
prices[, 2:(length(syms)+1)] <- NA
colnames(prices)[2:(length(syms)+1)] <- syms
for (i in syms) {
  datesi <- allt$ref.date[as.character(allt$ticker)==i]
  pricesi <- allt$price.adjusted[as.character(allt$ticker)==i]
  prices[match(datesi, dates), i] <- pricesi
  prices[, i] <- zoo(prices[, i], order.by = prices$dates)
  prices[, i] <- na.kalman(prices[, i])
}
prices2 <- prices
prices[,1] <- NULL
 
co <- rainbow(length(syms), end=0.75)
png("1-relative-prices.png", 700, 400, type="cairo")
plot(prices[,1]/mean(prices[,1]), col=co[1], lwd=1.5, main="Prices divided by their in-sample mean", xlab="Time", ylab="Ratio")
for (i in 2:length(syms)) {
  lines(prices[,i]/mean(prices[,i]), col=co[i])
}
legend("topleft", syms, col=co, lwd=1)
dev.off()
 
rs <- data.frame(rep(NA, nrow(prices)-1))
for (i in syms) {
  rs[, i] <- diff(log(prices[, i]))
}
rs[,1] <- NULL
 
w0 <- rep(1/length(syms), length(syms))
names(w0) <- syms
rmean <- zoo(as.numeric(as.matrix(rs) %*% w0), order.by = time(rs[,1]))
 
a <- table(as.numeric(format(time(rmean), "%Y")))
daysperyear <- mean(a[1:(length(a)-1)]) # Because the last year is incomplete
 
rsfull <- rs # Full sample
rs <- rs[time(rs[,1]) <= end.ins, ] # Training sample
rsnew <- rsfull[time(rsfull[,1]) > end.ins, ] # Test sample
 
meanrs <- colMeans(rs)
covrs <- cov(rs)
 
musigma <- function(w) {
  mu <- sum(meanrs * w)*daysperyear
  s <- sqrt((t(w) %*% covrs %*% w) * daysperyear)
  return(c(mu=mu, sigma=s, ratio=mu/s))
}
 
m0 <- musigma(w0)
 
num.workers <- if (.Platform$OS.type=="windows") 1 else detectCores()
MC <- 108000
iters <- 100
cc <- 150 # How many colours
cols <- rainbow(cc, end = 0.75)
sharpecuts <- seq(-0.1, 1.5, length.out = cc)
 
png("2-portfolios.png", 700, 600, type="cairo")
 
# First iteration
cat("Calculating", MC, "portfolios for the first time.\n")
ws <- halton(MC, length(syms), init=TRUE)
ws <- ws/rowSums(ws)
colnames(ws) <- syms
res <- mclapply(1:MC, function(i) musigma(ws[i, ]), mc.cores = num.workers)
res <- as.data.frame(matrix(unlist(res, use.names = FALSE), nrow=MC, byrow = TRUE))
colnames(res) <- c("mean", "sd", "ratio")
prcmu0 <- ecdf(res$mean)(m0[1]) # Which quantile is the mean
prcsd0 <- ecdf(res$sd)(m0[2]) # Which quantile is the sd
qmplus <- quantile(res$mean, prcmu0+0.005)
qmminus <- quantile(res$mean, prcmu0-0.005)
qsplus <- quantile(res$sd, prcsd0+0.005)
qsminus <- quantile(res$sd, prcsd0-0.005)
meanvicinity <- which(res$mean < qmplus & res$mean > qmminus)
sdvicinity <- which(res$sd < qsplus & res$sd > qsminus)
 
plot(res[,2:1], pch=16, col=cols[as.numeric(cut(res$ratio, breaks = sharpecuts))],
     xlim=range(res$sd)+c(-0.03, 0.03), ylim=pmax(range(res$mean)+c(-0.03, 0.03), c(0, -10)),
     main="Risk and return for various portfolios", xlab="Annualised volatility", ylab="Annualised log-return")
 
sharpes <- mus <- sds <- array(NA, dim=c(iters+1, length(syms)))
sharpes[1, ] <- ws[which.max(res$ratio), ]
mus[1, ] <- ws[which(res$mean == max(res$mean[sdvicinity])), ]
sds[1, ] <- ws[which(res$sd == min(res$sd[meanvicinity])), ]
pairsm <- pairss <- array(NA, dim=c(iters+1, 3))
pairsm[1, ] <- res$mean[c(which.max(res$ratio), which(res$mean == max(res$mean[sdvicinity])), which(res$sd == min(res$sd[meanvicinity])))]
pairss[1, ] <- res$sd[c(which.max(res$ratio), which(res$mean == max(res$mean[sdvicinity])), which(res$sd == min(res$sd[meanvicinity])))]
 
for (i in 1:iters) {
  cat("Calculating", MC, "portfolios, iteration", i, "out of", iters, "\n")
  # ws <- sobol(MC, length(syms), scrambling = 1, init=FALSE) # There is a bug in randtoolbox: it starts returning numbers outside [0, 1]
  ws <- halton(MC, length(syms), init=FALSE)
  ws <- ws/rowSums(ws)
  colnames(ws) <- syms
  res <- mclapply(1:MC, function(i) musigma(ws[i, ]), mc.cores = num.workers)
  res <- as.data.frame(matrix(unlist(res, use.names = FALSE), nrow=MC, byrow = TRUE))
  colnames(res) <- c("mean", "sd", "ratio")
  meanvicinity <- which(res$mean < qmplus & res$mean > qmminus)
  sdvicinity <- which(res$sd < qsplus & res$sd > qsminus)
  maxsh <- which.max(res$ratio)
  maxm <- which(res$mean == max(res$mean[sdvicinity]))
  minsd <- which(res$sd == min(res$sd[meanvicinity]))
  sharpes[i+1, ] <- ws[maxsh, ]
  mus[i+1, ] <- ws[maxm, ]
  sds[i+1, ] <- ws[minsd, ]
  pairsm[i+1, ] <- res$mean[c(maxsh, maxm, minsd)]
  pairss[i+1, ] <- res$sd[c(maxsh, maxm, minsd)]
  points(res[,2:1], pch=16, col=cols[as.numeric(cut(res$ratio, breaks = sharpecuts))])
  rm(res); gc()
}
 
points(as.numeric(pairss), as.numeric(pairsm))
 
sharpems <- apply(sharpes, 1, musigma)
wmaxsharpe <- sharpes[which.max(sharpems[3,]), ]
meanms <- apply(mus, 1, musigma)
wmaxmean <- mus[which.max(meanms[1,]), ] # Maximum returns while SD is the same
sdms <- apply(sds, 1, musigma)
wminsd <- sds[which.min(sdms[2,]), ] # Minimum SD while mean is the same
 
m1 <- musigma(wmaxsharpe)
m2 <- musigma(wmaxmean)
m3 <- musigma(wminsd)
finalstats <- cbind(rbind(w0, wmaxsharpe, wmaxmean, wminsd), rbind(m0, m1, m2, m3))
rownames(finalstats) <- c("Equal", "MaxSharpe", "MaxReturn", "MinRisk")
colnames(finalstats)[(length(syms)+1):(length(syms)+3)] <- c("RInS", "VInS", "SRInS")
 
 
points(m0[2], m0[1], cex=3, pch=16)
points(c(m1[2], m2[2], m3[2]), c(m1[1], m2[1], m3[1]), cex=2, pch=16)
dev.off()
 
# And now checking these revenues on the test sample
r0 <- zoo(as.numeric(as.matrix(rsnew) %*% w0), order.by = time(rsnew[,1]))
r1 <- zoo(as.numeric(as.matrix(rsnew) %*% wmaxsharpe), order.by = time(rsnew[,1]))
r2 <- zoo(as.numeric(as.matrix(rsnew) %*% wmaxmean), order.by = time(rsnew[,1]))
r3 <- zoo(as.numeric(as.matrix(rsnew) %*% wminsd), order.by = time(rsnew[,1]))
 
png("3-test-sample-returns.png", 700, 400, type="cairo")
plot(r0, main="Out-of-sample returns of weighted portfolios", ylab="Return", xlab="Time")
lines(r1, col=2)
lines(r2, col=3)
lines(r3, col=4)
legend("topright", c("Equal", "Sharpe", "Max. return", "Min. volatility"), col=1:4, lwd=1)
dev.off()
 
png("4-test-sample.png", 700, 400, type="cairo")
plot(exp(cumsum(r0)), ylim=c(0.6, 1.4), main="Portfolio growth over the last year in the test sample", ylab="Relative growth", xlab="Time")
abline(h=1)
lines(exp(cumsum(r1)), col=2)
lines(exp(cumsum(r2)), col=3)
lines(exp(cumsum(r3)), col=4)
legend("bottomleft", c("Equal", "Sharpe", "Max. return", "Min. volatility"), col=1:4, lwd=1)
dev.off()
 
finals <- cbind(r0, r1, r2, r3)
meanf <- colMeans(finals)*daysperyear
sdf <- apply(finals, 2, sd)*sqrt(daysperyear)
finalsoos <- cbind(meanf, sdf, meanf/sdf)
colnames(finalsoos) <- c("ROOS", "VOOS", "SRROS")
finalstats <- cbind(finalstats, finalsoos)
colnames(finalstats)[1:length(syms)] <- symnames
 
library(stargazer)
stargazer(t(finalstats), summary = FALSE, type="html", out="a.html")
 
write.csv(prices2, "stocks.csv", row.names = FALSE)

This is how the 13 U.S. stocks behaved.

Relative prices of top U.S. stocks

This is how 10 million random portfolios behaved (black dots indicate 100 best portfolios under the three criteria; big black dot for the equally weighted portfolio).

Return and risk of various U.S. portfolios

These are portfolio weights for U.S. companies with their performance.

Equal MaxSharpe MaxReturn MinRisk
Google 0.077 0.011 0.004 0.052
Apple 0.077 0.132 0.128 0.105
J.P. Morgan 0.077 0.133 0.099 0.028
Berkshire Hathaway 0.077 0.019 0.040 0.170
Wells Fargo 0.077 0.003 0.008 0.025
McDonald's 0.077 0.164 0.148 0.192
Bank of America 0.077 0.030 0.034 0.036
Microsoft 0.077 0.078 0.180 0.022
Amazon 0.077 0.088 0.161 0.005
Johnson & Johnson 0.077 0.138 0.064 0.138
Citigroup 0.077 0.003 0.003 0.009
Altria group 0.077 0.198 0.131 0.142
General Electric 0.077 0.003 0.001 0.078
Return in-sample 0.145 0.179 0.190 0.145
Risk in-sample 0.142 0.128 0.141 0.119
Sharpe ratio in-sample 1.023 1.395 1.345 1.217
Return in 2017.07–2018.06 0.135 0.163 0.241 0.053
Risk in 2017.07–2018.06 0.144 0.134 0.146 0.130
Sharpe ratio in 2017.07–2018.06 0.941 1.218 1.657 0.404

This is how these four U.S. portfolios look out of sample.

Returns for different optimal U.S. portfolios out-of-sample

If one had invested in these four portfolios a year before, this is what one would have harvested by now.

Cumulative growth of one's investments into U.S. stocks


And this is the result for top 9 Russian companies.

This is how the 9 Russian stocks behaved.

Relative prices of top Russian stocks

This is how 10 million random portfolios behaved (black dots indicate 100 best portfolios under the three criteria; big black dot for the equally weighted portfolio).

Return and risk of various Russian portfolios

Note that some portfolios can actually yield negative returns, and the volatility is generally higher!

These are portfolio weights for Russian companies with their performance.

Equal MaxSharpe MaxReturn MinRisk
SurgutNefteGas 0.111 0.005 0.009 0.175
RosNeft 0.111 0.162 0.031 0.269
GazProm 0.111 0.012 0.068 0.067
LukOil 0.111 0.018 0.083 0.037
SberBank 0.111 0.038 0.031 0.014
TransNeft 0.111 0.248 0.016 0.213
Norilsk Nickel 0.111 0.079 0.381 0.044
Alrosa 0.111 0.413 0.358 0.049
Yandex 0.111 0.026 0.022 0.131
Return in-sample 0.055 0.227 0.172 0.055
Risk in-sample 0.245 0.218 0.245 0.218
Sharpe ratio in-sample 0.223 1.042 0.701 0.250
Return in 2017.07–2018.06 0.202 0.164 0.223 0.164
Risk in 2017.07–2018.06 0.179 0.170 0.215 0.153
Sharpe ratio in 2017.07–2018.06 1.129 0.965 1.039 1.067

This is how these four Russian portfolios look out of sample.

Returns for different optimal Russian portfolios out-of-sample

If one had invested in these four portfolios a year before, this is what one would have harvested by now.

Cumulative growth of one's investments into Russian stocks

So actually if one used this investing strategy based on the information available in mid-2017, by now, depending on the criterion (Sharpe ratio or maximum profitability) they would have earned 14–22% on Russian stocks or 15–17% on American stocks!

TODO: constrained optimisation using the best stochastic solution as a starting point.

[2019-05-31] Update: for high-dimensional problems (hundreds of stocks), many more MC trials are required! One might also want to change the Halton sequence to a random uniform generator (like ws <- matrix(runif(length(syms)*MC), ncol = length(syms))) or a Latin hypercube sampler (lhs).

<...>

Why are you still here? Why are you not trading yet?

Using Twitter to model voters’ attitude towards presidential candidates: a case for Russia in 2018

NB. Any opinions that might be critical of the political system of Russia are justified by the fact that the author is a citizen of Russia himself, and retains the right to give opinions, or commentary, or judgement.

Introduction

Twitter has been known to serve as a news publication machine that drives the attention of mass media and shapes the preferences of voters; a good example of this would be the presidential elections in the USA in 2016 (Kapko, 2016). Presidential candidates were using this medium to communicate to voters directly and express their opinions outside presidential debates. It has been speculated that predictions made based on Twitter data can perform better than those obtained by robocalls or by examining the distribution of lawn signs.

Sometimes, as in the case with Donald Trump, using free-wheeling speech style can be really beneficial because it reassures voters that their candidates speaks what (s)he thinks. Candidates' tweets are news per se, and news agencies can easily get the low-hanging fruit by just re-publishing them and giving their own interpretation. Besides that, Twitter provides an opportunity to engage into a direct dialogue between the candidates and voters via the feature of directed tweets. Celebrities and influential people often respond directly to the candidates, which drives the evolution of voters' views. It has been noted that often it is clear who is writing the tweet: the candidate him/herself or his/her staffers.

As good as it may sound, there are obvious drawbacks to Twitter. First and foremost, it seems almost impossible to fit a meaningful argument within the limit of 140 (now 280) characters. As a result, it causes oversimplification of discourse. Another point is debatable: as soon as something is posted on Twitter, it is impossible to fully withdraw the tweet if someone else has noticed it. This can be both a drawback and a strong point: as soon as lewd material is posted, the poster bears full responsibility for it, and experiences all potential backlash.

Twitter as a predictor for election results has been rigorously analysed by researchers in the past. In fact, such analyses were applied to the 2010 Australian elections (Burgess and Bruns, 2012), 2010 Swedish elections (Larsson and Moe, 2012), 2011 Dutch elections (Sang and Bos, 2012), 2012 USA elections (Wang et al., 2012), 2015 Venezuelan elections (Yang, Macdonald and Ounis, 2017), 2016 USA elections (Bessi and Ferrara, 2016) etc. Some of them are using convolutional neural networks, some of them are applying clever tools to avoid misclassification, but some of them rely on a simple count of words from positive and negative dictionaries. We are going to adopt this simple approach in an attempt to analyse Russian data and evaluate the balance of powers (or lack thereof) in Russia.

In this work, I am using big data (70 000 tweets in total) and big data processing tools (Google Translate API) to investigate the relative popularity of candidates in the Russian presidential election of 2018. The methodology will be similar to that of Bhatia (2015).

Methodology

In this work, I develop a simple method to analyse the sentiment towards election candidates in arbitrary languages using Google Translate. It should be noted that as the moment of writing, Google is offering a free $300 trial for its Cloud API, so the author made use of his free access to the machine translation services. The workflow is as follows:

  1. Using R and twitteR package, obtain a large number (10 000) tweets mentioning the candidate;
  2. Break the tweets into small batches (100 tweets per batch separated with triple paragraphs) and write scripts for translation;
  3. Using Google Cloud Translation API, run the translation scripts and obtain the server response as a JSON object;
  4. Input the JSON object into R, break it into separate tweets again and clean up the text;
  5. Count the positive and negative words in those tweets using dictionaries and analyse the distributions.

Obviously, the main drawback of the simple dictionary match is the lack of context. There are examples of why this might return wrong results. Example 1: “Mr. X is going to put an end to corruption!” This author has a positive attitude towards Mr. X, but due to the presence of a negative word (“corruption”), is going to be counted as negative. Example 2: “@newsagency Yes, of course, Mr. Y is going to improve the situation, yeah, like, thousand times!” This is an example of irony; the attitude of this tweet is negative and sceptical, but the simple count will find a positive match (“improve”) and count this tweet as positive. Example 3: “Mr. Z is not honest!”

To counter these three points, an argument should be made. The error from the first category does not change much in terms of emotional attitude; if candidate X is viewed as a corruption fighter, the writer implies that the situation in the country is grave (hence the need to fight), and this negative score should be counted towards the general level of negativity, not just the negativity associated with candidate X. If candidates Y and Z also claim they are going to fight corruption, then all these negative votes will contribute to a fixed proportion of false negatives (say, n%), and the ratio of negative to positive votes should not change. This implies that ratios, not absolute values should be examined. The error from the second category is so semantically complicated that is it improbable that any advanced methods (including machine learning and neural networks) based on pure text analysis are going to detect irony and correctly interpret it in a negative sense. However, the proportion of such tweets is thought to be low, and an advanced user-based analysis is suggested in this case: if a user is generally opposing candidate Y (or is pessimistic in general), then any tweets with positive score about candidate Y should be re-weighted or discarded as “implausible”. The error from the third category, however, can be remedied with the use of multi-layered neural network or contextual analysis algorithms with look-ahead and look-behind: if the word from the “positive” list is preceded by “not”, “no”, “barely” or followed by “my a**” or other similar modifiers, the match should be counted for the “negative” list, and vice versa.

Implementation

Searching for tweets

Assumption 1. During the period from 2018-01-01 to 2018-03-18, any tweet containing the surname of a candidate is meant to be about the candidate, not his namesake.

This assumption allows us to identify tweets about the presidential candidates by a simple search using the name of the candidate. The candidates seem to have quite unique surnames, and almost all media covering people with their surnames are covering namely those people.

Assumption 2. There is no bias in Twitter search results; i. e., Twitter is not being in favour of any candidate, or is not filtering negative search results, or is not affected by spam bots who praise candidates.

This assumption is crucial for the unbiasedness of our results; so far, there have been no signs of Twitter censorship in Russia. However, there might be a population bias: people who are using Twitter can turn out to be younger and more progressive, and therefore, more opposing to the current political leader, Putin. It means that whereas most tweets about Putin on Twitter can come from a minority of political opposition, the majority of Putin voters are budget-depending elderly people and beneficiaries who live in faraway regions without internet access or any desire to express their opinions publicly.

api_key <- "..."
api_secret <- "..."
access_token <- "..."
access_token_secret <- "..."
 
library(twitteR)
setup_twitter_oauth(api_key, api_secret, access_token, access_token_secret)
 
find1 <- "Путин"
find2 <- "Бабурин"
find3 <- "Грудинин"
find4 <- "Жириновский"
find5 <- "Собчак"
find6 <- "Сурайкин"
find7 <- "Титов"
find8 <- "Явлинский"
find9 <- "Навальный"
finds <- c(find1, find2, find3, find4, find5, find6, find7, find8, find9)
 
number <- 10000
 
tweet1 <- searchTwitter(find1, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet2 <- searchTwitter(find2, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet3 <- searchTwitter(find3, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet4 <- searchTwitter(find4, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet5 <- searchTwitter(find5, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet6 <- searchTwitter(find6, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet7 <- searchTwitter(find7, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet8 <- searchTwitter(find8, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
tweet9 <- searchTwitter(find9, n=number, lang="ru", since="2018-01-01", until="2018-03-18")
 
save(tweet1, tweet2, tweet3, tweet4, tweet5, tweet6, tweet7, tweet8, tweet9, file="tweets.RData")

Breaking the tweets for translation

A simple snippet of code is suggested below to translate the long and heavy objects on Google Translate. As we know, Google charges API users on a character basis (as of writing, $20 per 1 million characters), with virtually no limit for translation (if the proper quotas are set to be infinite). However, just in order to be sure that the translations are made in time and no information is lost, we broke the messages into blocks of 100 tweets.

byhow <- 100
 
for (cand in 1:9) {
  texts <- get(paste0("tweet", cand))
  blocks <- ceiling(length(get(paste0("tweet", cand)))/byhow)
  for (i in 1:blocks) {
    sink(paste0("scripts", cand, "/scr", i, ".sh"))
    fr <- byhow*(i-1)+1
    to <- byhow*(i-1)+byhow
    cat('
    curl -s -X POST -H "Content-Type: application/json" \\
    -H "Authorization: Bearer "$(gcloud auth print-access-token) \\
    --data "{
      \'q\': \'
      ')
      if (to > length(texts)) to <- length(texts)
      for (j in fr:to) {
        cat("\n\n\n")
        a <- texts[[j]]$text
        a <- gsub("['\"`]", "", a)
        cat(a)
      }
      cat('
      \',
      \'source\': \'ru\',
      \'target\': \'en\',
      \'format\': \'text\'
    }" "https://translation.googleapis.com/language/translate/v2"
 
    ')
    sink()
  }
 
}

Therefore, for each of the nine candidates, it will write out up to 100 script files into a corresponding folder.

Using Google Translate API

Once the script files, titled scr1.sh, ..., scr100.sh are ready, all that remains is run them in series with a time-out (since Google API has a character-per-minute quota).

#!/bin/bash
for i in {1..100}
  do
  echo "Processing tweet file $i of 100"
  ./scr$i.sh > out$i.txt
  sleep 12
done

This script must be run in every folder (1–9), and the translation results will be stored in files out1.txt, ..., out100.txt (or up to the number of tweets divided by 100 if the candidate was not very popular).

Merging the objects into tweet lists

A simple for loop (for 9 candidates, from 1 to the number of output files) is required to merge the data back.

library(jsonlite)
 
alltrans <- list()
 
for (cand in 1:9) {
  lf <- list.files(paste0("./scripts", cand, "/"))
  lf <- lf[grep("out", lf)]
  howmany <- length(lf)
 
  alltrans[[cand]] <- list()
 
  for (i in 1:howmany) {
    a <- fromJSON(readLines(paste0("scripts", cand, "/out", i, ".txt")))
    a <- a[[1]]$translations$translatedText
    aa <- strsplit(a, "\n\n\n")[[1]]
    aa <- gsub("@ +", "@", aa)
    aa <- gsub("^\n +", "", aa)
    aa <- aa[aa!=""]
 
    alltrans[[cand]][[i]] <- aa
    if (! i%%10) print(paste0("Candidate ", cand, ": ", i, " out of ", howmany))
  }
}
 
clean <- function(t) {
  t <- gsub('[[:punct:]]', '', t)
  t <- gsub('[[:cntrl:]]', '', t) 
  t <- gsub('\\d+', '', t)
  t <- gsub('[[:digit:]]', '', t)
  t <- gsub('@\\w+', '', t)
  t <- gsub('http\\w+', '', t)
  t <- gsub("^\\s+|\\s+$", "", t)
  t <- gsub('"', '', t)
  t <- strsplit(t,"[ \t]+")
  t <- lapply(t, tolower)
  return(t)
}
 
alltrans2 <- lapply(alltrans, unlist)
 
alltrans3 <- lapply(alltrans2, clean)

As we see, this clean-up included removing punctuation, control characters, URLs, and digits from the tweets, and then splitting them at spaces or tabs. For matching, conversion to lowercase was used.

Analysing sentiments

First and foremost, let us look at the distributions of positive and negative tweet counts for various candidates.

positive <- scan('positive-words.txt', what='character', comment.char=';')
negative <- scan('negative-words.txt', what='character', comment.char=';')
 
posscore <- function(tweet) {
pos.match <- match(tweet, positive)
pos.match <- !is.na(pos.match)
pos.score <- sum(pos.match)
return(pos.score)
}
 
negscore <- function(tweet) {
neg.match <- match(tweet, negative)
neg.match <- !is.na(neg.match)
neg.score <- sum(neg.match)
return(neg.score)
}
 
posscores <- list()
negscores <- list()
for (i in 1:9) {
posscores[[i]] <- unlist(lapply(alltrans3[[i]], posscore))
negscores[[i]] <- unlist(lapply(alltrans3[[i]], negscore))
}
 
poscounts <- lapply(posscores, table)
negcounts <- lapply(negscores, table)
 
cols <- rainbow(10)
 
cairo_pdf("counts1.pdf", 7, 5)
plot(NULL, NULL, xlim=c(0, 6.7), ylim=c(0, 7516), main = "Positive tweet counts for candidates", ylab = "Tweets", xlab="Positive words in a tweet")
for (j in 1:9) {
for (i in 1:length(poscounts[[j]])) {
lines(c(i-1 + (j-1)/10, i-1 + (j-1)/10), c(0, poscounts[[j]][i]), lwd=5, col=cols[j])
}
}
legend("topright", finds, lwd=5, col=cols)
dev.off()
 
cairo_pdf("counts2.pdf", 7, 5)
plot(NULL, NULL, xlim=c(0, 6.7), ylim=c(0, 7516), main = "Negative tweet counts for candidates", ylab = "Tweets", xlab="Negative words in a tweet")
for (j in 1:9) {
for (i in 1:length(negcounts[[j]])) {
lines(c(i-1 + (j-1)/10, i-1 + (j-1)/10), c(0, negcounts[[j]][i]), lwd=5, col=cols[j])
}
}
legend("topright", finds, lwd=5, col=cols)
dev.off()

Distribution of positive word counts

Distribution of negative word counts

Fig. 1. Distribution of tweets

Since there is no real political competition to Putin (position 1) among the present candidates (positions 2–8), let us pay some attention to his closest political rival who was denied access to the elections, Navalny (position 9). As we can see from Fig. 1, the number of tweets related to Putin and containing 1–2 positive words is smaller than that related to Navalny. At the very same time, paradoxically, the number of tweets related to Putin and containing 1–2 negative words is also smaller than that of Navalny.

totpos <- numeric(9)
totneg <- numeric(9)
for (i in 1:9) {
totpos[i] <- sum(poscounts[[i]] * as.numeric(names(poscounts[[i]])))
totneg[i] <- sum(negcounts[[i]] * as.numeric(names(negcounts[[i]])))
}
 
# Ordering candidates not alpabetically but as they are in the database
a <- factor(finds)
finds <- factor(a, levels(a)[c(5, 1, 2, 3, 6, 7, 8, 9, 4)])
 
ratios <- totpos/totneg
ratiosdf <- data.frame(ratio=ratios, name=finds)
 
 
library(ggplot2)
p <- ggplot(data=ratiosdf, aes(x=name, y=ratio, fill=name)) + geom_bar(stat="identity") + ylab("Positive/negative ratio") +
theme_minimal() + coord_flip() + scale_x_discrete(limits = rev(levels(ratiosdf$name))) + theme(legend.position="none", axis.title.y=element_blank())
cairo_pdf("ratios.pdf", 8, 7)
plot(p)
dev.off()

Ratio of positive to negative tweet counts

Fig. 2. Tweet sentiment ratios for candidates

It is unsurprising that that the ratio for Navalny, who is known as a fighter against corruption, is lower than that of Putin, despite the fact that there is a substantial number of opposing bloggers, due to the fact that the word “corruption” is negative per se. In general, only one candidate (Grudinin) has this ratio greater or equal to one. This means that tweets containing the names of political leaders are negative in tone in general.

Now it it time to see which words were the top positive and top negative for each candidate.

poswords <- function(tweets) {
pmatch <- match(t, positive)
posw <- positive[pmatch]
posw <- posw[!is.na(posw)]
return(posw)
}
 
negwords <- function(tweets) {
pmatch <- match(t, negative)
negw <- negative[pmatch]
negw <- negw[!is.na(negw)]
return(negw)
}
 
pdata <- list()
ndata <- list()
for (j in 1:9) {
pdata[[j]] <- ""
ndata[[j]] <- ""
for (t in alltrans3[[j]]) {
pdata[[j]] <- c(poswords(t), pdata[[j]])
ndata[[j]] <- c(negwords(t), ndata[[j]])
}
}
 
pfreqs <- lapply(pdata, table)
nfreqs <- lapply(ndata, table)
 
dpfreq <- list()
dnfreq <- list()
for (j in 1:9) {
thresh <- length(alltrans2[[j]])/200
a <- pfreqs[[j]]
a <- a[a>thresh]
b <- nfreqs[[j]]
b <- b[b>thresh]
dpfreq[[j]] <- data.frame(word=as.character(names(a)), freq=as.numeric(a))
dpfreq[[j]]$word <- as.character(dpfreq[[j]]$word)
dnfreq[[j]] <- data.frame(word=as.character(names(b)), freq=as.numeric(b))
dnfreq[[j]]$word <- as.character(dnfreq[[j]]$word)
}
 
for (j in 1:9) {
p1 <- ggplot(dpfreq[[j]], aes(word, freq)) + geom_bar(stat="identity", fill="lightblue") + theme_bw() +
geom_text(aes(word, freq, label=freq), size=4) + labs(x="Major Positive Words", y="Frequency of Occurence",
title=paste0("Major Positive Words and Occurence for \n", finds[j], ", n = ", length(alltrans2[[j]]))) + 
theme(axis.text.x=element_text(angle=90))
p2 <- ggplot(dnfreq[[j]], aes(word, freq)) + geom_bar(stat="identity", fill="pink") + theme_bw() +
geom_text(aes(word, freq, label=freq), size=4) + labs(x="Major Negative Words", y="Frequency of Occurence",
title=paste0("Major Negative Words and Occurence for \n", finds[j], ", n = ", length(alltrans2[[j]]))) + 
theme(axis.text.x=element_text(angle=90))
cairo_pdf(paste0("posw", j, ".pdf"),  2+nrow(dpfreq[[j]])/8, 5)
plot(p1)
dev.off()
cairo_pdf(paste0("negw", j, ".pdf"),  2+nrow(dnfreq[[j]])/8, 5)
plot(p2)
dev.off()
}

Positive words for Putin
Negative words for Putin

Positive words for Baburin
Negative words for Baburin

Positive words for Grudinin
Negative words for Grudinin

Positive words for Zhirinovsky
Negative words for Zhirinovsky

Positive words for Sobchak
Negative words for Sobchak

Positive words for Surajkin
Negative words for Surajkin

Positive words for Titov
Negative words for Titov

Positive words for Yavlinsky
Negative words for Yavlinsky

Positive words for Yavlinsky
Negative words for Navalny

Fig. 3. Most frequent positive and negative words for candidates

We defined “most frequent” words as words whose count exceeded the total number of tweets divided by 200. There is just one small problem: such generic words as “like” and “right”, “well” etc. are on this table. We can remedy that by manually constructing a table of most popular meaningful positive and negative words for candidates (Table 1).

Candidate 1st + word 2nd + word 1st − word 2nd − word
Putin support 275 popular 252 problems 224 protests 153
Baburin important 30 right 27 curses 340 shit 235
Grudinin interesting 416 respect 266 dirty 232 poor 160
Zhirinovsky benefit 225 entertain 206 curses 339 whore 333
Sobchak useful 165 truthful 149 whore 571 beg 410
Surajkin free 86 regard 61 provocation 348 curses 339
Titov works 51 ready 30 shit 284 scandal 237
Yavlinsky liberty 91 great 83 sorry 339 shit 294
Navalny top 232 good 181 strike 355 boycott 352

Table 1. Two most frequent meaningful words for candidates

Conclusions

We can make several conclusions. First, several candidates (Baburin, Surajkin, Titov, Yavlinsky) are so unpopular that they were not mentioned 10 000 times over the course of 2.5 months since the beginning of 2018, which implies that they were added to the list of approved candidates to create the illusion of choice; they cannot compete against the obvious leader of the race. Next, on average, people tweet about politicians when they are discontent with something; it is reflected in the average negative score of tweets, except for Putin and Grudinin. Some politicians (Zhirinovsky, Sobchak) are viewed as politically promiscuous, which is reflected in the word frequency analysis. This all creates the illusion of choice presence, whereas most of the population of Russia does not have access to Twitter, or are not politically engaged, or are dependent on the subsidies issued by the present government, or cannot fully explore the majority of opposing opinions because they are carefully filtered out of the central and official media.

Bibliography

Bessi, Alessandro and Emilio Ferrara (2016). ‘Social bots distort the 2016 US Presidential election online discussion’.

Bhatia, Aankur (2015). Twitter Sentiment Analysis Tutorial. Rpubs.

Burgess, Jean and Axel Bruns (2012). ‘(Not) the Twitter election: the dynamics of the #ausvotes conversation in relation to the Australian media ecology’. Journalism Practice 6.3, pp. 384–402.

Kapko, Matt (2016). Twitter’s impact on 2016 presidential election is unmistakable. CIO.

Larsson, Anders Olof and Hallvard Moe (2012). ‘Studying political microblogging: Twitter users in the 2010 Swedish election campaign’. New Media & Society 14.5, pp. 729–747.

Sang, Erik Tjong Kim and Johan Bos (2012). ‘Predicting the 2011 dutch senate election results with twitter’. Proceedings of the workshop on semantic analysis in social media. Association for Computational Linguistics, pp. 53–60.

Wang, Hao et al. (2012). ‘A system for real-time twitter sentiment analysis of 2012 us presidential election cycle’. Proceedings of the ACL 2012 System Demonstrations. Association for Computational Linguistics, pp. 115–120.

Yang, Xiao, Craig Macdonald and Iadh Ounis (2017). ‘Using word embeddings in twitter election classification’. Information Retrieval Journal, pp. 1–25.

О вреде дробления выборки

Предположим, изначально у нас было одно полное облако точек. Мы разделили его на три непересекающихся подмножества (обозначены цветами). В каждом из цветных облаков есть значимая взаимосвязь между переменными \(X\) и \(Y\), которую обнаруживают любые регрессионные методы (приведены три: классический МНК, квантильная регрессия при \(q=0{,}5\) и локально линейная регрессия первой степени).

Ложные значимые эффекты на подвыборках

Однако если взглянуть на изначальное облако из всех точек всех трёх цветов, то эффекта никакого не будет (эта ничтожно слабая и малозначимая зависимость показана толстой чёрной линией). Получается парадокс: если индивид относится к зелёной подвыборке, то в ней влияние \(X\) на \(Y\) есть. То же верно и для всех остальных подвыборок. Везде эффект положительный! Однако одновременно для всех точек никакого эффекта нет.

Зависимая переменная: Y
Вся выборка Синие Зелёные Красные
X 0.013 0.839*** 0.482*** 0.598***
(0.030) (0.024) (0.048) (0.046)
Константа 0.004 0.013 -1.202*** 1.245***
(0.031) (0.019) (0.055) (0.053)
Наблюдений 1000 386 307 307
R2 0.0002 0.756 0.252 0.361
Adjusted R2 -0.001 0.755 0.250 0.359
Residual Std. Error 0.981 (df = 998) 0.373 (df = 384) 0.667 (df = 305) 0.640 (df = 305)
F Statistic 0.198 (df = 1; 998) 1,187.609*** (df = 1; 384) 102.762*** (df = 1; 305) 172.266*** (df = 1; 305)
Note: *p<0.1; **p<0.05; ***p<0.01

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

Код для репликации в R (отчего-то плагин глючит и склеивает всё в строку):

n <- 1000; set.seed(100); x <- rnorm(n); y <- rnorm(n); th <- 0.7; g1 <- (y < x + th) & (y > x - th);
g2 <- (y <= x - th); g3 <- (y >= x + th);
m <- lm(y~x); m1 <- lm(y~x, subset=g1); m2 <- lm(y~x, subset=g2); m3 <- lm(y~x, subset=g3)

seed не должен быть ручным

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

Правило 2. Если вам нужно получить несколько независимых случайных величин, не устанавливайте один и тот же seed. Для воспроизводимости принципиально установить его вообще. Если в исследовании требуются три случайные переменные: \(X_1\), \(X_2\), \(X_3\) — и чтобы каждая имела распределение \(\mathcal{U}[\underline{x}\vphantom{x}_i;\overline{x}\vphantom{x}_i]\), ни в коем случае не пишите так:

set seed 100
series K=randgen(u,10,30)
set seed 100
series L=randgen(u,20,40)
set seed 100
series M=randgen(u,-50,50)

Это приведёт к вырождению процесса и к тому, что между случайными величинами будет абсолютная мультиколлинеарность (они будут равны друг другу с точностью до линейного преобразования), а коэффициент корреляции будет равен 1. Это ещё хуже, чем компьютер RANDU, у которого «каждое число случайно само по себе, но не гарантируется того же для большего их количества», а единичная корреляция следует из соотношения \(x_{k+2}=6x_{k+1}-9x_{k}\).

degen-random

Хороший эконометрист всегда задаёт большие и как можно более отличные друг от друга сиды и записывает их в надёжное место. Отличный эконометрист генерирует случайные сиды и записывает начальный сид для генератора сидов; отличный эконометрист всегда должен мыслить на один уровень случайности выше.

set seed 100
series K1=randgen(u,10,30)
set seed 200
series L1=randgen(u,20,40)
set seed 300
series M1=randgen(u,-50,50)

A truly random 3-dimensional uniform distribution

Не надо думать, что зиккурат-алгоритм или алгоритм Бокса—Мюллера вывезет. Не вывезет!

Если вам необходимо прогнать цикл из 10 000 симуляций, установите сид для первой симуляции и не трогайте его, пока не закончится десятитысячная. У современных генераторов достаточно длинные периоды, чтобы не зациклиться даже на миллионе циклов. Алгоритм KISS («Keep it Simple Stupid», или «Не усложняй, придурок») обладает периодом, бо́льшим на 40 порядков, чем число атомов во Вселенной. Если вам необходимо воспроизвести 6 234-ю итерацию, не надо задавать вручную все 10 000, а затем выбирать 6 234-й, нет! Достаточно после каждого цикла в некоторую переменную записывать состояние датчика случайных чисел (random number generator state, или начальное значение, которое он выбрал себе сам для новой итерации), а затем для воспроизведения 6234-й итерации взять seed из полученной переменной. Не обесчещивайте славное имя Джорджа Марсальи, посмотрите в руководство пользователя или исходный код, найдите переменную, в которую записывается состояние датчика, и копируйте её значения в летописи.

Вторая причина, по которой не надо устанавливать seed вручную, — это плохая «случайность» чисел, получаемых с ручным сидом. Если вы получили 50 000 значений случайной переменной и хотите получить ещё 50 000, надо быть волшебником, чтобы при выставлении второго сида не наткнуться на начальное значение, которое уже проскакивало в первом цикле итераций, иначе части первой и второй последовательностей псевдослучайных чисел могут совпасть! Это смертный приговор специалисту по временным рядам или исследователю систем бинарных уравнений, где у ошибок-инноваций не должно быть корреляции со всеми предыдущими значениями случайной переменной. А у вас стоит идеальная машина, изготовленная по технологическому процессу 22 нм (диаметр среднего атома всего лишь в сто раз меньше), куда подаётся ток совершенно определённой силы, поэтому положение каждого электрона в машине детерминировано, и это не даёт возможности получить истинно случайные числа. Вот что вы будете делать с абсолютной корреляцией и единичными корнями, которые возникли только потому, что вы своим указующим перстом установили seed, совпавший с тем, который использовался для 49 998-го значения? Л’Экюйера и Симара на вас не хватает (TestU01, Pierre L’Ecuyer, Richard Simard).

Википедисты шутят — типографы смеются

Необычный пример злоупотребления «чёрточками», с которым довольно трудно что-либо поделать. Пример из статьи «Брана» (по состоянию на 05.09.13):

Membrane punctuation

Ещё раз посмотрим на это. «Инстантон — (−1)-брана». Тире, минус, дефис. Скажи, если бы они все были одинаковой длины, как бы ты их различил?

Не будь портных, — скажи: как различил бы ты служебные ведомства?

К. Прутков.

Задача по теории вероятностей

Задача. Передача «Что? Где? Когда?». На столе перед знатоками лежат по кругу тринадцать писем от телезрителей. Известно, что после того, как письмо отыграло, его снимают со стола, а на его месте остаётся стрелка, по которой осуществляется переход к следующему по часовой стрелке письму. Распорядитель крутит волчок, который выбирает письмо номер ноль для первого раунда. Оно убирается, на его месте остаётся стрелка, и далее по часовой стрелке письма от него нумеруются как первое, второе, третье, ..., двенадцатое.

  • Рассчитайте вероятность того, что в игре будет участвовать письмо с номером \(n\) для \(n=1,\ldots,12\).
  • Какова вероятность того, что в течение следующих двух ходов выпадет письмо 4, если уже выпали письма номер 1 и 2?
  • Какова вероятность того, что письмо 1 будет разыграно последним?

Нетрудная математическая задачка

Задача. Найти такие многочлены \(u(x)\), \(v(x)\), что
\[x^m u(x) + (1-x)^n v(x) =1. \]


Решение. Вспомним разложение бинома и запишем его сразу в форме суммы:
\[ 1 = \bigl( x + (1-x) \bigr)^{n+m-1} = \sum\limits_{i=0}^{n+m-1} C^i_{m+n-1} x^i (1-x)^{m+n-1-i} = 1. \]

Разобьём эту сумму от \(i=0\) до \(i=m+n-1\) на две суммы: от \(i=0\) до \(i=m-1\) и от \(i=m\) до \(i=m+n-1\). Из суммы \(\sum\limits_{i=0}^{i=m-1}\) можно будет вынести за знак суммирования \((1-x)^n\), а из суммы \(\sum\limits_{i=m}^{i=m+n-1}\) отлично выносится \(x^m\): эти выражения присутствовали в исходном вопросе. Для удобства сумму, из которой выносится \(x^m\), будем делать не по \(i\), а по \(j\), а сумму с вынесенным \((1-x)\) сделаем по \(k\). Следует пояснить, что \(\sum\limits_{j=m}^{j=m+n-1}\) — это то же, что и \(\sum\limits_{j=0}^{j=n-1}\). Тогда
\[\sum\limits_{i=0}^{n+m-1} C^i_{m+n-1} x^i (1-x)^{m+n-1-i} = \]
\[ = \sum\limits_{j=m}^{m+n-1} C^j_{m+n-1} x^j (1-x)^{m+n-1-j} + \sum\limits_{k=0}^{m-1} C^k_{m+n-1} x^k (1-x)^{m+n-1-k} = \]
\[ = x^m \sum\limits_{j=0}^{n-1} C^{j+m}_{m+n-1} x^j (1-x)^{n-1-j} + (1-x)^n \sum\limits_{k=0}^{m-1} C^k_{m+n-1} x^k (1-x)^{m-1-k} = \]
\[ = x^m u(x) + (1-x)^n v(x) \equiv 1. \]

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

[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

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

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

Заключённые и телефоны-1

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

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

Перед началом тюремного срока заключённые могут встретиться и обсудить свою стратегию, какой цвет следующего телефона сообщать тюремщику и когда заявлять о полном обзвоне всех заключённых. После этого они не смогут менять стратегию, так как камеры находятся в разных местах и не пропускают звука.

Тюремщик звонит с переменным интервалом и неравномерно, поэтому время между двумя звонками всегда разное, и способа измерения времени не существует. Заключённые обладают хорошей памятью, то есть каждый помнит, на какой телефон сколько раз ему звонили. О какой стратегии необходимо договориться заключённым, чтобы выйти из тюрьмы через максимально короткое время?

Аттолисы (атто-лисы)

В зарубежной терминологии есть чудесный термин «atto-fox» \((\mathrm{fox}\times10^{-18})\), которым обозначается одна из проблем, существующих в модели Лотки—Вольтерры (Lotka—Volterra model), описывающей динамику популяции хищников и жертв. Если быть кратким, чего я категорически не люблю, то в этой модели жертвы (травоядные) обеспечены бесконечным числом ресурса и плодятся экспоненциально, хищники обеспечены только жертвами, но если жертв слишком мало, то хищники вымирают, и от этого возобновляется численность жертв. Модель описывает цикл двух переменных — численности жертв и численности хищников.

Проблема «аттолис» заключается в том, что похожая модель, описывавшая динамику бешенства лис в Великобритании, допускала непрерывные значения переменной. Так, в этой модели в определённые промежутки времени численность лис может составлять \(10^{-18}\) лисы на квадратный километр (при изначально целочисленной популяции, например миллион лис). С таким размахом значений (24 порядка!) жди проблем!

В модели, рассмотренной в работе «On the Spatial Spread of Rabies among Foxes» (J. D. Murray; E. A. Stanley; D. L. Brown — Proceedings of the Royal Society of London. Series B. Biological Sciences, Vol. 229, No. 1255 (Nov. 22, 1986), 111–150.) Динамика численности лис \(S\) при отсутствии бешенства описывается следующей моделью:
\[ \frac{\mathrm{d}S}{\mathrm{d}T} = (a-b)\left( 1-\frac{S}{K}S \right) \]
\(T\) — время, \(a\) — рождаемость, \(b\) — смертность, \(K\) — ёмкость среды.

Пусть бешенство передаётся от бешеных лис здоровым лисам. Тогда хищниками, пожирающими жертв, становятся бешеные лисы, а жертвами — обычные лисы. Скорость заражения на одну лису равна \(\beta R\), где \(R\) — численность бешеных (Rabid) лис.

Заражённые лисы становятся заразными со скоростью \(\sigma\), где \(\1/sigma\) — инкубационный период.

[Дописать!]

Модель падающего тела при измерении расстояния, часть 1

В школе часто задают задачку: рассчитать глубину колодца, если в наличии есть только камень и секундомер. При этом предлагается воспользоваться простейшей формулой для измерения координаты тела, падающего равноускоренно без нулевой скорости:
\[ h = \frac{gt^2}{2}, \]

где \(h\) — путь (глубина колодца), \(g\) — ускорение свободного падения, \(t\) — время на секундомере.

При этом упускается важнейшая деталь. Скорость звука на этой планете очень невысока. Естественно, как только камень по параболе упадёт на дно, звук со дна начнёт лететь равномерно и прямолинейно, безо всякого замедления. Ушей слушателя он достигнет через время, равно отношению глубины ущелья к скорости звука. Если глубина невелика, то тогда всё нормально. Однако при глубине в 123 метра (или высоте, что то же самое; 123 метра — это высота 33-этажного MCI Centre в Лос-Анджелесе... или любого другого здания подобной этажности) звуку потребуется более трети секунды, чтобы достичь ушей слушателя. Как видно из иллюстраций, дальше — больше. Расстояние квадратично зависит от времени, поэтому доля возврата звука обратно в общем измеренном времени неуклонно растёт.

1s

5s

15s

42s

Если бы мы в модели с конечной скоростью звука просверлили в центре Эвереста отверстие-тоннель до земли и бросили туда камень, то он бы падал примерно 42 секунды, а обратно звук бы дошёл за половину этого времени (!). Радостные, мы бы стояли с секундомером, получили бы примерно минуту, подставили бы её в формулу \(h = \frac{gt^2}{2}\) и вычислили бы высоту Эвереста, которая бы отличалась от истинной в бо́льшую сторону... всего в несколько раз (иными словами, мы бы определили высоту горы с точностью до константы).

В формуле расчёта расстояния требуется ввести корректировку на скорость звука. Если этого не сделать, с каждой следующей секундой относительная погрешность измерения будет увеличиваться на 1,8–2 процентных пункта. Поэтому скорректируем наш ответ с учётом знания того, что скорость звука конечна и равна некой \(c\).

Звук летит равномерно прямолинейно без ускорения. Человек с секундомером измеряет сумму двух временных промежутков: время полёта камня (\(t_1\)) и время возврата звука (\(t_2\)). То есть мы знаем \(t = t_1+t_2\), а в формулу высоты, с которой падает тело, надо подставлять \(t=t_1\): \(h = \frac{gt_1^2}{2}\). Более того, \(h = \frac{gt_1^2}{2} = ct_2\), \(t=t_1+t_2\). Выразим \(t_1\) через \(t\):
\[ h = \frac{gt^2_1}{2} = ct_2 \quad \leadsto \quad t_2 = \frac{gt_1^2}{2c} \]
\[t = t_1 + \frac{gt_1^2}{2c} \quad \leadsto \quad \frac{g}{2c}t_1^2 + t_1 - t = 0\]
\[D = b^2-4ac = 1 + \frac{2gt}{c} \quad \leadsto \quad t_1 = \frac{c}{g}\left( \sqrt{1 + \frac{2gt}{c}} -1 \right)\]

Скорость звука (о ней будет сказано ниже) нам известна (≈340 м/с), ускорение свободного падения тоже известно (9,80665 м/с²). Поэтому итоговая формула, позволяющая определить высоту объекта по времени падения камня, который стучит о дно, будет иметь вид:
\[ h = \frac{gt_1^2}{2} = \frac{c^2}{2g}\left( \sqrt{1 + \frac{2gt}{c}} -1 \right)^2 = c\left( t - \frac{c}{g} \left( \sqrt{1 + \frac{2gt}{c}} +1 \right) \right).\]

Легко убедиться в том, что при очень большой скорости звука время \(t_1\) будет стремиться ко времени \(t\):
\[ \lim\limits_{c\rightarrow\infty} \frac{c}{g}\left( \sqrt{1 + \frac{2gt}{c}} -1 \right) = \lim\limits_{c\rightarrow\infty} \frac{c}{g}\frac{\left( \sqrt{1 + \frac{2gt}{c}} -1 \right)\left( \sqrt{1 + \frac{2gt}{c}} +1 \right)}{\left( \sqrt{1 + \frac{2gt}{c}} +1 \right)} = \]
\[= \lim\limits_{c\rightarrow\infty} \frac{c}{g}\frac{2gt}{c\left( \sqrt{1 + \frac{2gt}{c}} +1 \right)} = \lim\limits_{c\rightarrow\infty} \frac{2t}{ \sqrt{1 + \left\{\frac{2gt}{c_{\rightarrow\infty}} \right\}}_{\rightarrow0} +1 } = \lim\limits_{c\rightarrow\infty} \frac{2t}{\sqrt1 + 1} = t.\]

Вашему вниманию предлагаются графики, демонстрирующие отклонение результатов подсчёта по «школьной» и по скорректированной формуле.

diff1

diff2

diff3

Всем хороша данная модель, только сопротивления воздуха она не учитывает. Во второй части данной статьи я введу более совершенную модель, в которой будет учтена температура воздуха, его плотность, а также масса и аэродинамические свойства падающего объекта.

Сущность математики как Вселенной

диалог двух современных греков с замашками на монументальность

Анакрит. Диоклез, вы такой мечтатель! Знаете, люди (в том числе и вы) обречены на увлечение астрологией, если такого ещё не было.

Диоклез. Я? Обречён? Ни за что. Аз есмь сциентист. И физическим тяготением, слабым, сильным и электромагнитным взаимодействиями, сходящимися к одному в М-теории, всё буду объяснять. А не звёздами!

Анакрит. Подлинный сциентизм всегда возвращается, но никогда не бывает постоянным. Ему всегда сопутствует какая-нибудь лженаучная мистика. Сейчас у вас крайняя форма иррациональности — математика головного мозга. При должном лечении это переходит в более адекватную и реалистичную астрологию. Вы знакомы с биографией Юлиана Отступника? А со взглядами Кампанеллы?

Диоклез. «Астро» — это звезда. Вы хотите сказать, что я начну видеть мистику в расположении звёзд, которые движутся по инерции, а каждое следующее положение зависит от текущего, поэтому ничего изменить нельзя?

Анакрит. Верно! Только ведь науке неважно, что первично: инерция двигает звёзды, или звёзды порождают инерцию. Уравнение не изменится. Остаётся свобода трактовки. А вам, эгоцентричному символисту, астрологическое толкование окружающих явлений ох как пойдёт. Тем более что подлинная астрология — это не гороскопчики для вульгарных дамочек, а манускрипты на древних языках, перегрузка мозговых центров и математический анализ.

Диоклез. Звёзды не могут порождать инерцию. Импульс сидит в атомах. Потенциальная энергия сообщается телу, а не возникает в нём: за волосы вы себя из болота не вытянете. А эзотерика мне действительно пойдёт. Но продешевляться не желаю.

Анакрит. Кто сказал, что импульс сидит в атомах? Это утверждение с точки зрения позитивизма бессмысленно. Есть наблюдение, что за определённым приложением энергии следует определённое движение тела, и для описания этого процесс был введён термин «импульс». Это техническое понятие, вспомогательное и самостоятельной природы не имеющее. А дальше можно чему угодно его приписать — звёздам, атомам, Богу; это уже не наука. Тем более что в физике вроде как не признаётся уже принцип локальности. Хотя я этого не принимаю, как и индетерминизма.

Диоклез. Кинетическая энергия — это волна, распределённая в пространстве-времени. Возмущение гравитационного поля. Всё это можно измерить и определить.

Анакрит. Но опять же, «волна» — вспомогательное понятие. Непосредственно фиксируются изменения в строении вещества, для их объяснения вводится понятие волны, в своём роде техническое. Тем более что наиболее известный вид материи на грани вещества и волны — свет — является непосредственным продуктом деятельности звёзд.

Диоклез. Естественно, принцип локальности неверен. Взаимодействие известных сил может ослабляться экспоненциально, сила на расстоянии нескольких световых лет может быть обратна числу Грэма, но нулю она никогда не равна.

Анакрит. Ах, как вы красиво выразились: «обратно числу Грэма».

Диоклез. Свет — это не продукт деятельности звёзд, чёрт побери! Нет у них никакой деятельности. Он просто порождается наличием зарядов и их движением. Какую-то часть мы видим. Какую-то можем измерить. Энергия света и масса Вселенной связаны уравнением Эйнштейна, которое и дефект масс объясняет, и материю и энергию отождествляет. Допустите наличие тёмной материи — и всё у вас сойдётся в ноль. Рассчитайте искажение пространства-времени в сверхплотном сверхкомпакте — сможете потом определить, где будет рассеяна энергия, где будет масса, как тела будут двигаться и с какой силой каждое взаимодействие будет проявляться. И знаете что? Всё это уже было определено в момент зарождения Вселенной. Всё остальное — увы, мой друг, инерция.

Анакрит. Что значит «не продукт деятельности»? Бóльшая часть Вселенной что по массе, что по пространству — звездообразные структуры и вакуум (тёмную энергию пока не берём в расчёт). Свет преимущественно исходит из звёзд.

Кстати, всё определено! Что бы там ни говорили нынешние индетерминисты, людишки без философского мышления, не понимающие, что экспериментально опровергнуть детерминизм никак невозможно.

Диоклез. Да, свет исходит из звёзд. Но те просто зародились и продолжают излучать энергию по инерции. Звезда не имеет мозга или производственной функции, которая что-то с параметрами делает. Она просто теряет свет — вот и всё.

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

Анакрит. Опять метафизика без эмпирического смысла. Материального мозга нет, естественно. А «теряет свет», «производит свет», «дарует свет» — синонимичные высказывания с позитивистской точки зрения.

Диоклез. Да, это синонимы. Но я не об этом. Я о том, что звёзды не порождают инерцию, но являются её рабами.

Я предлагаю ради издевательства создать какую-нибудь свободную и открытую интернет-энциклопедию с названием навроде «Трещины на камнях» или «Каталог камешков менее 1 см». И подробно описывать, где камень лежал в таком-то месяце, куда переместился. Мол, мы начали исчислять Вселенную с этого, получаем такие данные об устройстве мира. Начинаем с малого. Создать тысячи записей, фотографий, биографий камней. Для этого надо всего лишь создать хороший словарь, шаблонизатор — пожалуйста, каталог материи Вселенной (часть первая: камни Земли менее 1 см) готов.

Анакрит. Верно. Но это наша задача как учёных — описание фактов движения. Вопрос в том, на каком языке описывать. А это как раз отлично ложится на неоплатоническую концепцию! И тёмная энергия хорошо впишется. Теория трёх солнц Юлиана вам известна?

Диоклез. О трёх солнцах знаю. Высшее Солнце является неточным аналогом моей математики, которая является божеством, как вы когда-то выразились. Вселенское ничто и объективность, то, что только отражается, но что справедливо всегда и присно — это моя математика, христианский Бог и абсолютное Солнце Юлиана.


Анакрит. Помню карикатуру «Факты: журнал для неопозитивистов». И там «\(A = A\)», «\(5 > 3\)» и так далее — на каждой странице.

Раз так грезите идеей точного моделирования Вселенной, которую я оставил года три назад и против которой могу привести сильные аргументы, вот вам задача. Гипотеза: ни в какой возможной вселенной не может существовать машина, вычисляющая точно состояние этой вселенной в некоторый будущий момент времени и при этом не равная по размеру этой вселенной (то есть сама вселенная моделью самой себя считаться не может). Доказать или опровергнуть — подумайте над задачей. Очевидно, что гипотеза верна, но требуется доказательство.

Диоклез. Логарифмы и законы природы — это у меня нечто исчислимое, у религиозников — веруемое, у Юлиана — недоступное. Но я лучше всех их, потому что я предлагаю это всё предельно аппроксимировать, провести эксперименты в бесконечном числе вселенных, и тогда мы получим вот этот самый абсолютный закон всего существующего.

Анакрит. Математика — это второе Солнце, Царь-солнце Юлиана. Высшее недоступно.

Диоклез. Оно доступно. Мы же только открываем математику, раскрываем квантовые дуальности, коробку Шрёдингера, а если будет бесконечное число таких откровений, то тогда высшее станет доступным. Проблема: высшее — это и есть вся эта плеяда бесконечных вселенных. Это будет полностью познанный мир разума, восприятия. Восприятие нетождественно реальности. Разум нетождественен восприятию. Восприятье есть реальность.

Анакрит. Что? Что же такое разум, если не восприятие? Как по мне, это лишняя категория: есть только восприятие, но нет воспринимающего. Позвольте, с чего восприятие есть реальность? То есть невоспринимаемое не есть реальность, как у Беркли? Вы субъективный идеалист?

Диоклез. Под разумом я подразумеваю законы. Только мы их несовершенными органами воспринимаем.

Анакрит. Какое отношение законы имеют к разуму? Вы как раз постоянно утверждали, что законы и без разума прекрасно существуют.

Диоклез. И реальность несовершенна. Орбитали квантуемы, а не непрерывны. Спины рациональные. Я имею в виду, что вся совокупность законов Вселенной — это бесконечный разум. А я хочу быть частицей с комплексным спином!

Анакрит. «Хочу!» — это что за субъективщина?! При чём тут «несовершенство»? Это несовершенство восприятия и его языка. Это хорошо сформулировал Докинз: человек развивался преимущественно на мезоуровне, поэтому в его мозгу нет частей, отвечающих за восприятие микро- и макроявлений.

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

Анакрит. Ага, понял вас. Чисто неоплатоническое словоупотребление; нус (примерно равно разуму) у Плотина близок к чему-то вроде совокупности законов, идее, математике нашего мира. Я употребляю разум в смысле восприятия и стремления нас к этим законам, а не в смысле самих законов.

Диоклез. А я хочу, потому что в идеальной Вселенной — в этой самой бесконечности — в совершенном Солнце — я смогу быть кем угодно! И с дуальностью, и на орбитали номер корень из пи. И порядок будет непрерывный. И поставлю себе на стол 1,5 монитора, и они будут все работать! А сам я буду находиться в минус трёх с четвертью местах одновременно.

Анакрит. Вы как-то слишком привязаны к материальному состоянию мира. Давайте определимся, что есть для вас математика. Напишите про каждое последующее определение, что думаете о нём:

1. Математика — в первую очередь вспомогательный инструмент, помогающий нам неточно приблизить наши представления о материальной реальности к ней самой.

2. Математика — совокупность сущностных законов, эссенция нашего мира, которую мы познаём через исследование этого мира (поиск единой теории).

3. Математика — набор абсолютных законов, стоящий вне и выше нашего мира и постигаемый разумом, ценный сам по себе.

Диоклез. Второе и третье верны одновременно. Как-то вы по-юлиановски написали...

Анакрит. В первом случае математика рукотворна и нужна для мира, во втором нерукотворна и мир нужен ей, в третьем нерукотворна и ни в ком не нуждается. Да и первое верно, просто это немного разные математики.

Диоклез. Всё, что стоит выше, реализуется в частной выборке. Мы наблюдаем что-то ещё в одной вселенной — приближаем среднее. И так определяем математику, которая нерукотворна. Мир ей не нужен. Просто получилось, что мы есть, поэтому можем попытаться её определить. И в бесконечности у нас это получится. Только мир математике не нужен.

Анакрит. Тогда смотрите. Я утверждаю, что формула \(e=mc^2\) относится ко второй математике, 2 яблока + 2 яблока = 4 яблока — к первой, а относительно третьей математики никакое высказывание не будет отличаться от другого, что-то вроде \(0 = 0\).

Не перескакивайте грани бесконечности! Вы меняете одно выборочное среднее на другое. Ни на какой выборке истинное среднее вы не сможете найти.

Диоклез. Это субъективный вопрос мира — попытаться что-то определить, и эту математику в том числе. Над ней ничего нет. Она неизменна. Но в бесконечности полностью вычислима. Тогда уже исчезнет смысл существования, Вселенная познает сама себя, всё кончится.

Вы очень хорошо всё охарактеризовали, только не вполне точно. Дело в том, что \(e=mc^2\), \(-e=-mc^2\), сложим соответствующие части, получим ноль. В сумме ноль равен нулю. Но это такой сложный нейтральный механизм!.. То, что \(x = 1\times x\), а \(0=x\times 0\), — это далеко не тривиально! И там столько механизмов, определяющих ноль... Их бесконечно много! Но все они познаются за бесконечный промежуток времени.

Анакрит. Мне кажется, уравнение — не лучшее, что физика даст математике. Я хотел бы видеть единую теорию не в виде уравнения, а в виде метода, который, имея на входе координаты и импульсы в одной точке времени, выдавал бы их в другой точке.

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

Диоклез. Тогда уже появятся механизмы, позволяющие задействовать части бесконечности, чтобы континуализировать и кватернионировать измерения, перпендикуляры, и мы сможем запустить механику части бесконечности. Например, четверти бесконечности. То есть спустя четверть бесконечности времени. И тогда Вселенная на четверть реализует себя.

Анакрит. «Мы сможем»?

Диоклез. Мы сможем! Механизм бесконечной генерализации механизмов может лежать в другой четверти, в другом орте. Но этому мехзанизму ещё что-то нужно.

Анакрит. Когда «мы сможем»? Я полагаю, через четверть бесконечности лет?

Диоклез. Да, через четверть бесконечности. Так. Метод должен быть максимально формализован и квантован. Должна появиться какая-то надматематика — генерализация вашей второй математики.

Анакрит. Мне одно не нравится: вы слишком уверенно пишете. Пишете о принципиальной ограниченности познания, а тон уверенный. Непроявление уважения к ней, или правда не чувствуете нерушимой грани между рассуждениями и Тем?

Диоклез. Не перебивайте! Надматематика будет комплексно относиться... Даже не знаю. Нет в текущей математике такой операции, которая бы описывала взаимоотношение генерализации средней математики и абсолютной. Какое-то... подобие, что ли.

Анакрит. Знаете, Александр Богданов свою науку тектологию, общую теорию организации и предшественницу кибернетики, подавал как генерализацию математики.

Диоклез. Просмотрел-с. Недостаточно общо. Дело в том, что математика и физика — это одно и то же. Познание ограничено только в дискретном и конечном промежутке. В бесконечности дисперсия среднего будет равна нулю, и средние значения отразят истинные, будут им тождественны.

Анакрит. Тектологией вы тоже обречены увлечься, потому пока не настаиваю.

Диоклез. Алгебра бесконечномерна. И с ней переплетены и ей тождественны физические законы мироустройства. Поэтому моё желание — стать кем-нибудь в −3,24 местах и одновременно на 1,5 мониторах печатать с отрицательной скоростью. Желание моё полностью законно как реализация частного исхода этой бесконечной системы бесконечностей, которая подарит нам возможность воспроизводить все те механизмы, которые мы откроем.

Анакрит. Выборочное среднее на бесконечности наблюдений не тождественно истинному среднему. Выборочное остаётся случайной величиной. Бесконечность может быть бесконечностью ошибок.

Диоклез. Бесконечная выборка — это генеральная совокупность, и потому среднее по ней истинно.

Анакрит. А, в смысле, что вы имеете в виду бесконечную выборку, являющуюся генеральной совокупностью. Тогда конечно. Но не любая бесконечная выборка — генеральная совокупность.

Диоклез. Вы согласились, что бесконечная выборка — это генеральная совокупность, мол, конечно, и тут же себя отрицаете. Вот правда: вся математика — это генеральная совокупность. Если вы бесконечное число раз измерите бесконечное число вселенных, то тогда найдёте истинные параметры. И через бесконечное число лет, как я уже говорил, произойдёт реализация открытых законов Вселенной, которые в неё изначально заложены... Нет, не заложены, просто из которых она состоит.

Анакрит. Генеральная совокупность может быть бесконечной выборкой, но не всякая бесконечная выборка — генеральная совокупность.

Лучше ответьте мне на такой вопрос. Монета выпадает решкой с вероятностью 0,5. Вы кидаете бесконечное число раз, и всё время выпадает решка. Вполне возможная ситуация?

Диоклез. Да, значит, во Вселенной есть какая-то хитрая дуальность, которая генерализует бесконечность на более широкий класс бесконечностей.

Анакрит. Удивительно... Это самое удивительное в пифагорейцах — в упор не видеть, что нет наук враждебнее, чем физика и математика, и упорно скреплять их мистикой, не имеющей отношения ни к той, ни к другой науке.

Диоклез. В пифагорейцах самое враждебное — это не ощипывать белого петуха через перекладину. (Здесь Диоклез высмеивает религиозные верования пифагорейцев. — Прим. авт.)

И в этой бесконечности у нас бесконечность решек, в другой — бесконечность орлов. И все их линейные комбинации. Под словом «бесконечность» я имею в виду бесконечно итерированную и генерализованную совокупность бесконечностей. Просто по определению монета — это 0,5 орла и 0,5 решки. Это определение справедливой монеты. И в этой самой бесконечной бесконечности, предельно полной, то есть в этой всей математике, всё, что распадается на две половины этой математики, является монетками.

Анакрит. Да, я понял. Вы об удивительном довольно явлении — выборке, тождественной генеральной совокупности. Стремление Души к Нусу (по Плотину).

Диоклез. То есть если есть прямая, являющаяся монеткой, то два луча — это орёл и решка. И понятно, что решка — это одни решки, то есть \(A=A\). Допустим, весь правый луч — решки, весь левый — орлы. Мы смотрим на прямую с левого торца и видим одних орлов в точках. Но дело в том, что эта бесконечность неполна. Это другой порядок малости. Аль забыли, что «о» малое бывает малых порядков? Так вот, им обратные величины бывают разных бесконечных порядков. Предельно бесконечный порядок — это и есть монетка, состоящая из двух половин. Это её определение.

Анакрит. Правда, эта чудо-выборка перестаёт быть выборкой. Но вы называете её выборкой потому, что вам приятно думать, что ваши жалкие земные конечные выборки имеют с ней некую общую природу, отличаясь лишь количественно. Вы отличие конечного и бесконечного мыслите количественным, не качественным.

Диоклез. Всё имеет единую природу бесконечной Математики, которая аппроксимируема и реализуема в бесконенчности!

Анакрит. Да, это я понял. Бесконечность бесконечностей.

Диоклез. Поэтому качественные и количественные различия — это одно и то же. Некорректно дифференцировать качественные и количественные изменения. Если величина меняется, то это уже изменение. Если из существующего набора возникает что-то другое, то это изменение как количественное, так и качественное. Даже репродукция объектов — это заполнение всей бесконечности пустотой; вот смысл размножения! То есть исчезновение, появление или изменение — это всё разные стороны одного и того же явления, именуемого нетождественностью.

Анакрит. Я бы дал такое определение: различие между двумя величинами является количественным, если одну можно получить из другой через конечное число арифметический действий. Тогда я прав, но это очевидно подогнанное определение.

Диоклез. Грубо. И определение у вас никчёмное. Мир бесконечно делим на бесконечно малые части, поэтому любое изменение — это уже бесконечность изменений. Вам нечем ответить!

Анакрит. Вы отрицаете качественные различия, я даю им определения в рамках признаваемых вами понятий. Спор терминологический, а не сущностный. Зачем вы пересказали апорию Зенона? Она тут при чём? «Мир бесконечно делим...» Ещё хуже, не к месту переоткрываете велосипед!

Диоклез. Если у нас есть бесконечность, то можно на неё поделить. Качественных различий не бывает. И количественных тоже не бывает. Есть просто тождественные состояния вселенной (пять яблок в этой и пять яблок в той). А есть нетождественные (пять и четыре соответственно). Индикатор различия есть единственно допустимый критерий оценки величины. С вашей точки зрения различие между пятью яблоками и четырьмя яблоками количественное. Между двумя яблоками и двумя грушами — качественное. Вопрос номер один: какова природа различий между двумя яблоками и одной грушей? Вопрос второй: какова природа различий между пятью яблоками и пятью яблоками? Качественные или количественные? Пока подумайте над этими вопросами.

Я вот о чём думаю. Если есть действие — написать письмо — то должно быть антидействие. Нет, не «не написать». И не отмотать время назад. Что-то, совсем проделывающее инверсию... Я это чувствую, но даже не знаю, как это описать. Сесть на стул. Обратное действие — это... как бы сотворить в одновременно существующей параллельной вселенной такое действие, чтобы их суммирование приводило к тому, что в результирующей вселенной никто никуда не садился. Сесть антиматерийно в обратной инверсии со всеми отрицательными показателями.

Анакрит. Теперь домножьте действие «сесть на стул» на \(i\) и опишите.

Диоклез. Понимаете, это сложнее. Обратное действие — это что-то, что очень естественно. Даже простое сложение, то бишь суперимпозицию двух посажений, надо сначала ввести. Смотрите. Я знаю, как получить ответ на ваш вопрос. Надо сначала ввести сложение во вселенной, потом получается отрицание, потом умножение (тоже через обратное действие), потенциирование. Главное правило — чтобы в сумме был ноль. Поэтому сложение, умножение и потенциирование тождественны тем, что определены через ноль.

Нет разницы между сложением, умножением на минус один и взятием корня. Любое действие во вселенной — это такое явление, суперимпозиция которого на своё противоявление даёт ноль. Под суперимпозицией я имею в виду наложение, одновременное сочетание, слияние воедино.

Анакрит. Вы мало того что признаёте альтернативные вселенные, так ещё и считаете нашу подчинённой некоей их совокупности, а не образующей её...

Диоклез. Понимаете, мы, как частный случай подчинения Великому Единому Закону Нуля, являемся только какой-то реализацией, которая за бесконечное число времени в бесконечном количестве ипостасей станет генеральной совокупностью.

Анакрит. Вы читали «Гамлета»? А «Розенкранц и Гильденстерн мертвы»?

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

Анакрит. Ох, ну вы так свободно переходите от математики к физике! Я так не умею.

Диоклез. Потому что всё возможно! И всё противоположное возможно! А в сумме — ноль.

Анакрит. Для меня ведь это бессмысленно.

Диоклез. Всё — это положительная часть Вселенной. Всё противоположное — это тёмная материя, которую мы не можем засечь. В сумме — ноль.

Анакрит. Относительно физики, то есть Вселенной, вы всё верно сказали, тут есть тождественность и нетождественность. А в математике количественно различие между одним и двумя, качественно между одним и бесконечностью.

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

Анакрит. Вы помните моё доказательство небытия Бога? Там утверждается, что понятия «существующее», «возможное» и «необходимое» тождественны (и главным является «существующее»), как и их антонимы.

Диоклез. Помню. Совершенно нормальная аксиоматика.

Анакрит. Это аксиоматика физикалистическая. Противоположная вашей.

Диоклез. Естественно. Есть материя, антиматерия, их линейные комбинации со множителями комплексными, кватернионными — какими угодно. Главное у меня вот что: ноль содержит сам себя. Это очень важно. Дело в том, что моя идеальная Бесконечность, то есть Математика-Ноль, всё это содержит. Половина решек (материя) и орлов (антиматерия) могут взяться в любой дикой комбинации из настолько генерализованных чисел, что представить нельзя, то есть чисел, генерализованных степени бесконечности Вселенной. Она содержит всё и саму себя. Это очень важно и нетривиально! Запомните, она сама содержит себя, антисебя, а также все их линейные и нелинейные комбинации — при этом всё это есть она сама. Если быть точнее, то линейные комбинации с коэффициентами, определёнными в бесконечной степени общо, степени общности Вселенной-Математики-Нуля. Любое осознание — ничто по сравнению с тем фактом, что ноль содержит ноль, при этом оставаясь ничего не содержащим нулём!

Анакрит. Возможна ли вселенная, в которой хоть один атом отличен от атома нашей? Очевидно, из линейной комбинации атомов нашей вселенной получить, вообразить такую можно. Так вот, эта вселенная, которую мы вообразим, — она возможна?

Диоклез. Возможна. Но мы не можем её получить средствами, доступными в нашей вселенной.

Анакрит. Да, это понятно. Вопрос второй: эта гипотетическая вселенная существует?

Диоклез. Да. Непременно существует.

Анакрит. Вопрос третий. Эта вселенная не просто есть, но обязана быть? То есть это не случайность, но фатальность, абсолютный факт?

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

Анакрит. Понимаю! Блестящее описание, кстати говоря. Прежде чем перейти дальше, один вопрос, который для меня является техническим, для вас может показаться издевательским, но, умоляю, задумайтесь о нём всерьёз на миг. Не может ли быть такого, что вы правы во всём, кроме того, что всё сводится к нулю? Самое интересное, что у вас, как и у меня, тождественны существование, возможность и необходимость. Но у меня главенствует существование, у вас — возможность. Смысл выходит противоположный. Существует ли вселенная, в которой утверждаемое вами неверно?

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

Анакрит. Верно ли, что, какое бы утверждение я ни написал в кавычках, существует бесчисленное множество вселенных, в которых оно будет истинным?

Диоклез. Да. И бесконечное множество вселенных, в которых будет не только истинным, но и не ложным, и не промежуточным, а истинное отношение будет задаваться оператором, который всё итерирует во всех направлениях бесконечности. Но и так понятно, что сумма всех направлений — это ноль.

Даже то, что есть плюс и минус, операция и обратная ей, — это несправедливо. Это как шахматы на двоих. Должен быть плюс, минус, что-то среднее, что-то полуторное, минус среднее и так далее. Понимаете, всё это будет иметь смысл только при бесконечном обобщении. А в бесконечном обобщении существует бесконечное обобщение. Я пока не могу придумать шахматы с 2,5 цветами на 2,5 человек, но на двоих (частный случай) сыграть можно.

Вот вам пара вопросов на подумать. Раз. Что такое две с половиной точки на плоскости? Два. Что такое минус точка? Чем она отличается от минус куба? Вопрос со звёздочкой: плоскость с точкой, принадлежащей этой плоскости, — это больше, чем плоскость, или столько же? Является ли точка частным случаем нуля как ничего? А сумма точки и минус точки?

Анакрит. Две точки на плоскости — прямая. Три — вся плоскость. Две с половиной образуют нечто между ними. Это полуплоскость — плоскость, ограниченная с одного конца. На ней две обычные точки и одна, лежащая на прямой, за которой бездна. А минус точка — это то, что относится к минус кубу так же, как точка к кубу. Большего сказать не могу.

Диоклез. Первый ответ неубедителен. Надо сначала специфицировать понятие «часть точки», а потом уже отвечать. 2,5 точки — это количество точек, размерность Хаусдорфа линейной оболочки (евклидовой) которых равна 1,5.

Анакрит. В последнем вопросе я, подобно Лейбницу, признаю реальность бесконечно малых. Точка и плоскость — это больше, чем плоскость, точка больше, чем ноль. Сумма точки и минус точки — не то ноль, не то пустое пространство. Скорее, пустое пространство. С размерностью той из точек, размерность которых была больше. Я не хотел оперировать понятием дробной размерности, поскольку не понимаю его, и попытался переизобрести его.

Диоклез. Постойте, плоскость — это множество точек. Это бесконечность размерности два. Точка \(A\) уже есть в этой плоскости. Поэтому плоскость — это \(\infty^2\) точек, просто одну из которых мы назвали и выделили.

Анакрит. Позвольте! Про «часть точки» я ничего не говорил!

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

Анакрит. А это смотря какую точку прибавляем! Я говорил о случае, когда мы складываем плоскость и точку, не лежащую на этой плоскости (то ли просто лежащую в другой точке этого n-мерного пространства, то ли вообще из другой вселенной). Тогда получается нечто большее.

Диоклез. Тогда, конечно, верно.

Анакрит. Размерность Хаусдорфа согласуется с нашими обычными представлениями о размерности в тех случаях, когда эти обычные представления есть.

Диоклез. Сумма точки и минус точки — это чистый арифметический нульмерный ноль. Сумма прямой и минус прямой — это чистый арифметический одномерный ноль. Поразмышляйте перед сном о разнице между нульмерным и одномерным нулями, и станете причастны к Моей аксиоматике.

Анакрит. Арифметический? В арифметике есть измерения? О, да вы же о координатах!

Диоклез. Я ошибся. Надо будет переспецифицировать нули. Арифметический, геометрический одномерный, нульмерный... Надо столько нулей придумать! Будет забавно, скажем, если у меня получится 17 видов нулей. Должно же быть ноль. Значит, в другой вселенной −17 нулей, всё остальное — линейная комбинация.

Анакрит. Подумайте о метавселенных побольше. Вот вам подарок:

Послание Единице

Всяк безумец, что ищет измены,
На тебя домножать будет зря:
Так, арабы изрезали вены,
Мысля корень из минус тебя.

Или в степень твою возведенье;
Ну да что; в степень степени (sic!):
Да и то — чёрт, а не изменение,
Даже не измененье, а пшик.

Небо яростно я вопрошаю,
Мозга ткани дымится покрой:
Не сказать то, что ты составная;
Как же можешь не быть ты простой?

Бог один, но различные лица,
Иль три бога, и каждый един?
\(i\) в квадрате, видать, единица,
Раз один в минус первой — один.

Ты противен мне, кол похоронный;
О, суккуб! Я дымлюсь, будто торф!
Раз система уж позиционна,
Единица — всегда автоморф.

Диоклез. Это вы?

Анакрит. Ага, под антисократический диалог. Надоело, не умею доводить стихотворения до конца. Не говоря о прочем; надо бы впихнуть стрелки Грэма и остальные свойства числа, да лень. Я склонен считать ноль натуральным числом. Хотя есть аргументы и против того. А вы?

Диоклез. Вопрос — что такое натуральное число?

Анакрит. Так вопрос как раз в том, какое определение принять. Если мы постепенно расширяем множества: пустое, ноль, натуральные, целые, рациональные, действительные — то с нулём. Если вспоминаем истоки понятия, то в счёте нуля нет, в самом деле.

Диоклез. В счёте ноль есть. Ноль яблок.

Анакрит. Греки единицу-то не всегда числом читали, что уж говорить о нуле.

Диоклез. Греки вообще любили извращения, и белокурыми рабами они не брезговали.

Анакрит. Ах, то есть вы математиколюб, но не признаёте греков. А кого? Немцев? Или хотя бы французов?

Диоклез. Греки вообще слишком ортодоксальными были. Если бы Вейерштрасс порядка не навёл, математика бы стала ещё одной социологией. В евклидовом пространстве настолько всё банально и небесконечно!

Анакрит. Конструктивистскую математику не любите, любите актуальную бесконечность?

Диоклез. Видимо, да.

Анакрит. Омерзительно! Не прельщают вас точные и изящные формулировки. Сравните однозначное «для каждого простого числа \(P\) существует простое число, большее \(P\)» — и невнятный лепет «множество простых чисел бесконечно». Абсурд!

Диоклез. Мало того, у этого множества есть мощность!

Анакрит. Пф-ф. Лженаука, ересь, мистика. Никакой стройности.

Диоклез. Наоборот, я за вейерштрассову стройность. Всё можно параметризовать, вычислить контрольную сумму.

Анакрит. Немцы слишком техничны, их системам недостаёт единства, удобства и понятности. При чём тут параметризация? Да и основа зыбкая. Аксиоматика должна строиться на реальном. А актуальная бесконечность иллюзорна.

Диоклез. При том, что у бесконечностей должна быть шкала абсолютной мощности. \(\aleph_0\), \(\aleph_1\) и так далее.

Анакрит. Я думал об абсолютной мощности года 2 назад...

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

Анакрит. Купил и хотел подарить вам книгу «Гениальность и помешательство» Ломброзо, да пообещал другому. Для вас, впрочем, другая вещь есть.

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

Анакрит. «Я не гений» — весьма провокационное, эгоцентричное и надменное высказывание. Даже от вас подобного дерзновения не ждал. У вас, кроме материи и антиматерии, есть другие возможные вселенные и породившая их и нашу математика.

Диоклез. Почему это дерзновение?!

Анакрит. Есть интеллектуальная последовательность: животное (вне интеллекта), посредственность (нулевой интеллект), рационалист (конечный интеллект), талант (стремятся в бесконечность, но не достигающий её интеллект) и гений — достигший бесконечности и обратившийся в безумие животного интеллект. Вы не животное: интеллект может быть оценён; не посредственность: обладаете им; не рационалист: явная претензия на контакт с бесконечностью. Казалось бы, талант — но талант стремится к гению; вы утверждаете, что, напротив, от бесконечности скрываетесь, не являясь потому талантом, как и гением. Претензия на бытие-над-схваткой: будто бы вас нельзя классифицировать, вы свободно регулируете собственный разум, становясь кем угодно по личной прихоти. Претензия на нечто большее, чем гений, безумец и даже животное. Фантастическая наглость!

Диоклез. А что такое «интеллект, обратившийся в безумие животного»?

Анакрит. Это сошедший с ума гений.

Диоклез. А почему не просто в безумие, а в безумие животного?

Анакрит. Потому что все животные безумны.

Диоклез. То есть бесконечность, ставшая нулём?

Анакрит. Нет. Ноль — это посредственность, животное — вне чисел. Переменная, достигшая бесконечности, а затем ставшая пустым множеством. Талант есть обещание гения. И талант во времени является отрезком; можно сколь угодно долго быть талантом. А гений моментален, гений становится таковым в некоторый миг, в предыдущий и последующий гением не являясь. Собственно, когда талант становится гением, — это не сила, а срыв, прекращение. Вы утверждаете, что способны быть просто талантом, что гениальность вам не грозит. Иммунитет.

Диоклез. Она грозит. Просто я это отрицаю, чтобы не лишиться рассудка и не выродиться в null space.

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

Анакрит, не в силах продолжить мысль после истощившей его беседы, под звонкий смех бородатого Диоклеза засыпает.