Options valuation

.OPTIONS.md

justin@justin-XPS-13-9360:~/work$ python3 black_scholes.py 
call: 4.75942
put : 0.80860
justin@justin-XPS-13-9360:~/work$ python3 binomial.py 
call: 4.75982
put : 0.80899
justin@justin-XPS-13-9360:~/work$ python3 monte_carlo.py 
call: 4.75913
put : 0.80862

binomial.py

import math

def option_value(fn, S, X, sigma, R, T, n=1000):
    dt=T/n
    a=math.exp(R*dt)
    u=math.exp(sigma*(dt**0.5))
    d=1/u
    p=(a-d)/(u-d)
    q=[fn(S*u**i*d**(n-i), X)
       for i in range(n+1)]            
    for i in range(n-1, -1, -1):
        q=[math.exp(-R*dt)*(p*q[j+1]+(1-p)*q[j])
           for j in range(i+1)]
    return q[0]

def call_value(S, X, sigma, R, T):
    return option_value(lambda S, X: max(0, S-X),
                        S, X, sigma, R, T)

def put_value(S, X, sigma, R, T):
    return option_value(lambda S, X: max(0, X-S),
                        S, X, sigma, R, T)

if __name__=="__main__":
    print ("call: %.5f" % call_value(42, 40, 0.2, 0.1, 0.5))
    print ("put : %.5f" % put_value(42, 40, 0.2, 0.1, 0.5))

black_scholes.py

import scipy.stats as ss

import math

CDF=ss.norm.cdf

def option_value(fn, S, X, sigma, R, T):
    sx=math.log(S/X)
    st=sigma*(T**0.5)
    d1=(sx+(R+sigma**2/2)*T)/st
    d2=(sx+(R-sigma**2/2)*T)/st
    df=math.exp(-R*T)
    return fn(S, X, d1, d2, df)

def call_value(S, X, sigma, R, T):
    return option_value(lambda S, X, d1, d2, df: S*CDF(d1)-X*df*CDF(d2),
                        S, X, sigma, R, T)

def put_value(S, X, sigma, R, T):
    return option_value(lambda S, X, d1, d2, df: X*df*CDF(-d2)-S*CDF(-d1),
                        S, X, sigma, R, T)
    
if __name__=="__main__":
    print ("call: %.5f" % call_value(42, 40, 0.2, 0.1, 0.5))
    print ("put : %.5f" % put_value(42, 40, 0.2, 0.1, 0.5))

monte_carlo.py

import numpy.random as npr

import numpy as np

def init_deviates(paths, var):
    q=npr.normal(size=paths)*var
    m=np.mean(q)
    return q-m    

def option_value(fn, S, X, sigma, R, T, seed=12345, paths=1000000):
    npr.seed((seed, seed))
    q=init_deviates(paths, sigma*(T**0.5))
    mu=(R-0.5*sigma**2)*T
    tv=S*np.exp(q+mu)
    payoff=fn(tv, X)
    df=np.exp(-R*T)
    return np.mean(payoff)*df

def call_value(S, X, sigma, R, T):
    return option_value(lambda S, X: np.maximum(0, S-X),
                        S, X, sigma, R, T)

def put_value(S, X, sigma, R, T): 
    return option_value(lambda S, X: np.maximum(0, X-S),
                        S, X, sigma, R, T)
    
if __name__=="__main__":
    print ("call: %.5f" % call_value(42, 40, 0.2, 0.1, 0.5))
    print ("put : %.5f" % put_value(42, 40, 0.2, 0.1, 0.5)) 

requirements.txt

numpy
scipy