Фильтр Калмана, руководство полного идиота

Итак, что же такое "фильтр Калмана"?

Предположим у вас есть некая величина X изменяющаяся во времени.

У вас есть некий прибор, измеряющий величину X с некоторой погрешностью и выдающий результат измерения Y.

Вы хотите оценить, чему на самом деле равна величина X - "отфильтровать погрешность, шум".

Например, в случае акций Xk может быть истинной ценой акции в момент времени k (она определяется доходами компании, курсом валют, ожиданиями инвесторов итп).

У нас есть только один способ измерить эту величину - пойти на биржу и посмотреть котировки. Но котировки на бирже подвержены влиянию массы факторов, цены колеблются в широком диапазоне и измеренная величина Yk может отличаться от Xx. Алгоритм Калмана позволяет нам оценить истинное значение Xk на основе истории изменения измеряемой величины Yk.

При этом, естественно, делается ряд предположений о связи между Y и X и характере изменения X.

Предположение номер 1

Xk = Ak * Xk-1 + Bk * Uk + Wk-1

Т.е. величина X в момент времени k представляет собой линейную комбинацию величины X в прошлый момент времени, управляющего сигнала U и шума W.

Что такое управляющий сигнал? Вернемся к нашему примеру с акциями. Предположим, что ФРС США активно печатает деньги, например 600 млрд долларов, эти деньги льются широкой рекой на рынок и неизбежно двигают цены акций вверх. Это и есть управляющий сигнал. В реальных задачах часто управляющий сигнал U отсутствует.

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

Предположение номер 2

Yk = Hk * Xk + Vk

Т.е. наша измеряемая величина Y это линейная комбинация реальной величины X и погрешности измерения V. При этом, делается дополнительное предположение что погрешность измерения тоже имеет нормальное распределение со своей дисперсией dV.

H - матрица измерителя. Очень часто на практике H = 1.

Что мы только что сделали?

Мы построили модель, описывающую нашу систему, при этом, заметьте, пока никакого фильтра Калмана не упоминалось.

Что мы хотим сделать?

Мы хотим получить некое уравнение

Xk = F (Xk-1, k-2..., Yk, k-1..., dW, dV)

позволяющее нам оценивать истинное значение X исходя из исторических оценок и измеряемых величин.

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

 

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

В конце концов, получается выражение для оценки X:

Xk = Kk * Yk + (1 - Kk) * Xk-1

Где Kk - это Калмановский коэффициент усиления. Именно эту величину, в процессе работы и определяет алгоритм фильтра Калмана. При этом фильтр Калмана находит оптимальную величину коэффициента усиления при заданных дисперсиях dW и dV. Чем лучше вы оцените эти дисперсии, тем лучше будет работать фильтр Калмана.

А теперь, представьте, что X и Y это не одномерные величины, а вектор. A, B и H - матрицы. Погрешность и шум V и W тоже векторы. Оказывается, что при этом, все уравнения будут работать точно так же, и алгоритм фильтра Калмана легко обобщается на многомерные величины.

 

А теперь, пример на R

Фильтр Калмана можно найти в нескольких модулях R, для примера я буду использовать реализацию из модуля dlm.

В данном примере мы применим фильтр Калмана к значениям индекса S&P500.

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

Итак, пример:


library(dlm)
library(quantmod)

getSymbols("^GSPC", from="2010-09-01")
y <- as.ts(Cl(GSPC))
x <- dlmFilter(y, dlmModPoly(1, dV=0.1, dW=0.01))

plot(y, type='l')
lines(x$f[-1], col = "red")
#Добавим EMA для сравнения
lines(EMA(y), col = "green")

Получилась картинка:

Для сравнения, на график, зеленым цветом, было добавлено 10-дневное экспоненциальное скользящее среднее