diff --git a/linearmodels/asset_pricing/model.py b/linearmodels/asset_pricing/model.py index 743094ebf0..fd4f2bfddb 100644 --- a/linearmodels/asset_pricing/model.py +++ b/linearmodels/asset_pricing/model.py @@ -25,6 +25,7 @@ LinearFactorModelResults, ) from linearmodels.iv.data import IVData, IVDataLike +from linearmodels.shared.alignment import check_alignment, row_labels from linearmodels.shared.exceptions import missing_warning from linearmodels.shared.hypotheses import WaldTestStatistic from linearmodels.shared.linalg import has_constant @@ -82,6 +83,11 @@ class _FactorModelBase: def __init__(self, portfolios: IVDataLike, factors: IVDataLike): self.portfolios = IVData(portfolios, var_name="portfolio") self.factors = IVData(factors, var_name="factor") + check_alignment( + row_labels(self.portfolios), + {"factors": self.factors}, + reference="portfolios", + ) self._name = self.__class__.__name__ self._formula: str | None = None self._validate_data() @@ -214,6 +220,18 @@ class TradedFactorModel(_FactorModelBase): tested using the estimated values :math:`\hat{\alpha}_i`. Risk premia, :math:`\lambda_i` are estimated using the sample averages of the factors, which must be excess returns on traded portfolios. + + Rows are paired by position, so that the first row of ``factors`` is the + factor return in the period of the first row of ``portfolios``. If both + ``portfolios`` and ``factors`` carry labels (a pandas Series or DataFrame + with an index other than the default 0, 1, ..., n-1, or an xarray + DataArray with coordinates), ``factors`` must use the same index, in the + same order, as ``portfolios``, so that a return cannot be attached to the + wrong period. A ``ValueError`` is raised if the index differs. Either + reindex ``factors`` to the index of ``portfolios``, or call + ``to_numpy()`` to pair its values by position. NumPy arrays and pandas + objects with a default index do not carry labels and are always paired by + position. """ def __init__(self, portfolios: IVDataLike, factors: IVDataLike): @@ -523,6 +541,18 @@ class LinearFactorModel(_LinearFactorModelBase): The model is tested using the estimated values :math:`\hat{\alpha}_i=\hat{\eta}_i`. + + Rows are paired by position, so that the first row of ``factors`` is the + factor return in the period of the first row of ``portfolios``. If both + ``portfolios`` and ``factors`` carry labels (a pandas Series or DataFrame + with an index other than the default 0, 1, ..., n-1, or an xarray + DataArray with coordinates), ``factors`` must use the same index, in the + same order, as ``portfolios``, so that a return cannot be attached to the + wrong period. A ``ValueError`` is raised if the index differs. Either + reindex ``factors`` to the index of ``portfolios``, or call + ``to_numpy()`` to pair its values by position. NumPy arrays and pandas + objects with a default index do not carry labels and are always paired by + position. """ def __init__( @@ -832,6 +862,18 @@ class LinearFactorModelGMM(_LinearFactorModelBase): The model is tested using the optimized objective function using the usual GMM J statistic. + + Rows are paired by position, so that the first row of ``factors`` is the + factor return in the period of the first row of ``portfolios``. If both + ``portfolios`` and ``factors`` carry labels (a pandas Series or DataFrame + with an index other than the default 0, 1, ..., n-1, or an xarray + DataArray with coordinates), ``factors`` must use the same index, in the + same order, as ``portfolios``, so that a return cannot be attached to the + wrong period. A ``ValueError`` is raised if the index differs. Either + reindex ``factors`` to the index of ``portfolios``, or call + ``to_numpy()`` to pair its values by position. NumPy arrays and pandas + objects with a default index do not carry labels and are always paired by + position. """ def __init__( diff --git a/linearmodels/iv/absorbing.py b/linearmodels/iv/absorbing.py index 9f001ea952..f55500ce2a 100644 --- a/linearmodels/iv/absorbing.py +++ b/linearmodels/iv/absorbing.py @@ -55,6 +55,7 @@ not_absorbed, preconditioner, ) +from linearmodels.shared.alignment import check_alignment, is_default_index, row_labels from linearmodels.shared.exceptions import missing_warning from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic from linearmodels.shared.utility import DataFrameWrapper, SeriesWrapper @@ -304,6 +305,11 @@ class Interaction: For each variable in `cont`, computes the interaction of the variable and the cartesian product of the categories. + If both ``cat`` and ``cont`` carry labels (a pandas Series or DataFrame + with an index other than the default 0, 1, ..., n-1, or an xarray + DataArray with coordinates), they must use the same index, in the same + order. + Examples -------- >>> import numpy as np @@ -369,6 +375,9 @@ def _check_data(self) -> None: self._cont_data = IVData(cont, "cont", nobs=nobs, convert_dummies=False) if self._cat_data.shape[1] == self._cont_data.shape[1] == 0: raise ValueError("Both cat and cont are empty arrays") + check_alignment( + row_labels(self._cat_data), {"cont": self._cont_data}, reference="cat" + ) cat_data = self._cat_data.pandas convert = [ col for col in cat_data if not (isinstance(cat_data[col], CategoricalDtype)) @@ -391,7 +400,11 @@ def cont(self) -> DataFrame: @property def isnull(self) -> Series: - return self.cat.isnull().any(axis=1) | self.cont.isnull().any(axis=1) + cat_null = self.cat.isnull().any(axis=1) + cont_null = self.cont.isnull().any(axis=1) + # Rows are paired by position since cat and cont need not have labels + index = cont_null.index if is_default_index(cat_null.index) else cat_null.index + return Series(cat_null.to_numpy() | cont_null.to_numpy(), index=index) def drop(self, locs: linearmodels.typing.data.BoolArray) -> None: self._cat_data.drop(locs) @@ -660,6 +673,18 @@ class AbsorbingLS: estimates a particular restricted set of parameters that captures the effect of non-redundant components in z. + Rows are paired by position. An input that carries labels (a pandas + Series or DataFrame with an index other than the default 0, 1, ..., n-1, + or an xarray DataArray with coordinates) must use the same index, in the + same order, as ``dependent``, so that a value cannot be attached to the + wrong observation. A ``ValueError`` is raised if the index of ``exog``, + ``absorb``, ``interactions`` or ``weights`` differs from the index of + ``dependent``. Either reindex the input to the index of ``dependent``, or + call ``to_numpy()`` to pair its values by position. NumPy arrays and + pandas objects with a default index do not carry labels and are always + paired by position. The same rule applies to ``clusters`` passed to + ``fit``. + See also -------- Interaction @@ -725,6 +750,7 @@ def __init__( self._absorbed_exog: DataFrame | None = None self._check_shape() + self._check_alignment() self._original_index = self._dependent.pandas.index self._drop_locs = self._drop_missing() self._columns = self._exog.cols @@ -786,6 +812,19 @@ def _check_shape(self) -> None: "observations".format(str(interact)) ) + def _check_alignment(self) -> None: + """Check that inputs with labels use the labels of dependent""" + expected = row_labels(self._dependent) + check_alignment(expected, {"exog": self._exog, "weights": self._weights}) + interactions: list[tuple[str, Interaction]] = [("absorb", self._absorb_inter)] + single = isinstance(self._interactions, (DataFrame, Interaction)) + for i, interact in enumerate(self._interaction_list): + name = "interactions" if single else f"interactions[{i}]" + interactions.append((name, interact)) + for name, interact in interactions: + check_alignment(expected, {name: interact.cat}) + check_alignment(expected, {name: interact.cont}) + @property def absorbed_dependent(self) -> DataFrame: """ @@ -1094,6 +1133,9 @@ def fit( stacklevel=2, ) absorb_options = dict(lsmr_options.items()) + check_alignment( + row_labels(self._dependent), {"clusters": cov_config.get("clusters")} + ) if self._absorbed_dependent is None: self._first_time_fit(use_cache, absorb_options, method) diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 5643020c30..4023db7c7c 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -4,7 +4,7 @@ from __future__ import annotations -from collections.abc import Callable +from collections.abc import Callable, Mapping from typing import Any, TypeVar, Union, cast import warnings @@ -50,6 +50,7 @@ OneWayClusteredWeightMatrix, ) from linearmodels.iv.results import IVGMMResults, IVResults, OLSResults +from linearmodels.shared.alignment import check_alignment, row_labels from linearmodels.shared.exceptions import IndexWarning, missing_warning from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic from linearmodels.shared.linalg import has_constant, inv_sqrth @@ -181,7 +182,17 @@ def __init__( self._original_index = self.dependent.pandas.index if weights is None: weights = ones(self.dependent.shape) - weights = IVData(weights).ndarray + weights_data = IVData(weights) + check_alignment( + row_labels(self.dependent), + { + "exog": self.exog, + "endog": self.endog, + "instruments": self.instruments, + "weights": weights_data, + }, + ) + weights = weights_data.ndarray if npany(weights <= 0): raise ValueError("weights must be strictly positive.") weights = weights / nanmean(weights) @@ -379,6 +390,12 @@ def _drop_missing(self) -> linearmodels.typing.data.BoolArray: missing_warning(missing, stacklevel=4) return missing + def _check_clusters(self, config: Mapping[str, Any]) -> None: + """Check that clusters with labels use the labels of dependent""" + check_alignment( + row_labels(self.dependent), {"clusters": config.get("clusters")} + ) + def wresids( self, params: linearmodels.typing.data.Float64Array ) -> linearmodels.typing.data.Float64Array: @@ -677,6 +694,7 @@ def fit( linearmodels.iv.covariance.KernelCovariance linearmodels.iv.covariance.ClusteredCovariance """ + self._check_clusters(cov_config) wy, wx, wz = self._wy, self._wx, self._wz kappa = self._kappa @@ -766,6 +784,18 @@ class IVLIML(_IVLSModelBase): \kappa-\alpha/(n-n_{instr}) + Rows are paired by position. An input that carries labels (a pandas + Series or DataFrame with an index other than the default 0, 1, ..., n-1, + or an xarray DataArray with coordinates) must use the same index, in the + same order, as ``dependent``, so that a value cannot be attached to the + wrong observation. A ``ValueError`` is raised if the index of ``exog``, + ``endog``, ``instruments`` or ``weights`` differs from the index of + ``dependent``. Either reindex the input to the index of ``dependent``, or + call ``to_numpy()`` to pair its values by position. NumPy arrays and + pandas objects with a default index do not carry labels and are always + paired by position. The same rule applies to ``clusters`` passed to + ``fit``. + .. todo:: * VCV: bootstrap @@ -887,6 +917,18 @@ class IV2SLS(_IVLSModelBase): The 2SLS estimator is a special case of a k-class estimator with :math:`\kappa=1`, + Rows are paired by position. An input that carries labels (a pandas + Series or DataFrame with an index other than the default 0, 1, ..., n-1, + or an xarray DataArray with coordinates) must use the same index, in the + same order, as ``dependent``, so that a value cannot be attached to the + wrong observation. A ``ValueError`` is raised if the index of ``exog``, + ``endog``, ``instruments`` or ``weights`` differs from the index of + ``dependent``. Either reindex the input to the index of ``dependent``, or + call ``to_numpy()`` to pair its values by position. NumPy arrays and + pandas objects with a default index do not carry labels and are always + paired by position. The same rule applies to ``clusters`` passed to + ``fit``. + .. todo:: * VCV: bootstrap @@ -1023,6 +1065,7 @@ def __init__( super().__init__(dependent, exog, endog, instruments, weights=weights) self._method = "IV-GMM" + self._check_clusters(weight_config) weight_matrix_estimator = WEIGHT_MATRICES[weight_type] self._weight = weight_matrix_estimator(**weight_config) self._weight_type = weight_type @@ -1140,6 +1183,18 @@ class IVGMM(_IVGMMBase): where :math:`W` is a positive definite weight matrix and :math:`Z` contains both the exogenous regressors and the instruments. + Rows are paired by position. An input that carries labels (a pandas + Series or DataFrame with an index other than the default 0, 1, ..., n-1, + or an xarray DataArray with coordinates) must use the same index, in the + same order, as ``dependent``, so that a value cannot be attached to the + wrong observation. A ``ValueError`` is raised if the index of ``exog``, + ``endog``, ``instruments`` or ``weights`` differs from the index of + ``dependent``. Either reindex the input to the index of ``dependent``, or + call ``to_numpy()`` to pair its values by position. NumPy arrays and + pandas objects with a default index do not carry labels and are always + paired by position. The same rule applies to ``clusters`` passed to + ``fit`` or in ``weight_config``. + .. todo:: * VCV: bootstrap @@ -1160,16 +1215,15 @@ def __init__( weight_type: str = "robust", **weight_config: Any, ): - super().__init__(dependent, exog, endog, instruments, weights=weights) - self._method = "IV-GMM" - - weight_matrix_estimator = WEIGHT_MATRICES[weight_type] - self._weight = weight_matrix_estimator(**weight_config) - self._weight_type = weight_type - self._weight_config = self._weight.config - # The arguments needed to create another weight estimator of this - # type. config does not contain all of them, e.g., optimal_bw - self._weight_kwargs = dict(weight_config) + super().__init__( + dependent, + exog, + endog, + instruments, + weights=weights, + weight_type=weight_type, + **weight_config, + ) @staticmethod def from_formula( @@ -1319,6 +1373,7 @@ def fit( -------- linearmodels.iv.gmm.IVGMMCovariance """ + self._check_clusters(cov_config) wy, wx, wz = self._wy, self._wx, self._wz nobs = wy.shape[0] weight_matrix = self._weight.weight_matrix @@ -1423,6 +1478,18 @@ class IVGMMCUE(_IVGMMBase): where :math:`W(\beta)` is a weight matrix that depends on :math:`\beta` through :math:`\epsilon_i = y_i - x_i\beta`. + Rows are paired by position. An input that carries labels (a pandas + Series or DataFrame with an index other than the default 0, 1, ..., n-1, + or an xarray DataArray with coordinates) must use the same index, in the + same order, as ``dependent``, so that a value cannot be attached to the + wrong observation. A ``ValueError`` is raised if the index of ``exog``, + ``endog``, ``instruments`` or ``weights`` differs from the index of + ``dependent``. Either reindex the input to the index of ``dependent``, or + call ``to_numpy()`` to pair its values by position. NumPy arrays and + pandas objects with a default index do not carry labels and are always + paired by position. The same rule applies to ``clusters`` passed to + ``fit`` or in ``weight_config``. + See Also -------- IV2SLS, IVLIML, IVGMM @@ -1686,6 +1753,7 @@ def fit( -------- linearmodels.iv.gmm.IVGMMCovariance """ + self._check_clusters(cov_config) wy, wx, wz = self._wy, self._wx, self._wz weight_matrix = self._weight.weight_matrix if starting is None: @@ -1745,6 +1813,17 @@ class _OLS(IVLIML): Uses IV2SLS internally by setting endog and instruments to None. Uses IVLIML with kappa=0 to estimate OLS models. + Rows are paired by position. An input that carries labels (a pandas + Series or DataFrame with an index other than the default 0, 1, ..., n-1, + or an xarray DataArray with coordinates) must use the same index, in the + same order, as ``dependent``, so that a value cannot be attached to the + wrong observation. A ``ValueError`` is raised if the index of ``exog`` or + ``weights`` differs from the index of ``dependent``. Either reindex the + input to the index of ``dependent``, or call ``to_numpy()`` to pair its + values by position. NumPy arrays and pandas objects with a default index + do not carry labels and are always paired by position. The same rule + applies to ``clusters`` passed to ``fit``. + See Also -------- statsmodels.regression.linear_model.OLS, @@ -1959,6 +2038,18 @@ class IVJIVE(_IVModelBase): variance, particularly when the instruments are weak (see Davidson and MacKinnon, 2006). + Rows are paired by position. An input that carries labels (a pandas + Series or DataFrame with an index other than the default 0, 1, ..., n-1, + or an xarray DataArray with coordinates) must use the same index, in the + same order, as ``dependent``, so that a value cannot be attached to the + wrong observation. A ``ValueError`` is raised if the index of ``exog``, + ``endog``, ``instruments`` or ``weights`` differs from the index of + ``dependent``. Either reindex the input to the index of ``dependent``, or + call ``to_numpy()`` to pair its values by position. NumPy arrays and + pandas objects with a default index do not carry labels and are always + paired by position. The same rule applies to ``clusters`` passed to + ``fit``. + See Also -------- IV2SLS, IVLIML, IVGMM, IVGMMCUE @@ -2090,6 +2181,7 @@ def fit( -------- linearmodels.iv.gmm.IVGMMCovariance """ + self._check_clusters(cov_config) wy, wx, wz = self._wy, self._wx, self._wz # The leverages and the first-stage predictions from an orthonormal diff --git a/linearmodels/panel/model.py b/linearmodels/panel/model.py index 10c8f8d37a..6a6ec6db67 100644 --- a/linearmodels/panel/model.py +++ b/linearmodels/panel/model.py @@ -11,15 +11,7 @@ from formulaic.parser.algos.tokenize import tokenize from formulaic.utils.context import capture_context import numpy as np -from pandas import ( - Categorical, - DataFrame, - Index, - MultiIndex, - RangeIndex, - Series, - get_dummies, -) +from pandas import Categorical, DataFrame, Index, MultiIndex, Series, get_dummies from scipy.linalg import lstsq as sp_lstsq from scipy.sparse.linalg import lsmr @@ -48,6 +40,13 @@ in_2core_graph, not_absorbed, ) +from linearmodels.shared.alignment import ( + check_row_alignment, + differs, + is_labelled, + labels_1d, + misaligned, +) from linearmodels.shared.exceptions import ( InferenceUnavailableWarning, MemoryWarning, @@ -297,72 +296,14 @@ class AmbiguityError(Exception): with a 2-level MultiIndex containing the entity- and time-ids.""" -def _is_default_index(index: Index) -> bool: - """Whether an index only numbers the rows, i.e., it is 0, 1, ..., n-1""" - return not isinstance(index, MultiIndex) and index.equals(RangeIndex(len(index))) - - def _is_labelled(value: object) -> bool: - """Whether the rows of an input are identified by labels. - - pandas, xarray and PanelData inputs carry labels. Raw arrays do not, and - are paired with the other inputs by position. The same is true of a pandas - Series or DataFrame that has a default index (and default columns if it - does not have a MultiIndex) since there are no labels to compare. - """ - if isinstance(value, Series): - return not _is_default_index(value.index) - if isinstance(value, DataFrame): - if isinstance(value.index, MultiIndex): - return True - return not (_is_default_index(value.index) and _is_default_index(value.columns)) - if isinstance(value, PanelData): - return True - try: - from xarray import DataArray # noqa: PLC0415 - except ImportError: - return False - return isinstance(value, DataArray) - + """Whether the rows of a panel input are identified by labels. -def _labels_1d(value: Series | DataFrame | Any) -> Index: - """Row labels of a Series, DataFrame or one-dimensional DataArray""" - if isinstance(value, (Series, DataFrame)): - return value.index - return value.get_index(value.dims[0]) - - -def _differs(labels: Index, expected: Index) -> bool: - """Whether labels, other than a default index, are not the expected labels""" - return not _is_default_index(labels) and not labels.equals(expected) - - -def _misaligned(name: str, what: str) -> ValueError: - return ValueError( - f"{what} of {name} does not match dependent. Rows are paired by " - "position, so an input with labels (pandas or xarray) must use the " - "same labels in the same order as dependent. Use a NumPy array to pair " - "the values by position." - ) - - -def _check_row_alignment(name: str, labels: Index, expected: Index) -> None: - """Reject labelled inputs whose row index is not the dependent index. - - Estimation pairs rows by position. An input with the same values in a - different order would attach those values to the wrong observation. - - Inputs with a different number of rows are not checked here since these - are reported using a length-specific error message. A default index has - no labels and is paired by position. + The columns of a DataFrame that has a default index also carry labels, since + they identify the entities when the frame is time by entity. See + :func:`linearmodels.shared.alignment.is_labelled`. """ - if ( - len(labels) != len(expected) - or _is_default_index(labels) - or labels.equals(expected) - ): - return - raise _misaligned(name, "The row index") + return is_labelled(value, columns=True) class _PanelModelBase: @@ -403,7 +344,7 @@ def __init__( self.exog = PanelData(exog, "Exog") self._dependent_labelled = _is_labelled(dependent) if self._dependent_labelled and _is_labelled(exog): - _check_row_alignment("exog", self.exog.index, self.dependent.index) + check_row_alignment("exog", self.exog.index, self.dependent.index) self._original_shape = self.dependent.shape self._constant = False self._formula: str | None = None @@ -466,7 +407,7 @@ def reformat_clusters( "and time periods as the model data." ) if self._dependent_labelled and _is_labelled(clusters): - _check_row_alignment("clusters", clusters_pd.index, self._original_index) + check_row_alignment("clusters", clusters_pd.index, self._original_index) clusters_pd.drop(~self.not_null) return clusters_pd.copy() @@ -503,7 +444,7 @@ def _adapt_weights(self, weights: PanelDataLike | None) -> PanelData: panel = PanelData(weights) if isinstance(weights, PanelData): if labelled: - _check_row_alignment("weights", panel.index, index) + check_row_alignment("weights", panel.index, index) return panel return self._grid_to_rows(panel, labelled) @@ -521,23 +462,23 @@ def _adapt_weights(self, weights: PanelDataLike | None) -> PanelData: ): panel = PanelData(DataFrame(weights)) if labelled: - _check_row_alignment("weights", panel.index, index) + check_row_alignment("weights", panel.index, index) return panel elif weights.shape[0] == nobs and weights.ndim == 1: - if labelled and _differs(_labels_1d(weights), times): - raise _misaligned("weights", "The time index") + if labelled and differs(labels_1d(weights), times): + raise misaligned("weights", "The time index") weights_arr = np.asarray(weights)[:, None] weights_arr = weights_arr @ np.ones((1, nentity)) grid = DataFrame(weights_arr, index=times, columns=entities) elif weights.shape[0] == nentity and weights.ndim == 1: - if labelled and _differs(_labels_1d(weights), entities): - raise _misaligned("weights", "The entity index") + if labelled and differs(labels_1d(weights), entities): + raise misaligned("weights", "The entity index") weights_arr = np.asarray(weights)[None, :] weights_arr = np.ones((nobs, 1)) @ weights_arr grid = DataFrame(weights_arr, index=times, columns=entities) elif weights.shape[0] == nentity * nobs and weights.size: if labelled: - _check_row_alignment("weights", _labels_1d(weights), index) + check_row_alignment("weights", labels_1d(weights), index) frame = self.dependent.dataframe.copy() # TODO: Fix this for pandas 3 frame.iloc[:, :] = np.asarray(weights)[:, None] @@ -575,7 +516,7 @@ def _grid_to_rows(self, grid: PanelData, labelled: bool) -> PanelData: if len(rows) != len(grid_index) or len(index) != len(grid_index): return grid if labelled and not rows.index.equals(grid_index): - raise _misaligned("weights", "The time and entity index") + raise misaligned("weights", "The time and entity index") rows = rows.set_axis(grid_index, axis=0) if not grid_index.equals(index): rows = rows.reindex(index) @@ -1491,7 +1432,7 @@ def _validate_effects( "entities and time periods as dependent." ) if labelled: - _check_row_alignment("other_effects", effects.index, self._original_index) + check_row_alignment("other_effects", effects.index, self._original_index) num_effects = effects.nvar if num_effects + self.entity_effects + self.time_effects > 2: diff --git a/linearmodels/shared/alignment.py b/linearmodels/shared/alignment.py new file mode 100644 index 0000000000..14f5f2a477 --- /dev/null +++ b/linearmodels/shared/alignment.py @@ -0,0 +1,261 @@ +""" +Pairing the rows of inputs that carry labels + +Estimators pair the rows of their inputs by position: the first row of the +regressors belongs to the first row of the dependent variable, and so on. An +input that carries labels (a pandas Series or DataFrame with an index, an +xarray DataArray, or the data containers built from these) must use the same +labels, in the same order, as the dependent variable. Otherwise the values +would be silently attached to the wrong observations. + +Inputs without labels, such as NumPy arrays and pandas objects that only have +a default index (0, 1, ..., n-1), are paired by position and only need to have +the right shape. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from typing import Any + +from pandas import DataFrame, Index, MultiIndex, RangeIndex, Series + + +def is_default_index(index: Index) -> bool: + """ + Whether an index only numbers the rows, i.e., it is 0, 1, ..., n-1. + + Parameters + ---------- + index : Index + The index to check. + + Returns + ------- + bool + True if the index is a default index, which does not have labels. + """ + return not isinstance(index, MultiIndex) and index.equals(RangeIndex(len(index))) + + +def is_labelled(value: object, *, columns: bool = False) -> bool: + """ + Whether the rows of an input are identified by labels. + + pandas, xarray, PanelData and IVData inputs can carry labels. Raw arrays do + not, and are paired with the other inputs by position. The same is true of + a pandas Series or DataFrame that has a default index, since there are no + labels to compare. + + Parameters + ---------- + value : object + The input to check. + columns : bool, optional + Whether the columns of a DataFrame also carry labels, so that a + DataFrame with a default index and non-default columns is labelled. + This is the case when the columns identify the entities of a panel. The + default is False, which is appropriate when the columns are variables. + + Returns + ------- + bool + True if the rows are identified by labels. + """ + from linearmodels.iv.data import IVData # noqa: PLC0415 + from linearmodels.panel.data import PanelData # noqa: PLC0415 + + if isinstance(value, Series): + return not is_default_index(value.index) + if isinstance(value, DataFrame): + return not is_default_index(value.index) or ( + columns and not is_default_index(value.columns) + ) + if isinstance(value, IVData): + # IVData stores the object it was built from + return is_labelled(value.original) + return isinstance(value, PanelData) or _is_data_array(value) + + +def _is_data_array(value: object) -> bool: + """Whether value is an xarray DataArray, if xarray is installed.""" + try: + from xarray import DataArray # noqa: PLC0415 + except ImportError: + return False + return isinstance(value, DataArray) + + +def labels_1d(value: Any) -> Index: + """ + Row labels of a Series, DataFrame, DataArray, IVData or PanelData. + + Parameters + ---------- + value : Series, DataFrame, DataArray, IVData or PanelData + The input. A DataArray must have named dimensions, and the labels are + taken from the first dimension. + + Returns + ------- + Index + The labels of the rows. + """ + from linearmodels.iv.data import IVData # noqa: PLC0415 + from linearmodels.panel.data import PanelData # noqa: PLC0415 + + if isinstance(value, (Series, DataFrame)): + return value.index + if isinstance(value, IVData): + return value.pandas.index + if isinstance(value, PanelData): + return value.index + return value.get_index(value.dims[0]) + + +def row_labels(value: object) -> Index | None: + """ + Row labels of an input, if it carries labels. + + Parameters + ---------- + value : object + The input, or None. + + Returns + ------- + Index or None + The labels of the rows, or None if the rows are not identified by + labels. Only the index of the rows is considered, so a DataFrame + with a default index is not labelled even if its columns are named. + """ + if not is_labelled(value): + return None + labels = labels_1d(value) + return None if is_default_index(labels) else labels + + +def differs(labels: Index, expected: Index) -> bool: + """ + Whether labels, other than a default index, are not the expected labels. + + Parameters + ---------- + labels : Index + The labels to check. + expected : Index + The expected labels. + + Returns + ------- + bool + True if labels is not a default index and is not equal to expected. + """ + return not is_default_index(labels) and not labels.equals(expected) + + +def misaligned(name: str, what: str, reference: str = "dependent") -> ValueError: + """ + Error to raise when the labels of an input are not those of the reference. + + Parameters + ---------- + name : str + Description of the input that does not match, which is used in the + message. + what : str + The part of the labels that does not match, e.g., "The row index". + reference : str, optional + Description of the input that sets the labels. The default is + "dependent". + + Returns + ------- + ValueError + The error, which explains how to resolve the mismatch. + """ + return ValueError( + f"{what} of {name} does not match {reference}. Rows are paired by " + "position, so an input with labels (pandas or xarray) must use the " + f"same labels in the same order as {reference}. Either reindex {name} " + f"to the index of {reference}, or call .to_numpy() on it to pair the " + "values by position." + ) + + +def check_row_alignment( + name: str, labels: Index, expected: Index, reference: str = "dependent" +) -> None: + """ + Reject labelled inputs whose row index is not the index of the reference. + + Estimation pairs rows by position. An input with the same values in a + different order would attach those values to the wrong observation. + + Inputs with a different number of rows are not checked here since these + are reported using a length-specific error message. A default index has + no labels and is paired by position. + + Parameters + ---------- + name : str + Description of the input, which is used in the message. + labels : Index + Labels of the rows of the input. + expected : Index + Labels of the rows of the reference. + reference : str, optional + Description of the input that sets the labels. The default is + "dependent". + + Raises + ------ + ValueError + If the labels differ from the expected labels. + """ + if ( + len(labels) != len(expected) + or is_default_index(labels) + or labels.equals(expected) + ): + return + raise misaligned(name, "The row index", reference) + + +def check_alignment( + expected: Index | None, + inputs: Mapping[str, object], + *, + reference: str = "dependent", +) -> None: + """ + Reject labelled inputs whose row index is not the expected index. + + Parameters + ---------- + expected : Index or None + Labels of the rows of the reference, usually dependent. If None, the + reference does not carry labels and so there is nothing to compare to: + all inputs are paired by position. + inputs : dict of str to object + The inputs to check, keyed by the name to use in the message. Inputs + that are None or do not carry labels are not checked. + reference : str, optional + Description of the input that sets the labels. The default is + "dependent". + + Raises + ------ + ValueError + If the labels of an input differ from the expected labels. + + See Also + -------- + row_labels : Row labels of an input, if it carries labels. + """ + if expected is None: + return + for name, value in inputs.items(): + labels = row_labels(value) + if labels is not None: + check_row_alignment(name, labels, expected, reference) diff --git a/linearmodels/system/model.py b/linearmodels/system/model.py index e836c99dab..ad5f6a01f7 100644 --- a/linearmodels/system/model.py +++ b/linearmodels/system/model.py @@ -27,6 +27,11 @@ from linearmodels.iv._utility import IVFormulaParser from linearmodels.iv.data import IVData +from linearmodels.shared.alignment import ( + check_alignment, + check_row_alignment, + row_labels, +) from linearmodels.shared.exceptions import missing_warning from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic from linearmodels.shared.linalg import has_constant @@ -485,6 +490,7 @@ def _validate_data(self) -> None: "{eq}".format(eq=self._eq_labels[i]) ) + self._check_alignment() self._drop_missing() self._common_exog = len(set(ids)) == 1 if self._common_exog: @@ -559,6 +565,49 @@ def _validate_data(self) -> None: ) self._constant_loc = constant_loc + def _check_alignment(self) -> None: + """Check that inputs with labels use the labels of dependent""" + nobs = [dep.shape[0] for dep in self._dependent] + for i, label in enumerate(self._eq_labels): + if nobs[i] != nobs[0]: + raise ValueError( + "All equations must have the same number of observations. " + f"The dependent variable in equation {label!r} has {nobs[i]} " + f"observations, but in equation {self._eq_labels[0]!r} it has " + f"{nobs[0]}." + ) + reference: Index | None = None + reference_name = "" + for i, label in enumerate(self._eq_labels): + dependent_name = f"dependent in equation {label!r}" + expected = row_labels(self._dependent[i]) + if expected is not None: + if reference is None: + reference, reference_name = expected, dependent_name + else: + check_row_alignment( + dependent_name, expected, reference, reference_name + ) + inputs = { + f"exog in equation {label!r}": self._exog[i], + f"endog in equation {label!r}": self._endog[i], + f"instruments in equation {label!r}": self._instr[i], + f"weights in equation {label!r}": self._weights[i], + } + check_alignment(expected, inputs, reference=dependent_name) + + def _check_clusters(self, config: Mapping[str, Any]) -> None: + """Check that clusters with labels use the labels of dependent""" + for i, label in enumerate(self._eq_labels): + expected = row_labels(self._dependent[i]) + if expected is not None: + check_alignment( + expected, + {"clusters": config.get("clusters")}, + reference=f"dependent in equation {label!r}", + ) + return + def _drop_missing(self) -> None: k = len(self._dependent) nobs = self._dependent[0].shape[0] @@ -1308,6 +1357,7 @@ def fit( linearmodels.system.covariance.KernelCovariance linearmodels.system.covariance.ClusteredCovariance """ + self._check_clusters(cov_config) if method is None: method = ( "ols" if (self._common_exog and self._constraints is None) else "gls" @@ -1463,6 +1513,19 @@ class IV3SLS(_LSSystemModelBase): \Omega = \Sigma \otimes I_N where :math:`\Sigma` is the covariance matrix of the residuals. + + Rows are paired by position, both within an equation and across the + equations. An input that carries labels (a pandas Series or DataFrame + with an index other than the default 0, 1, ..., n-1, or an xarray + DataArray with coordinates) must use the same index, in the same order, + as the dependent variable of its equation, so that a value cannot be + attached to the wrong observation. The dependent variables that carry + labels must also use the same index in every equation. A ``ValueError`` + is raised if an index differs. Either reindex the input to the index of + the dependent variable, or call ``to_numpy()`` to pair its values by + position. NumPy arrays and pandas objects with a default index do not + carry labels and are always paired by position. + The same rule applies to ``clusters`` passed to ``fit``. """ def __init__( @@ -1678,6 +1741,19 @@ class SUR(_LSSystemModelBase): SUR is a special case of 3SLS where there are no endogenous regressors and no instruments. + + Rows are paired by position, both within an equation and across the + equations. An input that carries labels (a pandas Series or DataFrame + with an index other than the default 0, 1, ..., n-1, or an xarray + DataArray with coordinates) must use the same index, in the same order, + as the dependent variable of its equation, so that a value cannot be + attached to the wrong observation. The dependent variables that carry + labels must also use the same index in every equation. A ``ValueError`` + is raised if an index differs. Either reindex the input to the index of + the dependent variable, or call ``to_numpy()`` to pair its values by + position. NumPy arrays and pandas objects with a default index do not + carry labels and are always paired by position. + The same rule applies to ``clusters`` passed to ``fit``. """ def __init__( @@ -1901,6 +1977,18 @@ class IVSystemGMM(_SystemModelBase): \hat{\beta}_{GMM} & = (X'ZW^{-1}Z'X)^{-1}X'ZW^{-1}Z'Y \\ where :math:`W` is a positive definite weighting matrix. + + Rows are paired by position, both within an equation and across the + equations. An input that carries labels (a pandas Series or DataFrame + with an index other than the default 0, 1, ..., n-1, or an xarray + DataArray with coordinates) must use the same index, in the same order, + as the dependent variable of its equation, so that a value cannot be + attached to the wrong observation. The dependent variables that carry + labels must also use the same index in every equation. A ``ValueError`` + is raised if an index differs. Either reindex the input to the index of + the dependent variable, or call ``to_numpy()`` to pair its values by + position. NumPy arrays and pandas objects with a default index do not + carry labels and are always paired by position. """ def __init__( diff --git a/linearmodels/tests/asset_pricing/test_input_alignment.py b/linearmodels/tests/asset_pricing/test_input_alignment.py new file mode 100644 index 0000000000..5910978958 --- /dev/null +++ b/linearmodels/tests/asset_pricing/test_input_alignment.py @@ -0,0 +1,324 @@ +""" +Factors that carry labels (pandas, xarray) have to use the labels of portfolios. + +The factor models pair the rows of portfolios and factors by position: the +return on the test portfolios in a period is paired with the value of the +factors in the same row. Factors whose rows are in a different order from +portfolios would silently be attached to the wrong periods. These tests check +that this is rejected, that inputs that carry no labels (NumPy arrays and +pandas objects with a default index) are paired by position, and that inputs +that use the labels of portfolios are accepted. +""" + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.asset_pricing import ( + LinearFactorModel, + LinearFactorModelGMM, + TradedFactorModel, +) +from linearmodels.iv.data import IVData +from linearmodels.shared.exceptions import MissingValueWarning + +NOBS = 120 +NFACTOR = 2 +NPORTFOLIO = 6 +MODELS = [TradedFactorModel, LinearFactorModel, LinearFactorModelGMM] +MODEL_IDS = [m.__name__ for m in MODELS] +FITTED = [TradedFactorModel, LinearFactorModel] +FITTED_IDS = [m.__name__ for m in FITTED] +INDEX_STYLES = ["datetime", "str", "int"] + + +def make_index(style, nobs=NOBS): + rng = np.random.default_rng(77) + if style == "str": + return pd.Index([f"t{i}" for i in range(nobs)], name="period") + if style == "datetime": + return pd.date_range("2000-01-31", periods=nobs, freq="D", name="date") + return pd.Index(rng.permutation(10 * nobs)[:nobs] + 5, name="id") + + +class Data: + """Labelled returns. portfolios and factors share an index""" + + def __init__(self, style="datetime", nobs=NOBS, seed=2): + rng = np.random.default_rng(seed) + self.index = idx = make_index(style, nobs) + self.factors = pd.DataFrame( + 0.5 + rng.standard_normal((nobs, NFACTOR)), + index=idx, + columns=[f"f{i}" for i in range(NFACTOR)], + ) + loadings = rng.chisquare(2, (NFACTOR, NPORTFOLIO)) / 2 + noise = rng.standard_normal((nobs, NPORTFOLIO)) + self.portfolios = pd.DataFrame( + self.factors.to_numpy() @ loadings + noise, + index=idx, + columns=[f"p{i}" for i in range(NPORTFOLIO)], + ) + + def frame(self): + return pd.concat([self.portfolios, self.factors], axis=1) + + +@pytest.fixture(params=INDEX_STYLES) +def data(request): + return Data(request.param) + + +@pytest.fixture +def simple(): + return Data("datetime") + + +def reverse(value): + """The same labels in the opposite order""" + return value.iloc[::-1] + + +def relabel(value): + """The same values, with the labels attached to the wrong rows""" + out = value.copy() + out.index = value.index[np.roll(np.arange(len(value)), 1)] + return out + + +CORRUPTIONS = {"reverse": reverse, "relabel": relabel} + + +def arrays(mod): + return [mod.portfolios.ndarray, mod.factors.ndarray] + + +def assert_same_arrays(left, right): + for a, b in zip(arrays(left), arrays(right), strict=True): + assert_allclose(a, b) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +def test_misaligned_factors_rejected(data, cls, corruption): + cls(data.portfolios, data.factors) + bad = CORRUPTIONS[corruption](data.factors) + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(data.portfolios, bad) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_message_names_input_and_remedy(simple, cls): + with pytest.raises(ValueError, match="does not match") as exc: + cls(simple.portfolios, reverse(simple.factors)) + message = str(exc.value) + assert "factors" in message + assert "does not match portfolios" in message + assert "in the same order as portfolios" in message + assert "reindex" in message + assert ".to_numpy()" in message + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +def test_misaligned_portfolios_rejected(simple, cls, corruption): + # portfolios sets the labels, so it is factors that do not match + bad = CORRUPTIONS[corruption](simple.portfolios) + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(bad, simple.factors) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_series_factor(simple, cls): + # A single factor is a Series. The portfolios are used with their labels + factor = simple.factors["f0"] + portfolios = simple.portfolios.iloc[:, :3] + cls(portfolios, factor) + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(portfolios, reverse(factor)) + mod = cls(portfolios, factor.to_numpy()[::-1]) + assert_allclose(mod.factors.ndarray, factor.to_numpy()[::-1, None]) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_same_labels_match_numpy(data, cls): + labelled = cls(data.portfolios, data.factors) + raw = cls(data.portfolios.to_numpy(), data.factors.to_numpy()) + assert_same_arrays(labelled, raw) + + +@pytest.mark.parametrize("cls", FITTED, ids=FITTED_IDS) +def test_same_labels_match_numpy_results(simple, cls): + labelled = cls(simple.portfolios, simple.factors).fit() + raw = cls(simple.portfolios.to_numpy(), simple.factors.to_numpy()).fit() + assert_allclose(labelled.risk_premia, raw.risk_premia) + assert_allclose(labelled.risk_premia_se, raw.risk_premia_se) + assert_allclose(labelled.alphas, raw.alphas) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_numpy_factors_are_paired_by_position(simple, cls): + values = simple.factors.to_numpy()[::-1] + mod = cls(simple.portfolios, values) + raw = cls(simple.portfolios.to_numpy(), values) + assert_same_arrays(mod, raw) + assert_allclose(mod.factors.ndarray, simple.factors.to_numpy()[::-1]) + + +@pytest.mark.parametrize("cls", FITTED, ids=FITTED_IDS) +def test_numpy_factors_are_paired_by_position_results(simple, cls): + values = simple.factors.to_numpy()[::-1] + res = cls(simple.portfolios, values).fit() + raw = cls(simple.portfolios.to_numpy(), values).fit() + assert_allclose(res.risk_premia, raw.risk_premia) + assert_allclose(res.alphas, raw.alphas) + labelled = cls(simple.portfolios, simple.factors).fit() + assert not np.allclose(res.alphas, labelled.alphas) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_numpy_portfolios_are_paired_by_position(simple, cls): + # Without labels in portfolios there is nothing to compare to + values = simple.portfolios.to_numpy() + for corrupt in CORRUPTIONS.values(): + mod = cls(values, corrupt(simple.factors)) + assert_allclose(mod.factors.ndarray, corrupt(simple.factors).to_numpy()) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_default_index_is_paired_by_position(simple, cls): + # pandas objects with a default index have no labels. This includes a + # DataFrame that has named columns + values = simple.factors.to_numpy()[::-1] + factors = pd.DataFrame(values, columns=simple.factors.columns) + assert isinstance(factors.index, pd.RangeIndex) + mod = cls(simple.portfolios, factors) + assert_allclose(mod.factors.ndarray, values) + # Without labels in portfolios there is nothing to compare to + portfolios = pd.DataFrame(simple.portfolios.to_numpy(), columns=list("abcdef")) + mod = cls(portfolios, reverse(simple.factors)) + assert_allclose(mod.factors.ndarray, reverse(simple.factors).to_numpy()) + # The labels of a reversed range are not 0, 1, ... + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(simple.portfolios, factors.iloc[::-1]) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_everything_reordered_together_is_allowed(simple, cls): + expected = cls(simple.portfolios, simple.factors) + reordered = cls(reverse(simple.portfolios), reverse(simple.factors)) + for a, b in zip(arrays(expected), arrays(reordered), strict=True): + assert_allclose(a, b[::-1]) + + +@pytest.mark.parametrize("cls", FITTED, ids=FITTED_IDS) +def test_everything_reordered_together_results(simple, cls): + expected = cls(simple.portfolios, simple.factors).fit() + reordered = cls(reverse(simple.portfolios), reverse(simple.factors)).fit() + assert_allclose(expected.risk_premia, reordered.risk_premia) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_index_names_are_ignored(simple, cls): + factors = simple.factors.rename_axis("another name") + mod = cls(simple.portfolios, factors) + assert_allclose(mod.factors.ndarray, simple.factors.to_numpy()) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_different_length_reports_the_length(simple, cls): + with pytest.raises(ValueError, match="number of observations in portfolios"): + cls(simple.portfolios, simple.factors.iloc[:-1]) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_missing_values_do_not_change_the_check(simple, cls): + portfolios = simple.portfolios.copy() + portfolios.iloc[5, 2] = np.nan + with pytest.warns(MissingValueWarning): + mod = cls(portfolios, simple.factors) + assert mod.portfolios.shape[0] == NOBS - 1 + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(portfolios, reverse(simple.factors)) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_ivdata_inputs(simple, cls): + cls(IVData(simple.portfolios), IVData(simple.factors)) + with pytest.raises(ValueError, match="row index of factors"): + cls(simple.portfolios, IVData(reverse(simple.factors))) + with pytest.raises(ValueError, match="row index of factors"): + cls(IVData(simple.portfolios), reverse(simple.factors)) + # IVData built from an array has no labels + mod = cls(simple.portfolios, IVData(simple.factors.to_numpy()[::-1])) + assert_allclose(mod.factors.ndarray, simple.factors.to_numpy()[::-1]) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_xarray_inputs(simple, cls): + xr = pytest.importorskip("xarray") + + def to_xarray(frame, coords=True): + obs = {"obs": frame.index.to_numpy()} if coords else {} + return xr.DataArray( + frame.to_numpy(), dims=["obs", "var"], coords={**obs, "var": list(frame)} + ) + + factors = to_xarray(simple.factors) + mod = cls(simple.portfolios, factors) + assert_allclose(mod.factors.ndarray, simple.factors.to_numpy()) + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(simple.portfolios, factors.isel(obs=slice(None, None, -1))) + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(to_xarray(simple.portfolios), reverse(simple.factors)) + cls(to_xarray(simple.portfolios), factors) + # No coordinates, no labels + plain = to_xarray(reverse(simple.factors), coords=False) + mod = cls(simple.portfolios, plain) + assert_allclose(mod.factors.ndarray, reverse(simple.factors).to_numpy()) + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_options_do_not_change_the_check(simple, cls): + kwargs = {"risk_free": True} + if cls is LinearFactorModel: + kwargs["sigma"] = np.eye(NPORTFOLIO) + if cls is TradedFactorModel: + kwargs = {} + cls(simple.portfolios, simple.factors, **kwargs) + with pytest.raises(ValueError, match="row index of factors does not match"): + cls(simple.portfolios, reverse(simple.factors), **kwargs) + + +# Formulas use a single DataFrame + + +@pytest.mark.parametrize("cls", MODELS, ids=MODEL_IDS) +def test_from_formula(simple, cls): + frame = simple.frame() + assert not isinstance(frame.index, pd.RangeIndex) + formula = "p0 + p1 + p2 + p3 ~ f0 + f1" + mod = cls.from_formula(formula, frame) + raw = cls(simple.portfolios.iloc[:, :4].to_numpy(), simple.factors.to_numpy()) + assert_same_arrays(mod, raw) + # Portfolios in a separate frame + portfolios = simple.portfolios + mod = cls.from_formula("f0 + f1", frame, portfolios=portfolios) + assert_same_arrays(mod, cls(portfolios.to_numpy(), simple.factors.to_numpy())) + for corrupt in CORRUPTIONS.values(): + with pytest.raises(ValueError, match="row index of factors does not match"): + cls.from_formula("f0 + f1", frame, portfolios=corrupt(portfolios)) + # The data has a default index, so there are no labels + default = frame.reset_index(drop=True) + cls.from_formula( + "f0 + f1", default, portfolios=reverse(default[portfolios.columns]) + ) + + +@pytest.mark.parametrize("cls", FITTED, ids=FITTED_IDS) +def test_from_formula_matches_arrays(simple, cls): + frame = simple.frame() + res = cls.from_formula("p0 + p1 + p2 + p3 ~ f0 + f1", frame).fit() + raw = cls(simple.portfolios.iloc[:, :4].to_numpy(), simple.factors.to_numpy()).fit() + assert_allclose(res.risk_premia, raw.risk_premia) diff --git a/linearmodels/tests/iv/test_absorbing.py b/linearmodels/tests/iv/test_absorbing.py index ecb0203c05..b5dacb8d2d 100644 --- a/linearmodels/tests/iv/test_absorbing.py +++ b/linearmodels/tests/iv/test_absorbing.py @@ -762,3 +762,72 @@ def test_lsmr_annihilate_empty(): y = np.empty((1000, 0)) y_out = lsmr_annihilate(x, y) assert y_out.shape == y.shape + + +@pytest.mark.parametrize( + "labels", ["both", "cat", "cont", "neither", "cat-only", "cont-only"] +) +def test_interaction_without_labels_pairs_by_position(labels): + # Only cat or cont is labelled, or one is missing. Missing flags are paired + # by position and do not depend on the labels + rng = np.random.default_rng(0) + nobs = 30 + index = pd.Index([f"r{i}" for i in range(nobs)]) + cat = pd.DataFrame({"g": pd.Categorical(rng.integers(0, 4, nobs))}, index=index) + cont = pd.DataFrame({"c": rng.standard_normal(nobs)}, index=index) + cat.iloc[5, 0] = np.nan + cont.iloc[11, 0] = np.nan + cats = {"both": cat, "cat": cat, "cont": cat.to_numpy(), "neither": cat.to_numpy()} + conts = { + "both": cont, + "cat": cont.to_numpy(), + "cont": cont, + "neither": cont.to_numpy(), + } + if labels == "cat-only": + interaction = Interaction(cat=cat) + expected = np.zeros(nobs, dtype=bool) + expected[5] = True + elif labels == "cont-only": + interaction = Interaction(cont=cont) + expected = np.zeros(nobs, dtype=bool) + expected[11] = True + else: + interaction = Interaction(cat=cats[labels], cont=conts[labels]) + expected = np.zeros(nobs, dtype=bool) + expected[[5, 11]] = True + isnull = interaction.isnull + assert isnull.shape == (nobs,) + assert_array_equal(isnull.to_numpy(), expected) + interaction.drop(expected) + assert interaction.cat.shape[0] == interaction.cont.shape[0] + + +@pytest.mark.parametrize("labels", ["cat", "cont", "cat-only", "cont-only"]) +def test_interactions_without_labels(labels): + # AbsorbingLS pairs a labelled dependent variable with interactions that + # only partly have labels by position + rng = np.random.default_rng(1) + nobs = 60 + index = pd.Index([f"r{i}" for i in range(nobs)]) + y = pd.Series(rng.standard_normal(nobs), index=index) + x = pd.DataFrame({"x": rng.standard_normal(nobs)}, index=index) + cat = pd.DataFrame({"g": pd.Categorical(rng.integers(0, 4, nobs))}, index=index) + cont = pd.DataFrame({"c": rng.standard_normal(nobs)}, index=index) + interactions = { + "cat": Interaction(cat=cat.to_numpy(), cont=cont), + "cont": Interaction(cat=cat, cont=cont.to_numpy()), + "cat-only": Interaction(cat=cat), + "cont-only": Interaction(cont=cont), + } + mod = AbsorbingLS(y, x, interactions=interactions[labels]) + res = mod.fit() + raw = AbsorbingLS( + y.to_numpy(), + x.to_numpy(), + interactions=Interaction( + cat=None if labels == "cont-only" else cat.to_numpy(), + cont=None if labels == "cat-only" else cont.to_numpy(), + ), + ).fit() + assert_allclose(res.params, raw.params) diff --git a/linearmodels/tests/iv/test_input_alignment.py b/linearmodels/tests/iv/test_input_alignment.py new file mode 100644 index 0000000000..f5963ef1c1 --- /dev/null +++ b/linearmodels/tests/iv/test_input_alignment.py @@ -0,0 +1,851 @@ +""" +Inputs that carry labels (pandas, xarray) have to use the labels of dependent. + +The IV and absorbing estimators pair the rows of their inputs by position. A +labelled input whose rows are in a different order from dependent would +silently be attached to the wrong observations. These tests check that this +is rejected, that inputs that carry no labels (NumPy arrays and pandas objects +with a default index) are paired by position, and that inputs that use the +labels of dependent are accepted. +""" + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels import OLS +from linearmodels.iv import ( + IV2SLS, + IVGMM, + IVGMMCUE, + IVJIVE, + IVLIML, + AbsorbingLS, + Interaction, +) +from linearmodels.iv.data import IVData +from linearmodels.iv.model import _OLS +from linearmodels.shared.exceptions import MissingValueWarning + +NOBS = 120 +ESTIMATORS = [IV2SLS, IVLIML, IVGMM, IVGMMCUE, IVJIVE] +ESTIMATOR_IDS = [e.__name__ for e in ESTIMATORS] +INDEX_STYLES = ["str", "datetime", "int", "multi"] +INPUTS = ["exog", "endog", "instruments", "weights"] +PARTS = ["dependent", "exog", "endog", "instruments"] +# Estimators that are fast to fit +FAST = [IV2SLS, IVLIML, IVGMM, IVJIVE] +FAST_IDS = [e.__name__ for e in FAST] +GMM = [IVGMM, IVGMMCUE] +GMM_IDS = [e.__name__ for e in GMM] +WITH_CLUSTERS = [*FAST, IVGMMCUE] +WITH_CLUSTERS_IDS = [e.__name__ for e in WITH_CLUSTERS] + + +def make_index(style, nobs=NOBS): + rng = np.random.default_rng(99) + if style == "str": + return pd.Index([f"r{i}" for i in range(nobs)], name="obs") + if style == "datetime": + return pd.date_range("2000-01-01", periods=nobs, freq="D", name="date") + if style == "int": + return pd.Index(rng.permutation(10 * nobs)[:nobs] + 7, name="id") + return pd.MultiIndex.from_product( + [np.arange(nobs // 10), np.arange(10)], names=["firm", "year"] + ) + + +class Data: + """Labelled IV data. The attributes are pandas objects with one index.""" + + def __init__(self, style="str", nobs=NOBS, seed=0): + rng = np.random.default_rng(seed) + self.index = idx = make_index(style, nobs) + self.instruments = pd.DataFrame( + rng.standard_normal((nobs, 2)), index=idx, columns=["z1", "z2"] + ) + z = self.instruments + self.endog = pd.DataFrame( + {"x": z.z1 + z.z2 + rng.standard_normal(nobs)}, index=idx + ) + self.exog = pd.DataFrame( + {"const": 1.0, "w": rng.standard_normal(nobs)}, index=idx + ) + self.dependent = pd.Series( + 1 + self.endog.x + self.exog.w + rng.standard_normal(nobs), + index=idx, + name="y", + ) + self.weights = pd.Series( + np.exp(rng.standard_normal(nobs)), index=idx, name="wt" + ) + self.clusters = pd.Series(rng.integers(0, 12, nobs), index=idx, name="cl") + + def parts(self): + return {name: getattr(self, name) for name in PARTS} + + def raw(self): + """The arguments of an IV model as NumPy arrays""" + return [value.to_numpy() for value in self.parts().values()] + + +@pytest.fixture(params=INDEX_STYLES) +def data(request): + return Data(request.param) + + +@pytest.fixture +def simple(): + return Data("str") + + +def reverse(value): + """The same labels in the opposite order""" + return value.iloc[::-1] + + +def relabel(value): + """The same values, with the labels attached to the wrong rows""" + out = value.copy() + out.index = value.index[np.roll(np.arange(len(value)), 1)] + return out + + +CORRUPTIONS = {"reverse": reverse, "relabel": relabel} + + +def build(cls, data, **replace): + """Model using data, where some inputs are replaced""" + weights = replace.pop("weights", None) + parts = {**data.parts(), **replace} + return cls(*parts.values(), weights=weights) + + +def arrays(mod): + """The arrays used to estimate a model""" + return [ + mod.dependent.ndarray, + mod.exog.ndarray, + mod.endog.ndarray, + mod.instruments.ndarray, + mod.weights.ndarray, + ] + + +def assert_same_arrays(left, right): + for a, b in zip(arrays(left), arrays(right), strict=True): + assert_allclose(a, b) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +@pytest.mark.parametrize("name", INPUTS) +def test_misaligned_input_rejected(data, cls, name, corruption): + # The model is accepted when the input uses the labels of dependent + build(cls, data, **{name: getattr(data, name)}) + bad = CORRUPTIONS[corruption](getattr(data, name)) + with pytest.raises(ValueError, match=f"row index of {name} does not match"): + build(cls, data, **{name: bad}) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_message_names_input_and_remedy(simple, cls): + with pytest.raises(ValueError, match="does not match") as exc: + build(cls, simple, instruments=reverse(simple.instruments)) + message = str(exc.value) + assert "instruments" in message + assert "reindex" in message + assert ".to_numpy()" in message + assert "in the same order as dependent" in message + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_dependent_reordered_rejected(simple, cls): + # It is the other inputs that do not match dependent + for corrupt in CORRUPTIONS.values(): + with pytest.raises(ValueError, match="row index of exog does not match"): + build(cls, simple, dependent=corrupt(simple.dependent)) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_everything_reordered_together_is_allowed(simple, cls): + expected = build(cls, simple, weights=simple.weights) + flipped = {name: reverse(value) for name, value in simple.parts().items()} + reordered = cls(*flipped.values(), weights=reverse(simple.weights)) + # The same observations in the opposite order + for a, b in zip(arrays(expected), arrays(reordered), strict=True): + assert_allclose(a, b[::-1]) + if cls is not IVGMMCUE: + assert_allclose(expected.fit().params, reordered.fit().params, rtol=1e-8) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_same_labels_match_numpy(data, cls): + labelled = build(cls, data, weights=data.weights) + raw = cls(*data.raw(), weights=data.weights.to_numpy()) + assert_same_arrays(labelled, raw) + + +@pytest.mark.parametrize("cls", FAST, ids=FAST_IDS) +def test_same_labels_match_numpy_results(simple, cls): + labelled = build(cls, simple, weights=simple.weights).fit() + raw = cls(*simple.raw(), weights=simple.weights.to_numpy()).fit() + assert_allclose(labelled.params, raw.params) + assert_allclose(labelled.std_errors, raw.std_errors) + assert_allclose(labelled.resids, raw.resids) + + +def numpy_equivalent(cls, data, name, values): + """Model where the labels of all inputs have been removed""" + raw = {part: value.to_numpy() for part, value in data.parts().items()} + if name == "weights": + return cls(*raw.values(), weights=values) + raw[name] = values + return cls(*raw.values()) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +@pytest.mark.parametrize("name", INPUTS) +def test_numpy_input_is_paired_by_position(simple, cls, name): + # A reordered array has no labels, and so it is paired by position + original = getattr(simple, name) + values = original.to_numpy()[::-1] + mod = build(cls, simple, **{name: values}) + assert_same_arrays(mod, numpy_equivalent(cls, simple, name, values)) + if name != "weights": + # The values were paired with the opposite rows + assert_allclose( + getattr(mod, name).ndarray, np.asarray(original).reshape(NOBS, -1)[::-1] + ) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +@pytest.mark.parametrize("name", INPUTS) +def test_default_index_is_paired_by_position(simple, cls, name): + # pandas objects that have a default index carry no labels. This is also + # the case for a DataFrame that has named columns + original = getattr(simple, name) + values = original.to_numpy()[::-1] + if isinstance(original, pd.DataFrame): + framed = pd.DataFrame(values, columns=original.columns) + else: + framed = pd.Series(values, name=original.name) + assert isinstance(framed.index, pd.RangeIndex) + mod = build(cls, simple, **{name: framed}) + assert_same_arrays(mod, numpy_equivalent(cls, simple, name, values)) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_reversed_default_index_has_labels(simple, cls): + # The labels of a reversed range are not 0, 1, ... + exog = pd.DataFrame(simple.exog.to_numpy(), columns=simple.exog.columns) + build(cls, simple, exog=exog) + reversed_exog = exog.iloc[::-1] + assert isinstance(reversed_exog.index, pd.RangeIndex) + with pytest.raises(ValueError, match="row index of exog does not match"): + build(cls, simple, exog=reversed_exog) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_unlabelled_dependent_is_not_checked(simple, cls): + # There are no labels to compare to, so the inputs are paired by position + y = simple.dependent.to_numpy() + for corrupt in CORRUPTIONS.values(): + mod = cls( + y, + corrupt(simple.exog), + corrupt(simple.endog), + corrupt(simple.instruments), + weights=corrupt(simple.weights), + ) + assert_allclose(mod.exog.ndarray, corrupt(simple.exog).to_numpy()) + assert_allclose(mod.dependent.ndarray, y[:, None]) + # The same holds if dependent has a default index + cls(pd.Series(y), reverse(simple.exog), simple.endog, simple.instruments) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_index_names_are_ignored(simple, cls): + exog = simple.exog.rename_axis("another name") + weights = simple.weights.rename_axis(None) + mod = build(cls, simple, exog=exog, weights=weights) + assert_allclose(mod.exog.ndarray, simple.exog.to_numpy()) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_duplicate_labels(cls): + # Duplicated labels are equal if they are in the same order + rng = np.random.default_rng(3) + index = pd.Index(np.repeat(["a", "b", "c", "d"], NOBS // 4)) + exog = pd.DataFrame({"c": 1.0, "w": rng.standard_normal(NOBS)}, index=index) + z = pd.DataFrame(rng.standard_normal((NOBS, 2)), index=index) + x = pd.Series(z.sum(axis=1) + rng.standard_normal(NOBS), index=index, name="x") + y = pd.Series(x + rng.standard_normal(NOBS), index=index, name="y") + cls(y, exog, x, z) + with pytest.raises(ValueError, match="row index of instruments"): + cls(y, exog, x, z.iloc[::-1]) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_different_length_reports_the_length(simple, cls): + with pytest.raises(ValueError, match="Array required to have"): + build(cls, simple, exog=simple.exog.iloc[:-1]) + with pytest.raises(ValueError, match="Array required to have"): + build(cls, simple, weights=reverse(simple.weights).iloc[:-1]) + + +def test_ols(simple): + # linearmodels.OLS is the same as _OLS + assert OLS is _OLS + mod = IV2SLS(simple.dependent, simple.exog, None, None) + assert mod._method == "OLS" + with pytest.raises(ValueError, match="row index of exog does not match"): + IV2SLS(simple.dependent, reverse(simple.exog), None, None) + with pytest.raises(ValueError, match="row index of weights does not match"): + IV2SLS( + simple.dependent, simple.exog, None, None, weights=reverse(simple.weights) + ) + with pytest.raises(ValueError, match="row index of exog does not match"): + OLS(simple.dependent, reverse(simple.exog)) + with pytest.raises(ValueError, match="row index of weights does not match"): + OLS(simple.dependent, simple.exog, weights=relabel(simple.weights)) + res = OLS(simple.dependent, simple.exog).fit( + cov_type="clustered", clusters=simple.clusters + ) + raw = OLS(simple.dependent, simple.exog).fit( + cov_type="clustered", clusters=simple.clusters.to_numpy() + ) + assert_allclose(res.std_errors, raw.std_errors) + with pytest.raises(ValueError, match="row index of clusters does not match"): + OLS(simple.dependent, simple.exog).fit( + cov_type="clustered", clusters=reverse(simple.clusters) + ) + labelled = _OLS(simple.dependent, simple.exog, weights=simple.weights) + raw = _OLS( + simple.dependent.to_numpy(), + simple.exog.to_numpy(), + weights=simple.weights.to_numpy(), + ) + assert_same_arrays(labelled, raw) + + +def test_inputs_that_are_not_present_are_not_compared(simple): + IV2SLS(simple.dependent, None, simple.endog, simple.instruments) + with pytest.raises(ValueError, match="row index of instruments"): + IV2SLS(simple.dependent, None, simple.endog, reverse(simple.instruments)) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_ivdata_inputs(simple, cls): + build(cls, simple, instruments=IVData(simple.instruments)) + with pytest.raises(ValueError, match="row index of instruments"): + build(cls, simple, instruments=IVData(reverse(simple.instruments))) + # IVData built from an array carries no labels + build(cls, simple, instruments=IVData(simple.instruments.to_numpy()[::-1])) + with pytest.raises(ValueError, match="row index of exog"): + build( + cls, + simple, + dependent=IVData(simple.dependent), + exog=IVData(reverse(simple.exog)), + ) + with pytest.raises(ValueError, match="row index of weights"): + build(cls, simple, weights=IVData(reverse(simple.weights))) + # IVData built from a pandas object with a default index carries no labels + default = IVData(pd.DataFrame(simple.exog.to_numpy()[::-1], columns=["a", "b"])) + build(cls, simple, exog=default) + + +def to_xarray(frame, coords=True): + xr = pytest.importorskip("xarray") + frame = frame.to_frame() if isinstance(frame, pd.Series) else frame + obs = {"obs": frame.index.to_numpy()} if coords else {} + return xr.DataArray( + frame.to_numpy(), dims=["obs", "var"], coords={**obs, "var": list(frame)} + ) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_xarray_inputs(simple, cls): + instruments = to_xarray(simple.instruments) + build(cls, simple, instruments=instruments) + # Reversing the rows of the DataArray reverses the coordinates + reversed_rows = instruments.isel(obs=slice(None, None, -1)) + with pytest.raises(ValueError, match="row index of instruments"): + build(cls, simple, instruments=reversed_rows) + with pytest.raises(ValueError, match="row index of exog"): + build( + cls, + simple, + dependent=to_xarray(simple.dependent), + exog=to_xarray(reverse(simple.exog)), + ) + # A DataArray without coordinates has no labels + no_coords = to_xarray(reverse(simple.instruments), coords=False) + mod = build(cls, simple, instruments=no_coords) + assert_allclose(mod.instruments.ndarray, simple.instruments.to_numpy()[::-1]) + # An xarray dependent sets the labels + with pytest.raises(ValueError, match="row index of endog"): + build( + cls, + simple, + dependent=to_xarray(simple.dependent), + endog=reverse(simple.endog), + ) + build(cls, simple, dependent=to_xarray(simple.dependent), endog=simple.endog) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_missing_values_do_not_change_the_check(simple, cls): + exog = simple.exog.copy() + exog.iloc[3, 1] = np.nan + weights = simple.weights.copy() + weights.iloc[10] = np.nan + with pytest.warns(MissingValueWarning): + mod = build(cls, simple, exog=exog, weights=weights) + assert mod.notnull.sum() == NOBS - 2 + with pytest.raises(ValueError, match="row index of instruments"): + build( + cls, + simple, + exog=exog, + weights=weights, + instruments=reverse(simple.instruments), + ) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_columns_are_not_compared(simple, cls): + # The columns are variables. Only the rows are paired + exog = simple.exog.copy() + exog["group"] = pd.Categorical(np.arange(NOBS) % 3) + mod = build(cls, simple, exog=exog) + assert mod.exog.shape[1] == 4 + with pytest.raises(ValueError, match="row index of exog"): + build(cls, simple, exog=reverse(exog)) + + +# Clusters + + +@pytest.mark.parametrize("cls", WITH_CLUSTERS, ids=WITH_CLUSTERS_IDS) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +def test_clusters_in_fit(simple, cls, corruption): + mod = build(cls, simple) + res = mod.fit(cov_type="clustered", clusters=simple.clusters) + raw = mod.fit(cov_type="clustered", clusters=simple.clusters.to_numpy()) + assert_allclose(res.std_errors, raw.std_errors) + bad = CORRUPTIONS[corruption](simple.clusters) + with pytest.raises(ValueError, match="row index of clusters does not match"): + mod.fit(cov_type="clustered", clusters=bad) + + +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +def test_clusters_index_styles(data, corruption): + mod = IV2SLS(*data.parts().values()) + res = mod.fit(cov_type="clustered", clusters=data.clusters) + raw = mod.fit(cov_type="clustered", clusters=data.clusters.to_numpy()) + assert_allclose(res.std_errors, raw.std_errors) + bad = CORRUPTIONS[corruption](data.clusters) + with pytest.raises(ValueError, match="row index of clusters does not match"): + mod.fit(cov_type="clustered", clusters=bad) + + +@pytest.mark.parametrize("cls", WITH_CLUSTERS, ids=WITH_CLUSTERS_IDS) +def test_clusters_numpy_is_paired_by_position(simple, cls): + mod = build(cls, simple) + values = simple.clusters.to_numpy()[::-1] + raw = cls(*simple.raw()).fit(cov_type="clustered", clusters=values) + res = mod.fit(cov_type="clustered", clusters=values) + assert_allclose(res.std_errors, raw.std_errors) + # A default index carries no labels + frame = pd.DataFrame({"cl": values}) + res = mod.fit(cov_type="clustered", clusters=frame) + assert_allclose(res.std_errors, raw.std_errors) + res = mod.fit(cov_type="clustered", clusters=pd.Series(values)) + assert_allclose(res.std_errors, raw.std_errors) + + +@pytest.mark.parametrize("cls", [IV2SLS, IVLIML]) +def test_two_way_clusters(simple, cls): + rng = np.random.default_rng(11) + second = pd.Series(rng.integers(0, 7, NOBS), index=simple.index, name="cl2") + clusters = pd.concat([simple.clusters, second], axis=1) + mod = build(cls, simple) + res = mod.fit(cov_type="clustered", clusters=clusters) + raw = mod.fit(cov_type="clustered", clusters=clusters.to_numpy()) + assert_allclose(res.cov, raw.cov) + with pytest.raises(ValueError, match="row index of clusters"): + mod.fit(cov_type="clustered", clusters=reverse(clusters)) + with pytest.raises(ValueError, match="row index of clusters"): + mod.fit(cov_type="clustered", clusters=relabel(clusters)) + # The columns are names, and do not need to match anything + unnamed = pd.DataFrame(clusters.to_numpy(), columns=["a", "b"]) + res = mod.fit(cov_type="clustered", clusters=unnamed) + assert_allclose(res.cov, raw.cov) + + +@pytest.mark.parametrize("cls", GMM, ids=GMM_IDS) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +def test_clusters_in_weight_config(simple, cls, corruption): + clustered = {"weight_type": "clustered", "clusters": simple.clusters} + labelled = cls(*simple.parts().values(), **clustered) + raw_clusters = {**clustered, "clusters": simple.clusters.to_numpy()} + raw = cls(*simple.raw(), **raw_clusters) + assert_allclose(labelled.fit().params, raw.fit().params) + bad = CORRUPTIONS[corruption](simple.clusters) + with pytest.raises(ValueError, match="row index of clusters does not match"): + cls(*simple.parts().values(), **{**clustered, "clusters": bad}) + # Arrays have no labels and are paired by position + cls(*simple.parts().values(), **{**clustered, "clusters": bad.to_numpy()}) + # As are pandas objects with a default index + cls( + *simple.parts().values(), + **{**clustered, "clusters": pd.Series(bad.to_numpy())}, + ) + + +@pytest.mark.parametrize("cls", GMM, ids=GMM_IDS) +def test_clusters_in_weight_config_and_fit(simple, cls): + clustered = {"weight_type": "clustered", "clusters": simple.clusters} + mod = cls(*simple.parts().values(), **clustered) + mod.fit(cov_type="clustered", clusters=simple.clusters) + with pytest.raises(ValueError, match="row index of clusters does not match"): + mod.fit(cov_type="clustered", clusters=reverse(simple.clusters)) + + +@pytest.mark.parametrize("cls", GMM, ids=GMM_IDS) +def test_other_weight_config_is_unchanged(simple, cls): + mod = build(cls, simple) + assert mod._weight_type == "robust" + mod = cls(*simple.parts().values(), weight_type="kernel", kernel="parzen") + assert mod._weight_type == "kernel" + assert mod._weight_config["kernel"] == "parzen" + mod = cls(*simple.parts().values(), weight_type="unadjusted", debiased=True) + assert mod._weight_config["debiased"] + assert mod._weight_kwargs["debiased"] + assert IVGMM(*simple.parts().values())._method == "IV-GMM" + + +@pytest.mark.parametrize("cls", WITH_CLUSTERS, ids=WITH_CLUSTERS_IDS) +def test_clusters_and_missing_values(simple, cls): + # Missing values are dropped when the model is created, and so clusters + # have to be provided for the observations that are used + exog = simple.exog.copy() + exog.iloc[7, 1] = np.nan + with pytest.warns(MissingValueWarning): + mod = build(cls, simple, exog=exog) + kept = simple.clusters[mod.notnull] + assert len(kept) == NOBS - 1 + res = mod.fit(cov_type="clustered", clusters=kept) + raw = mod.fit(cov_type="clustered", clusters=kept.to_numpy()) + assert_allclose(res.std_errors, raw.std_errors) + with pytest.raises(ValueError, match="row index of clusters"): + mod.fit(cov_type="clustered", clusters=reverse(kept)) + # Using all of the clusters is a length error, as before + with pytest.raises(ValueError, match="clusters has the wrong nobs"): + mod.fit(cov_type="clustered", clusters=simple.clusters) + + +@pytest.mark.parametrize("cls", GMM, ids=GMM_IDS) +def test_weight_clusters_and_missing_values(simple, cls): + exog = simple.exog.copy() + exog.iloc[7, 1] = np.nan + mask = np.ones(NOBS, dtype=bool) + mask[7] = False + kept = simple.clusters[mask] + parts = {**simple.parts(), "exog": exog} + with pytest.warns(MissingValueWarning): + mod = cls(*parts.values(), weight_type="clustered", clusters=kept) + assert mod.fit().nobs == NOBS - 1 + with ( + pytest.warns(MissingValueWarning), + pytest.raises(ValueError, match="row index of clusters"), + ): + cls(*parts.values(), weight_type="clustered", clusters=reverse(kept)) + + +# Formulas use a single DataFrame + + +FORMULA = "y ~ 1 + w + [x ~ z1 + z2]" + + +@pytest.fixture +def frame(simple): + return pd.concat( + [ + simple.dependent, + simple.exog[["w"]], + simple.endog, + simple.instruments, + simple.weights, + ], + axis=1, + ) + + +@pytest.mark.parametrize("cls", ESTIMATORS, ids=ESTIMATOR_IDS) +def test_from_formula(frame, cls): + assert not isinstance(frame.index, pd.RangeIndex) + plain = cls.from_formula(FORMULA, frame) + assert plain.exog.shape == (NOBS, 2) + weighted = cls.from_formula(FORMULA, frame, weights=frame["wt"]) + raw = cls.from_formula(FORMULA, frame, weights=frame["wt"].to_numpy()) + assert_same_arrays(weighted, raw) + for corrupt in CORRUPTIONS.values(): + with pytest.raises(ValueError, match="row index of weights does not match"): + cls.from_formula(FORMULA, frame, weights=corrupt(frame["wt"])) + # The data has a default index, so there are no labels + default = frame.reset_index(drop=True) + cls.from_formula(FORMULA, default, weights=default["wt"]) + cls.from_formula(FORMULA, default, weights=default["wt"].to_numpy()[::-1]) + + +@pytest.mark.parametrize("cls", FAST, ids=FAST_IDS) +def test_from_formula_matches_arrays(frame, cls): + formula = cls.from_formula(FORMULA, frame, weights=frame["wt"]).fit() + exog = np.c_[np.ones(NOBS), frame["w"].to_numpy()] + arrays_ = cls( + frame["y"].to_numpy(), + exog, + frame[["x"]].to_numpy(), + frame[["z1", "z2"]].to_numpy(), + weights=frame["wt"].to_numpy(), + ).fit() + assert_allclose(formula.params, arrays_.params) + + +# AbsorbingLS + + +class Absorbing: + def __init__(self, style="str", nobs=NOBS, seed=5): + rng = np.random.default_rng(seed) + self.index = idx = make_index(style, nobs) + self.exog = pd.DataFrame( + {"const": 1.0, "w": rng.standard_normal(nobs)}, index=idx + ) + self.cats = pd.DataFrame( + { + "g": pd.Categorical(rng.integers(0, 6, nobs)), + "h": pd.Categorical(rng.integers(0, 4, nobs)), + }, + index=idx, + ) + self.cont = pd.DataFrame( + rng.standard_normal((nobs, 2)), index=idx, columns=["c0", "c1"] + ) + self.dependent = pd.Series( + self.exog.w + self.cats["g"].cat.codes + rng.standard_normal(nobs), + index=idx, + name="y", + ) + self.weights = pd.Series( + np.exp(rng.standard_normal(nobs)), index=idx, name="wt" + ) + self.clusters = pd.Series(rng.integers(0, 12, nobs), index=idx, name="cl") + + def model(self, **kwargs): + return AbsorbingLS(self.dependent, self.exog, absorb=self.cats, **kwargs) + + +@pytest.fixture(params=INDEX_STYLES) +def absorbing(request): + return Absorbing(request.param) + + +@pytest.fixture +def abs_simple(): + return Absorbing("str") + + +def test_absorbing_misaligned_rejected(absorbing): + y, x, a, w = ( + absorbing.dependent, + absorbing.exog, + absorbing.cats, + absorbing.weights, + ) + AbsorbingLS(y, x, absorb=a, weights=w) + for corrupt in CORRUPTIONS.values(): + with pytest.raises(ValueError, match="row index of absorb does not match"): + AbsorbingLS(y, x, absorb=corrupt(a)) + with pytest.raises(ValueError, match="row index of exog does not match"): + AbsorbingLS(y, corrupt(x), absorb=a) + with pytest.raises(ValueError, match="row index of weights does not match"): + AbsorbingLS(y, x, absorb=a, weights=corrupt(w)) + with pytest.raises(ValueError, match="row index of exog does not match"): + AbsorbingLS(corrupt(y), x, absorb=a) + + +def test_absorbing_message(abs_simple): + with pytest.raises(ValueError, match="does not match") as exc: + AbsorbingLS( + abs_simple.dependent, abs_simple.exog, absorb=reverse(abs_simple.cats) + ) + message = str(exc.value) + assert "absorb" in message + assert "reindex" in message + assert ".to_numpy()" in message + + +def test_absorbing_continuous_absorb(abs_simple): + frame = pd.concat([abs_simple.cats, abs_simple.cont], axis=1) + AbsorbingLS(abs_simple.dependent, abs_simple.exog, absorb=frame) + with pytest.raises(ValueError, match="row index of absorb"): + AbsorbingLS(abs_simple.dependent, abs_simple.exog, absorb=reverse(frame)) + + +def test_absorbing_interaction_in_absorb(abs_simple): + y, x = abs_simple.dependent, abs_simple.exog + cat, cont = abs_simple.cats, abs_simple.cont + AbsorbingLS(y, x, absorb=Interaction(cat=cat, cont=cont)) + AbsorbingLS(y, x, absorb=Interaction(cat=cat)) + AbsorbingLS(y, x, absorb=Interaction(cont=cont)) + with pytest.raises(ValueError, match="row index of absorb"): + AbsorbingLS(y, x, absorb=Interaction(cat=reverse(cat), cont=reverse(cont))) + with pytest.raises(ValueError, match="row index of absorb"): + AbsorbingLS(y, x, absorb=Interaction(cat=reverse(cat))) + with pytest.raises(ValueError, match="row index of absorb"): + AbsorbingLS(y, x, absorb=Interaction(cont=reverse(cont))) + # Arrays have no labels + AbsorbingLS(y, x, absorb=Interaction(cat=cat.to_numpy()[::-1])) + + +def test_absorbing_interactions(abs_simple): + y, x = abs_simple.dependent, abs_simple.exog + cat, cont = abs_simple.cats, abs_simple.cont + frame = pd.concat([cat[["g"]], cont[["c0"]]], axis=1) + iact = Interaction(cat=cat[["h"]], cont=cont[["c1"]]) + flipped = Interaction(cat=reverse(cat[["h"]]), cont=reverse(cont[["c1"]])) + AbsorbingLS(y, x, interactions=frame) + AbsorbingLS(y, x, interactions=iact) + AbsorbingLS(y, x, interactions=[frame, iact]) + with pytest.raises(ValueError, match="row index of interactions does not match"): + AbsorbingLS(y, x, interactions=reverse(frame)) + with pytest.raises(ValueError, match="row index of interactions does not match"): + AbsorbingLS(y, x, interactions=flipped) + with pytest.raises(ValueError, match=r"row index of interactions\[0\]"): + AbsorbingLS(y, x, interactions=[reverse(frame), iact]) + with pytest.raises(ValueError, match=r"row index of interactions\[1\]"): + AbsorbingLS(y, x, interactions=[frame, relabel(frame)]) + with pytest.raises(ValueError, match=r"row index of interactions\[1\]"): + AbsorbingLS(y, x, interactions=[frame, flipped]) + # Arrays have no labels + arrays_ = Interaction(cat=cat[["h"]].to_numpy()[::-1], cont=cont[["c1"]]) + AbsorbingLS(y, x, interactions=[frame, arrays_]) + + +def test_interaction_cat_and_cont_must_match(abs_simple): + cat, cont = abs_simple.cats, abs_simple.cont + Interaction(cat=cat, cont=cont) + with pytest.raises(ValueError, match="row index of cont does not match cat"): + Interaction(cat=cat, cont=reverse(cont)) + with pytest.raises(ValueError, match="row index of cont does not match cat"): + Interaction(cat=cat, cont=relabel(cont)) + # There are no labels to compare + Interaction(cat=cat, cont=cont.to_numpy()[::-1]) + Interaction(cat=cat.to_numpy()[::-1], cont=cont) + default = pd.DataFrame(cont.to_numpy()[::-1], columns=cont.columns) + Interaction(cat=cat, cont=default) + + +def test_absorbing_numpy_is_paired_by_position(abs_simple): + y, x = abs_simple.dependent, abs_simple.exog + cat, w = abs_simple.cats, abs_simple.weights + unlabelled_cat = cat.reset_index(drop=True) + # exog as an array in the opposite order + flipped = x.to_numpy()[::-1] + res = AbsorbingLS(y, flipped, absorb=cat).fit() + raw = AbsorbingLS(y.to_numpy(), flipped, absorb=unlabelled_cat).fit() + assert_allclose(res.params, raw.params) + # weights as an array in the opposite order + weights = w.to_numpy()[::-1] + res = AbsorbingLS(y, x, absorb=cat, weights=weights).fit() + raw = AbsorbingLS( + y.to_numpy(), x.to_numpy(), absorb=unlabelled_cat, weights=weights + ).fit() + assert_allclose(res.params, raw.params) + # pandas objects with a default index + res = AbsorbingLS( + y, + pd.DataFrame(flipped, columns=x.columns), + absorb=cat, + weights=pd.Series(weights), + ).fit() + raw = AbsorbingLS(y.to_numpy(), flipped, absorb=unlabelled_cat, weights=weights) + assert_allclose(res.params, raw.fit().params) + + +def test_absorbing_same_labels_match_numpy(absorbing): + y, x, a, w = ( + absorbing.dependent, + absorbing.exog, + absorbing.cats, + absorbing.weights, + ) + res = AbsorbingLS(y, x, absorb=a, weights=w).fit() + raw = AbsorbingLS( + y.to_numpy(), + x.to_numpy(), + absorb=a.reset_index(drop=True), + weights=w.to_numpy(), + ).fit() + assert_allclose(res.params, raw.params) + assert_allclose(res.std_errors, raw.std_errors) + + +def test_absorbing_unlabelled_dependent_is_not_checked(abs_simple): + y = abs_simple.dependent.to_numpy() + AbsorbingLS( + y, + reverse(abs_simple.exog), + absorb=reverse(abs_simple.cats), + weights=reverse(abs_simple.weights), + ) + + +def test_absorbing_clusters(abs_simple): + mod = abs_simple.model() + res = mod.fit(cov_type="clustered", clusters=abs_simple.clusters) + raw = mod.fit(cov_type="clustered", clusters=abs_simple.clusters.to_numpy()) + assert_allclose(res.std_errors, raw.std_errors) + for corrupt in CORRUPTIONS.values(): + with pytest.raises(ValueError, match="row index of clusters does not match"): + mod.fit(cov_type="clustered", clusters=corrupt(abs_simple.clusters)) + values = abs_simple.clusters.to_numpy()[::-1] + res = mod.fit(cov_type="clustered", clusters=values) + raw = mod.fit(cov_type="clustered", clusters=pd.Series(values)) + assert_allclose(res.std_errors, raw.std_errors) + # The clusters are checked before the effects are absorbed + fresh = abs_simple.model() + with pytest.raises(ValueError, match="row index of clusters"): + fresh.fit(cov_type="clustered", clusters=reverse(abs_simple.clusters)) + assert fresh._absorbed_dependent is None + + +def test_absorbing_missing_values(abs_simple): + exog = abs_simple.exog.copy() + exog.iloc[4, 1] = np.nan + with pytest.warns(MissingValueWarning): + mod = AbsorbingLS(abs_simple.dependent, exog, absorb=abs_simple.cats) + kept = abs_simple.clusters[~np.asarray(exog.isnull().any(axis=1))] + assert mod.dependent.shape[0] == NOBS - 1 + # clusters are for the observations that are used + with pytest.raises(ValueError, match="row index of clusters"): + mod.fit(cov_type="clustered", clusters=reverse(kept)) + with pytest.raises(ValueError, match="row index of absorb"): + AbsorbingLS(abs_simple.dependent, exog, absorb=reverse(abs_simple.cats)) + + +def test_absorbing_xarray(abs_simple): + exog = to_xarray(abs_simple.exog) + AbsorbingLS(abs_simple.dependent, exog, absorb=abs_simple.cats) + with pytest.raises(ValueError, match="row index of exog"): + AbsorbingLS( + abs_simple.dependent, + exog.isel(obs=slice(None, None, -1)), + absorb=abs_simple.cats, + ) diff --git a/linearmodels/tests/shared/test_alignment.py b/linearmodels/tests/shared/test_alignment.py new file mode 100644 index 0000000000..48483911d8 --- /dev/null +++ b/linearmodels/tests/shared/test_alignment.py @@ -0,0 +1,217 @@ +import numpy as np +import pandas as pd +import pytest + +from linearmodels.iv.data import IVData +from linearmodels.panel.data import PanelData +from linearmodels.shared.alignment import ( + check_alignment, + check_row_alignment, + differs, + is_default_index, + is_labelled, + labels_1d, + misaligned, + row_labels, +) + +N = 6 + + +@pytest.fixture +def index(): + return pd.Index(list("abcdef"), name="obs") + + +@pytest.fixture +def series(index): + return pd.Series(np.arange(N, dtype=float), index=index, name="s") + + +@pytest.fixture +def frame(index): + return pd.DataFrame( + np.arange(2.0 * N).reshape(N, 2), index=index, columns=["a", "b"] + ) + + +@pytest.fixture +def default_frame(): + return pd.DataFrame(np.arange(2.0 * N).reshape(N, 2), columns=["a", "b"]) + + +@pytest.fixture +def multi_index(): + return pd.MultiIndex.from_product([["x", "y"], [1, 2, 3]], names=["e", "t"]) + + +def test_is_default_index(): + assert is_default_index(pd.RangeIndex(N)) + assert is_default_index(pd.Index(np.arange(N))) + assert is_default_index(pd.RangeIndex(0)) + assert is_default_index(pd.RangeIndex(N, name="named")) + assert not is_default_index(pd.RangeIndex(N)[::-1]) + assert not is_default_index(pd.RangeIndex(1, N + 1)) + assert not is_default_index(pd.Index(np.arange(N)[::-1])) + assert not is_default_index(pd.Index(list("abcdef"))) + assert not is_default_index(pd.date_range("2000-01-01", periods=N)) + assert not is_default_index(pd.MultiIndex.from_product([[0, 1], [0, 1, 2]])) + + +def test_is_labelled_pandas(series, frame, default_frame, multi_index): + assert is_labelled(series) + assert not is_labelled(series.reset_index(drop=True)) + assert is_labelled(series.iloc[::-1].reset_index(drop=True)) is False + assert is_labelled(frame) + assert not is_labelled(default_frame) + assert is_labelled(pd.Series(np.arange(N), index=multi_index)) + assert is_labelled(pd.DataFrame(np.ones((N, 2)), index=multi_index)) + # Reversed range + assert is_labelled(default_frame.iloc[::-1]) + assert is_labelled(pd.Series(np.arange(N)).iloc[::-1]) + + +def test_is_labelled_columns(default_frame): + # Columns are variables, unless they identify entities + assert not is_labelled(default_frame) + assert not is_labelled(default_frame, columns=False) + assert is_labelled(default_frame, columns=True) + unnamed = pd.DataFrame(np.ones((N, 2))) + assert not is_labelled(unnamed, columns=True) + named_rows = default_frame.set_axis(list("abcdef"), axis=0) + assert is_labelled(named_rows) + assert is_labelled(named_rows, columns=True) + + +def test_is_labelled_other_types(series, multi_index): + assert not is_labelled(None) + assert not is_labelled(np.arange(N)) + assert not is_labelled(np.ones((N, 2))) + assert not is_labelled(list(range(N))) + assert not is_labelled(1.0) + assert is_labelled(PanelData(pd.Series(np.arange(N), index=multi_index))) + assert is_labelled(IVData(series)) + assert not is_labelled(IVData(series.reset_index(drop=True))) + assert not is_labelled(IVData(series.to_numpy())) + assert not is_labelled(IVData(None, nobs=N)) + # IVData built from IVData + assert is_labelled(IVData(IVData(series))) + assert not is_labelled(IVData(IVData(series.to_numpy()))) + + +def test_is_labelled_xarray(frame): + xr = pytest.importorskip("xarray") + labelled = xr.DataArray(frame.to_numpy(), dims=["obs", "var"]) + # Any DataArray can carry labels. row_labels looks at the coordinates + assert is_labelled(labelled) + assert row_labels(labelled) is None + coords = {"obs": list(frame.index), "var": list(frame.columns)} + labelled = xr.DataArray(frame.to_numpy(), dims=["obs", "var"], coords=coords) + assert row_labels(labelled).equals(frame.index) + assert is_labelled(IVData(labelled)) + + +def test_row_labels(series, frame, default_frame, index): + assert row_labels(series).equals(index) + assert row_labels(frame).equals(index) + assert row_labels(series.reset_index(drop=True)) is None + assert row_labels(default_frame) is None + assert row_labels(None) is None + assert row_labels(series.to_numpy()) is None + assert row_labels([1, 2, 3]) is None + assert row_labels(default_frame.iloc[::-1]).equals(pd.RangeIndex(N)[::-1]) + + +def test_row_labels_ivdata(frame, index): + ivd = IVData(frame) + assert row_labels(ivd).equals(index) + assert row_labels(IVData(frame.to_numpy())) is None + assert row_labels(IVData(frame.reset_index(drop=True))) is None + # The rows that are used are the ones that are not dropped + drop = np.zeros(N, dtype=bool) + drop[[1, 4]] = True + ivd.drop(drop) + assert row_labels(ivd).equals(index[~drop]) + # An array that has had rows dropped has no labels, even though its + # remaining row numbers are not 0, 1, ... + raw = IVData(frame.to_numpy()) + raw.drop(drop) + assert row_labels(raw) is None + assert row_labels(IVData(ivd)).equals(index[~drop]) + + +def test_labels_1d(series, frame, index, multi_index): + assert labels_1d(series).equals(index) + assert labels_1d(frame).equals(index) + assert labels_1d(IVData(frame)).equals(index) + panel = PanelData(pd.Series(np.arange(N), index=multi_index)) + assert labels_1d(panel).equals(multi_index) + + +def test_labels_1d_xarray(frame, index): + xr = pytest.importorskip("xarray") + coords = {"obs": list(frame.index), "var": list(frame.columns)} + array = xr.DataArray(frame.to_numpy(), dims=["obs", "var"], coords=coords) + assert labels_1d(array).equals(index) + one_d = xr.DataArray( + frame.to_numpy()[:, 0], dims=["obs"], coords={"obs": list(frame.index)} + ) + assert labels_1d(one_d).equals(index) + + +def test_differs(index): + assert not differs(index, index) + assert not differs(index.copy(), index) + assert differs(index[::-1], index) + assert differs(pd.Index(list("abcdeg")), index) + # A default index has no labels + assert not differs(pd.RangeIndex(N), index) + assert differs(pd.RangeIndex(N)[::-1], index) + # Names are not labels + assert not differs(index.rename("another"), index) + + +def test_misaligned_message(): + exc = misaligned("exog", "The row index") + assert isinstance(exc, ValueError) + message = str(exc) + assert message.startswith("The row index of exog does not match dependent.") + assert "in the same order as dependent" in message + assert "reindex exog to the index of dependent" in message + assert ".to_numpy()" in message + exc = misaligned("weights", "The time index", "portfolios") + assert str(exc).startswith("The time index of weights does not match portfolios.") + assert "to the index of portfolios" in str(exc) + + +def test_check_row_alignment(index): + check_row_alignment("exog", index, index) + check_row_alignment("exog", index.copy(), index) + with pytest.raises(ValueError, match="row index of exog does not match dependent"): + check_row_alignment("exog", index[::-1], index) + with pytest.raises(ValueError, match="row index of x does not match reference"): + check_row_alignment("x", index[::-1], index, "reference") + # Default index and a different number of rows are not compared + check_row_alignment("exog", pd.RangeIndex(N), index) + check_row_alignment("exog", index[:-1], index) + check_row_alignment("exog", index[::-1][:-1], index) + + +def test_check_alignment(series, frame, default_frame, index): + inputs = {"exog": frame, "weights": series, "clusters": None} + check_alignment(index, inputs) + check_alignment(None, {"exog": frame.iloc[::-1]}) + check_alignment(index, {}) + check_alignment(index, {"exog": default_frame.iloc[:, :1], "w": np.ones(N)}) + with pytest.raises(ValueError, match="row index of exog does not match dependent"): + check_alignment(index, {"exog": frame.iloc[::-1]}) + with pytest.raises(ValueError, match="row index of w does not match the dependent"): + check_alignment( + index, + {"exog": frame, "w": series.iloc[::-1]}, + reference="the dependent", + ) + # The first input that does not match is reported + bad = {"a": frame, "b": frame.iloc[::-1], "c": series.iloc[::-1]} + with pytest.raises(ValueError, match="row index of b"): + check_alignment(index, bad) diff --git a/linearmodels/tests/system/test_input_alignment.py b/linearmodels/tests/system/test_input_alignment.py new file mode 100644 index 0000000000..ff700958a3 --- /dev/null +++ b/linearmodels/tests/system/test_input_alignment.py @@ -0,0 +1,606 @@ +""" +Inputs that carry labels (pandas, xarray) have to use the labels of dependent. + +The system estimators pair the rows of their inputs by position, both within +an equation and across the equations: the residuals of all equations in the +same row are assumed to be the errors of the same observation. A labelled input +whose rows are in a different order from dependent would silently be attached +to the wrong observations. These tests check that this is rejected, that inputs +that carry no labels (NumPy arrays and pandas objects with a default index) are +paired by position, and that inputs that use the labels of dependent are +accepted. +""" + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.iv.data import IVData +from linearmodels.shared.exceptions import MissingValueWarning +from linearmodels.system import IV3SLS, SUR, IVSystemGMM + +NOBS = 100 +LABELS = ["eq1", "eq2", "eq3"] +INDEX_STYLES = ["str", "datetime", "int", "multi"] +# (model, kind of equations, inputs of each equation that can be misaligned) +SUR_INPUTS = ["exog", "weights"] +IV_INPUTS = ["exog", "endog", "instruments", "weights"] +MODELS = [ + pytest.param(SUR, "sur", SUR_INPUTS, id="SUR"), + pytest.param(IV3SLS, "iv", IV_INPUTS, id="IV3SLS"), + pytest.param(IVSystemGMM, "iv", IV_INPUTS, id="IVSystemGMM"), +] +LS_MODELS = [ + pytest.param(SUR, "sur", id="SUR"), + pytest.param(IV3SLS, "iv", id="IV3SLS"), +] + + +def make_index(style, nobs=NOBS): + rng = np.random.default_rng(99) + if style == "str": + return pd.Index([f"r{i}" for i in range(nobs)], name="obs") + if style == "datetime": + return pd.date_range("2000-01-01", periods=nobs, freq="D", name="date") + if style == "int": + return pd.Index(rng.permutation(10 * nobs)[:nobs] + 7, name="id") + return pd.MultiIndex.from_product( + [np.arange(nobs // 10), np.arange(10)], names=["firm", "year"] + ) + + +class System: + """Labelled data for a system. All of the pandas objects share an index""" + + def __init__(self, style="str", nobs=NOBS, seed=4): + rng = np.random.default_rng(seed) + self.index = idx = make_index(style, nobs) + self.nobs = nobs + common = rng.standard_normal(nobs) + self.data = {} + for k, label in enumerate(LABELS): + instruments = pd.DataFrame( + rng.standard_normal((nobs, 2)), + index=idx, + columns=[f"z{k}_a", f"z{k}_b"], + ) + endog = pd.DataFrame( + {f"e{k}": instruments.iloc[:, 0] + instruments.iloc[:, 1]}, index=idx + ) + endog += rng.standard_normal((nobs, 1)) + exog = pd.DataFrame( + {"const": 1.0, f"x{k}": rng.standard_normal(nobs)}, index=idx + ) + noise = rng.standard_normal(nobs) + common + dependent = pd.Series( + 1 + exog.iloc[:, 1] + endog.iloc[:, 0] + noise, index=idx, name=f"y{k}" + ) + weights = pd.Series( + np.exp(rng.standard_normal(nobs)), index=idx, name=f"w{k}" + ) + self.data[label] = { + "dependent": dependent, + "exog": exog, + "endog": endog, + "instruments": instruments, + "weights": weights, + } + self.clusters = pd.Series(rng.integers(0, 10, nobs), index=idx, name="cl") + + def frame(self): + """A single frame with all of the variables""" + parts = [value for eq in self.data.values() for value in eq.values()] + frame = pd.concat(parts, axis=1) + return frame.loc[:, ~frame.columns.duplicated()] + + +@pytest.fixture(params=INDEX_STYLES) +def system(request): + return System(request.param) + + +@pytest.fixture +def simple(): + return System("str") + + +def reverse(value): + """The same labels in the opposite order""" + return value.iloc[::-1] + + +def relabel(value): + """The same values, with the labels attached to the wrong rows""" + out = value.copy() + out.index = value.index[np.roll(np.arange(len(value)), 1)] + return out + + +CORRUPTIONS = {"reverse": reverse, "relabel": relabel} + + +def equations(system, kind, *, weights=False, replace=None, form="dict", raw=False): + """ + Equations of a system + + kind is "sur" or "iv". replace is a dictionary with keys (label, input) and + values that replace the input in the equation. raw produces NumPy arrays. + """ + names = ["dependent", "exog"] + (["endog", "instruments"] if kind == "iv" else []) + if weights: + names.append("weights") + replace = {} if replace is None else replace + out = {} + for label in LABELS: + values = {} + for name in names: + value = replace.get((label, name), system.data[label][name]) + if raw and (label, name) not in replace: + value = value.to_numpy() + values[name] = value + out[label] = tuple(values.values()) if form == "tuple" else values + return out + + +def arrays(mod): + """Arrays used to estimate the model""" + return [*mod._y, *mod._x, *mod._z, *mod._w] + + +def assert_same_arrays(left, right): + for a, b in zip(arrays(left), arrays(right), strict=True): + assert_allclose(a, b) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +@pytest.mark.parametrize("form", ["dict", "tuple"]) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +@pytest.mark.parametrize("label", LABELS) +def test_misaligned_input_rejected(system, cls, kind, inputs, label, corruption, form): + # The model is accepted when all of the inputs use the labels of dependent + cls(equations(system, kind, weights=True, form=form)) + for name in inputs: + original = system.data[label][name] + bad = CORRUPTIONS[corruption](original) + replace = {(label, name): bad} + eqns = equations(system, kind, weights=True, replace=replace, form=form) + match = ( + f"row index of {name} in equation '{label}' does not match " + f"dependent in equation '{label}'" + ) + with pytest.raises(ValueError, match=match): + cls(eqns) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_message_names_input_and_remedy(simple, cls, kind, inputs): + replace = {("eq2", "exog"): reverse(simple.data["eq2"]["exog"])} + with pytest.raises(ValueError, match="does not match") as exc: + cls(equations(simple, kind, replace=replace)) + message = str(exc.value) + assert "exog in equation 'eq2'" in message + assert "dependent in equation 'eq2'" in message + assert "reindex" in message + assert ".to_numpy()" in message + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +@pytest.mark.parametrize("label", LABELS[1:]) +def test_equations_must_share_the_index(simple, cls, kind, inputs, label, corruption): + # Only dependent is different from the first equation + corrupt = CORRUPTIONS[corruption] + replace = {(label, "dependent"): corrupt(simple.data[label]["dependent"])} + match = ( + f"row index of dependent in equation '{label}' does not match " + "dependent in equation 'eq1'" + ) + with pytest.raises(ValueError, match=match): + cls(equations(simple, kind, replace=replace)) + # An equation where everything has the order of its own dependent variable + replace = { + (label, name): corrupt(simple.data[label][name]) + for name in ["dependent", "exog", "endog", "instruments", "weights"] + if (kind == "iv" or name in ("dependent", "exog", "weights")) + } + eqns = equations(simple, kind, weights=True, replace=replace) + with pytest.raises(ValueError, match=match): + cls(eqns) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_first_equation_is_the_reference(simple, cls, kind, inputs): + # eq1 is reversed, and the other equations do not match it + replace = { + ("eq1", name): reverse(simple.data["eq1"][name]) + for name in ["dependent", "exog", "endog", "instruments"] + if (kind == "iv" or name in ("dependent", "exog")) + } + with pytest.raises(ValueError, match="dependent in equation 'eq1'"): + cls(equations(simple, kind, replace=replace)) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +@pytest.mark.parametrize("name", ["exog", "weights"]) +def test_numpy_input_is_paired_by_position(simple, cls, kind, inputs, name): + # An array has no labels, and so it is paired by position + values = simple.data["eq2"][name].to_numpy()[::-1] + replace = {("eq2", name): values} + mod = cls(equations(simple, kind, weights=True, replace=replace)) + raw_replace = {("eq2", name): values} + raw = cls(equations(simple, kind, weights=True, replace=raw_replace, raw=True)) + assert_same_arrays(mod, raw) + # The values of eq2 are those of the opposite rows + if name == "exog": + assert_allclose(mod._exog[1].ndarray, values) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_default_index_is_paired_by_position(simple, cls, kind, inputs): + exog = simple.data["eq2"]["exog"] + flipped = pd.DataFrame(exog.to_numpy()[::-1], columns=exog.columns) + assert isinstance(flipped.index, pd.RangeIndex) + weights = pd.Series(simple.data["eq3"]["weights"].to_numpy()[::-1]) + replace = {("eq2", "exog"): flipped, ("eq3", "weights"): weights} + mod = cls(equations(simple, kind, weights=True, replace=replace)) + raw_replace = { + ("eq2", "exog"): exog.to_numpy()[::-1], + ("eq3", "weights"): weights.to_numpy(), + } + raw = cls(equations(simple, kind, weights=True, replace=raw_replace, raw=True)) + assert_same_arrays(mod, raw) + # A reversed range is not a default index + replace = {("eq2", "exog"): flipped.iloc[::-1]} + with pytest.raises(ValueError, match="row index of exog in equation 'eq2'"): + cls(equations(simple, kind, replace=replace)) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +@pytest.mark.parametrize("form", ["dict", "tuple"]) +def test_same_labels_match_numpy(system, cls, kind, inputs, form): + labelled = cls(equations(system, kind, weights=True, form=form)) + raw = cls(equations(system, kind, weights=True, form=form, raw=True)) + assert_same_arrays(labelled, raw) + + +@pytest.mark.parametrize(("cls", "kind"), LS_MODELS) +def test_same_labels_match_numpy_results(simple, cls, kind): + labelled = cls(equations(simple, kind, weights=True)).fit() + raw = cls(equations(simple, kind, weights=True, raw=True)).fit() + assert_allclose(labelled.params, raw.params) + assert_allclose(labelled.std_errors, raw.std_errors) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_everything_reordered_together_is_allowed(simple, cls, kind, inputs): + expected = cls(equations(simple, kind, weights=True)) + replace = { + (label, name): reverse(simple.data[label][name]) + for label in LABELS + for name in ["dependent", "exog", "endog", "instruments", "weights"] + if (kind == "iv" or name in ("dependent", "exog", "weights")) + } + reordered = cls(equations(simple, kind, weights=True, replace=replace)) + for a, b in zip(arrays(expected), arrays(reordered), strict=True): + assert_allclose(a, b[::-1]) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_unlabelled_dependent_is_not_checked(simple, cls, kind, inputs): + # There are no labels to compare to, so the inputs are paired by position + replace = {("eq1", "dependent"): simple.data["eq1"]["dependent"].to_numpy()} + for name in inputs: + replace[("eq1", name)] = reverse(simple.data["eq1"][name]) + mod = cls(equations(simple, kind, weights=True, replace=replace)) + assert_allclose( + mod._exog[0].ndarray, reverse(simple.data["eq1"]["exog"]).to_numpy() + ) + # The other equations are compared to their own dependent + replace = {("eq2", name): reverse(simple.data["eq2"][name]) for name in inputs} + replace[("eq2", "dependent")] = simple.data["eq2"]["dependent"].to_numpy() + cls(equations(simple, kind, weights=True, replace=replace)) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_inputs_are_compared_to_the_dependent_of_the_equation( + simple, cls, kind, inputs +): + # eq1 does not have labels. eq2 and eq3 are consistent with themselves + arrays_ = {("eq1", "dependent"): simple.data["eq1"]["dependent"].to_numpy()} + cls(equations(simple, kind, replace=arrays_)) + # eq2 is not consistent with itself + replace = {**arrays_, ("eq2", "exog"): reverse(simple.data["eq2"]["exog"])} + with pytest.raises(ValueError, match="row index of exog in equation 'eq2'"): + cls(equations(simple, kind, replace=replace)) + # eq3 dependent does not have the labels of eq2 dependent, which is the first + # dependent with labels + replace = { + **arrays_, + ("eq3", "dependent"): reverse(simple.data["eq3"]["dependent"]), + } + match = ( + "row index of dependent in equation 'eq3' does not match " + "dependent in equation 'eq2'" + ) + with pytest.raises(ValueError, match=match): + cls(equations(simple, kind, replace=replace)) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_index_names_are_ignored(simple, cls, kind, inputs): + exog = simple.data["eq2"]["exog"].rename_axis("another name") + mod = cls(equations(simple, kind, replace={("eq2", "exog"): exog})) + assert_allclose(mod._exog[1].ndarray, simple.data["eq2"]["exog"].to_numpy()) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_different_lengths_are_not_labels(simple, cls, kind, inputs): + # The number of observations in an equation is reported as before + exog = simple.data["eq2"]["exog"].iloc[:-1] + with pytest.raises(ValueError, match="Array required to have"): + cls(equations(simple, kind, replace={("eq2", "exog"): exog})) + # An equation with fewer observations than the others + names = ["dependent", "exog"] + (["endog", "instruments"] if kind == "iv" else []) + shorter = {("eq2", name): simple.data["eq2"][name].iloc[:-1] for name in names} + match = ( + "All equations must have the same number of observations. The dependent " + "variable in equation 'eq2' has 99 observations, but in equation 'eq1' " + "it has 100." + ) + with pytest.raises(ValueError, match=match): + cls(equations(simple, kind, replace=shorter)) + # The same is true without labels + values = {key: value.to_numpy() for key, value in shorter.items()} + with pytest.raises(ValueError, match=match): + cls(equations(simple, kind, replace=values)) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_missing_values_do_not_change_the_check(simple, cls, kind, inputs): + exog = simple.data["eq2"]["exog"].copy() + exog.iloc[3, 1] = np.nan + weights = simple.data["eq3"]["weights"].copy() + weights.iloc[10] = np.nan + replace = {("eq2", "exog"): exog, ("eq3", "weights"): weights} + with pytest.warns(MissingValueWarning): + mod = cls(equations(simple, kind, weights=True, replace=replace)) + assert mod._y[0].shape[0] == NOBS - 2 + replace[("eq1", "exog")] = reverse(simple.data["eq1"]["exog"]) + with pytest.raises(ValueError, match="row index of exog in equation 'eq1'"): + cls(equations(simple, kind, weights=True, replace=replace)) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_ivdata_inputs(simple, cls, kind, inputs): + exog = simple.data["eq2"]["exog"] + cls(equations(simple, kind, replace={("eq2", "exog"): IVData(exog)})) + with pytest.raises(ValueError, match="row index of exog in equation 'eq2'"): + cls(equations(simple, kind, replace={("eq2", "exog"): IVData(reverse(exog))})) + cls( + equations( + simple, kind, replace={("eq2", "exog"): IVData(exog.to_numpy()[::-1])} + ) + ) + + +@pytest.mark.parametrize(("cls", "kind", "inputs"), MODELS) +def test_xarray_inputs(simple, cls, kind, inputs): + xr = pytest.importorskip("xarray") + + def to_xarray(frame): + return xr.DataArray( + frame.to_numpy(), + dims=["obs", "var"], + coords={"obs": frame.index.to_numpy(), "var": list(frame)}, + ) + + exog = simple.data["eq2"]["exog"] + replace = {("eq2", "exog"): to_xarray(exog)} + mod = cls(equations(simple, kind, replace=replace)) + assert_allclose(mod._exog[1].ndarray, exog.to_numpy()) + flipped = to_xarray(exog).isel(obs=slice(None, None, -1)) + with pytest.raises(ValueError, match="row index of exog in equation 'eq2'"): + cls(equations(simple, kind, replace={("eq2", "exog"): flipped})) + # An xarray dependent sets the labels + dep = to_xarray(simple.data["eq1"]["dependent"].to_frame()) + replace = { + ("eq1", "dependent"): dep, + ("eq3", "exog"): reverse(simple.data["eq3"]["exog"]), + } + with pytest.raises(ValueError, match="row index of exog in equation 'eq3'"): + cls(equations(simple, kind, replace=replace)) + + +# Multivariate interfaces + + +class Multivariate: + def __init__(self, style="str"): + system = System(style) + names = LABELS + self.dependent = pd.concat( + [system.data[label]["dependent"] for label in names], axis=1 + ) + self.exog = system.data["eq1"]["exog"] + self.endog = system.data["eq1"]["endog"] + self.instruments = system.data["eq1"]["instruments"] + + +@pytest.fixture +def multi(): + return Multivariate("str") + + +def test_multivariate_ls(multi): + mod = SUR.multivariate_ls(multi.dependent, multi.exog) + raw = SUR.multivariate_ls(multi.dependent.to_numpy(), multi.exog.to_numpy()) + assert_same_arrays(mod, raw) + with pytest.raises(ValueError, match="row index of exog in equation 'y0'"): + SUR.multivariate_ls(multi.dependent, reverse(multi.exog)) + # An array has no labels + mod = SUR.multivariate_ls(multi.dependent, multi.exog.to_numpy()[::-1]) + assert_allclose(mod._exog[0].ndarray, multi.exog.to_numpy()[::-1]) + with pytest.raises(ValueError, match="row index of exog"): + SUR.multivariate_ls(reverse(multi.dependent), multi.exog) + + +def test_multivariate_iv(multi): + mod = IV3SLS.multivariate_iv( + multi.dependent, multi.exog, multi.endog, multi.instruments + ) + raw = IV3SLS.multivariate_iv( + multi.dependent.to_numpy(), + multi.exog.to_numpy(), + multi.endog.to_numpy(), + multi.instruments.to_numpy(), + ) + assert_same_arrays(mod, raw) + for name in ["exog", "endog", "instruments"]: + args = { + "exog": multi.exog, + "endog": multi.endog, + "instruments": multi.instruments, + name: reverse(getattr(multi, name)), + } + with pytest.raises(ValueError, match=f"row index of {name} in equation 'y0'"): + IV3SLS.multivariate_iv(multi.dependent, **args) + # Arrays have no labels + args[name] = getattr(multi, name).to_numpy()[::-1] + IV3SLS.multivariate_iv(multi.dependent, **args) + + +# Formulas use a single DataFrame + + +@pytest.mark.parametrize("cls", [SUR, IV3SLS, IVSystemGMM]) +def test_from_formula(simple, cls): + frame = simple.frame() + assert not isinstance(frame.index, pd.RangeIndex) + if cls is SUR: + formula = {"eq1": "y0 ~ 1 + x0", "eq2": "y1 ~ 1 + x1"} + else: + formula = { + "eq1": "y0 ~ 1 + x0 + [e0 ~ z0_a + z0_b]", + "eq2": "y1 ~ 1 + x1 + [e1 ~ z1_a + z1_b]", + } + plain = cls.from_formula(formula, frame) + weights = {"eq1": frame["w0"], "eq2": frame["w1"]} + weighted = cls.from_formula(formula, frame, weights=weights) + raw_weights = {key: value.to_numpy() for key, value in weights.items()} + raw = cls.from_formula(formula, frame, weights=raw_weights) + assert_same_arrays(weighted, raw) + assert len(plain._y) == 2 + for corrupt in CORRUPTIONS.values(): + bad = {"eq1": frame["w0"], "eq2": corrupt(frame["w1"])} + match = "row index of weights in equation 'eq2' does not match" + with pytest.raises(ValueError, match=match): + cls.from_formula(formula, frame, weights=bad) + # The same equations in the string format + if cls is SUR: + text = "{y0 ~ 1 + x0} {y1 ~ 1 + x1}" + else: + text = "{y0 ~ 1 + x0 + [e0 ~ z0_a + z0_b]} {y1 ~ 1 + x1 + [e1 ~ z1_a + z1_b]}" + string = cls.from_formula( + text, frame, weights={"y0": frame["w0"], "y1": frame["w1"]} + ) + assert len(string._y) == 2 + bad = {"y0": frame["w0"], "y1": reverse(frame["w1"])} + with pytest.raises(ValueError, match="row index of weights in equation 'y1'"): + cls.from_formula(text, frame, weights=bad) + + +@pytest.mark.parametrize(("cls", "kind"), LS_MODELS) +def test_from_formula_matches_equations(simple, cls, kind): + frame = simple.frame() + if kind == "sur": + formula = {"eq1": "y0 ~ 1 + x0", "eq2": "y1 ~ 1 + x1"} + else: + formula = { + "eq1": "y0 ~ 1 + x0 + [e0 ~ z0_a + z0_b]", + "eq2": "y1 ~ 1 + x1 + [e1 ~ z1_a + z1_b]", + } + formula_res = cls.from_formula(formula, frame).fit() + eqns = equations(simple, kind) + del eqns["eq3"] + eqn_res = cls(eqns).fit() + assert_allclose(formula_res.params, eqn_res.params) + + +# Clusters + + +@pytest.mark.parametrize(("cls", "kind"), LS_MODELS) +@pytest.mark.parametrize("corruption", list(CORRUPTIONS)) +def test_clusters_in_fit(simple, cls, kind, corruption): + mod = cls(equations(simple, kind)) + res = mod.fit(cov_type="clustered", clusters=simple.clusters) + raw = mod.fit(cov_type="clustered", clusters=simple.clusters.to_numpy()) + assert_allclose(res.std_errors, raw.std_errors) + bad = CORRUPTIONS[corruption](simple.clusters) + match = "row index of clusters does not match dependent in equation 'eq1'" + with pytest.raises(ValueError, match=match): + mod.fit(cov_type="clustered", clusters=bad) + + +@pytest.mark.parametrize(("cls", "kind"), LS_MODELS) +def test_clusters_numpy_is_paired_by_position(simple, cls, kind): + mod = cls(equations(simple, kind)) + values = simple.clusters.to_numpy()[::-1] + raw = mod.fit(cov_type="clustered", clusters=values) + # A default index has no labels + for default in (pd.Series(values), pd.DataFrame({"cl": values})): + res = mod.fit(cov_type="clustered", clusters=default) + assert_allclose(res.std_errors, raw.std_errors) + + +@pytest.mark.parametrize(("cls", "kind"), LS_MODELS) +def test_two_way_clusters(simple, cls, kind): + rng = np.random.default_rng(11) + second = pd.Series(rng.integers(0, 7, NOBS), index=simple.index, name="cl2") + clusters = pd.concat([simple.clusters, second], axis=1) + mod = cls(equations(simple, kind)) + res = mod.fit(cov_type="clustered", clusters=clusters) + raw = mod.fit(cov_type="clustered", clusters=clusters.to_numpy()) + assert_allclose(res.cov, raw.cov) + with pytest.raises(ValueError, match="row index of clusters"): + mod.fit(cov_type="clustered", clusters=reverse(clusters)) + unnamed = pd.DataFrame(clusters.to_numpy(), columns=["a", "b"]) + res = mod.fit(cov_type="clustered", clusters=unnamed) + assert_allclose(res.cov, raw.cov) + + +@pytest.mark.parametrize(("cls", "kind"), LS_MODELS) +def test_clusters_use_the_first_labelled_dependent(simple, cls, kind): + replace = {("eq1", "dependent"): simple.data["eq1"]["dependent"].to_numpy()} + mod = cls(equations(simple, kind, replace=replace)) + mod.fit(cov_type="clustered", clusters=simple.clusters) + match = "row index of clusters does not match dependent in equation 'eq2'" + with pytest.raises(ValueError, match=match): + mod.fit(cov_type="clustered", clusters=reverse(simple.clusters)) + # No labels in dependent. Clusters are paired by position + replace = { + (label, "dependent"): simple.data[label]["dependent"].to_numpy() + for label in LABELS + } + mod = cls(equations(simple, kind, replace=replace)) + mod.fit(cov_type="clustered", clusters=reverse(simple.clusters)) + + +@pytest.mark.parametrize(("cls", "kind"), LS_MODELS) +def test_clusters_and_missing_values(simple, cls, kind): + exog = simple.data["eq2"]["exog"].copy() + exog.iloc[7, 1] = np.nan + with pytest.warns(MissingValueWarning): + mod = cls(equations(simple, kind, replace={("eq2", "exog"): exog})) + mask = np.ones(NOBS, dtype=bool) + mask[7] = False + kept = simple.clusters[mask] + res = mod.fit(cov_type="clustered", clusters=kept) + raw = mod.fit(cov_type="clustered", clusters=kept.to_numpy()) + assert_allclose(res.std_errors, raw.std_errors) + with pytest.raises(ValueError, match="row index of clusters"): + mod.fit(cov_type="clustered", clusters=reverse(kept)) + # All of the clusters is a length error, as before + with pytest.raises(ValueError, match="clusters must be an ndarray"): + mod.fit(cov_type="clustered", clusters=simple.clusters)