American Binomial Pricer (with Dividends)

Prices an American call on a dividend-paying stock with a Cox–Ross–Rubinstein binomial tree. Dividends are handled by subtracting their present value from the spot at \(t=0\) (the “escrow” model), then re-adding the present value of any future dividends at each node so the early-exercise check uses the with-dividend stock price. Uses decimal.Decimal with 2000 digits of precision so the up/down factors don’t drift when steps get large.

Run it from the command line:

python american.py S K sigma r T dt [d1_amt d1_t] [d2_amt d2_t] ...

import sys
import math
from decimal import Decimal, getcontext

getcontext().prec = 2000
ctx = getcontext()


initialStockPrice = Decimal(sys.argv[1])
strikePrice       = Decimal(sys.argv[2])
volatility        = Decimal(sys.argv[3])
interestRate      = Decimal(sys.argv[4])
totalTime         = Decimal(sys.argv[5])
timeStep          = Decimal(sys.argv[6])


dividends = []
for i in range(7, len(sys.argv), 2):
    dividends.append((Decimal(sys.argv[i]), Decimal(sys.argv[i+1])))

initialStockDiv = initialStockPrice
for dividend in dividends:
    initialStockDiv -= dividend[0] * ctx.exp(-interestRate * (dividend[1]))

steps = int((totalTime / timeStep))


sqrt_dt = timeStep.sqrt()
up   = ctx.exp(volatility * sqrt_dt)
down = ctx.exp(-volatility * sqrt_dt)

p = (ctx.exp(interestRate * timeStep) - down) / (up - down)

one = Decimal(1)

numFinalNodes = steps + 1
topFinalNode = initialStockDiv * (up ** steps)

currPrice = []
prevPrice = []
currNodes = []
prevNodes = []
for j in range(numFinalNodes, 0, -1):
    if j == numFinalNodes:
        for i in range(j):
            if i == 0:
                nodeBelow = topFinalNode
            else:
                nodeBelow = down * (currNode / up)

            prevNodes.append(max(nodeBelow - strikePrice, Decimal(0)))
            prevPrice.append(nodeBelow)
            currNode = nodeBelow
        continue

    for i in range(j):
        eurOptionValue = ctx.exp(-interestRate * timeStep) * (
            p * prevNodes[i] + (one - p) * prevNodes[i + 1]
        )
        currStockPrice = prevPrice[i] / up
        currPrice.append(currStockPrice)

        currTime = (j - 1) * timeStep
        for dividend in dividends:
            if currTime < dividend[1]:
                currStockPrice = currStockPrice + dividend[0] * ctx.exp(
                    -interestRate * (dividend[1] - currTime)
                )

        currOption = currStockPrice - strikePrice
        currNodes.append(max(eurOptionValue, currOption, 0))

    prevNodes = currNodes
    prevPrice = currPrice
    currNodes = []
    currPrice = []

rounded = prevNodes[0].quantize(Decimal("0.000000000000001"))
print(rounded)