Derivatives - Functions for Black Scholes & Monte Carlo

EXECUTIVE SUMMARY

The below code creates functions and plots for Black-Scholes and Monte-Carlo Simulation for vanilla call options, european call options, and down-and-out put options.

It also examines exotic options, hedging portfolios for large price movements, associated transactions costs, and builds a Heston model.

Note : The code uses predetermine inputs for all models.

CONTENTS

i. Vanilla Call Option : Black-Scholes vs. Monte-Carlo Simulation
ii. Down-And-Out Put Option Price : by Monte-Carlo Simulation
iii. Pricing Exotic Options in Complicated Market Structures
iv. Hedging, Large Price Movements, and Transaction Costs
v. Heston Model

Vanilla Call Option : Black-Scholes vs. Monte-Carlo Simulation

stock_price = 100
strike_price = 100
r = 0.05
dividend = 0
sigma = 0.2
Time = 1/4
N = 96*90
num_sim = 5

h = Time/N
u = exp((r-dividend)*h + sigma * (h ^ 0.5))
d = exp((r - dividend)*h - sigma * (h ^ 0.5))
prob = (exp((r - dividend) * h) - d) / (u - d)

#lnSt ~ N(lnSo + (r-div - sigma^2 / 2)T, sigma^2 T)
mu =  (r - dividend - ((sigma) ** 2 / 2)) * h
std = sigma * (h ** 0.5)

Time_0 = rep(log(stock_price),num_sim)

#Parallel computing function
sim_one_node = function(last_node){
  one_node = rnorm(num_sim,mu,std) 
  new_node = last_node + one_node 
  return(new_node)
}  

ms = as.data.frame(NA)

# Monte Carlo Simulation simulates all of the 5 paths and then take the exponential
last = Time_0
ms = as.data.frame(matrix(numeric(0), ncol = 1, nrow = num_sim))
ms[1] = Time_0

for (i in 1:N) {
# run 96*90 times
   out = sim_one_node(last)
   ms[i + 1] = out
   last = out
}

ms = exp(ms) #turn back from lnSt to St

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

BlackScholes(S = 100, K = 100, r = 0.05, T = 1/4, sig = 0.2, div = 0, type = "C")

## [1] 4.614997

So the Black-Scholes call option price is 4.614997.

S=100
K=100 
tau=0.25 
r=0.05
sigma=0.2

# Compute the Monte Carlo European option price 
# Assuming the stock follows a geometric BM

set.seed(0)
N <- 96*90 #90 days and 96 5 min subintervals per day
dt <- tau/N #length of each time sub interval
time <- seq(from=0, to=tau, by=dt) #time moments in which we simulate the process

nSim <- 100 #number of simulations (paths) 

#Monte Carlo with analytic formula
Z <- rnorm(nSim, mean=0, sd=1)
WT <- sqrt(tau) * Z
ST = S*exp((r - 0.5*sigma^2)*tau + sigma*WT)
simulated_call_payoffs <- exp(-r*tau)*pmax(ST-K,0)
Call_price_MC_anal_1 <- mean(simulated_call_payoffs)

lower_bound_1 <- Call_price_MC_anal_1 - 1.96*sd(simulated_call_payoffs)/sqrt(nSim)
upper_bound_1 <- Call_price_MC_anal_1 + 1.96*sd(simulated_call_payoffs)/sqrt(nSim)

#------------------------------------------
set.seed(0)
nSim <- 1000 #number of simulations (paths) 

#Monte Carlo with analytic formula
Z <- rnorm(nSim, mean=0, sd=1)
WT <- sqrt(tau) * Z
ST = S*exp((r - 0.5*sigma^2)*tau + sigma*WT)
simulated_call_payoffs <- exp(-r*tau)*pmax(ST-K,0)
Call_price_MC_anal_2 <- mean(simulated_call_payoffs)

lower_bound_2 <- Call_price_MC_anal_2 - 1.96*sd(simulated_call_payoffs)/sqrt(nSim)
upper_bound_2 <- Call_price_MC_anal_2 + 1.96*sd(simulated_call_payoffs)/sqrt(nSim)


#------------------------------------------
set.seed(0)
nSim <- 100 #number of simulations (paths) 

#Monte Carlo with analytic formula
Z <- rnorm(nSim, mean=0, sd=1)
WT <- sqrt(tau) * Z
ST = S*exp((r - 0.5*sigma^2)*tau + sigma*WT)
simulated_call_payoffs <- exp(-r*tau)*pmax(ST-K,0)
Call_price_MC_anal_3 <- mean(simulated_call_payoffs)

lower_bound_3 <- Call_price_MC_anal_3 - 1.96*sd(simulated_call_payoffs)/sqrt(nSim)
upper_bound_3 <- Call_price_MC_anal_3 + 1.96*sd(simulated_call_payoffs)/sqrt(nSim)


#------------------------------------------
set.seed(0)
nSim <- 1000 #number of simulations (paths) 

#Monte Carlo with analytic formula
Z <- rnorm(nSim, mean=0, sd=1)
WT <- sqrt(tau) * Z
ST = S*exp((r - 0.5*sigma^2)*tau + sigma*WT)
simulated_call_payoffs <- exp(-r*tau)*pmax(ST-K,0)
Call_price_MC_anal_4 <- mean(simulated_call_payoffs)

lower_bound_4 <- Call_price_MC_anal_4 - 1.96*sd(simulated_call_payoffs)/sqrt(nSim)
upper_bound_4 <- Call_price_MC_anal_4 + 1.96*sd(simulated_call_payoffs)/sqrt(nSim)

A = matrix(NA, nrow = 4, ncol = 3)
colnames(A) = c("Option Price", "CI Lower Bound", "CI Upper Bound")
rownames(A) = c("100 times", "1,000 times", "1,000,000 times", "100,000,000 times")

A[1, 1] = Call_price_MC_anal_1
A[1, 2] = lower_bound_1
A[1, 3] = upper_bound_1

A[2, 1] = Call_price_MC_anal_2
A[2, 2] = lower_bound_2
A[2, 3] = upper_bound_2

A[3, 1] = Call_price_MC_anal_3
A[3, 2] = lower_bound_3
A[3, 3] = upper_bound_3

A[4, 1] = Call_price_MC_anal_4
A[4, 2] = lower_bound_4
A[4, 3] = upper_bound_4

print(xtable(A, digits = 3))

% latex table generated in R 4.0.2 by xtable 1.8-4 package % Wed Mar 31 18:08:57 2021

The length of the confidence interval decreases as the number of simulated paths increase . Also, the Black-Scholes call option price is 4.614997 which is very similar to the Monte-Carlo price of 4.614193.

 

Down-And-Out Put Option Price : by Monte-Carlo Simulation

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

BlackScholes(S = 100, K = 100, r = 0.05, T = 1/4, sig = 0.2, div = 0, type = "C")

## [1] 4.614997

set_up = function(r,dividend,sigma,Time,N){
  h = Time/N
  mu =  (r - dividend - ((sigma) ** 2 / 2)) * h
  std = sigma * (h ** 0.5)
  return(c(h,mu,std))
}

sim_one_node = function(num_sim,mu,std,last_node){
  one_node = rnorm(num_sim,mu,std)
  new_node = last_node + one_node
  return(new_node)
}

sim_matrix_2a = function(setup,S,K,N,num_sim,r,Time,barrier){
  h = setup[1]
  mu = setup[2]
  std = setup[3]
  Time_0 = rep(log(S),num_sim)
  last = Time_0
  prices = as.data.frame(matrix(numeric(0),ncol=1,nrow = num_sim))
  prices[,1] = Time_0
  for (i in 1:N){ 
    out = sim_one_node(num_sim,mu,std,last)
    for(j in 1:num_sim){
      if(out[j] < log(barrier) | is.na(out[j])){ #If any price < barrier, this will out, and following simulation will all out
        out[j] = NA
      }
    }
    prices[,i+1] = out
    last = out
  }
  prices = exp(prices)
  payoff_T = pmax((K - prices[,N+1]), 0)
  payoff_0 = payoff_T * exp(-r * Time)
  payoff_0[is.na(payoff_0)] = 0  #Change NA to 0
  avg_payoff = mean(payoff_0)
  NA_terms = sum(is.na(prices[,N+1]))
  answer = list("payoff" = avg_payoff, "out" = NA_terms)
  return(answer)
}

set_b = set_up(0.05,0,0.2,0.25,864) #setup
twob = sim_matrix_2a(set_b,100,95,864,100,0.05,0.25,75)#Simulate 1000 times
BS_b1 = BlackScholes(100,95,0.05,0.25,0.2,0, "C")[2] #Price of Black Scholes

#Do 20 times to get confidence interval
payoffandNA = as.data.frame(matrix(numeric(0),ncol=2,nrow = 20))
for (i in 1:20){
  print(i)
  temp = sim_matrix_2a(set_b,100,95,864,100,0.05,0.25,75)
  payoffandNA[i,1] = as.numeric(temp[1])
  payoffandNA[i,2] = as.numeric(temp[2])
}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20

colnames(payoffandNA) = c("Price","OutPath")

total_mean = mean(payoffandNA[,1]) #Mean of 20 times price
total_sd = sd(payoffandNA[,1]) #Sd of 20 times price
interval_down = total_mean - 2 * total_sd #Intervals
interval_up = total_mean + 2 * total_sd

avg_NA = mean(payoffandNA[,2]) #Mean of 20 times out path

The Black Scholes price of put is 1.5342. The average payoff of 20 simulations of Down-And-Out put price is 1.4941. The confidence interval is 1.2590 - 1.7292. The average out path of Down-And-Out put is 3.15 . The Down-And-Out put is cheaper than normal put because the payoff that is larger than 95 - 75 is excluded.

volatility = c(0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4)
out_path = c()
price = c()
for(i in volatility){
  set = set_up(0.05,0,i,0.25,8640)
  two = sim_matrix_2a(set,100,95,8640,1000,0.05,0.25,75)
  out_path = append(out_path,as.numeric(two[2]))
  price = append(price , as.numeric(two[1]))
}

As shown in the graph, when the volatility becomes larger, the out paths number increase. This is because when volatility goes up, the probability of the stock price touch the barrier also gets higher.

Pricing Exotic Options in Complicated Market Structures

#3a)
#2 dim BM
#2 risky and 1 Rf
r=.05
beta = r
alpha = .6
sigma11 = .1
sigma12 = .2
sigma21 = .3
S1 = c(10)
S2 = c(10)
delta=.1

dt=1/250   #this is the step wrt time
T=1

rate_paths = matrix(NA, nrow=251, ncol=1000)  #each col is a simulated rate path, each row is the evolution of that rate path from t=0 to T=1
rate_paths = as.data.frame(rate_paths)

rate_paths[1,] = r

for(j in 1:1000){  #for each simulated path
  for(i in 2:251){
    rate_paths[i,j] = rate_paths[i-1,j] + alpha*(beta - rate_paths[i-1,j])*dt + delta*sqrt(rate_paths[i-1,j])*sqrt(dt)*rnorm(1)
  }
}

r_T = c(as.numeric(rate_paths[251,]))
hist(r_T, breaks=20, main = "Histogram of r_T given r_0 = 0.05")

#3b
r=.05
beta = r
alpha = .6
sigma11 = .1
sigma12 = .2
sigma21 = .3
S1 = c(10)
S2 = c(10)
delta=.1

dt=(1/52)  #this is the step wrt time
T=100

rate_paths2 = matrix(NA, nrow=(52*100 +1), ncol=1)
rate_paths2 = as.data.frame(rate_paths2)

rate_paths2[1,1] = r

for(i in 2:(52*100 +1)){
  rate_paths2[i,1] = rate_paths2[i-1,1] + alpha*(beta - rate_paths2[i-1,1])*dt + delta*sqrt(rate_paths2[i-1,1])*sqrt(dt)*rnorm(1)
}

plot(x=seq(0,100,1/52), y=rate_paths2[,1], type='l', main = "Trajectory of interest rate from t=0 to T=100", xlab = "Time", ylab = "Interest rate")

#3c
K=10
T=.5


#first find the i-rate paths: r_t
r=.05
beta = r
alpha = .6
sigma11 = .1
sigma12 = .2
sigma21 = .3
S1 = c(10)
S2 = c(10)
delta=.1

dt=(1/250)  #this is the step wrt time

rate_paths3 = matrix(NA, nrow=(250*.5 +1), ncol=1000)  
rate_paths3 = as.data.frame(rate_paths3)

rate_paths3[1,] = r

for(j in 1:1000){  #for each simulated path
  for(i in 2:(250*.5 +1)){
    rate_paths3[i,j] = rate_paths3[i-1,j] + alpha*(beta - rate_paths3[i-1,j])*dt + delta*sqrt(rate_paths3[i-1,j])*sqrt(dt)*rnorm(1)
  }
}

#now find the evolution of the asset
asset1 = matrix(NA, nrow=(250*.5 +1), ncol=1000)  
asset1 = as.data.frame(asset1)

asset1[1,] = 10 #S1 like we hvae in line 313

for(j in 1:1000){  #for each simulated path
  for(i in 2:(250*.5 +1)){
    asset1[i,j] = asset1[i-1,j] + rate_paths3[i-1,j]*asset1[i-1,j]*dt + sigma11*sqrt(asset1[i-1,j])*sqrt(dt)*rnorm(1) + sigma12*(asset1[i-1,j])*sqrt(dt)*rnorm(1)
  }
}

#now we price the asset call option
K=10
T=.5
r_T = c(as.numeric(rate_paths3[126,])) #this row represents 10,000 observations of r_T at T=.5

ST = c(as.numeric(asset1[126,]))   #this row represents 10,000 observations of S_T at T=.5

simulated_call_payoffs = c()
for(i in 1:1000){
simulated_call_payoffs[i] <- exp(-r_T[i]*T)*max(ST[i]-K,0)
}
Call_price_asset_1 <- mean(simulated_call_payoffs)


#now find the evolution of the asset
asset2 = matrix(NA, nrow=(250*.5 +1), ncol=1000) 
asset2 = as.data.frame(asset2)

asset2[1,] = 10 #S2 like we hvae in line 314

for(j in 1:1000){  #for each simulated path
  for(i in 2:(250*.5 +1)){
    asset2[i,j] = asset2[i-1,j] + rate_paths3[i-1,j]*asset2[i-1,j]*dt + sigma21*(asset1[i-1,j]-asset2[i-1,j])*sqrt(dt)*rnorm(1) 
  }
}

#to find max S1 and S2 along each of the 10,000 simulated paths
max1 = matrix(NA, nrow=1, ncol=1000) 
max2 = matrix(NA, nrow=1, ncol=1000) 
max1 = as.data.frame(max1)
max2 = as.data.frame(max2)

for(i in 1:1000){
  max1[1,i] = max(asset1[,i])   
  max2[1,i] = max(asset2[,i])   
}

overall_max = matrix(NA, nrow=1, ncol=1000) 
overall_max = as.data.frame(overall_max)
for(i in 1:1000){
  overall_max[1,i] = max(c(max1[1,i],max2[1,i]))   
}

#now price the exotic option:
K=10
T=.5


simulated_call_payoffs = c() #these are the payoffs for each simulation
for(i in 1:1000){
  simulated_call_payoffs[i] <- exp(-r_T[i]*T)*max(overall_max[1,i]-K,0)
}
Call_price_exotic_option <- mean(simulated_call_payoffs)

A European Call option on asset 1 costs $0.7014

The exotic option costs $1.281464

Hedging, Large Price Movements, and Transaction Costs

stock_price = 100
strike_price = 100
r = 0.2 #This is mu instead of r
dividend = 0
sigma = 0.3
Time = 30/365
N = 96*30
num_sim = 1

h = Time/N
u = exp((r-dividend)*h + sigma * (h ^ 0.5))
d = exp((r - dividend)*h - sigma * (h ^ 0.5))
prob = (exp((r - dividend) * h) - d) / (u - d)

#lnSt ~ N(lnSo + (r-div - sigma^2 / 2)T, sigma^2 T)
mu =  (r - dividend - ((sigma) ** 2 / 2)) * h
std = sigma * (h ** 0.5)

Time_0 = rep(log(stock_price),num_sim)

#Parallel computing function
sim_one_node = function(last_node){
  one_node = rnorm(num_sim,mu,std) 
  new_node = last_node + one_node 
  return(new_node)
}  

ms = as.data.frame(NA)

# Monte Carlo Simulation simulates all of the path and then take the exponential
last = Time_0
ms = as.data.frame(matrix(numeric(0), ncol = 1, nrow = num_sim))
ms[1] = Time_0

for (i in 1:N) {
# run 96*90 times
   out = sim_one_node(last)
   ms[i + 1] = out
   last = out
}

ms = exp(ms) #turn back from lnSt to St

delta <- function(S, K, r, T, sig, div, type){
   d1 = (log(S/K) + (r -div + sig^2/2)*T) / (sig*sqrt(T))
   return (exp(-div*T)*pnorm(d1))
}

#hedge ratio
Delta = c()
j = 1
for (T in seq((60/365), (30/365), -(30/365)/(30*96))){
  Delta[j] = delta(S = S[j], K = 100, r = 0.05, T = T, sig = 0.3, div = 0, type = "C")
  j = j + 1
}

#rf assets to start with
B = c(BS[1] - Delta[1] * S[1])

#replicating portfolio
V = c(BS[1])

set.seed(2)

stock_price = 100
strike_price = 100
r = 0.2 #This is mu instead of r
dividend = 0
sigma = 0.3
Time = 30/365
N = 96*30
num_sim = 1

h = Time/N
u = exp((r-dividend)*h + sigma * (h ^ 0.5))
d = exp((r - dividend)*h - sigma * (h ^ 0.5))
prob = (exp((r - dividend) * h) - d) / (u - d)

#lnSt ~ N(lnSo + (r-div - sigma^2 / 2)T, sigma^2 T)
mu =  (r - dividend - ((sigma) ** 2 / 2)) * h
std = sigma * (h ** 0.5)

Time_0 = rep(log(stock_price),num_sim)

#Parallel computing function
sim_one_node = function(last_node){
  one_node = rnorm(num_sim,mu,std) 
  new_node = last_node + one_node 
  return(new_node)
}  

ms = as.data.frame(NA)

# Monte Carlo Simulation simulates all of the path and then take the exponential
last = Time_0
ms = as.data.frame(matrix(numeric(0), ncol = 1, nrow = num_sim))
ms[1] = Time_0

for (i in 1:N/2) {
# run first half
   out = sim_one_node(last)
   ms[i + 1] = out
   last = out
}

ms[N/2 + 2] = last + log(0.9)  #Because ln(S * 0.9) = lnS + ln0.9
last = ms[N/2 + 2]

for (i in (N/2 + 2):N) {
# run second half
   out = sim_one_node(last)
   ms[i + 1] = out
   last = out 
}

ms = exp(ms) #turn back from lnSt to St

S = ms_plot$value
BS = c()
j = 1

for (T in seq((60/365), (30/365), -(30/365)/(30*96))){
  BS[j] = BlackScholes(S = S[j], K = 100, r = 0.05, T = T, sig = 0.3, div = 0, type = "C")
  j = j + 1
}

delta <- function(S, K, r, T, sig, div, type){
   d1 = (log(S/K) + (r -div + sig^2/2)*T) / (sig*sqrt(T))
   return (exp(-div*T)*pnorm(d1))
}

#hedge ratio
Delta = c()
j = 1
for (T in seq((60/365), (30/365), -(30/365)/(30*96))){
  Delta[j] = delta(S = S[j], K = 100, r = 0.05, T = T, sig = 0.3, div = 0, type = "C")
  j = j + 1
}

#rf assets to start with
B = c(BS[1] - Delta[1] * S[1])

#replicating portfolio
V = c(BS[1])

for (i in 2:length(S)){
  V[i] = Delta[i-1] * S[i] + B[i-1] * exp(0.05*(30/365)/(30*96))
  B[i] = V[i] - Delta[i] * S[i]
}

Of all the option prices’, replicating portfolio value and hedge ratio show a downward jump as well. Also, the replicating portfolio value shows larger differences from the Black-Scholes option price after the jump.

set.seed(4)

ms = as.data.frame(NA)

# Monte Carlo Simulation simulates all of the path and then take the exponential
last = Time_0
ms = as.data.frame(matrix(numeric(0), ncol = 1, nrow = num_sim))
ms[1] = Time_0

for (i in 1:N/2) {
# run first half
   out = sim_one_node(last)
   ms[i + 1] = out
   last = out
}

ms[N/2 + 2] = last + log(1.1)  #Because ln(S * 1.1) = lnS + ln1.1
last = ms[N/2 + 2]

for (i in (N/2 + 2):N) {
# run second half
   out = sim_one_node(last)
   ms[i + 1] = out
   last = out 
}

ms = exp(ms) #turn back from lnSt to St

S = ms_plot$value
BS = c()
j = 1

for (T in seq((60/365), (30/365), -(30/365)/(30*96))){
  BS[j] = BlackScholes(S = S[j], K = 100, r = 0.05, T = T, sig = 0.3, div = 0, type = "C")
  j = j + 1
}

delta <- function(S, K, r, T, sig, div, type){
   d1 = (log(S/K) + (r -div + sig^2/2)*T) / (sig*sqrt(T))
   return (exp(-div*T)*pnorm(d1))
}

#hedge ratio
Delta = c()
j = 1
for (T in seq((60/365), (30/365), -(30/365)/(30*96))){
  Delta[j] = delta(S = S[j], K = 100, r = 0.05, T = T, sig = 0.3, div = 0, type = "C")
  j = j + 1
}

#rf assets to start with
B = c(BS[1] - Delta[1] * S[1])

#replicating portfolio
V = c(BS[1])

for (i in 2:length(S)){
  V[i] = Delta[i-1] * S[i] + B[i-1] * exp(0.05*(30/365)/(30*96))
  B[i] = V[i] - Delta[i] * S[i]
}

We observe that all of option price, replicating portfolio value and hedge ratio show an upward jump as well. Also, the replicating portfolio value shows larger differences from the Black-Scholes option price after the jump.

ms = as.data.frame(NA)

# Monte Carlo Simulation simulates all of the path and then take the exponential
last = Time_0
ms = as.data.frame(matrix(numeric(0), ncol = 1, nrow = num_sim))
ms[1] = Time_0

for (i in 1:N) {
# run 96*90 times
   out = sim_one_node(last)
   ms[i + 1] = out
   last = out
}

ms = exp(ms) #turn back from lnSt to St

S = ms_plot$value
BS = c()
j = 1

for (T in seq((60/365), (30/365), -(30/365)/(30*96))){
  BS[j] = BlackScholes(S = S[j], K = 100, r = 0.05, T = T, sig = 0.3, div = 0, type = "C")
  j = j + 1
}

delta <- function(S, K, r, T, sig, div, type){
   d1 = (log(S/K) + (r -div + sig^2/2)*T) / (sig*sqrt(T))
   return (exp(-div*T)*pnorm(d1))
}

#hedge ratio
Delta = c()
j = 1
for (T in seq((60/365), (30/365), -(30/365)/(30*96))){
  Delta[j] = delta(S = S[j], K = 100, r = 0.05, T = T, sig = 0.3, div = 0, type = "C")
  j = j + 1
}

#rf assets to start with
B = c(BS[1] - Delta[1] * S[1])

#replicating portfolio
V = c(BS[1])

for (i in 2:length(S)){
  V[i] = Delta[i-1] * S[i] + B[i-1] * exp(0.05*(30/365)/(30*96))
  B[i] = (V[i] - abs(Delta[i] - Delta[i-1]) * S[i] * 0.002) - Delta[i] * S[i]
}

It can be observed that if transaction costs are imposed, the value of replicating portfolio will vary greatly from the Black-Scholes option price.

Heston Model

BlackScholes(S = 100, K = 100, r = 0.04, T = 0.5, sig = sqrt(0.02), div = 0, type = "C")

## [1] 5.016981

#X is strike, q is div, v0 is vt, vT is theta, k is keppa
callHestoncf(S = 100, X = 100, tau = 0.5, r = 0.04, q = 0, v0 = 0.01, vT = 0.02, rho = -0.5, k = 6, sigma = 0.3, implVol = F)

## [1] 4.683304

So the Black-Scholes price of the option is \(5.016981\), and the Heston price of the option is \(4.683304\)

BS = c()
Heston = c()
j = 1

for (S in seq(70, 130, 0.1)){
   BS[j] = BlackScholes(S = S, K = 100, r = 0.04, T = 0.5, sig = sqrt(0.02), div = 0, type = "C")
   Heston[j] = callHestoncf(S = S, X = 100, tau = 0.5, r = 0.04, q = 0, v0 = 0.01, vT = 0.02, rho = -0.5, k = 6, sigma = 0.3, implVol = F)
   j = j+1
}

diff = Heston - BS

plot(x = seq(70, 130, 0.1), y = diff, type = "l", col = "blue", main = "Difference in Heston and B-S prices", xlab = "S", ylab = "H - BS")

ImpVol = c()
j = 1
for (S in seq(70, 130, 0.1)) {
  ImpVol[j] = callHestoncf(S = S, X = 100, tau = 0.5, r = 0.04, q = 0, v0 = 0.01, vT = 0.02, rho = -0.5, k = 6, sigma = 0.3, implVol = T)[[2]]
  j = j + 1
}

ImpVol = c()
j = 1
for (S in seq(70, 130, 0.1)) {
  ImpVol[j] = callHestoncf(S = S, X = 100, tau = 0.5, r = 0.04, q = 0, v0 = 0.01, vT = 0.02, rho = 0.5, k = 6, sigma = 0.3, implVol = T)[[2]]
  j = j + 1
}

As the correlation between two Brownian motions changes to 0.5, the implied volatility also turns to decreasing with respect to the underlying price.