#!/usr/bin/env python3 """ Geometric Mean Period Test Paper 3, Corollary 14.2 "Geometric mean period": T_nm^osc ~ 2*pi*sqrt(tau_n * tau_m) Tests whether predicted oscillation periods from autocorrelation times match actual spectral peak periods from cross-spectral analysis. Sub-tests: 3a: Spearman(T_predicted, T_actual) > 0.3 at p<0.05 3b: Median |T_predicted - T_actual| / T_actual < 0.50 Lean theorem: geometricMeanPeriod_pos in CESProofs/Dynamics/BusinessCycles.lean (proves positivity; oscillation_spectrum is schematic True := trivial) Review Notes (2026-03-12): Theory match: CORRECT — Paper 3, Corollary 14.2. Formula T_nm = 2π√(τ_n·τ_m) matches code exactly. No free parameters. I/O audit: CORRECT — τ from AR(1) autocorrelation of HP-filtered IP (input to formula), T_actual from cross-spectral peak (output to compare against). Data-signal alignment: ADEQUATE. AR(1) τ is a reasonable proxy for sector adjustment timescale. Cross-spectral peak is direct measure of co-oscillation. N=91 pairs from 14 Oberfield-Raval sectors. HP filter λ=129600 (monthly std). Frequency cutoff at 1/240 months (20 years) prevents long-period artifacts. Code quality: CLEAN (2026-03-12). Added _ALL_SUBTESTS, _cache, figures(). Confidence: HIGH — Spearman=0.428 (p<0.001) and median error 29% across 91 pairs with no fitted parameters. """ import warnings warnings.filterwarnings('ignore') import numpy as np import pandas as pd from scipy import stats from scipy.signal import welch, csd 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, hp_filter from common.plotting import setup_figure, setup_subplots, save_figure OR_SIGMA = { 'IPG336S': 0.50, 'IPG331S': 0.55, 'IPG327S': 0.60, 'IPG333S': 0.65, 'IPG321S': 0.70, 'IPG332S': 0.70, 'IPG335S': 0.70, 'IPG322S': 0.85, 'IPG337S': 0.85, 'IPG325S': 0.90, 'IPG326S': 0.95, 'IPG339S': 1.10, 'IPG311S': 1.15, 'IPG334S': 1.40, } OR_NAMES = { 'IPG336S': 'Transport Equip', 'IPG331S': 'Primary Metals', 'IPG327S': 'Nonmetallic Min', 'IPG333S': 'Machinery', 'IPG321S': 'Wood', 'IPG332S': 'Fabricated Metals', 'IPG335S': 'Electrical Equip', 'IPG322S': 'Paper', 'IPG337S': 'Furniture', 'IPG325S': 'Chemicals', 'IPG326S': 'Plastics/Rubber', 'IPG339S': 'Miscellaneous', 'IPG311S': 'Food', 'IPG334S': 'Computer/Elec', } _ALL_SUBTESTS = [ 'geometric_mean_period_spearman', 'geometric_mean_period_median_error', 'geometric_mean_period_overall', ] _cache = {} def main(): results = [] # Load data ip_series = {} for sid in OR_SIGMA: try: s = load_fred(sid) if len(s) >= 60: ip_series[sid] = s except FileNotFoundError: pass if len(ip_series) < 4: results.append(log_result("geometric_mean_period", np.nan, np.nan, "AMBIGUOUS", {"reason": "too few sectors"}, description="Tests whether predicted oscillation periods T_nm = 2*pi*sqrt(tau_n * tau_m) match actual cross-spectral peaks. Too few FRED IP sectors loaded.")) return results # Align to common range common_start = max(s.index.min() for s in ip_series.values()) common_end = min(s.index.max() for s in ip_series.values()) # HP-filter and estimate tau cycles = {} taus = {} for sid, s in ip_series.items(): sa = s[(s.index >= common_start) & (s.index <= common_end)].copy() if len(sa) < 60: continue log_ip = np.log(sa.values.astype(float)) cycle, trend = hp_filter(log_ip, lamb=129600) cycles[sid] = cycle phi_1 = np.corrcoef(cycle[:-1], cycle[1:])[0, 1] if 0.01 < phi_1 < 0.999: taus[sid] = -1.0 / np.log(phi_1) valid_sids = sorted(taus.keys()) if len(valid_sids) < 3: results.append(log_result("geometric_mean_period", np.nan, np.nan, "AMBIGUOUS", {"reason": "too few valid sectors"}, description="Tests geometric mean period prediction. Too few sectors with valid autocorrelation time estimates.")) return results # Pairwise cross-spectral analysis pair_results = [] for i, sid_n in enumerate(valid_sids): for sid_m in valid_sids[i + 1:]: T_pred = 2.0 * np.pi * np.sqrt(taus[sid_n] * taus[sid_m]) x = cycles[sid_n] y = cycles[sid_m] min_len = min(len(x), len(y)) x, y = x[:min_len], y[:min_len] if min_len < 48: continue nperseg = min_len nfft = 4 * nperseg freqs, Pxy = csd(x, y, fs=1.0, nperseg=nperseg, nfft=nfft, window='hann', noverlap=0) _, Pxx = welch(x, fs=1.0, nperseg=nperseg, nfft=nfft, window='hann', noverlap=0) _, Pyy = welch(y, fs=1.0, nperseg=nperseg, nfft=nfft, window='hann', noverlap=0) mask = freqs > 1.0 / 240.0 if mask.sum() == 0: continue freqs_f = freqs[mask] Pxy_f = Pxy[mask] Pxx_f = Pxx[mask] Pyy_f = Pyy[mask] denom = Pxx_f * Pyy_f denom[denom == 0] = 1e-30 coherence = np.abs(Pxy_f) ** 2 / denom max_coh = float(np.max(coherence)) abs_csd_vals = np.abs(Pxy_f) peak_idx = int(np.argmax(abs_csd_vals)) f_peak = float(freqs_f[peak_idx]) if f_peak <= 0: continue T_actual_peak = min(1.0 / f_peak, 240.0) pair_results.append({ 'T_pred': T_pred, 'T_peak': T_actual_peak, 'max_coherence': max_coh, 'sigma_n': OR_SIGMA.get(sid_n, 1.0), 'sigma_m': OR_SIGMA.get(sid_m, 1.0), }) df = pd.DataFrame(pair_results) if len(df) == 0: results.append(log_result("geometric_mean_period", np.nan, np.nan, "AMBIGUOUS", {"reason": "no valid pairs"}, description="Tests geometric mean period prediction. No valid sector pairs for cross-spectral analysis.")) return results df['rel_err'] = np.abs(df['T_pred'] - df['T_peak']) / df['T_peak'] # Coherence filter df_test = df[df['max_coherence'] > 0.3] if len(df[df['max_coherence'] > 0.3]) >= 5 else df # 3a: Spearman valid_p = df_test['T_peak'].notna() if valid_p.sum() >= 5: rho_s, p_s = stats.spearmanr(df_test.loc[valid_p, 'T_pred'], df_test.loc[valid_p, 'T_peak']) else: rho_s, p_s = np.nan, np.nan if rho_s > 0.3 and p_s < 0.05: verdict_3a = "CONSISTENT" elif rho_s > 0.2: verdict_3a = "AMBIGUOUS" else: verdict_3a = "INCONSISTENT" results.append(log_result("geometric_mean_period_spearman", rho_s, p_s, verdict_3a, {"n_pairs": len(df_test)}, description="Tests rank correlation between predicted oscillation periods T_pred = 2*pi*sqrt(tau_n * tau_m) and actual cross-spectral peak periods from Welch CSD analysis. Uses Spearman correlation across all sector pairs with coherence > 0.3. Data: HP-filtered FRED manufacturing IP subsectors. Prediction: Spearman rho > 0.3 at p < 0.05.")) # 3b: Median relative error valid_re = df_test['rel_err'].notna() med_re = df_test.loc[valid_re, 'rel_err'].median() if valid_re.sum() > 0 else 1.0 verdict_3b = "CONSISTENT" if med_re < 0.50 else "INCONSISTENT" results.append(log_result("geometric_mean_period_median_error", med_re, np.nan, verdict_3b, description="Tests calibration accuracy: median relative error |T_pred - T_actual| / T_actual across all sector pairs. Prediction: median relative error < 0.50, indicating the geometric mean formula provides order-of-magnitude accuracy.")) # Overall if verdict_3a == "CONSISTENT" and verdict_3b == "CONSISTENT": overall = "CONSISTENT" elif "AMBIGUOUS" in verdict_3a or rho_s > 0.2: overall = "AMBIGUOUS" else: overall = "INCONSISTENT" results.append(log_result("geometric_mean_period_overall", rho_s, p_s, overall, description="Combined verdict for geometric mean period test: requires both Spearman rank correlation and median relative error sub-tests to pass.")) _cache['results'] = results _cache['df_test'] = df_test _cache['df_all'] = df _cache['taus'] = taus _cache['valid_sids'] = valid_sids _cache['rho_s'] = rho_s _cache['p_s'] = p_s _cache['med_re'] = med_re return results def figures(): """Generate figures for geometric mean period subtests.""" if not _cache: main() paths = [] import matplotlib.pyplot as plt df_test = _cache.get('df_test') if df_test is None or len(df_test) == 0: return paths rho_s = _cache.get('rho_s', np.nan) p_s = _cache.get('p_s', np.nan) med_re = _cache.get('med_re', np.nan) # Figure 1: T_pred vs T_actual scatter (covers spearman + overall) fig, ax = setup_figure(width=8, height=6) ax.scatter(df_test['T_pred'], df_test['T_peak'], c='steelblue', s=40, alpha=0.6, edgecolor='k', linewidth=0.5, zorder=5) lims = [0, max(df_test['T_pred'].max(), df_test['T_peak'].max()) * 1.1] ax.plot(lims, lims, 'r--', alpha=0.5, label='Perfect prediction') ax.set_xlabel('Predicted period T = 2π√(τ_n·τ_m) [months]') ax.set_ylabel('Actual cross-spectral peak period [months]') ax.set_title(f'Geometric Mean Period: Predicted vs Actual\n' f'Spearman ρ={rho_s:.3f}, p={p_s:.2e}, N={len(df_test)} pairs') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) ax.set_xlim(lims) ax.set_ylim(lims) paths.extend(save_figure(fig, 'geometric_mean_period_spearman', subfolder='dynamics')) # Figure 2: Relative error histogram (covers median_error) fig, ax = setup_figure(width=8, height=5) rel_errs = df_test['rel_err'].dropna() ax.hist(rel_errs, bins=20, color='steelblue', edgecolor='k', alpha=0.7) ax.axvline(med_re, color='red', linestyle='--', linewidth=2, label=f'Median = {med_re:.3f}') ax.axvline(0.50, color='orange', linestyle=':', linewidth=1.5, label='Threshold = 0.50') ax.set_xlabel('Relative error |T_pred − T_actual| / T_actual') ax.set_ylabel('Count (sector pairs)') ax.set_title(f'Geometric Mean Period: Calibration Accuracy\n' f'Median relative error = {med_re:.1%}, N={len(rel_errs)} pairs') ax.legend(fontsize=9) ax.grid(True, alpha=0.3, axis='y') paths.extend(save_figure(fig, 'geometric_mean_period_median_error', subfolder='dynamics')) return paths if __name__ == '__main__': res = main() for r in res: print(f"{r['name']}: {r['verdict']} (stat={r['statistic']}, p={r['p_value']})")

p value	2.273e-5
statistic	0.4284

Geometric Mean Period Overall

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps