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$
  • 2
    $\begingroup$ It's bee a while but I think you should be dividing by vega, not gamma. $\endgroup$ – brian May 22 '14 at 19:00
  • 1
    $\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$ – Michael Le Barbier Grünewald May 23 '14 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$ – user2686641 May 23 '14 at 15:04
-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, privacy policy and cookie policy

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