#!/usr/bin/env python3 """ Procurement PoIC Test -- Paper 2, Myerson Mechanism Design ============================================================ Tests the Price of Incentive Compatibility (PoIC) prediction: the markup (award/ceiling) in federal procurement should be higher for differentiated product categories (low rho) than for commodity categories (high rho). Sub-tests: A. Kruskal-Wallis: size-adjusted markup differs across categories B. Monotonicity: size-adjusted median markup LOW > HIGH C. OLS: substitutability coefficient negative controlling for log(ceiling) (HC1 robust SEs) C2. OLS with explicit contract-type control (is_cost_type dummy) D. Spearman rank correlation: substitutability vs adjusted median markup E. Category dummies OLS — reproduces Paper 5 Table Column (4) F. Fine-grained 3-digit NAICS OLS — reclassifies electronics/aerospace as LOW G. Bootstrap CI for LOW-HIGH adjusted markup gap Data: procurement_contracts.csv from USAspending.gov (1,200 contracts) Review Notes (2026-03-11) ------------------------- Theory match: CORRECT. CES effective curvature predicts information rents increase with input complementarity (lower rho → higher PoIC). Data-signal alignment: ADEQUATE after controlling for contract size. CRITICAL: Simpson's paradox. The raw bivariate regression has the WRONG sign (beta = +0.018) because large contracts cluster in differentiated sectors AND large contracts have lower award/ceiling ratios (cost-type contracts). Controlling for log(ceiling) flips the sign to beta = -0.057 (p = 0.0002). This is documented in Paper 5 Table (Column 1 vs Column 3). Contract type confound: fixed-price contracts (commodity goods) have award/ceiling ≈ 1.0 by construction; cost-type contracts (professional services) have ratio < 1.0 because the ceiling is a budget cap, not a price. The log(ceiling) control captures this. Strict monotonicity (LOW > MEDIUM > HIGH) fails even after adjustment because MEDIUM > LOW in the paper's own data. The key theoretical prediction is LOW > HIGH (complementary vs commodity), which holds. NAICS 2-digit limitation: ALL 393 LOW contracts come from NAICS 54 (Professional Services). The paper mentions electronics and aircraft as LOW substitutability, but these fall under NAICS 33 (Manufacturing, classified MEDIUM at 2-digit). The test effectively compares "professional services vs commodity sectors." OLS coefficient discrepancy: test gets beta=-0.057 (R²=0.051) while paper reports beta=-0.038 (R²=0.099) for Column (3). The paper's "Full Controls" likely include additional variables beyond log(ceiling). Both coefficients are negative and significant. No Lean theorem exists for PoIC — purely empirical prediction from Paper 2. """ import warnings import numpy as np import pandas as pd from scipy import stats from common.config import COMPILED_DIR, RESULTS_DIR from common.data import load_compiled from common.plotting import setup_figure, setup_subplots, save_figure from common.stats import verdict, log_result warnings.filterwarnings('ignore') # NAICS 2-digit substitutability classification NAICS_SUBSTITUTABILITY = { '11': 'HIGH', '21': 'HIGH', '22': 'HIGH', '42': 'HIGH', '44': 'HIGH', '45': 'HIGH', '23': 'MEDIUM', '31': 'MEDIUM', '32': 'MEDIUM', '33': 'MEDIUM', '48': 'MEDIUM', '49': 'MEDIUM', '51': 'MEDIUM', '52': 'MEDIUM', '53': 'MEDIUM', '56': 'MEDIUM', '61': 'MEDIUM', '71': 'MEDIUM', '72': 'MEDIUM', '54': 'LOW', '55': 'LOW', '62': 'LOW', '81': 'LOW', '92': 'LOW', } SUBST_NUMERIC = {'HIGH': 3, 'MEDIUM': 2, 'LOW': 1} # 3-digit NAICS overrides for finer-grained classification # Paper describes LOW as "electronics, aircraft, professional services" # At 2-digit, electronics (334) and aerospace (336) fall under NAICS 33 = MEDIUM # Petroleum products (324) are commodity → HIGH NAICS3_OVERRIDES = { '334': 'LOW', # Computer & Electronic Product Manufacturing '336': 'LOW', # Transportation Equipment Manufacturing (aerospace, ships) '324': 'HIGH', # Petroleum & Coal Products Manufacturing (commodity) } _ALL_SUBTESTS = ['procurement_poic_kruskal_wallis', 'procurement_poic_monotonicity', 'procurement_poic_ols_subst', 'procurement_poic_ols_contract_type', 'procurement_poic_spearman', 'procurement_poic_category_dummies', 'procurement_poic_fine_naics', 'procurement_poic_monotonicity_boot'] def load_procurement_data(): """Load USAspending procurement data from compiled directory.""" try: df = load_compiled('procurement_contracts') except FileNotFoundError: return None # Compute markup if 'markup' not in df.columns: if 'current_total_value_of_award' in df.columns and 'potential_total_value_of_award' in df.columns: df['award_value'] = pd.to_numeric(df['current_total_value_of_award'], errors='coerce') df['ceiling_value'] = pd.to_numeric(df['potential_total_value_of_award'], errors='coerce') elif 'award_value' not in df.columns or 'ceiling_value' not in df.columns: return None df = df[(df['ceiling_value'] > 1000) & (df['award_value'] > 0)].copy() df['markup'] = df['award_value'] / df['ceiling_value'] df = df[(df['markup'] > 0.01) & (df['markup'] < 10)].copy() # Classify NAICS if 'subst_category' not in df.columns: naics_col = 'naics_code' if 'naics_code' in df.columns else 'naics' if naics_col not in df.columns: return None df['naics2'] = df[naics_col].astype(str).str[:2] df['subst_category'] = df['naics2'].map(NAICS_SUBSTITUTABILITY) df = df.dropna(subset=['subst_category']).copy() df['subst_numeric'] = df['subst_category'].map(SUBST_NUMERIC) # Fine-grained 3-digit classification (overrides 2-digit for specific subsectors) naics_col = 'naics_code' if 'naics_code' in df.columns else 'naics' df['naics3'] = df[naics_col].astype(str).str[:3] df['subst_fine'] = df['subst_category'].copy() for n3, cat in NAICS3_OVERRIDES.items(): mask = df['naics3'] == n3 df.loc[mask, 'subst_fine'] = cat df['subst_fine_numeric'] = df['subst_fine'].map(SUBST_NUMERIC) # Classify contract pricing type (fixed-price vs cost-type) if 'type_of_contract_pricing' in df.columns: fp_codes = {'J', 'K', 'L', 'A'} # Firm FP, FP EPA, FP Incentive, FP Redet cost_codes = {'R', 'S', 'U', 'V'} # Cost Plus Award/No/Fixed/Incentive Fee df['is_cost_type'] = df['type_of_contract_pricing'].isin(cost_codes).astype(int) # Compute log(ceiling) for size control and size-adjusted markup if 'ceiling_value' in df.columns: df['log_ceiling'] = np.log(df['ceiling_value']) # Residualize markup on log(ceiling) to remove Simpson's paradox try: import statsmodels.api as sm Xc = sm.add_constant(df['log_ceiling']) mc = sm.OLS(df['markup'], Xc).fit() df['adj_markup'] = mc.resid + mc.params['const'] except Exception: df['adj_markup'] = df['markup'] else: df['adj_markup'] = df['markup'] return df def test_kruskal_wallis(df): """Test A: Size-adjusted markup differs across substitutability categories.""" desc = ("Tests whether size-adjusted procurement markups differ significantly " "across substitutability categories (LOW/MEDIUM/HIGH) AND in the " "predicted direction (LOW median > HIGH median). Kruskal-Wallis H-test " "on USAspending contract markups residualized on log(ceiling).") if df is None: return log_result('procurement_poic_kruskal_wallis', np.nan, np.nan, 'SKIPPED', description=desc) groups = {c: df.loc[df['subst_category'] == c, 'adj_markup'].values for c in ['LOW', 'MEDIUM', 'HIGH']} groups_list = [g for g in groups.values() if len(g) >= 10] if len(groups_list) < 3: return log_result('procurement_poic_kruskal_wallis', np.nan, np.nan, 'SKIPPED', description=desc) H, p_kw = stats.kruskal(*groups_list) # Check direction: LOW median should be > HIGH median low_med = np.median(groups['LOW']) if len(groups['LOW']) >= 10 else np.nan high_med = np.median(groups['HIGH']) if len(groups['HIGH']) >= 10 else np.nan sign_ok = low_med > high_med verd = verdict(p_kw, sign_ok, threshold=0.10) return log_result('procurement_poic_kruskal_wallis', H, p_kw, verd, details={'n_groups': 3, 'adj_med_LOW': round(low_med, 4), 'adj_med_HIGH': round(high_med, 4)}, description=desc) def test_monotonicity(df): """Test B: Size-adjusted median markup LOW > HIGH (key prediction). Strict monotonicity (LOW > MEDIUM > HIGH) fails even in the paper's own data (MEDIUM > LOW after controls). The key theoretical prediction is LOW > HIGH: complementary goods have higher information rents than commodity goods. Full monotonicity is reported but not required. """ desc = ("Tests the PoIC directional prediction: size-adjusted median " "markup for LOW substitutability (complementary) should exceed " "HIGH substitutability (commodity). Uses markups residualized on " "log(ceiling) to remove the Simpson's paradox from contract type. " "CONSISTENT if LOW > HIGH. Full monotonicity (LOW > MEDIUM > HIGH) " "is reported but not required — the paper's own data shows MEDIUM > LOW.") if df is None: return log_result('procurement_poic_monotonicity', np.nan, np.nan, 'SKIPPED', description=desc) medians = {} for cat in ['LOW', 'MEDIUM', 'HIGH']: vals = df.loc[df['subst_category'] == cat, 'adj_markup'] if len(vals) < 10: return log_result('procurement_poic_monotonicity', np.nan, np.nan, 'SKIPPED', description=desc) medians[cat] = vals.median() low_gt_high = medians['LOW'] > medians['HIGH'] full_monotonic = medians['LOW'] > medians['MEDIUM'] > medians['HIGH'] gap = medians['LOW'] - medians['HIGH'] if low_gt_high: verd = 'CONSISTENT' else: verd = 'INCONSISTENT' return log_result('procurement_poic_monotonicity', gap, np.nan, verd, details={'adj_LOW': round(medians['LOW'], 4), 'adj_MEDIUM': round(medians['MEDIUM'], 4), 'adj_HIGH': round(medians['HIGH'], 4), 'full_monotonic': full_monotonic}, description=desc) def test_ols_substitutability(df): """Test C: OLS — substitutability predicts lower markup, controlling for log(ceiling). Matches Paper 5 Table Column (3). Without the log(ceiling) control, the bivariate regression has a POSITIVE coefficient (Simpson's paradox) because large cost-type contracts cluster in differentiated sectors. """ desc = ("Tests that higher substitutability predicts lower procurement markup, " "controlling for log(ceiling) to resolve Simpson's paradox. OLS with " "HC1 robust SEs: markup ~ subst_numeric + log(ceiling). Inspired by " "Paper 5 Table Column (3), but yields beta=-0.057 vs paper's -0.038 " "(paper may include additional controls not shown). " "Expects negative substitutability coefficient.") if df is None: return log_result('procurement_poic_ols_subst', np.nan, np.nan, 'SKIPPED', description=desc) try: import statsmodels.api as sm X = sm.add_constant(df[['subst_numeric', 'log_ceiling']]) y = df['markup'] model = sm.OLS(y, X).fit(cov_type='HC1') b = model.params['subst_numeric'] p = model.pvalues['subst_numeric'] sign_ok = b < 0 # higher substitutability => lower markup verd = verdict(p, sign_ok, threshold=0.10) return log_result('procurement_poic_ols_subst', b, p, verd, details={'n': int(model.nobs), 'r2': model.rsquared, 'beta_log_ceiling': float(model.params['log_ceiling'])}, description=desc) except Exception: return log_result('procurement_poic_ols_subst', np.nan, np.nan, 'SKIPPED', description=desc) def test_ols_contract_type(df): """Test C2: OLS with explicit contract-type control. Adds is_cost_type dummy to the OLS specification. Cost-type contracts (award/ceiling < 1.0 by construction) cluster in differentiated sectors, creating Simpson's paradox. This is a cleaner control than log(ceiling) alone because it directly captures the pricing mechanism confound. """ desc = ("Tests substitutability effect on markup with explicit contract-type control. " "OLS with HC1 robust SEs: markup ~ subst_numeric + log(ceiling) + is_cost_type. " "Cost-type contracts have award/ceiling < 1.0 by construction (ceiling = budget cap). " "Expects negative substitutability coefficient after controlling for both " "contract size and pricing mechanism.") if df is None or 'is_cost_type' not in df.columns: return log_result('procurement_poic_ols_contract_type', np.nan, np.nan, 'SKIPPED', description=desc) try: import statsmodels.api as sm X = sm.add_constant(df[['subst_numeric', 'log_ceiling', 'is_cost_type']]) y = df['markup'] model = sm.OLS(y, X).fit(cov_type='HC1') b = model.params['subst_numeric'] p = model.pvalues['subst_numeric'] sign_ok = b < 0 verd = verdict(p, sign_ok, threshold=0.10) return log_result('procurement_poic_ols_contract_type', b, p, verd, details={'n': int(model.nobs), 'r2': model.rsquared, 'beta_log_ceiling': float(model.params['log_ceiling']), 'beta_cost_type': float(model.params['is_cost_type']), 'p_cost_type': float(model.pvalues['is_cost_type'])}, description=desc) except Exception: return log_result('procurement_poic_ols_contract_type', np.nan, np.nan, 'SKIPPED', description=desc) def test_spearman_sector(df): """Test D: Spearman rank correlation at sector level (EXPLORATORY). With only 9 NAICS2 sectors in the data (1 LOW, 5 MEDIUM, 3 HIGH), this test is severely underpowered. Capped at AMBIGUOUS regardless of result because N < 12. """ desc = ("EXPLORATORY: Spearman rank correlation between NAICS 2-digit " "substitutability score and sector median size-adjusted markup. " "Capped at AMBIGUOUS due to only 9 sectors (1 LOW). Negative " "correlation means more substitutable sectors have lower markups.") if df is None: return log_result('procurement_poic_spearman', np.nan, np.nan, 'SKIPPED', description=desc) sector_stats = df.groupby(['naics2', 'subst_category']).agg( median_markup=('adj_markup', 'median'), n=('adj_markup', 'count'), ).reset_index() sector_stats = sector_stats[sector_stats['n'] >= 20] sector_stats['subst_num'] = sector_stats['subst_category'].map(SUBST_NUMERIC) if len(sector_stats) < 5: return log_result('procurement_poic_spearman', np.nan, np.nan, 'SKIPPED', description=desc) rho, p = stats.spearmanr(sector_stats['subst_num'], sector_stats['median_markup']) sign_ok = rho < 0 # more substitutable => lower markup # Cap at AMBIGUOUS: only ~9 sectors with 1 LOW — underpowered if len(sector_stats) < 12: verd = 'AMBIGUOUS' else: verd = verdict(p, sign_ok, threshold=0.10) return log_result('procurement_poic_spearman', rho, p, verd, details={'n_sectors': len(sector_stats)}, description=desc) def test_category_dummies(df): """Test E: Category dummy OLS — reproduces Paper 5 Table Column (4). Uses LOW/HIGH category dummies (ref=HIGH) with log(ceiling) control. This is the paper's strongest specification (R²=0.227) and directly estimates the LOW-HIGH gap in percentage points. """ desc = ("Reproduces Paper 5 Table Column (4): OLS with category dummies " "(LOW/MEDIUM vs HIGH reference) and log(ceiling) control. HC1 robust SEs. " "markup ~ LOW + MEDIUM + log(ceiling). Paper reports LOW=+23.3pp (p<10^-14), " "MEDIUM=+37.3pp, R²=0.227. Directly estimates the PoIC category gaps.") if df is None: return log_result('procurement_poic_category_dummies', np.nan, np.nan, 'SKIPPED', description=desc) try: import statsmodels.api as sm # Create dummy variables manually to avoid dtype issues X_df = pd.DataFrame({ 'const': 1.0, 'log_ceiling': df['log_ceiling'].astype(float), 'cat_LOW': (df['subst_category'] == 'LOW').astype(float), 'cat_MEDIUM': (df['subst_category'] == 'MEDIUM').astype(float), }) y = df['markup'].astype(float) model = sm.OLS(y, X_df).fit(cov_type='HC1') b_low = model.params.get('cat_LOW', np.nan) p_low = model.pvalues.get('cat_LOW', np.nan) b_med = model.params.get('cat_MEDIUM', np.nan) sign_ok = b_low > 0 # LOW should have higher markup than HIGH (reference) verd = verdict(p_low, sign_ok, threshold=0.10) _cache['cat_dummy_model'] = model return log_result('procurement_poic_category_dummies', float(b_low), float(p_low), verd, details={ 'n': int(model.nobs), 'r2': float(model.rsquared), 'beta_LOW': float(b_low), 'beta_MEDIUM': float(b_med), 'beta_log_ceiling': float(model.params['log_ceiling']), 'gap_LOW_HIGH_pp': float(b_low * 100), }, description=desc) except Exception as e: return log_result('procurement_poic_category_dummies', np.nan, np.nan, 'SKIPPED', description=desc, details={'error': str(e)}) def test_fine_naics(df): """Test F: OLS with fine-grained 3-digit NAICS classification. Reclassifies manufacturing subsectors using 3-digit NAICS: 334 (Electronics Mfg) → LOW 336 (Aerospace/Defense Mfg) → LOW 324 (Petroleum Products) → HIGH This matches the paper's description of LOW as "electronics, aircraft, professional services" and diversifies LOW from 1 sector to 3. """ desc = ("Tests PoIC with fine-grained 3-digit NAICS substitutability classification. " "Reclassifies NAICS 334 (electronics) and 336 (aerospace) from MEDIUM to LOW, " "and 324 (petroleum) from MEDIUM to HIGH, matching the paper's sector descriptions. " "OLS with HC1 robust SEs: markup ~ subst_fine_numeric + log(ceiling). " "Expects negative substitutability coefficient. LOW now has 3 sectors " "(professional services + electronics + aerospace) instead of 1.") if df is None or 'subst_fine_numeric' not in df.columns: return log_result('procurement_poic_fine_naics', np.nan, np.nan, 'SKIPPED', description=desc) try: import statsmodels.api as sm X = sm.add_constant(df[['subst_fine_numeric', 'log_ceiling']]) y = df['markup'] model = sm.OLS(y, X).fit(cov_type='HC1') b = model.params['subst_fine_numeric'] p = model.pvalues['subst_fine_numeric'] sign_ok = b < 0 verd = verdict(p, sign_ok, threshold=0.10) # Compute category stats fine_counts = df['subst_fine'].value_counts().to_dict() fine_n3 = df.groupby('subst_fine')['naics3'].nunique().to_dict() _cache['fine_naics_model'] = model _cache['fine_naics_df'] = df return log_result('procurement_poic_fine_naics', float(b), float(p), verd, details={ 'n': int(model.nobs), 'r2': float(model.rsquared), 'beta_subst': float(b), 'beta_log_ceiling': float(model.params['log_ceiling']), 'category_counts': fine_counts, 'n_sectors_per_cat': fine_n3, }, description=desc) except Exception as e: return log_result('procurement_poic_fine_naics', np.nan, np.nan, 'SKIPPED', description=desc, details={'error': str(e)}) def test_monotonicity_bootstrap(df): """Test G: Bootstrap CI for the LOW-HIGH adjusted markup gap. The monotonicity test (Test B) reports the gap without uncertainty. This adds a bootstrap 95% CI to quantify whether the LOW > HIGH direction is robust to sampling variability. """ desc = ("Bootstrap 95% CI for the LOW-HIGH size-adjusted markup gap. " "Resamples contracts within each category 2000 times and computes " "median(adj_markup_LOW) - median(adj_markup_HIGH). Tests whether " "the gap is significantly positive. Data: USAspending contracts, " "markups residualized on log(ceiling).") if df is None: return log_result('procurement_poic_monotonicity_boot', np.nan, np.nan, 'SKIPPED', description=desc) low_vals = df.loc[df['subst_category'] == 'LOW', 'adj_markup'].values high_vals = df.loc[df['subst_category'] == 'HIGH', 'adj_markup'].values if len(low_vals) < 20 or len(high_vals) < 20: return log_result('procurement_poic_monotonicity_boot', np.nan, np.nan, 'SKIPPED', description=desc) rng = np.random.RandomState(42) n_boot = 2000 gaps = np.empty(n_boot) for i in range(n_boot): low_sample = rng.choice(low_vals, size=len(low_vals), replace=True) high_sample = rng.choice(high_vals, size=len(high_vals), replace=True) gaps[i] = np.median(low_sample) - np.median(high_sample) observed_gap = np.median(low_vals) - np.median(high_vals) ci_low = np.percentile(gaps, 2.5) ci_high = np.percentile(gaps, 97.5) # Bootstrap p-value: fraction of bootstrap gaps <= 0 p_boot = np.mean(gaps <= 0) sign_ok = observed_gap > 0 if ci_low > 0 and sign_ok: verd = 'CONSISTENT' elif sign_ok: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' _cache['boot_gaps'] = gaps _cache['boot_observed'] = observed_gap return log_result('procurement_poic_monotonicity_boot', float(observed_gap), float(p_boot), verd, details={ 'observed_gap_pp': float(observed_gap * 100), 'ci_95': [float(ci_low), float(ci_high)], 'ci_95_pp': [float(ci_low * 100), float(ci_high * 100)], 'n_boot': n_boot, 'n_low': len(low_vals), 'n_high': len(high_vals), 'pct_positive': float(np.mean(gaps > 0) * 100), }, description=desc) _cache = {} def main(): results = [] df = load_procurement_data() results.append(test_kruskal_wallis(df)) results.append(test_monotonicity(df)) results.append(test_ols_substitutability(df)) results.append(test_ols_contract_type(df)) results.append(test_spearman_sector(df)) results.append(test_category_dummies(df)) results.append(test_fine_naics(df)) results.append(test_monotonicity_bootstrap(df)) _cache['results'] = results _cache['df'] = df return results def figures(): """Generate procurement PoIC figures — one per subtest.""" if not _cache: main() paths = [] df = _cache.get('df') if df is None or len(df) == 0: return paths colors_cat = {'LOW': '#d73027', 'MEDIUM': '#fee08b', 'HIGH': '#1a9850'} # Figure 1: Box plot — Kruskal-Wallis subtest fig, ax = setup_figure() groups = [] labels = [] colors_box = [] for cat in ['LOW', 'MEDIUM', 'HIGH']: vals = df.loc[df['subst_category'] == cat, 'adj_markup'].values if len(vals) >= 10: groups.append(vals) labels.append(cat) colors_box.append(colors_cat[cat]) if len(groups) >= 2: bp = ax.boxplot(groups, tick_labels=labels, patch_artist=True, showmeans=True, meanline=True, meanprops=dict(color='red', linewidth=2), showfliers=False) for patch, color in zip(bp['boxes'], colors_box): patch.set_facecolor(color) patch.set_alpha(0.7) ax.set_xlabel('Substitutability Category') ax.set_ylabel('Size-Adjusted Markup (residualized on log ceiling)') ax.set_title('Size-Adjusted Markup by Substitutability Category') ax.axhline(1.0, color='gray', linestyle='--', alpha=0.5) paths += save_figure(fig, 'procurement_poic_boxplot', subfolder='applications') # Figure 2: Bar chart of medians — Monotonicity subtest fig, ax = setup_figure() medians = {} for cat in ['LOW', 'MEDIUM', 'HIGH']: vals = df.loc[df['subst_category'] == cat, 'adj_markup'] if len(vals) >= 10: medians[cat] = vals.median() if medians: cats = list(medians.keys()) vals = [medians[c] for c in cats] bars = ax.bar(cats, vals, color=[colors_cat[c] for c in cats], edgecolor='black', alpha=0.7) ax.set_ylabel('Median Size-Adjusted Markup') ax.set_xlabel('Substitutability Category') ax.set_title('Median Markup: LOW > HIGH (PoIC Prediction)') ax.axhline(1.0, color='gray', linestyle='--', alpha=0.5) for bar, v in zip(bars, vals): ax.text(bar.get_x() + bar.get_width()/2, v + 0.005, f'{v:.3f}', ha='center', va='bottom', fontsize=10) paths += save_figure(fig, 'procurement_poic_competition', subfolder='applications') # Figure 3: OLS coefficients — OLS subtest try: import statsmodels.api as sm X = sm.add_constant(df[['subst_numeric', 'log_ceiling']]) model = sm.OLS(df['markup'], X).fit(cov_type='HC1') fig, ax = setup_figure() names = ['subst_numeric', 'log_ceiling'] coefs = [model.params[n] for n in names] cis = [model.conf_int().loc[n] for n in names] yerr = [[c - ci[0] for c, ci in zip(coefs, cis)], [ci[1] - c for c, ci in zip(coefs, cis)]] colors_coef = ['#d73027' if c < 0 else '#1a9850' for c in coefs] ax.barh(names, coefs, xerr=yerr, color=colors_coef, alpha=0.7, edgecolor='black', capsize=5) ax.axvline(0, color='black', linewidth=0.8) ax.set_xlabel('OLS Coefficient (HC1 robust)') ax.set_title(f'OLS: markup ~ subst + log(ceiling) (N={int(model.nobs)})') for i, (n, c) in enumerate(zip(names, coefs)): p = model.pvalues[n] ax.text(c, i, f' {c:.4f} (p={p:.4f})', va='center', fontsize=9) paths += save_figure(fig, 'procurement_poic_coefficients', subfolder='applications') except Exception: pass # Figure 4: OLS with contract-type control — OLS contract type subtest try: import statsmodels.api as sm if 'is_cost_type' in df.columns: X2 = sm.add_constant(df[['subst_numeric', 'log_ceiling', 'is_cost_type']]) model2 = sm.OLS(df['markup'], X2).fit(cov_type='HC1') fig, ax = setup_figure() names2 = ['subst_numeric', 'log_ceiling', 'is_cost_type'] coefs2 = [model2.params[n] for n in names2] cis2 = [model2.conf_int().loc[n] for n in names2] yerr2 = [[c - ci[0] for c, ci in zip(coefs2, cis2)], [ci[1] - c for c, ci in zip(coefs2, cis2)]] colors_coef2 = ['#d73027' if c < 0 else '#1a9850' for c in coefs2] ax.barh(names2, coefs2, xerr=yerr2, color=colors_coef2, alpha=0.7, edgecolor='black', capsize=5) ax.axvline(0, color='black', linewidth=0.8) ax.set_xlabel('OLS Coefficient (HC1 robust)') ax.set_title(f'OLS with Contract Type Control (N={int(model2.nobs)})') for i, (n, c) in enumerate(zip(names2, coefs2)): p = model2.pvalues[n] ax.text(c, i, f' {c:.4f} (p={p:.4f})', va='center', fontsize=9) paths += save_figure(fig, 'procurement_poic_contract_type', subfolder='applications') except Exception: pass # Figure 5: Category dummies coefficient plot cat_model = _cache.get('cat_dummy_model') if cat_model is not None: try: fig, ax = setup_figure() names = ['cat_LOW', 'cat_MEDIUM', 'log_ceiling'] coefs = [cat_model.params[n] for n in names] cis = [cat_model.conf_int().loc[n] for n in names] yerr = [[c - ci[0] for c, ci in zip(coefs, cis)], [ci[1] - c for c, ci in zip(coefs, cis)]] labels = ['LOW (vs HIGH)', 'MEDIUM (vs HIGH)', 'log(ceiling)'] colors_c = ['#d73027' if 'LOW' in n else '#fee08b' if 'MED' in n else '#1a9850' for n in names] ax.barh(labels, coefs, xerr=yerr, color=colors_c, alpha=0.7, edgecolor='black', capsize=5) ax.axvline(0, color='black', linewidth=0.8) ax.set_xlabel('OLS Coefficient (HC1 robust)') ax.set_title(f'Category Dummies (Paper 5 Column 4)\n' f'LOW = +{coefs[0]*100:.1f}pp, R² = {cat_model.rsquared:.3f}, ' f'N = {int(cat_model.nobs)}') for i, (l, c) in enumerate(zip(labels, coefs)): p = cat_model.pvalues[names[i]] ax.text(c, i, f' {c:.3f} (p={p:.1e})', va='center', fontsize=9) paths += save_figure(fig, 'procurement_poic_category_dummies', subfolder='applications') except Exception: pass # Figure 6: Fine-grained NAICS comparison fine_model = _cache.get('fine_naics_model') fine_df = _cache.get('fine_naics_df') if fine_df is not None: fig, (ax1, ax2) = setup_subplots(nrows=1, ncols=2, width=14, height=5) # Left: box plot by fine category groups_f = [] labels_f = [] colors_f = [] for cat in ['LOW', 'MEDIUM', 'HIGH']: vals = fine_df.loc[fine_df['subst_fine'] == cat, 'adj_markup'].values if len(vals) >= 10: groups_f.append(vals) labels_f.append(f'{cat}\n(n={len(vals)})') colors_f.append(colors_cat[cat]) if len(groups_f) >= 2: bp = ax1.boxplot(groups_f, tick_labels=labels_f, patch_artist=True, showmeans=True, meanline=True, showfliers=False, meanprops=dict(color='red', linewidth=2)) for patch, color in zip(bp['boxes'], colors_f): patch.set_facecolor(color) patch.set_alpha(0.7) ax1.set_ylabel('Size-Adjusted Markup') ax1.set_title('Fine-Grained NAICS (3-digit overrides)') ax1.axhline(1.0, color='gray', linestyle='--', alpha=0.5) # Right: OLS coefficients if fine_model is not None: names_fm = ['subst_fine_numeric', 'log_ceiling'] coefs_fm = [fine_model.params[n] for n in names_fm] cis_fm = [fine_model.conf_int().loc[n] for n in names_fm] yerr_fm = [[c - ci[0] for c, ci in zip(coefs_fm, cis_fm)], [ci[1] - c for c, ci in zip(coefs_fm, cis_fm)]] colors_fm = ['#d73027' if c < 0 else '#1a9850' for c in coefs_fm] ax2.barh(names_fm, coefs_fm, xerr=yerr_fm, color=colors_fm, alpha=0.7, edgecolor='black', capsize=5) ax2.axvline(0, color='black', linewidth=0.8) ax2.set_xlabel('OLS Coefficient (HC1 robust)') ax2.set_title(f'Fine NAICS OLS (R² = {fine_model.rsquared:.3f})') for i, (n, c) in enumerate(zip(names_fm, coefs_fm)): p = fine_model.pvalues[n] ax2.text(c, i, f' {c:.4f} (p={p:.4f})', va='center', fontsize=9) fig.tight_layout() paths += save_figure(fig, 'procurement_poic_fine_naics', subfolder='applications') # Figure 7: Bootstrap gap distribution boot_gaps = _cache.get('boot_gaps') boot_obs = _cache.get('boot_observed') if boot_gaps is not None: fig, ax = setup_figure() ax.hist(boot_gaps * 100, bins=40, color='steelblue', alpha=0.7, edgecolor='black') ax.axvline(0, color='black', linestyle='-', linewidth=1.5) ax.axvline(boot_obs * 100, color='red', linestyle='--', linewidth=2, label=f'Observed = {boot_obs*100:.1f}pp') ci_lo = np.percentile(boot_gaps, 2.5) * 100 ci_hi = np.percentile(boot_gaps, 97.5) * 100 ax.axvspan(ci_lo, ci_hi, alpha=0.15, color='red', label=f'95% CI [{ci_lo:.1f}, {ci_hi:.1f}]pp') ax.set_xlabel('LOW - HIGH median adjusted markup gap (pp)') ax.set_ylabel('Bootstrap count') pct_pos = np.mean(boot_gaps > 0) * 100 ax.set_title(f'Bootstrap: LOW > HIGH gap = {boot_obs*100:.1f}pp\n' f'95% CI [{ci_lo:.1f}, {ci_hi:.1f}]pp, {pct_pos:.0f}% positive') ax.legend(fontsize=9) paths += save_figure(fig, 'procurement_poic_bootstrap', subfolder='applications') # Figure 8: Sector scatter — Spearman subtest if 'naics2' in df.columns: sector = df.groupby(['naics2', 'subst_category']).agg( median_markup=('adj_markup', 'median'), n=('adj_markup', 'count')).reset_index() sector = sector[sector['n'] >= 20] sector['subst_num'] = sector['subst_category'].map(SUBST_NUMERIC) if len(sector) >= 5: fig, ax = setup_figure() for cat in ['LOW', 'MEDIUM', 'HIGH']: sub = sector[sector['subst_category'] == cat] ax.scatter(sub['subst_num'], sub['median_markup'], c=colors_cat.get(cat, 'gray'), edgecolors='black', linewidth=0.5, s=80, alpha=0.7, label=cat, zorder=3) for _, row in sub.iterrows(): ax.annotate(row['naics2'], (row['subst_num'], row['median_markup']), textcoords='offset points', xytext=(5, 5), fontsize=8) rho, p = stats.spearmanr(sector['subst_num'], sector['median_markup']) ax.set_xlabel('Substitutability Score (1=LOW, 2=MEDIUM, 3=HIGH)') ax.set_ylabel('Sector Median Size-Adjusted Markup') ax.set_title(f'Sector Spearman: rho={rho:.3f}, p={p:.3f} (N={len(sector)})') ax.legend(fontsize=9) paths += save_figure(fig, 'procurement_poic_scatter', subfolder='applications') 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']})")

p value	1.562e-52
n groups	3
statistic	238.5772
adj med LOW	1.3980
adj med HIGH	1.3447

Procurement Poic Kruskal Wallis

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps