#!/usr/bin/env python3 """ Test: EMD Timescale Hierarchy — Paper 5 (Empirical Tests) Blind discovery of economic timescale hierarchy using Ensemble EMD. Tests six predictions: P1: N_eff = 4-5 distinct timescale bands P2: Adjacent timescale ratios r* ~ 2 P3: Periods match classical business cycles P4: Sector N_eff negatively correlated with sigma P5: Adjacent IMF energy > non-adjacent (nearest-neighbor coupling) P6: Contraction amplitude > expansion amplitude (relaxation asymmetry) Data: FRED INDPRO, manufacturing subsectors, USREC. Lean backing: - Paper 4, Theorem thm:topology(iv): nearest-neighbor coupling topology is exact as r* → ∞, O(1/r*) approximation at finite r*. - Paper 4, Section on timescale separation: geometric progression r* ~ 2. Review Notes (2026-03-12): Theory match: CORRECT for P1-P3, P5-P6. PARTIAL for P4. P1-P3: EEMD is a blind, data-driven decomposition that discovers the timescale hierarchy predicted by CES theory. N_eff=5, r*≈2.3, and 4 classical cycle matches are strong results. P4: Only 1 of 7 sectors has different N_eff (Computer/Electronics=6, all others=5). Kendall tau = -0.534 is driven by a single observation. P5: Theory predicts adjacent coupling > non-adjacent, with non-adjacent suppressed by O(1/r*). At r*≈2, non-adjacent coupling is ~50% of adjacent. Test shows adj_mean=0.208 vs nonadj_mean=-0.044. P6: Business-cycle IMFs have higher amplitude during recessions (ratio=1.30), consistent with faster contraction dynamics on the CES landscape. I/O audit: CORRECT. All variables are outputs (IP growth decomposition). Data-signal alignment: STRONG. FRED INDPRO is 100+ years monthly. Code quality: MINOR ISSUES (fixed 2026-03-12). - P4 verdict ignored p-value (CONSISTENT on sign alone, p=0.13) - P5 had no p-value (Mann-Whitney test added) - EEMD is stochastic: different runs produce slightly different modes. IMF 4 (6.25yr) is borderline at energy=0.018 vs threshold 0.02. Confidence in verdict: MEDIUM-HIGH. P1-P3 are the strongest results. P4 is one-point driven. P5-P6 lack p-values for rigorous significance testing. """ 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 from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure try: from scipy import stats from scipy.signal import hilbert HAS_SCIPY = True except ImportError: HAS_SCIPY = False try: from PyEMD import EEMD HAS_EMD = True except ImportError: HAS_EMD = False EEMD_TRIALS = 200 EEMD_NOISE_WIDTH = 0.2 ENERGY_THRESHOLD = 0.02 SECTORS = { 'IPG334S': ('Computer/Electronics', 0.40), 'IPG335S': ('Electrical Equipment', 0.55), 'IPG333S': ('Machinery', 0.65), 'IPG325S': ('Chemicals', 0.70), 'IPG336S': ('Transport Equipment', 0.75), 'IPG331S': ('Primary Metals', 0.85), 'IPG311A2S': ('Food/Beverage', 0.90), } CLASSICAL_CYCLES = { 'Kitchin': (2, 4), 'Business': (4, 8), 'Juglar': (8, 12), 'Kuznets': (15, 25), 'Kondratiev': (40, 60), } def _to_growth(series): vals = series.values.astype(float) vals = vals[vals > 0] growth = np.diff(np.log(vals)) idx = series.index[1:len(growth)+1] return pd.Series(growth, index=idx, name=series.name) def _run_eemd(signal, trials=EEMD_TRIALS, noise_width=EEMD_NOISE_WIDTH): eemd = EEMD(trials=trials, noise_width=noise_width) eemd.noise_seed(42) return eemd.eemd(signal) def _characterize_imfs(imfs, dt_months=1.0): total_energy = sum(np.var(imf) for imf in imfs) if total_energy == 0: return [] results = [] for k, imf in enumerate(imfs): energy_frac = np.var(imf) / total_energy analytic = hilbert(imf) inst_phase = np.unwrap(np.angle(analytic)) inst_freq = np.diff(inst_phase) / (2.0 * np.pi * dt_months) inst_freq = inst_freq[inst_freq > 0] if len(inst_freq) > 10: inst_period = 1.0 / inst_freq max_period = len(imf) * dt_months inst_period = inst_period[(inst_period > 0) & (inst_period < max_period)] char_period = np.median(inst_period) if len(inst_period) > 0 else np.nan else: char_period = np.nan results.append({ 'imf_index': k, 'energy_fraction': energy_frac, 'char_period_months': char_period, 'char_period_years': char_period / 12.0 if not np.isnan(char_period) else np.nan, 'significant': energy_frac > ENERGY_THRESHOLD, }) return results _cache = {} def main(): results = [] np.random.seed(42) if not HAS_SCIPY or not HAS_EMD: results.append(log_result( "emd_hierarchy", np.nan, np.nan, "AMBIGUOUS", details={"reason": "PyEMD or scipy not installed"}, description="Skipped: PyEMD and scipy required for Ensemble Empirical Mode Decomposition." )) return results # Load INDPRO try: indpro_raw = load_fred('INDPRO') except FileNotFoundError: results.append(log_result( "emd_hierarchy", np.nan, np.nan, "AMBIGUOUS", details={"reason": "INDPRO not found"}, description="Skipped: FRED INDPRO (industrial production index) data not available for EMD decomposition." )) return results indpro_growth = _to_growth(indpro_raw) signal = indpro_growth.values.astype(float) # Run EEMD imfs = _run_eemd(signal) imf_info = _characterize_imfs(imfs) n_eff = sum(1 for i in imf_info if i['significant']) sig_periods = sorted([i['char_period_years'] for i in imf_info if i['significant'] and not np.isnan(i['char_period_years'])]) # P1: N_eff count p1_pass = 3 <= n_eff <= 7 results.append(log_result( "emd_hierarchy_neff_count", stat=float(n_eff), p=np.nan, verd="CONSISTENT" if p1_pass else "AMBIGUOUS", description="Tests the CES hierarchy prediction of N_eff = 4-5 distinct timescale bands. " "Method: count EEMD modes with energy fraction > 2% from INDPRO monthly " "log-growth (1919-present). Consistent if 3 <= N_eff <= 7." )) # P2: r* ~ 2 ratios = [] r_star = np.nan max_ratio = np.nan if len(sig_periods) >= 2: ratios = [sig_periods[i+1] / sig_periods[i] for i in range(len(sig_periods)-1) if sig_periods[i] > 0] r_star = np.median(ratios) if ratios else np.nan max_ratio = max(ratios) if ratios else np.nan p2_pass = 1.5 <= r_star <= 3.5 if not np.isnan(r_star) else False # Flag broken chain: if any individual ratio exceeds 5, the hierarchy # has a gap (likely a borderline mode dropped by the energy threshold) chain_note = "" if not np.isnan(max_ratio) and max_ratio > 5.0: chain_note = (f" NOTE: max individual ratio = {max_ratio:.1f} indicates a gap " f"in the hierarchy (a borderline mode at energy ~{ENERGY_THRESHOLD:.0%} " f"may have been excluded). Median r* is robust to this outlier.") results.append(log_result( "emd_hierarchy_adjacent_ratio", stat=float(r_star) if not np.isnan(r_star) else np.nan, p=np.nan, verd="CONSISTENT" if p2_pass else "AMBIGUOUS", details={"ratios": [float(r) for r in ratios], "max_ratio": float(max_ratio)}, description="Tests geometric timescale spacing r* ~ 2. Method: compute ratio of " "adjacent significant IMF characteristic periods (Hilbert instantaneous " "frequency median), report median ratio. Data: EEMD of INDPRO log-growth. " f"Consistent if 1.5 <= median r* <= 3.5.{chain_note}" )) else: results.append(log_result( "emd_hierarchy_adjacent_ratio", np.nan, np.nan, "AMBIGUOUS", details={"reason": "< 2 significant IMFs"}, description="Skipped: fewer than 2 significant IMFs extracted, cannot compute adjacent period ratios." )) # P3: Classical cycle match all_periods = [(i['char_period_years'], i['significant']) for i in imf_info if not np.isnan(i['char_period_years']) and i['energy_fraction'] > 0.01] matched = set() for period, _ in all_periods: for cycle_name, (lo, hi) in CLASSICAL_CYCLES.items(): if lo <= period <= hi: matched.add(cycle_name) break n_matched = len(matched) results.append(log_result( "emd_hierarchy_classical_cycle_match", stat=float(n_matched), p=np.nan, verd="CONSISTENT" if n_matched >= 2 else "AMBIGUOUS", details={"matched_cycles": list(matched)}, description="Tests whether EEMD modes match classical business cycle bands without prior " "knowledge of their existence. Method: check if IMF characteristic periods " "(energy > 1%) fall within Kitchin (2-4yr), Business (4-8yr), Juglar (8-12yr), " "Kuznets (15-25yr), Kondratiev (40-60yr) ranges. Consistent if >= 2 matched." )) # P4: Cross-sector N_eff vs sigma sector_neffs = [] sector_sigmas = [] for sid, (name, sigma) in SECTORS.items(): try: s = load_fred(sid) except FileNotFoundError: continue growth = _to_growth(s).values.astype(float) if len(growth) < 120: continue try: s_imfs = _run_eemd(growth, trials=100, noise_width=0.2) s_info = _characterize_imfs(s_imfs) s_neff = sum(1 for i in s_info if i['significant']) sector_neffs.append(s_neff) sector_sigmas.append(sigma) except Exception: continue if len(sector_neffs) >= 4: tau_stat, tau_p = stats.kendalltau(sector_sigmas, sector_neffs) # Require both correct sign AND significance for CONSISTENT sign_ok = tau_stat < 0 sig_ok = tau_p < 0.10 if sign_ok and sig_ok: p4_verd = "CONSISTENT" elif sign_ok: p4_verd = "AMBIGUOUS" else: p4_verd = "INCONSISTENT" n_unique = len(set(sector_neffs)) results.append(log_result( "emd_hierarchy_sector_neff_sigma", stat=float(tau_stat), p=float(tau_p), verd=p4_verd, details={"n_unique_neff": n_unique, "neffs": sector_neffs, "sigmas": sector_sigmas}, description="Tests whether more substitutable sectors (higher sigma) have fewer " "significant timescale bands (lower N_eff). Method: Kendall tau of " "sigma_OR vs EEMD N_eff across 7 FRED manufacturing subsectors. " "H0: tau >= 0. Prediction: tau < 0. " f"LIMITATION: {7-n_unique} of 7 sectors share the same N_eff — the " f"correlation is driven by {n_unique-1 if n_unique > 1 else 0} " f"sector(s) with different N_eff." )) else: results.append(log_result( "emd_hierarchy_sector_neff_sigma", np.nan, np.nan, "AMBIGUOUS", details={"reason": f"Only {len(sector_neffs)} sectors"}, description="Skipped: insufficient FRED sector IP series available for cross-sector N_eff vs sigma test." )) # P5: Nearest-neighbor coupling sig_idx = [i['imf_index'] for i in imf_info if i['significant']] adj_corrs = [] nonadj_corrs = [] if len(sig_idx) >= 3: roll_window = 24 energies = {} for idx in sig_idx: imf = imfs[idx] energies[idx] = pd.Series(imf).rolling(roll_window).var().dropna().values for ii, idx_i in enumerate(sig_idx): for jj, idx_j in enumerate(sig_idx): if jj <= ii: continue min_len = min(len(energies[idx_i]), len(energies[idx_j])) if min_len < 30: continue r, _ = stats.spearmanr(energies[idx_i][:min_len], energies[idx_j][:min_len]) if abs(ii - jj) == 1: adj_corrs.append(r) else: nonadj_corrs.append(r) if adj_corrs and nonadj_corrs: adj_mean = np.mean(adj_corrs) nonadj_mean = np.mean(nonadj_corrs) # Mann-Whitney test: are adjacent correlations higher? mw_stat, mw_p_two = stats.mannwhitneyu( adj_corrs, nonadj_corrs, alternative='greater') results.append(log_result( "emd_hierarchy_adjacent_coupling", stat=float(adj_mean - nonadj_mean), p=float(mw_p_two), verd="CONSISTENT" if adj_mean > nonadj_mean else "AMBIGUOUS", details={"adj_mean": float(adj_mean), "nonadj_mean": float(nonadj_mean), "adj_corrs": [float(x) for x in adj_corrs], "nonadj_corrs": [float(x) for x in nonadj_corrs], "mann_whitney_U": float(mw_stat)}, description="Tests nearest-neighbor coupling (Paper 4 Theorem thm:topology(iv)): " "rolling 24-month energy correlations (Spearman) between adjacent " "EEMD modes should exceed non-adjacent pairs. Method: Mann-Whitney U " "test, H0: adj <= nonadj (one-sided). At r*~2, non-adjacent coupling " "is O(1/r*) ~ 50% of adjacent, so the gap should be positive but modest." )) else: results.append(log_result( "emd_hierarchy_adjacent_coupling", np.nan, np.nan, "AMBIGUOUS", details={"reason": "Insufficient pairs"}, description="Skipped: not enough adjacent/non-adjacent IMF pairs to compare energy coupling strength." )) else: results.append(log_result( "emd_hierarchy_adjacent_coupling", np.nan, np.nan, "AMBIGUOUS", details={"reason": f"Only {len(sig_idx)} significant IMFs"}, description="Skipped: fewer than 3 significant IMFs, cannot test adjacent vs non-adjacent coupling." )) # P6: Relaxation asymmetry try: usrec = load_fred('USREC') common_start = max(indpro_growth.index[0], usrec.index[0]) common_end = min(indpro_growth.index[-1], usrec.index[-1]) usrec_trim = usrec.loc[common_start:common_end] growth_trim = indpro_growth.loc[common_start:common_end] recession_aligned = usrec_trim.reindex(growth_trim.index, method='ffill').fillna(0) is_recession = recession_aligned.values == 1.0 offset = indpro_growth.index.get_loc(common_start) trim_len = len(growth_trim) bc_indices = [i['imf_index'] for i in imf_info if not np.isnan(i['char_period_years']) and 2.0 <= i['char_period_years'] <= 8.0 and i['energy_fraction'] > 0.01] amp_ratios = [] for idx in bc_indices: imf = imfs[idx] imf_trim = imf[offset:offset + trim_len] if len(imf_trim) != trim_len: continue n_rec = is_recession.sum() n_exp = (~is_recession).sum() if n_rec > 10 and n_exp > 10: rms_rec = np.sqrt(np.mean(imf_trim[is_recession] ** 2)) rms_exp = np.sqrt(np.mean(imf_trim[~is_recession] ** 2)) if rms_exp > 0: amp_ratios.append(rms_rec / rms_exp) if amp_ratios: mean_ratio = np.mean(amp_ratios) results.append(log_result( "emd_hierarchy_relaxation_asymmetry", stat=float(mean_ratio), p=np.nan, verd="CONSISTENT" if mean_ratio > 1.0 else "AMBIGUOUS", details={"amp_ratios": [float(r) for r in amp_ratios], "bc_imf_indices": bc_indices}, description="Tests relaxation asymmetry: RMS amplitude of business-cycle EEMD " "modes (2-8yr period) should be larger during NBER recessions than " "expansions, reflecting faster contraction dynamics on the CES " "potential landscape. Data: INDPRO EEMD modes aligned with USREC. " "Consistent if recession/expansion amplitude ratio > 1." )) else: results.append(log_result( "emd_hierarchy_relaxation_asymmetry", np.nan, np.nan, "AMBIGUOUS", details={"reason": "No business-cycle IMFs found"}, description="Skipped: no EEMD modes with 2-8 year characteristic period found for asymmetry test." )) except FileNotFoundError: results.append(log_result( "emd_hierarchy_relaxation_asymmetry", np.nan, np.nan, "AMBIGUOUS", details={"reason": "USREC not found"}, description="Skipped: FRED USREC recession indicator not available for amplitude asymmetry test." )) # Populate cache for figures() _cache['imfs'] = imfs _cache['imf_info'] = imf_info _cache['sig_periods'] = sig_periods _cache['n_eff'] = n_eff _cache['indpro_growth'] = indpro_growth _cache['ratios'] = ratios _cache['r_star'] = r_star _cache['max_ratio'] = max_ratio _cache['sig_idx'] = sig_idx _cache['sector_neffs'] = sector_neffs _cache['sector_sigmas'] = sector_sigmas _cache['sector_names'] = [name for sid, (name, sigma) in SECTORS.items()][:len(sector_neffs)] _cache['adj_corrs'] = adj_corrs _cache['nonadj_corrs'] = nonadj_corrs _cache['classical_cycles'] = CLASSICAL_CYCLES _cache['all_periods'] = all_periods _cache['matched'] = matched return results def figures(): """Generate figures from cached main() results.""" if not _cache: main() paths = [] import matplotlib.pyplot as plt imfs = _cache.get('imfs') imf_info = _cache.get('imf_info', []) sig_periods = _cache.get('sig_periods', []) indpro_growth = _cache.get('indpro_growth') sector_neffs = _cache.get('sector_neffs', []) sector_sigmas = _cache.get('sector_sigmas', []) classical_cycles = _cache.get('classical_cycles', {}) all_periods = _cache.get('all_periods', []) if imfs is None: return [] fig, axes = setup_subplots(nrows=3, ncols=2, width=14, height=14) # (a) INDPRO IMFs stacked ax = axes[0, 0] sig_idx = _cache.get('sig_idx', []) for k, idx in enumerate(sig_idx): offset_val = k * 0.02 ax.plot(imfs[idx] + offset_val, linewidth=0.6, label=f'IMF {idx}') ax.set_title('(a) Significant INDPRO IMFs') ax.set_xlabel('Observation') ax.set_ylabel('Amplitude (stacked)') ax.legend(fontsize=7, ncol=2) # (b) r* histogram (adjacent timescale ratios) ax = axes[0, 1] ratios = _cache.get('ratios', []) if ratios: ax.hist(ratios, bins=max(5, len(ratios)), edgecolor='black', alpha=0.7) r_star = _cache.get('r_star', np.nan) if not np.isnan(r_star): ax.axvline(r_star, color='red', linestyle='--', label=f'median r*={r_star:.2f}') ax.axvline(2.0, color='green', linestyle=':', alpha=0.7, label='predicted r*=2') ax.set_title('(b) Adjacent Timescale Ratios') ax.set_xlabel('r* = T_{k+1} / T_k') ax.set_ylabel('Count') ax.legend(fontsize=8) # (c) N_eff vs sigma scatter ax = axes[1, 0] if sector_neffs and sector_sigmas: ax.scatter(sector_sigmas, sector_neffs, s=60, edgecolors='black', zorder=3) sector_names = _cache.get('sector_names', []) for k in range(min(len(sector_names), len(sector_sigmas))): ax.annotate(sector_names[k], (sector_sigmas[k], sector_neffs[k]), fontsize=6, ha='left', va='bottom') z = np.polyfit(sector_sigmas, sector_neffs, 1) x_line = np.linspace(min(sector_sigmas), max(sector_sigmas), 50) ax.plot(x_line, np.polyval(z, x_line), 'r--', alpha=0.6) ax.set_title('(c) Sector N_eff vs sigma') ax.set_xlabel('Substitution elasticity sigma') ax.set_ylabel('N_eff (significant IMFs)') # (d) Adjacent vs non-adjacent coupling bars ax = axes[1, 1] adj_corrs = _cache.get('adj_corrs', []) nonadj_corrs = _cache.get('nonadj_corrs', []) if adj_corrs or nonadj_corrs: means = [np.mean(adj_corrs) if adj_corrs else 0, np.mean(nonadj_corrs) if nonadj_corrs else 0] stds = [np.std(adj_corrs) if len(adj_corrs) > 1 else 0, np.std(nonadj_corrs) if len(nonadj_corrs) > 1 else 0] ax.bar(['Adjacent', 'Non-adjacent'], means, yerr=stds, capsize=5, color=['steelblue', 'salmon'], edgecolor='black') ax.set_title('(d) Energy Coupling: Adjacent vs Non-adjacent') ax.set_ylabel('Mean Spearman correlation') # (e) IMF periods vs classical cycle bands ax = axes[2, 0] if all_periods: periods_y = [p for p, _ in all_periods] energies_plot = [info['energy_fraction'] for info in imf_info if not np.isnan(info['char_period_years']) and info['energy_fraction'] > 0.01] ax.scatter(periods_y, energies_plot[:len(periods_y)], s=60, edgecolors='black', zorder=3) colors = ['#FFD700', '#87CEEB', '#90EE90', '#DDA0DD', '#FFA07A'] for k, (name, (lo, hi)) in enumerate(classical_cycles.items()): ax.axvspan(lo, hi, alpha=0.15, color=colors[k % len(colors)], label=name) ax.set_title('(e) IMF Periods vs Classical Cycles') ax.set_xlabel('Period (years)') ax.set_ylabel('Energy fraction') ax.set_xscale('log') ax.legend(fontsize=7, loc='upper left') # (f) Summary statistics table ax = axes[2, 1] ax.axis('off') n_eff = _cache.get('n_eff', '?') r_star_val = _cache.get('r_star', np.nan) matched = _cache.get('matched', set()) table_data = [ ['N_eff', str(n_eff)], ['r* (median)', f'{r_star_val:.2f}' if not np.isnan(r_star_val) else 'N/A'], ['Matched cycles', ', '.join(sorted(matched)) if matched else 'None'], ['N sectors', str(len(sector_neffs))], ] table = ax.table(cellText=table_data, colLabels=['Metric', 'Value'], loc='center', cellLoc='center') table.auto_set_font_size(False) table.set_fontsize(10) table.scale(1.0, 1.5) ax.set_title('(f) Summary Statistics') fig.suptitle('EMD Timescale Hierarchy', fontsize=14, y=1.01) fig.tight_layout() # Save as multi-panel figure covering all subtests paths.extend(save_figure(fig, 'emd_hierarchy', subfolder='hierarchy')) return paths if __name__ == '__main__': res = main() for r in res: print(f" {r['name']:45s} {r['verdict']:15s} stat={r['statistic']}")

Emd Hierarchy Adjacent Ratio

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources