#Parameters
T <- 1
N <- 100
S0 <- 50
sigma <- 0.15
intress <- 0.1
alfa <- 40
K <- 40
dt <- T/N
W <- sigma*sqrt(dt)

viimanes61m <- function(tyyp) {
  V <- c()
  S <- Spuu(N)
  A <- ApuuHW(N)
  if (tyyp == 1) {
    for (i in 1:length(S)) {
      for (j in 1:length(A)) {
        Vv <- max((A[j] - K), 0)
        V <- c(V, Vv)
      }
    }
  } else if (tyyp == 2) {
    for (i in 1:length(S)) {
      for (j in 1:length(A)) {
        Vv <- max((K - A[j]), 0)
        V <- c(V, Vv)
      }
    }
  } else if (tyyp == 3) {
    for (i in 1:length(S)) {
      for (j in 1:length(A)) {
        Vv <- max((S[i] - A[j]), 0)
        V <- c(V, Vv)
      }
    }
  } else if (tyyp == 4) {
    for (i in 1:length(S)) {
      for (j in 1:length(A)) {
        Vv <- max((A[j] - S[i]), 0)
        V <- c(V, Vv)
      }
    }
  }
  return(V)
}

Spuu <- function(N) {
  S <- c()
  if ((N + 1) %% 2 == 0) {
    for (i in seq(-N, -1, by=2)) {
      Ss <- S0 * exp(i * W)
      S <- c(S, Ss)
    }
    for (i in seq(1, N+1, by=2)) {
      Ss <- S0 * exp(i * W)
      S <- c(S, Ss)
    }
  } else {
    for (i in seq(-N, N+1, by=2)) {
      Ss <- S0 * exp(i * W)
      S <- c(S, Ss)
    }
  }
  return(S)
}

ApuuHW <- function(N) {
  A <- c()
  Y <- alfa * sqrt(0.25 / T) * sigma^2 * dt
  #Y <- alfa * sqrt(0.25 / T) * sigma * dt
  Amax <- S0 * (1 - exp(sigma * sqrt(dt) * (N + 1))) / ((N + 1) * (1 - exp(sigma * sqrt(dt))))
  Amin <- S0 * (1 - exp(-sigma * sqrt(dt) * (N + 1))) / ((N + 1) * (1 - exp(-sigma * sqrt(dt))))
  kmin <- floor((log(Amin) - log(S0)) / Y) - 1
  kmax <- floor((log(Amax) - log(S0)) / Y) + 1
  if (N == 0) {
    A <- c(A, S0)
  } else {
    for (i in seq(kmin, kmax, by=1)) {
      c <- S0 * exp(i * Y)
      A <- c(A, c)
    }
  }
  kmax-kmin
  return(A)
}


sammhaaval <- function(tyyp) {
  C <- c()
  W <- sigma * sqrt(dt)
  p <- (exp(intress * dt) - exp(-sigma * sqrt(dt))) /
    (exp(sigma * sqrt(dt)) - exp(-sigma * sqrt(dt)))
  V <- viimanes61m(tyyp)
  
  for (h in seq(N-1, 0, by=-1)) {
    C <- c()
    S <- Spuu(h)
    A <- ApuuHW(h)
    Ahiljem <- ApuuHW(h + 1)
    
    for (j in 1:length(S)) {
      for (k in 1:length(A)) {
        Aminus <- ((h + 1) * A[k] + (S[j] * exp(-W))) / (h + 2)
        Aplus <- ((h + 1) * A[k] + (S[j] * exp(W))) / (h + 2)
        
        
        a <- 1
        while (a <= length(Ahiljem) && Ahiljem[a] <= Aminus) {
          a <- a + 1
        }
        a <- min(a, length(Ahiljem))
        
        if (a > 1 && Ahiljem[a-1] == Aminus) {
          Aminusfloor <- Ahiljem[a-1]
          Aminusceil <- Ahiljem[a]
          Vminusfloor <- V[(j-1)*length(Ahiljem) + (a-1)]
          Vminusceil <- V[(j-1)*length(Ahiljem) + a]
        } else {
          if (a == 1) {
            Aminusfloor <- Ahiljem[1]
            Aminusceil <- Ahiljem[2]
            Vminusfloor <- V[(j-1)*length(Ahiljem) + 1]
            Vminusceil <- V[(j-1)*length(Ahiljem) + 2]
          } else {
            Aminusfloor <- Ahiljem[a-1]
            Aminusceil <- Ahiljem[a]
            Vminusfloor <- V[(j-1)*length(Ahiljem) + (a-1)]
            Vminusceil <- V[(j-1)*length(Ahiljem) + a]
          }
        }
        
        
        b <- 1
        while (b <= length(Ahiljem) && Ahiljem[b] <= Aplus) {
          b <- b + 1
        }
        b <- min(b, length(Ahiljem))
        
        if (b > 1 && Ahiljem[b-1] == Aplus) {
          Aplusfloor <- Ahiljem[b-1]
          Aplusceil <- Ahiljem[b]
          Vplusfloor <- V[j*length(Ahiljem) + (b-1)]
          Vplusceil <- V[j*length(Ahiljem) + b]
        } else {
          if (b == 1) {
            Aplusfloor <- Ahiljem[1]
            Aplusceil <- Ahiljem[2]
            Vplusfloor <- V[j*length(Ahiljem) + 1]
            Vplusceil <- V[j*length(Ahiljem) + 2]
          } else {
            Aplusfloor <- Ahiljem[b-1]
            Aplusceil <- Ahiljem[b]
            Vplusfloor <- V[j*length(Ahiljem) + (b-1)]
            Vplusceil <- V[j*length(Ahiljem) + b]
          }
        }
        
        eplus <- (Aplus - Aplusfloor) / (Aplusceil - Aplusfloor)
        emiinus <- (Aminus - Aminusfloor) / (Aminusceil - Aminusfloor)
        VV <- exp(-intress*dt)*(p*(eplus*Vplusceil + (1-eplus)*Vplusfloor) + 
                                      (1-p)*(emiinus*Vminusceil + (1-emiinus)*Vminusfloor))
        C <- c(C, VV)
      }
    }
    V <- C
  }
  return(V[1])  #return the option price
}


sammhaaval(1)