program define ueve, eclass
version 10

syntax varlist(min=2 numeric) [if] [in] [fweight aweight pweight iweight/] , GRoup(varname) ///
     [ewald eve b2sls Level(integer $S_level)]

marksample touse
local depvar: word 1 of `varlist'
local regs: list varlist - depvar
tempname b V g d _group

// Drop groups with 1 observation, display warning message with number of these groups

sort `group'
qui egen `g'=count(`depvar') if `touse', by(`group')
qui egen `d'=count(`depvar') if `g'==1
if `d' !=. {
	di as error "Warning: " `d' " groups with one observation are dropped"
}
qui replace `g'=. if `g'==1
markout `touse' `g'
qui egen `_group'=group(`group') if `touse'
markout `touse' `_group'

// Create weights (fweight aweight pweight iweight) for Mata
   tempvar normwt
   if `"`exp'"' != "" {
      qui gen double `normwt' = `exp' if `touse'
      if "`weight'" == "aweight" | "`weight'" == "pweight" {
         summ `normwt' if `touse', mean
         di as text "(sum of wgt is " %12.4e r(sum) ")"
         qui replace `normwt' = r(N)*`normwt'/r(sum)
      }
      local wtexp `"[`weight' = `exp']"'
      summ `normwt' if `touse', mean
      if "`weight'" == "iweight" {
         local normN = trunc(r(sum))
      }
      else {
         local normN = r(sum)
      }
      markout `touse' `normwt'
   }
   else {
      qui gen double `normwt' = 1 if `touse'
      qui count if `touse'
      local normN = r(N)
   }

// Display chosen estimation method, assign it to local scalar

if "`ewald'" !="" {
  display _newline as text "EWALD regression"
  local emethod=1
}
else if "`eve'" !="" {
  display _newline as text "EVE regression"
  local emethod=2
}
else {
  display _newline as text "UEVE regression"
  local emethod=0
}

// execute program computing estimator and its covariance matrix in Mata

mata: UEVE("`depvar'","`regs'","`_group'","`normwt'","`touse'","`emethod'")

// return results and report statistics

ereturn clear

matrix `b' = r(b)
matrix `V' = r(V)
local vnames `regs' _cons

matname `V' `vnames'
matname `b' `vnames', c(.)

local k = colsof(`b')
local N = r(N)
local df = `N' - `k' - 1

ereturn post `b' `V', dof(`df') obs(`N') depname(`depvar') esample(`touse')
ereturn local depvar = "`depvar'"
ereturn local indepvars = "`regs'"
ereturn scalar N = r(N)
ereturn scalar G = r(G)
ereturn scalar r2 = r(R2)
ereturn local cmd "ueve"

display _col(53) as text "Number of obs    = " as result %7.0f e(N)
display _col(53) as text "Number of groups = " as result %7.0f e(G)
display _col(53) as text "R-squared        = " as result %7.4f e(r2) _newline
ereturn display, level(`level')

end

// main Mata program computing the estimator and its covariance matrix

mata

real matrix UEVE (string scalar depvar, string scalar regs, string scalar group, string scalar normwt, string scalar touse, string scalar emethod)
{
real scalar ng, N, G, K, K1
real vector cons, lg, b
real matrix gwyx, _gwyx, Xg, Pg, Mg, Sg, sg, sxx, sxy, syy, s1ng, V

st_view(gwyx=., ., (tokens(group), tokens(normwt), tokens(depvar), tokens(regs)), touse)

N=rows(gwyx)
G=gwyx[rows(gwyx),1]
K1=cols(gwyx)
K=K1-2

Sg=J(K,K,0)
sg=J(K,1,0)
sxx=J(K,K,0)
sxy=J(K,1,0)
syy=J(1,1,0)
s1ng=0
Ybar=.
Xbar=J(1,K,.)

// compute variance of the sampling errors

for (i=1; i<=G; i++) {
	_gwyx=select(gwyx, gwyx[.,1]:==i)
	Yg = _gwyx[.,3]
	ng = rows(Yg)
	cons = J(ng,1,1)
	Xg = _gwyx[.,4..K1], cons
	lg = J(ng,1,1)
	w = sqrt(_gwyx[.,2])
	Xg=w:*Xg
	Yg=w:*Yg
	lg=w:*lg

	Pg = lg*invsym(lg'*lg)*lg'
	sxx = sxx + Xg'*Pg*Xg
	sxy = sxy + Xg'*Pg*Yg
	syy = syy + Yg'*Pg*Yg
	s1ng = s1ng + 1/ng

	Pg=.
	Mg=-lg*invsym(lg'*lg)*lg'
	dMg=J(ng,1,1)-(-diagonal(Mg))
	_diag(Mg,dMg)
	Sg = Sg + (1/(ng-1))*Xg'*Mg*Xg
	sg = sg + (1/(ng-1))*Xg'*Mg*Yg

	// create matrices of averages for R2
	ybar=colsum(Yg)/colsum(w)
	Ybar=Ybar\ybar
	xbar=colsum(Xg)/colsum(w)
	Xbar=Xbar\xbar
}

Mxx=(1/G)*sxx
Mxy=(1/G)*sxy
Myy=(1/G)*syy
Shat=(1/G)*Sg
shat=(1/G)*sg

// choose estimation method

// UEVE
if (emethod=="0") {
   a=(G-K-1)/G
}
// EWALD
if (emethod=="1") {
   a=0
}
// EVE
if (emethod=="2") {
   a=1
}

// Compute estimator

b=qrinv(Mxx-a*Shat)*(Mxy-a*shat)

// Compute covariance matrix of estimator

Omega=Mxx-a*Shat
rho=Myy-b'*Omega*b
A=Mxx*(Myy-b'*Omega*b+b'*Shat*b-2*shat'b)+(shat-Shat*b)*(shat-Shat*b)'
B=(1/G)*s1ng*(Shat*(rho+b'*Shat*b-2*shat'b)+Shat*(shat-Shat*b)*(shat-Shat*b)')
V=(1/G)*qrinv(Omega)*(A+(a^2)*B)*qrinv(Omega)
_makesymmetric(V)

// Compute R-squared using group means

Ybar=Ybar[2..(G+1),.]
Xbar=Xbar[2..(G+1),.]
i = J(G,1,1)
R2= 1 - (Ybar-Xbar*b)'(Ybar-Xbar*b)/(Ybar'Ybar-(i'Ybar)^2/G)

st_matrix("r(b)",b')
st_matrix("r(V)",V)
st_numscalar("r(N)",N)
st_numscalar("r(G)",G)
st_numscalar("r(R2)",R2)

}

end