随机微分方程的数值仿真

随机微分方程(Stochastic differential equation, SDE)在很多领域有广泛的应用，例如物理、工程、金融等。

只有少数几类随机微分方程具有解析解，在绝大多数情形下，我们需要进行数值近似。

考虑一个一般性的一维随机微分方程：

其中 和 为已知函数， 为一个布朗运动。

最简单的仿真方法是Euler方法。

假如我们需要对于 在时间区间 内进行仿真。

我们可以首先离散时间，即将时间区间 离散成 份，每一份的长度为 。

考察 在第 份时间区间 内的变化：

当时间步长 较小的情况下， 和 在内变化不大，它们分别可以用左端点的值来近似，即：

把这两个近似代入以上积分，我们就有：

我们可以根据这个近似来对 进行仿真。

令近似的仿真值为

满足迭代关系：

注意到 是一列独立的正态分布随机变量，且：

其中 为一列独立的标准正态分布的随机变量。

所以算法如下：

• 设 。
• 依次对于 ：生成标准正态随机变量 ，然后设：

该方法的误差的弱收敛阶为 ，即对于光滑函数 :

其中 为不依赖于 的常数。

其强收敛阶要低一些，为 。即存在常数 使得：

通过Milstein方法可以把强收敛阶提升为 。其与Euler方法不同的是将 近似为：

因此：

通过伊藤公式，我们可以求得：

利用

我们就得到了Milstein方法的递推公式：

Milstein方法与Euler方法的不同点就是最后一项。

这两种方法的代码实现和结果如下所示。

# 利用Euler和Milstein方法对随机微分方程进行仿真
# dX_t = mu(t, X_t)dt + sigma(t, X_t) dW_t

# 以Black-Schole模型为例

import numpy as np

# 参数
r, sig, x0 = 0.02, 0.3, 1.0

def mu(t, x): # mu函数
    return r*x

def sigma(t, x): # sigma函数
    return sig*x

def sigma_x(t, x): # sigma函数关于x的偏导数
    return sig

def Euler(T, m, n): 
    # n: sample path的数目
    # m: 时间离散的步数
    # T: 仿真的终止时间
    
    h = T/m # 时间步长
    
    Z = np.random.normal(0.0, 1.0, (n,m)) # 标准正态分布变量
    
    X = np.zeros((n,m+1)) # 存储n条路径的矩阵
    
    X[:,0] = x0 # 初始化
    
    for i in range(m):
        # 迭代：矩阵化实现
        X[:,i+1] = X[:,i] + mu(i*h,X[:,i])*h \
        + sigma(i*h,X[:,i])*np.sqrt(h)*Z[:,i]
    
    return X

def Milstein(T, m, n): 
    # n: sample path的数目
    # m: 时间离散的步数
    # T: 仿真的终止时间
    
    h = T/m # 时间步长
    
    Z = np.random.normal(0.0, 1.0, (n,m)) # 标准正态分布变量
    
    X = np.zeros((n,m+1)) # 存储n条路径的矩阵
    
    X[:,0] = x0 # 初始化
    
    for i in range(m):
        # 迭代：矩阵化实现
        X[:,i+1] = X[:,i] + mu(i*h,X[:,i])*h \
        + sigma(i*h,X[:,i])*np.sqrt(h)*Z[:,i] \
        + 0.5*sigma_x(i*h,X[:,i])*sigma(i*h,X[:,i])*h*(Z[:,i]**2-1)
    
    return X

# 实验
T, m, n = 1.0, 1000, 10

# 用两种方法模拟
X_Euler = Euler(T, m, n)
X_Mil = Milstein(T, m, n)

import matplotlib.pyplot as plt

# 时间格点
t_grid = np.linspace(0,T,m+1)

# 展示Euler Scheme的结果
for i in range(n):
    plt.plot(t_grid, X_Euler[i,:])

plt.xlabel('$t$')
plt.ylabel('$X_t$')
plt.title('Euler Scheme for BS')
plt.savefig('Euler_Scheme_BS.jpg')
plt.show()

# 展示Milstein Scheme的结果
for i in range(n):
    plt.plot(t_grid, X_Mil[i,:])

plt.xlabel('$t$')
plt.ylabel('$X_t$')
plt.title('Milstein Scheme for BS')
plt.savefig('Milstein_Scheme_BS.jpg')
plt.show()