Извлечение сигналов из моделей RNN, ARIMA и Prophet для прогнозирования с помощью Catboost

Исследования мощных моделей временных рядов надежны. Доступно множество вариантов, выходящих далеко за рамки классических методов, таких как ARIMA. В последнее время рекуррентные нейронные сети и LSTM стали предметом интереса многих исследователей. Судя по количеству загрузок, перечисленных на PePy, модель Prophet, вероятно, является наиболее широко используемой моделью прогнозистами временных рядов из-за ее низкого порога входа. Какой из этих вариантов лучше всего подходит для вашей серии?

Возможно, ответ заключается в том, что вы должны попробовать их все и объединить их сильные стороны. Обычный метод для этого известен как штабелирование. Популярная библиотека машинного обучения scikit-learn предлагает StackingRegressor, который можно использовать для задач временных рядов. Я ранее продемонстрировал, как его использовать именно для такого случая.

Однако у StackingRegressor есть ограничение; он принимает только другие классы моделей scikit-learn и API. Таким образом, такая модель, как ARIMA, недоступная в scikit-learn, или глубокие сети в tensorflow, становится недоступной для размещения в стеке. Есть решение. В этом посте я покажу, как расширить методы стекирования для временных рядов с помощью пакета scalecast и блокнота Jupyter. Используемый набор данных доступен в открытом доступе на GitHub. Необходимы следующие требования:

pip install --upgrade scalecast
conda install tensorflow
conda install shap
conda install -c conda-forge cmdstanpy
pip install prophet

Обзор набора данных

Набор данных ежечасный и разделен на набор поездов (700 наблюдений) и тестовый набор (48 наблюдений). В моей записной книжке используется серия H1, но модифицировать ее для использования любой почасовой серии из набора данных M4 несложно. Вот как можно прочитать данные и сохранить их в объекте Forecaster:

import pandas as pd
import numpy as np
from scalecast.Forecaster import Forecaster
from scalecast.util import metrics
import matplotlib.pyplot as plt
import seaborn as sns

def read_data(idx = 'H1', cis = True, metrics = ['smape']):
    info = pd.read_csv(
        'M4-info.csv',
        index_col=0,
        parse_dates=['StartingDate'],
        dayfirst=True,
    )
    train = pd.read_csv(
        f'Hourly-train.csv',
        index_col=0,
    ).loc[idx]
    test = pd.read_csv(
        f'Hourly-test.csv',
        index_col=0,
    ).loc[idx]
    y = train.values
    sd = info.loc[idx,'StartingDate']
    fcst_horizon = info.loc[idx,'Horizon']
    cd = pd.date_range(
        start = sd,
        freq = 'H',
        periods = len(y),
    )
    f = Forecaster(
        y = y, # observed values
        current_dates = cd, # current dates
        future_dates = fcst_horizon, # forecast length
        test_length = fcst_horizon, # test-set length
        cis = cis, # whether to evaluate intervals for each model
        metrics = metrics, # what metrics to evaluate
    )
    
    return f, test.values

f, test_set = read_data()
f # display the Forecaster object

Вот сюжет сериала:

Прикладные модели

Прежде чем мы начнем размещать модели, нам нужно сгенерировать на их основе прогнозы. Я решил использовать наивную модель в качестве эталона для моделей ARIMA, LSTM и Prophet. В последующих разделах я расскажу, почему я сделал тот или иной выбор, но можно было бы принять и другие решения, не менее интересные или даже более интересные.

Наивный тест

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

f.set_estimator('naive')
f.manual_forecast(seasonal=True)

АРИМА

Авторегрессионное интегрированное скользящее среднее — это популярный и простой метод временных рядов, который использует задержки и ошибки ряда для линейного прогнозирования его будущего. С помощью исследовательского анализа данных (который вы можете увидеть в связанной записной книжке) я определил, что этот ряд не был стационарным и что он был сильно сезонным. В конечном итоге я решил применить сезонную модель ARIMA порядка (5,1,4) x (1,1,1,24).

f.set_estimator('arima')
f.manual_forecast(
    order = (5,1,4),
    seasonal_order = (1,1,1,24),
    call_me = 'manual_arima',
)

ЛСТМ

Если ARIMA находится на более простой стороне моделей временных рядов, LSTM является одним из более продвинутых методов. Это метод глубокого обучения со многими параметрами, включая механизм внимания, который находит долгосрочные и краткосрочные закономерности в последовательных данных, что теоретически делает его идеальным выбором для временных рядов. Настроить эту модель с помощью tensorflow сложно, но не так уж и плохо в scalecast (см. эту статью). Я применил две модели LSTM: одну с функцией активации Tanh и одну с ReLu.

f.set_estimator('rnn')
f.manual_forecast(
    lags = 48,
    layers_struct=[
        ('LSTM',{'units':100,'activation':'tanh'}),
        ('LSTM',{'units':100,'activation':'tanh'}),
        ('LSTM',{'units':100,'activation':'tanh'}),
    ],
    optimizer = 'Adam',
    epochs = 15,
    plot_loss = True,
    validation_split=0.2,
    call_me = 'rnn_tanh_activation',
)

f.manual_forecast(
    lags = 48,
    layers_struct=[
        ('LSTM',{'units':100,'activation':'relu'}),
        ('LSTM',{'units':100,'activation':'relu'}),
        ('LSTM',{'units':100,'activation':'relu'}),
    ],
    optimizer = 'Adam',
    epochs = 15,
    plot_loss = True,
    validation_split=0.2,
    call_me = 'rnn_relu_activation',
)

Пророк

Несмотря на свою чрезвычайную популярность, модель Prophet в последнее время подвергается оклевете. Есть утверждения, что его точность не впечатляет, в основном из-за того, насколько нереалистично он экстраполирует тенденции и что он не учитывает локальные закономерности посредством авторегрессионного моделирования. Тем не менее, он делает некоторые вещи уникальными. Во-первых, он автоматически применяет праздничные эффекты к моделям. Также учитывается несколько видов сезонности. Он делает все это с минимальной спецификацией, требуемой от пользователя. Мне нравится идея использовать его в качестве сигнала, даже если он не идеален для создания точечных прогнозов.

f.set_estimator('prophet')
f.manual_forecast()

Сравните результаты

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

results = f.export(determine_best_by='TestSetSMAPE')
ms = results['model_summaries']
ms[
    [
        'ModelNickname',
        'TestSetLength',
        'TestSetSMAPE',
        'InSampleSMAPE',
    ]
]

Хорошей новостью здесь является то, что каждая модель превзошла наивный метод. Модель ARIMA показала лучшие результаты с процентной ошибкой 4,7%, за ней следует Prophet. Давайте посмотрим, как построены все прогнозы относительно проверочного набора:

f.plot(order_by="TestSetSMAPE",ci=True)
plt.show()

Все эти модели показали хорошие результаты в этой серии, без больших расхождений между ними. Складываем их!

Модели стекирования

Каждой модели суммирования требуется окончательный оценщик, который будет фильтровать различные оценки других моделей, чтобы создать новый набор прогнозов. Мы объединим наши ранее изученные результаты с оценщиком расширенного дерева под названием Catboost. Catboost — это мощная процедура, которая, как мы надеемся, позволит воплотить в жизнь лучшие сигналы каждой из уже примененных моделей.

f.add_signals(
    f.history.keys(), # add signals from all previously evaluated models
)
f.add_ar_terms(48)
f.set_estimator('catboost')

Приведенный выше код добавляет прогнозы каждой из оцененных моделей в объект Forecaster. Он называет эти предсказания «сигналами». Они обрабатываются так же, как и любые другие ковариаты, хранящиеся в этом же объекте. Мы также добавляем лаги последних 48 рядов в качестве дополнительных регрессоров, которые модель Catboost может использовать для прогнозирования. Давайте теперь вызовем три модели Catboost: 1, которая использует все доступные сигналы и запаздывания последовательности, одна, которая использует только сигналы, и третья, которая использует только запаздывания.

f.manual_forecast(
    Xvars='all',
    call_me='catboost_all_reg',
    verbose = False,
)
f.manual_forecast(
    Xvars=[x for x in f.get_regressor_names() if x.startswith('AR')], 
    call_me = 'catboost_lags_only',
    verbose = False,
)
f.manual_forecast(
    Xvars=[x for x in f.get_regressor_names() if not x.startswith('AR')], 
    call_me = 'catboost_signals_only',
    verbose = False,
)

Давайте подключимся к этому тестовому набору, который мы держали отдельно от объекта Forecaster в начале анализа, и сравним результаты всех примененных моделей. На этот раз мы рассмотрим два показателя: SMAPE и среднюю абсолютную масштабированную ошибку (MASE). Это две метрики, которые использовались в реальном соревновании M4.

test_results = pd.DataFrame(index = f.history.keys(),columns = ['smape','mase'])
for k, v in f.history.items():
    test_results.loc[k,['smape','mase']] = [
        metrics.smape(test_set,v['Forecast']),
        metrics.mase(test_set,v['Forecast'],m=24,obs=f.y),
    ]
    
test_results.sort_values('smape')

Объединив сигналы из разных классов моделей, мы создали два оценщика, которые превосходят другие: модель Catboost, которая была обучена с использованием всех сигналов и лагов 48 серий, за которой следует модель Catboost, которая только использовала сигналы. Оба они получили ошибку около 2,8% от выборки. Мы можем видеть эти две модели, построенные с фактическими значениями из тестового набора.

fig, ax = plt.subplots(figsize=(12,6))
f.plot(
    models = ['catboost_all_reg','catboost_signals_only'],
    ci=True,
    ax = ax
)
sns.lineplot(
    x = f.future_dates, 
    y = test_set, 
    ax = ax,
    label = 'held out actuals',
    color = 'darkblue',
    alpha = .75,
)
plt.show()

Какие сигналы были наиболее важными?

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

f.export_feature_importance('catboost_all_reg')

На приведенном выше снимке экрана показаны только первые несколько наиболее важных предикторов, но из этого мы можем видеть, что сигнал ARIMA был наиболее важным, за ним следовал первый лаг ряда, а затем сигнал пророка. Модели RNN также получили более высокие оценки, чем многие включенные лаги. Это может быть хорошей отправной точкой, если мы хотим обучить более экономичную модель в будущем.

Заключение

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