%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, M), where S is the stock price and M the running maximum
# Compute range first. (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,S0)] )
for n in range(N):
    dom[n+1] = set( (u*s, max(u*s,m)) for (s,m) in dom[n] )
    dom[n+1].update( (d*s, m) for (s,m) in dom[n] )

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

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

{(10, 10): 0.64}

f[1]

{(9.00, 10): 0.47, (11.00, 11.00): 0.80}

f[5]

{(16.11, 16.11): 1,
 (10.78, 11.00): 0.69,
 (13.18, 13.18): 1,
 (13.18, 13.18): 1,
 (13.18, 13.31): 1,
 (10.78, 11.98): 0.69,
 (7.22, 11.00): 0.19,
 (8.82, 10): 0.40,
 (10.78, 12.10): 1,
 (8.82, 10): 0.40,
 (8.82, 10.89): 0.40,
 (7.22, 10): 0.19,
 (8.82, 11.00): 0.40,
 (8.82, 11.00): 0.40,
 (10.78, 13.31): 1,
 (5.90, 10): 0.07,
 (10.78, 10.78): 0.69,
 (10.78, 10.89): 0.69,
 (8.82, 12.10): 1,
 (13.18, 14.64): 1}