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

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

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

В 2017 г. специалисты компании Facebook объявили о разработанном ими новом пакете для прогнозирования временных рядов - Prophet ("пророк"). Prophet решает указанные выше недостатки forecast и других подобных пакетов и позволяет в автоматическом режиме создавать модели, обладающие высокой точностью предсказаний. В следующих нескольких сообщениях мы рассмотрим основные возможности Prophet и особенности работы с ним.

Как Prophet работает

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

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

Оценивание параметров подгоняемой модели выполняется с использованием принципов байесовской статистики (либо методом нахождения апостериорного максимума (MAP), либо путем полного байесовского вывода). Для этого применяется платформа вероятностного программирования Stan. Prophet представляет собой ни что иное, как удобный интерфейс для работы с этой платформой из среды R (имеется также аналогичная библиотека для Python - fbprophet).

Инсталляция и документация

Prophet распространяется бесплатно по лицензии MIT. Этот пакет легко установить стандартным образом из CRAN (на Windows-машинах сначала нужно будет установить Rtools):

install.packages("prophet")

Если Вы работаете с компьютером Mac под управлением OS X, не забудьте добавить аргумент type = "source":

install.packages("prophet", type = "source")

Подробную документацию по работе с Prophet можно найти на Github-сайте проекта.

Данные, используемые в примерах

Для иллюстрации особенностей работы с Prophet воспользуемся историческими данными по стоимости биткоина на момент закрытия торгов (в долларах США). Подобные данные можно найти в удобной для скрейпинга табличной форме на нескольких сайтах, например CoinMarketCap (на момент написания этой статьи доступные данные охватывали период с 2013-04-28 по 2019-08-24 включительно):

require(rvest) # пакет для скрейпинга
require(dplyr)
require(lubridate)
require(ggplot2)
require(readr)
require(prophet)

# Собираем данные с сайта CoinMarketCap (за период с 2013-04-28 по 2019-08-24):
base_url <- "https://coinmarketcap.com/"
currency <- "currencies/bitcoin/historical-data/"
period <- "?start=20130428&end=20190825"
page <- paste0(base_url, currency, period)

dat <- read_html(page) %>% 
    html_table() %>% 
    .[[3]] %>% 
    rename(y = "Close**", ds = Date) %>% 
    select(y, ds) %>% 
    mutate(ds = mdy(ds)) %>% 
    mutate(y = gsub(",", "", y) %>% as.numeric(.))

# Представим собранные данные графически:
dat %>% 
    ggplot(., aes(ds, y)) +
    geom_line() + 
    theme_minimal()

Рисунок 1

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

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

# Удаляем данные, которые старше 2016-01-01 (чтобы сделать ряд немного
# короче и ускорить вычисления), и логарифмируем отклик (логарифмирование 
# поможет нам "бороться" с большой дисперсией в данных):
clean_data <- dat %>% 
    filter(ds >= as.Date("2016-01-01")) %>% 
    mutate(y = log(y))

# Разбиваем данные на обучающую выборку (все наблюдения, за исключением последних
# 90 дней) и проверочную выборку (последние 90 дней):
train_df <- clean_data %>% 
    arrange(ds) %>% 
    slice(1:(n() - 90))

test_df <- clean_data %>% 
    arrange(ds) %>% 
    tail(90)

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

Обучающие данные, подготовленные описанным выше способом, выглядят так:

Рисунок 2

Базовая модель

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

M0 <- prophet(train_df)

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

future_df <- make_future_dataframe(M0, periods = 90)
forecast_M0 <- predict(M0, future_df)

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

head(forecast_M0[, c("yhat", "yhat_lower", "yhat_upper")])
#     yhat yhat_lower yhat_upper
# 6.134964   6.033923   6.239148
# 6.134025   6.030052   6.231834
# 6.127710   6.017977   6.222559
# 6.124329   6.010517   6.229038
# 6.119314   6.016283   6.215383
# 6.111687   6.009227   6.211752

Таблицу forecast_M0 и объект M0 далее можно подать на функцию plot(), чтобы изобразить подогнанную модель и прогнозные значения на графике:

plot(M0, forecast_M0)

Рисунок 3

Черные точки на рис. 3 соответствуют значениям отклика из обучающей выборки. Сплошная голубая линия - это предсказанные моделью значения отклика, а огибающая эту линию светло-голубая "лента" - это 80% доверительные интервалы соответствующих предсказанных значений. Прогнозные значения у на следующие 90 дней видны в правой части графика.

Таблицу forecast_M0 и объект M0 можно также подать на функцию prophet_plot_components(),  что позволит изобразить отдельные компоненты модели :

prophet_plot_components(M0, forecast_M0)

Рисунок 4. Вверху - тренд; посредине - недельная сезонность; внизу - годовая сезонность.

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

Что дальше?

Для создания модели M0 потребовалась всего одна строка кода и процесс подгонки занял пару секунд. Учитывая, что анализируемый нами временной ряд - не самый простой для моделирования, полученная модель в целом хорошо передает свойства этого ряда. Возможность быстро и удобно создавать такие "хорошие" модели является основным преимуществом Prophet.

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

Подпишитесь на рассылку, чтобы не пропустить ничего интересного!