Implied Volatility in R ( Newton Raphson )

Just how fast is Newton Raphson method. We will not be able to understand the power unless we compare it with other methods
I am trying to find the Implied Volatility of an option using 2 methods

1) BISECTION METHOD

This involves providing a limit for the variable in question and adjusting the lower and higher limit based on the output of the current iteration

We have an ordinary Black Scholes equation
BS = function(S, K, T, r, sig, type="C"){
  d1 = (log(S/K) + (r + sig^2/2)*T) / (sig*sqrt(T))
  d2 = d1 - sig*sqrt(T)
  if(type=="C"){
    value = S*pnorm(d1) - K*exp(-r*T)*pnorm(d2)
  }
  if(type=="P"){
    value = K*exp(-r*T)*pnorm(-d2) - S*pnorm(-d1)
  }
  return(value)
}

The Bisection method will try to reduce the err term in the code below

err  = BS(S,K,T,r,sig,type) - market

The function implied.vol will take these inputs and minimise the err term upto 1000 iterations

implied.vol =

  function(S, K, T, r, market, type){

    sig = 0.20

    sig.up = 1

    sig.down = 0.001

    count = 0

    err = BS(S, K, T, r, sig, type) - market 

    

    ## repeat until error is sufficiently small or counter hits 1000

    while(abs(err) > 0.00001 && count<1000 div="">

      if(err < 0){

        sig.down = sig

        sig = (sig.up + sig)/2

      }else{

        sig.up = sig

        sig = (sig.down + sig)/2

      }

      err = BS(S, K, T, r, sig, type) - market

      count = count + 1

    }

    print(c("Counter is ",count),sep='')

    ## return NA if counter hit 1000

    if(count==1000){

      return(NA)

    }else{

      return(sig)

    }

  }

I clocked the execution and took the following AAPL option

Oct 2 Expiry as on 28th September with a Strike Price of 110

startTime = Sys.time()

implied.vol(112.44,110.00,4/252,.0031,3.50,"C")

print(Sys.time() - startTime)

0.02500105 secs 17 iterations

2) NEWTON RAPHSON

We will recalculate the option price using Newton Raphson

NR = function(f,tolerance=.0000001,start=1,iteration=100)

{

  deltax=.0000001

  counter=1

  current=start

  arrayofSolutions=numeric(iteration)

  while(counter <=iteration)

  {

    df.dx=(f(start+deltax) -f(start))/deltax

    current=start - (f(start)/df.dx)

    #print(current)

    arrayofSolutions[counter]=current

    counter = counter + 1

    if(abs(current-start) < tolerance) break

    start=current

  }

  return(arrayofSolutions[1:counter-1])

}

BSVolatility = function(volatility)

{

  BS(112.44,110.00,4/252,.0031,volatility,"C") - 3.5

}

I clocked the timings

startTime = Sys.time()

NR(BSVolatility,.000001,1,100)

print(Sys.time() - startTime)

0.01400089 sec 4 iterations