Comparar texto

Encontre a diferença entre dois arquivos de texto

Real-time diff

Unified diff

Collapse lines

Highlight change

Syntax highlighting

Ferramentas

Diffchecker Desktop The most secure way to run Diffchecker. Get the Diffchecker Desktop app: your diffs never leave your computer!Get Desktop

Comparing sensitive data, confidential files or internal emails?

Most legal and privacy policies prohibit uploading sensitive data online. Diffchecker Desktop ensures your confidential information never leaves your computer. Work offline and compare documents securely.

Download Diffchecker Desktop Learn more

Untitled diff

Created 6 years agoDiff never expires

81 removals

745 lines

77 additions

738 lines

*! ivreghdfe 1.0.0 07jul2018

*! this just adds absorb() to this code:

*! ivreg2 4.1.10 9Feb2016

*! authors cfb & mes

*! see end of file for version comments

* Variable naming:

* lhs = LHS endogenous

* endo = X1, RHS endogenous (instrumented) = #K1

* inexog = X2 = Z2 = included exogenous (instruments) = #K2 = #L2

* exexog = Z1 = excluded exogenous (instruments) = #L1

* iv = {inexog exexog} = all instruments

* rhs = {endo inexog} = RHS regressors

* no 0 or 1 at end of varlist means original varlist but after expansion of FV and TS vars

* 0 at the end of the name means the varlist after duplicates removed and collinearities/omitteds marked

* 1 means the same as 0 but after omitted vars dropped and extraneous FV operators "o", "b" and "n" removed.

* 0, 1 etc. also apply to _ct variables that are counts of these varlists

* dofminus is large-sample adjustment (e.g., #fixed effects)

* sdofminus is small-sample adjustment (e.g., #partialled-out regressors)

if c(version) < 12 & c(version) >= 9 {

* livreg2 Mata library.

* Ensure Mata library is indexed if new install.

* Not needed for Stata 12+ since ssc.ado does this when installing.

capture mata: mata drop m_calckw()

capture mata: mata drop m_omega()

capture mata: mata drop ms_vcvorthog()

capture mata: mata drop s_vkernel()

capture mata: mata drop s_cdsy()

mata: mata mlib index

}

*********************************************************************************

***************************** PARENT IVREG2 *************************************

****************** FORKS TO EXTERNAL IVREG2S IF CALLER < 11 *********************

*********************************************************************************

* Parent program, forks to versions as appropriate after version call

* Requires byable(onecall)

program define ivreg2, eclass byable(onecall) /* properties(svyj) */ sortpreserve

program define ivreghdfe, eclass byable(onecall) /* properties(svyj) */ sortpreserve

local lversion 04.1.10

* Minimum of version 8 required for parent program (earliest ivreg2 is ivreg28)

version 8

* Replay = no arguments before comma

if replay() {

* Call to ivreg2 will either be for version, in which case there should be no other arguments,

* or a postestimation call, in which case control should pass to main program.

syntax [, VERsion * ]

if "`version'"~="" & "`options'"=="" {

* Call to ivreg2 is for version

di in gr "`lversion'"

ereturn clear

ereturn local version `lversion'

exit

}

else if "`version'"~="" & "`options'"~="" {

* Improper use of version option

di as err "invalid syntax - cannot combine version with other options"

exit 198

}

else {

* Postestimation call, so put `options' macro (i.e. *) back into `0' macro with preceding comma

local 0 `", `options'"'

}

* replay can't be combined with by

if replay() & _by() {

di as err "invalid syntax - cannot use by with replay"

exit 601

}

* Handling of by. ivreg2x programs are byable(recall), so must set prefix for them.

if _by() {

local BY `"by `_byvars'`_byrc0':"'

}

* If calling version is < 11, pass control to earlier version

* Note that this means calls from version 11.0 will not go to legacy version

* but will fail requirement of version 11.2 in main code.

if _caller() < 11 {

local ver = round(_caller())

local ivreg2cmd ivreg2`ver'

* If replay, change e(cmd) macro to name of legacy ivreg2 before calling it, then change back

* Note by not allowed with replay; caught above so prefix not needed here.

if replay() {

ereturn local cmd "`ivreg2cmd'"

`ivreg2cmd' `0'

ereturn local cmd "ivreg2"

}

else {

* If not replay, call legacy ivreg2 and then add macros

`BY' `ivreg2cmd' `0'

ereturn local cmd "ivreg2"

ereturn local cmd "ivreghdfe"

ereturn local ivreg2cmd "`ivreg2cmd'"

ereturn local version `lversion'

ereturn local predict ivreg2_p

}

exit

}

// Version is 11 or above.

// Pass control to current estimation program ivreg211.

if replay() {

ivreg211 `0'

}

// If not replay, call ivreg211 and then add macros

else {

`BY' ivreg211 `0'

ereturn local cmd "ivreg2"

cap mata: mata drop HDFE // prefix ivreg211 call with capture?

ereturn local ivreg2cmd "ivreg2"

cap mata: mata drop hdfe_residuals

ereturn local cmd "ivreghdfe"

ereturn local ivreg2cmd "ivreghdfe"

ereturn local version `lversion'

ereturn local predict ivreg2_p

if (e(N_hdfe)!= .) ereturn local predict reghdfe_p

}

end

*********************************************************************************

*************************** END PARENT IVREG2 ***********************************

*********************************************************************************

********************* EXIT IF STATA VERSION < 11 ********************************

* When do file is loaded, exit here if Stata version calling program is < 11.

* Prevents loading of rest of program file (could cause earlier Statas to crash).

if c(stata_version) < 11 {

exit

}

******************** END EXIT IF STATA VERSION < 11 *****************************

*********************************************************************************

***************** BEGIN MAIN IVREG2 ESTIMATION CODE *****************************

*********************************************************************************

* Main estimation program

program define ivreg211, eclass byable(recall) sortpreserve

version 11.2

local ranktestversion 01.3.02

local ivreg2cmd "ivreg211" // actual command name

local ivreg2name "ivreg2" // name used in command line and for default naming of equations etc.

local ranktestcmd "ranktest"

if replay() {

syntax [, ///

FIRST FFIRST RF SFIRST ///

dropfirst droprf dropsfirst ///

Level(integer $S_level) ///

NOHEader NOFOoter ///

EForm(string) PLUS ///

NOOMITTED vsquish noemptycells ///

baselevels allbaselevels ///

VERsion ///

]

if "`version'" != "" & "`first'`ffirst'`rf'`noheader'`nofooter'`dropfirst'`droprf'`eform'`plus'" != "" {

di as err "option version not allowed"

error 198

}

if "`version'" != "" {

di in gr "`lversion'"

ereturn clear

ereturn local version `lversion'

exit

}

if `"`e(cmd)'"' != "ivreg2" {

if `"`e(cmd)'"' != "ivreghdfe" {

error 301

}

// Set display options

local dispopt eform(`eform') `noomitted' `vsquish' `noemptycells' `baselevels' `allbaselevels'

// On replay, set flag so saved eqns aren't dropped

if "`e(firsteqs)'" != "" & "`dropfirst'" == "" {

local savefirst "savefirst"

}

if "`e(rfeq)'" != "" & "`droprf'" == "" {

local saverf "saverf"

}

if "`e(sfirsteq)'" != "" & "`dropsfirst'" == "" {

local savesfirst "savesfirst"

}

// On replay, re-display collinearities and duplicates messages

DispCollinDups

}

else {

// MAIN CODE BLOCK

local cmdline `ivreg2name' `0' // `0' rather than `*' in case of any "s in string

// Confirm ranktest is installed (necessary component).

checkversion_ranktest `ranktestcmd' `ranktestversion'

// Start parsing

syntax [anything(name=0)] [if] [in] [aw fw pw iw/] [, ///

NOID NOCOLLIN ///

FIRST FFIRST SAVEFIRST SAVEFPrefix(name) ///

RF SAVERF SAVERFPrefix(name) ///

SFIRST SAVESFIRST SAVESFPrefix(name) ///

SMall NOConstant ///

Robust CLuster(varlist) kiefer dkraay(integer 0) ///

BW(string) kernel(string) center ///

GMM GMM2s CUE ///

LIML COVIV FULLER(real 0) Kclass(real 0) ///

ORTHOG(string) ENDOGtest(string) REDundant(string) ///

PARTIAL(string) FWL(string) ///

Absorb(string) ///

Level(integer $S_level) ///

NOHEader NOFOoter NOOUTput ///

bvclean NOOMITTED omitted vsquish noemptycells ///

baselevels allbaselevels ///

title(string) subtitle(string) ///

DEPname(string) EForm(string) PLUS ///

Tvar(varname) Ivar(varname) ///

B0(string) SMATRIX(string) WMATRIX(string) ///

sw psd0 psda useqr ///

dofminus(integer 0) sdofminus(integer 0) ///

NOPARTIALSMALL ///

fvall fvsep ///

]

if (`"`absorb'"' != "") {

// ensure reghdfe is installed

capture which reghdfe

if c(rc) {

di as err "Error: must have reghdfe installed"

di as err "To install, from within Stata type " _c

//di in smcl "{stata ssc install reghdfe :ssc install reghdfe}"

loc url "https://github.com/sergiocorreia/reghdfe/raw/version-4/src/"

di in smcl `"{stata `"net install reghdfe, from(`url')"':net install reghdfe, from(`url')}"'

exit 601

}

cap reghdfe, requirements // ensure ftools, moresyntax, etc. are installed

if (c(rc)) reghdfe, requirements

reghdfe, check // ensure the .mlib exists

// absorb implies...

local small small

local noconstant noconstant

local nopartialsmall

}

// Parse after clearing any sreturn macros (can be left behind in Stata 11)

sreturn clear

ivparse `0', ivreg2name(`ivreg2name') /// needed for some options

partial(`partial') ///

fwl(`fwl') /// legacy option

orthog(`orthog') ///

endogtest(`endogtest') ///

redundant(`redundant') ///

depname(`depname') ///

`robust' ///

cluster(`cluster') ///

bw(`bw') ///

kernel(`kernel') ///

dkraay(`dkraay') ///

`center' ///

`kiefer' ///

`sw' ///

`noconstant' ///

tvar(`tvar') ///

ivar(`ivar') ///

`gmm2s' ///

`gmm' /// legacy option, produces error message

`cue' ///

`liml' ///

fuller(`fuller') ///

kclass(`kclass') ///

b0(`b0') ///

wmatrix(`wmatrix') ///

`noid' ///

`savefirst' ///

savefprefix(`savefprefix') ///

`saverf' ///

saverfprefix(`saverfprefix') ///

`savesfirst' ///

savesfprefix(`savesfprefix') ///

dofminus(`dofminus') ///

`psd0' ///

`psda' ///

`nocollin' ///

`useqr' ///

`bvclean' ///

eform(`eform') ///

`noomitted' ///

`vsquish' ///

`noemptycells' ///

`baselevels' ///

`allbaselevels'

// varlists are unexpanded; may be empty

local lhs `s(lhs)'

local depname `s(depname)'

local endo `s(endo)'

local inexog `s(inexog)'

local exexog `s(exexog)'

local partial `s(partial)'

local cons =s(cons)

local partialcons =s(partialcons)

local tvar `s(tvar)'

local ivar `s(ivar)'

local tdelta `s(tdelta)'

local tsops =s(tsops)

local fvops =s(fvops)

local robust `s(robust)'

local cluster `s(cluster)'

local bw =`s(bw)' // arrives as string but return now as number

local bwopt `s(bwopt)'

local kernel `s(kernel)' // also used as flag for HAC estimation

local center =`s(center)' // arrives as string but now boolean

local kclassopt `s(kclassopt)'

local fulleropt `s(fulleropt)'

local liml `s(liml)'

local noid `s(noid)' // can also be triggered by b0(.) option

local useqr =`s(useqr)' // arrives as string but now boolean; nocollin=>useqr

local savefirst `s(savefirst)'

local savefprefix `s(savefprefix)'

local saverf `s(saverf)'

local saverfprefix `s(saverfprefix)'

local savesfirst `s(savesfirst)'

local savesfprefix `s(savesfprefix)'

local psd `s(psd)' // triggered by psd0 or psda

local dofmopt `s(dofmopt)'

local bvclean =`s(bvclean)' // arrives as string but return now as boolean

local dispopt `s(dispopt)'

// Can now tsset; sortpreserve will restore sort after exit

if `tsops' | "`kernel'"~="" {

cap tsset // restores sort if tsset or xtset but sort disrupted

if _rc>0 {

tsset `ivar' `tvar'

}

***********************************************************

// Weights

// fweight and aweight accepted as is

// iweight not allowed with robust or gmm and requires a trap below when used with summarize

// pweight is equivalent to aweight + robust

// Since we subsequently work with wvar, tsrevar of weight vars in weight `exp' not needed.

tempvar wvar

if "`weight'" == "fweight" | "`weight'"=="aweight" {

local wtexp `"[`weight'=`exp']"'

qui gen double `wvar'=`exp'

}

if "`weight'" == "fweight" & "`kernel'" !="" {

di in red "fweights not allowed (data are -tsset-)"

exit 101

}

if "`weight'" == "fweight" & "`sw'" != "" {

di in red "fweights currently not supported with -sw- option"

exit 101

}

if "`weight'" == "iweight" {

if "`robust'`cluster'`gmm2s'`kernel'" !="" {

di in red "iweights not allowed with robust or gmm"

exit 101

}

else {

local wtexp `"[`weight'=`exp']"'

qui gen double `wvar'=`exp'

}

if "`weight'" == "pweight" {

local wtexp `"[aweight=`exp']"'

qui gen double `wvar'=`exp'

local robust "robust"

}

if "`weight'" == "" {

* If no weights, define neutral weight variable

qui gen byte `wvar'=1

}

********************************************************************************

// markout sample

// include `tvar' to limit sample to where tvar is available, but only if TS operators used

marksample touse

if `tsops' {

markout `touse' `lhs' `inexog' `exexog' `endo' `cluster' `tvar', strok

}

else {

markout `touse' `lhs' `inexog' `exexog' `endo' `cluster', strok

}

* Create HDFE object and update touse

if (`"`absorb'"' != "") {

* fixed_effects(absvars | , touse, wtype, wtvar, dropsing, verbose)

mata: HDFE = fixed_effects("`absorb'", "`touse'", "`weight'", "`wvar'")

_assert ("`s(options)'"==""), msg("unsupported options in {bf:absorb()}: `s(options)'")

mata: HDFE.clustervars = HDFE.base_clustervars = tokens("`cluster'")

mata: HDFE.num_clusters = length(HDFE.clustervars)

mata: HDFE.estimate_dof() // compute degrees-of-freedom

mata: HDFE.save_touse("`touse'", 1) // 1 = overwrite

}

********************************************************************************

// weight factor and sample size

// Every time a weight is used, must multiply by scalar wf ("weight factor")

// wf=1 for no weights, fw and iw, wf = scalar that normalizes sum to be N if aw or pw

sum `wvar' if `touse' `wtexp', meanonly

// Weight statement

if "`weight'" ~= "" {

di in gr "(sum of wgt is " %14.4e `r(sum_w)' ")"

}

if "`weight'"=="" | "`weight'"=="fweight" | "`weight'"=="iweight" {

// Effective number of observations is sum of weight variable.

// If weight is "", weight var must be column of ones and N is number of rows

local wf=1

local N=r(sum_w)

}

else if "`weight'"=="aweight" | "`weight'"=="pweight" {

local wf=r(N)/r(sum_w)

local N=r(N)

}

else {

// Should never reach here

di as err "ivreg2 error - misspecified weights"

exit 198

}

if `N'==0 {

di as err "no observations"

exit 2000

}

***************************************************************

// Time-series data

// tindex used by Mata code so that ts operators work correctly

tempvar tindex

qui gen `tindex'=1 if `touse'

qui replace `tindex'=sum(`tindex') if `touse'

if `tsops' | "`kernel'"~="" {

// Report gaps in data

tsreport if `touse', panel

if `r(N_gaps)' != 0 {

di as text "Warning: time variable " as res "`tvar'" as text " has " ///

as res "`r(N_gaps)'" as text " gap(s) in relevant range"

}

// Set local macro T and check that bw < (T-1)

sum `tvar' if `touse', meanonly

local T = r(max)-r(min) + 1

local T1 = `T' - 1

if (`bw' > (`T1'/`tdelta')) {

di as err "invalid bandwidth in option bw() - cannot exceed timespan of data"

exit 198

}

// kiefer VCV = kernel(tru) bw(T) and no robust with tsset data

if "`kiefer'" ~= "" {

local bw =`T'

}

*********** Column of ones for constant set up here **************

if "`noconstant'"=="" {

// If macro not created, automatically omitted.

tempvar ones

qui gen byte `ones' = 1 if `touse'

}

************* Varlists, FV varlists, duplicates *****************

// Varlists come in 4 versions, e.g., for inexog:

// (a) inexog = full list of original expanded vnames; may have duplicates

// (b) inexog0 = as with inexog with duplicates removed but RETAINING base/omitted/etc. varnames

// (c) inexog1 = as with inexog0 but WITHOUT base/omitted/etc.

// (d) fv_inexog1 = corresponding list with temp vars minus base/omitted/etc., duplicates, collinearities etc.

// Varlists (c) and (d) are definitive, i.e., have the variables actually used in the estimation.

// Create consistent expanded varlists.

// "Consistent" means base vars for FVs must be consistent

// hence default rhs=endo+inexog is expanded as one.

// fvall: overrides, endo+inexog+exexog expanded as one

// fvsep: overrides, endo, inexog and exexog expanded separately

// NB: expanding endo+inexog+exexog is dangerous because

// fvexpand can zap a list in case of overlap

// e.g. fvexpand mpg + i(1/4).rep78 + i5.rep78

// => mpg 1b.rep78 2.rep78 3.rep78 4.rep78 5.rep78

// but fvexpand mpg + i.rep78 + i5.rep78

// => mpg 5.rep78

CheckDupsCollin, ///

lhs(`lhs') ///

endo(`endo') ///

inexog(`inexog') ///

exexog(`exexog') ///

partial(`partial') ///

orthog(`orthog') ///

endogtest(`endogtest') ///

redundant(`redundant') ///

touse(`touse') ///

wvar(`wvar') ///

wf(`wf') ///

`noconstant' ///

`nocollin' ///

`fvall' ///

`fvsep'

// Replace basic varlists and create "0" versions of varlists

foreach vl in lhs endo inexog exexog partial orthog endogtest redundant {

local `vl' `s(`vl')'

local `vl'0 `s(`vl'0)'

}

local dups `s(dups)'

local collin `s(collin)'

local ecollin `s(ecollin)'

// Create "1" and fv versions of varlists

foreach vl in lhs endo inexog exexog partial orthog endogtest redundant {

foreach var of local `vl'0 { // var-by-var so that fvrevar doesn't decide on base etc.

_ms_parse_parts `var'

if ~`r(omit)' { // create temp var only if not omitted

fvrevar `var' if `touse'

local `vl'1 ``vl'1' `var'

local fv_`vl'1 `fv_`vl'1' `r(varlist)'

}

local `vl'1 : list retokenize `vl'1

local fv_`vl'1 : list retokenize fv_`vl'1

}

// Check that LHS expanded to a single variable

local wrongvars_ct : word count `lhs'

if `wrongvars_ct' > 1 {

di as err "multiple dependent variables specified: `lhs'"

error 198

}

// Check that option varlists are compatible with main varlists

// orthog()

local wrongvars : list orthog1 - inexog1

local wrongvars : list wrongvars - exexog1

local wrongvars_ct : word count `wrongvars'

if `wrongvars_ct' {

di as err "Error: `wrongvars' listed in orthog() but does not appear as exogenous."

error 198

}

// endog()

local wrongvars : list endogtest1 - endo1

local wrongvars_ct : word count `wrongvars'

if `wrongvars_ct' {

di as err "Error: `wrongvars' listed in endog() but does not appear as endogenous."

error 198

}

// redundant()

local wrongvars : list redundant1 - exexog1

local wrongvars_ct : word count `wrongvars'

if `wrongvars_ct' {

di as err "Error: `wrongvars' listed in redundant() but does not appear as exogenous."

error 198

}

// And create allnames macros

local allnames `lhs' `endo' `inexog' `exexog'

local allnames0 `lhs0' `endo0' `inexog0' `exexog0'

local allnames1 `lhs1' `endo1' `inexog1' `exexog1'

local fv_allnames1 `fv_lhs1' `fv_endo1' `fv_inexog1' `fv_exexog1'

// *************** Partial-out block ************** //

// `partial' has all to be partialled out except for constant

if "`partial1'" != "" | `partialcons'==1 {

if "`partial1'" != "" | `partialcons'==1 | "`absorb'" != "" {

preserve

// Remove partial0 from inexog0.

// Remove partial1 from inexog1.

local inexog0 : list inexog0 - partial0

local inexog1 : list inexog1 - partial1

local fv_inexog1 : list fv_inexog1 - fv_partial1

// Check that cluster, weight, tvar or ivar variables won't be transformed

// Use allnames1 (expanded varlist)

if "`cluster'"~="" {

local pvarcheck : list cluster in allnames1

if `pvarcheck' {

di in r "Error: cannot use cluster variable `cluster' as dependent variable, regressor or IV"

di in r " in combination with -partial- option."

error 198

}

if "`tvar'"~="" {

local pvarcheck : list tvar in allnames1

if `pvarcheck' {

di in r "Error: cannot use time variable `tvar' as dependent variable, regressor or IV"

di in r " in combination with -partial- option."

error 198

}

if "`ivar'"~="" {

local pvarcheck : list ivar in allnames1

if `pvarcheck' {

di in r "Error: cannot use panel variable `ivar' as dependent variable, regressor or IV"

di in r " in combination with -partial- option."

error 198

}

if "`wtexp'"~="" {

tokenize `exp', parse("*/()+-^&|~")

local wvartokens `*'

local nwvarnames : list allnames1 - wvartokens

local wvarnames : list allnames1 - nwvarnames

if "`wvarnames'"~="" {

di in r "Error: cannot use weight variables as dependent variable, regressor or IV"

di in r " in combination with -partial- option."

error 198

}

// Partial out

// But first replace everything with doubles

recast double `fv_lhs1' `fv_endo1' `fv_inexog1' `fv_exexog1' `fv_partial1'

mata: s_partial ("`fv_lhs1'", ///

"`fv_endo1'", ///

"`fv_inexog1'", ///

"`fv_exexog1'", ///

"`fv_partial1'", ///

"`touse'", ///

"`weight'", ///

"`wvar'", ///

`wf', ///

`N', ///

`cons')

local partial_ct : word count `partial1'

if ("`absorb'" != "") {

local hdfe_varlist `fv_lhs1' `fv_endo1' `fv_inexog1' `fv_exexog1' `fv_partial1'

mata: st_store(HDFE.sample, tokens("`hdfe_varlist'"), HDFE.partial_out(tokens("`hdfe_varlist'")))

mata: st_local("absorb_ct", strofreal(HDFE.df_a))

assert `absorb_ct'`' != .

if (`absorb_ct'==0) local absorb_ct 1 // adjustment to match ivreg2 and old reghdfe (happens if absvar is nested in cluster)

local partial_ct 0

local partialcons `absorb_ct'

}

if ("`partial1'" != "" | `partialcons'==1) {

mata: s_partial ("`fv_lhs1'", ///

"`fv_endo1'", ///

"`fv_inexog1'", ///

"`fv_exexog1'", ///

"`fv_partial1'", ///

"`touse'", ///

"`weight'", ///

"`wvar'", ///

`wf', ///

`N', ///

`cons')

local partial_ct : word count `partial1'

}

// Constant is partialled out, unless nocons already specified in the first place

capture drop `ones'

local ones ""

if "`noconstant'" == "" {

// partial_ct used for small-sample adjustment to regression F-stat

local partial_ct = `partial_ct' + 1

local noconstant "noconstant"

local cons 0

}

if ("`absorb'" != "") {

local partial_ct = `partial_ct' + `absorb_ct'

}

else {

// Set count of partial vars to zero if option not used

local partial_ct 0

local partialcons 0

}

// Add partial_ct to small dof adjustment sdofminus

if "`nopartialsmall'"=="" {

local sdofminus = `sdofminus'+`partial_ct'

}

*********************************************

local rhs0 `endo0' `inexog0' // needed for display of omitted/base/etc.

local rhs1 `endo1' `inexog1'

local insts1 `exexog1' `inexog1'

local fv_insts1 `fv_exexog1' `fv_inexog1'

local fv_rhs1 `fv_endo1' `fv_inexog1'

local rhs0_ct : word count `rhs0' // needed for display of omitted/base/etc.

local rhs1_ct : word count `fv_rhs1'

local iv1_ct : word count `fv_insts1'

local endo1_ct : word count `fv_endo1'

local exex1_ct : word count `fv_exexog1'

local endoexex1_c : word count `fv_endo1' `fv_exexog1'

local inexog1_ct : word count `fv_inexog1'

// Counts modified to include constant if appropriate

local rhs1_ct = `rhs1_ct' + `cons'

local rhs0_ct = `rhs0_ct' + `cons' // needed for display of omitted/base/etc.

local iv1_ct = `iv1_ct' + `cons'

// Column/row names for matrices b, V, S, etc.

local cnb0 `endo0' `inexog0' // including omitted

local cnb1 `endo1' `inexog1' // excluding omitted

local cnZ0 `exexog0' `inexog0' // excluding omitted

local cnZ1 `exexog1' `inexog1' // excluding omitted

if `cons' {

local cnb0 "`cnb0' _cons"

local cnb1 "`cnb1' _cons"

local cnZ0 "`cnZ0' _cons"

local cnZ1 "`cnZ1' _cons"

}

*********************************************

// Remaining checks: variable counts, col/row names of b0, smatrix, wmatrix

CheckMisc, ///

rhs1_ct(`rhs1_ct') ///

iv1_ct(`iv1_ct') ///

bvector(`b0') ///

smatrix(`smatrix') ///

wmatrix(`wmatrix') ///

cnb1(`cnb1') ///

cnZ1(`cnZ1')

if "`b0'"~="" {

tempname b0 // so we can overwrite without changing original user matrix

mat `b0' = r(b0)

}

if "`smatrix'"~="" {

tempname S0

mat `S0' = r(S0)

}

if "`wmatrix'"~="" {

tempname wmatrix // so we can overwrite without changing original user matrix

mat `wmatrix' = r(W0)

}

*************** Commonly used matrices ****************

tempname YY yy yyc

tempname XX X1X1 X2X2 X1Z X1Z1 XZ Xy

tempname ZZ Z1Z1 Z2Z2 Z1Z2 Z1X2 Zy ZY Z2y Z2Y

tempname XXinv X2X2inv ZZinv XPZXinv

tempname rankxx rankzz condxx condzz

// use fv_ varlists

mata: s_crossprods ("`fv_lhs1'", ///

"`fv_endo1'", ///

"`fv_inexog1' `ones'", ///

"`fv_exexog1'", ///

"`touse'", ///

"`weight'", ///

"`wvar'", ///

`wf', ///

`N')

mat `XX' =r(XX)

mat `X1X1' =r(X1X1)

mat `X1Z' =r(X1Z)

mat `ZZ' =r(ZZ)

mat `Z2Z2' =r(Z2Z2)

mat `Z1Z2' =r(Z1Z2)

mat `XZ' =r(XZ)

mat `Xy' =r(Xy)

mat `Zy' =r(Zy)

mat `YY' =r(YY)

scalar `yy' =r(yy)

scalar `yyc' =r(yyc)

mat `ZY' =r(ZY)

mat `Z2y' =r(Z2y)

mat `Z2Y' =r(Z2Y)

mat `XXinv' =r(XXinv)

mat `ZZinv' =r(ZZinv)

mat `XPZXinv' =r(XPZXinv)

scalar `condxx' =r(condxx)

scalar `condzz' =r(condzz)

scalar `rankzz' = rowsof(`ZZinv') - diag0cnt(`ZZinv')

scalar `rankxx' = rowsof(`XXinv') - diag0cnt(`XXinv')

local overid = `rankzz' - `rankxx'

********** CLUSTER SETUP **********************************************

* Mata code requires data are sorted on (1) the first var cluster if there

* is only one cluster var; (2) on the 3rd and then 1st if two-way clustering,

* unless (3) two-way clustering is combined with kernel option, in which case

* the data are tsset and sorted on panel id (first cluster variable) and time

* id (second cluster variable).

* Second cluster var is optional and requires an identifier numbered 1..N_clust2,

* unless combined with kernel option, in which case it's the time variable.

* Third cluster var is the intersection of 1 and 2, unless combined with kernel

* opt, in which case it's unnecessary.

* Sorting on "cluster3 cluster1" means that in Mata, panelsetup works for

* both, since cluster1 nests cluster3.

* Note that it is possible to cluster on time but not panel, in which case

* cluster1 is time, cluster2 is empty and data are sorted on panel-time.

* Note also that if data are sorted here but happen to be tsset, will need

* to be re-tsset after estimation code concludes.

// No cluster options or only 1-way clustering

// but for Mata and other purposes, set N_clust vars =0

local N_clust=0

local N_clust1=0

local N_clust2=0

if "`cluster'"!="" {

local clopt "cluster(`cluster')"

tokenize `cluster'

local cluster1 "`1'"

local cluster2 "`2'"

if "`kernel'"~="" {

* kernel requires

Diffs salvos

Texto original

Open file

*! ivreg2 4.1.10  9Feb2016
*! authors cfb & mes
*! see end of file for version comments

*  Variable naming:
*  lhs = LHS endogenous
*  endo = X1, RHS endogenous (instrumented) = #K1
*  inexog = X2 = Z2 = included exogenous (instruments) = #K2 = #L2
*  exexog = Z1 = excluded exogenous (instruments) = #L1
*  iv = {inexog exexog} = all instruments
*  rhs = {endo inexog} = RHS regressors
*  no 0 or 1 at end of varlist means original varlist but after expansion of FV and TS vars
*  0 at the end of the name means the varlist after duplicates removed and collinearities/omitteds marked
*  1 means the same as 0 but after omitted vars dropped and extraneous FV operators "o", "b" and "n" removed.
*  0, 1 etc. also apply to _ct variables that are counts of these varlists
*  dofminus is large-sample adjustment (e.g., #fixed effects)
*  sdofminus is small-sample adjustment (e.g., #partialled-out regressors)

if c(version) < 12 & c(version) >= 9 {
* livreg2 Mata library.
* Ensure Mata library is indexed if new install.
* Not needed for Stata 12+ since ssc.ado does this when installing.
	capture mata: mata drop m_calckw()
	capture mata: mata drop m_omega()
	capture mata: mata drop ms_vcvorthog()
	capture mata: mata drop s_vkernel()
	capture mata: mata drop s_cdsy()
	mata: mata mlib index
}

*********************************************************************************
***************************** PARENT IVREG2 *************************************
****************** FORKS TO EXTERNAL IVREG2S IF CALLER < 11 *********************
*********************************************************************************

* Parent program, forks to versions as appropriate after version call
* Requires byable(onecall)
program define ivreg2, eclass byable(onecall) /* properties(svyj) */ sortpreserve
	local lversion 04.1.10

* Minimum of version 8 required for parent program (earliest ivreg2 is ivreg28)
	version 8

* Replay = no arguments before comma	
	if replay() {
* Call to ivreg2 will either be for version, in which case there should be no other arguments,
* or a postestimation call, in which case control should pass to main program.
		syntax [, VERsion * ]
		if "`version'"~="" & "`options'"=="" {
* Call to ivreg2 is for version
			di in gr "`lversion'"
			ereturn clear
			ereturn local version `lversion'
			exit	
		}
		else if "`version'"~="" & "`options'"~="" {
* Improper use of version option
di as err "invalid syntax - cannot combine version with other options"
			exit 198
		}
		else {
* Postestimation call, so put `options' macro (i.e. *) back into `0' macro with preceding comma
			local 0 `", `options'"'
		}
	}

* replay can't be combined with by
	if replay() & _by() {
di as err "invalid syntax - cannot use by with replay"
		exit 601
	}

* Handling of by. ivreg2x programs are byable(recall), so must set prefix for them.
	if _by() {
		local BY `"by `_byvars'`_byrc0':"'
	}

* If calling version is < 11, pass control to earlier version
* Note that this means calls from version 11.0 will not go to legacy version
* but will fail requirement of version 11.2 in main code.
	if _caller() < 11 {
		local ver = round(_caller())
		local ivreg2cmd ivreg2`ver'
* If replay, change e(cmd) macro to name of legacy ivreg2 before calling it, then change back
* Note by not allowed with replay; caught above so prefix not needed here.
		if replay() {
			ereturn local cmd "`ivreg2cmd'"
			`ivreg2cmd' `0'
			ereturn local cmd "ivreg2"

}
		else {
* If not replay, call legacy ivreg2 and then add macros
			`BY' `ivreg2cmd' `0'
			ereturn local cmd "ivreg2"
			ereturn local ivreg2cmd "`ivreg2cmd'"
			ereturn local version `lversion'
			ereturn local predict ivreg2_p
		}
		exit
	}

end
*********************************************************************************
*************************** END PARENT IVREG2 ***********************************
*********************************************************************************

********************* EXIT IF STATA VERSION < 11 ********************************

* When do file is loaded, exit here if Stata version calling program is < 11.
* Prevents loading of rest of program file (could cause earlier Statas to crash).

if c(stata_version) < 11 {
	exit
}

******************** END EXIT IF STATA VERSION < 11 *****************************

*********************************************************************************
***************** BEGIN MAIN IVREG2 ESTIMATION CODE *****************************
*********************************************************************************

* Main estimation program
program define ivreg211, eclass byable(recall) sortpreserve
	version 11.2
	local ranktestversion 01.3.02

local ivreg2cmd "ivreg211"			//  actual command name
	local ivreg2name "ivreg2"			//  name used in command line and for default naming of equations etc.
	local ranktestcmd "ranktest"

if replay() {
		syntax [, 												///
				FIRST FFIRST RF SFIRST							///
				dropfirst droprf dropsfirst						///
				Level(integer $S_level)							///
				NOHEader NOFOoter								///
				EForm(string) PLUS								///
				NOOMITTED vsquish noemptycells					///
				baselevels allbaselevels 						///
				VERsion											///
				]
		if "`version'" != "" & "`first'`ffirst'`rf'`noheader'`nofooter'`dropfirst'`droprf'`eform'`plus'" != "" {
			di as err "option version not allowed"
			error 198
		}
		if "`version'" != "" {
			di in gr "`lversion'"
			ereturn clear
			ereturn local version `lversion'
			exit
		}
		if `"`e(cmd)'"' != "ivreg2"  {
			error 301
		}
// Set display options
		local dispopt	eform(`eform') `noomitted' `vsquish' `noemptycells' `baselevels' `allbaselevels'

// On replay, set flag so saved eqns aren't dropped
		if "`e(firsteqs)'" != "" & "`dropfirst'" == "" {
			local savefirst "savefirst"
		}
		if "`e(rfeq)'" != "" & "`droprf'" == "" {
			local saverf "saverf"
		}
		if "`e(sfirsteq)'" != "" & "`dropsfirst'" == "" {
			local savesfirst "savesfirst"
		}
// On replay, re-display collinearities and duplicates messages
		DispCollinDups
	}
	else {
// MAIN CODE BLOCK
		local cmdline `ivreg2name' `0'		//  `0' rather than `*' in case of any "s in string

//  Confirm ranktest is installed (necessary component).
		checkversion_ranktest `ranktestcmd' `ranktestversion'

// Start parsing
		syntax [anything(name=0)] [if] [in] [aw fw pw iw/] [,		///
				NOID NOCOLLIN										///
				FIRST FFIRST SAVEFIRST SAVEFPrefix(name)			///
				RF SAVERF SAVERFPrefix(name)						///
				SFIRST SAVESFIRST SAVESFPrefix(name)				///
				SMall NOConstant 									///
				Robust CLuster(varlist) kiefer dkraay(integer 0)	///
				BW(string) kernel(string) center					///
				GMM GMM2s CUE										///
				LIML COVIV FULLER(real 0) Kclass(real 0)			///
				ORTHOG(string) ENDOGtest(string) REDundant(string)	///
				PARTIAL(string) FWL(string)							///
				Level(integer $S_level)								///
				NOHEader NOFOoter NOOUTput							///
				bvclean NOOMITTED omitted vsquish noemptycells		///
				baselevels allbaselevels 							///
				title(string) subtitle(string)						///
				DEPname(string) EForm(string) PLUS					///
				Tvar(varname) Ivar(varname)							///
				B0(string) SMATRIX(string) WMATRIX(string)			///
				sw psd0 psda useqr									///
				dofminus(integer 0) sdofminus(integer 0)			///
				NOPARTIALSMALL										///
				fvall fvsep											///
				]

// Parse after clearing any sreturn macros (can be left behind in Stata 11)
		sreturn clear
		ivparse `0',	ivreg2name(`ivreg2name')		///  needed for some options
						partial(`partial')				///
						fwl(`fwl')						///  legacy option
						orthog(`orthog')				///
						endogtest(`endogtest')			///
						redundant(`redundant')			///
						depname(`depname')				///
						`robust'						///
						cluster(`cluster')				///
						bw(`bw')						///
						kernel(`kernel')				///
						dkraay(`dkraay')				///
						`center'						///
						`kiefer'						///
						`sw'							///
						`noconstant'					///
						tvar(`tvar')					///
						ivar(`ivar')					///
						`gmm2s'							///
						`gmm'							///  legacy option, produces error message
						`cue'							///
						`liml'							///
						fuller(`fuller')				///
						kclass(`kclass')				///
						b0(`b0')						///
						wmatrix(`wmatrix')				///
						`noid'							///
						`savefirst'						///
						savefprefix(`savefprefix')		///
						`saverf'						///
						saverfprefix(`saverfprefix')	///
						`savesfirst'					///
						savesfprefix(`savesfprefix')	///
						dofminus(`dofminus')			///
						`psd0'							///
						`psda'							///
						`nocollin'						///
						`useqr'							///
						`bvclean'						///
						eform(`eform')					///
						`noomitted'						///
						`vsquish'						///
						`noemptycells'					///
						`baselevels'					///
						`allbaselevels'

// varlists are unexpanded; may be empty
		local lhs			`s(lhs)'
		local depname		`s(depname)'
		local endo			`s(endo)'
		local inexog		`s(inexog)'
		local exexog		`s(exexog)'
		local partial		`s(partial)'
		local cons			=s(cons)
		local partialcons	=s(partialcons)
		local tvar			`s(tvar)'
		local ivar			`s(ivar)'
		local tdelta		`s(tdelta)'
		local tsops			=s(tsops)
		local fvops			=s(fvops)
		local robust		`s(robust)'
		local cluster		`s(cluster)'
		local bw			=`s(bw)'				//  arrives as string but return now as number
		local bwopt			`s(bwopt)'
		local kernel		`s(kernel)'				//  also used as flag for HAC estimation
		local center		=`s(center)'			//  arrives as string but now boolean
		local kclassopt		`s(kclassopt)'
		local fulleropt		`s(fulleropt)'
		local liml			`s(liml)'
		local noid			`s(noid)'				//  can also be triggered by b0(.) option
		local useqr			=`s(useqr)'				//  arrives as string but now boolean; nocollin=>useqr
		local savefirst		`s(savefirst)'
		local savefprefix	`s(savefprefix)'
		local saverf		`s(saverf)'
		local saverfprefix	`s(saverfprefix)'
		local savesfirst	`s(savesfirst)'
		local savesfprefix	`s(savesfprefix)'
		local psd			`s(psd)'				//  triggered by psd0 or psda
		local dofmopt		`s(dofmopt)'
		local bvclean		=`s(bvclean)'			//  arrives as string but return now as boolean
		local dispopt		`s(dispopt)'

// Can now tsset; sortpreserve will restore sort after exit
		if `tsops' | "`kernel'"~="" {
			cap tsset								//  restores sort if tsset or xtset but sort disrupted
			if _rc>0 {
				tsset `ivar' `tvar'
			}
		}

***********************************************************

// Weights
// fweight and aweight accepted as is
// iweight not allowed with robust or gmm and requires a trap below when used with summarize
// pweight is equivalent to aweight + robust
// Since we subsequently work with wvar, tsrevar of weight vars in weight `exp' not needed.

tempvar wvar
		if "`weight'" == "fweight" | "`weight'"=="aweight" {
			local wtexp `"[`weight'=`exp']"'
			qui gen double `wvar'=`exp'
		}
		if "`weight'" == "fweight" & "`kernel'" !="" {
			di in red "fweights not allowed (data are -tsset-)"
			exit 101
		}
		if "`weight'" == "fweight" & "`sw'" != "" {
			di in red "fweights currently not supported with -sw- option"
			exit 101
		}
		if "`weight'" == "iweight" {
			if "`robust'`cluster'`gmm2s'`kernel'" !="" {
				di in red "iweights not allowed with robust or gmm"
				exit 101
			}
			else {
				local wtexp `"[`weight'=`exp']"'
				qui gen double `wvar'=`exp'
			}
		}
		if "`weight'" == "pweight" {
			local wtexp `"[aweight=`exp']"'
			qui gen double `wvar'=`exp'
			local robust "robust"
		}
		if "`weight'" == "" {
* If no weights, define neutral weight variable
			qui gen byte `wvar'=1
		}

********************************************************************************
// markout sample
// include `tvar' to limit sample to where tvar is available, but only if TS operators used
		marksample touse
		if `tsops' {
			markout `touse' `lhs' `inexog' `exexog' `endo' `cluster' `tvar', strok
		}
		else {
			markout `touse' `lhs' `inexog' `exexog' `endo' `cluster', strok
		}

********************************************************************************
// weight factor and sample size
// Every time a weight is used, must multiply by scalar wf ("weight factor")
// wf=1 for no weights, fw and iw, wf = scalar that normalizes sum to be N if aw or pw

sum `wvar' if `touse' `wtexp', meanonly
// Weight statement
		if "`weight'" ~= "" {
di in gr "(sum of wgt is " %14.4e `r(sum_w)' ")"
		}
		if "`weight'"=="" | "`weight'"=="fweight" | "`weight'"=="iweight" {
// Effective number of observations is sum of weight variable.
// If weight is "", weight var must be column of ones and N is number of rows
			local wf=1
			local N=r(sum_w)
		}
		else if "`weight'"=="aweight" | "`weight'"=="pweight" {
			local wf=r(N)/r(sum_w)
			local N=r(N)
		}
		else {
// Should never reach here
di as err "ivreg2 error - misspecified weights"
			exit 198
		}
		if `N'==0 {
di as err "no observations"
			exit 2000
		}

***************************************************************
// Time-series data
// tindex used by Mata code so that ts operators work correctly

tempvar tindex
		qui gen `tindex'=1 if `touse'
		qui replace `tindex'=sum(`tindex') if `touse'

if `tsops' | "`kernel'"~="" {
// Report gaps in data
			tsreport if `touse', panel
			if `r(N_gaps)' != 0 {
di as text "Warning: time variable " as res "`tvar'" as text " has "		///
				as res "`r(N_gaps)'" as text " gap(s) in relevant range"
			}
// Set local macro T and check that bw < (T-1)
			sum `tvar' if `touse', meanonly
			local T = r(max)-r(min) + 1
			local T1 = `T' - 1
			if (`bw' > (`T1'/`tdelta')) {
di as err "invalid bandwidth in option bw() - cannot exceed timespan of data"
				exit 198
			}
		}

// kiefer VCV = kernel(tru) bw(T) and no robust with tsset data
		if "`kiefer'" ~= "" {
			local bw	=`T'
		}

*********** Column of ones for constant set up here **************

if "`noconstant'"=="" {
// If macro not created, automatically omitted.
			tempvar ones
			qui gen byte `ones' = 1 if `touse'
		}

************* Varlists, FV varlists, duplicates *****************
// Varlists come in 4 versions, e.g., for inexog:
// (a) inexog = full list of original expanded vnames; may have duplicates
// (b) inexog0 = as with inexog with duplicates removed but RETAINING base/omitted/etc. varnames
// (c) inexog1 = as with inexog0 but WITHOUT base/omitted/etc.
// (d) fv_inexog1 = corresponding list with temp vars minus base/omitted/etc., duplicates, collinearities etc.
// Varlists (c) and (d) are definitive, i.e., have the variables actually used in the estimation.

// Create consistent expanded varlists.
// "Consistent" means base vars for FVs must be consistent
// hence default rhs=endo+inexog is expanded as one.
// fvall: overrides, endo+inexog+exexog expanded as one
// fvsep: overrides, endo, inexog and exexog expanded separately
// NB: expanding endo+inexog+exexog is dangerous because
//     fvexpand can zap a list in case of overlap
//     e.g. fvexpand mpg + i(1/4).rep78 + i5.rep78 
//             => mpg 1b.rep78 2.rep78 3.rep78 4.rep78 5.rep78
//     but fvexpand mpg + i.rep78 + i5.rep78
//				=> mpg 5.rep78

CheckDupsCollin,							///
							lhs(`lhs')				///
							endo(`endo')			///
							inexog(`inexog')		///
							exexog(`exexog')		///
							partial(`partial')		///
							orthog(`orthog')		///
							endogtest(`endogtest')	///
							redundant(`redundant')	///
							touse(`touse')			///
							wvar(`wvar')			///
							wf(`wf')				///
							`noconstant'			///
							`nocollin'				///
							`fvall'					///
							`fvsep'

// Replace basic varlists and create "0" versions of varlists
		foreach vl in lhs endo inexog exexog partial orthog endogtest redundant {
			local `vl'		`s(`vl')'
			local `vl'0		`s(`vl'0)'
		}
		local dups		`s(dups)'
		local collin	`s(collin)'
		local ecollin	`s(ecollin)'

// Create "1" and fv versions of varlists
		foreach vl in lhs endo inexog exexog partial orthog endogtest redundant {
			foreach var of local `vl'0 {						//  var-by-var so that fvrevar doesn't decide on base etc.
				_ms_parse_parts `var'
				if ~`r(omit)' {									//  create temp var only if not omitted
					fvrevar `var'	if `touse'
					local `vl'1		``vl'1' `var'
					local fv_`vl'1	`fv_`vl'1' `r(varlist)'
				}
			}
			local `vl'1		: list retokenize `vl'1
			local fv_`vl'1	: list retokenize fv_`vl'1
		}

// Check that LHS expanded to a single variable
		local wrongvars_ct	: word count `lhs'
		if `wrongvars_ct' > 1 {
di as err "multiple dependent variables specified: `lhs'"
			error 198
		}

// Check that option varlists are compatible with main varlists
// orthog()
		local wrongvars		: list orthog1 - inexog1
		local wrongvars		: list wrongvars - exexog1
		local wrongvars_ct	: word count `wrongvars'
		if `wrongvars_ct' {
di as err "Error: `wrongvars' listed in orthog() but does not appear as exogenous." 
			error 198
		}
// endog()
		local wrongvars		: list endogtest1 - endo1
		local wrongvars_ct	: word count `wrongvars'
		if `wrongvars_ct' {
di as err "Error: `wrongvars' listed in endog() but does not appear as endogenous." 
			error 198
		}
// redundant()
		local wrongvars		: list redundant1 - exexog1
		local wrongvars_ct	: word count `wrongvars'
		if `wrongvars_ct' {
di as err "Error: `wrongvars' listed in redundant() but does not appear as exogenous." 
			error 198
		}

// And create allnames macros
		local allnames		`lhs' `endo' `inexog' `exexog'
		local allnames0		`lhs0' `endo0' `inexog0' `exexog0'
		local allnames1		`lhs1' `endo1' `inexog1' `exexog1'
		local fv_allnames1	`fv_lhs1' `fv_endo1' `fv_inexog1' `fv_exexog1'

// *************** Partial-out block ************** //

// `partial' has all to be partialled out except for constant
		if "`partial1'" != "" | `partialcons'==1 {
			preserve

// Remove partial0 from inexog0.
// Remove partial1 from inexog1.
			local inexog0		: list inexog0 - partial0
			local inexog1		: list inexog1 - partial1
			local fv_inexog1	: list fv_inexog1 - fv_partial1

// Check that cluster, weight, tvar or ivar variables won't be transformed
// Use allnames1 (expanded varlist)
			if "`cluster'"~="" {
				local pvarcheck : list cluster in allnames1
				if `pvarcheck' {
di in r "Error: cannot use cluster variable `cluster' as dependent variable, regressor or IV"
di in r "       in combination with -partial- option." 
				error 198
				}
			}
			if "`tvar'"~="" {
				local pvarcheck : list tvar in allnames1
				if `pvarcheck' {
di in r "Error: cannot use time variable `tvar' as dependent variable, regressor or IV"
di in r "       in combination with -partial- option." 
				error 198
				}
			}
			if "`ivar'"~="" {
				local pvarcheck : list ivar in allnames1
				if `pvarcheck' {
di in r "Error: cannot use panel variable `ivar' as dependent variable, regressor or IV"
di in r "       in combination with -partial- option." 
				error 198
				}
			}
			if "`wtexp'"~="" {
				tokenize `exp', parse("*/()+-^&|~")
				local wvartokens `*'
				local nwvarnames : list allnames1 - wvartokens
				local wvarnames  : list allnames1 - nwvarnames
				if "`wvarnames'"~="" {
di in r "Error: cannot use weight variables as dependent variable, regressor or IV"
di in r "       in combination with -partial- option." 
				error 198
				}
			}
// Partial out
// But first replace everything with doubles
			recast double `fv_lhs1' `fv_endo1' `fv_inexog1' `fv_exexog1' `fv_partial1'
			mata: s_partial	("`fv_lhs1'",			///
							"`fv_endo1'",			///
							"`fv_inexog1'",			///
							"`fv_exexog1'",			///
							"`fv_partial1'",		///
							"`touse'",				///
							"`weight'",				///
							"`wvar'",				///
							`wf',					///
							`N',					///
							`cons')

local partial_ct : word count `partial1'
// Constant is partialled out, unless nocons already specified in the first place
			capture drop `ones'
			local ones ""
			if "`noconstant'" == "" {
// partial_ct used for small-sample adjustment to regression F-stat
				local partial_ct = `partial_ct' + 1
				local noconstant "noconstant"
				local cons 0
			}
		}
		else {
// Set count of partial vars to zero if option not used
			local partial_ct 0
			local partialcons 0
		}
// Add partial_ct to small dof adjustment sdofminus
		if "`nopartialsmall'"=="" {
			local sdofminus = `sdofminus'+`partial_ct'
		}

*********************************************

local rhs0			`endo0' `inexog0'				//  needed for display of omitted/base/etc.
		local rhs1			`endo1' `inexog1'
		local insts1		`exexog1' `inexog1'
		local fv_insts1		`fv_exexog1' `fv_inexog1'
		local fv_rhs1		`fv_endo1' `fv_inexog1'
		local rhs0_ct		: word count `rhs0'				//  needed for display of omitted/base/etc.
		local rhs1_ct		: word count `fv_rhs1'
		local iv1_ct		: word count `fv_insts1'
		local endo1_ct		: word count `fv_endo1'
		local exex1_ct		: word count `fv_exexog1'
		local endoexex1_c	: word count `fv_endo1' `fv_exexog1'
		local inexog1_ct	: word count `fv_inexog1'

// Counts modified to include constant if appropriate
		local rhs1_ct = `rhs1_ct' + `cons'
		local rhs0_ct = `rhs0_ct' + `cons'					//  needed for display of omitted/base/etc.
		local iv1_ct  = `iv1_ct' + `cons'

// Column/row names for matrices b, V, S, etc.
		local cnb0	`endo0' `inexog0'				//  including omitted
		local cnb1	`endo1' `inexog1'				//  excluding omitted
		local cnZ0	`exexog0' `inexog0'				//  excluding omitted
		local cnZ1	`exexog1' `inexog1'				//  excluding omitted
		if `cons' {
			local cnb0	"`cnb0' _cons"
			local cnb1	"`cnb1' _cons"
			local cnZ0	"`cnZ0' _cons"
			local cnZ1	"`cnZ1' _cons"
		}

*********************************************
// Remaining checks: variable counts, col/row names of b0, smatrix, wmatrix
		CheckMisc,							///
					rhs1_ct(`rhs1_ct')		///
					iv1_ct(`iv1_ct')		///
					bvector(`b0')			///
					smatrix(`smatrix')		///
					wmatrix(`wmatrix')		///
					cnb1(`cnb1')			///
					cnZ1(`cnZ1')

if "`b0'"~="" {
			tempname b0						//  so we can overwrite without changing original user matrix
			mat `b0'		= r(b0)
		}
		if "`smatrix'"~="" {
			tempname S0
			mat `S0'		= r(S0)
		}
		if "`wmatrix'"~="" {
			tempname wmatrix				//  so we can overwrite without changing original user matrix
			mat `wmatrix'	= r(W0)
		}
		
*************** Commonly used matrices ****************
		tempname YY yy yyc
		tempname XX X1X1 X2X2 X1Z X1Z1 XZ Xy
		tempname ZZ Z1Z1 Z2Z2 Z1Z2 Z1X2 Zy ZY Z2y Z2Y
		tempname XXinv X2X2inv ZZinv XPZXinv
		tempname rankxx rankzz condxx condzz

// use fv_ varlists
		mata: s_crossprods	("`fv_lhs1'",			///
							"`fv_endo1'",			///
							"`fv_inexog1' `ones'",	///
							"`fv_exexog1'",			///
							"`touse'",				///
							"`weight'",				///
							"`wvar'",				///
							`wf',					///
							`N')
		mat `XX'		=r(XX)
		mat `X1X1'		=r(X1X1)
		mat `X1Z'		=r(X1Z)
		mat `ZZ'		=r(ZZ)
		mat `Z2Z2'		=r(Z2Z2)
		mat `Z1Z2'		=r(Z1Z2)
		mat `XZ'		=r(XZ)
		mat `Xy'		=r(Xy)
		mat `Zy'		=r(Zy)
		mat `YY'		=r(YY)
		scalar `yy'		=r(yy)
		scalar `yyc'	=r(yyc)
		mat `ZY'		=r(ZY)
		mat `Z2y'		=r(Z2y)
		mat `Z2Y'		=r(Z2Y)
		mat `XXinv'		=r(XXinv)
		mat `ZZinv'		=r(ZZinv)
		mat `XPZXinv'	=r(XPZXinv)
		scalar `condxx'	=r(condxx)
		scalar `condzz'	=r(condzz)

scalar `rankzz'	= rowsof(`ZZinv') - diag0cnt(`ZZinv')
		scalar `rankxx'	= rowsof(`XXinv') - diag0cnt(`XXinv')
		local overid	= `rankzz' - `rankxx'

********** CLUSTER SETUP **********************************************

* Mata code requires data are sorted on (1) the first var cluster if there
* is only one cluster var; (2) on the 3rd and then 1st if two-way clustering,
* unless (3) two-way clustering is combined with kernel option, in which case
* the data are tsset and sorted on panel id (first cluster variable) and time
* id (second cluster variable).
* Second cluster var is optional and requires an identifier numbered 1..N_clust2,
* unless combined with kernel option, in which case it's the time variable.
* Third cluster var is the intersection of 1 and 2, unless combined with kernel
* opt, in which case it's unnecessary.
* Sorting on "cluster3 cluster1" means that in Mata, panelsetup works for
* both, since cluster1 nests cluster3.
* Note that it is possible to cluster on time but not panel, in which case
* cluster1 is time, cluster2 is empty and data are sorted on panel-time.
* Note also that if data are sorted here but happen to be tsset, will need
* to be re-tsset after estimation code concludes.

// No cluster options or only 1-way clustering
// but for Mata and other purposes, set N_clust vars =0
		local N_clust=0
		local N_clust1=0
		local N_clust2=0
		if "`cluster'"!="" {
			local clopt "cluster(`cluster')"
			tokenize `cluster'
			local cluster1 "`1'"
			local cluster2 "`2'"
			if "`kernel'"~="" {
* kernel requires either that cluster1 is time var and cluster2 is empty
* or that cluster1 is panel var and cluster2 is time var.
* Either way, data must be tsset and sorted for panel data.
				if "`cluster2'"~="" {
* Allow backwards order
					if "`cluster1'"=="`tvar'" & "`cluster2'"=="`ivar'" {
						local cluster1 "`2'"
						local cluster2 "`1'"
					}
					if "`cluster1'"~="`ivar'" | "`cluster2'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset panel & time vars."
di as err "       tsset panel var=`ivar'; tsset time var=`tvar'; cluster vars=`cluster1',`cluster2'"
						exit 198
					}
				}
				else {
					if "`cluster1'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset time variable."
di as err "       tsset time var=`tvar'; cluster var=`cluster1'"
						exit 198
					}
				}
			}
* Simple way to get quick count of 1st cluster variable without disrupting sort
* clusterid1 is numbered 1.._Nclust1.
			tempvar clusterid1
			qui egen `clusterid1'=group(`cluster1') if `touse'
			sum `clusterid1' if `touse', meanonly
			if "`cluster2'"=="" {
				local N_clust=r(max)
				local N_clust1=`N_clust'
				if "`kernel'"=="" {
* Single level of clustering and no kernel-robust, so sort on single cluster var.
* kernel-robust already sorted via tsset.
					sort `cluster1'
				}
			}
			else {
				local N_clust1=r(max)
				if "`kernel'"=="" {
					tempvar clusterid2 clusterid3
* New cluster id vars are numbered 1..N_clust2 and 1..N_clust3
					qui egen `clusterid2'=group(`cluster2') if `touse'
					qui egen `clusterid3'=group(`cluster1' `cluster2') if `touse'
* Two levels of clustering and no kernel-robust, so sort on cluster3/nested in/cluster1
* kernel-robust already sorted via tsset.
					sort `clusterid3' `cluster1'
					sum `clusterid2' if `touse', meanonly
					local N_clust2=r(max)
				}
				else {
* Need to create this only to count the number of clusters
					tempvar clusterid2
					qui egen `clusterid2'=group(`cluster2') if `touse'
					sum `clusterid2' if `touse', meanonly
					local N_clust2=r(max)
* Now replace with original variable
					local clusterid2 `cluster2'
				}

local N_clust=min(`N_clust1',`N_clust2')

}		// end 2-way cluster block
		}		// end cluster block

************************************************************************************************

tempname b W S V beta lambda j jp rss mss rmse sigmasq rankV rankS
		tempname arubin arubinp arubin_lin arubin_linp
		tempname r2 r2_a r2u r2c F Fp Fdf2 ivest

tempvar resid
		qui gen double `resid'=.
		
*******************************************************************************************
* LIML
*******************************************************************************************

if "`liml'`kclassopt'"~="" {

mata: s_liml(	"`ZZ'",								///
							"`XX'",								///
							"`XZ'",								///
							"`Zy'",								///
							"`Z2Z2'",							///
							"`YY'",								///
							"`ZY'",								///
							"`Z2Y'",							///
							"`Xy'",								///
							"`ZZinv'",							///
							"`fv_lhs1'", 						///
							"`fv_lhs1' `fv_endo1'",				///
							"`resid'",							///
							"`fv_endo1' `fv_inexog1' `ones'",	///
							"`fv_endo1'",						///
							"`fv_exexog1' `fv_inexog1' `ones'",	///
							"`fv_exexog1'",						///
							"`fv_inexog1' `ones'",				///
							`fuller',							///
							`kclass',							///
							"`coviv'",							///
							"`touse'",							///
							"`weight'",							///
							"`wvar'",							///
							`wf',								///
							`N',								///
							"`robust'",							///
							"`clusterid1'",						///
							"`clusterid2'",						///
							"`clusterid3'",						///
							`bw',								///
							"`kernel'",							///
							"`sw'",								///
							"`psd'",							///
							"`ivar'",							///
							"`tvar'",							///
							"`tindex'",							///
							`tdelta',							///
							`center',							///
							`dofminus',							///
							`useqr')

mat `b'=r(beta)
			mat `S'=r(S)
			mat `V'=r(V)
			scalar `lambda'=r(lambda)
			local kclass=r(kclass)
			scalar `j'=r(j)
			scalar `rss'=r(rss)
			scalar `sigmasq'=r(sigmasq)
			scalar `rankV'=r(rankV)
			scalar `rankS'=r(rankS)

scalar `arubin'=(`N'-`dofminus')*ln(`lambda')
			scalar `arubin_lin'=(`N'-`dofminus')*(`lambda'-1)
			
// collinearities can cause LIML to generate (spurious) OLS results
			if "`nocollin'"~="" & `kclass'<1e-8 {
di as err "warning: k=1 in LIML estimation; results equivalent to OLS;"
di as err "         may be caused by collinearities"
			}
		}

*******************************************************************************************
* OLS, IV and 2SGMM.  Also enter to get CUE starting values.
************************************************************************************************

if "`liml'`kclassopt'`b0'"=="" {

* Call to s_gmm1s to do 1st-step GMM.
* If W or S supplied, calculates GMM beta and residuals
* If none of the above supplied, calculates GMM beta using default IV weighting matrix and residuals
* Block not entered if b0 is provided.

* 1-step GMM is efficient and V/J/Sargan can be returned if:
*   - estimator is IV, W is known and S can be calculated from 1st-step residuals
*   - S is provided (and W is NOT) so W=inv(S) and beta can be calculated using W
* 1-step GMM is inefficient if:
*   - non-iid VCE is requested
*   - W is provided

local effic1s =		(													///
										"`gmm2s'`robust'`cluster'`kernel'"==""		///
									|	("`smatrix'"~="" & "`wmatrix'"=="")			///
								)

// use fv_ varlists
			mata: s_gmm1s(	"`ZZ'",									///
							"`XX'",									///
							"`XZ'",									///
							"`Zy'",									///
							"`ZZinv'",								///
							"`fv_lhs1'", 							///
							"`resid'",								///
							"`fv_endo1' `fv_inexog1' `ones'",		///
							"`fv_exexog1' `fv_inexog1' `ones'",		///
							"`touse'",								///
							"`weight'",								///
							"`wvar'",								///
							`wf',									///
							`N',									///
							"`wmatrix'",							///
							"`S0'",									///
							`dofminus',								///
							`effic1s',								///
							`overid',								///
							`useqr')
			mat `b'=r(beta)
			mat `W'=r(W)

* If 1st-step is efficient, save remaining results and we're done
			if `effic1s' {
				mat `V'=r(V)
				mat `S'=r(S)
				scalar `j'=r(j)
				scalar `rss'=r(rss)
				scalar `sigmasq'=r(sigmasq)
				scalar `rankV'=r(rankV)
				scalar `rankS'=r(rankS)
			}
			else {
* ...we're not done - do inefficient or 2-step efficient GMM

* Pick up matrix left by s_gmm1s(.)
				tempname QXZ_W_QZX
				mat `QXZ_W_QZX'=r(QXZ_W_QZX)

* Block calls s_omega to get cov matrix of orthog conditions, if not supplied
				if "`smatrix'"~="" {
					mat `S'=`S0'
				}
				else {
	
* NB: xtivreg2 calls ivreg2 with data sorted on ivar and optionally tvar.
*     Stock-Watson adjustment -sw- assumes data are sorted on ivar.  Checked at start of ivreg2.
	
* call abw code if bw() is defined and bw(auto) selected
					if `bw' != 0 {
						if `bw' == -1 {
							tempvar abwtouse
							gen byte `abwtouse' = (`resid' < .) 
							abw `resid' `exexog1' `inexog1' `abwtouse', 	/*
								*/	tindex(`tindex') nobs(`N') tobs(`T') noconstant kernel(`kernel')
							local bw `r(abw)'
							local bwopt "bw(`bw')"
							local bwchoice "`r(bwchoice)'"
						}
					}
* S covariance matrix of orthogonality conditions
// use fv_ varlists
					mata: s_omega(	"`ZZ'",									///
									"`resid'",								///
									"`fv_exexog1' `fv_inexog1' `ones'",		///
									"`touse'",								///
									"`weight'",								///
									"`wvar'",								///
									`wf',									///
									`N',									///
									"`robust'",								///
									"`clusterid1'",							///
									"`clusterid2'",							///
									"`clusterid3'",							///
									`bw',									///
									"`kernel'",								///
									"`sw'",									///
									"`psd'",								///
									"`ivar'",								///
									"`tvar'",								///
									"`tindex'",								///
									`tdelta',								///
									`center',								///
									`dofminus')
					mat `S'=r(S)
				}
	
* By this point: `b' has 1st-step inefficient beta
*                `resid' has resids from the above beta
*                `S' has vcv of orthog conditions using either `resid' or user-supplied `S0'
*                 `QXZ_W_QZX' was calculated in s_gmm1s(.) for use in s_iegmm(.)

* Inefficient IV.  S, W and b were already calculated above.
				if "`gmm2s'"=="" & "`robust'`cluster'`kernel'"~="" {
					mata: s_iegmm(	"`ZZ'",								///
									"`XX'",								///
									"`XZ'",								///
									"`Zy'",								///
									"`QXZ_W_QZX'",						///
									"`fv_lhs1'", 						///
									"`resid'",							///
									"`fv_endo1' `fv_inexog1' `ones'",	///
									"`fv_exexog1' `fv_inexog1' `ones'",	///
									"`touse'",							///
									"`weight'",							///
									"`wvar'",							///
									`wf',								///
									`N',								///
									"`W'",								///
									"`S'",								///
									"`b'",								///
									`dofminus',							///
									`overid',							///
									`useqr')
					}

* 2-step efficient GMM.  S calculated above, b and W will be updated.
				if "`gmm2s'"~="" {
					mata: s_egmm(	"`ZZ'",								///
									"`XX'",								///
									"`XZ'",								///
									"`Zy'",								///
									"`ZZinv'",							///
									"`fv_lhs1'", 						///
									"`resid'",							///
									"`fv_endo1' `fv_inexog1' `ones'",	///
									"`fv_exexog1' `fv_inexog1' `ones'",	///
									"`touse'",							///
									"`weight'",							///
									"`wvar'",							///
									`wf',								///
									`N',								///
									"`S'",								///
									`dofminus',							///
									`overid',							///
									`useqr')
					mat `b'=r(beta)
					mat `W'=r(W)
				}
	
				mat `V'=r(V)
				scalar `j'=r(j)
				scalar `rss'=r(rss)
				scalar `sigmasq'=r(sigmasq)
				scalar `rankV'=r(rankV)
				scalar `rankS'=r(rankS)
			}
* Finished with non-CUE/LIML block
		}

***************************************************************************************
* Block for cue gmm
*******************************************************************************************
		if "`cue'`b0'" != "" {

* s_gmmcue is passed initial b from IV/2-step GMM block above
* OR user-supplied b0 for evaluation of CUE obj function at b0
			mata: s_gmmcue(	"`ZZ'",								///
							"`XZ'",								///
							"`fv_lhs1'", 						///
							"`resid'",							///
							"`fv_endo1' `fv_inexog1' `ones'",	///
							"`fv_exexog1' `fv_inexog1' `ones'",	///
							"`touse'",							///
							"`weight'",							///
							"`wvar'",							///
							`wf',								///
							`N',								///
							"`robust'",							///
							"`clusterid1'",						///
							"`clusterid2'",						///
							"`clusterid3'",						///
							`bw',								///
							"`kernel'",							///
							"`sw'",								///
							"`psd'",							///
							"`ivar'",							///
							"`tvar'",							///
							"`tindex'",							///
							`tdelta',							///
							"`b'",								///
							"`b0'",								///
							`center',							///
							`dofminus',							///
							`useqr')

mat `b'=r(beta)
			mat `S'=r(S)
			mat `W'=r(W)
			mat `V'=r(V)
			scalar `j'=r(j)
			scalar `rss'=r(rss)
			scalar `sigmasq'=r(sigmasq)
			scalar `rankV'=r(rankV)
			scalar `rankS'=r(rankS)

}

****************************************************************
* Done with estimation blocks
****************************************************************

mat colnames `b' = `cnb1'
		mat colnames `V' = `cnb1'
		mat rownames `V' = `cnb1'
		mat colnames `S' = `cnZ1'
		mat rownames `S' = `cnZ1'
* No W matrix for LIML or kclass
		capture mat colnames `W' = `cnZ1'
		capture mat rownames `W' = `cnZ1'

*******************************************************************************************
* RSS, counts, dofs, F-stat, small-sample corrections
*******************************************************************************************

// rankxx = rhs1_ct except if nocollin
// rankzz = iv1_ct except if nocollin
// nocollin means count may exceed rank (because of dropped vars), so rank #s foolproof
		scalar `rmse'=sqrt(`sigmasq')
		if "`noconstant'"=="" {
			scalar `mss'=`yyc' - `rss'
		}
		else {
			scalar `mss'=`yy' - `rss'
		}

local Fdf1 = `rankxx' - `cons'
		local df_m = `rankxx' - `cons' + (`sdofminus'-`partialcons')

* Residual dof
		if "`cluster'"=="" {
* Use int(`N') because of non-integer N with iweights, and also because of
* possible numeric imprecision with N returned by above.
			local df_r = int(`N') - `rankxx' - `dofminus' - `sdofminus'
		}
		else {
* To match Stata, subtract 1
			local df_r = `N_clust' - 1
		}

* Sargan-Hansen J dof and p-value
* df=0 doesn't guarantee j=0 since can be call to get value of CUE obj fn
		local jdf = `rankzz' - `rankxx'
		if `jdf' == 0 & "`b0'"=="" {
			scalar `j' = 0
		}
		else {
			scalar `jp' = chiprob(`jdf',`j')
		}
		if "`liml'"~="" {
			scalar `arubinp' = chiprob(`jdf',`arubin')
			scalar `arubin_linp' = chiprob(`jdf',`arubin_lin')
		}

* Small sample corrections for var-cov matrix.
* If robust, the finite sample correction is N/(N-K), and with no small
* we change this to 1 (a la Davidson & MacKinnon 1993, p. 554, HC0).
* If cluster, the finite sample correction is (N-1)/(N-K)*M/(M-1), and with no small
* we change this to 1 (a la Wooldridge 2002, p. 193), where M=number of clusters.

if "`small'" != "" {
			if "`cluster'"=="" {
				matrix `V'=`V'*(`N'-`dofminus')/(`N'-`rankxx'-`dofminus'-`sdofminus')
			}
			else {
				matrix `V'=`V'*(`N'-1)/(`N'-`rankxx'-`sdofminus')		///
							* `N_clust'/(`N_clust'-1)
			}
			scalar `sigmasq'=`rss'/(`N'-`rankxx'-`dofminus'-`sdofminus')
			scalar `rmse'=sqrt(`sigmasq')
		}

scalar `r2u'=1-`rss'/`yy'
		scalar `r2c'=1-`rss'/`yyc'
		if "`noconstant'"=="" {
			scalar `r2'=`r2c'
			scalar `r2_a'=1-(1-`r2')*(`N'-1)/(`N'-`rankxx'-`dofminus'-`sdofminus')
		}
		else {
			scalar `r2'=`r2u'
			scalar `r2_a'=1-(1-`r2')*`N'/(`N'-`rankxx'-`dofminus'-`sdofminus')
		}
* `N' is rounded down to nearest integer if iweights are used.
* If aw, pw or fw, should already be integer but use round in case of numerical imprecision.
		local N=int(`N')

* Fstat
* To get it to match Stata's, must post separately with dofs and then do F stat by hand
*   in case weights generate non-integer obs and dofs
* Create copies so they can be posted
		tempname FB FV
		mat `FB'=`b'
		mat `FV'=`V'
		capture ereturn post `FB' `FV'
* If the cov matrix wasn't positive definite, the post fails with error code 506
		local rc = _rc
		if `rc' != 506 {
* Strip out omitted/base/etc. vars from RHS list
			ivreg2_fvstrip `rhs1', dropomit
			capture test `r(varlist)'
			if "`small'" == "" {
				if "`cluster'"=="" {
					capture scalar `F' = r(chi2)/`Fdf1' * `df_r'/(`N'-`dofminus')
				}
				else {
* sdofminus used here so that F-stat matches test stat from regression with no partial and small
					capture scalar `F' =	r(chi2)/`Fdf1' * 					///
											(`N_clust'-1)/`N_clust' *			///
											(`N'-`rankxx'-`sdofminus')/(`N'-1)
				}
			}
			else {
				capture scalar `F' = r(chi2)/`Fdf1'
			}
			capture scalar `Fp'=Ftail(`Fdf1',`df_r',`F')
			capture scalar `Fdf2'=`df_r'
		}

* If j==. or vcv wasn't full rank, then vcv problems and F is meaningless
		if `j' == . | `rc'==506 {
			scalar `F' = .
			scalar `Fp' = .
		}

* End of counts, dofs, F-stat, small sample corrections

********************************************************************************************
* Reduced form and first stage regression options
*******************************************************************************************
* Relies on proper count of (non-collinear) IVs generated earlier.
* Note that nocons option + constant in instrument list means first-stage
* regressions are reported with nocons option.  First-stage F-stat therefore
* correctly includes the constant as an explanatory variable.

if "`sfirst'`savesfirst'`rf'`saverf'`first'`ffirst'`savefirst'" != "" & (`endo1_ct' > 0) {

* Restore original order if changed for mata code above
			capture tsset

local sdofmopt = "sdofminus(`sdofminus')"
// Need to create Stata placeholders for Mata code so that Stata time-series operators can work on them
// fres1 is Nx1
// endo1_hat is NxK1
// fsresall is Nx(K1+1) (used for full system)
			tempname fsres1
			qui gen double `fsres1'=.
			local fsresall `fsres1'
			foreach x of local fv_endo1 {
				tempname fsres
				qui gen double `fsres'=.
				local fsresall "`fsresall' `fsres'"
			}

// mata code requires sorting on cluster 3 / cluster 1 (if 2-way) or cluster 1 (if one-way)
			if "`cluster'"!="" {
					sort `clusterid3' `cluster1'
			}
			mata: s_ffirst(	"`ZZ'",							///
							"`XX'",							///
							"`XZ'",							///
							"`ZY'",							///
							"`ZZinv'",						///
							"`XXinv'",						///
							"`XPZXinv'",					///
							"`Z2Z2'",						///
							"`Z1Z2'",						///
							"`Z2y'",						///
							"`fsres1'",						/// Nx1
							"`fsresall'",					/// Nx(K1+1)
							"`fv_lhs1'",					///
							"`fv_endo1'",					///
							"`fv_inexog1' `ones'",			///
							"`fv_exexog1'",					///
							"`touse'",						///
							"`weight'",						///
							"`wvar'",						///
							`wf',							///
							`N',							///
							`N_clust',						///
							"`robust'",						///
							"`clusterid1'",					///
							"`clusterid2'",					///
							"`clusterid3'",					///
							`bw',							///
							"`kernel'",						///
							"`sw'",							///
							"`psd'",						///
							"`ivar'",						///
							"`tvar'",						///
							"`tindex'",						///
							`tdelta',						///
							`center',						///
							`dofminus',						///
							`sdofminus')

tempname firstmat firstb firstv firsts
			mat `firstmat' = r(firstmat)
			mat rowname `firstmat' =	rmse sheapr2 pr2 F df df_r pvalue			///
										SWF SWFdf1 SWFdf2 SWFp SWchi2 SWchi2p SWr2	///
										APF APFdf1 APFdf2 APFp APchi2 APchi2p APr2
			mat colname `firstmat' = `endo1'
			mat `firstb'	= r(b)
			mat `firstv'	= r(V)
			mat `firsts'	= r(S)
			local archi2	=r(archi2)
			local archi2p	=r(archi2p)
			local arf		=r(arf)
			local arfp		=r(arfp)
			local ardf		=r(ardf)
			local ardf_r	=r(ardf_r)
			local sstat		=r(sstat)
			local sstatdf	=r(sstatdf)
			local sstatp	=r(sstatp)
			local rmse_rf	=r(rmse_rf)
			
* Restore original order if changed for mata code above
			capture tsset
// System of first-stage/reduced form eqns
			if "`sfirst'`savesfirst'" ~= "" {
				PostFirstRF if `touse',					///
							bmat(`firstb')				///
							vmat(`firstv')				///
							smat(`firsts')				///
							firstmat(`firstmat')		///
							lhs1(`lhs1')				///
							endo1(`endo1')				///
							znames0(`cnZ0')				///
							znames1(`cnZ1')				///
							bvclean(`bvclean')			///
							fvops(`fvops')				///
							partial_ct(`partial_ct')	///
							`robust'					///
							cluster(`cluster')			///
							cluster1(`cluster1')		///
							cluster2(`cluster2')		///
							nc(`N_clust')				///
							nc1(`N_clust1')				///
							nc2(`N_clust2')				///
							kernel(`kernel')			///
							bw(`bw')					///
							ivar(`ivar')				///
							tvar(`tvar')				///
							obs(`N')					///
							iv1_ct(`iv1_ct')			///
							cons(`cons')				///
							partialcons(`partialcons')	///
							dofminus(`dofminus')		///
							sdofminus(`sdofminus')
				local sfirsteq "`savesfprefix'sfirst_`lhs1'"
				local sfirsteq : subinstr local sfirsteq "." "_"
				capture est store `sfirsteq', title("System of first-stage/reduced form regressions")
				if _rc > 0 {
di
di in ye "Unable to store system of first-stage reduced form regressions."
di
				}
			}

// RF regression
			if "`rf'`saverf'" ~= "" {
				PostFirstRF if `touse',					///
							rf							/// extract RF regression as saved result
							rmse_rf(`rmse_rf')			/// provide RMSE for posting
							bmat(`firstb')				///
							vmat(`firstv')				///
							smat(`firsts')				///
							firstmat(`firstmat')		///
							lhs1(`lhs1')				///
							endo1(`endo1')				///
							znames0(`cnZ0')				///
							znames1(`cnZ1')				///
							bvclean(`bvclean')			///
							fvops(`fvops')				///
							partial_ct(`partial_ct')	///
							`robust'					///
							cluster(`cluster')			///
							cluster1(`cluster1')		///
							cluster2(`cluster2')		///
							nc(`N_clust')				///
							nc1(`N_clust1')				///
							nc2(`N_clust2')				///
							kernel(`kernel')			///
							bw(`bw')					///
							ivar(`ivar')				///
							tvar(`tvar')				///
							obs(`N')					///
							iv1_ct(`iv1_ct')			///
							cons(`cons')				///
							partialcons(`partialcons')	///
							dofminus(`dofminus')		///
							sdofminus(`sdofminus')
				local rfeq		"`saverfprefix'`lhs1'"
				local rfeq		: subinstr local rfeq "." "_"
				capture est store `rfeq', title("Reduced-form regression: `lhs'")
				if _rc > 0 {
di
di in ye "Unable to store reduced form regression of `lhs1'."
di
				}
			}

// Individual first-stage equations
			if "`first'`savefirst'" ~= "" {
				foreach vn in `endo1' {
	
					PostFirstRF if `touse',					///
								first(`vn')					/// extract first-stage regression
								bmat(`firstb')				///
								vmat(`firstv')				///
								smat(`firsts')				///
								firstmat(`firstmat')		///
								lhs1(`lhs1')				///
								endo1(`endo1')				///
								znames0(`cnZ0')				///
								znames1(`cnZ1')				///
								bvclean(`bvclean')			///
								fvops(`fvops')				///
								partial_ct(`partial_ct')	///
								`robust'					///
								cluster(`cluster')			///
								cluster1(`cluster1')		///
								cluster2(`cluster2')		///
								nc(`N_clust')				///
								nc1(`N_clust1')				///
								nc2(`N_clust2')				///
								kernel(`kernel')			///
								bw(`bw')					///
								ivar(`ivar')				///
								tvar(`tvar')				///
								obs(`N')					///
								iv1_ct(`iv1_ct')			///
								cons(`cons')				///
								partialcons(`partialcons')	///
								dofminus(`dofminus')		///
								sdofminus(`sdofminus')
					local eqname "`savefprefix'`vn'"
					local eqname : subinstr local eqname "." "_"
					capture est store `eqname', title("First-stage regression: `vn'")
					if _rc == 0 {
						local firsteqs "`firsteqs' `eqname'"
					}
					else {
di
di in ye "Unable to store first-stage regression of `vn'."
di
					}
				}
			}
		}
* End of RF and first-stage regression code

*******************************************************************************************
* Re-tsset if necessary
************************************************************************************************

capture tsset

*******************************************************************************************
* orthog option: C statistic (difference of Sargan statistics)
*******************************************************************************************
* Requires j dof from above
		if "`orthog'"!="" {
			tempname cj cstat cstatp
* Initialize cstat
			scalar `cstat' = 0
* Remove orthog from inexog and put in endo
* Remove orthog from exexog
			local cexexog1	: list fv_exexog1 - fv_orthog1
			local cinexog1	: list fv_inexog1 - fv_orthog1
			local cendo1	: list fv_inexog1 - cinexog1
			local cendo1	`fv_endo1' `cendo1'
			local clist_ct	: word count `orthog1'

* If robust, HAC/AC or GMM (but not LIML or IV), create optimal weighting matrix to pass to ivreg2
*   by extracting the submatrix from the full S and then inverting.
*   This guarantees the C stat will be non-negative.  See Hayashi (2000), p. 220. 
* Calculate C statistic with recursive call to ivreg2
* Collinearities may cause problems, hence -capture-.
* smatrix works generally, including homoskedastic case with Sargan stat
			capture _estimates hold `ivest', restore
			if _rc==1000 {
di as err "ivreg2 internal error - no room to save temporary estimation results"
di as err "Please drop one or more estimation results using -estimates drop-"
				exit 1000
			}
* clopt is omitted because it requires calculation of numbers of clusters, which is done
* only when S matrix is calculated
* S matrix has final varnames, but need to call ivreg2 with temp vars
* so must rename cols/rows of S
			tempname fv_S
			mat `fv_S'=`S'
			if `cons' {
				mat colnames `fv_S' = `fv_exexog1' `fv_inexog1' _cons
				mat rownames `fv_S' = `fv_exexog1' `fv_inexog1' _cons
			}
			else {
				mat colnames `fv_S' =  `fv_exexog1' `fv_inexog1'
				mat rownames `fv_S' =  `fv_exexog1' `fv_inexog1'
			}
			capture `ivreg2cmd'		`fv_lhs1'				///
									`cinexog1'				///
									(`cendo1'=`cexexog1')	///
									if `touse'				///
									`wtexp',				///
									`noconstant'			///
									`options'				///
									`small'					///
									`robust'				///
									`gmm2s'					///
									`bwopt'					///
									`kernopt'				///
									`dofmopt'				///
									`sw'					///
									`psd'					///
									smatrix("`fv_S'")		///
									noid					///
									nocollin
			local rc = _rc
			if `rc' == 481 {
				scalar `cstat' = 0
				local cstatdf = 0
			}
			else {
				scalar `cj'=e(j)
				local cjdf=e(jdf)
				scalar `cstat' = `j' - `cj'
				local cstatdf  = `jdf' - `cjdf'
			}
			_estimates unhold `ivest'
			scalar `cstatp'= chiprob(`cstatdf',`cstat')
* Collinearities may cause C-stat dof to differ from the number of variables in orthog()
* If so, set cstat=0
			if `cstatdf' != `clist_ct' {
				scalar `cstat' = 0
			}
		}
* End of orthog block

*******************************************************************************************
* Endog option
*******************************************************************************************
* Uses recursive call with orthog
		if "`endogtest'"!="" {
			tempname estat estatp
* Initialize estat
			scalar `estat' = 0
* Remove endogtest vars from endo and put in inexog
			local eendo1	: list fv_endo1 - fv_endogtest1
			local einexog1	`fv_inexog1' `fv_endogtest1'
			local elist_ct	: word count `endogtest1'

* Recursive call to ivreg2 using orthog option to obtain endogeneity test statistic
* Collinearities may cause problems, hence -capture-.
			capture {
				capture _estimates hold `ivest', restore
				if _rc==1000 {
di as err "ivreg2 internal error - no room to save temporary estimation results"
di as err "Please drop one or more estimation results using -estimates drop-"
					exit 1000
				}
				capture `ivreg2cmd'		`fv_lhs1'				///
										`einexog1'				///
										(`eendo1'=`fv_exexog1')	///
										if `touse'				///
										`wtexp', 				///
										`noconstant'			///
										`robust'				///
										`clopt'					///
										`gmm2s'					///
										`liml'					///
										`bwopt'					///
										`kernopt'				///
										`small'					///
										`dofmopt'				///
										`sw'					///
										`psd'					///
										`options'				///
										orthog(`fv_endogtest1')	///
										noid					///
										nocollin
				local rc = _rc
				if `rc' == 481 {
					scalar `estat' = 0
					local estatdf = 0
					}
				else {
					scalar `estat'=e(cstat)
					local  estatdf=e(cstatdf)
					scalar `estatp'=e(cstatp)
				}
				_estimates unhold `ivest'
* Collinearities may cause endog stat dof to differ from the number of variables in endog()
* If so, set estat=0
				if `estatdf' != `elist_ct' {
					scalar `estat' = 0
				}
			}
* End of endogeneity test block
		}

*******************************************************************************************
* Rank identification and redundancy block
*******************************************************************************************
		if `endo1_ct' > 0 & "`noid'"=="" {

// id=underidentification statistic, wid=weak identification statistic
			tempname idrkstat widrkstat iddf idp
			tempname ccf cdf rkf cceval cdeval cd
			tempname idstat widstat

// UNDERIDENTIFICATION
// Anderson canon corr underidentification statistic if homo, rk stat if not
// Need only id stat for testing full rank=(#cols-1)
// ranktest can exit with error if not full rank
// May not exit with error if e.g. ranktest (x y) (x w),
// i.e. collinearity across lists, so need to catch that.
// If no collinearity, can use iv1_ct and rhs1_ct etc.
			cap `ranktestcmd'							///
								(`fv_endo1')			///
								(`fv_exexog1')			///
								`wtexp'					///
								 if `touse',			///
								 partial(`fv_inexog1')	///
								 full					///
								 `noconstant'			///
								 `robust'				///
								 `clopt'				///
								 `bwopt'				///
								 `kernopt'
// Canonical correlations returned in r(ccorr), sorted in descending order.
// If largest = 1, collinearities so enter error block.
			local rkerror		= _rc>0 | r(chi2)==.
			if ~`rkerror' {
				local rkerror	= el(r(ccorr),1,1)==1
			}
			if `rkerror' {
di as err "warning: -ranktest- error in calculating underidentification test statistics;"
di as err "         may be caused by collinearities"
				scalar `idstat'	= .
				local iddf		= .
				scalar `idp'	= .
				scalar `cd'		= .
				scalar `cdf'	= .
			}
			else {
				if "`cluster'"=="" {
					scalar `idstat'=r(chi2)/r(N)*(`N'-`dofminus')
				}
				else {
// No dofminus adjustment needed for cluster-robust
					scalar `idstat'=r(chi2)
				}
				mat `cceval'=r(ccorr)
				mat `cdeval' = J(1,`endo1_ct',.)
				forval i=1/`endo1_ct' {
					mat `cceval'[1,`i'] = (`cceval'[1,`i'])^2
					mat `cdeval'[1,`i'] = `cceval'[1,`i'] / (1 - `cceval'[1,`i'])
				}
				local iddf = `iv1_ct' - (`rhs1_ct'-1)
				scalar `idp' = chiprob(`iddf',`idstat')
// Cragg-Donald F statistic.
// Under homoskedasticity, Wald cd eigenvalue = cc/(1-cc) Anderson canon corr eigenvalue.
				scalar `cd'=`cdeval'[1,`endo1_ct']
				scalar `cdf'=`cd'*(`N'-`sdofminus'-`iv1_ct'-`dofminus')/`exex1_ct'
			}	// end underidentification stat

// WEAK IDENTIFICATION
// Weak id statistic is Cragg-Donald F stat, rk Wald F stat if not
// ranktest exits with error if not full rank so can use iv1_ct and rhs1_ct etc.
			if "`robust'`cluster'`kernel'"=="" {
				scalar `widstat'=`cdf'
			}
			else {
// Need only test of full rank
				cap `ranktestcmd'						///
								(`fv_endo1')			///
								(`fv_exexog1')			///
								`wtexp'					///
								if `touse',				///
								partial(`fv_inexog1')	///
								full					///
								wald					///
								`noconstant'			///
								`robust'				///
								`clopt'					///
								`bwopt'					///
								`kernopt'
// Canonical correlations returned in r(ccorr), sorted in descending order.
// If largest = 1, collinearities so enter error block.
				local rkerror		= _rc>0 | r(chi2)==.
				if ~`rkerror' {
					local rkerror	= el(r(ccorr),1,1)==1
				}
				if `rkerror' {
di as err "warning: -ranktest- error in calculating weak identification test statistics;"
di as err "         may be caused by collinearities"
					scalar `rkf'		= .
					scalar `widstat'	= .
				}
				else {
// sdofminus used here so that F-stat matches test stat from regression with no partial
					if "`cluster'"=="" {
						scalar `rkf'=r(chi2)/r(N)*(`N'-`iv1_ct'-`sdofminus'-`dofminus')/`exex1_ct'
					}
					else {
						scalar `rkf' =	r(chi2)/(`N'-1) *				///
										(`N'-`iv1_ct'-`sdofminus') *	///
										(`N_clust'-1)/`N_clust' /		///
										`exex1_ct'
					}
					scalar `widstat'=`rkf'
				}
			}	// end weak-identification stat
		}	// end under- and weak-identification stats

* LM redundancy test
		if `endo1_ct' > 0 & "`redundant'" ~= "" & "`noid'"=="" {
* Use K-P rk statistics and LM version of test
* Statistic is the rank of the matrix of Z_1B*X_2, where Z_1B are the possibly redundant
* instruments and X_1 are the endogenous regressors; both have X_2 (exogenous regressors)
* and Z_1A (maintained excluded instruments) partialled out.  LM test of rank is
* is numerically equivalent to estimation of set of RF regressions and performing
* standard LM test of possibly redundant instruments.

local rexexog1		: list fv_exexog1 - fv_redundant1
			local redlist_ct	: word count `redundant1'
* LM version requires only -nullrank- rk statistics so would not need -all- option
			tempname rkmatrix
			qui `ranktestcmd'										///
								(`fv_endo1')						///
								(`fv_redundant1')					///
								`wtexp'								///
								if `touse',							///
								partial(`fv_inexog1' `rexexog1')	///
								null								///
								`noconstant'						///
								`robust'							///
								`clopt'								///
								`bwopt'								///
								`kernopt'
			mat `rkmatrix'=r(rkmatrix)
			tempname redstat redp
* dof adjustment needed because it doesn't use the adjusted S
			if "`cluster'"=="" {
				scalar `redstat' = `rkmatrix'[1,1]/r(N)*(`N'-`dofminus')
			}
			else {
* No dofminus adjustment needed for cluster-robust
				scalar `redstat' = `rkmatrix'[1,1]
			}
			local reddf = `endo1_ct'*`redlist_ct'
			scalar `redp' = chiprob(`reddf',`redstat')
		}

* End of identification stats block

*******************************************************************************************
* Error-checking block
*******************************************************************************************

* Check if adequate number of observations
		if `N' <= `iv1_ct' {
di in r "Error: number of observations must be greater than number of instruments"
di in r "       including constant."
			error 2001
		}

* Check if robust VCV matrix is of full rank
		if ("`gmm2s'`robust'`cluster'`kernel'" != "") & (`rankS' < `iv1_ct') {
* Robust covariance matrix not of full rank means either a singleton dummy or too few
*   clusters (in which case the indiv SEs are OK but no F stat or 2-step GMM is possible),
*   or there are too many AC/HAC-lags, or the HAC covariance estimator
*   isn't positive definite (possible with truncated and Tukey-Hanning kernels)
*   or nocollin option has been used.
* Previous versions of ivreg2 exited if 2-step GMM but beta and VCV may be OK.
* Continue but J, F, and C stat (if present) all possibly meaningless.
* Set j = missing so that problem can be reported in output.
			scalar `j' = .
			if "`orthog'"!="" {
				scalar `cstat' = .
			}
			if "`endogtest'"!="" {
				scalar `estat' = .
			}
		}

* End of error-checking block

**********************************************************************************************
* Post and display results.
*******************************************************************************************

// rankV = rhs1_ct except if nocollin
// rankS = iv1_ct except if nocollin
// nocollin means count may exceed rank (because of dropped vars), so rank #s foolproof

// Add back in omitted vars from "0" varlists unless bvclean requested
// or unless there are no omitted regressors that need adding back in.
		if ~`bvclean' & (`rhs0_ct' > `rhs1_ct') {
			AddOmitted, bmat(`b') vmat(`V') cnb0(`cnb0') cnb1(`cnb1')
			mat `b' = r(b)
			mat `V' = r(V)
// build fv info (base, empty, etc.) unless there was partialling out
			if `fvops' & ~`partial_ct' {
				local bfv "buildfvinfo"
			}
		}

*******************************************************************************************

// restore data if preserved for partial option
		if `partial_ct' {
			restore
		}

if "`small'"!="" {
			local NminusK = `N'-`rankxx'-`sdofminus'
			capture ereturn post `b' `V', dep(`depname') obs(`N') esample(`touse') dof(`NminusK') `bfv'
		}
		else {
			capture ereturn post `b' `V', dep(`depname') obs(`N') esample(`touse') `bfv'
		}

local rc = _rc
		if `rc' == 504 {
di in red "Error: estimated variance-covariance matrix has missing values"
			exit 504
		}
		if `rc' == 506 {
di in red "Error: estimated variance-covariance matrix not positive-definite"
			exit 506
		}
		if `rc' > 0 {
di in red "Error: estimation failed - could not post estimation results"
			exit `rc'
		}

local mok	=1													//  default - margins OK
		local mok	= `mok' & ~`partial_ct'								//  but not if partialling out
		local mok	= `mok' & ~(`fvops' & `bvclean')					//  nor if there are FVs and the base vars are not in e(b)
		if `mok' & `endo1_ct' {											//  margins can be used, endog regressors
			ereturn local marginsnotok	"Residuals SCores"				//  same as official -ivregress-
			ereturn local marginsok		"XB default"
		}
		else if `mok' & ~`endo1_ct' {									//  margins can be used, no endog regressors
			ereturn local marginsok		"XB default"					//  same as official -regress'
		}
		else {															//  don't allow margins
			ereturn local marginsnotok	"Residuals SCores XB default"
		}

// Original varlists without removed duplicates, collinears, etc.
// "0" varlists after removing duplicates and reclassifying vars, and including omitteds, FV base vars, etc.
// "1" varlists without omitted, FV base vars, and partialled-out vars
		ereturn local ecollin	`ecollin'
		ereturn local collin	`collin'
		ereturn local dups		`dups'
		ereturn local partial1	`partial1'
		ereturn local partial	`partial'
		ereturn local inexog1	`inexog1'
		ereturn local inexog0	`inexog0'
		ereturn local inexog	`inexog'
		ereturn local exexog1	`exexog1'
		ereturn local exexog0	`exexog0'
		ereturn local exexog	`exexog'
		ereturn local insts1 	`exexog1' `inexog1'
		ereturn local insts0 	`exexog0' `inexog0'
		ereturn local insts 	`exexog' `inexog'
		ereturn local instd1	`endo1'
		ereturn local instd0	`endo0'
		ereturn local instd		`endo'
		ereturn local depvar1	`lhs1'
		ereturn local depvar0	`lhs0'
		ereturn local depvar	`lhs'

ereturn scalar inexog_ct	=`inexog1_ct'
		ereturn scalar exexog_ct	=`exex1_ct'
		ereturn scalar endog_ct		=`endo1_ct'
		ereturn scalar partial_ct	=`partial_ct'

if "`smatrix'" == "" {
			ereturn matrix S `S'
		}
		else {
			ereturn matrix S `S0'					//  it's a copy so original won't be zapped
		}

* No weighting matrix defined for LIML and kclass
		if "`wmatrix'"=="" & "`liml'`kclassopt'"=="" {
			ereturn matrix W `W'
		}
		else if "`liml'`kclassopt'"=="" {
			ereturn matrix W `wmatrix'				//  it's a copy so original won't be zapped
		}

if "`kernel'"!="" {
			ereturn local kernel "`kernel'"
			ereturn scalar bw=`bw'
			ereturn local tvar "`tvar'"
			if "`ivar'" ~= "" {
				ereturn local ivar "`ivar'"
			}
			if "`bwchoice'" ~= "" {
				ereturn local bwchoice "`bwchoice'"
			}
		}

if "`small'"!="" {
			ereturn scalar df_r=`df_r'
			ereturn local small "small"
		}
		if "`nopartialsmall'"=="" {
			ereturn local partialsmall "small"
		}

if "`robust'" != "" {
			local vce "robust"
		}
		if "`cluster1'" != "" {
			if "`cluster2'"=="" {
				local vce "`vce' cluster"
			}
			else {
				local vce "`vce' two-way cluster"
			}
		}
		if "`kernel'" != "" {
			if "`robust'" != "" {
				local vce "`vce' hac"
			}
			else {
				local vce "`vce' ac"
			}
			local vce "`vce' `kernel' bw=`bw'"
		}
		if "`sw'" != "" {
			local vce "`vce' sw"
		}
		if "`psd'" != "" {
			local vce "`vce' `psd'"
		}
		local vce : list clean vce
		local vce = lower("`vce'")
		ereturn local vce `vce'

if "`cluster'"!="" {
			ereturn scalar N_clust=`N_clust'
			ereturn local clustvar `cluster'
		}
		if "`cluster2'"!="" {
			ereturn scalar N_clust1=`N_clust1'
			ereturn scalar N_clust2=`N_clust2'
			ereturn local clustvar1 `cluster1'
			ereturn local clustvar2 `cluster2'
		}

if "`robust'`cluster'" != "" {
			ereturn local vcetype "Robust"
		}

ereturn scalar df_m=`df_m'
		ereturn scalar sdofminus=`sdofminus'
		ereturn scalar dofminus=`dofminus'
		ereturn scalar center=`center'
		ereturn scalar r2=`r2'
		ereturn scalar rmse=`rmse'
		ereturn scalar rss=`rss'
		ereturn scalar mss=`mss'
		ereturn scalar r2_a=`r2_a'
		ereturn scalar F=`F'
		ereturn scalar Fp=`Fp'
		ereturn scalar Fdf1=`Fdf1'
		ereturn scalar Fdf2=`Fdf2'
		ereturn scalar yy=`yy'
		ereturn scalar yyc=`yyc'
		ereturn scalar r2u=`r2u'
		ereturn scalar r2c=`r2c'
		ereturn scalar condzz=`condzz'
		ereturn scalar condxx=`condxx'
		ereturn scalar rankzz=`rankzz'
		ereturn scalar rankxx=`rankxx'
		ereturn scalar rankS=`rankS'
		ereturn scalar rankV=`rankV'
		ereturn scalar ll = -0.5 * (`N'*ln(2*_pi) + `N'*ln(`rss'/`N') + `N')

* Always save J.  Also save as Sargan if homoskedastic; save A-R if LIML.
		ereturn scalar j=`j'
		ereturn scalar jdf=`jdf'
		if `j' != 0 & `j' != . {
			ereturn scalar jp=`jp'
		}
		if ("`robust'`cluster'"=="") {
			ereturn scalar sargan=`j'
			ereturn scalar sargandf=`jdf'
			if `j' != 0  & `j' != . {
				ereturn scalar sarganp=`jp'
			}
		}
		if "`liml'"!="" {
			ereturn scalar arubin=`arubin'
			ereturn scalar arubin_lin=`arubin_lin'
			if `j' != 0  & `j' != . {
				ereturn scalar arubinp=`arubinp'
				ereturn scalar arubin_linp=`arubin_linp'
			}
			ereturn scalar arubindf=`jdf'
		}

if "`orthog'"!="" {
			ereturn scalar cstat=`cstat'
			if `cstat'!=0  & `cstat' != . {
				ereturn scalar cstatp=`cstatp'
				ereturn scalar cstatdf=`cstatdf'
				ereturn local clist `orthog1'
			}
		}

if "`endogtest'"!="" {
			ereturn scalar estat=`estat'
			if `estat'!=0  & `estat' != . {
				ereturn scalar estatp=`estatp'
				ereturn scalar estatdf=`estatdf'
				ereturn local elist `endogtest1'
			}
		}

if `endo1_ct' > 0 & "`noid'"=="" {
			ereturn scalar idstat=`idstat'
			ereturn scalar iddf=`iddf'
			ereturn scalar idp=`idp'
			ereturn scalar cd=`cd'
			ereturn scalar widstat=`widstat'
			ereturn scalar cdf=`cdf'
			capture ereturn matrix ccev=`cceval'
			capture ereturn matrix cdev `cdeval'
			capture ereturn scalar rkf=`rkf'
		}

if "`redundant'"!="" & "`noid'"=="" {
			ereturn scalar redstat=`redstat'
			ereturn scalar redp=`redp'
			ereturn scalar reddf=`reddf'
			ereturn local  redlist `redundant1'
		}

if "`first'`ffirst'`savefirst'`sfirst'`savesfirst'" != "" & `endo1_ct'>0 {
// Capture here because firstmat may be empty if mvs encountered in 1st stage regressions
			capture ereturn matrix first `firstmat'
			ereturn scalar  arf=`arf'
			ereturn scalar  arfp=`arfp'
			ereturn scalar  archi2=`archi2'
			ereturn scalar  archi2p=`archi2p'
			ereturn scalar  ardf=`ardf'
			ereturn scalar  ardf_r=`ardf_r'
			ereturn scalar  sstat=`sstat'
			ereturn scalar  sstatp=`sstatp'
			ereturn scalar  sstatdf=`sstatdf'
		}
// not saved if empty
		ereturn local   firsteqs `firsteqs'
		ereturn local   rfeq `rfeq'
		ereturn local   sfirsteq `sfirsteq'

if "`liml'"!="" {
			ereturn local model "liml"
			ereturn scalar kclass=`kclass'
			ereturn scalar lambda=`lambda'
			if `fuller' > 0 & `fuller' < . {
				ereturn scalar fuller=`fuller'
			}
		}
		else if "`kclassopt'" != "" {
			ereturn local model "kclass"
			ereturn scalar kclass=`kclass'
		}
		else if "`gmm2s'`cue'`b0'`wmatrix'"=="" {
			if "`endo1'" == "" {
				ereturn local model "ols"
			}
			else {
				ereturn local model "iv"
			}
		}
		else if "`cue'`b0'"~="" {
				ereturn local model "cue"
			}
		else if "`gmm2s'"~="" {
			ereturn local model "gmm2s"
		}
		else if "`wmatrix'"~="" {
				ereturn local model "gmmw"
		}
		else {
* Should never enter here
			ereturn local model "unknown"
		}

if "`weight'" != "" { 
			ereturn local wexp "=`exp'"
			ereturn local wtype `weight'
		}
		ereturn local cmd			`ivreg2cmd'
		ereturn local cmdline		`cmdline'
		ereturn local version		`lversion'
		ereturn scalar nocollin		=("`nocollin'"~="")
		ereturn scalar partialcons	=`partialcons'
		ereturn scalar cons			=`cons'

ereturn local predict "`ivreg2cmd'_p"
		
		if "`e(model)'"=="gmm2s" & "`wmatrix'"=="" {
			local title2 "2-Step GMM estimation"
		}
		else if "`e(model)'"=="gmm2s" & "`wmatrix'"~="" {
			local title2 "2-Step GMM estimation with user-supplied first-step weighting matrix"
		}
		else if "`e(model)'"=="gmmw" {
			local title2 "GMM estimation with user-supplied weighting matrix"
		}
		else if "`e(model)'"=="cue" & "`b0'"=="" {
			local title2 "CUE estimation"
		}
		else if "`e(model)'"=="cue" & "`b0'"~="" {
			local title2 "CUE evaluated at user-supplied parameter vector"
		}
		else if "`e(model)'"=="ols" {
			local title2 "OLS estimation"
		}
		else if "`e(model)'"=="iv" {
			local title2 "IV (2SLS) estimation"
		}
		else if "`e(model)'"=="liml" {
			local title2 "LIML estimation"
		}
		else if "`e(model)'"=="kclass" {
			local title2 "k-class estimation"
		}
		else {
* Should never reach here
			local title2 "unknown estimation"
		}
		if "`e(vcetype)'" == "Robust" {
			local hacsubtitle1 "heteroskedasticity"
		}
		if "`e(kernel)'"!="" & "`e(clustvar)'"=="" {
			local hacsubtitle3 "autocorrelation"
		}
		if "`kiefer'"!="" {
			local hacsubtitle3 "within-cluster autocorrelation (Kiefer)"
		}
		if "`e(clustvar)'"!="" {
			if "`e(clustvar2)'"=="" {
				local hacsubtitle3 "clustering on `e(clustvar)'"
			}
			else {
				local hacsubtitle3 "clustering on `e(clustvar1)' and `e(clustvar2)'"
			}
			if "`e(kernel)'" != "" {
				local hacsubtitle4 "and kernel-robust to common correlated disturbances (Driscoll-Kraay)"
			}
		}
		if "`hacsubtitle1'"~="" & "`hacsubtitle3'" ~= "" {
			local hacsubtitle2 " and "
		}
		if "`title'"=="" {
			ereturn local title "`title1'`title2'"
		}
		else {
			ereturn local title "`title'"
		}
		if "`subtitle'"~="" {
			ereturn local subtitle "`subtitle'"
		}
		local hacsubtitle "`hacsubtitle1'`hacsubtitle2'`hacsubtitle3'"
		if "`b0'"~="" {
			ereturn local hacsubtitleB "Estimates based on supplied parameter vector"
		}
		else if "`hacsubtitle'"~="" & "`gmm2s'`cue'"~="" {
			ereturn local hacsubtitleB "Estimates efficient for arbitrary `hacsubtitle'"
		}
		else if "`wmatrix'"~="" {
			ereturn local hacsubtitleB "Efficiency of estimates dependent on weighting matrix"
		}
		else {
			ereturn local hacsubtitleB "Estimates efficient for homoskedasticity only"
		}
		if "`hacsubtitle'"~="" {
			ereturn local hacsubtitleV "Statistics robust to `hacsubtitle'"
		}
		else {
			ereturn local hacsubtitleV "Statistics consistent for homoskedasticity only"
		}
		if "`hacsubtitle4'"~="" {
			ereturn local hacsubtitleV2 "`hacsubtitle4'"
		}
		if "`sw'"~="" {
			ereturn local hacsubtitleV "Stock-Watson heteroskedastic-robust statistics (BETA VERSION)"
		}
	}

*******************************************************************************************
* Display results unless ivreg2 called just to generate stats or nooutput option

if "`nooutput'" == "" {
	
// Display supplementary first-stage/RF results
		if "`savesfirst'`saverf'`savefirst'" != "" {
			DispStored `"`savesfirst'"' `"`saverf'"' `"`savefirst'"'
		}
		if "`rf'" != "" {
			local eqname "`e(rfeq)'"
			tempname ivest
			_estimates hold `ivest', copy
			capture estimates restore `eqname'
			if _rc != 0 {
di
di in ye "Unable to display stored reduced form estimation."
di
			}
			else {
				DispSFirst "rf" `"`plus'"' `"`level'"' `"`nofooter'"' `"`ivreg2name'"' "`dispopt'"
			}
			_estimates unhold `ivest'
		}
		if "`first'" != "" {
			DispFirst `"`ivreg2name'"'
		}
		if "`sfirst'"!="" {
			local eqname "`e(sfirsteq)'"
			tempname ivest
			_estimates hold `ivest', copy
			capture estimates restore `eqname'
			if _rc != 0 {
di
di in ye "Unable to display stored first-stage/reduced form estimations."
di
			}
			else {
				DispSFirst "sfirst" `"`plus'"' `"`level'"' `"`nofooter'"' `"`ivreg2name'"' "`dispopt'"
			}
			_estimates unhold `ivest'
		}
		if "`first'`ffirst'`sfirst'" != "" {
			DispFFirst `"`ivreg2name'"'
		}

// Display main output.  Can be standard ivreg2, or first-stage-type results
		if "`e(model)'"=="first" | "`e(model)'"=="rf" | "`e(model)'"=="sfirst" {
			DispSFirst "`e(model)'" `"`plus'"' `"`level'"' `"`nofooter'"' `"`ivreg2name'"' "`dispopt'"
		}
		else {
			DispMain `"`noheader'"' `"`plus'"' `"`level'"' `"`nofooter'"' `"`ivreg2name'"' "`dispopt'"
		}
	}

// Drop first stage estimations unless explicitly saved or if replay
	if "`savefirst'" == "" {
		local firsteqs "`e(firsteqs)'"
		foreach eqname of local firsteqs {
			capture estimates drop `eqname'
		}
		ereturn local firsteqs
	}
// Drop reduced form estimation unless explicitly saved or if replay
	if "`saverf'" == "" {
		local eqname "`e(rfeq)'"
		capture estimates drop `eqname'
		ereturn local rfeq
	}
// Drop first stage/reduced form estimation unless explicitly saved or if replay
	if "`savesfirst'" == "" {
		local eqname "`e(sfirsteq)'"
		capture estimates drop `eqname'
		ereturn local sfirsteq
	}

end

*******************************************************************************************
* SUBROUTINES
*******************************************************************************************

// ************* Display system of or single first-stage and/or RF estimations ************ //

program define DispSFirst, eclass
	args model plus level nofooter helpfile dispopt
	version 11.2

di
	if "`model'"=="first" {
di in gr "First-stage regression of `e(depvar)':"
	}
	else if "`model'"=="rf" {
		local strlen = length("`e(depvar)'")+25
di in gr "Reduced-form regression: `e(depvar)'"
di in smcl in gr "{hline `strlen'}"
	}
	else if "`model'"=="sfirst" {
di in gr "System of first-stage/reduced-form regressions:"
di in smcl in gr "{hline 47}"
	}

// Display coefficients etc.
// Header info
	if "`e(hacsubtitleV)'" ~= "" {
di in gr _n "`e(hacsubtitleV)'"
	}
	if "`e(hacsubtitleV2)'" ~= "" {
di in gr "`e(hacsubtitleV2)'"
	}
di in gr "Number of obs = " _col(31) in ye %8.0f e(N)
	if "`e(kernel)'"!="" {
di in gr "  kernel=`e(kernel)'; bandwidth=" `e(bw)'
		if "`e(bwchoice)'"!="" {
di in gr "  `e(bwchoice)'"
		}
di in gr "  time variable (t):  " in ye e(tvar)
		if "`e(ivar)'" != "" {
di in gr "  group variable (i): " in ye e(ivar)
		}
	}
	if "`e(clustvar)'"!="" {
		if "`e(clustvar2)'"=="" {
			local N_clust `e(N_clust)'
			local clustvar `e(clustvar)'
		}
		else {
			local N_clust `e(N_clust1)'
			local clustvar `e(clustvar1)'
		}
di in gr "Number of clusters (`clustvar') = " _col(33) in ye %6.0f `N_clust'
	}
	if "`e(clustvar2)'"!="" {
di in gr "Number of clusters (" "`e(clustvar2)'" ") = " _col(33) in ye %6.0f e(N_clust2)
	}

// Unfortunate but necessary hack here: to suppress message about cluster adjustment of
//   standard error, clear e(clustvar) and then reset it after display
	local cluster `e(clustvar)'
	ereturn local clustvar

// Display output
	ereturn display, `plus' level(`level') `dispopt'
	ereturn local clustvar `cluster'

end

// ************* Display main estimation outpout ************** //

program define DispMain, eclass
	args noheader plus level nofooter helpfile dispopt
	version 11.2
* Prepare for problem resulting from rank(S) being insufficient
* Results from insuff number of clusters, too many lags in HAC,
*   to calculate robust S matrix, HAC matrix not PD, singleton dummy,
*   and indicated by missing value for j stat
* Macro `rprob' is either 1 (problem) or 0 (no problem)
	capture local rprob ("`e(j)'"==".")

if "`noheader'"=="" {
		if "`e(title)'" ~= "" {
di in gr _n "`e(title)'"
			local tlen=length("`e(title)'")
di in gr "{hline `tlen'}"
		}
		if "`e(subtitle)'" ~= "" {
di in gr "`e(subtitle)'"
		}
		if "`e(model)'"=="liml" | "`e(model)'"=="kclass" {
di in gr "k               =" %7.5f `e(kclass)'
		}
		if "`e(model)'"=="liml" {
di in gr "lambda          =" %7.5f `e(lambda)'
		}
		if e(fuller) > 0 & e(fuller) < . {
di in gr "Fuller parameter=" %-5.0f `e(fuller)'
		}
		if "`e(hacsubtitleB)'" ~= "" {
di in gr _n "`e(hacsubtitleB)'" _c
		}
		if "`e(hacsubtitleV)'" ~= "" {
di in gr _n "`e(hacsubtitleV)'"
		}
		if "`e(hacsubtitleV2)'" ~= "" {
di in gr "`e(hacsubtitleV2)'"
		}
		if "`e(kernel)'"!="" {
di in gr "  kernel=`e(kernel)'; bandwidth=" `e(bw)'
			if "`e(bwchoice)'"!="" {
di in gr "  `e(bwchoice)'"
			}
di in gr "  time variable (t):  " in ye e(tvar)
			if "`e(ivar)'" != "" {
di in gr "  group variable (i): " in ye e(ivar)
			}
		}
		di
		if "`e(clustvar)'"!="" {
			if "`e(clustvar2)'"=="" {
				local N_clust `e(N_clust)'
				local clustvar `e(clustvar)'
			}
			else {
				local N_clust `e(N_clust1)'
				local clustvar `e(clustvar1)'
			}
di in gr "Number of clusters (`clustvar') = " _col(33) in ye %6.0f `N_clust' _continue
		}
di in gr _col(55) "Number of obs = " in ye %8.0f e(N)
		if "`e(clustvar2)'"!="" {
di in gr "Number of clusters (" "`e(clustvar2)'" ") = " _col(33) in ye %6.0f e(N_clust2) _continue
		}
di in gr _c _col(55) "F(" %3.0f e(Fdf1) "," %6.0f e(Fdf2) ") = "
		if e(F) < 99999 {
di in ye %8.2f e(F)
		}
		else {
di in ye %8.2e e(F)
		}
di in gr _col(55) "Prob > F      = " in ye %8.4f e(Fp)

di in gr "Total (centered) SS     = " in ye %12.0g e(yyc) _continue
di in gr _col(55) "Centered R2   = " in ye %8.4f e(r2c)
di in gr "Total (uncentered) SS   = " in ye %12.0g e(yy) _continue
di in gr _col(55) "Uncentered R2 = " in ye %8.4f e(r2u)
di in gr "Residual SS             = " in ye %12.0g e(rss) _continue
di in gr _col(55) "Root MSE      = " in ye %8.4g e(rmse)
di
	}

* Display coefficients etc.
* Unfortunate but necessary hack here: to suppress message about cluster adjustment of
*   standard error, clear e(clustvar) and then reset it after display
	local cluster `e(clustvar)'
	ereturn local clustvar
	ereturn display, `plus' level(`level') `dispopt'
	ereturn local clustvar `cluster'

* Display 1st footer with identification stats
* Footer not displayed if -nofooter- option or if pure OLS, i.e., model="ols" and Sargan-Hansen=0
	if ~("`nofooter'"~="" | (e(model)=="ols" & (e(sargan)==0 | e(j)==0))) {

* Under ID test
		if "`e(instd)'"~="" & "`e(idstat)'"~="" {
di in smcl _c "{help `helpfile'##idtest:Underidentification test}"
			if "`e(vcetype)'`e(kernel)'"=="" {
di in gr _c " (Anderson canon. corr. LM statistic):"
			}
			else {
di in gr _c " (Kleibergen-Paap rk LM statistic):"
			}
di in ye _col(71) %8.3f e(idstat)
di in gr _col(52) "Chi-sq(" in ye e(iddf) /* 
	       			*/  in gr ") P-val =  " in ye _col(73) %6.4f e(idp)
* IV redundancy statistic
			if "`e(redlist)'"!="" {
di in gr "-redundant- option:"
di in smcl _c "{help `helpfile'##redtest:IV redundancy test}"
di in gr _c " (LM test of redundancy of specified instruments):"
di in ye _col(71) %8.3f e(redstat)
di in gr _col(52) "Chi-sq(" in ye e(reddf) /* 
	       			*/  in gr ") P-val =  " in ye _col(73) %6.4f e(redp)
di in gr "Instruments tested: " _c
					Disp `e(redlist)', _col(23)
			}
di in smcl in gr "{hline 78}"
		}
* Report Cragg-Donald statistic
		if "`e(instd)'"~="" & "`e(idstat)'"~="" {
di in smcl _c "{help `helpfile'##widtest:Weak identification test}"
di in gr " (Cragg-Donald Wald F statistic):" in ye _col(71) %8.3f e(cdf)
			if "`e(vcetype)'`e(kernel)'"~="" {
di in gr "                         (Kleibergen-Paap rk Wald F statistic):" in ye _col(71) %8.3f e(widstat)
			}
di in gr _c "Stock-Yogo weak ID test critical values:"
			Disp_cdsy, model(`e(model)') k2(`e(exexog_ct)') nendog(`e(endog_ct)') fuller("`e(fuller)'") col1(42) col2(73)
			if `r(cdmissing)' {
				di in gr _col(64) "<not available>"
			}
			else {
				di in gr "Source: Stock-Yogo (2005).  Reproduced by permission."
				if "`e(vcetype)'`e(kernel)'"~="" {
di in gr "NB: Critical values are for Cragg-Donald F statistic and i.i.d. errors."
				}
			}
			di in smcl in gr "{hline 78}"
		}

* Report either (a) Sargan-Hansen-C stats, or (b) robust covariance matrix problem
* e(model)="gmmw" means user-supplied weighting matrix and Hansen J using 2nd-step resids reported
		if `rprob' == 0 {
* Display overid statistic
			if "`e(vcetype)'" == "Robust" | "`e(model)'" == "gmmw" {
				if "`e(instd)'" != "" {
di in smcl _c "{help `helpfile'##overidtests:Hansen J statistic}"
di in gr _c " (overidentification test of all instruments):"
				}
				else {
di in smcl _c "{help `helpfile'##overidtests:Hansen J statistic}"
di in gr _c " (Lagrange multiplier test of excluded instruments):"
				}
			}
			else {
				if "`e(instd)'" != "" {
di in smcl _c "{help `helpfile'##overidtests:Sargan statistic}"
di in gr _c " (overidentification test of all instruments):"
				}
				else {
di in smcl _c "{help `helpfile'##overidtests:Sargan statistic}"
di in gr _c " (Lagrange multiplier test of excluded instruments):"
				}
			}
di in ye _col(71) %8.3f e(j)
			if e(jdf) {
di in gr _col(52) "Chi-sq(" in ye e(jdf) /* 
	       			*/  in gr ") P-val =  " in ye _col(73) %6.4f e(jp)
			}
			else {
di in gr _col(50) "(equation exactly identified)"
			}

* Display orthog option: C statistic (difference of Sargan statistics)
			if e(cstat) != . {
* If C-stat = 0 then warn, otherwise output
				if e(cstat) > 0  {
di in gr "-orthog- option:"
					if "`e(vcetype)'" == "Robust" {
di in gr _c "Hansen J statistic (eqn. excluding suspect orthog. conditions): "
					}
					else {
di in gr _c "Sargan statistic (eqn. excluding suspect orthogonality conditions):"
					}
di in ye _col(71) %8.3f e(j)-e(cstat)
di in gr _col(52) "Chi-sq(" in ye e(jdf)-e(cstatdf) in gr ") P-val =  " /*
				*/ in ye _col(73) %6.4f chiprob(e(jdf)-e(cstatdf),e(j)-e(cstat))
di in smcl _c "{help `helpfile'##ctest:C statistic}"
di in gr _c " (exogeneity/orthogonality of suspect instruments): "
di in ye _col(71) %8.3f e(cstat)
di in gr _col(52) "Chi-sq(" in ye e(cstatdf) in gr ") P-val =  " /*
				*/ in ye _col(73) %6.4f e(cstatp)
di in gr "Instruments tested:  " _c
					Disp `e(clist)', _col(23)
				}
				if e(cstat) == 0 {
di in gr _n "Collinearity/identification problems in eqn. excl. suspect orthog. conditions:"
di in gr "  C statistic not calculated for -orthog- option"
				}
			}
		}
		else {
* Problem exists with robust VCV - notify and list possible causes
di in r "Warning: estimated covariance matrix of moment conditions not of full rank."
			if e(j)==. {
di in r "         overidentification statistic not reported, and standard errors and"
			}
di in r "         model tests should be interpreted with caution."
di in r "Possible causes:"
			if e(nocollin) {
di in r "         collinearities in regressors or instruments (with -nocollin- option)"
			}
			if "`e(N_clust)'" != "" {
di in r "         number of clusters insufficient to calculate robust covariance matrix"
			}
			if "`e(kernel)'" != "" {
di in r "         covariance matrix of moment conditions not positive definite"
di in r "         covariance matrix uses too many lags"
			}
di in r "         singleton dummy variable (dummy with one 1 and N-1 0s or vice versa)"
di in r in smcl _c "{help `helpfile'##partial:partial}"
di in r " option may address problem."
		}

* Display endog option: endogeneity test statistic
		if e(estat) != . {
* If stat = 0 then warn, otherwise output
			if e(estat) > 0  {
di in gr "-endog- option:"
di in smcl _c "{help `helpfile'##endogtest:Endogeneity test}"
di in gr _c " of endogenous regressors: "
di in ye _col(71) %8.3f e(estat)
di in gr _col(52) "Chi-sq(" in ye e(estatdf) /* 
       			*/  in gr ") P-val =  " in ye _col(73) %6.4f e(estatp)
di in gr "Regressors tested:  " _c
				Disp `e(elist)', _col(23)
			}
			if e(estat) == 0 {
di in gr _n "Collinearity/identification problems in restricted equation:"
di in gr "  Endogeneity test statistic not calculated for -endog- option"
			}
		}

di in smcl in gr "{hline 78}"
* Display AR overid statistic if LIML and not robust
		if "`e(model)'" == "liml" & "`e(vcetype)'" ~= "Robust" & "`e(kernel)'" == "" {
			if "`e(instd)'" != "" {
di in smcl _c "{help `helpfile'##liml:Anderson-Rubin statistic}"
di in gr _c " (overidentification test of all instruments):"
				}
				else {
di in smcl _c "{help `helpfile'##liml:Anderson-Rubin statistic}"
di in gr _c " (LR test of excluded instruments):"
				}
di in ye _col(72) %7.3f e(arubin)
			if e(arubindf) {
di in gr _col(52) "Chi-sq(" in ye e(arubindf) /* 
	       			*/  in gr ") P-val =  " in ye _col(73) %6.4f e(arubinp)
			}
			else {
di in gr _col(50) "(equation exactly identified)"
			}
			di in smcl in gr "{hline 78}"
		}
	}

* Display 2nd footer with variable lists
	if "`nofooter'"=="" {

* Warn about dropped instruments if any
* Can happen with nocollin option and rank(S) < cols(S)
		if colsof(e(S)) > e(rankzz) {
di in gr "Collinearities detected among instruments: " _c
di in gr colsof(e(S))-e(rankzz) " instrument(s) dropped"
		}

if "`e(collin)'`e(dups)'" != "" | e(partial_ct) {
* If collinearities, duplicates or partial, abbreviated varlists saved with a 1 at the end
			local one "1"
		}
		if e(endog_ct) {
			di in gr "Instrumented:" _c
			Disp `e(instd`one')', _col(23)
		}
		if e(inexog_ct) {
			di in gr "Included instruments:" _c
			Disp `e(inexog`one')', _col(23)
		}
		if e(exexog_ct) {
			di in gr "Excluded instruments:" _c
			Disp `e(exexog`one')', _col(23)
		}
		if e(partial_ct) {
			if e(partialcons) {
				local partial "`e(partial`one')' _cons"
			}
			else {
				local partial "`e(partial`one')'"
			}
di in smcl _c "{help `helpfile'##partial:Partialled-out}"
			di in gr ":" _c
			Disp `partial', _col(23)
			if "`e(partialsmall)'"=="" {
di in gr _col(23) "nb: total SS, model F and R2s are after partialling-out;"
di in gr _col(23) "    any {help `helpfile'##s_small:small-sample adjustments} do not include"
di in gr _col(23) "    partialled-out variables in regressor count K"
			}
			else {
di in gr _col(23) "nb: total SS, model F and R2s are after partialling-out;"
di in gr _col(23) "    any {help `helpfile'##s_small:small-sample adjustments} include partialled-out"
di in gr _col(23) "    variables in regressor count K"
			}
		}
		if "`e(dups)'" != "" {
			di in gr "Duplicates:" _c
			Disp `e(dups)', _col(23)
		}
		if "`e(collin)'" != "" {
			di in gr "Dropped collinear:" _c
			Disp `e(collin)', _col(23)
		}
		if "`e(ecollin)'" != "" {
			di in gr "Reclassified as exog:" _c
			Disp `e(ecollin)', _col(23)
		}
		di in smcl in gr "{hline 78}"
	}
end

**************************************************************************************

// ************ Display collinearity and duplicates warning messages ************ //

program define DispCollinDups
	version 11.2
	if "`e(dups)'" != "" {
di in gr "Warning - duplicate variables detected"
di in gr "Duplicates:" _c
		Disp `e(dups)', _col(16)
	}
	if "`e(collin)'" != "" {
di in gr "Warning - collinearities detected"
di in gr "Vars dropped:" _c
		Disp `e(collin)', _col(16)
	}
end

// ************* Display all first-stage estimations ************ //

program define DispFirst
	version 11.2
	args helpfile
	tempname firstmat ivest sheapr2 pr2 F df df_r pvalue
	tempname SWF SWFdf1 SWFdf2 SWFp SWr2

mat `firstmat'=e(first)
	if `firstmat'[1,1] == . {
di
di in ye "Unable to display first-stage estimates; macro e(first) is missing"
		exit
	}
di in gr _newline "First-stage regressions"
di in smcl in gr "{hline 23}"
di
	local endo1 : colnames(`firstmat')
	local nrvars : word count `endo1'
	local firsteqs "`e(firsteqs)'"
	local nreqs : word count `firsteqs'
	if `nreqs' < `nrvars' {
di in ye "Unable to display all first-stage regressions."
di in ye "There may be insufficient room to store results using -estimates store-,"
di in ye "or names of endogenous regressors may be too long to store the results."
di in ye "Try dropping one or more estimation results using -estimates drop-,"
di in ye "using the -savefprefix- option, or using shorter variable names."
di
	}
	local robust "`e(vcetype)'"
	local cluster "`e(clustvar)'"
	local kernel "`e(kernel)'"
	foreach eqname of local firsteqs {
		_estimates hold `ivest'
		capture estimates restore `eqname'
		if _rc != 0 {
di
di in ye "Unable to list stored estimation `eqname'."
di in ye "There may be insufficient room to store results using -estimates store-,"
di in ye "or names of endogenous regressors may be too long to store the results."
di in ye "Try dropping one or more estimation results using -estimates drop-,"
di in ye "using the -savefprefix- option, or using shorter variable names."
di
		}
		else {
			local vn "`e(depvar)'"
			estimates replay `eqname', noheader
			mat `sheapr2' =`firstmat'["sheapr2","`vn'"]
			mat `pr2'     =`firstmat'["pr2","`vn'"]
			mat `F'       =`firstmat'["F","`vn'"]
			mat `df'      =`firstmat'["df","`vn'"]
			mat `df_r'    =`firstmat'["df_r","`vn'"]
			mat `pvalue'  =`firstmat'["pvalue","`vn'"]
			mat `SWF'     =`firstmat'["SWF","`vn'"]
			mat `SWFdf1'  =`firstmat'["SWFdf1","`vn'"]
			mat `SWFdf2'  =`firstmat'["SWFdf2","`vn'"]
			mat `SWFp'    =`firstmat'["SWFp","`vn'"]
			mat `SWr2'    =`firstmat'["SWr2","`vn'"]

di in gr "F test of excluded instruments:"
di in gr "  F(" %3.0f `df'[1,1] "," %6.0f `df_r'[1,1] ") = " in ye %8.2f `F'[1,1]
di in gr "  Prob > F      = " in ye %8.4f `pvalue'[1,1]

di in smcl "{help `helpfile'##swstats:Sanderson-Windmeijer multivariate F test of excluded instruments:}"
di in gr "  F(" %3.0f `SWFdf1'[1,1] "," %6.0f `SWFdf2'[1,1] ") = " in ye %8.2f `SWF'[1,1]
di in gr "  Prob > F      = " in ye %8.4f `SWFp'[1,1]

di
		}
		_estimates unhold `ivest'
	}
end

// ************* Display list of stored first-stage and RF estimations ************ //

program define DispStored
	args savesfirst saverf savefirst
	version 11.2

if "`savesfirst'" != "" {
		local eqlist "`e(sfirsteq)'"
	}
	if "`saverf'" != "" {
		local eqlist "`eqlist' `e(rfeq)'"
	}
	if "`savefirst'" != "" {
		local eqlist "`eqlist' `e(firsteqs)'"
	}
	local eqlist : list retokenize eqlist

di in gr _newline "Stored estimation results"
di in smcl in gr "{hline 25}" _c
	capture estimates dir `eqlist'
	if "`eqlist'" != "" & _rc == 0 {
// Estimates exist and can be listed
		estimates dir `eqlist'
	}
	else if "`eqlist'" != "" & _rc != 0 {
di
di in ye "Unable to list stored estimations."
di
	}
end

// ************* Display summary first-stage and ID test results ************ //

program define DispFFirst
	version 11.2
	args helpfile
	tempname firstmat
	tempname sheapr2 pr2 F df df_r pvalue
	tempname SWF SWFdf1 SWFdf2 SWFp SWchi2 SWchi2p SWr2
	mat `firstmat'=e(first)
	if `firstmat'[1,1] == . {
di
di in ye "Unable to display summary of first-stage estimates; macro e(first) is missing"
		exit
	}
	local endo   : colnames(`firstmat')
	local nrvars : word count `endo'
	local robust   "`e(vcetype)'"
	local cluster  "`e(clustvar)'"
	local kernel   "`e(kernel)'"
	local efirsteqs  "`e(firsteqs)'"

mat `df'      =`firstmat'["df",1]
	mat `df_r'    =`firstmat'["df_r",1]
	mat `SWFdf1'  =`firstmat'["SWFdf1",1]
	mat `SWFdf2'  =`firstmat'["SWFdf2",1]

di
di in gr _newline "Summary results for first-stage regressions"
di in smcl in gr "{hline 43}"
di

di _c in smcl _col(44) "{help `helpfile'##swstats:(Underid)}"
di    in smcl _col(65) "{help `helpfile'##swstats:(Weak id)}"

di _c in gr "Variable     |"
di _c in smcl _col(16) "{help `helpfile'##swstats:F}" in gr "("
di _c in ye _col(17) %3.0f `df'[1,1] in gr "," in ye %6.0f `df_r'[1,1] in gr ")  P-val"
di _c in gr _col(37) "|"
di _c in smcl _col(39) "{help `helpfile'##swstats:SW Chi-sq}" in gr "("
di _c in ye %3.0f `SWFdf1'[1,1] in gr ") P-val"
di _c in gr _col(60) "|"
di _c in smcl _col(62) "{help `helpfile'##swstats:SW F}" in gr "("
di    in ye _col(67) %3.0f `SWFdf1'[1,1] in gr "," in ye %6.0f `SWFdf2'[1,1] in gr ")"

local i = 1
	foreach vn of local endo {
	
		mat `sheapr2' =`firstmat'["sheapr2","`vn'"]
		mat `pr2'     =`firstmat'["pr2","`vn'"]
		mat `F'       =`firstmat'["F","`vn'"]
		mat `df'      =`firstmat'["df","`vn'"]
		mat `df_r'    =`firstmat'["df_r","`vn'"]
		mat `pvalue'  =`firstmat'["pvalue","`vn'"]
		mat `SWF'     =`firstmat'["SWF","`vn'"]
		mat `SWFdf1'  =`firstmat'["SWFdf1","`vn'"]
		mat `SWFdf2'  =`firstmat'["SWFdf2","`vn'"]
		mat `SWFp'    =`firstmat'["SWFp","`vn'"]
		mat `SWchi2'  =`firstmat'["SWchi2","`vn'"]
		mat `SWchi2p' =`firstmat'["SWchi2p","`vn'"]
		mat `SWr2'    =`firstmat'["SWr2","`vn'"]

local vnlen : length local vn
		if `vnlen' > 12 {
			local vn : piece 1 12 of "`vn'"
		}
di _c in y %-12s "`vn'" _col(14) in gr "|" _col(18) in y %8.2f `F'[1,1]
di _c _col(28) in y %8.4f  `pvalue'[1,1]
di _c _col(37) in g "|" _col(42) in y %8.2f `SWchi2'[1,1] _col(51) in y %8.4f  `SWchi2p'[1,1]
di    _col(60) in g "|" _col(65) in y %8.2f `SWF'[1,1] 
		local i = `i' + 1
	}
di

if "`robust'`cluster'" != "" {
		if "`cluster'" != "" {
			local rtype "cluster-robust"
		}
		else if "`kernel'" != "" {
			local rtype "heteroskedasticity and autocorrelation-robust"
		}
		else {
			local rtype "heteroskedasticity-robust"
		}
	}
	else if "`kernel'" != "" {
			local rtype "autocorrelation-robust"
	}
	if "`robust'`cluster'`kernel'" != "" {
di in gr "NB: first-stage test statistics `rtype'"
di
	}

local k2 = `SWFdf1'[1,1]
di in gr "Stock-Yogo weak ID F test critical values for single endogenous regressor:"
	Disp_cdsy, model(`e(model)') k2(`e(exexog_ct)') nendog(1) fuller("`e(fuller)'") col1(36) col2(67)
	if `r(cdmissing)' {
		di in gr _col(64) "<not available>"
	}
	else {
		di in gr "Source: Stock-Yogo (2005).  Reproduced by permission."
		if "`e(model)'"=="iv" & "`e(vcetype)'`e(kernel)'"=="" {
di in gr "NB: Critical values are for Sanderson-Windmeijer F statistic."
		}
		else {
di in gr "NB: Critical values are for i.i.d. errors only."
		}
di
	}

* Check that SW chi-sq and F denominator are correct and = underid test dof
	if e(iddf)~=`SWFdf1'[1,1] {
di in red "Warning: Error in calculating first-stage id statistics above;"
di in red "         dof of SW statistics is " `SWFdf1'[1,1] ", should be L-(K-1)=`e(iddf)'."
	}

tempname iddf idstat idp widstat cdf rkf
	scalar `iddf'=e(iddf)
	scalar `idstat'=e(idstat)
	scalar `idp'=e(idp)
	scalar `widstat'=e(widstat)
	scalar `cdf'=e(cdf)
	capture scalar `rkf'=e(rkf)
di in smcl "{help `helpfile'##idtest:Underidentification test}"
di in gr "Ho: matrix of reduced form coefficients has rank=K1-1 (underidentified)"
di in gr "Ha: matrix has rank=K1 (identified)"
	if "`robust'`kernel'"=="" {
di in ye "Anderson canon. corr. LM statistic" _c
	}
	else {
di in ye "Kleibergen-Paap rk LM statistic" _c
	}
di in gr _col(42) "Chi-sq(" in ye `iddf' in gr ")=" %-7.2f in ye `idstat' /*
	*/ _col(61) in gr "P-val=" %6.4f in ye `idp'

di
di in smcl "{help `helpfile'##widtest:Weak identification test}"
di in gr "Ho: equation is weakly identified"
di in ye "Cragg-Donald Wald F statistic" _col(65) %8.2f `cdf'
	if "`robust'`kernel'"~="" {
di in ye "Kleibergen-Paap Wald rk F statistic"  _col(65) %8.2f `rkf'
	}
di

di in gr "Stock-Yogo weak ID test critical values for K1=`e(endog_ct)' and L1=`e(exexog_ct)':"
	Disp_cdsy, model(`e(model)') k2(`e(exexog_ct)') nendog(`e(endog_ct)') fuller("`e(fuller)'") col1(36) col2(67)
	if `r(cdmissing)' {
		di in gr _col(64) "<not available>"
	}
	else {
		di in gr "Source: Stock-Yogo (2005).  Reproduced by permission."
		if "`e(vcetype)'`e(kernel)'"~="" {
di in gr "NB: Critical values are for Cragg-Donald F statistic and i.i.d. errors."
		}
	}
di

tempname arf arfp archi2 archi2p ardf ardf_r
	tempname sstat sstatp sstatdf
di in smcl "{help `helpfile'##wirobust:Weak-instrument-robust inference}"
di in gr "Tests of joint significance of endogenous regressors B1 in main equation"
di in gr "Ho: B1=0 and orthogonality conditions are valid"
* Needs to be small so that adjusted dof is reflected in F stat
	scalar `arf'=e(arf)
	scalar `arfp'=e(arfp)
	scalar `archi2'=e(archi2)
	scalar `archi2p'=e(archi2p)
	scalar `ardf'=e(ardf)
	scalar `ardf_r'=e(ardf_r)
	scalar `sstat'=e(sstat)
	scalar `sstatp'=e(sstatp)
	scalar `sstatdf'=e(sstatdf)
di in ye _c "Anderson-Rubin Wald test"
di in gr _col(36) "F(" in ye `ardf' in gr "," in ye `ardf_r' in gr ")=" /*
		*/	_col(49) in ye %7.2f `arf'    _col(61) in gr "P-val=" in ye %6.4f `arfp'
di in ye _c "Anderson-Rubin Wald test"
di in gr _col(36) "Chi-sq(" in ye `ardf' in gr ")=" /*
		*/	_col(49) in ye %7.2f `archi2' _col(61) in gr "P-val=" in ye %6.4f `archi2p'
di in ye _c "Stock-Wright LM S statistic"
di in gr _col(36) "Chi-sq(" in ye `sstatdf' in gr ")=" /*
		*/	_col(49) in ye %7.2f `sstat' _col(61) in gr "P-val=" in ye %6.4f `sstatp'
di
	if "`robust'`cluster'`kernel'" != "" {
di in gr "NB: Underidentification, weak identification and weak-identification-robust"
di in gr "    test statistics `rtype'"
di
	}

if "`cluster'" != "" & "`e(clustvar2)'"=="" {
di in gr "Number of clusters             N_clust  = " in ye %10.0f e(N_clust)
	}
	else if "`e(clustvar2)'" ~= "" {
di in gr "Number of clusters (1)         N_clust1 = " in ye %10.0f e(N_clust1)
di in gr "Number of clusters (2)         N_clust2 = " in ye %10.0f e(N_clust2)
	}
di in gr "Number of observations               N  = " in ye %10.0f e(N)
di in gr "Number of regressors                 K  = " in ye %10.0f e(rankxx)
di in gr "Number of endogenous regressors      K1 = " in ye %10.0f e(endog_ct)
di in gr "Number of instruments                L  = " in ye %10.0f e(rankzz)
di in gr "Number of excluded instruments       L1 = " in ye %10.0f e(ardf)
	if "`e(partial)'" != "" {
di in gr "Number of partialled-out regressors/IVs = " in ye %10.0f e(partial_ct)
di in gr "NB: K & L do not included partialled-out variables"
	}

end

// ************* Post first-stage and/or RF estimations ************ //

program define PostFirstRF, eclass
	version 11.2
	syntax [if]							///
		 [ ,							///
		 		first(string)			/// can be fv
		 		rf						/// omit first(.) and rf => post system of eqns
		 		rmse_rf(real 0)			///
				bmat(name)				///
				vmat(name)				///
				smat(name)				///
				firstmat(name)			///
				lhs1(string)			/// can be fv
				endo1(string)			///
				znames0(string)			///
				znames1(string)			///
				bvclean(integer 0)		///
				fvops(integer 0)		///
				partial_ct(integer 0)	///
				robust					///
				cluster(string)			///
				cluster1(string)		///
				cluster2(string)		///
				nc(integer 0)			///
				nc1(integer 0)			///
				nc2(integer 0)			///
				kernel(string)			///
				bw(real 0)				///
				ivar(name)				///
				tvar(name)				///
				obs(integer 0)			///
				iv1_ct(integer 0)		///
				cons(integer 0)			///
				partialcons(integer 0)	///
				dofminus(integer 0)		///
				sdofminus(integer 0)	///
			]

// renaming/copying
	local N			= `obs'
	local N_clust	= `nc'
	local N_clust1	= `nc1'
	local N_clust2	= `nc2'
	tempname b V S
	mat `b' = `bmat'
	mat `V' = `vmat'
	mat `S' = `smat'

marksample touse

mat colname `b' = `lhs1' `endo1'
	mat rowname `b' = `znames1'
	mat `b' = vec(`b')
	mat `b' = `b''
	mat colname `V' = `: colfullnames `b''
	mat rowname `V' = `: colfullnames `b''
	mat colname `S' = `: colfullnames `b''
	mat rowname `S' = `: colfullnames `b''
	
	if "`cluster'"=="" {
		matrix `V'=`V'*(`N'-`dofminus')/(`N'-`iv1_ct'-`dofminus'-`sdofminus')
	}
	else {
		matrix `V'=`V'*(`N'-1)/(`N'-`iv1_ct'-`sdofminus')		///
					* `N_clust'/(`N_clust'-1)
	}

// If RF or first-stage estimation required, extract it
// also set macros for model and depvar
	if "`rf'`first'"~="" {
		if "`rf'"~="" {											// RF
			local vnum		= 0
			local model		rf
			local depvar	`lhs1'
			local rmse		= `rmse_rf'
		}
		else {													// first-stage
			local vnum		: list posof "`first'" in endo1
			local vnum		= `vnum'
			local model 	first
			local depvar	`first'
			local rmse		= el(`firstmat', rownumb(`firstmat',"rmse"), colnumb(`firstmat',"`first'"))
		}
		local c0		= 1 + `vnum'*`iv1_ct'
		local c1		= (`vnum'+1)*`iv1_ct'
		mat `b'			= `b'[1,`c0'..`c1']
		mat `V'			= `V'[`c0'..`c1',`c0'..`c1']
		mat `S'			= `S'[`c0'..`c1',`c0'..`c1']
		mat coleq `b'	= ""
		mat coleq `V'	= ""
		mat roweq `V'	= ""
		mat coleq `S'	= ""
		mat roweq `S'	= ""
	}
	else {
		local model		sfirst
		local eqlist	`lhs1' `endo1'
	}

// reinsert omitteds etc. unless requested not to
// eqlist empty unless first-stage/rf system
	if ~`bvclean' {
		AddOmitted, bmat(`b') vmat(`V') cnb0(`znames0') cnb1(`znames1') eqlist(`eqlist')
		mat `b' = r(b)
		mat `V' = r(V)
// build fv info (base, empty, etc.) unless there was partialling out
		if `fvops' & ~`partial_ct' {
			local bfv "buildfvinfo"
		}
	}

local dof	= `N' - `iv1_ct' - `dofminus' - `sdofminus'
	ereturn post `b' `V', obs(`obs') esample(`touse') dof(`dof') depname(`depvar') `bfv'

// saved RF/first-stage equation scalars
	if "`rf'`first'"~="" {
		ereturn scalar rmse = `rmse'
		ereturn scalar df_r = `dof'
		ereturn scalar df_m = `iv1_ct' - `cons' + `sdofminus' - `partialcons'
	}
	ereturn scalar k_eq	= `: word count `endo1''
	ereturn local cmd	ivreg2
	ereturn local model	`model'
	ereturn matrix S	`S'

if "`kernel'"!="" {
		ereturn local kernel "`kernel'"
		ereturn scalar bw=`bw'
		ereturn local tvar "`tvar'"
		if "`ivar'" ~= "" {
			ereturn local ivar "`ivar'"
		}
	}

local vce : list clean vce
	local vce = lower("`vce'")
	ereturn local vce `vce'

if "`cluster'"!="" {
		ereturn scalar N_clust=`N_clust'
		ereturn local clustvar `cluster'
	}
	if "`cluster2'"!="" {
		ereturn scalar N_clust1=`N_clust1'
		ereturn scalar N_clust2=`N_clust2'
		ereturn local clustvar1 `cluster1'
		ereturn local clustvar2 `cluster2'
	}

if "`robust'`cluster'" != "" {
		ereturn local vcetype "Robust"
	}

// Assemble output titles
	if "`e(vcetype)'" == "Robust" {
		local hacsubtitle1 "heteroskedasticity"
	}
	if "`e(kernel)'"!="" & "`e(clustvar)'"=="" {
		local hacsubtitle3 "autocorrelation"
	}
	if "`kiefer'"!="" {
		local hacsubtitle3 "within-cluster autocorrelation (Kiefer)"
	}
	if "`e(clustvar)'"!="" {
		if "`e(clustvar2)'"=="" {
			local hacsubtitle3 "clustering on `e(clustvar)'"
		}
		else {
			local hacsubtitle3 "clustering on `e(clustvar1)' and `e(clustvar2)'"
		}
		if "`e(kernel)'" != "" {
			local hacsubtitle4 "and kernel-robust to common correlated disturbances (Driscoll-Kraay)"
		}
	}
	if "`hacsubtitle1'"~="" & "`hacsubtitle3'" ~= "" {
		local hacsubtitle2 " and "
	}
	local hacsubtitle "`hacsubtitle1'`hacsubtitle2'`hacsubtitle3'"
	if "`hacsubtitle'"~="" {
		ereturn local hacsubtitleV "Statistics robust to `hacsubtitle'"
	}
	else {
		ereturn local hacsubtitleV "Statistics consistent for homoskedasticity only"
	}
	if "`hacsubtitle4'"~="" {
		ereturn local hacsubtitleV2 "`hacsubtitle4'"
	}
	if "`sw'"~="" {
		ereturn local hacsubtitleV "Stock-Watson heteroskedastic-robust statistics (BETA VERSION)"
	}

end

**************************************************************************************
program define IsStop, sclass
				/* sic, must do tests one-at-a-time, 
				 * 0, may be very large */
	version 11.2
	if `"`0'"' == "[" {		
		sret local stop 1
		exit
	}
	if `"`0'"' == "," {
		sret local stop 1
		exit
	}
	if `"`0'"' == "if" {
		sret local stop 1
		exit
	}
* per official ivreg 5.1.3
	if substr(`"`0'"',1,3) == "if(" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "in" {
		sret local stop 1
		exit
	}
	if `"`0'"' == "" {
		sret local stop 1
		exit
	}
	else	sret local stop 0
end

// ************* Display list of variables ************ //

program define Disp 
	version 11.2
	syntax [anything] [, _col(integer 15) ]
	local maxlen = 80-`_col'
	local len = 0
	local first = 1
	foreach vn in `anything' {
* Don't display if base or omitted variable
		_ms_parse_parts `vn'
		if ~`r(omit)' {
			local vnlen		: length local vn
			if `len'+`vnlen' > `maxlen' {
				di
				local first = 1
				local len = `vnlen'
			}
			else {
				local len = `len'+`vnlen'+1
			}
			if `first' {
				local first = 0
				di in gr _col(`_col') "`vn'" _c
				}
			else {
				di in gr " `vn'" _c
			}
		}
	}
* Finish with a newline
	di
end

// *********** Display Cragg-Donald/Stock-Yogo critical values etc. ******** //

program define Disp_cdsy, rclass
	version 11.2
	syntax , col1(integer) col2(integer) model(string) k2(integer) nendog(integer) [ fuller(string) ]
	local cdmissing=1
	if "`model'"=="iv" | "`model'"=="gmm2s" | "`model'"=="gmmw" {
		cdsy, type(ivbias5) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') " 5% maximal IV relative bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(ivbias10) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "10% maximal IV relative bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(ivbias20) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "20% maximal IV relative bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(ivbias30) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "30% maximal IV relative bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(ivsize10) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "10% maximal IV size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(ivsize15) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "15% maximal IV size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(ivsize20) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "20% maximal IV size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(ivsize25) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "25% maximal IV size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
	}
	else if ("`model'"=="liml" & "`fuller'"=="") | "`model'"=="cue" {
		cdsy, type(limlsize10) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "10% maximal LIML size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(limlsize15) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "15% maximal LIML size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(limlsize20) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "20% maximal LIML size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(limlsize25) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "25% maximal LIML size" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
	}
	else if ("`model'"=="liml" & "`fuller'"~="") {
		cdsy, type(fullrel5) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') " 5% maximal Fuller rel. bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(fullrel10) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "10% maximal Fuller rel. bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(fullrel20) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "20% maximal Fuller rel. bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(fullrel30) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "30% maximal Fuller rel. bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(fullmax5) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') " 5% Fuller maximum bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(fullmax10) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "10% Fuller maximum bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(fullmax20) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "20% Fuller maximum bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		cdsy, type(fullmax30) k2(`k2') nendog(`nendog')
		if "`r(cv)'"~="." {
			di in gr _col(`col1') "30% Fuller maximum bias" in ye _col(`col2') %6.2f r(cv)
			local cdmissing=0
		}
		di in gr "NB: Critical values based on Fuller parameter=1"
	}	
	return scalar cdmissing	=`cdmissing'
end

program define cdsy, rclass
	version 11.2
	syntax , type(string) k2(integer) nendog(integer)

* type() can be ivbias5   (k2<=100, nendog<=3)
*               ivbias10  (ditto)
*               ivbias20  (ditto)
*               ivbias30  (ditto)
*               ivsize10  (k2<=100, nendog<=2)
*               ivsize15  (ditto)
*               ivsize20  (ditto)
*               ivsize25  (ditto)
*               fullrel5  (ditto)
*               fullrel10 (ditto)
*               fullrel20 (ditto)
*               fullrel30 (ditto)
*               fullmax5  (ditto)
*               fullmax10 (ditto)
*               fullmax20 (ditto)
*               fullmax30 (ditto)
*               limlsize10 (ditto)
*               limlsize15 (ditto)
*               limlsize20 (ditto)
*               limlsize25 (ditto)

tempname temp cv

* Initialize critical value as MV
	scalar `cv'=.

if "`type'"=="ivbias5" {
		mata: s_cdsy("`temp'", 1)
		if `k2'<=100 & `nendog'<=3 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="ivbias10" {
 		mata: s_cdsy("`temp'", 2)
		if `k2'<=100 & `nendog'<=3 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="ivbias20" {
	 	mata: s_cdsy("`temp'", 3)
		if `k2'<=100 & `nendog'<=3 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="ivbias30" {
	 	mata: s_cdsy("`temp'", 4)
		if `k2'<=100 & `nendog'<=3 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="ivsize10" {
	 	mata: s_cdsy("`temp'", 5)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="ivsize15" {
	 	mata: s_cdsy("`temp'", 6)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="ivsize20" {
	 	mata: s_cdsy("`temp'", 7)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="ivsize25" {
	 	mata: s_cdsy("`temp'", 8)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}
	
	if "`type'"=="fullrel5" {
	 	mata: s_cdsy("`temp'", 9)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}
	
	if "`type'"=="fullrel10" {
	 	mata: s_cdsy("`temp'", 10)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}
	
	if "`type'"=="fullrel20" {
	 	mata: s_cdsy("`temp'", 11)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="fullrel30" {
	 	mata: s_cdsy("`temp'", 12)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="fullmax5" {
	 	mata: s_cdsy("`temp'", 13)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}
	
	if "`type'"=="fullmax10" {
	 	mata: s_cdsy("`temp'", 14)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}
	
	if "`type'"=="fullmax20" {
	 	mata: s_cdsy("`temp'", 15)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="fullmax30" {
	 	mata: s_cdsy("`temp'", 16)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="limlsize10" {
	 	mata: s_cdsy("`temp'", 17)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="limlsize15" {
	 	mata: s_cdsy("`temp'", 18)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="limlsize20" {
	 	mata: s_cdsy("`temp'", 19)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}

if "`type'"=="limlsize25" {
	 	mata: s_cdsy("`temp'", 20)
		if `k2'<=100 & `nendog'<=2 {
			scalar `cv'=`temp'[`k2',`nendog']
		}
	}
	return scalar cv=`cv'
end

// ***************************** Parse ivreg2 arguments **************** //

program define ivparse, sclass
	version 11.2
		syntax [anything(name=0)]			///
			[ ,								///
				ivreg2name(name)			///
				partial(string)				///  as string because may have nonvariable in list
				fwl(string)					///  legacy option
				orthog(varlist fv ts)		///
				endogtest(varlist fv ts)	///
				redundant(varlist fv ts)	///
				depname(string)				///
				robust						///
				cluster(varlist fv ts)		///
				bw(string)					/// as string because may have noninteger option "auto"
				kernel(string)				///
				dkraay(integer 0)			///
				sw							///
				kiefer						///
				center						///
				NOCONSTANT					///
				tvar(varname)				///
				ivar(varname)				///
				gmm2s						///
				gmm							///
				cue							///
				liml						///
				fuller(real 0)				///
				kclass(real 0)				///
				b0(string)					///
				wmatrix(string)				///
				NOID						///
				savefirst					///
				savefprefix(name)			///
				saverf						///
				saverfprefix(name)			///
				savesfirst					///
				savesfprefix(name)			///
				psd0						///
				psda						///
				dofminus(integer 0)			///
				NOCOLLIN					///
				useqr						///
				bvclean						///
				eform(string)				///
				NOOMITTED					///
				vsquish						///
				noemptycells				///
				baselevels					///
				allbaselevels 				///
			]

// TS and FV opts based on option varlists
		local tsops		= ("`s(tsops)'"=="true")
		local fvops		= ("`s(fvops)'"=="true")
// useful boolean
		local cons		=("`noconstant'"=="")

local n 0
		gettoken lhs 0 : 0, parse(" ,[") match(paren)
		IsStop `lhs'
		while `s(stop)'==0 {
			if "`paren'"=="(" {
				local ++n
				if `n'>1 { 
di as err `"syntax is "(all instrumented variables = instrument variables)""'
					exit 198
				}
				gettoken p lhs : lhs, parse(" =")
				while "`p'"!="=" {
					if "`p'"=="" {
di as err `"syntax is "(all instrumented variables = instrument variables)""'
di as er `"the equal sign "=" is required"'
						exit 198
					}
					local endo `endo' `p'
					gettoken p lhs : lhs, parse(" =")
				}
				local exexog `lhs'
			}
			else {
				local inexog `inexog' `lhs'
			}
			gettoken lhs 0 : 0, parse(" ,[") match(paren)
			IsStop `lhs'
		}
// lhs attached to front of inexog
		gettoken lhs inexog	: inexog
		local endo			: list retokenize endo
		local inexog		: list retokenize inexog
		local exexog		: list retokenize exexog
// If depname not provided (default) name is lhs variable
		if "`depname'"=="" {
			local depname `lhs'
		}

// partial, including legacy FWL option
		local partial		`partial' `fwl'
// Need to nonvars "_cons" from list if present
// Also set `partialcons' local to 0/1
// Need word option so that varnames with cons in them aren't zapped
		local partial		: subinstr local partial "_cons" "", all count(local partialcons) word
		local partial		: list retokenize partial
		if "`partial'"=="_all" {
			local partial	`inexog'
		}
// constant always partialled out if present in regression and other inexog are being partialled out
// (incompatibilities caught in error-check section below)
		if "`partial'"~="" {
			local partialcons	= (`cons' | `partialcons')
		}

// detect if TS or FV operators used in main varlists
// clear any extraneous sreturn macros first
		sreturn clear
		local 0				`lhs' `inexog' `endo' `exexog' `partial'
		syntax				varlist(fv ts)
		local tsops			= ("`s(tsops)'"=="true") | `tsops'
		local fvops			= ("`s(fvops)'"=="true") | `fvops'

// TS operators not allowed with cluster, ivar or tvar.  Captured in -syntax-.
		if "`tvar'" == "" {
			local tvar		`_dta[_TStvar]'
		}
		if "`ivar'" == "" {
			local ivar		`_dta[_TSpanel]'
		}
		if "`_dta[_TSdelta]'" == "" {
			local tdelta	1
		}
		else {												//  use evaluator since _dta[_TSdelta] can
			local tdelta	= `_dta[_TSdelta]'				//  be stored as e.g. +1.0000000000000X+000
		}

sreturn local lhs			`lhs'
		sreturn local depname		`depname'
		sreturn local endo			`endo'
		sreturn local inexog		`inexog'
		sreturn local exexog 		`exexog'
		sreturn local partial		`partial'
		sreturn local cons			=`cons'
		sreturn local partialcons	=`partialcons'
		sreturn local tsops			=`tsops'
		sreturn local fvops			=`fvops'
		sreturn local tvar			`tvar'
		sreturn local ivar			`ivar'
		sreturn local tdelta		`tdelta'
		sreturn local noid			`noid'			//  can be overriden below
		sreturn local liml			`liml'			//  can be overriden below

//convert to boolean
		sreturn local useqr			=("`useqr'" ~= "")

// Cluster and SW imply robust
		if "`cluster'`sw'"~="" {
			local robust	"robust"
		}

// HAC estimation.

// First dkraay(bw): special case of HAC with clustering
// on time-series var in a panel + kernel-robust
		if `dkraay' {
			if "`bw'" == "" {
				local bw	`dkraay'
			}
			if "`cluster'" == "" {
				local cluster	`tvar'
			}
		}
// If bw is omitted, default `bw' is 0.
// bw(.) can be number or "auto" hence arrives as string, but is returned as number
// bw=-1 returned if "auto"
// If bw or kernel supplied, check/set `kernel'.
// Macro `kernel' is also used for indicating HAC in use.
// If bw or kernel not supplied, set bw=0
		if "`bw'" == "" & "`kernel'" == "" {
			local bw	0
		}
		else {
// Check it's a valid kernel and replace with unabbreviated kernel name; check bw.
// s_vkernel is in livreg2 mlib.
			mata: s_vkernel("`kernel'", "`bw'", "`ivar'")
			local kernel	`r(kernel)'
			local bw		`r(bw)'			//  = -1 if bw(auto) option chosen
			local tsops		= 1
		}
// kiefer = kernel(tru) bw(T) and no robust
		if "`kiefer'" ~= "" & "`kernel'" == "" {
			local kernel "Truncated"
		}

// Done parsing VCE opts
		sreturn local bw		`bw'
		sreturn local kernel	`kernel'
		sreturn local robust	`robust'
		sreturn local cluster	`cluster'
		if `bw' {
			sreturn local bwopt "bw(`bw')"
			sreturn local kernopt "kernel(`kernel')"
		}
// center arrives as string but is returned as boolean
		sreturn local center	=("`center'"=="center")

// Fuller implies LIML
		if `fuller' != 0 {
			sreturn local liml			"liml"
			sreturn local fulleropt		"fuller(`fuller')"
		}

if `kclass' != 0 {
			sreturn local kclassopt		"kclass(`kclass')"
		}

// b0 implies noid.
		if "`b0'" ~= "" {
			sreturn local noid			"noid"
		}

// save first, rf
		if "`savefprefix'" != "" {						//  savefprefix implies savefirst
			local savefirst				"savefirst"
		}
		else {											//  default savefprefix is _ivreg2_
			local savefprefix			"_`ivreg2name'_"
		}
		sreturn local savefirst			`savefirst'
		sreturn local savefprefix		`savefprefix'
		if "`saverfprefix'" != "" {						//  saverfprefix implies saverf
			local saverf				"saverf"
		}
		else {											// default saverfprefix is _ivreg2_
			local saverfprefix			"_`ivreg2name'_"
		}
		sreturn local saverf			`saverf'
		sreturn local saverfprefix		`saverfprefix'
		if "`savesfprefix'" != "" {					//  savesfprefix implies savesfirst
			local savesfirst			"savesfirst"
		}
		else {											// default saverfprefix is _ivreg2_
			local savesfprefix			"_`ivreg2name'_"
		}
		sreturn local savesfirst		`savesfirst'
		sreturn local savesfprefix		`savesfprefix'

// Macro psd has either psd0, psda or is empty
		sreturn local psd		"`psd0'`psda'"

// dofminus
		if `dofminus' {
			sreturn local dofmopt	dofminus(`dofminus')
		}

// display options
		local dispopt			eform(`eform') `vsquish' `noomitted' `noemptycells' `baselevels' `allbaselevels'
// now boolean - indicates that omitted and/or base vars should NOT be added to VCV
// automatically triggered by partial
		local bvclean			= wordcount("`bvclean'") | wordcount("`partial'") | `partialcons'
		sreturn local bvclean	`bvclean'
		sreturn local dispopt	`dispopt'

// ************ ERROR CHECKS ************* //

if `partialcons' & ~`cons' {
di in r "Error: _cons listed in partial() but equation specifies -noconstant-." 
			exit 198
		}
		if `partialcons' > 1 {
// Just in case of multiple _cons
di in r "Error: _cons listed more than once in partial()." 
			exit 198
		}

// User-supplied tvar and ivar checked if consistent with tsset.
		if "`tvar'"!="`_dta[_TStvar]'" {
di as err "invalid tvar() option - data already -tsset-"
			exit 5
		}
		if "`ivar'"!="`_dta[_TSpanel]'" {
di as err "invalid ivar() option - data already -xtset-"
			exit 5
		}

// dkraay
		if `dkraay' {
			if "`ivar'" == "" | "`tvar'" == "" {
di as err "invalid use of dkraay option - must use tsset panel data"
				exit 5
			}
			if "`dkraay'" ~= "`bw'" {
di as err "cannot use dkraay(.) and bw(.) options together"
				exit 198
			}
			if "`cluster'" ~= "`tvar'" {
di as err "invalid use of dkraay option - must cluster on `tvar' (or omit cluster option)"
				exit 198
			}
		}

// kiefer VCV = kernel(tru) bw(T) and no robust with tsset data
		if "`kiefer'" ~= "" {
			if "`ivar'" == "" | "`tvar'" == "" {
di as err "invalid use of kiefer option - must use tsset panel data"
				exit 5
			}
			if	"`robust'" ~= "" {
di as err "incompatible options: kiefer and robust"
				exit 198
			}
			if	"`kernel'" ~= "" & "`kernel'" ~= "Truncated" {
di as err "incompatible options: kiefer and kernel(`kernel')"
				exit 198
			}
			if	(`bw'~=0) {
di as err "incompatible options: kiefer and bw"
				exit 198
			}
		}

// sw=Stock-Watson robust SEs
		if "`sw'" ~= "" & "`cluster'" ~= "" {
di as err "Stock-Watson robust SEs not supported with -cluster- option"
				exit 198
		}
		if "`sw'" ~= "" & "`kernel'" ~= "" {
di as err "Stock-Watson robust SEs not supported with -kernel- option"
				exit 198
		}
		if "`sw'" ~= "" & "`ivar'"=="" {
di as err "Must -xtset- or -tsset- data or specify -ivar- with -sw- option"
			exit 198
		}

// LIML/kclass incompatibilities
		if "`liml'`kclassopt'" != "" {
			if "`gmm2s'`cue'" != "" {
di as err "GMM estimation not available with LIML or k-class estimators"
			exit 198
			}
			if `fuller' < 0 {
di as err "invalid Fuller option"
			exit 198
			}
			if "`liml'" != "" & "`kclassopt'" != "" {
di as err "cannot use liml and kclass options together"
			exit 198
			}
			if `kclass' < 0 {
di as err "invalid k-class option"
				exit 198
				}			
		}

if "`gmm2s'" != "" & "`cue'" != "" {
di as err "incompatible options: 2-step efficient gmm and cue gmm"
			exit 198
		}

if "`gmm2s'`cue'" != "" & "`exexog'" == "" {
di as err "option `gmm2s'`cue' invalid: no excluded instruments specified"
			exit 102
		}

// Legacy gmm option
		if "`gmm'" ~= "" {
di as err "-gmm- is no longer a supported option; use -gmm2s- with the appropriate option"
di as res "      gmm             =  gmm2s robust"
di as res "      gmm robust      =  gmm2s robust"
di as res "      gmm bw()        =  gmm2s bw()"
di as res "      gmm robust bw() =  gmm2s robust bw()"
di as res "      gmm cluster()   =  gmm2s cluster()"
			exit 198
		}
		
// b0 incompatible options.
		if "`b0'" ~= "" & "`gmm2s'`cue'`liml'`wmatrix'" ~= "" {
di as err "incompatible options: -b0- and `gmm2s' `cue' `liml' `wmatrix'"
			exit 198
		}
		if "`b0'" ~= "" & `kclass' ~= 0 {
di as err "incompatible options: -b0- and kclass(`kclass')"
			exit 198
		}

if "`psd0'"~="" & "`psda'"~="" {
di as err "cannot use psd0 and psda options together"
			exit 198
		}
end

// *************** Check varlists for for duplicates and collinearities ***************** //

program define CheckDupsCollin, sclass
	version 11.2
		syntax 								///
			[ ,								///
				lhs(string)					///
				endo(string)				///
				inexog(string)				///
				exexog(string)				///
				partial(string)				///
				orthog(string)				///
				endogtest(string)			///
				redundant(string)			///
				touse(string)				///
				wvar(string)				///
				wf(real 0)					///
				NOCONSTANT					///
				NOCOLLIN					///
				fvall						///
				fvsep						///
			]

if "`fvall'`fvsep'"=="" {				//  default, expand RHS and exexog separately
			local rhs	`endo' `inexog'
			foreach vl in lhs rhs exexog {
				fvexpand ``vl'' if `touse'
				local `vl' `r(varlist)'
			}
			local allvars	`rhs' `exexog'
		}
		else if "`fvall'"~="" {					//  expand all 3 varlists as one
			fvexpand `lhs' if `touse'
			local lhs		`r(varlist)'
			fvexpand `endo' `inexog' `exexog' if `touse'
			local allvars	`r(varlist)'
		}
		else if "`fvsep'"~="" {					//  expand 3 varlists separately
			foreach vl in lhs endo inexog exexog {
				fvexpand ``vl'' if `touse'
				local `vl' `r(varlist)'
			}
			local allvars	`endo' `inexog' `exexog'
		}
		else {									//  shouldn't reach here
di as err "internal ivreg2 err: CheckDupsCollin"
			exit 198
		}

// Create dictionary: `allvars' is list with b/n/o etc., sallvars is stripped version
// NB: lhs is not in dictionary and won't need to recreate it
		ivreg2_fvstrip `allvars'
		local sallvars	`r(varlist)'

// Create consistent expanded varlists
// (1) expand; (2) strip (since base etc. may be wrong); (3) recreate using dictionary
// NB: matchnames will return unmatched original name if not found in 2nd arg varlist
		foreach vl in endo inexog exexog partial orthog endogtest redundant {
			fvexpand ``vl'' if `touse'
			ivreg2_fvstrip `r(varlist)'
			local stripped	`r(varlist)'							//  create stripped version of varlist
			matchnames "`stripped'" "`sallvars'" "`allvars'"		//  match using dictionary
			local `vl'		`r(names)'								//  new consistent varlist with correct b/n/o etc.
		}

// Check for duplicates of variables
// (1)  inexog > endo
// (2)  inexog > exexog
// (3)  endo + exexog = inexog, as if it were "perfectly predicted"
		local lhs0		`lhs'					//  create here
		local dupsen1	: list dups endo
		local dupsin1	: list dups inexog
		local dupsex1	: list dups exexog
		foreach vl in endo inexog exexog partial orthog endogtest redundant {
			local `vl'0	: list uniq `vl'
		}
// Remove inexog from endo
		local dupsen2	: list endo0 & inexog0
		local endo0		: list endo0 - inexog0
// Remove inexog from exexog
		local dupsex2	: list exexog0 & inexog0
		local exexog0	: list exexog0 - inexog0
// Remove endo from exexog
		local dupsex3	: list exexog0 & endo0
		local exexog0	: list exexog0 - endo0
		local dups "`dupsen1' `dupsex1' `dupsin1' `dupsen2' `dupsex2' `dupsex3'"
		local dups		: list uniq dups

// Collinearity checks

// Need variable counts for "0" varlists
// These do NOT include the constant
		local endo0_ct		: word count `endo0'
		local inexog0_ct	: word count `inexog0'
		local rhs0_ct		: word count `inexog0' `exexog0'
		local exexog0_ct	: word count `exexog'

if "`nocollin'" == "" {

// Needed for ivreg2_rmcollright2
			tempvar normwt
			qui gen double `normwt' = `wf' * `wvar' if `touse'

// Simple case: no endogenous regressors, only included and excluded exogenous
			if `endo0_ct'==0 {
// Call ivreg2_rmcollright2 on "0" versions of inexog and exexog
// noexpand since already expanded and don't want inconsistant expansion
// newonly since don't want base vars in collinear list
				qui ivreg2_rmcollright2 `inexog0' `exexog0' if `touse',		///
						normwt(`normwt') `noconstant' noexpand newonly
// ivreg2_rmcollright2 returns fulll varlist with omitteds marked as omitted,
// so just need to separate the inexog and exexog lists
				if `r(k_omitted)' {
					local collin	`collin' `r(omitted)'
					local inexog0	""
					local exexog0	""
					local nvarlist	`r(varlist)'
					local i			1
					while `i' <= `rhs0_ct' {
						local nvar		: word `i' of `nvarlist'
						if `i' <= `inexog0_ct' {
							local inexog0	`inexog0' `nvar'		//  first batch go into inexog0
						}
						else {
							local exexog0	`exexog0' `nvar'		//  remainder go into exexog0
						}
						local ++i
					}
					local inexog0	: list retokenize inexog0
					local exexog0	: list retokenize exexog0
				}
			}
// Not-simple case: endogenous regressors
			else {

// 1st pass through - remove intra-endo collinears
				qui ivreg2_rmcollright2 `endo0' if `touse',				///
						normwt(`normwt') `noconstant' noexpand newonly
// ivreg2_rmcollright2 returns fulll varlist with omitteds marked as omitted,
// so just need to separate the inexog and exexog lists
				if `r(k_omitted)' {
					local collin	`collin' `r(omitted)'
					local endo0		`r(varlist)'
				}

// 2nd pass through - good enough unless endog appear as colllinear
// noexpand since already expanded and don't want inconsistent expansion
// newonly since don't want base vars in collinear list
				qui ivreg2_rmcollright2 `inexog0' `exexog0' `endo0' if `touse',		///
						normwt(`normwt') `noconstant' noexpand newonly
				if `r(k_omitted)' {
// Check if any endo are in the collinears.
// If yes, reclassify as inexog, then
// 3rd pass through - and then proceed to process inexog and exexog as above
					local ecollin	`r(omitted)'
					local ecollin	: list ecollin - inexog0
					local ecollin	: list ecollin - exexog0
					if wordcount("`ecollin'") {
// Collinears in endo, so reclassify as inexog, redo counts and call ivreg2_rmcollright2 again
						local endo0			: list endo0 - ecollin
						local inexog0		`ecollin' `inexog0'
						local inexog0		: list retokenize inexog0
						local endo0_ct		: word count `endo0'
						local inexog0_ct	: word count `inexog0'
						local rhs0_ct		: word count `inexog0' `exexog0'
// noexpand since already expanded and don't want inconsistant expansion
// newonly since don't want base vars in collinear list
						qui ivreg2_rmcollright2 `inexog0' `exexog0' `endo0' if `touse',		///
								normwt(`normwt') `noconstant' noexpand newonly
					}
// Collinears in inexog or exexog
					local collin	`collin' `r(omitted)'
					local inexog0	""
					local exexog0	""
					local nvarlist	`r(varlist)'
					local i			1
					while `i' <= `rhs0_ct' {
						local nvar		: word `i' of `nvarlist'
						if `i' <= `inexog0_ct' {
							local inexog0	`inexog0' `nvar'
						}
						else {
							local exexog0	`exexog0' `nvar'
						}
						local ++i
					}
					local inexog0	: list retokenize inexog0
					local exexog0	: list retokenize exexog0
				}
			}

// Collinearity and duplicates warning messages, if necessary
			if "`dups'" != "" {
di in gr "Warning - duplicate variables detected"
di in gr "Duplicates:" _c
				Disp `dups', _col(21)
			}
			if "`ecollin'" != "" {
di in gr "Warning - endogenous variable(s) collinear with instruments"
di in gr "Vars now exogenous:" _c
				Disp `ecollin', _col(21)
			}
			if "`collin'" != "" {
di in gr "Warning - collinearities detected"
di in gr "Vars dropped:" _c
				Disp `collin', _col(21)
			}
		}

// Last step: process partial0 so that names with o/b/n etc. match inexog0
		if wordcount("`partial0'") {
			ivreg2_fvstrip `inexog0' if `touse'
			local sinexog0		`r(varlist)'							//  for inexog dictionary
			ivreg2_fvstrip `partial0' if `touse'
			local spartial0		`r(varlist)'							//  for partial dictionary
			matchnames "`spartial0'" "`sinexog0'" "`inexog0'"			//  match using dictionary
			local partial0	`r(names)'									//  new partial0 with matches
			local partialcheck		: list partial0 - inexog0			//  unmatched are still in partial0
			if ("`partialcheck'"~="") {									//  so catch them
di in r "Error: `partialcheck' listed in partial() but not in list of regressors." 
				error 198
			}
		}
// Completed duplicates and collinearity checks

foreach vl in lhs endo inexog exexog partial orthog endogtest redundant {
			sreturn local `vl'		``vl''
			sreturn local `vl'0		``vl'0'
		}
		sreturn local dups		`dups'
		sreturn local collin	`collin'
		sreturn local ecollin	`ecollin'

end

// ******************* Misc error checks *************************** //

program define CheckMisc, rclass
	version 11.2
		syntax 								///
			[ ,								///
				rhs1_ct(integer 0)			///
				iv1_ct(integer 0)			///
				bvector(name)				///
				smatrix(name)				///
				wmatrix(name)				///
				cnb1(string)				///
				cnZ1(string)				///
			]

// Check variable counts
		if `rhs1_ct' == 0 {
di as err "error: no regressors specified"
			exit 102
		}
		if `rhs1_ct' > `iv1_ct' {
di as err "equation not identified; must have at least as many instruments"
di as err "not in the regression as there are instrumented variables"
			exit 481
		}

// Check user-supplied b vector
		if "`bvector'" != "" {
			tempname b0
// Rearrange/select columns to mat IV matrix
			cap matsort `bvector' "`cnb1'"
			matrix `b0'=r(sorted)
			local scols = colsof(`b0')
			local bcols : word count `cnb1'
			if _rc ~= 0 | (`scols'~=`bcols') {
di as err "-b0- option error: supplied b0 columns do not match regressor list"
exit 198
			}
			return mat b0 = `b0'
		}

// Check user-supplied S matrix
		if "`smatrix'" != "" {
			tempname S0
// Check that smatrix is indeed a matrix
			cap mat S0 = `smatrix'
			if _rc ~= 0 {
di as err "invalid matrix `smatrix' in smatrix option"
exit _rc
			}
// Rearrange/select columns to mat IV matrix
			cap matsort `smatrix' "`cnZ1'"
			matrix `S0'=r(sorted)
			local srows = rowsof(`S0')
			local scols = colsof(`S0')
			local zcols : word count `cnZ1'
			if _rc ~= 0 | (`srows'~=`zcols') | (`scols'~=`zcols') {
di as err "-smatrix- option error: supplied matrix columns/rows do not match IV list"
exit 198
			}
			if issymmetric(`S0')==0 {
di as err "-smatrix- option error: supplied matrix is not symmetric"
exit 198
			}
			return mat S0 = `S0'
		}

// Check user-supplied W matrix
		if "`wmatrix'" != "" {
			tempname W0
// Check that wmatrix is indeed a matrix
			cap mat W0 = `wmatrix'
			if _rc ~= 0 {
di as err "invalid matrix `wmatrix' in wmatrix option"
exit _rc
			}
// Rearrange/select columns to mat IV matrix
			cap matsort `wmatrix' "`cnZ1'"
			matrix `W0'=r(sorted)
			local srows = rowsof(`W0')
			local scols = colsof(`W0')
			local zcols : word count `cnZ1'
			if _rc ~= 0 | (`srows'~=`zcols') | (`scols'~=`zcols') {
di as err "-wmatrix- option error: supplied matrix columns/rows do not match IV list"
exit 198
			}
			if issymmetric(`W0')==0 {
di as err "-smatrix- option error: supplied matrix is not symmetric"
exit 198
			}
			return mat W0 = `W0'
		}
end

*******************************************************************************
************************* misc utilities **************************************
*******************************************************************************

// internal version of ivreg2_fvstrip 1.01 ms 24march2015
// takes varlist with possible FVs and strips out b/n/o notation
// returns results in r(varnames)
// optionally also omits omittable FVs
// expand calls fvexpand either on full varlist
// or (with onebyone option) on elements of varlist

program define ivreg2_fvstrip, rclass
	version 11.2
	syntax [anything] [if] , [ dropomit expand onebyone NOIsily ]
	if "`expand'"~="" {												//  force call to fvexpand
		if "`onebyone'"=="" {
			fvexpand `anything' `if'								//  single call to fvexpand
			local anything `r(varlist)'
		}
		else {
			foreach vn of local anything {
				fvexpand `vn' `if'									//  call fvexpand on items one-by-one
				local newlist	`newlist' `r(varlist)'
			}
			local anything	: list clean newlist
		}
	}
	foreach vn of local anything {									//  loop through varnames
		if "`dropomit'"~="" {										//  check & include only if
			_ms_parse_parts `vn'									//  not omitted (b. or o.)
			if ~`r(omit)' {
				local unstripped	`unstripped' `vn'				//  add to list only if not omitted
			}
		}
		else {														//  add varname to list even if
			local unstripped		`unstripped' `vn'				//  could be omitted (b. or o.)
		}
	}
// Now create list with b/n/o stripped out
	foreach vn of local unstripped {
		local svn ""											//  initialize
		_ms_parse_parts `vn'
		if "`r(type)'"=="variable" & "`r(op)'"=="" {			//  simplest case - no change
			local svn	`vn'
		}
		else if "`r(type)'"=="variable" & "`r(op)'"=="o" {		//  next simplest case - o.varname => varname
			local svn	`r(name)'
		}
		else if "`r(type)'"=="variable" {						//  has other operators so strip o but leave .
			local op	`r(op)'
			local op	: subinstr local op "o" "", all
			local svn	`op'.`r(name)'
		}
		else if "`r(type)'"=="factor" {							//  simple factor variable
			local op	`r(op)'
			local op	: subinstr local op "b" "", all
			local op	: subinstr local op "n" "", all
			local op	: subinstr local op "o" "", all
			local svn	`op'.`r(name)'							//  operator + . + varname
		}
		else if"`r(type)'"=="interaction" {						//  multiple variables
			forvalues i=1/`r(k_names)' {
				local op	`r(op`i')'
				local op	: subinstr local op "b" "", all
				local op	: subinstr local op "n" "", all
				local op	: subinstr local op "o" "", all
				local opv	`op'.`r(name`i')'					//  operator + . + varname
				if `i'==1 {
					local svn	`opv'
				}
				else {
					local svn	`svn'#`opv'
				}
			}
		}
		else if "`r(type)'"=="product" {
			di as err "ivreg2_fvstrip error - type=product for `vn'"
			exit 198
		}
		else if "`r(type)'"=="error" {
			di as err "ivreg2_fvstrip error - type=error for `vn'"
			exit 198
		}
		else {
			di as err "ivreg2_fvstrip error - unknown type for `vn'"
			exit 198
		}
		local stripped `stripped' `svn'
	}
	local stripped	: list retokenize stripped						//  clean any extra spaces
	
	if "`noisily'"~="" {											//  for debugging etc.
di as result "`stripped'"
	}

return local varlist	`stripped'								//  return results in r(varlist)
end

// **************** Add omitted vars to b  and V matrices ****************** //

program define AddOmitted, rclass
	version 11.2
	syntax 								///
		[ ,								///
			bmat(name)					///
			vmat(name)					///
			cnb0(string)				///
			cnb1(string)				///
			eqlist(string)				/// if empty, single-equation b and V
		]

tempname newb newV
	local eq_ct		=max(1,wordcount("`eqlist'"))
	local rhs0_ct	: word count `cnb0'
	local rhs1_ct	: word count `cnb1'

foreach vn in `cnb1' {
		local cnum	: list posof "`vn'" in cnb0
		local cnumlist "`cnumlist' `cnum'"
	}
// cnumlist is the list of columns in the single-equation new big matrix in which
// the non-zero entries from the reduced matrix (bmat or vmat) will appear.
// E.g., if newb will be [mpg o.mpg2 _cons] then cnum = [1 3].

mata: s_AddOmitted(					///
						"`bmat'",		///
						"`vmat'",		///
						"`cnumlist'",	///
						`eq_ct',		///
						`rhs0_ct',		///
						`rhs1_ct')
	mat `newb' = r(b)
	mat `newV' = r(V)

if `eq_ct'==1 {
		local allnames	`cnb0'						//  simple single-eqn case
	}
	else {
		foreach eqname in `eqlist' {
			foreach vname in `cnb0' {
				local allnames	"`allnames' `eqname':`vname'"
			}
		}
	}
	mat colnames `newb'	= `allnames'
	mat rownames `newb' = y1
	mat colnames `newV' = `allnames'
	mat rownames `newV' = `allnames'

return matrix b		=`newb'
	return matrix V 	=`newV'
end

// ************* More misc utilities ************** //

program define matsort, rclass
	version 11.2
	args bvmat names
	tempname m1 m2
	foreach vn in `names' {
		mat `m1'=nullmat(`m1'), `bvmat'[1...,"`vn'"]
	}
	if rowsof(`m1')>1 {
		foreach vn in `names' {
			mat `m2'=nullmat(`m2') \ `m1'["`vn'",1...]
		}
		return matrix sorted =`m2'
	}
	else {
		return matrix sorted =`m1'
	}
end

program define matchnames, rclass
	version 11.2
	args	varnames namelist1 namelist2

local k1 : word count `namelist1'
	local k2 : word count `namelist2'

if `k1' ~= `k2' {
		di as err "namelist error"
		exit 198
	}
	foreach vn in `varnames' {
		local i : list posof `"`vn'"' in namelist1
		if `i' > 0 {
			local newname : word `i' of `namelist2'
		}
		else {
* Keep old name if not found in list
			local newname "`vn'"
		}
		local names "`names' `newname'"
	}
	local names	: list clean names
	return local names "`names'"
end

program define checkversion_ranktest
	version 11.2
	args ranktestcmd ranktestversion

* Check that -ranktest- is installed
		capture `ranktestcmd', version
		if _rc != 0 {
di as err "Error: must have ranktest version `ranktestversion' or greater installed"
di as err "To install, from within Stata type " _c
di in smcl "{stata ssc install ranktest :ssc install ranktest}"
			exit 601
		}
		local vernum "`r(version)'"
		if ("`vernum'" < "`ranktestversion'") | ("`vernum'" > "09.9.99") {
di as err "Error: must have ranktest version `ranktestversion' or greater installed"
di as err "Currently installed version is `vernum'"
di as err "To update, from within Stata type " _c
di in smcl "{stata ssc install ranktest, replace :ssc install ranktest, replace}"
			exit 601
		}
end

// ************ Replacement _rmcollright with tweaks ****************** //

program define ivreg2_rmcollright2, rclass
	version 11.2
	syntax	[ anything ]					///  anything so that FVs aren't reordered
			[if] [in]						///
			[, 								///
			NORMWT(varname)					///
			NOCONStant						///
			NOEXPAND						///
			newonly							///
			lindep							///
			]

// Empty varlist, leave early
	if "`anything'"=="" {
		return scalar k_omitted =0
		exit
	}

marksample touse
	markout `touse' `anything'

local cons		= ("`noconstant'"=="")
	local expand	= ("`noexpand'"=="")
	local newonly	= ("`newonly'"~="")
	local forcedrop	= ("`forcedrop'"~="")
	local lindep	= ("`lindep'"~="")
	local 0			`anything'
	sreturn clear								//  clear any extraneous sreturn macros
	syntax varlist(ts fv)
	local tsops		= ("`s(tsops)'"=="true")
	local fvops		= ("`s(fvops)'"=="true")

if `tsops' | `fvops' {
		if `expand' {
			fvexpand `anything' if `touse'
			local anything		`r(varlist)'
			fvrevar `anything' if `touse'
			local fv_anything	`r(varlist)'
		}
		else {
// already expanded and in set order
// loop through fvrevar so that it doesn't rebase or reorder
			foreach var in `anything' {
				fvrevar `var' if `touse'
				local fv_anything	`fv_anything' `r(varlist)'
			}
		}
	}
	else {
		local fv_anything	`anything'
	}

tempname wname
	if "`normwt'"=="" {
		qui gen byte `wname'=1 if `touse'
	}
	else {
		qui gen double `wname' = `normwt' if `touse'
	}

mata: s_rmcoll2("`fv_anything'", "`anything'", "`wname'", "`touse'", `cons', `lindep')
	
	foreach var in `r(omitted)' {
		di as text "note: `var' omitted because of collinearity"
	}
	
	local omitted	"`r(omitted)'"			//  has all omitted, both newly and previously omitted
	local k_omitted	=r(k_omitted)			//  but newly omitted not marked with omit operator o
	if `lindep' {
		tempname lindepmat
		mat `lindepmat' = r(lindep)
		mat rownames `lindepmat' = `anything'
		mat colnames `lindepmat' = `anything'
	}

// Modern Stata version, add omitted notation to newly-missing vars
	if `k_omitted' {
		foreach var in `omitted' {
			_ms_parse_parts `var'			//  check if already omitted
			if r(omit) {					//  already omitted
				local alreadyomitted	`alreadyomitted' `var'
			}
			else {							//  not already omitted
				ivreg2_rmc2_ms_put_omit `var'	//  add omitted omit operator o and replace in main varlist
				local ovar	`s(ospec)'
				local anything	: subinstr local anything "`var'" "`ovar'", word
			}
		}
		if `newonly' {						//  omitted list should contain only newly omitted
			local omitted	: list omitted - alreadyomitted
			local k_omitted	: word count `omitted'
		}
	}

// Return results
	return scalar k_omitted =`k_omitted'
	return local omitted	`omitted'
	return local varlist	`anything'
	if `lindep' {
		return mat lindep	`lindepmat'
	}

end

// Used by ivreg2_rmcollright2
// taken from later Stata - not available in Stata 11
// version 1.0.0  28apr2011
program ivreg2_rmc2_ms_put_omit, sclass
	version 11.2						//  added by MS
        args vn
        _ms_parse_parts `vn'
        if r(type) =="variable" {
                local name `r(name)'
                local ovar o.`name'
        }
        if r(type) == "factor" {
                if !r(base) {
                        local name `r(name)'
                        if "`r(ts_op)'" != "" {
                                local name `r(ts_op)'.`name'
                        }       
                        local ovar `r(level)'o.`name'
                }
                else {
                        local ovar `vn'
                }
        }
        else if r(type) == "interaction" {
                local k = r(k_names)

forval i = 1/`k' {
                        local name = r(name`i')
                        if "`r(ts_op`i')'" != "" {
                                local name `r(ts_op`i')'.`name'
                        }
                        if "`r(level`i')'" != "" {
                                if r(base`i') {
                                        local name `r(level`i')'b.`name'
                                }
                                else {
                                        local name `r(level`i')'o.`name'
                                }
                        }
                        else {
                                local name o.`name'
                        }
                        local spec `spec'`sharp'`name'
                        local sharp "#"
                }
                local ovar `spec'
                                
        }
        _msparse `ovar'
        sreturn local ospec `r(stripe)'
end

*******************************************************************************
**************** SUBROUTINES FOR KERNEL-ROBUST ********************************
*******************************************************************************

// abw wants a varlist of [ eps | Z | touse]
// where Z includes all instruments, included and excluded, with constant if
// present as the last column; eps are a suitable set of residuals; and touse
// marks the observations in the data matrix used to generate the residuals
// (e.g. e(sample) of the appropriate model).
// The Noconstant option indicates that no constant term exists in the Z matrix.
// kern is the name of the HAC kernel. -ivregress- only provides definitions
// for Bartlett (default), Parzen, quadratic spectral.

// returns the optimal bandwidth as local abw

// abw 1.0.1  CFB 30jun2007
// 1.0.1 : redefine kernel names (3 instances) to match ivreg2
// 1.1.0 : pass nobs and tobs to s_abw; abw bug fix and also handles gaps in data correctly

prog def abw, rclass
	version 11.2
	syntax varlist(ts), [ tindex(varname) nobs(integer 0) tobs(integer 0) NOConstant Kernel(string)]
// validate kernel 
	if "`kernel'" == "" {
		local kernel = "Bartlett"
	}
// cfb B102
	if !inlist("`kernel'", "Bartlett", "Parzen", "Quadratic Spectral") {
		di as err "Error: kernel `kernel' not compatible with bw(auto)"
		return scalar abw = 1
		return local bwchoice "Kernel `kernel' not compatible with bw(auto); bw=1 (default)"
		exit
	}
	else {
// set constant
		local cons 1 
		if "`noconstant'" != "" {
			local cons 0
		}
// deal with ts ops 
		tsrevar `varlist'
		local varlist1 `r(varlist)'
		mata: s_abw("`varlist1'", "`tindex'", `nobs', `tobs', `cons', "`kernel'")
		return scalar abw = `abw'
		return local bwchoice "Automatic bw selection according to Newey-West (1994)"
	}
end

*******************************************************************************
************** END SUBROUTINES FOR KERNEL-ROBUST ******************************
*******************************************************************************

*******************************************************************************
*************************** BEGIN MATA CODE ***********************************
*******************************************************************************

// capture in case calling under version < 11.2
capture version 11.2

mata:

// For reference:
// struct ms_vcvorthog {
// 	string scalar	ename, Znames, touse, weight, wvarname
// 	string scalar	robust, clustvarname, clustvarname2, clustvarname3, kernel
// 	string scalar	sw, psd, ivarname, tvarname, tindexname
// 	real scalar		wf, N, bw, tdelta, dofminus
//  real scalar		center
// 	real matrix		ZZ
// 	pointer matrix	e
// 	pointer matrix	Z
// 	pointer matrix	wvar
// }

void s_abw	(string scalar Zulist,
			string scalar tindexname,
			real scalar nobs,
			real scalar tobs,
			real scalar cons,
			string scalar kernel
			)
{

// nobs = number of observations = number of data points available = rows(uZ)
// tobs = time span of data = t_N - t_1 + 1
// nobs = tobs if no gaps in data
// nobs < tobs if there are gaps
// nobs used below when calculating means, e.g., covariances in sigmahat.
// tobs used below when time span of data is needed, e.g., mstar.

string rowvector Zunames, tov
	string scalar v, v2
	real matrix uZ
	real rowvector h
	real scalar lenzu, abw
						
// access the Stata variables in Zulist, honoring touse stored as last column
	Zunames = tokens(Zulist)
	lenzu=cols(Zunames)-1
	v = Zunames[|1\lenzu|]
	v2 = Zunames[lenzu+1]
	st_view(uZ,.,v,v2)
	tnow=st_data(., tindexname)

// assume constant in last col of uZ if it exists
// account for eps as the first column of uZ
	if (cons) {
		nrows1=cols(uZ)-2
		nrows2=1
	}
	else {
		nrows1=cols(uZ)-1
		nrows2=0
	}
// [R] ivregress p.42: referencing Newey-West 1994 REStud 61(4):631-653
// define h indicator rowvector
	h = J(nrows1,1,1) \ J(nrows2,1,0)
	
// calc mstar per p.43
// Hannan (1971, 296) & Priestley (1981, 58) per Newey-West p. 633
// corrected per Alistair Hall msg to Brian Poi 17jul2008
//	T = rows(uZ)
//	oneT = 1/T
	expo = 2/9
	q = 1
//	cgamma = 1.4117
	cgamma = 1.1447
	if(kernel == "Parzen") { 
		expo = 4/25
		q = 2
		cgamma = 2.6614
	}
// cfb B102
	if(kernel == "Quadratic Spectral") {
		expo = 2/25
		q = 2
		cgamma = 1.3221
	}
// per Newey-West p.639, Anderson (1971), Priestley (1981) may provide
// guidance on setting expo for other kernels
//	mstar = trunc(20 *(T/100)^expo)
// use time span of data (not number of obs)
	mstar = trunc(20 *(tobs/100)^expo)

// calc uZ matrix
	u = uZ[.,1]
	Z = uZ[|1,2 \.,.|]
		
// calc f vector: (u_i Z_i) * h
	f =  (u :* Z) * h

// approach allows for gaps in time series
	sigmahat = J(mstar+1,1,0)
	for(j=0;j<=mstar;j++) {
		lsj = "L"+strofreal(j)
		tlag=st_data(., lsj+"."+tindexname)
		tmatrix = tnow, tlag
		svar=(tnow:<.):*(tlag:<.)		// multiply column vectors of 1s and 0s
		tmatrix=select(tmatrix,svar)	// to get intersection, and replace tmatrix
										// now calculate autocovariance; divide by nobs
		sigmahat[j+1] = quadcross(f[tmatrix[.,1],.], f[tmatrix[.,2],.]) / nobs
	}

// calc shat(q), shat(0)
	shatq = 0
	shat0 = sigmahat[1]
	for(j=1;j<=mstar;j++) {
		shatq = shatq + 2 * sigmahat[j+1] * j^q
		shat0 = shat0 + 2 * sigmahat[j+1]
	}

// calc gammahat
	expon = 1/(2*q+1)
	gammahat = cgamma*( (shatq/shat0)^2 )^expon
// use time span of data tobs (not number of obs T)
		m = gammahat * tobs^expon

// calc opt lag
	if(kernel == "Bartlett" | kernel == "Parzen") {
		optlag = min((trunc(m),mstar))
	}
	else if(kernel == "Quadratic Spectral") {
		optlag = min((m,mstar))
	}
		
// if optlag is the optimal lag to be used, we need to add one to 
// specify bandwidth in ivreg2 terms
	abw = optlag + 1
	st_local("abw",strofreal(abw))
} // end program s_abw

// *********** s_rmcoll2 (replacement for Stata _rmcollright etc. **********

void s_rmcoll2(	string scalar fv_vnames,
				string scalar vnames,
				string scalar wname,
				string scalar touse,
				scalar cons,
				scalar lindep)
{
	st_view(X=., ., tokens(fv_vnames), touse)
	st_view(w=., ., tokens(wname), touse)
	st_view(mtouse=., ., tokens(touse), touse)

if (cons) {
		Xmean=mean(X,w)
		XX=quadcrossdev(X,Xmean, w, X,Xmean)
	}
	else {
		XX=quadcross(X, w, X)
	}

XXinv=invsym(XX, range(1,cols(X),1))

st_numscalar("r(k_omitted)", diag0cnt(XXinv))
	if (lindep) {
		st_matrix("r(lindep)", XX*XXinv)
	}
	smat = (diagonal(XXinv) :== 0)'
	vl=tokens(vnames)
	vl_drop = select(vl, smat)
	vl_keep = select(vl, (1 :- smat))

if (cols(vl_keep)>0) {
		st_global("r(varlist)", invtokens(vl_keep))
	}
	if (cols(vl_drop)>0) {
		st_global("r(omitted)", invtokens(vl_drop))
	}
}	// end program s_rmcoll2

// ************** Add omitted Mata utility ************************

void s_AddOmitted(	string scalar bname,
					string scalar vname,
					string scalar cnumlist,
					scalar eq_ct,
					scalar rhs0_ct,
					scalar rhs1_ct)

{
	b = st_matrix(bname)
	V = st_matrix(vname)
	cn = strtoreal(tokens(cnumlist))
// cnumlist is the list of columns in the single-equation new big matrix in which
// the non-zero entries from the reduced matrix (bmat or vmat) will appear.
// E.g., if newb will be [mpg o.mpg2 _cons] then cnum = [1 3].
	col_ct = eq_ct * rhs0_ct
	
	newb = J(1,col_ct,0)
	newV = J(col_ct,col_ct,0)
	
// Code needs to accommodate multi-equation case.  Since all equations will have
// same reduced and full list of vars, in the same order, can do this with Kronecker
// products etc.  Second term below is basically the offset for each equation.
	cn = (J(1,eq_ct,1) # cn) + ((range(0,eq_ct-1,1)' # J(1,rhs1_ct,1) ) * rhs0_ct)

// Insert the values from the reduced matrices into the right places in the big matrices.
	newb[1, cn] = b
	newV[cn, cn] = V
	
	st_matrix("r(b)", newb)
	st_matrix("r(V)", newV)
	
}

// ************** Partial out *************************************

void s_partial(		string scalar yname,
					string scalar X1names,
					string scalar X2names,
					string scalar Z1names,
					string scalar Pnames,
					string scalar touse,
					string scalar weight,
					string scalar wvarname,
					scalar wf,
					scalar N,
					scalar cons)

{

// All varnames should be basic form, no FV or TS operators etc.
//  y = dep var
// X1 = endog regressors
// X2 = exog regressors = included IVs
// Z1 = excluded instruments
// Z2 = included IVs = X2
// PZ = variables to partial out
// cons = 0 or 1

ytoken=tokens(yname)
	X1tokens=tokens(X1names)
	X2tokens=tokens(X2names)
	Z1tokens=tokens(Z1names)
	Ptokens=tokens(Pnames)
	Ytokens = (ytoken, X1tokens, X2tokens, Z1tokens)

st_view(wvar, ., st_tsrevar(wvarname), touse)
	st_view(Y, ., Ytokens, touse)
	st_view(P, ., Ptokens, touse)
	L = cols(P)

if (cons & L>0) {					//  Vars to partial out including constant
		Ymeans = mean(Y,wf*wvar)
		Pmeans = mean(P,wf*wvar)
		PY = quadcrossdev(P, Pmeans, wf*wvar, Y, Ymeans)
		PP = quadcrossdev(P, Pmeans, wf*wvar, P, Pmeans)
	}
	else if (!cons & L>0) {				//  Vars to partial out NOT including constant
		PY = quadcross(P, wf*wvar, Y)
		PP = quadcross(P, wf*wvar, P)
	}
	else {								//  Only constant to partial out = demean
		Ymeans = mean(Y,wf*wvar)
	}

//	Partial-out coeffs. Default Cholesky; use QR if not full rank and collinearities present.
//	Not necessary if no vars other than constant
	if (L>0) {
		b = cholqrsolve(PP, PY)
	}
//	Replace with residuals
	if (cons & L>0) {					//  Vars to partial out including constant
		Y[.,.] = (Y :- Ymeans) - (P :- Pmeans)*b
	}
	else if (!cons & L>0) {				//  Vars to partial out NOT including constant
		Y[.,.] = Y - P*b
	}
	else {								//  Only constant to partial out = demean
		Y[.,.] = (Y :- Ymeans)
	}

} // end program s_partial

// ************** Common cross-products *************************************

void s_crossprods(	string scalar yname,
					string scalar X1names,
					string scalar X2names,
					string scalar Z1names,
					string scalar touse,
					string scalar weight,
					string scalar wvarname,
					scalar wf,
					scalar N)

{

//  y = dep var
// X1 = endog regressors
// X2 = exog regressors = included IVs
// Z1 = excluded instruments
// Z2 = included IVs = X2

ytoken=tokens(yname)
	X1tokens=tokens(X1names)
	X2tokens=tokens(X2names)
	Z1tokens=tokens(Z1names)

Xtokens = (X1tokens, X2tokens)
	Ztokens = (Z1tokens, X2tokens)

K1=cols(X1tokens)
	K2=cols(X2tokens)
	K=K1+K2
	L1=cols(Z1tokens)
	L2=cols(X2tokens)
	L=L1+L2

st_view(wvar, ., st_tsrevar(wvarname), touse)
	st_view(A, ., st_tsrevar((ytoken, Xtokens, Z1tokens)), touse)

AA = quadcross(A, wf*wvar, A)

if (K>0) {
		XX = AA[(2::K+1),(2..K+1)]
		Xy  = AA[(2::K+1),1]
	}
	if (K1>0) {
		X1X1 = AA[(2::K1+1),(2..K1+1)]
	}

if (L1 > 0) {
		Z1Z1 = AA[(K+2::rows(AA)),(K+2..rows(AA))]
	}

if (L2 > 0) {
		Z2Z2 = AA[(K1+2::K+1), (K1+2::K+1)]
		Z2y  = AA[(K1+2::K+1), 1]
	}

if ((L1>0) & (L2>0)) {
		Z2Z1 = AA[(K1+2::K+1), (K+2::rows(AA))]
		ZZ2 = Z2Z1, Z2Z2
		ZZ1 = Z1Z1, Z2Z1'
		ZZ = ZZ1 \ ZZ2
	}
	else if (L1>0) {
		ZZ = Z1Z1
	}
	else {
// L1=0
		ZZ  = Z2Z2
		ZZ2 = Z2Z2
	}

if ((K1>0) & (L1>0)) {						// K1>0, L1>0
		X1Z1 = AA[(2::K1+1), (K+2::rows(AA))]
	}

if ((K1>0) & (L2>0)) {
		X1Z2 = AA[(2::K1+1), (K1+2::K+1)]
		if (L1>0) {								// K1>0, L1>0, L2>0
			X1Z = X1Z1, X1Z2
			XZ = X1Z \ ZZ2
		}
		else {									// K1>0, L1=0, L2>0
			XZ = X1Z2 \ ZZ2
			X1Z = X1Z2
		}
	}
	else if (K1>0) {							// K1>0, L2=0
		XZ = X1Z1
		X1Z= X1Z1
	}
	else if (L1>0) {							// K1=0, L2>0
		XZ = AA[(2::K+1),(K+2..rows(AA))], AA[(2::K+1),(2..K+1)]
	}
	else {										// K1=0, L2=0
		XZ = ZZ
	}

if ((L1>0) & (L2>0)) {
		Zy = AA[(K+2::rows(AA)), 1] \ AA[(K1+2::K+1), 1]
		ZY = AA[(K+2::rows(AA)), (1..K1+1)] \ AA[(K1+2::K+1), (1..K1+1)]
		Z2Y = AA[(K1+2::K+1), (1..K1+1)]
	}
	else if (L1>0) {
		Zy = AA[(K+2::rows(AA)), 1]
		ZY = AA[(K+2::rows(AA)), (1..K1+1)]
	}
	else if (L2>0) {
		Zy = AA[(K1+2::K+1), 1]
		ZY = AA[(K1+2::K+1), (1..K1+1)]
		Z2Y = ZY
	}
// Zy, ZY, Z2Y not created if L1=L2=0

YY  = AA[(1::K1+1), (1..K1+1)]
	yy  = AA[1,1]
	st_subview(y, A, ., 1)
	ym    = sum(wf*wvar:*y)/N
	yyc   = quadcrossdev(y, ym, wf*wvar, y, ym)

XXinv = invsym(XX)
	if (Xtokens==Ztokens) {
		ZZinv = XXinv
		XPZXinv = XXinv
	}
	else {
		ZZinv = invsym(ZZ)
		XPZX  = makesymmetric(XZ*ZZinv*XZ')
		XPZXinv=invsym(XPZX)
	}

// condition numbers
	condxx=cond(XX)
	condzz=cond(ZZ)

st_matrix("r(XX)", XX)
	st_matrix("r(X1X1)", X1X1)
	st_matrix("r(X1Z)", X1Z)
	st_matrix("r(ZZ)", ZZ)
	st_matrix("r(Z2Z2)", Z2Z2)
	st_matrix("r(Z1Z2)", Z2Z1')
	st_matrix("r(Z2y)",Z2y)
	st_matrix("r(XZ)", XZ)
	st_matrix("r(Xy)", Xy)
	st_matrix("r(Zy)", Zy)
	st_numscalar("r(yy)", yy)
	st_numscalar("r(yyc)", yyc)
	st_matrix("r(YY)", YY)
	st_matrix("r(ZY)", ZY)
	st_matrix("r(Z2Y)", Z2Y)
	st_matrix("r(XXinv)", XXinv)
	st_matrix("r(ZZinv)", ZZinv)
	st_matrix("r(XPZXinv)", XPZXinv)
	st_numscalar("r(condxx)",condxx)
	st_numscalar("r(condzz)",condzz)

} // end program s_crossprods

// *************** 1st step GMM ******************** //
// Can be either efficient or inefficient.
// Can be IV or other 1-step GMM estimator.

void s_gmm1s(	string scalar ZZmatrix,
				string scalar XXmatrix,
				string scalar XZmatrix,
				string scalar Zymatrix,
				string scalar ZZinvmatrix,
				string scalar yname,
				string scalar ename,
				string scalar Xnames, 
				string scalar Znames, 
				string scalar touse,
				string scalar weight,
				string scalar wvarname,
				scalar wf,
				scalar N,
				string scalar Wmatrix,
				string scalar Smatrix,
				scalar dofminus,
				scalar efficient,			//  flag to indicate that 1st-step GMM is efficient
				scalar overid,				//  not guaranteed to be right if nocollin option used!
				scalar useqr)				//  flag to force use of QR instead of Cholesky solver
{

Ztokens=tokens(Znames)
	Xtokens=tokens(Xnames)

st_view(Z,    ., st_tsrevar(Ztokens),  touse)
	st_view(X,    ., st_tsrevar(Xtokens),  touse)
	st_view(y,    ., st_tsrevar(yname),  touse)
	st_view(e,    ., ename,  touse)
	st_view(wvar, ., st_tsrevar(wvarname), touse)

// Our convention is that regressors are [endog   included exog]
// and instruments are                   [excluded exog  included exog]
// If a constant is included, it is the last column.

ZZ = st_matrix(ZZmatrix)
	XX = st_matrix(XXmatrix)
	XZ = st_matrix(XZmatrix)
	Zy = st_matrix(Zymatrix)
	ZZinv = st_matrix(ZZinvmatrix)

QZZ    = ZZ   / N
	QXX    = XX   / N
	QXZ    = XZ   / N
	QZy    = Zy   / N
	QZZinv = ZZinv*N
	
	useqr = (diag0cnt(QZZinv)>0) | useqr

// Weighting matrix supplied (and inefficient GMM)
	if (Wmatrix~="") {
		W = st_matrix(Wmatrix)
		useqr = (diag0cnt(W)>0) | useqr
	}
// Var-cov matrix of orthog conditions supplied
	if (Smatrix~="") {
		omega=st_matrix(Smatrix)
		useqr = (diag0cnt(omega)>0) | useqr
	}

if (efficient) {									//  Efficient 1-step GMM block: OLS, IV or provided S
		if ((Xtokens==Ztokens) & (Smatrix=="")) {		//  OLS

beta = cholqrsolve(QZZ, QZy, useqr)
			beta = beta'
			e[.,.] = y - X * beta'						//  update residuals
			ee = quadcross(e, wf*wvar, e)
			sigmasq=ee/(N-dofminus)
			omega = sigmasq * QZZ
			W = 1/sigmasq * QZZinv
			V = 1/N * sigmasq * QZZinv
			rankS = rows(omega) - diag0cnt(QZZinv)		//  inv(omega) is proportional to inv(QZZ)
			rankV = rows(V) - diag0cnt(V)				//  inv(V) is proportional to inv(QZZ)
		}
		else if (Smatrix=="") {							//  IV
			aux1 = cholqrsolve(QZZ, QXZ', useqr)
			aux2 = cholqrsolve(QZZ, QZy, useqr)
			aux3 = makesymmetric(QXZ * aux1)
			beta = cholqrsolve(aux3, QXZ * aux2, useqr)
			beta = beta'
			e[.,.] = y - X * beta'						//  update residuals
			ee = quadcross(e, wf*wvar, e)
			sigmasq = ee/(N-dofminus)
			omega = sigmasq * QZZ
			W = 1/sigmasq * QZZinv
			V = 1/N * sigmasq * invsym(aux3)
			rankS = rows(omega) - diag0cnt(QZZinv)		//  inv(omega) is proportional to inv(QZZ)
			rankV = rows(V) - diag0cnt(V)				//  V is proportional to inv(aux3)
		}
		else {											//  efficient GMM with provided S (=omega)
			aux1 = cholqrsolve(omega, QXZ', useqr)
			aux2 = cholqrsolve(omega, QZy, useqr)
			aux3 = makesymmetric(QXZ * aux1)
			beta = cholqrsolve(aux3, QXZ * aux2, useqr)
			beta = beta'
			e[.,.] = y - X * beta'						//  update residuals
			ee = quadcross(e, wf*wvar, e)
			sigmasq=ee/(N-dofminus)
			W = invsym(omega)
			V = 1/N * invsym(aux3)						//  Normalize by N
			rankS = rows(omega) - diag0cnt(W)			//  since W=inv(omega)
			rankV = rows(V) - diag0cnt(V)				//  since V is prop to inv(aux3)
		}
		if (overid) {									// J if overidentified
			Ze = quadcross(Z, wf*wvar, e)
			gbar = Ze / N
			aux4 = cholqrsolve(omega, gbar, useqr)
			j = N * gbar' * aux4
		}
		else {
			j=0
		}
		st_matrix("r(beta)", beta)
		st_matrix("r(V)", V)
		st_matrix("r(S)", omega)
		st_matrix("r(W)", W)
		st_numscalar("r(rss)", ee)
		st_numscalar("r(j)", j)
		st_numscalar("r(sigmasq)", sigmasq)
		st_numscalar("r(rankS)", rankS)
		st_numscalar("r(rankV)", rankV)
	}
	else {												//  inefficient 1st-step GMM; don't need V, S, j etc.
		if ((Xtokens==Ztokens) & (Wmatrix=="")) {		//  OLS
			beta = cholqrsolve(QZZ, QZy, useqr)
			beta = beta'
			e[.,.] = y - X * beta'						//  update residuals
			ee = quadcross(e, wf*wvar, e)
			sigmasq=ee/(N-dofminus)
			W = 1/sigmasq * QZZinv
			QXZ_W_QZX = 1/sigmasq * QZZ					//  b/c W incorporates sigma^2
		}
		else if (Wmatrix=="") {							//  IV
			aux1 = cholqrsolve(QZZ, QXZ', useqr)
			aux2 = cholqrsolve(QZZ, QZy, useqr)
			aux3 = makesymmetric(QXZ * aux1)
			beta = cholqrsolve(aux3, QXZ * aux2, useqr)
			beta = beta'
			e[.,.] = y - X * beta'						//  update residuals
			ee = quadcross(e, wf*wvar, e)
			sigmasq=ee/(N-dofminus)
			W = 1/sigmasq * QZZinv
			QXZ_W_QZX = 1/sigmasq * aux3				//  b/c IV weighting matrix incorporates sigma^2
		}
		else {											//  some other 1st step inefficient GMM with provided W
			QXZ_W_QZX = QXZ * W * QXZ'
			_makesymmetric(QXZ_W_QZX)
			beta = cholqrsolve(QXZ_W_QZX, QXZ * W * QZy, useqr)
			beta = beta'
			e[.,.] = y - X * beta'						//  update residuals
		}
		st_matrix("r(QXZ_W_QZX)", QXZ_W_QZX)
		st_matrix("r(beta)", beta)
		st_matrix("r(W)",W)								//  always return W
	}

} // end program s_gmm1s

// *************** efficient GMM ******************** //
// Uses inverse of provided S matrix as weighting matrix.
// IV won't be done here but code would work for it as a special case.

void s_egmm(	string scalar ZZmatrix,
				string scalar XXmatrix,
				string scalar XZmatrix,
				string scalar Zymatrix,
				string scalar ZZinvmatrix,
				string scalar yname,
				string scalar ename,
				string scalar Xnames, 
				string scalar Znames, 
				string scalar touse,
				string scalar weight,
				string scalar wvarname,
				scalar wf,
				scalar N,
				string scalar Smatrix,					//  always provided
				scalar dofminus,
				scalar overid,							//  not guaranteed to be right if -nocollin- used!
				scalar useqr)
{

Ztokens=tokens(Znames)
	Xtokens=tokens(Xnames)

// Our convention is that regressors are [endog   included exog]
// and instruments are                   [excluded exog  included exog]
// If a constant is included, it is the last column.

ZZ			= st_matrix(ZZmatrix)
	XX			= st_matrix(XXmatrix)
	XZ			= st_matrix(XZmatrix)
	Zy			= st_matrix(Zymatrix)
	ZZinv		= st_matrix(ZZinvmatrix)

QZZ    = ZZ   / N
	QXX    = XX   / N
	QXZ    = XZ   / N
	QZy    = Zy   / N
	QZZinv = ZZinv*N

// Var-cov matrix of orthog conditions supplied
	omega=st_matrix(Smatrix)
	W = invsym(omega)							//  Efficient GMM weighting matrix
	rankS = rows(omega) - diag0cnt(W)			//  since W=inv(omega)
	
	if (rankS<rows(omega)) {					//  omega not full rank; W=inv(omega) dubious, exit with error
errprintf("\nError: estimated covariance matrix of moment conditions not of full rank,")
errprintf("\n       and optimal GMM weighting matrix not unique.")
errprintf("\nPossible causes:")
errprintf("\n       singleton dummy variable (dummy with one 1 and N-1 0s or vice versa)")
errprintf("\n       {help ivreg2##partial:partial} option may address problem.\n")
		exit(506)
	}

aux1 = cholqrsolve(omega, QXZ', useqr)
	aux2 = cholqrsolve(omega, QZy, useqr)
	aux3 = makesymmetric(QXZ * aux1)
	beta = cholqrsolve(aux3, QXZ * aux2, useqr)
	beta = beta'

// The GMM estimator is "root-N consistent", and technically we do
// inference on sqrt(N)*beta.  By convention we work with beta,
// so we adjust the var-cov matrix instead:
	V = 1/N * invsym(aux3)
	rankV = rows(V) - diag0cnt(V)		//  since V is proportional to inv(aux3)

//	above equivalent to but more accurate than:
//	W = invsym(omega) or W = QZZinv
//	QXZ_W_QZX = QXZ * W * QXZ'
//	QXZ_W_QZXinv=invsym(makesymmetric(QXZ_W_QZX))
//	beta = (QXZ_W_QZXinv * QXZ * W * QZy)
//	V = QXZ_W_QZXinv

e[.,.] = y - X * beta'
	ee = quadcross(e, wf*wvar, e)
	sigmasq=ee/(N-dofminus)

// J if overidentified
	if (cols(Z) > cols(X)) {
		Ze = quadcross(Z, wf*wvar, e)
		gbar = Ze / N
		aux4 = cholqrsolve(omega, gbar, useqr)
		j = N * gbar' * aux4
	}
	else {
		j=0
	}

st_matrix("r(beta)", beta)
	st_matrix("r(V)", V)
	st_matrix("r(W)", W)
	st_numscalar("r(rss)", ee)
	st_numscalar("r(j)", j)
	st_numscalar("r(sigmasq)", sigmasq)
	st_numscalar("r(rankV)",rankV)
	st_numscalar("r(rankS)",rankS)

} // end program s_egmm

// *************** inefficient GMM ******************** //

void s_iegmm(	string scalar ZZmatrix,
				string scalar XXmatrix,
				string scalar XZmatrix,
				string scalar Zymatrix,
				string scalar QXZ_W_QZXmatrix,
				string scalar yname,
				string scalar ename,
				string scalar Xnames, 
				string scalar Znames, 
				string scalar touse,
				string scalar weight,
				string scalar wvarname,
				scalar wf,
				scalar N,
				string scalar Wmatrix,
				string scalar Smatrix,
				string scalar bname,
				scalar dofminus,
				scalar overid,
				scalar useqr)
{

Ztokens=tokens(Znames)
	Xtokens=tokens(Xnames)

// Our convention is that regressors are [endog   included exog]
// and instruments are                   [excluded exog  included exog]
// If a constant is included, it is the last column.

QZZ			= st_matrix(ZZmatrix) / N
	QXX			= st_matrix(XXmatrix) / N
	QXZ			= st_matrix(XZmatrix) / N
	QZy			= st_matrix(Zymatrix) / N
	QXZ_W_QZX	= st_matrix(QXZ_W_QZXmatrix)

useqr = (diag0cnt(QXZ_W_QZX)>0) | useqr

// beta is supplied
	beta = st_matrix(bname)

// Inefficient weighting matrix supplied
	W = st_matrix(Wmatrix)

// Var-cov matrix of orthog conditions supplied
	omega=st_matrix(Smatrix)

// Residuals are supplied
	ee = quadcross(e, wf*wvar, e)
	sigmasq=ee/(N-dofminus)

// Calculate V and J.

// V
// The GMM estimator is "root-N consistent", and technically we do
// inference on sqrt(N)*beta.  By convention we work with beta, so we adjust
// the var-cov matrix instead:
	aux5 = cholqrsolve(QXZ_W_QZX, QXZ * W, useqr)
	V = 1/N * aux5 * omega * aux5'
	_makesymmetric(V)

// alternative
//	QXZ_W_QZXinv=invsym(QXZ_W_QZX)
//	V = 1/N * QXZ_W_QZXinv * QXZ * W * omega * W * QXZ' * QXZ_W_QZXinv

rankV = rows(V) - diag0cnt(invsym(V))				//  need explicitly to calc rank
	rankS = rows(omega) - diag0cnt(invsym(omega))		//  need explicitly to calc rank

// J if overidentified
	if (overid) {
// Note that J requires efficient GMM residuals, which means do 2-step GMM to get them.
//		QXZ_W2s_QZX = QXZ * W2s * QXZ'
//		_makesymmetric(QXZ_W2s_QZX)
//		QXZ_W2s_QZXinv=invsym(QXZ_W2s_QZX)
//		beta2s = (QXZ_W2s_QZXinv * QXZ * W2s * QZy)
		aux1 = cholqrsolve(omega, QXZ', useqr)
		aux2 = cholqrsolve(omega, QZy, useqr)
		aux3s = makesymmetric(QXZ * aux1)
		beta2s = cholqrsolve(aux3s, QXZ * aux2, useqr)
		beta2s = beta2s'
		e2s = y - X * beta2s'
		Ze2s = quadcross(Z, wf*wvar, e2s)
		gbar = Ze2s / N
		aux4 = cholqrsolve(omega, gbar, useqr)
		j = N * gbar' * aux4
	}
	else {
		j=0
	}

st_matrix("r(V)", V)
	st_numscalar("r(j)", j)
	st_numscalar("r(rss)", ee)
	st_numscalar("r(sigmasq)", sigmasq)
	st_numscalar("r(rankV)",rankV)
	st_numscalar("r(rankS)",rankS)

} // end program s_iegmm

// *************** LIML ******************** //

void s_liml(	string scalar ZZmatrix,
				string scalar XXmatrix,
				string scalar XZmatrix,
				string scalar Zymatrix,
				string scalar Z2Z2matrix,
				string scalar YYmatrix,
				string scalar ZYmatrix,
				string scalar Z2Ymatrix,
				string scalar Xymatrix,
				string scalar ZZinvmatrix,
				string scalar yname,
				string scalar Ynames,
				string scalar ename,
				string scalar Xnames,
				string scalar X1names,
				string scalar Znames,
				string scalar Z1names,
				string scalar Z2names,
				scalar fuller,
				scalar kclass,
				string scalar coviv,
				string scalar touse,
				string scalar weight,
				string scalar wvarname,
				scalar wf,
				scalar N,
				string scalar robust,
				string scalar clustvarname,
				string scalar clustvarname2,
				string scalar clustvarname3,
				scalar bw,
				string scalar kernel,
				string scalar sw,
				string scalar psd,
				string scalar ivarname,
				string scalar tvarname,
				string scalar tindexname,
				scalar tdelta,
				scalar center,
				scalar dofminus,
				scalar useqr)

{
	struct ms_vcvorthog scalar vcvo

vcvo.ename			= ename
	vcvo.Znames			= Znames
	vcvo.touse			= touse
	vcvo.weight			= weight
	vcvo.wvarname		= wvarname
	vcvo.robust			= robust
	vcvo.clustvarname	= clustvarname
	vcvo.clustvarname2	= clustvarname2
	vcvo.clustvarname3	= clustvarname3
	vcvo.kernel			= kernel
	vcvo.sw				= sw
	vcvo.psd			= psd
	vcvo.ivarname		= ivarname
	vcvo.tvarname		= tvarname
	vcvo.tindexname		= tindexname
	vcvo.wf				= wf
	vcvo.N				= N
	vcvo.bw				= bw
	vcvo.tdelta			= tdelta
	vcvo.center			= center
	vcvo.dofminus		= dofminus
	vcvo.ZZ				= st_matrix(ZZmatrix)

// X1 = endog regressors
// X2 = exog regressors = included IVs
// Z1 = excluded instruments
// Z2 = included IVs = X2

Ytokens=tokens(Ynames)
	Ztokens=tokens(Znames)
	Z1tokens=tokens(Z1names)
	Z2tokens=tokens(Z2names)
	Xtokens=tokens(Xnames)
	X1tokens=tokens(X1names)

st_view(Z,     ., st_tsrevar(Ztokens),   touse)
	st_view(X,     ., st_tsrevar(Xtokens),   touse)
	st_view(y,     ., st_tsrevar(yname),     touse)
	st_view(e,     ., ename,                 touse)
	st_view(wvar,  ., st_tsrevar(wvarname),  touse)

vcvo.e		= &e
	vcvo.Z		= &Z
	vcvo.wvar	= &wvar

// Our convention is that regressors are [endog   included exog]
// and instruments are                   [excluded exog  included exog]
// If a constant is included, it is the last column.

QZZ		= st_matrix(ZZmatrix) / N
	QXX		= st_matrix(XXmatrix) / N
	QXZ		= st_matrix(XZmatrix) / N
	QZy		= st_matrix(Zymatrix) / N
	QZ2Z2	= st_matrix(Z2Z2matrix) / N
	QYY		= st_matrix(YYmatrix) / N
	QZY		= st_matrix(ZYmatrix) / N
	QZ2Y	= st_matrix(Z2Ymatrix) / N
	QXy		= st_matrix(Xymatrix) / N
	QZZinv	= st_matrix(ZZinvmatrix)*N

useqr = (diag0cnt(QZZ)>0) | useqr

// kclass=0 => LIML or Fuller LIML so calculate lambda
	if (kclass == 0) {
		aux1 = cholqrsolve(QZZ, QZY, useqr)
		QWW = QYY - QZY'*aux1
		_makesymmetric(QWW)
		if (cols(Z2tokens) > 0) {
			aux2 = cholqrsolve(QZ2Z2, QZ2Y, useqr)
			QWW1 = QYY - QZ2Y'*aux2
			_makesymmetric(QWW1)
		}
		else {
// Special case of no exogenous regressors
			QWW1 = QYY
		}
		M=matpowersym(QWW, -0.5)
		Eval=symeigenvalues(M*QWW1*M)
		lambda=rowmin(Eval)
	}

// Exactly identified but might not be exactly 1, so make it so
	if (cols(Z)==cols(X)) {
		lambda=1
	}

if (fuller > (N-cols(Z))) {
printf("\n{error:Error: invalid choice of Fuller LIML parameter.}\n")
		exit(error(3351))
	}
	else if (fuller > 0) {
		k = lambda - fuller/(N-cols(Z))
	}
	else if (kclass > 0) {
		k = kclass
	}
	else {
		k = lambda
	}

aux3 = cholqrsolve(QZZ, QXZ', useqr)
	QXhXh=(1-k)*QXX + k*QXZ*aux3
	_makesymmetric(QXhXh)
	aux4 = cholqrsolve(QZZ, QZy, useqr)
	aux5 = cholqrsolve(QXhXh, QXZ, useqr)
	aux6 = cholqrsolve(QXhXh, QXy, useqr)
	beta = aux6*(1-k) + k*aux5*aux4
	beta = beta'

e[.,.] = y - X * beta'
	ee = quadcross(e, wf*wvar, e)
	sigmasq = ee /(N-dofminus)

omega = m_omega(vcvo)

QXhXhinv=invsym(QXhXh)

if ((robust=="") & (clustvarname=="") & (kernel=="")) {
// Efficient LIML
		if (coviv=="") {
// Note dof correction is already in sigmasq
			V = 1/N * sigmasq * QXhXhinv
			rankV = rows(V) - diag0cnt(V)		//  since V is proportional to inv(QXhXh)
		}
		else {
			aux7 = makesymmetric(QXZ * aux3)
			V = 1/N * sigmasq * invsym(aux7)
			rankV = rows(V) - diag0cnt(V)		//  since V is proportional to inv(aux7)
		}
		rankS = rows(omega) - diag0cnt(invsym(omega))
		if (cols(Z)>cols(X)) {
			Ze = quadcross(Z, wf*wvar, e)
			gbar = Ze / N
			aux8 = cholqrsolve(omega, gbar, useqr)
			j = N * gbar' * aux8
		}
		else {
			j=0
		}
	}
	else {
// Inefficient LIML
		if (coviv=="") {
			aux9 = cholqrsolve(QZZ, aux5', useqr)
			V = 1/N * aux9' * omega * aux9
			_makesymmetric(V)
			rankV = rows(V) - diag0cnt(invsym(V))				// need explicitly to calc rank
			rankS = rows(omega) - diag0cnt(invsym(omega))		// need explicitly to calc rank
		}
		else {
			aux10 = QXZ * aux3
			_makesymmetric(aux10)
			aux11 = cholqrsolve(aux10, aux3', useqr)
			V = 1/N * aux11 * omega * aux11'
			_makesymmetric(V)
			rankV = rows(V) - diag0cnt(invsym(V))				// need explicitly to calc rank
			rankS = rows(omega) - diag0cnt(invsym(omega))		// need explicitly to calc rank
		}
		if (cols(Z)>cols(X)) {
			aux12 = cholqrsolve(omega, QXZ', useqr)
			aux13 = cholqrsolve(omega, QZy, useqr)
			aux14 = makesymmetric(QXZ * aux12)
			beta2s = cholqrsolve(aux14, QXZ * aux13, useqr)
			beta2s = beta2s'
			e2s = y - X * beta2s'
			Ze2s = quadcross(Z, wf*wvar, e2s)
			gbar = Ze2s / N
			aux15 = cholqrsolve(omega, gbar, useqr)
			j = N * gbar' * aux15
		}
		else {
			j=0
		}
	}
	_makesymmetric(V)

st_matrix("r(beta)", beta)
	st_matrix("r(S)", omega)
	st_matrix("r(V)", V)
	st_numscalar("r(lambda)", lambda)
	st_numscalar("r(kclass)", k)
	st_numscalar("r(j)", j)
	st_numscalar("r(rss)", ee)
	st_numscalar("r(sigmasq)", sigmasq)
	st_numscalar("r(rankV)",rankV)
	st_numscalar("r(rankS)",rankS)

} // end program s_liml

// *************** CUE ******************** //

void s_gmmcue(	string scalar ZZmatrix,
				string scalar XZmatrix,
				string scalar yname,
				string scalar ename,
				string scalar Xnames, 
				string scalar Znames, 
				string scalar touse,
				string scalar weight,
				string scalar wvarname,
				scalar wf,
				scalar N,
				string scalar robust,
				string scalar clustvarname,
				string scalar clustvarname2,
				string scalar clustvarname3,
				scalar bw,
				string scalar kernel,
				string scalar sw,
				string scalar psd,
				string scalar ivarname,
				string scalar tvarname,
				string scalar tindexname,
				scalar tdelta,
				string scalar bname,
				string scalar b0name,
				scalar center,
				scalar dofminus,
				scalar useqr)

{

struct ms_vcvorthog scalar vcvo

Ztokens=tokens(Znames)
	Xtokens=tokens(Xnames)

st_view(Z,    ., st_tsrevar(Ztokens),  touse)
	st_view(X,    ., st_tsrevar(Xtokens),  touse)
	st_view(y,    ., st_tsrevar(yname),    touse)
	st_view(e,    ., ename,                touse)
	st_view(wvar, ., st_tsrevar(wvarname), touse)

// Pointers to views
	vcvo.e		= &e
	vcvo.Z		= &Z
	vcvo.wvar	= &wvar
	py			= &y
	pX			= &X

if (b0name=="") {

// CUE beta not supplied, so calculate/optimize

// Our convention is that regressors are [endog   included exog]
// and instruments are                   [excluded exog  included exog]
// If a constant is included, it is the last column.

// CUE is preceded by IV or 2-step GMM to get starting values.
// Stata convention is that parameter vectors are row vectors, and optimizers
// require this, so must conform to this in what follows.

beta_init = st_matrix(bname)

// What follows is how to set out an optimization in Stata.  First, initialize
// the optimization structure in the variable S.  Then tell Mata where the
// objective function is, that it's a minimization, that it's a "d0" type of
// objective function (no analytical derivatives or Hessians), and that the
// initial values for the parameter vector are in beta_init.  Finally, optimize.
		S = optimize_init()

optimize_init_evaluator(S, &m_cuecrit())
		optimize_init_which(S, "min")
		optimize_init_evaluatortype(S, "d0")
		optimize_init_params(S, beta_init)
// CUE objective function takes 3 extra arguments: y, X and the structure with omega details
		optimize_init_argument(S, 1, py)
		optimize_init_argument(S, 2, pX)
		optimize_init_argument(S, 3, vcvo)
		optimize_init_argument(S, 4, useqr)

beta = optimize(S)

// The last evaluation of the GMM objective function is J.
		j = optimize_result_value(S)

// Call m_omega one last time to get CUE weighting matrix.
		e[.,.] = y - X * beta'
		omega = m_omega(vcvo)
	}
	else {
// CUE beta supplied, so obtain maximized GMM obj function at b0
		beta = st_matrix(b0name)
		e[.,.] = y - X * beta'
		omega = m_omega(vcvo)
//		W = invsym(omega)
		gbar = 1/N * quadcross(Z, wf*wvar, e)
		j = N * gbar' * cholsolve(omega, gbar, useqr)
//		j = N * gbar' * W * gbar
	}

// Bits and pieces
	QXZ = st_matrix(XZmatrix)/N

ee = quadcross(e, wf*wvar, e)
	sigmasq=ee/(N-dofminus)

//	QXZ_W_QZX = QXZ * W * QXZ'
//	_makesymmetric(QXZ_W_QZX)
//	QXZ_W_QZXinv=invsym(QXZ_W_QZX)
//	V = 1/N * QXZ_W_QZXinv
	aux1 = cholsolve(omega, QXZ')
	if (aux1[1,1]==.) {					//  omega not full rank; W=inv(omega) dubious, exit with error
errprintf("\nError: estimated covariance matrix of moment conditions not of full rank,")
errprintf("\n       and optimal GMM weighting matrix not unique.")
errprintf("\nPossible causes:")
errprintf("\n       collinearities in instruments (if -nocollin- option was used)")
errprintf("\n       singleton dummy variable (dummy with one 1 and N-1 0s or vice versa)")
errprintf("\n       {help ivreg2##partial:partial} option may address problem.\n")
		exit(506)
	}
	aux3 = makesymmetric(QXZ * aux1)
	V = 1/N * invsym(aux3)
	if (diag0cnt(V)) {					//  V not full rank, likely caused by collinearities;
										//  b dubious, exit with error
errprintf("\nError: estimated variance matrix of b not of full rank, and CUE estimates")
errprintf("\n       unreliable; may be caused by collinearities\n")
		exit(506)
	}
	W = invsym(omega)
	
	st_matrix("r(beta)", beta)
	st_matrix("r(S)", omega)
	st_matrix("r(W)", W)
	st_matrix("r(V)", V)
	st_numscalar("r(j)", j)
	st_numscalar("r(rss)", ee)
	st_numscalar("r(sigmasq)", sigmasq)

} // end program s_gmmcue

// CUE evaluator function.
// Handles only d0-type optimization; todo, g and H are just ignored.
// beta is the parameter set over which we optimize, and 
// J is the objective function to minimize.

void m_cuecrit(todo, beta, pointer py, pointer pX, struct ms_vcvorthog scalar vcvo, useqr, j, g, H)
{
	*vcvo.e[.,.] = *py - *pX * beta'

omega = m_omega(vcvo)

// Calculate gbar=Z'*e/N
	gbar = 1/vcvo.N * quadcross(*vcvo.Z, vcvo.wf*(*vcvo.wvar), *vcvo.e)
	aux1 = cholqrsolve(omega, gbar, useqr)
	j = vcvo.N * gbar' * aux1

// old method
//	W = invsym(omega)
//	j = vcvo.N * gbar' * W * gbar

} // end program CUE criterion function

// ************** ffirst-stage stats *************************************

void s_ffirst(	string scalar ZZmatrix,
				string scalar XXmatrix,
				string scalar XZmatrix,
				string scalar ZYmatrix,
				string scalar ZZinvmatrix,
				string scalar XXinvmatrix,
				string scalar XPZXinvmatrix,
				string scalar X2X2matrix,
				string scalar Z1X2matrix,
				string scalar X2ymatrix,
				string scalar ename,				//  Nx1
				string scalar ematnames,			//  Nx(K1+1)
				string scalar yname,
				string scalar X1names,
				string scalar X2names,
				string scalar Z1names, 
				string scalar touse,
				string scalar weight,
				string scalar wvarname,
				scalar wf,
				scalar N,
				scalar N_clust,
				string scalar robust,
				string scalar clustvarname,
				string scalar clustvarname2,
				string scalar clustvarname3,
				scalar bw,
				string scalar kernel,
				string scalar sw,
				string scalar psd,
				string scalar ivarname,
				string scalar tvarname,
				string scalar tindexname,
				scalar tdelta,
				scalar center,
				scalar dofminus,
				scalar sdofminus)

{

struct ms_vcvorthog scalar vcvo

vcvo.Znames			= Znames
	vcvo.touse			= touse
	vcvo.weight			= weight
	vcvo.wvarname		= wvarname
	vcvo.robust			= robust
	vcvo.clustvarname	= clustvarname
	vcvo.clustvarname2	= clustvarname2
	vcvo.clustvarname3	= clustvarname3
	vcvo.kernel			= kernel
	vcvo.sw				= sw
	vcvo.psd			= psd
	vcvo.ivarname		= ivarname
	vcvo.tvarname		= tvarname
	vcvo.tindexname		= tindexname
	vcvo.wf				= wf
	vcvo.N				= N
	vcvo.bw				= bw
	vcvo.tdelta			= tdelta
	vcvo.center			= center
	vcvo.dofminus		= dofminus
	vcvo.ZZ				= st_matrix(ZZmatrix)

// X1 = endog regressors
// X2 = exog regressors = included IVs
// Z1 = excluded instruments
// Z2 = included IVs = X2

Xnames = invtokens( (X1names, X2names), " ")
	Znames = invtokens( (Z1names, X2names), " ")

st_view(y,     ., st_tsrevar(tokens(yname)),      touse)
	st_view(X1,    ., st_tsrevar(tokens(X1names)),    touse)
	st_view(Z1,    ., st_tsrevar(tokens(Z1names)),    touse)
	st_view(X,     ., st_tsrevar(tokens(Xnames)),     touse)
	st_view(Z,     ., st_tsrevar(tokens(Znames)),     touse)
	st_view(e,     ., ename,                          touse)
	st_view(emat,  ., tokens(ematnames),              touse)
	st_view(wvar,  ., st_tsrevar(wvarname),           touse)

vcvo.wvar	= &wvar
	vcvo.Z		= &Z
	vcvo.Znames	= Znames
	vcvo.ZZ		= st_matrix(ZZmatrix)

if ("X2names"~="") {
		st_view(X2,    ., st_tsrevar(tokens(X2names)),  touse)
	}

K1=cols(X1)
	K2=cols(X2)
	K=K1+K2
	L1=cols(Z1)
	L2=cols(X2)
	L=L1+L2
	df = L1
	df_r = N-L

ZZinv	= st_matrix(ZZinvmatrix)
	XXinv	= st_matrix(XXinvmatrix)
	XPZXinv	= st_matrix(XPZXinvmatrix)
	QZZ		= st_matrix(ZZmatrix)  / N
	QXX		= st_matrix(XXmatrix)  / N
	QZX		= st_matrix(XZmatrix)' / N
	QZY		= st_matrix(ZYmatrix)  / N
	QZZinv	= ZZinv*N
	QXXinv	= XXinv*N
	QX2X2	= st_matrix(X2X2matrix) / N
	QZ1X2	= st_matrix(Z1X2matrix) / N
	QX2y	= st_matrix(X2ymatrix)  / N

sheaall = (diagonal(XXinv) :/ diagonal(XPZXinv))  // (X1, X2) in column vector
	sheaall = (sheaall[(1::K1), 1 ])'  // Just X1 in row vector

// Full system of reduced form (col 1) and first-stage regressions
	bz = cholsolve(QZZ, QZY)
	Yhat = Z*bz
	Xhat = Yhat[.,(2..(K1+1))], X2
// VCV for full system
	eall = (y, X1) - Yhat
	ee = quadcross(eall, wf*wvar, eall)
// sigmas have large-sample dofminus correction incorporated but no small dof corrections
	sigmasqall = ee / (N-dofminus)
// rmses have small dof corrections
	rmseall = sqrt( ee / (N-L-dofminus-sdofminus) )
// V has all the classical VCVs in block diagonals
	V = sigmasqall # ZZinv
// For Wald test of excluded instruments
	R = I(L1) , J(L1, L2, 0)
// For AP and SW stats
	QXhXh = quadcross(Xhat, wf*wvar, Xhat) / N
	QXhX1 = quadcross(Xhat, wf*wvar, X1 ) / N

//	VCV for system of first-stage eqns
//  Can be robust; even if not, has correct off-block-diagonal covariances
	vcvo.ename	= ematnames
	vcvo.e		= &emat
	emat[.,.]	= eall
	omegar		= m_omega(vcvo)
	Vr			= makesymmetric(I(K1+1)#QZZinv * omegar * I(K1+1)#QZZinv) / N
	
// AR statistics from RF (first column)
	Rb			= bz[ (1::L1), 1 ]
	RVR			= Vr[| 1,1 \ L1, L1 |]
	ARWald		= Rb' * cholsolve(RVR, Rb)
	ARF			= ARWald
	ARdf		= L1
	if (clustvarname=="") {
		ARdf2	= (N-dofminus-L-sdofminus)
		ARF		= ARWald / (N-dofminus) * ARdf2 / ARdf
	}
	else {
		ARdf2	= N_clust - 1
		ARF		= ARWald / (N-1) * (N-L-sdofminus) /(N_clust) * ARdf2 / ARdf
	}
	ARFp		= Ftail(ARdf, ARdf2, ARF)
	ARchi2		= ARWald
	ARchi2p		= chi2tail(ARdf, ARchi2)

// Stock-Wright LM S statistic
// Equivalent to J stat for model with coeff on endog=0 and with inexog partialled out
// = LM version of AR stat (matches weakiv)
	if (K2>0) {
		by		= cholsolve(QX2X2, QX2y)
		e[.,.]	= y-X2*by
	}
	else {
		e[.,.]	= y
	}
// initialize residual for VCV calc to be single Nx1 vector 
	vcvo.e		= &e
	vcvo.ename	= ename
// get VCV and sstat=J
	omega		= m_omega(vcvo)
	gbar		= 1/N * quadcross(Z, wf*wvar, e)
	sstat		= N * gbar' * cholsolve(omega, gbar)
	sstatdf		= L1
	sstatp		= chi2tail(sstatdf, sstat)

// Prepare to loop over X1s for F, SW and AP stats
// initialize matrix to save first-stage results
	firstmat=J(21,0,0)
// initialize residual for VCV calc to be single Nx1 vector 
	vcvo.e		= &e
	vcvo.ename	= ename
	
	for (i=1; i<=K1; i++) {

// RMSEs for first stage start in SECOND row/column (first has RF)
		rmse = rmseall[i+1,i+1]
// Shea partial R2
		shea = sheaall[1,i]
// first-stage coeffs for ith X1.
// (nb: first column is reduced form eqn for y)
		b=bz[., (i+1)]
// Classical Wald stat (chi2 here); also yields partial R2
// Since r is an L1 x 1 zero vector, can use Rb instead of (Rb-r)
		Rb = b[ (1::L1), . ]
		RVR = V[| 1+i*L,1+i*L \ i*L+L1, i*L+L1 |]
		Wald = Rb' * cholsolve(RVR, Rb)
// Wald stat has dofminus correction in it via sigmasq,
// so remove it to calculate partial R2
		pr2 = (Wald/(N-dofminus)) / (1 + (Wald/(N-dofminus)))

// Robustify F stat if necessary.
		if ((robust~="") | (clustvarname~="") | (kernel~="") | (sw~="")) {
			RVR = Vr[| 1+i*L,1+i*L \ i*L+L1, i*L+L1 |]
			Wald = Rb' * cholsolve(RVR, Rb)
		}
// small dof adjustment is effectively additional L2, e.g., partialled-out regressors
		df = L1
		if (clustvarname=="") {
			df_r = (N-dofminus-L-sdofminus)
			F = Wald / (N-dofminus) * df_r / df
		}
		else {
			df_r = N_clust - 1
			F = Wald / (N-1) * (N-L-sdofminus) * (N_clust - 1) / N_clust / df
		}
		pvalue = Ftail(df, df_r, F)

// If #endog=1, AP=SW=standard F stat
		if (K1==1) {
			Fdf1	= df
			Fdf2	= df_r
			SWF		= F
			SWFp	= pvalue
			SWchi2	= Wald
			SWchi2p	= chi2tail(Fdf1, SWchi2)
			SWr2	= pr2
			APF		= SWF
			APFp	= SWFp
			APchi2	= SWchi2
			APchi2p	= SWchi2p
			APr2	= SWr2
		}
		else {

// Angrist-Pischke and Sanderson-Windmeijer stats etc.
// select matrix needed for both; will select all but the endog regressor of interest
			selmat=J(1,K,1)
			selmat[1,i]=0		// don't select endog regressor of interest
	
// AP
// QXhXh is crossproduct of X1hats (fitted Xs) plus Z2s
// QXhX1 is crossproduct with X1s
// QXhXhi and QXhX1i remove the row/col for the endog regressor of interest
			QXhXhi = select(select(QXhXh,selmat)', selmat)
			QXhX1i = select(QXhX1[.,i], selmat')
// 1st step - in effect, 2nd stage of 2SLS using FITTED X1hats, and then get residuals e1
			b1=cholsolve(QXhXhi, QXhX1i)
			QXhXhinv = invsym(QXhXhi)	//  Need this for V
			b1=QXhXhinv*QXhX1i
			e1 = X1[.,i] - select(Xhat,selmat)*b1
// 2nd step - regress e1 on all Zs and test excluded ones
			QZe1 = quadcross(Z, wf*wvar, e1 ) / N
			b2=cholsolve(QZZ, QZe1)
			APe2 = e1 - Z*b2
			ee = quadcross(APe2, wf*wvar, APe2)
			sigmasq2 = ee / (N-dofminus)
// Classical V
			Vi = sigmasq2 * QZZinv / N
			APRb=b2[ (1::L1), .]
			APRVR = Vi[ (1::L1), (1..L1) ]
			APWald = APRb' * cholsolve(APRVR, APRb)
// Wald stat has dofminus correction in it via sigmasq,
// so remove it to calculate partial R2
			APr2 = (APWald/(N-dofminus)) / (1 + (APWald/(N-dofminus)))
	
// Now SW stat
// Uses same 2SLS coeffs as AP but resids use ACTUAL X1 (not fitted X1)
			e1 = X1[.,i] - select(X,selmat)*b1
// next step - regress e on all Zs and test excluded ones
			QZe1 = quadcross(Z, wf*wvar, e1 ) / N
			b2=cholsolve(QZZ, QZe1)
			SWe2 = e1 - Z*b2
			ee = quadcross(SWe2, wf*wvar, SWe2)
			sigmasq2 = ee / (N-dofminus)
			Vi = sigmasq2 * QZZinv / N
			SWRb=b2[ (1::L1), .]
			SWRVR = Vi[ (1::L1), (1..L1) ]
			SWWald = SWRb' * cholsolve(SWRVR, SWRb)
// Wald stat has dofminus correction in it via sigmasq,
// so remove it to calculate partial R2
			SWr2 = (SWWald/(N-dofminus)) / (1 + (SWWald/(N-dofminus)))
	
// Having calculated AP and SW R-sq based on non-robust Wald, now get robust Wald if needed.
			if ((robust~="") | (clustvarname~="") | (kernel~="") | (sw~="")) {
				e[.,1]=APe2
				omega=m_omega(vcvo)
				Vi = makesymmetric(QZZinv * omega * QZZinv) / N
				APRVR = Vi[ (1::L1), (1..L1) ]
				APWald = APRb' * cholsolve(APRVR, APRb)				//  re-use APRb
				e[.,1]=SWe2
				omega=m_omega(vcvo)
				Vi = makesymmetric(QZZinv * omega * QZZinv) / N
				SWRVR = Vi[ (1::L1), (1..L1) ]
				SWWald = SWRb' * cholsolve(SWRVR, SWRb)				//  re-use SWRb
			}
	
// small dof adjustment is effectively additional L2, e.g., partialled-out regressors
			Fdf1 = (L1-K1+1)
			if (clustvarname=="") {
				Fdf2 = (N-dofminus-L-sdofminus)
				APF = APWald / (N-dofminus) * Fdf2 / Fdf1
				SWF = SWWald / (N-dofminus) * Fdf2 / Fdf1
			}
			else {
				Fdf2 = N_clust - 1
				APF = APWald / (N-1) * (N-L-sdofminus) * (N_clust - 1) / N_clust / Fdf1
				SWF = SWWald / (N-1) * (N-L-sdofminus) * (N_clust - 1) / N_clust / Fdf1
			}
			APFp = Ftail(Fdf1, Fdf2, APF)
			APchi2 = APWald
			APchi2p = chi2tail(Fdf1, APchi2)
			SWFp = Ftail(Fdf1, Fdf2, SWF)
			SWchi2 = SWWald
			SWchi2p = chi2tail(Fdf1, SWchi2)
		}

// Assemble results
		firstmat = firstmat ,													///
					(rmse \ shea \ pr2 \ F \ df \ df_r \ pvalue 				///
					\ SWF \ Fdf1 \ Fdf2 \ SWFp \ SWchi2 \ SWchi2p \ SWr2		///
					\ APF \ Fdf1 \ Fdf2 \ APFp \ APchi2 \ APchi2p \ APr2)
	} // end of loop for an X1 variable

st_numscalar("r(rmse_rf)", rmseall[1,1])
	st_matrix("r(firstmat)", firstmat)
	st_matrix("r(b)", bz)
	st_matrix("r(V)", Vr)
	st_matrix("r(S)", omegar)
	st_numscalar("r(archi2)", ARchi2)
	st_numscalar("r(archi2p)", ARchi2p)
	st_numscalar("r(arf)", ARF)
	st_numscalar("r(arfp)", ARFp)
	st_numscalar("r(ardf)", ARdf)
	st_numscalar("r(ardf_r)", ARdf2)
	st_numscalar("r(sstat)",sstat)
	st_numscalar("r(sstatp)",sstatp)
	st_numscalar("r(sstatdf)",sstatdf)

} // end program s_ffirst

// **********************************************************************

void s_omega(
 					string scalar ZZmatrix,
					string scalar ename,
					string scalar Znames,
					string scalar touse,
					string scalar weight,
					string scalar wvarname,
					scalar wf,
					scalar N,
					string scalar robust,
					string scalar clustvarname,
					string scalar clustvarname2,
					string scalar clustvarname3,
					scalar bw,
					string scalar kernel,
					string scalar sw,
					string scalar psd,
					string scalar ivarname,
					string scalar tvarname,
					string scalar tindexname,
					scalar tdelta,
					scalar center,
					scalar dofminus)
{

struct ms_vcvorthog scalar vcvo

st_view(Z,      ., st_tsrevar(tokens(Znames)),  touse)
	st_view(wvar,   ., st_tsrevar(wvarname), touse)
	st_view(e,      ., vcvo.ename,  touse)
	
	vcvo.e		= &e
	vcvo.Z		= &Z
	vcvo.wvar	= &wvar

ZZ = st_matrix(ZZmatrix)

S=m_omega(vcvo)

st_matrix("r(S)", S)
} // end of s_omega program

// Mata utility for sequential use of solvers
// Default is cholesky;
// if that fails, use QR;
// if overridden, use QR.

function cholqrsolve (	numeric matrix A,
						numeric matrix B,
						| real scalar useqr)
{
	if (args()==2) useqr = 0
	
	real matrix C

if (!useqr) {
		C = cholsolve(A, B)
		if (C[1,1]==.) {
			C = qrsolve(A, B)
		}
	}
	else {
		C = qrsolve(A, B)
	}

return(C)

}

end			//  end Mata section

exit		//  exit before loading comments

********************************** VERSION COMMENTS **********************************
*  Initial version cloned from official ivreg version 5.0.9  19Dec2001
*  1.0.2:  add logic for reg3. Sargan test
*  1.0.3:  add prunelist to ensure that count of excluded exogeneous is correct 
*  1.0.4:  revise option to exog(), allow included exog to be specified as well
*  1.0.5:  switch from reg3 to regress, many options and output changes
*  1.0.6:  fixed treatment of nocons in Sargan and C-stat, and corrected problems
*          relating to use of nocons combined with a constant as an IV
*  1.0.7:  first option reports F-test of excluded exogenous; prunelist bug fix
*  1.0.8:  dropped prunelist and switched to housekeeping of variable lists
*  1.0.9:  added collinearity checks; C-stat calculated with recursive call;
*          added ffirst option to report only F-test of excluded exogenous
*          from 1st stage regressions
*  1.0.10: 1st stage regressions also report partial R2 of excluded exogenous
*  1.0.11: complete rewrite of collinearity approach - no longer uses calls to
*          _rmcoll, does not track specific variables dropped; prunelist removed
*  1.0.12: reorganised display code and saved results to enable -replay()-
*  1.0.13: -robust- and -cluster- now imply -small-
*  1.0.14: fixed hascons bug; removed ivreg predict fn (it didn't work); allowed
*          robust and cluster with z stats and correct dofs
*  1.0.15: implemented robust Sargan stat; changed to only F-stat, removed chi-sq;
*          removed exog option (only orthog works)
*  1.0.16: added clusterised Sargan stat; robust Sargan handles collinearities;
*          predict now works with standard SE options plus resids; fixed orthog()
*          so it accepts time series operators etc.
*  1.0.17: fixed handling of weights.  fw, aw, pw & iw all accepted.
*  1.0.18: fixed bug in robust Sargan code relating to time series variables.
*  1.0.19: fixed bugs in reporting ranks of X'X and Z'Z
*          fixed bug in reporting presence of constant
*  1.0.20: added GMM option and replaced robust Sargan with (equivalent) J;
*          added saved statistics of 1st stage regressions
*  1.0.21: added Cragg HOLS estimator, including allowing empty endog list;
*          -regress- syntax now not allowed; revised code searching for "_cons"
*  1.0.22: modified cluster output message; fixed bug in replay for Sargan/Hansen stat;
*          exactly identified Sargan/Hansen now exactly zero and p-value not saved as e();
*          cluster multiplier changed to 1 (from buggy multiplier), in keeping with
*          eg Wooldridge 2002 p. 193.
*  1.0.23: fixed orthog option to prevent abort when restricted equation is underid.
*  1.0.24: fixed bug if 1st stage regressions yielded missing values for saving in e().
*  1.0.25: Added Shea version of partial R2
*  1.0.26: Replaced Shea algorithm with Godfrey algorithm
*  1.0.27: Main call to regress is OLS form if OLS or HOLS is specified; error variance
*          in Sargan and C statistics use small-sample adjustment if -small- option is
*          specified; dfn of S matrix now correctly divided by sample size
*  1.0.28: HAC covariance estimation implemented
*          Symmetrize all matrices before calling syminv
*          Added hack to catch F stats that ought to be missing but actually have a
*          huge-but-not-missing value
*          Fixed dof of F-stat - was using rank of ZZ, should have used rank of XX (couldn't use df_r
*          because it isn't always saved.  This is because saving df_r triggers small stats
*          (t and F) even when -post- is called without dof() option, hence df_r saved only
*          with -small- option and hence a separate saved macro Fdf2 is needed.
*          Added rankS to saved macros
*          Fixed trap for "no regressors specified"
*          Added trap to catch gmm option with no excluded instruments
*          Allow OLS syntax (no endog or excluded IVs specified)
*          Fixed error messages and traps for rank-deficient robust cov matrix; includes
*          singleton dummy possibility
*          Capture error if posting estimated VCV that isn't pos def and report slightly
*          more informative error message
*          Checks 3 variable lists (endo, inexog, exexog) separately for collinearities
*          Added AC (autocorrelation-consistent but conditionally-homoskedastic) option
*          Sargan no longer has small-sample correction if -small- option
*          robust, cluster, AC, HAC all passed on to first-stage F-stat
*          bw must be < T
*  1.0.29  -orthog- also displays Hansen-Sargan of unrestricted equation
*          Fixed collinearity check to include nocons as well as hascons
*          Fixed small bug in Godfrey-Shea code - macros were global rather than local
*          Fixed larger bug in Godfrey-Shea code - was using mixture of sigma-squares from IV and OLS
*            with and without small-sample corrections
*          Added liml and kclass
*  1.0.30  Changed order of insts macro to match saved matrices S and W
*  2.0.00  Collinearities no longer -qui-
*          List of instruments tested in -orthog- option prettified
*  2.0.01  Fixed handling of nocons with no included exogenous, including LIML code
*  2.0.02  Allow C-test if unrestricted equation is just-identified.  Implemented by
*          saving Hansen-Sargan dof as = 0 in e() if just-identified.
*  2.0.03  Added score() option per latest revision to official ivreg
*  2.0.04  Changed score() option to pscore() per new official ivreg
*  2.0.05  Fixed est hold bug in first-stage regressions
*          Fixed F-stat finite sample adjustment with cluster option to match official Stata
*          Fixed F-stat so that it works with hascons (collinearity with constant is removed)
*          Fixed bug in F-stat code - wasn't handling failed posting of vcv
*          No longer allows/ignores nonsense options
*  2.0.06  Modified lsStop to sync with official ivreg 5.1.3
*  2.0.07a Working version of CUE option
*          Added sortpreserve, ivar and tvar options
*          Fixed smalls bug in calculation of T for AC/HAC - wasn't using the last ob
*          in QS kernel, and didn't take account of possible dropped observations
*  2.0.07b Fixed macro bug that truncated long varlists
*  2.0.07c Added dof option.
*          Changed display of RMSE so that more digits are displayed (was %8.1g)
*          Fixed small bug where cstat was local macro and should have been scalar
*          Fixed bug where C stat failed with cluster.  NB: wmatrix option and cluster are not compatible!
*  2.0.7d  Fixed bug in dof option
*  2.1.0   Added first-stage identification, weak instruments, and redundancy stats
*  2.1.01  Tidying up cue option checks, reporting of cue in output header, etc.
*  2.1.02  Used Poskitt-Skeels (2002) result that C-D eval = cceval / (1-cceval)
*  2.1.03  Added saved lists of separate included and excluded exogenous IVs
*  2.1.04  Added Anderson-Rubin test of signif of endog regressors
*  2.1.05  Fix minor bugs relating to cluster and new first-stage stats
*  2.1.06  Fix bug in cue: capture estimates hold without corresponding capture on estimates unhold
*  2.1.07  Minor fix to ereturn local wexp, promote to version 8.2
*  2.1.08  Added dofminus option, removed dof option.  Added A-R test p-values to e().
*          Minor bug fix to A-R chi2 test - was N chi2, should have been N-L chi2.
*          Changed output to remove potentially misleading refs to N-L etc.
*          Bug fix to rhs count - sometimes regressors could have exact zero coeffs
*          Bug fix related to cluster - if user omitted -robust-, orthog would use Sargan and not J
*          Changed output of Shea R2 to make clearer that F and p-values do not refer to it
*          Improved handling of collinearites to check across inexog, exexog and endo lists
*          Total weight statement moved to follow summ command
*          Added traps to catch errors if no room to save temporary estimations with _est hold
*          Added -savefirst- option. Removed -hascons-, now synonymous with -nocons-.
*  2.1.09  Fixes to dof option with cluster so it no longer mimics incorrect areg behavior
*          Local ivreg2cmd to allow testing under name ivreg2
*          If wmatrix supplied, used (previously not used if non-robust sargan stat generated)
*          Allowed OLS using (=) syntax (empty endo and exexog lists)
*          Clarified error message when S matrix is not of full rank
*          cdchi2p, ardf, ardf_r added to saved macros
*          first and ffirst replay() options; DispFirst and DispFFirst separately codes 1st stage output
*          Added savefprefix, macro with saved first-stage equation names.
*          Added version option.
*          Added check for duplicate variables to collinearity checks
*          Rewrote/simplified Godfrey-Shea partial r2 code
* 2.1.10   Added NOOUTput option
*          Fixed rf bug so that first does not trigger unnecessary saved rf
*          Fixed cue bug - was not starting with robust 2-step gmm if robust/cluster
* 2.1.11   Dropped incorrect/misleading dofminus adjustments in first-stage output summary
* 2.1.12   Collinearity check now checks across inexog/exexog/endog simultaneously
* 2.1.13   Added check to catch failed first-stage regressions
*          Fixed misleading failed C-stat message
* 2.1.14   Fixed mishandling of missing values in AC (non-robust) block
* 2.1.15   Fixed bug in RF - was ignoring weights
*          Added -endog- option
*          Save W matrix for all cases; ensured copy is posted with wmatrix option so original isn't zapped
*          Fixed cue bug - with robust, was entering IV block and overwriting correct VCV
* 2.1.16   Added -fwl- option
*          Saved S is now robust cov matrix of orthog conditions if robust, whereas W is possibly non-robust
*          weighting matrix used by estmator.  inv(S)=W if estimator is efficient GMM.
*          Removed pscore option (dropped by official ivreg).
*          Fixed bug where -post- would fail because of missing values in vcv
*          Remove hascons as synonym for nocons
*          OLS now outputs 2nd footer with variable lists
* 2.1.17   Reorganization of code
*          Added ll() macro
*          Fixed N bug where weights meant a non-integer ob count that was rounded down
*          Fixed -fwl- option so it correctly handles weights (must include when partialling-out)
*          smatrix option takes over from wmatrix option.  Consistent treatment of both.
*          Saved smatrix and wmatrix now differ in case of inefficient GMM.
*          Added title() and subtitle() options.
*          b0 option returns a value for the Sargan/J stat even if exactly id'd.
*          (Useful for S-stat = value of GMM objective function.)
*          HAC and AC now allowed with LIML and k-class.
*          Collinearity improvements: bug fixed because collinearity was mistakenly checked across
*          inexog/exexog/endog simultaneously; endog predicted exactly by IVs => reclassified as inexog;
*          _rmcollright enforces inexog>endo>exexog priority for collinearities, if Stata 9.2 or later.
*          K-class, LIML now report Sargan and J.  C-stat based on Sargan/J.  LIML reports AR if homosked.
*          nb: can always easily get a C-stat for LIML based on diff of two AR stats.
*          Always save Sargan-Hansen as e(j); also save as e(sargan) if homoskedastic.
*          Added Stock-Watson robust SEs options sw()
* 2.1.18   Added Cragg-Donald-Stock-Yogo weak ID statistic critical values to main output
*          Save exexog_ct, inexog_ct and endog_ct as macros
*          Stock-Watson robust SEs now assume ivar is group variable
*          Option -sw- is standard SW.  Option -swpsd- is PSD version a la page 6 point 10.
*          Added -noid- option.  Suppresses all first-stage and identification statistics.
*          Internal calls to ivreg2 use noid option.
*          Added hyperlinks to ivreg2.hlp and helpfile argument to display routines to enable this.
* 2.1.19   Added matrix rearrangement and checks for smatrix and wmatrix options
*          Recursive calls to cstat simplified - no matrix rearrangement or separate robust/nonrobust needed
*          Reintroduced weak ID stats to ffirst output
*          Added robust ID stats to ffirst output for case of single endogenous regressor
*          Fixed obscure bug in reporting 1st stage partial r2 - would report zero if no included exogenous vars
*          Removed "HOLS" in main output (misleading if, e.g., estimation is AC but not HAC)
*          Removed "ML" in main output if no endogenous regressors - now all ML is labelled LIML
*          model=gmm is now model=gmm2s; wmatrix estimation is model=gmm
*          wmatrix relates to gmm estimator; smatrix relates to gmm var-cov matrix; b0 behavior equiv to wmatrix
*          b0 option implies nooutput and noid options
*          Added nocollin option to skip collinearity checks
*          Fixed minor display bug in ffirst output for endog vars with varnames > 12 characters
*          Fixed bug in saved rf and first-stage results for vars with long varnames; uses permname
*          Fixed bug in model df - had counted RHS, now calculates rank(V) since latter may be rank-deficient
*          Rank of V now saved as macro rankV
*          fwl() now allows partialling-out of just constant with _cons
*          Added Stock-Wright S statistic (but adds overhead - calls preserve)
*          Properties now include svyj.
*          Noted only: fwl bug doesn't allow time-series operators.
* 2.1.20   Fixed Stock-Wright S stat bug - didn't allow time-series operators
* 2.1.21   Fixed Stock-Wright S stat to allow for no exog regressors cases
* 2.2.00   CUE partials out exog regressors, estimates endog coeffs, then exog regressors separately - faster
*          gmm2s becomes standard option, gmm supported as legacy option
* 2.2.01   Added explanatory messages if gmm2s used.
*          States if estimates efficient for/stats consistent for het, AC, etc.
*          Fixed small bug that prevented "{help `helpfile'##fwl:fwl}" from displaying when -capture-d.
*          Error message in footer about insuff rank of S changed to warning message with more informative message.
*          Fixed bug in CUE with weights.
* 2.2.02   Removed CUE partialling-out; still available with fwl
*          smatrix and wmatrix become documented options. e(model)="gmmw" means GMM with arbitrary W
* 2.2.03   Fixed bug in AC with aweights; was weighting zi'zi but not ei'ei.
* 2.2.04   Added abw code for bw(), removed properties(svyj)
* 2.2.05   Fixed bug in AC; need to clear variable vt1 at start of loop
*          If iweights, N (#obs with precision) rounded to nearest integer to mimic official Stata treatment
*          and therefore don't need N scalar at all - will be same as N
*          Saves fwl_ct as macro.
*          -ffirst- output, weak id stat, etc. now adjust for number of partialled-out variables.
*          Related changes: df_m, df_r include adjustments for partialled-out variables.
*          Option nofwlsmall introduced - suppresses above adjustments.  Undocumented in ivreg2.hlp.
*          Replaced ID tests based on canon corr with Kleibergen-Paap rk-based stats if not homoskedastic
*          Replaced LR ID test stats with LM test stats.
*          Checks that -ranktest- is installed.
* 2.2.06   Fixed bug with missing F df when cue called; updated required version of ranktest 
* 2.2.07   Modified redundancy test statistic to match standard regression-based LM tests
*          Change name of -fwl- option to -partial-.
*          Use of b0 means e(model)=CUE.  Added informative b0 option titles. b0 generates output but noid.
*          Removed check for integer bandwidth if auto option used.
* 2.2.08   Add -nocollin- to internal calls and to -ivreg2_cue- to speed performance.
* 2.2.09   Per msg from Brian Poi, Alastair Hall verifies that Newey-West cited constant of 1.1447
*          is correct. Corrected mata abw() function. Require -ranktest- 1.1.03.
* 2.2.10   Added Angrist-Pischke multivariate f stats.  Rewrite of first and ffirst output.
*          Added Cragg-Donald to weak ID output even when non-iid.
*          Fixed small bug in non-robust HAC code whereby extra obs could be used even if dep var missing.
*             (required addition of  L`tau'.(`s1resid') in creation of second touse variable)
*          Fixed bugs that zapped varnames with "_cons" in them
*          Changed tvar and ivar setup so that data must be tsset or xtset.
*          Fixed bug in redundancy test stat when called by xtivreg2+cluster - no dofminus adj needed in this case
*          Changed reporting so that gaps between panels are not reported as such.
*          Added check that weight variable is not transformed by partialling out.
*          Changed Stock-Wright S statistic so that it uses straight partialling-out of exog regressors
*            (had been, in effect, doing 2SGMM partialling-out)
*          Fixed bug where dropped collinear endogenous didn't get a warning or listing
*          Removed N*CDEV Wald chi-sq statistic from ffirst output (LM stat enough)
* 3.0.00   Fully rewritten and Mata-ized code.  Require min Stata 10.1 and ranktest 1.2.00.
*          Mata support for Stock-Watson SEs for fixed effects estimator; doesn't support fweights.
*          Changed handling of iweights yielding non-integer N so that (unlike official -regress-) all calcs
*          for RMSE etc. use non-integer N and N is rounded down only at the end.
*          Added support for Thompson/Cameron-Gelbach-Miller 2-level cluster-robust vcvs.
* 3.0.01   Now exits more gracefully if no regressors survive after collinearity checks
* 3.0.02   -capture- instead of -qui- before reduced form to suppress not-full-rank error warning
*          Modified Stock-Wright code to partial out all incl Xs first, to reduce possibility of not-full-rank
*          omega and missing sstat.  Added check within Stock-Wright code to catch not-full-rank omega.
*          Fixed bug where detailed first-stage stats with cluster were disrupted if data had been tsset
*          using a different variables.
*          Fixed bug that didn't allow regression on just a constant.
*          Added trap for no observations.
*          Added trap for auto bw with panel data - not allowed.
* 3.0.03   Fixed bug in m_omega that always used Stock-Watson spectral decomp to create invertible shat
*          instead of only when (undocumented) spsd option is called.
*          Fixed bug where, if matsize too small, exited with wrong error (mistakenly detected as collinearities)
*          Removed inefficient call to -ranktest- that unnecessarily requested stats for all ranks, not just full.
* 3.0.04   Fixed coding error in m_omega for cluster+kernel.  Was *vcvo.e[tmatrix[.,1]], should have been (*vcvo.e)[tmatrix[.,1]].
*          Fixed bug whereby clusters defined by strings were not handled correctly.
*          Updated ranktest version check
* 3.0.05   Added check to catch unwanted transformations of time or panel variables by partial option.
* 3.0.06   Fixed partial bug - partialcons macro saved =0 unless _cons explicitly in partial() varlist
* 3.0.07   kclass was defaulting to LIML - fixed.
*          Renamed spsd option to psda (a=abs) following Stock-Watson 2008. Added psd0 option following Politis 2007.
*          Fixed bug that would prevent RF and first-stage with cluster and TS operators if cluster code changed sort order.
*          Modified action if S matrix is not full rank and 2-step GMM chosen.  Now continue but report problem in footer
*          and do not report J stat etc.
* 3.0.08   Fixed cluster+bw; was not using all observations of all panel units if panel was unbalanced.
*          Fixed inconsequential bug in m_omega that caused kernel loop to be entered (with no impact) even if kernel==""
*          Fixed small bug that compared bw to T instead of (correctly) to T/delta when checking that bw can't be too long.
*          Added dkraay option = cluster on t var + kernel-robust
*          Added kiefer option = truncated kernel, bw=T (max), and no robust
*          Fixed minor reporting bug that reported time-series gaps in entire panel dataset rather than just portion touse-d.
*          Recoded bw and kernel checks into subroutine vkernel.  Allow non-integer bandwidth within check as in ranktest.
* 3.1.01   First ivreg2 version with accompanying Mata library (shared with -ranktest-).  Mata library includes
*          struct ms_vcvorthog, m_omega, m_calckw, s_vkernel.
*          Fixed bug in 2-way cluster code (now in m_omega in Mata library) - would crash if K>1 (relevant for -ranktest- only).
* 3.1.02   Converted cdsy to Mata code and moved to Mata library. Standardized spelling/caps/etc. of QS as "Quadratic Spectral".
* 3.1.03   Improved partialling out in s_sstat and s_ffirst: replaced qrsolve with invsym.
* 3.1.04   Fixed minor bug in s_crossprod - would crash with L1=0 K1>0, and also with K=0
* 3.1.05   Fixed minor bug in orthog - wasn't saving est results if eqn w/o suspect instruments did not execute properly
*          Fixed minor bug in s_cccollin() - didn't catch perverse case of K1>0 (endog regressors) and L1=0 (no excl IVs)
* 3.1.06   Spelling fix for Danielle kernel, correct error check for bw vs T-1
* 3.1.07   Fixed bug that would prevent save of e(sample) when partialling out just a constant
* 3.1.08   01Jan14. Fixed reporting bug with 2-way clustering and kernel-robust that would give wrong count for 2nd cluster variable.
* 3.1.09   13July14. _rmcollright under version control has serious bug for v10 and earlier. Replaced with canon corr approach.
*          Fixed obscure bug in estimation sample - was not using obs when tsset tvar is missing, even if TS operators not used.
*          Fixed bug in auto bw code so now ivreg2 and ivregress agree. Also, ivreg2 auto bw code handles gaps in TS correctly.
* 4.0.00   25Jan15. Promote to require Stata version 11.2
*          Rewrite of s_gmm1s, s_iegmm, s_egmm etc. to use matrix solvers rather than inversion.
*          rankS and rankV now calculated along with estimators; rankS now always saved.
*          Returned to use of _rmcollright to detect collinearities since bug was in Stata 10's _rmcollright and now not relevant.
*          Added reporting of collinearities and duplicates in replay mode.
*          Rewrite of legacy support for previous ivreg2x version.  Main program calls ivreg2x depending on _caller().
*          Estimation and replay moved to ivreg211 subroutine above.
* 4.0.01   8Feb15. Fixed bug in default name and command used used for saved first and RF equations
*          Fixed bug in saved command line (was ivreg211, should be ivreg2).
* 4.0.02   9Feb15. Changed forced exit at Stata <11 before continuing loading to forced exit pre-Mata code at Stata <9.
* 4.1.00   Substantial rewrite to allow factor variables. Now also accepts TS ops as well as FV ops in partial varlist.
*          Rewrite included code for dropped/collinear/reclassified.
*          Saved RF and 1st-stage estimations have "if e(sample)" instead of "if `touse'" in e(cmdline).
*          Rewrite of s_gmm1s etc. to use qrsolve if weighting matrix not full rank or cholsolve fails
*          Fixed bug in display subroutines that would display hyperlink to wrong (nonexistent) help file.
* 4.1.01   15Jun15. Fixed bug that did not allow dropped variables to be in partial(.) varlist.
*          Major rewrite of parsing code and collinearity/dropped/reclassified code.
*          Added support for display options noomitted, vsquish, noemptycells, baselevels, allbaselevels.
*          Changed from _rmcoll/_rmcollright/_rmcoll2list to internal ivreg2_rmcollright2
*          Changed failure of ranktest to obtain id stats to non-fatal so that estimation proceeds.
*          Removed recount via _rmcoll if noid option specified
*          Added partial(_all) option.
*          Improved checks of smatrix, wmatrix, b0 options
*          Rewrite of first-stage and reduced form code; rewrite of replay(.) functionality
*          Added option for displaying system of first-stage/reduced form eqns.
*          Replaced AP first-stage test stats with SW (Sanderson-Windmeijer) first-stage stats
*          Corrected S LM stat option; now calcuated in effect as J stat for case of no endog (i.e. b=0)
*          with inexog partialled out i.e. LM version of AR stat; now matches weakiv
*          Undocumented FV-related options: fvsep (expand endo, inexog, exexog separately) fvall (expand together)
* 4.1.02   17Jun15.  Fixed bug in collinearity check - was ignoring weights.
*          More informative error message if invalid matrix provided to smatrix(.) or wmatrix(.) options.
*          Caught error if depvar was FV or TS var that expanded to >1 variable.
* 4.1.03   18Jun15.  Fixed bug with robust + rf option.
* 4.1.04   18Jun15.  Fixed bug in AR stat with dofminus option + cluster (was subtracting dof, shouldn't).
* 4.1.05   18Jun15.  Added rmse, df_m, df_r to saved RF and first-stage equation results.
* 4.1.06   4July15.  Replaced mvreg with Mata code for partialling out (big speed gains with many vars).
*          Rewrote AddOmitted to avoid inefficient loop; replaced with Mata subscripting.
*          Failure of id stats because of collinearities triggers error message only; estimation continues.
*          Calculation of dofs etc. uses rankS and rankV instead of iv1_ct and rhs1_ct;
*          counts are therefore correct even in presence of collinearities and use of nocollin option.
*          nocollin options triggers use of QR instead of default Cholesky.
*          rankxx and rankzz now based on diag0cnt of (XX)^-1 and (ZZ)^-1.
*          CUE fails if either S or V not full rank; can happen if nocollin option used.
*          Added undocumented useqr option to force use of QR instead of Cholesky.
*          Misc other code tweaks to make results more robust to nocollin option.
* 4.1.07   12July15.  Fixed bugs in calculation of rank(V) (had miscounted in some cases if omega not full rank)
*          Changed calc of dofs etc. from rankS and rankV to rankzz and rankxx (had miscounted in some cases etc.).
*          Restored warning message for all exog regressors case if S not full rank.
* 4.1.08   27July15. Replaced wordcount(.) function with word count macro in AddOmitted;
*          AddOmitted called only if any omitted regressors to add.
*          Added center option for centering moments.
* 4.1.09   20Aug15. Expanded error message for failure to save first-stage estimations (var name too long).
*          Fixed bug when weighting used with new partial-out code (see 4.1.06 4July15).
*          Tweaked code so that if called under Stata version < 11, main ivreg2.ado is exited immediately after
*          loading parent ivreg2 program.  Removed automatic use of QR solver when nocollin option used.
*          Added saved condition numbers for XX and ZZ.
*          e(cmdline) now saves original string including any "s (i.e., saves `0' instead of `*').
* 4.1.10   Fixed bug with posting first-stage results if sort had been disrupted by Mata code.
*          Fixed bug which mean endog(.) and orthog(.) varlists weren't saved or displayed.

Texto alterado

Open file

*! ivreghdfe 1.0.0  07jul2018
*! this just adds absorb() to this code:
*! ivreg2 4.1.10  9Feb2016
*! authors cfb & mes
*! see end of file for version comments

* Parent program, forks to versions as appropriate after version call
* Requires byable(onecall)
program define ivreghdfe, eclass byable(onecall) /* properties(svyj) */ sortpreserve
	local lversion 04.1.10

* Minimum of version 8 required for parent program (earliest ivreg2 is ivreg28)
	version 8

* replay can't be combined with by
	if replay() & _by() {
di as err "invalid syntax - cannot use by with replay"
		exit 601
	}

* Handling of by. ivreg2x programs are byable(recall), so must set prefix for them.
	if _by() {
		local BY `"by `_byvars'`_byrc0':"'
	}

}
		else {
* If not replay, call legacy ivreg2 and then add macros
			`BY' `ivreg2cmd' `0'
			ereturn local cmd "ivreghdfe"
			ereturn local ivreg2cmd "`ivreg2cmd'"
			ereturn local version `lversion'
			ereturn local predict ivreg2_p
		}
		exit
	}

// Version is 11 or above.
// Pass control to current estimation program ivreg211.
	if replay() {
		ivreg211 `0'
	}
// If not replay, call ivreg211 and then add macros
	else {
		`BY' ivreg211 `0'
		cap mata: mata drop HDFE // prefix ivreg211 call with capture?
		cap mata: mata drop hdfe_residuals
		ereturn local cmd "ivreghdfe"
		ereturn local ivreg2cmd "ivreghdfe"
		ereturn local version `lversion'
		ereturn local predict ivreg2_p
		if (e(N_hdfe)!=  .) ereturn local predict reghdfe_p
	}

********************* EXIT IF STATA VERSION < 11 ********************************

* When do file is loaded, exit here if Stata version calling program is < 11.
* Prevents loading of rest of program file (could cause earlier Statas to crash).

if c(stata_version) < 11 {
	exit
}

******************** END EXIT IF STATA VERSION < 11 *****************************

* Main estimation program
program define ivreg211, eclass byable(recall) sortpreserve
	version 11.2
	local ranktestversion 01.3.02

local ivreg2cmd "ivreg211"			//  actual command name
	local ivreg2name "ivreg2"			//  name used in command line and for default naming of equations etc.
	local ranktestcmd "ranktest"

if replay() {
		syntax [, 												///
				FIRST FFIRST RF SFIRST							///
				dropfirst droprf dropsfirst						///
				Level(integer $S_level)							///
				NOHEader NOFOoter								///
				EForm(string) PLUS								///
				NOOMITTED vsquish noemptycells					///
				baselevels allbaselevels 						///
				VERsion											///
				]
		if "`version'" != "" & "`first'`ffirst'`rf'`noheader'`nofooter'`dropfirst'`droprf'`eform'`plus'" != "" {
			di as err "option version not allowed"
			error 198
		}
		if "`version'" != "" {
			di in gr "`lversion'"
			ereturn clear
			ereturn local version `lversion'
			exit
		}
		if `"`e(cmd)'"' != "ivreghdfe"  {
			error 301
		}
// Set display options
		local dispopt	eform(`eform') `noomitted' `vsquish' `noemptycells' `baselevels' `allbaselevels'

//  Confirm ranktest is installed (necessary component).
		checkversion_ranktest `ranktestcmd' `ranktestversion'

// Start parsing
		syntax [anything(name=0)] [if] [in] [aw fw pw iw/] [,		///
				NOID NOCOLLIN										///
				FIRST FFIRST SAVEFIRST SAVEFPrefix(name)			///
				RF SAVERF SAVERFPrefix(name)						///
				SFIRST SAVESFIRST SAVESFPrefix(name)				///
				SMall NOConstant 									///
				Robust CLuster(varlist) kiefer dkraay(integer 0)	///
				BW(string) kernel(string) center					///
				GMM GMM2s CUE										///
				LIML COVIV FULLER(real 0) Kclass(real 0)			///
				ORTHOG(string) ENDOGtest(string) REDundant(string)	///
				PARTIAL(string) FWL(string)							///
				Absorb(string)           							///
				Level(integer $S_level)								///
				NOHEader NOFOoter NOOUTput							///
				bvclean NOOMITTED omitted vsquish noemptycells		///
				baselevels allbaselevels 							///
				title(string) subtitle(string)						///
				DEPname(string) EForm(string) PLUS					///
				Tvar(varname) Ivar(varname)							///
				B0(string) SMATRIX(string) WMATRIX(string)			///
				sw psd0 psda useqr									///
				dofminus(integer 0) sdofminus(integer 0)			///
				NOPARTIALSMALL										///
				fvall fvsep											///
				]

if (`"`absorb'"' != "") {
		// ensure reghdfe is installed
		capture which reghdfe
		if c(rc) {
			di as err "Error: must have reghdfe installed"
			di as err "To install, from within Stata type " _c
			//di in smcl "{stata ssc install reghdfe :ssc install reghdfe}"
			loc url "https://github.com/sergiocorreia/reghdfe/raw/version-4/src/"
			di in smcl `"{stata `"net install reghdfe, from(`url')"':net install reghdfe, from(`url')}"'
			exit 601
		}
		cap reghdfe, requirements // ensure ftools, moresyntax, etc. are installed
		if (c(rc)) reghdfe, requirements
		
		reghdfe, check // ensure the .mlib exists

// absorb implies...
		local small small
		local noconstant noconstant
		local nopartialsmall
	}

***********************************************************

* Create HDFE object and update touse
		if (`"`absorb'"' != "") {
			* fixed_effects(absvars | , touse, wtype, wtvar, dropsing, verbose)
			mata: HDFE = fixed_effects("`absorb'", "`touse'", "`weight'", "`wvar'")
			_assert ("`s(options)'"==""), msg("unsupported options in {bf:absorb()}: `s(options)'")
			mata: HDFE.clustervars = HDFE.base_clustervars = tokens("`cluster'")
			mata: HDFE.num_clusters = length(HDFE.clustervars)
			mata: HDFE.estimate_dof() // compute degrees-of-freedom
			mata: HDFE.save_touse("`touse'", 1) // 1 = overwrite
		}

***************************************************************
// Time-series data
// tindex used by Mata code so that ts operators work correctly

tempvar tindex
		qui gen `tindex'=1 if `touse'
		qui replace `tindex'=sum(`tindex') if `touse'

// kiefer VCV = kernel(tru) bw(T) and no robust with tsset data
		if "`kiefer'" ~= "" {
			local bw	=`T'
		}

*********** Column of ones for constant set up here **************

if "`noconstant'"=="" {
// If macro not created, automatically omitted.
			tempvar ones
			qui gen byte `ones' = 1 if `touse'
		}

// Check that LHS expanded to a single variable
		local wrongvars_ct	: word count `lhs'
		if `wrongvars_ct' > 1 {
di as err "multiple dependent variables specified: `lhs'"
			error 198
		}

// *************** Partial-out block ************** //

// `partial' has all to be partialled out except for constant
		if "`partial1'" != "" | `partialcons'==1 | "`absorb'" != "" {
			preserve

if ("`absorb'" != "") {
				local hdfe_varlist `fv_lhs1' `fv_endo1' `fv_inexog1' `fv_exexog1' `fv_partial1'
				mata: st_store(HDFE.sample, tokens("`hdfe_varlist'"), HDFE.partial_out(tokens("`hdfe_varlist'")))
				mata: st_local("absorb_ct", strofreal(HDFE.df_a))
				assert `absorb_ct'`' != .
				if (`absorb_ct'==0) local absorb_ct 1 // adjustment to match ivreg2 and old reghdfe (happens if absvar is nested in cluster)
				local partial_ct 0
				local partialcons `absorb_ct'
			}

if ("`partial1'" != "" | `partialcons'==1) {
				mata: s_partial	("`fv_lhs1'",			///
								"`fv_endo1'",			///
								"`fv_inexog1'",			///
								"`fv_exexog1'",			///
								"`fv_partial1'",		///
								"`touse'",				///
								"`weight'",				///
								"`wvar'",				///
								`wf',					///
								`N',					///
								`cons')
				local partial_ct : word count `partial1'
			}

// Constant is partialled out, unless nocons already specified in the first place
			capture drop `ones'
			local ones ""
			if "`noconstant'" == "" {
// partial_ct used for small-sample adjustment to regression F-stat
				local partial_ct = `partial_ct' + 1
				local noconstant "noconstant"
				local cons 0
			}

if ("`absorb'" != "") {
				local partial_ct = `partial_ct' + `absorb_ct'
			}
		}
		else {
// Set count of partial vars to zero if option not used
			local partial_ct 0
			local partialcons 0
		}
// Add partial_ct to small dof adjustment sdofminus
		if "`nopartialsmall'"=="" {
			local sdofminus = `sdofminus'+`partial_ct'
		}

*********************************************

scalar `rankzz'	= rowsof(`ZZinv') - diag0cnt(`ZZinv')
		scalar `rankxx'	= rowsof(`XXinv') - diag0cnt(`XXinv')
		local overid	= `rankzz' - `rankxx'

********** CLUSTER SETUP **********************************************

local N_clust=min(`N_clust1',`N_clust2')

}		// end 2-way cluster block
		}		// end cluster block

************************************************************************************************

tempname b W S V beta lambda j jp rss mss rmse sigmasq rankV rankS
		tempname arubin arubinp arubin_lin arubin_linp
		tempname r2 r2_a r2u r2c F Fp Fdf2 ivest

if "`liml'`kclassopt'"~="" {

if "`liml'`kclassopt'`b0'"=="" {

local effic1s =		(													///
										"`gmm2s'`robust'`cluster'`kernel'"==""		///
									|	("`smatrix'"~="" & "`wmatrix'"=="")			///
								)

* Pick up matrix left by s_gmm1s(.)
				tempname QXZ_W_QZX
				mat `QXZ_W_QZX'=r(QXZ_W_QZX)

mat `b'=r(beta)
			mat `S'=r(S)
			mat `W'=r(W)
			mat `V'=r(V)
			scalar `j'=r(j)
			scalar `rss'=r(rss)
			scalar `sigmasq'=r(sigmasq)
			scalar `rankV'=r(rankV)
			scalar `rankS'=r(rankS)

}

****************************************************************
* Done with estimation blocks
****************************************************************

* Store residuals if requested
if (`"`absorb'"' != "") {
	mata: HDFE.residuals = (HDFE.residuals == "" & HDFE.save_any_fe) ? "__temp_reghdfe_resid__" : HDFE.residuals
	mata: st_local("residuals_name", HDFE.residuals)
	if ("`residuals_name'" != "") mata: hdfe_residuals = st_data(., "`resid'", "`touse'")
}

local Fdf1 = `rankxx' - `cons'
		local df_m = `rankxx' - `cons' + (`sdofminus'-`partialcons')

* If j==. or vcv wasn't full rank, then vcv problems and F is meaningless
		if `j' == . | `rc'==506 {
			scalar `F' = .
			scalar `Fp' = .
		}

* End of counts, dofs, F-stat, small sample corrections

if "`sfirst'`savesfirst'`rf'`saverf'`first'`ffirst'`savefirst'" != "" & (`endo1_ct' > 0) {

* Restore original order if changed for mata code above
			capture tsset

capture tsset

// id=underidentification statistic, wid=weak identification statistic
			tempname idrkstat widrkstat iddf idp
			tempname ccf cdf rkf cceval cdeval cd
			tempname idstat widstat

* End of identification stats block

* End of error-checking block

// rankV = rhs1_ct except if nocollin
// rankS = iv1_ct except if nocollin
// nocollin means count may exceed rank (because of dropped vars), so rank #s foolproof

*******************************************************************************************

// restore data if preserved for partial option
		if (`partial_ct' | "`absorb'"!="") {
			restore
		}

ereturn scalar inexog_ct	=`inexog1_ct'
		ereturn scalar exexog_ct	=`exex1_ct'
		ereturn scalar endog_ct		=`endo1_ct'
		ereturn scalar partial_ct	=`partial_ct'

if "`smatrix'" == "" {
			ereturn matrix S `S'
		}
		else {
			ereturn matrix S `S0'					//  it's a copy so original won't be zapped
		}

if "`small'"!="" {
			ereturn scalar df_r=`df_r'
			ereturn local small "small"
		}
		if "`nopartialsmall'"=="" {
			ereturn local partialsmall "small"
		}

if "`robust'`cluster'" != "" {
			ereturn local vcetype "Robust"
		}

if "`redundant'"!="" & "`noid'"=="" {
			ereturn scalar redstat=`redstat'
			ereturn scalar redp=`redp'
			ereturn scalar reddf=`reddf'
			ereturn local  redlist `redundant1'
		}

if ("`absorb'" != "") {
			mata: HDFE.post_footnote()
			assert e(N_hdfe) != .
			mata: st_local("residuals_name", HDFE.residuals)
			if ("`residuals_name'" != "") {
				mata: HDFE.save_variable(HDFE.residuals, hdfe_residuals, "Residuals")
				mata: st_global("e(resid)", HDFE.residuals)
				reghdfe_store_alphas
			}
		}
	}

*******************************************************************************************
* Display results unless ivreg2 called just to generate stats or nooutput option

end

*******************************************************************************************
* SUBROUTINES
*******************************************************************************************

// ************* Display system of or single first-stage and/or RF estimations ************ //

program define DispSFirst, eclass
	args model plus level nofooter helpfile dispopt
	version 11.2

// Display output
	ereturn display, `plus' level(`level') `dispopt'
	ereturn local clustvar `cluster'

end

// ************* Display main estimation outpout ************** //

* Display 2nd footer with variable lists
	if "`nofooter'"=="" {

if (e(N_hdfe)!=  .) reghdfe_footnote
	}
end

**************************************************************************************

// ************ Display collinearity and duplicates warning messages ************ //

// ************* Display all first-stage estimations ************ //

program define DispFirst
	version 11.2
	args helpfile
	tempname firstmat ivest sheapr2 pr2 F df df_r pvalue
	tempname SWF SWFdf1 SWFdf2 SWFp SWr2

di in gr "F test of excluded instruments:"
di in gr "  F(" %3.0f `df'[1,1] "," %6.0f `df_r'[1,1] ") = " in ye %8.2f `F'[1,1]
di in gr "  Prob > F      = " in ye %8.4f `pvalue'[1,1]

di
		}
		_estimates unhold `ivest'
	}
end

// ************* Display list of stored first-stage and RF estimations ************ //

program define DispStored
	args savesfirst saverf savefirst
	version 11.2

// ************* Display summary first-stage and ID test results ************ //

mat `df'      =`firstmat'["df",1]
	mat `df_r'    =`firstmat'["df_r",1]
	mat `SWFdf1'  =`firstmat'["SWFdf1",1]
	mat `SWFdf2'  =`firstmat'["SWFdf2",1]

di
di in gr _newline "Summary results for first-stage regressions"
di in smcl in gr "{hline 43}"
di

di _c in smcl _col(44) "{help `helpfile'##swstats:(Underid)}"
di    in smcl _col(65) "{help `helpfile'##swstats:(Weak id)}"

end

// ************* Post first-stage and/or RF estimations ************ //

// renaming/copying
	local N			= `obs'
	local N_clust	= `nc'
	local N_clust1	= `nc1'
	local N_clust2	= `nc2'
	tempname b V S
	mat `b' = `bmat'
	mat `V' = `vmat'
	mat `S' = `smat'

marksample touse

local dof	= `N' - `iv1_ct' - `dofminus' - `sdofminus'
	ereturn post `b' `V', obs(`obs') esample(`touse') dof(`dof') depname(`depvar') `bfv'

if "`kernel'"!="" {
		ereturn local kernel "`kernel'"
		ereturn scalar bw=`bw'
		ereturn local tvar "`tvar'"
		if "`ivar'" ~= "" {
			ereturn local ivar "`ivar'"
		}
	}

local vce : list clean vce
	local vce = lower("`vce'")
	ereturn local vce `vce'

if "`cluster'"!="" {
		ereturn scalar N_clust=`N_clust'
		ereturn local clustvar `cluster'
	}
	if "`cluster2'"!="" {
		ereturn scalar N_clust1=`N_clust1'
		ereturn scalar N_clust2=`N_clust2'
		ereturn local clustvar1 `cluster1'
		ereturn local clustvar2 `cluster2'
	}

if "`robust'`cluster'" != "" {
		ereturn local vcetype "Robust"
	}

if ("`absorb'" != "") {
		mata: HDFE.post_footnote()
		assert e(N_hdfe) != .
	}
end