*! boottest 4.4.11 12 April 2024 *! Copyright (C) 2015-23 David Roodman * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * You should have received a copy of the GNU General Public License * along with this program. If not, see . mata mata clear mata set matastrict on mata set mataoptimize on mata set matalnum off struct smatrix { real matrix M } struct ssmatrix { struct smatrix matrix M } struct structboottestClust { real scalar N, multiplier, even real colvector order real matrix info } struct structFE { real colvector is, wt, sqrtwt } // return pointer to chosen columns of a matrix, but don't duplicate data if return value is whole matrix pointer (real matrix) pcol(real matrix A, real vector p) return(length(p)==cols(A)? &A : &A[,p]) // right-multiply a data matrix by a matrix, with efficient handling of special cases like latter is identity pointer (real matrix) scalar pXB(real matrix X, real matrix M) { real scalar r r = rows(M) return (all(colsum(M) :== 1) & all(colsum(!M) :== r-1)? // is M 0's but for one 1 in each col, so it's just copying/reordering cols? (all(diagonal(M)) & rows(M)==cols(M)? // is M just the identity matrix? &X : &X[,colsum(M:*(1::r))]) : // reorder X &(X * M)) } // return [X1 X2]B where X1 or X2 may have 0 cols pointer(real matrix) scalar pX12B(real matrix X1, real matrix X2, real matrix B) return(cols(B)? (cols(X1)? (cols(X2)? &(*pXB(X1, B[|.,.\cols(X1),.|]) + *pXB(X2, B[|cols(X1)+1,.\.,.|])) : pXB(X1, B)) : pXB(X2, B)) : &J(rows(X1),0,0)) // return &X[|S|] with appropiate behavior if cols(X)=0 or S requests no rows (S[2,1]=0), presuming that in degenerate cases S does not specify columns // if retval = X, doesn't duplicate data // S should be 2x1 because the function is only for selecting rows pointer(real matrix) scalar pXS(real matrix X, real colvector S) return(cols(X)? (S[2,1]? ((S[1,1]==. | S[1,1]==1) & (S[2,1]==. | S[2,1]==rows(X))? &X : &X[|S,(.\.)|]) : &J(0,cols(X),0)) : &J(editmissing(S[2,1],rows(X))-editmissing(S[1,1],1)+1,0,0)) // do X[|S|] = Y while allowing X to have no cols and S to be a colvector void setXS(real matrix X, real colvector S, real matrix Y) if (cols(X)) X[|S,(.\.)|] = Y;; matrix fold(matrix X) return(uppertriangle(X) + lowertriangle(X,0)') // fold matrix diagonally; returns same values as a quad form, but runs faster because of all the 0's class boottestOLS { // class for performing OLS real scalar y1y1, LIML, Fuller, kappa, isDGP, kZ, kX1, kX, y1bary1bar, restricted real colvector invXXXy1par, X1y1, dbetadr, beta0, y1bar, Zperpy1, t1, t1Y, deltadddot, X2y1, Xy1bar, deltaX, deltaY, u1dddot real rowvector y1Y2, Yendog, y1barU2ddot real matrix Z, ZZ, XZ, XX, PXZ, R1invR1R1, R1perp, Rpar, RperpX, RRpar, RparX, RparY, RR1invR1R1, YY, AR, XAR, XinvXX, Rt1, invXX, invH, Deltadddot, Y2bar, perpRperpX, XY2, XU2ddot, U2ddotU2ddot, Y2Y2, Piddot, ZY2, V, U2ddot pointer(real colvector) scalar py1par, pXy1par, py1 pointer(real matrix) scalar pA, pZperp, pX1, pX1par, pR1AR1, pZR1, pZperpZperp, pinvZperpZperp, pX2, pY2 pointer (class boottest scalar) scalar parent struct smatrix rowvector WXAR, CT_XAR, beta, u1ddot, XinvHg, invMg, Xg, XXg private void new(), InitTestDenoms() private virtual void InitVars(), SetR(), Estimate(), MakeResiduals() real matrix _select(), perp() } class boottestARubin extends boottestOLS { private virtual void InitVars(), Estimate() } class boottestIVGMM extends boottestOLS { real matrix H_2SLS, ZR1ZR1, X2ZR1, ZR1Y2, X1ZR1, ZZR1, H_2SLSmZZ, X2jk, Y2jk, X1jk, Zjk, ZR1jk real colvector Zy1, y1jk real rowvector twoy1ZR1 real scalar y1pary1par pointer(real rowvector) scalar pZy1par pointer(real rowvector) scalar py1parY2 pointer(real matrix) scalar pRperp pointer(struct smatrix rowvector) scalar py1parY2g, pZy1parg, pXy1parg // jk stuff struct smatrix rowvector XY2g, ZY2g, XXg, XZg, YYg, Zy1g, X1y1g, X2y1g, y1Y2g, ZZg, invXXg, H_2SLSg, H_2SLSmZZg, ZR1Y2g, ZR1ZR1g, twoy1ZR1g, ZZR1g, X1ZR1g, X2ZR1g, invHg real colvector _Nobsg, y1pary1parg, y1y1g pointer(real colvector) scalar py1pary1parg, py1parjk private void new() private virtual void InitVars(), Estimate(), MakeResiduals(), MakeH() } class boottest { real scalar scoreBS, B, small, auxwttype, null, dirty, initialized, ML, Nobs, _Nobs, kZ, kY2, kX1, sumwt, NClustVar, haswt, REst, multiplier, smallsample, quietly, FEboot, NErrClustCombs, /// sqrt, LIML, Fuller, kappa, WRE, WREnonARubin, ptype, twotailed, df, df_r, ARubin, willplot, notplotted, NumH0s, p, NBootClustVar, NErrClustVar, BootClust, FEdfadj, jk, /// NFE, granular, granularjk, purerobust, subcluster, Nstar, BFeas, v_sd, level, ptol, MaxMatSize, Nw, enumerate, bootstrapt, q, interpolable, interpolating, interpolate_u, robust, kX2, kX, Nall, bootneqcap real matrix AR, v, CI, CT_WE, infoBootData, infoErrAll, JNcapNstar, statDenom, SuwtXA, numer0, deltadenom_b, _Jcap, YYstar_b, YPXYstar_b, numerw, ustar0, YbarYbar, XYbar, invXXXZbar, PXZbar, Zbar, invXXXX1par, PiddotRparY real colvector DistCDR, plotX, plotY, beta, ClustShare, WeightGrpStart, WeightGrpStop, confpeak, gridmin, gridmax, gridpoints, numersum, anchor, poles, invFEwt, sqrtwt, Sstary1y1, CTstarFEy1 real rowvector peak, betas, As string scalar obswttype, madjtype, seed pointer (real matrix) scalar pX2, pR1, pR, pID, pFEID, pY2, pX1, pinfoAllData, pinfoCapData, pIDAll, pnumer, pustar, pU2parddot pointer (real colvector) scalar pr1, pr, py1, pSc, pwt, pA, pDist, pIDBootData, pIDBootAll class boottestOLS scalar DGP, Repl pointer(class boottestOLS scalar) scalar pM struct structboottestClust rowvector Clust struct smatrix matrix denom, Kcd, denom0, Jcd0, CTUX, FillingT0, SstarUU, negSstarUMZperpX, ScapXX1par struct smatrix rowvector Kd, dudr, dnumerdr, IDCTCapstar, infoCTCapstar, deltadenom, Zyi, SstarUPX, YYstar, YPXYstar, ScapXYbar, ScapPXYbarZperp, CT_FEcapYbar, SstarUX, SstarUXinvXX, SstarinvZperpZperpZperpU, SstarZperpU, CTstarFEU, SstarUYbar, SCTcapUXinvXX, SstarUMZperp, ScapXX, ScapYX, SstarY2X, SstarXY2par, ScapXZperp, SstarZR1X, SstarXy1, SstarZX, SstarZZ, SstarZR1Zperp, SstarZperpy1, SstarZZperp, SstarZperpY2par, SstarY2Zperp, SstarX1parY2par, SstarXX1par, SstarX1pary1, SstarY2X1par, SstarZR1X1par, SstarZX1par, SstarY2parY2par, Sstary1Y2par, Sstary1Y2, SstarZR1Y2par, SstarZR1Y2, SstarZY2par, SstarZY2, SstarY2Y2par, SstarZR1y1, SstarZy1, SstarZR1ZR1, SstarZR1Z, SstarY2Y2, SallXU2RparY, SallXu1ddot, SallU2X, CTstarFEY2, CTstarFEX, CTcapFEX, CTstarFEZR1, CTstarFEZ pointer(struct smatrix rowvector) pSstarXX, pSstarXZperp, pSallXy1, pSallZX, pSallZR1X, pSallY2X, pSallXX struct ssmatrix rowvector ddenomdr, dJcddr struct ssmatrix matrix ddenomdr2 pointer(struct smatrix matrix) scalar pJcd struct structFE rowvector FEs void new(), setsqrt(), setX1(), setptype(), setdirty(), setY2(), setY(), setX2(), setobswt(), setsc(), setML(), setLIML(), setARubin(), setauxwttype(), setFuller(), setkappa(), setquietly(), setbeta(), setA(), setsmall(), sethascons(), setjk(), setscoreBS(), setB(), setnull(), setID(), setFEID(), setlevel(), setptol(), setrobust(), setR1(), setR(), setwillplot(), setgrid(), setmadjust(), setMaxMatSize(), setstattype(), close() private void MakeNumerAndJ(), _clustAccum(), MakeWREStats(), MakeInterpolables(), _MakeInterpolables(), MakeNonWREStats(), UpdateBootstrapcDenom(), Init(), plot(), MakeWildWeights(), boottest(), crosstabCapstarMinus(), InitWRE(), PrepWRE(), storeWtGrpResults(), NoNullUpdate() real matrix getplot(), getCI(), getV(), getv() real scalar getp(), getpadj(), getstat(), getdf(), getdf_r(), getreps(), getrepsFeas(), getNBootClust() real rowvector getpeak() real colvector getdist(), getb() private real scalar r_to_p(), search() private real matrix count_binary(), crosstabFE(), HessianFixedkappa(), threebyone() private real rowvector _HessianFixedkappa() private pointer(real matrix) scalar Filling(), partialFE() private static real matrix combs() private static real colvector stableorder() private real vector _selectindex() static void _st_view() } void boottestOLS::new() { LIML = Fuller = kappa = 0 Rpar = isDGP = 1 beta = smatrix() } void boottestIVGMM::new() { Fuller = 0 kappa = isDGP = 1 u1ddot = beta = smatrix() } // do select() but handle case that both args are scalar and second is 0 by interpreting second arg as rowvector and returning J(1,0,0) // if v = 0 (so can't tell if row or col vector), returns J(1, 0, 0) real matrix boottestOLS::_select(real matrix X, real rowvector v) return (rows(X)==1 & cols(X)==1 & v==0? J(1,0,0) : select(X,v)) real matrix boottestOLS::perp(real matrix A) { real matrix vec; real rowvector val; pragma unset vec; pragma unset val symeigensystem(A*invsym(A'A)*A', vec, val); _edittozero(val, 1000) return (_select(vec, !val)) } // R1 is constraints. R is attack surface for null; only needed when using FWL for WRE // for DGP regression, R1 is maintained constraints + null if imposed while R should have 0 rows // for replication regressions R1 is maintained constraints, R is null void boottestOLS::SetR(real matrix R1, | real matrix R) { real matrix _R, vec; real rowvector val pragma unset vec; pragma unset val restricted = rows(R1) > 0 if (restricted) { R1invR1R1 = invsym(R1 * R1') if (all(diagonal(R1invR1R1))==0) _error(111, "Null hypothesis or model constraints are inconsistent or redundant.") R1invR1R1 = R1 ' R1invR1R1 symeigensystem(R1invR1R1 * R1, vec, val); _edittozero(val, 1000) R1perp = _select(vec, !val) // eigenvectors orthogonal to span of R1; foundation for parameterizing subspace compatible with constraints } else R1invR1R1 = J(parent->kZ,0,0) // and R1perp = I if (kappa) { // prepare to reduce regression via FWL _R = R \ J(parent->kY2, parent->kX1, 0), I(parent->kY2) // rows to prevent partialling out of endogenous regressors if (restricted) _R = _R * R1perp symeigensystem(_R ' invsym(_R * _R') * _R, vec, val); _edittozero(val, 1000) Rpar = _select(vec, val); if (restricted) Rpar = R1perp * Rpar // par of RR₁perp if (isDGP | !parent->WREnonARubin) { // in WRE, Rperp is same DGP and Repl; might not be in WUE, but is arranged so by call to SetR() RperpX = _select(vec, !val) if (restricted) RperpX = R1perp * RperpX RperpX = *pXS(RperpX, .\parent->kX1) } RparX = *pXS(Rpar, .\parent->kX1) // part of Rpar that refers to X1 RparY = *pXS(Rpar, parent->kX1+1\.) // part of Rpar that refers to Y2 if (isDGP==0) { RRpar = R * Rpar RR1invR1R1 = R * R1invR1R1 } } } // stuff that can be done before r set, and depends only on exogenous variables, which are fixed throughout all bootstrap methods void boottestOLS::InitVars(pointer(real matrix) scalar pRperp) { // pRperp is for replication regression--no null imposed real matrix H, negR1AR1; real scalar g; real colvector S u1ddot = smatrix(1 + parent->jk) // for jackknife, need jk'd residuals but also non-jk'd residuals for original test stat py1par = parent->py1 H = cross(*parent->pX1, *parent->pX1) invH = invsym(H) pR1AR1 = rows(R1perp)? &( R1perp * invsym( R1perp ' H * R1perp) * R1perp') : &invH // for DGP regression X1y1 = cross(*parent->pX1, *py1par) beta0 = *pR1AR1 * X1y1 dbetadr = *pR1AR1 * H * R1invR1R1 - R1invR1R1 if (parent->jk) if (parent->purerobust) // set up for optimized "robust" jk (invMg = smatrix()).M = 1 :/ (1 :- rowsum(*parent->pX1 * fold(*pR1AR1) :* *parent->pX1)) // standard hc3 multipliers else { u1ddot[2].M = J(parent->Nobs, 1, 0) Xg = smatrix(parent->Nstar) if (parent->granularjk) invMg = Xg else XXg = XinvHg = Xg negR1AR1 = - *pR1AR1 for (g=parent->Nstar; g; g--) { S = parent->NClustVar? parent->infoBootData[g,1] \ parent->infoBootData[g,2] : g\g Xg[g].M = *pXS(*parent->pX1,S) if (parent->granularjk) { // for many small clusters, faster to compute jackknife errors vaia hc3-like formula invMg[g].M = Xg[g].M * negR1AR1 * Xg[g].M' _diag(invMg[g].M, diagonal(invMg[g].M) :+ 1) invMg[g].M = invsym(invMg[g].M) } else { XXg[g].M = cross(Xg[g].M, Xg[g].M) XinvHg[g].M = Xg[g].M * (rows(R1perp)? R1perp * invsym(R1perp ' (H - XXg[g].M) * R1perp) * R1perp' : invsym(H - XXg[g].M)) } } } pA = rows(*pRperp)? &(*pRperp * invsym(*pRperp ' H * *pRperp) * *pRperp') : &invH // for replication regression AR = *pA * *parent->pR' if (parent->scoreBS | parent->robust) XAR = *parent->pX1 * AR } void boottestARubin::InitVars(| pointer(real matrix) pRperp) { pragma unused pRperp real matrix H, X2X1, _X2X1, _X1, _X2, tmp; real colvector S, X1y1, X2y1; real scalar g u1ddot = smatrix(1 + parent->jk) // for jackknife, need jk'd residuals but also non-jk'd residuals for original test stat X2X1 = cross(*parent->pX2, *parent->pX1) H = cross(*parent->pX1, *parent->pX1), X2X1' \ X2X1, cross(*parent->pX2, *parent->pX2) pA = &invsym(H) pR1AR1 = rows(R1perp)? &(R1perp * invsym(R1perp ' H * R1perp) * R1perp') : pA X1y1 = cross(*parent->pX1, *parent->py1) X2y1 = cross(*parent->pX2, *parent->py1) beta0 = *pR1AR1 * (X1y1 \ X2y1) dbetadr = *pR1AR1 * (cross(*parent->pX1, *parent->pY2) \ cross(*parent->pX2, *parent->pY2)) if (parent->jk) if (parent->purerobust) { // set up for optimized "robust" jk tmp = fold(*pR1AR1) (invMg=smatrix()).M = parent->kX1? rowsum(*parent->pX1 * tmp[|.,.\parent->kX1,parent->kX1|] :* *parent->pX1) + rowsum(*parent->pX1 * tmp[|.,parent->kX1+1\parent->kX1,.|] :* *parent->pX2) + rowsum(*parent->pX2 * tmp[|parent->kX1+1,parent->kX1+1\.,.|] :* *parent->pX2) : rowsum(*parent->pX2 * tmp :* *parent->pX2) invMg.M = 1 :/ (1 :- invMg.M) // standard hc3 multipliers } else { u1ddot[2].M = J(parent->Nobs, 1, 0) Xg = smatrix(parent->Nstar) if (parent->granularjk) invMg = Xg else XXg = XinvHg = smatrix(parent->Nstar) for (g=parent->Nstar; g; g--) { S = parent->NClustVar? parent->infoBootData[g,1] \ parent->infoBootData[g,2] : g\g _X1 = *pXS(*parent->pX1, S) _X2 = *pXS(*parent->pX2, S) Xg[g].M = _X1, _X2 if (parent->granularjk) { // for many small clusters, faster to compute jackknife errors vaia hc3-like formula invMg[g].M = -(Xg[g].M * *pR1AR1 * Xg[g].M') _diag(invMg[g].M, diagonal(invMg[g].M) :+ 1) invMg[g].M = invsym(invMg[g].M) } else { _X2X1 = cross(_X2, _X1) XXg[g].M = cross(_X1, _X1), _X2X1' \ _X2X1, cross(_X2, _X2) XinvHg[g].M = Xg[g].M * (rows(R1perp)? R1perp * invsym(R1perp ' (H - XXg[g].M) * R1perp) * R1perp' : invsym(H - XXg[g].M)) } } } AR = *pA * *parent->pR' // for replication regression if (parent->scoreBS | parent->robust) XAR = *pX12B(*parent->pX1, *parent->pX2, AR) } void boottestIVGMM::InitVars(|pointer(real matrix) scalar pRperp) { real matrix ZperpZ, ZperpZR1, _invZperpZperp, X1g, X2g, Zg, Zperpg, ZR1g, Y2g, _X2X1, _X1X1, _X2X2, _X1Y2, _X2Y2, X2X1, _invZperpZperpZperpX1, _invZperpZperpZperpX2, _invZperpZperpZperpZ, _invZperpZperpZperpZR1, _invZperpZperpZperpY2, ZperpX1, ZperpX2, ZperpY2, X1Y2, X2X2, X1X1, X2Y2, X1Z, X2Z, tX1, tX2, tY2, tZ, tZR1 real colvector S, y1g, _invZperpZperpZperpy1, Zperpy1, ty1 real rowvector y1ZR1 real scalar g this.pRperp = pRperp if (!isDGP & parent->WREnonARubin) { pZperp = parent->DGP.pZperp pinvZperpZperp = parent->DGP.pinvZperpZperp } else { perpRperpX = perp(RperpX) pZperp = pXB(*parent->pX1, RperpX) pinvZperpZperp = &invsym(*(pZperpZperp = &cross(*pZperp, *pZperp))) pX1 = pXB(*parent->pX1, perpRperpX) // retained exogenous regressors in instrument set ZperpX1 = cross(*pZperp, *pX1) ZperpX2 = cross(*pZperp, *parent->pX2) ZperpY2 = cross(*pZperp, *parent->pY2) Zperpy1 = cross(*pZperp, *parent->py1) kX = (kX1 = cols(*pX1)) + parent->kX2 u1ddot = smatrix() } kZ = cols(Rpar) Z = *(pX1par = pXB(*parent->pX1, RparX)) + *parent->pY2 * RparY // Zpar, X1par (retained exogenous regressors in regressor set) pZR1 = pX12B(*parent->pX1, *parent->pY2, R1invR1R1) // Z*R1 ZperpZ = cross(*pZperp, Z ) ZperpZR1 = cross(*pZperp, *pZR1) if (parent->jk & isDGP & parent->WREnonARubin) { XY2g = ZY2g = XXg = XZg = YYg = Zy1g = X1y1g = X2y1g = y1Y2g = invHg = ZZg = invXXg = H_2SLSg = H_2SLSmZZg = ZR1Y2g = ZR1ZR1g = twoy1ZR1g = ZZR1g = X1ZR1g = X2ZR1g = smatrix(parent->Nstar) y1y1g = J(parent->Nstar, 1, 0) beta = smatrix(parent->Nstar + 1) u1ddot.M = u1dddot = J(parent->Nobs, 1, 0); U2ddot = J(parent->Nobs, parent->kY2, 0) if (restricted) { py1pary1parg = &J(parent->Nstar,1,0) py1parY2g = &smatrix(parent->Nstar) pZy1parg = &smatrix(parent->Nstar) pXy1parg = &smatrix(parent->Nstar) } Y2jk = J(parent->Nobs, parent->kY2, 0) X2jk = J(parent->Nobs, parent->kX2, 0) X1jk = J(parent->Nobs, kX1, 0) Zjk = J(parent->Nobs, kZ, 0) ZR1jk = J(parent->Nobs, cols(*pZR1), 0) y1jk = J(parent->Nobs, 1, 0) X2X1 = cross(*parent->pX2, *pX1) X1X1 = cross(*pX1, *pX1) X2X2 = cross(*parent->pX2, *parent->pX2) X2Y2 = cross(*parent->pX2, *parent->pY2) X2y1 = cross(*parent->pX2, *parent->py1) X1y1 = cross(*pX1, *parent->py1) Zy1 = cross(Z, *parent->py1) X1Y2 = cross(*pX1, *parent->pY2) X1Z = cross(*pX1, Z) X2Z = cross(*parent->pX2, Z) ZZ = cross(Z, Z) y1Y2 = cross(*parent->py1, *parent->pY2) ZY2 = cross(Z, *parent->pY2) y1y1 = cross(*parent->py1, *parent->py1) X1ZR1 = cross(*pX1, *pZR1) X2ZR1 = cross(*parent->pX2, *pZR1) ZZR1 = cross(Z, *pZR1) y1ZR1 = cross(*parent->py1, *pZR1) ZR1ZR1 = cross(*pZR1, *pZR1) ZR1Y2 = cross(*pZR1, *parent->pY2) for (g=parent->Nstar; g; g--) { S = parent->NClustVar? parent->infoBootData[g,]' : g\g X2g = *pXS(*parent->pX2, S) Y2g = *pXS(*parent->pY2, S) y1g = *pXS(*parent->py1, S) Zperpg = *pXS(*pZperp, S) X1g = *pXS(*pX1, S) Zg = *pXS(Z , S) ZR1g = *pXS(*pZR1, S) _invZperpZperp = invsym(*pZperpZperp - cross(Zperpg, Zperpg)) _invZperpZperpZperpX1 = _invZperpZperp * (ZperpX1 - cross(Zperpg, X1g)) _invZperpZperpZperpX2 = _invZperpZperp * (ZperpX2 - cross(Zperpg, X2g)) _invZperpZperpZperpZ = _invZperpZperp * (ZperpZ - cross(Zperpg, Zg)) _invZperpZperpZperpZR1 = _invZperpZperp * (ZperpZR1 - cross(Zperpg, ZR1g)) _invZperpZperpZperpY2 = _invZperpZperp * (ZperpY2 - cross(Zperpg, Y2g)) _invZperpZperpZperpy1 = _invZperpZperp * (Zperpy1 - cross(Zperpg, y1g)) X1g = X1g - Zperpg * _invZperpZperpZperpX1 // FWL-process X2g = X2g - Zperpg * _invZperpZperpZperpX2 Zg = Zg - Zperpg * _invZperpZperpZperpZ ZR1g = ZR1g - Zperpg * _invZperpZperpZperpZR1 Y2g = Y2g - Zperpg * _invZperpZperpZperpY2 y1g = y1g - Zperpg * _invZperpZperpZperpy1 setXS(X1jk , S, X1g ) // save partialled-out vars from each jk iteration, to compute residuals later setXS(X2jk , S, X2g ) setXS(Y2jk , S, Y2g ) setXS(y1jk , S, y1g ) setXS(Zjk , S, Zg ) setXS(ZR1jk, S, ZR1g) tX1 = ZperpX1 - (*pZperpZperp) * _invZperpZperpZperpX1 tX2 = ZperpX2 - (*pZperpZperp) * _invZperpZperpZperpX2 tZ = ZperpZ - (*pZperpZperp) * _invZperpZperpZperpZ tY2 = ZperpY2 - (*pZperpZperp) * _invZperpZperpZperpY2 ty1 = Zperpy1 - (*pZperpZperp) * _invZperpZperpZperpy1 _X2X1 = X2X1 - ZperpX2 ' _invZperpZperpZperpX1 - _invZperpZperpZperpX2 ' tX1 - cross(X2g, X1g) _X1X1 = X1X1 - ZperpX1 ' _invZperpZperpZperpX1 - _invZperpZperpZperpX1 ' tX1 - cross(X1g, X1g) _X2X2 = X2X2 - ZperpX2 ' _invZperpZperpZperpX2 - _invZperpZperpZperpX2 ' tX2 - cross(X2g, X2g) _X1Y2 = X1Y2 - ZperpX1 ' _invZperpZperpZperpY2 - _invZperpZperpZperpX1 ' tY2 - cross(X1g, Y2g) _X2Y2 = X2Y2 - ZperpX2 ' _invZperpZperpZperpY2 - _invZperpZperpZperpX2 ' tY2 - cross(X2g, Y2g) y1Y2g[g].M = y1Y2 - Zperpy1 ' _invZperpZperpZperpY2 - _invZperpZperpZperpy1 ' tY2 - cross(y1g, Y2g) X2y1g[g].M = X2y1 - ZperpX2 ' _invZperpZperpZperpy1 - _invZperpZperpZperpX2 ' ty1 - cross(X2g, y1g) X1y1g[g].M = X1y1 - ZperpX1 ' _invZperpZperpZperpy1 - _invZperpZperpZperpX1 ' ty1 - cross(X1g, y1g) y1y1g[g] = y1y1 - 2 * (Zperpy1 ' _invZperpZperpZperpy1) + _invZperpZperpZperpy1 ' (*pZperpZperp) * _invZperpZperpZperpy1 - cross(y1g, y1g) Zy1g [g].M = Zy1 - ZperpZ ' _invZperpZperpZperpy1 - _invZperpZperpZperpZ ' ty1 - cross(Zg , y1g) XZg [g].M = X1Z - ZperpX1 ' _invZperpZperpZperpZ - _invZperpZperpZperpX1 ' tZ - cross(X1g, Zg) \ X2Z - ZperpX2 ' _invZperpZperpZperpZ - _invZperpZperpZperpX2 ' tZ - cross(X2g, Zg) ZZg[g].M = ZZ - ZperpZ ' _invZperpZperpZperpZ - _invZperpZperpZperpZ ' tZ - cross(Zg , Zg) ZY2g[g].M = ZY2 - ZperpZ ' _invZperpZperpZperpY2 - _invZperpZperpZperpZ ' tY2 - cross(Zg, Y2g) XY2g[g].M = _X1Y2 \ _X2Y2 invXXg[g].M = invsym((XXg[g].M = _X1X1, _X2X1' \ _X2X1, _X2X2)) H_2SLSg[g].M = XZg[g].M ' invXXg[g].M * XZg[g].M if (kappa!=1 | LIML) H_2SLSmZZg[g].M = H_2SLSg[g].M - ZZg[g].M if (restricted) { tZR1 = ZperpZR1 - (*pZperpZperp) * _invZperpZperpZperpZR1 X1ZR1g[g].M = X1ZR1 - ZperpX1 ' _invZperpZperpZperpZR1 - _invZperpZperpZperpX1 ' tZR1 - cross(X1g , ZR1g) X2ZR1g[g].M = X2ZR1 - ZperpX2 ' _invZperpZperpZperpZR1 - _invZperpZperpZperpX2 ' tZR1 - cross(X2g , ZR1g) ZZR1g[g].M = ZZR1 - ZperpZ ' _invZperpZperpZperpZR1 - _invZperpZperpZperpZ ' tZR1 - cross(Zg , ZR1g) twoy1ZR1g[g].M = y1ZR1 - Zperpy1 ' _invZperpZperpZperpZR1 - _invZperpZperpZperpy1 ' tZR1 - cross(y1g , ZR1g) twoy1ZR1g[g].M = twoy1ZR1g[g].M + twoy1ZR1g[g].M ZR1ZR1g[g].M = ZR1ZR1 - ZperpZR1 ' _invZperpZperpZperpZR1 - _invZperpZperpZperpZR1 ' tZR1 - cross(ZR1g, ZR1g) ZR1Y2g[g].M = ZR1Y2 - ZperpZR1 ' _invZperpZperpZperpY2 - _invZperpZperpZperpZR1 ' tY2 - cross(ZR1g, Y2g ) } } if (!restricted) { py1parY2g = &y1Y2g pZy1parg = &Zy1g py1pary1parg = &y1y1g pXy1parg = &smatrix(parent->Nstar); for (g=parent->Nstar;g;g--) (*pXy1parg)[g].M = X1y1g[g].M \ X2y1g[g].M py1parjk = &y1jk } if (Fuller) // number of observations outside each bootstrapping cluster _Nobsg = parent->_Nobs :- (parent->obswttype=="fweight"? panelsum(*parent->pwt, parent->infoBootData) : (parent->infoBootData[,2] - parent->infoBootData[,1]) :+ 1) } if (!isDGP & parent->WREnonARubin) { pX1 = parent->DGP.pX1 pX2 = parent->DGP.pX2 invXX = parent->DGP.invXX pY2 = parent->DGP.pY2 py1 = parent->DGP.py1 y1Y2 = parent->DGP.y1Y2 X2y1 = parent->DGP.X2y1 X1y1 = parent->DGP.X1y1 y1y1 = parent->DGP.y1y1 } else { pX1 = &(*pX1 - *pZperp * (*pinvZperpZperp * ZperpX1)) // FWL-process X1 pX2 = &(*parent->pX2 - *pZperp * (*pinvZperpZperp * ZperpX2)) pY2 = &(*parent->pY2 - *pZperp * (*pinvZperpZperp * ZperpY2)) py1 = &(*parent->py1 - *pZperp * (*pinvZperpZperp * Zperpy1)) X2X1 = cross(*pX2, *pX1) invXX = invsym(XX = (cross(*pX1, *pX1), X2X1' \ X2X1, cross(*pX2, *pX2))) if (parent->scoreBS==0) Y2Y2 = cross(*pY2, *pY2) XY2 = cross(*pX1, *pY2) \ cross(*pX2, *pY2) y1Y2 = cross(*py1, *pY2) X2y1 = cross(*pX2, *py1) X1y1 = cross(*pX1, *py1) y1y1 = cross(*py1, *py1) } Z = Z - *pZperp * (*pinvZperpZperp * ZperpZ) // XXX when is this avoidable, like in jk? Zy1 = cross(Z, *py1) XZ = cross(*pX1, Z) \ cross(*pX2, Z) V = invXX * XZ ZZ = cross(Z,Z) if (parent->WREnonARubin) ZY2 = cross(Z, *pY2) if (restricted) { pZR1 = &(*pZR1 - *pZperp * (*pinvZperpZperp * ZperpZR1)) // XXX when is this avoidable, like in jk? X2ZR1 = cross(*pX2 , *pZR1) X1ZR1 = cross(*pX1 , *pZR1) ZZR1 = cross(Z , *pZR1) twoy1ZR1 = cross( *py1 , *pZR1); twoy1ZR1 = twoy1ZR1 + twoy1ZR1 ZR1ZR1 = cross(*pZR1, *pZR1) ZR1Y2 = cross(*pZR1, *pY2 ) } else { py1parY2 = &y1Y2 pZy1par = &Zy1 y1pary1par = y1y1 pXy1par = &(X1y1 \ X2y1) py1par = py1 t1 = J(parent->kZ,1,0) t1Y = J(parent->kY2,1,0) } if (isDGP) { if (parent->jk & LIML==0 & parent->WREnonARubin) for (g=parent->Nstar; g; g--) invHg[g].M = invsym(kappa==1? H_2SLSg[g].M : ZZg[g].M + kappa * H_2SLSmZZg[g].M) H_2SLS = V ' XZ // Hessian if (kappa!=1 | LIML) H_2SLSmZZ = H_2SLS - ZZ if (LIML==0) // DGP is LIML except possibly when getting confidence peak for A-R plot; but LIML=0 when exactly id'd, for then kappa=1 always and Hessian doesn't depend on r1 and can be computed now MakeH() } else { pX1par = &(*pX1par - *pZperp * (*pinvZperpZperp * cross(*pZperp, *pX1par))) Yendog = 1, colsum(RparY :!= 0) // columns of Y = [y1par Zpar] that are endogenous (normally all) if (parent->robust & parent->bootstrapt & (parent->granular | parent->jk)) { // for WRE replication regression, prepare for CRVE XinvXX = *pX12B(*pX1, *pX2, invXX) if ((parent->granular | parent->jk) & parent->WREnonARubin) PXZ = *pX12B(*pX1, *pX2, V) } } } // do most of estimation; for LIML r1 must be passed now in order to solve eigenvalue problem involving it // inconsistency: for replication regression of Anderson-Rubin, r1 refers to the *null*, not the maintained constraints, because that's what affects the endogenous variables // For WRE, should only be called once for the replication regressions, since for them r1 is the unchanging model constraint void boottestOLS::Estimate(real scalar _jk, real colvector r1) { beta.M = beta0 - dbetadr * r1 if (_jk & rows(R1perp)) t1 = R1invR1R1 * r1 } void boottestARubin::Estimate(real scalar _jk, real colvector r1) { beta.M = beta0 - dbetadr * r1 py1par = &(*parent->py1 - *parent->pY2 * r1) if (_jk & rows(R1perp)) t1 = R1invR1R1 * r1 } void boottestIVGMM::MakeH() { pointer(real matrix) scalar pH pH = kappa==1? &H_2SLS : &(ZZ + kappa * H_2SLSmZZ) invH = invsym(*pH) if (pRperp) { // for score bootstrap pA = cols(*pRperp)? &(*pRperp * invsym(*pRperp ' (*pH) * *pRperp) * *pRperp') : &invH AR = *pA * (Rpar ' (*parent->pR')) XAR = *pX12B(*pX1, *pX2, V * AR) } } void boottestIVGMM::Estimate(real scalar _jk, real colvector r1) { real rowvector val; real matrix vec; real scalar g, kappag; real colvector ZXinvXXXy1par, invXXXy1parg; real matrix YPXY pragma unset vec; pragma unset val if (restricted) { t1 = R1invR1R1 * r1 if (isDGP) py1par = &(*py1 - *pZR1 * r1) y1pary1par = y1y1 - twoy1ZR1 * r1 + r1 ' ZR1ZR1 * r1 py1parY2 = &(y1Y2 - r1 ' ZR1Y2) pZy1par = &( Zy1 - ZZR1 * r1) pXy1par = &(X1y1 - X1ZR1 * r1 \ X2y1 - X2ZR1 * r1) if (_jk) { t1Y = t1[|parent->kX1+1,.\.,.|] for (g=parent->Nstar; g; g--) { py1parjk = &(y1jk - ZR1jk * r1) (*py1pary1parg)[g] = y1y1g[g] - twoy1ZR1g[g].M * r1 + r1 ' ZR1ZR1g[g].M * r1 (*py1parY2g) [g].M = y1Y2g[g].M - r1 ' ZR1Y2g[g].M (*pZy1parg) [g].M = Zy1g[g].M - ZZR1g[g].M * r1 (*pXy1parg) [g].M = X1y1g[g].M - X1ZR1g[g].M * r1 \ X2y1g[g].M - X2ZR1g[g].M * r1 } } } invXXXy1par = invXX * *pXy1par YY = y1pary1par, *pZy1par' \ *pZy1par, ZZ if (isDGP) { ZXinvXXXy1par = XZ ' invXXXy1par if (LIML) { YPXY = invXXXy1par ' (*pXy1par) , ZXinvXXXy1par' \ ZXinvXXXy1par , H_2SLS eigensystemselecti(invsym(YY) * YPXY, rows(YY)\rows(YY), vec, val) kappa = 1/(1 - Re(val)) // sometimes a tiny imaginary component sneaks into val if (Fuller) kappa = kappa - Fuller / (parent->_Nobs - parent->kX) MakeH() } beta.M = invH * (kappa==1? ZXinvXXXy1par : kappa * (ZXinvXXXy1par - *pZy1par) + *pZy1par) if (_jk) for (g=parent->Nstar; g; g--) { ZXinvXXXy1par = XZg[g].M ' (invXXXy1parg = invXXg[g].M * (*pXy1parg)[g].M) YYg[g].M = (*py1pary1parg)[g], (*pZy1parg)[g].M' \ (*pZy1parg)[g].M, ZZg[g].M if (LIML) { YPXY = invXXXy1parg ' (*pXy1parg)[g].M , ZXinvXXXy1par' \ ZXinvXXXy1par , H_2SLSg[g].M eigensystemselecti(invsym(YYg[g].M) * YPXY, rows(YYg[g].M)\rows(YYg[g].M), vec, val) kappag = 1/(1 - Re(val)) if (Fuller) kappag = kappag - Fuller / (_Nobsg[g] - parent->kX) invHg[g].M = invsym(ZZg[g].M + kappag * H_2SLSmZZg[g].M) } else kappag = kappa beta[g+1].M = invHg[g].M * (kappag==1? ZXinvXXXy1par : kappag * (ZXinvXXXy1par - (*pZy1parg)[g].M) + (*pZy1parg)[g].M) } } else if (parent->WREnonARubin) // if not score bootstrap of IV/GMM... Rt1 = RR1invR1R1 * r1 } void boottestOLS::MakeResiduals(real scalar _jk) { real scalar g, m; real colvector S, u1ddotg, Xt1 u1ddot.M = *py1par - *pX12B(*parent->pX1, *parent->pX2, beta.M) if (_jk) { m = parent->small? sqrt((parent->Nstar - 1) / parent->Nstar) : 1 if (parent->purerobust) // somewhat faster handling classic hc3 if (rows(R1perp)) { Xt1 = *pX12B(*parent->pX1, *parent->pX2, t1) u1ddot[2].M = m * ((invMg.M :* (u1ddot.M + Xt1)) - Xt1) } else u1ddot[2].M = m * invMg.M :* u1ddot.M else if (parent->granularjk) { for (g=parent->Nstar; g; g--) { S = parent->NClustVar? parent->infoBootData[g,]' : g\g if (rows(R1perp)) { Xt1 = Xg[g].M * t1 u1ddot[2].M[|S|] = m * (cross(invMg[g].M, u1ddot.M[|S|] + Xt1) - Xt1) } else u1ddot[2].M[|S|] = m * cross(invMg[g].M, u1ddot.M[|S|]) } } else for (g=parent->Nstar; g; g--) { S = parent->infoBootData[g,]' u1ddotg = u1ddot.M[|S|] u1ddot[2].M[|S|] = m * (u1ddotg + XinvHg[g].M * (rows(R1perp)? cross(Xg[g].M, u1ddotg) + XXg[g].M * t1 : cross(Xg[g].M, u1ddotg))) } } } void boottestIVGMM::MakeResiduals(real scalar _jk) { pragma unused _jk real matrix tmp, Piddotjk; real colvector Xu, negXuinvuu, _beta, delta, S, _t1Y; real scalar uu; real scalar g if (parent->jk) for (g=parent->Nstar; g; g--) { S = parent->NClustVar? parent->infoBootData[g,1] \ parent->infoBootData[g,2] : g\g u1ddot.M[|S|] = *pXS(*py1parjk,S) - *pXS(Zjk,S) * beta[1+g].M } else if (parent->granular | parent->scoreBS) u1ddot.M = *py1par - Z * beta.M if (parent->scoreBS==0) { _beta = 1 \ -beta.M uu = _beta ' YY * _beta Xu = *pXy1par - XZ * beta.M // after DGP regression, compute Y2 residuals by regressing Y2 on X while controlling for y1 residuals, done through FWL negXuinvuu = Xu / -uu Piddot = invsym(XX + negXuinvuu * Xu') * (negXuinvuu * (*py1parY2 - beta.M ' ZY2) + XY2) U2ddotU2ddot = Y2Y2 - Piddot'XY2 - XY2'Piddot + Piddot'XX*Piddot delta = Rpar * beta.M + t1 - parent->Repl.t1 deltaX = *pXS(delta, .\parent->kX1) deltaY = *pXS(delta, parent->kX1+1\.) tmp = perpRperpX'deltaX deltadddot = Piddot * deltaY; if (kX1) deltadddot[|.\kX1|] = deltadddot[|.\kX1|] + tmp // (X_∥'X_∥)^-1 * X_∥'y1bar Xy1bar = XX * deltadddot y1bary1bar = deltadddot'Xy1bar XU2ddot = XY2 - XX * Piddot y1barU2ddot = deltadddot'XU2ddot if (parent->granular | parent->jk) { Y2bar = *pX12B(*pX1, *pX2, Piddot) y1bar = *pX1 * tmp + Y2bar * deltaY } if (parent->jk) { _t1Y = t1Y - parent->Repl.t1Y for (g=parent->Nstar; g; g--) { S = parent->NClustVar? parent->infoBootData[g,1] \ parent->infoBootData[g,2] : g\g _beta = 1 \ -beta[1+g].M uu = _beta ' YYg[g].M * _beta Xu = (*pXy1parg)[g].M - XZg[g].M * beta[1+g].M negXuinvuu = Xu / -uu Piddotjk = invsym(XXg[g].M + negXuinvuu * Xu') * (negXuinvuu * ((*py1parY2g)[g].M - beta[1+g].M ' ZY2g[g].M) + XY2g[g].M) U2ddot [|S,(.\.)|] = tmp = *pXS(Y2jk,S) - *pX12B(*pXS(X1jk,S), *pXS(X2jk,S), Piddotjk) u1dddot[|S |] = u1ddot.M[|S|] + tmp * (RparY * beta[1+g].M + _t1Y) } } else if (parent->granular) { U2ddot = *pY2 - Y2bar u1dddot = u1ddot.M + U2ddot * deltaY } } } // non-WRE stuff that only depends on r in A-R case, for test stat denominators in replication regressions // since the non-AR OLS code never creates an object for replication regresssions, in that case this is called on the DGP regression object // depends on results of Estimate() only when doing OLS-style bootstrap on an overidentified IV/GMM regression--score bootstrap or A-R. Then kappa from DGP LIML affects Hessian, pH. void boottestOLS::InitTestDenoms() { real scalar d if (parent->bootstrapt & (parent->scoreBS | parent->robust)) { if (parent->granular | parent->purerobust) { WXAR = smatrix(parent->df) for (d=parent->df;d;d--) WXAR[d].M = XAR[,d] } if (parent->NFE & parent->robust & (parent->FEboot | parent->scoreBS)==0 & parent->granular < parent->NErrClustCombs) { // make first factor of second term of (64) for c=∩ (c=1) CT_XAR = smatrix(parent->df) for (d=parent->df;d;d--) CT_XAR[d].M = parent->crosstabFE(XAR[,d], *parent->pinfoCapData) } } } // partial fixed effects out of a data matrix pointer(real matrix) scalar boottest::partialFE(pointer(real matrix) scalar pIn) { real matrix Out, tmp; real scalar i if (NFE & pIn) { Out = *pIn if (haswt) for (i=NFE;i;i--) { tmp = Out[FEs[i].is,] Out[FEs[i].is,] = tmp :- FEs[i].sqrtwt * cross(FEs[i].wt, tmp) } else for (i=NFE;i;i--) { tmp = Out[FEs[i].is,] Out[FEs[i].is,] = tmp :- cross(FEs[i].wt, tmp) } return(&Out) } return (pIn) } void boottest::new() { ARubin = LIML = Fuller = WRE = small = scoreBS = auxwttype = ML = initialized = quietly = sqrt = ptype = robust = NFE = FEboot = granular = NErrClustCombs = subcluster = B = BFeas = interpolating = jk = 0 twotailed = null = dirty = willplot = v_sd = notplotted = FEdfadj = bootstrapt = 1 level = 95 ptol = 1e-6 confpeak = MaxMatSize = . pY2 = pX1 = pX2 = py1 = pSc = pID = pFEID = pR1 = pR = pwt = &J(0,0,0) pr1 = pr = &J(0,1,0) pIDBootData = pIDBootAll = &. } // break loops in data structure topology to enable garbage collection void boottest::close() { DGP.parent = Repl.parent = NULL } void boottest::setdirty(real scalar _dirty, | real scalar noinitialize) { dirty = _dirty if (_dirty & noinitialize!=1) initialized = 0 } void boottest::setsqrt(real scalar _sqrt) { if (_sqrt < sqrt) { if (dirty==0) { pDist = &(*pDist :* *pDist) multiplier = multiplier * multiplier } } else setdirty(1) sqrt = _sqrt } void boottest::setptype(string scalar ptype) { real scalar p p = cross( (strtrim(strlower(ptype)) :== ("symmetric"\"equaltail"\"lower"\"upper")), 1::4 ) - 1 if (p<0) _error(198, "p-value type must be "symmetric", "equaltail", "lower", or "upper"."') this.ptype = p this.twotailed = p<=1 } void boottest::setstattype(string scalar stattype) { real scalar p p = cross( (strtrim(strlower(stattype)) :== ("c"\"t")), 1::2 ) - 1 if (p<0) _error(198, "statistic type must be "t" or "c"."') this.bootstrapt = p setdirty(1) } void boottest::setX1(real matrix X1) { this.pX1 = &X1; setdirty(1) } void boottest::setX2(real matrix X2) { this.pX2 = &X2; setdirty(1) } void boottest::setY(real matrix y1) { this.py1 = &y1; setdirty(1) } void boottest::setY2(real matrix Y2) { this.pY2 = &Y2; setdirty(1) } void boottest::setobswt(real matrix wt, string scalar obswttype) { this.pwt = &wt; this.obswttype = obswttype; setdirty(1) } void boottest::setsc(real matrix Sc) { this.pSc = &Sc setdirty(1) } void boottest::setML(real scalar ML) { this.ML = ML; setdirty(1) if (ML) setscoreBS(1) } void boottest::setLIML(real scalar LIML) { this.LIML = LIML; setdirty(1) } void boottest::setARubin(real scalar ARubin) { this.ARubin = ARubin; setdirty(1) } void boottest::setFuller(real scalar Fuller) { this.Fuller = Fuller; setdirty(1) } void boottest::setkappa(real scalar kappa) { // kappa as in k-class this.kappa = kappa; setdirty(1) } void boottest::setquietly(real scalar quietly ) this.quietly = quietly void boottest::setbeta(real colvector beta) { this.beta = beta; setdirty(1) } void boottest::setA(real matrix A) { this.pA = &A; setdirty(1) } void boottest::setsmall(real scalar small) { this.small = small; setdirty(1) } void boottest::setscoreBS (real scalar scoreBS) { this.scoreBS = scoreBS; setdirty(1) } void boottest::setjk (real scalar jk) { this.jk = jk; setdirty(1) } void boottest::setB(real scalar B) { this.B = B if (B==0) setscoreBS(1) setdirty(1) } void boottest::setnull (real scalar null) { this.null = null; setdirty(1) } void boottest::setID (real matrix ID, | real scalar NBootClustVar, real scalar NErrClustVar) { this.pID = &ID; this.NBootClustVar = editmissing(NBootClustVar,1); this.NErrClustVar=editmissing(NErrClustVar,editmissing(NBootClustVar,1)); setdirty(1) if (cols(ID)) this.robust = 1 } void boottest::setFEID(real matrix ID, real scalar NFE, | real scalar FEdfadj) { this.pFEID = &ID; this.NFE = NFE; this.FEdfadj = editmissing(FEdfadj,1); setdirty(1) } void boottest::setlevel(real scalar level) this.level = level void boottest::setptol(real scalar ptol) this.ptol = ptol void boottest::setrobust(real scalar robust) { this.robust = robust if (robust==0) setID(J(0,0,0), 1, 1) setdirty(1) } void boottest::setR1(real matrix R1, real matrix r1) { this.pR1 = &R1; this.pr1 = &r1; setdirty(1) } void boottest::setR(real matrix R, real colvector r) { this.pR = &R; this.pr = &r; q = rows(R); setdirty(1) // q can differ from df in ARubin test } void boottest::setwillplot(real scalar willplot) { this.willplot = willplot } void boottest::setgrid(real rowvector gridmin, real rowvector gridmax, real rowvector gridpoints) { this.gridmin = gridmin; this.gridmax = gridmax; this.gridpoints = gridpoints } void boottest::setmadjust(string scalar madjtype, real scalar NumH0s) { this.madjtype = strlower(madjtype) this.NumH0s = NumH0s if (this.madjtype != "bonferroni" & this.madjtype != "sidak" & this.madjtype != "") _error(198, "Multiple-hypothesis adjustment type must be "Bonferroni" or "Sidak"."') } void boottest::setauxwttype(string scalar auxwttype) { auxwttype = strlower(auxwttype) if (.==(this.auxwttype = auxwttype=="rademacher" ? 0 : (auxwttype=="mammen" ? 1 : (auxwttype=="webb" ? 2 : (auxwttype=="normal" ? 3 : (auxwttype=="gamma" ? 4 : .)))))) _error(198, "Wild type must be "Rademacher", "Mammen", "Webb", "Normal", or "Gamma"."') setdirty(1) } void boottest::setMaxMatSize(real scalar MaxMatSize) { this.MaxMatSize = MaxMatSize; setdirty(1) } real colvector boottest::getdist(| string scalar diststat) { pointer (real rowvector) scalar _pnumer if (dirty) boottest() if (diststat == "numer") { _pnumer = v_sd==1? pnumer : &(*pnumer / v_sd) _sort( DistCDR = (*_pnumer)[|2\.|]' /*:+ *pr*/ , 1) } else if (rows(DistCDR)==0) if (cols(*pDist) > 1) _sort( DistCDR = multiplier * (*pDist)[|2\.|]' , 1) else DistCDR = J(0,1,0) return(DistCDR) } // get p value. Robust to missing bootstrapped values interpreted as +infinity. real scalar boottest::getp(|real scalar classical) { real scalar tmp; real scalar _p if (dirty) boottest() tmp = (*pDist)[1] if (tmp == .) return (.) if (B & classical==.) if (sqrt & ptype != 3) { if (ptype==0) p = rowsum(-abs(tmp) :> -abs(*pDist)) / BFeas // symmetric p value; do so as not to count missing entries in *pDist else if (ptype==1) // equal-tail p value p = 2 * min((rowsum(tmp :> *pDist) , rowsum(-tmp:>- *pDist))) / BFeas else p = rowsum( tmp :> *pDist) / BFeas // lower-tailed p value } else p = rowsum(-tmp :> - *pDist) / BFeas // upper-tailed p value or p value based on squared stats else { tmp = tmp * multiplier _p = small? Ftail(df, df_r, sqrt? tmp*tmp : tmp) : chi2tail(df, sqrt? tmp*tmp : tmp) if (sqrt & twotailed==0) { _p = _p / 2 if ((ptype==3) == (tmp<0)) _p = 1 - _p } if (classical != .) return(_p) p = _p // only save as the official p if this was not requested by plot() for AR case } return(p) } // numerator for full-sample test stat real colvector boottest::getb() { if (dirty) boottest() return(v_sd == 1? (*pnumer)[,1] : (*pnumer)[,1] / v_sd) } // denominator for full-sample test stat real matrix boottest::getV() { if (dirty) boottest() return (statDenom / ((v_sd == 1? smallsample : v_sd * v_sd * smallsample) * (sqrt? multiplier*multiplier : multiplier) * df)) } // wild weights real matrix boottest::getv() return(v_sd==1? v[|.,2\.,.|] : v[|.,2\.,.|] / v_sd) // Return number of bootstrap replications with feasible results // Returns 0 if getp() not yet accessed, or doing non-bootstrapping tests real scalar boottest::getrepsFeas() return (BFeas) real scalar boottest::getNBootClust() return (Nstar) // return number of replications, possibly reduced to 2^G real scalar boottest::getreps() return (B) real scalar boottest::getpadj(|real scalar classical) { real scalar _p _p = dirty | classical != . ? getp(classical) : p if (madjtype=="bonferroni") return(min((1, NumH0s*_p))) if (madjtype=="sidak" ) return(1 - (1 - _p)^NumH0s) return(_p) } real scalar boottest::getstat() { if (dirty) boottest() return(multiplier * (*pDist)[1]) } real scalar boottest::getdf() { if (dirty) boottest() return(df) } real scalar boottest::getdf_r() { if (dirty) boottest() return(df_r) } real matrix boottest::getplot() { if (notplotted) plot() return((plotX,plotY)) } real rowvector boottest::getpeak() { // x and y values of confidence curve peak (at least in OLS & ARubin) if (notplotted) plot() return(peak) } real matrix boottest::getCI() { if (notplotted) plot() return(CI) } void boottest::_st_view(real matrix V, real scalar i, string rowvector j, string scalar selectvar) { if (favorspeed() | 1) V = length(tokens(j))? st_data(i, j, selectvar) : st_data(i, J(1,0,0), selectvar) else st_view(V, i, j, selectvar) } // helper for summing over clusterings while factoring in clustering-specific parity and small-sample adjustments // replace X with Y if c=1; otherwise add it void boottest::_clustAccum(real matrix X, real scalar c, real matrix Y) X = c == 1? (Clust.even? (Clust. multiplier != 1? Clust. multiplier * Y : Y) : (Clust. multiplier != 1? (-Clust. multiplier) * Y : -Y)) : (Clust[c].even? X + (Clust[c].multiplier != 1? Clust[c].multiplier * Y : Y) : X - (Clust[c].multiplier != 1? Clust[c].multiplier * Y : Y)) // efficiently store bootstrap results for one wild weight group in a matrix, one col per replication, handling Nw>1 case (matsizegb() option) void boottest::storeWtGrpResults(pointer(real matrix) scalar pdest, real scalar w, real matrix content) if (Nw==1) pdest = &content else (*pdest)[|., WeightGrpStart[w] \ ., WeightGrpStop[w]|] = content void boottest::Init() { // for efficiency when varying r repeatedly to make CI, do stuff once that doesn't depend on r real colvector sortID, o, _FEID, _sqrtwt pointer (real colvector) scalar pIDAllData, pIDCapData real rowvector ClustCols real matrix Combs, tmp, IDCap real scalar i, j, c, minN, _B, i_FE, sumFEwt, _df pragma unset pIDAllData; pragma unset pIDCapData Nobs = rows(*py1) NClustVar = cols(*pID) kX = (kX1 = cols(*pX1)) + (kX2 = cols(*pX2)) if (kX2 == 0) pX2 = &J(Nobs,0,0) if ((kY2 = cols(*pY2)) == 0) pY2 = &J(Nobs,0,0) kZ = kX1 + kY2 if (LIML & kX2 == kY2) { // exactly identified LIML = 2SLS kappa = 1 LIML = 0 } if ((REst = rows(*pR1)) == 0) { // base model contains no restrictions? pR1 = &J(0,kZ,0) pr1 = &J(0,1 ,0) } if (kappa==.) kappa = kX2>0 // if kappa in kappa-class estimation not specified, it's 0 or 1 for OLS or 2SLS WRE = (kappa & scoreBS==0) | ARubin WREnonARubin = WRE & ARubin==0 if (haswt = rows(*pwt)>1) { sumwt = sum(*pwt) sqrtwt = sqrt(*pwt) py1 = &(*py1 :* sqrtwt) if (kY2) pY2 = &(*pY2 :* sqrtwt) if (kX1) pX1 = &(*pX1 :* sqrtwt) if (kX2) pX2 = &(*pX2 :* sqrtwt) } else pwt = &(sumwt = 1) _Nobs = haswt & obswttype=="fweight"? sumwt : Nobs if (WREnonARubin) if (NClustVar) infoBootData = _panelsetup(*pID, 1..NBootClustVar, pIDBootData) else { pIDBootData = &(1::Nobs) pinfoCapData = &(infoBootData = *pIDBootData,*pIDBootData) // no clustering, so no collapsing by cluster } else if (NClustVar) infoBootData = _panelsetup(*pID, 1..min((NClustVar, NBootClustVar))) // bootstrap cluster grouping defs rel to original data else pinfoCapData = pinfoAllData = &(infoBootData = J(Nobs,0,0)) // causes no collapsing of data in _panelsum() calls, only multiplying by weights if any Nstar = rows(infoBootData) if (NClustVar > NBootClustVar) { // info for grouping by intersections of all bootstrap & clustering vars wrt data; used to speed crosstab UXAR wrt bootstrapping cluster & intersection of all error clusters pinfoAllData = &_panelsetup(*pID, 1..NClustVar, pIDAllData) if (subcluster & rows(infoBootData) != rows(*pinfoAllData)) { errprintf("\nboottest can only perform the subcluster bootstrap when the bootstrap clusters are nested within the (intersections of the) error clusters.\n") _error(499) } } else { pinfoAllData = &infoBootData // info for grouping by intersections of all bootstrap & clustering vars wrt data; used to speed crosstab UXAR wrt bootstrapping cluster & intersection of all error clusters if (WREnonARubin) pIDAllData = pIDBootData } subcluster = NClustVar > NErrClustVar bootneqcap = NBootClustVar!=NClustVar | subcluster // bootstrapping and intersection of error clusterings different? if (subcluster) // info for intersections of error clustering wrt data pinfoCapData = &_panelsetup(*pID, subcluster+1..NClustVar, pIDCapData) else { pinfoCapData = pinfoAllData // info for intersections of error clustering wrt data if (WREnonARubin) pIDCapData = pIDAllData } Nall = rows(*pinfoAllData) IDCap = rows(*pinfoCapData)==Nobs? *pID : (*pID)[(*pinfoCapData)[,1],] // version of ID matrix with one row for each all-error-cluster-var intersection instead of 1 row for each obs; gets resorted pIDAll = Nall==Nobs ? pID : &((*pID)[(*pinfoAllData)[,1],]) // version of ID matrix with one row for each all-bootstrap & error cluster-var intersection instead of 1 row for each obs if (bootstrapt) { if (NClustVar) { minN = . Combs = combs(NErrClustVar) // represent all error clustering combinations. First is intersection of all error clustering vars Clust = structboottestClust(rows(Combs)-1) // leave out no-cluster combination NErrClustCombs = length(Clust) BootClust = 2^(NClustVar - NBootClustVar) // location of bootstrap clustering within list of cluster combinations for (c=1; c<=NErrClustCombs; c++) { // for each error clustering combination ClustCols = subcluster :+ _selectindex(Combs[c,]) Clust[c].even = mod(cols(ClustCols),2) if (c == 1) if (subcluster) { IDCap = IDCap[ Clust.order = stableorder(IDCap, ClustCols), ] Clust.info = _panelsetup(IDCap, ClustCols) } else Clust.info = J(rows(*pinfoAllData),0,0) // causes no collapsing of data in _panelsum() calls else { if (any(Combs[|c, min(_selectindex(Combs[c,] :!= Combs[c-1,])) \ c,.|])) // if this sort ordering same as last to some point and missing thereafter, no need to re-sort IDCap = IDCap[ Clust[c].order = stableorder(IDCap, ClustCols), ] Clust[c].info = _panelsetup(IDCap, ClustCols) } Clust[c].N = rows(Clust[c].info) if (small) { Clust[c].multiplier = Clust[c].N/(Clust[c].N-1) if (Clust[c].N < minN) minN = Clust[c].N } else Clust[c].multiplier = 1 } if (scoreBS | WREnonARubin==0) ClustShare = haswt? *_panelsum(*pwt, *pinfoCapData)/sumwt : ((*pinfoCapData)[,2]-(*pinfoCapData)[,1]:+ 1)/Nobs // share of observations by group } else { // if no clustering, cast "robust" as clustering by observation Clust = structboottestClust() Clust.multiplier = small? _Nobs / (_Nobs - 1) : 1 Clust.even = 1 Clust.N = Nobs Clust.info = J(Nobs, 0, 0) // signals _panelsum not to aggregate NErrClustCombs = 1 if (scoreBS | WREnonARubin==0) ClustShare = haswt? *pwt/sumwt : 1/_Nobs } } else minN = rows(infoBootData) purerobust = robust & (scoreBS | subcluster)==0 & Nstar==Nobs // do we ever error- *and* bootstrap-cluster by individual? granular = WREnonARubin? 2*Nobs*B*(2*Nstar+1) < Nstar*(Nstar*Nobs+Nall*B*(Nstar+1)) : !jk & robust & scoreBS==0 & (purerobust | (Nall+Nstar)*kZ*B + (Nall-Nstar)*B + kZ*B < Nall*kZ*kZ + Nobs*kZ + Nall * Nstar*kZ + Nall*Nstar) if (jk & !WREnonARubin) granularjk = kZ^3 + Nstar * (Nobs/Nstar*kZ^2 + (Nobs/Nstar)^2*kZ + (Nobs/Nstar)^2 + (Nobs/Nstar)^3) < Nstar * (kZ^2*Nobs/Nstar + kZ^3 + 2*kZ*(kZ + Nobs/Nstar)) if (bootstrapt) { if (robust & purerobust==0) { if (subcluster | granular) infoErrAll = _panelsetup(*pIDAll, subcluster+1..NClustVar) // info for error clusters wrt data collapsed to intersections of all bootstrapping & error clusters; used to speed crosstab UXAR wrt bootstrapping cluster & intersection of all error clusterings if ((scoreBS & B) | (WREnonARubin & granular==0 & !jk & bootstrapt)) JNcapNstar = J(Clust.N, Nstar, 0) } if (WREnonARubin & robust & bootstrapt & granular==0) { if (cols(*pIDAllData) == 0) { (void) _panelsetup(*pID, 1..NClustVar, pIDAllData) (void) _panelsetup(*pID, subcluster+1..NClustVar, pIDCapData) } IDCTCapstar = infoCTCapstar = smatrix(Nstar) for (i=Nstar;i;i--) { tmp = (*pIDAllData)[|infoBootData[i,]'|] // ID numbers w.r.t. intersection of all bootstrap/error clusterings contained in bootstrap cluster i infoCTCapstar[i].M = (*pinfoAllData)[tmp[1]::tmp[rows(tmp)],] // for each of those ID's, panel info for the all-bootstrap/error-clusterings data row groupings IDCTCapstar[i].M = (*pIDCapData)[infoCTCapstar[i].M[,1]] // ID numbers of those groupings w.r.t. the all-error-clusterings grouping } } } if (NFE) { // identify FE groups sortID = (*pFEID)[o = stableorder(*pFEID, 1)] i_FE = 1; FEboot = B>0 & NClustVar; j = Nobs; _FEID = J(Nobs, 1, 1) invFEwt = J(NFE,1,0) FEs = structFE(NFE) for (i=Nobs-1;i;i--) { if (sortID[i] != sortID[i+1]) { FEs[i_FE].is = o[|i+1\j|] if (haswt) { _sqrtwt = sqrtwt[FEs[i_FE].is] FEs[i_FE].wt = _sqrtwt / (sumFEwt = colsum((*pwt)[FEs[i_FE].is])) FEs[i_FE].sqrtwt = _sqrtwt } else FEs[i_FE].wt = J(j-i, 1, 1/(sumFEwt = j-i)) if ((B & robust & granular < NErrClustVar) | (WREnonARubin & robust & granular & bootstrapt)) invFEwt[i_FE] = 1 / sumFEwt j = i if (!WREnonARubin & FEboot) { // are all of this FE's obs in same bootstrapping cluster? (But no need to check if B=0 for then CT_WE in 2nd term of (62) orthogonal to v = col of 1's) tmp = (*pID)[FEs[i_FE].is, 1..NBootClustVar] FEboot = all(tmp :== tmp[1,]) } ++i_FE } _FEID[o[i]] = i_FE } FEs[NFE].is = o[|.\j|] if (haswt) { _sqrtwt = sqrtwt[FEs[NFE].is] FEs[NFE].wt = _sqrtwt / (sumFEwt = colsum((*pwt)[FEs[NFE].is])) FEs[NFE].sqrtwt = _sqrtwt } else FEs[NFE].wt = J(j,1,1/(sumFEwt = j)) if (robust & ((B & granular < NErrClustVar) | (WREnonARubin & granular & bootstrapt))) invFEwt[NFE] = 1 / sumFEwt if (FEboot) { // are all of this FE's obs in same bootstrapping cluster? tmp = (*pID)[FEs[NFE].is, 1..NBootClustVar] FEboot = all(tmp :== tmp[1,]) } pFEID = &_FEID // ordinal fixed effect ID pX1 = partialFE(pX1) pX2 = partialFE(pX2) py1 = partialFE(py1) pY2 = partialFE(pY2) } if (B & robust & granular & purerobust==0 & bootstrapt & WREnonARubin==0) if (NFE & FEboot==0) (void) _panelsetup(*pID , 1..NBootClustVar, pIDBootData) else (void) _panelsetup(*pIDAll, 1..NBootClustVar, pIDBootAll ) if (enumerate = (B & auxwttype==0 & Nstar*ln(2) < ln(B)+1e-6)) // generate full Rademacher set? MaxMatSize = . Nw = MaxMatSize == .? 1 : ceil((B+1) * max((rows(*pIDBootData), rows(*pIDBootAll), Nstar)) * 8 / MaxMatSize / 1.0X+1E) // 1.0X+1E = giga(byte) if (Nw == 1) { MakeWildWeights(B, 1) // make all wild weights, once if (enumerate) B = cols(v) - 1 // replications reduced to 2^G WeightGrpStart = 1 } else { seed = rseed() _B = ceil((B+1) / Nw) Nw = ceil((B+1) / _B) WeightGrpStart = (0::Nw-1) * _B :+ 1 (WeightGrpStop = (1::Nw ) * _B )[Nw] = B+1 } if (ML) df = rows(*pR) else { if (ARubin) { pR = &(J(kX2,kX1,0), I(kX2)) // attack surface is all endog vars pR1 = kX1 & rows(*pR1)? &((*pR1)[|.,.\.,kX1|], J(rows(*pR1),kX2,0)) : &J(0, kX, 0) // and convert model constraints from referring to X1, Y2 to X1, X2 } df = rows(*pR) if (WRE==0 & kappa==0) { // regular OLS DGP.parent = Repl.parent = &this DGP.LIML = this.LIML; DGP.Fuller = this.Fuller; DGP.kappa = this.kappa DGP.SetR(null? *pR1 \ *pR : *pR1) // DGP constraints: model constraints + null if imposed Repl.SetR(*pR1) // model constraints only DGP.InitVars(&Repl.R1perp) DGP.InitTestDenoms() pM = &DGP // estimator object from which to get A, AR, XAR } else if (ARubin) { if (willplot) { // for plotting/CI purposes get original point estimate since not normally generated DGP = boottestIVGMM(); DGP.parent = &this DGP.LIML = this.LIML; DGP.Fuller = this.Fuller; DGP.kappa = this.kappa DGP.SetR(*pR1, J(0,kZ,0)) // no-null model DGP.InitVars() DGP.Estimate(0, *pr1) confpeak = DGP.beta.M // estimated coordinate of confidence peak } DGP = boottestARubin(); DGP.parent = &this DGP.SetR(*pR1) DGP.InitVars() DGP.InitTestDenoms() pM = &DGP // estimator object from which to get A, AR, XAR kZ = kX } else if (WREnonARubin) { DGP = boottestIVGMM(); DGP.parent = &this DGP.LIML = 1 if (null) DGP.SetR(*pR1 \ *pR, J(0,kZ,0)) // DGP constraints: model constraints + imposed null else DGP.SetR(*pR1 , *pR) // when null not imposed, keep it in the attack surface, though not used there, so Zperp is same in DGP and Repl DGP.InitVars() Repl = boottestIVGMM() Repl.parent = &this Repl.isDGP = 0 Repl.LIML = this.LIML; Repl.Fuller = this.Fuller; Repl.kappa = this.kappa Repl.SetR(*pR1, *pR) Repl.InitVars() Repl.Estimate(0, *pr1) if (null==0) { // if not imposing null, then DGP constraints, kappa, Hessian, etc. do not vary with r and can be set now DGP.Estimate(jk, *pr1) DGP.MakeResiduals(0) } InitWRE() } else { // the score bootstrap for IV/GMM uses a IV/GMM DGP but then masquerades as an OLS test because most factors are fixed during the bootstrap. To conform, need DGP and Repl objects with different R, R1, one with FWL, one not DGP = boottestIVGMM(); DGP.parent = &this DGP.LIML = this.LIML; DGP.Fuller = this.Fuller; DGP.kappa = this.kappa DGP.SetR(null? *pR1 \ *pR : *pR1, J(0,kZ,0)) // DGP constraints: model constraints + null if imposed DGP.InitVars() Repl = boottestIVGMM(); Repl.parent = &this Repl.LIML = this.LIML; Repl.Fuller = this.Fuller; Repl.kappa = this.kappa Repl.SetR(*pR1, I(kZ)) // process replication restraints = model constraints only Repl.InitVars(&Repl.R1perp) Repl.Estimate(0, *pr1) // bit inefficient to estimate in both objects, but maintains the conformity Repl.InitTestDenoms() pM = &Repl // estimator object from which to get A, AR, XAR; DGP follows WRE convention of using FWL, Repl follows OLS convention of not; scoreBS for IV/GMM mixes the two if (null==0) { // if not imposing null, then DGP constraints, kappa, Hessian, etc. do not vary with r and can be set now DGP.Estimate(jk, *pr1) DGP.MakeResiduals(0) } } } if (bootstrapt) { denom = smatrix(df,df) if (WREnonARubin==0 & robust) { if (B) Kd = smatrix(df) Kcd = smatrix(NErrClustCombs, df) pJcd = B? &smatrix(NErrClustCombs, df) : &Kcd // if B = 0, Kcd will be multiplied by v, which is all 1's, and will constitute Jcd } } if (df==1) setsqrt(1) // work with t/z stats instead of F/chi2 if (small) { _df = _Nobs - kZ + cols(Repl.R1invR1R1) - FEdfadj df_r = NClustVar? minN - 1 : _df multiplier = (smallsample = _df / (_Nobs - robust)) / df // divide by # of constraints because F stat is so defined } else multiplier = smallsample = 1 if ((robust | ML)==0) multiplier = multiplier * _Nobs // will turn sum of squared errors in denom of t/z into mean if (sqrt) multiplier = sqrt(multiplier) if (bootstrapt & (WREnonARubin | df > 1+robust | MaxMatSize<.)) // unless nonWRE or df=1 or splitting weight matrix, code will create Dist element-by-element, so pre-allocate vector now pDist = &J(1, B+1, .) if (Nw>1 | WREnonARubin | (null==0 & df<=2)) pnumer = &J(df, B+1, .) if (WREnonARubin==0) { poles = anchor = J(0,0,0) interpolate_u = !(robust | ML) if (interpolable = bootstrapt & B & null & Nw==1 & (kappa==0 | ARubin)) { dnumerdr = smatrix(q) if (interpolate_u) dudr = dnumerdr if (robust) { ddenomdr = dJcddr = ssmatrix(q) ddenomdr2 = ssmatrix(q, q) for (i=q;i;i--) { ddenomdr[i].M = smatrix(df,df) for (j=i;j;j--) ddenomdr2[i,j].M = ddenomdr[i].M dJcddr[i].M = smatrix(NErrClustCombs, df) } } } } } // main routine void boottest::boottest() { real scalar w if (initialized==0) Init() else if (null==0) { NoNullUpdate() return } if (Nw > 1) { rseed(seed) MakeWildWeights(WeightGrpStop[1] - 1, 1) } if (WREnonARubin) PrepWRE() else MakeInterpolables() // make stuff that depends linearly on r, possibly by interpolating, for first weight group for (w=1; w<=Nw; w++) { // do group 1 first because it includes col 1, which is all that might need updating in constructing CI in WCU if (w > 1) MakeWildWeights(WeightGrpStop[w] - WeightGrpStart[w] + 1, 0) if (WREnonARubin) MakeWREStats(w) else MakeNonWREStats(w) } if (bootstrapt==0) UpdateBootstrapcDenom() BFeas = (*pDist)[1]==.? 0 : rownonmissing(*pDist) - 1 DistCDR = J(0,0,0) setdirty(0) initialized = 1 } // if not imposing null and we have returned to boottest(), then df=1 or 2; we're plotting or finding CI, and only test stat, not distribution, changes with r void boottest::NoNullUpdate() { if (WREnonARubin) (*pnumer)[,1] = Repl.RRpar * betas[,1] - *pr else if (ARubin) { DGP.Estimate(0, *pr) (*pnumer)[,1] = v_sd * DGP.Rpar * DGP.beta.M[|kX1+1\.|] // coefficients on excluded instruments in ARubin OLS } else (*pnumer)[,1] = v_sd * (*pR * (ML? beta : pM->Rpar * pM->beta.M) - *pr) // Analytical Wald numerator; if imposing null then numer[,1] already equals this. If not, then it's 0 before this (*pDist)[1] = df==1? (*pnumer)[1] / sqrt(statDenom) : (*pnumer)[,1] ' invsym(statDenom) * (*pnumer)[,1] } // compute bootstrap-c denominator from all bootstrap numerators // for 1-dimensional hypothesis, can ignore denominators (Young 2022, eq. 16), except that user might request actual distribution void boottest::UpdateBootstrapcDenom() { real colvector numer1 numer1 = (*pnumer)[,1] numersum = rowsum(*pnumer) - numer1 statDenom = (*pnumer * *pnumer' - numer1 * numer1' - numersum * numersum' / B) / B pDist = sqrt? &(*pnumer:/sqrt(statDenom)) : &colsum(*pnumer :* invsym(statDenom) * *pnumer) } // draw wild weight matrix of width _B. If first=1, insert column of 1s at front. For non-Anderson-Rubin WRE, subtract 1 from all weights void boottest::MakeWildWeights(real scalar _B, real scalar first) { real scalar m m = WREnonARubin & jk & small? sqrt(1 - 1 / Nstar) : 1 // finite-sample adjuster for JK regressions if (_B) { // in scoretest or waldtest WRE, still make v a col of 1's if (enumerate) v = J(Nstar,1,1), count_binary(Nstar, -m, m) // complete Rademacher set else if (auxwttype==3) v = rnormal(Nstar, _B+first, 0, m) // normal weights else if (auxwttype==4) v = m * (rgamma(Nstar, _B+first, 4, .5) :- 2) // Gamma weights else if (auxwttype==2) if (WREnonARubin) v = sqrt((2*m) * ceil(runiform(Nstar, _B+first) * 3)) :* ((runiform(Nstar, _B+first):>=.5):-.5) // Webb weights else { v = sqrt( ceil(runiform(Nstar, _B+first) * 3)) :* ((runiform(Nstar, _B+first):>=.5):-.5) // Webb weights, divided by sqrt(2) v_sd = 1.6a09e667f3bcdX-001 /*sqrt(.5)*/ } else if (auxwttype) if (WREnonARubin) v = ( rdiscrete(Nstar, _B+first, 1.727c9716ffb76X-001 \ 1.1b06d1d200914X-002 /*.5+sqrt(.05) \ .5-sqrt(.05)*/) :- 1.5 ) * (m * 1.1e3779b97f4a8X+001) /*sqrt(5)*/ :+ .5 // Mammen weights, minus 1 for convenience in WRE else { v = ( rdiscrete(Nstar, _B+first, 1.727c9716ffb76X-001 \ 1.1b06d1d200914X-002 /*.5+sqrt(.05) \ .5-sqrt(.05)*/) :- 1.5 ) :+ 1.c9f25c5bfedd9X-003 /*.5/sqrt(5)*/ // Mammen weights, divided by sqrt(5) v_sd = 1.c9f25c5bfedd9X-002 /*sqrt(.2)*/ } else if (WREnonARubin) v = jk? m :+ (runiform(Nstar, _B+first) :< .5) * (-2 * m) : (runiform(Nstar, _B+first) :< .5) * -2 :+ 1 // Rademacher weights; 1st exp simplifies to 2nd when jk=0 but is slower else { v = (runiform(Nstar, _B+first) :>= .5) :- .5 // Rademacher weights, divided by 2 v_sd = .5 } if (first) v[,1] = J(Nstar, 1, v_sd) // keep original residuals in first entry to compute base model stat } else v = J(0,1,0) // in places, cols(v) indicates B -- 1 for classical tests } void boottest::InitWRE() { // stuff done only once that knits together results from DGP and Repl regression prep real scalar i, g, j; real matrix X1g, X2g, Zperpg, X1parg, Zg, ZR1g, Y2g, S, Y2gi; real colvector y1g; real rowvector X1pargi if (LIML & Repl.kZ==1 & Nw==1) As = betas = J(1, B+1, 0) SstarUYbar = SstarinvZperpZperpZperpU = SstarZperpU = SstarUXinvXX = SstarUX = smatrix(Repl.kZ+1) SstarUU = smatrix(Repl.kZ+1, Repl.kZ+1) if (bootstrapt) { deltadenom_b = J(Repl.kZ, Repl.kZ, 0) Zyi = deltadenom = smatrix(Repl.kZ+1) SstarUMZperp = SstarUPX = SstarUX _Jcap = J(Clust.N, Repl.kZ, 0) if (granular==0) SCTcapUXinvXX = smatrix(Repl.kZ+1, Nstar) if (LIML | robust==0) { YYstar = YPXYstar = SstarUX YYstar_b = YPXYstar_b = J(Repl.kZ+1, Repl.kZ+1, 0) } if (NFE & FEboot==0 & (bootstrapt | kappa != 1 | LIML)) CTstarFEU = SstarUX negSstarUMZperpX = smatrix(Repl.kZ+1, Clust.N) } if (granular==0) { if (jk==0) { SstarZR1X = SstarXy1 = SstarZX = SstarY2X = SstarXY2par = smatrix(Nstar) if (kappa!=1 | LIML | bootstrapt) SstarZR1Zperp = SstarZperpy1 = SstarZZperp = SstarZperpY2par = SstarY2Zperp = smatrix(Nstar) if (kappa!=1 | LIML | robust==0 ) { SstarZR1Z = SstarZR1ZR1 = SstarX1parY2par = SstarXX1par = SstarX1pary1 = SstarY2X1par = SstarY2Y2 = SstarZR1X1par = SstarZX1par = SstarY2parY2par = Sstary1Y2par = Sstary1Y2 = SstarZR1Y2par = SstarZR1Y2 = SstarZY2par = SstarZY2 = SstarY2Y2par = SstarZR1y1 = SstarZy1 = SstarZZ = smatrix(Nstar) Sstary1y1 = J(Nstar, 1, 0) for (i=Repl.kZ; i>=0; i--) { SstarUYbar[i+1].M = J(Nstar, 1+Repl.kZ, 0) for (j=i; j>=0; j--) SstarUU[i+1,j+1].M = J(Nstar, 1, 0) } } if (!bootneqcap) { pSstarXX = &ScapXX pSstarXZperp = &ScapXZperp } for (i=Repl.kZ; i>=0; i--) SstarUX[i+1].M = J(DGP.kX, Nstar, 0) if (kappa!=1 | LIML | bootstrapt) for (i=Repl.kZ; i>=0; i--) SstarZperpU[i+1].M = J(cols(*DGP.pZperp), Nstar, 0) } if (robust & bootstrapt) { ScapXYbar = smatrix(Repl.kZ+1) if (jk) { CT_FEcapYbar = ScapPXYbarZperp = smatrix(Repl.kZ+1) FillingT0 = smatrix(Repl.kZ+1, Repl.kZ+1) // fixed component of groupwise term in sandwich filling } invXXXX1par = DGP.invXX * (cross(*DGP.pX1, *Repl.pX1par) \ cross(*DGP.pX2, *Repl.pX1par)) for (i=Repl.kZ; i>=0; i--) ScapXYbar[i+1].M = J(DGP.kX, Clust.N, 0) if (jk) for (i=Repl.kZ; i; i--) { ScapPXYbarZperp[i+1].M = J(Clust.N, cols(*DGP.pZperp), 0) for (j=i; j>=0; j--) FillingT0[i+1,j+1].M = J(Clust.N,1,0) } ScapXX = ScapXZperp = smatrix(Clust.N); ScapYX = smatrix(Repl.kZ + 1); ScapXX1par = smatrix(Clust.N, cols(Repl.Rpar)) for (i=Repl.kZ;i>=0;i--) ScapYX[i+1].M = J(Clust.N, DGP.kX, 0) for (g=Clust.N;g;g--) { S = (*pinfoCapData)[g,]' X1g = *pXS(*DGP.pX1, S) X2g = *pXS(*DGP.pX2, S) y1g = *pXS(*DGP.py1, S) Y2g = *pXS(*DGP.pY2, S) X1parg = *pXS(*Repl.pX1par,S) Zperpg = *pXS(*DGP.pZperp,S) ScapXX[g].M = cross(X1g, X2g); ScapXX[g].M = cross(X1g, X1g), ScapXX[g].M \ ScapXX[g].M', cross(X2g, X2g) ScapXZperp[g].M = cross(X1g,Zperpg) \ cross(X2g,Zperpg) for (i=cols(ScapXX1par);i;i--) { X1pargi = X1parg[,i] ScapXX1par[g,i].M = cross(X1g, X1pargi) \ cross(X2g, X1pargi) } ScapYX.M[g,] = cross(y1g, X1g) , cross(y1g, X2g) for (i=Repl.kZ;i;i--) { Y2gi = Y2g * Repl.RparY[,i] ScapYX[i+1].M[g,] = cross(Y2gi, X1g) , cross(Y2gi, X2g) } if (jk==0 & bootneqcap==0) { // if bootstrapping and intersection-of-error clusterings same, and using optimized non-jk code, then integrate by-boot and by-cap constructions for speed y1g = *pXS(*DGP.py1,S) Zg = *pXS(DGP.Z,S) ZR1g = *pXS(*DGP.pZR1,S) Y2g = *pXS(*DGP.pY2 ,S) SstarZR1X[g].M = cross(ZR1g, X1g) , cross(ZR1g, X2g) SstarXy1[g].M = cross(X1g, y1g) \ cross(X2g, y1g) SstarZX[g].M = cross(Zg, X1g) , cross(Zg, X2g) SstarXY2par[g].M = (SstarY2X[g].M = cross(Y2g, X1g) , cross(Y2g, X2g)) ' Repl.RparY if (kappa!=1 | LIML | bootstrapt) { SstarZR1Zperp[g].M = cross(ZR1g, Zperpg) SstarZperpy1[g].M = cross(Zperpg, y1g) SstarZZperp[g].M = cross(Zg, Zperpg) SstarZperpY2par[g].M = (SstarY2Zperp[g].M = cross(Y2g, Zperpg)) ' Repl.RparY } if (kappa!=1 | LIML | robust==0) { X1parg = *pXS(*Repl.pX1par,S) SstarX1parY2par[g].M = (SstarY2X1par[g].M = cross(Y2g, X1parg)) ' Repl.RparY SstarZX1par[g].M = cross(Zg, X1parg) SstarZR1X1par[g].M = cross(ZR1g, X1parg) SstarXX1par[g].M = cross(X1g, X1parg) \ cross(X2g, X1parg) SstarX1pary1[g].M = cross(X1parg, y1g) SstarY2parY2par[g].M = Repl.RparY ' cross(Y2g, Y2g) * Repl.RparY Sstary1Y2par[g].M = (Sstary1Y2[g].M = cross(y1g, Y2g)) * Repl.RparY SstarZR1Y2par[g].M = (SstarZR1Y2[g].M = cross(ZR1g, Y2g)) * Repl.RparY SstarZY2par[g].M = (SstarZY2[g].M = cross(Zg, Y2g)) * Repl.RparY SstarY2Y2par[g].M = cross(Y2g, Y2g) * Repl.RparY Sstary1y1[g] = cross(y1g, y1g) SstarZR1y1[g].M = cross(ZR1g, y1g) SstarZy1[g].M = cross(Zg, y1g) SstarZR1ZR1[g].M = cross(ZR1g, ZR1g) SstarZR1Z[g].M = cross(ZR1g, Zg) SstarZZ[g].M = cross(Zg, Zg) SstarY2Y2[g].M = cross(Y2g, Y2g) } } } if (jk==0 & robust & bootstrapt) { SallXu1ddot = SallXU2RparY = SallU2X = smatrix(Nall) if (bootneqcap) { pSallXy1 = &smatrix(Nall); pSallZX = &smatrix(Nall); pSallZR1X = &smatrix(Nall); pSallY2X = &smatrix(Nall); pSallXX = &smatrix(Nall) for (g=Nall;g;g--) { S = (*pinfoAllData)[g,]' X1g = *pXS(*DGP.pX1,S) X2g = *pXS(*DGP.pX2 ,S) y1g = *pXS(*DGP.py1,S) Zg = *pXS(DGP.Z,S) ZR1g = *pXS(*DGP.pZR1,S) Y2g = *pXS(*DGP.pY2 ,S) (*pSallXy1)[g].M = cross(X1g, y1g) \ cross(X2g, y1g) (*pSallZX)[g].M = cross(Zg, X1g), cross(Zg, X2g) (*pSallZR1X)[g].M = cross(ZR1g, X1g), cross(ZR1g, X2g) (*pSallY2X)[g].M = cross(Y2g, X1g), cross(Y2g, X2g) (*pSallXX)[g].M = cross(X1g, X2g); (*pSallXX)[g].M = cross(X1g, X1g), (*pSallXX)[g].M \ (*pSallXX)[g].M', cross(X2g, X2g) } } else { pSallXy1 = &SstarXy1 pSallZX = &SstarZX pSallZR1X = &SstarZR1X pSallY2X = &SstarY2X pSallXX = pSstarXX } } } if (jk==0 & (bootneqcap | !(robust & bootstrapt))) { pSstarXX = &smatrix(Nstar) if (kappa!=1 | LIML | bootstrapt) pSstarXZperp = &smatrix(Nstar) for (g=Nstar;g;g--) { S = infoBootData[g,]' X1g = *pXS(*DGP.pX1,S) X2g = *pXS(*DGP.pX2,S) (*pSstarXX)[g].M = cross(X1g, X2g); (*pSstarXX)[g].M = cross(X1g, X1g), (*pSstarXX)[g].M \ (*pSstarXX)[g].M', cross(X2g, X2g) y1g = *pXS(*DGP.py1,S) Zg = *pXS(DGP.Z,S) ZR1g = *pXS(*DGP.pZR1,S) Y2g = *pXS(*DGP.pY2 ,S) SstarZR1X[g].M = cross(ZR1g, X1g) , cross(ZR1g, X2g) SstarXy1[g].M = cross(X1g, y1g) \ cross(X2g, y1g) SstarZX[g].M = cross(Zg, X1g) , cross(Zg, X2g) SstarXY2par[g].M = (SstarY2X[g].M = cross(Y2g, X1g) , cross(Y2g, X2g)) ' Repl.RparY if (kappa!=1 | LIML | bootstrapt) { Zperpg = *pXS(*DGP.pZperp, S) (*pSstarXZperp)[g].M = cross(X1g, Zperpg) \ cross(X2g, Zperpg) SstarZR1Zperp[g].M = cross(ZR1g, Zperpg) SstarZperpy1[g].M = cross(Zperpg, y1g) SstarZZperp[g].M = cross(Zg, Zperpg) SstarZperpY2par[g].M = (SstarY2Zperp[g].M = cross(Y2g, Zperpg)) ' Repl.RparY } if (kappa!=1 | LIML | robust==0) { X1parg = *pXS(*Repl.pX1par,S) SstarX1parY2par[g].M = (SstarY2X1par[g].M = cross(Y2g, X1parg)) ' Repl.RparY SstarZX1par[g].M = cross(Zg, X1parg) SstarZR1X1par[g].M = cross(ZR1g, X1parg) SstarXX1par[g].M = cross(X1g, X1parg) \ cross(X2g, X1parg) SstarX1pary1[g].M = cross(X1parg, y1g) SstarY2parY2par[g].M = cross(Y2g, Y2g) Sstary1Y2par[g].M = (Sstary1Y2[g].M = cross(y1g, Y2g)) * Repl.RparY SstarZR1Y2par[g].M = (SstarZR1Y2[g].M = cross(ZR1g, Y2g)) * Repl.RparY SstarZY2par[g].M = (SstarZY2[g].M = cross(Zg, Y2g)) * Repl.RparY SstarY2Y2par[g].M = cross(Y2g, Y2g) * Repl.RparY Sstary1y1[g] = cross(y1g, y1g) SstarZR1y1[g].M = cross(ZR1g, y1g) SstarZy1[g].M = cross(Zg, y1g) SstarZR1ZR1[g].M = cross(ZR1g, ZR1g) SstarZR1Z[g].M = cross(ZR1g, Zg) SstarZZ[g].M = cross(Zg, Zg) SstarY2Y2[g].M = cross(Y2g, Y2g) } } } if (jk==0 & robust & bootstrapt & NFE & FEboot==0) { CTstarFEY2 = smatrix(kY2); for (i=kY2;i;i--) CTstarFEY2[i].M = crosstabFE(*pcol(*DGP.pY2,i), infoBootData) CTstarFEX = smatrix(DGP.kX); for (i=DGP.kX1;i ;i--) CTstarFEX[i].M = crosstabFE(*pcol(*DGP.pX1,i ), infoBootData) for (i=DGP.kX;i>DGP.kX1;i--) CTstarFEX[i].M = crosstabFE(*pcol(*DGP.pX2 ,i-DGP.kX1), infoBootData) CTstarFEy1 = crosstabFE(*DGP.py1, infoBootData) CTstarFEZR1 = smatrix(cols(*DGP.pZR1)); for (i=cols(*DGP.pZR1);i;i--) CTstarFEZR1[i].M = crosstabFE(*pcol(*DGP.pZR1,i), infoBootData) CTstarFEZ = smatrix(DGP.kZ); for (i=DGP.kZ;i;i--) CTstarFEZ[i].M = crosstabFE(*pcol(DGP.Z,i), infoBootData) CTcapFEX = smatrix(DGP.kX); for (i=DGP.kX1;i ;i--) CTcapFEX[i].M = invFEwt :* crosstabFE(*pcol(*DGP.pX1,i ), *pinfoCapData) for (i=DGP.kX;i>DGP.kX1;i--) CTcapFEX[i].M = invFEwt :* crosstabFE(*pcol(*DGP.pX2 ,i-DGP.kX1), *pinfoCapData) } } } void boottest::PrepWRE() { real scalar i, j, g, h; pointer (real colvector) scalar pu; real matrix ZU2ddotpar, deltadddot, tmp; real colvector PiddotRparYi, F1_0, _r, Pidddot, PiddotdeltaY if (null) { DGP.Estimate(jk, _r = *pr1 \ *pr) DGP.MakeResiduals(0) } else _r = *pr1 // slightly inefficient to reassign every time XYbar = DGP.Xy1bar, (invXXXZbar = Repl.XZ - DGP.XU2ddot * Repl.RparY) invXXXZbar = Repl.invXX * invXXXZbar ZU2ddotpar = (Repl.ZY2 - Repl.XZ'DGP.Piddot) * Repl.RparY YbarYbar = DGP.deltadddot'Repl.XZ - DGP.y1barU2ddot * Repl.RparY YbarYbar = DGP.y1bary1bar, YbarYbar \ YbarYbar' , (Repl.ZZ - ZU2ddotpar' - ZU2ddotpar + Repl.RparY ' DGP.U2ddotU2ddot * Repl.RparY) if (granular==0) { PiddotRparY = DGP.Piddot * Repl.RparY PiddotdeltaY = DGP.Piddot * DGP.deltaY } if (robust & bootstrapt) { // for WRE replication regression, prepare for CRVE if (granular | (NFE & FEboot==0)) PXZbar = *pX12B(*Repl.pX1, *Repl.pX2, invXXXZbar) if (granular==0) { deltadddot = (DGP.perpRperpX'DGP.deltaX \ J(kX2, 1, 0)) + PiddotdeltaY for (g=Clust.N;g;g--) ScapXYbar.M[,g] = ScapXX[g].M * deltadddot // S_cap(M_Zperp*y1 :* P_(MZperpX)]) for (i=Repl.kZ; i; i--) { PiddotRparYi = PiddotRparY[,i] for (g=Clust.N;g;g--) ScapXYbar[i+1].M[,g] = ScapXX1par[g,i].M + ScapXX[g].M ' PiddotRparYi // S_cap(M_Zperp[Z or y1] :* P_(MZperpX)]) } if (jk) { for (i=Repl.kZ; i; i--) { PiddotRparYi = PiddotRparY[,i] F1_0 = invXXXX1par[,i] + PiddotRparYi for (g=Clust.N;g;g--) { ScapPXYbarZperp[i+1].M[g,] = F1_0 ' ScapXZperp[g].M // S_cap(M_Zperp*y1 :* P_(MZperpX)]) for (j=i; j; j--) FillingT0[i+1,j+1].M[g] = (ScapXX1par[g,j].M + ScapXX[g].M * PiddotRparY[,j]) ' F1_0 FillingT0 [i+1, 1].M[g] = DGP.deltadddot ' ScapXX[g].M * F1_0 } if (NFE & FEboot==0) CT_FEcapYbar[i+1].M = crosstabFE(*pcol(PXZbar,i), *pinfoCapData) :* invFEwt } } } } if (granular | jk) { // construct objects in a way that computes residuals in intermediate step (O(N)) pU2parddot = pXB(DGP.U2ddot, Repl.RparY) Zbar = *Repl.pX1par + DGP.Y2bar * Repl.RparY for (i=Repl.kZ; i>=0; i--) { // precompute various clusterwise sums pu = i? pcol(*pU2parddot,i) : &DGP.u1dddot // S_star(u :* X), S_star(u :* Zperp) for residuals u for each endog var; store transposed SstarUX [i+1].M = *_panelsum2(*Repl.pX1, *Repl.pX2, *pu, infoBootData)' SstarUXinvXX [i+1].M = Repl.invXX * SstarUX[i+1].M if (kappa!=1 | LIML | bootstrapt) { if (Repl.Yendog[i+1]) SstarinvZperpZperpZperpU[i+1].M = *Repl.pinvZperpZperp * (SstarZperpU[i+1].M = *_panelsum(*Repl.pZperp, *pu, infoBootData)') if (NFE & FEboot==0) CTstarFEU[i+1].M = crosstabFE(*pu, infoBootData) } if (kappa!=1 | LIML | robust==0) { SstarUYbar[i+1].M = *_panelsum2(DGP.y1bar, Zbar, *pu, infoBootData) for (j=i; j>=0; j--) SstarUU[i+1,j+1].M = *_panelsum(j? *pcol(*pU2parddot,j) : DGP.u1dddot, *pu, infoBootData) } if (robust & bootstrapt & granular==0 & Repl.Yendog[i+1]) { SstarUPX[i+1].M = Repl.XinvXX * SstarUX[i+1].M SstarUMZperp[i+1].M = *Repl.pZperp * SstarinvZperpZperpZperpU[i+1].M if (Nobs == Nstar) // subtract "crosstab" of observation by cap-group of u _diag(SstarUMZperp[i+1].M, diagonal(SstarUMZperp[i+1].M) - (i? *pcol(*pU2parddot,i) : DGP.u1dddot)) // case: bootstrapping by observation else if (i) for (g=Nobs;g;g--) SstarUMZperp[i+1].M[g,(*pIDBootData)[g]] = SstarUMZperp[i+1].M[g,(*pIDBootData)[g]] - (*pU2parddot)[g,i] else for (g=Nobs;g;g--) SstarUMZperp[i+1].M[g,(*pIDBootData)[g]] = SstarUMZperp[i+1].M[g,(*pIDBootData)[g]] - DGP.u1dddot[g] if (NFE & FEboot==0) SstarUMZperp[i+1].M = SstarUMZperp[i+1].M + (invFEwt :* CTstarFEU[i+1].M)[*pFEID,] // CT_(*,FE) (U ̈_(∥j) ) (S_FE S_FE^' )^(-1) S_FE } if (robust & bootstrapt & granular==0 & Repl.Yendog[i+1] & !(NClustVar == NBootClustVar & !subcluster)) // Within each bootstrap cluster, groupwise sum by all-error-cluster-intersections of u:*X and u:*Zperp (and times invXX or invZperpZperp) for (g=Nstar;g;g--) SCTcapUXinvXX[i+1,g].M = *_panelsum(Repl.XinvXX, *pu, infoCTCapstar[g].M) } } else { // for coarse clustering with no jackknife, construct objects while minimizing O(N) operations if (kappa!=1 | LIML | robust==0) Pidddot = (DGP.perpRperpX'DGP.deltaX \ J(kX2, 1, 0)) + PiddotdeltaY if (robust & bootstrapt) for (g=Nall;g;g--) { SallXU2RparY[g].M = (SallU2X[g].M = (*pSallY2X)[g].M - DGP.Piddot ' (*pSallXX)[g].M) ' Repl.RparY SallXu1ddot[g].M = (*pSallXy1)[g].M - (*pSallZX)[g].M ' DGP.beta.M + SallU2X[g].M ' DGP.deltaY - (*pSallZR1X)[g].M ' _r } for (g=Nstar;g;g--) SstarUX.M[,g] = SstarXy1[g].M - SstarZR1X[g].M ' _r - SstarZX[g].M ' DGP.beta.M + SstarY2X[g].M ' DGP.deltaY - (*pSstarXX)[g].M ' PiddotdeltaY for (i=Repl.kZ; i>=0; i--) { // precompute various clusterwise sums // S_star(u :* X), S_star(u :* Zperp) for residuals u for each endog var; store transposed if (i) { PiddotRparYi = PiddotRparY[,i] for (g=Nstar;g;g--) SstarUX[i+1].M[,g] = SstarXY2par[g].M[,i] - (*pSstarXX)[g].M * PiddotRparYi } SstarUXinvXX[i+1].M = Repl.invXX * SstarUX[i+1].M if (kappa!=1 | LIML | bootstrapt) { if (i) { if (Repl.Yendog[i+1]) for (g=Nstar;g;g--) SstarZperpU[i+1].M[,g] = SstarZperpY2par[g].M[,i] - (*pSstarXZperp)[g].M ' PiddotRparYi } else for (g=Nstar;g;g--) SstarZperpU.M[,g] = SstarZperpy1[g].M - SstarZR1Zperp[g].M ' _r - SstarZZperp[g].M ' DGP.beta.M + (SstarY2Zperp[g].M - DGP.Piddot ' (*pSstarXZperp)[g].M) ' DGP.deltaY SstarinvZperpZperpZperpU[i+1].M = *Repl.pinvZperpZperp * SstarZperpU[i+1].M if (NFE & FEboot==0) CTstarFEU[i+1].M = i? threebyone(CTstarFEY2, Repl.RparY[,i]) :- threebyone(CTstarFEX, PiddotRparYi) : CTstarFEy1 :- threebyone(CTstarFEZR1, _r) :- threebyone(CTstarFEZ, DGP.beta.M) :+ threebyone(CTstarFEY2, DGP.deltaY) :- threebyone(CTstarFEX, DGP.Piddot * DGP.deltaY) } if (kappa!=1 | LIML | robust==0) { if (i) for (g=Nstar;g;g--) { SstarUYbar[i+1].M[ g,1 ] = SstarUX[i+1].M[,g] ' Pidddot SstarUYbar[i+1].M[|g,2\g,.|] = SstarUX[i+1].M[,g] ' PiddotRparY + SstarX1parY2par[g].M[i,] - PiddotRparYi ' SstarXX1par[g].M } else for (g=Nstar;g;g--) { SstarUYbar[i+1].M[g,1] = SstarUX.M[,g] ' Pidddot SstarUYbar[i+1].M[|g,2\g,.|] = SstarUX.M[,g] ' PiddotRparY + SstarX1pary1[g].M' - _r ' SstarZR1X1par[g].M - DGP.beta.M ' SstarZX1par[g].M + DGP.deltaY ' SstarY2X1par[g].M - PiddotdeltaY ' SstarXX1par[g].M } if (i) for (j=i; j>=0; j--) if (j) for (g=Nstar;g;g--) SstarUU[i+1,j+1].M[g] = SstarY2parY2par[g].M[i,j] - (SstarXY2par[g].M[,i] + SstarUX[i+1].M[,g]) ' PiddotRparY[,j] else for (g=Nstar;g;g--) SstarUU[i+1,1].M[g] = Sstary1Y2par[g].M[,i] - _r ' SstarZR1Y2par[g].M[,i] - DGP.beta.M ' SstarZY2par[g].M[,i] + DGP.deltaY ' SstarY2Y2par[g].M[,i] - PiddotdeltaY ' SstarXY2par[g].M[,i] - SstarUX.M[,g] ' PiddotRparYi else for (g=Nstar;g;g--) SstarUU.M[g] = Sstary1y1[g] + _r ' SstarZR1ZR1[g].M * _r + DGP.beta.M ' SstarZZ[g].M * DGP.beta.M + DGP.deltaY ' SstarY2Y2[g].M * DGP.deltaY - PiddotdeltaY ' (*pSstarXX)[g].M * PiddotdeltaY + 2 * (_r ' (SstarZR1Z[g].M * DGP.beta.M - SstarZR1y1[g].M - SstarZR1Y2[g].M * DGP.deltaY) + Sstary1Y2[g].M * DGP.deltaY - PiddotdeltaY ' SstarUX.M[,g] - DGP.beta.M ' (SstarZy1[g].M + SstarZY2[g].M * DGP.deltaY)) } if (robust & bootstrapt & Repl.Yendog[i+1]) { for (g=Clust.N;g;g--) negSstarUMZperpX[i+1,g].M = ScapXZperp[g].M * SstarinvZperpZperpZperpU[i+1].M // S_* diag⁡(U ̈_(∥j) ) Z_⊥ (Z_⊥^' Z_⊥ )^(-1) Z_(⊥g)^' X_(∥g) // - S_* diag⁡(U ̈_(∥j) ) I_g^' X_(∥g) if (subcluster) // crosstab c,c* is wide for (g=Clust.N;g;g--) for (h=infoErrAll[g,1]; h<=infoErrAll[g,2]; h++) negSstarUMZperpX[i+1,g].M[,h] = negSstarUMZperpX[i+1,g].M[,h] - (i? SallXU2RparY[h].M[,i] : SallXu1ddot[h].M) else if (NClustVar == NBootClustVar) // crosstab c,c* is square for (g=Nstar;g;g--) negSstarUMZperpX[i+1,g].M[,g] = negSstarUMZperpX[i+1,g].M[,g] - (i? SallXU2RparY[g].M[,i] : SallXu1ddot[g].M) else // crosstab c,c* is tall for (h=Nstar;h;h--) for (g=Clust[BootClust].info[h,1]; g<=Clust[BootClust].info[h,2]; g++) negSstarUMZperpX[i+1,g].M[,h] = negSstarUMZperpX[i+1,g].M[,h] - (i? SallXU2RparY[g].M[,i] : SallXu1ddot[g].M) if (NFE & FEboot==0) for (h=DGP.kX;h;h--) { tmp = CTcapFEX[h].M ' CTstarFEU[i+1].M for (g=Clust.N;g;g--) negSstarUMZperpX[i+1,g].M[h,] = negSstarUMZperpX[i+1,g].M[h,] + tmp[g,] // CT_(*,FE) (U ̈_(∥j) ) (S_FE S_FE^' )^(-1) S_FE } } } } } // For WRE, and with reference to Y = [y1 Z], given 0-based column indexes within it, ind1, ind2, return all bootstrap realizations of Y[,ind1]'((1-kappa)*M_Zperp-kappa*P_Xpar)*Y[,ind2] for kappa constant across replications // ind1 can be a rowvector // (only really the Hessian when we narrow Y to Z) real matrix boottest::HessianFixedkappa(real rowvector ind1, real scalar ind2, real scalar kappa, real scalar _jk) { real matrix retval; real scalar i if (cols(ind1) > 1) { retval = J(cols(ind1),cols(v),0) for (i=cols(ind1);i;i--) retval[i,] = _HessianFixedkappa(ind1[i], ind2, kappa, _jk) return(retval) } return(_HessianFixedkappa(ind1, ind2, kappa, _jk)) } real rowvector boottest::_HessianFixedkappa(real scalar ind1, real scalar ind2, real scalar kappa, real scalar _jk) { real matrix retval, _retval; pointer (real colvector) scalar pT1L, pT1R if (kappa) { pT1L = pcol(XYbar,ind1+1) // X_∥^' Y_(∥i) if (Repl.Yendog[ind1+1]) pT1L = &(*pT1L :+ SstarUX[ind1+1].M * v) pT1R = ind2? pcol(invXXXZbar,ind2) : &DGP.deltadddot if (Repl.Yendog[ind2+1]) pT1R = &(*pT1R :+ SstarUXinvXX[ind2+1].M * v) // right-side linear term retval = colsum(*pT1L :* *pT1R) // multiply in the left-side linear term } if (kappa != 1) { _retval = YbarYbar[ind1+1,ind2+1] :+ (SstarUYbar[ind2+1].M[,ind1+1] + SstarUYbar[ind1+1].M[,ind2+1]) ' v + colsum(v :* (ind1 <= ind2? SstarUU[ind2+1, ind1+1].M : SstarUU[ind1+1, ind2+1].M) :* v) if (Repl.Yendog[ind1+1] & Repl.Yendog[ind2+1]) _retval = _retval - colsum((SstarinvZperpZperpZperpU[ind1+1].M * v) :* (SstarZperpU[ind2+1].M * v)) if (NFE & FEboot==0) _retval = _retval - colsum(CTstarFEU[ind1+1].M * v :* (invFEwt :* CTstarFEU[ind2+1].M) * v) retval = kappa? kappa*retval:+(1-kappa)*_retval : _retval } if (_jk) { if (kappa) retval[,1] = cross(ind1? *pcol(Repl.XZ, ind1) : *Repl.pXy1par, ind2? *pcol(Repl.V , ind2) : Repl.invXXXy1par) if (kappa != 1) retval[,1] = kappa? kappa * retval[,1] + (1 - kappa) * Repl.YY[ind1+1,ind2+1] : Repl.YY[ind1+1,ind2+1] } return(cols(retval)>1? retval : J(1,cols(v),retval)) // if both vars exogenous, term is same for all b; this duplication is a bit inefficient, but only arises when exog vars involved in null } // Workhorse for WRE CRVE sandwich filling // Given column index ind1 and a matrix betas of all the bootstrap estimates, return all bootstrap realizations of P_X * Z[,ind1]_g ' u\hat_1g^* // for all groups in the intersection of all error clusterings // return value has one row per cap cluster, one col per bootstrap replication pointer(real matrix) scalar boottest::Filling(real scalar ind1, real matrix betas, real scalar _jk) { real scalar i, ind2; real matrix retval, CTstarFEUv, T1, F1, F1beta, F2, SstarUXv, SstarUZperpinvZperpZperp_v; pointer (real matrix) scalar pbetav; pointer (real colvector) pPXYstar; real rowvector _beta, SstarUMZperp_ind2_i; real colvector S pragma unset retval if (granular) { // create pieces of each N x B matrix one at a time rather than whole thing at once retval = J(Clust.N, cols(v), 0) SstarUXv = SstarUX[ind1+1].M * v for (ind2=0; ind2<=Repl.kZ; ind2++) { if (ind2) { _beta = betas[ind2,] pbetav = &(v :* _beta) } else pbetav = &v if (Repl.Yendog[ind2+1]) SstarUZperpinvZperpZperp_v = SstarinvZperpZperpZperpU[ind2+1].M * *pbetav if (NFE & FEboot==0) CTstarFEUv = (invFEwt :* CTstarFEU[ind2+1].M) * *pbetav for (i=Clust.N;i;i--) { S = (*pinfoCapData)[i,]' pPXYstar = &PXZbar[|S,(ind1\ind1)|] if (Repl.Yendog[ind1+1]) pPXYstar = &(*pPXYstar :+ *pXS(Repl.XinvXX,S) * SstarUXv) if (Repl.Yendog[ind2+1]) { SstarUMZperp_ind2_i = *pXS(*Repl.pZperp,S) * SstarUZperpinvZperpZperp_v - (ind2? (*pU2parddot)[|S, (ind2\ind2)|] :* (*pbetav)[*pXS(*pIDBootData,S),] : DGP.u1dddot [|S |] :* v [*pXS(*pIDBootData,S),]) if (NFE & FEboot==0) SstarUMZperp_ind2_i = SstarUMZperp_ind2_i + CTstarFEUv[(*pFEID)[|S|],] // CT_(*,FE) (U ̈_(∥j) ) (S_FE S_FE^' )^(-1) S_FE if (ind2) retval[i,] = retval[i,] - colsum(*pPXYstar :* (Zbar[|S,(ind2\ind2)|] * _beta :- SstarUMZperp_ind2_i)) else retval[i,] = colsum(*pPXYstar :* (DGP.y1bar[|S|] :- SstarUMZperp_ind2_i)) } else retval[i,] = retval[i,] - colsum(*pPXYstar :* (Zbar[|S,(ind2\ind2)|] * _beta)) } } } else if (jk) // coarse error clustering with O(N) operations for (ind2=0; ind2<=Repl.kZ; ind2++) { if (ind2) { _beta = -betas[ind2,] pbetav = &(v :* _beta) } else pbetav = &v T1 = Repl.Yendog[ind1+1]? ScapXYbar[ind2+1].M ' SstarUXinvXX[ind1+1].M : J(0,0,0) // S_∩ (Ybar_(∥j)) (X_∥^' X_∥ )^(-1) [S_* (U ̈_(∥i):*X_∥ )]^' if (Repl.Yendog[ind2+1]) { // add CT_(∩,*) (P_(X_∥ ) Y_(∥i):*U ̈_(∥j) ) if (NClustVar == NBootClustVar & !subcluster) // simple case of one clustering: full crosstab is diagonal if (cols(T1)) _diag(T1, diagonal(T1) + SstarUXinvXX[ind2+1].M ' (*pcol(XYbar,ind1+1))) else T1 = diag(SstarUXinvXX[ind2+1].M ' (*pcol(XYbar,ind1+1))) else { if (Repl.Yendog[ind1+1]==0) T1 = JNcapNstar for (i=Nstar;i;i--) T1[IDCTCapstar[i].M, i] = T1[IDCTCapstar[i].M, i] + SCTcapUXinvXX[ind2+1,i].M * *pcol(XYbar,ind1+1) } if (cols(*Repl.pZperp)) // subtract S_∩ (P_(X_∥ ) Y_(∥i):*Z_⊥ ) (Z_⊥^' Z_⊥ )^(-1) [S_* (U ̈_(∥j):*Z_⊥ )]^' T1 = T1 :- ScapPXYbarZperp[ind1+1].M * SstarinvZperpZperpZperpU[ind2+1].M if (NFE & FEboot==0) T1 = T1 :- CT_FEcapYbar[ind1+1].M ' CTstarFEU[ind2+1].M } if (ind2) { retval = retval + FillingT0[max((ind1,ind2))+1,min((ind1,ind2))+1].M * _beta if (cols(T1)) retval = retval + T1 * *pbetav // - x*beta components } else retval = FillingT0[ind1+1,1].M :+ T1 * v // y component if (Repl.Yendog[ind1+1] & Repl.Yendog[ind2+1]) for (i=Clust.N;i;i--) { S = (*pinfoCapData)[i,]', (.\.) retval[i,] = retval[i,] - colsum(v :* cross(SstarUPX[ind1+1].M[|S|], SstarUMZperp[ind2+1].M[|S|]) * *pbetav) } } else { // coarse error clustering without O(N) operations F1 = invXXXX1par[,ind1] + PiddotRparY[,ind1]; if (Repl.Yendog[ind1+1]) F1 = F1 :+ SstarUXinvXX[ind1+1].M * v retval = J(Clust.N, cols(v), 0) for (i=Clust.N;i;i--) retval[i,] = colsum(F1 :* (ScapXYbar.M[,i] :- negSstarUMZperpX[1,i].M * v)) for (ind2=Repl.kZ; ind2; ind2--) { F1beta = cols(F1)==1? F1 * betas[ind2,] : F1 :* betas[ind2,] for (i=Clust.N;i;i--) { F2 = ScapXYbar[ind2+1].M[,i]; if (Repl.Yendog[ind2+1]) F2 = F2 :- negSstarUMZperpX[ind2+1,i].M * v retval[i,] = retval[i,] - colsum(F1beta :* F2) } } } if (_jk) retval[,1] = *_panelsum(*pcol(Repl.PXZ, ind1), *Repl.py1par - Repl.Z * betas[,1], *pinfoCapData) return(&retval) } void boottest::MakeWREStats(real scalar w) { real scalar c, b, i, _jk real colvector numer_b real rowvector numerw, val, YY11, YY12, YY22, YPXY11, YPXY12, YPXY22, x11, x12, x21, x22, kappas, YY12YPXY12 real matrix deltanumer, Jcap, J_b, Jcaps, vec struct smatrix rowvector A pragma unset vec; pragma unset val _jk = jk & w==1 if (Repl.kZ == 1) { // optimized code for 1 coefficient in bootstrap regression if (LIML) { YY11 = HessianFixedkappa(0, 0, 0, _jk) // kappa=0 => Y*MZperp*Y YY12 = HessianFixedkappa(0, 1, 0, _jk) YY22 = HessianFixedkappa(1, 1, 0, _jk) YPXY11 = HessianFixedkappa(0, 0, 1, _jk) // kappa=1 => Y*PXpar*Y YPXY12 = HessianFixedkappa(0, 1, 1, _jk) YPXY22 = HessianFixedkappa(1, 1, 1, _jk) YY12YPXY12 = YY12 :* YPXY12 x11 = YY22 :* YPXY11 - YY12YPXY12 // elements of YYstar^-1 * YPXYstar up to factor of det(YYstar) x12 = YY22 :* YPXY12 - YY12 :* YPXY22 x21 = YY11 :* YPXY12 - YY12 :* YPXY11 x22 = YY11 :* YPXY22 - YY12YPXY12 kappas = .5 * (x11 + x22); kappas = 1 :/ (1 :- (kappas - sqrt(kappas:*kappas - x11:*x22 + x12:*x21)) :/ (YY11 :* YY22 - YY12 :* YY12)) // solve quadratic equation for smaller eignenvalue; last term is det(YYstar) if (Fuller) kappas = kappas :- Fuller / (_Nobs - kX) betas = (kappas :* (YPXY12 - YY12) + YY12) :/ (As = kappas :* (YPXY22 - YY22) + YY22) } else { As = HessianFixedkappa(1, 1, kappa, _jk) betas = HessianFixedkappa(1, 0, kappa, _jk) :/ As } if (null) numerw = betas :+ (Repl.Rt1 - *pr) / Repl.RRpar else { numerw = betas :- betas[1] if (w==1) numerw[1] = betas[1] + (Repl.Rt1 - *pr) / Repl.RRpar } storeWtGrpResults(pnumer, w, numerw) if (bootstrapt) { if (robust) { Jcaps = *Filling(1, betas, _jk) :/ As for (c=1; c<=NErrClustCombs; c++) { // sum sandwich over error clusterings if (NClustVar != 1 & rows(Clust[c].order)) Jcaps = Jcaps[Clust[c].order,] Jcap = *_panelsum(Jcaps, Clust[c].info) _clustAccum(denom.M, c, colsum(Jcap:*Jcap)) } } else denom.M = (HessianFixedkappa(0,0,0, _jk) - 2 * betas :* HessianFixedkappa(0, 1, 0, _jk) + betas:*betas :* HessianFixedkappa(1, 1, 0, _jk)) :/ As // classical error variance storeWtGrpResults(pDist, w, sqrt? numerw:/sqrt(denom.M) : numerw :* numerw :/ denom.M) denom.M = Repl.RRpar * Repl.RRpar * denom.M[1] // not a bug } } else { // WRE for >1 coefficient in bootstrap regression betas = J(Repl.kZ, cols(v), 0) A = smatrix(cols(v)) if (LIML) { for (i=Repl.kZ;i>=0;i--) { YYstar [i+1].M = HessianFixedkappa(0..i , i, 0, _jk) // kappa=0 => Y*MZperp*Y YPXYstar[i+1].M = HessianFixedkappa(i..Repl.kZ, i, 1, _jk) // kappa=1 => Y*PXpar*Y } for (b=cols(v); b; b--) { for (i=Repl.kZ;i>=0;i--) { YYstar_b [|. ,i+1\ i+1,i+1|] = YYstar [i+1].M[,b] // fill uppper triangle, which is all that invsym() looks at YPXYstar_b[|i+1,i+1\Repl.kZ+1,i+1|] = YPXYstar[i+1].M[,b] // fill lower triangle to prepare for _makesymmetric() } _makesymmetric(YPXYstar_b) eigensystemselecti(invsym(YYstar_b) * YPXYstar_b, Repl.kZ+1\Repl.kZ+1, vec, val) kappa = 1/(1 - Re(val)) // sometimes a tiny imaginary component sneaks into val if (Fuller) kappa = kappa - Fuller / (_Nobs - kX) betas[,b] = (A[b].M = invsym(kappa*YPXYstar_b[|2,2\.,.|] + (1-kappa)*YYstar_b[|2,2\.,.|])) * (kappa*YPXYstar_b[|2,1\.,1|] + (1-kappa)*YYstar_b[|1,2\1,.|]') } } else { deltanumer = HessianFixedkappa(1..Repl.kZ, 0, kappa, _jk) for (i=Repl.kZ;i;i--) deltadenom[i].M = HessianFixedkappa(1..i, i, kappa, _jk) for (b=cols(v); b; b--) { for (i=Repl.kZ;i;i--) deltadenom_b[|.,i\i,i|] = deltadenom[i].M[,b] // fill uppper triangle, which is all that invsym() looks at betas[,b] = (A[b].M = invsym(deltadenom_b)) * deltanumer[,b] } } if (robust) for(i=Repl.kZ;i;i--) Zyi[i].M = *Filling(i, betas, _jk) else for (i=Repl.kZ;i>=0;i--) YYstar[i+1].M = HessianFixedkappa(i..Repl.kZ, i, 0, _jk) // kappa=0 => Y*MZperp*Y for (b=cols(v); b; b--) { numer_b = null | w==1 & b==1? (Repl.RRpar * betas[,b] + Repl.Rt1) - *pr : Repl.RRpar * (betas[,b] - betas[,1]) if (bootstrapt) { if (robust) { // Compute denominator for this WRE test stat for(i=Repl.kZ;i;i--) _Jcap[,i] = Zyi[i].M[,b] Jcap = _Jcap * (A[b].M * Repl.RRpar') for (c=1; c<=NErrClustCombs; c++) { if (NClustVar != 1 & rows(Clust[c].order)) Jcap = Jcap[Clust[c].order,] J_b = *_panelsum(Jcap, Clust[c].info) _clustAccum(denom.M, c, cross(J_b,J_b)) } } else { // non-robust for (i=Repl.kZ;i>=0;i--) YYstar_b[|i+1,i+1\Repl.kZ+1,i+1|] = YYstar[i+1].M[,b] // fill lower triangle for makesymmetric() denom.M = (Repl.RRpar * A[b].M * Repl.RRpar') * ((-1 \ betas[,b]) ' makesymmetric(YYstar_b) * (-1 \ betas[,b])) // 2nd half is sig2 of errors } (*pDist)[b+WeightGrpStart[w]-1] = sqrt? numer_b/sqrt(denom.M) : cross(numer_b, invsym(denom.M) * numer_b) // hand-code for 2-dimensional? } (*pnumer)[,b+WeightGrpStart[w]-1] = numer_b // slight inefficiency: in usual bootstrap-t case, only need to save numerators in numer if getdist("numer") is coming because of svmat(numer) } } if (w==1 & bootstrapt) statDenom = denom.M // original-sample denominator } // // WCR/WCU // // For non-WRE, construct stuff that depends linearly or quadratically on r, possibly by interpolation void boottest::MakeInterpolables() { real scalar h1, h2, d1, d2, c; real matrix tmp; real colvector Delta, newPole if (interpolable) { if (rows(anchor)==0) { // first call? save current r as permanent anchor for interpolation _MakeInterpolables(anchor = *pr) numer0 = *pnumer if (interpolate_u) ustar0 = *pustar if (robust) Jcd0 = *pJcd return } if (rows(poles)) // been here at least twice? interpolate unless current r stretches range > 2X in some dimension(s) newPole = abs(*pr - anchor) :> 2 * abs(poles) else { // second call: from anchor make set of orthogonal poles, which equal anchor except in one dimension poles = *pr - anchor if (robust) // grab quadratic denominator from *previous* (1st) evaluation denom0 = denom newPole = J(q,1,1) // all poles new } if (any(newPole)) { // prep interpolation for (h1=1;h1<=q;h1++) if (newPole[h1]) { poles[h1] = (*pr)[h1] - anchor[h1] (tmp = anchor)[h1] = (*pr)[h1] // if q>1 this creates anchor points that are not graphed, an inefficiency. But simpler to make the deviations from 1st point orthogonal _MakeInterpolables(tmp) // calculate linear stuff at new anchor dnumerdr[h1].M = (*pnumer - numer0) / poles[h1] if (interpolate_u) dudr[h1].M = (*pustar - ustar0) / poles[h1] if (robust) // df > 1 for an ARubin test with >1 instruments. for (d1=1;d1<=df;d1++) { for (c=1;c<=NErrClustCombs;c++) { dJcddr[h1].M[c,d1].M = ((*pJcd)[c,d1].M - Jcd0[c,d1].M) / poles[h1] for (d2=1;d2<=d1;d2++) { tmp = colsum(Jcd0[c,d1].M :* dJcddr[h1].M[c,d2].M) if (d1 != d2) tmp = tmp + colsum(Jcd0[c,d2].M :* dJcddr[h1].M[c,d1].M) // for diagonal items, faster to just double after the c loop _clustAccum(ddenomdr[h1].M[d1,d2].M, c, tmp) } } ddenomdr[h1].M[d1,d1].M = ddenomdr[h1].M[d1,d1].M + ddenomdr[h1].M[d1,d1].M // double diagonal terms } } if (robust) // quadratic interaction terms for (h1=1;h1<=q;h1++) for (h2=h1;h2;h2--) if (newPole[h1] | newPole[h2]) for (d1=df;d1;d1--) for (d2=d1;d2;d2--) for (c=1;c<=NErrClustCombs;c++) _clustAccum(ddenomdr2[h1,h2].M[d1,d2].M, c, colsum(dJcddr[h1].M[c,d1].M :* dJcddr[h2].M[c,d2].M)) Delta = poles interpolating = 1 if (q==2) { // in this case we haven't yet actually computed interpolables at *pr, so interpolate them numerw = numer0 + dnumerdr.M * Delta[1] + dnumerdr[2].M * Delta[2] if (interpolate_u) { pustar = &(ustar0 + dudr.M * Delta[1] + dudr [2].M * Delta[2]) } } } else { // routine linear interpolation if the anchors not moved Delta = *pr - anchor numerw = numer0 + dnumerdr.M * Delta[1]; if (q > 1) numerw = numerw + dnumerdr[2].M * Delta[2] if (interpolate_u) { pustar = &(ustar0 + dudr.M * Delta[1]); if (q > 1) pustar = &(*pustar + dudr[2].M * Delta[2]) } } if (robust) // even if an anchor was just moved, and linear components just computed from scratch, do the quadratic interpolation now, from the updated linear factors if (q==1) for (d1=df;d1;d1--) for (d2=d1;d2;d2--) denom[d1,d2].M = denom0[d1,d2].M + ddenomdr.M[d1,d2].M * Delta + ddenomdr2.M[d1,d2].M * (Delta * Delta) else // q==2 for (d1=df;d1;d1--) for (d2=d1;d2;d2--) denom[d1,d2].M = denom0[d1,d2].M + ddenomdr[1].M[d1,d2].M * Delta[1] + ddenomdr[2].M[d1,d2].M * Delta[2] + ddenomdr2[1,1].M[d1,d2].M * (Delta[1] * Delta[1]) + ddenomdr2[2,1].M[d1,d2].M * (Delta[1] * Delta[2]) + ddenomdr2[2,2].M[d1,d2].M * (Delta[2] * Delta[2]) } else // non-interpolable cases _MakeInterpolables(*pr) } // Construct stuff that depends linearly or quadratically on r and doesn't depend on v. No interpolation. void boottest::_MakeInterpolables(real colvector r) { real scalar d, c, _jk; real colvector uXAR; pointer (real matrix) scalar pustarXAR, ptmp if (ML) uXAR = *pSc * (AR = *pA * *pR') else { if (ARubin) DGP.Estimate(jk, r) else if (kappa) { if (null) { // in score bootstrap for IV/GMM, if imposing null, then DGP constraints, kappa, Hessian, etc. do vary with r and must be set now DGP.Estimate(jk, *pr1 \ r) DGP.InitTestDenoms() } } else // regular OLS DGP.Estimate(jk, null? *pr1 \ r : *pr1) DGP.MakeResiduals(jk) } for (_jk=jk; _jk>=0; _jk--) { // for jackknife, run 2nd time on full original sample to get test stat if (!ML & (scoreBS | (robust & granular * !(jk & !_jk) < NErrClustCombs))) uXAR = DGP.u1ddot[1+_jk].M :* pM->XAR SuwtXA = scoreBS? (B? (NClustVar? *_panelsum(uXAR, infoBootData) : uXAR ) : colsum(uXAR)') : *DGP.pA * *_panelsum2(*pX1, *pX2, DGP.u1ddot[1+_jk].M, infoBootData)' // same calc as in score BS but broken apart to grab intermediate stuff, and assuming residuals defined; X2 empty except in Anderson-Rubin if (robust & granular < NErrClustCombs & bootstrapt) { pustarXAR = _panelsum(uXAR, *pinfoAllData) // collapse data to all-boot & error-cluster-var intersections. If no collapsing needed, _panelsum() will still fold in any weights if (B) { if (scoreBS) for (d=df;d;d--) Kd[d].M = JNcapNstar // inefficient, but not optimizing for the score bootstrap else for (d=df;d;d--) Kd[d].M = *_panelsum2(*pX1, *pX2, *pcol(DGP.XAR,d), *pinfoCapData) * SuwtXA // final term in (64), for c=intersection of all error clusters if (NFE & FEboot==0) CT_WE = crosstabFE(DGP.u1ddot[1+_jk].M, infoBootData) for (d=df;d;d--) { // subtract crosstab of u:*XAR wrt bootstrapping cluster combo and all-cluster-var intersections crosstabCapstarMinus(Kd[d].M, *pcol(*pustarXAR,d)) if (NFE & FEboot==0) Kd[d].M = Kd[d].M + pM->CT_XAR[d].M ' (invFEwt :* CT_WE) // middle term of (64) if (scoreBS) Kd[d].M = Kd[d].M - ClustShare * colsum(Kd[d].M) // recenter } for (c=1+granular; c<=NErrClustCombs; c++) { if (rows(Clust[c].order)) for (d=df;d;d--) Kd[d].M = Kd[d].M[Clust[c].order,] for (d=df;d;d--) Kcd[c,d].M = *_panelsum(Kd[d].M, Clust[c].info) } } else { // B = 0. In this case, only 1st term of (64) is non-zero after multiplying by v* (= all 1's), and it is then a one-way sum by c if (scoreBS) pustarXAR = &(*pustarXAR :- ClustShare * colsum(*pustarXAR)) // recenter if OLS for (c=1; c<=NErrClustCombs; c++) { if (rows(Clust[c].order)) pustarXAR = &((*pustarXAR)[Clust[c].order,]) ptmp = _panelsum(*pustarXAR, Clust[c].info) for (d=df;d;d--) Kcd[c,d].M = *pcol(*ptmp,d) } } } MakeNumerAndJ(1, _jk, r) // compute J = K * v; if Nw > 1, then this is for 1st group; if interpolating, it is only group, and may be needed now to prep interpolation } } // compute stuff depending linearly on v, needed to prep for interpolation // _jk = 0 for non-jackknife, or for non-jk, full original-sample estimation of test stat for jk void boottest::MakeNumerAndJ(real scalar w, real scalar _jk, | real colvector r) { // called to *prepare* interpolation, or when w>1, in which case there is no interpolation real scalar c, d; real matrix _v; real matrix betadev if (jk & !_jk) { if (!ARubin & null) (*pnumer)[,1] = v_sd * ( // full original-sample test stat numerator scoreBS? (B? colsum(SuwtXA)' : SuwtXA ) : (robust==0 | granular | purerobust ? *pR * (betadev = rowsum(SuwtXA)) : rowsum(*pR * SuwtXA))) } else { numerw = scoreBS? (B? cross(SuwtXA, v) : SuwtXA * v_sd ) : (robust==0 | granular | purerobust? *pR * (betadev = SuwtXA * v) : (*pR * SuwtXA) * v) if (w==1) { if ( ARubin) numerw[,1] = v_sd * DGP.Rpar * DGP.beta.M[|kX1+1\.|] // coefficients on excluded instruments in ARubin OLS else if (null==0) numerw[,1] = v_sd * (*pR * (ML? beta : pM->Rpar * pM->beta.M) - r) // Analytical Wald numerator; if imposing null then numer[,1] already equals this. If not, then it's 0 before this. } storeWtGrpResults(pnumer, w, numerw) } if (interpolate_u) { pustar = B? &(v :* DGP.u1ddot[1+_jk].M) : &DGP.u1ddot[1+_jk].M if (scoreBS) pustar = &(*pustar :- colsum(*pustar) * ClustShare) // Center variance if interpolated else pustar = &(*pustar - *pX12B(*pX1, *pX2, betadev)) // residuals of wild bootstrap regression are the wildized residuals after partialling out X (or XS) (Kline & Santos eq (11)) } if (B & robust & bootstrapt) if (jk & !_jk) // fill first col of K's with values for original, full-sample stat for (c=NErrClustCombs; c; c--) for (d=df;d;d--) ((*pJcd)[c,d].M)[,1] = rowsum(Kcd[c,d].M) * v_sd else { if (granular | purerobust) // optimized treatment when bootstrapping by many/small groups if (purerobust & !interpolable) pustar = &(*partialFE(&(DGP.u1ddot[1+_jk].M :* v)) - *pX12B(*pX1, *pX2, betadev)) // but avoid this idiosyncratic calculation and do everything through Jcd when interpolating, for clean code else { // clusters small but not all singletons if (NFE & FEboot==0) { pustar = partialFE(&(DGP.u1ddot.M :* v[*pIDBootData,])) for (d=df;d;d--) (*pJcd)[1,d].M = *_panelsum(*pustar, pM->WXAR[d].M, *pinfoCapData) - *_panelsum2(*pX1, *pX2, pM->WXAR[d].M, *pinfoCapData) * betadev } else { _v = v[*pIDBootAll,] for (d=df;d;d--) (*pJcd)[1,d].M = *_panelsum(*_panelsum(DGP.u1ddot[1+_jk].M, pM->WXAR[d].M, *pinfoAllData) :* _v, infoErrAll) - *_panelsum2(*pX1, *pX2, pM->WXAR[d].M, *pinfoCapData) * betadev } } for (c=NErrClustCombs; c>granular; c--) for (d=df;d;d--) (*pJcd)[c,d].M = Kcd[c,d].M * v } } void boottest::MakeNonWREStats(real scalar w) { real scalar i, c, j, k; real matrix ustar2, tmp, invdenom; real colvector numer_k, ustar_k; pointer (real matrix) scalar pAR; real rowvector t1, t2, t12 if (w > 1) MakeNumerAndJ(w, jk) if (bootstrapt == 0) return if (robust) { if (interpolating==0) { // these quadratic computations needed to *prepare* for interpolation but are superseded by interpolation once it is going if (purerobust & !interpolable) ustar2 = *pustar :* *pustar for (i=df;i;i--) for (j=i;j;j--) { if (c = purerobust & !interpolable) _clustAccum(denom[i,j].M, c, cross(pM->WXAR[i].M, pM->WXAR[j].M, ustar2)) // c=purerobust not a bug for (c++; c<=NErrClustCombs; c++) _clustAccum(denom[i,j].M, c, colsum((*pJcd)[c,i].M :* (*pJcd)[c,j].M)) // (60) } } if (df == 1) { storeWtGrpResults(pDist, w, numerw :/ sqrt(denom.M)) if (w==1) statDenom = denom.M[1] // original-sample denominator } else if (df == 2) { // hand-code 2D numer'inv(denom)*numer t1 = numerw[1,]; t2 = numerw[2,]; t12 = t1:*t2 storeWtGrpResults(pDist, w, (t1:*t1:*denom[2,2].M - (t12+t12):*denom[2,1].M + t2:*t2:*denom[1,1].M) :/ (denom[1,1].M:*denom[2,2].M - denom[2,1].M:*denom[2,1].M)) if (w==1) statDenom = denom[1,1].M[1], denom[2,1].M[1] \ denom[2,1].M[1], denom[2,2].M[1] // original-sample denominator } else { // build each replication's denominator from vectors that hold values for each position in denominator, all replications tmp = J(df,df,0) for (k=cols(v); k; k--) { for (i=df;i;i--) for (j=i;j;j--) tmp[j,i] = denom[i,j].M[k] // fill upper triangle, which is all invsym() looks at numer_k = numerw[,k] (*pDist)[k+WeightGrpStart[w]-1] = numer_k ' invsym(tmp) * numer_k // in degenerate cases, cross() would turn cross(.,.) into 0 } if (w==1) statDenom = makesymmetric(tmp') // original-sample denominator } } else { // non-robust pAR = ML? &AR : &pM->AR if (df == 1) { // optimize for one null constraint denom.M = *pR * *pAR denom.M = ML? denom.M * (v_sd * v_sd) : denom.M :* colsum(*pustar :* *pustar) storeWtGrpResults(pDist, w, numerw :/ sqrt(denom.M)) if (w==1) statDenom = denom.M[1] // original-sample denominator } else { denom.M = *pR * *pAR if (ML) { denom.M = denom.M * (v_sd * v_sd) for (k=cols(v); k; k--) { numer_k = numerw[,k] (*pDist)[k+WeightGrpStart[w]-1] = cross(numer_k, invsym(denom.M) * numer_k) } if (w==1) statDenom = denom.M // original-sample denominator } else { invdenom = invsym(denom.M) for (k=cols(v); k; k--) { numer_k = numerw[,k] ustar_k = (*pustar)[,k] (*pDist)[k+WeightGrpStart[w]-1] = cross(numer_k, invdenom * numer_k) / (tmp = cross(ustar_k, ustar_k)) } if (w==1) statDenom = denom.M * tmp // original-sample denominator } } } } // like panelsetup() but can group on multiple columns, like sort(), and faster. But doesn't take minobs, maxobs arguments. // Does take optional third argument, a matrix in which to store standardized ID variable, starting from 1 real matrix _panelsetup(real matrix X, real rowvector cols, | pointer(real colvector) scalar pID) { real matrix info; real scalar i, N; real scalar p; real rowvector tmp, id N = rows(X) info = J(N, 2, N); if (args()>2) pID = &J(N, 1, 1) info[1,1] = p = 1 id = X[1, cols] for (i=2; i<=N; i++) { if ((tmp=X[i,cols]) != id) { info[ p,2] = i - 1 info[++p,1] = i id = tmp } if (args()>2) (*pID)[i] = p } return (info[|.,.\p,.|]) } // Do panelsum() except that a single missing value in X doesn't make all results missing and // efficiently handles case when all groups have one row. pointer(real matrix) scalar _panelsum(real matrix X, real matrix arg2, | real matrix arg3) { if (args()==2) { if (rows(arg2)==0 | rows(arg2)==rows(X)) return(&X) } else if (rows(arg3)==0 | rows(arg3)==rows(X)) return(arg2==1? &X : &(X :* arg2)) // if no collapsing called for, still fold in provided weights return (cols(arg3)? &panelsum(X, arg2, arg3) : &panelsum(X, arg2)) } // concatenation of two _panelsum's pointer(real matrix) scalar _panelsum2(real matrix X1, real matrix X2, real matrix arg2, | real matrix arg3) return(args()==2? &(*_panelsum(X1,arg2),*_panelsum(X2,arg2)) : &(*_panelsum(X1,arg2,arg3),*_panelsum(X2,arg2,arg3))) // given a vector, return indices of the non-zero elements, like selectindex() function added in Stata 13 // if v = 0 (so can't tell if row or col vector), returns J(1, 0, 0) real vector boottest::_selectindex(real vector v) { real scalar rows if (v==0) return(J(1,0,0)) rows = rows(v) return(select(rows>1? 1::rows : 1..cols(v), v)) } // Return matrix that counts from 0 to 2^N-1 in binary, one column for each number, one row for each binary digit // except use provided lo and hi values for 0 and 1 real matrix boottest::count_binary(real scalar N, real scalar lo, real scalar hi) { real matrix tmp if (N<=1) return (lo , hi) tmp = count_binary(N-1, lo, hi) return (J(1, cols(tmp), lo), J(1, cols(tmp), hi) \ tmp, tmp) } // cross-tab sum of a column vector w.r.t. given panel info and fixed-effect var // one row per FE, one col per other grouping real matrix boottest::crosstabFE(real colvector v, real matrix info) { real matrix retval; real scalar i, j, tmp; real colvector _FEID, _v, vw retval = J(NFE, rows(info), 0) if (haswt) { vw = v :* sqrtwt if (cols(info) & rows(info)!=rows(v)) for (i=cols(retval);i;i--) { _FEID = panelsubmatrix(*pFEID, i, info) _v = panelsubmatrix(vw , i, info) for (j=rows(_FEID);j;j--) { tmp = _FEID[j] retval[tmp,i] = retval[tmp,i] + _v[j] } } else // "robust" case, no clustering for (i=cols(retval);i;i--) retval[(*pFEID)[i],i] = vw[i] } else if (cols(info) & rows(info)!=rows(v)) for (i=cols(retval);i;i--) { _FEID = panelsubmatrix(*pFEID, i, info) _v = panelsubmatrix(v , i, info) for (j=rows(_FEID);j;j--) { tmp = _FEID[j] retval[tmp,i] = retval[tmp,i] + _v[j] } } else // "robust" case, no clustering for (i=cols(retval);i;i--) retval[(*pFEID)[i],i] = v[i] return(retval) } // Multiply an M x N x L array, stored as an L-vector of M x N matrices, by an L-vector, returning M x N result real matrix boottest::threebyone(struct smatrix vector A, real colvector b) { real scalar i; real matrix retval if (length(b)) { retval = A.M * b[1] for (i=length(A);i>1;i--) retval = retval + A[i].M * b[i] return(retval) } else return(0) // can't determine dimensinons in this case } // subtract crosstab of v wrt bootstrapping cluster and all-cluster-var intersections from M // M should have one row for each all-cluster-var (including bootstrap cluster) intersection and one col for each bootstrap cluster // *** v needs to have been panelsum'd with pinfoAllData void boottest::crosstabCapstarMinus(real matrix M, real colvector v) { real colvector tmp; real scalar i if (subcluster) // crosstab c,c* is wide for (i=Clust.N;i;i--) { tmp = infoErrAll[i,]' M[|(i\i), tmp|] = M[|(i\i), tmp|] - v[|tmp|]' } else if (NClustVar == NBootClustVar) // crosstab c,c* is square _diag(M, diagonal(M) - v) else // crosstab c,c* is tall for (i=Nstar;i;i--) { tmp = Clust[BootClust].info[i,]' M[|tmp, (i\i)|] = M[|tmp, (i\i)|] - v[|tmp|] } } // given a pre-configured boottest linear model with one-degree null imposed, compute p value for given value of r // used with optimize() to construct confidence intervals // performs no error checking real scalar boottest::r_to_p(real colvector r) { pr = &r setdirty(1, 1) // set dirty = 1, but leave initialized=0, which we want when only changing r return (getpadj()) } // Chandrupatla 1997, "A new hybrid quadratic-bisection algorithm for finding the zero of a nonlinear function without using derivatives" // x1, x2 must bracket the true value, with f1=f(x1) and f2=f(x2). alpha is target p value. real scalar boottest::search(real scalar alpha, real scalar f1, real scalar x1, real scalar f2, real scalar x2) { real scalar x, t, fx, phi1, phi1_2, xi1, x3, f3 t = 0.5 do { fx = r_to_p(x = x1 + t * (x2 - x1)) if (fx>f1 == fx>f2) // violation of monotonicity because of precision problems? That's as good as it gets. return(x) if (fx xi1 | xi1 > phi1 + phi1 - phi1_2) t = 0.5 else { t = ((f3 - alpha) / (f1 - f2) + (x3 - x1) / ((x2 - x1) * (f3 - f1)) * (f2 - alpha)) * (f1 - alpha) / (f3 - f2) if (t < 0.000001) t = 0.000001 else if (t > 0.999999) t = 0.999999 } } while (1) } // derive wild bootstrap-based CI, for case of linear model with one-degree null imposed // also prepare data to plot confidence curve or surface void boottest::plot() { real scalar tmp, alpha, _quietly, c, d, i, j, halfwidth, p_lo, p_hi, p_confpeak, Phi; real colvector lo, hi; pointer (real colvector) scalar _pr _quietly = quietly; _pr = pr setquietly(1) alpha = 1 - level*.01 _editmissing(gridpoints, 25) boottest() Phi = invnormal(alpha/2) if (ARubin) halfwidth = abs(confpeak * sqrt((abs((*pDist)[1]) * multiplier) / df) / Phi) // abs(confpeak) * invnormal(getpadj(1)/2) / invnormal(alpha/2) else { halfwidth = (-1.5 * Phi) * sqrt(diagonal(getV())) confpeak = getb() + *pr } if (q==2) { // 2D plot lo = hi = J(2, 1, .) for(d=q;d;d--) { lo[d] = editmissing(gridmin[d], confpeak[d] - halfwidth[d]) hi[d] = editmissing(gridmax[d], confpeak[d] + halfwidth[d]) stata("_natscale " + strofreal(lo[d]) + " " + strofreal(hi[d]) + " 4") // using Stata logic for setting graph bounds ensures good-looking contour plot if (gridmin[d]==.) { stata("local min = r(min)") // for some reason st_global("r(min)") doesn't work lo[d] = strtoreal(st_local("min")) } if (gridmax[d]==.) { stata("local max = r(max)") hi[d] = strtoreal(st_local("max")) } } plotX = (rangen(lo[1], hi[1], gridpoints[1]) # J(gridpoints[2],1,1)), (J(gridpoints[1],1,1) # rangen(lo[2], hi[2], gridpoints[2])) plotY = J(rows(plotX), 1, .) } else { // 1D plot if (alpha<=0) alpha = .05 // if level=100, no CI constructed, but we need a reasonable alpha to choose graphing bounds if (alpha > 0 & cols(v)-1 <= 1/alpha-1e6) { setquietly(_quietly) if (quietly==0) errprintf("\nError: need at least %g replications to resolve a %g%% two-sided confidence interval.\n", ceil(1/alpha), level) return } if (gridmin[1]==. | gridmax[1]==.) { if (B) { lo = editmissing(gridmin[1], confpeak - halfwidth) // initial guess based on classical distribution hi = editmissing(gridmax[1], confpeak + halfwidth) } else { tmp = sqrt(statDenom[1,1]) * (small? invttail(df_r, alpha/2) : -invnormal(alpha/2)) lo = editmissing(gridmin[1], confpeak - tmp) hi = editmissing(gridmax[1], confpeak + tmp) if (scoreBS & (null | willplot)==0) { // if doing simple Wald test with no graph, we're done CI = lo, hi return } } if (abs(lo - *pr) > abs(hi - *pr)) { // brute force way to ensure that first trial bound tested is the farther one from *pr, for better interpolation if (gridmin[1]==. & ptype!=2) // unless lower-tailed p value, try at most 10 times to bracket confidence set by symmetrically widening for (i=10; i & -(p_lo=r_to_p(lo)) < -alpha; i--) { tmp = hi - lo lo = lo - tmp if (gridmax[1]==. & twotailed) hi = hi + tmp // maintain rough symmetry unless user specified upper bound } if (gridmax[1]==. & ptype!=3) // ditto for high side for (i=10; i & -(p_hi=r_to_p(hi)) < -alpha; i--) { tmp = hi - lo if (gridmin[1]==. & twotailed) lo = lo - tmp hi = hi + tmp } } else { if (gridmax[1]==. & ptype!=3) // ditto for high side for (i=10; i & -(p_hi=r_to_p(hi)) < -alpha; i--) { tmp = hi - lo if (gridmin[1]==. & twotailed) lo = lo - tmp hi = hi + tmp } if (gridmin[1]==. & ptype!=2) // unless upper-tailed p value, try at most 10 times to bracket confidence set by symmetrically widening for (i=10; i & -(p_lo=r_to_p(lo)) < -alpha; i--) { tmp = hi - lo lo = lo - tmp if (gridmax[1]==. & twotailed) hi = hi + tmp // maintain rough symmetry unless user specified upper bound } } } else { // both grid bounds pre-specified lo = gridmin[1] hi = gridmax[1] } plotX = rangen(lo, hi, gridpoints[1]) plotY = J(rows(plotX), 1, .) plotY[1] = p_lo; plotY[rows(plotX)] = p_hi p_confpeak = WRE? . : (twotailed? 1 : .5) if (confpeak < lo) { // insert original point estimate into grid if (gridmin[1] == .) { plotX = confpeak \ plotX plotY = p_confpeak \ plotY c = 1 } } else if (confpeak > hi) { if (gridmax[1] == .) { plotX = plotX \ confpeak plotY = plotY \ p_confpeak c = gridpoints[1] + 1 } } else { c = floor((confpeak - lo)/(hi - lo)*(gridpoints[1] - 1)) + 2 plotX = plotX[|.\c-1|] \ confpeak \ plotX[|c\.|] plotY = plotY[|.\c-1|] \ p_confpeak \ plotY[|c\.|] } } // end 1D plot prep i = 1 do { if (plotY[i] == .) plotY[i] = r_to_p(plotX[i,]') } while (1 < (i = mod(i-2,rows(plotX))+1)) if (hasmissing(plotY)) CI = . , . else if (q==1 & level<100) { // find CI bounds CI = (plotY :> alpha) :/ (plotY :< .); CI = CI[|2\.|] - CI[|.\rows(plotX)-1|] lo = _selectindex(CI:== 1) hi = _selectindex(CI:==-1) if (rows(lo)==0 & rows(hi)==0) CI = . , . else { if (rows(lo)==0) lo = . else if (rows(hi)==0) hi = . else { if ( lo[1 ] > hi[1 ]) lo = . \ lo // non-rejection ranges that are not within grid range if (-lo[rows(lo)] < -hi[rows(hi)]) hi = hi \ . } CI = lo, hi for (i=rows(lo); i; i--) for (j=2; j; j--) if (CI[i,j]<.) CI[i,j] = search(alpha, plotY[CI[i,j]], plotX[CI[i,j]], plotY[CI[i,j]+1], plotX[CI[i,j]+1]) } } if (c < .) { // now that it's done helping graph look good, remove peak point from returned grid for evenness, for Bayesian sampling purposes peak = plotX[c], plotY[c] if (c==1) { plotX = plotX[|2\.|]; plotY = plotY[|2\.|] } else if (c==gridpoints[1]+1) { plotX = plotX[|.\gridpoints[1]|]; plotY = plotY[|.\gridpoints[1]|] } else { plotX = plotX[|.\c-1|] \ plotX[|c+1\.|]; plotY = plotY[|.\c-1|] \ plotY[|c+1\.|] } } setquietly(_quietly); pr = _pr; dirty = 1 // restore backups notplotted = 0 } // return matrix whose rows are all the subsets of a row of numbers. Nil is at bottom. real matrix boottest::combs(real scalar d) { real matrix retval; real scalar i retval = J(2^d, 0, 0) for (i=d;i;i--) retval = retval, J(2^(d-i),1,1) # (1\0) # J(2^(i-1),1,1) return (retval) } // Like Mata's order() but does a stable sort real colvector boottest::stableorder(real matrix X, real rowvector idx) return (order((X, (1::rows(X))), (idx,cols(X)+1))) // Stata interface void boottest_stata(string scalar statname, string scalar dfname, string scalar dfrname, string scalar pname, string scalar padjname, string scalar ciname, string scalar plotname, string scalar peakname, real scalar level, real scalar ptol, real scalar ML, real scalar LIML, real scalar Fuller, real scalar kappa, real scalar ARubin, real scalar null, real scalar scoreBS, real scalar jk, string scalar auxweighttype, string scalar ptype, string scalar statistic, string scalar madjtype, real scalar NumH0s, string scalar X1names, string scalar Y2names, string scalar Ynames, string scalar bname, string scalar Aname, string scalar X2names, string scalar samplename, string scalar scnames, real scalar robust, string scalar IDnames, real scalar NBootClustVar, real scalar NErrClustVar, string scalar FEname, real scalar NFE, real scalar FEdfadj, string scalar wtname, string scalar obswttype, string scalar R1name, string scalar r1name, string scalar Rname, string scalar rname, real scalar B, string scalar repsname, string scalar repsFeasname, real scalar small, string scalar diststat, string scalar distname, string scalar gridmin, string scalar gridmax, string scalar gridpoints, real scalar MaxMatSize, real scalar quietly, string scalar b0name, string scalar V0name, string scalar vname, string scalar NBootClustname) { real matrix X2, ID, FEID, sc, Y2, X1 real colvector wt, Y class boottest scalar M pragma unset ID; pragma unset wt; pragma unset Y2; pragma unset X1; pragma unset Y; pragma unset X2; pragma unset sc M._st_view(sc, ., scnames, samplename) M._st_view(Y , ., Ynames , samplename) M._st_view(X2, ., X2names, samplename) if (FEname != "" ) FEID = st_data(., FEname , samplename) if (IDnames != "") ID = st_data(., IDnames, samplename) if (wtname != "") wt = st_data(., wtname , samplename) // panelsum() doesn't like views as weights M.setMaxMatSize(MaxMatSize) M.setsc(sc) M.setML(ML) M.setY (Y) M.setX2(X2) M.setobswt(wt, obswttype) M.setID(ID, NBootClustVar, NErrClustVar) M.setFEID(FEID, NFE, FEdfadj) M.setR1(st_matrix(R1name), st_matrix(r1name)) M.setR (st_matrix(Rname ), st_matrix(rname )) M.setnull(null) M.setsmall(small) M.setrobust(robust) M.setjk(jk) M.setscoreBS(scoreBS) M.setauxwttype(auxweighttype) M.setptype(ptype) M.setstattype(statistic) M.setB(B) M.setLIML(LIML) M.setFuller(Fuller) M.setkappa(kappa) M.setARubin(ARubin) M.setgrid(strtoreal(tokens(gridmin)), strtoreal(tokens(gridmax)), strtoreal(tokens(gridpoints))) M.setmadjust(madjtype, NumH0s) M.setlevel(level) M.setptol(ptol) M.setquietly(quietly) M._st_view(Y2, ., Y2names, samplename); M.setY2(Y2) M._st_view(X1, ., X1names, samplename); M.setX1(X1) if (bname != "") M.setbeta(st_matrix(bname)') if (Aname != "") M.setA (st_matrix(Aname) ) M.setwillplot(plotname != "") if (plotname != "" | (level<100 & ciname != "")) { if (plotname != "") st_matrix(plotname, M.getplot()) if (cols(M.peak)) st_matrix(peakname, M.getpeak()) if (level<100 & ciname != "") st_matrix(ciname, M.getCI()) } st_numscalar(statname, M.getstat()) st_numscalar(pname , M.getp ()) st_numscalar(repsname, M.getreps()) st_numscalar(repsFeasname, M.getrepsFeas()) st_numscalar(NBootClustname, M.getNBootClust()) st_numscalar(padjname, M.getpadj()) st_numscalar(dfname , M.getdf ()) st_numscalar(dfrname , M.getdf_r()) st_matrix(b0name, M.getb()') st_matrix(V0name, M.getV()) if (distname != "") st_matrix(distname, M.getdist(diststat)) if (vname != "" & B) st_matrix(vname, M.getv()) M.close() } mata mlib create lboottest, dir("c(sysdir_plus)'l") replace mata mlib add lboottest *(), dir("c(sysdir_plus)'l") mata mlib index end