from FMNM.BS_pricer import BS_pricer
from FMNM.Parameters import Option_param
from FMNM.Processes import Diffusion_process
import numpy as np
import scipy.stats as ss
from scipy.integrate import quad
from functools import partial
import matplotlib.pyplot as plt
%matplotlib inline
European option¶
Under the Black-Scholes (BS) model, the best method to price a vanilla European option is to use the BS closed formula.
The BS formula for a call is:
$$ C(t,T,S,K,r,\sigma) = S \, N(d_1) - K e^{-r(T-t)} N(d_2) $$with
$$ d_1 = \frac{1}{\sigma \sqrt{T-t}} \biggl[ \log \biggl( \frac{S}{K} \biggr) + \biggl(r + \frac{\sigma^2}{2} \biggr) (T-t) \biggr] \quad \text{and} \quad d_2 = d_1 - \sigma \sqrt{T-t} $$where $N$ is the cumulative distribution function of a standard normal random variable.
The formula for a put is similar and can be found in the wiki page.
The value of an option can be also computed as the discounted expectation of a future payoff in this way:
$$\begin{aligned} C(s,K,t, T) &= e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}\biggl[ (S_T - K)^+ \bigg| S_t=s \biggr] \\ &= e^{-r(T-t)} \int_0^{\infty}\, (s' - K)^+ \, f(s'|s) \, ds' \\ &= e^{-r(T-t)} \int_K^{\infty}\, (s' - K) \, f(s'|s) \, ds' \end{aligned}$$where $f(s'|s)$ is the risk neutral transition probability of the process $\{S_u\}_{u\in [t,T]}$ starting at $S_t=s$. This is a log-normal density function
$$ f(s'|s) = \frac{1}{s' \sigma \sqrt{2\pi (T-t)}} \; e^{- \frac{ \biggl[\log(s') - \bigl(\log(s) + (r-\frac{1}{2} \sigma^2)(T-t) \bigr) \biggr]^2}{2\sigma^2 (T-t)}} $$Comment:¶
Usually in statistics random variables are indicated by capital letters, and non-random variables are indicated with small letters. However, it is very common to indicate the stock price in the Black-Scholes formula with a capital letter e.g wiki notation. I'm going to use a capital letters to indicate random variables, and small letters for non random. However there will be exceptions where it will be clear from the context e.g. strike K and maturity T are non random, or $S_0$ most of the time is a fixed variable.
Let us analyze better the previous formulas¶
\begin{align*} C(s,t,K,T) &= e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}\biggl[ (S_T - K)^+ \bigg| S_t=s \biggr] \\ &= e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}\biggl[ S_T \mathbb{1}_{S_T >K} \bigg| S_t=s \biggr] - e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}\biggl[ K \mathbb{1}_{S_T >K} \bigg| S_t=s \biggr] \\ &= e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}\biggl[ S_T \mathbb{1}_{S_T >K} \bigg| S_t=s \biggr] - e^{-r(T-t)} K \, \underbrace{\mathbb{Q}\biggl[ S_T >K \bigg| S_t=s \biggr]}_{N(d_2)} \\ \end{align*}Let us introduce the following change of measure (under the stock numeraire):
$$ \begin{aligned} \frac{d \tilde{\mathbb{Q}} }{ d \mathbb{Q} } &= \frac{S_T}{\mathbb{E}^\mathbb{Q}[S_T | S_t]} = \frac{S_T}{S_t e^{r(T-t)}} \\ &= \frac{S_t e^{(r -\frac{1}{2}\sigma^2)(T-t) + \sigma W_{T-t}} }{S_t e^{r(T-t)}} \\ &= e^{ -\frac{1}{2}\sigma^2(T-t) + \sigma W_{T-t} } \quad \text{(exponential martingale)} \end{aligned} $$By Girsanov theorem, under $\tilde{\mathbb{Q}}$ the driving Brownian motion has the new dynamics
$$ \tilde{W_t} = W_t - \sigma t $$and the corresponding stock dynamics becomes
$$\begin{aligned} \frac{dS_t}{S_t} &= r dt + \sigma dW_t \\ &= (r+\sigma^2) dt + \sigma d\tilde{W}_t \end{aligned}$$The first term is
$$ \begin{aligned} e^{-r(T-t)} \mathbb{E}^{\mathbb{Q}}\biggl[ S_T \mathbb{1}_{S_T >K} \bigg| S_t=s \biggr] =& e^{-r(T-t)} \mathbb{E}^{\tilde{\mathbb{Q}}} \biggl[ \frac{d \mathbb{Q} }{ d \tilde{\mathbb{Q}}} S_T \mathbb{1}_{S_T >K} \bigg| S_t=s \biggr] \\ =& e^{-r(T-t)} \mathbb{E}^{\tilde{\mathbb{Q}}} \biggl[ \frac{e^{r(T-t)}S_t}{S_T} S_T \mathbb{1}_{S_T > K} \bigg| S_t=s \biggr] \\ =& \; s \, \underbrace{\tilde{\mathbb{Q}} ( S_T > K | S_t=s)}_{N(d_1)} \end{aligned}$$We have just seen how to interpret the terms $N(d_1)$ and $N(d_2)$. These are the risk neutral probabilities of $S_T > K$ in the stock and money market numeraires respectively.
I implemented the BS closed formula in the class BS_pricer.
Let us consider the following set of parameters, that will be recurrent in all the next notebooks.
S0 = 100.0 # spot stock price
K = 100.0 # strike
T = 1.0 # maturity
r = 0.1 # risk free rate
sig = 0.2 # diffusion coefficient or volatility
call = BS_pricer.BlackScholes("call", S0, K, T, r, sig)
put = BS_pricer.BlackScholes("put", S0, K, T, r, sig)
print("Call price: ", call)
print("Put price: ", put)
Call price: 13.269676584660893 Put price: 3.753418388256833
print(call)
print(put + S0 - K * np.exp(-r * T))
13.269676584660893 13.269676584660886
Numerical integration¶
I want to play a bit with the different formulas written above.
Let us compute the option prices by integrating the log-normal density:
e_ret = np.log(S0) + (r - 0.5 * sig**2) * T # expected return of the log-price
vol = sig * np.sqrt(T) # standard deviation of the log-price
# log-normal density (defined above)
def log_normal(x, e_ret, vol):
return 1 / (x * vol * np.sqrt(2 * np.pi)) * np.exp(-((np.log(x) - e_ret) ** 2) / (2 * vol**2))
plt.figure(figsize=(8, 5))
x = np.linspace(30, 190, 100)
plt.plot(x, log_normal(x, e_ret, vol))
plt.title("Lognormal distribution, conditioned on $S_0=100$")
plt.xlabel("$S_T$")
plt.show()
The function log_normal(x, e_ret, vol) defined above, corresponds to the scipy.stats function:
ss.lognorm.pdf(x, vol, scale=np.exp(e_ret) ).
In the next calculation, I'm going to use the scipy function.
Let us perform the integration with the scipy.integrate function quad:
def integrand_LN(S, strike, e_ret, vol, payoff):
if payoff == "call":
return (S - strike) * ss.lognorm.pdf(S, vol, scale=np.exp(e_ret))
elif payoff == "put":
return (strike - S) * ss.lognorm.pdf(S, vol, scale=np.exp(e_ret))
call = quad(integrand_LN, K, np.inf, args=(K, e_ret, vol, "call"))[0] * np.exp(-r * T)
put = quad(integrand_LN, 0, K, args=(K, e_ret, vol, "put"))[0] * np.exp(-r * T)
print("Call price: {} \nPut price: {} ".format(call, put))
Call price: 13.269676584660926 Put price: 3.753418388256828
The put option payoff $(K-S_T)^+$ is positive for $S_T < K$.
- In the call case, the integration is from $K$ to $\infty$.
- In the put case, the integration is from $0$ to $K$.
What if we use the change of measure proposed above? In this way the integrations are simpler.
Let us compute $\tilde{\mathbb{Q}}( S_T > K )$ and $\mathbb{Q}( S_T > K )$.
# expected return of the log-price under the new measure
e_ret_1 = np.log(S0) + (r + 0.5 * sig**2) * T
Q1 = quad(lambda S: ss.lognorm.pdf(S, vol, scale=np.exp(e_ret_1)), K, np.inf)[0]
print("Risk neutral probability under stock numeraire,\n Q1 =", Q1)
Q2 = quad(lambda S: ss.lognorm.pdf(S, vol, scale=np.exp(e_ret)), K, np.inf)[0]
print("Risk neutral probability under money market numeraire,\n Q2 =", Q2)
print("BS call price: ", S0 * Q1 - K * np.exp(-r * T) * Q2)
Risk neutral probability under stock numeraire, Q1 = 0.7257468822499276 Risk neutral probability under money market numeraire, Q2 = 0.6554217416103069 BS call price: 13.26967658466257
It is quite common to compute the Black-Scholes formula using $N(d_1)$ and $N(d_2)$.
The reason is that the cumulative function of the standard Normal distribution is more accessible (I guess). In the BS_pricer class I used the function scipy.stats.norm.cdf.
For completeness, let me recall that if $X_T$ is a Normal random variable, then $S_T = S_0 e^{X_T}$ is Log-Normal. Therefore we have:
$$ \mathbb{Q}( S_T > K ) = \mathbb{Q}\biggl( S_0 e^{X_T} > K \biggr) = \mathbb{Q}\biggl( X_T > \log \frac{K}{S_0} \biggr). $$This permits to use the Normal cumulative function.
Monte Carlo method¶
I'm going to simulate the random variables:
$$ S_T^i = S_0 e^{(r -\frac{1}{2}\sigma^2)T + \sigma W_{T}^i} $$for $1 \leq i \leq N$.
Then use the approximation for a call option:
For a put option I use this payoff $(K - S_T )^+$ inside the expectation.
np.random.seed(seed=44) # seed for random number generation
N = 10000000 # Number of random variables
W = ss.norm.rvs((r - 0.5 * sig**2) * T, np.sqrt(T) * sig, N)
S_T = S0 * np.exp(W)
call = np.mean(np.exp(-r * T) * np.maximum(S_T - K, 0))
put = np.mean(np.exp(-r * T) * np.maximum(K - S_T, 0))
call_err = ss.sem(np.exp(-r * T) * np.maximum(S_T - K, 0)) # standard error
put_err = ss.sem(np.exp(-r * T) * np.maximum(K - S_T, 0)) # standard error
print("Call price: {}, with error: {}".format(call, call_err))
print("Put price: {}, with error: {}".format(put, put_err))
Call price: 13.263338006636662, with error: 0.005093638687208466 Put price: 3.7553894006350093, with error: 0.002214066662789331
BS_pricer¶
In the next notebooks I will present better the class BS_pricer. But now let's have a look at the prices obtained by different pricing methods:
# Creates the object with the parameters of the option
opt_param = Option_param(S0=100, K=100, T=1, exercise="European", payoff="call")
# Creates the object with the parameters of the process
diff_param = Diffusion_process(r=0.1, sig=0.2)
# Creates the pricer object
BS = BS_pricer(opt_param, diff_param)
BS.closed_formula()
13.269676584660893
BS.Fourier_inversion()
13.269676584660623
BS.MC(N=30000000, Err=True, Time=True)
# output is: price, standard error and execution time
(array([13.26753511]), array([0.00294085]), 0.679786205291748)
BS.mesh_plt() # PDE method
The PDE approach is the topic of the notebook 2.1.
Binomial tree¶
Of course I cannot forget about the Binomial model! This is a simple but very powerful numerical method!
I expect you to be familiar with this model. If not, have a look at the wiki page.
Although I said I expect you to know the model, I'm not expecting that you have already implemented it!
Therefore, here I present an efficient implementation:
N = 15000 # number of periods or number of time steps
payoff = "call" # payoff
dT = float(T) / N # Delta t
u = np.exp(sig * np.sqrt(dT)) # up factor
d = 1.0 / u # down factor
V = np.zeros(N + 1) # initialize the price vector
# price S_T at time T
S_T = np.array([(S0 * u**k * d ** (N - k)) for k in range(N + 1)])
a = np.exp(r * dT) # risk free compounded return
p = (a - d) / (u - d) # risk neutral up probability
q = 1.0 - p # risk neutral down probability
if payoff == "call":
V[:] = np.maximum(S_T - K, 0.0)
else:
V[:] = np.maximum(K - S_T, 0.0)
for i in range(N - 1, -1, -1):
# the price vector is overwritten at each step
V[:-1] = np.exp(-r * dT) * (p * V[1:] + q * V[:-1])
print("BS Tree Price: ", V[0])
BS Tree Price: 13.269537371978052
The stock price on each node of a binomial tree is $$ S_n^{(k)} = u^{k} d^{n-k} S_0 \quad \text{for} \quad k \in \{ 0, ..., n \} \quad \text{and} \quad n \in \{ 0, ..., N \}$$ where the price at $S_n^{(k)}$ is obtained starting from $S_0^{(0)}$ := $S_0$ and applying $k$ up factors and $(n-k)$ down factors.
In a (recombining) binomial tree at each time point $n$ there are $i+1$ nodes. The total number of nodes of a binomial tree is
\begin{align} \sum_{i=0}^N (1+i) &= \sum_{i=0}^N 1 + \sum_{i=0}^N i \\ &= (N+1) + N(N+1)/2 \\ &= (N+1)(N+2)/2. \end{align}Since $u = 1/d$, we can rewrite the formula above as $$ S_n^{(k)} = u^{2k - n}\, S_0 \quad \text{for} \quad k \in \{ 0, ..., n \} \quad \text{and} \quad n \in \{ 0, ..., N \}$$
In this way, we can decrease the number of operations and increase the speed of the program:
%%timeit -n 20 -r 10
S_T = np.array([(S0 * u**j * d ** (N - j)) for j in range(N + 1)]) # price S_T at time T
2.64 ms ± 37.9 µs per loop (mean ± std. dev. of 10 runs, 20 loops each)
%%timeit -n 20 -r 10
S_0N = S0 / u**N
S_T = np.array([S_0N * u ** (2 * j) for j in range(N + 1)]) # price S_T at time T
1.59 ms ± 24.1 µs per loop (mean ± std. dev. of 10 runs, 20 loops each)
%%timeit -n 20 -r 10
S_0N = S0 / u**N
S_T = np.array([S_0N * u ** (j) for j in range(0, 2 * N + 1, 2)]) # price S_T at time T
1.49 ms ± 25 µs per loop (mean ± std. dev. of 10 runs, 20 loops each)
The last approach for the computation of $S_T$ is the fastest.
Limits of the model¶
BS_sigma = partial(BS_pricer.BlackScholes, "call", S0, K, T, r) # binding the function
sigmas = np.linspace(0.01, 10, 1000)
plt.plot(sigmas, BS_sigma(sigmas))
plt.xlabel("sig")
plt.ylabel("price")
plt.title("Black-Scholes price as function of volatility")
plt.show()
The BS formula is an increasing function of the volatility.
However, for higher volatilities, the graph becomes almost flat!!
We can conclude that the model is reliable for volatilities in the range $0 - 400\%$.