In [1]:
%pylab inline
%precision 2
import sys, tqdm

Populating the interactive namespace from numpy and matplotlib
In [2]:
# Parameters
u = 1.1
d = .9
r = .01
N = 500

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

print( f'p={p:.2f}, q={q:.2f}' )

p=0.55, q=0.45

Note at time $n$ the range of $S_n$ is simply $\{S_0 u^k d^{n-k} | k \in \{0, \dots, n\}\}$. Using this fact, we don't need to compute the domain of $f_n$ numerically. While it might seem like a minor savings, it actually saves a lot! Computing the domain numerically involves round-off errors at every step. As a result some values that are $10^{-16}$ apart are often treated as distinct values. These cause the size of your domain to grow much faster. I timed runs of my code below (with $N=500$), and as you can see it runs in miliseconds. If you used the older code I provided, I'm guessing it will take much longer (and may even not run at all).

In [3]:
# Stock price at time n if k heads have been tossed.
S = lambda n, k: S0*u**k * d**(n-k)

def Rn(k, f):
    # Rollback operator. Discounted conditional expectation of the price
    # assuming k heads have been tossed.
    # f should be the price at the next time.
    return ( f[k+1]*p + f[k]*q )/(1+r)

def price_option( g ):
    # Computes the price of an American option with intrinsic value g(S_n),
    # and also the price of a European option with maturity N and payoff g(S_N)
    # We index the prices by the number of heads. That is, f[n][k] gives
    # the price at time n when there have been k heads.
    f_am = empty( N+1, dtype=object )
    f_eu = empty_like( f_am )
    f_am[N] = [ g(N, S(N,k) ) for k in range(N+1) ]
    f_eu[N] = [ g(N, S(N,k) ) for k in range(N+1) ]
    
    # tqdm just gives you a progress bar and times your loop.
    for n in tqdm.tqdm_notebook( range(N-1, -1, -1)):
        f_am[n] = [ max( g(n, S(n,k)), Rn(k, f_am[n+1]) ) for k in range(n+1) ]
        f_eu[n] = [ Rn(k, f_eu[n+1]) for k in range(n+1) ]
        
    return (f_am, f_eu)
In [4]:
S0  = 10
K = S0*(1+r)**N
(am_call, eu_call) = price_option( lambda n, s: max(s - K, 0) )
(am_put, eu_put) = price_option( lambda n, s: max(K-s, 0) )
(am_straddle, eu_straddle ) = price_option( lambda n, s: abs(s-K) )

Widget Javascript not detected.  It may not be installed or enabled properly.










    






Widget Javascript not detected.  It may not be installed or enabled properly.










    






 
 















    















    






Widget Javascript not detected.  It may not be installed or enabled properly.










    






 
 















    
























In [5]:



    


am_call[0][0], am_put[0][0], am_call[0][0] + am_put[0][0] - am_straddle[0][0]



    















    
Out[5]:








(7.33, 1437.73, 7.33)



















In [6]:



    


eu_call[0][0], eu_put[0][0], eu_call[0][0] + eu_put[0][0] - eu_straddle[0][0]



    















    
Out[6]:








(7.33, 7.33, -0.00)



















In [7]:



    


# Let's figure out something about the minimal optimal exercise time
# The first time the AFP equals the intrinsic value is σ^*

g_put = lambda n, s: max(K-s, 0) # Intrinsic value for put
def AFPdifference(n, f, g):
    return [f[n][k] - g(n, S(n,k)) for k in range(n+1)]



    













In [8]:



    


# Choose smaller numbers (for readability) and re-solve
N = 50
K = 10
(am_put, eu_put) = price_option( g_put )



    















    






Widget Javascript not detected.  It may not be installed or enabled properly.










    






 
 















    
























In [9]:



    


AFPdifference(0, am_put, g_put)



    















    
Out[9]:








[1.33]



















In [10]:



    


AFPdifference(1, am_put, g_put)



    















    
Out[10]:








[0.70, 1.05]



















In [11]:



    


AFPdifference(2, am_put, g_put)



    















    
Out[11]:








[0.26, 1.25, 0.82]



















In [12]:



    


AFPdifference(3, am_put, g_put)



    















    
Out[12]:








[0.02, 0.64, 1.06, 0.63]



















In [13]:



    


AFPdifference(4, am_put, g_put)



    















    
Out[13]:








[0.00, 0.23, 1.18, 0.83, 0.48]



















In [14]:



    


AFPdifference(5, am_put, g_put)



    















    
Out[14]:








[0.00, 0.01, 0.59, 1.08, 0.63, 0.35]





















With these numbers it looks like at time 4 you should exercise this option if the first four coins come up tails. Otherwise you should hold it longer.













In [15]:



    


AFPdifference(6, am_put, g_put)



    















    
Out[15]:








[0.00, 0.00, 0.19, 1.11, 0.84, 0.48, 0.26]





















At time 6 if you get 1 heads you should exercise the option. (If you got 0 heads you would have already exercised at time 4)













In [ ]: