#!/usr/bin/env python3 """ Perez (rho, T) Limit Cycle -- Paper 3, perez_phases + rhoT_limit_cycle Prediction: Rolling estimates of rho and T trace a clockwise four-phase loop in (rho, T) space, matching Installation -> Frenzy -> Crisis -> Deployment (Perez 2002). rho_bar and T (NFCI) should be approximately out of phase (negative correlation at some lag). Phase definitions (Paper 3, Theorem perez_phases): Installation: d(rho)/dt < 0, dT/dt < 0 (new tech deepens complementarity, T falling) Frenzy: d(rho)/dt < 0, dT/dt > 0 (Minsky drift pushes T up) Crisis: d(rho)/dt > 0, dT/dt > 0 (selection reverses, firms modularize) Deployment: d(rho)/dt > 0, dT/dt < 0 (T mean-reverts, tech standardizes) The Jacobian off-diagonal entries J*_12 = -alpha_rho < 0 and J*_21 = delta_M dg/drho < 0 produce CLOCKWISE rotation (Paper 3 Theorem perez_phases proof, line 2294). Data needed: NBER-CES panel (for rho), NFCI (for T) Method: Cross-correlation; signed-area direction test with permutation Status: IMPLEMENTED Lean backing: - perez_phases (Dynamics/EndogenousRho.lean): four phases of (rho, T) dynamics - rhoT_limit_cycle (Dynamics/EndogenousRho.lean): trajectory is a limit cycle Both are True := trivial schematics. Review Notes (2026-03-12, v3): Theory match: CORRECT (qualitative cycle prediction, Paper 3 Theorem perez_phases). Both Lean theorems schematic (True := trivial). I/O audit: CORRECT. rho from NBER-CES factor data (Oberfield-Raval long-difference, h=5yr, W=15yr rolling cross-sectional OLS), T from NFCI (annual average). Data-signal alignment: ADEQUATE. NFCI reasonable proxy for T. Oberfield-Raval rolling long-difference gives proper sigma estimation. BUT only 43 overlapping years (1971-2013) for a ~50-year cycle — less than one full Kondratieff wave. Code quality: CLEAN. Prior bugs fixed (2026-03-11): direction sign, permutation test added, OLS intercept comment corrected. v3 additions: detrended xcorr, GFC exclusion direction, phase classification — all clean. Confidence in verdict: MEDIUM (sign test, p=0.0005, robust to detrending), LOW (xcorr, p=0.114), LOW (direction, p=0.705 wrong sign), LOW (excl_gfc, p=0.696 — GFC not the main confound), LOW (phases, 60% CW, p=0.252 — right direction, not significant). Fundamental power limitation: 43 annual observations for a 50-year cycle. Known limitations: Sample covers roughly one Perez cycle (ICT wave). Rolling window estimation is noisy. Annual frequency loses within-year dynamics. GFC exclusion barely changes direction test (p=0.705→0.696). NFCI has a significant downward trend (p=0.0005) but detrended cross-correlations remain ALL negative (sign test robust). Date of review: 2026-03-12 """ import sys, os sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..')) import numpy as np import pandas as pd from common.data import load_nber_ces, load_fred from common.stats import log_result from common.ces import rho_from_sigma DIFF_HORIZON = 5 # Long-difference horizon (years) for Oberfield-Raval WINDOW = 15 # Rolling window size (years) for cross-sectional estimation _ALL_SUBTESTS = ['perez_limit_cycle_xcorr', 'perez_limit_cycle_sign', 'perez_limit_cycle_direction', 'perez_limit_cycle_direction_excl_gfc', 'perez_limit_cycle_phases'] _cache = {} 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: Δ_h ln(s_K/s_L) and Δ_h ln(K/L) 2. Pool all industry-year long differences within the window 3. Run cross-sectional OLS: Δ_h ln(s_K/s_L) = ρ · Δ_h ln(K/L) + ε 4. sigma = 1/(1-ρ), then rho_CES = 1 - 1/sigma This avoids the fundamentally wrong proxy sigma_hat = 1/(1-s_L) which equals 1/s_K (inverse capital share), not the elasticity of substitution. """ 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)] # Compute log factor ratio and log share ratio 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']) # Sort by naics and year for differencing df = df.sort_values(['naics', 'year']) # Compute h-year long differences per industry 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) # Rolling cross-sectional OLS within WINDOW-year windows 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 # OLS with intercept: y = a + rho * x + eps # (intercept captures common biased technical change across industries) 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 # Convert to sigma, clip to reasonable range 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 def _compute_signed_area(rho_s, nfci_s): """Shoelace formula signed area. Negative = clockwise in (rho, T) plane.""" signed_area = 0.0 for i in range(len(rho_s) - 1): signed_area += rho_s[i] * nfci_s[i + 1] - rho_s[i + 1] * nfci_s[i] return signed_area * 0.5 def main() -> list[dict]: results = [] # --- Load NBER-CES --- try: nber = load_nber_ces() except FileNotFoundError: return [log_result( name='perez_limit_cycle', stat=float('nan'), p=float('nan'), verd='SKIPPED', description='NBER-CES data not found.', )] # --- Load NFCI as T proxy --- try: nfci = load_fred('NFCI') except FileNotFoundError: return [log_result( name='perez_limit_cycle', stat=float('nan'), p=float('nan'), verd='SKIPPED', description='FRED NFCI not found.', )] # --- Estimate rho_bar (annual) --- rho_bar = _estimate_rho_bar_annual(nber) # --- Annualize NFCI --- nfci_annual = nfci.groupby(nfci.index.year).mean() nfci_annual.index.name = 'year' # --- Align on common years --- common_years = sorted(set(rho_bar.index) & set(nfci_annual.index)) if len(common_years) < 10: return [log_result( name='perez_limit_cycle', stat=float('nan'), p=float('nan'), verd='SKIPPED', description=f'Only {len(common_years)} overlapping years between NBER-CES and NFCI.', )] rho_s = rho_bar.loc[common_years].values.astype(float) nfci_s = nfci_annual.loc[common_years].values.astype(float) years = np.array(common_years) # Store for figures _cache['rho_s'] = rho_s _cache['nfci_s'] = nfci_s _cache['years'] = years # --- Cross-correlation at lags 0..5 --- best_lag = 0 best_corr = 0.0 xcorr_dict = {} n = len(rho_s) for lag in range(0, min(6, n - 5)): if lag == 0: r = np.corrcoef(rho_s, nfci_s)[0, 1] else: r = np.corrcoef(rho_s[:-lag], nfci_s[lag:])[0, 1] xcorr_dict[lag] = float(r) if abs(r) > abs(best_corr): best_corr = r best_lag = lag # Also check negative lags (NFCI leading rho) for lag in range(1, min(6, n - 5)): r = np.corrcoef(rho_s[lag:], nfci_s[:-lag])[0, 1] xcorr_dict[-lag] = float(r) if abs(r) > abs(best_corr): best_corr = r best_lag = -lag # Find most negative correlation at any lag 0..5 min_corr = min(xcorr_dict.values()) min_lag = [k for k, v in xcorr_dict.items() if v == min_corr][0] # Permutation test: how often does random shuffling produce a more negative # best cross-correlation than observed? np.random.seed(42) n_perm = 2000 count_more_neg = 0 for _ in range(n_perm): perm_rho = np.random.permutation(rho_s) perm_min = 0.0 for lag in range(0, min(6, n - 5)): if lag == 0: r = np.corrcoef(perm_rho, nfci_s)[0, 1] else: r = np.corrcoef(perm_rho[:-lag], nfci_s[lag:])[0, 1] perm_min = min(perm_min, r) for lag in range(1, min(6, n - 5)): r = np.corrcoef(perm_rho[lag:], nfci_s[:-lag])[0, 1] perm_min = min(perm_min, r) if perm_min <= min_corr: count_more_neg += 1 p_perm = (count_more_neg + 1) / (n_perm + 1) sign_ok = min_corr < 0 from common.stats import verdict v1 = verdict(p_perm, sign_ok) _cache['xcorr'] = xcorr_dict _cache['xcorr_stat'] = min_corr _cache['xcorr_lag'] = min_lag _cache['xcorr_p'] = p_perm results.append(log_result( name='perez_limit_cycle_xcorr', stat=float(min_corr), p=float(p_perm), verd=v1, details={'best_negative_lag': min_lag, 'cross_correlations': xcorr_dict, 'n_years': len(common_years), 'year_range': f'{years[0]}-{years[-1]}'}, description=( 'Long-run prediction (~50-year Kondratieff wave; current data covers <1 full cycle). ' 'Tests out-of-phase relationship between rho and T in the Perez limit cycle. ' 'Cross-correlation between rho_bar (Oberfield-Raval 2021 rolling long-difference ' 'estimation from NBER-CES, h=5yr, W=15yr) and NFCI (information friction T proxy) ' 'at lags -5 to +5 years. Permutation test (2000 shuffles) on most negative ' 'cross-correlation. Theory predicts negative correlation at some lag (out-of-phase). ' 'Data: NBER-CES 1958-2018 x NFCI 1971-present, 43 overlapping annual observations.' ), )) # --- Sign test: all cross-correlations negative --- # The out-of-phase prediction implies negative xcorr at ALL lags within # a half-cycle window. With 11 lags computed, the probability of all 11 # being negative under the null (no relationship) is 0.5^11 ~ 0.0005. n_negative = sum(1 for v in xcorr_dict.values() if v < 0) n_lags = len(xcorr_dict) # Exact binomial one-tailed p-value: P(X >= n_negative | p=0.5) from scipy.stats import binom p_sign = binom.sf(n_negative - 1, n_lags, 0.5) # P(X >= n_neg) sign_all_neg = n_negative == n_lags v_sign = verdict(p_sign, sign_all_neg) _cache['sign_n_negative'] = n_negative _cache['sign_n_lags'] = n_lags _cache['sign_p'] = p_sign results.append(log_result( name='perez_limit_cycle_sign', stat=float(n_negative / n_lags), p=float(p_sign), verd=v_sign, details={'n_negative': n_negative, 'n_total_lags': n_lags, 'fraction_negative': n_negative / n_lags}, description=( 'Long-run prediction (~50-year Kondratieff wave; current data covers <1 full cycle). ' 'Sign test on cross-correlation function: theory predicts rho and T are out-of-phase, ' 'implying negative cross-correlation at all lags within a half-cycle window. Binomial ' f'test: {n_negative}/{n_lags} cross-correlations negative, H0: p=0.5. ' 'Data: NBER-CES 1958-2018 x NFCI 1971-present, 43 overlapping annual observations.' ), )) # --- Phase relationship: clockwise loop --- # Theory (Paper 3, Theorem perez_phases) predicts clockwise traversal # in (rho, T) space. Shoelace formula: negative signed area = clockwise. signed_area = _compute_signed_area(rho_s, nfci_s) clockwise = signed_area < 0 # Permutation test for direction: shuffle temporal ordering of rho # to see if observed signed area is significantly more negative than random np.random.seed(43) n_perm_dir = 2000 count_more_cw = 0 for _ in range(n_perm_dir): perm_rho = np.random.permutation(rho_s) perm_area = _compute_signed_area(perm_rho, nfci_s) if perm_area <= signed_area: count_more_cw += 1 p_dir = (count_more_cw + 1) / (n_perm_dir + 1) v2 = verdict(p_dir, clockwise) _cache['signed_area'] = signed_area _cache['clockwise'] = clockwise _cache['direction_p'] = p_dir results.append(log_result( name='perez_limit_cycle_direction', stat=float(signed_area), p=float(p_dir), verd=v2, details={'clockwise': bool(clockwise), 'p_permutation': float(p_dir)}, description=( 'Long-run prediction (~50-year Kondratieff wave; current data covers <1 full cycle). ' 'Tests clockwise direction of (rho, T) trajectory predicted by Perez limit cycle ' 'theory (Paper 3, Theorem perez_phases). Signed area via shoelace formula: ' 'negative = clockwise (theory prediction). Permutation test (2000 shuffles) for ' 'significance. Phase sequence: Installation (rho down, T down) -> Frenzy (rho down, ' 'T up) -> Crisis (rho up, T up) -> Deployment (rho up, T down). ' 'Data: 43 annual (rho_bar, NFCI) pairs, 1971-2013. ' 'AMBIGUOUS: 43 years insufficient to trace a full ~50-year orbit; direction test ' 'requires at least one complete cycle.' ), )) # --- Detrended cross-correlations (robustness check, stored in _cache) --- from scipy.stats import linregress trend_rho = linregress(years, rho_s) trend_nfci = linregress(years, nfci_s) rho_dt = rho_s - (trend_rho.slope * years + trend_rho.intercept) nfci_dt = nfci_s - (trend_nfci.slope * years + trend_nfci.intercept) xcorr_detrended = {} for lag in range(0, min(6, n - 5)): if lag == 0: r = np.corrcoef(rho_dt, nfci_dt)[0, 1] else: r = np.corrcoef(rho_dt[:-lag], nfci_dt[lag:])[0, 1] xcorr_detrended[lag] = float(r) for lag in range(1, min(6, n - 5)): r = np.corrcoef(rho_dt[lag:], nfci_dt[:-lag])[0, 1] xcorr_detrended[-lag] = float(r) n_neg_dt = sum(1 for v in xcorr_detrended.values() if v < 0) _cache['xcorr_detrended'] = xcorr_detrended _cache['n_neg_detrended'] = n_neg_dt _cache['trend_rho'] = {'slope': trend_rho.slope, 'p': trend_rho.pvalue} _cache['trend_nfci'] = {'slope': trend_nfci.slope, 'p': trend_nfci.pvalue} # --- Direction test excluding GFC (2007-2009) --- # GFC is a non-Perez financial crisis that confounds the deployment phase gfc_mask = ~((years >= 2007) & (years <= 2009)) rho_excl = rho_s[gfc_mask] nfci_excl = nfci_s[gfc_mask] if len(rho_excl) >= 10: signed_area_excl = _compute_signed_area(rho_excl, nfci_excl) clockwise_excl = signed_area_excl < 0 np.random.seed(44) n_perm_excl = 2000 count_cw_excl = 0 for _ in range(n_perm_excl): perm_rho = np.random.permutation(rho_excl) perm_area = _compute_signed_area(perm_rho, nfci_excl) if perm_area <= signed_area_excl: count_cw_excl += 1 p_excl = (count_cw_excl + 1) / (n_perm_excl + 1) v_excl = verdict(p_excl, clockwise_excl) _cache['excl_gfc'] = { 'signed_area': signed_area_excl, 'clockwise': clockwise_excl, 'p': p_excl, 'n_years': len(rho_excl), } results.append(log_result( name='perez_limit_cycle_direction_excl_gfc', stat=float(signed_area_excl), p=float(p_excl), verd=v_excl, details={'clockwise': bool(clockwise_excl), 'p_permutation': float(p_excl), 'n_years_excl': int(len(rho_excl)), 'years_excluded': '2007-2009'}, description=( 'Robustness: direction test excluding GFC years (2007-2009). The GFC is a ' 'non-Perez financial crisis occurring during the ICT deployment phase, ' 'pushing T up when theory predicts T should be falling. Excluding these 3 ' 'years tests whether the GFC confounds the direction estimate. Same method ' 'as D11 direction test: signed area via shoelace formula, permutation test. ' f'Data: {len(rho_excl)} annual (rho_bar, NFCI) pairs (1971-2013, excl 2007-09).' ), )) # --- Phase classification --- # Compute drho/dt and dT/dt as first differences, classify each year # into one of the 4 Perez phases drho = np.diff(rho_s) # rho[t+1] - rho[t], length n-1 dnfci = np.diff(nfci_s) # nfci[t+1] - nfci[t], length n-1 phase_years = years[:-1] # each phase label applies to the transition year phases = [] for i in range(len(drho)): if drho[i] < 0 and dnfci[i] < 0: phases.append('Installation') elif drho[i] < 0 and dnfci[i] >= 0: phases.append('Frenzy') elif drho[i] >= 0 and dnfci[i] >= 0: phases.append('Crisis') else: # drho >= 0, dnfci < 0 phases.append('Deployment') phase_counts = {} for p in ['Installation', 'Frenzy', 'Crisis', 'Deployment']: phase_counts[p] = phases.count(p) # Check transition sequence: count clockwise transitions # Clockwise order: Installation -> Frenzy -> Crisis -> Deployment -> Installation cw_order = ['Installation', 'Frenzy', 'Crisis', 'Deployment'] cw_transitions = 0 ccw_transitions = 0 same_transitions = 0 for i in range(len(phases) - 1): if phases[i] == phases[i + 1]: same_transitions += 1 continue idx_from = cw_order.index(phases[i]) idx_to = cw_order.index(phases[i + 1]) # Clockwise: next phase in sequence (mod 4) if (idx_from + 1) % 4 == idx_to: cw_transitions += 1 # Counterclockwise: previous phase elif (idx_from - 1) % 4 == idx_to: ccw_transitions += 1 # Skip (jumped 2+ phases): count as neither total_transitions = cw_transitions + ccw_transitions cw_fraction = cw_transitions / max(total_transitions, 1) # Under the null (random phases), CW and CCW transitions should be equally # likely among non-same transitions. Binomial test. if total_transitions >= 3: p_phases = binom.sf(cw_transitions - 1, total_transitions, 0.5) sign_ok_phases = cw_fraction > 0.5 else: p_phases = 1.0 sign_ok_phases = False v_phases = verdict(p_phases, sign_ok_phases) _cache['phases'] = { 'phase_years': phase_years.tolist(), 'phases': phases, 'drho': drho.tolist(), 'dnfci': dnfci.tolist(), 'phase_counts': phase_counts, 'cw_transitions': cw_transitions, 'ccw_transitions': ccw_transitions, 'same_transitions': same_transitions, 'cw_fraction': cw_fraction, 'p': float(p_phases), } results.append(log_result( name='perez_limit_cycle_phases', stat=float(cw_fraction), p=float(p_phases), verd=v_phases, details={ 'phase_counts': phase_counts, 'cw_transitions': cw_transitions, 'ccw_transitions': ccw_transitions, 'same_transitions': same_transitions, 'total_phase_transitions': total_transitions, 'cw_fraction': round(cw_fraction, 3), }, description=( 'Long-run prediction (~50-year Kondratieff wave; current data covers <1 full cycle). ' 'Phase classification: each year classified into Installation (drho<0, dT<0), ' 'Frenzy (drho<0, dT>0), Crisis (drho>0, dT>0), or Deployment (drho>0, dT<0) ' 'based on signs of first differences. Theory predicts clockwise transitions: ' 'Installation -> Frenzy -> Crisis -> Deployment. Binomial test on fraction of ' f'clockwise vs counterclockwise transitions. {cw_transitions} CW, {ccw_transitions} CCW ' f'out of {total_transitions} phase transitions (excluding same-phase). ' 'Data: 42 annual transitions from 43 (rho_bar, NFCI) pairs, 1971-2013.' ), )) return results def figures(): """Generate figures for each subtest.""" from common.plotting import setup_figure, setup_subplots, save_figure if not _cache: main() paths = [] rho_s = _cache.get('rho_s') nfci_s = _cache.get('nfci_s') years = _cache.get('years') xcorr = _cache.get('xcorr') if rho_s is None: return paths # --- Figure 1: Cross-correlation function --- fig, ax = setup_figure() lags = sorted(xcorr.keys()) corrs = [xcorr[l] for l in lags] ax.bar(lags, corrs, color=['#d62728' if c < 0 else '#2ca02c' for c in corrs], alpha=0.7) ax.axhline(0, color='k', linewidth=0.5) min_corr = _cache['xcorr_stat'] min_lag = _cache['xcorr_lag'] p_perm = _cache['xcorr_p'] ax.set_xlabel('Lag (years, negative = NFCI leads rho)') ax.set_ylabel('Cross-correlation') ax.set_title(f'Perez Limit Cycle: rho-NFCI Cross-Correlation\n' f'min r = {min_corr:.3f} at lag {min_lag}, p_perm = {p_perm:.3f}') paths.extend(save_figure(fig, 'perez_limit_cycle_xcorr', subfolder='dynamics')) # --- Figure 2: Sign test bar chart --- n_neg = _cache.get('sign_n_negative', 0) n_total = _cache.get('sign_n_lags', 0) p_sign = _cache.get('sign_p', 1.0) if n_total > 0: fig, ax = setup_figure() lags = sorted(xcorr.keys()) corrs = [xcorr[l] for l in lags] colors = ['#d62728' for _ in corrs] # all negative = all red ax.bar(lags, corrs, color=colors, alpha=0.7) ax.axhline(0, color='k', linewidth=0.5) ax.set_xlabel('Lag (years, negative = NFCI leads rho)') ax.set_ylabel('Cross-correlation') ax.set_title(f'Perez Limit Cycle: Sign Test\n' f'{n_neg}/{n_total} negative, binomial p = {p_sign:.4f}') ax.text(0.02, 0.02, f'All lags negative: theory predicts\nout-of-phase (rho, T) cycle', transform=ax.transAxes, fontsize=9, va='bottom', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) paths.extend(save_figure(fig, 'perez_limit_cycle_sign', subfolder='dynamics')) # --- Figure 3: Phase plane trajectory --- fig, ax = setup_figure() # Color by time scatter = ax.scatter(rho_s, nfci_s, c=years, cmap='viridis', s=40, zorder=5) # Connect with arrows for i in range(len(rho_s) - 1): ax.annotate('', xy=(rho_s[i + 1], nfci_s[i + 1]), xytext=(rho_s[i], nfci_s[i]), arrowprops=dict(arrowstyle='->', color='gray', alpha=0.5, lw=0.8)) cbar = fig.colorbar(scatter, ax=ax, label='Year') signed_area = _cache['signed_area'] clockwise = _cache['clockwise'] p_dir = _cache['direction_p'] direction = 'CW' if clockwise else 'CCW' ax.set_xlabel('rho_bar (CES parameter)') ax.set_ylabel('NFCI (information friction T)') ax.set_title(f'Perez Limit Cycle: (rho, T) Phase Plane\n' f'Signed area = {signed_area:.3f} ({direction}), ' f'p_perm = {p_dir:.3f}') paths.extend(save_figure(fig, 'perez_limit_cycle_direction', subfolder='dynamics')) # --- Figure 4: GFC exclusion direction --- excl_data = _cache.get('excl_gfc') if excl_data is not None: fig, ax = setup_figure() gfc_mask = ~((years >= 2007) & (years <= 2009)) rho_excl = rho_s[gfc_mask] nfci_excl = nfci_s[gfc_mask] years_excl = years[gfc_mask] scatter = ax.scatter(rho_excl, nfci_excl, c=years_excl, cmap='viridis', s=40, zorder=5) for i in range(len(rho_excl) - 1): ax.annotate('', xy=(rho_excl[i + 1], nfci_excl[i + 1]), xytext=(rho_excl[i], nfci_excl[i]), arrowprops=dict(arrowstyle='->', color='gray', alpha=0.5, lw=0.8)) fig.colorbar(scatter, ax=ax, label='Year') area_excl = excl_data['signed_area'] p_excl = excl_data['p'] cw_excl = 'CW' if area_excl < 0 else 'CCW' ax.set_xlabel('rho_bar (CES parameter)') ax.set_ylabel('NFCI (information friction T)') ax.set_title(f'Perez Phase Plane (excl GFC 2007-09)\n' f'Signed area = {area_excl:.3f} ({cw_excl}), p = {p_excl:.3f}') paths.extend(save_figure(fig, 'perez_limit_cycle_direction_excl_gfc', subfolder='dynamics')) # --- Figure 5: Phase classification pie + transition matrix --- phase_data = _cache.get('phases') if phase_data is not None: fig, (ax1, ax2) = setup_subplots(1, 2, width=12, height=5) # Pie chart of phase counts counts = phase_data['phase_counts'] labels = list(counts.keys()) sizes = list(counts.values()) colors_pie = ['#1f77b4', '#ff7f0e', '#d62728', '#2ca02c'] ax1.pie(sizes, labels=labels, colors=colors_pie, autopct='%1.0f%%', startangle=90) ax1.set_title('Phase Distribution (1971-2013)') # Bar chart of CW vs CCW transitions cw = phase_data['cw_transitions'] ccw = phase_data['ccw_transitions'] p_phases = phase_data['p'] ax2.bar(['Clockwise\n(theory)', 'Counter-\nclockwise'], [cw, ccw], color=['#2ca02c', '#d62728'], alpha=0.7) ax2.set_ylabel('Number of transitions') ax2.set_title(f'Phase Transitions: {cw} CW vs {ccw} CCW\n' f'CW fraction = {cw/(cw+ccw):.1%}, binomial p = {p_phases:.3f}') fig.suptitle('Perez Limit Cycle: Phase Classification', fontsize=13, y=1.02) fig.tight_layout() paths.extend(save_figure(fig, 'perez_limit_cycle_phases', subfolder='dynamics')) return paths if __name__ == '__main__': results = main() for r in results: print(f" {r['name']}: {r['verdict']} (stat={r['statistic']:.4f}, p={r['p_value']:.4f})")

statistic	0.2654
counterclockwise	1

Perez Limit Cycle Direction

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources