Option pricing in Vietnam

library(DT)
library(tidyverse)
library(htmlwidgets)
library(xts)
library(ggthemes)
library(plotly)

## set working directory to the folder containing data
setwd("/Users/tantoankhoa10/Downloads/forecast")
## Import csv file into R 
df <- read.csv(file="excel_vic.csv", head=T,stringsAsFactors = F)
df %>% head(n=100) %>% datatable()

df <- df[,c(2:7)] # only select interested columns
df %>% head(n=100) %>% datatable()

date <- as.Date(as.character(df$X.DTYYYYMMDD.), format="%Y%m%d")
# transform date column to a suitable formate
class(date)

df <- cbind(date,df)
df %>% head(n=100) %>% datatable()

df <- df[,-2] # eliminate the second column
df %>% head(n=100) %>% datatable()

df <- df[order(df$date),] # arrange data from oldest to newest
df %>% head(n=100) %>% datatable()

names(df) <- paste(c("date", "Open", "High", "Low", "Close", "Volume"))
df %>% head(n=100) %>% datatable()

df<-subset(df,df$date >= "2014-12-31" )
df_plot=ggplot(df, aes(x=date, y= Close, group=1))+geom_point(size=0.2)+geom_line(col="blue")+labs(x="Trading dates", y= "Prices", title="Price movement of VIC since 31/12/2014")
ggplotly(df_plot)

df <- xts(df[,2:6],order.by = df[,1])
X<-df[,4]
X <- as.numeric(na.omit(X))
plot(X,type="l",col="blue", main="VIC prices since 31/12/2014",xlab="Number of trading dates", ylab="Price")

library(yuima)
Delta <- 1/252
gBm <- setModel(drift="mu*x", diffusion="sigma*x")
mod <- setYuima(model=gBm, data=setData(X, delta=Delta))
fit <- qmle(mod, start=list(mu=0.01, sigma=0.01),
            lower=list(mu=0, sigma=0),
            upper=list(mu=100, sigma=100))
summary(fit)

## The stock price is governed by the equation dS_t= 0.302 S_tdt+0.272 S_tdW_t

GBM=function(T,N) {
# T=1 expiry time
# N=100 number of simulation points
h=T/N # the timestep of the simulation
X=rep(0, (N+1)) # create space to store value of Brownian motion
X[1]=91 ## start at 0
for(i in 1:N) { X[i+1]=X[i]*exp((0.302-0.272^2/2)*h+0.272*sqrt(h)*rnorm(1))}
return(X)
}


GBMSamplepaths <- function(T,N,nt)
{h=T/N ## time step
#nt: number of samples
t=seq(0,T, by=h) # produce a sequence of time points
X=matrix(rep(0, length(t) *nt), nrow=nt)
X_median=matrix(rep(0, length(t)), nrow=1)
## each row of the matrix stores one sample
# #return(X)
for (i in 1:nt) {X[i,]= GBM(T=T,N=N)}
for (j in 1:length(t)){X_median[j]= median(X[,j])}

# # ##Plot
ymax=max(X); ymin=min(X) #bounds for simulated prices
plot(t,X[1,], t="l" ,ylim=c(ymin, ymax), col=1,ylab="Price P(t)",xlab="time t")
for(i in 2:nt){lines(t,X[i,], t="l",ylim=c(ymin, ymax), col="gray")}
lines(t, X_median, col="red")
return(X_median)}
GBMSamplepaths(T=0.5,N=252/2,nt=252*5)

## Suppose we would like to compute the price of a European call option written on VIC price
## With current price S0=91, strike price K =95, expiry T=0.5, r= 0.07, volatility sigma=0.272

BS_call <-  function (S0,K,T,r,sigma){
d1 <- (log(S0/K)+(r+0.5*sigma^2)*T)/(sigma*sqrt(T))
d2 <- d1- sigma*sqrt(T)
opt.val <- S0*pnorm(d1)-K*exp(-r*T)*pnorm(d2)
return(opt.val)}

BS_call (S0=91,K=95,T=0.5,r=0.07,sigma=0.272)

European_call_binomial=function( T, t,r, X,N,S0,sigma){
delta_t=(T-t)/N
u = exp(sigma*sqrt(delta_t))
d = exp(-sigma*sqrt(delta_t))
p = (exp(r * delta_t) - d)/(u - d) ### probability of up move
Df = exp(-r * (T-t)) ### discounting factor
Ce=0;
for (i in 1:(N+1)){Ce =Ce+choose(N,i-1)*p^{i-1}*(1-p)^{N+1-i}*max(S0*u^{i-1}*d^{N-i+1}-X,0)}
return(Df*Ce)}
European_call_binomial(T=0.5, t=0, r=0.07, X=95, N=1000, S0=91, sigma=0.272)