Black-Scholes Formula

OK real talk — the title was a little misleading. Let me just come out and say it up front.

I cannot derive Black-Scholes. lmao

If you’re here for that, close the tab lmao

What I am going to do is code up the Black-Scholes Model. That’s the whole point. That’s all this is.

So first, let me just throw the thing at you — the Black-Scholes PDE (Partial Differential Equation), the one that literally won the Nobel Prize in Economics.

That right there is the Black-Scholes partial differential equation. The Nobel one.

But here’s the thing — we’re not actually using that one directly.

What we will use is the option pricing formula you get when you take that PDE and grind through it with some seriously heavy-duty math. That result — the Black-Scholes option pricing equation — is honestly the kind of thing that could’ve taken home a Fields Medal, which is basically the Nobel Prize of math, not of economics.

Anyway, here it is:

OK so before we just start plugging numbers in, let me walk through the assumptions baked into the Black-Scholes model. Then we’ll actually try computing a call price and a put price using those formulas above.

i) There is no arbitrage opportunity

(no way to make a riskless profit)

······· ’no arbitrage opportunity'

ii) It is possible to borrow and lend any amount

(even fractional amounts) of cash at the riskless rate

······· you can borrow or lend any amount of cash

iii) It is possible to buy and sell any amount

(even fractional amounts) of stock

······· you can also buy or sell any amount of stock

iv) None of these transactions incur any fees or costs

······· no transaction costs, no commissions — perfectly frictionless market

v) The rate of return on the riskless asset is constant

······· the riskless rate is just a fixed number, always

vi) The stock price follows geometric Brownian motion,

and we assume its drift and volatility are both constant

······· stock price follows Brownian motion with constant drift and constant volatility

(What’s the right way to think about “drift” here — like, the general prevailing direction of the stock price? I’m honestly not totally sure.)

Oh and — one thing that’s apparently a huge deal about BS: you only need the volatility of the stock. Not the expected return of the stock. That’s the key insight, supposedly… heh

Anyway. Let’s just use it.

S = 800: k = 850: r = 0.05: T = 0.5: v = 0.2:

d_1 = (log( s / k ) + (r + 0.5*v*v)*T) / (v*sqr(T))
d_2 = d_1 - v*sqr(T)

with worksheetfunction

c = S*.NormSDist(d_1) - k*exp(-r*T)*.NormSDist(d_2)

p = k*exp(-r*T)*.NormSDist(-d_2) - S*.NormSDist(-d_1)

End with

debug.print "call option price is : " & c

debug.print "put option price is : " & p

End sub

The derivation behind the formula is brutal, sure. But actually coding it up from the formula? Not that bad at all.

Let me at least think through the $N(d_1)$ piece a little.

What $N(x)$ means is:

So it’s the probability that a value drawn from the standard normal $N(0,1)$ comes out below $x$.

And $N(d_1)$ is, weirdly enough:

So $N(d_1)$ is the cumulative probability of getting a value below $d_1$ from the standard normal distribution.

But then what? T_T Because I genuinely don’t know what $d_1$ and $d_2$ are actually supposed to represent T_T

For now I just have to hold onto it and remember it like this… T_T T_T T_T

(Heads up: this screenshot was generated before the d1 fix above — the original code was missing the parentheses around v*math.sqrt(T) (Python reads d1 / v*math.sqrt(T) as (d1/v)*sqrt(T), not d1/(v*sqrt(T))), so the numbers in this image are off by a factor of T. Trust the corrected code above, not the screenshot. T_T)

import numpy as np
import scipy as sp
import scipy.stats
import math
import matplotlib.pyplot as plt
import os
path = r'C:\Users\GD Park\Desktop\financial-engineering'
os.chdir(path)
def BS_Formula(S=800, k=850, r=0.05, T=0.5, v=0.2):
    d1 = math.log(S/k) + (r + (v*v)*0.5)*T
    d1 = d1 / (v*math.sqrt(T))
    d2 = d1 - v*math.sqrt(T)
    norm = sp.stats.norm(loc=0, scale=1)
    c = S*norm.cdf(d1) - k*math.exp(-r*T)*norm.cdf(d2)
    p = k*math.exp(-r*T)*norm.cdf(-d2) - S*norm.cdf(-d1)
    return c, p
SS = np.linspace(700, 900, 1000)
C = np.zeros(len(SS))
P = np.zeros(len(SS))
for i in range(len(SS)):
    C[i], P[i] = BS_Formula(S=SS[i])
plt.figure(figsize=(6, 4))
plt.plot(SS, C)
plt.xlabel('S')
plt.ylabel('Price')
plt.title('BS_Formula_c')
plt.savefig('BS_Formula_c', dpi=200)
plt.figure(figsize=(6, 4))
plt.plot(SS, P)
plt.xlabel('S')
plt.ylabel('Price')
plt.title('BS_Formula_p')
plt.savefig('BS_Formula_p', dpi=200)