BUG: Fix issues in IV

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]
 

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.