Logo
Fully annotated reference manual - version 1.8.12
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
blackscholes.cpp
Go to the documentation of this file.
1/*
2 Copyright (C) 2019 Quaternion Risk Management Ltd
3 All rights reserved.
4
5 This file is part of ORE, a free-software/open-source library
6 for transparent pricing and risk analysis - http://opensourcerisk.org
7
8 ORE is free software: you can redistribute it and/or modify it
9 under the terms of the Modified BSD License. You should have received a
10 copy of the license along with this program.
11 The license is also available online at <http://opensourcerisk.org>
12
13 This program is distributed on the basis that it will form a useful
14 contribution to risk analytics and model standardisation, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the license for more details.
17*/
18
20
23
24#include <ql/math/comparison.hpp>
25#include <ql/math/matrixutilities/choleskydecomposition.hpp>
26#include <ql/math/matrixutilities/pseudosqrt.hpp>
27#include <ql/math/matrixutilities/symmetricschurdecomposition.hpp>
28
29namespace ore {
30namespace data {
31
32using namespace QuantLib;
33
34BlackScholes::BlackScholes(const Size paths, const std::string& currency, const Handle<YieldTermStructure>& curve,
35 const std::string& index, const std::string& indexCurrency,
36 const Handle<BlackScholesModelWrapper>& model, const McParams& mcParams,
37 const std::set<Date>& simulationDates, const IborFallbackConfig& iborFallbackConfig,
38 const std::string& calibration, const std::vector<Real>& calibrationStrikes)
39 : BlackScholes(paths, {currency}, {curve}, {}, {}, {}, {index}, {indexCurrency}, model, {}, mcParams,
40 simulationDates, iborFallbackConfig, calibration, {{index, calibrationStrikes}}) {}
41
43 const Size paths, const std::vector<std::string>& currencies, const std::vector<Handle<YieldTermStructure>>& curves,
44 const std::vector<Handle<Quote>>& fxSpots,
45 const std::vector<std::pair<std::string, QuantLib::ext::shared_ptr<InterestRateIndex>>>& irIndices,
46 const std::vector<std::pair<std::string, QuantLib::ext::shared_ptr<ZeroInflationIndex>>>& infIndices,
47 const std::vector<std::string>& indices, const std::vector<std::string>& indexCurrencies,
48 const Handle<BlackScholesModelWrapper>& model,
49 const std::map<std::pair<std::string, std::string>, Handle<QuantExt::CorrelationTermStructure>>& correlations,
50 const McParams& mcParams, const std::set<Date>& simulationDates, const IborFallbackConfig& iborFallbackConfig,
51 const std::string& calibration, const std::map<std::string, std::vector<Real>>& calibrationStrikes)
52 : BlackScholesBase(paths, currencies, curves, fxSpots, irIndices, infIndices, indices, indexCurrencies, model,
53 correlations, mcParams, simulationDates, iborFallbackConfig),
54 calibration_(calibration), calibrationStrikes_(calibrationStrikes) {}
55
57
59
60 // nothing to do if we do not have any indices
61
62 if (indices_.empty())
63 return;
64
65 // init underlying path where we map a date to a randomvariable representing the path values
66
67 for (auto const& d : effectiveSimulationDates_) {
68 underlyingPaths_[d] = std::vector<RandomVariable>(model_->processes().size(), RandomVariable(size(), 0.0));
69 if (trainingSamples() != Null<Size>())
71 std::vector<RandomVariable>(model_->processes().size(), RandomVariable(trainingSamples(), 0.0));
72 }
73
74 // compile the correlation matrix
75
76 Matrix correlation = getCorrelation();
77
78 // determine calibration strikes
79
80 std::vector<Real> calibrationStrikes;
81 if (calibration_ == "ATM") {
82 calibrationStrikes.resize(indices_.size(), Null<Real>());
83 } else if (calibration_ == "Deal") {
84 for (Size i = 0; i < indices_.size(); ++i) {
85 auto f = calibrationStrikes_.find(indices_[i].name());
86 if (f != calibrationStrikes_.end() && !f->second.empty()) {
87 calibrationStrikes.push_back(f->second[0]);
88 TLOG("calibration strike for index '" << indices_[i] << "' is " << f->second[0]);
89 } else {
90 calibrationStrikes.push_back(Null<Real>());
91 TLOG("calibration strike for index '" << indices_[i] << "' is ATMF");
92 }
93 }
94 } else {
95 QL_FAIL("BlackScholes: calibration '" << calibration_ << "' not supported, expected ATM, Deal");
96 }
97
98 // set reference date values, if there are no future simulation dates we are done
99
100 for (Size l = 0; l < indices_.size(); ++l) {
101 underlyingPaths_[*effectiveSimulationDates_.begin()][l].setAll(model_->processes()[l]->x0());
102 if (trainingSamples() != Null<Size>()) {
103 underlyingPathsTraining_[*effectiveSimulationDates_.begin()][l].setAll(model_->processes()[l]->x0());
104 }
105 }
106
107 if (effectiveSimulationDates_.size() == 1)
108 return;
109
110 // compute drift and covariances of log spots to evolve the process; the covariance computation is done
111 // on the refined grid where we assume the volatilities to be constant
112
113 // used for dirf adjustment eq / com that are not in base ccy below
114 std::vector<Size> forCcyDaIndex(indices_.size(), Null<Size>());
115 for (Size j = 0; j < indices_.size(); ++j) {
116 if (!indices_[j].isFx()) {
117 // do we have an fx index with the desired currency?
118 for (Size jj = 0; jj < indices_.size(); ++jj) {
119 if (indices_[jj].isFx()) {
121 forCcyDaIndex[j] = jj;
122 }
123 }
124 }
125 }
126
127 std::vector<Array> drift(effectiveSimulationDates_.size() - 1, Array(indices_.size(), 0.0));
128 std::vector<Matrix> sqrtCov;
130 std::vector<Matrix>(effectiveSimulationDates_.size() - 1, Matrix(indices_.size(), indices_.size(), 0.0));
131 Array variance(indices_.size(), 0.0), discountRatio(indices_.size(), 1.0);
132 Size tidx = 1; // index in the refined time grid
133
134 for (Size i = 1; i < effectiveSimulationDates_.size(); ++i) {
135
136 // covariance
137
138 while (tidx <= positionInTimeGrid_[i]) {
139 Array d_variance(indices_.size());
140 for (Size j = 0; j < indices_.size(); ++j) {
141 Real tmp = model_->processes()[j]->blackVolatility()->blackVariance(
142 timeGrid_[tidx],
143 calibrationStrikes[j] == Null<Real>()
144 ? atmForward(model_->processes()[j]->x0(), model_->processes()[j]->riskFreeRate(),
145 model_->processes()[j]->dividendYield(), timeGrid_[tidx])
146 : calibrationStrikes[j]);
147 d_variance[j] = std::max(tmp - variance[j], 1E-20);
148 variance[j] = tmp;
149 }
150 for (Size j = 0; j < indices_.size(); ++j) {
151 covariance_[i - 1][j][j] += d_variance[j];
152 for (Size k = 0; k < j; ++k) {
153 Real tmp = correlation[k][j] * std::sqrt(d_variance[j] * d_variance[k]);
154 covariance_[i - 1][k][j] += tmp;
155 covariance_[i - 1][j][k] += tmp;
156 }
157 }
158 ++tidx;
159 }
160
161 // salvage covariance matrix using spectral method if not positive semidefinite
162
163 SymmetricSchurDecomposition jd(covariance_[i - 1]);
164
165 bool needsSalvaging = false;
166 for (Size k = 0; k < covariance_[i - 1].rows(); ++k) {
167 if (jd.eigenvalues()[k] < -1E-16)
168 needsSalvaging = true;
169 }
170
171 if (needsSalvaging) {
172 Matrix diagonal(covariance_[i - 1].rows(), covariance_[i - 1].rows(), 0.0);
173 for (Size k = 0; k < jd.eigenvalues().size(); ++k) {
174 diagonal[k][k] = std::sqrt(std::max<Real>(jd.eigenvalues()[k], 0.0));
175 }
176 covariance_[i - 1] = jd.eigenvectors() * diagonal * diagonal * transpose(jd.eigenvectors());
177 }
178
179 // compute the _unique_ pos semidefinite square root
180
181 sqrtCov.push_back(CholeskyDecomposition(covariance_[i - 1], true));
182
183 // drift
184
185 Date d = *std::next(effectiveSimulationDates_.begin(), i);
186 for (Size j = 0; j < indices_.size(); ++j) {
187 Real tmp = model_->processes()[j]->riskFreeRate()->discount(d) /
188 model_->processes()[j]->dividendYield()->discount(d);
189 drift[i - 1][j] = -std::log(tmp / discountRatio[j]) - 0.5 * covariance_[i - 1][j][j];
190 discountRatio[j] = tmp;
191
192 // drift adjustment for eq / com indices that are not in base ccy
193 if (forCcyDaIndex[j] != Null<Size>()) {
194 drift[i - 1][j] -= covariance_[i - 1][forCcyDaIndex[j]][j];
195 }
196 }
197 }
198
199 // evolve the process using correlated normal variates and set the underlying path values
200
205 drift, sqrtCov);
206
207 if (trainingSamples() != Null<Size>()) {
212 drift, sqrtCov);
213 }
214
215 // set additional results provided by this model
216
217 for (Size i = 0; i < indices_.size(); ++i) {
218 for (Size j = 0; j < i; ++j) {
219 additionalResults_["BlackScholes.Correlation_" + indices_[i].name() + "_" + indices_[j].name()] =
220 correlation(i, j);
221 }
222 }
223
224 for (Size i = 0; i < calibrationStrikes.size(); ++i) {
225 additionalResults_["BlackScholes.CalibrationStrike_" + indices_[i].name()] =
226 (calibrationStrikes[i] == Null<Real>() ? "ATMF" : std::to_string(calibrationStrikes[i]));
227 }
228
229 for (Size i = 0; i < indices_.size(); ++i) {
230 Size timeStep = 0;
231 for (auto const& d : effectiveSimulationDates_) {
232 Real t = timeGrid_[positionInTimeGrid_[timeStep]];
233 Real forward = atmForward(model_->processes()[i]->x0(), model_->processes()[i]->riskFreeRate(),
234 model_->processes()[i]->dividendYield(), t);
235 if (timeStep > 0) {
236 Real volatility = model_->processes()[i]->blackVolatility()->blackVol(
237 t, calibrationStrikes[i] == Null<Real>() ? forward : calibrationStrikes[i]);
238 additionalResults_["BlackScholes.Volatility_" + indices_[i].name() + "_" + ore::data::to_string(d)] =
239 volatility;
240 }
241 additionalResults_["BlackScholes.Forward_" + indices_[i].name() + "_" + ore::data::to_string(d)] = forward;
242 ++timeStep;
243 }
244 }
245}
246
247void BlackScholes::populatePathValues(const Size nSamples, std::map<Date, std::vector<RandomVariable>>& paths,
248 const QuantLib::ext::shared_ptr<MultiPathVariateGeneratorBase>& gen,
249 const std::vector<Array>& drift, const std::vector<Matrix>& sqrtCov) const {
250
251 std::vector<std::vector<RandomVariable*>> rvs(indices_.size(),
252 std::vector<RandomVariable*>(effectiveSimulationDates_.size() - 1));
253 auto date = effectiveSimulationDates_.begin();
254 for (Size i = 0; i < effectiveSimulationDates_.size() - 1; ++i) {
255 ++date;
256 for (Size j = 0; j < indices_.size(); ++j) {
257 rvs[j][i] = &paths[*date][j];
258 rvs[j][i]->expand();
259 }
260 }
261
262 Array logState(indices_.size()), logState0(indices_.size());
263 for (Size j = 0; j < indices_.size(); ++j) {
264 logState0[j] = std::log(model_->processes()[j]->x0());
265 }
266
267 for (Size path = 0; path < nSamples; ++path) {
268 auto seq = gen->next();
269 logState = logState0;
270 for (Size i = 0; i < effectiveSimulationDates_.size() - 1; ++i) {
271 for (Size j = 0; j < indices_.size(); ++j) {
272 for (Size k = 0; k < indices_.size(); ++k)
273 logState[j] += sqrtCov[i][j][k] * seq.value[i][k];
274 }
275 logState += drift[i];
276 for (Size j = 0; j < indices_.size(); ++j)
277 rvs[j][i]->data()[path] = std::exp(logState[j]);
278 }
279 }
280} // populatePathValues()
281
282namespace {
283struct comp {
284 comp(const std::string& indexInput) : indexInput_(indexInput) {}
285 template <typename T> bool operator()(const std::pair<IndexInfo, QuantLib::ext::shared_ptr<T>>& p) const {
286 return p.first.name() == indexInput_;
287 }
288 const std::string indexInput_;
289};
290} // namespace
291
292RandomVariable BlackScholes::getFutureBarrierProb(const std::string& index, const Date& obsdate1, const Date& obsdate2,
293 const RandomVariable& barrier, const bool above) const {
294
295 // get the underlying values at the start and end points of the period
296
297 RandomVariable v1 = eval(index, obsdate1, Null<Date>());
298 RandomVariable v2 = eval(index, obsdate2, Null<Date>());
299
300 // check the barrier at the two endpoints
301
302 Filter barrierHit(barrier.size(), false);
303 if (above) {
304 barrierHit = barrierHit || v1 >= barrier;
305 barrierHit = barrierHit || v2 >= barrier;
306 } else {
307 barrierHit = barrierHit || v1 <= barrier;
308 barrierHit = barrierHit || v2 <= barrier;
309 }
310
311 // for IR / INF indices we have to check each date, this is fast though, since the index values are deteministic
312 // here
313
314 auto ir = std::find_if(irIndices_.begin(), irIndices_.end(), comp(index));
315 auto inf = std::find_if(infIndices_.begin(), infIndices_.end(), comp(index));
316 if (ir != irIndices_.end() || inf != infIndices_.end()) {
317 Date d = obsdate1 + 1;
318 while (d < obsdate2) {
319 RandomVariable res;
320 if (ir != irIndices_.end())
321 res = getIrIndexValue(std::distance(irIndices_.begin(), ir), d);
322 else
323 res = getInfIndexValue(std::distance(infIndices_.begin(), inf), d);
324 if (res.initialised()) {
325 if (above) {
326 barrierHit = barrierHit || res >= barrier;
327 } else {
328 barrierHit = barrierHit || res <= barrier;
329 }
330 }
331 ++d;
332 }
333 }
334
335 // Convert the bool result we have obtained so far to 1 / 0
336
337 RandomVariable result(barrierHit, 1.0, 0.0);
338
339 // return the result if we have an IR / INF index
340
341 if (ir != irIndices_.end() || inf != infIndices_.end())
342 return result;
343
344 // Now handle the dynamic indices in this model
345
346 // first make sure that v1 is not a historical fixing to avoid differences to a daily MC simulation where
347 // the process starts to evolve at the market data spot (only matters if the two values differ, which they
348 // should not too much anyway)
349
350 if (obsdate1 == referenceDate()) {
351 v1 = eval(index, obsdate1, Null<Date>(), false, true);
352 }
353
354 IndexInfo indexInfo(index);
355
356 // if we have an FX index, we "normalise" the tag by GENERIC (since it does not matter for projections), this
357 // is the same as in ModelImpl::eval()
358
359 if (indexInfo.isFx())
360 indexInfo = IndexInfo("FX-GENERIC-" + indexInfo.fx()->sourceCurrency().code() + "-" +
361 indexInfo.fx()->targetCurrency().code());
362
363 // get the average volatility over the period for the underlying, this is an approximation that could be refined
364 // by sampling more points in the interval and doing the below calculation for each subinterval - this would
365 // be slower though, so probably in the end we want this to be configurable in the model settings
366
367 // We might have one or two indices contributing to the desired volatility, since FX indices might require a
368 // triangulation. We look for the indices ind1 and ind2 so that the index is the quotient of the two.
369
370 Size ind1 = Null<Size>(), ind2 = Null<Size>();
371
372 auto i = std::find(indices_.begin(), indices_.end(), indexInfo);
373 if (i != indices_.end()) {
374 // we find the index directly in the index list
375 ind1 = std::distance(indices_.begin(), i);
376 } else {
377 // we can maybe triangulate an FX index (again, this is similar to ModelImpl::eval())
378 // if not, we can only try something else for FX indices
379 QL_REQUIRE(indexInfo.isFx(), "BlackScholes::getFutureBarrierProb(): index " << index << " not handled");
380 // is it a trivial fx index (CCY-CCY) or can we triangulate it?
381 if (indexInfo.fx()->sourceCurrency() == indexInfo.fx()->targetCurrency()) {
382 // pseudo FX index FX-GENERIC-CCY-CCY, leave ind1 and ind2 both set to zero
383 } else {
384 // if not, we can only try something else for FX indices
385 QL_REQUIRE(indexInfo.isFx(), "BlackScholes::getFutureBarrierProb(): index " << index << " not handled");
386 if (indexInfo.fx()->sourceCurrency() == indexInfo.fx()->targetCurrency()) {
387 // nothing to do, leave ind1 and ind2 = null
388 } else {
389 for (Size i = 0; i < indexCurrencies_.size(); ++i) {
390 if (indices_[i].isFx()) {
391 if (indexInfo.fx()->sourceCurrency().code() == indexCurrencies_[i])
392 ind1 = i;
393 if (indexInfo.fx()->targetCurrency().code() == indexCurrencies_[i])
394 ind2 = i;
395 }
396 }
397 }
398 }
399 }
400
401 // accumulate the variance contributions over [obsdate1, obsdate2]
402
403 Real variance = 0.0;
404
405 for (Size i = 1; i < effectiveSimulationDates_.size(); ++i) {
406
407 Date d1 = *std::next(effectiveSimulationDates_.begin(), i - 1);
408 Date d2 = *std::next(effectiveSimulationDates_.begin(), i);
409
410 if (obsdate1 <= d1 && d2 <= obsdate2) {
411 if (ind1 != Null<Size>()) {
412 variance += covariance_[i - 1][ind1][ind1];
413 }
414 if (ind2 != Null<Size>()) {
415 variance += covariance_[i - 1][ind2][ind2];
416 }
417 if (ind1 != Null<Size>() && ind2 != Null<Size>()) {
418 variance -= 2.0 * covariance_[i - 1][ind1][ind2];
419 }
420 }
421 }
422
423 // compute the hit probability
424 // see e.g. formulas 2, 4 in Emmanuel Gobet, Advanced Monte Carlo methods for barrier and related exotic options
425
426 if (!QuantLib::close_enough(variance, 0.0)) {
427 RandomVariable eps = RandomVariable(barrier.size(), 1E-14);
428 RandomVariable hitProb = exp(RandomVariable(barrier.size(), -2.0 / variance) * log(v1 / max(barrier, eps)) *
429 log(v2 / max(barrier, eps)));
430 result = result + applyInverseFilter(hitProb, barrierHit);
431 }
432
433 return result;
434}
435
436} // namespace data
437} // namespace ore
const std::string indexInput_
black scholes model for n underlyings (fx, equity or commodity)
const std::vector< Real > calibrationStrikes_
RandomVariable getInfIndexValue(const Size indexNo, const Date &d, const Date &fwd=Null< Date >()) const override
void performCalculations() const override
Size size() const override
RandomVariable getIrIndexValue(const Size indexNo, const Date &d, const Date &fwd=Null< Date >()) const override
const Date & referenceDate() const override
std::map< Date, std::vector< RandomVariable > > underlyingPathsTraining_
std::map< Date, std::vector< RandomVariable > > underlyingPaths_
std::set< Date > effectiveSimulationDates_
const Handle< BlackScholesModelWrapper > model_
Size trainingSamples() const override
std::vector< Size > positionInTimeGrid_
void performCalculations() const override
RandomVariable getFutureBarrierProb(const std::string &index, const Date &obsdate1, const Date &obsdate2, const RandomVariable &barrier, const bool above) const override
void populatePathValues(const Size nSamples, std::map< Date, std::vector< RandomVariable > > &paths, const QuantLib::ext::shared_ptr< MultiPathVariateGeneratorBase > &gen, const std::vector< Array > &drift, const std::vector< Matrix > &sqrtCov) const
std::vector< Matrix > covariance_
const std::map< std::string, std::vector< Real > > calibrationStrikes_
const std::string calibration_
BlackScholes(const Size paths, const std::vector< std::string > &currencies, const std::vector< Handle< YieldTermStructure > > &curves, const std::vector< Handle< Quote > > &fxSpots, const std::vector< std::pair< std::string, QuantLib::ext::shared_ptr< InterestRateIndex > > > &irIndices, const std::vector< std::pair< std::string, QuantLib::ext::shared_ptr< ZeroInflationIndex > > > &infIndices, const std::vector< std::string > &indices, const std::vector< std::string > &indexCurrencies, const Handle< BlackScholesModelWrapper > &model, const std::map< std::pair< std::string, std::string >, Handle< QuantExt::CorrelationTermStructure > > &correlations, const McParams &mcParams, const std::set< Date > &simulationDates, const IborFallbackConfig &iborFallbackConfig=IborFallbackConfig::defaultConfig(), const std::string &calibration="ATM", const std::map< std::string, std::vector< Real > > &calibrationStrikes={})
bool isFx() const
Definition: utilities.hpp:100
QuantLib::ext::shared_ptr< FxIndex > fx() const
Definition: utilities.hpp:109
std::map< std::string, boost::any > additionalResults_
Definition: model.hpp:153
std::vector< std::pair< IndexInfo, QuantLib::ext::shared_ptr< ZeroInflationIndex > > > infIndices_
Definition: modelimpl.hpp:113
std::vector< IndexInfo > indices_
Definition: modelimpl.hpp:114
RandomVariable eval(const std::string &index, const Date &obsdate, const Date &fwddate, const bool returnMissingMissingAsNull=false, const bool ignoreTodaysFixing=false) const override
Definition: modelimpl.cpp:181
std::vector< std::pair< IndexInfo, QuantLib::ext::shared_ptr< InterestRateIndex > > > irIndices_
Definition: modelimpl.hpp:112
const std::vector< std::string > indexCurrencies_
Definition: modelimpl.hpp:108
@ data
Definition: log.hpp:77
#define TLOG(text)
Logging Macro (Level = Data)
Definition: log.hpp:556
Shared utilities for model building and calibration.
RandomVariable max(RandomVariable x, const RandomVariable &y)
RandomVariable exp(RandomVariable x)
RandomVariable variance(const RandomVariable &r)
RandomVariable log(RandomVariable x)
boost::shared_ptr< MultiPathVariateGeneratorBase > makeMultiPathVariateGenerator(const SequenceType s, const Size dimension, const Size timeSteps, const BigNatural seed, const SobolBrownianGenerator::Ordering ordering, const SobolRsg::DirectionIntegers directionIntegers)
RandomVariable applyInverseFilter(RandomVariable x, const Filter &f)
Real atmForward(const Real s0, const Handle< YieldTermStructure > &r, const Handle< YieldTermStructure > &q, const Real t)
helper function that computes the atm forward
Definition: utilities.cpp:483
std::string to_string(const LocationInfo &l)
Definition: ast.cpp:28
Serializable Credit Default Swap.
Definition: namespaces.docs:23
QuantExt::SequenceType sequenceType
Definition: model.hpp:54
QuantLib::SobolBrownianGenerator::Ordering sobolOrdering
Definition: model.hpp:59
QuantExt::SequenceType trainingSequenceType
Definition: model.hpp:55
QuantLib::SobolRsg::DirectionIntegers sobolDirectionIntegers
Definition: model.hpp:60
string conversion utilities
string name