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
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)
}
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 iterationWe 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
