文章

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

讨论如何在 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 对象,仅需要对 driftevolve 方法作必要的修改。

“干净”的方法

当然,“干净”的方法要改变当前接口:

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

实现

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#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

示例

下面模拟两条曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#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;
}

参考文献

  • 《金融工程中的蒙特卡罗方法》
本文由作者按照 CC BY 4.0 进行授权