*! version 1.0 15feb2007 Z. Sajaia program define fastgini, rclass syntax varname [if] [in] [pweight fweight], /// [bin(numlist integer min=1 max=1 >0) noCHeck jk Level(cilevel)] version 9.2 if missing("`check'") { tempvar goodx marksample touse quietly count if `varlist' <= 0 & `touse' local ct = r(N) if `ct' > 0 { display display as text "Warning: `varlist' has `ct' values <= 0. Not used in calculations" } quietly generate `goodx' = 1 if (`varlist' > 0) markout `touse' `goodx' `varlist' quietly count if `touse' local ct = r(N) if `ct'==0 error 2000 } if ~missing("`exp'") { tempvar w generate double `w' `exp' } local vlist "`varlist' `w'" tempname gini gini_jk se mse mata: mata_callfastgini() if missing("`jk'") display as text _n "Gini coefficient = " as result %9.7f `gini' else { display _newline display as text "Gini coefficient" _c display as text _col(56) "Number of obs" _col(70) "= " as res %7.0f `ct' display "" tempname Tab t p ll ul z .`Tab' = ._tab.new, col(7) lmargin(0) ignore(.b) // column 1 2 3 4 5 6 7 .`Tab'.width 13 |11 11 9 9 12 12 .`Tab'.titlefmt . %11s %12s %7s . %24s . .`Tab'.strfmt . %11s . . . . . .`Tab'.pad . 2 2 1 3 3 3 .`Tab'.numfmt . %9.0g %9.0g %8.2f %5.3f %9.0g %9.0g .`Tab'.strcolor . result . . . . . .`Tab'.sep, top .`Tab'.titles "`varlist'" /// 1 "Gini" /// 2 "Std. Err." /// 3 "t" /// 4 "P>|t|" /// 5 `"[`=strsubdp("`level'")'% Conf. Interval]"' "" // 6 7 scalar `t' = `gini'/`se' scalar `p' = 2*norm(-abs(`t')) scalar `z' = invnorm((100+`level')/200) scalar `ll' = `gini' - `se'*`z' scalar `ul' = `gini' + `se'*`z' .`Tab'.sep .`Tab'.row "" `gini' `se' `t' `p' `ll' `ul' .`Tab'.sep, bottom return scalar mse = `mse' return scalar gini_jk = `gini_jk' return scalar se = `se' } return scalar gini = `gini' end mata: void function mata_callfastgini() { RET = fastgini(st_data( .,tokens(st_local("vlist")), st_local("touse")), st_local("weight"), (st_local("jk")!=""), strtoreal(st_local("bin"))) st_numscalar(st_local("gini"), RET[1,1]) st_numscalar(st_local("se"), RET[2,1]) st_numscalar(st_local("gini_jk"), RET[3,1]) st_numscalar(st_local("mse"), RET[4,1]) } real matrix function fastgini(real matrix X, | string scalar weight, real scalar jk, real scalar bin) { colvector Xi, WX N = rows(X) if (missing(bin)) { // do 'regular' gini R = X[order(X,1),] M = N } else { // 'fast' gini, do aggregation M = bin MM = minmax(X[.,1]) bsize = (MM[1,2] - MM[1,1])/M R = ((1::M):*bsize:+MM[1, 1], J(M, 1, 0)) Xi = trunc((X[.,1]:-MM[1,1]):/bsize:-0.01):+1 if (cols(X)==1) { for (i=1; i<=N; ++i) { // no weights R[Xi[i], 2] = R[Xi[i], 2] + 1 } } else { for (i=1; i<=N; ++i) { R[Xi[i], 2] = R[Xi[i], 2] + X[i,2] } } } Xi=. C = cols(R) R = R \ J(1, C, .) if (C == 1) { sumW=N } else { WX =R[., 1]:*R[., 2] sumW=quadcolsum(R)[1,2] } sumWX=0 if (missing(jk)) { //----------------- JUST GINI --------------------------------- if (C == 1) { sumWX=quadcross(R, ((M::1):*2:-1) \ .)*0.5 sumX=quadcolsum(R) } else { sumX =0 for (i=1; i<=M; ++i) { sumWX = sumWX + R[i, 2]*(sumX+WX[i]*0.5) sumX = sumX + WX[i] } } g = 1- 2*sumWX/sumW/sumX RET = g } else { // ---------------------- GINI PLUS JACKKNIFE STANDARD ERRORS ------------------- SUM = J(M+1, 2, 0) // W x*W if (C == 1) { sumWX=quadcross(R, ((M::1):*2:-1) \ .)*0.5 SUM[., 1] = (0::M) for (i=1; i<=M; ++i) { SUM[i+1, 2] = SUM[i, 2] + R[i,1] } } else { for (i=1; i<=M; ++i) { sumWX = sumWX + R[i, 2]*(SUM[i, 2]+WX[i]*0.5) SUM[i+1, 1] = SUM[i, 1] + R[i, 2] SUM[i+1, 2] = SUM[i, 2] + WX[i] } } sumW=SUM[M+1, 1] sumX=SUM[M+1, 2] g = 1- 2*sumWX/sumW/sumX RET = g if (weight == "pweight") { G= 1:-2:*(sumWX:-(sumW:-R[.,2]:*0.5:-SUM[.,1]):*WX:-R[.,2]:*SUM[.,2]):/(sumW:-R[.,2]):/(sumX:-WX) MV=quadmeanvariance(G) V1=MV[2,1]+M/(M-1)*(MV[1,1]-g)^2 RET = RET \ sqrt(MV[2,1]*(M-1)^2/M) RET = RET \ MV[1,1] RET = RET \ sqrt(V1*(M-1)^2/M) } else { G= (- SUM[.,2]-R[.,1]:*(-SUM[.,1]:+(sumW-0.5)):+sumWX):/(R[.,1]:-SUM[M+1, 2]):*(2/(sumW - 1)):+1 if (C == 1) MV=quadmeanvariance(G) else MV=quadmeanvariance(G, R[.,2]) V1=MV[2,1]+sumW/(sumW-1)*(MV[1,1]-g)^2 RET = RET \ sqrt(MV[2,1]*(sumW-1)^2/sumW) RET = RET \ MV[1,1] RET = RET \ sqrt(V1*(sumW-1)^2/sumW) } } return(RET) } end