/*------------------------------------*/
/*twowayjack */
/*written by Matt Webb */
/*version 1 2024-06-11 */
/*------------------------------------*/
version 13

/*----------------------------------*/

/*two-way jack program */

*this makes three calls to summclust
*uses the max operator
*reports cv1 and cv3 multi-way standard errors
/*----------------------------------*/

cap program drop twowayjack
program define twowayjack, rclass

  syntax varlist(min=3), cluster(varlist) [fevar(varlist) SAMple(string)]

  local XTEMP = word("`varlist'",2)

  /*check if two-way?*/
  local numdims = wordcount("`cluster'")
  else if `numdims' > 2 {
    disp as error "twowayjack only supports two-way clustering at the moment."
    exit 198
  }

  else if `numdims' == 2 {

    local clust1  = word("`cluster'",1)

    local clust2  = word("`cluster'",2)

  }

  else if `numdims' == 1{
    disp as error "For a single dimension of clustering use summclust"
    exit 198
  }

  /*add if to sample*/
  if "`sample'" != "" {
      local SMPL = "sample(`sample')"
  }

  /*matrix to store output*/
  *summtable
  mata: summtable = J(3,6,.)
  mata: numclust = J(1,3,.)

  *first dimension
   cap quietly mwsummclust `varlist', cluster("`clust1'") fevar(`fevar') ///
   `SMPL' nog gstar

   mata: cv1g = cv1se*cv1se
   mata: cv3g = cv3se*cv3se

   mata: summtable[1,1::4] = clustsum[7,1::4]
   mata: summtable[1,5] = rows(betag)
   mata: summtable[1,6] = gstarzero

   mata: numclust[1,1] = rows(betag)

  *second dimension
   cap quietly mwsummclust `varlist', cluster("`clust2'") fevar(`fevar') ///
   `SMPL' nog gstar

   mata: cv1h = cv1se*cv1se
   mata: cv3h = cv3se*cv3se

   mata: summtable[2,1::4] = clustsum[7,1::4]
   mata: summtable[2,5] = rows(betag)
   mata: summtable[2,6] = gstarzero

   mata: numclust[1,2] = rows(betag)

  *intersection

  tempvar ints

  egen `ints'  = group(`clust1' `clust2')

  cap quietly mwsummclust `varlist', cluster(`ints') fevar(`fevar') ///
  `SMPL' nog gstar

   mata: cv1i = cv1se*cv1se
   mata: cv3i = cv3se*cv3se

   mata: summtable[3,1::4] = clustsum[7,1::4]
   mata: summtable[3,5] = rows(betag)
   mata: summtable[3,6] = gstarzero

   mata: numclust[1,3] = rows(betag)

  *max

  mata: cv1trad = cv1g + cv1h - cv1i
  mata: cv3trad = cv3g + cv3h - cv3i

  mata: cv1set= (cv1trad,cv1g,cv1h)
  mata: cv3set= (cv3trad,cv3g,cv3h)


  mata: cv1max = sqrt(max(cv1set))
  *mata: cv1set
  *mata cv1max

  mata: cv3max = sqrt(max(cv3set))
  *mata: cv3set
  *mata cv3max

  mata: dfmax = min(numclust)

  *regression output
  *regoutput
  mata: mwcvstuff = J(2,6,.)

  mata mwcvstuff[1,1] = cvstuff[1,1]
  mata mwcvstuff[1,2] = cv1max
  mata mwcvstuff[1,3] = cvstuff[1,1] / cv1max
  mata mwcvstuff[1,4] = ///
  2*min((ttail(dfmax-1,mwcvstuff[1,3]),1-ttail(dfmax-1,mwcvstuff[1,3])))
  mata mwcvstuff[1,5] = cvstuff[1,1] - invttail(rows(betag)-1,0.025) * cv1max
  mata mwcvstuff[1,6] = cvstuff[1,1] + invttail(rows(betag)-1,0.025) * cv1max

  mata mwcvstuff[2,1] = cvstuff[1,1]
  mata mwcvstuff[2,2] = cv3max
  mata mwcvstuff[2,3] = cvstuff[1,1] / cv3max
  mata mwcvstuff[2,4] = ///
  2*min((ttail(dfmax-1,mwcvstuff[2,3]),1-ttail(dfmax-1,mwcvstuff[2,3])))
  mata mwcvstuff[2,5] = cvstuff[1,1] - invttail(dfmax-1,0.025) * cv3max
  mata mwcvstuff[2,6] = cvstuff[1,1] + invttail(dfmax-1,0.025) * cv3max

  mata: st_matrix("mwcvstuff", mwcvstuff)

  *report results

  *matrix list mwcvstuff

  /*output display*/
  disp ""
  disp ""
  disp "TWOWAYJACK"
  disp "Reference: James G. MacKinnon, Morten Ø. Nielsen, and Matthew D. Webb"
  disp "           Jackknife inference with two-way clustering"
  disp " "
  disp "Two-way cluster jackknife variance estimation."
  disp "Cluster summary statistics for `XTEMP' when clustered by `clust1' and `clust2'."
  *disp "There are `numobs' observations within `G' `cluster' clusters."

  local RSPECCV "&|&|"
  matrix colnames mwcvstuff = "Coeff" "Sd. Err." "t-stat" "P value" ///
  CI-lower CI-upper
  matrix rownames mwcvstuff = "CV1max" "CV3max"

  matlist mwcvstuff, title(Regression Output) rowtitle(s.e.) ///
  cspec(& %-4s w8 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 & %9.6f w10 ///
  & %9.6f w10 &) rspec("`RSPECCV'")

  mata: st_matrix("summtable",summtable)

  matrix colnames summtable = Ng Leverage "Partial L." "beta no g" "G" "Gstar"
  matrix rownames summtable = "`clust1'" "`clust2'" "intersect"
  matlist summtable, title(Coefficients of Variation, G, and G*) ///
  rowtitle(dimension) cspec(& %-4s w10 | %7.4f w8 & %7.4f w9 & %7.4f w10 ///
  & %7.4f w10 & %6.0f w7 & %8.2f w9 &) rspec("&|&&|")

end


/*summclust*/

/*------------------------------------*/
/*summclust */
/*written by Matt Webb */
/*version 4.210 2023-06-06 */
/*------------------------------------*/
version 13

cap program drop mwsummclust
program define mwsummclust, 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.1 - added summclust code