Newton Raphson Algorithm in R for Implied volatility

Question

I have the Black Scholes formula for options pricing in R :

BS = function(Flag,St, K, D, r, Ti, sigma) {
  d1  =  (log(St/K) + (r - D + (sigma^2)/2)*Ti) / (sigma*sqrt(Ti))
  d2  =  d1- sigma*sqrt(Ti)
  
  if(Flag == "call") price = St*exp(-D*(Ti)) * pnorm(d1)  - K*exp(-r*Ti)*pnorm(d2)
  if(Flag != "call") price = K*exp(-r*Ti)*pnorm(-d2)-St*exp(-D*Ti)*pnorm(-d1)
  return(price)}
BS("call",St=505.15, K=500, D=0, r=0.033, Ti=33/250, sigma=0.2)
[1] 18.48827

The first derivative of the above formula with respect to sigma :

vega_BS = function(St, K, D, r, Ti, sigma){
  d1  =  (log(St/K) + (r - D + (sigma^2)/2)*Ti) / (sigma*sqrt(Ti))
  vega =  St * dnorm(d1) * sqrt(Ti)
  return(round(vega,4))
}

I want to calculate the implied volatility for a given market value 32.4.

Doing so :

sig_implied = function(St, K, D,r, Ti,sigma,Market) {
  root_find = function(sigma){
    BS("call",St, K,D,r, Ti, sigma) - Market}
  round(uniroot(root_find, c(0,1))$root,3)
}
Market = 32.4
sig_implied(St=505.15, K=500, r=0.033,D=0, Ti=33/250,sigma=0.2,Market=Market)
[1] 0.394

Now I want to implement the NR algorithm for the implied volatility calculation.The result of NR has to be near 0.394 but doing so is far from close :

ImpliedVolNewton = function(Market,Flag, St, K, Ti, r, D,sigma, tol=0.0001, maxiter = 100) {

  s = 0.3
  not_converged = Ti
  vega = vega_BS(St, K, D, r, Ti, sigma)
  i = 1
  while (not_converged & (i < maxiter)) {
    err = (Market - BS(Flag,St, K, D, r, Ti, sigma) ) 
    s  =  s + err/vega
    not_converged = (abs(err/vega) > tol)
    i = i + 1
  }
  s }

ImpliedVolNewton(Market=32.4,"call",St=505.15, K=500, Ti=33/250, r=0.033, D=0,sigma=0.2,tol=0.0001)
[1] 22.73685

What I am doing wrong here ?

Any help ?

Update edit Or even this does not work

implied_volatility = function(Market,Flag,St,K,Ti,r,D,sigma,tol=0.0001,max_iterations=100){
  sigma0 = sqrt(abs(log(St/K)+r*Ti)*(2/Ti))
  for(i in max_iterations){
    diff = BS(Flag,St,K,Ti,r,D,sigma)-Market
    if(abs(diff)<tol){
      break
    }
    Sigma = sigma0 -diff/vega_BS(St,K,r,D,Ti,sigma)
  }
  return(Sigma)
}
implied_volatility (Market=32.4,"call",St=505.15, K=500, Ti=33/250, r=0.033, D=0,sigma=0.2,tol=0.0001)

Answer 1

implied_volatility_NR= function(Market,Flag, St, K,D,sigma, Ti, r) {
  max_iterations = 100
  tolerance      = 0.0001
  
  # Manaster & Koehler Seed Value
  sigma  =  sqrt(abs(log(St/K)+r*Ti)*(2/Ti))
  
  for(i in 1: max_iterations){

    vega  = vega_BS(St,K,r,D,Ti,sigma)
    diff  = Market - BS(Flag,St, K, D, r, Ti, sigma)
    
    if (abs(diff) < tolerance | vega < tolerance) {
      return(sigma)
    }
    sigma = sigma + diff/vega 
  }
  return(sigma)
}    
>implied_volatility_NR(32.4,Flag="call", St=505.15, K=500,D=0,sigma=0.2, Ti=33/250, r=0.033)
[1] 0.3936608