diff --git a/CHANGELOG.md b/CHANGELOG.md index 6adeb3f13..fae068371 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,18 @@ 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.21.0] - 2026-10-07 + +### Added +- Allows the defined benefit pension replacement rate `alpha_db` to vary across lifetime income groups and over time. The parameter is now a `T x J` array (`number_dims` changed from 0 to 2 in `default_parameters.json`), extrapolated like the other `T x J` parameters, so a single value (`[[0.02]]`), one value per income group, or a time path can be entered. The rate applied to a benefit is the one in effect in the period in which the benefit is paid, so a change in the rate affects current as well as future retirees. Reform JSON files or scripts that previously set `alpha_db` as a scalar (e.g., `0.02`) must now use the nested list form (e.g., `[[0.02]]`). +- Adds tests in `test_pensions.py` of DB benefits that vary by group in the steady state, of an incomplete lifetime straddling a change in the rate along the time path, of the derivative of benefits with a time-varying rate, and of consistency between the cohort-by-cohort and full-time-path computations of DB benefits with a rate that varies over time and by group. + +### Updated +- `pensions.DB_amount` now takes `t` and `method` arguments (`DB_amount(w, e, n, t, j, method, p)`) so that it can select the replacement rate for the period in which each benefit is paid. `pensions.DB_1dim_loop`, `DB_2dim_loop`, and `DB_3dim_loop` take the replacement rate as an array indexed by age (or by period for the full time path) rather than a scalar. +- `pensions.deriv_DB` and `pensions.deriv_theta` now take `t` and `j` arguments. With a time-varying rate, the derivative of benefits with respect to labor supply at each working age is evaluated at the rate in effect when the household retires. `deriv_theta` is not yet used in the household first order conditions, so this has no effect on model solutions. +- Fixes the title and description of `alpha_db` in `default_parameters.json`, which referred to a defined contribution system. +- Updates the defined benefit pension equations in the theory documentation (`government.md`) to use $\alpha^{DB}_{j,t}$. + ## [0.20.4] - 2026-10-05 ### Added @@ -735,6 +747,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Any earlier versions of OG-USA can be found in the [`OG-Core`](https://github.com/PSLmodels/OG-Core) repository [release history](https://github.com/PSLmodels/OG-Core/releases) from [v.0.6.4](https://github.com/PSLmodels/OG-Core/releases/tag/v0.6.4) (Jul. 20, 2021) or earlier. +[0.21.0]: https://github.com/PSLmodels/OG-Core/compare/v0.20.4...v0.21.0 [0.20.4]: https://github.com/PSLmodels/OG-Core/compare/v0.20.3...v0.20.4 [0.20.3]: https://github.com/PSLmodels/OG-Core/compare/v0.20.2...v0.20.3 [0.20.2]: https://github.com/PSLmodels/OG-Core/compare/v0.20.1...v0.20.2 diff --git a/docs/book/content/intro/parameters.md b/docs/book/content/intro/parameters.md index 179919b10..459be08e8 100644 --- a/docs/book/content/intro/parameters.md +++ b/docs/book/content/intro/parameters.md @@ -341,7 +341,8 @@ _Valid Choices:_['r', 'Curr GDP', 'LR GDP'] #### `alpha_db` -_Description:_ Replacement rate under a defined contribution system. +_Description:_ Replacement rate per year of contribution under a defined benefits system. Set value for base year, click '+' to add value for next year. All future years not specified are set to last value entered. +_Notes:_ Can be specified to vary by type J and over year T. The rate applied is the one in effect in the period the benefit is paid. _Value Type:_ float _Valid Range:_ min = 0.0 and max = 1.0 _Out-of-Range Action:_ error diff --git a/docs/book/content/theory/government.md b/docs/book/content/theory/government.md index e58face16..a7b296666 100644 --- a/docs/book/content/theory/government.md +++ b/docs/book/content/theory/government.md @@ -402,13 +402,13 @@ The defined benefit system pension amount is given as: ```{math} :label: eqn:db_pension - pension_{j,s,t} = \biggl[\frac{\sum_{s=R-ny}^{R-1}w_{t}e_{j,s,t}n_{j,s,t}}{ny}\biggr]\times Cy \times \alpha_{DB} \quad \forall s > R + pension_{j,s,t} = \biggl[\frac{\sum_{s=R-ny}^{R-1}w_{t}e_{j,s,t}n_{j,s,t}}{ny}\biggr]\times Cy \times \alpha^{DB}_{j,t} \quad \forall s > R ``` where: - $ny$ are the number of years over which average earnings are calculated - $Cy$ are the number of years of contributions. In our model, there is no exit from the labor force, so workers will contribute for $R$ years, but $Cy$ could be some number less than $R$ if there is a maximum number of years of contributions one can accrue under the DB system. - - $\alpha_{DB}$ is the replacement rate per year of contribution. + - $\alpha^{DB}_{j,t}$ is the replacement rate per year of contribution. This rate may vary over time and across lifetime income groups. The rate applied to a benefit is the one in effect in the period in which the benefit is paid, so a change in the rate affects the benefits of current retirees as well as those of future retirees. Given this pension system and the fact that there is only variation in labor supply along the intensive margin (so we don't need to consider changes in $Cy$), the partial derivatives from the household section are given by: @@ -417,7 +417,7 @@ where: \frac{\partial \theta_{j,u,t+u-s}}{\partial n_{j,s,t}} = \begin{cases} 0 , & \text{if}\ s < R - Cy \\ - w_{t}e_{j,s}\alpha_{DB}\times \frac{Cy}{ny}, & \text{if}\ R - Cy <= s < R \\ + w_{t}e_{j,s}\alpha^{DB}_{j,t+R-s}\times \frac{Cy}{ny}, & \text{if}\ R - Cy <= s < R \\ 0, & \text{if}\ s \geq R \\ \end{cases} ``` diff --git a/ogcore/__init__.py b/ogcore/__init__.py index dcd045ef1..a272fd99f 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.20.4" +__version__ = "0.21.0" diff --git a/ogcore/default_parameters.json b/ogcore/default_parameters.json index c1bffd247..c88dc890d 100644 --- a/ogcore/default_parameters.json +++ b/ogcore/default_parameters.json @@ -4145,15 +4145,19 @@ } }, "alpha_db": { - "title": "Replacement rate under a defined contribution system.", - "description": "Replacement rate under a defined contribution system.", + "title": "Replacement rate under a defined benefits system.", + "description": "Replacement rate per year of contribution under a defined benefits system. Set value for base year, click '+' to add value for next year. All future years not specified are set to last value entered.", + "param_notation": "$\\alpha^{DB}_{j,t}$", "section_1": "Fiscal Policy Parameters", "section_2": "Government Pension Parameters", - "notes": "", + "notes": "Can be specified to vary by type J and over year T. The rate applied is the one in effect in the period the benefit is paid.", "type": "float", + "number_dims": 2, "value": [ { - "value": 0.0 + "value": [[ + 0.0 + ]] } ], "validators": { diff --git a/ogcore/parameters.py b/ogcore/parameters.py index a6c14cd2f..b9325c45a 100644 --- a/ogcore/parameters.py +++ b/ogcore/parameters.py @@ -233,6 +233,7 @@ def compute_default_params(self): "income_tax_filer", "wealth_tax_filer", "tau_payroll", + "alpha_db", ] for item in tp_param_list3: param_in = np.asarray(getattr(self, item)) diff --git a/ogcore/pensions.py b/ogcore/pensions.py index 6102e000a..881996a94 100644 --- a/ogcore/pensions.py +++ b/ogcore/pensions.py @@ -140,7 +140,7 @@ def pension_amount(r, w, n, Y, theta, t, j, shift, method, e, factor, p): if p.pension_system == "US-Style Social Security": pension = SS_amount(w, n, theta, t, j, shift, method, e, p) elif p.pension_system == "Defined Benefits": - pension = DB_amount(w, e, n, j, p) + pension = DB_amount(w, e, n, t, j, method, p) pension = pension * replacement_rate_adjustment( pension, t, j, method, p ) @@ -249,15 +249,19 @@ def SS_amount(w, n, theta, t, j, shift, method, e, p): return pension -def DB_amount(w, e, n, j, p): +def DB_amount(w, e, n, t, j, method, p): r""" Calculate public pension from a defined benefits system. .. math:: pension{j,s,t} = \biggl[\frac{\sum_{s=R-ny}^{R-1}w_{t}e_{j,s,t} - n_{j,s,t}}{ny}\biggr]\times Cy \times \alpha_{DB} + n_{j,s,t}}{ny}\biggr]\times Cy \times \alpha^{DB}_{j,t} \quad \forall s > R + The replacement rate :math:`\alpha^{DB}_{j,t}` may vary over time + and across lifetime income groups. The rate applied to a benefit is + the one in effect in the period the benefit is paid. + Args: w (array_like): real wage rate e (Numpy array): effective labor units @@ -265,7 +269,12 @@ def DB_amount(w, e, n, j, p): (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 + t (int): model period of the first element of ``n`` (for a + partial or full lifetime path) or of the first row of ``n`` + (for the full time path) + j (int): index of lifetime income group, None if all groups + method (str): adjusts calculation dimensions based on 'SS', + 'TPI', or 'TPI_scalar' p (OG-Core Specifications object): model parameters Returns: @@ -297,6 +306,12 @@ def DB_amount(w, e, n, j, p): # 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) + # This household is age S - per_rmn at time t, so its benefit + # at age u is paid in period t + u - (S - per_rmn) + idx = np.clip( + t - (p.S - per_rmn) + np.arange(p.S), 0, p.alpha_db.shape[0] - 1 + ) + alpha_db_u = p.alpha_db[idx, j] DB = np.zeros(p.S) DB = DB_1dim_loop( @@ -313,13 +328,21 @@ def DB_amount(w, e, n, j, p): L_inc_avg, DB, equiv_periods, - p.alpha_db, + alpha_db_u, equiv_yr_contrib, ) DB = DB[-per_rmn:] else: if np.ndim(n) == 1: + if method == "SS": + alpha_db_u = np.full(p.S, p.alpha_db[-1, j]) + elif method == "TPI_scalar": + alpha_db_u = np.full(p.S, p.alpha_db[0, j]) + else: + # full lifetime of a cohort that is age 0 at time t + idx = np.clip(t + np.arange(p.S), 0, p.alpha_db.shape[0] - 1) + alpha_db_u = p.alpha_db[idx, j] DB = np.zeros(p.S) DB = DB_1dim_loop( w, @@ -332,11 +355,19 @@ def DB_amount(w, e, n, j, p): L_inc_avg, DB, equiv_periods, - p.alpha_db, + alpha_db_u, equiv_yr_contrib, ) elif np.ndim(n) == 2: + # cross-section of all S x J households in a single period + if method == "SS": + alpha_db_j = p.alpha_db[-1, :] + elif method == "TPI_scalar": + alpha_db_j = p.alpha_db[0, :] + else: + alpha_db_j = p.alpha_db[t, :] + alpha_db_uj = np.tile(alpha_db_j.reshape(1, p.J), (p.S, 1)) DB = np.zeros((p.S, p.J)) L_inc_avg_sj = np.zeros((equiv_periods, p.J)) DB = DB_2dim_loop( @@ -350,7 +381,7 @@ def DB_amount(w, e, n, j, p): L_inc_avg, DB, equiv_periods, - p.alpha_db, + alpha_db_uj, equiv_yr_contrib, ) @@ -360,6 +391,10 @@ def DB_amount(w, e, n, j, p): 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 + # rate in effect in each period of the path, T x J + idx = np.clip(t + np.arange(T), 0, p.alpha_db.shape[0] - 1) + alpha_db_tj = np.ascontiguousarray(p.alpha_db[idx, :]) + # alpha_db_tj = p.alpha_db[t + np.arange(T), :] DB = DB_3dim_loop( w_path, e_ss, @@ -370,7 +405,7 @@ def DB_amount(w, e, n, j, p): p.J, float(g_y_arr[-1]), equiv_periods, - p.alpha_db, + alpha_db_tj, equiv_yr_contrib, ) @@ -565,7 +600,7 @@ def PS_amount(w, e, n, j, factor, p): return PS -def deriv_theta(r, w, e, Y, per_rmn, factor, p): +def deriv_theta(r, w, e, Y, per_rmn, factor, t, j, p): """ Change in pension benefits for another unit of labor supply for pension system selected @@ -577,6 +612,9 @@ def deriv_theta(r, w, e, Y, per_rmn, factor, p): Y (array_like): GDP per_rmn (int): number of periods remaining in the model factor (scalar): scaling factor converting model units to + t (int): model period + j (int): index of lifetime income group + p (OG-Core Specifications object): model parameters Returns: d_theta (Numpy array): change in pension benefits for another @@ -584,7 +622,7 @@ def deriv_theta(r, w, e, Y, per_rmn, factor, p): """ # TODO: Add SS here... if p.pension_system == "Defined Benefits": - d_theta = deriv_DB(w, e, per_rmn, p) + d_theta = deriv_DB(w, e, per_rmn, t, j, p) d_theta = d_theta[-per_rmn:] elif p.pension_system == "Notional Defined Contribution": d_theta = deriv_NDC(r, w, e, Y, per_rmn, p) @@ -649,7 +687,7 @@ def deriv_NDC(r, w, e, Y, per_rmn, p): return d_theta -def deriv_DB(w, e, per_rmn, p): +def deriv_DB(w, e, per_rmn, t, j, p): r""" Change in DB pension benefits for another unit of labor supply @@ -657,15 +695,23 @@ def deriv_DB(w, e, per_rmn, p): \frac{\partial \theta_{j,u,t+u-s}}{\partial n_{j,s,t}} = \begin{cases} 0 , & \text{if}\ s < R - Cy \\ - w_{t}e_{j,s}\alpha_{DB}\times \frac{Cy}{ny}, + w_{t}e_{j,s}\alpha^{DB}_{j,t+R-s}\times \frac{Cy}{ny}, & \text{if}\ R - Cy <= s < R \\ 0, & \text{if}\ s \geq R \\ \end{cases} + With a replacement rate that varies over time, the change in the + benefit paid at each age u depends on the rate in effect at that + time. This function returns a single derivative for each working + age s, evaluated at the rate in effect when the household retires + (period t + R - s), as an approximation. + Args: w (array_like): real wage rate e (Numpy array): effective labor units per_rmn (int): number of periods remaining in the model + t (int): model period + j (int): index of lifetime income group p (OG-Core Specifications object): model parameters Returns: @@ -679,6 +725,12 @@ def deriv_DB(w, e, per_rmn, p): if per_rmn < (p.S - S_ret + 1): d_theta = np.zeros(p.S) else: + # element s of the loop is a household of age S - per_rmn + s + # at time t, which retires in period t + S_ret - age + age = p.S - per_rmn + np.arange(per_rmn) + idx = np.clip(t + S_ret - age, 0, p.alpha_db.shape[0] - 1) + alpha_db_s = p.alpha_db[idx, j] + # alpha_db_s = p.alpha_db[t + S_ret - age, j] d_theta = deriv_DB_loop( w, e, @@ -686,7 +738,7 @@ def deriv_DB(w, e, per_rmn, p): S_ret, per_rmn, equiv_periods, - p.alpha_db, + alpha_db_s, equiv_yr_contrib, ) return d_theta @@ -868,7 +920,8 @@ def deriv_DB_loop( S_ret (int): retirement age per_rmn (int): number of periods remaining in the model avg_earn_num_years (int): number of years AIME is computed from - alpha_db (scalar): replacement rate + alpha_db (Numpy array): replacement rate for each of the + per_rmn remaining periods yr_contr (scalar): years of contribution Returns: @@ -880,7 +933,9 @@ def deriv_DB_loop( # print("Average earnings years: ", avg_earn_num_years) num_per_retire = S - S_ret for s in range(per_rmn): - d_theta[s] = w[s] * e[s] * alpha_db * (yr_contr / avg_earn_num_years) + d_theta[s] = ( + w[s] * e[s] * alpha_db[s] * (yr_contr / avg_earn_num_years) + ) d_theta[-num_per_retire:] = 0.0 return d_theta @@ -1083,7 +1138,8 @@ def DB_1dim_loop( L_inc_avg (scalar): average labor income DB (Numpy array): pension amount for each household avg_earn_num_years (int): number of years AIME is computed from - alpha_db (scalar): replacement rate + alpha_db (Numpy array): replacement rate in effect at each age + u, length S yr_contr (scalar): years of contribution Returns: @@ -1098,7 +1154,7 @@ def DB_1dim_loop( w[s] / np.exp(g_y[-1] * (u - s)) * e[s] * n[s] ) L_inc_avg = L_inc_avg_s.sum() / avg_earn_num_years - rep_rate = yr_contr * alpha_db + rep_rate = yr_contr * alpha_db[u] DB[u] = rep_rate * L_inc_avg return DB @@ -1133,7 +1189,8 @@ def DB_2dim_loop( L_inc_avg (scalar): average labor income DB (Numpy array): pension amount for each household avg_earn_num_years (int): number of years AIME is computed from - alpha_db (scalar): replacement rate + alpha_db (Numpy array): replacement rate in effect at each age + u for each group j, size S x J yr_contr (scalar): years of contribution Returns: @@ -1146,7 +1203,7 @@ def DB_2dim_loop( w[s] / np.exp(g_y * (u - s)) * e[s, :] * n[s, :] ) L_inc_avg = L_inc_avg_sj.sum(axis=0) / avg_earn_num_years - rep_rate = yr_contr * alpha_db + rep_rate = yr_contr * alpha_db[u, :] DB[u, :] = rep_rate * L_inc_avg return DB @@ -1192,7 +1249,9 @@ def DB_3dim_loop( 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 + alpha_db (Numpy array): replacement rate per year of + contribution in effect in each period t for each group j, + size T x J yr_contr (int): years of contribution Returns: @@ -1214,7 +1273,9 @@ def DB_3dim_loop( e_ss[s_idx, :] * n_hist ) DB[t, u, :] = ( - (L.sum(axis=0) / avg_earn_num_years) * yr_contr * alpha_db + (L.sum(axis=0) / avg_earn_num_years) + * yr_contr + * alpha_db[t, :] ) return DB diff --git a/pyproject.toml b/pyproject.toml index c81d06338..ca8ce2071 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "ogcore" -version = "0.20.4" +version = "0.21.0" authors = [ {name = "Jason DeBacker and Richard W. Evans"}, ] diff --git a/tests/test_pensions.py b/tests/test_pensions.py index 84d6b3624..f7bd9cca1 100644 --- a/tests/test_pensions.py +++ b/tests/test_pensions.py @@ -124,7 +124,7 @@ def test_replacement_rate_vals(n, w, factor, j, p_in, expected): p = Specifications() p.S = 7 -p.alpha_db = 0.2 +p.alpha_db = np.ones((p.T + p.S, p.J)) * 0.2 p.retire = 4 p.avg_earn_num_years = 50 p.yr_contrib = 55 @@ -153,7 +153,7 @@ def test_replacement_rate_vals(n, w, factor, j, p_in, expected): L_inc_avg, DB, equiv_periods, - p.alpha_db, + np.full(p.S, 0.2), # rate in effect at each age u equiv_contrib_periods, ) @@ -205,7 +205,7 @@ def test_DB_1dim_loop(args, DB_loop_expected): per_rmn = p.S p.avg_earn_num_years = 50 p.yr_contrib = 55 -p.alpha_db = 0.2 +p.alpha_db = np.ones((p.T + p.S, p.J)) * 0.2 p.g_y = np.ones(p.T) * 0.03 w = np.array([1.2, 1.1, 1.21, 1, 1.01, 0.99, 0.8]) e = np.array([1.1, 1.11, 0.9, 0.87, 0.87, 0.7, 0.6]) @@ -222,7 +222,7 @@ def test_DB_1dim_loop(args, DB_loop_expected): p.retire, per_rmn, equiv_periods, - p.alpha_db, + np.full(per_rmn, 0.2), # rate in effect at each remaining period equiv_contrib_periods, ) @@ -281,7 +281,7 @@ def test_deriv_PS_loop(args, deriv_PS_loop_expected): p.retire = 4 p.avg_earn_num_years = 50 p.yr_contrib = 55 -p.alpha_db = 0.2 +p.alpha_db = np.ones((p.T + p.S, p.J)) * 0.2 p.g_y = np.ones(p.T) * 0.03 n_ddb1 = np.array([0.4, 0.45, 0.4, 0.42, 0.3, 0.2, 0.2]) w_ddb1 = np.array([1.2, 1.1, 1.21, 1, 1.01, 0.99, 0.8]) @@ -289,7 +289,7 @@ def test_deriv_PS_loop(args, deriv_PS_loop_expected): per_rmn = n_ddb1.shape[0] d_theta_empty = np.zeros_like(per_rmn) deriv_DB_expected1 = np.array([0.352, 0.3256, 0.2904, 0.232, 0.0, 0.0, 0.0]) -args_ddb1 = (w_ddb1, e_ddb1, per_rmn, p) +args_ddb1 = (w_ddb1, e_ddb1, per_rmn, 0, 1, p) #############non-zero d_theta: case 2############ p2 = Specifications() @@ -297,7 +297,7 @@ def test_deriv_PS_loop(args, deriv_PS_loop_expected): p2.retire = 5 p2.last_career_yrs = 2 p2.yr_contrib = 55 -p2.alpha_db = 0.2 +p2.alpha_db = np.ones((p2.T + p2.S, p2.J)) * 0.2 p2.g_y = 0.03 n_ddb2 = np.array([0.45, 0.4, 0.42, 0.3, 0.2, 0.2]) w_ddb1 = np.array([1.1, 1.21, 1, 1.01, 0.99, 0.8]) @@ -305,22 +305,50 @@ def test_deriv_PS_loop(args, deriv_PS_loop_expected): per_rmn = n_ddb2.shape[0] d_theta_empty = np.zeros_like(per_rmn) deriv_DB_expected2 = np.array([0.4884, 0.4356, 0.348, 0.35148, 0.0, 0.0]) -args_ddb2 = (w_ddb1, e_ddb1, per_rmn, p2) +args_ddb2 = (w_ddb1, e_ddb1, per_rmn, 0, 1, p2) + +# Rate that varies over time: the derivative for a household of age s +# at time t is evaluated at the rate in effect when it retires, period +# t + R - s. With S=7, R=4 and the rate doubling from period 2 on, the +# household of age 0 at t=0 retires in period 4 (0.4), age 1 in period 3 +# (0.4), age 2 in period 2 (0.4), age 3 in period 1 (0.2). +p3 = copy.deepcopy(p) +p3.alpha_db = np.ones((p3.T + p3.S, p3.J)) * 0.2 +p3.alpha_db[2:, :] = 0.4 +deriv_DB_expected3 = deriv_DB_expected1 * np.array( + [2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0] +) +args_ddb3 = ( + np.array([1.2, 1.1, 1.21, 1, 1.01, 0.99, 0.8]), + np.array([1.1, 1.11, 0.9, 0.87, 0.87, 0.7, 0.6]), + 7, + 0, + 1, + p3, +) -test_data = [(args_ddb1, deriv_DB_expected1), (args_ddb2, deriv_DB_expected2)] +test_data = [ + (args_ddb1, deriv_DB_expected1), + (args_ddb2, deriv_DB_expected2), + (args_ddb3, deriv_DB_expected3), +] @pytest.mark.parametrize( "args,deriv_DB_expected", test_data, - ids=["non-zero d_theta: case 1", "non-zero d_theta: case 2"], + ids=[ + "non-zero d_theta: case 1", + "non-zero d_theta: case 2", + "time-varying alpha_db", + ], ) def test_deriv_DB(args, deriv_DB_expected): """ Test of the pensions.deriv_DB() function. """ - w, e, per_rmn, p = args - deriv_DB = pensions.deriv_DB(w, e, per_rmn, p) + w, e, per_rmn, t, j, p = args + deriv_DB = pensions.deriv_DB(w, e, per_rmn, t, j, p) assert np.allclose(deriv_DB, deriv_DB_expected) @@ -440,7 +468,7 @@ def test_deriv_NDC(args, d_NDC_expected): p.retire = 4 p.avg_earn_num_years = 50 p.yr_contrib = 55 -p.alpha_db = 0.2 +p.alpha_db = np.ones((p.T + p.S, p.J)) * 0.2 p.g_y = np.ones(p.T) * 0.03 w_db = np.array([1.2, 1.1, 1.21, 1.0, 1.01, 0.99, 0.8]) e_db = np.array([1.1, 1.11, 0.9, 0.87, 0.87, 0.7, 0.6]) @@ -647,14 +675,23 @@ def test_delta_ret_loop(args, dir_delta_ret_expected): p.retire = 4 p.avg_earn_num_years = 50 p.yr_contrib = 55 -p.alpha_db = 0.2 +p.alpha_db = np.ones((p.T + p.S, p.J)) * 0.2 p.g_y = np.ones(p.T) * 0.03 j = 1 w = np.array([1.2, 1.1, 1.21, 1.0, 1.01, 0.99, 0.8]) e = np.array([1.1, 1.11, 0.9, 0.87, 0.87, 0.7, 0.6]) n = np.array([0.4, 0.45, 0.4, 0.42, 0.3, 0.2, 0.2]) DB_expected1 = np.array([0, 0, 0, 0, 0.337864778, 0.327879365, 0.318189065]) -args1 = (w, e, n, j, p) +args1 = (w, e, n, None, j, "SS", p) + +# Same steady state, but with a replacement rate that differs across +# lifetime income groups: group j=1 keeps 0.2 and the other groups get +# a very different rate, so the result for j=1 is unchanged only if the +# group-specific rate is used. +p_J = copy.deepcopy(p) +p_J.alpha_db = np.ones((p_J.T + p_J.S, p_J.J)) * 0.9 +p_J.alpha_db[:, 1] = 0.2 +args1_J = (w, e, n, None, j, "SS", p_J) #####################Incomplete############ p2 = Specifications() @@ -662,7 +699,7 @@ def test_delta_ret_loop(args, dir_delta_ret_expected): p2.retire = 4 p.avg_earn_num_years = 50 p2.yr_contrib = 55 -p2.alpha_db = 0.2 +p2.alpha_db = np.ones((p2.T + p2.S, p2.J)) * 0.2 p2.g_y = np.ones(p2.T) * 0.03 j = 1 w2 = np.array([1.21, 1.0, 1.01, 0.99, 0.8]) @@ -692,24 +729,144 @@ def test_delta_ret_loop(args, dir_delta_ret_expected): ) p2.e = e2 DB_expected2 = np.array([0, 0, 0.30593337, 0.29689167, 0.2881172]) -args2 = (w2, e2, n2, j, p2) +args2 = (w2, e2, n2, 0, j, "TPI", p2) + +# Incomplete lifetime straddling a change in the replacement rate. This +# household is age 2 at t=0 (S=7, 5 periods remaining) so its benefits +# at ages 4, 5, 6 are paid in periods 2, 3, 4. The rate doubles from +# period 3 on, so the age-4 benefit is unchanged and the age-5 and +# age-6 benefits double: the rate in effect when the benefit is paid +# applies, not the rate at retirement. +p2_tv = copy.deepcopy(p2) +p2_tv.alpha_db = np.ones((p2_tv.T + p2_tv.S, p2_tv.J)) * 0.2 +p2_tv.alpha_db[3:, :] = 0.4 +DB_expected2_tv = DB_expected2 * np.array([1.0, 1.0, 1.0, 2.0, 2.0]) +args2_tv = (w2, e2, n2, 0, j, "TPI", p2_tv) -test_data = [(args1, DB_expected1), (args2, DB_expected2)] +test_data = [ + (args1, DB_expected1), + (args1_J, DB_expected1), + (args2, DB_expected2), + (args2_tv, DB_expected2_tv), +] @pytest.mark.parametrize( - "args,DB_expected", test_data, ids=["SS/Complete", "Incomplete"] + "args,DB_expected", + test_data, + ids=[ + "SS/Complete", + "SS/Complete, alpha_db varies by J", + "Incomplete", + "Incomplete, alpha_db varies over time", + ], ) def test_DB(args, DB_expected): """ - Test of the pensions.get_DB() function. + Test of the pensions.DB_amount() function. """ - w, e, n, j, p = args - DB = pensions.DB_amount(w, e, n, j, p) + w, e, n, t, j, method, p = args + DB = pensions.DB_amount(w, e, n, t, j, method, p) assert np.allclose(DB, DB_expected) +def test_DB_cross_section_varies_by_J(): + """ + In the steady state with all groups computed at once (S x J inputs), + each group's DB benefit should scale with its own replacement rate. + """ + p_xs = Specifications() + p_xs.S = 7 + p_xs.J = 3 + p_xs.retire = 4 + p_xs.avg_earn_num_years = 50 + p_xs.yr_contrib = 55 + p_xs.g_y = 0.03 + alpha_J = np.array([0.1, 0.2, 0.4]) + p_xs.alpha_db = np.tile(alpha_J.reshape(1, 3), (p_xs.T + p_xs.S, 1)) + w_xs = np.array([1.2, 1.1, 1.21, 1.0, 1.01, 0.99, 0.8]) + e_xs = np.tile( + np.array([1.1, 1.11, 0.9, 0.87, 0.87, 0.7, 0.6]).reshape(7, 1), + (1, 3), + ) + n_xs = np.tile( + np.array([0.4, 0.45, 0.4, 0.42, 0.3, 0.2, 0.2]).reshape(7, 1), (1, 3) + ) + DB = pensions.DB_amount(w_xs, e_xs, n_xs, None, None, "SS", p_xs) + assert DB.shape == (7, 3) + # identical earnings histories, so benefits differ only by the rate + assert np.allclose(DB[:, 1], DB_expected1) + assert np.allclose(DB[:, 0], DB_expected1 * 0.5) + assert np.allclose(DB[:, 2], DB_expected1 * 2.0) + + +def test_DB_time_path_consistent_with_cohort_solves(): + """ + With a replacement rate that varies over time and by group, the DB + benefits computed for the full time path at once (the T x S x J + branch used for aggregates) must agree with the benefits computed + cohort by cohort (the partial- and full-lifetime branches used in + the household solves), since both index the rate by the period in + which the benefit is paid. + """ + p_tp = Specifications() + p_tp.S = 7 + p_tp.J = 2 + p_tp.T = 12 + p_tp.retire = 4 + p_tp.avg_earn_num_years = 50 + p_tp.yr_contrib = 55 + p_tp.g_y = 0.03 + rng = np.random.RandomState(3) + p_tp.alpha_db = 0.1 + 0.3 * rng.rand(p_tp.T + p_tp.S, p_tp.J) + p_tp.e = np.tile( + np.array([1.1, 1.11, 0.9, 0.87, 0.87, 0.7, 0.6]).reshape(1, 7, 1), + (p_tp.T + p_tp.S, 1, 2), + ) + p_tp.e[:, :, 1] *= 1.5 + # constant wage and labor supply over time so that the pre-time-path + # history used by the two branches is the same + w_path = np.ones(p_tp.T) * 1.2 + n_ss = np.tile( + np.array([0.4, 0.45, 0.4, 0.42, 0.3, 0.2, 0.2]).reshape(7, 1), (1, 2) + ) + n_path = np.tile(n_ss.reshape(1, 7, 2), (p_tp.T, 1, 1)) + p_tp.n_preTP = n_ss + DB_3D = pensions.DB_amount(w_path, None, n_path, 0, None, "TPI", p_tp) + assert DB_3D.shape == (p_tp.T, p_tp.S, p_tp.J) + for j in range(p_tp.J): + # cohorts alive at t=0, of age a, with S - a periods remaining + for a in range(1, p_tp.S): + per_rmn = p_tp.S - a + DB_1D = pensions.DB_amount( + w_path[:per_rmn], + p_tp.e[-1, a:, j], + n_ss[a:, j], + 0, + j, + "TPI", + p_tp, + ) + path_vals = np.array( + [DB_3D[tt, a + tt, j] for tt in range(per_rmn)] + ) + assert np.allclose(DB_1D, path_vals), (j, a) + # cohorts born at t >= 0 with full lifetimes inside the path + for t0 in range(p_tp.T - p_tp.S + 1): + DB_1D = pensions.DB_amount( + w_path[t0 : t0 + p_tp.S], + p_tp.e[-1, :, j], + n_ss[:, j], + t0, + j, + "TPI", + p_tp, + ) + path_vals = np.array([DB_3D[t0 + u, u, j] for u in range(p_tp.S)]) + assert np.allclose(DB_1D, path_vals), (j, t0) + + ################pension benefit derivative: DB############ p = Specifications() p.pension_system = "Defined Benefits" @@ -718,7 +875,7 @@ def test_DB(args, DB_expected): per_rmn = p.S p.avg_earn_num_years = 50 p.yr_contrib = 55 -p.alpha_db = 0.2 +p.alpha_db = np.ones((p.T + p.S, p.J)) * 0.2 w_ddb = np.array([1.2, 1.1, 1.21, 1, 1.01, 0.99, 0.8]) e_ddb = np.array([1.1, 1.11, 0.9, 0.87, 0.87, 0.7, 0.6]) p.g_y = np.ones(p.T) * 0.03 @@ -726,7 +883,7 @@ def test_DB(args, DB_expected): r_ddb = np.ones(p.T) * 0.03 factor = 2 d_theta_expected_ddb = np.array([0.352, 0.3256, 0.2904, 0.232, 0.0, 0.0, 0.0]) -args_ddb = (r_ddb, w_ddb, e_ddb, Y, per_rmn, factor, p) +args_ddb = (r_ddb, w_ddb, e_ddb, Y, per_rmn, factor, 0, 1, p) ################pension benefit derivative: NDC############ p2 = Specifications() @@ -750,7 +907,7 @@ def test_DB(args, DB_expected): # d_theta_expected_dndc = np.array([0.757437326, 0.680222841, 0, 0]) Y = None factor = 2 -args_dndc = (r_dndc, w_dndc, e_dndc, Y, per_rmn, factor, p2) +args_dndc = (r_dndc, w_dndc, e_dndc, Y, per_rmn, factor, 0, 1, p2) ################pension benefit derivative: PS############ @@ -769,7 +926,7 @@ def test_DB(args, DB_expected): p3.omega_SS = omegas per_rmn_ps = 5 d_theta_expected_dps = np.array([0.0026136, 0.002088, 0, 0, 0]) -args_dps = (r_dps, w_dps, e_dps, Y, per_rmn_ps, factor, p3) +args_dps = (r_dps, w_dps, e_dps, Y, per_rmn_ps, factor, 0, 1, p3) test_data = [ (args_ddb, d_theta_expected_ddb), @@ -785,8 +942,8 @@ def test_deriv_theta(args, d_theta_expected): """ Test of pensions.deriv_theta """ - r, w, e, Y, per_rmn, factor, p = args - d_theta = pensions.deriv_theta(r, w, e, Y, per_rmn, factor, p) + r, w, e, Y, per_rmn, factor, t, j, p = args + d_theta = pensions.deriv_theta(r, w, e, Y, per_rmn, factor, t, j, p) assert np.allclose(d_theta, d_theta_expected) @@ -983,7 +1140,7 @@ def test_get_PS(args, PS_expected): [ ( "Defined Benefits", - {"alpha_db": 0.02, "yr_contrib": 35, "avg_earn_num_years": 40}, + {"alpha_db": [[0.02]], "yr_contrib": 35, "avg_earn_num_years": 40}, ), ("Points System", {"vpoint": 0.5}), # The NDC growth-rate settings (ndc_growth_rate, @@ -1026,7 +1183,7 @@ def test_pension_amount_with_real_specifications(system, updates): "system,updates", [ ("US-Style Social Security", {}), - ("Defined Benefits", {"alpha_db": 0.01, "yr_contrib": 40}), + ("Defined Benefits", {"alpha_db": [[0.01]], "yr_contrib": 40}), ( "Points System", {"vpoint": 0.4, "points_growth_rate": "LR GDP"}, @@ -1093,7 +1250,7 @@ def test_SS_solve_defined_benefits(tmp_path): p.update_specifications( { "pension_system": "Defined Benefits", - "alpha_db": 0.02, + "alpha_db": [[0.02]], "yr_contrib": 35, "avg_earn_num_years": 40, } diff --git a/uv.lock b/uv.lock index 5990da103..d556da9ff 100644 --- a/uv.lock +++ b/uv.lock @@ -1521,7 +1521,7 @@ wheels = [ [[package]] name = "ogcore" -version = "0.20.4" +version = "0.21.0" source = { editable = "." } dependencies = [ { name = "dask" },