from __future__ import division, print_function, absolute_import
import math
import numpy as np
import scipy.special
import mafipy.function
# ----------------------------------------------------------------------------
# Black scholes european call/put
# ----------------------------------------------------------------------------
def _is_d1_or_d2_infinity(underlying, strike, vol):
"""is_d1_or_d2_infinity
:param float underlying:
:param float strike:
:param float vol:
:return: check whether :math:`d_{1}` and :math:`d_{2}` is infinity or not.
:rtype: bool
"""
return (np.isclose(underlying, 0.0)
or strike < 0.0
or vol < 0.0)
[docs]def func_d1(underlying, strike, rate, maturity, vol):
"""func_d1
calculate :math:`d_{1}` in black scholes formula.
See :py:func:`black_scholes_call_formula`.
:param float underlying: underlying/strike must be non-negative.
:param float strike: underlying/strike must be non-negative.
:param float rate:
:param float maturity: must be non-negative.
:param float vol: must be non-negative.
:return: :math:`d_{1}`.
:rtype: float
"""
assert(underlying / strike >= 0.0)
assert(maturity >= 0.0)
assert(vol >= 0.0)
numerator = (
math.log(underlying / strike) + (rate + vol * vol * 0.5) * maturity)
denominator = vol * math.sqrt(maturity)
return numerator / denominator
[docs]def func_d2(underlying, strike, rate, maturity, vol):
"""func_d2
calculate :math:`d_{2}` in black scholes formula.
See :py:func:`black_scholes_call_formula`.
:param float underlying: underlying/strike must be non-negative.
:param float strike: underlying/strike must be non-negative.
:param float rate:
:param float maturity: must be non-negative.
:param float vol: must be non-negative.
:return: :math:`d_{2}`.
:rtype: float.
"""
assert(underlying / strike >= 0.0)
assert(maturity >= 0.0)
assert(vol >= 0.0)
numerator = (
math.log(underlying / strike) + (rate - vol * vol * 0.5) * maturity)
denominator = vol * math.sqrt(maturity)
return numerator / denominator
[docs]def d_fprime_by_strike(underlying, strike, rate, maturity, vol):
"""d_fprime_by_strike
derivative of :math:`d_{1}` with respect to :math:`K`
where :math:`K` is strike.
See :py:func:`func_d1`.
.. math::
\\frac{\partial }{\partial K} d_{1}(K)
= \\frac{K}{\sigma S \sqrt{T}}.
Obviously, derivative of :math:`d_{1}` and :math:`d_{2}` is same.
That is
.. math::
\\frac{\partial }{\partial K} d_{1}(K)
= \\frac{\partial }{\partial K} d_{2}(K).
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: must be non-negative.
:param float vol:
:return: value of derivative.
:rtype: float
"""
assert(maturity > 0.0)
return - 1.0 / (math.sqrt(maturity) * vol * strike)
[docs]def d_fhess_by_strike(underlying, strike, rate, maturity, vol):
"""d_fhess_by_strike
second derivative of :math:`d_{i}\ (i = 1, 2)` with respect to :math:`K`,
where :math:`K` is strike.
.. math::
\\frac{\partial^{2}}{\partial K^{2}} d_{1}(K)
= \\frac{1}{S \sigma \sqrt{T} },
where
:math:`S` is underlying,
:math:`\sigma` is vol,
:math:`T` is maturity.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity:
:param float vol:
:return: value of second derivative of :math:`d_{1}` or :math:`d_{2}`.
:rtype: float
"""
assert(maturity > 0.0)
return 1.0 / (math.sqrt(maturity) * vol * strike * strike)
def black_scholes_call_formula(underlying, strike, rate, maturity, vol):
"""black_scholes_call_formula
calculate well known black scholes formula for call option.
.. math::
c(S, K, r, T, \sigma)
:= S N(d_{1}) - K e^{-rT} N(d_{2}),
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is vol,
:math:`N(\cdot)` is standard normal distribution,
and :math:`d_{1}` and :math:`d_{2}` are defined as follows:
.. math::
\\begin{eqnarray}
d_{1}
& = &
\\frac{\ln(S/K) + (r + \sigma^{2}/2)T}{\sigma \sqrt{T}},
\\
d_{2}
& = &
\\frac{\ln(S/K) + (r - \sigma^{2}/2)T} {\sigma \sqrt{T}},
\end{eqnarray}
:param float underlying: value of underlying.
:param float strike: strike of call option.
:param float rate: risk free rate.
:param float maturity: year fraction to maturity.
:param float vol: volatility.
:return: call value.
:rtype: float
"""
d1 = func_d1(underlying, strike, rate, maturity, vol)
d2 = func_d2(underlying, strike, rate, maturity, vol)
return (underlying * scipy.special.ndtr(d1)
- strike * math.exp(-rate * maturity) * scipy.special.ndtr(d2))
def black_scholes_put_formula(underlying, strike, rate, maturity, vol):
"""black_scholes_put_formula
calculate well known black scholes formula for put option.
Here value of put option is calculated by put-call parity.
.. math::
\\begin{array}{cccl}
& e^{-rT}(S - K)
& = & c(S, K, r, T, \sigma) - p(S, K, r, T, \sigma)
\\\\
\iff & p(S, K, r, T, \sigma)
& = & c(S, K, r, T, \sigma) - e^{-rT}(S - K)
\end{array}
where
:math:`c(\cdot)` denotes value of call option,
:math:`p(\cdot)` denotes value of put option,
:math:`S` is value of underlying at today,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is vol.
:math:`c(\cdot)` is calculated
by :py:func:`black_scholes_call_formula`.
:param float underlying: value of underlying.
:param float strike: strike of put option.
:param float rate: risk free rate.
:param float maturity: year fraction to maturity.
:param float vol: volatility.
:return: put value.
:rtype: float
"""
call_value = black_scholes_call_formula(
underlying, strike, rate, maturity, vol)
discount = math.exp(-rate * maturity)
return call_value - (underlying - strike * discount)
[docs]def black_scholes_call_value(
underlying,
strike,
rate,
maturity,
vol,
today=0.0):
"""black_scholes_call_value
calculate call value in the case of today is not zero.
(`maturity` - `today`) is treated as time to expiry.
See :py:func:`black_scholes_call_formula`.
* case :math:`S > 0, K < 0`
* return :math:`S - e^{-rT} K`
* case :math:`S < 0, K > 0`
* return 0
* case :math:`S < 0, K < 0`
* return :math:`S - e^{-rT}K + E[(-(S - K))^{+}]`
* case :math:`T \le 0`
* return 0
:param float underlying:
:param float strike:
:param float rate:
:param float maturity:
:param float vol: volatility. This must be positive.
:param float today:
:return: call value.
:rtype: float
"""
assert(vol >= 0.0)
time = maturity - today
# option is expired
if time < 0.0 or np.isclose(time, 0.0):
return 0.0
elif np.isclose(underlying, 0.0):
return math.exp(-rate * time) * max(-strike, 0.0)
elif np.isclose(strike, 0.0) and underlying > 0.0:
return math.exp(-rate * today) * underlying
elif np.isclose(strike, 0.0) and underlying < 0.0:
return 0.0
# never below strike
elif strike < 0.0 and underlying > 0.0:
return underlying - math.exp(-rate * time) * strike
# never beyond strike
elif strike > 0.0 and underlying < 0.0:
return 0.0
elif underlying < 0.0:
# max(S - K, 0) = (S - K) + max(-(S - K), 0)
value = black_scholes_call_formula(
-underlying, -strike, rate, time, vol)
return (underlying - strike) + value
return black_scholes_call_formula(
underlying, strike, rate, time, vol)
[docs]def black_scholes_put_value(
underlying,
strike,
rate,
maturity,
vol,
today=0.0):
"""black_scholes_put_value
evaluates value of put option using put-call parity so that
this function calls :py:func:`black_scholes_call_value`.
See :py:func:`black_scholes_put_formula`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity:
:param float vol:
:param float today:
:return: put value.
:rtype: float
"""
time = maturity - today
# option is expired
if time < 0.0 or np.isclose(time, 0.0):
return 0.0
elif np.isclose(strike, 0.0) and underlying > 0.0:
return 0.0
elif np.isclose(strike, 0.0) and underlying < 0.0:
return underlying * math.exp(-rate * today)
call_value = black_scholes_call_value(
underlying, strike, rate, maturity, vol, today)
discount = math.exp(-rate * time)
return call_value - (underlying - strike * discount)
[docs]def black_scholes_call_value_fprime_by_strike(
underlying, strike, rate, maturity, vol):
"""black_scholes_call_value_fprime_by_strike
First derivative of value of call option with respect to strike
under black scholes model.
See :py:func:`black_scholes_call_formula`.
.. math::
\\frac{\partial }{\partial K} c(K; S, r, T, \sigma)
= - e^{-rT} \Phi(d_{1}(K))
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is vol,
:math:`d_{1}, d_{2}` is defined
in :py:func:`black_scholes_call_formula`,
:math:`\Phi(\cdot)` is c.d.f. of standard normal distribution,
:math:`\phi(\cdot)` is p.d.f. of standard normal distribution.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: must be non-negative.
:param float vol: volatility. must be non-negative.
:return: value of derivative.
:rtype: float
"""
norm = scipy.stats.norm
assert(maturity > 0.0)
d2 = func_d2(underlying, strike, rate, maturity, vol)
discount = math.exp(-rate * maturity)
return -discount * norm.cdf(d2)
[docs]def black_scholes_call_value_fhess_by_strike(
underlying, strike, rate, maturity, vol):
"""black_scholes_call_value_fhess_by_strike
Second derivative of value of call option with respect to strike
under black scholes model.
See :py:func:`black_scholes_call_formula`
and :py:func:`black_scholes_call_value_fprime_by_strike`.
.. math::
\\begin{array}{ccl}
\\frac{\partial^{2}}{\partial K^{2}} c(0, S; T, K)
& = &
-e^{-rT}
\phi(d_{2}(K)) d^{\prime}(K)
\end{array}
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is vol,
:math:`d_{1}, d_{2}` is defined
in :py:func:`black_scholes_call_formula`,
:math:`\Phi(\cdot)` is c.d.f. of standard normal distribution,
:math:`\phi(\cdot)` is p.d.f. of standard normal distribution.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: non-negative.
:param float vol: volatility. non-negative.
:return: value of second derivative.
:rtype: float.
"""
norm = scipy.stats.norm
# option is expired
if maturity < 0.0 or np.isclose(maturity, 0.0):
return 0.0
# never below strike
elif strike <= 0.0 and underlying > 0.0:
return 0.0
# never beyond strike
elif strike > 0.0 and underlying < 0.0:
return 0.0
elif underlying < 0.0 and strike < 0.0:
underlying = -underlying
strike = -strike
discount = math.exp(-rate * maturity)
d2 = func_d2(underlying, strike, rate, maturity, vol)
d_fprime = d_fprime_by_strike(underlying, strike, rate, maturity, vol)
d2_density = norm.pdf(d2)
return -discount * d2_density * d_fprime
[docs]def black_scholes_call_value_third_by_strike(
underlying, strike, rate, maturity, vol):
"""black_scholes_call_value_third_by_strike
Third derivative of value of call option with respect to strike
under black scholes model.
See :py:func:`black_scholes_call_formula`
and :py:func:`black_scholes_call_value_fprime_by_strike`,
and :py:func:`black_scholes_call_value_fhess_by_strike`.
.. math::
\\begin{array}{ccl}
\\frac{\partial^{3}}{\partial K^{3}} c(0, S; T, K)
& = &
-e^{-rT}
\left(
\phi^{\prime}(d_{2}(K))(d^{\prime}(K))^{2}
+ \phi(d_{2}(K))d^{\prime\prime}(K)
\\right)
\end{array}
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is vol,
:math:`d_{1}, d_{2}` is defined
in :py:func:`black_scholes_call_formula`,
:math:`\Phi(\cdot)` is c.d.f. of standard normal distribution,
:math:`\phi(\cdot)` is p.d.f. of standard normal distribution.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: non-negative.
:param float vol: volatility. non-negative.
:return: value of third derivative.
:rtype: float.
"""
norm = scipy.stats.norm
assert(vol > 0.0)
# option is expired
if maturity < 0.0 or np.isclose(maturity, 0.0):
return 0.0
discount = math.exp(-rate * maturity)
d2 = func_d2(underlying, strike, rate, maturity, vol)
d_fprime = d_fprime_by_strike(underlying, strike, rate, maturity, vol)
d_fhess = d_fhess_by_strike(underlying, strike, rate, maturity, vol)
d2_density = norm.pdf(d2)
d2_density_fprime = mafipy.function.norm_pdf_fprime(d2)
term1 = d2_density_fprime * d_fprime * d_fprime
term2 = d2_density * d_fhess
return -discount * (term1 + term2)
# ----------------------------------------------------------------------------
# Black scholes greeks
# ----------------------------------------------------------------------------
[docs]def black_scholes_call_delta(underlying, strike, rate, maturity, vol):
"""black_scholes_call_delta
calculates black scholes delta.
.. math::
\\frac{\partial}{\partial S} c(S, K, r, T, \sigma)
= \Phi(d_{1}(S))
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is volatility,
:math:`\Phi` is standard normal c.d.f,
:math:`d_{1}` is defined in
:py:func:`func_d1`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: if maturity <= 0, this function returns 0.
:param float vol: volatility. This must be positive.
:return: value of delta.
:rtype: float.
"""
assert(vol >= 0.0)
if maturity <= 0.0:
return 0.0
d1 = func_d1(underlying, strike, rate, maturity, vol)
return scipy.stats.norm.cdf(d1)
[docs]def black_scholes_call_gamma(underlying, strike, rate, maturity, vol):
"""black_scholes_call_gamma
calculates black scholes gamma.
.. math::
\\frac{\partial^{2}}{\partial S^{2}} c(S, K, r, T, \sigma)
= -\phi(d_{1}(S, K, r, T, \sigma))
\\frac{1}{S^{2}\sigma\sqrt{T}}
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is volatility,
:math:`\Phi` is standard normal c.d.f,
:math:`d_{1}` is defined in
:py:func:`func_d1`.
See :py:func:`black_scholes_call_value`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity:
if maturity is not positive, this function returns 0.0.
:param float vol: volatility. This must be positive.
:return: value of gamma.
:rtype: float.
"""
assert(vol >= 0.0)
if maturity <= 0.0:
return 0.0
d1 = func_d1(underlying, strike, rate, maturity, vol)
denominator = underlying * vol * math.sqrt(maturity)
return scipy.stats.norm.pdf(d1) / denominator
[docs]def black_scholes_call_vega(underlying, strike, rate, maturity, vol):
"""black_scholes_call_vega
calculates black scholes vega.
.. math::
\\frac{\partial}{\partial \sigma} c(S, K, r, T, \sigma)
= \sqrt{T}S\phi(d_{1}(S, K, r, T, \sigma))
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is volatility,
:math:`\phi` is standard normal p.d.f,
:math:`d_{1}` is defined in
:py:func:`func_d1`.
See :py:func:`black_scholes_call_value`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: if maturity <= 0.0, this function returns 0.
:param float vol: volatility. This must be positive.
:return: value of vega.
:rtype: float.
"""
assert(vol >= 0.0)
if maturity <= 0.0:
return 0.0
d1 = func_d1(underlying, strike, rate, maturity, vol)
return math.sqrt(maturity) * underlying * scipy.stats.norm.pdf(d1)
[docs]def black_scholes_call_volga(underlying, strike, rate, maturity, vol):
"""black_scholes_call_volg
calculates black scholes volga.
.. math::
\\frac{\partial^{2}}{\partial \sigma^{2}} c(S, K, r, T, \sigma)
S \phi^{\prime}(d_{1}(\sigma))
\\frac{
(\\frac{1}{2} \sigma^{2} - r)T
}{
\sigma^{2}
}
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is volatility,
:math:`\phi` is standard normal p.d.f,
:math:`d_{1}` is defined in
:py:func:`func_d1`.
See :py:func:`black_scholes_call_value`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: must be non-negative.
:param float vol: volatility. This must be positive.
:return: value of volga.
:rtype: float.
"""
assert(vol >= 0.0)
if maturity < 0.0:
return 0.0
d1 = func_d1(underlying, strike, rate, maturity, vol)
pdf_fprime = mafipy.function.norm_pdf_fprime(d1)
ln_moneyness = math.log(underlying / strike)
numerator = -ln_moneyness + (0.5 * vol * vol - rate) * maturity
factor = numerator / (vol * vol)
return underlying * pdf_fprime * factor
[docs]def black_scholes_call_theta(underlying, strike, rate, maturity, vol, today):
"""black_scholes_call_theta
calculates black scholes theta.
.. math::
\\frac{\partial}{\partial t} c(t, S, K, r, T, \sigma)
= - S * \phi(d_{1})
\left(
\\frac{\sigma}{2\sqrt{T - t}}
\\right)
- r e^{-r(T - t)} K \Phi(d_{2})
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is volatility,
:math:`\phi` is standard normal p.d.f,
:math:`d_{1}` is defined in
:py:func:`func_d1`.
See :py:func:`black_scholes_call_value`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: must be non-negative.
:param float vol: volatility. This must be positive.
:return: value of theta.
:rtype: float.
"""
assert(maturity >= 0.0)
assert(vol >= 0.0)
norm = scipy.stats.norm
time = maturity - today
d1 = func_d1(underlying, strike, rate, time, vol)
d2 = func_d2(underlying, strike, rate, time, vol)
term1 = underlying * norm.pdf(d1) * (vol / (2.0 * math.sqrt(time)))
term2 = rate * math.exp(-rate * time) * strike * norm.cdf(d2)
return - term1 - term2
[docs]def black_scholes_call_rho(underlying, strike, rate, maturity, vol, today):
"""black_scholes_call_rho
calculates black scholes rho.
.. math::
\\frac{\partial}{\partial t} c(t, S, K, r, T, \sigma)
= (T - t)
e^{-r (T - t)}
K \Phi(d_{2})
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is volatility,
:math:`\phi` is standard normal p.d.f,
:math:`d_{2}` is defined in
:py:func:`func_d2`.
See :py:func:`black_scholes_call_value`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: must be non-negative.
:param float vol: volatility. This must be positive.
:return: value of rho.
:rtype: float.
"""
assert(maturity >= 0.0)
assert(vol >= 0.0)
norm = scipy.stats.norm
time = maturity - today
d2 = func_d2(underlying, strike, rate, time, vol)
return time * math.exp(-rate * time) * strike * norm.cdf(d2)
[docs]def black_scholes_call_vega_fprime_by_strike(
underlying, strike, rate, maturity, vol):
"""black_scholes_call_vega_fprime_by_strike
calculates derivative of black scholes vega with respect to strike.
This is required for :py:func:`sabr_pdf`.
.. math::
\\frac{\partial}{\partial K}
\mathrm{Vega}{\mathrm{BSCall}}(S, K, r, T, \sigma)
=
S\phi^{\prime}(d_{1}(S, K, r, T, \sigma))
\\frac{1}{\sigma K}
where
:math:`S` is underlying,
:math:`K` is strike,
:math:`r` is rate,
:math:`T` is maturity,
:math:`\sigma` is volatility,
:math:`\phi` is standard normal p.d.f,
:math:`d_{1}` is defined in
:py:func:`func_d1`.
See :py:func:`black_scholes_call_value`.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity: if maturity <= 0.0, this function returns 0.
:param float vol: volatility. This must be positive.
:return: derivative of vega with respect to strike.
:rtype: float.
"""
assert(vol >= 0.0)
if maturity <= 0.0:
return 0.0
d1 = func_d1(underlying, strike, rate, maturity, vol)
density_fprime = mafipy.function.norm_pdf_fprime(d1)
return -underlying * density_fprime / (vol * strike)
# ----------------------------------------------------------------------------
# Black scholes distributions
# ----------------------------------------------------------------------------
[docs]def black_scholes_cdf(underlying, strike, rate, maturity, vol):
"""black_scholes_cdf
calculates value of c.d.f. of black scholes model.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity:
:param float vol: must be positive.
:return: value of p.d.f. of black scholes model.
:rtype: float.
"""
assert(vol > 0.0)
return (1.0
+ black_scholes_call_value_fprime_by_strike(
underlying,
strike,
rate,
maturity,
vol) * math.exp(rate * maturity))
[docs]def black_scholes_pdf(underlying, strike, rate, maturity, vol):
"""black_scholes_pdf
calculates value of p.d.f. of black scholes model.
:param float underlying:
:param float strike:
:param float rate:
:param float maturity:
:param float vol: must be positive.
:return: value of p.d.f. of black scholes model.
:rtype: float.
"""
assert(vol > 0.0)
return (black_scholes_call_value_fhess_by_strike(
underlying,
strike,
rate,
maturity,
vol) * math.exp(rate * maturity))