/*************************************************************************/
/* DASP  : Distributive Analysis Stata Package  (version 2.3)            */
/*************************************************************************/
/* Conceived and programmed by Dr. Araar Abdelkrim   (2006-2012)          */
/* Universite Laval, Quebec, Canada                                      */
/* email : aabd@ecn.ulaval.ca                                            */
/* Phone : 1 418 656 7507                                                */
/*************************************************************************/
/* module  : dinineq                                                       */
/*************************************************************************/

** EDITED BY SEAN HIGGINS:
**  07-03-2016 Replace inefficient while loop through observations
**			   Formatting and remove table of results
**			   Calculate and return p-value of test of difference = 0

#delim ;

/*****************************************************/
/* Density function      : fw=Hweight*Hsize          */
/*****************************************************/
cap program drop ceqdinineq_den;
program define ceqdinineq_den, rclass;
	args fw x xval;
	qui su `x' [aw=`fw'], detail;
	local tmp = (`r(p75)'-`r(p25)')/1.34;
	local tmp = (`tmp'<`r(sd)')*`tmp'+(`tmp'>=`r(sd)')*`r(sd)' ;
	local h   = 0.9*`tmp'*_N^(-1.0/5.0);
	tempvar s1 s2;
	gen `s1' = sum( `fw' *exp(-0.5* ( ((`xval'-`x')/`h')^2  )  ));
	gen `s2' = sum( `fw' );
	return scalar den = `s1'[_N]/( `h'* sqrt(2*c(pi)) * `s2'[_N] );
end;

/***************************************/
/* Quantile  & GLorenz                 */
/***************************************/
cap program drop ceqdinineq_qua;
program define ceqdinineq_qua, rclass sortpreserve;
	args fw yyy xval order;
	preserve;
	sort `yyy', stable;
	qui cap drop if `yyy'>=. | `fw'>=.;
	tempvar ww qp glp pc;
	qui gen `ww'=sum(`fw');
	qui gen `pc'=`ww'/`ww'[_N];
	qui gen `qp' = `yyy' ;
	qui gen `glp' double = sum(`fw'*`yyy')/`ww'[_N];
	qui sum `yyy' [aw=`fw'];

	/* added by Sean to replace inefficient while loop through observations */
	tempvar ordering ;
	gen `ordering' = _n;
	count if (`pc'<`xval');
	if r(N) { ; // i.e. positive number of obs with `pc'<`xval' ;
		summ `ordering' if (`pc'<`xval');
		local ar = r(max);
	} ;
	else { ;
		local ar = 0;
	};
	local i = `ar' + 1;

	if (`i'> 1) {;
		local qnt =`qp'[`ar'] +((`qp'[`i'] -`qp'[`ar']) /(`pc'[`i']-`pc'[`ar']))*(`pc'[`i']-`pc'[`ar']);
		local glor=`glp'[`ar']+((`glp'[`i']-`glp'[`ar'])/(`pc'[`i']-`pc'[`ar']))*(`pc'[`i']-`pc'[`ar']);
	};
	if (`i'==1) {;
		local qnt =(max(0,`qp'[`i'])/(`pc'[`i']))*(`pc'[`i']);
		local glor=(max(0,`glp'[`i'])/(`pc'[`i']))*(`pc'[`i']);
	};

	return scalar qnt  = `qnt';
	return scalar glor = `glor';
	restore;
end;

/*********************************/
/* Main programs                 */
/*********************************/
cap program drop ceqdinineq2;
program define ceqdinineq2, rclass ;
	version 9.2;
	syntax varlist (min=1 max=1)
		[,
			HSize(varname)
			HWeight(varname)
			HGroup(varname)
			p1(real 0.1) p2(real 0.2) p3(real 0.8) p4(real 0.9)
			index(string)
			GNumber(int -1)
			CI(real 5)  CONF(string) LEVEL(real 95) vab(int 0)
		];

	tokenize `varlist';
						qui drop if `1'>=. ;
	if ("`hsize'"!="")   qui drop if `hsize'>=.;
	if ("`hweight'"!="") qui drop if `hweight'>=.;
	tempvar  hs sw fw ;
	gen `sw'=1;
	gen `hs'=1;

	if ("`hsize'"!="")     qui replace `hs' = `hsize';
	tempvar _in;
	if ("`hgroup'" != "")  qui gen    `_in' = (`hgroup' == `gnumber');
	if ("`hgroup'" != "")  qui replace `hs' = `hs' * `_in';
	if ("`hweight'"!="")   qui replace `sw'=`hweight';


	local hweight="";
	cap qui svy: total `1';
	local hweight=`"`e(wvar)'"'; // " ;
	cap ereturn clear;
	tempvar fw;
	gen `fw'=`hs';
	if (`"`hweight'"'~="") qui replace `fw'=`fw'*`hweight'; // " ;

	tempvar vec_a vec_b;

	if ( "`index'"=="qr") {;
		ceqdinineq_qua `fw' `1' `p1';
		local q1=`r(qnt)';
		ceqdinineq_qua `fw' `1' `p2';
		local q2=`r(qnt)';
		local est = `q1'/`q2';
		ceqdinineq_den `fw' `1' `q1';
		local fq1=`r(den)';
		ceqdinineq_den `fw' `1' `q2';
		local fq2=`r(den)';
		gen `vec_a' = -`hs'*((`q1'>`1')-`p1')/`fq1' + `hs'*`q1';
		gen `vec_b' = -`hs'*((`q2'>`1')-`p2')/`fq2' + `hs'*`q2';
		qui svy: ratio `vec_a'/`vec_b';
		cap drop matrix _vv;
		matrix _vv=e(V);
		local std = el(_vv,1,1)^0.5;
		global ws1=`q1';
		global ws2=`q2';
	};

	if ( "`index'"=="sr") {;
		ceqdinineq_qua `fw' `1' `p1';
		local q1=`r(qnt)'; local g1=`r(glor)';
		ceqdinineq_qua `fw' `1' `p2';
		local q2=`r(qnt)'; local g2=`r(glor)';
		ceqdinineq_qua `fw' `1' `p3';
		local q3=`r(qnt)'; local g3=`r(glor)';
		ceqdinineq_qua `fw' `1' `p4';
		local q4=`r(qnt)'; local g4=`r(glor)';

		local est = (`g2'-`g1')/(`g4'-`g3');
		global ws1=(`g2'-`g1');
		global ws2=(`g4'-`g3');

		gen `vec_a' = `hs'*(`q2'*`p2'+(`1'-`q2')*(`q2'>`1')) - `hs'*(`q1'*`p1'+(`1'-`q1')*(`q1'>`1')) ;
		gen `vec_b' = `hs'*(`q4'*`p4'+(`1'-`q4')*(`q4'>`1')) - `hs'*(`q3'*`p3'+(`1'-`q3')*(`q3'>`1')) ;;
		qui svy: ratio `vec_a'/`vec_b';
		cap drop matrix _vv;
		matrix _vv=e(V);
		local std = el(_vv,1,1)^0.5;

	};

	if (`vab'==1) {; // note that if the two vars come from the same data set, `vab' was set to 1 by ceqdinineq;
		gen __va=`vec_a';
		gen __vb=`vec_b';
	};

	qui svydes;
	local fr=`r(N_units)'-`r(N_strata)';
	local lvl=(100-`level')/100;
	if ("`conf'"=="ts") local lvl=`lvl'/2;
	local tt=invttail(`fr',`lvl');

	return scalar est  = `est';
	return scalar std  = `std';
	return scalar lb   = `est' - `tt'*`std';
	return scalar ub   = `est' + `tt'*`std';
	return scalar df   = `fr';

end;

capture program drop ceqdinineq;
program define ceqdinineq, eclass;
	version 9.2;
	syntax  namelist(min=2 max=2)
		[,
			FILE1(string) FILE2(string)
			p1(real 0.1) p2(real 0.2) p3(real 0.8) p4(real 0.9)
			index(string)
			HSize1(string) HSize2(string)
			COND1(string)  COND2(string)
			type(string)
			LEVEL(real 95) CONF(string) TEST(string) DEC(int 6)
		];

	global indica=3;
	tokenize `namelist';
	if ("`conf'"=="")          local conf="ts";
	if ("`index'"=="") local index = "qr";
	preserve;
	local indep = ( (`"`file1'"'=="" & `"`file2'"'=="") | (`"`file1'"'==`"`file2'"')  ); // " ;
	local vab=0;
	if ( (`"`file1'"'=="" & `"`file2'"'=="") | (`"`file1'"'==`"`file2'"')  ) local vab=1;  // " ;
	if ("`file1'" !="") use `"`file1'"', replace;  // " ;
	tempvar cd1;
	tempvar ths1;
	qui gen `ths1'=1;

	if ( "`hsize1'"!="") qui replace `ths1'=`hsize1';

	if ( "`cond1'"!="") {;
		gen `cd1'=`cond1';
		qui replace `ths1'=`ths1'*`cd1';
		qui sum `cd1';
		if (`r(sum)'==0) {;
			dis as error " With the condition(s) of distribution_1, the number of observations is 0.";
			exit;
		};
	};

	qui svyset ;
	if ( "`r(settings)'"==", clear") qui svyset _n, vce(linearized);

	local hweight1="";
	cap qui svy: total `1';
	local hweight1=`"`e(wvar)'"';  // " ;
	cap ereturn clear;

	ceqdinineq2 `1' ,
		hweight(`hweight1') hsize(`ths1')
		p1(`p1') p2(`p2') p3(`p3') p4(`p4')
		index(`index')
		conf(`conf') level(`level')
		vab(`vab')
	;
	matrix _res_d1  =(`r(est)',`r(std)',`r(lb)',`r(ub)', `r(df)') ;

	if (`vab'==1) {; // note that if the two vars come from the same data set, `vab' was set to 1 above;
		tempvar va vb;
		qui gen `va'=__va;
		qui gen `vb'=__vb;
		qui drop __va __vb;
		local sv1=$ws1;
		local sv2=$ws2;
	};

	if ("`file2'" !="" & `vab'!=1) use `"`file2'"', replace; // " ;
	tempvar cd2 ths2;
	qui gen `ths2'=1;
	if ( "`hsize2'"!="") qui replace `ths2'=`hsize2';
	if ( "`cond2'"!="") {;
		gen `cd2'=`cond2';
		qui replace `ths2'=`ths2'*`cd2';
		qui sum `cd2';
		if (`r(sum)'==0) {;
			dis as error " With the condition(s) of distribution_2 the number of observations is 0.";
			exit;
		};
	};

	qui svyset ;
	if ( "`r(settings)'"==", clear") qui svyset _n, vce(linearized);
	local hweight2="";
	cap qui svy: total `2';
	local hweight2=`"`e(wvar)'"'; // " ;
	cap ereturn clear;

	ceqdinineq2 `2' ,
		hweight(`hweight2') hsize(`ths2')
		p1(`p1') p2(`p2') p3(`p3') p4(`p4')
		index(`index')
		conf(`conf') level(`level')
		vab(`vab')
	;
	if (`vab'==1) {;
		tempvar vc vd;
		qui gen `vc'=__va;
		qui gen `vd'=__vb;
		qui drop __va __vb;
		local sv3=$ws1;
		local sv4=$ws2;
	};

	matrix _res_d2 =(`r(est)',`r(std)',`r(lb)',`r(ub)', `r(df)');
	local dif = el(_res_d2,1,1)-el(_res_d1,1,1);
	local std = (el(_res_d1,1,2)^2+el(_res_d2,1,2)^2)^0.5;

	if (`vab'==1) {;
		qui svy: mean `va' `vb' `vc' `vd';
		cap drop matrix mat;
		matrix mat=e(V);
		cap drop matrix gra;
		matrix gra=
			(
				1/`sv2'			\
				-`sv1'/`sv2'^2	\
				-1/`sv4'		\
				`sv3'/`sv4'^2
			);
		cap matrix drop _zz;
		matrix _zz=gra'*mat*gra;
		local std= el(_zz,1,1)^0.5;
	};

	local sdif = `std';

	local ind = "Quantile ratio";

	local lr = "p1 = `p1'";
	local hr = "p2 = `p2'";

	if ("`index'"=="sr") {;
		local ind = "Share ratio";
		local lr = "p1 = `p1' / p2 = `p2'";
		local hr = "p3 = `p3' / p4 = `p4'";
	};

	/*
      tempname table;
	.`table'  = ._tab.new, col(5);
	.`table'.width  24|16 16 16 16 ;
	.`table'.strcolor . . yellow . . ;
	.`table'.numcolor yellow yellow . yellow yellow;
	.`table'.numfmt %16.0g  %16.`dec'f  %16.`dec'f %16.`dec'f %16.`dec'f ;
       di _n as text in white "{col 5}Difference: `ind' index of inequality";
       di as text     "{col 5}Lower  rank      :  `lr'";
       di as text     "{col 5}Higher rank      :  `hr'";
	*/

	local sdif = `std';

	if ("`conf'"!="ts") local lvl=(100-`level')/100;
	if ("`conf'"=="ts") local lvl = (1-(100-`level')/200);
	local zzz=invnorm(`lvl');
	local lb   = `dif' -  `zzz'*`std';
	local ub   = `dif' +  `zzz'*`std';
	matrix _res_di =(`dif',`std',`lb',`ub');

	/* added by Sean: calculate p-value */
	qui svydes;
	local fr=`r(N_units)'-`r(N_strata)';
	scalar _p = ttail(`fr',abs(`dif'/`sdif'));

	local est1=el(_res_d1,1,1);
	local est2=el(_res_d2,1,1);

	local std1=el(_res_d1,1,2);
	local std2=el(_res_d2,1,2);

	local df1=el(_res_d1,1,5);
	local df2=el(_res_d2,1,5);

	local tyind1="Dist1";
	local tyind2="Dist2";

	/*
	_dasp_dif_table `2' `2',
	name1("`tyind1'")          name2("`tyind2'")
	m1(`est1')            m2(`est2')
	s1(`std1')            s2(`std2')
	df1(`df1')            df2(`df2')
	dif(`dif') sdif(`sdif')
	level(`level') conf(`conf') indep(`indep') test(`test');
	*/

	restore;
	ereturn clear;
	ereturn matrix d1 = _res_d1;
	ereturn matrix d2 = _res_d2;
	ereturn matrix di = _res_di;
	ereturn scalar p  = _p;

	cap matrix drop _res_d1;
	cap matrix drop _res_d2;
	cap matrix drop _res_di;
	cap scalar drop _p;

end;