Software 15.6.1 (Method of moments and maximum likelihood for an income distribution). The gamma distribution is a standard model for the distribution of income (Salem and Mount 1974). Fit it to per-capita disposable income in the \(50\) countries of LifeCycleSavings, first by matching the mean and variance, then by maximum likelihood. The likelihood equations for the gamma reduce to \(\lambda = \kappa /\bar y\) and the single equation \(\log \kappa - \psi (\kappa ) = \log \bar y - \overline {\log y}\), where \(\psi \) is the digamma function, which uniroot solves:

> y <- LifeCycleSavings$dpi; n <- length(y)   # income per head, 50 countries
> c(n = n, mean = mean(y), sd = sd(y))
        n      mean        sd
  50.0000 1106.7584  990.8689
> mm <- c(k = mean(y)^2 / var(y), rate = mean(y) / var(y))   # two moments
> f <- function(k) log(k) - digamma(k) - (log(mean(y)) - mean(log(y)))
> k.ml <- uniroot(f, c(0.01, 100))$root          # ML equation for the shape
> rbind(MM = mm, ML = c(k.ml, k.ml / mean(y)))
          k        rate
MM 1.247594 0.001127251
ML 1.234195 0.001115144

Reading the output. The shape parameter is \(1.25\) by the method of moments and \(1.23\) by maximum likelihood. A shape near one means a distribution close to the exponential: many poor countries and a long tail of rich ones, the standard deviation of \(991\) dollars being nearly as large as the mean of \(1107\). The two estimators use different moments, the mean and variance for MM, the mean and the mean of \(\log y\) for ML, and they agree closely, which is already some evidence that the gamma model fits.

Stata: gmm with the two moment expressions and the options winitial(identity) and onestep gives MM, and mlexp gives ML. EViews: Object \(\to \) New Object \(\to \) LogL for a user-written likelihood, or the GMM method in Quick \(\to \) Estimate Equation.

Software 15.6.2 (Two-step GMM and the \(J\) test). A general two-step GMM function, using all three moments, the mean, the mean square and the mean of \(\log y\), for two parameters, so that \(\ell - k = 1\). The standard errors use the numerical derivative of \(\bar g_n\) for \(\hat G\). The same test is then applied to annual rainfall in \(70\) US cities:

> gmm <- function(g, p0) {                        # two-step efficient GMM
+   Q <- function(p, W) {
+     if (any(p <= 0)) return(1e10)
+     gb <- colMeans(g(p)); drop(t(gb) %*% W %*% gb)
+   }
+   W1 <- diag(1 / apply(g(p0), 2, var))          # step 1: scaled identity
+   p1 <- optim(p0, Q, W = W1, control = list(reltol = 1e-14))$par
+   S <- crossprod(scale(g(p1), scale = FALSE)) / nrow(g(p1))
+   o <- optim(p1, Q, W = solve(S), control = list(reltol = 1e-14))
+   G <- sapply(1:2, function(j) {                # numerical derivative
+     h <- replace(numeric(2), j, 1e-6 * o$par[j])
+     (colMeans(g(o$par + h)) - colMeans(g(o$par - h))) / (2 * h[j])
+   })
+   m <- nrow(g(p0)); V <- solve(t(G) %*% solve(S) %*% G) / m
+   list(est = o$par, se = sqrt(diag(V)), J = m * o$value)
+ }
> gam <- function(y) function(p) cbind(y - p[1] / p[2],
+   y^2 - p[1] * (p[1] + 1) / p[2]^2, log(y) - digamma(p[1]) + log(p[2]))
> inc <- gmm(gam(y), mm)
> rbind(est = inc$est, se = inc$se)
            k         rate
est 1.2286083 0.0011359740
se  0.1326713 0.0001592832
> c(J = inc$J, p = pchisq(inc$J, 1, lower.tail = FALSE))
        J         p
0.0564122 0.8122592
> r <- as.numeric(precip)                         # rainfall, 70 US cities
> rain <- gmm(gam(r), c(mean(r)^2 / var(r), mean(r) / var(r)))
> c(J = rain$J, p = pchisq(rain$J, 1, lower.tail = FALSE))
          J           p
9.361445185 0.002215974

Reading the output. For income, efficient GMM gives \(\hat \kappa = 1.229\) (s.e. \(0.133\)), between the MM and ML values. The \(J\) statistic, \(0.056\) on one degree of freedom, is far from significant: the three moments agree, and the gamma model is not rejected. For rainfall it is \(9.36\), \(p = 0.002\): the gamma model is rejected. The reason is visible in the data: the rainfall figures are skewed to the left, with a skewness coefficient of \(-0.29\), while every gamma distribution is skewed to the right. The mean, mean square and mean log cannot all be matched by any gamma.

Stata: gmm with three moment expressions and wmatrix(robust), then estat overid. EViews: the \(J\)-statistic is printed with every GMM estimation.

Software 15.6.3 (Efficient GMM versus 2SLS under heteroscedasticity). \(2000\) samples of size \(200\) from \(y = x + u\) with three strong instruments and a disturbance whose standard deviation is proportional to \(e^{z_1}\). Each sample is estimated by 2SLS and by two-step GMM with \(\hat S = \frac 1n\sum \hat u_i^2z_iz_i'\), using the closed form of Theorem 15.3.1(i). A second experiment adds \(0.5z_3\) to the disturbance, making the third instrument invalid:

> set.seed(7)
> sim <- function(n = 200, R = 2000, bad = 0) {
+   replicate(R, {
+     Z <- matrix(rnorm(n * 3), n); v <- rnorm(n)
+     x <- drop(Z %*% c(0.5, 0.5, 0.5)) + v
+     u <- (0.5 * v + rnorm(n)) * exp(Z[, 1])     # heteroscedastic in z1
+     u <- u + bad * Z[, 3]                       # bad > 0: z3 invalid
+     y <- x + u; X <- cbind(1, x); Zc <- cbind(1, Z)
+     lin <- function(W) solve(t(X) %*% Zc %*% W %*% t(Zc) %*% X,
+                              t(X) %*% Zc %*% W %*% t(Zc) %*% y)
+     b1 <- lin(solve(crossprod(Zc)))             # 2SLS
+     e <- drop(y - X %*% b1); S <- crossprod(Zc * e) / n
+     b2 <- lin(solve(S))                         # efficient GMM
+     gb <- colMeans(Zc * drop(y - X %*% b2))
+     c(tsls = b1[2], gmm = b2[2], J = n * drop(t(gb) %*% solve(S) %*% gb))
+   })
+ }
> s0 <- sim(); s1 <- sim(bad = 0.5)
> round(rbind(mean = rowMeans(s0[1:2, ]), sd = apply(s0[1:2, ], 1, sd)), 3)
      tsls   gmm
mean 1.018 1.006
sd   0.365 0.249
> c(size.J = mean(s0[3, ] > qchisq(0.95, 2)),
+   power.J = mean(s1[3, ] > qchisq(0.95, 2)))
 size.J power.J
 0.0420  0.4595

Reading the output. Both estimators are nearly unbiased for the true value \(1\), but the standard deviation of GMM, \(0.249\), is a third smaller than that of 2SLS, \(0.365\): the efficient weight down-weights the moment conditions involving \(z_1\), which are noisy because the variance of \(u\) depends on \(z_1\). Under homoscedasticity the two would coincide. The \(J\) test, on \(\ell - k = 2\) degrees of freedom, rejects \(4.2\%\) of the time when all the instruments are valid, close to its nominal \(5\%\), and \(46\%\) of the time when \(z_3\) is invalid: a useful test, but one that misses the invalid instrument more often than not in samples of this size.

Stata: ivregress gmm y (x = z1 z2 z3), then estat overid. EViews: gmm y c x @ c z1 z2 z3, with the White weighting matrix.