estimate implied volatility using newton-raphson in python

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 

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)

-
It's bee a while but I think you should be dividing by vega, not gamma. – brian May 22 '14 at 19:00
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. – Michael Grünewald May 23 '14 at 5:18
Great help. Brian you are absolutely correct. It is vega instead of gamma. Changing just that made it work. – user2686641 May 23 '14 at 15:04

 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