BUG: Raise a ValueError for a debiased factor model without residual …

linearmodels/asset_pricing/model.py

@@ -136,9 +136,20 @@ def _validate_data(self) -> None:
                 "equivalent and must not have rank\n"
                 "less than the dimension of the smaller shape."
             )
-        if has_constant(f)[0]:
+        rank = np.linalg.matrix_rank(f)
+        constant = has_constant(f, rank)[0]
+        if constant and rank < f.shape[1]:
+            # has_constant is also True if the columns are collinear. Collinear
+            # factors have a constant if a column is constant, or, for a matrix
+            # with more rows than columns, if a constant is in the span of the
+            # columns
+            nobs = f.shape[0]
+            aug_rank = np.linalg.matrix_rank(np.c_[np.ones((nobs, 1)), f])
+            constant_column = np.any((np.ptp(f, axis=0) == 0) & ~np.all(f == 0, axis=0))
+            constant = bool(constant_column or (nobs > f.shape[1] and aug_rank == rank))
+        if constant:
             raise ValueError("factors must not contain a constant or equivalent.")
-        if np.linalg.matrix_rank(f) < f.shape[1]:
+        if rank < f.shape[1]:
             raise ValueError(
                 "Model cannot be estimated. factors do not have full column rank."
             )
@@ -148,6 +159,18 @@ def _validate_data(self) -> None:
                 "more returns to\nestimate the model parameters."
             )
 
+    def _check_debiased(self, debiased: bool) -> None:
+        """Debiasing needs more observations than regressors in the regressions"""
+        nobs, nfactor = self.factors.shape
+        if debiased and nobs <= nfactor + 1:
+            raise ValueError(
+                "The covariance cannot be debiased when the number of observations "
+                f"({nobs}) is not larger than the number of factors plus one "
+                f"({nfactor + 1}), the number of regressors of each time-series "
+                "regression, since there are no residual degrees of freedom. Use "
+                "debiased=False or more observations."
+            )
+
     @property
     def formula(self) -> str | None:
         return self._formula
@@ -326,6 +349,7 @@ def fit(
         ``kernel``, one of "bartlett", "parzen" or "qs" (quadratic spectral)
         and ``bandwidth`` (a positive integer).
         """
+        self._check_debiased(debiased)
         p = self.portfolios.ndarray
         f = self.factors.ndarray
         nportfolio = p.shape[1]
@@ -459,11 +483,33 @@ def __init__(
                 self.portfolios.shape[1]
             )
         else:
-            self._sigma = np.asarray(sigma)
-            vals, vecs = np.linalg.eigh(sigma)
+            self._sigma = self._validate_sigma(sigma)
+            vals, vecs = np.linalg.eigh(self._sigma)
             self._sigma_m12 = vecs @ np.diag(1.0 / np.sqrt(vals)) @ vecs.T
             self._sigma_inv = np.linalg.inv(self._sigma).astype(float, copy=False)
 
+    def _validate_sigma(
+        self, sigma: linearmodels.typing.data.ArrayLike
+    ) -> linearmodels.typing.data.Float64Array:
+        """Check that sigma is a symmetric positive definite matrix of the right size"""
+        nportfolio = self.portfolios.shape[1]
+        sigma_arr = np.asarray(sigma, dtype=float)
+        if sigma_arr.shape != (nportfolio, nportfolio):
+            raise ValueError(
+                f"sigma must be a square array with shape ({nportfolio}, "
+                f"{nportfolio}), where {nportfolio} is the number of test portfolios. "
+                f"The shape of sigma is {sigma_arr.shape}."
+            )
+        if not np.all(np.isfinite(sigma_arr)):
+            raise ValueError("sigma must not contain missing or infinite values.")
+        scale = np.max(np.abs(sigma_arr))
+        if np.max(np.abs(sigma_arr - sigma_arr.T)) > 1e-8 * scale:
+            raise ValueError("sigma must be symmetric.")
+        vals = np.linalg.eigvalsh((sigma_arr + sigma_arr.T) / 2)
+        if not vals[0] > nportfolio * np.finfo(float).eps * vals[-1]:
+            raise ValueError("sigma must be positive definite.")
+        return sigma_arr
+
     def __str__(self) -> str:
         out = super().__str__()
         if np.any(self._sigma != np.eye(self.portfolios.shape[1])):
@@ -509,8 +555,12 @@ class LinearFactorModel(_LinearFactorModelBase):
         Flag indicating whether the risk-free rate should be estimated
         from returns along other risk premia.  If False, the returns are
         assumed to be excess returns using the correct risk-free rate.
-    sigma : array_like
-        Positive definite residual covariance (nportfolio by nportfolio)
+    sigma : array_like, optional
+        Positive definite residual covariance (nportfolio by nportfolio) that
+        weights the test portfolios in the second stage. A ``ValueError`` is
+        raised if it does not have this shape, contains missing or infinite
+        values, or is not symmetric and positive definite. If not provided, the
+        identity matrix is used.
 
     Notes
     -----
@@ -591,8 +641,12 @@ def from_formula(
             Flag indicating whether the risk-free rate should be estimated
             from returns along other risk premia.  If False, the returns are
             assumed to be excess returns using the correct risk-free rate.
-        sigma : array_like
-            Positive definite residual covariance (nportfolio by nportfolio)
+        sigma : array_like, optional
+            Positive definite residual covariance (nportfolio by nportfolio) that
+            weights the test portfolios in the second stage. A ``ValueError`` is
+            raised if it does not have this shape, contains missing or infinite
+            values, or is not symmetric and positive definite. If not provided,
+            the identity matrix is used.
 
         Returns
         -------
@@ -659,6 +713,7 @@ def fit(
         """
         nobs, nf, nport, nrf, s1, s2, s3 = self._boundaries()
         excess_returns = not self._risk_free
+        self._check_debiased(debiased)
         f = self.factors.ndarray
         p = self.portfolios.ndarray
         nport = p.shape[1]
@@ -1009,7 +1064,9 @@ def fit(
         mod = LinearFactorModel(
             self.portfolios, self.factors, risk_free=self._risk_free
         )
-        res = mod.fit()
+        # The starting values do not use the covariance, which cannot be debiased
+        # if there are no residual degrees of freedom in the first stage
+        res = mod.fit(debiased=False)
         betas_1d = np.asarray(res.betas).ravel()
         lam = np.asarray(res.risk_premia)
         mu = self.factors.ndarray.mean(0)

linearmodels/shared/io.py

@@ -10,6 +10,8 @@ def _str(v: float) -> str:
     """Preferred basic formatter"""
     if np.isnan(v):
         return "        "
+    if np.isinf(v):
+        return f"{v:8.4g}"
     av = abs(v)
     digits = 0
     if av != 0:

linearmodels/tests/asset_pricing/test_factor_collinearity.py

@@ -0,0 +1,96 @@
+"""
+Tests of the messages for factors that are collinear or contain a constant
+
+has_constant is also True for a matrix that does not have full column rank, so that
+collinear factors were reported as containing a constant and the message for factors
+without full column rank could not be reached.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from linearmodels.asset_pricing.model import (
+    LinearFactorModel,
+    LinearFactorModelGMM,
+    TradedFactorModel,
+)
+
+MODELS = [TradedFactorModel, LinearFactorModel, LinearFactorModelGMM]
+RANK = r"Model cannot be estimated\. factors do not have full column rank\."
+CONSTANT = r"factors must not contain a constant or equivalent\."
+
+
+@pytest.fixture(scope="module")
+def data():
+    rng = np.random.default_rng(0)
+    factors = rng.standard_normal((100, 3))
+    portfolios = factors @ rng.standard_normal((3, 6)) + rng.standard_normal((100, 6))
+    return portfolios, factors
+
+
+def collinear_factors(kind, factors):
+    x1, x2, x3 = factors.T
+    ones = np.ones_like(x1)
+    cases = {
+        "duplicate": np.c_[factors, x1],
+        "combination": np.c_[factors, x1 - 2 * x2],
+        "combination first": np.c_[x1 + x2, factors],
+        # A constant column, or a constant in the span, makes the factors contain one
+        "duplicate constant": np.c_[factors, ones, ones],
+        "constant in span": np.c_[x1, x2, 1 - x1 - x2],
+        "constant in span, collinear": np.c_[x1, x2, 1 - x1 - x2, x1 + x2],
+        "constant and duplicate": np.c_[factors, x1, ones],
+        "constant not one": np.c_[factors, 3 * ones],
+        "constant not one and duplicate": np.c_[factors, 3 * ones, x1],
+        "scaled constant in span": np.c_[x1, x2, 5 - 2 * x1 - 2 * x2],
+        "x3 only": x3[:, None],
+    }
+    return cases[kind]
+
+
+@pytest.mark.parametrize("model", MODELS)
+@pytest.mark.parametrize("kind", ["duplicate", "combination", "combination first"])
+@pytest.mark.parametrize("pandas", [False, True])
+def test_collinear_factors(data, model, kind, pandas):
+    portfolios, factors = data
+    collinear = collinear_factors(kind, factors)
+    if pandas:
+        collinear = pd.DataFrame(collinear)
+    with pytest.raises(ValueError, match=RANK):
+        model(portfolios, collinear)
+
+
+@pytest.mark.parametrize("model", MODELS)
+@pytest.mark.parametrize(
+    "kind",
+    [
+        "duplicate constant",
+        "constant in span",
+        "constant in span, collinear",
+        "constant and duplicate",
+        "constant not one",
+        "constant not one and duplicate",
+        "scaled constant in span",
+    ],
+)
+def test_constant_is_reported_as_a_constant(data, model, kind):
+    portfolios, factors = data
+    with pytest.raises(ValueError, match=CONSTANT):
+        model(portfolios, collinear_factors(kind, factors))
+
+
+@pytest.mark.parametrize("model", MODELS)
+def test_valid_factors(data, model):
+    portfolios, factors = data
+    mod = model(portfolios, factors)
+    assert mod.factors.shape == (100, 3)
+    mod = model(portfolios, collinear_factors("x3 only", factors))
+    assert mod.factors.shape == (100, 1)
+
+
+def test_fewer_observations_than_factors(data):
+    # A wide matrix does not have full column rank, whatever it contains
+    portfolios, factors = data
+    with pytest.raises(ValueError, match=RANK):
+        TradedFactorModel(portfolios[:2], factors[:2])

linearmodels/tests/asset_pricing/test_model.py

@@ -169,7 +169,7 @@ def test_errors(data):
         TradedFactorModel(p3, data.factors)
     with pytest.raises(ValueError, match=r"factors must not contain"):
         TradedFactorModel(data.portfolios, f)
-    with pytest.raises(ValueError, match=r"factors must not contain"):
+    with pytest.raises(ValueError, match=r"factors do not have full column rank"):
         TradedFactorModel(data.portfolios, f2)
     with pytest.raises(ValueError, match=r"Model cannot be estimated"):
         TradedFactorModel(p5, f5)

linearmodels/tests/asset_pricing/test_saturated_model.py

@@ -0,0 +1,69 @@
+"""
+Tests of models with as many observations as regressors in the time-series regressions
+
+With nobs equal to nfactor + 1 the regressions have no residual degrees of freedom,
+so the covariance cannot be debiased.
+"""
+
+import numpy as np
+import pytest
+
+from linearmodels.asset_pricing.model import (
+    LinearFactorModel,
+    LinearFactorModelGMM,
+    TradedFactorModel,
+)
+
+NFACTOR = 2
+NPORTFOLIO = 5
+
+
+def make_data(nobs):
+    rng = np.random.default_rng(0)
+    factors = rng.standard_normal((nobs, NFACTOR))
+    betas = rng.standard_normal((NFACTOR, NPORTFOLIO))
+    portfolios = factors @ betas + 0.1 * rng.standard_normal((nobs, NPORTFOLIO))
+    return portfolios, factors
+
+
+@pytest.mark.parametrize("model", [TradedFactorModel, LinearFactorModel])
+@pytest.mark.parametrize("fit_options", [{}, {"debiased": True}])
+def test_saturated_model_cannot_be_debiased(model, fit_options):
+    portfolios, factors = make_data(NFACTOR + 1)
+    mod = model(portfolios, factors)
+    match = r"cannot be debiased.*\(3\).*\(3\).*Use debiased=False"
+    with pytest.raises(ValueError, match=match):
+        mod.fit(**fit_options)
+    with pytest.raises(ValueError, match=match):
+        mod.fit(cov_type="kernel", bandwidth=1, **fit_options)
+
+
+@pytest.mark.parametrize("model", [TradedFactorModel, LinearFactorModel])
+def test_saturated_model_without_debiasing(model):
+    portfolios, factors = make_data(NFACTOR + 1)
+    res = model(portfolios, factors).fit(debiased=False)
+    assert res.nobs == NFACTOR + 1
+    assert np.all(np.isfinite(res.params))
+
+
+@pytest.mark.parametrize("model", [TradedFactorModel, LinearFactorModel])
+def test_one_more_observation_can_be_debiased(model):
+    portfolios, factors = make_data(NFACTOR + 2)
+    res = model(portfolios, factors).fit()
+    assert np.all(np.isfinite(res.cov))
+
+
+def test_gmm_starting_values_are_not_debiased(monkeypatch):
+    # The starting values come from LinearFactorModel, and do not use the covariance,
+    # so that they must not fail if the first stage has no degrees of freedom
+    portfolios, factors = make_data(40)
+    calls = []
+    fit = LinearFactorModel.fit
+
+    def record(self, *args, **kwargs):
+        calls.append(kwargs.get("debiased", True))
+        return fit(self, *args, **kwargs)
+
+    monkeypatch.setattr(LinearFactorModel, "fit", record)
+    LinearFactorModelGMM(portfolios, factors).fit(disp=0)
+    assert calls == [False]