# VaR Monte Carlo - Lo que usaras en el trabajo
# ✅ Este codigo se puede ejecutar directamente en Google Colab
import numpy as np
from scipy.stats import norm

def var_monte_carlo(position, mu, sigma, days=1, n_sims=100000, confidence=0.95):
    """
    Calcula VaR usando Monte Carlo.
    
    Args:
        position: Valor de la posicion en USD
        mu: Retorno esperado diario (tipicamente ~0 para VaR)
        sigma: Volatilidad diaria
        days: Horizonte (1 para diario, 10 para Basilea)
        n_sims: Numero de simulaciones (100k es standard)
        confidence: 0.95 o 0.99
    
    Returns:
        VaR en USD
    """
    # Simular retornos
    np.random.seed(42)  # Para reproducibilidad
    returns = np.random.normal(mu * days, sigma * np.sqrt(days), n_sims)
    
    # Calcular P&L
    pnl = position * returns
    
    # VaR es el percentil negativo
    var = -np.percentile(pnl, (1 - confidence) * 100)
    
    return var

# Ejemplo: $50M en AAPL, vol diaria 2.1%
var_95 = var_monte_carlo(
    position=50_000_000,
    mu=0,           # Asumimos drift 0 para VaR
    sigma=0.021,    # 2.1% vol diaria
    days=1,
    confidence=0.95
)

print(f"VaR 95% (1 dia): ${var_95:,.0f}")
# Output: VaR 95% (1 dia): $1,726,842

# Expected Shortfall - Codigo de produccion
# ✅ Este codigo se puede ejecutar directamente en Google Colab
import numpy as np
from scipy.stats import norm

def expected_shortfall(position, sigma, confidence=0.975):
    """
    ES parametrico (asumiendo normalidad).
    Basilea IV usa 97.5% de confianza.
    """
    z = norm.ppf(confidence)
    phi_z = norm.pdf(z)  # Densidad en el cuantil
    es = position * sigma * phi_z / (1 - confidence)
    return es

def es_historico(returns, position, confidence=0.975):
    """
    ES historico - promedio de perdidas peores que VaR.
    """
    pnl = position * returns
    var = -np.percentile(pnl, (1 - confidence) * 100)
    # Promedio de perdidas que exceden el VaR
    tail_losses = pnl[pnl < -var]
    es = -np.mean(tail_losses)
    return es

# Ejemplo
position = 100_000_000
sigma_daily = 0.015

es_param = expected_shortfall(position, sigma_daily)
print(f"ES 97.5% (parametrico): ${es_param:,.0f}")
# Output: ES 97.5% (parametrico): $3,709,823

# Black-Scholes Greeks - Codigo de produccion
# ✅ Este codigo se puede ejecutar directamente en Google Colab
import numpy as np
from scipy.stats import norm

def black_scholes_greeks(S, K, T, r, sigma, option_type='call'):
    """
    Calcula precio y Greeks de una opcion europea.
    
    Args:
        S: Precio spot del subyacente
        K: Strike price
        T: Tiempo al vencimiento (en años)
        r: Tasa libre de riesgo (anual)
        sigma: Volatilidad implicita (anual)
        option_type: 'call' o 'put'
    
    Returns:
        dict con precio y todas las Greeks
    """
    # d1 y d2
    d1 = (np.log(S/K) + (r + sigma**2/2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    
    # Precio
    if option_type == 'call':
        price = S*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)
        delta = norm.cdf(d1)
    else:  # put
        price = K*np.exp(-r*T)*norm.cdf(-d2) - S*norm.cdf(-d1)
        delta = norm.cdf(d1) - 1
    
    # Greeks (iguales para call y put excepto Delta y Rho)
    gamma = norm.pdf(d1) / (S * sigma * np.sqrt(T))
    vega = S * norm.pdf(d1) * np.sqrt(T) / 100  # Por 1% cambio en vol
    theta = -(S*norm.pdf(d1)*sigma)/(2*np.sqrt(T)) - r*K*np.exp(-r*T)*norm.cdf(d2)
    theta = theta / 365  # Por dia
    
    return {
        'price': price,
        'delta': delta,
        'gamma': gamma,
        'vega': vega,
        'theta': theta
    }

# Ejemplo: Tesla call option
greeks = black_scholes_greeks(
    S=248,      # TSLA spot price
    K=250,      # Strike
    T=30/365,   # 30 dias
    r=0.05,     # 5% risk-free rate
    sigma=0.45  # 45% implied vol (TSLA es volatil)
)

print("TSLA $250 Call (30 dias):")
print(f"  Precio:  ${greeks['price']:.2f}")
print(f"  Delta:   {greeks['delta']:.4f}")
print(f"  Gamma:   {greeks['gamma']:.4f}")
print(f"  Vega:    ${greeks['vega']:.4f} por 1% vol")
print(f"  Theta:   ${greeks['theta']:.4f} por dia")

# Output:
# TSLA $250 Call (30 dias):
#   Precio:  $12.47
#   Delta:   0.4892
#   Gamma:   0.0142
#   Vega:    $0.2834 por 1% vol
#   Theta:   $-0.1523 por dia

# EWMA Volatility - Metodo RiskMetrics
# ✅ Este codigo se puede ejecutar directamente en Google Colab
import numpy as np
import pandas as pd

def ewma_volatility(returns, lambda_param=0.94):
    """
    Calcula volatilidad EWMA (metodo RiskMetrics).
    
    Args:
        returns: Serie de retornos
        lambda_param: Factor de decaimiento (0.94 es standard)
    
    Returns:
        Serie de volatilidades EWMA
    """
    returns = np.array(returns)
    n = len(returns)
    
    # Inicializar con varianza muestral
    var = np.zeros(n)
    var[0] = returns[0]**2
    
    # Recursion EWMA
    for t in range(1, n):
        var[t] = lambda_param * var[t-1] + (1 - lambda_param) * returns[t-1]**2
    
    return np.sqrt(var)

# Ejemplo con datos reales de SPY
np.random.seed(42)
# Simular retornos con cluster de volatilidad (como en la vida real)
returns = np.concatenate([
    np.random.normal(0, 0.01, 50),    # Periodo tranquilo
    np.random.normal(0, 0.03, 30),    # Spike de volatilidad
    np.random.normal(0, 0.015, 40)   # Vuelve a la normalidad
])

vol_ewma = ewma_volatility(returns)
vol_simple = pd.Series(returns).rolling(20).std()

print(f"EWMA vol al final: {vol_ewma[-1]*100:.2f}% diaria")
print(f"Simple vol (20d):  {vol_simple.iloc[-1]*100:.2f}% diaria")

# EWMA reacciona mas rapido a cambios en volatilidad

# GARCH usando la libreria arch (standard en la industria)
# ✅ Este codigo se puede ejecutar en Google Colab
# Ejecuta primero: !pip install arch

from arch import arch_model
import pandas as pd
import numpy as np

# Simular retornos (en produccion usas datos reales de Bloomberg)
np.random.seed(42)
returns = np.random.normal(0, 0.01, 1000) * 100  # En porcentaje

# Fit GARCH(1,1)
model = arch_model(returns, vol='Garch', p=1, q=1)
result = model.fit(disp='off')

print(result.summary())

# Pronosticar volatilidad para los proximos 5 dias
forecast = result.forecast(horizon=5)
vol_forecast = np.sqrt(forecast.variance.iloc[-1])

print("\nPronostico de volatilidad (proximos 5 dias):")
for i, v in enumerate(vol_forecast, 1):
    print(f"  Dia {i}: {v:.2f}% diaria ({v*np.sqrt(252):.1f}% anual)")