目录
- 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 中的项做出一些特殊约定，常见的约定有：
- \(N(t)\) 是参数等于 \(\lambda\) 的 Poisson 过程；
- \(\log Y_j\) 服从正态分布 \(N(\mu_{jump}, \sigma_{jump}^{2})\)；
- \(\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})\) 的算法如下：
- 生成 \(Z \sim N(0,1)\)；
- 生成 \(N \sim \text{Poisson}(\lambda(t_{i+1}-t_i))\)，若 \(N=0\)，则令 \(M=0\)，并转到第 4 步；
- 生成 \(\log Y_1,\dots,\log Y_N\)，令 \(M = \log Y_1+\dots+\log Y_N\)；
- 令 \(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 提供的接口两方面：
- QuantLib 中约定
StochasticProcess派生出的子类仅能描述 SDE 的结构信息，也就是 SDE 的参数、漂移和扩散项的函数形式，子类不携带有关随机数生成的信息，所有随机数生成的相关信息均由 Monte Carlo 框架中其他组件控制； - 生成随机过程路径的核心函数是
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 方法作必要的修改。
“干净”的方法
当然，“干净”的方法要改变当前接口：
- 声明一个和
StochasticProcess1D平行的新类StochasticProcess1DJump，二者唯一的区别是evolve方法，在StochasticProcess1DJump中形式是evolve(Time t0, Real x0, Time dt, const Array &dw)； - 将
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 objectsReal jumpIntensity = 0.2; // lambdaReal 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;
}
参考文献
- 《金融工程中的蒙特卡罗方法》
