function [theta_upd,theta_pred,P_upd,P_pred]=FFBS_kfilter(Z,C,F,CSI,Gamma,Q,theta00,P00)

%  Kalman filter steps prepared for FFBS implementation
%  Innovations of the measurement equations are heteroskedastic
%  Coefficients of measurement equation are states to be filtered
%  Specification is:
%
%  Z_t     = c_t+ F_t * theta_t     + eps_t, eps_t iid multnormal(0,Csi_t)
%  theta_t =    Gamma * theta_{t-1} + v_t,     v_t iid multnormal(0,Q)
%
%  Z_t is nx1, c_t is nx1  , F_t  is nxm , theta_t is mx1 , Csi_t is nxn
%
%  sample line: [theta_upd,theta_pred,P_upd,P_pred]=...
%                       FFBS_kfilter(Z,C,F,CSI,Gamma,Q,theta00,P00)
%
%  INPUTS:   Z is nxT ,  F is nxmxT  , CSI is nxnxT, C is nxT
%            Gamma is mxm,    Q is mxm    (state equation)
%            theta00 is mx1 and P00 is mxm  (initial conditions)
%  OUTPUTS: theta_upd is mx(T+1),    P_upd is  mxmx(T+1)
%           [the first time instance (the +1) is the initial condition]
%           theta_pred is mxT,       P_pred is mxmxT
%  Author: Francesco Corsello, May/June 2014

[n,T]=size(Z);
[~,m,~]=size(F);

theta_upd=zeros(m,T+1);  P_upd=zeros(m,m,T+1);
theta_pred=zeros(m,T); P_pred=zeros(m,m,T);

theta_upd(:,1)=theta00; P_upd(:,:,1)=P00;

gammaidn_bool=isequal(Gamma,eye(m));

for t=1:T
%prediction step
if gammaidn_bool
    theta_pred(:,t)=theta_upd(:,t);
    P_pred(:,:,t)=P_upd(:,:,t)+Q;
else
    theta_pred(:,t)=Gamma*theta_upd(:,t);
    P_pred(:,:,t)=Gamma*P_upd(:,:,t)*Gamma'+Q;
end

Z_pred=C(:,t)+F(:,:,t)*theta_pred(:,t);
kkt=P_pred(:,:,t)*F(:,:,t)';
Kkt=F(:,:,t)*kkt+CSI(:,:,t);
Kt=kkt/Kkt;
%[size(kkt) size(Kkt) size(Kt)]

errZ=Z(:,t)-Z_pred;

%updating step
theta_upd(:,t+1)=theta_pred(:,t)+Kt*errZ;
P_upd(:,:,t+1)=(eye(m)-Kt*F(:,:,t))*P_pred(:,:,t);
end


end