Take the 2-minute tour ×
Quantitative Finance Stack Exchange is a question and answer site for finance professionals and academics. It's 100% free, no registration required.

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)

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

1 Answer 1

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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