/*
Changelog
v.1.0.0 -- 1 September 2023
(1) Original release.
v.2.0.0 -- 11 November 2023
(1) Option -endogenous()- is now used to specify the names of the endogenous variables if any of them are not the same as the auxiliary dependent variables and are instead functions of the auxiliary dependent variables.
(2) Original -endogenous- option renamed -auxiliary-; use this to specify the auxiliary dependent variables.
(3) Replay functionality: previous estimation output is replayed after inputting "ivcloglog".
(4) Modified output for improved clarity.
(5) Bug fixes and minor tweaks.
v.2.0.1 -- 17 November 2023
(1) Bug fix for when multiple endogenous variables are supplied with -endogenous()-
*/
*! Version 2.0.1 -- William Liu (刘å¨å»‰) -- 17 November 2023
*! ivcloglog -- Complementary log-log model with endogenous covariates (either continuous or functions of observed continuous covariates), instrumented via the control function approach (i.e., 2SRI)
*!
*! Syntax:
*! ivcloglog binary_outcome_var [controls] [if] [in], VHATname(string) AUXiliary(auxdep_vars = instruments[, NOCONstant]) [ENDOgenous(transformed_endo_vars) NOCONstant order(integer) vce(vcetype) NOGENerate DIFFicult_vce SHOWstages]
program ivcloglog, eclass byable(onecall)
version 11 // FVvars and -_rmcoll- option -touse()- not introduced until v11.0
qui ds // Get original varlist
local orig_varlist `r(varlist)'
if _by() {
local BY `"by `_byvars'`_byrc0':"'
}
`BY' _vce_parserun ivcloglog, noeqlist jkopts(eclass): `0'
if "`s(exit)'" != "" {
ereturn local cmdline `"ivcloglog `0'"'
exit
}
if replay() { // For replaying the results
if `"`e(cmd)'"' != "ivcloglog" {
error 301
}
else if _by() {
error 190
}
else {
ivcloglog_out `0'
exit
}
}
cap noi { // To allow us to drop the cfvars
`BY' Estimate `0'
}
local est_rc = _rc
ereturn local cmdline `"ivcloglog `0'"'
drop_cfvars `0' orig_varlist(`orig_varlist') // Drop powers of vhat if requested with -nogenerate- (if they exist)
if `est_rc' { // To restore the original behavior -- with -capture noisily-, the error message is displayed but not the return code
exit `est_rc' // Note: -exit- with a return code only displays the error code; -error- with a return code also displays the corresponding generic error message
}
end
program define Estimate, eclass byable(recall)
syntax varlist(numeric fv ts) [if] [in] ///
[fweight iweight pweight aweight], ///
[ ///
VHATname(string) ///
noCONstant ///
order(integer 1) ///
vce(string) ///
NOGENerate ///
DIFFicult_vce ///
SHOWstages ///
* ///
]
_get_diopts diopts options, `options' // Note: the second argument must be provided anyway when there are non-display options
if ("`weight'" != "") {
local wgt [`weight'`exp']
}
if ("`vhatname'" == "") {
local vhatname "vhat"
}
else if (strpos("`vhatname'", " ")) {
di as err "option {bf:vhatname()} should only contain one word"
exit 198
}
gettoken depvar exogvars : varlist
_fv_check_depvar `depvar', k(1) // Stops depvar being a factor variable
local exogvars : list uniq exogvars // Check for duplicates
local nocon2 `constant' // For better readability
* Parse -auxiliary(...)- option
parse_remaining_opt, exogvars(`exogvars') `options'
local options `s(options)'
local num_endo = `s(num_endo)'
local num_aux = `s(num_aux)'
local allinst `s(allinst)'
local endovars `s(endovars)'
local auxdepvars `s(auxdepvars)'
local nocon1 `s(nocons1)'
marksample touse
_vce_parse `touse', optlist(Robust) argoptlist(CLuster) : , vce(`vce') // Apart from parsing the -vce()- option, this will also check variables inside of it for missing values
local clustervar "`r(cluster)'"
markout `touse' `endovars' `allinst'
* Check for "perfect predictors" of the outcome and collinearity amongst the regressors; also expand factor variables
// Checking the exogvars once rather than twice would be more efficient. (Recall that they appear as instruments!)
// Must check the second stage first since the number of observations can decrease if "perfect predictors" are present
tempname omit1 omit2
di as text "Checking second stage for perfect predictors of the outcome and regressor collinearity..."
iden_check_and_fvexpand_2, lhs(`depvar') exog(`exogvars') endo(`endovars') touse(`touse') wgt(`wgt') `difficult_vce' `nocon2'
if ("`exogvars'" != "") {
local exogvars `r(exog)'
matrix `omit2' = r(exog_omit)
matrix `omit2' = `omit2', J(1, `=`num_aux'*`order' + `num_endo'', 1) // Add elements for the endovars and residual powers
}
else {
matrix `omit2' = J(1, `=`num_aux'*`order' + `num_endo'', 1)
}
di as text "Checking first stage for regressor collinearity..."
iden_check_and_fvexpand_1, allinst(`allinst') touse(`touse') wgt(`wgt') `difficult_vce' `nocon1'
local allinst `r(allinst)'
matrix `omit1' = r(allinst_omit)
* Initial values
** First stage
tempname F1 F2 F1_forone_eq F
forval i = 1/`num_aux' {
forval p = 1/`order' {
local resid`i'_`p' `vhatname'`i'_`p' // For better readability
qui gen double `resid`i'_`p'' = . // Previously, for some reason, initializing the higher powers with missing values didn't work when not using -from(...)-
local resid`i'_powers = "`resid`i'_powers' `resid`i'_`p''"
}
local endovar`i' : word `i' of `auxdepvars'
init_values_1 `endovar`i'' `allinst' if `touse' `wgt', `nocon1' predictvars(`resid`i'_powers') order(`order') `showstages' // Store powers of vhat for auxdep variable `i' in local macro `resid`i'_powers'
matrix `F1_forone_eq' = r(param_est)
matrix `F1' = nullmat(`F1'), `F1_forone_eq' // -nullmat()- relaxes the restriction that `A' must exist
local yparam1 "`r(yparam)'"
local yparam1_combined "`yparam1_combined' `yparam1'"
local resid_combo "`resid_combo' `resid`i'_1'" // vhats from all endo. variables
forval p = 1/`order' {
local resid_powers_combo "`resid_powers_combo' `resid`i'_`p''" // Powers of the vhats from all endo. variables
}
local inst_opt1 "`inst_opt1' instruments(`endovar`i'':`allinst', `nocon1')"
}
** Second stage
local 2ndstage_vars "`exogvars' `endovars' `resid_powers_combo'"
init_values_2 `depvar' `2ndstage_vars' if `touse' `wgt', `nocon2' `showstages'
local yparam2 "`r(yparam)'"
local ll = `r(ll)'
local ll_0 = `r(ll_0)'
local rc = `r(rc)'
local ic = `r(ic)'
local converged = `r(converged)'
matrix `F2' = r(param_est)
matrix `F' = `F1', `F2'
local inst_opt2 "instruments(`depvar': `2ndstage_vars', `nocon2')"
* Calculate VCE
// The undocumented option -iter(0)- prevents numerical iterations -- we only want the sandwich formula.
// -haslfderivatives- differs from -hasderivatives- in that you don't need to specify each parameter, only the parameter groupings.
// These are what I call the parameters that form linear combinations.
// They are specified using the curly braces in the -parameter()- option or through residual equation names.
tempname b V N k N_clust
* -difficult_vce-
if ("`difficult_vce'" != "") {
if ("`clustervar'" == "") {
di as err "only vcetype {bf:cluster} is allowed with option {bf:difficult_vce}"
exit 198
}
forval eq = 1/`num_aux' {
forval q = 1/`=`order'-1' {
local resid_powerm1s_combo "`vhatname'`eq'_`q' `resid_powerm1s_combo'"
}
}
mata: calc_vce("`depvar'", "`2ndstage_vars'", "`allinst'", ///
"`resid_powerm1s_combo'", ///
"`resid_combo'", `order', ///
"`F'", "`touse'", "`clustervar'", ///
"`nocon1'", "`nocon2'", "`V'", ///
"`omit1'", "`omit2'")
qui gmm gmm_ivcloglog if `touse', equations(`auxdepvars' `depvar') ///
parameters("`yparam1_combined' `yparam2'") ///
y1(`depvar') y2(`auxdepvars') /// // Wooldridge-style notation for this line; the 1 and 2 on other lines represent stages
vhatname(`vhatname') order(`order') ///
`inst_opt1' `inst_opt2' ///
winitial(identity) ///
hasderivatives fakederivatives ///
onestep iter(0) from(`F') ///
vce(`vce') ///
valueid("EE criterion") `diopts'
matrix `b' = e(b)
matrix colnames `V' = `yparam1_combined' `yparam2'
matrix rownames `V' = `yparam1_combined' `yparam2'
ereturn repost b = `b' V = `V'
}
* Directly use GMM if no -difficult_vce- option
else {
qui gmm gmm_ivcloglog if `touse', equations(`auxdepvars' `depvar') ///
parameters("`yparam1_combined' `yparam2'") ///
y1(`depvar') y2(`auxdepvars') /// // Wooldridge-style notation for this line; the 1 and 2 on other lines represent stages
vhatname(`vhatname') order(`order') ///
`inst_opt1' `inst_opt2' ///
winitial(identity) ///
haslfderivatives onestep iter(0) from(`F') ///
vce(`vce') ///
valueid("EE criterion") `diopts'
}
* -ereturn- stuff and display results
// Note: in versions < 14 of Stata, -gmm- can count the number of parameters incorrectly,
// failing to exclude omission-flagged variables.
matrix `b' = e(b)
matrix `V' = e(V)
scalar `N' = e(N)
scalar `k' = e(n_moments) // Equals number of parameters but is not incorrectly counted in versions 14 of Stata
scalar `N_clust' = e(N_clust)
local clustvar `e(clustvar)'
local e_vce "`e(vce)'" // Named this way as a reminder that this is different from `vce'
local vcetype "`e(vcetype)'"
ereturn post `b' `V' /*`wgt'*/, depname(`depvar') o(`=`N'') esample(`touse') buildfvinfo
// Macros (apart from those from -ereturn post-)
ereturn hidden local vcetype "`vcetype'"
ereturn local clustvar `clustvar'
ereturn local vce "`e_vce'"
if ("`nocon1'" == "") {
ereturn local inst `allinst' _con
}
else {
ereturn local inst `allinst'
}
ereturn local auxdep `auxdepvars'
ereturn local endog `endovars'
if ("`nocon2'" == "") {
ereturn local exog `exogvars' _con
}
else {
ereturn local exog `exogvars'
}
ereturn local depvar `depvar'
ereturn local cmd "ivcloglog"
// Scalars
ereturn scalar k = `k'
ereturn scalar endog_ct = `num_endo'
ereturn scalar auxdep_ct = `num_aux'
ereturn scalar N_clust = `N_clust'
foreach var in `exogvars' `endovars' {
local temp "`temp' [`depvar']`var'"
}
qui test `temp' // Wald test for second stage, excluding constant (if it exists) and control function
ereturn scalar chi2 = r(chi2)
ereturn scalar df_m = r(df)
ereturn scalar p = chiprob(r(df), r(chi2))
ereturn scalar ll = `ll'
ereturn scalar ll_0 = `ll_0'
_post_vce_rank
ereturn scalar ic = `ic' // From -cloglog-
ereturn scalar rc = `rc' // From -cloglog-
ereturn scalar converged = `converged' // From -cloglog-
ivcloglog_out, `diopts'
end
program parse_remaining_opt, sclass // Parses -auxiliary()- and also checks for remaining options (should be none)
syntax [, exogvars(varlist fv ts) AUXiliary(string) ENDOgenous(varlist fv ts) *]
* Check for remaining options (should be none)
if ("`options'"!="") {
di as err "option(s) {bf:`options'} not allowed"
exit 198
}
* Parse contents of -auxiliary()-
gettoken auxiliary temp : auxiliary, parse(",") // Delimiter is stored in `temp' (which is not str.) and is then deleted
local options = substr("`temp'", 2, .)
parse_subopt_noconstant, `options'
local nocons1 `r(constant)' // Don't directly store `options' in `nocons1' since the former can be abbreviated
gettoken auxdepvars temp : auxiliary, parse("=")
local auxdepvars : list uniq auxdepvars // Check for duplicates
local extinst = substr("`temp'", 2, .)
local extinst : list uniq extinst // Check for duplicates
local extinst : list extinst - exogvars // Check for exogvars being explicitly added as instruments
* Count and parse auxiliary dependent variables
local num_aux : list sizeof auxdepvars
forval i = 1/`num_aux' {
local temp : word `i' of `auxdepvars' // Format is auxiliary(auxdepvar_1 ... auxdepvar_k = ...)
_fv_check_depvar `temp', k(1) // Stops the dependent variables in the auxiliary equations being factor variables
}
local num_extinst : list sizeof extinst
if (`num_aux' > `num_extinst') {
di as err "equation not identified; must have at least as many instruments not in the regression as there are instrumented variables"
exit 481
}
if ("`endogenous'" == "") {
local endogenous "`auxdepvars'"
local num_endo = `num_aux' // If option -endogenous()- not supplied, assume that the auxdepvars are the endovars
}
else {
local endogenous : list uniq endogenous // Check for duplicates
local num_endo : word count `endogenous' // Should not have quotes otherwise will return 1
}
sreturn local allinst "`exogvars' `extinst'"
sreturn local endovars "`endogenous'"
sreturn local auxdepvars "`auxdepvars'"
sreturn local nocons1 "`nocons1'"
sreturn local num_endo `num_endo'
sreturn local num_aux `num_aux'
sreturn local options `options'
end
program define parse_subopt_noconstant, rclass
syntax [anything], [noCONstant * ]
if ("`options'"!="") {
di as err "suboption(s) {bf:`options'} not allowed in option {bf:auxiliary()}"
exit 198
}
return local constant "`constant'"
end
program define iden_check_and_fvexpand_1, rclass
syntax, allinst(varlist fv ts) touse(name) [wgt(name) DIFFicult_vce noCONstant]
* Auxiliary model
cap noi _rmcoll `allinst' if `touse', expand `constant' noskipline // Note: all instruments are shared between the auxiliary models
if _rc {
error _rc
}
local vlist "`r(varlist)'"
* Recording the position of variables flagged for omission
// Constants are never flagged for omission. Since we focus on the variables that are flagged here, we don't need to worry about constants.
tempname temp allinst_omit
// "`allinst'" can obviously never be an empty string, so we do not need to worry about the problem explained below
matrix `temp' = J(1, wordcount("`vlist'"), 1)
matrix colnames `temp' = `vlist'
_ms_omit_info `temp'
matrix `allinst_omit' = `temp' - r(omit)
return local allinst "`vlist'"
return matrix allinst_omit = `allinst_omit'
end
program define iden_check_and_fvexpand_2, rclass
syntax, [lhs(varname ts) exog(varlist fv ts) endo(varlist ts) touse(name) wgt(name) DIFFicult_vce noCONstant]
* Primary model
cap noi _rmcoll `lhs' `endo' `exog' if `touse', touse(`touse') logit expand `constant' noskipline
if _rc { // Using c(rc), which means `=_rc', also works
error _rc
}
// Note: putting the exogenous variables at the end instead prioritizes them for omission.
// -logit- is an undocumented option of -_rmcoll-. Like ivprobit.ado, I use it here too.
// I use -_rmcoll, logit- instead of -_rmdcoll- and -_rmcoll- separately because -_rmdcoll- does not allow the option -touse()-.
* Extract the exogenous controls post-flagging
local vlist "`r(varlist)'" // `r(varlist)' is technically not a macro but is usually (not always) treated as such by Stata. Verbose here for safety just in case.
gettoken lhs vlist : vlist
extract_varlist, small(`endo') big(`vlist')
local endo `r(small)'
local vlist `r(big)'
local exog : copy local vlist
* Disallow dropping endogenous variables
checkifdropped, vars(`endo')
* Recording the position of variables flagged for omission
// Constants are never flagged for omission. Since we focus on the variables that are flagged here, we don't need to worry about constants.
tempname temp exog_omit
// Annoyingly, void matrices are allowed in Mata but not Stata, so I need to work around this.
if ("`exog'" != "") {
matrix `temp' = J(1, wordcount("`exog'"), 1) // Annoyingly, -_ms_omit_info- only works with matrices
matrix colnames `temp' = `exog'
_ms_omit_info `temp'
matrix `exog_omit' = `temp' - r(omit) // r(omit) from -_ms_omit_info- uses 1 to represent omission; we want it the other way round for Mata
}
return local exog "`exog'"
if ("`exog'" != "") {
return matrix exog_omit = `exog_omit'
}
end
program define extract_varlist, rclass
syntax, small(string) big(string)
local n : list sizeof small
local small // Clear this local macro since it is not empty
forval i = 1/`n' {
gettoken v big : big
local small `small' `v'
}
return local small `small'
return local big `big'
end
program define checkifdropped // Adapted from ivprobit.ado
syntax, vars(string)
while "`vars'" != "" {
gettoken var vars : vars
_ms_parse_parts `var' // Parses the matrix stripe `var' and returns the token parts in r()
if r(omit) { // From _ms_parse_parts
di as err "tried to drop an endogenous variable, which is not allowed"
exit 498
}
}
end
program define init_values_1, rclass // First stage
syntax [anything][if][in] ///
[fweight iweight pweight aweight], ///
[ ///
noCONstant ///
predictvars(string) ///
order(integer 1) ///
SHOWstages ///
]
marksample touse
if ("`weight'" != "") {
local wgt [`weight'`exp']
}
tempname param_est
if ("`showstages'" == "") {
local qui quietly
}
`qui' regress `anything' if `touse' `wgt', `constant'
matrix `param_est' = e(b) // e(b) from -regress- DOES NOT have the format "depvar: ...", so we must change the format
reformat_params // (Renaming the columns is unnecessary because -from()- extracts positionally and ignores names.)
local yparam "`r(yparam)'"
tempvar prediction
qui predict double `prediction' if `touse', resid
forval p = 1/`order' {
local predictvar_`p' : word `p' of `predictvars'
qui replace `predictvar_`p'' = `prediction'^`p'
}
return matrix param_est = `param_est'
return local yparam "`yparam'"
end
program define init_values_2, rclass // Second stage
syntax [anything][if][in] ///
[fweight iweight pweight aweight], ///
[ ///
noCONstant ///
SHOWstages ///
]
marksample touse
if ("`weight'" != "") {
local wgt [`weight'`exp']
}
tempname param_est
if ("`showstages'" == "") {
local qui quietly
}
`qui' cloglog `anything' if `touse' `wgt', `constant'
matrix `param_est' = e(b) // e(b) from -cloglog- DOES have the format "depvar: ..."
local yparam: colfullnames e(b)
return matrix param_est = `param_est'
return local yparam "`yparam'"
return local ll = e(ll)
return local ll_0 = e(ll_0)
return local ic = e(ic)
return local rc = e(rc)
return local converged = e(converged)
end
program define reformat_params, rclass
syntax [anything], [subtract(string)]
local paramy: colfullnames e(b)
local yvar "`e(depvar)'"
foreach var in `paramy' {
local yparam "`yparam' `yvar':`var'"
}
return local yparam "`yparam'"
end
mata:
//% Subroutine to get VCE
void calc_vce(string scalar y, string scalar x, string scalar z, ///
string scalar resid_powerm1s_combo, ///
string scalar resid_combo, real scalar order, ///
string scalar point_est, string scalar touse, string scalar clustervar, ///
string scalar nocon1, string scalar nocon2, string scalar vcename, ///
string scalar omit1, string scalar omit2)
{ //% Type declarations (what types these objects *initially* are)
real matrix temp_mat, X, Z, R1, info, g1_c, g2_c, g_c, Omega, vhat_powerm1s_combo_rump, vhat_powerm1s_eq, H_OLS, G_11, G_21_eq, G_21, G_22, V
real rowvector X_omit, Z_omit, multiZ_omit, ZX_omit, temp_b, b, rho_rump, rho_eq
real colvector y1, Xb, expXb, r1, r2, cvar, d, d_tilde_eq
real scalar n, num_endo, k_multiZ, k_multiZ_orig, nc, k_X, k_ZX, eq, orderm1
//% Set up matrices (representing the data)
// Note: st_data ignores factor variables with any omitted levels.
// It treats all levels as omitted, even ones with the omission operator "o", setting all their columns to zero!
st_view(y1=., ., y, touse) // Used "y" rather than "y1" as argument to avoid conflict
st_view(temp_mat=., ., x, touse) // Need "temp" because identical arguments not allowed
X_omit = st_matrix(omit2)
st_select(X, temp_mat, X_omit)
st_view(temp_mat=., ., z, touse)
Z_omit = st_matrix(omit1)
st_select(Z, temp_mat, Z_omit)
st_view(R1=., ., resid_combo, touse) // Matrix of auxiliary model residuals
n = rows(X) // All variables should have the same number of observations, conditioning on `touse'
// Add vector of ones to X and Z if constant requested.
// cross()/quadcross() can automatically do this for you, but this way is more readable.
if (nocon1 == "") {
Z = Z, J(n, 1, 1)
Z_omit = (Z_omit, 1)
}
if (nocon2 == "") {
X = X, J(n, 1, 1)
X_omit = (X_omit, 1)
}
//% Set up vectors (representing parameter estimates)
// Need to transpose because the estimates are inputted as `F', a row vector
num_endo = cols(R1)
k_multiZ = cols(Z) * num_endo // Defined in a SUR style; each instrument is counted once per eq. This makes intuitive sense if you consider having different instruments for each eq.
k_multiZ_orig = cols(Z_omit) * num_endo
temp_b = st_matrix(point_est)[|k_multiZ_orig+1 \ .|]
st_select(b, temp_b, X_omit)
b = b'
// Reminder: range subscripts are faster than list subscripts for subsetting a contiguous matrix (larger than 1x1)
// Reminder: "..." is not allowed for Mata matrix/vector subsetting
//% Omega (meat in the sandwich formula)
// We need to average the moment functions within each cluster to get the monent functions
st_view(cvar=., ., clustervar, touse)
info = panelsetup(cvar, 1) // Note: the other data must be in the same order as clustervar (which is true by default)
nc = rows(info)
k_X = cols(X)
k_ZX = k_multiZ + k_X
// Primary model, Part 1
Xb = X*b
expXb = exp(Xb)
r2 = mm_cond(y1, exp(Xb - expXb) :/ (-expm1(-expXb)), -expXb) // Not true residuals but pseudo-residuals
// Make by-cluster superobservations and then we can get the "robust SE"-style Gram matrix
g1_c = J(nc, 0, 0)
for (eq=1; eq<=num_endo; eq++) {
r1 = R1[., eq] // Residuals corresponding to auxiliary equation #eq
g1_c = (g1_c, panelsum(Z, r1, info)) // "c" here means summed within cluster (with R1 as the weights)
}
g2_c = panelsum(X, r2, info)
g_c = (g1_c, g2_c)
Omega = quadcross(g_c, g_c)
//% G_inv (bread in the sandwich formula); G = (G_11, G_12 \ G_21, G_22)
//% G_11
H_OLS = -quadcross(Z, Z) // Minus sign in this Hessian comes from maximizing the negative of the OLS objective function
G_11 = I(num_endo) # H_OLS // Used Kronecker product operator; could probably be optimized by using a loop with blockdiag() instead
//% G_12 is just a zero matrix
//% G_21
d = mm_cond(y1, ///
r2 :* (expm1(Xb) + exp(-expXb)) :/ expm1(-expXb), ///
-expXb)
// d is a vector of values "C_i", defined in the gmm moment evaluator function.
// I changed the notation here to avoid confusion since "c" was already used above.
if (order > 1) {
st_view(vhat_powerm1s_combo_rump, ., resid_powerm1s_combo, touse) // This contain the rest when we partition off a subvector from the front (like -gettoken-)
}
orderm1 = order - 1
// rho_rump is defined similarly to the above. (Also, order*num_endo is the size of rho)
if (nocon2 == "") {
rho_rump = b[|k_X-order*num_endo \ k_X-1|]
}
else {
rho_rump = b[|k_X-order*num_endo+1 \ .|]
}
G_21 = J(k_X, 0, .) // Initialize as zero-column matrix so it can be used in the loop
for (eq=1; eq<=num_endo; eq++) {
if (order > 1) {
vhat_powerm1s_eq = vhat_powerm1s_combo_rump[|1, 1 \ ., orderm1|]
if (eq < num_endo) {
vhat_powerm1s_combo_rump = vhat_powerm1s_combo_rump[|1, order \ ., .|] // Only take the rest if that actually exists
}
vhat_powerm1s_eq = J(n, 1, 1), vhat_powerm1s_eq // v^(p-1), p = 1, ...
}
else {
vhat_powerm1s_eq = J(n, 1, 1)
}
rho_eq = rho_rump[|1 \ order|]
if (eq < num_endo) {
rho_rump = rho_rump[|order+1 \ .|] // Only take the rest if that actually exists
}
d_tilde_eq = -d :* (vhat_powerm1s_eq * rho_eq)
G_21_eq = quadcross(X, d_tilde_eq, Z)
G_21 = (G_21, G_21_eq)
}
//% G_22
G_22 = quadcross(X, d, X)
//% V
// I use -lusolve()- rather than actually inverting G because this is more efficient.
// In addition, I use 0 tolerance because we already checked for collinearity with -_rmcoll-.
// Changing the tolerance is important because the data matrix might pass the -_rmcoll- collinearity check but
// fail the built-in collinearity check in -lusolve()-.
V = lusolve_special(G_11, G_21, G_22, lusolve_special(G_11, G_21, G_22, Omega, 0)', 0) // Would be more efficient to reuse the L and U instead of doing the same LU decomposition twice
// This V is currently actually only the VCE for the non-omission-flagged variables
multiZ_omit = J(1, 0, .)
for (eq=1; eq<=num_endo; eq++) {
multiZ_omit = (multiZ_omit, Z_omit)
}
ZX_omit = (multiZ_omit, X_omit)
V = insert_col_cons(V, ZX_omit, 0)
V = insert_col_cons(V', ZX_omit, 0)
V = makesymmetric(V)
st_matrix(vcename, V)
}
//% Subroutine similar to -lusolve()- that exploits the block lower triangular nature of G
real matrix lusolve_special(real matrix A1, real matrix A3, real matrix A4, ///
real matrix B, real scalar tolerance)
{
//% Type declarations
real matrix B1, B2, B3, B4, X1, X2, X3, X4, X
real scalar m1, n1
//% Main code
m1 = rows(A1)
n1 = cols(A3)
B1 = B[|1, 1 \ m1, n1|]
B2 = B[|1, n1+1 \ m1, .|]
B3 = B[|m1+1, 1 \ ., n1|]
B4 = B[|m1+1, n1+1 \ ., .|]
X1 = lusolve(A1, B1, tolerance)
X2 = lusolve(A1, B2, tolerance)
X3 = lusolve(A4, B3-A3*X1, tolerance)
X4 = lusolve(A4, B4-A3*X2, tolerance)
X = (X1, X2 \ X3, X4)
return(X)
}
//% Subroutine that inserts columns of a constant in a matrix corresponding to the 0s in a specified vector
real matrix insert_col_cons(real matrix M, real matrix v, real scalar cons)
{
//% Type declarations
real matrix R, temp
real scalar end_num_cols, num_rows, i
//% Main code
end_num_cols = cols(v)
num_rows = rows(M)
R = J(num_rows, end_num_cols, cons) // Easier to start with matrix of constants with final dimensions than to "expand" the original matrix
temp = M, J(num_rows, 1, 0) // Hacky buffer to prevent the loop from breaking at the last column of M
for (i=1; i<=end_num_cols; i++) {
if (v[i] == 1) {
R[., i] = temp[., 1]
temp = temp[|1, 2 \ ., .|]
}
}
return(R)
}
end
program ivcloglog_out // Modified version of -GMMout- from gmm.ado
syntax [, level(cilevel) COEFLegend selegend *]
_get_diopts diopts, `coeflegend' `selegend' `options'
local EXPRESSION = 1
local PROGRAM = 2
local robust
local hac
local bs
if bsubstr("`e(vcetype)'", 1, 6) == "Robust" | "`e(clustvar)'" != "" {
local robust "yes"
}
if bsubstr(`"`e(vcetype)'"', 1, 9) == "Bootstrap" | ///
bsubstr(`"`e(vcetype)'"', 1, 9) == "Jackknife" {
local bs "yes"
}
local anyrobust 0
display ""
display "{bf:Two-step cloglog model with endogenous covariates}"
display ""
tempname left right
.`left' = {}
.`right' = {}
local C1 "_col(1)"
local C2 "_col(21)"
local C3 "_col(43)"
local C4 "_col(64)"
.`left'.Arrpush `C1' as text "2nd-stage Wald chi2" `C2' "=" as res %14.2f e(chi2)
.`left'.Arrpush `C1' as text "2nd-stage Wald DF" `C2' "=" as res %14.0fc e(df_m)
.`left'.Arrpush `C1' as text "Prob > chi2" `C2' "=" as res %14.4f e(p)
.`right'.Arrpush `C3' as text "Number of parameters" `C4' "=" as res %14.0fc e(k)
.`right'.Arrpush `C3' as text "Number of obs" `C4' "=" as res %14.0fc e(N)
.`right'.Arrpush `C3' as text "Log-likelihood" `C4' "=" as res %14.2f e(ll)
local nl = `.`left'.arrnels'
local nr = `.`right'.arrnels'
local k = max(`nl', `nr')
forvalues i = 1/`k' {
di as text `.`left'[`i']' as text `.`right'[`i']'
}
di
_coef_table, level(`level') `diopts' // Stealthily uses e(vcetype), e(N_clust), and e(clustvar)
di in smcl "{p 0 4 4}{txt}Instruments: {res}`e(inst)'{p_end}"
ml_footnote
end
program define drop_cfvars
syntax anything, orig_varlist(varlist) [NOGENerate *]
if "`nogenerate'" != "" {
keep `orig_varlist'
// Yes, -nogenerate- is a misnomer because we do actually generate the powers of vhat. But, to be fair, that's also true for Stata tempvars.
// We do this because it's easier than renaming the columns of ("column-restriping") the matrix holding the coefficient estimates,
// which is necessary to make the powers of vhat be called "`vhatname'`i'_`p'".
}
end