function [lnL,status,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(Y,X,A,H,F,Q,R,KsiInit,initP,MaxIter,Tolerance,StartPeriod,AllowUnitRoot,c)
% KalmanFilterHt: Calculates the value of the log-likelihood function (based on Gaussian errors and initial entry
%                 for the state vector) for the observed data using the state-space representation.
%
% USAGE:
%
%       [lnL,status,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(Y,X,A,H,F,Q,R,KsiInit,initP,MaxIter,Tolerance,StartPeriod,AllowUnitRoot,c)
%
% The measurement equation is given by:
%
%                      y(t) = A'*x(t) + H(t)'*ksi(t) + w(t)
%
% while the state equation is:
%
%                      ksi(t) = F*ksi(t-1) + v(t)
%
% while w(t) and v(t) are white noise with zero mean and covariance matrices R and Q, respectively. The vector x(t) is
% predetermined. Furthermore, ksi(1) is uncorrelated with v(t) and w(t).
%
% REQUIRED INPUTS:  Y (nxT) matrix with observed data.
%
%                   X (kxT) matrix with exogenous of predetermined data.
%
%                   A (kxn) matrix with parameter values on the k X variables in the measurement equation.
%
%                   H (rxnxT) 3D matrix with parameter values on the r state variables (ksi) in the measurement equation.
%
%                   F (rxr) matrix with parameter values on lagged state variables in the state equation.
%
%                   Q (rxr) positive semidefinite matrix with covariances for the innovations in the state equation.
%
%                   R (nxn) positive semidefinite matrix with covariances for the innovations in the measurement equation.
%
%                   KsiInit (rx1) vector with initial values for the state vector.
%
%                   initP, variable indicating method of determining initial value for P(1|0), the MSE matrix for
%                          ksi(1|0). If initP=1 (default), the P(1|0) is given by the unconditional covariance matrix of ksi,
%                          calculated through the vec operator. Since this operation can be excessively slow for large systems
%                          (when r^2 is a large number) the Doubling Algorthm can be chosen to approximate the unconditional
%                          variance. This algorithm is called if initP=2. Finally, if initP=3, then P(1|0) = c*eye(r*r).
%                   NOYE: If F has eigenvalues outside the unit circle, then initP is always set to 3!
%
%                   MaxIter (integer), the maximum number of iterations to use with the doubling algorithm
%                          divided by 100. Must be multiplied by 100 when calling DoublingAlgorithmLyapunov.
%
%                   Tolerance (integer), determines the tolerance value for the doubling algorithm.
%
%                   StartPeriod, positive integer indicating first observation to use when calculating the log-likelihood value. The
%                       default value is 1.
%
%                   AllowUnitRoot (boolean) which is 1 if unit roots should be allowed for 0 otherwise.
%
%                   c, positive real number indicating the variance of P(1|0) when initP=3. The default value is 1000.
%
% REQUIRED OUTPUTS: lnL, scalar value of the log-likelihood function.
%
%                   status (boolean) indicating the covariance stationarity of the state variables. If all the
%                          eigenvalues of F are inside the unit circle, then status=0 and the log-likelihood function
%                          is calculated. If not, then status=1 and initP is set to 3.
%
% OPTIONAL OUTPUTS: lnLt (1x(T-StartPeriod+1)) vector with the time t values of the log-densities that sum to the value lnL.
%
%                   Yhat (nx(T-StartPeriod+1)) matrix with the predictions [y(StartPeriod|StartPeriod-1) ... y(T|T-1).
%
%                   MSEY (nxnx(T-StartPeriod+1)) 3 dimensional matrix with MSEY(:,:,t-StartPeriod+1) = H'*P(t|t-1)*H + R, the mean
%                        squared error matrix for y(t|t-1).
%
%                   Ksitt1, (rx(T-StartPeriod+1)) matrix with the predictions [ksi(StartPeriod|StartPeriod-1) ... ksi(T|T-1)].
%
%                   Ptt1 (rxrx(T-StartPeriod+1)) 3 dimensional matrix with Ptt1(:,:,t-StartPeriod+1) = P(t|t-1)
%
%
% DOCUMENTATION: Warne, A. (2006), "YADA Manual - Mathematical Details".
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: July 4, 2008.
%                        This version: December 18, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 18-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% initializing
%
lnL = 0;
status = 0;
%
% check the covariance stationarity of the state vector.
%
if max(abs(eig(F)))>=1;
   if AllowUnitRoot==0;
      status = 1;
   elseif (AllowUnitRoot==1)&(max(abs(eig(F)))>1)
      status = 1;
   end;
   initP = 3;
end;
%
% Determine dimensions. Checking if dimensions match is not undertaken!
%
[n,T] = size(Y);
if isempty(X)==0;
   k = size(X,1);
else;
   k = 0;
end;
r = size(H,1);
%
% Compute initial values
%
if initP==1;
   P_old = reshape(inv(eye(r*r)-kron(F,F))*vec(Q),r,r);
elseif initP==2;
   %
   % The maximum number of iterations is set to MaxIter*100, while the convergence value is 1.0e-8
   %
   [P_old,stat] = DoublingAlgorithmLyapunov(F,Q,MaxIter*100,10^(-(Tolerance+1)));
else;
   if c<=0;
      c = 1000;
   end;
   P_old = c*eye(r);
end;
%
% Initialize temporary matrices
%
if nargout>2;
   Yhat_temp = zeros(n,T);
   MSEY_temp = zeros(n,n,T);
   Ksitt1_temp = zeros(r,T);
   Ksitt1_temp(:,1) = KsiInit;
   Ptt1_temp = zeros(r,r,T);
end;
ksi_old = KsiInit;
lnLt_temp = zeros(1,T);
if k>0;
   Ax = A'*X;
else;
   Ax = zeros(1,T);
end;
%
% Now we can compute the Kalman filtered series
%
lnL_const = -n*log(2*pi);
for t=1:T;
   inv_msey = inv((H(:,:,t)'*P_old*H(:,:,t))+R);
   y_pred = Ax(:,t) + (H(:,:,t)'*ksi_old);
   y_error = Y(:,t)-y_pred;
   %
   % calculate the Kalman gain matrix
   %
   kalman_gain = F*P_old*H(:,:,t)*inv_msey;
   lnLt_temp(1,t) = lnL_const + (log(det(inv_msey))) - (y_error'*inv_msey*y_error);
   %
   % store values if asked to do so
   %
   if nargout>2;
      Yhat_temp(:,t) = y_pred;
      MSEY_temp(:,:,t) = (H(:,:,t)'*P_old*H(:,:,t))+R;
      Ksitt1_temp(:,t) = ksi_old;
      Ptt1_temp(:,:,t) = P_old;
   end;
   %
   % now we update ksi_old and P_old
   %
   ksi_old = (F*ksi_old)+(kalman_gain*y_error);
   F_dev = F-(kalman_gain*H(:,:,t)');
   P_old = (F_dev*P_old*F_dev') + (kalman_gain*R*kalman_gain') + Q;
end;
lnL = (1/2)*sum(lnLt_temp(1,StartPeriod:T));
%
% Fix output to match StartPeriod
%
if nargout==2;
   return;
elseif nargout==3;
   lnLt = (1/2)*lnLt_temp(1,StartPeriod:T);
elseif nargout==4;
   lnLt = (1/2)*lnLt_temp(1,StartPeriod:T);
   Yhat = Yhat_temp(:,StartPeriod:T);
elseif nargout==5;
   lnLt = (1/2)*lnLt_temp(1,StartPeriod:T);
   Yhat = Yhat_temp(:,StartPeriod:T);
   MSEY = MSEY_temp(:,:,StartPeriod:T);
elseif nargout==6;
   lnLt = (1/2)*lnLt_temp(1,StartPeriod:T);
   Yhat = Yhat_temp(:,StartPeriod:T);
   MSEY = MSEY_temp(:,:,StartPeriod:T);
   Ksitt1 = Ksitt1_temp(:,StartPeriod:T);
elseif nargout==7;
   lnLt = (1/2)*lnLt_temp(1,StartPeriod:T);
   Yhat = Yhat_temp(:,StartPeriod:T);
   MSEY = MSEY_temp(:,:,StartPeriod:T);
   Ksitt1 = Ksitt1_temp(:,StartPeriod:T);
   Ptt1 = Ptt1_temp(:,:,StartPeriod:T);
end;

%
% end of KalmanFilterHt.m
%
