*-------------------------------------------------------------------------------
* Command: "opl_tb_cba" = "OPL cost benefit analysis threshold based"
*-------------------------------------------------------------------------------
*! version 1.1.0 14jul2026
program define opl_tb_cba, eclass
    version 18.0

    syntax varname(numeric) [if] [in],                       ///
        COST(varname numeric)                               ///
        SELECT(varlist numeric min=1 max=2)                 ///
        [ LAMBDA(real 1)                                    ///
          NPOINTS(integer 101)                              ///
          CUSTOM(numlist)                                   ///
          CUSTOMPOLICY(varname numeric)                     ///
          GENerate(name)                                    ///
          GRAPH                                             ///
          REPLACE ]

    marksample touse
    markout `touse' `varlist' `cost' `select'

    if "`custompolicy'" != "" {
        markout `touse' `custompolicy'
    }

    if `lambda' < 0 {
        di as error "lambda() must be nonnegative"
        exit 198
    }

    if `npoints' < 2 {
        di as error "npoints() must be at least 2"
        exit 198
    }

    local tau `varlist'
    local nselect : word count `select'
    local x1 : word 1 of `select'

    if `nselect' == 2 {
        local x2 : word 2 of `select'
    }

    quietly count if `touse'
    local N = r(N)

    if `N' == 0 {
        di as error "no observations"
        exit 2000
    }

    /*
    Validate user-supplied customized policy.

    The policy must be coded:
        1 = treated
        0 = not treated
    */
    if "`custompolicy'" != "" {
        quietly count if `touse' & !inlist(`custompolicy', 0, 1)

        if r(N) > 0 {
            di as error ///
                "custompolicy() must contain only 0 and 1 in the estimation sample"
            exit 459
        }
    }

    tempvar net z1 z2 pfb popt pcand
    quietly gen double `net' = ///
        `tau' - `lambda' * `cost' if `touse'

    /*
    Standardize selection variables to [0,1]
    */
    quietly summarize `x1' if `touse', meanonly

    local x1min = r(min)
    local x1max = r(max)

    if (`x1max' - `x1min') <= 0 {
        di as error ///
            "selection variable `x1' is constant in the estimation sample"
        exit 498
    }

    quietly gen double `z1' = ///
        (`x1' - `x1min') / (`x1max' - `x1min') if `touse'

    if `nselect' == 2 {

        quietly summarize `x2' if `touse', meanonly

        local x2min = r(min)
        local x2max = r(max)

        if (`x2max' - `x2min') <= 0 {
            di as error ///
                "selection variable `x2' is constant in the estimation sample"
            exit 498
        }

        quietly gen double `z2' = ///
            (`x2' - `x2min') / (`x2max' - `x2min') if `touse'
    }

    /*
    First-best policy:
    treat individual i iff

        tau_i - lambda * cost_i > 0
    */
    quietly gen byte `pfb' = (`net' > 0) if `touse'

    quietly opl_tb_cba__eval,                         ///
        policy(`pfb')                                     ///
        tau(`tau')                                        ///
        cost(`cost')                                      ///
        lambda(`lambda')                                  ///
        touse(`touse')

    local N_fb        = r(N_treated)
    local coverage_fb = r(coverage)
    local TB_fb       = r(total_benefit)
    local TC_fb       = r(total_cost)
    local Q_fb        = r(Q)
    local ATET_fb     = r(ATET)
    local avgcost_fb  = r(avg_cost)
    local avgsurp_fb  = r(avg_surplus)

    /*
    Initialize constrained optimum at treat-none policy.
    */
    quietly gen byte `popt'  = 0 if `touse'
    quietly gen byte `pcand' = 0 if `touse'

    local bestQ  = 0
    local bestg1 = .
    local bestg2 = .

    /*
    Search over threshold policies.
    */
    tempname GRID

    if `nselect' == 1 {

        matrix `GRID' = J(`npoints', 9, .)

        matrix colnames `GRID' =                       ///
            threshold1                                ///
            N_treated                                 ///
            coverage                                  ///
            total_benefit                             ///
            total_cost                                ///
            Q                                         ///
            ATET                                      ///
            avg_cost                                  ///
            avg_surplus

        local row = 0

        forvalues m = 0/`=`npoints' - 1' {

            local ++row
            local g1 = `m' / (`npoints' - 1)

            quietly replace `pcand' = ///
                (`z1' >= `g1') if `touse'

            quietly opl_tb_cba__eval,             ///
                policy(`pcand')                        ///
                tau(`tau')                            ///
                cost(`cost')                          ///
                lambda(`lambda')                      ///
                touse(`touse')

            matrix `GRID'[`row',1] = `g1'
            matrix `GRID'[`row',2] = r(N_treated)
            matrix `GRID'[`row',3] = r(coverage)
            matrix `GRID'[`row',4] = r(total_benefit)
            matrix `GRID'[`row',5] = r(total_cost)
            matrix `GRID'[`row',6] = r(Q)
            matrix `GRID'[`row',7] = r(ATET)
            matrix `GRID'[`row',8] = r(avg_cost)
            matrix `GRID'[`row',9] = r(avg_surplus)

            if r(Q) > `bestQ' {

                local bestQ  = r(Q)
                local bestg1 = `g1'

                quietly replace `popt' = `pcand' if `touse'
            }
        }
    }
    else {

        local npolicies = `npoints' * `npoints'

        matrix `GRID' = J(`npolicies', 10, .)

        matrix colnames `GRID' =                       ///
            threshold1                                ///
            threshold2                                ///
            N_treated                                 ///
            coverage                                  ///
            total_benefit                             ///
            total_cost                                ///
            Q                                         ///
            ATET                                      ///
            avg_cost                                  ///
            avg_surplus

        local row = 0

        forvalues m1 = 0/`=`npoints' - 1' {

            local g1 = `m1' / (`npoints' - 1)

            forvalues m2 = 0/`=`npoints' - 1' {

                local ++row
                local g2 = `m2' / (`npoints' - 1)

                quietly replace `pcand' = ///
                    (`z1' >= `g1' & `z2' >= `g2') if `touse'

                quietly opl_tb_cba__eval,         ///
                    policy(`pcand')                    ///
                    tau(`tau')                        ///
                    cost(`cost')                      ///
                    lambda(`lambda')                  ///
                    touse(`touse')

                matrix `GRID'[`row',1]  = `g1'
                matrix `GRID'[`row',2]  = `g2'
                matrix `GRID'[`row',3]  = r(N_treated)
                matrix `GRID'[`row',4]  = r(coverage)
                matrix `GRID'[`row',5]  = r(total_benefit)
                matrix `GRID'[`row',6]  = r(total_cost)
                matrix `GRID'[`row',7]  = r(Q)
                matrix `GRID'[`row',8]  = r(ATET)
                matrix `GRID'[`row',9]  = r(avg_cost)
                matrix `GRID'[`row',10] = r(avg_surplus)

                if r(Q) > `bestQ' {

                    local bestQ  = r(Q)
                    local bestg1 = `g1'
                    local bestg2 = `g2'

                    quietly replace `popt' = `pcand' if `touse'
                }
            }
        }
    }

    /*
    Evaluate optimal constrained threshold policy.
    */
    quietly opl_tb_cba__eval,                         ///
        policy(`popt')                                    ///
        tau(`tau')                                        ///
        cost(`cost')                                      ///
        lambda(`lambda')                                  ///
        touse(`touse')

    local N_opt        = r(N_treated)
    local coverage_opt = r(coverage)
    local TB_opt       = r(total_benefit)
    local TC_opt       = r(total_cost)
    local Q_opt        = r(Q)
    local ATET_opt     = r(ATET)
    local avgcost_opt  = r(avg_cost)
    local avgsurp_opt  = r(avg_surplus)

    /*
    Convert optimal standardized thresholds back to original scale.
    */
    local bestg1_lev = .
    local bestg2_lev = .

    if `bestg1' < . {
        local bestg1_lev = ///
            `x1min' + `bestg1' * (`x1max' - `x1min')
    }

    if `nselect' == 2 & `bestg2' < . {
        local bestg2_lev = ///
            `x2min' + `bestg2' * (`x2max' - `x2min')
    }

    /*
    Evaluate user-supplied customized policy.
    */
    local has_custompolicy = 0

    tempname CUSTOMPOLICYMAT

    if "`custompolicy'" != "" {

        local has_custompolicy = 1

        quietly opl_tb_cba__eval,                     ///
            policy(`custompolicy')                         ///
            tau(`tau')                                    ///
            cost(`cost')                                  ///
            lambda(`lambda')                              ///
            touse(`touse')

        local N_custompolicy        = r(N_treated)
        local coverage_custompolicy = r(coverage)
        local TB_custompolicy       = r(total_benefit)
        local TC_custompolicy       = r(total_cost)
        local Q_custompolicy        = r(Q)
        local ATET_custompolicy     = r(ATET)
        local avgcost_custompolicy  = r(avg_cost)
        local avgsurp_custompolicy  = r(avg_surplus)

        matrix `CUSTOMPOLICYMAT' = J(1, 8, .)

        matrix colnames `CUSTOMPOLICYMAT' =             ///
            N_treated                                   ///
            coverage                                    ///
            total_benefit                               ///
            total_cost                                  ///
            Q                                           ///
            ATET                                        ///
            avg_cost                                    ///
            avg_surplus

        matrix `CUSTOMPOLICYMAT'[1,1] = `N_custompolicy'
        matrix `CUSTOMPOLICYMAT'[1,2] = `coverage_custompolicy'
        matrix `CUSTOMPOLICYMAT'[1,3] = `TB_custompolicy'
        matrix `CUSTOMPOLICYMAT'[1,4] = `TC_custompolicy'
        matrix `CUSTOMPOLICYMAT'[1,5] = `Q_custompolicy'
        matrix `CUSTOMPOLICYMAT'[1,6] = `ATET_custompolicy'
        matrix `CUSTOMPOLICYMAT'[1,7] = `avgcost_custompolicy'
        matrix `CUSTOMPOLICYMAT'[1,8] = `avgsurp_custompolicy'
    }

    /*
    Custom threshold policy or policies specified through custom().
    */
    tempname CUSTOMMAT
    local ncustom = 0

    if "`custom'" != "" {

        local ncustom : word count `custom'

        if `nselect' == 2 & `ncustom' != 2 {
            di as error ///
                "with two selection variables, custom() must contain exactly two thresholds"
            exit 198
        }

        foreach g of numlist `custom' {

            if (`g' < 0 | `g' > 1) {
                di as error ///
                    "all custom thresholds must lie in [0,1]"
                exit 198
            }
        }

        if `nselect' == 1 {

            matrix `CUSTOMMAT' = J(`ncustom', 9, .)

            matrix colnames `CUSTOMMAT' =                ///
                threshold1                               ///
                N_treated                                ///
                coverage                                 ///
                total_benefit                            ///
                total_cost                               ///
                Q                                        ///
                ATET                                     ///
                avg_cost                                 ///
                avg_surplus

            local crow = 0

            foreach g1 of numlist `custom' {

                local ++crow

                quietly replace `pcand' = ///
                    (`z1' >= `g1') if `touse'

                quietly opl_tb_cba__eval,           ///
                    policy(`pcand')                      ///
                    tau(`tau')                          ///
                    cost(`cost')                        ///
                    lambda(`lambda')                    ///
                    touse(`touse')

                matrix `CUSTOMMAT'[`crow',1] = `g1'
                matrix `CUSTOMMAT'[`crow',2] = r(N_treated)
                matrix `CUSTOMMAT'[`crow',3] = r(coverage)
                matrix `CUSTOMMAT'[`crow',4] = r(total_benefit)
                matrix `CUSTOMMAT'[`crow',5] = r(total_cost)
                matrix `CUSTOMMAT'[`crow',6] = r(Q)
                matrix `CUSTOMMAT'[`crow',7] = r(ATET)
                matrix `CUSTOMMAT'[`crow',8] = r(avg_cost)
                matrix `CUSTOMMAT'[`crow',9] = r(avg_surplus)

                if "`generate'" != "" {

                    local cname = "`generate'_custom`crow'"

                    capture confirm new variable `cname'

                    if _rc {

                        if "`replace'" != "" {
                            capture drop `cname'
                        }
                        else {
                            di as error ///
                                "variable `cname' already exists; specify replace"
                            exit 110
                        }
                    }

                    quietly gen byte `cname' = ///
                        `pcand' if `touse'

                    label variable `cname' ///
                        "Custom threshold policy `crow'"
                }
            }
        }
        else {

            local cg1 : word 1 of `custom'
            local cg2 : word 2 of `custom'

            quietly replace `pcand' = ///
                (`z1' >= `cg1' & `z2' >= `cg2') if `touse'

            quietly opl_tb_cba__eval,               ///
                policy(`pcand')                          ///
                tau(`tau')                              ///
                cost(`cost')                            ///
                lambda(`lambda')                        ///
                touse(`touse')

            matrix `CUSTOMMAT' = J(1, 10, .)

            matrix colnames `CUSTOMMAT' =                ///
                threshold1                               ///
                threshold2                               ///
                N_treated                                ///
                coverage                                 ///
                total_benefit                            ///
                total_cost                               ///
                Q                                        ///
                ATET                                     ///
                avg_cost                                 ///
                avg_surplus

            matrix `CUSTOMMAT'[1,1]  = `cg1'
            matrix `CUSTOMMAT'[1,2]  = `cg2'
            matrix `CUSTOMMAT'[1,3]  = r(N_treated)
            matrix `CUSTOMMAT'[1,4]  = r(coverage)
            matrix `CUSTOMMAT'[1,5]  = r(total_benefit)
            matrix `CUSTOMMAT'[1,6]  = r(total_cost)
            matrix `CUSTOMMAT'[1,7]  = r(Q)
            matrix `CUSTOMMAT'[1,8]  = r(ATET)
            matrix `CUSTOMMAT'[1,9]  = r(avg_cost)
            matrix `CUSTOMMAT'[1,10] = r(avg_surplus)

            if "`generate'" != "" {

                local cname = "`generate'_custom"

                capture confirm new variable `cname'

                if _rc {

                    if "`replace'" != "" {
                        capture drop `cname'
                    }
                    else {
                        di as error ///
                            "variable `cname' already exists; specify replace"
                        exit 110
                    }
                }

                quietly gen byte `cname' = ///
                    `pcand' if `touse'

                label variable `cname' ///
                    "Custom two-threshold policy"
            }
        }
    }

    /*
    Generate requested variables.
    */
    if "`generate'" != "" {

        foreach suffix in fb opt surplus z1 {

            local vname = "`generate'_`suffix'"

            capture confirm new variable `vname'

            if _rc {

                if "`replace'" != "" {
                    capture drop `vname'
                }
                else {
                    di as error ///
                        "variable `vname' already exists; specify replace"
                    exit 110
                }
            }
        }

        if `nselect' == 2 {

            local vname = "`generate'_z2"

            capture confirm new variable `vname'

            if _rc {

                if "`replace'" != "" {
                    capture drop `vname'
                }
                else {
                    di as error ///
                        "variable `vname' already exists; specify replace"
                    exit 110
                }
            }
        }

        quietly gen byte `generate'_fb = ///
            `pfb' if `touse'

        quietly gen byte `generate'_opt = ///
            `popt' if `touse'

        quietly gen double `generate'_surplus = ///
            `net' if `touse'

        quietly gen double `generate'_z1 = ///
            `z1' if `touse'

        if `nselect' == 2 {
            quietly gen double `generate'_z2 = ///
                `z2' if `touse'
        }

        label variable `generate'_fb ///
            "First-best treatment policy"

        label variable `generate'_opt ///
            "Optimal constrained threshold policy"

        label variable `generate'_surplus ///
            "Individual weighted surplus: tau - lambda*cost"

        label variable `generate'_z1 ///
            "`x1' standardized to [0,1]"

        if `nselect' == 2 {
            label variable `generate'_z2 ///
                "`x2' standardized to [0,1]"
        }
    }

    /*
    Graph.
    */
    if "`graph'" != "" {

        if `nselect' == 1 {

            local xlopt ""

            if `bestg1' < . {
                local xlopt ///
                    "xline(`bestg1', lpattern(solid))"
            }

            quietly twoway                                          ///
                (scatter `net' `z1' if `touse' & `popt' == 0,       ///
                    msymbol(Oh))                                    ///
                (scatter `net' `z1' if `touse' & `popt' == 1,       ///
                    msymbol(O)),                                    ///
                `xlopt'                                             ///
                yline(0, lpattern(dash))                            ///
                xtitle("`x1' standardized to [0,1]")                ///
                ytitle("Individual weighted surplus")               ///
                title("Optimal threshold welfare policy")           ///
                legend(order(1 "Not selected" 2 "Selected"))
        }
        else {

            local lines ""

            if `bestg1' < . {
                local lines ///
                    "`lines' xline(`bestg1', lpattern(solid))"
            }

            if `bestg2' < . {
                local lines ///
                    "`lines' yline(`bestg2', lpattern(solid))"
            }

            quietly twoway                                          ///
                (scatter `z2' `z1' if `touse' & `popt' == 0,        ///
                    msymbol(Oh))                                    ///
                (scatter `z2' `z1' if `touse' & `popt' == 1,        ///
                    msymbol(O)),                                    ///
                `lines'                                             ///
                xtitle("`x1' standardized to [0,1]")                ///
                ytitle("`x2' standardized to [0,1]")                ///
                title("Optimal two-threshold welfare policy")       ///
                legend(order(1 "Not selected" 2 "Selected"))
        }
    }

    /*
    Post e()-returns.
    */
    tempname b V

    matrix `b' = J(1,1,0)
    matrix colnames `b' = objective

    matrix `V' = J(1,1,0)
    matrix rownames `V' = objective
    matrix colnames `V' = objective

    ereturn post `b' `V', esample(`touse')

    ereturn scalar N         = `N'
    ereturn scalar lambda    = `lambda'
    ereturn scalar nselect   = `nselect'
    ereturn scalar npoints   = `npoints'

    ereturn scalar npolicies = ///
        cond(`nselect' == 1, `npoints', `npoints' * `npoints')

    ereturn scalar x1_min = `x1min'
    ereturn scalar x1_max = `x1max'

    if `nselect' == 2 {
        ereturn scalar x2_min = `x2min'
        ereturn scalar x2_max = `x2max'
    }

    ereturn scalar threshold1     = `bestg1'
    ereturn scalar threshold1_lev = `bestg1_lev'

    if `nselect' == 2 {
        ereturn scalar threshold2     = `bestg2'
        ereturn scalar threshold2_lev = `bestg2_lev'
    }

    /*
    First-best returns.
    */
    ereturn scalar N_fb             = `N_fb'
    ereturn scalar coverage_fb      = `coverage_fb'
    ereturn scalar total_benefit_fb = `TB_fb'
    ereturn scalar total_cost_fb    = `TC_fb'
    ereturn scalar Q_fb             = `Q_fb'
    ereturn scalar ATET_fb          = `ATET_fb'
    ereturn scalar avg_cost_fb      = `avgcost_fb'
    ereturn scalar avg_surplus_fb   = `avgsurp_fb'

    /*
    Optimal threshold policy returns.
    */
    ereturn scalar N_opt             = `N_opt'
    ereturn scalar coverage_opt      = `coverage_opt'
    ereturn scalar total_benefit_opt = `TB_opt'
    ereturn scalar total_cost_opt    = `TC_opt'
    ereturn scalar Q_opt             = `Q_opt'
    ereturn scalar ATET_opt          = `ATET_opt'
    ereturn scalar avg_cost_opt      = `avgcost_opt'
    ereturn scalar avg_surplus_opt   = `avgsurp_opt'

    /*
    Welfare loss relative to first-best.
    */
    ereturn scalar Q_loss = `Q_fb' - `Q_opt'

    if `Q_fb' > 0 {
        ereturn scalar Q_ratio = `Q_opt' / `Q_fb'
    }
    else {
        ereturn scalar Q_ratio = .
    }

    /*
    User-supplied customized policy returns.
    */
    ereturn scalar has_custompolicy = `has_custompolicy'

    if "`custompolicy'" != "" {

        ereturn scalar N_custompolicy = ///
            `N_custompolicy'

        ereturn scalar coverage_custompolicy = ///
            `coverage_custompolicy'

        ereturn scalar total_benefit_custompolicy = ///
            `TB_custompolicy'

        ereturn scalar total_cost_custompolicy = ///
            `TC_custompolicy'

        ereturn scalar Q_custompolicy = ///
            `Q_custompolicy'

        ereturn scalar ATET_custompolicy = ///
            `ATET_custompolicy'

        ereturn scalar avg_cost_custompolicy = ///
            `avgcost_custompolicy'

        ereturn scalar avg_surplus_custompolicy = ///
            `avgsurp_custompolicy'

        ereturn scalar Q_loss_custompolicy = ///
            `Q_fb' - `Q_custompolicy'

        if `Q_fb' > 0 {
            ereturn scalar Q_ratio_custompolicy = ///
                `Q_custompolicy' / `Q_fb'
        }
        else {
            ereturn scalar Q_ratio_custompolicy = .
        }

        ereturn matrix custompolicy = `CUSTOMPOLICYMAT'
        ereturn local custompolicyvar "`custompolicy'"
    }

    ereturn matrix grid = `GRID'

    if "`custom'" != "" {
        ereturn matrix custom = `CUSTOMMAT'
    }

    ereturn local cmd "opl_tb_cba"
    ereturn local cmdline `"`0'"'
    ereturn local tauvar "`tau'"
    ereturn local costvar "`cost'"
    ereturn local selectvars "`select'"
    ereturn local policytype "threshold"
    ereturn local direction "upper"
    ereturn local generate "`generate'"

    /*
    Display main comparison.
    */
    di as text _newline ///
        "Welfare-maximizing threshold policy"

    di as text "{hline 78}"

    di as text ///
        "Objective: sum_i pi_i * (tau_i - lambda * cost_i)"

    di as text ///
        "lambda = " as result %9.4f `lambda'

    di as text "{hline 78}"

    di as text ///
        %30s "Statistic"            ///
        %20s "First best"           ///
        %20s "Constrained"

    di as text "{hline 78}"

    di as text ///
        %30s "Treated units"        ///
        as result %20.0f `N_fb'     ///
        %20.0f `N_opt'

    di as text ///
        %30s "Coverage"             ///
        as result %20.4f `coverage_fb' ///
        %20.4f `coverage_opt'

    di as text ///
        %30s "Total benefit"        ///
        as result %20.4f `TB_fb'    ///
        %20.4f `TB_opt'

    di as text ///
        %30s "Total cost"           ///
        as result %20.4f `TC_fb'    ///
        %20.4f `TC_opt'

    di as text ///
        %30s "Weighted surplus Q"   ///
        as result %20.4f `Q_fb'     ///
        %20.4f `Q_opt'

    di as text ///
        %30s "ATET"                 ///
        as result %20.4f `ATET_fb'  ///
        %20.4f `ATET_opt'

    di as text ///
        %30s "Average cost"         ///
        as result %20.4f `avgcost_fb' ///
        %20.4f `avgcost_opt'

    di as text ///
        %30s "Average surplus"      ///
        as result %20.4f `avgsurp_fb' ///
        %20.4f `avgsurp_opt'

    di as text "{hline 78}"

    /*
    Display user-supplied policy.
    */
    if "`custompolicy'" != "" {

        di as text _newline ///
            "User-supplied customized policy: " ///
            as result "`custompolicy'"

        di as text "{hline 58}"

        di as text ///
            %30s "Statistic"        ///
            %20s "Customized policy"

        di as text "{hline 58}"

        di as text ///
            %30s "Treated units"    ///
            as result %20.0f `N_custompolicy'

        di as text ///
            %30s "Coverage"         ///
            as result %20.4f `coverage_custompolicy'

        di as text ///
            %30s "Total benefit"    ///
            as result %20.4f `TB_custompolicy'

        di as text ///
            %30s "Total cost"       ///
            as result %20.4f `TC_custompolicy'

        di as text ///
            %30s "Weighted surplus Q" ///
            as result %20.4f `Q_custompolicy'

        di as text ///
            %30s "ATET"             ///
            as result %20.4f `ATET_custompolicy'

        di as text ///
            %30s "Average cost"     ///
            as result %20.4f `avgcost_custompolicy'

        di as text ///
            %30s "Average surplus"  ///
            as result %20.4f `avgsurp_custompolicy'

        di as text "{hline 58}"
    }

    /*
    Display optimal thresholds.
    */
    if `bestg1' < . {

        di as text ///
            "Optimal standardized threshold 1: " ///
            as result %7.4f `bestg1'

        di as text ///
            "Optimal original-scale threshold 1: " ///
            as result %10.4f `bestg1_lev'
    }
    else {

        di as text ///
            "Optimal constrained policy: " ///
            as result "treat none"
    }

    if `nselect' == 2 & `bestg2' < . {

        di as text ///
            "Optimal standardized threshold 2: " ///
            as result %7.4f `bestg2'

        di as text ///
            "Optimal original-scale threshold 2: " ///
            as result %10.4f `bestg2_lev'
    }
end


*! version 1.1.0 14jul2026
capture program drop opl_tb_cba__eval
program define opl_tb_cba__eval, rclass
    version 18.0

    syntax,                                                 ///
        POLICY(varname numeric)                             ///
        TAU(varname numeric)                                ///
        COST(varname numeric)                               ///
        LAMBDA(real)                                        ///
        TOUSE(varname numeric)

    quietly count if `touse'
    local N = r(N)

    quietly count if `touse' & `policy' == 1
    local Nt = r(N)

    if `Nt' == 0 {

        return scalar N_treated    = 0
        return scalar coverage     = 0
        return scalar total_benefit = 0
        return scalar total_cost    = 0
        return scalar Q             = 0
        return scalar ATET          = .
        return scalar avg_cost      = .
        return scalar avg_surplus   = .

        exit
    }

    quietly summarize `tau' if ///
        `touse' & `policy' == 1, meanonly

    local TB   = r(sum)
    local ATET = r(mean)

    quietly summarize `cost' if ///
        `touse' & `policy' == 1, meanonly

    local TC = r(sum)
    local AC = r(mean)

    local Q = `TB' - `lambda' * `TC'

    return scalar N_treated     = `Nt'
    return scalar coverage      = `Nt' / `N'
    return scalar total_benefit = `TB'
    return scalar total_cost    = `TC'
    return scalar Q             = `Q'
    return scalar ATET          = `ATET'
    return scalar avg_cost      = `AC'
    return scalar avg_surplus   = `Q' / `Nt'
end