##################################################################

# Shengyu ZHENG

# ESSEC Business School

# June 21, 2022

# Illustration of Barrier Options Pricing with R

##################################################################

# 1. Installing packages and loading libraries

if (!require(derivmkts)) install.packages('derivmkts')
if (!require(ggplot2)) install.packages('ggplot2')
if (!require(dplyr)) install.packages('dplyr')

library("derivmkts") #Functions and R Code to Accompany Derivatives Markets
library("ggplot2") #For creating graphs
library("dplyr") #For data wrangling 


# 2. Characteristics of the barrier option and market data

STRIKE = 50 # Strike level
BARRIER = 40 # Barrier level
interval = 0 # Interval between two observation times (in days) (0 for continuous monitoring, 1 for daily monitoring)
T = 1 # Maturity (in years)

S0 = 50 # Spot price of the underlying asset
DIV = 0 # Dividend yield (annual rate)
VOL = 0.2 # Volatility of underlying (annual rate)
Rf = 0.1 # Risk-free rate (annual rate)


# 3. Using built-in formulas for pricing barrier options

# Constructing a function that can calculate all types of barrier options
BarrierOptionPrice <- function(InOut,UpDown,CallPut,s0,k,h,vol,div,tt,r,interval=0){

  InOut <- ifelse(tolower(InOut)=="in",1,-1)
  UpDown <- ifelse(tolower(UpDown)=="up",1,-1)
  CallPut <- ifelse(tolower(CallPut)=="call",1,-1) 
  
  if(interval != 0){
  
    Barrier_Modi <- h*exp(UpDown*0.5826*VOL*sqrt(interval/253))
    return(
      calldownin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(InOut+1)*(-UpDown+1)*(CallPut+1)/8+
      callupin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(InOut+1)*(UpDown+1)*(CallPut+1)/8+
      calldownout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(-InOut+1)*(-UpDown+1)*(CallPut+1)/8+
      callupout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(-InOut+1)*(UpDown+1)*(CallPut+1)/8+
      putdownin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(InOut+1)*(-UpDown+1)*(-CallPut+1)/8+
      putupin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(InOut+1)*(UpDown+1)*(-CallPut+1)/8+
      putdownout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(-InOut+1)*(-UpDown+1)*(-CallPut+1)/8+
      putupout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=Barrier_Modi)*(-InOut+1)*(UpDown+1)*(-CallPut+1)/8
    )

    }

else{
return(
      calldownin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(InOut+1)*(-UpDown+1)*(CallPut+1)/8+
      callupin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(InOut+1)*(UpDown+1)*(CallPut+1)/8+
      calldownout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(-InOut+1)*(-UpDown+1)*(CallPut+1)/8+
      callupout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(-InOut+1)*(UpDown+1)*(CallPut+1)/8+
      putdownin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(InOut+1)*(-UpDown+1)*(-CallPut+1)/8+
      putupin(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(InOut+1)*(UpDown+1)*(-CallPut+1)/8+
      putdownout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(-InOut+1)*(-UpDown+1)*(-CallPut+1)/8+
      putupout(s=s0,k=k,v=vol,r=r,tt=tt,d=div,H=h)*(-InOut+1)*(UpDown+1)*(-CallPut+1)/8
      )
  }
}


# 3.1 Continuous monitoring

## With the combined function
BarrierOptionPrice(InOut="out",UpDown="down",CallPut = "put",s0=S0,k=STRIKE,h=BARRIER, vol=VOL,div=DIV,tt=T,r=Rf,interval=0)

## With the built-in function
putdownout(s=S0,k=STRIKE,r=Rf,tt=T,d=DIV,H=BARRIER,v=VOL)

# 3.2 Discrete monitoring

## With the combined function
BarrierOptionPrice(InOut="out",UpDown="down",CallPut = "put",s0=S0,k=STRIKE,h=BARRIER, vol=VOL,div=DIV,tt=T,r=Rf,interval=1)

## With the built-in function
Barrier_Modi <- BARRIER*exp(-0.5826*VOL*sqrt(1/253))
putdownout(s=S0,k=STRIKE,r=Rf,tt=T,d=DIV,H=Barrier_Modi,v=VOL)


# 4. The simulation method (example for a down-and-out put)

m <- round(T*253) # Number of Knock observation days

s <- simprice(s0=S0,v=VOL,r=Rf,tt=T,d=DIV,trials=6000,periods = round(T*253),jump=FALSE, long=TRUE) #simulate trajectories 

s <- s %>%
  mutate (DnOut = if_else(price<BARRIER,1,0)) #mark the days where prices trigger the barrier event

LastDay = subset(s,s$period == round(T*253)) #a new data frame storing only the price for the last day


# Visualization 

Avg = mean(LastDay$price)
StdUp1 = Avg+1*sd(LastDay$price)
StdUp2 = Avg+2*sd(LastDay$price)
StdDn1 = Avg-1*sd(LastDay$price)
StdDn2 = Avg-2*sd(LastDay$price)

ggplot(s,aes(x=period,y=price,color=trial))+geom_line()+
  geom_hline(yintercept = Avg,color="black")+
  geom_hline(yintercept = StdUp1,color="green",linetype="dashed")+
  geom_hline(yintercept = StdUp2,color="green",linetype="dashed")+
  geom_hline(yintercept = StdDn1,color="red",linetype="dashed")+
  geom_hline(yintercept = StdDn2,color="red",linetype="dashed")+
  geom_hline(yintercept = STRIKE,color="blue",linetype="dashed")+
  geom_hline(yintercept = BARRIER,color="brown",linetype="dashed")+
  theme(legend.position = "none")


# Mark the trajectories where there has been at least one knock-out event

DnOut = subset(s,s$period %in% seq(round(T*253),1,by=-round(T*253/m)))
DnOut = aggregate(DnOut$DnOut,by=list(category=DnOut$trial),FUN=sum) 
LastDay$DnOut <- DnOut$x
LastDay <- LastDay %>% 
  mutate (Payoff = if_else(DnOut>0 | price>STRIKE,0,STRIKE-price)) #knocked-out and out-of-the-money options are marked with 0 payoff

# The present value of the mean of all possible payoffs should be the estimate of the premium (theoretical price) for the option

mean(LastDay$Payoff)*exp(-Rf*1)


# 5. Greeks

# Greeks of a single parameter setting

greeks(putdownout(s=S0,k=STRIKE,v=VOL,r=Rf,tt=T,d=DIV,H=BARRIER))


# Evolution of greeks of a changing parameter setting

STRIKES = seq(45,70,by=1)
p.greeks = greeks(putdownout(s=S0,k=STRIKES,v=VOL,r=Rf,tt=T,d=DIV,H=BARRIER),long=TRUE)
ggplot(p.greeks,aes(x=k,y=value,color=funcname))+
  geom_line()+
  labs(x="Strikes")+
  facet_wrap(~ greek, scales = "free")

##################################################################

# References:

# Broadie, Mark, Paul Glasserman and Steven Kou. (1997). A Continuity Correction for Discrete Barrier Option. Mathematical Finance, 7:325-349.

# Haug, Espen. (1997). The Complete Guide to Option Pricing. London/New York: McGraw-Hill.

# Jason R. (2020). Option Pricing with derivmkts | Option Greeks, Monte Carlo Simulation & Binomial Pricing, https://www.youtube.com/watch?v=aehwNnwn5VE&t=903s 

# Wang, Bing and Wang, Ling. (2011). Pricing Barrier Options using Monte Carlo Methods  

##################################################################