Implied volatility from American options using python

I am currently trying to construct volatility surface from american option prices (using Cox-Ross-Rubinstein tree) in Python 2.7. Below you can find the code I came up with. Any corrections would be welcome, however my biggest issue is performance. Single implied vol takes around 0.03 seconds to calculate. I guess that is because of the loops I use for moving backwards on the tree. Does anybody have an idea how to redesign the code to avoid looping for every step on the tree?

import numpy as np
from scipy.optimize import newton

# test parameters

n_steps=200    #number of steps
time=2.0    #time to maturity
strike=100.0    # strike price
price=102.0     # current stock price
time_step=float(time/n_steps) #length of a time step

R=0.05    #interest rate
D=0.02    #dividend yield
df_base=[1.0]    #discount factor for Today e.g. DF(TOday,Today)
df_base.extend([np.exp(-R*time_step)]*n_steps)
DFs=np.cumprod(df_base)    #vector of discount factors for all the steps

forwards=price*np.cumprod([np.exp((R-D)*time_step)]*n_steps)    # final version will make it possible to enter a vector of forward prices

call_price=17.3285914343
put_price=10.5343901817

# end of test parameters

def AmerImpliedVol(opt_price, price, forwards, strike, time,DFs,callPutInd=1.,n_steps=200,scale=1.):
time_step=float(time/n_steps)
time_sqrt=np.sqrt(time_step)

horizontal=np.array(xrange(0,n_steps+1),dtype='float64')
horizontal.shape=(1,n_steps+1)

vertical=np.array(xrange(0,-2*n_steps-2,-2),dtype='float64')
vertical.shape=(n_steps+1,1)

forwards_previous=[price]
forwards_previous.extend(forwards[:-1])
exp_drift=np.array(forwards)/np.array(forwards_previous)

step_DFs=np.array(DFs[1:])/np.array(DFs[:-1])

tree_template=vertical+horizontal

#vol=0.25
def step_price(vol):
power_vol=np.exp(vol*time_sqrt)
down=1/power_vol
p_up=(exp_drift-down)/(power_vol-down)
p_down=1.-p_up

tree=np.power(power_vol,tree_template)*price

cp=callPutInd
def payoff(values):
return np.maximum(cp*(values-strike),0.)

payoffs=payoff(tree[:,n_steps])
current=payoffs

for i in xrange(n_steps,0,-1):
after_step=step_DFs[i-1]*(p_up[i-1]*current[:n_steps]+p_down[i-1]*current[1:])
payoffs=payoff(tree[:,i-1])
current[:-1]=np.maximum(after_step,payoffs[:-1])
#current[:-1]=after_step

return scale*current[0]-opt_price

#return brentq(step_price,0.01,10.0,rtol=1.0e-8, maxiter=50)
return newton(step_price,0.3)

print AmerImpliedVol(put_price,price,forwards,strike,time,DFs,-1.)
print AmerImpliedVol(call_price,price,forwards,strike,time,DFs,1.)


For European options simple implied volatility procedure (almost no optimization) takes around 0.0014 s per option. I understand that implying vols from american option prices is slower but a factor of 20 seems large.

I fine-tuned code a little bit:
- I removed computationally intensive power function and replaced it with numpy cumprod
- I took payoff calculations outside of the loop
- I decreased accuracy to 0.5e-04, as higher accuracy could be considered noise

Now it takes around 0.015s per option.

def AmerImpliedVol2(opt_price, price, forwards, strike, time,DFs,callPutInd=1.,n_steps=200,scale=1.):
time_step=float(time/n_steps)
time_sqrt=np.sqrt(time_step)

forwards_previous=[price]
forwards_previous.extend(forwards[:-1])
exp_drift=np.array(forwards)/np.array(forwards_previous)

step_DFs=np.array(DFs[1:])/np.array(DFs[:-1])

tree_template=vertical+horizontal

#vol=0.25
def step_price(vol):
power_vol=np.exp(vol*time_sqrt)
down=1/power_vol
p_up=(exp_drift-down)/(power_vol-down)
p_down=1.-p_up

horizontal=np.ones(shape=(1,n_steps+1))*power_vol
horizontal[0,0]=1.
horizontal=np.cumprod(horizontal,axis=1)

vertical=np.ones(shape=(n_steps+1,1))/(power_vol*power_vol)
vertical[0,0]=1.
vertical=np.cumprod(vertical,axis=0)

tree=price*(vertical*horizontal)

cp=callPutInd
def payoff(values):
return np.maximum(cp*(values-strike),0.)

payoffs=payoff(tree)
del tree
current=payoffs[:,n_steps]

for i in xrange(n_steps,0,-1):
after_step=step_DFs[i-1]*(p_up[i-1]*current[:i]+p_down[i-1]*current[1:])
current=np.maximum(after_step,payoffs[:i,i-1])
#print current
#current[:-1]=after_step

return scale*current[0]-opt_price

seed=0.3

return newton(step_price,seed,tol=0.5e-04)

• Have you found a way to further optimize your code? Commented Jul 5, 2017 at 21:04