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

program define smfit
	version 4.0

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

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

	local varlist "req ex max(1)"
	local if "opt"
	local in "opt"
	local options "TRace noLOG Stats a0(real 2) b0(real 0.9) q0(real 2)"
	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 a_se b_se q_se 

	if (`a0' < 0) | (`b0' < 0) | (`q0' < 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"
		noi di in blue " Used in calculations"
		}

	  /* Exclude obs with `inc' <0 from calculations 
	   Substitute starred lines if want to exclude `inc'=0 as well */

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

	}

	/*	Allow optional starting values 
		Cf. matrix c0 = (1.0, 5.0, 0.5)
	*/
	if `a0' <= 1 { local c0a = 1 }
	else { local c0a = ln(`a0') }
	if `b0' <= 1 {
		qui summ `inc' [w=`wi'] if `touse'
		local c0b = ln( _result(3) )
		}
	else {	local c0b = ln(`b0')  }
	if `q0' <= 1 { local c0q = 1 }
	else { local c0q = ln(`q0') }
	
	matrix c0 = (`c0a', `c0b', `c0q')

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

	ml begin
	ml func sm_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 Singh-Maddala distribution)

	global S_E_cmd = "smfit"

  }  /* end of else */

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

	ml mlout smfit, level(`level')

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

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

	global S_a = "`a'"
	global S_b = "`b'"
	global S_q = "`q'"

	/* Derive s.e.s for a,b,q 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 a_se = (`a'-1) * sqrt(`v'[1,1])
	local b_se = (`b'-1) * sqrt(`v'[2,2])
	local q_se = `q' * sqrt(`v'[3,3])
	
	di in gr "a = 1+exp(p1) = " _col(18) in ye %9.5f `a'  _c
	di in gr "; std. err. = " _col(15) in ye %9.5f `a_se' _c
	di in gr "; z = " in ye %9.5f `a'/`a_se'

	di in gr "b = 1+exp(p2) = " _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 "q =   exp(p3) = " _col(18) in ye %9.5f `q'  _c
	di in gr "; std. err. = "  _col(15) in ye %9.5f `q_se' _c
	di in gr "; z = " in ye %9.5f `q'/`q_se'

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

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

	if "`stats'" ~= "" {

	/* (Option) stats implied by estimates */

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

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

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

	di " "
	di in gr "Singh-Maddala 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