estimación de la volatilidad implícita mediante newton-raphson en python

Question

Estoy tratando de calcular la volatilidad implícita utilizando newton-raphson en python, pero el valor diverge en lugar de converger. ¿Qué pasa con el código?

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

de salida:

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)

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