*! version 5.0 24Nov2000
*original: version 5.0 14Jul1999
*
* changelog:
*
*  for October 16th version, added ability of user to exclude first stage fixed effects using "nofsfe" stata option
*

program define areg2
version 5.0

local varlist "require existing min(1)"
local if "optional"
local in "optional"
*local weight "fweight aweight pweight iweight"
*local exp "optional noprefix"
local options "EXOG(string) ENDOG(string) IV(string) CHOICE(string) FE(string) Absorb(string) NOFSFE "
parse "`*'"

if "`exog'" ~= "" {
	unabbrev `exog'
	local exog "$S_1"
			}

if "`endog'" ~= "" {
	unabbrev `endog'
	local endog "$S_1"
			}
if "`iv'" ~= "" {
	unabbrev `iv'
	local iv "$S_1"
			}

di _newline

tempvar use
quietly mark `use' `if' `in'
quietly markout `use' `exog' `endog' `iv' `varlist' `fe'

if "`if'" == "" { local if "if `use' == 1" }
 else { local if "`if' & `use' == 1" }

if "`weight'" ~= "" { local weight " [weight `exp'] " }

parse "`endog'", parse(" ") /* making total number of endog variables */
local C 0
while "`1'" ~= "" { 
  local C = `C' + 1
  macro shift
			}

local c 0			/* a counter for the endog vars		*/
while `c' < `C' {
  local c = `c' + 1		/* incrementing the counter */

if `c' == 1 { preserve }	/* nec only once, so do it on the first erhs var */
quietly keep if `use' == 1

quietly { /* first stage */
  tempvar resid`c' eresid`c' hat`c' pred`c'
  local endog`c' : word `c' of `endog'

 if "`nofsfe'" == "nofsfe" { /* user requested that we exclude fixed effects from the first stage */
  if "`choice'" ~= "" {local nchoice = ", `choice'"} /* necessary because of , issue in regression below */
  reg `endog`c'' `exog' `iv' `weight' `if' `in' `nchoice' /* note the use of "nchoice" so we can use quote choice" below */
  pred double  `hat`c'' `if'  /* the prediction, not including fe */
  local hatvars "`hatvars' `hat`c''"
  di "hatvars is `hatvars'"
 }

 if "`nofsfe'" == "" { /* user did not request we exclude fixed effects from the first stage */
  areg `endog`c'' `exog' `iv' `weight' `if' `in', a(`fe') `choice'
  pred double  `pred`c'' `if'  /* the prediction, not including fe */
  gen double  `resid`c'' = `endog`c'' - `pred`c'' `if' /* residuals including fe */
  sort `fe'
  egen double `eresid`c'' = mean(`resid`c''), by(`fe')  /* the estimated fe */
  gen double `hat`c'' = `pred`c'' - `eresid`c'' /* this is the fitted val of endog c */
  local hatvars "`hatvars' `hat`c''"
  di "hatvars is `hatvars'"
 }
} /*ends quietly block for first stage*/
} /*ends while loop for endogenous variables*/

quietly { /* second stage */
    areg `varlist' `exog' `hatvars' `weight' `if' `in', a(`fe') `choice'
    local aregse = $S_E_sse
    local df = $S_E_tdf
    local numobs = $S_E_nobs
	tempname b V
    matrix `b' = get(_b)
    matrix `V' = get(VCE)  
    matrix dir
    matrix list `b'

	/*
	   in case you're wondering: changing the names assoc with `b' and `v' is critical:
	   it gets us correct ses b/c it means in predicting out the sum of sq res below,
	   we are using the *actual*, not fitted values, of endog vars. the fact that initial
	   results are computed using fitted values is irrelevant, b/c we multiply by ratio
	   of correct st. err. of reg to wrong st err of reg, so the wrong one disappears.
	*/

    matrix colnames `b' = `exog' `endog' _cons
    matrix colnames `V' = `exog' `endog' _cons
    matrix rownames `V' = `exog' `endog' _cons
    matrix post `b' `V', depname(`varlist')
    matrix `b' = get(_b)
    matrix `V' = get(VCE)  
    matrix list `b'
		}

  di "2SLS using areg in each stage" _column(56) "Number of obs =   `numobs'"
  di "Fixed effects included for variable:" _column(38) "`fe'"
  di "Endogenous variables:" _column(38) "`endog'"
  di "Instrumental variables:" _column(38) "`iv'"
  di "Exogenous rhs variables:" _column(38) "`exog'" _newline
  di "Note: constant is meaningless."

quietly {
    tempvar res2 pred2 resid2 eresid2 sqres
    pred double  `pred2' `if' `in'  /* the prediction, not including fe */
    gen double  `resid2' = `varlist' - `pred2' `if' `in' /* the residuals including fe */
    sort `fe'
    egen double  `eresid2' = mean(`resid2'), by(`fe')  /* the estimated fe */
    gen double `res2' = `resid2' - `eresid2' `if' `in' /* the residuals, not including the fe */
    gen double `sqres' = `res2'*`res2' /* squared residual */
    sum `sqres' `if' `in', detail
    local s = _result(18) 	/* sum of squared errors div by model df */
    local r = `s' / `aregse' /* ratio of correct to incorrect sse */
}

restore 

 matrix `V' = `V'*`r'
 matrix colnames `V' = `exog' `endog' _cons
 matrix rownames `V' = `exog' `endog' _cons
 matrix post `b' `V', depname("`varlist'") dof(`df')
 matrix mlout
end