#=====================#
# DISCLAIMER
# This code is provided "as is" and "with all faults" for educational and research use only.
# The author makes no representations or warranties of any kind concerning its correctness or
# suitability for any particular purposes. Numbers and graphs in published papers were
# produced in Excel, not this code. If you find errors, please let me know. If you find
# errors, please let me know.
#=====================#
# Setting parameters
S <- 1.0 #stock price at time t
K <- 0.8 #strike price
b <- 0.5 #barrier - for b = K case, set b just slightly less than K, to avoid divisions by zero problems.
tau <- 25 #time to maturity T - t (in years)
r <- 0.015 #risk-free annual interest rate (convenient to set to zero)
q <- 0.01 #deferment rate (yield) (needs slightly different from r, to avoid division-by-zero problems in theta)
sigma <- 0.13 #annual volatility of the stock price (standard deviation)
set.seed(2411) #set the seed
N <- 6300
#number of sub intervals, 300 steps for 25 years (monthly). 6300 steps for daily.
dt <- tau/N #length of each time sub interval
nSim <-10 #number of simulations (paths)
#Check validity of inputs
stopifnot(b <= min(S,K), r!=q, K!=b)
#analytic prices
# z's as in the paper
z1 <- (log(S/K) + (r - q + 0.5*sigma^2)*tau)/(sigma*sqrt(tau))
z2 <- (log(b^2/(K*S)) + (r - q + 0.5*sigma^2)*tau)/(sigma*sqrt(tau))
z3 <- (log(S/b) + (r - q + 0.5*sigma^2)*tau)/(sigma*sqrt(tau))
z4 <- (log(b/S) + (r - q + 0.5*sigma^2)*tau)/(sigma*sqrt(tau))
theta <- 2*(r-q)/sigma^2
BS_call <- S*exp(-q*tau)*pnorm(z1) - K*exp(-r*(tau))*pnorm(z1-sigma*sqrt(tau)) #BS value call
BS_put <- -S*exp(-q*tau)*pnorm(-z1) + K*exp(-r*tau)*pnorm(-z1+sigma*sqrt(tau)) #BS value put
Analytic_barrier_call <- S*exp(-q*tau)*pnorm(z1) - K*exp(-r*(tau))*pnorm(z1-sigma*sqrt(tau))+
1/theta*(S*exp(-q*tau)*(b/S)^(1+theta)*pnorm(z2) - K*exp(-r*(tau))*(K/b)^(theta-1)*pnorm(z2-theta*sigma*sqrt(tau)))
Analytic_barrier_put <- K*exp(-r*tau)*pnorm(-z1+sigma*sqrt(tau)) - S*exp(-q*tau)*pnorm(-z1)-
b*exp(-r*tau)*pnorm(-z3+sigma*sqrt(tau)) + S*exp(-q*tau)*pnorm(-z3)+
1/theta*(b*exp(-r*tau)*pnorm(-z3+sigma*sqrt(tau)) - S*exp(-q*tau)*(b/S)^(1+theta)*(pnorm(z4)-pnorm(z2)) -
K*exp(-r*tau)*(K/b)^(theta-1)*pnorm(z2-theta*sigma*sqrt(tau)))
# Compute the Monte Carlo prices
dt <- tau/N #length of each time sub interval
Z <- matrix(rnorm(nSim*N, mean=0, sd=1),nrow = nSim, ncol = N) #standard normal sample of N elements
dW <- Z*sqrt(dt) #Brownian motion increments (nSim simulations) x (N increments)
W <- matrix(numeric(nSim*(N+1)), nrow = nSim, ncol = (N+1))
Y_T <- numeric(nSim)
Option_payoff_written_put <- numeric(nSim)
Final_forward_hedge_position <- numeric(nSim)
Final_forward_profit <- numeric(nSim)
Forward_hedging_error <- numeric(nSim)
BS_final_cash <-numeric(nSim)
BS_final_stock <-numeric(nSim)
BS_hedging_error <- numeric(nSim)
Thomas_final_cash <-numeric(nSim)
Thomas_final_stock <-numeric(nSim)
Thomas_hedging_error <- numeric(nSim)
for(k in 1:nSim){
X <- numeric(N)
B <- numeric(N)
B_running_max <- numeric(N)
Y <- numeric(N)
X[1] <- S
Y[1] <- S
B_running_max[1] = b/X[1]
#Forward hedge
Forward_hedge_delta <- 1 # hedge is static!
Forward_hedge_position <- numeric(N)
#Black-Scholes replication of barrier put
# NB NB In this program we are using Black-Scholes to hedge the BARRIER put - NOT the ordinary Black-Scholes put. This always works,
# because Black-Scholes works for ANY random GBM path, including one which happens to always turn upwards at say 0.5 (but without
# a barrier in place at 0.5). The RGBM path (ie with an actual barrier at 0.5) is observationally indistinguihsable from this. So
# Black-Scholes replication always works for the barrier put. (But it's unnecessarily expensive; my replication scheme is cheaper.)
BS_delta_written_put <- numeric(N)
BS_stock_position_before_trade <- numeric(N)
BS_stock_position_after_trade <- numeric(N)
BS_trade_size <-numeric(N)
BS_cash_account <-numeric(N)
BS_delta_written_put[1] <- -exp(-q*tau)*(pnorm(z1)-1) #positive number
BS_stock_position_after_trade[1] <- BS_delta_written_put[1] #positive number, we're buying shares to synthesise a written put
BS_cash_account[1] <- -S * BS_delta_written_put[1] - BS_put #negative number, borrowing the money to buy the asset (no interest
#cost if r = q). Also borrowing more, because a written put involves owning the asset at maturity,but with a loss from the higher strike.
#Perhaps best understood by considering replication of BOUGHT put, You take position S * Delta_Put in the asset (a negative number,
# i.e. a short). You invest V_Put - S * Delta_Put in the bond, where Delta_Put is a negative number, so the investment is greater
# than V_Put. This enables you to replicate payoffs of a bought put. So wouldn't pay more than V_BS_Put for a bought put, because
# if you had that amount of cash, you could replicate put payoffs. For WRITTEN put, flip the sign on both trades.
BS_trade_size[1] <- 0
#Thomas put replication
Thomas_delta_written_put <- numeric(N)
Thomas_stock_position_after_trade <- numeric(N)
Thomas_trade_size <-numeric(N)
Thomas_cash_account <-numeric(N)
Thomas_z1 <-numeric(N)
Thomas_z2 <-numeric(N)
Thomas_z3 <-numeric(N)
Thomas_z4 <-numeric(N)
Thomas_delta_written_put[1] <- -exp(-q*tau)*(pnorm(z1) - pnorm(z3) + (b/S)^(1+theta)*(pnorm(z4)-pnorm(z2)))
#using z1,...,z4 from top of program...OK at time 1.
Thomas_stock_position_after_trade[1] <- Thomas_delta_written_put[1] # same as above
Thomas_cash_account[1] <- -S * Thomas_delta_written_put[1] - Analytic_barrier_put #negative number, borrowing to buy the asset.
Thomas_trade_size[1] <- 0
# Then do the loop with hedging at every step. This is slow, but given the recursive nature of the definition of Y_t
# and the need to hedge at every step, . I can't see how to speed it up.
for(j in 2:N){
X[j] <- X[j-1]*exp((r - q -0.5*sigma^2)*dt + sigma*dW[k,j]) # CHANGE r-q here to a positive drift, to check hedging works.
B[j] <- b/X[j]
B_running_max[j] <- ifelse(B[j] > B_running_max[j-1], B[j], B_running_max[j-1])
#BMax[j] <- max(B[1:j])
# Max of b/X so far (separate line so I could print and check it works)
#Done as a running maxinum to avoid calculating MAX(B[1;j]) in every interation of the loop
Y[j] <- X[j]*max(1,B_running_max[j])
#First forward hedging
Forward_hedge_position[j] <- Y[j] * Forward_hedge_delta
# trivial hedge, which gets ALL the benefit of the interventions.
# The hedge neutralises the drift, but NOT the interventions - correct for a forward!
#Then for BS hedging
BS_delta_written_put[j] <- -exp(-q*tau*(1-j/N))*(pnorm((log(Y[j]/K) +(r- q +0.5*sigma^2)*(tau*(1-j/N)))/(sigma*sqrt(tau*(1-j/N))))-1)
BS_trade_size[j] = BS_delta_written_put[j] - BS_delta_written_put[j-1] #if delta as gone down, need to sell some
BS_stock_position_after_trade[j] <- BS_stock_position_after_trade[j-1] + BS_trade_size[j]
# trade restores position to fractional size delta, based on new asset price
#E.g. Selling some if price has gone up, which implies delta of a written put has fallen
# stock numeraire here, no need to keep track of actual stock price until end, except in valuing the trade in the cash account
BS_cash_account[j] <- BS_cash_account[j-1] * exp(r*dt) + BS_stock_position_after_trade[j-1] *Y[j] * q *dt - BS_trade_size[j] * Y[j]
#Pounds numeraire here, so need stock price
#Then for Thomas hedging
Thomas_z1[j] <-(log(Y[j]/K) +(r - q +0.5*sigma^2)*(tau*(1-j/N)))/(sigma*sqrt(tau*(1-j/N)))
Thomas_z2[j] <-(log(b^2/(K*Y[j])) +(r - q +0.5*sigma^2)*(tau*(1-j/N)))/(sigma*sqrt(tau*(1-j/N)))
Thomas_z3[j] <-(log(Y[j]/b) +(r - q + 0.5*sigma^2)*(tau*(1-j/N)))/(sigma*sqrt(tau*(1-j/N)))
Thomas_z4[j] <-(log(b/Y[j]) +(r - q + 0.5*sigma^2)*(tau*(1-j/N)))/(sigma*sqrt(tau*(1-j/N)))
if (j!=N) {
Thomas_delta_written_put[j] <- -exp(-q*tau*(1-j/N))*(pnorm(Thomas_z1[j]) - pnorm(Thomas_z3[j]) +
(b/Y[j])^(1+theta) * (pnorm(Thomas_z4[j]) - pnorm(Thomas_z2[j])))
} else {
Thomas_delta_written_put[j] <- ifelse(K > Y[j], 1, 0)
}
# Setting the final delta manually avoids the deep bug of z2, z3, z4 evaluating as 0/0 if Y_T[k]] = b
# and time remaining = 0. Could alternatively fix by setting N to N+1 in the z's, which would
# probably run quicker, but give very slight discrepancy against others' evaluations.
# Alternatively, excluding the NaN cases where Y_T[k] = b from summary statistics would be OK, since
# the option pays off zero in those cases anyway.
Thomas_trade_size[j] <- Thomas_delta_written_put[j] - Thomas_delta_written_put[j-1]
Thomas_stock_position_after_trade[j] <- Thomas_stock_position_after_trade[j-1] + Thomas_trade_size[j]
Thomas_cash_account[j] <- Thomas_cash_account[j-1] * exp(r*dt) +
Thomas_stock_position_after_trade[j-1] *Y[j] * q *dt - Thomas_trade_size[j] * Y[j]
}
Y_T[k] <- Y[j]
Option_payoff_written_put[k] <- - max(K - Y_T[k],0) # minus sign as written
#First for forward hedge
Final_forward_hedge_position[k] <- Forward_hedge_position[j]
Final_forward_profit[k] <- Y_T[k] - Final_forward_hedge_position[k] # alwqays zero!
Forward_hedging_error <- Final_forward_profit #always zero!
# Then for BS hedging
BS_final_cash[k] <- BS_cash_account[j]
# The replication recipe gives 0 stock and 0 cash where put not due to be exercised, or 1 stock and
# -K cash where it is; in the latter case, you have a net wealth (of Y_T - K), which is negative, correct
# for the writer of an exercised put. So the recipe has always generated the payoffs of a written put.
BS_final_stock[k] <- BS_stock_position_after_trade[j]
#Final stock will always be correct, either 1 or 0, forced by calculation of the final delta either 1 or 0,
BS_hedging_error[k] <- ((BS_final_cash[k]+ BS_final_stock[k]*Y_T[k]) - Option_payoff_written_put[k])
# ("proceeds of replication scheme" - "option payoff")
#Then for Thomas hedging
Thomas_final_cash[k] <- Thomas_cash_account[j]
Thomas_final_stock[k] <- Thomas_stock_position_after_trade[j]
Thomas_hedging_error[k] <- ((Thomas_final_cash[k] + Thomas_final_stock[k]*Y_T[k]) - Option_payoff_written_put[k])
}
payoff_expiry_call <-pmax(Y_T-K,0) #pmax preserve the dimension of the matrix, so apply the max function to each element
expected_payoff_call <- mean(payoff_expiry_call)
Monte_Carlo_call_price <- exp(-r*(tau))*expected_payoff_call
payoff_expiry_put <-pmax(K-Y_T,0) #pmax preserve the dimension of the matrix, so apply the max function to each element
expected_payoff_put <- mean(payoff_expiry_put)
Monte_Carlo_put_price <- exp(-r*(tau))*expected_payoff_put
#First for forward
Mean_forward_hedging_error <- mean(Forward_hedging_error)
# Forward hedging error is always zero - I include this only to show that the usual forward hedging recipe works
# perfectly. (People often think the barrier will change the forward price. Of course if you introduce a barrier
# where there was none before, that will raise both spot and forward **in £ numeraire**; but the spot-to-forward
# arbitrage is unaffected.
#Then means for BS hedging
Mean_BS_final_cash <-mean(BS_final_cash)
Max_BS_final_cash <- max(BS_final_cash)
Min_BS_final_cash <- min(BS_final_cash)
Mean_BS_hedging_error <- mean(BS_hedging_error)
#Then means for Thomas hedging
Thomas_mean_final_cash <-mean(Thomas_final_cash)
Thomas_max_final_cash <- max(Thomas_final_cash)
Thomas_min_final_cash <- min(Thomas_final_cash)
Thomas_mean_hedging_error <- mean(Thomas_hedging_error)
Ratio_Call <- Monte_Carlo_call_price/Analytic_barrier_call
Ratio_Put <- Monte_Carlo_put_price/Analytic_barrier_put
#mybins <-seq(-0.04, 0.04, 0.001) This can be used to make both histograms look same, with narrow bins.
hist(Forward_hedging_error)# always zero!
hist(BS_hedging_error)
hist(Thomas_hedging_error)
cat("Analytic B-Scholes Call :",round(BS_call, digits=5), " Analytic B-Scholes Put :",round(BS_put, digits=5))
cat("\nAnalytic Barrier Call :",round(Analytic_barrier_call, digits=5), " Analytic Barrier Put :",round(Analytic_barrier_put, digits=5))
cat("\nMonte-Carlo Barrier Call:",round(Monte_Carlo_call_price, digits=5), " Monte-Carlo Barrier Put :",round(Monte_Carlo_put_price, digits=5))
cat("\n MC/Analytic for Call:",round(Ratio_Call, digits=5), " MC/Analytic for Put:",round(Ratio_Put, digits=5))
cat("\nMean B-Scholes hedging error (for BARRIER put):",round(Mean_BS_hedging_error,digits=5), "Thomas ditto: ",round(Thomas_mean_hedging_error,digits=5))
#=====================#
#This program is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#GNU General Public License for more details.
#You should have received a copy of the GNU General Public License
#along with this program. If not, see .