Derivatives in code

How fast does it change?

The derivative of a function f at x is its instantaneous rate of change: the slope of the curve at that point. In finance almost every sensitivity is a derivative. An option's delta is how much its price moves per unit move of the underlying, its gamma is how fast delta itself changes, and a bond's duration is a derivative of price with respect to yield. Calculus gives exact formulas for simple functions (the derivative of x**2 is 2*x), but real pricing functions are often too complicated or only available as code. Then we estimate the derivative numerically.

# A function is just a Python function
def f(x):
    return x * x

# Slope between two nearby points approaches the true slope (2x = 6 at x = 3)
x = 3.0
for h in (1.0, 0.1, 0.001):
    slope = (f(x + h) - f(x)) / h
    print(h, round(slope, 4))
# 1.0 7.0
# 0.1 6.1
# 0.001 6.001

Finite differences: forward and central

The forward difference looks one small step h ahead: (f(x + h) - f(x)) / h. Its error shrinks in proportion to h. The central difference looks a step each way: (f(x + h) - f(x - h)) / (2 * h). Errors on the two sides cancel, so its error shrinks like h squared and it is exact for quadratics. The second derivative (gamma) combines three points: (f(x + h) - 2 * f(x) + f(x - h)) / h**2. Pass the function itself as an argument (f), a skill from the functions topics.

import math

def forward_diff(f, x, h=1e-6):
    return (f(x + h) - f(x)) / h

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def second_derivative(f, x, h=1e-4):
    return (f(x + h) - 2 * f(x) + f(x - h)) / h**2

print(round(central_diff(math.sin, 0.0), 6))            # 1.0  (cos 0)
print(round(central_diff(lambda x: x**3, 2.0), 4))      # 12.0 (3 x^2)
print(round(second_derivative(lambda x: x**3, 2.0), 3)) # 12.0 (6 x)

# A quadratic: central is exact even for a big step, forward is not
print(central_diff(lambda x: x * x, 3.0, h=0.5))  # 6.0
print(forward_diff(lambda x: x * x, 3.0, h=0.5))  # 6.5

Choosing the step size

A smaller h reduces the truncation error of the formula, but subtracting two nearly equal numbers loses digits: this rounding error grows as h shrinks, and for a second derivative it grows faster because we also divide by h**2. Too small an h (like 1e-12) can give garbage. A good compromise is about 1e-5 for first derivatives with central differences and about 1e-4 for second derivatives. In real pricing code, "bump and reprice" is exactly this: re-run the pricer with the spot bumped up and down by a small amount.

import math

def central_diff(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

# True derivative of exp at 0 is 1; watch the error as h shrinks
for h in (1e-1, 1e-3, 1e-5, 1e-8, 1e-12):
    error = abs(central_diff(math.exp, 0.0, h) - 1.0)
    print(f"h={h:g} error={error:.1e}")
# h=0.1 error=1.7e-03
# h=0.001 error=1.7e-07
# h=1e-05 error=1.2e-11
# h=1e-08 error=5.3e-10     (rounding error takes over)
# h=1e-12 error=3.3e-05     (much worse: too small)

Using derivatives: Newton's method

A derivative tells you which way to move to reach a root, a point where f(x) = 0. Newton's method starts from a guess x and repeatedly replaces it with x - f(x) / f'(x), following the tangent line down to zero. It converges very quickly near a root, which is why it is used to back out an implied volatility from an option price. Stop when abs(f(x)) is below a tolerance, give up after a maximum number of iterations, and raise an error if the derivative is zero or the loop does not converge, rather than returning a wrong number.

def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def newton_root(f, x0, tol=1e-10, max_iter=50):
    x = x0
    for _ in range(max_iter):
        if abs(f(x)) < tol:
            return x
        slope = central_diff(f, x)
        if slope == 0:
            raise ValueError("flat slope: Newton cannot continue")
        x = x - f(x) / slope
    raise ValueError("did not converge")

# The square root of 2 is the positive root of x^2 - 2
root = newton_root(lambda x: x * x - 2, 1.0)
print(round(root, 8))           # 1.41421356

# x^2 + 1 never reaches zero
try:
    newton_root(lambda x: x * x + 1, 1.0)
except ValueError:
    print("no root found")      # no root found

• Forward (f(x+h) - f(x)) / h; central (f(x+h) - f(x-h)) / (2h); second (f(x+h) - 2f(x) + f(x-h)) / h**2.
• Not too big, not too small: about 1e-5 for first and 1e-4 for second derivatives.
• Newton: x = x - f(x) / f'(x) until abs(f(x)) is tiny; raise an error if it cannot converge.

Next up: integrals in code, adding up many thin slices to find areas, totals and expected values.

Lab · what's graded

• Estimate a slope looking ahead with forward_diff25 pts
• Estimate a slope from both sides with central_diff25 pts
• Estimate the curvature (gamma) with second_derivative25 pts
• Find where a function hits zero with newton_root25 pts