Выбор прогнозной модели

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

Вообще все методы выбора моделей можно сгруппировать следующим образом:

  1. Точность аппроксимации;
  2. Перекрёстная проверка (Cross-validation);
  3. Экспертная оценка;

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

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

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

Однако обо всём по порядку.

Точность аппроксимации

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

Самым простым методом выбора в этой категории является выбор на основе скорректированного коэффициента детерминации (\(R^2\) adjusted). Может возникнуть вопрос, зачем же клепать новые коэффициенты, почему нельзя использовать простой \(R^2\) из параграфа «Оценка качества прогнозных моделей» ? Дело в том, что \(R^2\) не учитывает число параметров модели и поэтому всегда будет выше для моделей с большим числом переменных. Действительно, дополнительные параметры обычно дают любой модели большую гибкость, которая не обязательно приводит к реальному улучшению прогнозной точности, но при этом позволяет лучше описать данные по обучающей выборке. Поэтому-то скорректированный коэффициент детерминации и был разработан. Вот его формула:

\begin{equation} \label{eq:ms_R2adj}
R^2_{adj} = 1-\frac{(T-1)}{(T-k)}R^2 = 1-\frac{(T-1) SSE}{(T-k) TSS}
\end{equation}

Из формулы \eqref{eq:ms_R2adj} видно, что с ростом числа параметров \(k\), \(R^2_{adj}\) будет снижаться, а с ростом числа наблюдений \(T\) он будет приближаться к значению коэффициента \(R^2\). Его интерпретация при этом соответствует интерпретации \(R^2\), так что подробней на нём мы останавливаться не будем. Заметим, однако, что обычно \(R^2\) используется в регрессионном анализе. Выбрать модель экспоненциального сглаживания или какую-нибудь авторегрессию с помощью него затруднительно. Впрочем, даже регрессионные модели обычно выбирают на основе чего-нибудь более сложного, чем \(R^2_{adj}\).

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

В основе информационных критериев лежит идея о некой идеальной модели, лежащей глубоко в изучаемом нами процессе (та самая абстракция, о которой мы говорили в параграфе про системы и модели). Идея метода, предложенного впервые в 1974 году Hirotugu Akaike , заключается в том, что каждая модель, параметры которой мы оцениваем, находится на некотором расстоянии от идеальной модели. Расстояние это можно измерить с помощью функции правдоподобия (см. параграф про методы оценки), и называется оно «Дивергенция Кулбака-Либлера» («Kullback–Leibler divergence»). Рассчитываться оно для некоторой m-й модели должно примерно по следующей формуле:
\begin{equation} \label{eq:ms_KLD}
I_m = \ell^*-\ell_m ,
\end{equation}
где \(\ell^*\) — значение логарифмированной функции правдоподобия для идеальной модели, а \(\ell_m\) — значение логарифмированной функции правдоподобия для имеющейся m-й модели.

Особенностью значения \(I_m\) является то, что это число почти всегда будет положительным, а равно нулю лишь в том случае, если мы имеем дело с идеальной моделью. Поэтому знать \(\ell^*\) в принципе и не обязательно: чем больше значение \(\ell_m\), тем ближе модель находится к идеальной. Стало быть для выбора наилучшей модели достаточно рассчитать \(\ell_m\) для всех моделей в пуле и выбрать ту, у которой это значение максимально. Однако из-за ряда особенностей выборочного метода использовать непосредственно логарифмированную функцию правдоподобия \(\ell_m\) нельзя, нужно её скорректировать на число оцениваемых параметров \(k\). Akaike предложил рассчитывать «An Information Criterion», скорректированный на число параметров, по следующей формуле:
\begin{equation} \label{eq:ms_AIC}
AIC = 2k-2 \ln(L(\theta|Y)),
\end{equation}
где \(\theta\) — вектор параметров модели (длиной \(k\)), \(Y\) — вектор фактических значений. Двойка в формуле \eqref{eq:ms_AIC} была добавлена Akaike для того, чтобы распределение дисперсий получаемых оценок соответствовало \(\chi^2\). В принципе, она в формуле ненужна, но так уж сложилось… Если убрать из формулы \eqref{eq:ms_AIC} эту несчастную двойку, то получится выборочная оценка функции правдоподобия «\(-\ell_m\)». Обратите внимание на минус! В формуле \eqref{eq:ms_AIC} используется негативное значение логарифмированной функции правдоподобия, поэтому выбирать нужно ту модель, у которой AIC будет наименьшим. Вот такая загагулина!

Иногда AIC так же интерпретируют как показатель, который основывается на принципе бритвы Оккама: модели с большим числом параметров будут иметь большее значение \(k\), что в свою очередь приведёт к большему значению AIC. То есть, якобы, таким образом сохраняется принцип: «не стоит плодить сущностей сверх необходимых» — и предпочтение отдаётся модели, описывающей процесс лучше, но при этом с меньшим числом параметров. Однако это скорее просто приятное дополнение. Сам коэффициент, как видим, выводится из совсем других предпосылок.

Но это ещё не всё. В 1978 году Nariaki Sugiura предложил модификацию AIC для случая с нормальным распределением, скорректировав коэффициент на число наблюдений в выборке. Итоговый критерий называется AICc (скорректированный AIC) и рассчитывается по следующей формуле:
\begin{equation} \label{eq:ms_AICc}
AICc = AIC + 2\frac{k(k+1)}{T-k-1} = 2 k \frac{T}{T-k-1}-2 \ln(L(\theta|Y)),
\end{equation}
где \(T\) — число наблюдений в выборке.

По первой части формулы \eqref{eq:ms_AICc} можно сказать, что с ростом числа наблюдений \(T\), разница между AIC и AICc будет уменьшаться, так как будет расти знаменатель второго члена. Вторая же часть формулы показывает, в чём именно заключается разница между AIC и AICc. В связи с тем, что на малых выборках AICc должен быть точнее, а на больших разницы практически не будет, статистики и прогнозисты рекомендуют для выбора модели использовать именно AICc. Однако надо заметить, что AICc применим только для моделей с нормальным распределением остатков (лог-нормальное тоже подойдёт). Для других распределений, нужно выводить другие формулы для AICc.

Среди информационных критериев встречаются и другие. Например, есть BIC («Bayesian IC» aka SIC — «Scwarz IC»), который выводится с использованием Байесовской статистики и рассчитывается по формуле:
\begin{equation} \label{eq:ms_BIC}
BIC = k \ln(T)-2 \ln(L(\theta|Y)),
\end{equation}

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

Пример в R
Рассмотрим, как выбор моделей работает на примере функции ets() пакета forecast для R. Для этого построим несколько моделей и сравним значения их информационных критериев.

x <- rnorm(100,0,1)

ets(x,model="ANN")
ets(x,model="AAN")
ets(x,model="AAA")

В результате выполнения этого кода, вы увидите несколько значений информационных критериев. Как вы считаете, какой из этих моделей на основе AICc нужно отдать предпочтение? Уверен на 146%, что это первая модель.

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

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

Второе "но" -- выбирать можно только среди моделей одного и того же типа, сформулированных одинаково. Например, нельзя выбрать с помощью IC между ETS в представлении пространства состояний и ARIMA в классическом виде. Чтобы это было возможно, их нужно привести к общему виду (например, пространство состояний).

Третье "но" заключается в том, что сравнивать с помощью IC можно только модели с одинаковой зависимой переменной. Сравнить, например, модель \(y_t=f(x_t)\) с \(\ln y_t = g(x_t)\) с помощью информационных критериев нельзя. Здесь уже нужно либо прибегать к ухищрениям по предварительному преобразованию переменной (например \(y_t = \exp(g(x_t)) \)), либо использовать другую функцию распределения при расчёте функции правдоподобия (например, лог-нормальную вместо нормальной), либо для выбора использовать какие-нибудь другие методы.

Подробней об информационных критериях можно почитать в замечательной книге "Burnham & Anderson, 1998. Model Selection and Multimodel Inference" . Саму тему использования информационных критериев мы ещё затронем, но уже применительно к конкретным моделям.

Перекрёстная проверка

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

Одним из достаточно популярных методов является нечто под названием "k-fold cross-validation". Идея его проста до безобразия. Имеющаяся выборка разбивается на некоторое число частей \(k\), после чего по \(k-1\) частям модель оценивается, а по оставшейся одной -- проверяется её точность с помощью каких-нибудь коэффициентов оценки точности прогнозов из параграфа про оценку качества моделей. После этого вырванную часть возвращают на место, а из ряда вытаскивают следующую, модель опять оценивают по \(k-1\) и проверяют по той самой отложенной в сторону. В результате этого в распоряжении исследователя оказывается \(k\) рядов с прогнозными ошибками, по которым можно принять решение о том, насколько модель хороша по сравнению с другими. Проблема метода заключается в том, что в случае с временными рядами разделение выборки на такие части приводит к потере структуры ряда. Например, имея месячные данные по продажам футболок с Путиным за 2011 -- 2016 годы, нельзя адекватно оценить модель, если вырван кусок ряда данных за 2014 год. Вторая проблема метода состоит в том, что весь процесс прогнозирования и выбора модели с ним усложняется и замедляется в \(k\) раз. Однако метод может быть эффективно применён при построении моделей регрессий и некоторых видов нейронных сетей.

Ну, и, конечно же, возникает закономерный вопрос: что такое \(k\), и каким оно должно быть? Всё очень просто. Достаточно часто на практике берут \(k=10\). Почему? А почему бы и нет?! Красивое же число! Но, если серьёзно, то число должно быть таким, чтобы по наблюдениям в \(k-1\) частей ряда можно было оценить модель без серьёзных проблем (см. параграф про базовые принципы).

Частным случаем "k-fold" является метод "leave-one-out" (мой вульгарный перевод на русский с помощью переводчика "Prompt" -- "отпуск один внешний"). Идея сводится к тому, чтобы оценить выборку по \(T-1\) наблюдению, а оценить по оставшемуся одному. Этот процесс повторяется \(T\) раз: каждый раз убирается одна точка, а по оставшимся модель оценивается. Минусы у метода примерно такие же, как и "k-fold", но вычислений оказывается значительно больше, так как фактически теперь \(k=T\). Явным плюсом выступает большое число наблюдений в проверочной выборке. В этом случае можно оценивать прогнозную точность модели, используя любые понравившиеся исследователю статистические методы.

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

Метод фиксированной точки отсчёта ("Fixed origin") фактически просто соответствует ретропрогнозу: мы убираем последние \(h\) наблюдений из выборки, строим модель по первым \(T-h\), после чего проверяем точность модели с помощью коэффициентов из параграфа про оценку качества моделей по тем самым убранным наблюдениям. \(h\) в данном случае -- это срок прогнозирования. Метод прост и эффективен, имеет понятную идея, но вообще-то не идеален. Давайте поясню.

Предположим, что мы пытаемся дать прогноз количества туристов, въезжающих в Россию на 2015-й год на основе квартальных данных за 2000 -- 2014 годы. 2014-й год мы оставили для проверки (как того и требует метод фиксированной точки), сделали на него прогноз и поняли, что всё плохо, всё тлен и ничего не работает кроме простой средней или какого-нибудь метода "Naïve". Как же нам интерпретировать этот результат и что же делать, если нам всё равно надо выбрать модель? Для начала надо понять, что же произошло. А дело в том, что 2014-й год аномальный -- в это время произошли качественные изменения в политике и экономике России, которые ни одна математическая модель не уловит. Поэтому проверять и сравнивать по нему модели не корректно -- эту аномалию нужно учесть либо через какие-нибудь фиктивные переменные (о которых мы ещё когда-нибудь поговорим), либо используя экспертные методы. А для выбора модели можно помимо 2014-го года использовать ещё и 2012-й, 2013-й.

Собственно говоря так мы плавно и подошли к методу сдвигающейся точки отсчёта ("Rolling origin"). Идея метода чем-то похожа на "k-fold", но при этом сохраняется структура временного ряда. Для начала исследователь выбирает некий горизонт \(r\) для сдвигающейся точки (сколько наблюдений с конца ряда убрать для процедуры), после чего выполняет следующую последовательность действий:

  1. оценивает модель по \(T-r\) наблюдениям,
  2. делает прогноз на \(h\) наблюдений вперёд из точки \(T-r\),
  3. записывает ошибки полученных прогнозов,
  4. добавляет одно наблюдение \(T-r+1\) в обучающую выборку, в результате чего его горизонт \(r = r-1\) (то есть уменьшается на одно наблюдение),
  5. повторяет процедуру, начиная с п.1 до тех пор, пока либо не кончатся наблюдения, либо не лопнет терпение.

В результате такой процедуры в распоряжении прогнозиста оказывается не \(h\) ошибок, а значительно больше. Причём, количество ошибок зависит от того, когда процесс останавливается. Если это происходит в тот момент, когда в проверочной выборке осталось \(h-1\) наблюдение (то есть полноценный прогноз на весь горизонт уже не сделать), то набор ошибок составит \(h \cdot (r-h+1) \). Если же продолжать до самого конца ряда данных, с каждым шагом уменьшая горизонт прогноза, то ошибок окажется ещё больше, однако не все они будут одинаково полезными: в этом случае получится \(r\) одношаговых прогнозных ошибок, \(r-1\) ошибок на два шага вперёд и так далее.

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

Пример в R
Посмотрим как работает функция сдвигающейся точки ro() из пакета greybox, которая была написана мною совместно с Ивом Сагаэром (Yves Sagaert), аспирантом университета Гент (Бельгия). Сгенерируем какой-нибудь ряд и возьмём функцию экспоненциального сглаживания es() из пакета smooth.

x <- rnorm(100,0,1)

ourcall <- "es(data,model='ANN',h=h)"
ourvalue <- "forecast"

ro(x,h=5,origins=5,ourcall,ourvalue,co=TRUE)

Переменная "ourcall" содержит вызов функции es() со всеми желаемыми параметрами и показывает, куда в этом вызове вставить данные "data", а куда -- горизонт прогнозирования \(h\). "ourvalue" говорит функции, какие значения вернуть по результатам использования нашей функции. Параметр "co=TRUE" в вызове функции ro() говорит о том, что обучающая выборка должна иметь фиксированную величину, то есть остановиться нужно тогда, когда наблюдений в ней осталось меньше \(h\). Функция достаточно универсально, но требует времени для освоения.

В результате выполнения этого кода мы получим две матрицы размером 5 на 5: "forecast" будет содержать прогнозы по нашей модели, а "actuals" будет содержать соответствующие фактические значения. По этим двум матрицам можно рассчитать ошибки на разные горизонты и прийти к каким-нибудь выводам.

Заметим, что метод сдвигающейся точки может использоваться по-разному. Например, иногда модель не переоценивают при добавлении наблюдений в обучающую выборку. Или добавляют не одно, а несколько наблюдений. Помимо этого, никто не знает, каким должен быть размер выборки для сдвигающейся точки \(r\). Тут так же надо руководствоваться здравым смыслом и базовыми принципами прогнозирования. В общем, всё это никак не регламентировано и статистически не обосновано. Помогает лишь теория прогнозирования.

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

Теперь, когда мы разобрались с инструментарием прогнозиста, мы, наконец, можем перейти к самому интересному -- методам и моделям.