价格按百分比变动
算术布朗运动不论股价是 5 还是 500，都加上同样金额的噪声，而且可能降到零以下。股票的表现不同：1% 的变动在任何价格上都是 1%。几何布朗运动把随机性放在百分比变化上：dS = mu * S * dt + sigma * S * dW。这就是布莱克-斯科尔斯所用的模型，mu 是年化预期收益，sigma 是波动率，也就是市场报价的那个量。
import math
def gbm_step(s, mu, sigma, dt, z):
return s * math.exp((mu - 0.5 * sigma ** 2) * dt + sigma * math.sqrt(dt) * z)
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 0.0), 4)) # a calm day: only the (mu - sigma^2/2) drift acts
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 1.0), 4)) # a one-sigma up day
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, -50.0), 4)) # an absurd fall: still positive精确解与负二分之一 σ 平方
解这个方程得到 S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t))。为什么有这个奇怪的 - sigma**2 / 2？因为指数函数是凸的：由詹森不等式，噪声会抬高 exp(噪声) 的平均，减去 sigma**2 / 2 恰好抵消它，所以*期望*价格以速率 mu 增长。典型的*路径*增长更慢，速率为 mu - sigma**2 / 2，这就是股价中位数低于均值的原因。下一节的伊藤引理会推导这一项。
import numpy as np
def gbm_path(s0, mu, sigma, dt, z):
log_moves = (mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z
return s0 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))
rng = np.random.default_rng(1)
z = rng.standard_normal(252)
path = gbm_path(100.0, 0.08, 0.2, 1 / 252, z)
print(path[:4].round(3), round(float(path[-1]), 3))大量路径：均值与中位数
模拟数千条一年期路径，看它们停在哪里。期末价格的*均值*接近 S0 * exp(mu * T)，即期望价格的公式。*中位数*更低，为 S0 * exp((mu - sigma**2/2) * T)，波动率高时差距很大：多数路径收在均值之下，均值被少数极高的路径拉高。这种右偏的对数正态分布，就是长期股价的形状。
import numpy as np
rng = np.random.default_rng(5)
s0, mu, sigma, T, steps, paths = 100.0, 0.08, 0.4, 1.0, 252, 20000
dt = T / steps
z = rng.standard_normal((paths, steps))
final = s0 * np.exp(np.sum((mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z, axis=1))
print(round(float(final.mean()), 1), round(s0 * np.exp(mu * T), 1)) # close to each other: the expected price
print(round(float(np.median(final)), 1), round(s0 * np.exp((mu - 0.5 * sigma ** 2) * T), 1)) # the median is lower
print(round(float(np.mean(final < s0)), 3)) # the share of paths that lose money估计 mu 和 sigma
对每隔 dt 年采样的价格，对数收益率 r = diff(log(prices)) 是均值为 (mu - sigma**2/2) * dt、标准差为 sigma * sqrt(dt) 的独立正态变量。所以 sigma = std(r) / sqrt(dt)，mu = mean(r) / dt + sigma**2 / 2。波动率几百个点就估得不错；漂移则出了名地噪声大，因为它的标准误只随总*时间*缩小，而不是随观测数量缩小。
import numpy as np
def estimate_gbm(prices, dt):
r = np.diff(np.log(prices))
sigma = r.std(ddof=1) / np.sqrt(dt)
return r.mean() / dt + 0.5 * sigma ** 2, sigma
rng = np.random.default_rng(8)
dt = 1 / 252
for years in (1, 10, 100):
n = int(years * 252)
log_moves = (0.10 - 0.5 * 0.25 ** 2) * dt + 0.25 * np.sqrt(dt) * rng.standard_normal(n)
prices = 100 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))
mu, sigma = estimate_gbm(prices, dt)
print(years, "years:", round(float(mu), 3), round(float(sigma), 3)) # true mu 0.10 and sigma 0.25- GBM：
S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t))；百分比变化是随机的，价格始终为正。 - sigma**2 / 2抵消了exp的凸性：期望价格是S0 * exp(mu * t)，而中位路径以mu - sigma**2 / 2增长。- 由对数收益估计：
sigma = std(r) / sqrt(dt)，mu = mean(r) / dt + sigma**2 / 2；波动率估得准，漂移估得差。
下一节：伊藤引理，解释 - sigma**2 / 2 从何而来的链式法则。
实验：评分内容
- 用
gbm_step走精确的一步25 分 - 用
gbm_path模拟价格路径25 分 - 用
gbm_expected对所有未来取平均25 分 - 用
estimate_gbm反推参数25 分