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

mu = param(1);
rho1 = param(2);
rho2 = param(3);
sigw = param(4);
sigv = param(5);
alphav = 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 = alphav*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_OC = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
        sb = sbm1 + K_OC*(z-CC*sbm1);
        Pb = (eye(2)-K_OC*CC)*Pbm1;

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

            %update
            z=y(n,t);
            K_OC = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
            sb = sbm1 + K_OC*(z-CC*sbm1);
            Pb = (eye(2)-K_OC*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
    
    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);
    
    ftp1t = rho1*ftt   + rho2*ftm1t;
    ftp2t = rho1*ftp1t + rho2*ftt;
    ftp3t = rho1*ftp2t + rho2*ftp1t;
    
    ftp1tm1 = rho1*fttm1   + rho2*ftm1tm1;
    ftp2tm1 = rho1*ftp1tm1 + rho2*fttm1;
    ftp3tm1 = rho1*ftp2tm1 + rho2*ftp1tm1;
    
    fe0 = s - ftt;
    fe1 = s(2:end) - ftp1t(:,1:end-1);
    fe2 = s(3:end) - ftp2t(:,1:end-2);
    fe3 = s(4:end) - ftp3t(:,1:end-3);

    fr0 = ftt - fttm1;
    fr1 = ftp1t(:,1:end-1) - ftp1tm1(:,1:end-1);
    fr2 = ftp2t(:,1:end-2) - ftp2tm1(:,1:end-2);
    fr3 = ftp3t(:,1:end-3) - ftp3tm1(:,1:end-3);
   
    sqq = ((1+s./100).^0.25-1)*100;
    fttqq=((1+ftt./100).^0.25-1)*100;
    ftp1tqq=((1+ftp1t./100).^0.25-1)*100;
    ftp2tqq=((1+ftp2t./100).^0.25-1)*100;
    ftp3tqq=((1+ftp3t./100).^0.25-1)*100;
    
    fttm1qq=((1+fttm1./100).^0.25-1)*100;
    ftp1tm1qq=((1+ftp1tm1./100).^0.25-1)*100;
    ftp2tm1qq=((1+ftp2tm1./100).^0.25-1)*100;
    ftp3tm1qq=((1+ftp3tm1./100).^0.25-1)*100;
    
    syoy = ((1+sqq(4:end)/100) .* (1+sqq(3:end-1)/100) .* (1+sqq(2:end-2)/100) .* (1+sqq(1:end-3)/100) - 1)*100;
    fyoy = ((1+fttqq(:,1:end-3)/100) .* (1+ftp1tqq(:,1:end-3)/100) .* (1+ftp2tqq(:,1:end-3)/100) .* (1+ftp3tqq(:,1:end-3)/100) -1 )*100;
    fyoym1 = ((1+fttm1qq(:,1:end-3)/100) .* (1+ftp1tm1qq(:,1:end-3)/100) .* (1+ftp2tm1qq(:,1:end-3)/100) .* (1+ftp3tm1qq(:,1:end-3)/100) -1 ) *100;

    fryoy = fyoy - fyoym1;
    feyoy = syoy - fyoy;
    
    %BMGS
    tmp_bgms=nancov(feyoy(:),fryoy(:));
    bgms(simct) = tmp_bgms(1,2)/tmp_bgms(2,2);

    %CG
    tmp_cg=nancov(mean(feyoy),mean(fryoy));
    cg(simct) = tmp_cg(1,2)/tmp_cg(2,2);

    %Disagree
    disagree(simct,:) = [std(ftt(:)) std(ftp1t(:)) std(ftp2t(:)) std(ftp3t(:))];

    %Updating rules
    Y = reshape(ftt,[N*T,1]);
    X = [reshape(fttm1,[N*T,1]) reshape(repmat(s,[N,1]),[N*T,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,[N*T,1]);
    TFE =repelem(eye(T),N,1);
    X = [reshape(fttm1,[N*T,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
