01-随机微积分

随机微积分

预计学习时间：4-5 小时

难度：⭐⭐⭐⭐⭐（传统量化中最难的部分）

核心问题：当你的变量不是”确定的”，而是”随机的”，微积分还成立吗？需要怎么改？

---

概述

随机微积分是普通微积分的”随机版升级”。在金融里，股价、利率、汇率每时每刻都在随机变动，普通微积分的工具不够用了。

本节从最直观的”抛硬币”出发，一步步带你理解：

• 布朗运动：随机过程的”标准模型”
• 伊藤积分：对随机过程的”积分”
• 伊藤引理：随机微积分中的”链式法则”
• SDE：随机微分方程，描述资产价格的”运动方程”

---

一、从随机游走到布朗运动

1.1 随机游走：最简单的随机过程

白话解释：一个人站在原点，每一步随机向左或向右走 1 米，走 100 步后他在哪？

答案是：不确定。但我们可以知道他的”平均位置”和”波动范围”。

import numpy as np
import matplotlib.pyplot as plt
 
# ============================================================
# 实验：抛硬币随机游走
# 每一步：抛一枚均匀硬币，正面 +1，反面 -1
# ============================================================
 
np.random.seed(42)
n_steps = 500      # 走 500 步
n_paths = 10       # 模拟 10 条路径
 
fig, axes = plt.subplots(1, 2, figsize=(16, 5))
 
# 左图：少量路径，看清每条路径的形状
coin_flips = np.random.choice([1, -1], size=(n_paths, n_steps))
paths = np.cumsum(coin_flips, axis=1)
 
for i in range(n_paths):
    axes[0].plot(paths[i], alpha=0.7, lw=1.2)
 
axes[0].axhline(y=0, color='black', linestyle='--', alpha=0.3)
axes[0].set_title(f'随机游走：{n_paths} 条路径，每步 +/-1')
axes[0].set_xlabel('步数')
axes[0].set_ylabel('位置')
 
# 右图：大量路径的分布（终点位置的直方图）
n_paths_many = 50000
coin_flips_many = np.random.choice([1, -1], size=(n_paths_many, n_steps))
final_positions = np.sum(coin_flips_many, axis=1)
 
axes[1].hist(final_positions, bins=80, density=True, alpha=0.6, color='steelblue')
# 叠加正态分布拟合
mu, sigma = 0, np.sqrt(n_steps)  # 理论值：均值 0，标准差 sqrt(n)
x = np.linspace(-4*sigma, 4*sigma, 200)
from scipy.stats import norm
axes[1].plot(x, norm.pdf(x, mu, sigma), 'r-', lw=2, label=f'正态 N(0, {n_steps})')
axes[1].set_title(f'终点位置分布（{n_paths_many} 条路径）')
axes[1].legend()
 
plt.tight_layout()
plt.show()
 
print(f"终点位置的理论均值 = 0")
print(f"终点位置的理论标准差 = sqrt({n_steps}) = {np.sqrt(n_steps):.2f}")
print(f"终点位置的实际均值 = {final_positions.mean():.4f}")
print(f"终点位置的实际标准差 = {final_positions.std():.4f}")

关键发现：

1. 每条路径看起来像”锯齿线”，忽上忽下
2. 走 500 步后，终点位置近似服从正态分布
3. 标准差是 $\sqrt{n}$，不是 $n$——这意味着路径”扩散”的速度比线性慢

白话解读：走 $n$ 步后，你期望回到原点附近（均值=0），但偏离范围大约是 $\sqrt{n}$。步数翻 4 倍，偏离范围才翻 2 倍。

1.2 布朗运动 = 连续时间的随机游走

现在把”离散步”变成”连续时间”：不再是每秒走一步，而是每时每刻都在随机抖动。

布朗运动 $W_t$ 的三条核心性质：

性质	数学表述	白话解读
起点为零	$W_0=0$	一开始站在原点
增量独立	$W_{t+dt}-W_t$ 与过去无关	未来走势只取决于现在，不取决于历史
增量正态	$W_{t+dt}-W_t\sim N(0,dt)$	在极短时间 $dt$ 内，变化的幅度大约是 $\sqrt{dt}$

还有一个反直觉的性质：布朗运动的路径是连续的，但处处不可导。

白话解读：路径没有”断裂”，每一点都和前一点相连（连续）。但如果你试图在任何一点求”斜率”，你会发现斜率不存在，因为路径在无限放大后仍然是锯齿状的（不可导）。

1.3 用 Python 模拟布朗运动

import numpy as np
import matplotlib.pyplot as plt
 
# ============================================================
# 模拟标准布朗运动 W_t
# 离散近似：W(t+dt) = W(t) + sqrt(dt) * Z, Z ~ N(0,1)
# ============================================================
 
np.random.seed(42)
T = 1.0           # 总时间（比如 1 年）
n_steps = 10000   # 离散步数
dt = T / n_steps  # 每步时间间隔
n_paths = 1000    # 模拟 1000 条路径
 
# 生成所有路径：每一步加一个 N(0, dt) 的随机增量
increments = np.random.normal(0, np.sqrt(dt), size=(n_paths, n_steps))
brownian = np.cumsum(increments, axis=1)
 
# 在前面补上起点 W_0 = 0
brownian = np.hstack([np.zeros((n_paths, 1)), brownian])
time_points = np.linspace(0, T, n_steps + 1)
 
fig, axes = plt.subplots(2, 2, figsize=(16, 10))
 
# 1. 左上：展示几条典型路径
sample_paths = [0, 1, 2, 3, 4]
colors = ['#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f00']
for i, c in zip(sample_paths, colors):
    axes[0, 0].plot(time_points, brownian[i], color=c, alpha=0.8, lw=1)
axes[0, 0].set_title('布朗运动：5 条典型路径')
axes[0, 0].set_xlabel('时间 t')
axes[0, 0].set_ylabel('W_t')
 
# 2. 右上：所有路径的均值和标准差包络
mean_path = brownian.mean(axis=0)
std_path = brownian.std(axis=0)
axes[0, 1].plot(time_points, mean_path, 'k-', lw=2, label='均值')
axes[0, 1].fill_between(time_points, mean_path - 2*std_path,
                        mean_path + 2*std_path, alpha=0.2, color='blue',
                        label='2 倍标准差')
axes[0, 1].plot(time_points, 2*std_path, 'b--', lw=1, alpha=0.5)
axes[0, 1].plot(time_points, -2*std_path, 'b--', lw=1, alpha=0.5)
axes[0, 1].set_title(f'路径统计：均值 ≈ 0，标准差 ≈ √t')
axes[0, 1].legend()
 
# 3. 左下：某一时刻的分布
snapshot_time = 0.5  # t = 0.5
snapshot_idx = int(snapshot_time / dt)
snapshot_values = brownian[:, snapshot_idx]
 
axes[1, 0].hist(snapshot_values, bins=50, density=True, alpha=0.6, color='steelblue')
from scipy.stats import norm
x = np.linspace(snapshot_values.min(), snapshot_values.max(), 200)
axes[1, 0].plot(x, norm.pdf(x, 0, np.sqrt(snapshot_time)), 'r-', lw=2,
               label=f'N(0, {snapshot_time})')
axes[1, 0].set_title(f't={snapshot_time} 时 W_t 的分布')
axes[1, 0].legend()
 
# 4. 右下：路径数量 vs 标准差收敛
path_counts = [10, 50, 200, 1000]
stds_sampled = [brownian[:n, -1].std() for n in path_counts]
std_theoretical = np.sqrt(T)
 
axes[1, 1].bar(range(len(path_counts)), stds_sampled,
               color='steelblue', alpha=0.7, tick_label=[str(n) for n in path_counts])
axes[1, 1].axhline(y=std_theoretical, color='red', linestyle='--', lw=2,
                  label=f'理论值 sqrt({T}) = {std_theoretical:.2f}')
axes[1, 1].set_title('路径数量越多，样本标准差越接近理论值')
axes[1, 1].set_xlabel('路径数量')
axes[1, 1].set_ylabel('终点标准差')
axes[1, 1].legend()
 
plt.tight_layout()
plt.show()
 
print(f"理论：E[W_t] = 0, Var(W_t) = t")
print(f"t=1 时理论标准差 = {np.sqrt(T):.4f}")
print(f"t=1 时 {n_paths} 条路径的标准差 = {brownian[:, -1].std():.4f}")

白话解读：

• 布朗运动在任意时刻 $t$ 的值服从 $N(0,t)$：均值永远是 0，方差随时间线性增长
• 1000 条路径的样本标准差已经非常接近理论值 $\sqrt{t}$

---

二、伊藤积分

2.1 为什么普通积分不适用于随机过程？

回顾普通积分（黎曼积分）的核心思想：

把面积切成无数个细条，每个细条的宽度趋近于 0，高度由函数在该点的值决定。

公式：${\int }_0^{T}f(t)dt={\lim }_{n\to \infty }{\sum }_{i=0}^{n-1}f(t_i)(t_{i+1}-t_i)$

这里的关键假设是：$f(t_i)$ 在区间 $[t_i,t_{i+1}]$ 内是”确定”的。但对于布朗运动，这个假设不成立：

布朗运动的增量是随机的，你在区间的哪一点取值，结果不一样。

2.2 伊藤积分的白话解释

普通积分：把面积切成细条，每条的高度是确定的。

伊藤积分：把面积切成细条，每条的高度取决于你”从哪一端开始算”。

伊藤选择的是左端点（$f(t_i)$），而不是右端点或中点。这个选择有深刻的原因：

左端点意味着”你只能用截至目前已知的信息来做决策”，不能偷看未来。

这和金融交易完全一致：你只能基于当前已知信息做交易决策，不能预知下一刻的价格。

2.3 伊藤积分的数学定义

${\int }_0^{T}f(t)d{W}_t={\lim }_{n\to \infty }{\sum }_{i=0}^{n-1}f(t_i)[W_{t_{i+1}}-W_{t_i}]$

白话解读：在时刻 $t_i$，用已知的信息 $f(t_i)$，乘以下一步的布朗运动增量 $W_{t_{i+1}}-W_{t_i}$。累加所有这些乘积，取极限，就是伊藤积分。

2.4 伊藤引理（随机微积分的”链式法则”）

这是整个随机微积分中最重要的一个结果。

普通微积分的链式法则：

如果 $Y=f(X)$，那么 $dY=f^{\prime }(X)\cdot dX$。

伊藤引理说：如果 $X$ 是布朗运动 $W_t$，链式法则需要多加一项！

$df(W_t)=f^{\prime }(W_t)d{W}_t+\frac{1}{2}{f}^{\prime \prime }(W_t)dt$

白话解读：对随机过程做链式法则时，除了原来的那一项，还要多加一个”二阶修正项” $\frac{1}{2}{f}^{\prime \prime }dt$。

为什么多出来？ 因为布朗运动的路径”太粗糙”了（处处不可导），一阶近似不够用，必须把二阶项也考虑进去。

2.5 伊藤引理的乘法表（关键工具）

以下规则是推导所有随机微积分公式的基础：

乘积	结果	白话解读
$dt\cdot dt$	$=0$	无穷小时间的平方是”更无穷小”，可以忽略
$dt\cdot dW$	$=0$	确定性变化和随机变化相乘，量级太小
$dW\cdot dW$	$=dt$	随机变化的平方等于时间变化（核心！）
$dW\cdot dW\cdot dW$	$=0$	三次方以上都忽略

为什么 $dW\cdot dW=dt$？

直观理解：$dW$ 大约是 $\sqrt{dt}$ 的量级，所以 $dW\cdot dW\approx (\sqrt{dt}{)}^2=dt$。

严格地说：$E[(dW{)}^2]=dt$，且当 $dt\to 0$ 时，$(dW{)}^2-dt$ 的方差趋于 0，所以 $(dW{)}^2$ 在均方意义下收敛到 $dt$。

2.6 用伊藤引理做推导（完整示例）

问题：如果 $f(W_t)=W_t^2$，求 $d(W_t^2)$。

推导过程：

第一步：用伊藤引理，$f(x)=x^2$

$df(W_t)=f^{\prime }(W_t)d{W}_t+\frac{1}{2}{f}^{\prime \prime }(W_t)dt$

第二步：计算导数

$f^{\prime }(x)=2x,f^{\prime \prime }(x)=2$

第三步：代入

$d(W_t^2)=2W_{t}d{W}_t+\frac{1}{2}\cdot 2\cdot dt=2W_{t}d{W}_t+dt$

白话解读：布朗运动的平方变化，等于”两项布朗运动的乘积”加上”纯时间项”。

对比普通微积分：如果 $W_t$ 是普通函数，链式法则只会给出 $d(W_t^2)=2W_{t}d{W}_t$。多出来的 $dt$ 项就是因为路径太粗糙导致的修正。

再推导一个：$f(W_t)=e^{W_t}$，求 $d(e^{W_t})$。

第一步：$f(x)=e^x$，$f^{\prime }(x)=e^x$，$f^{\prime \prime }(x)=e^x$

第二步：代入伊藤引理

$d(e^{W_t})=e^{W_t}d{W}_t+\frac{1}{2}{e}^{W_t}dt=e^{W_t}\left(d{W}_t+\frac{1}{2}dt\right)$

白话解读：指数布朗运动的变化，等于自身乘以”布朗增量加半个时间增量”。

2.7 多维伊藤引理

如果有两个相关的随机过程 $X_t$ 和 $Y_t$：

$dX={\mu }_{X}dt+{\sigma }_{X}d{W}_1$ $dY={\mu }_{Y}dt+{\sigma }_{Y}d{W}_2$

其中 $d{W}_1\cdot d{W}_2=\rho dt$（$\rho$ 是两个布朗运动的相关系数）。

对于 $f(t,X,Y)$：

$df=\frac{\partial f}{\partial t}dt+\frac{\partial f}{\partial X}dX+\frac{\partial f}{\partial Y}dY+\frac{1}{2}\left(\frac{{\partial }^{2}f}{\partial X^2}{\sigma }_X^2+\frac{{\partial }^{2}f}{\partial Y^2}{\sigma }_Y^2+2\frac{{\partial }^{2}f}{\partial X\partial Y}{\sigma }_X{\sigma }_Y\rho \right)dt$

白话解读：和一维类似，但多了交叉项。两个随机过程的相关性 $\rho$ 会影响修正项的大小。

乘法表扩展：

乘积	结果
$d{W}_1\cdot d{W}_1$	$dt$
$d{W}_2\cdot d{W}_2$	$dt$
$d{W}_1\cdot d{W}_2$	$\rho dt$
$d{W}_i\cdot dt$	$0$

---

三、随机微分方程（SDE）

3.1 SDE 是什么？

白话解释：普通微分方程描述”确定性变量”的变化规则，SDE 描述”随机变量”的变化规则。

普通微分方程：$\frac{dX}{dt}=\mu X$（人口增长、复利计算）

随机微分方程：$dX=\mu Xdt+\sigma XdW$（股价变化）

SDE 的一般形式：

$d{X}_t=\mu (X_t,t)dt+\sigma (X_t,t)d{W}_t$

项	含义	白话解读
$\mu (X_t,t)$	漂移项 (drift)	确定性的”趋势”方向，类似风速
$\sigma (X_t,t)$	扩散项 (diffusion)	随机波动的”强度”，类似颠簸程度
$d{W}_t$	布朗运动增量	随机的”推动力”，不可预测

3.2 几何布朗运动（GBM）：股票价格的标准模型

GBM 是最常用的股票价格模型：

$d{S}_t=\mu S_{t}dt+\sigma S_{t}d{W}_t$

每一项的白话解释：

• $\mu$：股票的年化预期收益率（漂移率），比如 0.08 表示年化期望收益 8%
• $\sigma$：股票的年化波动率，比如 0.25 表示年化波动 25%
• $d{W}_t$：随机的”冲击”，不可预测
• $S_{t}dt$：收益和波动都与当前股价成正比（涨 10% 和跌 10% 与股价高低无关）

为什么叫”几何”布朗运动？ 因为它的解是对数正态的：$\ln (S_t/S_0)$ 服从正态分布，等价于 $S_t$ 服从对数正态分布。

3.3 GBM 的解析解推导

要求解 $dS=\mu Sdt+\sigma SdW$，令 $f(S)=\ln S$。

用伊藤引理：$f^{\prime }(S)=1/S$，$f^{\prime \prime }(S)=-1/S^2$

$d(\ln S)=\frac{1}{S}dS+\frac{1}{2}\left(-\frac{1}{S^2}\right)(dS{)}^2$

把 $dS=\mu Sdt+\sigma SdW$ 代入：

$(dS{)}^2={\sigma }^2{S}^{2}dt$（用乘法表：$dW\cdot dW=dt$，其余高阶项为 0）

$d(\ln S)=\frac{1}{S}(\mu Sdt+\sigma SdW)-\frac{1}{2S^2}{\sigma }^2{S}^{2}dt$ $d(\ln S)=\mu dt+\sigma dW-\frac{{\sigma }^2}{2}dt$ $d(\ln S)=\left(\mu -\frac{{\sigma }^2}{2}\right)dt+\sigma dW$

从 0 到 T 积分：

$\ln S_T-\ln S_0=\left(\mu -\frac{{\sigma }^2}{2}\right)T+\sigma W_T$

$S_T=S_0\exp \left[\left(\mu -\frac{{\sigma }^2}{2}\right)T+\sigma W_T\right]$

白话解读：

• $\mu -\frac{{\sigma }^2}{2}$：对数收益率的均值（注意多减了 ${\sigma }^2/2$！这是伊藤修正）
• $\sigma W_T$：对数收益率的随机部分，服从 $N(0,{\sigma }^{2}T)$
• 股价永远为正（指数函数的值域），这符合现实

为什么多减 ${\sigma }^2/2$？ 因为”波动会吃掉收益”。直观理解：跌 50% 需要涨 100% 才能回本，波动率越大，这种不对称效应越强，长期来看期望收益会打折扣。

3.4 用 Python 模拟 GBM

import numpy as np
import matplotlib.pyplot as plt
 
# ============================================================
# 模拟几何布朗运动（GBM）股价路径
# S_T = S_0 * exp((mu - sigma^2/2)*T + sigma*W_T)
# ============================================================
 
np.random.seed(42)
 
# ---- 参数设置 ----
S0 = 100.0       # 初始股价（元）
mu = 0.08        # 年化预期收益率 8%
sigma = 0.25     # 年化波动率 25%
T = 1.0          # 模拟 1 年
n_steps = 252    # 约 252 个交易日
dt = T / n_steps
n_paths = 5      # 先画 5 条路径看细节
 
# ---- 模拟路径 ----
time_points = np.linspace(0, T, n_steps + 1)
Z = np.random.normal(0, 1, size=(n_paths, n_steps))  # 标准正态增量
# W_T = sum(sqrt(dt) * Z_i)
W = np.cumsum(np.sqrt(dt) * Z, axis=1)
W = np.hstack([np.zeros((n_paths, 1)), W])  # 补上 W_0 = 0
 
# S_t = S_0 * exp((mu - sigma^2/2) * t + sigma * W_t)
drift = (mu - sigma**2 / 2) * time_points
S = S0 * np.exp(drift[None, :] + sigma * W)
 
# ---- 绘图 ----
fig, axes = plt.subplots(2, 2, figsize=(16, 10))
 
# 1. 左上：5 条典型 GBM 路径
colors = ['#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f00']
for i in range(n_paths):
    axes[0, 0].plot(time_points, S[i], color=colors[i], alpha=0.8, lw=1.5)
 
# 画理论期望值 E[S_t] = S_0 * exp(mu * t)
expected = S0 * np.exp(mu * time_points)
axes[0, 0].plot(time_points, expected, 'k--', lw=2, label=f'E[S_t] = S0*exp(μt)')
axes[0, 0].set_title(f'GBM：S0={S0}, μ={mu}, σ={sigma}')
axes[0, 0].set_xlabel('时间（年）')
axes[0, 0].set_ylabel('股价（元）')
axes[0, 0].legend()
 
# 2. 右上：不同波动率的影响
sigmas = [0.10, 0.25, 0.50, 0.80]
for sig in sigmas:
    Z_sigma = np.random.normal(0, 1, size=(n_paths, n_steps))
    W_sigma = np.cumsum(np.sqrt(dt) * Z_sigma, axis=1)
    W_sigma = np.hstack([np.zeros((n_paths, 1)), W_sigma])
    drift_sigma = (mu - sig**2 / 2) * time_points
    S_sigma = S0 * np.exp(drift_sigma[None, :] + sig * W_sigma)
    # 只画均值路径
    axes[0, 1].plot(time_points, S_sigma.mean(axis=0), lw=2, label=f'σ={sig}')
 
axes[0, 1].plot(time_points, expected, 'k--', lw=2, label='期望值')
axes[0, 1].set_title('不同波动率下的 GBM 路径（均值）')
axes[0, 1].legend()
 
# 3. 左下：1000 条路径的终点分布
n_paths_large = 10000
Z_large = np.random.normal(0, 1, size=(n_paths_large, n_steps))
W_large = np.cumsum(np.sqrt(dt) * Z_large, axis=1)
final_S = S0 * np.exp((mu - sigma**2/2)*T + sigma*W_large[:, -1])
 
axes[1, 0].hist(final_S, bins=60, density=True, alpha=0.6, color='steelblue')
axes[1, 0].axvline(x=S0*np.exp(mu*T), color='red', linestyle='--', lw=2,
                  label=f'E[S_T] = {S0*np.exp(mu*T):.1f}')
axes[1, 0].axvline(x=S0, color='gray', linestyle=':', lw=2, label=f'S_0 = {S0}')
axes[1, 0].set_title(f'1 年后股价分布（{n_paths_large} 条路径）')
axes[1, 0].set_xlabel('股价（元）')
axes[1, 0].legend()
 
# 4. 右下：对数收益率的正态性检验
log_returns = np.log(final_S / S0)
axes[1, 1].hist(log_returns, bins=60, density=True, alpha=0.6, color='darkorange')
 
from scipy.stats import norm
x_lr = np.linspace(log_returns.min(), log_returns.max(), 200)
axes[1, 1].plot(x_lr, norm.pdf(x_lr, (mu-sigma**2/2)*T, sigma*np.sqrt(T)),
               'r-', lw=2, label='理论 N(μ-σ²/2, σ²T)')
axes[1, 1].set_title('对数收益率的正态性验证')
axes[1, 1].set_xlabel('ln(S_T / S_0)')
axes[1, 1].legend()
 
plt.tight_layout()
plt.show()
 
print("=== GBM 模拟统计 ===")
print(f"初始价格 S_0 = {S0}")
print(f"理论 E[S_T] = S_0 * exp(μT) = {S0 * np.exp(mu*T):.2f}")
print(f"模拟 E[S_T] = {final_S.mean():.2f}")
print(f"理论中位数 = S_0 * exp((μ-σ²/2)T) = {S0 * np.exp((mu-sigma**2/2)*T):.2f}")
print(f"模拟中位数 = {np.median(final_S):.2f}")
print(f"注意：期望 > 中位数，这就是 σ²/2 修正的效果")

3.5 Ornstein-Uhlenbeck 过程：均值回归模型

不是所有随机过程都像股价那样”越走越远”。有些变量会围绕一个均值来回波动——比如利率、价差、温度。

OU 过程的 SDE：

$d{X}_t=\kappa (\theta -X_t)dt+\sigma d{W}_t$

每个参数的白话解释：

参数	含义	白话解读
$\theta$	长期均值	像一根”弹簧”的平衡位置
$\kappa$	回归速度	弹簧的”硬度”，$\kappa$ 越大回弹越快
$\sigma$	波动率	随机扰动的强度
$(\theta -X_t)$	偏离程度	离平衡位置越远，拉力越大

白话解读：OU 过程像一个被弹簧拉住的粒子。偏离均值越远，被拉回来的力量越大；同时在随机力量的扰动下不断振荡。

import numpy as np
import matplotlib.pyplot as plt
 
# ============================================================
# 模拟 Ornstein-Uhlenbeck（OU）过程
# dX = kappa * (theta - X) * dt + sigma * dW
# 离散近似：X(t+dt) = X(t) + kappa*(theta-X(t))*dt + sigma*sqrt(dt)*Z
# ============================================================
 
np.random.seed(42)
 
T = 5.0          # 模拟 5 年
n_steps = 5000   # 离散步数
dt = T / n_steps
time_points = np.linspace(0, T, n_steps + 1)
 
# ---- 参数设置 ----
kappa = 3.0      # 回归速度：约 3 意味着半衰期 ≈ ln2/3 ≈ 0.23 年
theta = 0.0      # 长期均值：价差的均衡位置
sigma = 0.5      # 波动率
 
# ---- 模拟多条路径 ----
n_paths = 5
X = np.zeros((n_paths, n_steps + 1))
X[:, 0] = 3.0    # 起始点偏离均值较远
 
for t in range(n_steps):
    Z = np.random.normal(0, 1, n_paths)
    X[:, t+1] = X[:, t] + kappa * (theta - X[:, t]) * dt + sigma * np.sqrt(dt) * Z
 
# ---- 绘图 ----
fig, axes = plt.subplots(1, 2, figsize=(16, 5))
 
# 左图：路径
colors = ['#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f00']
for i in range(n_paths):
    axes[0].plot(time_points, X[i], color=colors[i], alpha=0.7, lw=1)
axes[0].axhline(y=theta, color='black', linestyle='--', lw=2, label=f'长期均值 θ={theta}')
axes[0].set_title(f'OU 过程：κ={kappa}, θ={theta}, σ={sigma}')
axes[0].set_xlabel('时间（年）')
axes[0].set_ylabel('X_t')
axes[0].legend()
 
# 右图：不同回归速度的对比
kappas = [0.5, 2.0, 5.0, 15.0]
for kap in kappas:
    X_k = np.zeros(n_steps + 1)
    X_k[0] = 3.0
    for t in range(n_steps):
        Z = np.random.normal(0, 1)
        X_k[t+1] = X_k[t] + kap * (theta - X_k[t]) * dt + sigma * np.sqrt(dt) * Z
    half_life = np.log(2) / kap
    axes[1].plot(time_points, X_k, lw=1.5, label=f'κ={kap} (半衰期={half_life:.2f}年)')
 
axes[1].axhline(y=theta, color='black', linestyle='--', lw=2)
axes[1].set_title('不同回归速度的 OU 过程')
axes[1].set_xlabel('时间（年）')
axes[1].legend()
 
plt.tight_layout()
plt.show()
 
print("=== OU 过程统计 ===")
print(f"半衰期 = ln(2)/κ = {np.log(2)/kappa:.4f} 年 ≈ {np.log(2)/kappa*252:.1f} 个交易日")
print(f"稳态分布：X ~ N(θ, σ²/(2κ)) = N({theta}, {sigma**2/(2*kappa):.4f})")

---

四、进阶概念

4.1 Girsanov 定理：测度变换

白话解释：想象你戴着不同颜色的眼镜看同一个随机过程。

• P 测度（物理测度）：戴上”现实”眼镜，看到的布朗运动有向上或向下的趋势
• Q 测度（风险中性测度）：戴上”假设”眼镜，看到的世界里所有资产的期望收益都等于无风险利率

Girsanov 定理说的是：通过适当的”概率调整”（数学上叫 Radon-Nikodym 导数），你可以在 P 测度和 Q 测度之间切换。

数学表述（直觉版）：

$d{W}_t^P=d{W}_t^Q+\lambda dt$

其中 $\lambda$ 是”风险溢价”——从 Q 测度的视角看，P 测度下的布朗运动只是多了一个漂移项。

金融含义：定价时，我们不需要知道投资者对风险的厌恶程度（不需要知道 $\mu$），只需要把视角切换到 Q 测度，然后折现即可。这是风险中性定价的数学基础。

4.2 Feynman-Kac 公式：PDE 和 SDE 的桥梁

白话解释：这是一个”翻译器”——把一个随机微分方程的问题翻译成偏微分方程，或者反过来。

核心思想：

给定 SDE：$dX=\mu (X,t)dt+\sigma (X,t)dW$

定义 $u(x,t)=E^{x,t}[g(X_T)]$（从 $(x,t)$ 出发的条件期望）

那么 $u(x,t)$ 满足 PDE：

$\frac{\partial u}{\partial t}+\mu \frac{\partial u}{\partial x}+\frac{1}{2}{\sigma }^2\frac{{\partial }^{2}u}{\partial x^2}=0$

终端条件：$u(x,T)=g(x)$

金融含义：期权定价可以用两种方式做——蒙特卡洛模拟 SDE，或者数值求解 PDE。Feynman-Kac 告诉我们这两种方法等价。

4.3 停时：什么时候该停止

白话解释：停时是一个”自己决定什么时候停下”的时间点，但决定只能基于已知信息（不能偷看未来）。

正式定义：$\tau$ 是停时，如果对于所有 $t$，事件 $\{\tau \le t\}$ 只依赖于截至 $t$ 的信息。

金融应用：

• 美式期权：持有人可以在到期前的任何时间行权，最优行权时间就是一个停时
• 止损策略：设定一个止损线，触发就卖出，这也是停时
• 最优交易时机：在均值回归策略中，什么时候入场、什么时候出场

可选停时定理（重要结论）：

如果 $W_t$ 是布朗运动，$\tau$ 是满足某些条件的停时，那么：

$E[W_{\tau }]=0$

白话解读：无论你怎么选择停止时机（只要你不能预知未来），布朗运动的期望值永远是 0。这意味着你不可能通过”择时”来战胜公平的随机游走。

---

五、实战应用

5.1 用 GBM 模拟股价路径并分析

import numpy as np
import matplotlib.pyplot as plt
 
# ============================================================
# 实战：模拟 A 股风格的 GBM 路径
# ============================================================
 
np.random.seed(42)
 
# A 股典型参数
S0 = 10.0       # 起始价格 10 元
mu = 0.06       # 年化预期收益 6%
sigma = 0.30    # 年化波动率 30%（A股波动较大）
T = 2.0         # 模拟 2 年（约 500 个交易日）
n_steps = 500
dt = T / n_steps
n_paths = 1000
 
# 模拟
Z = np.random.normal(0, 1, size=(n_paths, n_steps))
log_returns = (mu - sigma**2/2) * dt + sigma * np.sqrt(dt) * Z
log_prices = np.cumsum(log_returns, axis=1)
log_prices = np.hstack([np.zeros((n_paths, 1)), log_prices])  # ln(S_0/S_0) = 0
prices = S0 * np.exp(log_prices)
 
time_points = np.linspace(0, T, n_steps + 1)
 
# ---- 分析 ----
print("=== GBM 股价模拟分析 ===")
print(f"参数：S0={S0}, μ={mu}, σ={sigma}, T={T}年")
print(f"模拟路径数：{n_paths}")
print()
 
# 终点统计
final_prices = prices[:, -1]
print(f"2 年后股价统计：")
print(f"  期望值（理论）= {S0 * np.exp(mu*T):.2f}")
print(f"  期望值（模拟）= {final_prices.mean():.2f}")
print(f"  中位数（理论）= {S0 * np.exp((mu-sigma**2/2)*T):.2f}")
print(f"  中位数（模拟）= {np.median(final_prices):.2f}")
print(f"  5% 分位数 = {np.percentile(final_prices, 5):.2f}")
print(f"  95% 分位数 = {np.percentile(final_prices, 95):.2f}")
print(f"  最大值 = {final_prices.max():.2f}")
print(f"  最小值 = {final_prices.min():.2f}")
 
# 亏损概率
loss_prob = (final_prices < S0).mean()
print(f"  亏损概率 = {loss_prob:.1%}")
print()
 
# 年化对数收益的均值和标准差
total_log_return = np.log(final_prices / S0)
print(f"2 年对数收益：均值 = {total_log_return.mean():.4f}, 标准差 = {total_log_return.std():.4f}")
print(f"理论：均值 = (μ-σ²/2)*T = {(mu-sigma**2/2)*T:.4f}, 标准差 = σ√T = {sigma*np.sqrt(T):.4f}")
 
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
 
# 几条典型路径
for i in range(5):
    axes[0].plot(time_points, prices[i], alpha=0.7, lw=1)
axes[0].plot(time_points, S0*np.exp(mu*time_points), 'k--', lw=2, label='期望值')
axes[0].set_title('A 股风格 GBM 路径')
axes[0].set_xlabel('时间（年）')
axes[0].set_ylabel('股价（元）')
axes[0].legend()
 
# 终点分布
axes[1].hist(final_prices, bins=50, density=True, alpha=0.6, color='steelblue')
axes[1].axvline(x=S0, color='red', linestyle='--', label=f'S_0={S0}')
axes[1].set_title('2 年后股价分布')
axes[1].set_xlabel('股价（元）')
axes[1].legend()
 
# 对数收益率正态性
axes[2].hist(total_log_return, bins=50, density=True, alpha=0.6, color='darkorange')
from scipy.stats import norm
x = np.linspace(total_log_return.min(), total_log_return.max(), 200)
axes[2].plot(x, norm.pdf(x, (mu-sigma**2/2)*T, sigma*np.sqrt(T)), 'r-', lw=2,
            label='理论正态分布')
axes[2].set_title('对数收益率分布验证')
axes[2].set_xlabel('ln(S_T/S_0)')
axes[2].legend()
 
plt.tight_layout()
plt.show()

5.2 用 OU 过程模拟均值回归价差

import numpy as np
import matplotlib.pyplot as plt
 
# ============================================================
# 实战：用 OU 过程模拟配对交易的价差
# ============================================================
 
np.random.seed(42)
 
# 配对交易场景：股票 A 和股票 B 的价差应该围绕 0 波动
T = 2.0          # 2 年
n_steps = 504    # 约 2 年的交易日
dt = T / n_steps
time_points = np.linspace(0, T, n_steps + 1)
 
# OU 参数
kappa = 5.0      # 快速回归（半衰期 ≈ 50 个交易日）
theta = 0.0      # 价差的长期均值 = 0
sigma = 0.3      # 价差的波动率
n_paths = 3
 
# 模拟价差路径
spread = np.zeros((n_paths, n_steps + 1))
spread[:, 0] = 0.0
 
for t in range(n_steps):
    Z = np.random.normal(0, 1, n_paths)
    spread[:, t+1] = spread[:, t] + kappa * (theta - spread[:, t]) * dt + sigma * np.sqrt(dt) * Z
 
# ---- 简单的交易信号：价差偏离超过 1 倍标准差时入场 ----
std_long_term = sigma / np.sqrt(2 * kappa)  # 稳态标准差
entry_threshold = 1.0 * std_long_term
exit_threshold = 0.2 * std_long_term
 
# 对第一条路径做信号分析
s = spread[0]
positions = np.zeros(n_steps + 1)  # +1 做多价差，-1 做空价差，0 空仓
pnl = np.zeros(n_steps + 1)
 
for t in range(1, n_steps + 1):
    if positions[t-1] == 0:
        # 空仓状态：价差低于阈值做空（预期回归），高于阈值做多
        if s[t] < -entry_threshold:
            positions[t] = 1   # 做多价差（价差太低，预期回升）
        elif s[t] > entry_threshold:
            positions[t] = -1  # 做空价差（价差太高，预期回落）
        else:
            positions[t] = 0
    else:
        # 持仓状态：价差回到 0 附近时平仓
        if abs(s[t]) < exit_threshold:
            positions[t] = 0
        else:
            positions[t] = positions[t-1]
 
    # 计算 PnL（假设价差变动 1 单位 = 盈利 1 单位）
    pnl[t] = pnl[t-1] + positions[t] * (s[t] - s[t-1])
 
# ---- 绘图 ----
fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True,
                         gridspec_kw={'height_ratios': [2, 1, 1]})
 
# 价差路径
for i in range(n_paths):
    axes[0].plot(time_points, spread[i], alpha=0.5, lw=0.8)
axes[0].plot(time_points, spread[0], color='#e41a1c', lw=1.5, label='主分析路径')
axes[0].axhline(y=entry_threshold, color='green', linestyle='--', label=f'入场线 ±{entry_threshold:.2f}')
axes[0].axhline(y=-entry_threshold, color='green', linestyle='--')
axes[0].axhline(y=theta, color='black', linestyle='-', lw=1.5, label=f'均值 θ={theta}')
axes[0].set_title(f'OU 过程价差模拟：κ={kappa}, σ={sigma}, 半衰期≈{np.log(2)/kappa*252:.0f}天')
axes[0].legend(loc='upper right')
axes[0].set_ylabel('价差')
 
# 持仓
axes[1].step(time_points, positions[0], where='post', color='blue', lw=1.5)
axes[1].set_ylabel('持仓')
axes[1].set_yticks([-1, 0, 1])
axes[1].set_yticklabels(['做空', '空仓', '做多'])
 
# PnL
axes[2].plot(time_points, pnl[0], color='darkorange', lw=1.5)
axes[2].fill_between(time_points, 0, pnl[0], where=pnl[0]>=0, alpha=0.3, color='green')
axes[2].fill_between(time_points, 0, pnl[0], where=pnl[0]<0, alpha=0.3, color='red')
axes[2].set_ylabel('累计 PnL')
axes[2].set_xlabel('时间（年）')
 
plt.tight_layout()
plt.show()
 
print(f"=== 配对交易信号统计 ===")
print(f"稳态标准差 = σ/√(2κ) = {std_long_term:.4f}")
print(f"入场阈值 = {entry_threshold:.4f}")
print(f"出场阈值 = {exit_threshold:.4f}")
print(f"总交易次数 = {np.sum(np.diff(positions[0]) != 0) // 2}")
print(f"最终 PnL = {pnl[-1]:.4f}")
print(f"最大回撤 = {np.min(pnl - np.maximum.accumulate(pnl)):.4f}")

---

六、本节要点回顾

概念	核心要点	金融应用
布朗运动	连续、不可导、增量独立	股价随机波动的数学模型
$dW\cdot dW=dt$	随机变化的”平方”等于时间变化	所有 SDE 推导的基础
伊藤引理	随机版链式法则，多一个 $\frac{1}{2}{f}^{\prime \prime }dt$ 项	从股价模型推导期权 PDE
GBM	$dS=\mu Sdt+\sigma SdW$	标准股价模型
OU 过程	$dX=\kappa (\theta -X)dt+\sigma dW$	均值回归、利率、价差
Girsanov	P 测度 ↔ Q 测度切换	风险中性定价的数学基础
Feynman-Kac	SDE ↔ PDE 互译	蒙特卡洛 vs 有限差分
停时	基于已知信息的停止规则	美式期权、止损策略

下一步：掌握了随机微积分这个工具箱后，我们就可以进入 02-无套利定价理论.md，学习如何用这些工具给金融产品定价。