#!/usr/bin/env python3 """ Firm-Level Predictions Test — CES Free Energy Framework ======================================================== Tests firm-level predictions from Paper 3 (Dynamics on CES Landscape): Test #4: FDT (Fluctuation-Dissipation) — sigma^2 = alpha + T*chi cross-section Test #6: Equicorrelation — compound symmetry in correlation matrix Test #7: Directed Coupling — IRF asymmetry correlates with supply chain Test #8: Cyclical rho — countercyclical CES parameter (correlation proxy) Test #8b: Structural rho — countercyclical CES parameter (factor share identity) Data: 17 FRED Industrial Production series, BEA supply chain data (embedded) Migrated from scripts/test_firm_predictions.py Review Notes (2026-03-12) ------------------------- Theory match: CORRECT — Paper 3, Theorem 3.1 (VRI) and Theorem Minsky trap σ²_n = T · χ_n proved in Lean (gibbs_vri_full, static_vri, T_measurement_vri) Cyclical rho: minsky_cross_coupling in Lean (∂f_T/∂ρ < 0) I/O audit: CORRECT — all variables correctly oriented Data-signal alignment: WEAK for VRI cross-section (N=17), WEAK for cyclical rho (correlation proxy confounds ρ with common-shock variation), STRONG for other 2 Code quality: CLEAN after fixes BUG FIX (2026-03-12): cyclical_rho used step=6 (90% overlap) creating severely autocorrelated rolling windows. Changed to step=window//2 (50% overlap). Result: CONSISTENT → AMBIGUOUS. Matches paper's own result (tau=-0.069, p=0.09). NEW (2026-03-12): structural_rho subtest added — factor share identity method. Also AMBIGUOUS (tau=-0.018, p=0.329). Confirms proxy wasn't the bottleneck. Confidence: MEDIUM — 2/5 CONSISTENT (equicorrelation, directed coupling); 3/5 AMBIGUOUS (VRI cross-section, cyclical rho, structural rho) Known limitations: - Cross-sectional VRI has N=17 (power failure, not theory failure) - Bivariate VAR chi estimates are noisy (ignores cross-sector transmission) - Aggregate countercyclicality of rho is genuinely weak (paper acknowledges) - Cross-industry and cross-country variation stronger (see U3 Minsky test) """ import warnings warnings.filterwarnings('ignore') import os import numpy as np import pandas as pd from scipy import stats from common.config import COMPILED_DIR from common.data import load_fred from common.stats import verdict, log_result from common.plotting import setup_figure, setup_subplots, save_figure try: import statsmodels.api as sm from statsmodels.tsa.api import VAR HAS_SM = True except ImportError: HAS_SM = False START_DATE = '1972-01-01' LEAF_SECTORS = [ 'IPG331S', 'IPG332S', 'IPG333S', 'IPG334S', 'IPG335S', 'IPG336S', 'IPG337S', 'IPG339S', 'IPG311A2S', 'IPG313A4S', 'IPG315A6S', 'IPG321S', 'IPG322S', 'IPG323S', 'IPG324S', 'IPG325S', 'IPG326S', ] DURABLE_SECTORS = [ 'IPG331S', 'IPG332S', 'IPG333S', 'IPG334S', 'IPG335S', 'IPG336S', 'IPG337S', 'IPG339S', ] LABELS = { 'INDPRO': 'Total IP', 'IPMAN': 'Manufacturing', 'IPDMAN': 'Durable Goods', 'IPNMAN': 'Nondurable Goods', 'IPG331S': 'Primary Metals', 'IPG332S': 'Fabricated Metals', 'IPG333S': 'Machinery', 'IPG334S': 'Computer/Elec', 'IPG335S': 'Electrical Equip', 'IPG336S': 'Transport Equip', 'IPG337S': 'Furniture', 'IPG339S': 'Misc Mfg', 'IPG311A2S': 'Food/Bev/Tobacco', 'IPG313A4S': 'Textiles', 'IPG315A6S': 'Apparel/Leather', 'IPG321S': 'Wood Products', 'IPG322S': 'Paper', 'IPG323S': 'Printing', 'IPG324S': 'Petroleum/Coal', 'IPG325S': 'Chemicals', 'IPG326S': 'Plastics/Rubber', } # BEA supply chain direction matrix for 8 durable sectors BEA_DIRECTION = np.array([ [ 0, +1, +1, +1, +1, 0, +1, +1], [ -1, 0, +1, +1, +1, 0, +1, +1], [ -1, -1, 0, -1, -1, 0, -1, -1], [ -1, -1, +1, 0, +1, 0, +1, +1], [ -1, -1, +1, -1, 0, 0, +1, +1], [ 0, 0, 0, 0, 0, 0, 0, 0], [ -1, -1, +1, -1, -1, 0, 0, +1], [ -1, -1, +1, -1, -1, 0, -1, 0], ], dtype=float) _ALL_SUBTESTS = [ 'firm_predictions_variance_response', 'firm_predictions_equicorrelation', 'firm_predictions_directed_coupling', 'firm_predictions_cyclical_rho', 'firm_predictions_structural_rho', ] def _load_growth_data(): """Load FRED IP series and compute growth rates.""" all_ids = set(LEAF_SECTORS + DURABLE_SECTORS + ['INDPRO', 'IPDMAN', 'IPMAN']) raw = {} growth = {} for sid in sorted(all_ids): try: s = load_fred(sid) s = s[s.index >= START_DATE] raw[sid] = s growth[sid] = np.log(s).diff().dropna() except (FileNotFoundError, Exception): continue leaf_df = pd.DataFrame({sid: growth[sid] for sid in LEAF_SECTORS if sid in growth}) leaf_df = leaf_df.dropna() dur_df = pd.DataFrame({sid: growth[sid] for sid in DURABLE_SECTORS if sid in growth}) dur_df = dur_df.dropna() return raw, growth, leaf_df, dur_df def test_fdt(growth, leaf_df): """Test #4: FDT — sigma^2 = alpha + T*chi across sectors (cross-section).""" desc = ("Tests the fluctuation-dissipation theorem: cross-sector variance should " "scale linearly with response to aggregate shocks (sigma^2 = alpha + T*chi). " "OLS of sector log-growth variance on VAR-estimated cumulative impulse response " "chi (12-month response of sector IP to 1-unit INDPRO shock) across 17 FRED IP " "leaf sectors. Prediction: positive slope T > 0. Known low-power test (N=17); " "Paper 3 reports R^2=0.009 for this specification.") if 'INDPRO' not in growth or not HAS_SM: return log_result('firm_predictions_variance_response', np.nan, np.nan, 'SKIPPED', description=desc) indpro_g = growth['INDPRO'] sectors = [s for s in LEAF_SECTORS if s in growth] sector_stats = [] for sid in sectors: combined = pd.DataFrame({'sector': growth[sid], 'indpro': indpro_g}).dropna() if len(combined) < 150: continue sigma2 = combined['sector'].var() try: var_data = combined[['sector', 'indpro']].values model = VAR(var_data) p = max(model.select_order(maxlags=8).bic, 1) var_fit = model.fit(p) irf = var_fit.irf(12) chi = np.sum(irf.irfs[1:13, 0, 1]) sector_stats.append({'sector': sid, 'sigma2': sigma2, 'chi': chi, 'label': LABELS.get(sid, sid)}) except Exception: continue if len(sector_stats) < 5: return log_result('firm_predictions_variance_response', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(sector_stats)} sectors'}, description=desc) sdf = pd.DataFrame(sector_stats) X = sm.add_constant(sdf['chi'].values) y = sdf['sigma2'].values ols = sm.OLS(y, X).fit() T_est = ols.params[1] p_slope = ols.pvalues[1] r2 = ols.rsquared r_sp, p_sp = stats.spearmanr(sdf['sigma2'], sdf['chi']) sign_ok = T_est > 0 if sign_ok and (p_slope < 0.10 or p_sp < 0.10): verd = 'CONSISTENT' elif sign_ok: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['fdt_data'] = sdf _cache['fdt_ols'] = {'T': T_est, 'intercept': ols.params[0], 'r2': r2} return log_result('firm_predictions_variance_response', T_est, p_slope, verd, details={'N': len(sdf), 'R2': r2, 'spearman_r': r_sp, 'spearman_p': p_sp}, description=desc) def test_equicorrelation(leaf_df): """Test #6: CES compound symmetry — equal off-diagonal correlations.""" if leaf_df.shape[1] < 5: return log_result('firm_predictions_equicorrelation', np.nan, np.nan, 'SKIPPED', description="Tests CES compound symmetry: off-diagonal correlations in the IP sector correlation matrix should be approximately equal. Computes coefficient of variation of off-diagonal correlations and Frobenius distance from equicorrelation matrix, benchmarked against a permutation null (500 shuffles).") n = leaf_df.shape[1] corr_matrix = leaf_df.corr().values off_diag = [] for i in range(n): for j in range(i + 1, n): off_diag.append(corr_matrix[i, j]) off_diag = np.array(off_diag) mean_corr = np.mean(off_diag) std_corr = np.std(off_diag) cv_corr = std_corr / abs(mean_corr) if abs(mean_corr) > 1e-10 else np.inf equicorr = np.full_like(corr_matrix, mean_corr) np.fill_diagonal(equicorr, 1.0) frob_ratio = (np.linalg.norm(corr_matrix - equicorr, 'fro') / np.linalg.norm(corr_matrix, 'fro')) # Permutation null with fixed seed for reproducibility rng = np.random.RandomState(42) null_cvs = [] data_vals = leaf_df.values.copy() for _ in range(500): shuffled = data_vals.copy() for col in range(n): rng.shuffle(shuffled[:, col]) perm_corr = np.corrcoef(shuffled, rowvar=False) perm_od = [] for i in range(n): for j in range(i + 1, n): perm_od.append(perm_corr[i, j]) perm_od = np.array(perm_od) pm = np.mean(perm_od) ps = np.std(perm_od) null_cvs.append(ps / abs(pm) if abs(pm) > 1e-10 else np.inf) null_cvs = np.array(null_cvs) null_median = np.median(null_cvs) if cv_corr < null_median and frob_ratio < 0.5: verd = 'CONSISTENT' elif cv_corr < null_median: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['equicorr_data'] = { 'corr_matrix': corr_matrix, 'labels': [LABELS.get(s, s) for s in leaf_df.columns], 'cv': cv_corr, 'mean_corr': mean_corr, 'frob_ratio': frob_ratio, 'null_median': null_median } return log_result('firm_predictions_equicorrelation', cv_corr, np.nan, verd, details={'mean_corr': mean_corr, 'cv': cv_corr, 'frob_ratio': frob_ratio, 'null_cv_median': null_median}, description="Tests CES compound symmetry: off-diagonal correlations in the 17-sector IP correlation matrix should be approximately equal. Computes coefficient of variation of off-diagonal correlations and Frobenius distance from the equicorrelation matrix, benchmarked against a 500-shuffle permutation null (seed=42). Prediction: CV below null median and Frobenius ratio below 0.5.") def test_directed_coupling(dur_df): """Test #7: IRF asymmetry correlates with BEA supply chain direction.""" desc = ("Tests whether VAR impulse response asymmetry aligns with BEA supply chain " "direction. Decomposes the cumulative MA(12) coefficient matrix into symmetric " "and antisymmetric parts, then uses Spearman rank correlation between " "antisymmetric off-diagonal elements and the BEA direction matrix for 8 durable " "goods sectors. Prediction: positive Spearman r, asymmetry ratio > 0.3, and " "sign agreement > 50%.") if not HAS_SM or dur_df.shape[1] < 3: return log_result('firm_predictions_directed_coupling', np.nan, np.nan, 'SKIPPED', description=desc) n_vars = dur_df.shape[1] try: model = VAR(dur_df.values) p = max(model.select_order(maxlags=6).bic, 1) var_fit = model.fit(p) ma_coefs = var_fit.ma_rep(12) L = np.sum(ma_coefs[1:13], axis=0) S = (L + L.T) / 2 A = (L.T - L) / 2 frob_L = np.linalg.norm(L, 'fro') frob_A = np.linalg.norm(A, 'fro') asym_ratio = frob_A / frob_L if frob_L > 0 else np.nan D = BEA_DIRECTION[:n_vars, :n_vars] od_a, od_d = [], [] for i in range(n_vars): for j in range(n_vars): if i != j: od_a.append(A[i, j]) od_d.append(D[i, j]) od_a = np.array(od_a) od_d = np.array(od_d) r_sp, p_sp = stats.spearmanr(od_a, od_d) directed_mask = od_d != 0 if directed_mask.sum() > 0: sign_agree = np.sum(np.sign(od_a[directed_mask]) == np.sign(od_d[directed_mask])) sign_frac = sign_agree / directed_mask.sum() else: sign_frac = 0.0 if asym_ratio > 0.3 and r_sp > 0 and sign_frac > 0.5: verd = 'CONSISTENT' elif asym_ratio > 0.3 and r_sp > 0: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['directed_data'] = { 'A': A, 'D': D, 'L': L, 'od_a': od_a, 'od_d': od_d, 'labels': [LABELS.get(s, s) for s in dur_df.columns], 'asym_ratio': asym_ratio, 'sign_frac': sign_frac, 'r_sp': r_sp, 'p_sp': p_sp } return log_result('firm_predictions_directed_coupling', r_sp, p_sp, verd, details={'asym_ratio': asym_ratio, 'sign_frac': sign_frac, 'N_pairs': len(od_a)}, description=desc) except Exception: return log_result('firm_predictions_directed_coupling', np.nan, np.nan, 'SKIPPED', description=desc) def test_cyclical_rho(growth, dur_df): """Test #8: rho is countercyclical — rises in contractions.""" if 'INDPRO' not in growth or dur_df.shape[1] < 3: return log_result('firm_predictions_cyclical_rho', np.nan, np.nan, 'SKIPPED', description="Tests whether the CES parameter rho is countercyclical: average pairwise correlation among durable goods sectors (a rho proxy) should rise when INDPRO growth falls. Uses Kendall tau between rolling 60-month average pairwise correlation and rolling INDPRO growth.") n_dur = dur_df.shape[1] T = len(dur_df) window = 60 rolling_corr = [] step = window // 2 # 50% overlap (standard; avoids autocorrelation inflation) for start in range(0, T - window + 1, step): end = start + window sub = dur_df.iloc[start:end] center_date = sub.index[window // 2] corr_mat = sub.corr().values off_diag = [] for i in range(n_dur): for j in range(i + 1, n_dur): off_diag.append(corr_mat[i, j]) rolling_corr.append({'center_date': center_date, 'avg_corr': np.mean(off_diag)}) proxy_df = pd.DataFrame(rolling_corr).set_index('center_date') indpro_cycle = growth['INDPRO'].rolling(60, min_periods=30).mean() common = proxy_df.index.intersection(indpro_cycle.dropna().index) if len(common) < 10: return log_result('firm_predictions_cyclical_rho', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(common)} aligned windows'}, description="Tests whether the CES parameter rho is countercyclical: average pairwise correlation among durable goods sectors (a rho proxy) should rise when INDPRO growth falls. Uses Kendall tau between rolling 60-month average pairwise correlation and rolling INDPRO growth.") proxy_aligned = proxy_df.loc[common, 'avg_corr'] cycle_aligned = indpro_cycle.reindex(common) tau, p_tau = stats.kendalltau(proxy_aligned.values, cycle_aligned.values) # Prediction: negative tau (correlation rises when growth falls = countercyclical rho) sign_ok = tau < 0 verd = verdict(p_tau, sign_ok, threshold=0.10) _cache['cyclical_data'] = { 'proxy_df': proxy_df, 'indpro_cycle': indpro_cycle, 'tau': tau, 'p_tau': p_tau, 'common': common } return log_result('firm_predictions_cyclical_rho', tau, p_tau, verd, details={'N_windows': len(common)}, description=( "Tests countercyclical rho (Paper 3, Thm Minsky trap): average " "pairwise correlation among 8 durable goods IP sectors (rho proxy) " "should rise when INDPRO growth falls. Kendall tau of rolling " "60-month average pairwise correlation vs rolling INDPRO growth, " "50% overlap windows. Prediction: negative tau. " "AMBIGUOUS: sign correct but not significant — matches paper's own " "aggregate result (tau=-0.069, p=0.09). Cross-industry variation " "is stronger (see U3 Minsky rho-shift test)." )) def test_structural_rho(): """Test #8b: Structural rho from factor shares is countercyclical. Paper 3, Remark on Solow residual decomposition (line 390): Rolling 10-year OLS of Δlog(s_L/s_K) on Δlog(k) gives slope = -ρ̂(t). Prediction: Kendall τ(ρ̂, output gap) < 0. Paper reports: τ = -0.069, p = 0.09. """ desc = ( "Structural rho countercyclicality (Paper 3, Remark Solow decomposition): " "estimates rho(t) from factor share identity Δlog(s_L/s_K) = -ρ·Δlog(k) " "via rolling 40-quarter OLS. Tests Kendall tau of rho_hat vs HP-filtered " "output gap (GDPC1, lambda=1600). Data: PRS85006173 (labor share), GPDIC1 " "(investment → perpetual inventory capital), HOANBS (hours). " "Prediction: negative tau (rho rises in contractions = countercyclical). " "Paper reports tau=-0.069, p=0.09." ) if not HAS_SM: return log_result('firm_predictions_structural_rho', np.nan, np.nan, 'SKIPPED', description=desc + " Requires statsmodels.") try: labor_share = load_fred('PRS85006173') # nonfarm business labor share gdp = load_fred('GDPC1') # real GDP investment = load_fred('GPDIC1') # real gross private domestic investment hours = load_fred('HOANBS') # hours of all persons except FileNotFoundError as e: return log_result('firm_predictions_structural_rho', np.nan, np.nan, 'SKIPPED', details={'reason': str(e)}, description=desc) # Align quarterly data df = pd.DataFrame({ 'sL': labor_share, 'gdp': gdp, 'invest': investment, 'hours': hours }).dropna() if len(df) < 80: return log_result('firm_predictions_structural_rho', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(df)} quarterly obs'}, description=desc) # Construct capital stock via perpetual inventory # K(t+1) = (1-delta)*K(t) + I(t), delta = 10%/year = 2.5%/quarter delta_q = 0.025 invest_vals = df['invest'].values # Initial capital: K(0) = I(0) / (g + delta) where g = avg investment growth g_invest = np.mean(np.diff(np.log(invest_vals[:40]))) K = np.zeros(len(df)) K[0] = invest_vals[0] / max(g_invest + delta_q, 0.01) for t in range(1, len(df)): K[t] = (1 - delta_q) * K[t - 1] + invest_vals[t - 1] # Capital-labor ratio k = K / L k = K / df['hours'].values # Labor vs capital share # PRS85006173 is an INDEX (2017=100), not a percent. Calibrate to actual # share: nonfarm business labor share was ~57% in 2017. sL_raw = df['sL'].values idx_2017 = np.mean(sL_raw[(df.index >= '2017') & (df.index < '2018')]) if idx_2017 <= 0: idx_2017 = 100.0 sL = 0.57 * sL_raw / idx_2017 sK = 1.0 - sL # Check for valid values valid = (sL > 0.01) & (sK > 0.01) & (k > 0) if valid.sum() < 80: return log_result('firm_predictions_structural_rho', np.nan, np.nan, 'SKIPPED', details={'reason': f'Too few valid share observations ({valid.sum()})'}, description=desc) # Compute Δlog(s_L/s_K) and Δlog(k) log_share_ratio = np.log(sL / sK) log_k = np.log(k) d_log_share = np.diff(log_share_ratio) d_log_k = np.diff(log_k) # Rolling 40-quarter OLS: Δlog(s_L/s_K) = const - ρ·Δlog(k) window = 40 rho_hat = [] rho_dates = [] dates = df.index[1:] # dates for differenced series for start in range(len(d_log_k) - window + 1): end = start + window y = d_log_share[start:end] x = d_log_k[start:end] X = sm.add_constant(x) try: ols = sm.OLS(y, X).fit() rho_est = -ols.params[1] # slope = -ρ rho_hat.append(rho_est) rho_dates.append(dates[start + window // 2]) except Exception: continue if len(rho_hat) < 20: return log_result('firm_predictions_structural_rho', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(rho_hat)} rolling estimates'}, description=desc) rho_series = pd.Series(rho_hat, index=rho_dates) # HP-filter GDP for output gap (lambda=1600 for quarterly) from statsmodels.tsa.filters.hp_filter import hpfilter gdp_aligned = df['gdp'].reindex(rho_series.index, method='nearest') gdp_aligned = gdp_aligned.dropna() common_idx = rho_series.index.intersection(gdp_aligned.index) if len(common_idx) < 20: return log_result('firm_predictions_structural_rho', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(common_idx)} aligned obs'}, description=desc) log_gdp = np.log(gdp_aligned.loc[common_idx]) _, trend = hpfilter(log_gdp.values, lamb=1600) output_gap = log_gdp.values - trend rho_aligned = rho_series.loc[common_idx].values # Kendall tau: ρ̂ vs output gap — prediction: negative tau, p_tau = stats.kendalltau(rho_aligned, output_gap) sign_ok = tau < 0 verd = verdict(p_tau / 2 if sign_ok else 1.0, sign_ok, threshold=0.10) _cache['structural_rho'] = { 'rho_series': rho_series, 'output_gap': pd.Series(output_gap, index=common_idx), 'tau': tau, 'p_tau': p_tau, 'rho_mean': float(np.mean(rho_hat)), 'rho_std': float(np.std(rho_hat)), } return log_result('firm_predictions_structural_rho', tau, p_tau / 2 if sign_ok else 1.0, verd, details={ 'N_estimates': len(rho_hat), 'rho_mean': round(float(np.mean(rho_hat)), 3), 'rho_std': round(float(np.std(rho_hat)), 3), 'kendall_tau': round(tau, 4), 'kendall_p_two': round(p_tau, 4), 'kendall_p_one': round(p_tau / 2 if sign_ok else 1.0, 4), }, description=desc) _cache = {} def main(): results = [] raw, growth, leaf_df, dur_df = _load_growth_data() if leaf_df.shape[1] < 5: return [log_result('firm_predictions_variance_response', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {leaf_df.shape[1]} leaf sectors loaded'}, description="Early exit for all firm prediction tests. Fewer than 5 FRED IP leaf sectors could be loaded.")] # Test #4: FDT cross-section results.append(test_fdt(growth, leaf_df)) # Test #6: Equicorrelation results.append(test_equicorrelation(leaf_df)) # Test #7: Directed Coupling results.append(test_directed_coupling(dur_df)) # Test #8: Cyclical rho (correlation proxy) results.append(test_cyclical_rho(growth, dur_df)) # Test #8b: Structural rho (factor share identity) results.append(test_structural_rho()) _cache['results'] = results _cache['growth'] = growth _cache['leaf_df'] = leaf_df _cache['dur_df'] = dur_df return results def figures(): """Generate figures for firm-level prediction tests.""" if not _cache: main() paths = [] results = _cache.get('results', []) leaf_df = _cache.get('leaf_df') dur_df = _cache.get('dur_df') growth = _cache.get('growth', {}) import matplotlib.pyplot as plt # Figure 1: Variance-Response scatter (sigma^2 vs chi) fdt_data = _cache.get('fdt_data') fdt_ols = _cache.get('fdt_ols') if fdt_data is not None and fdt_ols is not None: sdf = fdt_data fig, ax = setup_figure(width=8, height=6) ax.scatter(sdf['chi'], sdf['sigma2'], c='steelblue', s=60, alpha=0.7, zorder=5) for _, row in sdf.iterrows(): ax.annotate(row['label'][:8], (row['chi'], row['sigma2']), fontsize=5, alpha=0.7, xytext=(3, 3), textcoords='offset points') # OLS fit line chi_range = np.linspace(sdf['chi'].min(), sdf['chi'].max(), 100) y_hat = fdt_ols['intercept'] + fdt_ols['T'] * chi_range ax.plot(chi_range, y_hat, 'r--', alpha=0.6, linewidth=1.5, label=f'OLS: T={fdt_ols["T"]:.2e}, R²={fdt_ols["r2"]:.3f}') ax.set_xlabel('Susceptibility χ (cumulative IRF to INDPRO shock)') ax.set_ylabel('Sector variance σ²') ax.set_title('Variance-Response Identity: σ² = α + T·χ\n(Cross-Section, N=17 FRED IP Sectors)') ax.legend(fontsize=8) paths.extend(save_figure(fig, 'firm_predictions_variance_response', subfolder='entry_exit')) # Figure 2: Equicorrelation heatmap eq_data = _cache.get('equicorr_data') if eq_data is not None: corr_matrix = eq_data['corr_matrix'] n = corr_matrix.shape[0] fig, ax = setup_figure(width=9, height=8) labels = [l[:10] for l in eq_data['labels']] im = ax.imshow(corr_matrix, cmap='RdBu_r', vmin=-0.5, vmax=1.0, aspect='auto') ax.set_xticks(range(n)) ax.set_yticks(range(n)) ax.set_xticklabels(labels, rotation=90, fontsize=6) ax.set_yticklabels(labels, fontsize=6) plt.colorbar(im, ax=ax, shrink=0.8, label='Correlation') ax.set_title(f'CES Equicorrelation: IP Sector Correlation Matrix\n' f'(mean r={eq_data["mean_corr"]:.3f}, CV={eq_data["cv"]:.2f})') paths.extend(save_figure(fig, 'firm_predictions_equicorrelation', subfolder='entry_exit')) # Figure 4: Directed coupling — antisymmetric vs BEA direction dir_data = _cache.get('directed_data') if dir_data is not None: od_a = dir_data['od_a'] od_d = dir_data['od_d'] fig, ax = setup_figure(width=8, height=6) # Jitter the discrete BEA direction for visibility rng = np.random.RandomState(42) jitter = rng.uniform(-0.1, 0.1, len(od_d)) colors = ['red' if d < 0 else ('green' if d > 0 else 'gray') for d in od_d] ax.scatter(od_d + jitter, od_a, c=colors, s=40, alpha=0.6) ax.axhline(0, color='k', linewidth=0.5, alpha=0.3) ax.axvline(0, color='k', linewidth=0.5, alpha=0.3) ax.set_xlabel('BEA Supply Chain Direction (-1/0/+1)') ax.set_ylabel('Antisymmetric IRF Component A[i,j]') r_sp = dir_data['r_sp'] p_sp = dir_data['p_sp'] ax.set_title(f'Directed Coupling: IRF Asymmetry vs BEA Direction\n' f'(Spearman r={r_sp:.3f}, p={p_sp:.4f})') paths.extend(save_figure(fig, 'firm_predictions_directed_coupling', subfolder='entry_exit')) # Figure 5: Cyclical rho — rolling correlation vs business cycle cyc_data = _cache.get('cyclical_data') if cyc_data is not None and 'INDPRO' in growth: proxy_df = cyc_data['proxy_df'] indpro_cycle = cyc_data['indpro_cycle'] tau = cyc_data['tau'] p_tau = cyc_data['p_tau'] fig, ax = setup_figure(width=10, height=5) ax.plot(proxy_df.index, proxy_df['avg_corr'], 'b-', linewidth=1.2, label='Avg pairwise correlation (60m window)') ax.set_ylabel('Average pairwise correlation', color='blue') ax.tick_params(axis='y', labelcolor='blue') ax2 = ax.twinx() cycle = indpro_cycle.dropna() ax2.plot(cycle.index, cycle.values, 'r-', linewidth=0.8, alpha=0.6, label='INDPRO growth (60m avg)') ax2.set_ylabel('INDPRO growth rate', color='red') ax2.tick_params(axis='y', labelcolor='red') ax2.invert_yaxis() ax.set_title(f'Cyclical Rho: Correlation Proxy vs Business Cycle\n' f'(Kendall tau={tau:.3f}, p={p_tau:.4f})') ax.legend(loc='upper left', fontsize=7) ax2.legend(loc='upper right', fontsize=7) paths.extend(save_figure(fig, 'firm_predictions_cyclical_rho', subfolder='entry_exit')) # Figure 6: Structural rho — rho_hat(t) vs output gap struct_data = _cache.get('structural_rho') if struct_data is not None: rho_s = struct_data['rho_series'] gap_s = struct_data['output_gap'] tau = struct_data['tau'] p_tau = struct_data['p_tau'] fig, ax = setup_figure(width=10, height=5) ax.plot(rho_s.index, rho_s.values, 'b-', linewidth=1.2, label=f'Structural rho (40q rolling, mean={struct_data["rho_mean"]:.2f})') ax.set_ylabel('Estimated rho', color='blue') ax.tick_params(axis='y', labelcolor='blue') ax.axhline(struct_data['rho_mean'], color='blue', linestyle=':', alpha=0.3) ax2 = ax.twinx() ax2.plot(gap_s.index, gap_s.values, 'r-', linewidth=0.8, alpha=0.6, label='Output gap (HP-filtered)') ax2.set_ylabel('Output gap (log deviation)', color='red') ax2.tick_params(axis='y', labelcolor='red') ax2.axhline(0, color='red', linestyle=':', alpha=0.3) ax.set_title(f'Structural Rho vs Output Gap\n' f'(Kendall tau={tau:.3f}, p_one={p_tau/2 if tau<0 else 1:.4f})') ax.legend(loc='upper left', fontsize=7) ax2.legend(loc='upper right', fontsize=7) paths.extend(save_figure(fig, 'firm_predictions_structural_rho', subfolder='entry_exit')) return paths if __name__ == '__main__': results = main() for r in results: print(f"{r['name']}: {r['verdict']} " f"(stat={r['statistic']}, p={r['p_value']})")

cv	0.3625
mean corr	0.3757
statistic	0.3625
frob ratio	0.2889
null cv median	18.0738

Firm Predictions Equicorrelation

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps