#!/usr/bin/env python3 """ Test: CES Isotropy Overidentification (U1) Paper 3 (Dynamics on the CES Landscape) / Paper 5 (Empirical Tests) Theory (Corollary 10.2, eq. 18): T_parallel = CV2(output growth)_n [rho-free longitudinal T estimator] T_perp = K_eff,n * Var(s_nj) [transverse T from input shares] Isotropy: T_parallel = T_perp Test: log(T_parallel) = alpha + beta * log(K_eff * Var(s)) + eps H0: beta = 1 (isotropy), rank order correct (Spearman > 0) Four approaches: A) BEA full I-O (79 commodities x 71 industries) + BEA quarterly real VA B) Euler-VRI at 4-digit NAICS: Census ECNMATFUEL + FRED IP (materials-only) C) Overidentification with BEA I-O 2022 + BEA GDP quarterly D) BEA 2017 Benchmark detailed Use Table (402 industries × 385 commodities) Data: BEA Use Table: data/bea/io/Use_Summary_259_2017.csv BEA GDP-by-Industry: data/bea/gdp_industry/RealValueAdded_10_quarterly.csv BEA I-O 2022: data/compiled/bea_io_use_2022.csv BEA GDP quarterly: data/compiled/bea_gdp_industry_quarterly.csv BEA Benchmark 2017: data/compiled/bea_use_detail_2017.csv (from Excel) Census ECNMATFUEL: data/compiled/ecnmatfuel_*.csv FRED 4-digit IP: data/fred/IPG????N.csv Review Notes (2026-03-13) ------------------------- Theory match: CORRECT. Directly implements the VRI eigenspace decomposition (Paper 3, Corollary): T_parallel = CV^2(HP-filtered log output), T_perp = K * Var(s_nj). Isotropy predicts beta=1 in log-log regression. I/O audit: CORRECT. T_parallel from output VA, T_perp from input shares. Data-signal alignment: ADEQUATE to STRONG. - BEA Use Table input shares: DIRECT measurement of intermediate input purchases - BEA quarterly real VA: DIRECT measurement of sector output - Oberfield-Raval sigma: published, well-cited calibrations - Benchmark 2017 detailed table provides 49:1 Var(s) variation vs 1.4:1 in the summary table Code quality: CLEAN (after 2026-03-13 fixes) - All approaches now use HC0-robust standard errors - Manufacturing verdict fallthrough fixed (can now return INCONSISTENT) - Approach C bug fixed: was including tax/VA rows as inputs (T005, VABAS etc.), giving beta=2.9. After fix: beta=0.51, N expanded from 14 to 66. - Var(s) computation standardized to np.mean((s - 1/J)^2) Confidence: MEDIUM (rank ordering robust, beta < 1 in 4/5 subtests) """ import warnings warnings.filterwarnings("ignore") import csv import numpy as np import pandas as pd import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.config import (BEA_IO_DIR, BEA_GDP_DIR, COMPILED_DIR, FRED_DIR, FRED_API_KEY, DATA_DIR) from common.stats import log_result, hp_filter from common.ces import curvature_K, rho_from_sigma from common.plotting import setup_figure, setup_subplots, save_figure try: from scipy import stats as scipy_stats HAS_SCIPY = True except ImportError: HAS_SCIPY = False try: import requests HAS_REQUESTS = True except ImportError: HAS_REQUESTS = False # ── Sigma calibrations ─────────────────────────────────────────────────── # BEA industry code -> sigma (Oberfield-Raval + broader literature) SIGMA_MAP_BEA = { '111CA': 0.65, '113FF': 0.65, '211': 0.90, '212': 0.55, '213': 0.65, '22': 0.60, '23': 0.70, '311FT': 1.15, '313TT': 0.70, '315AL': 0.70, '321': 0.70, '322': 0.85, '323': 0.85, '324': 0.90, '325': 0.90, '326': 0.95, '327': 0.60, '331': 0.55, '332': 0.70, '333': 0.65, '334': 1.40, '335': 0.70, '3361MV': 0.50, '3364OT': 0.50, '337': 0.85, '339': 1.10, '42': 1.10, '441': 1.10, '445': 1.10, '452': 1.10, '4A0': 1.10, '481': 0.90, '482': 0.60, '483': 0.90, '484': 0.90, '485': 0.90, '486': 0.60, '487OS': 0.90, '493': 0.90, '511': 1.40, '512': 1.20, '513': 1.40, '514': 1.40, '521CI': 1.20, '523': 1.20, '524': 1.20, '525': 1.20, '532RL': 0.80, 'HS': 0.80, 'ORE': 0.80, '5411': 1.10, '5412OP': 1.10, '5415': 1.40, '55': 1.10, '561': 1.10, '562': 0.70, '61': 0.80, '621': 0.85, '622': 0.85, '623': 0.85, '624': 0.85, '711AS': 1.10, '713': 1.10, '721': 1.20, '722': 1.20, '81': 1.10, } GOVT_CODES = {'GFE', 'GFGD', 'GFGN', 'GSLE', 'GSLG'} NON_COMMODITY_PREFIXES = ('T', 'V', 'PROF') # OR 3-digit sigma for 4-digit Euler-VRI test OR_SIGMA_3DIG = { '311': 1.15, '312': 0.85, '313': 0.70, '314': 0.70, '315': 0.70, '316': 0.70, '321': 0.70, '322': 0.85, '323': 0.85, '324': 0.90, '325': 0.90, '326': 0.95, '327': 0.60, '331': 0.55, '332': 0.70, '333': 0.65, '334': 1.40, '335': 0.70, '336': 0.50, '337': 0.85, '339': 1.10, } # OR sectors for overidentification test OR_SECTORS = { '336': ('Transport Equip', 0.50, '3361MV'), '331': ('Primary Metals', 0.55, '331'), '327': ('Nonmetallic Min', 0.60, '327'), '333': ('Machinery', 0.65, '333'), '321': ('Wood', 0.70, '321'), '332': ('Fabricated Metals', 0.70, '332'), '335': ('Electrical Equip', 0.70, '335'), '322': ('Paper', 0.85, '322'), '337': ('Furniture', 0.85, '337'), '325': ('Chemicals', 0.90, '325'), '326': ('Plastics/Rubber', 0.95, '326'), '339': ('Miscellaneous', 1.10, '339'), '311': ('Food', 1.15, '311FT'), '334': ('Computer/Elec', 1.40, '334'), } # ── OLS with HC0/HC1 ───────────────────────────────────────────────────── def _ols_hc(y, x, hc='hc0'): """OLS of y on [1, x]. Returns (alpha, beta, se_beta, r2, resid).""" n = len(y) X = np.column_stack([np.ones(n), x]) coef, *_ = np.linalg.lstsq(X, y, rcond=None) alpha_hat, beta_hat = coef resid = y - X @ coef XtX_inv = np.linalg.inv(X.T @ X) if hc == 'hc1': k = X.shape[1] S = X.T @ np.diag(resid ** 2) @ X * n / (n - k) else: S = X.T @ np.diag(resid ** 2) @ X cov = XtX_inv @ S @ XtX_inv se_b = np.sqrt(cov[1, 1]) ss_res = float(resid @ resid) ss_tot = float(np.sum((y - y.mean()) ** 2)) r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0 return alpha_hat, beta_hat, se_b, r2, resid # ══════════════════════════════════════════════════════════════════════════ # APPROACH A: BEA full I-O + quarterly VA # ══════════════════════════════════════════════════════════════════════════ def _approach_a(): """BEA Use Table (71 industries) + quarterly real VA.""" results = [] # Load Use Table use_path = BEA_IO_DIR / 'Use_Summary_259_2017.csv' if not use_path.exists(): return results use = pd.read_csv(use_path) all_row_codes = use.RowCode.unique() commodity_codes = [c for c in all_row_codes if isinstance(c, str) and not any(c.startswith(p) for p in NON_COMMODITY_PREFIXES) and c not in ('F010', 'F02E', 'F02N', 'F02R', 'F02S', 'F030', 'F040', 'F06C', 'F06E', 'F06N', 'F06S', 'F07C', 'F07E', 'F07N', 'F07S', 'F10C', 'F10E', 'F10N', 'F10S', 'T001')] industry_codes = [c for c in use.ColCode.unique() if isinstance(c, str) and c not in GOVT_CODES and not c.startswith('F') and not c.startswith('T')] use_data = {} for ind in industry_codes: if ind not in SIGMA_MAP_BEA: continue mask = (use.ColCode == ind) & (use.RowCode.isin(commodity_codes)) ind_data = use[mask].copy() ind_data['val'] = pd.to_numeric(ind_data.DataValue, errors='coerce') ind_data = ind_data.dropna(subset=['val']) ind_data = ind_data[ind_data.val > 0] if len(ind_data) < 3: continue total = ind_data.val.sum() shares = (ind_data.val / total).values J = len(shares) mean_share = 1.0 / J var_s = np.mean((shares - mean_share) ** 2) sigma = SIGMA_MAP_BEA[ind] rho = rho_from_sigma(sigma) K_eff = curvature_K(rho, J) use_data[ind] = { 'sigma': sigma, 'rho': rho, 'J': J, 'K_eff': K_eff, 'var_s': var_s, 'K_var_s': K_eff * var_s, } # Load quarterly VA va_path = BEA_GDP_DIR / 'RealValueAdded_10_quarterly.csv' if not va_path.exists(): return results va = pd.read_csv(va_path) va['val'] = pd.to_numeric(va.DataValue, errors='coerce') va = va.dropna(subset=['val']) va = va[va.val > 0] va_data = {} for ind in va.Industry.unique(): ind_data = va[va.Industry == ind].sort_values(['Year', 'Quarter']) y = ind_data.val.values if len(y) < 20: continue log_y = np.log(y) cycle, trend = hp_filter(log_y, lamb=1600) T_par = np.var(cycle) if T_par > 0: va_data[ind] = {'T_parallel': T_par} # Match common = set(use_data.keys()) & set(va_data.keys()) - GOVT_CODES if len(common) < 10: return results industries = sorted(common) T_par = np.array([va_data[i]['T_parallel'] for i in industries]) K_var_s = np.array([use_data[i]['K_var_s'] for i in industries]) valid = (T_par > 0) & (K_var_s > 0) industries = [i for i, v in zip(industries, valid) if v] T_par = T_par[valid] K_var_s = K_var_s[valid] N = len(industries) if N < 10 or not HAS_SCIPY: return results log_T = np.log(T_par) log_Kv = np.log(K_var_s) alpha_hat, slope, std_err, r2, resid = _ols_hc(log_T, log_Kv) r_spearman, p_spearman = scipy_stats.spearmanr(log_Kv, log_T) beta_lo = slope - 1.96 * std_err beta_hi = slope + 1.96 * std_err if r_spearman > 0.5 and p_spearman < 0.01 and beta_lo <= 1.0 <= beta_hi: verd = 'CONSISTENT' elif r_spearman > 0.3 and p_spearman < 0.05: verd = 'AMBIGUOUS' elif r_spearman > 0 and p_spearman < 0.10: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' results.append(log_result( name='ces_isotropy_bea_full_io', stat=float(slope), p=float(p_spearman), verd=verd, description="Tests isotropy prediction T_parallel = T_perp across 66 BEA industries. " "OLS of log(HP-cycle variance) on log(K_eff * Var(input shares)) with H0: beta=1. " "Uses BEA Use Table input shares and quarterly real value added. HC0-robust standard errors.", details={'n': N, 'beta': float(slope), 'se': float(std_err), 'spearman_r': float(r_spearman), 'r2': float(r2), 'beta_ci': [float(beta_lo), float(beta_hi)]}, )) print(f" A: N={N}, beta={slope:.3f} (SE={std_err:.3f}), " f"Spearman r={r_spearman:.3f} (p={p_spearman:.4f}) -> {verd}") _cache['records_a'] = [{'code': i, 'T_par': float(t), 'K_var_s': float(k)} for i, t, k in zip(industries, T_par, K_var_s)] _cache['se_a'] = float(std_err) # Manufacturing subsample mfg_codes = [i for i in industries if i[0] == '3'] if len(mfg_codes) >= 8: mfg_idx = [industries.index(c) for c in mfg_codes] mfg_T = log_T[mfg_idx] mfg_Kv = log_Kv[mfg_idx] _, sl_m, se_m, r2_m, _ = _ols_hc(mfg_T, mfg_Kv) rs_m, ps_m = scipy_stats.spearmanr(mfg_Kv, mfg_T) ci_lo = sl_m - 1.96 * se_m ci_hi = sl_m + 1.96 * se_m if rs_m > 0.3 and ps_m < 0.10 and ci_lo <= 1.0 <= ci_hi: verd_m = 'CONSISTENT' elif rs_m > 0 and ps_m < 0.10: verd_m = 'AMBIGUOUS' else: verd_m = 'INCONSISTENT' results.append(log_result( name='ces_isotropy_bea_mfg_subsample', stat=float(sl_m), p=float(ps_m), verd=verd_m, description="Tests isotropy on the manufacturing subsample (NAICS 3xx) where Oberfield-Raval sigma " "calibrations are most reliable. Same HC0-robust OLS specification as the full BEA I-O test, restricted to manufacturing industries.", details={'n': len(mfg_codes), 'beta': float(sl_m), 'se': float(se_m), 'spearman_r': float(rs_m)}, )) print(f" A-mfg: N={len(mfg_codes)}, beta={sl_m:.3f}, " f"Spearman r={rs_m:.3f} -> {verd_m}") return results # ══════════════════════════════════════════════════════════════════════════ # APPROACH B: 4-digit Euler-VRI (Census ECNMATFUEL + FRED IP) # ══════════════════════════════════════════════════════════════════════════ def _approach_b(): """4-digit NAICS Euler-VRI test.""" results = [] # Load Census ECNMATFUEL 4-digit Var(s) cache cache_4dig = COMPILED_DIR / 'ecnmatfuel_4digit_var_s.csv' if not cache_4dig.exists(): return results io_detail = {} with open(cache_4dig) as f: reader = csv.DictReader(f) for row in reader: io_detail[row['naics4']] = { 'J': int(row['J']), 'var_s': float(row['var_s']), 'H': float(row['H']), } if not io_detail: return results # Load FRED 4-digit IP series HP_LAMBDA_MONTHLY = 129_600 def _compute_cv2(values): arr = np.asarray(values, dtype=float) arr = arr[~np.isnan(arr) & (arr > 0)] if len(arr) < 60: return None cycle, trend = hp_filter(np.log(arr), HP_LAMBDA_MONTHLY) return float(np.var(cycle)) fred_results = {} # 4-digit NAICS candidates naics4_list = sorted(set(k[:4] for k in io_detail.keys())) for naics4 in naics4_list: series_id = f'IPG{naics4}N' path = FRED_DIR / f'{series_id}.csv' if not path.exists(): continue try: df = pd.read_csv(path, header=0) df.columns = ['date', 'value'] if len(df.columns) >= 2 else df.columns df['value'] = pd.to_numeric(df.iloc[:, -1], errors='coerce') vals = df['value'].dropna().values if len(vals) < 60: continue cv2 = _compute_cv2(vals) if cv2 is not None: fred_results[naics4] = cv2 except Exception: continue # Match records = [] for naics4, cv2 in fred_results.items(): parent3 = naics4[:3] sigma = OR_SIGMA_3DIG.get(parent3) if sigma is None: continue rho = rho_from_sigma(sigma) io_info = io_detail.get(naics4) if io_info is None or io_info['var_s'] <= 0 or io_info['J'] < 2: continue K_n = curvature_K(rho, io_info['J']) rhs = K_n * io_info['var_s'] if rhs <= 0 or cv2 <= 0: continue records.append({ 'naics4': naics4, 'sigma': sigma, 'rho': rho, 'J': io_info['J'], 'var_s': io_info['var_s'], 'K': K_n, 'rhs': rhs, 'lhs': cv2, }) if len(records) < 5 or not HAS_SCIPY: return results lhs_arr = np.array([r['lhs'] for r in records]) rhs_arr = np.array([r['rhs'] for r in records]) log_lhs = np.log(lhs_arr) log_rhs = np.log(rhs_arr) alpha_hat, beta_hat, se_beta, r2, resid = _ols_hc(log_lhs, log_rhs) r_s, p_s = scipy_stats.spearmanr(log_rhs, log_lhs) n = len(records) beta_ok = abs(beta_hat - 1.0) <= 1.5 * se_beta rank_ok = p_s < 0.05 if beta_ok and rank_ok: verd = 'CONSISTENT' elif rank_ok: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' results.append(log_result( name='ces_isotropy_euler_vri_4digit', stat=float(beta_hat), p=float(p_s), verd=verd, description="Tests isotropy at 4-digit NAICS using materials-only input shares from Census ECNMATFUEL " "and FRED monthly IP for output volatility. HC0-robust OLS of log(CV^2) on log(K*Var(s)) with Spearman rank check.", details={'n': n, 'beta': float(beta_hat), 'se': float(se_beta), 'spearman_r': float(r_s), 'r2': float(r2)}, )) print(f" B: N={n}, beta={beta_hat:.3f} (SE={se_beta:.3f}), " f"Spearman r={r_s:.3f} (p={p_s:.4f}) -> {verd}") return results # ══════════════════════════════════════════════════════════════════════════ # APPROACH C: Overidentification (BEA I-O 2022 + BEA GDP quarterly) # ══════════════════════════════════════════════════════════════════════════ def _approach_c(): """Overidentification using BEA I-O 2022 and quarterly GDP. Tests isotropy with a different I-O year (2022 vs 2017 in Approaches A/D) as a temporal stability check. Uses overidentification: tests both alpha=0 and beta=1. """ results = [] # Load BEA I-O 2022 — filter to intermediate commodity inputs only io_path = COMPILED_DIR / 'bea_io_use_2022.csv' if not io_path.exists(): return results # Non-commodity row prefixes to exclude (taxes, value-added, totals, govt) _EXCLUDE_PREFIXES = ('T', 'V', 'F', 'G', 'U', 'O') sector_inputs = {} with open(io_path) as f: reader = csv.DictReader(f) for row in reader: if row.get('ColType', '').strip() != 'Industry': continue col_code = row.get('ColCode', '').strip() row_code = row.get('RowCode', '').strip() # Filter: keep only intermediate commodity inputs if not row_code or any(row_code.startswith(p) for p in _EXCLUDE_PREFIXES): continue try: val = float(str(row.get('DataValue', '0')).replace(',', '')) except (ValueError, TypeError): val = 0.0 if val <= 0: continue if col_code not in sector_inputs: sector_inputs[col_code] = {} sector_inputs[col_code][row_code] = val if not sector_inputs: return results # Load quarterly VA (same source as Approach A for consistency) va_path = BEA_GDP_DIR / 'RealValueAdded_10_quarterly.csv' if not va_path.exists(): return results va = pd.read_csv(va_path) va['val'] = pd.to_numeric(va.DataValue, errors='coerce') va = va.dropna(subset=['val']) va = va[va.val > 0] va_data = {} for ind in va.Industry.unique(): ind_data = va[va.Industry == ind].sort_values(['Year', 'Quarter']) y = ind_data.val.values if len(y) < 20: continue log_y = np.log(y) cycle, trend = hp_filter(log_y, lamb=1600) T_par = np.var(cycle) if T_par > 0: va_data[ind] = T_par # Match all sectors with sigma calibrations (not just OR_SECTORS) matched = [] for bea_code, sigma in SIGMA_MAP_BEA.items(): rho = rho_from_sigma(sigma) T_par = va_data.get(bea_code) if T_par is None: continue bea_inputs = sector_inputs.get(bea_code) if bea_inputs is None: continue vals = np.array(list(bea_inputs.values()), dtype=float) if len(vals) < 3: continue total = vals.sum() if total <= 0: continue shares = vals / total J_n = len(shares) K_n = curvature_K(rho, J_n) var_s = np.mean((shares - 1.0 / J_n) ** 2) K_var_s = K_n * var_s if K_var_s <= 0: continue matched.append({ 'code': bea_code, 'sigma': sigma, 'rho': rho, 'T_par': T_par, 'K_var_s': K_var_s, 'J': J_n, }) N = len(matched) if N < 10 or not HAS_SCIPY: return results log_lhs = np.array([np.log(r['T_par']) for r in matched]) log_rhs = np.array([np.log(r['K_var_s']) for r in matched]) alpha, beta, se_beta, r2, resid = _ols_hc(log_lhs, log_rhs) spearman_r, spearman_p = scipy_stats.spearmanr(log_rhs, log_lhs) # Overidentification: test alpha=0 via HC0-robust SE XtX_inv = np.linalg.inv(np.column_stack([np.ones(N), log_rhs]).T @ np.column_stack([np.ones(N), log_rhs])) S = np.column_stack([np.ones(N), log_rhs]).T @ np.diag(resid ** 2) @ np.column_stack([np.ones(N), log_rhs]) cov = XtX_inv @ S @ XtX_inv se_alpha = np.sqrt(cov[0, 0]) t_alpha = alpha / se_alpha if se_alpha > 0 else np.inf t_beta = (beta - 1.0) / se_beta if se_beta > 0 else np.inf alpha_ok = abs(t_alpha) < 2.0 beta_ok = abs(t_beta) < 2.0 spearman_ok = spearman_r > 0 and spearman_p < 0.10 if alpha_ok and beta_ok and spearman_ok: verd = 'CONSISTENT' elif spearman_ok or beta_ok: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' results.append(log_result( name='ces_isotropy_overidentification', stat=float(beta), p=float(spearman_p), verd=verd, description="Overidentification test using BEA I-O 2022 (different year from Approaches A/D). " "Tests alpha=0 (no intercept bias) and beta=1 (unit slope) with HC0-robust SEs. " "Filters to intermediate commodity inputs only (excludes taxes, VA, final demand).", details={'n': N, 'alpha': float(alpha), 'beta': float(beta), 'se_alpha': float(se_alpha), 'se_beta': float(se_beta), 'spearman_r': float(spearman_r), 'r2': float(r2), 't_alpha': float(t_alpha), 't_beta': float(t_beta)}, )) print(f" C: N={N}, alpha={alpha:.3f}, beta={beta:.3f} (SE={se_beta:.3f}), " f"Spearman r={spearman_r:.3f} (p={spearman_p:.4f}) -> {verd}") return results # ══════════════════════════════════════════════════════════════════════════ # APPROACH D: BEA Benchmark 2017 detailed Use Table (402 industries) # ══════════════════════════════════════════════════════════════════════════ # Map detailed NAICS 4-digit prefix to GDP summary code _PREFIX_TO_GDP = { # Agriculture / forestry '1111': '111CA', '1112': '111CA', '1119': '111CA', '1121': '111CA', '1122': '111CA', '1123': '111CA', '1124': '111CA', '1125': '111CA', '1129': '111CA', '112A': '111CA', '1130': '113FF', '1140': '113FF', '1150': '113FF', # Mining '2110': '211', '2121': '212', '2122': '212', '2123': '212', '2131': '213', # Utilities '2211': '22', '2212': '22', '2213': '22', # Construction '2303': '23', '2332': '23', '2333': '23', '2334': '23', '2361': '23', '2362': '23', '2371': '23', '2372': '23', '2373': '23', '2379': '23', '2381': '23', '2382': '23', '2383': '23', '2389': '23', # Wholesale trade '4231': '42', '4232': '42', '4233': '42', '4234': '42', '4235': '42', '4236': '42', '4237': '42', '4238': '42', '4239': '42', '4241': '42', '4242': '42', '4243': '42', '4244': '42', '4245': '42', '4246': '42', '4247': '42', '4248': '42', '4249': '42', } # 3-digit NAICS to GDP for manufacturing composites _N3_TO_GDP = { '311': '311FT', '312': '311FT', '313': '313TT', '314': '313TT', '315': '315AL', '316': '315AL', } def _detail_to_gdp(code): """Map a detailed 6-digit BEA code to its GDP summary parent.""" code = str(code) if not code[0].isdigit(): return None # Try 4-digit prefix first (non-manufacturing) p4 = code[:4] if p4 in _PREFIX_TO_GDP: return _PREFIX_TO_GDP[p4] # Manufacturing: 3-digit composites n3 = code[:3] if n3 in _N3_TO_GDP: return _N3_TO_GDP[n3] # Motor vehicles vs other transport if n3 == '336': return '3361MV' if p4 in ('3361', '3362', '3363') else '3364OT' # Other manufacturing: 3-digit → GDP code if n3.startswith('3'): return n3 return None def _approach_d(): """BEA 2017 Benchmark detailed Use Table (402 industries × 385 commodities). Uses the detailed commodity breakdown to compute Var(s_nj) with full input bundle resolution, aggregated to GDP summary codes for VA matching. This gives 49:1 cross-sectional variation in Var(s) vs 1.4:1 in the summary table. """ results = [] detail_path = COMPILED_DIR / 'bea_use_detail_2017.csv' if not detail_path.exists(): return results use = pd.read_csv(detail_path) # Aggregate detailed purchases to GDP summary level gdp_inputs = {} for _, row in use.iterrows(): gdp = _detail_to_gdp(row['industry_code']) if gdp is None: continue if gdp not in gdp_inputs: gdp_inputs[gdp] = {} com = row['commodity_code'] gdp_inputs[gdp][com] = gdp_inputs[gdp].get(com, 0) + row['value'] # Load quarterly VA (same as approach A) va_path = BEA_GDP_DIR / 'RealValueAdded_10_quarterly.csv' if not va_path.exists(): return results va = pd.read_csv(va_path) va['val'] = pd.to_numeric(va.DataValue, errors='coerce') va = va[va.val > 0].dropna(subset=['val']) va_data = {} for ind in va.Industry.unique(): ind_data = va[va.Industry == ind].sort_values(['Year', 'Quarter']) y = ind_data.val.values if len(y) < 20: continue log_y = np.log(y) cycle, trend = hp_filter(log_y, lamb=1600) T_par = np.var(cycle) if T_par > 0: va_data[ind] = T_par # Match GDP industries with sigma estimates records = [] for gdp_code, inputs in gdp_inputs.items(): if gdp_code not in SIGMA_MAP_BEA and gdp_code not in va_data: continue sigma = SIGMA_MAP_BEA.get(gdp_code) T_par = va_data.get(gdp_code) if sigma is None or T_par is None: continue vals = np.array([v for v in inputs.values() if v > 0]) if len(vals) < 3: continue total = vals.sum() shares = vals / total J = len(shares) rho = rho_from_sigma(sigma) K = curvature_K(rho, J) var_s = np.mean((shares - 1.0 / J) ** 2) K_var_s = K * var_s if K_var_s <= 0: continue records.append({ 'gdp_code': gdp_code, 'sigma': sigma, 'rho': rho, 'J': J, 'var_s': var_s, 'K': K, 'K_var_s': K_var_s, 'T_par': T_par, }) N = len(records) if N < 8 or not HAS_SCIPY: return results log_T = np.array([np.log(r['T_par']) for r in records]) log_Kv = np.array([np.log(r['K_var_s']) for r in records]) alpha_hat, slope, std_err, r2, resid = _ols_hc(log_T, log_Kv) r_spearman, p_spearman = scipy_stats.spearmanr(log_Kv, log_T) beta_lo = slope - 1.96 * std_err beta_hi = slope + 1.96 * std_err if r_spearman > 0.5 and p_spearman < 0.01 and beta_lo <= 1.0 <= beta_hi: verd = 'CONSISTENT' elif r_spearman > 0.3 and p_spearman < 0.05: verd = 'AMBIGUOUS' elif r_spearman > 0 and p_spearman < 0.10: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' desc = ( "Tests isotropy T_parallel = K * Var(s) using BEA 2017 Benchmark detailed " "Use Table (402 industries × 385 commodities) aggregated to GDP summary codes, " "with quarterly real VA for output volatility. HC0-robust standard errors. " "The detailed table provides 49:1 cross-sectional variation in Var(s) vs 1.4:1 in the summary table." ) results.append(log_result( name='ces_isotropy_benchmark_detail', stat=float(slope), p=float(p_spearman), verd=verd, description=desc, details={ 'n': N, 'beta': float(slope), 'se': float(std_err), 'spearman_r': float(r_spearman), 'r2': float(r2), 'beta_ci': [float(beta_lo), float(beta_hi)], 'var_s_range': [float(min(r['var_s'] for r in records)), float(max(r['var_s'] for r in records))], }, )) print(f" D: N={N}, beta={slope:.3f} (SE={std_err:.3f}), " f"Spearman r={r_spearman:.3f} (p={p_spearman:.4f}) -> {verd}") # Manufacturing subsample mfg = [r for r in records if r['gdp_code'][0] == '3'] if len(mfg) >= 8: mfg_T = np.array([np.log(r['T_par']) for r in mfg]) mfg_Kv = np.array([np.log(r['K_var_s']) for r in mfg]) _, sl_m, se_m, r2_m, _ = _ols_hc(mfg_T, mfg_Kv) rs_m, ps_m = scipy_stats.spearmanr(mfg_Kv, mfg_T) ci_lo = sl_m - 1.96 * se_m ci_hi = sl_m + 1.96 * se_m if rs_m > 0.3 and ps_m < 0.10 and ci_lo <= 1.0 <= ci_hi: verd_m = 'CONSISTENT' elif rs_m > 0 and ps_m < 0.10: verd_m = 'AMBIGUOUS' else: verd_m = 'INCONSISTENT' results.append(log_result( name='ces_isotropy_benchmark_mfg', stat=float(sl_m), p=float(ps_m), verd=verd_m, description="Manufacturing subsample of the Benchmark detailed Use Table test. " "HC0-robust OLS. Oberfield-Raval sigma calibrations are most reliable for NAICS 31-33.", details={'n': len(mfg), 'beta': float(sl_m), 'se': float(se_m), 'spearman_r': float(rs_m), 'beta_ci': [float(ci_lo), float(ci_hi)]}, )) print(f" D-mfg: N={len(mfg)}, beta={sl_m:.3f}, " f"Spearman r={rs_m:.3f} -> {verd_m}") _cache['records_d'] = records _cache['se_d'] = float(std_err) return results # ── Main ────────────────────────────────────────────────────────────────── _cache = {} def main(): """Run CES isotropy tests. Returns list of result dicts.""" print("CES Isotropy Overidentification Tests (U1)") print("=" * 60) all_results = [] # Approach A: BEA full I-O print("\n[A] BEA full I-O (66 industries)...") res_a = _approach_a() all_results.extend(res_a) if not res_a: print(" No results (data not available)") all_results.append(log_result( name='ces_isotropy_bea_full_io', stat=np.nan, p=np.nan, verd='SKIPPED', description="Tests isotropy prediction T_parallel = T_perp across 66 BEA industries. " "OLS of log(HP-cycle variance) on log(K_eff * Var(input shares)) with H0: beta=1. " "Uses BEA Use Table input shares and quarterly real value added.", details={'reason': 'BEA I-O or GDP data not found'}, )) # Approach B: 4-digit Euler-VRI print("\n[B] 4-digit NAICS Euler-VRI...") res_b = _approach_b() all_results.extend(res_b) if not res_b: print(" No results (data not available)") all_results.append(log_result( name='ces_isotropy_euler_vri_4digit', stat=np.nan, p=np.nan, verd='SKIPPED', description="Tests isotropy at 4-digit NAICS using materials-only input shares from Census ECNMATFUEL " "and FRED monthly IP for output volatility. HC0-robust OLS of log(CV^2) on log(K*Var(s)) with Spearman rank check.", details={'reason': 'Census ECNMATFUEL or FRED 4-digit data not found'}, )) # Approach C: Overidentification print("\n[C] Overidentification (BEA I-O 2022)...") res_c = _approach_c() all_results.extend(res_c) if not res_c: print(" No results (data not available)") all_results.append(log_result( name='ces_isotropy_overidentification', stat=np.nan, p=np.nan, verd='SKIPPED', description="Overidentification test using BEA I-O 2022 input shares and quarterly GDP by industry " "for 14 Oberfield-Raval sectors. Tests both alpha=0 (no intercept bias) and beta=1 (unit slope) " "via t-tests, plus Spearman rank correlation.", details={'reason': 'BEA I-O 2022 or GDP quarterly data not found'}, )) # Approach D: Benchmark detailed Use Table print("\n[D] Benchmark detailed Use Table (402 industries)...") res_d = _approach_d() all_results.extend(res_d) if not res_d: print(" No results (data not available)") all_results.append(log_result( name='ces_isotropy_benchmark_detail', stat=np.nan, p=np.nan, verd='SKIPPED', description="Tests isotropy using BEA 2017 Benchmark detailed Use Table " "(402 industries × 385 commodities). Run download_data.py to build cache.", details={'reason': 'bea_use_detail_2017.csv not found'}, )) print(f"\n Results: {len(all_results)} tests") return all_results def figures(): """Generate scatter plots for CES isotropy tests.""" if not _cache: main() paths = [] import matplotlib.pyplot as plt # --- Approach D: Benchmark detailed Use Table scatter --- records = _cache.get('records_d') if records and len(records) >= 8: log_T = np.array([np.log(r['T_par']) for r in records]) log_Kv = np.array([np.log(r['K_var_s']) for r in records]) labels = [r['gdp_code'] for r in records] is_mfg = [r['gdp_code'][0] == '3' for r in records] from scipy import stats as sp_stats slope, intercept, _, _, _ = sp_stats.linregress(log_Kv, log_T) r_s, p_s = sp_stats.spearmanr(log_Kv, log_T) fig, axes = setup_subplots(1, 2, width=14, height=6) # Panel A: Full scatter with regression line and 1:1 line ax = axes[0] colors = ['#2166ac' if m else '#b2182b' for m in is_mfg] ax.scatter(log_Kv, log_T, c=colors, s=40, alpha=0.7, edgecolors='k', lw=0.3, zorder=3) for i, lab in enumerate(labels): ax.annotate(lab, (log_Kv[i], log_T[i]), fontsize=5, alpha=0.7, xytext=(3, 3), textcoords='offset points') xs = np.linspace(log_Kv.min() - 0.5, log_Kv.max() + 0.5, 100) ax.plot(xs, intercept + slope * xs, 'k-', lw=1.5, label=f'OLS: $\\beta$={slope:.2f} (SE={_cache.get("se_d", 0):.2f})') ax.plot(xs, xs, 'r--', lw=1, alpha=0.5, label='Isotropy ($\\beta$=1)') ax.set_xlabel('log(K · Var(s)) [transverse proxy]') ax.set_ylabel('log(T$_{\\parallel}$) [HP-cycle variance]') ax.set_title(f'(a) CES Isotropy: Benchmark Detail (N={len(records)})') ax.legend(fontsize=8) ax.grid(True, alpha=0.2) # Add Spearman annotation ax.text(0.05, 0.95, f'Spearman r={r_s:.3f}\np={p_s:.4f}', transform=ax.transAxes, fontsize=8, va='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) # Panel B: T_par / K*Var(s) ratio by sector ax = axes[1] ratios = [r['T_par'] / r['K_var_s'] for r in records] sorted_idx = np.argsort(ratios) sorted_ratios = [ratios[i] for i in sorted_idx] sorted_labels = [labels[i] for i in sorted_idx] sorted_colors = [colors[i] for i in sorted_idx] ax.barh(range(len(sorted_ratios)), sorted_ratios, color=sorted_colors, alpha=0.7, edgecolor='k', lw=0.3) ax.set_yticks(range(len(sorted_labels))) ax.set_yticklabels(sorted_labels, fontsize=6) ax.axvline(1.0, color='red', ls='--', lw=1, alpha=0.5, label='Isotropy (ratio=1)') ax.axvline(np.median(ratios), color='grey', ls=':', lw=1, label=f'Median={np.median(ratios):.1f}') ax.set_xlabel('$T_{\\parallel}$ / ($K$ · Var(s))') ax.set_title('(b) T-Anisotropy Ratio by Sector') ax.set_xscale('log') ax.legend(fontsize=7, loc='lower right') ax.grid(True, alpha=0.2, axis='x') fig.tight_layout() paths.extend(save_figure(fig, 'ces_isotropy_benchmark', subfolder='foundations')) return paths if __name__ == '__main__': results = main() for r in results: print(f" {r['name']:50s} {r['verdict']:15s} stat={r['statistic']}")

CES Isotropy Euler VRI 4digit

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources