Калман, Матлаб, и State Space Models

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

Во-первых, здесь явно присутствуют тренд и цикл с порядком около 300 дней. Теперь закроем график. Сходим покурим. Придем, откроем снова и посмотрим уже на сами цифры. Дата отгрузки указана в формате 01.09.2009, поэтому открываем в Экселе и если переведем данные в формат [$-F800]dddd, mmmm dd, yyyy], так чтобы показывался день недели, то заметим, что обычно в субботу значения по отгрузке намного меньше чем в остальные дни недели. Для примера две недели приведены на Графике 1. В общем, грузчик дядя Вася в субботу уходит домой пораньше смотреть футбол, а нам из-за этого придется дополнительно учитывать присутствие микроцикла с сезонностью в 7 дней. Кстати, мы не будем переводить данные в недельный интервал как в предыдущей статье, а продолжим работать с дневными данными.

Здесь я попытаюсь привести минимум теории, детальнее почему и откуда берутся все формулы можно прочитать в книге или я отвечу в комментариях. Допустим, имеется некий временной ряд , в нашем случае — данные по отгрузке товара. Как правильно указал kuznetsovin в предыдущей статье, данные явно являются нестационарным рядом с порядком интегрирования равном 1, и процедура ARIMA требовала бы дифференцирования рядов. Но так как мы этого делать не хотим, то следуя Harvey [1993], допустим что данный ряд может быть представлен как модель с ненаблюдаемыми компонентами (Unobservable Сomponent model):

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

Тренд можно разнообразить, построив его в виде модели локального линейного тренда:

где (2) — собственно тренд и (3) — наклон тренда, каждый со своими ошибками. Такая модель дает множество вариантов моделирования тренда, от random walk with drift до integrated random walk model. Многие эконометристы часто убирают ошибку в (2) чтобы получить чоткий гладкий тренд, а все случайные ошибки присвоить наклону тренда.

Стохастический цикл моделируется в виде суммы тригонометрических функций и их производных:

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

Ну и напоследок, недельный микроцикл с периодом s=7 также построим на основе суммы тригонометрических функций:

Теперь осталось все вышеприведенные формулы организовать в так называемый структурный формат подходящий для фильтра Калмана: Уравнение измерения (measurement equation) описывает данные которые мы наблюдаем:

и уравнение перехода (transition equation) — описывает динамику переменных, которые составляют , но мы их не видим (так называемые unobserved или латентные переменные):

В нашем случае, имеется 10 латентных переменных которые были определены в уравнениях (2)-(6):

Теперь переведем все вышеприведенные формулы в матричный формат. Матрица перехода в (7):

Неслабая такая матрица динамики латентных переменных в (8):

Вектор всех ошибок всех состояний (2)-(6):

оторый встраивается в уравнения динамики (8) с помощью матрицы .

И, наконец, матрица вариаций всех ошибок и пертурбаций в (8):

Ок, теперь мы готовы курить Калман. Об алгоритме фильтра уже неоднократно писали, разве что у нас немного изменены имена и фамилии обозначения переменных. Поэтому особо останавливаться на теории не будем, кратенько только так, минут на сорок. Итак, есть видимая переменная , и набор невидимых , для которых мы придумали динамику и структуру в (8)-(14), и которые мы и стараемся просчитать в попугаях с помощью фильтра Калмана. Так как латентные состояния невидимы, модель может быть не особо верна, да и обязательно пристутствуют разные ошибки измерения, то саму мы вряд ли встретим, а только ее математическое ожидание в момент времени t на основе наблюдаемых данных 1. t-1: , которое обладает дисперсией .

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

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

В дополнение к просчитыванию вектора состояний, нам еще интересно найти отдельные параметры модели, такие, например, как вариация ошибок, частота цикла, и параметр затухания цикла. Обзовем вектор искомых параметров как . Тогда если и распределены по Гауссу, то логарифмическая функция правдоподобия (Log-likelihood function) параметров для наших данных будет:

Так как и то результат из фильтра Калмана может быть использован чтобы вычислить функцию правдоподобия на основе ошибок фильтра:

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

Теперь пару слов о еще одной wunderwaffe — Калмановском Сглаживании (Kalman smoother), [Durbin, J., and Koopman, 2003], которая вроде на Хабре еще не упоминалось. В общем, идея похожа, только фильтр Калмана просчитывает каждое следуещее значение в момент времени t на основе предыдущих данных 1..t-1: . А Kalman Smoother немного читит, и, предполагая что у нас уже есть все данные, дает возможность уточнить смотря на весь временной ряд , то-есть, говоря простыми словами, он вычисляет . Это хорошо работает когда у нас уже есть все наблюдения, и интересуют не будущие оценки, а более лучше одеваться точные значения латентных переменных из которых складывается динамика ряда.

Сглаживание Калмана представляет собой обратную рекурсию:

где вектор сглаженных значений состояний будет иметь наименьшую дисперсию . Рекурсия сглаживания начинается в момент времени t=N, и запускается задом наперед задавая нулевые значения вектору и его дисперсии . Значения ошибки прогноза , ее дисперсия и Kalman gain берутся из фильтра, который прогоняется в первый заход.

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

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