BUG: Invert the covariance of the pricing errors on its N - K - nrf d… by bashtage · Pull Request #805 · bashtage/linearmodels

bashtage · 2026-10-12T22:09:09Z

…imensions

The J-statistic of LinearFactorModel was alpha' pinv(V) alpha, with V the
estimated covariance of the pricing errors and the default tolerance of pinv.
The second stage estimates make the pricing errors orthogonal to the K + nrf
columns of sigma^-1 [1 betas] (nrf is 1 if the risk-free rate is estimated), so
V has rank N - K - nrf in theory and the statistic has N - K - nrf degrees of
freedom. pinv also inverts the K + nrf directions in which V is null in theory.
Their eigenvalues are small but not zero in a sample, and the statistic

depended on the order of the rows with risk_free=True, where one eigenvalue is
numerically zero (1e-14 of the largest) and was inverted: 7.329 to 7.364 over
12 orderings of the same data, 6.8639 now for all of them;
was too large and rejected a correct model too often. In 200 simulated
samples with 250 observations, 25 portfolios and 3 factors the old statistic
was larger by a median of 4% and by up to 30%, and the 5% test rejected 12.5%
of the samples (10.0% now; 7.5% and 5.5% with 1,000 observations);
can be far from the statistic of the pricing errors when there are few
portfolios. For the five portfolios and three factors of the Fama-French
example of the documentation (819 months) it was 100.6 and is 41.6
(32.6 and 11.7 with risk_free=True).

The statistic is now the Wald statistic of the pricing errors in the space in
which they lie, (Q'a)'(Q'VQ)^-1(Q'a), where the columns of Q are an orthonormal
basis of the space that is orthogonal to the columns of sigma^-1 [1 betas]. This
is alpha' V+ alpha for the generalized inverse V+ = Q (Q'VQ)^-1 Q' of V, which
has rank N - K - nrf and the null space that the restrictions imply. It is the
Moore-Penrose pseudo-inverse when V is exactly singular in those directions.
Inverting the N - K - nrf largest eigenvalues of V instead agrees with it when V
is close to singular in those directions (to 0.2% in the simulations above), but
not otherwise: with sigma given (GLS) it rejected 2.3% of 300 samples of a
correct model at the nominal 5% (250 observations, 10 portfolios, 2 factors,
8.7% for this statistic and 11.0% for pinv), and with the Fama-French example it
gives 13.4.

An InvalidTestStatistic, which the summary can show, is returned if Q'VQ is
singular, as it is when there are too few observations for the portfolios. The
old statistic gave numbers such as -4.4e15 then. With as many portfolios as risk
premia the statistic is exactly 0.

The three tests of test_linear_factor_model.py that re-implemented the pinv
formula now compute the statistic from the restrictions with an independent
null space (scipy.linalg.null_space), and no longer need a tolerance of 10% for
the risk-free models. test_j_statistic.py checks the statistic for the
combinations of risk_free, sigma, cov_type and debiased, the invariance to the
order of the rows, the just identified model, and the helper for covariances
that are, and are not, singular in the restricted directions.

The docstrings of LinearFactorModel and of the J-statistic of the results, and
the mathematical detail of the documentation, describe the statistic and the
change.

…imensions The J-statistic of LinearFactorModel was alpha' pinv(V) alpha, with V the estimated covariance of the pricing errors and the default tolerance of pinv. The second stage estimates make the pricing errors orthogonal to the K + nrf columns of sigma^-1 [1 betas] (nrf is 1 if the risk-free rate is estimated), so V has rank N - K - nrf in theory and the statistic has N - K - nrf degrees of freedom. pinv also inverts the K + nrf directions in which V is null in theory. Their eigenvalues are small but not zero in a sample, and the statistic - depended on the order of the rows with risk_free=True, where one eigenvalue is numerically zero (1e-14 of the largest) and was inverted: 7.329 to 7.364 over 12 orderings of the same data, 6.8639 now for all of them; - was too large and rejected a correct model too often. In 200 simulated samples with 250 observations, 25 portfolios and 3 factors the old statistic was larger by a median of 4% and by up to 30%, and the 5% test rejected 12.5% of the samples (10.0% now; 7.5% and 5.5% with 1,000 observations); - can be far from the statistic of the pricing errors when there are few portfolios. For the five portfolios and three factors of the Fama-French example of the documentation (819 months) it was 100.6 and is 41.6 (32.6 and 11.7 with risk_free=True). The statistic is now the Wald statistic of the pricing errors in the space in which they lie, (Q'a)'(Q'VQ)^-1(Q'a), where the columns of Q are an orthonormal basis of the space that is orthogonal to the columns of sigma^-1 [1 betas]. This is alpha' V+ alpha for the generalized inverse V+ = Q (Q'VQ)^-1 Q' of V, which has rank N - K - nrf and the null space that the restrictions imply. It is the Moore-Penrose pseudo-inverse when V is exactly singular in those directions. Inverting the N - K - nrf largest eigenvalues of V instead agrees with it when V is close to singular in those directions (to 0.2% in the simulations above), but not otherwise: with sigma given (GLS) it rejected 2.3% of 300 samples of a correct model at the nominal 5% (250 observations, 10 portfolios, 2 factors, 8.7% for this statistic and 11.0% for pinv), and with the Fama-French example it gives 13.4. An InvalidTestStatistic, which the summary can show, is returned if Q'VQ is singular, as it is when there are too few observations for the portfolios. The old statistic gave numbers such as -4.4e15 then. With as many portfolios as risk premia the statistic is exactly 0. The three tests of test_linear_factor_model.py that re-implemented the pinv formula now compute the statistic from the restrictions with an independent null space (scipy.linalg.null_space), and no longer need a tolerance of 10% for the risk-free models. test_j_statistic.py checks the statistic for the combinations of risk_free, sigma, cov_type and debiased, the invariance to the order of the rows, the just identified model, and the helper for covariances that are, and are not, singular in the restricted directions. The docstrings of LinearFactorModel and of the J-statistic of the results, and the mathematical detail of the documentation, describe the statistic and the change.

…s of the GMM factor model The standard errors of the alphas of LinearFactorModelGMM were the square roots of the diagonal of S / T, with S = g'g / T the uncentered second moment of the moment conditions. The alphas are the averages of some of those moment conditions at the estimates, and S / T is their covariance at the true parameters. It ignores that the parameters were estimated, and did not depend on cov_type, center, kernel, bandwidth or debiased. The standard errors were too large, and the t-statistics of the alphas too small. With efficient GMM the covariance of the average moment conditions at the estimates is (S - G (G'S^-1 G)^-1 G') / (T - df), where S is the covariance of the moment conditions, G the Jacobian, and df the number of factors if the covariance is debiased and 0 otherwise. The standard errors of the alphas are now the square roots of its diagonal, with the S, G and scale that the covariance of the parameters, (G'S^-1 G)^-1 / (T - df), uses, so that the robust and kernel estimators, center and debiased apply to them. The covariance of the alphas has rank N - K - nrf, the degrees of freedom of the J-statistic. The covariance of the parameters, the estimates, the risk premia and the J-statistic do not change. - In 60 simulated samples with 400 observations, 6 portfolios and 2 factors and no pricing errors, the sampling standard deviation of the alphas is 0.00047 to 0.00136. The mean standard error was 0.00143 to 0.00145 for all six alphas (3 times too large for the portfolio with the smallest standard deviation), and is 0.00049 to 0.00133. - In a simulation with a portfolio that has a large beta, the ratio of the mean standard error to the sampling standard deviation of the alphas of the four portfolios was 1.03, 1.27, 1.03 and 2.60, and is 1.01, 1.22, 0.97 and 1.03 (50 samples with 600 observations). - For the five portfolios and three factors of the Fama-French example of the documentation (819 months), the standard errors of the alphas of the two-step estimator were 0.00123, 0.00081, 0.00060, 0.00130 and 0.00085, and are 0.00018, 0.00003, 0.00011, 0.00004 and 0.00003. The t-statistics were -0.66, 0.24, -0.16, -0.10 and -0.13, and are -4.53, 5.58, -0.88, -2.93 and -3.53. If the model is just identified, with as many portfolios as risk premia, the alphas are zero and so is their variance. The standard errors of the alphas were about 0.001 and are exactly 0 then, and their t-statistics are not finite. _str, which formats the numbers of the summaries, raised an OverflowError for an infinite value, so that full_summary failed for these models. It now formats infinite values, and the t-statistics are computed without a warning. HeteroskedasticCovariance, and so KernelCovariance, has a new moment_cov property that returns this covariance. It is zero if the Jacobian is square. test_gmm_alpha_se.py compares the standard errors and the t-statistics with an independent computation of M S M', with M = I - G (G'S^-1 G)^-1 G'S^-1, for risk_free, the two covariance estimators, center and debiased, checks the rank of the covariance, and runs the simulation above. test_covariance.py checks moment_cov. The docstrings of LinearFactorModelGMM and of the standard errors of the results, and the mathematical detail of the documentation, describe the standard errors and the change.

…covariance With debiased=True, the default, the covariance of all the parameters of TradedFactorModel was divided by nobs - nfactor - 1, the number of observations less the number of regressors of the time-series regressions. The risk premia are the sample means of the factors, so their covariance should use the degrees of freedom of a mean, nobs - 1, and risk_premia_se does: it is computed from a separate covariance of the demeaned factors with a divisor of nobs - 1. The block of cov for the risk premia was therefore larger than risk_premia_se squared, by (nobs - 1) / (nobs - nfactor - 1): 5.4% with 60 observations and 3 factors, and 0.4% in the Fama-French example of the documentation (819 months, 3 factors). The risk premia block of cov now uses nobs - 1 as well, so that it equals risk_premia_se squared, and is the variance of a sample mean, s^2 / nobs with s^2 the unbiased variance, for the robust covariance. The alphas and the loadings still use nobs - nfactor - 1, which is what the covariance of the OLS estimates uses, and nothing changes with debiased=False. The covariance between the risk premia and the other parameters is divided by the geometric mean of the two divisors, sqrt((nobs - nfactor - 1)(nobs - 1)). This rescales the standard errors of each group by the square root of its own factor, so that the covariance matrix remains positive semi-definite. The estimates, risk_premia_se and the J-statistic do not change. test_traded_risk_premia_cov.py compares the risk premia block with the variance of a sample mean for the robust and the kernel covariance, with and without debiasing, checks the scaling of each block against the covariance with debiased=False, and that the matrix is symmetric and positive semi-definite. test_model.py, which compared the whole covariance with a divisor of nobs - nfactor - 1, applies the new scale to the risk premia. The docstrings of TradedFactorModel.fit and of risk_premia_se, and the mathematical detail of the documentation, describe the scaling and the change.

InvalidTestStatistic is used by the J-statistic helper only, through the module instead of by name. There is no change in behavior. It leaves the lines that import from linearmodels.shared.exceptions and from linearmodels.shared.hypotheses unchanged, which another change edits.

The example in the introduction of the asset pricing documentation estimates a LinearFactorModel with the kernel covariance on six Fama-French portfolios and three factors. The J-statistic that it prints is 30.291 and was 39.109, with the statistic of the previous commits. The R-squared, the risk premia, their standard errors and the other values of the output are unchanged.

… test The test checks that the summary of a model with a zero standard error is computed without a warning from the division. Only the RuntimeWarning of a division matters, and other warnings, such as deprecations of other packages, should not fail the test.

…le notebook The notebook of the asset pricing examples says that the risk-free rate has little impact on the J-statistic of the two-step estimator, and that the model is still firmly rejected with a kernel covariance. With the J-statistic of the previous commits the statistic of the three factor model with nine portfolios is 24.9 (p-value 0.0004) and falls to 19.7 with the risk-free rate (0.0014) and to 14.6 with the kernel covariance (0.0240), and it was 36.4, 36.0 and 25.8. The two sentences only say that the model is rejected.

…uction The introduction wrote the J-statistics of the factor models without the inverse of the covariance of the pricing errors. It now has the inverse for the traded factor model, and the generalized inverse of rank P - K - nrf for the two-step estimator, which has P - K - nrf degrees of freedom when the risk-free rate is estimated.

…ortfolios The J-statistic of LinearFactorModel is the Wald statistic of the pricing errors in the space that is orthogonal to the columns of sigma^-1 [1 betas]. It depends on that space and not on the basis Q that is used to compute it, so that reordering the test portfolios, and sigma in the same way, reorders the pricing errors, their covariance and the restrictions in the same way and leaves the statistic unchanged. The tests check this for the robust and kernel covariance estimators, with and without GLS weights and with and without an estimated risk-free rate, and that the helper gives the same statistic for any basis of the space. The existing test covers the order of the observations.

The docstring of LinearFactorModel now says what the J-statistic is and why it is computed this way: * the second stage estimates make the pricing errors orthogonal to the K + nrf columns of sigma^-1 [1 betas], so they lie in a space with N - K - nrf dimensions, and the covariance of the pricing errors has that rank in theory; * the statistic is the Wald statistic of the coordinates Q'a of the pricing errors in that space, (Q'a)'(Q'VQ)^-1(Q'a), where Q has orthonormal columns that span the space that is orthogonal to the columns of sigma^-1 [1 betas]; * the estimated covariance is not exactly singular in a sample, so that its pseudo-inverse would either invert the directions that are null in theory or need a tolerance to decide which eigenvalues are zero. The projection uses the rank that the model implies and needs no tolerance; * the statistic depends on the space and not on the basis Q of the space, so that it does not depend on the order of the test portfolios (with sigma reordered in the same way). It does not depend on the order of the observations for the robust covariance estimator, but the kernel estimators use that order; * the statistic is an InvalidTestStatistic if Q'VQ is numerically singular, and it is 0 with 0 degrees of freedom if the model is just identified. The j_statistic property of the results and the helper that computes the statistic describe it in the same terms.

moment_cov, the covariance of the average moment conditions at the estimates, was a new public property of HeteroskedasticCovariance and KernelCovariance. It is only needed for the standard errors of the alphas of LinearFactorModelGMM. It is now _moment_cov, so that 7.1 does not add a public name for it. The estimates, the standard errors and the tests of the values are unchanged. A new test checks that the name is not public.

…unding The alphas of a LinearFactorModel that is just identified, with as many portfolios as risk premia, are zero and so are their variances. Rounding makes those variances tiny numbers of either sign, about -3e-19 in most samples. The standard error of a negative variance was NaN, with "RuntimeWarning: invalid value encountered in sqrt" in std_errors and in full_summary. The standard errors of the asset pricing results are now computed from the diagonal of the covariance with a helper. A variance that is negative by no more than 1e-10 times the largest variance of the covariance is zero up to rounding and its standard error is 0. In 1,124 simulated just identified models the most negative variance was 2e-13 times the largest, and no variance of the over identified models was negative. A variance that is more negative is not hidden: its standard error is NaN, with the RuntimeWarning, as before. The tolerance is relative, so that it does not depend on the units of the returns. It applies to std_errors and risk_premia_se, and to the alphas and the loadings of LinearFactorModelGMM. The tests cover the helper at the boundary of the tolerance, for several scales, the standard errors of just identified models with the robust and kernel covariance estimators, variances that are negative by rounding and by a lot in std_errors, risk_premia_se and the alphas of the GMM model, and that the standard errors of the other variances are unchanged.

… of 0 The alphas of a model that is just identified are zero by construction, up to rounding of about 1e-17, and their standard errors are 0: exactly for LinearFactorModelGMM and, since the standard errors of the variances that are zero up to rounding are 0, for most samples of LinearFactorModel. The t-statistic of an alpha was then plus or minus infinity and its p-value 0, which reads as a highly significant pricing error: alpha 3.886e-18 0.0000 inf 0.0000 3.886e-18 3.886e-18 The test of a parameter that is known is not defined. The t-statistic of a parameter with a standard error of 0 is now NaN, and so is its p-value. The summaries show blanks, as they do for any NaN. The risk premia t-statistics use the same rule. The t-statistics and p-values of the parameters with a standard error that is not 0 are unchanged.

codecov · 2026-10-13T01:22:49Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.66%. Comparing base (9fa3a90) to head (43f1b41).
⚠️ Report is 101 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #805      +/-   ##
==========================================
+ Coverage   99.60%   99.66%   +0.06%     
==========================================
  Files         136      165      +29     
  Lines       26057    31180    +5123     
  Branches     1873     2113     +240     
==========================================
+ Hits        25954    31077    +5123     
- Misses         38       40       +2     
+ Partials       65       63       -2

Flag	Coverage Δ
adder	`99.66% <100.00%> (+0.06%)`	⬆️
subtractor	`99.66% <100.00%> (+0.06%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

test_j_statistic_textbook also asserted that the new statistic is smaller than alpha' pinv(V) alpha, the statistic of 7.0. With risk_free=True the covariance of the alphas has an eigenvalue that is zero in theory and of the size of rounding errors in practice. Whether np.linalg.pinv keeps it depends on the rounding of the LAPACK build. On the Linux jobs with the newest NumPy (python314_latest, python314_no_binary, python315_latest) it falls below the tolerance of pinv, the two statistics agree to rounding, and the strict inequality failed for the kernel covariance without sigma. The inequality describes the old pseudo-inverse, not the statistic. The test still checks the value against the textbook computation, the degrees of freedom and the p-value.

CodeQL reported py/side-effect-in-assert for the call to np.delete inside the assert of test_std_errors_not_changed. np.delete returns a copy, so there was no side effect, but the call is now made before the assert.

bashtage added 13 commits October 12, 2026 20:34

github-advanced-security Bot found potential problems Oct 13, 2026

View reviewed changes

Comment thread linearmodels/tests/asset_pricing/test_std_errors_rounding.py Fixed

bashtage added 2 commits October 13, 2026 07:23

TST: Move a call out of an assert in the standard error rounding test

43f1b41

CodeQL reported py/side-effect-in-assert for the call to np.delete inside the assert of test_std_errors_not_changed. np.delete returns a copy, so there was no side effect, but the call is now made before the assert.

bashtage merged commit 181683d into main Oct 13, 2026
24 checks passed

bashtage deleted the lm71-fix-factor-model-inference branch October 13, 2026 08:05

Navigation Menu

Repository navigation

BUG: Invert the covariance of the pricing errors on its N - K - nrf d… - #805

BUG: Invert the covariance of the pricing errors on its N - K - nrf d…#805
bashtage merged 15 commits into
mainfrom
lm71-fix-factor-model-inference

bashtage commented Oct 12, 2026

Uh oh!

codecov Bot commented Oct 13, 2026 •

edited

Loading

Uh oh!

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

bashtage commented Oct 12, 2026

Uh oh!

codecov Bot commented Oct 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

codecov Bot commented Oct 13, 2026 •

edited

Loading