## 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.

## 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’).

## Binomial distribution in q

Recently I was using SciPy’s `scipy.stats.binom.pmf(x,n,p)`. I though it would be great if I could have such a function in q. So a simple idea is to construct a binomial tree with probabilities attached. Recalling that a Pascal triangle is generated using `n{0+':x,0}\1`, I modified it to get:

```q)pmf:{[n;p]n{(0,y*1-x)+x*y,0}[p]/1#1f}
q)pmf[6;0.3]
0.000729 0.010206 0.059535 0.18522 0.324135 0.302526 0.117649
q)sum pmf[1000;0.3]
1f
```

What is great about this method is that it is stable. Compared to SciPy 0.7.0, it was more accurate too (it is a known issue that older SciPy has buggy binom.pmf):

```>>> scipy.stats.binom.pmf(range(0,41),40,0.3)[-5:]
>>> array([3.33066907e-15, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.11022302e-16])
q)-5#pmf[40;0.7]
3.293487e-15 1.52594e-16 5.162955e-18 1.134715e-19 1.215767e-21
```

Unfortunately this method is too slow for large n. For large n, we need more sophisticated methods. For the interested reader, take a look at Catherine Loader’s Fast and Accurate Computation of Binomial Probabilities paper and an implementation of a binomial distribution in Boost

Posted in coding, k&q, python. Tags: , , . 1 Comment »