Untitled

import numpy as np
from datetime import datetime

def generate_ig(mu, lambda_, size=None):
    """Generate inverse Gaussian samples with mean mu and shape lambda_."""
    if size is None:
        size = 1
    if isinstance(size, int):
        size = (size,)
    Z = np.random.normal(0, 1, size=size)
    term1 = (mu**2 * Z**2) / (2 * lambda_)
    term2 = (mu / (2 * lambda_)) * np.sqrt(4 * mu * lambda_ * Z**2 + mu**2 * Z**4)
    V = mu + term1 - term2
    U = np.random.uniform(0, 1, size=size)
    mask = U <= mu / (mu + V)
    return np.where(mask, V, mu**2 / V)

def simulate_gh_bridge(start_date: str, start_price: float, 
                       end_date: str, end_price: float, 
                       current_date: str, n_paths: int, 
                       expected_return: float, volatility: float, 
                       alpha: float = 1.0) -> np.ndarray:
    """
    Simulate bridge paths with generalized hyperbolic increments.
    
    Parameters:
    - start_date (str): 'YYYY-MM-DD'
    - start_price (float): Starting price
    - end_date (str): 'YYYY-MM-DD'
    - end_price (float): Ending price
    - current_date (str): 'YYYY-MM-DD', between start and end
    - n_paths (int): Number of paths
    - expected_return (float): Annual expected return
    - volatility (float): Annual volatility
    - alpha (float): GH tail parameter (default 1.0)
    
    Returns:
    - np.ndarray: Shape (n_steps, n_paths) with simulated prices
    """
    # Date handling
    start_dt = datetime.strptime(start_date, '%Y-%m-%d')
    end_dt = datetime.strptime(end_date, '%Y-%m-%d')
    current_dt = datetime.strptime(current_date, '%Y-%m-%d')
    if not (start_dt <= current_dt <= end_dt):
        raise ValueError("Current date must be between start and end dates")
    
    # Time setup
    total_time = (end_dt - start_dt).days / 365.25
    n_steps = (end_dt - start_dt).days + 1
    dt = total_time / (n_steps - 1) if n_steps > 1 else 0
    
    # Parameters
    drift = (expected_return - 0.5 * volatility**2) * dt
    delta = alpha * volatility**2 * dt
    mu_ig = delta / alpha
    lambda_ig = delta * alpha  # Controls kurtosis
    
    # Generate increments
    V = generate_ig(mu_ig, lambda_ig, size=(n_steps-1, n_paths))
    Z = np.random.normal(0, 1, (n_steps-1, n_paths))
    increments = np.sqrt(V) * Z
    
    # Cumulative process
    M = np.vstack([np.zeros((1, n_paths)), np.cumsum(increments, axis=0)])
    t = np.linspace(0, total_time, n_steps)
    bridge_factor = t / total_time if total_time > 0 else np.zeros(n_steps)
    
    # Log prices
    log_S = (np.log(start_price) + drift * t[:, None] + M - 
             bridge_factor[:, None] * M[-1, :] + 
             bridge_factor[:, None] * (np.log(end_price / start_price) - drift * total_time))
    
    # Prices
    paths = np.exp(log_S)
    paths[0, :] = start_price  # Ensure exact start
    paths[-1, :] = end_price   # Ensure exact end
    
    return paths

# Example
if __name__ == "__main__":
    paths = simulate_gh_bridge(
        start_date="2025-01-01", start_price=100.0,
        end_date="2025-03-04", end_price=110.0,
        current_date="2025-02-01", n_paths=5,
        expected_return=0.05, volatility=0.2, alpha=1.0
    )
    print(f"Paths shape: {paths.shape}")
    print(f"Sample path:\n{paths[:, 0]}")