/*------------------------------------*/
/*summclust */
/*written by Matt Webb */
/*version 4.210 2023-06-06 */
/*------------------------------------*/
version 13
cap program drop summclust
program define summclust, rclass
syntax varlist(min=3), cluster(varname) [fevar(varlist) absorb(varname) gstar Rho(real 0.1234567) SAMple(string) TABle ADDmeans JACKknife REGtable NOGraph]
local numvars = wordcount("`varlist'")
local yvar = word("`varlist'",1)
local CTEMP ""
local start = 1 + `numvars' + 1
/*split indepvars into variables of interest and controls*/
local XTEMP = word("`varlist'",2)
forvalues i = 3/`numvars' {
local wordi = word("`varlist'",`i')
local CTEMP = "`CTEMP' `wordi' "
}
local FEVAR "`fevar'"
local SMPL "`sample'"
local svars = "`addmeans'"
/*add if to sample*/
if "`sample'" != "" {
local SMPL = "if `sample'"
}
/*check rho argument*/
if inrange(`rho',0.00001,1)==0 {
disp as error "The value for rho, `rho', is invalid. rho must be between 0 and 1."
exit 198
}
else {
/*turns gstar on if rho option is specified*/
if inrange(`rho',0.1234,0.1235) == 0 {
local gstarind = 1
}
else {
local gstarind = 0
}
}
/*flag to calculate g-star*/
if "`gstar'" != "" {
local gstarind = 1
}
/*flag for absorb used properly*/
local nocv3 = 0
/*flags for singular clusters*/
local singclust = 0
local singclustfull = 0
local colsdrop = 0
local countfevar = 0
/*declare tempvars*/
tempvar ytilde
tempvar temp_sd
local fevars " "
/*xvar labels for regtable*/
local name_xvar = "`CTEMP'"
local catvars = " "
foreach catvar in `FEVAR' {
local catvars = " `catvars' " + "i.`catvar'"
}
local absvars = " "
foreach avar in `absorb' {
local absvars = " `absvars' " + "i.`avar'"
}
tempvar temp_sample
mark `temp_sample'
markout `temp_sample' `yvar' `XTEMP' `CTEMP' `catvars' `absvars' `cluster',strok
/*impose the "sample option" restriction*/
if "`sample'" != "" {
tempvar temp_samp
qui gen `temp_samp' = 0
qui replace `temp_samp' = 1 `SMPL'
qui replace `temp_sample' = `temp_sample'*`temp_samp'
}
/*for fevar*/
if `"`fevar'"' != "" {
*regtable names labels
local names_fe = ""
local tempvars = " "
local j=0
local q = 0
foreach fvar in `FEVAR' {
local q = `q' + 1
local countfevar = `countfevar' + 1
qui levelsof `fvar' if `temp_sample'==1, local(fevarlevel)
local jstart = `j'+1
foreach flevel in `fevarlevel' {
local j = `j'+1
local tj = "t_`j'"
tempvar `tj'
qui gen ``tj'' = `fvar'==`flevel'
/*create labels for these variables for regtable*/
local name_tj = "`fvar'_`flevel'"
local names_fe = "`names_fe' `name_tj'"
} /*end of flevel*/
/*drop the last fe for the second and greater fe*/
local jend = `j'-1
if `q'==1 & `"`absorb'"' == "" {
local fevars = "`' `t_1'-`t_`j' ' "
local names_fe = subinstr("`names_fe'", "`name_tj'"," ", .)
}
else if `q'==1 & `"`absorb'"' != "" {
local fevars = "`' `t_1'-`t_`jend' ' "
local names_fe = subinstr("`names_fe'", "`name_tj'"," ", .)
}
else if `q' > 1 {
local fevars = "`fevars' `t_`jstart'' - `t_`jend'' "
local names_fe = subinstr("`names_fe'", "`name_tj'"," ", .)
}
} /*end of fvar*/
*local fevars = "`' `t_1'-`t_`j' ' "
local XVAR = "`CTEMP'"
local CTEMP = "`CTEMP' `fevars'"
} /*end of fevar*/
/*sort the data by the clustering variable*/
/*make a sequential-numeric cluster indicator*/
tempvar temp_indexa
qui egen `temp_indexa' = group(`cluster') if `temp_sample'==1
qui summ `temp_indexa'
local G = r(max)
qui sort `temp_indexa'
qui putmata CVAR=`temp_indexa' if `temp_sample'==1, replace
mata: numobs = rows(CVAR)
mata: info = panelsetup(CVAR,1)
mata: ng = info[.,2] - info[.,1] :+ 1
mata: xtilde= J(numobs,0,.)
local xtilde = " "
local XVAR = "`CTEMP'"
/*dump everything into mata here, then demean by absorb*/
qui gen `ytilde' = `yvar' if `temp_sample'==1
foreach xvar in `XTEMP' {
qui putmata `xvar' = `xvar' if `temp_sample'==1, replace
mata: xtilde = xtilde, `xvar'
local xtilde = "`xtilde'" + " ``varname''"
}
if `"`absorb'"' != "" {
/* want to check that cluster is constant within absorb */
qui putmata absorb =`absorb' if `temp_sample'==1, replace
qui putmata cabsorb = `temp_indexa' if `temp_sample'==1, replace
mata: ainfo = panelsetup(absorb,1)
mata: ang = ainfo[.,2] - ainfo[.,1] :+ 1
mata: absave = panelsum(cabsorb,ainfo) :/ ang
mata: absorb2 = cabsorb:*cabsorb
mata: absorb2ave = panelsum(absorb2, ainfo) :/ ang
mata absconstant = floatround(absorb2ave)!=floatround( absave:*absave)
mata: st_numscalar("nocv3", absconstant)
mata: st_local("nocv3", strofreal(absconstant))
/*other way*/
if "`absorb'" == "`cluster'" {
local nogood = 0
}
else {
/*check for invariance of cluster variable*/
cap drop `temp_sd'
qui bysort `absorb': egen `temp_sd' = sd(`cluster')
qui summ `temp_sd'
local sd_mean = r(mean)
if inrange(`sd_mean',-0.0001,0.0001)==1{
local nogood = 0
}
else{
local nogood = 1
disp as error "Cluster variable not constant within absorb variable. Use fevar instead."
exit 198
}
}
}
/*maybe make this for when there is neither absorb and fevar*/
if `"`absorb'"' == "" & `"`fevar'"' == "" {
/*add constant*/
tempvar tempcons
qui gen `tempcons' = 1
local xtilde = "`xtilde'" + " `tempcons' "
mata: ones = J(rows(xtilde),1,1)
mata: xtilde = xtilde, ones
}
/*put everything in mata*/
sort `temp_indexa'
qui putmata Y = `ytilde' if `temp_sample'==1, replace
qui putmata XOTHER = (`CTEMP' ) if `temp_sample'==1, replace
mata: X = xtilde, XOTHER
mata: k = cols(X)
mata: st_numscalar("k", k)
mata: st_local("k", strofreal(k))
mata: numobs = rows(X)
mata: st_numscalar("numobs", numobs)
mata: st_local("numobs", strofreal(numobs))
mata: sk = rank(X)
mata: st_numscalar("sk", sk)
mata: st_local("sk", strofreal(sk))
local coeffnum = `sk'
mata: iota =J(numobs,1,1)
/*for absorb option*/
if `"`absorb'"' != "" {
/*need to partition everything by absorb variable*/
mata: abinfo = panelsetup(absorb,1)
mata : abng = abinfo[.,2] - abinfo[.,1] :+ 1
mata: numab = rows(abinfo)
mata: st_numscalar("H", numab)
mata: st_local("H", strofreal(numab))
local coeffnum = `sk' + `H'
/*demean y */
mata: ymean = panelsum(Y,abinfo) :/ abng
mata: Ydm = J(0,1,0)
forvalues g=1/`H' {
mata: Ydm`g' = panelsubmatrix(Y,`g',abinfo)
mata: Ydm`g' = Ydm`g' :- ymean[`g',1]
mata: Ydm = Ydm \ Ydm`g'
}
mata Y = Ydm
/*demean x*/
mata: xmean = panelsum(X, abinfo) :/ abng
mata: Xdm = J(0,cols(X),0)
forvalues g=1/`H' {
mata: Xdm`g' = panelsubmatrix(X,`g',abinfo)
mata: Xdm`g' = Xdm`g' :- xmean[`g',.]
mata: Xdm = Xdm \ Xdm`g'
}
mata X = Xdm
} /*end of absorb if*/
/*count number of x variables*/
mata: xnum = cols(xtilde)
/*create all the submatrices*/
forvalues g=1/`G' {
mata: X`g' = panelsubmatrix(X,`g',info)
mata: Y`g' = panelsubmatrix(Y,`g',info)
mata: iota`g' = panelsubmatrix(iota,`g',info)
mata: X`g'pX`g' = quadcross(X`g',X`g')
mata: X`g'pY`g' = quadcross(X`g',Y`g')
}
/*summ over the sub-blocks*/
mata: XpX = X1pX1
mata: XpY = X1pY1
forvalues g=2/`G' {
mata: XpX = XpX + X`g'pX`g'
mata: XpY = XpY + X`g'pY`g'
}
mata: XpXi = invsym(XpX)
mata: beta = XpXi*XpY
/*generate scores*/
forvalues g=1/`G' {
mata: W`g' = Y`g' - X`g' *beta
mata: Z`g' = X`g''* W`g'
mata: Z`g'Z`g'p = Z`g'*Z`g''
}
mata: ZZp = Z1Z1p
forvalues g=2/`G'{
mata: ZZp = ZZp + Z`g'Z`g'p
}
mata: crve1 = XpXi * ZZp * XpXi
local gfact = (`G'*(`numobs'-1))/((`G'-1)*(`numobs'-`coeffnum'))
mata: crve1 = `gfact':* crve1
mata: betatemp = beta[1,1]
mata: st_numscalar("BETATEMP", betatemp)
mata: st_local("BETATEMP", strofreal(betatemp))
mata: cv1se = sqrt(crve1[1,1])
mata: st_numscalar("setemp", cv1se)
mata: st_local("setemp", strofreal(cv1se))
local cv1t = `BETATEMP'/`setemp'
local cv1p = 2*min(ttail(`G'-1, `cv1t'), 1-ttail(`G'-1, `cv1t'))
local cv1lci = `BETATEMP' - invttail(`G'-1,0.025) * `setemp'
local cv1uci = `BETATEMP' + invttail(`G'-1,0.025) * `setemp'
/*beta-g, Hg, Lg*/
mata: betags = J(k,0,.)
mata: betadropgs = J(k,0,.)
mata: Leverage = J(`G',1,.)
mata: cv3block = J(k,`G',.)
/*string for problem clusters*/
local probclust = ""
local probclustfull = ""
forvalues g=1/`G'{
mata: left = XpX - X`g'pX`g'
mata: right = XpY - X`g'pY`g'
mata: leftinv = invsym(left)
mata: beta`g' = leftinv *right
mata: betags = betags, beta`g'
mata: H2`g' = X`g''*X`g'*XpXi
mata: L`g' = trace(H2`g')
mata: Leverage[`g',1] = L`g'
/*cv3*/
mata: betatemp = (beta :- beta`g')
mata: cv3block[.,`g'] = betatemp
/*check for singularities for all elements*/
mata anyzerosmall = beta`g':<= 0.00000001
mata anyzerobig = beta`g':>= -0.00000001
mata anyzero = anyzerosmall :* anyzerobig
mata countzero = colsum(anyzero)
mata: st_numscalar("countzero", countzero)
mata: st_local("countzero", strofreal(countzero))
if `countzero' > 0 {
mata: betadrop`g' = J(rows(beta`g'),1,0)
local singclustfull = `singclustfull' + 1
local probclustfull = "`probclustfull'" + " `g'"
}
else {
mata: betadrop`g' = beta`g'
}
mata: betadropgs = betadropgs, betadrop`g'
} /*end of g loop*/
local singfrac = (`singclustfull'/ `G')*100
if `singclustfull' >=1 {
mata: betadropgs = select(betadropgs, betadropgs[1.,]:!=0)
mata gdrop = cols(betadropgs)
mata: st_numscalar("gdrop", gdrop)
mata: st_local("gdrop", strofreal(gdrop))
}
/*jackknife*/
if "`jackknife'" != "" {
mata: betagbar = mean(betags')
mata: betadiff = betags' :- betagbar
mata: betadiffsum = betadiff[1,.]'* betadiff[1,.]
forvalues g = 2 /`G'{
mata: betadiffsum = betadiffsum + betadiff[`g',.]'* betadiff[`g',.]
}
mata: cv3j = sqrt( ((`G'-1)/`G') :* betadiffsum)
mata: cv3jse = cv3j[1,1]
mata: st_numscalar("cv3jse", cv3jse)
mata: st_local("cv3jse", strofreal(cv3jse))
local cv3jt = `BETATEMP'/`cv3jse'
local cv3jp = 2*min(ttail(`G'-1, `cv3jt'), 1-ttail(`G'-1, `cv3jt'))
local cv3jlci = `BETATEMP' - invttail(`G'-1,0.025) * `cv3jse'
local cv3juci = `BETATEMP' + invttail(`G'-1,0.025) * `cv3jse'
/*dropped singular version*/
if `singclustfull' >= 1 {
mata: betagdropbar = mean(betadropgs')
mata: betadiffdrop = betadropgs' :- betagdropbar
mata: betadiffsumdrop = betadiffdrop[1,.]'* betadiffdrop[1,.]
forvalues g = 2 /`gdrop'{
mata: betadiffsumdrop = betadiffsumdrop + betadiffdrop[`g',.]'* betadiffdrop[`g',.]
}
mata: cv3jdrop = sqrt( ((`gdrop'-1)/`gdrop') :* betadiffsumdrop)
mata: cv3jsedrop = cv3jdrop[1,1]
mata: st_numscalar("cv3jsedrop", cv3jsedrop)
mata: st_local("cv3jsedrop", strofreal(cv3jsedrop))
local cv3jtdrop = `BETATEMP'/`cv3jsedrop'
local cv3jpdrop = 2*min(ttail(`gdrop'-1, `cv3jtdrop'), 1-ttail(`gdrop'-1, `cv3jtdrop'))
local cv3jlcidrop = `BETATEMP' - invttail(`gdrop'-1,0.025) * `cv3jsedrop'
local cv3jucidrop = `BETATEMP' + invttail(`gdrop'-1,0.025) * `cv3jsedrop'
} /*end of singclustfull*/
} /*end of jackknife*/
/*partial leverage*/
qui reg ``varname'' ``tempcons'' `CTEMP' if `temp_sample'==1, noc
local xtt = "xtt"
tempvar xtt
qui predict `xtt' if `temp_sample'==1, resid
qui putmata xtt = `xtt' if `temp_sample'==1, replace
mata: xttxtt = cross(xtt,xtt)
mata: pLeverage = J(`G',1,.)
mata: cv3sum = J(k,k,0)
mata: cv3sumdrop = J(k,k,0)
forvalues g=1/`G'{
mata: xtt`g' = panelsubmatrix(xtt,`g',info)
mata: pL`g' = cross(xtt`g',xtt`g') / xttxtt
mata: pLeverage[`g',1] = pL`g'
mata: cv3sum = cv3sum + cv3block[.,`g']*cv3block[.,`g']'
}
mata: cv3 = ((`G'-1)/`G'):*cv3sum
mata: cv3se = sqrt(cv3[1,1])
mata: st_numscalar("cv3se", cv3se)
mata: st_local("cv3se", strofreal(cv3se))
local cv3t = `BETATEMP'/`cv3se'
local cv3p = 2*min(ttail(`G'-1, `cv3t'), 1-ttail(`G'-1, `cv3t'))
local cv3lci = `BETATEMP' - invttail(`G'-1,0.025) * `cv3se'
local cv3uci = `BETATEMP' + invttail(`G'-1,0.025) * `cv3se'
local SUMVAR "ng Leverage pLeverage betag"
/*matrices to store output*/
mata: clustsum= J(7,4,.)
mata: bonus = J(6,4,.)
/*dropped singular version*/
if `singclustfull' >=1 {
forvalues g = 1/`gdrop' {
mata: betatempdrop = (beta :- betadropgs[.,`g'])
mata: cv3sumdrop = cv3sumdrop + betatempdrop*betatempdrop'
}
mata: cv3drop = ((`gdrop'-1)/`gdrop'):*cv3sumdrop
mata: cv3sedrop = sqrt(cv3drop[1,1])
mata: st_numscalar("cv3sedrop", cv3sedrop)
mata: st_local("cv3sedrop", strofreal(cv3sedrop))
local cv3tdrop = `BETATEMP'/`cv3sedrop'
local cv3pdrop = 2*min(ttail(`gdrop'-1, `cv3tdrop'), 1-ttail(`gdrop'-1, `cv3tdrop'))
local cv3lcidrop = `BETATEMP' - invttail(`gdrop'-1,0.025) * `cv3sedrop'
local cv3ucidrop = `BETATEMP' + invttail(`gdrop'-1,0.025) * `cv3sedrop'
/*statistics adding dropped betas*/
mata: colsdrop = cols(betadropgs)
mata: st_numscalar("colsdrop", colsdrop)
mata: st_local("colsdrop", strofreal(colsdrop))
if `colsdrop' > 0 {
local SUMVAR "ng Leverage pLeverage betag betadropg"
mata: clustsum= J(7,5,.)
mata: bonus = J(6,5,.)
mata: betadropg = betadropgs[1,.]'
}
} /*end of singclustfull */
mata: betag = betags[1,.]'
local s = 0
foreach svar in `SUMVAR' {
local s = `s' + 1
local tempsvar = "temp_`svar'"
tempvar `tempsvar'
qui getmata ``tempsvar'' = `svar', force
qui summ ``tempsvar'', det
/*min*/
local min = r(min)
mata: clustsum[1,`s'] = `min'
/*q1*/
local q1 = r(p25)
mata: clustsum[2,`s'] = `q1'
/*median*/
local median = r(p50)
mata: clustsum[3,`s'] = `median'
/*mean*/
local mean = r(mean)
mata: clustsum[4,`s'] = `mean'
local graph_`s' = `mean'
/*q3*/
local q1 = r(p75)
mata: clustsum[5,`s'] = `q1'
/*max*/
local max = r(max)
mata: clustsum[6,`s'] = `max'
/*coeff of variation*/
mata: meandiff = `svar' :- `mean'
mata: meandiff2 = meandiff:*meandiff
mata: meandiffsum = colsum(meandiff2)
mata: denom= 1/((`G'-1)*(`mean'^2))
mata: scalvar = denom*meandiffsum
mata: scalvar = sqrt(abs(scalvar))
mata: clustsum[7,`s'] = scalvar
/*svar options*/
if "`svars'"!="" {
/*harmonic mean*/
mata: sinv = 1:/`svar'
mata: meaninv = mean(sinv)
mata: harmonic = 1/meaninv
mata: bonus[1,`s'] = harmonic
mata: harmrat = harmonic/`mean'
mata: bonus[2,`s'] = harmrat
/*geometric mean*/
mata: logvar = log(`svar')
mata: meanlog = mean(logvar)
mata: geo = exp(meanlog)
mata: bonus[3,`s'] = geo
mata: georat = geo/`mean'
mata: bonus[4,`s'] = georat
/*quadratic mean*/
mata: squares = `svar':*`svar'
mata: sqmean = mean(squares)
mata: qdmean = sqrt(sqmean)
mata: bonus[5,`s'] = qdmean
mata: quadrat = qdmean/`mean'
mata: bonus[6,`s'] = quadrat
} /*end of svars if*/
mata: st_matrix("`svar'", `svar')
} /*end of sum var loop*/
/*replace harmonic and geometric values for betanog*/
mata: bonus[1,4] = .
mata: bonus[2,4] = .
mata: bonus[3,4] = .
mata: bonus[4,4] = .
mata: st_matrix("clustsum", clustsum)
matrix rownames clustsum = min q1 median mean q3 max coefvar
if `singclustfull' >=1 & `colsdrop' >= 1 {
matrix colnames clustsum = Ng Leverage "Partial L." "all beta no g" "kept beta no g"
mata: bonus[1,5] = .
mata: bonus[2,5] = .
mata: bonus[3,5] = .
mata: bonus[4,5] = .
}
else {
matrix colnames clustsum = Ng Leverage "Partial L." "beta no g"
}
/*rowlabels for cluster by cluster matrix*/
qui levelsof `cluster' if `temp_sample'==1, local(CLUSTERNAMES)
mata: cnames = J(`G',1,".")
local f = 0
foreach cnames in `CLUSTERNAMES'{
local f = `f' + 1
mata: cnames[`f',1] = "`cnames'"
local clustlabel`f' = "`cnames'"
}
/****************************/
/*gstar option*/
/***************************/
/*check if fixed effects invariant within clusters*/
if `gstarind' == 1 {
if `"`fevar'"' != "" {
foreach fvar in `FEVAR' {
if "`fvar'" == "`cluster'" {
local gstarind = 2
}
else {
cap drop `temp_sd'
qui bysort `fvar': egen `temp_sd' = sd(`cluster')
qui summ `temp_sd'
local sd_mean = r(mean)
if inrange(`sd_mean',-0.0001,0.0001)==1{
local gstarind = 2
}
}
} /*end of fvar*/
} /*end of fevar if*/
if `"`absorb'"' != "" {
if "`absorb'" == "`cluster'" {
local gstarind = 2
}
else {
/*check for invariance of cluster variable*/
cap drop `temp_sd'
qui bysort `absorb': egen `temp_sd' = sd(`cluster')
qui summ `temp_sd'
local sd_mean = r(mean)
if inrange(`sd_mean',-0.0001,0.0001)==1{
local gstarind = 2
}
}
} /*end of absorb if*/
} /*end of if 1*/
if `gstarind' == 1 {
forvalues g= 1/ `G'{
mata: w`g' = XpXi[.,1]
mata: gamzero`g' = w`g'' * X`g'pX`g' * w`g'
mata: gamone`g' = (iota`g'' * X`g' * w`g')' * (iota`g'' * X`g' * w`g')
mata: gamrho`g' = `rho' * gamone`g' + (1-`rho')*gamzero`g'
}
mata: gamzt = 0
mata: gamrt = 0
mata: gamot = 0
forvalues g = 1 / `G'{
mata: gamzt = gamzt + gamzero`g'
mata: gamot = gamot + gamone`g'
mata: gamrt = gamrt + gamrho`g'
}
mata: gambarz = gamzt/ `G'
mata: gambaro = gamot / `G'
mata: gambarr = gamrt / `G'
mata: gammazero = 0
mata: gammaone = 0
mata: gammarho = 0
forvalues g = 1 / `G'{
mata: tempone = ((gamone`g' - gambaro)/gambaro)^2
mata: tempzero = ((gamzero`g' - gambarz)/gambarz)^2
mata: temprho = ((gamrho`g'- gambarr)/gambarr)^2
mata: gammazero = gammazero + tempzero
mata: gammaone = gammaone + tempone
mata: gammarho = gammarho + temprho
}
mata: gammazero = gammazero / `G'
mata: gammaone = gammaone / `G'
mata: gammarho = gammarho / `G'
mata: gstarone = `G'/(1 + gammaone)
mata: gstarzero = `G'/(1 + gammazero)
mata: gstarrho = `G'/(1 + gammarho)
mata: st_numscalar("gstarzero", gstarzero)
mata: st_local("gstarzero", strofreal(gstarzero))
mata: st_numscalar("gstarone", gstarone)
mata: st_local("gstarone", strofreal(gstarone))
mata: st_numscalar("gstarrho", gstarrho)
mata: st_local("gstarrho", strofreal(gstarrho))
}
else if `gstarind' == 2 {
forvalues g= 1/ `G'{
mata: w`g' = XpXi[.,1]
mata: gamzero`g' = w`g'' * X`g'pX`g' * w`g'
}
mata: gamzt = 0
forvalues g = 1 / `G'{
mata: gamzt = gamzt + gamzero`g'
}
mata: gambarz = gamzt/ `G'
mata: gammazero = 0
forvalues g = 1 / `G'{
mata: tempzero = ((gamzero`g' - gambarz)/gambarz)^2
mata: gammazero = gammazero + tempzero
}
mata: gammazero = gammazero / `G'
mata: gstarzero = `G'/(1 + gammazero)
mata: st_numscalar("gstarzero", gstarzero)
mata: st_local("gstarzero", strofreal(gstarzero))
}
/*output display*/
disp ""
disp ""
disp "SUMMCLUST - MacKinnon, Nielsen, and Webb"
disp " "
disp "Cluster summary statistics for `XTEMP' when clustered by `cluster'."
disp "There are `numobs' observations within `G' `cluster' clusters."
/*singular clusters warnings*/
if `singclustfull' >= 1 & `colsdrop' >= 1{
mata problems = J(0,1,".")
local clustnames = ""
foreach problem in `probclustfull' {
mata: problems = problems \ cnames[`problem',1]
local clustnames = "`clustnames' " + "`clustlabel`problem'' "
}
local singfrac = round(`singfrac',0.01)
disp " "
disp as error "***************************************************************************"
disp "WARNING - Elements of beta undefined when certain cluster(s) are omitted"
disp "***************************************************************************"
disp "Two standard errors are calculated: "
disp "The first standard error uses a generalized inverse."
disp "The second standard error drops the singularities."
disp "***************************************************************************"
disp "There are `singclustfull' problem clusters, out of `G' clusters."
disp "The problematic `cluster' cluster(s) are: `clustnames'"
if `singfrac' < 20 {
disp "As only" %7.2f `singfrac' " % of subsamples are singular, dropping them is preferred."
}
else {
disp "As" %7.2f `singfrac' "% of subsamples are singular, the generalized inverse is preferred."
}
disp " "
disp " "
} /*end of singclust if*/
else if `singclustfull' >= 1 & `colsdrop' == 0{
mata problems = J(0,1,".")
local clustnames = ""
foreach problem in `probclustfull' {
mata: problems = problems \ cnames[`problem',1]
local clustnames = "`clustnames' " + "`clustlabel`problem'' "
}
local singfrac = round(`singfrac',0.01)
disp " "
disp as error "***************************************************************************"
disp "WARNING - Elements of beta undefined when certain cluster(s) are omitted"
disp "***************************************************************************"
disp "The standard error uses a generalized inverse."
disp "***************************************************************************"
disp "There are `singclustfull' problem clusters, out of `G' clusters."
disp "All clusters are problematic."
disp " "
disp " "
} /*end of singclust else if*/
if "`jackknife'"!="" {
matrix define cvstuff = J(3,6,.)
matrix rownames cvstuff = CV1 CV3 CV3J
local RSPECCV "&|&&|"
}
else {
matrix define cvstuff = J(2,6,.)
matrix rownames cvstuff = CV1 CV3
local RSPECCV "&|&|"
}
/*check if absorb used improperly*/
if `nocv3'==1 {
matrix cvstuff[1,1] = `BETATEMP'
matrix cvstuff[1,2] = `setemp'
matrix cvstuff[1,3] = `cv1t'
matrix cvstuff[1,4] = `cv1p'
matrix cvstuff[1,5] = `cv1lci'
matrix cvstuff[1,6] = `cv1uci'
matrix cvstuff[2,1] = `BETATEMP'
disp" "
disp as error "*********************************************************************"
disp "WARNING - `cluster' is not constant within the absorb variable, `absorb' ".
disp "Beta no g, and CV3 are not computed correctly in this case. "
disp "Use the fevar option for `absorb' instead, if the above are desired."
disp "*********************************************************************"
disp" "
/*also replace the beta-no-g, partial leverage, and leverage*/
forvalues i = 1/ 7 {
matrix clustsum[`i',2] = .
matrix clustsum[`i',3] = .
matrix clustsum[`i',4] = .
}
} /*end of nocv3 */
else {
matrix cvstuff[1,1] = `BETATEMP'
matrix cvstuff[1,2] = `setemp'
matrix cvstuff[1,3] = `cv1t'
matrix cvstuff[1,4] = `cv1p'
matrix cvstuff[1,5] = `cv1lci'
matrix cvstuff[1,6] = `cv1uci'
matrix cvstuff[2,1] = `BETATEMP'
matrix cvstuff[2,2] = `cv3se'
matrix cvstuff[2,3] = `cv3t'
matrix cvstuff[2,4] = `cv3p'
matrix cvstuff[2,5] = `cv3lci'
matrix cvstuff[2,6] = `cv3uci'
if "`jackknife'"!="" {
matrix cvstuff[3,1] = `BETATEMP'
matrix cvstuff[3,2] = `cv3jse'
matrix cvstuff[3,3] = `cv3jt'
matrix cvstuff[3,4] = `cv3jp'
matrix cvstuff[3,5] = `cv3jlci'
matrix cvstuff[3,6] = `cv3juci'
if `nocv3'== 1 {
forvalues j= 2 /6 {
matrix cvstuff[3,`j'] = .
}
}
}
} /*end else*/
/*singular clusters*/
if `singclustfull' >=1 {
if "`jackknife'" != "" {
matrix define cvstuffdrop = J(2,6,.)
matrix rownames cvstuffdrop = CV3 CV3J
local RSPECCVDROP "&|&|"
matrix cvstuffdrop[2,1] = `BETATEMP'
matrix cvstuffdrop[2,2] = `cv3jsedrop'
matrix cvstuffdrop[2,3] = `cv3jtdrop'
matrix cvstuffdrop[2,4] = `cv3jpdrop'
matrix cvstuffdrop[2,5] = `cv3jlcidrop'
matrix cvstuffdrop[2,6] = `cv3jucidrop'
}
else {
matrix define cvstuffdrop = J(1,6,.)
matrix rownames cvstuffdrop = CV3
local RSPECCVDROP "|&|"
}
matrix cvstuffdrop[1,1] = `BETATEMP'
matrix cvstuffdrop[1,2] = `cv3sedrop'
matrix cvstuffdrop[1,3] = `cv3tdrop'
matrix cvstuffdrop[1,4] = `cv3pdrop'
matrix cvstuffdrop[1,5] = `cv3lcidrop'
matrix cvstuffdrop[1,6] = `cv3ucidrop'
} /* end of singclustfull if */
matrix colnames cvstuff = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper
matlist cvstuff, title(Regression Output) rowtitle(s.e.) ///
cspec(& %-4s w6 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 & %9.6f w10 & %9.6f w10 &) ///
rspec("`RSPECCV'")
if `nocv3'== 1 {
disp as error "WARNING - CV3 estimates not available as absorb option improperly specified."
}
else {
if `singclustfull' >=1 & `colsdrop' >= 1 {
matrix colnames cvstuffdrop = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper
matlist cvstuffdrop, title(Regression Output -- Dropping Singular Omit-One-Cluster Subsamples) rowtitle(s.e.) ///
cspec(& %-4s w6 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 & %9.6f w10 & %9.6f w10 &) ///
rspec("`RSPECCVDROP'")
}
}
mata: cvstuff=st_matrix("cvstuff")
if `singclustfull' >=1 & `colsdrop' >= 1 {
matlist clustsum , title("Cluster Variability") rowtitle(Statistic) ///
cspec(& %-10s | %8.2f o2 & %9.6f o2 & %9.6f w10 & %9.6f o2 & %9.6f o2 &) ///
rspec(&-&&&&&-&)
}
else {
matlist clustsum , title("Cluster Variability") rowtitle(Statistic) ///
cspec(& %-10s | %8.2f o4 & %9.6f o4 & %9.6f w10 & %9.6f o4 &) ///
rspec(&-&&&&&-&)
}
/*error for no-beta-g*/
if `nocv3'==1 {
disp as error "Some results not available as absorb option is improperly specified."
disp as text ""
}
/*-------------------*/
/*optional output*/
/*-------------------*/
/*gstar*/
if `gstarind'==1 {
disp " "
disp "Effective Number of Clusters"
disp "-----------------------------"
disp "G*(0) = " %6.3f `gstarzero'
if inrange(`rho',0.1234,0.1235)==0 {
disp "G*(`rho') = " %6.3f `gstarrho'
}
disp "G*(1) = " %6.3f `gstarone'
disp "-----------------------------"
}
else if `gstarind'==2 {
disp " "
disp "Effective Number of Clusters"
disp "-----------------------------"
disp "G*(0) = " %6.3f `gstarzero'
disp "-----------------------------"
disp as error "G*(rho) and G*(1) are not available."
disp as error "There are fixed effects at the cluster or subcluster level."
}
/*additional summary statistics*/
if "`svars'"!=""{
mata: st_matrix("bonus", bonus)
matrix rownames bonus = "Harmonic Mean" "Harmonic Ratio" "Geometric Mean" "Geometric Ratio" "Quadratic Mean" "Quadratic Ratio"
if `singclustfull' >= 1 & `colsdrop' >= 1 {
matrix colnames bonus = Ng Leverage "Partial L." "all beta no g" "kept beta no g"
matlist bonus, title("Alternative Sample Means and Ratios to Arithmetic Mean") ///
cspec(& %-10s w15 | %11.3f o2 & %9.6f w10 & %9.6f w10 & %9.6f o2 & %9.6f o2 &) rspec(&|&&&&&|)
}
else {
matrix colnames bonus = Ng Leverage "Partial L." "beta no g"
matlist bonus, title("Alternative Sample Means and Ratios to Arithmetic Mean") ///
cspec(& %-10s w15 | %11.3f o4 & %9.6f w10 & %9.6f w10 & %9.6f o4 &) rspec(&|&&&&&|)
}
} /*end of svars*/
mata: scall = (ng, Leverage, pLeverage, betag)
/*display std errors for all coefficients*/
if "`regtable'" != "" {
/*the cv3J matrix is the square root already, CV3 is not*/
if "`jackknife'" == "" {
mata fullcv3 = diagonal(cv3)
mata: fullcv3 = sqrt(fullcv3)
mata: regresstab = J(rows(fullcv3),6,.)
mata: regresstab[.,1] =beta
mata: regresstab[.,2] = fullcv3
mata: tstats = beta :/ fullcv3
mata: regresstab[.,3] = tstats
mata: pvals = J(rows(fullcv3),2,.)
forvalues i = 1 / `sk' {
mata: pvals[`i',1] = ttail(`G'-1,regresstab[`i',3])
mata: pvals[`i',2] = 1 - ttail(`G'-1,regresstab[`i',3])
mata: regresstab[`i',4] = 2*min(pvals[`i',.])
}
mata: regresstab[.,5] = regresstab[.,1] - invttail(`G'-1,.025)* regresstab[.,2]
mata: regresstab[.,6] = regresstab[.,1] + invttail(`G'-1,.025)* regresstab[.,2]
mata: st_matrix("regresstab",regresstab)
/*dimensions of this matrix*/
mata: rtrows = rows(regresstab)
mata: st_numscalar("rtrows", rtrows)
mata: st_local("rtrows", strofreal(rtrows))
matrix colnames regresstab = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper
/*include constant row when absorb & fevar options not used */
if `"`absorb'"' == "" & `"`fevar'"' == "" {
matrix rownames regresstab = "`XTEMP' " "_cons" `name_xvar' `names_fe'
}
else {
matrix rownames regresstab = "`XTEMP' " `name_xvar' `names_fe'
}
if `rtrows' <= 50 {
local rspecreg: display "-" _dup(`rtrows') "&" "-"
matlist regresstab, title(Linear Regression -- CV3) rowtitle(`yvar') ///
cspec(& %-12s w6 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 & %9.6f w10 & %9.6f w10 &) ///
rspec(`rspecreg')
}
else{
disp "Linear Regression -- CV3"
matrix list regresstab
}
mata: st_matrix("cv3",cv3)
return matrix cv3 = cv3
} /*end of no jackknife if*/
if "`jackknife'" != "" {
mata fullcv3 = diagonal(cv3j)
mata: cv3jsq = cv3j:*cv3j
mata: st_matrix("cv3j",cv3jsq)
return matrix cv3j = cv3j
mata: regresstab = J(rows(fullcv3),6,.)
mata: regresstab[.,1] =beta
mata: regresstab[.,2] = fullcv3
mata: tstats = beta :/ fullcv3
mata: regresstab[.,3] = tstats
mata: pvals = J(rows(fullcv3),2,.)
forvalues i = 1 / `sk' {
mata: pvals[`i',1] = ttail(`G'-1,regresstab[`i',3])
mata: pvals[`i',2] = 1 - ttail(`G'-1,regresstab[`i',3])
mata: regresstab[`i',4] = 2*min(pvals[`i',.])
}
mata: regresstab[.,5] = regresstab[.,1] - invttail(`G'-1,.025)* regresstab[.,2]
mata: regresstab[.,6] = regresstab[.,1] + invttail(`G'-1,.025)* regresstab[.,2]
mata: st_matrix("regresstab",regresstab)
matrix colnames regresstab = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper
/*dimensions of this matrix*/
mata: rtrows = rows(regresstab)
mata: st_numscalar("rtrows", rtrows)
mata: st_local("rtrows", strofreal(rtrows))
/*include constant row when absorb & fevar options not used */
if `"`absorb'"' == "" & `"`fevar'"' == "" {
matrix rownames regresstab = "`XTEMP'" "_cons" `name_xvar' `names_fe'
}
else {
matrix rownames regresstab = "`XTEMP' " `name_xvar' `names_fe'
}
if `rtrows' <= 50 {
local rspecreg: display "-" _dup(`rtrows') "&" "-"
matlist regresstab, title(Linear Regression -- CV3) rowtitle(`yvar') ///
cspec(& %-12s w6 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 & %9.6f w10 & %9.6f w10 &) ///
rspec(`rspecreg')
}
else{
disp "Linear Regression -- CV3"
matrix list regresstab
}
} /*end of yes jackknife if */
/*singular clusters*/
if `singclustfull' >= 1 {
/*the cv3J matrix is the square root already, CV3 is not*/
if "`jackknife'" == "" {
mata fullcv3drop = diagonal(cv3drop)
mata: fullcv3drop = sqrt(fullcv3drop)
mata: regresstabdrop = J(rows(fullcv3drop),6,.)
mata: regresstabdrop[.,1] =beta
mata: regresstabdrop[.,2] = fullcv3drop
mata: tstatsdrop = beta :/ fullcv3drop
mata: regresstabdrop[.,3] = tstatsdrop
mata: pvalsdrop = J(rows(fullcv3drop),2,.)
forvalues i = 1 / `sk' {
mata: pvalsdrop[`i',1] = ttail(`gdrop'-1,regresstabdrop[`i',3])
mata: pvalsdrop[`i',2] = 1 - ttail(`gdrop'-1,regresstabdrop[`i',3])
mata: regresstabdrop[`i',4] = 2*min(pvalsdrop[`i',.])
}
mata: regresstabdrop[.,5] = regresstabdrop[.,1] - invttail(`gdrop'-1,.025)* regresstabdrop[.,2]
mata: regresstabdrop[.,6] = regresstabdrop[.,1] + invttail(`gdrop'-1,.025)* regresstabdrop[.,2]
mata: st_matrix("regresstabdrop",regresstabdrop)
/*dimensions of this matrix*/
mata: rtrows = rows(regresstabdrop)
mata: st_numscalar("rtrows", rtrows)
mata: st_local("rtrows", strofreal(rtrows))
matrix colnames regresstabdrop = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper
/*include constant row when absorb & fevar options not used */
if `"`absorb'"' == "" & `"`fevar'"' == "" {
matrix rownames regresstabdrop = "`XTEMP' " "_cons" `name_xvar' `names_fe'
}
else {
matrix rownames regresstabdrop = "`XTEMP' " `name_xvar' `names_fe'
}
if `rtrows' <= 50 {
local rspecreg: display "-" _dup(`rtrows') "&" "-"
matlist
regresstabdrop, title(Linear Regression -- CV3 -- Omitting Singular Omit-One-Cluster Subsamples) rowtitle(`yvar') ///
cspec(& %-12s w6 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 & %9.6f w10 & %9.6f w10 &) ///
rspec(`rspecreg')
}
else{
disp "Linear Regression -- CV3 -- Omitting Singular Omit-One-Cluster Subsamples"
matrix list regresstabdrop
}
mata: st_matrix("cv3drop",cv3drop)
return matrix cv3drop = cv3drop
} /* end of no jackknife if */
if "`jackknife'" != "" {
mata fullcv3drop = diagonal(cv3jdrop)
mata: cv3jsqdrop = cv3jdrop:*cv3jdrop
mata: st_matrix("cv3jdrop",cv3jsqdrop)
return matrix cv3jdrop = cv3jdrop
mata: regresstabdrop = J(rows(fullcv3),6,.)
mata: regresstabdrop[.,1] =beta
mata: regresstabdrop[.,2] = fullcv3drop
mata: tstatsdrop = beta :/ fullcv3drop
mata: regresstabdrop[.,3] = tstatsdrop
mata: pvalsdrop = J(rows(fullcv3drop),2,.)
forvalues i = 1 / `sk' {
mata: pvalsdrop[`i',1] = ttail(`G'-1,regresstabdrop[`i',3])
mata: pvalsdrop[`i',2] = 1 - ttail(`G'-1,regresstabdrop[`i',3])
mata: regresstabdrop[`i',4] = 2*min(pvals[`i',.])
}
mata: regresstabdrop[.,5] = regresstabdrop[.,1] - invttail(`G'-1,.025)* regresstabdrop[.,2]
mata: regresstabdrop[.,6] = regresstabdrop[.,1] + invttail(`G'-1,.025)* regresstabdrop[.,2]
mata: st_matrix("regresstabdrop",regresstabdrop)
matrix colnames regresstabdrop = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper
/*dimensions of this matrix*/
mata: rtrows = rows(regresstabdrop)
mata: st_numscalar("rtrows", rtrows)
mata: st_local("rtrows", strofreal(rtrows))
/*include constant row when absorb & fevar options not used */
if `"`absorb'"' == "" & `"`fevar'"' == "" {
matrix rownames regresstabdrop = "`XTEMP' " "_cons" `name_xvar' `names_fe'
}
else {
matrix rownames regresstabdrop = "`XTEMP' " `name_xvar' `names_fe'
}
if `rtrows' <= 50 {
local rspecreg: display "-" _dup(`rtrows') "&" "-"
matlist
regresstabdrop, title(Linear Regression -- CV3 -- Omitting Singular Omit-One-Cluster Subsamples) rowtitle(`yvar') ///
cspec(& %-12s w6 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 & %9.6f w10 & %9.6f w10 &) ///
rspec(`rspecreg')
}
else{
disp "Linear Regression -- CV3 -- Omitting Singular Omit-One-Cluster Subsamples"
matrix list regresstabdrop
}
}/*end of yes jackknife if*/
} /*end of singclustfull*/
} /*end of regtable*/
if "`table'"!="" {
if `G'< 53 {
mata: st_matrix("scall", scall)
matrix rownames scall = `CLUSTERNAMES'
matrix colnames scall = Ng Leverage "Partial L." "beta no g"
local TITLE "Cluster by Cluster Statistics"
local Gm1 = `G'- 1
local rspec : display "rspec(&-"_dup(`Gm1') "&" "-)"
matlist scall, title("`TITLE'") rowtitle("`cluster'") ///
cspec(& %-10s | %8.0f o4 & %9.6f o4 & %9.6f w10 & %9.6f o4 &) `rspec'
return matrix scall = scall
}
else {
disp " "
disp "Cluster by Cluster Statistics"
disp as error "The number of clusters, `G' exceeds 52, results reported as matrix. "
disp as text ""
mata: scall
}
} /*end of table*/
/*stuff to return*/
return matrix ng = ng
return matrix lever = Leverage
return matrix partlev = pLeverage
return matrix betanog = betag
/*nograph - default is to draw a graph*/
if "`nograph'" == "" {
tempvar ng
qui getmata `ng' = ng, force
label variable `ng' "Observations Per Cluster"
tempvar lever
qui getmata `lever' = Leverage, force
label variable `lever' "Leverage"
tempvar plever
qui getmata `plever' = pLeverage, force
label variable `plever' "Partial Leverage"
tempvar betanog
qui getmata `betanog' = betag, force
label variable `betanog' "Omit-One-Cluster Coefficients"
qui cap graph drop name summclust_temp_ln.gph
qui cap graph drop name summclust_temp_pn.gph
qui cap graph drop name summclust_temp_lb.gph
qui cap graph drop name summclust_temp_pb.gph
qui cap erase summclust_temp_ln.gph
qui cap erase summclust_temp_pn.gph
qui cap erase summclust_temp_lb.gph
qui cap erase summclust_temp_pb.gph
set graphics off
qui twoway scatter `lever' `ng', saving(summclust_temp_ln) scheme(s1mono) xline( `graph_1' ) yline( `graph_2' )
qui twoway scatter `plever' `ng', saving(summclust_temp_pn) scheme(s1mono) xline( `graph_1' ) yline(`graph_3' )
qui twoway scatter `lever' `betanog', saving(summclust_temp_lb) scheme(s1mono) xline( `graph_4' ) yline( `graph_2' )
qui twoway scatter `plever' `betanog', saving(summclust_temp_pb) scheme(s1mono) xline(`graph_4' ) yline( `graph_3' )
set graphics on
graph combine summclust_temp_ln.gph summclust_temp_pn.gph summclust_temp_lb.gph summclust_temp_pb.gph, title("Cluster Specific Statistics For `G' `cluster' Clusters") scheme(s1mono)
}
else if "`nograph'" != "" {
disp "Graph option suppressed."
}
/*additional stuff to return*/
/*scalars*/
return local beta = `BETATEMP'
/*cv1*/
return local cv1se = `setemp'
return local cv1t = `cv1t'
return local cv1p = `cv1p'
return local cv1lci = `cv1lci'
return local cv1uci = `cv1uci'
/*cv3*/
return local cv3se = `cv3se'
return local cv3t = `cv3t'
return local cv3p = `cv3p'
return local cv3lci = `cv3lci'
return local cv3uci = `cv3uci'
/*cv3drop*/
if `singclustfull' >= 1 {
return local cv3sedrop = `cv3sedrop'
return local cv3tdrop = `cv3tdrop'
return local cv3pdrop = `cv3pdrop'
return local cv3lcidrop = `cv3lcidrop'
return local cv3ucidrop = `cv3ucidrop'
}
/*jackknife*/
if "`jackknife'"!="" {
return local beta = `BETATEMP'
return local cv3Jse = `cv3jse'
return local cv3Jt = `cv3jt'
return local cv3Jp = `cv3jp'
return local cv3Jlci = `cv3jlci'
return local cv3Juci = `cv3juci'
if `singclustfull' >= 1 {
return local cv3sedrop = `cv3jsedrop'
return local cv3tdrop = `cv3jtdrop'
return local cv3pdrop = `cv3jpdrop'
return local cv3lcidrop = `cv3jlcidrop'
return local cv3ucidrop = `cv3jucidrop'
}
}
/*gstar*/
if `gstarind'==1 {
return local gstarzero = `gstarzero'
return local gstarone = `gstarone'
return local gstarrho = `gstarrho'
}
else if `gstarind'==2 {
return local gstarzero = `gstarzero'
}
end
/*--------------------------------------*/
/* Change Log */
/*--------------------------------------*/
*1.004 - mata cv table, CV3J addition, absorb check, return stuff
*1.005 - sample option, make gstar locals conditional
*1.006 - gstar1 caution, large table as matrix
*1.007 - geometric mean
*1.008 - absorb-jack conflict bug corrected
*1.009 - added version number for ssc submission
*1.010 - check for singularities and string cluster variable correction
*2.000 - replace globals with locals, tempvars, absorb-treat conflict resolution
*2.001 - additional global replacements, gstar formatting
*2.002 - fixed constant issue, d.o.f. for cv1
*2.003 - small formatting changes to tables
*3.000 - added regtable and nograph options, fixed bugs
*3.001 - fixed absorb errors
*4.000 - changed handling for singularities
*4.001 - changed absorb calculations to panelsum
*4.002 - debugged absorb and regtable errors
*4.003 - debugged marksample errors
*4.100 - debugged additional absorb error
*4.200 - debugged singular flag error
*4.201 - debugged singular flag absorb error
*4.202 - additional returned values
*4.203 - "singular subsample" terminology used in output
*4.204 - several fixes
*4.208 - cleaned up singular clusters warnings
*4.210 - remove Lg from output when absorb not properly used