1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2004 Ferdinando Ametrano
5 Copyright (C) 2007 StatPro Italia srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include <ql/exercise.hpp>
22#include <ql/math/distributions/poissondistribution.hpp>
23#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
24#include <ql/pricingengines/vanilla/jumpdiffusionengine.hpp>
25#include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
26#include <ql/termstructures/yield/flatforward.hpp>
27#include <ql/utilities/dataformatters.hpp>
28#include <utility>
29
30namespace QuantLib {
31
32 JumpDiffusionEngine::JumpDiffusionEngine(ext::shared_ptr<Merton76Process> process,
33 Real relativeAccuracy,
34 Size maxIterations)
35 : process_(std::move(process)), relativeAccuracy_(relativeAccuracy),
36 maxIterations_(maxIterations) {
37 registerWith(h: process_);
38 }
39
40
41 void JumpDiffusionEngine::calculate() const {
42
43 Real jumpSquareVol = process_->logJumpVolatility()->value()
44 * process_->logJumpVolatility()->value();
45 Real muPlusHalfSquareVol = process_->logMeanJump()->value()
46 + 0.5*jumpSquareVol;
47 // mean jump size
48 Real k = std::exp(x: muPlusHalfSquareVol) - 1.0;
49 Real lambda = (k+1.0) * process_->jumpIntensity()->value();
50
51 ext::shared_ptr<StrikedTypePayoff> payoff =
52 ext::dynamic_pointer_cast<StrikedTypePayoff>(r: arguments_.payoff);
53 QL_REQUIRE(payoff, "non-striked payoff given");
54
55 Real variance =
56 process_->blackVolatility()->blackVariance(
57 d: arguments_.exercise->lastDate(),
58 strike: payoff->strike());
59
60 DayCounter voldc = process_->blackVolatility()->dayCounter();
61 Calendar volcal = process_->blackVolatility()->calendar();
62 Date volRefDate = process_->blackVolatility()->referenceDate();
63 Time t = voldc.yearFraction(d1: volRefDate,
64 d2: arguments_.exercise->lastDate());
65 Rate riskFreeRate = -std::log(x: process_->riskFreeRate()->discount(
66 d: arguments_.exercise->lastDate()))/t;
67 Date rateRefDate = process_->riskFreeRate()->referenceDate();
68
69 PoissonDistribution p(lambda*t);
70
71 Handle<Quote> stateVariable = process_->stateVariable();
72 Handle<YieldTermStructure> dividendTS = process_->dividendYield();
73 RelinkableHandle<YieldTermStructure> riskFreeTS(
74 *process_->riskFreeRate());
75 RelinkableHandle<BlackVolTermStructure> volTS(
76 *process_->blackVolatility());
77
78 ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
79 new GeneralizedBlackScholesProcess(stateVariable, dividendTS,
80 riskFreeTS, volTS));
81
82 AnalyticEuropeanEngine baseEngine(bsProcess);
83
84 auto* baseArguments = dynamic_cast<VanillaOption::arguments*>(baseEngine.getArguments());
85
86 baseArguments->payoff = arguments_.payoff;
87 baseArguments->exercise = arguments_.exercise;
88
89 baseArguments->validate();
90
91 const auto* baseResults =
92 dynamic_cast<const VanillaOption::results*>(baseEngine.getResults());
93
94 results_.value = 0.0;
95 results_.delta = 0.0;
96 results_.gamma = 0.0;
97 results_.theta = 0.0;
98 results_.vega = 0.0;
99 results_.rho = 0.0;
100 results_.dividendRho = 0.0;
101
102 Real r, v, weight, lastContribution = 1.0;
103 Size i;
104 Real theta_correction;
105 // Haug arbitrary criterium is:
106 //for (i=0; i<11; i++) {
107 for (i=0; (lastContribution>relativeAccuracy_ && i<maxIterations_)
108 || i < Size(lambda*t); i++) {
109
110 // constant vol/rate assumption. It should be relaxed
111 v = std::sqrt(x: (variance + i*jumpSquareVol)/t);
112 r = riskFreeRate - process_->jumpIntensity()->value()*k
113 + i*muPlusHalfSquareVol/t;
114 riskFreeTS.linkTo(h: ext::shared_ptr<YieldTermStructure>(new
115 FlatForward(rateRefDate, r, voldc)));
116 volTS.linkTo(h: ext::shared_ptr<BlackVolTermStructure>(new
117 BlackConstantVol(rateRefDate, volcal, v, voldc)));
118
119 baseArguments->validate();
120 baseEngine.calculate();
121
122 weight = p(Size(i));
123 results_.value += weight * baseResults->value;
124 results_.delta += weight * baseResults->delta;
125 results_.gamma += weight * baseResults->gamma;
126 results_.vega += weight * (std::sqrt(x: variance/t)/v)*
127 baseResults->vega;
128 // theta modified
129 theta_correction = baseResults->vega*((i*jumpSquareVol)/
130 (2.0*v*t*t)) +
131 baseResults->rho*i*muPlusHalfSquareVol/(t*t);
132 results_.theta += weight *(baseResults->theta + theta_correction +
133 lambda*baseResults->value);
134 if(i != 0){
135 results_.theta -= (p(Size(i-1))*lambda* baseResults->value);
136 }
137 //end theta calculation
138 results_.rho += weight * baseResults->rho;
139 results_.dividendRho += weight * baseResults->dividendRho;
140
141 lastContribution = std::fabs(x: baseResults->value /
142 (std::fabs(x: results_.value)>QL_EPSILON ? results_.value : 1.0));
143
144 lastContribution = std::max<Real>(a: lastContribution,
145 b: std::fabs(x: baseResults->delta /
146 (std::fabs(x: results_.delta)>QL_EPSILON ? results_.delta : 1.0)));
147
148 lastContribution = std::max<Real>(a: lastContribution,
149 b: std::fabs(x: baseResults->gamma /
150 (std::fabs(x: results_.gamma)>QL_EPSILON ? results_.gamma : 1.0)));
151
152 lastContribution = std::max<Real>(a: lastContribution,
153 b: std::fabs(x: baseResults->theta /
154 (std::fabs(x: results_.theta)>QL_EPSILON ? results_.theta : 1.0)));
155
156 lastContribution = std::max<Real>(a: lastContribution,
157 b: std::fabs(x: baseResults->vega /
158 (std::fabs(x: results_.vega)>QL_EPSILON ? results_.vega : 1.0)));
159
160 lastContribution = std::max<Real>(a: lastContribution,
161 b: std::fabs(x: baseResults->rho /
162 (std::fabs(x: results_.rho)>QL_EPSILON ? results_.rho : 1.0)));
163
164 lastContribution = std::max<Real>(a: lastContribution,
165 b: std::fabs(x: baseResults->dividendRho /
166 (std::fabs(x: results_.dividendRho)>QL_EPSILON ?
167 results_.dividendRho : 1.0)));
168
169 lastContribution *= weight;
170 }
171 QL_ENSURE(i<maxIterations_,
172 i << " iterations have been not enough to reach "
173 << "the required " << relativeAccuracy_
174 << " accuracy. The " << io::ordinal(i)
175 << " addendum was " << lastContribution
176 << " while the running sum was " << results_.value);
177 }
178
179}
180
181

source code of quantlib/ql/pricingengines/vanilla/jumpdiffusionengine.cpp