American Binomial Model in Python

Having written about pricing American-style options on a binomial tree in q, I thought it would be instructive to do the same in Python and NumPy. Here is the code:

import functools as ft
import numpy as np

def BPTree(n, S, u, d):
  r = [np.array([S])]
  for i in range(n):
    r.append(np.concatenate((r[-1][:1]*u, r[-1]*d)))
  return r

def GBM(R, P, S, T, r, b, v, n):
  t = float(T)/n
  u = np.exp(v * np.sqrt(t))
  d = 1./u
  p = (np.exp(b * t) - d)/(u - d)
  ptree = BPTree(n, S, u, d)[::-1]
  R_ = ft.partial(R, np.exp(-r*t), p)
  return reduce(R_, map(P, ptree))[0]

def American(D, p, a, b): return np.maximum(b, D*(a[:-1]*p + a[1:]*(1-p)))
def VP(S, K): return np.maximum(K - S, 0)
ABM = ft.partial(GBM, American)

There is a minor deviation from the q code: we are allowing d to be specified in BPTree. But otherwise, they are doing the same thing. Performance (as measured in ipython) isn’t too far-off either:

In [1]: from binomial import *
In [2]: %timeit ABM(ft.partial(VP,K=102.0), 100.0, 1.0, 0.08, 0.08, 0.2, 1000)
10 loops, best of 3: 38.4 ms per loop
In [3]: ABM(ft.partial(VP,K=102.0), 100.0, 1.0, 0.08, 0.08, 0.2, 1000)
Out[3]: 6.2215001602514555

Note the similarity between the q and Python code. The similarity is a result of using NumPy and functools which enabled Python to perform array-oriented computation and partial function application. We did use a loop in BPTree as Python/NumPy does not seem to have the same “scan” operation as q. I suppose we could have created a numpy.ufunc to use accumulate()… but the loop felt cleaner and more Pythonic.

American Binomial Model in q

American-style options may be exercised prior to expiry. To price them on a binomial tree, the full tree has to be constructed. The prices are then computed at each node through a backward reduction process starting from the prices at maturity. At each node, we take the maximum of the payoff at that node versus the price implied by the binomial model.

Constructing a binomial price tree is relatively easy in q:

BPTree:{[n;S;u] n{(x*y 0),y%x}[u]\1#S}  / binomial price tree

where n is the depth of the tree, S is the current price and u is the scale of the up-move. We simply let the scale of the down-move be 1/u and hence the y%x in the expression.

The general binomial model can then be implemented as follow:

GBM:{[R;P;S;T;r;b;v;n]                  / General Binomial Model (CRR)
 t:T%n;                                 / time interval
 u:exp v*sqrt t;                        / up; down is 1/u
 p:(exp[b*t]-1%u)%(u-1%u);              / probability of up
 ptree:reverse BPTree[n;S;u];           / reverse binomial price tree
 first R[exp[neg r*t];p] over P ptree }

where R is a reduction function, P is the payoff function, S is the current price, T is the time to maturity, r is the risk-free rate, b is the cost of carry, v is the volatility and n is the depth of the tree. For American and European options, the reduction functions may be expressed as:

American:{[D;p;a;b] max(b;D*(-1_a*p)+1_a*1-p)}
European:{[D;p;a;b] D*(-1_a*p)+1_a*1-p}

where D is the discount factor and p is the probability of an up-move. Consequently, we can express the American and European binomial models simply as:


Testing the code on an American vanilla put option (strike = 102; price = 100; time to maturity = 1 year; risk-free rate = cost of carry = 8%; volatility = 20%; depth of tree = 1000):

q) VP:{[S;K]max(K-S;0)}    / vanilla put: max(K-S,0)
q) \t show ABM[VP[;102];100;1;0.08;0.08;0.2;1000]

It took 31 ms to compute. Pretty nice for so little code.

Note: It turns out that this implementation of EBM is faster than the one in the previous post. The reason is that I avoided using the expensive xexp function this time round. Otherwise, the previous implementation should be faster since it only computes the payoffs at maturity and not the intermediate nodes.

European Binomial Model in q and Python

In the previous post, we created a binomial probability mass function (pmf). We can use that to easily evaluate European-style options:

EBM:{[P;S;K;T;r;b;v;n]                 / European Binomial Model (CRR)
 t:T%n;                                / time interval
 u:exp v*sqrt t;                       / up
 d:1%u;                                / down
 p:(exp[b*t]-d)%(u-d);                 / probability of up
 ns:til n+1;                           / 0, 1, 2, ..., n
 us:u xexp ns;                         / u**0, u**1, ...
 ds:d xexp ns;                         / d**0, d**1, ...
 Ss:S*ds*reverse us;                   / prices at tree leaves
 ps:pmf[n;p];                          / probabilities at tree leaves
 exp[neg r*T]*sum P[Ss;K]*ps }

Note that P is the payoff, S is the current price, K is the strike price, T is the time to maturity, r is the risk-free rate, v is the volatility, b is the cost of carry and n is the depth of the binomial tree. The Python version using NumPy and SciPy actually looks quite similar:

def EuropeanBinomialModel(P, S, K, T, r, b, v, n):
  n = int(n)
  t = float(T)/n                              # time interval
  u = np.exp(v * np.sqrt(t))                  # up
  d = 1/u                                     # down
  p = (np.exp(b*t)-d)/(u-d)                   # probability of up
  ns = np.arange(0, n+1, 1)                   # 0, 1, 2, ..., n
  us = u**ns                                  # u**0, u**1, ...
  ds = d**ns                                  # d**0, d**1, ...
  Ss = S*us*ds[::-1]                          # prices at leaves
  ps = binom_pmf(ns, n, p)                    # probabilities at leaves
  return np.exp(-r*T) * np.sum(P(Ss,K) * ps)

As we can see, both code has no explicit loops. This is possible in Python as NumPy and SciPy are array-oriented. NumPy and SciPy’s idea of “broadcasting” has some similarity with k/q’s concept of “atomic functions” (definition: a function f of any number of arguments is atomic if f is identical to f’).