*! paretofit.ado, 1.1.0, Stephen P. Jenkins, 2015-10-15
*! paretofit.ado, 1.0.0, Stephen P. Jenkins & Philippe Van Kerm, 2007-04-06

*----------------------------------------------------------------------
*                                       Revision history :
* 1.1.0: added more output from stats option; fixed bug in variance calculation
*----------------------------------------------------------------------
*----------------------------------------------------------------------
*                                       Description :
*
* Fit Pareto distribution by ML to unit record data.
* Code largely borrowed from -smfit- (for fitting Singh-Maddala distributions).
* See similar commands -dagumfit-, -gb2fit-, -lognfit-.
*
*----------------------------------------------------------------------


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

program define paretofit, eclass byable(onecall) properties(ml_score svyr)
	version 9.2
	if replay() {
		if "`e(cmd)'" != "paretofit" {
			noi di as error "results for {cmd:paretofit} 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) sortpreserve

	syntax varlist(max=1) [if] [in] [aw fw pw iw] [,  ///
		Avar(varlist numeric) ///
		X0(string) ME(string) ///
		CDF(namelist max=1) PDF(namelist max=1) ///
		POORfrac(real 0) STats ///
		Robust Cluster(varname) From(real 0)  ///
		Level(cilevel)    ///
		noLOG * ]

	local title "ML fit of Pareto distribution"

	local inc "`varlist'"

	local na : word count `avar'

	if ("`avar'"~="") {
		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
	}

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

	local option0 `options'

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

	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'
	markout `touse' `cluster' , strok
	mlopts mlopts, `options'

	if ! `:length local x0'  {
		qui su `inc' if `touse', meanonly
		loc x0 = r(min)
	}

	if  (`x0' < 0)  {
		di as error "x0 must be positive"
		exit 198
	}
	if  (`poorfrac' != 0) & (`poorfrac' < `x0')  {
		di as error "poorfrac value must be higher than x0 (here `x0')"
		exit 198
	}

	if ! `:length local me'  {
		loc t = `x0'
	}
	else local t = `me'
	if  `t' < `x0'   {
		di as error "ME threshold cannot be less than x0 (here `x0')"
		exit 198
	}


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

	qui count if `touse'
	if (r(N) == 0) 	error 2000
	else			loc nobs = r(N)

	if (`from' == 0) 	local init ""
	else 				local init "init(`from' , copy)"


	global S_mlinc "`inc'"
	global S_x0 "`x0'"

	`log' ml model lf paretofit_ll (a: `avar')  	///
		`wgt' if `touse' , maximize					///
		collinear title(`title') `robust' 	///
		search(on) `clopt' `level' `mlopts' `stdopts' `modopts' `init'

	eret local cmd "paretofit"
	eret local depvar "`inc'"
	eret scalar x0 = `x0'
	eret scalar t = `t'

	if ("`avar'"!="") 	eret scalar nocov = 0
	else {
		tempname e
		mat `e' = e(b)
		local a = `e'[1,1]
		eret scalar ba = `a'
		eret scalar nocov = 1

		/* Estimated Pareto c.d.f. */
		if "`cdf'" ~= "" {
			qui ge `cdf' = 1 - (`inc'/`x0')^(-`a') if `touse'
		 	eret local cdfvar "`cdf'"
		}
		/* Estimated Pareto p.d.f. */
		if "`pdf'" ~= "" {
		 	qui ge `pdf' = `a' * (`inc')^(-(`a'+1)) * (`x0')^(`a') if `touse'
			eret local pdfvar "`pdf'"
		}
		/* Summary statistics and Lorenz ordinates */
		if "`stats'" ~= "" {
			eret scalar mode = `x0'

			quietly nlcom 	(alpha: (_b[a:_cons]) ) 			///
			(beta: ( _b[a:_cons]/(_b[a:_cons] - 1) ) ) 			///
			(gini: ( 1/ (2*_b[a:_cons] - 1 ) )  )				///
			(mld: ( log(_b[a:_cons]/(_b[a:_cons] - 1)) - (1/_b[a:_cons])  ) )	///
			(theil: (  (1/(_b[a:_cons] - 1)) - log(_b[a:_cons]/(_b[a:_cons] - 1)) ) ) ///
			(mean:  ( `x0'*_b[a:_cons]/(_b[a:_cons] - 1)  ) )			///
			(me: ( `t'/(_b[a:_cons] - 1)  ) ), noheader

			tempname b V
			mat `b' = r(b)
			mat `V' = r(V)
			ereturn scalar alpha = `b'[1,1]
			ereturn scalar se_alpha = sqrt( `V'[1,1] )
			ereturn scalar beta = `b'[1,2]
			ereturn scalar se_beta = sqrt( `V'[2,2] )
			ereturn scalar gini = `b'[1,3]
			ereturn scalar se_gini = sqrt( `V'[3,3] )
			ereturn scalar ge0 = `b'[1,4]
			ereturn scalar se_ge0 = sqrt( `V'[4,4] )
			ereturn scalar ge1 = `b'[1,5]
			ereturn scalar se_ge1 = sqrt( `V'[5,5] )
			ereturn scalar mean = `b'[1,6]
			ereturn scalar se_mean = sqrt( `V'[6,6] )
			ereturn scalar me = `b'[1,7]
			ereturn scalar se_me = sqrt( `V'[7,7] )
				// SPJ changed `a' >=2 to `a' > 2 (see K & K)
			if (`a' > 2) {

				quietly nlcom 	(ge2: ( 1/ ( 2*_b[a:_cons]*(_b[a:_cons]-2) )) ) 	 ///
				 (var: ( (`x0'^2)*_b[a:_cons]/( (_b[a:_cons]-2)*(_b[a:_cons]-1)^2) )  )   ///
				 (sd: (  sqrt((`x0'^2)*_b[a:_cons]/( (_b[a:_cons]-2)*(_b[a:_cons]-1)^2) )) ) , noheader

				tempname c W
				mat `c' = r(b)
				mat `W' = r(V)
				ereturn scalar ge2 = `c'[1,1]
				ereturn scalar se_ge2 = sqrt( `V'[1,1] )
				ereturn scalar var = `c'[1,2]
				ereturn scalar se_var = sqrt( `V'[2,2] )
				ereturn scalar sd = `c'[1,3]
				ereturn scalar se_sd = sqrt( `V'[3,3] )
			}

			local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99"
			foreach x of local ptile {
				eret scalar p`x' = `x0'*(1-`x'/100)^(-(1/`a'))
				eret scalar Lp`x' = 1 - (1-`x'/100)^(1 - (1/`a'))
			}
		}
		/* Fraction with income below given level */
		if ("`poorfrac'" ~= "" & `poorfrac' > 0) {
			eret scalar poorfrac = 1 - ( `x0' / `poorfrac' )^`a'
			eret scalar pline = `poorfrac'
			local pfrac "poorfrac(`poorfrac')"
		}
	}

	Display, `level' `pfrac'

end

program define Display
	syntax [,Level(cilevel) POORfrac(real 0) *]
	local diopts "`options'"
	ml display, level(`level') `diopts'
	if `level' < 10 | `level' > 99 {
		local level = 95
	}

	if 	"`e(mean)'" ~= ""  {
		local sname "alpha beta mean me gini ge0 ge1"
		local i = 1
		foreach j of local sname {
				local n: word `i' of `sname'
				di as txt "   `n'"  _col(14) "{c |} " as res %10.6f e(`j') _c
				di as res _col(28) %9.6f e(se_`j')  _col(34)  _c
				di %9.3f as result e(`j')/e(se_`j') _col(44) _c
				di %9.3f as result 2*(1-normal(`e(`j')'/`e(se_`j')')) _c
				di _col(58) %9.6g as result `e(`j')'+invnormal((100-`level')/200)*`e(se_`j')' _c
				di _col(70) %9.6g as result `e(`j')'-invnormal((100-`level')/200)*`e(se_`j')'
				local i = `i' + 1
		}
			if `e(alpha)' > 2 {
				local vname = "ge2 var sd"
				local k = 1
				foreach m of local vname {
					local n: word `k' of `vname'
					di as txt "   `n'"  _col(14) "{c |} " as res %10.6f e(`m') _c
					di as res _col(28) %9.6f e(se_`m')  _col(34)  _c
					di %9.3f as result e(`m')/e(se_`m') _col(44) _c
					di %9.3f as result 2*(1-normal(`e(`m')'/`e(se_`m')')) _c
					di _col(58) %9.6g as result `e(`m')'+invnormal((100-`level')/200)*`e(se_`m')' _c
					di _col(70) %9.6g as result `e(`m')'-invnormal((100-`level')/200)*`e(se_`m')'
					local k = `k' + 1
			}

		}

		di as text "{hline 13}{c BT}{hline 65}"
	}

	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.5g 1 - ( `e(x0)' / `poorfrac' )^`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.5g `e(p1)' _col(20) %9.5g `e(Lp1)'
		di as txt  " 5%" _col(6) as res %9.5g `e(p5)' _col(20) %9.5g `e(Lp5)'
		di as txt  "10%" _col(6) as res %9.5g `e(p10)' _col(20) %9.5g `e(Lp10)'
		di as txt  "20%" _col(6) as res %9.5g `e(p20)' _col(20) %9.5g `e(Lp20)'
		di as txt  "25%" _col(6) as res %9.5g `e(p25)' _col(20) %9.5g `e(Lp25)'
		di as txt  "30%" _col(6) as res %9.5g `e(p30)' _col(20) %9.5g `e(Lp30)'
		di as txt  "40%" _col(6) as res %9.5g `e(p40)' _col(20) %9.5g `e(Lp40)' _c
		di as txt  _col(30) "Mode" _col(42) as res %9.5g `e(mode)'
		di as txt  "50%" _col(6) as res %9.5g `e(p50)' _col(20) %9.5g `e(Lp50)' _c
		di as txt  _col(30) "Mean" _col(42) as res %9.5g `e(mean)'
		di as txt  "60%" _col(6) as res %9.5g `e(p60)' _col(20) %9.5g `e(Lp60)' _c
		di as txt  _col(30) "Std. Dev." _col(42) as res %9.5g `e(sd)'
		di as txt  "70%" _col(6) as res %9.5g `e(p70)' _col(20) %9.5g `e(Lp70)'
		di as txt  "75%" _col(6) as res %9.5g `e(p75)' _col(20) %9.5g `e(Lp75)' _c
		di as txt  _col(30) "Variance" _col(42) as res %9.5g `e(var)'
		di as txt  "80%" _col(6) as res %9.5g `e(p80)' _col(20) %9.5g `e(Lp80)' _c
		di as txt  _col(30) "Half CV^2" _col(42) as res %9.5g `e(i2)'
		di as txt  "90%" _col(6) as res %9.5g `e(p90)' _col(20) %9.5g `e(Lp90)' _c
		di as txt  _col(30) "Gini coeff." _col(42) as res %9.5g `e(gini)'
		di as txt  "95%" _col(6) as res %9.5g `e(p95)' _col(20) %9.5g `e(Lp95)' _c
		di as txt  _col(30) "p90/p10" _col(42) as res %9.5g `e(p90)'/`e(p10)'
		di as txt  "99%" _col(6) as res %9.5g `e(p99)' _col(20) %9.5g `e(Lp99)' _c
		di as txt  _col(30) "p75/p25" _col(42) as res %9.5g `e(p75)'/`e(p25)'
		di as txt "{hline 60}"


	}


end



exit