Single-sector default, heterogeneous CES epsilon, and UK-calibrated structural parameters

examples/run_oguk.py

@@ -30,17 +30,27 @@
 )
 
 
-def run_steady_state(age_specific: str = "pooled"):
+def run_steady_state(age_specific: str = "pooled", multi_sector: bool = False):
     """Run baseline and reform steady state, print results."""
-    print(f"Solving baseline steady state (age_specific='{age_specific}')...")
+    sector_label = "8-sector" if multi_sector else "1-sector"
+    print(
+        f"Solving baseline steady state (age_specific='{age_specific}', {sector_label})..."
+    )
     t0 = time.time()
-    baseline = solve_steady_state(start_year=2026, age_specific=age_specific)
+    baseline = solve_steady_state(
+        start_year=2026, age_specific=age_specific, multi_sector=multi_sector
+    )
     print(f"  Done in {time.time() - t0:.1f}s")
 
-    print(f"Solving reform steady state (age_specific='{age_specific}')...")
+    print(
+        f"Solving reform steady state (age_specific='{age_specific}', {sector_label})..."
+    )
     t0 = time.time()
     reform = solve_steady_state(
-        start_year=2026, policy=REFORM, age_specific=age_specific
+        start_year=2026,
+        policy=REFORM,
+        age_specific=age_specific,
+        multi_sector=multi_sector,
     )
     print(f"  Done in {time.time() - t0:.1f}s")
 
@@ -87,17 +97,21 @@ def run_steady_state(age_specific: str = "pooled"):
     print("=" * 60)
 
 
-def run_tpi():
+def run_tpi(multi_sector: bool = False):
     """Run baseline and reform transition paths, print results."""
     from dask.distributed import Client
 
-    print("Running baseline + reform transition paths...")
+    sector_label = "8-sector" if multi_sector else "1-sector"
+    print(f"Running baseline + reform transition paths ({sector_label})...")
     print("(This solves SS + TPI for both scenarios — may take a while.)")
     client = Client(n_workers=2, threads_per_worker=1, memory_limit="2GB")
     t0 = time.time()
     try:
         base_tp, reform_tp = run_transition_path(
-            start_year=2026, policy=REFORM, client=client
+            start_year=2026,
+            policy=REFORM,
+            client=client,
+            multi_sector=multi_sector,
         )
     finally:
         client.close()
@@ -137,12 +151,13 @@ def run_tpi():
 def main():
     mode = sys.argv[1] if len(sys.argv) > 1 else "ss"
     age_specific = sys.argv[2] if len(sys.argv) > 2 else "pooled"
+    multi_sector = "multi-sector" in sys.argv
     if mode == "ss":
-        run_steady_state(age_specific=age_specific)
+        run_steady_state(age_specific=age_specific, multi_sector=multi_sector)
     elif mode == "tpi":
-        run_tpi()
+        run_tpi(multi_sector=multi_sector)
     else:
-        print(f"Usage: {sys.argv[0]} [ss|tpi] [pooled|brackets|each]")
+        print(f"Usage: {sys.argv[0]} [ss|tpi] [pooled|brackets|each] [multi-sector]")
         sys.exit(1)
 
 

oguk/api.py

@@ -127,7 +127,7 @@ class Config:
     tax_revenue: np.ndarray = Field(description="Total tax revenue")
     debt: np.ndarray = Field(description="Government debt")
     g_y_annual: float = Field(
-        default=0.010,
+        default=0.011,
         description="Annual productivity growth rate used in this run",
     )
     # Per-industry arrays (T × M)
@@ -696,16 +696,47 @@ def _build_specs(
     max_iter: int = 250,
     age_specific: str = "pooled",
     param_overrides: dict | None = None,
+    multi_sector: bool = False,
 ):
     """Build a calibrated Specifications object (internal).
 
+    UK-calibrated structural parameters (set in oguk_default_parameters.json):
+
+        frisch = 0.35
+            Blundell, MaCurdy & Meghir (1999) "Labor Supply Models: Unobserved
+            Heterogeneity, Nonparticipation and Dynamics", ch. 27 in Handbook
+            of Labor Economics. Central UK estimate ~0.3-0.4.
+            https://doi.org/10.1016/S1573-4463(99)30039-9
+
+        g_y_annual = 0.011
+            OBR "Economic and Fiscal Outlook" (Nov 2025), Table 1.1: potential
+            output growth ~1.1% per year.
+            https://obr.uk/efo/economic-and-fiscal-outlook-november-2025/
+
+        delta_annual = 0.065
+            ONS "Capital stocks and fixed capital consumption" (CAPSTK):
+            ratio of consumption of fixed capital to net capital stock, whole
+            economy excluding dwellings, 2022 (~6-7%).
+            https://www.ons.gov.uk/economy/nationalaccounts/uksectoraccounts/datasets/capitalstocksandfixedcapitalconsumption
+
+        beta_annual = 0.965
+            Calibrated to match UK household saving ratio (~6-9%).
+            ONS "Households saving ratio" (NRJS).
+            https://www.ons.gov.uk/economy/grossdomesticproductgdp/timeseries/dgd8
+
+        world_int_rate_annual = 0.02
+            UK 10-year gilt real yield, Bank of England yield curve data.
+            https://www.bankofengland.co.uk/statistics/yield-curves
+
     Args:
         param_overrides: Optional dict of OG-Core parameter names to values
             that are applied after all other configuration. Use this to shock
             structural parameters (e.g. ``{"g_y_annual": 0.011}``). Overrides
             are applied last so they take precedence over defaults and
             calibration outputs. Demographic / tax-function keys are not
             permitted here — those are set by calibrate().
+        multi_sector: If True, use 8-sector industry calibration (M=8).
+            If False (default), use a single representative sector (M=1).
     """
     from ogcore.parameters import Specifications
 
@@ -738,20 +769,21 @@ def _build_specs(
     ]:
         defaults.pop(key, None)
 
-    # Strip single-industry defaults that will be replaced by industry_params
-    for key in [
-        "gamma",
-        "gamma_g",
-        "epsilon",
-        "Z",
-        "cit_rate",
-        "io_matrix",
-        "alpha_c",
-        "delta_tau_annual",
-        "inv_tax_credit",
-        "tau_c",
-    ]:
-        defaults.pop(key, None)
+    if multi_sector:
+        # Strip single-industry defaults that will be replaced by industry_params
+        for key in [
+            "gamma",
+            "gamma_g",
+            "epsilon",
+            "Z",
+            "cit_rate",
+            "io_matrix",
+            "alpha_c",
+            "delta_tau_annual",
+            "inv_tax_credit",
+            "tau_c",
+        ]:
+            defaults.pop(key, None)
 
     p = Specifications(
         baseline=baseline, output_base=output_base, baseline_dir=baseline_dir
@@ -789,10 +821,14 @@ def _build_specs(
         }
     )
 
-    # Apply 8-sector industry calibration
-    defaults.update(get_industry_params())
-    # Levenberg-Marquardt is more robust than the default 'hybr' for M>1
-    defaults["SS_root_method"] = "lm"
+    if multi_sector:
+        # Apply 8-sector industry calibration
+        defaults.update(get_industry_params())
+        # hybr (Powell hybrid) for heterogeneous CES; LM gets stuck at ~1e-5
+        defaults["SS_root_method"] = "hybr"
+        # Relax tolerances for heterogeneous CES
+        defaults["mindist_SS"] = 1e-4
+        defaults["RC_SS"] = 1e-4
 
     if param_overrides:
         defaults.update(param_overrides)
@@ -850,6 +886,7 @@ def solve_steady_state(
     max_iter: int = 250,
     age_specific: str = "pooled",
     param_overrides: dict | None = None,
+    multi_sector: bool = False,
 ) -> SteadyStateResult:
     """Solve for steady state equilibrium.
 
@@ -863,6 +900,8 @@ def solve_steady_state(
             "each"     — separate function per individual age (80)
         param_overrides: Optional dict of OG-Core parameter names to values
             for structural shocks (e.g. ``{"g_y_annual": 0.011}``).
+        multi_sector: If True, use 8-sector industry calibration (M=8).
+            If False (default), use a single representative sector (M=1).
 
     Returns:
         SteadyStateResult with equilibrium values
@@ -878,6 +917,7 @@ def solve_steady_state(
             max_iter=max_iter,
             age_specific=age_specific,
             param_overrides=param_overrides,
+            multi_sector=multi_sector,
         )
         ss = run_SS(p, client=None)
         return _ss_dict_to_result(ss)
@@ -889,6 +929,7 @@ def run_transition_path(
     client=None,
     age_specific: str = "pooled",
     param_overrides: dict | None = None,
+    multi_sector: bool = False,
 ) -> tuple[TransitionPathResult, TransitionPathResult | None]:
     """Run baseline (and optionally reform) transition path.
 
@@ -906,6 +947,8 @@ def run_transition_path(
             "each"     — separate function per individual age (80)
         param_overrides: Optional dict of OG-Core parameter names to values
             for structural shocks (e.g. ``{"Z": [[1.004]]}``).
+        multi_sector: If True, use 8-sector industry calibration (M=8).
+            If False (default), use a single representative sector (M=1).
 
     Returns:
         (baseline_tp, reform_tp) — reform_tp is None if no policy/overrides
@@ -930,6 +973,7 @@ def run_transition_path(
             base_dir,
             baseline=True,
             age_specific=age_specific,
+            multi_sector=multi_sector,
         )
 
         # Solve SS first to auto-calibrate alpha_G.
@@ -963,6 +1007,7 @@ def run_transition_path(
                 baseline=False,
                 age_specific=age_specific,
                 param_overrides=param_overrides,
+                multi_sector=multi_sector,
             )
 
             ss_reform = SS.run_SS(p_reform, client=client)

oguk/industry_params.py

@@ -256,14 +256,36 @@ def _sector_gva() -> np.ndarray:
     0.0,  # Manufacturing
 ]
 
+# --- Annual depreciation rates by sector ---
+# Source: ONS "Capital stocks and fixed capital consumption" (CAPSTK),
+# ratio of consumption of fixed capital to net capital stock by SIC 2007
+# section, 2022.
+# https://www.ons.gov.uk/economy/nationalaccounts/uksectoraccounts/datasets/capitalstocksandfixedcapitalconsumption
+# Energy and manufacturing have higher depreciation (machinery-heavy).
+# Real estate has lower depreciation (long-lived structures dominate).
+_DELTA_ANNUAL = [
+    0.08,  # Energy       — rigs, turbines depreciate fast
+    0.06,  # Construction — equipment + vehicles
+    0.07,  # Trade & Transport — vehicles, fit-outs
+    0.10,  # Info & Finance — IT equipment, short-lived assets
+    0.03,  # Real Estate  — dominated by long-lived structures
+    0.08,  # Business Svcs— IT equipment, office fit-outs
+    0.05,  # Public & Other — mix of long/short-lived assets
+    0.09,  # Manufacturing— machinery, plant, equipment
+]
 
-def _sector_tfp() -> list:
+
+def _sector_tfp(epsilon=None, gamma=None) -> list:
     """Solow-residual TFP by sector, normalised so the GVA-weighted mean = 1.
 
-    Computes Z_m = GVA_m / (K_m^gamma_m * L_m^(1 - gamma_m)) for each sector,
-    then rescales so that sum(gva_share_m * Z_m) = 1.0.  This ensures the
-    aggregate economy matches the baseline calibration while allowing
-    cross-sector productivity differences.
+    When epsilon is all 1.0 (Cobb-Douglas), computes:
+        Z_m = GVA_m / (K_m^gamma_m * L_m^(1 - gamma_m))
+
+    When epsilon differs from 1.0 (CES), computes the CES residual:
+        Z_m = GVA_m / [gamma_m^(1/eps) * K_m^((eps-1)/eps)
+               + (1-gamma_m)^(1/eps) * L_m^((eps-1)/eps)]^(eps/(eps-1))
+
+    Then rescales so that sum(gva_share_m * Z_m) = 1.0.
 
     Sources:
         GVA: _GVA_BY_SIC_SECTION (ONS Blue Book 2024)
@@ -273,9 +295,31 @@ def _sector_tfp() -> list:
     gva = _sector_gva()
     capital = np.array(_CAPITAL_STOCK, dtype=float)
     labour = np.array(_WORKFORCE_JOBS, dtype=float)
-    gamma = np.array(_GAMMA, dtype=float)
-
-    z_raw = gva / (capital**gamma * labour ** (1 - gamma))
+    if gamma is None:
+        gamma = np.array(_GAMMA, dtype=float)
+    else:
+        gamma = np.array(gamma, dtype=float)
+    if epsilon is None:
+        epsilon = np.ones(M, dtype=float)
+    else:
+        epsilon = np.array(epsilon, dtype=float)
+
+    z_raw = np.empty(M, dtype=float)
+    for m in range(M):
+        eps = epsilon[m]
+        g = gamma[m]
+        K = capital[m]
+        L = labour[m]
+        if eps == 1.0:
+            # Cobb-Douglas
+            z_raw[m] = gva[m] / (K**g * L ** (1 - g))
+        else:
+            # CES: Y = Z * [gamma^(1/eps)*K^((eps-1)/eps) + (1-gamma)^(1/eps)*L^((eps-1)/eps)]^(eps/(eps-1))
+            ces_aggregate = (
+                g ** (1 / eps) * K ** ((eps - 1) / eps)
+                + (1 - g) ** (1 / eps) * L ** ((eps - 1) / eps)
+            ) ** (eps / (eps - 1))
+            z_raw[m] = gva[m] / ces_aggregate
 
     # Normalise so GVA-weighted mean equals 1
     gva_shares = gva / gva.sum()
@@ -291,14 +335,12 @@ def get_industry_params() -> dict:
     All arrays are in list form for JSON compatibility.
 
     Gamma is shrunk toward the aggregate UK mean for solver stability.
-    Epsilon is set to 1.0 (Cobb-Douglas) because OG-Core's TPI solver
-    does not converge with heterogeneous CES elasticities (produces NaN
-    at iteration 1). The raw literature values are preserved in _EPSILON
-    for future use when OG-Core TPI is fixed.
+    Epsilon uses calibrated CES elasticities from the literature.
+    Depreciation rates are sector-specific from ONS capital stock data.
 
     Shrinkage:
       - gamma: 40% shrinkage toward 0.35 (aggregate UK capital share)
-      - epsilon: Cobb-Douglas (1.0) for all sectors (TPI stability)
+      - epsilon: 50% shrinkage toward 1.0 (Cobb-Douglas) applied in _EPSILON values
     """
     gva = _sector_gva()
     gva_shares = gva / gva.sum()
@@ -309,18 +351,20 @@ def get_industry_params() -> dict:
     # Shrink gamma 40% toward aggregate mean (0.35) for solver stability
     gamma_shrunk = [0.35 + 0.6 * (g - 0.35) for g in _GAMMA]
 
+    epsilon = list(_EPSILON)
+
     return {
         "M": M,
         "I": NUM_CONSUMPTION_GOODS,
         "gamma": gamma_shrunk,
         "gamma_g": [0.0] * M,
-        "epsilon": [1.0] * M,  # Cobb-Douglas; heterogeneous epsilon breaks TPI
-        "Z": [_sector_tfp()],
+        "epsilon": epsilon,  # calibrated CES elasticities by sector
+        "Z": [_sector_tfp(epsilon=epsilon, gamma=gamma_shrunk)],
         "cit_rate": [[0.27] * M],
         "io_matrix": _IO_MATRIX,
         "alpha_c": alpha_c.tolist(),
         "c_min": list(_C_MIN),
-        "delta_tau_annual": [[0.05] * M],
+        "delta_tau_annual": [list(_DELTA_ANNUAL)],  # sector-specific tax depreciation
         "inv_tax_credit": [[0.0] * M],
         "tau_c": [[0.19] * NUM_CONSUMPTION_GOODS],
     }