Why don't my cluster-robust (panel-robust) standard errors match those in Stata? (solved)

I'm trying to write my own code for cluster-robust (AKA panel-robust, AKA heteroskedasticity and serial-correlation-consistent) standard errors, so that I can make a couple of small extensions. But I can't get my results to match Stata's (in which this procedure is routine), so I am probably missing some detail. Would appreciate pointers on the code, below.

The model is $$y_{ig} = \alpha_g + X'_{ig}\beta + \epsilon_{ig}$$ The groupwise mean is the subtracted from every term above, getting rid of the $\alpha$, and leaving de-meaned group-varying variables: $$ydm_{ig} = Xdm'_{ig}\beta + \epsilon dm_{ig}$$

The variance of $\beta$ is estimated as

$$\left(\displaystyle\sum_G Xdm_g'Xdm_g\right)^{-1} \displaystyle\sum_G Xdm_g' \epsilon dm_g \epsilon dm_g' Xdm_g \left(\displaystyle\sum_G Xdm_g'Xdm_g\right)^{-1}$$

Standard errors are the square root of the diagonal of this matrix, inflated by $G/(G-1) \times (N-1)/(N-K)$.

This is all textbook econometrics. Standard in stata, and code exists to do it in R. (this, for example)

I want my code to not rely on plm, lmtest or sandwich.

The following script is supposed to implement the math above, simulating a panel dataset in which outcomes are autocorrelated and groups have different forms of heteroskedasticity. It gives the right point estimates but the standard errors are bigger than those given by stata:

rm(list=ls())
set.seed(999)
N = 1000
G = 200
x = rnorm(N)
z = rexp(N)
obs = N/G
fe = data.frame(ID=1:G,fe = rnorm(G)+5)
fe = data.frame(ID = rep(fe$ID,obs),fe = rep(fe$fe,obs))
t = rep(1,G); for (i in 2:obs) {t = c(t,rep(i,G))}
data = data.frame(y=x+z+fe$fe+x*rnorm(N,mean=0,sd=fe$fe*runif(N)), ID=fe$ID,x,z,t) write.csv(data,"testdata.csv") demean = function(var,ID){ dat = data.frame(var,ID) library(doBy) means = summaryBy(var~ID,data=dat,fun=mean) d = data.frame(var,ID) a = merge(d,means,by="ID") adm = a[,2]-a[,3] adm } xdm = demean(data$x,data$ID) ydm = demean(data$y,data$ID) zdm = demean(data$z,data$ID) mdm = lm(ydm~xdm+zdm-1) summary(mdm) e = mdm$resid

lpm = cbind(xdm,zdm)

tofu = matrix(0,ncol(lpm),ncol(lpm))
K = mdm$rank for (i in 1:G){ X = lpm[data$ID==unique(data$ID)[i],] bread = t(X)%*%X +bread r = e[data$ID==unique(data\$ID)[i]]
tofu = t(X) %*% r %*% t(r) %*% X + tofu
}
se = G/(G-1)*(N-1)/(N-K)*sqrt(diag(vcv))
summary(mdm)
se


Using that testdata.csv file, I can compare to Stata:

clear all
insheet using "testdata.csv"
xtset id t
xtreg y x z,fe cluster(id)


My R-code:

1> se
xdm        zdm
0.16946120 0.08793485


Stata's output:

             |               Robust
y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
x |   1.223739   .1992449     6.14   0.000     .8308371    1.616642
z |    .876592   .0960943     9.12   0.000     .6870981    1.066086


Mine are consistently more optimistic than Stata's. Anybody help me figure out why? I feel like it is probably some small bug, but I can't pinpoint where.

EDIT: The answer to the question is here and here. Stata doesn't run a textbook within'' estimator, rather it adds the averages back into each variable, which is what allows for the estimation of a constant term (which has the interpretation of the mean of the fixed effects, and which allows for prediction. When I alter my code to copy those procedures, my SE's are equivalent up to the 3rd or 4th decimal point.

• for reproducibility, you may want to provide the set.seed() in R – Metrics Aug 7 '13 at 11:42
• A complete side-issue, and I don't use R, but isn't there an easier way to "demean" (horrible word here)? Your customised function seems to do stuff every time it is called, but my gut is that R should allow easier ways to do it. – Nick Cox Aug 7 '13 at 11:56
• @Nick Cox: probably. I'll optimize the code once I get it functional. I started trying to make a function to demean a whole data frame, handle missing data, etc, but pared back my goals when I realized that I wans't getting the correct results in the first place. – generic_user Aug 7 '13 at 11:59

STATA computes a slight variation of the within estimator. However, it is possible to compute clustered standard errors in R that are identical to those in STATA. The following post gives you a function and the code that replicated STATA's clustered standard errors in R:

Clustered standard errors and robust standard errors