function momvec = SMM_model(xval,shocks,datamomvec,N,T,burn)
%% Set some parameters
H=8;%number of horizons 
Ns=2;%number of state variables in LGSS model
%% Read in xval
mu=2.4;
rho=xval(1);
sigw=xval(2);
sige=xval(3);
sigv=xval(4);
q=xval(5);
p=xval(6);
%% Generate DGP and signals
s=zeros(1,T+burn);
s(1,1)=mu;
for t=2:T+burn
    s(1,t) = (1-rho)*mu + rho*s(1,t-1) + sigw*shocks.w(1,t);
end
x=s+sige*shocks.e;
xm1=x;
s=s(1,burn+1:end);
x=x(1,burn+1:end);
xm1=xm1(1,burn:end-1);
y=repmat(s,[N,1]) + sigv*shocks.v;
qdraw=shocks.qdraw<q;
pdraw=shocks.pdraw<p;
%% Run Kalman filter
MU = [mu*(1-rho);0];
A = [rho 0; 1 0]; 
B = [sigw;0];
C = eye(Ns);
D = [sigv 0; 0 sige];

v_s = sigw^2/(1-rho^2);
s0=mu*ones(2,1); 
P0=[v_s rho*v_s; rho*v_s v_s] ;

%arrays that store KF predictions
stt=NaN(N,T,H);
sttm1=NaN(N,T,H);

% Loop for Kalman filter
for n=1:N
    
	t=1;
	%predict
    sbm1 = MU + A*s0;
    Pbm1 = A*P0*A' + B*B';
    %update
    z=[y(n,t);xm1(t)];
    K = Pbm1*C'/(C*Pbm1*C'+D*D');
    sb = sbm1 + K*(z-C*sbm1);
    Pb = (eye(Ns)-K*C)*Pbm1;
    %store forecasts
    sttm1(n,t,1) = sbm1(1);
    stt(n,t,1) = sb(1);
    for h=2:H
        sttm1(n,t,h) = (1-rho)*mu + rho*sttm1(n,t,h-1);
        stt(n,t,h) = (1-rho)*mu + rho*stt(n,t,h-1);
    end
    
    %Now do the rest
	for t=2:T
        %predict
        sbm1 = MU + A*sb;
	    Pbm1 = A*Pb*A' + B*B';
       
        %update
        z=[y(n,t);xm1(1,t)];
        K = Pbm1*C'/(C*Pbm1*C'+D*D');
        sb = sbm1 + K*(z-C*sbm1);
        Pb = (eye(Ns)-K*C)*Pbm1;
            
        %store forecasts 
        sttm1(n,t,1) = sbm1(1);
        stt(n,t,1) = sb(1);
        for h=2:H
            sttm1(n,t,h) = (1-rho)*mu + rho*sttm1(n,t,h-1);
            stt(n,t,h) = (1-rho)*mu + rho*stt(n,t,h-1);
        end

	end %t
        
end %n
%% Now build reported forecasts, updates
xhat = NaN(N,T,H);
xhat1a=NaN(N,T);
xhat2a=NaN(N,T);

%initialize
xhat(:,1:4,:) = stt(:,1:4,:);
xhat1a(:,4) = mean(stt(:,1,1:4),3);
xhat2a(:,4) = mean(stt(:,1,5:8),3);

for n=1:N
    for t=5:T
        if mod(t,4)==1
            xhataimplied=mean(stt(n,t,1:4));
            xhatimplied2a=mean(stt(n,t,5:8));
        end
        if mod(t,4)==2
            xhataimplied=(x(1,t-1) + sum(stt(n,t,1:3)) )/4;
            xhatimplied2a=mean(stt(n,t,4:7));
        end
        if mod(t,4)==3
            xhataimplied=(x(1,t-2) + x(1,t-1) + sum(stt(n,t,1:2)) )/4;
            xhatimplied2a=mean(stt(n,t,3:6));
        end
        if mod(t,4)==0
            xhataimplied=(x(1,t-3) + x(1,t-2) + x(1,t-1) + stt(n,t,1) )/4;
            xhatimplied2a=mean(stt(n,t,2:5));
        end
        
        if qdraw(n,t)==0 && pdraw(n,t)==0 %Case 1: No updating
            
            xhat(n,t,1:H-1) = xhat(n,t-1,2:H);
            if mod(t,4)==0
                xhat1a(n,t) = xhat2a(n,t-1);
                xhat2a(n,t) = NaN;
            else
                xhat1a(n,t) = xhat1a(n,t-1);
                xhat2a(n,t) = xhat2a(n,t-1);
            end
            
        elseif qdraw(n,t)==1 && pdraw(n,t)==0 %Case 2: Update Q not A
            if mod(t,4)==1
                xhat(n,t,1:4) = stt(n,t,1:4) + (xhat1a(n,t-1) - xhataimplied); 
                xhat(n,t,5) = stt(n,t,5) + (xhat2a(n,t-1) - xhatimplied2a);
            end
            if mod(t,4)==2
                xhat(n,t,1:3) = stt(n,t,1:3) + (xhat1a(n,t-1) - xhataimplied); 
                xhat(n,t,4:5) = stt(n,t,4:5) + (xhat2a(n,t-1) - xhatimplied2a); 
            end
            if mod(t,4)==3
                xhat(n,t,1:2) = stt(n,t,1:2) + (xhat1a(n,t-1) - xhataimplied); 
                xhat(n,t,3:5) = stt(n,t,3:5) + (xhat2a(n,t-1) - xhatimplied2a); 
            end
            if mod(t,4)==0
                xhat(n,t,1) = stt(n,t,1) + (xhat1a(n,t-1) - xhataimplied); 
                xhat(n,t,2:5) = stt(n,t,2:5) + (xhat2a(n,t-1) - xhatimplied2a); 
            end
            
            if mod(t,4)==0
                xhat1a(n,t) = xhat2a(n,t-1);
                xhat2a(n,t) = NaN;
            else
                xhat1a(n,t) = xhat1a(n,t-1);
                xhat2a(n,t) = xhat2a(n,t-1);
            end
            
        elseif qdraw(n,t)==0 && pdraw(n,t)==1 %Case 3: Update A not Q
            xhat(n,t,1:H-1) = xhat(n,t-1,2:H);
            if mod(t,4)==1 %if in Q1, yesterday's FE doesn't figure into A today
                xhat1a(n,t) = xhat1a(n,t-1);
            else %else if in Q2-Q4, yesterday's FE does figure into A today
                xhat1a(n,t) = xhat1a(n,t-1) + (1/4)*(xm1(1,t) - xhat(n,t-1,1));
            end
            xhat2a(n,t) = xhat2a(n,t-1);
    
        elseif qdraw(n,t)==1 && pdraw(n,t)==1 %Case 4: Update both
            xhat(n,t,:) = stt(n,t,:);
            xhat1a(n,t) = xhataimplied;
            xhat2a(n,t) = xhatimplied2a;  
        end
        
    end %t
end %n
%% Construct moments
%clearing the first initialization periods 
%and getting the timing right across vars
fr = xhat(:,5:end,1) - xhat(:,4:end-1,2);
lfe = repmat(x(4:end-1),[N,1]) - xhat(:,4:end-1,1);
feq = repmat(x(5:end),[N,1]) - xhat(:,5:end,1);
xar = movmean(x,[3,0]);
fey = repmat(xar(5:end),[N,1]) - xhat1a(:,5:end);

xhat=xhat(:,5:end,1);
xhat1a=xhat1a(:,5:end);

%computing target moments
cov1=nancov(xhat,xhat1a);
cov2=nancov(fr,lfe);
momvec=zeros(length(datamomvec),1);
momvec(1) = cov1(1,1);%C(xhat,xhat1a)_[1,1]
momvec(2) = cov1(1,2);%C(xhat,xhat1a)_[1,2]
momvec(3) = cov1(2,2);%C(xhat,xhat1a)_[2,2]
momvec(4) = cov2(1,1);%C(fr,lfe)_[1,1]
momvec(5) = cov2(1,2);%C(fr,lfe)_[1,2]
momvec(6) = cov2(2,2);%C(fr,lfe)_[2,2]
momvec(7) = nanmean(feq.^2,'all');%MSEQ
momvec(8) = nanmean(fey.^2,'all');%MSEY

end   