Software 16.7.1 (Pooled, fixed effects and first differences). A simulated panel of \(300\) workers over \(5\) years. Log wages depend on training hours \(x\) (true effect \(0.5\)), on a time-invariant characteristic \(z\), and on unobserved ability \(\alpha _i\), which also raises training:

> set.seed(8); N <- 300; T <- 5                 # 300 workers, 5 years
> id <- rep(1:N, each = T)
> a <- rnorm(N)                                  # unobserved ability
> x <- 0.8 * a[id] + rnorm(N * T)                # training, higher for able
> z <- rbinom(N, 1, 0.5)[id]                     # time-invariant (e.g. sex)
> y <- 1 + 0.5 * x + 0.3 * z + a[id] + rnorm(N * T)   # log wage
> dm <- function(v) v - ave(v, id)               # within transformation
> pooled <- lm(y ~ x + z)
> lsdv <- lm(y ~ x + factor(id))
> within <- lm(dm(y) ~ dm(x) - 1)
> fd <- lm(diff(y)[-seq(T, N * T, T)] ~ diff(x)[-seq(T, N * T, T)])
> c(pooled = coef(pooled)[["x"]], LSDV = coef(lsdv)[["x"]],
+   within = coef(within)[[1]], FD = coef(fd)[[2]])
   pooled      LSDV    within        FD
0.9663257 0.5174025 0.5174025 0.5328128

Reading the output. Pooled least squares estimates the return to training as \(0.97\), nearly double the truth, because training proxies for ability. The dummy-variable regression, with \(300\) worker intercepts, and the within regression on demeaned data give the identical \(0.517\), as Theorem 16.2.1 says; first differences gives \(0.533\). Both are close to \(0.5\). The index -seq(T, N * T, T) drops the differences that would straddle two workers.

Stata: xtset id year, then xtreg y x, fe (and regress D.(y x) for first differences). EViews: a Panel workfile, with Cross-section: Fixed in the Panel Options tab.

Software 16.7.2 (Random effects, Mundlak and Hausman). The variance components by the Swamy–Arora method, the random effects estimator by quasi-demeaning, and the two tests:

> s2e <- sum(resid(within)^2) / (N * T - N - 1)  # sigma_eps^2 (with df fix)
> xb <- tapply(x, id, mean); yb <- tapply(y, id, mean)
> zb <- tapply(z, id, mean)
> between <- lm(yb ~ xb + zb)
> s2a <- sum(resid(between)^2) / (N - 3) - s2e / T     # sigma_alpha^2
> theta <- 1 - sqrt(s2e / (s2e + T * s2a))
> qd <- function(v) v - theta * ave(v, id)       # quasi-demeaning
> re <- lm(qd(y) ~ I(rep(1 - theta, N * T)) + qd(x) + qd(z) - 1)
> c(s2e = s2e, s2a = s2a, theta = theta, RE = coef(re)[[2]])
      s2e       s2a     theta        RE
1.0204634 0.1943769 0.2843246 0.8187059
> mundlak <- lm(y ~ x + z + ave(x, id))          # add the unit means of x
> coef(summary(mundlak))[c("x", "ave(x, id)"), 1:3]
            Estimate Std. Error  t value
x          0.5174025 0.03080082 16.79834
ave(x, id) 0.9245214 0.04420127 20.91617
> V.fe <- sum(resid(within)^2) / (N * T - N - 1) / sum(dm(x)^2)
> V.re <- vcov(re)[2, 2]
> H <- (coef(within)[[1]] - coef(re)[[2]])^2 / (V.fe - V.re)
> c(se.FE = sqrt(V.fe), se.RE = sqrt(V.re), Hausman = H)
       se.FE        se.RE      Hausman
  0.02824600   0.02574678 672.77197414

Reading the output. The estimated \(\theta = 0.28\) is small, so random effects is close to pooled least squares, and gives \(0.82\): inconsistent, because ability is correlated with training. The estimated \(\hat \sigma _\alpha ^2 = 0.19\) is far below the true \(1\): part of the ability effect has been attributed to training in the between regression, so the variance components are themselves contaminated. Mundlak’s regression reproduces the fixed effects coefficient exactly, \(0.5174\), as Theorem 16.4.1 says, and the coefficient of the worker mean of training, \(0.92\) with \(t = 20.9\), measures the correlation between training and ability. The Hausman statistic, \(673\) on one degree of freedom, rejects random effects overwhelmingly.

Stata: xtreg y x z, re, then hausman fe re after storing both. EViews: Cross-section: Random, and View \(\to \) Fixed/Random Effects Testing \(\to \) Correlated Random Effects - Hausman Test.

Software 16.7.3 (Clustered standard errors). For the pooled regression, the conventional and cluster-robust standard errors, clustering by worker. rowsum adds the scores \(x_{it}\hat u_{it}\) within each worker, giving the cluster sums \(X_i'\hat u_i\):

> cl <- function(fit, g) {                       # cluster-robust variance
+   X <- model.matrix(fit); e <- resid(fit); B <- solve(crossprod(X))
+   M <- crossprod(rowsum(X * e, g)); B %*% M %*% B
+ }
> round(rbind(conventional = sqrt(diag(vcov(pooled))),
+             clustered = sqrt(diag(cl(pooled, id)))), 4)
             (Intercept)      x      z
conventional      0.0444 0.0251 0.0648
clustered         0.0582 0.0316 0.0883

Reading the output. The clustered standard errors are \(26\) to \(36\%\) larger than the conventional ones. The ratio is largest for \(z\), which is constant within workers (\(\rho _x = 1\)), as the Moulton factor predicts: the intraclass correlation of the pooled residuals is \(0.21\), and with \(T = 5\) the factor is \(1 + 4 \times 0.21 = 1.85\), a standard error ratio of \(1.36\), exactly the \(0.0883/0.0648 = 1.36\) observed.

Stata: regress y x z, vce(cluster id). EViews: Options \(\to \) Coefficient covariance \(\to \) Cluster robust.

Software 16.7.4 (Worker effects in productivity data). The Machines data in the recommended package nlme record productivity scores for six workers, each tested three times on each of three machines. The question is how the machines differ; the workers differ too:

> options(show.signif.stars = FALSE)
> data(Machines, package = "nlme")
> m <- data.frame(Machines); table(m$Worker)

6 2 4 1 3 5
9 9 9 9 9 9
> f0 <- lm(score ~ Machine, m)                   # ignores the workers
> f1 <- lm(score ~ Machine + Worker, m)          # worker fixed effects
> round(rbind(pooled = coef(summary(f0))[2:3, 2],
+             FE = coef(summary(f1))[2:3, 2],
+             clustered = sqrt(diag(cl(f0, m$Worker)))[2:3]), 3)
          MachineB MachineC
pooled       1.925    1.925
FE           1.054    1.054
clustered    2.210    1.406
> coef(f1)[2:3]
 MachineB  MachineC
 7.966667 13.916667
> anova(f1)
Analysis of Variance Table

Response: score
          Df  Sum Sq Mean Sq F value    Pr(>F)
Machine    2 1755.26  877.63  87.798 < 2.2e-16
Worker     5 1241.89  248.38  24.848 4.867e-12
Residuals 46  459.82   10.00
> ms <- anova(f1)[["Mean Sq"]]
> c(s2.worker = (ms[2] - ms[3]) / 9, s2.eps = ms[3],
+   ICC = (ms[2] - ms[3]) / 9 / ((ms[2] - ms[3]) / 9 + ms[3]))
 s2.worker     s2.eps        ICC
26.4869984  9.9960145  0.7260091

Reading the output. The machine effects are the same whether or not the worker effects are included, \(7.97\) and \(13.92\) points above machine A, because the design is balanced and the worker dummies are orthogonal to the machine dummies. But their standard errors fall from \(1.93\) to \(1.05\) when worker effects are included, because worker differences are removed from the residual. Most of the variation not explained by the machines is between workers: the estimated intraclass correlation is \(0.73\). In this orthogonal design random effects gives the same machine effects as fixed effects. The clustered standard errors are larger still, \(2.21\) and \(1.41\), and differ between the two machines, because they also allow for workers who are better on some machines than others, the worker-by-machine interaction; but with only six clusters they are themselves very imprecise.

Stata: xtreg score i.machine, i(worker) fe or mixed score i.machine || worker:. EViews: a panel workfile with worker as the cross-section identifier.

Software 16.7.5 (Nickell bias and the Anderson–Hsiao estimator). \(500\) samples of \(N = 500\) units with \(y_{it} = 0.5y_{i,t-1} + \alpha _i + \varepsilon _{it}\), using five periods in the within regression; Anderson–Hsiao uses \(y_{i,t-2}\) as the instrument for \(\Delta y_{i,t-1}\):

> set.seed(9)
> dyn <- function(N = 500, T = 6, rho = 0.5) {
+   a <- rnorm(N); y <- matrix(0, N, T + 50)
+   for (t in 2:(T + 50)) y[, t] <- rho * y[, t - 1] + a + rnorm(N)
+   y <- y[, -(1:50)]                            # drop the burn-in
+   Y <- y[, 2:T]; L <- y[, 1:(T - 1)]
+   fe <- sum(dm2(L) * dm2(Y)) / sum(dm2(L)^2)
+   dY <- Y[, -1] - Y[, -(T - 1)]; dL <- L[, -1] - L[, -(T - 1)]
+   zz <- L[, -(T - 1)]                          # y(t-2) instruments dy(t-1)
+   ah <- sum(zz * dY) / sum(zz * dL)
+   c(FE = fe, AH = ah)
+ }
> dm2 <- function(M) M - rowMeans(M)
> r <- replicate(500, dyn())
> round(rbind(mean = rowMeans(r), sd = apply(r, 1, sd)), 3)
        FE    AH
mean 0.167 0.491
sd   0.022 0.066
> nickell <- function(rho, T) {                # large-N bias, T periods
+   A <- 1 - (1 - rho^T) / (T * (1 - rho))
+   -(1 + rho) / (T - 1) * A / (1 - 2 * rho * A / ((1 - rho) * (T - 1)))
+ }
> c(plim.FE = 0.5 + nickell(0.5, 5), leading.term = 0.5 - 1.5 / 4)
     plim.FE leading.term
   0.1689189    0.1250000

Reading the output. Fixed effects averages \(0.167\) instead of \(0.5\), with a small standard deviation of \(0.022\): it is precisely wrong. Nickell’s formula gives the probability limit \(0.169\), almost exactly the simulated mean, while its large-\(T\) approximation, \(0.125\), is too crude for \(T = 5\). Anderson–Hsiao averages \(0.491\), close to the truth, at the cost of a standard deviation three times as large.

Stata: xtabond y or xtdpdsys y for Arellano–Bond and Blundell–Bond. EViews: Quick \(\to \) Estimate Equation, method GMM/DPD, with the dynamic panel wizard.