diff --git a/linearmodels/iv/common.py b/linearmodels/iv/common.py index 6317fae137..e329141d29 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,23 +247,41 @@ 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 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 ---------- - .. [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) @@ -319,11 +338,16 @@ 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 ---------- - .. [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): @@ -419,8 +443,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) @@ -459,17 +492,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 @@ -518,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 -------- @@ -532,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 @@ -594,16 +631,33 @@ 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 ---------- - .. [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] + return _kleibergen_paap(endog, instr, exog, cov_type, cov_config, wald=False)[0] def kleibergen_paap_f( @@ -684,19 +738,29 @@ 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 ---------- - .. [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/iv/model.py b/linearmodels/iv/model.py index 5643020c30..945e44e528 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, ) @@ -36,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, @@ -64,6 +68,7 @@ "IV2SLS", "IVGMM", "IVGMMCUE", + "IVJIVE", "IVLIML", "WEIGHT_MATRICES", "_OLS", @@ -1070,7 +1075,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 = ( @@ -1879,6 +1894,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"): @@ -2027,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) @@ -2079,9 +2116,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/iv/results.py b/linearmodels/iv/results.py index c02e1cbcb8..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,16 +762,24 @@ 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()) @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: @@ -808,8 +823,17 @@ 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._kleibergen_paap[0] + return self._rk_lm @property def kleibergen_paap_f(self) -> float: @@ -839,8 +863,13 @@ 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._kleibergen_paap[1] + return self._rk_wald_f @cached_property def diagnostics(self) -> DataFrame: @@ -1656,7 +1685,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 +1736,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/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/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(): 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..d6c5bd9672 --- /dev/null +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -0,0 +1,744 @@ +""" +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 + it describes. +* IVGMMResults.centered_j_stat: the units of the data, continuously updating + GMM, weights, clusters, kernels and singular covariances. +""" + +import doctest +import inspect +import pickle +import re +import warnings + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.datasets import wage +from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE, IVJIVE, IVLIML, compare +from linearmodels.iv.common import ( + cragg_donald, + cragg_donald_f, + 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 + + +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 + + +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 +# ---------------------------------------------------------------------------- + + +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 + + +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 +# ---------------------------------------------------------------------------- + +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 + + +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 + + +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 + # 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 + + +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 +# ---------------------------------------------------------------------------- + + +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) 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%"]