Black Scholes Implementation

Full Script

import math from scipy.stats import norm def call(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t) d2 = d1 - sigma * math.sqrt(t) u = norm.cdf(d1) * S - norm.cdf(d2) * K * math.exp(-r * t) return float(u); def put(S, K, sigma, r, t): #based on put-call parity result = K * math.exp(-1*r*t) - S + call(S, K, sigma, r, t) return result; def call_delta(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t) return float(norm.cdf(d1)) def putDelta(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t) return -1*float(norm.cdf(-1*d1)) def callRho(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t) d2 = d1 - sigma * math.sqrt(t) val = K*t*math.exp(-1*r*t)*float(norm.cdf(d2)) return val; def putRho(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t) d2 = d1 - sigma * math.sqrt(t) val = -1*K*t*math.exp(-1*r*t)*float(norm.cdf(-1*d2)) return val; def callVega(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t); d2 = d1 - sigma * math.sqrt(t); return S*Nprime(d1)*math.sqrt(t); def Nprime(x): return math.exp(-1*x*x/2)/math.sqrt(2*math.pi) def putVega(S, K, sigma, r, t): return callVega(S, K, sigma, r, t) def callTheta(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t) d2 = d1 - sigma * math.sqrt(t) factor1 = -1*S*sigma*(1/(2*math.sqrt(t)))*Nprime(d1) factor2 = -1*r*K*math.exp(-1*r*t)*float(norm.cdf(d2)) return factor1 + factor2 def putTheta(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t) d2 = d1 - sigma * math.sqrt(t) factor1 = -1*S*sigma*(1/(2*math.sqrt(t)))*Nprime(-1*d1) factor2 = r*K*math.exp(-1*r*t)*float(norm.cdf(-1*d2)) return factor1 + factor2 def callGamma(S, K, sigma, r, t): d1 = (1 / (sigma * math.sqrt(t))) * (math.log(S / K) + (r + (sigma * sigma) / 2) * t); return Nprime(d1)/(sigma*math.sqrt(t)*S) def putGamma(S, K, sigma, r, t): return callGamma(S, K, sigma, r, t)