#!/usr/bin/env python3 """ Crisis Sequence Ordering Test Tests Proposition 3.4: K^2 < K for K in (0,1) implies correlation robustness (proportional to K^2) degrades FASTER than superadditivity (proportional to K) as K_eff -> 0. Observable prediction: cross-sector correlations INCREASE before aggregate OUTPUT DECLINES in pre-crisis periods. Three subtests: 1. crisis_sequence_ppi_correlations: PPI commodity input correlations vs INDPRO (all 6 NBER recessions, Newey-West corrected) 2. crisis_sequence_demand_only: Same but excluding supply-shock recessions (1973-74, 1980) to address commodity confounding 3. crisis_sequence_financial: Financial market correlations (BAA10Y, T10Y2Y, DGS1, DCOILWTICO) vs INDPRO — more direct proxy for K^2 channel Data: - INPUTS: 12 FRED PPI commodity groups (monthly, 1970-2026) - FINANCIAL: BAA10Y, T10Y2Y, DGS1, DCOILWTICO (monthly, 1986-2026) - OUTPUT: INDPRO (Industrial Production Index, monthly) Review Notes (2026-03-11): Theory match: PARTIAL. The Lean theorem (crisis_sequence_ordering) proves a SENSITIVITY ordering ((1-x)^2 < (1-x)), not a TEMPORAL ordering. The temporal implication is plausible (the faster-degrading channel hits threshold first) but not directly proven. Paper 3 Testable Prediction #3 (line 1706) makes the temporal prediction explicit. I/O audit: CORRECT. PPI → input correlation proxy; financial series → correlation robustness proxy; INDPRO → output decline proxy. Data-signal alignment: ADEQUATE (upgraded from WEAK after adding demand-only and financial subtests). - PPI correlations: REASONABLE but confounded by supply shocks - Financial correlations: REASONABLE — more direct proxy for K^2 channel - Demand-only subsample: addresses supply-shock confounding - Newey-West HAC: addresses serial correlation (replaces Bartlett) Code quality: CLEAN (bugs fixed from prior evaluation). Confidence in verdict: MEDIUM. Date of review: 2026-03-11 """ import warnings warnings.filterwarnings('ignore') import math import numpy as np import pandas as pd from scipy import stats from itertools import combinations from statsmodels.regression.linear_model import OLS from statsmodels.tools import add_constant import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.data import load_fred from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure # PPI commodity groups (production INPUTS) PPI_SECTORS = { 'WPU01': 'Farm Products', 'WPU02': 'Processed Foods', 'WPU03': 'Textiles & Apparel', 'WPU06': 'Chemicals', 'WPU07': 'Rubber & Plastics', 'WPU08': 'Lumber & Wood', 'WPU09': 'Pulp & Paper', 'WPU10': 'Metals', 'WPU11': 'Machinery', 'WPU12': 'Furniture', 'WPU13': 'Nonmetallic Minerals', 'WPU14': 'Transport Equipment', } # Financial series for correlation robustness proxy FINANCIAL_SERIES = { 'BAA10Y': 'Corporate Credit Spread', 'T10Y2Y': 'Yield Curve Slope', 'DGS1': '1-Year Treasury Yield', 'DCOILWTICO': 'Oil Price', } # NBER recessions (peak month) — COVID excluded (exogenous shock) RECESSIONS_ALL = [ ('1973-74 Oil+Overheating', '1973-11', 'supply'), ('1980 Volcker', '1980-01', 'supply'), ('1981-82 Volcker II', '1981-07', 'demand'), ('1990-91 S&L/Kuwait', '1990-07', 'demand'), ('2001 Dot-com', '2001-03', 'demand'), ('2007-09 GFC', '2007-12', 'demand'), ] RECESSIONS_DEMAND = [(n, p, t) for n, p, t in RECESSIONS_ALL if t == 'demand'] def _compute_rolling_correlation(growth, window=12): """Rolling mean pairwise correlation across sectors.""" pairs = list(combinations(growth.columns, 2)) corr_vals = [] for i in range(window - 1, len(growth)): block = growth.iloc[i - window + 1:i + 1] cm = block.corr() vals = [cm.loc[a, b] for a, b in pairs] corr_vals.append(np.nanmean(vals)) return pd.Series(corr_vals, index=growth.index[window - 1:]) def _newey_west_p(x, y, max_lag=12): """OLS regression y ~ x with Newey-West HAC standard errors. Returns (beta, t_stat, p_value, r_squared, n). """ mask = np.isfinite(x) & np.isfinite(y) x, y = x[mask], y[mask] n = len(x) if n < 20: return np.nan, np.nan, np.nan, np.nan, n X = add_constant(x) model = OLS(y, X).fit(cov_type='HAC', cov_kwds={'maxlags': max_lag}) beta = model.params[1] t_stat = model.tvalues[1] p_val = model.pvalues[1] r2 = model.rsquared return beta, t_stat, p_val, r2, n def _cross_corr_nw(c_series, y_series, lags, nw_maxlag=12): """Cross-correlation at multiple lags using Newey-West for p-values.""" common = c_series.index.intersection(y_series.index) c_full = c_series.loc[common].values y_neg = (-y_series.loc[common]).values results = {} best_r, best_lag, best_p = -1, 0, 1.0 for lag in lags: if lag == 0: cv, yv = c_full, y_neg else: cv, yv = c_full[:-lag], y_neg[lag:] r, p_naive = stats.pearsonr(cv, yv) beta, t_nw, p_nw, r2, n = _newey_west_p(cv, yv, max_lag=nw_maxlag) results[lag] = {'r': r, 'p_naive': p_naive, 'p_nw': p_nw, 't_nw': t_nw, 'n': n} if r > best_r: best_r, best_lag, best_p = r, lag, p_nw return results, best_r, best_lag, best_p def _threshold_crossing(C_input, Y_output, recessions): """Test threshold crossing ordering for a list of recessions.""" n_correct, n_total = 0, 0 details = [] for name, peak_str, _ in recessions: peak = pd.Timestamp(peak_str) start = peak - pd.DateOffset(months=24) baseline_start = peak - pd.DateOffset(months=84) C_base = C_input[(C_input.index >= baseline_start) & (C_input.index < start)] if len(C_base) < 24: continue thresh = C_base.mean() + 0.5 * C_base.std() Cw = C_input[(C_input.index >= start) & (C_input.index <= peak)] Yw = Y_output[(Y_output.index >= start) & (Y_output.index <= peak)] c_breaks = Cw[Cw > thresh] c_break = c_breaks.index[0] if len(c_breaks) > 0 else None y_breaks = Yw[Yw < 0] y_break = y_breaks.index[0] if len(y_breaks) > 0 else None if c_break is not None and y_break is not None: n_total += 1 lead_months = (y_break - c_break).days / 30.4 if lead_months > 0: n_correct += 1 details.append({'recession': name, 'lead_months': round(lead_months, 1)}) elif c_break is not None: n_total += 1 n_correct += 1 details.append({'recession': name, 'lead_months': 'corr only'}) elif y_break is not None: n_total += 1 details.append({'recession': name, 'lead_months': 'output only'}) p_binom = stats.binom.sf(n_correct - 1, n_total, 0.5) if n_total > 0 else np.nan return n_correct, n_total, p_binom, details def _pre_recession_trends(C_input, Y_output, recessions): """Test pre-recession rising correlation + falling output trends.""" n_both, n_trend = 0, 0 details = [] for name, peak_str, _ in recessions: peak = pd.Timestamp(peak_str) start = peak - pd.DateOffset(months=18) Cw = C_input[(C_input.index >= start) & (C_input.index <= peak)] Yw = Y_output[(Y_output.index >= start) & (Y_output.index <= peak)] if len(Cw) < 6: continue n_trend += 1 rc, _ = stats.spearmanr(np.arange(len(Cw)), Cw.values) ry, _ = stats.spearmanr(np.arange(len(Yw)), Yw.values) both_ok = rc > 0 and ry < 0 if both_ok: n_both += 1 details.append({'recession': name, 'corr_r': round(rc, 3), 'output_r': round(ry, 3), 'match': both_ok}) return n_both, n_trend, details def _combine_verdict(best_r, best_lag, p_nw, n_correct, n_total, n_both, n_trend): """Combine three sub-assessments into overall verdict.""" assessments = [] # Cross-correlation (Newey-West) sign_ok = best_r > 0 and best_lag > 0 if sign_ok and p_nw < 0.10: assessments.append('CONSISTENT') elif sign_ok: assessments.append('AMBIGUOUS') else: assessments.append('INCONSISTENT') # Threshold crossing if n_total > 0 and n_correct / n_total >= 0.6: assessments.append('CONSISTENT') elif n_total > 0 and n_correct / n_total >= 0.5: assessments.append('AMBIGUOUS') else: assessments.append('INCONSISTENT') # Pre-recession trends if n_trend > 0 and n_both / n_trend >= 0.5: assessments.append('CONSISTENT') elif n_trend > 0 and n_both / n_trend >= 1.0 / 3.0: assessments.append('AMBIGUOUS') else: assessments.append('INCONSISTENT') n_c = sum(1 for a in assessments if a == 'CONSISTENT') n_i = sum(1 for a in assessments if a == 'INCONSISTENT') if n_c >= 2: overall = 'CONSISTENT' elif n_i >= 2: overall = 'INCONSISTENT' else: overall = 'AMBIGUOUS' return overall, assessments _cache = {} def main(): results = [] # ── Load PPI data ── ppi_dfs = {} for code in PPI_SECTORS: try: ppi_dfs[code] = load_fred(code) except FileNotFoundError: pass if len(ppi_dfs) < 6: results.append(log_result("crisis_sequence_ppi_correlations", np.nan, np.nan, "AMBIGUOUS", {"reason": f"Only {len(ppi_dfs)} PPI series"}, description="Insufficient PPI series for crisis sequence ordering test.")) return results ppi_panel = pd.DataFrame(ppi_dfs).dropna() ppi_panel = ppi_panel[ppi_panel.index >= '1970-01-01'] # ── Load output ── try: indpro = load_fred('INDPRO') except FileNotFoundError: results.append(log_result("crisis_sequence_ppi_correlations", np.nan, np.nan, "AMBIGUOUS", {"reason": "INDPRO not available"}, description="INDPRO not available for crisis sequence ordering test.")) return results # ── Compute PPI input correlations and output growth ── ppi_inflation = np.log(ppi_panel).diff(12).dropna() C_ppi = _compute_rolling_correlation(ppi_inflation, window=12) C_ppi.name = 'ppi_input_corr' Y_output = (np.log(indpro) - np.log(indpro.shift(12))).dropna() Y_output.name = 'output_growth' lag_range = list(range(0, 19, 3)) # ════════════════════════════════════════════════════════════ # Subtest 1: PPI correlations — all recessions, Newey-West # ════════════════════════════════════════════════════════════ xcorr_ppi, best_r_ppi, best_lag_ppi, best_p_ppi = _cross_corr_nw( C_ppi, Y_output, lag_range, nw_maxlag=12) nc_ppi, nt_ppi, pb_ppi, thresh_ppi = _threshold_crossing( C_ppi, Y_output, RECESSIONS_ALL) nb_ppi, ntr_ppi, trend_ppi = _pre_recession_trends( C_ppi, Y_output, RECESSIONS_ALL) verd_ppi, subs_ppi = _combine_verdict( best_r_ppi, best_lag_ppi, best_p_ppi, nc_ppi, nt_ppi, nb_ppi, ntr_ppi) results.append(log_result( "crisis_sequence_ppi_correlations", stat=best_r_ppi, p=best_p_ppi, verd=verd_ppi, details={ "best_lag_months": best_lag_ppi, "best_r": round(best_r_ppi, 4), "p_newey_west": round(best_p_ppi, 4) if np.isfinite(best_p_ppi) else None, "p_naive": float(f"{xcorr_ppi[best_lag_ppi]['p_naive']:.2e}") if best_lag_ppi in xcorr_ppi else None, "threshold_correct": f"{nc_ppi}/{nt_ppi}", "p_binom": round(pb_ppi, 4) if np.isfinite(pb_ppi) else None, "trend_correct": f"{nb_ppi}/{ntr_ppi}", "sub_assessments": subs_ppi, "n_ppi_sectors": len(ppi_dfs), "xcorr_profile": {str(k): round(v['r'], 4) for k, v in xcorr_ppi.items()}, }, description=( "Tests K^2 < K degradation ordering: PPI commodity input correlations " "should rise before INDPRO output declines. Method: (1) Pearson cross-" "correlation at lags 0-18m with Newey-West HAC p-values (maxlag=12); " "(2) threshold crossing ordering across 6 NBER recessions (1973-2007); " "(3) pre-recession 18-month Spearman trend. Prediction: positive r at " "positive lag. CAUTION: Supply-shock recessions (1973, 1980) confound — " "commodity price co-movement mechanically precedes output decline." ), )) _cache['C_ppi'] = C_ppi _cache['Y_output'] = Y_output _cache['xcorr_ppi'] = xcorr_ppi _cache['best_lag_ppi'] = best_lag_ppi _cache['best_r_ppi'] = best_r_ppi _cache['best_p_ppi'] = best_p_ppi _cache['thresh_ppi'] = thresh_ppi # ════════════════════════════════════════════════════════════ # Subtest 2: PPI correlations — demand-only recessions # ════════════════════════════════════════════════════════════ nc_dem, nt_dem, pb_dem, thresh_dem = _threshold_crossing( C_ppi, Y_output, RECESSIONS_DEMAND) nb_dem, ntr_dem, trend_dem = _pre_recession_trends( C_ppi, Y_output, RECESSIONS_DEMAND) # Cross-correlation still uses full time series (same lag structure) # but verdict combines with demand-only threshold/trend verd_dem, subs_dem = _combine_verdict( best_r_ppi, best_lag_ppi, best_p_ppi, nc_dem, nt_dem, nb_dem, ntr_dem) results.append(log_result( "crisis_sequence_demand_only", stat=best_r_ppi, p=best_p_ppi, verd=verd_dem, details={ "best_lag_months": best_lag_ppi, "best_r": round(best_r_ppi, 4), "p_newey_west": round(best_p_ppi, 4) if np.isfinite(best_p_ppi) else None, "threshold_correct": f"{nc_dem}/{nt_dem}", "p_binom": round(pb_dem, 4) if np.isfinite(pb_dem) else None, "trend_correct": f"{nb_dem}/{ntr_dem}", "sub_assessments": subs_dem, "recessions_used": [n for n, _, _ in RECESSIONS_DEMAND], }, description=( "Robustness check excluding supply-shock recessions (1973-74 Oil, 1980 " "Volcker) where commodity price co-movement mechanically causes both rising " "PPI correlations and output decline. Tests only demand-driven recessions " "(1981-82, 1990-91, 2001, GFC) where the K^2 < K degradation mechanism " "applies without supply-shock confounding. Same methodology as subtest 1." ), )) _cache['thresh_dem'] = thresh_dem # ════════════════════════════════════════════════════════════ # Subtest 3: Financial correlations — direct K^2 proxy # ════════════════════════════════════════════════════════════ fin_dfs = {} for code in FINANCIAL_SERIES: try: s = load_fred(code) # Resample daily to monthly s_monthly = s.resample('MS').mean().dropna() fin_dfs[code] = s_monthly except FileNotFoundError: pass if len(fin_dfs) >= 3: fin_panel = pd.DataFrame(fin_dfs).dropna() fin_panel = fin_panel[fin_panel.index >= '1986-01-01'] # Monthly changes (first differences for spreads, log-returns for oil) fin_changes = pd.DataFrame(index=fin_panel.index[1:]) for col in fin_panel.columns: if col == 'DCOILWTICO': fin_changes[col] = np.log(fin_panel[col]).diff().iloc[1:] else: fin_changes[col] = fin_panel[col].diff().iloc[1:] fin_changes = fin_changes.dropna() C_fin = _compute_rolling_correlation(fin_changes, window=12) C_fin.name = 'financial_corr' # Financial-era recessions only (post-1986) RECESSIONS_FIN = [(n, p, t) for n, p, t in RECESSIONS_ALL if pd.Timestamp(p).year >= 1986] xcorr_fin, best_r_fin, best_lag_fin, best_p_fin = _cross_corr_nw( C_fin, Y_output, lag_range, nw_maxlag=12) nc_fin, nt_fin, pb_fin, thresh_fin = _threshold_crossing( C_fin, Y_output, RECESSIONS_FIN) nb_fin, ntr_fin, trend_fin = _pre_recession_trends( C_fin, Y_output, RECESSIONS_FIN) verd_fin, subs_fin = _combine_verdict( best_r_fin, best_lag_fin, best_p_fin, nc_fin, nt_fin, nb_fin, ntr_fin) results.append(log_result( "crisis_sequence_financial", stat=best_r_fin, p=best_p_fin, verd=verd_fin, details={ "best_lag_months": best_lag_fin, "best_r": round(best_r_fin, 4), "p_newey_west": round(best_p_fin, 4) if np.isfinite(best_p_fin) else None, "threshold_correct": f"{nc_fin}/{nt_fin}", "p_binom": round(pb_fin, 4) if np.isfinite(pb_fin) else None, "trend_correct": f"{nb_fin}/{ntr_fin}", "sub_assessments": subs_fin, "n_financial_series": len(fin_dfs), "financial_series": list(fin_dfs.keys()), "recessions_used": [n for n, _, _ in RECESSIONS_FIN], "xcorr_profile": {str(k): round(v['r'], 4) for k, v in xcorr_fin.items()}, }, description=( "Tests K^2 < K degradation using financial market correlations as a " "more direct proxy for the correlation robustness channel. Computes " "rolling 12-month pairwise correlation among BAA10Y (credit spread), " "T10Y2Y (yield curve), DGS1 (short-rate), and DCOILWTICO (oil) monthly " "changes. Tests whether rising financial correlations lead INDPRO output " "declines. Newey-West HAC p-values. Post-1986 sample (4 recessions). " "Financial correlations are a more direct measure of diversification " "failure (K^2 channel) than PPI commodity correlations." ), )) _cache['C_fin'] = C_fin _cache['xcorr_fin'] = xcorr_fin _cache['best_lag_fin'] = best_lag_fin _cache['best_r_fin'] = best_r_fin _cache['best_p_fin'] = best_p_fin _cache['thresh_fin'] = thresh_fin else: results.append(log_result( "crisis_sequence_financial", stat=np.nan, p=np.nan, verd="AMBIGUOUS", details={"reason": f"Only {len(fin_dfs)} financial series available"}, description="Tests K^2 < K ordering with financial correlations. Insufficient financial series available." )) _cache['results'] = results return results _ALL_SUBTESTS = [ 'crisis_sequence_ppi_correlations', 'crisis_sequence_demand_only', 'crisis_sequence_financial', ] def figures(): """Generate figures for each crisis sequence subtest.""" if not _cache: main() paths = [] C_ppi = _cache.get('C_ppi') Y_output = _cache.get('Y_output') xcorr_ppi = _cache.get('xcorr_ppi', {}) best_lag_ppi = _cache.get('best_lag_ppi', 0) best_r_ppi = _cache.get('best_r_ppi', 0) best_p_ppi = _cache.get('best_p_ppi', 1.0) if C_ppi is None or Y_output is None: return paths # ── Figure 1: PPI correlations (covers subtests 1 & 2) ── fig, axes = setup_subplots(1, 2, width=14, height=5) # Panel A: Time series ax = axes[0] common = C_ppi.index.intersection(Y_output.index) ax.plot(C_ppi.loc[common].index, C_ppi.loc[common].values, 'b-', linewidth=1, alpha=0.7, label='Mean PPI pairwise corr') ax2 = ax.twinx() ax2.plot(Y_output.loc[common].index, Y_output.loc[common].values * 100, 'r-', linewidth=1, alpha=0.7, label='IP 12m growth (%)') ax2.set_ylabel('IP 12m growth (%)', color='r') ax2.tick_params(axis='y', labelcolor='r') for name, peak_str, rtype in RECESSIONS_ALL: peak = pd.Timestamp(peak_str) color = '#ff9999' if rtype == 'supply' else '#cccccc' ax.axvspan(peak, peak + pd.DateOffset(months=12), alpha=0.2, color=color) ax.set_ylabel('Mean pairwise PPI correlation', color='b') ax.tick_params(axis='y', labelcolor='b') ax.set_title('PPI Input Correlations vs Output Growth\n' '(pink = supply-shock recessions)') ax.legend(loc='upper left', fontsize=8) ax2.legend(loc='upper right', fontsize=8) # Panel B: Cross-correlation function ax = axes[1] lags = sorted(xcorr_ppi.keys()) rs = [xcorr_ppi[l]['r'] for l in lags] colors = ['#e74c3c' if l == best_lag_ppi else '#2c3e50' for l in lags] ax.bar(lags, rs, width=2.5, color=colors, edgecolor='black', linewidth=0.5) ax.axhline(0, color='gray', linewidth=0.5) ax.set_xlabel('Lag (months, C_input leads)') ax.set_ylabel('Pearson r (C_input vs -Y_output)') p_str = f'{best_p_ppi:.4f}' if np.isfinite(best_p_ppi) else 'N/A' ax.set_title(f'Cross-Correlation: r={best_r_ppi:.3f} at lag={best_lag_ppi}m\n' f'(Newey-West p={p_str})') fig.tight_layout() paths.extend(save_figure(fig, 'crisis_sequence_ppi_correlations', subfolder='dynamics')) # ── Figure 2: Financial correlations ── C_fin = _cache.get('C_fin') xcorr_fin = _cache.get('xcorr_fin', {}) best_lag_fin = _cache.get('best_lag_fin', 0) best_r_fin = _cache.get('best_r_fin', 0) best_p_fin = _cache.get('best_p_fin', 1.0) if C_fin is not None: fig, axes = setup_subplots(1, 2, width=14, height=5) # Panel A: Time series ax = axes[0] common_fin = C_fin.index.intersection(Y_output.index) ax.plot(C_fin.loc[common_fin].index, C_fin.loc[common_fin].values, 'b-', linewidth=1, alpha=0.7, label='Mean financial pairwise corr') ax2 = ax.twinx() ax2.plot(Y_output.loc[common_fin].index, Y_output.loc[common_fin].values * 100, 'r-', linewidth=1, alpha=0.7, label='IP 12m growth (%)') ax2.set_ylabel('IP 12m growth (%)', color='r') ax2.tick_params(axis='y', labelcolor='r') for name, peak_str, rtype in RECESSIONS_ALL: peak = pd.Timestamp(peak_str) if peak.year >= 1986: ax.axvspan(peak, peak + pd.DateOffset(months=12), alpha=0.15, color='gray') ax.set_ylabel('Mean pairwise financial correlation', color='b') ax.tick_params(axis='y', labelcolor='b') ax.set_title('Financial Correlations vs Output Growth\n' '(BAA10Y, T10Y2Y, DGS1, DCOILWTICO)') ax.legend(loc='upper left', fontsize=8) ax2.legend(loc='upper right', fontsize=8) # Panel B: Cross-correlation function ax = axes[1] lags_fin = sorted(xcorr_fin.keys()) rs_fin = [xcorr_fin[l]['r'] for l in lags_fin] colors_fin = ['#e74c3c' if l == best_lag_fin else '#2c3e50' for l in lags_fin] ax.bar(lags_fin, rs_fin, width=2.5, color=colors_fin, edgecolor='black', linewidth=0.5) ax.axhline(0, color='gray', linewidth=0.5) ax.set_xlabel('Lag (months, C_financial leads)') ax.set_ylabel('Pearson r (C_financial vs -Y_output)') p_str_fin = f'{best_p_fin:.4f}' if np.isfinite(best_p_fin) else 'N/A' ax.set_title(f'Cross-Correlation: r={best_r_fin:.3f} at lag={best_lag_fin}m\n' f'(Newey-West p={p_str_fin})') fig.tight_layout() paths.extend(save_figure(fig, 'crisis_sequence_financial', subfolder='dynamics')) return paths if __name__ == '__main__': res = main() for r in res: p_str = f'{r["p_value"]:.4f}' if np.isfinite(r['p_value']) else 'N/A' print(f"{r['name']}: {r['verdict']} (stat={r['statistic']:.4f}, p={p_str})")

p value	1.373e-6
statistic	0.1908
n ppi sectors	12
best lag months	6
threshold correct	3/4
trend both correct	2/6

Crisis Sequence Ordering

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps