Временные ряды: от ARIMA до градиентного бустинга

Почти всё, что мы строили в этом треке, опиралось на одно тихое допущение: наблюдения независимы и одинаково распределены (i.i.d.). Именно оно оправдывает случайное перемешивание перед разбиением на train/test, именно оно делает кросс-валидацию несмещённой оценкой риска, именно на нём стоят теоремы обобщения из статьи про переобучение и регуляризацию.

Во временных рядах это допущение неверно дважды. Во-первых, наблюдения зависимы: продажи сегодня коррелируют с продажами вчера. Во-вторых, распределение меняется во времени: инфляция, рост базы пользователей, смена продукта. Поэтому наивный перенос ML-рецептов на ряды даёт катастрофу конкретного, узнаваемого вида — модель показывает R² = 0.98 на кросс-валидации и разваливается в первую же неделю продакшена.

Эта статья про то, как прогнозировать честно. Мы пройдём путь от классической статистики (ARIMA, ETS), которая моделирует процесс порождения ряда, до современного признакового подхода, который сводит прогнозирование к обычной регрессии и отдаёт её градиентному бустингу — именно этот подход выигрывает большинство прикладных задач и почти все соревнования.

Постановка задачи

Временной ряд — последовательность наблюдений y₁, y₂, …, y_T, упорядоченных по времени и (обычно) равномерно расположенных с шагом Δt: час, день, неделя, месяц.

Задача прогнозирования: по истории до момента T предсказать значения на горизонте h:

ŷ_{T+1}, ŷ_{T+2}, …, ŷ_{T+h}  =  f(y₁…y_T, X₁…X_{T+h})

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

Несколько развилок, которые нужно осознать до выбора алгоритма:

Последняя строка важнее всех остальных. Прогноз почти никогда не является конечным продуктом: он питает решение — сколько заказать товара, сколько поставить смен, сколько зарезервировать мощности. Стоимость недопрогноза и перепрогноза почти всегда асимметрична, и это должно попадать в функцию потерь, а не в устные договорённости.

Эволюция подходов

Полезный вывод из истории: сложность не выигрывает автоматически. Соревнования M3 и M4 (Makridakis et al., The M4 Competition: Results, findings, conclusion) показали, что чистые нейросети регулярно проигрывали простому экспоненциальному сглаживанию. Перелом наступил в M5, где данных было много, рядов десятки тысяч, и глобальные модели на градиентном бустинге получили достаточно материала, чтобы учиться поперёк рядов.

Компоненты ряда и декомпозиция

Классическая оптика смотрит на ряд как на сумму (или произведение) нескольких компонент:

• Тренд T(t) — медленное направленное изменение уровня.
• Сезонность S(t) — колебания фиксированного периода: недельный профиль, годовой профиль, суточный профиль.
• Цикл — колебания непостоянного периода (деловой цикл); часто сливают с трендом.
• Остаток R(t) — то, что не объяснено: шум плюс аномалии.

Две базовые модели:

аддитивная:        y(t) = T(t) + S(t) + R(t)
мультипликативная: y(t) = T(t) · S(t) · R(t)

Выбирать просто: если амплитуда сезонных колебаний растёт вместе с уровнем ряда (декабрьский всплеск продаж в абсолютных рублях больше, когда бизнес крупнее) — модель мультипликативная. Логарифмирование переводит мультипликативную в аддитивную:

log y = log T + log S + log R

Это первая причина, по которой ряды почти всегда стоит попробовать в логах. Вторая — стабилизация дисперсии (частный случай преобразования Бокса–Кокса). Третья — прогноз в логах автоматически неотрицателен после обратного преобразования. Цена: exp(E[log y]) ≠ E[y], при обратном переходе нужна поправка Даффи (·exp(σ²/2) в предположении логнормальности), иначе прогноз систематически занижен.

STL (Seasonal-Trend decomposition using Loess, Cleveland et al., 1990) — рабочая лошадка декомпозиции. В отличие от наивного скользящего среднего, STL допускает медленное изменение сезонного профиля, устойчив к выбросам (robust-режим) и работает с любым периодом.

import pandas as pd
from statsmodels.tsa.seasonal import STL

# ряд с датой в индексе и явно заданной частотой — без этого statsmodels не поймёт период
s = pd.read_csv("sales.csv", parse_dates=["date"], index_col="date").asfreq("D")["y"]

stl = STL(s, period=7, seasonal=13, robust=True)   # period=7: недельная сезонность
res = stl.fit()

trend, season, resid = res.trend, res.seasonal, res.resid

# практическое применение №1: детекция аномалий по остатку через устойчивый z-score
mad = (resid - resid.median()).abs().median() * 1.4826    # робастная оценка σ
anomalies = resid[(resid - resid.median()).abs() > 4 * mad]

# практическое применение №2: сезонная корректировка ряда перед моделированием
deseasonalized = s - season

Практический смысл декомпозиции — не красивая картинка, а три вещи: (1) детекция аномалий по остатку, (2) сезонная корректировка перед подачей в модель, (3) визуальная проверка гипотез до того, как вы потратили день на подбор гиперпараметров.

Стационарность: почему без неё классика не работает

Ряд строго стационарен, если совместное распределение любого набора наблюдений не меняется при сдвиге по времени. На практике пользуются слабой стационарностью: постоянное матожидание, постоянная дисперсия и автоковариация, зависящая только от лага, а не от абсолютного времени.

Зачем это нужно? Модели семейства ARMA оценивают конечный набор параметров по одной реализации процесса. Если распределение дрейфует, «средний» параметр не описывает ничего. Ещё опаснее эффект ложной регрессии (spurious regression): два независимых случайных блуждания дают R² около 0.7 и «значимые» коэффициенты просто потому, что оба содержат стохастический тренд (Granger & Newbold, 1974). Знаменитые графики «корреляция потребления маргарина и разводов в Мэне» — ровно этот артефакт.

Как проверять

Два теста с противоположными нулевыми гипотезами — и применять их надо в паре:

from statsmodels.tsa.stattools import adfuller, kpss

def stationarity_report(x, name=""):
    adf_p = adfuller(x.dropna(), autolag="AIC")[1]
    kpss_p = kpss(x.dropna(), regression="c", nlags="auto")[1]
    # ADF: H0 = единичный корень; KPSS: H0 = стационарность
    if adf_p < 0.05 and kpss_p > 0.05:
        verdict = "стационарен"
    elif adf_p >= 0.05 and kpss_p <= 0.05:
        verdict = "нестационарен — нужно дифференцирование"
    else:
        verdict = "тесты противоречат: возможен тренд-стационарный ряд или мало данных"
    print(f"{name}: ADF p={adf_p:.4f}, KPSS p={kpss_p:.4f} -> {verdict}")

stationarity_report(s, "исходный")
stationarity_report(s.diff(), "после d=1")
stationarity_report(s.diff().diff(7), "после d=1, D=1 (период 7)")

Противоречие тестов — не баг, а полезный сигнал: ряд может быть стационарен вокруг детерминированного тренда (тогда лечится вычитанием тренда, а не дифференцированием) или наоборот.

Дифференцирование

∇y_t = y_t − y_{t−1} убирает линейный тренд, ∇² — квадратичный, сезонная разность ∇_m y_t = y_t − y_{t−m} убирает сезонность периода m.

Главная ошибка новичка — передифференцировать. Каждая разность добавляет шум и MA-компоненту, а лишняя разность делает ряд «переразностным» с искусственной отрицательной автокорреляцией на лаге 1. Правило: d почти никогда не больше 2, D почти никогда не больше 1. В pmdarima/statsmodels порядок подбирают тестами (ndiffs, nsdiffs) вместо ручного перебора.

ACF и PACF — язык, на котором ряд говорит о себе

Автокорреляционная функция ρ(k) = corr(y_t, y_{t−k}) показывает суммарную связь с лагом k, включая опосредованную. Частная автокорреляция φ(k) — связь с лагом k после исключения влияния всех промежуточных лагов.

Аналогия: если A влияет на B, а B на C, то ACF увидит связь A–C, а PACF её обнулит, показав только прямое влияние. Это позволяет различать AR и MA:

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, ax = plt.subplots(2, 1, figsize=(11, 6))
plot_acf(s.diff().dropna(), lags=40, ax=ax[0])    # синяя зона = доверительный коридор ±1.96/√n
plot_pacf(s.diff().dropna(), lags=40, ax=ax[1], method="ywm")

Значимость определяется выходом за коридор ±1.96/√n. При 40 лагах примерно два всплеска окажутся «значимыми» случайно — это множественное тестирование, не приписывайте смысл каждому.

Базовые модели, которые обязательно нужно побить

Прежде чем обучать что-либо, зафиксируйте baseline. Без него невозможно понять, что вы вообще делаете полезного.

import numpy as np

def naive(y, h):                     # прогноз = последнее значение
    return np.repeat(y[-1], h)

def seasonal_naive(y, h, m):         # прогноз = значение из прошлого периода
    return np.array([y[-m + (i % m)] for i in range(h)])

def drift(y, h):                     # линейная экстраполяция первой и последней точки
    slope = (y[-1] - y[0]) / (len(y) - 1)
    return y[-1] + slope * np.arange(1, h + 1)

Сезонный наивный прогноз на рядах с сильной недельной или годовой сезонностью — на удивление сильный соперник. На дневных розничных данных он часто даёт MAPE 15–20%, и если ваш LSTM показывает 18%, работа была проделана впустую. В M4 значительная доля участников не побила комбинацию простых методов.

Метрики: почему MAPE врёт

Разбор метрик регрессии есть в статье про оценку моделей, но у рядов своя специфика.

Про MAPE стоит сказать отдельно, потому что это метрика по умолчанию в большинстве бизнес-отчётов. Проблемы:

1. Деление на ноль. Товар не продавался в понедельник — MAPE не определён. Интермиттентный спрос (много нулей) делает её бессмысленной целиком.
2. Асимметрия. Занижение прогноза ограничено 100% ошибки, завышение — не ограничено ничем. Оптимизация MAPE систематически толкает прогноз вниз.
3. Несопоставимость. Ошибка 20% на товаре с продажами 3 шт./день и на товаре с 3000 шт./день — это разные по стоимости ошибки, а метрика их уравнивает.

MASE (Hyndman & Koehler, 2006) лишена этих проблем: она нормирована на ошибку сезонного наивного прогноза, безразмерна, определена при нулях и напрямую интерпретируется — MASE < 1 значит «лучше наивного».

def mase(y_true, y_pred, y_train, m=1):
    """Mean Absolute Scaled Error. m — период сезонности (1 = обычный naive)."""
    scale = np.mean(np.abs(y_train[m:] - y_train[:-m]))   # ошибка наивного на обучении
    if scale == 0:
        return np.nan                                     # константный ряд
    return np.mean(np.abs(y_true - y_pred)) / scale


def pinball_loss(y_true, y_pred, q):
    """Квантильная (pinball) потеря: асимметричный штраф для q-квантиля."""
    d = y_true - y_pred
    return np.mean(np.maximum(q * d, (q - 1) * d))

Pinball loss — то, что нужно, когда стоимость ошибки асимметрична. Дефицит товара стоит вам маржи и лояльности, избыток — заморозки оборотных средств и списания. Если недостача обходится втрое дороже излишка, вам нужен не средний прогноз, а квантиль уровня q = 3/(3+1) = 0.75. Это классическая задача газетчика (newsvendor), и она решается прямо на уровне функции потерь: LGBMRegressor(objective="quantile", alpha=0.75).

Честная валидация: rolling origin

Обычная k-fold кросс-валидация на рядах — прямая утечка из будущего. Модель, обученная на данных за март, предсказывающая февраль, «знает» и общий уровень, и произошедшие события. Оценка получится оптимистичной в разы.

Правильная схема — бэктест с плавающей точкой отсечения (rolling origin, walk-forward):

Три детали, которые обычно упускают:

1. Gap (эмбарго). Если вы прогнозируете на 14 дней вперёд, между концом обучения и началом теста нужен разрыв в 14 дней — иначе модель тестируется в условиях, недоступных в проде, где последние две недели фактов ещё не пришли. Отдельная тонкость: лаговые признаки, построенные по краю обучающего окна, «протекают» в тест, если gap не заложен.
2. Expanding vs sliding. Расширяющееся окно использует всю историю — хорошо при стабильном процессе. Скользящее окно фиксированной длины лучше при дрейфе распределения (сменился ассортимент, поменялась цена) и делает время обучения предсказуемым.
3. Число фолдов и агрегация. Одна точка отсечения — это одно наблюдение метрики с огромной дисперсией. Нужно 5–10 фолдов, и агрегировать не только среднее, но и разброс.

import numpy as np
import pandas as pd

def rolling_origin_splits(n, initial, horizon, step, gap=0, max_train=None):
    """Генератор индексов (train_idx, test_idx) для бэктеста.

    initial   — минимальный размер обучения
    horizon   — длина тестового окна = горизонт прогноза
    step      — на сколько сдвигается точка отсечения между фолдами
    gap       — эмбарго между train и test (обычно = horizon)
    max_train — если задан, окно скользящее фиксированной длины
    """
    cutoff = initial
    while cutoff + gap + horizon <= n:
        start = 0 if max_train is None else max(0, cutoff - max_train)
        yield np.arange(start, cutoff), np.arange(cutoff + gap, cutoff + gap + horizon)
        cutoff += step


def backtest(df, model_factory, feature_cols, target="y", **split_kw):
    """Возвращает метрики по фолдам. model_factory() создаёт свежую необученную модель."""
    scores = []
    for k, (tr, te) in enumerate(rolling_origin_splits(len(df), **split_kw)):
        model = model_factory()
        model.fit(df.iloc[tr][feature_cols], df.iloc[tr][target])
        pred = model.predict(df.iloc[te][feature_cols])
        scores.append({
            "fold": k,
            "mae": np.mean(np.abs(df.iloc[te][target] - pred)),
            "mase": mase(df.iloc[te][target].values, pred, df.iloc[tr][target].values, m=7),
        })
    out = pd.DataFrame(scores)
    print(out, "\nMASE: %.3f ± %.3f" % (out.mase.mean(), out.mase.std()))
    return out

В scikit-learn есть готовый TimeSeriesSplit, но у него нет параметра gap до версии 0.24 и нет ограничения на длину обучающего окна в старых версиях — проверяйте, что именно вы получаете. Общая логика утечек разобрана в статье про данные и признаки; во временных рядах утечка — не риск, а состояние по умолчанию, которое надо активно предотвращать.

Диагностический цикл Бокса–Дженкинса

Обратите внимание на шаг с проверкой остатков. Если в остатках ARIMA осталась автокорреляция (тест Люнга–Бокса отвергает гипотезу о белом шуме), модель не выжала из ряда всю структуру — и, что важнее, её доверительные интервалы неверны.

ARIMA и SARIMAX

AR(p): текущее значение — линейная комбинация p прошлых значений.

y_t = c + φ₁y_{t−1} + … + φ_p y_{t−p} + ε_t

MA(q): текущее значение — линейная комбинация q прошлых ошибок (шоков).

y_t = c + ε_t + θ₁ε_{t−1} + … + θ_q ε_{t−q}

Разница содержательная. AR — это инерция: высокие продажи вчера тянут за собой сегодняшние. MA — это затухающее эхо разовых шоков: сломался склад, всплеск ошибки прогноза, и последствия расходятся ещё несколько дней.

ARIMA(p, d, q) = ARMA, применённая к ряду после d дифференцирований. SARIMA(p,d,q)(P,D,Q)_m добавляет те же компоненты на сезонных лагах. SARIMAX — плюс экзогенные регрессоры.

import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX

# экзогенные признаки должны быть известны и на горизонте прогноза
exog = pd.DataFrame({
    "promo":   promo_flag,          # запланированные акции — известны заранее
    "holiday": holiday_flag,        # производственный календарь — известен заранее
}, index=s.index)

model = SARIMAX(
    s,
    exog=exog,
    order=(1, 1, 1),                # (p, d, q)
    seasonal_order=(1, 1, 1, 7),    # (P, D, Q, m) — недельная сезонность
    enforce_stationarity=False,     # оставить оценке свободу; проверять остатки вручную
    enforce_invertibility=False,
)
fit = model.fit(disp=False)
print(fit.summary())

# диагностика остатков: они обязаны быть белым шумом
from statsmodels.stats.diagnostic import acorr_ljungbox
lb = acorr_ljungbox(fit.resid, lags=[7, 14, 21], return_df=True)
print(lb)                            # p-value > 0.05 => автокорреляция не обнаружена

fc = fit.get_forecast(steps=14, exog=future_exog)
point = fc.predicted_mean
ci80 = fc.conf_int(alpha=0.20)       # 80%-й интервал

Автоподбор порядков. Перебор по сетке с критерием AICc — стандарт де-факто; алгоритм Хиндмана–Хандакара из пакета forecast реализован в Python как pmdarima.auto_arima (документация).

import pmdarima as pm

auto = pm.auto_arima(
    s, exogenous=exog, m=7,
    seasonal=True, d=None, D=None,      # порядки разностей определить тестами
    stepwise=True,                      # пошаговый поиск вместо полного перебора
    information_criterion="aicc",
    suppress_warnings=True, error_action="ignore",
)
print(auto.summary())

Сложность. Оценка методом максимального правдоподобия через фильтр Калмана стоит O(T · k³) на итерацию, где k — размерность вектора состояния (порядка p + q + m·(P+Q)). Сезонный член m = 365 для дневных данных с годовой сезонностью раздувает k до сотен — SARIMA на дневных данных с годовым периодом практически неприменима. Именно для этого случая используют Fourier-члены как экзогенные регрессоры (ARIMA с гармониками): K пар синус/косинус вместо m сезонных параметров.

# годовая сезонность через гармоники: 4 пары вместо 365 сезонных лагов
import numpy as np
t = np.arange(len(s))
fourier = pd.DataFrame({
    f"{fn.__name__}_{k}": fn(2 * np.pi * k * t / 365.25)
    for k in range(1, 5) for fn in (np.sin, np.cos)
}, index=s.index)

Когда ARIMA всё ещё лучший выбор: один ряд, короткая история (50–200 точек), нужны статистически обоснованные доверительные интервалы, требуется объяснимость перед регулятором. При тысячах рядов и богатых признаках она проигрывает.

ETS и экспоненциальное сглаживание

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

простое сглаживание:  ℓ_t = α·y_t + (1 − α)·ℓ_{t−1},   ŷ_{t+h} = ℓ_t

Метод Хольта добавляет компоненту тренда b_t, Хольта–Уинтерса — сезонную s_t. Обобщение — таксономия ETS (Error, Trend, Seasonal), где каждая компонента бывает аддитивной, мультипликативной, демпфированной или отсутствующей: 30 комбинаций, из которых модель выбирается по AICc автоматически.

from statsmodels.tsa.holtwinters import ExponentialSmoothing

hw = ExponentialSmoothing(
    s,
    trend="add", damped_trend=True,   # демпфирование почти всегда улучшает длинный горизонт
    seasonal="add", seasonal_periods=7,
    initialization_method="estimated",
).fit()
pred = hw.forecast(14)

Демпфированный тренд (damped_trend=True) — недооценённый практический приём. Линейный тренд, экстраполированный на 90 дней вперёд, даёт абсурдные значения; демпфирование делает его затухающим и почти всегда снижает ошибку на длинном горизонте. В M3 демпфированный Хольт был одним из лучших методов в целом.

Theta-метод — победитель M3, эквивалентный простому сглаживанию с дрейфом, реализован как statsmodels.tsa.forecasting.theta.ThetaModel. Это отличный «сильный baseline» стоимостью в три строки.

Prophet: инструмент аналитика

Prophet от Meta — не статистическая модель ряда, а аддитивная регрессия по времени:

y(t) = g(t) + s(t) + h(t) + ε

где g — кусочно-линейный (или логистический) тренд с автоматически найденными точками перелома, s — сезонности через ряды Фурье, h — эффекты праздников с окнами до и после.

Сильные стороны: устойчивость к пропускам, встроенный календарь праздников по странам, интерпретируемые компоненты, работа «из коробки» без подбора. Слабые: это по сути сглаженная кривая по времени, она плохо использует автокорреляцию малых лагов, склонна переобучаться на точках перелома тренда и в бенчмарках регулярно проигрывает как ETS, так и бустингу (критический разбор Хиндмана). Разумная позиция: Prophet — быстрый способ получить приличный прогноз и красивую разбивку по компонентам для аналитического отчёта, но не финальная продакшн-модель для тысяч рядов.

Признаковый подход: сводим прогноз к регрессии

Ключевая идея, на которой стоит весь современный прикладной форкастинг: превратить ряд в обычную табличную задачу, где строка — это «объект в момент времени», а признаки построены только из информации, доступной на момент прогноза. Дальше работает любой табличный регрессор, и на практике это градиентный бустинг.

Что даёт этот подход:

• естественная работа с сотнями экзогенных признаков (цена, промо, погода, остатки на складе);
• глобальные модели — одна модель на все ряды сразу, обучающаяся переносить закономерности между похожими товарами; это спасает короткие ряды и новинки;
• нелинейности и взаимодействия без ручной спецификации;
• готовая инфраструктура ML: та же валидация, тот же мониторинг, тот же MLOps.

Инженерия признаков

import pandas as pd
import numpy as np

def make_features(df, horizon, lags=(7, 14, 28, 365), windows=(7, 28, 91)):
    """df: колонки [date, series_id, y, price, promo]. horizon — на сколько шагов вперёд.

    ГЛАВНОЕ ПРАВИЛО: любой признак, построенный из y, обязан отставать
    минимум на horizon шагов — в момент прогноза свежих фактов ещё нет.
    """
    df = df.sort_values(["series_id", "date"]).copy()
    g = df.groupby("series_id")["y"]

    # 1. лаги: прямая память ряда
    for L in lags:
        df[f"lag_{L}"] = g.shift(L)

    # 2. скользящие статистики — ОБЯЗАТЕЛЬНО от сдвинутого ряда, иначе утечка текущего y
    shifted = g.shift(horizon)
    for w in windows:
        r = shifted.rolling(w)
        df[f"roll_mean_{w}"] = r.mean()
        df[f"roll_std_{w}"] = r.std()
        df[f"roll_max_{w}"] = r.max()
    # экспоненциальное сглаживание как признак «текущего уровня»
    df["ewm_28"] = shifted.ewm(span=28, adjust=False).mean()

    # 3. календарь — известен на любую дату в будущем
    d = df["date"].dt
    df["dow"], df["dom"], df["month"], df["woy"] = d.dayofweek, d.day, d.month, d.isocalendar().week
    df["is_weekend"] = (df["dow"] >= 5).astype(int)
    df["days_to_month_end"] = d.days_in_month - d.day

    # 4. циклическое кодирование: декабрь и январь должны быть рядом
    df["dow_sin"] = np.sin(2 * np.pi * df["dow"] / 7)
    df["dow_cos"] = np.cos(2 * np.pi * df["dow"] / 7)
    df["moy_sin"] = np.sin(2 * np.pi * df["month"] / 12)
    df["moy_cos"] = np.cos(2 * np.pi * df["month"] / 12)

    # 5. события: расстояние до праздника важнее бинарного флага
    df["days_to_holiday"] = days_to_nearest_holiday(df["date"])

    # 6. экзогенные, известные на будущее
    df["price_rel"] = df["price"] / df.groupby("series_id")["price"].transform("median")
    df["promo_next_7"] = df.groupby("series_id")["promo"].transform(
        lambda x: x.shift(-1).rolling(7).max())      # будущее промо — это ЗАПЛАНИРОВАННОЕ промо

    # 7. счётчики интермиттентности
    df["days_since_sale"] = g.transform(
        lambda x: x.groupby((x > 0).cumsum()).cumcount())
    return df

Пункт 2 — самая частая и самая дорогая ошибка в проде. rolling(7).mean() без предварительного shift(horizon) включает в признак сегодняшнее значение целевой переменной. На бэктесте это даёт MAE лучше в 3–5 раз, в проде — обычный уровень наивного прогноза и неприятный разговор с заказчиком.

Обучение бустинга

import lightgbm as lgb

FEATURES = [c for c in df.columns if c not in ("date", "y", "series_id")]

params = dict(
    objective="tweedie",        # для спроса с нулями и правым хвостом; см. ниже
    tweedie_variance_power=1.1,
    learning_rate=0.03,
    num_leaves=127,
    min_data_in_leaf=100,
    feature_fraction=0.8,
    bagging_fraction=0.8, bagging_freq=1,
    lambda_l2=1.0,
    n_estimators=3000,
)

train = df[df.date < cutoff]
valid = df[(df.date >= cutoff) & (df.date < cutoff + pd.Timedelta(days=28))]

model = lgb.LGBMRegressor(**params)
model.fit(
    train[FEATURES], train["y"],
    eval_set=[(valid[FEATURES], valid["y"])],
    eval_metric="l1",
    categorical_feature=["series_id", "dow", "month"],
    callbacks=[lgb.early_stopping(200), lgb.log_evaluation(100)],
)

Про objective: для розничного спроса с массой нулей и редкими крупными продажами распределение Твиди (tweedie, p ≈ 1.1–1.3) заметно лучше квадратичной потери — это была одна из ключевых находок победителей M5. Для симметричных непрерывных рядов подойдёт l1 (оптимизирует медиану, устойчива к выбросам) или l2.

Осторожно с трендом. Деревья не экстраполируют: если ряд растёт, а в обучении не было значений выше 1000, модель никогда не предскажет 1200. Три лекарства: (1) прогнозировать не уровень, а отношение y_t / roll_mean_28 или разность; (2) вычесть тренд заранее и добавить обратно после; (3) гибрид — линейная модель на тренд плюс бустинг на остаток. Ограничение фундаментальное и в проде проявляется через несколько месяцев после запуска, когда бизнес вырос за пределы обучающего диапазона.

Рекурсивная или прямая стратегия

Практический компромисс, который используют чаще всего, — прямая стратегия с группировкой горизонтов: одна модель на дни 1–7, другая на 8–14, третья на 15–28. Плюс вариант «одна модель с признаком horizon» — дешевле по инфраструктуре, но требует, чтобы модель сама выучила зависимость точности от горизонта.

Сложность. Обучение LightGBM: O(n · f · log n) при построении гистограмм, где n — число строк (ряды × даты), f — число признаков; память O(n · f) в бинаризованном виде (обычно 1 байт на признак-значение). На панели 30 000 товаров × 1800 дней = 54 млн строк × 60 признаков это ~3–4 ГБ и десятки минут на CPU. Рекурсивный прогноз на h шагов требует h последовательных пересчётов признаков — на больших панелях это узкое место инференса, поэтому прямая стратегия в проде часто выигрывает не точностью, а временем ответа.

Интервалы предсказания

Точечный прогноз без интервала почти бесполезен для принятия решения: страховой запас определяется именно верхним квантилем.

Три рабочих способа:

1. Аналитические интервалы ARIMA/ETS — выводятся из предположения о нормальности остатков. Известно, что они систематически слишком узкие: не учитывают неопределённость самих оценённых параметров и ошибку спецификации модели.
2. Квантильная регрессия. Обучить бустинг с objective="quantile" для нужных уровней (0.1, 0.5, 0.9). Просто, работает, но квантили могут пересекаться — их надо сортировать постобработкой.
3. Conformal prediction. Модельно-агностичный подход с гарантией покрытия: берём абсолютные ошибки на калибровочном отрезке, берём их (1−α)-квантиль, строим интервал ŷ ± q. Для рядов существует адаптивная версия EnbPI (Xu & Xie, 2021), корректирующая ширину под наблюдаемое покрытие.

def conformal_interval(model, X_cal, y_cal, X_new, alpha=0.1):
    """Простейший split-conformal: гарантия покрытия 1-alpha при обмениваемости.
    Для рядов калибровочная выборка берётся ПОЗЖЕ обучающей (out-of-time)."""
    residuals = np.abs(y_cal - model.predict(X_cal))
    n = len(residuals)
    q = np.quantile(residuals, np.ceil((n + 1) * (1 - alpha)) / n, method="higher")
    pred = model.predict(X_new)
    return pred - q, pred, pred + q

Проверять интервалы надо не по «выглядит правдоподобно», а по эмпирическому покрытию: доля фактов, попавших внутрь 80%-го интервала, должна быть около 80%. Если 55% — интервалы не работают, и построенный на них страховой запас систематически недостаточен.

Иерархические ряды и согласование

Реальный прогноз почти всегда иерархичен: товар → категория → магазин → регион → компания. Наивное решение — прогнозировать каждый уровень отдельно — даёт несогласованные числа: сумма прогнозов по товарам не равна прогнозу по категории. Финансовый директор такое не примет.

Подходы:

• Bottom-up — прогнозировать нижний уровень и суммировать. Согласовано по построению, но нижние ряды самые шумные.
• Top-down — прогнозировать агрегат и распределять по историческим долям. Стабильно наверху, теряет специфику отдельных товаров.
• Optimal reconciliation (MinT) — прогнозировать все уровни независимо, затем спроецировать вектор прогнозов на пространство согласованных решений с минимальной дисперсией ошибки (Wickramasuriya et al., 2019). На практике даёт лучший результат и реализован в hierarchicalforecast от Nixtla.

Отдельная категория — интермиттентный спрос: запчасти, редкие SKU, где большинство дней ноль. Классика здесь — метод Кростона и его вариант SBA, разделяющие ряд на «размер продажи» и «интервал между продажами». В ML-подходе то же самое достигается признаком days_since_sale и Твиди-целью.

Выбор подхода

Практический порядок действий, который экономит недели:

1. Построить график ряда и посмотреть на него глазами. Серьёзно — половина проблем (сдвиг единиц измерения, дыры в данных, смена системы учёта) видна сразу.
2. Посчитать seasonal_naive и зафиксировать его MASE как планку.
3. Прогнать ETS/Theta через statsforecast — это минуты и часто уже достаточная точность.
4. Если рядов много и есть экзогенные факторы — собрать признаковый датасет и обучить LightGBM.
5. Нейросети — только если пункт 4 упёрся в потолок и есть десятки тысяч рядов.

Про пункт 5: современные глобальные архитектуры (DeepAR, N-BEATS, Temporal Fusion Transformer, PatchTST) действительно выигрывают на больших однородных панелях с длинным контекстом, особенно когда нужен вероятностный прогноз распределения целиком. Устройство рекуррентных сетей и внимания разбирается в треке про нейронные сети; здесь важно понимать, что порог входа по данным и инфраструктуре высокий, а выигрыш над хорошо настроенным бустингом обычно измеряется единицами процентов.

Типичные ошибки

1. Перемешивание перед разбиением. Random split вместо временного — оценка завышена в разы.
2. Скользящие статистики без сдвига. rolling(7).mean() включает текущее y.
3. Признаки, недоступные в момент прогноза. Фактическая погода, фактические остатки на складе, курс на день прогноза. Правило: для каждого признака ответьте, во сколько именно он станет известен.
4. Ретроспективно пересчитанные данные. Цены и справочники в хранилище обновляются задним числом; на бэктесте вы видите исправленную версию, в проде — сырую. Лечится snapshot-таблицами с valid_from/valid_to.
5. Нет gap между train и test. Модель тестируется в невозможных условиях.
6. MAPE на рядах с нулями. Метрика не определена или взрывается.
7. Оптимизация MAE там, где стоимость ошибки асимметрична. Нужен квантиль, а не среднее.
8. Игнорирование смены режима. COVID, ребрендинг, смена ассортимента. Обучение на данных до разрыва без индикатора режима — гарантированное смещение.
9. Экстраполяция трендов деревьями. Разобрано выше: прогнозируйте отношения или детрендируйте.
10. Один фолд бэктеста. Дисперсия метрики на одном отсечении огромна; решение принимается по шуму.
11. Дырки во времени, интерпретированные как нули. Отсутствие строки за день ≠ ноль продаж: возможно, магазин был закрыт или ETL упал. Реиндексация по полному календарю с явным различением «ноль» и «нет данных» обязательна.
12. Праздники как бинарный флаг. 31 декабря и 25 декабря — разные эффекты; важны также дни до и после.

Жизненный цикл прогнозной системы в проде

Что из этого действительно важно на практике:

Fallback обязателен. Прогнозная система с внешними зависимостями (погода, партнёрские данные) рано или поздно останется без входа. Автоматический откат на сезонный наивный прогноз с пометкой в метаданных лучше, чем отсутствие числа или прогноз на устаревших признаках.

Ежедневный мониторинг — по горизонтам отдельно. Ошибка на дне 1 и дне 28 живёт своей жизнью; агрегированная метрика скроет деградацию короткого горизонта. Отслеживать надо MASE по горизонту, эмпирическое покрытие интервалов и смещение (mean error) — систематический перекос вверх или вниз опаснее шума, потому что напрямую превращается в лишний склад или дефицит.

Регламент переобучения. Три опции: по расписанию (просто, предсказуемо), по триггеру деградации (экономно, но нужен надёжный детектор), гибрид. Практика: еженедельное переобучение по расписанию плюс схема «чемпион–претендент», где новая модель выходит в прод только если побила текущую на свежем out-of-time окне.

Реестр признаков и точки во времени. Отдельный слой (feature store) с временными метками доступности каждого признака — единственный способ гарантировать, что обучение и инференс видят одни и те же данные. Детали инфраструктуры — в статье про MLOps.

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

Полезные библиотеки

Мини-итог

• Временные ряды нарушают i.i.d.: любое разбиение, кросс-валидация и построение признаков должны уважать стрелку времени.
• Классика (ARIMA, ETS) моделирует процесс порождения ряда: сильна на одиночных рядах с короткой историей, даёт статистически обоснованные интервалы, плохо масштабируется на длинные сезонности и панели.
• Признаковый подход сводит прогноз к регрессии; глобальный градиентный бустинг на лагах, скользящих статистиках и календаре — рабочая лошадка индустрии и победитель M5.
• Главный источник фиктивно хороших результатов — утечка из будущего: несдвинутые скользящие статистики, признаки, недоступные в момент прогноза, ретроспективно исправленные данные.
• Метрика должна отражать стоимость решения: MASE для сравнения поперёк рядов, pinball loss при асимметричной стоимости ошибки, MAPE — почти никогда.
• Бэктест на нескольких точках отсечения с gap, обязательное сравнение с seasonal_naive, мониторинг по горизонтам и fallback в проде.

Источники

• Rob J. Hyndman, George Athanasopoulos. Forecasting: Principles and Practice, 3rd ed. — бесплатный онлайн-учебник, лучший вход в тему.
• George Box, Gwilym Jenkins et al. Time Series Analysis: Forecasting and Control — первоисточник методологии ARIMA.
• Makridakis, Spiliotis, Assimakopoulos. The M4 Competition: 100 000 time series and 61 forecasting methods, 2020.
• Makridakis et al. The M5 competition: Background, organization, and implementation, 2022.
• Hyndman, Koehler. Another look at measures of forecast accuracy, 2006 — обоснование MASE.
• Cleveland et al. STL: A Seasonal-Trend Decomposition Procedure Based on Loess, 1990.
• Wickramasuriya, Athanasopoulos, Hyndman. Optimal forecast reconciliation (MinT), 2019.
• Lim et al. Temporal Fusion Transformers for Interpretable Multi-horizon Time Series Forecasting, 2019.
• Oreshkin et al. N-BEATS: Neural basis expansion analysis, 2019.
• Xu, Xie. Conformal prediction interval for dynamic time-series (EnbPI), 2021.

Что дальше

Во всех задачах этого трека модель делала предсказание и на этом её работа заканчивалась: прогноз не менял мир, в котором он был сделан. В следующей статье — Обучение с подкреплением — мы перейдём к постановке, где агент действует, среда отвечает, и каждое действие меняет распределение будущих данных. Марковские процессы принятия решений, Q-обучение, policy gradient и то, почему исследование против эксплуатации оказывается центральным компромиссом.