Here are a few examples on pricing securities with a fixed maturity time. This code is mainly for illustrative purposes, and not optimized to run well.

%pylab inline
%precision 3

Populating the interactive namespace from numpy and matplotlib

'%.3f'

# Fix a few parameters for all examples.
u = 1.05
d = .9
S0 = 10
r = .03
N = 90

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

# Price a European Call option with strike K
K = 100
g = lambda s: max(s-K, 0) # Payoff of the call

# Choose our state process Y = S. First find Range(Y_n) with a forward in time loop.
# Then find f_n(Y_n) with a backward in time loop.

# In python range is a pre-defined function, so we use dom to denote the Range(Y_n)
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 )
        dom[n+1].add( u*s )
        
# Now dom[n] is Range(S_n). Let's compute f by backward induction.
# It's convenient to define a function that computes the AFP at time n, given the AFP at time n+1
# This is often called the "rollback operator" in numerical libraries
def Rn(n, s):
    return (f[n+1][u*s]*p + f[n+1][d*s]*q) / (1+r)

# Now run a loop backwards to compute f[n] starting from n = N
f = empty( N+1, dtype=object )
f[N] = { s:g(s) for s in dom[N] }
for n in range(N-1, -1, -1):
    f[n] = { s: Rn(n, s) for s in dom[n] }

# That's it, now print a few values of the AFP

print( f'f[0]={f[0]}' )
print( f'f[1]={f[1]}' ) # Should print a list of prices, depending on the stock price

f[0]={10: 3.576314189888773}
f[1]={9.0: 2.626151496736848, 10.5: 3.8462885569467575}

Look simple? Well, there are a few traps we shouldn't fall for. For instance, we know range(S_n) has n+1 elements in it. Let's see how many elements are in our computed $\operatorname{range}(S_n)$ for $n=30$

len(dom[30])

283

Woops. I got 283 on my system. I should have gotten 31! What went wrong?

Turns out it's floating point errors. Mathematically $u(u^m d^n) = u^{m+1} d^n$ for any $u, d \in \mathbb R$. On computers this is false, due to round off errors. Let's check. (Note in Python $u^m$ is written as u**m.)

u**2 * d**5 == u*(u * d**5)

False

That gave false on my system, even though mathematically the identity is true...

How do you avoid this? There are a few ways. Here's a suggestion for one workaround: Instead of storing the floating point number $u^i d^j$ inside dom[n], why don't we just store the integers $(i, j)$. (Or even just $i$, since $j$ can be worked out from $i$ and $n$.) This will avoid all floating point errors. See if you can do this (or find a different approach that avoids these floating point errors).

If you can't find a way to do this, no worries: the floating error you make is super small (about $10^{-17}$). So the price of your security will likely be correct to a few decimal places. The real cost is the run time! Due to round off errors, dom[n] will grow much faster than it should. So if we wanted $N = 1000$, say, this algorithm may not finish as quickly as you would like. Try setting $N = 1000$ above for instance. (Good code will handle that instantly; bad code may not finish...)

# Lets price an up and out call option, with the same strike price and upper barrier 500
# If the stock price ever exceeds or equals U, the option is worthless. If not, its a regular call.
K = 100
U = 500

We had an algorithm to price this in class. Choose the state process $Y = (S, M)$, where $M$ is the running maximum, and write the price in the form $V_n = f_n( S_n, M_n)$ and use the recurrence relation $$ f_{n}(s,m) = \frac{1}{1+r} \bigl( \tilde p f_{n+1}(us, m \vee (us) ) + \tilde q f_{n+1}(ds, m \vee (ds) \bigr) $$ To code this up, note that if $m \geq U$, the price is always $0$. So We only need to find the price when $m < U$. We also only need to keep track of the stock price when it is below $U$. Turns out, we won't even have to keep track of the range of $M$, and can simplify the above to $$ f_n(s) = \frac{1}{1+r} \bigl( \tilde p f_{n+1}(us) 1_{us < U} + \tilde q f_{n+1}(ds) 1_{ds < U} \bigr) $$

# Note: I'm going to use the same bad code as above. It's a bit easier to read, but it is still
# susseptable to floating point errors, and will choke if N is too large.


# Find the range(Y_n) first.
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]:
        if d*s < U: 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,
f = empty( N+1, dtype=object )
f[N] = { s:max(s - K, 0) for s in dom[N] } # Payoff if M < U
def Rn(n, s):
    # Rollback operator
    return (  p*(f[n+1][u*s] if u*s < U else 0 )
            + q*(f[n+1][d*s] if d*s < U else 0 )
           ) / (1+r)
for n in range(N-1, -1, -1):
    f[n] = { s: Rn(n, s) for s in dom[n] }
    
print( f'f[0]={f[0]}' )
print( f'f[1]={f[1]}' ) # Should print a list of prices, depending on the stock price

f[0]={10: 3.534795277316776}
f[1]={9.0: 2.616497274172393, 10.5: 3.798430191246108}