function [smW,smPU,sxi] = smoothin(F,P0,P1,xi0,xi1,T)
for i=1:T
    J(:,:,i) = P1(:,:,i)*F'*inv(P0(:,:,i+1));
end

% smoothed series %
    sxi(:,T) = xi1(:,T);
    sxi(:,T-1) = xi1(:,T)+J(:,:,T-1)*(sxi(:,T)-xi0(:,T));
    
for i=2:T-1
    sxi(:,T-i) = xi1(:,T-i)+J(:,:,T-i)*(sxi(:,T+1-i)-xi0(:,T+1-i));
end

% corresponding mean squared errors %
 SP(:,:,T) = P1(:,:,T);
 SP(:,:,T-1) = P0(:,:,T-1) + J(:,:,T-1)*(P1(:,:,T)-P0(:,:,T))*J(:,:,T-1)';
for i=2:T-1
    SP(:,:,T-i) = P1(:,:,T-i)+J(:,:,T-i)*(SP(:,:,T+1-i)-P0(:,:,T+1-i))*J(:,:,T-i)';
end

% smoothed inflows and continuation probabilities
for t=1:T
smW(t,1) = sxi(1,t)';
smW(t,2) = sxi(2,t)';
smPU(t,1) = exp(-exp(sxi(3,t)));
smPU(t,2) = exp(-exp(sxi(4,t)));
end
end
