正向随机微分方程的经典数值格式模拟

前言
一：随机微分方程（SDE）的相关性质
二：方程的数值解格式
三：实验模拟
四：总结

前言

随机微分方程的发展分支可以说到现在非常广泛了，从19世纪的布朗运动发现到上世纪40年代的伊藤积分再到正向伊藤随机微分方程解的存在唯一性定理，后又跨越到上世纪90年代的倒向随机微分方程解的存在唯一性定理发现，最后到近现在的带跳、带马尔科夫、带延迟等一系列正倒向随机微方程解的存在唯一性定理等发现，可以说体系越发广泛，这也是顺应时代发展的必然趋势。

本篇文章将数值模拟经典下的正向随机微分方程数值格式的模拟，由于不是存粹数学研究，有些假设和定理和推导会简单概括，重点是计算机编程模拟，也算是学习随机微分方程数值法的入门介绍。

一：随机微分方程（SDE）的相关性质

简介

随机过程 XtX_tXt 被称作一个扩散过程，如果它满足下面的随机微分方程：

Xt=X0+∫t0ta(s,Xs)ds+∫t0tb(s,Xs)dWs,t∈[t0,T],X_t=X_0+\int_{t_0}^{t}a(s,X_s)ds+\int_{t_0}^{t}b(s,X_s)dW_s,t\in[t_0,T],Xt=X0+∫t0ta(s,Xs)ds+∫t0tb(s,Xs)dWs,t∈[t0,T], 其中 as=a(s,Xs)a_s=a(s,X_s)as=a(s,Xs) 和 bs=b(s,Xs)b_s=b(s,X_s)bs=b(s,Xs) 为 [0,T]×Rm[0,T]\times R^m[0,T]×Rm 上的可测函数， $W_t $为概率空间 (Ω,F,P)(\Omega,F,P)(Ω,F,P) 上的 d−d -d−维标准维纳过程。但为了此随机模型在实际中得以应用，保证上述扩散过程得以存在，我们将对此扩散过程给与一定条件，此SDESDESDE的解才会存在且唯一！

SDE解的存在唯一性定理

令 T>0T>0T>0 ，且对于扩散过程和漂移过程分别有

a(.,.):[0,T]×Rm→Rm,b(.,.):[0,T]×Rm→Rm×da(.,.):[0,T]\times R^m\rightarrow R^m,b(.,.):[0,T]\times R^m\rightarrow R^{m\times d}a(.,.):[0,T]×Rm→Rm,b(.,.):[0,T]×Rm→Rm×d 为可测函数，

满足 ∣a(s,x)+b(s,x)∣≤C(1+∣x∣),x∈Rm,s∈[t0,T]\left| a(s,x)+b(s,x) \right|\leq C(1+|x|),x\in R^m,s\in[t_0,T]∣a(s,x)+b(s,x)∣≤C(1+∣x∣),x∈Rm,s∈[t0,T],

其中 C,LC,LC,L 为有界常数，令 Φ\PhiΦ 为独立于由 Ws,s≥0W_s,s\geq0Ws,s≥0 产生的 σ\sigmaσ 代数 FFF ，使得

E[∣Z2∣]<∞E\left[ |Z^2| \right]<\inftyE[∣Z2∣]<∞。

则上述随机微分方程存在唯一的连续解，其中 X0∈ΦX_0\in\PhiX0∈Φ ,并且我们还知道

XtX_tXt 适应于 Φ\PhiΦ 和 Ws,s≤tW_s,s\leq tWs,s≤t 产生的域流 FtΦ且E[∫0T∣Xt∣2dt]<∞F_t^\Phi 且 E\left[ \int_{0}^{T}|X_t|^2dt \right]<\inftyFtΦ且E[∫0T∣Xt∣2dt]<∞。

结合伊藤公式的一般形式的随机算子表示

对于 m 维扩散过程和 d 维的维纳过程，且

φ(t,X)∈C2,1([0,T]×Rm)\varphi(t,X)\in C^{2,1}\left( [0,T]\times R^m \right)φ(t,X)∈C2,1([0,T]×Rm) ，定义算子

L0φ,Ljφ从[0,T]×Rm→RmL^0\varphi,L^j\varphi 从 [0,T]\times R^m\rightarrow R^mL0φ,Ljφ从[0,T]×Rm→Rm :

L0φ=∂φ∂t+∑k=1m∂φ∂xkak+12∑q=1d∑l,k=1q∂2φφxkφxlbl,qbk,q,LjφL^0\varphi=\frac{\partial\varphi}{\partial t}+\sum_{k=1}^{m}{\frac{\partial\varphi}{\partial x_k}a_k+\frac{1}{2}\sum_{q=1}^{d}{\sum_{l,k=1}^{q}{\frac{\partial^2\varphi}{\varphi x_k\varphi x_l}b_{l,q}b_{k,q}}}},L^j\varphiL0φ=∂t∂φ+∑k=1m∂xk∂φak+21∑q=1d∑l,k=1qφxkφxl∂2φbl,qbk,q,Ljφ

=∑k=1m∂φφxkbk,j=\sum_{k=1}^{m}{\frac{\partial \varphi}{\varphi x_k}b_{k,j}}=∑k=1mφxk∂φbk,j。

二：方程的数值解格式

伊藤展开

使用上述随机算子后，对于含任ddd 为维纳过程的 φ(t,X)∈C2,1([0,T]×Rm)\varphi(t,X)\in C^{2,1}\left( [0,T]\times R^m \right)φ(t,X)∈C2,1([0,T]×Rm) ，并且在时间 tnt_ntn 处展开，由伊藤公式我们有如下积分形式为：

φ(tn+1,Xtn+1)=φ(tn,Xtn)+∫tntn+1L0φ(s,Xs)ds+\varphi(t_{n+1},X_{t_{n+1}})=\varphi(t_{n},X_{t_n})+\int_{t_n}^{t_{n+1}}L^0\varphi(s,X_s)ds+φ(tn+1,Xtn+1)=φ(tn,Xtn)+∫tntn+1L0φ(s,Xs)ds+

∑j=1d∫tntn+1Ljφ(s,Xs)dWsj\sum_{j=1}^{d}{\int_{t_n}^{t_{n+1}}L^j\varphi(s,X_s)dW_s^j}∑j=1d∫tntn+1Ljφ(s,Xs)dWsj,

为了简化模型，我们对上述的只含一维维纳过程的 \varphi(.,.) 继续使用伊藤公式可得：

φ(tn+1,Xtn+1)=φ(tn,Xtn)+\varphi(t_{n+1},X_{t_{n+1}})=\varphi(t_{n},X_{t_n})+φ(tn+1,Xtn+1)=φ(tn,Xtn)+

∫tntn+1[L0φ(s,Xs)+∫tnsL0L0φ(τ,Xτ)+∫tnsL1L0φ(τ,Xτ)dWτ]ds+\int_{t_n}^{t_{n+1}}\left[ L^0\varphi(s,X_{s})+\int_{t_n}^{s}L^0L^0\varphi(\tau,X_{\tau})+\int_{t_n}^{s}L^1L^0\varphi(\tau,X_{\tau})dW_\tau\right]ds +∫tntn+1[L0φ(s,Xs)+∫tnsL0L0φ(τ,Xτ)+∫tnsL1L0φ(τ,Xτ)dWτ]ds+

∫tntn+1[L1φ(s,Xs)+∫tnsL0L1φ(τ,Xτ)dτ+∫tnsL1L1φ(τ,Xτ)dWτ]dWs\int_{t_n}^{t_{n+1}}\left[ L^1\varphi(s,X_s)+\int_{t_n}^{s}L^0L^1\varphi(\tau,X_\tau)d\tau +\int_{t_n}^{s}L^1L^1\varphi(\tau,X_\tau)dW_\tau\right]dW_s∫tntn+1[L1φ(s,Xs)+∫tnsL0L1φ(τ,Xτ)dτ+∫tnsL1L1φ(τ,Xτ)dWτ]dWs.

PS:只要上述 φ(.,.)\varphi(.,.)φ(.,.) 函数足够光滑，我们可以继续展开，可以得到更多的复杂高阶数值格式！

收敛阶的定义

给定时间分割 0=t0<t1<...<tN−1<tN=T0=t_0<t_1<...<t_{N-1}<t_N=T0=t0<t1<...<tN−1<tN=T ，假设 Δt=tn+1−tn\Delta t =t_{n+1}-t_nΔt=tn+1−tn ,

ΔWn=Wtn+1−Wtn(n=0,1,2,...,N)\Delta W_n=W_{t_{n+1}}-W_{t_n}(n=0,1,2,...,N)ΔWn=Wtn+1−Wtn(n=0,1,2,...,N),

且存在光滑函数 G∈Cp2(β+1)G\in C_p^{2(\beta+1)}G∈Cp2(β+1) 和正常数 CCC 满足不等式

∣E[G(Xt)−G(XN)]∣≤C(1+E[∣X0∣b])(Δt)β\left| E[G(X_t)-G(X^N)] \right|\leq C\left( 1+E[|X_0|^b] \right)(\Delta t)^\beta∣∣E[G(Xt)−G(XN)]∣∣≤C(1+E[∣X0∣b])(Δt)β ,

其中 β,b∈N\beta,b\in Nβ,b∈N ,则 XNX^NXN 以阶数 β\betaβ 弱收敛于 XTX_TXT. 当然有弱收敛就有强收敛，当我们不以整个轨迹误差取期望再绝对值，取而代之的是以每个点的绝对误差取期望后，就是我们的强收敛误差要求，高阶强收敛格式更加复杂，条件更加苛刻，阶数可以是0.5的倍数关系。

经典数值格式介绍

设初始条件 X0X_0X0 已知，对于 0≤n≤N−10\leq n\leq N-10≤n≤N−1 ，对于 mmm 维的伊藤过程的第 kkk 个分量，我们有如下格式：
1：欧拉格式（弱1.0阶，强0.5阶）：

Xkn+1=Xkn+akΔt+bkΔWnX^{n+1}_k=X_k^n+a_k\Delta t+b_k\Delta W_nXkn+1=Xkn+akΔt+bkΔWn

2：米尔斯坦格式（强弱都是1.0阶）：

Xkn+1=Xkn+akΔt+bkΔWn+12L1bk((ΔWn)2−Δt)X^{n+1}_k=X_k^n+a_k\Delta t+b_k\Delta W_n+\frac{1}{2}L^1b_k\left( (\Delta W_n)^2-\Delta t \right)Xkn+1=Xkn+akΔt+bkΔWn+21L1bk((ΔWn)2−Δt)

3：弱二阶格式：

Xkn+1=Xkn+akΔt+bkΔWn+12L1bk((ΔWn)2−Δt)+(bk+12Δt(L0bk+L1ak)ΔWn)+12L0ak(Δt)2X^{n+1}_k=X_k^n+a_k\Delta t+b_k\Delta W_n+\frac{1}{2}L^1b_k\left( (\Delta W_n)^2-\Delta t \right)+\left( b_k+\frac{1}{2}\Delta t(L^0b_k+L^1a_k)\Delta W_n \right)+\frac{1}{2}L^0a_k(\Delta t)^2Xkn+1=Xkn+akΔt+bkΔWn+21L1bk((ΔWn)2−Δt)+(bk+21Δt(L0bk+L1ak)ΔWn)+21L0ak(Δt)2

三：实验模拟

我们以几何布朗运动 Xt=aXtdt+bXtdWtX_t=aX_tdt+bX_tdW_tXt=aXtdt+bXtdWt ，只要简单的对 ln(Xt)ln(X_t)ln(Xt) 使用伊藤公式，就可得此模型真解为Xt=X0exp[(a−12b2)t+bWt]X_t=X_0exp[ (a-\frac{1}{2}b^2)t+bW_t ]Xt=X0exp[(a−21b2)t+bWt] ，所以在实验中能清楚地比较真解和数值解的误差结果，正向随机微分方程经典模型还有Ornstein−UhlenbeckOrnstein-UhlenbeckOrnstein−Uhlenbeck过程 Xt=a(μ−Xt)dt+bdWtX_t=a(\mu-X_t)dt+bdW_tXt=a(μ−Xt)dt+bdWt 等，也是存在真解的，但实验内容是跟几何布朗运动完全一致的。

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['font.sans-serif']=['SimHei']##中文乱码问题！
plt.rcParams['axes.unicode_minus']=False#横坐标负号显示问题！def schemes(a,b,names):logtimels = []logerrorls = []CRls = []msg = '''1：exp函数2：sin函数3：原函数的3次方4：原函数'''print(msg)##定义初始化内容choose = int(input('请选择g(x)函数：'))echoes = 100time_ls = [8,16,32,64,128]error_ls = []value_matx = np.zeros((echoes,len(time_ls)))Euler = 'y + a * dt * y + b * y * dwt_temp'Milstein = Euler + ' + 0.5 * b * b * y * (dwt_temp * dwt_temp - dt)'weak_2_order = Milstein + ' + 0.5 * a * a * y * dt * dt + 0.5 * dt * (b * a * y + a * b * y) * dwt_temp'schemes_ls = [Euler,Milstein,weak_2_order]for s,n in zip(schemes_ls,names):for e in range(echoes):for i in range(len(time_ls)):x0 = 1##SDE真解初始值y = 1##数值格式初始值T = 1N = time_ls[i]dt = 1 / Nyt = []dwt = np.random.normal(0,1,N) * np.sqrt(dt)Wt = np.cumsum(dwt)###真解WT = Wt[N - 1]vtemp = x0 * np.exp((a - 0.5 * b * b) * T + b * WT)xT = 0if choose == 1:xT = np.exp(vtemp)if choose == 2:xT = np.sin(vtemp)if choose == 3:xT = np.power(vtemp,3)if choose == 4:xT = vtemp##数值格式解for j in range(N):dwt_temp = dwt[j]y = eval(s)if choose == 1:yt.append(np.exp(y))if choose == 2:yt.append(np.sin(y))if choose == 3:yt.append(np.power(y,3))if choose == 4:yt.append(y)error_temp = yt[N - 1] - xT ##真解与数值解末端值作差error_ls.append(error_temp)value_matx[e,:] = error_lserror_ls = []error_fin = np.abs(sum(value_matx) / echoes)##弱收敛最终求期望log_time = np.log(1 / np.array(time_ls))log_error_fin = np.log(error_fin)CR = np.polyfit(log_time,log_error_fin,1)##一次多项式拟合print('\033[1;36m{0:*^80}\033[0m'.format('计算结果'))print('格式%s的相对误差为：%s'%(n,np.round(error_fin,6)))print('格式%s的弱收敛阶为：%s'%(n,round(CR[0],6)))logerrorls.append(np.round(log_error_fin,6))logtimels.append(log_time)CRls.append(round(CR[0],6))return logtimels,logerrorls,CRlsmu = 1##漂移系数
sigm = 0.05##扩散系数
scheme_names = ['Euler', 'Milstein', 'weak_2_order']
R = schemes(mu,sigm,scheme_names)

图1

图2

图3

marker_ls = ['-*','-o','-^']
def make_figure():plt.figure(figsize=(8, 6))for i,j,k,l,m in zip(R[0],R[1],R[2],scheme_names,marker_ls):plt.plot(i,j,m,label=l)plt.plot(np.array([-4,-3]),np.array([-4,-3]) + 1,'--',label='slope=1')plt.plot(np.array([-4, -3]), 2 * np.array([-4, -3]) - 1.5,'--',label='slope=2')plt.xlabel('stepsize logN')plt.ylabel('log(|E(Y_T-X^N)|)')plt.legend()plt.show()
make_figure()

图4

四：总结

经典的正向随机微分方程数值格式和收敛阶的模拟大致如上，但最开始也提到过，由正向衍生出来的各种方程在如今是非常多的。与正向相对性的倒向随机微分方程也是一个大类研究方向，今后如果有机会，作者还会发布些关于倒向随机微分方程的数值模拟，感兴趣的朋友希望点赞关注加收藏啦！