*! version 1.2.2 Stephen P. Jenkins, July 2012 *! fix typo in post-estimation prediction of variance *! (Thanks to Michal Brzezinski for spotting the typo) *! version 1.2.1 Stephen P. Jenkins, March 2007 *! Fit GB2 distribution by ML to unit record data /*------------------------------------------------ playback request */ program define gb2fit, eclass byable(onecall) version 8.2 if replay() { if "`e(cmd)'" != "gb2fit" { noi di as error "results for gb2fit 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) CDF(namelist max=1) PDF(namelist max=1) POORfrac(real 0) /// Robust Cluster(varname) SVY STats From(string) /// Level(integer $S_level) /// 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' 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 txt "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 txt "Warning: {res:`inc'} has `ct' values = 0;" _c noi di as text " not used in calculations" } replace `touse' = 0 if `inc' <= 0 } qui count if `touse' if r(N) == 0 { error 2000 } if "`from'" != "" { local b0 "`from'" } global S_mlinc "`inc'" `log' ml model lf gb2fit_ll (a: `avar') (b: `bvar') (p: `pvar') (q: `qvar') /// `wgt' if `touse' , maximize /// collinear title(`title') `robust' `svy' init(`b0') /// search(on) `clopt' `level' `mlopts' `stdopts' `modopts' eret local cmd "gb2fit" eret local depvar "`inc'" tempname b ba bb bp bq mat `b' = e(b) mat `ba' = `b'[1,"a:"] mat `bb' = `b'[1,"b:"] mat `bp' = `b'[1,"p:"] mat `bq' = `b'[1,"q:"] eret matrix b_a = `ba' eret matrix b_b = `bb' eret matrix b_p = `bp' eret matrix b_q = `bq' eret scalar length_b_a = 1 + `na' eret scalar length_b_b = 1 + `nb' eret scalar length_b_p = 1 + `np' eret scalar length_b_q = 1 + `nq' if ("`avar'"!="" | "`bvar'"!="" | "`pvar'"!="" | "`qvar'"!="") { eret scalar nocov = 0 } if "`avar'"=="" & "`bvar'"=="" & "`pvar'"=="" & "`qvar'"=="" { tempname e mat `e' = e(b) local a = `e'[1,1] local b = `e'[1,2] local p = `e'[1,3] local q = `e'[1,4] eret scalar ba = `a' eret scalar bb = `b' eret scalar bp = `p' eret scalar bq = `q' eret scalar nocov = 1 /* 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'))*exp(lngamma(`q')) /// *( (1 +(`inc'/`b')^`a')^(`p'+`q') ) /// ) if `touse' eret local pdfvar "`pdf'" } /* selected quantiles predicted from GB2 model */ /* Lorenz curve ordinates at selected quantiles */ if "`stats'" ~= "" { eret scalar mean = `b'*exp(lngamma(`p'+1/`a')) /// *exp(lngamma(`q'-1/`a')) /// /( exp(lngamma(`p'))*exp(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')) /// *exp(lngamma(`q'-2/`a')) /// /( exp(lngamma(`p'))*exp(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 { local 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') ) } } /* 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' } } 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'" ml display, level(`level') `diopts' if `level' < 10 | `level' > 99 { local level = 95 } 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}" } end