0.5(Euler-Maruyama), 1(Milstein), 和1.5 阶强Stochastic Differential Equation格式总结

本文的目的在于提供0.5，1和1.5阶强SDE数值格式的推导和内容，所有推导基于 $I t o - T a y l or$ 展开，由于国内外网站缺少关于强SDE数值阶的总结，笔者在此特作总结，为使用SDE数值格式的读者提供帮助，本文需要读者预先已经知道了有关Brown-motion的基本性质和随机微分的基本性质，否则格式推导会看不懂。如果想要实际应用直接看数值格式即可。

1. SDE强收敛阶判断定理

本文考虑如下SDE:
$dX=b(x)dt+\sigma(x)dW.$
使用数值SDE在时间 $t$ 下以 $x$ 为初值的近似解定义为 $\overline{X}_{t,x}(t+h)$ ,其中 $t+h\geq t$ .真解为: $X_{t,x}(t+h)$ .若满足如下定理:
定理1:
若存在 $p_1,p_2$ ,满足:
$|X_{t,x}(t+h)-\overline{X}_{t,x}(t+h)|\leq k(1+|x|^2)h^{p_1}.$
$|X_{t,x}(t+h)-\overline{X}_{t,x}(t+h)|^2]^{\frac{1}{2}}\leq k(1+|x|^2)^{\frac{1}{2}}h^{p_2}.$
且满足 $p_1\geq p_2+\frac{1}{2}$ , $p_2\geq\frac{1}{2}$ .则可认为强收敛阶为 $p_2-\frac{1}{2}$ .
该定理的证明笔者在此不做推导，感兴趣的读者自行阅读其他文献或教材，这一般都是有的，这里我们直接利用它进行以下的 $I t o - T a y l or$ 展开。

2. $I t o - T a y l or$ 展开

我们考虑:
$dx=b(x)dt+\sigma(x)dW.$
若有一个非常光滑的 $f (x)$ ,根据 $I t o$ 公式，会有:
$df(x)=[b(x)f_x(x)+\frac{1}{2} \sigma^2 f_{xx}(x)]dt+[\sigma(x)f_x(x)]dW.$
这两个微分算子可以被定义为:
$L^0=b(x)\frac{d}{dx}+\frac{1}{2}\sigma^2(x)\frac{d^2}{dx^2}$
$L^1=\sigma(x)\frac{d}{dx}$
则对两边取积分会有:
$f(x_t)=f(x_{t0})+\int_{t_0}^t L^0f(x_s)ds+\int_{t_0}^tL^1f(x_s)dW_s.$
注意到:
$\int_{t_0}^t L^0f(x_s)ds=L^0 f(x_{t_0})(t-t_0)+\int_{t_0}^t[L^0f(x_s)-L^0f(x_{t_0})].$
$\int_{t_0}^t L^1f(x_s)dWs=L^1 f(x_{t_0})(W_t-W_{t_0})+\int_{t_0}^t[L^1f(x_s)-L^1f(x_{t_0})]dWs.$
注意到:
$\int_{t_0}^t[L^0f(x_s)-L^0f(x_{t_0})]=\int_{t_0}^t\int_{t_0}^s d[L^0f(x_{s_1})]ds$
$\int_{t_0}^t[L^1f(x_s)-L^1f(x_{t_0})]=\int_{t_0}^t\int_{t_0}^s d[L^1f(x_{s_1})]dWs$
使用 $I t o$ 公式我们可以知道:
$d[L^0f(x)]=[L^0L^0f]dt+[L^1L^0f]dW.$
$d[L^1f(x)]=[L^0L^1f]dt+[L^1L^1f]dW.$
我们便可以得到:
$f(x_t)=f(x_{t0})+L^0 f(x_{t_0})(t-t_0)+L^1 f(x_{t_0})(W_t-W_{t_0})+\overline{R}$
其中:
$\overline{R}=\int_{t_0}^t\int_{t_0}^s L^0L^0f(s_1)ds_1ds+\int_{t_0}^t\int_{t_0}^s L^1L^0f(s_1)dWs_1ds+\int_{t_0}^t\int_{t_0}^s L^0L^1f(s_1)ds_1dWs+\int_{t_0}^t\int_{t_0}^s L^1L^1f(s_1)dWs_1dWs$
当然可以继续展开下去。
我们在这里只关心 $f (x) = x$ 的情况，因此这会导出本文的主题0.5(Euler-Maruyama), 1(Milstein), 和1.5 阶强Stochastic Differential Equation格式。

3、强格式1: 0.5阶强SDE (Euler-Maruyama)

令 $f (x) = x$ ,显然地,我们会得到: $L^0f=L^0x=b(x)$ , $L^1f=L^1x=\sigma(x)$ .
则根据：
$f(x_t)=f(x_{t0})+L^0 f(x_{t_0})(t-t_0)+L^1 f(x_{t_0})(W_t-W_{t_0})+\overline{R}$
代入 $f (x) = x$ 可得:
$x_t=x_{t_0}+b(x_{t_0})(t-t_0)+\sigma(x_{t_0})(W_t-W_{t_0})+\overline{R}$
此时 $\overline{R}$ 为:
$\overline{R}=\int_{t_0}^t\int_{t_0}^s L^0b(x_{s_1})ds_1ds+\int_{t_0}^t\int_{t_0}^s L^1b(x_{s_1})dWs_1ds+\int_{t_0}^t\int_{t_0}^s L^0\sigma(x_{s_1})ds_1dWs+\int_{t_0}^t\int_{t_0}^s L^1\sigma(x_{s_1})dWs_1dWs$
称下式为Euler-Maruyama格式:
$x_t=x_{t_0}+b(x_{t_0})(t-t_0)+\sigma(x_{t_0})(W_t-W_{t_0})$ 需要用到如下随机微分性质: $E[dWs^2]=ds$
根据定理和随机微分的性质可知余项 $\overline{R}$ 中满足: $p_1=2$ (具有 $d W s$ 会使得阶数为0), $p_2=1$ (一个 $d s$ 贡献一个阶,一个 $d W s$ 贡献0.5个阶), 满足定理格式，此时强收敛阶为 $p_2-0.5=0.5$ 阶,则该格式具有0.5阶强收敛性。

4、强格式2/3: 1阶强SDE (Milstein)和1.5阶强SDE.

注意到余项:
$\overline{R}=\int_{t_0}^t\int_{t_0}^s L^0b(x_{s_1})ds_1ds+\int_{t_0}^t\int_{t_0}^s L^1b(x_{s_1})dWs_1ds+\int_{t_0}^t\int_{t_0}^s L^0\sigma(x_{s_1})ds_1dWs+\int_{t_0}^t\int_{t_0}^s L^1\sigma(x_{s_1})dWs_1dWs$
我们想要提升强格式阶，那么显然此时若能够提升 $p_2$ 的阶就可以了，注意到限制了 $p_2$ 的阶为 $\overline{R}$ 中这一项:
$G=\int_{t_0}^t\int_{t_0}^s L^1\sigma(x_{s_1})dWs_1dWs$ 那么，我们需要对其进行再一次的展开，此时根据 $I t o$ 公式，我们可以将它展开为:
$G=L^1\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^s dWs_1dWs+\int_{t_0}^t\int_{t_0}^s\int _{t_0}^{s_1}L^0L^1\sigma(x_{s_2})d{s_2}dWs_1dWs+\int_{t_0}^t\int_{t_0}^s\int _{t_0}^{s_1}L^1L^1\sigma(x_{s_2})d{Ws_2}dWs_1dWs$
此时可以得到Milstein格式为:
$x_t=x_{t_0}+b(x_{t_0})(t-t_0)+\sigma(x_{t_0})(W_t-W_{t_0})+L^1\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^s dWs_1dWs$ 这即:
$x_t=x_{t_0}+b(x_{t_0})(t-t_0)+\sigma(x_{t_0})(W_t-W_{t_0})+\sigma^{'}(x_{t_0})\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^s dWs_1dWs$ 其中:
$\int_{t_0}^t\int_{t_0}^s dWs_1dWs=\frac{1}{2}(W_t-W_{t_0})^2-\frac{1}{2}(t-t_0)$ 这即可以得到最终格式为:
$x_t=x_{t_0}+b(x_{t_0})(t-t_0)+\sigma(x_{t_0})(W_t-W_{t_0})+\frac{1}{2}\sigma^{'}(x_{t_0})\sigma(x_{t_0})[(W_t-W_{t_0})^2-(t-t_0)]$
此时注意到 $p_1=2$ 不变，因为提升了 $p_2=\frac{3}{2}$ ，满足定理要求，因此该格式收敛阶为1.下面继续提升，注意到此时若想提升强收敛阶，很明显提升 $p_2$ 是不够的了，因为此时 $p_1$ 也需要被提升，但是很明显的是在 $\overline{R}$ 中有 $p_2=\frac{3}{2}$ 的项两个:
$\int_{t_0}^t\int_{t_0}^s L^0\sigma(x_{s_1})ds_1dWs$ 和
$\int_{t_0}^t\int_{t_0}^s L^1b(x_{s_1})dWs_1ds$ 则他们也需要被包含在内才可以提升该格式的收敛阶，因此在 $M i l s t e in$ 的基础上，我们继续延拓格式 $T$ 为:
$T = M i l s t e in + K$
$K=L^0\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^sds_1dWs+L^1b(x_{t_0})\int_{t_0}^t\int_{t_0}^s dWs_1ds+L^0b(x_{t_0})\int_{t_0}^t\int_{t_0}^s ds_1ds+(L^1)^2\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^s\int _{t_0}^{s_1}d{Ws_2}dWs_1dWs$ 我们仅考虑从 $n$ 到 $n + 1$ 会发生什么，令 $\triangle t=t_{n+1}-t_n$ , $\triangle W=W_{t_{n+1}}-W_{t_n}$
$(L^1)^2\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^s\int _{t_0}^{s_1}d{Ws_2}dWs_1dWs=(L^1)^2\sigma(x_{t_0})[\frac{1}{6}\triangle W^2-\frac{1}{2}\triangle t]\triangle W$ 这一点使用 $I t o$ 公式很容易证明。
另一方面令:
$\triangle Z=\int_{t_n}^{t_{n+1}}\int_{t_n}^s dWs_1ds=\int_{t_n}^{t_{n+1}}[t_{n+1}-s] dW_s$ 根据随机微分的性质,显然
$\triangle Z \sim N(0,\int_{t_n}^{t_{n+1}}[t_{n+1}-s]^2 ds)=N(0,\frac{1}{3} \triangle t^3)$ 注意到:
$E[\triangle Z \triangle W]=\int_{t_n}^{t_{n+1}}[t_{n+1}-s] ds=\frac{1}{2} \triangle t^2$ 因此:
$E[(\triangle Z-\frac{1}{2} \triangle W)\triangle W]=0$ 这样找到了独立的样本，且注意到此时根据随机微分的性质可得到:
$(\triangle Z-\frac{1}{2} \triangle W)～N(0,\frac{1}{12}\triangle t^3)$ 注意到:
$\int_{t_0}^t\int_{t_0}^sds_1dWs=\triangle t \triangle W-\triangle Z$ 那么显然地,我们可以根据如下方法来给出1.5阶格式,若我们有 $G_1$ , $G_2$ 两个独立的同分布正态分布变量
$G_1～N(0,1),G_2～N(0,1)$ 那么根据我们的推导
$K=L^0\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^sds_1dWs+L^1b(x_{t_0})\int_{t_0}^t\int_{t_0}^s dWs_1ds+L^0b(x_{t_0})\int_{t_0}^t\int_{t_0}^s ds_1ds+(L^1)^2\sigma(x_{t_0})\int_{t_0}^t\int_{t_0}^s\int _{t_0}^{s_1}d{Ws_2}dWs_1dWs$
$K=L^0\sigma(x_{t_0})[\triangle t \triangle W-\triangle Z]+L^1b(x_{t_0})\triangle Z+L^0b(x_{t_0})\int_{t_0}^t\int_{t_0}^s ds_1ds+(L^1)^2\sigma(x_{t_0})[\frac{1}{6}\triangle W^2-\frac{1}{2}\triangle t]\triangle W$ 其中:
$\triangle W= \sqrt{\triangle t} G_1.$
$\triangle Z= \frac{1}{2} \triangle t\triangle W+\frac{1}{2 \sqrt{3}} \triangle t^{\frac{3}{2}} G_2$ 此时 $K$ 中所有成分可求，这样的1.5阶格式为:

$x_t=x_{t_0}+b(x_{t_0})(t-t_0)+\sigma(x_{t_0})(W_t-W_{t_0})+\frac{1}{2}\sigma^{'}(x_{t_0})\sigma(x_{t_0})[(W_t-W_{t_0})^2-(t-t_0)]+K$
$K=L^0\sigma(x_{t_0})[\triangle t \triangle W-\triangle Z]+L^1b(x_{t_0})\triangle Z+L^0b(x_{t_0})\int_{t_0}^t\int_{t_0}^s ds_1ds+(L^1)^2\sigma(x_{t_0})[\frac{1}{6}\triangle W^2-\frac{1}{2}\triangle t]\triangle W$ 当然，读者可以自行推导更高阶的 $I t o - T a y l or$ 展开，笔者在这里不过多介绍。