diff --git a/linearmodels/asset_pricing/model.py b/linearmodels/asset_pricing/model.py index aaf0f5c296..fd99047975 100644 --- a/linearmodels/asset_pricing/model.py +++ b/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) diff --git a/linearmodels/shared/io.py b/linearmodels/shared/io.py index 8e296c837a..28a62b7df1 100644 --- a/linearmodels/shared/io.py +++ b/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: diff --git a/linearmodels/tests/asset_pricing/test_factor_collinearity.py b/linearmodels/tests/asset_pricing/test_factor_collinearity.py new file mode 100644 index 0000000000..18cda77b81 --- /dev/null +++ b/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]) diff --git a/linearmodels/tests/asset_pricing/test_model.py b/linearmodels/tests/asset_pricing/test_model.py index 18ed297794..de937a0cd8 100644 --- a/linearmodels/tests/asset_pricing/test_model.py +++ b/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) diff --git a/linearmodels/tests/asset_pricing/test_saturated_model.py b/linearmodels/tests/asset_pricing/test_saturated_model.py new file mode 100644 index 0000000000..6a6fe45026 --- /dev/null +++ b/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] diff --git a/linearmodels/tests/asset_pricing/test_sigma_validation.py b/linearmodels/tests/asset_pricing/test_sigma_validation.py new file mode 100644 index 0000000000..c5f9d5182b --- /dev/null +++ b/linearmodels/tests/asset_pricing/test_sigma_validation.py @@ -0,0 +1,132 @@ +""" +Tests of the validation of sigma, the covariance used for GLS by LinearFactorModel + +sigma must be a symmetric positive definite array with one row and one column for each +test portfolio. +""" + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.asset_pricing.model import LinearFactorModel + +NOBS = 200 +NPORTFOLIO = 6 + + +@pytest.fixture(scope="module") +def data(): + rng = np.random.default_rng(0) + factors = rng.standard_normal((NOBS, 2)) + betas = rng.standard_normal((2, NPORTFOLIO)) + portfolios = factors @ betas + rng.standard_normal((NOBS, NPORTFOLIO)) + return portfolios, factors + + +@pytest.mark.parametrize( + "shape", [(5, 5), (7, 7), (6, 5), (5, 6), (6, 6, 1), (36,), ()] +) +def test_wrong_shape(data, shape): + portfolios, factors = data + sigma = np.ones(shape) if shape else np.array(1.0) + with pytest.raises(ValueError, match=r"sigma must be a square array with shape"): + LinearFactorModel(portfolios, factors, sigma=sigma) + + +def test_wrong_shape_message(data): + portfolios, factors = data + match = r"shape \(6, 6\), where 6 is the number of test portfolios.*\(5, 5\)" + with pytest.raises(ValueError, match=match): + LinearFactorModel(portfolios, factors, sigma=np.eye(5)) + + +@pytest.mark.parametrize("bad", [np.nan, np.inf, -np.inf]) +def test_not_finite(data, bad): + portfolios, factors = data + sigma = np.eye(NPORTFOLIO) + sigma[1, 2] = sigma[2, 1] = bad + with pytest.raises(ValueError, match=r"sigma must not contain missing or infinite"): + LinearFactorModel(portfolios, factors, sigma=sigma) + + +def test_not_symmetric(data): + portfolios, factors = data + sigma = np.eye(NPORTFOLIO) + sigma[1, 2] = 0.2 + with pytest.raises(ValueError, match=r"sigma must be symmetric"): + LinearFactorModel(portfolios, factors, sigma=sigma) + + +@pytest.mark.parametrize("kind", ["negative", "indefinite", "singular", "zero"]) +def test_not_positive_definite(data, kind): + portfolios, factors = data + if kind == "negative": + sigma = -np.eye(NPORTFOLIO) + elif kind == "indefinite": + sigma = np.eye(NPORTFOLIO) + sigma[0, 0] = -1.0 + elif kind == "singular": + sigma = np.eye(NPORTFOLIO) + sigma[0, 1] = sigma[1, 0] = 1.0 + else: + sigma = np.zeros((NPORTFOLIO, NPORTFOLIO)) + with pytest.raises(ValueError, match=r"sigma must be positive definite"): + LinearFactorModel(portfolios, factors, sigma=sigma) + + +def test_sample_covariance_with_too_few_observations(): + # The covariance of 6 portfolios from 4 observations is singular + rng = np.random.default_rng(1) + portfolios = rng.standard_normal((4, NPORTFOLIO)) + factors = rng.standard_normal((4, 1)) + with pytest.raises(ValueError, match=r"sigma must be positive definite"): + LinearFactorModel(portfolios, factors, sigma=np.cov(portfolios.T)) + + +@pytest.mark.parametrize("kind", ["array", "frame", "list", "ints"]) +def test_valid_sigma(data, kind): + portfolios, factors = data + sigma = np.cov(portfolios.T) + if kind == "frame": + sigma = pd.DataFrame(sigma) + elif kind == "list": + sigma = sigma.tolist() + elif kind == "ints": + sigma = np.diag(np.arange(1, NPORTFOLIO + 1)) + mod = LinearFactorModel(portfolios, factors, sigma=sigma) + assert "using GLS" in str(mod) + res = mod.fit() + sigma_arr = np.asarray(sigma, dtype=float) + # The second stage is GLS with the sigma that was given + betas = np.asarray(res.betas) + lam = np.linalg.solve( + betas.T @ np.linalg.solve(sigma_arr, betas), + betas.T @ np.linalg.solve(sigma_arr, portfolios.mean(0)), + ) + assert_allclose(res.risk_premia, lam) + + +def test_symmetric_up_to_rounding(data): + # A product of matrices is symmetric up to rounding + portfolios, factors = data + rng = np.random.default_rng(2) + a = rng.standard_normal((NPORTFOLIO, NPORTFOLIO)) + sigma = a @ np.diag(np.arange(1.0, NPORTFOLIO + 1)) @ a.T + sigma = sigma + 1e-14 * rng.standard_normal((NPORTFOLIO, NPORTFOLIO)) + res = LinearFactorModel(portfolios, factors, sigma=sigma).fit() + assert np.all(np.isfinite(res.cov)) + + +def test_from_formula_validates_sigma(data): + portfolios, factors = data + frame = pd.DataFrame( + np.c_[portfolios, factors], + columns=[f"p{i}" for i in range(NPORTFOLIO)] + ["f0", "f1"], + ) + formula = " + ".join(f"p{i}" for i in range(NPORTFOLIO)) + " ~ f0 + f1" + with pytest.raises(ValueError, match=r"sigma must be a square array"): + LinearFactorModel.from_formula(formula, frame, sigma=np.eye(3)) + mod = LinearFactorModel.from_formula(formula, frame, sigma=np.eye(NPORTFOLIO)) + assert mod.portfolios.shape[1] == NPORTFOLIO diff --git a/linearmodels/tests/shared/test_io.py b/linearmodels/tests/shared/test_io.py new file mode 100644 index 0000000000..9e6a6eb8d6 --- /dev/null +++ b/linearmodels/tests/shared/test_io.py @@ -0,0 +1,61 @@ +""" +Tests of the number formatters of the summaries +""" + +import numpy as np +import pytest + +from linearmodels.shared.io import _str, pval_format +from linearmodels.system import SUR + + +@pytest.mark.parametrize("value", [np.inf, np.float64("inf")], ids=["float", "numpy"]) +def test_str_positive_infinity(value): + assert _str(value) == " inf" + + +@pytest.mark.parametrize("value", [-np.inf, np.float64("-inf")], ids=["float", "numpy"]) +def test_str_negative_infinity(value): + assert _str(value) == " -inf" + + +def test_str_has_the_width_of_the_other_numbers(): + assert len(_str(np.inf)) == len(_str(-np.inf)) == len(_str(1e-5)) == 8 + assert len(_str(np.nan)) == 8 + + +@pytest.mark.parametrize( + ("value", "expected"), + [ + (np.nan, " "), + (0.0, "0.0000"), + (1.23456, "1.2346"), + (-0.5, "-0.5000"), + (12345.678, "1.235e+04"), + (123456.789, "1.235e+05"), + (1e-5, " 1e-05"), + ], +) +def test_str_finite_values_are_not_changed(value, expected): + assert _str(value) == expected + + +def test_pval_format(): + assert pval_format(np.nan) == " " + assert pval_format(0.04567) == "0.0457" + + +def test_summary_with_a_zero_standard_error(): + # An exact fit has a zero standard error and an infinite t-statistic + rng = np.random.default_rng(0) + x = rng.standard_normal((50, 1)) + equations = { + "a": {"dependent": 2 * x[:, 0], "exog": x}, + "b": {"dependent": x[:, 0] + rng.standard_normal(50), "exog": x}, + } + res = SUR(equations).fit(method="ols") + assert np.asarray(res.std_errors)[0] == 0.0 + assert np.isinf(res.tstats.iloc[0]) + summary = str(res.summary) + assert "inf" in summary + assert "Equation: a" in summary