# R-version of Pataki code

v0 <- 0.04;
S0 <- 1;
T <- 1;
vbar <- 0.04;
lambda <- 1.15;
eta <- 0.39;
rho <- -0.64;


HestonMC <- function(v0, S0, T, vbar, lambda, eta, rho, N, n, AK, vneg)
{
# v0: Initial volatility
# S0: Initial stock price
# T: Time frame for model
# vbar: Long term mean volatility
# lambda:
# eta:
# rho: Correlation between two B.M.’s
# N: Number of Monte Carlo paths
# n: Number of time steps per M.C. path
# AK: Array of strikes to evaluate call values for
# vneg: What to do if volatility goes negative:
# 0: absorbing (0 -> v)
# 1: reflecting (-v -> v)
# 2: Roger Lord
# 3: Milstein + Roger Lord

  dt = T/n;
  
  rho2m1 <- sqrt(1 - rho^2);
  negCount <- 0;
  # We use a vertical array, one element per M.C. path
  S <- rep(0,N); v <- rep(1,N)*v0;
  for (i in 1:n)
  {
    # Two sets of correlated normal random vars.
    # As alternative to antithetic variates, force mean<-0 and sd<-1
    W1 <- rnorm(N);
    W2 <- rnorm(N);
    W1 <- W1 - mean(W1); W1 <- W1/sd(W1);
    W2 <- W2 - mean(W2); W2 <- W2/sd(W2);
    # Now W1 and W2 are forced to have mean<-0 and sd<-1
    W2 <- rho*W1 + rho2m1*W2;
    # Stock and volatility SDE discretization
    if (vneg == 0){
      # Absorbing condition
      sqvdt <- sqrt(v*dt);
      S <- S  -v/2*dt + sqvdt * W1;
      v <- v + lambda*(vbar - v)* dt + eta * sqvdt * W2;
      negCount <- negCount + length(v[v < 0]);
      v[v < 0] <- 0;}
    if (vneg == 1){
      # Reflecting condition
      # Code to be filled in by you!
                    }
    if (vneg == 2){
      # Roger Lord
      # Code to be filled in by you!
      }
    if (vneg == 3){
      # Milstein + Roger Lord
      # Code to be filled in by you!
      }
  }
  negCount <- negCount / (n*N);
  M <- length(AK);
  AV <- numeric(M); AVdev <- numeric(M);
  BSV <- numeric(M); BSVH <- numeric(M); BSVL <- numeric(M);
  
  S <- S0*exp(S);
  # Evaluate mean call value for each path
  for (i in 1:M)
  {
    K <- AK[i];
    V <- (S>K)*(S - K); # Boundary condition for European call
    AV[i] <- mean(V);
    AVdev[i] <- 2 * sd(V) / sqrt(N);
  }
  return(data.frame(AV,AVdev));
}

AK <- c(0.7,1,1.4);
HestonMC(v0, S0, T, vbar, lambda, eta, rho, 500, 1000, AK, 0)