mata:mata clear
version 10.1
mata: mata set matastrict on
mata:
// m_mygmm2s 1.0.0 MES/CFB 11aug2008
void m_mygmm2s(string scalar yname, 
               string scalar endognames, 
               string scalar inexognames, 
               string scalar exexognames, 
               string scalar touse, 
               string scalar robust)

{
	real matrix Y, X1, X2, Z1, X, Z, QZZ, QZX, W, omega, V
	real vector cons, beta_iv, beta_gmm, e, gbar
	real scalar K, L, N, j
	
// Use st_tsrevar in case any variables use Stata's time-series operators.

	st_view(Y, ., st_tsrevar(tokens(yname)), touse)
	st_view(X1, ., st_tsrevar(tokens(endognames)), touse)
	st_view(X2, ., st_tsrevar(tokens(inexognames)), touse)
	st_view(Z1, ., st_tsrevar(tokens(exexognames)), touse)

// Our convention is that regressors are [endog   included exog]
// and instruments are                   [excluded exog  included exog]
// Constant is added by default and is the last column.
	cons = J(rows(X2), 1, 1)
	X2 = X2, cons
	X = X1, X2
	Z = Z1, X2

	K = cols(X)
	L = cols(Z)
	N = rows(Y)

	QZZ = 1/N * quadcross(Z, Z)
	QZX = 1/N * quadcross(Z, X)
		
// First step of 2-step feasible efficient GMM: IV (2SLS).  Weighting matrix
// is inv of Z'Z (or QZZ).
	W = invsym(QZZ)
	beta_iv = (invsym(X'Z * W * Z'X) * X'Z * W * Z'Y)
// By convention, Stata parameter vectors are row vectors
	beta_iv = beta_iv'
// Use first-step residuals to calculate optimal weighting matrix for 2-step FEGMM
	omega = m_myomega(beta_iv, Y, X, Z, robust)
// Second step of 2-step feasible efficient GMM: IV (2SLS).  Weighting matrix
// is inv of Z'Z (or QZZ).
	W = invsym(omega)
	beta_gmm = (invsym(X'Z * W * Z'X) * X'Z * W * Z'Y)
// By convention, Stata parameter vectors are row vectors
	beta_gmm = beta_gmm'

// Sargan-Hansen J statistic: first we calculate the second-step residuals
	e = Y - X * beta_gmm'
// Calculate gbar = 1/N * Z'*e
	gbar = 1/N * quadcross(Z, e)
	j = N * gbar' * W * gbar

// Sandwich var-cov matrix (no finite-sample correction)
// Reduces to classical var-cov matrix if Omega is not robust form.
// But 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(QZX' * W * QZX)

// Easiest way of returning results to Stata: as r-class macros.
	st_matrix("r(beta)", beta_gmm)
	st_matrix("r(V)", V)
	st_matrix("r(omega)", omega)
	st_numscalar("r(j)", j)
	st_numscalar("r(N)", N)
	st_numscalar("r(L)", L)
	st_numscalar("r(K)", K)
}
end

mata: mata mosave m_mygmm2s(), dir(PERSONAL) replace