Session 6: Cross-Sectional Analysis

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 statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan, het_white
from datetime import timedelta
from dataclasses import dataclass
from typing import List, Dict, Optional, Tuple
import warnings
warnings.filterwarnings('ignore')

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

# Extended sample with firm characteristics for cross-sectional analysis
# Using Q2-Q3 2023 earnings announcements with additional firm data
EVENTS = [
    # Tech - Large Cap
    {'ticker': 'AAPL', 'date': '2023-08-03', 'name': 'Apple', 'sector': 'Tech', 'eps_surprise': 0.05},
    {'ticker': 'MSFT', 'date': '2023-07-25', 'name': 'Microsoft', 'sector': 'Tech', 'eps_surprise': 0.08},
    {'ticker': 'GOOGL', 'date': '2023-07-25', 'name': 'Alphabet', 'sector': 'Tech', 'eps_surprise': 0.12},
    {'ticker': 'AMZN', 'date': '2023-08-03', 'name': 'Amazon', 'sector': 'Tech', 'eps_surprise': 0.45},
    {'ticker': 'META', 'date': '2023-07-26', 'name': 'Meta', 'sector': 'Tech', 'eps_surprise': 0.22},
    {'ticker': 'NVDA', 'date': '2023-08-23', 'name': 'Nvidia', 'sector': 'Tech', 'eps_surprise': 0.35},
    # Tech - Smaller
    {'ticker': 'AMD', 'date': '2023-08-01', 'name': 'AMD', 'sector': 'Tech', 'eps_surprise': -0.02},
    {'ticker': 'INTC', 'date': '2023-07-27', 'name': 'Intel', 'sector': 'Tech', 'eps_surprise': 0.18},
    {'ticker': 'CRM', 'date': '2023-08-30', 'name': 'Salesforce', 'sector': 'Tech', 'eps_surprise': 0.15},
    {'ticker': 'ADBE', 'date': '2023-09-14', 'name': 'Adobe', 'sector': 'Tech', 'eps_surprise': 0.06},
    # Finance
    {'ticker': 'JPM', 'date': '2023-07-14', 'name': 'JPMorgan', 'sector': 'Finance', 'eps_surprise': 0.10},
    {'ticker': 'BAC', 'date': '2023-07-18', 'name': 'Bank of America', 'sector': 'Finance', 'eps_surprise': 0.03},
    {'ticker': 'WFC', 'date': '2023-07-14', 'name': 'Wells Fargo', 'sector': 'Finance', 'eps_surprise': 0.08},
    {'ticker': 'GS', 'date': '2023-07-19', 'name': 'Goldman Sachs', 'sector': 'Finance', 'eps_surprise': -0.15},
    {'ticker': 'MS', 'date': '2023-07-18', 'name': 'Morgan Stanley', 'sector': 'Finance', 'eps_surprise': -0.05},
    # Healthcare
    {'ticker': 'JNJ', 'date': '2023-07-20', 'name': 'Johnson & Johnson', 'sector': 'Healthcare', 'eps_surprise': 0.04},
    {'ticker': 'UNH', 'date': '2023-07-14', 'name': 'UnitedHealth', 'sector': 'Healthcare', 'eps_surprise': 0.02},
    {'ticker': 'PFE', 'date': '2023-08-01', 'name': 'Pfizer', 'sector': 'Healthcare', 'eps_surprise': -0.08},
    {'ticker': 'MRK', 'date': '2023-07-27', 'name': 'Merck', 'sector': 'Healthcare', 'eps_surprise': 0.12},
    {'ticker': 'ABBV', 'date': '2023-07-27', 'name': 'AbbVie', 'sector': 'Healthcare', 'eps_surprise': 0.05},
    # Consumer
    {'ticker': 'PG', 'date': '2023-07-28', 'name': 'Procter & Gamble', 'sector': 'Consumer', 'eps_surprise': 0.03},
    {'ticker': 'KO', 'date': '2023-07-26', 'name': 'Coca-Cola', 'sector': 'Consumer', 'eps_surprise': 0.05},
    {'ticker': 'PEP', 'date': '2023-07-13', 'name': 'PepsiCo', 'sector': 'Consumer', 'eps_surprise': 0.07},
    {'ticker': 'WMT', 'date': '2023-08-17', 'name': 'Walmart', 'sector': 'Consumer', 'eps_surprise': 0.09},
    {'ticker': 'COST', 'date': '2023-09-26', 'name': 'Costco', 'sector': 'Consumer', 'eps_surprise': 0.04},
]

EST_WINDOW, GAP, PRE, POST = 120, 10, 5, 5
print(f"Sample: {len(EVENTS)} earnings announcements across {len(set(e['sector'] for e in EVENTS))} sectors")

@dataclass
class EventResult:
    """Container for event study results with firm characteristics."""
    ticker: str
    name: str
    sector: str
    event_date: pd.Timestamp
    eps_surprise: float
    market_cap: float
    beta: float
    sigma: float
    pre_event_return: float
    event_data: pd.DataFrame
    car_var: float  # Variance of CAR for WLS


def get_firm_characteristics(ticker: str) -> Dict:
    """Fetch firm characteristics from Yahoo Finance."""
    try:
        stock = yf.Ticker(ticker)
        info = stock.info
        return {
            'market_cap': info.get('marketCap', np.nan),
            'beta': info.get('beta', np.nan),
        }
    except:
        return {'market_cap': np.nan, 'beta': np.nan}


def process_event(event: Dict, est_window: int, gap: int, pre: int, post: int) -> Optional[EventResult]:
    """Process single event and extract characteristics."""
    try:
        ticker = event['ticker']
        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)
        
        # Split windows
        est_end = -(gap + pre)
        est_data = df[(df['event_time'] >= est_end - est_window) & (df['event_time'] < est_end)]
        evt_data = df[(df['event_time'] >= -pre) & (df['event_time'] <= post)].copy()
        
        # Market model estimation
        y, x = est_data['stock_ret'].values, est_data['market_ret'].values
        X = sm.add_constant(x)
        ols = sm.OLS(y, X).fit()
        
        alpha, beta = ols.params[0], ols.params[1]
        sigma = np.std(ols.resid, ddof=2)
        
        # Abnormal returns
        evt_data['expected'] = alpha + beta * evt_data['market_ret']
        evt_data['AR'] = evt_data['stock_ret'] - evt_data['expected']
        
        # Pre-event cumulative return (momentum proxy)
        pre_event_data = df[(df['event_time'] >= -60) & (df['event_time'] < -pre)]
        pre_event_return = (1 + pre_event_data['stock_ret']).prod() - 1 if len(pre_event_data) > 0 else 0
        
        # CAR variance for WLS
        L = len(evt_data[(evt_data['event_time'] >= -1) & (evt_data['event_time'] <= 1)])
        car_var = L * sigma**2
        
        # Get firm characteristics
        chars = get_firm_characteristics(ticker)
        
        return EventResult(
            ticker=ticker,
            name=event['name'],
            sector=event['sector'],
            event_date=event_dt,
            eps_surprise=event['eps_surprise'],
            market_cap=chars['market_cap'],
            beta=beta,
            sigma=sigma,
            pre_event_return=pre_event_return,
            event_data=evt_data,
            car_var=car_var
        )
    except Exception as e:
        print(f"{event['ticker']}: FAILED - {e}")
        return None


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

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


# Create cross-sectional dataset
def create_cross_sectional_data(event_results: List[EventResult], tau1: int, tau2: int) -> pd.DataFrame:
    """Create dataset for cross-sectional regression."""
    data = []
    for r in event_results:
        car = calculate_car(r.event_data, tau1, tau2)
        data.append({
            'ticker': r.ticker,
            'name': r.name,
            'sector': r.sector,
            'CAR': car,
            'eps_surprise': r.eps_surprise,
            'log_mcap': np.log(r.market_cap) if r.market_cap and r.market_cap > 0 else np.nan,
            'beta': r.beta,
            'sigma': r.sigma,
            'momentum': r.pre_event_return,
            'car_var': r.car_var,
            'car_se': np.sqrt(r.car_var)
        })
    
    df = pd.DataFrame(data)
    df = df.dropna()
    return df


# Create dataset for CAR(-1,+1)
cs_data = create_cross_sectional_data(event_results, -1, 1)
print(f"\nCross-sectional dataset: {len(cs_data)} observations")
print(cs_data[['ticker', 'sector', 'CAR', 'eps_surprise', 'log_mcap']].head(10).to_string())

# Summary statistics by sector
print("Summary Statistics by Sector:")
print("="*70)

sector_stats = cs_data.groupby('sector').agg({
    'CAR': ['count', 'mean', 'median', 'std'],
    'eps_surprise': 'mean',
    'log_mcap': 'mean'
}).round(4)

sector_stats.columns = ['N', 'Mean CAR', 'Median CAR', 'Std CAR', 'Avg EPS Surp', 'Avg Log MCap']
sector_stats['Mean CAR'] = sector_stats['Mean CAR'].apply(lambda x: f"{x*100:.2f}%")
sector_stats['Median CAR'] = sector_stats['Median CAR'].apply(lambda x: f"{x*100:.2f}%")
sector_stats['Std CAR'] = sector_stats['Std CAR'].apply(lambda x: f"{x*100:.2f}%")
print(sector_stats.to_string())

# Visualize relationships
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# CAR vs EPS Surprise
ax1 = axes[0, 0]
colors = {'Tech': 'blue', 'Finance': 'green', 'Healthcare': 'red', 'Consumer': 'purple'}
for sector in cs_data['sector'].unique():
    mask = cs_data['sector'] == sector
    ax1.scatter(cs_data.loc[mask, 'eps_surprise'], cs_data.loc[mask, 'CAR']*100, 
                label=sector, alpha=0.7, s=80, c=colors.get(sector, 'gray'))

# Add regression line
z = np.polyfit(cs_data['eps_surprise'], cs_data['CAR']*100, 1)
p = np.poly1d(z)
x_line = np.linspace(cs_data['eps_surprise'].min(), cs_data['eps_surprise'].max(), 100)
ax1.plot(x_line, p(x_line), 'k--', alpha=0.7, label='OLS fit')

ax1.set_xlabel('EPS Surprise', fontsize=11)
ax1.set_ylabel('CAR (-1,+1) %', fontsize=11)
ax1.set_title('CAR vs Earnings Surprise', fontsize=12)
ax1.axhline(0, color='gray', linestyle='-', linewidth=0.5)
ax1.axvline(0, color='gray', linestyle='-', linewidth=0.5)
ax1.legend(fontsize=9)

# CAR vs Market Cap
ax2 = axes[0, 1]
ax2.scatter(cs_data['log_mcap'], cs_data['CAR']*100, alpha=0.7, s=80)
z = np.polyfit(cs_data['log_mcap'], cs_data['CAR']*100, 1)
p = np.poly1d(z)
x_line = np.linspace(cs_data['log_mcap'].min(), cs_data['log_mcap'].max(), 100)
ax2.plot(x_line, p(x_line), 'r--', alpha=0.7)
ax2.set_xlabel('Log Market Cap', fontsize=11)
ax2.set_ylabel('CAR (-1,+1) %', fontsize=11)
ax2.set_title('CAR vs Firm Size', fontsize=12)
ax2.axhline(0, color='gray', linestyle='-', linewidth=0.5)

# CAR distribution by sector
ax3 = axes[1, 0]
sector_order = cs_data.groupby('sector')['CAR'].mean().sort_values().index
cs_data['sector_cat'] = pd.Categorical(cs_data['sector'], categories=sector_order, ordered=True)
bp = ax3.boxplot([cs_data[cs_data['sector'] == s]['CAR']*100 for s in sector_order],
                  labels=sector_order, patch_artist=True)
for patch, sector in zip(bp['boxes'], sector_order):
    patch.set_facecolor(colors.get(sector, 'gray'))
    patch.set_alpha(0.6)
ax3.axhline(0, color='red', linestyle='--', alpha=0.7)
ax3.set_ylabel('CAR (-1,+1) %', fontsize=11)
ax3.set_title('CAR Distribution by Sector', fontsize=12)

# CAR vs Pre-event Momentum
ax4 = axes[1, 1]
ax4.scatter(cs_data['momentum']*100, cs_data['CAR']*100, alpha=0.7, s=80)
z = np.polyfit(cs_data['momentum']*100, cs_data['CAR']*100, 1)
p = np.poly1d(z)
x_line = np.linspace(cs_data['momentum'].min()*100, cs_data['momentum'].max()*100, 100)
ax4.plot(x_line, p(x_line), 'r--', alpha=0.7)
ax4.set_xlabel('Pre-Event 60-day Return (%)', fontsize=11)
ax4.set_ylabel('CAR (-1,+1) %', fontsize=11)
ax4.set_title('CAR vs Momentum', fontsize=12)
ax4.axhline(0, color='gray', linestyle='-', linewidth=0.5)
ax4.axvline(0, color='gray', linestyle='-', linewidth=0.5)

plt.tight_layout()
plt.show()

# Correlation matrix
corr_vars = ['CAR', 'eps_surprise', 'log_mcap', 'beta', 'sigma', 'momentum']
corr_matrix = cs_data[corr_vars].corr()

fig, ax = plt.subplots(figsize=(8, 6))
sns.heatmap(corr_matrix, annot=True, fmt='.2f', cmap='RdBu_r', center=0,
            vmin=-1, vmax=1, ax=ax, square=True)
ax.set_title('Correlation Matrix of Cross-Sectional Variables', fontsize=12)
plt.tight_layout()
plt.show()

def run_cross_sectional_regression(data: pd.DataFrame, 
                                    y_var: str, 
                                    x_vars: List[str],
                                    add_constant: bool = True) -> sm.regression.linear_model.RegressionResultsWrapper:
    """Run OLS cross-sectional regression."""
    y = data[y_var]
    X = data[x_vars]
    
    if add_constant:
        X = sm.add_constant(X)
    
    model = sm.OLS(y, X).fit()
    return model


# Model 1: Simple - EPS Surprise only
print("="*80)
print("MODEL 1: CAR on EPS Surprise")
print("="*80)

model1 = run_cross_sectional_regression(cs_data, 'CAR', ['eps_surprise'])
print(model1.summary())

# Model 2: Add firm characteristics
print("\n" + "="*80)
print("MODEL 2: CAR on EPS Surprise + Firm Characteristics")
print("="*80)

model2 = run_cross_sectional_regression(cs_data, 'CAR', 
                                         ['eps_surprise', 'log_mcap', 'momentum'])
print(model2.summary())

# Model 3: Add sector fixed effects
print("\n" + "="*80)
print("MODEL 3: With Sector Fixed Effects")
print("="*80)

# Create sector dummies
cs_data_fe = cs_data.copy()
sector_dummies = pd.get_dummies(cs_data_fe['sector'], prefix='sector', drop_first=True)
cs_data_fe = pd.concat([cs_data_fe, sector_dummies], axis=1)

x_vars_fe = ['eps_surprise', 'log_mcap', 'momentum'] + list(sector_dummies.columns)
model3 = run_cross_sectional_regression(cs_data_fe, 'CAR', x_vars_fe)
print(model3.summary())

 ['eps_surprise', 'log_mcap', 'momentum'] + list(sector_dummies.columns)

# Test for heteroskedasticity
print("Heteroskedasticity Tests:")
print("="*60)

# Get residuals from Model 2
y = cs_data['CAR']
X = sm.add_constant(cs_data[['eps_surprise', 'log_mcap', 'momentum']])

# Breusch-Pagan test
bp_lm, bp_pval, bp_f, bp_fpval = het_breuschpagan(model2.resid, X)
print(f"\nBreusch-Pagan Test:")
print(f"  LM statistic: {bp_lm:.4f}")
print(f"  p-value: {bp_pval:.4f}")
print(f"  Conclusion: {'Reject H0 - Heteroskedasticity present' if bp_pval < 0.05 else 'Fail to reject H0'}")

# White test
white_lm, white_pval, white_f, white_fpval = het_white(model2.resid, X)
print(f"\nWhite Test:")
print(f"  LM statistic: {white_lm:.4f}")
print(f"  p-value: {white_pval:.4f}")
print(f"  Conclusion: {'Reject H0 - Heteroskedasticity present' if white_pval < 0.05 else 'Fail to reject H0'}")

# Visualize heteroskedasticity
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Residuals vs Fitted
ax1 = axes[0]
ax1.scatter(model2.fittedvalues*100, model2.resid*100, alpha=0.7)
ax1.axhline(0, color='red', linestyle='--')
ax1.set_xlabel('Fitted Values (%)', fontsize=11)
ax1.set_ylabel('Residuals (%)', fontsize=11)
ax1.set_title('Residuals vs Fitted Values', fontsize=12)

# Residuals vs Sigma (proxy for heteroskedasticity)
ax2 = axes[1]
ax2.scatter(cs_data['sigma']*100, np.abs(model2.resid)*100, alpha=0.7)
z = np.polyfit(cs_data['sigma']*100, np.abs(model2.resid)*100, 1)
p = np.poly1d(z)
x_line = np.linspace(cs_data['sigma'].min()*100, cs_data['sigma'].max()*100, 100)
ax2.plot(x_line, p(x_line), 'r--', alpha=0.7)
ax2.set_xlabel('Firm Volatility (σ) %', fontsize=11)
ax2.set_ylabel('|Residuals| (%)', fontsize=11)
ax2.set_title('Absolute Residuals vs Volatility', fontsize=12)

plt.tight_layout()
plt.show()

print("\nNote: If residuals spread out with volatility, we have heteroskedasticity.")

# Compare OLS with Robust Standard Errors
print("\n" + "="*80)
print("COMPARISON: OLS vs Robust Standard Errors")
print("="*80)

# OLS with robust SEs (HC1 = White's heteroskedasticity-consistent)
model2_robust = run_cross_sectional_regression(cs_data, 'CAR', 
                                                ['eps_surprise', 'log_mcap', 'momentum'])
model2_robust_hc = model2_robust.get_robustcov_results(cov_type='HC1')

# Create comparison table
comparison = pd.DataFrame({
    'Variable': model2.params.index,
    'Coef': model2.params.values,
    'OLS SE': model2.bse.values,
    'OLS t-stat': model2.tvalues.values,
    'Robust SE': model2_robust_hc.bse,
    'Robust t-stat': model2_robust_hc.tvalues,
})

comparison['SE Change'] = (comparison['Robust SE'] / comparison['OLS SE'] - 1) * 100

print("\nCoefficient Comparison:")
for col in ['Coef', 'OLS SE', 'Robust SE']:
    comparison[col] = comparison[col].apply(lambda x: f"{x:.6f}")
for col in ['OLS t-stat', 'Robust t-stat']:
    comparison[col] = comparison[col].apply(lambda x: f"{x:.3f}")
comparison['SE Change'] = comparison['SE Change'].apply(lambda x: f"{x:+.1f}%")

print(comparison.to_string(index=False))
print("\nNote: Positive SE Change means robust SEs are larger (more conservative).")

def run_wls_regression(data: pd.DataFrame, 
                       y_var: str, 
                       x_vars: List[str],
                       weight_var: str) -> sm.regression.linear_model.RegressionResultsWrapper:
    """Run Weighted Least Squares regression."""
    y = data[y_var]
    X = sm.add_constant(data[x_vars])
    
    # Weights = 1/variance (inverse variance weighting)
    weights = 1 / data[weight_var]
    
    model = sm.WLS(y, X, weights=weights).fit()
    return model


# Run WLS
print("="*80)
print("WEIGHTED LEAST SQUARES (Inverse Variance Weighting)")
print("="*80)

model_wls = run_wls_regression(cs_data, 'CAR', 
                                ['eps_surprise', 'log_mcap', 'momentum'],
                                'car_var')
print(model_wls.summary())

# Compare OLS, Robust OLS, and WLS
print("\n" + "="*80)
print("COMPARISON: OLS vs Robust OLS vs WLS")
print("="*80)

# Create comprehensive comparison
vars_to_compare = ['const', 'eps_surprise', 'log_mcap', 'momentum']

comparison_full = pd.DataFrame({
    'Variable': vars_to_compare,
    'OLS Coef': [model2.params[v] for v in vars_to_compare],
    'OLS t': [model2.tvalues[v] for v in vars_to_compare],
    'Robust t': [model2_robust_hc.tvalues[v] for v in range(len(vars_to_compare))],
    'WLS Coef': [model_wls.params[v] for v in vars_to_compare],
    'WLS t': [model_wls.tvalues[v] for v in vars_to_compare],
})

def sig(t, df=20):
    p = 2 * (1 - stats.t.cdf(abs(t), df))
    if p < 0.01: return '***'
    elif p < 0.05: return '**'
    elif p < 0.10: return '*'
    return ''

print("\n" + "-"*80)
print(f"{'Variable':<15} {'OLS Coef':>12} {'OLS t':>10} {'Robust t':>10} {'WLS Coef':>12} {'WLS t':>10}")
print("-"*80)
for _, row in comparison_full.iterrows():
    print(f"{row['Variable']:<15} {row['OLS Coef']:>12.6f} {row['OLS t']:>8.2f}{sig(row['OLS t']):2s} "
          f"{row['Robust t']:>8.2f}{sig(row['Robust t']):2s} {row['WLS Coef']:>12.6f} {row['WLS t']:>8.2f}{sig(row['WLS t']):2s}")
print("-"*80)
print(f"{'R-squared':<15} {model2.rsquared:>12.4f} {'':<22} {model_wls.rsquared:>12.4f}")
print(f"{'N':<15} {int(model2.nobs):>12} {'':<22} {int(model_wls.nobs):>12}")
print("\nSignificance: *** p<0.01, ** p<0.05, * p<0.10")

def calculate_vif(data: pd.DataFrame, variables: List[str]) -> pd.DataFrame:
    """Calculate Variance Inflation Factors."""
    X = data[variables]
    X = sm.add_constant(X)
    
    vif_data = []
    for i, var in enumerate(X.columns):
        if var == 'const':
            continue
        vif = variance_inflation_factor(X.values, i)
        vif_data.append({'Variable': var, 'VIF': vif})
    
    return pd.DataFrame(vif_data)


# Calculate VIFs
print("Multicollinearity Diagnostics:")
print("="*50)

vif_df = calculate_vif(cs_data, ['eps_surprise', 'log_mcap', 'momentum', 'beta', 'sigma'])
vif_df['VIF'] = vif_df['VIF'].apply(lambda x: f"{x:.2f}")
vif_df['Concern'] = vif_df['VIF'].apply(lambda x: 'High' if float(x) > 10 else 'Moderate' if float(x) > 5 else 'Low')

print("\nVariance Inflation Factors:")
print(vif_df.to_string(index=False))
print("\nRule of thumb: VIF > 10 indicates serious multicollinearity")

# Create interaction term
cs_data['surprise_x_size'] = cs_data['eps_surprise'] * cs_data['log_mcap']

# Model with interaction
print("="*80)
print("MODEL WITH INTERACTION: Does Surprise Effect Vary by Size?")
print("="*80)

model_interact = run_cross_sectional_regression(
    cs_data, 'CAR', 
    ['eps_surprise', 'log_mcap', 'surprise_x_size', 'momentum']
)
print(model_interact.summary())

# Visualize interaction effect
fig, ax = plt.subplots(figsize=(10, 6))

# Split by size (median)
median_mcap = cs_data['log_mcap'].median()
small = cs_data[cs_data['log_mcap'] < median_mcap]
large = cs_data[cs_data['log_mcap'] >= median_mcap]

# Scatter plots
ax.scatter(small['eps_surprise'], small['CAR']*100, color='blue', alpha=0.7, 
           s=80, label=f'Small firms (n={len(small)})')
ax.scatter(large['eps_surprise'], large['CAR']*100, color='red', alpha=0.7, 
           s=80, label=f'Large firms (n={len(large)})')

# Separate regression lines
for data, color, label in [(small, 'blue', 'Small'), (large, 'red', 'Large')]:
    if len(data) > 2:
        z = np.polyfit(data['eps_surprise'], data['CAR']*100, 1)
        p = np.poly1d(z)
        x_line = np.linspace(data['eps_surprise'].min(), data['eps_surprise'].max(), 100)
        ax.plot(x_line, p(x_line), f'{color[0]}--', linewidth=2, alpha=0.7)

ax.axhline(0, color='gray', linestyle='-', linewidth=0.5)
ax.axvline(0, color='gray', linestyle='-', linewidth=0.5)
ax.set_xlabel('EPS Surprise', fontsize=12)
ax.set_ylabel('CAR (-1,+1) %', fontsize=12)
ax.set_title('Earnings Response: Small vs Large Firms', fontsize=14)
ax.legend()

plt.tight_layout()
plt.show()

# Report slope differences
if len(small) > 2 and len(large) > 2:
    slope_small = np.polyfit(small['eps_surprise'], small['CAR'], 1)[0]
    slope_large = np.polyfit(large['eps_surprise'], large['CAR'], 1)[0]
    print(f"\nEarnings Response Coefficient:")
    print(f"  Small firms: {slope_small:.4f} ({slope_small*100:.2f}% CAR per unit surprise)")
    print(f"  Large firms: {slope_large:.4f} ({slope_large*100:.2f}% CAR per unit surprise)")

# Subsample: Positive vs Negative Surprises
print("="*80)
print("SUBSAMPLE ANALYSIS: Positive vs Negative Surprises")
print("="*80)

pos_surprise = cs_data[cs_data['eps_surprise'] > 0]
neg_surprise = cs_data[cs_data['eps_surprise'] <= 0]

print(f"\nPositive surprises: n={len(pos_surprise)}, Mean CAR={pos_surprise['CAR'].mean()*100:.2f}%")
print(f"Negative surprises: n={len(neg_surprise)}, Mean CAR={neg_surprise['CAR'].mean()*100:.2f}%")

# Run regressions on subsamples
if len(pos_surprise) > 5:
    model_pos = run_cross_sectional_regression(pos_surprise, 'CAR', ['eps_surprise', 'log_mcap'])
    print(f"\nPositive Surprise Subsample:")
    print(f"  EPS coefficient: {model_pos.params['eps_surprise']:.4f} (t={model_pos.tvalues['eps_surprise']:.2f})")

if len(neg_surprise) > 5:
    model_neg = run_cross_sectional_regression(neg_surprise, 'CAR', ['eps_surprise', 'log_mcap'])
    print(f"\nNegative Surprise Subsample:")
    print(f"  EPS coefficient: {model_neg.params['eps_surprise']:.4f} (t={model_neg.tvalues['eps_surprise']:.2f})")

# Subsample: By Sector
print("\n" + "="*80)
print("SUBSAMPLE ANALYSIS: By Sector")
print("="*80)

sector_results = []
for sector in cs_data['sector'].unique():
    sector_data = cs_data[cs_data['sector'] == sector]
    
    if len(sector_data) >= 4:
        try:
            model = run_cross_sectional_regression(sector_data, 'CAR', ['eps_surprise'])
            sector_results.append({
                'Sector': sector,
                'N': len(sector_data),
                'Mean CAR': sector_data['CAR'].mean(),
                'EPS Coef': model.params['eps_surprise'],
                't-stat': model.tvalues['eps_surprise'],
                'R-squared': model.rsquared
            })
        except:
            pass

sector_df = pd.DataFrame(sector_results)
sector_df['Mean CAR'] = sector_df['Mean CAR'].apply(lambda x: f"{x*100:.2f}%")
sector_df['EPS Coef'] = sector_df['EPS Coef'].apply(lambda x: f"{x:.4f}")
sector_df['t-stat'] = sector_df['t-stat'].apply(lambda x: f"{x:.2f}")
sector_df['R-squared'] = sector_df['R-squared'].apply(lambda x: f"{x:.3f}")

print("\nEarnings Response by Sector:")
print(sector_df.to_string(index=False))

# Create SCAR variable
cs_data['SCAR'] = cs_data['CAR'] / cs_data['car_se']

print("="*80)
print("SCAR REGRESSION (Standardized Dependent Variable)")
print("="*80)

model_scar = run_cross_sectional_regression(cs_data, 'SCAR', 
                                             ['eps_surprise', 'log_mcap', 'momentum'])
print(model_scar.summary())

# Compare CAR vs SCAR regressions
print("\n" + "="*70)
print("COMPARISON: CAR vs SCAR as Dependent Variable")
print("="*70)

print("\nNote: SCAR regression gives t-statistics directly comparable to")
print("standardized test statistics from Sessions 4-5.")
print(f"\nMean SCAR in sample: {cs_data['SCAR'].mean():.3f}")
print(f"Std SCAR: {cs_data['SCAR'].std():.3f}")
print(f"\nIf SCAR ~ N(0,1) under null, mean significantly > 0 indicates positive abnormal returns.")

# Demonstrate using pre-determined variables
print("="*80)
print("ADDRESSING ENDOGENEITY: Using Pre-Determined Variables")
print("="*80)

print("\nVariables in our model and their timing:")
print("-" * 60)
print("  eps_surprise:  Realized at event (potentially endogenous)")
print("  log_mcap:      Measured before event (pre-determined)")
print("  momentum:      60-day pre-event return (pre-determined)")
print("  beta:          Estimated from pre-event data (pre-determined)")
print("  sigma:         Estimated from pre-event data (pre-determined)")

# Model using only pre-determined variables
print("\nModel with only pre-determined variables:")
model_predet = run_cross_sectional_regression(cs_data, 'CAR', 
                                               ['log_mcap', 'momentum', 'beta', 'sigma'])
print(f"R-squared: {model_predet.rsquared:.4f}")
print(f"\nNote: Low R-squared without eps_surprise suggests surprise")
print("      is the main driver, as expected for earnings announcements.")

def create_regression_table(models: List, model_names: List[str], 
                            var_order: List[str]) -> pd.DataFrame:
    """Create publication-style regression table."""
    rows = []
    
    for var in var_order:
        coef_row = {'Variable': var}
        se_row = {'Variable': ''}
        
        for model, name in zip(models, model_names):
            if var in model.params.index:
                coef = model.params[var]
                se = model.bse[var]
                t = model.tvalues[var]
                p = model.pvalues[var]
                
                sig = '***' if p < 0.01 else '**' if p < 0.05 else '*' if p < 0.10 else ''
                coef_row[name] = f"{coef:.4f}{sig}"
                se_row[name] = f"({se:.4f})"
            else:
                coef_row[name] = ''
                se_row[name] = ''
        
        rows.append(coef_row)
        rows.append(se_row)
    
    # Add statistics
    rows.append({'Variable': ''})
    
    r2_row = {'Variable': 'R-squared'}
    n_row = {'Variable': 'N'}
    
    for model, name in zip(models, model_names):
        r2_row[name] = f"{model.rsquared:.3f}"
        n_row[name] = f"{int(model.nobs)}"
    
    rows.append(r2_row)
    rows.append(n_row)
    
    return pd.DataFrame(rows)


# Create table
print("\n" + "="*90)
print("TABLE: Cross-Sectional Determinants of Earnings Announcement Returns")
print("="*90)
print("Dependent Variable: CAR(-1,+1)")
print("-"*90)

models = [model1, model2, model_wls, model_interact]
names = ['(1)', '(2) OLS', '(3) WLS', '(4) Interact']
var_order = ['const', 'eps_surprise', 'log_mcap', 'momentum', 'surprise_x_size']

table = create_regression_table(models, names, var_order)
print(table.to_string(index=False))
print("-"*90)
print("Standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.10")
print("WLS uses inverse CAR variance as weights.")

# Exercise 3: Clustered Standard Errors
print("Exercise 3: Clustered Standard Errors")
print("="*60)

# Cluster by sector
y = cs_data['CAR']
X = sm.add_constant(cs_data[['eps_surprise', 'log_mcap', 'momentum']])

# OLS with clustered SEs by sector
model_cluster = sm.OLS(y, X).fit(cov_type='cluster', 
                                  cov_kwds={'groups': cs_data['sector']})

print("\nOLS with Standard Errors Clustered by Sector:")
print("-"*60)
print(f"{'Variable':<15} {'Coef':>12} {'Clustered SE':>14} {'t-stat':>10}")
print("-"*60)
for var in model_cluster.params.index:
    coef = model_cluster.params[var]
    se = model_cluster.bse[var]
    t = model_cluster.tvalues[var]
    print(f"{var:<15} {coef:>12.6f} {se:>14.6f} {t:>10.2f}")