function [moments, moments_sd, disagree, K_DE] = simulated_moments_DE(param,N,T,burn,Nsim,shocks)

mu     = param(1);
rho1   = param(2);
rho2   = param(3);
sigw   = param(4);
sigv   = param(5);
varphi = param(6);

bgms = NaN(Nsim,1);
cg   = NaN(Nsim,1);
disagree = NaN(Nsim,4);
for simct=1:Nsim

    s = zeros(1,T+burn);
    for t=3:T+burn
        s(1,t) = rho1*s(1,t-1) + rho2*s(1,t-2) + sigw*shocks.w(simct,1,t);
    end
    y = s + sigv*squeeze(shocks.v(simct,:,:));

    AA = [rho1 rho2; 1 0]; 
    BB = [sigw;0];
    CC = [1 0];
    DD = sigv;
    s0 = zeros(2,1); 
    P0 = eye(2);

    % Run Kalman filter
    ftt=NaN(N,T+burn);
    fttm1=NaN(N,T+burn);
    ftm1t=NaN(N,T+burn);
    ftm1tm1=NaN(N,T+burn);

    for n=1:N

        t=1;
        %predict
        sbm1 = AA*s0;
        Pbm1 = AA*P0*AA' + BB*BB';
        %update
        z=y(n,t);
        K_RE = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
        sb = sbm1 + K_RE*(z-CC*sbm1);
        Pb = (eye(2)-K_RE*CC)*Pbm1;

        for t=2:T+burn
            
            %predict
             sbm1 = AA*sb;
             Pbm1 = AA*Pb*AA' + BB*BB';

            %update
            z=y(n,t);
            K_RE = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
            sb = sbm1 + K_RE*(z-CC*sbm1);
            Pb = (eye(2)-K_RE*CC)*Pbm1;

            ftt(n,t) = sb(1);
            fttm1(n,t) = sbm1(1);
            ftm1t(n,t) = sb(2);
            ftm1tm1(n,t) = sbm1(2);
            
        end %t 

    end %n
    
    K_DE = K_RE*(1+varphi);
    
    s = s(burn+1:end);
    ftt = ftt(:,burn+1:end);
    fttm1=fttm1(:,burn+1:end);
    ftm1t=ftm1t(:,burn+1:end);
    ftm1tm1=ftm1tm1(:,burn+1:end);
    
    sqq = ((1+s./100).^0.25-1)*100;
    fttqq=((1+ftt./100).^0.25-1)*100;
    fttm1qq=((1+fttm1./100).^0.25-1)*100;
    ftm1tqq = ((1+ftm1t./100).^0.25-1)*100;
    ftm1tm1qq = ((1+ftm1tm1./100).^0.25-1)*100;
    
    ftt_DE = (1+varphi)*ftt - varphi*fttm1;
    ftm1t_DE = (1+varphi)*ftm1t - varphi*ftm1tm1;
    ftm1tm1_DE = ftt_DE(:,1:end-1);
    fttm1_DE = rho1*ftt_DE(:,1:end-1) + rho2*ftm1t_DE(:,1:end-1);
    ftt_DE = ftt_DE(:,2:end);
    ftm1t_DE = ftm1t_DE(:,2:end);
    
    fttqq_DE = (1+varphi)*fttqq - varphi*fttm1qq;
    ftm1tqq_DE = (1+varphi)*ftm1tqq - varphi*ftm1tm1qq;
    ftm1tm1qq_DE = fttqq_DE(:,1:end-1);
    fttm1qq_DE = rho1*fttqq_DE(:,1:end-1) + rho2*ftm1tqq_DE(:,1:end-1);
    fttqq_DE = fttqq_DE(:,2:end);
    ftm1tqq_DE = ftm1tqq_DE(:,2:end);

    ftp1t_DE = rho1*ftt_DE   + rho2*ftm1t_DE;
    ftp2t_DE = rho1*ftp1t_DE + rho2*ftt_DE;
    ftp3t_DE = rho1*ftp2t_DE + rho2*ftp1t_DE;

    ftp1tm1_DE = rho1*fttm1_DE   + rho2*ftm1tm1_DE;
    ftp2tm1_DE = rho1*ftp1tm1_DE + rho2*fttm1_DE;
    ftp3tm1_DE = rho1*ftp2tm1_DE + rho2*ftp1tm1_DE;

    ftp1tqq_DE = rho1*fttqq_DE   + rho2*ftm1tqq_DE;
    ftp2tqq_DE = rho1*ftp1tqq_DE + rho2*fttqq_DE;
    ftp3tqq_DE = rho1*ftp2tqq_DE + rho2*ftp1tqq_DE;
    
    ftp1tm1qq_DE = rho1*fttm1qq_DE   + rho2*ftm1tm1qq_DE;
    ftp2tm1qq_DE = rho1*ftp1tm1qq_DE + rho2*fttm1qq_DE;
    ftp3tm1qq_DE = rho1*ftp2tm1qq_DE + rho2*ftp1tm1qq_DE;
    
    syoy = ((1+sqq(5:end)/100) .* (1+sqq(4:end-1)/100) .* (1+sqq(3:end-2)/100) .* (1+sqq(2:end-3)/100) - 1)*100;
    fyoy_DE = ((1+fttqq_DE(:,1:end-3)/100) .* (1+ftp1tqq_DE(:,1:end-3)/100) .* (1+ftp2tqq_DE(:,1:end-3)/100) .* (1+ftp3tqq_DE(:,1:end-3)/100) -1 )*100;
    fyoym1_DE = ((1+fttm1qq_DE(:,1:end-3)/100) .* (1+ftp1tm1qq_DE(:,1:end-3)/100) .* (1+ftp2tm1qq_DE(:,1:end-3)/100) .* (1+ftp3tm1qq_DE(:,1:end-3)/100) -1 ) *100;
    
    feyoy = syoy - fyoy_DE;
    fryoy = fyoy_DE - fyoym1_DE;
    
    fe0 = s(2:end) - ftt_DE;
    fe1 = s(3:end) - ftp1t_DE(:,1:end-1);
    fe2 = s(4:end) - ftp2t_DE(:,1:end-2);
    fe3 = s(5:end) - ftp3t_DE(:,1:end-3);

    fr0 = ftt_DE - fttm1_DE;
    fr1 = ftp1t_DE(:,1:end-1) - ftp1tm1_DE(:,1:end-1);
    fr2 = ftp2t_DE(:,1:end-2) - ftp2tm1_DE(:,1:end-2);
    fr3 = ftp3t_DE(:,1:end-3) - ftp3tm1_DE(:,1:end-3);

    tmp_bgms=nancov(feyoy(:),fryoy(:));
    bgms(simct) = tmp_bgms(1,2)/tmp_bgms(2,2);

    tmp_cg=nancov(mean(feyoy),mean(fryoy));
    cg(simct) = tmp_cg(1,2)/tmp_cg(2,2);
    
    disagree(simct,:) = [std(ftt_DE(:)) std(ftp1t_DE(:)) std(ftp2t_DE(:)) std(ftp3t_DE(:))];
    
    %Updating rules
    Y = reshape(ftt_DE,[N*(T-1),1]);
    X = [reshape(fttm1_DE,[N*(T-1),1]) reshape(repmat(s(2:end),[N,1]),[N*(T-1),1])];    
    tmp_update=mvregress(X,Y);
    beta1_update(simct) = tmp_update(1);
    beta2_update(simct) = tmp_update(2);
    beta3_update(simct) = 0;

    %With TFE    
    Y = reshape(ftt_DE,[N*(T-1),1]);
    TFE =repelem(eye(T-1),N,1);
    X = [reshape(fttm1_DE,[N*(T-1),1]) TFE];
    tmp_update=mvregress(X,Y);
    beta1_update_tfe(simct) = tmp_update(1);
    beta2_update_tfe(simct) = 0;
    beta3_update_tfe(simct) = 0;

end %simct

moments    = [mean(bgms) mean(cg)  mean(beta1_update) mean(beta2_update) mean(beta1_update_tfe)];
moments_sd = [std(bgms)  std(cg)   std(beta1_update)  std(beta2_update)  std(beta1_update_tfe) ];
disagree = mean(disagree,1);

end