*! version 11.0.0  13Feb2018
program define mymean11, eclass

    version 15.1

    syntax varlist(numeric) [if] [in]
    marksample touse
    tempname b V N

    mata: mymean_work("`varlist'", "`touse'", "`b'", "`V'", "`N'")

    matrix colnames `b'  = `varlist'
    matrix colnames `V'  = `varlist'
    matrix rownames `V'  = `varlist'

    ereturn post `b' `V', esample(`touse')
    ereturn scalar   N   = `N'
    ereturn scalar df_r  = `N'-1
    ereturn display
end        

mata:
void mymean_work(string scalar vlist,          ///
    string scalar touse, string scalar bname,  ///
    string scalar vname, string scalar nname )
{
    real matrix X, E, V, samp
    real vector b
    real scalar n

    st_view(samp=., ., touse)
    st_view(X=., ., vlist)
    
    MyAve(X, samp, b, n)
    MyV(X, samp, b, V)

    st_matrix(bname, b)
    st_matrix(vname, V)
    st_numscalar(nname, n)
}

void MyAve(real matrix X, real vector samp, b, n)
{
    real scalar   r, c, i, j    

    r = rows(X)
    c = cols(X)
    b = J(1, c, 0)
    n = 0
    for(i=1; i<=r; i++) {
        if (samp[i]==1) {
            ++n
            for(j=1; j<=c; j++) {
                b[1,j] = b[1,j] + X[i,j]
            }
        }
    }
    b = (1/n)*b
}

void MyV(real matrix X, real vector samp, real matrix b, V)
{
    real scalar r, c, i, j, j2, n

    r = rows(X)
    c = cols(X)

    V = J(c, c, 0)
    e = J(1, c, 0)
    n = 0
    if (rows(b)!=1 | cols(b)!=c) {
        printf("{err}sample-average vector and VCE are not conformable\n")
        exit(error(503))
    }
    for(i=1; i<=r; i++) {
        if (samp[i]==1) {
            ++n
            for(j=1; j<=c; j++) {
                e[1,j] = X[i,j] - b[1,j]        
            }
            for(j=1; j<=c; j++) {
                for(j2=1; j2<=j; j2++) {
                    V[j, j2] = V[j, j2] + e[1,j]*e[1,j2] 
                }
            }
        }
    }
    for(j=1; j<=c; j++) {
        for(j2=j+1; j2<=c; j2++) {
            V[j, j2] = V[j2, j] 
        }
    }

    V = (1/n)*(1/(n-1))*V
}

end