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

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

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

Для решения описанной проблемы при работе с временными рядами применяют несколько модификаций перекрестной проверки (подробнее см., например, Hyndman & Athanasopoulos 2019, Tashman 2000, а также документацию по пакету caret). В пакете Prophet, в частности, реализован т.н. метод "имитированных исторических прогнозов" (Simulated Historical Forecasts, SHF). Рассмотрим, что этот метод собой представляет, и как им пользоваться для выбора оптимальной модели.

Метод имитированных исторических прогнозов

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

Рис. 1

Метод SHF (рис. 2) пытается сымитировать описанную выше процедуру. В пределах отрезка с исходными обучающими данными выбирают \(K\) "точек отсчета" ("cut-off points" по терминологии Prophet), на основе которых формируются блоки данных для выполнения перекрестной проверки: все исторические наблюдения, предшествующие \(k\)-й точке отсчета (а также сама эта точка), образуют обучающие данные для подгонки соответствующей модели, а \(H\) исторических наблюдений, следующих за точкой отсчета, образуют прогнозный горизонт. Расстояние между точками отсчета называется "периодом" (period) и по умолчанию составляет \(H / 2\). Обучающие наблюдения в первом из \(K\) блоков образуют т.н. "initial period" ("начальный отрезок"). В Prophet длина этого отрезка по умолчанию составляет \(3 \times H\), однако исследователь при желании может ее изменить.

Каждый раз после подгонки модели на обучающих данных из \(k\)-го блока рассчитываются предсказания для прогнозного горизонта того же блока, что позволяет оценить качество прогноза с помощью подходящей случаю метрики (например, RMSE; см. также ниже). Значения этой метрики, усредненные по каждой дате прогнозных горизонтов каждого блока, в итоге дают оценку качества предсказаний, которую можно ожидать от модели, построенной по всем исходным обучающим данным. Это в свою очередь позволяет сравнить несколько альтернативных моделей и выбрать оптимальную.

Рис. 2

Таким образом, SHF близок к классическому способу выполнения перекрестной проверки моделей временных рядов по методу "скользящей точки отсчета" ("rolling origin evaluation", Tashman 2000), однако отличается от последнего тем, что для оценивания качества предсказаний используется меньше блоков и прогнозные горизонты этих блоков удалены друг от друга на некоторое расстояние. С одной стороны, это является преимуществом метода SHF, поскольку он требует меньше вычислений и получаемые оценки качества предсказаний не так сильно коррелируют друг с другом, как в случае с перекрестной проверкой по методу "скользящей точки отсчета". С другой стороны, при небольшом количестве блоков \(K\) оценка качества модели может оказаться ненадежной. С решением последней проблемы отчасти помогает увеличение отрезка исходных исторических данных (поскольку чем он длиннее, тем больше блоков можно в него "вместить"). Однако следует помнить, что модели, построенные на длинных временных рядах, часто демонстрируют низкое качество прогнозов, поскольку параметры таких моделей задать труднее и возрастает риск переобучения.

Код для воспроизведения примеров

Все примеры, приведенные в сообщениях из этой серии, можно воспроизвести с помощью кода, который хранится в Github-репозитории ranalytics/intro_to_prophet. Описание используемых в примерах данных по стоимости биткоина представлено в первом сообщении.

Функция cross_validation()

В пакете Prophet перекрестная проверка по методу имитированных исторических прогнозов выполняется с помощью функции cross_validation(), которая имеет следующие аргументы:

• model - модельный объект;
• horizon - длина прогнозного горизонта в каждом блоке данных, используемом для выполнения перекрестной проверки;
• units - название единицы измерения времени (например, "days", "hours", "secs" - см. справочный файл по базовой функции difftime());
• initial - длина начального отрезка с обучающими данным в первом блоке.

Функция cross_validation() возвращает таблицу c истинными (y) и оцененными (yhat) значениями моделируемой переменной, а также доверительными границами предсказанных значений (yhat_lower и yhat_upper) для каждой точки отсчета (cutoff) и каждой даты соответствующего прогнозного периода (ds):

library(prophet)
library(dplyr)

# Загружаем ранее созданные и сохраненные модели:
load("./workspaces/intro_to_prophet.RData")

# Пример использования функции cross_validation():
M3_cv <- cross_validation(M3, 
                          initial = 730, # длина первого блока (2 года)
                          period = 90,   # расстояние между точками отсчета
                          horizon = 90,  # длина прогнозного горизонта (90 дней)
                          units = "days")

head(M3_cv)

## # A tibble: 6 x 6
##   ds                      y  yhat yhat_lower yhat_upper cutoff             
##   <dttm>              <dbl> <dbl>      <dbl>      <dbl> <dttm>             
## 1 2018-03-03 00:00:00  9.35  9.23       9.12       9.35 2018-03-02 00:00:00
## 2 2018-03-04 00:00:00  9.35  9.21       9.09       9.33 2018-03-02 00:00:00
## 3 2018-03-05 00:00:00  9.36  9.21       9.09       9.31 2018-03-02 00:00:00
## 4 2018-03-06 00:00:00  9.29  9.19       9.07       9.30 2018-03-02 00:00:00
## 5 2018-03-07 00:00:00  9.21  9.18       9.07       9.30 2018-03-02 00:00:00
## 6 2018-03-08 00:00:00  9.15  9.17       9.05       9.28 2018-03-02 00:00:00

Функция performance_metrics()

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

• Среднеквадратичная ошибка (mean squared error, MSE): \[MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]
• Квадратный корень из среднеквадратичной ошибки (root mean squared error, RMSE): \[RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}\]
• Средняя абсолютная ошибка (mean absolute error, MAE): \[MAE = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i|\]
• Средняя абсолютная удельная ошибка (mean absolute percentage error, MAPE): \[MAPE =\frac{1}{n} \sum_{i=1}^{n} \bigl| \frac{y_i - \hat{y}_i}{y_i} \bigr| \]
• "Покрытие" (coverage): доля истинных значений моделируемой переменной, которые находятся в пределах доверительных границ прогноза.

В приведенных выше формулах \(y_i\) и \(\hat{y}_i\) - истинное и предсказанное значения моделируемой переменной соответственно, а \(n\) - количество наблюдений.

Функция performance_metrics() имеет следующие аргументы:

• df - таблица, полученная с помощью функции cross_validation();
• metrics - вектор с названиями метрик качества модели (по умолчанию этот аргумент принимает значение NULL, что приводит к расчету всех перечисленных выше метрик, т.е. c("mse", "rmse", "mae", "mape", "coverage"));
• rolling_window - размер "скользящего окна", в пределах которого происходит усреднение каждой метрики (по умолчанию принимает значение 0.1, т.е. 10% от длины прогнозного горизонта; см. пояснения ниже).

На рис. 3 представлено схематичное изображение того, как функция performance_metrics() рассчитывает метрики качества модели для случая, когда перекрестная проверка выполняется на 5 блоках данных с длиной прогнозного горизонта \(H = 100\) и аргументом rolling_window = 0.1.

Рис. 3

В приведенном ниже примере функция performance_metrics() применена для расчета среднеквадратичной ошибки прогноза модели построенной нами ранее модели M3:

performance_metrics(M3_cv,
                    metrics = "mse",
                    rolling_window = 0.1) %>% 
    head()

##   horizon        mse
## 1  9 days 0.06424673
## 2 10 days 0.07572021
## 3 11 days 0.08695892
## 4 12 days 0.09707747
## 5 13 days 0.10990773
## 6 14 days 0.12260315

Как видно из полученного результата, первое усредненное значение MSE приходится на 9-й день прогнозного горизонта, поскольку длина этого горизонта для модели M3 составляет 90 дней (см. выше код для создания таблицы M3_cv), а 10% от 90 (размер скользящего окна, задаваемый аргументом rolling_window) составляет 9.

Если аргументу rolling_window присвоить значение 0, то функция performance_metrics рассчитает запрашиваемые метрики качества для каждой даты прогнозного горизонта (т.е. размер скользящего окна в данном случае фактически равен 1):

performance_metrics(M3_cv,
                    metrics = "mse",
                    rolling_window = 0) %>% 
    head()

##   horizon        mse
## 1  1 days 0.03595728
## 2  2 days 0.03941046
## 3  3 days 0.05343681
## 4  4 days 0.05035319
## 5  5 days 0.05454100
## 6  6 days 0.07256616

Если же аргументу rolling_window присвоить значение 1, то запрашиваемые метрики качества будут усреднены по всему прогнозному горизонту:

performance_metrics(M3_cv,
                    metrics = "mse",
                    rolling_window = 1) %>% 
    head()

##    horizon        mse
## 1  90 days   0.192834

Функция plot_cross_validation_metric()

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

• df_cv - таблица, полученная с помощью функции cross_validation();
• metric - название метрики;
• rolling_window - размер "скользящего окна", в пределах которого происходит усреднение метрики (см. рис. 3 и описание функции performance_metrics()).

Подобно другим функциям для визуализации данных в пакете Prophet, plot_cross_validation_metric() возвращает объект класса ggplot (рис. 4):

# Визуализация метрики качества предсказаний, полученной по результатам 
# перекрестной проверки:
plot_cross_validation_metric(M3_cv,
                             metric = "mse",
                             rolling_window = 0.1)

Рис. 4

На рис. 4 приведены оценки MSE для каждой из дат прогнозного горизонта (\(H = 90\)) каждого из \(K = 5\) блоков данных, участвовавших в перекрестной проверке. Голубая линия соответствует усредненным значениям в пределах каждого скользящего окна размером в 9 наблюдений (см. рис. 3). Судя по большому разбросу полученных оценок MSE, качество модели M3 желает оставлять лучшего.

Пример выбора оптимальной модели с помощью перекрестной проверки

Рассмотрим теперь, как описанную выше методологию выполнения перекрестной проверки можно применить для выбора оптимальной модели из нескольких альтернативных. Предположим, что перед нами стоит задача выбрать оптимальную модель стоимости биткоина из построенных ранее моделей M4, M5 и M12. Для оценивания качества этих моделей воспользуемся двумя метриками - MAPE и покрытие (см. выше).  Для упрощения примера предположим также, что нас интересует качество предсказаний в целом для 90-дневного прогнозного горизонта (т.е. нам неинтересны отдельные даты этого горизонта). Рассчитаем обе метрики качества для каждой из моделей-кандидатов:

# Пример выбора оптимальной модели:
M4_cv <- cross_validation(M4, 
                          initial = 730,
                          period = 180,
                          horizon = 90,
                          units = "days")
M5_cv <- cross_validation(M5, 
                          initial = 730,
                          period = 180,
                          horizon = 90,
                          units = "days")
M12_cv <- cross_validation(M12, 
                           initial = 730,
                           period = 180,
                           horizon = 90,
                           units = "days")
M4_perf <- performance_metrics(M4_cv,
                               metrics = c("mape", "coverage"),
                               rolling_window = 1)
M5_perf <- performance_metrics(M5_cv,
                               metrics = c("mape", "coverage"),
                               rolling_window = 1)
M12_perf <- performance_metrics(M12_cv,
                                metrics = c("mape", "coverage"),
                                rolling_window = 1)

M4_perf
##   horizon       mape  coverage
## 1 90 days 0.03214045 0.2703704

M5_perf
##   horizon       mape  coverage
## 1 90 days 0.03714357 0.4185185

M12_perf
##   horizon       mape  coverage
## 1 90 days 0.02347762 0.6444444

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

Рис. 5

Вспомним, что до сих пор мы строили все модели по обучающим данным из таблицы train_df. Однако у нас есть и тестовый набор данных - test_df. Посмотрим, как выбранная нами оптимальная модель M12 сработает на этой тестовой выборке. На рис. 6 представлены обучающие данные (черные точки; для удобства показаны наблюдения только за 2019 г.) и истинные значения стоимости биткоина в прогнозном периоде (красные точки). Голубая сплошная линия на этом графике соответствует предсказанным моделью значениям, а светло-голубая полоса вокруг нее - 80%-ной доверительной области предсказанных значений:

plot(M12, forecast_M12) + 
    coord_cartesian(xlim = c(as.POSIXct("2019-01-01"),
                             as.POSIXct("2019-08-24"))) +
    geom_point(data = test_df, aes(as.POSIXct(ds), y), col = "red")

Рис. 6

Хотя выбранная нами в качестве оптимальной модель M12 не смогла правильно предсказать некоторые локальные колебания стоимости биткоина в прогнозном периоде, в целом она дала неплохой результат - большинство истинных значений стоимости оказалось в пределах 80%-ной доверительной полосы. Следует подчеркнуть, однако, что это всего лишь пример использования перекрестной проверки для сравнения качества предсказаний нескольких моделей - совершенно точно не стоит разрабатывать стратегию торговли биткоином на основе приведенных здесь результатов ;)

***

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

Другие статьи из серии:

• Введение
• Параметры моделей
• Эффекты "праздников"
• Сезонные компоненты
• Добавление предикторов