#!/usr/bin/env python3 """ Empirical Tests for Potts Crystallization (Paper 3c) ===================================================== Tests first-order regime shift (fold bifurcation, hysteresis) using Census Business Dynamics Statistics: 1. Fold threshold crossings — recessions push more sectors past the collapse fold (firm growth < -2%). Mann-Whitney U test. Tests Theorem 3c.2 (fold bifurcation). 2. Fold synchronization — fraction of sectors with negative firm growth jumps during recessions. Mann-Whitney U test. Tests Proposition 3.5 (synchronized fold crossing). 3. Entry rate hysteresis — entry rate recoveries are slower than entry rate declines. Tests asymmetric adjustment from the proved hysteresis theorem (Theorem 3c.3: cost_up < cost_down). Mann-Whitney U test on |negative Δentry| vs positive Δentry. 4. Job destruction asymmetry — job destruction rate spikes (positive changes) are larger than normalizations (negative changes). Second facet of hysteresis via job flow dynamics. Mann-Whitney U test. 5. Herfindahl structural breaks — CUSUM test for breaks in sector concentration. Data: census_bds_naics.csv (19 NAICS sectors, 1978-2023) Uses firms, estabs_entry_rate, estabs_exit_rate, job_destruction_rate columns. Lean theorems: - fold_bifurcation (CESProofs/Dynamics/RegimeShift.lean) — schematic - hysteresis_width_positive — real proof - hysteresis_wider_for_complements — real proof - irreversibility_of_adoption — real proof Review Notes (2026-03-12): Rewritten from v1 which used only firm counts. V1 had bugs: (1) RECESSION_YEARS missing 1980-82, 1990-91; (2) transition clustering 2pp threshold produced stat=Infinity; (3) entropy jumps p=1.0 (secular trend, not regime shifts). V2 adds Mann-Whitney p-values to fold tests, replaces uninformative entropy/clustering tests with hysteresis tests that directly test the proved Theorem 3c.3, and uses BDS entry/exit/job flow data which v1 completely ignored. Theory match: CORRECT for subtests 1-4 (fold + hysteresis theorems). I/O audit: CORRECT — firm counts, entry/exit rates, job destruction rates are all input-side measures. Data-signal alignment: STRONG (subtests 1-2), ADEQUATE (3-4), ADEQUATE (5). Code quality: CLEAN after rewrite. Confidence: MEDIUM — 4 of 5 subtests have formal p-values. """ import numpy as np 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.plotting import setup_figure, setup_subplots, save_figure from common.stats import log_result RECESSION_YEARS = {1980, 1981, 1982, 1990, 1991, 2001, 2008, 2009, 2020} _ALL_SUBTESTS = [ 'potts_crystallization_fold_crossings', 'potts_crystallization_fold_synchronization', 'potts_crystallization_entry_hysteresis', 'potts_crystallization_destruction_asymmetry', 'potts_crystallization_hhi_breaks', ] _cache = {} def _load_bds(): """Load BDS data as DataFrame.""" try: return load_compiled('census_bds_naics') except FileNotFoundError: return None def _firm_growth_by_year(df): """Compute year-over-year firm growth rates by sector. Returns dict: {year: np.array of growth rates across sectors}. """ years = sorted(df['year'].unique()) sectors = sorted(df['naics'].unique()) result = {} for i in range(1, len(years)): yr, yr_prev = years[i], years[i - 1] cur = df[df['year'] == yr].set_index('naics') prev = df[df['year'] == yr_prev].set_index('naics') rates = [] for sec in sectors: if sec in cur.index and sec in prev.index: f_c = float(cur.loc[sec, 'firms']) f_p = float(prev.loc[sec, 'firms']) if f_p > 0: rates.append((f_c - f_p) / f_p) if len(rates) >= 5: result[yr] = np.array(rates) return result # ── Test 1: Fold Threshold Crossings (Mann-Whitney) ── def _test_fold_crossings(df): """Sectors with firm growth < -2%: recession vs normal years. Mann-Whitney U test, one-sided.""" growth = _firm_growth_by_year(df) FOLD_THRESHOLD = -0.02 rec_counts, nor_counts = [], [] for yr, rates in growth.items(): n_below = int(np.sum(rates < FOLD_THRESHOLD)) if yr in RECESSION_YEARS: rec_counts.append(n_below) else: nor_counts.append(n_below) if len(rec_counts) < 3 or len(nor_counts) < 3: return log_result('potts_crystallization_fold_crossings', np.nan, np.nan, 'SKIPPED', description="Insufficient recession/normal data.") u_stat, p_mw = stats.mannwhitneyu(rec_counts, nor_counts, alternative='greater') mean_rec = float(np.mean(rec_counts)) mean_nor = float(np.mean(nor_counts)) ratio = mean_rec / mean_nor if mean_nor > 0 else float('inf') if p_mw < 0.05: verd = 'CONSISTENT' elif p_mw < 0.10: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['fold_crossings'] = { 'recession_counts': rec_counts, 'normal_counts': nor_counts, 'ratio': ratio, 'p': p_mw, } return log_result('potts_crystallization_fold_crossings', ratio, p_mw, verd, description=( f"Tests Theorem 3c.2 fold bifurcation: recessions should push " f"more sectors past the collapse fold. Mann-Whitney U test " f"(one-sided): count of sectors with firm growth < -2% in " f"recession years ({len(rec_counts)}) vs normal years " f"({len(nor_counts)}). Data: Census BDS 19 NAICS sectors, " f"1978-2023. Recession mean: {mean_rec:.2f}; normal: " f"{mean_nor:.2f} (ratio {ratio:.2f}). p = {p_mw:.4f}." ), details={'mean_recession': round(mean_rec, 2), 'mean_normal': round(mean_nor, 2), 'ratio': round(ratio, 2), 'U': float(u_stat), 'n_recession_years': len(rec_counts), 'n_normal_years': len(nor_counts)}) # ── Test 2: Fold Synchronization (Mann-Whitney) ── def _test_fold_synchronization(df): """Fraction of sectors declining: recession vs normal years. Mann-Whitney U test, one-sided.""" growth = _firm_growth_by_year(df) rec_frac, nor_frac = [], [] for yr, rates in growth.items(): frac = float(np.mean(rates < 0)) if yr in RECESSION_YEARS: rec_frac.append(frac) else: nor_frac.append(frac) if len(rec_frac) < 3 or len(nor_frac) < 3: return log_result('potts_crystallization_fold_synchronization', np.nan, np.nan, 'SKIPPED', description="Insufficient recession/normal data.") u_stat, p_mw = stats.mannwhitneyu(rec_frac, nor_frac, alternative='greater') mean_rec = float(np.mean(rec_frac)) mean_nor = float(np.mean(nor_frac)) ratio = mean_rec / mean_nor if mean_nor > 0 else float('inf') if p_mw < 0.05: verd = 'CONSISTENT' elif p_mw < 0.10: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['synchronization'] = { 'recession_frac': rec_frac, 'normal_frac': nor_frac, 'ratio': ratio, 'p': p_mw, } return log_result('potts_crystallization_fold_synchronization', ratio, p_mw, verd, description=( f"Tests Proposition 3.5: fraction of sectors with negative " f"firm growth should jump during recessions via synchronized " f"fold crossing. Mann-Whitney U test (one-sided): fraction " f"declining in recession ({len(rec_frac)} yr) vs normal " f"({len(nor_frac)} yr). Data: Census BDS 19 sectors, " f"1978-2023. Recession: {mean_rec:.1%}; normal: " f"{mean_nor:.1%} (ratio {ratio:.2f}). p = {p_mw:.4f}." ), details={'recession_mean': round(mean_rec, 4), 'normal_mean': round(mean_nor, 4), 'ratio': round(ratio, 2), 'U': float(u_stat), 'p': round(p_mw, 4)}) # ── Test 3: Entry Rate Hysteresis ── def _test_entry_hysteresis(df): """Tests Theorem 3c.3 (hysteresis): entry rate drops should be larger than entry rate recoveries, because entry requires coordination (cost_up < cost_down → entry is harder than exit). Method: Compare |negative Δentry_rate| vs positive Δentry_rate using Mann-Whitney U test. Hysteresis predicts negative changes are larger in magnitude than positive changes.""" years = sorted(df['year'].unique()) # Compute economy-wide mean entry rate per year entry_rates = np.array([ float(df[df['year'] == yr]['estabs_entry_rate'].mean()) for yr in years ]) d_entry = np.diff(entry_rates) neg = d_entry[d_entry < 0] pos = d_entry[d_entry > 0] if len(neg) < 3 or len(pos) < 3: return log_result('potts_crystallization_entry_hysteresis', np.nan, np.nan, 'SKIPPED', description="Insufficient entry rate change data.") u_stat, p_mw = stats.mannwhitneyu(np.abs(neg), pos, alternative='greater') mean_neg = float(np.mean(np.abs(neg))) mean_pos = float(np.mean(pos)) asym_ratio = mean_neg / mean_pos if mean_pos > 0 else float('inf') if p_mw < 0.05: verd = 'CONSISTENT' elif p_mw < 0.10: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['entry_hysteresis'] = { 'years': years[1:], 'd_entry': d_entry, 'neg_magnitudes': np.abs(neg), 'pos_magnitudes': pos, 'p': p_mw, 'asym_ratio': asym_ratio, } return log_result('potts_crystallization_entry_hysteresis', asym_ratio, p_mw, verd, description=( f"Tests Theorem 3c.3 hysteresis (proved): entry rate " f"should decline faster than it recovers because entry " f"requires coordination (cost_up < cost_down). Mann-Whitney " f"U test: |negative Δentry_rate| ({len(neg)} episodes) vs " f"positive Δentry_rate ({len(pos)} episodes). Data: Census " f"BDS establishment entry rates across 19 sectors, 1978-2023. " f"Mean |decline|: {mean_neg:.3f}; mean recovery: {mean_pos:.3f} " f"(ratio {asym_ratio:.2f}). p = {p_mw:.4f}." ), details={'mean_neg_magnitude': round(mean_neg, 3), 'mean_pos_magnitude': round(mean_pos, 3), 'asymmetry_ratio': round(asym_ratio, 2), 'U': float(u_stat), 'n_declines': len(neg), 'n_recoveries': len(pos)}) # ── Test 4: Job Destruction Asymmetry ── def _test_destruction_asymmetry(df): """Tests hysteresis via job flows: destruction rate spikes should be larger than destruction rate normalizations. Method: Compare positive Δjob_destruction_rate (crisis onset) vs |negative Δjob_destruction_rate| (recovery) using Mann-Whitney U.""" years = sorted(df['year'].unique()) jd_rates = np.array([ float(df[df['year'] == yr]['job_destruction_rate'].mean()) for yr in years ]) d_jd = np.diff(jd_rates) pos = d_jd[d_jd > 0] # spikes (crisis onset) neg = d_jd[d_jd < 0] # normalizations (recovery) if len(pos) < 3 or len(neg) < 3: return log_result('potts_crystallization_destruction_asymmetry', np.nan, np.nan, 'SKIPPED', description="Insufficient job destruction data.") u_stat, p_mw = stats.mannwhitneyu(pos, np.abs(neg), alternative='greater') mean_spike = float(np.mean(pos)) mean_norm = float(np.mean(np.abs(neg))) asym_ratio = mean_spike / mean_norm if mean_norm > 0 else float('inf') if p_mw < 0.05: verd = 'CONSISTENT' elif p_mw < 0.10: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['destruction_asymmetry'] = { 'years': years[1:], 'd_jd': d_jd, 'spike_magnitudes': pos, 'norm_magnitudes': np.abs(neg), 'p': p_mw, 'asym_ratio': asym_ratio, } return log_result('potts_crystallization_destruction_asymmetry', asym_ratio, p_mw, verd, description=( f"Tests hysteresis via job flows: destruction rate spikes " f"(positive Δjd) should exceed normalizations (|negative Δjd|). " f"Mann-Whitney U test (one-sided). Data: Census BDS job " f"destruction rates across 19 sectors, 1978-2023. " f"Mean spike: {mean_spike:.3f}; mean normalization: " f"{mean_norm:.3f} (ratio {asym_ratio:.2f}). p = {p_mw:.4f}." ), details={'mean_spike': round(mean_spike, 3), 'mean_normalization': round(mean_norm, 3), 'asymmetry_ratio': round(asym_ratio, 2), 'U': float(u_stat), 'n_spikes': len(pos), 'n_normalizations': len(neg)}) # ── Test 5: Herfindahl Structural Breaks ── def _test_hhi_breaks(df): """CUSUM test for structural breaks in sector concentration. Prediction: |CUSUM| > 1.96.""" years = sorted(df['year'].unique()) sectors = sorted(df['naics'].unique()) hhi_series = [] for yr in years: sub = df[df['year'] == yr].set_index('naics') firms = np.array([float(sub.loc[s, 'firms']) if s in sub.index else 0 for s in sectors]) total = np.sum(firms) shares = firms / total if total > 0 else np.zeros(len(sectors)) hhi_series.append(float(np.sum(shares ** 2))) hhi_vals = np.array(hhi_series) mean_hhi = np.mean(hhi_vals) std_hhi = np.std(hhi_vals) cusum = np.cumsum(hhi_vals - mean_hhi) cusum_norm = cusum / (std_hhi * np.sqrt(len(hhi_vals))) max_idx = int(np.argmax(np.abs(cusum_norm))) max_cusum = float(abs(cusum_norm[max_idx])) break_year = int(years[max_idx]) verd = 'CONSISTENT' if max_cusum > 1.96 else 'INCONSISTENT' _cache['hhi'] = { 'years': [int(y) for y in years], 'hhi': hhi_vals, 'cusum': cusum_norm, 'break_year': break_year, 'max_cusum': max_cusum, } return log_result('potts_crystallization_hhi_breaks', max_cusum, np.nan, verd, description=( f"Tests for structural breaks in sector concentration via " f"CUSUM test on Herfindahl index of NAICS firm-count shares. " f"|CUSUM| > 1.96 indicates significant break. Data: Census " f"BDS 19 sectors, 1978-2023. Break at {break_year} " f"(CUSUM={max_cusum:.2f}). Note: may reflect secular " f"composition shifts (services growth) rather than " f"Potts crystallization." ), details={'break_year': break_year, 'max_cusum': round(max_cusum, 3)}) # ── Main ────────────────────────────────────────────────────────── def main(): """Run all Potts crystallization tests. Returns list of result dicts.""" df = _load_bds() if df is None: return [log_result('potts_crystallization', np.nan, np.nan, 'SKIPPED', description="Census BDS NAICS data file not found.", details={'reason': 'BDS data not found'})] results = [ _test_fold_crossings(df), _test_fold_synchronization(df), _test_entry_hysteresis(df), _test_destruction_asymmetry(df), _test_hhi_breaks(df), ] _cache['results'] = results return results def figures(): """Generate figures for Potts crystallization subtests.""" if not _cache: main() paths = [] # --- Figure 1: Fold crossings (box plot with p-value) --- d = _cache.get('fold_crossings') if d is not None: fig, ax = setup_figure() data = [d['normal_counts'], d['recession_counts']] bp = ax.boxplot(data, labels=['Normal Years', 'Recession Years'], patch_artist=True) bp['boxes'][0].set_facecolor('lightblue') bp['boxes'][1].set_facecolor('salmon') ax.set_ylabel('Sectors with Firm Growth < -2%') ax.set_title('Fold Threshold Crossings: Recession vs Normal') ax.text(0.98, 0.95, f'MW p = {d["p"]:.4f}\nratio = {d["ratio"]:.2f}', transform=ax.transAxes, ha='right', va='top', fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) paths.extend(save_figure(fig, 'potts_crystallization_fold_crossings', subfolder='curvature_roles')) # --- Figure 2: Fold synchronization (box plot with p-value) --- d = _cache.get('synchronization') if d is not None: fig, ax = setup_figure() data = [d['normal_frac'], d['recession_frac']] bp = ax.boxplot(data, labels=['Normal Years', 'Recession Years'], patch_artist=True) bp['boxes'][0].set_facecolor('lightblue') bp['boxes'][1].set_facecolor('salmon') ax.set_ylabel('Fraction of Sectors Declining') ax.set_title('Fold Synchronization: Recession vs Normal') ax.text(0.98, 0.95, f'MW p = {d["p"]:.4f}\nratio = {d["ratio"]:.2f}', transform=ax.transAxes, ha='right', va='top', fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) paths.extend(save_figure(fig, 'potts_crystallization_synchronization', subfolder='curvature_roles')) # --- Figure 3: Entry rate hysteresis --- d = _cache.get('entry_hysteresis') if d is not None: fig, axes = setup_subplots(nrows=2, ncols=1, width=10, height=8) # Top: entry rate changes time series colors = ['salmon' if v < 0 else 'steelblue' for v in d['d_entry']] axes[0].bar(d['years'], d['d_entry'], color=colors, alpha=0.7) axes[0].axhline(0, color='gray', linewidth=0.8) axes[0].set_ylabel('Δ Entry Rate (pp)') axes[0].set_title('Establishment Entry Rate Changes') # Bottom: comparison of magnitudes bp = axes[1].boxplot([d['pos_magnitudes'], d['neg_magnitudes']], labels=['Recoveries (+Δ)', 'Declines (|−Δ|)'], patch_artist=True) bp['boxes'][0].set_facecolor('lightblue') bp['boxes'][1].set_facecolor('salmon') axes[1].set_ylabel('Magnitude (pp)') axes[1].set_title(f'Hysteresis: Declines > Recoveries (p={d["p"]:.4f})') fig.tight_layout() paths.extend(save_figure(fig, 'potts_crystallization_entry_hysteresis', subfolder='curvature_roles')) # --- Figure 4: Job destruction asymmetry --- d = _cache.get('destruction_asymmetry') if d is not None: fig, axes = setup_subplots(nrows=2, ncols=1, width=10, height=8) colors = ['salmon' if v > 0 else 'steelblue' for v in d['d_jd']] axes[0].bar(d['years'], d['d_jd'], color=colors, alpha=0.7) axes[0].axhline(0, color='gray', linewidth=0.8) axes[0].set_ylabel('Δ Job Destruction Rate (pp)') axes[0].set_title('Job Destruction Rate Changes') bp = axes[1].boxplot([d['spike_magnitudes'], d['norm_magnitudes']], labels=['Spikes (+Δjd)', 'Normalizations (|−Δjd|)'], patch_artist=True) bp['boxes'][0].set_facecolor('salmon') bp['boxes'][1].set_facecolor('lightblue') axes[1].set_ylabel('Magnitude (pp)') axes[1].set_title(f'Hysteresis: Spikes > Normalizations (p={d["p"]:.4f})') fig.tight_layout() paths.extend(save_figure(fig, 'potts_crystallization_destruction_asymmetry', subfolder='curvature_roles')) # --- Figure 5: HHI CUSUM --- d = _cache.get('hhi') if d is not None: fig, axes = setup_subplots(nrows=2, ncols=1, width=10, height=8) axes[0].plot(d['years'], d['hhi'], 'o-', color='steelblue', markersize=4) axes[0].set_ylabel('Herfindahl Index') axes[0].set_title('NAICS Sector Concentration') axes[1].plot(d['years'], d['cusum'], 'o-', color='darkgreen', markersize=4) axes[1].axhline(1.96, color='red', linestyle='--', label='1.96') axes[1].axhline(-1.96, color='red', linestyle='--', label='-1.96') axes[1].axvline(d['break_year'], color='orange', linewidth=2, label=f'Break: {d["break_year"]}') axes[1].set_xlabel('Year') axes[1].set_ylabel('Normalized CUSUM') axes[1].set_title(f'CUSUM Test (max={d["max_cusum"]:.2f})') axes[1].legend(fontsize=8) fig.tight_layout() paths.extend(save_figure(fig, 'potts_crystallization_hhi', subfolder='curvature_roles')) 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']})")

max cusum	2.7541
statistic	2.7541
break year	1992

Potts Crystallization Hhi Breaks

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps