Integrals in code

Adding up thin slices

An integral measures the area under a curve: add up the value of f over many tiny steps. In finance integrals are everywhere in disguise: the cumulative profit of a position over time, the present value of a continuous cash flow, an expected payoff as the average of the payoff weighted by a probability density, or the total exposure of a portfolio over a period. Few real integrands have a tidy formula, so we approximate the area with slices. Split the interval [a, b] into n equal slices of width h = (b - a) / n and add up the slice areas.

# Area under x^2 on [0, 1] is exactly 1/3. Left rectangles, n slices:
def left_riemann(f, a, b, n):
    h = (b - a) / n
    return sum(f(a + i * h) for i in range(n)) * h

square = lambda x: x * x
for n in (10, 100, 1000):
    print(n, round(left_riemann(square, 0, 1, n), 6))
# 10 0.285
# 100 0.32835
# 1000 0.332834

Better slices: trapezoids and Simpson

Left rectangles are crude: each slice pretends the curve is flat. The trapezoid rule joins the two ends of each slice with a straight line, so the area of a slice is h * (f(left) + f(right)) / 2; adding them, the inner points count twice and the two end points once: h * (f(a) / 2 + f(a + h) + ... + f(b - h) + f(b) / 2). It is exact for straight lines. Simpson's rule fits a parabola through every three points: h / 3 * (f(a) + 4 * f(odd points) + 2 * f(even inner points) + f(b)), which needs an even n, and is exact for cubics. For smooth functions Simpson is usually far more accurate for the same work.

def trapezoid(f, a, b, n):
    h = (b - a) / n
    return h * (f(a) / 2 + sum(f(a + i * h) for i in range(1, n)) + f(b) / 2)

def simpson(f, a, b, n):
    if n % 2:
        raise ValueError("n must be even")
    h = (b - a) / n
    odd = sum(f(a + i * h) for i in range(1, n, 2))
    even = sum(f(a + i * h) for i in range(2, n, 2))
    return h / 3 * (f(a) + 4 * odd + 2 * even + f(b))

square = lambda x: x * x
for n in (10, 100, 1000):
    print(n, round(trapezoid(square, 0, 1, n), 6), round(simpson(square, 0, 1, n), 6))
# 10 0.335 0.333333
# 100 0.33335 0.333333
# 1000 0.333334 0.333333

print(simpson(lambda x: x**3, 0, 2, 4))   # 4.0  (exact for cubics)

Integrating data you already have

Often there is no function, only measurements at regular times: a list of exposures sampled every day. The trapezoid rule works directly on such samples. With values ys spaced dt apart, each neighbouring pair forms a trapezoid of area dt * (ys[i] + ys[i + 1]) / 2, and the sum approximates the area under the sampled curve, for example the total exposure-days. With fewer than two samples there is no interval, so the area is 0.0.

def trapezoid_samples(ys, dt):
    if len(ys) < 2:
        return 0.0
    return dt * (ys[0] / 2 + sum(ys[1:-1]) + ys[-1] / 2)

# Exposure sampled daily: 0, 2, 4, 2, 0 million
exposure = [0, 2, 4, 2, 0]
print(trapezoid_samples(exposure, 1.0))   # 8.0  exposure-days
print(trapezoid_samples([5.0], 1.0))      # 0.0

import math

# An expected value is an integral: E[X^2] for a standard normal is 1
def density(x):
    return math.exp(-x * x / 2) / math.sqrt(2 * math.pi)

def simpson(f, a, b, n):
    h = (b - a) / n
    odd = sum(f(a + i * h) for i in range(1, n, 2))
    even = sum(f(a + i * h) for i in range(2, n, 2))
    return h / 3 * (f(a) + 4 * odd + 2 * even + f(b))

second_moment = simpson(lambda x: x * x * density(x), -8, 8, 400)
print(round(second_moment, 6))                   # 1.0

# Probability that X is within one standard deviation: about 68%
print(round(simpson(density, -1, 1, 100), 4))    # 0.6827

• Slices of width h = (b - a) / n; left rectangles h * sum(f(a + i*h)) for i < n.
• Trapezoid: ends count half; Simpson: weights 1, 4, 2, 4, ..., 4, 1 times h / 3, n even.
• For sampled data use dt * (first / 2 + middle values + last / 2); test on known integrals.

Next up: vectors and dot products, the building block of portfolio returns, correlations and, later, NumPy.

Lab · what's graded

• Approximate an area with rectangles using left_riemann25 pts
• Do better with trapezoids using trapezoid25 pts
• Fit parabolas with simpson25 pts
• Integrate measured data with trapezoid_samples25 pts