In this report, we are going to implement theory and exercices from John Hull (2009): Fundamentals of Futures and Options Markets (7th edition)

1 The Lognormal Distribution (p.316)

Black-Scholes assumes that \(log(S_T)\) are normally distributed

\[log(S_T) \sim N(log(S_0) + (\mu- \frac{\sigma^2}{2})T , \sigma\sqrt{T} )\]

\[\frac{dS}{S} \sim N(\mu\delta t , \sigma\sqrt{\delta t} )\]

where \(S_0\) is the current stock price

1.1 Example 1 (p.317):

Consider a stock with \(S_0\) = 40, expected return = 16% p.a. and volatility = 20% p.a.

The probability distribution of \(S_T\) in six months is:

S0 <- 40
mu <- 16/100 
sd <- 20/100
T <- 0.5

(mean <- log(S0) + (mu-sd^2/2)*T)
## [1] 3.76
(var   <- sd^2*T)
## [1] 0.02
(stdev <- sqrt(var))
## [1] 0.141
# CI of ln(S_T) 
(CI <- c(mean-1.96*stdev, mean+1.96*stdev))
## [1] 3.48 4.04
# Plot distribution
x <- seq(3,5,by = .01)
dens <- dnorm(x,mean = mean, sd = stdev)
plot(x     , dens , type='l', main= "distribution of ln(S_t)")
i <- x >= 3.48 & x <= 4.04
polygon(c(3.48,x[i],4.04), c(0,dens[i],0), col="lightblue") 

# This implies CI of S_T = 
(CI2 <- exp(CI))
## [1] 32.5 56.6
plot(exp(x), dens , type='l', main= "distribution of S_t")
i <- exp(x) >= 32.5 & exp(x) <= 56.6
polygon(c(32.5,exp(x[i]),56.6), c(0,dens[i],0), col="lightblue")

1.2 Example 2 (p.318)

Consider a stock with - expected return = 17% p.a. - volatility = 20% p.a.

The pdf for the rate of return (continuously compounded) over one year is normal

(mu <- 0.17 - 0.2^2/2)
## [1] 0.15
(sd <- 0.20) # percent
## [1] 0.2
# CI for returns = 
c(mu-1.96*sd, mu+1.96*sd)
## [1] -0.242  0.542
x <- seq(-0.5,0.8,by = .01)
dens <- dnorm(x,mean = mu, sd = sd)
plot(x     , dens , type='l', main= "distribution of ln(S_t)")
i <- x >= -0.242 & x <= 0.542
polygon(c(-0.242,x[i],0.542), c(0,dens[i],0), col="lightblue")

2 Pricing functions

According to the Black-Scholes model, the following functions give the price of plain vanilla european call and put options respectively:

\[c=S_0e^{-qT}N(d_1) - Ke^{-rT}N(d_2)\]

\[p=Ke^{-rT}N(-d_2) - S_0e^{-rT}N(-d_1)\] where

\[d_1 = \frac{ln(S_0/K) + (r-q+\sigma^2/2)T }{ \sigma\sqrt{T} }\]

\[d_2 = \frac{ln(S_0/K) + (r-q-\sigma^2/2)T }{ \sigma\sqrt{T}} = d_1 - \sigma\sqrt{T}\]

When the stock price \(S_0\) becomes large, \(d_1\) and \(d_2\) become very large, and

  • \(N(d_1)\) and \(N(d_2)\) are both close to \(1\)
  • Call almost certain to be exercised (price close to \(c = S_0 - Ke^{-rT}\) )
  • \(N(-d_1)\) and \(N(-d_2)\) are both close to \(0\)
  • Put price close to \(0\)

When the stock price \(S_0\) becomes small, \(d_1\) and \(d_2\) become very large negative, and

  • \(N(d_1)\) and \(N(d_2)\) are both close to \(0\)
  • Call price becomes \(0\)
  • \(N(-d_1)\) and \(N(-d_2)\) are both close to \(1\)
  • put price becomes \(Ke^{-rT} - S_0\)

Here are the functions implemented in R:

Call <- function(S, K, r, T, sigma) {
  d1  <-  (log(S/K) + (r + sigma^2/2)*T) / (sigma*sqrt(T))
  d2  <-  d1 - sigma*sqrt(T)
  S * pnorm(d1)  - K*exp(-r*T)*pnorm(d2)
}

Put <- function(S, K, r, T, sigma) {
  d1  <-  (log(S/K) + (r + sigma^2/2)*T) / (sigma*sqrt(T))
  d2  <-  d1 - sigma*sqrt(T)
  -S * pnorm(-d1) + K*exp(-r*T)*pnorm(-d2)
}

Now, we have to set some parameters:

Variable Description Example
S Spot Price 100
K Strike Price 70
r risk-free interest rate 0.05
T Time to Maturity 0.5 => 6 months
sigma Implied annual volatility 0.16 => 1% daily volatility 
# Parameters
S     <- 100
K     <- 70
r     <- 0
T     <- 1       # seq(10/52,1/52,-1/52)
sigma <- c(.16,.21,.3, .4,.5,.6,.75,1,1.06)

Some intermediate calculations:

d1  <-  (log(S/K) + (r + sigma^2/2)*T) / (sigma*sqrt(T))
d2  <-  d1 - sigma*sqrt(T)
p   <-  cbind(S, K, T, sigma,  p = -S * pnorm(-d1) + K*exp(-r*T)*pnorm(-d2))  
c   <-  cbind(S, K, T, sigma,  c =  S * pnorm(d1)  - K*exp(-r*T)*pnorm(d2))

Probability of exercise:

pnorm(d2) %>% round(2)
## [1] 0.98 0.94 0.85 0.76 0.68 0.62 0.54 0.44 0.42

Expected value of payoff at maturity:

K*pnorm(d2) %>% round(2)
## [1] 68.6 65.8 59.5 53.2 47.6 43.4 37.8 30.8 29.4

3 Plotting \(\theta\)

\(\theta\) is the time value of an option. Here we are looking at

cprct <- c[,c(1,5)] / matrix(rep(c[1,c(1,5)],each=length(c[,5]) ),ncol = 2) * 100
ts.plot(cprct, col = c(1,4), main = "Evolution of Stocks and Calls in %")
legend('topleft', c(paste("Stock +", tail(cprct,1)[1]-100, "%") , 
                    paste("Long Call +", round(tail(cprct,1)[2])-100, "%")), 
       col=c(1,4), lty=1)

pprct <- p[,c(1,5)] / matrix(rep(p[1,c(1,5)],each=length(c[,5]) ),ncol = 2) * 100
ts.plot(pprct, col=1:2, main = "Evolution of Stocks and Puts in %")

4 Put-Call Parity

c[,5] + K*exp(-r*T)
## [1] 100 100 101 103 106 108 112 119 121
p[,5] + S %>% round(1)
## [1] 100 100 101 103 106 108 112 119 121

5 Find Implicit Volatility

S <- 40
C <- 1.50
VperSigma <- Call(S,K,r,T,seq(0.01,0.50,0.005))
IV <- which.min(abs(VperSigma-C))/2
print(paste("Implicit Volatility is around", IV, "%)"))
## [1] "Implicit Volatility is around 47.5 %)"

Parameters:

mu <- 0
sd <- 0.2

x <- rlnorm(10000, mu,sd)
plot(density(x), xlim=c(0,4)) ; abline(v=c(1, 1.5), lty=2)

dlnorm(1.2, mu,sd,log = T)   # pdf
## [1] 0.0927
plnorm(1.54, mu,sd)          # CDF
## [1] 0.985
qlnorm(0.95, mu,sd)          # Quantile Function
## [1] 1.39

Plotting:

par(mfrow = c(1, 2))
curve(dnorm, from = -3, to = 3, n = 1000, main = "Normal Probability Density Function")
text(-2, 0.3, expression(f(x) == paste(frac(1, sqrt(2 * pi * sigma^2)),
                                       " ", e^{
                                         frac(-(x - mu)^2, 2 * sigma^2)
                                       })), cex = 1.2)
x <- dnorm(seq(-3, 3, 0.001))
plot(seq(-3, 3, 0.001), cumsum(x)/sum(x), type = "l", col = "blue",
     xlab = "x", main = "Normal Cumulative Distribution Function")
text(-1.5, 0.7, expression(phi(x) == paste(frac(1, sqrt(2 * pi)),
                                           " ", integral(e^(-t^2/2) * dt, -infinity, x))), cex = 1.2)