Work-around for Stata crash when number of fixed effects is very large: require # of FE as input, and don't represent them as linked list.

Author: droodman

boottest.ado

@@ -1,4 +1,4 @@
-*! boottest 2.1.8 3 August 2018
+*! boottest 2.1.9 4 August 2018
 *! Copyright (C) 2015-18 David Roodman
 
 * This program is free software: you can redistribute it and/or modify
@@ -33,7 +33,7 @@ program define _boottest, rclass sortpreserve
 	version 11
 
 	mata st_local("StataVersion", boottestStataVersion()); st_local("CodeVersion", boottestVersion())
-	if `StataVersion' != c(stata_version) | "`CodeVersion'" < "02.01.00" {
+	if `StataVersion' != c(stata_version) | "`CodeVersion'" < "02.01.09" {
 		cap findfile "lboottest.mlib"
 		while !_rc {
 			erase "`r(fn)'"
@@ -227,7 +227,8 @@ program define _boottest, rclass sortpreserve
 	local scoreBS = "`boottype'"=="score"
 
 	local FEname = cond(inlist("`cmd'","xtreg","xtivreg","xtivreg2"), "`e(ivar)'", "`e(absvar)'`e(absvars)'")
-
+	local NFE    = cond(inlist("`cmd'","xtreg","xtivreg","xtivreg2"),   e(N_g)   , 0`e(k_absorb)'`e(K1)')
+	
 	if `"`seed'"'!="" set seed `seed'
 
 	tempname p padj se stat df df_r hold C C0 CC0 b V keepC keepW repsname repsFeasname
@@ -322,7 +323,7 @@ program define _boottest, rclass sortpreserve
 
 		mata _boottestp = J(`cons',1,`k') \ order(tokens("`colnames'")', 1)[invorder(order(tokens("`Xnames_exog' `Xnames_endog'")', 1))] \ `k'+1
 
-		if "`FEname'"!="" & `cons' {
+		if 0`NFE' & `cons' {
 			mata _boottestp = _boottestp[|2\.|]
 			local cons 0
 			local _cons
@@ -414,7 +415,7 @@ program define _boottest, rclass sortpreserve
 		_estimates unhold `hold'
 		if r(k_autoCns) mat `C0' = `C0'[r(k_autoCns)+1...,1...]
 
-		if "`FEname'" != "" {
+		if `NFE' {
 			mata st_local("rc", strofreal(any(0:!=select(st_matrix("`C0'"),st_matrixcolstripe("`C0'")[,2]':=="_cons"))))
 
 			if `rc' {
@@ -577,7 +578,7 @@ program define _boottest, rclass sortpreserve
 		mata boottest_stata("`stat'", "`df'", "`df_r'", "`p'", "`padj'", "`cimat'", "`plotmat'", "`peakmat'", `level', `ML', `LIML', 0`fuller', `K', `ar', `null', `scoreBS', "`weighttype'", "`ptype'", ///
 												"`madjust'", `N_h0s', "`Xnames_exog'", "`Xnames_endog'", 0`cons', ///
 												"`Ynames'", "`b'", "`V'", "`W'", "`ZExclnames'", "`hold'", "`scnames'", `hasrobust', "`allclustvars'", `:word count `bootcluster'', `:word count `clustvars'', ///
-												"`FEname'", "`wtname'", "`wtype'", "`C'", "`C0'", `reps', "`repsname'", "`repsFeasname'", `small', "`svmat'", "`dist'", ///
+												"`FEname'", 0`NFE', "`wtname'", "`wtype'", "`C'", "`C0'", `reps', "`repsname'", "`repsFeasname'", `small', "`svmat'", "`dist'", ///
 												`gridmin', `gridmax', `gridpoints', `matsizegb')
 		_estimates unhold `hold'
 
@@ -693,6 +694,7 @@ program define _boottest, rclass sortpreserve
 end
 
 * Version history
+* 2.1.9 Work-around for Stata crash when number of fixed effects is very large: require # of FE as input, and don't represent them as linked list.
 * 2.1.8 Fixed 2.1.6 crash with FE. Fixed parsing of matsizegb() option.
 * 2.1.6 Changed to faster computation for WCR with many clusters, but not quite WB. In multiway only applies to intersection of all clusters.
 * 2.1.4 Fixed CI scaling issue introduced in 2.0.5 that affects scoretest, waldtest

boottest.mata

@@ -1,4 +1,4 @@
-*! boottest 2.1.8 3 August 2018
+*! boottest 2.1.9 4 August 2018
 *! Copyright (C) 2015-18 David Roodman
 
 * This program is free software: you can redistribute it and/or modify
@@ -20,7 +20,7 @@ mata set mataoptimize on
 mata set matalnum off
 
 string scalar boottestStataVersion() return("`c(stata_version)'")
-string scalar      boottestVersion() return("02.01.00")
+string scalar      boottestVersion() return("02.01.09")
 
 struct smatrix {
 	real matrix M
@@ -68,7 +68,7 @@ class boottestModel {
 	pointer (class AnalyticalModel scalar) scalar pM_Repl, pM
 	struct smatrix matrix denom
 	struct smatrix colvector Kcd, XEndstar, XExXEndstar, ZExclXEndstar, XZi, eZi, euZVR0
-	pointer (struct structFE scalar) scalar FEs
+	struct structFE rowvector FEs
 	pointer (real matrix) matrix pQ
 	pointer (struct boottest_clust scalar) scalar pBootClust
 
@@ -386,8 +386,8 @@ void boottestModel::setID      (real matrix ID, | real scalar NBootClustVar, rea
 	this.pID = &ID; this.NBootClustVar = editmissing(NBootClustVar,1); this.NErrClust=editmissing(NErrClust,1); setdirty(1)
 	if (cols(ID)) this.robust = 1
 }
-void boottestModel::setFEID      (real matrix ID) {
-	this.pFEID = &ID; setdirty(1)
+void boottestModel::setFEID(real matrix ID, real scalar NFE) {
+	this.pFEID = &ID; this.NFE = NFE; setdirty(1)
 }
 void boottestModel::setlevel  (real scalar level  )
 	this.level = level
@@ -602,7 +602,7 @@ void boottestModel::boottest() {
 	real colvector rAll, sortID, o, _FEID
 	real rowvector val, ClustCols
 	real matrix RAll, L, LAll, vec, Combs, t, IDErr
-	real scalar i, j, c, minN, sumN, _reps, g
+	real scalar i, j, c, minN, sumN, _reps, g, i_FE
 	pointer (real matrix) scalar _pR0, pIDAll
 	class AnalyticalModel scalar M_WRE
 	pragma unset vec; pragma unset val; pragma unused M_WRE
@@ -713,37 +713,36 @@ void boottestModel::boottest() {
 				J_ClustN_NBootClust = J(Clust.N, NBootClust, 0)
 		}
 
-		if (cols(*pFEID)) { // fixed effect prep
+		if (NFE) {
 			sortID = (*pFEID)[o = order(*pFEID, 1)]
-			NFE = 1; FEboot = reps>0; j = Nobs; _FEID = wtFE = J(Nobs, 1, 1)
+			i_FE = 1; FEboot = reps>0; j = Nobs; _FEID = wtFE = J(Nobs, 1, 1)
+			FEs = structFE(NFE)
 			for (i=Nobs-1;i;i--) {
 				if (sortID[i] != sortID[i+1]) {
-					NFE++
-					next = FEs; (FEs = &(structFE()))->next = next  // add new FE to linked list
-					FEs->is = o[|i+1\j|]
+					FEs[i_FE].is = o[|i+1\j|]
 					if (weights) {
-						FEs->wt  = (*pwt)[FEs->is]
-						FEs->wt = FEs->wt / colsum(FEs->wt)
+						t  = (*pwt)[FEs[i_FE].is]
+						FEs[i_FE].wt = t / colsum(t)
 					} else
-						FEs->wt = J(j-i,1,1/(j-i))
-					wtFE[FEs->is] = FEs->wt
+						FEs[i_FE].wt = J(j-i,1,1/(j-i))
+					wtFE[FEs[i_FE].is] = FEs[i_FE].wt
 					j = i
 					
 					if (FEboot & NClustVar) { // are all of this FE's obs in same bootstrapping cluster?
-						t = (*pID)[FEs->is, 1..NBootClustVar]
+						t = (*pID)[FEs[i_FE].is, 1..NBootClustVar]
 						FEboot = all(t :== t[1,])
 					}
+					++i_FE
 				}
-				_FEID[o[i]] = NFE
+				_FEID[o[i]] = i_FE
 			}
-			next = FEs; (FEs = &(structFE()))->next = next
-			FEs->is = FEs->is = o[|.\j|]
+			FEs[NFE].is = FEs[NFE].is = o[|.\j|]
 			if (weights) {
-				FEs->wt  = (*pwt)[FEs->is]
-				FEs->wt = FEs->wt / colsum(FEs->wt)
+				t  = (*pwt)[FEs[NFE].is]
+				FEs[NFE].wt = t / colsum(t)
 			} else
-				FEs->wt = J(j-i,1,1/(j-i))
-			wtFE[FEs->is] = FEs->wt
+				FEs[NFE].wt = J(j-i,1,1/(j-i))
+			wtFE[FEs[NFE].is] = FEs[NFE].wt
 			pFEID = &_FEID // ordinal fixed effect ID
 			if (robust & !(scoreBS | FEboot) & granular < NErrClustCombs)
 				infoBootAll = _panelsetup(*pIDAll, 1..NBootClustVar) // info for bootstrapping clusters wrt data collapsed to intersections of all bootstrapping & error clusters
@@ -1094,9 +1093,12 @@ real scalar boottestModel::MakeNonWRENumers(real scalar thisWeightGrpStart, real
 
 
 
+
+
+
+
 void boottestModel::MakeNonWREStats(real scalar thisWeightGrpStart, real scalar thisWeightGrpStop) {
-	real scalar d, i, c, j, l; real matrix eu, eueu, t; real colvector numer_l
-	pointer (real matrix) scalar peZVR0, pVR0
+	real scalar d, i, c, j, l; real matrix eu, eueu, t; real colvector numer_l; pointer (real matrix) scalar peZVR0, pVR0
 
 	if (robust) {
 		if (granular < NErrClustCombs & thisWeightGrpStart==1) { // if pure robust and no multi-way clustering, initialized stuff that depends on r0 but not u
@@ -1169,7 +1171,7 @@ void boottestModel::MakeNonWREStats(real scalar thisWeightGrpStart, real scalar
 					denom[i,j].M = colsum(u :* QQ * u)
 				}
 		else { // alternative core computational loop, avoiding computing Q'Q which has cubic time cost in numbers of bootstrapping clusters
-			if (granular) // prep optimized treatment when bootstrapping by small groups
+			if (granular) // prep optimized treatment when bootstrapping by many/small groups
 				if (purerobust) {
 					eu = *M_DGP.partialFE(&(pM->e :* u)) - *pX * betadev
 					eueu = eu:*eu
@@ -1332,13 +1334,12 @@ real matrix boottestModel::count_binary(real scalar N, real scalar lo, real scal
 
 // partial fixed effects out of a data matrix
 pointer(real matrix) scalar AnalyticalModel::partialFE(pointer(real matrix) scalar pIn) {
+	real matrix Out, t; real scalar i
 	if (parent->NFE & pIn!=NULL) {
-		real matrix Out, t; pointer (struct structFE scalar) scalar thisFE
-		thisFE = parent->FEs; Out = J(rows(*pIn), cols(*pIn), .)
-		while (thisFE != NULL) {
-			t = (*pIn)[thisFE->is,]
-			Out[thisFE->is,] = t :- cross(thisFE->wt, t)
-			thisFE = thisFE->next
+		Out = *pIn
+		for (i=parent->NFE;i;i--) {
+			t = Out[parent->FEs[i].is,]
+			Out[parent->FEs[i].is,] = t :- cross(parent->FEs[i].wt, t)
 		}
 		return(&Out)
 	}
@@ -1516,7 +1517,7 @@ void boottest_stata(string scalar statname, string scalar dfname, string scalar
 	real scalar K, real scalar AR, real scalar null, real scalar scoreBS, string scalar weighttype, string scalar ptype, string scalar madjtype, real scalar NumH0s,
 	string scalar XExnames, string scalar XEndnames, real scalar hascons, string scalar Ynames, string scalar bname, string scalar Vname, string scalar Wname, 
 	string scalar ZExclnames, string scalar samplename, string scalar scnames, real scalar robust, string scalar IDnames, real scalar NBootClustVar, real scalar NErrClust, 
-	string scalar FEname, string scalar wtname, string scalar wttype, string scalar Cname, string scalar C0name, real scalar reps, string scalar repsname, string scalar repsFeasname, 
+	string scalar FEname, real scalar NFE, string scalar wtname, string scalar wttype, string scalar Cname, string scalar C0name, real scalar reps, string scalar repsname, string scalar repsFeasname, 
 	real scalar small, string scalar diststat, string scalar distname, real scalar gridmin, real scalar gridmax, real scalar gridpoints, real scalar MaxMatSize) {
 
 	real matrix C, R, C0, R0, ZExcl, ID, FEID, sc, XEnd, XEx
@@ -1548,7 +1549,7 @@ void boottest_stata(string scalar statname, string scalar dfname, string scalar
 	M.setZExcl(ZExcl)
 	M.setwt (wt)
 	M.setID(ID, NBootClustVar, NErrClust)
-	M.setFEID(FEID)
+	M.setFEID(FEID, NFE)
 	M.setR (R , r )
 	M.setR0(R0, r0)
 	M.setnull(null)

crosstabs.do

@@ -0,0 +1,104 @@
+cap mata mata drop CT1()
+cap mata mata drop CT2()
+
+mata
+mata clear
+mata set matastrict on
+mata set mataoptimize on
+mata set matalnum off
+
+real matrix CT1(real colvector X, real colvector ID1, real colvector ID2) {
+	transmorphic CT, inds1, inds2, loc1, loc2; real scalar i, i1, i2, ind1, ind2; real colvector key; real matrix retval
+		
+	CT = asarray_create("real", 2); asarray_notfound(CT, 0)
+	inds1 = asarray_create("real", 1)
+	inds2 = asarray_create("real", 1)
+	
+	for (i=rows(X);i;i--) {
+		key = (ind1 = ID1[i]), (ind2  = ID2[i])
+		asarray(CT, key, asarray(CT, key) + X[i])
+		asarray(inds1, ID1[i], 1) // add to catalogs of observed instances
+		asarray(inds2, ID2[i], 1)
+	}
+
+	retval = J(asarray_elements(inds1), asarray_elements(inds2), 0)
+	i1 = 0
+	for (loc1=asarray_first(inds1); loc1!=NULL; loc1=asarray_next(inds1, loc1)) {
+		ind1 = asarray_key(inds1, loc1)
+		i1++
+		i2 = 0
+		for (loc2=asarray_first(inds2); loc2!=NULL; loc2=asarray_next(inds2, loc2)) {
+			ind2 = asarray_key(inds2, loc2)
+			retval[i1,++i2] = asarray(CT, (ind1,ind2))
+		}
+	}
+	return(retval)
+}
+
+real matrix CT2(real colvector X, real colvector ID1, real colvector ID2) {
+	class Factor scalar F
+	F = _factor((ID1,ID2))
+	return(rowshape(aggregate_sum(F, F.sort(X), ., ""), rows(uniqrows(F.keys[,1]))))
+}
+
+real matrix CT3(real colvector X, real colvector ID1, real colvector ID2) {
+	real scalar i, j, t, N, N1, N2; real matrix retval; real colvector sortID1, sortID2, sortX, _ID2, _ID2i, Xi, o1, o2; real matrix info
+	
+	N = rows(X)
+	
+	sortID1 = ID1[o1 = order(ID1, 1)]
+	sortID2 = ID2[o1]
+	sortX   = X[o1]
+
+	N1 = rows(info = _panelsetup(sortID1, 1))
+	
+	sortID2 = sortID2[o2 = order(sortID2, 1)]
+	
+	N2 = 1; _ID2 = J(N, 1, 1)
+	for (i=N-1;i;i--) {
+		if (sortID2[i] != sortID2[i+1])
+			N2++
+		_ID2[o2[i]] = N2 // ordinal version of ID2 expressed with respect to ID1-sorted data
+	}
+
+	retval = J(N1, N2, 0)
+	for (i=N1;i;i--) {
+		_ID2i = panelsubmatrix(_ID2, i, info)
+		Xi    = panelsubmatrix(sortX    , i, info)
+		for (j=rows(_ID2i);j;j--) {
+			t = _ID2i[j] 
+			retval[i,t] = retval[i,t] + Xi[j]
+		}
+	}
+	return(retval)
+}
+
+real matrix CT4(real colvector X, real colvector ID1, real colvector ID2) {
+	real colvector ID, o; real matrix info
+	ID = ID1, ID2
+	ID = ID[o = order(ID, (1,2)),]
+	info = _panelsetup(ID, (1,2))
+	return(1)
+//	return(_panelsum(X[o], _panelsetup(ID, (1,2))))
+}
+
+
+rseed(9821374)
+N = 1000000
+N1 = 5
+N2 = 5
+X = runiform(N,1)
+ID1 = rdiscrete(N,1,J(N1,1,1/N1)) :+ 5
+ID2 = rdiscrete(N,1,J(N2,1,1/N2)) :+ 10
+//CT1(X,ID1,ID2)
+CT2(X,ID1,ID2)
+CT3(X,ID1,ID2)
+CT4(X,ID1,ID2)
+
+timer_clear()
+//timer_on(1);(void) CT1(X,ID1,ID2);timer_off(1)
+timer_on(2);for(i=1;i;i--) (void) CT2(X,ID1,ID2);timer_off(2)
+timer_on(3);for(i=1;i;i--) (void) CT3(X,ID1,ID2);timer_off(3)
+timer_on(4);for(i=1;i;i--) (void) CT4(X,ID1,ID2);timer_off(4)
+timer()
+end