Quantitative Analysis of Tail Risk

import numpy as np

def calculate_bkm_skew(strikes, prices, is_call, S, r, T, F):
    """
    Calculate BKM model-free skewness from option chain.
    
    Parameters:
    -----------
    strikes : array-like
        Strike prices (sorted ascending)
    prices : array-like
        Option prices (OTM only, bid-ask midpoints)
    is_call : array-like
        Boolean array indicating if option is a call
    S : float
        Current spot price
    r : float
        Risk-free rate (annualized)
    T : float
        Time to maturity (years)
    F : float
        Forward price = S * exp(r*T)
    
    Returns:
    --------
    float : CBOE SKEW index value
    """
    # 1. Setup Data
    strikes = np.array(strikes)
    prices = np.array(prices)
    
    # 2. Calculate Strike Intervals (dK) - Trapezoidal Rule
    dK = np.zeros(len(strikes))
    if len(strikes) > 1:
        # Interior points: average of adjacent intervals
        dK[1:-1] = (strikes[2:] - strikes[:-2]) / 2
        # Boundary points
        dK[0] = strikes[1] - strikes[0]
        dK[-1] = strikes[-1] - strikes[-2]
    
    # 3. Initialize Moments
    V = 0  # Variance (2nd moment)
    W = 0  # Cubic (3rd moment - Skewness)
    X = 0  # Quartic (4th moment - Kurtosis)
    
    # 4. Summation Loop (Numerical Integration)
    for i, K in enumerate(strikes):
        # Filter for OTM only (puts below forward, calls above)
        if (K < F and is_call[i]) or (K > F and not is_call[i]):
            continue
        
        # Log-moneyness: k = ln(K/F)
        k_val = np.log(K / F)
        
        # BKM Weighting Functions (from Bakshi-Kapadia-Madan 1999)
        # Variance weight
        v_weight = (2 * (1 - k_val)) / (K**2)
        # Skewness weight (cubic)
        w_weight = (6 * k_val - 3 * (k_val**2)) / (K**2)
        # Kurtosis weight (quartic)
        x_weight = (12 * (k_val**2) - 4 * (k_val**3)) / (K**2)
        
        # Add weighted contributions
        V += v_weight * prices[i] * dK[i]
        W += w_weight * prices[i] * dK[i]
        X += x_weight * prices[i] * dK[i]
    
    # 5. Discount to Present Value
    df = np.exp(r * T)
    V *= df
    W *= df
    X *= df
    
    # 6. Mean Adjustment Correction
    # Account for drift in risk-neutral measure
    mu = np.exp(r * T) - 1 - V/2 - W/6 - X/24
    
    # 7. Calculate Raw Skewness (3rd standardized moment)
    # Skew = E[(X-μ)³] / σ³
    variance = V - mu**2
    if variance <= 0:
        return np.nan  # Invalid: negative variance
    
    skew_raw = (W - 3*mu*V + 2*mu**3) / (variance**1.5)
    
    # 8. Transform to CBOE Index Scale
    # SKEW = 100 - 10*S (so negative skew → values > 100)
    skew_index = 100 - 10 * skew_raw
    
    return skew_index


# Example Usage
def example_calculation():
    """
    Example: Calculate SKEW from sample option chain
    """
    # Sample data (simplified)
    strikes = np.array([4800, 4850, 4900, 4950, 5000, 5050, 5100, 5150, 5200])
    prices = np.array([5.2, 8.5, 13.2, 20.1, 30.5, 22.3, 15.8, 10.2, 6.5])
    is_call = np.array([False, False, False, False, False, True, True, True, True])
    
    S = 5000  # Spot price
    r = 0.05  # 5% risk-free rate
    T = 30/365  # 30 days to expiration
    F = S * np.exp(r * T)  # Forward price
    
    skew = calculate_bkm_skew(strikes, prices, is_call, S, r, T, F)
    print(f"CBOE SKEW Index: {skew:.2f}")
    
    # Interpretation
    if skew > 135:
        print("⚠️  Extreme tail risk pricing")
    elif skew > 120:
        print("📊 Elevated skew - market hedging")
    else:
        print("✅ Normal distribution approximation")

example_calculation()