#!/usr/bin/env python3 """ Estimation Paradox (Theorem 7 of "CES as Theorem") ==================================================== Tests that Fisher information I(rho) = Var_P[log x_j] = 0 at symmetry. Countries with more balanced sectors should yield LESS precise sigma estimates. Uses OECD STAN sectoral data for cross-sectional estimation: log(s_j) = alpha + (1-sigma)*log(P_j) + eps per country-year. Review Notes (2026-03-12): Theory match: CORRECT — tests Theorem thm:paradox (Paper 1) I(rho) = Var_P[log x_j] = 0 at symmetry. Lean: fisherInfoRho_zero_at_symmetry (real proof, not schematic). I/O audit: CORRECT — nominal shares for CES regression, deflators as prices, Var[log real VA] as Fisher proxy. Sign conventions all correct. Data-signal alignment: ADEQUATE — 775 country-year obs, 38 OECD countries, 1995-2019. Fisher proxy (Var[log real VA]) is direct. Balance (1/CV) is indirect. Code quality: CLEAN — MSE control correctly isolates information channel. Confidence: MEDIUM — Tests 1 and 3 CONSISTENT with clustered SEs (p<0.001); test 2 AMBIGUOUS (clustered p=0.343, SE inflated 3.6x by within-country correlation). Country FE reverses the sign for balance (test 3) and var_log_shares (test 2), but Fisher proxy (test 1) is robust to FE (slope -0.35, t=-2.64, p=0.009). Known limitations: - MSE confound: Spearman(balance, MSE) = -0.645. MSE control resolves this but only partially in the within-country dimension. - Country FE reversal: within-country balance changes show OPPOSITE sign for tests 2-3. The pooled result is driven by between-country variation. - No year FE: common year shocks could affect both balance and SE. The companion MLE test (test_estimation_paradox_mle) avoids the MSE confound entirely by using the escort-distribution MLE whose SE follows the Cramer-Rao bound directly. Migrated from scripts/test_estimation_paradox.py """ import numpy as np import pandas as pd from scipy import stats from common.data import load_oecd_stan from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure # Non-overlapping ISIC Rev.4 sectors 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): """Parse OECD STAN columns into usable format (handles raw and pre-processed).""" 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 estimate_sigma_cross_sectional(nom_panels, defl_panels, real_panels): rows = [] for country in nom_panels: if country not in defl_panels or country not in real_panels: continue nom, defl, real = nom_panels[country], defl_panels[country], real_panels[country] common_sectors = sorted(set(nom.columns) & set(defl.columns) & set(real.columns)) common_years = sorted(set(nom.index) & set(defl.index) & set(real.index)) if len(common_sectors) < 8: continue for year in common_years: nom_row = nom.loc[year, common_sectors] defl_row = defl.loc[year, common_sectors] real_row = real.loc[year, common_sectors] valid = (nom_row > 0) & (defl_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) dv = defl_row[sv].values.astype(float) rv = real_row[sv].values.astype(float) J = len(sv) shares = nv / nv.sum() log_s = np.log(shares) log_p = np.log(dv) x = log_p - log_p.mean() y = log_s - log_s.mean() xx = np.sum(x**2) if xx < 1e-15: continue beta = np.sum(x * y) / xx resid = y - beta * x dof = J - 2 if dof <= 0: continue s2 = np.sum(resid**2) / dof se_beta = np.sqrt(s2 / xx) fisher_proxy = np.var(np.log(rv)) cv_shares = np.std(shares) / np.mean(shares) balance = 1.0 / cv_shares if cv_shares > 0 else np.nan rows.append({ 'country': country, 'year': year, 'se_sigma': se_beta, 'log_se': np.log(se_beta) if se_beta > 0 else np.nan, '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, 'var_log_shares': np.var(log_s), 'log_var_log_s': np.log(np.var(log_s)) if np.var(log_s) > 0 else np.nan, 'log_MSE': np.log(s2) if s2 > 0 else np.nan, 'n_sectors': J, }) return pd.DataFrame(rows) _cache = {} def main(): try: stan_raw = load_oecd_stan() except FileNotFoundError: return [log_result('estimation_paradox', np.nan, np.nan, 'SKIPPED', description="Cross-sectional estimation paradox tests require OECD STAN sectoral value added data.")] stan = parse_stan(stan_raw) nom_va = extract_b1g(stan, 'V') real_va = extract_b1g(stan, 'L') deflators = extract_b1g(stan, 'D') nom_panels = build_panels(nom_va, "Nominal") real_panels = build_panels(real_va, "Real") defl_panels = build_panels(deflators, "Deflator") cs_df = estimate_sigma_cross_sectional(nom_panels, defl_panels, real_panels) cs_df = cs_df.dropna(subset=['log_se', 'log_fisher', 'log_balance', 'log_MSE']) if len(cs_df) < 20: return [log_result('estimation_paradox', np.nan, np.nan, 'SKIPPED', description="Cross-sectional estimation paradox tests require at least 20 valid country-year observations.")] results = [] n_obs = len(cs_df) countries = cs_df['country'].values unique_countries = np.unique(countries) n_clusters = len(unique_countries) def ols_with_mse_control(x_col, y_col='log_se'): """Multiple regression y ~ 1 + x + log(MSE) with clustered SEs by country. Returns (slope, clustered_p, details_dict) where details_dict contains OLS SE, clustered SE, and country-FE robustness results. """ X = np.column_stack([np.ones(n_obs), cs_df[x_col].values, cs_df['log_MSE'].values]) y = cs_df[y_col].values betas = np.linalg.lstsq(X, y, rcond=None)[0] resids = y - X @ betas n, k = n_obs, X.shape[1] # OLS SE bread = np.linalg.inv(X.T @ X) ols_se = np.sqrt(np.diag(np.sum(resids**2) / (n - 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 / (n - 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: within-demean x_col, log_MSE, y_col fe_df = cs_df[[x_col, 'log_MSE', y_col, 'country']].copy() for col in [x_col, 'log_MSE', y_col]: fe_df[col] = fe_df.groupby('country')[col].transform(lambda s: s - s.mean()) X_fe = np.column_stack([np.ones(n), fe_df[x_col].values, fe_df['log_MSE'].values]) y_fe = fe_df[y_col].values betas_fe = np.linalg.lstsq(X_fe, y_fe, rcond=None)[0] resids_fe = y_fe - X_fe @ betas_fe dof_fe = n - k - n_clusters + 1 bread_fe = np.linalg.inv(X_fe.T @ X_fe) se_fe = np.sqrt(np.diag(np.sum(resids_fe**2) / max(dof_fe, 1) * bread_fe)) t_fe = betas_fe[1] / se_fe[1] if se_fe[1] > 0 else 0.0 p_fe = 2 * (1 - stats.t.cdf(abs(t_fe), max(dof_fe, 1))) detail = { 'n_obs': n, 'n_clusters': n_clusters, 'ols_se': float(ols_se[1]), 'clustered_se': float(cl_se[1]), 'clustered_t': float(t_cl), 'fe_slope': float(betas_fe[1]), 'fe_t': float(t_fe), 'fe_p': float(p_fe), } return betas[1], p_cl, detail # Test 1: log(SE) ~ log(Fisher) + log(MSE) slope1, p1, d1 = ols_with_mse_control('log_fisher') v1 = 'CONSISTENT' if slope1 < 0 and p1 < 0.05 else ('AMBIGUOUS' if slope1 < 0 else 'INCONSISTENT') results.append(log_result('estimation_paradox_fisher_to_se', slope1, p1, v1, details=d1, description="Tests whether higher Fisher information (Var[log real VA]) yields lower SE of sigma. " "Multiple regression of log(SE) on log(Fisher proxy) + log(MSE) across country-year " "cross-sectional sigma estimates from OECD STAN data (N=%d, %d countries). " "Controls for residual variance to isolate the information channel. " "SEs clustered by country. FE robustness: slope=%.3f, t=%.2f." % (d1['n_obs'], d1['n_clusters'], d1['fe_slope'], d1['fe_t']))) # Test 2: log(SE) ~ log(Var[log shares]) + log(MSE) slope2, p2, d2 = ols_with_mse_control('log_var_log_s') v2 = 'CONSISTENT' if slope2 < 0 and p2 < 0.05 else ('AMBIGUOUS' if slope2 < 0 else 'INCONSISTENT') results.append(log_result('estimation_paradox_var_logs_to_se', slope2, p2, v2, details=d2, description="Tests the estimation paradox via share dispersion: higher Var[log shares] should reduce SE of sigma. " "Multiple regression of log(SE) on log(Var[log shares]) + log(MSE); " "negative slope confirms that more dispersed shares improve sigma identification " "after controlling for residual variance. SEs clustered by country. " "FE robustness: slope=%.3f, t=%.2f (sign reverses under FE)." % (d2['fe_slope'], d2['fe_t']))) # Test 3: log(SE) ~ log(Balance) + log(MSE) slope3, p3, d3 = ols_with_mse_control('log_balance') v3 = 'CONSISTENT' if slope3 > 0 and p3 < 0.05 else ('AMBIGUOUS' if slope3 > 0 else 'INCONSISTENT') results.append(log_result('estimation_paradox_balance_to_se', slope3, p3, v3, details=d3, description="Tests the paradox directly: more balanced sectors (lower CV of shares) should yield higher SE of sigma. " "Multiple regression of log(SE) on log(1/CV) + log(MSE); positive slope confirms " "that balanced allocation impairs sigma estimation after controlling for residual variance. " "SEs clustered by country. " "FE robustness: slope=%.3f, t=%.2f (sign reverses under FE)." % (d3['fe_slope'], d3['fe_t']))) _cache['cs_df'] = cs_df return results def _partial_residuals(cs_df, x_col, y_col='log_se', control='log_MSE'): """Compute partial residuals: residualize both y and x on control.""" c = cs_df[control].values C = np.column_stack([np.ones(len(c)), c]) C_inv = np.linalg.lstsq(C, np.eye(len(c)), rcond=None)[0] y = cs_df[y_col].values y_resid = y - C @ (np.linalg.lstsq(C, y, rcond=None)[0]) x = cs_df[x_col].values x_resid = x - C @ (np.linalg.lstsq(C, x, rcond=None)[0]) return x_resid, y_resid def figures(): """Generate estimation paradox partial residual plots.""" if not _cache: main() paths = [] cs_df = _cache.get('cs_df') if cs_df is None or len(cs_df) < 20: return paths fig, axes = setup_subplots(1, 3, width=16, height=5) panels = [ ('log_fisher', 'steelblue', '(a) SE vs Fisher info | MSE', 'log(Fisher proxy) | log(MSE)', -1), ('log_var_log_s', 'forestgreen', '(b) SE vs share variance | MSE', 'log(Var[log shares]) | log(MSE)', -1), ('log_balance', 'darkorange', '(c) SE vs balance | MSE', 'log(Balance) | log(MSE)', 1), ] for ax, (x_col, color, title, xlabel, expected_sign) in zip(axes, panels): x_r, y_r = _partial_residuals(cs_df, x_col) ax.scatter(x_r, y_r, s=10, alpha=0.4, c=color) sl, il, _, _, _ = stats.linregress(x_r, y_r) xr = np.linspace(x_r.min(), x_r.max(), 50) ax.plot(xr, il + sl * xr, 'r-', lw=1.5, label=f'partial slope={sl:.3f}') ax.set_xlabel(xlabel) ax.set_ylabel('log(SE) | log(MSE)') ax.set_title(title) ax.legend(fontsize=8) ax.grid(True, alpha=0.3) fig.suptitle('Estimation Paradox: Partial Residual Plots (controlling for MSE)', fontsize=12, y=1.02) paths.extend(save_figure(fig, 'estimation_paradox', 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	6.274e-4
statistic	-0.6468

Estimation Paradox Var Logs To Se

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps