#!/usr/bin/env python3 """ CES vs. Translog Convergence Horse Race ======================================== Tests predictions of the CES stability theorem (Paper 1, Thm 3.2): P1: rho preserved across aggregation levels P2: Translog interaction terms decay monotonically with aggregation P3: Log-log slope of interaction ratio vs agg factor is negative P4: Crossover — CES wins AIC/BIC at macro when J is large enough P5: FA-CES crossover — factor-augmented CES closes the crossover gap P6: Random re-grouping — random blocks give steeper slope than NAICS hierarchy Data: NBER-CES Manufacturing Industry Database (1958-2018, 2012 NAICS). The crossover test uses KLEM (Capital, Labor, Materials, Energy) with J=4 factors, where translog has 16 parameters vs CES's 6. This creates an AIC penalty gap of 20, making the crossover testable. The original J=2 test (K, L only) had a penalty gap of just 6, insufficient for crossover. Lean backing: - aggregation_invariance_of_rho: rho preserved exactly under aggregation (P1) - convergence_rate_cubic: dominant non-CES mode (m=3) decays under aggregation (P2) - aggregation_fixed_point: CES is the unique fixed point of the aggregation operator Review Notes (2026-03-11): Theory match: CORRECTED. The theorem's k^{-1/2} rate assumes random block averaging (CLT: independent draws within each block). NAICS hierarchical grouping creates highly correlated blocks (steel grouped with metals, not with food), reducing the effective block size. The -0.5 exponent is an UPPER BOUND for random grouping. For correlated grouping, the exponent should be between -0.5 and 0. The key testable prediction is that the interaction ratio decays monotonically (slope < 0), not the exact rate. Data-signal alignment: STRONG. NBER-CES provides consistent capital, labor, and value-added series for manufacturing at all NAICS levels (6->2 digit). Translog interaction ratio is a direct proxy for non-CES deviation amplitude. Known limitations: Only N=5 aggregation levels (5 data points for slope regression). Bootstrap CI with N=5 is of limited reliability. NAICS hierarchy is fixed; cannot test random re-groupings. Review Notes — Crossover subtest (2026-03-11, revised): The AIC crossover (P4) is an asymptotic inference from the convergence theorem. With J=2 inputs (K,L), translog's 3 extra parameters have a trivial AIC penalty (2*3=6), so even tiny non-CES deviations beat CES. Solution: use KLEM (J=4) where translog has 10 extra second-order terms (AIC penalty = 20), making it much harder for translog to justify its flexibility at macro levels where non-CES deviations have decayed. BIC (penalty = log(n)*k) provides an even stricter test. Migrated from scripts/test_convergence_horse_race.py """ import warnings import numpy as np import pandas as pd from scipy.optimize import minimize from scipy import stats from common.data import load_nber_ces from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure warnings.filterwarnings('ignore') YEAR_MIN = 1977 YEAR_MAX = 2018 MIN_OBS = 15 def _info_criteria(ss_res, n, k_params): """Compute AIC and BIC from residual sum of squares.""" aic = n * np.log(ss_res / n + 1e-30) + 2 * k_params bic = n * np.log(ss_res / n + 1e-30) + np.log(n) * k_params return aic, bic def fit_translog_2(ln_y, ln_k, ln_l): n = len(ln_y) t = np.arange(n, dtype=float) t = (t - t.mean()) / (t.std() + 1e-10) X = np.column_stack([np.ones(n), t, ln_k, ln_l, 0.5*ln_k**2, 0.5*ln_l**2, ln_k*ln_l]) b = np.linalg.lstsq(X, ln_y, rcond=None)[0] ss_res = np.sum((ln_y - X @ b)**2) ss_tot = np.sum((ln_y - ln_y.mean())**2) R2 = 1 - ss_res/ss_tot if ss_tot > 1e-20 else np.nan k_params = 7 aic, bic = _info_criteria(ss_res, n, k_params) alpha_norm = np.sqrt(b[2]**2 + b[3]**2) beta_norm = np.sqrt(b[4]**2 + b[5]**2 + b[6]**2) ir = beta_norm / alpha_norm if alpha_norm > 1e-10 else np.nan return {'R2': R2, 'aic': aic, 'bic': bic, 'ir': ir} def fit_ces_2(ln_y, ln_k, ln_l): n = len(ln_y) t = np.arange(n, dtype=float) t = (t - t.mean()) / (t.std() + 1e-10) def sse(params): ln_A, gamma, w_logit, rho = params w = 1.0 / (1.0 + np.exp(-np.clip(w_logit, -10, 10))) rho_c = np.clip(rho, -5, 0.99) if abs(rho_c) < 1e-6: pred = ln_A + gamma*t + w*ln_k + (1-w)*ln_l else: inner = w*np.exp(rho_c*ln_k) + (1-w)*np.exp(rho_c*ln_l) inner = np.maximum(inner, 1e-30) pred = ln_A + gamma*t + (1.0/rho_c)*np.log(inner) return np.sum((ln_y - pred)**2) best = None m = np.mean(ln_y) for rho0 in [-3.0, -1.5, -0.5, -0.1, 0.01, 0.3, 0.6, 0.9]: for w0 in [-1.0, 0.0, 1.0]: try: res = minimize(sse, [m, 0.0, w0, rho0], method='Nelder-Mead', options={'maxiter': 10000, 'xatol': 1e-10, 'fatol': 1e-13}) if best is None or res.fun < best.fun: best = res except Exception: pass if best is None: return {'R2': np.nan, 'aic': np.nan, 'bic': np.nan, 'rho': np.nan} ss_res = best.fun ss_tot = np.sum((ln_y - ln_y.mean())**2) R2 = 1 - ss_res/ss_tot if ss_tot > 1e-20 else np.nan k_params = 4 aic, bic = _info_criteria(ss_res, n, k_params) rho = np.clip(best.x[3], -5, 0.99) return {'R2': R2, 'aic': aic, 'bic': bic, 'rho': rho} def fit_translog_J(ln_y, ln_inputs): """Fit J-factor translog: J linear + J(J+1)/2 quadratic + intercept + trend.""" n, J = ln_inputs.shape t = np.arange(n, dtype=float) t = (t - t.mean()) / (t.std() + 1e-10) # Build regressor matrix: intercept, trend, J linear, J(J+1)/2 quadratic cols = [np.ones(n), t] for j in range(J): cols.append(ln_inputs[:, j]) for j in range(J): for k in range(j, J): if j == k: cols.append(0.5 * ln_inputs[:, j] ** 2) else: cols.append(ln_inputs[:, j] * ln_inputs[:, k]) X = np.column_stack(cols) k_params = X.shape[1] # 2 + J + J(J+1)/2 b = np.linalg.lstsq(X, ln_y, rcond=None)[0] ss_res = np.sum((ln_y - X @ b) ** 2) ss_tot = np.sum((ln_y - ln_y.mean()) ** 2) R2 = 1 - ss_res / ss_tot if ss_tot > 1e-20 else np.nan aic, bic = _info_criteria(ss_res, n, k_params) # Interaction ratio: ||quadratic coeffs|| / ||linear coeffs|| alpha_norm = np.linalg.norm(b[2:2 + J]) beta_norm = np.linalg.norm(b[2 + J:]) ir = beta_norm / alpha_norm if alpha_norm > 1e-10 else np.nan return {'R2': R2, 'aic': aic, 'bic': bic, 'ir': ir, 'k_params': k_params} def fit_ces_J(ln_y, ln_inputs): """Fit J-factor CES: ln_A, gamma, (J-1) weight logits, rho.""" n, J = ln_inputs.shape t = np.arange(n, dtype=float) t = (t - t.mean()) / (t.std() + 1e-10) def sse(params): ln_A = params[0] gamma = params[1] rho = np.clip(params[2], -5, 0.99) # Softmax for weights logits = np.zeros(J) logits[:J - 1] = params[3:] w = np.exp(logits - logits.max()) w = w / w.sum() if abs(rho) < 1e-6: pred = ln_A + gamma * t + np.sum(w * ln_inputs, axis=1) else: inner = np.sum(w * np.exp(rho * ln_inputs), axis=1) inner = np.maximum(inner, 1e-30) pred = ln_A + gamma * t + (1.0 / rho) * np.log(inner) return np.sum((ln_y - pred) ** 2) best = None m = np.mean(ln_y) for rho0 in [-2.0, -0.5, -0.1, 0.01, 0.3, 0.6]: init = [m, 0.0, rho0] + [0.0] * (J - 1) try: res = minimize(sse, init, method='Nelder-Mead', options={'maxiter': 20000, 'xatol': 1e-10, 'fatol': 1e-13}) if best is None or res.fun < best.fun: best = res except Exception: pass if best is None: return {'R2': np.nan, 'aic': np.nan, 'bic': np.nan, 'rho': np.nan} ss_res = best.fun ss_tot = np.sum((ln_y - ln_y.mean()) ** 2) R2 = 1 - ss_res / ss_tot if ss_tot > 1e-20 else np.nan k_params = 3 + (J - 1) # ln_A, gamma, rho, (J-1) weights aic, bic = _info_criteria(ss_res, n, k_params) rho = np.clip(best.x[2], -5, 0.99) return {'R2': R2, 'aic': aic, 'bic': bic, 'rho': rho, 'k_params': k_params} def fit_fa_ces_J(ln_y, ln_inputs): """Fit J-factor factor-augmented CES with input-specific technical change. Model: ln Y = ln_A + (1/rho) * ln[sum_j w_j * exp(rho * (ln x_j + gamma_j * t))] Parameters: ln_A, rho, (J-1) weight logits, J augmentation rates = 2J+1 total. For KLEM (J=4): 9 params vs translog's 16. """ n, J = ln_inputs.shape t = np.arange(n, dtype=float) t = (t - t.mean()) / (t.std() + 1e-10) def sse(params): ln_A = params[0] rho = np.clip(params[1], -5, 0.99) logits = np.zeros(J) logits[:J - 1] = params[2:J + 1] logits -= logits.max() w = np.exp(logits) w /= w.sum() gammas = params[J + 1:2 * J + 1] aug = ln_inputs + np.outer(t, gammas) if abs(rho) < 1e-6: pred = ln_A + aug @ w else: inner = np.exp(rho * aug) @ w inner = np.maximum(inner, 1e-30) pred = ln_A + np.log(inner) / rho diff = ln_y - pred return diff @ diff best = None m = np.mean(ln_y) bounds = [(None, None), (-5.0, 0.99)] + [(None, None)] * (J - 1) + [(None, None)] * J for rho0 in [-1.0, 0.0, 0.5]: init = [m, rho0] + [0.0] * (J - 1) + [0.0] * J try: res = minimize(sse, init, method='L-BFGS-B', bounds=bounds, options={'maxiter': 2000, 'ftol': 1e-12}) if best is None or res.fun < best.fun: best = res except Exception: pass if best is None: return {'R2': np.nan, 'aic': np.nan, 'bic': np.nan, 'rho': np.nan} ss_res = best.fun ss_tot = np.sum((ln_y - ln_y.mean()) ** 2) R2 = 1 - ss_res / ss_tot if ss_tot > 1e-20 else np.nan k_params = 2 * J + 1 aic, bic = _info_criteria(ss_res, n, k_params) rho = np.clip(best.x[1], -5, 0.99) return {'R2': R2, 'aic': aic, 'bic': bic, 'rho': rho, 'k_params': k_params} def aggregate(df_in, level_col, vars_to_sum): agg_dict = {v: (v, 'sum') for v in vars_to_sum} agg_dict['n_sub'] = ('n6', 'nunique') return df_in.groupby([level_col, 'year']).agg(**agg_dict).reset_index().rename( columns={level_col: 'code'}) def estimate_level(agg_df, use_klem=False): """Estimate CES vs translog at each industry within an aggregation level. Parameters ---------- agg_df : DataFrame with columns code, year, and input variables use_klem : if True, use 4-factor KLEM (gross output); else 2-factor KL (value added) """ results = [] for code in agg_df['code'].unique(): sub = agg_df[agg_df['code'] == code].sort_values('year') if len(sub) < MIN_OBS: continue if use_klem: ln_y = np.log(sub['vship'].values) ln_inputs = np.column_stack([ np.log(sub['cap'].values), np.log(sub['prodh'].values), np.log(sub['real_mat'].values), np.log(sub['real_energy'].values), ]) tl = fit_translog_J(ln_y, ln_inputs) ces = fit_ces_J(ln_y, ln_inputs) else: ln_y = np.log(sub['vadd'].values) ln_k = np.log(sub['cap'].values) ln_l = np.log(sub['prodh'].values) tl = fit_translog_2(ln_y, ln_k, ln_l) ces = fit_ces_2(ln_y, ln_k, ln_l) results.append({ 'code': code, 'R2_ces': ces['R2'], 'aic_ces': ces['aic'], 'bic_ces': ces['bic'], 'R2_tl': tl['R2'], 'aic_tl': tl['aic'], 'bic_tl': tl['bic'], 'ir': tl['ir'], 'rho': ces['rho'], 'daic': tl['aic'] - ces['aic'], 'dbic': tl['bic'] - ces['bic'], }) return pd.DataFrame(results) def _random_regrouping(df, n_reps=100): """Test interaction ratio decay with random block averaging (2-factor K,L). Instead of NAICS hierarchy (correlated blocks), randomly assign 6-digit industries into blocks. Theory predicts slope = -0.5 for i.i.d. blocks (CLT). """ k_values = [1, 2, 4, 8, 16, 32, 64] good_codes = [] for code in df['n6'].unique(): sub = df[df['n6'] == code] if len(sub) >= MIN_OBS: good_codes.append(code) good_codes = sorted(good_codes) n_ind = len(good_codes) all_years = sorted(df['year'].unique()) n_years = len(all_years) year_to_idx = {y: i for i, y in enumerate(all_years)} # Pre-compute: (n_ind, n_years, 3) array for vadd, cap, prodh data_matrix = np.zeros((n_ind, n_years, 3)) valid_matrix = np.zeros((n_ind, n_years), dtype=bool) for i, code in enumerate(good_codes): sub = df[df['n6'] == code].sort_values('year') idxs = [year_to_idx[y] for y in sub['year'].values] data_matrix[i, idxs, 0] = sub['vadd'].values data_matrix[i, idxs, 1] = sub['cap'].values data_matrix[i, idxs, 2] = sub['prodh'].values valid_matrix[i, idxs] = True rng = np.random.RandomState(42) results = [] for k in k_values: n_groups = n_ind // k if n_groups < 3: continue if k == 1: # Baseline: individual industries, one "rep" irs = [] for i in range(n_ind): valid = valid_matrix[i] vals = data_matrix[i, valid] if len(vals) < MIN_OBS or np.any(vals <= 0): continue tl = fit_translog_2(np.log(vals[:, 0]), np.log(vals[:, 1]), np.log(vals[:, 2])) if np.isfinite(tl['ir']): irs.append(tl['ir']) if irs: results.append({'k': 1, 'ir_median': np.median(irs), 'ir_lo': np.median(irs), 'ir_hi': np.median(irs)}) continue ir_medians = [] for _ in range(n_reps): perm = rng.permutation(n_ind) irs = [] for g in range(n_groups): indices = perm[g * k:(g + 1) * k] group_valid = valid_matrix[indices].all(axis=0) if group_valid.sum() < MIN_OBS: continue agg = data_matrix[indices].sum(axis=0) vals = agg[group_valid] if np.any(vals <= 0): continue tl = fit_translog_2(np.log(vals[:, 0]), np.log(vals[:, 1]), np.log(vals[:, 2])) if np.isfinite(tl['ir']): irs.append(tl['ir']) if len(irs) >= 3: ir_medians.append(np.median(irs)) if ir_medians: results.append({ 'k': k, 'ir_median': np.median(ir_medians), 'ir_lo': np.percentile(ir_medians, 2.5), 'ir_hi': np.percentile(ir_medians, 97.5), }) return pd.DataFrame(results) _cache = {} def main(): try: df = load_nber_ces() except FileNotFoundError: return [log_result('convergence_horse_race', np.nan, np.nan, 'SKIPPED', description="CES vs translog convergence horse race using NBER-CES manufacturing data. " "Tests predictions of the CES stability theorem across NAICS aggregation levels.")] df = df[(df['year'] >= YEAR_MIN) & (df['year'] <= YEAR_MAX)].copy() df['naics'] = df['naics'].astype(str).str.strip() for v in ['vadd', 'cap', 'prodh']: df = df[df[v] > 0] # --- 2-factor (K, L) analysis --- df['n6'] = df['naics'].str[:6] df['n5'] = df['naics'].str[:5] df['n4'] = df['naics'].str[:4] df['n3'] = df['naics'].str[:3] df['n2'] = df['naics'].str[:2] n6_count = df['n6'].nunique() levels = [(6, 'n6'), (5, 'n5'), (4, 'n4'), (3, 'n3'), (2, 'n2')] summary_rows = [] for digits, col in levels: agg = aggregate(df, col, ['vadd', 'cap', 'prodh']) ind_results = estimate_level(agg) if len(ind_results) == 0: continue n_codes = agg['code'].nunique() k_agg = n6_count / n_codes summary_rows.append({ 'digits': digits, 'n_estimated': len(ind_results), 'k_agg': k_agg, 'ir_median': ind_results['ir'].median(), 'rho_median': ind_results['rho'].median(), 'rho_iqr': ind_results['rho'].quantile(0.75) - ind_results['rho'].quantile(0.25), 'pct_ces_wins_aic': (ind_results['daic'] > 0).mean() * 100, 'pct_ces_wins_bic': (ind_results['dbic'] > 0).mean() * 100, 'daic_median': ind_results['daic'].median(), 'dbic_median': ind_results['dbic'].median(), }) if len(summary_rows) < 3: return [log_result('convergence_horse_race', np.nan, np.nan, 'SKIPPED', description="CES vs translog convergence horse race using NBER-CES manufacturing data. " "Tests predictions of the CES stability theorem across NAICS aggregation levels.")] summary = pd.DataFrame(summary_rows) results = [] # P1: rho preserved rho_vals = summary['rho_median'].values rho_cv = np.std(rho_vals) / abs(np.mean(rho_vals)) if abs(np.mean(rho_vals)) > 0.01 else np.nan v1 = 'CONSISTENT' if (np.isfinite(rho_cv) and rho_cv < 0.15) else 'AMBIGUOUS' results.append(log_result('convergence_horse_race_rho_preserved', rho_cv, np.nan, v1, description="Tests whether CES rho is preserved across NAICS aggregation levels (6-digit to 2-digit). " "Coefficient of variation of median rho across levels; CV < 0.15 is consistent with aggregation invariance.")) # P3: convergence rate (slope < 0) log_k = np.log(summary['k_agg'].values) log_ir = np.log(summary['ir_median'].values) slope, intercept, r_value, p_value, std_err = stats.linregress(log_k, log_ir) n_pts = len(log_k) # One-sided test: slope < 0 (theory predicts negative) t_stat_neg = slope / std_err if std_err > 0 else np.nan p_negative = stats.t.cdf(t_stat_neg, n_pts - 2) if n_pts > 2 else np.nan # Bootstrap CI np.random.seed(42) slopes_boot = [] for _ in range(5000): idx = np.random.choice(n_pts, n_pts, replace=True) if len(set(idx)) < 2: continue s, _, _, _, _ = stats.linregress(log_k[idx], log_ir[idx]) slopes_boot.append(s) slopes_boot = np.array(slopes_boot) ci_lo, ci_hi = np.percentile(slopes_boot, [2.5, 97.5]) if slope < 0 and p_negative < 0.05 and -0.5 <= slope: v5 = 'CONSISTENT' elif slope < 0 and p_negative < 0.10: v5 = 'AMBIGUOUS' else: v5 = 'INCONSISTENT' results.append(log_result('convergence_horse_race_slope', slope, p_negative, v5, description="Tests whether the log-log slope of translog interaction ratio vs aggregation " "factor is negative (non-CES deviations decay under aggregation). " "Theory predicts slope = -0.5 for random block averaging (Thm 3.2); NAICS " "hierarchical grouping creates correlated blocks, yielding a flatter slope " "between -0.5 and 0. The key prediction is monotone decay (slope < 0), not " "the exact rate. OLS regression on NBER-CES data across 5 NAICS levels.", details={'ci_lo': ci_lo, 'ci_hi': ci_hi, 'theoretical_upper_bound': -0.5, 'r_squared': r_value**2})) _cache['summary'] = summary _cache['log_k'] = log_k _cache['log_ir'] = log_ir _cache['slope'] = slope _cache['intercept'] = intercept _cache['ci_lo'] = ci_lo _cache['ci_hi'] = ci_hi # --- P4: KLEM crossover (4-factor: K, L, M, E) --- # With J=4, translog has 16 params vs CES's 6 — AIC penalty gap = 20. # This makes the crossover testable: translog's overhead is large enough # that at macro levels (where non-CES deviations have decayed) CES wins. # Use gross output (vship) as Y, since KLEM includes intermediate inputs. df_klem = df.copy() # Deflate materials and energy to get real quantities for v in ['vship', 'matcost', 'energy', 'pimat', 'pien']: df_klem = df_klem[df_klem[v] > 0] df_klem['real_mat'] = df_klem['matcost'] / df_klem['pimat'] df_klem['real_energy'] = df_klem['energy'] / df_klem['pien'] klem_vars = ['vship', 'cap', 'prodh', 'real_mat', 'real_energy'] klem_summary_rows = [] for digits, col in levels: agg = aggregate(df_klem, col, klem_vars) ind_results = estimate_level(agg, use_klem=True) if len(ind_results) == 0: continue # Cache per-industry translog results for FA-CES reuse _cache.setdefault('klem_tl_cache', {})[digits] = ind_results[['code', 'aic_tl', 'bic_tl']].copy() n_codes = agg['code'].nunique() k_agg = n6_count / n_codes klem_summary_rows.append({ 'digits': digits, 'n_estimated': len(ind_results), 'k_agg': k_agg, 'ir_median': ind_results['ir'].median(), 'pct_ces_wins_aic': (ind_results['daic'] > 0).mean() * 100, 'pct_ces_wins_bic': (ind_results['dbic'] > 0).mean() * 100, 'daic_median': ind_results['daic'].median(), 'dbic_median': ind_results['dbic'].median(), }) if len(klem_summary_rows) >= 3: klem_summary = pd.DataFrame(klem_summary_rows) _cache['klem_summary'] = klem_summary # Primary crossover test: does CES win rate increase from micro to macro? pct_aic = klem_summary['pct_ces_wins_aic'].values pct_bic = klem_summary['pct_ces_wins_bic'].values aic_diff = pct_aic[-1] - pct_aic[0] bic_diff = pct_bic[-1] - pct_bic[0] # Use the stronger criterion (BIC, which penalizes more) # Crossover: CES wins >50% at macro AND <50% at micro bic_crossover = pct_bic[0] < 50 and pct_bic[-1] > 50 aic_crossover = pct_aic[0] < 50 and pct_aic[-1] > 50 # Direction: CES wins more at macro than micro bic_direction = bic_diff > 0 aic_direction = aic_diff > 0 if bic_crossover or aic_crossover: v4 = 'CONSISTENT' elif bic_direction or aic_direction: v4 = 'AMBIGUOUS' else: v4 = 'INCONSISTENT' results.append(log_result('convergence_horse_race_crossover', bic_diff, np.nan, v4, description="Tests for an information-criterion crossover using KLEM production " "(J=4: capital, labor, materials, energy). With 4 factors, translog has 16 " "parameters vs CES's 6, creating an AIC penalty gap of 20 (BIC gap even larger). " "Theory: non-CES deviations decay under aggregation, so CES should win at macro " "levels where translog's extra flexibility is not worth the parameter cost. " "BIC helps at micro (14% CES wins vs 2% AIC) but translog still dominates at " "all aggregation levels. INCONSISTENT because CES win rate decreases with " "aggregation (wrong direction). Root cause: aggregation reduces noise " "proportionally to non-CES deviations in correlated NAICS hierarchies.", details={ 'J': 4, 'translog_params': 16, # 2 + J + J(J+1)/2 = 2+4+10 'ces_params': 6, # ln_A, gamma, rho, (J-1) weights 'pct_ces_wins_aic_micro': float(pct_aic[0]), 'pct_ces_wins_aic_macro': float(pct_aic[-1]), 'pct_ces_wins_bic_micro': float(pct_bic[0]), 'pct_ces_wins_bic_macro': float(pct_bic[-1]), 'daic_micro': float(klem_summary['daic_median'].iloc[0]), 'daic_macro': float(klem_summary['daic_median'].iloc[-1]), 'dbic_micro': float(klem_summary['dbic_median'].iloc[0]), 'dbic_macro': float(klem_summary['dbic_median'].iloc[-1]), })) else: results.append(log_result('convergence_horse_race_crossover', np.nan, np.nan, 'SKIPPED', description="KLEM crossover test skipped — insufficient KLEM data at 3+ aggregation levels.")) # --- P5: FA-CES crossover (factor-augmented CES vs translog, KLEM) --- # Reuses cached translog AIC/BIC from KLEM section to avoid redundant fits klem_tl_cache = _cache.get('klem_tl_cache', {}) if len(klem_summary_rows) >= 3: fa_summary_rows = [] for digits, col in levels: cached_tl = klem_tl_cache.get(digits) if cached_tl is None: continue tl_lookup = dict(zip(cached_tl['code'], zip(cached_tl['aic_tl'], cached_tl['bic_tl']))) agg = aggregate(df_klem, col, klem_vars) fa_res = [] for code in agg['code'].unique(): if code not in tl_lookup: continue tl_aic, tl_bic = tl_lookup[code] sub = agg[agg['code'] == code].sort_values('year') if len(sub) < MIN_OBS: continue ln_y = np.log(sub['vship'].values) ln_inputs = np.column_stack([ np.log(sub['cap'].values), np.log(sub['prodh'].values), np.log(sub['real_mat'].values), np.log(sub['real_energy'].values), ]) fa = fit_fa_ces_J(ln_y, ln_inputs) fa_res.append({ 'daic_fa': tl_aic - fa['aic'], 'dbic_fa': tl_bic - fa['bic'], }) if fa_res: fa_df = pd.DataFrame(fa_res) n_codes = agg['code'].nunique() k_agg = n6_count / n_codes fa_summary_rows.append({ 'digits': digits, 'k_agg': k_agg, 'n_estimated': len(fa_df), 'pct_fa_wins_aic': (fa_df['daic_fa'] > 0).mean() * 100, 'pct_fa_wins_bic': (fa_df['dbic_fa'] > 0).mean() * 100, 'daic_fa_median': fa_df['daic_fa'].median(), 'dbic_fa_median': fa_df['dbic_fa'].median(), }) if len(fa_summary_rows) >= 3: fa_summary = pd.DataFrame(fa_summary_rows) _cache['fa_summary'] = fa_summary pct_aic_fa = fa_summary['pct_fa_wins_aic'].values pct_bic_fa = fa_summary['pct_fa_wins_bic'].values bic_cross = pct_bic_fa[0] < 50 and pct_bic_fa[-1] > 50 aic_cross = pct_aic_fa[0] < 50 and pct_aic_fa[-1] > 50 bic_dir = pct_bic_fa[-1] > pct_bic_fa[0] aic_dir = pct_aic_fa[-1] > pct_aic_fa[0] if bic_cross or aic_cross: v_fa = 'CONSISTENT' elif bic_dir or aic_dir: v_fa = 'AMBIGUOUS' else: v_fa = 'INCONSISTENT' bic_diff_fa = pct_bic_fa[-1] - pct_bic_fa[0] results.append(log_result('convergence_horse_race_fa_crossover', bic_diff_fa, np.nan, v_fa, description="Tests whether factor-augmented CES (input-specific technical change) " "achieves the crossover that standard CES fails. FA-CES has 2J+1=9 parameters " "for KLEM vs translog's 16. If FA-CES captures the factor augmentation " "heterogeneity that standard CES misses, its win rate should increase from " "micro to macro levels.", details={ 'fa_params': 2 * 4 + 1, 'translog_params': 16, 'pct_fa_wins_aic_micro': float(pct_aic_fa[0]), 'pct_fa_wins_aic_macro': float(pct_aic_fa[-1]), 'pct_fa_wins_bic_micro': float(pct_bic_fa[0]), 'pct_fa_wins_bic_macro': float(pct_bic_fa[-1]), 'daic_fa_micro': float(fa_summary['daic_fa_median'].iloc[0]), 'daic_fa_macro': float(fa_summary['daic_fa_median'].iloc[-1]), 'dbic_fa_micro': float(fa_summary['dbic_fa_median'].iloc[0]), 'dbic_fa_macro': float(fa_summary['dbic_fa_median'].iloc[-1]), })) else: results.append(log_result('convergence_horse_race_fa_crossover', np.nan, np.nan, 'SKIPPED', description="FA-CES crossover test skipped — insufficient data.")) else: results.append(log_result('convergence_horse_race_fa_crossover', np.nan, np.nan, 'SKIPPED', description="FA-CES crossover test skipped — no KLEM data.")) # --- P6: Random re-grouping slope --- rand_df = _random_regrouping(df) if len(rand_df) >= 3: _cache['rand_regrouping'] = rand_df log_k_r = np.log(rand_df['k'].values) log_ir_r = np.log(rand_df['ir_median'].values) slope_r, intercept_r, r_r, p_r, se_r = stats.linregress(log_k_r, log_ir_r) # One-sided test: slope < 0 n_pts_r = len(log_k_r) t_stat_r = slope_r / se_r if se_r > 0 else np.nan p_neg_r = stats.t.cdf(t_stat_r, n_pts_r - 2) if n_pts_r > 2 else np.nan _cache['rand_slope'] = slope_r _cache['rand_intercept'] = intercept_r # Verdict: random slope should be (1) significantly negative and # (2) steeper than NAICS slope (less correlation => faster decay) naics_slope_val = _cache.get('slope', 0) if slope_r < 0 and p_neg_r < 0.05 and slope_r < naics_slope_val: v_rand = 'CONSISTENT' elif slope_r < 0 and p_neg_r < 0.10: v_rand = 'AMBIGUOUS' else: v_rand = 'INCONSISTENT' results.append(log_result('convergence_horse_race_random_slope', slope_r, p_neg_r, v_rand, description="Tests the k^{-1/2} prediction using random block averaging instead of " "NAICS hierarchy. Random blocks should give steeper slope than NAICS (less " "within-block correlation => faster CLT convergence). Theory predicts -0.5 " "for i.i.d. blocks; NAICS gives -0.149 due to hierarchical correlation. " "100 replications per block size (k=1,2,4,8,16,32,64).", details={ 'naics_slope': float(naics_slope_val), 'theoretical_slope': -0.5, 'r_squared': r_r ** 2, 'n_k_values': len(rand_df), })) else: results.append(log_result('convergence_horse_race_random_slope', np.nan, np.nan, 'SKIPPED', description="Random re-grouping test skipped — insufficient data.")) return results def figures(): """Generate per-subtest figures for convergence horse race.""" if not _cache: main() paths = [] summary = _cache.get('summary') if summary is None: return paths import matplotlib.pyplot as plt digits = summary['digits'].values # --- Figure 1: rho preserved across aggregation levels --- fig, ax = setup_figure(width=7, height=5) ax.plot(digits, summary['rho_median'].values, 'o-', c='darkorange', markersize=10, linewidth=2) lo = summary['rho_median'] - summary['rho_iqr'] / 2 hi = summary['rho_median'] + summary['rho_iqr'] / 2 ax.fill_between(digits, lo, hi, alpha=0.2, color='darkorange') rho_cv = np.std(summary['rho_median'].values) / abs(np.mean(summary['rho_median'].values)) ax.set_xlabel('NAICS digits (finer → coarser)', fontsize=11) ax.set_ylabel('Median CES rho', fontsize=11) ax.set_title(f'CES rho preserved across aggregation levels (CV={rho_cv:.3f})', fontsize=12) ax.invert_xaxis() ax.set_xticks(digits) paths.extend(save_figure(fig, 'convergence_horse_race_rho_preserved', subfolder='foundations')) # --- Figure 2: log-log convergence slope --- log_k = _cache['log_k'] log_ir = _cache['log_ir'] slope = _cache['slope'] intercept = _cache['intercept'] ci_lo = _cache['ci_lo'] ci_hi = _cache['ci_hi'] fig, ax = setup_figure(width=7, height=5) ax.scatter(log_k, log_ir, c='steelblue', edgecolors='k', s=80, zorder=3) for i, d in enumerate(digits): ax.annotate(f'{d}-digit', (log_k[i], log_ir[i]), textcoords='offset points', xytext=(8, 5), fontsize=9) xr = np.linspace(log_k.min() - 0.2, log_k.max() + 0.2, 50) ax.plot(xr, intercept + slope * xr, 'r-', linewidth=2, label=f'OLS: slope={slope:.3f} [{ci_lo:.3f}, {ci_hi:.3f}]') ax.plot(xr, intercept + (-0.5) * xr, 'k--', alpha=0.5, linewidth=1.5, label='Theory (random blocks): slope=-0.50') ax.set_xlabel('log(aggregation factor k)', fontsize=11) ax.set_ylabel('log(interaction ratio ||beta||/||alpha||)', fontsize=11) ax.set_title('Non-CES deviation decay under aggregation', fontsize=12) ax.legend(fontsize=9, loc='upper right') paths.extend(save_figure(fig, 'convergence_horse_race_slope', subfolder='foundations')) # --- Figure 3: KLEM crossover (4-factor AIC/BIC) --- klem_summary = _cache.get('klem_summary') if klem_summary is not None: klem_digits = klem_summary['digits'].values fig, axes = setup_subplots(1, 2, width=14, height=5) # (a) CES win rate: AIC vs BIC across levels ax = axes[0] x = np.arange(len(klem_digits)) w = 0.35 aic_vals = klem_summary['pct_ces_wins_aic'].values bic_vals = klem_summary['pct_ces_wins_bic'].values ax.bar(x - w/2, aic_vals, w, label='AIC', color='steelblue', edgecolor='k') ax.bar(x + w/2, bic_vals, w, label='BIC', color='darkorange', edgecolor='k') ax.axhline(50, color='k', linestyle='--', alpha=0.5, label='50% crossover') ax.set_xticks(x) ax.set_xticklabels([f'{d}-digit' for d in klem_digits]) ax.set_xlabel('NAICS level (finer -> coarser)', fontsize=11) ax.set_ylabel('% industries where CES wins', fontsize=11) ax.set_title('(a) KLEM (J=4): CES win rate by level', fontsize=11) ax.set_ylim(0, 100) ax.legend(fontsize=9) # (b) Median delta-IC by level ax = axes[1] ax.plot(klem_digits, klem_summary['daic_median'].values, 's-', c='steelblue', markersize=8, linewidth=2, label='delta-AIC') ax.plot(klem_digits, klem_summary['dbic_median'].values, 'o-', c='darkorange', markersize=8, linewidth=2, label='delta-BIC') ax.axhline(0, color='k', linestyle='--', alpha=0.5) ax.set_xlabel('NAICS digits', fontsize=11) ax.set_ylabel('Median delta-IC (translog - CES)', fontsize=11) ax.set_title('(b) KLEM: information criterion gap', fontsize=11) ax.invert_xaxis() ax.set_xticks(klem_digits) ax.legend(fontsize=9) # Shade regions for d_vals, color in [(klem_summary['daic_median'].values, 'steelblue'), (klem_summary['dbic_median'].values, 'darkorange')]: ax.fill_between(klem_digits, 0, d_vals, where=np.array(d_vals) > 0, alpha=0.1, color='forestgreen') ax.fill_between(klem_digits, 0, d_vals, where=np.array(d_vals) <= 0, alpha=0.1, color='tomato') fig.suptitle('KLEM Crossover: CES vs Translog (J=4, 16 vs 6 params)', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'convergence_horse_race_crossover', subfolder='foundations')) else: # Fallback: 2-factor AIC/BIC comparison fig, axes = setup_subplots(1, 2, width=13, height=5) ax = axes[0] x = np.arange(len(digits)) w = 0.35 aic_vals = summary['pct_ces_wins_aic'].values bic_vals = summary['pct_ces_wins_bic'].values ax.bar(x - w/2, aic_vals, w, label='AIC', color='steelblue', edgecolor='k') ax.bar(x + w/2, bic_vals, w, label='BIC', color='darkorange', edgecolor='k') ax.axhline(50, color='k', linestyle='--', alpha=0.5, label='50% crossover') ax.set_xticks(x) ax.set_xticklabels([f'{d}-digit' for d in digits]) ax.set_xlabel('NAICS level', fontsize=11) ax.set_ylabel('% industries where CES wins', fontsize=11) ax.set_title('(a) J=2: CES win rate (AIC vs BIC)', fontsize=11) ax.set_ylim(0, 100) ax.legend(fontsize=9) ax = axes[1] ax.plot(digits, summary['daic_median'].values, 's-', c='steelblue', markersize=8, linewidth=2, label='delta-AIC') ax.plot(digits, summary['dbic_median'].values, 'o-', c='darkorange', markersize=8, linewidth=2, label='delta-BIC') ax.axhline(0, color='k', linestyle='--', alpha=0.5) ax.set_xlabel('NAICS digits', fontsize=11) ax.set_ylabel('Median delta-IC (translog - CES)', fontsize=11) ax.set_title('(b) J=2: AIC and BIC gap', fontsize=11) ax.invert_xaxis() ax.set_xticks(digits) ax.legend(fontsize=9) fig.suptitle('Crossover Test (J=2): CES vs Translog', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'convergence_horse_race_crossover', subfolder='foundations')) # --- Figure 4: FA-CES crossover --- fa_summary = _cache.get('fa_summary') klem_summary_fig = _cache.get('klem_summary') if fa_summary is not None and klem_summary_fig is not None: fa_digits = fa_summary['digits'].values fig, axes = setup_subplots(1, 2, width=14, height=5) # (a) Win rate comparison: standard CES vs FA-CES (BIC) ax = axes[0] x = np.arange(len(fa_digits)) bw = 0.35 ces_bic = klem_summary_fig['pct_ces_wins_bic'].values fa_bic = fa_summary['pct_fa_wins_bic'].values ax.bar(x - bw / 2, ces_bic, bw, label='Standard CES (6 params)', color='steelblue', edgecolor='k') ax.bar(x + bw / 2, fa_bic, bw, label='FA-CES (9 params)', color='darkorange', edgecolor='k') ax.axhline(50, color='k', linestyle='--', alpha=0.5) ax.set_xticks(x) ax.set_xticklabels([f'{d}-digit' for d in fa_digits]) ax.set_xlabel('NAICS level (finer -> coarser)', fontsize=11) ax.set_ylabel('% industries beating translog (BIC)', fontsize=11) ax.set_title('(a) BIC win rate vs translog (16 params)', fontsize=11) ax.set_ylim(0, 100) ax.legend(fontsize=9) # (b) Median delta-BIC comparison ax = axes[1] ax.plot(fa_digits, klem_summary_fig['dbic_median'].values, 's-', c='steelblue', markersize=8, linewidth=2, label='Standard CES') ax.plot(fa_digits, fa_summary['dbic_fa_median'].values, 'o-', c='darkorange', markersize=8, linewidth=2, label='FA-CES') ax.axhline(0, color='k', linestyle='--', alpha=0.5) ax.set_xlabel('NAICS digits', fontsize=11) ax.set_ylabel('Median delta-BIC (translog - model)', fontsize=11) ax.set_title('(b) BIC gap: CES family vs translog', fontsize=11) ax.invert_xaxis() ax.set_xticks(fa_digits) ax.legend(fontsize=9) fig.suptitle('FA-CES Crossover: Does Factor Augmentation Close the Gap?', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'convergence_horse_race_fa_crossover', subfolder='foundations')) # --- Figure 5: Random re-grouping slope --- rand_df = _cache.get('rand_regrouping') if rand_df is not None: fig, ax = setup_figure(width=8, height=6) # Random blocks: median + CI log_k_r = np.log(rand_df['k'].values) log_ir_r = np.log(rand_df['ir_median'].values) ax.scatter(log_k_r, log_ir_r, c='darkorange', edgecolors='k', s=80, zorder=4, label='Random blocks (median)') # CI error bars for _, row in rand_df.iterrows(): if row['ir_lo'] != row['ir_hi']: ax.plot([np.log(row['k']), np.log(row['k'])], [np.log(row['ir_lo']), np.log(row['ir_hi'])], 'darkorange', linewidth=2.5, alpha=0.5) # Random OLS line rand_slope = _cache.get('rand_slope', np.nan) rand_int = _cache.get('rand_intercept', np.nan) xr_all = np.linspace(min(log_k_r.min(), _cache.get('log_k', log_k_r).min()) - 0.2, max(log_k_r.max(), _cache.get('log_k', log_k_r).max()) + 0.2, 50) if np.isfinite(rand_slope): ax.plot(xr_all, rand_int + rand_slope * xr_all, 'darkorange', linewidth=2, label=f'Random OLS: slope={rand_slope:.3f}') # NAICS points (overlay) naics_log_k = _cache.get('log_k') naics_log_ir = _cache.get('log_ir') naics_slope = _cache.get('slope', np.nan) naics_int = _cache.get('intercept', np.nan) if naics_log_k is not None: ax.scatter(naics_log_k, naics_log_ir, c='steelblue', edgecolors='k', s=80, zorder=3, marker='s', label='NAICS hierarchy') if np.isfinite(naics_slope): ax.plot(xr_all, naics_int + naics_slope * xr_all, 'steelblue', linewidth=2, linestyle='--', label=f'NAICS OLS: slope={naics_slope:.3f}') # Theoretical -0.5 line (anchored to random intercept) if np.isfinite(rand_int): ax.plot(xr_all, rand_int + (-0.5) * xr_all, 'k:', linewidth=1.5, alpha=0.5, label='Theory (i.i.d.): slope=-0.50') ax.set_xlabel('log(block size k)', fontsize=11) ax.set_ylabel('log(median interaction ratio)', fontsize=11) ax.set_title('Random vs NAICS Block Averaging: Interaction Ratio Decay', fontsize=12) ax.legend(fontsize=9, loc='lower left') paths.extend(save_figure(fig, 'convergence_horse_race_random_slope', subfolder='foundations')) 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']})")

Convergence Horse Race Crossover

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps