Allow defined benefits replacement rate to vary by T and J

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

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  

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}
   ```

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"

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": {

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

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,23 +249,32 @@ 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
         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
+        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,14 +612,17 @@ 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
             unit of labor supply
     """
     # 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,23 +687,31 @@ 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
 
     .. math::
         \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,14 +725,20 @@ 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,
             p.S,
             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