Calculate Implied Volatility of an option price given its market price

Calculate_Implied_Vol.py

from numpy import *
from scipy.stats import norm
from scipy.optimize import root, fsolve, newton
import logging
 
logging.basicConfig(format='%(asctime)s %(levelname)s %(message)s', level=logging.DEBUG)
 
 
def BlackScholesCall(S, K, T, sigma, r = 0., q = 0.):
        #logging.debug("S=" + str(S) + ",K=" +  str(K) + ",T=" +  str(T) + ",sigma=" +  str(sigma) + ",r=" + str(r) + ",q=" +  str(q) )
        d1 = (1/(sigma*sqrt(T)))*(log(S/K) + (r + sigma**2)*T)
        d2 = d1 - sigma*sqrt(T)
        return S*exp(-q*T)*norm.cdf(d1) - K*exp(-r*T)*norm.cdf(d2)
    
def BlackScholesPut(S, K, T, sigma, r = 0., q = 0.):
        d1 = (1/(sigma*sqrt(T)))*(log(S/K) + (r + sigma**2)*T)
        d2 = d1 - sigma*sqrt(T)
        return  (-S*exp(-q*T)*norm.cdf(-d1)) + (K*exp(-r*T)*norm.cdf(-d2))
 
def Delta(sigma, S, K, T,r, q, pricer):
        h = 1.e-3
        return (pricer(S + h, K, T, sigma, r, q) - pricer(S - h, K, T, sigma -h, r, q))/2*h
    
def Gamma(sigma, S, K, T,r, q, pricer):
        h = 1.e-2
        return (pricer(S + h, K, T, sigma , r, q) - 2* pricer(S , K, T, sigma -h, r, q) + pricer(S -h, K, T, sigma -h, r, q))/h**2
    
def Vega(sigma, S, K, T,r, q, pricer):
        h = 1.e-3
        return (pricer(S, K, T, sigma + h, r, q) - pricer(S, K, T, sigma -h, r, q))/2*h  


def Theta(sigma, S, K, T,r, q, pricer):
        h = 1.e-3
        return (pricer(S, K, T + h, sigma, r, q) - pricer(S, K, T -h, sigma, r, q))/2*h
    
option_err = lambda sigma, option_price, S, K, T, r, q, pricer: abs(pricer(S, K, T, sigma, r, q) - option_price)

def GetImpVol(sigma0, call_price, S, K, T, r, q, pricer, method ):
        dict = {"root":root, "fsolve": fsolve}        
        return dict[method]( option_err , sigma0, args = (call_price, S, K, T, r, q, pricer))  
    
#def GetImpVol_newton(sigma0, call_price, S, K, T, r, q, pricer):              
#        return newton( option_err , sigma0, fprime = Vega, args = (call_price, S, K, T, r, q, pricer))  
    
#vec_GetImpVol = vectorize(GetImpVol_newton)  

if __name__ == "__main__":
        no_opt = 2 
        S = 100.*ones(no_opt)
        K = linspace(80,120,no_opt)
        sigma = linspace(0.25,0.35,no_opt)
        T = 1.*ones(no_opt)
        r = 0.02*ones(no_opt)
        q = 0.01*ones(no_opt)
 
        call_prices = BlackScholesCall(S, K, T, sigma, r, q )
        put_prices = BlackScholesPut(S, K, T, sigma, r, q )
        logging.debug(call_prices)
        logging.debug(put_prices)
        
        logging.debug("Original Sigma: " + str(sigma))
        logging.debug("Call Price: " + str(call_prices))
 
        sigma0 = 0.1*ones(no_opt)
        
        sigma_result = GetImpVol(sigma0,call_prices,S,K,T,r,q,BlackScholesCall, "root")
        call_prices_result = BlackScholesCall(S, K, T, sigma_result.x, r, q  )
        logging.debug("Sigma Result: " + str(sigma_result.x))        
        logging.debug("Call Price Result: " + str(call_prices_result))
        
        sigma_result = GetImpVol(sigma0,call_prices,S,K,T,r,q,BlackScholesCall, "fsolve")
        call_prices_result = BlackScholesCall(S, K, T, sigma_result, r, q  )
        logging.debug("Sigma Result: " + str(sigma_result))
        logging.debug("Call Price Result: " + str(call_prices_result))