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

в 4:59, , рубрики: Без рубрики

Недавно kuznetsovin опубликовал пост об использовании Питона для анализа временных рядов в экономике. В качестве модели была выбрана «рабочая лошадка» эконометрики — ARIMA, пожалуй, одна из наиболее распространенных моделей для временных данных. В то же время, главный недостаток АRIMA-подобных моделей в том, что они не приспособлены для работы с нестационарными рядами. Например, если в данных присутствует тренд или сезонность, то математическое ожидание будет иметь разное значение в разных участках серии — Калман, Матлаб, и State Space Models, что не есть хорошо. Для избежания этого, АRIMA предполагает работать не с исходными данными, а с их разностью (так называемое дифференцирование — от «taking a difference»). Все бы хорошо, но тут возникают две проблемы — (а) мы возможно теряем значимую информацию беря разницу ряда, и (б) упускается возможность разложить ряд данных на составляющие компоненты — тренд, цикл, и т.п. Поэтому, в данной статье я хотел бы привести альтернативный метод анализа — State Space Modeling (SSM), в русском переводе — Модель Пространства Состояний.

Примечание

Для многих англоязычных названий в этой сфере перевод либо отсутствует, либо разнится у разных переводчиков, либо же коряв до полной невозможности использовать его в приличном обществе. Поэтому, многие имена и названия будут приведены в английском варианте. Для интересующихся — одна из лучших книг по моделированию SSM. Также оказалось, что в русскоязычном переводе хороших книг по теме практически нет. Как вариант — работа Александра Цыплакова, которая хоть и опубликована как отдельная статья, но является фактически копией книги в предельно сокращенном варианте.

Итак, приступим.

1. Начальный анализ данных

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

График 1

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

2. State Space Modes

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

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

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

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

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

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

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

Неслабая такая матрица динамики латентных переменных в (8):
Калман, Матлаб, и State Space Models

Вектор всех ошибок всех состояний (2)-(6):
Калман, Матлаб, и State Space Models
оторый встраивается в уравнения динамики (8) с помощью матрицы Калман, Матлаб, и State Space Models.

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

3. Фильтр и Сглаживатель Калмана

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

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

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

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

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

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

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

4. Регрессия

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

Все работает все примерно так:

  1. Главная программа otgr_ssm.m загружает данные, подготавливает структуру ssmopt, в которой записаны много ценных замечаний и важных переменных которые будут передаваться в разные места кода. Последние 70 значений по отгрузке (10 недель) отложим в сторону для сравнения с прогнозом, который обязательно построим в конце.
  2. Данные выгружаются в подпрограмму максимизации Log-likelihood function. Последняя процедура часто используется в эконометрике, поэтому была быстро на коленке написана в общем виде отдельной функцией mle_my.m чтобы не замусоривать основной код.
  3. так как искомые параметры включают дисперсии состояний, они должны быть строго положительны, что не всегда получается в ходе числовой оптимизации. Поэтому, сменим пешки на рюмашки и все входные значения дисперсий будут трансформированы как Калман, Матлаб, и State Space Models, а значение коэффициента затухания нормировано в диапазоне Калман, Матлаб, и State Space Models как Калман, Матлаб, и State Space Models. Ну и на выходе, сооствественно трансформированы обратно, чтобы получить правильные значения — Калман, Матлаб, и State Space Models для вариаций и Калман, Матлаб, и State Space Models для коэффициента затухания. Значение ssmopt.trans указывает, надо ли трансформировать (1) или нет (0), и в какую сторону — 'in', или 'out'. Все это происходит с помощью отдельной функции transform.m
  4. mle_my.m вызывает обьектную функцию f_obj.m, которая прогоняет фильтр Калмана (15)-(22) с помощью kfmy2.m, вычисляет значение Log-likelihood function (22), и, чтобы два раза не вставать, просчитывает Сглаживатель Калмана (23)-(27) используя ksmy2.m. Это все отправляется назад в mle_my.m, которая проверяет, достигли мы максимума функции (22) или нет. Если да, то все на выход, в пункт 4. Если нет, то повторяем шаги 2-3.
    Самый первый запуск фильтра начинается с предположения что Калман, Матлаб, и State Space Models. Тут можно поиграться, так как переменных много, и мы ведем поиск глобального максимума шестимерной функции, поэтому возможны многие локальные экстремумы. По-хорошему, можно было бы построить графики возможных значений лог-функции и посмотреть на возможные экстремумы.
  5. Выводим результаты — найденные параметры и графики. Заодно посчитаем предсказание для 70 дней и сравним его с реальными данными используя frcst.m

Собственно, код:

otgr_ssm.m

clear all;
clc;
close all;
format long;

%------------------- 1. Load and prepare data ------------------------------
load otgruzka.mat;
      
% Structure ssmopt contains several important records used throughout the code
ssmopt.frcst=70;                  % forecast length
yend=y(end-ssmopt.frcst+1:end);   % saved last observation for the forecast comparison
y=y(1:end-ssmopt.frcst);
ssmopt.N=length(y);              
ssmopt.trans=1;    			% transform estimated parameters to preserve positiveness of variances
ssmopt.sv=[5e+8;500;5e+8;5e+8;0.05;0.9];    % starting values 
ssmopt.mle='f_obj';			 				% name of the objective function for the maximization
ssmopt.sv=transform(ssmopt.sv,'in',ssmopt);
ssmopt.filter='kfmy2';					 % name of the function computing Kalman Filter
ssmopt.smooth='ksmy2';					 % name of the function computing Kalman Smoother

%------------------- 2. Estimate the model ------------------------------
result=mle_my(y,ssmopt);				 % call Maximum Likelihood maximization function
b=transform(result.b,'out',ssmopt);	     % transform parameters back

% recompute data based on the correct non-transformed parameters
ssmopt.trans=0; 
ssmopt.sv=b;
[LH,KF_out,Ksm_out] = feval(ssmopt.mle,b,y, ssmopt);

% Recover filtered states series - trend, cycle, and seasonal
a_trend=KF_out.Afilt(1,:);
a_cycle=KF_out.Afilt(3,:);
a_seas=KF_out.Afilt(5,:)+KF_out.Afilt(7,:)+KF_out.Afilt(9,:);
y_filt=a_trend+a_cycle+a_seas;	% build the estimated filtered series Y

% Recover smoothed states series - trend, cycle, and seasonal
a_trendsm=Ksm_out.Asm(1,:);
a_cyclesm=Ksm_out.Asm(3,:);
a_seassm=Ksm_out.Asm(5,:)+Ksm_out.Asm(7,:)+Ksm_out.Asm(9,:);
y_smooth=a_trendsm+a_cyclesm+a_seassm;	% build the estimated smoothed series Y

result=mle_my(y,ssmopt);		% find correct Hessian for non-transformed parameters
%------------------- 3. Compute estimation statistics ------------------------------
%Find standard errors, and p-values
se=sqrt(abs(diag(inv(result.hessian)/ssmopt.N)));         % s.e.(b)
tstat=b./se;                             				  % t-statistics  
pval=2*(1-tcdf(abs(tstat),ssmopt.N-length(ssmopt.sv)));   % p-value

% Display output
fprintf('Estimated parameters and p-values:n');
out=[b pval]
period=2*pi/b(end-1)

% Compute R-squared
resid=y-y_filt;                                % estimation errors
SSE=resid*resid';                              % Sum of Squared Errors   
SST=(y-mean(y))*(y-mean(y))';                  % Sum of Squares Total
R2=1-SSE/SST                                   % R-squared'

%------------------- 4. Make Forecast ------------------------------
af0=KF_out.Afilt(:,end-ssmopt.frcst);
[yf,af]=frcst(b,y,ssmopt, af0);

%------------------- 5. Plot results ------------------------------
%p=ssmopt.N;
p=600;
t=[1:1:p];
figure(1)	
plot(t,y(1:p),'k',t,y_filt(1:p),'b',t,y_smooth(1:p),'r--')
title('Original, Filtered, and Smoothed data')
legend('y(t)','y filtered','y smoothed');

figure(2)	
plot(t,y(1:p),'k',t,a_trend(1:p),'b',t,a_trendsm(1:p),'r--')
title('Original data, Filtered and Smoothed trend')
legend('y(t)','Filtered trend','Smoothed trend');

figure(3)   
plot(t,a_cycle(1:p),'b',t,a_cyclesm(1:p),'r--')
title('Filtered and Smoothed cycle')
legend('Filtered cycle','Smoothed cycle');

figure(4)	% Filtered + smoothed seasonal
plot(t,a_seas(1:p),'b',t,a_seassm(1:p),'r--')
title('Filtered and Smoothed weekly seasonal')
legend('Filtered seasonal','Smoothed seasonal');

t=[1:1:ssmopt.frcst];
figure(5)
plot(t,yend,'b',t,yf,'r--')
title('Original data and Forecast')
legend('Original data','Forecast');

RMSE=sqrt(sum((yend - yf).^2)/ssmopt.frcst)  % Root Mean Squared Error

mle_my.m

function result_mle=mle_my(y,mleopt);
warning off;

%---------------- 1. Set-up minimization options ----------------
options=optimset('fminunc');  
options=optimset('LargeScale', 'off' , ...
                 'HessUpdate', 'bfgs' , ...
                 'LineSearchType', 'quadcubic' , ...
                 'GradObj' , 'off'  , ...  
                 'Display','off' , ...                  
                 'MaxIter' ,  1000  , ...
                 'TolX',  1e-12 , ...
                 'TolFun' , 1e-12, ... 
                 'DerivativeCheck' , 'off' , ...
                 'Diagnostics' , 'off' , ... 
                 'MaxFunEvals', 1000);

%---------------- 2. Run minimization ----------------
[b,fval,exitflag,output,grad,hessian]=fminunc(mleopt.mle,mleopt.sv,options,y,mleopt);
warning on;

result_mle.b=b;
result_mle.fval=fval;
result_mle.output=output;
result_mle.hessian=hessian;

f_obj.m

function [obj,KF_out,Ksm_out]=f_obj_tr(b,y,ssmopt);
   
%---------------- 1. Recover parameters ------------------------------------
b=transform(b,'out',ssmopt);
s2_irr=b(1);
s2_tr=b(2);
s2_cyc=b(3);
s2_seas=b(4);
freq=b(5);
rho=b(6);

%----------------  2. Build the model  ------------------------------------
ssmopt.ssmodel.states=10;
ssmopt.ssmodel.Z=[1 0 1 0 1 0 1 0 1 0];

T1 = [1 1 0 0; 0 1 0 0; 0 0 rho*cos(freq) rho*sin(freq); 0 0 -rho*sin(freq) rho*cos(freq)];
T2=[cos(2*pi/7) sin(2*pi/7) 0 0 0 0;... 
   -sin(2*pi/7) cos(2*pi/7) 0 0 0 0;...
	 0 0 cos(4*pi/7) sin(4*pi/7) 0 0;...
	 0 0 -sin(4*pi/7) cos(4*pi/7) 0 0;...
	 0 0 0 0 cos(6*pi/7) sin(6*pi/7);...
	 0 0 0 0 -sin(6*pi/7) cos(6*pi/7)];

ssmopt.ssmodel.T = [T1 zeros(rows(T1),cols(T2));zeros(rows(T2),cols(T1)) T2];
ssmopt.ssmodel.R=eye(10); ssmopt.ssmodel.R(1,1)=0;

H=s2_irr;

Q=zeros(10,10);
Q(2,2)=s2_tr; Q(3,3)=s2_cyc; Q(4,4)=s2_cyc; 
Q(5,5)=s2_seas; Q(6,6)=s2_seas; Q(7,7)=s2_seas; Q(8,8)=s2_seas;
Q(9,9)=s2_seas; Q(10,10)=s2_seas;

%----------------  3. Suggest starting conditions for the states ------------------------
a0=[y(1);0;0;0;0;0;0;0;0;0];
P0=eye(ssmopt.ssmodel.states)*1e+10;

%---------------- 4. Run Kalman filter ------------------------
KF_out = feval(ssmopt.filter,y, ssmopt, Q, H, a0, P0);
obj=KF_out.LH;

%---------------- 5. Run Kalman smoother ------------------------
if nargout>2
	ssmopt.ssmodel.num_etas=3;				 % number of the state variances (required for Kalman smoother)
	Ksm_out = feval(ssmopt.smooth,KF_out, ssmopt);
end

kfmy2.m

% Kalman filter 
%     y[t] = Z*alpha[t] + eps,  eps ~ N(0,H).
%     alpha[t] = T*alpha[t-1] + R*eta,  eta ~ N(0,Q).
%     v[t]=y[t]-E(y[t]) = y[t]-Z*a[t] 
%     F[t]=var(v[t])

function KF_out = kfmy_koop(y, ssmopt, Q, H, a, P);

N=ssmopt.N;
m=ssmopt.ssmodel.states;
%---------------- 1. Recover parameters and prepare matrices ----------------
T=ssmopt.ssmodel.T;
Z=ssmopt.ssmodel.Z;
R=ssmopt.ssmodel.R;
KF_out.Vmat=zeros(1,N);  KF_out.Fmat=zeros(1,N);
KF_out.Afilt=zeros(m,N); KF_out.Pfilt=zeros(m,m,N);
KF_out.Kmat=zeros(m,N);  KF_out.Lmat=zeros(m,m,N);
LHmat=zeros(1,N);

if ~isfield(ssmopt,'exactcheck');
		ssmopt.exactcheck=1;       % set exact filter initialization by default
end;

%---------------- 2. Set default starting values for a and P , if none provided  ----------------
Pinf=eye(m);   
if nargin < 6						  
	if ssmopt.exactcheck==1
		P=zeros(m,m);
	else
		P=eye(m)*1000000000;
	end	
end  
  
if nargin < 5
    a=[y(1); zeros(m-1,1)];
end
  
KF_out.Afilt(:,1)= a;
KF_out.Pfilt(:,:,1) = P;
d=0;

%---------------- 3. Exact Filtering ----------------
if ssmopt.exactcheck==1
	evals=10;      % number of time steps to evaluate until Pinf converges to zero
	KF_out.exact.F1=zeros(1,evals);
	KF_out.exact.F2=zeros(1,evals);
	KF_out.exact.L1=zeros(m,m,evals);
	KF_out.exact.Pinf=zeros(m,m,evals); KF_out.exact.Pinf(:,:,1)=Pinf;
	for i=1:evals
		if sum(sum(Pinf))<1e-20;
			d=i-1;                  % time point after which Pinf-->0, and after which we may start regular Kalman filter
			break;
		else
			if  sum(Pinf*Z')>0		% Pinf is not singular
				F1=inv(Z*Pinf*Z');	F2=-F1*(Z*P*Z'+H)*F1;
				K=T*Pinf*Z'*F1;		K1=T*(P*Z'*F1+Pinf*Z'*F2);
				L=T-K*Z;			L1=-K1*Z;
				P=T*Pinf*L1' + T*P*L' + R*Q*R';
				Pinf=T*Pinf*L';
			else
				F1=Z*P*Z'+H;		F2=F1;
				K=T*P*Z'/F1;
				L=T-K*Z;			L1=L;	
				P=T*P*L' + R*Q*R';
				Pinf=T*Pinf*T';
			end
			v=y(i) - Z*a;
			a=T*a+K*v;
		
			%save filtered estimates
			KF_out.Afilt(:,i+1)=a;  KF_out.Pfilt(:,:,i+1)=P;
			KF_out.Vmat(i)=v;       KF_out.Fmat(i)=F1;
			KF_out.Kmat(:,i)=K;     KF_out.Lmat(:,:,i)=L;      
			LHmat(i) = -0.5*(log(2*pi*F1) + v^2/F1);
			
			%save exact values for smoother
			KF_out.exact.F1(i)=F1;
			KF_out.exact.F2(i)=F2;
			KF_out.exact.L1(:,:,i)=L1;
			KF_out.exact.Pinf(:,:,i+1)=Pinf;
		end
	end
end

%---------------- 4.  Regular Filtering ----------------
for i=d+1:N
  v=y(i) - Z*a;
  f=Z*P*Z' + H;
  
  K=T*P*Z'/f;
  L=T-K*Z;
  
  a=T*a+K*v;
  P=T*P*L'+R*Q*R';

  if i<N
    KF_out.Afilt(:,i+1)=a;     KF_out.Pfilt(:,:,i+1)=P;
  end  
  
  KF_out.Vmat(i)=v;       KF_out.Fmat(i)=f;
  KF_out.Kmat(:,i)=K;     KF_out.Lmat(:,:,i)=L;      
  LHmat(i) = -0.5*(log(2*pi*f) + v^2/f);

end

%---------------- 5. Prepare output ----------------
KF_out.LH=-sum(LHmat);
KF_out.Q=Q;
KF_out.H=H;
KF_out.exact.d=d;
end

 

ksmy2.m

function [Ksm_out, Kdism_out] = ksmy2(KF_out, ssmopt); 

[m,N]=size(KF_out.Afilt);
meta=ssmopt.ssmodel.num_etas;

%---------------- 1. Recover filtered matrices ----------------
Fmat=KF_out.Fmat;
Vmat=KF_out.Vmat;
Pfilt=KF_out.Pfilt;
Afilt=KF_out.Afilt;
Q=KF_out.Q;
H=KF_out.H;

%---------------- 2. Recover Model structure ----------------
Z=ssmopt.ssmodel.Z;
T=ssmopt.ssmodel.T;
R=ssmopt.ssmodel.R;
Asm=zeros(m,N);
Psm=zeros(m,m,N);
rmat=zeros(m,N);
Nmat=zeros(m,m,N);
Eps=zeros(1,N);
Eta=zeros(meta,N);
Kmat=KF_out.Kmat;
Lmat=KF_out.Lmat;

if ~isfield(KF_out,'exact');		
	KF_out.exact.d=0;
end

d=KF_out.exact.d;

if KF_out.exact.d>0
	L1=KF_out.exact.L1;
	F1=KF_out.exact.F1;
	F2=KF_out.exact.F2;
	Pinf=KF_out.exact.Pinf;
end
%---------------- 3. Regular Smoothing for t=N..d+1 observations ----------------
for i=N:-1:d+1
	r=Z'/Fmat(i)*Vmat(i) + Lmat(:,:,i)'*rmat(:,i);
	N=Z'/Fmat(i)*Z + Lmat(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i);
	Asm(:,i)=Afilt(:,i) + Pfilt(:,:,i)*r;
	Psm(:,:,i)=Pfilt(:,:,i)-Pfilt(:,:,i)*N*Pfilt(:,:,i);
  
  if i>1
    rmat(:,i-1)=r;
    Nmat(:,:,i-1)=N;
  end  
  
  if nargout>1
    Eps(i)=H*(1/(Fmat(i))*Vmat(i)-Kmat(:,i)'*rmat(:,i));
    Eta(:,i)=Q*R'*rmat(:,i);
  end  
end
%---------------- 4. Exact Smoothing for t=d..1 observations ----------------
if KF_out.exact.d>0
	r1=zeros(m,1);	N1=zeros(m,m); N2=zeros(m,m);
	for i=d:-1:1
		if  sum(Pinf(:,:,i)*Z')>0		%cond(Pinf)<1e+12	% Pinf is not singular
			r1=Z'*F1(i)*Vmat(i) + Lmat(:,:,i)'*r1 + L1(:,:,i)'*rmat(:,i);
			N2=Z'*F2(i)*Z + Lmat(:,:,i)'*N2*Lmat(:,:,i) + Lmat(:,:,i)'*N1*L1(:,:,i) + L1(:,:,i)'*N1*Lmat(:,:,i) + L1(:,:,i)'*Nmat(:,:,i)*L1(:,:,i);
			N1=Z'*F1(i)*Z + Lmat(:,:,i)'*N1*Lmat(:,:,i) + L1(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i);

			r=Lmat(:,:,i)'*r1;
			N=Lmat(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i);
			
			if nargout>1
				Eps(i)=-H*Kmat(:,i)'*rmat(:,i);
				Eta(:,i)=Q*R'*rmat(:,i);
			end  
  
		else				% Pinf is singular
			r1=T'*rmat(:,i);
			N2=T'*N2*T;
			N1=T'*N1*Lmat(:,:,i);
			
			r=Z'/(Fmat(i))*Vmat(i) + Lmat(:,:,i)'*rmat(:,i);
			N=Z'/(Fmat(i))*Z + Lmat(:,:,i)'*Nmat(:,:,i)*Lmat(:,:,i);
			
			if nargout>1
				Eps(i)=H*(1/Fmat(i)*Vmat(i) - Kmat(:,i)'*rmat(:,i));
				Eta(:,i)=Q*R'*rmat(:,i);
			end  
		end		
		
		if i>1
			rmat(:,i-1)=r;
			Nmat(:,:,i-1)=N;
		end
		Asm(:,i)=Afilt(:,i) + Pfilt(:,:,i)*r + Pinf(:,:,i)*r1;
		Psm(:,:,i)=Pfilt(:,:,i)-Pfilt(:,:,i)*N*Pfilt(:,:,i) - (Pinf(:,:,i)*N1*Pfilt(:,:,i))' - Pinf(:,:,i)*N1*Pfilt(:,:,i) -  Pinf(:,:,i)*N2*Pinf(:,:,i);

	end
end
%---------------- 5. Prepare output ----------------
Ksm_out.Asm=Asm;
Ksm_out.Psm=Psm;
Ksm_out.Kmat=Kmat;
Ksm_out.Lmat=Lmat;
Ksm_out.Nmat=Nmat;
Ksm_out.rmat=rmat;

Kdism_out.Eps=Eps;
Kdism_out.Eta=Eta;

transform.m

function b=transform(b,howto,ssmopt);
k=length(b);

if strcmp(howto,'in') % in-transformation

  if ssmopt.trans==0  % no transformation
     b=b;
  end;

  if ssmopt.trans==1  % transformation to preserve the positiveness of variances
     b(1:k-1,:)=log(b(1:k-1,:));      
     b(k)=log(1/b(k)-1);
  end;

else   % out-transformation

  if ssmopt.trans==0  % no transformation
     b=b;
  end;

  if ssmopt.trans==1
    b(1:k-1,:)=exp(b(1:k-1,:));
    b(k)=1/(1+exp(b(k)));   
  end;
end

5. Результаты

(p-values найденных оценок указаны в скобках)

Дисперсия ошибки наблюдаемого ряда Калман, Матлаб, и State Space Models 1.77E+009 (0.00)
Дисперсия ошибки тренда Калман, Матлаб, и State Space Models 348.73 (0.00)
Дисперсия цикла Калман, Матлаб, и State Space Models 6.07E+008 (0.00)
Дисперсия ошибки сезонной компоненты Калман, Матлаб, и State Space Models 3.91E+006 (0.00)
Частота цикла Калман, Матлаб, и State Space Models 3.91E+006 (0.00)
Период цикла (в днях) 362.6 (0.00)
Коэффициент затухания циклаКалман, Матлаб, и State Space Models 0.891 (0.00)
R-квадрат регрессии 0.78

6. Графики

Для удобства восприятия, показаны графики только для первых 600 дней, для всего ряда спрятаны под спойлеры.
a. Исходные, отфильтрованные и сглаженные данные
Калман, Матлаб, и State Space Models

Весь ряд

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

b. Исходные данные, отфильтрованный тренд, сглаженный тренд
Как видно, фильтр Калмана пытается угадать тренд на основе предыдущих значений, и поэтому колеблется вместе с линией партии, но немного запаздывая, пытаясь угадать куда дальше направятся наши данные. Сглаживатель Калмана же «видит» весь ряд, и поэтому тренд выглядит намного ровнее и спокойнее:
Калман, Матлаб, и State Space Models

Весь ряд

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

c. Отфильтрованный и сглаженный цикл
Как видно из таблицы с результатами, средняя длина цикла составляет порядка 362 дня, или практически один год (кто бы удивлялся). Также видно как в самом начале фильтр начинает калиброваться и совершенно не попадает в данные, ведь мы задали начальные значения латентных переменных равными нулю и специально очень большой дисперсией порядка 1e+10. Но обычно достоаточно несколько (2-5) первых попыток чтобы попасть в ритм. Кстати, в этой работе использовалась точная инициализация фильтра (exact initialization), которая помогает отфильтрованным значениям быстрее сбежаться с данными.
Калман, Матлаб, и State Space Models

Весь ряд

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

d. Отфильтрованный и сглаженный недельный сезонный фактор
Постепенно растет количество отгрузок, увеличивается и дневная волатильность — лишь мгновение ты наверху и стремительно падаешь вниз.
Калман, Матлаб, и State Space Models

Весь ряд

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

6. Прогноз

На основе полученных параметров и используя последние значения отфильтрованных состояних, строим прогноз на 70 дней (10 недель) и сравниваем с существующими данными. В целом прогноз выглядит неплохо, вот что фильтр животворящий делает. Особенно радует угаданная волатильность по дням недели. Если присмотреться ко всей длине построенного прогноза и включить воображение, то видно еще и как прогноз прогибается под годовой цикла отгрузки. Единственный момент где прогноз не сработал — с 26 по 32 день. Но там явно случилось почти недельное падение отгрузки, так же как и резкий скачок сразу за ней, которые вряд-ли возможно было предвидеть, так как подобный случай встречался только единожды в самом начале серий, в общем, а что упал, так то от помутненья.
Калман, Матлаб, и State Space Models
RMSE прогноза — 1.112e+005, на случай если мы захотим сравнить модель.

Ну вот и все.

Примечание

Не следует воспринимать State Space Models как что-то стоящее осторонь в эконометрике. Наоборот, они представляют собой обобщенный вариант для многих более специфических моделей. Например, MA(1) процесс
Калман, Матлаб, и State Space Models
легко представить в форме SSM:
Калман, Матлаб, и State Space Models
Или же ARMA(2,1) процесс:
Калман, Матлаб, и State Space Models
Упакованный в формат SSM:
Калман, Матлаб, и State Space Models

Литература по теме

  • Durbin, J., and Koopman, S.J. «Time Series Analysis by State Space Methods». Oxford: Oxford University Press, 2001.
  • Durbin, J., and Koopman, S.J. «A simple and efficient simulation smoother for state space time series analysis» Biometrika vol. 89, issue 3, 2002.
  • Harvey, A.C., and Jaeger, A. “Detrending, stylised facts and the business cycle.” Journal of Applied Econometrics (8), 1993.
  • Harvey, A.C. «Forecasting, Structural Time Series Models and the Kalman Filter». Cambridge: Cambridge University Press, 1989.

Автор: werwooolf

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js