program define diff_ge

    /*************************************************************/
    /* computes differences of some generalized entropy measures */
    /* and bootstrap confidence intervals from potentially unba- */
    /* lanced and dependent samples 1 and 2                      */
    /*                                                           */
    /* syntax: diff_ge <w1> <inc1> <c1> <w2> <inc2> <c2>         */
    /*         where  <w1> and <w2> are weights                  */
    /*                <inc1> and <inc2> are income variables     */
    /*                <c1> and <c2> are sample inclusion dummies */
    /*                                                           */
    /* incomes of non-participants in a wave (i.e. c1=0 or c2=0) */
    /* must be set to strictly positive values                   */
    /*                                                           */
    /* if no weights are needed, generate two variables with all */
    /* values set to 1 and use them as <w1> and <w2>             */
    /*                                                           */
    /* for details see Biewen (2002): Bootstrap inference for    */
    /* inequality, mobility and poverty measurement, Journal     */
    /* of Econometrics, Vol. 108, pp. 317 – 342                  */
    /*                                                           */
    /* copyright Martin Biewen 2001                              */
    /*************************************************************/

    version 6



    /* parse syntax */

    syntax varlist(min=6 max=6)



    /* set number of bootstrap replications and confidence level */

    tempname reps conf nobs

    scalar `reps' = 100

    scalar `conf' = 0.95

    scalar `nobs' = _N



    /* check if incomes > 0 */

    quietly sum `2'

    if r(min) <= 0 {

        di in red _newline "incomes must be strictly positive"
        di in red "incomes of non-participants must be set to strictly positive values"

        exit
        }

    quietly sum `5'

    if r(min) <= 0 {

        di in red _newline "incomes must be strictly positive"
        di in red "incomes of non-participants must be set to strictly positive values"

        exit
        }


    /* compute estimates */

    ge_diff `1' `2' `3' `4' `5' `6'


    tempname dgem1 dmld dtheil dge2

    scalar `dgem1' = r(dgem1)

    scalar `dmld' = r(dmld)

    scalar `dtheil' = r(dtheil)

    scalar `dge2' = r(dge2)



    /* preserve and initialize post */

    preserve

    tempname postnam

    tempfile boot

    quietly postfile `postnam' bgem1 bmld btheil bge2 using `boot', replace



    /* bootstrap loop */

    local i = 1

    while `i' <= `reps' {

        restore, preserve

        bsample

        ge_diff `1' `2' `3' `4' `5' `6'

        post `postnam' sqrt(`nobs') * (r(dgem1) - `dgem1') /*
                        */ sqrt(`nobs') * (r(dmld) - `dmld') /*
                        */ sqrt(`nobs') * (r(dtheil) - `dtheil') /*
                        */ sqrt(`nobs') * (r(dge2) - `dge2')

        local i = `i' + 1

        }

    postclose `postnam'



    /* generate lower and upper confidence bounds */

    use `boot', clear


    sort bgem1

    tempname lbgem1 ubgem1

    scalar `lbgem1' = `dgem1' - bgem1[`reps' - int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')

    scalar `ubgem1' = `dgem1' - bgem1[int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')


    sort bmld

    tempname lbmld ubmld

    scalar `lbmld' = `dmld' - bmld[`reps' - int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')

    scalar `ubmld' = `dmld' - bmld[int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')


    sort btheil

    tempname lbtheil ubtheil

    scalar `lbtheil' = `dtheil' - btheil[`reps' - int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')

    scalar `ubtheil' = `dtheil' - btheil[int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')


    sort bge2

    tempname lbge2 ubge2

    scalar `lbge2' = `dge2' - bge2[`reps' - int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')

    scalar `ubge2' = `dge2' - bge2[int((-(`conf'/2) + 0.5) * `reps')] / sqrt(`nobs')


    /* print out results */

    di _newline "------------------------------------------------------------------"
    di          "Generalized entropy differences and bootstrap confidence intervals"
    di          "------------------------------------------------------------------"
    di          "Confidence level = " %3.2f `conf' ",               " %4.0f `reps' " bootstrap replications
    di          "------------------------------------------------------------------"
    di          "GE^-1 = " in green %7.6f `dgem1' in yellow   ",     confidence interval " in green "[" %7.6f `lbgem1' " ; " %7.6f `ubgem1' "]"
    di          "MLD   = " in green %7.6f `dmld'  in yellow   ",     confidence interval " in green "[" %7.6f `lbmld' " ; " %7.6f `ubmld' "]"
    di          "THEIL = " in green %7.6f `dtheil' in yellow  ",     confidence interval " in green "[" %7.6f `lbtheil' " ; " %7.6f `ubtheil' "]"
    di          "GE^2  = " in green %7.6f `dge2'   in yellow  ",     confidence interval " in green "[" %7.6f `lbge2' " ; " %7.6f `ubge2' "]"
    di          "------------------------------------------------------------------"





    restore

end



program define ge_diff, rclass


    /* returns differences r(dgem1) r(dmld) r(dtheil) r(dge2)    */


quietly {


    /* period 1 */

    /* calculate products */

    tempvar w x c cw cwx cw1dx cwx2 cwlogx cwxlogx

    gen `w' = `1'

    gen `x' = `2'

    gen `c' = `3'

    gen `cw' = `c' * `w'

    gen `cwx' = `c' * `w' * `x'

    gen `cw1dx' = `c' * `w' / `x'

    gen `cwx2' = `c' * `w' * `x' * `x'

    gen `cwlogx' = `c' * `w' * log(`x')

    gen `cwxlogx' = `c' * `w' * `x' * log(`x')


    /* calculate conditional moments*/

    tempname c1 m0 m1 mm1 m2 t01 t11

    sum `c'

    scalar `c1' = r(mean)

    sum `cw'

    scalar `m0' = r(mean)/`c1'

    sum `cwx'

    scalar `m1' = r(mean)/`c1'

    sum `cw1dx'

    scalar `mm1' = r(mean)/`c1'

    sum `cwx2'

    scalar `m2' = r(mean)/`c1'

    sum `cwlogx'

    scalar `t01' = r(mean)/`c1'

    sum `cwxlogx'

    scalar `t11' = r(mean)/`c1'



    /* calculate indices for period 1 */

    tempname gem11 mld1 theil1 ge21

    scalar `gem11' = 0.5 * (`m1' * `mm1' /(`m0' * `m0') - 1)

    scalar `mld1' = - `t01' /`m0' + log(`m1'/`m0')

    scalar `theil1' = `t11' /`m1' - log(`m1'/`m0')

    scalar `ge21' = 0.5 * (`m0' * `m2' /(`m1' * `m1') -1)



    /* period 2 */

    /* calculate products */

    tempvar w x c cw cwx cw1dx cwx2 cwlogx cwxlogx

    gen `w' = `4'

    gen `x' = `5'

    gen `c' = `6'

    gen `cw' = `c' * `w'

    gen `cwx' = `c' * `w' * `x'

    gen `cw1dx' = `c' * `w' / `x'

    gen `cwx2' = `c' * `w' * `x' * `x'

    gen `cwlogx' = `c' * `w' * log(`x')

    gen `cwxlogx' = `c' * `w' * `x' * log(`x')


    /* calculate conditional moments*/

    tempname c1 m0 m1 mm1 m2 t01 t11

    sum `c'

    scalar `c1' = r(mean)

    sum `cw'

    scalar `m0' = r(mean)/`c1'

    sum `cwx'

    scalar `m1' = r(mean)/`c1'

    sum `cw1dx'

    scalar `mm1' = r(mean)/`c1'

    sum `cwx2'

    scalar `m2' = r(mean)/`c1'

    sum `cwlogx'

    scalar `t01' = r(mean)/`c1'

    sum `cwxlogx'

    scalar `t11' = r(mean)/`c1'



    /* calculate indices for period 2 */

    tempname gem12 mld2 theil2 ge22

    scalar `gem12' = 0.5 * (`m1' * `mm1' /(`m0' * `m0') - 1)

    scalar `mld2' = - `t01' /`m0' + log(`m1'/`m0')

    scalar `theil2' = `t11' /`m1' - log(`m1'/`m0')

    scalar `ge22' = 0.5 * (`m0' * `m2' /(`m1' * `m1') -1)



    /* return results */

    return scalar dgem1 = `gem12' - `gem11'

    return scalar dmld = `mld2' - `mld1'

    return scalar dtheil = `theil2' - `theil1'

    return scalar dge2 = `ge22' - `ge21'

    }

end