#!/usr/bin/env python3 """ Test: Directed Technical Change under CES ========================================== Tests whether technical change is directed toward the relatively scarce factor, as predicted by the CES model when sigma < 2 (Acemoglu 2002). Five sub-tests: 1. Time-series rho estimation from national data 2. Cross-state rho estimation from BEA state data 3. Consistency (Wald test: time-series rho ~ cross-state rho) 4. sigma < 2 threshold (bootstrap confidence interval) 5. Capital deepening direction (K/L ratio increasing over time) Data: BEA state GDP, compensation, GOS, employment; national GDI and capital stock. Review Notes (2026-03-12): - Theory match: CORRECT — subtests 1,3,4 match Lean theorems exactly; subtest 5 now tests the DTC directional prediction (Spearman(dlog(K/L), dlog(w_K/w_L)) < 0); subtest 2 is known-weak proxy (GDP/worker != K/L) - I/O audit: CORRECT for subtests 1,3,4,5; WEAK for subtest 2 - Data-signal alignment: STRONG for time-series and DTC direction (BEA NIPA direct), WEAK for cross-state (no state-level K data) - Code quality: CLEAN — real data figures, proper _cache, correct figure registry - Confidence: HIGH — 3 informative CONSISTENT subtests (ts_rho, sigma<2, DTC direction), 1 expected AMBIGUOUS (cross-state), 1 low-power CONSISTENT (consistency) """ import numpy as np import pandas as pd from scipy import stats import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) from common.data import load_compiled from common.ces import sigma_from_rho, rho_from_sigma from common.stats import verdict, log_result from common.plotting import setup_figure, setup_subplots, save_figure def load_state_data(): """Load BEA state-level data.""" gdp = load_compiled('state_gdp_2022') comp = load_compiled('state_compensation_2022') gos = load_compiled('state_gos_2022') emp = load_compiled('state_employment_2022') return gdp, comp, gos, emp def load_national_ts(): """Load national time series.""" gdi = load_compiled('us_gdi_timeseries') cap = load_compiled('us_capital_stock_timeseries') emp = load_compiled('us_employment_timeseries') return gdi, cap, emp def merge_state_data(gdp, comp, gos, emp): """Merge state datasets.""" state_col = None for col in ['state', 'State', 'GeoName', 'geoname', 'GeoFips']: if col in gdp.columns: state_col = col break if state_col is None: state_col = gdp.columns[0] def val_col(df): numeric = df.select_dtypes(include=[np.number]).columns return numeric[-1] if len(numeric) > 0 else df.columns[-1] merged = pd.DataFrame() merged['state'] = gdp[state_col] merged['gdp'] = pd.to_numeric(gdp[val_col(gdp)], errors='coerce') merged['compensation'] = pd.to_numeric(comp[val_col(comp)], errors='coerce') merged['gos'] = pd.to_numeric(gos[val_col(gos)], errors='coerce') merged['employment'] = pd.to_numeric(emp[val_col(emp)], errors='coerce') merged = merged.dropna() merged['labor_share'] = merged['compensation'] / merged['gdp'] merged['capital_share'] = merged['gos'] / merged['gdp'] return merged def extract_ts(df): """Extract year and value arrays from a dataframe.""" year_col = [c for c in df.columns if 'year' in c.lower()] if year_col: years = df[year_col[0]].values val_col = [c for c in df.columns if c != year_col[0]] vals = pd.to_numeric(df[val_col[0]], errors='coerce').values else: years = np.arange(len(df)) vals = pd.to_numeric(df.iloc[:, -1], errors='coerce').values return years, vals # ── Sub-tests ───────────────────────────────────────────────────── def test_time_series_rho(gdi, cap, emp): """Test 1: Estimate rho from national time series. Levels regression: log(s_L/s_K) = c + rho * log(L/K) where s_L = Compensation/GDI, s_K = (NOS+CFC)/GDI. """ # Merge on Year df = gdi.merge(cap, on='Year').merge(emp, on='Year') # Factor shares df['s_L'] = df['Compensation'] / df['GDI'] df['s_K'] = (df['NetOperatingSurplus'] + df['CFC']) / df['GDI'] df['LK_ratio'] = df['TotalEmployment'] / df['NonresidentialCapitalStock'] df = df[(df['s_L'] > 0) & (df['s_K'] > 0) & (df['LK_ratio'] > 0)].copy() if len(df) < 10: return log_result('directed_technical_change_ts_rho', 0.0, 1.0, 'SKIPPED', description="Estimates rho from US national time series using OLS on log(s_L/s_K) ~ log(L/K). Tests whether the CES elasticity of substitution sigma = 1/(1-rho) falls in a reasonable range (0, 3).") log_share_ratio = np.log(df['s_L'].values / df['s_K'].values) log_LK = np.log(df['LK_ratio'].values) slope, intercept, r_val, p_val, se = stats.linregress(log_LK, log_share_ratio) rho_est = slope sigma_est = sigma_from_rho(rho_est) sign_correct = p_val < 0.05 and 0.0 < sigma_est < 3.0 verd = verdict(p_val, sign_correct) _cache['ts_x'] = log_LK _cache['ts_y'] = log_share_ratio _cache['ts_slope'] = slope _cache['ts_intercept'] = intercept return log_result( 'directed_technical_change_ts_rho', rho_est, p_val, verd, {'rho': rho_est, 'sigma': sigma_est, 'r_squared': r_val ** 2, 'slope': slope, 'se': se}, description="Estimates rho from US national time series using OLS on log(s_L/s_K) ~ log(L/K). Tests whether the CES elasticity of substitution sigma = 1/(1-rho) falls in a reasonable range (0, 3)." ) def test_cross_state_rho(state_df): """Test 2: Estimate rho from cross-state variation.""" if state_df is None or len(state_df) < 10: return log_result('directed_technical_change_xs_rho', 0.0, 1.0, 'SKIPPED', description="Estimates rho from cross-state variation using OLS on log(s_L/s_K) ~ log(GDP/worker). Uses BEA state GDP, compensation, and GOS data to test whether sigma falls in (0, 3).") df = state_df.copy() df = df[(df['labor_share'] > 0) & (df['capital_share'] > 0)] # Use GDP/employment as independent proxy for capital intensity # (NOT gos/compensation which is algebraically -1 * LHS) df = df[df['employment'] > 0] df['gdp_per_worker'] = df['gdp'] / df['employment'] df = df[df['gdp_per_worker'] > 0] if len(df) < 10: return log_result('directed_technical_change_xs_rho', 0.0, 1.0, 'SKIPPED', description="Estimates rho from cross-state variation using OLS on log(s_L/s_K) ~ log(GDP/worker). Uses BEA state GDP, compensation, and GOS data to test whether sigma falls in (0, 3).") log_share_ratio = np.log(df['labor_share'] / df['capital_share']).values log_kl = np.log(df['gdp_per_worker']).values slope, intercept, r_val, p_val, se = stats.linregress(log_kl, log_share_ratio) rho_est = slope sigma_est = sigma_from_rho(rho_est) sign_correct = 0.0 < sigma_est < 3.0 verd = verdict(p_val, sign_correct) _cache['xs_x'] = log_kl _cache['xs_y'] = log_share_ratio _cache['xs_slope'] = slope _cache['xs_intercept'] = intercept return log_result( 'directed_technical_change_xs_rho', rho_est, p_val, verd, {'rho': rho_est, 'sigma': sigma_est, 'slope': slope, 'se': se, 'r_squared': r_val ** 2, 'n_states': len(df)}, description="Estimates rho from cross-state variation using OLS on log(s_L/s_K) ~ log(GDP/worker). Uses BEA state GDP, compensation, and GOS data to test whether sigma falls in (0, 3)." ) def test_consistency(ts_result, xs_result): """Test 3: Wald test — time-series rho ~ cross-state rho.""" rho_ts = ts_result.get('rho', None) rho_xs = xs_result.get('rho', None) se_ts = ts_result.get('se', None) se_xs = xs_result.get('se', None) if any(v is None for v in [rho_ts, rho_xs, se_ts, se_xs]): return log_result('directed_technical_change_consistency', 0.0, 1.0, 'SKIPPED', description="Wald test for equality of time-series and cross-state rho estimates. Consistent (high p-value) means both identification strategies yield the same CES parameter.") if se_ts <= 0 or se_xs <= 0: return log_result('directed_technical_change_consistency', 0.0, 1.0, 'SKIPPED', description="Wald test for equality of time-series and cross-state rho estimates. Consistent (high p-value) means both identification strategies yield the same CES parameter.") # Wald test: (rho_ts - rho_xs) / sqrt(se_ts^2 + se_xs^2) diff = rho_ts - rho_xs se_diff = np.sqrt(se_ts ** 2 + se_xs ** 2) z_stat = diff / se_diff p_val = 2.0 * stats.norm.sf(abs(z_stat)) # Consistent if cannot reject equality (high p-value) sign_correct = p_val > 0.05 verd = 'CONSISTENT' if p_val > 0.10 else ( 'AMBIGUOUS' if p_val > 0.01 else 'INCONSISTENT') return log_result( 'directed_technical_change_consistency', z_stat, p_val, verd, {'rho_ts': rho_ts, 'rho_xs': rho_xs, 'diff': diff}, description="Wald test for equality of time-series and cross-state rho estimates. Computes z = (rho_ts - rho_xs) / sqrt(se_ts^2 + se_xs^2). Consistent (high p-value) means both identification strategies yield the same CES parameter." ) def test_sigma_threshold(): """Test 4: sigma < 2 (bootstrap CI from time-series data). Uses the time-series identification (log(s_L/s_K) ~ log(L/K)) which has direct K and L data, unlike the cross-state proxy (GDP/worker). Lean: dtcPremium_threshold_sigma_2. """ skip_desc = ( "Bootstrap test that sigma < 2, the Acemoglu (2002) threshold for " "directed technical change toward the scarce factor. Uses 10000 bootstrap " "resamples of time-series rho from log(s_L/s_K) ~ log(L/K) on BEA NIPA " "data; p-value is fraction of bootstraps with sigma >= 2." ) ts_x = _cache.get('ts_x') ts_y = _cache.get('ts_y') if ts_x is None or len(ts_x) < 10: return log_result('directed_technical_change_sigma_threshold', 0.0, 1.0, 'SKIPPED', description=skip_desc) # Bootstrap rho estimates from time-series data n_boot = 10000 np.random.seed(42) sigma_boots = [] for _ in range(n_boot): idx = np.random.choice(len(ts_x), len(ts_x), replace=True) try: slope, _, _, _, _ = stats.linregress(ts_x[idx], ts_y[idx]) sigma_b = sigma_from_rho(slope) if 0 < sigma_b < 20: sigma_boots.append(sigma_b) except Exception: continue if len(sigma_boots) < 100: return log_result('directed_technical_change_sigma_threshold', 0.0, 1.0, 'SKIPPED', description=skip_desc) sigma_boots = np.array(sigma_boots) ci_lower = np.percentile(sigma_boots, 2.5) ci_upper = np.percentile(sigma_boots, 97.5) median_sigma = np.median(sigma_boots) _cache['sigma_boots'] = sigma_boots # Test: upper CI < 2.0 sign_correct = ci_upper < 2.0 # Approximate p-value: fraction of bootstraps with sigma >= 2 p_val = np.mean(sigma_boots >= 2.0) verd = verdict(p_val if sign_correct else 1.0 - p_val, sign_correct) return log_result( 'directed_technical_change_sigma_threshold', median_sigma, p_val, verd, {'median_sigma': median_sigma, 'ci_lower': ci_lower, 'ci_upper': ci_upper, 'n_boot': len(sigma_boots)}, description=skip_desc ) def test_capital_deepening(gdi, cap, emp): """Test 5: Capital deepening lowers capital's relative price (DTC direction). When sigma in (1,2), the DTC premium exponent is sigma-2 < 0, so increasing K/L should LOWER w_K/w_L. Tests via Spearman(dlog(K/L), dlog(w_K/w_L)) < 0. Lean: dtcPremium_complement_regime. """ skip_desc = ( "Tests the DTC directional prediction: when sigma < 2, capital deepening " "lowers capital's relative price. Spearman correlation of dlog(K/L) vs " "dlog(w_K/w_L) in first differences. Uses BEA NIPA factor shares, capital " "stock, and employment. Prediction: negative correlation (price effect " "dominates market size effect)." ) df = gdi.merge(cap, on='Year').merge(emp, on='Year') df['s_L'] = df['Compensation'] / df['GDI'] df['s_K'] = (df['NetOperatingSurplus'] + df['CFC']) / df['GDI'] df['K'] = df['NonresidentialCapitalStock'] df['L'] = df['TotalEmployment'] df = df[(df['s_L'] > 0) & (df['s_K'] > 0) & (df['K'] > 0) & (df['L'] > 0)].copy() if len(df) < 10: return log_result('directed_technical_change_capital_deepening', 0.0, 1.0, 'SKIPPED', description=skip_desc) kl = (df['K'] / df['L']).values w_K = (df['s_K'] * df['GDI'] / df['K']).values w_L = (df['s_L'] * df['GDI'] / df['L']).values # First differences of logs d_log_kl = np.diff(np.log(kl)) d_log_wkwl = np.diff(np.log(w_K / w_L)) rho_s, p_s = stats.spearmanr(d_log_kl, d_log_wkwl) # Theory predicts negative correlation (sigma < 2 => price effect dominates) sign_correct = rho_s < 0 verd = verdict(p_s, sign_correct) _cache['deepening_d_log_kl'] = d_log_kl _cache['deepening_d_log_wkwl'] = d_log_wkwl return log_result( 'directed_technical_change_capital_deepening', rho_s, p_s, verd, {'spearman_r': rho_s, 'n_diffs': len(d_log_kl)}, description=skip_desc ) # ── Main ────────────────────────────────────────────────────────── _cache = {} def main(): """Run all directed technical change tests. Returns list of result dicts.""" results = [] state_df = None try: gdp, comp, gos, emp = load_state_data() state_df = merge_state_data(gdp, comp, gos, emp) except FileNotFoundError as e: print(f" State data loading failed: {e}") try: gdi, cap, emp_ts = load_national_ts() except FileNotFoundError as e: print(f" National data loading failed: {e}") gdi, cap, emp_ts = None, None, None # Test 1 if gdi is not None: r1 = test_time_series_rho(gdi, cap, emp_ts) else: r1 = log_result('directed_technical_change_ts_rho', 0.0, 1.0, 'SKIPPED', description="Estimates rho from US national time series using OLS on log(s_L/s_K) ~ log(L/K). Tests whether the CES elasticity of substitution sigma = 1/(1-rho) falls in a reasonable range (0, 3).") results.append(r1) # Test 2 r2 = test_cross_state_rho(state_df) results.append(r2) # Test 3: consistency results.append(test_consistency(r1, r2)) # Test 4: sigma < 2 (uses ts_x/ts_y from _cache, so must run after test 1) results.append(test_sigma_threshold()) # Test 5: capital deepening direction if gdi is not None: results.append(test_capital_deepening(gdi, cap, emp_ts)) else: results.append(log_result('directed_technical_change_capital_deepening', 0.0, 1.0, 'SKIPPED', description="Tests the DTC directional prediction: when sigma < 2, capital deepening " "lowers capital's relative price. Spearman correlation of dlog(K/L) vs " "dlog(w_K/w_L) in first differences. Uses BEA NIPA factor shares, capital " "stock, and employment. Prediction: negative correlation.")) _cache['results'] = results _cache['state_df'] = state_df return results def figures(): """Generate directed technical change figures with actual data plots.""" if not _cache: main() paths = [] results = _cache.get('results', []) if not results: return paths # Figure 1: Time-series rho scatter with regression line r1 = next((r for r in results if r['name'] == 'directed_technical_change_ts_rho'), None) ts_x = _cache.get('ts_x') ts_y = _cache.get('ts_y') if r1 and r1['verdict'] != 'SKIPPED' and ts_x is not None: fig, ax = setup_figure(8, 6) ax.scatter(ts_x, ts_y, alpha=0.7, s=30, zorder=3) x_line = np.linspace(ts_x.min(), ts_x.max(), 100) ax.plot(x_line, _cache['ts_slope'] * x_line + _cache['ts_intercept'], 'r-', linewidth=2, label=f'OLS: rho={r1["rho"]:.3f} (sigma={r1["sigma"]:.2f})') ax.set_xlabel('log(L/K)') ax.set_ylabel('log(s_L/s_K)') ax.set_title(f'Time-Series rho: R2={r1["r_squared"]:.3f}, p={r1["p_value"]:.1e}') ax.legend(fontsize=10) ax.grid(True, alpha=0.3) paths += save_figure(fig, 'dtc_timeseries_rho', subfolder='macro') # Figure 2: Cross-state rho scatter r2 = next((r for r in results if r['name'] == 'directed_technical_change_xs_rho'), None) xs_x = _cache.get('xs_x') xs_y = _cache.get('xs_y') if r2 and r2['verdict'] != 'SKIPPED' and xs_x is not None: fig, ax = setup_figure(8, 6) ax.scatter(xs_x, xs_y, alpha=0.6, s=30, zorder=3) x_line = np.linspace(xs_x.min(), xs_x.max(), 100) ax.plot(x_line, _cache['xs_slope'] * x_line + _cache['xs_intercept'], 'r-', linewidth=2, label=f'OLS: rho={r2["rho"]:.3f} (R2={r2["r_squared"]:.4f})') ax.set_xlabel('log(GDP/worker)') ax.set_ylabel('log(s_L/s_K)') ax.set_title(f'Cross-State rho: N={r2["n_states"]} states, p={r2["p_value"]:.3f}') ax.legend(fontsize=10) ax.grid(True, alpha=0.3) paths += save_figure(fig, 'dtc_crossstate_rho', subfolder='macro') # Figure 3: Sigma bootstrap histogram r4 = next((r for r in results if r['name'] == 'directed_technical_change_sigma_threshold'), None) sigma_boots = _cache.get('sigma_boots') if r4 and r4['verdict'] != 'SKIPPED' and sigma_boots is not None: fig, ax = setup_figure(8, 5) ax.hist(sigma_boots, bins=40, density=True, alpha=0.7, color='steelblue', edgecolor='white') ax.axvline(2, color='black', linewidth=2, linestyle='--', label='Acemoglu threshold: sigma=2') ax.axvline(r4['ci_lower'], color='red', linewidth=1.5, linestyle=':') ax.axvline(r4['ci_upper'], color='red', linewidth=1.5, linestyle=':', label=f'95% CI: [{r4["ci_lower"]:.3f}, {r4["ci_upper"]:.3f}]') ax.axvline(r4['median_sigma'], color='red', linewidth=2, label=f'Median: {r4["median_sigma"]:.3f}') ax.set_xlabel('sigma = 1/(1-rho)') ax.set_ylabel('Density') ax.set_title(f'Bootstrap: sigma < 2 (p={r4["p_value"]:.3f}, N={r4["n_boot"]})') ax.legend(fontsize=9) ax.grid(True, alpha=0.3) paths += save_figure(fig, 'dtc_sigma_bootstrap', subfolder='macro') # Figure 4: Capital deepening direction scatter d_kl = _cache.get('deepening_d_log_kl') d_wkwl = _cache.get('deepening_d_log_wkwl') r5 = next((r for r in results if r['name'] == 'directed_technical_change_capital_deepening'), None) if r5 and r5['verdict'] != 'SKIPPED' and d_kl is not None: fig, ax = setup_figure(8, 6) ax.scatter(d_kl, d_wkwl, alpha=0.7, s=40, zorder=3) # Add regression line from scipy.stats import linregress sl, ic, _, _, _ = linregress(d_kl, d_wkwl) x_line = np.linspace(d_kl.min(), d_kl.max(), 100) ax.plot(x_line, sl * x_line + ic, 'r-', linewidth=2) ax.axhline(0, color='gray', linewidth=0.5, linestyle='-') ax.axvline(0, color='gray', linewidth=0.5, linestyle='-') ax.set_xlabel(r'$\Delta\log(K/L)$') ax.set_ylabel(r'$\Delta\log(w_K/w_L)$') rho_s = r5.get('spearman_r', r5['statistic']) ax.set_title(f'DTC Direction: Spearman r={rho_s:.3f}, p={r5["p_value"]:.4f}\n' f'Negative = price effect dominates (sigma < 2)') ax.grid(True, alpha=0.3) paths += save_figure(fig, 'dtc_capital_deepening', subfolder='macro') return paths if __name__ == '__main__': results = main() for r in results: print(f" {r['name']}: {r['verdict']} " f"(stat={r['statistic']:.4f}, p={r['p_value']:.4f})")

se	0.0187
rho	0.1885
sigma	1.2323
slope	0.1885
p value	1.846e-10
r squared	0.7958
statistic	0.1885

Directed Technical Change Ts Rho

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps