#Not parallel
MCPrice <- function(x = 1, tv = 0, Tv = 1, rv = 1, sigma = 1,
                    Mv = 1000, fn){
    h <- function(m) {
        u <- rnorm(m/2)
        tmp <- c(x * exp((rv - 0.5 * sigma^2) * (Tv - tv) + sigma *
          sqrt(Tv - tv) * u), x * exp((rv - 0.5 * sigma^2) * (Tv -
          tv) + sigma * sqrt(Tv - tv) * (-u)))
        mean(sapply(tmp, function(xx) fn(xx)))
    }
    p <- h(Mv)
    p * exp(-rv * (Tv - tv))
}

S0<-100
K<-110
r<-0.05
Tval<-1/4
sigma<-0.25

f<-function(x){return(max(0,x-K))}

K<-110
set.seed(123)
M<-50000
MCPrice(x=S0,tv=0,Tv=Tval,rv=r,sigma,Mv=M,fn=f)

###########################################################################
f<-function(x){return(max(0,x-110))}

set.seed(123)
M<-5e6
system.time(Price0<-MCPrice(x=S0,tv=0,Tv=Tval,rv=r,sigma,Mv=M,fn=f))
print(Price0)

MCPrice.par <- function(x = 1, tv = 0, Tv = 1, rv = 1, sigma = 1,
                    Mv = 1000, fn){
        require(foreach)
        h <- function(m) {
          u <- rnorm(m/2)
          tmp <- c(x * exp((rv - 0.5 * sigma^2) * (Tv - tv) + sigma *
            sqrt(Tv - tv) * u), x * exp((rv - 0.5 * sigma^2) * (Tv -
            tv) + sigma * sqrt(Tv - tv) * (-u)))
            mean(sapply(tmp, function(xx) fn(xx)))
        }
        nodes <- getDoParWorkers()
        print(paste('Nodes being used: ',nodes))
        p <- foreach(m = rep(Mv/nodes, nodes),.combine = "c") %dopar% h(m)
        p <- mean(p)
        p * exp(-rv * (Tv - tv))
}
set.seed(123)
M<-5e6
system.time(Price<-MCPrice.par(x=S0,tv=0,Tv=Tval,rv=r,sigma,Mv=M,fn=f))
print(Price)
###########################################################################
#install.packages('snowfall')
#install.packages('doSNOW')

require(snowfall)
sfInit(parallel = TRUE, cpus = 6)
cl <- sfGetCluster()
clusterSetupRNG(cl, seed = rep(123, 2))
require(foreach)
require(doSNOW)
registerDoSNOW(cl)

getDoParWorkers()
M <- 1e7
system.time(MCPrice.par(x = S0, tv = 0, Tv = T, rv = r, sigma, Mv = M, fn = f))
system.time(MCPrice(x = S0, tv = 0, Tv = T, rv = r, sigma, Mv = M, fn = f))


set.seed(123)
m <- c(10, 50, 100, 150, 200, 250, 500, 1000, 10000)
p1 <- NULL
err <- NULL
nM <- length(m)
repl <- 100
mat <- matrix(, repl, nM)
for (k in 1:nM) {
  tmp <- numeric(repl)
  for (i in 1:repl) tmp[i] <- MCPrice(x = S0, tv = 0, Tv = T, rv = r, sigma, Mv = m[k], fn = f)
  mat[, k] <- tmp
  p1 <- c(p1, mean(tmp))
  err <- c(err, sd(tmp))
}
colnames(mat) <- m

#install.packages('fOptions')
require(fOptions)
p0 <- GBSOption(TypeFlag = "c", S = S0, X = K, Time = T, r = r, b = r, sigma = sigma)@price
minP <- min(p1 - err)
maxP <- max(p1 + err)
plot(m, p1, type = "n", ylim = c(minP, maxP), axes = F, ylab = "MC price", xlab = "MC replications")
lines(m, p1 + err, col = "blue")
lines(m, p1 - err, col = "blue")
axis(2, p0, "B&S price")
axis(1, m)
boxplot(mat, add = TRUE, at = m, boxwex = 15, col = "orange", axes = F)
points(m, p1, col = "blue", lwd = 3, lty = 3)
abline(h = p0, lty = 2, col = "red", lwd = 3)
sfStop()
registerDoSEQ()