diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index abdaa4b548..5b960984c5 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -282,13 +282,23 @@ def __init__( self._sqrt_weights = w = sqrt(self.weights.ndarray) # dependent variable self._y = self.dependent.ndarray - self._wy = self._y * w # model regressors self._x = c_[self.exog.ndarray, self.endog.ndarray] - self._wx = self._x * w # first-stage regressors self._z = c_[self.exog.ndarray, self.instruments.ndarray] - self._wz = self._z * w + if npall(w == 1.0): + # All of the weights are 1, as when weights is None or constant, + # and so scaling by the square root of the weights does not + # change the data. The weighted data is the data, and avoiding + # the copies substantially reduces the memory used when the + # data is large (GH 607). + self._wy = self._y + self._wx = self._x + self._wz = self._z + else: + self._wy = self._y * w + self._wx = self._x * w + self._wz = self._z * w self._has_constant = False self._regressor_is_exog = array( diff --git a/linearmodels/tests/iv/test_model.py b/linearmodels/tests/iv/test_model.py index 85df0410eb..445eaa53cb 100644 --- a/linearmodels/tests/iv/test_model.py +++ b/linearmodels/tests/iv/test_model.py @@ -17,6 +17,7 @@ from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE, IVLIML from linearmodels.iv.model import _OLS from linearmodels.iv.results import compare +from linearmodels.shared.exceptions import MissingValueWarning from linearmodels.shared.hypotheses import WaldTestStatistic from linearmodels.shared.utility import AttrDict @@ -508,3 +509,202 @@ def test_first_state_f(data): indiv = _OLS(endog[:, i], rhs).fit(cov_type="unadjusted") stat = indiv.wald_test(null_rest, null_value) assert_allclose(diag.iloc[i].loc["f.stat"], stat.stat / nz) + + +def test_weights_all_one_not_materialized(data): + # GH 607: weighted copies of the data are not stored when all of the + # weights are 1, which is the case when weights is None or constant + mod = IV2SLS(data.dep, data.exog, data.endog, data.instr) + assert mod._wy is mod._y + assert mod._wx is mod._x + assert mod._wz is mod._z + + weights = np.full(data.nobs, 2.5) + mod = IV2SLS(data.dep, data.exog, data.endog, data.instr, weights=weights) + assert mod._wy is mod._y + assert mod._wx is mod._x + assert mod._wz is mod._z + + +def test_weights_nonconstant_materialized(data): + rs = np.random.RandomState(12345) + weights = 0.5 + rs.random_sample(data.nobs) + mod = IV2SLS(data.dep, data.exog, data.endog, data.instr, weights=weights) + assert mod._wz is not mod._z + root_w = np.sqrt(mod.weights.ndarray) + assert_allclose(mod._wy, mod._y * root_w) + assert_allclose(mod._wx, mod._x * root_w) + assert_allclose(mod._wz, mod._z * root_w) + + +def test_weights_all_one_results_unchanged(data): + # Estimation with no weights or constant weights is identical to + # explicitly weighting by the square root of the normalized weights + res = IV2SLS(data.dep, data.exog, data.endog, data.instr).fit() + weights = np.full(data.nobs, 3.0) + res_const = IV2SLS( + data.dep, data.exog, data.endog, data.instr, weights=weights + ).fit() + assert_allclose(res_const.params, res.params, rtol=0, atol=0) + assert_allclose(res_const.cov, res.cov, rtol=0, atol=0) + + +def test_weights_float32_not_upcast_by_weighting(): + # GH 607: float32 data with float32 weights produces the same estimates + # as the float64 data, up to the precision of the float32 inputs, and + # constant float32 weights do not materialize weighted copies + rs = np.random.RandomState(12345) + n = 500 + dep = rs.standard_normal((n, 1)) + exog = rs.standard_normal((n, 2)) + endog = rs.standard_normal((n, 1)) + instr = rs.standard_normal((n, 2)) + weights = 0.5 + rs.random_sample(n) + + res64 = IV2SLS(dep, exog, endog, instr, weights=weights).fit() + res32 = IV2SLS( + dep.astype(np.float32), + exog.astype(np.float32), + endog.astype(np.float32), + instr.astype(np.float32), + weights=weights.astype(np.float32), + ).fit() + assert_allclose(res32.params, res64.params, rtol=1e-5, atol=1e-6) + + mod32 = IV2SLS( + dep.astype(np.float32), + exog.astype(np.float32), + endog.astype(np.float32), + instr.astype(np.float32), + weights=np.ones(n, dtype=np.float32), + ) + assert mod32._wz is mod32._z + + +def test_weights_float32_all_one_fit(): + # GH 607: float32 data with all-one float32 weights are fitted through + # the unmaterialized path, and the estimates agree with the estimates + # of the same data stored as float64 + rs = np.random.RandomState(12345) + n = 500 + dep = rs.standard_normal((n, 1)).astype(np.float32) + exog = rs.standard_normal((n, 2)).astype(np.float32) + endog = rs.standard_normal((n, 1)).astype(np.float32) + instr = rs.standard_normal((n, 2)).astype(np.float32) + weights = np.ones(n, dtype=np.float32) + + mod = IV2SLS(dep, exog, endog, instr, weights=weights) + assert mod._wy is mod._y + assert mod._wx is mod._x + assert mod._wz is mod._z + res = mod.fit() + res64 = IV2SLS( + dep.astype(np.float64), + exog.astype(np.float64), + endog.astype(np.float64), + instr.astype(np.float64), + ).fit() + assert_allclose(res.params, res64.params, rtol=1e-6, atol=1e-8) + assert_allclose(res.cov, res64.cov, rtol=1e-5, atol=1e-9) + + +def test_weights_integer_inputs_fit(): + # GH 607: integer data with constant integer weights are fitted through + # the unmaterialized path, and the estimates are identical to the + # estimates of the same data stored as float64 + rs = np.random.RandomState(12345) + n = 500 + dep = np.round(10 * rs.standard_normal((n, 1))).astype(np.int64) + exog = np.round(10 * rs.standard_normal((n, 2))).astype(np.int64) + endog = np.round(10 * rs.standard_normal((n, 1))).astype(np.int64) + instr = np.round(10 * rs.standard_normal((n, 2))).astype(np.int64) + weights = np.full(n, 2) + + mod = IV2SLS(dep, exog, endog, instr, weights=weights) + assert mod._wy is mod._y + assert mod._wx is mod._x + assert mod._wz is mod._z + res = mod.fit() + res64 = IV2SLS( + dep.astype(np.float64), + exog.astype(np.float64), + endog.astype(np.float64), + instr.astype(np.float64), + weights=np.full(n, 2.0), + ).fit() + assert_allclose(res.params, res64.params, rtol=0, atol=0) + assert_allclose(res.cov, res64.cov, rtol=0, atol=0) + + +def test_weights_all_one_matches_materialized_reference(data): + # GH 607: the estimates from the unmaterialized path agree with a + # reference computed directly from explicitly materialized weighted + # data, which is what the model stored before the copies were skipped + y = np.asarray(data.dep, dtype=np.float64) + x = np.c_[ + np.asarray(data.exog, dtype=np.float64), + np.asarray(data.endog, dtype=np.float64), + ] + z = np.c_[ + np.asarray(data.exog, dtype=np.float64), + np.asarray(data.instr, dtype=np.float64), + ] + n = y.shape[0] + root_w = np.ones((n, 1)) + wy, wx, wz = y * root_w, x * root_w, z * root_w + pz = wz @ pinv(wz.T @ wz) @ wz.T + x_pz_x_inv = pinv(wx.T @ pz @ wx) + params_ref = x_pz_x_inv @ (wx.T @ pz @ wy) + resids_ref = wy - wx @ params_ref + s2 = float(np.squeeze(resids_ref.T @ resids_ref)) / n + cov_ref = s2 * x_pz_x_inv + + res = IV2SLS(data.dep, data.exog, data.endog, data.instr).fit( + cov_type="unadjusted", debiased=False + ) + assert_allclose(res.params, params_ref.ravel(), rtol=1e-8, atol=1e-10) + assert_allclose( + np.asarray(res.resids).ravel(), resids_ref.ravel(), rtol=1e-8, atol=1e-10 + ) + assert_allclose(res.cov, cov_ref, rtol=1e-8, atol=1e-12) + + +def test_weights_single_non_one_materialized(data): + # GH 607: a single weight that is not 1 takes the materializing path, + # and the stored data is scaled by the square root of the weights + weights = np.ones(data.nobs) + weights[0] = 2.0 + mod = IV2SLS(data.dep, data.exog, data.endog, data.instr, weights=weights) + assert mod._wy is not mod._y + assert mod._wx is not mod._x + assert mod._wz is not mod._z + root_w = np.sqrt(mod.weights.ndarray) + assert_allclose(mod._wy, mod._y * root_w) + assert_allclose(mod._wx, mod._x * root_w) + assert_allclose(mod._wz, mod._z * root_w) + + +def test_weights_nan_weight_drops_observation(): + # GH 607: a NaN weight is a missing value, so the observation is + # dropped, as it was before the weighted copies were skipped, and the + # remaining all-one weights use the unmaterialized path + rs = np.random.RandomState(12345) + n = 200 + dep = rs.standard_normal((n, 1)) + exog = np.column_stack([np.ones(n), rs.standard_normal(n)]) + endog = rs.standard_normal((n, 1)) + instr = rs.standard_normal((n, 2)) + weights = np.ones(n) + weights[3] = np.nan + + with pytest.warns(MissingValueWarning): + mod = IV2SLS(dep, exog, endog, instr, weights=weights) + res = mod.fit() + keep = np.ones(n, dtype=bool) + keep[3] = False + res_ref = IV2SLS(dep[keep], exog[keep], endog[keep], instr[keep]).fit() + assert res.nobs == n - 1 + assert_allclose(res.params, res_ref.params, rtol=0, atol=0) + assert mod._wy is mod._y + assert mod._wx is mod._x + assert mod._wz is mod._z