Обнаружение критических точек на фондовом рынке

. Это делается с помощью матрицы векторов сезонности для каждого значения t из исторических и будущих данных. Например, для N = 10 , а сезонная компонента имеет вид . При этом , а увеличение N (порядка ряда Фурье) позволяет поймать более часто меняющиеся сезонные схемы, хоть и повышает риск переобучения.

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

Рассмотренный алгоритм в паре с оффлайн CPD-методами применяется для решения первой из четырех задач, поставленных во введении.

2.3 Онлайн CPD

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

, тревога (alarm), если , где - параметр скольжения (drift), - порог (threshold), - разница между последними двумя наблюдаемыми значениями. Другими словами, алгоритм считает сумму положительных изменений в наблюдаемой последовательности и сравнивает её с порогом, при превышении которого выдается сигнал тревоги, а сумма обнуляется. Параметр скольжения нужен для того, чтобы избежать сигналов тревоги при отсутствии реальных изменений. Само уравнение носило смысл так называемого «фильтра низких частот». Первоначальная форма CUSUM-теста реагировала только на изменения в среднем наблюдаемой последовательности в большую сторону, однако позднее появились различные расширения этого алгоритма, учитывающие изменения среднего в обе стороны, изменения дисперсии, хранящие более длинную историю наблюдений для и т. д.

В данной работе применяется следующая версия CUSUM-алгоритма:

Приведенный алгоритм зависит от параметров drift и threshold, которые необходимо правильно подобрать. В книге «Adaptive Filtering and Change Detection» F. Gustafsson «Adaptive Filtering and Change Detection», 2000 предлагается следующая процедура подбора этих параметров:

Начать следует с очень большого порога;

Параметр скольжения следует взять равным половине ожидаемого изменения, либо подобрать его таким, чтобы g = 0 выполнялось более чем в половине случаев;

Затем нужно установить порог таким, чтобы достигнуть требуемого количества ложных сигналов тревоги или требуемой задержки определения change point`а;

Если критическую точку необходимо определять быстрее, следует уменьшить параметр скольжения;

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

Если наблюдается множество критических точек, не связанных с реальными изменениями, следует увеличить параметр скольжения.

Помимо онлайн определения change point`ов, используемый CUSUM-алгоритм включает в себя ретроспективную процедуру вычисления начальной и конечной точки изменения, так что можно наглядно увидеть процесс изменения и момент, в который алгоритм дал сигнал тревоги.

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

Предполагается, что можно вычислить прогнозное распределение длины текущего сегмента с учетом нового наблюдения, опираясь на апостериорное распределение текущего наблюдаемого сегмента: . Апостериорное распределение находится рекурсивно из совместного распределения длины текущего сегмента и известных данных:



.

Стоит заметить, что прогнозное распределение зависит только от текущих данных . Таким образом, можно сгенерировать так называемый рекурсивный «алгоритм передачи сообщений» для совместного распределения длины текущего сегмента и известных данных, опираясь на два вычисления:

Априорное распределение при известном ;

Предиктивное распределение с учетом нового наблюдения при известных наблюдениях со времен последней критической точки.

Условное априорное распределение принимает ненулевые значения только в двух случаях: если длина текущего сегмента продолжает увеличиваться и , либо если возникает критическая точка и

. То есть

,

где - так называемая функция опасности (hazard function), определяемая как

.

В случае, если - дискретное экспоненциальное распределение (оно же геометрическое) с параметром , функция опасности принимает константный вид .

Алгоритм должен определять не только рекуррентное соотношение, но и так называемые «условия инициализации». Рассматривается два случая:

Критическая точка возникла до первого наблюдения .

Наблюдается недавнее подмножество данных, как, например, в случае мониторинга климатических изменений. Тогда априорная функция распределения вероятностей первоначального сегмента определяется так называемой нормализованной «функцией выживания» (survival function):

,

где Z - подходящая нормализирующая константа, .

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

2.4 LSTM нейросеть

Для задачи краткосрочного прогнозирования в данном исследовании используется LSTM нейросеть. Этот алгоритм является представителем класса рекуррентных нейронных сетей, главное отличие которых от обычных нейросетей заключается в том, что они «имеют память». Формально, они просто сохраняют информацию на каждом шаге и передают её на следующий. Более четко этот процесс можно увидеть на рисунке 6:

Рис. 6. Развернутая рекуррентная нейросеть.

Источник: http://colah.github.io/posts/2015-08-Understanding-LSTMs/

Здесь A - повторяющийся модуль, - входные данные, - выходные данные. Такое свойство делает алгоритм теоретически подходящим для краткосрочных предсказаний финансовых рядов.

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

Рис. 7. Структура повторяющегося модуля LSTM.

Источник: http://colah.github.io/posts/2015-08-Understanding-LSTMs/

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

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

[0, 1] с помощью сигмоид-функции , тем самым определяя, какую долю информации для каждого компонента необходимо запомнить. Как видно на рисунке 7, повторяющийся модуль включает три штуки таких «ворот» для защиты и контроля конвейерной ленты. Данные проходят через все слои по порядку слева направо.

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

Рис. 8. Процесс «забывания» ненужной информации.

Следующий шаг заключается в том, чтобы решить, какую новую информацию следует добавить на «конвейер». Этот процесс определяется сразу двумя следующими по порядку слоями. Сначала сигмоидный слой, называющийся «input gate layer», решает, какие значения следует обновить. Затем тангенс-слой (tanh) создает вектор потенциальных кандидатов , после чего выходные значения обоих слоев будут объединены для обновления состояния узла. В случае рассмотренного выше примера - необходимо добавить род нового подлежащего взамен забытого старого. Данный шаг можно увидеть на рисунке 9:

Рис. 9. Подготовка к обновлению состояния узла.

Источник: http://colah.github.io/posts/2015-08-Understanding-LSTMs/

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



Рис. 10. Процесс обновления состояния узла.

Источник: http://colah.github.io/posts/2015-08-Understanding-LSTMs/

Заключительный шаг - определить, какую информацию отправить «на выход». Эта информация основана на текущем состоянии узла, однако её еще необходимо отфильтровать. Сначала сигмоидный слой определяет, какие части текущей информации на конвейере необходимо отправить дальше. Затем данные проходят через тангенс-слой, который преобразует их в формат отрезка [-1, 1] с помощью гиперболической тангенс-функции . Наконец, результаты обоих слоев перемножаются между собой и дальше отправляется только необходимая информация. Возвращаясь к примеру с предсказанием слова, если алгоритм получил данные о подлежащем, он может решить отправить дальше информацию, влияющую на форму глагола (например, число), если следующее слово ожидается глаголом. Заключительный шаг представлен на рисунке 11:

Рис. 11. Заключительный этап повторяющегося модуля.

Источник: http://colah.github.io/posts/2015-08-Understanding-LSTMs/

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

2.5 Математическая модель броуновского движения

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

,

где S - цена акции, - коэффициент скольжения, - коэффициент рассеивания,

- броуновское движение. Под броуновским движением в данном случае понимается его математическая модель, в теории случайных процессов называемая винеровским процессом https://ru.wikipedia.org/wiki/Винеровский_процесс. Рассуждая в терминах прибыли , коэффициент скольжения отражает среднюю прибыль за какой-то период времени, а коэффициент рассеивания - стандартное отклонение этой же величины. Стохастическая часть представленного выше уравнения заключена в . Каждое броуновское приращение считается путем умножения стандартной нормальной случайной величины на квадратный корень приращения времени , т. е. . Кумулятивная сумма броуновских приращений является так называемым дискретизированным броуновским путем . Для рассмотренного стохастического дифференциального уравнения движение цены акции с начальной ценой можно описать моделью геометрического броуновского движения . В данном исследовании такая модель используется для подсчета вероятности достижения определенного значения цены до указанного момента времени.

Рассмотрим стандартный винеровский процесс , начинающийся в точке , и некоторое значение . Момент времени - точка, в которой впервые процесс принял значение , т. е.

. Некоторые из таких процессов достигнут довольно быстро, другие уйдут в отрицательные значения и достигнут только через какое-то долгое время. Поэтому необходимо определить распределение вероятностей для . При этом такое распределение будем искать при условии, что для какого-то заранее определенного t:

.

Второе слагаемое в этом уравнении равно 0, так как оно отражает «пустое» событие. Таким образом, . Теперь представим, что после момента наш процесс «начался заново», тогда, согласно принципу отражения винеровских процессов https://en.wikipedia.org/wiki/Reflection_principle_(Wiener_process)

, или

.

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

Здесь важно заметить, что на самом деле - случайное, а не фиксированное время, оно различается для каждого сгенерированного процесса. Однако, при определенных условиях, случайное время может вести себя так же, как и фиксированное. В этом случае это обуславливается тем, что стандартный винеровский процесс обладает сильными марковскими свойствами https://en.wikipedia.org/wiki/Markov_property (то есть условное распределение вероятностей будущих состояний процесса зависит только от настоящего момента времени и не зависит от всего что ему предшествовало), а - марковское время, соответствующее событию первого достижения значения .

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

Глава 3. Результаты исследования

3.1 Предсказывание критических точек с помощью комбинации оффлайн CPD и алгоритма Prophet

Представим следующую ситуацию: некий агент наблюдает колебания индекса в течение какого-то относительно долгого периода времени (в интервале от года до двух), после чего решается открыть позицию, опираясь на статистический прогноз. Но что, если впереди неожиданный кризис? Сможет ли алгоритм предсказать правильное направление тренда, чтобы можно было получить прибыль, открыв короткую позицию? Чтобы ответить на этот вопрос, попробуем предсказать три известных случая обвала фондовых рынков: крах индекса Moex 2008 года и крахи индекса DJIA 1929 и 1987 годов.

Процедуру настраивания алгоритма прогноза проведем в два этапа: сначала определим лучшие параметры модели Prophet с помощью процедуры поиска по сетке (grid search); затем добавим в модель критические точки, которые найдет оффлайн CPD-алгоритм. Заметим, что в Prophet уже имеется встроенная процедура поиска change point`ов, которая по умолчанию определяет 25 критических точек на заданной последовательности. Поэтому на втором этапе решение о добавлении «своих» change point`ов принимается исходя из сравнения качества моделей с разными вариантами критических точек.

Начнем с описания настраиваемых параметров модели:

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

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

Порядок ряда Фурье - чувствительность модели к сезонным изменениям (чем выше, тем чувствительнее). Рассматриваются варианты fo = {4, 6, 8, 10, 12};

Длина сезонного цикла. Рассматриваются варианты cycle_length =

{50, 100, 150, 200, 250, 300, 350};

Моделирование методом Монте-Карло по схеме марковской цепи (MCMC) с заданным количеством выборок, либо методом оценки апостериорного максимума https://ru.wikipedia.org/wiki/Оценка_апостериорного_максимума (MAP) - по сути дополненный метод максимального правдоподобия. Первый вариант показал себя хуже, поскольку тяжело подобрать оптимальное количество выборок, так что используется MAP-оценивание.

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

Сначала рассмотрим ситуацию краха индекса Moex 2008 года. В качестве обучающей выборки взят период с 9 января 2007 года по 1 июля 2008 года, всего 368 торговых дней. Горизонт прогнозирования - до 30 апреля 2009 года, то есть следующие 209 торговых дней. В таблице 1 можно увидеть топ-5 моделей с линейным трендом. В первом столбце находятся значения параметров, отвечающие за модель сезонности, порядок ряда Фурье и длину сезонного цикла соответственно. Во втором - ошибка прогноза.

Таблица 1. Источник: расчеты автора.

Trend = Linear
Model params	MAE
['multiplicative', 4, 250]	244.267
['multiplicative', 12, 250]	252.0496
['multiplicative', 6, 250]	252.4453
['additive', 10, 250]	290.8859
['additive', 12, 250]	308.8785

Несмотря на то, что наилучший результат показал набор параметров [`multiplicative', 4, 250], для дальнейшей работы используется вторая по качеству модель. Это сделано по трем причинам: во-первых, 4-й порядок рядов делает кривую слишком гладкой; во-вторых, разница показателей качества несущественна; в-третьих, модель [`multiplicative', 12, 250] лучше всего (из линейных) показала себя на «докризисном» прогнозе, который можно увидеть на рисунке 12:

Рис. 12. Докризисный прогноз модели линейного тренда.

Источник: расчеты автора.

Тренировочная выборка здесь - период с 9 января 2007 по 28 декабря 2007, горизонт прогнозирования - до 1 августа 2008. Зеленой звездочкой обозначена дата «входа на рынок», что в данном случае является последней датой тренировочной выборки. Как видно на графике, обучающие данные очень похожие на данные с явной сезонностью, чем и объясняется такой прогноз модели. Однако, модель хорошо ловит «первичный» убывающий тренд, что, по сути, и является целью. Полагаться на точные долгосрочные прогнозы фондового рынка все равно нельзя, так что дальнейшие ошибки прогноза для рассматриваемой задачи не так важны, если «первичный» тренд предсказан правильно.

Теперь проверим эту модель на «кризисном» прогнозе. Локальные экстремумы прогноза будем считать критическими точками, так как они сигнализируют об изменении направления тренда. Результат можно увидеть на рисунке 13:

Рис. 13. Прогноз обвала индекса Moex 2008 года.

Источник: расчеты автора.

Как видим, тренд предсказан очень неплохо. Если бы агент открыл короткую позицию, опираясь на этот прогноз, и закрывал её в локальных минимумах, его прибыль могла бы составить 11.3%, 62.9% или 62.4% соответственно. Ранее упоминалось, что доверять долгосрочному прогнозу нельзя, но даже с учетом рискофобии агента, он остается в прибыли.

Попробуем еще улучшить качество модели, добавив в неё критические точки, определенные рассмотренными в главе 2 оффлайн CPD-алгоритмами. Средняя абсолютная ошибка модели с автоматически определенными Prophet change point`ами составляет 252.0496, попытаемся улучшить этот показатель.

В следующих таблицах приведены результаты добавления в модель критических точек, определенных всеми возможных комбинациями рассмотренных cost и search-функций.

Таблица 2. Источник: расчеты автора.

Pelt (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	50	68	250.6947
normal	5	44	249.4407
ar	5500	14	218.3051
mahalanobis	0.006	49	250.6947
rbf	0.1	56	248.4888

Dynp (pre-set n_chpts)
cost	n_chpts	MAE
l2	65	255.6115
normal	71	254.4073
ar	14	218.3051
mahalanobis	65	255.6115
rbf	65	248.9117

Таблица 3. Источник: расчеты автора.

Таблица 4. Источник: расчеты автора.

Window (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	0	3	881.2166
normal	0	3	885.784
ar	0	3	885.784
mahalanobis	0	3	881.2166
rbf	0	3	885.784

Таблица 5. Источник: расчеты автора.

Binseg (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	50	68	250.6947
normal	6	42	235.1821
ar	4000	16	231.5926
mahalanobis	0.006	68	250.6947
rbf	0.1	56	248.4888

Таблица 6. Источник: расчеты автора.

BotUp (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	50	68	250.6947
normal	5	44	249.4407
ar	3500	31	257.7799
mahalanobis	0.006	68	250.6947
rbf	0.06	65	255.1771

Параметры для каждого случая были подобраны вручную. Алгоритм может очень чувствительно реагировать на минимальную разницу в количестве и расположении критических точек, при этом вариантов различных наборов параметров для использованных CPD-процедур огромное множество, так что представленные результаты могут быть не оптимальными.

Стоит заметить, что одни и те же cost-функции могут выдавать одинаковый результат в различных алгоритмах, но это работает не всегда (даже при идентичных параметрах два «оптимальных» алгоритма могут выдавать разные результаты). Результаты работы функций издержек и (mahalanobis) для всех search-функций одинаковы, так как без указания спецификации матрицы M они ничем не отличаются. Процедура «Window-sliding» оказалась для рассматриваемых данных худшей, для всех остальных же удалось подобрать параметры, улучшающие изначальный показатель MAE. Лучшей комбинацией оказалась авторегрессионная функция издержек для 14 change point`ов при одном из оптимальных алгоритмов (Pelt или Dynp он же Opt). Распределение этих критических точек проиллюстрировано на рисунке 14.

Рис. 14. Наилучшие критические точки для тренировочной выборки краха Moex 2008.

Источник: расчеты автора.

Заменим определенные Prophet`ом change point`ы на найденные:

Рис. 15. Скорректированный прогноз обвала Moex 2008. Источник: расчеты автора.

Как видим, кривая местами стала ближе к истинным значениям и пропал небольшой скачек вверх в начале прогноза, однако глобально ничего не поменялось. Один из оставшихся локальных минимумов все же немного сдвинулся, теперь возможная прибыль могла бы составлять 64.5% либо 62.4%.

Посмотрим теперь как покажет себя модель логистического тренда. Рассмотрим следующие значения пропускной способности для известных и будущих данных соответственно: 2000 и 2000, так как выше этой отметки значения не поднимались; 2500 и 2000 в надежде поймать убывающий в будущем тренд; 2000 и 1600 чтобы получить возможно хорошую, но пессимистичную на момент входа на рынок модель. Топ-5 кандидатов для каждого случая можно увидеть в таблицах ниже.

Таблица 7. Источник: расчеты автора.

Trend = Logistic, cap=(2000, 2000)
Model params	MAE
['additive', 12, 200]	868.5953
['multiplicative', 12, 200]	880.5548
['multiplicative', 10, 200]	881.3598
['additive', 10, 200]	881.8598
['multiplicative', 8, 200]	882.3476

Таблица 8. Источник: расчеты автора.

Trend = Logistic, cap=(2500, 2000)
Model params	MAE
['multiplicative', 8, 200]	509.0105
['additive', 12, 200]	510.8004
['multiplicative', 10, 200]	511
['additive', 10, 200]	511.9039
['additive', 6, 200]	511.917

Таблица 9. Источник: расчеты автора.

Trend = Logistic, cap=(2000, 1600)
Model params	MAE
['additive', 12, 200]	553.1646
['additive', 8, 200]	560.5562
['additive', 4, 200]	560.7597
['additive', 6, 200]	560.8176
['additive', 10, 200]	561.0437

Как видим, модели логистического тренда подходят для рассматриваемой задачи гораздо хуже линейных. Лучшую из рассмотренных можно увидеть на рисунке 16:

Рис. 16. Прогноз кризиса Moex 2008 модели логистического тренда.

Источник: расчеты автора.

Модель не ловит «правильный» убывающий тренд. Для двух случаев краха DJIA будем рассматривать только линейные модели, однако стоит заметить, что хоть логистические модели и не подходят для «кризисных» прогнозов, на «докризисном» они показывают себя довольно неплохо. На рисунке 17 представлен «докризисный» прогноз модели логистического тренда с параметрами [`additive', 4, 200] и пропускной способностью 2500 и 2000 для известных и будущих данных соответственно. Критические точки для неё определены встроенной процедурой Prophet.



Рис. 17. Докризисный прогноз модели логистического тренда.

Источник: расчеты автора.

Рассмотрим теперь ситуацию краха индекса DJIA 1929 года, также известную как «Великую депрессию». В качестве обучающей выборки взят период я 1 июня 1928 года по 1 октября 1929 года, всего 400 торговых дней. Горизонт прогнозирования - до 31 декабря 1929 года, то есть следующие 68 торговых дней. В таблице 10 представлены топ-5 моделей с линейной функцией тренда:

Таблица 10. Источник: расчеты автора.

Trend = Linear
Model params	MAE
['multiplicative', 12, 300]	28.46617
['additive', 10, 300]	29.89035
['additive', 12, 300]	30.65023
['multiplicative', 10, 300]	32.12636
['multiplicative', 8, 300]	40.62104

Заметим, что здесь лучшая модель имеет те же параметры функции сезонности и порядка ряда Фурье, что и в предыдущем случае. Отличается лишь длина цикла, что логично, на разных рынках в разные времена длины циклов скорее всего должны отличаться. Далее сразу приведем различные варианты критических точек, после чего опишем лучшую модель.

Таблица 11. Источник: расчеты автора.

Pelt (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	25	48	27.97674
normal	12	28	27.01974
ar	100	16	27.861
mahalanobis	0.01	48	27.97674
rbf	0.01	59	27.32042

Таблица 12. Источник: расчеты автора.

Dynp (pre-set n_chpts)
cost	n_chpts	MAE
l2	41	28.00685
normal	28	26.9708
ar	16	27.861
mahalanobis	42	28.07607
rbf	59	27.32042

Таблица 13. Источник: расчеты автора.

Window (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	0	3	90.66405
normal	0	2	131.1072
ar	0	4	115.63
mahalanobis	0	3	90.66405
rbf	0	3	98.81838

Таблица 14. Источник: расчеты автора.

Binseg (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	22	51	28.15341
normal	12	29	26.79046
ar	50	18	28.97014
mahalanobis	0.01	48	28.16428
rbf	0.01	58	27.31028

Таблица 15. Источник: расчеты автора.

BotUp (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	25	48	27.97674
normal	12	29	26.95657
ar	50	33	28.25292
mahalanobis	0.01	48	27.97674
rbf	0.01	59	27.32042

Примечательно, что в данном случае лучшие результаты показал приближенный метод «Binseg» с функцией издержек (normal) для 29 критических точек (но это не значит, что он лучше определил сами положения точек, значит лишь что такое расположение лучше сказалось на прогнозе модели). Распределение этих точек можно увидеть на рисунке 18:

Рис. 18. Наилучшие критические точки для тренировочной выборки краха DJIA 1929.

Источник: расчеты автора.

Теперь посмотрим, как выглядит прогноз с учетом замены change point`ов, автоматически определенных встроенной процедурой Prophet, на найденные выше.

Рис. 19. Прогноз кризиса DJIA 1929 модели линейного тренда.

Источник: расчеты автора.

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

Последняя рассматриваемая в этом разделе ситуация - крах DJIA 1987 года, он же «Черный понедельник». Тренировочная выборка здесь - период со 2 января 1986 года по 1 октября 1987 года, всего 443 торговых дня. Горизонт прогнозирования - до 31 декабря 1987 года, то есть следующие 64 торговых дня. Топ-5 моделей с линейной функцией тренда представлены в таблице 16:

Таблица 16. Источник: расчеты автора.

Trend = Linear
Model params	MAE
['multiplicative', 6, 200]	479.1701
['multiplicative', 10, 200]	479.189
['multiplicative', 4, 200]	483.5446
['multiplicative', 8, 200]	485.5137
['multiplicative', 12, 200]	489.4003

Здесь, как и во всех предыдущих аналогичных таблицах, топ-5 моделей имеют одинаковую длину цикла. На самом деле, весь топ-n моделей как в этом, так и в каждом предыдущем рассмотренном случае занимают модели с различными комбинациями параметров для наилучшей длины цикла. Это говорит о том, что длина цикла - важнейший параметр модели для нашей задачи. Также отметим, что модель [`multiplicative', 12, 200], хорошо подошедшая для первых двух ситуаций, хоть и не является лучшей в данном случае, но не сильно проигрывает. Рассмотрим варианты критических точек:

Таблица 17. Источник: расчеты автора.

Pelt (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	700	56	479.4236
normal	15	31	494.2902
ar	1000	63	442.6816
mahalanobis	0	89	450.4048
rbf	0	88	443.3144

Таблица 18. Источник: расчеты автора.

Dynp (pre-set n_chpts)
cost	n_chpts	MAE
l2	56	481.7908
normal	31	494.2902
ar	63	442.6816
mahalanobis	88	443.0938
rbf	87	438.2808

Таблица 19. Источник: расчеты автора.

Window (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	0	3	670.3462
normal	0	3	765.1443
ar	0	2	811.3956
mahalanobis	0	3	670.3462
rbf	0	3	671.0094

Таблица 20. Источник: расчеты автора.

Binseg (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	700	56	479.4236
normal	8	39	499.638
ar	1000	33	456.8591
mahalanobis	0	89	450.4048
rbf	0	87	437.0308

Таблица 21. Источник: расчеты автора.

BotUp (n_chpts based on penalty)
cost	penalty	n_chpts	MAE
l2	700	56	479.4236
normal	11	39	498.646
ar	1500	50	454.0521
mahalanobis	0	89	450.4048
rbf	0	87	437.0308

В данном случае лучший результат у методов «Binseg» и «BotUp» для функции издержек и 87 критических точек. Их распределение можно увидеть на рисунке 20.

Рассмотрев различные варианты оффлайн CPD-алгоритмов для всех трех случаев, можно сделать следующие выводы:

Алгоритм «Window» для данной задачи не подходит;

Все остальные алгоритмы показывают приблизительно одинаковые результаты (в терминах качества прогноза Prophet);

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

Рис. 20. Наилучшие критические точки для тренировочной выборки краха DJIA 1987.

Прогноз в данном случае выглядит следующим образом:

Рис. 21. Прогноз кризиса DJIA 1987 модели линейного тренда.

Здесь, в отличие от предыдущих двух случаев, нельзя сказать, что факт кризиса был предсказан. Возможно, это не получилось потому, что обвал DJIA 1987 года был «слишком резким». Однако, убывающий тренд вновь угадан, а открытие короткой позиции и закрытие её в предсказанном локальном минимуме обернулось бы для агента прибылью в 26.3%.

Напоследок рассмотрим модель [`multiplicative', 12, 200] на данных кризиса DJIA 1987 года, которая хорошо подошла для первых двух случаев (но с другой длиной цикла). Как уже отмечалось, более низкий порядок ряда Фурье делает кривую более гладкой, что иногда может немного снизить ошибку прогноза (как в третьем случае), однако гладкие кривые не свойственны колебаниям на фондовом рынке. Критические точки для модели

[`multiplicative', 12, 200] оставим определенными встроенной процедурой, так как нас в данном случае интересует форма прогнозной кривой, а не минимизация средней абсолютной ошибки (которая все равно существенно не изменится). Результат можно увидеть на рисунке 22.

Рис. 22. Прогноз краха DJIA 1987 моделью с линейным трендом и параметрами [`multiplicative', 12, 200].

Источник: расчеты автора.

Как и ожидалось, тренд сохранился, а колебания стали более чувствительными. Закрывая короткую позицию в локальных минимумах, агент мог бы получить прибыль в размере 25.7%, 33% либо 24%. Можно сказать, что мультипликативная сезонность и порядок ряда Фурье, равный 12 подходит для всех трех рассмотренных случаев.

На этом завершается решение первой из четырех поставленных во введении задач.

3.2 Применение онлайн CPD-алгоритмов

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

Итак, рассмотрим ситуацию краха индекса Moex 2008 года. Допустим, мы ожидаем серьезных изменений, поэтому параметр порогового значения установим равным 400. Коэффициент скольжения в данном случае выбран равным 1. Результат можно у видеть на рисунке 23.

Рис. 23. Результат применения Online CUSUM (threshold=400, drift=1) для кризиса Moex 2008 года.

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

Теперь снизим пороговое значение в два раза: threshold = 200. Коэффициент скольжения для него подобран равным 5.

Рис. 24. Результат применения Online CUSUM (threshold=200, drift=5) для кризиса Moex 2008 года.

Источник: расчеты автора.

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

Перейдем к ситуации кризиса DJIA 1929 года. Первый набор параметров возьмем следующим: threshold = 100, drift = 0. Результаты на рисунке 25.

Рис. 25. Результат применения Online CUSUM (threshold=100, drift=0) для кризиса DJIA 1929 года.

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

Снизим пороговое значение до 50, коэффициент скольжение оставим равным нулю. Результат работы алгоритма в этом случае можно увидеть на рисунке 26.



Рис. 26. Результат применения Online CUSUM (threshold=50, drift=0) для кризиса DJIA 1929 года.

«Пре-кризисная» точка вновь осталась на своем месте, а предшествующее ей глобальное изменение теперь определяется раньше. Добавились еще несколько локальных точек, определяющих изменения либо около середины, либо чуть ближе к концу. Как и в предыдущей рассмотренной кризисной ситуации, результат применения онлайн CUSUM-алгоритма для кризиса DJIA 1929 года вполне можно считать удовлетворительным.

Наконец, рассмотрим ситуацию «Черного понедельника». Для первого случая выбраны параметры порогового значения, равного 400, и коэффициента скольжения, равного 0. Результат на рисунке 27.



Рис. 27. Результат применения Online CUSUM (threshold=400, drift=0) для кризиса DJIA 1987 года.

Источник: расчеты автора.

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

Вновь снизим параметр порогового значения в два раза: threshold = 200. Параметр скольжения для него возьмем равным 3. Результат можно увидеть на рисунке 28.



Рис. 28. Результат применения Online CUSUM (threshold=200, drift=3) для кризиса DJIA 1987 года.

Источник: расчеты автора.

Впервые «пре-кризисная» тока сдвинулась - теперь будущее падение определяется раньше. Предшествующие ей положительные изменения разбиты на два сегмента, первый из которых определен раньше середины, второй чуть позже. «Пост-кризисные» локальные изменения определены в момент, близкий или совпадающий с их завершением. В целом же результаты алгоритма и в этом случае скорее удовлетворительные, чем нет. Заметим также, что и при решении первой задачи исследования, и для «CUSUM-части» второй, данные о кризисе Moex 2008 и DJIA 1929 лучше «поддаются» используемым алгоритмам, нежели данные «Черного понедельника».

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

Будем полагать, что длина текущего сегмента имеет дискретное экспоненциальное (оно же геометрическое) распределение с параметром . Такое предположение взято из статьи Ryan P. Adams, David J. C. MacKay «Bayesian Online Changepoint Detection», 2007, где аналогичный алгоритм применялся к данным о колебаниях индекса DJIA за период 1972-1975 годов. Само значение настраивается отдельно для каждой рассмотренной выборки.

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

.

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

Таким образом, для случая краха Moex 2008 года, возьмем стартовые коэффициенты (также из статьи выше), (как ожидаемое среднее значение), а для остальных двух коэффициентов рассмотрим три случая: , . Эти значения получены путем выражения из уравнения для scale при предположении, что стандартное отклонение будет равно 200 и фиксированных . Параметр опытным путем установлен равным 100. Сразу отметим, что вероятность быть критической точкой для первого наблюдения последовательности всегда будет равна единице, чтобы длину первоначального сегмента начинать считать с нуля.

Каждый из описанных выше случаев рассмотрим по порядку. Результат применения алгоритма для первого из них можно увидеть на рисунке 29.



Рис. 29. Результат применения Bayesian Online CPD (, , ) для кризиса Moex 2008 года.

Источник: расчеты автора.

Красными точками здесь и далее отмечены значения, которые являются change point`ами с вероятностью . Представленное на графике распределение вероятностей считалось апостериорно для каждой точки, после появления одного нового наблюдения. Всего здесь 12 критических точек (включая первое наблюдение, которое, на самом деле, ею не является): с 9 по 12 января 2007 года, с 21 по 24 января 2008, 13 и 14 мая 2008, 25 июля 2008 и 19 сентября 2008. Подряд идущие точки, конечно, смущают, и, скорее всего, обусловлены шумом, а не реальными изменениями параметров распределения. Зато две последние определены довольно неплохо, причем первая из них очень близка к аналогичной «кризисной» точке, определенной CUSUM-алгоритмом (здесь она определена на день позже), а вторая - прямо в момент начала локального падения.

Посмотрим, как изменится ситуация для второго набора стартовых коэффициентов. Результат на рисунке 30:



Рис. 30. Результат применения Bayesian Online CPD (, , ) для кризиса Moex 2008 года.

Источник: расчеты автора.

Здесь 5 change point`ов не считая начальной точки последовательности:

11 января 2007, 18 июля 2007, 21 января 2008, 25 июля 2008 и 28 июля 2008. Подряд идущих точек теперь нет, лишь две последние расположены близко друг к другу. При этом «кризисная» точка осталась на своем месте, а вот для последующих данных алгоритм ничего не выявил, определив все последующие наблюдения в один сегмент.

Наконец, рассмотрим третий набор стартовых коэффициентов. Результаты можно увидеть на рисунке 31.



Рис. 31. Результат применения Bayesian Online CPD (, , ) для кризиса Moex 2008 года.

Источник: расчеты автора.

Как видим, результат, по сравнению с предыдущим случаем, изменился на одну точку. Теперь вместо 18 июля 2007 появилось 1 июля 2008.

Если рассматривать результаты с точки зрения предсказания кризиса, можно сказать, что все три протестированных набора параметров подходят. Однако, если учитывать реакцию алгоритма на более локальные изменения, лучше всего подходит первый набор параметров, хоть и содержит в себе достаточно большое количество false alarm`ов. Скорее всего это можно исправить, рассмотрев большое множество наборов стартовых коэффициентов, либо поменяв используемое априорное распределение наблюдений (например, на гамма-распределение), однако это выходит за рамки данного исследования.

Теперь применим алгоритм к данным «Великой депрессии». Значения

и оставим как в «лучшей» Moex-модели, ожидаемое среднее установим равным 300, а ожидаемое стандартное отклонение - 50, откуда получим . Параметр здесь подобран равным 50. Результат на рисунке 32.

Рис. 32. Результат применения Bayesian Online CPD (, , ) для кризиса DJIA 1929 года.

Источник: расчеты автора.

Алгоритм определил 10 критических точек, не считая начальной точки последовательности: 4 и 5 июля 1986, 20 августа 1986, 10 и 11 октября 1986, 22 мая 1929, 24 октября и с 28 по 30 октября 1929. Как видим, проблема подряд идущих точек осталась, хотя здесь их меньше, чем в случае Moex-выборки. Теперь с большей уверенностью можно сказать, что это проблема стартовых коэффициентов модели. В целом же алгоритм показывает себя неплохо - большинство определенных критических точек соответствуют видимым изменениям последовательности, а «кризисная» точка вновь определена очень близко к той, что ранее определил CUSUM (и вновь на день позже).

Напоследок рассмотрим ситуацию «Черного понедельника». Значение оставим не тронутым, коэффициент в этот раз возьмем равным 0.1, ожидаемое среднее - 2000, а ожидаемое стандартное отклонение - 400, откуда . Параметр . Результат можно увидеть на рисунке 33.

Рис. 33. Результат применения Bayesian Online CPD (, , ) для кризиса DJIA 1987 года.

Источник: расчеты автора.

Здесь алгоритм определил 4 критические точки, не считая начальной точки последовательности: 19 июля 1987, 16, 19 и 20 октября 1987. При этом на графике распределений можно увидеть, что есть еще как минимум две точки, скорее всего соответствующие локальным изменениям, однако, определенная для них вероятность немного не дотягивает до ...