*! version 3.0.0 SPJ, October 2007 Add left-truncation option
*! version 2.1.0 Stephen P. Jenkins, 11 April 2007
*!	Revise Gini calc; add GE indices; etc
*! version 2.0.0 Stephen P. Jenkins, 28 May 2006
*! version 1.2.0 Stephen P. Jenkins, April 2004
*! Fit GB2 distribution by ML to unit record data (log parameter metric)


/*------------------------------------------------ playback request */

program define gb2lfit, eclass byable(onecall)
	version 8.2
	if replay() {
		if "`e(cmd)'" != "gb2lfit" {
			noi di as error "results for gb2lfit not found"
			exit 301
		}
		if _by() {
			error 190
		}
		Display `0'
		exit `rc'
	}
	if _by() {
		by `_byvars'`_byrc0': Estimate `0'
	}
	else	Estimate `0'
end

/*------------------------------------------------ estimation */

program define Estimate, eclass byable(recall)

	syntax varlist(max=1) [if] [in] [aw fw pw iw] [,  ///
		Avar(varlist numeric) Bvar(varlist numeric) Pvar(varlist numeric) Qvar(varlist numeric) ///
		ABPQ(varlist numeric) CENSvar(varname) From(string) ///
		CDF(namelist max=1) PDF(namelist max=1) POORfrac(real 0) ///
		Robust Cluster(varname) SVY  STats  ///
		Level(integer $S_level) CGINI EPSilon(real 1e-10) EXtras ///
		IGINI NIPS(integer 5000) LEFTtr(varname) ///
		noLOG  * ]

	local title "ML fit of GB2 distribution"

	local inc "`varlist'"

	if "`abpq'" ~= "" & (("`avar'"~="")|("`bvar'"~="")|("`pvar'"~="")|("`qvar'"~="")) {
		di as error "Cannot use abpq(.) option in conjunction with avar(.), bvar(.), pvar(.), qvar(.) options"
		exit 198
	}

	if "`abpq'" ~= "" {
		local avar "`abpq'"
		local bvar "`abpq'"
		local pvar "`abpq'"
		local qvar "`abpq'"
	}

	local na : word count `avar'
	local nb : word count `bvar'
	local np : word count `pvar'
	local nq : word count `qvar'

	if ("`avar'"~="")|("`bvar'"~="")|("`pvar'"~="")|("`qvar'"~="") {
		if ("`stats'"~= "") {
			noi di as error "stats option not available for model with covariates"
			exit 198
		}
		if ("`poorfrac'"~="") & `poorfrac' > 0   {
			noi di as error "poorfrac option not available for model with covariates"
			exit 198
		}
		if ("`pdf'"~="") {
			noi di as error "pdf option not available for model with covariates"
			exit 198
		}
		if ("`cdf'"~="") {
			noi di as error "cdf option not available for model with covariates"
			exit 198
		}
	}

	if "`cdf'" ~= "" {
		confirm new variable `cdf'
	}
	if "`pdf'" ~= "" {
		confirm new variable `pdf'
	}

	if  (`poorfrac' < 0)  {
		di as error "poorfrac value must be positive"
		exit 198
	}


	local option0 `options'

	local wtype `weight'
	local wtexp `"`exp'"'
	if "`weight'" != "" {
		local wgt `"[`weight'`exp']"'
	}

	if "`weight'" == "pweight" & "`svy'" == "" {
		di as error "To use pweights, -svyset- your data, and use -svy- option"
		exit
	}

	if "`weight'" == "pweight" | "`cluster'" != "" {
	        local robust "robust"
	}

	if "`cluster'" ! = "" {
		local clopt "cluster(`cluster')"
	}

	if "`level'" != "" {
		local level "level(`level')"
	}

        local log = cond("`log'" == "", "noisily", "quietly")

	marksample touse
	markout `touse' `varlist' `avar' `bvar' `pvar' `qvar' `cluster' `censvar'
	if "`svy'" ~= "" {
		svymarkout `touse'
		svyopts modopts diopts options, `option0'

		eret local diopts "`diopts'"
	}
	mlopts mlopts, `options'


	set more off

quietly {

	count if `inc' < 0 & `touse'
	local ct =  r(N)
	if `ct' > 0 {
		noi di " "
		noi di as text "Warning: {res:`inc'} has `ct' values < 0;" _c
		noi di as text " not used in calculations"
	  }

	count if `inc' == 0 & `touse'
	local ct =  r(N)
	if `ct' > 0 {
		noi di " "
		noi di as text "Warning: {res:`inc'} has `ct' values = 0;" _c
		noi di as text " not used in calculations"
	}

	replace `touse' = 0 if `inc' <= 0


        count if `touse'
        if r(N) == 0 {
		error 2000
	}

	if "`censvar'" != "" {
		capture quietly assert (`censvar' == 1 | `censvar' == 0) if `touse'
		if _rc != 0 {
			di as error "censvar must be 0/1 variable"
			exit 198
		}
	}

	if "`lefttr'" != "" replace `touse' = 0 if `inc' < `lefttr'
	global S_mlz "`lefttr'"


	if "`from'" != ""  {
		local b0 "`from'"
	}

	global S_mlinc "`inc'"
	global S_mlcens "`censvar'"

	`log' ml model lf gb2lfit_ll (ln_a: `avar') (ln_b: `bvar') (ln_p: `pvar') (ln_q: `qvar') 	///
		`wgt' if `touse' , maximize 						///
		collinear title(`title') `robust' `svy' init(`b0')		 		///
		search(on) `clopt' `level' `mlopts' `stdopts' `modopts'

	eret local cmd "gb2lfit"
	eret local depvar "`inc'"
	if "`censvar'" != "" eret local censvar "`censvar'"
	if "`lefttr'" != "" eret local lefttr "`lefttr'"


	tempname b ba bb bp bq
	mat `b' = e(b)
	mat `ba' = `b'[1,"ln_a:"]
	mat `bb' = `b'[1,"ln_b:"]
	mat `bp' = `b'[1,"ln_p:"]
	mat `bq' = `b'[1,"ln_q:"]

	eret matrix b_lna = `ba'
	eret matrix b_lnb = `bb'
	eret matrix b_lnp = `bp'
	eret matrix b_lnq = `bq'
	eret scalar length_b_lna = 1 + `na'
	eret scalar length_b_lnb = 1 + `nb'
	eret scalar length_b_lnp = 1 + `np'
	eret scalar length_b_lnq = 1 + `nq'

	if ("`avar'"!="" | "`bvar'"!="" | "`pvar'"!=""  | "`qvar'"!="") {
		eret scalar nocov = 0
	}

	if "`avar'"=="" & "`bvar'"=="" & "`pvar'"=="" & "`qvar'"=="" {

		tempname e v

		mat `e' = e(b)
		local a = exp(`e'[1,1])
		local b = exp(`e'[1,2])
		local p = exp(`e'[1,3])
		local q = exp(`e'[1,4])

		eret scalar ba = `a'
		eret scalar bb = `b'
		eret scalar bp = `p'
		eret scalar bq = `q'

		mat `v' = e(V)
		eret scalar se_a = `a' * sqrt(`v'[1,1])
		eret scalar se_b = `b' * sqrt(`v'[2,2])
		eret scalar se_p = `p' * sqrt(`v'[3,3])
		eret scalar se_q = `q' * sqrt(`v'[4,4])

		eret scalar nocov = 1
		eret local svy "`svy'"

			/* Estimated GB2 c.d.f. */

		if "`cdf'" ~= "" {
			qui ge `cdf' = ibeta(`p',`q', (`inc'/`b')^`a'/(1+(`inc'/`b')^`a') ) if `touse'
		 	eret local cdfvar "`cdf'"
		}


			/* Estimated GB2 p.d.f. */

		if "`pdf'" ~= "" {
		 	qui ge `pdf' = (`a'*(`inc')^(`a'*`p'-1))*exp(lngamma(`p'+`q')) / (		  ///
	 				 (`b'^(`a'*`p'))*exp(lngamma(`p') + lngamma(`q'))  ///
					  *( (1 +(`inc'/`b')^`a')^(`p'+`q') ) ///
					) if `touse'
			eret local pdfvar "`pdf'"
		}


			/* Fraction with income below given level */

		if "`poorfrac'" ~= "" & `poorfrac' > 0 {
			eret scalar poorfrac = ibeta(`p',`q', (`poorfrac'/`b')^`a'/(1+(`poorfrac'/`b')^`a') )
			eret scalar pline = `poorfrac'
		}


			/* selected quantiles predicted from GB2 model */
			/* Lorenz curve ordinates at selected quantiles */

		if "`stats'" ~= "" {
			eret scalar mean = `b'*exp(				///
						lngamma(`p'+1/`a') 		///
						+ lngamma(`q'-1/`a') 		///
						- lngamma(`p')			///
						- lngamma(`q') 			///
						)
			eret scalar mode = cond(`a'*`p'>1,`b'*(((`a'*`p'-1)/(`a'*`q'+1))^(1/`a')),0,.)
			eret scalar var = `b'*`b'*exp(				///
						lngamma(`p'+2/`a') 		///
					   	+ lngamma(`q'-2/`a')		///
					        - lngamma(`p')			///
						- lngamma(`q')			///
						) - (`e(mean)'*`e(mean)')
			eret scalar sd = sqrt(`e(var)')
			eret scalar i2 = .5*`e(var)'/(`e(mean)'*`e(mean)')
			eret scalar gini = .
			// Gini coeff is function of generalized hypergeometric function!!!
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname ib
				scalar `ib' = invibeta(`p',`q',`x'/100)
				eret scalar p`x' =  `b'* (`ib'/(1-`ib'))^(1/`a')
				eret scalar Lp`x' = ibeta(`p'+1/`a',`q'-1/`a',(`e(p`x')'/`b')^`a'/(1+(`e(p`x')'/`b')^`a') )
			}
			eret scalar p90p10 = `e(p90)'/`e(p10)'
			eret scalar p75p25 = `e(p75)'/`e(p25)'

		}

			/* Selected GE inequality indices & other stats */
		if "`extras'" != "" {
			nlcom (aXp:    ( exp([ln_a]_cons) * exp([ln_p]_cons)) ) ///
		          (aXq:    ( exp([ln_a]_cons) * exp([ln_q]_cons) ) )	///
			(mean: exp([ln_b]_cons)*exp(				///
			  lngamma( exp([ln_p]_cons)+1/exp([ln_a]_cons))  	///
			+ lngamma( exp([ln_q]_cons)-1/exp([ln_a]_cons))  	///
			- lngamma( exp([ln_p]_cons))  				///
			- lngamma( exp([ln_q]_cons)) 				///
			))							///
			(var: (exp([ln_b]_cons))^2*exp(				///
			  lngamma( exp([ln_p]_cons)+2/exp([ln_a]_cons))  	///
			+ lngamma( exp([ln_q]_cons)-2/exp([ln_a]_cons))  	///
			- lngamma( exp([ln_p]_cons))  				///
			- lngamma( exp([ln_q]_cons)) 				///
			)							///
			- ( exp([ln_b]_cons)*exp(				///
			  lngamma( exp([ln_p]_cons)+1/exp([ln_a]_cons))  	///
			+ lngamma( exp([ln_q]_cons)-1/exp([ln_a]_cons))  	///
			- lngamma( exp([ln_p]_cons))  				///
			- lngamma( exp([ln_q]_cons)) ))^2			///
			)							///
			(sd: sqrt( (exp([ln_b]_cons))^2*exp(			///
			  lngamma( exp([ln_p]_cons)+2/exp([ln_a]_cons))  	///
			+ lngamma( exp([ln_q]_cons)-2/exp([ln_a]_cons))  	///
			- lngamma( exp([ln_p]_cons))  				///
			- lngamma( exp([ln_q]_cons)) 				///
			)							///
			- ( exp([ln_b]_cons)*exp(				///
			  lngamma( exp([ln_p]_cons)+1/exp([ln_a]_cons))  	///
			+ lngamma( exp([ln_q]_cons)-1/exp([ln_a]_cons))  	///
			- lngamma( exp([ln_p]_cons))  				///
			- lngamma( exp([ln_q]_cons)) ))^2 )			///
			)							///
			(GEm1: -.5 + exp( lngamma(exp([ln_p]_cons)-1/exp([ln_a]_cons)) ///
				+ lngamma(exp([ln_q]_cons)+1/exp([ln_a]_cons)) 	///
				+ lngamma(exp([ln_p]_cons)+1/exp([ln_a]_cons))	///
				+ lngamma(exp([ln_q]_cons)-1/exp([ln_a]_cons))	///
				- 2*lngamma(exp([ln_p]_cons))			///
				- 2*lngamma(exp([ln_q]_cons))  - ln(2) ) 	///
			)							///
			(GE0: lngamma( exp([ln_p]_cons)+1/exp([ln_a]_cons))	///
			+ lngamma( exp([ln_q]_cons)-1/exp([ln_a]_cons))  	///
			- lngamma( exp([ln_p]_cons))  				///
			- lngamma( exp([ln_q]_cons)) 				///
			- digamma(exp([ln_p]_cons))/exp([ln_a]_cons)		///
			+ digamma(exp([ln_q]_cons))/exp([ln_a]_cons)		///
			)							///
			(GE1: lngamma( exp([ln_p]_cons))			///
			+ lngamma( exp([ln_q]_cons)) 				///
			- lngamma( exp([ln_p]_cons)+1/exp([ln_a]_cons))		///
			- lngamma( exp([ln_q]_cons)-1/exp([ln_a]_cons))  	///
			+ digamma(exp([ln_p]_cons)+1/exp([ln_a]_cons))/exp([ln_a]_cons)	///
			- digamma(exp([ln_q]_cons)-1/exp([ln_a]_cons))/exp([ln_a]_cons)	///
			)							///
			(GE2: -.5 + exp( lngamma(exp([ln_p]_cons)) + lngamma(exp([ln_q]_cons)) ///
				+ lngamma(exp([ln_p]_cons)+2/exp([ln_a]_cons))	///
				+ lngamma(exp([ln_q]_cons)-2/exp([ln_a]_cons))	///
				- 2*lngamma(exp([ln_p]_cons)+1/exp([ln_a]_cons))	///
				- 2*lngamma(exp([ln_q]_cons)-1/exp([ln_a]_cons))  - ln(2) ) ///
			)

			tempname bGE VGE
			matrix `VGE' = r(V)
			matrix	`bGE' = r(b)
			ereturn matrix V_GE = `VGE', copy
			ereturn matrix b_GE = `bGE', copy

			ereturn scalar aXp = `bGE'[1,1]
			ereturn scalar aXq = `bGE'[1,2]
			ereturn scalar cmean = `bGE'[1,3]
			ereturn scalar cvar = `bGE'[1,4]
			ereturn scalar csd = `bGE'[1,5]
			ereturn scalar GEm1 = `bGE'[1,6]
			ereturn scalar GE0 = `bGE'[1,7]
			ereturn scalar GE1 = `bGE'[1,8]
			ereturn scalar GE2 = `bGE'[1,9]

			ereturn scalar se_aXp = sqrt(`VGE'[1,1])
			ereturn scalar se_aXq = sqrt(`VGE'[2,2])
			ereturn scalar se_cmean = sqrt(`VGE'[3,3])
			ereturn scalar se_cvar = sqrt(`VGE'[4,4])
			ereturn scalar se_csd = sqrt(`VGE'[5,5])
			ereturn scalar se_GEm1 = sqrt(`VGE'[6,6])
			ereturn scalar se_GE0 = sqrt(`VGE'[7,7])
			ereturn scalar se_GE1 = sqrt(`VGE'[8,8])
			ereturn scalar se_GE2 = sqrt(`VGE'[9,9])

			// estimated quantiles, top shares, etc.

			* overall
			quantcalc `a' `b' `p' `q'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				eret scalar cp`x' = `p`x'est'
				eret scalar cLp`x' = `Lp`x'est'
			}
			eret scalar cp75p25 = `p75p25est'
			eret scalar cp90p10 = `p90p10est'
			eret scalar cshtop01 = `shtop01est'
			eret scalar cshtop05 = `shtop05est'
			eret scalar cshtop10 = `shtop10est'

			* Recalculations for numerical differentiation

			tempname a1 a2 b1 b2 p1 p2 q1 q2
			scalar `a1' = `a' - `a'/100
			scalar `a2' = `a' + `a'/100
			scalar `b1' = `b' - `b'/100
			scalar `b2' = `b' + `b'/100
			scalar `p1' = `p' - `p'/100
			scalar `p2' = `p' + `p'/100
			scalar `q1' = `q' - `q'/100
			scalar `q2' = `q' + `q'/100

			quantcalc `a1' `b' `p' `q'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'a1 cLp`x'a1
				scalar `cp`x'a1' = `p`x'est'
				scalar `cLp`x'a1' = `Lp`x'est'
			}
			tempname cp75p25a1 cp90p10a1 cshtop01a1 cshtop05a1 cshtop10a1
			scalar `cp75p25a1' = `p75p25est'
			scalar `cp90p10a1' = `p90p10est'
			scalar `cshtop01a1' = `shtop01est'
			scalar `cshtop05a1' = `shtop05est'
			scalar `cshtop10a1' = `shtop10est'

			quantcalc `a2' `b' `p' `q'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'a2 cLp`x'a2
				scalar `cp`x'a2' = `p`x'est'
				scalar `cLp`x'a2' = `Lp`x'est'
			}
			tempname cp75p25a2 cp90p10a2 cshtop01a2 cshtop05a2 cshtop10a2
			scalar `cp75p25a2' = `p75p25est'
			scalar `cp90p10a2' = `p90p10est'
			scalar `cshtop01a2' = `shtop01est'
			scalar `cshtop05a2' = `shtop05est'
			scalar `cshtop10a2' = `shtop10est'

				// partial multiplied by parameter since logparameter metric
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname aprimep`x' aprimeLp`x'
				scalar `aprimep`x'' = `a' * (`cp`x'a1' - `cp`x'a2')/(`a'/50)
				scalar `aprimeLp`x'' = `a' * (`cLp`x'a1' - `cLp`x'a2')/(`a'/50)
			}
			tempname aprimep75p25 aprimep90p10 aprimetop01 aprimetop05 aprimetop10
			scalar `aprimep75p25' = `a' * (`cp75p25a1' - `cp75p25a2')/(`a'/50)
			scalar `aprimep90p10' = `a' * (`cp90p10a1' - `cp90p10a2')/(`a'/50)
			scalar `aprimetop01' = `a' * (`cshtop01a1' - `cshtop01a2')/(`a'/50)
			scalar `aprimetop05' = `a' * (`cshtop05a1' - `cshtop05a2')/(`a'/50)
			scalar `aprimetop10' = `a' * (`cshtop10a1' - `cshtop10a2')/(`a'/50)

			quantcalc `a' `b1' `p' `q'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'b1 cLp`x'b1
				scalar `cp`x'b1' = `p`x'est'
				scalar `cLp`x'b1' = `Lp`x'est'
			}
			tempname cp75p25b1 cp90p10b1 cshtop01b1 cshtop05b1 cshtop10b1
			scalar `cp75p25b1' = `p75p25est'
			scalar `cp90p10b1' = `p90p10est'
			scalar `cshtop01b1' = `shtop01est'
			scalar `cshtop05b1' = `shtop05est'
			scalar `cshtop10b1' = `shtop10est'

			quantcalc `a' `b2' `p' `q'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'b2 cLp`x'b2
				scalar `cp`x'b2' = `p`x'est'
				scalar `cLp`x'b2' = `Lp`x'est'
			}
			tempname cp75p25b2 cp90p10b2 cshtop01b2 cshtop05b2 cshtop10b2
			scalar `cp75p25b2' = `p75p25est'
			scalar `cp90p10b2' = `p90p10est'
			scalar `cshtop01b2' = `shtop01est'
			scalar `cshtop05b2' = `shtop05est'
			scalar `cshtop10b2' = `shtop10est'

				// partial multiplied by parameter since logparameter metric
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname bprimep`x' bprimeLp`x'
				scalar `bprimep`x'' = `b' * (`cp`x'b1' - `cp`x'b2')/(`b'/50)
				scalar `bprimeLp`x'' = `b' * (`cLp`x'b1' - `cLp`x'b2')/(`b'/50)
			}
			tempname bprimep75p25 bprimep90p10 bprimetop01 bprimetop05 bprimetop10
			scalar `bprimep75p25' = `b' * (`cp75p25b1' - `cp75p25b2')/(`b'/50)
			scalar `bprimep90p10' = `b' * (`cp90p10b1' - `cp90p10b2')/(`b'/50)
			scalar `bprimetop01' = `b' * (`cshtop01b1' - `cshtop01b2')/(`b'/50)
			scalar `bprimetop05' = `b' * (`cshtop05b1' - `cshtop05b2')/(`b'/50)
			scalar `bprimetop10' = `b' * (`cshtop10b1' - `cshtop10b2')/(`b'/50)

			quantcalc `a' `b' `p1' `q'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'p1 cLp`x'p1
				scalar `cp`x'p1' = `p`x'est'
				scalar `cLp`x'p1' = `Lp`x'est'
			}
			tempname cp75p25p1 cp90p10p1 cshtop01p1 cshtop05p1 cshtop10p1
			scalar `cp75p25p1' = `p75p25est'
			scalar `cp90p10p1' = `p90p10est'
			scalar `cshtop01p1' = `shtop01est'
			scalar `cshtop05p1' = `shtop05est'
			scalar `cshtop10p1' = `shtop10est'

			quantcalc `a' `b' `p2' `q'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'p2 cLp`x'p2
				scalar `cp`x'p2' = `p`x'est'
				scalar `cLp`x'p2' = `Lp`x'est'
			}
			tempname cp75p25p2 cp90p10p2 cshtop01p2 cshtop05p2 cshtop10p2
			scalar `cp75p25p2' = `p75p25est'
			scalar `cp90p10p2' = `p90p10est'
			scalar `cshtop01p2' = `shtop01est'
			scalar `cshtop05p2' = `shtop05est'
			scalar `cshtop10p2' = `shtop10est'

				// partial multiplied by parameter since logparameter metric
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname pprimep`x' pprimeLp`x'
				scalar `pprimep`x'' = `p' * (`cp`x'p1' - `cp`x'p2')/(`p'/50)
				scalar `pprimeLp`x'' = `p' * (`cLp`x'p1' - `cLp`x'p2')/(`p'/50)
			}
			tempname pprimep75p25 pprimep90p10 pprimetop01 pprimetop05 pprimetop10
			scalar `pprimep75p25' = `p' * (`cp75p25p1' - `cp75p25p2')/(`p'/50)
			scalar `pprimep90p10' = `p' * (`cp90p10p1' - `cp90p10p2')/(`p'/50)
			scalar `pprimetop01' = `p' * (`cshtop01p1' - `cshtop01p2')/(`p'/50)
			scalar `pprimetop05' = `p' * (`cshtop05p1' - `cshtop05p2')/(`p'/50)
			scalar `pprimetop10' = `p' * (`cshtop10p1' - `cshtop10p2')/(`p'/50)

			quantcalc `a' `b' `p' `q1'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'q1 cLp`x'q1
				scalar `cp`x'q1' = `p`x'est'
				scalar `cLp`x'q1' = `Lp`x'est'
			}
			tempname cp75p25q1 cp90p10q1 cshtop01q1 cshtop05q1 cshtop10q1
			scalar `cp75p25q1' = `p75p25est'
			scalar `cp90p10q1' = `p90p10est'
			scalar `cshtop01q1' = `shtop01est'
			scalar `cshtop05q1' = `shtop05est'
			scalar `cshtop10q1' = `shtop10est'

			quantcalc `a' `b' `p' `q2'
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname cp`x'q2 cLp`x'q2
				scalar `cp`x'q2' = `p`x'est'
				scalar `cLp`x'q2' = `Lp`x'est'
			}
			tempname cp75p25q2 cp90p10q2 cshtop01q2 cshtop05q2 cshtop10q2
			scalar `cp75p25q2' = `p75p25est'
			scalar `cp90p10q2' = `p90p10est'
			scalar `cshtop01q2' = `shtop01est'
			scalar `cshtop05q2' = `shtop05est'
			scalar `cshtop10q2' = `shtop10est'

				// partial multiplied by parameter since logparameter metric
			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				tempname qprimep`x' qprimeLp`x'
				scalar `qprimep`x'' = `q' * (`cp`x'q1' - `cp`x'q2')/(`q'/50)
				scalar `qprimeLp`x'' = `q' * (`cLp`x'q1' - `cLp`x'q2')/(`q'/50)
			}
			tempname qprimep75p25 qprimep90p10 qprimetop01 qprimetop05 qprimetop10
			scalar `qprimep75p25' = `q' * (`cp75p25q1' - `cp75p25q2')/(`q'/50)
			scalar `qprimep90p10' = `q' * (`cp90p10q1' - `cp90p10q2')/(`q'/50)
			scalar `qprimetop01' = `q' * (`cshtop01q1' - `cshtop01q2')/(`q'/50)
			scalar `qprimetop05' = `q' * (`cshtop05q1' - `cshtop05q2')/(`q'/50)
			scalar `qprimetop10' = `q' * (`cshtop10q1' - `cshtop10q2')/(`q'/50)

			tempname prime V var
			mat V = e(V)

			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				mat `prime' = [`aprimep`x'', `bprimep`x'', `pprimep`x'', `qprimep`x'']
				mat varg = `prime'*V*`prime''
				eret scalar se_cp`x' = sqrt(varg[1,1])
				mat `prime' = [`aprimeLp`x'', `bprimeLp`x'', `pprimeLp`x'', `qprimeLp`x'']
				mat varg = `prime'*V*`prime''
				eret scalar se_cLp`x' = sqrt(varg[1,1])
			}

			mat `prime' = [`aprimep75p25', `bprimep75p25', `pprimep75p25', `qprimep75p25']
			mat varg = `prime'*V*`prime''
			eret scalar se_cp75p25 = sqrt(varg[1,1])

			mat `prime' = [`aprimep90p10', `bprimep90p10', `pprimep90p10', `qprimep90p10']
			mat varg = `prime'*V*`prime''
			eret scalar se_cp90p10 = sqrt(varg[1,1])

			mat `prime' = [`aprimetop01', `bprimetop01', `pprimetop01', `qprimetop01']
			mat varg = `prime'*V*`prime''
			eret scalar se_cshtop01 = sqrt(varg[1,1])

			mat `prime' = [`aprimetop05', `bprimetop05', `pprimetop05', `qprimetop05']
			mat varg = `prime'*V*`prime''
			eret scalar se_cshtop05 = sqrt(varg[1,1])

			mat `prime' = [`aprimetop10', `bprimetop10', `pprimetop10', `qprimetop10']
			mat varg = `prime'*V*`prime''
			eret scalar se_cshtop10 = sqrt(varg[1,1])

		}


		if "`cgini'" != "" {

			capture assert [ln_a]_cons + [ln_q]_cons > 0
			if _rc != 0 {
				noi di " "
				noi di as error "a*q <= 1. 3F2 function is non-convergent. Gini not calculated."
				local jump 1
			}
		}

		if "`cgini'" != "" & "`jump'" == "" {

			noi di " "
			noi di as txt "Calculating Gini using 3F2 functions - be patient"
			noi di as txt "	Convergence tolerance for 3F2 function = " as res `epsilon'

			// estimated Gini
			* overall Gini
			ginicalc `a' `b' `p' `q' `epsilon'
			eret scalar cgini = `giniest'
			eret scalar num_iter = `num_iter'

			tempname a1 a2 p1 p2 q1 q2
			tempname cginia1 cginia2 cginip1 cginip2 cginiq1 cginiq2
			tempname aprime pprime qprime

			* Ginis used for numerical differentiation
			scalar `a1' = `a' - `a'/100
			scalar `a2' = `a' + `a'/100
			ginicalc `a1' `b' `p' `q'  `epsilon'
			scalar `cginia1' = `giniest'
			ginicalc `a2' `b' `p' `q'  `epsilon'
			scalar `cginia2' = `giniest'
				// partial multiplied by parameter since logparameter metric
			scalar `aprime' = `a' * (`cginia1' - `cginia2')/(`a'/50)

			scalar `p1' = `p' - `p'/100
			scalar `p2' = `p' + `p'/100
			ginicalc `a' `b' `p1' `q'  `epsilon'
			scalar `cginip1' = `giniest'
			ginicalc `a' `b' `p2' `q'  `epsilon'
			scalar `cginip2' = `giniest'
				// partial multiplied by parameter since logparameter metric
			scalar `pprime' = `p' * (`cginip1' - `cginip2')/(`p'/50)

			scalar `q1' = `q' - `q'/100
			scalar `q2' = `q' + `q'/100
			ginicalc `a' `b' `p' `q1'  `epsilon'
			scalar `cginiq1' = `giniest'
			ginicalc `a' `b' `p' `q2'  `epsilon'
			scalar `cginiq2' = `giniest'
				// partial multiplied by parameter since logparameter metric
			scalar `qprime' = `q' * (`cginiq1' - `cginiq2')/(`q'/50)

			tempname gprime V varg


			mat `gprime' = [`aprime', 0, `pprime', `qprime']
			mat V = e(V)
			mat varg = `gprime'*V*`gprime''
			eret scalar se_cgini = sqrt(varg[1,1])

		}


		if "`igini'" != "" {


			noi di " "
			noi di as txt "Calculating Gini by numerical integration. # integration points = " as res `nips'


			// estimated Gini

			tempvar s
			range `s' 0 1 `nips'
			eret scalar nips = `nips'
			global S_s "`s'"

			* overall Gini
			giniint `a' `b' `p' `q'  `nips'
			eret scalar igini = `intgini'

			tempname a1 a2 p1 p2 q1 q2
			tempname iginia1 iginia2 iginip1 iginip2 iginiq1 iginiq2
			tempname aprime pprime qprime

			* Ginis used for numerical differentiation
			scalar `a1' = `a' - `a'/100
			scalar `a2' = `a' + `a'/100
			giniint `a1' `b' `p' `q' `nips'
			scalar `iginia1' = `intgini'
			giniint `a2' `b' `p' `q' `nips'
			scalar `iginia2' = `intgini'
				// partial multiplied by parameter since logparameter metric
			scalar `aprime' = `a' * (`iginia1' - `iginia2')/(`a'/50)

			scalar `p1' = `p' - `p'/100
			scalar `p2' = `p' + `p'/100
			giniint `a' `b' `p1' `q'  `nips'
			scalar `iginip1' = `intgini'
			giniint `a' `b' `p2' `q'  `nips'
			scalar `iginip2' = `intgini'
				// partial multiplied by parameter since logparameter metric
			scalar `pprime' = `p' * (`iginip1' - `iginip2')/(`p'/50)

			scalar `q1' = `q' - `q'/100
			scalar `q2' = `q' + `q'/100
			giniint `a' `b' `p' `q1'  `nips'
			scalar `iginiq1' = `intgini'
			giniint `a' `b' `p' `q2'  `nips'
			scalar `iginiq2' = `intgini'
				// partial multiplied by parameter since logparameter metric
			scalar `qprime' = `q' * (`iginiq1' - `iginiq2')/(`q'/50)

			tempname gprime V varg


			mat `gprime' = [`aprime', 0, `pprime', `qprime']
			mat V = e(V)
			mat varg = `gprime'*V*`gprime''
			eret scalar se_igini = sqrt(varg[1,1])

					// get rid of global now
			global S_s

		}


	} // end block executed if no covariates


} // end -quietly block-


	if "`poorfrac'" ~= "" & `poorfrac' > 0 {
		local pfrac "poorfrac(`poorfrac')"
	}


	Display, `level' `pfrac' `diopts'

end




program define Display

	syntax [,Level(int $S_level) POORfrac(real 0)  *]
	local diopts "`options'"
	if `level' < 10 | `level' > 99 {
		local level = 95
	}

	local plus "plus"
	if "`e(svy)'" ~= ""  {
		local plus ""
	}

	if `e(nocov)' == 0 {
		ml display, level(`level') `diopts'
	}

	if `e(nocov)' == 1 {
		ml display, level(`level') `diopts' `plus'
		_diparm ln_a, level(`level') exp prob label("Parameters a")
		di in smcl in gr "{hline 13}{c +}{hline 64}"
		_diparm ln_b, level(`level') exp prob label("b")
		di in smcl in gr "{hline 13}{c +}{hline 64}"
		_diparm ln_p, level(`level') exp prob label("p")
		di in smcl in gr "{hline 13}{c +}{hline 64}"
		_diparm ln_q, level(`level') exp prob label("q")
		di in smcl in gr "{hline 13}{c BT}{hline 64}"
	}

	if  (`poorfrac' < 0)  {
		di as error "poorfrac value must be positive"
		exit
	}

	if `poorfrac' > 0 & (`e(nocov)' == 0) {
		di as error "poorfrac option not available (model was specified with covariates)"
		exit
	}

	if `poorfrac' > 0 & (`e(nocov)' == 1) {
		di " "
		di "Fraction with {res: `e(depvar)'} < `poorfrac' = " as res %9.5f ibeta(`e(bp)',`e(bq)', (`poorfrac'/`e(bb)')^`e(ba)'/(1+(`poorfrac'/`e(bb)')^`e(ba)') )
		di " "
	}


	if "`e(mean)'" != ""  {

		di as txt "{hline 60}"
		di as txt  _col(6) "Quantiles" _col(22) "Cumulative shares of"
		di as txt  _col(22) "total `e(depvar)' (Lorenz ordinates)"
		di as txt "{hline 60}"
		di as txt  " 1%" _col(6) as res %9.5f `e(p1)' _col(20) %9.5f `e(Lp1)'
		di as txt  " 5%" _col(6) as res %9.5f `e(p5)' _col(20) %9.5f `e(Lp5)'
		di as txt  "10%" _col(6) as res %9.5f `e(p10)' _col(20) %9.5f `e(Lp10)'
		di as txt  "20%" _col(6) as res %9.5f `e(p20)' _col(20) %9.5f `e(Lp20)'
		di as txt  "25%" _col(6) as res %9.5f `e(p25)' _col(20) %9.5f `e(Lp25)'
		di as txt  "30%" _col(6) as res %9.5f `e(p30)' _col(20) %9.5f `e(Lp30)'
		di as txt  "40%" _col(6) as res %9.5f `e(p40)' _col(20) %9.5f `e(Lp40)' _c
		di as txt  _col(30) "Mode" _col(42) as res %9.5f `e(mode)'
		di as txt  "50%" _col(6) as res %9.5f `e(p50)' _col(20) %9.5f `e(Lp50)' _c
		di as txt  _col(30) "Mean" _col(42) as res %9.5f `e(mean)'
		di as txt  "60%" _col(6) as res %9.5f `e(p60)' _col(20) %9.5f `e(Lp60)' _c
		di as txt  _col(30) "Std. Dev." _col(42) as res %9.5f `e(sd)'
		di as txt  "70%" _col(6) as res %9.5f `e(p70)' _col(20) %9.5f `e(Lp70)'
		di as txt  "75%" _col(6) as res %9.5f `e(p75)' _col(20) %9.5f `e(Lp75)' _c
		di as txt  _col(30) "Variance" _col(42) as res %9.5f `e(var)'
		di as txt  "80%" _col(6) as res %9.5f `e(p80)' _col(20) %9.5f `e(Lp80)' _c
		di as txt  _col(30) "Half CV^2" _col(42) as res %9.5f `e(i2)'
		di as txt  "90%" _col(6) as res %9.5f `e(p90)' _col(20) %9.5f `e(Lp90)'
*		di as txt  "90%" _col(6) as res %9.5f `e(p90)' _col(20) %9.5f `e(Lp90)' _c
*		di as txt  _col(30) "Gini coeff." _col(42) as res %9.5f `e(gini)'
		di as txt  "95%" _col(6) as res %9.5f `e(p95)' _col(20) %9.5f `e(Lp95)' _c
		di as txt  _col(30) "p90/p10" _col(42) as res %9.5f `e(p90)'/`e(p10)'
		di as txt  "99%" _col(6) as res %9.5f `e(p99)' _col(20) %9.5f `e(Lp99)' _c
		di as txt  _col(30) "p75/p25" _col(42) as res %9.5f `e(p75)'/`e(p25)'
		di as txt "{hline 60}"
	}

	if "`e(cp10)'" != "" {

		di as text "{hline 9}{c TT}{hline 65}"
		di as text "Statistic"_col(10) as text "{c |}   Estimate " _c
		di as text _col(25)  "Std. Err." _col(39) "z"  _c
		di as text _col(46) "P>|z|"  _col(56) "[`level'% Conf. Interval]"
		di as text "{hline 9}{c +}{hline 65}"

		local stats "aXp aXq cmean cvar csd"
		local stats "`stats' cp1 cp5 cp10 cp20 cp25 cp30 cp40 cp50 cp60 cp70 cp75 cp80 cp90 cp95 cp99"
		local stats "`stats' cLp1 cLp5 cLp10 cLp20 cLp25 cLp30 cLp40 cLp50 cLp60 cLp70 cLp75 cLp80 cLp90 cLp95 cLp99"
		local stats "`stats' cshtop01 cshtop05 cshtop10"
		local stats "`stats' cp75p25 cp90p10 GEm1 GE0 GE1 GE2"
		local sname "a*p a*q mean var sd"
		local sname "`sname' p1 p5 p10 p20 p25 p30 p40 p50 p60 p70 p75 p80 p90 p95 p99"
		local sname "`sname' Lp1 Lp5 Lp10 Lp20 Lp25 Lp30 Lp40 Lp50 Lp60 Lp70 Lp75 Lp80 Lp90 Lp95 Lp99"
		local sname "`sname' shTop01pc shTop05pc shTop10pc p75:p25 p90:p10"
		local sname "`sname' GE(-1) GE(0) GE(1) GE(2)"

		if ("`e(cgini)'" != "") & ("`e(igini)'" == "") {
			local stats "`stats' cgini"
			local sname "`sname' cGini"
		}

		if ("`e(cgini)'" == "") & ("`e(igini)'" != "") {
			local stats "`stats' igini"
			local sname "`sname' iGini"
		}

		if ("`e(cgini)'" != "") & ("`e(igini)'" != "") {
			local stats "`stats' cgini igini"
			local sname "`sname' cGini iGini"
		}

		local i = 1
		foreach j of local stats {
			local n: word `i' of `sname'
			di as txt "`n'"  _col(10) "{c |} " as res %10.6f e(`j') _c
			di as res _col(24) %9.6f e(se_`j')  _col(31)  _c
			di %9.3f as result e(`j')/e(se_`j') _col(40) _c
			di %9.3f as result 2*(1-normal(`e(`j')'/`e(se_`j')')) _c
			di _col(56) %9.6g as result `e(`j')'+invnormal((100-`level')/200)*`e(se_`j')' _c
			di _col(67) %9.6g as result `e(`j')'-invnormal((100-`level')/200)*`e(se_`j')'
			local i = `i' + 1
		}
		di as text "{hline 9}{c BT}{hline 65}"
	}

	if "`e(cp10)'" == "" & "`e(cgini)'" != "" & "`e(igini)'" == "" {

		di as text "{hline 9}{c TT}{hline 65}"
		di as text "Statistic"_col(10) as text "{c |}   Estimate " _c
		di as text _col(25)  "Std. Err." _col(39) "z"  _c
		di as text _col(46) "P>|z|"  _col(56) "[`level'% Conf. Interval]"
		di as text "{hline 9}{c +}{hline 65}"
		di as text " cGini" _col(10) "{c |} "  as result %10.7f e(cgini) _c
		di as result _col(24) e(se_cgini) _col(31)  _c
		di %9.3f as result e(cgini)/e(se_cgini) _col(40) _c
		di %9.3f as result 2*(1-normal(`e(cgini)'/`e(se_cgini)')) _c
		di _col(56) %9.0g as result e(cgini)+invnormal((100-`level')/200)*e(se_cgini) _c
		di _col(67) %9.0g as result e(cgini)-invnormal((100-`level')/200)*e(se_cgini)
		di as text "{hline 9}{c BT}{hline 65}"
	}

	if "`e(cp10)'" == "" & "`e(cgini)'" == "" & "`e(igini)'" != "" {

		di as text "{hline 9}{c TT}{hline 65}"
		di as text "Statistic"_col(10) as text "{c |}   Estimate " _c
		di as text _col(25)  "Std. Err." _col(39) "z"  _c
		di as text _col(46) "P>|z|"  _col(56) "[`level'% Conf. Interval]"
		di as text "{hline 9}{c +}{hline 65}"
		di as text " iGini" _col(10) "{c |} "  as result %10.7f e(igini) _c
		di as result _col(24) e(se_igini) _col(31)  _c
		di %9.3f as result e(igini)/e(se_igini) _col(40) _c
		di %9.3f as result 2*(1-normal(`e(igini)'/`e(se_igini)')) _c
		di _col(56) %9.0g as result e(igini)+invnormal((100-`level')/200)*e(se_igini) _c
		di _col(67) %9.0g as result e(igini)-invnormal((100-`level')/200)*e(se_igini)
		di as text "{hline 9}{c BT}{hline 65}"
	}

	if "`e(cp10)'" == "" & "`e(cgini)'" != "" & "`e(igini)'" != "" {

		di as text "{hline 9}{c TT}{hline 65}"
		di as text "Statistic"_col(10) as text "{c |}   Estimate " _c
		di as text _col(25)  "Std. Err." _col(39) "z"  _c
		di as text _col(46) "P>|z|"  _col(56) "[`level'% Conf. Interval]"
		di as text "{hline 9}{c +}{hline 65}"
		di as text " cGini" _col(10) "{c |} "  as result %10.7f e(cgini) _c
		di as result _col(24) e(se_cgini) _col(31)  _c
		di %9.3f as result e(cgini)/e(se_cgini) _col(40) _c
		di %9.3f as result 2*(1-normal(`e(cgini)'/`e(se_cgini)')) _c
		di _col(56) %9.0g as result e(cgini)+invnormal((100-`level')/200)*e(se_cgini) _c
		di _col(67) %9.0g as result e(cgini)-invnormal((100-`level')/200)*e(se_cgini)
		di as text " iGini" _col(10) "{c |} "  as result %10.7f e(igini) _c
		di as result _col(24) e(se_igini) _col(31)  _c
		di %9.3f as result e(igini)/e(se_igini) _col(40) _c
		di %9.3f as result 2*(1-normal(`e(igini)'/`e(se_igini)')) _c
		di _col(56) %9.0g as result e(igini)+invnormal((100-`level')/200)*e(se_igini) _c
		di _col(67) %9.0g as result e(igini)-invnormal((100-`level')/200)*e(se_igini)
		di as text "{hline 9}{c BT}{hline 65}"
	}


end

program define ginicalc

	args ca cb cp cq eps

	// estimated Gini a la McDonald (Econometrica 1984) formula
	//   	= B * [ (1/p)* 3F2A - (1/(p + 1/a))* 3F2B]

	tempname B f3F2A f3F2B sumA sumB sumAlast sumBlast ginicalc

	scalar `B' = exp( lngamma(2*`cq' - (1/`ca')) + lngamma(2*`cp' + (1/`ca'))  ///
			- lngamma(2*`cp' + 2*`cq') + 2*lngamma(`cp' + `cq')	///
			- lngamma(`cp') - lngamma(`cq') 			///
			- lngamma(`cp' + (1/`ca')) - lngamma(`cq' - (1/`ca'))	///
			)

	scalar `sumA' = 0
	scalar `sumB' = 0
	scalar `sumAlast' = .
	scalar `sumBlast' = .

		// 3F2 functions are infinite series;
		//    calculate to pre-specified # decimal places
	local i = 1
	while ( abs(`sumA' - `sumAlast' ) > `eps' ) | ( abs(`sumB' - `sumBlast' ) > `eps' ) {

		scalar `sumAlast' = `sumA'
		scalar `sumA' = `sumA' + 				///
				exp( lngamma(1 + `i') 			///
				+ lngamma(`cp' + `cq' + `i')  		///
				+ lngamma(2*`cp' + (1/`ca') + `i') 	///
				- lngamma(`cp' + 1 + `i')		///
				- lngamma(2*`cp' + 2*`cq' + `i')	///
				- lnfact(`i')				///
				)
		scalar `sumBlast' = `sumB'
		scalar `sumB' = `sumB' + 				///
				exp( lngamma(1 + `i') 			///
				+ lngamma(`cp' + `cq' + `i')  		///
				+ lngamma(2*`cp' + (1/`ca') + `i') 	///
				- lngamma(`cp' + (1/`ca') + 1 + `i')	///
				- lngamma(2*`cp' + 2*`cq' + `i')	///
				- lnfact(`i')				///
				)

		local i = `i' + 1
	}

	c_local num_iter = `i'

	scalar `f3F2A' = 1 +  exp( lngamma(`cp' + 1) + lngamma(2*`cp' + 2*`cq') ///
				- lngamma(`cp' + `cq') - lngamma(2*`cp' + (1/`ca')) ///
				) * `sumA'

	scalar `f3F2B' = 1 +  exp( lngamma(`cp' + (1/`ca') + 1) + lngamma(2*`cp' + 2*`cq') ///
				- lngamma(`cp' + `cq') - lngamma(2*`cp' + (1/`ca')) ///
				) * `sumB'

	scalar `ginicalc' = `B' * ( `f3F2A'/`cp' - `f3F2B'/(`cp' + (1/`ca')) )

	c_local giniest = `ginicalc'

end


program define quantcalc
	args ca cb cp cq

	local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
	foreach x of local ptile {
		tempname ib p`x' Lp`x'
		scalar `ib' = invibeta(`cp',`cq',`x'/100)
		scalar `p`x'' =  `cb'* (`ib'/(1-`ib'))^(1/`ca')
		scalar `Lp`x'' = ibeta(`cp'+1/`ca',`cq'-1/`ca',(`p`x''/`cb')^`ca'/(1+(`p`x''/`cb')^`ca') )

		c_local p`x'est = `p`x''
		c_local Lp`x'est = `Lp`x''

	}

	tempname p90p10 p75p25 shtop01 shtop05 shtop10

	scalar `p90p10' = `p90'/`p10'
	scalar `p75p25' = `p75'/`p25'
	scalar `shtop01' = 1 - `Lp99'
	scalar `shtop05' = 1 - `Lp95'
	scalar `shtop10' = 1 - `Lp90'
	c_local p90p10est = `p90p10'
	c_local p75p25est = `p75p25'
	c_local shtop01est = `shtop01'
	c_local shtop05est = `shtop05'
	c_local shtop10est = `shtop10'

end




program define giniint

	args ca cb cp cq cnips

	// estimated Gini using Stata -intreg-  Gini = int(0,1)  { 2*(s-L(s) } ds
	//	Calculate over range of p = F(y), i.e. [0,1]
	// # integration points in interval is set by NIPS options

	tempvar z L y
	ge double `z' =  `cb'* ( invibeta(`cp',`cq',$S_s) / ( 1 - invibeta(`cp',`cq', $S_s) ))^(1/`ca')
	ge double `L' = ibeta(`cp'+1/`ca',`cq'-1/`ca',(`z'/`cb')^`ca'/(1+ (`z'/`cb')^`ca') )
	replace `L' = 0 in 1
	replace `L' = 1 in `cnips'

	ge double `y' = 2*($S_s - `L')
	integ `y' $S_s

	c_local intgini = r(integral)

end