diff --git a/linearmodels/panel/covariance.py b/linearmodels/panel/covariance.py
index 880451024d..701234971c 100644
--- a/linearmodels/panel/covariance.py
+++ b/linearmodels/panel/covariance.py
@@ -105,6 +105,13 @@ def __init__(
         self._nobs_eff = self._nobs - extra_df
         if debiased:
             self._nobs_eff -= self._nvar
+        if self._nobs_eff <= 0:
+            regressors = f"{self._nvar} regressors and " if debiased else ""
+            raise ValueError(
+                "The covariance of the model cannot be estimated since the model has "
+                f"no residual degrees of freedom. It has {self._nobs} observations "
+                f"and uses {regressors}{extra_df} degrees of freedom for the effects."
+            )
         self._scale = self._nobs / self._nobs_eff
         self._name = "Unadjusted"
 
diff --git a/linearmodels/panel/data.py b/linearmodels/panel/data.py
index aecca87421..a39f629729 100644
--- a/linearmodels/panel/data.py
+++ b/linearmodels/panel/data.py
@@ -442,6 +442,9 @@ def _demean_both(self, weights: PanelData | None) -> PanelData:
         e = self.demean(group, weights=weights)
         assert isinstance(e, PanelData)
         d = self.dummies(dummy, drop_first=True)
+        if d.shape[1] == 0:
+            # A single entity or period: there are no dummy variables to remove
+            return e
         d.index = e.index
         d_pd = PanelData(d).demean(group, weights=weights)
         d_arr = d_pd.values2d
diff --git a/linearmodels/panel/model.py b/linearmodels/panel/model.py
index a63ca58938..102b16cc10 100644
--- a/linearmodels/panel/model.py
+++ b/linearmodels/panel/model.py
@@ -523,9 +523,10 @@ def _grid_to_rows(self, grid: PanelData, labelled: bool) -> PanelData:
         return PanelData(rows)
 
     def _check_exog_rank(self) -> int:
-        if not self._check_rank:
-            return self.exog.shape[1]
         x = cast("linearmodels.typing.data.Float64Array", self.exog.values2d)
+        if not self._check_rank:
+            # The shape of the PanelData is (variables, periods, entities)
+            return x.shape[1]
         _, _, rank_of_x, _ = _lstsq(x, np.ones(x.shape[0]))
         if rank_of_x < x.shape[1]:
             raise ValueError(
@@ -574,8 +575,9 @@ def _validate_data(self) -> None:
             self._not_null = cast("BoolArray", np.asarray(~missing))
 
         w_df = self.weights.dataframe
-        if np.any(np.asarray(w_df) <= 0):
-            raise ValueError("weights must be strictly positive.")
+        w_values = np.asarray(w_df)
+        if np.any(w_values <= 0) or np.any(np.isinf(w_values)):
+            raise ValueError("weights must be strictly positive and finite.")
         w_df = w_df / w_df.mean()
         self.weights = PanelData(w_df)
         rank_of_x = self._check_exog_rank()
@@ -607,6 +609,10 @@ def _f_statistic(
         weps_const = y
         num_df = x.shape[1]
         name = "Model F-statistic (homoskedastic)"
+        if df_resid <= 0:
+            return InvalidTestStatistic(
+                "The model has no residual degrees of freedom", name=name
+            )
         if self.has_constant:
             if num_df == 1:
                 return InvalidTestStatistic("Model contains only a constant", name=name)
@@ -2205,17 +2211,29 @@ def fit(
                 df_num -= 1
             weps_pooled = wy - wx @ _lstsq(wx, wy, rcond=None)[0]
             resid_ss_pooled = float(np.squeeze(weps_pooled.T @ weps_pooled))
-            num = (resid_ss_pooled - resid_ss) / df_num
-
-            denom = resid_ss / df_denom
-            stat = num / denom
-            f_pooled = WaldTestStatistic(
-                stat,
-                "Effects are zero",
-                df_num,
-                df_denom=df_denom,
-                name="Pooled F-statistic",
+            name = "Pooled F-statistic"
+            f_pooled: (
+                WaldTestStatistic | InapplicableTestStatistic | InvalidTestStatistic
             )
+            if df_num <= 0:
+                # A single entity (or period) with a constant: the effect is the constant
+                f_pooled = InapplicableTestStatistic(
+                    reason="The effects do not use any degrees of freedom", name=name
+                )
+            elif df_denom <= 0:
+                f_pooled = InvalidTestStatistic(
+                    "The model has no residual degrees of freedom", name=name
+                )
+            else:
+                num = (resid_ss_pooled - resid_ss) / df_num
+                denom = resid_ss / df_denom
+                f_pooled = WaldTestStatistic(
+                    num / denom,
+                    "Effects are zero",
+                    df_num,
+                    df_denom=df_denom,
+                    name=name,
+                )
             res.update(f_pooled=f_pooled)
             effects = DataFrame(
                 eps_effects - eps,
@@ -2421,10 +2439,20 @@ def fit(
 
         wx = root_w * x
         wy = root_w * y
-        params = _lstsq(wx, wy, rcond=None)[0]
+        params, _, rank, _ = _lstsq(wx, wy, rcond=None)
+        if self._check_rank and rank < x.shape[1]:
+            raise ValueError(
+                "The regressors are collinear in the means of the entities, which "
+                f"are the {y.shape[0]} observations of the regression, so that the "
+                "parameters are not identified. This happens if there are fewer "
+                "entities than regressors, or if a regressor is a combination of "
+                "the others only in the means, e.g., the dummy variables of the "
+                "periods and the constant in a balanced panel. If you wish to "
+                "proceed, you can set check_rank=False."
+            )
 
         df_resid = y.shape[0] - x.shape[1]
-        df_model = (x.shape[1],)
+        df_model = x.shape[1]
         nobs = y.shape[0]
         cov_config = self._setup_clusters(cov_config)
         extra_df = 0
@@ -3311,17 +3339,16 @@ def _mr(ex: linearmodels.typing.data.Float64Array) -> int:
         if valid_blocks.sum() < exog.shape[1]:
 
             warnings.warn(
-                "The number of time-series observation available to estimate "
-                "cross-sectional\nregressions, {}, is less than the number of "
-                "parameters in the model. Parameter\ninference is not "
-                "available.".format(valid_blocks.sum()),
+                "The number of cross-sectional regressions that can be estimated, "
+                "{}, is less than the number of parameters in the model. Parameter "
+                "inference is not available.".format(valid_blocks.sum()),
                 InferenceUnavailableWarning,
                 stacklevel=3,
             )
         elif valid_blocks.sum() < valid_blocks.shape[0]:
 
             warnings.warn(
-                "{} of the time-series regressions cannot be estimated due to "
+                "{} of the cross-sectional regressions cannot be estimated due to "
                 "deficient rank.".format(valid_blocks.shape[0] - valid_blocks.sum()),
                 MissingValueWarning,
                 stacklevel=3,
diff --git a/linearmodels/panel/utility.py b/linearmodels/panel/utility.py
index 6f2e47cbd4..8e0470aa34 100644
--- a/linearmodels/panel/utility.py
+++ b/linearmodels/panel/utility.py
@@ -642,6 +642,17 @@ class PanelModelData(NamedTuple):
     clusters: DataFrame
 
 
+def _integers(
+    rng: numpy.random.RandomState | numpy.random.Generator,
+    high: int,
+    size: tuple[int, ...],
+) -> linearmodels.typing.data.IntArray:
+    """Random integers in [0, high) from a RandomState or a Generator"""
+    if isinstance(rng, numpy.random.Generator):
+        return rng.integers(0, high, size)
+    return rng.randint(0, high, size)
+
+
 def generate_panel_data(
     nentity: int = 971,
     ntime: int = 7,
@@ -650,7 +661,7 @@ def generate_panel_data(
     missing: float = 0,
     other_effects: int = 2,
     ncats: int | list[int] = 4,
-    rng: numpy.random.RandomState | None = None,
+    rng: numpy.random.RandomState | numpy.random.Generator | None = None,
 ) -> PanelModelData:
     """
     Simulate panel data for testing
@@ -673,9 +684,9 @@ def generate_panel_data(
         The number of categories to use in other_effects and variance
         clusters. If list-like, then it must have as many elements
         as other_effects.
-    rng : RandomState
-        A NumPy RandomState instance. If not provided, one is initialized
-        using a fixed seed.
+    rng : RandomState or Generator
+        A NumPy RandomState or Generator instance. If not provided, a
+        RandomState is initialized using a fixed seed.
 
     Returns
     -------
@@ -721,12 +732,11 @@ def generate_panel_data(
         _c = []
         for i in range(other_effects):
             nc = ncats[i]
-            _c.append(rng.randint(0, nc, (1, t, n)))
+            _c.append(_integers(rng, nc, (1, t, n)))
         c = np.concatenate(_c, 0)
 
     vcats = [f"varcat.{i}" for i in range(2)]
-    vc2 = np.ones((2, t, 1)) @ rng.randint(0, n // 2, (2, 1, n))
-    vc1 = vc2[[0]]
+    vc2 = np.ones((2, t, 1)) @ _integers(rng, n // 2, (2, 1, n))
 
     if const:
         x[0] = 1.0
@@ -760,14 +770,9 @@ def generate_panel_data(
         c, items=cats, major_axis=time, minor_axis=entities, swap=True
     )
     other_eff = c_df.reindex(index)
-    vc1_df = panel_to_frame(
-        vc1, items=vcats[:1], major_axis=time, minor_axis=entities, swap=True
-    )
-    vc1_df = vc1_df.reindex(index)
-    vc2_df = panel_to_frame(
+    clusters = panel_to_frame(
         vc2, items=vcats, major_axis=time, minor_axis=entities, swap=True
     )
-    vc2_df = vc2_df.reindex(index)
-    clusters = concat([vc1_df, vc2_df], sort=False)
+    clusters = clusters.reindex(index)
     data = concat([y_df, x_df], axis=1, sort=False)
     return PanelModelData(data, w_df, other_eff, clusters)
diff --git a/linearmodels/tests/panel/test_between_means.py b/linearmodels/tests/panel/test_between_means.py
new file mode 100644
index 0000000000..de58f5a23c
--- /dev/null
+++ b/linearmodels/tests/panel/test_between_means.py
@@ -0,0 +1,60 @@
+"""
+BetweenOLS where the means of the regressors are collinear, and df_model
+
+The regression of BetweenOLS uses the means of the entities. The regressors can be
+collinear in the means when they are not collinear in the data: there are fewer
+entities than regressors, or the dummy variables of the periods are the constant 1/T in
+the means of a balanced panel and so they are collinear with the constant. The
+covariance raised numpy's "Singular matrix" error when it was used (and the parameters
+were not identified), and the error has to say what is wrong.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from linearmodels.datasets import wage_panel
+from linearmodels.panel import BetweenOLS
+
+
+def make_panel(nentity, nperiod, nvar, seed=0):
+    rg = np.random.default_rng(seed)
+    index = pd.MultiIndex.from_product(
+        [np.arange(nentity), np.arange(nperiod)], names=["entity", "time"]
+    )
+    nobs = nentity * nperiod
+    x = pd.DataFrame(
+        rg.standard_normal((nobs, nvar)),
+        index=index,
+        columns=[f"x{i}" for i in range(nvar)],
+    )
+    y = pd.Series(rg.standard_normal(nobs), index=index, name="y")
+    return y, x
+
+
+def test_between_with_more_regressors_than_entities():
+    y, x = make_panel(2, 3, 3)
+    with pytest.raises(ValueError, match="collinear in the means of the entities"):
+        BetweenOLS(y, x).fit()
+    # The rank is not checked if the user does not want it
+    res = BetweenOLS(y, x, check_rank=False).fit(debiased=False)
+    assert res.params.shape == (3,)
+
+
+def test_between_with_time_dummies_in_a_balanced_panel():
+    data = wage_panel.load()
+    year = pd.Categorical(data.year)
+    data = data.set_index(["nr", "year"])
+    data["year"] = year
+    data["const"] = 1.0
+    with pytest.raises(ValueError, match="collinear in the means of the entities"):
+        BetweenOLS(data.lwage, data[["const", "expersq", "union", "year"]]).fit()
+    res = BetweenOLS(data.lwage, data[["const", "expersq", "union"]]).fit()
+    assert res.df_resid == 545 - 3
+
+
+def test_between_df_model_is_the_number_of_regressors():
+    y, x = make_panel(10, 3, 2)
+    res = BetweenOLS(y, x).fit()
+    assert res.df_model == 2
+    assert isinstance(res.df_model, int)
diff --git a/linearmodels/tests/panel/test_both_effects_degenerate.py b/linearmodels/tests/panel/test_both_effects_degenerate.py
new file mode 100644
index 0000000000..3785765494
--- /dev/null
+++ b/linearmodels/tests/panel/test_both_effects_degenerate.py
@@ -0,0 +1,49 @@
+"""
+Entity and time effects in a panel with a single entity or a single period
+
+The effects remove all of the variation of the regressors. The error has to say that the
+effects absorbed the variables. The default algorithm that removes both effects
+raised "No objects to concatenate" from pandas, since it uses the dummy variables of the
+shorter dimension, and there are none.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from linearmodels.panel import PanelOLS
+from linearmodels.panel.data import PanelData
+from linearmodels.panel.utility import AbsorbingEffectError
+
+
+def make_panel(nentity, nperiod, nvar, seed=0):
+    rg = np.random.default_rng(seed)
+    index = pd.MultiIndex.from_product(
+        [np.arange(nentity), np.arange(nperiod)], names=["entity", "time"]
+    )
+    nobs = nentity * nperiod
+    x = pd.DataFrame(
+        rg.standard_normal((nobs, nvar)),
+        index=index,
+        columns=[f"x{i}" for i in range(nvar)],
+    )
+    y = pd.Series(rg.standard_normal(nobs), index=index, name="y")
+    return y, x
+
+
+@pytest.mark.parametrize("low_memory", [None, False, True])
+@pytest.mark.parametrize("shape", [(1, 6), (6, 1)])
+def test_entity_and_time_effects_absorb_a_single_entity_or_period(shape, low_memory):
+    y, x = make_panel(*shape, 2)
+    mod = PanelOLS(y, x, entity_effects=True, time_effects=True)
+    with pytest.raises(AbsorbingEffectError, match="absorbed"):
+        mod.fit(low_memory=low_memory)
+
+
+@pytest.mark.parametrize("shape", [(1, 6), (6, 1)])
+def test_demean_both_of_a_single_entity_or_period(shape):
+    # There is nothing to remove besides the means of the group with the observations
+    _, x = make_panel(*shape, 2)
+    demeaned = PanelData(x).demean("both")
+    assert demeaned.dataframe.shape == x.shape
+    np.testing.assert_allclose(demeaned.values2d, 0.0, atol=1e-14)
diff --git a/linearmodels/tests/panel/test_check_rank_false.py b/linearmodels/tests/panel/test_check_rank_false.py
new file mode 100644
index 0000000000..4d35c85886
--- /dev/null
+++ b/linearmodels/tests/panel/test_check_rank_false.py
@@ -0,0 +1,124 @@
+"""
+The constant is detected whether or not the rank of the regressors is checked
+
+check_rank=False skips the computation of the rank of the regressors, and the models
+then use the number of regressors, which is the rank if they are not collinear. They
+used the number of periods of the panel instead. The detection of a constant that is
+not a column, e.g., dummy variables that add to one, compares the rank of the
+regressors with the rank of the regressors and a constant, so that it decided that
+there was no constant in a panel with 5 periods and 2 dummy variables, and that there
+was a constant in a panel with 2 periods and 3 regressors. The R-squared and the
+F-statistic then treated the models as if they had, or did not have, a constant.
+"""
+
+import numpy as np
+from numpy.testing import assert_allclose
+import pandas as pd
+import pytest
+
+from linearmodels.panel import (
+    BetweenOLS,
+    FamaMacBeth,
+    FirstDifferenceOLS,
+    PanelOLS,
+    PooledOLS,
+    RandomEffects,
+)
+
+MODELS = ["pooled", "between", "random-effects", "fama-macbeth", "panel-ols"]
+
+# Two cross-sections cannot give the standard errors of 3 parameters
+pytestmark = pytest.mark.filterwarnings(
+    "ignore:The number of cross-sectional regressions:"
+    "linearmodels.shared.exceptions.InferenceUnavailableWarning"
+)
+
+
+def build(model, y, x, check_rank):
+    if model == "panel-ols":
+        return PanelOLS(y, x, entity_effects=True, check_rank=check_rank)
+    models = {
+        "pooled": PooledOLS,
+        "between": BetweenOLS,
+        "random-effects": RandomEffects,
+        "fama-macbeth": FamaMacBeth,
+        "first-difference": FirstDifferenceOLS,
+    }
+    return models[model](y, x, check_rank=check_rank)
+
+
+def make_index(nentity, nperiod):
+    return pd.MultiIndex.from_product(
+        [np.arange(nentity), np.arange(nperiod)], names=["entity", "time"]
+    )
+
+
+def implicit_constant_panel():
+    """Two dummy variables that add to one, and 5 periods"""
+    rg = np.random.default_rng(0)
+    index = make_index(40, 5)
+    group = rg.integers(0, 2, len(index))
+    x = pd.DataFrame(
+        {"a": (group == 0) * 1.0, "b": (group == 1) * 1.0},
+        index=index,
+    )
+    y = pd.Series(rg.standard_normal(len(index)) + group, index=index, name="y")
+    return y, x
+
+
+def no_constant_panel():
+    """Three regressors with a mean of zero, and 2 periods"""
+    rg = np.random.default_rng(1)
+    index = make_index(40, 2)
+    x = pd.DataFrame(
+        rg.standard_normal((len(index), 3)), index=index, columns=["a", "b", "c"]
+    )
+    y = pd.Series(rg.standard_normal(len(index)) + 1, index=index, name="y")
+    return y, x
+
+
+def r_squared(y, x, centered):
+    yv, xv = y.to_numpy(), x.to_numpy()
+    resid = yv - xv @ np.linalg.lstsq(xv, yv, rcond=None)[0]
+    total = yv - yv.mean() if centered else yv
+    return 1 - resid @ resid / (total @ total)
+
+
+@pytest.mark.parametrize("model", MODELS)
+def test_implicit_constant(model):
+    y, x = implicit_constant_panel()
+    checked = build(model, y, x, True)
+    skipped = build(model, y, x, False)
+    assert checked.has_constant
+    assert skipped.has_constant
+    res_checked, res_skipped = checked.fit(), skipped.fit()
+    assert_allclose(res_skipped.params, res_checked.params)
+    assert_allclose(res_skipped.rsquared, res_checked.rsquared)
+    assert_allclose(res_skipped.f_statistic.stat, res_checked.f_statistic.stat)
+
+
+@pytest.mark.parametrize("model", MODELS + ["first-difference"])
+def test_no_constant_with_more_regressors_than_periods(model):
+    y, x = no_constant_panel()
+    checked = build(model, y, x, True)
+    skipped = build(model, y, x, False)
+    assert not checked.has_constant
+    assert not skipped.has_constant
+    res_checked, res_skipped = checked.fit(), skipped.fit()
+    assert_allclose(res_skipped.params, res_checked.params)
+    assert_allclose(res_skipped.rsquared, res_checked.rsquared)
+    assert_allclose(res_skipped.f_statistic.stat, res_checked.f_statistic.stat)
+
+
+@pytest.mark.parametrize("check_rank", [True, False])
+def test_pooled_r_squared_of_an_implicit_constant_is_centered(check_rank):
+    y, x = implicit_constant_panel()
+    res = PooledOLS(y, x, check_rank=check_rank).fit()
+    assert_allclose(res.rsquared, r_squared(y, x, True))
+
+
+@pytest.mark.parametrize("check_rank", [True, False])
+def test_pooled_r_squared_without_a_constant_is_not_centered(check_rank):
+    y, x = no_constant_panel()
+    res = PooledOLS(y, x, check_rank=check_rank).fit()
+    assert_allclose(res.rsquared, r_squared(y, x, False))
diff --git a/linearmodels/tests/panel/test_fama_macbeth.py b/linearmodels/tests/panel/test_fama_macbeth.py
index 6da1ab917c..16ba6342f7 100644
--- a/linearmodels/tests/panel/test_fama_macbeth.py
+++ b/linearmodels/tests/panel/test_fama_macbeth.py
@@ -129,9 +129,14 @@ def test_block_size_warnings():
     )
     y = pd.DataFrame(y, index=idx, columns=["y"])
     x = pd.DataFrame(x, index=idx, columns=["x1", "x2", "x3"])
-    with pytest.warns(MissingValueWarning):
+    # The periods are cross-sectional regressions, and not time-series regressions
+    with pytest.warns(
+        MissingValueWarning, match=r"1 of the cross-sectional regressions"
+    ):
         FamaMacBeth(y.iloc[:11], x.iloc[:11])
-    with pytest.warns(InferenceUnavailableWarning):
+    with pytest.warns(
+        InferenceUnavailableWarning, match=r"cross-sectional regressions that can be"
+    ):
         FamaMacBeth(y.iloc[::4], x.iloc[::4])
 
 
@@ -338,7 +343,9 @@ def test_inference_n40_t10():
 def test_inference_dropped_periods(bandwidth, debiased):
     # T is the number of periods that are estimated, 21 of the 25 periods
     y, x = dropped_periods_panel()
-    with pytest.warns(MissingValueWarning, match=r"4 of the time-series regressions"):
+    with pytest.warns(
+        MissingValueWarning, match=r"4 of the cross-sectional regressions"
+    ):
         mod = FamaMacBeth(y, x)
     cov_type = "kernel" if bandwidth else "unadjusted"
     res = mod.fit(
diff --git a/linearmodels/tests/panel/test_formula.py b/linearmodels/tests/panel/test_formula.py
index 7897ae6ec7..6573ab5636 100644
--- a/linearmodels/tests/panel/test_formula.py
+++ b/linearmodels/tests/panel/test_formula.py
@@ -334,6 +334,16 @@ def _mr(z):
         # FMB has tighter checks that cannot be disabled
         return
 
+    if model == FirstDifferenceOLS:
+        # The check for a constant reads regressors that are not full rank as having an
+        # implicit constant, which first differences do not allow. The rank that it
+        # used when the rank was not checked was the number of periods, not the rank
+        with pytest.raises(ValueError, match=r"Constants are not allowed"):
+            model.from_formula(formula, joined, check_rank=False)
+        with pytest.raises(ValueError, match=r"Constants are not allowed"):
+            model(y, x, check_rank=False)
+        return
+
     mod = model.from_formula(formula, joined, check_rank=False)
     res = mod.fit()
     assert isinstance(res.summary.as_text(), str)
diff --git a/linearmodels/tests/panel/test_no_residual_df.py b/linearmodels/tests/panel/test_no_residual_df.py
new file mode 100644
index 0000000000..abe21026ba
--- /dev/null
+++ b/linearmodels/tests/panel/test_no_residual_df.py
@@ -0,0 +1,151 @@
+"""
+Models without residual degrees of freedom
+
+A model that has as many parameters as observations, after the effects are removed,
+raised the ZeroDivisionError of the scale of the covariance, n / (n - k), or of an
+F-statistic. The covariance of such a model is not defined, and the error has to say
+so. The F-statistic of the model is not defined either, which is the result if the
+covariance is not debiased. The F-statistic that tests that the effects are zero is not
+defined if the effects do not use any degree of freedom: a single entity with a
+constant, or a single period with time effects.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from linearmodels.panel import (
+    BetweenOLS,
+    FamaMacBeth,
+    FirstDifferenceOLS,
+    PanelOLS,
+    PooledOLS,
+)
+from linearmodels.shared.exceptions import InferenceUnavailableWarning
+from linearmodels.shared.hypotheses import (
+    InapplicableTestStatistic,
+    InvalidTestStatistic,
+)
+
+NO_DF = "no residual degrees of freedom"
+
+
+def make_panel(nentity, nperiod, nvar, seed=0):
+    rg = np.random.default_rng(seed)
+    index = pd.MultiIndex.from_product(
+        [np.arange(nentity), np.arange(nperiod)], names=["entity", "time"]
+    )
+    nobs = nentity * nperiod
+    x = pd.DataFrame(
+        rg.standard_normal((nobs, nvar)),
+        index=index,
+        columns=[f"x{i}" for i in range(nvar)],
+    )
+    y = pd.Series(rg.standard_normal(nobs), index=index, name="y")
+    return y, x
+
+
+def saturated_pooled():
+    # 2 observations, a constant and a regressor
+    y, x = make_panel(1, 2, 1)
+    x.insert(0, "const", 1.0)
+    return y, x
+
+
+def saturated_between():
+    # 2 entities, a constant and a regressor in the regression of the means
+    y, x = make_panel(2, 3, 1)
+    x.insert(0, "const", 1.0)
+    return y, x
+
+
+# Each model has as many parameters as observations once the effects are removed
+SATURATED = {
+    "pooled": lambda: PooledOLS(*saturated_pooled()),
+    # 2 entities with 2 periods, 2 regressors and 2 effects
+    "entity-effects": lambda: PanelOLS(*make_panel(2, 2, 2), entity_effects=True),
+    "time-effects": lambda: PanelOLS(*make_panel(2, 2, 2), time_effects=True),
+    "between": lambda: BetweenOLS(*saturated_between()),
+    # 1 entity with 3 periods has 2 differences, and there are 2 regressors
+    "first-difference": lambda: FirstDifferenceOLS(*make_panel(1, 3, 2)),
+}
+# BetweenOLS does not have the kernel covariances
+COVARIANCES = [
+    (model, cov_type)
+    for model in SATURATED
+    for cov_type in (["robust"] if model == "between" else ["robust", "kernel", "ac"])
+]
+
+
+@pytest.mark.parametrize("model", list(SATURATED))
+def test_saturated_model_raises_value_error(model):
+    mod = SATURATED[model]()
+    with pytest.raises(ValueError, match=NO_DF) as exc_info:
+        mod.fit()
+    assert exc_info.type is ValueError
+
+
+@pytest.mark.parametrize(("model", "cov_type"), COVARIANCES)
+def test_saturated_model_with_other_covariances(model, cov_type):
+    with pytest.raises(ValueError, match=NO_DF):
+        SATURATED[model]().fit(cov_type=cov_type)
+
+
+@pytest.mark.parametrize("model", list(SATURATED))
+def test_saturated_model_that_is_not_debiased(model):
+    # The estimates and the covariance do not need the degrees of freedom. The
+    # F-statistic does, and it is not defined
+    res = SATURATED[model]().fit(debiased=False)
+    assert np.all(np.isfinite(res.params))
+    assert isinstance(res.f_statistic, InvalidTestStatistic)
+    assert np.isnan(res.f_statistic.stat)
+    assert NO_DF in str(res.f_statistic)
+    assert str(res.summary)
+
+
+def test_famamacbeth_with_as_many_observations_as_parameters():
+    # One cross-section that has as many entities as regressors
+    y, x = make_panel(2, 1, 2)
+    with pytest.warns(InferenceUnavailableWarning):
+        mod = FamaMacBeth(y, x)
+    with pytest.raises(ValueError, match=NO_DF):
+        mod.fit()
+
+
+def single_entity(const):
+    y, x = make_panel(1, 8, 2)
+    if const:
+        x.insert(0, "const", 1.0)
+    return y, x
+
+
+def single_period(const):
+    y, x = make_panel(8, 1, 2)
+    if const:
+        x.insert(0, "const", 1.0)
+    return y, x
+
+
+@pytest.mark.parametrize("const", [False, True])
+def test_single_entity_with_entity_effects(const):
+    # The effect is the constant: the slopes are those of the regression with a constant
+    y, x = single_entity(const)
+    res = PanelOLS(y, x, entity_effects=True).fit()
+    slopes = ["x0", "x1"]
+    pooled_x = x[slopes].copy()
+    pooled_x.insert(0, "const", 1.0)
+    pooled = PooledOLS(y, pooled_x).fit()
+    np.testing.assert_allclose(res.params[slopes], pooled.params[slopes])
+    # The effects do not use any degree of freedom: there is nothing to test
+    assert isinstance(res.f_pooled, InapplicableTestStatistic)
+    assert "effects" in str(res.f_pooled)
+    assert np.isnan(res.f_pooled.stat)
+    assert str(res.summary)
+
+
+@pytest.mark.parametrize("const", [False, True])
+def test_single_period_with_time_effects(const):
+    y, x = single_period(const)
+    res = PanelOLS(y, x, time_effects=True).fit()
+    assert isinstance(res.f_pooled, InapplicableTestStatistic)
+    assert np.all(np.isfinite(res.params))
diff --git a/linearmodels/tests/panel/test_utility.py b/linearmodels/tests/panel/test_utility.py
index 8613a473e8..723fcfa127 100644
--- a/linearmodels/tests/panel/test_utility.py
+++ b/linearmodels/tests/panel/test_utility.py
@@ -5,6 +5,7 @@
 import pandas as pd
 import pytest
 
+from linearmodels.panel import PanelOLS
 from linearmodels.panel.utility import (
     AbsorbingEffectError,
     PanelModelData,
@@ -224,7 +225,9 @@ class SubArray(np.ndarray):
 @pytest.mark.parametrize("const", [True, False])
 @pytest.mark.parametrize("other_effects", [0, 1, 2])
 @pytest.mark.parametrize("cat_list", [True, False])
-@pytest.mark.parametrize("random_state", [None, np.random.RandomState(0)])
+@pytest.mark.parametrize(
+    "random_state", [None, np.random.RandomState(0), np.random.default_rng(0)]
+)
 def test_generate_panel_data(missing, const, other_effects, cat_list, random_state):
     if cat_list:
         ncats: list[int] | int = [13] * other_effects
@@ -232,7 +235,11 @@ def test_generate_panel_data(missing, const, other_effects, cat_list, random_sta
         ncats = 21
 
     dataset = generate_panel_data(
-        missing=missing, const=const, other_effects=other_effects, ncats=ncats
+        missing=missing,
+        const=const,
+        other_effects=other_effects,
+        ncats=ncats,
+        rng=random_state,
     )
 
     assert isinstance(dataset, PanelModelData)
@@ -244,6 +251,81 @@ def test_generate_panel_data(missing, const, other_effects, cat_list, random_sta
     assert dataset.other_effects.shape == (dataset.data.shape[0], other_effects)
 
 
+@pytest.mark.parametrize("make_rng", [np.random.RandomState, np.random.default_rng])
+def test_generate_panel_data_rng(make_rng):
+    # A Generator has integers and not randint. The same seed gives the same data
+    options = {"nentity": 30, "ntime": 5, "other_effects": 2, "ncats": [3, 4]}
+    first = generate_panel_data(rng=make_rng(1234), **options)
+    second = generate_panel_data(rng=make_rng(1234), **options)
+    other = generate_panel_data(rng=make_rng(4321), **options)
+    for result in (first, second, other):
+        assert result.data.shape == (150, 6)
+        assert result.other_effects.shape == (150, 2)
+    pd.testing.assert_frame_equal(first.data, second.data)
+    pd.testing.assert_frame_equal(first.other_effects, second.other_effects)
+    pd.testing.assert_frame_equal(first.clusters, second.clusters)
+    assert not first.data.equals(other.data)
+    # The categories of the effects are integers in the ranges
+    effects = first.other_effects
+    assert set(effects["cat.0"].unique()) <= set(range(3))
+    assert set(effects["cat.1"].unique()) <= set(range(4))
+
+
+def test_generate_panel_data_random_state_is_unchanged():
+    # The draws of a RandomState do not depend on Generators being accepted
+    result = generate_panel_data(
+        nentity=6,
+        ntime=3,
+        nexog=2,
+        other_effects=2,
+        ncats=[3, 5],
+        rng=np.random.RandomState(7),
+    )
+    effects = result.other_effects.to_numpy().T
+    assert_array_equal(
+        effects,
+        [
+            [1, 0, 1, 2, 1, 0, 1, 1, 0, 2, 0, 2, 2, 1, 1, 0, 1, 1],
+            [2, 1, 1, 0, 3, 1, 2, 1, 3, 0, 0, 2, 2, 0, 0, 1, 1, 2],
+        ],
+    )
+
+
+def test_generate_panel_data_clusters_random_state_is_unchanged():
+    # The clusters are the draws that the second column of the clusters had
+    result = generate_panel_data(
+        nentity=6,
+        ntime=3,
+        nexog=2,
+        other_effects=2,
+        ncats=[3, 5],
+        rng=np.random.RandomState(7),
+    )
+    first_periods = result.clusters.iloc[::3]
+    assert_array_equal(first_periods["varcat.0"], [2, 2, 0, 2, 0, 0])
+    assert_array_equal(first_periods["varcat.1"], [1, 0, 1, 0, 2, 2])
+
+
+@pytest.mark.parametrize("missing", [0, 0.1])
+def test_generate_panel_data_clusters(missing):
+    # The clusters have one row for each observation, and two variables
+    result = generate_panel_data(nentity=20, ntime=5, missing=missing)
+    assert result.clusters.shape == (result.data.shape[0], 2)
+    assert result.clusters.index.equals(result.data.index)
+    assert not result.clusters.isnull().any().any()
+    assert result.clusters.to_numpy().max() < 10
+    assert list(result.clusters.columns) == ["varcat.0", "varcat.1"]
+    # They are the clusters of the entities: the same in all periods of an entity
+    assert (result.clusters.groupby(level=0).nunique() == 1).all().all()
+    # And can be used in a covariance
+    data = result.data.dropna()
+    clusters = result.clusters.loc[data.index]
+    res = PanelOLS(data["y"], data[["x0", "x1"]], entity_effects=True).fit(
+        cov_type="clustered", clusters=clusters
+    )
+    assert np.all(np.isfinite(res.params))
+
+
 def test_not_absorbed_const():
     x = np.random.standard_normal((200, 3))
     x[:, 0] = 0
diff --git a/linearmodels/tests/panel/test_weights_validation.py b/linearmodels/tests/panel/test_weights_validation.py
new file mode 100644
index 0000000000..e3e4c74eb4
--- /dev/null
+++ b/linearmodels/tests/panel/test_weights_validation.py
@@ -0,0 +1,80 @@
+"""
+The weights of the panel models must be finite and strictly positive
+
+The weights are scaled to have a mean of one. Zero and negative weights were checked, but
+an infinite weight was scaled to NaN and the other weights to zero, and the models failed
+with errors from NumPy, SciPy and pandas that do not say what is wrong ("array must not
+contain infs or NaNs", "SVD did not converge", an AbsorbingEffectError, and "operands could not
+be broadcast together" for FirstDifferenceOLS). Missing weights are dropped with their
+observations, as the other missing values are, and have to give the model of the data
+without them.
+"""
+
+import numpy as np
+from numpy.testing import assert_allclose
+import pandas as pd
+import pytest
+
+from linearmodels.panel import (
+    BetweenOLS,
+    FamaMacBeth,
+    FirstDifferenceOLS,
+    PanelOLS,
+    PooledOLS,
+    RandomEffects,
+)
+from linearmodels.shared.exceptions import MissingValueWarning
+
+MODELS = {
+    "pooled": lambda y, x, w: PooledOLS(y, x, weights=w),
+    "panel-entity-effects": lambda y, x, w: PanelOLS(
+        y, x, weights=w, entity_effects=True
+    ),
+    "panel-time-effects": lambda y, x, w: PanelOLS(y, x, weights=w, time_effects=True),
+    "between": lambda y, x, w: BetweenOLS(y, x, weights=w),
+    "first-difference": lambda y, x, w: FirstDifferenceOLS(y, x, weights=w),
+    "random-effects": lambda y, x, w: RandomEffects(y, x, weights=w),
+    "fama-macbeth": lambda y, x, w: FamaMacBeth(y, x, weights=w),
+}
+
+
+def make_panel(seed=0):
+    rg = np.random.default_rng(seed)
+    index = pd.MultiIndex.from_product(
+        [np.arange(30), np.arange(6)], names=["entity", "time"]
+    )
+    x = pd.DataFrame(
+        rg.standard_normal((index.shape[0], 2)), index=index, columns=["a", "b"]
+    )
+    y = pd.Series(
+        x @ [0.5, -0.3] + rg.standard_normal(index.shape[0]), index=index, name="y"
+    )
+    w = pd.Series(rg.chisquare(4, index.shape[0]) / 4 + 0.2, index=index, name="w")
+    return y, x, w
+
+
+@pytest.mark.parametrize("bad", [np.inf, -np.inf, 0.0, -1.0])
+@pytest.mark.parametrize("model", list(MODELS))
+def test_weights_must_be_finite_and_strictly_positive(model, bad):
+    y, x, w = make_panel()
+    w.iloc[7] = bad
+    with pytest.raises(
+        ValueError, match="weights must be strictly positive and finite"
+    ):
+        MODELS[model](y, x, w)
+
+
+@pytest.mark.parametrize("model", list(MODELS))
+def test_missing_weights_are_dropped_with_their_observations(model):
+    y, x, w = make_panel()
+    missing = [3, 17, 18, 100]
+    w_missing = w.copy()
+    w_missing.iloc[missing] = np.nan
+    with pytest.warns(MissingValueWarning):
+        res = MODELS[model](y, x, w_missing).fit()
+    keep = np.ones(y.shape[0], dtype=bool)
+    keep[missing] = False
+    expected = MODELS[model](y[keep], x[keep], w[keep]).fit()
+    assert res.nobs == expected.nobs
+    assert_allclose(res.params, expected.params)
+    assert_allclose(res.std_errors, expected.std_errors)