From 1667dc3fb288df3585db87cf45e01b4f1da86969 Mon Sep 17 00:00:00 2001
From: saidamir <amamatov@yahoo.com>
Date: Mon, 17 Jun 2024 08:57:11 -0400
Subject: [PATCH 1/9] Added solver to feols and created new test file for it

---
 pyfixest/estimation/feols_.py | 11 +++++--
 tests/test_solver.py          | 61 +++++++++++++++++++++++++++++++++++
 2 files changed, 69 insertions(+), 3 deletions(-)
 create mode 100644 tests/test_solver.py

diff --git a/pyfixest/estimation/feols_.py b/pyfixest/estimation/feols_.py
index df7e9046..3050d811 100644
--- a/pyfixest/estimation/feols_.py
+++ b/pyfixest/estimation/feols_.py
@@ -173,6 +173,7 @@ def __init__(
         coefnames: list[str],
         weights_name: Optional[str],
         weights_type: Optional[str],
+        solver: str = "np.linalg.solve",
     ) -> None:
         self._method = "feols"
         self._is_iv = False
@@ -193,7 +194,7 @@ def __init__(
             self._X = X
 
         self.get_nobs()
-
+        self._solver = solver
         _feols_input_checks(Y, X, weights)
 
         if self._X.shape[1] == 0:
@@ -302,12 +303,16 @@ def get_fit(self) -> None:
         _X = self._X
         _Y = self._Y
         _Z = self._Z
-
+        _solver = self._solver
         self._tZX = _Z.T @ _X
         self._tZy = _Z.T @ _Y
 
         # self._tZXinv = np.linalg.inv(self._tZX)
-        self._beta_hat = np.linalg.solve(self._tZX, self._tZy).flatten()
+        if _solver == "np.linalg.lstsq":
+            self._beta_hat, _, _, _ = np.linalg.lstsq(self._tZX, self._tZy, rcond=None)
+        else:  # Default to np.linalg.solve
+            self._beta_hat = np.linalg.solve(self._tZX, self._tZy).flatten()
+
         # self._beta_hat, _, _, _ = lstsq(self._tZX, self._tZy, lapack_driver='gelsy')
 
         # self._beta_hat = (self._tZXinv @ self._tZy).flatten()
diff --git a/tests/test_solver.py b/tests/test_solver.py
new file mode 100644
index 00000000..619e9db6
--- /dev/null
+++ b/tests/test_solver.py
@@ -0,0 +1,61 @@
+import unittest
+
+import numpy as np
+
+from pyfixest import feols
+
+
+class TestSolvers(unittest.TestCase):
+    """A test case class for testing the solvers in the feols class."""
+
+    def setUp(self):
+        """
+        Set up the test environment.
+
+        This method is called before each test case is executed.
+        It initializes the test data
+        by generating random arrays for X, Y, and Z.
+        """
+        np.random.seed(42)  # Ensure reproducibility
+        self._X = np.random.rand(100, 5)
+        self._Y = np.random.rand(100, 1)
+        self._Z = np.random.rand(
+            100, 5
+        )  # Assuming Z is used similarly to X in your context
+
+    def test_solver_equivalence(self):
+        """
+        Test the equivalence of different solvers for the feols class.
+
+        This method initializes an object with test data and compares the results
+        obtained from two different solvers: np.linalg.lstsq
+        and np.linalg.solve. It asserts that
+        the results are identical or very close within a tolerance.
+
+        Returns
+        -------
+            None
+        """
+        obj = feols()  # Replace with your actual class initialization
+        obj._X = self._X
+        obj._Y = self._Y
+        obj._Z = self._Z
+
+        # Use np.linalg.lstsq solver
+        obj._solver = "np.linalg.lstsq"
+        obj.get_fit()
+        beta_hat_lstsq = obj._beta_hat.copy()
+
+        # Use np.linalg.solve solver
+        obj._solver = "np.linalg.solve"
+        obj.get_fit()
+        beta_hat_solve = obj._beta_hat.copy()
+
+        # Assert that the results are identical or very close within a tolerance
+        np.testing.assert_array_almost_equal(
+            beta_hat_lstsq, beta_hat_solve, decimal=6, err_msg="Solvers' results differ"
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()

From 88808b9828a710d3ef571d16f0421e9eba17beb1 Mon Sep 17 00:00:00 2001
From: saidamir <amamatov@yahoo.com>
Date: Wed, 19 Jun 2024 22:13:20 -0400
Subject: [PATCH 2/9] Updated test file and feos class with new solver

---
 pyfixest/estimation/feols_.py |   5 ++
 tests/test_solver.py          | 111 ++++++++++++++++------------------
 2 files changed, 56 insertions(+), 60 deletions(-)

diff --git a/pyfixest/estimation/feols_.py b/pyfixest/estimation/feols_.py
index 3050d811..b9957f2e 100644
--- a/pyfixest/estimation/feols_.py
+++ b/pyfixest/estimation/feols_.py
@@ -55,6 +55,9 @@ class Feols:
     weights_type : Optional[str]
         Type of the weights variable. Either "aweights" for analytic weights or
         "fweights" for frequency weights.
+     solver : str, optional.
+        The solver to use for the regression. Can be either "np.linalg.solve" or
+        "np.linalg.lstsq". Defaults to "np.linalg.solve".
 
     Attributes
     ----------
@@ -79,6 +82,8 @@ class Feols:
         Indices of collinear variables.
     _Z : np.ndarray
         Alias for the _X array, used for calculations.
+    _solver: str
+        The solver used for the regression.
     _weights : np.ndarray
         Array of weights for each observation.
     _N : int
diff --git a/tests/test_solver.py b/tests/test_solver.py
index 619e9db6..a2ba6fdb 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -1,61 +1,52 @@
-import unittest
-
 import numpy as np
-
-from pyfixest import feols
-
-
-class TestSolvers(unittest.TestCase):
-    """A test case class for testing the solvers in the feols class."""
-
-    def setUp(self):
-        """
-        Set up the test environment.
-
-        This method is called before each test case is executed.
-        It initializes the test data
-        by generating random arrays for X, Y, and Z.
-        """
-        np.random.seed(42)  # Ensure reproducibility
-        self._X = np.random.rand(100, 5)
-        self._Y = np.random.rand(100, 1)
-        self._Z = np.random.rand(
-            100, 5
-        )  # Assuming Z is used similarly to X in your context
-
-    def test_solver_equivalence(self):
-        """
-        Test the equivalence of different solvers for the feols class.
-
-        This method initializes an object with test data and compares the results
-        obtained from two different solvers: np.linalg.lstsq
-        and np.linalg.solve. It asserts that
-        the results are identical or very close within a tolerance.
-
-        Returns
-        -------
-            None
-        """
-        obj = feols()  # Replace with your actual class initialization
-        obj._X = self._X
-        obj._Y = self._Y
-        obj._Z = self._Z
-
-        # Use np.linalg.lstsq solver
-        obj._solver = "np.linalg.lstsq"
-        obj.get_fit()
-        beta_hat_lstsq = obj._beta_hat.copy()
-
-        # Use np.linalg.solve solver
-        obj._solver = "np.linalg.solve"
-        obj.get_fit()
-        beta_hat_solve = obj._beta_hat.copy()
-
-        # Assert that the results are identical or very close within a tolerance
-        np.testing.assert_array_almost_equal(
-            beta_hat_lstsq, beta_hat_solve, decimal=6, err_msg="Solvers' results differ"
-        )
-
-
-if __name__ == "__main__":
-    unittest.main()
+import pytest
+
+import pyfixest as pf
+from pyfixest.estimation.feols_ import Feols
+
+
+@pytest.fixture()
+def data():
+    data = pf.get_data()
+    Y = data["Y"]
+    X = data[["X1", "X2"]]
+    Z = data[["Z1"]]
+    weights = data["weights"]  # Define the variable "weights"
+    return X, Y, Z, weights
+
+
+def test_solver_equivalence(data):
+    """
+    Test the equivalence of different solvers for the feols class.
+    This function initializes an object with test data and compares the results
+    obtained from two different solvers: np.linalg.lstsq
+    and np.linalg.solve. It asserts that
+    the results are identical or very close within a tolerance.
+    """
+    X, Y, Z, weights = data
+
+    obj = Feols(
+        X=X,
+        Y=Y,
+        weights=weights,
+        collin_tol=1e-08,
+        coefnames=["X1", "X2"],
+        weights_name="weights",
+        weights_type=None,
+        solver="np.linalg.solve",
+    )
+
+    # Use np.linalg.lstsq solver
+    obj._solver = "np.linalg.lstsq"
+    obj.get_fit()
+    beta_hat_lstsq = obj._beta_hat.copy()
+
+    # Use np.linalg.solve solver
+    obj._solver = "np.linalg.solve"
+    obj.get_fit()
+    beta_hat_solve = obj._beta_hat.copy()
+
+    # Assert that the results are identical or very close within a tolerance
+    np.testing.assert_array_almost_equal(
+        beta_hat_lstsq, beta_hat_solve, decimal=6, err_msg="Solvers' results differ"
+    )

From 3ca68c48afeaf05e7a0a35c7e6513e40e84174c9 Mon Sep 17 00:00:00 2001
From: saidamir <amamatov@yahoo.com>
Date: Sun, 30 Jun 2024 11:15:07 -0400
Subject: [PATCH 3/9] Adding solver to feiv and fepois, and upding test_solver
 file

---
 pyfixest/estimation/feiv_.py   | 31 ++++++++++++-
 pyfixest/estimation/fepois_.py | 31 ++++++++++++-
 tests/test_solver.py           | 82 ++++++++++++++++++++++++++++++++--
 3 files changed, 139 insertions(+), 5 deletions(-)

diff --git a/pyfixest/estimation/feiv_.py b/pyfixest/estimation/feiv_.py
index 1603bca6..39130e02 100644
--- a/pyfixest/estimation/feiv_.py
+++ b/pyfixest/estimation/feiv_.py
@@ -94,6 +94,7 @@ def __init__(
         collin_tol: float,
         weights_name: Optional[str],
         weights_type: Optional[str],
+        solver: str = "np.linalg.solve",
     ) -> None:
         super().__init__(
             Y=Y,
@@ -103,6 +104,7 @@ def __init__(
             collin_tol=collin_tol,
             weights_name=weights_name,
             weights_type=weights_type,
+            solver=solver,
         )
 
         if self._has_weights:
@@ -127,6 +129,32 @@ def __init__(
         self._support_iid_inference = True
         self._supports_cluster_causal_variance = False
 
+    def solve_ols(self, tZX, tZY, solver):
+        """
+        Solve the ordinary least squares problem using the specified solver.
+
+        Parameters
+        ----------
+        tZX (array-like): The design matrix.
+        tZY (array-like): The response variable.
+        solver (str): The solver to use. Supported solvers are "np.linalg.lstsq"
+        and "np.linalg.solve".
+
+        Returns
+        -------
+        array-like: The solution to the ordinary least squares problem.
+
+        Raises
+        ------
+        ValueError: If the specified solver is not supported.
+        """
+        if solver == "np.linalg.lstsq":
+            return np.linalg.lstsq(tZX, tZY, rcond=None)[0]
+        elif solver == "np.linalg.solve":
+            return np.linalg.solve(tZX, tZY).flatten()
+        else:
+            raise ValueError(f"Solver {solver} not supported.")
+
     def get_fit(self) -> None:
         """Fit a IV model using a 2SLS estimator."""
         _X = self._X
@@ -141,8 +169,9 @@ def get_fit(self) -> None:
         H = self._tXZ @ self._tZZinv
         A = H @ self._tZX
         B = H @ self._tZy
+        solver = self._solver
 
-        self._beta_hat = np.linalg.solve(A, B).flatten()
+        self._beta_hat = self.solve_ols(A, B, solver)
 
         self._Y_hat_link = self._X @ self._beta_hat
         self._u_hat = self._Y.flatten() - self._Y_hat_link.flatten()
diff --git a/pyfixest/estimation/fepois_.py b/pyfixest/estimation/fepois_.py
index 4b53fe14..60f28166 100644
--- a/pyfixest/estimation/fepois_.py
+++ b/pyfixest/estimation/fepois_.py
@@ -68,6 +68,7 @@ def __init__(
         maxiter: int = 25,
         tol: float = 1e-08,
         fixef_tol: float = 1e-08,
+        solver: str = "np.linalg.solve",
         weights_name: Optional[str] = None,
         weights_type: Optional[str] = None,
     ):
@@ -79,6 +80,7 @@ def __init__(
             collin_tol=collin_tol,
             weights_name=weights_name,
             weights_type=weights_type,
+            solver=solver,
         )
 
         # input checks
@@ -113,6 +115,32 @@ def __init__(
         self.deviance = None
         self._Xbeta = np.array([])
 
+    def solve_ols(self, tZX, tZY, solver):
+        """
+        Solve the ordinary least squares problem using the specified solver.
+
+        Parameters
+        ----------
+        tZX (array-like): The design matrix.
+        tZY (array-like): The response variable.
+        solver (str): The solver to use. Supported solvers are "np.linalg.lstsq"
+        and "np.linalg.solve".
+
+        Returns
+        -------
+        array-like: The solution to the ordinary least squares problem.
+
+        Raises
+        ------
+        ValueError: If the specified solver is not supported.
+        """
+        if solver == "np.linalg.lstsq":
+            return np.linalg.lstsq(tZX, tZY, rcond=None)[0]
+        elif solver == "np.linalg.solve":
+            return np.linalg.solve(tZX, tZY).flatten()
+        else:
+            raise ValueError(f"Solver {solver} not supported.")
+
     def get_fit(self) -> None:
         """
         Fit a Poisson Regression Model via Iterated Weighted Least Squares (IWLS).
@@ -151,6 +179,7 @@ def get_fit(self) -> None:
         _iwls_maxiter = 25
         _tol = self.tol
         _fixef_tol = self.fixef_tol
+        solver = self._solver
 
         def compute_deviance(_Y: np.ndarray, mu: np.ndarray):
             with warnings.catch_warnings():
@@ -215,7 +244,7 @@ def compute_deviance(_Y: np.ndarray, mu: np.ndarray):
             XWX = WX.transpose() @ WX
             XWZ = WX.transpose() @ WZ
 
-            delta_new = np.linalg.solve(XWX, XWZ)  # eq (10), delta_new -> reg_z
+            delta_new = self.solve_ols(XWX, XWZ, solver)  # eq (10), delta_new -> reg_z
             resid = Z_resid - X_resid @ delta_new
 
             mu_old = mu.copy()
diff --git a/tests/test_solver.py b/tests/test_solver.py
index a2ba6fdb..b06f3717 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -2,19 +2,95 @@
 import pytest
 
 import pyfixest as pf
+from pyfixest.estimation.feiv_ import Feiv
 from pyfixest.estimation.feols_ import Feols
+from pyfixest.estimation.fepois_ import Fepois
 
 
 @pytest.fixture()
 def data():
     data = pf.get_data()
-    Y = data["Y"]
+    Y = data["Y"].to_numpy()
     X = data[["X1", "X2"]]
-    Z = data[["Z1"]]
-    weights = data["weights"]  # Define the variable "weights"
+    Z = data[["Z1"]].to_numpy()
+    weights = data["weights"].to_numpy()  # Define the variable "weights"
     return X, Y, Z, weights
 
 
+def test_solver_fepois(data):
+    """
+    Test the equivalence of different solvers for the fepois class.
+    This function initializes an object with test data and compares the results
+    obtained from two different solvers: np.linalg.lstsq
+    and np.linalg.solve. It asserts that
+    the results are identical or very close within a tolerance.
+    """
+    X, Y, Z, weights = data
+    obj = Fepois(
+        X=X,
+        Y=Y,
+        weights=weights,
+        coefnames=["X1", "X2"],
+        collin_tol=1e-08,
+        weights_name="weights",
+        weights_type=None,
+        solver="np.linalg.solve",
+        fe=None,
+        drop_singletons=False,
+    )
+    # Use np.linalg.lstsq solver
+    obj._solver = "np.linalg.lstsq"
+    obj.get_fit()
+    beta_hat_lstsq = obj._beta_hat.copy()
+    # Use np.linalg.solve solver
+    obj._solver = "np.linalg.solve"
+    obj.get_fit()
+    beta_hat_solve = obj._beta_hat.copy()
+    # Assert that the results are identical or very close within a tolerance
+    np.testing.assert_array_almost_equal(
+        beta_hat_lstsq, beta_hat_solve, decimal=6, err_msg="Solvers' results differ"
+    )
+
+
+def test_solver_feiv(data):
+    """
+    Test the equivalence of different solvers for the feiv class.
+    This function initializes an object with test data and compares the results
+    obtained from two different solvers: np.linalg.lstsq
+    and np.linalg.solve. It asserts that
+    the results are identical or very close within a tolerance.
+    """
+    X, Y, Z, weights = data
+
+    obj = Feiv(
+        Y=Y,
+        X=X,
+        Z=Z,
+        weights=weights,
+        coefnames_x=["X1", "X2"],
+        collin_tol=1e-08,
+        weights_name="weights",
+        weights_type=None,
+        solver="np.linalg.solve",
+        coefnames_z=["Z1"],
+    )
+
+    # Use np.linalg.lstsq solver
+    obj._solver = "np.linalg.lstsq"
+    obj.get_fit()
+    beta_hat_lstsq = obj._beta_hat.copy()
+
+    # Use np.linalg.solve solver
+    obj._solver = "np.linalg.solve"
+    obj.get_fit()
+    beta_hat_solve = obj._beta_hat.copy()
+
+    # Assert that the results are identical or very close within a tolerance
+    np.testing.assert_array_almost_equal(
+        beta_hat_lstsq, beta_hat_solve, decimal=6, err_msg="Solvers' results differ"
+    )
+
+
 def test_solver_equivalence(data):
     """
     Test the equivalence of different solvers for the feols class.

From 1bea8f088fd74e2532cac510d783553a864d3cec Mon Sep 17 00:00:00 2001
From: Alexander Fischer <alexander-fischer1801@t-online.de>
Date: Sun, 30 Jun 2024 21:55:40 +0200
Subject: [PATCH 4/9] update tests

---
 tests/test_solver.py | 27 +++++++++++++++++++--------
 1 file changed, 19 insertions(+), 8 deletions(-)

diff --git a/tests/test_solver.py b/tests/test_solver.py
index b06f3717..826f436a 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -8,16 +8,26 @@
 
 
 @pytest.fixture()
-def data():
+def data_feols():
     data = pf.get_data()
     Y = data["Y"].to_numpy()
-    X = data[["X1", "X2"]]
+    X = data[["X1", "X2"]].to_numpy()
     Z = data[["Z1"]].to_numpy()
     weights = data["weights"].to_numpy()  # Define the variable "weights"
     return X, Y, Z, weights
 
 
-def test_solver_fepois(data):
+@pytest.fixture()
+def data_fepois():
+    data = pf.get_data(model="Fepois")
+    Y = data["Y"].to_numpy()
+    X = data[["X1", "X2"]].to_numpy()
+    Z = data[["Z1"]].to_numpy()
+    weights = data["weights"].to_numpy()  # Define the variable "weights"
+    return X, Y, Z, weights
+
+
+def test_solver_fepois(data_fepois):
     """
     Test the equivalence of different solvers for the fepois class.
     This function initializes an object with test data and compares the results
@@ -25,7 +35,7 @@ def test_solver_fepois(data):
     and np.linalg.solve. It asserts that
     the results are identical or very close within a tolerance.
     """
-    X, Y, Z, weights = data
+    X, Y, Z, weights = data_fepois
     obj = Fepois(
         X=X,
         Y=Y,
@@ -52,7 +62,7 @@ def test_solver_fepois(data):
     )
 
 
-def test_solver_feiv(data):
+def test_solver_feiv(data_feols):
     """
     Test the equivalence of different solvers for the feiv class.
     This function initializes an object with test data and compares the results
@@ -60,12 +70,13 @@ def test_solver_feiv(data):
     and np.linalg.solve. It asserts that
     the results are identical or very close within a tolerance.
     """
-    X, Y, Z, weights = data
+    X, Y, Z, weights = data_feols
 
     obj = Feiv(
         Y=Y,
         X=X,
         Z=Z,
+        endogvar=X[:, 0],
         weights=weights,
         coefnames_x=["X1", "X2"],
         collin_tol=1e-08,
@@ -91,7 +102,7 @@ def test_solver_feiv(data):
     )
 
 
-def test_solver_equivalence(data):
+def test_solver_feols(data_feols):
     """
     Test the equivalence of different solvers for the feols class.
     This function initializes an object with test data and compares the results
@@ -99,7 +110,7 @@ def test_solver_equivalence(data):
     and np.linalg.solve. It asserts that
     the results are identical or very close within a tolerance.
     """
-    X, Y, Z, weights = data
+    X, Y, Z, weights = data_feols
 
     obj = Feols(
         X=X,

From 945f611a151866f5ee3c187093e41cfb4a385f1f Mon Sep 17 00:00:00 2001
From: Alexander Fischer <alexander-fischer1801@t-online.de>
Date: Sun, 30 Jun 2024 21:56:31 +0200
Subject: [PATCH 5/9] move solve_ols to Feols class, delete in other methods;
 add docstrings

---
 .coverage                      | Bin 53248 -> 53248 bytes
 pyfixest/estimation/feiv_.py   |  34 ++++----------------------------
 pyfixest/estimation/feols_.py  |  35 ++++++++++++++++++++++++++++-----
 pyfixest/estimation/fepois_.py |  33 +++++--------------------------
 4 files changed, 39 insertions(+), 63 deletions(-)

diff --git a/.coverage b/.coverage
index cb311f7f41fb7da06eaa4ff2efd26a05b933c0c6..dc0b802595222df324e5f93ee87a4cf6b6fc263b 100644
GIT binary patch
literal 53248
zcmeI5eQX=$8Nly+7oX4e`E}wXPMWs44Q-*&1`1=eg|$i3KnrY8RzPgg+>3o_j`<__
z&S?UI#%-a9iD~qoZ4zQknzlbCP2C?LjZLG{KKR(ihbC>4rj2Y8(6O=4f?zPR=XvkW
zKDP<t+({MPbBdq$?%wx#?s<Oi*S)sS`SkAVEXUNRD%FzV=$8s<K@^1%T^9sFg_i~|
ze^X#b_75Nx1N$rOs>1Ls*J$i!AsKx@V4u^jXB(0aX&+6zkSwY9CmvP|H4hoI5dk7V
z1c<=@H-Y-@q>}CL7Z2X$7?VZQsTz5+>Pl_w#_`cT<NBV_?bnU#?m7MZxDMNf4f?2F
zt?bul%&I<R6;0hL7p%PDSmk}XGi~~%YqnXyi4M-eREGtJxjydYtO95`<~}enQ?*J)
zbzc9xIlrD4@UG^VbB=ogjF{HGGM3QK_te+vRddR$n&rG{yI!1c71qVq9Xb$KvYCu{
z#N;s<)ja$qTB3my%}{ntfuVe*T0rk7t429LZQAQ~qY+`g0)EapZnV`(=_0417p*c6
zvus<ARVnM{oSCmVW?|U?*CRB419-t`z}8+2Xs@g3T`NGRsl1G)SKL#c0_VyVA8S()
zFAd7v=W?wy3CM4&m1eMVK*PPTQ&~xKC0(cD1^d41W|_Kd#@kYzK60#4SF2XGZm*e*
z%3rTHHP`bn%9Zsx1_3SLBZ+8sb4C=PRho15?M1lX${97M!jE%M`MHfCl7~BEN_KEi
zJaRj4D2OLFX*qUdM{Q{|jngYM9n4?3-C!?iZ6w%>ibcIL#k&k@>u%*hYYz6Ai+mVH
z6+4<<@d|>cXSKDVv*3e1S*%R9wl&^jY!f;Nj0JwNsaiRYUB$*2@|9V$YV0$Iy<B%#
zINJ>gn$AWwnjIfpsp&MDKyG7WH(u|LD%pVn@u13^kyl$@qO>$2zb2O%M*XakY1p-@
zKajlFAn9u~M@hercsUTx1|^yu9at$GuUu}Uw-c)?vXbrT5$kS`z;Ov*9W6oNnPwL_
zWg>Y0?#~#`w0^qpcp4lZ>WF5q=vm2ebRoA9W%STzBi!X9{hWRY27P=<EiU7krpwow
z^BZQ3B8(<8hE;_w3uCZ51L#<=U_fDGvQl$+`^7`vJaO)6`s*5K9{cQ})H+@qy7vtg
z+qPg@#c9jl;-A%rrVSfMA$v%ltW=7oQFiAB(AZ_eH7#4`KAb))ZEM=F=h4Q#*&*(e
zKe^yAzXv*<-@jm|R5Hp1Z@$85*tH5{KXy*Mjys5eqCNzj(uOX2icR0Bmy08oZ-00O
zn4r~?QE<CER!YF){%P|xbM9=k!*EQ;Dw#SK;m+`ewcoJyeAR?Zm`m3UN$=GRu(om=
zQQ|5&hPQI{K~c$Gc%fLYd7aT6XL2@7>LtT>+dGW3G?G&?o9*R#7;xckC)iGqT>u9j
zx>vx#&8@*nmx5o)yVr@Y!v(&`-WK4GHX=X-hyW2F0z`la5CI}U1c(3;AOb|-+#?`~
z5mCYQe}w%(V1ENa8xbG^M1Tko0U|&IhyW2F0z`la5CJ0a!6cB3NNbq?6_1OVC=F!2
zF9BRWeA)0Nafm9yP6+G-d;NnMAbCv$hyW2F0z`la5CI}U1c(3;AOb{y2>1k&(i+i!
z7a$Us22}4&0Q~%aP<d8hKVkQ<?Wwm@kERN#0qs}XgW66lntUpGPjXx0Wa3z&km!s*
z6F(9kQ<u~y)d$qwDvSL*_NCZOu?(owMg)ie5g-CYfC!v50%L=+P*03)nYhV@XD26&
zqB&<w?5N~xC9~|<6EpKu)|_cO&7GQK743<fIhQwQ_~X#<2<6CHFwzmi2v~$)OhEOC
ze0~-zY+fS?^C1irtipt|f64^aEvqHr&XB6b88864qXXc$5>nZm1vEL^sTJli_EjNc
zN7tM2V*vf(0{|nJ*r2<qPZIWr^xqBr<5@|t!)Th5wS92K?Y$7M8bYyZ&Qz+-#H?l4
zjG}c1R@wHR_H5WSco?|q+8FH)S#_QtxtvwDoLmmv7|(zS79x~pwM~?2MF*Zuo&htv
zyTHs)n`SBnvzUXfP_^cAZdXIkuIU7m{q32&3VxWYRy1xfs?ehmb{cGjsjxtlQ>HbG
z{_J94Bip7wd;s|*@Uqa<jue>E+A}o;P06y+%BTib!gO4B*eaN_Ilc!bu1re8fe?PU
z7JTiI1UL|8$aB>Lt?%kMSjdEMAW$!}@T_+ukTDeug&D#FhN|$>BDMoG#$xRonKCOy
zUMp8C?T6HK3cR^K3Wm;Wk4?dJ%sh<G@M}2mgFx_m9lX_bfWbCfUBQGOg(#zyaT%<%
zIbZ~8epCW8VH&SltyW-i4w%>)0TW^fKsf6J^f!y)=G%ZOe*QnGJSVW9s!y?d)i0_y
zs_ED-Vh_e{W@Cxh6K3k2Sa<5Nl$lzc_^$Tr<Vd0`{xj{Uwkv*Ld`Dc+6y>?(y~(G~
z+6@m~m<SL7B0vO)01*&`9dbPp=0?K1k8N51501+9jxa`=H+$atf9)1&K8yjpo$=QH
zYp#&)46Ev||5tC9>q=PV*7g5D*x1qaX8dqhxvGr-zzARe_ivW=hxH$H;i`R?OLizt
zZ~dPglj~|2#bwt2y=~eEuK#<&R^93A|L&{gItvp@v)Y2||ICeYeW+bC%dh{tu9fTk
zZJIoF{offDb|6ZB{hz)@u4mizr*-|$#^t)!rm5!oKNYt2-C--Z{?|Ss9SGxxYk^DX
zKFP4%-WxZY?>;4V%Jobb2Lkodvi^^U9l`^K8tZ@c42(3^|FJU&sp%BZRc?^$=e5aZ
z+4X<)3|d`t{oip014dv4Er)Hqr?3B|ZE{@<0|;lGfD4h$(tMcr1Ouw{`+s!XNCb!g
z5g-CYfCvx)B0vO)01+SpMBtnzAi=i=H2nNuWbX*@M;j3!0z`la5CI}U1c(3;AOb{y
z2oM1xa1IiXBu&N7|9cZ(5!k=jpV@EO^XwS=2K!8ETdG5QO?yfEp7yY2X>#(_<lV`w
z$-!iA;yEZx8xbG^M1Tko0U|&IhyW2F0{>G2m-Wh`TI{_;Sa=ESo;=WE{P>a*IF=W;
zxktwjh~qsSqFRsl%zkYt5-76FZmMz4-5bGU{Kw_6iApzIJ>C6bc})g(uFO1^aWmeO
z+|-lu)-K4tyX&tnvzGJzfw&s>;8Aq^U+5!(;5i);bU4-pQApp`iSF*{96rv^!D-=y
zdx^VnUl96GoD}6Wrk;O`{9;(F83P=Mem<o9ay$j$1{Lc|gRhsQF9?M2pC#lP$usog
zD)#$#Z}mu?7<82vw)#d^N_C8c|F3uCEvAKtF32NGjf;Sb3wRK_GzFKRPCd4XXN@eq
zh?f>lucYG`hnKOpUjwC&X|pN>vDmPg9uf4|1xd)dIXU|llz-%zC%HAy16LNFEZ-9C
zTO9oUL%0438;=l+CZKF8aj-WI>2Un<w>|_;_X*-&%oP^dCKb+C)i<PL|Mcw~Z-(%W
z3*$IX7T|Kpze7!O-N?bd7_d8haqNdLj*L9T%>~8D15X57DprLK`Cg2{l7i?`#EM?P
zj1d&PwD7#^l+Mp^aY}z^W(WfESZ3o@#Ls(e9QV`c5Z^F`78x56{@vnmdg0v$cCV<L
zpr5FOuSIy_=7*)@!NN`jTUA<l^_8XNdk2<k-QH|MZ)7YkiSZ~@PjB>7NkNFN?STEQ
z9b1-84O3bYHVJPkstozv^3k0~52qy9T`Ns_TO=BR<5VP)5M@<eCC3ixg8TV@FPj$F
z+w3HJjlIlXWWQlAuxHqh*bmsF>^tmR>=E`=c9eaA9cFj2Syq9<v=IR!Km>>Y5g-CY
zfCvx)B0vO)01+Sp%M!pPyx7~zlb#-)ba(S4li^8M7f(7nd6G`^gfX6^QasT#o+Ohz
zNhEj@=kE@Ps>+XIF`g(2Pohztbae1UmU$vcJc&ePcpHF!{txei&_)D^01+SpM1Tko
z0U|&IhyW2F0z`laoWlh0`F~pfpTiYQp%DQhKm>>Y5g-CYfCvx)B0vO)01*HJ`27DM
zGx7WXZ?J!`zpz)>AKCBV8Gv81XW7&2$LtCAI6Mb%jD4GZlYN~%#J&vA0zAO(V~5z?
zY@RvrJiu+BP8$&*0z`la5CI}U1c(3;AOb{y2oM1xAPvGV?MrJ>tU<9F#Q=&`DEd+K
zp~#}>MbU$z8$||17m7|4X%q}a3WbIui6Vg_jzUEdL!qFEqUb;&qmWQUP~h+X{}+iY
BIHCXm

literal 53248
zcmeI5eUKbQ6~KFT_G|X5_r7zv+}%hCmyi3%<p_kJ$t4ggD3CxR5(8}JcIR@F?CxA<
zXD@ern7t51Q3RtFWf34&iIim-Q>BR(id-2i0+tdAtU?8Hpd=`$C;<b+tnbb2-tHxL
zyBq(3#5Z;OdZwpe_j|v7-Lt*BGq>iND>YM<`*kCsm~sybGeKazvdkEh;1z^de)GVN
zGk*Z7P_#eYuEaW*{yxaBV*&SFjNcSo#XADKf*1Io2_&RD{SSB)DGDcOBLYN#2oQn)
zZvvy&1UwCmjl%ZZOr<ZbnuZcp4O<2;yLQF0wJYSc%a&icLbmhd**+Pzjt+U5Z0K9%
zLDi7^wYVy4$(R;ZOf9)lHV4#v>2yks;e`(J;8KSLJ1u?Oi)b;>GS!Vx#h{@j6k|l*
ztd6u>1#+#K>ab~NKqacSF^MJQ**Wz&vZ3~?hMJ73Dcg$KT5OJQ&h6WLo`#wl;SSaE
zNim}E=bz*aWE5Pvq8}=X>P8H$?=zHSbU;nD%gTg@Q620YHf?VWJ<(<Aa$HMVZYEQj
zsp&~s9af`hQ;kjOU^_y~Zvd~@8>njoI<&#o%(YcSr=UEAX0CWnxgU&6>UpdMMQj?B
zEt`|+L?0mEluis{=YWRYu&Ga{Ii0Sl=bWv$Z>H$Wg1?ixGmjj5)YeKTwJm9NqVwD3
zf@ZtrMpAE=(Furv^Zo9IOKJoLQK=56w#4E6Rzyjgx^*0Z&X0712<)ixdKy|=g*&!b
zfr5A<eVUn?*pVg$%|!Mzp@ZkAj~g5%lY<0DQ9Lf|{nn5{ZQ1S|i00xk)0Kyz7&y>m
zJy#Gc{XnA)jRhOzzPR2uIo7O*NvSYEU@piP>rW@6I8;*TLsTD94P~R+nLBGw3m<eK
zLFjZ#?uHeu(}m7N2t>Lkc73CbZcjr~ldxT~!YJ2Uxx_Omgz`N(#Wb2fl~9#b+Q?5N
zr6$R|Mq!rBkCEINaHqrLZdledojbX5k?!1%chuwbG}P4zqxOivc?n;xNlsu+VF;XZ
zA*kH^LB$-9Pq!Vj!SFVhyW!Hh=?q5`BHbw6x370tCg0I$(NDmnk1uhO$(U0x`BUmh
z$B+_-*<?`B3>dO72isQw84Jb~D6I79Y14{dJhZ|SuRR5O?FtIVdG=6hj#V6n_i8<r
z(%`a+m#x&|e6HL!prmjXO0~&-x*k`Rq<w7wjVn^nrk0W|8&1C}EiKfM8o@d?4z*b}
z<u5Kc&98+)mmgoSlSnAZSnhg-mtnhBnEP>XVmtOE0*Z1Q49XM?(NpToSDHJEGgiL-
z&Ou;;wI-C9J>0QVifY_Cpq^HpeKopVF;!Dbs4^C@jIjc1tCEtVh6*R)S~{mqEENXW
zTaj**xJvfoTe;Cz!P78ro-mru4MuyOiKO77o>1~;=MEi{g5;Ewg?L#76E3{lDUK6l
zSA+u(?G|uyvwN`5rjRd{d#~f0gA06tzt7;0HX=X-hyW2F0z`la5CI}U1c(3;AOb|-
z+#?_g4#9)#e+Pe=@qd7zjR+6{B0vO)01+SpM1Tko0U|&IhyW2dI|&3FVl&VG#A7}e
z#HNPaj{q+2T-e#;1FszXHO61#Z=anGlGQ|j2oM1xKm>>Y5g-CYfCvx)B0vO)K%PKA
zY!>pL0yunPla%`q0KfllkzQi_IsP5Kg)ipL&@V$@3at&b29E@P9Q=In+Tg6f%Yh#R
zwgfu-@Awb+@Aj|vU*LPo_k?ekPxpP?7Y1e8hyW2F0z`la5CI}U1U@7L=C_D!W7mM5
zP`frMYS-XMzc#F<%)$;Fz{4)^M)<C1bO`RBo!2a~8%t@#v{)Cs7Nnr?{8?qH$HIjs
zu==J7s*ewX9u%H4Q)KHZumy_FZ4}x1az)Lp{VJ60o*^==Tv>HUO`4IEnU0NsN>78x
z5~WnqrWQ|iMbzP_I%qvs0F)Ngi)?c#B|{z54YO-VOQn^#b~D<xpsq5VR2puiHf){l
z+A{v++KEK8q-I7U(9FUbk!>kuTA|ar66v@J_iYb?+VW~A8<i^6(_?Bp0>Ncy!x1~~
zu<6BBBD<x+Dn0>EB}>PZ)ryf)4Wtt;L#N2>epMU7_AcQf8>zHCYjj)%4r{1-VMt_p
zg{u9~ms$$7x`So9!;Y1hIux<?K&2xfGNn{Qwo1PEe1Dl?XBRK(?~8pROP6Y(2;30d
zc|KuTk5ne`ixdrbut~}UcY3^)EA^{-+`@LTr*abujY4a^ZU{<Oh0euPQ;ov7hKK9G
z^M&B?A67rRD(~l*3eS{CqSiuZ<zuL*E4sxp;h`FaZophp7_FTSkwr^6h|_FQ`ME-w
z$*`y_e*fPh9cKJ@{2|}{zAE1m-!b2y?;tk<uLKN!sdSjTLI*;32KIzL9c&JC`u`o0
zgGc?3`|t7J5Zvc)3$6)%;X`t<BQuBq5g-CYfCvx)B0vNf>lQZ(rQSIf-bv@y|1Is}
zjU_eky(YY@^?!4@>alR433$^ytCH&T>;I;?;<`#~fqG}o6W5n5np^)jE)cb{Whbxy
zXIv;IN~+{mT$9%S4L#!Kl1fvo|LZH%DPI5AmGkHH^?&U};+B%86*{eW{a<sLb5yEW
zZ|e1b^`+u1l~!@;`oF3iog%aI>;Ld#airq*Ca?ebLQ$_&wXpsVmFo^WR*KjE!A?;r
z*^sSbNeK2pxnY-!#=;8Hzer4%Y@i6-r1ihAT;LZePOSf>Gf<jX|9j8C#6qJ$!?Vme
z>Z;VaDcAq*GwA2S`rma1W2k7s?JO4_r?3CTE-_lt!Bf`%j=AFclAqEQnM2S2qqhk}
zfCvx)B0vO)01+SpM1Tko0U|&I&S?T7{7fK--~S8z1cN`?hyW2F0z`la5CI}U1c(3;
zAOb{y2oQmDkbo!#CH(%s-v2qq-{P<HKk`F-FFXfeP3Xeld%;(Ne+cdg-V<COI2L#=
zFcA1?z~w*T-v@<hBLYN#2oM1xKm>>Y5g-CY;C~>nu-++1@%o!t=2g%>wrx@sV_8qp
zv2%Q>eKfXBSW)K^q)}hp(1Te=QIRQjLlb%SUN0WwzbA)H@YF)<;o9?@%{8#|v6}rg
z_KD-7o#v!-c{MiYtp58EZf8CEz>@c0*!I#lmfz}Nng6tAZ-s2h3IFn3W(PyE2NK3s
zJ0J<eOBk<VW3EkASjXC`&N0ls;rc&iPy7X&IWcq_NOpX4BPPdQd>d`Z=m{5N@tSdH
zPwo)ep@?n$J--@zaDQ;Nz&1E?J`6Onm7uV4=k?-@%-0}u*da(QVQ^<jcy7Us%*pJk
z@MAaKz3u#wwjNRCXv-)0ODER>(>G@x$-dr-t!AJ7(GteSvKQUCM0hiM?V4Ej#KZSz
zdG_)g<U`xctLHLdw*9}|F;fULhs>P4*@>n?EOwvBcF&MlxF-Z=421SCvQEC16+i)q
zoVMA$Z2iCZj4&$*N{fR-jTt8!pQtt5%VaO}I_C!9>~(>mccHG`Pd$P0lKIJ3n<Yq|
zNPfmWW1RhL*E;FX*o4>Zhq596_Ie+roxVrjeN~9ff4)Q5fM?k_UnD`kA-yB+{TJ4Q
z+gQPm6G#MdtRw@?1>OWZH(;J)`z$X2+_A{}i^ILW2QjO7JCt3bd~@5r1NWUgx^&f%
zuV-B^?D!3I-|;<f-OIe04*4sOE^P*TeXJL?GTE!%i$iBawJ4Y{{=--Io<zs>RtuPw
zySreRbzCG4=3=ra3wL`bj^&9xmdUg+_Quw=KlQX@rr_;mjzdce=fNCy%yzsu6ZbRl
zl3DLV=ivd6*z0n$UU2L1M>6fQkcs*Q48CQK@oIUM2g9qy^Vg%<scpk5@>_;wpMHAw
zyHof#o~>K`y~}sX;jqsQJyq|%I>4B_#RdCIT#M~!`t1Iel^6Z`$&4zri&?ga9rs91
zIA81BxpL=@kO;diVt;OPxE*jDatL6mB+Ybsx691_{=c3NF#bM2#{b5T@WcFh{tSPL
zAK*Xd`}iLIFyGC;#dq>A^Bw#)KE!n>OdAm(0z`la5CI}U1c(3;AOb{y2oM1xFeL$8
z!VC5FR#I1ICAGCyQd46k)zwx~Rb?gNu$6FbC83a&1cO!*2v~{VZzVqK(*Z$}tRt`2
zN<1Dbal5U=<+2i|(@I3qN*oR+d<<ZH|3B4-Qd2rXnM8mH5CI}U1c(3;AOb{y2oM1x
zKm>@uSxErz|EKl;S=n`DG!Y;IM1Tko0U|&IhyW2F0z`la5P^yb;QjwCd=viu{~i8M
z{suqFU*RvpJpj-0-|{E<<NPuH2;2v-mw%ssmw%h@;`hS60C({(^4t06_y{-Qet=CC
zH%0}C01+SpM1Tko0U|&IhyW2F0z`laobv?47I@sg*o<NpiY62@Q8c2MfuaFLJ&HOM
uwJ2&(RHLXu5k|pLgir)g1W@=<_)thFyeK><+$dZqoG3&T4ixzO|NjDFp32_<

diff --git a/pyfixest/estimation/feiv_.py b/pyfixest/estimation/feiv_.py
index a5a32173..1f7fa803 100644
--- a/pyfixest/estimation/feiv_.py
+++ b/pyfixest/estimation/feiv_.py
@@ -33,6 +33,8 @@ class Feiv(Feols):
         Names of the coefficients of Z.
     collin_tol : float
         Tolerance for collinearity check.
+    solver: str, default is 'np.linalg.solve'
+        Solver to use for the estimation. Alternative is 'np.linalg.lstsq'.
     weights_name : Optional[str]
         Name of the weights variable.
     weights_type : Optional[str]
@@ -140,38 +142,13 @@ def __init__(
         self._supports_cluster_causal_variance = False
         self._endogvar = endogvar
 
-    def solve_ols(self, tZX, tZY, solver):
-        """
-        Solve the ordinary least squares problem using the specified solver.
-
-        Parameters
-        ----------
-        tZX (array-like): The design matrix.
-        tZY (array-like): The response variable.
-        solver (str): The solver to use. Supported solvers are "np.linalg.lstsq"
-        and "np.linalg.solve".
-
-        Returns
-        -------
-        array-like: The solution to the ordinary least squares problem.
-
-        Raises
-        ------
-        ValueError: If the specified solver is not supported.
-        """
-        if solver == "np.linalg.lstsq":
-            return np.linalg.lstsq(tZX, tZY, rcond=None)[0]
-        elif solver == "np.linalg.solve":
-            return np.linalg.solve(tZX, tZY).flatten()
-        else:
-            raise ValueError(f"Solver {solver} not supported.")
-
     def get_fit(self) -> None:
         """Fit a IV model using a 2SLS estimator."""
         _X = self._X
         _Z = self._Z
         _Y = self._Y
         _endogvar = self._endogvar
+        _solver = self._solver
 
         #  Start First Stage
 
@@ -203,10 +180,7 @@ def get_fit(self) -> None:
         H = self._tXZ @ self._tZZinv
         A = H @ self._tZX
         B = H @ self._tZy
-        solver = self._solver
-
-        self._beta_hat = self.solve_ols(A, B, solver)
-
+        self._beta_hat = self.solve_ols(A, B, _solver)
 
         # Predicted values and residuals
         self._Y_hat_link = self._X @ self._beta_hat
diff --git a/pyfixest/estimation/feols_.py b/pyfixest/estimation/feols_.py
index caed86cf..bb1dac71 100644
--- a/pyfixest/estimation/feols_.py
+++ b/pyfixest/estimation/feols_.py
@@ -62,7 +62,7 @@ class Feols:
     weights_type : Optional[str]
         Type of the weights variable. Either "aweights" for analytic weights or
         "fweights" for frequency weights.
-     solver : str, optional.
+    solver : str, optional.
         The solver to use for the regression. Can be either "np.linalg.solve" or
         "np.linalg.lstsq". Defaults to "np.linalg.solve".
 
@@ -171,7 +171,8 @@ class Feols:
         Adjusted R-squared value of the model.
     _adj_r2_within : float
         Adjusted R-squared value computed on demeaned dependent variable.
-
+    _solver: str
+        The solver used to fit the normal equation.
     """
 
     def __init__(
@@ -301,6 +302,32 @@ def __init__(
         self.summary = functools.partial(_tmp, models=[self])
         self.summary.__doc__ = _tmp.__doc__
 
+    def solve_ols(self, tZX: np.ndarray, tZY: np.ndarray, solver: str):
+        """
+        Solve the ordinary least squares problem using the specified solver.
+
+        Parameters
+        ----------
+        tZX (array-like): Z'X.
+        tZY (array-like): Z'Y.
+        solver (str): The solver to use. Supported solvers are "np.linalg.lstsq"
+        and "np.linalg.solve".
+
+        Returns
+        -------
+        array-like: The solution to the ordinary least squares problem.
+
+        Raises
+        ------
+        ValueError: If the specified solver is not supported.
+        """
+        if solver == "np.linalg.lstsq":
+            return np.linalg.lstsq(tZX, tZY, rcond=None)[0]
+        elif solver == "np.linalg.solve":
+            return np.linalg.solve(tZX, tZY).flatten()
+        else:
+            raise ValueError(f"Solver {solver} not supported.")
+
     def get_fit(self) -> None:
         """
         Fit an OLS model.
@@ -322,9 +349,7 @@ def get_fit(self) -> None:
         else:  # Default to np.linalg.solve
             self._beta_hat = np.linalg.solve(self._tZX, self._tZy).flatten()
 
-        # self._beta_hat, _, _, _ = lstsq(self._tZX, self._tZy, lapack_driver='gelsy')
-
-        # self._beta_hat = (self._tZXinv @ self._tZy).flatten()
+        self._beta_hat = self.solve_ols(self._tZX, self._tZy, _solver)
 
         self._Y_hat_link = self._X @ self._beta_hat
         self._u_hat = self._Y.flatten() - self._Y_hat_link.flatten()
diff --git a/pyfixest/estimation/fepois_.py b/pyfixest/estimation/fepois_.py
index 60f28166..3f6d11da 100644
--- a/pyfixest/estimation/fepois_.py
+++ b/pyfixest/estimation/fepois_.py
@@ -48,8 +48,11 @@ class Fepois(Feols):
         Maximum number of iterations for the IRLS algorithm.
     tol : Optional[float], default=1e-08
         Tolerance level for the convergence of the IRLS algorithm.
+    solver: str, default is 'np.linalg.solve'
+        Solver to use for the estimation. Alternative is 'np.linalg.lstsq'.
     fixef_tol: float, default = 1e-08.
         Tolerance level for the convergence of the demeaning algorithm.
+    solver:
     weights_name : Optional[str]
         Name of the weights variable.
     weights_type : Optional[str]
@@ -115,32 +118,6 @@ def __init__(
         self.deviance = None
         self._Xbeta = np.array([])
 
-    def solve_ols(self, tZX, tZY, solver):
-        """
-        Solve the ordinary least squares problem using the specified solver.
-
-        Parameters
-        ----------
-        tZX (array-like): The design matrix.
-        tZY (array-like): The response variable.
-        solver (str): The solver to use. Supported solvers are "np.linalg.lstsq"
-        and "np.linalg.solve".
-
-        Returns
-        -------
-        array-like: The solution to the ordinary least squares problem.
-
-        Raises
-        ------
-        ValueError: If the specified solver is not supported.
-        """
-        if solver == "np.linalg.lstsq":
-            return np.linalg.lstsq(tZX, tZY, rcond=None)[0]
-        elif solver == "np.linalg.solve":
-            return np.linalg.solve(tZX, tZY).flatten()
-        else:
-            raise ValueError(f"Solver {solver} not supported.")
-
     def get_fit(self) -> None:
         """
         Fit a Poisson Regression Model via Iterated Weighted Least Squares (IWLS).
@@ -179,7 +156,7 @@ def get_fit(self) -> None:
         _iwls_maxiter = 25
         _tol = self.tol
         _fixef_tol = self.fixef_tol
-        solver = self._solver
+        _solver = self._solver
 
         def compute_deviance(_Y: np.ndarray, mu: np.ndarray):
             with warnings.catch_warnings():
@@ -244,7 +221,7 @@ def compute_deviance(_Y: np.ndarray, mu: np.ndarray):
             XWX = WX.transpose() @ WX
             XWZ = WX.transpose() @ WZ
 
-            delta_new = self.solve_ols(XWX, XWZ, solver)  # eq (10), delta_new -> reg_z
+            delta_new = self.solve_ols(XWX, XWZ, _solver)  # eq (10), delta_new -> reg_z
             resid = Z_resid - X_resid @ delta_new
 
             mu_old = mu.copy()

From ed5c7fde7421cb3164ac7dbff40a7722d80c92ad Mon Sep 17 00:00:00 2001
From: Alexander Fischer <alexander-fischer1801@t-online.de>
Date: Sun, 30 Jun 2024 22:22:07 +0200
Subject: [PATCH 6/9] fix tests except for Attribute error Fepois has no
 attribute _N

---
 .coverage                     | Bin 53248 -> 53248 bytes
 pyfixest/estimation/feols_.py |   8 +-------
 tests/test_solver.py          |  27 +++++++++++++++------------
 3 files changed, 16 insertions(+), 19 deletions(-)

diff --git a/.coverage b/.coverage
index dc0b802595222df324e5f93ee87a4cf6b6fc263b..27469ab0cbae0b2e3e3351b16bfacd8744a48697 100644
GIT binary patch
delta 372
zcmZozz}&Eac>`O6UM~axXa3jxxA`ydZ{T0TAH`?S$IAPb_bKmA-g&$!JpXy_^Az(K
z^N8?pb06XE-7F{&$-Q|-cMzk3oB#_Wr=5TnGXos_{a-)t$Jbv!!Y6${!!wzsH`U65
z52&b!Z<_^JmBau0^<SooyklZuc)$qP1!5e4aDfsG3<(Sj96?-@=k<y+9pIb1y?0gk
z_S5@+J@)<odA|DRlmD3*K5+bq%q;^cMFkEF3}6FTfWitODFG0{zyTt@$b(1-8_Z((
zD?iaug2|z7qNA7<SQ}Uv$yqEqFb;_3aA05(X!r*c0xM=&@E@pK31lB5gO(^r&VqrV
aLDztRp@@OuM;v2{J10n?!se*{WexyKwN9D<

delta 340
zcmZozz}&Eac>`O6-V6r*&-}0XZ}VT^-@w0wKZ?(ukCpc=?^E8Lyz_W7d0BW~@zn9y
z@yPRtai8a&u~|?cg?n>vcMzkLga8X8r;&gbGXos_{a+7cGfn>AlPYb@2b9j>+hzfl
zclduFB*nl4A|}`ON==^EE6(J=FnN3Ls&ICYAPlfDFem`o2N+=j2%3?B1H=&k;sy}I
z0YrdRGh6_h4q<~u8NSp{bd+TJ$vDwb)EcY}EDSUq1_T&j91yL*z{ugi@C(H$0uBF|
cf!dWo7J}4(@B<(Q0hs#~K%9onQT@vt0DW*a#sB~S

diff --git a/pyfixest/estimation/feols_.py b/pyfixest/estimation/feols_.py
index bb1dac71..0fa7beb7 100644
--- a/pyfixest/estimation/feols_.py
+++ b/pyfixest/estimation/feols_.py
@@ -322,7 +322,7 @@ def solve_ols(self, tZX: np.ndarray, tZY: np.ndarray, solver: str):
         ValueError: If the specified solver is not supported.
         """
         if solver == "np.linalg.lstsq":
-            return np.linalg.lstsq(tZX, tZY, rcond=None)[0]
+            return np.linalg.lstsq(tZX, tZY, rcond=None)[0].flatten()
         elif solver == "np.linalg.solve":
             return np.linalg.solve(tZX, tZY).flatten()
         else:
@@ -343,12 +343,6 @@ def get_fit(self) -> None:
         self._tZX = _Z.T @ _X
         self._tZy = _Z.T @ _Y
 
-        # self._tZXinv = np.linalg.inv(self._tZX)
-        if _solver == "np.linalg.lstsq":
-            self._beta_hat, _, _, _ = np.linalg.lstsq(self._tZX, self._tZy, rcond=None)
-        else:  # Default to np.linalg.solve
-            self._beta_hat = np.linalg.solve(self._tZX, self._tZy).flatten()
-
         self._beta_hat = self.solve_ols(self._tZX, self._tZy, _solver)
 
         self._Y_hat_link = self._X @ self._beta_hat
diff --git a/tests/test_solver.py b/tests/test_solver.py
index 826f436a..3fca3c15 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -9,22 +9,25 @@
 
 @pytest.fixture()
 def data_feols():
-    data = pf.get_data()
-    Y = data["Y"].to_numpy()
+    data = pf.get_data().dropna()
+    Y = data["Y"].to_numpy().reshape((-1, 1))
     X = data[["X1", "X2"]].to_numpy()
-    Z = data[["Z1"]].to_numpy()
-    weights = data["weights"].to_numpy()  # Define the variable "weights"
+    Z = data[["Z1", "X2"]].to_numpy()
+    weights = (
+        data["weights"].to_numpy().reshape((-1, 1))
+    )  # Define the variable "weights"
     return X, Y, Z, weights
 
 
 @pytest.fixture()
 def data_fepois():
-    data = pf.get_data(model="Fepois")
-    Y = data["Y"].to_numpy()
+    data = pf.get_data(model="Fepois").dropna()
+    Y = data["Y"].to_numpy().reshape((-1, 1))
     X = data[["X1", "X2"]].to_numpy()
-    Z = data[["Z1"]].to_numpy()
-    weights = data["weights"].to_numpy()  # Define the variable "weights"
-    return X, Y, Z, weights
+    weights = (
+        data["weights"].to_numpy().reshape((-1, 1))
+    )  # Define the variable "weights"
+    return X, Y, weights
 
 
 def test_solver_fepois(data_fepois):
@@ -35,7 +38,7 @@ def test_solver_fepois(data_fepois):
     and np.linalg.solve. It asserts that
     the results are identical or very close within a tolerance.
     """
-    X, Y, Z, weights = data_fepois
+    X, Y, weights = data_fepois
     obj = Fepois(
         X=X,
         Y=Y,
@@ -76,7 +79,7 @@ def test_solver_feiv(data_feols):
         Y=Y,
         X=X,
         Z=Z,
-        endogvar=X[:, 0],
+        endogvar=X[:, 0].reshape((-1, 1)),
         weights=weights,
         coefnames_x=["X1", "X2"],
         collin_tol=1e-08,
@@ -110,7 +113,7 @@ def test_solver_feols(data_feols):
     and np.linalg.solve. It asserts that
     the results are identical or very close within a tolerance.
     """
-    X, Y, Z, weights = data_feols
+    X, Y, _, weights = data_feols
 
     obj = Feols(
         X=X,

From 0c75437d8617e394918fb71c1bfd35d0ce0c38cc Mon Sep 17 00:00:00 2001
From: Alexander Fischer <alexander-fischer1801@t-online.de>
Date: Tue, 2 Jul 2024 07:29:51 +0200
Subject: [PATCH 7/9] Update fepois_.py

add reshape to delta_new
---
 pyfixest/estimation/fepois_.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyfixest/estimation/fepois_.py b/pyfixest/estimation/fepois_.py
index 3f6d11da..1db646cf 100644
--- a/pyfixest/estimation/fepois_.py
+++ b/pyfixest/estimation/fepois_.py
@@ -221,7 +221,7 @@ def compute_deviance(_Y: np.ndarray, mu: np.ndarray):
             XWX = WX.transpose() @ WX
             XWZ = WX.transpose() @ WZ
 
-            delta_new = self.solve_ols(XWX, XWZ, _solver)  # eq (10), delta_new -> reg_z
+            delta_new = self.solve_ols(XWX, XWZ, _solver).reshape((-1,1))  # eq (10), delta_new -> reg_z
             resid = Z_resid - X_resid @ delta_new
 
             mu_old = mu.copy()

From e7d980e59d4cfdec7c038c86d62783c3f036fbac Mon Sep 17 00:00:00 2001
From: Alexander Fischer <alexander-fischer1801@t-online.de>
Date: Tue, 2 Jul 2024 21:55:55 +0200
Subject: [PATCH 8/9] fix tests

---
 .coverage                      | Bin 53248 -> 53248 bytes
 pyfixest/estimation/fepois_.py |   4 +++-
 tests/test_solver.py           |  20 +++++++++++++++++---
 3 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/.coverage b/.coverage
index 27469ab0cbae0b2e3e3351b16bfacd8744a48697..153b61681dac2c386ad3833811fd700daeaf44ee 100644
GIT binary patch
delta 279
zcmZozz}&Eac>`O6ULOPhXa3jxxA`ydZ{T0TAI0au$Hx1P_Zjal-ub+#ybL@KcuIIo
zctm-4xQ}x8Z59-W;@;fb{gILX|G#gH^*kaMa<XIEFYUfPxwm&SyEq3cBWH`=WcNND
zm#sG$8884N0|!V<0Einv3<nUwpa3E+)PqO}8_Z((Qa{mAlIbVoL`P9;ur{zT3q%q`
z3VeWZK(qn_BZmXSFLs1HP({N(W}tQ@kcBJ^D?&i(9su#S0}Ko&7#SW^vrdpz1@ajh
KHb?a@a{vHh3pxY<

delta 279
zcmZozz}&Eac>`O6UM~axXa3jxxA`ydZ{T0TAH`?S$IAPb_bKmA-g&$!JpXy_^Az(K
z^N8?pb06XE-7F{&$-TL+`y(SiBLgxxFuAXHGrI%_D<fx%!DNp<9oNX*GGql1wgUqL
z6NtqE#0nsW0El4V01;p0K_r9?W-<JgpXeyT<WM)!QOpXg4J?dgF^dk21EM(`7?=ba
x{=tO6idh!?2kKS=*~iGBB?^+WU|?v_HDF*UVqo|Y$C%>I2~w!AIjVn|0|2RjE;IlD

diff --git a/pyfixest/estimation/fepois_.py b/pyfixest/estimation/fepois_.py
index 3f6d11da..9ba44092 100644
--- a/pyfixest/estimation/fepois_.py
+++ b/pyfixest/estimation/fepois_.py
@@ -221,7 +221,9 @@ def compute_deviance(_Y: np.ndarray, mu: np.ndarray):
             XWX = WX.transpose() @ WX
             XWZ = WX.transpose() @ WZ
 
-            delta_new = self.solve_ols(XWX, XWZ, _solver)  # eq (10), delta_new -> reg_z
+            delta_new = self.solve_ols(XWX, XWZ, _solver).reshape(
+                (-1, 1)
+            )  # eq (10), delta_new -> reg_z
             resid = Z_resid - X_resid @ delta_new
 
             mu_old = mu.copy()
diff --git a/tests/test_solver.py b/tests/test_solver.py
index 3fca3c15..fd20ff2a 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -46,7 +46,7 @@ def test_solver_fepois(data_fepois):
         coefnames=["X1", "X2"],
         collin_tol=1e-08,
         weights_name="weights",
-        weights_type=None,
+        weights_type="aweights",
         solver="np.linalg.solve",
         fe=None,
         drop_singletons=False,
@@ -55,7 +55,21 @@ def test_solver_fepois(data_fepois):
     obj._solver = "np.linalg.lstsq"
     obj.get_fit()
     beta_hat_lstsq = obj._beta_hat.copy()
+
     # Use np.linalg.solve solver
+    X, Y, weights = data_fepois
+    obj = Fepois(
+        X=X,
+        Y=Y,
+        weights=weights,
+        coefnames=["X1", "X2"],
+        collin_tol=1e-08,
+        weights_name="weights",
+        weights_type="aweights",
+        solver="np.linalg.solve",
+        fe=None,
+        drop_singletons=False,
+    )
     obj._solver = "np.linalg.solve"
     obj.get_fit()
     beta_hat_solve = obj._beta_hat.copy()
@@ -84,7 +98,7 @@ def test_solver_feiv(data_feols):
         coefnames_x=["X1", "X2"],
         collin_tol=1e-08,
         weights_name="weights",
-        weights_type=None,
+        weights_type="aweights",
         solver="np.linalg.solve",
         coefnames_z=["Z1"],
     )
@@ -122,7 +136,7 @@ def test_solver_feols(data_feols):
         collin_tol=1e-08,
         coefnames=["X1", "X2"],
         weights_name="weights",
-        weights_type=None,
+        weights_type="aweights",
         solver="np.linalg.solve",
     )
 

From f5c2eaee69b36b5156400f730ab2ddef77609e19 Mon Sep 17 00:00:00 2001
From: Alexander Fischer <alexander-fischer1801@t-online.de>
Date: Tue, 2 Jul 2024 22:03:04 +0200
Subject: [PATCH 9/9] add scipy.sparse.linalg.lsqr solver

---
 .coverage                     | Bin 53248 -> 53248 bytes
 pyfixest/estimation/feols_.py |   8 +++++---
 tests/test_solver.py          |   9 +++++++++
 3 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/.coverage b/.coverage
index 153b61681dac2c386ad3833811fd700daeaf44ee..1658c0a5e6ae67d93dfaae938c5c9bb7c37ca102 100644
GIT binary patch
delta 297
zcmZozz}&Eac>`O6UM~axXa3jxxA`ydZ{T0TAH`?S$IAPb_bKmA-g&$!JpXy_^Az(K
z^N8?pb06XE-7F{&$-TL+`xm45>Oxts_j&9LAix148W<TAJTxcI?_I<z!NJPN*<!%4
z@4kb>jmds}D(=b;uBSoNf|SAm3xfiP?EoT}Km;R*XaEro{}F5upMl}Rf5wT95=>9n
zCpwB+D}aSSgaDky;BY_~$^i*BFt8{vFussS$b)S9#Shf3>;M!LU|0|&2vWz)z;J_!
XnL&Y_q2ax#L)jFt2BytX{mUEx+!Q&?

delta 295
zcmZozz}&Eac>`O6ULOPhXa3jxxA`ydZ{T0TAI0au$Hx1P_Zjal-ub+#ybL@KcuIIo
zctm-4xQ}x8Z59-W;@;fb{fkjNJEr~8?%RwEAixA75*Qdbg19En>s`bv&cVva*`mj>
z@4kb>iOIfwDsCJ<w%&xO1u2CCMg|TLTL6d~Knw>E!Jq&lF4Tia2ph~|_)<U7QIhE=
z<3vYMYp^!3FbhNyL<)R>aX_>J10#n6!!LG(JWxf$KW3nIC6I+I3@buF>K*{`wgU_d
XCm0zXRI^TyRR!`H8a7AuFLM9@-zz!*

diff --git a/pyfixest/estimation/feols_.py b/pyfixest/estimation/feols_.py
index 0fa7beb7..0ec111ff 100644
--- a/pyfixest/estimation/feols_.py
+++ b/pyfixest/estimation/feols_.py
@@ -8,7 +8,7 @@
 import pandas as pd
 from formulaic import Formula
 from scipy.sparse import csr_matrix
-from scipy.sparse.linalg import spsolve
+from scipy.sparse.linalg import lsqr, spsolve
 from scipy.stats import f, norm, t
 
 from pyfixest.errors import VcovTypeNotSupportedError
@@ -310,8 +310,8 @@ def solve_ols(self, tZX: np.ndarray, tZY: np.ndarray, solver: str):
         ----------
         tZX (array-like): Z'X.
         tZY (array-like): Z'Y.
-        solver (str): The solver to use. Supported solvers are "np.linalg.lstsq"
-        and "np.linalg.solve".
+        solver (str): The solver to use. Supported solvers are "np.linalg.lstsq",
+        "np.linalg.solve", and "scipy.sparse.linalg.lsqr".
 
         Returns
         -------
@@ -325,6 +325,8 @@ def solve_ols(self, tZX: np.ndarray, tZY: np.ndarray, solver: str):
             return np.linalg.lstsq(tZX, tZY, rcond=None)[0].flatten()
         elif solver == "np.linalg.solve":
             return np.linalg.solve(tZX, tZY).flatten()
+        elif solver == "scipy.sparse.linalg.lsqr":
+            return lsqr(tZX, tZY)[0].flatten()
         else:
             raise ValueError(f"Solver {solver} not supported.")
 
diff --git a/tests/test_solver.py b/tests/test_solver.py
index fd20ff2a..14f69f3f 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -150,7 +150,16 @@ def test_solver_feols(data_feols):
     obj.get_fit()
     beta_hat_solve = obj._beta_hat.copy()
 
+    # scipy.sparse.linalg.lsqr solver
+    obj._solver = "scipy.sparse.linalg.lsqr"
+    obj.get_fit()
+    beta_hat_lsqr = obj._beta_hat.copy()
+
     # Assert that the results are identical or very close within a tolerance
     np.testing.assert_array_almost_equal(
         beta_hat_lstsq, beta_hat_solve, decimal=6, err_msg="Solvers' results differ"
     )
+
+    np.testing.assert_array_almost_equal(
+        beta_hat_lstsq, beta_hat_lsqr, decimal=6, err_msg="Solvers' results differ"
+    )