From 2210df45959881958b023928d3c3ea9a76d0e911 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:11:58 +0100 Subject: [PATCH 1/8] BUG: Do not treat badly scaled moment covariances as singular in centered_j_stat IVGMMResults.centered_j_stat returned an InvalidTestStatistic, saying that the covariance of the moment conditions is singular, for models in which it is not. It used matrix_rank on the covariance, and the tolerance of matrix_rank is relative to the largest eigenvalue. When the variances of two moment conditions differ by a factor of about 1e16, which is the case if two instruments are measured in units that differ by a factor of 1e8, the smaller one is below the tolerance and the covariance is reported as singular. For a model with four instruments, one of which is multiplied by a constant, the J statistic is 7.6196 at every scale. The centered statistic is 7.906 for constants up to 1e7 and is invalid from 1e8 on, and the same for constants down to 1e-8. This occurred for the unadjusted, robust, kernel and clustered estimators, with and without weights, and when the scale was that of an exogenous regressor, which is also an instrument. The rank is now that of the correlation matrix of the moment conditions, which does not depend on the units. A covariance whose variances are not positive or that is not finite is also singular. The statistic is computed as before, and all of the existing reference values are unchanged. A covariance that is singular is still found, e.g., when there are too few clusters, with any scale. The tests that are added check the definition of the statistic for IVGMM and IVGMMCUE, without using the weight matrix classes, with weights, one-way clusters, a kernel and the small sample adjustment. The reference values from R do not cover continuously updating GMM with weights or a kernel, or weights with clusters. They also check continuously updating GMM with and without centering, and that the statistic is invalid when there are as many clusters as moment conditions. 32 tests of the units of the data fail without the change. --- linearmodels/iv/model.py | 15 +- linearmodels/iv/results.py | 8 +- .../tests/iv/test_new_feature_hardening.py | 251 ++++++++++++++++++ 3 files changed, 272 insertions(+), 2 deletions(-) create mode 100644 linearmodels/tests/iv/test_new_feature_hardening.py diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 5643020c30..550a8885dd 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -17,15 +17,18 @@ average, c_, column_stack, + diag, dtype, eye, float64, + isfinite, isscalar, logical_not, nan, nanmean, ndarray, ones, + outer, sqrt, squeeze, ) @@ -1070,7 +1073,17 @@ def _centered_j_statistic( weight_matrix_estimator = WEIGHT_MATRICES[self._weight_type] centered = weight_matrix_estimator(**{**self._weight_kwargs, "center": True}) s_c = centered.weight_matrix(x, z, eps) - if matrix_rank(s_c, hermitian=True) < ninstr: + # The rank is that of the correlation matrix of the moment conditions, + # so that it does not depend on the units in which the instruments are + # measured. The tolerance of the rank test is relative to the largest + # eigenvalue, and so it treats the covariance as singular if the + # variances of the moment conditions differ by a factor of about 1e16 + variances = diag(s_c) + singular = not (npall(isfinite(s_c)) and npall(variances > 0)) + if not singular: + std = sqrt(variances) + singular = matrix_rank(s_c / outer(std, std), hermitian=True) < ninstr + if singular: # The centered cluster sums add to zero, so a clustered estimator # has rank at most one less than the number of clusters reason = ( diff --git a/linearmodels/iv/results.py b/linearmodels/iv/results.py index c02e1cbcb8..6ba2912d71 100644 --- a/linearmodels/iv/results.py +++ b/linearmodels/iv/results.py @@ -1656,7 +1656,7 @@ def centered_j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: Returns ------- - WaldTestStatistic + WaldTestStatistic or InvalidTestStatistic J statistic test of overidentifying restrictions that uses a covariance of the moment conditions that has been centered. An InvalidTestStatistic is returned if that covariance is @@ -1707,6 +1707,12 @@ def centered_j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: It also has no effect in :class:`~linearmodels.iv.model.IVGMMCUE` models unless ``center=False`` was set. + In a weighted model the moment conditions are + :math:`w_i z_i \hat{\epsilon}_i`, where :math:`w_i` are the weights of + the model, as in :attr:`j_stat`. The statistic and the test of whether + the covariance is singular do not depend on the units in which the + variables are measured. + References ---------- Hall, A. R. (2000). Covariance matrix estimation and the power of the diff --git a/linearmodels/tests/iv/test_new_feature_hardening.py b/linearmodels/tests/iv/test_new_feature_hardening.py new file mode 100644 index 0000000000..2d0ac0b954 --- /dev/null +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -0,0 +1,251 @@ +""" +Checks of the instrumental variable features that were added after version 7.0 + +* IVGMMResults.centered_j_stat: the units of the data, continuously updating + GMM, weights, clusters, kernels and singular covariances. +""" + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.iv import IVGMM, IVGMMCUE +from linearmodels.iv.results import IVGMMResults +from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic + + +def simulate(seed, nobs=300, ninstr=4, nendog=1, violation=0.0, weighted=False): + """ + Heteroskedastic IV data with a constant and an exogenous regressor + + The last instrument is invalid if violation is not 0. The weights are None + unless weighted. + """ + rng = np.random.default_rng(seed) + z = rng.standard_normal((nobs, ninstr)) + v = rng.standard_normal((nobs, nendog)) + u = (1 + 0.5 * np.abs(z[:, 0])) * (0.5 * v.sum(1) + rng.standard_normal(nobs)) + endog = z @ rng.uniform(0.3, 0.8, (ninstr, nendog)) + v + exog = np.column_stack([np.ones(nobs), rng.standard_normal(nobs)]) + y = exog @ [1.0, 0.5] + endog.sum(1) + violation * z[:, -1] + u + w = rng.uniform(0.3, 3.0, nobs) if weighted else None + return y, exog, endog, z, w + + +# ---------------------------------------------------------------------------- +# IVGMMResults.centered_j_stat +# ---------------------------------------------------------------------------- + + +def moment_conditions(res, y, exog, endog, z, w): + """Moment conditions of the weighted model at the estimates""" + sw = np.ones(y.shape[0]) if w is None else np.sqrt(w / w.mean()) + x = np.column_stack([exog, endog]) * sw[:, None] + instruments = np.column_stack([exog, z]) * sw[:, None] + eps = y * sw - x @ res.params.to_numpy() + return instruments * eps[:, None] + + +def centered_covariance( + g, weight_type, nvar, *, debiased=False, clusters=None, bandwidth=None +): + """The covariance of the moments about their mean, from the definition""" + nobs = g.shape[0] + gc = g - g.mean(0) + if weight_type == "robust": + s = gc.T @ gc / nobs + scale = nobs / (nobs - nvar) + elif weight_type == "kernel": + s = gc.T @ gc + for lag in range(1, bandwidth + 1): + gamma = gc[lag:].T @ gc[:-lag] + s = s + (1 - lag / (bandwidth + 1)) * (gamma + gamma.T) + s = s / nobs + scale = nobs / (nobs - nvar) + else: + sums = pd.DataFrame(gc).groupby(clusters).sum().to_numpy() + s = sums.T @ sums / nobs + ngroups = sums.shape[0] + scale = (nobs - 1) / (nobs - nvar) * ngroups / (ngroups - 1) + return s * scale if debiased else s + + +def weight_options(weight_type, nobs): + if weight_type == "kernel": + return {"kernel": "bartlett", "bandwidth": 3} + if weight_type == "clustered": + return {"clusters": np.arange(nobs) // 10} + return {} + + +def expected_centered_j(res, data, weight_type, debiased, options): + y, exog, endog, z, w = data + g = moment_conditions(res, y, exog, endog, z, w) + nvar = exog.shape[1] + endog.shape[1] + s_c = centered_covariance( + g, + weight_type, + nvar, + debiased=debiased, + clusters=options.get("clusters"), + bandwidth=options.get("bandwidth"), + ) + gbar = g.mean(0) + return g.shape[0] * gbar @ np.linalg.solve(s_c, gbar) + + +@pytest.mark.parametrize("debiased", [False, True]) +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize("weight_type", ["robust", "kernel", "clustered"]) +def test_centered_j_gmm_matches_the_definition(weight_type, weighted, debiased): + data = simulate(20, violation=0.3, weighted=weighted) + y, exog, endog, z, w = data + options = weight_options(weight_type, y.shape[0]) + mod = IVGMM( + y, + exog, + endog, + z, + weights=w, + weight_type=weight_type, + debiased=debiased, + **options, + ) + res = mod.fit() + expected = expected_centered_j(res, data, weight_type, debiased, options) + assert type(res.centered_j_stat) is WaldTestStatistic + assert res.centered_j_stat.df == 3 + assert_allclose(res.centered_j_stat.stat, expected, rtol=1e-9) + # The violation of the moment conditions makes the statistics differ + assert not np.isclose(res.centered_j_stat.stat, res.j_stat.stat, rtol=1e-3) + + +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize("weight_type", ["robust", "kernel", "clustered"]) +def test_centered_j_cue_matches_the_definition(weight_type, weighted): + data = simulate(21, violation=0.3, weighted=weighted) + y, exog, endog, z, w = data + options = weight_options(weight_type, y.shape[0]) + mod = IVGMMCUE(y, exog, endog, z, weights=w, weight_type=weight_type, **options) + res = mod.fit() + assert isinstance(res, IVGMMResults) + expected = expected_centered_j(res, data, weight_type, False, options) + assert_allclose(res.centered_j_stat.stat, expected, rtol=1e-8) + # The CUE centers by default, and so the weight matrix of the objective is + # the covariance of the statistic + assert res.weight_config["center"] is True + assert_allclose(res.centered_j_stat.stat, res.j_stat.stat, rtol=1e-8) + assert res.centered_j_stat.df == res.j_stat.df == 3 + + +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize("weight_type", ["robust", "kernel", "clustered"]) +def test_centered_j_cue_without_centering(weight_type, weighted): + data = simulate(22, violation=0.3, weighted=weighted) + y, exog, endog, z, w = data + options = weight_options(weight_type, y.shape[0]) + mod = IVGMMCUE( + y, exog, endog, z, weights=w, weight_type=weight_type, center=False, **options + ) + res = mod.fit() + expected = expected_centered_j(res, data, weight_type, False, options) + assert res.weight_config["center"] is False + assert_allclose(res.centered_j_stat.stat, expected, rtol=1e-8) + assert abs(res.centered_j_stat.stat - res.j_stat.stat) > 1e-3 + if weight_type == "robust": + # S = S_c + gbar gbar' at the estimates, as in iterated GMM + nobs = y.shape[0] + j = res.j_stat.stat + assert_allclose(res.centered_j_stat.stat, j / (1 - j / nobs), rtol=1e-8) + + +@pytest.mark.parametrize("weighted", [False, True]) +def test_centered_j_cue_with_too_few_clusters(weighted): + # The centered cluster sums add to zero, so that the covariance has at most + # nclusters - 1 non-zero eigenvalues, and the model has 5 moment conditions. + # The uncentered covariance of the CUE has nclusters + y, exog, endog, z, w = simulate(23, ninstr=3, violation=0.3, weighted=weighted) + outcomes = {} + for nclusters in (5, 6): + mod = IVGMMCUE( + y, + exog, + endog, + z, + weights=w, + weight_type="clustered", + clusters=np.arange(y.shape[0]) % nclusters, + center=False, + ) + res = mod.fit() + outcomes[nclusters] = res + assert np.isfinite(res.j_stat.stat) + assert res.j_stat.df == 2 + invalid = outcomes[5].centered_j_stat + assert isinstance(invalid, InvalidTestStatistic) + assert np.isnan(invalid.stat) + assert np.isnan(invalid.pval) + assert invalid.critical_values is None + assert "singular" in str(invalid) + assert "Centered J-test" in str(invalid) + valid = outcomes[6].centered_j_stat + assert type(valid) is WaldTestStatistic + assert valid.df == 2 + assert np.isfinite(valid.stat) + + +@pytest.mark.parametrize("scale", [1e9, 1e-9]) +@pytest.mark.parametrize("column", ["instrument", "exog"]) +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize("weight_type", ["robust", "unadjusted", "kernel", "clustered"]) +def test_centered_j_does_not_depend_on_the_units_of_the_data( + weight_type, weighted, column, scale +): + # Measuring a variable in other units multiplies a moment condition by a + # constant. The J statistic does not change, and the covariance of the + # moment conditions is not singular. The test of its rank used the largest + # eigenvalue of the covariance and so rejected any covariance for which the + # variances of the moment conditions differed by a factor of 1e16 + y, exog, endog, z, w = simulate(24, violation=0.3, weighted=weighted) + options = weight_options(weight_type, y.shape[0]) + + def fit(exog_in, z_in): + mod = IVGMM( + y, exog_in, endog, z_in, weights=w, weight_type=weight_type, **options + ) + return mod.fit() + + base = fit(exog, z) + exog_in, z_in = exog.copy(), z.copy() + if column == "instrument": + z_in[:, 1] *= scale + else: + exog_in[:, 1] *= scale + res = fit(exog_in, z_in) + assert type(base.centered_j_stat) is WaldTestStatistic + assert_allclose(res.j_stat.stat, base.j_stat.stat, rtol=1e-6) + assert type(res.centered_j_stat) is WaldTestStatistic + assert res.centered_j_stat.df == base.centered_j_stat.df + assert_allclose(res.centered_j_stat.stat, base.centered_j_stat.stat, rtol=1e-6) + assert_allclose(res.centered_j_stat.pval, base.centered_j_stat.pval, rtol=1e-4) + + +def test_centered_j_singular_covariance_is_still_found_with_other_units(): + # The check of the rank does not accept a covariance that is singular. There + # are 5 moment conditions and 5 clusters, so that the covariance that is + # centered has rank 4 + y, exog, endog, z, _ = simulate(25, ninstr=3, violation=0.3) + z = z.copy() + z[:, 1] *= 1e9 + mod = IVGMM( + y, + exog, + endog, + z, + weight_type="clustered", + clusters=np.arange(y.shape[0]) % 5, + ) + res = mod.fit() + assert isinstance(res.centered_j_stat, InvalidTestStatistic) + assert "singular" in str(res.centered_j_stat) From 36d2bffd4111b507b1e45d54f80f984ac3173d98 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:33:34 +0100 Subject: [PATCH 2/8] BUG: Select the kernel bandwidth of IVJIVE using the data The covariance of IVJIVE is computed with the weight matrix classes of the GMM models. With cov_type="kernel" and no bandwidth they use a bandwidth of nobs - 2, which is not a consistent estimator, and the standard errors are too small. IV2SLS and IVLIML select the bandwidth using the data, and the documentation of IVJIVE says that its covariance estimators are the same as in the other IV models. For one simulated data set with 500 observations, 4 instruments and independent errors, the standard error of the coefficient on the endogenous variable is 0.050 with the default, 0.065 for IV2SLS and 0.066 for the robust covariance estimator of IVJIVE. Coverage of nominal 95% confidence intervals for the same coefficient, using 400 simulated data sets with 4 instruments, so that the standard error of the frequencies is 0.012: nobs errors IV2SLS, kernel IVJIVE, nobs - 2 IVJIVE, selected 250 iid 0.915 0.693 0.927 500 iid 0.925 0.667 0.943 500 AR(1), 0.5 0.932 0.703 0.943 1000 AR(1), 0.5 0.938 0.620 0.943 The bandwidth is now selected with the option of the kernel weight matrix that selects it using the scores, unless it is given. The bandwidth that is used is in cov_config of the results, as before, and a bandwidth of nobs - 2 can be used by passing it. IVGMM and IVGMMCUE are not changed. They document the default of nobs - 2, but it has the same effect there: the coverage of IVGMM is about 0.66 in the simulation of the second row. The first-stage tests of the rank of their results use that bandwidth when the covariance is a kernel. Two other options of the IVJIVE covariance estimators gave errors that were poor, and now give the errors of the other IV models: * An unknown kernel raised a KeyError. It now raises the ValueError of IV2SLS and IVLIML that the kernel is unknown, whether or not the bandwidth is given. * Two-way clusters, which IV2SLS and IVLIML support, raised a ValueError from NumPy about the dimensions of arrays. IVJIVE supports one-way clustering, and it now raises a ValueError that says so. A 2-d array with one column is accepted, as in the other models. The tests of the default bandwidth of IVJIVE in test_jive.py required nobs - 2, and are updated. The new tests check that the bandwidth is selected, with weights and for each kernel, that the first-stage tests of the rank use it, that it depends on the data, the coverage of the intervals, and the errors above. --- linearmodels/iv/model.py | 26 ++++++- linearmodels/tests/iv/test_jive.py | 117 +++++++++++++++++++++++++++-- 2 files changed, 134 insertions(+), 9 deletions(-) diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 550a8885dd..4116cf54c4 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -39,6 +39,7 @@ from linearmodels.iv._utility import IVFormulaParser from linearmodels.iv.common import f_statistic, find_constant from linearmodels.iv.covariance import ( + KERNEL_LOOKUP, ClusteredCovariance, HeteroskedasticCovariance, HomoskedasticCovariance, @@ -1892,6 +1893,21 @@ def __init__( ) -> None: # The weight matrix has no effect in a just identified model w = eye(x.shape[1]) + if cov_type == "kernel": + kernel = cov_config.get("kernel", "bartlett") + if kernel not in KERNEL_LOOKUP: + raise ValueError(f"Unknown kernel: {kernel}") + # Select the bandwidth using the data unless it is given, as in + # IV2SLS and IVLIML. The default of the weight matrix estimator is + # nobs - 2, which is not a consistent covariance estimator + cov_config = {"optimal_bw": True, **cov_config} + elif cov_type == "clustered": + clusters = cov_config.get("clusters") + if clusters is not None and asarray(clusters).squeeze().ndim > 1: + raise ValueError( + "IVJIVE supports one-way clustering only. clusters must be " + "a 1-d array, or a 2-d array with a single column." + ) super().__init__(x, y, x_loo, params, w, cov_type, debiased, **cov_config) self._name = "JIVE Covariance" if cov_type in ("unadjusted", "homoskedastic"): @@ -2092,9 +2108,13 @@ def fit( ----- The covariance estimators are the same as in the other IV models with the leave-one-out first-stage predictions as the instruments. The - additional parameters depend on the estimator: ``kernel`` and - ``bandwidth`` for "kernel" and ``clusters`` for "clustered". The - defaults are used if none are provided. + additional parameters depend on the estimator: ``kernel`` (the default + is "bartlett") and ``bandwidth`` for "kernel" and ``clusters`` for + "clustered". If ``bandwidth`` is not provided it is selected using the + data, as in :class:`~linearmodels.iv.model.IV2SLS`, and the bandwidth + that was used is in ``cov_config`` of the results. "clustered" is + one-way clustering, and ``clusters`` has one element for each + observation. Since JIVE is not a k-class estimator, ``kappa`` is None in the results. diff --git a/linearmodels/tests/iv/test_jive.py b/linearmodels/tests/iv/test_jive.py index ece7abb1d3..080d332643 100644 --- a/linearmodels/tests/iv/test_jive.py +++ b/linearmodels/tests/iv/test_jive.py @@ -491,16 +491,115 @@ def test_cov_type_aliases(alias): @pytest.mark.parametrize("kernel", ["bartlett", "parzen", "qs"]) def test_kernel_options(kernel): y, exog, endog, z = simulate(29) - default = IVJIVE(y, exog, endog, z).fit(cov_type="kernel", kernel=kernel) - again = IVJIVE(y, exog, endog, z).fit( - cov_type="kernel", kernel=kernel, bandwidth=y.shape[0] - 2 - ) + mod = IVJIVE(y, exog, endog, z) + default = mod.fit(cov_type="kernel", kernel=kernel) + # As in the other IV models, the bandwidth is selected using the data if it + # is not given. It is not nobs - 2, which is not a consistent estimator + chosen = default.cov_config["bandwidth"] + assert 1 <= chosen < y.shape[0] / 4 + again = mod.fit(cov_type="kernel", kernel=kernel, bandwidth=chosen) assert_allclose(default.cov, again.cov, rtol=1e-12) - short = IVJIVE(y, exog, endog, z).fit(cov_type="kernel", kernel=kernel, bandwidth=2) + long = mod.fit(cov_type="kernel", kernel=kernel, bandwidth=y.shape[0] - 2) + assert not np.allclose(default.cov, long.cov) + short = mod.fit(cov_type="kernel", kernel=kernel, bandwidth=2) assert not np.allclose(default.cov, short.cov) assert short.cov_config["kernel"] == kernel +def test_kernel_bandwidth_is_selected_using_the_data_in_a_weighted_model(): + # A bandwidth of nobs - 2 is a covariance estimator that is not consistent, + # and it understates the standard errors + y, exog, endog, z = simulate(6, nobs=500) + w = np.random.default_rng(6).uniform(0.3, 2.0, y.shape[0]) + mod = IVJIVE(y, exog, endog, z, weights=w) + res = mod.fit(cov_type="kernel") + bandwidth = res.cov_config["bandwidth"] + assert 1 <= bandwidth < 40 + assert res.cov_config["kernel"] == "bartlett" + # It is the bandwidth that the model would use if it was given + same = mod.fit(cov_type="kernel", bandwidth=bandwidth) + assert_allclose(res.cov, same.cov, rtol=1e-12) + assert same.cov_config["bandwidth"] == bandwidth + # nobs - 2 is the default of IVGMM, and gives smaller standard errors + long = mod.fit(cov_type="kernel", bandwidth=y.shape[0] - 2) + assert long.cov_config["bandwidth"] == y.shape[0] - 2 + assert (long.std_errors < 0.9 * same.std_errors).any() + # Each kernel selects a bandwidth + for kernel in ("bartlett", "parzen", "qs"): + chosen = mod.fit(cov_type="kernel", kernel=kernel).cov_config["bandwidth"] + assert 1 <= chosen < y.shape[0] / 4 + + +def test_first_stage_uses_the_bandwidth_that_was_selected(): + # As for IV2SLS, the first stage uses the bandwidth of the model, which is + # the one that was selected when it was not given + y, exog, endog, z = simulate(6, nobs=500) + mod = IVJIVE(y, exog, endog, z) + res = mod.fit(cov_type="kernel") + bandwidth = res.cov_config["bandwidth"] + explicit = mod.fit(cov_type="kernel", bandwidth=bandwidth).first_stage + first_stage = res.first_stage + assert np.isfinite(first_stage.kleibergen_paap.stat) + assert_allclose(first_stage.kleibergen_paap.stat, explicit.kleibergen_paap.stat) + assert_allclose(first_stage.kleibergen_paap_f, explicit.kleibergen_paap_f) + long = mod.fit(cov_type="kernel", bandwidth=y.shape[0] - 2).first_stage + assert not np.isclose(long.kleibergen_paap.stat, first_stage.kleibergen_paap.stat) + + +def test_kernel_bandwidth_depends_on_the_data(): + # Autocorrelated errors need a longer bandwidth than independent errors + rng = np.random.default_rng(7) + nobs = 800 + z = rng.standard_normal((nobs, 3)) + v = rng.standard_normal(nobs) + chosen = {} + for rho in (0.0, 0.8): + u = np.zeros(nobs) + innovations = rng.standard_normal(nobs) + 0.5 * v + for t in range(1, nobs): + u[t] = rho * u[t - 1] + innovations[t] + x = z @ [0.6, 0.4, 0.3] + v + y = 1 + x + u + res = IVJIVE(y, np.ones((nobs, 1)), x[:, None], z).fit(cov_type="kernel") + chosen[rho] = res.cov_config["bandwidth"] + assert chosen[0.8] > chosen[0.0] + assert chosen[0.0] < 20 + + +def test_kernel_confidence_intervals_cover(): + # With a bandwidth of nobs - 2, the coverage of nominal 95% intervals was + # about 0.65. A fixed seed is used. The standard error of a coverage + # frequency is about 0.016 for 200 replications + rng = np.random.default_rng(8) + nobs, reps = 250, 200 + covered = np.zeros(reps, dtype=bool) + for rep in range(reps): + z = rng.standard_normal((nobs, 4)) + v = rng.standard_normal(nobs) + x = z @ [0.5, 0.4, 0.3, 0.2] + v + y = 1 + x + 0.5 * v + rng.standard_normal(nobs) + res = IVJIVE(y, np.ones((nobs, 1)), x[:, None], z).fit(cov_type="kernel") + lower, upper = res.conf_int().iloc[-1] + covered[rep] = lower <= 1.0 <= upper + assert 0.88 < covered.mean() < 0.995 + + +def test_two_way_clusters_are_not_supported(): + y, exog, endog, z = simulate(9) + nobs = y.shape[0] + two_way = np.column_stack([np.arange(nobs) // 10, np.arange(nobs) % 7]) + # The other models do support two-way clustering + IV2SLS(y, exog, endog, z).fit(cov_type="clustered", clusters=two_way) + with pytest.raises(ValueError, match="one-way clustering"): + IVJIVE(y, exog, endog, z).fit(cov_type="clustered", clusters=two_way) + # A column is one-way clustering + column = IVJIVE(y, exog, endog, z).fit( + cov_type="clustered", clusters=two_way[:, :1] + ) + flat = IVJIVE(y, exog, endog, z).fit(cov_type="clustered", clusters=two_way[:, 0]) + assert_allclose(column.cov, flat.cov, rtol=1e-12) + + def test_cov_config_is_not_shared_between_fits(): y, exog, endog, z = simulate(30) mod = IVJIVE(y, exog, endog, z) @@ -538,9 +637,15 @@ def test_invalid_cov_type(): def test_unknown_kernel(): + # The error is the one of the other IV models, whether or not the bandwidth + # is selected using the data y, exog, endog, z = simulate(8) - with pytest.raises(KeyError): + with pytest.raises(ValueError, match="Unknown kernel: unknown"): IVJIVE(y, exog, endog, z).fit(cov_type="kernel", kernel="unknown") + with pytest.raises(ValueError, match="Unknown kernel: unknown"): + IVJIVE(y, exog, endog, z).fit(cov_type="kernel", kernel="unknown", bandwidth=3) + with pytest.raises(ValueError, match="Unknown kernel: unknown"): + IV2SLS(y, exog, endog, z).fit(cov_type="kernel", kernel="unknown") def test_options_for_a_different_covariance_estimator_are_not_ignored(): From 5a9e69ba47e1726ca72ce8f36b5ed4518c83f79d Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:33:44 +0100 Subject: [PATCH 3/8] DOC: Remove the unreferenced footnotes from the rank test docstrings The docstrings of cragg_donald, cragg_donald_f, kleibergen_paap and kleibergen_paap_f list their references as footnotes, ".. [1] ...", that the text never cites. Sphinx warns "Footnote [n] is not referenced" for each of them: 3, 1, 2 and 3 warnings for the four functions, which are 9 of the warnings of the documentation build in CI. A build of these docstrings alone with Sphinx 9.1 gives the same 9 warnings at the same lines before the change and none after. The references are a plain list now, as in the docstrings of centered_j_stat and IVJIVE, and the text cites the sources that it uses by author and year: Cragg and Donald (1993) for the statistic, Anderson (1951) for the test of the rank that uses the canonical correlations, and Stock and Yogo (2005) for the critical values. The test finds the footnotes that a docstring defines and does not cite, and the citations that it does not define, for these functions and for the same properties of FirstStageResults. It fails for the four functions without the change, and it keeps the references of each function. --- linearmodels/iv/common.py | 64 ++++++++------- .../tests/iv/test_new_feature_hardening.py | 81 ++++++++++++++++++- 2 files changed, 114 insertions(+), 31 deletions(-) diff --git a/linearmodels/iv/common.py b/linearmodels/iv/common.py index 6317fae137..9af8d910c1 100644 --- a/linearmodels/iv/common.py +++ b/linearmodels/iv/common.py @@ -217,7 +217,7 @@ def cragg_donald( :math:`\chi^2` approximation is poor, and 2SLS can be badly biased and its tests badly sized even when the null is rejected. - The p-value is asymptotic. Stock-Yogo (2005) critical values for weak + The p-value is asymptotic. Stock and Yogo (2005) critical values for weak instruments apply to the F form under i.i.d. errors and are not provided here. Choosing whether to use an instrumental variable estimator based on the outcome of a test like this one distorts the @@ -235,7 +235,8 @@ def cragg_donald( Let :math:`X = Z\Pi + V`, where :math:`\Pi` is the :math:`k \times m` matrix of first-stage coefficients on the excluded instruments, after the exogenous regressors have been partialled out. The null hypothesis - is :math:`\mathrm{rank}(\Pi) < m`. The statistic is + is :math:`\mathrm{rank}(\Pi) < m`. The statistic is the minimum + eigenvalue statistic of Cragg and Donald (1993), .. math:: @@ -246,7 +247,8 @@ def cragg_donald( instruments, and :math:`c` is the number of exogenous regressors. :math:`r^2_{\min}/(1 - r^2_{\min})` is the smallest eigenvalue of :math:`(X'M_ZX)^{-1}X'P_ZX`. The statistic is asymptotically - :math:`\chi^2_{k-m+1}` under the null. + :math:`\chi^2_{k-m+1}` under the null. Anderson (1951) derived a test of + the rank that is a function of the canonical correlations. The statistic is **not** the first-stage F-statistic. It equals :math:`k` times the classical F-statistic when there is a single endogenous @@ -254,15 +256,17 @@ def cragg_donald( References ---------- - .. [1] Cragg, J. G., & Donald, S. G. (1993). Testing identifiability - and specification in instrumental variable models. Econometric - Theory, 9(2), 222-240. - .. [2] Anderson, T. W. (1951). Estimating linear restrictions on - regression coefficients for multivariate normal distributions. - Annals of Mathematical Statistics, 22(3), 327-351. - .. [3] Stock, J. H., & Yogo, M. (2005). Testing for weak instruments in - linear IV regression. In Identification and Inference for Econometric - Models. Cambridge University Press. + Cragg, J. G., & Donald, S. G. (1993). Testing identifiability and + specification in instrumental variable models. *Econometric Theory*, 9(2), + 222-240. + + Anderson, T. W. (1951). Estimating linear restrictions on regression + coefficients for multivariate normal distributions. *Annals of + Mathematical Statistics*, 22(3), 327-351. + + Stock, J. H., & Yogo, M. (2005). Testing for weak instruments in linear IV + regression. In *Identification and Inference for Econometric Models*. + Cambridge University Press. """ name = "Cragg-Donald Test" dec = _first_stage_rank(endog, instr, exog, name) @@ -321,9 +325,9 @@ def cragg_donald_f( References ---------- - .. [1] Stock, J. H., & Yogo, M. (2005). Testing for weak instruments in - linear IV regression. In Identification and Inference for Econometric - Models. Cambridge University Press. + Stock, J. H., & Yogo, M. (2005). Testing for weak instruments in linear IV + regression. In *Identification and Inference for Econometric Models*. + Cambridge University Press. """ test = cragg_donald(endog, instr, exog) if isinstance(test, InvalidTestStatistic): @@ -596,12 +600,12 @@ def kleibergen_paap( References ---------- - .. [1] Kleibergen, F., & Paap, R. (2006). Generalized reduced rank tests - using the singular value decomposition. Journal of Econometrics, - 133(1), 97-126. - .. [2] Kleibergen, F., & Schaffer, M. E. (2007). ranktest: Stata module - to test the rank of a matrix using the Kleibergen-Paap rk statistic. - Statistical Software Components S456865, Boston College. + Kleibergen, F., & Paap, R. (2006). Generalized reduced rank tests using the + singular value decomposition. *Journal of Econometrics*, 133(1), 97-126. + + Kleibergen, F., & Schaffer, M. E. (2007). ranktest: Stata module to test + the rank of a matrix using the Kleibergen-Paap rk statistic. *Statistical + Software Components* S456865, Boston College. """ return _kleibergen_paap(endog, instr, exog, cov_type, cov_config)[0] @@ -688,15 +692,15 @@ def kleibergen_paap_f( References ---------- - .. [1] Kleibergen, F., & Paap, R. (2006). Generalized reduced rank tests - using the singular value decomposition. Journal of Econometrics, - 133(1), 97-126. - .. [2] Olea, J. L. M., & Pflueger, C. (2013). A robust test for weak - instruments. Journal of Business & Economic Statistics, 31(3), - 358-369. - .. [3] Baum, C. F., Schaffer, M. E., & Stillman, S. (2007). Enhanced - routines for instrumental variables/GMM estimation and testing. Stata - Journal, 7(4), 465-506. + Kleibergen, F., & Paap, R. (2006). Generalized reduced rank tests using the + singular value decomposition. *Journal of Econometrics*, 133(1), 97-126. + + Olea, J. L. M., & Pflueger, C. (2013). A robust test for weak instruments. + *Journal of Business & Economic Statistics*, 31(3), 358-369. + + Baum, C. F., Schaffer, M. E., & Stillman, S. (2007). Enhanced routines for + instrumental variables/GMM estimation and testing. *Stata Journal*, 7(4), + 465-506. """ return _kleibergen_paap(endog, instr, exog, cov_type, cov_config)[1] diff --git a/linearmodels/tests/iv/test_new_feature_hardening.py b/linearmodels/tests/iv/test_new_feature_hardening.py index 2d0ac0b954..62898d336c 100644 --- a/linearmodels/tests/iv/test_new_feature_hardening.py +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -1,17 +1,30 @@ """ Checks of the instrumental variable features that were added after version 7.0 +* The tests of the rank of the first-stage coefficient matrix, which are + cragg_donald, cragg_donald_f, kleibergen_paap and kleibergen_paap_f: that the + documentation is well-formed, and the behavior with weights and clusters that + it describes. * IVGMMResults.centered_j_stat: the units of the data, continuously updating GMM, weights, clusters, kernels and singular covariances. """ +import inspect +import re + import numpy as np from numpy.testing import assert_allclose import pandas as pd import pytest from linearmodels.iv import IVGMM, IVGMMCUE -from linearmodels.iv.results import IVGMMResults +from linearmodels.iv.common import ( + cragg_donald, + cragg_donald_f, + kleibergen_paap, + kleibergen_paap_f, +) +from linearmodels.iv.results import FirstStageResults, IVGMMResults from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic @@ -33,6 +46,72 @@ def simulate(seed, nobs=300, ninstr=4, nendog=1, violation=0.0, weighted=False): return y, exog, endog, z, w +# ---------------------------------------------------------------------------- +# Tests of the rank of the first-stage coefficient matrix +# ---------------------------------------------------------------------------- + +DOCUMENTED = { + "cragg_donald": cragg_donald, + "cragg_donald_f": cragg_donald_f, + "kleibergen_paap": kleibergen_paap, + "kleibergen_paap_f": kleibergen_paap_f, + "FirstStageResults.cragg_donald": FirstStageResults.cragg_donald, + "FirstStageResults.cragg_donald_f": FirstStageResults.cragg_donald_f, + "FirstStageResults.kleibergen_paap": FirstStageResults.kleibergen_paap, + "FirstStageResults.kleibergen_paap_f": FirstStageResults.kleibergen_paap_f, +} + + +def footnote_problems(doc): + """Footnotes that are defined and not cited, and citations not defined""" + defined = set(re.findall(r"^\s*\.\. \[(\w+)\]", doc, flags=re.MULTILINE)) + cited = set(re.findall(r"\[(\w+)\]_", doc)) + return sorted(defined - cited), sorted(cited - defined) + + +def test_footnote_check_finds_problems(): + doc = "Text [1]_\n\nReferences\n----------\n.. [1] A\n.. [2] B\n" + assert footnote_problems(doc) == (["2"], []) + assert footnote_problems("A citation [3]_") == ([], ["3"]) + assert footnote_problems("References\n----------\nA\n\nB\n") == ([], []) + + +@pytest.mark.parametrize("name", list(DOCUMENTED)) +def test_footnotes_are_cited(name): + # Sphinx warns "Footnote [n] is not referenced" for each of them. Use + # author-year references, as other docstrings do, or cite the footnotes + doc = inspect.getdoc(DOCUMENTED[name]) + assert footnote_problems(doc) == ([], []) + + +SOURCES = { + "cragg_donald": [ + "Cragg, J. G., & Donald, S. G. (1993)", + "Anderson, T. W. (1951)", + "Stock, J. H., & Yogo, M. (2005)", + ], + "cragg_donald_f": ["Stock, J. H., & Yogo, M. (2005)"], + "kleibergen_paap": [ + "Kleibergen, F., & Paap, R. (2006)", + "Kleibergen, F., & Schaffer, M. E. (2007)", + ], + "kleibergen_paap_f": [ + "Kleibergen, F., & Paap, R. (2006)", + "Olea, J. L. M., & Pflueger, C. (2013)", + "Baum, C. F., Schaffer, M. E., & Stillman, S. (2007)", + ], +} + + +@pytest.mark.parametrize("name", list(SOURCES)) +def test_sources_are_listed(name): + # Removing the footnote markup keeps the references + doc = inspect.getdoc(DOCUMENTED[name]) + references = doc.split("References\n----------\n")[1] + for source in SOURCES[name]: + assert source in references + + # ---------------------------------------------------------------------------- # IVGMMResults.centered_j_stat # ---------------------------------------------------------------------------- From 361c443c62887b63dce7b3e0a5de3c2747e45a83 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:33:55 +0100 Subject: [PATCH 4/8] BUG: Do not compute the rk Wald covariance when only the LM statistic is requested kleibergen_paap computed the rk Wald statistic as well as the LM statistic, and discarded it. When the covariance of the Wald statistic was rank deficient, it warned that "the covariance of the Kleibergen-Paap statistic is rank deficient", about a statistic that was not requested. The scores of the Wald statistic are first-stage residuals, which sum to zero, and so their cluster sums have a covariance whose rank is one less than the number of clusters. The LM statistic needs ninstr - nendog + 1 clusters, and the Wald statistic needs one more. With exactly ninstr - nendog + 1 clusters the LM statistic is well defined, and kleibergen_paap and FirstStageResults.kleibergen_paap warned anyway. The rk Wald covariance is now only computed for the F statistic, by kleibergen_paap_f and FirstStageResults.kleibergen_paap_f, which warn if it is rank deficient as before. The values of all of the statistics are unchanged, and the LM statistic avoids about half of the work. The tests check, for several numbers of instruments and endogenous regressors, that the LM statistic does not warn with ninstr - nendog + 1 clusters and warns with one fewer, that the F statistic warns with ninstr - nendog + 1 clusters and not with one more, and the same through the results of a weighted model. They fail without the change. --- linearmodels/iv/common.py | 37 +++++++---- linearmodels/iv/results.py | 16 +++-- .../tests/iv/test_new_feature_hardening.py | 61 ++++++++++++++++++- 3 files changed, 94 insertions(+), 20 deletions(-) diff --git a/linearmodels/iv/common.py b/linearmodels/iv/common.py index 9af8d910c1..6c23750c61 100644 --- a/linearmodels/iv/common.py +++ b/linearmodels/iv/common.py @@ -423,8 +423,17 @@ def _kleibergen_paap( exog: linearmodels.typing.data.Float64Array, cov_type: str, cov_config: dict[str, Any] | None, + *, + wald: bool = True, ) -> tuple[WaldTestStatistic | InvalidTestStatistic, float]: - """Compute the rk LM test and the rk Wald F-statistic""" + """ + Compute the rk LM test and the rk Wald F-statistic + + The F-statistic is NaN if wald is False. The covariance of the scores of + the Wald statistic is not computed in that case. Its rank can be less than + the rank of the covariance of the LM statistic, and so it can warn when + the LM statistic is well-defined. + """ name = "Kleibergen-Paap rk LM Test" family = _cov_family(cov_type) config = {} if cov_config is None else dict(cov_config) @@ -463,17 +472,19 @@ def _kleibergen_paap( return singular, nan lm_chi2, lm_df = lm_result lm_test = WaldTestStatistic(lm_chi2, _RANK_NULL, lm_df, name=name) - - # The Wald version uses the first-stage residuals to form the scores - residuals = dec.x - dec.z @ dec.pihat - wald_result = _rk_chi2(dec, residuals, family, config) - dof = dec.nobs - dec.ninstr - dec.nexog - if family == "clustered": - nclusters = unique(config["clusters"]).shape[0] - scale = dof * (nclusters - 1) / nclusters / (dec.nobs - 1) - else: - scale = dof / dec.nobs - f = nan if wald_result is None else wald_result[0] * scale / dec.ninstr + f = nan + if wald: + # The Wald version uses the first-stage residuals to form the scores + residuals = dec.x - dec.z @ dec.pihat + wald_result = _rk_chi2(dec, residuals, family, config) + dof = dec.nobs - dec.ninstr - dec.nexog + if family == "clustered": + nclusters = unique(config["clusters"]).shape[0] + scale = dof * (nclusters - 1) / nclusters / (dec.nobs - 1) + else: + scale = dof / dec.nobs + if wald_result is not None: + f = wald_result[0] * scale / dec.ninstr return lm_test, f @@ -607,7 +618,7 @@ def kleibergen_paap( the rank of a matrix using the Kleibergen-Paap rk statistic. *Statistical Software Components* S456865, Boston College. """ - return _kleibergen_paap(endog, instr, exog, cov_type, cov_config)[0] + return _kleibergen_paap(endog, instr, exog, cov_type, cov_config, wald=False)[0] def kleibergen_paap_f( diff --git a/linearmodels/iv/results.py b/linearmodels/iv/results.py index 6ba2912d71..681c067cd3 100644 --- a/linearmodels/iv/results.py +++ b/linearmodels/iv/results.py @@ -759,12 +759,16 @@ def cragg_donald_f(self) -> float: return cragg_donald_f(*self._weighted_arrays()) @cached_property - def _kleibergen_paap( - self, - ) -> tuple[WaldTestStatistic | InvalidTestStatistic, float]: + def _rk_lm(self) -> WaldTestStatistic | InvalidTestStatistic: + return _kleibergen_paap( + *self._weighted_arrays(), self._cov_type, self._cov_config, wald=False + )[0] + + @cached_property + def _rk_wald_f(self) -> float: return _kleibergen_paap( *self._weighted_arrays(), self._cov_type, self._cov_config - ) + )[1] @property def kleibergen_paap(self) -> WaldTestStatistic | InvalidTestStatistic: @@ -809,7 +813,7 @@ def kleibergen_paap(self) -> WaldTestStatistic | InvalidTestStatistic: Stata's ``ivreg2``. See :func:`linearmodels.iv.common.kleibergen_paap` for the definition and for conventions that differ from Stata. """ - return self._kleibergen_paap[0] + return self._rk_lm @property def kleibergen_paap_f(self) -> float: @@ -840,7 +844,7 @@ def kleibergen_paap_f(self) -> float: "Kleibergen-Paap rk Wald F statistic". See :func:`linearmodels.iv.common.kleibergen_paap_f` for the definition. """ - return self._kleibergen_paap[1] + return self._rk_wald_f @cached_property def diagnostics(self) -> DataFrame: diff --git a/linearmodels/tests/iv/test_new_feature_hardening.py b/linearmodels/tests/iv/test_new_feature_hardening.py index 62898d336c..babf8c7889 100644 --- a/linearmodels/tests/iv/test_new_feature_hardening.py +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -11,13 +11,14 @@ import inspect import re +import warnings import numpy as np from numpy.testing import assert_allclose import pandas as pd import pytest -from linearmodels.iv import IVGMM, IVGMMCUE +from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE from linearmodels.iv.common import ( cragg_donald, cragg_donald_f, @@ -112,6 +113,64 @@ def test_sources_are_listed(name): assert source in references +def rank_design(seed, nobs=300, ninstr=4, nendog=2, weighted=False): + rng = np.random.default_rng(seed) + exog = np.column_stack([np.ones(nobs), rng.standard_normal(nobs)]) + z = rng.standard_normal((nobs, ninstr)) + pi = 0.4 * rng.standard_normal((ninstr, nendog)) + scale = np.exp(0.5 * z[:, [0]]) + x = z @ pi + 0.5 * exog[:, [1]] + scale * rng.standard_normal((nobs, nendog)) + w = rng.uniform(0.3, 3.0, nobs) if weighted else None + return x, z, exog, w + + +@pytest.mark.parametrize(("ninstr", "nendog"), [(4, 2), (5, 1), (6, 2), (5, 3)]) +def test_number_of_clusters_needed(ninstr, nendog): + # The covariance of the LM statistic has dimension ninstr - nendog + 1, and + # that of the Wald F statistic has scores with a mean of zero, so that it + # needs one more cluster. The LM statistic is computed without the Wald + # covariance, and so does not warn if it is full rank + x, z, exog, _ = rank_design(11, nobs=240, ninstr=ninstr, nendog=nendog) + df = ninstr - nendog + 1 + + def config(nclusters): + return {"clusters": np.arange(x.shape[0]) % nclusters} + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*rank deficient.*") + lm = kleibergen_paap(x, z, exog, "clustered", config(df)) + wald = kleibergen_paap_f(x, z, exog, "clustered", config(df + 1)) + assert lm.df == df + assert np.isfinite(lm.stat) + assert np.isfinite(wald) + with pytest.warns(UserWarning, match="rank deficient"): + short = kleibergen_paap_f(x, z, exog, "clustered", config(df)) + assert np.isfinite(short) + with pytest.warns(UserWarning, match="rank deficient"): + lm_short = kleibergen_paap(x, z, exog, "clustered", config(df - 1)) + assert lm_short.df == df - 1 + + +def test_number_of_clusters_needed_through_the_model(): + # The covariance of the model has to be full rank, which needs more clusters + y, exog, endog, z, w = simulate(13, nobs=240, ninstr=6, nendog=2, weighted=True) + df = 5 + + def first_stage(nclusters): + clusters = np.arange(y.shape[0]) % nclusters + mod = IV2SLS(y, exog, endog, z, weights=w) + return mod.fit(cov_type="clustered", clusters=clusters).first_stage + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*rank deficient.*") + assert first_stage(df).kleibergen_paap.df == df + assert np.isfinite(first_stage(df + 1).kleibergen_paap_f) + with pytest.warns(UserWarning, match="rank deficient"): + _ = first_stage(df).kleibergen_paap_f + with pytest.warns(UserWarning, match="rank deficient"): + _ = first_stage(df - 1).kleibergen_paap + + # ---------------------------------------------------------------------------- # IVGMMResults.centered_j_stat # ---------------------------------------------------------------------------- From 3150bf899f97c1c15cc67054b5ed7556fb83741a Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:33:57 +0100 Subject: [PATCH 5/8] DOC: Document the weights and clusters of the first-stage rank tests Document what the first-stage tests of the rank do with weights and clusters, in the docstrings of cragg_donald, cragg_donald_f, kleibergen_paap, kleibergen_paap_f and the corresponding properties of FirstStageResults. The behavior has not changed, and each statement is checked by a test. * Weights. The functions use the weighted data: each array has to be multiplied by the square root of the observation weights, including the constant column of exog, which is then the root of the weights and not a column of ones. The sample size is the number of rows. The statistics do not change if all of the weights are multiplied by a constant, and the results of the IV models use weights that are scaled to have mean 1. A test checks the scale of the weights for the functions, using a scale from 1e-6 to 1e6, and for the models, that the constant has to be weighted, and that the statistics of a weighted model are the ones of the functions. * Clusters. cragg_donald does not use a covariance estimator, and so does not depend on the clusters of the model, which is checked by comparing the four covariance estimators. kleibergen_paap sums the weighted scores within the clusters in cov_config, with one element for each row, and does not support two-way clusters. The LM statistic needs at least ninstr - nendog + 1 clusters and the Wald F statistic needs one more, and they warn if there are fewer. The previous commit makes the count correct for the LM statistic. * The first-stage tests of the rank of the results are the same for IV2SLS, IVLIML, IVGMM, IVGMMCUE and IVJIVE, with and without weights, and for the unadjusted, robust and clustered covariance estimators. The warnings of kleibergen_paap now also say that the covariance is formed as a dense matrix with (ninstr * nendog)**2 elements, so that the memory grows with the square of the number of instruments, which matters for the many-instrument models that JIVE is used for. The Raises section includes clusters with the wrong number of elements. The sources of kleibergen_paap_f are cited in the text. The new tests pass without this change, since they check the behavior. --- linearmodels/iv/common.py | 69 ++++++++++-- linearmodels/iv/results.py | 25 +++++ .../tests/iv/test_new_feature_hardening.py | 106 +++++++++++++++++- 3 files changed, 189 insertions(+), 11 deletions(-) diff --git a/linearmodels/iv/common.py b/linearmodels/iv/common.py index 6c23750c61..e329141d29 100644 --- a/linearmodels/iv/common.py +++ b/linearmodels/iv/common.py @@ -254,6 +254,21 @@ def cragg_donald( times the classical F-statistic when there is a single endogenous regressor. The F form, ``CD / k``, is reported by :func:`cragg_donald_f`. + **Weights.** In a weighted model, each of the arrays has to be multiplied + by the square root of the observation weights, including the constant + column of ``exog``. The constant is then a column of the square roots of + the weights and not a column of ones, since the first-stage regressions + are weighted. The sample size :math:`n` is the number of rows, and not the + sum of the weights. The statistic does not change if all of the weights are + multiplied by the same constant. The results of the IV models do this + using weights that have been scaled to have mean 1. + + **Clusters.** The statistic does not use a covariance estimator, and so + does not depend on clusters. It is the same for a model that is estimated + with a clustered, kernel or robust covariance estimator as it is for one + with an unadjusted covariance estimator. See :func:`kleibergen_paap` for a + statistic that does. + References ---------- Cragg, J. G., & Donald, S. G. (1993). Testing identifiability and @@ -323,6 +338,11 @@ def cragg_donald_f( critical values for. It is reported by Stata's ``ivreg2`` as the "Cragg-Donald Wald F statistic". + The weights are used as in :func:`cragg_donald`. Each array, including the + constant of ``exog``, is multiplied by the square root of the weights in a + weighted model, and the statistic does not depend on the scale of the + weights or on clusters. + References ---------- Stock, J. H., & Yogo, M. (2005). Testing for weak instruments in linear IV @@ -533,7 +553,8 @@ def kleibergen_paap( ------ ValueError If ``cov_type`` is not supported, or is "clustered" and - ``cov_config`` does not contain ``clusters``. + ``cov_config`` does not contain ``clusters`` or ``clusters`` does not + have an element for each row of the arrays. Warnings -------- @@ -547,14 +568,15 @@ def kleibergen_paap( instruments and first-stage errors are clustered or persistent. The covariance has dimension :math:`km \times km` and is poorly estimated - unless :math:`n` is large relative to :math:`km`. With clustering it needs - many clusters, and with few clusters the test can be very conservative. - A kernel bandwidth that is too short for the dependence does not remove - the over-rejection. If there are too few clusters or observations the - covariance is rank deficient. A warning is issued, a generalized inverse - is used and the degrees of freedom are reduced, but the result should not - be relied on. The chi-square approximation is also less reliable with - many instruments. + unless :math:`n` is large relative to :math:`km`. It is formed as a dense + matrix, and so the memory that is needed grows with :math:`(km)^2`. With + clustering it needs many clusters, and with few clusters the test can be + very conservative. A kernel bandwidth that is too short for the dependence + does not remove the over-rejection. If there are too few clusters or + observations the covariance is rank deficient. A warning is issued, a + generalized inverse is used and the degrees of freedom are reduced, but + the result should not be relied on. The chi-square approximation is also + less reliable with many instruments. The null is that the model is underidentified. Rejecting it does not show that the instruments are strong, so this is not a test for weak @@ -609,6 +631,23 @@ def kleibergen_paap( :math:`1 - j/(b+1)`, so a Stata bandwidth of :math:`b+1` corresponds to a ``bandwidth`` of :math:`b`. + **Weights.** In a weighted model, each of the arrays has to be multiplied + by the square root of the observation weights, including the constant + column of ``exog``, as for :func:`cragg_donald`. The scores are then those + of the weighted first-stage regressions, and the covariance estimator is + applied to them. The statistic does not change if all of the weights are + multiplied by the same constant. + + **Clusters.** With "clustered", the scores are summed within the clusters + in ``cov_config["clusters"]``. ``clusters`` has one element for each row + of the arrays, in the order of the rows. It is a 1-d array, or a 2-d array + with a single column. Two-way clustering is not supported. The LM + statistic is a quadratic form in a covariance with :math:`k-m+1` rows, + whose rank is at most the number of clusters, and so it needs at least + :math:`k-m+1` clusters. The scores of the Wald F statistic of + :func:`kleibergen_paap_f` sum to zero, and so it needs one more cluster. + A warning is issued if there are fewer clusters. + References ---------- Kleibergen, F., & Paap, R. (2006). Generalized reduced rank tests using the @@ -699,7 +738,17 @@ def kleibergen_paap_f( where :math:`M` is the number of clusters. With ``cov_type`` equal to "unadjusted" or "homoskedastic" the statistic is the Cragg-Donald F statistic from :func:`cragg_donald_f`. These are the definitions used by - Stata's ``ivreg2`` for the "Kleibergen-Paap rk Wald F statistic". + Stata's ``ivreg2`` for the "Kleibergen-Paap rk Wald F statistic" (Baum, + Schaffer and Stillman, 2007). + + The weights and clusters are used as in :func:`kleibergen_paap`. In a + weighted model, each of the arrays, including the constant of ``exog``, has + to be multiplied by the square root of the weights, and the statistic does + not depend on their scale. The clusters have one element for each row of + the arrays. If there are :math:`m` endogenous regressors, the scores of the + Wald statistic sum to zero, so that at least :math:`k-m+2` clusters are + needed, which is one more than for :func:`kleibergen_paap`. A warning is + issued if there are fewer clusters. References ---------- diff --git a/linearmodels/iv/results.py b/linearmodels/iv/results.py index 681c067cd3..4652f52a8a 100644 --- a/linearmodels/iv/results.py +++ b/linearmodels/iv/results.py @@ -728,6 +728,13 @@ def cragg_donald(self) -> WaldTestStatistic | InvalidTestStatistic: It equals ``ninstr`` times the classical F-statistic when there is one endogenous regressor. See :func:`linearmodels.iv.common.cragg_donald` for the definition. + + In a weighted model the statistic is computed from the data multiplied + by the square root of the weights of the model, which are scaled to + have mean 1. It does not depend on the scale of the weights. It does not + use a covariance estimator, and so it does not depend on the clusters + or the kernel of the model, and is the same for every covariance + estimator. """ return cragg_donald(*self._weighted_arrays()) @@ -755,6 +762,10 @@ def cragg_donald_f(self) -> float: This is the statistic that Stock and Yogo (2005) tabulate critical values for, and Stata's ``ivreg2`` reports it as the "Cragg-Donald Wald F statistic". See :func:`linearmodels.iv.common.cragg_donald_f`. + + The weights are used as for :attr:`cragg_donald`, and the statistic + does not depend on the scale of the weights or on the clusters of the + model. """ return cragg_donald_f(*self._weighted_arrays()) @@ -812,6 +823,15 @@ def kleibergen_paap(self) -> WaldTestStatistic | InvalidTestStatistic: corresponds to the "Kleibergen-Paap rk LM statistic" reported by Stata's ``ivreg2``. See :func:`linearmodels.iv.common.kleibergen_paap` for the definition and for conventions that differ from Stata. + + In a weighted model the scores are those of the weighted first-stage + regressions, and the statistic does not depend on the scale of the + weights. A clustered covariance estimator uses the ``clusters`` that + were passed to ``fit``. There must be at least ``ninstr - nendog + 1`` + clusters, or ``ninstr - nendog + 2`` for :attr:`kleibergen_paap_f`, and + a warning is issued if there are fewer. A kernel covariance estimator + uses the kernel and the bandwidth of the model, including a bandwidth + that the model selected. """ return self._rk_lm @@ -843,6 +863,11 @@ def kleibergen_paap_f(self) -> float: is the Cragg-Donald F statistic. Stata's ``ivreg2`` reports it as the "Kleibergen-Paap rk Wald F statistic". See :func:`linearmodels.iv.common.kleibergen_paap_f` for the definition. + + The weights, clusters and kernel are used as for + :attr:`kleibergen_paap`. With a clustered covariance estimator there + must be at least ``ninstr - nendog + 2`` clusters, and a warning is + issued if there are fewer. """ return self._rk_wald_f diff --git a/linearmodels/tests/iv/test_new_feature_hardening.py b/linearmodels/tests/iv/test_new_feature_hardening.py index babf8c7889..2a50dceedd 100644 --- a/linearmodels/tests/iv/test_new_feature_hardening.py +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -18,7 +18,7 @@ import pandas as pd import pytest -from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE +from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE, IVJIVE, IVLIML from linearmodels.iv.common import ( cragg_donald, cragg_donald_f, @@ -124,6 +124,91 @@ def rank_design(seed, nobs=300, ninstr=4, nendog=2, weighted=False): return x, z, exog, w +COV_CASES = [ + ("unadjusted", None), + ("robust", None), + ("clustered", {"clusters": np.arange(300) // 10}), + ("kernel", {"bandwidth": 4}), +] + + +def all_rank_statistics(endog, instr, exog): + out = [cragg_donald(endog, instr, exog).stat, cragg_donald_f(endog, instr, exog)] + for cov_type, config in COV_CASES: + out.append(kleibergen_paap(endog, instr, exog, cov_type, config).stat) + out.append(kleibergen_paap_f(endog, instr, exog, cov_type, config)) + return np.array(out) + + +def test_rank_statistics_do_not_depend_on_the_scale_of_the_weights(): + # The arrays are the data multiplied by the square root of the weights + x, z, exog, w = rank_design(1, weighted=True) + root = np.sqrt(w)[:, None] + base = all_rank_statistics(root * x, root * z, root * exog) + assert np.isfinite(base).all() + for scale in (1e-6, 1 / 37.5, 37.5, 1e6): + factor = np.sqrt(scale) + arrays = [factor * root * a for a in (x, z, exog)] + assert_allclose(all_rank_statistics(*arrays), base, rtol=1e-8) + + +def test_first_stage_does_not_depend_on_the_scale_of_the_weights(): + y, exog, endog, z, w = simulate(2, weighted=True) + options = {"cov_type": "clustered", "clusters": np.arange(y.shape[0]) // 10} + base = IV2SLS(y, exog, endog, z, weights=w).fit(**options).first_stage + scaled = IV2SLS(y, exog, endog, z, weights=37.5 * w).fit(**options).first_stage + ones = IV2SLS(y, exog, endog, z, weights=np.ones(y.shape[0])).fit(**options) + plain = IV2SLS(y, exog, endog, z).fit(**options) + for name in ("kleibergen_paap_f", "cragg_donald_f"): + assert_allclose(getattr(scaled, name), getattr(base, name), rtol=1e-10) + assert_allclose( + getattr(ones.first_stage, name), getattr(plain.first_stage, name) + ) + assert_allclose(scaled.kleibergen_paap.stat, base.kleibergen_paap.stat) + # The weights matter + assert not np.isclose(base.kleibergen_paap_f, plain.first_stage.kleibergen_paap_f) + + +def test_the_constant_of_a_weighted_model_is_the_root_of_the_weights(): + # The functions use the weighted arrays, including the constant + y, exog, endog, z, w = simulate(3, weighted=True) + first_stage = ( + IV2SLS(y, exog, endog, z, weights=w).fit(cov_type="robust").first_stage + ) + root = np.sqrt(w / w.mean())[:, None] + stat = kleibergen_paap(root * endog, root * z, root * exog, "robust").stat + assert_allclose(stat, first_stage.kleibergen_paap.stat, rtol=1e-10) + assert_allclose( + cragg_donald_f(root * endog, root * z, root * exog), + first_stage.cragg_donald_f, + rtol=1e-10, + ) + # A column of ones is not the constant of the weighted model + unit_constant = exog.copy() + unit_constant[:, 0] = 1.0 + wrong = kleibergen_paap(root * endog, root * z, unit_constant, "robust").stat + assert not np.isclose(wrong, first_stage.kleibergen_paap.stat, rtol=1e-3) + + +@pytest.mark.parametrize("estimator", [IV2SLS, IVLIML, IVGMM, IVGMMCUE, IVJIVE]) +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize("cov_type", ["unadjusted", "robust", "clustered"]) +def test_first_stage_rank_tests_do_not_depend_on_the_estimator( + estimator, cov_type, weighted +): + y, exog, endog, z, w = simulate(12, ninstr=4, nendog=2, weighted=weighted) + options = {} + if cov_type == "clustered": + options["clusters"] = np.arange(y.shape[0]) // 10 + reference = IV2SLS(y, exog, endog, z, weights=w).fit(cov_type=cov_type, **options) + res = estimator(y, exog, endog, z, weights=w).fit(cov_type=cov_type, **options) + expected, actual = reference.first_stage, res.first_stage + assert_allclose(actual.kleibergen_paap.stat, expected.kleibergen_paap.stat) + assert_allclose(actual.kleibergen_paap_f, expected.kleibergen_paap_f) + assert_allclose(actual.cragg_donald.stat, expected.cragg_donald.stat) + assert_allclose(actual.cragg_donald_f, expected.cragg_donald_f) + + @pytest.mark.parametrize(("ninstr", "nendog"), [(4, 2), (5, 1), (6, 2), (5, 3)]) def test_number_of_clusters_needed(ninstr, nendog): # The covariance of the LM statistic has dimension ninstr - nendog + 1, and @@ -171,6 +256,25 @@ def first_stage(nclusters): _ = first_stage(df - 1).kleibergen_paap +def test_cragg_donald_ignores_the_covariance_estimator_and_clusters(): + y, exog, endog, z, w = simulate(15, weighted=True) + mod = IV2SLS(y, exog, endog, z, weights=w) + clusters = np.arange(y.shape[0]) // 10 + fits = [ + mod.fit(cov_type="unadjusted"), + mod.fit(cov_type="robust"), + mod.fit(cov_type="clustered", clusters=clusters), + mod.fit(cov_type="kernel", bandwidth=4), + ] + cragg = [fit.first_stage.cragg_donald.stat for fit in fits] + cragg_f = [fit.first_stage.cragg_donald_f for fit in fits] + assert_allclose(cragg, cragg[0]) + assert_allclose(cragg_f, cragg_f[0]) + # and the Kleibergen-Paap statistic uses it + kp = [fit.first_stage.kleibergen_paap.stat for fit in fits] + assert len({round(value, 6) for value in kp}) == len(kp) + + # ---------------------------------------------------------------------------- # IVGMMResults.centered_j_stat # ---------------------------------------------------------------------------- From b71f8445a0ae4110a24cb2ebc4a838fa21d1875a Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:34:06 +0100 Subject: [PATCH 6/8] DOC: Fix the IVJIVE.from_formula example and add IVJIVE to model.__all__ The example in the docstring of IVJIVE.from_formula is np.log(wage) ~ 1 + exper + exper ** 2 + brthord + [educ ~ sibs] In a formula, a ** 2 is the power of a set of terms, which is a again, and not its square. The model that the example creates has no squared term. Its parameters are Intercept, exper, brthord and educ. The example now uses I(exper ** 2), and the notes say why. They also say that the formula expressions have to be valid and that a constant has to be included with "1 +", as the other from_formula docstrings do. The examples of IVLIML, IV2SLS, IVGMM and IVGMMCUE have the same formula, and are not changed here. IVJIVE is also not in __all__ of linearmodels.iv.model, which lists the other estimators, and is now. The tests run the examples in the docstrings of IVJIVE and IVJIVE.from_formula, check that the formula in the example has the squared term, as a parameter of the model, and that IVJIVE is in __all__. --- linearmodels/iv/model.py | 10 +++- .../tests/iv/test_new_feature_hardening.py | 51 +++++++++++++++++++ 2 files changed, 60 insertions(+), 1 deletion(-) diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 4116cf54c4..945e44e528 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -68,6 +68,7 @@ "IV2SLS", "IVGMM", "IVGMMCUE", + "IVJIVE", "IVLIML", "WEIGHT_MATRICES", "_OLS", @@ -2056,13 +2057,20 @@ def from_formula( block of the form [endog ~ instruments] which is used to indicate the list of endogenous variables and instruments. The general structure of an IV formula is `dependent ~ exog [endog ~ instruments]` + and it must be the case that the formula expressions constructed from + blocks `dependent ~ exog endog` and `dependent ~ exog instruments` are + both valid formulas. + + A constant must be explicitly included using "1 +" if required. Use + `I(exper ** 2)` for the square of a variable, since `exper ** 2` is + the power of a set of terms in a formula, which is `exper`. Examples -------- >>> from linearmodels.datasets import wage >>> from linearmodels.iv import IVJIVE >>> data = wage.load() - >>> formula = "np.log(wage) ~ 1 + exper + exper ** 2 + brthord + [educ ~ sibs]" + >>> formula = "np.log(wage) ~ 1 + exper + I(exper ** 2) + brthord + [educ ~ sibs]" >>> mod = IVJIVE.from_formula(formula, data) """ parser = IVFormulaParser(formula, data) diff --git a/linearmodels/tests/iv/test_new_feature_hardening.py b/linearmodels/tests/iv/test_new_feature_hardening.py index 2a50dceedd..9823cfd4b3 100644 --- a/linearmodels/tests/iv/test_new_feature_hardening.py +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -1,6 +1,9 @@ """ Checks of the instrumental variable features that were added after version 7.0 +* IVJIVE: the docstring examples, formulas with weights and missing values, + prediction, comparison with the other estimators, the summary and the + first-stage tests of the rank. * The tests of the rank of the first-stage coefficient matrix, which are cragg_donald, cragg_donald_f, kleibergen_paap and kleibergen_paap_f: that the documentation is well-formed, and the behavior with weights and clusters that @@ -9,6 +12,7 @@ GMM, weights, clusters, kernels and singular covariances. """ +import doctest import inspect import re import warnings @@ -18,6 +22,7 @@ import pandas as pd import pytest +from linearmodels.datasets import wage from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE, IVJIVE, IVLIML from linearmodels.iv.common import ( cragg_donald, @@ -25,7 +30,9 @@ kleibergen_paap, kleibergen_paap_f, ) +import linearmodels.iv.model from linearmodels.iv.results import FirstStageResults, IVGMMResults +from linearmodels.shared.exceptions import MissingValueWarning from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic @@ -47,6 +54,50 @@ def simulate(seed, nobs=300, ninstr=4, nendog=1, violation=0.0, weighted=False): return y, exog, endog, z, w +# ---------------------------------------------------------------------------- +# IVJIVE +# ---------------------------------------------------------------------------- + + +def run_docstring_examples(obj, name): + """Run the examples in the docstrings of an object and of its members""" + finder = doctest.DocTestFinder(exclude_empty=True) + failures = [] + with_examples = set() + for test in finder.find(obj, name, globs={}): + runner = doctest.DocTestRunner(optionflags=doctest.ELLIPSIS) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", MissingValueWarning) + runner.run(test, out=failures.append) + if test.examples: + with_examples.add(test.name) + return with_examples, failures + + +def test_jive_docstring_examples_run(): + with_examples, failures = run_docstring_examples(IVJIVE, "IVJIVE") + assert not failures, "".join(failures) + assert {"IVJIVE", "IVJIVE.from_formula"} <= with_examples + + +def test_jive_from_formula_example_has_every_term(): + # exper ** 2 is the power of a set of terms in a formula, which is exper, and + # not its square. Without I() the example silently omits the term that it + # names + doc = IVJIVE.from_formula.__doc__ + formula = re.search(r'formula = "(.*)"', doc).group(1) + assert "exper ** 2" not in formula.replace("I(exper ** 2)", "") + with warnings.catch_warnings(): + warnings.simplefilter("ignore", MissingValueWarning) + mod = IVJIVE.from_formula(formula, wage.load()) + assert "I(exper ** 2)" in list(mod.fit().params.index) + + +def test_jive_is_in_the_model_module_all(): + assert "IVJIVE" in linearmodels.iv.model.__all__ + assert linearmodels.iv.model.IVJIVE is IVJIVE + + # ---------------------------------------------------------------------------- # Tests of the rank of the first-stage coefficient matrix # ---------------------------------------------------------------------------- From e5aa78e64dadd6d501f22e8eac30978e000916b1 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:34:09 +0100 Subject: [PATCH 7/8] TST: Test IVJIVE with formulas, weights, prediction, comparison and pickling Add tests of the parts of IVJIVE and of the first-stage tests of the rank that the tests of the features did not check. They pass without any change to the code, since they found no defect in the behavior that they test. * from_formula with weights, when an exogenous variable has missing values, is the model of the data without those rows, with the weights of those rows also dropped, and the weights are used. Only the exogenous variables of a formula can have missing values. * predict with data, which does not need the dependent variable or the instruments, with exog and endog, and in sample, for a model with weights. * compare of IVJIVE with IV2SLS, IVLIML and IVGMM, with weights, from a dict or a list, with each precision and stars, and for models from formulas. * The summary with weights and the small sample adjustment, for each of the covariance estimators. * The first-stage tests of the rank of IVJIVE equal those of IV2SLS for each covariance estimator with and without weights, and the first stage, including its summary, works for each estimator and covariance estimator, with weights. * The results of IV2SLS, IVGMM, IVGMMCUE and IVJIVE can be pickled, including the first stage and the centered J statistic. --- .../tests/iv/test_new_feature_hardening.py | 202 +++++++++++++++++- 1 file changed, 201 insertions(+), 1 deletion(-) diff --git a/linearmodels/tests/iv/test_new_feature_hardening.py b/linearmodels/tests/iv/test_new_feature_hardening.py index 9823cfd4b3..d6c5bd9672 100644 --- a/linearmodels/tests/iv/test_new_feature_hardening.py +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -14,6 +14,7 @@ import doctest import inspect +import pickle import re import warnings @@ -23,7 +24,7 @@ import pytest from linearmodels.datasets import wage -from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE, IVJIVE, IVLIML +from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE, IVJIVE, IVLIML, compare from linearmodels.iv.common import ( cragg_donald, cragg_donald_f, @@ -54,6 +55,27 @@ def simulate(seed, nobs=300, ninstr=4, nendog=1, violation=0.0, weighted=False): return y, exog, endog, z, w +def simulated_frame(seed, nobs=200, weighted=True): + """The same type of data in a DataFrame, with names for a formula""" + y, exog, endog, z, w = simulate(seed, nobs=nobs, ninstr=3, weighted=weighted) + data = pd.DataFrame( + { + "y": y, + "x1": exog[:, 1], + "x2": endog[:, 0], + "z1": z[:, 0], + "z2": z[:, 1], + "z3": z[:, 2], + } + ) + if weighted: + data["w"] = w + return data + + +FORMULA = "y ~ 1 + x1 + [x2 ~ z1 + z2 + z3]" + + # ---------------------------------------------------------------------------- # IVJIVE # ---------------------------------------------------------------------------- @@ -98,6 +120,184 @@ def test_jive_is_in_the_model_module_all(): assert linearmodels.iv.model.IVJIVE is IVJIVE +def test_jive_from_formula_weights_and_missing_values(): + # Only the exogenous regressors of a formula can contain missing values + data = simulated_frame(1) + data.loc[[4, 77, 10], "x1"] = np.nan + with pytest.warns(MissingValueWarning, match="missing values"): + mod = IVJIVE.from_formula(FORMULA, data, weights=data.w) + res = mod.fit() + keep = data.notnull().all(axis=1) + clean = data[keep] + expected = IVJIVE.from_formula(FORMULA, clean, weights=clean.w).fit() + unweighted = IVJIVE.from_formula(FORMULA, clean).fit() + assert res.nobs == keep.sum() == 197 + assert list(res.params.index) == ["Intercept", "x1", "x2"] + assert_allclose(res.params, expected.params, rtol=1e-10) + assert_allclose(res.cov, expected.cov, rtol=1e-10) + assert_allclose(res.resids, expected.resids, rtol=1e-10) + # The weights are used + assert not np.allclose(res.params, unweighted.params, rtol=1e-3) + assert mod.formula == FORMULA + + +def test_jive_predict_with_formula_data_and_weights(): + data = simulated_frame(2) + res = IVJIVE.from_formula(FORMULA, data, weights=data.w).fit(cov_type="robust") + new = data.iloc[:7] + x_new = np.column_stack([np.ones(7), new.x1, new.x2]) + expected = x_new @ res.params.to_numpy() + pred = res.predict(data=new) + assert list(pred.columns) == ["predictions"] + assert pred.index.equals(new.index) + assert_allclose(pred.to_numpy().squeeze(), expected, rtol=1e-12) + # The dependent variable and the instruments are not needed + only_regressors = res.predict(data=new[["x1", "x2"]]) + assert_allclose(only_regressors.to_numpy().squeeze(), expected, rtol=1e-12) + # In sample predictions are the fitted values, whatever the weights + in_sample = res.predict(data=data) + assert_allclose(in_sample.to_numpy(), res.fitted_values.to_numpy(), rtol=1e-10) + assert_allclose(res.predict().to_numpy(), res.fitted_values.to_numpy()) + # Values of the regressors are an alternative to data + exog_new = new[["x1"]].assign(const=1.0)[["const", "x1"]] + arrays = res.predict(exog_new, new[["x2"]]) + assert_allclose(arrays.to_numpy().squeeze(), expected, rtol=1e-12) + with pytest.raises(ValueError, match="but not both"): + res.predict(new[["x1"]], new[["x2"]], data=new) + + +def test_jive_compare_with_other_estimators_and_weights(): + y, exog, endog, z, w = simulate(3, weighted=True) + results = { + "JIVE": IVJIVE(y, exog, endog, z, weights=w).fit(), + "2SLS": IV2SLS(y, exog, endog, z, weights=w).fit(), + "LIML": IVLIML(y, exog, endog, z, weights=w).fit(), + "GMM": IVGMM(y, exog, endog, z, weights=w).fit(), + } + comparison = compare(results) + assert list(comparison.params.columns) == list(results) + for name, res in results.items(): + assert_allclose(comparison.params[name].to_numpy(), res.params.to_numpy()) + assert_allclose(comparison.tstats[name].to_numpy(), res.tstats.to_numpy()) + assert_allclose(comparison.rsquared[name], res.rsquared) + assert comparison.nobs[name] == res.nobs + methods = comparison.estimator_method + assert methods.to_dict() == { + "JIVE": "IV-JIVE", + "2SLS": "IV-2SLS", + "LIML": "IV-LIML", + "GMM": "IV-GMM", + } + assert set(comparison.cov_estimator) == {"robust"} + text = str(comparison.summary) + for label in ("JIVE", "2SLS", "LIML", "GMM", "IV-JIVE", "IV-GMM"): + assert label in text + # A list of results, stars and the other precisions + as_list = compare(list(results.values()), precision="std_errors", stars=True) + assert_allclose(as_list.std_errors.iloc[:, 0], results["JIVE"].std_errors) + assert "***" in str(as_list.summary) + + +def test_jive_compare_formula_models(): + data = simulated_frame(4) + results = { + "JIVE": IVJIVE.from_formula(FORMULA, data, weights=data.w).fit(), + "2SLS": IV2SLS.from_formula(FORMULA, data, weights=data.w).fit(), + } + comparison = compare(results) + assert list(comparison.params.index) == ["Intercept", "x1", "x2"] + assert_allclose(comparison.params["JIVE"], results["JIVE"].params) + + +@pytest.mark.parametrize("cov_type", ["unadjusted", "robust", "kernel", "clustered"]) +def test_jive_summary_with_weights_and_small_sample_adjustment(cov_type): + y, exog, endog, z, w = simulate(5, weighted=True) + options = {} + if cov_type == "kernel": + options["bandwidth"] = 3 + elif cov_type == "clustered": + options["clusters"] = np.arange(y.shape[0]) // 12 + res = IVJIVE(y, exog, endog, z, weights=w).fit( + cov_type=cov_type, debiased=True, **options + ) + text = str(res.summary) + assert "IV-JIVE" in text + assert "Debiased: True" in text + assert "JIVE Covariance" in text + # The model F-statistic of a debiased model is an F + assert res.f_statistic.dist_name.startswith("F(") + assert res.f_statistic.dist_name in text + assert "Endogenous: endog" in text + assert "Instruments: instruments.0" in text + assert max(len(line) for line in text.split("\n")) <= 78 + + +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize( + ("cov_type", "options"), + [ + ("unadjusted", {}), + ("robust", {}), + ("clustered", {"clusters": np.arange(300) // 10}), + ("kernel", {"kernel": "bartlett", "bandwidth": 4}), + ], +) +def test_jive_first_stage_rank_tests_match_2sls(cov_type, options, weighted): + # The first stage does not depend on the estimator, only on the covariance + # estimator, the weights and the clusters of the model + y, exog, endog, z, w = simulate(10, ninstr=5, nendog=2, weighted=weighted) + jive = IVJIVE(y, exog, endog, z, weights=w).fit(cov_type=cov_type, **options) + tsls = IV2SLS(y, exog, endog, z, weights=w).fit(cov_type=cov_type, **options) + stats_jive, stats_tsls = jive.first_stage, tsls.first_stage + assert type(stats_jive.kleibergen_paap) is WaldTestStatistic + assert stats_jive.kleibergen_paap.df == 4 + for name in ("kleibergen_paap_f", "cragg_donald_f"): + assert_allclose(getattr(stats_jive, name), getattr(stats_tsls, name)) + for name in ("kleibergen_paap", "cragg_donald"): + assert_allclose(getattr(stats_jive, name).stat, getattr(stats_tsls, name).stat) + + +@pytest.mark.parametrize("estimator", [IV2SLS, IVLIML, IVGMM, IVGMMCUE, IVJIVE]) +@pytest.mark.parametrize("cov_type", ["unadjusted", "robust", "kernel", "clustered"]) +def test_first_stage_of_every_estimator_and_covariance(estimator, cov_type): + # The covariance configuration of the results is used by the first stage, + # and differs by the estimator + y, exog, endog, z, w = simulate(30, nobs=200, ninstr=3, weighted=True) + options = {} + if cov_type == "clustered": + options["clusters"] = np.arange(200) // 10 + elif cov_type == "kernel": + options["bandwidth"] = 3 + res = estimator(y, exog, endog, z, weights=w).fit(cov_type=cov_type, **options) + first_stage = res.first_stage + assert "First Stage Estimation Results" in str(first_stage) + assert first_stage.diagnostics.shape == (1, 6) + assert type(first_stage.kleibergen_paap) is WaldTestStatistic + assert first_stage.kleibergen_paap.df == 3 + assert np.isfinite(first_stage.kleibergen_paap_f) + assert np.isfinite(first_stage.cragg_donald_f) + + +@pytest.mark.parametrize("estimator", [IV2SLS, IVGMM, IVGMMCUE, IVJIVE]) +def test_results_can_be_pickled(estimator): + y, exog, endog, z, w = simulate(31, nobs=200, ninstr=3, weighted=True) + clusters = np.arange(200) // 10 + res = estimator(y, exog, endog, z, weights=w).fit( + cov_type="clustered", clusters=clusters + ) + clone = pickle.loads(pickle.dumps(res)) + assert_allclose(clone.params, res.params) + assert_allclose(clone.cov, res.cov) + assert str(clone.summary) == str(res.summary) + first_stage, original = clone.first_stage, res.first_stage + assert_allclose(first_stage.kleibergen_paap.stat, original.kleibergen_paap.stat) + assert_allclose(first_stage.kleibergen_paap_f, original.kleibergen_paap_f) + assert_allclose(first_stage.cragg_donald_f, original.cragg_donald_f) + if isinstance(res, IVGMMResults): + assert_allclose(clone.centered_j_stat.stat, res.centered_j_stat.stat) + assert str(clone.centered_j_stat) == str(res.centered_j_stat) + + # ---------------------------------------------------------------------------- # Tests of the rank of the first-stage coefficient matrix # ---------------------------------------------------------------------------- From aed232b52879730b228f367455bd668e392f306a Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 20:34:09 +0100 Subject: [PATCH 8/8] TST: Test NormalTestStatistic and document its tails NormalTestStatistic, which pesaran_cd returns, had tests of a p-value, the critical values and the str and repr of a two-sided statistic. Add tests of the rest, which pass without a change to the code, since they found no defect: * The critical values are the quantiles of the normal distribution for tests of size 10%, 5% and 1%, two-sided and one-sided, and the p-value at a critical value is the size of the test. The one-sided critical values are smaller. The two-sided p-value is the p-value of a chi-squared statistic with 1 degree of freedom, and is symmetric. The one-sided p-values of a statistic and its negative add to 1. * The str and repr of a one-sided statistic, with and without a name, of a small p-value, and the p-values and critical values for infinite and NaN statistics, and for statistics that are NumPy scalars, 0-d arrays or ints. * NormalTestStatistic, WaldTestStatistic with and without a denominator degree of freedom, InvalidTestStatistic and InapplicableTestStatistic can be pickled with each protocol, copied and deep copied, and give the same strings, p-values and critical values, and the frozen distribution still works. This is needed to send results to another process. The docstring said that p-values are two-sided if two_sided is True, and did not say what the one-sided p-value and critical values are. They are for a test that rejects when the statistic is large: the p-value is the upper tail probability, and the critical values are the 0.9, 0.95 and 0.99 quantiles, which means that the p-value of a large negative statistic is close to 1. The str of a one-sided statistic does not say that it is one-sided, as before. --- linearmodels/shared/hypotheses.py | 20 +- .../shared/test_normal_test_statistic.py | 195 ++++++++++++++++++ 2 files changed, 211 insertions(+), 4 deletions(-) create mode 100644 linearmodels/tests/shared/test_normal_test_statistic.py diff --git a/linearmodels/shared/hypotheses.py b/linearmodels/shared/hypotheses.py index c65b0ca449..d1df946744 100644 --- a/linearmodels/shared/hypotheses.py +++ b/linearmodels/shared/hypotheses.py @@ -100,7 +100,7 @@ def __repr__(self) -> str: class NormalTestStatistic: - """ + r""" Test statistic holder for tests with an asymptotic standard normal law. Parameters @@ -109,11 +109,23 @@ class NormalTestStatistic: The test statistic. null : str A statement of the test's null hypothesis. - name : str + name : str, optional Name of test. - two_sided : bool + two_sided : bool, optional Flag indicating whether p-values and critical values are computed - using a two-sided test. + using a two-sided test. The default is True. + + Notes + ----- + A two-sided test rejects the null if the absolute value of the statistic + :math:`s` is large. The p-value is :math:`2\Pr(Z > |s|)`, where :math:`Z` + has a standard normal distribution, and the critical values for tests of + size 10%, 5% and 1% are the 0.95, 0.975 and 0.995 quantiles of the standard + normal. + + A one-sided test rejects the null if the statistic is large. The p-value is + :math:`\Pr(Z > s)`, which is close to 1 if the statistic is large and + negative, and the critical values are the 0.9, 0.95 and 0.99 quantiles. """ def __init__( diff --git a/linearmodels/tests/shared/test_normal_test_statistic.py b/linearmodels/tests/shared/test_normal_test_statistic.py new file mode 100644 index 0000000000..1fa347c47b --- /dev/null +++ b/linearmodels/tests/shared/test_normal_test_statistic.py @@ -0,0 +1,195 @@ +""" +Tests of NormalTestStatistic and of the other holders of test statistics + +The p-values and the critical values are checked against each other and against +the normal distribution and the chi-squared distribution with one degree of +freedom, and the objects are checked to survive being pickled and copied, which +is needed to send results to other processes. +""" + +import copy +import pickle + +import numpy as np +from numpy.testing import assert_allclose +import pytest +from scipy import stats + +from linearmodels.shared.hypotheses import ( + InapplicableTestStatistic, + InvalidTestStatistic, + NormalTestStatistic, + WaldTestStatistic, +) + +SIZES = {"10%": 0.10, "5%": 0.05, "1%": 0.01} +PROTOCOLS = range(pickle.HIGHEST_PROTOCOL + 1) + + +@pytest.mark.parametrize("two_sided", [True, False]) +def test_critical_values_are_quantiles_of_the_normal(two_sided): + crit = NormalTestStatistic(0.0, "_NULL_", two_sided=two_sided).critical_values + assert list(crit) == list(SIZES) + # The test rejects when the statistic exceeds the critical value in + # absolute value if two-sided and when it is large if one-sided + tail = np.array(list(SIZES.values())) / (2 if two_sided else 1) + assert_allclose(list(crit.values()), stats.norm.isf(tail), rtol=1e-12) + assert_allclose(list(crit.values()), stats.norm.ppf(1 - tail), rtol=1e-12) + # Smaller sizes have larger critical values + assert list(crit.values()) == sorted(crit.values()) + + +@pytest.mark.parametrize("two_sided", [True, False]) +@pytest.mark.parametrize("label", list(SIZES)) +def test_p_value_at_the_critical_value_is_the_size(label, two_sided): + crit = NormalTestStatistic(0.0, "_NULL_", two_sided=two_sided).critical_values + at_critical_value = NormalTestStatistic(crit[label], "_NULL_", two_sided=two_sided) + assert_allclose(at_critical_value.pval, SIZES[label], rtol=1e-10) + # Statistics beyond the critical value reject + for sign in (1, -1) if two_sided else (1,): + beyond = NormalTestStatistic( + sign * (crit[label] + 0.01), "_NULL_", two_sided=two_sided + ) + assert beyond.pval < SIZES[label] + + +def test_one_sided_critical_values_are_smaller_than_two_sided(): + one = NormalTestStatistic(0.0, "_NULL_", two_sided=False).critical_values + two = NormalTestStatistic(0.0, "_NULL_", two_sided=True).critical_values + for label in SIZES: + assert one[label] < two[label] + assert_allclose(one["10%"], stats.norm.isf(0.1)) + assert_allclose(two["5%"], stats.norm.isf(0.025)) + + +@pytest.mark.parametrize("stat", [-4.0, -1.3, 0.0, 0.2, 1.3, 4.0]) +def test_two_sided_p_value_is_a_chi2_p_value_and_symmetric(stat): + ts = NormalTestStatistic(stat, "_NULL_") + wald = WaldTestStatistic(stat**2, "_NULL_", 1) + assert_allclose(ts.pval, wald.pval, rtol=1e-10) + assert_allclose(NormalTestStatistic(-stat, "_NULL_").pval, ts.pval, rtol=1e-12) + assert 0 < ts.pval <= 1 + + +@pytest.mark.parametrize("stat", [-4.0, -1.3, 0.0, 0.2, 1.3, 4.0]) +def test_one_sided_p_value_is_for_a_large_statistic(stat): + ts = NormalTestStatistic(stat, "_NULL_", two_sided=False) + assert_allclose(ts.pval, stats.norm.sf(stat), rtol=1e-12) + flipped = NormalTestStatistic(-stat, "_NULL_", two_sided=False) + assert_allclose(ts.pval + flipped.pval, 1.0, rtol=1e-12) + larger = NormalTestStatistic(stat + 0.5, "_NULL_", two_sided=False) + assert larger.pval < ts.pval + + +def test_nonfinite_statistics(): + ts = NormalTestStatistic(np.nan, "_NULL_") + assert np.isnan(ts.pval) + assert "Statistic: nan" in str(ts) + assert "P-value: nan" in str(ts) + # The critical values do not depend on the statistic + assert_allclose( + list(ts.critical_values.values()), + list(NormalTestStatistic(1.0, "_NULL_").critical_values.values()), + ) + assert NormalTestStatistic(np.inf, "_NULL_").pval == 0.0 + assert NormalTestStatistic(-np.inf, "_NULL_").pval == 0.0 + assert NormalTestStatistic(np.inf, "_NULL_", two_sided=False).pval == 0.0 + assert NormalTestStatistic(-np.inf, "_NULL_", two_sided=False).pval == 1.0 + + +@pytest.mark.parametrize( + "stat", + [np.float64(2.0), np.array(2.0), 2, 2.0], + ids=["numpy-scalar", "0-d-array", "int", "float"], +) +def test_numpy_and_integer_statistics(stat): + ts = NormalTestStatistic(stat, "_NULL_", name="_NAME_") + assert_allclose(ts.pval, 2 * stats.norm.sf(2.0), rtol=1e-12) + assert "Statistic: 2.0000" in str(ts) + + +def test_str_and_repr_of_one_sided_statistic(): + ts = NormalTestStatistic(1.5, "_NULL_", name="_NAME_", two_sided=False) + expected = ( + "_NAME_\nH0: _NULL_\nStatistic: 1.5000\nP-value: 0.0668\nDistributed: N(0,1)" + ) + assert str(ts) == expected + assert repr(ts) == expected + "\nNormalTestStatistic, id: " + hex(id(ts)) + unnamed = NormalTestStatistic(-1.5, "_NULL_", two_sided=False) + assert str(unnamed) == ( + "H0: _NULL_\nStatistic: -1.5000\nP-value: 0.9332\nDistributed: N(0,1)" + ) + + +def test_str_and_repr_of_two_sided_statistic_with_a_small_p_value(): + ts = NormalTestStatistic(-9.0, "_NULL_", name="_NAME_") + assert "P-value: 0.0000\n" in str(ts) + assert str(ts).startswith("_NAME_\nH0: _NULL_\nStatistic: -9.0000\n") + assert repr(ts).startswith(str(ts)) + assert repr(ts).endswith(f"NormalTestStatistic, id: {hex(id(ts))}") + + +CASES = { + "normal": NormalTestStatistic(1.5, "_NULL_", name="_NAME_"), + "normal-one-sided": NormalTestStatistic( + -1.5, "_NULL_", name="_NAME_", two_sided=False + ), + "normal-nan": NormalTestStatistic(np.nan, "_NULL_"), + "wald-chi2": WaldTestStatistic(3.0, "_NULL_", 2, name="_NAME_"), + "wald-f": WaldTestStatistic(3.0, "_NULL_", 2, 40), + "invalid": InvalidTestStatistic("_REASON_", name="_NAME_"), + "inapplicable": InapplicableTestStatistic(reason="_REASON_", name="_NAME_"), + "inapplicable-default": InapplicableTestStatistic(), +} + + +def assert_same_statistic(clone, original): + assert type(clone) is type(original) + assert clone.dist_name == original.dist_name + if hasattr(original, "df"): + assert clone.df == original.df + assert clone.df_denom == original.df_denom + assert str(clone) == str(original) + assert_allclose(clone.stat, original.stat) + assert_allclose(clone.pval, original.pval) + assert clone.null == original.null + if original.critical_values is None: + assert clone.critical_values is None + else: + assert clone.critical_values.keys() == original.critical_values.keys() + assert_allclose( + list(clone.critical_values.values()), + list(original.critical_values.values()), + ) + # The distribution that is used to compute the p-value is usable + assert_allclose( + clone.dist.sf(1.25), original.dist.sf(1.25), rtol=1e-12, equal_nan=True + ) + + +@pytest.mark.parametrize("protocol", PROTOCOLS) +@pytest.mark.parametrize("name", list(CASES)) +def test_pickle_round_trip(name, protocol): + original = CASES[name] + clone = pickle.loads(pickle.dumps(original, protocol=protocol)) + assert clone is not original + assert_same_statistic(clone, original) + # The hex id in the repr is the id of the object + assert repr(clone).endswith(hex(id(clone))) + + +@pytest.mark.parametrize("name", list(CASES)) +def test_copy_and_deepcopy(name): + original = CASES[name] + assert_same_statistic(copy.copy(original), original) + assert_same_statistic(copy.deepcopy(original), original) + + +def test_pickle_preserves_the_side_of_the_test(): + two_sided = pickle.loads(pickle.dumps(NormalTestStatistic(1.5, "_NULL_"))) + one_sided = pickle.loads( + pickle.dumps(NormalTestStatistic(1.5, "_NULL_", two_sided=False)) + ) + assert_allclose(two_sided.pval, 2 * stats.norm.sf(1.5)) + assert_allclose(one_sided.pval, stats.norm.sf(1.5)) + assert one_sided.critical_values["5%"] < two_sided.critical_values["5%"]