*! version 1.0.0 Stephen P. Jenkins, Dec 1997    STB-48 sg106
*! Fit Dagum distribution by ML to unit record data

program define dagumfit
	version 4.0

/*------------------------------------------------ playback request */
  local options "LEvel(integer $S_level) "
  if "`*'" == "" | substr("`1'",1,1) == "," {
	if "$S_E_cmd" ~= "dagumfit" {
	error 301
	}
  parse "`*'"
  }

/*------------------------------------------------ estimation */
  else {

	local varlist "req ex max(1)"
	local if "opt"
	local in "opt"
	local options "TRace noLOG Stats"
	local options "`options' b0(real 54.59815) d0(real 1.1051709)"
	local options "`options' h0(real 442414)"
	local options "`options' LEvel(integer $S_level)"
	local options "`options' cdf(string) pdf(string)"
	local weight "aweight fweight"
	parse "`*'"
	parse "`varlist'", parse (" ")
	local inc "`1'"

	tempvar wi touse badinc mysamp
	tempname c v f V b_se d_se h_se 

	if (`b0' < 0) | (`d0' < 0) | (`h0' < 0) {
          di in r "Initial values must be > 0"
	  exit 198
	}	

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

	if "`weight'" == "" { qui ge `wi' = 1}
	else {qui ge `wi' `exp'}

	mark `touse' `if' `in'
	markout `touse' `varlist' 
	set more 1

	quietly {

	  count if `inc' < 0 & `touse'
	  local ct = _result(1)
	  if `ct' > 0 {
		noi di " "
		noi di in blue "Warning: `inc' has `ct' values < 0." _c
		noi di in blue " Not used in calculations"
		}

	  count if `inc' == 0 & `touse'
	  local ct = _result(1)
	  if `ct' > 0 {
		noi di " "
		noi di in blue "Warning: `inc' has `ct' values = 0." _c
		noi di in blue " Not used in calculations"
		}

	  /* Exclude obs with `inc' <= 0 from calculations  */

	  ge `badinc' = 0
	  replace `badinc' =. if `inc' <= 0
  	  markout `touse'  `badinc'

	}

	/*	Allow optional starting values 
		Cf. matrix c0 = (4, .1, 13)
	*/
	
	if `b0' <= 1 { local c0b = 1 }
	else { local c0b = ln(`b0') }

	if `d0' <= 1 { local c0d = 1 }
	else { local c0d = ln(`d0') }

	if `h0' <= 1 { local c0h = 13 }
	else { local c0h = ln(`h0') }

	/*
	if `h0' <= 1 {
		qui summ `inc' [w=`wi'] if `touse'
		local c0h = ln( (_result(3))^`d0' - 1 ) 
		}
	else {	local c0h = ln(`h0')  }
	*/

	matrix c0 = (`c0b', `c0d', `c0h')

*	matrix c0 = (4,.1,13)

	matrix colnames c0 = p1:_cons p2:_cons p3:_cons

	ml begin
	ml func dagum_ll
	ml method deriv0

	eq p1 :          
	eq p2 :          
	eq p3 :          

	ml model c = p1 p2 p3, depv(000) from(c0)

	global S_mldepn "`inc'"
	global S_mlwgt "`wi'"

	qui ml sample `mysamp' [w=`wi'] if `touse'

	if "`log'" != "" {
		local prel "quietly "
		}

	`prel' ml maximize `f' `V', `trace'  

	ml post mlnorm, title(ML fit of Dagum distribution)

	global S_E_cmd = "dagumfit"

  }  /* end of else */

	if `level' < 10 | `level' > 99 {
		local level = 95
		}

	ml mlout dagumfit, level(`level')

	matrix `c' = get(_b)
	matrix `v' = get(VCE)

	local b = exp(`c'[1,1]) 
	local d = exp(`c'[1,2]) 
	local h = exp(`c'[1,3]) + 1

	global S_b = "`b'"
	global S_d = "`d'"
	global S_h = "`h'"

	/* Derive s.e.s for b,d,h using delta method  [Greene, 3/e, p. 278]
	Given J functions of estimate vector f(b)
	JxK matrix C = df(b)/db'
	Est.Asym.Var[f(b)] = C[covar(b)]C'   JxK,KxK,KxJ
	... in the present case, C has only diagonal elements, so can treat
	as in univariate situation: cf Kendall & Stuart 5/e Vol. 1, p.324: 
		s.e.[g(x)] = (dg/dx)*s.e.[x]  
	*/

	local b_se = (`b') * sqrt(`v'[1,1])
	local d_se = (`d') * sqrt(`v'[2,2])
	local h_se = (`h'-1) * sqrt(`v'[3,3])
	
	di in gr "b = exp(p1)   = " _col(18) in ye %9.5f `b'  _c
	di in gr "; std. err. = " _col(15) in ye %9.5f `b_se' _c
	di in gr "; z = " in ye %9.5f `b'/`b_se'

	di in gr "d = exp(p2)   = " _col(18) in ye %9.5f `d'  _c
	di in gr "; std. err. = "  _col(15) in ye %9.5f `d_se' _c
	di in gr "; z = " in ye %9.5f `d'/`d_se'

	di in gr "h = 1+exp(p3) = " _col(18) in ye %9.5f `h'  _c
	di in gr "; std. err. = "  _col(15) in ye %9.5f `h_se' _c
	di in gr "; z = " in ye %9.5f `h'/`h_se'


/*	S.E. FORMULAE TO BE CHECKED
	* Derive s.e.s for a,b,q using delta method [Greene, 3/e, p. 278]

	local  z1 = `h'- 1

	matrix `del' = (`b',`d',`z1')
	matrix list `del'

	matrix `delv' = `del'*`v'
	matrix `delvdel' = `delv'*`del''
	matrix list `delvdel'

	local b_se = sqrt(`delvdel'[1,1])
	local d_se = sqrt(`delvdel'[2,2])
	local h_se = sqrt(`delvdel'[3,3])

	local b_t = `a'/`a_se'
	local d_t = `b'/`b_se'
	local h_t = `h'/`h_se'

	di " "
	di "s.e. estimates maybe not yet correct"
	di "b_se b_t d_se d_t h_se h_t   "
	di `b_se' " " `b_t' " " `d_se' " " `d_t' " " `h_se' " " `h_t'
*/

		/* Estimated Dagum c.d.f. */
	if "`cdf'" ~= "" {
		qui ge `cdf' = ( 1 + `h'*`inc'^(-`d') )^(-`b') if `touse'
	}

		/* Estimated Dagum p.d.f. */
	if "`pdf'" ~= "" {
		qui ge `pdf' = (`b'*`d'*`h')* (`inc')^(-`d'-1) /*
		*/ / (1 + `h'*`inc'^(-`d'))^(`b'+1) if `touse'
	}

	if "`stats'" ~= "" {

	/* (Option) stats implied by estimates */

	local mean = (`b'*`h'^(1/`d')) * exp(lngamma(1-1/`d')) /*
		*/ *exp(lngamma(`b'+1/`d'))/exp(lngamma(1+`b')) 
	local var =  (`b'*`h'^(2/`d')) * exp(lngamma(1-2/`d')) /*
		*/ *exp(lngamma(`b'+2/`d'))/exp(lngamma(1+`b')) /*
	        */ 	- (`mean'*`mean')
	local sd = sqrt(`var')
	local i2 = .5*`var'/(`mean'*`mean')
	local gini = -1 + ( (exp(lngamma(`b'))/exp(lngamma(2*`b')) ) /* 
	     */ / ( exp(lngamma(`b'+1/`d'))/exp(lngamma(2*`b'+1/`d')) ) )

	/* percentiles predicted from Dagum model */
	local x01 = `h'^(1/`d') * ( (.01)^(-1/`b') - 1 )^(-1/`d') 
	local x05 = `h'^(1/`d') * ( (.05)^(-1/`b') - 1 )^(-1/`d') 
	local x10 = `h'^(1/`d') * ( (.10)^(-1/`b') - 1 )^(-1/`d') 
	local x20 = `h'^(1/`d') * ( (.20)^(-1/`b') - 1 )^(-1/`d') 
	local x25 = `h'^(1/`d') * ( (.25)^(-1/`b') - 1 )^(-1/`d') 
	local x30 = `h'^(1/`d') * ( (.30)^(-1/`b') - 1 )^(-1/`d') 
	local x40 = `h'^(1/`d') * ( (.40)^(-1/`b') - 1 )^(-1/`d') 
	local x50 = `h'^(1/`d') * ( (.50)^(-1/`b') - 1 )^(-1/`d') 
	local x60 = `h'^(1/`d') * ( (.60)^(-1/`b') - 1 )^(-1/`d') 
	local x70 = `h'^(1/`d') * ( (.70)^(-1/`b') - 1 )^(-1/`d') 
	local x75 = `h'^(1/`d') * ( (.75)^(-1/`b') - 1 )^(-1/`d') 
	local x80 = `h'^(1/`d') * ( (.80)^(-1/`b') - 1 )^(-1/`d') 
	local x90 = `h'^(1/`d') * ( (.90)^(-1/`b') - 1 )^(-1/`d') 
	local x95 = `h'^(1/`d') * ( (.95)^(-1/`b') - 1 )^(-1/`d') 
	local x99 = `h'^(1/`d') * ( (.99)^(-1/`b') - 1 )^(-1/`d') 

	/* Lorenz curve ordinates at selected percentiles */
	local Lp01  = ibeta( `b'+1/`d',1-1/`d',(.01)^(1/`b') ) 
	local Lp05  = ibeta( `b'+1/`d',1-1/`d',(.05)^(1/`b') ) 
	local Lp10  = ibeta( `b'+1/`d',1-1/`d',(.10)^(1/`b') ) 
	local Lp20  = ibeta( `b'+1/`d',1-1/`d',(.20)^(1/`b') ) 
	local Lp25  = ibeta( `b'+1/`d',1-1/`d',(.25)^(1/`b') ) 
	local Lp30  = ibeta( `b'+1/`d',1-1/`d',(.30)^(1/`b') ) 
	local Lp40  = ibeta( `b'+1/`d',1-1/`d',(.40)^(1/`b') ) 
	local Lp50  = ibeta( `b'+1/`d',1-1/`d',(.50)^(1/`b') ) 
	local Lp60  = ibeta( `b'+1/`d',1-1/`d',(.60)^(1/`b') ) 
	local Lp70  = ibeta( `b'+1/`d',1-1/`d',(.70)^(1/`b') ) 
	local Lp75  = ibeta( `b'+1/`d',1-1/`d',(.75)^(1/`b') ) 
	local Lp80  = ibeta( `b'+1/`d',1-1/`d',(.80)^(1/`b') ) 
	local Lp90  = ibeta( `b'+1/`d',1-1/`d',(.90)^(1/`b') ) 
	local Lp95  = ibeta( `b'+1/`d',1-1/`d',(.95)^(1/`b') ) 
	local Lp99  = ibeta( `b'+1/`d',1-1/`d',(.99)^(1/`b') ) 

	di " "
	di in gr "Dagum model estimates for distribution of `inc'"
	di in gr _dup(60) "-"
	di in gr _col(5) "Percentiles" _col(17) "Cumulative shares of" 
	di in gr _col(17) "total `inc' (Lorenz ordinates)"
	di " "
	di in gr " 1%" _col(6) in ye %9.5f `x01' _col(20) %9.5f `Lp01'
	di in gr " 5%" _col(6) in ye %9.5f `x05' _col(20) %9.5f `Lp05'
	di in gr "10%" _col(6) in ye %9.5f `x10' _col(20) %9.5f `Lp10'
	di in gr "20%" _col(6) in ye %9.5f `x20' _col(20) %9.5f `Lp20'
	di in gr "25%" _col(6) in ye %9.5f `x25' _col(20) %9.5f `Lp25'
	di in gr "30%" _col(6) in ye %9.5f `x30' _col(20) %9.5f `Lp30'
	di in gr "40%" _col(6) in ye %9.5f `x40' _col(20) %9.5f `Lp40'
	di in gr "50%" _col(6) in ye %9.5f `x50' _col(20) %9.5f `Lp50' _c
	di in gr _col(10) "Mean" _col(22) in ye %9.5f `mean'
	di in gr "60%" _col(6) in ye %9.5f `x60' _col(20) %9.5f `Lp60' _c
	di in gr _col(10) "Std. Dev." _col(22) in ye %9.5f `sd'
	di in gr "70%" _col(6) in ye %9.5f `x70' _col(20) %9.5f `Lp70'
	di in gr "75%" _col(6) in ye %9.5f `x75' _col(20) %9.5f `Lp75' _c
	di in gr _col(10) "Variance" _col(22) in ye %9.5f `var'
	di in gr "80%" _col(6) in ye %9.5f `x80' _col(20) %9.5f `Lp80' _c
	di in gr _col(10) "Half CV^2" _col(22) in ye %9.5f `i2'
	di in gr "90%" _col(6) in ye %9.5f `x90' _col(20) %9.5f `Lp90' _c
	di in gr _col(10) "Gini coeff." _col(22) in ye %9.5f `gini'
	di in gr "95%" _col(6) in ye %9.5f `x95' _col(20) %9.5f `Lp95' _c
	di in gr _col(10) "p90/p10" _col(22) in ye %9.5f `x90'/`x10'
	di in gr "99%" _col(6) in ye %9.5f `x99' _col(20) %9.5f `Lp99' _c
	di in gr _col(10) "p75/p25" _col(22) in ye %9.5f `x75'/`x25'

	} /* end of Stats if block */


end