QuantLib 金融計算——高級話題之模擬跳擴散過程

時間  2019-12-08
標籤 quantlib 金融 計算 高級 話題 模擬 擴散 過程 简体版
原文   原文鏈接

目錄ios

若是未作特別說明，文中的程序都是 C++11 代碼。算法

QuantLib 金融計算——高級話題之模擬跳擴散過程

跳擴散過程

1976 年，Merton 最先在衍生品訂價中引入並分析了跳擴散過程，正由於如此 QuantLib 中和跳擴散相關的隨機過程類稱之爲 Merton76Process，一個通常的跳擴散過程能夠由下面的 SDE 描述，app

\[ \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。dom

在應用於衍生品訂價時，須要對上述 SDE 中的項作出一些特殊約定，常見的約定有：ide

  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} \]ui

在 \(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\)

那麼spa

\[ \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;
}

參考文獻

  1. 《金融工程中的蒙特卡羅方法》
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程