#!/usr/bin/env python3 """ Compound Symmetry Test: CES Covariance Eigenstructure ===================================================== Paper 3 (Dynamics on the CES Potential Landscape), Section "Permutation Invariance" Paper 3b (Dynamics under General Weights), Section "Conservation Laws under General Weights" Lean theorems (equal-weight / compound symmetry): compound_symmetry_cov, compound_symmetry_cov_perp, compound_symmetry_decomp, compound_symmetry_market_eigenvalue, compound_symmetry_idio_eigenvalue, compound_symmetry_trace (all in CESProofs/Dynamics/ConservationLaws.lean, proved, 0 sorry) Lean theorems (general-weight extension): gen_hessian_equal_marginals, genCesHessianQF, secular_equation, generalCurvatureK, K_reduction_equal_weights (CESProofs/Foundations/GeneralWeights.lean) Mathematical prediction: EQUAL WEIGHTS: Under permutation-symmetric CES, the equilibrium covariance matrix takes compound symmetry form: Sigma = (s^2 - g)*I + g*1*1^T Eigenvalues: 1 market mode s^2 + (J-1)*g on 1, J-1 equal idiosyncratic modes s^2 - g on 1-perp. All off-diagonal correlations equal. GENERAL WEIGHTS (Paper 1 Sec 10, Paper 3b): When weights a_j differ, permutation invariance is broken. The Hessian becomes a rank-1 perturbation of diag(1/p_j) where p_j = a_j^sigma are cost-minimizing shares: H = alpha * 11^T - beta * diag(1/p_j) Eigenvalues on 1-perp are determined by the secular equation: sum_j 1/(w_j - mu) = 0, w_j = a_j^{-sigma} interlacing the diagonal entries. The key structural feature — 1 dominant eigenvalue (rank-1 perturbation) — is preserved under all weight vectors. Equicorrelation is broken: correlations become weight-dependent. Sub-tests: 6a: Eigenvalue dominance — lambda_1 / sum(lambda) >> 1/J (rank-1 structure) 6b: Equicorrelation — CV of off-diagonal correlations (equal-weight test) 6c: Implied parameters — g > 0 (positive common covariance) and s^2 > g 6d: Marchenko-Pastur — bulk eigenvalues vs noise (tests equal-weight special case) 6e: Sub-period stability — structure holds in both halves of sample Review Notes: - Theory match: CORRECT — tests both equal-weight (compound symmetry) and structural (rank-1 dominance) predictions - I/O audit: CORRECT — uses FRED IP (output), theory predicts covariance of sector outputs under CES aggregation - Data-signal alignment: STRONG — FRED manufacturing IP is direct measure, J=18 sectors, T=648 months - Code quality: BUGS FOUND (fixed 2026-03-11) — no p-values, single log_result, missing _ALL_SUBTESTS - Confidence: HIGH - Date of review: 2026-03-11 """ import numpy as np import pandas as pd from common.data import load_fred from common.stats import log_result, verdict from common.plotting import setup_figure, setup_subplots, save_figure # Manufacturing IP series (3-digit NAICS subsectors) MFG_SECTORS = [ 'IPG311A2S', 'IPG313A4S', 'IPG315A6S', 'IPG321S', 'IPG322S', 'IPG323S', 'IPG324S', 'IPG325S', 'IPG326S', 'IPG327S', 'IPG331S', 'IPG332S', 'IPG333S', 'IPG334S', 'IPG335S', 'IPG336S', 'IPG337S', 'IPG339S', ] _ALL_SUBTESTS = [ 'compound_symmetry_eigenvalue_dominance', 'compound_symmetry_equicorrelation', 'compound_symmetry_implied_parameters', 'compound_symmetry_marchenko_pastur', 'compound_symmetry_stability', 'compound_symmetry_variance_heterogeneity', 'compound_symmetry_variance_share', 'compound_symmetry_broad_eigenvalue', 'compound_symmetry_broad_variance_share', 'compound_symmetry_correlation_size', 'compound_symmetry_eigenmode_diagonality', 'compound_symmetry_secular_eigenvalues', 'compound_symmetry_mfg_stability', 'compound_symmetry_mfg_eigenmode_diag', ] # FRED IP series → BEA 3-digit NAICS mapping for VA share lookup # BEA splits transportation (336) into 3361MV + 3364OT; we sum them. FRED_TO_BEA = { 'IPG311A2S': ['311FT'], 'IPG313A4S': ['313TT'], 'IPG315A6S': ['315AL'], 'IPG321S': ['321'], 'IPG322S': ['322'], 'IPG323S': ['323'], 'IPG324S': ['324'], 'IPG325S': ['325'], 'IPG326S': ['326'], 'IPG327S': ['327'], 'IPG331S': ['331'], 'IPG332S': ['332'], 'IPG333S': ['333'], 'IPG334S': ['334'], 'IPG335S': ['335'], 'IPG336S': ['3361MV', '3364OT'], # Combined transportation equipment 'IPG337S': ['337'], 'IPG339S': ['339'], } # BEA aggregate codes to exclude (leaves only) BEA_AGGREGATES = { 'GDP', 'PVT', 'PGOOD', 'PSERV', '31G', '31ND', '33DG', '44RT', '48TW', '51', '52', '53', '54', '6', '7', '11', '21', '56', 'FIRE', 'PROF', 'G', 'GF', 'GFE', 'GFG', 'GFGD', 'GFGN', 'GSL', 'GSLE', 'GSLG', 'HS', 'ORE', 'NABI', } def _load_bea_broad_data(): """Load BEA quarterly VA for all leaf industries, return (growth_df, shares_dict). Returns quarterly log-growth rates and average VA shares for ~65 leaf industries. Excludes aggregates, government, imputed sectors, and industries with non-positive VA. """ try: df = pd.read_csv('data/compiled/bea_real_va_quarterly.csv') except FileNotFoundError: return None, None df['DataValue'] = pd.to_numeric(df['DataValue'], errors='coerce') leaves = [c for c in df['Industry'].unique() if c not in BEA_AGGREGATES] leaf_df = df[df['Industry'].isin(leaves)].copy() # Pivot to wide format pivot = leaf_df.pivot_table( index=['Year', 'Quarter'], columns='Industry', values='DataValue', aggfunc='first') pivot = pivot.sort_index() # Drop industries with any non-positive values (can't take log) for col in pivot.columns: if (pivot[col] <= 0).any(): pivot = pivot.drop(columns=[col]) # Drop industries with missing data pivot = pivot.dropna(axis=1) if pivot.shape[1] < 10: return None, None # Compute log-growth log_va = np.log(pivot) growth = log_va.diff().iloc[1:] # Compute average VA shares avg_va = pivot.mean() total_va = avg_va.sum() shares = {ind: avg_va[ind] / total_va for ind in avg_va.index} return growth, shares def _load_bea_va_shares(): """Load BEA quarterly VA data and compute average manufacturing sector shares.""" try: df = pd.read_csv('data/compiled/bea_real_va_quarterly.csv') except FileNotFoundError: return None leaf_codes = set() for bea_list in FRED_TO_BEA.values(): leaf_codes.update(bea_list) mfg = df[df['Industry'].isin(leaf_codes)].copy() mfg['DataValue'] = pd.to_numeric(mfg['DataValue'], errors='coerce') avg_va = mfg.groupby('Industry')['DataValue'].mean() # Map to FRED series shares = {} for fred_id, bea_codes in FRED_TO_BEA.items(): va = sum(avg_va.get(c, 0) for c in bea_codes) shares[fred_id] = va total = sum(shares.values()) if total <= 0: return None return {k: v / total for k, v in shares.items()} def _compute_corr_cv(growth_df): """Compute CV of off-diagonal correlations from growth data.""" cov_matrix = growth_df.cov().values diag = np.sqrt(np.diag(cov_matrix)) corr = cov_matrix / np.outer(diag, diag) mask = np.triu(np.ones_like(corr, dtype=bool), k=1) off_diag = corr[mask] mean_od = np.mean(off_diag) if abs(mean_od) < 1e-10: return float('inf'), off_diag, corr cv = np.std(off_diag, ddof=1) / abs(mean_od) return cv, off_diag, corr def _frobenius_ratio(corr): """Frobenius distance ratio: ||Sigma - Sigma_eq||_F / ||Sigma||_F.""" mask = np.triu(np.ones_like(corr, dtype=bool), k=1) mean_od = np.mean(corr[mask]) # Equicorrelation matrix: 1 on diagonal, mean_od off-diagonal eq_corr = np.full_like(corr, mean_od) np.fill_diagonal(eq_corr, 1.0) diff = corr - eq_corr return np.linalg.norm(diff, 'fro') / np.linalg.norm(corr, 'fro') def _permutation_cv_test(growth_df, observed_cv, n_perm=2000): """ Permutation test for CV of off-diagonal correlations. Null: independently permute each column's time indices. Under no common structure, pairwise correlations are near zero and their CV is very large. """ rng = np.random.RandomState(42) null_cvs = [] vals = growth_df.values.copy() J = vals.shape[1] for _ in range(n_perm): # Independently permute each sector's time series perm_vals = np.column_stack([rng.permutation(vals[:, j]) for j in range(J)]) perm_df = pd.DataFrame(perm_vals) cv_null, _, _ = _compute_corr_cv(perm_df) null_cvs.append(cv_null) null_cvs = np.array(null_cvs) # One-sided: we want observed CV to be LOWER than null p_value = max(np.mean(null_cvs <= observed_cv), 1 / n_perm) return p_value, null_cvs def period_stats(g_sub): """Compute eigenvalue dominance and corr CV for a sub-period.""" cov_sub = g_sub.cov().values evals = np.linalg.eigh(cov_sub)[0] evals = np.sort(evals)[::-1] l1_share = evals[0] / np.sum(evals) cv, _, _ = _compute_corr_cv(g_sub) return l1_share, cv, len(g_sub) _cache = {} def main(): results = [] # Load data ip_data = {} for sid in MFG_SECTORS: try: ip_data[sid] = load_fred(sid) except FileNotFoundError: pass skip_desc = ( "Tests whether manufacturing IP covariance has compound symmetry structure " "(1 large eigenvalue + J-1 equal bulk eigenvalues) as predicted by equal-weight CES. " "Five sub-criteria on eigenvalue spectrum, correlation homogeneity, and sub-period " "stability from 18 FRED IP sectors." ) if len(ip_data) < 5: return [log_result(name, np.nan, np.nan, 'SKIPPED', description=skip_desc) for name in _ALL_SUBTESTS] ip_df = pd.DataFrame(ip_data).dropna() J = ip_df.shape[1] if J < 5: return [log_result(name, np.nan, np.nan, 'SKIPPED', description=skip_desc) for name in _ALL_SUBTESTS] # Compute log-growth log_ip = np.log(ip_df) growth = log_ip.diff().iloc[1:] T = len(growth) # Covariance and eigendecomposition cov_matrix = growth.cov().values eigenvalues = np.linalg.eigh(cov_matrix)[0] eigenvalues = np.sort(eigenvalues)[::-1] total_var = np.sum(eigenvalues) lam1 = eigenvalues[0] bulk = eigenvalues[1:] bulk_mean = np.mean(bulk) # Correlation matrix analysis corr_cv, off_diag_corr, corr_matrix = _compute_corr_cv(growth) frob_ratio = _frobenius_ratio(corr_matrix) # Implied compound symmetry parameters g_implied = (lam1 - bulk_mean) / J s2_implied = bulk_mean + g_implied # Marchenko-Pastur comparison q = J / T sigma2_mp = bulk_mean mp_upper = sigma2_mp * (1 + np.sqrt(q)) ** 2 fraction_below_mp = np.sum(bulk <= mp_upper) / len(bulk) # Sub-period stability midpoint = growth.index[len(growth) // 2] g_early = growth.loc[:midpoint] g_late = growth.loc[midpoint:] l1_early, cv_early, _ = period_stats(g_early) l1_late, cv_late, _ = period_stats(g_late) # ---- Test 6a: Eigenvalue dominance ---- lam1_share = lam1 / total_var # Permutation test: independently permute each column to destroy cross-correlations. # Under H0 (no common factor), lambda1_share ~ 1/J + noise. rng = np.random.RandomState(42) n_perm = 2000 null_shares = [] vals = growth.values.copy() for _ in range(n_perm): perm_vals = np.column_stack([rng.permutation(vals[:, j]) for j in range(J)]) perm_cov = np.cov(perm_vals.T) perm_evals = np.sort(np.linalg.eigh(perm_cov)[0])[::-1] null_shares.append(perm_evals[0] / np.sum(perm_evals)) null_shares = np.array(null_shares) # One-sided: observed should be HIGHER than null (dominant eigenvalue) p_6a = max(np.mean(null_shares >= lam1_share), 1 / n_perm) verd_6a = verdict(p_6a, lam1_share > np.median(null_shares)) results.append(log_result( 'compound_symmetry_eigenvalue_dominance', lam1_share, p_6a, verd_6a, description=( "Tests that the first eigenvalue of the manufacturing IP covariance matrix dominates " "(lambda_1/sum >> 1/J), as predicted by compound symmetry (1 market mode + J-1 equal " "idiosyncratic modes). Permutation test (2000 draws, columns independently permuted) " "provides null distribution under no common structure. " "Data: 18 FRED manufacturing IP series, monthly log-growth rates." ), details={ 'lam1_share': round(lam1_share, 4), 'null_share_median': round(float(np.median(null_shares)), 4), 'null_share_max': round(float(np.max(null_shares)), 4), 'J': J, 'T': T, } )) # ---- Test 6b: Equicorrelation (CV of off-diagonal correlations) ---- p_6b, null_cvs = _permutation_cv_test(growth, corr_cv, n_perm=2000) # sign_correct: CV should be LOW (less variable off-diagonal corr = more equicorrelated) verd_6b = verdict(p_6b, corr_cv < np.median(null_cvs)) results.append(log_result( 'compound_symmetry_equicorrelation', corr_cv, p_6b, verd_6b, description=( "Tests equicorrelation: the CV of off-diagonal correlations should be low under " "compound symmetry (CV=0 for perfect equicorrelation). Permutation test (2000 draws, " "each column independently permuted) provides null distribution. H0: observed CV is " "consistent with no common structure. Data: 18 FRED manufacturing IP series. " "Also reports Frobenius ratio ||Sigma - Sigma_eq||_F / ||Sigma||_F." ), details={ 'corr_cv': round(corr_cv, 4), 'null_cv_median': round(float(np.median(null_cvs)), 4), 'null_cv_percentile': round(float(np.mean(null_cvs <= corr_cv) * 100), 1), 'frobenius_ratio': round(frob_ratio, 4), 'mean_off_diag_corr': round(float(np.mean(off_diag_corr)), 4), 'J': J, 'T': T, } )) # ---- Test 6c: Implied parameters ---- # g > 0 means positive common covariance; s^2 > g means positive idiosyncratic variance # Bootstrap test n_boot = 2000 boot_g = [] boot_s2 = [] for _ in range(n_boot): idx = rng.choice(T, size=T, replace=True) boot_cov = np.cov(growth.values[idx].T) boot_evals = np.sort(np.linalg.eigh(boot_cov)[0])[::-1] b_lam1 = boot_evals[0] b_bulk_mean = np.mean(boot_evals[1:]) b_g = (b_lam1 - b_bulk_mean) / J b_s2 = b_bulk_mean + b_g boot_g.append(b_g) boot_s2.append(b_s2) boot_g = np.array(boot_g) boot_s2 = np.array(boot_s2) # p for g > 0: fraction of bootstrap where g <= 0 p_g = max(np.mean(boot_g <= 0), 1 / n_boot) # p for s^2 > g: fraction where s^2 <= g p_s2g = max(np.mean(boot_s2 <= boot_g), 1 / n_boot) p_6c = max(p_g, p_s2g) # both must hold sign_6c = (g_implied > 0) and (s2_implied > g_implied) verd_6c = verdict(p_6c, sign_6c) results.append(log_result( 'compound_symmetry_implied_parameters', g_implied, p_6c, verd_6c, description=( "Tests that implied compound symmetry parameters are physically valid: g > 0 (positive " "common covariance) and s^2 > g (positive idiosyncratic variance). Under compound " "symmetry Sigma = (s^2-g)I + g*11^T, g = (lambda_1 - bulk_mean)/J and s^2 = bulk_mean " "+ g. Bootstrap (2000 resamples) tests H0: g <= 0 or s^2 <= g. Data: 18 FRED " "manufacturing IP series." ), details={ 'g_implied': round(g_implied, 6), 's2_implied': round(s2_implied, 6), 'g_boot_ci_95': [round(float(np.percentile(boot_g, 2.5)), 6), round(float(np.percentile(boot_g, 97.5)), 6)], 'J': J, 'T': T, } )) # ---- Test 6d: Marchenko-Pastur comparison ---- # Under compound symmetry, bulk eigenvalues should be equal (= s^2 - g). # The MP law gives the expected distribution for a sample covariance with # i.i.d. noise. If bulk eigenvalues fall within MP bounds, they're consistent # with noise around a single common factor. # Test: what fraction of bulk eigenvalues fall below MP upper bound? # Bootstrap for p-value boot_frac = [] for _ in range(n_boot): idx = rng.choice(T, size=T, replace=True) boot_cov = np.cov(growth.values[idx].T) boot_evals = np.sort(np.linalg.eigh(boot_cov)[0])[::-1] b_bulk = boot_evals[1:] b_bulk_mean = np.mean(b_bulk) b_mp_upper = b_bulk_mean * (1 + np.sqrt(J / T)) ** 2 b_frac = np.sum(b_bulk <= b_mp_upper) / len(b_bulk) boot_frac.append(b_frac) boot_frac = np.array(boot_frac) # p-value: fraction of bootstrap where fraction_below < 0.80 p_6d = max(np.mean(boot_frac < 0.80), 1 / n_boot) verd_6d = verdict(p_6d, fraction_below_mp >= 0.80) results.append(log_result( 'compound_symmetry_marchenko_pastur', fraction_below_mp, p_6d, verd_6d, description=( "Tests whether bulk eigenvalues (lambda_2 through lambda_J) are consistent with " "noise around a single common factor. Under EQUAL-weight compound symmetry, all bulk " "eigenvalues equal s^2-g. Under GENERAL weights (Paper 3b), the secular equation " "predicts unequal bulk eigenvalues interlacing the weight-dependent diagonal entries " "w_j = a_j^{-sigma}. Some eigenvalues above MP is expected for heterogeneous sector " "weights. Threshold: >= 80% below sigma^2*(1+sqrt(J/T))^2. Bootstrap (2000). " "Data: 18 FRED manufacturing IP series." ), details={ 'fraction_below_mp': round(fraction_below_mp, 4), 'mp_upper': round(mp_upper, 6), 'bulk_mean': round(bulk_mean, 6), 'q_ratio': round(q, 4), 'n_bulk_above_mp': int(np.sum(bulk > mp_upper)), 'J': J, 'T': T, } )) # ---- Test 6e: Sub-period stability ---- # Both halves should show eigenvalue dominance and low CV stable = (l1_early > 0.30 and l1_late > 0.30 and cv_early < 0.70 and cv_late < 0.70) # Bootstrap: resample each half and check stability boot_stable = [] n_early = len(g_early) n_late = len(g_late) for _ in range(n_boot): idx_e = rng.choice(n_early, size=n_early, replace=True) idx_l = rng.choice(n_late, size=n_late, replace=True) boot_e = pd.DataFrame(g_early.values[idx_e], columns=g_early.columns) boot_l = pd.DataFrame(g_late.values[idx_l], columns=g_late.columns) l1_e, cv_e, _ = period_stats(boot_e) l1_l, cv_l, _ = period_stats(boot_l) boot_stable.append(l1_e > 0.30 and l1_l > 0.30 and cv_e < 0.70 and cv_l < 0.70) boot_stable = np.array(boot_stable) # p-value: fraction of bootstrap where stability fails p_6e = max(1 - np.mean(boot_stable), 1 / n_boot) verd_6e = verdict(p_6e, stable) results.append(log_result( 'compound_symmetry_stability', float(stable), p_6e, verd_6e, description=( "Tests sub-period stability of compound symmetry: eigenvalue dominance (lambda_1 share " "> 0.30) and low correlation CV (< 0.70) should hold in both early and late halves of " "the sample. Bootstrap (2000 resamples per half) estimates stability probability. " "Data: 18 FRED manufacturing IP series split at midpoint." ), details={ 'l1_early': round(l1_early, 4), 'l1_late': round(l1_late, 4), 'cv_early': round(cv_early, 4), 'cv_late': round(cv_late, 4), 'boot_stability_rate': round(float(np.mean(boot_stable)), 4), 'midpoint': str(midpoint), 'J': J, 'T_early': n_early, 'T_late': n_late, } )) # ---- Test 6f: Variance heterogeneity (general-weight prediction) ---- # Under compound symmetry (equal weights), all sector variances are equal. # Under general weights, the CES Hessian diagonal is -(1-rho)/(c*p_j), # so Var(g_j) ∝ p_j = a_j^sigma — sectors with different weights have # different variances. Bartlett's test for homogeneity of variances. from scipy.stats import bartlett as bartlett_test sector_vars = np.array([growth[col].var() for col in growth.columns]) var_cv = np.std(sector_vars, ddof=1) / np.mean(sector_vars) if np.mean(sector_vars) > 0 else 0 bartlett_stat, bartlett_p = bartlett_test(*[growth[col].values for col in growth.columns]) bartlett_p = max(bartlett_p, 1e-15) # avoid exact zero # Under general weights, we WANT to reject H0 (variances are unequal). # sign_correct = True if variances ARE heterogeneous (CV > threshold) # The p-value from Bartlett is for H0: all variances equal. # Small p → reject H0 → supports general weights. verd_6f = verdict(bartlett_p, True) # sign is always correct (we want rejection) results.append(log_result( 'compound_symmetry_variance_heterogeneity', var_cv, bartlett_p, verd_6f, description=( "Tests the general-weight CES prediction (Paper 3b) that sector growth variances " "should differ. Under equal-weight compound symmetry, Var(g_j) = s^2 for all j. " "Under general weights, the Hessian diagonal H_jj = -(1-rho)/(c*p_j) where " "p_j = a_j^sigma, so Var(g_j) varies with CES weight (Lean: genCesHessianQF). " "Bartlett's test for homogeneity of variances: H0: all Var(g_j) equal. " "Rejection supports general weights over pure compound symmetry. " "Data: 18 FRED manufacturing IP series, monthly log-growth rates." ), details={ 'var_cv': round(var_cv, 4), 'bartlett_stat': round(bartlett_stat, 2), 'bartlett_p': bartlett_p, 'var_min': round(float(np.min(sector_vars)), 6), 'var_max': round(float(np.max(sector_vars)), 6), 'var_ratio': round(float(np.max(sector_vars) / np.min(sector_vars)), 2), 'J': J, 'T': T, } )) # ---- Test 6g: Variance-share correlation (general-weight quantitative test) ---- # Under general weights, the CES Hessian gives level variance Var(x_j) ∝ p_j, # but equilibrium level mu_j ∝ p_j too (cost-minimizing input allocation). # So for LOG-GROWTH: Var(dlog x_j) = Var(x_j)/mu_j^2 ∝ p_j/p_j^2 = 1/p_j. # Growth rate variance is INVERSELY proportional to CES weight. # Use BEA VA shares as proxy for CES weights a_j. # Prediction: Spearman(VA_share, Var(g_j)) < 0 (NEGATIVE — larger sectors # have lower growth volatility because their mean level dominates). va_shares = _load_bea_va_shares() if va_shares is not None: # Match sectors that are in both FRED IP data and BEA VA matched_shares = [] matched_vars = [] matched_labels = [] for col in growth.columns: if col in va_shares and va_shares[col] > 0: matched_shares.append(va_shares[col]) matched_vars.append(growth[col].var()) matched_labels.append(col) if len(matched_shares) >= 8: from scipy.stats import spearmanr r_sv, p_sv = spearmanr(matched_shares, matched_vars) verd_6g = verdict(p_sv, r_sv < 0) # negative = correct sign _cache['variance_share_data'] = { 'shares': matched_shares, 'vars': matched_vars, 'labels': matched_labels, 'r': r_sv, 'p': p_sv, } else: r_sv, p_sv = np.nan, np.nan verd_6g = 'SKIPPED' else: r_sv, p_sv = np.nan, np.nan verd_6g = 'SKIPPED' matched_shares = [] matched_vars = [] n_matched = len(matched_shares) if va_shares is not None else 0 herfindahl = sum(s**2 for s in matched_shares) if matched_shares else np.nan results.append(log_result( 'compound_symmetry_variance_share', r_sv, p_sv, verd_6g, description=( "Tests the general-weight CES prediction (Paper 3b, Lean: genCesHessianQF, " "secular_equation) that sectors with larger CES weights have LOWER log-growth " "variance: Var(dlog x_j) proportional to 1/p_j = a_j^{-sigma}. Level variance " "scales as p_j but mean level also scales as p_j, so growth rate variance = " "Var/mu^2 scales inversely. Spearman correlation of BEA value-added shares " "(proxy for a_j) vs growth variance. Prediction: r < 0 (negative). " "Data: 18 FRED manufacturing IP series matched to BEA quarterly VA shares " "(2005-2019). Herfindahl index of VA shares measures weight concentration." ), details={ 'spearman_r': round(r_sv, 4) if not np.isnan(r_sv) else None, 'n_matched': n_matched, 'herfindahl': round(herfindahl, 4) if not np.isnan(herfindahl) else None, 'equal_weight_H': round(1 / J, 4) if J > 0 else None, 'J': J, 'T': T, } )) # ---- Tests 6h-6i: Broad economy (BEA ~65 leaf industries) ---- broad_growth, broad_shares = _load_bea_broad_data() broad_skip_desc = ( "Tests compound symmetry predictions across ~65 BEA leaf industries (all sectors, " "quarterly VA 2005-2023). Extends the 18-sector FRED manufacturing tests to the " "full economy with 3.5x more industries for greater statistical power." ) if broad_growth is not None and broad_shares is not None: J_broad = broad_growth.shape[1] T_broad = broad_growth.shape[0] # ---- Test 6h: Broad eigenvalue dominance ---- broad_cov = broad_growth.cov().values broad_evals = np.sort(np.linalg.eigh(broad_cov)[0])[::-1] broad_lam1_share = broad_evals[0] / np.sum(broad_evals) # Permutation test rng_broad = np.random.RandomState(123) n_perm_broad = 2000 broad_null_shares = [] broad_vals = broad_growth.values.copy() for _ in range(n_perm_broad): pv = np.column_stack([rng_broad.permutation(broad_vals[:, j]) for j in range(J_broad)]) pc = np.cov(pv.T) pe = np.sort(np.linalg.eigh(pc)[0])[::-1] broad_null_shares.append(pe[0] / np.sum(pe)) broad_null_shares = np.array(broad_null_shares) p_6h = max(np.mean(broad_null_shares >= broad_lam1_share), 1 / n_perm_broad) verd_6h = verdict(p_6h, broad_lam1_share > np.median(broad_null_shares)) _cache['broad_lam1_share'] = broad_lam1_share _cache['broad_null_shares'] = broad_null_shares _cache['broad_evals'] = broad_evals _cache['J_broad'] = J_broad _cache['T_broad'] = T_broad results.append(log_result( 'compound_symmetry_broad_eigenvalue', broad_lam1_share, p_6h, verd_6h, description=( "Tests eigenvalue dominance (rank-1 structure) across the full US economy, " "not just manufacturing. Under CES aggregation, a single market mode should " "dominate the covariance spectrum regardless of sector count. Permutation test " "(2000 draws, columns independently permuted) against no-structure null. " "Data: ~65 BEA leaf industries, quarterly real VA log-growth (2005-2023). " f"J/T={J_broad}/{T_broad}={J_broad/T_broad:.2f} — high ratio makes " "dominance harder to detect, so finding it is more impressive." ), details={ 'lam1_share': round(broad_lam1_share, 4), 'null_median': round(float(np.median(broad_null_shares)), 4), 'null_max': round(float(np.max(broad_null_shares)), 4), 'J': J_broad, 'T': T_broad, 'JT_ratio': round(J_broad / T_broad, 3), } )) # ---- Test 6i: Broad variance-share Spearman ---- broad_sector_vars = {col: broad_growth[col].var() for col in broad_growth.columns} share_list = [] var_list = [] label_list = [] for ind in broad_growth.columns: if ind in broad_shares and broad_shares[ind] > 0: share_list.append(broad_shares[ind]) var_list.append(broad_sector_vars[ind]) label_list.append(ind) if len(share_list) >= 10: from scipy.stats import spearmanr r_broad, p_broad = spearmanr(share_list, var_list) verd_6i = verdict(p_broad, r_broad < 0) _cache['broad_variance_share_data'] = { 'shares': share_list, 'vars': var_list, 'labels': label_list, 'r': r_broad, 'p': p_broad, } else: r_broad, p_broad = np.nan, np.nan verd_6i = 'SKIPPED' n_broad = len(share_list) herf_broad = sum(s**2 for s in share_list) if share_list else np.nan results.append(log_result( 'compound_symmetry_broad_variance_share', r_broad, p_broad, verd_6i, description=( "Tests the general-weight CES prediction Var(dlog x_j) proportional to 1/p_j across " f"{n_broad} BEA leaf industries spanning the full US economy (not just manufacturing). " "This is the broad-economy version of test 6g, with ~3.5x more sectors for greater " "statistical power. Spearman correlation of average VA share (proxy for CES weight) " "vs quarterly log-growth variance. Prediction: r < 0 (larger sectors less volatile). " "Data: BEA quarterly real VA (2005-2023), all private non-aggregate industries." ), details={ 'spearman_r': round(r_broad, 4) if not np.isnan(r_broad) else None, 'n_industries': n_broad, 'herfindahl': round(herf_broad, 4) if not np.isnan(herf_broad) else None, 'J_broad': J_broad, 'T_broad': T_broad, } )) # ---- Test 6j: Pairwise correlation vs sector size ---- # Under rank-1 Hessian (genCesHessianQF), the common factor loading is the # same for all sectors, but idiosyncratic variance ∝ 1/p_j. So larger sectors # (higher p_j, lower idiosyncratic variance) have higher average pairwise # correlation. Test: Spearman(share_j, avg_cor_j) > 0. from scipy.stats import spearmanr broad_corr = broad_growth.corr().values shares_arr = np.array([broad_shares[col] for col in broad_growth.columns]) avg_corr_per_sector = [] for j in range(J_broad): off_diag = [broad_corr[j, k] for k in range(J_broad) if k != j] avg_corr_per_sector.append(np.mean(off_diag)) avg_corr_per_sector = np.array(avg_corr_per_sector) r_6j, p_6j = spearmanr(shares_arr, avg_corr_per_sector) verd_6j = verdict(p_6j, r_6j > 0) # positive = correct sign _cache['correlation_size_data'] = { 'shares': shares_arr.tolist(), 'avg_corr': avg_corr_per_sector.tolist(), 'labels': list(broad_growth.columns), 'r': r_6j, 'p': p_6j, } results.append(log_result( 'compound_symmetry_correlation_size', r_6j, p_6j, verd_6j, description=( "Tests the general-weight CES prediction that larger sectors have higher " "average pairwise correlations. Under the rank-1 Hessian (Lean: genCesHessianQF), " "all sectors share the same common factor loading, but idiosyncratic variance " "is 1/p_j. Larger sectors (higher p_j) have less idiosyncratic noise relative " "to the common factor, so their correlations with other sectors are higher. " "Spearman of VA share vs row-mean of correlation matrix (excluding diagonal). " f"Prediction: r > 0 (positive). Data: {J_broad} BEA industries, quarterly VA." ), details={ 'spearman_r': round(r_6j, 4), 'n_industries': J_broad, 'avg_corr_range': [round(float(avg_corr_per_sector.min()), 4), round(float(avg_corr_per_sector.max()), 4)], } )) # ---- Test 6k: Eigenmode diagonality ---- # Paper 3b Theorem weighted_vri: the covariance matrix is DIAGONAL in the # eigenmode basis of diag(w) restricted to 1⊥, where w_j = 1/p_j. # Test: rotate observed covariance by predicted eigenvectors, measure # off-diagonal energy. Should be lower than permutation null. w_arr = 1.0 / shares_arr # Project onto 1⊥: P = I - (1/J)11^T P_proj = np.eye(J_broad) - np.ones((J_broad, J_broad)) / J_broad # Restricted diag(w) W_restricted = P_proj @ np.diag(w_arr) @ P_proj # Eigendecompose — first J-1 eigenvectors (last has eigenvalue ~0) w_evals, w_evecs = np.linalg.eigh(W_restricted) # Sort descending by absolute eigenvalue w_order = np.argsort(np.abs(w_evals))[::-1] w_evals = w_evals[w_order] w_evecs = w_evecs[:, w_order] pred_evecs = w_evecs[:, :J_broad - 1] # drop the ~0 eigenvalue mode # Rotate observed covariance by predicted eigenvectors Sigma_rotated = pred_evecs.T @ broad_cov @ pred_evecs # Off-diagonal energy ratio diag_mask = np.eye(J_broad - 1, dtype=bool) off_diag_energy = np.sum(Sigma_rotated[~diag_mask] ** 2) total_energy = np.sum(Sigma_rotated ** 2) obs_offdiag_ratio = np.sqrt(off_diag_energy / total_energy) # Permutation test: permute share labels, recompute predicted eigenvectors rng_diag = np.random.RandomState(456) n_perm_diag = 2000 null_offdiag_ratios = [] for _ in range(n_perm_diag): perm_w = rng_diag.permutation(w_arr) pW = P_proj @ np.diag(perm_w) @ P_proj pe, pv = np.linalg.eigh(pW) po = np.argsort(np.abs(pe))[::-1] pv = pv[:, po][:, :J_broad - 1] Sr = pv.T @ broad_cov @ pv oe = np.sum(Sr[~diag_mask] ** 2) te = np.sum(Sr ** 2) null_offdiag_ratios.append(np.sqrt(oe / te)) null_offdiag_ratios = np.array(null_offdiag_ratios) # One-sided: observed ratio should be LOWER than null p_6k = max(np.mean(null_offdiag_ratios <= obs_offdiag_ratio), 1 / n_perm_diag) verd_6k = verdict(p_6k, obs_offdiag_ratio < np.median(null_offdiag_ratios)) _cache['eigenmode_diag_data'] = { 'obs_ratio': obs_offdiag_ratio, 'null_ratios': null_offdiag_ratios, } results.append(log_result( 'compound_symmetry_eigenmode_diagonality', obs_offdiag_ratio, p_6k, verd_6k, description=( "Tests the weighted VRI prediction (Paper 3b, Theorem weighted_vri, " "Lean: genCesHessianQF): the covariance matrix should be approximately " "DIAGONAL when rotated into the eigenmode basis predicted by sector weights. " "Eigenvectors of diag(1/VA_share) restricted to 1-perp provide the predicted " "basis. Off-diagonal Frobenius energy ratio measures diagonality. Permutation " "test (2000 draws, share labels shuffled) provides null. Observed ratio should " f"be lower than null. Data: {J_broad} BEA industries, quarterly VA." ), details={ 'obs_offdiag_ratio': round(obs_offdiag_ratio, 4), 'null_median': round(float(np.median(null_offdiag_ratios)), 4), 'null_percentile': round(float(np.mean( null_offdiag_ratios <= obs_offdiag_ratio) * 100), 1), 'J': J_broad, } )) # ---- Test 6l: Secular equation eigenvalue ordering ---- # The secular equation Σ 1/(w_j - μ) = 0 predicts the eigenvalues of the # Hessian on 1⊥. Covariance eigenvalues ∝ 1/μ_k. Compare predicted ordering # against observed eigenvalues of the covariance restricted to 1⊥. from scipy.optimize import brentq w_sorted = np.sort(w_arr) # Add small jitter to break exact ties for i in range(1, len(w_sorted)): if w_sorted[i] <= w_sorted[i - 1]: w_sorted[i] = w_sorted[i - 1] + 1e-8 def secular_fn(mu, w_vals): return np.sum(1.0 / (w_vals - mu)) secular_roots = [] for k in range(len(w_sorted) - 1): lo = w_sorted[k] + 1e-10 hi = w_sorted[k + 1] - 1e-10 if hi <= lo: continue try: root = brentq(secular_fn, lo, hi, args=(w_sorted,)) secular_roots.append(root) except (ValueError, RuntimeError): pass # Restrict observed covariance to 1⊥ and get eigenvalues Sigma_restricted = P_proj @ broad_cov @ P_proj obs_evals_full = np.sort(np.linalg.eigh(Sigma_restricted)[0])[::-1] obs_evals_1perp = obs_evals_full[:J_broad - 1] # drop the ~0 eigenvalue n_roots = min(len(secular_roots), len(obs_evals_1perp)) if n_roots >= 10: # Predicted covariance eigenvalues ∝ 1/μ_k, sorted descending pred_cov_evals = 1.0 / np.array(secular_roots[:n_roots]) pred_cov_sorted = np.sort(pred_cov_evals)[::-1] obs_sorted = obs_evals_1perp[:n_roots] r_6l, p_6l = spearmanr(pred_cov_sorted, obs_sorted) verd_6l = verdict(p_6l, r_6l > 0) # positive = correct ordering _cache['secular_eigenvalue_data'] = { 'pred_evals': pred_cov_sorted.tolist(), 'obs_evals': obs_sorted.tolist(), 'n_roots': n_roots, 'r': r_6l, 'p': p_6l, } else: r_6l, p_6l = np.nan, np.nan verd_6l = 'SKIPPED' results.append(log_result( 'compound_symmetry_secular_eigenvalues', r_6l, p_6l, verd_6l, description=( "Tests the secular equation prediction (Lean: secular_equation, Paper 1 " "Proposition 10.2): eigenvalues of diag(w_j) restricted to 1-perp, where " "w_j = 1/VA_share_j, determine the Hessian eigenmode spectrum. Covariance " "eigenvalues are proportional to 1/mu_k (inverse secular roots). Spearman rank " "correlation of predicted vs observed eigenvalue ordering on 1-perp. " f"Prediction: r > 0 (positive). Data: {J_broad} BEA industries, {n_roots} " "secular roots solved via Brent's method." ), details={ 'spearman_r': round(r_6l, 4) if not np.isnan(r_6l) else None, 'n_roots': n_roots if n_roots >= 10 else 0, 'J': J_broad, } )) else: results.append(log_result( 'compound_symmetry_broad_eigenvalue', np.nan, np.nan, 'SKIPPED', description=broad_skip_desc)) results.append(log_result( 'compound_symmetry_broad_variance_share', np.nan, np.nan, 'SKIPPED', description=broad_skip_desc)) for name in ['compound_symmetry_correlation_size', 'compound_symmetry_eigenmode_diagonality', 'compound_symmetry_secular_eigenvalues']: results.append(log_result(name, np.nan, np.nan, 'SKIPPED', description=broad_skip_desc)) # ---- Test 6m: Manufacturing monthly stability (FRED, J=18, T~648) ---- # Same test as 6e but with monthly data (T~324 per half vs ~41 quarterly). # Much higher power: T/J ratio goes from ~2 to ~18 per half. mid_m = len(growth) // 2 g_early_m = growth.iloc[:mid_m] g_late_m = growth.iloc[mid_m:] l1_early_m, cv_early_m, n_early_m = period_stats(g_early_m) l1_late_m, cv_late_m, n_late_m = period_stats(g_late_m) stable_m = (l1_early_m > 0.30 and l1_late_m > 0.30 and cv_early_m < 0.70 and cv_late_m < 0.70) # Bootstrap boot_stable_m = [] rng_m = np.random.RandomState(789) for _ in range(n_boot): idx_e = rng_m.choice(n_early_m, size=n_early_m, replace=True) idx_l = rng_m.choice(n_late_m, size=n_late_m, replace=True) boot_e = pd.DataFrame(g_early_m.values[idx_e], columns=g_early_m.columns) boot_l = pd.DataFrame(g_late_m.values[idx_l], columns=g_late_m.columns) l1_e, cv_e, _ = period_stats(boot_e) l1_l, cv_l, _ = period_stats(boot_l) boot_stable_m.append(l1_e > 0.30 and l1_l > 0.30 and cv_e < 0.70 and cv_l < 0.70) boot_stable_m = np.array(boot_stable_m) p_6m = max(1 - np.mean(boot_stable_m), 1 / n_boot) verd_6m = verdict(p_6m, stable_m) _cache['mfg_stability_data'] = { 'l1_early': l1_early_m, 'l1_late': l1_late_m, 'cv_early': cv_early_m, 'cv_late': cv_late_m, 'boot_rate': float(np.mean(boot_stable_m)), } results.append(log_result( 'compound_symmetry_mfg_stability', float(stable_m), p_6m, verd_6m, description=( "Tests sub-period stability using FRED monthly IP data (higher power than " "quarterly 6e). Splits ~648 months into two halves (~324 each). Each half " "must show eigenvalue dominance (lambda_1 share > 0.30) and low correlation " "CV (< 0.70). Bootstrap (2000 resamples per half). T/J ratio ~18 per half " "vs ~2 in quarterly test. Data: 18 FRED manufacturing IP series, monthly." ), details={ 'l1_early': round(l1_early_m, 4), 'l1_late': round(l1_late_m, 4), 'cv_early': round(cv_early_m, 4), 'cv_late': round(cv_late_m, 4), 'boot_stability_rate': round(float(np.mean(boot_stable_m)), 4), 'T_early': n_early_m, 'T_late': n_late_m, 'J': J, } )) # ---- Test 6n: Manufacturing eigenmode diagonality (FRED monthly, J/T=0.028) ---- # Same test as 6k but using monthly FRED data (J/T=0.028 vs 0.80 for BEA quarterly). # Much better eigenmode estimation due to low J/T ratio. if va_shares is not None: share_arr_m = np.array([va_shares.get(sid, 0) for sid in growth.columns]) valid_m = share_arr_m > 0 if np.sum(valid_m) >= 5: # Use only sectors with valid VA shares cols_m = [c for c, v in zip(growth.columns, valid_m) if v] growth_m = growth[cols_m] share_m = share_arr_m[valid_m] J_m = len(cols_m) T_m = len(growth_m) # Compute w_j = 1/share_j (inverse share ~ 1/p_j) w_m = 1.0 / share_m # Project diag(w) onto 1-perp ones_m = np.ones(J_m) P_m = np.eye(J_m) - np.outer(ones_m, ones_m) / J_m Pw_m = P_m @ np.diag(w_m) @ P_m # Eigenvectors of projected weight matrix w_evals_m, w_evecs_m = np.linalg.eigh(Pw_m) w_order_m = np.argsort(np.abs(w_evals_m))[::-1] w_evecs_m = w_evecs_m[:, w_order_m] pred_evecs_m = w_evecs_m[:, :J_m - 1] # drop ~0 eigenvalue mode # Observed covariance rotated by predicted eigenvectors cov_m = growth_m.cov().values Sigma_rot_m = pred_evecs_m.T @ cov_m @ pred_evecs_m diag_mask_m = np.eye(J_m - 1, dtype=bool) off_energy_m = np.sum(Sigma_rot_m[~diag_mask_m] ** 2) total_energy_m = np.sum(Sigma_rot_m ** 2) obs_ratio_m = np.sqrt(off_energy_m / total_energy_m) # Permutation test: shuffle share labels rng_n = np.random.RandomState(101) n_perm_n = 2000 null_ratios_m = [] for _ in range(n_perm_n): perm_w = rng_n.permutation(w_m) pW = P_m @ np.diag(perm_w) @ P_m pe, pv = np.linalg.eigh(pW) po = np.argsort(np.abs(pe))[::-1] pv = pv[:, po][:, :J_m - 1] Sr = pv.T @ cov_m @ pv oe = np.sum(Sr[~diag_mask_m] ** 2) te = np.sum(Sr ** 2) null_ratios_m.append(np.sqrt(oe / te)) null_ratios_m = np.array(null_ratios_m) # One-sided: observed should be LOWER than null p_6n = max(np.mean(null_ratios_m <= obs_ratio_m), 1 / n_perm_n) verd_6n = verdict(p_6n, obs_ratio_m < np.median(null_ratios_m)) _cache['mfg_eigenmode_data'] = { 'obs_ratio': obs_ratio_m, 'null_ratios': null_ratios_m, 'J': J_m, 'T': T_m, } results.append(log_result( 'compound_symmetry_mfg_eigenmode_diag', obs_ratio_m, p_6n, verd_6n, description=( "Tests weighted VRI diagonality using FRED monthly IP (higher power than " "quarterly BEA 6k). J/T ratio = 0.028 (vs 0.80 in quarterly test), giving " "much more reliable eigenmode estimation. Eigenvectors of diag(1/VA_share) " "restricted to 1-perp predict the covariance eigenmode basis. Off-diagonal " "Frobenius energy ratio measures diagonality. Permutation test (2000 draws, " f"share labels shuffled). Data: {J_m} FRED sectors, {T_m} monthly obs." ), details={ 'obs_offdiag_ratio': round(obs_ratio_m, 4), 'null_median': round(float(np.median(null_ratios_m)), 4), 'null_percentile': round(float(np.mean( null_ratios_m <= obs_ratio_m) * 100), 1), 'J': J_m, 'T': T_m, 'JT_ratio': round(J_m / T_m, 4), } )) else: results.append(log_result( 'compound_symmetry_mfg_eigenmode_diag', np.nan, np.nan, 'SKIPPED', description="Insufficient VA share data for FRED eigenmode diagonality test.")) else: results.append(log_result( 'compound_symmetry_mfg_eigenmode_diag', np.nan, np.nan, 'SKIPPED', description="BEA VA shares not available for FRED eigenmode diagonality test.")) # Store cache for figures _cache['eigenvalues'] = eigenvalues _cache['off_diag_corr'] = off_diag_corr _cache['corr_matrix'] = corr_matrix _cache['lam1_share'] = lam1_share _cache['null_shares'] = null_shares _cache['corr_cv'] = corr_cv _cache['frob_ratio'] = frob_ratio _cache['J'] = J _cache['T'] = T _cache['bulk'] = bulk _cache['bulk_mean'] = bulk_mean _cache['mp_upper'] = mp_upper _cache['null_cvs'] = null_cvs _cache['null_shares_6a'] = null_shares _cache['g_implied'] = g_implied _cache['s2_implied'] = s2_implied _cache['l1_early'] = l1_early _cache['l1_late'] = l1_late _cache['cv_early'] = cv_early _cache['cv_late'] = cv_late _cache['fraction_below_mp'] = fraction_below_mp _cache['sector_vars'] = sector_vars _cache['sector_labels'] = list(growth.columns) _cache['bartlett_stat'] = bartlett_stat _cache['bartlett_p'] = bartlett_p return results def figures(): """Generate per-subtest figures for compound symmetry test.""" if not _cache: main() paths = [] eigenvalues = _cache.get('eigenvalues') if eigenvalues is None: return paths import matplotlib.pyplot as plt J = _cache['J'] bulk = _cache['bulk'] mp_upper = _cache['mp_upper'] bulk_mean = _cache['bulk_mean'] off_diag_corr = _cache['off_diag_corr'] # --- Figure 1: Eigenvalue dominance (6a) --- fig, axes = setup_subplots(1, 2, width=12, height=5) lam1_share = _cache['lam1_share'] # (a) Eigenvalue spectrum ax = axes[0] ax.bar(range(1, len(eigenvalues) + 1), eigenvalues, color='steelblue', edgecolor='k') ax.axhline(mp_upper, color='r', linestyle='--', label=f'MP upper = {mp_upper:.4f}') ax.axhline(bulk_mean, color='orange', linestyle=':', label=f'Bulk mean = {bulk_mean:.4f}') ax.set_xlabel('Eigenvalue index') ax.set_ylabel('Eigenvalue') ax.set_title(f'(a) Eigenvalue Spectrum (J={J})') ax.legend(fontsize=8) ax.grid(True, alpha=0.3) # (b) Permutation null for lambda1 share ax = axes[1] null_shares = _cache.get('null_shares', np.array([])) if len(null_shares) > 0: ax.hist(null_shares, bins=50, color='lightcoral', edgecolor='k', alpha=0.7) ax.axvline(lam1_share, color='blue', linewidth=2, label=f'Observed = {lam1_share:.3f}') ax.set_xlabel('λ₁ / Σλ') ax.set_ylabel('Count') ax.set_title('(b) Permutation null for λ₁ share') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) fig.suptitle(f'Eigenvalue Dominance (λ₁ share={lam1_share:.3f})', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'compound_symmetry_eigenvalue_dominance', subfolder='foundations')) # --- Figure 2: Equicorrelation (6b) --- fig, axes = setup_subplots(1, 2, width=12, height=5) # (a) Off-diagonal correlation histogram ax = axes[0] ax.hist(off_diag_corr, bins=30, color='forestgreen', edgecolor='k', alpha=0.7) mean_od = np.mean(off_diag_corr) ax.axvline(mean_od, color='r', linestyle='--', label=f'mean = {mean_od:.3f}') ax.set_xlabel('Off-diagonal correlation') ax.set_ylabel('Count') ax.set_title(f'(a) Off-diagonal correlations (CV={_cache["corr_cv"]:.3f})') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) # (b) Permutation null distribution ax = axes[1] null_cvs = _cache['null_cvs'] ax.hist(null_cvs, bins=50, color='lightcoral', edgecolor='k', alpha=0.7) ax.axvline(_cache['corr_cv'], color='blue', linewidth=2, label=f'Observed CV = {_cache["corr_cv"]:.3f}') ax.set_xlabel('CV of off-diagonal correlations') ax.set_ylabel('Count') ax.set_title('(b) Permutation null distribution') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) fig.suptitle('Equicorrelation Test', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'compound_symmetry_equicorrelation', subfolder='foundations')) # --- Figure 3: Implied parameters (6c) --- fig, ax = setup_figure() g = _cache['g_implied'] s2 = _cache['s2_implied'] bars = ax.bar(['g (common)', 's² (total)', 's²−g (idio)'], [g, s2, s2 - g], color=['steelblue', 'orange', 'forestgreen'], edgecolor='k') ax.axhline(0, color='k', linewidth=0.5) ax.set_ylabel('Variance component') ax.set_title(f'Implied Compound Symmetry Parameters (g={g:.5f}, s²={s2:.5f})') for bar, val in zip(bars, [g, s2, s2 - g]): ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.00001, f'{val:.5f}', ha='center', va='bottom', fontsize=9) paths.extend(save_figure(fig, 'compound_symmetry_implied_parameters', subfolder='foundations')) # --- Figure 4: Marchenko-Pastur (6d) --- fig, ax = setup_figure() ax.bar(range(2, len(bulk) + 2), bulk, color='steelblue', edgecolor='k', alpha=0.7) ax.axhline(mp_upper, color='r', linestyle='--', linewidth=2, label=f'MP upper = {mp_upper:.5f}') ax.axhline(bulk_mean, color='orange', linestyle=':', linewidth=2, label=f'Bulk mean = {bulk_mean:.5f}') n_above = int(np.sum(bulk > mp_upper)) frac = _cache['fraction_below_mp'] ax.set_xlabel('Eigenvalue index (2 through J)') ax.set_ylabel('Eigenvalue') ax.set_title(f'Bulk Eigenvalues vs MP Bound ({n_above}/{len(bulk)} above, ' f'{frac:.1%} below)') ax.legend(fontsize=9) paths.extend(save_figure(fig, 'compound_symmetry_marchenko_pastur', subfolder='foundations')) # --- Figure 5: Sub-period stability (6e) --- fig, axes = setup_subplots(1, 2, width=12, height=5) for ax, label, l1, cv in [(axes[0], 'Early', _cache['l1_early'], _cache['cv_early']), (axes[1], 'Late', _cache['l1_late'], _cache['cv_late'])]: ax.bar(['λ₁ share', 'Corr CV'], [l1, cv], color=['steelblue', 'forestgreen'], edgecolor='k') ax.axhline(0.30, color='r', linestyle='--', alpha=0.5, label='λ₁ threshold (0.30)') ax.axhline(0.70, color='orange', linestyle='--', alpha=0.5, label='CV threshold (0.70)') ax.set_ylim(0, max(l1, cv, 0.75) * 1.2) ax.set_title(f'{label} half: λ₁={l1:.3f}, CV={cv:.3f}') ax.legend(fontsize=8) ax.grid(True, alpha=0.3) fig.suptitle('Sub-Period Stability', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'compound_symmetry_stability', subfolder='foundations')) # --- Figure 6: Variance heterogeneity (6f) --- sector_vars = _cache.get('sector_vars') if sector_vars is not None: fig, ax = setup_figure(width=10, height=5) labels = _cache.get('sector_labels', [f'S{i}' for i in range(len(sector_vars))]) short_labels = [s.replace('IPG', '').replace('S', '') for s in labels] order = np.argsort(sector_vars)[::-1] ax.bar(range(len(sector_vars)), sector_vars[order], color='steelblue', edgecolor='k', alpha=0.7) ax.set_xticks(range(len(sector_vars))) ax.set_xticklabels([short_labels[i] for i in order], rotation=45, ha='right', fontsize=8) mean_var = np.mean(sector_vars) ax.axhline(mean_var, color='r', linestyle='--', label=f'Mean variance = {mean_var:.5f}') bp = _cache.get('bartlett_p', np.nan) ax.set_xlabel('NAICS Sector') ax.set_ylabel('Variance of monthly log-growth') ax.set_title(f'Sector Growth Variance (Bartlett p={bp:.2e}, ' f'ratio max/min={np.max(sector_vars)/np.min(sector_vars):.1f}x)') ax.legend(fontsize=9) paths.extend(save_figure(fig, 'compound_symmetry_variance_heterogeneity', subfolder='foundations')) # --- Figure 7: Variance-share Spearman (6g) --- vs_data = _cache.get('variance_share_data') if vs_data is not None: fig, ax = setup_figure() shares = np.array(vs_data['shares']) * 100 # to percent vars_ = np.array(vs_data['vars']) ax.scatter(shares, vars_, color='steelblue', edgecolor='k', s=60, zorder=5) for i, lbl in enumerate(vs_data['labels']): short = lbl.replace('IPG', '').replace('S', '') ax.annotate(short, (shares[i], vars_[i]), fontsize=7, xytext=(4, 4), textcoords='offset points') # Fit line z = np.polyfit(shares, vars_, 1) x_fit = np.linspace(min(shares), max(shares), 100) ax.plot(x_fit, np.polyval(z, x_fit), 'r--', alpha=0.7) r_val = vs_data['r'] p_val = vs_data['p'] ax.set_xlabel('BEA Value-Added Share (%)') ax.set_ylabel('Variance of monthly log-growth') ax.set_title(f'Variance vs VA Share (Spearman r={r_val:.3f}, p={p_val:.3f})') ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'compound_symmetry_variance_share', subfolder='foundations')) # --- Figure 8: Broad eigenvalue dominance (6h) --- broad_null = _cache.get('broad_null_shares') if broad_null is not None: fig, axes = setup_subplots(1, 2, width=12, height=5) broad_evals = _cache['broad_evals'] J_b = _cache['J_broad'] # (a) Eigenvalue spectrum ax = axes[0] ax.bar(range(1, min(31, len(broad_evals) + 1)), broad_evals[:30], color='steelblue', edgecolor='k', alpha=0.7) ax.set_xlabel('Eigenvalue index (top 30)') ax.set_ylabel('Eigenvalue') ax.set_title(f'(a) BEA {J_b}-Industry Eigenvalue Spectrum') ax.grid(True, alpha=0.3) # (b) Permutation null ax = axes[1] broad_lam1 = _cache['broad_lam1_share'] ax.hist(broad_null, bins=50, color='lightcoral', edgecolor='k', alpha=0.7) ax.axvline(broad_lam1, color='blue', linewidth=2, label=f'Observed = {broad_lam1:.3f}') ax.set_xlabel('Lambda-1 / sum(lambda)') ax.set_ylabel('Count') ax.set_title('(b) Permutation null distribution') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) fig.suptitle(f'Broad Economy Eigenvalue Dominance ({J_b} industries)', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'compound_symmetry_broad_eigenvalue', subfolder='foundations')) # --- Figure 9: Broad variance-share (6i) --- bvs_data = _cache.get('broad_variance_share_data') if bvs_data is not None: fig, ax = setup_figure(width=10, height=6) shares_pct = np.array(bvs_data['shares']) * 100 vars_b = np.array(bvs_data['vars']) ax.scatter(shares_pct, vars_b, color='steelblue', edgecolor='k', s=40, zorder=5) for i, lbl in enumerate(bvs_data['labels']): ax.annotate(lbl, (shares_pct[i], vars_b[i]), fontsize=5, xytext=(3, 3), textcoords='offset points') z = np.polyfit(shares_pct, vars_b, 1) x_fit = np.linspace(min(shares_pct), max(shares_pct), 100) ax.plot(x_fit, np.polyval(z, x_fit), 'r--', alpha=0.7) r_b = bvs_data['r'] p_b = bvs_data['p'] ax.set_xlabel('BEA Value-Added Share (%)') ax.set_ylabel('Variance of quarterly log-growth') ax.set_title(f'Broad Economy: Variance vs VA Share ' f'(N={len(bvs_data["labels"])}, Spearman r={r_b:.3f}, p={p_b:.4f})') ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'compound_symmetry_broad_variance_share', subfolder='foundations')) # --- Figure 10: Correlation vs sector size (6j) --- cs_data = _cache.get('correlation_size_data') if cs_data is not None: fig, ax = setup_figure(width=10, height=6) sh = np.array(cs_data['shares']) * 100 ac = np.array(cs_data['avg_corr']) ax.scatter(sh, ac, color='steelblue', edgecolor='k', s=40, zorder=5) for i, lbl in enumerate(cs_data['labels']): ax.annotate(lbl, (sh[i], ac[i]), fontsize=5, xytext=(3, 3), textcoords='offset points') z = np.polyfit(sh, ac, 1) x_fit = np.linspace(min(sh), max(sh), 100) ax.plot(x_fit, np.polyval(z, x_fit), 'r--', alpha=0.7) ax.set_xlabel('BEA Value-Added Share (%)') ax.set_ylabel('Average pairwise correlation') ax.set_title(f'Correlation vs Sector Size ' f'(N={len(cs_data["labels"])}, ' f'Spearman r={cs_data["r"]:.3f}, p={cs_data["p"]:.4f})') ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'compound_symmetry_correlation_size', subfolder='foundations')) # --- Figure 11: Eigenmode diagonality (6k) --- em_data = _cache.get('eigenmode_diag_data') if em_data is not None: fig, ax = setup_figure() null_r = em_data['null_ratios'] obs_r = em_data['obs_ratio'] ax.hist(null_r, bins=50, color='lightcoral', edgecolor='k', alpha=0.7, label='Permutation null') ax.axvline(obs_r, color='blue', linewidth=2, label=f'Observed = {obs_r:.3f}') ax.set_xlabel('Off-diagonal Frobenius energy ratio') ax.set_ylabel('Count') pctile = np.mean(null_r <= obs_r) * 100 ax.set_title(f'Eigenmode Diagonality Test ' f'(observed at {pctile:.1f}th percentile)') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'compound_symmetry_eigenmode_diagonality', subfolder='foundations')) # --- Figure 12: Secular equation eigenvalues (6l) --- se_data = _cache.get('secular_eigenvalue_data') if se_data is not None: fig, ax = setup_figure() pred = np.array(se_data['pred_evals']) obs = np.array(se_data['obs_evals']) # Normalize both to [0, 1] for visual comparison pred_norm = pred / pred.max() if pred.max() > 0 else pred obs_norm = obs / obs.max() if obs.max() > 0 else obs n = len(pred_norm) ax.scatter(range(1, n + 1), pred_norm, color='steelblue', s=30, label='Predicted (1/secular roots)', zorder=5) ax.scatter(range(1, n + 1), obs_norm, color='orange', s=30, marker='^', label='Observed (covariance)', zorder=5) ax.set_xlabel('Eigenvalue rank (descending)') ax.set_ylabel('Normalized eigenvalue') ax.set_title(f'Secular Equation Eigenvalue Ordering ' f'(N={n}, Spearman r={se_data["r"]:.3f}, p={se_data["p"]:.4f})') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'compound_symmetry_secular_eigenvalues', subfolder='foundations')) # --- Figure 13: Manufacturing monthly stability (6m) --- ms_data = _cache.get('mfg_stability_data') if ms_data is not None: fig, axes = setup_subplots(1, 2, width=12, height=5) ax = axes[0] labels_m = ['Early', 'Late'] l1_vals = [ms_data['l1_early'], ms_data['l1_late']] colors_m = ['steelblue' if v > 0.30 else 'salmon' for v in l1_vals] ax.bar(labels_m, l1_vals, color=colors_m, edgecolor='k') ax.axhline(0.30, color='r', linestyle='--', label='Threshold (0.30)') ax.set_ylabel('Lambda-1 share') ax.set_title('(a) Eigenvalue Dominance by Half') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) ax = axes[1] cv_vals = [ms_data['cv_early'], ms_data['cv_late']] colors_cv = ['steelblue' if v < 0.70 else 'salmon' for v in cv_vals] ax.bar(labels_m, cv_vals, color=colors_cv, edgecolor='k') ax.axhline(0.70, color='r', linestyle='--', label='Threshold (0.70)') ax.set_ylabel('Correlation CV') ax.set_title('(b) Correlation CV by Half') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) fig.suptitle(f'Monthly Manufacturing Stability ' f'(boot rate={ms_data["boot_rate"]:.3f})', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'compound_symmetry_mfg_stability', subfolder='foundations')) # --- Figure 14: Manufacturing eigenmode diagonality (6n) --- me_data = _cache.get('mfg_eigenmode_data') if me_data is not None: fig, ax = setup_figure() null_r = me_data['null_ratios'] obs_r = me_data['obs_ratio'] ax.hist(null_r, bins=50, color='lightcoral', edgecolor='k', alpha=0.7, label='Permutation null') ax.axvline(obs_r, color='blue', linewidth=2, label=f'Observed = {obs_r:.3f}') ax.set_xlabel('Off-diagonal Frobenius energy ratio') ax.set_ylabel('Count') pctile = np.mean(null_r <= obs_r) * 100 jt = me_data['J'] / me_data['T'] ax.set_title(f'Mfg Eigenmode Diagonality (J/T={jt:.3f}, ' f'{pctile:.1f}th percentile)') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths.extend(save_figure(fig, 'compound_symmetry_mfg_eigenmode_diag', 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']})")

J	18
T	648
n pass	4
corr cv	0.3468
pass 6a	True
pass 6b	True
pass 6c	True
pass 6d	False
pass 6e	True
g implied	1.847e-4
statistic	0.5012

Compound Symmetry

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources