function [Yact,Yhat,Yf1,Yf0,Ycq] = predictfnc_RE(datmat,xvec)
%Read in data
Y=datmat;
%Rearrange dimensions to get NxTx3-dimensional Y

%draw shocks
rng(2501);
shocks.v=randn(1,1000);
shocks.w=randn(1,1000);

N=size(Y,1);    %number of forecasters
Nper=size(Y,2); %number of periods

%within forcaster-date settings:
Nz=3; %number of observations
Ns=7; %number of state variables

%reshape the data
fe=reshape(Y(:,:,1),[N*Nper,1]);
fr=reshape(Y(:,:,2),[N*Nper,1]);
lfe=reshape(Y(:,:,3),[N*Nper,1]);
yit=reshape([fe fr lfe]',[Nz*N,Nper]); 

%% Run a filter to recover k1 k2
rho1=xvec(1);
rho2=xvec(2);
sigw=xvec(3);
sigv=xvec(4);

T=length(shocks.w);
s=zeros(1,T);
for t=3:T
    s(1,t) = rho1*s(1,t-1) + rho2*s(1,t-2) + sigw*shocks.w(1,t);
end
y = s + sigv*shocks.v;
AA = [rho1 rho2; 1 0]; 
BB = [sigw;0];
CC = [1 0];
DD = sigv;
s0 = zeros(2,1); 
P0 = eye(2);

% Run Kalman filter to get k1 k2
t=1;
    %predict
    sbm1 = AA*s0;
    Pbm1 = AA*P0*AA' + BB*BB';
    %update
    z=y(1,t);
    K = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
    sb = sbm1 + K*(z-CC*sbm1);
    Pb = (eye(2)-K*CC)*Pbm1;
for t=2:T
    %predict
     sbm1 = AA*sb;
     Pbm1 = AA*Pb*AA' + BB*BB';
       
    %update
    z=y(1,t);
    K = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
    sb = sbm1 + K*(z-CC*sbm1);
    Pb = (eye(2)-K*CC)*Pbm1;
end 
k1=K(1,1); 
k2=K(2,1);

%% Next run through another filter to construct fitted values
a = [rho1 rho2 0 0 0 0 0;
    1 0 0 0 0 0 0;
    0 0 0 0 0 0 0;
    k1*rho1 k1*rho2 0 (1-k1)*rho1 (1-k1)*rho2 0 0;
    k2*rho1 k2*rho2 0 1-k2*rho1 -k2*rho2 0 0;
    0 0 0 rho1 rho2 0 0;
    0 0 0 1 0 0 0];

b = [0 1 0;
    0 0 0;
    1 0 0;
    0 k1 k1;
    0 k2 k2;
    0 0 0;
    0 0 0];

% FE, FR, FEMN -- 1Q
c = [rho1 rho2 1 -rho1 -rho2 0 0;
    0 0 0 rho1 rho2 -rho1 -rho2;
    0 1 0 0 0 0 -1];

cf0 = [0 0 0 1 0 0 0; 0 0 0 0 0 0 0; 0 0 0 0 0 0 0];
cf1 = [0 0 0 rho1 rho2 0 0; 0 0 0 0 0 0 0; 0 0 0 0 0 0 0];
ccq = [1 0 0 -1 0 0 0;
    0 0 0 1 0 -1 0;
    0 1 0 0 0 0 -1];


%define relevant matrices (kron products are taken such that they are
%consistent with how yit is shaped)
T = kron(speye(N),a); 
B = kron(speye(N),b); 
Q = kron(speye(N),[sigw^2 0 0; 0 sigw^2 0; 0 0 sigv^2]);
W = (kron(speye(N),c));
Wf0  = (kron(speye(N),cf0));
Wf1 = (kron(speye(N),cf1));
Wcq = (kron(speye(N),ccq));

%Kalman Filter
P0=speye(N*Ns); 
F={}; K={}; 
Pttm1={}; Ptt={};
yttm1={}; ytt={};
y0=sparse(zeros(Ns*N,1)); 
Yact={}; Yhat={};

for t=1:Nper
	yt=yit(:,t);
	Mt=speye(Nz*N);
	Mt(isnan(yt),:)=[]; 
	yt(isnan(yt))=[];
	Wt=Mt*W;
    Wtf0 = Mt*Wf0;
    Wtf1 = Mt*Wf1;
    Wtcq = Mt*Wcq;
    
	%Run predict/update
    if t==1
        yttm1{t} = T*y0;
        Pttm1{t} = T*P0*T' + B*Q*B'; 
    
        F{t} = Wt*Pttm1{t}*Wt';
        K{t} = Pttm1{t}*Wt'/F{t};
    
        ytt{t} = yttm1{t} + K{t}*(yt - Wt*yttm1{t});
        Ptt{t} = Pttm1{t} - K{t}*Wt*Pttm1{t};
        
        Yact{t} = yt;
        Yhat{t} = Wt*yttm1{t};
        Ycq{t} = Wtcq*yttm1{t};
        Yf0{t} = Wtf0*yttm1{t};
        Yf1{t} = Wtf1*yttm1{t};
        
    else
        yttm1{t} = T*ytt{t-1};
        Pttm1{t} = T*Ptt{t-1}*T' + B*Q*B'; 
    
        F{t} = Wt*Pttm1{t}*Wt';
        K{t} = Pttm1{t}*Wt'/F{t};
    
        ytt{t} = yttm1{t} + K{t}*(yt - Wt*yttm1{t});
        Ptt{t} = Pttm1{t} - K{t}*Wt*Pttm1{t};
        
        Yact{t} = yt;
        Yhat{t} = Wt*yttm1{t};
        Ycq{t} = Wtcq*yttm1{t};
        Yf0{t} = Wtf0*yttm1{t};
        Yf1{t} = Wtf1*yttm1{t};
    end
end %t

end
