#!/usr/bin/env python3 """ Test: Leontief Spectral Stability — Paper 5 (Empirical Tests), Section: Leontief Hierarchy Discovery Tests hierarchy depth and spectral spacing using the off-diagonal structure of the Leontief inverse matrix L = (I-A)^{-1}. Tests: 1. Off-diagonal SVD: r* = sigma_1/sigma_2 in [1.5, 4.0] 2. Hierarchy depth: N_eff (80% energy) in [3, 8] 3. Rolling stability: CV of N_eff < 0.3 (skipped — no rolling data) 4. EMD cross-validation: |N_eff_leontief - N_eff_EMD| <= 2 5. Random matrix null: economy more concentrated than random Data: BEA Leontief inverse (71x71), EMD results. Theory: Paper 5, Result (res:leontief_levels): "the effective rank of L_off (80% cumulative energy) is N_eff = 5, identical to the EMD estimate." The hierarchy discoverable from static I-O structure should match the hierarchy discoverable from dynamic timescale decomposition (EMD). Lean: CESProofs/Hierarchy/ModuliSpace.lean — hierarchy_depth_finite; CESProofs/Hierarchy/SpectralHierarchy.lean — spectralGapRatio, timescale_ratio_eq_rate_ratio. Review Notes (2026-03-11): Theory match: CORRECT — all subtests directly implement Paper 5 claims I/O audit: CORRECT — BEA Leontief inverse is the total requirements matrix; EMD N_eff from INDPRO EEMD decomposition Data-signal: STRONG — both measurements are direct (no proxies) Code quality: MINOR ISSUES — Test 3 skipped (no rolling_neff data); random null p=0.0 is technically p<0.005 (200 permutations) Confidence: HIGH — all 4 running subtests CONSISTENT with strong margins """ import warnings warnings.filterwarnings("ignore") import numpy as np import pandas as pd import csv import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.config import COMPILED_DIR, DATA_DIR from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure def _load_leontief_matrix(): """Load the 71x71 Leontief inverse matrix.""" path = COMPILED_DIR / 'bea_leontief_full.csv' if not path.exists(): raise FileNotFoundError(f"Leontief matrix not found at {path}") with open(path) as f: reader = csv.reader(f) header = next(reader) rows = [] for row in reader: vals = [] for v in row[1:]: try: vals.append(float(v)) except ValueError: vals.append(0.0) rows.append(vals) return np.array(rows), header[1:] def _offdiag_svd(L): L_offdiag = L - np.diag(np.diag(L)) sv = np.linalg.svd(L_offdiag, compute_uv=False) return L_offdiag, np.sort(sv)[::-1] def _load_csv_dicts(fname): path = COMPILED_DIR / fname if not path.exists(): return None with open(path) as f: reader = csv.DictReader(f) rows = [] for row in reader: d = {} for k, v in row.items(): try: d[k] = float(v.replace(",", "")) except (ValueError, AttributeError): d[k] = v rows.append(d) return rows _cache = {} def main(): results = [] try: L, codes = _load_leontief_matrix() except FileNotFoundError as e: results.append(log_result( "leontief_spectral", np.nan, np.nan, "AMBIGUOUS", details={"reason": str(e)}, description="Skipped: BEA 71x71 Leontief inverse matrix not available for spectral analysis." )) return results L_offdiag, sv = _offdiag_svd(L) _cache['sv'] = sv # --- Test 1: Off-diagonal SVD spectrum --- r_star = sv[0] / sv[1] if sv[1] > 0 else float('inf') if 1.5 <= r_star <= 4.0: verd = "CONSISTENT" elif r_star > 4.0: verd = "AMBIGUOUS" else: verd = "INCONSISTENT" results.append(log_result( "leontief_spectral_offdiag_svd", stat=float(r_star), p=np.nan, verd=verd, details={"sigma_1": float(sv[0]), "sigma_2": float(sv[1])}, description="Tests whether the off-diagonal SVD spectrum of the BEA Leontief inverse shows a " "leading singular value ratio r* = sigma_1/sigma_2 in [1.5, 4.0], consistent with " "hierarchical layer structure rather than flat random coupling." )) # --- Test 2: Hierarchy depth N_eff --- energy_total = np.sum(sv ** 2) cum_energy = np.cumsum(sv ** 2) / energy_total n_eff = int(np.searchsorted(cum_energy, 0.80) + 1) _cache['cum_energy'] = cum_energy _cache['n_eff'] = n_eff verd = "CONSISTENT" if 3 <= n_eff <= 8 else "INCONSISTENT" results.append(log_result( "leontief_spectral_hierarchy_depth", stat=float(n_eff), p=np.nan, verd=verd, description="Counts the effective hierarchy depth N_eff as the number of singular values needed " "to capture 80% of off-diagonal energy in the Leontief inverse. " "Consistent if 3 <= N_eff <= 8, matching the CES-predicted hierarchy depth." )) # --- Test 3: Rolling stability --- rolling_data = _load_csv_dicts("sector_layers/rolling_neff.csv") if rolling_data is not None: sectors = {} for r in rolling_data: sec = r.get("sector_id") or r.get("name", "unknown") neff_val = r.get("n_eff") if neff_val is not None and not (isinstance(neff_val, float) and np.isnan(neff_val)): sectors.setdefault(sec, []).append(float(neff_val)) cvs = [] for sec, vals in sectors.items(): if len(vals) >= 3: mean_v = np.mean(vals) if mean_v > 0: cvs.append(np.std(vals) / mean_v) if cvs: overall_cv = np.mean(cvs) if overall_cv < 0.3: verd = "CONSISTENT" elif overall_cv < 0.5: verd = "AMBIGUOUS" else: verd = "INCONSISTENT" results.append(log_result( "leontief_spectral_rolling_stability", stat=float(overall_cv), p=np.nan, verd=verd, description="Tests whether the effective hierarchy depth N_eff is stable over time by computing " "the coefficient of variation (CV) of rolling N_eff across sectors. " "Consistent if CV < 0.3, confirming structural rather than transient hierarchy." )) # --- Test 4: EMD cross-validation --- emd_data = _load_csv_dicts("emd_hierarchy_results.csv") if emd_data is not None: emd_neff = None for r in emd_data: if "Neff" in str(r.get("test", "")): emd_neff = float(r["value"]) break if emd_neff is not None: diff = abs(n_eff - emd_neff) verd = "CONSISTENT" if diff <= 2 else "INCONSISTENT" results.append(log_result( "leontief_spectral_emd_crossval", stat=float(diff), p=np.nan, verd=verd, details={"leontief_neff": n_eff, "emd_neff": float(emd_neff)}, description="Cross-validates Leontief-derived N_eff against EMD-derived N_eff from time-series " "decomposition. Two independent methods should agree on hierarchy depth. " "Consistent if |N_eff_leontief - N_eff_EMD| <= 2." )) # --- Test 5: Random matrix null --- np.random.seed(42) n = L.shape[0] offdiag_mask = L_offdiag > 1e-10 density = np.mean(offdiag_mask) offdiag_vals = L_offdiag[offdiag_mask] val_mean = np.mean(offdiag_vals) val_std = np.std(offdiag_vals) actual_energy_5 = float(cum_energy[min(4, len(cum_energy)-1)]) n_random = 200 random_energy_5 = [] for _ in range(n_random): R = np.zeros((n, n)) mask = np.random.random((n, n)) < density np.fill_diagonal(mask, False) R[mask] = np.abs(np.random.normal(val_mean, val_std, size=np.sum(mask))) sv_r = np.linalg.svd(R, compute_uv=False) sv_r = np.sort(sv_r)[::-1] cum_e = np.cumsum(sv_r ** 2) / np.sum(sv_r ** 2) random_energy_5.append(float(cum_e[min(4, len(cum_e)-1)])) random_energy_5 = np.array(random_energy_5) _cache['random_energy_5'] = random_energy_5 _cache['actual_energy_5'] = actual_energy_5 p_value = float(np.mean(random_energy_5 >= actual_energy_5)) verd = "CONSISTENT" if p_value < 0.05 else "INCONSISTENT" results.append(log_result( "leontief_spectral_random_null", stat=float(actual_energy_5), p=p_value, verd=verd, details={"random_mean": float(np.mean(random_energy_5))}, description="Random matrix null test: checks whether the top-5 singular values of the Leontief " "off-diagonal capture more energy than expected from 200 random matrices matched on " "density and value distribution. Consistent if p < 0.05 (economy is more concentrated)." )) _cache['results'] = results return results def figures(): """Generate figures for each subtest.""" if not _cache: main() paths = [] results = _cache.get('results', []) sv = _cache.get('sv') cum_energy = _cache.get('cum_energy') n_eff = _cache.get('n_eff') # Figure 1: SVD spectrum (covers offdiag_svd + hierarchy_depth) if sv is not None and cum_energy is not None: fig, (ax1, ax2) = setup_subplots(1, 2, width=14, height=5) # Left: singular value bar chart n_show = min(15, len(sv)) ax1.bar(range(1, n_show + 1), sv[:n_show], color='steelblue', alpha=0.8) r_star = sv[0] / sv[1] if sv[1] > 0 else float('inf') ax1.set_xlabel('Singular value index') ax1.set_ylabel('Singular value') ax1.set_title(f'Off-diagonal SVD Spectrum (BEA Leontief 71x71)\n' f'r* = σ₁/σ₂ = {r_star:.2f}') ax1.axhline(sv[0], color='red', linestyle='--', alpha=0.3) ax1.axhline(sv[1], color='orange', linestyle='--', alpha=0.3) ax1.annotate(f'σ₁ = {sv[0]:.3f}', (1, sv[0]), xytext=(3, sv[0] + 0.05), fontsize=8, color='red') ax1.annotate(f'σ₂ = {sv[1]:.3f}', (2, sv[1]), xytext=(3.5, sv[1] + 0.05), fontsize=8, color='orange') # Right: cumulative energy curve n_show_cum = min(20, len(cum_energy)) ax2.plot(range(1, n_show_cum + 1), cum_energy[:n_show_cum], 'o-', color='steelblue') ax2.axhline(0.80, color='red', linestyle='--', label='80% threshold') ax2.axvline(n_eff, color='green', linestyle='--', alpha=0.5, label=f'N_eff = {n_eff}') ax2.set_xlabel('Number of singular values') ax2.set_ylabel('Cumulative energy fraction') ax2.set_title(f'Hierarchy Depth: N_eff = {n_eff}\n(80% cumulative energy)') ax2.legend(fontsize=9) ax2.set_ylim(0, 1.05) ax2.grid(True, alpha=0.3) fig.tight_layout() paths.extend(save_figure(fig, 'leontief_spectral_svd_depth', subfolder='hierarchy')) # Figure 2: EMD cross-validation comparison r_emd = next((r for r in results if r['name'] == 'leontief_spectral_emd_crossval'), None) if r_emd and r_emd['verdict'] != 'SKIPPED': leon_neff = r_emd.get('leontief_neff', n_eff) emd_neff = r_emd.get('emd_neff', 0) fig, ax = setup_figure(8, 5) bars = ax.bar(['Leontief SVD\n(BEA I-O matrix)', 'EMD\n(INDPRO time series)'], [leon_neff, emd_neff], color=['steelblue', 'darkorange'], width=0.4) for bar, val in zip(bars, [leon_neff, emd_neff]): ax.text(bar.get_x() + bar.get_width()/2, val + 0.1, f'{val:.0f}', ha='center', va='bottom', fontweight='bold', fontsize=14) ax.set_ylabel('N_eff (hierarchy depth)') diff = abs(leon_neff - emd_neff) ax.set_title(f'EMD Cross-Validation: Two Methods Agree\n' f'|{leon_neff} - {emd_neff:.0f}| = {diff:.0f} ≤ 2') ax.set_ylim(0, max(leon_neff, emd_neff) + 2) ax.axhspan(leon_neff - 2, leon_neff + 2, color='green', alpha=0.1, label='±2 agreement band') ax.legend(fontsize=9) paths.extend(save_figure(fig, 'leontief_spectral_emd_crossval', subfolder='hierarchy')) # Figure 3: Random matrix null random_e5 = _cache.get('random_energy_5') actual_e5 = _cache.get('actual_energy_5') if random_e5 is not None and actual_e5 is not None: fig, ax = setup_figure(8, 5) ax.hist(random_e5, bins=25, color='gray', alpha=0.7, edgecolor='black', label='Random matrices (N=200)') ax.axvline(actual_e5, color='red', linewidth=2, linestyle='--', label=f'Actual economy = {actual_e5:.3f}') p_val = float(np.mean(random_e5 >= actual_e5)) ax.set_xlabel('Top-5 SV cumulative energy fraction') ax.set_ylabel('Count') ax.set_title(f'Random Matrix Null Test\n' f'Economy ({actual_e5:.3f}) vs Random (mean {np.mean(random_e5):.3f}), ' f'p < {max(p_val, 1/201):.3f}') ax.legend(fontsize=9) paths.extend(save_figure(fig, 'leontief_spectral_random_null', subfolder='hierarchy')) return paths if __name__ == '__main__': res = main() for r in res: print(f" {r['name']:50s} {r['verdict']:15s} stat={r['statistic']}")

p value	0
statistic	0.8181
random mean	0.6904

Leontief Spectral Random Null

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps