Binomial Tree Option Pricing › Complete Implementation

import numpy as np


class BinomialTree:
    """
    A comprehensive binomial tree implementation for option pricing.
    """

    def __init__(
        self, S0, K, T, r, sigma, N, option_type="call", american=False
    ):
        self.S0 = S0
        self.K = K
        self.T = T
        self.r = r
        self.sigma = sigma
        self.N = N
        self.option_type = option_type
        self.american = american

        # CRR parameters
        self.dt = T / N
        self.u = np.exp(sigma * np.sqrt(self.dt))
        self.d = 1 / self.u
        self.p = (np.exp(r * self.dt) - self.d) / (self.u - self.d)
        self.disc = np.exp(-r * self.dt)

        # Build tree
        self._build_tree()

    def _build_tree(self):
        """Construct stock price and option value trees."""
        N = self.N

        # Stock prices
        self.stock_tree = np.zeros((N + 1, N + 1))
        for i in range(N + 1):
            for j in range(i + 1):
                self.stock_tree[i, j] = (
                    self.S0 * (self.u**j) * (self.d ** (i - j))
                )

        # Option values
        self.option_tree = np.zeros((N + 1, N + 1))
        self.exercise_tree = np.zeros((N + 1, N + 1), dtype=bool)

        # Terminal payoffs
        for j in range(N + 1):
            S = self.stock_tree[N, j]
            if self.option_type == "call":
                self.option_tree[N, j] = max(S - self.K, 0)
            else:
                self.option_tree[N, j] = max(self.K - S, 0)

        # Backward induction
        for i in range(N - 1, -1, -1):
            for j in range(i + 1):
                # Continuation value
                cont = self.disc * (
                    self.p * self.option_tree[i + 1, j + 1]
                    + (1 - self.p) * self.option_tree[i + 1, j]
                )

                if self.american:
                    # Exercise value
                    S = self.stock_tree[i, j]
                    if self.option_type == "call":
                        exercise = max(S - self.K, 0)
                    else:
                        exercise = max(self.K - S, 0)

                    if exercise > cont:
                        self.option_tree[i, j] = exercise
                        self.exercise_tree[i, j] = True
                    else:
                        self.option_tree[i, j] = cont
                else:
                    self.option_tree[i, j] = cont

    @property
    def price(self):
        """Return the option price."""
        return self.option_tree[0, 0]

    @property
    def delta(self):
        """Calculate delta using the first step of the tree."""
        if self.N < 1:
            return None
        return (self.option_tree[1, 1] - self.option_tree[1, 0]) / (
            self.stock_tree[1, 1] - self.stock_tree[1, 0]
        )

    @property
    def gamma(self):
        """Calculate gamma using the second step of the tree."""
        if self.N < 2:
            return None
        # Delta at upper node after one step
        delta_up = (self.option_tree[2, 2] - self.option_tree[2, 1]) / (
            self.stock_tree[2, 2] - self.stock_tree[2, 1]
        )
        # Delta at lower node after one step
        delta_down = (self.option_tree[2, 1] - self.option_tree[2, 0]) / (
            self.stock_tree[2, 1] - self.stock_tree[2, 0]
        )
        # Average stock price change
        h = 0.5 * (self.stock_tree[2, 2] - self.stock_tree[2, 0])
        return (delta_up - delta_down) / h

    @property
    def theta(self):
        """Calculate theta using the second step of the tree."""
        if self.N < 2:
            return None
        # Option value at t=2dt for middle node
        f_2dt = self.option_tree[2, 1]
        # Per-year theta
        return (f_2dt - self.price) / (2 * self.dt)

# Demonstrate the BinomialTree class
S0, K, T, r, sigma = 100, 100, 1, 0.05, 0.2
N = 100

# Price European and American puts
euro_put = BinomialTree(S0, K, T, r, sigma, N, "put", american=False)
amer_put = BinomialTree(S0, K, T, r, sigma, N, "put", american=True)
premium = amer_put.price - euro_put.price