capture program drop locmtest4
program locmtest, rclass
	version 13

	
	_iv_parse `0'
	
	local y `s(lhs)'
	local s `s(endog)'
	local x `s(exog)'
	local z `s(inst)'
	local 0 `s(zero)'
	
	syntax [if] [, GRaph COEFficients]
	

	clear mata

	
	preserve   		

	display " "
	display " =================================================== "
	display "  Output for the Lochner & Moretti (2015) Wald Test  "
	display " =================================================== "
	display " "
	display " Output Variable y: " "`y'"
	display " Endogenous Variable s: " "`s'"
	display " Instruments z: " "`z'"
	display " "


	****************************
	* Eliminate Missing Values *
	****************************
	
	marksample touse
	qui reg `y' `s' `z' `x' if `touse' 
	quietly keep if e(sample)
	local n = e(N)
	display " Number of observations = " `n'


	*******************************************************
	* Hausman Test for Endogeneity (Auxiliary Regression) *
	*******************************************************


															
	* Step 1: Regress Endogenous variable on Exogenous variables
	quietly xi: reg `s' `z' `x'
	
	* Step 2: Obtain residuals
	quietly predict s_res, residuals


	* Step 3: Regress including Residuals
	quietly xi: reg `y' `s' s_res `x'

	* Step 4: Get the F-Statistic 	
	quietly test s_res 

	local HAux = r(F)
	local pHAux = r(p)

	drop s_res  

	************************
	* Lochner-Moretti Test *
	************************
	
	********************
	* Generate Dummies *
	********************
	quietly tab `s',  matrow(CAT)
	local ns = rowsof(CAT)
	display " Number of Categories of Endogenous variable is: " "`ns'"

	
	local nd = `ns'-1
	display " Number of Dummies is: " "`nd'"

	mata: C = J(`nd',1,.)
	
	forvalues k=2/`ns'{
		local a = CAT[`k',1]
		generate D`a' = (`s'>=`a')
		mata: C[`k'-1,1] = `a'
	}

	* Number of Regressors *
	local nr =0	

	* Get the Number of Excluded Instruments *
	
	local ni = 0
	foreach var in `z'{
  		local ni = `ni' + 1
  	}

	display " "
	display " The number of Excluded Instruments is: " `ni'
	*****************
	* Project Out x *
	*****************

	foreach var in `y' `z'{
		quietly   xi: reg `var' `x'

		predict `var'tilde, residuals
		quietly keep if e(sample)
	}
	
	******************
	* Build Matrix Z *
	******************
	mata: Z = J(`n',`ni',0)
	
	foreach var in `z'{
		mata: st_view(`var'tilde=., ., "`var'tilde")
	}
	
	local pi = 0
	foreach var in `z'{
		local pi = `pi' + 1
		mata: Z[.,`pi'] = `var'tilde
		mata: mata drop `var'tilde
	}

	
	foreach var in `z'{
	drop `var'tilde
	}

	mata: st_view(y = ., ., "`y'tilde")

	quietly   xi: reg `s' `x'
		local nx = e(rank)

	predict `s'tilde, residuals
	quietly keep if e(sample)
	
	forvalues k=2/`ns'{
		local a = CAT[`k',1]

		quietly   xi: reg D`a' `x'
  
		predict Dtilde`a', residuals
	}
	*****************************
	* Build Matrix of Dummies D *
	*****************************
	quietly tab `s',  matrow(CAT)

	forvalues k = 2/`ns'{
		local a = CAT[`k',1]
	      mata: st_view(D`a' = ., ., "Dtilde`a'")
      }

	mata: D = J(`n',`nd',0)

	forvalues k=1/`nd'{
		local l = `k' + 1
		local a = CAT[`l',1]
		mata: D[.,`k'] = D`a'		
		mata: mata drop D`a'
	}

	drop Dtilde*
	***********************************
	* Build Matrices S, y, X1, X2, Z2 *
	***********************************
	mata: st_view(S = ., ., "`s'tilde")
	* Rename Matrices *
	mata: mata rename D X1 
	mata: mata rename S X2 
	mata: mata rename Z Z2 
		
	*************************************
	* Get OLS of Non-Linear Equation: B *
	*************************************
	
	mata: B1 = invsym(X1'X1)*X1'y

	mata: EPS =  y - X1*B1
	mata: V1 = (1/(`n'-`nx'-1))*(EPS'*EPS)*invsym(X1'X1)

	mata: B = J(`nd',1,0)
	mata: VB = J(`nd',1,0)

	forvalues k=1/`nd'{
		mata: B[`k'] = B1[`k']
		mata: VB[`k'] = sqrt(V1[`k',`k'])
	}



	************************
	* IV of Linear Eqn: BL *
	************************


	mata: BL1 = invsym(X2'Z2*invsym(Z2'Z2)*Z2'X2)*X2'Z2*invsym(Z2'Z2)*Z2'y

	mata: v =  y - X2*BL1
	
	mata: BL = BL1[1]


	*********************
	* OLS of Linear Eqn *
	*********************

	mata: BL1ols = invsym(X2'X2)*X2'y

	mata: BLols = BL1ols[1]
	mata: EPSL = y - X2*BL1ols
	mata: sigmaepsL = (1/`n')*(EPSL'EPSL)
	mata: SIGMAEPSL = sigmaepsL*invsym(X2'X2)
	mata: SD_BLols = sqrt(SIGMAEPSL[1,1])


	** Drop Unnecessary Matrices **

	mata: mata drop BL1ols EPSL sigmaepsL SIGMAEPSL y

	*************************
	* Computing OLS Weights *
	*************************
	mata: Wols = J(`nd',1,0)
	mata: VWols = J(`nd',1,0)
	
	forvalues k = 1/`nd'{
		local l = `k' + 1
		mata: Wols[`k'] = invsym(X2'X2)*X2'X1[.,`k']
		mata: VWols[`k'] = sqrt((((X1[.,`k']-X2*Wols[`k'])'*(X1[.,`k']-X2*Wols[`k']))/(`n'-`nx' -1))*invsym(X2'X2))
	}


	**************************
	* Computing 2SLS Weights *
	**************************

	mata: W = J(`nd',1,0)
	mata: VW = J(`nd',1,0)

	forvalues k=1/`nd'{
		local l = `k' + 1
		local a = CAT[`l',1]
		mata: W[`k'] = invsym(X2'Z2*invsym(Z2'Z2)*Z2'X2)*X2'Z2*invsym(Z2'Z2)*Z2'X1[.,`k']
		mata: VW[`k'] = sqrt((((X1[.,`k']-X2*W[`k'])'*(X1[.,`k']-X2*W[`k']))/(`n'-`nx' - 1))*invsym(X2'Z2*invsym(Z2'Z2)*Z2'X2))
		mata: PSI`a'= X1[.,`k'] - X2*W[`k']
	}
	
	
	***********************
	* Vector of Estimates *
	*    ( Without X )    *
	***********************

	mata: THETA = (B1 \ BL1 \ W)


	******************************
	* Build Matrices A, UZ and V *
	******************************

	mata: UZ = X1:*EPS
	mata: A11 = invsym(X1'X1)

	* Release Memory Space *
	mata: mata drop X1 EPS

	
	mata: UZ = (UZ, Z2:*v)
	mata: mata drop v

	forvalues k=1/`nd'{
		local l = `k' + 1
		local a = CAT[`l',1]
		mata: UZ = (UZ, Z2:*PSI`a')
		mata: mata drop PSI`a'
	}
	
	mata: CAPLAMBDA = UZ'UZ
	mata: mata drop UZ


	mata: A22T = invsym(X2'Z2*invsym(Z2'Z2)*Z2'X2)*X2'Z2*invsym(Z2'Z2)
	mata: mata drop X2 Z2

	mata: A22 = I(`nd'+1)#A22T

	mata: A = blockdiag(A11,A22)
	mata: mata drop A11 A22 A22T

	mata: V = A*CAPLAMBDA*A'
	mata: mata drop A CAPLAMBDA


	*********************************
	* Standard Errors of Estimators *
	*********************************

	mata: SD = sqrt(diagonal(V))
	
	mata: SDadj = sqrt(diagonal(V))*sqrt(`n'/(`n'-(1+`nr')))


	******************************	
	* Compute the Test Statistic *
	******************************

	mata: T = BL - W'B
	
	mata: OX = J(`nr',1,0)
	
	mata: Q = (B, J(`nd',1,1)#OX')

	mata: G = (-W \ OX \ 1 \ OX)

	forvalues k=1/`nd'{
		mata: G = (G \ -Q[`k',.]')
	}

	mata: WN = (T^2)/(G'V*G)
	

	*********************************************
	********** O U P U T   T A B L E S **********
	*********************************************

	*********************************
	* Getting Results for the Table *
	*********************************

	* Std Error of IV *
	local sdiv = `nd' + `nr' + 1
	mata: SD_BLiv = SD[`sdiv',1]

	* Std Error of Reweigthed OLS *
	mata: G2 = (W \ OX \ 0 \ OX)
	forvalues k=1/`nd'{
		mata: G2 = (G2 \ Q[`k',.]')
	}
	mata: SD_WBols = sqrt(G2'V*G2)

	* Difference BLiv - WBols *
	mata: SD_T = sqrt(G'V*G)

	* p-value Lochner & Moretti Wald Test *
	mata: st_numscalar("r(wm)", WN)
	local pWN = chi2tail(1,r(wm))
	mata: WNpvalue = `pWN'

	* Durbin-Wu-Hausman Test *
	mata: pHA = `pHAux'
	mata: HA = `HAux'

	*******************************
	* [(IV-OLS)/StdDev(IV-OLS)]^2 *
	*******************************

	mata: NaiveTest = ((BL-BLols)/(SD_BLiv-SD_BLols))^2
	mata: st_numscalar("r(NT)",NaiveTest)
	local NT = r(NT)
	local pNT =	chi2tail(1,`NT')
	mata: NTpvalue = `pNT'

	*** TABLES ***
	
	display " "
	display " Estimated Coefficients "
	display " "
	mata: printf("{txt}{space 13}{c |}      Coef.    Std. Err.\n")
	mata: printf("{hline 13}{c +}{hline 24}\n")
	mata: printf("{txt}%12s {c |} {res}%10.0g  %10.0g\n", "OLS", BLols, SD_BLols)
	mata: printf("{txt}%12s {c |} {res}%10.0g  %10.0g\n", "IV", BL, SD_BLiv)
	mata: printf("{txt}%12s {c |} {res}%10.0g  %10.0g\n", "RWOLS", B'W, SD_WBols)
	mata: printf("{hline 13}{c +}{hline 24}\n")
	
	display " "
	display " Estimated Test Statistics "
	mata: printf("{txt}{space 13}{c |}      Test     p-value\n")
	mata: printf("{hline 13}{c +}{hline 24}\n")
	mata: printf("{txt}%12s {c |} {res}%10.0g  %10.0g\n", "LM-Wald", WN, WNpvalue)
	mata: printf("{txt}%12s {c |} {res}%10.0g  %10.0g\n", "Naive Wald", NaiveTest, NTpvalue)
	mata: printf("{txt}%12s {c |} {res}%10.0g  %10.0g\n", "DWH Test", HA, pHA)	
	mata: printf("{hline 13}{c +}{hline 24}\n")

	display " "
	display " NOTES: "
	display " "
	display " RWOLS = Reweighted OLS using TSLS Weights "
	display " "
	display " LM-Wald = Lochner-Moretti Wald Test"
	display " "
	display " Naive Wald = [ (IV - OLS) / SD(IV-OLS) ]^2 "	
	display " "
	display " DWH Test: Durbin-Wu-Hausman Test (Augmented Regression) "


clear results
	*****************
	* Saved Results *
	*****************	
	
	* Matrix of regression Coefficients *
	*************************************
	mata: st_matrix("e(B)", B)
	mata: st_matrix("e(VB)",VB)

	mata: st_matrix("e(W)", W)
	mata: st_matrix("e(VW)", VW)

	mata: st_matrix("e(Wols)", Wols)
	mata: st_matrix("e(VWols)", VWols)
	
	* OLS of Linear Eqn *
	*********************
	mata: st_numscalar("e(BLols)",BLols)
	mata: st_numscalar("e(SDBLols)",SD_BLols)

	* 2SLS of Linear Eqn *
	**********************
	mata: st_numscalar("e(BLiv)",BL)
	mata: st_numscalar("e(SDBLiv)",SD_BLiv)

	* Difference BLiv - BLols *
	***************************
	mata: st_numscalar("e(DIVOLS)",BL-BLols)
	mata: st_numscalar("e(SDDIVOLS)",SD_BLiv-SD_BLols)

	* Reweigthed OLS *
	******************	
	mata: st_numscalar("e(WBols)",B'W)
	mata: st_numscalar("e(SDWBols)",SD_WBols)	

	* Difference BLiv - WBols *
	***************************
	mata: st_numscalar("e(T)",T)
	mata: st_numscalar("e(SDT)",SD_T)

	* LM Wald Test *
	****************
	mata: st_numscalar("e(wm)", WN)
	mata: st_numscalar("e(pwm)", WNpvalue)
	
	* Naive Wald *
	**************
	mata: st_numscalar("e(nw)", NaiveTest)
	mata: st_numscalar("e(pnw)", NTpvalue)
	
	* DWH Test *
	************
	mata: st_numscalar("e(dwh)", HA)
	mata: st_numscalar("e(pdwh)", pHA)
	

	**************************
	* Print out Coefficients *
	**************************
if ("`coefficients'" == "coefficients"){
	display " "
	display " Estimated Coefficients: "
	mata: RES1 = (B,VB, W, VW, Wols, VWols)
	mata: st_matrix ("RES", RES1)
	local matrows = ""
	forvalues k=2/`ns'{
		local a = CAT[`k',1]
		local matrows = "`matrows'" + " " + "`a' "
	}
	matrix colnames RES = "B" "seB" "W2SLS" "seW2sls" "Wols" "seWols"
	matrix rownames RES = `matrows'
	matlist RES , aligncolnames(c) 

}
	**********************
	*    Create Graph    *
	**********************	
if ("`graph'" == "graph") {
qui{

local xvar : variable label `s'

**Convert matrix to stata variables
getmata B, replace force
getmata W, replace force
getmata Wols, replace force
getmata C, replace force

}
twoway (bar B C, yaxis(1) yscale(alt axis(1)) color(gs10))  /*
	 */ (line W C, yaxis(2) color(gs0) ) /*
	*/ (line Wols C, yaxis(2)  lpattern(-) yscale(alt axis(2)) color(gs0) ) /*
	*/ , ytitle("Weights", axis(2)) ytitle("OLS Coefficents", axis(1) ) /*
	*/ xtitle("`xvar'")	legend(label(1 "OLS Coefficients") label(2 "2SLS Weights") /*
	*/ label(3 "OLS Weights") )
}
end