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