Прогнозирование временных рядов с использованием модели TBATS

Мультисезонность Модель прогнозирования временных рядов с кодами Python и R

Введение

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

TBATS — это метод прогнозирования для моделирования данных временных рядов. Основная цель этого - прогнозировать временные ряды со сложными сезонными закономерностями с использованием экспоненциального сглаживания.

TBATS: Tригонометрическая сезонность, преобразование Bокса-Кокса, ошибки ARMA, Trend и T strong>Sсезонные компоненты.

Обзор

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

Разделите данные на две части (скажем, train_data и test_data). Train_data используется для обучения модели и подгонки модели к данным. Обученная модель оценивается с использованием test_data.
Предоставьте модели информацию о продолжительности сезона (например, если присутствуют почасовые данные, модель можно строить еженедельно для всех 24*7 часов в неделю).
Подгоните модель к train_data, используя train_data к модели.
Прогнозируйте модель вперед на определенный период времени, для которого вы хотите сделать прогноз.

КАК TBATS ВЫБИРАЕТ ОКОНЧАТЕЛЬНУЮ МОДЕЛЬ

TBATS рассмотрит различные альтернативы и подойдет для нескольких моделей. Он будет рассматривать модели:

с преобразованием Бокса-Кокса и без преобразования Бокса-Кокса.
с учетом тренда и без тренда.
с демпфированием тренда и без демпфирования тренда.
с процессом ARIMA(p,q) и без ARMA(p,q), используемым для моделирования остатков.
несезонная модель.
различное количество гармоник, используемых для моделирования сезонных эффектов

Окончательная модель будет выбрана с использованием информационного критерия Акаике (AIC).

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

Существующая реализация

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

Новая реализация

Новая реализация TBATS в Python доступна на GitHub. В оставшейся части статьи мы приведем пример использования и сравним производительность этой реализации с другими методами.

Монтаж

Из ПиПИ:

pip install tbats

Импорт через:

from tbats import BATS, TBATS

Минимальный рабочий пример:

from tbats import TBATS
import numpy as np
# required on windows for multi-processing,
# see https://docs.python.org/2/library/multiprocessing.html#windows
if __name__ == '__main__':
    np.random.seed(2342)
    t = np.array(range(0, 160))
    y = 5 * np.sin(t * 2 * np.pi / 7) + 2 * np.cos(t * 2 * np.pi / 30.5) + \
        ((t / 20) ** 1.5 + np.random.normal(size=160) * t / 50) + 10
    
    # Create estimator
    estimator = TBATS(seasonal_periods=[14, 30.5])
    
    # Fit model
    fitted_model = estimator.fit(y)
    
    # Forecast 14 steps ahead
    y_forecasted = fitted_model.forecast(steps=14)
    
    # Summarize fitted model
    print(fitted_model.summary())

Чтение сведений о модели

# Time series analysis
print(fitted_model.y_hat) # in sample prediction
print(fitted_model.resid) # in sample residuals
print(fitted_model.aic)
# Reading model parameters
print(fitted_model.params.alpha)
print(fitted_model.params.beta)
print(fitted_model.params.x0)
print(fitted_model.params.components.use_box_cox)
print(fitted_model.params.components.seasonal_harmonics)

Поиск неисправностей

BATS и TBATS пробуют множество моделей под капотом и могут показаться медленными при подгонке слишком длинных временных рядов. Чтобы ускорить его, вы можете начать с ограниченного пространства поиска модели. Рекомендуется запускать его без преобразования Бокса-Кокса и моделирования ошибок ARMA, которые являются самыми медленными элементами модели:

# Create estimator
estimator = TBATS(
    seasonal_periods=[14, 30.5],
    use_arma_errors=False,  # shall try only models without ARMA
    use_box_cox=False  # will not use Box-Cox
)
fitted_model = estimator.fit(y)

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

estimator = TBATS(
    seasonal_periods=[14, 30.5],
    n_jobs=1
)
fitted_model = estimator.fit(y)

Сравнительные тесты пакета прогнозов R. Те, которые НЕ ЗАПУСКАЮТСЯ с тестовой командой по умолчанию, вам необходимо установить пакет прогноза R:

python setup.py test_r

Сравнение с реализацией R

Реализация Python должна быть максимально эквивалентна реализации R в пакете прогнозов.

Некоторые примеры и набор данных временных рядов

Чтобы протестировать методы прогнозирования, нам нужны некоторые данные временных рядов. Давайте воспользуемся временными рядами из Вызова прогнозирования спроса на товары в магазине Kaggle. Это игровой вызов, и набор, скорее всего, искусственный (см. комментарии в ядрах и обсуждениях). Данные там содержат ежедневные продажи 50 товаров в 10 магазинах за период 5 лет (всего 500 различных временных рядов). Для нашей цели нам нужен только один временной ряд, поэтому я произвольно возьму продажи товара 1 в магазине 1.

import pandas as pddf = pd.read_csv('kaggle_sales.csv')
df = df[(df['store'] == 1) & (df['item'] == 1)] # item 1 in store 1
df = df.set_index('date')
y = df['sales']y_to_train = y.iloc[:(len(y)-365)]
y_to_test = y.iloc[(len(y)-365):] # last year for testing

Рис. 1. Ежедневные продажи товара 1 в магазине 1.

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

Модель TBATS

Чтобы начать прогнозирование, нам нужно установить пакет tbats и подогнать модель. Единственное, что мы должны вручную предоставить модели, — это продолжительность сезона:

from tbats import TBATS, BATS# Fit the model
estimator = TBATS(seasonal_periods=(7, 365.25))
model = estimator.fit(y_to_train)# Forecast 365 days ahead
y_forecast = model.forecast(steps=365)

Вы могли заметить, что продолжительность годового сезона не является целым числом. Он равен 365,25 для учета високосных лет, с чем может справиться TBATS.

TBATS, кажется, проделал неплохую работу по моделированию обоих сезонных эффектов:

Рис. 2: данные о продажах за 3 последних года. TBATS моделирует ежегодные сезонные эффекты.

Рис. 3: данные за 12 недель. Еженедельный сезонный эффект также моделируется TBATS.

Если мы заглянем под капот и просмотрим параметры модели, мы обнаружим, что 3 сезонные гармоники используются для моделирования недельной модели, а 11 гармоник используются для моделирования годовой модели. TBATS решил использовать преобразование Бокса-Кокса с лямбдой 0,234955. Тренд не моделируется, и ARMA не используется для моделирования остатков, поскольку p, q равно 0.

Use Box-Cox: True
Use trend: False
Use damped trend: False
Seasonal periods: [  7.   365.25]
Seasonal harmonics [ 3 11]
ARMA errors (p, q): (0, 0)
Box-Cox Lambda 0.234955
Smoothing (Alpha): 0.015789

Модель SARIMA с недельной сезонностью

Сравним ТБАТС с другим широко применяемым и широко известным методом: САРИМА. SARIMA доказала, что предоставляет самые современные решения для прогнозирования временных рядов. К сожалению, у него есть два основных недостатка: (1) можно смоделировать только один сезонный эффект, (2) продолжительность сезона не должна быть слишком большой.

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

from pmdarima import auto_arima
arima_model = auto_arima(y_to_train, seasonal=True, m=7)
y_arima_forecast = arima_model.predict(n_periods=365)

Auto arima выбрала модель SARIMA(0, 1, 1)x(1, 0, 1, 7). Поскольку ожидаемая годовая картина не моделируется (см. рис. 4).

Рис. 4: SARIMA моделирует только недельные паттерны. Сравните с рис. 2.

SARIMAX с условиями Фурье

Можно применить трюк [4], чтобы использовать экзогенные переменные в SARIMAX для моделирования дополнительных сезонностей с членами Фурье.

Мы продолжим моделировать недельную модель с помощью сезонной части SARIMA. Для годового сезонного шаблона мы будем использовать вышеупомянутый прием. Я сравнил несколько вариантов количества членов Фурье, и 2 дает наиболее точные прогнозы. Поэтому мы будем использовать 2 члена Фурье в качестве экзогенных переменных.

# prepare Fourier terms
exog = pd.DataFrame({'date': y.index})
exog = exog.set_index(pd.PeriodIndex(exog['date'], freq='D'))
exog['sin365'] = np.sin(2 * np.pi * exog.index.dayofyear / 365.25)
exog['cos365'] = np.cos(2 * np.pi * exog.index.dayofyear / 365.25)
exog['sin365_2'] = np.sin(4 * np.pi * exog.index.dayofyear / 365.25)
exog['cos365_2'] = np.cos(4 * np.pi * exog.index.dayofyear / 365.25)
exog = exog.drop(columns=['date'])
exog_to_train = exog.iloc[:(len(y)-365)]
exog_to_test = exog.iloc[(len(y)-365):]# Fit model
arima_exog_model = auto_arima(y=y_to_train, exogenous=exog_to_train, seasonal=True, m=7)# Forecast
y_arima_exog_forecast = arima_exog_model.predict(n_periods=365, exogenous=exog_to_test)

С помощью членов Фурье SARIMAX может моделировать обе сезонные модели (рис. 5).

Рис. 5: SARIMAX с моделями членов Фурье как для недельных, так и для годовых моделей

Сравнение моделей

До сих пор мы рассматривали только участки. Давайте сравним производительность моделей, используя прогноз на 365 дней вперед. Мы будем использовать среднюю абсолютную ошибку в качестве нашей метрики:

ТБАТ: 3,8577
САРИМА: 7,2249
SARIMAX с 2 членами Фурье: 3,9045

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

Преимущества

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

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

Недостатки

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

В отличие от SARIMAX, BATS и TBATS не позволяют добавлять в модель экзогенные переменные для улучшения прогнозов. По словам Роба Дж. Хайндмана, включить их и не попасть в вопросы предсказуемости не так уж и просто.

Модель TBATS (модель пространства состояний экспоненциального сглаживания с преобразованием Бокса-Кокса, ошибками ARMA, трендовыми и сезонными компонентами)

Соответствует модели TBATS, примененной к y, как описано в De Livera, Hyndman & Snyder (2011). Параллельная обработка используется по умолчанию для ускорения вычислений.

tbats(
  y,
  use.box.cox = NULL,
  use.trend = NULL,
  use.damped.trend = NULL,
  seasonal.periods = NULL,
  use.arma.errors = TRUE,
  use.parallel = length(y) > 1000,
  num.cores = 2,
  bc.lower = 0,
  bc.upper = 1,
  biasadj = FALSE,
  model = NULL,
  ...
)

Заключение

TBATS позволяет пользователям легко обрабатывать данные с несколькими сезонными моделями. Эта модель предпочтительнее, когда сезонность меняется во времени.