#!/usr/bin/env python3 """ Tests for Crisis Duration Scaling Tests the two-term wave duration decomposition from Paper 4 (eq:wave_duration): tau_wave = DeltaPhi/eps_drift + pi/sqrt(|a|*eps_drift) where the first term is the drift phase and the second is the fold passage (canard bifurcation delay, Theorem thm:canard). Sub-tests: 1. Wave duration vs alpha (Spearman, N=12) 2. Two-term NLS: tau = c1*DeltaPhi/alpha + c2*pi/sqrt(alpha) 3. FRED IP recession recovery by sector (EXPLORATORY) 4. Combined assessment Geometric dimension d assignments follow GEOMETRIC_DIM_BY_PRODUCT in tests/applications/test_tech_cycle_predictions.py (d=2 for planar/lithographic processes like semiconductors, solar PV, NAND; d=1 for linear/assembly processes). The d=1 vs d=2 duration comparison is handled by test_tech_cycle_predictions.py (tech_cycle_geometric_dimension subtest), not duplicated here. Review Notes (2026-03-11, v3): Theory match: CORRECT for Test 2. The two-term decomposition directly tests eq:wave_duration from Paper 4. The NLS fit estimates both drift (c1) and fold passage (c2) coefficients, with c2 > 0 confirming the canard contribution. Lean axiom: transition_duration_scaling (CESProofs/Hierarchy/TransitionDynamics.lean). I/O audit: CORRECT. alpha (input) from published experience curves, wave duration (output) from historical records, DeltaPhi = ln(cost_reduction_ratio) from wave-specific cost data. Data-signal alignment: ADEQUATE for Tests 1-2 (N=12 technologies with literature alpha values and wave-specific cost reductions). WEAK for Test 3 (FRED IP recession recovery maps structural prediction to cyclical dynamics). Confidence in verdict: MEDIUM for Test 2 (structural model), LOW for others. """ 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_compiled from common.stats import log_result, spearman_test from common.plotting import setup_figure, setup_subplots, save_figure # NBER recession peaks and troughs RECESSIONS = [ ("1973-11", "1975-03"), ("1981-07", "1982-11"), ("1990-07", "1991-03"), ("2001-03", "2001-11"), ("2007-12", "2009-06"), ("2020-02", "2020-04"), ] # FRED IP series -> sector info # Geometric dimension d follows GEOMETRIC_DIM_BY_PRODUCT in # tests/applications/test_tech_cycle_predictions.py: # d=2 for planar/lithographic (semiconductors), d=1 for linear/assembly (all others) # Alpha = d * alpha_0 where alpha_0 = 0.12 (Wright 1936; Nagy et al. 2013). ALPHA_0 = 0.12 FRED_SECTORS = { "IPG331S": {"name": "Primary metals", "d": 1, "alpha": 1 * ALPHA_0}, "IPG332S": {"name": "Fabricated metals", "d": 1, "alpha": 1 * ALPHA_0}, "IPG333S": {"name": "Machinery", "d": 1, "alpha": 1 * ALPHA_0}, "IPG334S": {"name": "Computer & electronic", "d": 2, "alpha": 2 * ALPHA_0}, "IPG335S": {"name": "Electrical equipment", "d": 1, "alpha": 1 * ALPHA_0}, "IPG3361T3S": {"name": "Motor vehicles", "d": 1, "alpha": 1 * ALPHA_0}, "IPG3364T9S": {"name": "Other transportation", "d": 1, "alpha": 1 * ALPHA_0}, "IPG337S": {"name": "Furniture", "d": 1, "alpha": 1 * ALPHA_0}, "IPG339S": {"name": "Misc manufacturing", "d": 1, "alpha": 1 * ALPHA_0}, "IPG325S": {"name": "Chemicals", "d": 1, "alpha": 1 * ALPHA_0}, "IPG322S": {"name": "Paper", "d": 1, "alpha": 1 * ALPHA_0}, "IPG324S": {"name": "Petroleum & coal", "d": 1, "alpha": 1 * ALPHA_0}, "IPG326S": {"name": "Plastics & rubber", "d": 1, "alpha": 1 * ALPHA_0}, } def _date_to_months(date_str): parts = date_str.split("-") return int(parts[0]) * 12 + int(parts[1]) def _measure_recession_duration(ip_data, peak_date, trough_date): """Measure peak-to-recovery duration for a sector during a recession.""" peak_m = _date_to_months(peak_date) trough_m = _date_to_months(trough_date) # Find IP at recession peak peak_ip = None for date_str, val in ip_data: m = _date_to_months(date_str) if abs(m - peak_m) <= 1: peak_ip = val break if peak_ip is None or peak_ip <= 0: return None # Find trough IP trough_ip = float('inf') for date_str, val in ip_data: m = _date_to_months(date_str) if peak_m <= m <= trough_m + 6: if val < trough_ip: trough_ip = val if trough_ip == float('inf'): return None depth_pct = (peak_ip - trough_ip) / peak_ip # Find recovery recovery_m = None for date_str, val in ip_data: m = _date_to_months(date_str) if m > trough_m and val >= peak_ip: recovery_m = m break if recovery_m is None: last_date, last_val = ip_data[-1] recovery_m = _date_to_months(last_date) if last_val < peak_ip: return None p2t = trough_m - peak_m p2r = recovery_m - peak_m return (p2t, p2r, depth_pct) def _load_ip_as_tuples(series_id): """Load a FRED IP series as list of (date_str, value) tuples.""" try: s = load_fred(series_id) return [(d.strftime("%Y-%m-%d"), v) for d, v in s.items()] except FileNotFoundError: return None _cache = {} def main(): results = [] # Load wave data (expanded to N=12 with cost_reduction_ratio for DeltaPhi) wave_data = None try: wave_data = load_compiled("wave_learning_rates") except FileNotFoundError: pass # --- Test 1: Wave duration vs alpha (Spearman) --- if wave_data is not None: data_points = [] names_list = [] for _, r in wave_data.iterrows(): alpha_val = r.get("alpha_lit", None) dur = r.get("wave_duration_yr", None) if alpha_val is not None and dur is not None: alpha_val = float(alpha_val) dur = float(dur) if alpha_val > 0 and dur > 0: data_points.append((alpha_val, dur)) names_list.append(str(r.get('era', ''))) if len(data_points) >= 3: alphas = [d[0] for d in data_points] durs = [d[1] for d in data_points] rho_s, p_s = spearman_test(alphas, durs) if rho_s < -0.5 and p_s < 0.10: verd = "CONSISTENT" elif rho_s < 0: verd = "AMBIGUOUS" else: verd = "INCONSISTENT" _cache['wave_vs_alpha'] = { 'alphas': alphas, 'durs': durs, 'rho': rho_s, 'p': p_s, 'names': names_list } results.append(log_result("crisis_duration_wave_vs_alpha", rho_s, p_s, verd, {"n_waves": len(data_points)}, description="Tests whether industries with higher learning rates (alpha) " "have shorter technology wave durations. Spearman correlation between " "literature alpha values and wave duration in years from compiled wave " f"data (N={len(data_points)}). Prediction: rho_s < 0.")) else: results.append(log_result("crisis_duration_wave_vs_alpha", np.nan, np.nan, "AMBIGUOUS", {"reason": "insufficient data"}, description="Tests wave duration vs learning rate alpha. Insufficient data.")) else: results.append(log_result("crisis_duration_wave_vs_alpha", np.nan, np.nan, "AMBIGUOUS", {"reason": "wave_learning_rates.csv not found"}, description="Tests wave duration vs learning rate alpha. Data not found.")) # --- Test 2: Two-term NLS model (KEY TEST) --- # Theory (Paper 4, eq:wave_duration): # tau_wave = DeltaPhi/eps_drift + pi/sqrt(|a|*eps_drift) # where DeltaPhi = ln(c0/c*) is the potential of discovery and eps_drift ~ alpha. # We fit: tau = c1 * (DeltaPhi/alpha) + c2 * (pi/sqrt(alpha)) # and test: c1 > 0 (drift present), c2 > 0 (fold passage present), R^2 reasonable. if wave_data is not None: nls_points = [] for _, r in wave_data.iterrows(): alpha_val = r.get("alpha_lit", None) dur = r.get("wave_duration_yr", None) cr = r.get("cost_reduction_ratio", None) if alpha_val is not None and dur is not None and cr is not None: alpha_val = float(alpha_val) dur = float(dur) cr = float(cr) if alpha_val > 0 and dur > 0 and cr > 1: dphi = np.log(cr) nls_points.append({ 'era': str(r.get('era', '')), 'alpha': alpha_val, 'dur': dur, 'dphi': dphi, 'drift': dphi / alpha_val, 'fold': np.pi / np.sqrt(alpha_val) }) if len(nls_points) >= 5: from scipy.optimize import curve_fit durs_arr = np.array([p['dur'] for p in nls_points]) drift_arr = np.array([p['drift'] for p in nls_points]) fold_arr = np.array([p['fold'] for p in nls_points]) def two_term(x, c1, c2): return c1 * x[0] + c2 * x[1] X = np.vstack([drift_arr, fold_arr]) try: popt, pcov = curve_fit(two_term, X, durs_arr, p0=[0.5, 1.0]) c1, c2 = popt perr = np.sqrt(np.diag(pcov)) residuals = durs_arr - two_term(X, *popt) ss_res = np.sum(residuals**2) ss_tot = np.sum((durs_arr - durs_arr.mean())**2) r_sq = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan rmse = np.sqrt(np.mean(residuals**2)) # t-statistics for coefficients t_c1 = c1 / perr[0] if perr[0] > 0 else np.nan t_c2 = c2 / perr[1] if perr[1] > 0 else np.nan # Verdict: key question is whether both duration components exist. # c1 > 0 means drift phase (DeltaPhi/alpha) contributes. # c2 > 0 means fold passage (pi/sqrt(alpha)) contributes — this # is the canard bifurcation prediction under test. both_positive = c1 > 0 and c2 > 0 c2_significant = t_c2 > 2.0 # fold passage significant at ~5% c1_significant = t_c1 > 2.0 # drift also significant if both_positive and c2_significant: verd = "CONSISTENT" elif c2 <= 0: verd = "INCONSISTENT" else: verd = "AMBIGUOUS" predicted = two_term(X, *popt) _cache['two_term_nls'] = { 'points': nls_points, 'c1': c1, 'c2': c2, 'c1_se': perr[0], 'c2_se': perr[1], 'r_sq': r_sq, 'rmse': rmse, 'predicted': predicted.tolist(), 'actual': durs_arr.tolist(), 'residuals': residuals.tolist(), } results.append(log_result("crisis_duration_two_term_nls", r_sq, np.nan, verd, {"c1_drift": round(c1, 3), "c1_se": round(perr[0], 3), "c2_fold": round(c2, 3), "c2_se": round(perr[1], 3), "t_c1": round(t_c1, 2), "t_c2": round(t_c2, 2), "r_squared": round(r_sq, 3), "rmse_years": round(rmse, 1), "n_technologies": len(nls_points)}, description="Tests the two-term wave duration decomposition from Paper 4 " "eq:wave_duration: tau = c1*(DeltaPhi/alpha) + c2*(pi/sqrt(alpha)). " "DeltaPhi = ln(cost_reduction_ratio) is the potential of discovery; " "alpha is the Wright's Law learning rate. NLS fit on N compiled " "technology waves. Prediction: c1 > 0 (drift phase present), " "c2 > 0 (fold passage / canard bifurcation delay present). " "c2 > 0 is the key test of the canard contribution.")) except Exception as e: results.append(log_result("crisis_duration_two_term_nls", np.nan, np.nan, "AMBIGUOUS", {"reason": f"NLS fit failed: {e}"}, description="Two-term NLS fit of wave duration decomposition. Fit failed.")) else: results.append(log_result("crisis_duration_two_term_nls", np.nan, np.nan, "AMBIGUOUS", {"reason": f"need >= 5 technologies with cost data, have {len(nls_points)}"}, description="Two-term NLS fit. Insufficient technologies with cost reduction data.")) else: results.append(log_result("crisis_duration_two_term_nls", np.nan, np.nan, "AMBIGUOUS", {"reason": "wave_learning_rates.csv not found"}, description="Two-term NLS fit. Data not found.")) # --- Test 3: FRED IP recession recovery (EXPLORATORY) --- # Maps structural transition prediction to cyclical recession recovery. # Capped at AMBIGUOUS because the canard bifurcation theory predicts # structural technology transition durations, not business cycle recovery. sector_avg_p2r = {} for series_id, info in FRED_SECTORS.items(): ip_data = _load_ip_as_tuples(series_id) if ip_data is None: continue durations = [] for peak, trough in RECESSIONS: result = _measure_recession_duration(ip_data, peak, trough) if result is not None: durations.append(result[1]) # peak-to-recovery if durations: sector_avg_p2r[series_id] = np.mean(durations) if len(sector_avg_p2r) >= 4: alphas = [FRED_SECTORS[s]["alpha"] for s in sector_avg_p2r] p2rs = list(sector_avg_p2r.values()) n_distinct = len(set(alphas)) if n_distinct < 3: d1_p2r = [p for s, p in sector_avg_p2r.items() if FRED_SECTORS[s]["d"] == 1] d2_p2r = [p for s, p in sector_avg_p2r.items() if FRED_SECTORS[s]["d"] == 2] if d1_p2r and d2_p2r: mean_d1 = np.mean(d1_p2r) mean_d2 = np.mean(d2_p2r) _cache['ip_vs_alpha'] = {'d1_p2r': d1_p2r, 'd2_p2r': d2_p2r, 'mean_d1': mean_d1, 'mean_d2': mean_d2} verd = "AMBIGUOUS" # capped regardless of direction results.append(log_result("crisis_duration_ip_recovery", mean_d1 - mean_d2, np.nan, verd, {"n_sectors": len(sector_avg_p2r), "n_distinct_alpha": n_distinct, "mean_d1_p2r": round(mean_d1, 1), "mean_d2_p2r": round(mean_d2, 1), "exploratory": True}, description="Exploratory: compares mean recession recovery duration " "for d=1 vs d=2 manufacturing sectors in FRED IP data. " "d assignments follow GEOMETRIC_DIM_BY_PRODUCT in " "test_tech_cycle_predictions.py. Capped at AMBIGUOUS because the canard " "bifurcation prediction applies to structural technology transitions " "(decades), not cyclical recession recovery (months).")) else: results.append(log_result("crisis_duration_ip_recovery", np.nan, np.nan, "AMBIGUOUS", {"reason": "need both d=1 and d=2", "exploratory": True}, description="Exploratory: recession recovery vs learning rate. " "Need both d=1 and d=2 sectors.")) else: rho_s, p_s = spearman_test(alphas, p2rs) verd = "AMBIGUOUS" # capped — exploratory results.append(log_result("crisis_duration_ip_recovery", rho_s, p_s, verd, {"n_sectors": len(sector_avg_p2r), "n_distinct_alpha": n_distinct, "exploratory": True}, description="Exploratory: Spearman(alpha, peak-to-recovery) across FRED IP sectors. " "Capped at AMBIGUOUS.")) else: results.append(log_result("crisis_duration_ip_recovery", np.nan, np.nan, "AMBIGUOUS", {"reason": "insufficient FRED data", "exploratory": True}, description="Exploratory: recession duration vs learning rate. Insufficient data.")) # --- Test 4: Combined --- core_results = [r for r in results if not r.get("exploratory", False)] n_core_consistent = sum(1 for r in core_results if r["verdict"] == "CONSISTENT") n_core_total = len(core_results) n_total = len(results) n_consistent = sum(1 for r in results if r["verdict"] == "CONSISTENT") n_inconsistent = sum(1 for r in results if r["verdict"] == "INCONSISTENT") if n_core_consistent >= 2: combined_verd = "CONSISTENT" elif n_inconsistent > 0: combined_verd = "INCONSISTENT" elif n_core_consistent >= 1: combined_verd = "AMBIGUOUS" else: combined_verd = "AMBIGUOUS" results.append(log_result("crisis_duration_combined", float(n_consistent), np.nan, combined_verd, {"n_consistent": n_consistent, "n_total": n_total, "n_core_consistent": n_core_consistent, "n_core_total": n_core_total}, description=f"Combined verdict across crisis duration sub-tests. " f"Core tests (wave_vs_alpha, two_term_nls) directly test the canard " f"bifurcation duration prediction. {n_consistent} of {n_total} CONSISTENT.")) _cache['results'] = results return results def figures(): """Generate figures for each crisis duration subtest.""" if not _cache: main() paths = [] # --- Figure 1: Wave duration vs alpha scatter --- data = _cache.get('wave_vs_alpha') if data: fig, ax = setup_figure() ax.scatter(data['alphas'], data['durs'], s=80, zorder=5) for i, name in enumerate(data.get('names', [])): if i < len(data['alphas']): ax.annotate(name, (data['alphas'][i], data['durs'][i]), textcoords="offset points", xytext=(5, 5), fontsize=7) ax.set_xlabel('Learning rate (alpha)') ax.set_ylabel('Wave duration (years)') ax.set_title(f'Wave Duration vs Learning Rate\n' f'Spearman rho={data["rho"]:.3f}, p={data["p"]:.3f}, N={len(data["alphas"])}') paths.extend(save_figure(fig, 'crisis_duration_wave_vs_alpha', subfolder='dynamics')) # --- Figure 2: Two-term NLS predicted vs actual --- data = _cache.get('two_term_nls') if data: fig, axes = setup_subplots(1, 2, width=14, height=5) # Panel A: Predicted vs Actual ax = axes[0] actual = data['actual'] predicted = data['predicted'] ax.scatter(actual, predicted, s=80, zorder=5, edgecolors='black') for i, pt in enumerate(data['points']): ax.annotate(pt['era'], (actual[i], predicted[i]), textcoords="offset points", xytext=(5, 5), fontsize=7) lims = [min(min(actual), min(predicted)) - 2, max(max(actual), max(predicted)) + 2] ax.plot(lims, lims, 'k--', alpha=0.5, label='Perfect fit') ax.set_xlim(lims) ax.set_ylim(lims) ax.set_xlabel('Actual wave duration (years)') ax.set_ylabel('Predicted: c1*(DeltaPhi/alpha) + c2*(pi/sqrt(alpha))') ax.set_title(f'Two-Term NLS: R²={data["r_sq"]:.3f}, RMSE={data["rmse"]:.1f}yr') ax.legend() # Panel B: Coefficient decomposition (stacked bar) ax = axes[1] eras = [pt['era'][:15] for pt in data['points']] drift_contrib = [data['c1'] * pt['drift'] for pt in data['points']] fold_contrib = [data['c2'] * pt['fold'] for pt in data['points']] x = range(len(eras)) ax.bar(x, drift_contrib, label=f'Drift: c1={data["c1"]:.2f}', color='#4878CF', alpha=0.7) ax.bar(x, fold_contrib, bottom=drift_contrib, label=f'Fold passage: c2={data["c2"]:.2f}', color='#E24A33', alpha=0.7) # Overlay actual durations ax.scatter(x, [pt['dur'] for pt in data['points']], marker='D', color='black', s=50, zorder=5, label='Actual') ax.set_xticks(x) ax.set_xticklabels(eras, rotation=45, ha='right', fontsize=7) ax.set_ylabel('Duration (years)') ax.set_title('Duration Decomposition: Drift + Fold Passage') ax.legend(fontsize=8) fig.tight_layout() paths.extend(save_figure(fig, 'crisis_duration_two_term_nls', subfolder='dynamics')) # --- Figure 3: IP recession recovery box plot --- data = _cache.get('ip_vs_alpha') if data: fig, ax = setup_figure(width=6, height=5) bp = ax.boxplot([data['d1_p2r'], data['d2_p2r']], labels=['d=1 (linear)', 'd=2 (planar)'], patch_artist=True) bp['boxes'][0].set_facecolor('#4878CF') bp['boxes'][1].set_facecolor('#E24A33') for box in bp['boxes']: box.set_alpha(0.5) ax.set_ylabel('Peak-to-recovery (months)') ax.set_title(f'Recession Recovery Duration by Geometric Dimension\n' f'd=1 mean={data["mean_d1"]:.1f}mo, d=2 mean={data["mean_d2"]:.1f}mo') paths.extend(save_figure(fig, 'crisis_duration_ip_recovery', subfolder='dynamics')) # --- Figure 4: Combined verdict summary --- res = _cache.get('results', []) if res: fig, ax = setup_figure(width=7, height=4) names = [r['name'].replace('crisis_duration_', '') for r in res] verdicts = [r['verdict'] for r in res] colors = {'CONSISTENT': '#2CA02C', 'AMBIGUOUS': '#FF7F0E', 'INCONSISTENT': '#D62728'} ax.barh(range(len(names)), [1]*len(names), color=[colors.get(v, '#999') for v in verdicts]) ax.set_yticks(range(len(names))) ax.set_yticklabels(names, fontsize=9) ax.set_xlim(0, 1.5) ax.set_xticks([]) for i, v in enumerate(verdicts): ax.text(0.5, i, v, ha='center', va='center', fontweight='bold', fontsize=10) ax.set_title('Crisis Duration: Subtest Verdicts') ax.invert_yaxis() paths.extend(save_figure(fig, 'crisis_duration_combined', 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']})")

mean d1	28
mean d2	26.7500
predicted	1.4142
statistic	1.0467

Crisis Duration D1 D2 Ratio

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps