MonteCarlo Parisian Option

require(fOptions)

sobolInnovations = function(mcSteps, pathLength, init, ...) {
  innovations = rnorm.sobol(mcSteps, pathLength, init, ...)
  innovations
}

wienerPath = function(eps) {
  path = (b-sigma*sigma/2)*delta.t + sigma*sqrt(delta.t)*eps
  path
}

ParisianPayoff = function(path) {
  ST = S*exp(sum(path))

  path = S*exp(cumsum(path))
  path = (path <= H) + 0
  path = (rle(path)$values[rle(path)$lengths >= k] == 1)

  if (sum(path) > 0) {payoff = 0}
  else {
    if (TypeFlag == "c") payoff = exp(-r*Time)*max(ST-X, 0)
    if (TypeFlag == "p") payoff = exp(-r*Time)*max(0, X-ST)
  }
  payoff
}

MonteCarloOption2 = function (delta.t, pathLength, k, mcSteps, mcLoops, init = TRUE,
                              innovations.gen, path.gen, payoff.calc, antithetic = TRUE,
                              standardization = FALSE, trace = TRUE, ...)
{
  delta.t <<- delta.t
  k <<- k
  if (trace)
    cat("\nMonte Carlo Simulation Path:\n\n")
  iteration = rep(0, length = mcLoops)
  cat("\nLoop:\t", "No\t")
  for (i in 1:mcLoops) {
    if (i > 1)
      init = FALSE
    eps = innovations.gen(mcSteps, pathLength, init = init, ...)
    if (antithetic)
      eps = rbind(eps, -eps)
    if (standardization)
      eps = (eps - mean(eps))/sqrt(var(as.vector(eps)))
    path = t(path.gen(eps))
    payoff = NULL

    # Wrong in MonteCarloOption: for (j in 1:dim(path)[1])
    for (j in 1:dim(path)[2]) payoff = c(payoff, payoff.calc(path[,j]))
    iteration[i] = mean(payoff)
    if (trace)
      cat("\nLoop:\t", i, "\t:", iteration[i], sum(iteration)/i)
  }
  if (trace)
    cat("\n")
  iteration
}

TypeFlag <<- "c"; S <<- 47; X <<- 50; H <<- 45
Time <<- 0.5; sigma <<- 0.4; r <<- 0.1; b <<- 0.1

mc = MonteCarloOption2(delta.t = 1/2000, pathLength = floor(Time/delta.t), k = 0.1*floor(Time/delta.t), mcSteps = 50,
                      mcLoops = 3000, init = TRUE, innovations.gen = sobolInnovations,
                      path.gen = wienerPath, payoff.calc = ParisianPayoff,
                      antithetic = TRUE, standardization = FALSE, trace = TRUE,
                      scrambling = 2, seed = 4711)

par(mfrow = c(1, 1))
mcPrice = cumsum(mc)/(1:length(mc))
plot(mcPrice, type = "l", main = "Parisian Option", xlab = "Monte Carlo Loops", ylab = "Option Price")