*! version 10.0.0  13Feb2018
program define mymean10, 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 
    real vector b
    real scalar n

    X = st_data(., vlist, touse)
    b = mean(X)
    E = (X :- b)
    n = rows(E)
    V = (1/n)*(1/(n-1))*quadcross(E,E)

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