*! version 0.0 21sep2021 Liyang Sun, lsun20@mit.edu capture program drop manyweakivtest program define manyweakivtest, eclass sortpreserve version 13 set more off _iv_parse 0' local depvar s(lhs)' local endog s(endog)' local covariates s(exog)' local instr s(inst)' local 0 s(zero)' syntax [if] [in] [aweight fweight], /// [NOConstant] // // * Mark sample (reflects the if/in conditions, and includes only nonmissing observations) // marksample touse // markout touse' by' xq' covariates', strok * Parse the dependent variable tempname h y yhat x xhat * dis "covariates'" qui regress depvar' covariates', noconstant' // partial out controls from Y (if empty, then partial out the constant term) qui predict double y', residual qui regress endog' covariates', noconstant' // partial out controls from X (if empty, then partial out the constant term) qui predict double x', residual if "covariates'" == "" & "noconstant'" != "" { local instr_partialed "instr'" // nothing to partial out } else if "covariates'" == "" & "noconstant'" == "" { local instr_partialed "" local k = 1 foreach z of varlist instr' { tempvar zk' // dis "z'" qui regress z', noconstant' qui predict double zk'', residual // partial out the constant local instr_partialed "instr_partialed' zk''" local k = k' + 1 } } else { qui mvreg instr' = covariates', noconstant' // partial out controls from Z local instr_partialed "" local k = 1 foreach z of varlist instr' { tempvar zk' // dis "z'" qui predict double zk'', residual equation(#k') // partial out controls from Z local instr_partialed "instr_partialed' zk''" local k = k' + 1 } } ** now the dep, endogenous varibale and instruments have controls partialled out ** first-stage regression qui regress x' instr_partialed', nocons // the constant term is already partialled out qui predict double h', hat // leverage Z_i'(Z'Z)^-1 Z_i qui predict double xhat', // predicted value Z\hat{\pi} ** reduced-form regression qui regress y' instr_partialed', nocons // the constant term is already partialled out qui predict double yhat', // predicted value Z\hat{\delta} ** move to mata for matrix calculation mata: Sigma_fun("instr_partialed'","yhat'","y'","xhat'","x'","`h'") // dis "The analytical solution to the jackknife AR test inversion are:" // tempname roots Sigma1 // matrix list r(roots) // matrix list r(Sigma1) // ereturn clear // ereturn matrix r r(roots)' // ereturn matrix S r(Sigma1)' end mata: void Sigma_fun( string scalar Z_name, /// string scalar Yhat_name, /// string scalar Y_name, /// string scalar Xhat_name, /// string scalar X_name, /// string scalar H_name /// ) { Z = st_data(.,Z_name) Yhat = st_data(.,Yhat_name) Y = st_data(.,Y_name) Xhat = st_data(.,Xhat_name) X = st_data(.,X_name) H = st_data(.,H_name) N = rows(Z) K = cols(Z) H_diag = diag(H) M_diag = J(N,1,1)-H ZZ_inv = qrinv(Z'*Z) ZZZ_inv = Z*ZZ_inv XMX = X:*X - X:*Xhat YMY = Y:*Y - Y:*Yhat YMX = Y:*X - Y:*Xhat XMY = X:*Y - X:*Yhat XMYYMX = -YMX-XMY Sigma1_hh = 0 Sigma1_gg_0 = 0 Sigma1_gg_1 = 0 Sigma1_gg_2 = 0 Sigma1_gg_3 = 0 Sigma1_gg_4 = 0 for (i=1; i<=N; i++) { if (mod(i,10000) == 0) { printf("Finished %g observations\n",i) } Zi = Z[i,] Pi = ZZZ_inv * Zi' PP = Pi:^2 PP_off = (PP):/((1-H[i])*M_diag + PP) PP_off[i]=0 Sigma1_hh = Sigma1_hh + XMX[i]*PP_off'*XMX Sigma1_gg_0 = Sigma1_gg_0 + YMY[i]*PP_off'*YMY; Sigma1_gg_1 = Sigma1_gg_1 + 2 * YMY[i]*PP_off'*XMYYMX; Sigma1_gg_2 = Sigma1_gg_2 + 2 * XMX[i]*PP_off'*YMY + XMYYMX[i]*PP_off'*XMYYMX; Sigma1_gg_3 = Sigma1_gg_3 + 2 * XMX[i]*PP_off'*XMYYMX; Sigma1_gg_4 = Sigma1_gg_4 + XMX[i]*PP_off'*XMX; } Sigma1 = (Sigma1_gg_0, Sigma1_gg_1, Sigma1_gg_2, Sigma1_gg_3, Sigma1_gg_4) // Sigma1 XPX = X'*Xhat - X'*H_diag*X //a2 YPY = Y'*Yhat - Y'*H_diag*Y //a0 YPX = Y'*Xhat - Y'*H_diag*X //a1 // TODO: numerator is a quadratic inequality // printf("Ybar' P Ybar=%g, %g, %g, and determinant is %g\n",YPY, YPX, XPX, YPX^2 - 4*YPY*XPX) cnum = polyroots((YPY, -2*YPX, XPX)/sqrt(K)) // normalized by \sqrt(K) // cnum cnum_real = J(1,0,.); cnum_nreal = 0; // count the number of real roots for (i=1; i<=2; i++) { if (isrealvalues(cnum[i]) == 1) { cnum_real = (cnum_real, Re(cnum[i])) cnum_real = sort(cnum_real,1) cnum_nreal++ } } cnum_real = sort(cnum_real',1)' // Solving the quartic inequality where c_j is coef on jth order term crit2 = 1.64^2 // TODO: add in user-specified critical value c0 = 2*crit2*Sigma1_gg_0-YPY^2 c1 = 2*crit2*Sigma1_gg_1+4*YPY*YPX c2 = 2*crit2*Sigma1_gg_2-2*YPY*XPX-4*YPX^2 c3 = 2*crit2*Sigma1_gg_3+4*YPX*XPX c4 = 2*crit2*Sigma1_gg_4-XPX^2 c = polyroots((c0, c1, c2, c3, c4)/K) // normalized by K creal = J(1,0,.); cnreal = 0; // count the number of real roots // printf("Coefficients are %g, %g, %g, %g, %g\n",c0, c1, c2, c3, c4) // c for (i=1; i<=4; i++) { if (isrealvalues(c[i]) == 1) { creal = (creal, Re(c[i])) cnreal++ } creal = sort(creal',1)' st_matrix("r(roots)", creal) } printf("The analytical solution to the jackknife AR test inversion are:\n") if (c4 <0 & XPX >0) { printf("Bounded interval, union of the following intervals\n") if (cnum_nreal == 2) { printf("[%9.0g , %9.0g]\n",cnum_real[1],cnum_real[2]) } // cnum_real // bounded interval from quadratic numerator (can be empty tho) // creal // bounded interval from quartic inequality printf("[%9.0g , %9.0g]\n",creal[1],creal[2]) } if (c4 <0 & XPX <0) { printf("Unbounded interval,union of the following intervals\n") if (cnum_nreal == 2) { printf("[-inf, %g ]U[ %g ,+inf]\n",cnum_real[1],cnum_real[2]) } // unbounded interval from quadratic numerator // creal printf("[%9.0g , %9.0g]\n",creal[1],creal[2]) } if (c4 >0 & XPX >0) { printf("Unbounded interval,union of the following intervals\n") // cnum_real // bounded interval from quadratic numerator if (cnum_nreal == 2) { printf("[%9.0g , %9.0g]\n",cnum_real[1],cnum_real[2]) } // unbounded interval from quartic numerator if (cnreal == 0) { printf("[-inf,+inf]\n") } if (cnreal == 2) { printf("[-inf,%g]U[%g,+inf]\n",creal[1],creal[2]) } if (cnreal == 4) { printf("[-inf,%g]U[%g,%g]U[%g,+inf]\n",creal[1],creal[2],creal[3],creal[4]) } } if (c4 >0 & XPX <0) { printf("Unbounded interval,union of the following intervals\n") if (cnum_nreal == 2) { printf("[-inf, %g ]U[ %g ,+inf]",cnum_real[1],cnum_real[2]) } // unbounded interval from quadratic numerator if (cnreal == 0) { printf("[-inf,+inf]\n") } if (cnreal == 2) { printf("[-inf,%g]U[%g,+inf]\n",creal[1],creal[2]) } if (cnreal == 4) { printf("[-inf,%g]U[%g,%g]U[%g,+inf]\n",creal[1],creal[2],creal[3],creal[4]) } } st_matrix("r(Sigma1)", Sigma1) } end