/*
The Inequality Key Figures disseminated on the LIS website are computed using the R programming language. When attempting to replicate these figures in Stata, small discrepancies may occasionally arise. These differences are due to variations in numerical precision between the two software environments, as well as dataset-specific characteristics. For example, in datasets with very small subsamples (such as poorsm – rs19), results may be highly sensitive to whether a few households are classified as poor or not. In other datasets with many identical values, even minor precision differences can affect poverty classification and lead to noticeable mismatches.

Although such discrepancies are very rare, users seeking full replication of the published figures are advised to use R.

For further questions, please contact the LIS User Support team at usersupport@lisdatacenter.org
*/




// package = stata
// project = lis

// To select specific datasets, other than an entire country series, use option 'ccyy()' instead of 'iso2()' --> example: ccyy(de22 it16)


	// 1) Load data 

lissyuse, iso2(lu) ///
		  hvars(hid hwgt nhhmem dhi dname year) // (data is at household-level)

levelsof dname, local(levels)
foreach ccyy of local levels {
	
**
preserve
**
di "`ccyy'"
keep if dname == "`ccyy'"


	// 2) Data preparation 

*===================================
* Filter out missing observations 
*===================================
 
drop if missing(dhi)
 
*===================================
* Create person weight (data is at household-level)
*===================================
 
generate double pwt = hwgt * nhhmem 
 
*===================================
* Bottom and top coding / outlier detection
*===================================

generate double dhi_log = log(dhi) 
replace dhi_log = 0 if dhi_log == . & dhi != . // keep negatives and 0 

sort dhi_log hwgt
qui percentils dhi_log [aw = hwgt], p(25 75)
local p75 = e(Perc_75)
local p25 = e(Perc_25) 
gen double iqr = `p75' - `p25'  // interquartile range
gen double upper_bound = `p75'  + (iqr * 3)  // upper bound for extreme values
gen double lower_bound = `p25'  - (iqr * 3)  // lower bound for extreme values
replace dhi=exp(upper_bound) if dhi>exp(upper_bound)  // top code income at upper bound for extreme values
replace dhi=exp(lower_bound) if dhi<exp(lower_bound)  // bottom code income at lower bound for extreme values

*===================================
* Equivalise income (square root)
*===================================
generate double ey = ( dhi / (nhhmem^0.5) )


	// 3) Compute Key Figures 

*===================================
* Gini Coefficient
*===================================

quietly sum ey [w = pwt] , de
local cumsum_w = r(sum_w)
local mean_ey = r(mean) 
sort ey
egen double gini = sum( (2 * sum(pwt) + 1 - pwt) * pwt * (ey - `mean_ey') / (`cumsum_w'^2 * `mean_ey') )
dis "Gini Coefficient = " gini

*===================================
* Atkinson Coefficients
*===================================

* epsilon equal to 0.5	
local epsilon "0.5"
quietly sum ey [aw=pwt] 
quietly gen double new0p5 = sum(pwt * ((ey/r(mean))^(1-`epsilon')))    
quietly sum ey [aw=pwt] , de 
dis "Atkinson Coefficient (e=0.5) = " 1-(new0p5[_N]/r(sum_w))^(1/(1-`epsilon'))

* epsilon equal to 1				
quietly sum ey [aw=pwt] 
quietly gen double new1 = exp(sum((pwt*log((ey/r(mean))^(1/r(sum_w))))))
dis "Atkinson Coefficient (e=1) = " 1-(new1[_N])
 
*===================================
* Percentile Ratios
*===================================
 
sort ey pwt 
qui percentils ey [aw = pwt], p(90 10)
local p90 = e(Perc_90)
local p10 = e(Perc_10) 
generate double ratio_p90_p10 = `p90'/`p10'
dis "Percentile Ratio (90/10) = " ratio_p90_p10 

sort ey pwt 
qui percentils ey [aw = pwt], p(90 50)
local p90 = e(Perc_90)
local p50 = e(Perc_50) 
generate double ratio_p90_p50 = `p90'/`p50'
dis "Percentile Ratio (90/50) = " ratio_p90_p50 

sort ey pwt 
qui percentils ey [aw = pwt], p(80 20)
local p80 = e(Perc_80)
local p20 = e(Perc_20) 
generate double ratio_p80_p20 = `p80'/`p20'
dis "Percentile Ratio (80/20) = " ratio_p80_p20 

**
restore
**

}