ENH: Add Cragg-Donald test of reduced rank for first-stage regression (#622)

linearmodels/iv/common.py

@@ -5,6 +5,7 @@
 
 from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic
 import linearmodels.typing.data
+from linearmodels.iv._utility import annihilate, proj
 
 
 def find_constant(x: linearmodels.typing.data.Float64Array) -> int | None:
@@ -26,6 +27,97 @@ def find_constant(x: linearmodels.typing.data.Float64Array) -> int | None:
         return None
 
 
+def cragg_donald(
+    endog: linearmodels.typing.data.Float64Array,
+    instr: linearmodels.typing.data.Float64Array,
+    exog: linearmodels.typing.data.Float64Array,
+) -> WaldTestStatistic | InvalidTestStatistic:
+    r"""
+    Cragg-Donald test of reduced rank for the first-stage regression
+
+    Parameters
+    ----------
+    endog : ndarray
+        Weighted endogenous regressor array (nobs, nendog)
+    instr : ndarray
+        Weighted instrument array (nobs, ninstr)
+    exog : ndarray
+        Weighted exogenous regressor array (nobs, nexog), partialled out
+        before testing. Include a constant column here if the model has one.
+
+    Returns
+    -------
+    WaldTestStatistic
+        Test statistic, distributed chi2(ninstr - nendog + 1) under the
+        null that the first-stage coefficient matrix does not have full
+        column rank. Returns an InvalidTestStatistic if there are fewer
+        instruments than endogenous regressors.
+
+    Notes
+    -----
+    Let :math:`X = Z\Pi + V`, where :math:`\Pi \in \mathbb{R}^{k \times m}`.
+    The null hypothesis is :math:`\mathrm{rank}(\Pi) < m`. The test
+    statistic is
+
+    .. math::
+
+        \mathrm{CD} = (n - k - c) \cdot \lambda_{\min}\left(
+        (X^T M_Z X)^{-1} X^T P_Z X\right)
+
+    where :math:`P_Z` is the projection onto the column space of the
+    (control-partialled) instruments, :math:`M_Z = I - P_Z`, and
+    :math:`\lambda_{\min}` is the smallest eigenvalue. With a single
+    endogenous regressor, this statistic reduces to the standard
+    first-stage F-statistic.
+
+    The reported p-value uses the asymptotic chi2 distribution (Anderson
+    1951) rather than Stock-Yogo (2005) finite-sample critical values,
+    which require choosing a tolerance for maximal size distortion or
+    worst-case bias and so are not implemented here.
+
+    References
+    ----------
+    .. [1] Cragg, J. G., & Donald, S. G. (1993). Testing identifiability
+       and specification in instrumental variable models. Econometric
+       Theory, 9(2), 222-240.
+    .. [2] Anderson, T. W. (1951). Estimating linear restrictions on
+       regression coefficients for multivariate normal distributions.
+       Annals of Mathematical Statistics, 22(3), 327-351.
+    """
+    import scipy.linalg
+
+    n, k = instr.shape
+    m = endog.shape[1]
+    name = "Cragg-Donald Test"
+    if m == 0:
+        return InvalidTestStatistic(
+            "Model contains no endogenous regressors; the Cragg-Donald "
+            "statistic is not defined.",
+            name=name,
+        )
+    if k < m:
+        return InvalidTestStatistic(
+            "Number of instruments is less than the number of endogenous "
+            "regressors; the Cragg-Donald statistic is not defined.",
+            name=name,
+        )
+
+    n_controls = exog.shape[1]
+    x = annihilate(endog, exog) if n_controls > 0 else endog
+    z = annihilate(instr, exog) if n_controls > 0 else instr
+
+    x_proj = proj(x, z)
+    signal = x.T @ x_proj
+    noise = x.T @ (x - x_proj)
+
+    lambda_min = scipy.linalg.eigh(signal, noise, eigvals_only=True)[0]
+    statistic = (n - k - n_controls) * lambda_min
+
+    df = k - m + 1
+    null = "Instruments jointly identify the endogenous regressors (full rank)"
+    return WaldTestStatistic(statistic, null, df, name=name)
+
+
 def f_statistic(
     params: linearmodels.typing.data.Float64Array,
     cov: linearmodels.typing.data.Float64Array,

linearmodels/iv/results.py

@@ -675,6 +675,34 @@ def __init__(
         self._cov_type = cov_type
         self._cov_config = cov_config
 
+    @cached_property
+    def cragg_donald(self) -> WaldTestStatistic | InvalidTestStatistic:
+        """
+        Cragg-Donald test of reduced rank for the joint first-stage regression
+
+        Returns
+        -------
+        WaldTestStatistic
+            Test statistic for the null that the instruments do not
+            jointly identify all endogenous regressors. See
+            linearmodels.iv.common.cragg_donald for details.
+
+        Notes
+        -----
+        Unlike the per-variable F-statistics in ``diagnostics``, this test
+        accounts for correlation among the fitted values of the endogenous
+        regressors, and so correctly flags cases where instruments are
+        individually strong but cannot jointly distinguish between
+        correlated endogenous variables.
+        """
+        from linearmodels.iv.common import cragg_donald
+
+        w = sqrt(self.weights.ndarray)
+        endog = w * self.endog.ndarray.astype(float, copy=False)
+        instr = w * self.instr.ndarray.astype(float, copy=False)
+        exog = w * self.exog.ndarray.astype(float, copy=False)
+        return cragg_donald(endog, instr, exog)
+
     @cached_property
     def diagnostics(self) -> DataFrame:
         """

linearmodels/tests/iv/test_postestimation.py

@@ -155,3 +155,54 @@ def test_linear_restriction(data):
     formula_str = " = ".join(formula_dict.keys()) + " = 0"
     ts2 = res.wald_test(formula=formula_str)
     assert_allclose(ts.stat, ts2.stat)
+
+
+def test_cragg_donald(data):
+    res = IV2SLS(data.dep, data.exog, data.endog, data.instr).fit(cov_type="unadjusted")
+    cd = res.first_stage.cragg_donald
+    # Reference value computed via mlondschien/ivmodels rank_test on the
+    # same GH issue (statsmodels/linearmodels#622) example data
+    assert cd.df == data.instr.shape[1] - data.endog.shape[1] + 1
+
+
+def test_cragg_donald_degenerate():
+    rng = np.random.default_rng(0)
+    n = 200
+    z = rng.normal(size=(n, 1))  # 1 instrument
+    x = rng.normal(size=(n, 2))  # 2 endogenous -- k < m
+    exog_empty = np.empty((n, 0))
+
+    from linearmodels.iv.common import cragg_donald
+
+    result = cragg_donald(x, z, exog_empty)
+    assert np.isnan(result.pval)
+
+
+def test_cragg_donald_known_value():
+    # Reproduces the exact example from GH issue #622, validated against
+    # mlondschien/ivmodels's independent rank_test implementation
+    # (statistic=0.8939161043879634, p_value=0.6395707363012899)
+    rng = np.random.default_rng(0)
+    n = 1000
+    z = rng.normal(size=(n, 3))
+    h = rng.normal(size=(n, 3))
+    x = z @ np.ones((3, 2)) + h @ np.array([[1, 0], [0, -1], [0, 0]])
+    y = h @ np.array([1, 1, 0.1])
+
+    from linearmodels.iv.common import cragg_donald
+
+    exog_empty = np.empty((n, 0))
+    result = cragg_donald(x, z, exog_empty)
+    assert_allclose(result.stat, 0.8939161043879634, rtol=1e-6)
+    assert_allclose(result.pval, 0.6395707363012899, rtol=1e-6)
+
+
+def test_cragg_donald_no_endog():
+    from linearmodels.iv.common import cragg_donald
+
+    n = 100
+    empty_endog = np.empty((n, 0))
+    z = np.random.default_rng(0).normal(size=(n, 3))
+    exog_empty = np.empty((n, 0))
+    result = cragg_donald(empty_endog, z, exog_empty)
+    assert np.isnan(result.pval)