With the programs below, I price European calls and puts using several options pricing models. This was a definitely a pedagogical exercise: at some point: you have to admit that pricing in 0.005 seconds is just good enough regardless of the intricacies of the model. Regardless, I go through multiple so as to compare runtime and implement in Julia.
Let's consider a stock with underlying price $S$ with volatility $σ$ and an option with strike price $K$ and maturity $T$. The risk-free rate is $r$.
S = 100.0 # Underlying price
σ = 0.2 # Volatility
K = 100.0 # Strike price
T = 1.0 # Maturity
r = 0.05; # Risk-free rate
The cop-out.
We assume that $S$ follows a geometric Brownian motion so that $$dS = μS dt + σS dW,$$ where $μ$ is the drift, $σ$ is the volatility, and $W$ is a standard Brownian motion.
Everything happens from this stochastic differential equation. Since we use the risk-neutral measure for pricing, we replace the actual drift $μ$ with $r$. We then discretize the continuous-time geometric Brownian motion into $N$ time steps, in the form $$S(t+dt) = S(t) * \exp((r-\frac{1}{2}σ^2)dt + σ \sqrt{dt} * \epsilon),$$ where $\epsilon$ is a standard normal random variable.
Recall how option payoff works: a call makes $\max(S_T - K, 0)$ and and a put makes $\max(K - S_T, 0)$, where $S_T$ is the price at $T$. Option price is then given by $$p = \exp(-rT) * E[\text{Payoff}].$$
How do we determine $E[\text{Payoff}]$, you ask? Easy. We run $M$ price paths, each with $N$ time steps, noting each resulting price and the resultant option payoff. Then we just average and discount. This is incredibly simple to implement and very accurate, at the cost of computation.
include("mc.jl");
N = 1000 # Time steps
M = 100000 # Paths
MC.call(S, K, T, r, σ, N, M)
MC.put(S, K, T, r, σ, N, M)
print("Call option: ")
@time MCcall_price = MC.call(S, K, T, r, σ, N, M)
MCcall_memory = @allocated MC.call(S, K, T, r, σ, N, M)
print("Put option: ")
@time MCput_price = MC.put(S, K, T, r, σ, N, M)
MCput_memory = @allocated MC.put(S, K, T, r, σ, N, M)
println("Call price using Monte-Carlo Simulation: ", MCcall_price)
println("Put price using Monte-Carlo Simulation: ", MCput_price)
Call option: 0.557122 seconds (14 allocations: 781.609 KiB) Put option: 0.642294 seconds (14 allocations: 781.609 KiB) Call price using Monte-Carlo Simulation: 10.427693444670505 Put price using Monte-Carlo Simulation: 5.593117474536605
The supposed industry standard.
Go read Black and Scholes, 1973 "The Pricing of Options and Corporate Liabilities" if you really want.
We again assume $S$ follows a geometric Brownian motion with constant drift and volatility (that SDE will show up quite a lot).
The formula for the price of a call option $C$ is $$C = S * N(d_+) - K * \exp(-rT) * N(d_-),$$ and the formula for the price of a put option $P$ is $$P = K * \exp(-rT) * N(-d_-) - S * N(-d_+),$$ where $N(x)$ is the cumulative distribution function of the standard normal distribution and $$d_+=\frac{1}{\sigma\sqrt{T}}(\ln(\frac{S}{K})+(r+\frac{\sigma^2}{2})T),$$ $$d_-=d_+ - \sigma\sqrt{T}.$$
include("bs.jl");
BS.call(S, K, T, r, σ)
BS.put(S, K, T, r, σ)
print("Call option: ")
@time BScall_price = BS.call(S, K, T, r, σ)
BScall_memory = @allocated BS.call(S, K, T, r, σ)
print("Put option: ")
@time BSput_price = BS.put(S, K, T, r, σ)
BSput_memory = @allocated BS.put(S, K, T, r, σ)
println("Call price using Black-Scholes: ", BScall_price)
println("Put price using Black-Scholes: ", BSput_price)
Call option: 0.000003 seconds (2 allocations: 32 bytes) Put option: 0.000004 seconds (2 allocations: 32 bytes) Call price using Black-Scholes: 10.450583572185565 Put price using Black-Scholes: 5.573526022256971
Discrete Black-Scholes.
This formalization is from Cox, Ross, and Rubinstein 1979 "Option pricing: A simplified approach."
Over a small time interval, $S$ can either move up by a factor $u$ or down by a factor $d$, with probabilities $p$ and $1-p$ respectively. As the number of time steps increases, this process converges to the continuous-time geometric Brownian motion assumed above.
Key parameters in addition to the ones above are $\delta t=\frac{T}{N}$, the length of each time step, $u=\exp(\sigma*\sqrt{\delta t})$, $d=1/u$, and $p=\frac{\exp(r*\delta t)-d}{u-d}$. Note that by construction, $u*d=1$. Essentially, we construct a tree, work backwards, and then discount.
include("crr.jl");
N = 1000 # Time steps
CRR.call(S, K, T, r, σ, N)
CRR.put(S, K, T, r, σ, N)
print("Call option: ")
@time CRRcall_price = CRR.call(S, K, T, r, σ, N)
CRRcall_memory = @allocated CRR.call(S, K, T, r, σ, N)
print("Put option: ")
@time CRRput_price = CRR.put(S, K, T, r, σ, N)
CRRput_memory = @allocated CRR.put(S, K, T, r, σ, N)
println("Call price using Cox-Ross-Rubinstein: ", CRRcall_price)
println("Put price using Cox-Ross-Rubinstein: ", CRRput_price)
Call option: 0.000222 seconds (4 allocations: 16.281 KiB) Put option: 0.000226 seconds (4 allocations: 16.281 KiB) Call price using Cox-Ross-Rubinstein: 10.44858410376464 Put price using Cox-Ross-Rubinstein: 5.571526553833646
include("jr.jl");
N = 1000 # Time steps
JR.call(S, K, T, r, σ, N)
JR.put(S, K, T, r, σ, N)
print("Call option: ")
@time JRcall_price = JR.call(S, K, T, r, σ, N)
JRcall_memory = @allocated JR.call(S, K, T, r, σ, N)
print("Put option: ")
@time JRput_price = JR.put(S, K, T, r, σ, N)
JRput_memory = @allocated JR.put(S, K, T, r, σ, N)
println("Call price using Jarrow-Rudd: ", JRcall_price)
println("Put price using Jarrow-Rudd: ", JRput_price)
Call option: 0.000218 seconds (4 allocations: 16.281 KiB) Put option: 0.000220 seconds (4 allocations: 16.281 KiB) Call price using Jarrow-Rudd: 10.452179348628865 Put price using Jarrow-Rudd: 5.57513513189215
include("lr.jl");
N = 1000 # Time steps
LR.call(S, K, T, r, σ, N)
LR.put(S, K, T, r, σ, N)
print("Call option: ")
@time LRcall_price = LR.call(S, K, T, r, σ, N)
LRcall_memory = @allocated LR.call(S, K, T, r, σ, N)
print("Put option: ")
@time LRput_price = LR.put(S, K, T, r, σ, N)
LRput_memory = @allocated LR.put(S, K, T, r, σ, N)
println("Call price using Leisen-Reimer: ", LRcall_price)
println("Put price using Leisen-Reimer: ", LRput_price)
Call option: 0.000220 seconds (4 allocations: 16.281 KiB) Put option: 0.000277 seconds (4 allocations: 16.281 KiB) Call price using Leisen-Reimer: 10.44858410376464 Put price using Leisen-Reimer: 5.571526553833646