import math
import scipy
import scipy.stats
import numpy as np

result = 0.0
cdf = 0.0

"""

*** This module includes the blackScholes implementation as well as the Greeks.  We will need one of the Greeks to 
implement our Implied Volatility module.

This is an implementation to calculate Cummulative Density Function (CDF) and Option Price (in Black-Scholes Model)
using the Pseudocode from the Primer of Mathematics for Financial Engineering by Dan Stefanica. I will check to see how
the results compare to the ones generated by Matlab's Financial Instrument ToolBox at a later time. But for now this is just an
implementation of the code from the above book. Black-Scholes is an analytical formula and good for valuing instruments
that cannot be exercised before maturity such as European Call or Put. To value options such as American Put or call
that can be exercised prior to maturity, you may have to use binomial models or trees. Check out Financial Engineering
and Risk Management on Coursera, the course gave a good treatment of these topics (Thank goodness for these free courses,
you can just learn on a continous basis (as time permits) and don't have to pay a dime). I will post codes for building
binomial lattices later (but it may take a long while).

COMPARISON OF CODE VERSUS MATLAB FINANCIAL INSTRUMENT TOOLBOX
You can use Matlab to price options as well by calling like so:
[Call, Put] = blsprice(Price, Strike, Rate, Time, Volatility, Yield)
So if you replace the parameters above with the following values you get the following prices which converges to the
ones generated by this python module:
[Call, Put] = blsprice(42, 40, 0.05, 0.5, 0.3, 0.03)

Call = 4.7053

Put = 2.3430


If you run this code in python, like so
(blackScholes(0,42.0,40.0,0.5,0.3,0.05,0.03,'call'))
(blackScholes(0,42.0,40.0,0.5,0.3,0.05,0.03,'put'))

You get same results as above:
The value of a call with the given parameters is : 4.7053

The value of a put with the given parameters is : 2.3430



"""
def cumm_dens_function(t):
    z = abs(t)
    y = 1/(1 + 0.2316419*z)
    a1 = 0.319381530
    a2 = -0.356563782
    a3 = 1.781477937
    a4 = -1.821255978
    a5 = 1.330274429
    
    cdf = 1 - math.exp(-t**2/2)*(a1*y + a2*y**2 + a3*y**3 + a4*y**4+a5*y**5)/math.sqrt(2*math.pi)
    
    if t > 0:
        result = cdf
    else:
        result = 1 - cdf
        
    return result

"""
This is an alternative implementation to calculate CDF using scipy.
This implementation is more straight forward.
"""
def cumm_dens_function_scipy(t):
    return scipy.stats.norm.cdf(t)

"""
Black-Scholes for European Call and Put
"""

def blackScholes(t,S,K,T,a,r,q,Type):
    #t = beginning time
    #S = Spot Price
    #K = Strike
    #T = Maturity
    #a = volatility
    #r = constant interest rate
    #q = continous dividend rate of the underlying asset
    #Type can be call or put - enter as strings e.g. 'call' or 'put'
    call = 0.0
    put = 0.0
    d1 =0.0
    d2 = 0.0
    S = float(S)
    K = float(K)
    
    d1 = (np.log(S/K) +(r - q + a**2/2)*(T - t))/(a*math.sqrt(T-t))
    d2 = d1 - (a*math.sqrt(T-t))
    
    call = S*math.exp(-q*(T-t))*cumm_dens_function_scipy(d1) - K*math.exp(-r*(T-t))*cumm_dens_function_scipy(d2)
    put = K*math.exp(-r*(T-t))*cumm_dens_function_scipy(-d2) - S*math.exp(-q*(T-t))*cumm_dens_function_scipy(-d1)
    
    if Type == 'call':
        return call
    else:
        return put

#Sample Calls:
#(blackScholes(0,42.0,40.0,0.5,0.3,0.05,0.03,'call'))
#(blackScholes(0,42.0,40.0,0.5,0.3,0.05,0.03,'put'))
        

def deltaCallPut(t,S,K,T,a,r,q,Type):
    #t = beginning time
    #S = Spot Price
    #K = Strike
    #T = Maturity
    #a = volatility
    #r = constant interest rate
    #q = continous dividend rate of the underlying asset
    #Type can be call or put - enter as strings e.g. 'call' or 'put'
    d1 = (np.log(S/K) +(r - q + a**2/2)*(T - t))/(a*math.sqrt(T-t))
    d2 = d1 - (a*math.sqrt(T-t))
    
    if Type == 'call':
        return math.exp(-q*(T-t))*cumm_dens_function_scipy(d1)
    else:
        return -1*math.exp(-q*(T-t))*cumm_dens_function_scipy(-d1)

#Sample Calls:
#deltaCallPut(0,42.0,40.0,0.5,0.3,0.05,0.03,'call')    
#deltaCallPut(0,42.0,40.0,0.5,0.3,0.05,0.03,'put')

def gammaCallPut(t,S,K,T,a,r,q):
    #t = beginning time
    #S = Spot Price
    #K = Strike
    #T = Maturity
    #a = volatility
    #r = constant interest rate
    #q = continous dividend rate of the underlying asset
    #Type can be call or put - enter as strings e.g. 'call' or 'put'
    d1 = (np.log(S/K) +(r - q + a**2/2)*(T - t))/(a*math.sqrt(T-t))
    d2 = d1 - (a*math.sqrt(T-t))
    return math.exp(-q*(T-t))/(S*a*math.sqrt(T-t)) * (1.0/(math.sqrt(2*math.pi))*math.exp(-d1**2/2.0))

#Sample Calls:
#gammaCallPut(0,42.0,40.0,0.5,0.3,0.05,0.03)

def vegaCallPut(t,S,K,T,a,r,q):
    #t = beginning time
    #S = Spot Price
    #K = Strike
    #T = Maturity
    #a = volatility
    #r = constant interest rate
    #q = continous dividend rate of the underlying asset
    #Type can be call or put - enter as strings e.g. 'call' or 'put'
    d1 = (np.log(S/K) +(r - q + a**2/2)*(T - t))/(a*math.sqrt(T-t))
    return S * math.exp(-q*(T-t)) * math.sqrt(T-t) * (1.0/(math.sqrt(2*math.pi))*math.exp(-d1**2/2.0))
    
#Sample Calls:
#vegaCallPut(0,42.0,40.0,0.5,0.3,0.05,0.03)

def thetaCallPut(t,S,K,T,a,r,q,Type):
    #t = beginning time
    #S = Spot Price
    #K = Strike
    #T = Maturity
    #a = volatility
    #r = constant interest rate
    #q = continous dividend rate of the underlying asset
    #Type can be call or put - enter as strings e.g. 'call' or 'put'
    d1 = (np.log(S/K) +(r - q + a**2/2)*(T - t))/(a*math.sqrt(T-t))
    d2 = d1 - (a*math.sqrt(T-t))
    
    if Type == 'call':
        return ((-1*S*a*math.exp(-q*(T-t)))/(2.0*math.sqrt(2*math.pi*(T-t))) * math.exp(-d1**2/2.0)) + q * S * math.exp(-q*(T-t)) * cumm_dens_function_scipy(d1) - r*K*math.exp(-r*(T-t))*cumm_dens_function_scipy(d2)
    else:
        return ((-1*S*a*math.exp(-q*(T-t)))/(2*math.sqrt(2*math.pi*(T-t))) * math.exp(-d1**2/2.0)) - q * S * math.exp(-q*(T-t)) * cumm_dens_function_scipy(-d1) + r*K*math.exp(-r*(T-t))*cumm_dens_function_scipy(-d2)
    
#Sample Calls:
#thetaCallPut(0,42.0,40.0,0.5,0.3,0.05,0.03,'call')
#thetaCallPut(0,42.0,40.0,0.5,0.3,0.05,0.03,'put')

def rhoCallPut(t,S,K,T,a,r,q,Type):
    #t = beginning time
    #S = Spot Price
    #K = Strike
    #T = Maturity
    #a = volatility
    #r = constant interest rate
    #q = continous dividend rate of the underlying asset
    #Type can be call or put - enter as strings e.g. 'call' or 'put'
    d1 = (np.log(S/K) +(r - q + a**2/2)*(T - t))/(a*math.sqrt(T-t))
    d2 = d1 - (a*math.sqrt(T-t))
    
    if Type == 'call':
        return K * (T-t) * math.exp(-r*(T-t)) * cumm_dens_function_scipy(d2)
    else:
        return -1 * K * (T-t) * math.exp(-r*(T-t)) * cumm_dens_function_scipy(-d2)
 

