/*-----------------------------*/
/*ADO file*/
/* logitjack  - logit jack(knife) estimation */
/* implements procedures found in this paper: */
/* logit-jack */
/* written by Matt Webb */
/* version 3.3 */
/* date 05/29/24 */
/*-----------------------------*/


capture program drop logitjack
program define logitjack, rclass
	syntax varlist(min = 2) , CLUSter(varname) [fevar(varlist) VARiables(real 1) BOOTstrap reps(real 999)  JACKknife  SAMple(string) NOnull ]

	version 13

	local SMPL "`sample'"

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

	/*flag for linearized or boot*/
	local flagboot = 0
	if "`bootstrap'" != "" {
		local flagboot = 1
	}
	if "`nonull'" != "" {
		local flagboot = 1
	}
	if `reps' != 999 {
		local flagboot = 1
	}

	local flaglinear = 0
	if "`linear'"!=""{
		local flaglinear = 1
	}
	local flagboth = `flagboot'+`flaglinear'

	mata cvstuff = J(1,6,.)
	mata linstuff = J(1,6,.)
	mata bootstuff = J(2,4,.)
	mata: bootci = J(4,4,.)

	/*fevar setup*/
	local FEVAR "`fevar'"
	local countfevar = 0
	local catvars = " "
	foreach catvar in `FEVAR' {
	    local catvars = " `catvars' " +  "i.`catvar'"
	}

	/*sort the data by the clustering variable*/
		/*make a sequential-numeric cluster indicator*/

	tempvar temp_sample
	mark `temp_sample'
	markout `temp_sample' `varlist' `absorb' `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'
	}

	/*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: st_numscalar("numobs", numobs)
	mata: st_local("numobs", strofreal(numobs))

	mata: info = panelsetup(CVAR,1)
	mata: ng = info[.,2] - info[.,1] :+ 1
	mata: G = rows(ng)

	if `variables' == 1{
		local x = word("`varlist'",2)
		local y = word("`varlist'",1)
	}

	*to impose null
	if "`nonull'" == "" {
		local newvarlist = subinstr("`varlist'","`x'","",1)
	}

	/*when fevar specified*/
	if "`fevar'" != ""{

		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 */
			local jend = `j'-1

			local fevars = "`fevars' `t_`jstart'' - `t_`jend''  "

		} /*end of fvar*/

		local allti = "`fevars'"

		qui logit `varlist' `allti',  cluster(`cluster')

		local beta = _b[`x']
		mata: beta = `beta'

		local se =  _se[`x']
		mata: cv1se = `se'

		local yhat = "yhat"
		tempvar yhat
		qui predict `yhat' if `temp_sample'==1

		qui putmata yhat = `yhat' if `temp_sample'==1, replace

		if "`nonull'" == ""{

			qui logit `newvarlist' `allti',  cluster(`cluster')

			local yhatr = "yhatr"
			tempvar yhatr
			qui predict `yhatr' if `temp_sample'==1

			qui putmata yhatr = `yhatr' if `temp_sample'==1, replace

			/*dump everything to mata*/
			qui putmata ALLR = (`newvarlist' `allti') if `temp_sample'==1, replace

			mata: ones = J(rows(ALLR),1,1)

			mata ALLR = ALLR, ones

			mata numall = cols(ALLR)
			mata XR = ALLR[.,2::numall]

		}

		/*dump everything to mata*/
		qui putmata ALL = (`varlist' `allti') if `temp_sample'==1, replace

		mata: ones = J(rows(ALL),1,1)

		mata ALL = ALL, ones

		mata numall = cols(ALL)
		mata Y = ALL[.,1]
		mata X = ALL[.,2::numall]

		mata: cv1se = sqrt(cv1se*cv1se* (((rows(X)-1)/(rows(X)-cols(X)))))

		mata: st_numscalar("cv1se", cv1se)
		mata: st_local("cv1se", strofreal(cv1se))

		* correct the factor from Stata's output
		mata: tcv1 = beta/cv1se
		mata: st_numscalar("tcv1", tcv1)
		mata: st_local("tcv1", strofreal(tcv1))

		local p1 = 2*min(ttail(`G'-1,`tcv1'), 1 -ttail(`G'-1,`tcv1'))

		local cv1lci = `beta' - invttail(`G'-1,0.025) * `cv1se'
		local cv1uci = `beta' + invttail(`G'-1,0.025) * `cv1se'

		mata: cvstuff[1,1] = `beta'
		mata: cvstuff[1,2] = `cv1se'
		mata: cvstuff[1,3] = `tcv1'
		mata: cvstuff[1,4] = `p1'
		mata: cvstuff[1,5] = `cv1lci'
		mata: cvstuff[1,6] = `cv1uci'


		/*brute force jackknife*/

		if "`jackknife'" != "" {

			mata: betas= J(`G',1,.)

			forvalues g = 1 / `G' {

				 qui logit `varlist' `allti' if `temp_indexa'!=`g',  cluster(`cluster')
				 local betas = _b[`x']
				 mata: betas[`g',1]= `betas'

			} /*end of j*/

			mata: cv3diff = betas:-beta
			mata: cv3diffsqsum = (`G'-1)/(`G')* sum(cv3diff:*cv3diff)
			mata: cv3se= sqrt(cv3diffsqsum)
			mata: st_numscalar("cv3se", cv3se)
			mata: st_local("cv3se", strofreal(cv3se))

			local tcv3 = `beta'/`cv3se'
			mata: tcv3 = beta/cv1se

			local p3 = 2*min(ttail(`G'-1,`tcv3'), 1 -ttail(`G'-1,`tcv3'))

			local cv3lci = `beta' - invttail(`G'-1,0.025) * `cv3se'
			local cv3uci = `beta' + invttail(`G'-1,0.025) * `cv3se'

			mata: cvstuffjack = J(1,6,.)
			mata: cvstuffjack[1,1] = `beta'
			mata: cvstuffjack[1,2] = `cv3se'
			mata: cvstuffjack[1,3] = `tcv3'
			mata: cvstuffjack[1,4] = `p3'
			mata: cvstuffjack[1,5] = `cv3lci'
			mata: cvstuffjack[1,6] = `cv3uci'

		} /* end of jackknife if */

	}/*end of fevar*/

	/*when absorb not specified*/
	if "`fevar'" == ""{

		*estimate logit
		qui logit `varlist', cluster(`cluster')

		local beta = _b[`x']
		mata: beta = `beta'

		local se =  _se[`x']

		local yhat = "yhat"
		tempvar yhat
		qui predict `yhat' if `temp_sample'==1

		qui putmata yhat = `yhat' if `temp_sample'==1, replace

		/*dump everything to mata*/
		qui putmata ALL = (`varlist') if `temp_sample'==1, replace

		*add a constant to X
		qui mata: ones = J(rows(ALL),1,1)

		mata numall = cols(ALL)
		mata Y = ALL[.,1]
		mata X = (ALL[.,2::numall], ones)

		mata: cv1se = `se'
		mata: cv1se = sqrt(cv1se*cv1se* (((rows(X)-1)/(rows(X)-cols(X)))))

		mata: st_numscalar("cv1se", cv1se)
		mata: st_local("cv1se", strofreal(cv1se))

		* correct the factor from Stata's output
		mata: tcv1 = beta/cv1se
		mata: st_numscalar("tcv1", tcv1)
		mata: st_local("tcv1", strofreal(tcv1))

		local p1 = 2*min(ttail(`G'-1,`tcv1'), 1 -ttail(`G'-1,`tcv1'))

		local cv1lci = `beta' - invttail(`G'-1,0.025) * `cv1se'
		local cv1uci = `beta' + invttail(`G'-1,0.025) * `cv1se'

		mata: cvstuff[1,1] = `beta'
		mata: cvstuff[1,2] = `cv1se'
		mata: cvstuff[1,3] = `tcv1'
		mata: cvstuff[1,4] = `p1'
		mata: cvstuff[1,5] = `cv1lci'
		mata: cvstuff[1,6] = `cv1uci'


		if "`nonull'" == "" {
			*estimate logit
			qui logit `newvarlist', cluster(`cluster')

			local yhatr = "yhatr"
			tempvar yhatr
			qui predict `yhatr' if `temp_sample'==1

			qui putmata yhatr = `yhatr' if `temp_sample'==1, replace

			/*dump everything to mata*/
			qui putmata ALLR = (`newvarlist') if `temp_sample'==1, replace

			*add a constant to X
			qui mata: ones = J(rows(ALLR),1,1)

			mata numall = cols(ALLR)
			mata XR = (ALLR[.,2::numall], ones)

		} /*end of nonull */

		if "`jackknife'" != ""  {

			*bruteforce jackknife
			mata: betas= J(`G',1,.)

			forvalues g = 1 / `G' {

				 qui logit `varlist' if `temp_indexa'!=`g',  cluster(`cluster')
				 local betas = _b[`x']
				 mata: betas[`g',1]= `betas'

			} /*end of g*/

			mata: cv3diff = betas:-beta
			mata: cv3diffsqsum = (`G'-1)/(`G')* sum(cv3diff:*cv3diff)
			mata: cv3se= sqrt(cv3diffsqsum)
			mata: st_numscalar("cv3se", cv3se)
			mata: st_local("cv3se", strofreal(cv3se))

			local tcv3 = `beta'/`cv3se'
			mata: tcv3 = beta/cv1se

			local p3 = 2*min(ttail(`G'-1,`tcv3'), 1 -ttail(`G'-1,`tcv3'))

			local cv3lci = `beta' - invttail(`G'-1,0.025) * `cv3se'
			local cv3uci = `beta' + invttail(`G'-1,0.025) * `cv3se'

			mata: cvstuffjack = J(1,6,.)
			mata: cvstuffjack[1,1] = `beta'
			mata: cvstuffjack[1,2] = `cv3se'
			mata: cvstuffjack[1,3] = `tcv3'
			mata: cvstuffjack[1,4] = `p3'
			mata: cvstuffjack[1,5] = `cv3lci'
			mata: cvstuffjack[1,6] = `cv3uci'

		} /* end of jackknife if */

	} /*end of no-fevar*/

		mata: scores = (Y :- yhat):*X
		mata: score_g =  panelsum(scores, info)

		if "`nonull'" == "" {

			mata: scoresr = (Y :- yhatr):*XR
			mata: score_rg =  panelsum(scoresr, info)
			mata: scoresr2 = (Y :- yhatr):*X
			mata: score_r2g =  panelsum(scoresr2, info)
		}

		forvalues g = 1 / `G' {

			mata: firstrow = info[`g',1]
			mata: st_numscalar("firstrow", firstrow)
			mata: st_local("firstrow", strofreal(firstrow))
			mata: lastrow = info[`g',2]
			mata: st_numscalar("lastrow", lastrow)
			mata: st_local("lastrow", strofreal(lastrow))

			local f = `firstrow'
			mata: infomat_g`g' = yhat[`f',1]*(1 - yhat[`f',1]):*X[`f',.]'X[`f',.]

			if "`nonull'" == "" {
				mata: infomat_rg`g' = yhatr[`f',1]*(1 - yhatr[`f',1]):*XR[`f',.]'XR[`f',.]

				mata: infomat_r2g`g' = yhatr[`f',1]*(1 - yhatr[`f',1]):*X[`f',.]'X[`f',.]
			}

			local next = `firstrow'+1
			forvalues i = `next'/`lastrow'{

				mata: infomat_g`g' = infomat_g`g' + yhat[`i',1]*(1 - yhat[`i',1]):*X[`i',.]'X[`i',.]

				if "`nonull'" == "" {

					mata: infomat_rg`g' = infomat_rg`g' + yhatr[`i',1]*(1 - yhatr[`i',1]):*XR[`i',.]'XR[`i',.]

					mata: infomat_r2g`g' = infomat_r2g`g' + yhatr[`i',1]*(1 - yhatr[`i',1]):*X[`i',.]'X[`i',.]

				}

			} /*end of i*/

		}/*end of g*/

		mata: score_all = colsum(score_g)
		mata: infomat_all = infomat_g1

		forvalues g = 2/`G'{
			mata: infomat_all = infomat_all + infomat_g`g'
		}

		mata: invinfomat_all = invsym(infomat_all)

		mata: linbeta = invinfomat_all*score_all'
		mata: linbetax = linbeta[1,1]

		mata: st_numscalar("linbetax", linbetax)
		mata: st_local("linbetax", strofreal(linbetax))

		*if "`nonull'" != ""{
			mata: score_us = J(`G',cols(X),.)
		*}
		if "`nonull'" == ""{
			mata: score_rall = colsum(score_rg)
			mata: score_r2all = colsum(score_r2g)

			mata: infomat_rall = infomat_rg1
			mata: infomat_r2all = infomat_r2g1

			forvalues g = 2/`G'{
				mata: infomat_rall = infomat_rall + infomat_rg`g'
				mata: infomat_r2all = infomat_r2all + infomat_r2g`g'
			}

			mata: invinfomat_rall = invsym(infomat_rall)
			mata: invinfomat_r2all = invsym(infomat_r2all)

			mata: nm1 =cols(infomat_all)
			mata: st_numscalar("nm1", nm1)
			mata: st_local("nm1", strofreal(nm1))

			mata: linbetar = invinfomat_rall*score_rall'

			mata: score_rs = J(`G',cols(X)-1,.)
			mata: score_r2s = J(`G',cols(X),.)

		} /*end of nonull*/

		/*-----------------------*/
		/*cv3l*/
		/*-----------------------*/

		mata: cv3lsum = J(rows(linbeta),rows(linbeta),0)

		mata: betaog = J(`G',1,.)

		forvalues g = 1 / `G' {

			*unrestricted
			mata: beta_o`g' = invsym(infomat_all - infomat_g`g')*(score_all'-score_g[`g',.]')

			mata: cv3lsum = cv3lsum + (beta_o`g':-linbeta)*(beta_o`g'':-linbeta')

			*transform the scores

				*unrestricted - equation 33

				mata: score_us[`g',.] = score_g[`g',.] - (infomat_g`g'* beta_o`g')'

			if "`nonull'" == ""{

				*restricted  - equation 34
				mata: beta_ro`g' = invsym(infomat_rall - infomat_rg`g')*(score_rall'-score_rg[`g',.]')
				mata: score_rs[`g',.] = score_rg[`g',.] - (infomat_rg`g'* beta_ro`g')'

				mata: score_r2s[`g',.] = score_r2g[`g',.] - (infomat_r2g`g'[.,2::cols(infomat_r2g`g')]* beta_ro`g')'

			}

			mata: betaog[`g',1] = 	beta_o`g'[1,1]

		} /*end of g*/

		mata: cv3l = ((`G'-1)/`G'):*cv3lsum

		mata: cv3lse = sqrt(cv3l[1,1])
		mata: st_numscalar("cv3lse", cv3lse)
		mata: st_local("cv3lse", strofreal(cv3lse))

		mata: tcv3l = `beta' / cv3lse

		mata: st_numscalar("tcv3l", tcv3l)
		mata: st_local("tcv3l", strofreal(tcv3l))

		local p3l = 2*min(ttail(`G'-1,`tcv3l'), 1 -ttail(`G'-1,`tcv3l'))

		local cv3llci = `beta' - invttail(`G'-1,0.025) * `cv3lse'
		local cv3luci = `beta' + invttail(`G'-1,0.025) * `cv3lse'

		mata: linstuff[1,1] = `beta'
		mata: linstuff[1,2] = cv3lse
		mata: linstuff[1,3] = tcv3l
		mata: linstuff[1,4] = `p3l'
		mata: linstuff[1,5] = `cv3llci'
		mata: linstuff[1,6] = `cv3luci'

		/*summary statistics*/

		if "`jackknife'" == "" {
			local SUMVAR "ng betalg "

			mata: clustsum= J(7,2,.)
		}

		if "`jackknife'" != "" {
			local SUMVAR "ng betalg betas"

			mata: clustsum= J(7,3,.)
		}

		local s = 0

		mata: betalg = betaog:+beta

		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


		  } /*end of svars if*/

		 	mata: st_matrix("clustsum", clustsum)
			matrix rownames clustsum = min q1 median mean q3 max coefvar

			if "`jackknife'" == "" {
				matrix colnames clustsum = Ng "Lin beta no g"
			}
			else {

				matrix colnames clustsum = Ng "Lin beta no g" "beta no g"
			}



	*} /*end of flag*/

	/*----------------*/
	/* bootstrap */
	/*----------------*/

	if `flagboot' == 1 {

		mata: w_sum = J(cols(score_g),cols(score_g),0)

		forvalues g = 1 / `G'{
			mata: w_g = score_g[`g',.]' - infomat_g`g' * linbeta
			mata: w_sum = w_sum + w_g*w_g'
		}

		*variance
		mata: var_cv1 = ((`G')/(`G'-1)):* invinfomat_all * w_sum * invinfomat_all

		mata: t = linbeta[1,1]/sqrt(var_cv1[1,1])

		mata: tboots = J(`reps',1,.)
		mata: tbootss = J(`reps',1,.)
		mata: betabootc = J(`reps',1,.)
		mata: betaboots = J(`reps',1,.)

		mata: tbootnewc = J(`reps',1,.)
		mata: tbootnews = J(`reps',1,.)


		forvalues b = 1 / `reps' {

			*switch to six-point if G is less than 12

			if `G' > 12 {

				*draw a G-vector of rademacher weights
				mata  e = runiform(`G',1)

				mata v = 2:* (e :> 0.5) :- 1

				local BOOTSTRING "P-values calculated with `reps' replications and Rademacher weights."
			}
			if `G' <= 12 {

				*draw a G-vector of six point weights
				mata  e = runiform(`G',1)
				mata: st_matrix("e", e)

				matrix vh = J(`G',1,.)

				forvalues h = 1 / `G' {

					if inrange(e[`h',1],0/6,6/6) {
						matrix vh[`h',1] = -sqrt(3/2)
					}
					if inrange( e[`h',1],1/6,2/6){
						matrix vh[`h',1] = -sqrt(2/2)
					}
					if inrange( e[`h',1],2/6,3/6) {
						matrix vh[`h',1] = -sqrt(1/2)
					}
					if inrange( e[`h',1],3/6,4/6) {
						matrix vh[`h',1] =  sqrt(1/2)
					}
					if inrange( e[`h',1],4/6,5/6) {
						matrix vh[`h',1] =  sqrt(2/2)
					}
					if inrange( e[`h',1],5/6,6/6) {
						matrix vh[`h',1] =  sqrt(3/2)
					}

				}
				mata v = st_matrix("vh")

				local BOOTSTRING "P-values calculated with `reps' replications and Webb weights."


			} /*end of six point*/

			/*unrestricted*/
			if "`nonull'" != "" {

				******************************
				/*classic - C*/

				*transform the scores, then sum them
				mata: sum_score_b = colsum(v:*score_g)

				*bootstrap beta
				mata beta_boot = invinfomat_all * sum_score_b'
				mata: betabootc[`b',1] = beta_boot[1,1]

				*bootstrap empirical scores
				mata: w_sum = J(cols(score_g),cols(score_g),0)
				forvalues g = 1 / `G'{
					mata: w_g = v[`g',1]:*score_g[`g',.]' - infomat_g`g' * beta_boot
					mata: w_sum = w_sum + w_g*w_g'
				}

				*bootstrap variance
				 mata: var_boot = ((`G')/(`G'-1))*((rows(X)-1)/(rows(X)-cols(X))):* invinfomat_all * w_sum * invinfomat_all

				*bootstrap t
				mata: tboots[`b',.] = beta_boot[1,1]/sqrt(var_boot[1,1])

				********************
				/*scores - S*/
				*transform the scores, then sum them
				mata: sum_score_sb = colsum(v:*score_us)

				*bootstrap beta
				mata beta_boot = invinfomat_all * sum_score_sb'
				mata: betaboots[`b',1] = beta_boot[1,1]

				*bootstrap empirical scores
				mata: w_sum = J(cols(score_g),cols(score_g),0)
				forvalues g = 1 / `G'{
					mata: w_g = v[`g',1]:*score_us[`g',.]' - infomat_g`g' * beta_boot
					mata: w_sum = w_sum + w_g*w_g'
				}

				*bootstrap variance
				 mata: var_boot = ((`G')/(`G'-1))*((rows(X)-1)/(rows(X)-cols(X))):* invinfomat_all * w_sum * invinfomat_all

				*bootstrap t
				mata: tbootss[`b',.] = beta_boot[1,1]/sqrt(var_boot[1,1])

			} /* end of unrestricted */

			********************
			/*restricted*/
			if "`nonull'" == "" {

				/*classic - C*/
				*transform the scores, then sum them
				mata: sum_score_b = colsum(v:*score_r2g)

				*bootstrap beta
				mata beta_boot = invinfomat_r2all * sum_score_b'

				*beta - for variance
					mata: sum_score_b = colsum(v:*score_r2g)
					mata beta_boot_var = invinfomat_r2all * sum_score_b'

				*bootstrap empirical scores
				mata: w_sum = J(cols(score_r2g),cols(score_r2g),0)
				forvalues g = 1 / `G'{
					mata: w_g = v[`g',1]:*score_r2g[`g',.]' - infomat_r2g`g' * beta_boot_var
					mata: w_sum = w_sum + w_g*w_g'
				}

				*bootstrap variance
				 mata: var_boot = ((`G')/(`G'-1))*((rows(X)-1)/(rows(X)-cols(X))):* invinfomat_r2all * w_sum * invinfomat_r2all

				*bootstrap t
				mata: tbootnewc[`b',.] = beta_boot[1,1]/sqrt(var_boot[1,1])

				************************
				/*score - S*/
				*transform the scores, then sum them
				mata: sum_score_sb = colsum(v:*score_r2s)

				*bootstrap beta
				mata beta_boot = invinfomat_r2all * sum_score_sb'

				*beta for variance
					mata: sum_score_sb = colsum(v:*score_r2s)
					mata beta_boot_var = invinfomat_r2all * sum_score_sb'

				*bootstrap empirical scores
				mata: w_sum = J(cols(score_r2g),cols(score_r2g),0)
				forvalues g = 1 / `G'{
					mata: w_g = v[`g',1]:*score_r2s[`g',.]' - infomat_r2g`g' * beta_boot_var
					mata: w_sum = w_sum + w_g*w_g'
				}

				*bootstrap variance
				 *mata: var_boot = ((`G'-1)):* invinfomat_rall * w_sum * invinfomat_rall
				 mata: var_boot = ((`G')/(`G'-1))*((rows(X)-1)/(rows(X)-cols(X))):* invinfomat_r2all * w_sum * invinfomat_r2all

				*bootstrap t
				mata: tbootnews[`b',.] = beta_boot[1,1]/sqrt(var_boot[1,1])


			} /* end of restricted */

		} /*end of bootstrap*/

		*bootstrap p

		*classic
		mata: bootc_p = mean(abs(tcv1[1,1]):<abs(tboots))

		mata: st_numscalar("bootc_p", bootc_p)
		mata: st_local("bootc_p", strofreal(bootc_p))

		*scores
		mata: boots_p = mean(abs(tcv1[1,1]):<abs(tbootss))

		mata: st_numscalar("boots_p", boots_p)
		mata: st_local("boots_p", strofreal(boots_p))

		*new
		mata: bootnews_p = mean(abs(tcv1[1,1]):<abs(tbootnews))
		mata: bootnewc_p = mean(abs(tcv1[1,1]):<abs(tbootnewc))

		mata: bootstuff[1,1] = `beta'
		mata: bootstuff[1,2] = `cv1se'
		mata: bootstuff[1,3] = `tcv1'

		mata: bootstuff[2,1] = `beta'
		mata: bootstuff[2,2] = `cv1se'
		mata: bootstuff[2,3] = `tcv1'

		if  "`nonull'" != "" {
			mata: bootstuff[1,4] = bootc_p
			mata: bootstuff[2,4] = boots_p
		}

		if  "`nonull'" == "" {
			mata: bootstuff[1,4] = bootnewc_p
			mata: bootstuff[2,4] = bootnews_p
		}


		/*calculate bootstrap confidence intervals*/
		if "`nonull'" != "" {

			*bootstrap quantiles
			mata: qlower = ceil(0.025*rows(tboots))
			mata: qupper = ceil(0.975*rows(tboots))

			mata: sortc = sort(tboots,1)
			mata: sorts = sort(tbootss,1)

			mata: tlowerc = sortc[qlower,1]
			mata: tupperc = sortc[qupper,1]

			mata: tlowers = sorts[qlower,1]
			mata: tuppers = sorts[qupper,1]


			mata: cilowc = `beta' - abs(tlowerc)*`cv1se'
			mata: ciupc = `beta' + abs(tupperc)*`cv1se'

			mata: cilows = `beta' - abs(tlowers)*`cv1se'
			mata: ciups = `beta' + abs(tuppers)*`cv1se'


			*second sets of CIs
			mata: bootsec = sqrt(1/(rows(betabootc)-1)*sum((betabootc:-mean(betabootc)):*(betabootc:-mean(betabootc))))

			mata: bootses = sqrt(1/(rows(betaboots)-1)*sum((betaboots:-mean(betaboots)):*(betaboots:-mean(betaboots))))

			local critval = invttail(`G'-1,0.025)

			mata: cilowbootc = `beta' - `critval'*bootsec
			mata: ciupbootc = `beta' + `critval'*bootsec

			mata: cilowboots = `beta' - `critval'*bootses
			mata: ciupboots = `beta' + `critval'*bootses

			*"CLASSIC-CV1-se"
			mata: bootci[1,1] = `beta'

			mata: bootci[1,2] = `cv1se'

			mata: bootci[1,3] = cilowc

			mata: bootci[1,4] = ciupc

			*"CLASSIC-boot-se"
			mata: bootci[2,1] = `beta'

			mata: bootci[2,2] = bootsec

			mata: bootci[2,3] = cilowbootc

			mata: bootci[2,4] = ciupbootc

			*"SCORE-CV1-se"
			mata: bootci[3,1] = `beta'

			mata: bootci[3,2] =  `cv1se'

			mata: bootci[3,3] = cilows

			mata: bootci[3,4] = ciups

			*"SCORE-boot-se"
			mata: bootci[4,1] = `beta'

			mata: bootci[4,2] = bootses

			mata: bootci[4,3] = cilowboots

			mata: bootci[4,4] = ciupboots

		}


	} /*end of boots if*/


	/*output display*/
		disp ""
		disp ""
		disp "LOGITJACK - MacKinnon, Nielsen, and Webb"
		disp " "
		disp "Jackknife cluster statistics for binary response models."
		disp "Estimates for `x' when clustered by `cluster'."
		disp "There are `numobs' observations within `G' `cluster' clusters."

	if "`jackknife'"!="" {

		mata: cvtable  = cvstuff \ cvstuffjack \ linstuff


		mata: st_matrix("cvtable",cvtable)
		matrix rownames cvtable = CV1 CV3 CV3L
		local RSPECCV "&|&&|"

		matrix colnames cvtable = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper

		matlist cvtable, title(Logistic 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 "`jackknife'"=="" {

		mata: cvtable  = cvstuff \ linstuff


		mata: st_matrix("cvtable",cvtable)
		matrix rownames cvtable = CV1  CV3L
		local RSPECCV "&|&|"

		matrix colnames cvtable = "Coeff" "Sd. Err." "t-stat" "P value" CI-lower CI-upper

		matlist cvtable, title(Logistic 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'")
	}

		/*summary statistic output*/
		if "`jackknife'"=="" {

			matlist clustsum , title("Cluster Variability") rowtitle(Statistic) ///
			cspec(& %-10s | %8.2f o2 & %9.6f  w14 & ) ///
			rspec(&-&&&&&-&)
		}
		else  {
			matlist clustsum , title("Cluster Variability") rowtitle(Statistic) ///
			cspec(& %-10s | %8.2f o4 & %9.6f  w14 & %9.6f  w10 &) ///
			rspec(&-&&&&&-&)
		}



	if `flagboot'== 1{

		if "`nonull'" == "" {
			local BOOTTITLE "Restricted Bootstrapped Linearized Regression Output"
			local BOOTROW "WCLR"
		}

		if "`nonull'" != ""{
			local BOOTTITLE "Unrestricted Bootstrapped Linearized Regression Output"
			local BOOTROW "WCLU"
		}

		mata: st_matrix("bootstuff",bootstuff)

		matrix rownames bootstuff = CLASSIC SCORE
		local RSPECBOOT "&|&|"

		matrix colnames bootstuff = "Coeff" "Sd. Err." "t-stat" "P value"

		matlist bootstuff, title(`BOOTTITLE') rowtitle(`BOOTROW') ///
			cspec(& %-10s w6 | %9.6f w10 & %9.6f & %6.4f w7 & %6.4f w7 &) ///
			rspec("`RSPECBOOT'")
			disp "`BOOTSTRING'"



		if "`nonull'" != ""{

			/* bootstrap confidence intervals  */

			mata: st_matrix("bootci",bootci)

			matrix rownames bootci = "CLASSIC-CV1-se" "CLASSIC-WB-se" "SCORE-CV1-se" "SCORE-WB-se"
			local RSPECBOOT "&|&|&|"

			matrix colnames bootci = "Coeff" "std.er." "WCLU CI-low" "WCLU CI-up"

			local BOOTCITITLE "Unrestricted Bootstrapped Confidence Intervals"


			matlist bootci, title(`BOOTCITITLE') rowtitle(`BOOTROW') ///
				cspec(& %-14s w14 | %9.6f w10 & %9.6f & %6.4f w15 & %6.4f w15 &) ///
				rspec("`RSPECBOOT'")

		} /*end of no null*/

	} /*end of flag boot*/

end

/*----------------------------------------*/
/*version history*/
*v1 - absorb cv3
*v1.1 - bootstraps
*v1.2 - nonull
*v1.3 - cleaned up display output
*v2 - confidence intervals
*v2.1 - correct absorb error
*v2.2 - added six point, changed var-factor issue
*v3 - corrected WCLR and more CIs
*v3.1 - changes to output and CI formula
*v3.2 - minor bug fixes
*v3.3 - added summary statistics