In my previous post Options and Volatility Smile , we used Black-Scholes formula to derive Implied Volatility from given option strike and tenor. Most of the options traded on exchanges are American (with a few index options being European) and can be exercised at any time prior to expiration. Whether it is optimal to exercise an option early depends on whether the stock pays dividend or the level of interest rates and is a very complex subject. What I want to focus on is using Binomial model to price an American option. I summarize below Binomial model theory from the excellent book Derivative Markets 3rd Ed. by Robert McDonald

There are many flavors of the Binomial model but they all have following steps in common:

- Simulate future prices of underlying stock at various points in time until expiry
- Calculate the payoff of the option at expiry
- Discount the payoff back to today to get the option price

#### 1. Simulate future prices of the underlying

Assume the price of a stock is $S_0$ and over a small time period $\Delta t$, the price could either be $S_0u$ or $S_0d$ where u is a factor by which price rises and d is a factor by which price drops. The stock is assumed to follow a random walk, also assume that $p$ is the probability of the stock price rising and $(1-p)$ is the probability of it falling. There are many ways to approach the values of $u$,$d$ and $p$ and the various Binomial models differ in the ways these three parameters are calculated.

In Cox-Ross-Rubinstein (CRR) model, $u = \frac{1}{d}$ is assumed. Since we have 3 unknowns, 2 more equations are needed and those come from risk neutral pricing assumption. Over a small $\Delta t$ the expected return of the stock is

$$pu + (1-p)d = e^{r \Delta t}$$

and the expected variance of the returns is

$$pu^2 + (1-p)d^2 – (e^{r \Delta t})^2 = \sigma ^2 \Delta t$$

Solving for $u$, $d$ and $p$, we get

$$p = \frac{e^{r\Delta t} – d}{u-d}$$

$$u = e^{\sigma \sqrt{\Delta t}}$$

$$d = e^{-\sigma\sqrt{\Delta t}}$$

The CRR model generates a following tree as we simulate multi step stock price movements, this a recombining tree centered around $S_0$.

2. Calculating payoffs at expiry

In this step, we calculate the payoff at each node that corresponds to expiry.

For a put option, $payoff = max(K – S_N, 0)$

For a call option, $payoff = max(S_N – K, 0)$

where $N$ is node at expiry with a stock price $S_N$ and $K$ is the strike.

#### 3. Discounting the payoffs

In this step, we discount the payoffs at expiry back to today using *backward induction* where we start at expiry node and step backwards through time calculating option value at each node of the tree.

For American put, $V_n = max(K – S_n, e^{-r \Delta t} (p V_u + (1-p) V_d))$

For American call, $V_n = max(S_n – K, e^{-r \Delta t} (p V_u + (1-p) V_d))$

$V_n$ is the option value at node n

$S_n$ is the stock price at node n

$r$ is risk free interest rate

$\Delta t$ is time step

$V_u$ is the option value from the upper node at n+1

$V_d$ is the option value from the lower node at n+1

All the variants of Binomial model, including CRR, converge to Black-Scholes in the limit $\Delta t \to 0$ but the rate of convergence is different. The variant of Binomial model that I would like to use is called Leisen-Reimer Model which converges much faster. Please see the original paper for formulas and a C++ implementation at Volopta.com which I have ported to Python in the next section.

The python code is going to look very similar to Options and Volatility Smile post except that we will swap out Black-Scholes framework with Leisen-Reimer model. We will also use the same option chain data AAPL_BBG_vols.csv

##### A note on Python code

*I usually do not write code like below, I am purposely avoiding using any classes so that the focus remains on the objective which is to understand the technique.*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 |
import math import numpy as np import matplotlib as mpl import pandas as pd import matplotlib.pyplot as plt from scipy.integrate import quad from scipy.optimize import brentq from scipy.interpolate import interp2d from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm def BinomialLeisenReimer(S0, K, r, T, N, sigma, div): is_call = False odd_N = N if (N % 2 == 1) else (N + 1) d1 = (math.log(S0 / K) + ((r - div) + (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T)) d2 = (math.log(S0 / K) + ((r - div) - (sigma ** 2) / 2.) * T) / (sigma * math.sqrt(T)) dt = T / odd_N r_n = math.exp((r - div) * dt) method = 2. Term1 = math.pow((d1 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0) pp = 0.5 + math.copysign(1, d1) * 0.5 * math.sqrt(1 - math.exp(-Term1)) Term2 = math.pow((d2 / (N + 1.0 / 3.0 - (1 - method) * 0.1 / (N + 1))), 2) * (N + 1.0 / 6.0) p = 0.5 + math.copysign(1, d2) * 0.5 * math.sqrt(1 - math.exp(-Term2)) u = r_n * pp / p d = (r_n - p * u) / (1 - p) qu = p qd = 1 - p #Initialize price tree STs = [np.array([S0])] # Simulate the possible stock prices path for i in range(N): prev_branches = STs[-1] st = np.concatenate((prev_branches * u, [prev_branches[-1] * d])) STs.append(st) # Add nodes at each time step #Initialize payoff tree payoffs = np.maximum(0, (STs[N]-K) if is_call else (K-STs[N])) def check_early_exercise(payoffs, node): early_ex_payoff = (STs[node] - K) if is_call else (K - STs[node]) return np.maximum(payoffs, early_ex_payoff) #Begin tree traversal for i in reversed(range(N)): # The payoffs from NOT exercising the option payoffs = (payoffs[:-1] * qu + payoffs[1:] * qd) * df # Payoffs from exercising, for American options early_ex_payoff = (STs[i] - K) if is_call else (K - STs[i]) payoffs = np.maximum(payoffs, early_ex_payoff) return payoffs[0] headers = ['Date','Strike','Call_Bid','Call_Ask','Call','Maturity','Put_Bid','Put_Ask','Put'] dtypes = {'Date': 'str', 'Strike': 'float', 'Call_Bid': 'float', 'Call_Ask': 'float', 'Call':'float','Maturity':'str','Put_Bid':'float','Put_Ask':'float','Put':'float'} parse_dates = ['Date', 'Maturity'] data = pd.read_csv('AAPL_BBG_vols.csv',skiprows=1,header=None, names=headers, dtype=dtypes, parse_dates=parse_dates) data['LR_Imp_Vol'] = 0.0 data['TTM'] = 0.0 r = 0.0248 #risk-free rate S0 = 153.97 # spot price as of 5/11/2017 div = 0.0182 for i in data.index: t = data['Date'][i] T = data['Maturity'][i] K = data['Strike'][i] Put = data['Put'][i] sigma_guess = 0.2 time_to_maturity = (T - t).days/365.0 data.loc[i, 'TTM'] = time_to_maturity def func_LR(sigma): return BinomialLeisenReimer(S0, K, r, time_to_maturity, ((T - t).days) - 1, sigma, div) - Put lr_imp_vol = brentq(func_LR, 0.03, 1.0) data.loc[i, 'LR_Imp_Vol'] = lr_imp_vol ttm = data['TTM'].tolist() strikes = data['Strike'].tolist() imp_vol = data['LR_Imp_Vol'].tolist() f = interp2d(strikes,ttm,imp_vol, kind='cubic') plot_strikes = np.linspace(data['Strike'].min(), data['Strike'].max(),25) plot_ttm = np.linspace(0, data['TTM'].max(), 25) fig = plt.figure() ax = fig.gca(projection='3d') X, Y = np.meshgrid(plot_strikes, plot_ttm) Z = np.array([f(x,y) for xr, yr in zip(X, Y) for x, y in zip(xr,yr) ]).reshape(len(X), len(X[0])) surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0.1) ax.set_xlabel('Strikes') ax.set_ylabel('Time-to-Maturity') ax.set_zlabel('Implied Volatility') fig.colorbar(surf, shrink=0.5, aspect=5) |