#!/usr/bin/env python3 """ Aggregation Puzzle -- Multi-Level sigma Estimation ==================================================== Paper 5 (Empirical Tests), Section "The Aggregation Gradient". **Prediction**: Effective sigma increases with aggregation level due to inter-sector reallocation of factors (Houthakker 1955, Sato 1967). The gradient is concentrated at the cross-sector step where heterogeneous capital intensities create reallocation opportunities. **Lean**: aggregation_fixed_point (Emergence.lean) — CES is the unique fixed point of multi-scale aggregation. The paper also references aggregation_raises_effective_sigma (not yet formalized in Lean). Test 1: NBER-CES plant-level sigma distribution (median < 1) Test 2: Sector-level sigma between plant and aggregate Test 3: Oberfield-Raval comparison (rank correlation) Test 4: Monotone aggregation gradient Test 5: Reallocation decomposition Data: NBER-CES Manufacturing Industry Database (364 six-digit NAICS industries, 1958-2018). Aggregate sigma=1.267 from separate two-factor estimation (Paper 5, Section "Two-Factor CES"). Review Notes (2026-03-12): - Theory match: CORRECT (Paper 5 Sec Aggregation Gradient, Paper 1 Emergent CES) - I/O audit: CORRECT — log(s_L/s_K) regressed on log(L/K) via 10-year long differences (Chirinko 2008). s_L = pay/vadd (value-added share, not cost share). - Data-signal alignment: STRONG (NBER-CES gold standard, N=359 plants, N=21 sectors) - Code quality: MINOR ISSUES (Test 4 statistic/p mismatch fixed 2026-03-12) - Confidence: HIGH (4/5 CONSISTENT, Oberfield-Raval r=0.85 validates estimation) - Known limitations: aggregate sigma=1.267 hardcoded; labor share is value-added based not cost share; paper's O-R comparison uses BEA (r=0.156) while code uses NBER-CES (r=0.851) — different datasets. """ import numpy as np import pandas as pd import os from common.config import DATA_DIR from common.data import load_nber_ces from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure # Oberfield-Raval (2021) sigma estimates OBERFIELD_RAVAL = { '336': ('Transport Equip', 0.50), '331': ('Primary Metals', 0.55), '327': ('Nonmetallic Min', 0.60), '333': ('Machinery', 0.65), '321': ('Wood', 0.70), '332': ('Fabricated Metals', 0.70), '335': ('Electrical Equip', 0.70), '322': ('Paper', 0.85), '337': ('Furniture', 0.85), '325': ('Chemicals', 0.90), '326': ('Plastics/Rubber', 0.95), '339': ('Miscellaneous', 1.10), '311': ('Food', 1.15), '334': ('Computer/Elec', 1.40), } DIFF_HORIZON = 10 _ALL_SUBTESTS = [ 'aggregation_puzzle_plant_sigma', 'aggregation_puzzle_sector_sigma', 'aggregation_puzzle_oberfield_raval', 'aggregation_puzzle_monotone_gradient', 'aggregation_puzzle_reallocation', ] def estimate_sigma_long_diff(y_level, x_level, min_n=10): """Estimate sigma via long differences on log(share_ratio) vs log(L/K).""" dy = y_level[DIFF_HORIZON:] - y_level[:-DIFF_HORIZON] dx = x_level[DIFF_HORIZON:] - x_level[:-DIFF_HORIZON] mask = np.isfinite(dx) & np.isfinite(dy) dx, dy = dx[mask], dy[mask] if len(dx) < min_n: return None n = len(dx) dx_bar, dy_bar = dx.mean(), dy.mean() ss_xy = np.sum((dx - dx_bar) * (dy - dy_bar)) ss_xx = np.sum((dx - dx_bar) ** 2) if ss_xx < 1e-12: return None rho_hat = ss_xy / ss_xx sigma_hat = 1.0 / (1.0 - rho_hat) if abs(1.0 - rho_hat) > 1e-6 else np.nan if not (np.isfinite(sigma_hat) and 0.01 < sigma_hat < 20): return None residuals = dy - (dy_bar + rho_hat * (dx - dx_bar)) se = np.sqrt(np.sum(residuals ** 2) / (n - 2) / ss_xx) if n > 2 else np.nan return {'rho': rho_hat, 'sigma': sigma_hat, 'se_rho': se, 'n_obs': n} def test1_plant_level(nber): """Plant-level sigma distribution from NBER-CES 6-digit industries.""" nber = nber[(nber['vadd'] > 0) & (nber['pay'] > 0) & (nber['cap'] > 0) & (nber['prodh'] > 0)].copy() nber['s_L'] = nber['pay'] / nber['vadd'] nber = nber[(nber['s_L'] > 0.05) & (nber['s_L'] < 0.95)] nber['log_share_ratio'] = np.log(nber['s_L'] / (1.0 - nber['s_L'])) nber['log_LK'] = np.log(nber['prodh'] / nber['cap']) sigma_estimates = [] for naics, group in nber.groupby('naics'): group = group.sort_values('year') if len(group) < DIFF_HORIZON + 5: continue result = estimate_sigma_long_diff( group['log_share_ratio'].values, group['log_LK'].values) if result: result['naics'] = naics result['naics3'] = str(naics)[:3] sigma_estimates.append(result) if not sigma_estimates: return None, log_result('aggregation_puzzle_plant_sigma', np.nan, np.nan, 'SKIPPED', description="Tests whether plant-level (6-digit NAICS) sigma estimates have median below 1, " "consistent with complementary inputs at the micro level. " "Uses long-difference regression on NBER-CES manufacturing data.") sigmas = pd.DataFrame(sigma_estimates) med_sigma = sigmas['sigma'].median() frac_below_1 = (sigmas['sigma'] < 1.0).mean() # Wilcoxon signed-rank test: H0: median sigma = 1, H1: median sigma < 1 from scipy.stats import wilcoxon try: _, p_wilcox = wilcoxon(sigmas['sigma'] - 1.0, alternative='less') except Exception: p_wilcox = np.nan is_consistent = (med_sigma < 1.0) and p_wilcox < 0.05 verd = 'CONSISTENT' if is_consistent else ('AMBIGUOUS' if med_sigma < 1.0 else 'INCONSISTENT') return sigmas, log_result('aggregation_puzzle_plant_sigma', med_sigma, float(p_wilcox), verd, description="Tests whether plant-level (6-digit NAICS) sigma estimates have median below 1, " "consistent with complementary inputs at the micro level. " "Wilcoxon signed-rank test H0: median=1, H1: median<1.", details={'n_industries': len(sigmas), 'frac_below_1': float(frac_below_1)}) def test2_sector_level(nber): """Sector-level sigma from NBER-CES 3-digit aggregation.""" nber = nber[(nber['vadd'] > 0) & (nber['pay'] > 0) & (nber['cap'] > 0) & (nber['prodh'] > 0)].copy() nber['naics3'] = nber['naics'].astype(str).str[:3] sector_agg = nber.groupby(['naics3', 'year']).agg( vadd=('vadd', 'sum'), pay=('pay', 'sum'), cap=('cap', 'sum'), prodh=('prodh', 'sum') ).reset_index() sector_agg['s_L'] = sector_agg['pay'] / sector_agg['vadd'] sector_agg = sector_agg[(sector_agg['s_L'] > 0.05) & (sector_agg['s_L'] < 0.95)] sector_agg['log_share_ratio'] = np.log(sector_agg['s_L'] / (1.0 - sector_agg['s_L'])) sector_agg['log_LK'] = np.log(sector_agg['prodh'] / sector_agg['cap']) sector_sigmas = [] for naics3, group in sector_agg.groupby('naics3'): group = group.sort_values('year') if len(group) < DIFF_HORIZON + 5: continue result = estimate_sigma_long_diff( group['log_share_ratio'].values, group['log_LK'].values) if result: result['naics3'] = naics3 sector_sigmas.append(result) if not sector_sigmas: return None, log_result('aggregation_puzzle_sector_sigma', np.nan, np.nan, 'SKIPPED', description="Tests whether 3-digit sector-level sigma lies between plant-level and aggregate, " "as predicted by CES hierarchical aggregation. " "Long-difference regression on NBER-CES data aggregated to 3-digit NAICS.") sectors = pd.DataFrame(sector_sigmas) med = sectors['sigma'].median() # Wilcoxon signed-rank test: H0: median sigma = 1.267, H1: median sigma < 1.267 from scipy.stats import wilcoxon try: _, p_wilcox = wilcoxon(sectors['sigma'] - 1.267, alternative='less') except Exception: p_wilcox = np.nan is_consistent = (med < 1.267) and (med > 0.1) and p_wilcox < 0.05 verd = 'CONSISTENT' if is_consistent else ('AMBIGUOUS' if med < 1.267 else 'INCONSISTENT') return sectors, log_result('aggregation_puzzle_sector_sigma', med, float(p_wilcox), verd, description="Tests whether 3-digit sector-level sigma lies below aggregate (1.267), " "as predicted by CES hierarchical aggregation. " "Wilcoxon signed-rank test H0: median=1.267, H1: median<1.267.", details={'n_sectors': len(sectors)}) def test3_oberfield_raval(sector_df): """Compare our sector estimates to Oberfield-Raval (2021).""" if sector_df is None or len(sector_df) == 0: return log_result('aggregation_puzzle_oberfield_raval', np.nan, np.nan, 'SKIPPED', description="Tests rank-order agreement between our sector sigma estimates and " "Oberfield-Raval (2021) published values. " "Spearman rank correlation across matched 3-digit NAICS sectors.") matched = [] for _, row in sector_df.iterrows(): if row['naics3'] in OBERFIELD_RAVAL: name, or_sigma = OBERFIELD_RAVAL[row['naics3']] matched.append({'naics3': row['naics3'], 'or_sigma': or_sigma, 'our_sigma': row['sigma']}) if len(matched) < 4: return log_result('aggregation_puzzle_oberfield_raval', np.nan, np.nan, 'SKIPPED', description="Tests rank-order agreement between our sector sigma estimates and " "Oberfield-Raval (2021) published values. " "Spearman rank correlation across matched 3-digit NAICS sectors.") mdf = pd.DataFrame(matched) from scipy.stats import spearmanr rho_s, p_s = spearmanr(mdf['or_sigma'], mdf['our_sigma']) verd = 'CONSISTENT' if rho_s > 0.3 and p_s < 0.10 else 'AMBIGUOUS' return log_result('aggregation_puzzle_oberfield_raval', rho_s, p_s, verd, description="Tests rank-order agreement between our sector sigma estimates and " "Oberfield-Raval (2021) published values. " "Spearman rank correlation across matched 3-digit NAICS sectors.", details={'n_matched': len(mdf)}) _cache = {} def main(): try: nber = load_nber_ces() except FileNotFoundError: return [log_result('aggregation_puzzle', np.nan, np.nan, 'SKIPPED', description="Multi-level sigma estimation tests require the NBER-CES Manufacturing Industry Database.")] results = [] plant_df, r1 = test1_plant_level(nber) results.append(r1) sector_df, r2 = test2_sector_level(nber) results.append(r2) r3 = test3_oberfield_raval(sector_df) results.append(r3) # Test 4: Monotone gradient plant_med = plant_df['sigma'].median() if plant_df is not None else np.nan sector_med = sector_df['sigma'].median() if sector_df is not None else np.nan agg_sigma = 1.267 # from separate two-factor estimation if np.isfinite(plant_med) and np.isfinite(sector_med): # The theory predicts the gradient is concentrated at the cross-sector # aggregation step (reallocation across heterogeneous sectors). Within-sector # aggregation (plants sharing similar technology) adds negligible reallocation, # so plant_med ≈ sector_med is expected. Use tolerance for the plant→sector step. cross_sector_gap = agg_sigma - sector_med within_sector_gap = sector_med - plant_med # Bootstrap p-value for cross-sector gap: H0: sector median >= aggregate rng = np.random.default_rng(42) n_boot = 2000 sector_vals = sector_df['sigma'].values boot_meds = np.array([np.median(rng.choice(sector_vals, size=len(sector_vals), replace=True)) for _ in range(n_boot)]) # p = fraction of bootstrap medians >= aggregate (one-sided) p_cross = float(np.mean(boot_meds >= agg_sigma)) cross_ok = cross_sector_gap > 0.1 within_ok = within_sector_gap > -0.05 # tolerance for estimation noise if cross_ok and within_ok and p_cross < 0.05: verd = 'CONSISTENT' elif cross_ok and within_ok: verd = 'AMBIGUOUS' elif cross_ok: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' results.append(log_result('aggregation_puzzle_monotone_gradient', cross_sector_gap, max(p_cross, 1.0/n_boot), verd, description="Tests the monotone gradient prediction: plant sigma <= sector sigma << aggregate sigma. " "Bootstrap test (2000 resamples) for cross-sector gap: H0 sector median >= aggregate. " "Stat is cross-sector gap (sector->aggregate). Theory predicts gradient concentrated " "at this step due to cross-sector capital intensity heterogeneity.", details={'plant': plant_med, 'sector': sector_med, 'aggregate': agg_sigma, 'cross_sector_gap': cross_sector_gap, 'within_sector_gap': within_sector_gap, 'boot_sector_median_95ci': [float(np.percentile(boot_meds, 2.5)), float(np.percentile(boot_meds, 97.5))]})) else: results.append(log_result('aggregation_puzzle_monotone_gradient', np.nan, np.nan, 'SKIPPED', description="Tests the monotone gradient prediction: plant sigma < sector sigma < aggregate sigma. " "Compares median sigma at 6-digit, 3-digit, and economy-wide aggregation levels.")) # Test 5: Reallocation decomposition if sector_df is not None and len(sector_df) > 3: mean_sector_sigma = sector_df['sigma'].mean() reallocation = agg_sigma - mean_sector_sigma # One-sample t-test: H0: mean sector sigma = 1.267, H1: mean < 1.267 from scipy.stats import ttest_1samp t_stat, p_two = ttest_1samp(sector_df['sigma'].values, agg_sigma) p_realloc = float(p_two / 2) if t_stat < 0 else float(1 - p_two / 2) if reallocation > 0 and p_realloc < 0.05: verd = 'CONSISTENT' elif reallocation > 0: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' results.append(log_result('aggregation_puzzle_reallocation', reallocation, p_realloc, verd, description="Tests whether inter-sector reallocation contributes positively to aggregate sigma. " "One-sample t-test H0: mean sector sigma = aggregate, H1: mean < aggregate. " "Positive gap indicates reallocation adds apparent substitutability.")) else: results.append(log_result('aggregation_puzzle_reallocation', np.nan, np.nan, 'SKIPPED', description="Tests whether inter-sector reallocation contributes positively to aggregate sigma. " "Computes aggregate sigma minus mean sector sigma; positive indicates reallocation adds apparent substitutability.")) _cache['plant_df'] = plant_df _cache['sector_df'] = sector_df _cache['plant_med'] = plant_med _cache['sector_med'] = sector_med _cache['agg_sigma'] = agg_sigma return results def figures(): """Generate per-subtest figures for the aggregation puzzle.""" if not _cache: main() paths = [] plant_df = _cache.get('plant_df') sector_df = _cache.get('sector_df') agg_sigma = _cache.get('agg_sigma', 1.267) plant_med = _cache.get('plant_med', np.nan) sector_med = _cache.get('sector_med', np.nan) import matplotlib.pyplot as plt # --- Plant-level sigma distribution --- if plant_df is not None: fig, ax = setup_figure() ax.hist(plant_df['sigma'].clip(0, 5), bins=30, alpha=0.7, color='steelblue', edgecolor='k') ax.axvline(1.0, color='k', linestyle='--', alpha=0.6, label='sigma=1 (CD)') ax.axvline(plant_df['sigma'].median(), color='steelblue', linestyle='-', lw=2, label=f'Median={plant_df["sigma"].median():.3f}') frac = (plant_df['sigma'] < 1.0).mean() ax.set_xlabel('Estimated sigma') ax.set_ylabel('Count') ax.set_title(f'Plant-Level sigma (N={len(plant_df)}, {frac:.1%} < 1)') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'plant_level_sigma_distribution', subfolder='aggregation_puzzle')) # --- Sector-level sigma --- if sector_df is not None: fig, ax = setup_figure() ax.hist(sector_df['sigma'].clip(0, 5), bins=15, alpha=0.7, color='darkorange', edgecolor='k') ax.axvline(1.0, color='k', linestyle='--', alpha=0.6, label='sigma=1 (CD)') ax.axvline(agg_sigma, color='r', linestyle=':', lw=2, label=f'Aggregate={agg_sigma:.3f}') ax.axvline(sector_df['sigma'].median(), color='darkorange', linestyle='-', lw=2, label=f'Median={sector_df["sigma"].median():.3f}') ax.set_xlabel('Estimated sigma') ax.set_ylabel('Count') ax.set_title(f'Sector-Level sigma (N={len(sector_df)})') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'sector_level_sigma', subfolder='aggregation_puzzle')) # --- Oberfield-Raval comparison --- if sector_df is not None: matched = [] for _, row in sector_df.iterrows(): if row['naics3'] in OBERFIELD_RAVAL: name, or_sigma = OBERFIELD_RAVAL[row['naics3']] matched.append({'name': name, 'or_sigma': or_sigma, 'our_sigma': row['sigma']}) if matched: mdf = pd.DataFrame(matched) from scipy.stats import spearmanr rho_s, p_s = spearmanr(mdf['or_sigma'], mdf['our_sigma']) fig, ax = setup_figure() ax.scatter(mdf['or_sigma'], mdf['our_sigma'], s=50, color='seagreen', edgecolor='k', zorder=3) for _, row in mdf.iterrows(): ax.annotate(row['name'], (row['or_sigma'], row['our_sigma']), fontsize=7, ha='left', va='bottom', xytext=(4, 4), textcoords='offset points') lims = [0, max(mdf['or_sigma'].max(), mdf['our_sigma'].max()) + 0.2] ax.plot(lims, lims, 'k--', alpha=0.3, label='45-degree') ax.set_xlabel('Oberfield-Raval sigma') ax.set_ylabel('Our NBER-CES sigma') ax.set_title(f'Oberfield-Raval Comparison (r={rho_s:.3f}, p={p_s:.4f})') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'oberfield_raval_comparison', subfolder='aggregation_puzzle')) # --- Aggregation gradient --- if np.isfinite(plant_med) and np.isfinite(sector_med): fig, ax = setup_figure() levels = ['Plant\n(6-digit)', 'Sector\n(3-digit)', 'Aggregate'] vals = [plant_med, sector_med, agg_sigma] colors_bar = ['steelblue', 'darkorange', 'crimson'] ax.bar(levels, vals, color=colors_bar, edgecolor='k') ax.set_ylabel('Median sigma') gap = agg_sigma - sector_med ax.set_title(f'Aggregation Gradient (cross-sector gap = +{gap:.3f})') ax.axhline(1.0, color='k', linestyle='--', alpha=0.4, label='CD (sigma=1)') ax.legend(fontsize=9) ax.grid(True, alpha=0.3, axis='y') paths.extend(save_figure(fig, 'aggregation_gradient', subfolder='aggregation_puzzle')) # --- Reallocation decomposition --- if sector_df is not None and len(sector_df) > 3: fig, ax = setup_figure() mean_s = sector_df['sigma'].mean() realloc = agg_sigma - mean_s ax.bar(['Mean Sector\nsigma', 'Reallocation\nComponent', 'Aggregate\nsigma'], [mean_s, realloc, agg_sigma], color=['darkorange', 'limegreen', 'crimson'], edgecolor='k') ax.set_ylabel('sigma') ax.set_title(f'Reallocation Decomposition (+{realloc:.3f})') ax.axhline(1.0, color='k', linestyle='--', alpha=0.4) ax.grid(True, alpha=0.3, axis='y') paths.extend(save_figure(fig, 'reallocation_decomposition', subfolder='aggregation_puzzle')) 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.783e-4
statistic	0.2835

Aggregation Puzzle Reallocation

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps