Stop counting scrap twice in the emission factors (#1293)

bedrock/transform/eeio/__tests__/test_useeio_b_steps.py

@@ -60,6 +60,12 @@ def test_useeio_b_adjust_divide_transform_steps(
         lambda original_year, target_year: ratio,
     )
     monkeypatch.setattr(dc, "derive_Vnorm_scrap_corrected", lambda: vnorm)
+    # Gross output includes scrap, so this branch divides by 1 - s (#1293).
+    monkeypatch.setattr(
+        dc,
+        "scrap_share_of_output",
+        lambda: pd.Series([0.1, 0.0], index=["I1", "I2"]),
+    )
 
     try:
         out = dc.derive_B_via_vnorm()
@@ -69,10 +75,10 @@ def test_useeio_b_adjust_divide_transform_steps(
 
     # Step 1: adjusted x = x_nominal / ratio = [125, 160]
     # Step 2: Bi = E / adjusted x = [0.08, 0.125]
-    # Step 3: B = Bi @ I = Bi
+    # Step 3: B = Bi @ (I / (1 - s)), with s = [0.1, 0]
     expected = pd.DataFrame(
         [
-            [0.08, 0.125],  # CO2
+            [0.08 / 0.9, 0.125],  # CO2
             [0.0, 0.0],  # CH4
             [0.0, 0.0],  # N2O
             [0.0, 0.0],  # HFCs

bedrock/transform/eeio/usa_matrices.py

@@ -308,12 +308,53 @@ def derive_q() -> pd.Series[float]:
     return q
 
 
+def scrap_share_of_output() -> pd.Series[float]:
+    """Each industry's scrap (``S00401``) over its output **including** scrap.
+
+    That is USEEIO's ``s``. The Cornerstone Make has no scrap, so
+    ``derive_Vnorm_scrap_corrected`` needs no adjustment where ``x`` is its row
+    sums. Where ``x`` is BEA gross output, which includes scrap, B divides the
+    market shares by ``1 - s`` (``derive_B_via_vnorm``, #1293).
+
+    A BEA industry split into several Cornerstone industries (the waste
+    children of ``562000``, the electricity children) shares its scrap across
+    them **by output**, not by the correspondence's equal weights. So every
+    child gets its parent's scrap share, as in USEEIO's waste disaggregation.
+    """
+    V = derive_V()
+    x = compute_x(V=V)
+    corresp = industry_corresp()
+    by_output = corresp.gt(0).mul(x.reindex(corresp.index).fillna(0.0), axis=0)
+    totals = by_output.sum(axis=0)
+    # A BEA industry whose children have no output keeps the correspondence.
+    weights = by_output.div(totals.where(totals > 0), axis=1).fillna(corresp)
+    scrap = weights @ load_detail_V_usa().loc[:, 'S00401']
+    if get_usa_config().implement_electricity_disaggregation:
+        parent = float(scrap.get(ELECTRICITY_AGGREGATE_SECTOR, 0.0))
+        scrap = scrap.drop(ELECTRICITY_AGGREGATE_SECTOR, errors='ignore')
+        children = [c for c in ELECTRICITY_DISAGG_SECTORS if c in x.index]
+        child_x = float(x.reindex(children).sum())
+        for code in children:
+            scrap.loc[code] = parent * float(x[code]) / child_x if child_x else 0.0
+    scrap = scrap.reindex(V.index, fill_value=0.0)
+    return (scrap / (x.reindex(V.index, fill_value=0.0) + scrap)).fillna(0.0)
+
+
 @functools.cache
 def derive_Vnorm_scrap_corrected(
     apply_inflation: bool | None = None,
     target_year: int = 0,
 ) -> pd.DataFrame:
-    """Scrap-corrected V norm. Inflation is applied via ``derive_V``.
+    """Market shares ``V q̂⁻¹``, net of scrap. Inflation is applied via ``derive_V``.
+
+    USEEIO divides market shares by ``1 - s``, with ``s`` scrap's share of an
+    industry's output **including** scrap, so that the industry's inputs and
+    emissions land on its other products. The Cornerstone Make table has no
+    scrap: ``S00401`` maps with weight 0, so ``x`` is already output without
+    scrap, and ``U / x`` and ``E / x`` already land on the other products.
+    Dividing by ``1 - s`` again would count scrap twice (#1293). The plain
+    market shares give exactly USEEIO's A and B:
+    ``U x_incl⁻¹ V_n / (1 - scrap / x_incl) = U x_excl⁻¹ V_n``.
 
     When ``target_year`` is not positive and inflation is on, uses
     ``USAConfig.model_base_year``.
@@ -327,27 +368,8 @@ def derive_Vnorm_scrap_corrected(
 
     V = derive_V(use_inflation, effective_year)
     q = compute_q(V=V)
-    x = compute_x(V=V)
-
-    Vnorm = compute_Vnorm_matrix(V=V, q=q)
 
-    scrap_detail = load_detail_V_usa().loc[:, 'S00401']
-    scrap_fraction = industry_corresp() @ scrap_detail
-    if get_usa_config().implement_electricity_disaggregation:
-        parent_scrap = float(scrap_fraction.get(ELECTRICITY_AGGREGATE_SECTOR, 0.0))
-        scrap_fraction = scrap_fraction.drop(
-            ELECTRICITY_AGGREGATE_SECTOR, errors='ignore'
-        )
-        for code in ELECTRICITY_DISAGG_SECTORS:
-            scrap_fraction.loc[code] = parent_scrap
-    # The delivery sector has no scrap: the reindex below gives it zero.
-    scrap_fraction = scrap_fraction.reindex(V.index, fill_value=0.0)
-    x_aligned = x.reindex(V.index, fill_value=0.0)
-
-    V_scrap_corrected = Vnorm.divide(
-        (1.0 - (scrap_fraction / x_aligned).fillna(0.0)), axis=0
-    )
-    V_scrap_corrected = V_scrap_corrected.reindex(
+    V_scrap_corrected = compute_Vnorm_matrix(V=V, q=q).reindex(
         index=V.index, columns=V.columns, fill_value=0.0
     )
     validate_cornerstone(V_scrap_corrected, 'V')
@@ -834,6 +856,12 @@ def derive_B_via_vnorm() -> pd.DataFrame:
     else:
         x = derive_x_after_redefinition() if cfg.use_ghg_year_x_in_B else derive_x()
     Vnorm = derive_Vnorm_scrap_corrected()
+    if cfg.usa_detail_io_source != 'nowcast' and (
+        cfg.deflate_x_to_detail_io_year_for_B or cfg.use_ghg_year_x_in_B
+    ):
+        # BEA gross output includes scrap, so the scrap correction applies here,
+        # as in USEEIO: divide by 1 - s with s over output including scrap (#1293).
+        Vnorm = Vnorm.divide(1.0 - scrap_share_of_output(), axis=0)
     Bi = E.divide(x, axis=1).fillna(0.0)
     B = Bi @ Vnorm
     validate_cornerstone(B, 'B')

bedrock/utils/economic/inflation_helpers_usa.py

@@ -423,12 +423,10 @@ def get_vnorm_adjusted_commodity_price_ratio(
     aligned = industry_ratio.reindex(Vnorm.index, fill_value=1.0)
 
     # Normalize V_norm columns to sum to 1 so the dot-product is a true weighted
-    # *average*, not a weighted *sum*. Scrap correction in
-    # `derive_Vnorm_scrap_corrected` is a row-axis scaling that
-    # leaves column sums slightly > 1 — using the un-renormalized form would
-    # produce a non-1 ratio at year=year (off by the column-sum drift, ~5–7%).
-    # Done locally here rather than upstream because A/B matrix consumers want
-    # the un-renormalized scrap-corrected V_norm.
+    # *average*, not a weighted *sum*. `derive_Vnorm_scrap_corrected` returns the
+    # plain market shares (#1293), whose columns already sum to 1, so this is a
+    # guard: it keeps the ratio at 1 for year = year if the shares ever change.
+    # B's gross-output path scales them by 1 / (1 - s) separately.
     column_sums = Vnorm.sum(axis=0)
     weights = Vnorm.divide(column_sums.where(column_sums > 1e-9, 1.0), axis=1)
     commodity_ratio = aligned @ weights

bedrock/utils/validation/__tests__/test_eeio_diagnostics.py

@@ -417,8 +417,8 @@ def test_compare_output_and_L_y(
 @pytest.mark.xfail(
     reason=(
         'Make and Use disagree on commodity output. On MUT v1.0.0_8365fa4 '
-        '(2024), with the F09000 sales in y: 334118 (+4.4%, imports above '
-        'use) and S00402 ($4M on q = 0).'
+        '(2024) one sector fails: 334118 (+4.4%), whose imports exceed its '
+        'use. 2017-2023 pass.'
     ),
 )
 def test_domestic_L_y_from_use_final_demand_equals_q(
@@ -569,6 +569,21 @@ def test_domestic_lci_equals_e(canonical_usa_config: None) -> None:
     )
 
 
+@pytest.mark.eeio_integration
+def test_b_times_q_conserves_emissions(canonical_usa_config: None) -> None:
+    """``B @ q`` sums to E: commodity factors neither add nor lose emissions.
+
+    It failed by 4.0 Mt when the scrap correction divided market shares by
+    ``1 - s`` on an ``x`` that already excluded scrap (#1293). The LCI≈E test
+    above can't see that, because it puts the same market shares on both sides.
+    """
+    B = derive_B()
+    q = derive_Aq_scaled().scaled_q
+    assert float((B.reindex(columns=q.index, fill_value=0.0) @ q).sum()) == (
+        pytest.approx(float(derive_E_usa().values.sum()), rel=1e-9)
+    )
+
+
 @pytest.mark.eeio_integration
 @pytest.mark.xfail(
     reason=(
@@ -600,3 +615,20 @@ def test_total_lci_equals_e(canonical_usa_config: None) -> None:
         f'total LCI≈E failed ({len(result.failing_sectors)} cells): '
         f'{result.failing_sectors[:20]}'
     )
+
+
+def test_output_vs_leontief_skips_the_accounting_sectors() -> None:
+    """S00402 has no output, so q = L @ y isn't checked for it (#1270)."""
+    sectors = ['1111A0', 'S00402']
+    L = pd.DataFrame(np.eye(2), index=sectors, columns=sectors)
+    q = pd.Series([100.0, 0.0], index=sectors)
+
+    off_on_used_goods = compare_output_vs_leontief_x_demand(
+        output=q, L=L, y=pd.Series([100.0, -6.0e6], index=sectors)
+    )
+    assert off_on_used_goods.passed
+
+    off_on_oilseeds = compare_output_vs_leontief_x_demand(
+        output=q, L=L, y=pd.Series([103.0, 0.0], index=sectors)
+    )
+    assert not off_on_oilseeds.passed

bedrock/utils/validation/eeio_diagnostics.py

@@ -28,6 +28,12 @@
 # (e.g. S00402 used goods) compare |diff| against atol instead of (diff / q).
 _VALIDATE_RESULT_ATOL = 1e-4
 
+# BEA's accounting rows: scrap, used and secondhand goods, noncomparable
+# imports and the rest-of-world adjustment. They have no output of their own
+# (S00402's q is 0), so q = L @ y can't hold for them and isn't checked, as in
+# USEEIO 2. Of the four, only S00402 is a Cornerstone commodity today (#1270).
+ACCOUNTING_SECTORS: frozenset[str] = frozenset({'S00300', 'S00401', 'S00402', 'S00900'})
+
 
 @dc.dataclass
 class DiagnosticResult:
@@ -347,7 +353,7 @@ def compare_output_vs_leontief_x_demand(
     Compares the total sector output (commodity or industry) against
     the model result calculation of L @ y.
     Pass/fail: ``|diff| <= tolerance * |value| + atol`` for all sectors (see
-    ``validate_result``).
+    ``validate_result``) except the ``ACCOUNTING_SECTORS``, which aren't checked.
 
     Parameters
     ----------
@@ -387,8 +393,14 @@ def compare_output_vs_leontief_x_demand(
     output_check = backcompute_q_from_L_and_y(L=L, y=y)
     name = "compare output and L * y"
 
+    # L @ y uses every sector; only the comparison leaves the accounting rows out.
+    checked = output.index.difference(list(ACCOUNTING_SECTORS), sort=False)
     d_result = validate_result(
-        name, output, output_check, tolerance=tolerance, include_details=include_details
+        name,
+        output.loc[checked],
+        output_check.loc[checked],
+        tolerance=tolerance,
+        include_details=include_details,
     )
 
     return d_result