function [thetas_sample,theta0_sample]=FFBS_sampling(theta_upd,theta_pred,P_upd,P_pred,Gamma)

%  Backward sampling steps prepared for FFBS implementation
%  see Carter and   Kohn (1994)
%
%  Start to sample from multnormal ( theta_(T|T), P(T|T) )
%  and then backward each theta(t) for t=T-1:-1:0 from
%  MultiNormal( theta(t|t+1), P(t|t+1) )  where
%  theta(t|t+1)= theta(t|t)+ P(t|t)*Gamma'*(P(t+1|t))^(-1) * (theta(t+1)-Gamma*theta(t|t))
%  P(t|t+1)   =  P(t|t)- P(t|t) Gamma' P(t+1|t)^(-1) Gamma P(t|t)
%
%  line: [theta_sample]=FFBS_sampling(theta_upd,theta_pred,P_upd,P_pred)
%
%  INPUTS: 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
%           Gamma is mxm   (if no Gamma is supplied, Gamma=identity matrix)
%  OUTPUTS: thetas_sample is mxT
%           theta0_sample is mx1
%  Author: Francesco Corsello, May/June 2014
[m,T]=size(theta_pred);
if nargin < 5
    Gamma=eye(m);
    gammaidn_bool=true;
else
    gammaidn_bool=false;
end


theta_sample=zeros(m,T+1);

theta_sample(:,T+1)=transpose(mvnrnd( theta_upd(:,T+1)' ,...
    correct_precise_eig(make_sim_precise( P_upd(:,:,T+1)))));
for t=T:-1:1
    if gammaidn_bool
        P_upd_invpred=P_upd(:,:,t)/P_pred(:,:,t);
        theta_smooth=theta_upd(:,t)+P_upd_invpred*...
            (theta_sample(:,t+1)-theta_upd(:,t));
        P_smooth=(eye(m)- P_upd_invpred )* P_upd(:,:,t) ;
    else
        theta_smooth=theta_upd(:,t)+...
            (  P_upd(:,:,t)* Gamma' ) *  (P_pred(:,:,t)\...
            (theta_sample(:,t+1)-Gamma*theta_upd(:,t)));
        P_smooth=P_upd(:,:,t)*(eye(m)-  Gamma'*  (P_pred(:,:,t)\Gamma)  *   P_upd(:,:,t));
    end
    
    try
        
        P_smooth_prec=make_sim_precise(P_smooth);
        sampled=sample_multinorm(  theta_smooth' , P_smooth_prec      );
        theta_sample(:,t)=   sampled';
        
    catch
        theta_sample(:,t)=transpose(sample_multinorm(  theta_smooth' , ...
            correct_precise_eig(make_sim_precise(P_smooth)  )));
    end

end

thetas_sample=theta_sample(:,2:T+1);
theta0_sample=theta_sample(:,1);
