Raise if mx < k. Handle mx = k separately.

CHANGELOG.rst

@@ -17,6 +17,10 @@ Changelog
   :func:`~ivmodels.tests.lagrange_multiplier.inverse_lagrange_multiplier_test` via the
   `kwargs`.
 
+- :class:`~ivmodels.models.kclass.KClass` now raises if `kappa >= 1` (as for the LIML
+  and TSLS estimators) and the number of instruments is less than the number of
+  endogenous regressors.
+
 0.4.0 - 2024-08-08
 ------------------
 

ivmodels/models/kclass.py

@@ -366,6 +366,18 @@
         n, k = Z.shape
         mx, mc = X.shape[1], C.shape[1]
 
+        if (
+            (
+                isinstance(self.kappa, str)
+                and self.kappa.lower() in {"tsls", "2sls", "liml"}
+            )
+            or (isinstance(self.kappa, (int, float)) and self.kappa >= 1)
+        ) and k < mx:
+            raise ValueError(
+                f"Need at least as many instruments (got {k}) as endogenous regressors "
+                f"(got {mx})."
+            )
+
         # Including an intercept is equivalent to replacing y <- M_1 y, X <- M_1 X,
         # C <- M_1 C, Z <- M_1 Z, where M_1 = Id - 1/n 1 1^T is the projection
         # orthogonal to np.ones(n). I.e., center the columns of all design matrices.
@@ -418,21 +430,25 @@
                 self.kappa_ = 0
             else:
                 self.fuller_alpha_ = self._fuller_alpha(self.kappa)
-                # If C!=None, we compute the ar_min as if we removed C from all design
-                # matrices. I.e., we replace Z <- M_C Z, X <- M_C X, y <- M_C y.
-                # We also exclude the columns in X coming from C.
-                X_proj_C, y_proj_C = proj(C, X[:, :mx], y)
-                # Here ar_min = lambdamin (
-                #   (X y)^T M_{[Z, C]} (X y)^{-1} (X y)^T P_{M_C Z} (X y)
-                # ).
-                # Thus X_proj <- P_[M_C Z] X = P_[Z, C] X - P_C X = X_proj - X_proj_C
-                # and X <- M_C X = X - X_proj_C. Some for y.
-                self.ar_min_ = self.ar_min(
-                    X[:, :mx] - X_proj_C,
-                    y - y_proj_C,
-                    X_proj=X_proj[:, :mx] - X_proj_C,
-                    y_proj=y_proj - y_proj_C,
-                )
+
+                if mx == k:
+                    self.ar_min_ = 0
+                else:
+                    # If C!=None, we compute the ar_min as if we removed C from all design
+                    # matrices. I.e., we replace Z <- M_C Z, X <- M_C X, y <- M_C y.
+                    # We also exclude the columns in X coming from C.
+                    X_proj_C, y_proj_C = proj(C, X[:, :mx], y)
+                    # Here ar_min = lambdamin (
+                    #   (X y)^T M_{[Z, C]} (X y)^{-1} (X y)^T P_{M_C Z} (X y)
+                    # ).
+                    # Thus X_proj <- P_[M_C Z] X = P_[Z, C] X - P_C X = X_proj - X_proj_C
+                    # and X <- M_C X = X - X_proj_C. Some for y.
+                    self.ar_min_ = self.ar_min(
+                        X[:, :mx] - X_proj_C,
+                        y - y_proj_C,
+                        X_proj=X_proj[:, :mx] - X_proj_C,
+                        y_proj=y_proj - y_proj_C,
+                    )
                 self.kappa_liml_ = 1 + self.ar_min_
                 self.kappa_ = self.kappa_liml_ - self.fuller_alpha_ / (
                     n - k - mc - self.fit_intercept

ivmodels/tests/likelihood_ratio.py

@@ -76,6 +76,12 @@
     n, k = Z.shape
     mx, mw, mc, md = X.shape[1], W.shape[1], C.shape[1], D.shape[1]
 
+    if k < mx + mw:
+        raise ValueError(
+            "The number of instruments must be at least the number of endogenous "
+            "regressors."
+        )
+
     if fit_intercept:
         C = np.hstack([np.ones((n, 1)), C])
 
@@ -91,11 +97,14 @@
 
     XWy_proj = np.hstack([X_proj, W_proj, y_proj.reshape(-1, 1)])
 
-    ar_min = _characteristic_roots(
-        a=oproj(D, XWy).T @ XWy_proj,
-        b=XWy.T @ (XWy - XWy_proj),
-        subset_by_index=[0, 0],
-    )[0]
+    if k == mx + mw:
+        ar_min = 0
+    else:
+        ar_min = _characteristic_roots(
+            a=oproj(D, XWy).T @ XWy_proj,
+            b=XWy.T @ (XWy - XWy_proj),
+            subset_by_index=[0, 0],
+        )[0]
 
     if md > 0:
         X = np.hstack([X, D])
@@ -189,13 +198,16 @@
     XWy_proj = np.concatenate([XW_proj, y_proj.reshape(-1, 1)], axis=1)
     XWy = np.concatenate([XW, y.reshape(-1, 1)], axis=1)
 
-    kappa_liml = np.real(
-        _characteristic_roots(
-            a=oproj(D, XWy).T @ XWy_proj,
-            b=XWy.T @ (XWy - XWy_proj),
-            subset_by_index=[0, 0],
-        )[0]
-    )
+    if k == mx + mw:
+        kappa_liml = 0
+    else:
+        kappa_liml = np.real(
+            _characteristic_roots(
+                a=oproj(D, XWy).T @ XWy_proj,
+                b=XWy.T @ (XWy - XWy_proj),
+                subset_by_index=[0, 0],
+            )[0]
+        )
 
     dfd = n - k - mc - md - fit_intercept
     quantile = scipy.stats.chi2.ppf(1 - alpha, df=mx + md) + dfd * kappa_liml

tests/tests/test_anderson_rubin.py

@@ -70,7 +70,7 @@ def test_more_powerful_sAR_critical_value_function_integrates_to_one(k, hat_kapp
 
 @pytest.mark.parametrize("alpha", [0.1, 0.05, 0.01])
 @pytest.mark.parametrize(
-    "n, k, mx, mw", [(100, 2, 1, 0), (100, 2, 2, 1), (100, 20, 5, 5)]
+    "n, k, mx, mw", [(100, 2, 2, 0), (100, 3, 2, 1), (100, 20, 5, 5)]
 )
 def test_inverse_anderson_rubin_confidence_set_alternative_formulation(
     alpha, n, k, mx, mw