Implied volatility from American options using python

Question

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)