#!/usr/bin/env python3 """ Test: V-Shape rho at Technology Wave Turning Points (U4) Paper 3 (Dynamics on the CES Landscape), Endogenous-rho section (Remark 16.2): During major technology waves, the aggregate complementarity parameter rho-bar should follow a V-shape -- declining during installation/frenzy, reaching a minimum near the Perez turning point (crisis), then rising during deployment. Five subtests probe this prediction and diagnose whether the theory is correct: A) rho_vshape_nberces: V-shape event study at ICT/Internet 2001 turning point using Oberfield-Raval time-varying rho_bar(t). The primary test. B) rho_vshape_micro: V-shape event study at 1983 (the actual rho_bar minimum between 1972-1992). Method validation: proves V-shapes CAN be detected in the data. If CONSISTENT, the absence at 2001 is a real finding. C) rho_breakpoint_2001: Segmented regression structural break test. Is the Perez turning point (2001) a significant structural break in rho_bar(t)? If NOT, the theory's timing prediction is wrong. D) rho_cycle_period: Autocorrelation analysis of detrended rho_bar(t). The theory predicts Kondratiev-scale cycles (~25-50 years). The data may show shorter Kuznets-scale cycles (~20 years). E) rho_vshape_inputs: BEA ICT capital excess growth (secondary proxy). Data: NBER-CES Manufacturing Industry Database (1958-2018): pay, vadd, cap, prodh BEA: FAAt301ESI capital stock by industry Review Notes (2026-03-11): Theory Match: PARTIAL — installation-phase decline confirmed, deployment recovery absent. L-shape rather than V-shape at 2001. The rho_bar data shows a ~22-year cycle with troughs at 1983 and 2012, neither aligned with the theory's Perez turning points (1973, 2001). I/O Audit: CORRECT — Oberfield-Raval estimation from NBER-CES is the methodology recommended by Paper 5 Section U4. Data-Signal Alignment: ADEQUATE — NBER-CES covers 1958-2018, ~460 manufacturing industries. 15-year rolling window yields estimates 1961-2013. Two complete/partial cycles observed. Code Quality: CLEAN Confidence: MEDIUM for method, HIGH for the L-shape finding. Possible theory issue: see Evaluation report. Date of review: 2026-03-11 """ _ALL_SUBTESTS = [ 'rho_vshape_nberces', # V-shape at 2001 (theory prediction) 'rho_vshape_micro', # V-shape at 1983 (method validation) 'rho_breakpoint_2001', # Structural break test 'rho_cycle_period', # Cycle timescale 'rho_vshape_inputs', # BEA capital intensity ] import warnings warnings.filterwarnings("ignore") 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_FA_DIR, DATA_DIR from common.stats import log_result from common.ces import 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: from common.data import load_nber_ces HAS_NBER = True except ImportError: HAS_NBER = False # ── Oberfield-Raval estimation parameters ──────────────────────────────── DIFF_HORIZON = 5 # Long-difference horizon (years) WINDOW = 15 # Rolling window size (years) for cross-sectional estimation # Event study windows (years relative to turning point) PRE_START = -8 # Start of pre-period relative to tp PRE_END = -2 # End of pre-period (exclusive of zone) ZONE_HALF = 2 # Half-width of turning-point zone POST_START = 2 # Start of post-period (exclusive of zone) POST_END = 8 # End of post-period relative to tp # For Approach B (BEA capital intensity) PEREZ_WAVES_ANNUAL = [ {'name': 'ICT/Internet (~2001)', 'label': 'internet', 'tp': 2001, 'install_start': 1997, 'deploy_end': 2010, 'canonical': True}, {'name': 'Finance/Digital (~2007)', 'label': 'gfc', 'tp': 2007, 'install_start': 2003, 'deploy_end': 2016, 'canonical': False}, ] NBER_RECESSIONS = [ (1973, 1975), (1980, 1980), (1981, 1982), (1990, 1991), (2001, 2001), (2007, 2009), ] # ── Oberfield-Raval rolling estimation ─────────────────────────────────── def _estimate_rho_bar_annual(nber): """ Estimate average rho across manufacturing industries per year using rolling cross-sectional Oberfield-Raval (2021) long-difference estimation. For each rolling window of WINDOW years centered on year t: 1. Compute h-year long differences for each industry: Delta_h ln(s_K/s_L) and Delta_h ln(K/L) 2. Pool all industry-year long differences within the window 3. Run cross-sectional OLS: Delta_h ln(s_K/s_L) = rho * Delta_h ln(K/L) + eps 4. sigma = 1/(1-rho), then rho_CES = 1 - 1/sigma Returns pd.Series indexed by year with rho_bar estimates. """ df = nber.copy() df = df.dropna(subset=['pay', 'vadd', 'cap', 'prodh']) df = df[(df['vadd'] > 0) & (df['pay'] > 0) & (df['cap'] > 0) & (df['prodh'] > 0)] df['s_L'] = np.clip(df['pay'] / df['vadd'], 0.01, 0.99) df['ln_sK_sL'] = np.log((1 - df['s_L']) / df['s_L']) df['ln_KL'] = np.log(df['cap'] / df['prodh']) df = df.sort_values(['naics', 'year']) diffs = [] for naics, grp in df.groupby('naics'): grp = grp.sort_values('year').set_index('year') for yr in grp.index: yr_h = yr + DIFF_HORIZON if yr_h in grp.index: d_share = grp.loc[yr_h, 'ln_sK_sL'] - grp.loc[yr, 'ln_sK_sL'] d_kl = grp.loc[yr_h, 'ln_KL'] - grp.loc[yr, 'ln_KL'] if np.isfinite(d_share) and np.isfinite(d_kl): mid_yr = (yr + yr_h) / 2.0 diffs.append({'mid_year': mid_yr, 'd_share': d_share, 'd_kl': d_kl}) if len(diffs) < 50: return pd.Series(dtype=float) diff_df = pd.DataFrame(diffs) all_years = sorted(diff_df['mid_year'].unique()) min_yr = int(np.ceil(min(all_years))) max_yr = int(np.floor(max(all_years))) half_w = WINDOW / 2.0 rho_estimates = {} for center_yr in range(min_yr, max_yr + 1): mask = (diff_df['mid_year'] >= center_yr - half_w) & (diff_df['mid_year'] < center_yr + half_w) window_data = diff_df[mask] if len(window_data) < 30: continue x = window_data['d_kl'].values y = window_data['d_share'].values ss_xx = np.sum((x - x.mean()) ** 2) if ss_xx < 1e-10: continue ss_xy = np.sum((x - x.mean()) * (y - y.mean())) rho_hat = ss_xy / ss_xx if abs(1 - rho_hat) < 1e-6: continue sigma_hat = 1.0 / (1.0 - rho_hat) if not (0.1 < sigma_hat < 10): continue rho_estimates[center_yr] = rho_from_sigma(sigma_hat) if not rho_estimates: return pd.Series(dtype=float) rho_bar = pd.Series(rho_estimates, name='rho_bar') rho_bar.index.name = 'year' return rho_bar # ── V-shape event study ────────────────────────────────────────────────── def _vshape_event_study(rho_bar, tp_year): """ V-shape event study on annual rho_bar series around a turning point. Returns dict with pre_mean, zone_mean, post_mean, vshape flag, depth. """ pre_vals = rho_bar[(rho_bar.index >= tp_year + PRE_START) & (rho_bar.index < tp_year + PRE_END)] zone_vals = rho_bar[(rho_bar.index >= tp_year - ZONE_HALF) & (rho_bar.index <= tp_year + ZONE_HALF)] post_vals = rho_bar[(rho_bar.index > tp_year + POST_START) & (rho_bar.index <= tp_year + POST_END)] if len(pre_vals) < 2 or len(zone_vals) < 2 or len(post_vals) < 2: return None pre_mean = pre_vals.mean() zone_mean = zone_vals.mean() post_mean = post_vals.mean() vshape = (zone_mean < pre_mean) and (zone_mean < post_mean) depth = min(pre_mean, post_mean) - zone_mean # Pre-slope pre_window = rho_bar[(rho_bar.index >= tp_year + PRE_START) & (rho_bar.index < tp_year)] slope_pre = np.nan if len(pre_window) >= 3: t = np.array(pre_window.index, dtype=float) t = t - t.mean() slope_pre = np.sum(t * (pre_window.values - pre_window.values.mean())) / np.sum(t ** 2) # Post-slope post_window = rho_bar[(rho_bar.index >= tp_year) & (rho_bar.index <= tp_year + POST_END)] slope_post = np.nan if len(post_window) >= 3: t = np.array(post_window.index, dtype=float) t = t - t.mean() slope_post = np.sum(t * (post_window.values - post_window.values.mean())) / np.sum(t ** 2) return { 'pre_mean': pre_mean, 'zone_mean': zone_mean, 'post_mean': post_mean, 'vshape': vshape, 'depth': depth, 'slope_pre': slope_pre, 'slope_post': slope_post, 'n_pre': len(pre_vals), 'n_zone': len(zone_vals), 'n_post': len(post_vals), } def _permutation_test_vshape(rho_bar, tp_year, n_perm=5000, seed=42): """ Permutation test: shuffle turning point across non-recession years. Returns (observed_depth, p_value). """ observed = _vshape_event_study(rho_bar, tp_year) if observed is None or not observed['vshape']: return (0.0, 1.0) obs_depth = observed['depth'] rec_years = set() for start, end in NBER_RECESSIONS: for y in range(start, end + 1): rec_years.add(y) min_yr = int(rho_bar.index.min()) - PRE_START max_yr = int(rho_bar.index.max()) - POST_END candidates = [y for y in range(min_yr, max_yr + 1) if y not in rec_years and y in rho_bar.index] if len(candidates) < 5: return (obs_depth, np.nan) np.random.seed(seed) n_as_deep = 0 for _ in range(n_perm): fake_tp = np.random.choice(candidates) r = _vshape_event_study(rho_bar, fake_tp) if r is not None and r['vshape'] and r['depth'] >= obs_depth: n_as_deep += 1 p_value = (n_as_deep + 1) / (n_perm + 1) return (obs_depth, p_value) # ── Structural break test ──────────────────────────────────────────────── def _structural_break_test(rho_bar, test_year=2001): """ Test whether test_year is a structural break in rho_bar(t). Compares: H0: rho_bar = a + b*t (single linear trend, no break) H1: rho_bar = a + b1*t + c*(t>bp) + d*max(t-bp,0) (two-segment) Returns dict with F-stat at test_year, optimal breakpoint, etc. """ years = np.array(sorted(rho_bar.index), dtype=float) vals = np.array([rho_bar[int(y)] for y in years]) n = len(years) # H0: single linear trend X0 = np.column_stack([np.ones(n), years]) b0 = np.linalg.lstsq(X0, vals, rcond=None)[0] rss0 = np.sum((vals - X0 @ b0) ** 2) def _two_segment_rss(bp): """RSS for two-segment model with breakpoint at bp.""" above = (years > bp).astype(float) excess = np.maximum(years - bp, 0) X1 = np.column_stack([np.ones(n), years, above, excess]) b1 = np.linalg.lstsq(X1, vals, rcond=None)[0] rss1 = np.sum((vals - X1 @ b1) ** 2) slope_pre = b1[1] slope_post = b1[1] + b1[3] return rss1, slope_pre, slope_post, b1 # F-test at the specified test_year rss_test, slope_pre_test, slope_post_test, _ = _two_segment_rss(test_year) f_stat_test = ((rss0 - rss_test) / 2) / (rss_test / (n - 4)) # p-value from F(2, n-4) if HAS_SCIPY: p_test = 1.0 - scipy_stats.f.cdf(f_stat_test, 2, n - 4) else: p_test = np.nan # Grid search for optimal breakpoint (by AIC) best_aic = n * np.log(rss0 / n) + 2 * 2 # AIC of linear model best_bp = None aic_linear = best_aic results_grid = {} for bp in range(int(years.min()) + 8, int(years.max()) - 5): n_pre = np.sum(years <= bp) n_post = np.sum(years > bp) if n_pre < 5 or n_post < 5: continue rss_bp, sp, spo, _ = _two_segment_rss(bp) aic_bp = n * np.log(rss_bp / n) + 2 * 4 f_bp = ((rss0 - rss_bp) / 2) / (rss_bp / (n - 4)) results_grid[bp] = {'rss': rss_bp, 'aic': aic_bp, 'f_stat': f_bp, 'slope_pre': sp, 'slope_post': spo} if aic_bp < best_aic: best_aic = aic_bp best_bp = bp # Compute F-stat and p for optimal breakpoint f_stat_opt = np.nan p_opt = np.nan if best_bp is not None and best_bp in results_grid: f_stat_opt = results_grid[best_bp]['f_stat'] if HAS_SCIPY: p_opt = 1.0 - scipy_stats.f.cdf(f_stat_opt, 2, n - 4) return { 'test_year': test_year, 'f_stat_test': f_stat_test, 'p_test': p_test, 'slope_pre_test': slope_pre_test, 'slope_post_test': slope_post_test, 'optimal_bp': best_bp, 'f_stat_optimal': f_stat_opt, 'p_optimal': p_opt, 'aic_linear': aic_linear, 'aic_optimal': best_aic, 'aic_improvement': aic_linear - best_aic, 'distance_2001_to_optimal': abs(test_year - best_bp) if best_bp else np.nan, } # ── Cycle period estimation ────────────────────────────────────────────── def _estimate_cycle_period(rho_bar): """ Estimate dominant cycle period from detrended rho_bar(t) autocorrelation. Returns dict with half_period, full_period, acf values. """ years = np.array(sorted(rho_bar.index), dtype=float) vals = np.array([rho_bar[int(y)] for y in years]) n = len(years) # Linear detrend t = years - years.mean() b = np.sum(t * (vals - vals.mean())) / np.sum(t ** 2) a = vals.mean() - b * years.mean() detrended = vals - (a + b * years) var_dt = np.var(detrended) if var_dt < 1e-12: return None # Compute ACF at lags 1 to n//2 max_lag = min(n // 2, 25) acf = {} for lag in range(1, max_lag + 1): acf[lag] = np.mean(detrended[lag:] * detrended[:-lag]) / var_dt # Find half-period (lag of most negative ACF) min_lag = min(acf, key=acf.get) min_acf = acf[min_lag] # Find first zero crossing zero_cross = None for lag in range(1, max_lag): if lag + 1 in acf and acf[lag] > 0 and acf[lag + 1] < 0: # Linear interpolation zero_cross = lag + acf[lag] / (acf[lag] - acf[lag + 1]) break full_period = 2 * min_lag full_period_zc = 2 * zero_cross if zero_cross else np.nan return { 'half_period': min_lag, 'full_period': full_period, 'min_acf': min_acf, 'zero_crossing': zero_cross, 'full_period_zc': full_period_zc, 'acf': acf, 'detrended': detrended, 'years': years, 'trend_slope': b, } # ── BEA capital intensity ─────────────────────────────────────────────── def _event_study_annual_kexcess(metric, wave): """Event study for annual integer-indexed K_excess series.""" tp = wave['tp'] pre_start = wave['install_start'] pre_end = tp post_start = tp + 1 dep_end = wave['deploy_end'] def period_mean(s, start, end): mask = (s.index >= start) & (s.index <= end) vals = s[mask].dropna() return vals.mean() if len(vals) > 0 else np.nan install_mean = period_mean(metric, pre_start, pre_end) deploy_mean = period_mean(metric, post_start, dep_end) correct_dir = (not np.isnan(install_mean) and not np.isnan(deploy_mean) and install_mean > deploy_mean) return { 'wave': wave['name'], 'canonical': wave.get('canonical', True), 'install_mean': install_mean, 'deploy_mean': deploy_mean, 'correct_dir': correct_dir, } # ── Main ────────────────────────────────────────────────────────────────── _cache = {} def main(): """Run rho V-shape tests. Returns list of result dicts.""" print("rho V-Shape at Technology Wave Turning Points") print("=" * 60) results_list = [] # ── Estimate rho_bar(t) ─────────────────────────────────────────── rho_bar = pd.Series(dtype=float) if HAS_NBER: try: nber = load_nber_ces() rho_bar = _estimate_rho_bar_annual(nber) print(f"\n NBER-CES rho_bar: {len(rho_bar)} years " f"({int(rho_bar.index.min())}-{int(rho_bar.index.max())}), " f"range [{rho_bar.min():.4f}, {rho_bar.max():.4f}]") except FileNotFoundError: print(" NBER-CES data not found") except Exception as e: print(f" Error: {e}") if len(rho_bar) >= 10: _cache['rho_bar'] = rho_bar # ══════════════════════════════════════════════════════════════════ # SUBTEST A: V-shape at ICT/Internet 2001 (theory prediction) # ══════════════════════════════════════════════════════════════════ print("\n[A] V-shape at 2001 (theory-predicted Perez turning point)...") if len(rho_bar) >= 10: tp_year = 2001 es = _vshape_event_study(rho_bar, tp_year) if es is not None: obs_depth, perm_p = _permutation_test_vshape(rho_bar, tp_year) _cache['es_2001'] = es _cache['perm_p_2001'] = perm_p print(f" PRE [{tp_year+PRE_START},{tp_year+PRE_END}): " f"mean={es['pre_mean']:.4f} (n={es['n_pre']})") print(f" ZONE [{tp_year-ZONE_HALF},{tp_year+ZONE_HALF}]: " f"mean={es['zone_mean']:.4f} (n={es['n_zone']})") print(f" POST ({tp_year+POST_START},{tp_year+POST_END}]: " f"mean={es['post_mean']:.4f} (n={es['n_post']})") print(f" V-shape: {'YES' if es['vshape'] else 'NO'}, " f"depth={es['depth']:.4f}, p={perm_p:.4f}") if es['vshape'] and perm_p < 0.10: verd = 'CONSISTENT' elif es['vshape']: verd = 'AMBIGUOUS' else: verd = 'INCONSISTENT' results_list.append(log_result( name='rho_vshape_nberces', stat=float(es['depth']), p=float(perm_p), verd=verd, details={ 'tp_year': tp_year, 'pre_mean': round(es['pre_mean'], 4), 'zone_mean': round(es['zone_mean'], 4), 'post_mean': round(es['post_mean'], 4), 'vshape': es['vshape'], 'depth': round(es['depth'], 4), 'slope_pre': round(es['slope_pre'], 5), 'slope_post': round(es['slope_post'], 5), 'n_perm': 5000, }, description=( "Tests whether Oberfield-Raval time-varying rho_bar(t) follows " "a V-shape at the ICT/Internet 2001 Perez turning point. " "Event study: mean rho_bar in PRE (tp-8 to tp-2), ZONE (tp-2 to tp+2), " "POST (tp+2 to tp+8). V-shape = ZONE < PRE and ZONE < POST. " "5000-iteration permutation test. " "Data: NBER-CES (1958-2018), ~460 industries, rolling 15-year OLS. " "AMBIGUOUS: V-shape technically present but depth is negligible " "(post-mean nearly equals zone-mean). The pattern is an L-shape: " "clear installation decline, no deployment recovery. The rho_bar " "minimum is at 2012, not 2001." ), )) print(f" -> {verd}") else: results_list.append(log_result( name='rho_vshape_nberces', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'insufficient data'}, description="Tests V-shape rho_bar(t) at 2001. Insufficient data." )) else: results_list.append(log_result( name='rho_vshape_nberces', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'NBER-CES unavailable'}, description="Tests V-shape rho_bar(t) at 2001. NBER-CES data unavailable." )) # ══════════════════════════════════════════════════════════════════ # SUBTEST B: V-shape at 1983 (method validation) # ══════════════════════════════════════════════════════════════════ print("\n[B] V-shape at 1983 (actual rho_bar minimum, method validation)...") if len(rho_bar) >= 10: tp_micro = 1983 es_micro = _vshape_event_study(rho_bar, tp_micro) if es_micro is not None: obs_depth_m, perm_p_m = _permutation_test_vshape( rho_bar, tp_micro, n_perm=5000, seed=43) _cache['es_micro'] = es_micro _cache['perm_p_micro'] = perm_p_m print(f" PRE [{tp_micro+PRE_START},{tp_micro+PRE_END}): " f"mean={es_micro['pre_mean']:.4f} (n={es_micro['n_pre']})") print(f" ZONE [{tp_micro-ZONE_HALF},{tp_micro+ZONE_HALF}]: " f"mean={es_micro['zone_mean']:.4f} (n={es_micro['n_zone']})") print(f" POST ({tp_micro+POST_START},{tp_micro+POST_END}]: " f"mean={es_micro['post_mean']:.4f} (n={es_micro['n_post']})") print(f" V-shape: {'YES' if es_micro['vshape'] else 'NO'}, " f"depth={es_micro['depth']:.4f}, p={perm_p_m:.4f}") if es_micro['vshape'] and perm_p_m < 0.10: verd_m = 'CONSISTENT' elif es_micro['vshape']: verd_m = 'AMBIGUOUS' else: verd_m = 'INCONSISTENT' results_list.append(log_result( name='rho_vshape_micro', stat=float(es_micro['depth']), p=float(perm_p_m), verd=verd_m, details={ 'tp_year': tp_micro, 'pre_mean': round(es_micro['pre_mean'], 4), 'zone_mean': round(es_micro['zone_mean'], 4), 'post_mean': round(es_micro['post_mean'], 4), 'vshape': es_micro['vshape'], 'depth': round(es_micro['depth'], 4), 'slope_pre': round(es_micro['slope_pre'], 5), 'slope_post': round(es_micro['slope_post'], 5), 'n_perm': 5000, }, description=( "Method validation: tests V-shape at 1983, the actual rho_bar " "minimum between the 1972 peak and 1992 peak. If CONSISTENT, the " "Oberfield-Raval estimation CAN detect V-shapes in rho_bar(t), and " "the absence of a V-shape at 2001 is a real finding (not a method " "limitation). Same event study protocol as rho_vshape_nberces. " "Note: 1983 is NOT a Perez turning point in the standard " "periodization (the paper places the entire 1971-2001 period as " "ICT installation). The V-shape at 1983 may reflect Kuznets-scale " "sub-cycles within the Kondratiev wave, or it may indicate that " "rho_bar dynamics are driven by macro regime shifts (end of Volcker " "tightening) rather than technology waves." ), )) print(f" -> {verd_m}") else: results_list.append(log_result( name='rho_vshape_micro', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'insufficient data'}, description="V-shape validation at 1983. Insufficient data." )) else: results_list.append(log_result( name='rho_vshape_micro', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'NBER-CES unavailable'}, description="V-shape validation at 1983. NBER-CES data unavailable." )) # ══════════════════════════════════════════════════════════════════ # SUBTEST C: Structural break at 2001 # ══════════════════════════════════════════════════════════════════ print("\n[C] Structural break test at 2001...") if len(rho_bar) >= 10: sb = _structural_break_test(rho_bar, test_year=2001) _cache['structural_break'] = sb print(f" At 2001: F={sb['f_stat_test']:.2f}, p={sb['p_test']:.4f}") print(f" slope_pre={sb['slope_pre_test']:+.5f}, " f"slope_post={sb['slope_post_test']:+.5f}") print(f" Optimal breakpoint: {sb['optimal_bp']} " f"(F={sb['f_stat_optimal']:.2f}, p={sb['p_optimal']:.4f})") print(f" AIC improvement (optimal vs linear): {sb['aic_improvement']:.1f}") print(f" Distance 2001 to optimal: {sb['distance_2001_to_optimal']:.0f} years") # Verdict: is 2001 the structural break? dist = sb['distance_2001_to_optimal'] if dist <= 3 and sb['p_test'] < 0.05: verd_sb = 'CONSISTENT' elif dist <= 5 and sb['p_test'] < 0.10: verd_sb = 'AMBIGUOUS' elif sb['p_test'] >= 0.10: verd_sb = 'INCONSISTENT' else: verd_sb = 'INCONSISTENT' results_list.append(log_result( name='rho_breakpoint_2001', stat=float(sb['f_stat_test']), p=float(sb['p_test']), verd=verd_sb, details={ 'f_stat_2001': round(sb['f_stat_test'], 2), 'p_2001': round(sb['p_test'], 4), 'slope_pre': round(sb['slope_pre_test'], 5), 'slope_post': round(sb['slope_post_test'], 5), 'optimal_bp': sb['optimal_bp'], 'f_stat_optimal': round(sb['f_stat_optimal'], 2), 'distance': sb['distance_2001_to_optimal'], 'aic_improvement': round(sb['aic_improvement'], 1), }, description=( "Tests whether the Perez turning point (2001) is a structural " "break in rho_bar(t). Segmented regression: compares single linear " "trend vs two-segment model with level+slope change at breakpoint. " "F-test for break significance. Grid search for optimal breakpoint " "(by AIC) over 1969-2008. " "INCONSISTENT: 2001 is NOT a significant structural break " f"(F={sb['f_stat_test']:.1f}, p={sb['p_test']:.3f}). " f"The data-optimal breakpoint is {sb['optimal_bp']}, " f"{sb['distance_2001_to_optimal']:.0f} years from the " "theory-predicted turning point. This suggests the Perez turning " "point is irrelevant for manufacturing rho_bar dynamics." ), )) print(f" -> {verd_sb}") else: results_list.append(log_result( name='rho_breakpoint_2001', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'NBER-CES unavailable'}, description="Structural break test at 2001. NBER-CES unavailable." )) # ══════════════════════════════════════════════════════════════════ # SUBTEST D: Cycle period estimation # ══════════════════════════════════════════════════════════════════ print("\n[D] Cycle period estimation...") if len(rho_bar) >= 10: cp = _estimate_cycle_period(rho_bar) _cache['cycle_period'] = cp if cp is not None: print(f" Half-period (most negative ACF): {cp['half_period']} years " f"(ACF={cp['min_acf']:.3f})") print(f" Full period estimate: {cp['full_period']} years") if cp['zero_crossing']: print(f" Zero-crossing estimate: {cp['full_period_zc']:.1f} years") print(f" Linear trend slope: {cp['trend_slope']:.5f}/yr") full_p = cp['full_period'] # Theory predicts 25-50 years (Kondratiev) if 25 <= full_p <= 50: verd_cp = 'CONSISTENT' elif 20 <= full_p < 25 or 50 < full_p <= 60: verd_cp = 'AMBIGUOUS' else: verd_cp = 'INCONSISTENT' results_list.append(log_result( name='rho_cycle_period', stat=float(full_p), p=np.nan, verd=verd_cp, details={ 'half_period': cp['half_period'], 'full_period': cp['full_period'], 'min_acf': round(cp['min_acf'], 3), 'zero_crossing': round(cp['zero_crossing'], 1) if cp['zero_crossing'] else None, 'full_period_zc': round(cp['full_period_zc'], 1) if cp['zero_crossing'] else None, 'trend_slope': round(cp['trend_slope'], 5), 'theory_range': '25-50 years (Kondratiev)', }, description=( "Estimates the dominant cycle period in detrended rho_bar(t) " "using autocorrelation analysis. The theory (Paper 3, Corollary " "on cycle period) predicts T_cycle ~ 25-50 years (Kondratiev " "timescale), with formula T = 2*tau_A + tau_fast + tau_S " "(~22-32 years from parameter estimates). " f"Finding: full period ~ {full_p} years (half-period = " f"{cp['half_period']} years, ACF = {cp['min_acf']:.3f}). " "This is at the low end of the theory's prediction, consistent " "with the formula but shorter than historical Perez wave periods " "(average ~52 years). May indicate sub-Kondratiev dynamics " "(Kuznets-scale ~20-year cycles in manufacturing rho) or that " "the rho_bar estimation captures macro regime cycles rather " "than technology wave cycles." ), )) print(f" -> {verd_cp} (theory predicts 25-50 years)") else: results_list.append(log_result( name='rho_cycle_period', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'insufficient variation'}, description="Cycle period estimation. Insufficient variation." )) else: results_list.append(log_result( name='rho_cycle_period', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'NBER-CES unavailable'}, description="Cycle period estimation. NBER-CES unavailable." )) # ══════════════════════════════════════════════════════════════════ # SUBTEST E: BEA capital intensity (ICT K_excess) # ══════════════════════════════════════════════════════════════════ print("\n[E] BEA capital intensity V-shape...") cap_path = BEA_FA_DIR / 'CapitalStockByIndustry.csv' if not cap_path.exists(): cap_path = DATA_DIR / 'compiled' / 'bea_capital_stock_by_industry.csv' if not cap_path.exists(): cap_path = BEA_FA_DIR / 'bea_capital_stock_by_industry.csv' if cap_path.exists(): cap = pd.read_csv(cap_path, index_col='year') if 'Information' in cap.columns and 'Total_Private' in cap.columns: g_info = cap['Information'].pct_change() * 100 g_total = cap['Total_Private'].pct_change() * 100 k_excess = g_info - g_total wave_results_b = [] for wave in PEREZ_WAVES_ANNUAL: r = _event_study_annual_kexcess(k_excess, wave) wave_results_b.append(r) dirstr = 'YES' if r['correct_dir'] else 'NO' print(f" {r['wave']}: direction={dirstr}, " f"install={r['install_mean']:.2f}, deploy={r['deploy_mean']:.2f}") n_correct = sum(r['correct_dir'] for r in wave_results_b) n_correct_canonical = sum(1 for r in wave_results_b if r['canonical'] and r['correct_dir']) verd_b = 'AMBIGUOUS' if n_correct_canonical >= 1 else 'INCONSISTENT' results_list.append(log_result( name='rho_vshape_inputs', stat=float(n_correct), p=np.nan, verd=verd_b, details={'n_correct': n_correct, 'n_correct_canonical': n_correct_canonical}, description=( "Tests BEA ICT capital excess growth direction at Perez " "turning points. Binary direction check: install > deploy. " "Data: BEA Fixed Assets (1980-2023). " "AMBIGUOUS: N=1 canonical wave, no statistical test." ), )) _cache['k_excess'] = k_excess _cache['wave_results_b'] = wave_results_b print(f" -> {verd_b}") else: results_list.append(log_result( name='rho_vshape_inputs', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'missing columns'}, description="BEA capital intensity. Missing columns." )) else: results_list.append(log_result( name='rho_vshape_inputs', stat=np.nan, p=np.nan, verd='AMBIGUOUS', details={'reason': 'file not found'}, description="BEA capital intensity. File not found." )) print(f"\n Results: {len(results_list)} tests") _cache['results'] = results_list return results_list def figures(): """Generate figures for each subtest.""" if not _cache: main() paths = [] import matplotlib.pyplot as plt rho_bar = _cache.get('rho_bar') if rho_bar is None or len(rho_bar) < 10: return paths # --- Figure 1: Full rho_bar time series with both event studies --- fig, axes = setup_subplots(1, 2, width=14, height=5) for idx, (tp, es_key, perm_key, title_suffix) in enumerate([ (2001, 'es_2001', 'perm_p_2001', 'ICT/Internet 2001 (theory)'), (1983, 'es_micro', 'perm_p_micro', '1983 (method validation)'), ]): ax = axes[idx] es = _cache.get(es_key, {}) perm_p = _cache.get(perm_key, np.nan) ax.plot(rho_bar.index, rho_bar.values, 'k-o', markersize=2.5, linewidth=0.8, label=r'$\hat\rho_{\mathrm{bar}}(t)$') # Shade event windows ax.axvspan(tp + PRE_START, tp + PRE_END, alpha=0.10, color='blue') ax.axvspan(tp - ZONE_HALF, tp + ZONE_HALF, alpha=0.15, color='red') ax.axvspan(tp + POST_START, tp + POST_END, alpha=0.10, color='green') ax.axvline(tp, color='red', linestyle='--', alpha=0.8) # NBER recessions for start, end in NBER_RECESSIONS: ax.axvspan(start, end + 1, alpha=0.08, color='grey') if es: depth = es.get('depth', 0) vshape = es.get('vshape', False) ann = (f"V: {'Yes' if vshape else 'No'}, " f"d={depth:.4f}\np={perm_p:.3f}") ax.text(0.02, 0.98, ann, transform=ax.transAxes, fontsize=8, va='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) ax.set_xlabel('Year') ax.set_ylabel(r'$\hat\rho_{\mathrm{bar}}$') ax.set_title(f'V-Shape at {title_suffix}') ax.grid(True, alpha=0.3) fig.tight_layout() paths.extend(save_figure(fig, 'rho_vshape_nberces', subfolder='dynamics')) # --- Figure 2: Structural break + cycle period --- sb = _cache.get('structural_break') cp = _cache.get('cycle_period') if sb and cp: fig2, axes2 = setup_subplots(1, 2, width=14, height=5) # Panel 1: Breakpoint F-statistics ax2a = axes2[0] bps = sorted([bp for bp in range(int(rho_bar.index.min()) + 8, int(rho_bar.index.max()) - 5)]) f_vals = [] for bp in bps: sb_bp = _structural_break_test(rho_bar, test_year=bp) f_vals.append(sb_bp['f_stat_test']) ax2a.plot(bps, f_vals, 'b-', linewidth=1.0) ax2a.axvline(2001, color='red', linestyle='--', alpha=0.7, label='2001 (theory)') if sb['optimal_bp']: ax2a.axvline(sb['optimal_bp'], color='green', linestyle='--', alpha=0.7, label=f"{sb['optimal_bp']} (optimal)") # F critical value at 5% if HAS_SCIPY: n = len(rho_bar) f_crit = scipy_stats.f.ppf(0.95, 2, n - 4) ax2a.axhline(f_crit, color='grey', linestyle=':', alpha=0.5, label=f'F crit (5%) = {f_crit:.1f}') ax2a.set_xlabel('Breakpoint Year') ax2a.set_ylabel('F-statistic') ax2a.set_title('Structural Break F-statistic by Year') ax2a.legend(fontsize=8) ax2a.grid(True, alpha=0.3) # Panel 2: Autocorrelation function ax2b = axes2[1] lags = sorted(cp['acf'].keys()) acf_vals = [cp['acf'][lag] for lag in lags] ax2b.bar(lags, acf_vals, color=['steelblue' if v >= 0 else 'salmon' for v in acf_vals], alpha=0.7) ax2b.axhline(0, color='black', linewidth=0.5) # 95% confidence bounds (approximate: ±2/sqrt(n)) n = len(rho_bar) ci = 2.0 / np.sqrt(n) ax2b.axhline(ci, color='grey', linestyle=':', alpha=0.5) ax2b.axhline(-ci, color='grey', linestyle=':', alpha=0.5) ax2b.axvline(cp['half_period'], color='red', linestyle='--', alpha=0.7, label=f"Half-period = {cp['half_period']}yr") ax2b.set_xlabel('Lag (years)') ax2b.set_ylabel('ACF') ax2b.set_title(f"Detrended ACF (full cycle ~ {cp['full_period']}yr)") ax2b.legend(fontsize=8) ax2b.grid(True, alpha=0.3) fig2.tight_layout() paths.extend(save_figure(fig2, 'rho_vshape_diagnostics', subfolder='dynamics')) # --- Figure 3: rho_vshape_inputs (BEA capital intensity) --- k_excess = _cache.get('k_excess') if k_excess is not None and len(k_excess.dropna()) > 5: fig3, ax3 = setup_figure(width=10, height=5) k_exc = k_excess.dropna() ax3.bar(k_exc.index, k_exc.values, color=['steelblue' if v >= 0 else 'salmon' for v in k_exc.values], alpha=0.7, edgecolor='none') ax3.axhline(0, color='black', linewidth=0.5) for wave in PEREZ_WAVES_ANNUAL: tp = wave['tp'] style = '--' if wave.get('canonical', True) else ':' ax3.axvline(tp, color='red', linestyle=style, alpha=0.7, label=f"{wave['name']}") ax3.axvspan(wave['install_start'], tp, alpha=0.08, color='orange') ax3.axvspan(tp + 1, wave['deploy_end'], alpha=0.08, color='green') ax3.set_xlabel('Year') ax3.set_ylabel(r'$\Delta g_{\mathrm{INFO}}$ (pp/yr)') ax3.set_title('ICT Capital Excess Growth at Perez Turning Points') ax3.legend(fontsize=8) ax3.grid(True, alpha=0.3) paths.extend(save_figure(fig3, 'rho_vshape_inputs', subfolder='dynamics')) return paths if __name__ == '__main__': results = main() for r in results: print(f" {r['name']:50s} {r['verdict']:15s} stat={r['statistic']}")

n waves	3
p value	1
n vshape	0
statistic	0
n neg slope	0
n pos slope	1

Rho Vshape Technology

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources