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

Статья написана по мотивам доклада Синимекс DATA Meetup, так что те, кому больше нравится видео-формат, могут посмотреть запись:

Запись доклада с митапа

А мы двинемся дальше в текстовом формате :)

Модельная задача

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

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

Сэмпл данных отгрузок может выглядеть примерно так

record_id

brik_id

product_id

ds

units

amount

1

1

1

2019-01-01

83

6972

2

1

1

2019-02-01

27

2268

3

1

1

2019-03-01

23

1932

4

1

1

2019-04-01

18

1512

5

1

1

2019-05-01

11

924

6

1

1

2019-06-01

5

420

7

1

1

2019-07-01

4

336

8

1

1

2019-08-01

7

588

9

1

1

2019-09-01

29

2436

10

1

1

2019-10-01

34

2856

Пример истории продаж и бейзлайн-прогноза одного продукта в территориальной единице
Пример истории продаж и бейзлайн-прогноза одного продукта в территориальной единице

Иерархии

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

Чаще всего используемая у клиента территориальная иерархия не сильно отличается от административного деления (федеральные округа (брик-1) > субъекты федерации (брик-2) > города > районы, деревни, ...), однако в нее могут вноситься важные для бизнеса изменения (например, Москву или Санкт-Петербург могут вытаскивать на уровень брик-1 наравне с федеральными округами из-за важности для бизнеса и т.п.).

Пример построения территориальной иерархии с 4-мя уровнями бриков
Пример построения территориальной иерархии с 4-мя уровнями бриков

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

Пример продуктовой иерархии: SKU (форма выпуска) > Бренд > Total. Названия препаратов и иллюстрации сгенерированы ChatGPT
Пример продуктовой иерархии: SKU (форма выпуска) > Бренд > Total. Названия препаратов и иллюстрации сгенерированы ChatGPT

Зачем об этом думать? DS-специалисты решают не просто интересную математическую задачу прогнозирования, а закрывают какую-то потребность бизнеса. У прогнозов есть конкретный потребитель, и для того чтобы понимать, что ему нужно, необходимо осознавать, на каком уровне иерархии он работает. Например, если потребителем наших предсказаний является генеральный директор, то самые низкие уровни предсказания по территориальной иерархии (продажи в небольших городах или районах городов) его в принципе не волнуют - его интересует, как его бизнес в следующем месяце будет вести себя по всей линейке препаратов и по всей территории суммарно. Регионального менеджера, который отвечает, например, за развитие бизнеса в республике Бурятия, могут интересовать точные прогнозы именно на этом уровне, а уровни выше, скорее всего, волновать не будут.

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

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

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

Реконсиляция

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

Рассмотрим игрушечную иерархию временных рядов в виде дерева из восьми узлов и трех уровней:

Пример из учебника https://otexts.com/fpp3/hts.html
Пример из учебника https://otexts.com/fpp3/hts.html

Отношения между временными рядами в подобной иерархии для любого отсчета t можно описать в виде системы линейных уравнений:

y_t = y_{A,t} + y_{B,t} \\ y_{A,t} = y_{AA,t} + y_{AB,t} + y_{AC,t} \\ y_{B,t} = y_{BA,t} + y_{BB,t} \\ y_t = y_{AA,t} + y_{AB,t} + y_{AC,t} + y_{BA,t} + y_{BB,t}

Запишем матричный вид полученной СЛАУ:

\textbf{y}_t = \textbf{S}\textbf{b}_t
Вектор-столбец y состоит из отсчетов всех временных рядов в иерархии (узлов дерева), вектор-столбец b - из отсчетов самых базовых временных рядов (листья дерева). Нижний блок матрицы является единичной матрицей, т.к. выражает тривиальные уравнения вида Ybb,t = Ybb,t
Вектор-столбец y состоит из отсчетов всех временных рядов в иерархии (узлов дерева), вектор-столбец b - из отсчетов самых базовых временных рядов (листья дерева). Нижний блок матрицы является единичной матрицей, т.к. выражает тривиальные уравнения вида Ybb,t = Ybb,t

Матрица S размера (n × m) (количество узлов × количество листьев) называется матрицей суммирования (summing matrix) и отвечает на вопрос «как получить отсчеты временных рядов по всей иерархии, если мы знаем отсчеты временных рядов на самом низком уровне иерархии?» Вид матрицы суммирования зависит только от дерева иерархии временных рядов и полностью им определяется. Построив матрицу суммирования, для получения отсчетов рядов по всем узлам иерархии нам достаточно иметь отсчеты для рядов-листьев.

Теперь допустим, что мы построили прогнозные модели в каждом узле дерева иерархии и в каждом узле имеем базовые несогласованные прогнозы. Для получения согласованных прогнозов мы теперь можем "спроецировать" информацию из всех узлов в листовые узлы (т.е. уточнить прогнозы в листьях за счет прогнозов в узлах), а после этого собрать уже согласованные прогнозы вверх по иерархии с помощью матрицы суммирования. Для того чтобы формализовать "проецирование", введем матрицу G размера (m × n) (количество листьев × количество узлов), которую назовем матрицей отображения (mapping matrix):

\tilde{\textbf{b}}_h = \textbf{G} \hat{\textbf{y}}_h
Здесь h - горизонт прогнозирования, y с крышкой - базовые прогнозы, y с тильдой - новые прогнозы в листьях
Здесь h - горизонт прогнозирования, y с крышкой - базовые прогнозы, y с тильдой - новые прогнозы в листьях

Конкретный вид матрицы G определяет алгоритм реконсиляции, а общий вид произвольного алгоритма линейной реконсиляции можно отобразить в таком виде:

\tilde{\textbf{y}}_h = \textbf{S} \textbf{G} \hat{\textbf{y}}_h

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

Top-Down

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

Подход Top-Down предусматривает строительство одной прогнозной модели на самом верхнем уровне какой-то выбранной нами иерархии (например, total по всей РФ), после чего этот прогноз через коэффициенты (подобранные экспертно или эвристически) распределяется между нижними уровнями.

Сам подход можно выразить простейшими формулами:

\tilde{y}_{i, h} = \alpha_i \hat{y}_{total, h} \\ \sum_i \alpha_i = 1

А в терминах нашего формализма и игрушечного примера он будет выглядеть так:

Матрица S не меняется. Единственный ненулевой столбец в матрице G - первый, т.к. для реконсиляции мы используем только прогноз самого верхнего уровня и раскладываем его в листья с помощью коэффициентов альфа.
Матрица S не меняется. Единственный ненулевой столбец в матрице G - первый, т.к. для реконсиляции мы используем только прогноз самого верхнего уровня и раскладываем его в листья с помощью коэффициентов альфа.

Какие у такого подхода есть плюсы?

  • Использование всего одной прогнозной модели на самом высоком уровне иерархии позволяет выделить много ресурсов на ее тюнинг и добиться хорошего качества.

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

Минусы логичным образом вытекают из плюсов:

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

  • По мере «распиливания» прогнозов по иерархии, дисперсия прогнозов на каждом отдельном уровне имеет тенденцию к увеличению.

  • Проблема определения оптимальных коэффициентов: есть несколько методов, каждый со своими плюсами и минусами. Например, можно использовать исторические данные по доле каждого брика в общих продажах (метод Historical Average), однако в таком случае эти коэффициенты будут всегда запаздывать, если на рынке есть динамика. Можно также построить какие-либо дополнительные прогнозные модели, которые будут прогнозировать эту долю (метод Forecast Proportions), что уже сложнее и требует дополнительных моделей.

Bottom-Up

Второй простой подход - полная альтернатива Top-Down - заключается в построении моделей только в листьях иерархии. Матрица G в таком случае будет совсем тривиальной, т.к. никакого "проецирования" информации и не происходит:

Матрица G просто "копирует" прогнозы в листьях
Матрица G просто "копирует" прогнозы в листьях

Плюсы такого подхода:

  • Поскольку модели строятся на самом низком уровне, мы можем "нативным" способом добавить в них какие-нибудь локальные признаки.

  • Сам подход в некотором смысле проще, чем Top-Down - не требуется дополнительно подбирать коэффициенты, нужно просто построить достаточно качественные модели прогнозирования в листьях.

Минусы:

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

  • Вместо одной модели у нас появляется m моделей, и заниматься уникальным тюнингом каждой отдельной модели становится накладно по ресурсам. Например, можно обучить одну модель на всю Россию, но уже сложнее обучить 86 моделей на каждый субъект Федерации. Мы можем прийти к десяткам тысяч моделей, и в каждую такую модель мы уже не можем вложить много времени, что ограничивает возможность улучшать их метрики.

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

Нужно учитывать, что метрики такого подхода при агрегации далее наверх могут вести к более низкому качеству прогнозов, чем Top-Down, и наоборот. Необходимо играть с балансом, определяя, что важнее в конкретной задаче - точные предсказания агрегированных величин или отдельных листьев в иерархии.

Также стоит упомянуть про подход Middle-Out - в нем мы строим прогнозные модели на каком-то промежуточном уровне иерархии, а потом от него двигаемся вверх как Bottom-Up, а вниз - как Top-Down. Матрица G для такого подхода будет смесью матриц GBU и GTD.

Multi-level reconсilation

Рассмотренные ранее методы называются single-level реконсиляцией, так как используют прогнозы только с какого-то одного уровня иерархии. Однако в пределе у нас может быть столько прогнозных моделей, сколько узлов в дереве иерархии - в таком случае нам потребуется multi-level реконсиляция.

Используя ранее введённый формализм, можно выражать разные алгоритмы реконсиляции в виде разных оптимизационных задач на поиск оптимальной матрицы G и различных дополнительных предпосылок или ограничений на такие задачи. Посмотрим на два подобных метода - MinTrace и ERM.

MinTrace (Minimum Trace)

MinTrace - это не так давно предложенный метод (статья опубликована в 2018г.). За подробным выводом отсылаю читателя к оригинальной статье "Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization" (Wickramasuriya, S. L., Athanasopoulos, G., & Hyndman, R. J., ссылка на DOI). Я же приведу только основные предпосылки и результаты из статьи.

Первым делом потребуем несмещенность вектора итоговых согласованных предсказаний, то есть

\mathbb{E}[\tilde{\textbf{y}}_h] = \textbf{y}_{T+h}

Из этого требования в статье выводится ограничение на матрицу G - она должна соответствовать уравнению

\textbf{S}\textbf{G}\textbf{S} = \textbf{S}

Интересное следствие данного ограничения - в статье показывается, что никакие Top-Down методы не удовлетворяют этому условию, а значит, Top-Down реконсиляция всегда дает смещенные прогнозы.

Будем искать оптимальную матрицу G в смысле минимизации среднеквадратичных ошибок согласованных прогнозов (т.е. минимизируем след матрицы ковариации ошибок, откуда и берется название метода). Матрицу ковариации согласованных прогнозов можно выразить через матрицы S, G и матрицу ковариации базовых прогнозов:

\textbf{G}_{opt} = \arg\min  \ Tr(\textbf{V}_h) \\ \textbf{V}_h = Var[\textbf{y}_{t+h} - \tilde{\textbf{y}}_h] = \textbf{S} \textbf{G} \textbf{W}_h \textbf{G}^T \textbf{S}^T \\ \textbf{W}_h = Var[\textbf{y}_{t+h} - \hat{\textbf{y}}_h]

Отсюда решение для G можно получить в аналитической форме (closed-form solution):

\textbf{G}_{opt} = (\textbf{S}^T \textbf{W}_h^{-1}\textbf{S})^{-1} \textbf{S}^T \textbf{W}_h^{-1}

Матрицу ковариации ошибок базовых прогнозов Wh нужно отдельно оценить, для этого в статье приводится несколько подходов к ее аппроксимации в разных случаях. Из дополнительного - у метода есть модификация Non-Negative MinTrace, которая может быть полезна в случае неотрицательного таргета.

ERM (Empirical Risk Minimization)

В математическом плане этот метод проще предыдущего - в нем не ставится никаких ограничений на матрицу G, а просто предлагается искать ее через регрессию. Метод описан в статье "Regularized Regression for Hierarchical Forecasting Without Unbiasedness Conditions" (ссылка на DOI).

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

\textbf{G}_{opt} = \arg\min ||Y - \textbf{S}\textbf{G}\hat{Y}||_2^2

где Y - это матрица, составленная из таргета на валидационной выборке, Y с крышкой - матрица из базовых прогнозов на этой выборке.

В случае обычной задачи регрессии решение также можно получить в аналитическом виде:

\textbf{G}_{opt} = (\textbf{S}^T\textbf{S})^{-1}Y^T\hat{Y}(\hat{Y}\hat{Y})^{-1}

Как и в обычной регрессии, можно добавить в задачу регуляризацию, например, L1 относительно Bottom-Up подхода:

\textbf{G}_{opt} = \arg\min ||Y - \textbf{S}\textbf{G}\hat{Y}||_2^2 + \lambda || \textbf{G} - \textbf{G}_{BU} ||_1

что может быть полезно, если возникает какая-то численная нестабильность при обращении матриц. Решения в аналитическом виде в таком случае уже не будет, но можно решить задачу, например, градиентным спуском или другими численными итеративными методами.

Плюсы и минусы

Каких плюсов можно ожидать от multi-level реконсиляции?

  • Очевидный плюс - можно получить согласованный прогноз, который в среднем по всей иерархии будет лучше, чем Top-Down, Bottom-Up или их смеси, а также каких-либо эвристических подходов, так как мы используем и агрегируем некоторым оптимальным образом информацию от прогнозных моделей сразу со многих уровней иерархий.

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

Минусы:

  • Нужно строить еще больше базовых прогнозных моделей в узлах иерархии - в пределе n узлов всего дерева, которое может быть достаточно большим. Это приводит к тому, что мы перестаем тюнить какие-то отдельные модели и выходим в обобщенные пайплайны, которые строят, например, 10 000 моделей по всем временным рядам, считают усредненные метрики и т.д. Это более сложная работа, и всегда есть риск что-то упустить.

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

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

Имплементации

Последний немаловажный вопрос - где можно найти готовые реализации упомянутых алгоритмов? Во многих фреймворках для прогнозирования временных рядов (Darts, ETNA, Greykite, Merlion) есть реализации разных методов реконсиляции (где-то только TopDown/BottomUp, где-то есть более продвинутые типа MinTrace). Я хотел бы обратить внимание на библиотеку HierarchicalForecast от Nixtla. Во-первых, в ней реализован самый широкий набор алгоритмов по сравнению с другими известными мне фреймворками, а во-вторых, ей удобно пользоваться именно как библиотекой - получить базовые предсказания с помощью любых других библиотек или фреймворков, а потом воспользоваться только алгоритмом реконсиляции.

На этом у меня всё, спасибо за внимание и точных вам прогнозов! :)

Комментарии (0)