Add --use-tob flag for TOB revenue calibration

Adds calibration targets for Taxation of Benefits (TOB) revenue:
- OASDI TOB: tob_revenue_oasdi (tier 1 taxation)
- HI TOB: tob_revenue_medicare_hi (tier 2 taxation)

Changes:
- Add load_oasdi_tob_projections() and load_hi_tob_projections() to ssa_data.py
- Add TOB constraint parameters to calibration.py
- Add --use-tob CLI flag to run_household_projection.py
- Clean hi_tob_billions_nominal_usd column in social_security_aux.csv
- Update README.md documentation

Closes #459

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <noreply@anthropic.com>

Author: baogorek

changelog_entry.yaml

@@ -0,0 +1,4 @@
+- bump: minor
+  changes:
+    added:
+    - Added --use-tob flag for TOB (Taxation of Benefits) revenue calibration targeting OASDI and HI trust fund revenue

policyengine_us_data/datasets/cps/long_term/README.md

@@ -6,8 +6,8 @@
 Run projections using `run_household_projection.py`:
 
 ```bash
-# Recommended: GREG with all three constraint types
-python run_household_projection.py 2100 --greg --use-ss --use-payroll --save-h5
+# Recommended: GREG with all constraint types
+python run_household_projection.py 2100 --greg --use-ss --use-payroll --use-tob --save-h5
 
 # IPF with only age distribution constraints (faster, less accurate)
 python run_household_projection.py 2050
@@ -21,6 +21,7 @@ python run_household_projection.py 2100 --greg --use-ss
 - `--greg`: Use GREG calibration instead of IPF
 - `--use-ss`: Include Social Security benefit totals as calibration target (requires `--greg`)
 - `--use-payroll`: Include taxable payroll totals as calibration target (requires `--greg`)
+- `--use-tob`: Include TOB (Taxation of Benefits) revenue as calibration target (requires `--greg`)
 - `--save-h5`: Save year-specific .h5 files to `./projected_datasets/` directory
 
 **Estimated runtime:** ~2 minutes/year without `--save-h5`, ~3 minutes/year with `--save-h5`
@@ -36,7 +37,7 @@ python run_household_projection.py 2100 --greg --use-ss
 
 **GREG (Generalized Regression Estimator)**
 - Solves for weights matching multiple constraints simultaneously
-- Can enforce age distribution + Social Security benefits + taxable payroll
+- Can enforce age distribution + Social Security benefits + taxable payroll + TOB revenue
 - One-shot solution using `samplics` package
 - **Recommended** for accurate long-term projections
 
@@ -57,14 +58,20 @@ python run_household_projection.py 2100 --greg --use-ss
    - Calculated as: `taxable_earnings_for_social_security` + `social_security_taxable_self_employment_income`
    - Source: SSA Trustee Report 2024 (`social_security_aux.csv`)
 
+4. **TOB Revenue** (`--use-tob`, GREG only)
+   - Taxation of Benefits revenue for OASDI and Medicare HI trust funds
+   - OASDI: `tob_revenue_oasdi` (tier 1 taxation, 0-50% of benefits)
+   - HI: `tob_revenue_medicare_hi` (tier 2 taxation, 50-85% of benefits)
+   - Source: SSA Trustee Report 2024 (`social_security_aux.csv`)
+
 ---
 
 ### Data Sources
 
 All data from **SSA 2024 Trustee Report**:
 
 - `SSPopJul_TR2024.csv` - Population projections 2025-2100 by single year of age
-- `social_security_aux.csv` - OASDI costs and taxable payroll projections 2025-2100
+- `social_security_aux.csv` - OASDI costs, taxable payroll, and TOB revenue projections 2025-2100
   - Extracted from `SingleYearTRTables_TR2025.xlsx` Table VI.G9 using `extract_ssa_costs.py`
 
 Files located in: `policyengine_us_data/storage/`

policyengine_us_data/datasets/cps/long_term/calibration.py

@@ -85,6 +85,10 @@ def calibrate_greg(
     payroll_target=None,
     h6_income_values=None,
     h6_revenue_target=None,
+    oasdi_tob_values=None,
+    oasdi_tob_target=None,
+    hi_tob_values=None,
+    hi_tob_target=None,
     n_ages=86,
 ):
     """
@@ -101,6 +105,10 @@ def calibrate_greg(
         payroll_target: Optional taxable payroll target total
         h6_income_values: Optional H6 reform income values per household
         h6_revenue_target: Optional H6 reform total revenue impact target
+        oasdi_tob_values: Optional OASDI TOB revenue values per household
+        oasdi_tob_target: Optional OASDI TOB revenue target total
+        hi_tob_values: Optional HI TOB revenue values per household
+        hi_tob_target: Optional HI TOB revenue target total
         n_ages: Number of age groups
 
     Returns:
@@ -116,6 +124,8 @@ def calibrate_greg(
         (ss_values is not None and ss_target is not None)
         or (payroll_values is not None and payroll_target is not None)
         or (h6_income_values is not None and h6_revenue_target is not None)
+        or (oasdi_tob_values is not None and oasdi_tob_target is not None)
+        or (hi_tob_values is not None and hi_tob_target is not None)
     )
 
     if needs_aux_df:
@@ -135,6 +145,14 @@ def calibrate_greg(
             aux_df["h6_revenue"] = h6_income_values
             controls["h6_revenue"] = h6_revenue_target
 
+        if oasdi_tob_values is not None and oasdi_tob_target is not None:
+            aux_df["oasdi_tob"] = oasdi_tob_values
+            controls["oasdi_tob"] = oasdi_tob_target
+
+        if hi_tob_values is not None and hi_tob_target is not None:
+            aux_df["hi_tob"] = hi_tob_values
+            controls["hi_tob"] = hi_tob_target
+
         aux_vars = aux_df
     else:
         aux_vars = X
@@ -160,6 +178,10 @@ def calibrate_weights(
     payroll_target=None,
     h6_income_values=None,
     h6_revenue_target=None,
+    oasdi_tob_values=None,
+    oasdi_tob_target=None,
+    hi_tob_values=None,
+    hi_tob_target=None,
     n_ages=86,
     max_iters=100,
     tol=1e-6,
@@ -180,6 +202,10 @@ def calibrate_weights(
         payroll_target: Optional payroll target (for GREG with payroll)
         h6_income_values: Optional H6 reform income values per household
         h6_revenue_target: Optional H6 reform total revenue impact target
+        oasdi_tob_values: Optional OASDI TOB revenue values per household
+        oasdi_tob_target: Optional OASDI TOB revenue target total
+        hi_tob_values: Optional HI TOB revenue values per household
+        hi_tob_target: Optional HI TOB revenue target total
         n_ages: Number of age groups
         max_iters: Max iterations for IPF
         tol: Convergence tolerance for IPF
@@ -204,6 +230,10 @@ def calibrate_weights(
                 payroll_target,
                 h6_income_values,
                 h6_revenue_target,
+                oasdi_tob_values,
+                oasdi_tob_target,
+                hi_tob_values,
+                hi_tob_target,
                 n_ages,
             )
         except Exception as e:

policyengine_us_data/datasets/cps/long_term/run_household_projection.py

@@ -3,19 +3,20 @@
 
 
 Usage:
-    python run_household_projection.py [START_YEAR] [END_YEAR] [--greg] [--use-ss] [--use-payroll] [--use-h6-reform] [--save-h5]
+    python run_household_projection.py [START_YEAR] [END_YEAR] [--greg] [--use-ss] [--use-payroll] [--use-h6-reform] [--use-tob] [--save-h5]
 
     START_YEAR: Optional starting year (default: 2025)
     END_YEAR: Optional ending year (default: 2035)
     --greg: Use GREG calibration instead of IPF (optional)
     --use-ss: Include Social Security benefit totals as calibration target (requires --greg)
     --use-payroll: Include taxable payroll totals as calibration target (requires --greg)
     --use-h6-reform: Include H6 reform income impact ratio as calibration target (requires --greg)
+    --use-tob: Include TOB (Taxation of Benefits) revenue as calibration target (requires --greg)
     --save-h5: Save year-specific .h5 files with calibrated weights to ./projected_datasets/
 
 Examples:
     python run_household_projection.py 2045 2045 --greg --use-ss  # single year
-    python run_household_projection.py 2025 2100 --greg --use-ss --use-payroll --use-h6-reform --save-h5
+    python run_household_projection.py 2025 2100 --greg --use-ss --use-payroll --use-tob --save-h5
 """
 
 import sys
@@ -256,6 +257,16 @@ def create_h6_reform():
         USE_GREG = True
     from ssa_data import load_h6_income_rate_change
 
+USE_TOB = "--use-tob" in sys.argv
+if USE_TOB:
+    sys.argv.remove("--use-tob")
+    if not USE_GREG:
+        print(
+            "Warning: --use-tob requires --greg, enabling GREG automatically"
+        )
+        USE_GREG = True
+    from ssa_data import load_oasdi_tob_projections, load_hi_tob_projections
+
 SAVE_H5 = "--save-h5" in sys.argv
 if SAVE_H5:
     sys.argv.remove("--save-h5")
@@ -286,6 +297,8 @@ def create_h6_reform():
     print(f"  Including taxable payroll constraint: Yes")
 if USE_H6_REFORM:
     print(f"  Including H6 reform income impact constraint: Yes")
+if USE_TOB:
+    print(f"  Including TOB revenue constraint: Yes")
 if SAVE_H5:
     print(f"  Saving year-specific .h5 files: Yes (to {OUTPUT_DIR}/)")
     os.makedirs(OUTPUT_DIR, exist_ok=True)
@@ -464,6 +477,33 @@ def create_h6_reform():
             del reform_sim
             gc.collect()
 
+    oasdi_tob_values = None
+    oasdi_tob_target = None
+    hi_tob_values = None
+    hi_tob_target = None
+    if USE_TOB:
+        oasdi_tob_hh = sim.calculate(
+            "tob_revenue_oasdi", period=year, map_to="household"
+        )
+        oasdi_tob_values = oasdi_tob_hh.values
+        oasdi_tob_target = load_oasdi_tob_projections(year)
+
+        hi_tob_hh = sim.calculate(
+            "tob_revenue_medicare_hi", period=year, map_to="household"
+        )
+        hi_tob_values = hi_tob_hh.values
+        hi_tob_target = load_hi_tob_projections(year)
+
+        if year in display_years:
+            oasdi_baseline = np.sum(oasdi_tob_values * baseline_weights)
+            hi_baseline = np.sum(hi_tob_values * baseline_weights)
+            print(
+                f"  [DEBUG {year}] OASDI TOB baseline: ${oasdi_baseline/1e9:.1f}B, target: ${oasdi_tob_target/1e9:.1f}B"
+            )
+            print(
+                f"  [DEBUG {year}] HI TOB baseline: ${hi_baseline/1e9:.1f}B, target: ${hi_tob_target/1e9:.1f}B"
+            )
+
     y_target = target_matrix[:, year_idx]
 
     w_new, iterations = calibrate_weights(
@@ -478,13 +518,19 @@ def create_h6_reform():
         payroll_target=payroll_target,
         h6_income_values=h6_income_values,
         h6_revenue_target=h6_revenue_target,
+        oasdi_tob_values=oasdi_tob_values,
+        oasdi_tob_target=oasdi_tob_target,
+        hi_tob_values=hi_tob_values,
+        hi_tob_target=hi_tob_target,
         n_ages=n_ages,
         max_iters=100,
         tol=1e-6,
         verbose=False,
     )
 
-    if year in display_years and (USE_SS or USE_PAYROLL or USE_H6_REFORM):
+    if year in display_years and (
+        USE_SS or USE_PAYROLL or USE_H6_REFORM or USE_TOB
+    ):
         if USE_SS:
             ss_achieved = np.sum(ss_values * w_new)
             print(
@@ -507,6 +553,15 @@ def create_h6_reform():
             print(
                 f"  [DEBUG {year}] H6 achieved revenue: ${h6_revenue_achieved/1e9:.3f}B (error: {error_pct:.1f}%)"
             )
+        if USE_TOB:
+            oasdi_achieved = np.sum(oasdi_tob_values * w_new)
+            hi_achieved = np.sum(hi_tob_values * w_new)
+            print(
+                f"  [DEBUG {year}] OASDI TOB achieved: ${oasdi_achieved/1e9:.1f}B (error: {(oasdi_achieved - oasdi_tob_target)/oasdi_tob_target*100:.1f}%)"
+            )
+            print(
+                f"  [DEBUG {year}] HI TOB achieved: ${hi_achieved/1e9:.1f}B (error: {(hi_achieved - hi_tob_target)/hi_tob_target*100:.1f}%)"
+            )
 
     weights_matrix[:, year_idx] = w_new
     baseline_weights_matrix[:, year_idx] = baseline_weights

policyengine_us_data/datasets/cps/long_term/ssa_data.py

@@ -89,3 +89,39 @@ def load_h6_income_rate_change(year):
     row = df[df["year"] == year]
     # CSV stores as percentage (e.g., -0.18), convert to decimal
     return row["h6_income_rate_change"].values[0] / 100
+
+
+def load_oasdi_tob_projections(year):
+    """
+    Load OASDI TOB (Taxation of Benefits) revenue target for a given year.
+
+    Args:
+        year: Year to load OASDI TOB revenue for
+
+    Returns:
+        Total OASDI TOB revenue in nominal dollars
+    """
+    csv_path = STORAGE_FOLDER / "social_security_aux.csv"
+    df = pd.read_csv(csv_path)
+
+    row = df[df["year"] == year]
+    nominal_billions = row["oadsi_tob_billions_nominal_usd"].values[0]
+    return nominal_billions * 1e9
+
+
+def load_hi_tob_projections(year):
+    """
+    Load HI (Medicare) TOB revenue target for a given year.
+
+    Args:
+        year: Year to load HI TOB revenue for
+
+    Returns:
+        Total HI TOB revenue in nominal dollars
+    """
+    csv_path = STORAGE_FOLDER / "social_security_aux.csv"
+    df = pd.read_csv(csv_path)
+
+    row = df[df["year"] == year]
+    nominal_billions = row["hi_tob_billions_nominal_usd"].values[0]
+    return nominal_billions * 1e9