diff --git a/CHANGELOG.md b/CHANGELOG.md index 69d9fb503..e60d69cc5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,55 @@ 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). +## [0.17.0] - 2026-07-16 12:00:00 + +### Bug Fixes + +- Fixes the Defined Benefits and Points System pension paths so they run + against a real `Specifications` object (Issues #1014 and #1075): + `p.retire` (an array since PR #433) was passed as the scalar `S_ret` + into the numba loops, the scalar `p.g_y` was indexed as an array inside + the loops, the scalar steady-state wage was indexed as a path, and the + time-varying 3-D `e` matrix (PR #895) was sliced as 2-D. The systems + use the steady-state retirement age and earnings profile for now; full + time variation remains open in Issue #1014. The same coercions are + applied to the NDC path, but the NDC system additionally requires + growth-rate settings (`ndc_growth_rate`, `dir_growth_rate`) that are + not yet parameters in `default_parameters.json`, so it still cannot + run against a real `Specifications` object (see Issue #1169). +- Wires the pre-time-path inputs the DB/NDC/PS benefit formulas need + into the TPI: labor supplied before the time path begins comes from + the model's initial condition (the same baseline object that + initializes wealth), and pre-time-path wages are anchored to the + period-0 wage of the current path. Adds the full time-path (T x S x J) + evaluation of Defined Benefits amounts (`DB_3dim_loop`) used when the + TPI computes aggregate revenues, built from each cohort's own wage and + labor history so aggregates are consistent with household behavior. + With these changes a country model using the Defined Benefits system + solves both the steady state and the transition path. +- Adds regression tests that call `pension_amount` with a real + `Specifications` object per pension system (the existing tests + pre-scalarized the inputs and so never exercised the real interface) + and a local-marked steady-state solve test with the Defined Benefits + system. See PR [#1167](https://github.com/PSLmodels/OG-Core/pull/1167). +- Fixes the tax-liability revenue calculation using the first year's tax + noncompliance and filer rates for every period, even when those rates + are set to vary over time (Issue #1168): households responded to the + changing rates while government revenue, the budget, and debt stayed on + the first-year values. `income_tax_liab` now slices the rates over the + transition path. See PR [#1174](https://github.com/PSLmodels/OG-Core/pull/1174). +- Fixes the steady-state `etr_ss` and `mtry_ss` diagnostics using the + labor-income tax noncompliance rate in place of the capital-income rate + (`capital_noncompliance_rate_2D` in `SS.py`, Issue #1170). This affects + only the post-solve SS diagnostics -- the solution itself already used + the correct rates -- and is a no-op when the two rates are equal (the + default). See PR [#1171](https://github.com/PSLmodels/OG-Core/pull/1171). +- Fixes `alpha_FA` (direct foreign aid, added in 0.16.0) not being + extended over the model time path: it was never registered in + `tp_param_list`, so a multi-year path stayed short and broke the + transition solve. It now extends to `T + S` with the last value carried + forward, like `alpha_G`/`alpha_T`/`alpha_I`. See PR [#1166](https://github.com/PSLmodels/OG-Core/pull/1166). + ## [0.16.4] - 2026-07-02 12:00:00 ### Added diff --git a/ogcore/TPI.py b/ogcore/TPI.py index 4dfed4d24..faec9be5d 100644 --- a/ogcore/TPI.py +++ b/ogcore/TPI.py @@ -164,6 +164,13 @@ def get_initial_SS_values(p): B0 = aggr.get_B(ss_baseline_vars["b_sp1"], p, "SS", True) initial_b = ss_baseline_vars["b_sp1"] * (ss_baseline_vars["B"] / B0) initial_n = ss_baseline_vars["n"] + # The DB/NDC/PS pension formulas need the labor supplied before the + # time path begins by cohorts alive at t=0. Use the model's initial + # labor condition (the same baseline object that initializes wealth); + # pre-time-path wages are anchored to the period-0 wage of the + # current path inside pensions.py. See Issue #1014 for a fully + # history-consistent treatment. + p.n_preTP = initial_n Ybaseline = None TRbaseline = None diff --git a/ogcore/__init__.py b/ogcore/__init__.py index 5a471f5e0..abe887769 100644 --- a/ogcore/__init__.py +++ b/ogcore/__init__.py @@ -21,4 +21,4 @@ from ogcore.txfunc import * # noqa: F403 from ogcore.utils import * # noqa: F403 -__version__ = "0.16.4" +__version__ = "0.17.0" diff --git a/ogcore/pensions.py b/ogcore/pensions.py index 0727ba8c2..5e2d9acc0 100644 --- a/ogcore/pensions.py +++ b/ogcore/pensions.py @@ -213,7 +213,10 @@ def DB_amount(w, e, n, j, p): Args: w (array_like): real wage rate e (Numpy array): effective labor units - n (Numpy array): labor supply + n (Numpy array): labor supply; may be a partial lifetime path + (length < S, for cohorts alive when the time path begins), + a steady-state vector (S,) or matrix (S, J), or the full + time path (T, S, J) j (int): index of lifetime income group p (OG-Core Specifications object): model parameters @@ -226,23 +229,38 @@ def DB_amount(w, e, n, j, p): equiv_periods = int(round((p.S / 80.0) * p.avg_earn_num_years)) - 1 equiv_yr_contrib = int(round((p.S / 80.0) * p.yr_contrib)) - 1 L_inc_avg_s = np.zeros(equiv_periods) + # Retirement age is potentially time varying (PR #433), but the DB + # system supports a single retirement age for now; use the + # steady-state value (see Issue #1014). g_y is a scalar parameter, + # but the numba loops index it, so pass it as a 1-element array. + S_ret = int(np.asarray(p.retire).flat[-1]) + g_y_arr = np.atleast_1d(np.asarray(p.g_y, dtype=float)).ravel() + # The wage may arrive as a scalar (steady state, or one remaining + # period in the TPI); the loops index it as a path, so broadcast it + # to the length of the labor-supply vector. + if np.ndim(w) == 0: + w = np.full(n.shape[0], float(w)) if n.shape[0] < p.S: per_rmn = n.shape[0] - # TODO: think about how to handle setting w_preTP and n_preTP - # TODO: will need to update how the e matrix is handled here - # and else where to allow for it to be time varying - w_S = np.append((p.w_preTP * np.ones(p.S))[:(-per_rmn)], w) + # Pre-time-path wages: the recent past is anchored to the + # period-0 wage of the current path (trend growth is handled + # by the de-trending inside the benefit formula); labor comes + # from the model's initial condition n_preTP (see Issue #1014) + w_S = np.append((w[0] * np.ones(p.S))[:(-per_rmn)], w) n_S = np.append(p.n_preTP[:(-per_rmn), j], n) DB = np.zeros(p.S) DB = DB_1dim_loop( w_S, - p.e[:, j], + # TODO: will need to update how the e matrix is handled + # here and elsewhere to allow for it to be time varying + # (see Issue #1014); use the steady-state profile for now + p.e[-1, :, j] if np.ndim(p.e) == 3 else p.e[:, j], n_S, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, L_inc_avg_s, L_inc_avg, DB, @@ -259,9 +277,9 @@ def DB_amount(w, e, n, j, p): w, e, n, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, L_inc_avg_s, L_inc_avg, DB, @@ -277,9 +295,9 @@ def DB_amount(w, e, n, j, p): w, e, n, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, L_inc_avg_sj, L_inc_avg, DB, @@ -288,6 +306,26 @@ def DB_amount(w, e, n, j, p): equiv_yr_contrib, ) + elif np.ndim(n) == 3: + T = n.shape[0] + w_path = np.squeeze(np.asarray(w)) + if w_path.ndim == 0: + w_path = np.full(T, float(w_path)) + e_ss = p.e[-1] if np.ndim(p.e) == 3 else p.e + DB = DB_3dim_loop( + w_path, + e_ss, + n, + p.n_preTP, + S_ret, + p.S, + p.J, + float(g_y_arr[-1]), + equiv_periods, + p.alpha_db, + equiv_yr_contrib, + ) + return DB @@ -315,22 +353,34 @@ def NDC_amount(w, e, n, r, Y, j, p): """ g_ndc_amount = g_ndc(r, Y, p) delta_ret_amount = delta_ret(r, Y, p) + # Single retirement age for now (see Issue #1014); g_y is a scalar + # but the numba loops index it; the steady-state wage is a scalar + # but the loops index it as a path. + S_ret = int(np.asarray(p.retire).flat[-1]) + g_y_arr = np.atleast_1d(np.asarray(p.g_y, dtype=float)).ravel() + if np.ndim(w) == 0: + w = np.full(n.shape[0], float(w)) if n.shape[0] < p.S: per_rmn = n.shape[0] - w_S = np.append((p.w_preTP * np.ones(p.S))[:(-per_rmn)], w) + # Pre-time-path wages: the recent past is anchored to the + # period-0 wage of the current path (trend growth is handled + # by the de-trending inside the benefit formula); labor comes + # from the model's initial condition n_preTP (see Issue #1014) + w_S = np.append((w[0] * np.ones(p.S))[:(-per_rmn)], w) n_S = np.append(p.n_preTP[:(-per_rmn), j], n) - NDC_s = np.zeros(p.retire) + NDC_s = np.zeros(S_ret) NDC = np.zeros(p.S) NDC = NDC_1dim_loop( w_S, - p.e[:, j], + # steady-state earnings profile for now (see Issue #1014) + p.e[-1, :, j] if np.ndim(p.e) == 3 else p.e[:, j], n_S, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, p.tau_p, g_ndc_amount, delta_ret_amount, @@ -341,15 +391,15 @@ def NDC_amount(w, e, n, r, Y, j, p): else: if np.ndim(n) == 1: - NDC_s = np.zeros(p.retire) + NDC_s = np.zeros(S_ret) NDC = np.zeros(p.S) NDC = NDC_1dim_loop( w, e, n, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, p.tau_p, g_ndc_amount, delta_ret_amount, @@ -357,15 +407,15 @@ def NDC_amount(w, e, n, r, Y, j, p): NDC, ) elif np.ndim(n) == 2: - NDC_sj = np.zeros((p.retire, p.J)) + NDC_sj = np.zeros((S_ret, p.J)) NDC = np.zeros((p.S, p.J)) NDC = NDC_2dim_loop( w, e, n, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, p.tau_p, g_ndc_amount, delta_ret_amount, @@ -397,19 +447,32 @@ def PS_amount(w, e, n, j, factor, p): PS (Numpy array): pension amount for each household """ + # Single retirement age for now (see Issue #1014); g_y is + # a scalar but the numba loops index it; the steady-state wage is a + # scalar but the loops index it as a path. + S_ret = int(np.asarray(p.retire).flat[-1]) + g_y_arr = np.atleast_1d(np.asarray(p.g_y, dtype=float)).ravel() + if np.ndim(w) == 0: + w = np.full(n.shape[0], float(w)) + if n.shape[0] < p.S: per_rmn = n.shape[0] - w_S = np.append((p.w_preTP * np.ones(p.S))[:(-per_rmn)], w) + # Pre-time-path wages: the recent past is anchored to the + # period-0 wage of the current path (trend growth is handled + # by the de-trending inside the benefit formula); labor comes + # from the model's initial condition n_preTP (see Issue #1014) + w_S = np.append((w[0] * np.ones(p.S))[:(-per_rmn)], w) n_S = np.append(p.n_preTP[:(-per_rmn), j], n) - L_inc_avg_s = np.zeros(p.retire) + L_inc_avg_s = np.zeros(S_ret) PS = np.zeros(p.S) PS = PS_1dim_loop( w_S, - p.e[:, j], + # steady-state earnings profile for now (see Issue #1014) + p.e[-1, :, j] if np.ndim(p.e) == 3 else p.e[:, j], n_S, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, p.vpoint, factor, L_inc_avg_s, @@ -419,15 +482,15 @@ def PS_amount(w, e, n, j, factor, p): else: if np.ndim(n) == 1: - L_inc_avg_s = np.zeros(p.retire) + L_inc_avg_s = np.zeros(S_ret) PS = np.zeros(p.S) PS = PS_1dim_loop( w, e, n, - p.retire, + S_ret, p.S, - p.g_y, + g_y_arr, p.vpoint, factor, L_inc_avg_s, @@ -435,16 +498,16 @@ def PS_amount(w, e, n, j, factor, p): ) elif np.ndim(n) == 2: - L_inc_avg_sj = np.zeros((p.retire, p.J)) + L_inc_avg_sj = np.zeros((S_ret, p.J)) PS = np.zeros((p.S, p.J)) PS = PS_2dim_loop( w, e, n, - p.retire, + S_ret, p.S, p.J, - p.g_y, + g_y_arr, p.vpoint, factor, L_inc_avg_sj, @@ -513,9 +576,11 @@ def deriv_NDC(r, w, e, Y, per_rmn, p): d_theta (Numpy array): change in NDC pension benefits for another unit of labor supply """ + # Single retirement age for now (see Issue #1014) + S_ret = int(np.asarray(p.retire).flat[-1]) if per_rmn == 1: d_theta = 0 - elif per_rmn < (p.S - p.retire + 1): + elif per_rmn < (p.S - S_ret + 1): d_theta = np.zeros(per_rmn) else: d_theta_empty = np.zeros(per_rmn) @@ -526,7 +591,7 @@ def deriv_NDC(r, w, e, Y, per_rmn, p): e, per_rmn, p.S, - p.retire, + S_ret, p.tau_p, g_ndc_amount, delta_ret_amount, @@ -561,14 +626,16 @@ def deriv_DB(w, e, per_rmn, p): """ equiv_periods = int(round((p.S / 80.0) * p.avg_earn_num_years)) - 1 equiv_yr_contrib = int(round((p.S / 80.0) * p.yr_contrib)) - 1 - if per_rmn < (p.S - p.retire + 1): + # Single retirement age for now (see Issue #1014) + S_ret = int(np.asarray(p.retire).flat[-1]) + if per_rmn < (p.S - S_ret + 1): d_theta = np.zeros(p.S) else: d_theta = deriv_DB_loop( w, e, p.S, - p.retire, + S_ret, per_rmn, equiv_periods, p.alpha_db, @@ -602,12 +669,14 @@ def deriv_PS(w, e, per_rmn, factor, p): """ - if per_rmn < (p.S - p.retire + 1): + # Single retirement age for now (see Issue #1014) + S_ret = int(np.asarray(p.retire).flat[-1]) + if per_rmn < (p.S - S_ret + 1): d_theta = np.zeros(p.S) else: d_theta_empty = np.zeros(p.S) d_theta = deriv_PS_loop( - w, e, p.S, p.retire, per_rmn, d_theta_empty, p.vpoint, factor + w, e, p.S, S_ret, per_rmn, d_theta_empty, p.vpoint, factor ) d_theta = d_theta[-per_rmn:] @@ -718,10 +787,12 @@ def delta_ret(r, Y, p): """ surv_rates = 1 - p.mort_rates_SS - dir_delta_s_empty = np.zeros(p.S - p.retire + 1) + # Single retirement age for now (see Issue #1014) + S_ret = int(np.asarray(p.retire).flat[-1]) + dir_delta_s_empty = np.zeros(p.S - S_ret + 1) g_dir_value = g_dir(r, Y, p.g_y, p.g_n, p.dir_growth_rate) dir_delta = delta_ret_loop( - p.S, p.retire, surv_rates, g_dir_value, dir_delta_s_empty + p.S, S_ret, surv_rates, g_dir_value, dir_delta_s_empty ) delta_ret = 1 / (dir_delta + p.indR - p.k_ret) @@ -1027,6 +1098,73 @@ def DB_2dim_loop( @numba.jit(nopython=True) +def DB_3dim_loop( + w_path, + e_ss, + n, + n_preTP, + S_ret, + S, + J, + g_y, + avg_earn_num_years, + alpha_db, + yr_contr, +): + r""" + Calculate public pension from a defined benefits system over the + full time path. + + Used when the TPI solution evaluates taxes for all periods at once + to compute aggregates. Each retiree's benefit is computed from their + own cohort's wage and labor history: the wage entering the average + for a household aged u at time t, earned at age s, is the wage at + time t - (u - s). Histories that predate the time path use the + period-0 wage (trend growth is handled by the de-trending in the + benefit formula) and the model's initial labor supply, so this + reproduces exactly what the per-cohort household solves compute and + household behavior and aggregates are consistent. The inner work is + vectorized rather than looped. + + Args: + w_path (Numpy array): real wage rate path, length T + e_ss (Numpy array): effective labor units, size SxJ + n (Numpy array): labor supply, size TxSxJ + n_preTP (Numpy array): pre-time-path labor supply, size SxJ + S_ret (int): retirement age + S (int): number of periods in the model + J (int): number of lifetime income groups + g_y (scalar): growth rate of technology + avg_earn_num_years (int): number of years earnings are averaged + over + alpha_db (scalar): replacement rate per year of contribution + yr_contr (int): years of contribution + + Returns: + DB (Numpy array): pension amount for each household, size TxSxJ + """ + T = n.shape[0] + s_idx = np.arange(S_ret - avg_earn_num_years, S_ret) + DB = np.zeros((T, S, J)) + for t in range(T): + for u in range(S_ret, S): + tau = t - (u - s_idx) + w_hist = w_path[np.clip(tau, 0, T - 1)] + n_hist = np.where( + (tau >= 0)[:, None], + n[np.clip(tau, 0, T - 1), s_idx, :], + n_preTP[s_idx, :], + ) + L = (w_hist[:, None] / np.exp(g_y * (u - s_idx))[:, None]) * ( + e_ss[s_idx, :] * n_hist + ) + DB[t, u, :] = ( + (L.sum(axis=0) / avg_earn_num_years) * yr_contr * alpha_db + ) + + return DB + + def NDC_1dim_loop(w, e, n, S_ret, S, g_y, tau_p, g_ndc, delta_ret, NDC_s, NDC): """ Calculate public pension from a notional defined contribution diff --git a/pyproject.toml b/pyproject.toml index 7e3a57e1f..728ed3f05 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "ogcore" -version = "0.16.4" +version = "0.17.0" authors = [ {name = "Jason DeBacker and Richard W. Evans"}, ] diff --git a/tests/test_pensions.py b/tests/test_pensions.py index d7b3302aa..5c675a485 100644 --- a/tests/test_pensions.py +++ b/tests/test_pensions.py @@ -817,9 +817,13 @@ def test_deriv_theta(args, d_theta_expected): p2.k_ret = 0.4615 p2.mort_rates_SS = np.array([0.01, 0.05, 0.3, 0.4, 1]) p2.e = e -NDC_expected2 = np.array([0, 0.25185784, 0.24441432]) +# Pre-time-path wages are anchored to the period-0 wage of the passed +# path (w[0] = 1.1), not the old w_preTP attribute. Earlier expected +# values under the old convention, kept for reference: +# NDC_expected2 = np.array([0, 0.25185784, 0.24441432]) # TODO: why move from numbers below to those above ? Diff in numpy rounding?? # NDC_expected2 = np.array([0, 0.251721214, 0.244281728]) +NDC_expected2 = np.array([0, 0.25705613, 0.24945897]) args2 = (w, e, n, r, None, j, p2) test_data = [(args1, NDC_expected1), (args2, NDC_expected2)] @@ -935,7 +939,8 @@ def test_PS_1dim_loop(args, PS_loop_expected): ] ) p2.e = e2 -PS_expected2 = np.array([0, 0, 0.003585952, 0.003479971, 0.003377123]) +# Pre-time-path wages anchored to the period-0 wage (w2[0] = 1.21) +PS_expected2 = np.array([0, 0, 0.003852672, 0.003738809, 0.00362831]) args2 = (w2, e2, n2, j, factor, p2) test_data = [(args1, PS_expected1), (args2, PS_expected2)] @@ -952,3 +957,76 @@ def test_get_PS(args, PS_expected): PS = pensions.PS_amount(w, e, n, j, factor, p) print("PS inside of the test", PS) assert np.allclose(PS, PS_expected) + + +# Regression tests for the array/scalar bugs hit when the non-US pension +# systems run against a real Specifications object (Issues #1014 and #1075): +# p.retire is an array (time-varying since PR #433) but was passed as the +# scalar S_ret into the numba loops; p.g_y is a scalar but the loops index +# it; and in the SS the wage is a scalar but the loops index it as a path. +# The mock-parameter tests above pre-scalarize these inputs, so they never +# exercised the real interface. +@pytest.mark.parametrize( + "system,updates", + [ + ( + "Defined Benefits", + {"alpha_db": 0.02, "yr_contrib": 35, "avg_earn_num_years": 40}, + ), + ("Points System", {"vpoint": 0.5}), + ], + ids=["DB", "PS"], +) +# The Notional Defined Contribution system is excluded above: its +# growth-rate settings (ndc_growth_rate, dir_growth_rate) are not yet +# parameters in default_parameters.json, so it cannot run against a real +# Specifications object at all. See Issue #1169. +def test_pension_amount_with_real_specifications(system, updates): + """ + pension_amount must accept a real Specifications object in the SS: + retire is an array, g_y a scalar, and the steady-state wage a scalar. + """ + p = Specifications() + p.update_specifications(dict(updates, pension_system=system)) + j = 0 + w = 1.2 # steady-state wage is a scalar + r, Y, factor = 0.05, 1.0, 100000.0 + n = 0.4 * np.ones(p.S) + e = p.e[-1, :, j] + theta = np.zeros(p.J) + pension = pensions.pension_amount( + r, w, n, Y, theta, None, j, False, "SS", e, factor, p + ) + S_ret = int(p.retire[-1]) + assert pension.shape == (p.S,) + assert np.all(pension[:S_ret] == 0.0) + assert np.all(pension[S_ret:] > 0.0) + + +@pytest.mark.local +def test_SS_solve_defined_benefits(tmp_path): + """ + The steady state must solve with the Defined Benefits pension system + and produce positive aggregate pension outlays (Issue #1014 asked for + model-run tests of each pension system; the NDC system cannot run yet + -- see above -- and the analogous Points System run should be added + with the resolution of that issue). + """ + from ogcore.execute import runner + from ogcore import utils + + p = Specifications(baseline=True, num_workers=1) + p.update_specifications( + { + "pension_system": "Defined Benefits", + "alpha_db": 0.02, + "yr_contrib": 35, + "avg_earn_num_years": 40, + } + ) + p.baseline_dir = p.output_base = str(tmp_path) + runner(p, time_path=False, client=None) + ss = utils.safe_read_pickle(str(tmp_path / "SS" / "SS_vars.pkl")) + Y = np.asarray(ss["Y"]).sum() + assert ss["agg_pension_outlays"] > 0.0 + assert 0.0 < ss["agg_pension_outlays"] / Y < 1.0