Move from eta formulation to kappa formulation.

ivmodels/kclass.py

@@ -176,8 +176,8 @@ def _fuller_alpha(self, kappa):
         else:
             return 0.0
 
-    def _eta_liml(self, X, y, Z=None, X_proj=None, y_proj=None):
-        """Compute the eta parameter of the LIML estimator.
+    def _kappa_liml(self, X, y, Z=None, X_proj=None, y_proj=None):
+        """Compute the kappa parameter of the LIML estimator.
 
         Either ``Z`` or both ``X_proj`` and ``y_proj`` must be specified.
 
@@ -196,7 +196,7 @@ def _eta_liml(self, X, y, Z=None, X_proj=None, y_proj=None):
 
         Returns
         -------
-        eta_liml: float
+        kappa_liml: float
             Smallest eigenvalue of ``((X y)^T (X y))^{-1} (X y)^T P_Z (X y)``, where
             ``P_Z`` is the projection matrix onto the subspace spanned by Z.
         """
@@ -207,8 +207,8 @@ def _eta_liml(self, X, y, Z=None, X_proj=None, y_proj=None):
 
         Xy = np.concatenate([X, y.reshape(-1, 1)], axis=1)
         Xy_proj = np.concatenate([X_proj, y_proj.reshape(-1, 1)], axis=1)
-        W = np.linalg.solve(Xy.T @ Xy, Xy.T @ Xy_proj)
-        return min(np.linalg.eigvals(W))
+        W = np.linalg.solve((Xy - Xy_proj).T @ Xy, Xy.T @ Xy_proj)
+        return 1 + min(np.linalg.eigvals(W))
 
     def _solve_normal_equations(self, X, y, X_proj, y_proj, alpha=0):
         if alpha != 0:
@@ -258,11 +258,12 @@ def fit(self, X, y, Z=None, *args, **kwargs):
             if self.kappa.lower() in {"tsls", "2sls"}:
                 self.kappa_ = 1
             elif self.kappa.lower() == "ols":
-                self.kappa = 0
+                self.kappa_ = 0
             else:
                 self.fuller_alpha_ = self._fuller_alpha(self.kappa)
-                self.eta_liml_ = self._eta_liml(X, y, X_proj=X_proj, y_proj=y_proj)
-                self.kappa_ = 1 / (1 - self.eta_liml_) - self.fuller_alpha_ / (n - q)
+                self.kappa_liml_ = self._kappa_liml(X, y, X_proj=X_proj, y_proj=y_proj)
+                self.kappa_ = self.kappa_liml_ - self.fuller_alpha_ / (n - q)
+
         else:
             self.kappa_ = self.kappa
 
@@ -307,10 +308,10 @@ class KClass(KClassMixin, GeneralizedLinearRegressor):
        (1 - \\kappa) \\| y - X \\beta \\|_2^2 + \\kappa \\|P_Z (y - X \\beta) \\|_2^2
        \\\\
        &= (X^T (\\kappa P_Z + (1 - \\kappa) \\mathrm{Id}) X)^{-1} X^T
-       (\\kappa P_Z + (1 - \\kappa) \\mathrm{Id}) X) y
+       (\\kappa P_Z + (1 - \\kappa) \\mathrm{Id}) X) y,
 
-    where :math:`P_Z` is the projection matrix onto the subspace spanned by :math:`Z`
-    and :math:`\\mathrm{Id}` is the identity matrix.
+    where :math:`P_Z = Z (Z^T Z)^{-1} Z^T` is the projection matrix onto the subspace
+    spanned by :math:`Z` and :math:`\\mathrm{Id}` is the identity matrix.
     This includes the the ordinary least-squares (OLS) estimator (:math:`\\kappa = 0`),
     the two-stage least-squares (2SLS) estimator
     (:math:`\\kappa = 1`), the limited information maximum likelihood (LIML) estimator
@@ -322,10 +323,11 @@ class KClass(KClassMixin, GeneralizedLinearRegressor):
     ----------
     kappa: float or {fuller(a), liml}
         The kappa parameter of the k-class estimator. If float, then kappa must be in
-        :math:`[0, \\hat\\kappa_\\mathrm{LIML} := 1 / (1 - \\eta_\\mathrm{LIML})
-        \\geq 1]`, where :math:`\\eta_\\mathrm{LIML}` is the smallest eigenvalue of
-        the matrix :math:`((X \\ \\ y)^T (X \\ \\ y))^{-1} (X \\ \\ y)^T P_Z (X \\ y)`
-        and :math:`P_Z` is the projection matrix onto the subspace spanned by :math:`Z`.
+        :math:`[0, \\hat\\kappa_\\mathrm{LIML}]]`, where
+        :math:`\\kappa_\\mathrm{LIML} \\geq 1` is 1 plus the smallest eigenvalue of the
+        matrix :math:`((X \\ \\ y)^T M_Z (X \\ \\ y))^{-1} (X \\ \\ y)^T P_Z (X \\ y)`.
+        and :math:`P_Z` is the projection matrix onto the subspace spanned by :math:`Z`
+        and :math:`M_Z = Id - P_Z`.
         If string, then must be one of ``"liml"``, ``"fuller"``, or ``"fuller(a)"``,
         where ``a`` is numeric. If ``kappa="liml"``, then
         :math:`\\kappa = \\hat\\kappa_\\mathrm{LIML}` is used. If ``kappa="fuller(a)"``,
@@ -361,9 +363,9 @@ class KClass(KClassMixin, GeneralizedLinearRegressor):
     fuller_alpha_: float
         If ``kappa`` is one of ``{"fuller", "fuller(a)", "liml"}`` for some numeric
         value ``a``, the alpha parameter of the Fuller estimator.
-    eta_liml_: float
+    kappa_liml_: float
         If ``kappa`` is one of ``{"fuller", "fuller(a)", "liml"}`` for some numeric
-        value ``a``, the eta parameter of the LIML estimator.
+        value ``a``, the kappa parameter of the LIML estimator.
 
     References
     ----------

ivmodels/tests.py

@@ -621,8 +621,8 @@ def bounded_inverse_anderson_rubin(Z, X):
 
     X_proj = proj(Z, X)
 
-    W = np.linalg.solve(X.T @ X, X.T @ X_proj)
-    eta_min = min(np.real(np.linalg.eigvals(W)))
+    W = np.linalg.solve(X.T @ (X - X_proj), X.T @ X_proj)
+    kappa = min(np.real(np.linalg.eigvals(W)))
 
-    cdf = scipy.stats.f.cdf((n - q) / q * eta_min / (1 - eta_min), dfn=q, dfd=n - q)
+    cdf = scipy.stats.f.cdf((n - q) / q * kappa, dfn=q, dfd=n - q)
     return 1 - cdf

tests/test_k_class.py

@@ -28,7 +28,7 @@ def test__fuller_alpha(kappa, expected):
 
 
 @pytest.mark.parametrize("n, p, q, u", [(100, 2, 2, 1), (100, 2, 5, 2)])
-def test__eta_liml_same_with_shortcut(n, p, q, u):
+def test__kappa_liml_same_with_shortcut(n, p, q, u):
     Z, X, y = simulate_gaussian_iv(n, p, q, u)
 
     X = X - X.mean(axis=0)
@@ -39,24 +39,24 @@ def test__eta_liml_same_with_shortcut(n, p, q, u):
 
     k_class = KClass(kappa="liml")
 
-    lbda = k_class.fit(X, y, Z).eta_liml_
-    assert np.allclose(lbda, k_class._eta_liml(X, y, X_proj=X_proj, y_proj=y_proj))
-    assert np.allclose(lbda, k_class._eta_liml(X, y, Z=Z, X_proj=X_proj))
-    assert np.allclose(lbda, k_class._eta_liml(X, y, Z=Z, y_proj=y_proj))
-    assert np.allclose(lbda, k_class._eta_liml(X, y, Z=Z))
+    lbda = k_class.fit(X, y, Z).kappa_liml_
+    assert np.allclose(lbda, k_class._kappa_liml(X, y, X_proj=X_proj, y_proj=y_proj))
+    assert np.allclose(lbda, k_class._kappa_liml(X, y, Z=Z, X_proj=X_proj))
+    assert np.allclose(lbda, k_class._kappa_liml(X, y, Z=Z, y_proj=y_proj))
+    assert np.allclose(lbda, k_class._kappa_liml(X, y, Z=Z))
 
 
 @pytest.mark.parametrize("n, p, q, u", [(100, 2, 2, 1), (100, 2, 5, 2)])
-def test__eta_liml_positive(n, p, q, u):
+def test__kappa_liml_positive(n, p, q, u):
     Z, X, y = simulate_gaussian_iv(n, p, q, u)
 
     k_class = KClass(kappa="liml")
     k_class.fit(X, y.flatten(), Z)
 
     if Z.shape[1] > X.shape[1]:
-        assert k_class.eta_liml_ > 0
+        assert k_class.kappa_liml_ > 1
     else:
-        assert np.allclose(k_class.eta_liml_, 0)
+        assert np.allclose(k_class.kappa_liml_, 1)
 
 
 @pytest.mark.parametrize("n, p, q, u", [(100, 2, 2, 1), (100, 2, 5, 2)])
@@ -110,7 +110,7 @@ def test_liml_equal_to_tsls_in_just_identified_setting(n, p, q, u):
 
 
 @pytest.mark.parametrize("n, p, q, u", [(100, 2, 2, 1), (100, 4, 4, 2)])
-def test_anderson_rubin_at_liml_is_equal_to_eta_liml(n, p, q, u):
+def test_anderson_rubin_at_liml_is_equal_to_kappa_liml(n, p, q, u):
     Z, X, y = simulate_gaussian_iv(n, p, q, u)
     y = y.flatten()
 
@@ -119,7 +119,7 @@ def test_anderson_rubin_at_liml_is_equal_to_eta_liml(n, p, q, u):
 
     assert np.allclose(
         anderson_rubin_test(Z, X, y, liml.coef_)[0],
-        (n - q) / q * liml.eta_liml_ / (1 - liml.eta_liml_),
+        (n - q) / q * (liml.kappa_liml_ - 1),
         atol=1e-8,
     )
 
@@ -158,7 +158,6 @@ def test_fuller_bias_and_mse(n, beta, Pi, gamma, delta):
 
     kappas = ["liml", "fuller(1)", "fuller(4)", 0]  # 0 is for OLS
     results = {kappa: np.zeros(shape=(n_iterations, p)) for kappa in kappas}
-    limls = np.zeros(shape=n_iterations)
 
     for seed in range(n_iterations):
         rng = np.random.RandomState(seed)
@@ -169,7 +168,6 @@ def test_fuller_bias_and_mse(n, beta, Pi, gamma, delta):
 
         for kappa in kappas:
             results[kappa][seed, :] = KClass(kappa=kappa).fit(X, y, Z).coef_
-        limls[seed] = KClass(kappa="liml").fit(X, y, Z).eta_liml_
 
     mses = {k: np.mean((v - beta.flatten()) ** 2) for k, v in results.items()}
     # biases = {k: np.mean(v - beta.flatten(), axis=0) for k, v in results.items()}