***************************************************************************************************;
#delimit;
** define the program that calculates the Pendakur, Pendakur and Woodcock (2008) Severity Index;
**    cite Pendakur, Krishna, Ravi Pendakur and Simon D. Woodcock, 2008, A Representation Index:        ;
**         Measuring the Representation of Minorities in the Income Distribution, Berkeley Electronic   ;
**         Journal of Economic Analysis and Policy: ADVANCES.  						; 
** this program uses UNIQUE;
** -- VARLIST is dependent variable followed by covariates;
** -- can use weighted data with EWEIGHT option (with UNWEIGHTED set to 0);
** -- QTILE selects the quantile, and ABOVE specifies above measures (with below measures as default);
** -- GROUPS is a single variable taking a different value for each group;
** -- ANCHOR is a scalar specifying the group number of the anchoring distribution;
** --   (0 indicates the population as the anchoring distribution);
** -- TOL gives the convergence tolerance for the iterated WLS expectile estimator;
program define severity, sortpreserve byable(recall);
version 10;
syntax varlist [if] [in], [etile(real 0.5)] Groups(varlist numeric) [above] [anchor (integer 0)] [tol (real 0.00001)] [EWeight (varlist numeric)] [UNweighted (integer 1)];
    marksample touse;
    preserve;
    tokenize `varlist';
    local depvar `1';
        mac shift;
        local xvar `*';
    unique `groups' if `touse';
    scalar ng=r(sum);
    di "severity index will be computed on " ng " groups using group " `anchor' " as the anchor";
    if "`above'"=="above" di "calculating severity ABOVE " `etile';
    else di "calculating severity BELOW " `etile' " expectile";   
    matrix rep = J(ng,1,0);     * this holds severity index;
    * get point estimates;
    di "computing expectiles ... ";
    g anchor_group=`groups';
    if `anchor'==0 {;
	quietly replace anchor_group=0;
    };
    if `unweighted'==1 {;
	quietly g aw=1;
    };
    else {;
	quietly g aw=`eweight';
    };
    * pre-estimate;
    quietly reg `varlist' [aweight=aw] if ((`touse')&(anchor_group==`anchor'));
    predict ehat, residual;
    g ehat_old=ehat;
    g edif=abs(ehat-ehat_old);
    g ew=abs(`etile'-(ehat<0))*aw;
    scalar e_test=10;
    * iterate to convergence;
    while e_test>`tol' {;
	    quietly replace ehat_old=ehat;
	    quietly drop ehat;
	    quietly reg `varlist' [aweight=ew] if ((`touse')&(anchor_group==`anchor'));
	    quietly predict ehat, residual;
	    quietly replace ew=abs(`etile'-(ehat<0))*aw;
	    quietly replace edif=abs(ehat-ehat_old);
	    quietly summ edif;
	    scalar e_test=r(max);
    };
    quietly predict yhat ; 
    if "`above'"=="above" {;
        quietly g abs_yhat_minus_y=abs(yhat-`depvar') if ((`touse')&(`depvar'~=.)&(yhat~=.));
	quietly g trunc_yhat_minus_y=(yhat<`depvar')*abs(yhat-`depvar') if ((`touse')&(`depvar'~=.)&(yhat~=.));
   };
    else {;
        quietly g abs_yhat_minus_y=abs(yhat-`depvar') if ((`touse')&(`depvar'~=.)&(yhat~=.));
	quietly g trunc_yhat_minus_y=(yhat>`depvar')*abs(yhat-`depvar') if ((`touse')&(`depvar'~=.)&(yhat~=.));
    };
    local gg=1;
    di "severity of the following group numbers";
    levelsof `groups', local(levels);
    di "using anchor group " `anchor';
    foreach g of local levels {;
        quietly summ trunc_yhat_minus_y [aweight=aw] if ((`touse') & (`groups'==`g')&(`depvar'~=.)&(yhat~=.)) ; 
	scalar num=r(sum); 
        quietly summ abs_yhat_minus_y [aweight=aw] if ((`touse') & (`groups'==`g')&(`depvar'~=.)&(yhat~=.)) ; 
	scalar den=r(sum); 
        matrix rep[`gg',1] = num/den;
        local gg = `gg' + 1;
    };
    quietly drop yhat anchor_group aw ew ehat ehat_old abs_yhat_minus_y trunc_yhat_minus_y;
    matrix rep=rep';
    matrix list rep;
restore;
end;