*! version 1.0.4  21dec2020  Ben Jann

local rc 0
capt findfile lmoremata.mlib
if _rc {
    di as error "the {bf:moremata} package is required; type {stata ssc install moremata, replace}"
    error 499
}

program robstat, eclass properties(svyb svyj)
    version 11
    if replay() {
        Display `0'
        exit
    }
    local version : di "version " string(_caller()) ":"
    Parse_opts `0' // returns lhs, options, statistics, jbtest, jbtest2, jbwald, cluster
    `version' _vce_parserun robstat, mark(over) bootopts(`cluster') ///
        jkopts(`cluster') wtypes(pw iw fw) noeqlist: ///
        `lhs', nose `statistics' `options'
    if "`s(exit)'" != "" {
        ereturn local cmdline `"robstat `0'"'
        if "`jbtest'"!="" {
            JBtest, `jbtest2' `jbwald' display
        }
        exit
    }
    Estimate `0' // returns diopts
    ereturn local cmdline `"robstat `0'"'
    if "`jbtest'"!="" {
        JBtest, `jbtest2' `jbwald'
    }
    Display, `diopts'
    if `"`e(IFvars)'"'!="" {
        local IFnote "(influence functions stored in: `e(IFvars)')"
        local i 0
        while (1) {
            local ++i
            local line: piece `i' 80 of "`IFnote'"
            if "`line'"=="" continue, break
            di as txt "`line'"
        }
        di
    }
end

program Parse_opts
    _parse comma lhs 0 : 0
    syntax [, nose CLuster(passthru) Statistics(str) swap ///
        JBtest JBtest2(str) WALD  * ]
    if `"`jbtest2'"'!="" local jbtest jbtest
    if "`jbtest'"!="" & `"`statistics'"'!="" {
        di as err "statistics() not allowed with jbtest"
        exit 198
    }
    if "`jbtest'"!="" & "`swap'"!="" {
        di as err "swap not allowed with jbtest"
        exit 198
    }
    if "`wald'"!="" & "`jbtest'"=="" {
        di as err "wald only allowed if jbtest is specified"
        exit 198
    }
    if "`jbtest'"!="" JBtest_args `jbtest2' // sets local statistics
    c_local lhs `"`lhs'"'
    c_local cluster `cluster'
    c_local jbtest `jbtest'
    if "`jbtest2'"!="" {
        c_local jbtest2 jbtest2(`jbtest2')
    }
    c_local jbwald `wald'
    if "`statistics'"!="" {
        c_local statistics statistics(`statistics')
    }
    c_local options `swap' `options'
end

program Display
    syntax [, noHEader noTABle Level(passthru) CILog * ]
    if "`cilog'"!="" {
        if c(stata_version)<15 {
            di as err "option cilog only allowed in Stata 15 or newer"
            exit 198
        }
    }
    _get_diopts diopts, `options'
    if `"`e(cmd)'"'!="robstat" {
        di as err "last robstat results not found"
        exit 301
    }
    if `"`level'"'=="" local level level(`e(level)')
    if `"`header'"'=="" {
        if c(stata_version)>=12 {
            nobreak {
                Display_ereturn local cmd "total" // mimick header of -total-
                capture noisily break {
                    _coef_table_header
                }
                Display_ereturn local cmd "robstat"
                if _rc exit _rc
            }
        }
        else _coef_table_header
        if `"`e(over)'"'!="" {
            _svy_summarize_legend
        }
        else di ""
    }
    if `"`table'"'=="" {
        if "`cilog'"!="" {
            mata: st_local("cilog", "cilog" * anyof(st_matrix("e(class)"),2))
        }
        nobreak {
            local depvar `"`e(depvar)'"'
            local modified 0
            if e(N_stats)==1 {
                if e(N_vars)>1 | `"`e(over)'"'!="" {
                    Display_ereturn local depvar `"`e(statistics)'"'
                    local modified 1
                }
            }
            capture noisily break {
                Display_table "`cilog'" `"`level'"' `"`options'"'
            }
            if `modified' {
                Display_ereturn local depvar `"`depvar'"'
            }
            if _rc exit _rc
        }
    }
    capt confirm matrix e(jbtest)
    if _rc==0 JBtest_display
end
prog Display_ereturn, eclass
    ereturn `0'
end
program Display_table
    args cilog level options
    if c(stata_version)>=15 {
        if "`cilog'"!="" {
            qui _coef_table, `level'
            tempname CI
            mata: robstat_cilog()
            _coef_table, nopvalue cimat(`CI') `level' `options'
            mat drop `CI'
            di as txt `"(`cilognote')"'
        }
        else {
            eret di, nopvalue `level' `options'
        }
    }
    else if c(stata_version)>=14 {
        eret di, nopvalue `level' `options'
    }
    else if c(stata_version)>=13 {
        quietly update
        if r(inst_ado)>=d(26jun2014) {
            eret di, nopvalue `level' `options'
        }
        else {
            _coef_table, cionly `level' `options'
        }
    }
    else if c(stata_version) >= 12 {
        _coef_table, cionly `level' `options'
    }
    else {
        eret di, `level' `options'
    }
end

program JBtest, eclass
    syntax [, jbtest2(str) WALD display ]
    JBtest_args `jbtest2'
    local vce   "`e(vce)'"
    local wtype "`e(vce)'"
    if "`vce'"=="" {
        if "`wtype'"!="" & "`wtype'"!="fweight" local wald wald
    }
    else if "`vce'"!="analytic" local wald wald
    local ntests: list sizeof jbtests
    local eqs: coleq e(b)
    local eqs: list uniq eqs
    if "`wald'"!="" {
        if e(df_r)<. {
            local jbtype F
            local jbtitle "Normality Tests (Wald F; df_r = `e(df_r)')"
        }
        else {
            local jbtype chi2
            local jbtitle "Normality Tests (Wald chi2)"
        }
        tempname jbtest tmp
        foreach eq of local eqs {
            local rown
            matrix `tmp' = J(`ntests', 3,.)
            local m 0
            foreach t of local jbtests {
                local ++m
                if "`t'"=="jbera" {
                    local rown `rown' JB
                    local exp (_b[`eq':skewness]=0) (_b[`eq':kurtosis]=3)
                }
                else if "`t'"=="moors" {
                    local rown `rown' MOORS
                    local exp (_b[`eq':SK25]=0) (_b[`eq':QW25]=1.23)
                }
                else if "`t'"=="mc" {
                    local rown `rown' MC
                    local exp (_b[`eq':MC]=0)
                }
                else if "`t'"=="lmc" {
                    local rown `rown' LMC
                    local exp (_b[`eq':LMC]=0.199)
                }
                else if "`t'"=="rmc" {
                    local rown `rown' RMC
                    local exp (_b[`eq':RMC]=0.199)
                }
                else if "`t'"=="lr" {
                    local rown `rown' LR
                    local exp (_b[`eq':LMC]=0.199) (_b[`eq':RMC]=0.199)
                }
                else if "`t'"=="mcl" {
                    local rown `rown' MC-L
                    local exp (_b[`eq':MC]=0) (_b[`eq':LMC]=0.199)
                }
                else if "`t'"=="mcr" {
                    local rown `rown' MC-R
                    local exp (_b[`eq':MC]=0) (_b[`eq':RMC]=0.199)
                }
                else if "`t'"=="mclr" {
                    local rown `rown' MC-LR
                    local exp (_b[`eq':MC]=0) (_b[`eq':LMC]=0.199) (_b[`eq':RMC]=0.199)
                }
                qui test `exp'
                matrix `tmp'[`m',1] = r(`jbtype'), r(df), r(p)
            }
            mat rown `tmp' = `rown'
            mat roweq `tmp' = `eq'
            mat `jbtest' = nullmat(`jbtest') \ `tmp'
        }
        mat coln `jbtest' = `jbtype' df "Prob>`jbtype'"
    }
    else {
        local nover = e(N_over)
        local jbtype chi2
        local jbtitle "Normality Tests"
        tempname jbtest tmp D V
        local k 0
        foreach eq of local eqs {
            local ++k
            local N = el(e(_N), mod(`k'-1,`nover')+1, 1)
            local rown
            matrix `tmp' = J(`ntests', 3,.)
            local m 0
            foreach t of local jbtests {
                local ++m
                if "`t'"=="jbera" {
                    local rown `rown' JB
                    local df 2
                    matrix `D' = (_b[`eq':skewness], _b[`eq':kurtosis]-3)'
                    matrix `V' = (6, 0) \ (0, 24)
                }
                else if "`t'"=="moors" {
                    local rown `rown' MOORS
                    local df 2
                    matrix `D' = (_b[`eq':SK25], _b[`eq':QW25]-1.23)'
                    matrix `V' = (1.84, 0) \ (0, 3.14)
                }
                else if "`t'"=="mc" {
                    local rown `rown' MC
                    local df 1
                    matrix `D' = (_b[`eq':MC])
                    matrix `V' = (1.25)
                }
                else if "`t'"=="lmc" {
                    local rown `rown' LMC
                    local df 1
                    matrix `D' = (_b[`eq':LMC]-0.199)
                    matrix `V' = (2.62)
                }
                else if "`t'"=="rmc" {
                    local rown `rown' RMC
                    local df 1
                    matrix `D' = (_b[`eq':RMC]-0.199)
                    matrix `V' = (2.62)
                }
                else if "`t'"=="lr" {
                    local rown `rown' LR
                    local df 2
                    matrix `D' = (_b[`eq':LMC]-0.199, _b[`eq':RMC]-0.199)'
                    matrix `V' = (2.62, -.0123) \ (-.0123, 2.62)
                }
                else if "`t'"=="mcl" {
                    local rown `rown' MC-L
                    local df 2
                    matrix `D' = (_b[`eq':MC], _b[`eq':LMC]-0.199)'
                    matrix `V' = (1.25, .323) \ (.323, 2.62)
                }
                else if "`t'"=="mcr" {
                    local rown `rown' MC-R
                    local df 2
                    matrix `D' = (_b[`eq':MC], _b[`eq':RMC]-0.199)'
                    matrix `V' = (1.25, .323) \ (.323, 2.62)
                }
                else if "`t'"=="mclr" {
                    local rown `rown' MC-LR
                    local df 3
                    matrix `D' = (_b[`eq':MC], _b[`eq':LMC]-0.199, _b[`eq':RMC]-0.199)'
                    matrix `V' = (1.25, .323, -.323) \  ///
                                 (.323, 2.62, -.0123) \ ///
                                 (-.323, -.0123, 2.62)
                }
                matrix `tmp'[`m',1] = `N'*`D''*invsym(`V')*`D', `df'
                matrix `tmp'[`m',3] = chi2tail(`df', `tmp'[`m',1])
            }
            mat rown `tmp' = `rown'
            mat roweq `tmp' = `eq'
            mat `jbtest' = nullmat(`jbtest') \ `tmp'
        }
        mat coln `jbtest' = `jbtype' df "Prob>`jbtype'"

    }
    ereturn matrix jbtest = `jbtest'
    ereturn local jbwald `wald'
    ereturn local jbtype `jbtype'
    ereturn local jbtitle "`jbtitle'"
    if "`display'"!="" JBtest_display
end
program JBtest_args
    if strtrim(`"`0'"')=="all" {
        local 0 jbera moors mclr mcl mcr mc lr lmc rmc
    }
    local jbtests
    foreach s of local 0 {
        local t = lower(`"`s'"')
        local l  = strlen(`"`t'"')
        // jbera
        if "`t'"==substr("jbera", 1, max(2,`l')) {
            local jbtests `jbtests' jbera
            local s1 skewness kurtosis
        }
        else if "`t'"==substr("moors", 1, max(2,`l')) {
            local jbtests `jbtests' moors
            local s2 sk25 qw25
        }
        else if "`t'"=="mc" {
            local jbtests `jbtests' mc
            local s3 mc
        }
        else if "`t'"==substr("lmc", 1, max(1,`l')) {
            local jbtests `jbtests' lmc
            local s4 lmc
        }
        else if "`t'"==substr("rmc", 1, max(1,`l')) {
            local jbtests `jbtests' rmc
            local s5 rmc
        }
        else if "`t'"=="lr" {
            local jbtests `jbtests' lr
            local s4 lmc
            local s5 rmc
        }
        else if inlist("`t'", "mcl", "mc-l") {
            local jbtests `jbtests' mcl
            local s3 mc
            local s4 lmc
        }
        else if inlist("`t'", "mcr", "mc-r") {
            local jbtests `jbtests' mcr
            local s3 mc
            local s5 rmc
        }
        else if inlist("`t'", "mclr", "mc-lr") {
            local jbtests `jbtests' mclr
            local s3 mc
            local s4 lmc
            local s5 rmc
        }
        else {
            di as err `"`s' not allowed in jbtest()"'
            exit 198
        }
    }
    if "`jbtests'"=="" {
        local jbtests jbera moors mclr
        local statistics skewness kurtosis sk25 qw25 mc lmc rmc
    }
    else {
        local statistics `s1' `s2' `s3' `s4' `s5'
    }
    c_local statistics `statistics'
    c_local jbtests `jbtests'
end

program JBtest_display
    local rspec
    local k_eq = e(k_eq)
    local njb  = rowsof(e(jbtest)) / `k_eq' - 1
    forv i = 1/`k_eq' {
        forv j = 1/`njb' {
            local rspec `rspec'&
        }
        local rspec `rspec'-
    }
    matlist e(jbtest), rspec(--`rspec') cspec(&%12s|%10.2f&%5.0g&%10.4f&) ///
        title(`e(jbtitle)')
end

program Estimate, eclass
    // syntax
    syntax varlist(numeric) [if] [in] [pw iw fw/], [          ///
        Statistics(str asis) swap                             ///
        over(varname numeric) Total                           ///
        vce(str) CLuster(varname) svy SVY2(str) nose          ///
        GENerate GENerate2(name) replace                      ///
        /// normality test (just for parsing)
        JBtest JBtest2(str) WALD                              ///
        /// optimization options for M-estimation
        TOLerance(real 1e-10) ITERate(integer `c(maxiter)')   ///
        /// kernel density estimation options
        Kernel(name) bw(name) Adaptive(int 2) n(int 512)      ///
        /// display options
        Level(cilevel) noHEader noTABle CILog *               ///
        ]
    if "`cilog'"!="" {
        if c(stata_version)<15 {
            di as err "option cilog only allowed in Stata 15 or newer"
            exit 198
        }
    }
    if `"`generate2'"'!=""   local generate generate
    else if "`generate'"!="" local generate2 _IF_
    if `"`jbtest'`jbtest2'"'!="" {
        JBtest_args `jbtest2' // sets local statistics
    }
    if `tolerance'<=0 {
            di as err  "tolerance() must be positive"
            exit 198
    }
    if `iterate'<=0 {
            di as err  "iterate() must be positive"
            exit 198
    }
    local estopts tolerance(`tolerance') iterate(`iterate')
    if "`kernel'"=="" local kernel epan2
    if "`bw'"=="" local bw dpi
    if `adaptive'<0 {
        di as err "adaptive() must be positive"
        exit 198
    }
    if `n'<=2 {
        di as err "n() must be > 2"
        exit 198
    }
    local estopts `estopts' kernel(`kernel') bw(`bw') adaptive(`adaptive') n(`n')

    // varlist
    local varlist: list uniq varlist // remove repeated varnames
    local ndepv: list sizeof varlist

    // statistics
    Parse_stats `statistics'    // returns stats, statslbl
    local nstats: list sizeof stats

    // over(), total
    if "`over'"=="" {
        if "`total'"!="" {
            di as err "total only allowed if over() is specified"
            exit 198
        }
    }

    // vce
    local vce0 `"`vce'"'
    if `"`cluster'"'!="" {
        if "`se'"!="" {
            di as err "cluster() and nose not both allowed"
            exit 198
        }
        if `"`svy'`svy2'"'!="" {
            di as err "cluster() and svy() not both allowed"
            exit 198
        }
        if `"`vce'"'!="" {
            di as err "cluster() and vce() not both allowed"
            exit 198
        }
        local vce `"cluster `cluster'"'
        local cluster
    }
    if `"`vce'"'!="" {
        if "`se'"!="" {
            di as err "vce() and nose not both allowed"
            exit 198
        }
        if `"`svy'`svy2'"'!="" {
            di as err "vce() and svy() not both allowed"
            exit 198
        }
        gettoken vce clustvar : vce
        if `"`vce'"'=="analytic" & `"`clustvar'"'=="" {
            local vce analytic
        }
        else if `:list sizeof clustvar'==1 & ///
            substr("cluster", 1, max(2, strlen(`"`vce'"')))==`"`vce'"' {
            local vce cluster
            local vcetype Robust
            gettoken clustvar : clustvar
            local clustopt cluster(`clustvar')
        }
        else {
            di as err "invalid vce()"
            exit 198
        }
    }
    else if "`se'"=="" {
        local vce analytic
    }

    // svy
    if `"`svy'`svy2'"'!="" {
        if "`weight'"!="" {
            di as err "weights not allowed with svy; supply weights to {help svyset}"
            exit 101
        }
        if "`se'"!="" {
            di as err "svy() and nose not both allowed"
            exit 198
        }
        local svy svy
        if `"`svy2'"'!="" {
            local svy2 `"subpop(`svy2')"'
        }
    }

    // display options
    local levelopt level(`level')
    _get_diopts diopts, `options'
    c_local diopts `header' `table' `levelopt' `cilog' `diopts'

    // sample and weights
    marksample touse
    markout `touse' `clustvar' `over'
    if "`svy'"!="" {
        if c(stata_version)>=14 {
            tempvar subpop exp
            _svy_setup `touse' `subpop' `exp', svy `svy2'
            local weight `"`r(wtype)'"'
            if "`weight'"=="" {
                drop `exp'
                local exp
            }
        }
        else {
            tempvar subpop
            _svy_setup `touse' `subpop', svy `svy2'
            local weight `"`r(wtype)'"'
            local exp `"`r(wvar)'"'
        }
        if `"`r(vce)'"'!="linearized" {
            di as err "option svy is only allowed if VCE is set to linearized; " ///
                `"use the {helpb svy} prefix command for `r(vce)' survey estimation"'
            exit 498
        }
        local svy_posts `"`r(poststrata)'"'
        local svy_postw `"`r(postweight)'"'
        if `"`svy_posts'"'!="" {
            if "`weight'"!="" local wexp `"[`weight' = `exp']"'
            tempvar exp
            svygen post double `exp' `wexp' if `touse', ///
                posts(`svy_posts') postw(`svy_postw')
            if "`weight'"=="" local weight pweight
        }
    }
    else local subpop `touse'
    if "`over'"!="" {
        capt assert ((`over'==floor(`over')) & (`over'>=0)) if `subpop'
        if _rc {
            di as err "variable in over() must be integer and nonnegative"
            exit 452
        }
        qui levelsof `over' if `subpop', local(overvals)
        local N_over: list sizeof overvals
        if "`total'"!="" {
            local ++N_over
        }
        local over_labels
        foreach overval of local overvals {
            local over_labels `"`over_labels' `"`: label (`over') `overval''"'"'
        }
        local over_labels: list clean over_labels
    }
    else {
        local total total
        local N_over 1
    }
    if "`weight'"!="" {
        capt confirm variable `exp'
        if _rc {
            local wexp `"= `exp'"'
            tempvar exp
            qui gen double `exp' `wexp'
        }
        local wgt `"[`weight' = `exp']"'
        if "`weight'"=="pweight" {
            local swgt `"[aw = `exp']"'
        }
        else {
            local swgt `"`wgt'"'
        }
    }
    else {
        local wgt
        local swgt
    }
    if "`svy'"!="" & `"`exp'"'!="" {
        su `touse' if `touse', meanonly
    }
    else {
        su `touse' `swgt' if `touse', meanonly
    }
    local N = r(N)
    if `N'==0 error 2000

    // check whether variance estimation is supported
    if "`se'"=="" {
        if "`weight'"=="fweight" {
            if "`vce'"=="analytic" & `"`vce0'"'=="" {
                di as txt "(variance estimation not supported with fweights)"
                local se nose
            }
            else if "`svy'`vce'"!="" {
                di as err "variance estimation not supported with fweights"
                exit 498
            }
        }
    }

    // generate: check whether variables already exist
    if "`se'"=="" & "`generate'"!="" & "`replace'"=="" {
        foreach v of local varlist {
            foreach overval in `overvals' `total' {
                foreach stat of local stats {
                    Make_IF_Name `generate2' `stat' `v' `ndepv' "`over'" ///
                        `overval' `stat' `nstats' // returns vname
                    confirm new variable `vname'
                }
            }
        }
    }

    // compute statistics
    tempname b btmp class cltmp aux auxtmp _N
    local uvars
    if `nstats'>1 {
        local k_eq   = `ndepv' * `N_over'
        local k_coef = `nstats'
    }
    else {
        if "`over'"!="" {
            local k_eq   = `ndepv'
            local k_coef = `N_over'
        }
        else {
            local k_eq   = 1
            local k_coef = `ndepv'
        }
    }
    if "`over'"!="" {
        mat `_N' = J(`N_over', 1, .)
    }
    else {
        mat `_N' = J(1, 1, .)
    }
    mat rown `_N' = `overvals' `total'
    mat coln `_N' = "N"
    foreach v of local varlist {
        local k 0
        foreach overval in `overvals' `total' {
            local ++k
            // select obs
            if `"`overval'"'=="total" local touse1 (1)
            else                      local touse1 (`over'==`overval')
            // equation label
            if `k_eq'>1 {
                if `ndepv'==1        local eq "`overval'"
                else if "`over'"=="" local eq "`v'"
                else if `nstats'==1  local eq "`v'"
                else                 local eq "`v'_`overval'"
            }
            else if "`over'"!="" & `nstats'>1 local eq "`overval'"
            else local eq
            // count obs
            if "`svy'"!="" & `"`exp'"'!="" {
                su `touse' if `touse' & `subpop' & `touse1', meanonly
            }
            else {
                su `touse' `swgt' if `touse' & `subpop' & `touse1', meanonly
            }
            mat `_N'[`k', 1] = r(N)
            // compute stats
            mat `btmp' = J(1, `nstats', .)
            mat coln `btmp' = `statslbl'
            mat `cltmp' = `btmp'
            if `nstats'==1 {
                if "`over'"!="" {
                    mat coln `btmp' = `overval'
                }
                else if `ndepv'>1 {
                    mat coln `btmp' = `v'
                }
            }
            mat coleq `btmp' = `eq'
            mat `auxtmp' = `btmp'
            local j 0
            foreach stat of local stats {
                local u
                if "`se'"=="" {
                    tempvar u
                    qui gen double `u' = 0 if `touse'
                    local uvars `uvars' `u'
                }
                local ++j
                Estimate_`stat' `wgt' if `touse' & `subpop' & `touse1', ///
                    v(`v') u(`u') `estopts'
                mat `btmp'[1, `j']   = r(b)
                mat `cltmp'[1, `j']  = r(class)
                mat `auxtmp'[1, `j'] = r(k)
            }
            mat `b'     = nullmat(`b'), `btmp'
            mat `class' = nullmat(`class'), `cltmp'
            mat `aux'   = nullmat(`aux'), `auxtmp'
        }
    }

    // compute standard errors
    if "`se'"=="" {
        tempname V
        if "`svy'"!="" {
            //qui svy, `svy2': total `uvars' if `touse'
            qui svy, `svy2': mean `uvars' if `touse'
        }
        else {
            //qui total `uvars' `wgt' if `touse', `clustopt'
            qui mean `uvars' `wgt' if `touse', `clustopt'
        }
        matrix `V' = e(V)
        local rank = e(rank)
        local df_r = e(df_r)
        if "`vce'"=="cluster" {
            local N_clust = e(N_clust)
        }
        if "`svy'"!="" {
            local svy_scalars
            foreach l in N_sub N_strata N_strata_omit singleton census N_pop ///
                         N_subpop N_psu N_poststrata stages {
                if e(`l')<. {
                    local svy_`l' = e(`l')
                    local svy_scalars `svy_scalars' `l'
                }
            }
            local svy_macros
            foreach l in prefix wtype wvar wexp singleunit strata psu fpc ///
                         poststrata postweight vce vcetype mse subpop adjust {
                local svy_`l' `"`e(`l')'"'
                local svy_macros `svy_macros' `l'
            }
            forv l=1/`svy_stages' {
                local svy_su`l' `"`e(su`l')'"'
                local svy_fpc`l' `"`e(fpc`l')'"'
                local svy_weight`l' `"`e(weight`l')'"'
                local svy_strata`l' `"`e(strata`l')'"'
                local svy_macros `svy_macros' su`l' fpc`l' weight`l' strata`l'
            }
            local svy_matrices
            foreach l in V_srs V_srssub V_srswr V_srssubwr _N_strata_single ///
                _N_strata_certain _N_strata _N_postsum _N_postsize {
                capt confirm matrix e(`l')
                if _rc==0 {
                    tempname svy_`l'
                    mat `svy_`l'' = e(`l')
                    local svy_matrices `svy_matrices' `l'
                }
            }
        }
    }

    // post results
    if "`se'"=="" {
        local coln: colfullnames `b'
        mat coln `V' = `coln'
        mat rown `V' = `coln'
    }
    else if "`weight'"=="iweight" { // add empty e(V) for svy
        tempname V
        mat `V' = `b'' * `b' * 0
    }
    if "`swap'"!="" {   // flip equations and coefficients
        mata: _robstat_swap_eq_and_coefs(`k_eq', `k_coef')
    }
    eret post `b' `V' `wgt', obs(`N') esample(`touse')
    eret local cmd "robstat"
    eret local depvar "`varlist'"
    if "`svy'"!="" eret local title "Survey: Robust Statistics"
    else           eret local title "Robust Statistics"
    eret scalar k_eq = `k_eq'
    eret scalar N_stats = `nstats'
    eret local statistics "`statslbl'"
    eret scalar N_vars = `ndepv'
    if "`over'"!="" {
        eret local total "`total'"
        eret local over_labels `"`over_labels'"'
        eret local over_namelist `"`overvals'"'
        eret local over "`over'"
    }
    eret scalar N_over = `N_over'
    eret mat _N = `_N'
    eret mat class = `class'
    eret mat aux = `aux'
    if "`se'"=="" {
        eret local vcetype "`vcetype'"
        eret local vce "`vce'"
        eret scalar rank = `rank'
        eret scalar df_r = `df_r'
        eret scalar level = `level'
        if "`vce'"=="cluster" {
            eret local clustvar "`clustvar'"
            eret scalar N_clust = `N_clust'
        }
        if "`svy'"!="" {
            foreach l of local svy_scalars {
                eret scalar `l' = `svy_`l''
            }
            foreach l of local svy_macros {
                eret local `l' `"`svy_`l''"'
            }
            foreach l of local svy_matrices {
                if substr("`l'", 1, 1)=="V" {
                    mat coln `svy_`l'' = `coln'
                    mat rown `svy_`l'' = `coln'
                }
                eret matrix `l' = `svy_`l''
            }
        }
    }

    // generate: rename uvars
    if "`se'"=="" & "`generate'"!="" {
        local IFvars
        local i 0
        foreach v of local varlist {
            foreach overval in `overvals' `total' {
                foreach stat of local stats {
                    Make_IF_Name `generate2' `stat' `v' `ndepv' "`over'" ///
                        `overval' `stat' `nstats' // returns vname
                    capt confirm new variable `vname'
                    if _rc drop `vname'
                    local ++i
                    local uvar: word `i' of `uvars'
                    rename `uvar' `vname'
                    local IFvars `IFvars' `vname'
                }
            }
        }
        eret local IFvars `IFvars'
    }
    Return_clear // for some reason -return clear- does not delete r(S)
end
program Return_clear, rclass
    local x
end

program Make_IF_Name
    args vname stat v ndepv over overval stat nstats
    local stat: subinstr local stat ":" ""
    local uscore
    if `ndepv'>1 {
        local vname `vname'`v'
        local uscore _
    }
    if "`over'"!="" {
        local vname `vname'`uscore'`overval'
        local uscore _
    }
    if `nstats'>1 | "`uscore'"=="" {
        local vname `vname'`uscore'`stat'
    }
    c_local vname `vname'
end

program Parse_stats
    local stats
    local statslbl
    foreach s of local 0 {
        local ok = regexm(lower(`"`s'"'), "^([a-z]+)([0-9]*)$")
        if `ok' {
            local s1 = regexs(1)
            local s2 = regexs(2)
        }
        else {
            di as err `"`s' not allowed in statistics()"'
            exit 198
        }
        local l = strlen(`"`s1'"')
        // mean
        if "`s1'"==substr("mean", 1, max(1,`l')) & `"`s2'"'=="" {
            local stats     `stats'   mean
            local statslbl `statslbl' mean
            continue
        }
        // alpha-trimmed mean
        if "`s1'"==substr("alpha", 1, max(1,`l')) {
            if "`s2'"=="" local s2 5                // default is alpha5
            Parse_stats_confirm_num 1 49 `s2' `s'   // range is [1,49]
            local stats    `stats'    alpha:`s2'
            local statslbl `statslbl' alpha`s2'
            continue
        }
        // median
        if "`s1'"==substr("median", 1, max(3,`l')) & `"`s2'"'=="" {
            local stats     `stats'   median
            local statslbl `statslbl' median
            continue
        }
        // Hodges-Lehmann
        if "`s1'"=="hl" & `"`s2'"'=="" {
            local stats    `stats'    hl
            local statslbl `statslbl' HL
            continue
        }
        if "`s1'"=="hlnaive" & `"`s2'"'=="" {
            local stats    `stats'    hlnaive
            local statslbl `statslbl' HL(naive)
            continue
        }
        // Huber M
        if "`s1'"==substr("huber", 1, max(1,`l')) {
            if "`s2'"=="" local s2 95                // default is huber95
            Parse_stats_confirm_num 64 99 `s2' `s'   // range is [64,99]
            local stats    `stats'    huber:`s2'
            local statslbl `statslbl' Huber`s2'
            continue
        }
        // Biweight M
        if "`s1'"==substr("biweight", 1, max(2,`l')) {
            if "`s2'"=="" local s2 95                // default is biweight95
            Parse_stats_confirm_num 1 99 `s2' `s'    // range is [1,99]
            local stats    `stats'    biweight:`s2'
            local statslbl `statslbl' Biweight`s2'
            continue
        }
        // standard deviation
        if "`s1'"=="sd" & `"`s2'"'=="" {
            local stats     `stats'   sd
            local statslbl `statslbl' SD
            continue
        }
        // IQR
        if "`s1'"=="iqr" & `"`s2'"'=="" {
            local stats    `stats'    iqr
            local statslbl `statslbl' IQR
            continue
        }
        // IQRc
        if "`s1'"=="iqrc" & `"`s2'"'=="" {
            local stats    `stats'    iqrc
            local statslbl `statslbl' IQRc
            continue
        }
        // MAD
        if "`s1'"=="mad" & `"`s2'"'=="" {
            local stats    `stats'    mad
            local statslbl `statslbl' MAD
            continue
        }
        // MADN
        if "`s1'"=="madn" & `"`s2'"'=="" {
            local stats    `stats'    madn
            local statslbl `statslbl' MADN
            continue
        }
        // Qn coefficient
        if "`s1'"==substr("qn", 1, max(1,`l')) & `"`s2'"'=="" {
            local stats    `stats'    qn
            local statslbl `statslbl' Qn
            continue
        }
        if "`s1'"=="qnnaive" & `"`s2'"'=="" {
            local stats    `stats'    qnnaive
            local statslbl `statslbl' Qn(naive)
            continue
        }
        // M estimate of scale
        if "`s1'"=="s" {
            if "`s2'"=="" local s2 50               // default is s50
            Parse_stats_confirm_num 1 50 `s2' `s'   // range is [1,50]
            local stats    `stats'    s:`s2'
            local statslbl `statslbl' S`s2'
            continue
        }
        // skewness
        if "`s1'"==substr("skewness", 1, max(3,`l')) & `"`s2'"'=="" {
            local stats    `stats'    skewness
            local statslbl `statslbl' skewness
            continue
        }
        // SK(p)
        if "`s1'"=="sk" {
            if "`s2'"=="" local s2 25               // default is sk25
            Parse_stats_confirm_num 1 49 `s2' `s'   // range is [1,49]
            local stats    `stats'    sk:`s2'
            local statslbl `statslbl' SK`s2'
            continue
        }
        // medcouple
        if "`s1'"=="mc" & `"`s2'"'=="" {
            local stats    `stats'    mc
            local statslbl `statslbl' MC
            continue
        }
        if "`s1'"=="mcnaive" & `"`s2'"'=="" {
            local stats    `stats'    mcnaive
            local statslbl `statslbl' MC(naive)
            continue
        }
        // kurtosis
        if "`s1'"==substr("kurtosis", 1, max(1,`l')) & `"`s2'"'=="" {
            local stats    `stats'    kurtosis
            local statslbl `statslbl' kurtosis
            continue
        }
        // LQW(p)
        if "`s1'"=="qw" {
            if "`s2'"=="" local s2 25               // default is qw25
            Parse_stats_confirm_num 1 49 `s2' `s'   // range is [1,49]
            local stats    `stats'    qw:`s2'
            local statslbl `statslbl' QW`s2'
            continue
        }
        // LQW(p)
        if "`s1'"=="lqw" {
            if "`s2'"=="" local s2 25               // default is lqw25
            Parse_stats_confirm_num 1 49 `s2' `s'   // range is [1,49]
            local stats    `stats'    lqw:`s2'
            local statslbl `statslbl' LQW`s2'
            continue
        }
        // RQW(p)
        if "`s1'"=="rqw" {
            if "`s2'"=="" local s2 25               // default is rqw25
            Parse_stats_confirm_num 1 49 `s2' `s'   // range is [1,49]
            local stats    `stats'    rqw:`s2'
            local statslbl `statslbl' RQW`s2'
            continue
        }
        // LMC
        if "`s1'"=="lmc" & `"`s2'"'=="" {
            local stats    `stats'    lmc
            local statslbl `statslbl' LMC
            continue
        }
        if "`s1'"=="lmcnaive" & `"`s2'"'=="" {
            local stats    `stats'    lmcnaive
            local statslbl `statslbl' LMC(naive)
            continue
        }
        // RMC
        if "`s1'"=="rmc" & `"`s2'"'=="" {
            local stats    `stats'    rmc
            local statslbl `statslbl' RMC
            continue
        }
        if "`s1'"=="rmcnaive" & `"`s2'"'=="" {
            local stats    `stats'    rmcnaive
            local statslbl `statslbl' RMC(naive)
            continue
        }
        di as err `"`s' not allowed in statistics()"'
        exit 198
    }
    if "`stats'"=="" {
        local stats mean
        local statslbl mean
    }
    c_local stats `stats'
    c_local statslbl `statslbl'
end
program Parse_stats_confirm_num
    args l u n s
    capt numlist "`n'", min(1) max(1) range(>=`l' <=`u')
    if _rc {
        di as err `"`s' not allowed in statistics()"'
        exit 198
    }
end

program Estimate_mean, rclass
    syntax if [pw iw fw], v(str) [ u(str) * ]
    if "`weight'"=="pweight" local swgt [aw`exp']
    else if "`weight'"!=""   local swgt [`weight'`exp']
    su `v' `swgt' `if', meanonly
    tempname b
    scalar `b' = r(mean)
    if "`u'"!="" {
        qui replace `u' = (`v' - `b') `if'
    }
    return scalar b = `b'
    return scalar class = 1
end

program Estimate_alpha, rclass sortpreserve
    syntax anything(name=p) if [pw iw fw], v(str) [ u(str) * ]
    if "`weight'"=="pweight" local swgt  [aw`exp']
    else if "`weight'"!=""   local swgt  [`weight'`exp']
    local p = substr("`p'", 2, .) // strip leading colon
    marksample touse
    tempvar W
    if `"`exp'"'!="" qui gen double `W' `exp' if `touse'
    else             qui gen double `W' = 1 if `touse'
    sort `touse' `v' `W'
    qui by `touse': replace `W' = sum(`W') if `touse'
    qui by `touse': replace `W' = `W'*100/`W'[_N] if `touse'
    tempvar select
    qui by `touse': gen byte `select' = ///
        -1 + (`W'>`p') + (`W'[_n-1]>=(100-`p') & _n>1) if `touse'
    su `v' `swgt' if `touse' & `select'==0, meanonly
    tempname b
    scalar `b' = r(mean)
    if "`u'"!="" {
        tempname qlo qup plo pup
        scalar `qlo' = r(min)
        scalar `qup' = r(max)
        tempvar zlo zup
        qui gen byte `zlo' = (`select'==-1) if `touse'
        su `zlo' `swgt' if `touse', meanonly
        scalar `plo' = r(mean)
        qui gen byte `zup' = (`select'==1) if `touse'
        su `zup' `swgt' if `touse', meanonly
        scalar `pup' = r(mean)
        qui replace `u' = (`v'*(`select'==0) + `qlo'*(`zlo'-`plo') ///
                      + `qup'*(`zup'-`pup')) / (1-`plo'-`pup') - `b' if `touse'
    }
    return scalar b = `b'
    return scalar class = 1
end

program Estimate_median, rclass
    syntax if [pw iw fw/], v(str) [ u(str) ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    marksample touse
    _pctile `v' `pcwgt' if `touse', percentiles(50)
    tempname b
    scalar `b' = r(r1)
    if "`u'"!="" {
        if "`weight'"=="pweight" local swgt  [aw = `exp']
        else if "`weight'"!=""   local swgt  [`weight' = `exp']
        tempname d
        mata: robstat_kdens_s("`b'", "`d'")
        tempvar z
        qui generate byte `z' = (`v'<=`b') if `touse'
        summarize `z' `swgt' if `touse', meanonly
        qui replace `u' = (r(mean) - `z') / `d' if `touse'
    }
    return scalar b = `b'
    return scalar class = 1
end

program Estimate_hl, rclass
    syntax if [pw iw fw/], v(str) [ u(str) naive ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="pweight" local swgt  [aw = `exp']
    else if "`weight'"!=""   local swgt  [`weight' = `exp']
    marksample touse
    tempname b
    mata: robstat_estimate_pairwise("hl") // sets scalar b
    if "`u'"!="" {
        tempvar v1 F1 dvar
        qui generate double `v1'   = 2*`b' - `v' if `touse'
        qui generate double `F1'   = .
        qui generate double `dvar' = .
        mata: robstat_relrank("`v1'", "`F1'")
        tempname mF1
        summarize `F1' `swgt' if `touse', meanonly
        scalar `mF1' = r(mean)
        mata: robstat_kdens_v("`v1'", "`dvar'")
        summarize `dvar' `swgt' if `touse', meanonly
        qui replace `u' = (`mF1' - `F1') / r(mean) if `touse'
    }
    return scalar b = `b'
    return scalar class = 1
end
program Estimate_hlnaive
    Estimate_hl `0' naive
end

program Estimate_huber, rclass
    syntax anything(name=p) if [pw iw fw/], v(str) [ u(str) ///
        tolerance(str) iterate(str) ///
        kernel(str) bw(str) adaptive(str) n(str) ]
    local p = substr("`p'", 2, .) // strip leading colon
    marksample touse
    tempname b k
    mata: robstat_estimate_m(`p', "huber", `tolerance', `iterate')
        // returns result in b and k, fills in u
    return scalar b = `b'
    return scalar k = `k'
    return scalar class = 1
end

program Estimate_biweight, rclass
    syntax anything(name=p) if [pw iw fw/], v(str) [ u(str) ///
        tolerance(str) iterate(str) ///
        kernel(str) bw(str) adaptive(str) n(str) ]
    local p = substr("`p'", 2, .) // strip leading colon
    marksample touse
    tempname b k
    mata: robstat_estimate_m(`p', "biweight", `tolerance', `iterate')
        // returns result in b and k, fills in u
    return scalar b = `b'
    return scalar k = `k'
    return scalar class = 1
end

program Estimate_sd, rclass
    syntax if [pw iw fw], v(str) [ u(str) * ]
    if "`weight'"=="pweight" local swgt  [aw`exp']
    else if "`weight'"!=""   local swgt  [`weight'`exp']
    qui su `v' `swgt' `if'
    tempname mean sd
    scalar `mean' = r(mean)
    scalar `sd' = r(sd)
    if (`sd'>=.) scalar `sd' = 0  // can happen if n=1
    if "`u'"!="" {
        if (`sd'!=0) {
            tempname c
            if "`weight'"=="iweight" scalar `c' = r(sum_w) / (r(sum_w) - 1)
            else                     scalar `c' = r(N) / (r(N) - 1)
            qui replace `u' = 1/(2*`sd') * (`c'*(`v'-`mean')^2 - `sd'^2) `if'
        }
    }
    return scalar b = `sd'
    return scalar class = 2
end

program Estimate_iqr, rclass
    syntax if [pw iw fw/], v(str) [ u(str) corr ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    marksample touse
    tempname q1 q3 b d
    scalar `d' = cond("`corr'"!="", 1/(invnormal(0.75)-invnormal(0.25)), 1)
    _pctile `v' `pcwgt' if `touse', percentiles(25 75)
    scalar `q1' = r(r1)
    scalar `q3' = r(r2)
    scalar `b' = (`q3' - `q1') * `d'
    if "`u'"!="" {
        if "`weight'"=="pweight" local swgt  [aw = `exp']
        else if "`weight'"!=""   local swgt  [`weight' = `exp']
        tempvar  z1 z3
        tempname p1 p3
        qui generate byte `z1' = (`v'<`q1') if `touse'
        summarize `z1' `swgt' if `touse', meanonly
        scalar `p1' = r(mean)
        qui generate byte `z3' = (`v'<`q3') if `touse'
        summarize `z3' `swgt' if `touse', meanonly
        scalar `p3' = r(mean)
        tempname d1 d3
        mata: robstat_kdens_s(("`q1'", "`q3'"), ("`d1'", "`d3'"))
        qui replace `u' = `d' * ((`p3'-`z3')/`d3' - (`p1'-`z1')/`d1') if `touse'
    }
    return scalar b = `b'
    return scalar class = 2
end
program Estimate_iqrc
    Estimate_iqr `0' corr
end

program Estimate_mad, rclass
    syntax if [pw iw fw/], v(str) [ u(str) corr ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    marksample touse
    tempname b med mad d
    scalar `d' = cond("`corr'"!="", 1/invnormal(0.75), 1)
    _pctile `v' `pcwgt' if `touse', percentiles(50)
    scalar `med' = r(r1)
    tempvar tmp
    qui generate double `tmp' = abs(`v' - `med') if `touse'
    _pctile `tmp' `pcwgt' if `touse', percentiles(50)
    scalar `mad' = r(r1)
    scalar `b' = `d' * `mad'
    if "`u'"!="" {
        if "`weight'"=="pweight" local swgt [aw = `exp']
        else if "`weight'"!=""   local swgt [`weight' = `exp']
        tempvar  z1 z2
        tempname p1 p2
        qui generate byte `z1' = (`tmp' <= (`b'/`d'))  if `touse'
        summarize `z1' `swgt' if `touse', meanonly
        scalar `p1' = r(mean)
        qui generate byte `z2' = (`v' <= `med')  if `touse'
        summarize `z2' `swgt' if `touse', meanonly
        scalar `p2' = r(mean)
        tempname q1 q3 d1 d2 d3
        scalar `q1' = `med' - `mad'
        scalar `q3' = `med' + `mad'
        mata: robstat_kdens_s(("`q1'", "`med'", "`q3'"), ///
                              ("`d1'", "`d2'", "`d3'"))
        qui replace `u' = `d' * ((`p1'-`z1')/(`d1'+`d3') + ///
            (`p2'-`z2')*((`d1'-`d3')/(`d2'*(`d1'+`d3')))) if `touse'
    }
    return scalar b = `b'
    return scalar class = 2
end
program Estimate_madn
    Estimate_mad `0' corr
end

program Estimate_qn, rclass
    syntax if [pw iw fw/], v(str) [ u(str) naive ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="pweight" local swgt  [aw = `exp']
    else if "`weight'"!=""   local swgt  [`weight' = `exp']
    marksample touse
    tempname b d
    scalar `d' = 1/(sqrt(2) * invnormal(5/8))
    mata: robstat_estimate_pairwise("qn") // sets scalar b
    if "`u'"!="" {
        tempvar v1 v2 F1 F2 dvar
        qui generate double `v1'   = `v' + `b'/`d' if `touse'
        qui generate double `v2'   = `v' - `b'/`d' if `touse'
        qui generate double `F1'   = .
        qui generate double `F2'   = .
        qui generate double `dvar' = .
        mata: robstat_relrank(("`v1'", "`v2'"), ("`F1'", "`F2'"))
        tempname p1 p2
        summarize `F1' `swgt' if `touse', meanonly
        scalar `p1' = r(mean)
        summarize `F2' `swgt' if `touse', meanonly
        scalar `p2' = r(mean)
        mata: robstat_kdens_v("`v1'", "`dvar'")
        summarize `dvar' `swgt' if `touse', meanonly
        qui replace `u' = `d' * ((`p1'-`p2') - (`F1'-`F2')) / r(mean) if `touse'
    }
    return scalar b = `b'
    return scalar class = 2
end
program Estimate_qnnaive
    Estimate_qn `0' naive
end
//     if (n<=9) {
//         d = d * (0.399, 0.994, 0.512, 0.844, 0.611, 0.857, 0.669, 0.872)[n-1]
//     }
//     else {
//         if (mod(n,2)) d = d * n/(n + 1/4)
//         else          d = d * n/(n + 3.8)
//     }

program Estimate_s, rclass
    syntax anything(name=p) if [pw iw fw/], v(str) [ u(str) ///
        tolerance(str) iterate(str) ///
        kernel(str) bw(str) adaptive(str) n(str) ]
    local p = substr("`p'", 2, .) // strip leading colon
    marksample touse
    tempname b k
    mata: robstat_estimate_m(`p', "scale", `tolerance', `iterate')
        // returns result in b and k, fills in u
    return scalar b = `b'
    return scalar k = `k'
    return scalar class = 2
end

program Estimate_skewness, rclass
    syntax if [pw iw fw], v(str) [ u(str) * ]
    if inlist("`weight'","pweight","iweight") local swgt  [aw`exp']
    else if "`weight'"!=""                    local swgt  [`weight'`exp']
    qui su `v' `swgt' `if', detail
    tempname b
    scalar `b' = r(skewness)
    if (`b'>=.) scalar `b' = 0  // can happen if n=1
    else if "`u'"!="" {
        tempvar z
        qui generate double `z' = (`v'-r(mean))/r(sd) `if'
        qui replace `u' = `z'^3 - 3*`z' - `b' - `b'*(3/2)*(`z'^2-1) `if'
    }
    return scalar b = `b'
    return scalar class = 3
end

program Estimate_sk, rclass
    syntax anything(name=p1) if [pw iw fw/], v(str) [ u(str) ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    local p1 = substr("`p1'", 2, .) // strip leading colon
    local p3 = 100 - `p1'
    marksample touse
    _pctile `v' `pcwgt' if `touse', percentiles(`p1' 50 `p3')
    tempname q1 q2 q3 b
    scalar `q1' = r(r1)
    scalar `q2' = r(r2)
    scalar `q3' = r(r3)
    scalar `b' = (`q1' + `q3' - 2*`q2') / (`q3' - `q1')
    if (`b'>=.) scalar `b' = 0  // can happen if n=1
    else if "`u'"!="" {
        if "`weight'"=="pweight" local swgt [aw = `exp']
        else if "`weight'"!=""   local swgt [`weight' = `exp']
        forv i=1/3 {
            tempvar  z`i'
            tempname P`i'
            qui gen byte `z`i'' = (`v' <= `q`i'')  if `touse'
            summarize `z`i'' `swgt' if `touse', meanonly
            scalar `P`i'' = r(mean)
        }
        tempname d1 d2 d3
        mata: robstat_kdens_s(("`q1'", "`q2'", "`q3'"), ("`d1'", "`d2'", "`d3'"))
        qui replace `u' = 2 * ///
            ((`q3'-`q2')*((`P1'-`z1')/`d1' - (`P2'-`z2')/`d2')  ///
           - (`q2'-`q1')*((`P2'-`z2')/`d2' - (`P3'-`z3')/`d3')) ///
            / (`q3' - `q1')^2 if `touse'
    }
    return scalar b = `b'
    return scalar class = 3
end

program Estimate_mc, rclass
    syntax if [pw iw fw/], V0(str) [ u(str) naive ///
        kernel(str) bw(str) adaptive(str) n(str) lmc rmc * ]
    if "`weight'"=="pweight" local swgt  [aw = `exp']
    else if "`weight'"!=""   local swgt  [`weight' = `exp']
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    marksample touse
    tempname v b q med
    if "`lmc'"!="" {
        _pctile `v0' `pcwgt' if `touse', percentiles(50)
        scalar `med' = r(r1)
        local touse0 `touse'
        tempvar touse
        qui gen byte `touse' = `touse0' & (`v0'<`med')
        qui count if `touse'
        if r(N)<1 {
            return scalar b = 0
            exit
        }
        _pctile `v0' `pcwgt' if `touse', percentiles(50)
        scalar `q'   = r(r1)
        qui generate double `v' = (`q' - `v0') if `touse'
        mata: robstat_estimate_pairwise("mc") // sets scalar b
        local touse `touse0'
    }
    else if "`rmc'"!="" {
        _pctile `v0' `pcwgt' if `touse', percentiles(50)
        scalar `med' = r(r1)
        local touse0 `touse'
        tempvar touse
        qui gen byte `touse' = `touse0' & (`v0'>`med')
        qui count if `touse'
        if r(N)<1 {
            return scalar b = 0
            exit
        }
        _pctile `v0' `pcwgt' if `touse', percentiles(50)
        scalar `q'   = r(r1)
        qui generate double `v' = `v0' - `q' if `touse'
        mata: robstat_estimate_pairwise("mc") // sets scalar b
        local touse `touse0'
    }
    else {
        _pctile `v0' `pcwgt' if `touse', percentiles(50)
        scalar `q' = r(r1)
        qui generate double `v' = `v0' - `q' if `touse'
        mata: robstat_estimate_pairwise("mc") // sets scalar b
    }
    if "`u'"!="" {
        if (abs(`b')==1) {
            di as txt "(mc = " `b' ": cannot compute standard errors)"
            return scalar b = `b'
            exit
        }
        tempvar g1 g2 F1 F2 dg
        tempname mc Idg dq dH gmed Fmed
        if "`lmc'"!="" {
            scalar `mc' = -`b'
            local v `v0'
        }
        else if "`rmc'"!="" {
            scalar `mc' = -`b'
            scalar `q' = -`q'
            scalar `med' = -`med'
            drop `v'
            qui generate `v' = -`v0' if `touse'
        }
        else {
            scalar `mc' = `b'
            scalar `med' = .
            local v `v0'
        }
        qui generate `g1'   = (`v'*(`mc'-1) + 2*`q') / (`mc'+1) if `touse'
        qui generate `g2'   = (`v'*(`mc'+1) - 2*`q') / (`mc'-1) if `touse'
        qui generate `F1'   = .
        qui generate `F2'   = .
        mata: robstat_relrank(("`g1'", "`g2'"), ("`F1'", "`F2'"))
        qui generate `dg' = .
        mata: robstat_kdens_v("`g1'", "`dg'")
        qui replace `dg' = `dg' * (`v'>=`q' & `v'<=`med') if `touse'
        summarize `dg' `swgt' if `touse', meanonly
        scalar `Idg' = r(mean)
        mata: robstat_kdens_s("`q'", "`dq'")
        if "`lmc'`rmc'"!="" {
            qui replace `g1' = 32 * `dg' * ((`q'-`v')/(`mc'+1)^2) if `touse'
        }
        else {
            qui replace `g1' = 8 * `dg' * ((`v'-`q')/(`mc'+1)^2) if `touse'
        }
        summarize `g1' `swgt' if `touse', meanonly
        scalar `dH' = r(mean)
        if (`dH'==0) {
            //di as txt "(mc: cannot estimate standard errors)"
        }
        else if "`lmc'`rmc'"!="" {
            scalar `gmed' = (`med'*(`mc'-1) + 2*`q') / (`mc'+1)
            mata: robstat_relrank_s("`gmed'", "`Fmed'")
            qui replace `u' = 1/`dH' * (1 ///
                - 16*`F1'*((`v'>`q' & `v'<`med')) ///
                - 16*(`F2'-.25)*(`v'>`gmed' & `v'<`q') ///
                - 4*(`v'<`gmed') ///
                - 8*sign(`v'-`med')*`Fmed' ///
                + (.25-(`v'<`q'))*(4 - 32*`Idg'/(`dq'*(`mc'+1)))) if `touse'
        }
        else {
            qui replace `u' = 1/`dH' * ///
                (1 - 4*`F1'*(`v'>`q') - 4*(`F2'-.5)*(`v'<`q') ///
                + sign(`v'-`q')*(1 - 4*`Idg'/(`dq'*(`mc'+1)))) if `touse'
        }
        capt assert (`u'<.) if `touse'
        if _rc {
            di as txt "(mc: IF = .; cannot compute standard errors)"
            qui replace `u' = 0 if `touse'
        }
        else {
            // make sure that IF is centered at zero
            summarize `u' `swgt' if `touse', meanonly
            qui replace `u' = `u' - r(mean) if `touse'
        }
    }
    return scalar b = `b'
    return scalar class = 3
end
program Estimate_mcnaive
    Estimate_mc `0' naive
end

program Estimate_kurtosis, rclass
    syntax if [pw iw fw], v(str) [ u(str) * ]
    if inlist("`weight'","pweight","iweight") local swgt  [aw`exp']
    else if "`weight'"!=""                    local swgt  [`weight'`exp']
    qui su `v' `swgt' `if', detail
    tempname b
    scalar `b' = r(kurtosis)
    if (`b'>=.) scalar `b' = 0  // can happen if n=1
    else if "`u'"!="" {
        tempvar z
        qui generate double `z' = (`v'-r(mean))/r(sd) `if'
        qui replace `u' = (`z'^2 - `b')^2 - `b'*(`b'-1) - 4*r(skewness)*`z' `if'
    }
    return scalar b = `b'
    return scalar class = 4
end

program Estimate_qw, rclass
    syntax anything(name=p) if [pw iw fw/], v(str) [ u(str) ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    local p = substr("`p'", 2, .) // strip leading colon
    local p1 = `p'/2
    local p2 = `p'          // or should this be 25?
    local p3 = 50 - `p'/2
    local p4 = 50 + `p'/2
    local p5 = 100 - `p'    // or should this be 75?
    local p6 = 100 - `p'/2
    marksample touse
    _pctile `v' `pcwgt' if `touse', percentiles(`p1' `p2' `p3' `p4' `p5' `p6')
    tempname q1 q2 q3 q4 q5 q6 b
    scalar `q1' = r(r1)
    scalar `q2' = r(r2)
    scalar `q3' = r(r3)
    scalar `q4' = r(r4)
    scalar `q5' = r(r5)
    scalar `q6' = r(r6)
    scalar `b' = (`q6' - `q4' + `q3' - `q1')/ (`q5' - `q2')
    if (`q1'==`q3') scalar `b' = 0
    else if "`u'"!="" {
        if "`weight'"=="pweight" local swgt [aw = `exp']
        else if "`weight'"!=""   local swgt [`weight' = `exp']
        forv i=1/6 {
            tempvar  z`i'
            tempname P`i'
            qui gen byte `z`i'' = (`v' <= `q`i'')  if `touse'
            summarize `z`i'' `swgt' if `touse', meanonly
            scalar `P`i'' = r(mean)
        }
        tempname d1 d2 d3 d4 d5 d6
        mata: robstat_kdens_s(  ///
            ("`q1'", "`q2'", "`q3'", "`q4'", "`q5'", "`q6'"), ///
            ("`d1'", "`d2'", "`d3'", "`d4'", "`d5'", "`d6'"))
        qui replace `u' = ///
            ((`q5'-`q2')*((`P6'-`z6')/`d6' - (`P4'-`z4')/`d4'  ///
                        + (`P3'-`z3')/`d3' - (`P1'-`z1')/`d1') ///
            - (`q6' - `q4' + `q3' - `q1') * ///
              ((`P5'-`z5')/`d5' - (`P2'-`z2')/`d2')) ///
            / (`q5' - `q2')^2 if `touse'
    }
    return scalar b = `b'
    return scalar class = 4
end

program Estimate_lqw, rclass
    syntax anything(name=p) if [pw iw fw/], v(str) [ u(str) ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    local p = substr("`p'", 2, .) // strip leading colon
    local p1 = `p'/2
    local p2 = 25
    local p3 = 50 - `p'/2
    marksample touse
    _pctile `v' `pcwgt' if `touse', percentiles(`p1' `p2' `p3')
    tempname q1 q2 q3 b
    scalar `q1' = r(r1)
    scalar `q2' = r(r2)
    scalar `q3' = r(r3)
    scalar `b' = - (`q1' + `q3' - 2*`q2') / (`q3' - `q1')
    if (`q1'==`q3') scalar `b' = 0
    else if "`u'"!="" {
        if "`weight'"=="pweight" local swgt [aw = `exp']
        else if "`weight'"!=""   local swgt [`weight' = `exp']
        forv i=1/3 {
            tempvar  z`i'
            tempname P`i'
            qui gen byte `z`i'' = (`v' <= `q`i'')  if `touse'
            summarize `z`i'' `swgt' if `touse', meanonly
            scalar `P`i'' = r(mean)
        }
        tempname d1 d2 d3
        mata: robstat_kdens_s(("`q1'", "`q2'", "`q3'"), ("`d1'", "`d2'", "`d3'"))
        qui replace `u' = 2 * ///
            ((`q2'-`q1')*((`P2'-`z2')/`d2' - (`P3'-`z3')/`d3')  ///
           - (`q3'-`q2')*((`P1'-`z1')/`d1' - (`P2'-`z2')/`d2')) ///
            / (`q3' - `q1')^2 if `touse'
    }
    return scalar b = `b'
    return scalar class = 4
end

program Estimate_rqw, rclass
    syntax anything(name=p) if [pw iw fw/], v(str) [ u(str) ///
        kernel(str) bw(str) adaptive(str) n(str) * ]
    if "`weight'"=="iweight" local pcwgt [aw = `exp']
    else if "`weight'"!=""   local pcwgt [`weight' = `exp']
    local p = substr("`p'", 2, .) // strip leading colon
    local p1 = 50 + `p'/2
    local p2 = 75
    local p3 = 100 - `p'/2
    marksample touse
    _pctile `v' `pcwgt' if `touse', percentiles(`p1' `p2' `p3')
    tempname q1 q2 q3 b
    scalar `q1' = r(r1)
    scalar `q2' = r(r2)
    scalar `q3' = r(r3)
    scalar `b' = (`q1' + `q3' - 2*`q2') / (`q3' - `q1')
    if (`q1'==`q3') scalar `b' = 0
    else if "`u'"!="" {
        if "`weight'"=="pweight" local swgt [aw = `exp']
        else if "`weight'"!=""   local swgt [`weight' = `exp']
        forv i=1/3 {
            tempvar  z`i'
            tempname P`i'
            qui gen byte `z`i'' = (`v' <= `q`i'')  if `touse'
            summarize `z`i'' `swgt' if `touse', meanonly
            scalar `P`i'' = r(mean)
        }
        tempname d1 d2 d3
        mata: robstat_kdens_s(("`q1'", "`q2'", "`q3'"), ("`d1'", "`d2'", "`d3'"))
        qui replace `u' = 2 * ///
            ((`q3'-`q2')*((`P1'-`z1')/`d1' - (`P2'-`z2')/`d2')  ///
           - (`q2'-`q1')*((`P2'-`z2')/`d2' - (`P3'-`z3')/`d3')) ///
            / (`q3' - `q1')^2 if `touse'
    }
    return scalar b = `b'
    return scalar class = 4
end

program Estimate_lmc
    Estimate_mc `0' lmc
end
program Estimate_lmcnaive
    Estimate_mc `0' naive lmc
end

program Estimate_rmc
    Estimate_mc `0' rmc
end
program Estimate_rmcnaive
    Estimate_mc `0' naive rmc
end

version 11
mata mata set matastrict on
mata:

void robstat_cilog()
{
    real scalar      i, n
    string scalar    note
    real rowvector   b, se, crit, ll, ul, type
    real matrix      table
    string vector    cstripe, rstripe

    table   = st_matrix("r(table)")
    cstripe = st_matrixcolstripe("r(table)")[,2]
    rstripe = st_matrixrowstripe("r(table)")[,2]
    type    = st_matrix("e(class)")
    b = se = ll = ul = crit = J(1, cols(table), .)
    n = rows(rstripe)
    for (i=1;i<=n;i++) {
        if      (rstripe[i]=="b")    b    = table[i,]
        else if (rstripe[i]=="se")   se   = table[i,]
        else if (rstripe[i]=="ll")   ll   = table[i,]
        else if (rstripe[i]=="ul")   ul   = table[i,]
        else if (rstripe[i]=="crit") crit = table[i,]
    }
    n = rows(cstripe)
    for (i=1;i<=n;i++) {
        if (type[i]==2) { // Scale statistic
            se[i] = se[i] :/ b[i]
            b[i]  = ln(b[i])
            ll[i] = exp(b[i] - crit[i] * se[i])
            ul[i] = exp(b[i] + crit[i] * se[i])
        }
    }
    st_matrix(st_local("CI"), ll \ ul)
    note = "log-transformed confidence interval"
    if (allof(type, 2)) {
        if (length(type)>1) note = note + "s"
    }
    else note = note + "s for scale statistics"
    st_local("cilognote", note)
}

// helper program to rearrange b and V
void _robstat_swap_eq_and_coefs(real scalar keq, real scalar kcoef)
{
    real scalar    i, j, k
    real colvector p
    string matrix  mstripe

    // anything to do?
    if (keq==1 & (st_local("over")=="" | st_local("nstats")=="1")) return
    // permutation vector
    p = J(keq*kcoef, 1, .)
    i = 0
    for (j=1;j<=kcoef;j++) {
        for (k=1;k<=keq;k++) {
            p[++i] = j + kcoef*(k-1)
        }
    }
    st_local("k_eq", st_local("k_coef"))
    st_local("k_coef", strofreal(keq))
    // rearrange labels
    mstripe = st_matrixcolstripe(st_local("b"))[p,(2,1)]
    // rearrange b
    st_replacematrix(st_local("b"), st_matrix(st_local("b"))[1,p])
    st_matrixcolstripe(st_local("b"), mstripe)
    // rearrange V
    if (st_local("V")!="") {
        st_replacematrix(st_local("V"), st_matrix(st_local("V"))[p,p])
        st_matrixcolstripe(st_local("V"), mstripe)
        st_matrixrowstripe(st_local("V"), mstripe)
    }
}

// M-estimation
void robstat_estimate_m(
    real scalar p,
    string scalar obj,
    real scalar tol,
    real scalar iter)
{
    real scalar    b, k, med, pw
    real colvector x, w, u
    transmorphic   S

    // data
    x = st_data(., st_local("v"), st_local("touse"))
    if (st_local("weight")!="") {
        w = st_data(., st_local("exp"), st_local("touse"))
        if (st_local("weight")!="fweight") w = w :/ quadsum(w) * rows(w) // normalize weights
    }
    else  w = 1
    // estimation
    med = mm_median(x, w)
    if (obj=="scale") {
        S = mm_mscale(x, w, p, ., med, 0, tol, iter)
        b = mm_mscale_b(S)
        k = mm_mscale_k(S)
    }
    else {
        S = mm_mloc(x, w, p, obj, med, ., 0, tol, iter)
        b = mm_mloc_b(S)
        k = mm_mloc_k(S)
    }
    st_numscalar(st_local("b"), b)
    st_numscalar(st_local("k"), k)
    // IFs
    if (st_local("u")=="") return
    pw = (st_local("weight")=="pweight")
    if (obj=="scale")
         u = _robstat_mscale_IF(x, w, pw, k, b, med, mm_mscale_delta(S))
    else u = _robstat_mloc_IF(x, w, pw, obj, k, b, mm_mloc_s(S), med)
    st_store(., st_local("u"), st_local("touse"), u)
}

real colvector _robstat_mloc_IF(
    real colvector x,
    real colvector w,
    real scalar    pw,
    string scalar  obj,
    real scalar    k,
    real scalar    b,
    real scalar    s,
    real scalar    med)
{
    real scalar    mad
    real colvector d, z, phi, psi, zmed, zmad

    mad = s * invnormal(0.75)
    d   = robstat_kdens(x, w, pw, (med-mad, med, med+mad)')
    z   = (x :- b) / s
    if (obj=="biweight") {
        phi = mm_biweight_phi(z, k)
        psi = mm_biweight_psi(z, k)
    }
    else {
        phi = mm_huber_phi(z, k)
        psi = mm_huber_psi(z, k)
    }
    zmad = (abs(x :- med) :<= mad)
    zmed = (x :<= med)
    return((s*psi - (mean(phi:*z, w) / ((d[1]+d[3]) * invnormal(0.75))) *
            ((mean(zmad, w) :- zmad) +
            ((d[1]-d[3])/d[2]) * (mean(zmed, w) :- zmed))) / mean(phi, w))
}

real colvector _robstat_mscale_IF(
    real colvector x,
    real colvector w,
    real scalar    pw,
    real scalar    k,
    real scalar    s,
    real scalar    med,
    real scalar    delta)
{
    real colvector d, z, psi, zmed

    d    = robstat_kdens(x, w, pw, med)
    z    = (x :- med) / s
    psi  = mm_biweight_psi(z, k)
    zmed = (x :<= med)
    return((s * (mm_biweight_rho(z, k) :- delta)
        - mean(psi, w)/d * (mean(zmed, w) :- zmed))
        / mean(psi:*z, w))
}

// Pairwise estimators
void robstat_estimate_pairwise(string scalar s)
{
    real scalar    fw, naive
    real colvector x, w

    fw    = (st_local("weight")=="fweight")
    naive = (st_local("naive")!="")
    x     = st_data(., st_local("v"), st_local("touse"))
    if (st_local("weight")=="") w = 1
    else w = st_data(., st_local("exp"), st_local("touse"))
    if      (s=="hl") st_numscalar(st_local("b"), mm_hl(x, w, fw, naive))
    else if (s=="qn") st_numscalar(st_local("b"), mm_qn(x, w, fw, naive))
    else if (s=="mc") st_numscalar(st_local("b"), mm_mc(x, w, fw, naive))
}

// density estimations
void robstat_kdens_v(string scalar atvar, string scalar dvar)
{   // store density estimate in variable
    real scalar    pw
    string scalar  touse, wtype
    real colvector x, w, at

    // data
    touse = st_local("touse")
    x  = st_data(., st_local("v"), touse)
    wtype = st_local("weight")
    if (wtype!="") {
        w = st_data(., st_local("exp"), touse)
        if (wtype!="fweight") w = w :/ quadsum(w) * rows(w) // normalize weights
    }
    else w = 1
    pw = wtype=="pweight"
    // density estimate
    at = st_data(., atvar, touse)
    st_store(., dvar, touse, robstat_kdens(x, w, pw, at))
}

void robstat_kdens_s(string rowvector in, string rowvector out)
{   // store density estimate in scalars
    real scalar    i, pw
    string scalar  touse, wtype
    real colvector x, w, at, d

    // data
    touse = st_local("touse")
    x  = st_data(., st_local("v"), touse)
    wtype = st_local("weight")
    if (wtype!="") {
        w = st_data(., st_local("exp"), touse)
        if (wtype!="fweight") w = w :/ quadsum(w) * rows(w) // normalize weights
    }
    else w = 1
    pw = wtype=="pweight"
    // density estimate
    at = J(cols(in), 1, .)
    for (i=1; i<=cols(in); i++) at[i] = st_numscalar(in[i])
    d = robstat_kdens(x, w, pw, at)
    for (i=1; i<=cols(out); i++) st_numscalar(out[i], d[i])
}

real colvector robstat_kdens(
    real colvector x,
    real colvector w,
    real scalar    pw,
    real colvector at)
{
    class mm_density scalar D

    if (mm_isconstant(x)) return(x[1]:==at)
    D.data(x, w, pw)
    D.kernel(st_local("kernel"), strtoreal(st_local("adaptive")))
    D.bw(st_local("bw"), 1, 2, 1)
    D.n(strtoreal(st_local("n")))
    return(D.d(at))
}

// relative ranks
void robstat_relrank(string rowvector in, string rowvector out)
{
    string scalar  touse
    real colvector x, w

    touse = st_local("touse")
    x = st_data(., st_local("v"), touse)
    w = (st_local("weight")!="" ? st_data(., st_local("exp"), touse) : 1)
    st_store(., out, touse,
        mm_relrank(x, w, st_data(., in, touse), "midpoints"!=""))
}

void robstat_relrank_s(string rowvector in, string rowvector out)
{
    string scalar  touse
    real scalar    i
    real colvector x, w, at, r

    touse = st_local("touse")
    x = st_data(., st_local("v"), touse)
    w = (st_local("weight")!="" ? st_data(., st_local("exp"), touse) : 1)
    at = J(cols(in), 1, .)
    for (i=1; i<=cols(in); i++) at[i] = st_numscalar(in[i])
    r = mm_relrank(x, w, at, "midpoints"!="")
    for (i=1; i<=cols(out); i++) st_numscalar(out[i], r[i])
}

end

exit