35 beta_(
beta), jumpIntensity_(jumpIntensity), eta_(eta)
37 QL_REQUIRE(eps > 0.0 && eps < 1.0, "eps > 0.0 and eps < 1.0
");
38 QL_REQUIRE(steps > 1, "minimum number of steps is two
");
40 const Real start = 0.0;
41 const Real end = 1.0-eps;
42 const Real dx = (end-start)/(steps-1);
43 const Real scale = 1/(1-std::exp(-beta/jumpIntensity));
45 for (Size i=0; i < steps; ++i) {
46 const Real p = start + i*dx;
47 locations_[i] = scale*(-1.0/eta*std::log(1.0-p));
50 for (Size i=0; i < steps-1; ++i) {
51 dminus_[i+1] = dplus_[i] = locations_[i+1]-locations_[i];
53 dplus_.back() = dminus_.front() = Null<Real>();
57 Real ExponentialJump1dMesher::jumpSizeDensity(Real x, Time t) const {
58 const Real a = 1.0-jumpIntensity_/beta_;
59 const Real norm = 1.0-std::exp(-jumpIntensity_*t);
61 = std::exp(GammaFunction().logValue(1.0-jumpIntensity_/beta_));
62 return jumpIntensity_*gammaValue/norm
63 *( incompleteGammaFunction(a, x*std::exp(beta_*t)*eta_)
64 -incompleteGammaFunction(a, x*eta_))
65 *std::pow(eta_, jumpIntensity_/beta_)
66 /(beta_*std::pow(x, a));
69 Real ExponentialJump1dMesher::jumpSizeDensity(Real x) const {
70 const Real a = 1.0-jumpIntensity_/beta_;
72 = std::exp(GammaFunction().logValue(jumpIntensity_/beta_));
73 return std::exp(-x*eta_)*std::pow(x, -a) * std::pow(eta_, 1.0-a)
77 Real ExponentialJump1dMesher::jumpSizeDistribution(Real x, Time t) const {
78 const Real xmin = std::min(x, 1.0e-100);
80 return GaussLobattoIntegral(1000000, 1e-12)(
81 [&](Real _x){ return jumpSizeDensity(_x, t); },
82 xmin, std::max(x, xmin));
85 Real ExponentialJump1dMesher::jumpSizeDistribution(Real x) const {
86 const Real a = jumpIntensity_/beta_;
87 const Real xmin = std::min(x, QL_EPSILON);
89 = std::exp(GammaFunction().logValue(jumpIntensity_/beta_));
92 (std::pow(xmin, a)/a - std::pow(xmin, a+1)/(a+1))/gammaValue;
94 return lowerEps + GaussLobattoIntegral(10000, 1e-12)(
95 [&](Real _x){ return jumpSizeDensity(_x); },
96 xmin/eta_, std::max(x, xmin/eta_));
ExponentialJump1dMesher(Size steps, Real beta, Real jumpIntensity, Real eta, Real eps=1e-3)
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
mesher for a exponential jump mesher with high mean reversion rate and low jump intensity
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
std::size_t Size
size of a container
Incomplete Gamma function.