用代码理解积分

把细条加起来

积分衡量曲线下方的面积：把 f 在许多微小步长上的值加起来。在金融里，积分无处不在，只是披着别的外衣：持仓随时间的累计盈亏、连续现金流的现值、以概率密度加权的收益期望，或者一段时期内投资组合的总敞口。很少有真实的被积函数能写成整齐的公式，所以我们用细条来近似面积。把区间 [a, b] 分成 n 个宽度为 h = (b - a) / n 的等宽细条，再把面积加起来。

# 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

更好的切法：梯形与辛普森

左矩形法很粗糙：每个细条都假装曲线是平的。梯形法用直线连接每个细条的两端，所以一个细条的面积是 h * (f(左) + f(右)) / 2；相加时内部的点算两次，两个端点各算一次：h * (f(a) / 2 + f(a + h) + ... + f(b - h) + f(b) / 2)。它对直线是精确的。辛普森法对每三个点拟合一条抛物线：h / 3 * (f(a) + 4 * f(奇数点) + 2 * f(内部偶数点) + f(b))，要求 n 为偶数，对三次函数是精确的。对光滑函数，在同样的工作量下辛普森通常精确得多。

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)

对已有数据积分

很多时候没有函数，只有按固定时间间隔的测量值：比如每天采样的敞口列表。梯形法可以直接用于这类样本。对间隔为 dt 的数值 ys，每对相邻值构成一个面积为 dt * (ys[i] + ys[i + 1]) / 2 的梯形，总和近似采样曲线下的面积，例如总的敞口天数。少于两个样本时没有区间，面积为 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

• 宽度 h = (b - a) / n 的细条；左矩形法为对 i < n 求 h * sum(f(a + i*h))。
• 梯形法：两端算一半；辛普森法：权重 1、4、2、4 ... 4、1 乘以 h / 3，n 为偶数。
• 对采样数据用 dt * (首项 / 2 + 中间各值 + 末项 / 2)；用已知积分来检验。

下一节：向量与点积，投资组合收益、相关性乃至之后 NumPy 的基本构件。

实验：评分内容

• 用 left_riemann 以矩形近似面积25 分
• 用 trapezoid 以梯形做得更好25 分
• 用 simpson 拟合抛物线25 分
• 用 trapezoid_samples 对测量数据积分25 分