Valuing American Options using Tilley algorithm

Question

Hey I want to implement Tilley's algorithm (Valuing American Options in a Path Simulation Model by JA Tilley, 1993) to price american options. Where can I find implementation of this method in any programming language (Python preferable)? Is this method popular or rather not? I am not good programmer so some ready sources could be helpful.

I tried to do it but my result is wrong. Can anyone look at my code (in R)? I think that something must be wrong with step 6 or $z$ matrix.

S0 = 40
r = 0.07
sigma = 0.3
K = 45
n_steps_per_year = 252
dt = 1 / n_steps_per_year
T = 3
n_steps = n_steps_per_year * T
R = n_paths
Q = 70
P = 72
n_paths = P * Q
d = exp(-r * dt)

N = matrix(rnorm(n_paths * n_steps, mean = 0, sd = 1), n_paths, n_steps)

paths_S = matrix(nrow = n_paths, ncol = n_steps + 1, S0)  

for(i in 1:n_paths){
  for(j in 1:n_steps){
    paths_S[i, j + 1] = paths_S[i, j] * exp((r - 0.5 * sigma ^ 2) * dt + sigma * sqrt(dt) * N[i, j])
  }
}


I = matrix(nrow = n_paths, ncol = n_steps + 1)  # Intristic value matrix
I[, n_steps + 1] = sapply(K - paths_S[, n_steps + 1], max, 0)  
V = matrix(nrow = n_paths, ncol = n_steps + 1)  # option value matrix
V[, n_steps + 1] = I[, n_steps + 1] 
y = matrix(nrow = n_paths, ncol = n_steps + 1)  
z = matrix(nrow = n_paths, ncol = n_steps + 1)
H = matrix(nrow = n_paths, ncol = n_steps + 1)  # Holding value
  
# At moment T
for(k in 1:n_paths){
  z[k, n_steps + 1] = ifelse(I[k, n_steps + 1] > 0, 1, 0)   
}

### Backward-Induction ###
paths_S = as.data.frame(paths_S)

for(t in n_steps:2){
  # sorting
  o = order(paths_S[, t], decreasing = TRUE)
  paths_S = paths_S[o, ]
  H = H[o, ]
  y = y[o, ]
  V = V[o, ] 
  I = I[o, ]
  # Intristic value
  I[, t] = sapply(K - paths_S[, t], max, 0)

  # holding value 
  for(k in 1:n_paths){
    H[k, t] = d / P * sum(V[(floor((k - 1) / P) * P + 1):(floor((k - 1) / P) * P + P), t + 1])
  }
  
  x = as.integer(I[, t] > H[, t])
  
  # Step 6 #
  if(max(x) == 0){
    k_star = n_paths + 1
  }
  else{
    k_star = which.max(with(rle(x), rep(values * (lengths > c(tail(lengths, - 1), 0)), lengths)))
  }
  #  
  # k_star = which.max(with(rle(x), rep(values * (lengths > c(tail(lengths, - 1), 0)), lengths)))
  for(k in 1:n_paths){
    y[k, t] = ifelse(k >= k_star, 1, 0)
  }
  
  # value
  for(k in 1:n_paths){
    V[k, t] = ifelse(y[k, t] == 1, I[k, t], H[k, t])
  }
}

# z matrix
for(k in 1:n_paths){
  for(t in 2:(n_steps)){
    z[k, t] = ifelse(y[k, t] == 1 & max(y[k, 2:(t - 1)]) == 0, 1, 0)
  }
}


sum = 0
for(k in 1:n_paths){
  for(t in 2:(n_steps + 1)){
    sum = sum + z[k, t] * exp(-r * (t - 1) * dt) * I[k, t]
  }
}

price = sum / n_paths
price