*! 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 z`k'
// dis "`z'"
qui regress `z', `noconstant'
qui predict double `z`k'', residual // partial out the constant
local instr_partialed "`instr_partialed' `z`k''"
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 z`k'
// dis "`z'"
qui predict double `z`k'', residual equation(#`k') // partial out controls from Z
local instr_partialed "`instr_partialed' `z`k''"
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