some_option_pricing_exercise.py

# -*- coding: utf-8 -*-
import math
import numpy
import matplotlib.pyplot
import scipy.stats

def d1(S0, K, r, sigma, T):
   return (math.log(S0 / K) + (r + sigma**2 / 2) * T) / (sigma * math.sqrt(T))

def d2(S0, K, r, sigma, T):
   return (math.log(S0 / K) + (r - sigma**2 / 2) * T) / (sigma * math.sqrt(T))

def BlackScholesEuropeanCallPrice(S0, K, r, sigma, T):
   return S0 * scipy.stats.norm.cdf(d1(S0, K, r, sigma, T)) - K * math.exp(-r * T) * scipy.stats.norm.cdf(d2(S0, K, r, sigma, T))

def BlackScholesEuropeanPutPrice(S0, K, r, sigma, T):
   return K * math.exp(-r * T) * scipy.stats.norm.cdf(-d2(S0, K, r, sigma, T)) - S0 * scipy.stats.norm.cdf(-d1(S0, K, r, sigma, T))

S0 = 100.0
K = 110.0
r = 0.03
sigma = 0.1
T = 0.5

print ("S0\tstock price at time 0:", S0)
print ("K\tstrike price:", K)
print ("r\tcontinuously compounded risk-free rate:", r)
print ("sigma\tvolatility of the stock price per year:", sigma)
print ("T\ttime to maturity in trading years:", T)

c_BS = BlackScholesEuropeanCallPrice(S0, K, r, sigma, T)
p_BS = BlackScholesEuropeanPutPrice(S0, K, r, sigma, T)

print
print ("c_BS\tBlack-Scholes European call price:", c_BS)
print ("p_BS\tBlack-Scholes European put price:", p_BS)

#picture
"""
S0s = numpy.arange(90.0, 120.0, 0.5)
cs_1D = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=1.0 / 252.0) for S0 in S0s]
cs_3M = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=0.25) for S0 in S0s]
cs_6M = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=0.5) for S0 in S0s]
cs_1Y = [BlackScholesEuropeanCallPrice(S0, K, r, sigma, T=1.0) for S0 in S0s]
matplotlib.pyplot.plot(S0s, cs_1D, 'r-', lw=2)
matplotlib.pyplot.plot(S0s, cs_3M, 'y-', lw=2)
matplotlib.pyplot.plot(S0s, cs_6M, 'b-', lw=2)
matplotlib.pyplot.plot(S0s, cs_1Y, 'g-', lw=2)
matplotlib.pyplot.xlabel("$S_0$")
matplotlib.pyplot.ylabel("$c$")
matplotlib.pyplot.title("European call price versus $S_0$, $T = \\frac{1}{252}, \\frac{1}{4}, \\frac{1}{2}, 1$")
matplotlib.pyplot.grid(True)
matplotlib.pyplot.show()
"""

#unusable part
def BinomialTreeEuropeanCallPrice(S0, K, r, sigma, T, N=30):
   deltaT = T / N

   u = math.exp(sigma * math.sqrt(deltaT))
   d = 1.0 / u
   # Initialise our f_{i,j} tree with zeros
   fs = [[0.0 for j in range(i + 1)] for i in range(N + 1)]

   a = math.exp(r * deltaT)

   p = (a - d) / (u - d)
   oneMinusP = 1.0 - p

   # Compute the leaves, f_{N, j}
   for j in range(i+1):
       fs[N][j] = max(S0 * u**j * d**(N - j) - K, 0.0)

   for i in range(N-1, -1, -1):
       for j in range(i + 1):
           fs[i][j] = math.exp(-r * deltaT) * (p * fs[i + 1][j + 1] + oneMinusP * fs[i + 1][j])

   # print fs

   return fs[0][0]

c_BT = BinomialTreeEuropeanCallPrice(S0, K, r, sigma, T)
print ("c_BT\tBinomial tree European call price:", c_BT)
#unusable part

def MonteCarloEuropeanCallPrice(S0, K, r, sigma, T, N=100, pathCount=5000):
   deltaT = T / N

   fs = [0.0 for i in range(pathCount)]

   for i in range(pathCount):
       S = S0
       epsilons = scipy.stats.norm.rvs(size=N)
       for j in range(N):
           S += r * S * deltaT + sigma * S * epsilons[j] * math.sqrt(deltaT)
       fs[i] = max(S - K, 0.0)

   return numpy.mean(fs)

c_MC = MonteCarloEuropeanCallPrice(S0, K, r, sigma, T)
print ("c_MC\tMonte Carlo European call price:", c_MC)

This is a nearly perfect option pricing program. It didn't take me much time to write it at the beginning, but I combined multiple pricing into one program, and I can show the difference between different pricing methods in a graphical way.
The only bug is that the code of the binomial tree model seems to have some problems and cannot be priced correctly.