#!/usr/bin/env python3 """ Tsallis Tail Distribution Test ============================== Paper 2 (CES Potential), §11.6 Manufacturing Tail Distributions Tests whether q-exponential (Tsallis) distributions fit manufacturing IP absolute log-returns better than exponential (Shannon) distributions. Per-sector test (tsallis_tails_qexp): Fits exponential and q-exponential MLEs to absolute monthly IP log-returns. Likelihood ratio test LR ~ chi2(1) per sector. CONSISTENT if majority of 17 sectors show significant q-exp superiority. GARCH robustness (tsallis_tails_garch_robustness): Re-runs LR test on GARCH(1,1) standardized residuals to confirm heavy tails are not a volatility clustering artifact. Data: 17 leaf FRED Manufacturing IP series (monthly, cached) Note: The cross-sector q̂ vs ρ̂ regression (formerly tsallis_tails_q_rho_regression) was dropped after evaluation (2026-03-11). The q = ρ identification applies to equilibrium allocation distributions, not time-series output returns. The proxy gap is unbridgeable with current data (q̂ ≈ 1.0 vs ρ ∈ 0.54–0.85). See evaluations/Tsallis Tails Q Rho Regression.md for full analysis. Review Notes (2026-03-11) ------------------------- Theory match: PARTIAL — q = ρ identification is a theorem (Lean: tsallis_uniqueness, schematic). Theory derives q-exponential for equilibrium allocations, not output returns. Within-sector LR test is valid (tests q ≠ 1, not q = ρ). I/O audit: CORRECT for within-sector test (tests heavy tails vs exponential). Data-signal: ADEQUATE — ~600 monthly obs per sector, absolute log-returns are stationary. GARCH control addresses volatility clustering confound. Code quality: CLEAN after fixes. Confidence: MEDIUM (qexp and garch robustness). """ import warnings warnings.filterwarnings('ignore') import numpy as np import pandas as pd from scipy import stats from scipy.optimize import minimize from common.config import COMPILED_DIR, RESULTS_DIR, FIGURES_DIR from common.data import load_fred from common.stats import log_result from common.plotting import setup_figure, setup_subplots, save_figure START_DATE = '1972-01-01' LEAF_SERIES = { 'IPG311A2S': 'Food/Bev/Tobacco', 'IPG313A4S': 'Textiles', 'IPG315A6S': 'Apparel/Leather', 'IPG321S': 'Wood Products', 'IPG322S': 'Paper', 'IPG323S': 'Printing', 'IPG324S': 'Petroleum/Coal', 'IPG325S': 'Chemicals', 'IPG326S': 'Plastics/Rubber', 'IPG331S': 'Primary Metals', 'IPG332S': 'Fabricated Metals', 'IPG333S': 'Machinery', 'IPG334S': 'Computer/Electronics', 'IPG335S': 'Electrical Equipment', 'IPG336S': 'Transport Equipment', 'IPG337S': 'Furniture', 'IPG339S': 'Misc Manufacturing', } def q_exp_loglik(params, data): """Negative log-likelihood for q-exponential distribution.""" q, beta = params if beta <= 0 or q >= 2.0 or q <= 0.0: return 1e12 n = len(data) if abs(q - 1.0) < 1e-10: return -(n * np.log(beta) - beta * np.sum(data)) if q < 2: log_C = np.log(beta) + np.log(2 - q) else: return 1e12 bracket = 1.0 - (1.0 - q) * beta * data if np.any(bracket <= 0): return 1e12 log_pdf = (1.0 / (1.0 - q)) * np.log(bracket) total = n * log_C + np.sum(log_pdf) return -total def fit_exponential(data): """MLE for exponential distribution.""" lam = 1.0 / np.mean(data) ll = len(data) * np.log(lam) - lam * np.sum(data) return lam, ll def fit_q_exponential(data): """MLE for q-exponential distribution.""" lam0 = 1.0 / np.mean(data) best_nll = 1e12 best_params = (1.0, lam0) for q0 in np.arange(0.3, 1.8, 0.1): try: res = minimize(q_exp_loglik, x0=[q0, lam0], args=(data,), method='Nelder-Mead', options={'maxiter': 2000, 'xatol': 1e-6, 'fatol': 1e-6}) if res.fun < best_nll and 0 < res.x[0] < 2 and res.x[1] > 0: best_nll = res.fun best_params = res.x except Exception: continue try: res2 = minimize(q_exp_loglik, x0=best_params, args=(data,), method='L-BFGS-B', bounds=[(0.01, 1.99), (0.01, None)], options={'maxiter': 1000}) if res2.fun < best_nll: best_nll = res2.fun best_params = res2.x except Exception: pass q_hat, beta_hat = best_params ll = -best_nll return q_hat, beta_hat, ll def load_all_returns(): """Load all leaf series and compute absolute log-returns.""" returns = {} for sid, label in LEAF_SERIES.items(): try: s = load_fred(sid) s = s[s.index >= START_DATE] except FileNotFoundError: continue lr = np.log(s / s.shift(1)).dropna() abs_lr = lr.abs() abs_lr = abs_lr[abs_lr > 0] returns[sid] = abs_lr return returns _cache = {} def main(): returns = load_all_returns() if len(returns) < 10: return [log_result('tsallis_tails_qexp', np.nan, np.nan, 'SKIPPED', details={'reason': f'Only {len(returns)} series loaded'}, description="Tests whether q-exponential (Tsallis) distributions fit manufacturing " "IP absolute log-returns better than exponential (Shannon) distributions, " "via likelihood ratio tests across 17 FRED leaf manufacturing IP series.", )] results_list = [] sector_results = [] for sid in sorted(returns.keys()): data = returns[sid].values lam_hat, ll_exp = fit_exponential(data) q_hat, beta_hat, ll_qexp = fit_q_exponential(data) lr_stat = max(2 * (ll_qexp - ll_exp), 0) lr_pval = 1 - stats.chi2.cdf(lr_stat, df=1) sector_results.append({ 'series': sid, 'q_hat': q_hat, 'lr_stat': lr_stat, 'lr_pval': lr_pval, }) df = pd.DataFrame(sector_results) n_sig = (df['lr_pval'] < 0.05).sum() mean_q = df['q_hat'].mean() # q-exponential significantly better in majority of sectors majority_sig = n_sig > len(df) / 2 if majority_sig: verd_qexp = 'CONSISTENT' elif n_sig >= 5: verd_qexp = 'AMBIGUOUS' else: verd_qexp = 'INCONSISTENT' results_list.append(log_result( 'tsallis_tails_qexp', float(n_sig) / len(df), np.nan, verd_qexp, details={ 'n_sectors': len(df), 'n_significant_05': int(n_sig), 'mean_q_hat': mean_q, }, description="Tests whether q-exponential (Tsallis) fits manufacturing IP absolute log-returns " "better than exponential (Shannon). LR test: 2*(ll_qexp - ll_exp) ~ chi2(1), per " "sector. Data: 17 FRED leaf IP series, monthly |log(IP_t/IP_{t-1})| from 1972. " "Prediction: majority of sectors show significant q-exp superiority (heavy tails).", )) # GARCH robustness garch_results = [] try: from arch import arch_model garch_n_sig = 0 garch_total = 0 for sid in sorted(returns.keys()): try: s = load_fred(sid) s = s[s.index >= START_DATE] lr = np.log(s / s.shift(1)).dropna() lr = lr[lr != 0] lr_scaled = lr * 100 am = arch_model(lr_scaled, vol='Garch', p=1, q=1, dist='normal') res = am.fit(disp='off', show_warning=False) std_resid = res.std_resid.dropna() abs_resid = np.abs(std_resid.values) abs_resid = abs_resid[abs_resid > 0] _, ll_exp_g = fit_exponential(abs_resid) _, _, ll_qexp_g = fit_q_exponential(abs_resid) lr_stat_g = max(2 * (ll_qexp_g - ll_exp_g), 0) lr_pval_g = 1 - stats.chi2.cdf(lr_stat_g, df=1) garch_total += 1 if lr_pval_g < 0.05: garch_n_sig += 1 garch_results.append({'series': sid, 'lr_stat': lr_stat_g, 'lr_pval': lr_pval_g}) except Exception: continue if garch_total >= 5: garch_frac = garch_n_sig / garch_total verd_garch = 'CONSISTENT' if garch_frac > 0.5 else 'AMBIGUOUS' results_list.append(log_result( 'tsallis_tails_garch_robustness', garch_frac, np.nan, verd_garch, details={'n_sig': garch_n_sig, 'n_total': garch_total}, description="GARCH robustness: re-runs the q-exponential vs exponential LR test on " "GARCH(1,1) standardized residuals instead of raw returns. If heavy tails " "persist after volatility clustering is removed, it confirms the CES-Tsallis " "prediction is not an artifact of conditional heteroskedasticity.", )) except ImportError: pass _cache['sector_df'] = df _cache['returns'] = returns _cache['results'] = results_list _cache['garch_results'] = garch_results return results_list def figures(): """Generate Tsallis tails figures: QQ plots + LR bar chart, GARCH comparison.""" if not _cache: main() paths = [] df = _cache.get('sector_df') returns = _cache.get('returns', {}) if df is None or len(df) == 0: return paths fig, axes = setup_subplots(1, 3, width=15, height=5) # (a) QQ plot for best sector (highest LR stat) ax = axes[0] best_row = df.loc[df['lr_stat'].idxmax()] sid = best_row['series'] if sid in returns: data = np.sort(returns[sid].values) n = len(data) theoretical_q = stats.expon.ppf(np.arange(1, n + 1) / (n + 1), scale=np.mean(data)) ax.scatter(theoretical_q, data, s=4, alpha=0.5, color='steelblue') lim = max(data.max(), theoretical_q.max()) * 1.05 ax.plot([0, lim], [0, lim], 'r--', lw=1, label='1:1 line') ax.legend(fontsize=8) ax.set_xlabel('Exponential quantiles') ax.set_ylabel('Empirical quantiles') ax.set_title(f'(a) Best sector: {LEAF_SERIES.get(sid, sid)}') ax.grid(True, alpha=0.3) # (b) QQ plot for worst sector (lowest LR stat) ax = axes[1] worst_row = df.loc[df['lr_stat'].idxmin()] sid2 = worst_row['series'] if sid2 in returns: data2 = np.sort(returns[sid2].values) n2 = len(data2) theoretical_q2 = stats.expon.ppf(np.arange(1, n2 + 1) / (n2 + 1), scale=np.mean(data2)) ax.scatter(theoretical_q2, data2, s=4, alpha=0.5, color='darkorange') lim2 = max(data2.max(), theoretical_q2.max()) * 1.05 ax.plot([0, lim2], [0, lim2], 'r--', lw=1, label='1:1 line') ax.legend(fontsize=8) ax.set_xlabel('Exponential quantiles') ax.set_ylabel('Empirical quantiles') ax.set_title(f'(b) Worst sector: {LEAF_SERIES.get(sid2, sid2)}') ax.grid(True, alpha=0.3) # (c) LR bar chart ax = axes[2] sorted_df = df.sort_values('lr_stat', ascending=True) colors = ['steelblue' if p < 0.05 else 'lightgrey' for p in sorted_df['lr_pval']] y_pos = range(len(sorted_df)) ax.barh(y_pos, sorted_df['lr_stat'], color=colors, edgecolor='grey', lw=0.5) ax.set_yticks(y_pos) ax.set_yticklabels([LEAF_SERIES.get(s, s) for s in sorted_df['series']], fontsize=7) crit = stats.chi2.ppf(0.95, df=1) ax.axvline(crit, color='red', ls='--', lw=1, label=f'chi2(1) 5% = {crit:.2f}') ax.set_xlabel('Likelihood ratio statistic') ax.set_title('(c) LR test: q-exponential vs exponential') ax.legend(fontsize=8) fig.tight_layout() paths = save_figure(fig, 'tsallis_tails_test', subfolder='potential') # --- Figure 2: GARCH robustness comparison --- garch_res = _cache.get('garch_results', []) if garch_res: garch_df = pd.DataFrame(garch_res) raw_df = df[['series', 'lr_stat', 'lr_pval']].rename( columns={'lr_stat': 'raw_lr', 'lr_pval': 'raw_pval'}) merged = raw_df.merge(garch_df.rename( columns={'lr_stat': 'garch_lr', 'lr_pval': 'garch_pval'}), on='series') merged = merged.sort_values('raw_lr', ascending=True) merged['label'] = [LEAF_SERIES.get(s, s) for s in merged['series']] fig2, ax2 = setup_figure() y_pos = np.arange(len(merged)) bar_h = 0.35 ax2.barh(y_pos + bar_h / 2, merged['raw_lr'], bar_h, color='steelblue', alpha=0.8, label='Raw returns') ax2.barh(y_pos - bar_h / 2, merged['garch_lr'], bar_h, color='darkorange', alpha=0.8, label='GARCH residuals') crit = stats.chi2.ppf(0.95, df=1) ax2.axvline(crit, color='red', ls='--', lw=1, label=f'5% critical = {crit:.2f}') ax2.set_yticks(y_pos) ax2.set_yticklabels(merged['label'], fontsize=7) ax2.set_xlabel('Likelihood ratio statistic') ax2.set_title('GARCH Robustness: Raw vs Standardized Residuals') ax2.legend(fontsize=8) ax2.grid(True, alpha=0.3, axis='x') fig2.tight_layout() paths.extend(save_figure(fig2, 'tsallis_tails_garch', subfolder='potential')) return paths if __name__ == '__main__': results = main() for r in results: print(f"{r['name']}: {r['verdict']} " f"(stat={r['statistic']}, p={r['p_value']})")

N	17
R2	0.0323
slope	-0.2985
p value	0.4898
statistic	-0.2985

Tsallis Tails Q Rho Regression

Description

Figures

Statistics

Script

▸ Source Code

Reproduction Steps

Data Sources