*! version 1.2.0 Stephen P. Jenkins April 2004 *! predict program to accompany -smfit- for predictions when covariates program define smpred, eclass version 8.2 if "`e(cmd)'" != "smfit" { di in red "smfit was not the last estimation command" exit 301 } syntax [if] [in] [, Aval(string) Bval(string) Qval(string) /// ABQval(string) CDF(namelist max=1) PDF(namelist max=1) /// POORfrac(real 0) ] if "`abqval'" ~= "" & (("`aval'"~="")|("`bval'"~="")|("`qval'"~="")) { di as error "Cannot use abqval(.) option in conjunction with aval(.), bval(.), qval(.) options" exit 198 } if "`abqval'" ~= "" { local aval "`abqval'" local bval "`abqval'" local qval "`abqval'" } if ("`if'" ~= "" | "`in'" ~= "") & ("`cdf'" == "") & ("`pdf'" == "") { di as txt "Warning: if and in have effect only if cdf or pdf options specified" } if "`cdf'" ~= "" { confirm new variable `cdf' } if "`pdf'" ~= "" { confirm new variable `pdf' } if "`poorfrac'" ~= "" { capture assert `poorfrac' < 0 if _rc == 0 { di as error "poorfrac value must be positive" exit 198 } } if `e(nocov)' ~= 1 & ("`aval'" == ""|"`bval'" == ""|"`qval'" == "") { di as error "Model used covariates: specify covariate values" exit } if `e(nocov)' == 1 { local a = `e(ba)' local b = `e(bb)' local q = `e(bq)' } else { tempname maval mbval mqval mat `maval' = (`aval') mat `mbval' = (`bval') mat `mqval' = (`qval') if colsof(`maval') ~= `e(length_b_a)' { di as error "# aval() elements ~= # covariates in a: eqn" exit } if colsof(`mbval') ~= `e(length_b_b)' { di as error "# bval() elements ~= # covariates in b: eqn" exit } if colsof(`mqval') ~= `e(length_b_q)' { di as error "# qval() elements ~= # covariates in q: eqn" exit } tempname ba bb bq mat `ba' = e(b_a)' mat `bb' = e(b_b)' mat `bq' = e(b_q)' mat `ba' = `maval'*`ba' mat `bb' = `mbval'*`bb' mat `bq' = `mqval'*`bq' local a = `ba'[1,1] local b = `bb'[1,1] local q = `bq'[1,1] } eret scalar bba = `a' eret scalar bbb = `b' eret scalar bbq = `q' eret scalar mean = `b'*exp(lngamma(1+1/`a')) /// *exp(lngamma(`q'-1/`a'))/exp(lngamma(`q')) eret scalar mode = cond(`a'>1,`b'*(((`a'-1)/(`a'*`q'+1))^(1/`a')),0,.) eret scalar var = `b'*`b'*exp(lngamma(1+2/`a')) /// *exp(lngamma(`q'-2/`a'))/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 = 1-(exp(lngamma(`q'))*exp(lngamma(2*`q'-1/`a')) /// / (exp(lngamma(`q'-1/`a'))*exp(lngamma(2*`q')))) local ptile "1 5 10 20 25 30 40 50 60 70 75 80 90 95 99" /* selected quantiles predicted from S-M model */ /* Lorenz curve ordinates at selected quantiles */ foreach x of local ptile { eret scalar p`x' = `b' * ( (1-`x'/100)^(-1/`q') - 1 )^(1/`a') eret scalar Lp`x' = ibeta(1+1/`a',`q' -1/`a',1-(1-`x'/100)^(1/`q')) } 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)' _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}" /* Fraction with income below given level */ if "`poorfrac'" ~= "" & `poorfrac' > 0 { eret scalar poorfrac = 1 - ( 1/(1 + (`poorfrac'/`b')^`a')^`q' ) eret scalar pline = `poorfrac' } if "`e(poorfrac)'" ~= "" { di " " di "Fraction with `e(depvar)' < `e(pline)' = " as res %9.5f e(poorfrac) di " " } marksample touse markout `touse' /* Estimated S-M c.d.f. */ if "`cdf'" ~= "" { qui ge `cdf' = 1 - ( 1/(1 + (`e(depvar)'/`b')^`a')^`q' ) if `touse' eret local cdfvar "`cdf'" } /* Estimated S-M p.d.f. */ if "`pdf'" ~= "" { qui ge `pdf' = (`a'*`q'/`b')* /// ((1 + (`e(depvar)'/`b')^`a')^-(`q'+1)) /// * ((`e(depvar)'/`b')^(`a'-1)) if `touse' eret local pdfvar "`pdf'" } end