#!/usr/bin/env python3 """ Within-Industry Productivity Dispersion ~ K -- WP1 Prediction 1.7 Prediction: Industries with higher complementarity (lower rho, higher K) show wider productivity dispersion. Cross-sector regression of log(90/10 productivity ratio) on log(K) should yield positive slope. Theory (Paper 2, Proposition: Productivity Dispersion): Var[log F] ~ K^4 / (4(J-1)^2 c^4) * Var[T_i] => sigma(log F) ~ K^2 => log(90/10 ratio) ~ K^2 (not K^4!) The 90/10 ratio relates to the std dev (K^2), not the variance (K^4). A log-log regression log(ratio) ~ beta*log(K) is an approximation since the true relationship is polynomial, not power-law. Lean: dispersion_increases_with_K (real proof, K^2 scaling) Data needed: NBER-CES manufacturing panel Method: OLS of log(90/10 ratio) on log(|K|) across 4-digit NAICS groups Status: IMPLEMENTED Review Notes (2026-03-12): Theory match: PARTIAL — Theory predicts Var(log F) ~ K^4 but the test uses log(90/10 ratio) which scales as sigma ~ K^2. The log-log regression is a functional form approximation (theory gives polynomial K^2, not power-law K^beta). The description "beta ~ 4" was incorrect; beta ~ 0.89 is consistent with the K^2 polynomial through the log-log lens over the empirical K range. I/O audit: CORRECT — vadd/emp = labor productivity (output per worker), pay/vadd = labor share, both standard. K computed from crude sigma estimate via 1/(1-s_L) which is approximate. Data-signal alignment: ADEQUATE — NBER-CES is DIRECT for manufacturing productivity. sigma = 1/(1-s_L) is a WEAK proxy for sigma (not the standard CES-consistent estimator). N=62 groups adequate. Code quality: MINOR ISSUES — Two-tailed p when one-tailed appropriate; panel OLS lacks clustered SEs (t=36 is inflated); |K| loses sign info. Confidence in verdict: MEDIUM — All 3 subtests CONSISTENT. Qualitative prediction (more K => more dispersion) confirmed. Quantitative K^4 claim was incorrectly stated (should be K^2 for sigma, K^4 for variance). """ import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) import numpy as np import pandas as pd from common.data import load_nber_ces from common.stats import log_result, verdict from common.ces import curvature_K, rho_from_sigma from common.plotting import setup_figure, setup_subplots, save_figure _cache = {} def main() -> list[dict]: try: df = load_nber_ces() except FileNotFoundError: return [ log_result( name='productivity_dispersion', stat=float('nan'), p=float('nan'), verd='SKIPPED', description='NBER-CES data not found.', ) ] needed = ['naics', 'year', 'vadd', 'emp', 'pay'] if not all(c in df.columns for c in needed): return [ log_result( name='productivity_dispersion', stat=float('nan'), p=float('nan'), verd='SKIPPED', description=f'Missing columns. Have: {list(df.columns[:20])}', ) ] df = df.dropna(subset=needed) df = df[(df['emp'] > 0) & (df['vadd'] > 0) & (df['pay'] > 0)].copy() df['naics'] = df['naics'].astype(str).str.strip() df['naics4'] = df['naics'].str[:4] # Labor productivity df['lp'] = df['vadd'] / df['emp'] df = df[df['lp'] > 0].copy() # Labor share for sigma estimation df['s_L'] = df['pay'] / df['vadd'] # Clip to sensible range df['s_L'] = df['s_L'].clip(0.05, 0.95) results = [] # For each 4-digit group and year, compute: # - 90/10 ratio of labor productivity # - J = number of 6-digit industries # - sigma_hat = 1/(1 - s_L) using median labor share # - K from sigma_hat and J records = [] for (naics4, year), gdf in df.groupby(['naics4', 'year']): J = gdf['naics'].nunique() if J < 3: continue lp_vals = gdf['lp'].values p90 = np.percentile(lp_vals, 90) p10 = np.percentile(lp_vals, 10) if p10 <= 0: continue ratio_90_10 = p90 / p10 # Estimate sigma from median labor share med_sL = gdf['s_L'].median() # sigma = 1/(1 - s_L) is a crude proxy # For s_L close to 1, sigma -> inf (Cobb-Douglas); for s_L small, sigma ~ 1 sigma_hat = 1.0 / (1.0 - med_sL) if sigma_hat <= 0 or not np.isfinite(sigma_hat): continue rho_hat = rho_from_sigma(sigma_hat) K = curvature_K(rho_hat, J) records.append({ 'naics4': naics4, 'year': year, 'ratio_90_10': ratio_90_10, 'J': J, 'sigma_hat': sigma_hat, 'rho_hat': rho_hat, 'K': K, }) if len(records) < 10: return [ log_result( name='productivity_dispersion', stat=float('nan'), p=float('nan'), verd='SKIPPED', description=f'Only {len(records)} group-year obs with J>=3.', ) ] rdf = pd.DataFrame(records) # Average across years for each 4-digit group for stability avg = rdf.groupby('naics4').agg( ratio_90_10=('ratio_90_10', 'median'), K=('K', 'median'), J=('J', 'median'), n_years=('year', 'count'), ).reset_index() # Use |K| since some groups may have K < 0 (substitutes regime) avg['absK'] = avg['K'].abs() avg = avg[avg['absK'] > 1e-6].copy() avg = avg[avg['ratio_90_10'] > 1.0].copy() # ratio must exceed 1 if len(avg) < 10: return [ log_result( name='productivity_dispersion', stat=float('nan'), p=float('nan'), verd='SKIPPED', description=f'Only {len(avg)} groups after filtering.', ) ] from scipy.stats import t as t_dist # OLS: log(90/10 ratio) ~ beta * log(|K|) # Theory predicts sigma(logF) ~ K^2, so log(ratio) ~ c*K^2 (polynomial) # Log-log regression is an approximation; we test beta > 0 y = np.log(avg['ratio_90_10'].values) x = np.log(avg['absK'].values) # Add constant for OLS X = np.column_stack([np.ones(len(x)), x]) try: beta, residuals, rank, sv = np.linalg.lstsq(X, y, rcond=None) except np.linalg.LinAlgError: return [ log_result( name='productivity_dispersion', stat=float('nan'), p=float('nan'), verd='SKIPPED', description='OLS failed (singular matrix).', ) ] slope = beta[1] n = len(y) k = 2 # intercept + slope # Compute standard error and p-value y_hat = X @ beta resid = y - y_hat rss = np.sum(resid ** 2) se_sq = rss / (n - k) XtX_inv = np.linalg.inv(X.T @ X) se_beta = np.sqrt(se_sq * XtX_inv[1, 1]) if se_beta > 0: t_stat = slope / se_beta # One-tailed: theory predicts positive slope p_val = t_dist.sf(t_stat, df=n - k) if slope > 0 else 1 - t_dist.sf(-t_stat, df=n - k) else: t_stat = float('nan') p_val = float('nan') sign_ok = slope > 0 verd1 = verdict(p_val, sign_ok) # R-squared ss_tot = np.sum((y - y.mean()) ** 2) r_sq = 1 - rss / ss_tot if ss_tot > 0 else float('nan') results.append( log_result( name='productivity_dispersion_K4', stat=float(slope), p=float(p_val), verd=verd1, details={ 'n_groups': n, 't_stat': float(t_stat), 'r_squared': float(r_sq), 'intercept': float(beta[0]), 'note': 'Theory: Var(logF) ~ K^4 => sigma ~ K^2. ' 'Log-log regression approximates polynomial K^2 ' 'relationship; beta < 2 expected due to functional ' 'form mismatch (polynomial vs power-law).', }, description=( "Tests Prediction 1.7: higher CES curvature K amplifies " "within-industry productivity dispersion (Lean: " "dispersion_increases_with_K). OLS: log(90/10 labor " "productivity ratio) ~ beta * log(|K|) across 62 four-digit " "NAICS groups, one-tailed. Data: NBER-CES vadd/emp for " "productivity, pay/vadd for labor share -> sigma -> K. " "Prediction: beta > 0 (p < 0.05)." ), ) ) # Test 2: Spearman rank correlation as nonparametric check from common.stats import spearman_test rho_s, p_s = spearman_test(avg['absK'].values, avg['ratio_90_10'].values) sign_ok_s = rho_s > 0 verd2 = verdict(p_s, sign_ok_s) results.append( log_result( name='productivity_dispersion_spearman', stat=float(rho_s), p=float(p_s), verd=verd2, details={'n_groups': len(avg)}, description=( "Nonparametric check: Spearman rank correlation of |K| vs " "90/10 labor productivity ratio across 62 four-digit NAICS " "groups. Data: NBER-CES manufacturing. Predicted: positive " "(higher K => wider dispersion)." ), ) ) # Test 3: Panel version (group-year level, not averaged) panel = rdf[rdf['K'].abs() > 1e-6].copy() panel = panel[panel['ratio_90_10'] > 1.0].copy() if len(panel) >= 30: y_p = np.log(panel['ratio_90_10'].values) x_p = np.log(panel['K'].abs().values) X_p = np.column_stack([np.ones(len(x_p)), x_p]) beta_p, _, _, _ = np.linalg.lstsq(X_p, y_p, rcond=None) slope_p = beta_p[1] y_hat_p = X_p @ beta_p resid_p = y_p - y_hat_p rss_p = np.sum(resid_p ** 2) n_p = len(y_p) se_sq_p = rss_p / (n_p - 2) XtX_inv_p = np.linalg.inv(X_p.T @ X_p) se_p = np.sqrt(se_sq_p * XtX_inv_p[1, 1]) if se_p > 0: t_p = slope_p / se_p # One-tailed: theory predicts positive slope p_p = t_dist.sf(t_p, df=n_p - 2) if slope_p > 0 else 1 - t_dist.sf(-t_p, df=n_p - 2) else: t_p = float('nan') p_p = float('nan') verd3 = verdict(p_p, slope_p > 0) results.append( log_result( name='productivity_dispersion_panel', stat=float(slope_p), p=float(p_p), verd=verd3, details={ 'n_obs': n_p, 't_stat': float(t_p), 'note': 'Unclustered SEs; t-stat may be inflated by ' 'within-group correlation.', }, description=( "Panel OLS (group-year level): log(90/10 labor " "productivity ratio) ~ beta * log(|K|), one-tailed. " "Data: NBER-CES 4-digit NAICS x year. Note: unclustered " "SEs (t-stat overstated). Prediction: beta > 0." ), ) ) # Test 4: Polynomial functional form — theory predicts log(ratio) = a + c*K^2 # This is the DIRECT theoretical form (not the log-log approximation) K2_vals = avg['absK'].values ** 2 y_poly = np.log(avg['ratio_90_10'].values) X_poly = np.column_stack([np.ones(len(K2_vals)), K2_vals]) try: beta_poly, _, _, _ = np.linalg.lstsq(X_poly, y_poly, rcond=None) slope_poly = beta_poly[1] n_poly = len(y_poly) y_hat_poly = X_poly @ beta_poly resid_poly = y_poly - y_hat_poly rss_poly = np.sum(resid_poly ** 2) ss_tot_poly = np.sum((y_poly - y_poly.mean()) ** 2) r_sq_poly = 1 - rss_poly / ss_tot_poly if ss_tot_poly > 0 else float('nan') se_sq_poly = rss_poly / (n_poly - 2) XtX_inv_poly = np.linalg.inv(X_poly.T @ X_poly) se_poly = np.sqrt(se_sq_poly * XtX_inv_poly[1, 1]) if se_poly > 0: t_poly = slope_poly / se_poly p_poly = t_dist.sf(t_poly, df=n_poly - 2) if slope_poly > 0 else 1 - t_dist.sf(-t_poly, df=n_poly - 2) else: t_poly = float('nan') p_poly = float('nan') # Compare AIC/BIC: polynomial vs log-log from common.stats import compute_aic_bic aic_loglog, bic_loglog = compute_aic_bic(rss, 2, n) aic_poly, bic_poly = compute_aic_bic(rss_poly, 2, n_poly) verd_poly = verdict(p_poly, slope_poly > 0) results.append( log_result( name='productivity_dispersion_polynomial', stat=float(slope_poly), p=float(p_poly), verd=verd_poly, details={ 'n_groups': n_poly, 't_stat': float(t_poly), 'r_squared': float(r_sq_poly), 'intercept': float(beta_poly[0]), 'aic_polynomial': float(aic_poly), 'bic_polynomial': float(bic_poly), 'aic_loglog': float(aic_loglog), 'bic_loglog': float(bic_loglog), 'note': 'Direct theoretical form: log(ratio) = a + c*K^2. ' 'Compare AIC/BIC with log-log form.', }, description=( "Tests the DIRECT theoretical functional form: " "log(90/10 ratio) = a + c*K^2 (from Var(logF) ~ K^4 => " "sigma ~ K^2). OLS across 4-digit NAICS groups, one-tailed. " "Data: NBER-CES. Compared to log-log via AIC/BIC. " "Prediction: c > 0 (p < 0.05)." ), ) ) _cache['polynomial'] = { 'K2': K2_vals, 'y': y_poly, 'slope': slope_poly, 'intercept': beta_poly[0], 'r_sq': r_sq_poly, 'p': p_poly, 'aic_poly': aic_poly, 'bic_poly': bic_poly, 'aic_loglog': aic_loglog, 'bic_loglog': bic_loglog, } except (np.linalg.LinAlgError, ValueError): pass # Test 5: Panel with clustered SEs (by naics4) if len(panel) >= 30: y_pc = np.log(panel['ratio_90_10'].values) x_pc = np.log(panel['K'].abs().values) X_pc = np.column_stack([np.ones(len(x_pc)), x_pc]) beta_pc, _, _, _ = np.linalg.lstsq(X_pc, y_pc, rcond=None) slope_pc = beta_pc[1] y_hat_pc = X_pc @ beta_pc resid_pc = y_pc - y_hat_pc n_pc = len(y_pc) # Cluster-robust SEs (CR1) by naics4 groups = panel['naics4'].values unique_groups = np.unique(groups) G = len(unique_groups) XtX_inv_pc = np.linalg.inv(X_pc.T @ X_pc) meat = np.zeros((2, 2)) for g in unique_groups: mask = groups == g Xg = X_pc[mask] eg = resid_pc[mask] score_g = Xg.T @ eg # (2,) meat += np.outer(score_g, score_g) # Small-sample correction: G/(G-1) * (n-1)/(n-k) correction = (G / (G - 1)) * ((n_pc - 1) / (n_pc - 2)) V_cluster = correction * XtX_inv_pc @ meat @ XtX_inv_pc se_cluster = np.sqrt(V_cluster[1, 1]) if se_cluster > 0: t_cluster = slope_pc / se_cluster p_cluster = t_dist.sf(t_cluster, df=G - 1) if slope_pc > 0 else 1 - t_dist.sf(-t_cluster, df=G - 1) else: t_cluster = float('nan') p_cluster = float('nan') verd_cluster = verdict(p_cluster, slope_pc > 0) results.append( log_result( name='productivity_dispersion_panel_clustered', stat=float(slope_pc), p=float(p_cluster), verd=verd_cluster, details={ 'n_obs': n_pc, 'n_clusters': int(G), 't_stat_clustered': float(t_cluster), 't_stat_unclustered': float(slope_pc / np.sqrt((np.sum(resid_pc ** 2) / (n_pc - 2)) * XtX_inv_pc[1, 1])), 'note': 'Clustered SEs by naics4 group. df = G-1.', }, description=( "Panel OLS with cluster-robust SEs: log(90/10 labor " "productivity ratio) ~ beta * log(|K|), clustered by " "4-digit NAICS group, one-tailed. Data: NBER-CES " "group-year panel. Corrects for within-group correlation " "that inflates unclustered t-stats. Prediction: beta > 0." ), ) ) # Test 6: Complements vs Substitutes split comp = avg[avg['K'] > 0].copy() # complements: rho < 0, K > 0 subs = avg[avg['K'] < 0].copy() # substitutes: rho > 0, K < 0 if len(comp) >= 5 and len(subs) >= 5: from scipy.stats import mannwhitneyu # Theory: complements should show MORE dispersion per unit |K| # Simple check: compare median dispersion ratio comp_med = comp['ratio_90_10'].median() subs_med = subs['ratio_90_10'].median() # Mann-Whitney: complements have higher 90/10 ratio? u_stat, p_mw_two = mannwhitneyu( comp['ratio_90_10'].values, subs['ratio_90_10'].values, alternative='greater' ) sign_ok_cs = comp_med > subs_med verd_cs = verdict(p_mw_two, sign_ok_cs) results.append( log_result( name='productivity_dispersion_comp_vs_subs', stat=float(comp_med / subs_med if subs_med > 0 else float('nan')), p=float(p_mw_two), verd=verd_cs, details={ 'n_complements': len(comp), 'n_substitutes': len(subs), 'median_ratio_complements': float(comp_med), 'median_ratio_substitutes': float(subs_med), 'mann_whitney_U': float(u_stat), }, description=( "Complements vs substitutes: compares 90/10 productivity " "ratio between complement sectors (K > 0, rho < 0) and " "substitute sectors (K < 0, rho > 0). Mann-Whitney U, " "one-tailed. Data: NBER-CES 4-digit NAICS groups. " "Theory predicts complementary sectors show MORE " "dispersion. Prediction: complement median > substitute." ), ) ) _cache['comp_vs_subs'] = { 'comp_ratios': comp['ratio_90_10'].values, 'subs_ratios': subs['ratio_90_10'].values, 'comp_med': comp_med, 'subs_med': subs_med, 'p': p_mw_two, } elif len(comp) >= 5 or len(subs) >= 5: # Only one regime has enough data — note but skip results.append( log_result( name='productivity_dispersion_comp_vs_subs', stat=float('nan'), p=float('nan'), verd='SKIPPED', description=( f"Complements vs substitutes: only {len(comp)} complement " f"and {len(subs)} substitute groups (need >= 5 each)." ), ) ) # Cache for figures _cache['x_logK'] = x # log(|K|) for cross-section _cache['y_logR'] = y # log(90/10 ratio) for cross-section _cache['slope'] = slope _cache['intercept'] = beta[0] _cache['r_sq'] = r_sq _cache['p_val'] = p_val _cache['avg'] = avg _cache['rho_s'] = rho_s if 'rho_s' in dir() else None _cache['p_s'] = p_s if 'p_s' in dir() else None return results def figures(): """Generate figures for productivity dispersion subtests.""" if not _cache: main() paths = [] x = _cache.get('x_logK') y = _cache.get('y_logR') if x is not None and y is not None: slope = _cache.get('slope', 0) intercept = _cache.get('intercept', 0) r_sq = _cache.get('r_sq', 0) p_val = _cache.get('p_val', 1) # ── Figure 1: OLS scatter (covers K4 and Spearman subtests) ── fig, ax = setup_figure(width=8, height=6) ax.scatter(x, y, alpha=0.5, s=30, color='steelblue', edgecolor='black', linewidth=0.5) x_line = np.linspace(x.min(), x.max(), 100) ax.plot(x_line, intercept + slope * x_line, 'r-', linewidth=2, label=f'OLS: β={slope:.3f} (p={p_val:.1e})') ax.set_xlabel('log(|K|)') ax.set_ylabel('log(90/10 productivity ratio)') ax.set_title('Productivity Dispersion vs CES Curvature') rho_s = _cache.get('rho_s') p_s = _cache.get('p_s') note = f'R² = {r_sq:.3f}' if rho_s is not None: note += f'\nSpearman r = {rho_s:.3f}' ax.annotate(note, xy=(0.05, 0.95), xycoords='axes fraction', fontsize=10, va='top', bbox=dict(boxstyle='round', fc='wheat', alpha=0.5)) ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'productivity_dispersion_scatter', subfolder='foundations')) # ── Figure 2: Polynomial form (K^2) vs log-log comparison ── poly = _cache.get('polynomial') if poly is not None: fig, (ax1, ax2) = setup_subplots(ncols=2, width=12, height=5) K2 = poly['K2'] y_p = poly['y'] sl_p = poly['slope'] int_p = poly['intercept'] # Left: polynomial fit ax1.scatter(K2, y_p, alpha=0.5, s=30, color='steelblue', edgecolor='black', linewidth=0.5) K2_line = np.linspace(K2.min(), K2.max(), 100) ax1.plot(K2_line, int_p + sl_p * K2_line, 'r-', linewidth=2, label=f'c={sl_p:.3f} (p={poly["p"]:.1e})') ax1.set_xlabel('K²') ax1.set_ylabel('log(90/10 productivity ratio)') ax1.set_title('Polynomial Form (Theory)') ax1.annotate(f'R² = {poly["r_sq"]:.3f}\nAIC = {poly["aic_poly"]:.1f}', xy=(0.05, 0.95), xycoords='axes fraction', fontsize=10, va='top', bbox=dict(boxstyle='round', fc='wheat', alpha=0.5)) ax1.legend(fontsize=9) ax1.grid(True, alpha=0.3) # Right: log-log fit for comparison if x is not None and y is not None: ax2.scatter(x, y, alpha=0.5, s=30, color='steelblue', edgecolor='black', linewidth=0.5) ax2.plot(x_line, intercept + slope * x_line, 'r-', linewidth=2, label=f'β={slope:.3f} (p={p_val:.1e})') ax2.set_xlabel('log(|K|)') ax2.set_ylabel('log(90/10 productivity ratio)') ax2.set_title('Log-Log Form (Approximation)') ax2.annotate(f'R² = {r_sq:.3f}\nAIC = {poly["aic_loglog"]:.1f}', xy=(0.05, 0.95), xycoords='axes fraction', fontsize=10, va='top', bbox=dict(boxstyle='round', fc='wheat', alpha=0.5)) ax2.legend(fontsize=9) ax2.grid(True, alpha=0.3) fig.tight_layout() paths.extend(save_figure(fig, 'productivity_dispersion_polynomial', subfolder='foundations')) # ── Figure 3: Complements vs Substitutes box plot ── cs = _cache.get('comp_vs_subs') if cs is not None: fig, ax = setup_figure(width=6, height=5) data = [cs['comp_ratios'], cs['subs_ratios']] bp = ax.boxplot(data, labels=['Complements\n(K > 0)', 'Substitutes\n(K < 0)'], patch_artist=True, widths=0.5) bp['boxes'][0].set_facecolor('steelblue') bp['boxes'][0].set_alpha(0.6) bp['boxes'][1].set_facecolor('coral') bp['boxes'][1].set_alpha(0.6) ax.set_ylabel('90/10 Productivity Ratio') ax.set_title('Dispersion: Complements vs Substitutes') ax.annotate( f'Median ratio: {cs["comp_med"]:.2f} vs {cs["subs_med"]:.2f}\n' f'Mann-Whitney p = {cs["p"]:.3f}', xy=(0.05, 0.95), xycoords='axes fraction', fontsize=10, va='top', bbox=dict(boxstyle='round', fc='wheat', alpha=0.5)) ax.grid(True, alpha=0.3, axis='y') paths.extend(save_figure(fig, 'productivity_dispersion_comp_subs', subfolder='foundations')) return paths if __name__ == '__main__': for r in main(): print(f" {r['name']}: {r['verdict']} (stat={r['statistic']:.4f}, p={r['p_value']:.4f})")

n obs	3782
t stat	36.1896
p value	1.661e-246
statistic	0.7141

Productivity Dispersion Panel

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps