The Black Scholes model estimates the value of a European call or put option by using the following parameters:

S = Stock Price

K = Strike Price at Expiration 

r = Risk-free Interest Rate

T = Time to Expiration

sig = Volatility of the Underlying asset

Using R, we can write a function to compute the option price once we have the values of these 5 parameters. The BlackScholes function takes these parameters and returns the value of the call or put option.

BlackScholes <- function(S, K, r, T, sig, type){
  
  if(type=="C"){
  d1 <- (log(S/K) + (r + sig^2/2)*T) / (sig*sqrt(T))
  d2 <- d1 - sig*sqrt(T)
  
  value <- S*pnorm(d1) - K*exp(-r*T)*pnorm(d2)
  return(value)}
  
  if(type=="P"){
  d1 <- (log(S/K) + (r + sig^2/2)*T) / (sig*sqrt(T))
  d2 <- d1 - sig*sqrt(T)
  
  value <-  (K*exp(-r*T)*pnorm(-d2) - S*pnorm(-d1))
  return(value)}
}

This function first calculates the d1 and d2 parameters required for the Black-Scholes model  and then uses pnorm() command of R which simulates a cumulative normal distribution. In order to use the BlackScholes function to value a call and a put option, we can run the following lines:

call <- BlackScholes(110,100,0.04,1,0.2,"C")
put <- BlackScholes(110,100,0.04,1,0.2,"P")

Some basic assumption that this model has is that assumes that the volatility is constant over the life of the option. There are two ways to estimate the volatility that is plugged in the Black Scholes model. The first approach is by calculating the historical standard deviations of the stock returns and the second method uses the option market prices to find the Implied Volatility or the volatility that the market estimate for the stock.

As an example we can download data for AAPL stock from Quandl and compute the standard deviation of the returns for a 50 days period.

# Get the Close column(4) from the WIKI dataset from Quandl of the APPL stock 

data <- Quandl("WIKI/AAPL.4")

# Retrieve the first 50 values

recent_data <- data[1:50,]

# Use the arrange function from dplyr package to get old values at top. This would guarantee that the calculations are performed with the oldest values of the series at first.

recent_data <- arrange(recent_data, -row_number())

# Add a new column with the returns of prices. Input NA to the first value of the column in order to fit the column in the recent_data data frame.

recent_data$returns <- c('NA',round(diff(log(recent_data$Close)),3))

# Convert the column to numeric

recent_data$returns <- as.numeric(recent_data$returns)
 
# Calculate the standard deviation of the log returns

Standard_deviation <- sd(recent_data$returns,na.rm=TRUE)

The result of this is a standard deviation of the returns of 0.01807. If we are dealing with Time until Expiration and Interest Rate on a year basis we need to convert the standard deviation measure to year values. This can be done with the following formula using 250 trading days in a year:

annual_sigma <- standard_deviation * sqrt(250)

The annual standard deviation is 0.2857366.

The second approach to calculate volatility is to find the Implied Volatility value that the market estimates for a specific option. We can use the Black Scholes model to calculate the Implied Volatility by using known values of the stock price, strike price, time to expiration, interest rate and an array of standard deviation values that would allow us to find an option price which is the nearest price related to the option market price. 

Suppose that the parameters of the Black Scholes model are the following:

S <- 50

K <- 70

C <- 1.5 (Actual Call Price)

r <- 0

Time <- 1

sigma <- seq(0.01,0.5,0.005)

The variable CallValueperSigma would store all the option prices that are calculated with the different values of sigma along the sigma vector.

# The variable CallValueperSigma would store all the option prices that are 
# calculated with the different values of sigma along the sigma vector. 

CallValueperSigma <- BlackScholes(S,K,r,Time,sigma) 

The next step is to find which of these values is the nearest compared to the actual call price which is 1.5.

# IV is the the min value of the difference between each value of the  #CallValueperSigma  vector and the call price. 

IV <- which.min(abs(CallValueperSigma-C))/2 

#The “which” command would return the index of the value in the sigma vector that #correspond to the option price derived from the BlackScholes  function  that is #nearest respect  the call market price. 
#This value is divided by two in order to get the specific Implied Volatility value from the sigma vector, because sigma starts from 0 to 0.5.

The Implied Volatility value from  above is 32%. So we can conclude that the market assign a 32 % of Implied Volatility to this contract option.