#!/usr/bin/env python3 """ Empirical Tests for Sectoral Ramsey / Inverse-Elasticity Rule ============================================================= Lean: CESProofs/Macro/GrowthTax.lean, ramsey_inverse_elasticity; CESProofs/Macro/TaxStructure.lean, saezTopRate_decreasing_in_e. Tests five properties using US manufacturing data: 1. Sector elasticity aggregation (SITC3 → NAICS via 2-digit prefix) 2. Effective tax rate by sector (Taxes/VA from BEA composition) 3. Ramsey inverse-elasticity (Spearman rho of tau vs sigma; prediction: rho < 0) 4. Recession ordering (output drops correlated with sigma) 5. Revenue decomposition (low-sigma vs high-sigma effective tax rates) Data: broda_weinstein_us_sigmas.csv, bea_composition.csv Review Notes (2026-03-12) ------------------------- Theory match: CORRECT (with sign ambiguity). Lean theorems ramsey_inverse_elasticity and saezTopRate_decreasing_in_e are REAL PROOFS (not schematic). Both prove τ*(e₂) < τ*(e₁) when e₂ > e₁. CRITICAL DISTINCTION: This is a NORMATIVE prediction (what taxes SHOULD be under Ramsey optimality), not a POSITIVE prediction (what taxes ARE). US manufacturing tax policy is set by political economy, historical accident, and industry lobbying — not by optimal tax theory. An INCONSISTENT result means the government doesn't follow the Ramsey rule, NOT that the theorem is wrong. Tests 3-5 are capped at AMBIGUOUS when statistically insignificant. SIGN AMBIGUITY: The mapping from CES σ to behavioral elasticity e is non-trivial. Demand-side (BW trade σ): higher σ → more elastic → negative σ-τ prediction. CES supply-side (curvature K ∝ 1/σ): higher K → more distortion → positive σ-τ prediction. Paper 2 Remark rem:ramsey-ces has a text inconsistency: says "lower behavioral elasticity" but the logic requires "higher effective behavioral elasticity." The observed ρ_s=+0.300 is positive (consistent with supply-side), but completely insignificant at N=9. I/O audit: MIXED — both key proxies (trade σ for behavioral e, BEA tax/VA for optimal τ) are WEAK. Trade sigma measures substitutability between imported varieties; e measures responsiveness of the tax base to rate changes. BEA tax/VA includes property taxes and subsidies, not narrowly-defined Ramsey commodity taxes. Data-signal alignment: WEAK. N=9 manufacturing sectors (critical |rho_s| ~ 0.643 for p<0.05). Recession ordering confounded by durability. Capital intensity confounds tax rates (primary metals highest at 3.99% due to property taxes on mills). Code quality: CLEAN (after 2026-03-11 SITC concordance fix). Confidence in verdict: LOW. Normative prediction with weak proxies and N=9. Labor Elasticity tests (5/5 CONSISTENT) provide much stronger evidence for the same Lean theorem. Date of review: 2026-03-12 """ import numpy as np from scipy import stats as sp_stats import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.data import load_compiled from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure # SITC2 → NAICS concordance # Corrected 2026-03-11: removed SITC 61 (leather→NAICS 316), 62 (rubber→326), # 66 (non-metallic minerals→327) — these map to NAICS sectors outside our 9-sector # panel and were contaminating NAICS 331 (Primary metals). Fixed SITC 68 # (non-ferrous metals) from 332→331. SITC_TO_NAICS = { "63": "337", "64": "339", "65": "339", "67": "331", "68": "331", "69": "332", "71": "333", "72": "333", "73": "333", "74": "333", "75": "334", "76": "334", "77": "335", "78": "3361MV", "79": "3364OT", "81": "337", "82": "337", "83": "339", "84": "339", "85": "339", "87": "334", "88": "334", "89": "339", } NAICS_NAMES = { "331": "Primary metals", "332": "Fabricated metals", "333": "Machinery", "334": "Computer & electronic", "335": "Electrical equipment", "3361MV": "Motor vehicles", "3364OT": "Other transportation", "337": "Furniture", "339": "Miscellaneous mfg", } def _spearman(x, y): """Compute Spearman rank correlation with p-value via scipy.""" rho, p = sp_stats.spearmanr(x, y) return rho, p def _permutation_p(x, y, observed_rho, n_perm=5000, seed=45, alternative='less'): """Permutation test for Spearman correlation. alternative='less': p = fraction of permuted rho <= observed (for rho < 0 prediction) alternative='greater': p = fraction of permuted rho >= observed (for rho > 0 prediction) """ np.random.seed(seed) x_arr = np.array(x) count = 0 for _ in range(n_perm): perm_x = np.random.permutation(x_arr) perm_rho, _ = sp_stats.spearmanr(perm_x, y) if alternative == 'less': if perm_rho <= observed_rho: count += 1 else: if perm_rho >= observed_rho: count += 1 return count / n_perm def _aggregate_sigmas(bw): """Aggregate Broda-Weinstein SITC3 → NAICS sector median sigmas.""" sector_sigmas = {naics: [] for naics in NAICS_NAMES} for _, row in bw.iterrows(): raw = row['sitc3_code'] sitc_code = str(int(raw)) if isinstance(raw, float) else str(raw).strip() sigma = float(row['sigma']) if sigma <= 0 or sigma > 100: continue prefix = sitc_code[:2] if prefix in SITC_TO_NAICS: sector_sigmas[SITC_TO_NAICS[prefix]].append(sigma) result = {} for naics in NAICS_NAMES: sigs = sector_sigmas[naics] if sigs: result[naics] = { 'median': np.median(sigs), 'mean': np.mean(sigs), 'n': len(sigs), } return result def _get_bea_data(bea): """Parse BEA composition data into {(industry, year): {desc: value}}.""" data = {} for _, row in bea.iterrows(): ind = row['industry'] if isinstance(ind, float): ind = str(int(ind)) else: ind = str(ind).strip() yr = int(row['year']) desc = str(row['desc']).strip().strip('"') val = float(row['value']) key = (ind, yr) if key not in data: data[key] = {} data[key][desc] = val return data # ── Sub-tests ───────────────────────────────────────────────────── def test_sector_elasticities(sector_sigmas): """Test 1: Meaningful variation in σ across sectors.""" valid_medians = [v['median'] for v in sector_sigmas.values()] if len(valid_medians) < 5: return log_result('sectoral_ramsey_sector_elasticities', np.nan, np.nan, 'AMBIGUOUS', details={'n_sectors': len(valid_medians)}, description="Tests that trade elasticities sigma vary meaningfully across manufacturing sectors (CV > 10%). Aggregates Broda-Weinstein SITC3-level sigmas to NAICS sectors via concordance table.") cv = np.std(valid_medians) / np.mean(valid_medians) verd = 'CONSISTENT' if cv > 0.1 else 'AMBIGUOUS' return log_result('sectoral_ramsey_sector_elasticities', cv, np.nan, verd, details={'n_sectors': len(valid_medians), 'cv': cv}, description="Tests that trade elasticities sigma vary meaningfully across manufacturing sectors (CV > 10%). Aggregates Broda-Weinstein SITC3-level sigmas to NAICS sectors via concordance table.") def test_sector_tax_rates(bea_data): """Test 2: Tax rates vary meaningfully across sectors.""" sector_taxes = {} for naics in NAICS_NAMES: ratios = [] for yr in range(1997, 2025): key = (naics, yr) if key in bea_data: d = bea_data[key] taxes = d.get('Taxes on production and imports less subsidies', 0) va = d.get('Value added', 0) if va > 0: ratios.append(taxes / va) if ratios: sector_taxes[naics] = np.mean(ratios) if len(sector_taxes) < 5: return log_result('sectoral_ramsey_tax_rates', np.nan, np.nan, 'AMBIGUOUS', description="Tests that effective tax rates (taxes on production less subsidies / value added) vary meaningfully across manufacturing sectors (CV > 10%). Uses BEA composition data averaged over 1997-2024."), sector_taxes taus = list(sector_taxes.values()) cv = np.std(taus) / np.mean(taus) verd = 'CONSISTENT' if cv > 0.1 else 'AMBIGUOUS' return log_result('sectoral_ramsey_tax_rates', cv, np.nan, verd, details={'n_sectors': len(sector_taxes), 'cv': cv}, description="Tests that effective tax rates (taxes on production less subsidies / value added) vary meaningfully across manufacturing sectors (CV > 10%). Uses BEA composition data averaged over 1997-2024."), sector_taxes def test_ramsey_inverse_elasticity(sector_sigmas, sector_taxes): """Test 3: Spearman rho(sigma, tau) < 0 — higher elasticity -> lower tax. NORMATIVE test: the Ramsey rule predicts what taxes SHOULD be under optimal policy. US manufacturing taxes are set by political economy, not Ramsey optimization. An insignificant result means we cannot detect whether policy follows the normative rule — capped at AMBIGUOUS. """ desc = ("NORMATIVE: Tests the Ramsey inverse-elasticity rule (proved: " "saezTopRate_decreasing_in_e). Spearman rank correlation between " "sector median sigma (Broda-Weinstein SITC3→NAICS, N=9) and " "effective tax rate (BEA taxes on production less subsidies / VA, " "1997-2024 mean). Prediction: rho < 0 (higher demand elasticity → " "lower optimal tax). Capped at AMBIGUOUS because (1) normative " "prediction — policy need not follow Ramsey optimum; (2) sign is " "ambiguous — CES production-side curvature K predicts positive " "sigma-tau correlation (Paper 2 Remark rem:ramsey-ces).") paired_sigmas, paired_taus = [], [] for naics in NAICS_NAMES: if naics in sector_sigmas and naics in sector_taxes: paired_sigmas.append(sector_sigmas[naics]['median']) paired_taus.append(sector_taxes[naics]) if len(paired_sigmas) < 4: return log_result('sectoral_ramsey_inverse_elasticity', np.nan, np.nan, 'SKIPPED', description=desc) rho_s, p_scipy = _spearman(paired_sigmas, paired_taus) p_perm = _permutation_p(paired_sigmas, paired_taus, rho_s, n_perm=5000, seed=45, alternative='less') if rho_s < 0 and p_perm < 0.10: verd = 'CONSISTENT' elif rho_s < 0: # Right sign but not significant verd = 'AMBIGUOUS' else: # Wrong sign — but this is normative, so cap at AMBIGUOUS # (policy doesn't follow the rule, not evidence against theorem) verd = 'AMBIGUOUS' return log_result('sectoral_ramsey_inverse_elasticity', rho_s, p_perm, verd, details={'spearman_rho': round(rho_s, 4), 'p_scipy': round(p_scipy, 4), 'p_perm': round(p_perm, 4), 'n': len(paired_sigmas)}, description=desc) def test_recession_ordering(sector_sigmas, bea_data): """Test 4: Higher sigma -> smaller recession drops (rho > 0). Confounded by durability: motor vehicles and primary metals are both high-sigma AND highly cyclical due to durable-goods demand volatility, not substitutability. Only 3 recessions with very different causes (tech bubble, GFC, pandemic). Capped at AMBIGUOUS when insignificant given these design limitations. """ desc = ("Tests Proposition recession_ordering: higher-sigma sectors should " "experience smaller recession output drops. Spearman correlation " "between sector median sigma (Broda-Weinstein, N=9) and mean " "VA decline (BEA) across three recessions (2001, 2009, 2020). " "Prediction: rho > 0. Capped at AMBIGUOUS due to N=9, durability " "confound (motor vehicles/primary metals are high-sigma AND " "highly cyclical), and heterogeneous recession types.") # Value added by (industry, year) — desc == 'Value added' va_data = {} for (naics, yr), d in bea_data.items(): if 'Value added' in d: va_data[(naics, yr)] = d['Value added'] recessions = [(2000, 2001), (2007, 2009), (2019, 2020)] sector_drops = {naics: [] for naics in NAICS_NAMES} for peak_yr, trough_yr in recessions: for naics in NAICS_NAMES: va_peak = va_data.get((naics, peak_yr)) va_trough = va_data.get((naics, trough_yr)) if va_peak and va_trough and va_peak > 0: sector_drops[naics].append((va_trough - va_peak) / va_peak) paired_sigmas, paired_drops = [], [] for naics in NAICS_NAMES: if naics in sector_sigmas and sector_drops[naics]: paired_sigmas.append(sector_sigmas[naics]['median']) paired_drops.append(np.mean(sector_drops[naics])) if len(paired_sigmas) < 4: return log_result('sectoral_ramsey_recession_ordering', np.nan, np.nan, 'SKIPPED', description=desc) rho_s, p_scipy = _spearman(paired_sigmas, paired_drops) p_perm = _permutation_p(paired_sigmas, paired_drops, rho_s, n_perm=5000, seed=46, alternative='greater') if rho_s > 0 and p_perm < 0.10: verd = 'CONSISTENT' elif rho_s > 0: verd = 'AMBIGUOUS' else: # Wrong sign — cap at AMBIGUOUS due to design limitations # (durability confound, heterogeneous recessions, N=9) verd = 'AMBIGUOUS' return log_result('sectoral_ramsey_recession_ordering', rho_s, p_perm, verd, details={'spearman_rho': round(rho_s, 4), 'p_scipy': round(p_scipy, 4), 'p_perm': round(p_perm, 4), 'n': len(paired_sigmas)}, description=desc) def test_revenue_decomposition(sector_sigmas, sector_taxes): """Test 5: Low-sigma sectors have higher effective tax rates than high-sigma. Normative prediction, same caveats as Test 3. """ desc = ("NORMATIVE: Tests that low-sigma sectors bear higher effective tax " "rates than high-sigma sectors, as predicted by the Ramsey " "inverse-elasticity rule. Splits sectors at the median sigma and " "compares mean effective tax rates. Difference near zero indicates " "tax policy does not differentiate by substitutability.") paired = [] for naics in NAICS_NAMES: if naics in sector_sigmas and naics in sector_taxes: paired.append({ 'sigma': sector_sigmas[naics]['median'], 'tau': sector_taxes[naics], }) if len(paired) < 4: return log_result('sectoral_ramsey_revenue_decomposition', np.nan, np.nan, 'SKIPPED', description=desc) paired.sort(key=lambda x: x['sigma']) mid = len(paired) // 2 mean_tau_low = np.mean([p['tau'] for p in paired[:mid]]) mean_tau_high = np.mean([p['tau'] for p in paired[mid:]]) if mean_tau_low > mean_tau_high: verd = 'CONSISTENT' elif abs(mean_tau_low - mean_tau_high) < 0.01: verd = 'AMBIGUOUS' else: # Wrong direction but tiny difference — normative test, cap at AMBIGUOUS verd = 'AMBIGUOUS' return log_result('sectoral_ramsey_revenue_decomposition', mean_tau_low - mean_tau_high, np.nan, verd, details={'mean_tau_low_sigma': mean_tau_low, 'mean_tau_high_sigma': mean_tau_high}, description=desc) # ── Main ────────────────────────────────────────────────────────── _cache = {} def main(): """Run all sectoral Ramsey tests. Returns list of result dicts.""" try: bw = load_compiled('broda_weinstein_us_sigmas') bea = load_compiled('bea_composition') except FileNotFoundError as e: return [log_result('sectoral_ramsey', np.nan, np.nan, 'SKIPPED', details={'error': str(e)}, description="Master fallback: all sectoral Ramsey tests skipped due to data loading failure for Broda-Weinstein sigmas or BEA composition data.")] sector_sigmas = _aggregate_sigmas(bw) bea_data = _get_bea_data(bea) results = [] results.append(test_sector_elasticities(sector_sigmas)) tax_result, sector_taxes = test_sector_tax_rates(bea_data) results.append(tax_result) results.append(test_ramsey_inverse_elasticity(sector_sigmas, sector_taxes)) results.append(test_recession_ordering(sector_sigmas, bea_data)) results.append(test_revenue_decomposition(sector_sigmas, sector_taxes)) # Cache data for figures _cache['results'] = results _cache['sector_sigmas'] = sector_sigmas _cache['sector_taxes'] = sector_taxes _cache['bea_data'] = bea_data return results _ALL_SUBTESTS = [ 'sectoral_ramsey_sector_elasticities', 'sectoral_ramsey_tax_rates', 'sectoral_ramsey_inverse_elasticity', 'sectoral_ramsey_recession_ordering', 'sectoral_ramsey_revenue_decomposition', ] def figures(): """Generate per-subtest figures for sectoral Ramsey tests.""" if not _cache: main() paths = [] sector_sigmas = _cache.get('sector_sigmas', {}) sector_taxes = _cache.get('sector_taxes', {}) bea_data = _cache.get('bea_data', {}) if not sector_sigmas or not sector_taxes: return paths # --- Figure 1: Sector elasticities bar chart --- medians = {naics: sector_sigmas[naics]['median'] for naics in NAICS_NAMES if naics in sector_sigmas} if medians: fig, ax = setup_figure() names = [NAICS_NAMES[k] for k in medians] vals = list(medians.values()) bars = ax.barh(names, vals, color='steelblue', alpha=0.8) ax.set_xlabel('Median trade elasticity (sigma)') ax.set_title('Sector Elasticities (Broda-Weinstein)') cv = np.std(vals) / np.mean(vals) ax.axvline(np.mean(vals), color='red', ls='--', alpha=0.5, label=f'Mean={np.mean(vals):.1f}') ax.legend(fontsize=8) ax.text(0.98, 0.02, f'CV = {cv:.2f}', transform=ax.transAxes, ha='right', va='bottom', fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) fig.tight_layout() paths.extend(save_figure(fig, 'sectoral_ramsey_elasticities', subfolder='macro')) # --- Figure 2: Tax rates bar chart --- if sector_taxes: fig, ax = setup_figure() names = [NAICS_NAMES[k] for k in sector_taxes] vals = [v * 100 for v in sector_taxes.values()] ax.barh(names, vals, color='coral', alpha=0.8) ax.set_xlabel('Effective tax rate (%)') ax.set_title('Sector Tax Rates (BEA)') cv = np.std(vals) / np.mean(vals) ax.text(0.98, 0.02, f'CV = {cv:.2f}', transform=ax.transAxes, ha='right', va='bottom', fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) fig.tight_layout() paths.extend(save_figure(fig, 'sectoral_ramsey_tax_rates', subfolder='macro')) # --- Figure 3: Inverse elasticity scatter --- paired_sigmas, paired_taus, paired_names = [], [], [] for naics in NAICS_NAMES: if naics in sector_sigmas and naics in sector_taxes: paired_sigmas.append(sector_sigmas[naics]['median']) paired_taus.append(sector_taxes[naics] * 100) paired_names.append(NAICS_NAMES[naics]) if paired_sigmas: fig, ax = setup_figure() ax.scatter(paired_sigmas, paired_taus, c='steelblue', s=60, zorder=3) for name, s, t in zip(paired_names, paired_sigmas, paired_taus): ax.annotate(name, (s, t), fontsize=7, ha='left', va='bottom', xytext=(4, 4), textcoords='offset points') rho_s, p = _spearman(paired_sigmas, paired_taus) # Fit line z = np.polyfit(paired_sigmas, paired_taus, 1) xs = np.linspace(min(paired_sigmas), max(paired_sigmas), 50) ax.plot(xs, np.polyval(z, xs), 'r--', alpha=0.5) ax.set_xlabel('Median trade elasticity (sigma)') ax.set_ylabel('Effective tax rate (%)') ax.set_title('Ramsey Inverse-Elasticity Test') ax.text(0.98, 0.98, f'Spearman rho = {rho_s:.3f}\nPrediction: rho < 0\nN = {len(paired_sigmas)}', transform=ax.transAxes, ha='right', va='top', fontsize=9, bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8)) ax.grid(True, alpha=0.3) fig.tight_layout() paths.extend(save_figure(fig, 'sectoral_ramsey_inverse_elasticity', subfolder='macro')) # --- Figure 4: Recession ordering --- va_data = {} for (naics, yr), d in bea_data.items(): if 'Value added' in d: va_data[(naics, yr)] = d['Value added'] recessions = [(2000, 2001), (2007, 2009), (2019, 2020)] sector_drops = {naics: [] for naics in NAICS_NAMES} for peak_yr, trough_yr in recessions: for naics in NAICS_NAMES: va_peak = va_data.get((naics, peak_yr)) va_trough = va_data.get((naics, trough_yr)) if va_peak and va_trough and va_peak > 0: sector_drops[naics].append((va_trough - va_peak) / va_peak) paired_s, paired_d, paired_n = [], [], [] for naics in NAICS_NAMES: if naics in sector_sigmas and sector_drops[naics]: paired_s.append(sector_sigmas[naics]['median']) paired_d.append(np.mean(sector_drops[naics]) * 100) paired_n.append(NAICS_NAMES[naics]) if paired_s: fig, ax = setup_figure() ax.scatter(paired_s, paired_d, c='darkorange', s=60, zorder=3) for name, s, d in zip(paired_n, paired_s, paired_d): ax.annotate(name, (s, d), fontsize=7, ha='left', va='bottom', xytext=(4, 4), textcoords='offset points') rho_s, _ = _spearman(paired_s, paired_d) ax.set_xlabel('Median trade elasticity (sigma)') ax.set_ylabel('Mean recession VA drop (%)') ax.set_title('Recession Ordering vs Elasticity') ax.text(0.98, 0.02, f'Spearman rho = {rho_s:.3f}\nPrediction: rho > 0\nN = {len(paired_s)}', transform=ax.transAxes, ha='right', va='bottom', fontsize=9, bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8)) ax.grid(True, alpha=0.3) fig.tight_layout() paths.extend(save_figure(fig, 'sectoral_ramsey_recession', subfolder='macro')) # --- Figure 5: Revenue decomposition --- paired = [] for naics in NAICS_NAMES: if naics in sector_sigmas and naics in sector_taxes: paired.append({ 'name': NAICS_NAMES[naics], 'sigma': sector_sigmas[naics]['median'], 'tau': sector_taxes[naics] * 100, }) if len(paired) >= 4: paired.sort(key=lambda x: x['sigma']) mid = len(paired) // 2 fig, ax = setup_figure() colors = ['#d73027' if i < mid else '#1a9850' for i in range(len(paired))] ax.barh([p['name'] for p in paired], [p['tau'] for p in paired], color=colors, alpha=0.8) ax.set_xlabel('Effective tax rate (%)') ax.set_title('Revenue Decomposition: Low-sigma (red) vs High-sigma (green)') mean_lo = np.mean([p['tau'] for p in paired[:mid]]) mean_hi = np.mean([p['tau'] for p in paired[mid:]]) ax.text(0.98, 0.02, f'Mean low-sigma: {mean_lo:.2f}%\nMean high-sigma: {mean_hi:.2f}%\nDiff: {mean_lo - mean_hi:+.3f}pp', transform=ax.transAxes, ha='right', va='bottom', fontsize=9, bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8)) fig.tight_layout() paths.extend(save_figure(fig, 'sectoral_ramsey_revenue', subfolder='macro')) 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']})")

n	9
p perm	0.6934
p scipy	0.6368
p value	0.6934
statistic	-0.1833
spearman rho	-0.1833

Sectoral Ramsey Recession Ordering

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps