#!/usr/bin/env python3 """ Estimation Paradox via Shares-Only MLE ======================================= Matches observed nominal shares to escort probabilities: pi_j(rho) = x_j^rho / sum x_k^rho Theory (Paper 1, Theorem EST-2, fisherInfoRho_zero_at_symmetry): Fisher information I(rho) = Var_P[log x_j] = 0 at symmetric equilibrium. Countries with more balanced sectors should yield LESS precise rho estimates. Tests: A. log(SE) vs log(Var[log real VA]) — independent Fisher proxy, slope ~ -0.5 B. log(SE) vs log(Balance) — positive slope (balanced → less precise) C. Spearman rank correlation between Var[log real VA] and SE Review Notes (2026-03-13): - Theory match: CORRECT — tests Theorem EST-2 (Paper 1) I(rho) = Var_P[log x] = 0 at symmetry. Lean: fisherInfoRho_zero_at_symmetry (real proof, not schematic). - I/O audit: CORRECT (after fix) — Tests A/C originally used analytical escort Fisher info which is tautologically = 1/SE^2 (the Hessian of the escort LL equals -Var_P[log x] exactly). Fixed to use Var[log real VA] as independent proxy. - Data-signal alignment: ADEQUATE — 813 country-year obs, 38 OECD countries. - Code quality: BUGS FOUND — (1) Tests A/C were tautological (slope=-0.5 and Spearman=-1 by construction). Fixed with independent Fisher proxy. (2) Test B lacked clustered SEs. Fixed. - Confidence: MEDIUM — Test B robust with clustered SEs; FE reversal for within-country variation (same pattern as companion test). Migrated from scripts/test_estimation_paradox_mle.py """ import numpy as np import pandas as pd from scipy import stats, optimize from common.data import load_oecd_stan from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure TARGET_SECTORS = [ 'A', 'B', 'C10T12', 'C13T15', 'C16T18', 'C19', 'C20_21', 'C22_23', 'C24_25', 'C26T28', 'C29_30', 'C31T33', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M_N', 'O', 'P', 'Q' ] SECTOR_ALTERNATIVES = {'D': ['D_E'], 'E': ['D_E'], 'M_N': ['M', 'N']} def parse_stan(stan): stan = stan.copy() # Already pre-processed (from oecd_stan_va.csv) if 'country' in stan.columns and 'activity' in stan.columns: stan['year'] = pd.to_numeric(stan['year'], errors='coerce') stan['value'] = pd.to_numeric(stan['value'], errors='coerce') return stan.dropna(subset=['year', 'value']) # Raw OECD SDMX format stan['country'] = stan['REF_AREA: Reference area'].str.split(':').str[0].str.strip() stan['activity'] = stan['ACTIVITY: Economic activity'].str.split(':').str[0].str.strip() stan['measure'] = stan['MEASURE: Measure'].str.split(':').str[0].str.strip() stan['price_base'] = stan['PRICE_BASE: Price base'].str.split(':').str[0].str.strip() stan['year'] = pd.to_numeric(stan['TIME_PERIOD: Time period'], errors='coerce') stan['value'] = pd.to_numeric(stan['OBS_VALUE'], errors='coerce') return stan.dropna(subset=['year', 'value']) def extract_b1g(stan, price_base_code, year_min=1995, year_max=2019): mask = ((stan['measure'] == 'B1G') & (stan['price_base'] == price_base_code) & (stan['year'] >= year_min) & (stan['year'] <= year_max)) return stan[mask][['country', 'activity', 'year', 'value']].copy() def find_sector(activity, available): if activity in available: return activity for alt in SECTOR_ALTERNATIVES.get(activity, []): if alt in available: return alt return None def build_panels(va_data, label=""): available = set(va_data['activity'].unique()) mapped = list(set(find_sector(s, available) for s in TARGET_SECTORS) - {None}) va_f = va_data[va_data['activity'].isin(mapped)] panels = {} for country in va_f['country'].unique(): cdata = va_f[va_f['country'] == country] pivot = cdata.pivot_table(index='year', columns='activity', values='value', aggfunc='mean') if pivot.shape[1] >= 10 and pivot.shape[0] >= 5: pivot = pivot.dropna(axis=1, how='all').dropna(axis=0, how='any') if pivot.shape[1] >= 10 and pivot.shape[0] >= 5: panels[country] = pivot return panels def escort_log_likelihood(rho, log_x, obs_shares): log_weights = rho * log_x log_Z = np.logaddexp.reduce(log_weights) log_pi = log_weights - log_Z return np.sum(obs_shares * log_pi) def escort_fisher_info(rho, log_x): log_weights = rho * log_x log_Z = np.logaddexp.reduce(log_weights) P = np.exp(log_weights - log_Z) mean_log_x = np.sum(P * log_x) return np.sum(P * (log_x - mean_log_x)**2) def estimate_rho_mle(log_x, obs_shares): def neg_ll(rho): return -escort_log_likelihood(rho, log_x, obs_shares) try: result = optimize.minimize_scalar(neg_ll, bounds=(-5, 5), method='bounded') rho_hat = result.x except Exception: return None, None, None, False h = 1e-5 ll_plus = escort_log_likelihood(rho_hat + h, log_x, obs_shares) ll_center = escort_log_likelihood(rho_hat, log_x, obs_shares) ll_minus = escort_log_likelihood(rho_hat - h, log_x, obs_shares) hessian = (ll_plus - 2*ll_center + ll_minus) / h**2 if hessian >= 0: return rho_hat, None, None, False se_rho = 1.0 / np.sqrt(-hessian) fisher_info = escort_fisher_info(rho_hat, log_x) return rho_hat, se_rho, fisher_info, True _cache = {} def main(): try: stan_raw = load_oecd_stan() except FileNotFoundError: return [log_result('estimation_paradox_mle', np.nan, np.nan, 'SKIPPED', description="Estimation paradox MLE tests require OECD STAN sectoral value added data.")] stan = parse_stan(stan_raw) nominal_va = extract_b1g(stan, 'V') real_va = extract_b1g(stan, 'L') nom_panels = build_panels(nominal_va, "Nominal") real_panels = build_panels(real_va, "Real") rows = [] for country in sorted(nom_panels.keys()): if country not in real_panels: continue nom, real = nom_panels[country], real_panels[country] common_sectors = sorted(set(nom.columns) & set(real.columns)) common_years = sorted(set(nom.index) & set(real.index)) if len(common_sectors) < 8: continue for year in common_years: nom_row = nom.loc[year, common_sectors] real_row = real.loc[year, common_sectors] valid = (nom_row > 0) & (real_row > 0) if valid.sum() < 8: continue sv = [s for s, v in zip(common_sectors, valid) if v] nv = nom_row[sv].values.astype(float) rv = real_row[sv].values.astype(float) obs_shares = nv / nv.sum() log_x = np.log(rv) rho_hat, se_rho, fisher_info, converged = estimate_rho_mle(log_x, obs_shares) if rho_hat is None or se_rho is None or not converged: continue if se_rho <= 0 or fisher_info is None or fisher_info <= 0: continue cv_shares = np.std(obs_shares) / np.mean(obs_shares) balance = 1.0 / cv_shares if cv_shares > 0 else np.nan # Independent Fisher proxy: Var[log real VA] (not the analytical # escort Fisher info, which is tautologically = 1/SE^2) fisher_proxy = np.var(log_x) rows.append({ 'country': country, 'year': year, 'se_rho': se_rho, 'log_se': np.log(se_rho), 'fisher_proxy': fisher_proxy, 'log_fisher': np.log(fisher_proxy) if fisher_proxy > 0 else np.nan, 'balance': balance, 'log_balance': np.log(balance) if balance > 0 and np.isfinite(balance) else np.nan, }) df = pd.DataFrame(rows) df = df.dropna(subset=['log_se', 'log_fisher', 'log_balance']) if len(df) < 20: return [log_result('estimation_paradox_mle', np.nan, np.nan, 'SKIPPED', description="Estimation paradox MLE tests require at least 20 valid country-year observations.")] results = [] n_obs = len(df) countries = df['country'].values unique_countries = np.unique(countries) n_clusters = len(unique_countries) def _clustered_ols(x_col, y_col='log_se'): """OLS with clustered SEs by country. Returns slope, clustered_p, details.""" X = np.column_stack([np.ones(n_obs), df[x_col].values]) y = df[y_col].values betas = np.linalg.lstsq(X, y, rcond=None)[0] resids = y - X @ betas k = 2 bread = np.linalg.inv(X.T @ X) ols_se = np.sqrt(np.diag(np.sum(resids**2) / (n_obs - k) * bread)) # Clustered SE (Cameron-Gelbach-Miller sandwich) meat = np.zeros((k, k)) for c in unique_countries: mask = countries == c score = X[mask].T @ resids[mask] meat += np.outer(score, score) V_cl = bread @ meat @ bread * n_clusters / (n_clusters - 1) * n_obs / (n_obs - k) cl_se = np.sqrt(np.diag(V_cl)) t_cl = betas[1] / cl_se[1] p_cl = 2 * (1 - stats.t.cdf(abs(t_cl), n_clusters - 1)) # Country FE robustness fe_df = df[[x_col, y_col, 'country']].copy() for col in [x_col, y_col]: fe_df[col] = fe_df.groupby('country')[col].transform(lambda s: s - s.mean()) sl_fe, _, _, p_fe, _ = stats.linregress(fe_df[x_col], fe_df[y_col]) return betas[1], p_cl, { 'n_obs': n_obs, 'n_clusters': n_clusters, 'ols_se': float(ols_se[1]), 'clustered_se': float(cl_se[1]), 'fe_slope': float(sl_fe), 'fe_p': float(p_fe), } # Test A: log(SE) vs log(Fisher proxy) -- slope should be negative # Fisher proxy = Var[log real VA], an INDEPENDENT measure of input dispersion # (NOT the analytical escort Fisher info, which is = 1/SE^2 by construction) slope_a, p_a, d_a = _clustered_ols('log_fisher') va = 'CONSISTENT' if slope_a < 0 and p_a < 0.05 else ('AMBIGUOUS' if slope_a < 0 else 'INCONSISTENT') results.append(log_result('estimation_paradox_mle_cramer_rao_slope', slope_a, p_a, va, details=d_a, description="Tests the Cramer-Rao prediction: higher cross-sectional dispersion of real VA " "(independent Fisher proxy = Var[log real VA]) should yield lower SE of rho. " "OLS of log(SE) on log(Var[log real VA]) across %d country-year escort MLE estimates " "(%d countries). Clustered SEs by country. FE slope=%.3f, p=%.4f." % (d_a['n_obs'], d_a['n_clusters'], d_a['fe_slope'], d_a['fe_p']))) # Test B: log(SE) vs log(Balance) -- positive slope slope_b, p_b, d_b = _clustered_ols('log_balance') vb = 'CONSISTENT' if slope_b > 0 and p_b < 0.05 else ('AMBIGUOUS' if slope_b > 0 else 'INCONSISTENT') results.append(log_result('estimation_paradox_mle_balance_to_se', slope_b, p_b, vb, details=d_b, description="Tests the estimation paradox: more balanced sectoral allocation (lower CV of shares) " "should yield higher SE of rho. OLS of log(SE) on log(1/CV) across %d country-year " "escort MLE estimates (%d countries); positive slope confirms balanced allocation " "impairs rho estimation. Clustered SEs by country. FE slope=%.3f, p=%.4f." % (d_b['n_obs'], d_b['n_clusters'], d_b['fe_slope'], d_b['fe_p']))) # Test C: Spearman rank correlation (Fisher proxy vs SE) rho_sf, p_sf = stats.spearmanr(df['fisher_proxy'], df['se_rho']) vc = 'CONSISTENT' if rho_sf < 0 and p_sf < 0.05 else ('AMBIGUOUS' if rho_sf < 0 else 'INCONSISTENT') results.append(log_result('estimation_paradox_mle_spearman_fisher_se', rho_sf, p_sf, vc, description="Nonparametric test: Spearman rank correlation between Var[log real VA] " "(independent Fisher proxy) and SE of rho across country-year MLE estimates. " "Negative correlation confirms that higher cross-sectional input dispersion " "yields tighter rho estimates, as predicted by the estimation paradox.")) _cache['df'] = df return results def figures(): """Generate estimation paradox MLE figure.""" if not _cache: main() paths = [] df = _cache.get('df') if df is None or len(df) < 20: return paths fig, axes = setup_subplots(1, 3, width=16, height=5) # (a) log(SE) vs log(Fisher proxy = Var[log real VA]) ax = axes[0] ax.scatter(df['log_fisher'], df['log_se'], s=10, alpha=0.4, c='steelblue') from scipy import stats as sp_stats sl, il, _, _, _ = sp_stats.linregress(df['log_fisher'], df['log_se']) xr = np.linspace(df['log_fisher'].min(), df['log_fisher'].max(), 50) ax.plot(xr, il + sl * xr, 'r-', label=f'slope={sl:.3f}') ax.set_xlabel('log(Var[log real VA])') ax.set_ylabel('log(SE of rho)') ax.set_title('(a) SE vs Fisher proxy') ax.legend(fontsize=8) ax.grid(True, alpha=0.3) # (b) log(SE) vs log(Balance) ax = axes[1] ax.scatter(df['log_balance'], df['log_se'], s=10, alpha=0.4, c='darkorange') sl2, il2, _, _, _ = sp_stats.linregress(df['log_balance'], df['log_se']) xr2 = np.linspace(df['log_balance'].min(), df['log_balance'].max(), 50) ax.plot(xr2, il2 + sl2 * xr2, 'r-', label=f'slope={sl2:.3f}') ax.set_xlabel('log(Balance)') ax.set_ylabel('log(SE of rho)') ax.set_title('(b) SE vs sectoral balance') ax.legend(fontsize=8) ax.grid(True, alpha=0.3) # (c) Fisher proxy vs SE (levels, rank correlation) ax = axes[2] ax.scatter(df['fisher_proxy'], df['se_rho'], s=10, alpha=0.4, c='forestgreen') ax.set_xlabel('Var[log real VA]') ax.set_ylabel('SE of rho') ax.set_title('(c) Spearman rank: Fisher proxy vs SE') ax.grid(True, alpha=0.3) fig.suptitle('Estimation Paradox: MLE (Escort Distribution)', fontsize=13, y=1.02) paths.extend(save_figure(fig, 'estimation_paradox_mle', subfolder='foundations')) 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.546e-175
statistic	0.6572

Estimation Paradox MLE Balance To Se

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps