ENH: Add Hansen-Lee misspecification-robust J-test to IVGMM

linearmodels/iv/model.py

@@ -29,7 +29,7 @@
     sqrt,
     squeeze,
 )
-from numpy.linalg import eigvalsh, inv, matrix_rank, pinv
+from numpy.linalg import eigvalsh, inv, matrix_rank, pinv, solve
 from pandas import DataFrame, Series, concat
 from scipy.optimize import minimize
 
@@ -1044,10 +1044,44 @@ def _gmm_post_estimation(
             "weight_config": self._weight.config,
             "iterations": iters,
             "j_stat": self._j_statistic(params, weight_mat),
+            "centered_j_stat": self._centered_j_statistic(params),
         }
 
         return gmm_specific
 
+    def _centered_j_statistic(
+        self, params: linearmodels.typing.data.Float64Array
+    ) -> WaldTestStatistic | InvalidTestStatistic:
+        """
+        J statistic that uses a centered estimate of the moment covariance
+
+        The covariance of the moment conditions is estimated using the same
+        estimator as the model, i.e., the one selected by ``weight_type`` and
+        ``weight_config``, except that the mean of the moment conditions is
+        always subtracted, and it is evaluated at ``params``.
+        """
+        y, x, z = self._wy, self._wx, self._wz
+        nobs, nvar, ninstr = y.shape[0], x.shape[1], z.shape[1]
+        name = "Centered J-test of overidentifying restrictions"
+        eps = y - x @ params
+        g_bar = (z * eps).mean(0)
+
+        weight_matrix_estimator = WEIGHT_MATRICES[self._weight_type]
+        centered = weight_matrix_estimator(**{**self._weight_kwargs, "center": True})
+        s_c = centered.weight_matrix(x, z, eps)
+        if matrix_rank(s_c, hermitian=True) < ninstr:
+            # The centered cluster sums add to zero, so a clustered estimator
+            # has rank at most one less than the number of clusters
+            reason = (
+                "The centered covariance of the moment conditions is singular, "
+                "so the statistic is not defined. This occurs, for example, "
+                "when there are too few clusters."
+            )
+            return InvalidTestStatistic(reason, name=name)
+        stat = float(nobs * g_bar @ solve(s_c, g_bar))
+        null = "Expected moment conditions are equal to 0"
+        return WaldTestStatistic(stat, null, ninstr - nvar, name=name)
+
     def _j_statistic(
         self,
         params: linearmodels.typing.data.Float64Array,
@@ -1337,6 +1371,7 @@ def _gmm_post_estimation(
             "weight_config": self._weight.config,
             "iterations": iters,
             "j_stat": self._j_statistic(params, weight_mat),
+            "centered_j_stat": self._centered_j_statistic(params),
         }
 
         return gmm_specific

linearmodels/iv/results.py

@@ -1589,6 +1589,7 @@ def __init__(
         self._weight_config = results["weight_config"]
         self._iterations = results["iterations"]
         self._j_stat = results["j_stat"]
+        self._centered_j_stat = results["centered_j_stat"]
 
     @property
     def weight_matrix(self) -> linearmodels.typing.data.Float64Array:
@@ -1648,6 +1649,75 @@ def j_stat(self) -> InvalidTestStatistic | WaldTestStatistic:
         """
         return self._j_stat
 
+    @property
+    def centered_j_stat(self) -> InvalidTestStatistic | WaldTestStatistic:
+        r"""
+        J-test of overidentifying restrictions using a centered moment covariance
+
+        Returns
+        -------
+        WaldTestStatistic
+            J statistic test of overidentifying restrictions that uses a
+            covariance of the moment conditions that has been centered.
+            An InvalidTestStatistic is returned if that covariance is
+            singular, which can happen if there are too few clusters.
+
+        Notes
+        -----
+        The statistic is defined as
+
+        .. math ::
+
+          n \bar{g}'\hat{S}_c^{-1}\bar{g} \sim \chi^2_q
+
+        where :math:`\bar{g} = n^{-1}\sum \hat{\epsilon}_i z_i`,
+        :math:`z_i` includes both the exogenous variables and instruments,
+        :math:`\hat{\epsilon}_i` are the model residuals and the degree of
+        freedom is :math:`q = n_{instr} - n_{endog}`. :math:`\hat{S}_c` is the
+        estimator of the variance of :math:`\sqrt{n}\bar{g}` selected by
+        ``weight_type`` and ``weight_config``, applied to the centered moment
+        conditions :math:`g_i - \bar{g}` and evaluated at the final parameter
+        estimates. For ``weight_type="robust"`` it is
+
+        .. math ::
+
+          \hat{S}_c = n^{-1}\sum_{i=1}^{n}(g_i - \bar{g})(g_i - \bar{g})'
+
+        :math:`\hat{S}_c` and the :math:`W` of :attr:`j_stat` estimate the same
+        quantity when the model is correctly specified, and so the
+        statistics have the same asymptotic :math:`\chi^2_q` null
+        distribution. If the model is misspecified, :math:`E[g_i] \neq 0`,
+        both statistics are consistent tests and grow at rate :math:`n`, but
+        the uncentered estimator also contains :math:`\bar{g}\bar{g}'`. When
+        it is evaluated using weights from the final parameter estimates, as
+        in iterated GMM, this bounds the J statistic by :math:`n`, however
+        large the violation of the moment conditions, whereas the centered
+        statistic is not bounded. For ``weight_type="robust"`` and converged
+        iterated GMM, the two are related by
+        :math:`n \bar{g}'\hat{S}_c^{-1}\bar{g} = J / (1 - J / n)`, where
+        :math:`J` is the statistic in :attr:`j_stat`. This relationship is
+        :math:`J / (1 - s J / n)` if the estimator is scaled by
+        :math:`s = n / (n - k)` because ``debiased`` is ``True``. The
+        difference matters when :math:`J / n` is not small.
+
+        Centering has no effect when ``weight_type`` is "unadjusted",
+        since that estimator uses the residual variance about its mean, and
+        the statistic then coincides with :attr:`j_stat` up to the step of
+        the estimation at which the weight matrix is evaluated.
+        It also has no effect in :class:`~linearmodels.iv.model.IVGMMCUE`
+        models unless ``center=False`` was set.
+
+        References
+        ----------
+        Hall, A. R. (2000). Covariance matrix estimation and the power of the
+        overidentifying restrictions test. *Econometrica*, 68(6), 1517-1528.
+
+        Hansen, B. E. and Lee, S. (2021). Inference for iterated GMM under
+        misspecification. *Econometrica*, 89(3), 1419-1447. Theorem 1 and
+        Example 1 discuss the centered weight matrix.
+        """
+        return self._centered_j_stat
+
     def c_stat(self, variables: list[str] | str | None = None) -> WaldTestStatistic:
         r"""
         C-test of endogeneity

linearmodels/tests/iv/results/centered-j-package-check.R

@@ -0,0 +1,225 @@
+# Cross-check of centered-j-reference.csv against the R packages gmm and AER.
+#
+# Run from this directory:   Rscript centered-j-package-check.R
+# Requires:                  gmm (tested with 1.9.1) and AER (tested with 1.2.17)
+# Output:                    a table of the agreement, and an error if the
+#                            values in the file do not agree with the packages
+#
+# centered-j-reference.R, which produces the reference values, is written in
+# base R so that nothing is shared with linearmodels. This script is not part
+# of the tests and does not write anything. It checks that the reference values
+# are what established R implementations compute, to the precision that they
+# have, which is limited by their numerical optimizers and not by the
+# statistics.
+#
+#   * gmm::gmm with type="iterative" is the iterated GMM estimator. Its
+#     centeredVcov argument selects whether the covariance of the moment
+#     conditions is centered, and J is n times the objective. The J of a model
+#     whose estimator is not centered is gmm with centeredVcov=FALSE. The
+#     centered statistic is gmm::evalGmm with centeredVcov=TRUE evaluated at the
+#     estimates, with the weight matrix evaluated at the same estimates.
+#   * gmm::gmm with type="cue" is the continuously updating estimator.
+#   * vcov="iid" is the covariance of the moment conditions without any
+#     assumption on their variance, i.e., the heteroskedasticity robust
+#     estimator (evalGmm does not have vcov="MDS"), and vcov="HAC" with
+#     kernel="Bartlett" and prewhite=0 the kernel estimator. The
+#     bandwidth bw of gmm is the number of lags plus one, so that bw=5 is
+#     the 4 lags of linearmodels with weights 1 - j / 5.
+#   * AER::ivreg Sargan test is the J statistic with the unadjusted estimator
+#     of the covariance of the moments. It is the 2SLS J of the model, which
+#     is the iterated GMM estimator with that weight.
+#
+# What is not checked is what the packages do not offer: gmm does not have a
+# one-step weight matrix that is the same as linearmodels' (so the two-step
+# estimates are not compared, their first step is identity weighting), clustered
+# covariances, or the small-sample scale of debiased estimators. Those are
+# checked against theory by centered-j-reference.R.
+
+suppressMessages({
+  library(foreign)
+  library(gmm)
+  library(AER)
+})
+
+ref <- read.csv("centered-j-reference.csv", stringsAsFactors = FALSE)
+sim <- read.dta("simulated-data.dta")
+sim$const <- 1
+hs <- read.csv("housing.csv")
+hs$const <- 1
+reg <- model.matrix(~ region, hs)[, -1]
+hs <- cbind(hs, reg)
+mis <- read.csv("misspecified-data.csv")
+mis$const <- 1
+
+words <- function(s) strsplit(s, " ", fixed = TRUE)[[1]]
+
+# Data of a row of the reference table, weighted as linearmodels weights it
+model_data <- function(r) {
+  d <- switch(r$dataset, sim = sim, housing = hs, mis = mis)
+  sw <- if (r$weighted) sqrt(d$weights / mean(d$weights)) else rep(1, nrow(d))
+  exog <- words(r$exog)
+  endog <- words(r$endog)
+  instr <- words(r$instruments)
+  list(y = d[[r$dependent]] * sw,
+       X = as.matrix(d[, c(exog, endog), drop = FALSE]) * sw,
+       Z = as.matrix(d[, c(exog, instr), drop = FALSE]) * sw)
+}
+
+# The moment function of gmm, with the data in a single matrix [y, X, Z]
+moments <- function(k) {
+  function(theta, dat) {
+    e <- dat[, 1] - dat[, 2:(k + 1)] %*% theta
+    dat[, -(1:(k + 1)), drop = FALSE] * drop(e)
+  }
+}
+
+gmm_args <- function(r) {
+  if (r$weight_type == "kernel") {
+    list(vcov = "HAC", kernel = "Bartlett", bw = as.integer(r$bandwidth) + 1,
+         prewhite = 0)
+  } else {
+    list(vcov = "iid")
+  }
+}
+
+# Fit by gmm, returning J and the estimates
+fit_gmm <- function(r, md, centered, type) {
+  k <- ncol(md$X)
+  dat <- cbind(md$y, md$X, md$Z)
+  n <- nrow(dat)
+  g <- moments(k)
+  # 2SLS starting values
+  start <- drop(solve(crossprod(md$X, md$Z) %*% solve(crossprod(md$Z)) %*%
+                        crossprod(md$Z, md$X),
+                      crossprod(md$X, md$Z) %*% solve(crossprod(md$Z)) %*%
+                        crossprod(md$Z, md$y)))
+  run <- function(t0, type) {
+    args <- c(list(g = g, x = dat, t0 = t0, type = type, wmatrix = "optimal",
+                   centeredVcov = centered, crit = 1e-12),
+              gmm_args(r))
+    if (type == "iterative") args$itermax <- 1000
+    do.call(gmm, args)
+  }
+  fit <- run(start, type)
+  if (type == "cue") {
+    # The CUE objective is hard to minimize. Use the best of two starts, the
+    # 2SLS estimates and the iterated GMM estimates
+    other <- run(coef(run(start, "iterative")), "cue")
+    if (other$objective < fit$objective) fit <- other
+  }
+  list(b = unname(coef(fit)), J = unname(fit$objective * n), n = n, dat = dat, g = g)
+}
+
+# Objective of gmm times n at the parameters theta, with the weight matrix
+# evaluated at theta_w, and with or without centering of the covariance of the
+# moment conditions
+eval_at <- function(r, md, theta, theta_w, centered) {
+  k <- ncol(md$X)
+  dat <- cbind(md$y, md$X, md$Z)
+  args <- c(list(g = moments(k), x = dat, t0 = theta, tetw = theta_w,
+                 wmatrix = "optimal", centeredVcov = centered), gmm_args(r))
+  unname(do.call(evalGmm, args)$objective * nrow(dat))
+}
+
+# Sargan statistic of the 2SLS fit, which is the J of the unadjusted estimator
+sargan <- function(r, md) {
+  dat <- data.frame(y = md$y)
+  xcols <- colnames(md$X)
+  zcols <- colnames(md$Z)
+  d <- cbind(dat, as.data.frame(md$X), as.data.frame(md$Z[, setdiff(zcols, xcols), drop = FALSE]))
+  exog <- words(r$exog)
+  form <- as.formula(sprintf("y ~ 0 + %s | 0 + %s",
+                             paste(xcols, collapse = " + "),
+                             paste(zcols, collapse = " + ")))
+  fit <- ivreg(form, data = d)
+  unname(summary(fit, diagnostics = TRUE)$diagnostics["Sargan", "statistic"])
+}
+
+check <- function(label, got, want, rtol) {
+  err <- abs(got - want) / abs(want)
+  if (!is.finite(err) || err > rtol) {
+    stop(sprintf("%s: %.10g (package) vs %.10g (reference), relative difference %.2e",
+                 label, got, want, err))
+  }
+  err
+}
+
+part_a <- list()
+part_b <- list()
+for (i in seq_len(nrow(ref))) {
+  r <- ref[i, ]
+  if (r$weight_type == "clustered" || r$debiased || is.na(r$centered_j)) next
+  md <- model_data(r)
+  label <- sprintf("%s [%s]", r$id, r$method)
+
+  if (r$weight_type == "unadjusted") {
+    # 2SLS is the iterated GMM estimator with this weight, and its J is Sargan's
+    if (r$method == "two_step") next
+    s <- sargan(r, md)
+    part_a[[label]] <- c(j = check(paste(label, "J vs Sargan"), s, r$j, 1e-6),
+                         centered_j = check(paste(label, "centered J vs Sargan"), s,
+                                            r$centered_j, 1e-6))
+    next
+  }
+
+  # A. The statistics at the estimates of the reference, which are the same
+  # parameters for both, so that only the statistics are compared. J uses the
+  # weight matrix that was used to estimate the parameters
+  theta <- as.numeric(words(r$params))
+  theta_w <- as.numeric(words(r$weight_params))
+  part_a[[label]] <- c(
+    j = check(paste(label, "J"), eval_at(r, md, theta, theta_w, r$center), r$j, 1e-8),
+    centered_j = check(paste(label, "centered J"),
+                       eval_at(r, md, theta, theta, TRUE), r$centered_j, 1e-8)
+  )
+
+  # B. The estimators of the package. gmm does not have a two-step estimator
+  # that starts as linearmodels does, and its CUE is only compared for the
+  # robust estimator
+  if (r$method == "two_step" || (r$method == "cue" && r$weight_type != "robust")) next
+  fit <- fit_gmm(r, md, r$center, c(iterated = "iterative", cue = "cue")[[r$method]])
+  if (r$method == "cue") {
+    # The reference value is a minimum, so that the package cannot be lower.
+    # The numerical optimizer of the package does not always find the minimum,
+    # in which case there is nothing to compare
+    if (fit$J < r$j * (1 - 1e-6)) {
+      stop(sprintf("%s: the package finds J=%.8g, below the reference minimum %.8g",
+                   label, fit$J, r$j))
+    }
+    if (fit$J > r$j * (1 + 5e-4)) {
+      cat(sprintf("note: %s: gmm stops at J=%.6g, above the minimum %.6g; not compared
+",
+                  label, fit$J, r$j))
+      next
+    }
+  }
+  part_b[[label]] <- c(
+    j = check(paste(label, "J of the estimator"), fit$J, r$j, 5e-4),
+    params = max(abs(fit$b - theta) / pmax(abs(theta), 1e-2))
+  )
+  # The objective is flat in some directions, so that the estimates of the
+  # package are less precise than its J
+  if (part_b[[label]]["params"] >= 2e-2) {
+    stop(sprintf("%s: the estimates of the package differ by %.2e (relative) from the reference",
+                 label, part_b[[label]]["params"]))
+  }
+}
+
+cat("
+A. Statistics at the parameters of the reference (relative differences)
+")
+res_a <- do.call(rbind, part_a)
+print(signif(res_a, 3))
+cat("
+B. Estimators of the package (relative differences in J and in the parameters)
+")
+res_b <- do.call(rbind, part_b)
+print(signif(res_b, 3))
+cat(sprintf(paste0("
+A: %d models, largest relative difference %.2e in J and %.2e in ",
+                   "the centered J
+B: %d models, largest relative difference %.2e in J and ",
+                   "%.2e in the parameters
+"),
+            nrow(res_a), max(res_a[, "j"]), max(res_a[, "centered_j"]),
+            nrow(res_b), max(res_b[, "j"]), max(res_b[, "params"])))