moments.R

first_moment<-function(s_0,r,TotalTime)
{
  (s_0*(exp(r*TotalTime)-1))/(r*TotalTime)
}

second_moment<-function(s_0,r,s,TotalTime)
{
  part1 = (2*s_0^2)/(r*(r+s^2)*(2*r+s^2)*TotalTime^2)
  
  part2 = r*exp((2*r+s^2)*TotalTime)-(2*r+s^2)*exp(r*TotalTime)+r+s^2
  
  part1*part2
}

simulateGBM.R

simulateGBM=function(s_0,drift,s,Time,dt)
{
  time_points= seq(from=0,to=Time,by=dt)
  stock=seq(from=0,to=0,length=length(time_points))
  stock_two=seq(from=0,to=0,length=length(time_points))
  
  #for antithetic variables
  stock[1]= log(s_0)
  stock_two[1]= log(s_0)
  
  for(i in 2:length(time_points))
  {
    z=rnorm(1)
    stock[i]=stock[i-1]+(drift-0.5*s^2)*dt+s*sqrt(dt)*z
    stock_two[i]=stock_two[i-1]+(drift-0.5*s^2)*dt+s*sqrt(dt)*(rnorm(1))*(-z)
  }
  
  return(cbind(exp(stock),exp(stock_two)))
}

simulationStudy.R

# Parameters
s_0 = 100
s = 0.3
TotalTime = 1
delta = 0.0001
r = 0.09
time_points = seq(from=0,to=TotalTime,by=delta)

## Set up for the monte carlo simulation
total_paths = c(50,100,500,1000,5000,10000,20000,40000,100000)
strike_vector = c(0,50,90,100,110,150,200)
result = matrix(data = 0, nrow = length(total_paths), ncol = length(strike_vector))
std_dev = matrix(data = 0, nrow = length(total_paths), ncol = length(strike_vector))

lowr_bnd = matrix(data = 0, nrow = length(total_paths), ncol = length(strike_vector))
uppr_bnd = matrix(data = 0, nrow = length(total_paths), ncol = length(strike_vector))

for(path_number in 1:length(total_paths))
{
  no_paths= total_paths[path_number]
  paths= matrix(data = 0, nrow = length(time_points), ncol = no_paths)
  paths_two= matrix(data = 0, nrow = length(time_points), ncol = no_paths)
  
  for(i in 1:no_paths)
  {
    paths[,i]=simulateGBM(s_0,r,s,TotalTime,delta)[,1]
    paths_two[,i]=simulateGBM(s_0,r,s,TotalTime,delta)[,2]
    
  }
  
  for(m in 1:length(strike_vector))
  {
    Strike = strike_vector[m]
    one = mean(pmax(colMeans(paths)-Strike,0))*exp(-r*TotalTime)
    two = mean(pmax(colMeans(paths_two)-Strike,0))*exp(-r*TotalTime)
    result[path_number,m] = (one+two)/2
    
    sum = (pmax(colMeans(paths)-Strike,0)+pmax(colMeans(paths_two)-Strike,0)) *exp(-r*TotalTime)/2
    std_dev[path_number,m] = sd(sum)
    
    lowr_bnd[path_number,m]= result[path_number,m]-1.96*std_dev[path_number,m] /sqrt(no_paths)
    uppr_bnd[path_number,m]= result[path_number,m]+1.96*std_dev[path_number,m] /sqrt(no_paths)
  }
  
  print(paste("Path Index  :", path_number, " of ", length(total_paths)))
  print(paste("Number of paths :", no_paths))
}


### Example convergence profile of Monte Carlo
x_m = result[1:8,2]
lower = lowr_bnd[1:8,2]
upper = uppr_bnd[1:8,2]

xvals= total_paths[1:8]#seq(from=100,to=100*total_paths,by=100)
plot(xvals,x_m,type="l",ylim=c(47,52),xlab="Number of paths",ylab= "Asian Option value",
     main="Convergence Profile Of Arithmetic Asian Option")
lines(xvals,lower,type="l",col="red")
lines(xvals,upper,type="l",col="red")
legend("topright", legend=c("Asian Option Value","95% confidence interval"),
       col=c("black","red"), lty=1:2, cex=0.8)

###### Log normal approximation
m1 = first_moment(s_0,r,TotalTime)
m2 = second_moment(s_0,r,s,TotalTime)

mu_of_average = log(m1^2/sqrt(m2))
sig_of_average = sqrt(log(m2/m1^2))

ln_result= seq(from=0,to=0,length=length(strike_vector))

for(i in 1:length(ln_result))
{
  ln_result[i]=price_option(mu_of_average,r,TotalTime,sig_of_average,strike_vector[i])
}

### Percentage difference between the best Monte Carlo and Log Normal approximation
((result[8,]-ln_result)/result[8,])*100