05 - 随机微分方程基础¶

import numpy as np
import matplotlib.pyplot as plt

def simulate_wiener_process(T=1.0, N=1000, num_paths=5):
    """
    模拟维纳过程

    参数:
        T: 终止时间
        N: 时间步数
        num_paths: 模拟路径数

    返回:
        t: 时间点
        W: 维纳过程路径
    """
    dt = T / N
    t = np.linspace(0, T, N+1)

    # 生成增量
    dW = np.random.normal(0, np.sqrt(dt), (num_paths, N))

    # 累积求和
    W = np.zeros((num_paths, N+1))
    W[:, 1:] = np.cumsum(dW, axis=1)

    return t, W

# 模拟
np.random.seed(42)
t, W = simulate_wiener_process(T=1.0, N=1000, num_paths=5)

# 可视化
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
for i in range(5):
    plt.plot(t, W[i], alpha=0.7, label=f'Path {i+1}')
plt.xlabel('Time')
plt.ylabel('$W_t$')
plt.title('Wiener Process Sample Paths')
plt.legend()
plt.grid(True, alpha=0.3)

# 验证性质
plt.subplot(1, 2, 2)
# 计算多个路径的统计量
num_samples = 1000
_, W_many = simulate_wiener_process(T=1.0, N=100, num_paths=num_samples)

# 计算均值和方差
mean = np.mean(W_many, axis=0)
variance = np.var(W_many, axis=0)

plt.plot(t[::10], mean[::10], 'b-', label='Empirical Mean', linewidth=2)
plt.plot(t[::10], variance[::10], 'r-', label='Empirical Variance', linewidth=2)
plt.plot(t, np.zeros_like(t), 'b--', label='Theoretical Mean (0)', alpha=0.5)
plt.plot(t, t, 'r--', label='Theoretical Variance (t)', alpha=0.5)
plt.xlabel('Time')
plt.ylabel('Value')
plt.title('Statistical Properties')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('wiener_process.png', dpi=150)
plt.show()

print("维纳过程性质验证:")
print(f"终点均值: {np.mean(W[:, -1]):.4f} (理论: 0)")
print(f"终点方差: {np.var(W[:, -1]):.4f} (理论: 1.0)")

class ItoIntegrator:
    """
    Itô积分计算
    """
    def __init__(self, seed=None):
        if seed is not None:
            np.random.seed(seed)

    def integrate(self, f, T=1.0, N=10000, num_samples=1000):
        """
        计算Itô积分 ∫₀ᵀ f(t) dW_t

        参数:
            f: 被积函数 f(t)
            T: 积分上限
            N: 时间步数
            num_samples: 蒙特卡洛样本数

        返回:
            integral: 积分值的蒙特卡洛估计
            variance: 估计的方差
        """
        dt = T / N
        t = np.linspace(0, T, N+1)

        integrals = []

        for _ in range(num_samples):
            # 生成维纳过程
            dW = np.random.normal(0, np.sqrt(dt), N)

            # 计算Itô积分（左端点法则）
            f_values = f(t[:-1])  # 左端点
            integral = np.sum(f_values * dW)
            integrals.append(integral)

        return np.mean(integrals), np.var(integrals)

    def verify_ito_isometry(self, f, T=1.0, N=10000, num_samples=5000):
        """
        验证Itô等距
        """
        dt = T / N
        t = np.linspace(0, T, N+1)

        # 计算左边: E[(∫fdW)²]
        integrals = []
        for _ in range(num_samples):
            dW = np.random.normal(0, np.sqrt(dt), N)
            f_values = f(t[:-1])
            integral = np.sum(f_values * dW)
            integrals.append(integral**2)

        lhs = np.mean(integrals)

        # 计算右边: E[∫f²dt]
        f_squared = f(t[:-1])**2
        rhs = np.sum(f_squared) * dt

        return lhs, rhs

# 测试Itô积分
integrator = ItoIntegrator(seed=42)

# 测试函数: f(t) = 1
f1 = lambda t: np.ones_like(t)  # lambda匿名函数
mean1, var1 = integrator.integrate(f1, T=1.0, N=1000, num_samples=10000)
print(f"\nItô积分 ∫₀¹ 1 dW_t:")
print(f"均值: {mean1:.4f} (理论: 0)")
print(f"方差: {var1:.4f} (理论: 1.0)")

# 测试函数: f(t) = t
f2 = lambda t: t
mean2, var2 = integrator.integrate(f2, T=1.0, N=1000, num_samples=10000)
print(f"\nItô积分 ∫₀¹ t dW_t:")
print(f"均值: {mean2:.4f} (理论: 0)")
print(f"方差: {var2:.4f} (理论: 1/3 ≈ 0.333)")

# 验证Itô等距
lhs, rhs = integrator.verify_ito_isometry(f2, T=1.0, N=1000, num_samples=5000)
print(f"\nItô等距验证:")
print(f"左边 E[(∫fdW)²]: {lhs:.4f}")
print(f"右边 E[∫f²dt]: {rhs:.4f}")
print(f"相对误差: {abs(lhs - rhs) / rhs * 100:.2f}%")

def verify_ito_formula():
    """
    验证Itô公式: 对于f(X_t) = X_t², 其中dX_t = dW_t

    Itô公式给出: df = 2X_t dX_t + dt = 2W_t dW_t + dt

    因此: W_t² = t + 2∫₀ᵗ W_s dW_s
    """
    T = 1.0
    N = 10000
    dt = T / N
    num_samples = 1000

    results = []

    for _ in range(num_samples):
        # 生成维纳过程
        dW = np.random.normal(0, np.sqrt(dt), N)
        W = np.zeros(N+1)
        W[1:] = np.cumsum(dW)

        # 计算Itô积分 ∫W_s dW_s
        ito_integral = np.sum(W[:-1] * dW)

        # 验证: W_t² = t + 2∫W_s dW_s
        lhs = W[-1]**2  # [-1]负索引取最后元素
        rhs = T + 2 * ito_integral

        results.append(abs(lhs - rhs))

    mean_error = np.mean(results)
    print(f"Itô公式验证 (W_t² = t + 2∫W_s dW_s):")
    print(f"平均误差: {mean_error:.6f}")
    print(f"验证通过！" if mean_error < 0.1 else "误差较大")

verify_ito_formula()

class SDESimulator:
    """
    SDE模拟器
    """
    def __init__(self, seed=None):
        if seed is not None:
            np.random.seed(seed)

    def euler_maruyama(self, b, sigma, x0, T, N, num_paths=1):
        """
        Euler-Maruyama方法求解SDE

        参数:
            b: 漂移函数 b(t, x)
            sigma: 扩散函数 sigma(t, x)
            x0: 初始条件
            T: 终止时间
            N: 时间步数
            num_paths: 模拟路径数

        返回:
            t: 时间点
            X: 解路径
        """
        dt = T / N
        t = np.linspace(0, T, N+1)

        X = np.zeros((num_paths, N+1))
        X[:, 0] = x0

        for i in range(N):
            dW = np.random.normal(0, np.sqrt(dt), num_paths)
            X[:, i+1] = X[:, i] + b(t[i], X[:, i]) * dt + sigma(t[i], X[:, i]) * dW

        return t, X

    def milstein(self, b, sigma, sigma_derivative, x0, T, N, num_paths=1):
        """
        Milstein方法（更高精度）
        """
        dt = T / N
        t = np.linspace(0, T, N+1)

        X = np.zeros((num_paths, N+1))
        X[:, 0] = x0

        for i in range(N):
            dW = np.random.normal(0, np.sqrt(dt), num_paths)
            sigma_val = sigma(t[i], X[:, i])

            X[:, i+1] = (X[:, i] + b(t[i], X[:, i]) * dt +
                        sigma_val * dW +
                        0.5 * sigma_val * sigma_derivative(t[i], X[:, i]) * (dW**2 - dt))

        return t, X

# 测试：几何布朗运动
simulator = SDESimulator(seed=42)

# 参数
mu, sigma = 0.1, 0.3
x0, T, N = 1.0, 1.0, 1000

# 漂移和扩散函数
b = lambda t, x: mu * x
sigma_func = lambda t, x: sigma * x
sigma_derivative = lambda t, x: sigma

# Euler-Maruyama
t, X_em = simulator.euler_maruyama(b, sigma_func, x0, T, N, num_paths=100)

# Milstein
t, X_mil = simulator.milstein(b, sigma_func, sigma_derivative, x0, T, N, num_paths=100)

# 解析解
W_T = np.random.normal(0, np.sqrt(T), 100)
X_exact = x0 * np.exp((mu - 0.5 * sigma**2) * T + sigma * W_T)

print(f"\n几何布朗运动模拟:")
print(f"Euler-Maruyama终点均值: {np.mean(X_em[:, -1]):.4f}")
print(f"Milstein终点均值: {np.mean(X_mil[:, -1]):.4f}")
print(f"解析解终点均值: {np.mean(X_exact):.4f}")
print(f"理论值: {x0 * np.exp(mu * T):.4f}")