#!/usr/bin/env python3 """ Test: Horse Race and Out-of-Sample Forecasting — Paper 5 (Empirical Tests) HORSE RACE TESTS (HR1-HR3): HR1: Damping cancellation — CES transient decay vs "delayed positive effect" Fit exponential decay vs hump-shaped IRF to Basel III local projections HR2: Cycle amplitude — CES sigma vs sector size vs trend growth rate HR3: Timescale hierarchy — CES r*=2.19 vs Marchenko-Pastur random matrix null OUT-OF-SAMPLE FORECASTS (OOS1-OOS3): OOS1: Wright's Law solar PV — train on 1976-2015, forecast 2016-2023 OOS2: Damping cancellation — pre-Basel III IRF predicts post-Basel III non-effect OOS3: Semiconductor recovery — pre-2020 recession patterns predict 2020 recovery Data: BCL regulation panel, IMF FDI, FRED IP indices, solar PV (hardcoded). Review Notes (2026-03-12) — OOS2 (Basel III) evaluation ----------------------------------------------------------------- OOS2 Theory match: CORRECT — Paper 4 Proposition 6.1 damping cancellation. Regulatory tightening has transient but not persistent effect on financial development. Lean: damping_cancellation_algebraic (REAL PROOF, field_simp). OOS2 I/O audit: CORRECT — CSI change as treatment, FDI change as outcome. OOS2 Data-signal alignment: ADEQUATE — N=101 (48 treated, 53 control). Treatment threshold (CSI increase >= 2) captures 48 countries vs paper's 25 Basel III adopters. No pre-trend check. Two-sample t-test vs paper's full panel DID with country FE. OOS2 Code quality: MINOR ISSUES — pre-IRF uses shift(-h) (row-shift) instead of year-based lookup; doesn't affect verdict. No _cache writes or figure (FIXED 2026-03-12). OOS2 Confidence: MEDIUM — p=0.37 supports null (no persistent effect), consistent with paper's p=0.95. Different methodology explains p-value gap. Review Notes (2026-03-11) — focused on OOS3 (covid_recovery) ----------------------------------------------------------------- OOS3 Theory match: PARTIAL — alpha = d*alpha_0 predicts TECHNOLOGY TRANSITION duration (decades), not recession recovery (months). OOS3 I/O audit: CORRECT for forecast, MIXED for d comparison (N=1 for d=2) OOS3 Data-signal alignment: ADEQUATE (forecast) / WEAK (d comparison — N=1, survivorship bias) OOS3 Code quality: CLEAN after fixes (cache writes, figure, restructured verdict) OOS3 Confidence: MEDIUM — OOS forecast passes (-39% error), d comparison is auxiliary Known limitations: - d=2 has only one sector (IPG334S) — cross-sectional d comparison has no power - 2001 and 2020 recessions specifically targeted tech sector (confounders) - Non-recovering sectors excluded from d comparison (survivorship bias) """ import warnings warnings.filterwarnings("ignore") import numpy as np import pandas as pd import csv import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.config import COMPILED_DIR, FRED_DIR from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure try: from scipy import stats from scipy.optimize import curve_fit HAS_SCIPY = True except ImportError: HAS_SCIPY = False try: import statsmodels.api as sm HAS_SM = True except ImportError: HAS_SM = False # --- Data loading utilities --- def _load_csv_dicts(fname): path = COMPILED_DIR / fname if not path.exists(): return None with open(path) as f: reader = csv.DictReader(f) rows = [] for row in reader: d = {} for k, v in row.items(): try: d[k] = float(v.replace(",", "")) except (ValueError, AttributeError): d[k] = v rows.append(d) return rows def _load_fred_ip(series_id): """Load FRED IP monthly series -> [(date_str, value), ...].""" path = FRED_DIR / f"{series_id}.csv" if not path.exists(): return None data = [] with open(path) as f: reader = csv.reader(f) next(reader) for row in reader: if len(row) >= 2 and row[1].strip(): try: data.append((row[0].strip(), float(row[1]))) except ValueError: pass return data def _date_to_months(date_str): parts = date_str.split("-") return int(parts[0]) * 12 + int(parts[1]) 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_SECTORS = { "IPG331S": {"name": "Primary metals", "d": 1, "alpha": 0.12}, "IPG332S": {"name": "Fabricated metals", "d": 1, "alpha": 0.12}, "IPG333S": {"name": "Machinery", "d": 1, "alpha": 0.14}, "IPG334S": {"name": "Computer & electronic", "d": 2, "alpha": 0.24}, "IPG335S": {"name": "Electrical equipment", "d": 1.5, "alpha": 0.18}, "IPG3361T3S": {"name": "Motor vehicles", "d": 1, "alpha": 0.14}, "IPG3364T9S": {"name": "Other transportation", "d": 1, "alpha": 0.18}, "IPG337S": {"name": "Furniture", "d": 1, "alpha": 0.12}, "IPG339S": {"name": "Misc manufacturing", "d": 1, "alpha": 0.12}, "IPG325S": {"name": "Chemicals", "d": 1, "alpha": 0.13}, "IPG322S": {"name": "Paper", "d": 1, "alpha": 0.12}, "IPG324S": {"name": "Petroleum & coal", "d": 1, "alpha": 0.12}, "IPG326S": {"name": "Plastics & rubber", "d": 1, "alpha": 0.13}, } PREDICTED_MODEL = { 'capital_stringency_idx': 'A', 'activity_restrictions_idx': 'B', 'supervisory_power_idx': 'B', 'overall_restrictiveness_idx': 'B', } REGULATION_TYPE = { 'capital_stringency_idx': 'direct stock', 'activity_restrictions_idx': 'institutional', 'supervisory_power_idx': 'institutional', 'overall_restrictiveness_idx': 'composite/institutional', } OR_SECTORS = { '336': ('Transport Equip', 0.50, 'IPG336S'), '331': ('Primary Metals', 0.55, 'IPG331S'), '327': ('Nonmetallic Min', 0.60, 'IPG327S'), '333': ('Machinery', 0.65, 'IPG333S'), '321': ('Wood', 0.70, 'IPG321S'), '332': ('Fabricated Metals', 0.70, 'IPG332S'), '335': ('Electrical Equip', 0.70, 'IPG335S'), '322': ('Paper', 0.85, 'IPG322S'), '337': ('Furniture', 0.85, 'IPG337S'), '325': ('Chemicals', 0.90, 'IPG325S'), '326': ('Plastics/Rubber', 0.95, 'IPG326S'), '339': ('Miscellaneous', 1.10, 'IPG339S'), '311': ('Food', 1.15, 'IPG311S'), '334': ('Computer/Elec', 1.40, 'IPG334S'), } def _measure_recession(ip_data, peak_date, trough_date): peak_m = _date_to_months(peak_date) trough_m = _date_to_months(trough_date) peak_ip = None for ds, val in ip_data: m = _date_to_months(ds) if abs(m - peak_m) <= 1: peak_ip = val break if peak_ip is None or peak_ip <= 0: return None trough_ip = float('inf') for ds, val in ip_data: m = _date_to_months(ds) if peak_m <= m <= trough_m + 6: if val < trough_ip: trough_ip = val if trough_ip == float('inf'): return None depth = (peak_ip - trough_ip) / peak_ip recovery_m = None for ds, val in ip_data: m = _date_to_months(ds) if m > trough_m and val >= peak_ip: recovery_m = m break if recovery_m is None: return None p2r = recovery_m - peak_m pre_vals = [] for ds, val in ip_data: m = _date_to_months(ds) if peak_m - 36 <= m < peak_m: pre_vals.append(val) volatility = np.std(np.diff(pre_vals)) / np.mean(pre_vals) if len(pre_vals) > 6 else np.nan return {'p2r': p2r, 'depth': depth, 'peak_ip': peak_ip, 'volatility': volatility} # ====================================================================== # HR1: Damping cancellation horse race # ====================================================================== _ALL_SUBTESTS = [ 'horse_race_oos_irf_stratification', 'horse_race_oos_sigma_amplitude', 'horse_race_oos_eigenvalue_spacing', 'horse_race_oos_wrights_law', 'horse_race_oos_basel_iii', 'horse_race_oos_covid_recovery', ] _cache = {} def _hr1_damping_horse_race(): results = [] if not HAS_SM or not HAS_SCIPY: results.append(log_result( "horse_race_oos_damping", np.nan, np.nan, "AMBIGUOUS", details={"reason": "requires statsmodels and scipy"}, description="Skipped: statsmodels and scipy required for local projection IRF and curve fitting." )) return results bcl_path = COMPILED_DIR / 'bank_regulation_full_panel.csv' fdi_path = COMPILED_DIR / 'imf_financial_development.csv' if not bcl_path.exists() or not fdi_path.exists(): results.append(log_result( "horse_race_oos_damping", np.nan, np.nan, "AMBIGUOUS", details={"reason": "BCL or FDI data not found"}, description="Skipped: BCL regulation panel or IMF Financial Development data not available." )) return results bcl = pd.read_csv(bcl_path) fdi = pd.read_csv(fdi_path) # Prepare BCL shocks bcl_cols = ['country_code', 'year', 'wave', 'capital_stringency_idx', 'activity_restrictions_idx', 'supervisory_power_idx', 'overall_restrictiveness_idx'] bcl_cols = [c for c in bcl_cols if c in bcl.columns] bcl_sub = bcl[bcl_cols].copy().sort_values(['country_code', 'year']) for col in ['capital_stringency_idx', 'activity_restrictions_idx', 'supervisory_power_idx', 'overall_restrictiveness_idx']: if col in bcl_sub.columns: bcl_sub[f'd_{col}'] = bcl_sub.groupby('country_code')[col].diff() # FDI lookup fdi_col = None for c in ['financial_development_overall', 'fd_index', 'value']: if c in fdi.columns: fdi_col = c break if fdi_col is None: results.append(log_result( "horse_race_oos_damping", np.nan, np.nan, "AMBIGUOUS", details={"reason": "Cannot find FDI column"}, description="Skipped: financial development outcome column not found in IMF dataset." )) return results cc_col = 'country_code' if 'country_code' in fdi.columns else 'iso2' fdi_sub = fdi[[cc_col, 'year', fdi_col]].copy() fdi_sub.columns = ['country_code', 'year', 'fdi'] fdi_sub['year'] = fdi_sub['year'].astype(int) fdi_lookup = {} for _, row in fdi_sub.iterrows(): fdi_lookup[(row['country_code'], int(row['year']))] = row['fdi'] bcl_sub['fdi_base'] = bcl_sub.apply( lambda r: fdi_lookup.get((r['country_code'], int(r['year'])), np.nan), axis=1) bcl_sub = bcl_sub.dropna(subset=['fdi_base']) shock_vars = [c for c in bcl_sub.columns if c.startswith('d_') and 'idx' in c] n_correct = 0 n_shocks = 0 for shock in shock_vars: shock_name = shock.replace('d_', '').replace('_idx', '').replace('_', ' ').title() horizons = list(range(9)) betas = [] ses = [] for h in horizons: df_h = bcl_sub.copy() df_h['fdi_forward'] = df_h.apply( lambda r: fdi_lookup.get((r['country_code'], int(r['year']) + h), np.nan), axis=1) df_h['dfdi'] = df_h['fdi_forward'] - df_h['fdi_base'] df_h = df_h.dropna(subset=[shock, 'dfdi', 'fdi_base']) if len(df_h) < 30: betas.append(np.nan) ses.append(np.nan) continue for col in [shock, 'dfdi', 'fdi_base']: df_h[col] = df_h[col] - df_h.groupby('country_code')[col].transform('mean') X = sm.add_constant(df_h[[shock, 'fdi_base']]) y = df_h['dfdi'] try: model = sm.OLS(y, X).fit(cov_type='cluster', cov_kwds={'groups': bcl_sub.loc[df_h.index, 'country_code']}) betas.append(model.params[shock]) ses.append(model.bse[shock]) except Exception: betas.append(np.nan) ses.append(np.nan) betas = np.array(betas, dtype=float) ses = np.array(ses, dtype=float) valid = ~np.isnan(betas) h_valid = np.array(horizons)[valid] b_valid = betas[valid] # Fit Model A: monotone decay model_a_aic = np.inf try: def decay_model(h, a, lam): return a * np.exp(-lam * h) popt_a, _ = curve_fit(decay_model, h_valid, b_valid, p0=[b_valid[0] if len(b_valid) > 0 else 0, 0.3], maxfev=5000) resid_a = b_valid - decay_model(h_valid, *popt_a) ss_a = np.sum(resid_a**2) n = len(b_valid) model_a_aic = n * np.log(ss_a / n + 1e-20) + 2 * 2 except Exception: pass # Fit Model B: hump-shaped model_b_aic = np.inf try: def hump_model(h, a, lam): return a * h * np.exp(-lam * h) popt_b, _ = curve_fit(hump_model, h_valid, b_valid, p0=[0.01, 0.5], maxfev=5000) resid_b = b_valid - hump_model(h_valid, *popt_b) ss_b = np.sum(resid_b**2) n = len(b_valid) model_b_aic = n * np.log(ss_b / n + 1e-20) + 2 * 2 except Exception: pass # Fit Model C: biexponential cascade model_c_aic = np.inf try: def biexp_model(h, a, lam1, lam2): return a * (np.exp(-lam1 * h) - np.exp(-lam2 * h)) popt_c, _ = curve_fit(biexp_model, h_valid, b_valid, p0=[0.02, 0.3, 0.8], bounds=([0, 1e-4, 1e-4], [np.inf, 10, 10]), maxfev=10000) if popt_c[2] > popt_c[1]: resid_c = b_valid - biexp_model(h_valid, *popt_c) ss_c = np.sum(resid_c**2) n = len(b_valid) model_c_aic = n * np.log(ss_c / n + 1e-20) + 2 * 3 except Exception: pass aic_map = {'A': model_a_aic, 'B': model_b_aic, 'C': model_c_aic} best_model = min(aic_map, key=aic_map.get) shock_key = shock.replace('d_', '') predicted = PREDICTED_MODEL.get(shock_key, '?') match = (best_model == predicted) or (predicted == 'B' and best_model == 'C') if match: n_correct += 1 n_shocks += 1 # Cache per-shock IRF data for figures _cache.setdefault('hr1_shocks', []).append({ 'shock_name': shock_name, 'horizons': h_valid.tolist(), 'betas': b_valid.tolist(), 'ses': ses[valid].tolist(), 'best_model': best_model, 'predicted': predicted, 'match': match, 'aic_map': {k: float(v) for k, v in aic_map.items()}, 'popt_a': popt_a.tolist() if model_a_aic < np.inf else None, 'popt_b': popt_b.tolist() if model_b_aic < np.inf else None, 'popt_c': popt_c.tolist() if model_c_aic < np.inf else None, }) # Overall verdict based on stratification if n_shocks == 0: verd = "AMBIGUOUS" elif n_correct == n_shocks: verd = "CONSISTENT" elif n_correct >= n_shocks - 1: verd = "AMBIGUOUS" else: verd = "INCONSISTENT" results.append(log_result( "horse_race_oos_irf_stratification", stat=float(n_correct), p=np.nan, verd=verd, details={"n_correct": n_correct, "n_shocks": n_shocks}, description="HR1 horse race: fits three IRF shape models (A: monotone decay, B: hump-shaped, " "C: biexponential) to local projection estimates for each BCL regulatory shock. " "CES predicts Model A for direct-stock shocks (capital stringency) and Model B/C for " "institutional shocks. Consistent if AIC-best model matches prediction for all shocks." )) return results # ====================================================================== # HR2: Cycle amplitude horse race # ====================================================================== def _hr2_cycle_amplitude(): results = [] if not HAS_SCIPY: results.append(log_result( "HR2: Cycle amplitude", np.nan, np.nan, "AMBIGUOUS", details={"reason": "requires scipy"}, description="Skipped: scipy required for Spearman rank correlation in cycle amplitude analysis." )) return results records = [] for naics, (name, sigma, fred_id) in OR_SECTORS.items(): ip_data = _load_fred_ip(fred_id) if ip_data is None: alt_ids = {'IPG336S': 'IPG3361T3S', 'IPG311S': 'IPG311A2S'} if fred_id in alt_ids: ip_data = _load_fred_ip(alt_ids[fred_id]) if ip_data is None: continue values = np.array([v for _, v in ip_data if v > 0]) if len(values) < 120: continue log_growth = np.diff(np.log(values)) amplitude = np.std(log_growth) * np.sqrt(12) mean_level = np.mean(values) trend_growth = 12 * np.mean(log_growth) records.append({ 'naics': naics, 'name': name, 'sigma': sigma, 'amplitude': amplitude, 'mean_level': mean_level, 'trend_growth': trend_growth, }) if len(records) < 8: results.append(log_result( "HR2: Cycle amplitude", np.nan, np.nan, "AMBIGUOUS", details={"reason": f"Only {len(records)} matched sectors"}, description="Skipped: insufficient FRED manufacturing sector IP series available for amplitude analysis." )) return results sigma = np.array([r['sigma'] for r in records]) amplitude = np.array([r['amplitude'] for r in records]) size = np.log(np.array([r['mean_level'] for r in records])) growth = np.array([r['trend_growth'] for r in records]) _cache['hr2_records'] = records _cache['hr2_sigma'] = sigma _cache['hr2_amplitude'] = amplitude r_sigma, p_sigma = stats.spearmanr(sigma, amplitude) r_size, p_size = stats.spearmanr(size, amplitude) r_growth, p_growth = stats.spearmanr(growth, amplitude) ces_strongest = abs(r_sigma) >= max(abs(r_size), abs(r_growth)) ces_correct_sign = r_sigma < 0 if ces_correct_sign and ces_strongest and p_sigma < 0.10: verd = "CONSISTENT" elif ces_correct_sign and p_sigma < 0.10: verd = "AMBIGUOUS" elif ces_correct_sign: verd = "AMBIGUOUS" else: verd = "INCONSISTENT" results.append(log_result( "horse_race_oos_sigma_amplitude", stat=float(r_sigma), p=float(p_sigma), verd=verd, details={"r_size": float(r_size), "p_size": float(p_size), "r_growth": float(r_growth), "p_growth": float(p_growth)}, description="HR2 horse race: tests whether CES sigma is the strongest predictor of annualized cycle " "amplitude across 14 manufacturing sectors, competing against sector size and trend growth. " "Uses Spearman rank correlations on FRED IP data. Consistent if |r_sigma| > |r_size|, |r_growth| and r_sigma < 0." )) return results # ====================================================================== # HR3: Timescale Marchenko-Pastur null # ====================================================================== def _hr3_timescale_mp(): results = [] sector_keys = ["IPG334S", "IPG335S", "IPG333S", "IPG325S", "IPG336S", "IPG331S", "IPG311A2S"] series_data = {} for key in sector_keys: ip = _load_fred_ip(key) if ip is not None: series_data[key] = {d: v for d, v in ip} if len(series_data) < 5: results.append(log_result( "horse_race_oos_eigenvalue_spacing", np.nan, np.nan, "AMBIGUOUS", details={"reason": f"Only {len(series_data)} sectors available"}, description="Skipped: insufficient FRED sector IP series for eigenvalue spacing analysis." )) return results all_dates = None for key, data in series_data.items(): dates = set(data.keys()) all_dates = dates if all_dates is None else all_dates & dates common_dates = sorted(all_dates) keys = list(series_data.keys()) growth_matrix = [] for i in range(1, len(common_dates)): row = [] for key in keys: v0 = series_data[key].get(common_dates[i-1]) v1 = series_data[key].get(common_dates[i]) if v0 and v1 and v0 > 0 and v1 > 0: row.append(np.log(v1 / v0)) else: row.append(np.nan) growth_matrix.append(row) G = np.array(growth_matrix) mask = ~np.any(np.isnan(G), axis=1) G = G[mask] n, p = G.shape if n < 50: results.append(log_result( "horse_race_oos_eigenvalue_spacing", np.nan, np.nan, "AMBIGUOUS", details={"reason": "Insufficient observations"}, description="Skipped: too few common monthly observations across sectors for correlation eigenvalue analysis." )) return results corr = np.corrcoef(G.T) eigvals = np.sort(np.linalg.eigvalsh(corr))[::-1] sig_eig = eigvals[eigvals > 0.1] if len(sig_eig) < 2: results.append(log_result( "horse_race_oos_eigenvalue_spacing", np.nan, np.nan, "AMBIGUOUS", details={"reason": "Fewer than 2 significant eigenvalues"}, description="Skipped: fewer than 2 correlation eigenvalues exceed 0.1; cannot compute spacing ratios." )) return results obs_ratios = sig_eig[:-1] / sig_eig[1:] obs_median_ratio = float(np.median(obs_ratios)) # Monte Carlo np.random.seed(42) n_mc = 5000 mc_median_ratios = [] for _ in range(n_mc): X_rand = np.random.randn(n, p) C_rand = np.corrcoef(X_rand.T) eig_rand = np.sort(np.linalg.eigvalsh(C_rand))[::-1] sig_rand = eig_rand[eig_rand > 0.1] if len(sig_rand) >= 2: ratios_rand = sig_rand[:-1] / sig_rand[1:] mc_median_ratios.append(np.median(ratios_rand)) mc_median_ratios = np.array(mc_median_ratios) p_median = float(np.mean(mc_median_ratios >= obs_median_ratio)) if p_median < 0.05: verd = "CONSISTENT" elif p_median < 0.10: verd = "AMBIGUOUS" else: verd = "AMBIGUOUS" results.append(log_result( "horse_race_oos_eigenvalue_spacing", stat=obs_median_ratio, p=p_median, verd=verd, details={"mc_mean": float(np.mean(mc_median_ratios))}, description="HR3 Marchenko-Pastur null: tests whether the observed eigenvalue spacing ratio of the " "cross-sector correlation matrix exceeds what random matrices produce. Monte Carlo with " "5000 random matrices matching observed density. Consistent if p < 0.05 (structured hierarchy)." )) return results # ====================================================================== # OOS1: Wright's Law solar PV # ====================================================================== def _oos1_solar_pv(): results = [] if not HAS_SCIPY: results.append(log_result( "horse_race_oos_wrights_law", np.nan, np.nan, "AMBIGUOUS", details={"reason": "requires scipy"}, description="Skipped: scipy required for Wright's Law regression and forecasting." )) return results solar = [ (1976, 76.0, 0.001), (1980, 30.0, 0.01), (1985, 12.0, 0.1), (1990, 8.0, 0.5), (1995, 5.5, 0.8), (2000, 4.0, 1.5), (2005, 3.5, 5.0), (2008, 2.5, 16.0), (2010, 1.5, 40.0), (2012, 0.7, 100.0), (2014, 0.55, 180.0), (2016, 0.35, 300.0), (2018, 0.25, 500.0), (2020, 0.22, 750.0), (2022, 0.26, 1200.0), (2023, 0.22, 1500.0), ] years = np.array([s[0] for s in solar]) costs = np.array([s[1] for s in solar]) cum_gw = np.array([s[2] for s in solar]) train_mask = years <= 2015 test_mask = years > 2015 log_cost_train = np.log(costs[train_mask]) log_Q_train = np.log(cum_gw[train_mask]) log_cost_test = np.log(costs[test_mask]) log_Q_test = np.log(cum_gw[test_mask]) slope, intercept, r_value, p_value, se = stats.linregress(log_Q_train, log_cost_train) beta_train = -slope r2_train = r_value**2 log_cost_pred = intercept + slope * log_Q_test cost_pred = np.exp(log_cost_pred) cost_actual = costs[test_mask] errors = [(cost_pred[i] - cost_actual[i]) / cost_actual[i] * 100 for i in range(len(cost_actual))] mape = np.mean(np.abs(errors)) _cache['oos1_years'] = years _cache['oos1_costs'] = costs _cache['oos1_cum_gw'] = cum_gw _cache['oos1_train_mask'] = train_mask _cache['oos1_test_mask'] = test_mask _cache['oos1_slope'] = slope _cache['oos1_intercept'] = intercept _cache['oos1_cost_pred'] = cost_pred last_train_cost = costs[train_mask][-1] naive_errors = [(last_train_cost - c) / c * 100 for c in cost_actual] naive_mape = np.mean(np.abs(naive_errors)) if mape < 50 and mape < naive_mape: verd = "CONSISTENT" elif mape < 100: verd = "AMBIGUOUS" else: verd = "INCONSISTENT" results.append(log_result( "horse_race_oos_wrights_law", stat=float(mape), p=np.nan, verd=verd, details={"beta_train": float(beta_train), "r2_train": float(r2_train), "naive_mape": float(naive_mape)}, description="OOS1: trains Wright's Law (log-log cost vs cumulative production) on solar PV data " "through 2015, then forecasts 2016-2023. Reports MAPE vs naive (last-observation) forecast. " "Consistent if MAPE < 50% and beats naive, validating learning-curve extrapolation." )) return results # ====================================================================== # OOS2: Damping cancellation out-of-sample # ====================================================================== def _oos2_damping_oos(): results = [] if not HAS_SM or not HAS_SCIPY: results.append(log_result( "OOS2: Damping OOS", np.nan, np.nan, "AMBIGUOUS", details={"reason": "requires statsmodels and scipy"}, description="Skipped: statsmodels and scipy required for pre/post-Basel III damping OOS test." )) return results bcl_path = COMPILED_DIR / 'bank_regulation_full_panel.csv' fdi_path = COMPILED_DIR / 'imf_financial_development.csv' if not bcl_path.exists() or not fdi_path.exists(): results.append(log_result( "OOS2: Damping OOS", np.nan, np.nan, "AMBIGUOUS", details={"reason": "Data not found"}, description="Skipped: BCL regulation or IMF Financial Development data not available for OOS2." )) return results bcl = pd.read_csv(bcl_path) fdi = pd.read_csv(fdi_path) fdi_col = None for c in ['financial_development_overall', 'fd_index', 'value']: if c in fdi.columns: fdi_col = c break cc_col = 'country_code' if 'country_code' in fdi.columns else 'iso2' fdi_sub = fdi[[cc_col, 'year', fdi_col]].copy() fdi_sub.columns = ['country_code', 'year', 'fdi'] fdi_sub['year'] = fdi_sub['year'].astype(int) bcl_pre = bcl[bcl['year'] <= 2007].copy() bcl_pre.sort_values(['country_code', 'year'], inplace=True) shock = 'd_capital_stringency_idx' if 'capital_stringency_idx' in bcl_pre.columns: bcl_pre[shock] = bcl_pre.groupby('country_code')['capital_stringency_idx'].diff() else: results.append(log_result( "OOS2: Damping OOS", np.nan, np.nan, "AMBIGUOUS", details={"reason": "capital_stringency_idx not available"}, description="Skipped: capital stringency index column not available in BCL data for OOS2 test." )) return results pre = bcl_pre.merge(fdi_sub, on=['country_code', 'year'], how='inner') # Build FDI lookup for year-based horizon (not row-shift) fdi_lookup = {} for _, row in fdi_sub.iterrows(): fdi_lookup[(row['country_code'], int(row['year']))] = row['fdi'] pre_betas = [] for h in range(6): df_h = pre.copy() df_h['fdi_forward'] = df_h.apply( lambda r: fdi_lookup.get((r['country_code'], int(r['year']) + h), np.nan), axis=1) df_h['dfdi'] = df_h['fdi_forward'] - df_h['fdi'] df_h = df_h.dropna(subset=[shock, 'dfdi', 'fdi']) if len(df_h) < 20: pre_betas.append(np.nan) continue for col in [shock, 'dfdi', 'fdi']: df_h[col] = df_h[col] - df_h.groupby('country_code')[col].transform('mean') X = sm.add_constant(df_h[[shock, 'fdi']]) y = df_h['dfdi'] try: model = sm.OLS(y, X).fit() pre_betas.append(model.params[shock]) except Exception: pre_betas.append(np.nan) pre_betas = np.array(pre_betas, dtype=float) # Post-period DID did_estimate = np.nan p_did = np.nan if 'capital_stringency_idx' in bcl.columns: w3 = bcl[bcl['year'] == 2007][['country_code', 'capital_stringency_idx']].set_index('country_code') w5 = bcl[bcl['year'] == 2019][['country_code', 'capital_stringency_idx']].set_index('country_code') both = w3.join(w5, lsuffix='_w3', rsuffix='_w5').dropna() both['change'] = both['capital_stringency_idx_w5'] - both['capital_stringency_idx_w3'] treated = set(both[both['change'] >= 2].index) fdi_2007 = fdi_sub[fdi_sub['year'] == 2007].set_index('country_code')['fdi'] fdi_2019 = fdi_sub[fdi_sub['year'] == 2019].set_index('country_code')['fdi'] common = set(fdi_2007.index) & set(fdi_2019.index) & set(both.index) if len(common) > 10: did_data = [] for cc in common: did_data.append({ 'treated': int(cc in treated), 'dfdi': fdi_2019[cc] - fdi_2007[cc] }) did_df = pd.DataFrame(did_data) t_treated = did_df[did_df['treated'] == 1]['dfdi'] t_control = did_df[did_df['treated'] == 0]['dfdi'] did_estimate = float(t_treated.mean() - t_control.mean()) _, p_did = stats.ttest_ind(t_treated, t_control) p_did = float(p_did) if not np.isnan(p_did) and p_did > 0.10: verd = "CONSISTENT" elif not np.isnan(p_did): verd = "AMBIGUOUS" else: verd = "AMBIGUOUS" n_treated = int(t_treated.shape[0]) if 't_treated' in dir() else 0 n_control = int(t_control.shape[0]) if 't_control' in dir() else 0 _cache['oos2_data'] = { 'did_estimate': float(did_estimate) if not np.isnan(did_estimate) else None, 'p_did': float(p_did) if not np.isnan(p_did) else None, 'n_treated': n_treated, 'n_control': n_control, 'treated_mean': float(t_treated.mean()) if n_treated > 0 else None, 'control_mean': float(t_control.mean()) if n_control > 0 else None, 'treated_dfdi': t_treated.values.tolist() if n_treated > 0 else [], 'control_dfdi': t_control.values.tolist() if n_control > 0 else [], 'pre_betas': pre_betas.tolist(), } results.append(log_result( "horse_race_oos_basel_iii", stat=float(did_estimate) if not np.isnan(did_estimate) else np.nan, p=p_did if not np.isnan(p_did) else np.nan, verd=verd, details={"pre_mean_late_beta": float(np.nanmean(pre_betas[3:6])), "n_treated": n_treated, "n_control": n_control}, description="Tests damping cancellation prediction (Paper 4, Prop 6.1): regulatory tightening has " "transient but not persistent effect on financial development. DID: treated = countries with " "BCL capital stringency increase >= 2 between 2007-2019 (Basel III era). Two-sample t-test " "of FDI change (treated vs control). Consistent if p > 0.10 (null of no persistent effect " "cannot be rejected). Pre-2007 local projection IRF provides auxiliary evidence of decay." )) return results # ====================================================================== # OOS3: Semiconductor recovery # ====================================================================== def _oos3_semiconductor_recovery(): results = [] desc = ("OOS3: trains on median peak-to-recovery time from pre-2020 recessions for " "IPG334S (Computer & Electronic) to forecast COVID semiconductor recovery. " "Also reports geometric dimension context: d=2 sectors should recover faster " "than d=1 sectors across all recessions (theory: alpha = d*alpha_0). " "Verdict based on forecast error alone (<50% = CONSISTENT). " "d comparison is auxiliary (N=1 for d=2, confounded by semiconductor shortage in COVID).") target_series = "IPG334S" ip_target = _load_fred_ip(target_series) if ip_target is None: results.append(log_result( "horse_race_oos_covid_recovery", np.nan, np.nan, "SKIPPED", details={"reason": f"{target_series} not found"}, description=desc )) return results pre_2020 = RECESSIONS[:-1] covid = RECESSIONS[-1] train_recoveries = [] for peak, trough in pre_2020: result = _measure_recession(ip_target, peak, trough) if result is not None: train_recoveries.append(result) if len(train_recoveries) < 2: results.append(log_result( "horse_race_oos_covid_recovery", np.nan, np.nan, "SKIPPED", details={"reason": "Insufficient pre-2020 recessions"}, description=desc )) return results p2r_train = [r['p2r'] for r in train_recoveries] median_p2r = np.median(p2r_train) covid_result = _measure_recession(ip_target, covid[0], covid[1]) if covid_result is not None: actual_p2r = covid_result['p2r'] forecast_error = (actual_p2r - median_p2r) / median_p2r * 100 # All-recession d comparison (auxiliary evidence) all_recession_ratios = [] covid_ratio = np.nan all_ip = {} for series_id in FRED_SECTORS: ip = _load_fred_ip(series_id) if ip is not None: all_ip[series_id] = ip for rec_idx, (peak, trough) in enumerate(RECESSIONS): d1_p2r, d2_p2r = [], [] d1_no_recovery, d2_no_recovery = 0, 0 for series_id, info in FRED_SECTORS.items(): if series_id not in all_ip: continue rec = _measure_recession(all_ip[series_id], peak, trough) if rec is not None: if info['d'] == 1: d1_p2r.append(rec['p2r']) elif info['d'] == 2: d2_p2r.append(rec['p2r']) else: if info['d'] == 1: d1_no_recovery += 1 elif info['d'] == 2: d2_no_recovery += 1 if d1_p2r and d2_p2r: ratio = np.mean(d1_p2r) / np.mean(d2_p2r) all_recession_ratios.append({ 'peak': peak, 'ratio': ratio, 'd1_mean': np.mean(d1_p2r), 'd2_mean': np.mean(d2_p2r), 'd1_n': len(d1_p2r), 'd2_n': len(d2_p2r), 'd1_no_recovery': d1_no_recovery, }) if peak == covid[0]: covid_ratio = ratio n_correct = sum(1 for r in all_recession_ratios if r['ratio'] > 1.0) n_total = len(all_recession_ratios) pre_covid_correct = sum(1 for r in all_recession_ratios if r['peak'] != covid[0] and r['ratio'] > 1.0) pre_covid_total = sum(1 for r in all_recession_ratios if r['peak'] != covid[0]) # Verdict based on forecast error alone (the true OOS prediction) if abs(forecast_error) < 50: verd = "CONSISTENT" elif abs(forecast_error) < 100: verd = "AMBIGUOUS" else: verd = "INCONSISTENT" _cache['oos3_data'] = { 'actual_p2r': actual_p2r, 'median_p2r': median_p2r, 'forecast_error': forecast_error, 'p2r_train': p2r_train, 'all_recession_ratios': all_recession_ratios, 'covid_ratio': covid_ratio, } results.append(log_result( "horse_race_oos_covid_recovery", stat=float(actual_p2r), p=np.nan, verd=verd, details={ "forecast_p2r": float(median_p2r), "forecast_error_pct": float(forecast_error), "d1_d2_ratio_covid": float(covid_ratio) if not np.isnan(covid_ratio) else None, "d_correct_all": f"{n_correct}/{n_total}", "d_correct_pre_covid": f"{pre_covid_correct}/{pre_covid_total}", }, description=desc )) else: results.append(log_result( "horse_race_oos_covid_recovery", np.nan, np.nan, "SKIPPED", details={"reason": "COVID recovery not measurable"}, description=desc )) return results def main(): results = [] results.extend(_hr1_damping_horse_race()) results.extend(_hr2_cycle_amplitude()) results.extend(_hr3_timescale_mp()) results.extend(_oos1_solar_pv()) results.extend(_oos2_damping_oos()) results.extend(_oos3_semiconductor_recovery()) return results def figures(): """Generate figures from cached main() results.""" if not _cache: main() paths = [] import matplotlib.pyplot as plt # HR1 figure: IRF curves per shock type with model fits hr1_shocks = _cache.get('hr1_shocks', []) if hr1_shocks: n_shocks = min(len(hr1_shocks), 3) fig, axes = setup_subplots(nrows=1, ncols=n_shocks, width=5*n_shocks, height=4) if n_shocks == 1: axes = [axes] for k, shock_data in enumerate(hr1_shocks[:n_shocks]): ax = axes[k] h = np.array(shock_data['horizons']) b = np.array(shock_data['betas']) se = np.array(shock_data['ses']) ax.errorbar(h, b, yerr=1.96*se, fmt='ko', capsize=3, label='LP estimates') ax.axhline(0, color='gray', linestyle='-', alpha=0.5) h_fine = np.linspace(0, max(h), 100) best = shock_data['best_model'] if shock_data['popt_a'] is not None: a_p = shock_data['popt_a'] y_a = a_p[0] * np.exp(-a_p[1] * h_fine) style = '-' if best == 'A' else '--' ax.plot(h_fine, y_a, f'b{style}', alpha=0.7, label=f'A: decay') if shock_data['popt_b'] is not None: b_p = shock_data['popt_b'] y_b = b_p[0] * h_fine * np.exp(-b_p[1] * h_fine) style = '-' if best == 'B' else '--' ax.plot(h_fine, y_b, f'r{style}', alpha=0.7, label=f'B: hump') if shock_data['popt_c'] is not None: c_p = shock_data['popt_c'] y_c = c_p[0] * (np.exp(-c_p[1] * h_fine) - np.exp(-c_p[2] * h_fine)) style = '-' if best == 'C' else '--' ax.plot(h_fine, y_c, f'g{style}', alpha=0.7, label=f'C: biexp') match_str = 'MATCH' if shock_data['match'] else 'MISMATCH' ax.set_title(f'{shock_data["shock_name"]}\nbest={best} pred={shock_data["predicted"]} [{match_str}]', fontsize=9) ax.set_xlabel('Horizon (waves)') ax.set_ylabel('IRF coefficient') ax.legend(fontsize=7) fig.suptitle('HR1: Damping Cancellation IRF Stratification', fontsize=12, y=1.02) fig.tight_layout() paths.extend(save_figure(fig, 'horse_race_hr1_irf', subfolder='hierarchy')) # HR2 figure: Cycle amplitude vs sigma if 'hr2_sigma' in _cache and 'hr2_amplitude' in _cache: sigma = _cache['hr2_sigma'] amplitude = _cache['hr2_amplitude'] records = _cache['hr2_records'] fig, ax = setup_figure(width=7, height=5) ax.scatter(sigma, amplitude, s=60, edgecolors='black', zorder=3) for rec in records: ax.annotate(rec['name'], (rec['sigma'], rec['amplitude']), fontsize=7, ha='left', va='bottom') z = np.polyfit(sigma, amplitude, 1) x_line = np.linspace(min(sigma), max(sigma), 50) ax.plot(x_line, np.polyval(z, x_line), 'r--', alpha=0.6) ax.set_xlabel('CES substitution elasticity (sigma)') ax.set_ylabel('Annualized cycle amplitude') ax.set_title('HR2: Cycle Amplitude vs Sigma') paths.extend(save_figure(fig, 'horse_race_hr2_amplitude', subfolder='hierarchy')) # OOS1 figure: Wright's Law solar PV if 'oos1_cum_gw' in _cache: fig, ax = setup_figure(width=7, height=5) cum_gw = _cache['oos1_cum_gw'] costs = _cache['oos1_costs'] train_mask = _cache['oos1_train_mask'] test_mask = _cache['oos1_test_mask'] slope = _cache['oos1_slope'] intercept = _cache['oos1_intercept'] cost_pred = _cache['oos1_cost_pred'] ax.scatter(cum_gw[train_mask], costs[train_mask], c='steelblue', s=50, edgecolors='black', label='Train (<=2015)', zorder=3) ax.scatter(cum_gw[test_mask], costs[test_mask], c='orange', s=50, edgecolors='black', label='Test (>2015)', zorder=3) q_line = np.logspace(np.log10(cum_gw.min()), np.log10(cum_gw.max()), 100) ax.plot(q_line, np.exp(intercept + slope * np.log(q_line)), 'r--', label=f'Fit: slope={slope:.2f}') ax.scatter(cum_gw[test_mask], cost_pred, marker='x', c='red', s=40, label='Predicted', zorder=4) ax.set_xscale('log') ax.set_yscale('log') ax.set_xlabel('Cumulative installed capacity (GW)') ax.set_ylabel('Cost ($/W)') ax.set_title("OOS1: Wright's Law Solar PV Forecast") ax.legend(fontsize=8) paths.extend(save_figure(fig, 'horse_race_oos1_solar', subfolder='hierarchy')) # OOS2 figure: Basel III DID — treated vs control FDI change oos2 = _cache.get('oos2_data') if oos2 is not None and oos2.get('treated_dfdi') and oos2.get('control_dfdi'): fig, axes = setup_subplots(nrows=1, ncols=2, width=14, height=5) # Panel A: Box plots of FDI change ax = axes[0] treated_vals = oos2['treated_dfdi'] control_vals = oos2['control_dfdi'] bp = ax.boxplot([treated_vals, control_vals], labels=[f'Treated\n(CSI+≥2, N={oos2["n_treated"]})', f'Control\n(N={oos2["n_control"]})'], patch_artist=True) bp['boxes'][0].set_facecolor('coral') bp['boxes'][1].set_facecolor('steelblue') ax.axhline(0, color='gray', linestyle='--', alpha=0.5) ax.set_ylabel('FDI Change (2007→2019)') p_val = oos2['p_did'] did_est = oos2['did_estimate'] ax.set_title(f'OOS2: Basel III DID\nDID={did_est:.4f}, p={p_val:.3f}') # Panel B: Pre-2007 IRF decay ax = axes[1] pre_betas = oos2.get('pre_betas', []) if pre_betas: horizons = list(range(len(pre_betas))) valid = [(h, b) for h, b in zip(horizons, pre_betas) if not np.isnan(b)] if valid: hs, bs = zip(*valid) ax.bar(hs, bs, color='steelblue', edgecolor='black', width=0.6) ax.axhline(0, color='gray', linestyle='--', alpha=0.5) ax.set_xlabel('Horizon (years)') ax.set_ylabel('IRF coefficient (pre-2007)') ax.set_title('Pre-Basel III IRF: CSI → FDI') else: ax.text(0.5, 0.5, 'No valid pre-2007 IRF betas', transform=ax.transAxes, ha='center', va='center') else: ax.text(0.5, 0.5, 'No pre-2007 IRF data', transform=ax.transAxes, ha='center', va='center') fig.tight_layout() paths.extend(save_figure(fig, 'horse_race_oos2_basel_iii', subfolder='hierarchy')) # OOS3 figure: COVID recovery forecast + d comparison across recessions oos3 = _cache.get('oos3_data') if oos3 is not None: all_ratios = oos3.get('all_recession_ratios', []) fig, axes = setup_subplots(nrows=1, ncols=2, width=14, height=5) # Panel A: Forecast vs actual ax = axes[0] ax.bar(['Forecast\n(median pre-2020)', 'Actual\n(COVID)'], [oos3['median_p2r'], oos3['actual_p2r']], color=['steelblue', 'coral'], edgecolor='black', width=0.5) ax.set_ylabel('Peak-to-Recovery (months)') err = oos3['forecast_error'] ax.set_title(f'OOS3: IPG334S Recovery Forecast\nError: {err:+.1f}%') # Annotate training data for i, p2r in enumerate(oos3['p2r_train']): ax.scatter(0, p2r, c='gray', s=30, alpha=0.5, zorder=5, marker='_', linewidths=2) # Panel B: d1/d2 ratio across recessions ax = axes[1] if all_ratios: labels = [r['peak'][:4] for r in all_ratios] ratios = [r['ratio'] for r in all_ratios] colors = ['coral' if r['peak'] == '2020-02' else 'steelblue' for r in all_ratios] bars = ax.bar(labels, ratios, color=colors, edgecolor='black', width=0.6) ax.axhline(1.0, color='red', linestyle='--', alpha=0.6, label='d1/d2 = 1') ax.set_ylabel('d=1 / d=2 Recovery Ratio') ax.set_xlabel('Recession Peak Year') n_correct = sum(1 for r in ratios if r > 1.0) ax.set_title(f'd=2 Recovers Faster: {n_correct}/{len(ratios)} Recessions') ax.legend(fontsize=8) fig.tight_layout() paths.extend(save_figure(fig, 'horse_race_oos3_covid', subfolder='hierarchy')) return paths if __name__ == '__main__': res = main() for r in res: print(f" {r['name']:50s} {r['verdict']:15s} stat={r['statistic']}")

statistic	14
d1 d2 ratio	0.8061
forecast p2r	23
forecast error pct	-39.1304

Horse Race OOS Covid Recovery

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps