#!/usr/bin/env python3 """ Laffer Curve Test — CES Complementarity Predictions ==================================================== Tests CES prediction that sector friction sensitivity depends on complementarity: Test 1: Panel interaction — GDP growth response to tax changes varies by curvature K Test 2: TCJA SALT cap event study — high-tax states hit harder, with sector interaction Test 3: Monetary tightening sensitivity — IP response to FFR changes varies by input HHI Theory: Paper 2 (CES Potential), Section "Laffer Curve as CES Potential Result". CES curvature K = (1-rho)(J-1)/J amplifies welfare loss from friction. Sectors with higher input concentration (higher K proxy) suffer larger output losses from friction shocks. Lean: lafferPeak(eta) = 1/(2*eta), ramsey_inverse_elasticity (GrowthTax.lean). Data: Tests 1-2: FRED state GDP by industry (50 states x 6 sectors), Tax Foundation rates. Test 3: FRED IP subsector indices (18 mfg sectors, monthly 1972-2026), DFF (daily Fed Funds Rate → monthly), BEA I-O Use Table 2022 (input concentration HHI). Review Notes (2026-03-11, v3): Theory match: PARTIAL — CES curvature prediction: high-K sectors more friction-sensitive. Tests 1-2 have proxy failure (state income tax ≠ manufacturing friction). Test 3 uses monetary tightening as universal friction + BEA I-O input HHI as K proxy. I/O audit: CORRECT — FFR changes as friction (affects all sectors via cost of capital); IP indices as output (direct production measure); input HHI from BEA I-O as K proxy. Data-signal alignment: Test 3 STRONG — 18 sectors, monthly data, wide HHI variation (0.056 to 0.634), frequent rate changes. Tests 1-2 WEAK (sparse tax variation). Code quality: CLEAN. Confidence in verdict: Test 3 MEDIUM-HIGH; Tests 1-2 LOW-MEDIUM. Migrated from scripts/test_laffer_curve.py """ 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 HAS_SM = True except ImportError: HAS_SM = False # Sector complementarity — rho from Oberfield-Raval (2021) and Atalay (2017) # K = (1-rho) for large J approximation SECTOR_RHO = { 'Construction': {'fred_code': 'CONSTNQGSP', 'naics': '23', 'rho': -0.073}, 'Manufacturing': {'fred_code': 'MANNQGSP', 'naics': '31-33', 'rho': -0.299}, 'Healthcare': {'fred_code': 'HLTHSOCASSNQGSP', 'naics': '62', 'rho': -0.073}, 'Finance': {'fred_code': 'FININSNQGSP', 'naics': '52', 'rho': -0.073}, 'Information': {'fred_code': 'INFONQGSP', 'naics': '51', 'rho': -0.191}, 'ProfServices': {'fred_code': 'PROBUSNQGSP', 'naics': '54', 'rho': -0.073}, } STATES = [ 'AL', 'AK', 'AZ', 'AR', 'CA', 'CO', 'CT', 'DE', 'FL', 'GA', 'HI', 'ID', 'IL', 'IN', 'IA', 'KS', 'KY', 'LA', 'ME', 'MD', 'MA', 'MI', 'MN', 'MS', 'MO', 'MT', 'NE', 'NV', 'NH', 'NJ', 'NM', 'NY', 'NC', 'ND', 'OH', 'OK', 'OR', 'PA', 'RI', 'SC', 'SD', 'TN', 'TX', 'UT', 'VT', 'VA', 'WA', 'WV', 'WI', 'WY', ] NO_INCOME_TAX = {'AK', 'FL', 'NV', 'NH', 'SD', 'TN', 'TX', 'WA', 'WY'} TAX_SCHEDULE = { 'AL': [(2005, 5.0)], 'AZ': [(2005, 4.54), (2021, 2.98), (2023, 2.5)], 'AR': [(2005, 7.0), (2015, 6.9), (2019, 5.9), (2023, 4.4), (2024, 3.9)], 'CA': [(2005, 9.3), (2009, 10.3), (2012, 13.3)], 'CO': [(2005, 4.63), (2020, 4.55), (2022, 4.4)], 'CT': [(2005, 5.0), (2009, 6.5), (2011, 6.7), (2015, 6.99)], 'DE': [(2005, 5.95), (2010, 6.6)], 'GA': [(2005, 6.0), (2024, 5.49)], 'HI': [(2005, 8.25), (2009, 11.0)], 'ID': [(2005, 7.8), (2012, 6.925), (2018, 6.5), (2022, 6.0), (2023, 5.8)], 'IL': [(2005, 3.0), (2011, 5.0), (2015, 3.75), (2017, 4.95)], 'IN': [(2005, 3.4), (2015, 3.3), (2017, 3.23), (2023, 3.15), (2024, 3.05)], 'IA': [(2005, 8.98), (2023, 6.0), (2024, 5.7)], 'KS': [(2005, 6.45), (2013, 3.9), (2017, 5.7)], 'KY': [(2005, 6.0), (2018, 5.0), (2023, 4.5), (2024, 4.0)], 'LA': [(2005, 6.0), (2022, 4.25)], 'ME': [(2005, 8.5), (2016, 7.15)], 'MD': [(2005, 5.5), (2008, 6.25), (2012, 5.75)], 'MA': [(2005, 5.3), (2012, 5.15), (2014, 5.1), (2016, 5.05), (2020, 5.0)], 'MI': [(2005, 3.9), (2007, 4.35), (2012, 4.25)], 'MN': [(2005, 7.85), (2013, 9.85)], 'MS': [(2005, 5.0), (2024, 4.7)], 'MO': [(2005, 6.0), (2019, 5.4), (2020, 5.3), (2022, 4.95), (2024, 4.8)], 'MT': [(2005, 6.9), (2024, 5.9)], 'NE': [(2005, 6.84), (2024, 5.84)], 'NJ': [(2005, 8.97), (2018, 10.75)], 'NM': [(2005, 4.9), (2021, 5.9)], 'NY': [(2005, 6.85), (2009, 8.82), (2021, 10.9)], 'NC': [(2005, 8.25), (2014, 5.8), (2015, 5.75), (2017, 5.499), (2019, 5.25), (2022, 4.99), (2023, 4.75), (2024, 4.5)], 'ND': [(2005, 5.54), (2009, 4.86), (2013, 3.22), (2019, 2.9), (2024, 1.95)], 'OH': [(2005, 7.5), (2008, 5.925), (2014, 5.0), (2019, 4.797), (2021, 3.99), (2024, 3.5)], 'OK': [(2005, 6.65), (2008, 5.5), (2016, 5.0), (2022, 4.75)], 'OR': [(2005, 9.0), (2010, 9.9)], 'PA': [(2005, 3.07)], 'RI': [(2005, 9.9), (2011, 5.99)], 'SC': [(2005, 7.0), (2022, 6.5), (2024, 6.4)], 'UT': [(2005, 7.0), (2007, 5.0), (2023, 4.85), (2024, 4.65)], 'VT': [(2005, 9.5), (2007, 8.95), (2009, 8.75)], 'VA': [(2005, 5.75)], 'WV': [(2005, 6.5), (2024, 5.12)], 'WI': [(2005, 6.75), (2013, 7.65)], } def get_tax_rate(state, year): if state in NO_INCOME_TAX: return 0.0 schedule = TAX_SCHEDULE.get(state) if schedule is None: return np.nan rate = schedule[0][1] for yr, r in schedule: if yr <= year: rate = r else: break return rate def _load_state_gdp(): """Load state GDP by industry from FRED cache or compiled data.""" rows = [] for st in STATES: for sector, info in SECTOR_RHO.items(): fred_code = info['fred_code'] sid = f"{st}{fred_code}" rho = info['rho'] K = 1 - rho # curvature K ~ (1-rho) for large J try: s = load_fred(sid) if s is not None and len(s) > 0: for date, val in s.items(): rows.append({ 'state': st, 'sector': sector, 'year': date.year, 'quarter': date.quarter, 'gdp': val, 'rho': rho, 'K': K, }) except (FileNotFoundError, Exception): continue if not rows: path = COMPILED_DIR / 'laffer_state_gdp.csv' if path.exists(): return pd.read_csv(path) return None return pd.DataFrame(rows) # FRED IP series -> BEA I-O industry code mapping (18 manufacturing sectors) FRED_TO_BEA = { 'IPG311S': '311FT', # Food -> Food, beverage, tobacco 'IPG313A4S': '313TT', # Textile mills + products 'IPG315A6S': '315AL', # Apparel + leather 'IPG321S': '321', # Wood products 'IPG322S': '322', # Paper products 'IPG323S': '323', # Printing 'IPG324S': '324', # Petroleum and coal 'IPG325S': '325', # Chemical 'IPG326S': '326', # Plastics and rubber 'IPG327S': '327', # Nonmetallic mineral 'IPG331S': '331', # Primary metals 'IPG332S': '332', # Fabricated metals 'IPG333S': '333', # Machinery 'IPG334S': '334', # Computer and electronic 'IPG335S': '335', # Electrical equipment 'IPG336S': '336', # Transportation equipment (BEA: 3361MV + 3364OT) 'IPG337S': '337', # Furniture 'IPG339S': '339', # Miscellaneous manufacturing } # BEA I-O codes that combine into FRED 336 (transportation equipment) BEA_336_COMPONENTS = ['3361MV', '3364OT'] def _compute_input_hhi(): """Compute input concentration (HHI) from BEA I-O Use Table 2022. Returns dict mapping BEA industry code -> HHI (intermediate inputs only, excluding value-added components V001/V002/V003). """ io_path = COMPILED_DIR / 'bea_io_use_2022.csv' if not io_path.exists(): return None df = pd.read_csv(io_path) # Filter out totals and value-added rows totals = ['T001', 'T005', 'T006', 'T007', 'T008', 'T013', 'T014', 'T015', 'T016', 'T017', 'T018', 'T019', 'VABAS', 'VAPRO'] va_rows = ['V001', 'V002', 'V003'] flows = df[~df['ColCode'].isin(totals) & ~df['RowCode'].isin(totals) & ~df['RowCode'].isin(va_rows)].copy() flows['DataValue'] = pd.to_numeric(flows['DataValue'], errors='coerce') hhi_dict = {} bea_codes = set(FRED_TO_BEA.values()) for bea_code in bea_codes: if bea_code == '336': # Combine 3361MV + 3364OT ind_flows = flows[flows['ColCode'].isin(BEA_336_COMPONENTS)] ind_flows = ind_flows.groupby('RowCode')['DataValue'].sum().reset_index() else: ind_flows = flows[flows['ColCode'] == bea_code][['RowCode', 'DataValue']].copy() ind_flows = ind_flows.dropna(subset=['DataValue']) ind_flows = ind_flows[ind_flows['DataValue'] > 0] total = ind_flows['DataValue'].sum() if total <= 0: continue shares = ind_flows['DataValue'] / total hhi_dict[bea_code] = float((shares ** 2).sum()) return hhi_dict def _run_monetary_sensitivity_test(): """Test 3: Monetary Tightening Sensitivity by Input Concentration. Panel regression: IP_growth_it = alpha_i + gamma_t + beta1*dFFR_t + beta2*(dFFR_t * HHI_i) + eps_it CES prediction: beta2 < 0 (sectors with concentrated inputs lose more output when rates rise, because complementary production is more sensitive to friction — high HHI proxies high K). """ results = [] # 1. Compute input HHI from BEA I-O hhi_dict = _compute_input_hhi() if hhi_dict is None or len(hhi_dict) < 10: return results # 2. Load Fed Funds Rate (DFF, daily -> monthly average) try: ffr = load_fred('DFF') if ffr is None or len(ffr) < 100: return results except FileNotFoundError: return results ffr_monthly = ffr.resample('MS').mean() ffr_monthly.name = 'ffr' ffr_df = ffr_monthly.to_frame().reset_index() ffr_df.columns = ['date', 'ffr'] ffr_df['dffr'] = ffr_df['ffr'].diff() ffr_df = ffr_df.dropna(subset=['dffr']) # 3. Load IP subsector series and build panel rows = [] sector_names = {} bea_descriptions = { '311FT': 'Food/Bev/Tobacco', '313TT': 'Textiles', '315AL': 'Apparel', '321': 'Wood', '322': 'Paper', '323': 'Printing', '324': 'Petroleum/Coal', '325': 'Chemicals', '326': 'Plastics/Rubber', '327': 'Nonmetal Minerals', '331': 'Primary Metals', '332': 'Fab Metals', '333': 'Machinery', '334': 'Computers/Electronics', '335': 'Elec Equipment', '336': 'Transport Equipment', '337': 'Furniture', '339': 'Misc Mfg', } for fred_code, bea_code in FRED_TO_BEA.items(): if bea_code not in hhi_dict: continue try: s = load_fred(fred_code) if s is None or len(s) < 100: continue except FileNotFoundError: continue hhi = hhi_dict[bea_code] sector_names[bea_code] = bea_descriptions.get(bea_code, bea_code) ip_df = s.to_frame('ip').reset_index() ip_df.columns = ['date', 'ip'] ip_df['date'] = pd.to_datetime(ip_df['date']) ip_df = ip_df.set_index('date').resample('MS').mean().reset_index() ip_df['ip_growth'] = np.log(ip_df['ip']).diff() ip_df = ip_df.dropna(subset=['ip_growth']) for _, row in ip_df.iterrows(): rows.append({ 'date': row['date'], 'sector': bea_code, 'ip_growth': row['ip_growth'], 'hhi': hhi, }) if not rows: return results panel = pd.DataFrame(rows) panel = panel.merge(ffr_df[['date', 'dffr']], on='date', how='inner') # Also add 5-month lagged rate change (monetary transmission lag) ffr_df['dffr_L5'] = ffr_df['dffr'].shift(5) panel = panel.merge(ffr_df[['date', 'dffr_L5']], on='date', how='inner') panel['dffr_x_hhi'] = panel['dffr'] * panel['hhi'] panel['dffr_L5_x_hhi'] = panel['dffr_L5'] * panel['hhi'] panel['year'] = panel['date'].dt.year if len(panel) < 500 or not HAS_SM: return results try: from scipy.stats import linregress, spearmanr # ---- Sector FE, cluster-robust SEs by sector ---- # HHI is sector-level, so cluster at sector level. # 18 clusters is small — results are conservative check alongside # the cross-sectional Spearman. sec_dum = pd.get_dummies(panel['sector'], prefix='sec', drop_first=True, dtype=float) # Contemporaneous specification p_clean = panel.dropna(subset=['dffr', 'dffr_x_hhi', 'ip_growth']) X = pd.concat([ p_clean[['dffr', 'dffr_x_hhi']].reset_index(drop=True), pd.get_dummies(p_clean['sector'], prefix='sec', drop_first=True, dtype=float).reset_index(drop=True), ], axis=1) X = sm.add_constant(X) y = p_clean['ip_growth'].reset_index(drop=True) sectors_arr = p_clean['sector'].values res = sm.OLS(y.astype(float), X.astype(float)).fit( cov_type='cluster', cov_kwds={'groups': sectors_arr}) beta_dffr = res.params['dffr'] beta_int = res.params['dffr_x_hhi'] p_int = res.pvalues['dffr_x_hhi'] p_dffr = res.pvalues['dffr'] # Lag-5 specification (monetary transmission delay) p5 = panel.dropna(subset=['dffr_L5', 'dffr_L5_x_hhi', 'ip_growth']) X5 = pd.concat([ p5[['dffr_L5', 'dffr_L5_x_hhi']].reset_index(drop=True), pd.get_dummies(p5['sector'], prefix='sec', drop_first=True, dtype=float).reset_index(drop=True), ], axis=1) X5 = sm.add_constant(X5) y5 = p5['ip_growth'].reset_index(drop=True) res5 = sm.OLS(y5.astype(float), X5.astype(float)).fit( cov_type='cluster', cov_kwds={'groups': p5['sector'].values}) beta_int_L5 = res5.params['dffr_L5_x_hhi'] p_int_L5 = res5.pvalues['dffr_L5_x_hhi'] # Cross-sectional Spearman: HHI vs sector-level dFFR slope sector_slopes = {} sector_slopes_L5 = {} for bea_code in sorted(panel['sector'].unique()): sub = panel[panel['sector'] == bea_code].dropna(subset=['dffr']) if len(sub) > 50: sl = linregress(sub['dffr'].values, sub['ip_growth'].values) sector_slopes[bea_code] = { 'name': sector_names.get(bea_code, bea_code), 'hhi': float(hhi_dict.get(bea_code, np.nan)), 'slope': float(sl.slope), 'p': float(sl.pvalue), 'n': len(sub), } sub5 = panel[panel['sector'] == bea_code].dropna(subset=['dffr_L5']) if len(sub5) > 50: sl5 = linregress(sub5['dffr_L5'].values, sub5['ip_growth'].values) sector_slopes_L5[bea_code] = { 'name': sector_names.get(bea_code, bea_code), 'hhi': float(hhi_dict.get(bea_code, np.nan)), 'slope': float(sl5.slope), 'p': float(sl5.pvalue), 'n': len(sub5), } if len(sector_slopes) >= 5: hhis = [v['hhi'] for v in sector_slopes.values()] slopes = [v['slope'] for v in sector_slopes.values()] r_spear, p_spear = spearmanr(hhis, slopes) else: r_spear, p_spear = np.nan, np.nan if len(sector_slopes_L5) >= 5: hhis5 = [v['hhi'] for v in sector_slopes_L5.values()] slopes5 = [v['slope'] for v in sector_slopes_L5.values()] r_spear_L5, p_spear_L5 = spearmanr(hhis5, slopes5) else: r_spear_L5, p_spear_L5 = np.nan, np.nan # Verdict: use lag-5 Spearman (cross-sectional, avoids small-cluster # issue) as primary, contemporaneous panel as secondary # CES prediction: negative (high HHI = more sensitive) sign_ok = r_spear_L5 < 0 if not np.isnan(r_spear_L5) else beta_int < 0 p_primary = p_spear_L5 if not np.isnan(p_spear_L5) else p_int verd = verdict(p_primary, sign_ok, threshold=0.10) _cache['monetary_test'] = { 'beta_dffr': float(beta_dffr), 'p_dffr': float(p_dffr), 'beta_int': float(beta_int), 'p_int': float(p_int), 'beta_int_L5': float(beta_int_L5), 'p_int_L5': float(p_int_L5), 'N': int(res.nobs), 'R2': float(res.rsquared), 'n_sectors': len(panel['sector'].unique()), 'hhi_range': (float(panel['hhi'].min()), float(panel['hhi'].max())), 'sector_slopes': sector_slopes, 'sector_slopes_L5': sector_slopes_L5, 'spearman_r': float(r_spear), 'spearman_p': float(p_spear), 'spearman_r_L5': float(r_spear_L5), 'spearman_p_L5': float(p_spear_L5), 'date_range': (str(panel['date'].min().date()), str(panel['date'].max().date())), } n_sec = len(panel['sector'].unique()) hhi_lo = panel['hhi'].min() hhi_hi = panel['hhi'].max() results.append(log_result( 'laffer_curve_monetary_sensitivity', r_spear_L5, p_spear_L5, verd, details={ 'N': int(res.nobs), 'R2': float(res.rsquared), 'beta_dffr_x_hhi': float(beta_int), 'p_dffr_x_hhi_cluster': float(p_int), 'beta_dffr_L5_x_hhi': float(beta_int_L5), 'p_dffr_L5_x_hhi_cluster': float(p_int_L5), 'n_sectors': n_sec, 'hhi_range': f'{hhi_lo:.3f} to {hhi_hi:.3f}', 'spearman_L0': f'r={r_spear:.3f}, p={p_spear:.3f}', 'spearman_L5': f'r={r_spear_L5:.3f}, p={p_spear_L5:.3f}', }, description=( "Monetary tightening sensitivity by input concentration: tests whether " "sectors with concentrated inputs (high HHI from BEA I-O Use Table 2022) " "suffer larger output losses from Fed rate hikes, as CES curvature predicts. " "Two-step approach: (1) estimate each sector's dFFR→IP_growth slope, " "(2) Spearman correlation of HHI vs slope at 5-month lag (monetary " "transmission delay). Panel confirmation: IP_growth ~ dFFR + dFFR*HHI " f"with sector FE and sector-clustered SEs. {n_sec} FRED manufacturing IP " f"subsectors, monthly ({int(res.nobs)} obs, HHI {hhi_lo:.3f}-{hhi_hi:.3f}). " f"DFF (daily Fed Funds Rate → monthly avg) as friction. " f"Prediction: negative correlation (high HHI = more sensitive). " f"Lag-5 Spearman r={r_spear_L5:.3f}, p={p_spear_L5:.3f}; " f"panel beta(dFFR×HHI)={beta_int:.4f}, cluster p={p_int:.4f}." ), )) except Exception: pass return results _cache = {} def main(): results = [] if not HAS_SM: return [log_result('laffer_curve_tax_sensitivity', np.nan, np.nan, 'SKIPPED', details={'reason': 'statsmodels not installed'}, description="Tests CES curvature prediction for tax sensitivity.", )] gdp_df = _load_state_gdp() if gdp_df is None or len(gdp_df) == 0: return [log_result('laffer_curve_tax_sensitivity', np.nan, np.nan, 'SKIPPED', details={'reason': 'State GDP data not found'}, description="Tests CES curvature prediction for tax sensitivity.", )] # Annualize if quarterly if 'quarter' in gdp_df.columns: gdp_ann = gdp_df.groupby(['state', 'sector', 'year', 'rho', 'K'] )['gdp'].mean().reset_index() else: gdp_ann = gdp_df.copy() # Add tax rates gdp_ann['tax_rate'] = gdp_ann.apply( lambda r: get_tax_rate(r['state'], r['year']), axis=1) gdp_ann = gdp_ann.dropna(subset=['tax_rate', 'gdp']) gdp_ann = gdp_ann[gdp_ann['gdp'] > 0] if len(gdp_ann) < 100: return [log_result('laffer_curve_tax_sensitivity', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(gdp_ann)} obs'}, description="Tests CES curvature prediction for tax sensitivity.", )] gdp_ann['log_gdp'] = np.log(gdp_ann['gdp']) # ================================================================== # Test 1: Panel interaction — GDP growth ~ dtax * K + state FE + year FE # CES prediction: high K (more complementary) → more sensitive to tax friction # So dtax * K interaction should be negative (tax increase * high K = more damage) # ================================================================== gdp_ann_sorted = gdp_ann.sort_values(['state', 'sector', 'year']) gdp_ann_sorted['gdp_growth'] = gdp_ann_sorted.groupby( ['state', 'sector'])['log_gdp'].diff() gdp_ann_sorted['dtax'] = gdp_ann_sorted.groupby( ['state', 'sector'])['tax_rate'].diff() panel = gdp_ann_sorted.dropna(subset=['gdp_growth', 'dtax']).copy() panel['dtax_x_K'] = panel['dtax'] * panel['K'] if len(panel) >= 100: try: st_dum = pd.get_dummies(panel['state'], prefix='st', drop_first=True, dtype=float) yr_dum = pd.get_dummies(panel['year'], prefix='yr', drop_first=True, dtype=float) X = pd.concat([panel[['dtax', 'dtax_x_K']].reset_index(drop=True), st_dum.reset_index(drop=True), yr_dum.reset_index(drop=True)], axis=1) X = sm.add_constant(X) y = panel['gdp_growth'].reset_index(drop=True) res = sm.OLS(y.astype(float), X.astype(float)).fit(cov_type='HC1') beta_int = res.params['dtax_x_K'] p_int = res.pvalues['dtax_x_K'] beta_dtax = res.params['dtax'] n_nz = int((panel['dtax'] != 0).sum()) # CES prediction: interaction should be negative # (high K + tax increase = more output loss) sign_ok = beta_int < 0 verd = verdict(p_int, sign_ok, threshold=0.10) _cache['panel_test'] = { 'beta_int': beta_int, 'p_int': p_int, 'beta_dtax': beta_dtax, 'N': int(res.nobs), 'R2': res.rsquared, 'n_nonzero_dtax': n_nz, } results.append(log_result( 'laffer_curve_tax_sensitivity', beta_int, p_int, verd, details={'N': int(res.nobs), 'R2': res.rsquared, 'beta_dtax': beta_dtax, 'beta_dtax_x_K': beta_int, 'n_nonzero_dtax': n_nz, 'rho_source': 'Oberfield-Raval 2021 / Atalay 2017'}, description=( "CES curvature tax sensitivity: GDP growth ~ delta_tax + delta_tax * K " "with state and year FE (HC1 robust SEs). Tests whether high-K " "(complementary) sectors suffer larger output losses from tax changes, " "as predicted by CES potential (Paper 2, line 795: 'marginal deadweight " "loss is larger where the Laffer decline is steeper'). " "K = 1 - rho from Oberfield-Raval (2021) / Atalay (2017): Manufacturing " "K=1.30, Information K=1.19, others K=1.07. FRED state GDP by industry " f"(50 states, 6 sectors, {int(res.nobs)} obs, {n_nz} non-zero dtax). " "Prediction: beta(dtax*K) < 0 (high K = more sensitive to tax friction)." ), )) except Exception: pass # ================================================================== # Test 2: TCJA SALT cap event study (2018) # The 2017 TCJA capped SALT deductions at $10K, effectively raising # tax burden for high-tax states. Tests whether the SALT shock's # effect on sector GDP varies by curvature K. # ================================================================== # Treatment: states with pre-TCJA tax rate >= 7% (top quartile) # Control: states with pre-TCJA tax rate <= 3% (bottom quartile, incl. no-income-tax) # Pre: 2015-2017, Post: 2018-2019 (exclude 2020 COVID) pre_years = [2015, 2016, 2017] post_years = [2018, 2019] # Pre-TCJA tax rates (as of 2017) state_tax_2017 = {st: get_tax_rate(st, 2017) for st in STATES} high_tax_states = {st for st, rate in state_tax_2017.items() if not np.isnan(rate) and rate >= 7.0} low_tax_states = {st for st, rate in state_tax_2017.items() if not np.isnan(rate) and rate <= 3.0} tcja_data = gdp_ann[ (gdp_ann['year'].isin(pre_years + post_years)) & (gdp_ann['state'].isin(high_tax_states | low_tax_states)) ].copy() if len(tcja_data) >= 50: tcja_data['post'] = (tcja_data['year'] >= 2018).astype(float) tcja_data['treated'] = tcja_data['state'].isin(high_tax_states).astype(float) tcja_data['did'] = tcja_data['post'] * tcja_data['treated'] tcja_data['did_x_K'] = tcja_data['did'] * tcja_data['K'] try: # DID with sector curvature interaction: # log_gdp ~ treated + post + did + did*K + sector FE + state FE sec_dum = pd.get_dummies(tcja_data['sector'], prefix='sec', drop_first=True, dtype=float) st_dum2 = pd.get_dummies(tcja_data['state'], prefix='st', drop_first=True, dtype=float) yr_dum2 = pd.get_dummies(tcja_data['year'], prefix='yr', drop_first=True, dtype=float) X2 = pd.concat([ tcja_data[['treated', 'post', 'did', 'did_x_K']].reset_index(drop=True), sec_dum.reset_index(drop=True), st_dum2.reset_index(drop=True), yr_dum2.reset_index(drop=True), ], axis=1) X2 = sm.add_constant(X2) y2 = tcja_data['log_gdp'].reset_index(drop=True) res2 = sm.OLS(y2.astype(float), X2.astype(float)).fit(cov_type='HC1') beta_did = res2.params['did'] p_did = res2.pvalues['did'] beta_did_K = res2.params['did_x_K'] p_did_K = res2.pvalues['did_x_K'] # CES prediction: did_x_K should be negative # (SALT shock hurts high-K sectors more in treated states) sign_ok_tcja = beta_did_K < 0 verd_tcja = verdict(p_did_K, sign_ok_tcja, threshold=0.10) n_treated = int(tcja_data['treated'].sum()) n_control = int((1 - tcja_data['treated']).sum()) _cache['tcja_test'] = { 'beta_did': beta_did, 'p_did': p_did, 'beta_did_K': beta_did_K, 'p_did_K': p_did_K, 'N': int(res2.nobs), 'R2': res2.rsquared, 'n_treated': n_treated, 'n_control': n_control, 'high_tax_states': sorted(high_tax_states), 'low_tax_states': sorted(low_tax_states), } results.append(log_result( 'laffer_curve_tcja_salt', beta_did_K, p_did_K, verd_tcja, details={ 'N': int(res2.nobs), 'R2': res2.rsquared, 'beta_did': beta_did, 'p_did': p_did, 'beta_did_K': beta_did_K, 'p_did_K': p_did_K, 'n_high_tax_states': len(high_tax_states), 'n_low_tax_states': len(low_tax_states), }, description=( "TCJA SALT cap event study: the 2017 TCJA capped SALT deductions at $10K, " "effectively raising tax burden for high-tax states. DID: log(GDP) ~ " "treated * post * K with state, sector, and year FE (HC1 robust SEs). " f"Treatment: {len(high_tax_states)} states with 2017 rate >= 7%. " f"Control: {len(low_tax_states)} states with rate <= 3%. " "Pre: 2015-2017, Post: 2018-2019 (excluding COVID 2020). " "CES prediction: DID*K < 0 (high-K sectors lose more in treated states " "because curvature amplifies misallocation from tax friction)." ), )) except Exception: pass # ================================================================== # Test 3: Monetary Tightening Sensitivity by Input Concentration # Uses Fed Funds Rate changes as universal friction + BEA I-O input # HHI as K proxy + FRED IP subsectors as output. Solves proxy failure # in Tests 1-2: monetary policy affects all sectors via cost of capital. # ================================================================== monetary_results = _run_monetary_sensitivity_test() if monetary_results: results.extend(monetary_results) if not results: results.append(log_result('laffer_curve_tax_sensitivity', np.nan, np.nan, 'SKIPPED', description="Fallback when no subtests could run.", )) _cache['gdp_ann'] = gdp_ann _cache['panel'] = panel if 'panel' in dir() else None _cache['results'] = results return results def figures(): """Generate Laffer curve dashboard: panel sensitivity and TCJA event study.""" if not _cache: main() paths = [] gdp_ann = _cache.get('gdp_ann') if gdp_ann is None or len(gdp_ann) == 0: return paths import matplotlib.pyplot as plt # Determine layout based on which tests ran has_panel = 'panel_test' in _cache has_tcja = 'tcja_test' in _cache has_monetary = 'monetary_test' in _cache n_panels = has_panel + has_tcja + has_monetary if n_panels == 0: return paths fig, axes = setup_subplots(1, max(n_panels, 1), width=7 * max(n_panels, 1), height=5) if n_panels == 1: axes = [axes] ax_idx = 0 # --- Panel interaction: tax change sensitivity by sector --- if has_panel: ax = axes[ax_idx] ax_idx += 1 gdp_sorted = gdp_ann.sort_values(['state', 'sector', 'year']).copy() gdp_sorted['gdp_growth'] = gdp_sorted.groupby(['state', 'sector'])['log_gdp'].diff() gdp_sorted['dtax'] = gdp_sorted.groupby(['state', 'sector'])['tax_rate'].diff() nz = gdp_sorted.dropna(subset=['gdp_growth', 'dtax']) nz = nz[nz['dtax'] != 0] # Color by K value k_vals = sorted(nz['K'].unique()) cmap = plt.cm.viridis for K_val in k_vals: sub = nz[nz['K'] == K_val] sector_name = sub['sector'].iloc[0] if len(sub) > 0 else '?' color = cmap((K_val - min(k_vals)) / (max(k_vals) - min(k_vals) + 0.01)) ax.scatter(sub['dtax'], sub['gdp_growth'], s=15, alpha=0.4, color=color, label=f'{sector_name} (K={K_val:.2f})') if len(sub) > 5: slope, intercept, _, _, _ = stats.linregress(sub['dtax'], sub['gdp_growth']) x_fit = np.linspace(sub['dtax'].min(), sub['dtax'].max(), 50) ax.plot(x_fit, intercept + slope * x_fit, color=color, alpha=0.7, lw=2) pt = _cache.get('panel_test', {}) beta = pt.get('beta_int', np.nan) p_int = pt.get('p_int', np.nan) ax.axhline(0, color='black', lw=0.5) ax.axvline(0, color='black', lw=0.5) ax.set_xlabel('Tax rate change (pp)') ax.set_ylabel('GDP growth (log diff)') ax.set_title(f'Panel: dtax*K interaction beta={beta:.4f}, p={p_int:.3f}') ax.legend(fontsize=7, loc='upper right') ax.grid(True, alpha=0.3) # --- TCJA SALT cap event study --- if has_tcja: ax = axes[ax_idx] ax_idx += 1 tcja = _cache.get('tcja_test', {}) high_tax = set(tcja.get('high_tax_states', [])) low_tax = set(tcja.get('low_tax_states', [])) pre_years = [2015, 2016, 2017] post_years = [2018, 2019] # Compute average GDP growth by treatment group and K tcja_data = gdp_ann[ (gdp_ann['year'].isin(pre_years + post_years)) & (gdp_ann['state'].isin(high_tax | low_tax)) ].copy() tcja_data['treated'] = tcja_data['state'].isin(high_tax) tcja_data['period'] = tcja_data['year'].apply( lambda y: 'Pre (2015-17)' if y <= 2017 else 'Post (2018-19)') # Bar chart: DID by sector sectors = sorted(tcja_data['sector'].unique()) did_by_sector = [] for sector in sectors: s = tcja_data[tcja_data['sector'] == sector] K_val = s['K'].iloc[0] treat_pre = s[(s['treated']) & (s['year'] <= 2017)]['log_gdp'].mean() treat_post = s[(s['treated']) & (s['year'] >= 2018)]['log_gdp'].mean() ctrl_pre = s[(~s['treated']) & (s['year'] <= 2017)]['log_gdp'].mean() ctrl_post = s[(~s['treated']) & (s['year'] >= 2018)]['log_gdp'].mean() did_est = (treat_post - treat_pre) - (ctrl_post - ctrl_pre) did_by_sector.append({'sector': sector, 'K': K_val, 'did': did_est}) did_df = pd.DataFrame(did_by_sector).sort_values('K') colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(did_df))) bars = ax.bar(range(len(did_df)), did_df['did'], color=colors, edgecolor='black', lw=0.5) ax.set_xticks(range(len(did_df))) ax.set_xticklabels([f"{r['sector']}\nK={r['K']:.2f}" for _, r in did_df.iterrows()], fontsize=7, rotation=30, ha='right') ax.axhline(0, color='black', lw=0.5) ax.set_ylabel('DID estimate (log GDP)') beta_did_K = tcja.get('beta_did_K', np.nan) p_did_K = tcja.get('p_did_K', np.nan) ax.set_title(f'TCJA SALT Cap DID by Sector\n' f'DID*K beta={beta_did_K:.4f}, p={p_did_K:.3f}') ax.grid(True, alpha=0.3, axis='y') # --- Monetary tightening: HHI vs sector FFR sensitivity (lag 5) --- if has_monetary: ax = axes[ax_idx] ax_idx += 1 mt = _cache['monetary_test'] # Use lag-5 slopes as primary (monetary transmission delay) sector_slopes = mt.get('sector_slopes_L5', mt.get('sector_slopes', {})) if sector_slopes: hhis = [v['hhi'] for v in sector_slopes.values()] slopes = [v['slope'] for v in sector_slopes.values()] names = [v['name'] for v in sector_slopes.values()] ax.scatter(hhis, slopes, s=40, alpha=0.7, c='steelblue', edgecolors='black', lw=0.5) for i, name in enumerate(names): ax.annotate(name, (hhis[i], slopes[i]), fontsize=5.5, xytext=(3, 3), textcoords='offset points') if len(hhis) >= 3: sl, intercept, _, _, _ = stats.linregress(hhis, slopes) x_fit = np.linspace(min(hhis), max(hhis), 50) ax.plot(x_fit, intercept + sl * x_fit, 'r-', lw=1.5, alpha=0.7) ax.axhline(0, color='black', lw=0.5) ax.set_xlabel('Input Concentration (HHI)') ax.set_ylabel('dFFR(t-5) → IP growth slope') r_sp = mt.get('spearman_r_L5', np.nan) p_sp = mt.get('spearman_p_L5', np.nan) ax.set_title(f'Monetary Sensitivity (lag 5m)\n' f'Spearman r={r_sp:.3f}, p={p_sp:.3f}') ax.grid(True, alpha=0.3) fig.tight_layout() # Save combined figure for all subtests paths = save_figure(fig, 'laffer_curve_tax_sensitivity', subfolder='potential') # Also save a separate figure for the monetary test if it exists if has_monetary: fig2, ax2 = setup_figure() mt = _cache['monetary_test'] sector_slopes = mt.get('sector_slopes_L5', mt.get('sector_slopes', {})) if sector_slopes: hhis = [v['hhi'] for v in sector_slopes.values()] slopes = [v['slope'] for v in sector_slopes.values()] names = [v['name'] for v in sector_slopes.values()] ax2.scatter(hhis, slopes, s=50, alpha=0.7, c='steelblue', edgecolors='black', lw=0.5) for i, name in enumerate(names): ax2.annotate(name, (hhis[i], slopes[i]), fontsize=6.5, xytext=(4, 4), textcoords='offset points') if len(hhis) >= 3: sl, intercept, _, _, _ = stats.linregress(hhis, slopes) x_fit = np.linspace(min(hhis), max(hhis), 50) ax2.plot(x_fit, intercept + sl * x_fit, 'r-', lw=1.5, alpha=0.7) ax2.axhline(0, color='black', lw=0.5) ax2.set_xlabel('Input Concentration (HHI from BEA I-O)') ax2.set_ylabel('FFR Change (t-5) → IP Growth Slope') r_sp = mt.get('spearman_r_L5', mt.get('spearman_r', np.nan)) p_sp = mt.get('spearman_p_L5', mt.get('spearman_p', np.nan)) beta = mt.get('beta_int', np.nan) p_int = mt.get('p_int', np.nan) n_sec = mt.get('n_sectors', 0) N = mt.get('N', 0) ax2.set_title(f'Monetary Tightening Sensitivity by Input Concentration\n' f'Spearman r={r_sp:.3f}, p={p_sp:.3f} (lag 5m) | ' f'Panel β={beta:.4f}, cluster p={p_int:.4f}\n' f'{n_sec} sectors, {N} monthly obs') ax2.grid(True, alpha=0.3) fig2.tight_layout() paths.extend(save_figure(fig2, 'laffer_curve_monetary_sensitivity', subfolder='potential')) 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']})")

N	5981
R2	0.1636
p value	0.6782
statistic	0.0043

Laffer Curve Tax Sensitivity

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources