binomial_pricing.py

# -*- coding: utf-8 -*-
"""Binomial pricing.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1tzeboCjpLzrc-W6bI72nBMProrJxPcDw
"""

import pandas as pd
import pandas_datareader.data as web
import numpy as np
import datetime
import math
import yfinance as yf
import matplotlib.pyplot as plt

def additive_binomial_tree(n):
    for i in range(n):
        x = [1 , 0, 1]
        for j in range(i):
            x.append(0)
            x.append(1)
        x = np.array(x) + i
        y = np.arange(-(i+1), i+2)[::-1]
        y = y + 13
        #y = np.arange(n-1-i, n+2+i)[::-1]
        plt.plot(x, y, 'bo-', color='black');
        plt.xlabel('n');
        plt.ylabel('G');
    c = np.linspace(0,n,10)
    d = 0*c
    plt.plot(c, d, color='red')
    plt.show();

additive_binomial_tree(14)

def multiplicative_binomial_tree(n, a):
    seq = 3
    for i in range(n):
        x = [1, 0, 1]
        y = [n/a**(i+1)]
        for j in range(i):
            x.append(0)
            x.append(1)
        for b in range(seq-1):
          y.append(y[b]*a)
        seq += 2
        x = np.array(x) + i
        y = np.array(y)[::-1]
        plt.plot(x, y, '-', color='black');

    c = np.linspace(0,n,10)
    d = 0*c
    plt.plot(c, d, color='red')
    plt.xlabel('n');
    plt.ylabel('S');
    plt.show();

multiplicative_binomial_tree(13, 1.2)

"""## Inputs"""

T = 1
n = 10
dt = T / n
r = 0.05 * T
m = 1+r
u = 1.1
d = 1/1.1
p = (1+r-d)/(u-d)
q = 1-p
s0 = 100
K = 100
option_type = 'P'

"""## Binomial tree"""

def binomial_tree(n, u, d, s0):
    prices = np.zeros((n+1, n+1))
    prices[0, 0] = s0
    for i in range(1, n+1):
        prices[0:i, i] = prices[0:i, i-1] * u
        prices[i, i] = prices[i-1, i-1] * d
    return prices

stock_prices = binomial_tree(n, u, d, s0)
stock_prices

def draw_bin_tree(n, u, d, s0):
    tree = binomial_tree(n, u, d, s0)
    y = []
    for i in range(n+1):
        for j in tree[:, i]:
            if j != 0:
                y.append((j, i))
    y = np.array(y)
    plt.scatter(y[:, 1], y[:, 0])

draw_bin_tree(n, u, d, s0)

"""## Option payoff at maturity"""

def option_maturity_payoffs(stock_prices, K, option_type):
    if option_type == 'C':
        payoffs = stock_prices[:, len(stock_prices)-1] - K
    else:
        payoffs = K - stock_prices[:, len(stock_prices)-1]
    for j in range(len(payoffs)):
        payoffs[j] = max(0, payoffs[j])
    return payoffs

maturity_payoff = option_maturity_payoffs(stock_prices, K, option_type)
maturity_payoff

"""### Backward recursion to determine option price at time 0"""

def backward_payoffs(n, r, p, stock_prices, maturity_payoff):
    back_option_payoffs = np.zeros((n+1, n+1))
    back_option_payoffs[:, n] = maturity_payoff
    for i in range(n-1, -1, -1):
        for j in range(n-1, -1, -1):
            if stock_prices[i][j] != 0:
                back_option_payoffs[i][j] = (1/m)*(back_option_payoffs[i][j+1]*p + back_option_payoffs[i+1][j+1]*q) 
    return back_option_payoffs

b_payoffs = backward_payoffs(n, r, p, stock_prices, maturity_payoff)
b_payoffs

def draw_payoffs(b_payoffs):
    tree = b_payoffs
    y = []
    for i in range(n+1):
        for j in range(len(tree[:, i])):
            if tree[j, i] != 0 or  i>=j:
                y.append((tree[j, i], i))
    y = np.array(y)
    plt.scatter(y[:, 1], y[:, 0])

draw_payoffs(b_payoffs)

"""### American put options"""

def put_anticipated_payoffs(b_payoffs, K, stock_prices):
    intrinsic_value = np.zeros((n+1, n+1))
    intrinsic_value[:, n] = b_payoffs[:, n]
    for i in range(n-1, -1, -1):
        for j in range(n-1, -1, -1):
            if b_payoffs[i][j] != 0:
                intrinsic_value[i][j] = max(b_payoffs[i][j], K - stock_prices[i][j])
    intrinsic_value[0][0] = (1/m)*(p*intrinsic_value[1][0]+q*intrinsic_value[1][1])
    return intrinsic_value

b_payoffs = put_anticipated_payoffs(b_payoffs, K, stock_prices)
b_payoffs

draw_payoffs(b_payoffs)

"""## Practical usage"""

r = np.log(1+r)
m = np.exp(r*dt)
sigma = 0.2
u = np.exp(sigma*np.sqrt(dt))
d = 1/u
p = (m - d)/(u-d)
q = 1-p

def combination(n, k):
    return math.factorial(n) / (math.factorial(k)*math.factorial(n - k))
    
def binomial_f(k, n, p):
    return combination(n, k)*(p**k)*(1-p)**(n-k)

def payoff(n, p, u, d, S0, K, option_type):
    pay = 0
    for k in range(n+1):
        if option_type == 'call':
            pay += binomial_f(k, n, p) * max(0, S0*(u**k)*(d**(n-k)) - K)
        else:
            pay += binomial_f(k, n, p) * max(0, K - (u**k)*(d**(n-k))*S0)
    return np.exp(-r*T) * pay

price = payoff(n, p, u, d, s0, K, 'call')
price

def binomial_formula(S0, K , T, r, sigma, N, option_type):
    dt = T/N
    u = np.exp(sigma * np.sqrt(dt))
    d = np.exp(-sigma * np.sqrt(dt))
    p = (  np.exp(r*dt) - d )  /  (  u - d )
    return payoff(N, p, u, d, S0, K, option_type)

#def binom_EU1(S0, K , T, r, sigma, N, type_ = 'call'):
#    dt = T/N
#    u = np.exp(sigma * np.sqrt(dt))
#    d = np.exp(-sigma * np.sqrt(dt))
#    p = (  np.exp(r*dt) - d )  /  (  u - d )
#    value = 0 
#    for i in range(N+1):
#        node_prob = combination(N, i)*p**i*(1-p)**(N-i)
#        ST = S0*(u)**i*(d)**(N-i)
#        if type_ == 'call':
#            value += max(ST-K,0) * node_prob
#        elif type_ == 'put':
#            value += max(K-ST, 0)*node_prob
#
#    return value*np.exp(-r*T)
    
def get_data(symbol, n):
    obj = yf.Ticker(symbol)
    expiry_dates = obj.options
    options = obj.option_chain(expiry_dates[n])
    df = options.calls
    df.reset_index(inplace=True)
    df['Time'] = (datetime.datetime.strptime(expiry_dates[n], '%Y-%m-%d') - datetime.datetime.now()).days
    df['expiration'] = datetime.datetime.strptime(expiry_dates[n], '%Y-%m-%d')
    df['mid_price'] = (df.bid + df.ask) / 2
    return df

df = get_data('AAPL', 1)

prices = [] 

for row in df.itertuples():
    price = binomial_formula(165, row.strike, row.Time / 255, 0.03, 0.018 * np.sqrt(row.Time), 100, 'call')
    prices.append(price)
    
df['Price'] = prices
df['error'] = df.mid_price - df.Price 

plt.plot(df.strike, df.mid_price, label= 'Mid Price')
plt.plot(df.strike, df.Price, label = 'Calculated Price')
plt.xlabel('Strike')
plt.ylabel('Call Value')
plt.legend()

df.head()

plt.plot(df.strike, df.error);

prices = []
for i in range(1, 100):
  prices.append(binomial_formula(100, 100, 1, 0.05, 0.2, i, 'call'))

plt.plot(prices)

from scipy.stats import norm

N = norm.cdf
def BS_CALL(S, K, T, r, sigma):
    d1 = (np.log(S/K) + (r + sigma**2/2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * N(d1) - K * np.exp(-r*T)* N(d2)

prices1 = []
for i in range(1, 100):
  prices1.append(BS_CALL(100, 100, 1, 0.05, 0.2))

prices = []
for i in range(1, 100):
  prices.append(binomial_formula(100, 100, 1, 0.05, 0.2, i, 'call'))

plt.plot(prices)
plt.plot(prices1)

scarti = []
for i in range(1, 99):  
  scarti.append(np.array(prices)[i] - np.array(prices)[i-1])
plt.plot(scarti)

stds = []
N = 100
for j in range(2, N):
  prices = []
  for i in range(1, j):
    prices.append(binomial_formula(100, 100, 1, 0.05, 0.2, i, 'call'))

  somma = 0
  for i in range (j-1):
    somma += (prices[i] - sum(prices)/len(prices))**2
  variance = (1/len(prices)) * somma
  stds.append(variance ** 1/2)

plt.plot(stds)
plt.xlim(1, 100)