#Parameters
S0 <- 100 #Initial stock price
K <- 90 #Strike price
r <- 0.09 #Risk-free rate
sigma <- 0.5 # Volatility
T <- 1 #Time to maturity
n <- 100000 #Number of simulations
N <- 120 # Number of periods
b <- 0  #Include S0

# Arithmetic Asian Option -------------------------------------------------

asian_option_mc <- function(S0, K, r, sigma, T, n, N, option_type="call", b) {
  dt <- T/N  #Time step
  
  #Simulating asset paths
  payoff <- numeric(n)
  for (i in 1:n) {
    S <- numeric(N+1)
    S[1] <- S0
    for (j in 2:(N+1)) {
      Z <- rnorm(1)
      S[j] <- S[j-1] * exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*Z)
    }
    
    #Discrete arithmetic average
    if (b == 0) {
      A <- mean(S)
    } else {
      A <- mean(S[-1])
    }
    
    # Compute payoff
    if (option_type == "call") {
      payoff[i] <- pmax(A - K, 0) * exp(-r * T)
    } else {
      payoff[i] <- pmax(K - A, 0) * exp(-r * T)
    }
  }
  
  option_price <- mean(payoff)
  error <- sqrt(1/(n-1)*sum((payoff - option_price)^2))
  
  return(list(price=option_price, standard_error=error/sqrt(n)))
}

result <- asian_option_mc(S0, K, r, sigma, T, n, N, "call", 0)
print(result)


# Antithetic Asian Options ------------------------------------------------

asian_option_An <- function(S0, K, r, sigma, T, n, N, option_type="call", b) {
  dt <- T/N  #Time step
  
  payoff <- numeric(n)
  for (i in 1:(n)) {
    S1 <- numeric(N+1)
    S2 <- numeric(N+1)
    S1[1] <- S0
    S2[1] <- S0
    for (j in 2:(N+1)) {
      Z <- rnorm(1)
      S1[j] <- S1[j-1] * exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*Z)
      S2[j] <- S2[j-1] * exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*(-Z))
    }
    
    #Discrete arithmetic average
    if (b == 0) {
      A1 <- mean(S1)
      A2 <- mean(S2)
    } else {
      A1 <- mean(S1[-1])
      A2 <- mean(S2[-1])
    }
    
    #Payoff
    if (option_type == "call") {
      payoff[i] <- exp(-r*T)* (pmax(A1 - K, 0) + pmax(A2 - K, 0))/2
    } else {
      payoff[i] <- exp(-r*T) * (pmax(K - A1, 0) + pmax(K - A2, 0))/2
    }
  }
  
  option_price <- mean(payoff)
  #error <- sqrt(1/(n-1)*sum((payoff - option_price)^2))
  
  return(list(price=option_price, standard_error=sd(payoff)/sqrt(n)))
}

#Parameters
S0 <- 100      # Initial stock price
K <- 105       # Strike price
r <- 0.05     # Risk-free rate
sigma <- 0.05  # Volatility
T <- 1        # Time to maturity
n <- 100000   # Number of simulations
N <- 40       # Number of fixing dates
b <- 0        # Include S0 = 0


resultAn <- asian_option_An(S0, K, r, sigma, T, n, N, "call",0)
print(resultAn)


# Control variates Asian Options ------------------------------------------

asian_option_cv <- function(S0, K, r, sigma, T, n, N, option_type="call", b) {
  dt <- T/N  # Time step
  
  # Simulating asset paths
  payoff <- numeric(n)
  control_variate <- numeric(n)
  for (i in 1:n) {
    S <- numeric(N+1)
    S[1] <- S0
    for (j in 2:(N+1)) {
      Z <- rnorm(1)
      S[j] <- S[j-1] * exp((r-0.5*sigma^2)*dt + sigma*sqrt(dt)*Z)
    }
    
    #Discrete arithmetic average
    if (b == 0) {
      A <- mean(S)
    } else {
      A <- mean(S[-1])
    }
    
    #Geometric average (control)
    if (b == 0) {
      #G <- prod(S)^(1/(N+1))
      G <- exp(1/(N+1)*sum(log(S)))
    } else {
      #G <- prod(S[-1])^(1/N)
      G <- exp(1/(N)*sum(log(S[-1])))
    }
    
    #Black-Scholes formula for geometric Asian option (Control Variate)
    sigma_adj2=((2*N+1)*(N+1)^b)/(6*N^b*(N+1-b)) *sigma^2
    mu = (N+1)/(2*(N+1-b))*(r-sigma^2/2)+sigma_adj2/2
    d1=(log(S0/K) + (mu + 0.5*sigma_adj2*T))/(sqrt(sigma_adj2*T))
    d2=d1 - sqrt(sigma_adj2*T)
    
    #Payoff
    if (option_type == "call") {
      geo_price=exp(-r*T)*( S0*exp(mu*T)*pnorm(d1)-K*pnorm(d2))
      payoff[i] <- pmax(A - K, 0) * exp(-r * T)
      control_variate[i] <- pmax(G - K, 0) * exp(-r * T)
      
    } else {
      geo_price=exp(-r*T)*(K*pnorm(-d2) - S0*exp(mu*T)*pnorm(-d1))
      payoff[i] <- pmax(K - A, 0) * exp(-r * T)
      control_variate[i] <- pmax(K - G, 0) * exp(-r * T)
    }
  }
  
  #Control variate optimal coefficient cv
  cv <- cov(control_variate, payoff) / var(control_variate)
  
  #Adjusted payoff
  adj_p <- payoff - cv * (control_variate - geo_price)
  option_price <- mean(adj_p)
  error <- sqrt(1/(n-1)*sum((adj_p - option_price)^2))
  
  return(list(price=option_price, standard_error=error/sqrt(n)))
}

#Parameters
S0 <- 50      # Initial stock price
K <- 60       # Strike price
r <- 0.1     # Risk-free rate
sigma <- 0.3  # Volatility
T <- 1        # Time to maturity
n <- 100000   # Number of simulations
N <- 100       # Number of fixing dates
b <- 0        # Include S0 = 0


result_cv <- asian_option_cv(S0, K, r, sigma, T, n, N, "put", 0)
print(result_cv)


#asian_option_mc(S0, K, r, sigma, T, n, N, "put")
#asian_option_An(S0, K, r, sigma, T, n, N, "put")
#asian_option_cv(S0, K, r, sigma, T, n, N, "put")











# Floating Strike Arithmetic Asian Option ------------------------------------

asian_floating_option_mc <- function(S0, lambda, r, sigma, T, n, m, option_type="call", b) {
  dt <- T/m  # Time step
  
  payoff <- numeric(n)
  for (i in 1:n) {
    S <- numeric(m + 1)
    S[1] <- S0
    for (j in 2:(m + 1)) {
      Z <- rnorm(1)
      S[j] <- S[j - 1] * exp((r - 0.5 * sigma^2) * dt + sigma * sqrt(dt) * Z)
    }
    
    # Discrete arithmetic average (with or without S0 based on b)
    if (b == 0) {
      A <- mean(S)
    } else {
      A <- mean(S[-1])
    }
    
    ST <- S[m + 1]  #for floating
    
    # Compute floating strike payoff
    if (option_type == "call") {
      payoff[i] <- pmax(lambda*ST - A, 0) * exp(-r * T)
    } else {
      payoff[i] <- pmax(A - lambda*ST, 0) * exp(-r * T)
    }
  }
  
  option_price <- mean(payoff)
  error <- sqrt(1 / (n - 1) * sum((payoff - option_price)^2))
  
  return(list(price = option_price, standard_error = error / sqrt(n)))
}

result <- asian_floating_option_mc(S0, lambda, r, sigma, T, n, m, "put", 0)
print(result)



# Antithetic floating Asian Options ------------------------------------------------

asian_option_An <- function(S0, lambda, r, sigma, T, n, m, option_type="call", b) {
  dt <- T/m  #Time step
  
  payoff <- numeric(n)
  for (i in 1:(n)) {
    S1 <- numeric(m+1)
    S2 <- numeric(m+1)
    S1[1] <- S0
    S2[1] <- S0
    for (j in 2:(m+1)) {
      Z <- rnorm(1)
      S1[j] <- S1[j-1] * exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*Z)
      S2[j] <- S2[j-1] * exp((r-0.5*sigma^2)*dt+sigma*sqrt(dt)*(-Z))
    }
    
    #Discrete arithmetic average
    if (b == 0) {
      A1 <- mean(S1)
      A2 <- mean(S2)
    } else {
      A1 <- mean(S1[-1])
      A2 <- mean(S2[-1])
    }
    ST1 <- S1[m + 1]  #for floating
    ST2 <- S2[m + 1]
    
    #Payoff
    if (option_type == "call") {
      payoff[i] <- exp(-r*T)* (pmax(lambda*ST1 - A1, 0) + pmax(lambda*ST2 - A2, 0))/2
    } else {
      payoff[i] <- exp(-r*T) * (pmax(A1 - lambda*ST1, 0) + pmax(A2 - lambda*ST2, 0))/2
    }
  }
  
  option_price <- mean(payoff)
  error <- sqrt(1/(n-1)*sum((payoff - option_price)^2))
  
  return(list(price=option_price, standard_error=error/sqrt(n)))
}


#Parameters

S0 <- 100
lambda <- 1
r <- 0.09
sigma <- 0.5
T <- 1
n <- 100000
m <- 120
b <- 0

resultAn <- asian_option_An(S0, lambda, r, sigma, T, n, m, "put",0)
print(resultAn)