g

gqzhang

V1

2022/12/06阅读:23主题:默认主题

随机微分方程的数值仿真

随机微分方程(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.020.31.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.01.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.01.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.0100010

# 用两种方法模拟
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()

分类:

数学

标签:

数学

作者介绍

g
gqzhang
V1