Adv Calc HW 4 | Python Fiddle

# -*- coding: utf-8 -*-
"""
Created on Tue Feb 24 21:23:39 2015

@author: Jawndalf
"""

# -*- coding: utf-8 -*-
"""
Created on Tue Feb 24 13:26:32 2015

@author: Jawndalf
"""

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return np.exp((-x**2.0)/2.0)

def simp_rule(a, b, n):
    h = (b-a)/float(n)
    total = f(a) / 6.0 + f(b) / 6.0
    
    for i in range (1, n):
        total = total + f(a + i*h) / 3.0
        
    for i in range (1, n+1):
        total = total + 2 * f(a + (i-0.5) * h) / 3.0
        
    return total * h

def cum_distrib(x, n):
    total = 0.5
    return total + 1/(np.sqrt(2.0 * np.pi)) * simp_rule(0.0, x, n)

def cum_distrib_approxnn(t):
    z = np.abs(t)
    y = 1.0 / (1 + 0.2316419 * z)
    a1 = 0.319381530
    a2 = -0.356563782
    a3 = 1.781477937
    a4 = -1.821255978
    a5 = 1.330274429
    
    m = 1 - np.exp((-t**2)/2.0) * (a1*y+ a2*y**2 + a3*y**3 + a4*y**4 + a5*y**5) / np.sqrt(2*np.pi)
    if t > 0:
        nn = m
    else:
        nn = 1-m
    
    return nn

def cum_distrib_tol(tol, x):
    tolerance = float(tol)
    n = 4
    old_int = cum_distrib(x,n)
    n = 2 * n
    new_int = cum_distrib(x,n)
    
    while abs(new_int - old_int) > tolerance:
        old_int = new_int
        n = 2* n
        new_int = cum_distrib(x,n)
    
    return new_int
"""return (n, new_int, abs(new_int - old_int) ) ==> copy this in if you want to see how many iterations it takes, and what the difference between the new and old estimates are"""

""" Black Scholes Theorem"""

def PutCallOptionPrice (S, K, r, q, vol, T, t):
    d1 = ( np.log(float(S) / K) + (r - q + (vol**2)/2) * (T - t) ) / (vol * np.sqrt(T - t))
    d2 = d1 - vol * np.sqrt(T - t)
    
    C = S* np.exp(-q* (T-t)) * cum_distrib_tol(10**-12, d1) - K * np.exp(-r *(T-t)) * cum_distrib_tol(10**-12, d2)
    P = K* np.exp(-r* (T-t)) * cum_distrib_tol(10**-12, -d2) - S * np.exp(-q *(T-t)) * cum_distrib_tol(10**-12, -d1)
    
    delta = np.exp(-q*(T-t)) * cum_distrib_tol(10**-12,d1)
    
    return (C, delta, cum_distrib_tol(10**-12,d1),cum_distrib_tol(10**-12,d2))
    
"""return (C, P, d1, d2)
return (C, P, d1, cum_distrib_tol(10**-12,-d1))"""

def PutCallOptionPrice_appx(S, K, r, q, vol, T, t):
    d1 = ( np.log(float(S) / K) + (r - q + (vol**2)/2) * (T - t) ) / (vol * np.sqrt(T - t))
    d2 = d1 - vol * np.sqrt(T - t)
    
    C = S* np.exp(-q* (T-t)) * cum_distrib_approxnn(d1) - K * np.exp(-r *(T-t)) * cum_distrib_approxnn(d2)
    P = K* np.exp(-r* (T-t)) * cum_distrib_approxnn(-d2) - S * np.exp(-q *(T-t)) * cum_distrib_approxnn(-d1)
    
    
    return (C, P, d1, d2,cum_distrib_approxnn(d1),cum_distrib_approxnn(d2))

print PutCallOptionPrice(50.0, 45.0, 0.03, 0.01, 0.3, 5/12.0, 0.0)

print PutCallOptionPrice_appx(40.0, 40.0, 0.05, 0.01, 0.2, 0.25, 0.0)

Python Fiddle

Python Cloud IDE