[pmat,beta_tt]=....
    getfactors(dataFx,T,ns,H,R,F,Q,MU,beta0,P00)

%Carter and Kohn algorithm to draw the factor
beta_tt=zeros(T,ns);          %will hold the filtered state variable
ptt=zeros(ns,ns,T);    % will hold its variance
% %%%%%%%%%%%Step 6a run Kalman Filter
i=1;
x=H;
%Prediction
beta10=MU+beta0*F';
p10=F*P00*F'+Q;
yhat=(x*(beta10)')';                                                
eta=dataFx(i,:)-yhat;
feta=(x*p10*x')+R;
ifeta=invpd(feta);
%updating
K=(p10*x')*ifeta;
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);
beta_tt(i,:)=beta11;
ptt(:,:,i)=p11;

for i=2:T

    %Prediction
beta10=MU+beta11*F';
p10=F*p11*F'+Q;
yhat=(x*(beta10)')';  

eta=dataFx(i,:)-yhat;
feta=(x*p10*x')+R;
ifeta=invpd(feta);
%updating
K=(p10*x')*ifeta;
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);
ptt(:,:,i)=p11;
beta_tt(i,:)=beta11;

end
% Backward recursion to calculate the mean and variance of the distribution of the state
%vector
beta2 = zeros(T,ns);   %this will hold the draw of the state variable
jv1=1:ns; %index of state variables to extract
jv=1:ns;
f=F(jv,:);
q=Q(jv,jv);
mu=MU(jv);   
i=T;  %period t
p00=ptt(jv1,jv1,i); 
beta2(i,:)=beta_tt(i,:);
beta2(i,jv1)=mvnrnd(beta_tt(i:i,jv1),p00,1);%beta_tt(i:i,jv1)+(wa(i:i,jv1)*cholx(p00));   %draw for beta in period t from N(beta_tt,ptt)
%periods t-1..to .1
for i=T-1:-1:1


    
pt=ptt(:,:,i); 
ifptfq=invpd(f*pt*f'+q);
bm=beta_tt(i:i,:)+(pt*f'*ifptfq*(beta2(i+1:i+1,jv)-mu-beta_tt(i,:)*f')')';  
pm=pt-pt*f'*ifptfq*f*pt;  
beta2(i,:)=bm;
beta2(i:i,jv1)=mvnrnd(bm(jv1),pm(jv1,jv1),1);%bm(jv1)+(wa(i:i,jv1)*cholx(pm(jv1,jv1)));  
end


pmat=beta2(:,jv1);   %update the factors