From 875b7461d9104383d7693dadfc0f51fef8c07198 Mon Sep 17 00:00:00 2001 From: marcelolafleur Date: Tue, 28 Jul 2026 10:42:05 -0400 Subject: [PATCH 1/2] Divide bequest pools by actual group populations, not birth shares MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With demographics varying across lifetime income groups (PR #1165), group population shares no longer equal lambdas: survivorship differs by group. The use_zeta=False branch of household.get_bq still divided each group's bequest pool by lambdas[j], so bequests received did not sum to bequests left (about +6.5% over-distribution in an OG-ZAF calibration with South African mortality gradients) and the steady state failed the aggregate resource constraint by about 1% of GDP. Divide by the group's actual population share instead — omega_SS[:, j] for the steady state, per-period omega sums along the time path — matching what the use_zeta=True branch already does. When demographics are common across groups the two coincide, so existing results are unchanged (all existing get_bq tests pass unmodified). Adds a conservation test with non-separable omega covering SS and TPI, by-j and pooled. Fixes #1186. --- CHANGELOG.md | 14 ++++++++++++++ ogcore/household.py | 17 +++++++++++++---- tests/test_household.py | 32 ++++++++++++++++++++++++++++++++ 3 files changed, 59 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 63f4fa4eb..da2a0d57e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,20 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Bug Fixes + +- Fixes Issue [#1186](https://github.com/PSLmodels/OG-Core/issues/1186): + with demographics that vary across lifetime income groups (PR #1165), the + `use_zeta = False` branch of `household.get_bq` divided each group's + bequest pool by its birth share (`lambdas[j]`) rather than its actual + population share, so bequests received did not sum to bequests left and + the steady-state resource constraint failed. Receipts are now divided by + the group's actual population (from `omega_SS` / `omega`), matching the + `use_zeta = True` branch. Results are unchanged when demographics are + common across groups. + ## [0.18.1] - 2026-07-22 12:00:00 ### Bug Fixes diff --git a/ogcore/household.py b/ogcore/household.py index 4b5c068d5..2b0c518a1 100644 --- a/ogcore/household.py +++ b/ogcore/household.py @@ -192,21 +192,30 @@ def get_bq(BQ, j, p, method): * utils.to_timepath_shape(BQ) ) / p.omega[:len_T, :, :] else: + # Divide each group's bequest pool by its actual population share + # rather than its birth share (lambdas): with demographic gradients + # across income groups (PR #1165), survivorship differs by group and + # the two are no longer equal. With common demographics they + # coincide, so results are unchanged in that case. if j is not None: if method == "SS": - bq = np.tile(BQ[j], p.S) / p.lambdas[j] + pop_j = p.omega_SS[:, j].sum() + bq = np.tile(BQ[j], p.S) / pop_j if method == "TPI": len_T = BQ.shape[0] + pop_j = p.omega[:len_T, :, j].sum(axis=1) bq = np.tile( - np.reshape(BQ[:, j] / p.lambdas[j], (len_T, 1)), (1, p.S) + np.reshape(BQ[:, j] / pop_j, (len_T, 1)), (1, p.S) ) else: if method == "SS": - BQ_per = BQ / np.squeeze(p.lambdas) + pop = p.omega_SS.sum(axis=0) + BQ_per = BQ / pop bq = np.tile(np.reshape(BQ_per, (1, p.J)), (p.S, 1)) if method == "TPI": len_T = BQ.shape[0] - BQ_per = BQ / p.lambdas.reshape(1, p.J) + pop = p.omega[:len_T, :, :].sum(axis=1) + BQ_per = BQ / pop bq = np.tile(np.reshape(BQ_per, (len_T, 1, p.J)), (1, p.S, 1)) return bq diff --git a/tests/test_household.py b/tests/test_household.py index 7694fd6f2..db33b51d3 100644 --- a/tests/test_household.py +++ b/tests/test_household.py @@ -179,6 +179,38 @@ def test_get_bq(BQ, j, p, method, expected): assert np.allclose(test_value, expected) +def test_get_bq_conserves_with_income_varying_demographics(): + """ + Bequests received must equal bequests left when group population + shares differ from lambdas, as they do with demographic gradients + across lifetime income groups (mortality varying by j). + """ + p = Specifications() + p.S = 3 + p.J = 2 + p.T = 3 + p.lambdas = np.array([0.6, 0.4]) + p.use_zeta = False + # non-separable omega: group 0 has relatively fewer old survivors + p.omega_SS = np.array([[0.18, 0.08], [0.15, 0.10], [0.22, 0.27]]) + assert not np.allclose(p.omega_SS.sum(axis=0), p.lambdas) + p.omega = np.tile(p.omega_SS.reshape((1, p.S, p.J)), (p.T, 1, 1)) + BQ = np.array([1.7, 3.1]) + # SS, all j + bq = household.get_bq(BQ, None, p, "SS") + received = (bq * p.omega_SS).sum() + assert np.allclose(received, BQ.sum()) + # SS, each j + for j in range(p.J): + bq_j = household.get_bq(BQ, j, p, "SS") + assert np.allclose((bq_j * p.omega_SS[:, j]).sum(), BQ[j]) + # TPI, all j + BQ_path = np.tile(BQ.reshape((1, p.J)), (p.T, 1)) + bq_path = household.get_bq(BQ_path, None, p, "TPI") + received_path = (bq_path * p.omega[: p.T]).sum(axis=(1, 2)) + assert np.allclose(received_path, BQ_path.sum(axis=1)) + + p1 = Specifications() p1.eta = np.tile( np.array([[0.1, 0.3], [0.15, 0.4], [0.05, 0.0]]).reshape(1, p2.S, p2.J), From 621b97b7d6d5406290e8826a7e63a0205d84eb74 Mon Sep 17 00:00:00 2001 From: marcelolafleur Date: Tue, 28 Jul 2026 16:36:28 -0400 Subject: [PATCH 2/2] Format code blocks in tests/BENCHMARK_README.md so ruff format --check passes --- tests/BENCHMARK_README.md | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/tests/BENCHMARK_README.md b/tests/BENCHMARK_README.md index 892f684fb..727b8b562 100644 --- a/tests/BENCHMARK_README.md +++ b/tests/BENCHMARK_README.md @@ -87,17 +87,17 @@ Each benchmark produces a `BenchmarkResult` with the following metrics: ```python @dataclass class BenchmarkResult: - test_name: str # Name of the test - platform: str # Operating system - scheduler: str # Dask scheduler used - num_workers: int # Number of workers - compute_time: float # Execution time in seconds - peak_memory_mb: float # Peak memory usage in MB - avg_memory_mb: float # Average memory usage in MB - data_size_mb: float # Input data size in MB - num_tasks: int # Number of parallel tasks - success: bool # Whether test succeeded - error_message: str # Error details if failed + test_name: str # Name of the test + platform: str # Operating system + scheduler: str # Dask scheduler used + num_workers: int # Number of workers + compute_time: float # Execution time in seconds + peak_memory_mb: float # Peak memory usage in MB + avg_memory_mb: float # Average memory usage in MB + data_size_mb: float # Input data size in MB + num_tasks: int # Number of parallel tasks + success: bool # Whether test succeeded + error_message: str # Error details if failed ``` ### Key Metrics to Monitor @@ -228,7 +228,7 @@ cluster = LocalCluster( n_workers=num_workers, threads_per_worker=2, processes=False, # Use threads, not processes - memory_limit='4GB', + memory_limit="4GB", ) client = Client(cluster) ```