1
$\begingroup$

I am trying to calculate the implied volatility using newton-raphson in python, but the value diverges instead of converge. What is wrong with the code?

s = stock price
k = strike
t = time to maturity
rf = risk free interest
cp = +/-1 call/put
price = option price

def newtonRap(cp, price, s, k, t, rf):
    v = sqrt(2*pi/t)*price/s
    print "initial volatility: ",v
    for i in range(1, 10):
        d1 = (log(s/k)+(rf+0.5*pow(v,2))*t)/(v*sqrt(t))
        d2 = d1 - v*sqrt(t)
        gamma = norm.pdf(d1)/(s*v*sqrt(t))
        price0 = cp*s*norm.cdf(cp*d1) - cp*k*exp(-rf*t)*norm.cdf(cp*d2)
        v = v - (price0 - price)/gamma
        print "price, gamma, volatility\n",(price0, gamma, v)
        if abs(price0 - price) < 1e-10 :
            break
    return v


v = newtonRap(cp=1, price = 1.52, s=23.95, k=24, t=71.0/365, rf=0.05)
print v

output:

initial volatility:  0.36069926906  
price, gamma, volatility  
(1.6055072570611344, 0.10385864414094476, -0.46260492259786345)  
price, gamma, volatility  
(-1.8488599102758396, -0.080851497020229368, -42.129859511995726)  
price, gamma, volatility  
(-23.767706818545953, -1.6137689013848907e-22, -1.5669967860233743e+23)  
price, gamma, volatility  
(-23.767706818545953, -0.0, -inf)  

RuntimeWarning: divide by zero encountered in double_scalars  
  v = v - (price0 - price)/gamma  
RuntimeWarning: invalid value encountered in double_scalars  
  d1 = (log(s/k)+(rf+0.5*pow(v,2))*t)/(v*sqrt(t))  
price, gamma, volatility  
(nan, nan, nan)  
price, gamma, volatility  
(nan, nan, nan)  
$\endgroup$
3
  • 2
    $\begingroup$ It's bee a while but I think you should be dividing by vega, not gamma. $\endgroup$
    – brian
    Commented May 22, 2014 at 19:00
  • 2
    $\begingroup$ Split your code in three functions, which you can test individually: the first function implements the Newton-Raphson method—test it on examples which are easier to understand—the second function implements the volatility function and the second its derivative. $\endgroup$ Commented May 23, 2014 at 5:18
  • 1
    $\begingroup$ Great help. Brian you are absolutely correct. It is vega instead of gamma. Changing just that made it work. $\endgroup$ Commented May 23, 2014 at 15:04

1 Answer 1

-1
$\begingroup$
def newtonRap(cp, price, s, k, t, rf):
    v = sqrt(2*pi/t)*price/s
    print "initial volatility: ",v
    for i in range(1, 100):
        d1 = (log(s/k)+(rf+0.5*pow(v,2))*t)/(v*sqrt(t))
        d2 = d1 - v*sqrt(t)
        vega = s*norm.pdf(d1)*sqrt(t)
        price0 = cp*s*norm.cdf(cp*d1) - cp*k*exp(-rf*t)*norm.cdf(cp*d2)
        v = v - (price0 - price)/vega
        print "price, vega, volatility\n",(price0, vega, v)
        if abs(price0 - price) < 1e-25 :
            break
    return v
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.