*! version 1.0.0
*! Dean Jolliffe and Bohdan Krushelnytskyy, CERGE UK
*! July 1999
*! Bootstrap standard errors for indices of inequality.

#delimit ;

program define ineqerr;
  version 5.0;
  local varlist "required existing";
  local weight "aweight fweight";
  local if "optional";
  local in "optional";
  local options "Reps(integer 100) psu(string) psuwt(string)";
  parse "`*'";

  quietly
  {;
    if `reps'<1 {; display in red "reps() must be positive"; exit; };
    if "`psuwt'" ~= "" & "`psu'"==""
    {; display in red "psuwt() can only be specified if psu() is specified"; exit; };

    parse "`varlist'", parse(" ");
    tempvar use weivar psuwei;
    generate byte `use'=1 `if' `in';
    if "`weight'"=="" {; generate `weivar'=1; }; else {; generate `weivar'`exp'; };
    if "`psuwt'" ~= "" {; generate `psuwei'=`psuwt'; local psuwt "`psuwei'"; };
    preserve;

    if "`psu'"==""
    {;
      while "`1'" ~= ""
      {;
        restore, preserve;
        local vlabel : variable label `1';
        if "`vlabel'"=="" {; local vlabel "(unlabeled)"; };
        local llength = 77-length("`1'")-length("`vlabel'");
        noisily display _newline "`1' " _dup(`llength') "-" " `vlabel'";

        drop if `use'~=1 | `1'==. | `weivar'==.;
        keep `1' `weivar';
        count if `1'<0;
        if _result(1)>0
        {;
          noisily display _result(1) " values < 0. Not used in calculations.";
          drop if `1'<0;
        };
        count if `1'==0;
        if _result(1)>0
        {;
          noisily display _result(1) " values = 0. Not used in calculations.";
          drop if `1'==0;
        };
        noisily display "(obs=" _N ")";

        tempname gini theil varlogs memfile;
        tempfile datfile resfile;
        save `datfile', replace;

        ineqmeas `1' `weivar';
        scalar `gini'=$gini;
        scalar `theil'=$theil;
        scalar `varlogs'=$varlogs;

        postfile `memfile' Gini Theil Varlogs using `resfile', replace;
        local i=1;
        while `i'<=`reps'
        {;
          use `datfile', clear;
          bsample;
          ineqmeas `1' `weivar';
          post `memfile' $gini $theil $varlogs;
          local i=`i'+1;
        };
        postclose `memfile';

        use `resfile', clear;
        char Gini[bstrap] `gini';
        char Theil[bstrap] `theil';
        char Varlogs[bstrap] `varlogs';
        noisily bstat Gini Theil Varlogs;
        macro shift;
      };
    };
    else
    {;
      while "`1'" ~= ""
      {;
        restore, preserve;
        local vlabel : variable label `1';
        if "`vlabel'"=="" {; local vlabel "(unlabeled)"; };
        local llength = 77-length("`1'")-length("`vlabel'");
        noisily display _newline "`1' " _dup(`llength') "-" " `vlabel'";

        drop if `use'~=1 | `1'==. | `psu'==. | `weivar'==.;
        if "`psuwt'" ~= "" {; drop if `psuwt'==.; };
        count if `1'<0;
        if _result(1)>0
        {;
          noisily display _result(1) " values < 0. Not used in calculations.";
          drop if `1'<0;
        };
        count if `1'==0;
        if _result(1)>0
        {;
          noisily display _result(1) " values = 0. Not used in calculations.";
          drop if `1'==0;
        };
        noisily display "(obs=" _N ")";

        tempname gini theil varlogs memfile;
        tempfile datfile uniqpsu resfile;
        tempvar bootvar bootwei count;

        ineqmeas `1' `weivar';
        scalar `gini'=$gini;
        scalar `theil'=$theil;
        scalar `varlogs'=$varlogs;

        * save only relevant variables from the original dataset to use for further merges;
        keep `1' `psu' `weivar' `psuwt';
        sort `psu';
        save `datfile', replace;

        * create the datafile with unique PSUs and Stage1 Weights;
        sort `psu';
        by `psu': generate num=_n;
        keep if num==1;
        keep `psu' `psuwt';
        save `uniqpsu', replace;

        postfile `memfile' Gini Theil Varlogs using `resfile', replace;
        local i=1;
        while `i'<=`reps'
        {;
          * select PSUs;
          use `uniqpsu', clear;
          if "`psuwt'" ~= ""
          {;
            tempname memdraw;
            tempfile psu_tmp resdraw;
            tempvar cum_wt;

            generate `cum_wt'=sum(`psuwt');
            generate `bootvar'=`cum_wt'[_N]*uniform();
            save `psu_tmp', replace;

            postfile `memdraw' `psu' using `resdraw', replace;
            local j=1;
            while `j'<=_N
            {;
              keep if `cum_wt'>`bootvar'[`j'];
              post `memdraw' `psu'[1];
              local j=`j'+1;
              use `psu_tmp', clear;
            };
            postclose `memdraw';

            use `resdraw', clear;
            sort `psu';
            merge `psu' using `uniqpsu';
            keep if _merge==3;
            collapse (count) `count'=`psuwt', by(`psu');
          };
          else
          {;
            generate `bootvar'=`psu'[int(_N*uniform())+1];
            collapse (count) `count'=`psu', by(`bootvar');
            rename `bootvar' `psu';
          };

          * merge with the original dataset;
          merge `psu' using `datfile';
          keep if _merge==3;
          expand `count';
          keep `1' `psu' `weivar';

          * randomly select households from each PSU;
          sort `psu';
          by `psu': generate `count'=int(_N*uniform())+1;
          by `psu': generate `bootvar'=`1'[`count'];
          by `psu': generate `bootwei'=`weivar'[`count'];

          * post the inequality measures;
          ineqmeas `bootvar' `bootwei';
          post `memfile' $gini $theil $varlogs;
          local i=`i'+1;
        };
        postclose `memfile';

        use `resfile', clear;
        char Gini[bstrap] `gini';
        char Theil[bstrap] `theil';
        char Varlogs[bstrap] `varlogs';
        noisily bstat Gini Theil Varlogs;
        macro shift;
      };
    };
  macro drop gini theil varlogs;
  };
end;

program define ineqmeas;
  version 5.0;
  quietly
  {;
    gsort - `1';

    summarize `1' [aweight=`2'];
    tempname numind mean;
    scalar `numind'=_result(2);
    scalar `mean'=_result(3);

    tempvar rho newvar;
    generate `rho'=`2'[_n-1];
    replace `rho'=1 in 1;
    replace `rho'=sum(`rho');

    generate `newvar'=`2'*`1'*(`rho'+(`2'-1)/2);
    summarize `newvar';
    global gini=1+1/`numind'-2/(`mean'*`numind'^2)*_result(18);  /* Myles, STB-48 */

    replace `newvar'=`1'/`mean'*ln(`1'/`mean');
    summarize `newvar' [aweight=`2'];
    global theil=_result(3);                                     /* Deaton */

    replace `newvar'=ln(`1');
    summarize `newvar' [aweight=`2'];
    scalar `mean'=_result(3);
    replace `newvar'=(ln(`1')-`mean')^2;
    summarize `newvar' [aweight=`2'];
    global varlogs=_result(3);
  };
end;