
% Function that computes sample moments from SCF

function mhat_scf = mhat_scf(ln_r1,ln_r2,ln_a1,ln_a2)

    % Matrix to extract covariance from variance-covariance matrix
    M           = [0,0;1,0];

    % Sample moments from SCF
    cov_r1r2    = sum(cov(ln_r1,ln_r2).*M,'all');
    cov_r1a1    = sum(cov(ln_r1,ln_a1).*M,'all');
    cov_r1a2    = sum(cov(ln_r1,ln_a2).*M,'all');
    cov_r2a1    = sum(cov(ln_r2,ln_a1).*M,'all');
    cov_r2a2    = sum(cov(ln_r2,ln_a2).*M,'all');
    cov_a1a2    = sum(cov(ln_a1,ln_a2).*M,'all');
    clear M

    % Vector of sample moments from SCF
    mhat_scf    = [cov_r1r2;cov_r1a1;cov_r1a2;cov_r2a1;cov_r2a2;cov_a1a2];
