QuantLib 金融计算——模拟跳扩散过程

由于版本问题，代码可能与最新版不兼容。

如果未做特别说明，文中的程序都是 C++11 代码。

QuantLib 金融计算——模拟跳扩散过程

跳扩散过程

1976 年，Merton 最早在衍生品定价中引入并分析了跳扩散过程，正因为如此 QuantLib 中和跳扩散相关的随机过程类称之为 Merton76Process，一个一般的跳扩散过程可以由下面的 SDE 描述，

\[\frac{dS(t)}{S(t-)} = \mu dt + \sigma dW(t) + dJ(t)\\ J(t) = \sum_{j=1}^{N(t)}(Y_j-1)\]

其中 $Y_j$ 是随机变量，$N(t)$ 是计数过程。$dJ(t)$ 表示 $J(t)$ 在 $t$ 时刻发生的跳跃幅度，如果发生一次跳，则幅度为 $Y_j-1$；如果没有发生跳，则幅度为 0。

在应用于衍生品定价时，需要对上述 SDE 中的项做出一些特殊约定，常见的约定有：

1. $N(t)$ 是参数等于 $\lambda$ 的 Poisson 过程；
2. $\log Y_j$ 服从正态分布 $N(\mu_{jump}, \sigma_{jump}^{2})$；
3. $\mu = r - \lambda m$，其中 $m = E[Y_j] - 1$，$r$ 代表无风险利率

模拟算法

令 $X(t) = \log S(t)$，那么

\[\begin{aligned} X(t_{i+1}) = & X(t_i) + \left(\mu - \frac12 \sigma^2 \right)(t_{i+1} - t_i)\\ & +\sigma[W(t_{i+1}) - W(t_i)] + \sum_{j = N(t_i)+1}^{N(t_{i+1})}\log Y_j \end{aligned}\]

在 $X(t_i)$ 的基础上模拟 $X(t_{i+1})$ 的算法如下：

1. 生成 $Z \sim N(0,1)$；
2. 生成 $N \sim \text{Poisson}(\lambda(t_{i+1}-t_i))$，若 $N=0$，则令 $M=0$，并转到第 4 步；
3. 生成 $\log Y_1,\dots,\log Y_N$，令 $M = \log Y_1+\dots+\log Y_N$；
4. 令 $X(t_{i+1}) = X(t_i) + \left(\mu - \frac12 \sigma^2 \right)(t_{i+1} - t_i) +\sigma \sqrt{t_{i+1} - t_i} Z + M$

那么

\[\begin{aligned} S_{t_{i+1}} =& S_{t_i} e^{(r-\lambda m -\frac{1}{2}\sigma^{2})\Delta t+ \sigma Z \sqrt{\Delta t} + M}\\ =& S_{t_i} e^{(r-\frac{1}{2}\sigma^{2})\Delta t+ \sigma Z \sqrt{\Delta t}}e^{(-\lambda m)\Delta t+M} \end{aligned}\]

其中，$\Delta t = t_{i+1} - t_i$，而 $e^{(-\lambda m)\Delta t+M}$ 是跳扩散相对于一般 Black-Scholes-Merton 过程的修正项。

面临的问题

目前 QuantLib 中提供的 Merton76Process 类只能配合“解析定价引擎”使用，本身不具备模拟随机过程路径的功能。究其原因，问题出在 QuantLib 的编码约定和 StochasticProcess1D 提供的接口两方面：

1. QuantLib 中约定 StochasticProcess 派生出的子类仅能描述 SDE 的结构信息，也就是 SDE 的参数、漂移和扩散项的函数形式，子类不携带有关随机数生成的信息，所有随机数生成的相关信息均由 Monte Carlo 框架中其他组件控制；
2. 生成随机过程路径的核心函数是 evolve 方法，StochasticProcess1D 提供的接口是 evolve(Time t0, Real x0, Time dt, Real dw)。如果 Merton76Process 按约定实现 evolve 方法的话，形式必须是 evolve(Time t0, Real x0, Time dt, const Array &dw)，因为模拟跳需要额外的随机性，所以 dw 必须是一个 Array。很明显，不匹配。

“脏”的方法

在不改变当前接口的前提下，若要实现模拟跳扩散过程，需要用比较“脏”一点儿的手段，即打破约定，让随机过程类携带一个随机数发生器，为模拟跳提供额外的随机性。

具体来说，需要声明一个 Merton76Process 的派生类，该类携带一个高斯随机数发生器。因为从数学上来讲跳扩散过程推广自一般 Black-Scholes-Merton 过程，添加了一个修正项，所以遵循“适配器模式”（或“装饰器模式”）的思想，绝大部分计算可以委托给一个 BlackScholesMertonProcess 对象，仅需要对 drift 和 evolve 方法作必要的修改。

“干净”的方法

当然，“干净”的方法要改变当前接口：

1. 声明一个和 StochasticProcess1D 平行的新类 StochasticProcess1DJump，二者唯一的区别是 evolve 方法，在 StochasticProcess1DJump 中形式是 evolve(Time t0, Real x0, Time dt, const Array &dw)；
2. 将 Merton76Process 改成继承自 StochasticProcess1DJump。

实现

下面的代码实现了前面提到的“脏”的方法，因为随机数发生器的种类有很多，且没有基类提供统一的接口，所以使用了模板技术让类可以接受不同类型的随机数发生器。同时，许多计算被委托给了一个 BlackScholesMertonProcess 对象。

#ifndef MERTON76JUMPDIFFUSIONPROCESS_HPP
#define MERTON76JUMPDIFFUSIONPROCESS_HPP

#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/math/distributions/poissondistribution.hpp>
#include <ql/math/randomnumbers/boxmullergaussianrng.hpp>
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
#include <ql/processes/blackscholesprocess.hpp>
#include <ql/processes/merton76process.hpp>

namespace QuantLib
{
template<typename GAUSS_RNG>
class Merton76JumpDiffusionProcess : public Merton76Process
{
  public:
    Merton76JumpDiffusionProcess(const Handle<Quote>& stateVariable,
                                 const Handle<YieldTermStructure>& dividendTS,
                                 const Handle<YieldTermStructure>& riskFreeTS,
                                 const Handle<BlackVolTermStructure>& blackVolTS,
                                 const Handle<Quote>& jumpInt,
                                 const Handle<Quote>& logJMean,
                                 const Handle<Quote>& logJVol,
                                 const GAUSS_RNG& gauss_rng,
                                 const ext::shared_ptr<discretization>& disc =
                                     ext::shared_ptr<discretization>(new EulerDiscretization))
      : Merton76Process(
            stateVariable,
            dividendTS,
            riskFreeTS,
            blackVolTS,
            jumpInt,
            logJMean,
            logJVol,
            disc)
      , blackProcess_(
            new BlackScholesMertonProcess(
                stateVariable,
                dividendTS,
                riskFreeTS,
                blackVolTS,
                disc))
      , gauss_rng_(gauss_rng)
    {
    }
    virtual ~Merton76JumpDiffusionProcess() {}

    Real x0() const
    {
        return blackProcess_->x0();
    }
    Time time(const Date& d) const
    {
        return blackProcess_->time(d);
    }
    Real diffusion(Time t,
                   Real x) const
    {
        return blackProcess_->diffusion(t, x);
    }
    Real apply(Real x0,
               Real dx) const
    {
        return blackProcess_->apply(x0, dx);
    }
    Size factors() const
    {
        return 1;
    }
    Real drift(Time t,
               Real x) const
    {
        Real lambda_ = Merton76Process::jumpIntensity()->value();
        Real delta_ = Merton76Process::logJumpVolatility()->value();
        Real nu_ = Merton76Process::logMeanJump()->value();
        Real m_ = std::exp(nu_ + 0.5 * delta_ * delta_) - 1;

        return blackProcess_->drift(t, x) - lambda_ * m_;
    }
    Real evolve(Time t0,
                Real x0,
                Time dt,
                Real dw) const;

  private:
    const CumulativeNormalDistribution cumNormalDist_;
    ext::shared_ptr<GeneralizedBlackScholesProcess> blackProcess_;
    GAUSS_RNG gauss_rng_;
};

template<typename GAUSS_RNG>
Real Merton76JumpDiffusionProcess<GAUSS_RNG>::evolve(Time t0,
                                                     Real x0,
                                                     Time dt,
                                                     Real dw) const
{
    Real lambda_ = Merton76Process::jumpIntensity()->value();
    Real delta_ = Merton76Process::logJumpVolatility()->value();
    Real nu_ = Merton76Process::logMeanJump()->value();
    Real m_ = std::exp(nu_ + 0.5 * delta_ * delta_) - 1;

    Real p = cumNormalDist_(gauss_rng_.next().value);
    if (p < 0.0)
        p = 0.0;
    else if (p >= 1.0)
        p = 1.0 - QL_EPSILON;

    Real j = gauss_rng_.next().value;
    const Real n = InverseCumulativePoisson(lambda_ * dt)(p);
    Real retVal = blackProcess_->evolve(
        t0, x0, dt, dw);
    retVal *=
        std::exp(-lambda_ * m_ * dt + nu_ * n + delta_ * std::sqrt(n) * j);

    return retVal;
}
}
#endif // MERTON76JUMPDIFFUSIONPROCESS_HPP

示例

下面模拟两条曲线

#include <iostream>

#include <ql/math/randomnumbers/boxmullergaussianrng.hpp>
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
#include <ql/processes/blackscholesprocess.hpp>
#include <ql/quotes/simplequote.hpp>
#include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
#include <ql/termstructures/yield/flatforward.hpp>
#include <ql/time/calendars/target.hpp>
#include <ql/time/date.hpp>
#include <ql/time/daycounters/actualactual.hpp>

#include "Merton76JumpDiffusionProcess.hpp"

int main() {
    using namespace std;
    using namespace QuantLib;

    Date refDate = Date(27, Mar, 2019);
    Rate riskFreeRate = 0.03;
    Rate dividendRate = 0.01;
    Real spot = 52.0;
    Rate vol = 0.2;
    Calendar cal = TARGET();
    DayCounter dc = ActualActual();

    ext::shared_ptr<YieldTermStructure> rdStruct(
        new FlatForward(refDate, riskFreeRate, dc));
    ext::shared_ptr<YieldTermStructure> rqStruct(
        new FlatForward(refDate, dividendRate, dc));
    Handle<YieldTermStructure> rdHandle(rdStruct);
    Handle<YieldTermStructure> rqHandle(rqStruct);

    ext::shared_ptr<SimpleQuote> spotQuote(new SimpleQuote(spot));
    Handle<Quote> spotHandle(spotQuote);

    ext::shared_ptr<BlackVolTermStructure> volQuote(
        new BlackConstantVol(refDate, cal, vol, dc));
    Handle<BlackVolTermStructure> volHandle(volQuote);

    // Specify the jump intensity, jump mean and
    // jump volatility objects

    Real jumpIntensity = 0.2;    // lambda
    Real jumpVolatility = 0.3;
    Real jumpMean = 0.0;

    ext::shared_ptr<SimpleQuote> jumpInt(new SimpleQuote(jumpIntensity));
    ext::shared_ptr<SimpleQuote> jumpVol(new SimpleQuote(jumpVolatility));
    ext::shared_ptr<SimpleQuote> jumpMn(new SimpleQuote(jumpMean));

    Handle<Quote> jumpI(jumpInt), jumpV(jumpVol), jumpM(jumpMn);

    ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
        new BlackScholesMertonProcess(
            spotHandle, rqHandle, rdHandle, volHandle));

    unsigned long seed = 12324u;
    MersenneTwisterUniformRng unifMt(seed);
    MersenneTwisterUniformRng unifMtJ(25u);

    typedef BoxMullerGaussianRng<MersenneTwisterUniformRng> GAUSS;

    GAUSS bmGauss(unifMt);
    GAUSS jGauss(unifMtJ);

    ext::shared_ptr<Merton76JumpDiffusionProcess<GAUSS>> mtProcess(
        new Merton76JumpDiffusionProcess<GAUSS>(
            spotHandle, rqHandle, rdHandle, volHandle,
            jumpI, jumpM, jumpV, jGauss));

    Time dt = 0.004, t = 0.0;
    Real x = spotQuote->value();
    Real y = spotQuote->value();
    Real dw;
    Size numVals = 250;

    std::cout << "Time, Jump, NoJump" << std::endl;
    std::cout << t << ", " << x << ", " << y << std::endl;

    for (Size j = 1; j <= numVals; ++j) {
        dw = bmGauss.next().value;
        x = mtProcess->evolve(t, x, dt, dw);
        y = bsmProcess->evolve(t, y, dt, dw);
        std::cout << t + dt << ", " << x << ", " << y << std::endl;
        t += dt;
    }

    return EXIT_SUCCESS;
}

参考文献

• 《金融工程中的蒙特卡罗方法》