Instantly share code, notes, and snippets.

@WilsonMongwe WilsonMongwe/moments.R
Last active Dec 14, 2017

Embed
What would you like to do?
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=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)))
}
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment