ivreg2.ado

*! ivreg2 4.1.11  22Nov2019
*! 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.11

* local to store Stata version of calling program
	local caller = _caller()

* 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
	}

// 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 {
		// use to separate main args from options
		syntax [anything] [if] [in] [aw fw pw iw] [, * ]
		// append caller(.) to options
		`BY' ivreg211 `anything' `if' `in' [`weight' `exp'], `options' caller(`caller')
//		`BY' ivreg211 `0'
		ereturn local cmd "ivreg2"
		ereturn local ivreg2cmd "ivreg2"
		ereturn local version `lversion'
		ereturn local predict ivreg2_p
		ereturn local cmdline ivreg2 `0'		//  `0' rather than `*' in case of any "s in string
	}

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 ivreg2cmd "ivreg211"			//  actual command name
	local ivreg2name "ivreg2"			//  name used in command line and for default naming of equations etc.

	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											///
				caller(real 0)									///
				]
		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

// 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											///
				caller(real 0)										///
				]

//  Confirm ranktest is installed (necessary component).
		checkversion_ranktest `caller'
		local ranktestcmd `r(ranktestcmd)'

// 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'",								///