随机微分方程的模拟（一）

无论是资产定价还是风险管理，都要涉及到SDE的蒙特卡洛模拟。

本文仅讨论一维SDE的数值模拟。

对大多数微分方程，我们无法求出解析解，所以为了确定这个方程所对应的函数，就需要进行数值模拟。随机微分方程也是一样的，对于SDE：
$d X_{t}=\mu\left(X_{t}, t\right) d t+\sigma\left(X_{t}, t\right) d W_{t}$ 大多数情况下，要求出随机过程 $X_t$ 在某个时刻的分布，还得靠随机模拟，模拟出若干条路径，来得到近似的分布。模拟的方法，主要就是两种，Euler–Maruyama method 和 Milstein method。

Euler 方法是最直接的方法，扩散过程（diffusion process）是一个马尔可夫过程，因为 $W_{t+\Delta t}-W_t\sim \mathcal{N}(0, \Delta t)$ ，自然就有：
$\small X_{t+\Delta t} \simeq X_t+\mu\left(X_{t}, t\right) \Delta t+\sigma\left(X_{t}, t\right) Z \sqrt{\Delta t}$

Milstein 方法会更加精确一些： $X_{t+\Delta t}\simeq X_t+\mu\left(X_{t}, t\right) \Delta t+\sigma\left(X_{t}, t\right) Z_{t} \sqrt{t}+\frac{1}{2} \sigma\left(X_{t}, t\right) \frac{\partial}{\partial x} \sigma\left(X_{t}, t\right)\left[Z^{2} \Delta t-(\Delta t)^{2}\right]$

其中 $Z\sim \mathcal{N}(0, 1)$ 。

Milstein 方法的初衷是让估计：
$\int_t^{t+\Delta t}\sigma(X_u, u)dW_u=\sigma(X_t, t)\int_t^{t+\Delta t}dW_u$ 变得稍微精确一点。所以，对 $\sigma(X_t, t)$ 应用伊藤公式，求出 $d [\sigma(X_t, t)]$ 关于 $dW_t$ 和 $dt$ 的关系，再积分得到 $\sigma(X_t, t)$ 的表达式，最后估计 $\int_t^{t+\Delta t}\sigma(X_u, u)dW_u$ 。这是得到米尔斯坦方法的思路，沿着这个思路下去就是随机泰勒展开！

注意到 Milstein 方法需要 $\sigma(X_t, t)$ 是二阶可微的。

米尔斯坦方法对应到几何布朗运动 ${d} X_{t}=\mu X {d} t+\sigma X d W_{t}$ 就是：
$X_{t+\Delta t}=X_t + \mu X_t \Delta t + \sigma X_t \Delta W_t+\frac{1}{2}\sigma^2\left[(\Delta W_t)^2-\Delta t \right]$ 它其实是在 Euler 方法的基础上增加了更高阶的小量。

SDE的模拟本质上我认为还是在做数值积分。

order of convergence

Milstein 方法比 Euler 方法要精确，那么这种精确程度该如何来度量呢？
这就要引入收敛阶数的概念，而收敛性分为两种，强收敛和弱收敛。

对于随机过程 $X_t$ ，我们以 $\Delta t$ 构造等距节点，使用某种方法插值模拟出的随机过程记为 $Y^{\Delta}_t$ 。

强收敛

对于任意的时间 $T$ ，如果存在不依赖于 $\Delta t$ 的常数 $C$ 和 $\delta$ ，使得：
${E}\left(\left|{X}_{{T}}-{Y}_{{T}}^{{\Delta}}\right|\right) \leq {C} {\Delta t}^{\gamma}, \quad \forall \Delta t \in (0, \delta)$ 我们就说这个方法是 $\gamma$ 阶强收敛的！

弱收敛

对于任意的时间 $T$ 和可测函数 $g$ ，如果存在不依赖于 $\Delta t$ 的常数 $C$ 和 $\delta$ ，使得：
$|E(g(X_T)) - E(g(Y_T^{\Delta}))| \leq C {\Delta t}^\beta, \quad \forall \Delta t \in (0, \delta)$ 我们就说这个方法是 $\beta$ 阶弱收敛的！

收敛阶数不仅与模拟方法有关，还与这个随机微分方程的性质有关。 但是一般来说，我们认为 Euler 方法的强收敛阶数是 0.5，Milstein 方法是1；Euler 方法的弱收敛阶数是1。

stochastic taylor expansion

其实 Euler 方法和 Milstein 方法都可以从随机泰勒展开中得到，并且可以类似的分析它们的收敛阶数。

首先考虑一个自治（autonomous）的ODE：
$\frac{d}{d t} X(t)=a[X(t)]$ 令 $f$ 是 $X_t$ 的函数，那么根据链导法则：
$\frac{d}{d t} f[X(t)]=a[X(t)] \frac{\partial}{\partial X} f[X(t)]$ 定义算子 $\mathcal{L} \equiv a(X) \frac{\partial}{\partial X}$ ，上式可以写成积分的形式：
$f[X(t)]=f\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{t} \mathcal{L} f\left[X\left(\tau_{1}\right)\right] d \tau_{1}$ 对 $f\left[X\left(\tau_{1}\right)\right]$ 也可以写成上式一般的展开形式，继而有：
$\begin{aligned} f[X(t)] &=f\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{t} \mathcal{L}\left[f\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{\tau_{1}} \mathcal{L} f\left[X\left(\tau_{2}\right)\right] d \tau_{2}\right] d \tau_{1} \\ &=f\left[X\left(t_{0}\right)\right]+\mathcal{L} f\left[X\left(t_{0}\right)\right] \int_{t_{0}}^{t} d \tau_{1}+\int_{t_{0}}^{t} \int_{t_{0}}^{\tau_{1}} \mathcal{L}^{2} f\left[X\left(\tau_{2}\right)\right] d \tau_{2} d \tau_{1} \\ &=f\left[X\left(t_{0}\right)\right]+\mathcal{L} f\left[X\left(t_{0}\right)\right]\left(t-t_{0}\right)+\int_{t_{0}}^{t} \int_{t_{0}}^{\tau_{1}} \mathcal{L}^{2} f\left[X\left(\tau_{2}\right)\right] d \tau_{2} d \tau_{1} \end{aligned}$ 因为：
$\begin{aligned} \mathcal{L}^{2} f\left[X\left(\tau_{2}\right)\right] &=\mathcal{L}^{2}\left[f\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{\tau_{2}} \mathcal{L} f\left[X\left(\tau_{3}\right)\right] d \tau_{3}\right] \\ &=\mathcal{L}^{2} f\left[X\left(t_{0}\right)\right]+\underbrace{\mathcal{L}^{2} \int_{t_{0}}^{\tau_{2}} \mathcal{L} f\left[X\left(\tau_{3}\right)\right] d \tau_{3}}_{O\left(\mathcal{L}^{3}\right)} \end{aligned}$ 上面其实就是我们熟悉的Taylor展开：
$f[X(t)]=f\left[X\left(t_{0}\right)\right]+\mathcal{L} f\left[X\left(t_{0}\right)\right]\left(t-t_{0}\right)+\frac{1}{2} \mathcal{L}^{2} f\left[X\left(t_{0}\right)\right]\left(t-t_{0}\right)^{2}+O\left(\mathcal{L}^{3}\right)$

接下来我们把上述方法应用到自治的SDE中。定义算子：
$\begin{array}{l} \mathcal{L}^{0} \equiv a \frac{\partial}{\partial X}+\frac{1}{2} b^{2}[X] \frac{\partial^{2}}{\partial X^{2}} \\ \mathcal{L}^{1} \equiv b[X] \frac{\partial}{\partial X} \end{array}$ 那么伊藤公式
$d X(t)=a[X(t)] d t+b[X(t)] d W(t)$ 可以写成：
$f[X(t)]=f\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{t} \mathcal{L}^{0} f[X(s)] d s+\int_{t_{0}}^{t} \mathcal{L}^{1} f[X(s)] d W(s)$ 接下来，

令 $f(x)=x$ ，得到：
$\Large X(t)=X\left(t_{0}\right)+\int_{t_{0}}^{t} a[X(s)] d s+\int_{t_{0}}^{t} b[X(s)] d W(s) \tag{1}$
令 $f(x)=a(x)$ ，得到：
$\Large a[X(t)]=a\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{t} \mathcal{L}^{0} a[X(s)] d s+\int_{t_{0}}^{t} \mathcal{L}^{1} a[X(s)] d W(s) \tag{2}$
令 $f(x)=b(x)$ ，得到：
$\Large b[X(t)]=b\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{t} \mathcal{L}^{0} b[X(s)] d s+\int_{t_{0}}^{t} \mathcal{L}^{1} b[X(s)] d W(s) \tag{3}$

现在我们要做的是，把伊藤公式应用到伊藤公式中，来得到一个比伊藤公式更加精确的公式！

把(2)和(3)代入(1)：

$\begin{aligned} X(t)=& X\left(t_{0}\right)+\int_{t_{0}}^{t}\left\{a\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{s_{1}} \mathcal{L}^{0} a\left[X\left(s_{2}\right)\right] d s_{2}+\int_{t_{0}}^{s_{1}} \mathcal{L}^{1} a\left[X\left(s_{2}\right)\right] d W\left(s_{2}\right)\right\} d s_{1} \\ &+\int_{t_{0}}^{t}\left\{b\left[X\left(t_{0}\right)\right]+\int_{t_{0}}^{s_{1}} \mathcal{L}^{0} b\left[X\left(s_{2}\right)\right] d s_{2}+\int_{t_{0}}^{s_{1}} \mathcal{L}^{1} b\left[X\left(s_{2}\right)\right] d W\left(s_{2}\right)\right\} d W\left(s_{1}\right) \end{aligned}$ 于是：
$\large X(t)=X\left(t_{0}\right)+a\left[X\left(t_{0}\right)\right] \int_{t_{0}}^{t} d s_{1}+b\left[X\left(t_{0}\right)\right] \int_{t_{0}}^{t} d W\left(s_{1}\right)+R \tag{4}$ 其中余项： $\large \begin{aligned} R \equiv & \int_{t_{0}}^{t} \int_{t_{0}}^{s_{1}} \mathcal{L}^{0} a\left[X\left(s_{2}\right)\right] d s_{2} d s_{1}+\int_{t_{0}}^{t} \int_{t_{0}}^{s_{1}} \mathcal{L}^{1} a\left[X\left(s_{2}\right)\right] d W\left(s_{2}\right) d s_{1} \\ &+\int_{t_{0}}^{t} \int_{t_{0}}^{s_{1}} \mathcal{L}^{0} b\left[X\left(s_{2}\right)\right] d s_{2} d W\left(s_{1}\right)+\int_{t_{0}}^{t} \int_{t_{0}}^{s_{1}} \mathcal{L}^{1} b\left[X\left(s_{2}\right)\right] d W\left(s_{2}\right) d W\left(s_{1}\right) \end{aligned} \tag{5}$ 我们上面所做的工作就是重复应用伊藤公式获取更高阶的近似！

注意到(4)式右边前三项就是 Euler 近似！所以我们可以通过分析余项 $R$ 来得到 Euler 方法的余项阶数。

分析一下(5)式中余项 $R$ 各项的阶数，如果把 $R$ 的最后一项考虑进去，我们就能得到 Milstein 方法。
注意到： $\small \int_{t_{0}}^{t} \int_{t_{0}}^{s_{1}} d W\left(s_{2}\right) d W\left(s_{1}\right)=\frac{1}{2}\left[W(t)-W\left(t_{0}\right)\right]^{2}-\frac{1}{2}\left(t-t_{0}\right)$

如果我们把余项 $R$ 的后三项都考虑进去，就能得到比 Milstein 方法更精确的数值方法！

其它的模拟方法

隐式方法

一般我们用的都是显示的 Euler 方法，但其实也可以有类似于 ODE 里的隐式方法：
$X_{t+\Delta t} \simeq X_t+\mu\left(X_{t+\Delta t}, t + \Delta t \right) \Delta t+\sigma\left(X_{t}, t\right) Z \sqrt{\Delta t}$ 我感觉一般常研究的SDE的稳定性都还不错。

Runge-Kutta Methods

Milstein 方法需要计算一次导数，我们也可以不去直接计算导数，而是去近似它。
$\left\{\begin{array}{l} \tilde{X}_{i}=X_{i}+a\left(X_{i}\right) \Delta t+b\left(X_{i}\right) \sqrt{\Delta t} \\ X_{i+1}=X_{i}+a\left(X_{i}\right) \Delta t+b\left(X_{i}\right) \Delta W_{i}+\frac{1}{2 \sqrt{\Delta t}}\left[b\left(\tilde{X}_{i}\right)-b\left(X_{i}\right)\right]\left[\left(\Delta W_{i}\right)^{2}-\Delta t\right] \end{array}\right.$

Simplified Methods

布朗运动的本质是随机游走的极限，当 $\Delta t$ 足够小的时候，用离散的随机游走去代替布朗运动，能起到节省计算时间的效果。如果目的只是估计 $E(X_T)$ ，那么 Simplified Methods 的效果也不错。

参考资料：
https://www.mimuw.edu.pl/~apalczew/CFP_lecture5.pdf