pricing.py

import numpy as np
from scipy.stats import norm
from scipy.optimize import brentq

def black_scholes(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    if option_type == "call":
        return S*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)
    else:
        return K*np.exp(-r*T)*norm.cdf(-d2) - S*norm.cdf(-d1)

def binomial_tree(S, K, T, r, sigma, N=100, option_type="call", american=False):
    dt = T/N
    u = np.exp(sigma*np.sqrt(dt))
    d = 1/u
    p = (np.exp(r*dt) - d) / (u - d)
    discount = np.exp(-r*dt)
    ST = np.array([S * u**j * d**(N-j) for j in range(N+1)])
    if option_type == "call":
        V = np.maximum(ST - K, 0)
    else:
        V = np.maximum(K - ST, 0)
    for i in range(N-1, -1, -1):
        ST = np.array([S * u**j * d**(i-j) for j in range(i+1)])
        V = discount * (p*V[1:] + (1-p)*V[:-1])
        if american:
            if option_type == "call":
                V = np.maximum(V, ST - K)
            else:
                V = np.maximum(V, K - ST)
    return V[0]

def monte_carlo(S, K, T, r, sigma, n_simulations=100000, option_type="call"):
    np.random.seed(42)
    Z = np.random.standard_normal(n_simulations)
    ST = S * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
    if option_type == "call":
        payoff = np.maximum(ST - K, 0)
    else:
        payoff = np.maximum(K - ST, 0)
    price = np.exp(-r*T) * np.mean(payoff)
    std_error = np.exp(-r*T) * np.std(payoff) / np.sqrt(n_simulations)
    return price, std_error

def asian_option(S, K, T, r, sigma, n_simulations=50000, n_steps=252, option_type="call"):
    np.random.seed(42)
    dt = T/n_steps
    paths = np.zeros((n_simulations, n_steps+1))
    paths[:, 0] = S
    for t in range(1, n_steps+1):
        Z = np.random.standard_normal(n_simulations)
        paths[:, t] = paths[:, t-1] * np.exp((r - 0.5*sigma**2)*dt + sigma*np.sqrt(dt)*Z)
    avg_price = paths[:, 1:].mean(axis=1)
    if option_type == "call":
        payoff = np.maximum(avg_price - K, 0)
    else:
        payoff = np.maximum(K - avg_price, 0)
    return np.exp(-r*T) * np.mean(payoff)

def implied_volatility(market_price, S, K, T, r, option_type="call"):
    try:
        iv = brentq(lambda sigma: black_scholes(S, K, T, r, sigma, option_type) - market_price,
                    0.001, 10.0, xtol=1e-6)
        return iv
    except:
        return np.nan