
% Function that evaluates model moments at the estimated parameter vector

function m = estmoms(thetahat,k,mhat_scf,sa1,mhat_scf_param4)

    % Extract estimated parameters or sample moments from SCF
    if k==20
        mu_r1       = thetahat(1);
        mu_r2       = thetahat(2);
        mu_a1       = thetahat(3);
        mu_a2       = thetahat(4);
        var_r1      = thetahat(5);
        var_r2      = thetahat(6);
        var_a1      = thetahat(7);
        var_a2      = thetahat(8);
        cov_r1r2    = thetahat(9);
        cov_r1a1    = thetahat(10);
        cov_r1a2    = thetahat(11);
        cov_r2a1    = thetahat(12);
        cov_r2a2    = thetahat(13);
        cov_a1a2    = thetahat(14);
        pi_i1       = thetahat(15);
        pi_i2       = thetahat(16);
        pi_c1       = thetahat(17);
        pi_c2       = thetahat(18);
        u_i         = thetahat(19);
        u_c         = thetahat(20);
    elseif k==18
        mu_r1       = thetahat(1);
        mu_r2       = thetahat(2);
        mu_a1       = thetahat(3);
        mu_a2       = thetahat(4);
        var_r1      = thetahat(5);
        var_r2      = thetahat(6);
        var_a1      = thetahat(7);
        var_a2      = thetahat(8);
        cov_r1r2    = thetahat(9);
        cov_r1a1    = thetahat(10);
        cov_r1a2    = thetahat(11);
        cov_r2a1    = thetahat(12);
        cov_r2a2    = thetahat(13);
        cov_a1a2    = thetahat(14);
        pi_i1       = thetahat(15);
        pi_i2       = thetahat(16);
        pi_c1       = thetahat(17);
        pi_c2       = thetahat(18);
        u_i         = 0;
        u_c         = 0;
    elseif k==19
        mu_r1       = thetahat(1);
        mu_r2       = thetahat(2);
        mu_a1       = thetahat(3);
        mu_a2       = thetahat(4);
        var_r1      = thetahat(5);
        var_r2      = thetahat(6);
        var_a1      = thetahat(7);
        var_a2      = thetahat(8);
        cov_r1r2    = thetahat(9);
        cov_r1a1    = thetahat(10);
        cov_r1a2    = thetahat(11);
        cov_r2a1    = thetahat(12);
        cov_r2a2    = thetahat(13);
        cov_a1a2    = thetahat(14);
        pi_i1       = thetahat(15);
        pi_i2       = thetahat(16);
        pi_c1       = thetahat(17);
        pi_c2       = thetahat(18);
        u_i         = 0;
        u_c         = 0;
        sa1         = thetahat(19);
    elseif k==12 
        mu_r1       = thetahat(1);
        mu_r2       = thetahat(2);
        mu_a1       = thetahat(3);
        mu_a2       = thetahat(4);
        var_r1      = thetahat(5);
        var_r2      = thetahat(6);
        var_a1      = thetahat(7);
        var_a2      = thetahat(8);
        cov_r1r2    = mhat_scf(1);
        cov_r1a1    = mhat_scf(2);
        cov_r1a2    = mhat_scf(3);
        cov_r2a1    = mhat_scf(4);
        cov_r2a2    = mhat_scf(5);
        cov_a1a2    = mhat_scf(6);
        pi_i1       = thetahat(9);
        pi_i2       = thetahat(10);
        pi_c1       = thetahat(11);
        pi_c2       = thetahat(12);
        u_i         = 0;
        u_c         = 0;
   elseif k==4 
        mu_r1       = mhat_scf_param4(1);
        mu_r2       = mhat_scf_param4(2);
        mu_a1       = mhat_scf_param4(3);
        mu_a2       = mhat_scf_param4(4);
        var_r1      = mhat_scf_param4(5);
        var_r2      = mhat_scf_param4(6);
        var_a1      = mhat_scf_param4(7);
        var_a2      = mhat_scf_param4(8);
        cov_r1r2    = mhat_scf(1);
        cov_r1a1    = mhat_scf(2);
        cov_r1a2    = mhat_scf(3);
        cov_r2a1    = mhat_scf(4);
        cov_r2a2    = mhat_scf(5);
        cov_a1a2    = mhat_scf(6);
        pi_i1       = thetahat(1);
        pi_i2       = thetahat(2);
        pi_c1       = thetahat(3);
        pi_c2       = thetahat(4);
        u_i         = 0;
        u_c         = 0;    
    end
    
    % Evaluate model moments at estimated parameter vector
    mu_y1       = mu_r1+mu_a1;
    mu_y2       = mu_r2+mu_a2;
    mu_a        = sa1*mu_a1+(1-sa1)*mu_a2;
    mu_ri       = pi_i1*mu_r1+pi_i2*mu_r2+u_i;
    mu_rc       = pi_c1*mu_r1+pi_c2*mu_r2+u_c;
    var_y1      = var_r1+var_a1+2*cov_r1a1;
    var_y2      = var_r2+var_a2+2*cov_r2a2;
    var_a       = (sa1^2)*var_a1+((1-sa1)^2)*var_a2+2*sa1*(1-sa1)*cov_a1a2;
    var_ri      = (pi_i1^2)*var_r1+(pi_i2^2)*var_r2+2*pi_i1*pi_i2*cov_r1r2;
    var_rc      = (pi_c1^2)*var_r1+(pi_c2^2)*var_r2+2*pi_c1*pi_c2*cov_r1r2;
    cov_y1y2    = cov_r1r2+cov_r1a2+cov_r2a1+cov_a1a2;
    cov_y1a     = sa1*cov_r1a1+sa1*var_a1+...
                (1-sa1)*cov_r1a2+(1-sa1)*cov_a1a2;
    cov_y1ri    = pi_i1*var_r1+pi_i1*cov_r1a1+...
                pi_i2*cov_r1r2+pi_i2*cov_r2a1;
    cov_y1rc    = pi_c1*var_r1+pi_c1*cov_r1a1+...
                pi_c2*cov_r1r2+pi_c2*cov_r2a1;
    cov_y2a     = sa1*cov_r2a1+sa1*cov_a1a2+...
                (1-sa1)*cov_r2a2+(1-sa1)*var_a2;
    cov_y2ri    = pi_i1*cov_r1r2+pi_i1*cov_r1a2+...
                pi_i2*var_r2+pi_i2*cov_r2a2;
    cov_y2rc    = pi_c1*cov_r1r2+pi_c1*cov_r1a2+...
                pi_c2*var_r2+pi_c2*cov_r2a2;
    cov_ari     = pi_i1*sa1*cov_r1a1+pi_i1*(1-sa1)*cov_r1a2+...
                pi_i2*sa1*cov_r2a1+pi_i2*(1-sa1)*cov_r2a2;
    cov_arc     = pi_c1*sa1*cov_r1a1+pi_c1*(1-sa1)*cov_r1a2+...
                pi_c2*sa1*cov_r2a1+pi_c2*(1-sa1)*cov_r2a2;
    cov_rirc    = pi_i1*pi_c1*var_r1+(pi_i1*pi_c2+pi_i2*pi_c1)*cov_r1r2+...
                pi_i2*pi_c2*var_r2;
            
    % Vector of model moments evaluated at the estimated parameter vector
    m           = [mu_y1;mu_y2;mu_a;mu_ri;mu_rc;var_y1;cov_y1y2;...
                cov_y1a;cov_y1ri;cov_y1rc;var_y2;cov_y2a;cov_y2ri;...
                cov_y2rc;var_a;cov_ari;cov_arc;var_ri;cov_rirc;var_rc];
