*! version 12.0.0  16Jan2016
program define myregress12, eclass sortpreserve
    version 14.1

    syntax varlist(numeric ts fv) [if] [in] [, noCONStant vce(string) ]
    marksample touse

    _vce_parse `touse' , optlist(Robust) argoptlist(CLuster) : , vce(`vce')
    local vce        "`r(vce)'"
    local clustervar "`r(cluster)'"
    if "`vce'" == "robust" | "`vce'" == "cluster" {
        local vcetype "Robust"
    }
    if "`clustervar'" != "" {
        capture confirm numeric variable `clustervar'
        if _rc {
            display in red "invalid vce() option"
            display in red "cluster variable {bf:`clustervar'} is " ///
                "string variable instead of a numeric variable"
            exit(198)
        }
        sort `clustervar'
    }

    gettoken depvar indepvars : varlist
    _fv_check_depvar `depvar'

    fvexpand `indepvars'
    local cnames `r(varlist)'

    tempname b V N rank df_r

    mata: mywork("`depvar'", "`cnames'", "`touse'", "`constant'",    ///
       "`vce'", "`clustervar'",                                  ///
       "`b'", "`V'", "`N'", "`rank'", "`df_r'")

    if "`constant'" == "" {
        local cnames `cnames' _cons
    }

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

    ereturn post `b' `V', esample(`touse') buildfvinfo
    ereturn scalar N        = `N'
    ereturn scalar rank     = `rank'
    ereturn scalar df_r     = `df_r'
    ereturn local  vce      "`vce'"
    ereturn local  vcetype  "`vcetype'"
    ereturn local  clustvar "`clustervar'"
    ereturn local  cmd      "myregress12"

    ereturn display

end

mata:

void mywork( string scalar depvar,  string scalar indepvars,
             string scalar touse,   string scalar constant,
             string scalar vcetype, string scalar clustervar,
             string scalar bname,   string scalar Vname,
             string scalar nname,   string scalar rname,
             string scalar dfrname)
{

    real vector    y, b, e, e2, cvar, ei
    real matrix    X, XpXi, M, info, xi
    real scalar    n, p, k, nc, i, dfr

    y    = st_data(., depvar, touse)
    X    = st_data(., indepvars, touse)
    n    = rows(X)

    if (constant == "") {
        X    = X,J(n,1,1)
    }

    XpXi = quadcross(X, X)
    XpXi = invsym(XpXi)
    b    = XpXi*quadcross(X, y)
    e    = y - X*b
    e2   = e:^2
    p    = cols(X)
    k    = p - diag0cnt(XpXi)
    if (vcetype == "robust") {
        M    = quadcross(X, e2, X)
        dfr  = n - k
        V    = (n/dfr)*XpXi*M*XpXi
    }
    else if (vcetype == "cluster") {
        cvar = st_data(., clustervar, touse)
        info = panelsetup(cvar, 1)
        nc   = rows(info)
        M    = J(k, k, 0)
        dfr  = nc - 1
        for(i=1; i<=nc; i++) {
            xi = panelsubmatrix(X,i,info)
            ei = panelsubmatrix(e,i,info)
            M  = M + xi'*(ei*ei')*xi
        }
        V    = ((n-1)/(n-k))*(nc/(nc-1))*XpXi*M*XpXi
    }
    else {                 // vcetype must IID
        dfr  = n - k
        V    = (quadsum(e2)/dfr)*XpXi
    }

    st_matrix(bname, b')
    st_matrix(Vname, V)
    st_numscalar(nname, n)
    st_numscalar(rname, k)
    st_numscalar(dfrname, dfr)

}

end