# R Monte Carlo Simulation Price Path Converging Volatility #number of simulation runs nSims = 1000 S = 100 r = 0 q = 0 volatility = 0.25 #drift term mu = r - q #every trading day increment dt = 1 / 365 vol_vec = vector() for (tDays in 1:365) { # standard normal distribution random number z=rnorm(tDays*nSims, mean=0, sd=1) # generate log-normal return matrix return_matrix=matrix(exp((mu-0.5*volatility^2) *dt+volatility*sqrt(dt) *z), ncol=nSims) # return value: price path matrix path=rbind(matrix(rep(S,nSims),ncol=nSims), S*apply(return_matrix,2,cumprod)) # calculate the volatility of path vol=mean(apply(path, 2, function(x){365*mean( (diff(log(x),n=1)-0.5*(volatility*dt)^2 )^2 ) } ) ) vol_vec=c(vol_vec, vol) } plot(vol_vec, type="l", col="blue") abline(h=volatility, col="black", pch=22, lty=2)