%pylab inline
%precision 2

Populating the interactive namespace from numpy and matplotlib

'%.2f'

# Parameters
u = 1.1
d = .9
S0 = 10
r = .01
N = 20

# Risk neutral probabilities
p = (1 + r - d) / (u-d)
q = (u - 1 -r) / (u - d)

# Up rebate option: Pays A the first time stock price crosses U
A = 1
U = 1.2 * S0

# Choose the state process Y = S, where S
# Compute range first. We only need the stock prices when S < U.
# (This requires a forward loop from 0 to N)
#
# Beware -- for more complicated models you may have to worry about round off errors
dom = empty( N+1, dtype=object )
dom[0] = set( [S0] ) # Assumes S0 < U
for n in range(N):
    dom[n+1] = set()
    for s in dom[n]:
        dom[n+1].add( d*s )
        if u*s < U: dom[n+1].add( u*s )

# Compute price. f[n][]s] gives the price at time n when stock price is s and the option has not paid
f = empty( N+1, dtype=object )
f[N] = { s:0 for s in dom[N] }
def Rn(n, s):
    # Rollback operator
    return ( (f[n+1][u*s] if u*s < U else A )*p + f[n+1][d*s]*q )/(1+r)
for n in range(N-1, -1, -1):
    f[n] = { s: Rn(n, s) for s in dom[n] }

# Print a few values of the arbitrage free price
f[0]

{10: 0.64}

f[1]

{9.00: 0.47, 11.00: 0.80}

f[5]

{5.90: 0.07, 7.22: 0.19, 8.82: 0.40, 8.82: 0.40, 10.78: 0.69}