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

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

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

Функция add_regressor()

Эта функция имеет следующие аргументы:
  • m - модельный объект;
  • name - название добавляемого предиктора;
  • prior.scale - параметр, используемый для регуляризации эффекта добавляемого предиктора (в этой статье не рассматривается, однако см., например, здесь);
  • standardize - позволяет стандартизовать значения добавляемого предиктора перед подгонкой модели (по умолчанию принимает значение "auto" - в этом случае предиктор будет стандартизован, если только он не является индикаторной переменной со значениями 1 и 0; другие возможные значения этого аргумента - TRUE и FALSE);
  • mode - необязательный параметр, определяющий характер сезонности добавляемого предиктора (по умолчанию равен m$seasonality.mode).
Все предикторы, добавляемые в модель с помощью функции add_regressor(), должны присутствовать в таблице с обучающими данными. Аналогично, для расчета предсказаний будущие значения каждого предиктора должны присутствовать в таблице с датами, задающими прогнозный отрезок времени. Последнее обстоятельство делает расчет прогноза с использованием количественных предикторов и многих качественных переменных проблематичным: будущие значения таких предикторов (в отличие, например, от дат официальных праздников и других регулярных событий) обычно исследователю не известны. Как правило, для решения этой проблемы сначала строят отдельные модели временных рядов для каждого предиктора, а затем используют предсказанные с помощью таких моделей будущие значения предикторов для получения искомого прогноза. Конечно, такой подход существенно усложняет весь процесс моделирования (поскольку возникает необходимость построения надежных моделей для отдельных предикторов), однако в большинстве случаев это - лучшее, что можно сделать.  Именно такой стратегией мы воспользуемся и в описанном ниже примере.

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

В качестве примера добавим в модель стоимости биткоина цены на акции трех компаний - Amazon, Google и Facebook. Стоит подчеркнуть, что выбор цен на акции в качестве предикторов, равно как и выбор именно этих компаний, ни с чем не связаны - просто такие данные легко получить из публично доступных источников. Для сбора этих данных воспользуемся пакетом tidyquant:

library(tidyquant)
library(prophet)
library(dplyr)
library(ggplot2)

# Собираем данные по цене акций Amazon, Google и Facebook на момент
# закрытия торгов, переименовываем переменные для удобства и
# логарифмируем значения цены:
AMZN <- tq_get(x = "AMZN", get = "stock.prices", 
               from = "2016-01-01", to = "2019-05-26") %>% 
    rename(ds = date, amzn = close) %>% 
    select(ds, amzn) %>% 
    mutate(amzn = log(amzn))

GOOG <- tq_get(x = "GOOG", get = "stock.prices", 
               from = "2016-01-01", to = "2019-05-26") %>% 
    rename(ds = date, goog = close) %>% 
    select(ds, goog) %>% 
    mutate(goog = log(goog))

FB <- tq_get(x = "FB", get = "stock.prices", 
             from = "2016-01-01", to = "2019-05-26") %>% 
    rename(ds = date, fb = close) %>% 
    select(ds, fb) %>% 
    mutate(fb = log(fb))

# Добавляем данные по цене акций в таблицу с обучающими данными: 
train_df <- left_join(train_df, AMZN) %>% 
    left_join(., GOOG) %>% 
    left_join(FB)

Обратите внимание: собранные с помощью функции tq_get() данные относятся к периоду с 2016-01-01 по 2019-05-26 (включительно), что соответствует периоду обучающих данных (таблица train_df). Важно также понимать, что на момент написания этой статьи (5 октября 2019 г.) интересующий нас прогнозный период (2019-05-27 - 2019-08-24) уже стал историей и, соответственно, цены на акции для этого периода были известны. Но, конечно же, мы не можем использовать эти ставшие историей данные для прогнозирования стоимости биткоина на период с 2019-05-27 по 2019-08-24: в реальной ситуации при расчете прогноза подобное "заглядывание в будущее" было бы невозможно.

Есть одна небольшая проблема с собранными нами данными по цене акций Amazon, Google и Facebook: все три переменные имеют пропущенные значения, что связано с прекращением торгов в выходные дни и во время государственных праздников США.

train_df

# A tibble: 1,242 x 5
       y ds          amzn  goog    fb
   <dbl> <date>     <dbl> <dbl> <dbl>
 1  6.07 2016-01-01  NA    NA    NA   
 2  6.07 2016-01-02  NA    NA    NA   
 3  6.06 2016-01-03  NA    NA    NA   
 4  6.07 2016-01-04  6.46  6.61  4.63
 5  6.07 2016-01-05  6.45  6.61  4.63
 6  6.06 2016-01-06  6.45  6.61  4.63
 7  6.13 2016-01-07  6.41  6.59  4.58
 8  6.12 2016-01-08  6.41  6.57  4.58
 9  6.10 2016-01-09  NA    NA    NA   
10  6.10 2016-01-10  NA    NA    NA   
# ... with 1,232 more rows

Хотя Prophet-модели мало чувствительны к наличию пропущенных значений в зависимой переменной (Taylor & Letham 2017), пропущенные значения в предикторах недопустимы. Для восстановления NA-значений в переменных amzn, goog и fb мы воспользуемся простой и хорошо подходящей для этого случая стратегией LOСF ("last observation carried forward"), которая заключается в том, что каждое пропущенное значение заменяется на последнее предшествующее ему непропущенное значение (или на ближайшее следующее за ним непропущенное значение, если ряд начинается с NA, как в нашем случае). Одну из реализаций метода LOCF можно найти в функции na_locf() из пакета imputeTS:

library(imputeTS)

# Восстанавливаем пропущенные значения цены акций по методу LOCF:
train_df <- train_df %>% 
    mutate(amzn = na_locf(amzn),
           goog = na_locf(goog),
           fb = na_locf(fb))

train_df

# A tibble: 1,242 x 5
       y ds          amzn  goog    fb
   <dbl> <date>     <dbl> <dbl> <dbl>
 1  6.07 2016-01-01  6.46  6.61  4.63
 2  6.07 2016-01-02  6.46  6.61  4.63
 3  6.06 2016-01-03  6.46  6.61  4.63
 4  6.07 2016-01-04  6.46  6.61  4.63
 5  6.07 2016-01-05  6.45  6.61  4.63
 6  6.06 2016-01-06  6.45  6.61  4.63
 7  6.13 2016-01-07  6.41  6.59  4.58
 8  6.12 2016-01-08  6.41  6.57  4.58
 9  6.10 2016-01-09  6.41  6.57  4.58
10  6.10 2016-01-10  6.41  6.57  4.58
# ... with 1,232 more rows

Как отмечено выше, цены на акции Amazon, Google и Facebook выбраны в качестве предикторов исключительно из удобства. Тем не менее, все три переменные демонстрируют определенную корреляцию со стоимостью биткоина и, соответственно, вполне подходят для моделирования:

Рис. 1

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

# Строим отдельные модели для каждого предиктора (параметры моделей
# подобраны после предварительного экспериментирования):
m_amzn <- prophet(rename(AMZN, y = amzn),
                  n.changepoints = 15, 
                  changepoint.range = 0.9,
                  weekly.seasonality = FALSE)
m_goog <- prophet(rename(GOOG, y = goog),
                  n.changepoints = 15, 
                  changepoint.range = 0.9,
                  weekly.seasonality = FALSE)
m_fb <- prophet(rename(FB, y = fb),
                n.changepoints = 15, 
                changepoint.range = 0.9,
                weekly.seasonality = FALSE)

# Рассчитываем прогнозные значения для цен на акции:
future_df_shares <- future_df %>% select(ds)
amzn_forecast <- predict(m_amzn, future_df_shares)
goog_forecast <- predict(m_goog, future_df_shares)
fb_forecast <- predict(m_fb, future_df_shares)

future_df_shares <- future_df_shares %>% 
    mutate(amzn = amzn_forecast$yhat,
           goog = goog_forecast$yhat,
           fb = fb_forecast$yhat)

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

Рис. 2

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


# Модель стоимости биткоина, включающая три предиктора:
M16 <- prophet(n.changepoints = 15, changepoint.range = 0.9)
M16 <- add_regressor(M16, 'amzn')
M16 <- add_regressor(M16, 'goog')
M16 <- add_regressor(M16, 'fb')
M16 <- fit.prophet(M16, train_df)

forecast_M16 <- predict(M16, future_df_shares)

# Визуализация полученной модели и прогнозных значений:
plot(M16, forecast_M16)

Рис. 3

На рис. 4 представлены оцененные компоненты модели. Обратите внимание, что эффекты всех трех предикторов объединены в один компонент (см. "extra_regressors_additive" на графике внизу):

# Визуализация компонентов модели:
prophet_plot_components(M16, forecast_M16)

Рис. 4

С помощью рассмотренной ранее функции plot_forecast_component() мы можем также изобразить эффекты отдельных предикторов (рис. 5):

library(gridExtra)

amzn_comp <- plot_forecast_component(M16, forecast_M16, name = "amzn") +
    ggtitle("Amazon")
goog_comp <- plot_forecast_component(M16, forecast_M16, name = "goog") +
    ggtitle("Google")
fb_comp <- plot_forecast_component(M16, forecast_M16, name = "fb") +
    ggtitle("Facebook")
grid.arrange(amzn_comp, goog_comp, fb_comp, ncol = 3)

Рис. 5

Согласно полученной модели, эффект переменной amzn оказался сильнее (примерный диапазон значений от -0.2 до 0.2), чем эффекты goog (от -0.005 до 0.005) и fb (от -0.06 до 0.06). Мы можем сделать такое заключение благодаря тому, что эффекты всех трех предикторов представлены на одной шкале и поэтому сравнимы (значения всех этих переменных перед подгонкой модели были стандартизованы  - см. выше описание аргумента standardize функции add_regressor()). Интересно также, что характер изменения эффекта amzn во времени почти зеркально отличается от такового у переменных goog и fb. Я оставлю читателю возможность подумать над интерпретацией этого наблюдения.

Как насчет качественных предикторов?

В рассмотренном нами примере все добавленные в модель дополнительные переменные были количественными. В случае с качественными предикторами у исследователя есть два варианта на выбор: либо представить значения таких переменных в виде "праздников", либо воспользоваться описанной выше более общей функцией add_regressor(). В обоих случаях качественные переменные необходимо преобразовать в индикаторные (т.е. создать отдельные столбцы для каждого уровня переменной со значениями 1/0).


***
В следующем сообщении из этой серии мы рассмотрим способ выбора оптимальной модели временного ряда с использованием перекрестной проверки.

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

4 Комментарии

В.Шитиков написал(а)…
Прекрасно рассказано, дорогой Сергей Эдуардович!
Но вот такой вопрос: Прогнозируется не вещественная переменная, а вероятность какого-то события. Например, ежедневно фиксировались частота покупок определенного товара из покупателей некой фокус-группы и нужно оценить динамику вероятности продаж на какой-то период. Есть ли какие-то принципиальные трудности замены отклика - натуральной переменной на частоты произошедших событий при использовании функций пакета? Всего наилучшего.
Sergey Mastitsky написал(а)…
Здравствуйте, Владимир Кириллович!

Как таковой возможности моделировать пропорции в Prophet, к сожалению, нету. Но возможны следующие подходы.

1. Если очень хочется использовать Prophet, то:
1a: во избежание возможных ситуаций, когда предсказания будут бессмысленными (например, отрицательные значения), можно попробовать моделировать не исходные пропорции, а их логарифм (с последующим экспоненцированием предсказанных значений для перехода на оригинальную шкалу)
1b: построить две модели - одну для количества единиц интересующего товара, а вторую - для общего количества всех остальных товаров. Далее на основе полученных прогнозов по каждой модели можно было бы рассчитать долю для интересующего товара на каждую дату прогнозного периода. Я бы пошел именно этим путем, поскольку он дал бы бизнесу более полную картину: и возможные величины продаж в абсолютном выражении, и долю интересующего товара.

2. Воспользоваться другими подходами, в частности GLARMA (https://cran.r-project.org/web/packages/glarma/vignettes/glarma.pdf)
Анонимно написал(а)…
Добрый день! Не планируете рассмотреть метод анализа временных рядов Singular spectrum analysis (SSA), что реализован в пакете Rssa?
Sergey Mastitsky написал(а)…
Спасибо за вопрос! Рассмотрение SSA не планировал, но добавлю в свой список тем, которые потенциально интересны читателям блога.
Новые Старые