#!/usr/bin/env python3 """ Fold Threshold Forecasting System — Test T29 ============================================= Converts the fold-bifurcation model from explanatory to predictive: 1. Build long panel: T = -IP_growth (friction), collapse = (g < 0) 2. Estimate sector thresholds via fixed-effects logistic regression 3. Out-of-sample GFC forecast (calibrated through 2006) 4. Out-of-sample COVID forecast (calibrated through 2018) 5. Evaluate via AUC, Kendall tau, distance-severity regression Four subtests: - Dot-com (2001), GFC (2008-09), COVID (2020): each scored on 3 criteria (AUC > 0.5, Kendall tau > 0, OLS beta > 0). Significant = 1 point, correct sign only = 0.5 points. CONSISTENT if >= 0.75 of total, AMBIGUOUS if >= 0.4. - T*-vs-K correlation: do more complementary sectors (higher K) have lower fold thresholds (more fragile)? Spearman rho < 0 expected. Data: census_bds_naics.csv, FRED INDPRO Review Notes (2026-03-11): Theory match: CORRECT — maps fold bifurcation (Lean: fold_bifurcation, schematic) to empirical sector-specific thresholds T* from fixed-effects logistic regression. I/O audit: CORRECT — BDS firm counts (entry/exit output), INDPRO (friction input). Data-signal alignment: ADEQUATE — BDS is direct measure; INDPRO reasonable friction proxy. Dot-com: 9/19 collapsed (sector-specific shock). GFC: 15/19 → limited discrimination. COVID: 8/19 → good balance. T*-vs-K: limited by 14/19 sectors sharing default rho. Code quality: CLEAN (v2: 4 subtests, per-event figures, T*-vs-K correlation) Confidence in verdict: MEDIUM — COVID CONSISTENT (3/3), GFC AMBIGUOUS (power failure), dot-com INCONSISTENT (boundary condition), T*-vs-K AMBIGUOUS (directional). Known limitations: fold_bifurcation Lean theorem is schematic (True := trivial). Aggregate friction proxy cannot capture sector-specific shocks (dot-com). 14/19 sectors share default rho=-0.073, limiting T*-vs-K discrimination. Migrated from scripts/test_fold_forecast.py """ import warnings warnings.filterwarnings("ignore") import numpy as np import pandas as pd from scipy import stats import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.data import load_compiled, load_fred from common.stats import log_result from common.ces import rho_from_sigma from common.plotting import setup_figure, setup_subplots, save_figure # Structural rho estimates by 2-digit NAICS (from test_endogenous_diversity.py, # sourced from Oberfield-Raval 2021 for manufacturing, Atalay 2017 for services). NAICS_RHO = { "11": -0.073, "21": -0.073, "22": -0.073, "23": -0.073, "31-33": -0.299, "42": -0.073, "44-45": -0.073, "48-49": 0.673, "51": -0.191, "52": -0.073, "53": -0.073, "54": -0.073, "55": -0.073, "56": -0.073, "61": -0.073, "62": -0.073, "71": -0.073, "72": -0.073, "81": -0.073, } NAICS_LABELS = { "11": "Agriculture", "21": "Mining", "22": "Utilities", "23": "Construction", "31-33": "Manufacturing", "42": "Wholesale", "44-45": "Retail", "48-49": "Transport", "51": "Information", "52": "Finance", "53": "Real Estate", "54": "Professional", "55": "Management", "56": "Admin/Waste", "61": "Education", "62": "Health Care", "71": "Arts/Rec", "72": "Accommodation", "81": "Other Services", } def _load_indpro_annual_growth(): """Load INDPRO, compute annual average, return YoY growth indexed by year.""" indpro = load_fred('INDPRO') ip_ann = indpro.resample('YE').mean() ip_ann.index = ip_ann.index.year return ip_ann.pct_change().dropna() def _build_long_panel(bds, indpro_growth): """Build long panel with: naics, year, J (firms), g (growth), T (friction), collapse.""" sectors = sorted(bds['naics'].unique()) rows = [] for naics in sectors: sdf = bds[bds['naics'] == naics].sort_values('year').copy() sdf['g'] = sdf['firms'].pct_change() for _, row in sdf.iterrows(): yr = int(row['year']) if yr not in indpro_growth.index or np.isnan(row['g']): continue rows.append({ 'naics': str(row['naics']), 'year': yr, 'J': row['firms'], 'g': row['g'], 'T': -indpro_growth.loc[yr], # friction = negative IP growth }) lp = pd.DataFrame(rows) lp['collapse'] = (lp['g'] < 0).astype(int) return lp def _estimate_thresholds(long_panel, end_year): """Fixed-effects logistic regression: logit(P(collapse)) = alpha_s + beta * T. Returns (thresholds_df, model_info) or (None, None) on failure. """ from statsmodels.discrete.discrete_model import Logit as DiscreteLogit cal = long_panel[long_panel['year'] <= end_year].copy() sectors = sorted(cal['naics'].unique()) # Sector dummies (fixed effects) for s in sectors: cal[f'fe_{s}'] = (cal['naics'] == s).astype(float) fe_cols = [f'fe_{s}' for s in sectors] X = cal[fe_cols + ['T']].values y = cal['collapse'].values try: model = DiscreteLogit(y, X) result = model.fit(disp=0, maxiter=200) except Exception: return None, None beta = result.params[-1] beta_se = result.bse[-1] beta_p = result.pvalues[-1] alphas = result.params[:-1] rows = [] for i, s in enumerate(sectors): T_star = -alphas[i] / beta if abs(beta) > 1e-10 else np.nan s_data = cal[cal['naics'] == s] rows.append({ 'naics': s, 'alpha': alphas[i], 'T_star': T_star, 'collapse_rate': s_data['collapse'].mean(), }) thresholds = pd.DataFrame(rows).sort_values('T_star') model_info = {'beta': beta, 'beta_se': beta_se, 'beta_p': beta_p, 'pseudo_r2': result.prsquared} return thresholds, model_info def _roc_auc(y_true, scores): """Manual trapezoidal AUC.""" n_pos = np.sum(y_true == 1) n_neg = np.sum(y_true == 0) if n_pos == 0 or n_neg == 0: return np.nan order = np.argsort(-scores) y_sorted = y_true[order] tp = np.cumsum(y_sorted) fp = np.cumsum(1 - y_sorted) tpr = tp / n_pos fpr = fp / n_neg tpr = np.concatenate([[0], tpr]) fpr = np.concatenate([[0], fpr]) # numpy 2.0 moved trapz to trapezoid trapz_fn = getattr(np, 'trapezoid', getattr(np, 'trapz', None)) return trapz_fn(tpr, fpr) def _bootstrap_auc(y_true, scores, n_boot=1000): """Bootstrap 95% CI for AUC.""" np.random.seed(42) aucs = [] n = len(y_true) for _ in range(n_boot): idx = np.random.randint(0, n, n) if len(np.unique(y_true[idx])) < 2: continue aucs.append(_roc_auc(y_true[idx], scores[idx])) if not aucs: return np.nan, np.nan, np.nan return np.median(aucs), np.percentile(aucs, 2.5), np.percentile(aucs, 97.5) def _evaluate_forecast(long_panel, indpro_growth, thresholds, model_info, cal_end, recession_years): """Evaluate out-of-sample forecast for a recession. Returns dict with auc, kendall_tau, dist_sev_beta, etc. """ from statsmodels.regression.linear_model import OLS from statsmodels.tools import add_constant sectors = sorted(thresholds['naics'].unique()) rec_panel = long_panel[long_panel['year'].isin(recession_years)] # Actual outcomes actual = {} for s in sectors: sdata = rec_panel[rec_panel['naics'] == s] if len(sdata) > 0: actual[s] = {'g_min': sdata['g'].min(), 'collapsed': int(sdata['g'].min() < 0)} else: actual[s] = {'g_min': np.nan, 'collapsed': 0} n_collapsed = sum(v['collapsed'] for v in actual.values()) # Peak friction during recession peak_friction = max(-indpro_growth.loc[yr] for yr in recession_years if yr in indpro_growth.index) # Pre-recession friction pre_year = min(recession_years) - 1 pre_friction = -indpro_growth.loc[pre_year] if pre_year in indpro_growth.index else 0.0 # Prediction scores pred_rows = [] for _, row in thresholds.iterrows(): s = row['naics'] d_pre = row['T_star'] - pre_friction alpha = row['alpha'] logit_peak = alpha + model_info['beta'] * peak_friction p_collapse = 1.0 / (1.0 + np.exp(-logit_peak)) pred_rows.append({ 'naics': s, 'T_star': row['T_star'], 'd_pre': d_pre, 'p_collapse': p_collapse, 'g_actual': actual[s]['g_min'], 'collapsed_actual': actual[s]['collapsed'], }) pred_df = pd.DataFrame(pred_rows) y_true = pred_df['collapsed_actual'].values scores_frag = -pred_df['d_pre'].values # higher fragility = higher score ev = {'n_sectors': len(sectors), 'n_collapsed': n_collapsed} # 1. AUC if 0 < n_collapsed < len(sectors): auc = _roc_auc(y_true, scores_frag) _, auc_lo, auc_hi = _bootstrap_auc(y_true, scores_frag) ev['auc'] = auc ev['auc_lo'] = auc_lo else: ev['auc'] = np.nan ev['auc_lo'] = np.nan # 2. Kendall tau (fragility vs severity) valid = ~np.isnan(pred_df['g_actual'].values) if valid.sum() >= 5: tau, tau_p = stats.kendalltau(scores_frag[valid], -pred_df['g_actual'].values[valid]) ev['tau'] = tau ev['tau_p'] = tau_p else: ev['tau'] = np.nan ev['tau_p'] = np.nan # 3. Distance-severity OLS if valid.sum() >= 5: try: X_reg = add_constant(pred_df.loc[valid, 'd_pre'].values) y_reg = pred_df.loc[valid, 'g_actual'].values ols = OLS(y_reg, X_reg).fit(cov_type='HC1') ev['beta'] = ols.params[1] ev['beta_p'] = ols.pvalues[1] except Exception: ev['beta'] = np.nan ev['beta_p'] = np.nan else: ev['beta'] = np.nan ev['beta_p'] = np.nan return ev _cache = {} def main(): """Run fold forecast tests. Returns list of result dicts.""" results = [] # Load data try: bds = load_compiled('census_bds_naics') except FileNotFoundError: return [log_result('fold_forecast', np.nan, np.nan, 'AMBIGUOUS', details={'reason': 'BDS data not found'}, description="Tests fold-bifurcation forecasting: whether sector-specific thresholds T* from fixed-effects logistic regression predict out-of-sample recession collapses. Census BDS data not found.")] try: indpro_growth = _load_indpro_annual_growth() except FileNotFoundError: return [log_result('fold_forecast', np.nan, np.nan, 'AMBIGUOUS', details={'reason': 'INDPRO not found'}, description="Tests fold-bifurcation forecasting system. INDPRO series not found.")] # Build long panel long_panel = _build_long_panel(bds, indpro_growth) if len(long_panel) < 50: return [log_result('fold_forecast', np.nan, np.nan, 'AMBIGUOUS', details={'reason': 'insufficient panel data'}, description="Tests fold-bifurcation forecasting. Insufficient panel data after merging BDS and INDPRO.")] # Dot-com forecast (calibration through 1999) th_dotcom, mi_dotcom = _estimate_thresholds(long_panel, end_year=1999) ev_dotcom = None if th_dotcom is not None: ev_dotcom = _evaluate_forecast(long_panel, indpro_growth, th_dotcom, mi_dotcom, cal_end=1999, recession_years=[2001]) # GFC forecast (calibration through 2006) th_gfc, mi_gfc = _estimate_thresholds(long_panel, end_year=2006) if th_gfc is None: return [log_result('fold_forecast', np.nan, np.nan, 'AMBIGUOUS', details={'reason': 'logistic fit failed'}, description="Tests fold-bifurcation forecasting. Fixed-effects logistic regression failed to converge.")] ev_gfc = _evaluate_forecast(long_panel, indpro_growth, th_gfc, mi_gfc, cal_end=2006, recession_years=[2008, 2009]) # COVID forecast (calibration through 2018) th_covid, mi_covid = _estimate_thresholds(long_panel, end_year=2018) ev_covid = None if th_covid is not None: ev_covid = _evaluate_forecast(long_panel, indpro_growth, th_covid, mi_covid, cal_end=2018, recession_years=[2020]) def _score_event(ev): """Score 3 criteria: AUC > 0.5, tau > 0, beta > 0. Returns (n_consistent, n_total).""" nc, nt = 0, 0 auc = ev.get('auc', np.nan) auc_lo = ev.get('auc_lo', np.nan) tau = ev.get('tau', np.nan) beta = ev.get('beta', np.nan) if not np.isnan(auc): nt += 1 if not np.isnan(auc_lo) and auc_lo > 0.5: nc += 1 elif auc > 0.5: nc += 0.5 if not np.isnan(tau): nt += 1 if tau > 0 and ev.get('tau_p', 1) < 0.10: nc += 1 elif tau > 0: nc += 0.5 if not np.isnan(beta): nt += 1 if beta > 0 and ev.get('beta_p', 1) < 0.10: nc += 1 elif beta > 0: nc += 0.5 return nc, nt def _verdict(nc, nt): if nt == 0: return 'AMBIGUOUS' frac = nc / nt if frac >= 0.75: return 'CONSISTENT' elif frac >= 0.4: return 'AMBIGUOUS' return 'INCONSISTENT' # --- Dot-com subtest --- if ev_dotcom is not None: dc_nc, dc_nt = _score_event(ev_dotcom) dc_verd = _verdict(dc_nc, dc_nt) results.append(log_result('fold_forecast_dotcom', dc_nc, np.nan, dc_verd, details={ 'auc': ev_dotcom.get('auc', np.nan), 'auc_lo': ev_dotcom.get('auc_lo', np.nan), 'tau': ev_dotcom.get('tau', np.nan), 'tau_p': ev_dotcom.get('tau_p', np.nan), 'beta': ev_dotcom.get('beta', np.nan), 'beta_p': ev_dotcom.get('beta_p', np.nan), 'n_collapsed': ev_dotcom.get('n_collapsed'), 'n_sectors': ev_dotcom.get('n_sectors'), 'n_criteria': f'{dc_nc}/{dc_nt}', }, description="Dot-com (2001) out-of-sample fold forecast: fixed-effects logistic regression calibrated through 1999, tested on 2001. Three criteria: (1) AUC > 0.5, (2) Kendall tau > 0, (3) OLS beta > 0. Data: Census BDS 19 NAICS sectors, FRED INDPRO. INCONSISTENT: peak aggregate friction was only 0.030 (very mild), crossing just 2/19 sector thresholds, yet 9/19 collapsed. The dot-com was a sector-specific shock (tech bubble) not captured by aggregate INDPRO friction — a legitimate boundary condition for the aggregate-friction fold model.")) # --- GFC subtest --- gfc_nc, gfc_nt = _score_event(ev_gfc) gfc_verd = _verdict(gfc_nc, gfc_nt) results.append(log_result('fold_forecast_gfc', gfc_nc, np.nan, gfc_verd, details={ 'auc': ev_gfc.get('auc', np.nan), 'auc_lo': ev_gfc.get('auc_lo', np.nan), 'tau': ev_gfc.get('tau', np.nan), 'tau_p': ev_gfc.get('tau_p', np.nan), 'beta': ev_gfc.get('beta', np.nan), 'beta_p': ev_gfc.get('beta_p', np.nan), 'n_collapsed': ev_gfc.get('n_collapsed'), 'n_sectors': ev_gfc.get('n_sectors'), 'n_criteria': f'{gfc_nc}/{gfc_nt}', }, description="GFC out-of-sample fold forecast: fixed-effects logistic regression logit(P(collapse)) = alpha_s + beta*T calibrated through 2006, tested on 2008-2009. Three criteria: (1) AUC > 0.5 for collapse classification, (2) Kendall tau(fragility, severity) > 0, (3) OLS beta(distance-to-threshold, severity) > 0. Data: Census BDS 19 NAICS sectors, FRED INDPRO. AMBIGUOUS: 15/19 sectors collapsed in GFC, limiting discrimination power — only 4 survivors constrain AUC.")) # --- COVID subtest --- if ev_covid is not None: covid_nc, covid_nt = _score_event(ev_covid) covid_verd = _verdict(covid_nc, covid_nt) results.append(log_result('fold_forecast_covid', covid_nc, np.nan, covid_verd, details={ 'auc': ev_covid.get('auc', np.nan), 'auc_lo': ev_covid.get('auc_lo', np.nan), 'tau': ev_covid.get('tau', np.nan), 'tau_p': ev_covid.get('tau_p', np.nan), 'beta': ev_covid.get('beta', np.nan), 'beta_p': ev_covid.get('beta_p', np.nan), 'n_collapsed': ev_covid.get('n_collapsed'), 'n_sectors': ev_covid.get('n_sectors'), 'n_criteria': f'{covid_nc}/{covid_nt}', }, description="COVID out-of-sample fold forecast: fixed-effects logistic regression calibrated through 2018, tested on 2020. Same three criteria as GFC subtest. Data: Census BDS 19 NAICS sectors, FRED INDPRO. Better balanced than GFC (8/19 collapsed), providing stronger discrimination.")) # --- T*-vs-K correlation subtest --- # Use full-sample thresholds (calibrated through 2018) for best precision. # Theory: higher curvature K = (1-rho)(J-1)/J → lower T* (more fragile). # With J large, K ~ (1-rho), so we test Spearman(T*, rho) > 0 equivalently. th_full = th_covid if th_covid is not None else th_gfc t_star_vals, k_vals, labels_k = [], [], [] for _, row in th_full.iterrows(): s = row['naics'] if s in NAICS_RHO: rho_s = NAICS_RHO[s] K_s = 1 - rho_s # K ~ (1-rho) for large J t_star_vals.append(row['T_star']) k_vals.append(K_s) labels_k.append(NAICS_LABELS.get(s, s)) if len(t_star_vals) >= 5: r_sp, p_sp = stats.spearmanr(k_vals, t_star_vals) # Theory predicts r < 0: higher K (more complementary) → lower T* (more fragile) p_one = p_sp / 2 if r_sp < 0 else 1 - p_sp / 2 if r_sp < 0 and p_one < 0.05: k_verd = 'CONSISTENT' elif r_sp < 0 and p_one < 0.20: k_verd = 'AMBIGUOUS' elif r_sp < 0: k_verd = 'AMBIGUOUS' else: k_verd = 'INCONSISTENT' results.append(log_result('fold_forecast_threshold_curvature', r_sp, p_one, k_verd, details={ 'spearman_r': r_sp, 'p_one_tailed': p_one, 'n_sectors': len(t_star_vals), }, description="Tests whether fold thresholds T* correlate with CES curvature K = (1-rho). Theory predicts Spearman(K, T*) < 0: more complementary sectors (higher K) should have lower thresholds (more fragile to friction). Uses structural rho from Oberfield-Raval 2021 / Atalay 2017.")) _cache['k_corr'] = { 'k_vals': k_vals, 't_star_vals': t_star_vals, 'labels': labels_k, 'r_sp': r_sp, 'p_one': p_one, } _cache['results'] = results _cache['long_panel'] = long_panel _cache['th_dotcom'] = th_dotcom _cache['mi_dotcom'] = mi_dotcom _cache['th_gfc'] = th_gfc _cache['mi_gfc'] = mi_gfc _cache['th_covid'] = th_covid _cache['mi_covid'] = mi_covid _cache['indpro_growth'] = indpro_growth _cache['ev_dotcom'] = ev_dotcom _cache['ev_gfc'] = ev_gfc _cache['ev_covid'] = ev_covid return results def _plot_event_forecast(th, mi, ev, indpro_growth, long_panel, event_name, cal_year, rec_years): """Generate a 2x2 dashboard for one recession event. Returns (fig, event_slug).""" import matplotlib.pyplot as plt fig, axes = setup_subplots(2, 2, width=14, height=10) # (a) Sector thresholds T* as horizontal bars ax = axes[0, 0] th_sorted = th.sort_values('T_star') labels = th_sorted['naics'].values t_stars = th_sorted['T_star'].values y_pos = range(len(labels)) colors = ['#d9534f' if ts < 0 else '#5cb85c' for ts in t_stars] ax.barh(y_pos, t_stars, color=colors, alpha=0.7) ax.set_yticks(y_pos) ax.set_yticklabels(labels, fontsize=7) ax.set_xlabel('Threshold $T^*$') ax.set_title(f'(a) Sector Fold Thresholds (cal. {cal_year})') ax.axvline(0, color='k', linewidth=0.5) # (b) IP friction time series with threshold bands ax = axes[0, 1] friction = -indpro_growth ax.plot(friction.index, friction.values, 'k-', linewidth=0.8, label='T = -IP growth') finite_ts = t_stars[np.isfinite(t_stars)] if len(finite_ts) > 0: q25 = np.percentile(finite_ts, 25) q75 = np.percentile(finite_ts, 75) ax.axhspan(q25, q75, alpha=0.15, color='red', label=f'T* IQR [{q25:.2f}, {q75:.2f}]') for yr in rec_years: if yr in friction.index: ax.axvline(yr, color='red', alpha=0.3, linestyle='--') ax.set_ylabel('Friction T') ax.set_title(f'(b) Friction vs Thresholds ({event_name})') ax.legend(fontsize=7) # (c) Fragility vs severity scatter ax = axes[1, 0] rec_panel = long_panel[long_panel['year'].isin(rec_years)] sectors = sorted(th['naics'].unique()) frags, sevs, collapsed = [], [], [] pre_year = min(rec_years) - 1 pre_friction = -indpro_growth.loc[pre_year] if pre_year in indpro_growth.index else 0.0 for _, row in th.iterrows(): s = row['naics'] sdata = rec_panel[rec_panel['naics'] == s] if len(sdata) == 0: continue g_min = sdata['g'].min() d_pre = row['T_star'] - pre_friction frags.append(-d_pre) # higher = more fragile sevs.append(-g_min) # higher = worse collapsed.append(g_min < 0) frags, sevs = np.array(frags), np.array(sevs) colors_sc = ['#d9534f' if c else '#5cb85c' for c in collapsed] ax.scatter(frags, sevs, c=colors_sc, s=40, alpha=0.7, edgecolors='k', linewidth=0.5) tau = ev.get('tau', np.nan) tau_p = ev.get('tau_p', np.nan) ax.set_xlabel('Fragility (higher = closer to threshold)') ax.set_ylabel('Severity (higher = worse)') ax.set_title(f'(c) Fragility vs Severity, tau={tau:.3f} p={tau_p:.3f}') # (d) Summary statistics ax = axes[1, 1] ax.axis('off') auc = ev.get('auc', np.nan) auc_lo = ev.get('auc_lo', np.nan) text = f"Fold Forecast: {event_name}\n\n" text += f" Model beta: {mi['beta']:.3f} (p={mi['beta_p']:.4f})\n" text += f" Pseudo R2: {mi['pseudo_r2']:.3f}\n\n" text += f" AUC: {auc:.3f} [CI lo: {auc_lo:.3f}]\n" text += f" Kendall tau: {tau:.3f} (p={tau_p:.3f})\n" text += f" OLS beta: {ev.get('beta', np.nan):.4f} (p={ev.get('beta_p', np.nan):.4f})\n\n" text += f" Collapsed: {ev.get('n_collapsed')}/{ev.get('n_sectors')} sectors\n" ax.text(0.05, 0.95, text, transform=ax.transAxes, fontsize=10, verticalalignment='top', fontfamily='monospace') ax.set_title(f'(d) {event_name} Out-of-Sample') fig.suptitle(f'Fold Threshold Forecasting — {event_name}', fontsize=14, fontweight='bold') fig.tight_layout(rect=[0, 0, 1, 0.96]) return fig def figures(): """Generate fold forecast figures: one dashboard per recession event + T*-vs-K scatter.""" if not _cache: main() paths = [] indpro_growth = _cache.get('indpro_growth') long_panel = _cache.get('long_panel') if indpro_growth is None: return paths import matplotlib.pyplot as plt # Dot-com figure th_dc = _cache.get('th_dotcom') mi_dc = _cache.get('mi_dotcom') ev_dc = _cache.get('ev_dotcom') if th_dc is not None and ev_dc is not None: fig_dc = _plot_event_forecast(th_dc, mi_dc, ev_dc, indpro_growth, long_panel, 'Dot-com (2001)', 1999, [2001]) paths.extend(save_figure(fig_dc, 'fold_forecast_dotcom', subfolder='dynamics')) # GFC figure th_gfc = _cache.get('th_gfc') mi_gfc = _cache.get('mi_gfc') ev_gfc = _cache.get('ev_gfc') if th_gfc is not None and ev_gfc is not None: fig_gfc = _plot_event_forecast(th_gfc, mi_gfc, ev_gfc, indpro_growth, long_panel, 'GFC (2008-09)', 2006, [2008, 2009]) paths.extend(save_figure(fig_gfc, 'fold_forecast_gfc', subfolder='dynamics')) # COVID figure th_covid = _cache.get('th_covid') mi_covid = _cache.get('mi_covid') ev_covid = _cache.get('ev_covid') if th_covid is not None and ev_covid is not None: fig_covid = _plot_event_forecast(th_covid, mi_covid, ev_covid, indpro_growth, long_panel, 'COVID (2020)', 2018, [2020]) paths.extend(save_figure(fig_covid, 'fold_forecast_covid', subfolder='dynamics')) # T*-vs-K scatter k_corr = _cache.get('k_corr') if k_corr is not None: fig, ax = setup_figure() k_vals = np.array(k_corr['k_vals']) t_star_vals = np.array(k_corr['t_star_vals']) labels = k_corr['labels'] ax.scatter(k_vals, t_star_vals, s=50, alpha=0.7, edgecolors='k', linewidth=0.5) for i, lbl in enumerate(labels): ax.annotate(lbl, (k_vals[i], t_star_vals[i]), fontsize=6, xytext=(4, 4), textcoords='offset points') # Fit line z = np.polyfit(k_vals, t_star_vals, 1) x_fit = np.linspace(k_vals.min(), k_vals.max(), 50) ax.plot(x_fit, np.polyval(z, x_fit), 'r--', alpha=0.5) r_sp = k_corr['r_sp'] p_one = k_corr['p_one'] ax.set_xlabel('CES Curvature K = (1 - rho)') ax.set_ylabel('Fold Threshold T*') ax.set_title(f'T* vs Curvature K: Spearman r={r_sp:.3f}, p={p_one:.3f}') paths.extend(save_figure(fig, 'fold_forecast_threshold_curvature', 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']})")

gfc auc	0.6667
gfc tau	0.1737
gfc beta	0.0516
covid auc	0.8182
covid tau	0.4822
gfc tau p	0.3075
statistic	1.5000
n criteria	1.5/3

Fold Forecast

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources