function [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(Y,Yhat,Ksitt1,Ptt1,H,F,R)
% StateSmootherHt: Calculates smoothed and updated values of the state variables from the Kalman filter, as well as
%                  the MSE matrices for these predictors.
%
% USAGE:
%
%       [Ksitt,Ptt,KsitT,PtT,Ksit1t] = StateSmootherHt(Y,Yhat,Ksitt1,Ptt1,H,F,R)
%
% The measurement equation is given by:
%
%                      y(t) = A'*x(t) + H'*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).
%
% The updated and smooth predictions along with their MSE matrices are calculated using the formulas from
% Hamilton (1994, Chapter 13). They are also given in Warne (2006).
%
%
% REQUIRED INPUTS:  Y (nxT) matrix with observed data.
%
%                   Yhat (nxT) matrix with the forecasted data (1 step ahead forecasts).
%
%                   Ksitt1 (rxT) matrix with the forecasted state vector (1 step ahead).
%
%                   Ptt1 (rxrxT) 3 dimensional matrix with the MSEs for the forecasted state vector.
%
%                   H (rxnxT) 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.
%
%                   R (nxn) positive semidefinite matrix with covariances for the innovations in the measurement equation.
%
% REQUIRED OUTPUTS: Ksitt (rxT) matrix with the updated predictions of the state vector.
%
%                   Ptt (rxrxT) 3 dimensional matrix with the MSEs of the updated predictions of the state vector.
%
%                   KsitT (rxT) matrix with the smooth predictions of the state vector.
%
%                   PtT (rxrxT) 3 dimensional matrix with the MSEs of the smooth predictions of the state vector.
%
%                   Ksit1t (rxT-1) matrix with 1-step smooth predictions of ksi(t-1) conditional on time t data.
%
%
% 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.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[r,T] = size(Ksitt1);
n = size(R,1);
%
% Initializa the output
%
Ksitt = zeros(r,T);
Ptt = zeros(r,r,T);
KsitT = zeros(r,T);
PtT = zeros(r,r,T);
Ksit1t = zeros(r,T-1);
%
% we loop from period T and backwards to t=1
%
t = T;
while t>=1;
   inv_msey = inv((H(:,:,t)'*Ptt1(:,:,t)*H(:,:,t))+R);
   %
   % calculate the updated state prediction and its MSE
   %
   Ksitt(:,t) = Ksitt1(:,t)+(Ptt1(:,:,t)*H(:,:,t)*inv_msey*(Y(:,t)-Yhat(:,t)));
   Ptt(:,:,t) = Ptt1(:,:,t)-(Ptt1(:,:,t)*H(:,:,t)*inv_msey*H(:,:,t)'*Ptt1(:,:,t));
   if t<T;
      if rank(Ptt1(:,:,t+1))==length(Ptt1(:,:,t+1));
         m = eye(r);
      else;
         m = SchurDecomposition(Ptt1(:,:,t+1));
      end;
      %
      % Now we roll back the smoother
      %
%      J_t = m*m'*Ptt(:,:,t)*m*m'*F'*m*inv(m'*Ptt1(:,:,t+1)*m)*m';
      J_t = Ptt(:,:,t)*F'*m*inv(m'*Ptt1(:,:,t+1)*m)*m';
      KsitT(:,t) = Ksitt(:,t)+(J_t*(KsitT(:,t+1)-Ksitt1(:,t+1)));
      PtT(:,:,t) = Ptt(:,:,t)+(J_t*(PtT(:,:,t+1)-Ptt1(:,:,t+1))*J_t');
      Ksit1t(:,t) = Ksitt(:,t)+(J_t*(Ksitt(:,t+1)-Ksitt1(:,t+1)));
   else;
      %
      % Updates and smooth predictions and MSEs have the same values for t=T
      %
      KsitT(:,T) = Ksitt(:,T);
      PtT(:,:,T) = Ptt(:,:,T);
   end;
   t = t-1;
end;

%
% end of StateSmootherHt.m
%
