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

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

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;

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

% 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

 

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;

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

• 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.