#!/usr/bin/env python3 """ Test: Wage Complementarity and CES Factor Shares =================================================== Paper 5 (Empirical Tests), Section "Wage-Complementarity Link". **Prediction**: CES curvature controls adjustment speed — low-sigma sectors (strong complements) should have stickier wages and more stable factor shares than high-sigma sectors. **Lean**: factorShare_eq_escort (Economics.lean) — factor shares are CES escort distributions; laborShare_increases_with_kl_ratio_complements (TwoFactorCES.lean) — labor share monotonicity. Five sub-tests: 1. Compensation volatility positively correlated with sigma 2. Wage stickiness gradient (low-sigma sectors have stickier wages) 3. Wage cyclicality (high-sigma → comp growth more correlated with VA growth) 4. Factor share stability (low-sigma → lower CV of labor share s_L = comp/VA) 5. Combined assessment Data: bea_composition.csv (9 NAICS mfg sectors, 1997-2024), broda_weinstein_us_sigmas.csv (2716 product sigmas via SITC2→NAICS). Review Notes (2026-03-12): - Theory match: CORRECT (Paper 5, Sec Wage-Complementarity Link) - I/O audit: CORRECT after fixes (Tests 3 & 4 had bugs — see Change Log) - Data-signal alignment: WEAK (N=9, all sigma>1 substitute regime) - Code quality: BUGS FOUND and FIXED (Tests 1&3 identical, Test 4 wrong variable) - Confidence: LOW (N=9 underpowered; all sectors in substitute regime) - Known limitations: Broda-Weinstein sigmas are trade elasticities, not production elasticities. No sectors with sigma<1 (complement regime) where CES stabilization is strongest. """ import numpy as np import pandas as pd from scipy import 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 verdict, log_result from common.plotting import setup_figure, setup_subplots, save_figure # ── SITC2 → NAICS concordance (manufacturing only) ────────────── SITC2_TO_NAICS = { "61": "331", "62": "331", "63": "337", "64": "339", "65": "339", "66": "331", "67": "331", "68": "332", "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", } _ALL_SUBTESTS = [ 'wage_complementarity_volatility', 'wage_complementarity_stickiness', 'wage_complementarity_cyclicality', 'wage_complementarity_factor_share_stability', 'wage_complementarity_combined', ] def load_data(): """Load BEA composition and Broda-Weinstein elasticities.""" bea = load_compiled('bea_composition') bw = load_compiled('broda_weinstein_us_sigmas') return bea, bw def aggregate_sigmas(bw): """Aggregate Broda-Weinstein product-level sigmas to NAICS sectors via SITC2.""" sector_sigmas = {naics: [] for naics in NAICS_NAMES} for _, row in bw.iterrows(): sitc_code = str(int(row['sitc3_code'])) if isinstance(row['sitc3_code'], float) else str(row['sitc3_code']).strip() sigma = row['sigma'] if sigma <= 0 or sigma > 100: continue prefix = sitc_code[:2] if prefix in SITC2_TO_NAICS: sector_sigmas[SITC2_TO_NAICS[prefix]].append(sigma) result = {} for naics in NAICS_NAMES: sigs = sector_sigmas[naics] if len(sigs) > 0: result[naics] = np.median(sigs) return result def get_sector_time_series(bea, desc_filter='Compensation'): """Extract sector-level time series from BEA composition data (long format).""" if 'desc' not in bea.columns or 'industry' not in bea.columns: return None comp = bea[bea['desc'].str.contains(desc_filter, case=False, na=False)] if comp.empty: return None sectors = {} for ind, grp in comp.groupby('industry'): s = grp.set_index('year')['value'].sort_index().dropna() if len(s) >= 3: sectors[str(ind)] = s return sectors if sectors else None # ── Sub-tests ───────────────────────────────────────────────────── def test_volatility_vs_sigma(bea, bw): """Test 1: Compensation volatility negatively correlated with sigma.""" if bea is None or bw is None: return log_result('wage_complementarity_volatility', 0.0, 1.0, 'SKIPPED', description="Tests whether compensation growth volatility is correlated with sector-level sigma. Uses Spearman rank correlation between Broda-Weinstein median sigma and standard deviation of log-differenced BEA compensation by NAICS sector.") sigma_dict = aggregate_sigmas(bw) sectors = get_sector_time_series(bea) if sectors is None or len(sigma_dict) < 3: return log_result('wage_complementarity_volatility', 0.0, 1.0, 'SKIPPED', description="Tests whether compensation growth volatility is correlated with sector-level sigma. Uses Spearman rank correlation between Broda-Weinstein median sigma and standard deviation of log-differenced BEA compensation by NAICS sector.") pairs = [] for naics, sigma in sigma_dict.items(): if naics in sectors: vals = sectors[naics] growth = np.diff(np.log(np.maximum(vals.values, 1e-10))) if len(growth) >= 3: vol = np.std(growth) pairs.append((sigma, vol, naics)) if len(pairs) < 4: return log_result('wage_complementarity_volatility', 0.0, 1.0, 'SKIPPED', description="Tests that high-sigma (more substitutable) sectors have more volatile compensation growth. Spearman rank correlation of Broda-Weinstein median sigma vs std(log-diff compensation) across 9 NAICS manufacturing sectors (BEA, 1997-2024). Prediction: r > 0.") sigmas = [p[0] for p in pairs] vols = [p[1] for p in pairs] r, p_val = stats.spearmanr(sigmas, vols) sign_correct = r > 0 # higher sigma → more volatile verd = verdict(p_val, sign_correct) _cache['volatility'] = {'sigmas': sigmas, 'vols': vols, 'names': [NAICS_NAMES.get(p[2], p[2]) for p in pairs], 'r': r, 'p': p_val} return log_result( 'wage_complementarity_volatility', r, p_val, verd, {'spearman_r': r, 'n_sectors': len(pairs)}, description="Tests that high-sigma (more substitutable) sectors have more volatile compensation growth. Spearman rank correlation of Broda-Weinstein median sigma vs std(log-diff compensation) across 9 NAICS manufacturing sectors (BEA, 1997-2024). Prediction: r > 0." ) def test_wage_stickiness(bea, bw): """Test 2: Wage stickiness gradient — low-sigma more sticky.""" if bea is None or bw is None: return log_result('wage_complementarity_stickiness', 0.0, 1.0, 'SKIPPED', description="Tests that low-sigma (more complementary) sectors have stickier wages, measured by higher AR(1) persistence in compensation growth. Uses Spearman rank correlation between sigma and AR(1) coefficient; prediction is r < 0.") sigma_dict = aggregate_sigmas(bw) sectors = get_sector_time_series(bea) if sectors is None or len(sigma_dict) < 3: return log_result('wage_complementarity_stickiness', 0.0, 1.0, 'SKIPPED', description="Tests that low-sigma (more complementary) sectors have stickier wages, measured by higher AR(1) persistence in compensation growth. Uses Spearman rank correlation between sigma and AR(1) coefficient; prediction is r < 0.") pairs = [] for naics, sigma in sigma_dict.items(): if naics in sectors: vals = sectors[naics] growth = np.diff(np.log(np.maximum(vals.values, 1e-10))) if len(growth) >= 5: ar1 = np.corrcoef(growth[:-1], growth[1:])[0, 1] if not np.isnan(ar1): pairs.append((sigma, ar1, naics)) if len(pairs) < 4: return log_result('wage_complementarity_stickiness', 0.0, 1.0, 'SKIPPED', description="Tests that low-sigma (more complementary) sectors have stickier wages, measured by higher AR(1) persistence in compensation growth. Spearman rank correlation of Broda-Weinstein sigma vs AR(1) of log-diff compensation across 9 NAICS sectors (BEA, 1997-2024). Prediction: r < 0.") sigmas = [p[0] for p in pairs] ar1s = [p[1] for p in pairs] r, p_val = stats.spearmanr(sigmas, ar1s) sign_correct = r < 0 # higher sigma → less sticky → lower AR(1) verd = verdict(p_val, sign_correct) _cache['stickiness'] = {'sigmas': sigmas, 'ar1s': ar1s, 'names': [NAICS_NAMES.get(p[2], p[2]) for p in pairs], 'r': r, 'p': p_val} return log_result( 'wage_complementarity_stickiness', r, p_val, verd, {'spearman_r': r, 'n_sectors': len(pairs)}, description="Tests that low-sigma (more complementary) sectors have stickier wages, measured by higher AR(1) persistence in compensation growth. Spearman rank correlation of Broda-Weinstein sigma vs AR(1) of log-diff compensation across 9 NAICS sectors (BEA, 1997-2024). Prediction: r < 0." ) def test_wage_cyclicality(bea, bw): """Test 3: High-sigma sectors have comp growth more correlated with VA growth.""" desc = ("Tests that high-sigma sectors have compensation growth more " "correlated with VA growth (more cyclical wages). Computes " "Pearson corr(log-diff comp, log-diff VA) per sector, then " "Spearman rank correlation of that cyclicality measure vs " "Broda-Weinstein sigma across 9 NAICS sectors (BEA, 1997-2024). " "Prediction: r > 0.") if bea is None or bw is None: return log_result('wage_complementarity_cyclicality', 0.0, 1.0, 'SKIPPED', description=desc) sigma_dict = aggregate_sigmas(bw) comp_sectors = get_sector_time_series(bea, 'Compensation') va_sectors = get_sector_time_series(bea, 'Value added') if comp_sectors is None or va_sectors is None: return log_result('wage_complementarity_cyclicality', 0.0, 1.0, 'SKIPPED', description=desc) pairs = [] for naics, sigma in sigma_dict.items(): if naics in comp_sectors and naics in va_sectors: comp = comp_sectors[naics] va = va_sectors[naics] # Align years common_years = comp.index.intersection(va.index) if len(common_years) < 5: continue comp_aligned = comp[common_years].sort_index() va_aligned = va[common_years].sort_index() comp_g = np.diff(np.log(np.maximum(comp_aligned.values, 1e-10))) va_g = np.diff(np.log(np.maximum(va_aligned.values, 1e-10))) if len(comp_g) >= 4: cyclicality = np.corrcoef(comp_g, va_g)[0, 1] if not np.isnan(cyclicality): pairs.append((sigma, cyclicality, naics)) if len(pairs) < 4: return log_result('wage_complementarity_cyclicality', 0.0, 1.0, 'SKIPPED', description=desc) sigmas = [p[0] for p in pairs] cycl = [p[1] for p in pairs] r, p_val = stats.spearmanr(sigmas, cycl) sign_correct = r > 0 # higher sigma → comp more correlated with VA verd = verdict(p_val, sign_correct) _cache['cyclicality'] = {'sigmas': sigmas, 'cycl': cycl, 'names': [NAICS_NAMES.get(p[2], p[2]) for p in pairs], 'r': r, 'p': p_val} return log_result( 'wage_complementarity_cyclicality', r, p_val, verd, {'spearman_r': r, 'n_sectors': len(pairs)}, description=desc ) def test_factor_share_stability(bea, bw): """Test 4: High-sigma sectors have less stable labor shares (s_L = comp/VA).""" desc = ("Tests that high-sigma sectors have less stable labor shares " "s_L = comp/VA over time. Computes CV of annual labor share per " "sector, then Spearman rank correlation vs Broda-Weinstein sigma " "across 9 NAICS sectors (BEA, 1997-2024). Prediction: r > 0 " "(higher sigma → more volatile factor shares).") if bea is None or bw is None: return log_result('wage_complementarity_factor_share_stability', 0.0, 1.0, 'SKIPPED', description=desc) sigma_dict = aggregate_sigmas(bw) comp_sectors = get_sector_time_series(bea, 'Compensation') va_sectors = get_sector_time_series(bea, 'Value added') if comp_sectors is None or va_sectors is None: return log_result('wage_complementarity_factor_share_stability', 0.0, 1.0, 'SKIPPED', description=desc) pairs = [] for naics, sigma in sigma_dict.items(): if naics in comp_sectors and naics in va_sectors: comp = comp_sectors[naics] va = va_sectors[naics] common_years = comp.index.intersection(va.index) if len(common_years) < 3: continue comp_aligned = comp[common_years].sort_index() va_aligned = va[common_years].sort_index() # Labor share = compensation / value added labor_share = comp_aligned / va_aligned.replace(0, np.nan) labor_share = labor_share.dropna() if len(labor_share) >= 3 and labor_share.mean() > 0: cv = labor_share.std() / labor_share.mean() pairs.append((sigma, cv, naics)) if len(pairs) < 4: return log_result('wage_complementarity_factor_share_stability', 0.0, 1.0, 'SKIPPED', description=desc) sigmas = [p[0] for p in pairs] cvs = [p[1] for p in pairs] r, p_val = stats.spearmanr(sigmas, cvs) sign_correct = r > 0 # higher sigma → less stable factor shares verd = verdict(p_val, sign_correct) _cache['factor_share'] = {'sigmas': sigmas, 'cvs': cvs, 'names': [NAICS_NAMES.get(p[2], p[2]) for p in pairs], 'r': r, 'p': p_val} return log_result( 'wage_complementarity_factor_share_stability', r, p_val, verd, {'spearman_r': r, 'n_sectors': len(pairs)}, description=desc ) def test_combined_assessment(results_so_far): """Test 5: Combined assessment of wage-complementarity predictions.""" consistent_count = sum( 1 for r in results_so_far if r['verdict'] == 'CONSISTENT') total = len(results_so_far) skipped = sum(1 for r in results_so_far if r['verdict'] == 'SKIPPED') active = total - skipped if active == 0: return log_result('wage_complementarity_combined', 0.0, 1.0, 'SKIPPED', description="Combined assessment of all four wage-complementarity subtests. Uses a binomial test on the fraction of CONSISTENT results (under H0: 25% consistent by chance). Consistent if majority of active subtests pass.") frac = consistent_count / active p_val = 1.0 - stats.binom.cdf(consistent_count - 1, active, 0.25) sign_correct = frac > 0.5 verd = verdict(p_val, sign_correct) return log_result( 'wage_complementarity_combined', frac, p_val, verd, {'consistent': consistent_count, 'active': active, 'skipped': skipped}, description="Combined assessment of all four wage-complementarity subtests (volatility, stickiness, cyclicality, factor share stability). Uses a binomial test on the fraction of CONSISTENT results under H0: 25% consistent by chance." ) # ── Main ────────────────────────────────────────────────────────── _cache = {} def main(): """Run all wage complementarity tests. Returns list of result dicts.""" bea, bw = None, None try: bea, bw = load_data() except FileNotFoundError as e: print(f" Data loading failed: {e}") results = [] results.append(test_volatility_vs_sigma(bea, bw)) results.append(test_wage_stickiness(bea, bw)) results.append(test_wage_cyclicality(bea, bw)) results.append(test_factor_share_stability(bea, bw)) results.append(test_combined_assessment(results)) _cache['results'] = results return results def figures(): """Generate per-subtest figures.""" if not _cache: main() paths = [] # --- Volatility scatter --- data = _cache.get('volatility') if data: fig, ax = setup_figure() ax.scatter(data['sigmas'], data['vols'], color='steelblue', s=60, zorder=3) for i, name in enumerate(data['names']): ax.annotate(name, (data['sigmas'][i], data['vols'][i]), fontsize=7, ha='left', va='bottom', xytext=(4, 4), textcoords='offset points') ax.set_xlabel('Broda-Weinstein σ (median)') ax.set_ylabel('Std(log Δ compensation)') ax.set_title(f"Compensation Volatility vs σ\nr={data['r']:.3f}, p={data['p']:.3f}") paths.extend(save_figure(fig, 'wage_complementarity_volatility', subfolder='macro')) # --- Stickiness scatter --- data = _cache.get('stickiness') if data: fig, ax = setup_figure() ax.scatter(data['sigmas'], data['ar1s'], color='darkorange', s=60, zorder=3) for i, name in enumerate(data['names']): ax.annotate(name, (data['sigmas'][i], data['ar1s'][i]), fontsize=7, ha='left', va='bottom', xytext=(4, 4), textcoords='offset points') ax.set_xlabel('Broda-Weinstein σ (median)') ax.set_ylabel('AR(1) of compensation growth') ax.set_title(f"Wage Stickiness vs σ\nr={data['r']:.3f}, p={data['p']:.3f}") ax.axhline(0, color='gray', ls='--', lw=0.8) paths.extend(save_figure(fig, 'wage_complementarity_stickiness', subfolder='macro')) # --- Cyclicality scatter --- data = _cache.get('cyclicality') if data: fig, ax = setup_figure() ax.scatter(data['sigmas'], data['cycl'], color='seagreen', s=60, zorder=3) for i, name in enumerate(data['names']): ax.annotate(name, (data['sigmas'][i], data['cycl'][i]), fontsize=7, ha='left', va='bottom', xytext=(4, 4), textcoords='offset points') ax.set_xlabel('Broda-Weinstein σ (median)') ax.set_ylabel('corr(Δlog comp, Δlog VA)') ax.set_title(f"Wage Cyclicality vs σ\nr={data['r']:.3f}, p={data['p']:.3f}") ax.axhline(0, color='gray', ls='--', lw=0.8) paths.extend(save_figure(fig, 'wage_complementarity_cyclicality', subfolder='macro')) # --- Factor share stability scatter --- data = _cache.get('factor_share') if data: fig, ax = setup_figure() ax.scatter(data['sigmas'], data['cvs'], color='crimson', s=60, zorder=3) for i, name in enumerate(data['names']): ax.annotate(name, (data['sigmas'][i], data['cvs'][i]), fontsize=7, ha='left', va='bottom', xytext=(4, 4), textcoords='offset points') ax.set_xlabel('Broda-Weinstein σ (median)') ax.set_ylabel('CV(labor share = comp/VA)') ax.set_title(f"Factor Share Stability vs σ\nr={data['r']:.3f}, p={data['p']:.3f}") paths.extend(save_figure(fig, 'wage_complementarity_factor_share', subfolder='macro')) # --- Combined summary bar --- results = _cache.get('results') if results: fig, ax = setup_figure() names = [r['name'].replace('wage_complementarity_', '') for r in results if r['name'] != 'wage_complementarity_combined'] verdicts = [r['verdict'] for r in results if r['name'] != 'wage_complementarity_combined'] colors = {'CONSISTENT': 'green', 'AMBIGUOUS': 'gold', 'INCONSISTENT': 'red', 'SKIPPED': 'gray'} bars = ax.barh(names, [r['statistic'] for r in results if r['name'] != 'wage_complementarity_combined'], color=[colors.get(v, 'gray') for v in verdicts]) ax.set_xlabel('Spearman r') ax.set_title('Wage Complementarity: Sub-test Statistics') ax.axvline(0, color='black', lw=0.8) for i, v in enumerate(verdicts): ax.annotate(v, (0, i), fontsize=8, ha='left', va='center', xytext=(5, 0), textcoords='offset points') paths.extend(save_figure(fig, 'wage_complementarity_combined', subfolder='macro')) return paths if __name__ == '__main__': results = main() for r in results: print(f" {r['name']}: {r['verdict']} " f"(stat={r['statistic']:.4f}, p={r['p_value']:.4f})")

p value	0.3085
n sectors	9
statistic	0.3833
spearman r	0.3833

Wage Complementarity Factor Share Stability

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources