#!/usr/bin/env python3 """ Baumol Convergence Test -- Paper 4, Hierarchical Ceiling Cascade ================================================================= Paper 4 (Hierarchical CES Architecture), Section "The Hierarchical Ceiling Cascade" (Proposition 5.1). **Prediction**: Each level's output is bounded by the level below it in the timescale hierarchy: F_n <= phi_n(F_{n-1})/(sigma_n * J). The long-run growth rate equals the improvement rate of the slowest sector. Baumol cost disease = CES complementarity constraint. **Lean**: ceiling_cascade (TransitionDynamics.lean), ceiling_bound, ceiling_nonneg. Sub-tests: 4a: Var(g_10yr)/Var(g_5yr) < 0.5 (i.i.d. null, not just < 1.0) 4b: Supply chain growth gaps shrink from 5yr to 10yr horizons 4c: Bottleneck upstream dominates downstream growth (Leontief-based) 4d: Kendall tau(growth_rank, upstreamness) < 0 Data: BEA real value added quarterly, BEA Use table (supply chain), BEA Leontief inverse (upstreamness). Review Notes (2026-03-12): - Theory match: CORRECT (Paper 4 Prop 5.1, ceiling cascade) - I/O audit: CORRECT — growth from real VA (output), upstreamness from Leontief - Data-signal alignment: ADEQUATE (64 leaf industries, 2005-2019) - Code quality: BUGS FOUND and FIXED (see Change Log) - Confidence: MEDIUM (3/4 CONSISTENT after fixes, bottleneck test redesigned) - Known limitations: 2005-2009 period includes GFC (inflates 5yr variance); bottleneck test uses Leontief inverse (industry-by-industry) instead of Use table (commodity-by-industry) to fix code mismatch. """ import warnings import numpy as np import pandas as pd from scipy import stats from common.config import COMPILED_DIR, RESULTS_DIR from common.plotting import setup_figure, setup_subplots, save_figure from common.stats import verdict, log_result warnings.filterwarnings('ignore') GOVT_HOUSING = {'GFE', 'GFGD', 'GFGN', 'GSLE', 'GSLG', 'HS', 'ORE'} VA_AGGREGATES = { 'GDP', 'PVT', 'G', 'GF', 'GFG', 'GSL', '11', '21', '31G', '31ND', '33DG', '44RT', '48TW', '51', '52', '53', '54', '56', '6', '62', '7', '71', '72', 'FIRE', 'NABI', 'PGOOD', 'PROF', 'PSERV', '531', } _ALL_SUBTESTS = [ 'baumol_convergence_dispersion_shrinks', 'baumol_convergence_supply_chain_gap', 'baumol_convergence_bottleneck_dominates', 'baumol_convergence_upstream_slower', ] def compute_growth(df, start_year, end_year): """Compute annualized log-growth between two periods.""" horizon = end_year - start_year start_vals = df[df['Year'] == start_year].groupby('Industry')['val'].mean() end_vals = df[df['Year'] == end_year].groupby('Industry')['val'].mean() common = start_vals.index.intersection(end_vals.index) s = start_vals[common]; e = end_vals[common] valid = (s > 0) & (e > 0) return (np.log(e[valid]) - np.log(s[valid])) / horizon def pair_gaps(g, pairs): """Compute |g_downstream - g_upstream| for each supply chain pair.""" gaps = [] for _, row in pairs.iterrows(): u, d = row['upstream'], row['downstream'] if u in g.index and d in g.index: gaps.append(abs(g[d] - g[u])) return np.array(gaps) _cache = {} def main(): results = [] # Load VA data try: va = pd.read_csv(COMPILED_DIR / 'bea_real_va_quarterly.csv') except FileNotFoundError: try: va = pd.read_csv(COMPILED_DIR.parent / 'bea_gdp_industry' / 'RealValueAdded_10_quarterly.csv') except FileNotFoundError: return [log_result('baumol_convergence_dispersion_shrinks', np.nan, np.nan, 'SKIPPED', details={'reason': 'VA data not found'}, description="Tests Baumol convergence: growth rate dispersion shrinks at longer horizons.")] va['val'] = pd.to_numeric(va.get('DataValue', va.get('val', pd.Series())), errors='coerce') if 'Year' not in va.columns and 'year' in va.columns: va = va.rename(columns={'year': 'Year'}) if 'Industry' not in va.columns and 'industry' in va.columns: va = va.rename(columns={'industry': 'Industry'}) # Load Use table for supply chain pairs try: use = pd.read_csv(COMPILED_DIR / 'bea_use_summary.csv') except FileNotFoundError: try: use = pd.read_csv(COMPILED_DIR.parent / 'bea_io' / 'Use_Summary_259_2017.csv') except FileNotFoundError: use = None # Load Leontief inverse for upstreamness and bottleneck try: L_mat = pd.read_csv(COMPILED_DIR / 'bea_leontief_full.csv', index_col=0) except FileNotFoundError: try: L_mat = pd.read_csv(COMPILED_DIR.parent / 'bea_leontief_full.csv', index_col=0) except FileNotFoundError: L_mat = None # Filter to leaf industries if use is not None: use['val'] = pd.to_numeric(use.get('DataValue', use.get('val', pd.Series())), errors='coerce') if 'RowCode' in use.columns and 'ColCode' in use.columns: use_ind = use[(use.get('RowType') == 'Commodity') & (use.get('ColType') == 'Industry')].copy() use_ind = use_ind[~use_ind['RowCode'].isin(GOVT_HOUSING) & ~use_ind['ColCode'].isin(GOVT_HOUSING)] use_industry_codes = set(use_ind['ColCode'].unique()) | set(use_ind['RowCode'].unique()) use_industry_codes -= GOVT_HOUSING else: use_ind = None use_industry_codes = set() else: use_ind = None use_industry_codes = set() leaf_codes = use_industry_codes - VA_AGGREGATES if use_industry_codes else set(va['Industry'].unique()) va_leaf = va[va['Industry'].isin(leaf_codes)].copy() # Compute growth rates g_5yr_early = compute_growth(va_leaf, 2005, 2009) g_5yr_late = compute_growth(va_leaf, 2015, 2019) g_10yr_A = compute_growth(va_leaf, 2005, 2014) g_10yr_B = compute_growth(va_leaf, 2010, 2019) g_15yr = compute_growth(va_leaf, 2005, 2019) # ── 4a: Var(g_10yr) < Var(g_5yr) with bootstrap p-value ── var_5yr_avg = (g_5yr_early.var() + g_5yr_late.var()) / 2 var_10yr_avg = (g_10yr_A.var() + g_10yr_B.var()) / 2 if var_5yr_avg > 0: ratio_10_5 = var_10yr_avg / var_5yr_avg # Under i.i.d. annual shocks, Var(annualized_g_T) = sigma^2/T # so the expected ratio Var(g_10)/Var(g_5) = 5/10 = 0.5 # Test: is observed ratio significantly below 0.5? iid_null = 0.5 rng = np.random.default_rng(42) n_boot = 2000 # Bootstrap: resample industries, compute ratio each time industries = list(set(g_5yr_early.index) & set(g_5yr_late.index) & set(g_10yr_A.index) & set(g_10yr_B.index)) boot_ratios = [] for _ in range(n_boot): idx = rng.choice(industries, size=len(industries), replace=True) v5 = (g_5yr_early[idx].var() + g_5yr_late[idx].var()) / 2 v10 = (g_10yr_A[idx].var() + g_10yr_B[idx].var()) / 2 if v5 > 0: boot_ratios.append(v10 / v5) boot_ratios = np.array(boot_ratios) # p-value: fraction of bootstrap ratios >= i.i.d. null (0.5) p_disp = float(np.mean(boot_ratios >= iid_null)) if len(boot_ratios) > 0 else np.nan p_disp = max(p_disp, 1.0 / n_boot) # floor at 1/n_boot sign_correct = ratio_10_5 < iid_null verd = verdict(p_disp, sign_correct) _cache['dispersion'] = { 'ratio': ratio_10_5, 'p': p_disp, 'iid_null': iid_null, 'var_5yr': var_5yr_avg, 'var_10yr': var_10yr_avg, 'boot_ratios': boot_ratios, 'n_industries': len(industries), } results.append(log_result('baumol_convergence_dispersion_shrinks', ratio_10_5, p_disp, verd, details={'var_5yr': var_5yr_avg, 'var_10yr': var_10yr_avg, 'iid_null': iid_null, 'n_industries': len(industries)}, description="Tests that cross-industry growth rate variance shrinks FASTER than " "mechanical time-averaging (i.i.d. null ratio = 0.5). Bootstrap (2000 resamples) " "of Var(g_10yr)/Var(g_5yr) across BEA leaf industries. Ratio < 0.5 indicates " "genuine Baumol convergence beyond mechanical averaging.")) else: results.append(log_result('baumol_convergence_dispersion_shrinks', np.nan, np.nan, 'SKIPPED', description="Tests that cross-industry growth rate variance shrinks faster than i.i.d. null.")) # Build supply chain pairs from Use table if use_ind is not None and 'RowCode' in use_ind.columns: total_inputs = use_ind.groupby('ColCode')['val'].sum() use_pairs = use_ind.groupby(['RowCode', 'ColCode'])['val'].sum().reset_index() use_pairs = use_pairs.merge(total_inputs.rename('total_input'), left_on='ColCode', right_index=True) use_pairs['share'] = use_pairs['val'] / use_pairs['total_input'] supply_pairs = use_pairs[(use_pairs['share'] > 0.05) & (use_pairs['RowCode'] != use_pairs['ColCode'])].copy() supply_pairs = supply_pairs.rename(columns={'RowCode': 'upstream', 'ColCode': 'downstream'}) else: supply_pairs = None # ── 4b: Supply chain gap shrinks ── if supply_pairs is not None and len(supply_pairs) > 0: gaps_5_e = pair_gaps(g_5yr_early, supply_pairs) gaps_5_l = pair_gaps(g_5yr_late, supply_pairs) gaps_10_A = pair_gaps(g_10yr_A, supply_pairs) gaps_10_B = pair_gaps(g_10yr_B, supply_pairs) gaps_5yr = (gaps_5_e.mean() + gaps_5_l.mean()) / 2 if len(gaps_5_e) > 0 and len(gaps_5_l) > 0 else 0 gaps_10yr = (gaps_10_A.mean() + gaps_10_B.mean()) / 2 if len(gaps_10_A) > 0 and len(gaps_10_B) > 0 else 0 if gaps_5yr > 0: gap_ratio = gaps_10yr / gaps_5yr # Permutation test: shuffle horizon labels rng = np.random.default_rng(43) n_perm = 2000 perm_ratios = [] # Pool all gap arrays and randomly assign to "5yr" and "10yr" all_gaps = np.concatenate([gaps_5_e, gaps_5_l, gaps_10_A, gaps_10_B]) n_5e, n_5l = len(gaps_5_e), len(gaps_5_l) n_10A, n_10B = len(gaps_10_A), len(gaps_10_B) for _ in range(n_perm): perm = rng.permutation(all_gaps) p5 = (perm[:n_5e].mean() + perm[n_5e:n_5e+n_5l].mean()) / 2 p10 = (perm[n_5e+n_5l:n_5e+n_5l+n_10A].mean() + perm[n_5e+n_5l+n_10A:].mean()) / 2 if p5 > 0: perm_ratios.append(p10 / p5) p_gap = float(np.mean(np.array(perm_ratios) <= gap_ratio)) if perm_ratios else np.nan p_gap = max(p_gap, 1.0 / n_perm) sign_correct = gap_ratio < 1.0 verd = verdict(p_gap, sign_correct) _cache['supply_chain'] = { 'ratio': gap_ratio, 'p': p_gap, 'gaps_5yr': gaps_5yr, 'gaps_10yr': gaps_10yr, } results.append(log_result('baumol_convergence_supply_chain_gap', gap_ratio, p_gap, verd, details={'gaps_5yr': gaps_5yr, 'gaps_10yr': gaps_10yr}, description="Tests that supply chain growth gaps shrink at longer horizons. " "Computes mean |g_downstream - g_upstream| for BEA Use Table pairs at 5yr " "and 10yr horizons. Permutation test (2000 shuffles) on gap ratio; " "ratio < 1 means convergence.")) else: results.append(log_result('baumol_convergence_supply_chain_gap', np.nan, np.nan, 'SKIPPED', description="Tests supply chain growth gap convergence.")) else: results.append(log_result('baumol_convergence_supply_chain_gap', np.nan, np.nan, 'SKIPPED', description="Tests supply chain growth gap convergence.")) # ── 4c: Bottleneck test (using Leontief inverse, not Use table) ── # Use Leontief inverse (industry-by-industry) for supply chain relationships # because Use table mixes commodity/industry codes causing poor overlap with VA data. if L_mat is not None and len(g_15yr) >= 5: L_cols = [c for c in L_mat.columns if c not in GOVT_HOUSING] L_rows = [r for r in L_mat.index if r not in GOVT_HOUSING] L_clean = L_mat.loc[L_rows, L_cols] # For each downstream industry d, find its top upstream suppliers # (industries with largest L[upstream, downstream] off-diagonal entries) common_industries = list(set(g_15yr.index) & set(L_clean.columns) & set(L_clean.index)) if len(common_industries) >= 5: g_min_leontief = {} g_avg_leontief = {} for d in common_industries: if d not in L_clean.columns: continue col = L_clean[d].drop(d, errors='ignore') # Top upstream suppliers (Leontief coefficient > 0.02) top_ups = col[col > 0.02].index ups_with_growth = [u for u in top_ups if u in g_15yr.index] if len(ups_with_growth) >= 2: # need at least 2 for min != avg g_min_leontief[d] = g_15yr[ups_with_growth].min() g_avg_leontief[d] = g_15yr[ups_with_growth].mean() if len(g_min_leontief) >= 5: g_min_s = pd.Series(g_min_leontief) g_avg_s = pd.Series(g_avg_leontief) common_d = g_15yr.index.intersection(g_min_s.index) r_min, p_min = stats.spearmanr(g_15yr[common_d], g_min_s[common_d]) r_avg, _ = stats.spearmanr(g_15yr[common_d], g_avg_s[common_d]) pass_4c = (r_min > 0) and (r_min > r_avg) verd = 'CONSISTENT' if pass_4c else 'AMBIGUOUS' _cache['bottleneck'] = { 'r_min': r_min, 'r_avg': r_avg, 'p_min': p_min, 'n_industries': len(common_d), 'g_down': g_15yr[common_d].values, 'g_min_up': g_min_s[common_d].values, 'g_avg_up': g_avg_s[common_d].values, } results.append(log_result('baumol_convergence_bottleneck_dominates', r_min, p_min, verd, details={'r_min': r_min, 'r_avg': r_avg, 'n_industries': len(common_d)}, description="Tests that slowest (bottleneck) upstream industry dominates " "downstream growth more than average upstream. Uses Leontief inverse " "(industry-by-industry) to identify upstream suppliers with coefficient > 0.02. " "Spearman r(downstream growth, min-upstream growth) vs r(downstream, avg-upstream). " "Prediction: r_min > r_avg > 0.")) else: results.append(log_result('baumol_convergence_bottleneck_dominates', np.nan, np.nan, 'SKIPPED', description="Bottleneck test requires >= 5 downstream industries with >= 2 upstream suppliers.")) else: results.append(log_result('baumol_convergence_bottleneck_dominates', np.nan, np.nan, 'SKIPPED', description="Bottleneck test skipped: insufficient industry overlap.")) else: results.append(log_result('baumol_convergence_bottleneck_dominates', np.nan, np.nan, 'SKIPPED', description="Bottleneck test requires Leontief inverse data.")) # ── 4d: Kendall tau(growth, upstreamness) < 0 ── if L_mat is not None: L_cols = [c for c in L_mat.columns if c not in GOVT_HOUSING] upstreamness = L_mat[L_cols].sum(axis=0) common = g_15yr.index.intersection(upstreamness.index) if len(common) >= 5: tau, p = stats.kendalltau(g_15yr[common], upstreamness[common]) p_one = p / 2 if tau < 0 else 1 - p / 2 pass_4d = (tau < 0) and (p_one < 0.10) verd = 'CONSISTENT' if pass_4d else 'AMBIGUOUS' _cache['upstream'] = { 'tau': tau, 'p': p_one, 'growth': g_15yr[common].values, 'upstreamness': upstreamness[common].values, 'names': list(common), 'n': len(common), } results.append(log_result('baumol_convergence_upstream_slower', tau, p_one, verd, details={'n_industries': len(common)}, description="Tests that more upstream industries grow slower, as predicted by " "hierarchical ceiling cascade (Paper 4, Prop 5.1). Kendall tau between " "15-year annualized log-growth and Leontief-inverse upstreamness (forward/" "backward linkage ratio) across BEA industries. Prediction: tau < 0.")) else: results.append(log_result('baumol_convergence_upstream_slower', np.nan, np.nan, 'SKIPPED', description="Upstream slower test skipped: insufficient industry overlap.")) else: results.append(log_result('baumol_convergence_upstream_slower', np.nan, np.nan, 'SKIPPED', description="Upstream slower test requires Leontief inverse data.")) _cache['results'] = results return results def figures(): """Generate per-subtest figures.""" if not _cache: main() paths = [] # --- Dispersion shrinks --- data = _cache.get('dispersion') if data and len(data.get('boot_ratios', [])) > 0: fig, ax = setup_figure() ax.hist(data['boot_ratios'], bins=40, alpha=0.7, color='steelblue', edgecolor='k') ax.axvline(data['ratio'], color='red', lw=2, label=f'Observed = {data["ratio"]:.3f}') ax.axvline(data['iid_null'], color='orange', lw=2, ls='--', label=f'i.i.d. null = {data["iid_null"]:.1f}') ax.set_xlabel('Var(g_10yr) / Var(g_5yr)') ax.set_ylabel('Bootstrap count') ax.set_title(f'Dispersion Convergence (p={data["p"]:.4f})') ax.legend(fontsize=9) paths.extend(save_figure(fig, 'baumol_dispersion_shrinks', subfolder='applications')) # --- Supply chain gap --- data = _cache.get('supply_chain') if data: fig, ax = setup_figure() ax.bar(['5-year\nhorizon', '10-year\nhorizon'], [data['gaps_5yr'], data['gaps_10yr']], color=['coral', 'steelblue'], edgecolor='k') ax.set_ylabel('Mean |g_downstream - g_upstream|') ax.set_title(f'Supply Chain Gap Convergence (ratio={data["ratio"]:.3f}, p={data["p"]:.4f})') ax.grid(True, alpha=0.3, axis='y') paths.extend(save_figure(fig, 'baumol_supply_chain_gap', subfolder='applications')) # --- Bottleneck --- data = _cache.get('bottleneck') if data: fig, axes = setup_subplots(1, 2, width=12, height=5) ax = axes[0] ax.scatter(data['g_min_up'], data['g_down'], s=40, color='crimson', edgecolor='k', alpha=0.7, label=f'min upstream (r={data["r_min"]:.3f})') ax.set_xlabel('Min upstream growth') ax.set_ylabel('Downstream growth') ax.set_title('Bottleneck (min upstream)') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) ax = axes[1] ax.scatter(data['g_avg_up'], data['g_down'], s=40, color='steelblue', edgecolor='k', alpha=0.7, label=f'avg upstream (r={data["r_avg"]:.3f})') ax.set_xlabel('Avg upstream growth') ax.set_ylabel('Downstream growth') ax.set_title('Average upstream') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) fig.suptitle(f'Bottleneck Dominance (N={data["n_industries"]})', fontsize=12) paths.extend(save_figure(fig, 'baumol_bottleneck', subfolder='applications')) # --- Upstream slower --- data = _cache.get('upstream') if data: fig, ax = setup_figure() ax.scatter(data['upstreamness'], data['growth'], s=40, color='seagreen', edgecolor='k', alpha=0.7) ax.set_xlabel('Upstreamness (forward/backward linkage)') ax.set_ylabel('15-year annualized log growth') ax.set_title(f'Upstream Slower (tau={data["tau"]:.3f}, p={data["p"]:.4f})') ax.grid(True, alpha=0.3) # Trend line z = np.polyfit(data['upstreamness'], data['growth'], 1) x_line = np.linspace(min(data['upstreamness']), max(data['upstreamness']), 100) ax.plot(x_line, np.polyval(z, x_line), 'r--', alpha=0.5) paths.extend(save_figure(fig, 'baumol_upstream_slower', subfolder='applications')) return paths if __name__ == '__main__': results = main() for r in results: print(f"{r['name']}: {r['verdict']} " f"(stat={r['statistic']}, p={r['p_value']})")

p value	4.232e-4
statistic	-0.2857

Baumol Convergence Upstream Slower

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources