Black-Scholes

Black-Scholes model We first create a function implementing the Black-Scholes formula and related Greeks:

BScall <- function(t=0,T,S,K,r,q=0,sigma,isPut=0) {
# t and T are measured in years; all parameters are annualized
# q is the continuous dividend yield
d1 <- (log(S/K)+(r-q+sigma^2/2)*(T-t))/(sigma*sqrt(T-t))
d2 <- d1-sigma*sqrt(T-t)
binary <- pnorm(-d2)*exp(-r*T)

# Call Delta at t
Delta <- exp(-q*(T-t))*pnorm(d1)
Gamma <- exp(-q*(T-t))*exp(-d1^2/2)/sqrt(2*pi)/S/sigma/sqrt(T-t)
Vega <- S*exp(-q*(T-t))/sqrt(2*pi)*exp(-d1^2/2)*sqrt(T-t)
Theta <- -S*exp(-q*(T-t))*sigma/sqrt(T-t)/2*dnorm(d1) - r*K*exp(-r*(T-t))*pnorm(d2) + 
    q*S*exp(-q*(T-t))*pnorm(d1)
Rho <- (T-t)*K*exp(-r*(T-t))*pnorm(d2)


# Black-Scholes formula for Calls
BSprice <- -K*exp(-r*(T-t))*pnorm(d2)+S*Delta

if (isPut==1) {
    Delta <- -exp(-q*(T-t))*pnorm(-d1)
    BSprice <- S*Delta+K*exp(-r*(T-t))*pnorm(-d2)
    Theta <- -S*exp(-q*(T-t))*sigma/sqrt(T-t)/2*dnorm(d1) + r*K*exp(-r*(T-t))*pnorm(-d2) - 
    q*S*exp(-q*(T-t))*pnorm(-d1)
    Rho <- -(T-t)*K*exp(-r*(T-t))*pnorm(-d2)
}
Bank <- BSprice-Delta*S

return (list(Delta=Delta,Gamma=Gamma,Theta=Theta,Vega=Vega,Rho=Rho,Price=BSprice,d1=d1,d2=d2,B=Bank))
}

Here is an example usage for a Call option (Call maturity of 20 weeks, initial share price \(S_0=49\) ):

BScall(t=0,T=20/52,S=49,K=50,r=0.05,q=0,sigma=0.2,isPut=0)
## $Delta
## [1] 0.5216047
## 
## $Gamma
## [1] 0.06554404
## 
## $Theta
## [1] -4.30533
## 
## $Vega
## [1] 12.10548
## 
## $Rho
## [1] 8.906962
## 
## $Price
## [1] 2.400527
## 
## $d1
## [1] 0.05418135
## 
## $d2
## [1] -0.06985338
## 
## $B
## [1] -23.1581

And another (in the context of Delta-hedging, where we change \(t,S_t\), keeping the rest of the inputs fixed):

BScall(t=1/52,T=20/52,S=50,K=50,r=0.05,q=0,sigma=0.2)
## $Delta
## [1] 0.5837767
## 
## $Gamma
## [1] 0.064538
## 
## $Theta
## [1] -4.542943
## 
## $Vega
## [1] 11.7906
## 
## $Rho
## [1] 9.617237
## 
## $Price
## [1] 2.867973
## 
## $d1
## [1] 0.2115647
## 
## $d2
## [1] 0.09067058
## 
## $B
## [1] -26.32086

Simulation of stock prices

We simulate paths of Brownian motion with drift, \(\{B(t)\}\).

To do so we repeatedly draw standard normal increments.

Below the drift is \(\mu=1\) , the volatility is \(\sigma =1\) and the initial condition is \(B(0)=0\)

We generate 2 independent paths with \(t\in [0,10]\).

dt <- 0.01
path_of_BM <- rep(0, 1001)
for (i in 2:1001)
   path_of_BM[i] <- path_of_BM[i-1] + 1*dt + 1*sqrt(dt)*rnorm(1)

# Plot the first path and set the plot parameters
plot(seq(0,10,by=dt), path_of_BM, type="l", lwd=2, xlim=c(0, 10),
     xlab='Time', ylab='Brownian Motion w/Drift')

# Simulate and plot the second path
path_of_BM_2 <- path_of_BM
for (i in 2:1001)
   path_of_BM_2[i] <- path_of_BM_2[i-1] + 1*dt + 2*sqrt(dt)*rnorm(1)
lines(seq(0,10,by=dt), path_of_BM_2, col="blue", lwd=2)

# Add a legend to the plot
legend("topleft", legend=c("Path 1", "Path 2"), col=c("black", "blue"), lty=1, lwd=2)

dt <- 0.01 sets the time step for the simulation to 0.01. This means that each time the loop runs, it updates the value of the Brownian motion at a time 0.01 seconds later than the previous step.

path_of_BM <- rep(0, 1001) initializes an array of length 1001 to store the simulated values of the Brownian motion. The first value is set to 0, since we start the process at time 0.

The first loop for (i in 2:1001) updates the values in path_of_BM using the formula for Brownian motion with drift: \(B(t+\Delta t) = B(t) + \mu\Delta t + \sigma\sqrt{\Delta t}Z\), where \(Z\) is a standard normal random variable. The sqrt(dt)*rnorm(1) term generates a single standard normal random variable, which is multiplied by the volatility \(\sigma=1\) and the square root of the time step \(\sqrt{dt}\) to get the increment in the Brownian motion.

The plot() function plots the first simulated path in black, using solid lines of width 2. The x-axis is labeled “Time” and the y-axis is labeled “Brownian Motion w/Drift”. The limits of the x-axis are set to [0, 10].

The second loop for (i in 2:1001) updates a second array path_of_BM_2 in the same way as path_of_BM, but with a different value of the volatility: \(\sigma=2\). This generates a second, independent path of Brownian motion with drift. The lines() function then plots this second path in blue, using the same x-axis as the first plot.

Now, for simulating stock prices, we need to switch to Geometric Brownian motion. Similar, except we need to exponentiate to go from normal to log-normal random variables.

First, fix T and sample from the log-normal distribution of \(S_T\):

# Set the random seed for reproducibility
set.seed(170)

# Define the parameters of the model
r <- 0.05  # risk-free rate
delt <- 0  # dividend yield
S0 <- 49  # initial stock price
sigm <- 0.2  # volatility
T <- 20/52  # time horizon in years (here, 20 weeks)

# Simulate the stock prices at time T for 10,000 paths using the Black-Scholes model
S_T <- S0*exp((r-delt-0.5*sigm^2)*T + sigm*sqrt(T)*rnorm(10000))

# Create a histogram of the simulated stock prices with 100 bins and store the counts and breaks
hs <- hist(S_T, 100, plot = FALSE)

# Plot the histogram and add a title with the mean and standard deviation of the simulated stock prices
hist(S_T, 100, main = paste("Mean =", prettyNum(mean(S_T)), "sd =", prettyNum(sd(S_T))))

The histogram created above can be used to estimate the price of an option using the Monte Carlo method and the risk-neutral pricing approach. Specifically, we can estimate the expected discounted payoff of a Call option with a strike price of K = 50 by taking the average of the sampled payoffs from the simulated stock prices. In this case, we are using the Q-probability, which is already incorporated into our sampling above, to estimate the option premium. Therefore, the estimated Call premium is given by the empirical average of the 1000 sampled payoffs.

exp(-r*T)*mean( pmax(S_T - 50, 0))
## [1] 2.398787

As can be seen, this estimate is quite good, and within 1 cent of the true answer above coming from the exact Black-Scholes formula.

Below we repeat with 2 additional Monte Carlo esimates that are of much worse quality because they only use 100 samples. The first estimate is more than 30 cents off.

S_T <- S0*exp( (r-delt-0.5*sigm^2)*T + sigm*sqrt(T)*rnorm(100))
exp(-r*T)*mean( pmax(S_T - 50, 0))
## [1] 2.752859
# Try it again
S_T <- S0*exp( (r-delt-0.5*sigm^2)*T + sigm*sqrt(T)*rnorm(100))
exp(-r*T)*mean( pmax(S_T - 50, 0))
## [1] 2.522091

##Paths of GBM We use the same logic as for BM with drift to sample three representative paths of Geometric Brownian Motion.

set.seed(1)
dt <- 1/520; M <- (T/dt)+1
path_of_S <- rep(S0, M)
for (i in 2:M)
   path_of_S[i] <- path_of_S[i-1] + (r-delt)*path_of_S[i-1]*dt + sigm*sqrt(dt)*path_of_S[i-1]*rnorm(1)
plot(seq(0,T,by=dt), path_of_S, xlab='Time', ylab=expression(S[t]), type="l", lwd=2, col="blue", ylim=c(40,60))
legend("topleft", legend=c("Path 1", "Path 2", "Path 3"), col=c("blue", "red", "green"), lwd=2)

for (i in 2:M)
   path_of_S[i] <- path_of_S[i-1] + (r-delt)*path_of_S[i-1]*dt + sigm*sqrt(dt)*path_of_S[i-1]*rnorm(1)
lines(seq(0,T,by=dt), path_of_S, lwd=2, col="red")

for (i in 2:M)
   path_of_S[i] <- path_of_S[i-1] + (r-delt)*path_of_S[i-1]*dt + sigm*sqrt(dt)*path_of_S[i-1]*rnorm(1)
lines(seq(0,T,by=dt), path_of_S, lwd=2, col="green")

The final plot below super-imposes a path of the GBM against the terminal distribution of \(S_T\) , linking together the above path-based simulations and the previous Monte Carlo-based option pricing algorithm.

par(mfrow=c(1,2))
for (i in 2:M)
   path_of_S[i] <- path_of_S[i-1] + (r-delt)*path_of_S[i-1]*dt + sigm*sqrt(dt)*path_of_S[i-1]*rnorm(1)

plot(seq(0,T,by=dt), path_of_S, xlab='Time', ylab=expression(S[t]), type="l", lwd=2, col="blue", ylim=c(40,60))
barplot(hs$counts, names=hs$mids, horiz=T,  width=0.5,xlab=expression(S[T]))