Session 5: Statistical Inference II - Non-Parametric Tests

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import yfinance as yf
import statsmodels.api as sm
from datetime import timedelta
from dataclasses import dataclass
from typing import List, Dict, Optional
import warnings
warnings.filterwarnings('ignore')

plt.style.use('seaborn-v0_8-whitegrid')
print("Libraries loaded!")

# Sample events
EVENTS = [
    {'ticker': 'AAPL', 'date': '2023-08-03', 'name': 'Apple Q3'},
    {'ticker': 'MSFT', 'date': '2023-07-25', 'name': 'Microsoft Q4'},
    {'ticker': 'GOOGL', 'date': '2023-07-25', 'name': 'Alphabet Q2'},
    {'ticker': 'AMZN', 'date': '2023-08-03', 'name': 'Amazon Q2'},
    {'ticker': 'META', 'date': '2023-07-26', 'name': 'Meta Q2'},
    {'ticker': 'NVDA', 'date': '2023-08-23', 'name': 'Nvidia Q2'},
    {'ticker': 'TSLA', 'date': '2023-07-19', 'name': 'Tesla Q2'},
    {'ticker': 'AMD', 'date': '2023-08-01', 'name': 'AMD Q2'},
    {'ticker': 'INTC', 'date': '2023-07-27', 'name': 'Intel Q2'},
    {'ticker': 'CRM', 'date': '2023-08-30', 'name': 'Salesforce Q2'},
]

EST_WINDOW, GAP, PRE, POST = 120, 10, 10, 10
WINDOWS = [(-1, +1), (0, 0), (-5, +5), (-10, +10)]

@dataclass
class EventResult:
    ticker: str
    model_sigma: float
    model_n_obs: int
    model_market_mean: float
    model_market_var: float
    event_data: pd.DataFrame
    estimation_data: pd.DataFrame

def process_event(ticker, event_date, est_window, gap, pre, post):
    """Download and process single event."""
    try:
        event_dt = pd.to_datetime(event_date)
        start = event_dt - timedelta(days=int((est_window + gap + pre) * 1.5))
        end = event_dt + timedelta(days=int(post * 2.5))
        
        stock = yf.download(ticker, start=start, end=end, progress=False)['Close']
        market = yf.download('^GSPC', start=start, end=end, progress=False)['Close']
        
        df = pd.DataFrame({'stock': stock.squeeze(), 'market': market.squeeze()})
        df['stock_ret'] = df['stock'].pct_change()
        df['market_ret'] = df['market'].pct_change()
        df = df.dropna()
        
        if event_dt not in df.index:
            idx = df.index.get_indexer([event_dt], method='nearest')[0]
            event_dt = df.index[idx]
        
        event_idx = df.index.get_loc(event_dt)
        df['event_time'] = range(-event_idx, len(df) - event_idx)
        
        est_end = -(gap + pre)
        est_data = df[(df['event_time'] >= est_end - est_window) & (df['event_time'] < est_end)].copy()
        evt_data = df[(df['event_time'] >= -pre) & (df['event_time'] <= post)].copy()
        
        # Market model
        y, x = est_data['stock_ret'].values, est_data['market_ret'].values
        X = sm.add_constant(x)
        ols = sm.OLS(y, X).fit()
        
        est_data['AR'] = est_data['stock_ret'] - (ols.params[0] + ols.params[1] * est_data['market_ret'])
        evt_data['AR'] = evt_data['stock_ret'] - (ols.params[0] + ols.params[1] * evt_data['market_ret'])
        
        return EventResult(
            ticker=ticker, model_sigma=np.std(ols.resid, ddof=2),
            model_n_obs=len(y), model_market_mean=np.mean(x),
            model_market_var=np.sum((x - np.mean(x))**2),
            event_data=evt_data, estimation_data=est_data
        )
    except Exception as e:
        print(f"{ticker}: FAILED - {e}")
        return None

print("Processing events...")
event_results = [r for e in EVENTS if (r := process_event(e['ticker'], e['date'], EST_WINDOW, GAP, PRE, POST))]
print(f"Processed {len(event_results)} events")

def calculate_car(event_data, tau1, tau2):
    mask = (event_data['event_time'] >= tau1) & (event_data['event_time'] <= tau2)
    return event_data.loc[mask, 'AR'].sum()

all_ars = np.concatenate([r.estimation_data['AR'].values for r in event_results])
cars_11 = np.array([calculate_car(r.event_data, -1, 1) for r in event_results])

print("Distribution of Estimation Period Residuals:")
print(f"  Skewness: {stats.skew(all_ars):.3f}")
print(f"  Kurtosis: {stats.kurtosis(all_ars):.3f} (excess)")
jb_stat, jb_p = stats.jarque_bera(all_ars)
print(f"  Jarque-Bera p-value: {jb_p:.4f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].hist(all_ars*100, bins=50, density=True, alpha=0.7, edgecolor='black')
axes[0].set_title('Estimation Period ARs')
axes[0].set_xlabel('AR (%)')
axes[1].hist(cars_11*100, bins=8, density=True, alpha=0.7, edgecolor='black')
axes[1].axvline(0, color='red', linestyle='--')
axes[1].set_title('CAR(-1,+1) Distribution')
axes[1].set_xlabel('CAR (%)')
plt.tight_layout()
plt.show()

def sign_test(event_results, tau1, tau2):
    cars = np.array([calculate_car(r.event_data, tau1, tau2) for r in event_results])
    N = len(cars)
    N_pos = np.sum(cars > 0)
    Z = (N_pos - 0.5 * N) / np.sqrt(0.25 * N)
    p_exact = stats.binomtest(N_pos, N, 0.5)
    return {'N': N, 'N_pos': N_pos, 'prop_pos': N_pos/N, 'Z': Z, 'p_value': p_exact, 'CAAR': np.mean(cars)}

print("Sign Test Results:")
print("="*60)
for tau1, tau2 in WINDOWS:
    r = sign_test(event_results, tau1, tau2)
    r['p_value'] = r['p_value'].pvalue
    sig = '***' if r['p_value'] < 0.01 else '**' if r['p_value'] < 0.05 else '*' if r['p_value'] < 0.10 else ''
    print(f"[{tau1:+d},{tau2:+d}]: {r['N_pos']}/{r['N']} positive ({r['prop_pos']*100:.0f}%), Z={r['Z']:+.3f}, p={r['p_value']:.4f} {sig}")

def generalized_sign_test(event_results, tau1, tau2):
    cars = np.array([calculate_car(r.event_data, tau1, tau2) for r in event_results])
    N = len(cars)
    
    # Estimate p_hat from estimation periods
    p_hat = np.mean([np.mean(r.estimation_data['AR'].values > 0) for r in event_results])
    
    N_pos = np.sum(cars > 0)
    expected = N * p_hat
    Z = (N_pos - expected) / np.sqrt(N * p_hat * (1 - p_hat))
    p_value = 2 * (1 - stats.norm.cdf(abs(Z)))
    
    return {'N': N, 'N_pos': N_pos, 'p_hat': p_hat, 'expected': expected, 'Z': Z, 'p_value': p_value}

print("\nGeneralized Sign Test Results:")
print("="*60)
p_hat = np.mean([np.mean(r.estimation_data['AR'].values > 0) for r in event_results])
print(f"Estimated baseline P(AR>0) = {p_hat:.3f}")
for tau1, tau2 in WINDOWS:
    r = generalized_sign_test(event_results, tau1, tau2)
    sig = '***' if r['p_value'] < 0.01 else '**' if r['p_value'] < 0.05 else '*' if r['p_value'] < 0.10 else ''
    print(f"[{tau1:+d},{tau2:+d}]: {r['N_pos']}/{r['N']} vs {r['expected']:.1f} expected, Z={r['Z']:+.3f}, p={r['p_value']:.4f} {sig}")

def corrado_rank_test(event_results, tau1, tau2):
    N = len(event_results)
    transformed_ranks = []
    all_K = []
    
    for r in event_results:
        all_ars = np.concatenate([r.estimation_data['AR'].values, r.event_data['AR'].values])
        T = len(all_ars)
        ranks = stats.rankdata(all_ars, method='average')
        K = ranks / (T + 1) - 0.5
        
        all_K.extend(K[:len(r.estimation_data)])
        
        evt_K = K[len(r.estimation_data):]
        mask = (r.event_data['event_time'].values >= tau1) & (r.event_data['event_time'].values <= tau2)
        transformed_ranks.append(np.mean(evt_K[mask]))
    
    K_bar = np.mean(transformed_ranks)
    S_K = np.std(all_K, ddof=1) / np.sqrt(N)
    Z = K_bar / S_K if S_K > 0 else 0
    p_value = 2 * (1 - stats.norm.cdf(abs(Z)))
    
    return {'N': N, 'mean_rank': K_bar, 'Z': Z, 'p_value': p_value}

print("\nCorrado (1989) Rank Test Results:")
print("="*60)
for tau1, tau2 in WINDOWS:
    r = corrado_rank_test(event_results, tau1, tau2)
    sig = '***' if r['p_value'] < 0.01 else '**' if r['p_value'] < 0.05 else '*' if r['p_value'] < 0.10 else ''
    print(f"[{tau1:+d},{tau2:+d}]: Mean Rank={r['mean_rank']:+.4f}, Z={r['Z']:+.3f}, p={r['p_value']:.4f} {sig}")

def wilcoxon_test(event_results, tau1, tau2):
    cars = np.array([calculate_car(r.event_data, tau1, tau2) for r in event_results])
    stat, p_value = stats.wilcoxon(cars, alternative='two-sided')
    return {'N': len(cars), 'median': np.median(cars), 'W': stat, 'p_value': p_value}

print("\nWilcoxon Signed-Rank Test Results:")
print("="*60)
for tau1, tau2 in WINDOWS:
    r = wilcoxon_test(event_results, tau1, tau2)
    sig = '***' if r['p_value'] < 0.01 else '**' if r['p_value'] < 0.05 else '*' if r['p_value'] < 0.10 else ''
    print(f"[{tau1:+d},{tau2:+d}]: Median={r['median']*100:+.3f}%, W={r['W']:.1f}, p={r['p_value']:.4f} {sig}")

def bootstrap_test(event_results, tau1, tau2, n_boot=5000):
    cars = np.array([calculate_car(r.event_data, tau1, tau2) for r in event_results])
    N = len(cars)
    caar = np.mean(cars)
    
    np.random.seed(42)
    boot_caars = [np.mean(cars[np.random.choice(N, N, replace=True)]) for _ in range(n_boot)]
    boot_caars = np.array(boot_caars)
    
    boot_se = np.std(boot_caars)
    ci_lower, ci_upper = np.percentile(boot_caars, [2.5, 97.5])
    p_value = np.mean(np.abs(boot_caars - caar) >= np.abs(caar))
    t_stat = caar / boot_se if boot_se > 0 else 0
    
    return {'N': N, 'CAAR': caar, 'boot_SE': boot_se, 't': t_stat, 
            'CI_lower': ci_lower, 'CI_upper': ci_upper, 'p_value': p_value, 'distribution': boot_caars}

print("\nBootstrap Test Results (B=5000):")
print("="*70)
boot_results = []
for tau1, tau2 in WINDOWS:
    r = bootstrap_test(event_results, tau1, tau2)
    boot_results.append(r)
    sig = '***' if r['p_value'] < 0.01 else '**' if r['p_value'] < 0.05 else '*' if r['p_value'] < 0.10 else ''
    print(f"[{tau1:+d},{tau2:+d}]: CAAR={r['CAAR']*100:+.3f}%, 95% CI=[{r['CI_lower']*100:+.2f}%,{r['CI_upper']*100:+.2f}%], t={r['t']:+.3f}, p={r['p_value']:.4f} {sig}")

# Visualize bootstrap distributions
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, ((tau1, tau2), r) in enumerate(zip(WINDOWS, boot_results)):
    ax = axes.flatten()[i]
    ax.hist(r['distribution']*100, bins=50, density=True, alpha=0.7, edgecolor='black')
    ax.axvline(r['CAAR']*100, color='red', lw=2, label=f'CAAR={r["CAAR"]*100:.2f}%')
    ax.axvline(r['CI_lower']*100, color='green', ls='--', label='95% CI')
    ax.axvline(r['CI_upper']*100, color='green', ls='--')
    ax.axvline(0, color='black', ls=':')
    ax.set_title(f'Window [{tau1:+d},{tau2:+d}]')
    ax.set_xlabel('CAAR (%)')
    ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

# Add parametric tests for comparison
def t_test(event_results, tau1, tau2):
    cars = np.array([calculate_car(r.event_data, tau1, tau2) for r in event_results])
    N = len(cars)
    t_stat = np.mean(cars) / (np.std(cars, ddof=1) / np.sqrt(N))
    p_value = 2 * (1 - stats.t.cdf(abs(t_stat), df=N-1))
    return {'t': t_stat, 'p': p_value}

print("\n" + "="*100)
print("COMPREHENSIVE TEST COMPARISON")
print("="*100)

for tau1, tau2 in WINDOWS:
    cars = [calculate_car(r.event_data, tau1, tau2) for r in event_results]
    print(f"\nWindow [{tau1:+d},{tau2:+d}]: CAAR={np.mean(cars)*100:+.3f}%, Median={np.median(cars)*100:+.3f}%")
    print("-"*80)
    
    tt = t_test(event_results, tau1, tau2)
    st = sign_test(event_results, tau1, tau2)
    gst = generalized_sign_test(event_results, tau1, tau2)
    cr = corrado_rank_test(event_results, tau1, tau2)
    wt = wilcoxon_test(event_results, tau1, tau2)
    bt = bootstrap_test(event_results, tau1, tau2, n_boot=2000)
    
    def sig(p): return '***' if p < 0.01 else '**' if p < 0.05 else '*' if p < 0.10 else ''
    
    results = [
        ('t-test (Param)', tt['t'], tt['p']),
        ('Sign Test', st['Z'], st['p_value'].pvalue),
        ('Gen. Sign Test', gst['Z'], gst['p_value']),
        ('Corrado Rank', cr['Z'], cr['p_value']),
        ('Wilcoxon', np.nan, wt['p_value']),
        ('Bootstrap', bt['t'], bt['p_value']),
    ]
    
    for name, stat, p in results:
        stat_str = f"{stat:+.3f}" if not np.isnan(stat) else "  N/A"
        print(f"  {name:18s}: stat={stat_str}, p={p:.4f} {sig(p)}")

# Compare p-values visually
fig, ax = plt.subplots(figsize=(12, 5))

test_names = ['t-test', 'Sign', 'GST', 'Corrado', 'Wilcoxon', 'Bootstrap']
x = np.arange(len(WINDOWS))
width = 0.12

for j, name in enumerate(test_names):
    p_vals = []
    for tau1, tau2 in WINDOWS:
        if name == 't-test':
            p_vals.append(t_test(event_results, tau1, tau2)['p'])
        elif name == 'Sign':
            p_vals.append(sign_test(event_results, tau1, tau2)['p_value'].pvalue)
        elif name == 'GST':
            p_vals.append(generalized_sign_test(event_results, tau1, tau2)['p_value'])
        elif name == 'Corrado':
            p_vals.append(corrado_rank_test(event_results, tau1, tau2)['p_value'])
        elif name == 'Wilcoxon':
            p_vals.append(wilcoxon_test(event_results, tau1, tau2)['p_value'])
        else:
            p_vals.append(bootstrap_test(event_results, tau1, tau2, 1000)['p_value'])
    ax.bar(x + j*width, p_vals, width, label=name, alpha=0.8)

ax.axhline(0.05, color='red', ls='--', label='5% level')
ax.axhline(0.10, color='orange', ls=':', label='10% level')
ax.set_xticks(x + 2.5*width)
ax.set_xticklabels([f'[{t1:+d},{t2:+d}]' for t1, t2 in WINDOWS])
ax.set_ylabel('p-value')
ax.set_title('Comparison of p-values Across Tests')
ax.legend(bbox_to_anchor=(1.02, 1))
ax.set_ylim(0, 1)
plt.tight_layout()
plt.show()

def monte_carlo_power(N_events, N_sims, true_car, skewness=0, alpha=0.05):
    """Monte Carlo for test size/power."""
    rejections = {'t': 0, 'sign': 0, 'wilcoxon': 0}
    
    for _ in range(N_sims):
        if skewness == 0:
            cars = np.random.normal(true_car, 0.05, N_events)
        else:
            cars = stats.skewnorm.rvs(skewness, loc=true_car, scale=0.05, size=N_events)
        
        _, t_p = stats.ttest_1samp(cars, 0)
        if t_p < alpha: rejections['t'] += 1
        
        sign_p = stats.binomtest(np.sum(cars > 0), N_events, 0.5).pvalue
        if sign_p < alpha: rejections['sign'] += 1
        
        try:
            _, wilc_p = stats.wilcoxon(cars)
            if wilc_p < alpha: rejections['wilcoxon'] += 1
        except: pass
    
    return {k: v/N_sims for k, v in rejections.items()}

print("Monte Carlo Analysis (N=30, 500 sims):")
print("="*50)

print("\n1. SIZE (CAR=0, Normal):")
size = monte_carlo_power(30, 500, 0, 0)
for t, r in size.items(): print(f"   {t}: {r*100:.1f}%")

print("\n2. SIZE (CAR=0, Skewed):")
size_skew = monte_carlo_power(30, 500, 0, 3)
for t, r in size_skew.items(): print(f"   {t}: {r*100:.1f}%")

print("\n3. POWER (CAR=1%):")
power = monte_carlo_power(30, 500, 0.01, 0)
for t, r in power.items(): print(f"   {t}: {r*100:.1f}%")

# Publication table
print("\n" + "="*90)
print("PUBLICATION TABLE")
print("="*90)

def sig(p): return '***' if p < 0.01 else '**' if p < 0.05 else '*' if p < 0.10 else ''

rows = []
for tau1, tau2 in WINDOWS:
    cars = [calculate_car(r.event_data, tau1, tau2) for r in event_results]
    tt = t_test(event_results, tau1, tau2)
    gst = generalized_sign_test(event_results, tau1, tau2)
    cr = corrado_rank_test(event_results, tau1, tau2)
    
    rows.append({
        'Window': f'[{tau1:+d},{tau2:+d}]',
        'N': len(cars),
        'CAAR': f"{np.mean(cars)*100:+.2f}%",
        'Median': f"{np.median(cars)*100:+.2f}%",
        '%Pos': f"{np.mean(np.array(cars)>0)*100:.0f}%",
        't-stat': f"{tt['t']:.2f}{sig(tt['p'])}",
        'GST': f"{gst['Z']:.2f}{sig(gst['p_value'])}",
        'Rank': f"{cr['Z']:.2f}{sig(cr['p_value'])}"
    })

print(pd.DataFrame(rows).to_string(index=False))
print("\nNotes: ***, **, * = 1%, 5%, 10% significance")