Skip to content

Instantly share code, notes, and snippets.

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

 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= log(s_0) stock_two= 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
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.