有限差分法

网格上的导数

导数是一个极限，但计算机只能用有限的步长。有限差分用小步长 h 上的斜率代替 f'(x)：前向 (f(x + h) - f(x)) / h，或中心 (f(x + h) - f(x - h)) / (2h)，后者更准确，因为两侧的误差相互抵消（误差按 h**2 缩小，而前向只按 h）。同样的思路给出二阶导数，也就是曲率：(f(x + h) - 2 f(x) + f(x - h)) / h**2。对作为现价函数的期权价格，它就是 gamma。

import math

def forward(f, x, h):
    return (f(x + h) - f(x)) / h

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

exact = math.cos(1.0)                               # the derivative of sin at 1
for h in (0.1, 0.01, 0.001):
    print(h, abs(forward(math.sin, 1.0, h) - exact), abs(central(math.sin, 1.0, h) - exact))

# Gamma of a call from the second difference of its Black-Scholes price
import math

def norm_cdf(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def bs_call(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))

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

price = lambda s: bs_call(s, 100, 0.05, 0.2, 1)
print(round(second_difference(price, 100, 0.5), 6))     # 0.018761
print(round(second_difference(price, 100, 0.01), 6))    # 0.018762

从导数到方程

这种方法的目的是求解偏微分方程。布莱克-斯科尔斯就是一个：期权价格随时间的变化取决于它对现价的曲率。经过变量替换，它变成热方程 du/dt = D * d2u/dx2：热量会扩散，温度曲线越弯，变化越快。把数值放在网格上，把两个导数都换成差分，方程就变成了根据邻点更新每个点的规则。

# The explicit scheme: a spike of heat spreads to its neighbours, step by step
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

grid = np.zeros(9)
grid[4] = 1.0
for step in range(4):
    print(step, np.round(grid, 3))
    grid = explicit_step(grid, 0.25)

代码中，alpha = D * dt / dx**2 把时间步长 dt 和网格间距 dx 联系起来。一次更新 new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2]) 用切片一次移动所有内部点，两个端点是边界条件，这里保持固定。这个格式很简单，陷阱也很简单。

稳定性：陷阱

显式格式只有在 alpha <= 0.5，即 D * dt / dx**2 <= 0.5 时才稳定。超过这个值，舍入误差每一步都被放大，数值会以锯齿状爆炸。这个条件很苛刻：网格间距减半，时间步长必须缩小到四分之一，成本迅速上升。所以专业代码常用隐式格式，它每一步求解一个小线性方程组，对任何 dt 都稳定；显式格式则最容易编写和理解。

# The same spike with alpha = 0.4 (stable) and alpha = 0.6 (unstable) after 60 steps
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

for alpha in (0.4, 0.6):
    grid = np.zeros(21)
    grid[10] = 1.0
    for _ in range(60):
        grid = explicit_step(grid, alpha)
    print(alpha, "largest value", float(np.abs(grid).max()))       # 0.4: 0.0556, 0.6: about 3.1e7

• 有限差分用网格上的值代替导数：中心差分 (f(x + h) - f(x - h)) / (2h) 的误差约为 h**2；二阶差分 (f(x + h) - 2 f(x) + f(x - h)) / h**2 给出曲率（gamma）。
• 偏微分方程变成网格上的更新规则：热方程的显式一步是 values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1])，用切片书写。
• 显式格式需要 D * dt / dx**2 <= 0.5 才稳定；隐式格式去掉了这一限制，代价是每步要解一个线性方程组。

下一节：协方差矩阵与相关矩阵，投资组合优化的第一步。

实验：评分内容

• 用 second_difference 估计曲率25 分
• 用 explicit_step 把网格推进一步25 分
• 用 is_stable 检查时间步长25 分
• 用 solve_heat 沿时间向前推进25 分