*! mpi
*! version 1.0.3 16jun2017
*! Author: Daniele Pacifico, daniele.pacifico@oecd.org
*! Author: Felix Pöge, felix.poege@ip.mpg.de

/* ---------------------------------------------------------------------------

Computing Main Indicators

    - If the by-option is specified, ONLY the by-version is returned.
    - Do not call independently.

--------------------------------------------------------------------------- */

program define _mpi_main, sortpreserve rclass

    syntax anything(equalok)

    quietly {

        // Read in all the to-be-used locals supplied to the program.
        foreach pair in `*' {
            local eq_pos = strpos("`pair'", "=")
            local before = substr("`pair'", 1, `eq_pos'-1)
            local after = substr("`pair'", `eq_pos'+1, .)
            local `before' = "`after'"
        }

        tempname res resV res_add resV_add

        /*
            If a by-option was specified, prepare the respective statements.
        */
        if "`by'" != "" {
            local over_statement = ", over(`by', nolabel)"
        }
        else {
            local over_statement = ""
        }

        /*
            Check if all individuals are non-poor, as this requires additional
            messages below.
        */
        cap assert `_H' == 0
        if _rc == 9 {
            local warn_no_poor = 0
        }
        else {
            error _rc
            local warn_no_poor = 1
        }

        local main_eval = "`_H' `_M0'"
        // "Average Deprivation Share Among Poor"
        local add_eval = "(`_M0' / `_H')"

        // If possible, compute poverty gaps
        if `use_thresholds' {
            // `gap_1': "Adjusted Multidimensional Deprivation Gap"
            // G = |g1|/|g0|
            // `gap_2' "Adjusted Squared Multidimensional Deprivation Gap"

            local add_eval = "`add_eval' (`gap_1' / `_M0') (`gap_2' / `_M0')"

            // M1 = H*A*G = M0*G (gap 1)
            // M2 = H*A*S(2) = M0*S(2) (gap 2)
            local main_eval = "`main_eval' `gap_1' `gap_2'"

            // Computing M_alpha
            foreach a of local alpha {
                if inlist(`a', 1, 2) {
                    continue
                }
                local a_name = subinstr("`a'", ".", "p", .)
                local a_name = subinstr("`a_name'", "-", "m", .)

                // `gap_`a_name'' "Adjusted Multidimensional Deprivation Gap with transformation alpha"
                local main_eval = "`main_eval' `gap_`a_name''"
                local add_eval = "`add_eval' (`gap_`a_name'' / `_M0')"
            } // end alpha iteration
        } // end use_thresholds

        /*
            Perform the evaluation
        */
        `svy' mean `main_eval' `weight_exp' `l_if' `over_statement'
        matrix `res' = e(b)
        matrix `resV' = e(V)
        local N = e(N)

        // The additional variables (A, G, S, ...) will not be by-decomposed.
        if "`by'" == "" {
            `svy' ratio `add_eval' `weight_exp' `l_if'
            matrix `res_add' = e(b)
            matrix `resV_add' = e(V)

            /*
                There (used to?) exist a bug in Stata where a ratio of (x / y) with y=0 returns invalid
                values in e(V). This code here will detect the invalid entries and replace them
                with zero.
                Then, the behavior is equivalent to a zero-variation entry.
                See:
                http://www.statalist.org/forums/forum/general-stata-discussion/general/1348095-create-or-test-for-1-ind-in-a-matrix
            */
            local cols = colsof(`resV_add')
            forvalues i = 1/`cols' {
                forvalues j = 1/`cols' {
                    // The entry is not missing if unmodified, but if
                    // used in a formula, it becomes missing.
                    if missing(1 * `resV_add'[`i', `j']) & !missing(`resV_add'[`i', `j']) {
                        matrix `resV_add'[`i', `j'] = 0
                    }
                }
            }

            // Checking the computation: M0 = H*A
            local H = `res'[1, 1]
            local M0 = `res'[1, 2]
            local A = `res_add'[1, 1]
            capture assert abs(`M0' - `H'*`A') < `assertion_sensitivity' | ///
                           (missing(`M0') | missing(`H'*`A'))
            if _rc == 9 {
                noi di as error "Computation M0 failed. M0=`M0' H=`H' A=`A'."
            }
            else {
                error _rc
            }

            if `use_thresholds' {
                local G = `res_add'[1, 2]
                local S = `res_add'[1, 3]
                local M1 = `res'[1, 3]
                local M2 = `res'[1, 4]
                local alt_M1 = `M0'*`G'
                cap assert abs(`M1' - `alt_M1') < `assertion_sensitivity' | ///
                           (missing(`alt_M1') | missing(`M1'))
                if _rc != 0 {
                    noi display as error "Computation M1 failed (`alt_M1' vs. " `M1' ")"
                }

                local alt_M2 = `M0'*`S'
                cap assert abs(`M2' - `alt_M2') < `assertion_sensitivity' | ///
                           (missing(`M2') | missing(`alt_M2'))
                if _rc != 0 {
                    noi display as error "Computation M2 failed (`alt_M2' vs. " `M2' ")"
                }
            }
        }


        /*
            Calculating proportional results in the case of there being a by-variable.
        */

        if "`by'" != "" {

            // Naming the matrix (for later)
            local levelnames = e(over_namelist)

            // Estimate the proportional contribution
            tempname res_perc res_percV

            // Ratios to evaluate in the end.
            local ratios = ""

            local lvl_idx = 0
            foreach i of local levelnames {
                local lvl_idx = `lvl_idx' + 1
                tempvar l_H_`lvl_idx'
                gen `l_H_`lvl_idx'' = `_H' * (`by' == `i')
                local ratios = "`ratios' (`l_H_`lvl_idx''/`_H')"
            }

            foreach i of local levelnames {
                local lvl_idx = `lvl_idx' + 1
                tempvar l_M0_`lvl_idx'
                gen `l_M0_`lvl_idx'' = `_M0' * (`by' == `i')
                local ratios = "`ratios' (`l_M0_`lvl_idx''/`_M0')"
            }

            // M_alpha
            if `use_thresholds' {
                foreach a of local alpha {
                    local a_name = subinstr("`a'", ".", "p", .)
                    local a_name = subinstr("`a_name'", "-", "m", .)
                    foreach i of local levelnames {
                        local lvl_idx = `lvl_idx' + 1
                        tempvar l_M`a_name'_`lvl_idx'
                        gen `l_M`a_name'_`lvl_idx'' = `gap_`a_name'' * (`by' == `i')
                        local ratios = "`ratios' (`l_M`a_name'_`lvl_idx'' / `gap_`a_name'')"
                    }
                }
            }

            `svy' ratio `ratios' `weight_exp' `l_if'
            matrix `res_perc' = e(b)
            matrix `res_percV' = e(V)

            /*
                There (used to?) exist a bug where a ratio of (x / y) with y=0 returns invalid
                values in e(V). This code here will detect the invalid entries and replace them
                with zero.
                Then, the behavior is equivalent to a zero-variation entry.
                See:
                http://www.statalist.org/forums/forum/general-stata-discussion/general/1348095-create-or-test-for-1-ind-in-a-matrix
            */
            local cols = colsof(`res_percV')
            forvalues i = 1/`cols' {
                forvalues j = 1/`cols' {
                    // The entry is not missing if unmodified, but if
                    // used in a formula, it becomes missing.
                    if missing(1 * `res_percV'[`i', `j']) & !missing(`res_percV'[`i', `j']) {
                        matrix `res_percV'[`i', `j'] = 0
                    }
                }
            }

        } // End by exists?

        /*
            Name the absolute results created above.
        */
        local name_main = ""
        local name_add = ""
        if "`by'" == "" {
            // Labels and output of estimate matrix.
            local by_n_1 = `"Main"'
            local by_n_2 = `"Additional"'
            if `use_thresholds' == 0 {
                local name_main = `"`name_main' "Main:H" "Main:M0""'
                local name_add = `"`name_add' "Additional:A""'
            }
            else {
                local name_main = `"`name_main' "Main:H" "Main:M0" "Main:M1" "Main:M(2)""'
                local name_add = `"`name_add' "Additional:A" "Additional:G" "Additional:S(2)""'

                foreach a of local alpha {
                    if inlist(`a', 1, 2) {
                        continue
                    }
                    local a_name = subinstr("`a'", ".", "p", .)
                    local a_name = subinstr("`a_name'", "-", "m", .)

                    // `gap_`a_name'' "Adjusted Multidimensional Deprivation Gap with transformation alpha"
                    local main_eval = "`main_eval' `gap_`a_name''"
                    local add_eval = "`add_eval' (`gap_`a_name'' / `_M0')"

                    local a_name = subinstr("`a'", ".", ",", .)
                    local name_main = `"`name_main' "Main:M(`a_name')""'
                    local name_add = `"`name_add' "Additional:S(`a_name')""'
                } // end alpha iteration
            }
        } // End no by.
        else {
            local levels = `"`levelnames'"'

            if `use_thresholds' == 0 {
                foreach by_level of local levels {
                    local name_main = `"`name_main' "H:`by'_`by_level'" "M0:`by'_`by_level'""'
                }
            }
            else {
                local names = "H M0"
                foreach a of local alpha {
                    local a_name = subinstr("`a'", ".", "p", .)
                    local a_name = subinstr("`a_name'", "-", "m", .)
                    local names = "`names' M(`a_name')"
                } // end alpha iteration

                foreach name of local names {
                    foreach by_level of local levels {
                        local name_main = `"`name_main' "`name':`by'_`by_level'""'
                    }
                }
            }
        }

        /*
            Apply the naming for the standard results.
        */
        matrix rownames `res' = "Estimate"
        matrix colnames `res' = `name_main'
        matrix colnames `resV' = `name_main'
        matrix rownames `resV' = `name_main'

        // Additional results are only computed in the no-by case.
        if "`by'" == "" {
            matrix rownames `res_add' = "Estimate"
            matrix colnames `res_add' = `name_add'
            matrix colnames `resV_add' = `name_add'
            matrix rownames `resV_add' = `name_add'
        }
        // Proportional contributions only matter in the by-case.
        else {
            matrix rownames `res_perc' = "Estimate"
            matrix colnames `res_perc' = `name_main'
            matrix colnames `res_percV' = `name_main'
            matrix rownames `res_percV' = `name_main'
        }

        /*
            Return results.
        */

        return matrix mpi_main = `res'
        return matrix mpi_main_V = `resV'
        if "`by'" == "" {
            return matrix mpi_add = `res_add'
            return matrix mpi_add_V = `resV_add'
        }
        else {
            return matrix mpi_perc = `res_perc'
            return matrix mpi_perc_V = `res_percV'
            return local over_namelist `levelnames'
        }
        return scalar N = `N'

        // Remove the ereturn results that this program created.
        ereturn clear
    } // End quietly

end

/* ---------------------------------------------------------------------------

Decomposition of M_alpha by Indicator

    - If the by-option is specified, ONLY the by-version is returned.
    - Do not call independently.
--------------------------------------------------------------------------- */

program define _mpi_domains, sortpreserve rclass

    syntax anything(equalok)

    quietly {
        // Read in all the to-be-used locals supplied to the program.
        foreach pair in `*' {
            local eq_pos = strpos("`pair'", "=")
            local before = substr("`pair'", 1, `eq_pos'-1)
            local after = substr("`pair'", `eq_pos'+1, .)
            local `before' = "`after'"
        }

        /*
            If a by-option was specified, prepare the respective statements.
        */
        if "`by'" != "" {
            local over_statement = ", over(`by', nolabel)"
        }
        else {
            local over_statement = ""
        }

        tempname m m_V

        local tot_m0_shares = 0
        foreach a of local alpha {
            local a_name = subinstr("`a'", ".", "p", .)
            local a_name = subinstr("`a_name'", "-", "m", .)
            local tot_m`a_name'_shares = 0
        }
        local prev_ind = 0

        // The ratios to finally evaluate.
        local eval_ratio = ""

        // Iterate over domains
        forvalues j = 1 / `ndom' {
            local nind`j' = wordcount("`d`j''")

            /*
                Iterate over indicators, create the to-be-estimated variables.
            */
            forvalues i = 1 / `nind`j'' {
                local w = word("`w`j''", `i')
                local ind = word("`d`j''", `i')

                tempname m_`j'_`i'
                generate `m_`j'_`i'' = `ispoor_`j'_`i'' * `_H'

                // Generate the M_alpha as a share of total M0 for this indicator
                local eval_ratio = "`eval_ratio' (`m_`j'_`i'' / `_M0')"

                if `use_thresholds' {
                    tempname gap_poor
                    local a_idx = 1
                    foreach a of local alpha {
                        local a_idx = `a_idx' + 1
                        local a_name = subinstr("`a'", ".", "p", .)
                        local a_name = subinstr("`a_name'", "-", "m", .)
                        // Generate Malpha as share of total Malpha
                        // = (|galpha(select)|/|g0(total)|)/Malpha(total)
                        tempname gap_poor_`j'_`i'_`a_name'
                        gen `gap_poor_`j'_`i'_`a_name'' = `gap_`j'_`i'_`a_name'' * `_H'

                        if `a' == 1 {
                            local eval_ratio = "`eval_ratio' (`gap_poor_`j'_`i'_`a_name'' / `gap_1')"
                        }
                        else if `a' == 2 {
                            local eval_ratio = "`eval_ratio' (`gap_poor_`j'_`i'_`a_name'' / `gap_2')"
                        }
                        else {
                            local eval_ratio = "`eval_ratio' (`gap_poor_`j'_`i'_`a_name'' / `gap_`a_name'')"
                        }
                    }
                }
            }
        } // End iterating over domains

        // Calculating ratios.
        `svy' ratio `eval_ratio' `weight_exp' `l_if' `over_statement'
        matrix `m' = e(b)
        matrix `m_V' = e(V)

        /*
            There (used to?) exist a bug in Stata where a ratio of (x / y) with y=0 returns invalid
            values in e(V). This code here will detect the invalid entries and replace them
            with zero.
            Then, the behavior is equivalent to a zero-variation entry.
            See:
            http://www.statalist.org/forums/forum/general-stata-discussion/general/1348095-create-or-test-for-1-ind-in-a-matrix
        */
        local cols = colsof(`m')
        forvalues i = 1/`cols' {
            forvalues j = 1/`cols' {
                // The entry is not missing if unmodified, but if
                // used in a formula, it becomes missing.
                if missing(1 * `m_V'[`i', `j']) & !missing(`m_V'[`i', `j']) {
                    matrix `m_V'[`i', `j'] = 0
                }
            }
        }

        // Naming the matrix (for later)
        local levelnames = e(over_namelist)

        /*
            The problem at this point is that the results should be
            contributions. The weights need to be included to make the
            contributions add up to one.
        */
        tempname mat_weights
        local n_levels = cond("`by'" == "", 1, wordcount(`"`=e(over_namelist)'"'))
        local size = (wordcount("`alpha'") * `use_thresholds' + 1) * `ni' * `n_levels'
        matrix `mat_weights' = J(`size', `size', 0)
        local idx = 0
        forvalues j = 1 / `ndom' {
            forvalues i = 1 / `nind`j'' {
                if `use_thresholds' {
                    local iterate = "0 `alpha'"
                }
                else {
                    local iterate = "0"
                }
                foreach a of local iterate {
                    forvalues level = 1 / `n_levels' {
                        local idx = `idx' + 1
                        matrix `mat_weights'[`idx', `idx'] = `=word("`w`j''", `i')'
                    }
                }
            }
        }
        matrix `m' = `m' * `mat_weights'
        matrix `m_V' = `mat_weights'' * `m_V' * `mat_weights'

        if "`by'" == "" {
            local levels = `" "" "'
        }
        else {
            local levels = `"`levelnames'"'
        }

        // Holds the names for the estimated matrices.
        local names = ""

        /*
            Iterate over domains/indicators, create the labels.
        */
        forvalues j = 1 / `ndom' {
            forvalues i = 1 / `nind`j'' {
                local ind = word("`d`j''", `i')

                /*
                    Create the names for each indicator within each domain
                */
                if `use_thresholds' {
                    local idx = 0
                    local addlabel_alpha = ""
                    foreach a in 0 `alpha' {
                        local idx = `idx' + 1
                        local a_name = subinstr("`a'", ".", ",", .)
                        if "`a'" != "0" {
                            local a_name = "(`a_name')"
                        }
                        foreach by_level of local levels {
                            // Labels and output of estimate matrix.
                            if `"`by_level'"' != "" {
                                local by_n = `"`by'_`by_level':"'
                            }
                            else {
                                local by_n = `""'
                            }
                            local names = `"`names' "`by_n'`ind'_M`a_name'""'
                        } // Iterating over by-levels (if they exist)
                    } // Iterate over alpha levels.
                }
                else {
                    foreach by_level of local levels {
                        // Labels and output of estimate matrix.
                        if `"`by_level'"' != "" {
                            local by_n = `"`by'_`by_level':"'
                        }
                        else {
                            local by_n = `""'
                        }
                        local names = `"`names' "`by_n'`ind'_M0""'
                    } // Iterating over by-levels (if they exist)
                }
            } // Iterating over indicators
        } // Iterating over domains.

        /*
            Put the names.
        */
        matrix colnames `m' = `names'
        matrix colnames `m_V' = `names'
        matrix rownames `m_V' = `names'

        // Erase ereturn results created within this program
        ereturn clear

        return matrix mpi_decomposed = `m'
        return matrix mpi_decomposed_V = `m_V'
        if "`by'" != "" {
            return local over_namelist `levelnames'
        }
    } // End quietly

end

/*

    Full mpi program.

*/

program mpi, sortpreserve eclass

    cap syntax anything [in] [if] [pw fw], Cutoff(real) ///
        [by(varname) Alpha(numlist) ///
        Level(real -1) CATegories(integer 20) ///
        Svy SUBpop(passthru) ///
        NOSummary NODecomposition ///
        POSTMain POSTAdditional POSTIndicators POSTDomains ///
        POSTBYMain POSTBYProportion POSTBYIndicators POSTBYDomains ///
        DEPRIVEDDummy(namelist min=1 max=1) DEPRIVEDScore(namelist min=1 max=1)]

    // Take care of error messages when calling mpi through bootstrap.
    if _rc == 100 {
        local cmd = e(cmd)
        if `"`cmd'"' == "mpi" {
            noi display as text "Replaying last e(b) / e(V) results as mpi was started without options."
            noi display as text "This is normal for {help bootstrap}. Otherwise, see {help mpi} on how to specify options."
            noi cap noi ereturn display
            if _rc == 1 {
                error 1
            }
            else if _rc != 0 {
                noi _coef_table
            }
            exit
        }
        else {
            error _rc
        }
    }
    else {
        syntax anything [in] [if] [pw fw], Cutoff(real) ///
                [by(varname) Alpha(numlist) ///
                Level(real -1) CATegories(integer 20) ///
                Svy SUBpop(passthru) ///
                NOSummary NODecomposition ///
                POSTMain POSTAdditional POSTIndicators POSTDomains ///
                POSTBYMain POSTBYProportion POSTBYIndicators POSTBYDomains ///
                DEPRIVEDDummy(namelist min=1 max=1) DEPRIVEDScore(namelist min=1 max=1)]
    }

    quietly {

        version 13.0

        // Check that the correct amount of post-statements was given: zero or one.
        local cnt_list = 0
        foreach a in postmain postadditional postindicators postdomains ///
                 postbymain postbyproportion postbyindicators postbydomains {
            if "``a''" != "" {
                local cnt_list = `cnt_list'+1
            }
        }
        if `cnt_list' > 1 {
            noi display as error "Error: Supplied more than one option to list results in e(b) and e(V)." ///
                " This is not allowed, specify at maximum one."
            error 198
        }
        // Choose postmain as default option.
        else if `cnt_list' == 0 {
            local postmain = "postmain"
        }

        // Check that, if given, the deprived dummy/score variable does not exist.
        if "`depriveddummy'" != "" {
            noi confirm new var `depriveddummy'
        }
        if "`deprivedscore'" != "" {
            noi confirm new var `deprivedscore'
        }

        // Globally required names present everywhere
        tempvar sweight dscore gap_1 gap_2 _H _M0

        local var_submit = `"sweight=`sweight'"'
        local var_submit = `"`var_submit' dscore=`dscore'"'
        local var_submit = `"`var_submit' gap_1=`gap_1'"'
        local var_submit = `"`var_submit' gap_2=`gap_2'"'
        local var_submit = `"`var_submit' _H=`_H'"'
        local var_submit = `"`var_submit' _M0=`_M0'"'

        // Mark sample
        marksample touse

        // Prepare weights
        if "`exp'" != "" {
            generate `sweight' `exp'
        }
        else {
            generate `sweight' = 1
        }
        if "`weight'`exp'" != "" {
            local weight_exp = "[`weight'`exp']"
        }
        else {
            local weight_exp = ""
        }

        // Prepare svy usage
        if "`svy'" != "" {
            // Standard If and In are not allowed with svy.
            if "`if'" != "" & "`in'" != "" {
                noi display as error "Error: The svy option cannot be combined with if and in. Use the extended svy option."
                error 101
            }
            else if "`if'" != "" {
                noi display as error "Error: The svy option cannot be combined with if. Use the extended svy option."
                error 101
            }
            else if "`in'" != "" {
                noi display as error "Error: The svy option cannot be combined with in. Use the extended svy option."
                error 101
            }
            else if `"`subpop'"' != "" {
                // Get content between braces for subpop clause.
                local firstbrace = strpos(`"`subpop'"', "(")
                local lastbrace = -1
                forvalues i = 1/`=strlen(`"`subpop'"')' {
                    if substr(`"`subpop'"', `i', `i'+1) == ")" {
                        local lastbrace = `i'
                    }
                }
                // If the subpop clause is a variable, change it into an if statement.
                local subpop_clause = substr(`"`subpop'"', `firstbrace'+1, `lastbrace'-`firstbrace'-1)
                capture confirm variable `clause'
                if _rc == 0 {
                    local clause = `"if `clause'"''
                }

                local subpop_clause = substr(`"`subpop'"', 1, `firstbrace') + `"`clause'"' + substr(`"`subpop'"', `lastbrace', .)
                local svy = `"svy, `subpop': "'
            }
            else {
                local svy = "svy:"
            }
            local l_if = ""
        }
        // Standard: no svy, if as in touse
        else {
            local svy = ""
            local l_if = "if `touse'"
        }

        if `"`subpop'"' != "" & "`svy'" == "" {
            noi display as error "Error: The subpop option requires the svy option to be specified."
            error 101
        }

        // Prepare CI level
        if `level' == -1 {
            local level = c(level)
        }

        // As long as no weight information is given, assume weighting to be equal.
        local equalw = 1

        // Maximum Variable Length - to be generated later
        local max_var_length = -1

        // Sensitivity for all assumptions
        local assertion_sensitivity = 0.00001

        // To compute:
        // Number of indicators
        local ni = 0
        // Number of domains
        local ndom = 0
        // nidX: Number of indicators in each domain
        // dX: Domain X
        // wX: Weights X
        // thresX: Thresholds X
        // Whether thresholds for the indicators are given
        local use_thresholds = 0

        while strpos("`anything'", ")") > 0 {
            local entry = trim(substr("`anything'", 1, strpos("`anything'", ")")))
            local anything = substr("`anything'", strpos("`anything'", ")") + 1, .)
            // Argument type should be "w" or "d"
            local argtype = lower(substr(trim("`entry'"), 1, 1))
            if !(inlist("`argtype'", "w", "weight", "d", "dom", "domain") ///
                 | inlist("`argtype'", "t", "thres", "threshold", "p", "poor")) {
                noi display as error "Only domains (d), weights (w) or thresholds (t) can be supplied. Found: '`argtype'"
            }
            // Get the number of weight / indicator
            local brstartpos = strpos("`entry'", "(")
            local idx = trim(substr("`entry'", 2, `brstartpos'-2))
            capture confirm number `idx'
            if _rc != 0 {
                noi display as error "Argument '`entry'' did not follow the proper form, after d/w should follow the index."
                error 197
            }
            else {
                error _rc
            }
            // Get content of brackets
            local inbr = substr("`entry'", `brstartpos' + 1, strpos("`entry'", ")") - `brstartpos' - 1)
            // Check that the variables exist if it is a domain
            if "`argtype'" == "d" | "`argtype'" == "dom" | "`argtype'" == "domain" {
                local count_ind = 0
                unab inbr : `inbr'
                foreach d of varlist `inbr' {
                    confirm var `d'
                    local count_ind = `count_ind' + 1
                    local max_var_length = max(`max_var_length', strlen("`d'"))
                }
                // Save the number of indicators for this domain
                local nid`idx' = `count_ind'
                local ni = `ni' + `nid`idx''
                local argtype = "d"
            }
            if "`argtype'" == "w" | "`argtype'" == "weight" {
                foreach w of local inbr {
                    confirm number `w'
                }
                local equalw = 0
                local argtype = "w"
            }
            if "`argtype'" == "t" | "`argtype'" == "thres" | "`argtype'" == "threshold" {
                /* Threshold definitions:
                An individual is poor if the value of the indicator is below the threshold.
                A person at the threshold is always not poor.
                For this, all the threshold arguments need to be specified without
                spaces in between.                                              */
                local inbr_new = ""
                foreach t of local inbr {
                    confirm number `t'
                    local inbr_new = "`inbr_new' `t'"
                }
                local use_thresholds = 1
                local argtype = "thres"
                local inbr = "`inbr_new'"
            }

            // Set weight / indicator
            local `argtype'`idx' = "`inbr'"
            local ndom = max(`idx', `ndom')
        }

        // Make sure that at least 1 domain is specified.
        if `ndom' == 0 {
            noi display as error "No deprivation domains specified"
            error 197
        }

        // If no thresholds are in use, every threshold is, by definition, 1.
        // >0.99 is used instead of 1 since the poverty is defined as being strictly
        // lower or higher than the cutoff, which would otherwise cause problems in
        // datasets without explicit thresholds.
        if `use_thresholds' == 0 {
            forvalues i = 1 / `ndom' {
                local thres`i' = ""
                forvalues j = 1/`nid`i'' {
                    local thres`i' = "0.99 `thres`i''"
                }
            }
        }

        // Make sure that no domains or weights are missing
        forvalues idx = 1/`ndom' {
            if ("`d`idx''" == "") {
                noi display as error "Deprivation domain `idx' not specified"
                error 197
            }
            if (`equalw' == 0) & ("`w`idx''" == "") {
                noi display as error "Weight `idx' not specified."
                error 197
            }
            if (`equalw' == 0) & (wordcount("`w`idx''") != `nid`idx'') {
                noi display as error "Vector w`idx'(`w`idx'') does not have the same number of elements as the number of indicators in domain d`idx'(`d`idx'')"
                error 197
            }
            // Check that all thresholds are given, always.
            if (wordcount("`thres`idx''") != `nid`idx'') {
                noi display as error "Vector thres`idx'(`thres`idx'') does not have the same number of elements as the number of indicators in domain d`idx'(`d`idx'')"
                error 197
            }
        }

        // Make sure that alpha does not include 0, 1 and 2, they are computed by default.
        local newalpha = ""
        foreach a of local alpha {
            if !inlist(`a', 0, 1, 2) {
                local newalpha = "`newalpha' `a'"
            }
        }
        local alpha = "1 2 `newalpha'"

        // Check if any observations are left
        tempname all_missing
        gen `all_missing' = 0
        forvalues i = 1 / `ndom' {
            foreach vari of varlist `d`i'' {
                replace `all_missing' = `all_missing' | missing(`vari')
            }
        }
        if "`by'" != "" {
            replace `all_missing' = 1 if missing(`by') & `touse'
        }
        count if `all_missing' & `touse'
        if r(N) == _N {
            // No observations
            error 2000
        }

        // Exclude missing values from the estimation sample
        count if `all_missing' == 1
        if r(N)>0 {
            noi display in gr "Note: Missing values encountered, excluding them."
            recode `all_missing' 1=.
            markout `touse' `sweight' `all_missing'
        }
        drop `all_missing'

        // Make sure that all variables are ordinal or all variables are real-valued.
        // Guess the type of each variable by the number of observations
        local error_message = ""
        local t_vars = ""
        forvalues i=1/4 {
            local msg_`i'_shown = 0
        }
        tempvar dist_obs
        forvalues i = 1 / `ndom' {
            local vartype`i' = ""
            foreach vari of varlist `d`i'' {
                bys `vari': gen `dist_obs' = (_n == 1)
                count if `dist_obs'
                local n_dist_obs = r(N)
                // Likely ordinal or even binary, see syntax for the origin of
                // the local
                if `n_dist_obs' <= `categories' {
                    // Set the expected variabe type for subsequent iterations if it is not set.
                    if "`t_vars'" == "" {
                        if (`n_dist_obs' == 2) {
                            local t_vars = "Binary"
                        }
                        else {
                            local t_vars = "Ordinal"
                        }
                    }
                    // The first detected variable was real-valued but now it looks binary.
                    else if ("`t_vars'" == "Real-valued") & (`n_dist_obs' <= 2) {
                        if !`msg_1_shown' {
                            local error_message = `"`error_message' ""' + ///
                                    `"Mixing binary (`vari') and real-valued indicators is not recommended: binary indicators automatically receive a higher weight during the calculation of the normalized poverty gaps. See {help mpi} for details." _n "'
                        }
                        local msg_1_shown = 1
                    }
                    // The first detected variable was real-valued but now it looks ordinal.
                    else if ("`t_vars'" == "Real-valued") & (`n_dist_obs' != 2) {
                        if !`msg_3_shown' {
                            local error_message = `"`error_message' ""' + ///
                                    `"Variable `vari' has little variation (less than `categories' categories) whereas other variables seem real-valued. "' + ///
                                    `"The statistical properties of AF indices with alpha>0 hold when the distance from the deprivation threshold can be measured correctly. "' + ///
                                    `"With ordinal and categorical variables the use of M0 (after dichotomisation into deprived and non-deprived states) is recommended.  See {help mpi} for details." _n"'
                        }
                        local msg_4_shown = 1
                    }
                    if `use_thresholds' {
                        if (`n_dist_obs' == 2) {
                            local vartype`i' = "`vartype`i'' Binary"
                        }
                        else {
                            local vartype`i' = "`vartype`i'' Ordinal"
                        }
                    }
                    else {
                        local vartype`i' = "`vartype`i'' Binary"
                    }
                }
                // The variable is likely real-valued
                else {
                    if "`t_vars'" == "" {
                        local t_vars = "Real-valued"
                    }
                    else if "`t_vars'" != "Real-valued" {
                        if !`msg_4_shown' {
                            local error_message = `"`error_message' "Variable `vari' is likely real-valued, but previous variables had little variation (less than `categories' categories). "' + ///
                                    `"The statistical properties of AF indices with alpha>0 hold when the distance from the deprivation threshold can be measured correctly. "' + ///
                                    `"With ordinal and categorical variables the use of M0 (after dichotomisation into deprived and non-deprived states) is recommended. See {help mpi} for details." _n"'
                        }
                        local msg_4_shown = 1
                    }
                    local vartype`i' = "`vartype`i'' Real-valued"
                }
                drop `dist_obs'
            }
        }

        // If equal weighting is specified, assign equal weights for every indicator.
        if `equalw' {
            local domshare = 1 / `ndom'
            forvalues i = 1 / `ndom' {
                local w`i' = ""
                local indshare = `domshare'/`nid`i''
                forvalues j = 1/`nid`i'' {
                    local w`i' = "`indshare' `w`i''"
                }
            }
        }

        confirm number `cutoff'

        // Check that relevant variables are numeric and get the number of missing values
        // Check that each variable contains only 0, 1 and missing entries if it is the
        // non-threshold case.
        forvalues j=1/`ndom' {
            foreach ind of varlist `d`j'' {
                capture confirm numeric variable `ind'
                if _rc != 0 {
                    noi display as error "Variable `ind' is not numeric."
                    exit 459
                }
                else {
                    error _rc
                }
                // If no thresholds are used, 0 is not poor and 1 is poor. Hence,
                // check that the variable contents match to the case.
                if `use_thresholds' == 0 {
                    capture assert inlist(`ind', 0, 1) | missing(`ind') | !`touse'
                    if _rc == 9 {
                        noi display as error "Variable `ind' contains values besides " ///
                            "0, 1 and missing. Supply indicators containing 1 for poor " ///
                            "and 0 for non-poor individuals or supply poverty thresholds."
                        exit 459
                    }
                    else {
                        error _rc
                    }
                }
            }
        }

        // Check if the supplied domainal weights add up to one:
        local totwcheck = 0
        forvalues j = 1/`ndom' {
            local countw = 1
            foreach ind of varlist `d`j'' {
                local w = word("`w`j''", `countw')
                local totwcheck = `totwcheck' + `w'
                local countw = `countw' + 1
            }
        }
        if abs(`totwcheck' - 1) > .01 {
            noi display in re "Total sum of weight is " _c
            noi display in ye `totwcheck' in re ", it should be 1."
            error 9
        }
        // If the difference between the weights is close to 1, scale them and the
        // cutoff to remove even the slight differences.
        else if abs(`totwcheck' - 1) > `assertion_sensitivity' {
            noi display in ye "Note: the total sum of weight was close to 1, the weights and the cutoff were rescaled to add up to exactly 1."
            forvalues j = 1/`ndom' {
                local new_weight = ""
                local countw = 1
                foreach ind of varlist `d`j'' {
                    local w = word("`w`j''", `countw')
                    local w  = `w' / `totwcheck'
                    local new_weight = "`new_weight' `w'"
                    local countw = `countw' + 1
                }
                local w`j' = "`new_weight'"
            }
            local cutoff = `cutoff' / `totwcheck'
        }

        if `totwcheck' < `cutoff' & "`noinitialcall'" == "" {
            noi display in re "Attention, cutoff is larger than the sum of weights. " _c
            noi display in ye "(C: `cutoff' W: 1). " _c
            noi display in re "No individual could possibly be poor."
        }

        // Create matrix of deprivation domains and overall deprivation score
        foreach a of local alpha {
            local a_name = subinstr("`a'", ".", "p", .)
            local a_name = subinstr("`a_name'", "-", "m", .)
            tempname gap_`a_name'
            generate `gap_`a_name'' = 0
            local var_submit = `"`var_submit' gap_`a_name'=`gap_`a_name''"'
        }
        generate `dscore' = 0
        forvalues j = 1/`ndom' {
            tempvar wdom`j'
            foreach a of local alpha {
                local a_name = subinstr("`a'", ".", "p", .)
                local a_name = subinstr("`a_name'", "-", "m", .)
                tempvar gap`j'_`a_name'
                generate `gap`j'_`a_name'' = 0
                local var_submit = `"`var_submit' gap`j'_`a_name'=`gap`j'_`a_name''"'
            }
            generate `wdom`j'' = 0
            c_local wdom`j' "`wdom`j''"
            local nind`j' = wordcount("`d`j''")
            forvalues i = 1 / `nind`j'' {
                local w = word("`w`j''", `i')
                local ind = word("`d`j''", `i')
                local threshold = word("`thres`j''", `i')
                // Poor if indicator < Threshold.
                // If no thresholds are defined: special case, threshold is at 0.99 and
                // above is poor (as the indicators will have value 1)
                // Always: Being at the threshold is NOT poor.
                local comparison = cond(`use_thresholds' == 0, 1, -1)
                tempvar ispoor_`j'_`i'
                generate `ispoor_`j'_`i'' = (`comparison'*(`ind' - `threshold') > 0) if !missing(`ind')
                local var_submit = `"`var_submit' ispoor_`j'_`i'=`ispoor_`j'_`i''"'
                replace `wdom`j'' = `wdom`j'' + `ispoor_`j'_`i''*`w'
                // Generate poverty gaps. Absolute value is to prevent negative values, which
                // can occur if the threshold is < 0.
                tempvar gap_`j'_`i'_1
                // Gap = [Poor]*[Achievement gap]/[Threshold]
                // Assuming that an individual is poor IF the indicator is higher
                // than the threshold. If the comparison is the other way round because of binary variables,
                // multiply by -1. (Done through *`comparison')
                generate `gap_`j'_`i'_1' = abs(`ispoor_`j'_`i''*`comparison'*(`ind'-`threshold') / `threshold')
                local var_submit = `"`var_submit' gap_`j'_`i'_1=`gap_`j'_`i'_1'"'
                replace `gap`j'_1' = `gap`j'_1' + `gap_`j'_`i'_1'*`w'
                foreach a of local alpha {
                    if `a' != 1 {
                        local a_name = subinstr("`a'", ".", "p", .)
                        local a_name = subinstr("`a_name'", "-", "m", .)
                        tempvar gap_`j'_`i'_`a_name'
                        generate double `gap_`j'_`i'_`a_name'' = `gap_`j'_`i'_1' ^ `a'
                        local var_submit = `"`var_submit' gap_`j'_`i'_`a_name'=`gap_`j'_`i'_`a_name''"'
                        cap assert `gap_`j'_`i'_`a_name'' >= 0
                        if _rc == 9 {
                            noi display as error "Problem in comparison, negative gap detected."
                        }
                        else {
                            error _rc
                        }
                        replace `gap`j'_`a_name'' = `gap`j'_`a_name'' + `gap_`j'_`i'_`a_name''*`w'
                    }
                }
            }
            replace `dscore' = `dscore' + `wdom`j''
            foreach a of local alpha {
                local a_name = subinstr("`a'", ".", "p", .)
                local a_name = subinstr("`a_name'", "-", "m", .)
                replace `gap_`a_name'' = `gap_`a_name'' + `gap`j'_`a_name''
            }
        }

        foreach a of local alpha {
            local a_name = subinstr("`a'", ".", "p", .)
            local a_name = subinstr("`a_name'", "-", "m", .)
            replace `gap_`a_name'' = 0 if `dscore' < `cutoff'
        }

        // Generate Multidimensional Deprivation Headcount (H)
        generate `_H' = (`dscore' >= `cutoff') if !missing(`dscore')

        // If the user wants a _H variable to be generated, do it.
        if "`depriveddummy'" != "" {
            generate `depriveddummy' = `_H' if `touse'
            label var `depriveddummy' "mpi: Dummy multidimensionally deprived"
            label define `depriveddummy'_deprived 0 "No" 1 "Yes", modify
            label values `depriveddummy' `depriveddummy'_deprived
        }

        // If the user wants a deprivation variable to be generated, do it.
        if "`deprivedscore'" != "" {
            generate `deprivedscore' = `dscore' if `touse'
            label var `deprivedscore' "mpi: Multidimensional deprivation score"
        }


        // Generate Adjusted Multidimensional Deprivation Headcount (M0)
        generate `_M0' = cond(`_H' > 0, `dscore', 0)

        /*
            Return all locals generated over the course of the program.
        */
        foreach loc_name in weight_exp svy equalw max_var_length assertion_sensitivity ///
                            ni ndom use_thresholds alpha level touse l_if svy {
            local var_submit = `"`var_submit' "`loc_name'=``loc_name''" "'
        }
        forvalues i = 1 / `ndom' {
            local var_submit = `"`var_submit' "thres`i'=`thres`i''""'
            local var_submit = `"`var_submit' "d`i'=`d`i''""'
            local var_submit = `"`var_submit' "w`i'=`w`i''""'
            local var_submit = `"`var_submit' "nind`i'=`nind`i''""'
            local var_submit = `"`var_submit' "vartype`i'=`vartype`i''""'
        }

        // A is a local for naming mata matrices, res for temporarily saving results.
        tempname res A

    /* ---------------------------------------------------------------------------

    Create Summary Table

    --------------------------------------------------------------------------- */
        // Establish total number of individuals
        if "`nosummary'" == "" {

            noi display as result _n _n "Summary of {cmd:mpi} indicators" _n "{hline 26}"

            local table_shift = max(5 + `max_var_length', 11)
            local table_1 = `table_shift' + 13 - (1-`use_thresholds')*4
            local table_2 = `table_1' + 10
            if `use_thresholds' {
                local table_3 = `table_2' + 10
            }
            else {
                local table_3 = `table_1' + 10
            }
            local linew = `table_3' + 7

            noi display as text "Indicator" _c
            noi display as text _col(`table_shift') "Type" _c
            noi display as text _col(`table_1') "Weight" _c
            if `use_thresholds' {
                noi display as text _col(`table_2') "Threshold" _c
            }
            noi display as text _col(`table_3') "Deprived" _c
            noi display as text _n "{hline `linew'}" _c

            forvalues j = 1 / `ndom' {
                local nind`j' = wordcount("`d`j''")
                noi display as result _n "Domain `j'" _c
                forvalues i = 1 / `nind`j'' {
                    // Extract weight, indicator and thresholds from
                    // the respective locals
                    local w = round(real(word("`w`j''", `i')), 0.01)
                    local ind = word("`d`j''", `i')
                    local threshold = word("`thres`j''", `i')
                    local vartype = word("`vartype`j''", `i')

                    // Calculate how many people are deprived.
                    `svy' proportion `ispoor_`j'_`i'' `weight_exp' `l_if'
                    matrix `res' = e(b)
                    local perc_poor = `res'[1, 2] * 100
                    if missing(`perc_poor') {
                        local perc_poor = 0
                    }

                    noi display as text _n _col(3) "`ind'" _c
                    noi display as text _col(`table_shift') "`vartype'" _c
                    noi display as result _col(`table_1') "`w'" _c
                    if `use_thresholds' {
                        noi display as result _col(`table_2') "`threshold'" _c
                    }
                    noi display as result _col(`table_3') %6.3f = `perc_poor' " %" _c
                }
            }
            noi display as text _n "{hline `linew'}"
            noi display as text "Deprived: Percentage of individuals whose indicator values are"
            noi display as text "below the threshold." _n _n

            if `"`error_message'"' != "" {
                noi display as result "Note: " `error_message'
            }
        }

    /* ---------------------------------------------------------------------------

    Computing Main Indicators

    --------------------------------------------------------------------------- */

        /*
            Check if all individuals are non-poor, as this requires additional
            messages below.
        */
        cap assert `_H' == 0
        if _rc == 9 {
            local warn_no_poor = 0
        }
        else {
            error _rc
            local warn_no_poor = 1
        }

        noi _mpi_main `"`var_submit'"'

        tempname hm_1 hm_2 hm_V_1 hm_V_2

        matrix `hm_1' = r(mpi_main)
        matrix `hm_V_1' = r(mpi_main_V)
        matrix `hm_2' = r(mpi_add)
        matrix `hm_V_2' = r(mpi_add_V)
        local N = r(N)

        /*
            Intermission: Display the results.
        */
        tempname disp_b disp_V tmp

        local hm_1_col = colsof(`hm_1')
        local hm_2_col = colsof(`hm_2')
        matrix `disp_b' = `hm_1', `hm_2'
        matrix `tmp' = J(`hm_1_col', `hm_2_col', 0)
        local names : colfullnames `hm_1'
        matrix rownames `tmp' = `names'
        local names : colfullnames `hm_2'
        matrix colnames `tmp' = `names'
        matrix `disp_V' = (`hm_V_1', `tmp')
        matrix `disp_V' = `disp_V' \ (`tmp'', `hm_V_2')
        ereturn post `disp_b' `disp_V', obs(`N')

        local n_string = "N = `N'"
        local len_n_string = strlen("`n_string'")
        noi display as text "Main results" _col(`=62-`len_n_string'') "`n_string'"

        // Stata 14 and newer versions of Stata 13.
        quietly capture ereturn display, level(`level') nopvalue cformat("%9.3f")
        if _rc == 1 {
            error 1
        }
        else if _rc != 0 {
            // Older versions.
            quietly capture _coef_table, cionly level(`level') cformat("%9.3f")
            if _rc == 1 {
                error 1
            }
            else if _rc != 0 {
                // Final fallback.
                tempname hm hm_se hm_se_ci
                matrix `hm' = e(b)
                matrix `hm' = `hm''
                matrix `hm_se' = e(V)
                matrix `hm_se' = vecdiag(`hm_se'')
                mata: `A' = st_matrix("`hm_se'")
                mata: `A' = sqrt(`A')
                mata: st_matrix("`hm_se'", `A')
                matrix `hm_se' = `hm_se''
                matrix `hm_se_ci' = J(`=rowsof(`hm')', 2, .)
                forvalues i = 1/`=rowsof(`hm')' {
                    matrix `hm_se_ci'[`i', 1] = `hm'[`i', 1] - invnormal(1 - (100-`level')/200) * `hm_se'[`i', 1]
                    matrix `hm_se_ci'[`i', 2] = `hm'[`i', 1] + invnormal(1 - (100-`level')/200) * `hm_se'[`i', 1]
                }
                matrix colnames `hm' = "Estimate"
                matrix colnames `hm_se' = "Std Err"
                matrix colnames `hm_se_ci' = "[`level'% Conf"  "Interval]"
                noi: matlist (`hm', `hm_se', `hm_se_ci'), bor(all) for(%9.3f) aligncolnames(r) rowtitle("")
            }
            else {
                noisily _coef_table, cionly level(`level') cformat("%9.3f")
            }
        }
        else {
            noisily ereturn display, level(`level') nopvalue cformat("%9.3f")
        }

        local add_space = 56
        if `use_thresholds' {
            local add_space = 58
        }
        noi display as text "Note:" _col(7) "Adjusted Multidimensional Headcount" _col(53) "M0" _col(`add_space') "= H*A"
        if `use_thresholds' {
            noi display as text _col(7) "Adjusted Poverty Gap" _col(53) "M1" _col(`add_space') "= H*A*G"
            if wordcount("`alpha'") > 2 {
                noi display as text _col(7) "Adjusted Foster-Greer-Thorbecke (FGT) Measure"_col(53) "M(a)" _col(`add_space') "= H*A*S(a)"
            }
            else {
                noi display as text _col(7) "Adjusted Foster-Greer-Thorbecke (FGT) Measure"_col(53) "M(2)" _col(`add_space') "= H*A*S(2)"
            }
        }

        if `warn_no_poor' == 1 {
            error _rc
            noi di as result "Note: No individual is multidimensionally deprived."
        }
        noi display as text ""

    /* ---------------------------------------------------------------------------

    Decomposition of M_alpha by Indicator

    --------------------------------------------------------------------------- */
        if "`nodecomposition'" == "" {
            if `ni' != 1 {
                // Variables for output
                tempname decomb_ind decomb_dom decomb_ind_V decomb_dom_V mat_transform

                // Variables for display
                tempname disp decomb_ind_disp decomb_dom_disp

                noi _mpi_domains `"`var_submit'"'

                matrix `decomb_ind' = r(mpi_decomposed)
                matrix `decomb_ind_V' = r(mpi_decomposed_V)

                local ncol = `=1 + wordcount("`alpha'")*`use_thresholds''

                mata: `A' = st_matrix("r(mpi_decomposed)")
                mata: cols = strtoreal(st_local("ncol"))
                mata: `A' = colshape(`A', cols)
                mata: st_matrix("`decomb_ind_disp'", `A')

                // Transform from the all-indicator view to the only-domain view
                matrix `mat_transform' = J(`=`ncol'*`ni'', `=`ncol'*`ndom'', 0)
                local idx = 0
                local count_ind = 0
                forvalues j = 1 / `ndom' {
                    forvalues i = 1 / `nind`j'' {
                        forvalues idx = 1 / `ncol' {
                            matrix `mat_transform'[`count_ind'*`ncol' + `idx', (`j'-1)*`ncol' + `idx'] = 1
                        }
                        local count_ind = `count_ind' + 1
                    }
                }

                matrix `decomb_dom' = `decomb_ind' * `mat_transform'
                matrix `decomb_dom_V' = `mat_transform'' * `decomb_ind_V' * `mat_transform'

                mata: `A' = colshape(st_matrix("`decomb_dom'"), cols)
                mata: st_matrix("`decomb_dom_disp'", `A')

                /*
                    Create the correct names.
                */
                if `use_thresholds' {
                    local idx = 0
                    local addlabel_alpha = ""
                    foreach a of local alpha {
                        local idx = `idx' + 1
                        local a_name = subinstr("`a'", ".", ",", .)
                        local addlabel_alpha = `"`addlabel_alpha' "M(`a_name')""'
                        local a_name = subinstr("`a'", ".", "p", .)
                        local a_name = subinstr("`a_name'", "-", "m", .)
                    }
                    local names = `""M0" `addlabel_alpha'"'
                }
                else {
                    local names = `""M0""'
                }
                matrix colnames `decomb_ind_disp' = `names'
                matrix colnames `decomb_dom_disp' = `names'

                local e_names = ""
                local lb_ind =  " "
                local lb_dom = " "
                forvalues j = 1 / `ndom' {
                    local lb_dom = `"`lb_dom' "domain `j'""'
                    local nind`j' = wordcount("`d`j''")

                    /*
                        Names for the domain-ereturn.
                    */
                    if `use_thresholds' {
                        local idx = 0
                        local addlabel_alpha = ""
                        foreach a of local alpha {
                            local idx = `idx' + 1
                            local a_name = subinstr("`a'", ".", ",", .)
                            local addlabel_alpha = `"`addlabel_alpha' "d_`j':M(`a_name')""'
                            local a_name = subinstr("`a'", ".", "p", .)
                            local a_name = subinstr("`a_name'", "-", "m", .)
                        }
                        local e_names = `"`e_names' "d_`j':M0" `addlabel_alpha'"'
                    }
                    else {
                        local e_names = `"`e_names' "d_`j':M0""'
                    }

                    /*
                        Iterate over indicators, create the to-be-estimated variables and the labels.
                    */
                    forvalues i = 1 / `nind`j'' {
                        local ind = word("`d`j''", `i')
                        local lb_ind = `"`lb_ind' "domain `j':`ind'""''
                    }
                }

                matrix rownames `decomb_ind_disp' = `lb_ind'
                matrix rownames `decomb_dom_disp' = `lb_dom'

                matrix colnames `decomb_dom' = `e_names'
                matrix colnames `decomb_dom_V' = `e_names'
                matrix rownames `decomb_dom_V' = `e_names'
            }

            /*
                Prepare for display and do consistency checks.
            */
            if `ni' != 1 {
                tempname check
                mata: `A' = colsum(st_matrix("`decomb_ind_disp'"))
                mata: st_matrix("`check'", `A')

                forvalues i = 1/`=colsof(`check')' {
                    cap assert abs(`check'[1, `i'] - 1) < `assertion_sensitivity' | ///
                               missing(`check'[1, `i']) | `warn_no_poor'
                    if _rc == 9 {
                        noi di as error "Decomposition by ind/dom failed in column `i'. (`=`check'[1, `i']' vs. 1)"
                    }
                    else {
                        error _rc
                    }
                }
                matrix `disp' = (`decomb_ind_disp' \ `check')
                matrix rownames `disp' = `lb_ind' Total
                noi matlist `disp', bor(all) for(%9.3f) lines(rowtotal) rowtitle(Indicator)  twidth(15)
                noi display in gr "Contribution of each indicator (%)"
                if `warn_no_poor' {
                    noi di as result "Note: No individual is multidimensionally deprived."
                }
                noi display as text ""
            }

            if (`ni' != `ndom') & (`ndom' != 1) {
                matrix `disp' = (`decomb_dom_disp' \ `check')
                matrix rownames `disp' = `lb_dom' Total
                noi matlist `disp', bor(all) for(%9.3f) lines(rowtotal) rowtitle(Domain)
                noi display in gr "Contribution of each domain (%)"
                if `warn_no_poor' {
                    noi di as result "Note: No individual is multidimensionally deprived."
                }
                noi display as text ""
            }
        }

    /* ---------------------------------------------------------------------------

        Decomposition of M_alpha by subgroups

    --------------------------------------------------------------------------- */

        if ("`by'" != "")   {

            unab by : `by'

            noi: display as result _n _n "Decomposition by subgroups"
            noi: display as result "{hline 32}"

            noi display in gr "MPI by: " in ye "`by'"

            // Display a warning message if single groups have no deprived individuals.
            local warn_no_poor = 0
            local warn_no_poor_groups = ""

            levelsof `by' `l_if', local(by_levels)
            foreach by_level of local by_levels {
                cap assert `_H' == 0 if `by' == `by_level'
                if _rc == 0 {
                    local warn_no_poor = `warn_no_poor' + 1
                    local warn_no_poor_groups = "`warn_no_poor_groups' `by_level'"
                }
            }
            foreach by_level of local by_levels {
                cap assert `_H' == 0 if `by' == "`by_level'"
                if _rc == 0 {
                    local warn_no_poor = `warn_no_poor' + 1
                    local warn_no_poor_groups = `"`warn_no_poor_groups' "`by_level'""'
                }
            }
            if `warn_no_poor' > 0 {
                if `warn_no_poor' > 1 {
                    local s = "s"
                }
                else {
                    local s = ""
                }
                noi display as result `"Note: No deprived individuals in by-group`s' `warn_no_poor_groups'."'
                local warn_no_poor = 1
            }

            // Get estimates
            noi _mpi_main `"`var_submit' by=`by'"'

            // Retreive the results
            local levelnames = r(over_namelist)

            // Variables for storange and returning
            tempname hm_by hm_by_V hm_by_perc hm_by_perc_V

            // Variables for display
            tempname hm_by_disp hm_by_perc_disp

            matrix `hm_by' = r(mpi_main)
            matrix `hm_by_V' = r(mpi_main_V)
            matrix `hm_by_perc' = r(mpi_perc)
            matrix `hm_by_perc_V' = r(mpi_perc_V)

            local ncol = `=wordcount("`levelnames'")'
            mata: `A' = st_matrix("`hm_by'")
            mata: cols = strtoreal(st_local("ncol"))
            mata: `A' = colshape(`A', cols)
            mata: st_matrix("`hm_by_disp'", `A')

            local ncol = `=wordcount("`levelnames'")'
            mata: `A' = st_matrix("`hm_by_perc'")
            mata: cols = strtoreal(st_local("ncol"))
            mata: `A' = colshape(`A', cols)
            mata: st_matrix("`hm_by_perc_disp'", `A')

            // Prepare subpopulation labels
            local levels ""
            // Check if the by-variable has a label.
            local by_lab: value label `by'
            if "`by_lab'" != "" {
                foreach i of local levelnames {
                    local label: label `by_lab' `i'
                    local levels = `"`levels' ":`label'""''
                }
            }
            else {
                foreach i of local levelnames {
                    local levels = `"`levels' ":`by'_`i'""''
                }
            }

            /*
                Apply the labels (from above)
            */
            if `use_thresholds' {
                local addlabel_alpha = ""
                foreach a of local alpha {
                    // Malpha
                    local a_name = subinstr("`a'", ".", "p", .)
                    local a_name = subinstr("`a_name'", "-", "m", .)
                    local a_name = subinstr("`a'", ".", ",", .)
                    local addlabel_alpha = `"`addlabel_alpha' "M(`a_name')""'
                }
            }


            // Printing the absolute indices by subgroups

            // Get subpopulation count for all _valid_ subgroups.
            // The correction is only necessary in the svy-case, because
            // proportion does not exclude observations when the other
            // covariates are missing by default.
            count if `touse' == 0
            local use_touse = (r(N) != 0 & "`svy'" != "")
            if `use_touse' {
                `svy' proportion `by' `weight_exp' `l_if', over(`touse')
            }
            else {
                `svy' proportion `by' `weight_exp' `l_if'
            }

            // Cut missing values from the result vector, if they exist. This is only
            // relevant for the svy case.
            tempname proportions
            mat `proportions' = e(b)

            // Remove results which are missing (exactly =0)
            // and such that were from the touse==0 scenario.
            local i = 1
            local real_i = 1
            local max = colsof(`proportions')
            while `i' < `max' {
                if `use_touse' {
                    local condition = ((`proportions'[1, `i'] == 0) | (mod(`real_i', 2)==1))
                }
                else {
                    local condition = (`proportions'[1, `i'] == 0)
                }
                if `condition' {
                    if `i' != 1 & `i' != `max' {
                        matrix `proportions' = `proportions'[1, 1..`=`i'-1'], `proportions'[1, `=`i'+1'..`max']
                    }
                    else if `i' == 1 {
                        matrix `proportions' = `proportions'[1, 2..`max']
                    }
                    else {
                        matrix `proportions' = `proportions'[1, 1..`=`max'-1']
                    }
                    local i = `i'-1
                    local max = `max'-1
                }
                local i = `i'+1
                local real_i = `real_i'+1
            }

            matrix `hm_by_disp' = `hm_by_disp' \ (`proportions')
            matrix `hm_by_disp' = `hm_by_disp' , ( (`hm_1')' \ 1 )

            matrix colnames `hm_by_disp' = `levels' "Total"
            if `use_thresholds' {
                matrix rownames `hm_by_disp' = "H" "M0" `addlabel_alpha' "pop share"
            }
            else {
                matrix rownames `hm_by_disp' = "H" "M0" `addlabel_alpha' "pop share"
            }
            */

            // Checking that the decomposition was correct
            local ncol = colsof(`hm_by_disp')
            local nrow = rowsof(`hm_by_disp')
            forvalues col = 1/`ncol' {
                local total = 0
                forvalues row = 1/`nrow' {
                    // Values * Population share
                    if !missing(`hm_by_disp'[`row', `col']) {
                        local total = `total' + `hm_by_disp'[`row', `col']*`hm_by_disp'[`row', `ncol'+1]
                    }
                    else {
                        local total = .
                    }
                }
                cap assert abs(`total' - `hm_by_disp'[`nrow'+1, `col']) < `assertion_sensitivity' | ///
                           missing(`total')
                if _rc == 9 {
                    noi display as error "By-Decomposition failed in col `col'."
                }
                else {
                    error _rc
                }
            }

            noi matlist `hm_by_disp', bor(all) for(%9.3f) lines(rctotal)
            noi display as text "Indices by subgroup (absolute)"
            if `warn_no_poor' {
                noi di as result "Note: In at least one category, no individual is multidimensionally deprived."
            }
            noi display as text ""

            // Printing the proportional contribution of subgroups to the indices
            tempname check disp
            mata: `A' = rowsum(st_matrix("`hm_by_perc_disp'"))
            mata: st_matrix("`check'", `A')
            matrix `disp' = `hm_by_perc_disp' , `check'

            local ncol = colsof(`disp')
            local nrow = rowsof(`disp')
            forvalues row = 1/`nrow' {
                forvalues col = 1/`ncol' {
                    // If the value is zero, (if there are no
                    // poor people in the subgroup), ignore it.
                    if `disp'[`row', `col'] == 0 {
                        matrix `disp'[`row', `col'] = .
                    }
                }
                cap assert abs(`disp'[`row', `ncol'] - 1) < `assertion_sensitivity' | missing(`disp'[`row', `ncol'])
                if _rc == 9 {
                    noi di as error "Computation of row sums in proportions decomposition failed. `rowsum' vs 1 in row `row'."
                }
                else {
                    error _rc
                }
            }

            matrix colnames `disp' = `levels' "Total"
            if `use_thresholds' {
                matrix rownames `disp' = "H" "M0" `addlabel_alpha'
            }
            else {
                matrix rownames `disp' = "H" "M0"
            }

            noi matlist `disp', bor(all) for(%9.3f) lines(coltotal)
            noi display as text "Contribution of subgroups to indices (%)"
            if `warn_no_poor' {
                noi di as result "Note: In at least one category, no individual is multidimensionally deprived."
            }
            noi display as text ""

    /* ---------------------------------------------------------------------------

    Decomposition of M_alpha by Subgroup AND Indicator

    --------------------------------------------------------------------------- */

            if "`nodecomposition'" == "" {
                if (`ni' != 1) {
                    noi _mpi_domains `"`var_submit' by=`by'"'

                    // Retreive the results
                    local levelnames = r(over_namelist)

                    // Variables for returning.
                    tempname hm_by_ind hm_by_dom hm_by_ind_V hm_by_dom_V

                    // Variables for display.
                    tempname hm_by_ind_disp hm_by_dom_disp

                    matrix `hm_by_ind' = r(mpi_decomposed)
                    matrix `hm_by_ind_V' = r(mpi_decomposed_V)

                    local ncol = `=wordcount("`levelnames'")'
                    mata: `A' = st_matrix("`hm_by_ind'")
                    mata: cols = strtoreal(st_local("ncol"))
                    mata: `A' = colshape(`A', cols)
                    mata: st_matrix("`hm_by_ind_disp'", `A')

                    // Transform from the all-indicator view to the only-domain view
                    local n_alpha = wordcount("`alpha'")*`use_thresholds' + 1
                    matrix `mat_transform' = J(`=`ncol'*`ni'*`n_alpha'', `=`ncol'*`ndom'*`n_alpha'', 0)
                    local count_ind = 0
                    forvalues j = 1 / `ndom' {
                        forvalues i = 1 / `nind`j'' {
                            forvalues idx = 1 / `=`ncol'*`n_alpha'' {
                                matrix `mat_transform'[`count_ind'*`ncol'*`n_alpha' + `idx', (`j'-1)*`ncol'*`n_alpha' + `idx'] = 1
                            }
                            local count_ind = `count_ind' + 1
                        }
                    }

                    matrix `hm_by_dom' = `hm_by_ind' * `mat_transform'
                    matrix `hm_by_dom_V' = `mat_transform'' * `hm_by_ind_V' * `mat_transform'

                    mata: `A' = colshape(st_matrix("`hm_by_dom'"), cols)
                    mata: st_matrix("`hm_by_dom_disp'", `A')

                    // The resulting matrices need to be brought to the final form.
                    // Number of columns = Number of levels of `by'
                    local ind_names = ""
                    local dom_names = ""
                    local e_dom_names = ""
                    forvalues j = 1/`ndom' {
                        forvalues i = 1/`nind`j'' {
                            local ind = word("`d`j''", `i')
                            local ind_names = `"`ind_names' "M0:`ind'""'
                            if `use_thresholds' {
                                foreach a of local alpha {
                                    local a_idx = `a_idx' + 1
                                    local a_name = subinstr("`a'", ".", ",", .)
                                    local ind_names = `"`ind_names' "M(`a_name'):`ind'""'
                                }
                            }
                        }
                        local dom_names = `"`dom_names' "M0:domain `j'""'
                        foreach i of local levelnames {
                            local e_dom_names = `"`e_dom_names' "`by'_`i':d`j'_M0""'
                        }

                        if `use_thresholds' {
                            foreach a of local alpha {
                                local a_idx = `a_idx' + 1
                                local a_name = subinstr("`a'", ".", ",", .)
                                local dom_names = `"`dom_names' "M(`a_name'):domain `j'""'
                                foreach i of local levelnames {
                                    local e_dom_names = `"`e_dom_names' "`by'_`i':d`j'_M(`a_name')""'
                                }
                            }
                        }
                    }
                    matrix rownames `hm_by_ind_disp' = `ind_names'
                    matrix colnames `hm_by_ind_disp' = `levels'

                    matrix colnames `hm_by_dom' = `e_dom_names'
                    matrix colnames `hm_by_dom_V' = `e_dom_names'
                    matrix rownames `hm_by_dom_V' = `e_dom_names'

                    matrix rownames `hm_by_dom_disp' = `dom_names'
                    matrix colnames `hm_by_dom_disp' = `levels'

                    // The domains are ordered as d1:M0 d1:M1 d2:M0 d2:M1 ...
                    // Reorder them as d1:M0 d2:M0 d1:M1 d2:M1
                    tempname hm_by_ind_new hm_by_dom_new
                    * tempname hm_by_ind_se_new hm_by_dom_se_new
                    local idx_ind = 0
                    if `use_thresholds' {
                        local nindices = wordcount("0 `alpha'")
                    }
                    else {
                        local nindices = 1
                    }
                    forvalues index = 1/`nindices' {
                        local n_prev_dom = 0
                        forvalues j = 1/`ndom' {
                            matrix `hm_by_dom_new' = nullmat(`hm_by_dom_new') \ `hm_by_dom_disp'[(`j'-1)*`nindices'+`index', 1...]
                            * matrix `hm_by_dom_se_new' = nullmat(`hm_by_dom_se_new') \ `hm_by_dom_se'[(`j'-1)*`nindices'+`index', 1...]
                            local nind`j' = wordcount("`d`j''")
                            forvalues i = 1/`nind`j'' {
                                // Do the same for the indicators
                                matrix `hm_by_ind_new' = nullmat(`hm_by_ind_new') \ `hm_by_ind_disp'[(`i'-1+`n_prev_dom')*`nindices'+`index', 1...]
                                * matrix `hm_by_ind_se_new' = nullmat(`hm_by_ind_se_new') \ `hm_by_ind'[(`i'-1+`n_prev_dom')*`nindices'+`index', 1...]
                                local n_ind = `n_ind' + 1
                            }
                            local n_prev_dom = `n_prev_dom' + `nind`j''
                        }
                    }
                    matrix `hm_by_ind_disp' = `hm_by_ind_new'
                    matrix `hm_by_dom_disp' = `hm_by_dom_new'

                    // Prepare the matrix of totals and their names to add them at
                    // the proper places in the display matrix
                    tempname total c
                    matrix `total' = J(1, colsof(`hm_by_ind_disp'), .)
                    matrix rownames `total' = "Total"
                    local total_names = `" "M0:Total" "'
                    foreach a of local alpha {
                        local a_idx = `a_idx' + 1
                        local a_name = subinstr("`a'", ".", ",", .)
                        local total_names = `"`total_names' "M(`a_name'):Total""'
                    }

                    // Select the individual indicators to reorganize the display matrix
                    // and add the 'total' column
                    cap matrix drop `disp'
                    local cspec = ""
                    local rspec = ""
                    forvalues index = 1/`nindices' {
                        matrix rownames `total' = `=word(`"`total_names'"', `index')'
                        local to = (`index')*`ni'
                        local to = cond(`to' > rowsof(`hm_by_ind_disp'), ., `to')
                        local from = (`index'-1)*`ni'+1
                        matrix `c' = `decomb_ind_disp'[1..`ni', `index']
                        matrix colnames `c' = "Total"
                        matrix `disp' = nullmat(`disp') \ ((`hm_by_ind_disp'[`from'..`to', 1...], `c') \ (`total', 1))
                        local rspec = "`rspec'-"
                        forvalues i = 1/`=`ni'-1' {
                            local rspec = "`rspec'&"
                        }
                        local rspec = "`rspec'-"
                        // Check that the total is really 1, by levels
                        forvalues j = 1/`=colsof(`hm_by_ind_disp')' {
                            local total_check = 0
                            local count_zero = 0    // Count the number of exact zeros.
                            forvalues i = `from'/`to' {
                                if `hm_by_ind_disp'[`i', `j'] == 0 {
                                    local count_zero = `count_zero' + 1
                                }
                                else {
                                    local total_check = `total_check' + `hm_by_ind_disp'[`i', `j']
                                }
                            }
                            // If the number of zero-entries is equal to the number of entries,
                            // the by-groups is completely missing. Make this visible in the output matrix.
                            if `count_zero' == `to'-`from'+1 {
                                local total_check = .
                                forvalues i = 1/`count_zero' {
                                    matrix `disp'[`=rowsof(`disp')-`i'', `j'] = .
                                }
                            }
                            cap assert abs(`total_check' - 1) < `assertion_sensitivity' | ///
                                       missing(`total_check')
                            if _rc == 9 {
                                noi di as error "Decomposition by by+ind failed, index=`index', j=`j'. (`total_check' vs. 1)"
                            }
                            else {
                                error _rc
                                matrix `disp'[`=rowsof(`disp')', `j'] = `total_check'
                            }
                            if `index' == 1 {
                                local cspec = "`cspec' %9.3f"
                                if `j' != `=colsof(`hm_by_ind_disp')' {
                                    local cspec = "`cspec' &"
                                }
                            }
                        }
                    }
                    noi matlist `disp', rowtitle("") cspec(| b %16s | `cspec' | %9.3f |) rspec(-`rspec'-)
                    noi display as text "Contribution of each indicator (%)"
                    if `warn_no_poor' {
                        noi di as result "Note: In at least one category, no individual is multidimensionally deprived."
                    }
                    noi display as text ""

                    // Displaying Domain decompositions
                    if (`ni' != `ndom')  & (`ndom' != 1) {
                        matrix drop `disp'
                        local cspec = ""
                        local rspec = ""
                        forvalues index = 1/`nindices' {
                            matrix rownames `total' = `=word(`"`total_names'"', `index')'
                            local to = ((`index')*`ndom')
                            local to = cond(`to' > rowsof(`hm_by_dom_disp'), ., `to')
                            local from = ((`index'-1)*`ndom'+1)
                            matrix `c' = `decomb_dom_disp'[1..`ndom', `index']
                            matrix colnames `c' = "Total"
                            matrix `disp' = nullmat(`disp') \ ((`hm_by_dom_disp'[`from'..`to', 1...], `c') \ (`total', 1))
                            local rspec = "`rspec'-"
                            forvalues i = 1/`=`ndom'-1' {
                                local rspec = "`rspec'&"
                            }
                            local rspec = "`rspec'-"
                            // Check that the total is really 1
                            forvalues j = 1/`=colsof(`hm_by_dom_disp')' {
                                local total_check = 0
                                local count_zero = 0    // Count the number of exact zeros.
                                forvalues i = `from'/`to' {
                                    if `hm_by_dom_disp'[`i', `j'] == 0 {
                                        local count_zero = `count_zero' + 1
                                    }
                                    else {
                                        local total_check = `total_check' + `hm_by_dom_disp'[`i', `j']
                                    }
                                }
                                // If the number of zero-entries is equal to the number of entries,
                                // the by-groups is completely missing. Make this visible in the output matrix.
                                if `count_zero' == `to'-`from'+1 {
                                    local total_check = .
                                    forvalues i = 1/`count_zero' {
                                        matrix `disp'[`=rowsof(`disp')-`i'', `j'] = .
                                    }
                                }
                                cap ass
                                cap assert abs(`total_check' - 1) < `assertion_sensitivity' | ///
                                           missing(`total_check')
                                if _rc == 9 {
                                    noi di as error "Decomposition by by+dom failed, index=`index', j=`j'. (`total_check' vs. 1)"
                                }
                                else {
                                    error _rc
                                    matrix `disp'[`=rowsof(`disp')', `j'] = `total_check'
                                }
                                if `index' == 1 {
                                    local cspec = "`cspec' %9.3f"
                                    if `j' != `=colsof(`hm_by_dom_disp')' {
                                        local cspec = "`cspec' &"
                                    }
                                }
                            }
                        }
                        noi matlist `disp', rowtitle("") cspec(| b %16s | `cspec' | %9.3f |) rspec(-`rspec'-)
                        noi display as text "Contribution of each domain (%)"
                        if `warn_no_poor' {
                            noi di as result "Note: In at least one category, no individual is multidimensionally deprived."
                        }
                        noi display as text ""
                    }
                }
            } // Should the decomposition by domains/indicators be computed ?
        } // End decomposition (by)

        ereturn clear

        // Posting the correct values in e(b) and e(V).
        // At the begin of mpi, there is a check that at most one
        // post statement is present.
        tempname mat_keep_1 mat_keep_2
        if "`postmain'" != "" {
            matrix `mat_keep_1' = `hm_1'
            matrix `mat_keep_2' = `hm_V_1'
            ereturn post `hm_1' `hm_V_1' `weight_exp', esample(`touse') obs(`N')
            matrix `hm_1' = `mat_keep_1'
            matrix `hm_V_1' = `mat_keep_2'
            ereturn local type_b_return = "Main results."
        }

        if "`postadditional'" != "" {
            matrix `mat_keep_1' = `hm_2'
            matrix `mat_keep_2' = `hm_V_2'
            ereturn post `hm_2' `hm_V_2' `weight_exp', esample(`touse') obs(`N')
            matrix `hm_2' = `mat_keep_1'
            matrix `hm_V_2' = `mat_keep_2'
            ereturn local type_b_return = "Additional results."
        }

        if "`postindicators'" != "" {
            if "`nodecomposition'" == "" & `ni' != 1 {
                matrix `mat_keep_1' = `decomb_ind'
                matrix `mat_keep_2' = `decomb_ind_V'
                ereturn post `decomb_ind' `decomb_ind_V' `weight_exp', esample(`touse') obs(`N')
                matrix `decomb_ind' = `mat_keep_1'
                matrix `decomb_ind_V' = `mat_keep_2'
                ereturn local type_b_return = "Main results decomposed by indicators."
            }
            else {
                noi display as error "Error: postindicators not possible in the context."
                error 198
            }
        }

        if "`postdomains'" != "" {
            if "`nodecomposition'" == "" & (`ni' != `ndom')  & (`ndom' != 1) {
                matrix `mat_keep_1' = `decomb_dom'
                matrix `mat_keep_2' = `decomb_dom_V'
                ereturn post `decomb_dom' `decomb_dom_V' `weight_exp', esample(`touse') obs(`N')
                matrix `decomb_dom' = `mat_keep_1'
                matrix `decomb_dom_V' = `mat_keep_2'
                ereturn local type_b_return = "Main results decomposed by domains."
            }
            else {
                noi display as error "Error: postdomains not possible in the context."
                error 198
            }
        }

        if "`postbymain'" != "" {
            if "`by'" != "" {
                matrix `mat_keep_1' = `hm_by'
                matrix `mat_keep_2' = `hm_by_V'
                ereturn post `hm_by' `hm_by_V' `weight_exp', esample(`touse') obs(`N')
                matrix `hm_by' = `mat_keep_1'
                matrix `hm_by_V' = `mat_keep_2'
                ereturn local type_b_return = "Main results decomposed by by-levels."
            }
            else {
                noi display as error "Error: postbymain not possible in the context. Specify by."
                error 198
            }
        }

        if "`postbyproportion'" != "" {
            if "`by'" != "" {
                matrix `mat_keep_1' = `hm_by_perc'
                matrix `mat_keep_2' = `hm_by_perc_V'
                ereturn post `hm_by_perc' `hm_by_perc_V' `weight_exp', esample(`touse') obs(`N')
                matrix `hm_by_perc' = `mat_keep_1'
                matrix `hm_by_perc_V' = `mat_keep_2'
                ereturn local type_b_return = "by proportion"
            }
            else {
                noi display as error "Error: postbyproportion not possible in the context. Specify by."
                error 198
            }
        }

        if "`postbyindicators'" != "" {
            if "`by'" != "" & "`nodecomposition'" == "" & `ni' != 1 {
                matrix `mat_keep_1' = `hm_by_ind'
                matrix `mat_keep_2' = `hm_by_ind_V'
                ereturn post `hm_by_ind' `hm_by_ind_V' `weight_exp', esample(`touse') obs(`N')
                matrix `hm_by_ind' = `mat_keep_1'
                matrix `hm_by_ind_V' = `mat_keep_2'
                ereturn local type_b_return = "Main results decomposed by indicators and by-levels."
            }
            else {
                noi display as error "Error: postbyindicators not possible in the context."
                error 198
            }
        }

        if "`postbydomains'" != "" {
            if "`by'" != "" & "`nodecomposition'" == "" & (`ni' != `ndom') & (`ndom' != 1) {
                matrix `mat_keep_1' = `hm_by_dom'
                matrix `mat_keep_2' = `hm_by_dom_V'
                ereturn post `hm_by_dom' `hm_by_dom_V' `weight_exp', esample(`touse') obs(`N')
                matrix `hm_by_dom' = `mat_keep_1'
                matrix `hm_by_dom_V' = `mat_keep_2'
                ereturn local type_b_return = "Main results decomposed by domains and by-levels."
            }
            else {
                noi display as error "Error: postbydomains not possible in the context."
                error 198
            }
        }

        // Return by-results
        if "`by'" != "" {
            if "`nodecomposition'" == "" {
                if `ni' != 1 {
                    ereturn matrix by_ind_V = `hm_by_ind_V'
                    ereturn matrix by_ind = `hm_by_ind'
                }
                if (`ni' != `ndom') & (`ndom' != 1) {
                    ereturn matrix by_dom_V = `hm_by_dom_V'
                    ereturn matrix by_dom = `hm_by_dom'
                }
            }
            ereturn matrix by_mpi_V = `hm_by_V'
            ereturn matrix by_mpi = `hm_by'
            ereturn matrix by_mpi_pc_V = `hm_by_perc_V'
            ereturn matrix by_mpi_pc = `hm_by_perc'
        }

        if "`nodecomposition'" == "" {
            if (`ni' != `ndom')  & (`ndom' != 1) {
                ereturn matrix dom_V = `decomb_dom_V'
                ereturn matrix dom = `decomb_dom'
            }
            if `ni' != 1 {
                ereturn matrix ind_V = `decomb_ind_V'
                ereturn matrix ind = `decomb_ind'
            }
        }

        ereturn matrix mpi_add_V = `hm_V_2'
        ereturn matrix mpi_add = `hm_2'
        ereturn matrix mpi_main_V = `hm_V_1'
        ereturn matrix mpi_main = `hm_1'

        ereturn local cmd = "mpi"
        ereturn local cmdline = `"mpi `*'"'

        return clear

        noi display in smcl "Type {stata ereturn list:ereturn list} to see the list of saved results and more" _n ///
                            "information on the estimation sample."
    }

end