function [L,varargout] = kalcvf(data, lead, a, F, b, H, var, varargin)
%KALCVF The Kalman filter
%
% State space model is defined as follows:
% z(t+1) = a+F*z(t)+eta(t) (state or transition equation)
% y(t) = b+H*z(t)+eps(t) (observation or measurement equation)
%
% [logl, >]=kalcvf(data, lead, a, F, b, H, var, )
% computes the one-step prediction and the filtered estimate, as well as their covariance matrices.
% The function uses forward recursions, and you can also use it to obtain k-step estimates.
%
% The inputs to the KALCVF function are as follows:
% data is a Ny x T matrix containing data (y(1), ... , y(T)). NaNs are treated as missing values.
% lead is the number of steps to forecast after the end of the data.
% a is an Nz x 1 vector for a time-invariant input vector in the transition equation.
% F is an Nz x Nz matrix for a time-invariant transition matrix in the transition equation.
% b is an Ny x 1 vector for a time-invariant input vector in the measurement equation.
% H is an Ny x Nz matrix for a time-invariant measurement matrix in the measurement equation.
% var is an (Ny+Nz) x (Ny+Nz) matrix for a time-invariant variance matrix for
% the error in the transition equation and the error in the measurement equation,
% that is, [eta(t)', eps(t)']'.
% z0 is an optional Nz x 1 initial state vector.
% vz0 is an optional Nz x Nz covariance matrix of an initial state vector.
%
% The KALCVF function returns the following output:
% logl is a value of the average log likelihood function of the SSM
% under assumption that observation noise eps(t) is normally distributed
% pred is an optional Nz x (T+lead) matrix containing one-step predicted state vectors.
% vpred is an optional Nz x Nz x (T+lead) matrix containing mean square errors of predicted state vectors.
% filt is an optional Nz x T matrix containing filtered state vectors.
% vfilt is an optional Nz x Nz x T matrix containing mean square errors of filtered state vectors.
%
% The initial state vector and its covariance matrix of the time invariant Kalman filters
% are computed under the stationarity condition:
% z0 = (I-F)\a
% vz0 = (I-kron(F,F))\V(:)
% where F and V are the time invariant transition matrix and the covariance matrix of transition equation noise,
% and vec(V) is an Nz^2 x 1 column vector that is constructed by the stacking Nz columns of matrix V.
% Note that all eigenvalues of the matrix F are inside the unit circle when the SSM is stationary.
% When the preceding formula cannot be applied, the initial state vector estimate is set to a
% and its covariance matrix is given by 1E6I. Optionally, you can specify initial values.
%
% This is a M-file for MATLAB.
% $Revision: 1.4 $ $Date: 09/30/2016 $
% Author: Iskander Karibzhanov (kais@bankofcanada.ca)
%==========================================================================
% Revision history:
%
% 03/19/2003 - algorithm and interface were adapted from SAS/IML KALCVF subroutine for use in MATLAB M file
% 09/17/2015 - added line P=(P+P')/2 to make sure P is numerically symmetric and does not explode over time.
% 09/30/2016 - NaNs in the data matrix are treated as missing values
%
%==========================================================================
T = size(data,2);
Nz = size(a,1);
Ny = size(b,1);
nin = nargin;
if nin~=7 && nin~=9
error('Seven or nine input arguments required.')
end
if nin==9
z = varargin{1};
P = varargin{2};
end
nout = nargout;
if nout~=1 && nout ~=3 && nout ~=5
error('One, three, or five output arguments required.')
end
% check input matrix dimensions
if size(data,1)~=Ny
error('data and b must have the same number of rows')
end
if size(a,2)~=1
error('a must be column vector')
end
if any(size(F)~=[Nz Nz])
error('F must be square')
end
if size(b,2)~=1
error('b must be column vector')
end
if any(size(H)~=[Ny Nz])
error('H must be Ny by Nz matrix')
end
if any(size(var)~=[(Ny+Nz) (Ny+Nz)])
error('var must be (Ny+Nz) by (Ny+Nz) matrix')
end
if nin==9 && any(size(z)~=[Nz 1])
error('z0 must be column vector of length Nz')
end
if nin==9 && any(size(P)~=[Nz Nz])
error('vz0 must be Nz by Nz matrix')
end
% V(t) and R(t) are variances of eta(t) and eps(t), respectively,
% and G(t) is a covariance of eta(t) and eps(t)
V = var(1:Nz,1:Nz);
R = var(Nz+1:end,Nz+1:end);
G = var(1:Nz,Nz+1:end);
if nin==7
e = eig(F);
if all(all(e*e'-eye(Nz)))
z = (eye(Nz)-F)\a;
P = reshape((eye(Nz^2)-kron(F,F))\V(:),Nz,Nz);
else
z = a;
P = eye(Nz)*1e6;
end
end
if nout>1
pred = zeros(Nz,T+lead);
vpred = zeros(Nz,Nz,T+lead);
end
if nout==5
filt = zeros(Nz,T);
vfilt = zeros(Nz,Nz,T);
end
L = 0;
if nout>1
pred(:,1) = z;
vpred(:,:,1) = P;
end
for t=1:T
i = ~isnan(data(:,t));
h = H(i,:); Pht = P*h';
e = data(i,t)-b(i,1)-h*z; % prediction error
D = h*Pht+R(i,i); % prediction error variance
if nout==5
filt(:,t) = z+Pht*(D\e);
vfilt(:,:,t) = P-Pht/D*Pht';
end
K = (F*Pht+G(:,i))/D;
z = a+F*z+K*e;
P = F*P*F'+V-K*D*K';
P = (P+P')/2;
if nout>1
pred(:,t+1) = z;
vpred(:,:,t+1) = P;
end
L = L-(log(det(D))+e'*D\e+sum(i)*log(2*pi))/2;
end
if lead>=1 && nout>1
for t=T+1:T+lead
z = F*z+a;
P = F*P*F'+V;
pred(:,t+1) = z;
vpred(:,:,t+1) = P;
end
end
if nout>1
varargout(1) = {pred};
varargout(2) = {vpred};
end
if nout==5
varargout(3) = {filt};
varargout(4) = {vfilt};
end