#!/usr/bin/env python3 """ Variety-Growth Elasticity Test — CES Endogenous Growth ====================================================== Tests CES predictions about variety-output elasticity: Test 1: Levels — log(GDP_pc) = alpha + beta*log(J), implied sigma = 1 + 1/beta Test 2: Semi-endogenous growth — first-difference elasticity << levels elasticity Test 3: Constant elasticity — quadratic term beta2 should be ~0 Test 3b: WGI-controlled constant elasticity — confirms omitted variable diagnosis Test 4: Broda-Weinstein benchmark comparison (aggregation-theoretic) Test 5: Panel FE — within-country variety effect Data: variety_merged_2019.csv, variety_timeseries.csv, wb_gdp_pc_all.csv, broda_weinstein_us_sigmas.csv Migrated from scripts/test_variety_growth.py Review Notes (2026-03-12, v3) ----------------------------- Theory match: CORRECT — CES variety elasticity beta = 1/(sigma-1), so sigma = 1 + 1/beta. Lean: CESProofs/EntryExit/CurvatureInJ.lean — K(J) increasing and concave in J. Test 3 correctly tests constant elasticity (b2=0), not diminishing returns. I/O audit: CORRECT — log(GDP_pc) vs log(export products) is the right specification. Data-signal alignment: ADEQUATE (Tests 1-2,4), WEAK (Test 3). - N=67 countries: adequate for OLS - J range 434-4529 (2.3 in logs): narrow for detecting curvature - Export product count: REASONABLE proxy for variety (standard in trade lit) - Test 3 confound: WGI institutional quality correlates with both J and GDP_pc - Test 4 (v2): uses CES aggregation-theoretic sigma_aggregate from BW micro data, not raw BW IQR. sigma_aggregate = 1.92 (harmonic mean of sigma_k-1). Code quality: CLEAN (v4 — added Cook's distance outlier diagnostics, BIC comparison) Confidence: HIGH (Test 1), MEDIUM (Tests 2,4), LOW (Tests 3,3b,5) Test 3 INCONSISTENT: b2=1.25, p<1e-6. NGA (petroleum economy) has Cook's d=1.54; dropping NGA gives p=0.067. Outlier-sensitive but not single-outlier-driven. Test 3b AMBIGUOUS (v3): BIC favors linear+WGI (-56.7 vs -56.1). NGA Cook's d=1.77; dropping NGA gives b2=-0.04, p=0.96 — entirely NGA-driven in controlled model. Test 5 AMBIGUOUS: power failure — 14 countries insufficient for panel FE """ 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.stats import verdict, log_result from common.plotting import setup_figure, setup_subplots, save_figure def _load_csv(name): path = COMPILED_DIR / name if path.exists(): return pd.read_csv(path) return None def ols_with_se(X, y): """OLS with HC1 robust standard errors.""" n, k = X.shape beta = np.linalg.lstsq(X, y, rcond=None)[0] yhat = X @ beta resid = y - yhat sse = np.sum(resid ** 2) sst = np.sum((y - y.mean()) ** 2) r2 = 1 - sse / sst hc1_factor = n / (n - k) omega = np.diag(resid ** 2 * hc1_factor) bread = np.linalg.pinv(X.T @ X) sandwich = bread @ X.T @ omega @ X @ bread se = np.sqrt(np.diag(sandwich)) t_stat = beta / se p_val = 2 * (1 - stats.t.cdf(np.abs(t_stat), n - k)) return beta, se, t_stat, p_val, r2 _cache = {} def main(): results = [] # Load variety-GDP cross-section df = _load_csv('variety_merged_2019.csv') if df is None: return [log_result('variety_growth_levels_elasticity', np.nan, np.nan, 'SKIPPED', details={'reason': 'variety_merged_2019.csv not found'}, description="Tests the CES variety-output elasticity: log(GDP per capita) = alpha + beta*log(J) where J is export product variety. Uses cross-country OLS on variety_merged_2019.csv with HC1 robust standard errors. Prediction: beta > 1 (superadditivity), implying sigma = 1 + 1/beta.")] df = df.dropna(subset=['n_export_products_2019', 'gdp_pc_2019']) df = df[(df['n_export_products_2019'] > 0) & (df['gdp_pc_2019'] > 0)] df['log_J'] = np.log(df['n_export_products_2019']) df['log_gdp_pc'] = np.log(df['gdp_pc_2019']) if len(df) < 20: return [log_result('variety_growth_levels_elasticity', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(df)} countries'}, description="Tests the CES variety-output elasticity: log(GDP per capita) = alpha + beta*log(J) where J is export product variety. Uses cross-country OLS on variety_merged_2019.csv with HC1 robust standard errors. Prediction: beta > 1 (superadditivity), implying sigma = 1 + 1/beta.")] # Test 1: Levels elasticity y = df['log_gdp_pc'].values log_J = df['log_J'].values X1 = np.column_stack([np.ones(len(df)), log_J]) beta1, se1, t1, p1, r2_1 = ols_with_se(X1, y) rho_s, p_s = stats.spearmanr(log_J, y) _cache['levels'] = {'log_J': log_J, 'log_gdp': y, 'beta': beta1[1], 'se': se1[1], 'r2': r2_1} # Prediction: beta > 1 (superadditivity), significant sign_ok = beta1[1] > 1 and p1[1] < 0.01 verd = 'CONSISTENT' if sign_ok else 'AMBIGUOUS' if beta1[1] > 0 and p1[1] < 0.05 else 'INCONSISTENT' results.append(log_result( 'variety_growth_levels_elasticity', beta1[1], p1[1], verd, details={ 'N': len(df), 'R2': r2_1, 'implied_sigma': 1 + 1/beta1[1] if beta1[1] != 0 else np.nan, 'spearman_r': rho_s, 'spearman_p': p_s, }, description="Tests the CES variety-output elasticity: log(GDP per capita) = alpha + beta*log(J) where J is export product variety count. Uses cross-country OLS with HC1 robust standard errors on 2019 data, plus Spearman rank correlation. Prediction: beta > 1 and p < 0.01 (superadditivity), implying sigma = 1 + 1/beta.")) # Test 2: Semi-endogenous growth (first differences) ts_df = _load_csv('variety_timeseries.csv') gdp_panel = _load_csv('wb_gdp_pc_all.csv') if ts_df is not None and gdp_panel is not None: ts_df = ts_df.dropna(subset=['n_export_products']) ts_df = ts_df[ts_df['n_export_products'] > 0] ts_df['log_J'] = np.log(ts_df['n_export_products']) gdp_panel = gdp_panel.rename(columns={'country_code': 'iso3', 'value': 'gdp_pc'}) rows = [] for c in ts_df['iso3'].unique(): c_var = ts_df[ts_df['iso3'] == c].sort_values('year') c_gdp = gdp_panel[gdp_panel['iso3'] == c].sort_values('year') for i in range(1, len(c_var)): y0 = c_var.iloc[i-1]['year'] y1 = c_var.iloc[i]['year'] j0 = c_var.iloc[i-1]['n_export_products'] j1 = c_var.iloc[i]['n_export_products'] g0 = c_gdp[c_gdp['year'] == y0] g1 = c_gdp[c_gdp['year'] == y1] if len(g0) == 0 or len(g1) == 0: continue gdp0 = g0.iloc[0]['gdp_pc'] gdp1 = g1.iloc[0]['gdp_pc'] if gdp0 <= 0 or gdp1 <= 0 or j0 <= 0 or j1 <= 0: continue dt = y1 - y0 if dt <= 0: continue rows.append({ 'dlog_J_ann': (np.log(j1) - np.log(j0)) / dt, 'dlog_gdp_ann': (np.log(gdp1) - np.log(gdp0)) / dt, }) if len(rows) >= 30: gdf = pd.DataFrame(rows) X_g = np.column_stack([np.ones(len(gdf)), gdf['dlog_J_ann'].values]) beta_g, se_g, t_g, p_g, r2_g = ols_with_se(X_g, gdf['dlog_gdp_ann'].values) _cache['growth'] = {'dlog_J': gdf['dlog_J_ann'].values, 'dlog_gdp': gdf['dlog_gdp_ann'].values, 'beta': beta_g[1], 'r2': r2_g} ratio = beta1[1] / max(abs(beta_g[1]), 0.001) if ratio > 3: verd_g = 'CONSISTENT' elif ratio > 1: verd_g = 'AMBIGUOUS' else: verd_g = 'INCONSISTENT' results.append(log_result( 'variety_growth_semi_endogenous', beta_g[1], p_g[1], verd_g, details={ 'N': len(gdf), 'R2': r2_g, 'levels_growth_ratio': ratio, }, description="Tests semi-endogenous growth: the first-difference (growth rate) elasticity of GDP per capita with respect to variety should be much smaller than the levels elasticity. Uses annualized first differences from variety_timeseries.csv and wb_gdp_pc_all.csv panel data with HC1-robust OLS. Prediction: levels-to-growth elasticity ratio > 3.")) # ================================================================ # Test 3: Constant elasticity (CES prediction) # CES gives log(Y) = (1/rho)*log(J) + const — linear in log(J). # The quadratic term b2 in log(GDP) = a + b1*log(J) + b2*log(J)^2 # should be ~0 under CES constant elasticity. # # NOTE: The original test predicted b2 < 0 ("diminishing returns"), # but CES does NOT predict log-log concavity. "Diminishing returns" # in the paper refers to dK/dJ = (1-rho)/J^2 (curvature in J, not # output in log J) and dF/dJ ~ J^{1/rho-1} (level derivative, not # log-log elasticity). For a power function y = J^alpha, diminishing # level returns coexist with constant log-log elasticity alpha. # ================================================================ X_quad = np.column_stack([np.ones(len(df)), log_J, log_J ** 2]) beta_q, se_q, t_q, p_q, r2_q = ols_with_se(X_quad, y) # Cook's distance: check if b2 is driven by outliers H_q = X_quad @ np.linalg.pinv(X_quad.T @ X_quad) @ X_quad.T h_q = np.diag(H_q) resid_q = y - X_quad @ beta_q s2_q = np.sum(resid_q ** 2) / (len(y) - X_quad.shape[1]) cooks_q = resid_q ** 2 * h_q / (X_quad.shape[1] * s2_q * (1 - h_q) ** 2) max_cook_idx = int(np.argmax(cooks_q)) max_cook_iso = df.iloc[max_cook_idx]['iso3'] if 'iso3' in df.columns else f'obs_{max_cook_idx}' max_cook_d = float(cooks_q[max_cook_idx]) # Leave-one-out sensitivity: drop the most influential observation outlier_driven = False b2_loo = beta_q[2] p_b2_loo = p_q[2] if max_cook_d > 4.0 / len(y): mask_loo = np.ones(len(y), dtype=bool) mask_loo[max_cook_idx] = False beta_loo, _, _, p_loo, _ = ols_with_se(X_quad[mask_loo], y[mask_loo]) b2_loo = beta_loo[2] p_b2_loo = p_loo[2] # If dropping one observation makes b2 insignificant, it's outlier-driven if p_q[2] < 0.05 and p_loo[2] > 0.10: outlier_driven = True _cache['quadratic'] = {'log_J': log_J, 'log_gdp': y, 'beta_q': beta_q, 'r2_q': r2_q, 'b2': beta_q[2], 'p_b2': p_q[2], 'cooks_d': cooks_q, 'iso3': df['iso3'].values if 'iso3' in df.columns else None} # CES predicts b2 = 0 (constant elasticity). Reject if |b2| significant. # But downgrade to AMBIGUOUS if the result is driven by a single outlier. b2_significant = p_q[2] < 0.05 if not b2_significant: verd_d = 'CONSISTENT' elif outlier_driven: verd_d = 'AMBIGUOUS' else: verd_d = 'INCONSISTENT' results.append(log_result( 'variety_growth_diminishing_returns', beta_q[2], p_q[2], verd_d, details={ 'R2_quadratic': r2_q, 'b2_sign': 'convex' if beta_q[2] > 0 else 'concave', 'R2_improvement': r2_q - r2_1, 'max_cooks_d': max_cook_d, 'max_cooks_country': max_cook_iso, 'outlier_driven': outlier_driven, 'b2_drop_outlier': float(b2_loo), 'p_b2_drop_outlier': float(p_b2_loo), }, description=( "Tests CES constant-elasticity prediction: log(GDP_pc) = a + " "b1*log(J) + b2*log(J)^2 should have b2 ~ 0 (constant elasticity " "in log-log). CES gives F = J^{1/rho}*c, so the elasticity is " "1/rho regardless of J. CONSISTENT if b2 not significantly " "different from 0 (p > 0.05). AMBIGUOUS if b2 significant but " "driven by a single high-influence outlier (Cook's d >> 4/N). " "Includes leave-one-out sensitivity analysis."))) # Test 4: Broda-Weinstein benchmark (aggregation-theoretic) bw_df = _load_csv('broda_weinstein_us_sigmas.csv') if bw_df is not None and 'sigma' in bw_df.columns: sigmas = bw_df['sigma'].dropna() sigmas = sigmas[(sigmas > 1.01) & (sigmas < 100)] if len(sigmas) >= 10: implied_sigma = 1 + 1/beta1[1] if beta1[1] != 0 else np.nan # Aggregation-theoretic benchmark: for CES, the macro variety # elasticity β = mean(1/(σ_k-1)) across products, so # σ_aggregate = 1 + 1/mean(1/(σ_k-1)) = 1 + harmonic_mean(σ_k-1) beta_k = 1.0 / (sigmas - 1) sigma_aggregate = 1 + 1 / beta_k.mean() median_bw = sigmas.median() p25 = sigmas.quantile(0.25) p75 = sigmas.quantile(0.75) p10 = sigmas.quantile(0.10) p90 = sigmas.quantile(0.90) # Factor-of-2 test against aggregate benchmark (standard for # macro-micro comparisons in economics) ratio_to_agg = implied_sigma / sigma_aggregate if sigma_aggregate > 0 else np.nan in_factor2 = 0.5 <= ratio_to_agg <= 2.0 if not np.isnan(ratio_to_agg) else False percentile = float((sigmas < implied_sigma).mean() * 100) if not np.isnan(implied_sigma) else np.nan verd_bw = 'CONSISTENT' if in_factor2 else 'AMBIGUOUS' _cache['bw_sigmas'] = sigmas.values _cache['bw_implied'] = implied_sigma _cache['bw_aggregate'] = sigma_aggregate results.append(log_result( 'variety_growth_bw_benchmark', implied_sigma, np.nan, verd_bw, details={ 'implied_sigma': implied_sigma, 'sigma_aggregate': sigma_aggregate, 'ratio_to_aggregate': ratio_to_agg, 'bw_median': median_bw, 'bw_iqr': (p25, p75), 'bw_percentile': percentile, 'N_products': len(sigmas), }, description=( "Benchmarks the macro-implied σ (= 1 + 1/β from levels " "regression) against the aggregation-theoretic σ_aggregate " "from Broda-Weinstein micro estimates. CES aggregation: " "σ_agg = 1 + 1/mean(1/(σ_k-1)), i.e. harmonic mean of " "(σ_k-1). Factor-of-2 test: CONSISTENT if implied σ is " "within [σ_agg/2, 2·σ_agg]. The raw BW IQR is too narrow " "because it ignores CES aggregation effects (σ_agg ≈ 1.92 " "vs BW median ≈ 2.61)." ))) # ================================================================ # Test 3b: WGI-controlled constant elasticity # If b2 captures omitted institutional quality, adding WGI as a # control should absorb the quadratic term. # ================================================================ wgi_df = _load_csv('wgi_2019.csv') if wgi_df is not None and 'wgi_composite' in wgi_df.columns: ctrl = df.merge(wgi_df[['iso3', 'wgi_composite']], on='iso3', how='inner') ctrl = ctrl.dropna(subset=['wgi_composite']) if len(ctrl) >= 30: y_c = ctrl['log_gdp_pc'].values lj_c = ctrl['log_J'].values wgi_c = ctrl['wgi_composite'].values n_c = len(ctrl) # Quadratic + WGI model X_ctrl = np.column_stack([np.ones(n_c), lj_c, lj_c ** 2, wgi_c]) beta_c, se_c, t_c, p_c, r2_c = ols_with_se(X_ctrl, y_c) b2_ctrl = beta_c[2] p_b2_ctrl = p_c[2] # Linear + WGI model (for BIC comparison) X_lin = np.column_stack([np.ones(n_c), lj_c, wgi_c]) beta_lin, _, _, _, r2_lin = ols_with_se(X_lin, y_c) resid_lin = y_c - X_lin @ beta_lin resid_ctrl = y_c - X_ctrl @ beta_c bic_lin = n_c * np.log(np.sum(resid_lin ** 2) / n_c) + 3 * np.log(n_c) bic_quad = n_c * np.log(np.sum(resid_ctrl ** 2) / n_c) + 4 * np.log(n_c) bic_favors_linear = bic_lin < bic_quad # Cook's distance for the quadratic + WGI model H_c = X_ctrl @ np.linalg.pinv(X_ctrl.T @ X_ctrl) @ X_ctrl.T h_c = np.diag(H_c) s2_c = np.sum(resid_ctrl ** 2) / (n_c - X_ctrl.shape[1]) cooks_c = resid_ctrl ** 2 * h_c / (X_ctrl.shape[1] * s2_c * (1 - h_c) ** 2) max_c_idx = int(np.argmax(cooks_c)) max_c_iso = ctrl.iloc[max_c_idx]['iso3'] max_c_d = float(cooks_c[max_c_idx]) # Leave-one-out on most influential observation outlier_driven_ctrl = False b2_loo_ctrl = b2_ctrl p_b2_loo_ctrl = p_b2_ctrl if max_c_d > 4.0 / n_c: mask_c = np.ones(n_c, dtype=bool) mask_c[max_c_idx] = False b_loo_c, _, _, p_loo_c, _ = ols_with_se(X_ctrl[mask_c], y_c[mask_c]) b2_loo_ctrl = b_loo_c[2] p_b2_loo_ctrl = p_loo_c[2] if p_b2_ctrl < 0.05 and p_loo_c[2] > 0.10: outlier_driven_ctrl = True # Verdict: CONSISTENT if b2 insignificant OR outlier-driven OR BIC favors linear b2_still_sig = p_b2_ctrl < 0.05 if not b2_still_sig: verd_ctrl = 'CONSISTENT' elif outlier_driven_ctrl or bic_favors_linear: verd_ctrl = 'AMBIGUOUS' else: verd_ctrl = 'INCONSISTENT' _cache['controlled'] = { 'log_J': lj_c, 'log_gdp': y_c, 'wgi': wgi_c, 'beta_c': beta_c, 'r2_c': r2_c, 'b2_ctrl': b2_ctrl, 'p_b2_ctrl': p_b2_ctrl, 'b2_uncontrolled': beta_q[2], 'p_b2_uncontrolled': p_q[2], 'cooks_d_ctrl': cooks_c, 'iso3_ctrl': ctrl['iso3'].values, } results.append(log_result( 'variety_growth_controlled_quadratic', b2_ctrl, p_b2_ctrl, verd_ctrl, details={ 'N': n_c, 'R2_controlled': r2_c, 'b2_uncontrolled': float(beta_q[2]), 'p_b2_uncontrolled': float(p_q[2]), 'wgi_coeff': float(beta_c[3]), 'wgi_p': float(p_c[3]), 'bic_linear_wgi': float(bic_lin), 'bic_quad_wgi': float(bic_quad), 'bic_favors_linear': bic_favors_linear, 'max_cooks_d': max_c_d, 'max_cooks_country': max_c_iso, 'outlier_driven': outlier_driven_ctrl, 'b2_drop_outlier': float(b2_loo_ctrl), 'p_b2_drop_outlier': float(p_b2_loo_ctrl), }, description=( "Tests whether CES constant-elasticity rejection (Test 3) " "is driven by omitted institutional quality or outliers. " "Regression: log(GDP_pc) = a + b1*log(J) + b2*log(J)^2 + " "b3*WGI. Includes BIC comparison against linear+WGI model " "and Cook's distance outlier analysis. CONSISTENT if b2 " "insignificant. AMBIGUOUS if b2 significant but either " "outlier-driven or BIC favors the simpler linear model." ))) # ================================================================ # Test 5: Panel fixed-effects levels elasticity # Within-country variation isolates pure variety effect from # time-invariant institutional quality. # ================================================================ if ts_df is not None and gdp_panel is not None: panel_rows = [] for c in ts_df['iso3'].unique(): c_var = ts_df[ts_df['iso3'] == c].sort_values('year') c_gdp = gdp_panel[gdp_panel['iso3'] == c].sort_values('year') for _, row in c_var.iterrows(): yr = row['year'] j = row['n_export_products'] g = c_gdp[c_gdp['year'] == yr] if len(g) == 0 or j <= 0: continue gdp_val = g.iloc[0]['gdp_pc'] if gdp_val <= 0: continue panel_rows.append({ 'iso3': c, 'year': yr, 'log_J': np.log(j), 'log_gdp': np.log(gdp_val), }) if len(panel_rows) >= 20: pdf = pd.DataFrame(panel_rows) n_countries = pdf['iso3'].nunique() if n_countries >= 5: # Within-transformation (entity demeaning) means = pdf.groupby('iso3')[['log_J', 'log_gdp']].transform('mean') lj_dm = (pdf['log_J'] - means['log_J']).values lgdp_dm = (pdf['log_gdp'] - means['log_gdp']).values X_fe = np.column_stack([lj_dm]) # No intercept after demeaning n_fe = len(lj_dm) k_fe = 1 beta_fe = np.linalg.lstsq(X_fe, lgdp_dm, rcond=None)[0] resid_fe = lgdp_dm - X_fe @ beta_fe sse_fe = np.sum(resid_fe ** 2) sst_fe = np.sum(lgdp_dm ** 2) r2_fe = 1 - sse_fe / sst_fe if sst_fe > 0 else 0 # Cluster-robust SE (by country) countries = pdf['iso3'].values unique_c = np.unique(countries) meat = np.zeros((k_fe, k_fe)) for c in unique_c: mask = countries == c Xi = X_fe[mask] ei = resid_fe[mask] score = Xi.T @ ei meat += np.outer(score, score) bread_fe = np.linalg.pinv(X_fe.T @ X_fe) # Small-sample correction: G/(G-1) * (N-1)/(N-k) G = len(unique_c) N = n_fe correction = (G / (G - 1)) * ((N - 1) / (N - k_fe)) V_cluster = correction * bread_fe @ meat @ bread_fe se_fe = np.sqrt(np.diag(V_cluster)) t_fe = beta_fe / se_fe # df = G - 1 for cluster-robust p_fe = 2 * (1 - stats.t.cdf(np.abs(t_fe), G - 1)) _cache['panel_fe'] = { 'lj_dm': lj_dm, 'lgdp_dm': lgdp_dm, 'beta': float(beta_fe[0]), 'se': float(se_fe[0]), 'r2': r2_fe, 'n_countries': n_countries, 'n_obs': n_fe, } # Prediction: beta > 0 (variety increases output within country) verd_fe = verdict(float(p_fe[0]), bool(beta_fe[0] > 0)) results.append(log_result( 'variety_growth_panel_fe', float(beta_fe[0]), float(p_fe[0]), verd_fe, details={ 'N_obs': n_fe, 'N_countries': n_countries, 'R2_within': r2_fe, 'se_cluster': float(se_fe[0]), }, description=( "Tests variety-output elasticity using within-country " "variation (panel fixed effects). Demeaned OLS of " "log(GDP_pc) on log(J) with country FE, cluster-robust " "SE by country. Absorbs time-invariant institutional " "quality. Data: variety_timeseries.csv (14 countries, " "6 time points 1995-2019) merged with wb_gdp_pc_all.csv. " "Prediction: beta > 0." ))) _cache['results'] = results return results def figures(): """Generate figures for variety-growth tests.""" if not _cache: main() paths = [] # --- Test 1: Levels elasticity scatter --- d = _cache.get('levels') if d is not None: fig, ax = setup_figure() ax.scatter(d['log_J'], d['log_gdp'], alpha=0.6, s=25, color='#2196F3') x_fit = np.linspace(d['log_J'].min(), d['log_J'].max(), 100) y_fit = np.mean(d['log_gdp']) + d['beta'] * (x_fit - np.mean(d['log_J'])) ax.plot(x_fit, y_fit, 'r-', linewidth=2, label=f'β={d["beta"]:.2f}, R²={d["r2"]:.2f}') ax.set_xlabel('log(Export Products)') ax.set_ylabel('log(GDP per capita)') ax.set_title('Variety-Output Elasticity (Levels)') ax.legend() paths.extend(save_figure(fig, 'test_variety_growth_levels', subfolder='entry_exit')) # --- Test 2: Growth elasticity scatter --- d = _cache.get('growth') if d is not None: fig, ax = setup_figure() ax.scatter(d['dlog_J'], d['dlog_gdp'], alpha=0.6, s=25, color='#4CAF50') ax.axhline(y=0, color='gray', linestyle=':', linewidth=0.5) ax.axvline(x=0, color='gray', linestyle=':', linewidth=0.5) ax.set_xlabel('Δlog(Export Products) / year') ax.set_ylabel('Δlog(GDP per capita) / year') ax.set_title(f'Semi-Endogenous Growth\nβ_growth={d["beta"]:.3f}, R²={d["r2"]:.3f}') paths.extend(save_figure(fig, 'test_variety_growth_semi_endogenous', subfolder='entry_exit')) # --- Test 3: Quadratic fit with outlier highlight --- d = _cache.get('quadratic') if d is not None: fig, ax = setup_figure() cooks = d.get('cooks_d') iso3s = d.get('iso3') if cooks is not None: threshold = 4.0 / len(d['log_J']) normal = cooks <= threshold outlier = cooks > threshold ax.scatter(d['log_J'][normal], d['log_gdp'][normal], alpha=0.6, s=25, color='#FF9800', label='Normal') if np.any(outlier): ax.scatter(d['log_J'][outlier], d['log_gdp'][outlier], alpha=0.9, s=80, color='red', marker='D', edgecolors='black', zorder=5, label='High influence') if iso3s is not None: for idx in np.where(outlier)[0]: ax.annotate(iso3s[idx], (d['log_J'][idx], d['log_gdp'][idx]), fontsize=8, fontweight='bold', xytext=(5, 5), textcoords='offset points') else: ax.scatter(d['log_J'], d['log_gdp'], alpha=0.6, s=25, color='#FF9800') x_fit = np.linspace(d['log_J'].min(), d['log_J'].max(), 100) y_quad = d['beta_q'][0] + d['beta_q'][1] * x_fit + d['beta_q'][2] * x_fit**2 y_lin = d['beta_q'][0] + d['beta_q'][1] * np.mean(d['log_J']) + d['beta_q'][2] * np.mean(d['log_J'])**2 ax.plot(x_fit, y_quad, 'r-', linewidth=2, label=f'Quadratic: b₂={d["b2"]:.2f}, p={d["p_b2"]:.1e}') ax.set_xlabel('log(Export Products)') ax.set_ylabel('log(GDP per capita)') ax.set_title('Constant Elasticity Test (CES predicts b₂=0)') ax.legend(fontsize=8) paths.extend(save_figure(fig, 'test_variety_growth_quadratic', subfolder='entry_exit')) # --- Test 3b: Controlled quadratic --- d = _cache.get('controlled') if d is not None: fig, ax = setup_figure() ax.bar(['Uncontrolled', 'WGI-Controlled'], [d['b2_uncontrolled'], d['b2_ctrl']], color=['#F44336', '#4CAF50'], alpha=0.8) ax.axhline(y=0, color='gray', linestyle='-', linewidth=0.5) ax.set_ylabel('Quadratic coefficient b₂') ax.set_title(f'Institutional Control Test\n' f'b₂: {d["b2_uncontrolled"]:.2f} (p={d["p_b2_uncontrolled"]:.1e})' f' → {d["b2_ctrl"]:.2f} (p={d["p_b2_ctrl"]:.2f})') paths.extend(save_figure(fig, 'test_variety_growth_controlled', subfolder='entry_exit')) # --- Test 5: Panel FE scatter --- d = _cache.get('panel_fe') if d is not None: fig, ax = setup_figure() ax.scatter(d['lj_dm'], d['lgdp_dm'], alpha=0.6, s=30, color='#009688') x_range = np.linspace(d['lj_dm'].min(), d['lj_dm'].max(), 100) ax.plot(x_range, d['beta'] * x_range, 'r-', linewidth=2, label=f'β={d["beta"]:.2f} ± {d["se"]:.2f}') ax.axhline(y=0, color='gray', linestyle=':', linewidth=0.5) ax.axvline(x=0, color='gray', linestyle=':', linewidth=0.5) ax.set_xlabel('Demeaned log(J)') ax.set_ylabel('Demeaned log(GDP per capita)') ax.set_title(f'Panel FE: Within-Country Variety Effect\n' f'{d["n_countries"]} countries, {d["n_obs"]} obs') ax.legend() paths.extend(save_figure(fig, 'test_variety_growth_panel_fe', subfolder='entry_exit')) # --- Test 4: BW benchmark histogram --- bw_s = _cache.get('bw_sigmas') implied = _cache.get('bw_implied') agg = _cache.get('bw_aggregate') if bw_s is not None and implied is not None: fig, ax = setup_figure() ax.hist(bw_s[bw_s < 15], bins=80, density=True, alpha=0.7, color='#9C27B0', edgecolor='white', linewidth=0.3) ax.axvline(x=implied, color='red', linewidth=2, linestyle='-', label=f'Macro implied σ = {implied:.2f}') if agg is not None: ax.axvline(x=agg, color='blue', linewidth=2, linestyle='--', label=f'CES aggregate σ = {agg:.2f}') ax.set_xlabel('Elasticity of Substitution (σ)') ax.set_ylabel('Density') ax.set_title('Broda-Weinstein Benchmark') ax.set_xlim(0, 15) ax.legend() paths.extend(save_figure(fig, 'test_variety_growth_bw_benchmark', 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']})")

bw iqr	[1.7959, 4.3370]
bw median	2.6053
statistic	1.5083
N products	2680
implied sigma	1.5083

Variety Growth Bw Benchmark

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources