"""Valuation of European call options
in Black-Scholes-Merton model
incl. vega function and implied volatility estimation
"""

#
# Author: Yehui He <yehui.he@hotmail.com>
# License: Apache-2.0 License

def bsm_call_value(S0, K, T, r, sigma):
    """Valuation of European call option in BSM model.
    Analytical formula.
    Parameters
    ==========
    S0 : float
        initial stock/index level
    K : float
        strike price
    T : float
    maturity date (in year fractions)
    r : float
        constant risk-free short rate
    sigma : float
        volatility factor in diffusion term

    Returns
    =======
    value : float
        present value of the European call option
    """
    from math import log, sqrt, exp
    from scipy import stats

    S0 = float(S0)
    d1 = (log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
    d2 = (log(S0 / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
    # stats.norm.cdf --> cumulative distribution function
    # for normal distribution
    value = (S0 * stats.norm.cdf(d1, 0.0, 1.0) -
             K * exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0))
    return value


def bsm_vega(S0, K, T, r, sigma):
    """Vega of European option in BSM model.

    Parameters
    ----------
    S0 : float
        initial stock/index level

    K : float
        Strick price

    T : float
        maturity date (in year fractions)

    r : float
        constant risk-free short rate

    sigma : float
        volatility factor in diffusion term

    Returns
    -------
    vega : float
        partial derivative of BSM formula with respect
        to sigma, i.e. vega
    """
    from math import log, sqrt
    from scipy import stats

    S0 = float(S0)

    d1 = (log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
    vega = S0 * stats.norm.pdf(d1, 0.0, 1.0) * sqrt(T)
    return vega


def bsm_call_imp_vol(S0, K, T, r, C0, sigma_est, it=100):
    """Implied volatility of European call option in BSM model.

    Parameters
    ==========
    S0: float
        initial stock/index level
    K: float
        strike price
    T: float
        maturity date (in year fractions)
    r: float
        constant risk-free short rate
    sigma_est: float
        estimate of impl. volatility
    it: integer
        number of iterations
    Returns
    =======
    simga_est: float
        numerically estimated implied volatility
    """
    for i in range(it):
        sigma_est -= ((bsm_call_value(S0, K, T, r, sigma_est) - C0) /
                      bsm_vega(S0, K, T, r, sigma_est))
    return sigma_est
