#!/usr/bin/env python3 """ Eigenvalue Early Warning Signals -- Paper 2, Prop critical_slowing + Remark early_warning Paper 3, early_warning_signals Theory (Paper 2, Section 6.5; Paper 3, Section 3.4): The CES potential Hessian eigenvalues on 1-perp are proportional to K_eff = K(1 - T/T*). As information friction T -> T*, the relaxation time tau = tau_0 / (1 - T/T*) diverges, producing three observable early warning signals (Remark, Paper 2): (a) Increased autocorrelation: rho(dt) = exp(-dt/tau) -> 1 (b) Increased variance: Var = T*tau/(2-q) diverges (c) Flickering near T* The covariance matrix Sigma = T * (H_eff)^{-1} concentrates variance into the first principal component as the smallest H_eff eigenvalue vanishes (proof sketch in early_warning_signals, Paper 3). Lean theorem: early_warning_signals in CESProofs/Dynamics/FluctuationResponse.lean Subtests: 1. US variance-stress correspondence: Spearman(explained_var, NFCI) > 0. Block bootstrap corrects for serial correlation from rolling window. 2. US recession vs expansion contrast (episode-level): explained_var is higher during NBER recession episodes. Block permutation test. 3. US pre-recession eigenvalue trend: explained_var increases in the 12 months before NBER peaks. Permutation test. 4. Cross-country GFC concentration surge: For each of ~35 OECD countries, sector correlation matrix eigenvalue concentration increases from pre-GFC (2001-2007) to GFC (2008-2013). Wilcoxon signed-rank test. This provides N=35 independent paired observations, overcoming the power limitation of the US-only subtests. 5. Cross-country dot-com placebo: Same method as subtest 4 but comparing pre-dot-com (1995-1998) to dot-com (1999-2002). The 2001 recession was sector-specific (tech), not a broad information friction increase, so eigenvalue concentration should NOT increase. CONSISTENT if p > 0.10. Data: FRED IP subsectors + NFCI + USREC (subtests 1-3); OECD STAN deflated VA by industry (subtests 4-5). Review Notes (2026-03-11): Complete rewrite from threshold-based test (lambda_2/lambda_1 < 0.3 in 6-18 months before peaks, which fired 1/6 recessions, p=0.918). The 36-month rolling correlation window creates extreme serial correlation (AC1=0.985, N_eff ~ 18 vs N_raw=612). Naive p-values inflated ~35x. Fixed with block bootstrap (Subtest 1) and episode-level block permutation (Subtest 2). US-only subtests are AMBIGUOUS (correct direction, insufficient power with 7 recessions). Subtest 4 (cross-country GFC) exploits the 2008 financial crisis hitting all OECD countries, providing 35 independent paired observations. The Wilcoxon p = 0.011 is CONSISTENT. The 2001 dot-com bust shows no signal (p=0.754) because it was sector-specific (tech), not a broad information friction increase — correctly discriminating between crisis types. """ import warnings warnings.filterwarnings("ignore") import numpy as np import pandas as pd import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.data import load_fred, load_fred_multi from common.stats import log_result IP_SUBSECTORS = [ 'IPG331S', 'IPG332S', 'IPG333S', 'IPG334S', 'IPG335S', 'IPG322S', 'IPG324S', 'IPG325S', 'IPG326S', 'IPG327S', 'IPG336S', 'IPG337S', 'IPG339S', ] NBER_PEAKS = ['1973-11', '1980-01', '1981-07', '1990-07', '2001-03', '2007-12', '2020-02'] ROLLING_WINDOW = 36 # months # OECD STAN non-overlapping ISIC Rev.4 sectors (shared with test_estimation_paradox) OECD_SECTORS = [ 'A', 'B', 'C10T12', 'C13T15', 'C16T18', 'C19', 'C20_21', 'C22_23', 'C24_25', 'C26T28', 'C29_30', 'C31T33', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M_N', 'O', 'P', 'Q' ] def rolling_eigen(df, window=36): """Compute rolling eigenvalues of correlation matrix.""" results = [] for i in range(window, len(df)): sub = df.iloc[i - window:i] sub_clean = sub.dropna(axis=1) if sub_clean.shape[1] < 3: continue corr = sub_clean.corr() eigvals = np.linalg.eigvalsh(corr.values) eigvals = np.sort(eigvals)[::-1] # descending results.append({ 'date': df.index[i], 'lambda1': eigvals[0], 'lambda2': eigvals[1] if len(eigvals) > 1 else np.nan, 'ratio': eigvals[1] / eigvals[0] if len(eigvals) > 1 and eigvals[0] > 0 else np.nan, 'explained_var': eigvals[0] / eigvals.sum() if eigvals.sum() > 0 else np.nan, }) return pd.DataFrame(results).set_index('date') def _explained_var_safe(growth_df): """Compute explained variance from correlation matrix, handling degeneracies.""" corr = growth_df.corr() m = corr.values if np.any(np.isnan(m)) or np.any(np.isinf(m)): return np.nan try: eigvals = np.linalg.eigvalsh(m) s = eigvals.sum() return float(max(eigvals) / s) if s > 0 else np.nan except np.linalg.LinAlgError: return np.nan def _crosscountry_eigenvalue_test(pivot, pre_range, crisis_range, min_sectors=12, min_years=4): """Compute per-country eigenvalue concentration in two periods. Returns dict mapping country -> {ev_pre, ev_crisis, delta, n_sectors} or empty dict if insufficient data. """ country_results = {} for c in sorted(pivot.index.get_level_values('country').unique()): cdata = pivot.loc[c].sort_index() avail = [s for s in OECD_SECTORS if s in cdata.columns and cdata[s].notna().sum() >= 15] if len(avail) < min_sectors: continue growth = cdata[avail].pct_change().iloc[1:] pre = growth.loc[pre_range[0]:pre_range[1]].dropna(axis=1, how='all') crisis = growth.loc[crisis_range[0]:crisis_range[1]].dropna(axis=1, how='all') common = sorted(set(pre.columns) & set(crisis.columns)) if len(common) < min_sectors: continue pre = pre[common].dropna() crisis = crisis[common].dropna() if len(pre) < min_years or len(crisis) < min_years: continue ev_pre = _explained_var_safe(pre) ev_crisis = _explained_var_safe(crisis) if np.isnan(ev_pre) or np.isnan(ev_crisis): continue country_results[c] = { 'ev_pre': ev_pre, 'ev_crisis': ev_crisis, 'delta': ev_crisis - ev_pre, 'n_sectors': len(common), } return country_results _cache = {} def main(): results = [] # ========================================================= # Load FRED IP data for US subtests (1-3) # ========================================================= try: ip_df = load_fred_multi(IP_SUBSECTORS) except FileNotFoundError: ip_df = None eigen_df = None if ip_df is not None and ip_df.shape[1] >= 4: log_ret = np.log(ip_df / ip_df.shift(1)).dropna() eigen_df = rolling_eigen(log_ret, window=ROLLING_WINDOW) if len(eigen_df) < 60: eigen_df = None # ========================================================= # Subtest 1: US Variance-stress correspondence # ========================================================= if eigen_df is not None: try: nfci = load_fred('NFCI') nfci_monthly = nfci.resample('MS').mean() except (FileNotFoundError, Exception): nfci_monthly = None if nfci_monthly is not None and len(nfci_monthly) > 0: from scipy.stats import spearmanr merged = eigen_df.join(nfci_monthly.rename('nfci'), how='inner').dropna() if len(merged) >= 30: rho_sp, p_naive = spearmanr(merged['explained_var'], merged['nfci']) # Block bootstrap for proper p-value rng = np.random.default_rng(42) n_boot = 5000 block_size = ROLLING_WINDOW n_blocks = len(merged) // block_size boot_rhos = [] for _ in range(n_boot): idxs = rng.integers(0, n_blocks, size=n_blocks) bi = np.concatenate([ np.arange(j * block_size, min((j + 1) * block_size, len(merged))) for j in idxs ]) bi = bi[bi < len(merged)] b = merged.iloc[bi] br, _ = spearmanr(b['explained_var'], b['nfci']) boot_rhos.append(br) boot_rhos = np.array(boot_rhos) p_block = float(np.mean(boot_rhos <= 0)) ci_lo = float(np.percentile(boot_rhos, 2.5)) ci_hi = float(np.percentile(boot_rhos, 97.5)) sign_ok = rho_sp > 0 if sign_ok and p_block < 0.05: v1 = 'CONSISTENT' elif sign_ok and p_block < 0.10: v1 = 'AMBIGUOUS' elif not sign_ok and p_block < 0.05: v1 = 'INCONSISTENT' else: v1 = 'AMBIGUOUS' ev = merged['explained_var'] acfs = [float(ev.autocorr(k)) for k in range(1, 2 * ROLLING_WINDOW)] n_eff_factor = 1 + 2 * sum(max(a, 0) for a in acfs) n_eff = int(len(merged) / n_eff_factor) results.append(log_result( 'eigenvalue_early_warning_stress', stat=float(rho_sp), p=float(p_block), verd=v1, details={ 'spearman_rho': round(float(rho_sp), 3), 'p_naive': float(p_naive), 'p_block_bootstrap': round(float(p_block), 4), 'boot_95ci': [round(ci_lo, 3), round(ci_hi, 3)], 'n_obs': len(merged), 'n_eff': n_eff, }, description=( f'US variance-stress: Spearman(explained_var, NFCI) = ' f'{rho_sp:.3f}. Block bootstrap p = {p_block:.4f} ' f'(naive p = {p_naive:.2e}; N_eff = {n_eff} vs N_raw = {len(merged)}). ' f'95% CI: [{ci_lo:.3f}, {ci_hi:.3f}]. ' f'Theory: rising T concentrates variance in PC1. ' f'CONSISTENT if rho > 0 and p < 0.05.' ), )) else: results.append(log_result( 'eigenvalue_early_warning_stress', stat=np.nan, p=np.nan, verd='SKIPPED', description=f'Insufficient NFCI overlap ({len(merged)} obs).', )) else: results.append(log_result( 'eigenvalue_early_warning_stress', stat=np.nan, p=np.nan, verd='SKIPPED', description='NFCI data not available.', )) else: results.append(log_result( 'eigenvalue_early_warning_stress', stat=np.nan, p=np.nan, verd='SKIPPED', description='FRED IP subsector data not available.', )) # ========================================================= # Subtest 2: US Recession vs expansion contrast (episode-level) # ========================================================= if eigen_df is not None: try: usrec = load_fred('USREC') usrec_monthly = usrec.resample('MS').last() except (FileNotFoundError, Exception): usrec_monthly = None if usrec_monthly is not None and len(usrec_monthly) > 0: merged_rec = eigen_df.join(usrec_monthly.rename('usrec'), how='inner').dropna() rec_months = merged_rec[merged_rec['usrec'] == 1].copy() exp_months = merged_rec[merged_rec['usrec'] == 0] if len(rec_months) >= 6 and len(exp_months) >= 30: date_diffs = rec_months.index.to_series().diff() episode_breaks = date_diffs > pd.Timedelta(days=120) rec_months['episode'] = episode_breaks.cumsum().values episode_means = rec_months.groupby('episode')['explained_var'].mean() n_episodes = len(episode_means) exp_mean = float(exp_months['explained_var'].mean()) rec_ep_mean = float(episode_means.mean()) obs_diff = rec_ep_mean - exp_mean n_above = int(sum(episode_means > exp_mean)) rng = np.random.default_rng(42) rec_lengths = rec_months.groupby('episode').size().values all_ev = merged_rec['explained_var'].values n_total = len(all_ev) perm_diffs = [] for _ in range(5000): fake_vals = [] for length in rec_lengths: start = rng.integers(0, max(1, n_total - length)) fake_vals.append(float(np.mean(all_ev[start:start + length]))) perm_diffs.append(np.mean(fake_vals) - np.mean(all_ev)) perm_diffs = np.array(perm_diffs) p_ep = float(np.mean(perm_diffs >= obs_diff)) sign_ok = obs_diff > 0 if sign_ok and p_ep < 0.05: v2 = 'CONSISTENT' elif sign_ok and p_ep < 0.10: v2 = 'AMBIGUOUS' elif not sign_ok and p_ep < 0.05: v2 = 'INCONSISTENT' else: v2 = 'AMBIGUOUS' results.append(log_result( 'eigenvalue_early_warning_recession', stat=float(obs_diff), p=float(p_ep), verd=v2, details={ 'rec_episode_mean': round(rec_ep_mean, 3), 'exp_mean': round(exp_mean, 3), 'diff': round(obs_diff, 3), 'n_episodes': n_episodes, 'n_above_expansion': n_above, }, description=( f'US recession episodes: mean explained_var = ' f'{rec_ep_mean:.3f} ({n_episodes} episodes, {n_above} above ' f'expansion mean {exp_mean:.3f}). ' f'Diff = {obs_diff:+.3f}. Block permutation p = {p_ep:.4f}. ' f'CONSISTENT if rec > exp and p < 0.05.' ), )) else: results.append(log_result( 'eigenvalue_early_warning_recession', stat=np.nan, p=np.nan, verd='SKIPPED', description='Insufficient US recession data.', )) else: results.append(log_result( 'eigenvalue_early_warning_recession', stat=np.nan, p=np.nan, verd='SKIPPED', description='USREC data not available.', )) else: results.append(log_result( 'eigenvalue_early_warning_recession', stat=np.nan, p=np.nan, verd='SKIPPED', description='FRED IP data not available.', )) # ========================================================= # Subtest 3: US Pre-recession eigenvalue trend # ========================================================= if eigen_df is not None: from scipy.stats import spearmanr as sp_corr testable_peaks = [] trend_correct = [] for peak_str in NBER_PEAKS: peak_date = pd.Timestamp(peak_str) window_start = peak_date - pd.DateOffset(months=12) if window_start < eigen_df.index.min(): continue if peak_date > eigen_df.index.max(): continue mask = (eigen_df.index >= window_start) & (eigen_df.index <= peak_date) w = eigen_df.loc[mask, 'explained_var'].dropna() if len(w) < 6: continue testable_peaks.append(peak_str) r_trend, _ = sp_corr(np.arange(len(w)), w.values) trend_correct.append(r_trend > 0) n_testable = len(testable_peaks) n_correct = sum(trend_correct) if n_testable > 0: fraction = n_correct / n_testable rng = np.random.default_rng(42) valid_dates = eigen_df.index[ (eigen_df.index >= eigen_df.index.min() + pd.DateOffset(months=12)) & (eigen_df.index <= eigen_df.index.max()) ] perm_fractions = [] for _ in range(5000): fake_peaks = rng.choice(valid_dates, size=n_testable, replace=False) fake_correct = 0 for fp in fake_peaks: w_start = fp - pd.DateOffset(months=12) mask = (eigen_df.index >= w_start) & (eigen_df.index <= fp) w = eigen_df.loc[mask, 'explained_var'].dropna() if len(w) >= 6: r_t, _ = sp_corr(np.arange(len(w)), w.values) if r_t > 0: fake_correct += 1 perm_fractions.append(fake_correct / n_testable) perm_fractions = np.array(perm_fractions) p_perm = float(np.mean(perm_fractions >= fraction)) sign_ok = fraction >= 0.5 if sign_ok and p_perm < 0.10: v3 = 'CONSISTENT' elif sign_ok: v3 = 'AMBIGUOUS' else: v3 = 'INCONSISTENT' results.append(log_result( 'eigenvalue_early_warning_trend', stat=float(fraction), p=float(p_perm), verd=v3, details={ 'n_testable': n_testable, 'n_correct_trend': n_correct, 'testable_peaks': testable_peaks, 'perm_mean': round(float(np.mean(perm_fractions)), 3), }, description=( f'US pre-recession trend: explained_var increases in 12 months ' f'before {n_correct}/{n_testable} NBER peaks. ' f'Permutation p = {p_perm:.3f}. ' f'CONSISTENT if >= 50% correct and p < 0.10.' ), )) else: results.append(log_result( 'eigenvalue_early_warning_trend', stat=np.nan, p=np.nan, verd='SKIPPED', description='No testable US recession peaks.', )) else: results.append(log_result( 'eigenvalue_early_warning_trend', stat=np.nan, p=np.nan, verd='SKIPPED', description='FRED IP data not available.', )) # ========================================================= # Subtests 4-5: Cross-country concentration surge tests # ========================================================= try: from common.data import load_oecd_stan stan = load_oecd_stan() except (FileNotFoundError, ImportError, Exception): stan = None pivot = None if stan is not None and len(stan) > 0: from scipy.stats import wilcoxon defl = stan[stan['price_base'] == 'D'] sector_data = defl[defl['activity'].isin(OECD_SECTORS)] pivot = sector_data.pivot_table( index=['country', 'year'], columns='activity', values='value', ) # --- Subtest 4: GFC surge (broad financial crisis → T↑ globally) --- if pivot is not None: gfc_results = _crosscountry_eigenvalue_test(pivot, (2001, 2007), (2008, 2013)) n_countries = len(gfc_results) if n_countries >= 15: pre_vals = np.array([r['ev_pre'] for r in gfc_results.values()]) gfc_vals = np.array([r['ev_crisis'] for r in gfc_results.values()]) deltas = gfc_vals - pre_vals n_up = int(np.sum(deltas > 0)) stat_w, p_w = wilcoxon(gfc_vals, pre_vals, alternative='greater') mean_delta = float(np.mean(deltas)) sign_ok = mean_delta > 0 if sign_ok and p_w < 0.05: v4 = 'CONSISTENT' elif sign_ok and p_w < 0.10: v4 = 'AMBIGUOUS' elif not sign_ok and p_w < 0.05: v4 = 'INCONSISTENT' else: v4 = 'AMBIGUOUS' _cache['gfc_country_results'] = gfc_results results.append(log_result( 'eigenvalue_early_warning_gfc_crosscountry', stat=float(mean_delta), p=float(p_w), verd=v4, details={ 'n_countries': n_countries, 'n_concentration_increased': n_up, 'mean_ev_pre': round(float(np.mean(pre_vals)), 3), 'mean_ev_gfc': round(float(np.mean(gfc_vals)), 3), 'mean_delta': round(mean_delta, 3), 'pre_range': [2001, 2007], 'gfc_range': [2008, 2013], 'top5_delta': [ {'country': c, 'delta': round(r['delta'], 3)} for c, r in sorted(gfc_results.items(), key=lambda x: -x[1]['delta'])[:5] ], }, description=( f'Cross-country GFC surge: sector eigenvalue concentration ' f'increases from pre-GFC (2001-2007) to ' f'GFC (2008-2013) in {n_up}/{n_countries} ' f'OECD countries. Mean delta = {mean_delta:+.3f}. ' f'Wilcoxon signed-rank p = {p_w:.4f}. ' f'Theory: the 2008 financial crisis raised information friction T ' f'globally, concentrating sector variance into PC1. ' f'Each country provides one independent paired observation, ' f'overcoming the power limitation of the US-only subtests. ' f'CONSISTENT if delta > 0 and p < 0.05.' ), )) else: results.append(log_result( 'eigenvalue_early_warning_gfc_crosscountry', stat=np.nan, p=np.nan, verd='SKIPPED', details={'n_countries_found': n_countries}, description=f'Insufficient OECD countries with sector data ({n_countries}).', )) else: results.append(log_result( 'eigenvalue_early_warning_gfc_crosscountry', stat=np.nan, p=np.nan, verd='SKIPPED', description='OECD STAN data not available.', )) # --- Subtest 5: Dot-com placebo (sector-specific tech bust → no broad T↑) --- if pivot is not None: dotcom_results = _crosscountry_eigenvalue_test(pivot, (1995, 1998), (1999, 2002)) n_dc = len(dotcom_results) if n_dc >= 15: pre_dc = np.array([r['ev_pre'] for r in dotcom_results.values()]) crisis_dc = np.array([r['ev_crisis'] for r in dotcom_results.values()]) deltas_dc = crisis_dc - pre_dc n_up_dc = int(np.sum(deltas_dc > 0)) # Two-sided test: we're testing the null (no effect expected) stat_dc, p_dc_two = wilcoxon(crisis_dc, pre_dc) # Also one-sided for comparison _, p_dc_greater = wilcoxon(crisis_dc, pre_dc, alternative='greater') mean_delta_dc = float(np.mean(deltas_dc)) # CONSISTENT if NOT significant (placebo should show no effect) if p_dc_greater > 0.10: v5 = 'CONSISTENT' elif p_dc_greater > 0.05: v5 = 'AMBIGUOUS' else: v5 = 'INCONSISTENT' _cache['dotcom_country_results'] = dotcom_results results.append(log_result( 'eigenvalue_early_warning_dotcom_placebo', stat=float(mean_delta_dc), p=float(p_dc_greater), verd=v5, details={ 'n_countries': n_dc, 'n_concentration_increased': n_up_dc, 'mean_ev_pre': round(float(np.mean(pre_dc)), 3), 'mean_ev_dotcom': round(float(np.mean(crisis_dc)), 3), 'mean_delta': round(mean_delta_dc, 3), 'pre_range': [1995, 1998], 'dotcom_range': [1999, 2002], 'p_two_sided': round(float(p_dc_two), 4), }, description=( f'Dot-com placebo: sector eigenvalue concentration from ' f'pre-dot-com (1995-1998) to dot-com (1999-2002) in ' f'{n_up_dc}/{n_dc} OECD countries. ' f'Mean delta = {mean_delta_dc:+.3f}. ' f'Wilcoxon one-sided p = {p_dc_greater:.4f}. ' f'Theory: the 2001 dot-com bust was sector-specific (tech), ' f'not a broad information friction increase. Unlike the GFC, ' f'it should NOT concentrate sector variance into PC1. ' f'Discriminates between broad crises (GFC: T↑ everywhere) ' f'and sector-specific shocks (dot-com: T↑ in tech only). ' f'CONSISTENT if p > 0.10 (no significant concentration increase).' ), )) else: results.append(log_result( 'eigenvalue_early_warning_dotcom_placebo', stat=np.nan, p=np.nan, verd='SKIPPED', details={'n_countries_found': n_dc}, description=f'Insufficient OECD countries for dot-com period ({n_dc}).', )) else: results.append(log_result( 'eigenvalue_early_warning_dotcom_placebo', stat=np.nan, p=np.nan, verd='SKIPPED', description='OECD STAN data not available.', )) _cache['eigen_df'] = eigen_df _cache['results'] = results return results def figures(): """Generate figures for eigenvalue early warning subtests.""" if not _cache: main() paths = [] results = _cache.get('results', []) eigen_df = _cache.get('eigen_df') from common.plotting import setup_figure, setup_subplots, save_figure # --- Figure 1: Stress correspondence (subtest 1) --- r1 = next((r for r in results if r['name'] == 'eigenvalue_early_warning_stress'), None) if r1 and r1['verdict'] != 'SKIPPED' and eigen_df is not None: try: nfci = load_fred('NFCI') nfci_monthly = nfci.resample('MS').mean() merged = eigen_df.join(nfci_monthly.rename('nfci'), how='inner').dropna() if len(merged) >= 30: fig, axes = setup_subplots(1, 2, width=14, height=5) # Panel 1: time series ax = axes[0] ax.plot(merged.index, merged['explained_var'], 'b-', alpha=0.7, label='Explained var (PC1)') ax2 = ax.twinx() ax2.plot(merged.index, merged['nfci'], 'r-', alpha=0.5, label='NFCI') ax2.set_ylabel('NFCI', color='red') ax.set_ylabel('Explained Variance (PC1)') ax.set_title('Eigenvalue Concentration vs Financial Stress') ax.legend(loc='upper left', fontsize=8) ax2.legend(loc='upper right', fontsize=8) # Panel 2: scatter ax = axes[1] ax.scatter(merged['nfci'], merged['explained_var'], alpha=0.3, s=10) rho_sp = r1.get('spearman_rho', r1['statistic']) p_block = r1.get('p_block_bootstrap', r1['p_value']) ax.set_xlabel('NFCI (financial stress)') ax.set_ylabel('Explained Variance (PC1)') ax.set_title(f'Spearman r={rho_sp:.3f}, block p={p_block:.3f}') ax.grid(True, alpha=0.3) fig.tight_layout() paths += save_figure(fig, 'eigenvalue_ew_stress', subfolder='potential') except Exception: pass # --- Figure 2: Recession vs expansion (subtest 2) --- r2 = next((r for r in results if r['name'] == 'eigenvalue_early_warning_recession'), None) if r2 and r2['verdict'] != 'SKIPPED': try: rec_mean = r2.get('rec_episode_mean', 0) exp_mean = r2.get('exp_mean', 0) n_ep = r2.get('n_episodes', 0) fig, ax = setup_figure(7, 5) bars = ax.bar(['Expansion', 'Recession'], [exp_mean, rec_mean], color=['steelblue', 'firebrick'], alpha=0.8) ax.set_ylabel('Mean Explained Variance (PC1)') diff = r2.get('diff', rec_mean - exp_mean) p_ep = r2['p_value'] ax.set_title(f'Recession vs Expansion: diff={diff:+.3f}, ' f'block perm p={p_ep:.3f} (N={n_ep} episodes)') ax.grid(True, alpha=0.3, axis='y') paths += save_figure(fig, 'eigenvalue_ew_recession', subfolder='potential') except Exception: pass # --- Figure 3: Pre-recession trend (subtest 3) --- r3 = next((r for r in results if r['name'] == 'eigenvalue_early_warning_trend'), None) if r3 and r3['verdict'] != 'SKIPPED' and eigen_df is not None: from scipy.stats import spearmanr peaks = r3.get('testable_peaks', []) if peaks: n_panels = min(len(peaks), 6) ncols = min(3, n_panels) nrows = (n_panels + ncols - 1) // ncols fig, axes = setup_subplots(nrows, ncols, width=4*ncols, height=3.5*nrows) axes_flat = axes.flat if hasattr(axes, 'flat') else [axes] for i, peak_str in enumerate(peaks[:6]): ax = axes_flat[i] peak_date = pd.Timestamp(peak_str) w_start = peak_date - pd.DateOffset(months=12) mask = (eigen_df.index >= w_start) & (eigen_df.index <= peak_date) w = eigen_df.loc[mask, 'explained_var'].dropna() if len(w) >= 6: t = np.arange(len(w)) ax.plot(t, w.values, 'b-o', markersize=3) r_t, _ = spearmanr(t, w.values) color = 'green' if r_t > 0 else 'red' ax.set_title(f'{peak_str}: r={r_t:.2f}', fontsize=9, color=color) ax.set_xlabel('Months before peak', fontsize=8) ax.set_ylabel('Expl. Var', fontsize=8) # Hide unused axes for j in range(n_panels, len(axes_flat)): axes_flat[j].set_visible(False) fig.suptitle(f'Pre-Recession Trend: {r3.get("n_correct_trend",0)}/{len(peaks)} ' f'increasing (p={r3["p_value"]:.3f})', fontsize=11) fig.tight_layout() paths += save_figure(fig, 'eigenvalue_ew_trend', subfolder='potential') # --- Figure 4: Cross-country GFC (subtest 4) --- r4 = next((r for r in results if r['name'] == 'eigenvalue_early_warning_gfc_crosscountry'), None) if r4 and r4['verdict'] != 'SKIPPED': n_c = r4.get('n_countries', 0) mean_pre = r4.get('mean_ev_pre', 0) mean_gfc = r4.get('mean_ev_gfc', 0) n_up = r4.get('n_concentration_increased', 0) top5 = r4.get('top5_delta', []) fig, axes = setup_subplots(1, 2, width=12, height=5) # Panel 1: paired comparison ax = axes[0] ax.bar(['Pre-GFC\n(2001-07)', 'GFC\n(2008-13)'], [mean_pre, mean_gfc], color=['steelblue', 'firebrick'], alpha=0.8) ax.set_ylabel('Mean Explained Variance (PC1)') ax.set_title(f'Cross-Country: {n_up}/{n_c} countries increased\n' f'Wilcoxon p={r4["p_value"]:.4f}') ax.grid(True, alpha=0.3, axis='y') # Panel 2: top delta countries if top5: countries = [d['country'] for d in top5] deltas = [d['delta'] for d in top5] ax = axes[1] ax.barh(countries[::-1], deltas[::-1], color='firebrick', alpha=0.7) ax.set_xlabel('Delta (GFC - Pre)') ax.set_title('Top 5 Countries by Concentration Increase') ax.grid(True, alpha=0.3, axis='x') fig.tight_layout() paths += save_figure(fig, 'eigenvalue_ew_gfc', subfolder='potential') # --- Figure 5: GFC vs dot-com comparison (subtests 4+5) --- r5 = next((r for r in results if r['name'] == 'eigenvalue_early_warning_dotcom_placebo'), None) if r4 and r4['verdict'] != 'SKIPPED' and r5 and r5['verdict'] != 'SKIPPED': fig, axes = setup_subplots(1, 2, width=12, height=5) # Panel 1: GFC delta distribution gfc_cr = _cache.get('gfc_country_results', {}) dc_cr = _cache.get('dotcom_country_results', {}) if gfc_cr and dc_cr: gfc_deltas = [r['delta'] for r in gfc_cr.values()] dc_deltas = [r['delta'] for r in dc_cr.values()] ax = axes[0] positions = [1, 2] bp = ax.boxplot([gfc_deltas, dc_deltas], positions=positions, widths=0.6, patch_artist=True) bp['boxes'][0].set_facecolor('firebrick') bp['boxes'][0].set_alpha(0.6) bp['boxes'][1].set_facecolor('steelblue') bp['boxes'][1].set_alpha(0.6) ax.axhline(0, color='black', linewidth=0.5, linestyle='--') ax.set_xticks(positions) ax.set_xticklabels(['GFC\n(2008-13)', 'Dot-com\n(1999-02)']) ax.set_ylabel('Delta Explained Variance') ax.set_title('Broad Crisis vs Sector-Specific Shock') ax.grid(True, alpha=0.3, axis='y') # Panel 2: side-by-side bar comparison ax = axes[1] gfc_p = r4['p_value'] dc_p = r5['p_value'] gfc_d = r4.get('mean_delta', 0) dc_d = r5.get('mean_delta', 0) bars = ax.bar(['GFC', 'Dot-com'], [gfc_d, dc_d], color=['firebrick', 'steelblue'], alpha=0.8) ax.axhline(0, color='black', linewidth=0.5, linestyle='--') ax.set_ylabel('Mean Delta (crisis - pre)') ax.set_title(f'GFC p={gfc_p:.4f} (broad T↑)\n' f'Dot-com p={dc_p:.3f} (sector-specific)') ax.grid(True, alpha=0.3, axis='y') fig.tight_layout() paths += save_figure(fig, 'eigenvalue_ew_dotcom_placebo', subfolder='potential') return paths if __name__ == '__main__': res = main() for r in res: print(f" {r['name']}: {r['verdict']} " f"(stat={r['statistic']}, p={r.get('p_value', 'n/a')})") if r.get('details'): for k, v in r['details'].items(): print(f" {k}: {v}")

diff	0.0380
p value	0.2338
exp mean	0.3390
statistic	0.0378
n episodes	7
rec episode mean	0.3770
n above expansion	4

Eigenvalue Early Warning Recession

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps