Divide bequest pools by actual group populations, not birth shares

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.19.0] - 2026-07-29 12:00:00
 
 ### Added

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
 

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),