program define tsls, eclass

	syntax [anything(name=0)] [, fe(varlist) areg robust small cluster(varlist) demean replace first]

	**************************************
	** tsls.ado --                      **
	**************************************
		* By Jake Robbins
		* Data should be demeaned -- otherwise a demean option is available

	********************************
	** Parsing the Syntax         **
	********************************

		local n 0

		gettoken lhs 0 : 0, parse(" ,[") match(paren)
		IsStop `lhs'
		if `s(stop)' {
			error 198
		}
		while `s(stop)'==0 {
			if "`paren'"=="(" {
				local n = `n' + 1
				if `n'>1 {
capture noi error 198
di in red `"syntax is "(all instrumented variables = instrument variables)""'
exit 198
				}
				gettoken p lhs : lhs, parse(" =")
				while "`p'"!="=" {
					if "`p'"=="" {
capture noi error 198
di in red `"syntax is "(all instrumented variables = instrument variables)""'
di in red `"the equal sign "=" is required"'
exit 198
					}
					local endo `endo' `p'
					gettoken p lhs : lhs, parse(" =")
				}
* To enable Cragg HOLS estimator, allow for empty endo list
				local temp_ct  : word count `endo'
				if `temp_ct' > 0 {
					tsunab endo : `endo'
				}
* To enable OLS estimator with (=) syntax, allow for empty exexog list
				local temp_ct  : word count `lhs'
				if `temp_ct' > 0 {
					tsunab exexog : `lhs'
				}
			}
			else {
				local inexog `inexog' `lhs'
			}
			gettoken lhs 0 : 0, parse(" ,[") match(paren)
			IsStop `lhs'
		}
		local 0 `"`lhs' `0'"'

		tsunab inexog : `inexog'
		tokenize `inexog'
		local lhs "`1'"
		local 1 " "
		local inexog `*'

	*************************************
	** Demean!                         **
	*************************************

		if "`demean'"!="" {

			if "`replace'"=="" preserve

			dropobs `lhs' `inexog' `endo' `exexog', fe(`fe') cluster(`cluster')

			demean `lhs' `inexog' `endo' `exexog', fe(`fe') cluster(`cluster')

		}

	*************************************
	** AREG Sample Restrict            **
	*************************************

		if "`areg'"!="" {

			preserve

			dropobs `lhs' `inexog' `endo' `exexog', fe(`fe') cluster(`cluster')

		}
	*************************************
	** Options Logic                   **
	*************************************

		if "`cluster'"!="" local cluster = "cluster(`cluster')"
		if "`cluster'"!="" local cluster2 = ",`cluster'"

		if "`robust'"!="" local robust = "`robust'"


	**************************************
	** Memory Saving FE Reg             **
	**************************************

		if "`areg'"!="" local adj = "a"
		if "`areg'"!="" local adj2 = "d"
		if "`areg'"!="" local absorb1 = ",absorb(`fe')"
		if "`areg'"!="" local absorb2 = "absorb(`fe')"

	***************************************
	** First Stage                       **
	***************************************

		* Need a first stage for each endogenous variable. This will create predicted variables, the names of which will be stored in the local endogenous_hat.
		* Also, if we are saving first stage results, need to create local variables to store results of the first stage F test


		local endog_count = 0

		local endogenous_hat ""

		foreach endog_var in `endo' {

			local endog_count = `endog_count' + 1

			* The 'adj' is for the areg version, as is 'absorb2'

				* If first option is used, display first stage regs

			local quiet = "quietly"
			if "`first'"!="" local quiet = ""
			`quiet' `adj'reg `endog_var' `exexog' `inexog', level(95) `absorb2' `cluster'

			if "`first'"!="" {

				test `exexog'
				global ftest_`endog_count' = r(F)


			}

			* The 'adj2' is also for the areg version
			predict double `endog_var'hat if e(sample), xb`adj2'

			* Create a string of all of the 'hat' variables

			local endogenous_hat = "`endogenous_hat' `endog_var'hat"


		}

		* Temporarily renaming the endogenous_hat variables. This is only because we want the names in the results matrix to match up with the original names.


		foreach endog_var in `endo' {

			rename `endog_var' `endog_var'real
			rename `endog_var'hat `endog_var'


		}

	***************************************
	** Second Stage                      **
	***************************************

		quietly `adj'reg `lhs' `endo' `inexog' , level(95) `absorb2'

		local sst_t = e(mss) + e(rss)

			* Areg FE Coefficients
			if "`areg'"!="" predict feval if e(sample), d

		***************************
		** Store E Class Returns **
		***************************

			scalar e_N=e(N)
*			scalar e_df_m = e(df_m)
*			scalar e_df_r = e(df_r)
*			scalar e_F = e(F)
*			scalar e_r2 = e(r2)
*			scalar e_rmse = e(rmse)
*			scalar e_mss = e(mss)
*			scalar e_rss = e(rss)
*			scalar e_r2_a = e(r2_a)
*			scalar e_ll = e(ll)
*			scalar e_ll_0 = e(ll_0)

		************************
		** Store COEF Results **
		************************

			matrix bmatrix = e(b)

		************************************
		** Save DOF From Second Stage Reg **
		************************************

			local dof = e(df_r)

	******************************
	** Adjust Standard Errors   **
	******************************

		***************
		** Step 1    **
		***************

			* Calculating the correct residuals, using the betas calculated in the second stage regression, but the using the "real" endogenous variables instead of the endogenous hat variables. As a trick to do this, we simply rename the "real" endogenous variables and then use the Stata predict command

			foreach endog_var in `endo' {

				rename `endog_var' `endog_var'hold
				rename `endog_var'real `endog_var'


			}

			* For AREG regression we can't get the residual directly -- we must calculate it by hand (y - xB). Here we get xb.

			if "`areg'"=="" predict double res, residual
			if "`areg'"!="" predict double pre1, xb


			* We now undue the renaming.

			foreach endog_var in `endo' {

				rename `endog_var' `endog_var'real
				rename  `endog_var'hold `endog_var'

			}

			* Areg Residuals
			if "`areg'"!="" gen res = `lhs' - pre1 - feval

		**************************
		** Step 2               **
		**************************

			* Trick to get clustered standard errors - Regress the residuals from the second stage on the endogenous hat variables


			* Degrees of Freedom adjustment - only necessary if there areg or there is no clustering
			local dof_alt = ""
			local dof_final = .
			local distinct = .

				if "`cluster'"=="" | "`areg'"!="" {

					capture distinct `fe'

					* The degree of freedom will be the same degree of freedom from the second stage, minus the number of distinct values of the fixed effect, plus one
					local dof_final = `dof' - r(ndistinct) + 1
					local dof_alt = "dof(`dof_final')"
					local distinct = r(ndistinct)

				}
			if "`areg'"!="" local dof_alt ""
			quietly `adj'reg res `endo' `inexog', level(95) `cluster' `robust' `absorb2' `dof_alt'

		**************************
		** Step 3               **
		**************************

			* Final adjustments to standard errors

			local N_clust = e(N_clust)

			local rss_t = e(rss)

			* AREG Adjustments to formulas. This is necessary because areg adjusts automatically the DOF for the absorbed variables, while regressions do not. These put the areg results in line with the demeaned regression results

			if "`areg'"!="" {
				local dof_final = `dof'
				local dof = `dof_final' + `distinct' - 1
			}


			* Correct for Non Small Sample Standard Errors if neccesary
			if "`areg'"=="" | "`areg'"!=""  {

				if "`small'"!="small" {

					if "`robust'"=="" {
						if "`cluster'"!="" mat vpr=e(V)*( (`dof')*(e(N_clust) - 1) ) / ( (e(N) - 1) * e(N_clust))
						if "`cluster'"=="" mat vpr=e(V)*(`dof_final' / (`dof_final' + ( e(N)-(`dof'+1) ) ) )


					}
					else {

						if "`cluster'"!="" mat vpr=e(V)*( (`dof')*(e(N_clust) - 1) ) / ( (e(N) - 1) * e(N_clust))
						if "`cluster'"=="" mat vpr=e(V)*(`dof_final' / (`dof_final' + ( e(N)-(`dof'+1) ) ) )*(`dof'/`dof_final')

					}

				}
				else {


					if "`robust'"=="" mat vpr=e(V)
					else if "`robust'"=="robust" & "`cluster'"=="" mat vpr=e(V)*(`dof'/`dof_final')
					else mat vpr=e(V)

				}

			}
			else {


				if "`small'"!="small" {

					if "`robust'"=="" {

						mat vpr=e(V)*`dof'/(`dof'+1)

					}
					else {

					}

				}
				else {

					if "`robust'"=="" {
						mat vpr=e(V)
					}
					else {

					}
				}




			}

			* AREG and XTIVREG2 differ when clustering is involved
			if "`areg'"!="" & ("`cluster'"!="" | "`robust'"!="") mat vpr = vpr*`dof_final'/`dof'


	*******************************
	** Step 4                    **
	*******************************

		* Display Output


			ereturn post bmatrix vpr, depname(`lhs')

			ereturn scalar N = e_N
			ereturn scalar r2= 1-(`rss_t')/(`sst_t')
			ereturn scalar N_clust = (`N_clust')
			ereturn scalar mss = `sst_t' - `rss_t'
			ereturn scalar rss = `rss_t'
			ereturn scalar df_rn = `dof'
			ereturn scalar rmse = sqrt(e(rss)/e(df_rn))
			ereturn scalar yyc = `sst_t'
			ereturn scalar yy = e(yyc)
			ereturn scalar r2u = e(r2)
			ereturn scalar r2c = e(r2u)
*			ereturn scalar df_m = e_df_m
*			ereturn scalar df_r = e_df_r
*		 	ereturn scalar F= e_F
*			ereturn scalar r2_a = e_r2_a
*			ereturn scalar ll = e_ll
*			ereturn scalar ll_0 = e_ll_0
*			ereturn local title "Fixed Effects Two Stage Least Squares Regression"
*			ereturn local method "Secret Sauce"
*			ereturn local depvar "`lhs'"
*			ereturn local cmd "xt2sls"


			* Returning First Stage Information
			if "`first'"!="" {
				forvalues i = 1/`endog_count' {

					ereturn scalar fstat_`i' = ${ftest_`i'}

				}
			}


			di

			di in gr "tsls, Fixed-Effects Two-Stage Instrumental Variable Results"

			if "`cluster'"!="" di in gr "Number of clusters (`clustvar') = " _col(33) in ye %6.0f e(N_clust) _continue

			di in gr _col(55) "Number of obs = " in ye %8.0f e(N)

*			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 .
			di

			ereturn display

	***************************
	** Dropping Variables    **
	***************************

		* Unrenaming

		foreach endog_var in `endo' {

			rename `endog_var' `endog_var'hat
			rename `endog_var'real `endog_var'

		}

		foreach var in `endogenous_hat' res {

			drop `var'

		}

		if "`areg'"!="" drop feval pre1

	*********************************
	** Restore Demean/Dropped Obs  **
	*********************************

		if "`replace'"=="" & "`demean'"!=""  restore

		if "`areg'"!="" restore


end



**************************************************************************************
program define IsStop, sclass
				/* sic, must do tests one-at-a-time,
				 * 0, may be very large */
	version 8.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

**************************************************************************************
program define demean

	syntax namelist(name=0), fe(varlist) [cluster(str)]

	sort `fe'

	foreach var in `0' {

		by `fe' : egen float `var'_m = mean(`var')

		quietly replace `var'_m = `var' - `var'_m

		drop `var'

		rename `var'_m `var'

	}

end

**************************************************************************************
program define dropobs

	* This program drops observations in which there are
	syntax namelist(name=0), fe(varlist) [cluster(str)]


	foreach var in `0' {

		quietly drop if `var'>=.

	}

	* Taking care of cluster now
	if "`cluster'"!="" {
		capture drop if `cluster' >=.
		capture drop if `cluster'==""
	}

	sort `fe'

		* Get rid of singleton groups -- this is to match up with xtivreg2.

			by `fe': egen temp_count = count(`fe')

			quietly drop if temp_count==1

			drop temp_count


end

* Fin