function [Yact,Yhat,Yf1,Yf0] = predictfnc_ar1_ME(datmat,xvec)
%Read in data
Y=datmat;
%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=6; %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]); 

%% Step 1: Run a filter to recover k1 k2

rho=xvec(1);
sigw=xvec(2);
sigv=xvec(3);
rhohat=xvec(4);

num = (rho-rhohat)^2;
denom = (1-rho^2);
Veps = sigw^2*(num/denom + 1);
sigeps=sqrt(Veps);

%% First, run through filter to get k1 and k2
T=length(shocks.w);
s=zeros(1,T);
for t=2:T
    s(1,t) = rho*s(1,t-1) + sigw*shocks.w(1,t);
end
y=s + sigv*shocks.v;

AA = rhohat; 
BB = sigeps;
CC = 1;
DD = sigv;
s0=zeros(1,1); 
P0=eye(1);

% Run Kalman filter to get k1 
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(1)-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(1)-K*CC)*Pbm1;
end 
k1=K(1,1); 

%% Next run through another filter to construct fitted values

a = [rho 0 0 0 0 0;
    1 0 0 0 0 0;
    0 0 0 0 0 0;
    k1*rho 0 0 (1-k1)*rhohat 0 0;
    0 0 0 rhohat 0 0;
    0 0 0 1 0 0];

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

c = [rho 0 1 -rhohat 0 0; 
    0 0 0 rhohat -rhohat 0; 
    0 1 0 0 0 -1];

cf0 = [0 0 0 1 0 0 ; 0 0 0 0 0 0; 0 0 0 0 0 0];
cf1 = [0 0 0 rhohat 0 0 ; 0 0 0 0 0 0; 0 0 0 0 0 0];

%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));

%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;
     
	%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};
        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};
        Yf0{t} = Wtf0*yttm1{t};
        Yf1{t} = Wtf1*yttm1{t};
    end
end %t

end