def delta(flag, S, K, t, r, sigma):
    """Return Black-Scholes delta of an option.
    
    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str      

    Example 17.1, page 355, Hull:
    
    >>> S = 49
    >>> K = 50 
    >>> r = .05
    >>> t = 0.3846
    >>> sigma = 0.2
    >>> flag = 'c'
    >>> delta_calc = delta(flag, S, K, t, r, sigma)
    >>> # 0.521601633972
    >>> delta_text_book = 0.522
    >>> abs(delta_calc - delta_text_book) < .01
    True
    """

    d_1 = d1(S, K, t, r, sigma)

    if flag == 'p':
        return N(d_1) - 1.0
    else:
        return N(d_1) def gamma(flag, S, K, t, r, sigma):
    """Return Black-Scholes gamma of an option.
    
    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str      
    
    Example 17.4, page 364, Hull:

    >>> S = 49
    >>> K = 50 
    >>> r = .05
    >>> t = 0.3846
    >>> sigma = 0.2
    >>> flag = 'c'
    >>> gamma_calc = gamma(flag, S, K, t, r, sigma)
    >>> # 0.0655453772525
    >>> gamma_text_book = 0.066
    >>> abs(gamma_calc - gamma_text_book) < .001
    True
    """

    d_1 = d1(S, K, t, r, sigma)
    v_squared = sigma**2
    return pdf(d_1)/(S*sigma*numpy.sqrt(t)) def theta(flag, S, K, t, r, sigma):
    """Return Black-Scholes theta of an option.
    
    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str      

    The text book analytical formula does not divide by 365,
    but in practice theta is defined as the change in price
    for each day change in t, hence we divide by 365.

    Example 17.2, page 359, Hull:

    >>> S = 49
    >>> K = 50 
    >>> r = .05
    >>> t = 0.3846
    >>> sigma = 0.2
    >>> flag = 'c'
    >>> annual_theta_calc = theta(flag, S, K, t, r, sigma) * 365
    >>> # -4.30538996455
    >>> annual_theta_text_book = -4.31
    >>> abs(annual_theta_calc - annual_theta_text_book) < .01
    True

    Using the same inputs with a put.
    >>> S = 49
    >>> K = 50 
    >>> r = .05
    >>> t = 0.3846
    >>> sigma = 0.2
    >>> flag = 'p'
    >>> annual_theta_calc = theta(flag, S, K, t, r, sigma) * 365
    >>> # -1.8530056722
    >>> annual_theta_reference = -1.8530056722
    >>> abs(annual_theta_calc - annual_theta_reference) < .000001
    True
    """
    
    two_sqrt_t = 2 * numpy.sqrt(t)

    D1 = d1(S, K, t, r, sigma)
    D2 = d2(S, K, t, r, sigma)

    first_term = (-S * pdf(D1) * sigma) / two_sqrt_t 

    if flag == 'c':

        second_term = r * K * numpy.exp(-r*t) * N(D2)
        return (first_term - second_term)/365.0
    
    if flag == 'p':
    
        second_term = r * K * numpy.exp(-r*t) * N(-D2)
        return (first_term + second_term)/365.0 def vega(flag, S, K, t, r, sigma):
    """Return Black-Scholes vega of an option.
    
    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str
    
    The text book analytical formula does not multiply by .01,
    but in practice vega is defined as the change in price
    for each 1 percent change in IV, hence we multiply by 0.01.

    Example 17.6, page 367, Hull:
    
    >>> S = 49
    >>> K = 50 
    >>> r = .05
    >>> t = 0.3846
    >>> sigma = 0.2
    >>> flag = 'c'
    >>> vega_calc = vega(flag, S, K, t, r, sigma)
    >>> # 0.121052427542
    >>> vega_text_book = 0.121
    >>> abs(vega_calc - vega_text_book) < .01
    True
    """

    d_1 = d1(S, K, t, r, sigma)
    return S * pdf(d_1) * numpy.sqrt(t) * 0.01 def rho(flag, S, K, t, r, sigma):
    """Return Black-Scholes rho of an option.
    
    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str      
    
    The text book analytical formula does not multiply by .01,
    but in practice rho is defined as the change in price
    for each 1 percent change in r, hence we multiply by 0.01.

    Example 17.7, page 368, Hull:

    >>> S = 49
    >>> K = 50 
    >>> r = .05
    >>> t = 0.3846
    >>> sigma = 0.2
    >>> flag = 'c'
    >>> rho_calc = rho(flag, S, K, t, r, sigma)
    >>> # 0.089065740988
    >>> rho_text_book = 0.0891
    >>> abs(rho_calc - rho_text_book) < .0001
    True
    """

    d_2 = d2(S, K, t, r, sigma)
    e_to_the_minus_rt = numpy.exp(-r*t)
    if flag == 'c':
        return t*K*e_to_the_minus_rt * N(d_2) * .01
    else:
        return -t*K*e_to_the_minus_rt * N(-d_2) * .01 def delta(flag, S, K, t, r, sigma):
    """Return Black-Scholes delta of an option.
    
    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str
    """
    
    b = r

    return numerical_delta(flag, S, K, t, r, sigma, b, f) def gamma(flag, S, K, t, r, sigma):
    """Return Black-Scholes gamma of an option.

    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str
    """

    b = r

    return numerical_gamma(flag, S, K, t, r, sigma, b, f) def theta(flag, S, K, t, r, sigma):
    """Return Black-Scholes theta of an option.

    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str      

    """    
    
    b = r

    return numerical_theta(flag, S, K, t, r, sigma, b, f) def vega(flag, S, K, t, r, sigma):
    """Return Black-Scholes vega of an option.

    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str
    """

    b = r

    return numerical_vega(flag, S, K, t, r, sigma, b, f) def rho(flag, S, K, t, r, sigma):
    """Return Black-Scholes rho of an option.

    :param S: underlying asset price
    :type S: float
    :param K: strike price
    :type K: float
    :param sigma: annualized standard deviation, or volatility
    :type sigma: float
    :param t: time to expiration in years
    :type t: float
    :param r: risk-free interest rate
    :type r: float
    :param flag: 'c' or 'p' for call or put.
    :type flag: str
    """    
    
    b = r

    return numerical_rho(flag, S, K, t, r, sigma, b, f)
from itertools import product
from py_vollib.black_scholes.greeks.analytical import delta
from py_vollib.black_scholes.greeks.analytical import gamma
from py_vollib.black_scholes.greeks.analytical import theta
from py_vollib.black_scholes.greeks.analytical import vega
from py_vollib.black_scholes.greeks.analytical import rho
from py_vollib.black_scholes.greeks.numerical import delta as ndelta
from py_vollib.black_scholes.greeks.numerical import gamma as ngamma
from py_vollib.black_scholes.greeks.numerical import theta as ntheta
from py_vollib.black_scholes.greeks.numerical import vega as nvega
from py_vollib.black_scholes.greeks.numerical import rho as nrho
import unittest
import numpy
# -----------------------------------------------------------------------------
# IMPORTS

# Standard library imports
from __future__ import print_function

# Related third party imports

# Local application/library specific imports


class TestBSGreeks(unittest.TestCase):
    def setUp(self):
        self.epsilon = 0.001
        self.flags = ['c', 'p']
        self.S = 100
        self.Ks = numpy.linspace(20, 200, 10)
        self.ts = numpy.linspace(0.01, 2, 10)
        self.rs = numpy.linspace(0, 0.2, 10)
        self.sigmas = numpy.linspace(0.1, 0.5, 10)
        self.arg_combinations = list(product(self.flags, [self.S], self.Ks, self.ts, self.rs, self.sigmas))

    def diff_mean(self, left, right):
        left_arr = numpy.array(left)
        right_arr = numpy.array(right)
        abs_diff = numpy.abs(left_arr - right_arr)
        return numpy.mean(abs_diff)

    def test_theta(self):
        thetas = [theta(flag, S, K, t, r, sigma) for flag, S, K, t, r, sigma in self.arg_combinations]
        nthetas = [ntheta(flag, S, K, t, r, sigma) for flag, S, K, t, r, sigma in self.arg_combinations]
        self.assertTrue(self.diff_mean(thetas, nthetas) < self.epsilon)

    def test_delta(self):
        deltas = [delta(flag, S, K, t, r, sigma) for flag, S, K, t, r, sigma in self.arg_combinations]
        ndeltas = [ndelta(flag, S, K, t, r, sigma) for flag, S, K, t, r, sigma in self.arg_combinations]
        self.assertTrue(self.diff_mean(deltas, ndeltas) < self.epsilon)

    def test_gamma(self):
        gammas = [gamma(flag, S, K, t, r, sigma) for flag, S, K, t, r, sigma in self.arg_combinations]
        ngammas = [ngamma(flag, S, K, t, r, sigma) for flag, S, K, t, r, sigma in self.arg_combinations]
        self.assertTrue(self.diff_mean(gammas, ngammas) < self.epsilon)

    def test_vega(self):

        vegas = [vega(flag, S, K, t, r, sigma) for flag, S, K, t, r, sigma in self.arg_combinations]