Monday, September 28, 2015

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

Newton-Raphson method (univariate)


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

No comments: