#!/usr/bin/env python3 """ Test: Natural Gradient (Slerp vs Linear) on GDP-by-Industry Shares =================================================================== Tests whether spherical linear interpolation (Slerp) on the simplex outperforms linear interpolation for GDP-by-industry share dynamics during recessions. The CES information geometry predicts Slerp is the natural geodesic on the share manifold (Amari 1998, Paper 1 InformationGeometry.lean). Lean theorem: `natural_gradient_agrees_at_equilibrium` — at symmetric equilibrium, Slerp and linear agree (Fisher-Rao distance = 0). The test checks whether they diverge during large share movements (recessions). Eleven sub-tests: T1: Paired Wilcoxon of per-quarter RSS (Slerp vs linear) during downturns T2: Paired Wilcoxon of per-quarter RSS (Slerp vs linear) during recoveries T3: Midpoint concentration (Slerp preserves simplex curvature) — binomial T4: Fisher-Rao distance orders recession severity — diagnostic T5: Pooled AIC comparison (delta-AIC and Akaike weight) T6: COVID exclusion sensitivity — T1 re-run without 2020 recession T7: O(theta^2) scaling — recession-only (N=5, underpowered) T8: Rolling Slerp vs Linear across ALL quarters (N=20 non-overlapping windows) T9: Rolling geometric scaling — ||slerp-linear||^2 vs theta (theory: slope~4) Note: T9 tests a mathematical identity (simplex geometry), not economic theory T10: Share volatility scaling — std(Δs_i) vs s_i^0.5 (Fisher-Rao metric test) THE KEY TEST: directly validates the Fisher-Rao metric on real data T11: Geodesic momentum — informative negative (no Slerp extrapolation advantage) Data: BEA quarterly GDP-by-industry value-added shares. Review Notes (2026-03-11): Theory match: CORRECT — Slerp is the Fisher-Rao geodesic under the Amari embedding xi_j = sqrt(s_j). The test checks whether real share dynamics follow this geodesic vs Euclidean (linear) interpolation. I/O audit: CORRECT — BEA value-added shares are the right data object. Industry shares lie on the simplex, which is the natural domain. Data-signal: WEAK — quarterly GDP shares change by <1% per quarter. Slerp and linear agree to O(theta^2) where theta ~ 0.03-0.13. BEA data starts 2005, so 2001 recession is missing. Effective N=1 recession for downturn RSS (GFC only; 2020 has 1 trivial quarter). Recovery phase has N=2 events (GFC + COVID, 8 quarters each). Code quality: BUGS FOUND — previous version had fabricated p-values: T1/T2 used p_val=min(ratio,1.0) (not a p-value), T4 used hardcoded p=0.05, T5 used hardcoded AIC-to-p mapping. All fixed with proper tests (paired Wilcoxon, diagnostic reframe, AIC weight). Confidence: LOW — fundamental power limitation: for small theta, Slerp ~ linear + O(theta^2), so detecting the difference requires either very large shocks or many more observations. Update (2026-03-11): Added T6 (COVID exclusion) and T7 (theta^2 scaling). T6 CONSISTENT (p=0.031): but only 1 recession (GFC) after COVID exclusion. T7 CONSISTENT (slope=0.83, n=5): underpowered; RSS-based confounds geometry with noise. Update (2026-03-11): Added T8 (rolling all quarters) and T9 (geometric scaling). T8 AMBIGUOUS: 10/20 windows favor Slerp (exactly 50%), no practical advantage. T9 CONSISTENT (slope=4.34, R^2=0.925, p<1e-10): geometric prediction confirmed! The theory IS geometrically correct (||slerp-linear||^2 ~ theta^4) but the effect is too small to improve predictions at typical quarterly share changes (theta~0.03). CAVEAT: T9 tests a mathematical identity — ||slerp-linear||^2 depends only on start/end shares, not actual intermediate data. It validates simplex geometry. Update (2026-03-11): Added T10 (share volatility scaling) and T11 (geodesic momentum). T10 CONSISTENT (slope=0.42, 95% CI contains 0.5, rejects 0.0 and 1.0): THE STRONGEST RESULT — directly validates Fisher-Rao metric on industry shares. Fisher-Rao predicts Var(Ds_i) proportional to s_i, so log(std) vs log(mean_s) has slope=0.5. T11 AMBIGUOUS (informative negative): no geodesic momentum in share dynamics. Slerp extrapolation performs identically to linear and worse than random walk. """ import numpy as np import pandas as pd from scipy import stats from scipy.special import rel_entr import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.config import BEA_API_KEY, RESULTS_DIR, FIGURES_DIR from common.plotting import setup_figure, setup_subplots, save_figure from common.stats import verdict, log_result # ── Geometric helpers ────────────────────────────────────────────── def bhattacharyya_distance(p, q): """Bhattacharyya distance between two probability vectors.""" p = np.asarray(p, dtype=float) q = np.asarray(q, dtype=float) p = p / p.sum() q = q / q.sum() bc = np.sum(np.sqrt(p * q)) bc = min(bc, 1.0) return -np.log(max(bc, 1e-30)) def fisher_rao_distance(p, q): """Fisher-Rao distance = 2 arccos(sum sqrt(p_i q_i)).""" p = np.asarray(p, dtype=float) q = np.asarray(q, dtype=float) p = p / p.sum() q = q / q.sum() bc = np.sum(np.sqrt(np.maximum(p, 1e-30) * np.maximum(q, 1e-30))) bc = min(bc, 1.0) return 2.0 * np.arccos(bc) def slerp(p, q, t): """Spherical linear interpolation on the positive orthant of the sphere.""" sp = np.sqrt(np.maximum(p, 1e-30)) sq = np.sqrt(np.maximum(q, 1e-30)) sp = sp / np.linalg.norm(sp) sq = sq / np.linalg.norm(sq) dot = np.clip(np.dot(sp, sq), -1.0, 1.0) omega = np.arccos(dot) if omega < 1e-10: result_sqrt = (1 - t) * sp + t * sq else: result_sqrt = (np.sin((1 - t) * omega) * sp + np.sin(t * omega) * sq) / np.sin(omega) result = result_sqrt ** 2 return result / result.sum() def linear_interp(p, q, t): """Standard linear interpolation on the simplex.""" result = (1 - t) * p + t * q return result / result.sum() def herfindahl(shares): """Herfindahl index of concentration.""" s = np.asarray(shares, dtype=float) s = s / s.sum() return np.sum(s ** 2) # ── BEA data loading ────────────────────────────────────────────── def fetch_bea_va_shares(): """ Load BEA quarterly GDP-by-industry value-added shares. Tries local cache first, then BEA API. """ from common.config import BEA_GDP_DIR, COMPILED_DIR # Try compiled cache (may be raw BEA format or pre-processed shares) cache = COMPILED_DIR / 'bea_quarterly_va_shares.csv' if cache.exists(): df = pd.read_csv(cache) # Check if already processed (date index, industry columns) if 'DataValue' not in df.columns and len(df) > 20: df = pd.read_csv(cache, index_col=0, parse_dates=True) return df # Raw BEA format: process into shares if 'DataValue' in df.columns and 'Industry' in df.columns: df['value'] = pd.to_numeric(df['DataValue'].astype(str).str.replace(',', ''), errors='coerce') df = df.dropna(subset=['value']) # Build quarter date — handle Roman numeral quarters roman_map = {'I': 1, 'II': 2, 'III': 3, 'IV': 4} df['qnum'] = df['Quarter'].map(roman_map).fillna(pd.to_numeric(df['Quarter'], errors='coerce')) df = df.dropna(subset=['qnum']) df['date'] = pd.to_datetime( df['Year'].astype(int).astype(str) + '-' + ((df['qnum'].astype(int) - 1) * 3 + 1).astype(str).str.zfill(2) + '-01') pivot = df.pivot_table(index='date', columns='Industry', values='value', aggfunc='first') drop_cols = [c for c in pivot.columns if c in ('ALL', 'TOTAL', '...') or 'total' in str(c).lower()] pivot = pivot.drop(columns=drop_cols, errors='ignore') shares = pivot.div(pivot.sum(axis=1), axis=0).dropna(how='all') if len(shares) > 20: return shares # Try BEA GDP directory for name in ['RealValueAdded_10_quarterly.csv', 'gdp_by_industry_quarterly.csv']: path = BEA_GDP_DIR / name if path.exists(): raw = pd.read_csv(path) return raw # Fetch from BEA API if key available if not BEA_API_KEY: raise FileNotFoundError( "BEA quarterly VA data not found and no API key available. " "Place bea_quarterly_va_shares.csv in data/compiled/" ) import requests url = "https://apps.bea.gov/api/data/" params = { 'UserID': BEA_API_KEY, 'method': 'GetData', 'DataSetName': 'GDPbyIndustry', 'TableID': '10', 'Frequency': 'Q', 'Year': ','.join(str(y) for y in range(2000, 2025)), 'Industry': 'ALL', 'ResultFormat': 'JSON' } resp = requests.get(url, params=params, timeout=60) resp.raise_for_status() data = resp.json() rows = data.get('BEAAPI', {}).get('Results', {}).get('Data', []) if not rows: raise ValueError("Empty BEA API response") records = [] for row in rows: try: val = float(row['DataValue'].replace(',', '')) records.append({ 'quarter': row['Quarter'], 'industry': row['Industry'], 'value': val }) except (ValueError, KeyError): continue df_raw = pd.DataFrame(records) pivot = df_raw.pivot_table(index='quarter', columns='industry', values='value', aggfunc='first') drop_cols = [c for c in pivot.columns if c in ('ALL', 'TOTAL', '...') or 'total' in str(c).lower()] pivot = pivot.drop(columns=drop_cols, errors='ignore') shares = pivot.div(pivot.sum(axis=1), axis=0) shares = shares.dropna(how='all') cache.parent.mkdir(parents=True, exist_ok=True) shares.to_csv(cache) return shares # ── NBER recession dates ────────────────────────────────────────── RECESSIONS = [ ('2001-03-01', '2001-11-01'), ('2007-12-01', '2009-06-01'), ('2020-02-01', '2020-04-01'), ] def get_recession_windows(shares_df, pre_quarters=4, post_quarters=8): """Extract pre/trough/post windows around each recession.""" windows = [] idx = pd.to_datetime(shares_df.index) for start_str, end_str in RECESSIONS: start = pd.Timestamp(start_str) end = pd.Timestamp(end_str) pre_mask = (idx >= start - pd.DateOffset(months=3 * pre_quarters)) & (idx < start) trough_mask = (idx >= start) & (idx <= end) post_mask = (idx > end) & (idx <= end + pd.DateOffset(months=3 * post_quarters)) pre = shares_df.loc[pre_mask] trough = shares_df.loc[trough_mask] post = shares_df.loc[post_mask] if len(pre) > 0 and len(trough) > 0: windows.append({ 'name': f"Recession {start.year}", 'pre': pre, 'trough': trough, 'post': post, 'start': start, 'end': end, }) return windows # ── Sub-tests ───────────────────────────────────────────────────── def _compute_paired_rss(windows, phase_key='trough'): """Compute per-quarter Slerp and linear RSS for a given phase. Returns arrays of per-quarter RSS: (slerp_rss_arr, linear_rss_arr, n_points). """ slerp_rss_list = [] linear_rss_list = [] for w in windows: pre_mean = w['pre'].mean(axis=0).values pre_mean = np.maximum(pre_mean, 1e-10) pre_mean = pre_mean / pre_mean.sum() phase = w[phase_key] if len(phase) == 0: continue if phase_key == 'trough': endpoint = phase.iloc[-1].values else: endpoint = phase.iloc[-1].values endpoint = np.maximum(endpoint, 1e-10) endpoint = endpoint / endpoint.sum() for i, (_, row) in enumerate(phase.iterrows()): actual = row.values actual = np.maximum(actual, 1e-10) actual = actual / actual.sum() t = (i + 1) / max(len(phase), 1) if phase_key == 'post': # Recovery: interpolate from trough to recovery end trough_end = w['trough'].iloc[-1].values trough_end = np.maximum(trough_end, 1e-10) trough_end = trough_end / trough_end.sum() pred_slerp = slerp(trough_end, endpoint, t) pred_linear = linear_interp(trough_end, endpoint, t) else: pred_slerp = slerp(pre_mean, endpoint, t) pred_linear = linear_interp(pre_mean, endpoint, t) slerp_rss_list.append(np.sum((actual - pred_slerp) ** 2)) linear_rss_list.append(np.sum((actual - pred_linear) ** 2)) return np.array(slerp_rss_list), np.array(linear_rss_list) def test_downturn_rss(windows): """T1: Paired Wilcoxon test — does Slerp have lower per-quarter RSS during downturns?""" slerp_arr, linear_arr = _compute_paired_rss(windows, 'trough') n = len(slerp_arr) t1_desc = ( "Tests whether Slerp (Fisher-Rao geodesic) achieves lower per-quarter RSS than " "linear interpolation during recession downturns. Paired Wilcoxon signed-rank test " "on per-quarter RSS differences (Slerp - Linear) across downturn quarters in " "NBER recessions (2007-09, 2020), using BEA quarterly GDP-by-industry " "value-added shares (starts 2005; 2001 recession not in data). " "H0: median(RSS_slerp - RSS_linear) >= 0, H1: < 0 (one-sided). " "Caveat: 2020 has only 1 trough quarter (trivially predicted), so effectively " "all evidence comes from the 2007-09 GFC (6 quarters). " "For small share changes (theta ~ 0.08), Slerp and linear agree to O(theta^2).") if n < 5: return log_result('natural_gradient_downturn_rss', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n} downturn quarters', 'n_points': n}, description=t1_desc) diff = slerp_arr - linear_arr ratio = slerp_arr.sum() / max(linear_arr.sum(), 1e-30) # Paired Wilcoxon signed-rank test (one-sided: Slerp < Linear) try: stat, p_twosided = stats.wilcoxon(slerp_arr, linear_arr, alternative='less') p_val = p_twosided # already one-sided with alternative='less' except Exception: stat, p_twosided = stats.wilcoxon(slerp_arr, linear_arr) p_val = p_twosided / 2 if diff.mean() < 0 else 1 - p_twosided / 2 sign_ok = ratio < 1.0 verd = verdict(p_val, sign_ok, threshold=0.10) return log_result( 'natural_gradient_downturn_rss', float(ratio), float(p_val), verd, details={ 'slerp_rss_total': float(slerp_arr.sum()), 'linear_rss_total': float(linear_arr.sum()), 'ratio': float(ratio), 'n_points': n, 'mean_diff': float(diff.mean()), 'median_diff': float(np.median(diff)), 'wilcoxon_stat': float(stat), }, description=t1_desc) def test_recovery_rss(windows): """T2: Paired Wilcoxon test — does Slerp have lower per-quarter RSS during recoveries?""" slerp_arr, linear_arr = _compute_paired_rss(windows, 'post') n = len(slerp_arr) t2_desc = ( "Tests whether Slerp achieves lower per-quarter RSS than linear interpolation " "during post-recession recoveries. Paired Wilcoxon signed-rank test on per-quarter " "RSS differences across 8 quarters after each NBER trough (2007-09, 2020), " "using BEA quarterly GDP-by-industry value-added shares. One-sided H1: Slerp < Linear. " "AMBIGUOUS: Slerp is slightly worse during recoveries (ratio > 1), suggesting " "recovery dynamics may be more linear than geodesic — possibly because recovery " "follows institutional/policy paths rather than manifold curvature.") if n < 5: return log_result('natural_gradient_recovery_rss', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n} recovery quarters', 'n_points': n}, description=t2_desc) diff = slerp_arr - linear_arr ratio = slerp_arr.sum() / max(linear_arr.sum(), 1e-30) try: stat, p_val = stats.wilcoxon(slerp_arr, linear_arr, alternative='less') except Exception: stat, p_twosided = stats.wilcoxon(slerp_arr, linear_arr) p_val = p_twosided / 2 if diff.mean() < 0 else 1 - p_twosided / 2 sign_ok = ratio < 1.0 verd = verdict(p_val, sign_ok, threshold=0.10) return log_result( 'natural_gradient_recovery_rss', float(ratio), float(p_val), verd, details={ 'slerp_rss_total': float(slerp_arr.sum()), 'linear_rss_total': float(linear_arr.sum()), 'ratio': float(ratio), 'n_points': n, 'mean_diff': float(diff.mean()), 'wilcoxon_stat': float(stat), }, description=t2_desc) def test_midpoint_concentration(windows): """T3: Binomial test — is Slerp midpoint more concentrated than linear midpoint?""" slerp_wins = 0 total = 0 hhi_pairs = [] for w in windows: pre_mean = w['pre'].mean(axis=0).values pre_mean = np.maximum(pre_mean, 1e-10) pre_mean = pre_mean / pre_mean.sum() if len(w['trough']) == 0: continue trough_mean = w['trough'].iloc[-1].values trough_mean = np.maximum(trough_mean, 1e-10) trough_mean = trough_mean / trough_mean.sum() mid_slerp = slerp(pre_mean, trough_mean, 0.5) mid_linear = linear_interp(pre_mean, trough_mean, 0.5) hhi_slerp = herfindahl(mid_slerp) hhi_linear = herfindahl(mid_linear) hhi_pairs.append((w['name'], hhi_slerp, hhi_linear)) if hhi_slerp > hhi_linear: slerp_wins += 1 total += 1 t3_desc = ( "Tests whether Slerp midpoints are more concentrated (higher HHI) than linear " "midpoints between pre-recession and trough industry shares. Binomial test across " "recession episodes (H0: P(Slerp more concentrated) = 0.5, one-sided). CES " "information geometry predicts Slerp preserves simplex curvature, yielding higher " "concentration at the midpoint. BEA quarterly GDP-by-industry shares.") if total == 0: return log_result('natural_gradient_midpoint_concentration', np.nan, np.nan, 'SKIPPED', details={'reason': 'No recession windows'}, description=t3_desc) frac = slerp_wins / total from scipy.stats import binomtest try: p_val = binomtest(slerp_wins, total, 0.5, alternative='greater').pvalue except Exception: p_val = 1.0 - stats.binom.cdf(slerp_wins - 1, total, 0.5) sign_correct = frac > 0.5 verd = verdict(p_val, sign_correct) return log_result( 'natural_gradient_midpoint_concentration', float(frac), float(p_val), verd, details={ 'slerp_wins': slerp_wins, 'total': total, 'hhi_pairs': [(n, round(s, 6), round(l, 6)) for n, s, l in hhi_pairs], }, description=t3_desc) def test_shock_magnitude(windows): """T4: Does Fisher-Rao distance order recession severity? Diagnostic.""" if len(windows) < 2: return log_result('natural_gradient_shock_ordering', np.nan, np.nan, 'SKIPPED', description="Diagnostic: Fisher-Rao distance vs recession severity. Insufficient data.") distances = [] for w in windows: pre_mean = w['pre'].mean(axis=0).values pre_mean = np.maximum(pre_mean, 1e-10) pre_mean = pre_mean / pre_mean.sum() if len(w['trough']) == 0: continue trough_mean = w['trough'].iloc[-1].values trough_mean = np.maximum(trough_mean, 1e-10) trough_mean = trough_mean / trough_mean.sum() fr = fisher_rao_distance(pre_mean, trough_mean) distances.append((w['name'], fr)) if len(distances) < 2: return log_result('natural_gradient_shock_ordering', np.nan, np.nan, 'SKIPPED', description="Diagnostic: Fisher-Rao distance vs recession severity. Insufficient data.") # Sort by distance distances_sorted = sorted(distances, key=lambda x: x[1], reverse=True) names = [d[0] for d in distances_sorted] vals = [d[1] for d in distances_sorted] # With N=2-3, can't do a proper test. Report as diagnostic. # Check if known severity ordering is matched: # COVID (large, short) > GFC (large, long) > 2001 (mild) ratio = vals[0] / max(vals[-1], 1e-10) # With only 2-3 recessions, assign AMBIGUOUS (no proper test possible) # Use Spearman if we have at least 3 points if len(distances) >= 3: # Known GDP decline ordering (approximate) known_severity = {'Recession 2020': 3, 'Recession 2007': 2, 'Recession 2001': 1} fr_vals = [d[1] for d in distances] sev_vals = [known_severity.get(d[0], 0) for d in distances] r_sp, p_sp = stats.spearmanr(fr_vals, sev_vals) sign_ok = r_sp > 0 verd = verdict(p_sp / 2 if r_sp > 0 else 1 - p_sp / 2, sign_ok, threshold=0.10) p_report = p_sp / 2 if r_sp > 0 else 1 - p_sp / 2 stat_report = r_sp else: verd = 'AMBIGUOUS' p_report = float('nan') stat_report = ratio return log_result( 'natural_gradient_shock_ordering', float(stat_report), float(p_report), verd, details={'distances': dict(distances), 'ratio_max_min': float(ratio), 'n_recessions': len(distances)}, description=( "Diagnostic: tests whether Fisher-Rao distance on the BEA industry share " "manifold correctly orders recession severity. Computes geodesic distance " "between pre-recession and trough shares for each NBER recession. Spearman " "correlation with known severity ranking if N >= 3; otherwise diagnostic only " "(N too small for inference). Fisher-Rao distance: 2*arccos(sum sqrt(p_i*q_i)).")) def test_pooled_aic(windows): """T5: Pooled AIC comparison with Akaike weight.""" from common.stats import compute_aic_bic slerp_rss = 0.0 linear_rss = 0.0 n_total = 0 J = 0 for w in windows: pre_mean = w['pre'].mean(axis=0).values pre_mean = np.maximum(pre_mean, 1e-10) pre_mean = pre_mean / pre_mean.sum() J = len(pre_mean) for phase_key in ['trough', 'post']: phase = w[phase_key] if len(phase) == 0: continue end_mean = phase.iloc[-1].values end_mean = np.maximum(end_mean, 1e-10) end_mean = end_mean / end_mean.sum() for i, (_, row) in enumerate(phase.iterrows()): actual = row.values actual = np.maximum(actual, 1e-10) actual = actual / actual.sum() t_param = (i + 1) / max(len(phase), 1) pred_s = slerp(pre_mean, end_mean, t_param) pred_l = linear_interp(pre_mean, end_mean, t_param) slerp_rss += np.sum((actual - pred_s) ** 2) linear_rss += np.sum((actual - pred_l) ** 2) n_total += 1 t5_desc = ( "Tests whether Slerp achieves lower AIC than linear interpolation pooled across " "all recession downturn and recovery phases. Computes AIC from RSS of Slerp vs " "linear predictions of BEA quarterly industry shares; reports delta-AIC and Akaike " "weight w_slerp = exp(-delta/2) / (1 + exp(-delta/2)). delta-AIC > 2 provides " "positive evidence; > 10 strong evidence (Burnham & Anderson 2002).") if n_total == 0: return log_result('natural_gradient_pooled_aic', np.nan, np.nan, 'SKIPPED', details={'reason': 'No data'}, description=t5_desc) T = n_total * J k_slerp = 2 # omega, t k_linear = 2 # slope, intercept aic_slerp, _ = compute_aic_bic(slerp_rss, k_slerp, T) aic_linear, _ = compute_aic_bic(linear_rss, k_linear, T) delta_aic = aic_linear - aic_slerp # positive = Slerp better # Akaike weight for Slerp if abs(delta_aic) < 500: # avoid overflow w_slerp = np.exp(delta_aic / 2) / (1 + np.exp(delta_aic / 2)) else: w_slerp = 1.0 if delta_aic > 0 else 0.0 sign_correct = delta_aic > 0 # Use delta-AIC > 2 as evidence threshold (Burnham & Anderson) if delta_aic > 10: verd = 'CONSISTENT' if sign_correct else 'INCONSISTENT' elif delta_aic > 2: verd = 'CONSISTENT' if sign_correct else 'AMBIGUOUS' else: verd = 'AMBIGUOUS' return log_result( 'natural_gradient_pooled_aic', float(delta_aic), float(w_slerp), verd, details={ 'aic_slerp': float(aic_slerp), 'aic_linear': float(aic_linear), 'delta_aic': float(delta_aic), 'akaike_weight_slerp': float(w_slerp), 'slerp_rss': float(slerp_rss), 'linear_rss': float(linear_rss), 'n_obs': n_total, 'n_dimensions': J, }, description=t5_desc) def test_downturn_no_covid(windows): """T6: COVID exclusion — re-run T1 Wilcoxon excluding 2020 recession.""" windows_no_covid = [w for w in windows if '2020' not in w['name']] t6_desc = ( "COVID exclusion sensitivity: re-runs the downturn RSS Wilcoxon test (T1) " "excluding the 2020 pandemic recession. Tests whether the Slerp advantage " "during downturns survives when COVID is removed. " "Paired Wilcoxon on BEA quarterly GDP-by-industry shares for 2007-09 " "recession only (2001 outside BEA data range, starts 2005). " "One-sided H1: Slerp < Linear. " "Caveat: reduces to single-event test (n=6 from GFC only).") slerp_arr, linear_arr = _compute_paired_rss(windows_no_covid, 'trough') n = len(slerp_arr) if n < 5: return log_result('natural_gradient_no_covid', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n} downturn quarters without COVID', 'n_points': n}, description=t6_desc) diff = slerp_arr - linear_arr ratio = slerp_arr.sum() / max(linear_arr.sum(), 1e-30) try: stat, p_val = stats.wilcoxon(slerp_arr, linear_arr, alternative='less') except Exception: stat, p_twosided = stats.wilcoxon(slerp_arr, linear_arr) p_val = p_twosided / 2 if diff.mean() < 0 else 1 - p_twosided / 2 sign_ok = ratio < 1.0 verd = verdict(p_val, sign_ok, threshold=0.10) return log_result( 'natural_gradient_no_covid', float(ratio), float(p_val), verd, details={ 'slerp_rss_total': float(slerp_arr.sum()), 'linear_rss_total': float(linear_arr.sum()), 'ratio': float(ratio), 'n_points': n, 'mean_diff': float(diff.mean()), 'median_diff': float(np.median(diff)), 'wilcoxon_stat': float(stat), 'recessions_used': [w['name'] for w in windows_no_covid], }, description=t6_desc) def _compute_per_quarter_theta_and_diff(windows, phase_key='trough'): """Compute per-quarter Fisher-Rao theta and |RSS_slerp - RSS_linear|. Returns arrays: (thetas, abs_diffs, slerp_rss, linear_rss). """ thetas = [] abs_diffs = [] slerp_list = [] linear_list = [] for w in windows: pre_mean = w['pre'].mean(axis=0).values pre_mean = np.maximum(pre_mean, 1e-10) pre_mean = pre_mean / pre_mean.sum() phase = w[phase_key] if len(phase) == 0: continue endpoint = phase.iloc[-1].values endpoint = np.maximum(endpoint, 1e-10) endpoint = endpoint / endpoint.sum() for i, (_, row) in enumerate(phase.iterrows()): actual = row.values actual = np.maximum(actual, 1e-10) actual = actual / actual.sum() t_param = (i + 1) / max(len(phase), 1) # Theta: Fisher-Rao distance from pre_mean to actual theta = fisher_rao_distance(pre_mean, actual) pred_slerp = slerp(pre_mean, endpoint, t_param) pred_linear = linear_interp(pre_mean, endpoint, t_param) rss_s = np.sum((actual - pred_slerp) ** 2) rss_l = np.sum((actual - pred_linear) ** 2) thetas.append(theta) abs_diffs.append(abs(rss_s - rss_l)) slerp_list.append(rss_s) linear_list.append(rss_l) return np.array(thetas), np.array(abs_diffs), np.array(slerp_list), np.array(linear_list) def test_theta_squared_scaling(windows): """T7: Tests whether |RSS_slerp - RSS_linear| scales as O(theta^2).""" thetas, abs_diffs, _, _ = _compute_per_quarter_theta_and_diff(windows, 'trough') t7_desc = ( "Tests the O(theta^2) scaling prediction: the Slerp-linear RSS difference " "should grow quadratically with Fisher-Rao distance theta on the share manifold. " "Computes per-quarter theta (geodesic distance from pre-recession shares to actual) " "and |RSS_slerp - RSS_linear|. OLS of log|diff| on log(theta) estimates the scaling " "exponent; theory predicts slope ~ 2. Spearman correlation tests monotonicity. " "BEA quarterly GDP-by-industry shares across NBER recessions.") # Filter out zero diffs (identical Slerp/linear) mask = (thetas > 1e-10) & (abs_diffs > 1e-20) thetas_f = thetas[mask] diffs_f = abs_diffs[mask] n = len(thetas_f) if n < 5: return log_result('natural_gradient_theta_scaling', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n} valid points', 'n_points': n}, description=t7_desc) # Log-log regression: log|diff| = a + b*log(theta) log_theta = np.log(thetas_f) log_diff = np.log(diffs_f) slope, intercept, r_val, p_val_ols, se = stats.linregress(log_theta, log_diff) # Spearman: does larger theta -> larger diff? r_sp, p_sp = stats.spearmanr(thetas_f, diffs_f) # Theory predicts slope ~ 2 (quadratic scaling) # Sign check: slope should be positive (larger theta -> larger diff) # AND slope should be near 2 sign_ok = slope > 0 and r_sp > 0 # Use Spearman p for verdict (more robust with small N) p_one = p_sp / 2 if r_sp > 0 else 1 - p_sp / 2 verd = verdict(p_one, sign_ok, threshold=0.10) _cache['theta_scaling'] = { 'thetas': thetas_f, 'diffs': diffs_f, 'log_theta': log_theta, 'log_diff': log_diff, 'slope': slope, 'intercept': intercept, } return log_result( 'natural_gradient_theta_scaling', float(slope), float(p_one), verd, details={ 'ols_slope': float(slope), 'ols_intercept': float(intercept), 'ols_r_squared': float(r_val ** 2), 'ols_se': float(se), 'ols_p': float(p_val_ols), 'spearman_r': float(r_sp), 'spearman_p': float(p_sp), 'n_points': n, 'theory_prediction': 'slope ~ 2 (quadratic)', 'theta_range': [float(thetas_f.min()), float(thetas_f.max())], }, description=t7_desc) # ── Rolling window tests (use ALL quarters, not just recessions) ── def _compute_rolling_window_stats(shares, h=4): """Compute per-window Slerp vs Linear RSS for non-overlapping h-quarter windows. Returns list of dicts with theta, slerp_rss, linear_rss, rss_diff per window. """ n_dates = len(shares) window_data = [] for w_idx in range(n_dates // h): t_start = w_idx * h t_end = t_start + h if t_end >= n_dates: break start = shares.iloc[t_start].values.copy() start = np.maximum(start, 1e-10) start = start / start.sum() end = shares.iloc[t_end].values.copy() end = np.maximum(end, 1e-10) end = end / end.sum() theta = fisher_rao_distance(start, end) slerp_rss = 0.0 linear_rss = 0.0 pred_dist = 0.0 for k in range(1, h): actual = shares.iloc[t_start + k].values.copy() actual = np.maximum(actual, 1e-10) actual = actual / actual.sum() t_param = k / h pred_s = slerp(start, end, t_param) pred_l = linear_interp(start, end, t_param) slerp_rss += np.sum((actual - pred_s) ** 2) linear_rss += np.sum((actual - pred_l) ** 2) pred_dist += np.sum((pred_s - pred_l) ** 2) window_data.append({ 'window': w_idx, 'date_start': shares.index[t_start], 'date_end': shares.index[t_end], 'theta': theta, 'slerp_rss': slerp_rss, 'linear_rss': linear_rss, 'rss_diff': slerp_rss - linear_rss, # negative = Slerp better 'pred_dist': pred_dist, # ||slerp - linear||^2, for geometric scaling test }) return window_data def test_rolling_all_quarters(shares): """T8: Rolling Slerp vs Linear across ALL quarters, not just recessions. Uses non-overlapping 4-quarter windows across the full BEA sample. Each window interpolates from s_t to s_{t+4} and predicts s_{t+1}, s_{t+2}, s_{t+3}. Paired Wilcoxon on per-window RSS differences. """ t8_desc = ( "Tests Slerp vs linear interpolation using ALL quarterly data, not just recession " "episodes. Non-overlapping 4-quarter (1-year) windows across the full BEA sample " "(2005-2025, N ~ 20 independent windows). For each window, predicts 3 intermediate " "quarters via both Slerp and linear interpolation between start and end shares. " "Paired Wilcoxon signed-rank test on per-window total RSS differences. " "H1: Slerp < Linear (one-sided). This test dramatically improves statistical power " "over the recession-only tests (N=1-2 events -> N ~ 20 independent windows).") window_data = _compute_rolling_window_stats(shares, h=4) _cache['rolling_window_data'] = window_data n = len(window_data) if n < 8: return log_result('natural_gradient_rolling_all', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n} windows', 'n_windows': n}, description=t8_desc) diffs = np.array([w['rss_diff'] for w in window_data]) slerp_total = sum(w['slerp_rss'] for w in window_data) linear_total = sum(w['linear_rss'] for w in window_data) ratio = slerp_total / max(linear_total, 1e-30) n_slerp_better = int(np.sum(diffs < 0)) try: stat, p_val = stats.wilcoxon(diffs, alternative='less') except Exception: stat, p_twosided = stats.wilcoxon(diffs) p_val = p_twosided / 2 if np.mean(diffs) < 0 else 1 - p_twosided / 2 sign_ok = ratio < 1.0 verd = verdict(p_val, sign_ok, threshold=0.10) return log_result( 'natural_gradient_rolling_all', float(ratio), float(p_val), verd, details={ 'n_windows': n, 'n_slerp_better': n_slerp_better, 'ratio': float(ratio), 'slerp_total': float(slerp_total), 'linear_total': float(linear_total), 'mean_theta': float(np.mean([w['theta'] for w in window_data])), 'median_theta': float(np.median([w['theta'] for w in window_data])), 'theta_range': f"{min(w['theta'] for w in window_data):.4f}-{max(w['theta'] for w in window_data):.4f}", 'wilcoxon_stat': float(stat), }, description=t8_desc) def test_rolling_theta_scaling(shares): """T9: Theta-squared scaling using rolling window prediction distances. Uses ||pred_slerp - pred_linear|| (the geometric distance between the two predictions) rather than |RSS_slerp - RSS_linear| (which confounds geometry with theta-dependent noise). Theory predicts ||slerp - linear|| ~ O(theta^2), so log(pred_dist) vs log(theta) should have slope ~ 4 (since pred_dist is squared distance: ||delta||^2 ~ theta^4). """ t9_desc = ( "Tests O(theta^2) scaling prediction using the GEOMETRIC distance between Slerp " "and linear predictions: ||slerp(t) - linear(t)||^2 vs theta (Fisher-Rao distance " "between window endpoints). Theory: slerp and linear diverge as O(theta^2), so " "squared divergence scales as theta^4. OLS of log(pred_dist) on log(theta) should " "yield slope ~ 4. Uses non-overlapping 4-quarter rolling windows across full BEA " "sample (N ~ 20), not just recession quarters. This avoids the confound in T7 " "where |RSS_diff| mixes geometric scaling with theta-dependent prediction noise.") window_data = _cache.get('rolling_window_data') if not window_data: window_data = _compute_rolling_window_stats(shares, h=4) _cache['rolling_window_data'] = window_data thetas = np.array([w['theta'] for w in window_data]) pred_dists = np.array([w['pred_dist'] for w in window_data]) mask = (thetas > 1e-6) & (pred_dists > 1e-30) thetas_f = thetas[mask] dists_f = pred_dists[mask] n = len(thetas_f) if n < 8: return log_result('natural_gradient_rolling_theta_scaling', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n} valid windows'}, description=t9_desc) log_theta = np.log(thetas_f) log_dist = np.log(dists_f) # Theory: ||slerp - linear||^2 ~ theta^4 (since ||delta|| ~ theta^2) slope, intercept, r_val, p_val_ols, se = stats.linregress(log_theta, log_dist) r_sp, p_sp = stats.spearmanr(thetas_f, dists_f) sign_ok = slope > 0 and r_sp > 0 p_one = p_sp / 2 if r_sp > 0 else 1 - p_sp / 2 verd = verdict(p_one, sign_ok, threshold=0.10) _cache['rolling_theta_scaling'] = { 'thetas': thetas_f, 'diffs': dists_f, 'log_theta': log_theta, 'log_diff': log_dist, 'slope': slope, 'intercept': intercept, } return log_result( 'natural_gradient_rolling_theta_scaling', float(slope), float(p_one), verd, details={ 'ols_slope': float(slope), 'ols_intercept': float(intercept), 'ols_r_squared': float(r_val ** 2), 'ols_se': float(se), 'ols_p': float(p_val_ols), 'spearman_r': float(r_sp), 'spearman_p': float(p_sp), 'n_windows': n, 'theory_prediction': 'slope ~ 4 (||delta||^2 ~ theta^4)', 'theta_range': [float(thetas_f.min()), float(thetas_f.max())], }, description=t9_desc) # ── Direct metric tests ────────────────────────────────────────── def test_share_volatility_scaling(shares): """T10: Share volatility scaling — THE KEY TEST for Fisher-Rao metric. The Fisher-Rao metric on the simplex predicts that the local variance of share changes scales linearly with the share level: Var(Δs_i) ∝ s_i. Equivalently, std(Δs_i) ∝ s_i^0.5, so log(std) vs log(mean_s) has slope 0.5. Alternative predictions: - Euclidean (flat metric): slope = 0 (volatility independent of level) - Log-normal: slope = 1.0 (volatility proportional to level) - Fisher-Rao: slope = 0.5 (volatility proportional to sqrt of level) This tests the METRIC ITSELF, not geodesic paths, making it the most direct test of CES information geometry on real economic data. """ t10_desc = ( "Tests the Fisher-Rao metric prediction that share change volatility scales as " "std(Delta_s_i) proportional to s_i^0.5 across industries. For each BEA industry, computes " "time-series std of quarterly share changes vs mean share level. OLS of " "log(std) on log(mean_share) estimates the scaling exponent. Fisher-Rao predicts " "slope = 0.5; Euclidean predicts 0; log-normal predicts 1.0. " "Data: BEA quarterly GDP-by-industry value-added shares (2005-2025), N ~ 75 industries. " "Prediction: slope in [0.3, 0.7] containing 0.5, rejecting both 0 and 1.") # Compute per-industry volatility statistics deltas = shares.diff().iloc[1:] # quarterly share changes mean_shares = shares.mean(axis=0) std_deltas = deltas.std(axis=0) # Filter: need positive mean share and positive std mask = (mean_shares > 1e-6) & (std_deltas > 1e-12) mean_s = mean_shares[mask].values std_d = std_deltas[mask].values n = len(mean_s) if n < 20: return log_result('natural_gradient_volatility_scaling', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n} valid industries', 'n_industries': n}, description=t10_desc) log_mean = np.log(mean_s) log_std = np.log(std_d) slope, intercept, r_val, p_val_ols, se = stats.linregress(log_mean, log_std) # Test against theory predictions # H0: slope = 0.5 (Fisher-Rao) t_stat_fr = (slope - 0.5) / se p_fr = 2 * stats.t.sf(abs(t_stat_fr), n - 2) # H0: slope = 0 (Euclidean) t_stat_euc = slope / se p_euc = 2 * stats.t.sf(abs(t_stat_euc), n - 2) # H0: slope = 1.0 (log-normal) t_stat_ln = (slope - 1.0) / se p_ln = 2 * stats.t.sf(abs(t_stat_ln), n - 2) # 95% CI ci_low = slope - 1.96 * se ci_high = slope + 1.96 * se # Verdict: CONSISTENT if slope is closer to 0.5 than to 0 or 1, and CI contains 0.5 sign_ok = (0.2 < slope < 0.8) # broadly in Fisher-Rao range contains_half = (ci_low <= 0.5 <= ci_high) rejects_zero = (p_euc < 0.05) rejects_one = (p_ln < 0.05) if contains_half and rejects_zero and rejects_one: verd = 'CONSISTENT' elif sign_ok and (rejects_zero or rejects_one): verd = 'CONSISTENT' elif sign_ok: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['volatility_scaling'] = { 'log_mean': log_mean, 'log_std': log_std, 'slope': slope, 'intercept': intercept, 'se': se, 'r_squared': r_val ** 2, 'n': n, 'ci_low': ci_low, 'ci_high': ci_high, } return log_result( 'natural_gradient_volatility_scaling', float(slope), float(p_fr), verd, details={ 'ols_slope': float(slope), 'ols_se': float(se), 'ols_r_squared': float(r_val ** 2), 'ols_p': float(p_val_ols), 'ci_95': [float(ci_low), float(ci_high)], 'n_industries': n, 'p_reject_euclidean': float(p_euc), 'p_reject_lognormal': float(p_ln), 'p_reject_fisher_rao': float(p_fr), 'theory_predictions': { 'fisher_rao': 0.5, 'euclidean': 0.0, 'lognormal': 1.0, }, }, description=t10_desc) def test_geodesic_momentum(shares): """T11: Geodesic momentum test — informative negative. Tests whether Slerp extrapolation from past share changes predicts future shares better than linear extrapolation or no-change (random walk). Theory would predict geodesic momentum if adjustment dynamics follow natural gradient flow. Result: NO geodesic momentum detected — share dynamics follow random walk, not geodesics. """ t11_desc = ( "Tests geodesic momentum: does Slerp extrapolation of past share changes predict " "future shares? For each quarter, extrapolates from s_{t-1} through s_t to predict " "s_{t+1} via both Slerp and linear extrapolation. Compares RSS against no-change " "benchmark (predict s_{t+1} = s_t). Wilcoxon test of Slerp vs Linear RSS. " "Data: BEA quarterly GDP-by-industry shares. " "RESULT: Informative negative — no geodesic momentum. Both extrapolation methods " "are worse than random walk (ratio ~ 1.9), and Slerp does not beat linear. " "Share dynamics lack geodesic inertia at quarterly frequency.") n_dates = len(shares) if n_dates < 10: return log_result('natural_gradient_geodesic_momentum', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {n_dates} dates'}, description=t11_desc) slerp_rss_list = [] linear_rss_list = [] nochange_rss_list = [] for t in range(1, n_dates - 1): prev = shares.iloc[t - 1].values.copy() prev = np.maximum(prev, 1e-10) prev = prev / prev.sum() curr = shares.iloc[t].values.copy() curr = np.maximum(curr, 1e-10) curr = curr / curr.sum() actual_next = shares.iloc[t + 1].values.copy() actual_next = np.maximum(actual_next, 1e-10) actual_next = actual_next / actual_next.sum() # Slerp extrapolation: continue geodesic from prev through curr to t=2 pred_slerp = slerp(prev, curr, 2.0) # Linear extrapolation: curr + (curr - prev), renormalized pred_linear = 2 * curr - prev pred_linear = np.maximum(pred_linear, 1e-10) pred_linear = pred_linear / pred_linear.sum() # No-change benchmark pred_nc = curr.copy() slerp_rss_list.append(np.sum((actual_next - pred_slerp) ** 2)) linear_rss_list.append(np.sum((actual_next - pred_linear) ** 2)) nochange_rss_list.append(np.sum((actual_next - pred_nc) ** 2)) slerp_arr = np.array(slerp_rss_list) linear_arr = np.array(linear_rss_list) nochange_arr = np.array(nochange_rss_list) n = len(slerp_arr) ratio_sl = slerp_arr.sum() / max(linear_arr.sum(), 1e-30) ratio_snc = slerp_arr.sum() / max(nochange_arr.sum(), 1e-30) n_slerp_beats_linear = int(np.sum(slerp_arr < linear_arr)) # Wilcoxon: is Slerp better than Linear? try: stat, p_val = stats.wilcoxon(slerp_arr, linear_arr, alternative='less') except Exception: stat, p_twosided = stats.wilcoxon(slerp_arr, linear_arr) p_val = p_twosided / 2 if slerp_arr.mean() < linear_arr.mean() else 1 - p_twosided / 2 sign_ok = ratio_sl < 1.0 verd = verdict(p_val, sign_ok, threshold=0.10) _cache['geodesic_momentum'] = { 'slerp_arr': slerp_arr, 'linear_arr': linear_arr, 'nochange_arr': nochange_arr, 'ratio_sl': ratio_sl, 'ratio_snc': ratio_snc, } return log_result( 'natural_gradient_geodesic_momentum', float(ratio_sl), float(p_val), verd, details={ 'n_quarters': n, 'slerp_linear_ratio': float(ratio_sl), 'slerp_nochange_ratio': float(ratio_snc), 'n_slerp_beats_linear': n_slerp_beats_linear, 'wilcoxon_stat': float(stat), 'interpretation': 'No geodesic momentum: both extrapolations worse than random walk', }, description=t11_desc) # ── Main ────────────────────────────────────────────────────────── _cache = {} def main(): """Run all natural gradient tests. Returns list of result dicts.""" results = [] try: shares = fetch_bea_va_shares() except (FileNotFoundError, Exception) as e: print(f" Data loading failed: {e}") return [log_result('natural_gradient', 0.0, 1.0, 'SKIPPED', description="Skipped: BEA quarterly GDP-by-industry value-added share data could not be loaded.", details={'error': str(e)})] # Ensure shares is a proper DataFrame with datetime index if isinstance(shares, pd.DataFrame): if not isinstance(shares.index, pd.DatetimeIndex): try: shares.index = pd.to_datetime(shares.index) except Exception: pass # Ensure numeric and normalized shares = shares.apply(pd.to_numeric, errors='coerce') shares = shares.dropna(how='all', axis=1).dropna(how='all', axis=0) shares = shares.div(shares.sum(axis=1), axis=0) windows = get_recession_windows(shares) if not windows: return [log_result('natural_gradient', 0.0, 1.0, 'SKIPPED', description="Skipped: no NBER recession windows found within the BEA data.", details={'error': 'No recession windows found in data'})] _cache['shares'] = shares _cache['windows'] = windows results.append(test_downturn_rss(windows)) results.append(test_recovery_rss(windows)) results.append(test_midpoint_concentration(windows)) results.append(test_shock_magnitude(windows)) results.append(test_pooled_aic(windows)) results.append(test_downturn_no_covid(windows)) results.append(test_theta_squared_scaling(windows)) # Rolling window tests — use ALL quarters, not just recessions results.append(test_rolling_all_quarters(shares)) results.append(test_rolling_theta_scaling(shares)) # Direct metric tests results.append(test_share_volatility_scaling(shares)) results.append(test_geodesic_momentum(shares)) _cache['results'] = results return results def figures(): """Generate figures for natural gradient subtests.""" if not _cache: main() paths = [] windows = _cache.get('windows') if not windows: return paths # Figure 1: Per-quarter RSS comparison (Slerp vs Linear) — downturn + recovery slerp_down, linear_down = _compute_paired_rss(windows, 'trough') slerp_rec, linear_rec = _compute_paired_rss(windows, 'post') if len(slerp_down) > 0 or len(slerp_rec) > 0: fig, (ax1, ax2) = setup_subplots(nrows=1, ncols=2, width=14, height=6) if len(slerp_down) > 0: x = np.arange(len(slerp_down)) ax1.bar(x - 0.15, slerp_down * 1e6, 0.3, color='steelblue', label='Slerp', alpha=0.8) ax1.bar(x + 0.15, linear_down * 1e6, 0.3, color='darkorange', label='Linear', alpha=0.8) ax1.set_xlabel('Quarter') ax1.set_ylabel('RSS ($\\times 10^{-6}$)') ax1.set_title('Downturn: Per-Quarter RSS') ax1.legend() if len(slerp_rec) > 0: x = np.arange(len(slerp_rec)) ax2.bar(x - 0.15, slerp_rec * 1e6, 0.3, color='steelblue', label='Slerp', alpha=0.8) ax2.bar(x + 0.15, linear_rec * 1e6, 0.3, color='darkorange', label='Linear', alpha=0.8) ax2.set_xlabel('Quarter') ax2.set_ylabel('RSS ($\\times 10^{-6}$)') ax2.set_title('Recovery: Per-Quarter RSS') ax2.legend() fig.suptitle('Natural Gradient: Slerp vs Linear RSS by Phase', fontsize=13) paths += save_figure(fig, 'natural_gradient_rss_comparison', subfolder='curvature_roles') # Figure 2: Fisher-Rao distances by recession fr_distances = [] for w in windows: pre_mean = w['pre'].mean(axis=0).values pre_mean = np.maximum(pre_mean, 1e-10) pre_mean = pre_mean / pre_mean.sum() if len(w['trough']) == 0: continue trough_mean = w['trough'].iloc[-1].values trough_mean = np.maximum(trough_mean, 1e-10) trough_mean = trough_mean / trough_mean.sum() fr = fisher_rao_distance(pre_mean, trough_mean) fr_distances.append((w['name'], fr)) if len(fr_distances) >= 2: fig, ax = setup_figure(width=8, height=5) names = [d[0] for d in fr_distances] vals = [d[1] for d in fr_distances] colors = ['#e74c3c' if 'COVID' in n or '2020' in n else '#3498db' for n in names] ax.bar(range(len(names)), vals, color=colors) ax.set_xticks(range(len(names))) ax.set_xticklabels(names) ax.set_ylabel('Fisher-Rao Distance') ax.set_title('Recession Severity on the Share Manifold') for i, v in enumerate(vals): ax.annotate(f'{v:.4f}', (i, v), textcoords='offset points', xytext=(0, 5), ha='center', fontsize=10) paths += save_figure(fig, 'natural_gradient_fisher_rao', subfolder='curvature_roles') # Figure 3: Theta-squared scaling (log-log plot) ts = _cache.get('theta_scaling') if ts is not None and len(ts.get('thetas', [])) >= 3: fig, ax = setup_figure(width=8, height=6) ax.scatter(ts['log_theta'], ts['log_diff'], alpha=0.7, color='steelblue', edgecolors='navy', s=50, zorder=3) # OLS fit line x_fit = np.linspace(ts['log_theta'].min(), ts['log_theta'].max(), 50) y_fit = ts['slope'] * x_fit + ts['intercept'] ax.plot(x_fit, y_fit, 'r-', linewidth=2, label=f'OLS: slope={ts["slope"]:.2f} (theory: 2.0)') # Reference slope=2 line through data centroid cx, cy = ts['log_theta'].mean(), ts['log_diff'].mean() y_ref = 2.0 * (x_fit - cx) + cy ax.plot(x_fit, y_ref, 'k--', linewidth=1, alpha=0.5, label='Reference slope=2') ax.set_xlabel(r'$\log(\theta)$ (Fisher-Rao distance)') ax.set_ylabel(r'$\log(|RSS_{slerp} - RSS_{linear}|)$') ax.set_title(r'Slerp-Linear Difference vs $\theta$: Testing $O(\theta^2)$ Scaling') ax.legend() paths += save_figure(fig, 'natural_gradient_theta_scaling', subfolder='curvature_roles') # Figure 4: Rolling window RSS ratios over time rolling_data = _cache.get('rolling_window_data') if rolling_data and len(rolling_data) >= 3: fig, (ax1, ax2) = setup_subplots(nrows=1, ncols=2, width=14, height=5) # Left: per-window RSS ratio over time dates = [w['date_start'] for w in rolling_data] ratios = [w['slerp_rss'] / max(w['linear_rss'], 1e-30) for w in rolling_data] colors = ['steelblue' if r < 1 else 'darkorange' for r in ratios] ax1.bar(range(len(ratios)), ratios, color=colors, alpha=0.8, edgecolor='black', linewidth=0.3) ax1.axhline(1.0, color='black', linestyle='--', alpha=0.5) ax1.set_xlabel('Window (4-quarter non-overlapping)') ax1.set_ylabel('RSS ratio (Slerp / Linear)') n_slerp = sum(1 for r in ratios if r < 1) ax1.set_title(f'Rolling Slerp vs Linear: {n_slerp}/{len(ratios)} windows favor Slerp') ax1.text(0.02, 0.98, 'Blue = Slerp better\nOrange = Linear better', transform=ax1.transAxes, va='top', fontsize=8, style='italic') # Right: theta vs RSS difference (rolling version of theta scaling) thetas_r = np.array([w['theta'] for w in rolling_data]) diffs_r = np.array([w['rss_diff'] for w in rolling_data]) colors_r = ['steelblue' if d < 0 else 'darkorange' for d in diffs_r] ax2.scatter(thetas_r, diffs_r * 1e9, c=colors_r, s=50, alpha=0.8, edgecolors='black', linewidth=0.3) ax2.axhline(0, color='black', linestyle='--', alpha=0.5) ax2.set_xlabel(r'$\theta$ (Fisher-Rao distance)') ax2.set_ylabel(r'RSS diff (Slerp $-$ Linear) $\times 10^{-9}$') ax2.set_title(r'Per-Window Slerp Advantage vs Share Distance $\theta$') fig.tight_layout() paths += save_figure(fig, 'natural_gradient_rolling_all', subfolder='curvature_roles') # Figure 5: Rolling theta scaling (log-log) rts = _cache.get('rolling_theta_scaling') if rts is not None and len(rts.get('thetas', [])) >= 3: fig, ax = setup_figure(width=8, height=6) ax.scatter(rts['log_theta'], rts['log_diff'], alpha=0.7, color='steelblue', edgecolors='navy', s=60, zorder=3) x_fit = np.linspace(rts['log_theta'].min(), rts['log_theta'].max(), 50) y_fit = rts['slope'] * x_fit + rts['intercept'] ax.plot(x_fit, y_fit, 'r-', linewidth=2, label=f'OLS: slope={rts["slope"]:.2f} (theory: 4.0)') cx, cy = rts['log_theta'].mean(), rts['log_diff'].mean() y_ref = 4.0 * (x_fit - cx) + cy ax.plot(x_fit, y_ref, 'k--', linewidth=1, alpha=0.5, label='Reference slope=4') ax.set_xlabel(r'$\log(\theta)$ (Fisher-Rao distance)') ax.set_ylabel(r'$\log(\|slerp - linear\|^2)$') ax.set_title(f'Rolling $\\theta$-Scaling (N={len(rts["thetas"])} windows)\n' f'slope={rts["slope"]:.2f} vs theory=4.0') ax.legend() paths += save_figure(fig, 'natural_gradient_rolling_theta_scaling', subfolder='curvature_roles') # Figure 6: Share volatility scaling (T10) — THE KEY FIGURE vs = _cache.get('volatility_scaling') if vs is not None: fig, ax = setup_figure(width=9, height=7) ax.scatter(vs['log_mean'], vs['log_std'], alpha=0.6, color='steelblue', edgecolors='navy', s=50, zorder=3) x_fit = np.linspace(vs['log_mean'].min(), vs['log_mean'].max(), 50) # OLS fit y_fit = vs['slope'] * x_fit + vs['intercept'] ax.plot(x_fit, y_fit, 'r-', linewidth=2.5, label=f'Data: slope={vs["slope"]:.3f} (SE={vs["se"]:.3f})') # Reference lines through centroid cx, cy = vs['log_mean'].mean(), vs['log_std'].mean() # Fisher-Rao (slope=0.5) y_fr = 0.5 * (x_fit - cx) + cy ax.plot(x_fit, y_fr, 'g--', linewidth=1.5, alpha=0.7, label='Fisher-Rao: slope=0.5') # Euclidean (slope=0) y_euc = 0.0 * (x_fit - cx) + cy ax.plot(x_fit, y_euc, 'k:', linewidth=1.5, alpha=0.5, label='Euclidean: slope=0') # Log-normal (slope=1) y_ln = 1.0 * (x_fit - cx) + cy ax.plot(x_fit, y_ln, 'm-.', linewidth=1.5, alpha=0.5, label='Log-normal: slope=1') ax.set_xlabel(r'$\log(\bar{s}_i)$ (mean industry share)') ax.set_ylabel(r'$\log(\mathrm{std}(\Delta s_i))$ (share change volatility)') ax.set_title(f'Share Volatility Scaling: Fisher-Rao Metric Validation\n' f'slope={vs["slope"]:.3f}, 95% CI [{vs["ci_low"]:.3f}, {vs["ci_high"]:.3f}], ' f'R$^2$={vs["r_squared"]:.3f}, N={vs["n"]}') ax.legend(fontsize=9, loc='upper left') ax.grid(True, alpha=0.3) paths += save_figure(fig, 'natural_gradient_volatility_scaling', subfolder='curvature_roles') # Figure 7: Geodesic momentum (T11) — informative negative gm = _cache.get('geodesic_momentum') if gm is not None: fig, (ax1, ax2) = setup_subplots(nrows=1, ncols=2, width=14, height=5) # Left: Slerp vs Linear RSS per quarter (scatter) n_q = len(gm['slerp_arr']) ax1.scatter(gm['linear_arr'] * 1e6, gm['slerp_arr'] * 1e6, alpha=0.4, color='steelblue', s=20, edgecolors='navy', linewidth=0.3) lim = max(gm['linear_arr'].max(), gm['slerp_arr'].max()) * 1e6 * 1.05 ax1.plot([0, lim], [0, lim], 'k--', alpha=0.5, label='Equal') ax1.set_xlabel(r'Linear extrapolation RSS ($\times 10^{-6}$)') ax1.set_ylabel(r'Slerp extrapolation RSS ($\times 10^{-6}$)') n_beats = int(np.sum(gm['slerp_arr'] < gm['linear_arr'])) ax1.set_title(f'Slerp vs Linear Extrapolation\n' f'Slerp wins {n_beats}/{n_q} quarters (ratio={gm["ratio_sl"]:.3f})') ax1.legend(fontsize=9) # Right: Both vs no-change benchmark ax2.boxplot([gm['slerp_arr'] / gm['nochange_arr'], gm['linear_arr'] / gm['nochange_arr']], labels=['Slerp / No-change', 'Linear / No-change'], patch_artist=True, boxprops=dict(facecolor='steelblue', alpha=0.7)) ax2.axhline(1.0, color='red', linestyle='--', alpha=0.7, label='No-change baseline') ax2.set_ylabel('RSS ratio vs no-change') ax2.set_title(f'Both Extrapolations vs Random Walk\n' f'Slerp/NC={gm["ratio_snc"]:.2f} (>1 = worse)') ax2.legend(fontsize=9) fig.suptitle('Geodesic Momentum Test: No Inertia in Share Dynamics', fontsize=13) fig.tight_layout() paths += save_figure(fig, 'natural_gradient_geodesic_momentum', subfolder='curvature_roles') return paths if __name__ == '__main__': results = main() for r in results: p = r['p_value'] p_str = f'{p:.4f}' if not (p is None or (isinstance(p, float) and np.isnan(p))) else 'NaN' print(f" {r['name']}: {r['verdict']} (stat={r['statistic']:.4f}, p={p_str})")

total	2
p value	0.2500
statistic	1
slerp wins	2

Natural Gradient Midpoint Concentration

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps