专栏名称: 风控斋

一种随机微分方程数值模拟的另类方法

风控斋 · 知乎专栏 · · 2017-08-20 00:08

正文

最近看到这么一个话题， Stochastic Calculus 在现实的 Quant 世界中到底可以干什么，惊动了不少大咖的关注。这让我想起了许多年前在设计distressed asset风险模型的时候曾经用到的一种数值模拟 SDE的方法。其实这种方法并不是什么新方法，因为其发表距今已经10多年了。这篇文章的作者之一现在是JPMorgan Rates QR的MD。后文为了叙述方便，我称这种方法为WA(Weak Approximation, 见原文标题）.

原文链接

与传统的SDE数值模拟方法，如大家熟知的Euler和Milstein方法相比，WA方法的思路非常不同，其物理的意义更浓郁一些，首先其出发点是SDE的Stratonovich形式，而非在金融数学中大家普遍使用的Ito形式。相对于后者，前者的优点是复合函数的微分满足传统微积分的链式法则(chain rule), 而不是Ito's Lemma. 这对于许多不太熟悉金融数学的人来说更加直观。此外，这种方法的关键思想是利用物理中类似路径积分算子分解的思想，将SDE在每一步的求解过程分为三步，每一步求解一个常微分方程。由于这种简化，每一步的常微分方程往往较易求解，甚至可以通过解析解求得而无须诉诸于离散化的数值方法。

这里为了直观起见，仅以一维的SDE为例，并且在叙述中采用大家更加熟悉的方式而非是抽象的算符形式。关于一般的情况，特别是多维SDE, 大家有兴趣可以参考原文。

考虑如下的一维标量SDE:

$d X_t = \mu(X_t, t) dt + \sigma (X_t, t) dW_t$

给定初始条件 $X(0)=X_0$ , 我们的目标是求其在t时间t的解

$X_t = X_0 + \int^{t}_{0} \mu(X_s, s) ds + \int^{t}_{0} \sigma(X_s,s) d W_{s}$

我们通常熟悉的两种方法是通过离散化上述SDE来求解该方程，即从t=0开始，每一步进行离散化更新 $X_{t}$ .

Euler 方法： $X(t+\Delta t) = X(t) + \mu(X_t,t) \Delta t + \sigma(X_t,t) Z_{t} \sqrt{t}$
Milstein 方法： $X(t+\Delta t) = X(t) + \mu(X_t,t) \Delta t + \sigma(X_t,t) Z_{t} \sqrt{t} + \frac{1}{2} \sigma (X_t, t) \frac{\partial}{\partial x} \sigma(X_t,t)\left[ Z^{2}_{t} \Delta t - (\Delta t)^{2} ] \right]$

这里 $Z_t \sim N(0,1)$ .

为了保证解的收敛性和稳定性，时间间隔 $\Delta t$ 往往需要取的较小。这样如果需要求得期限较长的解需要的部数就会很多。而 WA方法则可以有效地降低求解步数从而提高求解效率。

用WA求解的步骤如下：

将上述的SDE转化为Stratonovich形式： $d X_t = \mu^{'}(X_t, t) dt + \sigma (X_t, t) \circ dW_t \\ \mu^{'}(X_t, t) = \mu(X_t,t) - \frac{1}{2} \sigma(X_t, t) \frac{\partial}{\partial X_{t}} \sigma(X_t,t)$ （注意多出的一项与Milstein 的最后一项在形式上的相似性。）
将方程离散化，每一步分别求解3个常微分方程，其解可以分别用如下的3个积分来表示：

$X^{1}_{t, \Delta t} = X_{t} + \int^{t +\frac{1}{2}\Delta t}_{t} \mu^{'}(X_s,s)ds\\ X^{2}_{t, \Delta t} = X^{1}_{t, \Delta t} + \int^{t + \frac{1}{2}\Delta t + \sqrt{\Delta t} Z_{t} }_{t+\frac{1}{2} \Delta t} \sigma(X_{s}, s) ds \\ X_{t+\Delta t} = X^{2}_{t, \Delta t} + \int^{t + \Delta t}_{t+\frac{1}{2} \Delta t} \mu^{'}(X_s,s)ds$

这样就完成了一步的求解过程。

我们以一个具体的例子来说明这一求解过程。考虑大家都很熟悉的 Cox-Ingersoll-Ross (CIR)模型：

$d X_t = \kappa (\theta - X_t) dt + \sigma \sqrt{X_t} dW_t$

首先将其化为 Stratonovich形式 $d X_t = \kappa (\theta^{'} - X_t ) dt + \sigma \sqrt{X_t} \circ dW_t\\ \theta^{'} = \theta - \frac{ \sigma^{2}}{4\kappa}$
求解三个常微分方程，得到每步的更新规则（在第一和第三步保证解为证）。

$X^{1}_{t, \Delta t} = max(0, \theta^{'} + (X_{t} - \theta^{'}) e^{-\frac{\kappa \Delta t}{2}})\\ X^{2}_{t, \Delta t } = \frac{1}{4} (2 \sqrt{ X^{1}_{t, \Delta t} } + \sigma Z_{t} \sqrt{\Delta t})^2\\ X_{t+\Delta t} =max(0, \theta^{'} + (X^{2}_{t, \Delta t } - \theta^{'}) e^{-\frac{\kappa \Delta t}{2}})$

我们都知道CIR模型的理论解服从非中心卡方分布（non-central $\chi^{2}$ distribution)，其自由度和偏中心率为 $X_{0}, \kappa, \theta$ 的函数。这使得我们可以将理论分布与模拟分布进行比较。下图表明，用Milstein与WA模拟的SDE的解都服从非中心卡方分布。

模拟的参数如下： $\kappa = 0.1\\ \theta = 0.5 \\ \sigma = 0.2\\ X_{0} = 0.45\\ \Delta t = 0.01\\ t = 10$

模拟次数为10000次。

一种随机微分方程数值模拟的另类方法

正文

请到「今天看啥」查看全文