1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3 /*
4 Copyright (C) 2006 Ferdinando Ametrano
5 Copyright (C) 2006 Marco Bianchetti
6 Copyright (C) 2006 Cristina Duminuco
7 Copyright (C) 2006 StatPro Italia srl
8 Copyright (C) 2008 Mark Joshi
9 Copyright (C) 2012 Peter Caspers
10
11 This file is part of QuantLib, a free-software/open-source library
12 for financial quantitative analysts and developers - http://quantlib.org/
13
14 QuantLib is free software: you can redistribute it and/or modify it
15 under the terms of the QuantLib license. You should have received a
16 copy of the license along with this program; if not, please email
17 <quantlib-dev@lists.sf.net>. The license is also available online at
18 <http://quantlib.org/license.shtml>.
19
20 This program is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
22 FOR A PARTICULAR PURPOSE. See the license for more details.
23 */
24
25 #include "marketmodel.hpp"
26 #include "utilities.hpp"
27 #include <ql/models/marketmodels/accountingengine.hpp>
28 #include <ql/models/marketmodels/browniangenerators/mtbrowniangenerator.hpp>
29 #include <ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp>
30 #include <ql/models/marketmodels/callability/collectnodedata.hpp>
31 #include <ql/models/marketmodels/callability/lsstrategy.hpp>
32 #include <ql/models/marketmodels/callability/nothingexercisevalue.hpp>
33 #include <ql/models/marketmodels/callability/parametricexerciseadapter.hpp>
34 #include <ql/models/marketmodels/callability/swapbasissystem.hpp>
35 #include <ql/models/marketmodels/callability/swapratetrigger.hpp>
36 #include <ql/models/marketmodels/callability/triggeredswapexercise.hpp>
37 #include <ql/models/marketmodels/callability/upperboundengine.hpp>
38 #include <ql/models/marketmodels/curvestates/lmmcurvestate.hpp>
39 #include <ql/models/marketmodels/driftcomputation/lmmdriftcalculator.hpp>
40 #include <ql/models/marketmodels/evolvers/lognormalfwdrateeuler.hpp>
41 #include <ql/models/marketmodels/evolvers/lognormalfwdrateeulerconstrained.hpp>
42 #include <ql/models/marketmodels/evolvers/lognormalfwdrateipc.hpp>
43 #include <ql/models/marketmodels/evolvers/lognormalfwdrateballand.hpp>
44 #include <ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp>
45 #include <ql/models/marketmodels/evolvers/normalfwdratepc.hpp>
46 #include <ql/models/marketmodels/discounter.hpp>
47 #include <ql/models/marketmodels/models/abcdvol.hpp>
48 #include <ql/models/marketmodels/models/flatvol.hpp>
49 #include <ql/models/marketmodels/correlations/expcorrelations.hpp>
50 #include <ql/models/marketmodels/correlations/timehomogeneousforwardcorrelation.hpp>
51 #include <ql/models/marketmodels/products/multiproductcomposite.hpp>
52 #include <ql/models/marketmodels/products/multistep/callspecifiedmultiproduct.hpp>
53 #include <ql/models/marketmodels/products/multistep/exerciseadapter.hpp>
54 #include <ql/models/marketmodels/products/multistep/multistepcoinitialswaps.hpp>
55 #include <ql/models/marketmodels/products/multistep/multistepcoterminalswaps.hpp>
56 #include <ql/models/marketmodels/products/multistep/multistepcoterminalswaptions.hpp>
57 #include <ql/models/marketmodels/products/multistep/multistepperiodcapletswaptions.hpp>
58 #include <ql/models/marketmodels/products/multistep/multistepforwards.hpp>
59 #include <ql/models/marketmodels/products/multistep/multistepnothing.hpp>
60 #include <ql/models/marketmodels/products/multistep/multistepoptionlets.hpp>
61 #include <ql/models/marketmodels/products/multistep/multistepswap.hpp>
62 #include <ql/models/marketmodels/products/onestep/onestepforwards.hpp>
63 #include <ql/models/marketmodels/products/onestep/onestepoptionlets.hpp>
64 #include <ql/models/marketmodels/forwardforwardmappings.hpp>
65 #include <ql/models/marketmodels/proxygreekengine.hpp>
66 #include <ql/models/marketmodels/swapforwardmappings.hpp>
67 #include <ql/models/marketmodels/models/fwdperiodadapter.hpp>
68 #include <ql/models/marketmodels/models/fwdtocotswapadapter.hpp>
69 #include <ql/models/marketmodels/models/cotswaptofwdadapter.hpp>
70 #include <ql/models/marketmodels/utilities.hpp>
71 #include <ql/methods/montecarlo/genericlsregression.hpp>
72 #include <ql/legacy/libormarketmodels/lmlinexpcorrmodel.hpp>
73 #include <ql/legacy/libormarketmodels/lmextlinexpvolmodel.hpp>
74 #include <ql/time/schedule.hpp>
75 #include <ql/time/calendars/nullcalendar.hpp>
76 #include <ql/time/daycounters/simpledaycounter.hpp>
77 #include <ql/pricingengines/blackformula.hpp>
78 #include <ql/pricingengines/blackcalculator.hpp>
79 #include <ql/utilities/dataformatters.hpp>
80 #include <ql/math/integrals/segmentintegral.hpp>
81 #include <ql/math/statistics/convergencestatistics.hpp>
82 #include <ql/termstructures/volatility/abcd.hpp>
83 #include <ql/termstructures/volatility/abcdcalibration.hpp>
84 #include <ql/math/functional.hpp>
85 #include <ql/math/optimization/simplex.hpp>
86 #include <ql/quotes/simplequote.hpp>
87 #include <ql/auto_ptr.hpp>
88
89 #include <ql/models/marketmodels/products/pathwise/pathwiseproductcaplet.hpp>
90 #include <ql/models/marketmodels/products/pathwise/pathwiseproductswaption.hpp>
91
92 #include <ql/models/marketmodels/pathwiseaccountingengine.hpp>
93 #include <ql/models/marketmodels/pathwisegreeks/ratepseudorootjacobian.hpp>
94 #include <ql/models/marketmodels/pathwisegreeks/swaptionpseudojacobian.hpp>
95
96 #include <ql/models/marketmodels/models/pseudorootfacade.hpp>
97
98 #include <ql/models/marketmodels/pathwisegreeks/bumpinstrumentjacobian.hpp>
99
100 #include <ql/models/marketmodels/evolvers/volprocesses/squarerootandersen.hpp>
101
102 #include <ql/models/marketmodels/evolvers/svddfwdratepc.hpp>
103 #include <ql/processes/hestonprocess.hpp>
104 #include <ql/models/equity/hestonmodel.hpp>
105 #include <ql/time/daycounters/actualactual.hpp>
106 #include <ql/pricingengines/vanilla/analytichestonengine.hpp>
107
108 #include <ql/models/marketmodels/products/multistep/multistepinversefloater.hpp>
109 #include <ql/models/marketmodels/products/pathwise/pathwiseproductinversefloater.hpp>
110 #include <ql/models/marketmodels/products/multistep/multisteppathwisewrapper.hpp>
111
112 #include <boost/math/special_functions/fpclassify.hpp>
113 #include <boost/preprocessor/iteration/local.hpp>
114 #include <ql/functional.hpp>
115 #include <sstream>
116
117 using namespace QuantLib;
118 using namespace boost::unit_test_framework;
119
120 using std::fabs;
121 using std::sqrt;
122
123 #define BEGIN(x) (x+0)
124 #define END(x) (x+LENGTH(x))
125
126 namespace market_model_test {
127
128 Date todaysDate, startDate, endDate;
129 Schedule dates;
130 std::vector<Time> rateTimes, paymentTimes;
131 std::vector<Real> accruals;
132 Calendar calendar;
133 DayCounter dayCounter;
134 std::vector<Rate> todaysForwards, todaysCoterminalSwapRates;
135 Rate meanForward;
136 std::vector<Real> coterminalAnnuity;
137 Spread displacement;
138 std::vector<DiscountFactor> todaysDiscounts;
139 std::vector<Volatility> volatilities, blackVols, normalVols;
140 std::vector<Volatility> swaptionsVolatilities, swaptionsBlackVols;
141 Real a, b, c, d;
142 Real longTermCorrelation, beta;
143 Size measureOffset_;
144 unsigned long seed_;
145 Size paths_, trainingPaths_;
146 bool printReport_ = false;
147
148
149 // a simple structure to store some data which will be used during tests
150 struct SubProductExpectedValues {
SubProductExpectedValuesmarket_model_test::SubProductExpectedValues151 explicit SubProductExpectedValues(const std::string& descr)
152 : description(descr), testBias(false) {}
153 std::string description;
154 std::vector<Real> values;
155 bool testBias;
156 Real errorThreshold;
157 };
158
setup()159 void setup() {
160
161 // Times
162 calendar = NullCalendar();
163 todaysDate = Settings::instance().evaluationDate();
164 //startDate = todaysDate + 5*Years;
165 endDate = todaysDate + 5*Years;
166 dates =Schedule(todaysDate, endDate, Period(Semiannual),
167 calendar, Following, Following,
168 DateGeneration::Backward, false);
169 rateTimes = std::vector<Time>(dates.size()-1);
170 paymentTimes = std::vector<Time>(rateTimes.size()-1);
171 accruals = std::vector<Real>(rateTimes.size()-1);
172 dayCounter = SimpleDayCounter();
173 for (Size i=1; i<dates.size(); ++i)
174 rateTimes[i-1] = dayCounter.yearFraction(todaysDate, dates[i]);
175 std::copy(rateTimes.begin()+1, rateTimes.end(), paymentTimes.begin());
176 for (Size i=1; i<rateTimes.size(); ++i)
177 accruals[i-1] = rateTimes[i] - rateTimes[i-1];
178
179 // Rates & displacement
180 todaysForwards = std::vector<Rate>(paymentTimes.size());
181 displacement = 0.0;
182 meanForward=0.0;
183
184 for (Size i=0; i<todaysForwards.size(); ++i)
185 {
186 // FLOATING_POINT_EXCEPTION
187 todaysForwards[i] = 0.03 + 0.0010*i;
188 meanForward+= todaysForwards[i];
189 }
190 meanForward /= todaysForwards.size();
191
192
193
194 // Discounts
195 todaysDiscounts = std::vector<DiscountFactor>(rateTimes.size());
196 todaysDiscounts[0] = 0.95;
197 for (Size i=1; i<rateTimes.size(); ++i)
198 todaysDiscounts[i] = todaysDiscounts[i-1] /
199 (1.0+todaysForwards[i-1]*accruals[i-1]);
200
201 // Coterminal swap rates & annuities
202 Size N = todaysForwards.size();
203 todaysCoterminalSwapRates = std::vector<Rate>(N);
204 coterminalAnnuity = std::vector<Real>(N);
205 Real floatingLeg = 0.0;
206 for (Size i=1; i<=N; ++i) {
207 if (i==1) {
208 coterminalAnnuity[N-1] = accruals[N-1]*todaysDiscounts[N];
209 } else {
210 coterminalAnnuity[N-i] = coterminalAnnuity[N-i+1] +
211 accruals[N-i]*todaysDiscounts[N-i+1];
212 }
213 floatingLeg = todaysDiscounts[N-i]-todaysDiscounts[N];
214 todaysCoterminalSwapRates[N-i] = floatingLeg/coterminalAnnuity[N-i];
215 }
216
217 // Cap/Floor Volatilities
218 Volatility mktVols[] = {
219 0.15541283,
220 0.18719678,
221 0.20890740,
222 0.22318179,
223 0.23212717,
224 0.23731450,
225 0.23988649,
226 0.24066384,
227 0.24023111,
228 0.23900189,
229 0.23726699,
230 0.23522952,
231 0.23303022,
232 0.23076564,
233 0.22850101,
234 0.22627951,
235 0.22412881,
236 0.22206569,
237 0.22009939
238 };
239
240 a = -0.0597;
241 b = 0.1677;
242 c = 0.5403;
243 d = 0.1710;
244 volatilities = std::vector<Volatility>(todaysForwards.size());
245 blackVols = std::vector<Volatility>(todaysForwards.size());
246 normalVols = std::vector<Volatility>(todaysForwards.size());
247 for (Size i=0; i<std::min(LENGTH(mktVols),todaysForwards.size()); i++) {
248 volatilities[i] = todaysForwards[i]*mktVols[i]/
249 (todaysForwards[i]+displacement);
250 blackVols[i]= mktVols[i];
251 normalVols[i]= mktVols[i]*todaysForwards[i];
252 }
253
254 // Swaption volatility quick fix
255 swaptionsVolatilities = volatilities;
256
257 // Cap/Floor Correlation
258 longTermCorrelation = 0.5;
259 beta = 0.2;
260 measureOffset_ = 5;
261
262 // Monte Carlo
263 seed_ = 42;
264
265 #ifdef _DEBUG
266 paths_ = 127;
267 trainingPaths_ = 31;
268 #else
269 paths_ = 32767; //262144-1; //; // 2^15-1
270 trainingPaths_ = 8191; // 2^13-1
271 #endif
272 }
273
274 ext::shared_ptr<SequenceStatisticsInc>
simulate(const ext::shared_ptr<MarketModelEvolver> & evolver,const MarketModelMultiProduct & product)275 simulate(const ext::shared_ptr<MarketModelEvolver>& evolver,
276 const MarketModelMultiProduct& product) {
277 Size initialNumeraire = evolver->numeraires().front();
278 Real initialNumeraireValue = todaysDiscounts[initialNumeraire];
279
280 AccountingEngine engine(evolver, product, initialNumeraireValue);
281 ext::shared_ptr<SequenceStatisticsInc> stats(
282 new SequenceStatisticsInc(product.numberOfProducts()));
283 engine.multiplePathValues(*stats, paths_);
284 return stats;
285 }
286
287
marketModelTypeToString(MarketModelTest::MarketModelType type)288 std::string marketModelTypeToString(MarketModelTest::MarketModelType type) {
289 switch (type) {
290 case MarketModelTest::ExponentialCorrelationFlatVolatility:
291 return "Exp. Corr. Flat Vol.";
292 case MarketModelTest::ExponentialCorrelationAbcdVolatility:
293 return "Exp. Corr. Abcd Vol.";
294 //case CalibratedMM:
295 // return "CalibratedMarketModel";
296 default:
297 QL_FAIL("unknown MarketModelEvolver type");
298 }
299 }
300
makeMarketModel(bool logNormal,const EvolutionDescription & evolution,Size numberOfFactors,MarketModelTest::MarketModelType marketModelType,Spread forwardBump=0.0,Volatility volBump=0.0)301 ext::shared_ptr<MarketModel> makeMarketModel(
302 bool logNormal,
303 const EvolutionDescription& evolution,
304 Size numberOfFactors,
305 MarketModelTest::MarketModelType marketModelType,
306 Spread forwardBump = 0.0,
307 Volatility volBump = 0.0) {
308
309 std::vector<Time> fixingTimes(evolution.rateTimes());
310 fixingTimes.pop_back();
311 ext::shared_ptr<LmVolatilityModel> volModel(new
312 LmExtLinearExponentialVolModel(fixingTimes,0.5, 0.6, 0.1, 0.1));
313 ext::shared_ptr<LmCorrelationModel> corrModel(
314 new LmLinearExponentialCorrelationModel(evolution.numberOfRates(),
315 longTermCorrelation, beta));
316
317 std::vector<Rate> bumpedForwards(todaysForwards.size());
318 std::transform(todaysForwards.begin(), todaysForwards.end(),
319 bumpedForwards.begin(),
320 add<Rate>(forwardBump));
321
322 std::vector<Volatility> bumpedVols(volatilities.size());
323 if (logNormal)
324 std::transform(volatilities.begin(), volatilities.end(),
325 bumpedVols.begin(),
326 add<Volatility>(volBump));
327 else
328 std::transform(normalVols.begin(), normalVols.end(),
329 bumpedVols.begin(),
330 add<Volatility>(volBump));
331
332 Matrix correlations = exponentialCorrelations(evolution.rateTimes(),
333 longTermCorrelation,
334 beta);
335 ext::shared_ptr<PiecewiseConstantCorrelation> corr(new
336 TimeHomogeneousForwardCorrelation(correlations,
337 evolution.rateTimes()));
338 switch (marketModelType) {
339 case MarketModelTest::ExponentialCorrelationFlatVolatility:
340 return ext::shared_ptr<MarketModel>(new
341 FlatVol(bumpedVols,
342 corr,
343 evolution,
344 numberOfFactors,
345 bumpedForwards,
346 std::vector<Spread>(bumpedForwards.size(), displacement)));
347 case MarketModelTest::ExponentialCorrelationAbcdVolatility:
348 return ext::shared_ptr<MarketModel>(new
349 AbcdVol(0.0,0.0,1.0,1.0,
350 bumpedVols,
351 corr,
352 evolution,
353 numberOfFactors,
354 bumpedForwards,
355 std::vector<Spread>(bumpedForwards.size(), displacement)));
356 //case CalibratedMM:
357 // return ext::shared_ptr<MarketModel>(new
358 // CalibratedMarketModel(volModel, corrModel,
359 // evolution,
360 // numberOfFactors,
361 // bumpedForwards,
362 // displacement));
363 default:
364 QL_FAIL("unknown MarketModel type");
365 }
366 }
367
368 enum MeasureType { ProductSuggested, Terminal,
369 MoneyMarket, MoneyMarketPlus };
370
measureTypeToString(MeasureType type)371 std::string measureTypeToString(MeasureType type) {
372 switch (type) {
373 case ProductSuggested:
374 return "ProductSuggested measure";
375 case Terminal:
376 return "Terminal measure";
377 case MoneyMarket:
378 return "Money Market measure";
379 case MoneyMarketPlus:
380 return "Money Market Plus measure";
381 default:
382 QL_FAIL("unknown measure type");
383 }
384 }
385
makeMeasure(const MarketModelMultiProduct & product,MeasureType measureType)386 std::vector<Size> makeMeasure(const MarketModelMultiProduct& product,
387 MeasureType measureType) {
388 std::vector<Size> result;
389 const EvolutionDescription& evolution(product.evolution());
390 switch (measureType) {
391 case ProductSuggested:
392 result = product.suggestedNumeraires();
393 break;
394 case Terminal:
395 result = terminalMeasure(evolution);
396 if (!isInTerminalMeasure(evolution, result)) {
397 BOOST_ERROR("\nfailure in verifying Terminal measure:\n"
398 << to_stream(result));
399 }
400 break;
401 case MoneyMarket:
402 result = moneyMarketMeasure(evolution);
403 if (!isInMoneyMarketMeasure(evolution, result)) {
404 BOOST_ERROR("\nfailure in verifying MoneyMarket measure:\n"
405 << to_stream(result));
406 }
407 break;
408 case MoneyMarketPlus:
409 result = moneyMarketPlusMeasure(evolution, measureOffset_);
410 if (!isInMoneyMarketPlusMeasure(evolution, result, measureOffset_)) {
411 BOOST_ERROR("\nfailure in verifying MoneyMarketPlus(" <<
412 measureOffset_ << ") measure:\n" << to_stream(result));
413 }
414 break;
415 default:
416 QL_FAIL("unknown measure type");
417 }
418 checkCompatibility(evolution, result);
419 if (printReport_) {
420 BOOST_TEST_MESSAGE(" " << measureTypeToString(measureType) << ": " << to_stream(result));
421 }
422 return result;
423 }
424
425 enum EvolverType { Ipc, Balland, Pc, NormalPc};
426
evolverTypeToString(EvolverType type)427 std::string evolverTypeToString(EvolverType type) {
428 switch (type) {
429 case Ipc:
430 return "iterative predictor corrector";
431 case Balland:
432 return "Balland predictor corrector";
433 case Pc:
434 return "predictor corrector";
435 case NormalPc:
436 return "predictor corrector for normal case";
437 default:
438 QL_FAIL("unknown MarketModelEvolver type");
439 }
440 }
441
makeMarketModelEvolver(const ext::shared_ptr<MarketModel> & marketModel,const std::vector<Size> & numeraires,const BrownianGeneratorFactory & generatorFactory,EvolverType evolverType,Size initialStep=0)442 ext::shared_ptr<MarketModelEvolver> makeMarketModelEvolver(
443 const ext::shared_ptr<MarketModel>& marketModel,
444 const std::vector<Size>& numeraires,
445 const BrownianGeneratorFactory& generatorFactory,
446 EvolverType evolverType,
447 Size initialStep = 0) {
448 switch (evolverType) {
449 case Ipc:
450 return ext::shared_ptr<MarketModelEvolver>(
451 new LogNormalFwdRateIpc(marketModel, generatorFactory,
452 numeraires, initialStep));
453 case Balland:
454 return ext::shared_ptr<MarketModelEvolver>(
455 new LogNormalFwdRateBalland(marketModel, generatorFactory,
456 numeraires, initialStep));
457 case Pc:
458 return ext::shared_ptr<MarketModelEvolver>(
459 new LogNormalFwdRatePc(marketModel, generatorFactory,
460 numeraires, initialStep));
461 case NormalPc:
462 return ext::shared_ptr<MarketModelEvolver>(
463 new NormalFwdRatePc(marketModel, generatorFactory,
464 numeraires, initialStep));
465 default:
466 QL_FAIL("unknown MarketModelEvolver type");
467 }
468 }
469
checkMultiProductCompositeResults(const SequenceStatisticsInc & stats,const std::vector<SubProductExpectedValues> & subProductExpectedValues,const std::string & config)470 void checkMultiProductCompositeResults (const SequenceStatisticsInc& stats,
471 const std::vector<SubProductExpectedValues>& subProductExpectedValues,
472 const std::string& config) {
473
474 std::vector<Real> results = stats.mean();
475 std::vector<Real> errors = stats.errorEstimate();
476
477 // size check
478 Size nbOfResults = 0;
479 for (Size i=0; i<subProductExpectedValues.size(); ++i) {
480 for (Size j=0; j<subProductExpectedValues[i].values.size(); ++j)
481 ++nbOfResults;
482 }
483
484 if (nbOfResults != results.size())
485 BOOST_ERROR("mismatch between the size of the result and the \
486 number of results");
487 Size currentResultIndex = 0;
488
489 std::vector<Real> stdDevs;
490 std::vector<SubProductExpectedValues>::const_iterator subProductExpectedValue;
491 for (subProductExpectedValue = subProductExpectedValues.begin();
492 subProductExpectedValue != subProductExpectedValues.end();
493 ++subProductExpectedValue) {
494 Real minError = QL_MAX_REAL;
495 Real maxError = QL_MIN_REAL;
496 Real errorThreshold = subProductExpectedValue->errorThreshold;
497 for (Size j=0; j<subProductExpectedValue->values.size(); ++j) {
498 Real stdDev =
499 (results[currentResultIndex]-subProductExpectedValue->values[j])
500 /errors[currentResultIndex];
501 stdDevs.push_back(stdDev);
502 maxError = std::max(maxError, stdDev);
503 minError = std::min(minError, stdDev);
504 ++currentResultIndex;
505 }
506 bool isBiased = minError > 0.0 || maxError < 0.0;
507 if (printReport_
508 || (subProductExpectedValue->testBias && isBiased)
509 || std::max(-minError, maxError) > errorThreshold) {
510 BOOST_TEST_MESSAGE(config);
511 currentResultIndex = 0;
512 for (Size j=0; j<subProductExpectedValue->values.size(); ++j) {
513 BOOST_TEST_MESSAGE(io::ordinal(j+1)
514 << " " << subProductExpectedValue->description
515 << ": " << io::rate(results[currentResultIndex])
516 << "\t" << io::rate(subProductExpectedValue->values[j])
517 << "\t" << io::rate(errors[currentResultIndex])
518 << "; discrepancy = "
519 << stdDevs[currentResultIndex]
520 << "\n");
521 ++currentResultIndex;
522 }
523 BOOST_ERROR("test failed");
524 }
525 }
526 }
527
528
checkForwardsAndOptionlets(const SequenceStatisticsInc & stats,const std::vector<Rate> & forwardStrikes,const std::vector<ext::shared_ptr<StrikedTypePayoff>> & displacedPayoffs,const std::string & config)529 void checkForwardsAndOptionlets(const SequenceStatisticsInc& stats,
530 const std::vector<Rate>& forwardStrikes,
531 const std::vector<ext::shared_ptr<StrikedTypePayoff> >& displacedPayoffs,
532 const std::string& config) {
533
534 std::vector<Real> results = stats.mean();
535 std::vector<Real> errors = stats.errorEstimate();
536 std::vector<Real> stdDevs(todaysForwards.size());
537
538 Size N = todaysForwards.size();
539 std::vector<Rate> expectedForwards(N), expectedCaplets(N);
540 std::vector<Real> forwardStdDevs(N), capletStdDev(N);
541 Real minError = QL_MAX_REAL;
542 Real maxError = QL_MIN_REAL;
543 // forwards check
544 for (Size i=0; i<N; ++i) {
545 expectedForwards[i] = (todaysForwards[i]-forwardStrikes[i])
546 *accruals[i]*todaysDiscounts[i+1];
547 forwardStdDevs[i] = (results[i]-expectedForwards[i])/errors[i];
548 if (forwardStdDevs[i]>maxError)
549 maxError = forwardStdDevs[i];
550 else if (forwardStdDevs[i]<minError)
551 minError = forwardStdDevs[i];
552 Time expiry = rateTimes[i];
553 expectedCaplets[i] =
554 BlackCalculator(displacedPayoffs[i],
555 todaysForwards[i]+displacement,
556 volatilities[i]*std::sqrt(expiry),
557 todaysDiscounts[i+1]*accruals[i]).value();
558 capletStdDev[i] = (results[i+N]-expectedCaplets[i])/errors[i+N];
559 if (capletStdDev[i]>maxError)
560 maxError = capletStdDev[i];
561 else if (capletStdDev[i]<minError)
562 minError = capletStdDev[i];
563 }
564
565 Real errorThreshold = 2.50;
566 if ( printReport_ || minError > 0.0 || maxError < 0.0 ||
567 minError <-errorThreshold || maxError > errorThreshold) {
568 BOOST_TEST_MESSAGE(config);
569 Size i;
570 for (i=0; i<N; ++i) {
571 BOOST_TEST_MESSAGE(io::ordinal(i+1) << " forward: "
572 << io::rate(results[i])
573 << "\t" << io::rate(expectedForwards[i])
574 << "\t" << io::rate(errors[i])
575 << "; discrepancy = "
576 << forwardStdDevs[i]
577 << "\n");
578 }
579 for (i=0; i<N; ++i) {
580 BOOST_TEST_MESSAGE(
581 io::ordinal(i+1) << "\t"
582 << io::rate(results[i+N])
583 << " +- " << io::rate(errors[i+N])
584 << "\t" << io::rate(expectedCaplets[i])
585 << "\t" << io::rate(errors[i+N])
586 << "; discrepancy = "
587 << (results[i+N]-expectedCaplets[i])/(errors[i+N] == 0.0 ?
588 1.0 : errors[i+N])
589 << "\n");
590 }
591 BOOST_ERROR("test failed");
592 }
593 }
594
595
596
checkNormalForwardsAndOptionlets(const SequenceStatisticsInc & stats,const std::vector<Rate> & forwardStrikes,const std::vector<ext::shared_ptr<PlainVanillaPayoff>> & displacedPayoffs,const std::string & config)597 void checkNormalForwardsAndOptionlets(const SequenceStatisticsInc& stats,
598 const std::vector<Rate>& forwardStrikes,
599 const std::vector<ext::shared_ptr<PlainVanillaPayoff> >& displacedPayoffs,
600 const std::string& config) {
601
602 std::vector<Real> results = stats.mean();
603 std::vector<Real> errors = stats.errorEstimate();
604 std::vector<Real> stdDevs(todaysForwards.size());
605
606 Size N = todaysForwards.size();
607 std::vector<Rate> expectedForwards(N), expectedCaplets(N);
608 std::vector<Real> forwardStdDevs(N), capletStdDev(N);
609 Real minError = QL_MAX_REAL;
610 Real maxError = QL_MIN_REAL;
611 // forwards check
612 for (Size i=0; i<N; ++i) {
613 expectedForwards[i] = (todaysForwards[i]-forwardStrikes[i])
614 *accruals[i]*todaysDiscounts[i+1];
615 forwardStdDevs[i] = (results[i]-expectedForwards[i])/errors[i];
616 if (forwardStdDevs[i]>maxError)
617 maxError = forwardStdDevs[i];
618 else if (forwardStdDevs[i]<minError)
619 minError = forwardStdDevs[i];
620 Time expiry = rateTimes[i];
621 expectedCaplets[i] =
622 bachelierBlackFormula(displacedPayoffs[i],
623 todaysForwards[i]+displacement,
624 normalVols[i]*std::sqrt(expiry),
625 todaysDiscounts[i+1]*accruals[i]);
626 capletStdDev[i] = (results[i+N]-expectedCaplets[i])/errors[i+N];
627 if (capletStdDev[i]>maxError)
628 maxError = capletStdDev[i];
629 else if (capletStdDev[i]<minError)
630 minError = capletStdDev[i];
631 }
632
633 Real errorThreshold = 2.50;
634 if (minError > 0.0 || maxError < 0.0 ||
635 minError <-errorThreshold || maxError > errorThreshold) {
636 BOOST_TEST_MESSAGE(config);
637 Size i;
638 for (i=0; i<N; ++i) {
639 BOOST_TEST_MESSAGE(io::ordinal(i+1) << " forward: "
640 << io::rate(results[i])
641 << " +- " << io::rate(errors[i])
642 << "; expected: " << io::rate(expectedForwards[i])
643 << "; discrepancy = "
644 << forwardStdDevs[i]
645 << " standard errors");
646 }
647 for (i=0; i<N; ++i) {
648 BOOST_TEST_MESSAGE(
649 io::ordinal(i+1) << " caplet: "
650 << io::rate(results[i+N])
651 << " +- " << io::rate(errors[i+N])
652 << "; expected: " << io::rate(expectedCaplets[i])
653 << "; discrepancy = "
654 << (results[i+N]-expectedCaplets[i])/(errors[i+N] == 0.0 ?
655 1.0 : errors[i+N])
656 << " standard errors");
657 }
658 BOOST_ERROR("test failed");
659 }
660 }
661
662
663
checkCallableSwap(const SequenceStatisticsInc & stats,const std::string & config)664 void checkCallableSwap(const SequenceStatisticsInc& stats,
665 const std::string& config) {
666 Real payerNPV = stats.mean()[0];
667 Real receiverNPV = stats.mean()[1];
668 Real bermudanNPV = stats.mean()[2];
669 Real callableNPV = stats.mean()[3];
670 Real tolerance = 1.1e-15;
671 Real swapError = std::fabs(receiverNPV+payerNPV);
672 Real callableError = std::fabs(receiverNPV+bermudanNPV-callableNPV);
673
674 if (swapError>tolerance || bermudanNPV<0.0 ||
675 callableNPV<receiverNPV || callableError>tolerance)
676 BOOST_TEST_MESSAGE(config); // detailed error info below
677 if (swapError>tolerance)
678 BOOST_ERROR("agreement between payer and receiver swap failed:"
679 "\n payer swap: " << payerNPV <<
680 "\n receiver swap: " << receiverNPV <<
681 "\n error: " << swapError <<
682 "\n tolerance: " << tolerance);
683 if (bermudanNPV<0.0)
684 BOOST_ERROR("negative bermudan option value:"
685 "\n bermudan: " << bermudanNPV);
686 if (callableNPV<receiverNPV)
687 BOOST_ERROR("callable receiver less valuable than plain receiver:"
688 "\n receiver swap: " << receiverNPV <<
689 "\n callable: " << callableNPV);
690 if (callableError>tolerance)
691 BOOST_ERROR("agreement between receiver+bermudan and callable failed:"
692 "\n receiver swap: " << receiverNPV <<
693 "\n bermudan: " << bermudanNPV <<
694 "\n receiver+bermudan: " << receiverNPV+bermudanNPV <<
695 "\n callable: " << callableNPV <<
696 "\n error: " << callableError <<
697 "\n tolerance: " << tolerance);
698 if (printReport_) {
699 BOOST_TEST_MESSAGE(std::setprecision(2) <<
700 " payer swap: " << io::rate(payerNPV) << " +/- " << io::rate(stats.errorEstimate()[0]) <<
701 "\n receiver swap: " << io::rate(receiverNPV) << " +/- " << io::rate(stats.errorEstimate()[1]) <<
702 "\n bermudan: " << io::rate(bermudanNPV) << " +/- " << io::rate(stats.errorEstimate()[2]) <<
703 "\n receiver+bermudan: " << io::rate(receiverNPV+bermudanNPV) <<
704 "\n callable: " << io::rate(callableNPV) << " +/- " << io::rate(stats.errorEstimate()[3]));
705 }
706 }
707
708 }
709
710
testOneStepForwardsAndOptionlets()711 void MarketModelTest::testOneStepForwardsAndOptionlets() {
712
713 BOOST_TEST_MESSAGE("Testing exact repricing of "
714 "one-step forwards and optionlets "
715 "in a lognormal forward rate market model...");
716
717 using namespace market_model_test;
718
719 setup();
720
721 std::vector<Rate> forwardStrikes(todaysForwards.size());
722 std::vector<ext::shared_ptr<Payoff> > optionletPayoffs(todaysForwards.size());
723 std::vector<ext::shared_ptr<StrikedTypePayoff> >
724 displacedPayoffs(todaysForwards.size());
725 for (Size i=0; i<todaysForwards.size(); ++i) {
726 forwardStrikes[i] = todaysForwards[i] + 0.01;
727 optionletPayoffs[i] = ext::shared_ptr<Payoff>(new
728 PlainVanillaPayoff(Option::Call, todaysForwards[i]));
729 displacedPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
730 PlainVanillaPayoff(Option::Call, todaysForwards[i]+displacement));
731 }
732
733 OneStepForwards forwards(rateTimes, accruals,
734 paymentTimes, forwardStrikes);
735 OneStepOptionlets optionlets(rateTimes, accruals,
736 paymentTimes, optionletPayoffs);
737
738 MultiProductComposite product;
739 product.add(forwards);
740 product.add(optionlets);
741 product.finalize();
742
743 EvolutionDescription evolution = product.evolution();
744
745 MarketModelType marketModels[] = {
746 // CalibratedMM,
747 ExponentialCorrelationFlatVolatility,
748 ExponentialCorrelationAbcdVolatility };
749 for (Size j=0; j<LENGTH(marketModels); j++) {
750
751 // one step must be always full factors
752 Size testedFactors[] = { todaysForwards.size()};
753 for (Size m=0; m<LENGTH(testedFactors); ++m) {
754 Size factors = testedFactors[m];
755
756 // for one step product ProductSuggested is equal to Terminal
757 // for one step product MoneyMarketPlus is equal to Terminal
758 MeasureType measures[] = { MoneyMarket,
759 Terminal };
760 for (Size k=0; k<LENGTH(measures); k++) {
761 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
762
763 bool logNormal = true;
764 ext::shared_ptr<MarketModel> marketModel =
765 makeMarketModel(logNormal, evolution, factors, marketModels[j]);
766
767 EvolverType evolvers[] = { Pc, Balland, Ipc};
768 ext::shared_ptr<MarketModelEvolver> evolver;
769 Size stop =
770 isInTerminalMeasure(evolution, numeraires) ? 0 : 1;
771 for (Size i=0; i<LENGTH(evolvers)-stop; i++) {
772
773 for (Size n=0; n<1; n++) {
774 MTBrownianGeneratorFactory generatorFactory(seed_);
775 //SobolBrownianGeneratorFactory generatorFactory(
776 // SobolBrownianGenerator::Diagonal, seed_);
777
778 evolver = makeMarketModelEvolver(marketModel,
779 numeraires,
780 generatorFactory,
781 evolvers[i]);
782 std::ostringstream config;
783 config <<
784 marketModelTypeToString(marketModels[j]) << ", " <<
785 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
786 measureTypeToString(measures[k]) << ", " <<
787 evolverTypeToString(evolvers[i]) << ", " <<
788 "MT BGF";
789 if (printReport_)
790 BOOST_TEST_MESSAGE(" " << config.str());
791
792 ext::shared_ptr<SequenceStatisticsInc> stats =
793 simulate(evolver, product);
794 checkForwardsAndOptionlets(*stats,
795 forwardStrikes,
796 displacedPayoffs,
797 config.str());
798 }
799 }
800 }
801 }
802 }
803 }
804
testOneStepNormalForwardsAndOptionlets()805 void MarketModelTest::testOneStepNormalForwardsAndOptionlets() {
806
807 BOOST_TEST_MESSAGE("Testing exact repricing of "
808 "one-step forwards and optionlets "
809 "in a normal forward rate market model...");
810
811 using namespace market_model_test;
812
813 setup();
814
815 std::vector<Rate> forwardStrikes(todaysForwards.size());
816 std::vector<ext::shared_ptr<Payoff> > optionletPayoffs(todaysForwards.size());
817 std::vector<ext::shared_ptr<PlainVanillaPayoff> >
818 displacedPayoffs(todaysForwards.size());
819 for (Size i=0; i<todaysForwards.size(); ++i) {
820 forwardStrikes[i] = todaysForwards[i] + 0.01;
821 optionletPayoffs[i] = ext::shared_ptr<Payoff>(new
822 PlainVanillaPayoff(Option::Call, todaysForwards[i]));
823 displacedPayoffs[i] = ext::make_shared<PlainVanillaPayoff>(Option::Call, todaysForwards[i]+displacement);
824 }
825
826 OneStepForwards forwards(rateTimes, accruals,
827 paymentTimes, forwardStrikes);
828 OneStepOptionlets optionlets(rateTimes, accruals,
829 paymentTimes, optionletPayoffs);
830
831 MultiProductComposite product;
832 product.add(forwards);
833 product.add(optionlets);
834 product.finalize();
835
836 EvolutionDescription evolution = product.evolution();
837
838 MarketModelType marketModels[] = {
839 // CalibratedMM,
840 ExponentialCorrelationFlatVolatility,
841 ExponentialCorrelationAbcdVolatility };
842 for (Size j=0; j<LENGTH(marketModels); j++) {
843
844 // one step must be always full factors
845 Size testedFactors[] = { todaysForwards.size()};
846 for (Size m=0; m<LENGTH(testedFactors); ++m) {
847 Size factors = testedFactors[m];
848
849 // for one step product ProductSuggested is equal to Terminal
850 // for one step product MoneyMarketPlus is equal to Terminal
851 MeasureType measures[] = { MoneyMarket,
852 Terminal };
853 for (Size k=0; k<LENGTH(measures); k++) {
854 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
855
856 bool logNormal = false;
857 ext::shared_ptr<MarketModel> marketModel =
858 makeMarketModel(logNormal, evolution, factors, marketModels[j]);
859
860 EvolverType evolvers[] = { NormalPc};
861 ext::shared_ptr<MarketModelEvolver> evolver;
862 Size stop =
863 isInTerminalMeasure(evolution, numeraires) ? 0 : 1;
864 for (Size i=0; i<LENGTH(evolvers)-stop; i++) {
865
866 for (Size n=0; n<1; n++) {
867 MTBrownianGeneratorFactory generatorFactory(seed_);
868 //SobolBrownianGeneratorFactory generatorFactory(
869 // SobolBrownianGenerator::Diagonal, seed_);
870
871 evolver = makeMarketModelEvolver(marketModel,
872 numeraires,
873 generatorFactory,
874 evolvers[i]);
875 std::ostringstream config;
876 config <<
877 marketModelTypeToString(marketModels[j]) << ", " <<
878 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
879 measureTypeToString(measures[k]) << ", " <<
880 evolverTypeToString(evolvers[i]) << ", " <<
881 "MT BGF";
882 if (printReport_)
883 BOOST_TEST_MESSAGE(" " << config.str());
884
885 ext::shared_ptr<SequenceStatisticsInc> stats =
886 simulate(evolver, product);
887 checkNormalForwardsAndOptionlets(*stats,
888 forwardStrikes,
889 displacedPayoffs,
890 config.str());
891 }
892 }
893 }
894 }
895 }
896 }
897
testInverseFloater()898 void MarketModelTest::testInverseFloater()
899 {
900
901 BOOST_TEST_MESSAGE("Testing exact repricing of "
902 "inverse floater "
903 "in forward rate market model...");
904
905 using namespace market_model_test;
906
907 setup();
908
909
910 std::vector<Real> fixedStrikes(accruals.size(), 0.1);
911 std::vector<Real> fixedMultipliers(accruals.size(), 2.0);
912 std::vector<Real> floatingSpreads(accruals.size(), 0.002);
913 std::vector<Real> fixedAccruals(accruals);
914 std::vector<Real> floatingAccruals(accruals);
915
916 bool payer = true;
917
918
919 MultiStepInverseFloater product(
920 rateTimes,
921 fixedAccruals,
922 floatingAccruals,
923 fixedStrikes,
924 fixedMultipliers,
925 floatingSpreads,
926 paymentTimes,
927 payer);
928
929 MarketModelPathwiseInverseFloater productPathwise(
930 rateTimes,
931 fixedAccruals,
932 floatingAccruals,
933 fixedStrikes,
934 fixedMultipliers,
935 floatingSpreads,
936 paymentTimes,
937 payer);
938
939 MultiProductPathwiseWrapper productWrapped(productPathwise);
940
941 MultiProductComposite productComposite;
942 productComposite.add(product);
943 productComposite.add(productWrapped);
944 productComposite.finalize();
945
946
947
948
949 EvolutionDescription evolution = productComposite.evolution();
950
951 MarketModelType marketModels[] = {
952 // CalibratedMM,
953 ExponentialCorrelationFlatVolatility,
954 ExponentialCorrelationAbcdVolatility };
955 for (Size j=0; j<LENGTH(marketModels); j++)
956 {
957
958 Size testedFactors[] = { std::min<Size>(todaysForwards.size(),3)};
959 for (Size m=0; m<LENGTH(testedFactors); ++m)
960 {
961 Size factors = testedFactors[m];
962
963 MeasureType measures[] = { MoneyMarket};
964 for (Size k=0; k<LENGTH(measures); k++)
965 {
966 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
967
968 bool logNormal = false;
969 ext::shared_ptr<MarketModel> marketModel =
970 makeMarketModel(logNormal, evolution, factors, marketModels[j]);
971
972 EvolverType evolvers[] = {Pc};
973 ext::shared_ptr<MarketModelEvolver> evolver;
974
975 for (Size i=0; i<LENGTH(evolvers); i++)
976 {
977
978
979 MTBrownianGeneratorFactory generatorFactory(seed_);
980 //SobolBrownianGeneratorFactory generatorFactory(
981 // SobolBrownianGenerator::Diagonal, seed_);
982
983 evolver = makeMarketModelEvolver(marketModel,
984 numeraires,
985 generatorFactory,
986 evolvers[i]);
987 std::ostringstream config;
988 config <<
989 marketModelTypeToString(marketModels[j]) << ", " <<
990 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
991 measureTypeToString(measures[k]) << ", " <<
992 evolverTypeToString(evolvers[i]) << ", " <<
993 "MT BGF";
994 if (printReport_)
995 BOOST_TEST_MESSAGE(" " << config.str());
996
997 ext::shared_ptr<SequenceStatisticsInc> stats =
998 simulate(evolver, productComposite);
999
1000 std::vector<Real> modelVolatilities(accruals.size());
1001 for (Size i=0; i < accruals.size(); ++i)
1002 modelVolatilities[i] = sqrt(marketModel->totalCovariance(i)[i][i]);
1003
1004
1005
1006 Real truePrice =0.0;
1007
1008 for (Size i=0; i < accruals.size(); ++i)
1009 {
1010 Real floatingCouponPV = floatingAccruals[i] *(todaysForwards[i]+floatingSpreads[i])*todaysDiscounts[i+1];
1011 Real inverseCouponPV = 2*fixedAccruals[i] *todaysDiscounts[i+1]* blackFormula(Option::Put,
1012 fixedStrikes[i]/2.0,
1013 todaysForwards[i],
1014 modelVolatilities[i]);
1015
1016 truePrice += floatingCouponPV - inverseCouponPV;
1017 }
1018
1019
1020
1021
1022
1023 Real priceError = stats->mean()[0] - truePrice;
1024 Real priceSD = stats->errorEstimate()[0];
1025
1026 Real errorInSds = priceError/priceSD;
1027 if (fabs(errorInSds) > 4.0)
1028 BOOST_FAIL("Inverse floater product has price error equal to " <<errorInSds << " sds . Price " <<truePrice << " MC price " << stats->mean()[0] << " \n" );
1029
1030 Real numericalTolerance = 1E-12;
1031
1032 if (fabs(stats->mean()[0] - stats->mean()[1]) > numericalTolerance)
1033 BOOST_FAIL("Inverse floater and wrapper pathwise inverse floater do not agree:" << stats->mean()[0] << " " << stats->mean()[1] );
1034
1035
1036
1037
1038
1039 } // evolvers
1040 } // measures
1041 } // factors
1042 }
1043 }
1044
testMultiProductComposite(const MarketModelMultiProduct & product,const std::vector<market_model_test::SubProductExpectedValues> & subProductExpectedValues,const std::string & testDescription)1045 void testMultiProductComposite(const MarketModelMultiProduct& product,
1046 const std::vector<market_model_test::SubProductExpectedValues>& subProductExpectedValues,
1047 const std::string& testDescription) {
1048
1049 using namespace market_model_test;
1050
1051 BOOST_TEST_MESSAGE(
1052 "Testing exact repricing of "
1053 << testDescription
1054 << "in a lognormal forward rate market model...");
1055
1056 setup();
1057
1058 const EvolutionDescription& evolution = product.evolution();
1059
1060 MarketModelTest::MarketModelType marketModels[] = {
1061 // CalibratedMM,
1062 MarketModelTest::ExponentialCorrelationFlatVolatility,
1063 MarketModelTest::ExponentialCorrelationAbcdVolatility };
1064 for (Size j=0; j<LENGTH(marketModels); j++) {
1065
1066 Size testedFactors[] = { 4, 8,
1067 todaysForwards.size()};
1068 for (Size m=0; m<LENGTH(testedFactors); ++m) {
1069 Size factors = testedFactors[m];
1070
1071 // Composite's ProductSuggested is the Terminal one
1072 MeasureType measures[] = { // ProductSuggested,
1073 Terminal,
1074 MoneyMarketPlus,
1075 MoneyMarket};
1076 for (Size k=0; k<LENGTH(measures); k++) {
1077 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
1078
1079 bool logNormal = true;
1080 ext::shared_ptr<MarketModel> marketModel =
1081 makeMarketModel(logNormal, evolution, factors, marketModels[j]);
1082
1083
1084 EvolverType evolvers[] = { Pc, Balland, Ipc };
1085 ext::shared_ptr<MarketModelEvolver> evolver;
1086 Size stop =
1087 isInTerminalMeasure(evolution, numeraires) ? 0 : 1;
1088 for (Size i=0; i<LENGTH(evolvers)-stop; i++) {
1089
1090 for (Size n=0; n<1; n++) {
1091 //MTBrownianGeneratorFactory generatorFactory(seed_);
1092 SobolBrownianGeneratorFactory generatorFactory(
1093 SobolBrownianGenerator::Diagonal, seed_);
1094
1095 evolver = makeMarketModelEvolver(marketModel,
1096 numeraires,
1097 generatorFactory,
1098 evolvers[i]);
1099 std::ostringstream config;
1100 config <<
1101 marketModelTypeToString(marketModels[j]) << ", " <<
1102 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
1103 measureTypeToString(measures[k]) << ", " <<
1104 evolverTypeToString(evolvers[i]) << ", " <<
1105 "MT BGF";
1106 if (printReport_)
1107 BOOST_TEST_MESSAGE(" " << config.str());
1108
1109 ext::shared_ptr<SequenceStatisticsInc> stats =
1110 simulate(evolver, product);
1111 checkMultiProductCompositeResults(*stats,
1112 subProductExpectedValues,
1113 config.str());
1114 }
1115 }
1116 }
1117 }
1118 }
1119 }
1120
addForwards(MultiProductComposite & product,std::vector<market_model_test::SubProductExpectedValues> & subProductExpectedValues)1121 void addForwards(MultiProductComposite& product,
1122 std::vector<market_model_test::SubProductExpectedValues>& subProductExpectedValues) {
1123
1124 using namespace market_model_test;
1125
1126 // create forwards and add them to the product...
1127 std::vector<Rate> forwardStrikes(todaysForwards.size());
1128
1129 for (Size i=0; i<todaysForwards.size(); ++i)
1130 forwardStrikes[i] = todaysForwards[i] + 0.01;
1131
1132 MultiStepForwards forwards(rateTimes, accruals,
1133 paymentTimes, forwardStrikes);
1134 product.add(forwards);
1135
1136 // computing and storing expected values
1137 subProductExpectedValues.push_back(SubProductExpectedValues("Forward"));
1138 subProductExpectedValues.back().errorThreshold = 2.50;
1139 for (Size i=0; i<todaysForwards.size(); ++i) {
1140 subProductExpectedValues.back().values.push_back(
1141 (todaysForwards[i]-forwardStrikes[i])
1142 *accruals[i]*todaysDiscounts[i+1]);
1143 }
1144 }
1145
addOptionLets(MultiProductComposite & product,std::vector<market_model_test::SubProductExpectedValues> & subProductExpectedValues)1146 void addOptionLets(MultiProductComposite& product,
1147 std::vector<market_model_test::SubProductExpectedValues>& subProductExpectedValues) {
1148
1149 using namespace market_model_test;
1150
1151 // create the products...
1152 std::vector<ext::shared_ptr<Payoff> > optionletPayoffs(todaysForwards.size());
1153 std::vector<ext::shared_ptr<StrikedTypePayoff> >
1154 displacedPayoffs(todaysForwards.size());
1155
1156 for (Size i=0; i<todaysForwards.size(); ++i) {
1157 optionletPayoffs[i] = ext::shared_ptr<Payoff>(new
1158 PlainVanillaPayoff(Option::Call, todaysForwards[i]));
1159 //CashOrNothingPayoff(Option::Call, todaysForwards[i], 0.01));
1160 displacedPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
1161 PlainVanillaPayoff(Option::Call, todaysForwards[i]+displacement));
1162 //CashOrNothingPayoff(Option::Call, todaysForwards[i]+displacement, 0.01));
1163 }
1164
1165 MultiStepOptionlets optionlets(rateTimes, accruals,
1166 paymentTimes, optionletPayoffs);
1167 product.add(optionlets);
1168
1169 // computing and storing expected values
1170 subProductExpectedValues.push_back(SubProductExpectedValues("Caplet"));
1171 subProductExpectedValues.back().errorThreshold = 2.50;
1172 for (Size i=0; i<todaysForwards.size(); ++i) {
1173 subProductExpectedValues.back().values.push_back(
1174 BlackCalculator(displacedPayoffs[i],
1175 todaysForwards[i]+displacement,
1176 volatilities[i]*std::sqrt(rateTimes[i]),
1177 todaysDiscounts[i+1]*accruals[i]).value());
1178 }
1179 }
1180
1181
addCoinitialSwaps(MultiProductComposite & product,std::vector<market_model_test::SubProductExpectedValues> & subProductExpectedValues)1182 void addCoinitialSwaps(MultiProductComposite& product,
1183 std::vector<market_model_test::SubProductExpectedValues>& subProductExpectedValues) {
1184
1185 using namespace market_model_test;
1186
1187 // create the products...
1188 Real fixedRate = 0.04;
1189 MultiStepCoinitialSwaps multiStepCoinitialSwaps(rateTimes, accruals, accruals,
1190 paymentTimes, fixedRate);
1191 product.add(multiStepCoinitialSwaps);
1192 // computing and storing expected values
1193 subProductExpectedValues.push_back(SubProductExpectedValues("coinitial swap"));
1194 subProductExpectedValues.back().testBias = false;
1195 subProductExpectedValues.back().errorThreshold = 2.32;
1196 Real coinitialSwapValue = 0;
1197 for (Size i=0; i<todaysForwards.size(); ++i) {
1198 coinitialSwapValue += (todaysForwards[i]-fixedRate)
1199 *accruals[i]*todaysDiscounts[i+1];
1200 subProductExpectedValues.back().values.push_back(coinitialSwapValue);
1201 }
1202 }
1203
addCoterminalSwapsAndSwaptions(MultiProductComposite & product,std::vector<market_model_test::SubProductExpectedValues> & subProductExpectedValues)1204 void addCoterminalSwapsAndSwaptions(MultiProductComposite& product,
1205 std::vector<market_model_test::SubProductExpectedValues>& subProductExpectedValues) {
1206
1207 using namespace market_model_test;
1208
1209 Real fixedRate = 0.04;
1210 MultiStepCoterminalSwaps swaps(rateTimes, accruals, accruals,
1211 paymentTimes, fixedRate);
1212
1213 std::vector<ext::shared_ptr<StrikedTypePayoff> > payoffs(todaysForwards.size());
1214 for (Size i = 0; i < payoffs.size(); ++i)
1215 payoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
1216 PlainVanillaPayoff(Option::Call, todaysForwards[i]));
1217
1218 MultiStepCoterminalSwaptions swaptions(rateTimes,
1219 rateTimes, payoffs);
1220 product.add(swaps);
1221 product.add(swaptions);
1222
1223 subProductExpectedValues.push_back(SubProductExpectedValues("coterminal swap"));
1224 subProductExpectedValues.back().testBias = false;
1225 subProductExpectedValues.back().errorThreshold = 2.32;
1226 LMMCurveState curveState(rateTimes); // not the best way to detect errors in LMMCurveState...
1227 curveState.setOnForwardRates(todaysForwards);
1228 std::vector<Rate> atmRates = curveState.coterminalSwapRates();
1229 for (Size i=0; i<todaysForwards.size(); ++i) {
1230 Real expectedNPV = curveState.coterminalSwapAnnuity(0, i) * (atmRates[i]-fixedRate) *
1231 todaysDiscounts.front();
1232 subProductExpectedValues.back().values.push_back(expectedNPV);
1233 }
1234 // we clone the prooduct to be able to finalize it and call evolution function member on it
1235 MultiProductComposite productClone = product;
1236 productClone.finalize();
1237 subProductExpectedValues.push_back(SubProductExpectedValues("coterminal swaption"));
1238 subProductExpectedValues.back().testBias = false;
1239 subProductExpectedValues.back().errorThreshold = 2.32;
1240 const Spread displacement = 0;
1241 Matrix jacobian =
1242 SwapForwardMappings::coterminalSwapZedMatrix(
1243 curveState, displacement);
1244 bool logNormal = true;
1245
1246 EvolutionDescription evolution = productClone.evolution();
1247 Size factors = todaysForwards.size();
1248 MarketModelTest::MarketModelType marketModelType = MarketModelTest::ExponentialCorrelationFlatVolatility;
1249 ext::shared_ptr<MarketModel> marketModel =
1250 makeMarketModel(logNormal, evolution, factors, marketModelType);
1251 for (Size i=0; i<todaysForwards.size(); ++i) {
1252 const Matrix& forwardsCovariance = marketModel->totalCovariance(i);
1253 Matrix cotSwapsCovariance =
1254 jacobian * forwardsCovariance * transpose(jacobian);
1255 ext::shared_ptr<PlainVanillaPayoff> payoff(
1256 new PlainVanillaPayoff(Option::Call, todaysForwards[i]+displacement));
1257
1258 Real expectedSwaption =
1259 BlackCalculator(payoff,
1260 todaysCoterminalSwapRates[i]+displacement,
1261 std::sqrt(cotSwapsCovariance[i][i]),
1262 curveState.coterminalSwapAnnuity(0,i) *
1263 todaysDiscounts[0]).value();
1264 subProductExpectedValues.back().values.push_back(expectedSwaption);
1265 }
1266 }
1267
1268
testAllMultiStepProducts()1269 void MarketModelTest::testAllMultiStepProducts() {
1270 std::string testDescription = "all multi-step products ";
1271
1272 using namespace market_model_test;
1273
1274 setup();
1275
1276 MultiProductComposite product;
1277 std::vector<SubProductExpectedValues> subProductExpectedValues;
1278 addForwards(product, subProductExpectedValues);
1279 addOptionLets(product, subProductExpectedValues);
1280 addCoinitialSwaps(product, subProductExpectedValues);
1281 addCoterminalSwapsAndSwaptions(product, subProductExpectedValues);
1282 product.finalize();
1283 testMultiProductComposite(product, subProductExpectedValues,
1284 testDescription);
1285 }
1286
1287
testPeriodAdapter()1288 void MarketModelTest::testPeriodAdapter() {
1289
1290 BOOST_TEST_MESSAGE("Testing period-adaptation routines in LIBOR market model...");
1291
1292 using namespace market_model_test;
1293
1294 setup();
1295 LMMCurveState cs(rateTimes);
1296 cs.setOnForwardRates(todaysForwards);
1297
1298 Size period=2;
1299 Size offset =0;
1300
1301 LMMCurveState bigRateCS(
1302 ForwardForwardMappings::RestrictCurveState(cs,
1303 period,
1304 offset
1305 ));
1306
1307 std::vector<Time> swaptionPaymentTimes(bigRateCS.rateTimes());
1308 swaptionPaymentTimes.pop_back();
1309
1310 std::vector<Time> capletPaymentTimes(swaptionPaymentTimes);
1311
1312
1313 Size numberBigRates = bigRateCS.numberOfRates();
1314
1315 std::vector<ext::shared_ptr<StrikedTypePayoff> > optionletPayoffs(numberBigRates);
1316 std::vector<ext::shared_ptr<StrikedTypePayoff> > swaptionPayoffs(numberBigRates);
1317 std::vector<ext::shared_ptr<StrikedTypePayoff> > displacedOptionletPayoffs(numberBigRates);
1318 std::vector<ext::shared_ptr<StrikedTypePayoff> > displacedSwaptionPayoffs(numberBigRates);
1319
1320 for (Size i=0; i<numberBigRates; ++i)
1321 {
1322 optionletPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
1323 PlainVanillaPayoff(Option::Call, bigRateCS.forwardRate(i)));
1324 swaptionPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
1325 PlainVanillaPayoff(Option::Call, bigRateCS.coterminalSwapRate(i)));
1326 displacedOptionletPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
1327 PlainVanillaPayoff(Option::Call, bigRateCS.forwardRate(i)+displacement));
1328 displacedSwaptionPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
1329 PlainVanillaPayoff(Option::Call, bigRateCS.coterminalSwapRate(i)+displacement));
1330
1331 }
1332
1333 MultiStepPeriodCapletSwaptions theProduct(rateTimes,
1334 capletPaymentTimes,
1335 swaptionPaymentTimes,
1336 optionletPayoffs,
1337 swaptionPayoffs,
1338 period,
1339 offset);
1340
1341 const EvolutionDescription& evolution(theProduct.evolution());
1342
1343 bool logNormal = true;
1344 Size factors = 5;
1345
1346 ext::shared_ptr<MarketModel> originalModel =
1347 makeMarketModel(logNormal,
1348 evolution,
1349 factors,
1350 ExponentialCorrelationAbcdVolatility);
1351
1352 std::vector<Spread> newDisplacements;
1353
1354 ext::shared_ptr<MarketModel> adaptedforwardModel(new FwdPeriodAdapter(originalModel,
1355 period,
1356 offset,
1357 newDisplacements));
1358
1359 ext::shared_ptr<MarketModel> adaptedSwapModel(new
1360 FwdToCotSwapAdapter(adaptedforwardModel));
1361
1362 Matrix finalForwardCovariances(adaptedforwardModel->totalCovariance(adaptedforwardModel->numberOfSteps()-1));
1363 Matrix finalSwapCovariances(adaptedSwapModel->totalCovariance(adaptedSwapModel->numberOfSteps()-1));
1364
1365 std::vector<Volatility> adaptedForwardSds(adaptedforwardModel->numberOfRates());
1366 std::vector<Volatility> adaptedSwapSds(adaptedSwapModel->numberOfRates());
1367 std::vector<Real> approxCapletPrices(adaptedforwardModel->numberOfRates());
1368 std::vector<Real> approxSwaptionPrices(adaptedSwapModel->numberOfRates());
1369
1370 for (Size j=0; j < adaptedSwapModel->numberOfRates(); ++j)
1371 {
1372 adaptedForwardSds[j] = sqrt(finalForwardCovariances[j][j]);
1373 adaptedSwapSds[j] = sqrt(finalSwapCovariances[j][j]);
1374
1375 Real capletAnnuity = todaysDiscounts[0]*bigRateCS.discountRatio(j+1,0)
1376 *bigRateCS.rateTaus()[j];
1377
1378 approxCapletPrices[j] = BlackCalculator(displacedOptionletPayoffs[j],
1379 bigRateCS.forwardRate(j)+displacement,
1380 adaptedForwardSds[j],
1381 capletAnnuity).value();
1382
1383 Real swaptionAnnuity = todaysDiscounts[0]
1384 *bigRateCS.coterminalSwapAnnuity(0,j);
1385
1386 approxSwaptionPrices[j] = BlackCalculator(displacedSwaptionPayoffs[j],
1387 bigRateCS.coterminalSwapRate(j)+displacement,
1388 adaptedSwapSds[j],
1389 swaptionAnnuity).value();
1390 }
1391 SobolBrownianGeneratorFactory generatorFactory(
1392 SobolBrownianGenerator::Diagonal, seed_);
1393
1394
1395
1396 ext::shared_ptr<MarketModelEvolver> evolver = makeMarketModelEvolver(originalModel,
1397 theProduct.suggestedNumeraires(),
1398 generatorFactory,
1399 Pc);
1400
1401 ext::shared_ptr<SequenceStatisticsInc> stats =
1402 simulate(evolver, theProduct);
1403
1404 std::vector<Real> results = stats->mean();
1405 std::vector<Real> errors = stats->errorEstimate();
1406
1407 std::vector<Real> capletErrorsInSds(numberBigRates);
1408 std::vector<Real> swaptionErrorsInSds(numberBigRates);
1409
1410 if (2*numberBigRates != results.size())
1411 BOOST_ERROR("mismatch between the size of the result and the \
1412 number of results");
1413
1414 for (Size i=0; i < numberBigRates; ++i)
1415 {
1416 capletErrorsInSds[i]= (results[i]-approxCapletPrices[i])/errors[i];
1417 swaptionErrorsInSds[i]= (results[i+numberBigRates]-approxSwaptionPrices[i])/errors[i+numberBigRates];
1418 }
1419
1420 Real capletTolerance = 4;
1421 Real swaptionTolerance = 4;
1422
1423
1424 for (Size i=0; i < numberBigRates; ++i) {
1425 if (fabs(capletErrorsInSds[i]) > capletTolerance) {
1426 BOOST_FAIL(io::ordinal(i+1) << "caplet , approx price " <<
1427 approxCapletPrices[i] <<
1428 ", \t simulation price " << results[i] <<
1429 ", \t error in sds " << capletErrorsInSds[i]);
1430 }
1431 }
1432 for (Size i=0; i < numberBigRates; ++i) {
1433 if (fabs(swaptionErrorsInSds[i]) > swaptionTolerance) {
1434 BOOST_FAIL(io::ordinal(i+1) << "swaption, approx price " <<
1435 approxSwaptionPrices[i] <<
1436 ", \t simulation price " << results[i+numberBigRates] <<
1437 ", \t error in sds " << swaptionErrorsInSds[i]);
1438 }
1439 }
1440 }
testCallableSwapNaif()1441 void MarketModelTest::testCallableSwapNaif() {
1442
1443 BOOST_TEST_MESSAGE("Pricing callable swap with naif exercise strategy in a LIBOR market model...");
1444
1445 using namespace market_model_test;
1446
1447 setup();
1448
1449 Real fixedRate = 0.04;
1450
1451 // 0. a payer swap
1452 MultiStepSwap payerSwap(rateTimes, accruals, accruals, paymentTimes,
1453 fixedRate, true);
1454
1455 // 1. the equivalent receiver swap
1456 MultiStepSwap receiverSwap(rateTimes, accruals, accruals, paymentTimes,
1457 fixedRate, false);
1458
1459 //exercise schedule
1460 std::vector<Rate> exerciseTimes(rateTimes);
1461 exerciseTimes.pop_back();
1462 //std::vector<Rate> exerciseTimes;
1463 //for (Size i=2; i<rateTimes.size()-1; i+=2)
1464 // exerciseTimes.push_back(rateTimes[i]);
1465
1466 // naif exercise strategy
1467 std::vector<Rate> swapTriggers(exerciseTimes.size(), fixedRate);
1468 SwapRateTrigger naifStrategy(rateTimes, swapTriggers, exerciseTimes);
1469
1470 // Longstaff-Schwartz exercise strategy
1471 std::vector<std::vector<NodeData> > collectedData;
1472 std::vector<std::vector<Real> > basisCoefficients;
1473 NothingExerciseValue control(rateTimes);
1474 SwapBasisSystem basisSystem(rateTimes,exerciseTimes);
1475 NothingExerciseValue nullRebate(rateTimes);
1476
1477 CallSpecifiedMultiProduct dummyProduct =
1478 CallSpecifiedMultiProduct(receiverSwap, naifStrategy,
1479 ExerciseAdapter(nullRebate));
1480
1481 const EvolutionDescription& evolution = dummyProduct.evolution();
1482
1483 MarketModelType marketModels[] = {
1484 // CalibratedMM,
1485 ExponentialCorrelationFlatVolatility,
1486 ExponentialCorrelationAbcdVolatility };
1487 for (Size j=0; j<LENGTH(marketModels); j++) {
1488
1489 Size testedFactors[] = { 4, // 8,
1490 todaysForwards.size()};
1491 for (Size m=0; m<LENGTH(testedFactors); ++m) {
1492 Size factors = testedFactors[m];
1493
1494 // Composite's ProductSuggested is the Terminal one
1495 MeasureType measures[] = { // ProductSuggested,
1496 MoneyMarketPlus
1497 // MoneyMarket,
1498 // Terminal
1499 };
1500 for (Size k=0; k<LENGTH(measures); k++) {
1501 std::vector<Size> numeraires = makeMeasure(dummyProduct, measures[k]);
1502
1503 bool logNormal = true;
1504 ext::shared_ptr<MarketModel> marketModel =
1505 makeMarketModel(logNormal, evolution, factors, marketModels[j]);
1506
1507
1508 EvolverType evolvers[] = { Pc, Balland, Ipc };
1509 ext::shared_ptr<MarketModelEvolver> evolver;
1510 Size stop =
1511 isInTerminalMeasure(evolution, numeraires) ? 0 : 1;
1512 for (Size i=0; i<LENGTH(evolvers)-stop; i++) {
1513
1514 for (Size n=0; n<1; n++) {
1515 //MTBrownianGeneratorFactory generatorFactory(seed_);
1516 SobolBrownianGeneratorFactory generatorFactory(
1517 SobolBrownianGenerator::Diagonal, seed_);
1518
1519 evolver = makeMarketModelEvolver(marketModel,
1520 numeraires,
1521 generatorFactory,
1522 evolvers[i]);
1523 std::ostringstream config;
1524 config <<
1525 marketModelTypeToString(marketModels[j]) << ", " <<
1526 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
1527 measureTypeToString(measures[k]) << ", " <<
1528 evolverTypeToString(evolvers[i]) << ", " <<
1529 "MT BGF";
1530 if (printReport_)
1531 BOOST_TEST_MESSAGE(" " << config.str());
1532
1533 // use the naif strategy
1534
1535 // 2. bermudan swaption to enter into the payer swap
1536 CallSpecifiedMultiProduct bermudanProduct =
1537 CallSpecifiedMultiProduct(
1538 MultiStepNothing(evolution),
1539 naifStrategy, payerSwap);
1540
1541 // 3. callable receiver swap
1542 CallSpecifiedMultiProduct callableProduct =
1543 CallSpecifiedMultiProduct(
1544 receiverSwap, naifStrategy,
1545 ExerciseAdapter(nullRebate));
1546
1547 // lower bound: evolve all 4 products togheter
1548 MultiProductComposite allProducts;
1549 allProducts.add(payerSwap);
1550 allProducts.add(receiverSwap);
1551 allProducts.add(bermudanProduct);
1552 allProducts.add(callableProduct);
1553 allProducts.finalize();
1554
1555 ext::shared_ptr<SequenceStatisticsInc> stats =
1556 simulate(evolver, allProducts);
1557 checkCallableSwap(*stats, config.str());
1558
1559
1560 // upper bound
1561
1562 //MTBrownianGeneratorFactory uFactory(seed_+142);
1563 SobolBrownianGeneratorFactory uFactory(
1564 SobolBrownianGenerator::Diagonal, seed_+142);
1565 evolver = makeMarketModelEvolver(marketModel,
1566 numeraires,
1567 uFactory,
1568 evolvers[i]);
1569
1570 std::vector<ext::shared_ptr<MarketModelEvolver> >
1571 innerEvolvers;
1572
1573 std::valarray<bool> isExerciseTime =
1574 isInSubset(evolution.evolutionTimes(),
1575 naifStrategy.exerciseTimes());
1576 for (Size s=0; s < isExerciseTime.size(); ++s) {
1577 if (isExerciseTime[s]) {
1578 MTBrownianGeneratorFactory iFactory(seed_+s);
1579 ext::shared_ptr<MarketModelEvolver> e =
1580 makeMarketModelEvolver(marketModel,
1581 numeraires,
1582 iFactory,
1583 evolvers[i],
1584 s);
1585 innerEvolvers.push_back(e);
1586 }
1587 }
1588
1589 Size initialNumeraire = evolver->numeraires().front();
1590 Real initialNumeraireValue =
1591 todaysDiscounts[initialNumeraire];
1592
1593 UpperBoundEngine uEngine(evolver, innerEvolvers,
1594 receiverSwap, nullRebate,
1595 receiverSwap, nullRebate,
1596 naifStrategy,
1597 initialNumeraireValue);
1598 Statistics uStats;
1599 uEngine.multiplePathValues(uStats,255,256);
1600 Real delta = uStats.mean();
1601 Real deltaError = uStats.errorEstimate();
1602 if (printReport_)
1603 BOOST_TEST_MESSAGE(" upper bound delta: " << io::rate(delta) << " +- " << io::rate(deltaError));
1604
1605 }
1606 }
1607 }
1608 }
1609 }
1610 }
1611
testCallableSwapLS()1612 void MarketModelTest::testCallableSwapLS() {
1613
1614 BOOST_TEST_MESSAGE("Pricing callable swap with Longstaff-Schwartz exercise strategy in a LIBOR market model...");
1615
1616 using namespace market_model_test;
1617
1618 setup();
1619
1620 Real fixedRate = 0.04;
1621
1622 // 0. a payer swap
1623 MultiStepSwap payerSwap(rateTimes, accruals, accruals, paymentTimes,
1624 fixedRate, true);
1625
1626 // 1. the equivalent receiver swap
1627 MultiStepSwap receiverSwap(rateTimes, accruals, accruals, paymentTimes,
1628 fixedRate, false);
1629
1630 //exercise schedule
1631 std::vector<Rate> exerciseTimes(rateTimes);
1632 exerciseTimes.pop_back();
1633 //std::vector<Rate> exerciseTimes;
1634 //for (Size i=2; i<rateTimes.size()-1; i+=2)
1635 // exerciseTimes.push_back(rateTimes[i]);
1636
1637 // naif exercise strategy
1638 std::vector<Rate> swapTriggers(exerciseTimes.size(), fixedRate);
1639 SwapRateTrigger naifStrategy(rateTimes, swapTriggers, exerciseTimes);
1640
1641 // Longstaff-Schwartz exercise strategy
1642 std::vector<std::vector<NodeData> > collectedData;
1643 std::vector<std::vector<Real> > basisCoefficients;
1644 NothingExerciseValue control(rateTimes);
1645 SwapBasisSystem basisSystem(rateTimes,exerciseTimes);
1646 NothingExerciseValue nullRebate(rateTimes);
1647
1648 CallSpecifiedMultiProduct dummyProduct =
1649 CallSpecifiedMultiProduct(receiverSwap, naifStrategy,
1650 ExerciseAdapter(nullRebate));
1651
1652 const EvolutionDescription& evolution = dummyProduct.evolution();
1653
1654 MarketModelType marketModels[] = {
1655 // CalibratedMM,
1656 ExponentialCorrelationFlatVolatility,
1657 ExponentialCorrelationAbcdVolatility };
1658 for (Size j=0; j<LENGTH(marketModels); j++) {
1659
1660 Size testedFactors[] = { 4, // 8,
1661 todaysForwards.size()};
1662 for (Size m=0; m<LENGTH(testedFactors); ++m) {
1663 Size factors = testedFactors[m];
1664
1665 // Composite's ProductSuggested is the Terminal one
1666 MeasureType measures[] = { // ProductSuggested,
1667 // MoneyMarketPlus,
1668 MoneyMarket
1669 //Terminal
1670 };
1671 for (Size k=0; k<LENGTH(measures); k++) {
1672 std::vector<Size> numeraires = makeMeasure(dummyProduct, measures[k]);
1673
1674 bool logNormal = true;
1675 ext::shared_ptr<MarketModel> marketModel =
1676 makeMarketModel(logNormal, evolution, factors, marketModels[j]);
1677
1678
1679 EvolverType evolvers[] = { Pc, Balland, Ipc };
1680 ext::shared_ptr<MarketModelEvolver> evolver;
1681 Size stop =
1682 isInTerminalMeasure(evolution, numeraires) ? 0 : 1;
1683 for (Size i=0; i<LENGTH(evolvers)-stop; i++) {
1684
1685 for (Size n=0; n<1; n++) {
1686 //MTBrownianGeneratorFactory generatorFactory(seed_);
1687 SobolBrownianGeneratorFactory generatorFactory(
1688 SobolBrownianGenerator::Diagonal, seed_);
1689
1690 evolver = makeMarketModelEvolver(marketModel,
1691 numeraires,
1692 generatorFactory,
1693 evolvers[i]);
1694 std::ostringstream config;
1695 config <<
1696 marketModelTypeToString(marketModels[j]) << ", " <<
1697 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
1698 measureTypeToString(measures[k]) << ", " <<
1699 evolverTypeToString(evolvers[i]) << ", " <<
1700 "MT BGF";
1701 if (printReport_)
1702 BOOST_TEST_MESSAGE(" " << config.str());
1703
1704 // calculate the exercise strategy
1705 collectNodeData(*evolver,
1706 receiverSwap, basisSystem, nullRebate,
1707 control, trainingPaths_, collectedData);
1708 genericLongstaffSchwartzRegression(collectedData,
1709 basisCoefficients);
1710 LongstaffSchwartzExerciseStrategy exerciseStrategy(
1711 basisSystem, basisCoefficients,
1712 evolution, numeraires,
1713 nullRebate, control);
1714
1715 // 2. bermudan swaption to enter into the payer swap
1716 CallSpecifiedMultiProduct bermudanProduct =
1717 CallSpecifiedMultiProduct(
1718 MultiStepNothing(evolution),
1719 exerciseStrategy, payerSwap);
1720
1721 // 3. callable receiver swap
1722 CallSpecifiedMultiProduct callableProduct =
1723 CallSpecifiedMultiProduct(
1724 receiverSwap, exerciseStrategy,
1725 ExerciseAdapter(nullRebate));
1726
1727 // lower bound: evolve all 4 products togheter
1728 MultiProductComposite allProducts;
1729 allProducts.add(payerSwap);
1730 allProducts.add(receiverSwap);
1731 allProducts.add(bermudanProduct);
1732 allProducts.add(callableProduct);
1733 allProducts.finalize();
1734
1735 ext::shared_ptr<SequenceStatisticsInc> stats =
1736 simulate(evolver, allProducts);
1737 checkCallableSwap(*stats, config.str());
1738
1739
1740 // upper bound
1741
1742 //MTBrownianGeneratorFactory uFactory(seed_+142);
1743 SobolBrownianGeneratorFactory uFactory(
1744 SobolBrownianGenerator::Diagonal, seed_+142);
1745 evolver = makeMarketModelEvolver(marketModel,
1746 numeraires,
1747 uFactory,
1748 evolvers[i]);
1749
1750 std::vector<ext::shared_ptr<MarketModelEvolver> >
1751 innerEvolvers;
1752
1753 std::valarray<bool> isExerciseTime =
1754 isInSubset(evolution.evolutionTimes(),
1755 exerciseStrategy.exerciseTimes());
1756 for (Size s=0; s < isExerciseTime.size(); ++s) {
1757 if (isExerciseTime[s]) {
1758 MTBrownianGeneratorFactory iFactory(seed_+s);
1759 ext::shared_ptr<MarketModelEvolver> e =
1760 makeMarketModelEvolver(marketModel,
1761 numeraires,
1762 iFactory,
1763 evolvers[i],
1764 s);
1765 innerEvolvers.push_back(e);
1766 }
1767 }
1768
1769 Size initialNumeraire = evolver->numeraires().front();
1770 Real initialNumeraireValue =
1771 todaysDiscounts[initialNumeraire];
1772
1773 UpperBoundEngine uEngine(evolver, innerEvolvers,
1774 receiverSwap, nullRebate,
1775 receiverSwap, nullRebate,
1776 exerciseStrategy,
1777 initialNumeraireValue);
1778 Statistics uStats;
1779 uEngine.multiplePathValues(uStats,255,256);
1780 Real delta = uStats.mean();
1781 Real deltaError = uStats.errorEstimate();
1782 if (printReport_)
1783 BOOST_TEST_MESSAGE(" upper bound delta: " << io::rate(delta) << " +- " << io::rate(deltaError));
1784
1785 }
1786 }
1787 }
1788 }
1789 }
1790 }
1791
testCallableSwapAnderson(MarketModelType marketModelType,Size testedFactor)1792 void MarketModelTest::testCallableSwapAnderson(
1793 MarketModelType marketModelType, Size testedFactor) {
1794
1795 using namespace market_model_test;
1796
1797 BOOST_TEST_MESSAGE("Pricing callable swap with Anderson exercise "
1798 "strategy in a LIBOR market model for test factor "
1799 << testedFactor << " and model type "
1800 << marketModelTypeToString(marketModelType)
1801 << "...");
1802
1803 setup();
1804
1805 Real fixedRate = 0.04;
1806
1807 // 0. a payer swap
1808 MultiStepSwap payerSwap(rateTimes, accruals, accruals, paymentTimes,
1809 fixedRate, true);
1810
1811 // 1. the equivalent receiver swap
1812 MultiStepSwap receiverSwap(rateTimes, accruals, accruals, paymentTimes,
1813 fixedRate, false);
1814
1815 //exercise schedule
1816 std::vector<Rate> exerciseTimes(rateTimes);
1817 exerciseTimes.pop_back();
1818 //std::vector<Rate> exerciseTimes;
1819 //for (Size i=2; i<rateTimes.size()-1; i+=2)
1820 // exerciseTimes.push_back(rateTimes[i]);
1821
1822 // naif exercise strategy
1823 std::vector<Rate> swapTriggers(exerciseTimes.size(), fixedRate);
1824 SwapRateTrigger naifStrategy(rateTimes, swapTriggers, exerciseTimes);
1825
1826 // Anderson exercise strategy
1827 std::vector<std::vector<NodeData> > collectedData;
1828 std::vector<std::vector<Real> > parameters;
1829 NothingExerciseValue control(rateTimes);
1830 NothingExerciseValue nullRebate(rateTimes);
1831 TriggeredSwapExercise parametricForm(rateTimes, exerciseTimes,
1832 std::vector<Time>(exerciseTimes.size(),fixedRate));
1833
1834 CallSpecifiedMultiProduct dummyProduct =
1835 CallSpecifiedMultiProduct(receiverSwap, naifStrategy,
1836 ExerciseAdapter(nullRebate));
1837
1838 const EvolutionDescription& evolution = dummyProduct.evolution();
1839
1840 Size factors = testedFactor;
1841
1842 // Composite's ProductSuggested is the Terminal one
1843 MeasureType measures[] = { // ProductSuggested,
1844 // MoneyMarketPlus,
1845 // MoneyMarket,
1846 Terminal
1847 };
1848 for (Size k=0; k<LENGTH(measures); k++) {
1849 std::vector<Size> numeraires = makeMeasure(dummyProduct, measures[k]);
1850 bool logNormal = true;
1851 ext::shared_ptr<MarketModel> marketModel =
1852 makeMarketModel(logNormal, evolution, factors, marketModelType);
1853 EvolverType evolvers[] = { Pc, Balland, Ipc };
1854 ext::shared_ptr<MarketModelEvolver> evolver;
1855 Size stop =
1856 isInTerminalMeasure(evolution, numeraires) ? 0 : 1;
1857 for (Size i=0; i<LENGTH(evolvers)-stop; i++) {
1858 for (Size n=0; n<1; n++) {
1859 //MTBrownianGeneratorFactory generatorFactory(seed_);
1860 SobolBrownianGeneratorFactory generatorFactory(
1861 SobolBrownianGenerator::Diagonal, seed_);
1862 evolver = makeMarketModelEvolver(marketModel,
1863 numeraires,
1864 generatorFactory,
1865 evolvers[i]);
1866 std::ostringstream config;
1867 config <<
1868 marketModelTypeToString(marketModelType) << ", " <<
1869 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
1870 measureTypeToString(measures[k]) << ", " <<
1871 evolverTypeToString(evolvers[i]) << ", " <<
1872 "MT BGF";
1873 if (printReport_)
1874 BOOST_TEST_MESSAGE(" " << config.str());
1875 // 1. calculate the exercise strategy
1876 collectNodeData(*evolver,
1877 receiverSwap, parametricForm, nullRebate,
1878 control, trainingPaths_, collectedData);
1879 Simplex om(0.01);
1880 EndCriteria ec(1000, 100, 1e-8, 1e-16, 1e-8);
1881 Size initialNumeraire = evolver->numeraires().front();
1882 Real initialNumeraireValue = todaysDiscounts[initialNumeraire];
1883 Real firstPassValue = genericEarlyExerciseOptimization(
1884 collectedData, parametricForm, parameters, ec, om) *
1885 initialNumeraireValue;
1886 if (printReport_)
1887 BOOST_TEST_MESSAGE(" initial estimate: " << io::rate(firstPassValue));
1888 ParametricExerciseAdapter exerciseStrategy(parametricForm, parameters);
1889
1890 // 2. bermudan swaption to enter into the payer swap
1891 CallSpecifiedMultiProduct bermudanProduct =
1892 CallSpecifiedMultiProduct(
1893 MultiStepNothing(evolution),
1894 exerciseStrategy, payerSwap);
1895
1896 // 3. callable receiver swap
1897 CallSpecifiedMultiProduct callableProduct =
1898 CallSpecifiedMultiProduct(
1899 receiverSwap, exerciseStrategy,
1900 ExerciseAdapter(nullRebate));
1901 // lower bound: evolve all 4 products togheter
1902 MultiProductComposite allProducts;
1903 allProducts.add(payerSwap);
1904 allProducts.add(receiverSwap);
1905 allProducts.add(bermudanProduct);
1906 allProducts.add(callableProduct);
1907 allProducts.finalize();
1908 ext::shared_ptr<SequenceStatisticsInc> stats =
1909 simulate(evolver, allProducts);
1910 checkCallableSwap(*stats, config.str());
1911
1912 // upper bound
1913 //MTBrownianGeneratorFactory uFactory(seed_+142);
1914 SobolBrownianGeneratorFactory uFactory(
1915 SobolBrownianGenerator::Diagonal, seed_+142);
1916 evolver = makeMarketModelEvolver(marketModel,
1917 numeraires,
1918 uFactory,
1919 evolvers[i]);
1920 std::vector<ext::shared_ptr<MarketModelEvolver> >
1921 innerEvolvers;
1922 std::valarray<bool> isExerciseTime =
1923 isInSubset(evolution.evolutionTimes(),
1924 exerciseStrategy.exerciseTimes());
1925 for (Size s=0; s < isExerciseTime.size(); ++s) {
1926 if (isExerciseTime[s]) {
1927 MTBrownianGeneratorFactory iFactory(seed_+s);
1928 ext::shared_ptr<MarketModelEvolver> e =
1929 makeMarketModelEvolver(marketModel,
1930 numeraires,
1931 iFactory,
1932 evolvers[i],
1933 s);
1934 innerEvolvers.push_back(e);
1935 }
1936 }
1937 UpperBoundEngine uEngine(evolver, innerEvolvers,
1938 receiverSwap, nullRebate,
1939 receiverSwap, nullRebate,
1940 exerciseStrategy,
1941 initialNumeraireValue);
1942 Statistics uStats;
1943 uEngine.multiplePathValues(uStats,255,256);
1944 Real delta = uStats.mean();
1945 Real deltaError = uStats.errorEstimate();
1946 if (printReport_)
1947 BOOST_TEST_MESSAGE(" upper bound delta: " << io::rate(delta) << " +- " << io::rate(deltaError));
1948
1949 }
1950 }
1951 }
1952 }
1953
1954
1955
testGreeks()1956 void MarketModelTest::testGreeks() {
1957
1958 BOOST_TEST_MESSAGE("Testing caplet greeks in a lognormal forward rate market model using partial proxy simulation...");
1959
1960 using namespace market_model_test;
1961
1962 setup();
1963
1964 std::vector<ext::shared_ptr<Payoff> > payoffs(todaysForwards.size());
1965 std::vector<ext::shared_ptr<StrikedTypePayoff> >
1966 displacedPayoffs(todaysForwards.size());
1967 for (Size i=0; i<todaysForwards.size(); ++i) {
1968 payoffs[i] = ext::shared_ptr<Payoff>(new
1969 //PlainVanillaPayoff(Option::Call, todaysForwards[i]));
1970 CashOrNothingPayoff(Option::Call, todaysForwards[i], 0.01));
1971 displacedPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
1972 //PlainVanillaPayoff(Option::Call, todaysForwards[i]+displacement));
1973 CashOrNothingPayoff(Option::Call, todaysForwards[i]+displacement, 0.01));
1974 }
1975
1976 MultiStepOptionlets product(rateTimes, accruals,
1977 paymentTimes, payoffs);
1978
1979 const EvolutionDescription& evolution = product.evolution();
1980
1981 MarketModelType marketModels[] = {
1982 // CalibratedMM,
1983 // ExponentialCorrelationFlatVolatility,
1984 ExponentialCorrelationAbcdVolatility };
1985 for (Size j=0; j<LENGTH(marketModels); j++) {
1986
1987 Size testedFactors[] = { 4, 8, todaysForwards.size() };
1988 for (Size m=0; m<LENGTH(testedFactors); ++m) {
1989 Size factors = testedFactors[m];
1990
1991 MeasureType measures[] = { //MoneyMarketPlus,
1992 MoneyMarket//,
1993 //Terminal
1994 };
1995 for (Size k=0; k<LENGTH(measures); k++) {
1996 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
1997
1998 for (Size n=0; n<1; n++) {
1999 //MTBrownianGeneratorFactory generatorFactory(seed_);
2000 SobolBrownianGeneratorFactory generatorFactory(
2001 SobolBrownianGenerator::Diagonal,
2002 seed_);
2003
2004 bool logNormal = true;
2005 ext::shared_ptr<MarketModel> marketModel =
2006 makeMarketModel(logNormal, evolution, factors,
2007 marketModels[j]);
2008
2009 ext::shared_ptr<MarketModelEvolver> evolver(new
2010 LogNormalFwdRateEuler(marketModel,
2011 generatorFactory,
2012 numeraires));
2013 SequenceStatisticsInc stats(product.numberOfProducts());
2014
2015
2016 std::vector<Size> startIndexOfConstraint;
2017 std::vector<Size> endIndexOfConstraint;
2018
2019 for (Size i=0; i<evolution.evolutionTimes().size(); ++i) {
2020 startIndexOfConstraint.push_back(i);
2021 endIndexOfConstraint.push_back(i+1);
2022 }
2023
2024
2025 std::vector<
2026 std::vector<ext::shared_ptr<ConstrainedEvolver> > >
2027 constrainedEvolvers;
2028 std::vector<std::vector<std::vector<Real> > > diffWeights;
2029 std::vector<std::vector<SequenceStatisticsInc> > greekStats;
2030
2031 std::vector<ext::shared_ptr<ConstrainedEvolver> >
2032 deltaGammaEvolvers;
2033 std::vector<std::vector<Real> > deltaGammaWeights(
2034 2, std::vector<Real>(3));
2035 std::vector<SequenceStatisticsInc> deltaGammaStats(2,stats);
2036
2037
2038 Spread forwardBump = 1.0e-6;
2039 marketModel =
2040 makeMarketModel(logNormal, evolution, factors,
2041 marketModels[j], -forwardBump);
2042 deltaGammaEvolvers.push_back(
2043 ext::shared_ptr<ConstrainedEvolver>(new
2044 LogNormalFwdRateEulerConstrained(marketModel,
2045 generatorFactory,
2046 numeraires)));
2047 deltaGammaEvolvers.back()->setConstraintType(
2048 startIndexOfConstraint, endIndexOfConstraint);
2049 marketModel =
2050 makeMarketModel(logNormal, evolution, factors,
2051 marketModels[j], forwardBump);
2052 deltaGammaEvolvers.push_back(
2053 ext::shared_ptr<ConstrainedEvolver>(new
2054 LogNormalFwdRateEulerConstrained(marketModel,
2055 generatorFactory,
2056 numeraires)));
2057 deltaGammaEvolvers.back()->setConstraintType(
2058 startIndexOfConstraint, endIndexOfConstraint);
2059
2060 deltaGammaWeights[0][0] = 0.0;
2061 deltaGammaWeights[0][1] = -1.0/(2.0*forwardBump);
2062 deltaGammaWeights[0][2] = 1.0/(2.0*forwardBump);
2063
2064 deltaGammaWeights[1][0] = -2.0/(forwardBump*forwardBump);
2065 deltaGammaWeights[1][1] = 1.0/(forwardBump*forwardBump);
2066 deltaGammaWeights[1][2] = 1.0/(forwardBump*forwardBump);
2067
2068
2069 std::vector<ext::shared_ptr<ConstrainedEvolver> >
2070 vegaEvolvers;
2071 std::vector<std::vector<Real> > vegaWeights(
2072 1, std::vector<Real>(3));
2073 std::vector<SequenceStatisticsInc> vegaStats(1,stats);
2074
2075 Volatility volBump = 1.0e-4;
2076 marketModel =
2077 makeMarketModel(logNormal, evolution, factors,
2078 marketModels[j], 0.0, -volBump);
2079 vegaEvolvers.push_back(
2080 ext::shared_ptr<ConstrainedEvolver>(new
2081 LogNormalFwdRateEulerConstrained(marketModel,
2082 generatorFactory,
2083 numeraires)));
2084 vegaEvolvers.back()->setConstraintType(
2085 startIndexOfConstraint, endIndexOfConstraint);
2086 marketModel =
2087 makeMarketModel(logNormal, evolution, factors,
2088 marketModels[j], 0.0, volBump);
2089 vegaEvolvers.push_back(
2090 ext::shared_ptr<ConstrainedEvolver>(new
2091 LogNormalFwdRateEulerConstrained(marketModel,
2092 generatorFactory,
2093 numeraires)));
2094 vegaEvolvers.back()->setConstraintType(
2095 startIndexOfConstraint, endIndexOfConstraint);
2096
2097 vegaWeights[0][0] = 0.0;
2098 vegaWeights[0][1] = -1.0/(2.0*volBump);
2099 vegaWeights[0][2] = 1.0/(2.0*volBump);
2100
2101
2102
2103 constrainedEvolvers.push_back(deltaGammaEvolvers);
2104 diffWeights.push_back(deltaGammaWeights);
2105 greekStats.push_back(deltaGammaStats);
2106
2107 constrainedEvolvers.push_back(vegaEvolvers);
2108 diffWeights.push_back(vegaWeights);
2109 greekStats.push_back(vegaStats);
2110
2111 std::ostringstream config;
2112 config <<
2113 marketModelTypeToString(marketModels[j]) << ", " <<
2114 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
2115 measureTypeToString(measures[k]) << ", " <<
2116 "MT BGF";
2117 if (printReport_)
2118 BOOST_TEST_MESSAGE(" " << config.str());
2119
2120 Size initialNumeraire = evolver->numeraires().front();
2121 Real initialNumeraireValue =
2122 todaysDiscounts[initialNumeraire];
2123
2124 ProxyGreekEngine engine(evolver,
2125 constrainedEvolvers, diffWeights,
2126 startIndexOfConstraint,
2127 endIndexOfConstraint,
2128 product,
2129 initialNumeraireValue);
2130
2131 engine.multiplePathValues(stats, greekStats, paths_);
2132
2133 std::vector<Real> values = stats.mean();
2134 std::vector<Real> errors = stats.errorEstimate();
2135 std::vector<Real> deltas = greekStats[0][0].mean();
2136 std::vector<Real> deltaErrors = greekStats[0][0].errorEstimate();
2137 std::vector<Real> gammas = greekStats[0][1].mean();
2138 std::vector<Real> gammaErrors = greekStats[0][1].errorEstimate();
2139 std::vector<Real> vegas = greekStats[1][0].mean();
2140 std::vector<Real> vegaErrors = greekStats[1][0].errorEstimate();
2141
2142 std::vector<DiscountFactor> discPlus(todaysForwards.size()+1, todaysDiscounts[0]);
2143 std::vector<DiscountFactor> discMinus(todaysForwards.size()+1, todaysDiscounts[0]);
2144 std::vector<Rate> fwdPlus(todaysForwards.size());
2145 std::vector<Rate> fwdMinus(todaysForwards.size());
2146 std::vector<Rate> pricePlus(todaysForwards.size());
2147 std::vector<Rate> price0(todaysForwards.size());
2148 std::vector<Rate> priceMinus(todaysForwards.size());
2149 for (Size i=0; i<todaysForwards.size(); ++i) {
2150 Time tau = rateTimes[i+1]-rateTimes[i];
2151 fwdPlus[i]=todaysForwards[i]+forwardBump;
2152 fwdMinus[i]=todaysForwards[i]-forwardBump;
2153 discPlus[i+1]=discPlus[i]/(1.0+fwdPlus[i]*tau);
2154 discMinus[i+1]=discMinus[i]/(1.0+fwdMinus[i]*tau);
2155 pricePlus[i]=BlackCalculator(displacedPayoffs[i], fwdPlus[i],
2156 volatilities[i]*sqrt(rateTimes[i]),
2157 discPlus[i+1]*tau).value();
2158 price0[i]=BlackCalculator(displacedPayoffs[i], todaysForwards[i],
2159 volatilities[i]*sqrt(rateTimes[i]),
2160 todaysDiscounts[i+1]*tau).value();
2161 priceMinus[i]=BlackCalculator(displacedPayoffs[i], fwdMinus[i],
2162 volatilities[i]*sqrt(rateTimes[i]),
2163 discMinus[i+1]*tau).value();
2164 }
2165
2166 for (Size i=0; i<product.numberOfProducts(); ++i) {
2167 Real numDelta = (pricePlus[i]-priceMinus[i])/(2.0*forwardBump);
2168 Real numGamma = (pricePlus[i]-2*price0[i]+priceMinus[i])/(forwardBump*forwardBump);
2169 if (printReport_) {
2170 BOOST_TEST_MESSAGE(io::ordinal(i+1) << " caplet: "
2171 << "value = " << price0[i] << ", "
2172 << "delta = " << numDelta << ", "
2173 << "gamma = " << numGamma);
2174 BOOST_TEST_MESSAGE(io::ordinal(i+1) << " caplet: "
2175 << "value = " << values[i]
2176 << " +- " << errors[i]
2177 << " (" << (values[i]-price0[i])/errors[i] << " s.e.), "
2178 << "delta = " << deltas[i]
2179 << " +- " << deltaErrors[i]
2180 << " (" << (deltas[i]-numDelta)/deltaErrors[i] << " s.e.), "
2181 << "gamma = " << gammas[i]
2182 << " +- " << gammaErrors[i]
2183 << " (" << (gammas[i]-numGamma)/gammaErrors[i] << " s.e.), "
2184 << "vega = " << vegas[i]
2185 << " +- " << vegaErrors[i]);
2186 }
2187 }
2188 }
2189 }
2190 }
2191 }
2192 }
2193
2194 // pathwise deltas
2195
2196
testPathwiseGreeks()2197 void MarketModelTest::testPathwiseGreeks()
2198 {
2199
2200 BOOST_TEST_MESSAGE("Testing caplet deltas in a lognormal forward rate market model using pathwise method...");
2201
2202 using namespace market_model_test;
2203
2204 setup();
2205
2206
2207
2208 std::vector<ext::shared_ptr<Payoff> > payoffs(todaysForwards.size());
2209 std::vector<ext::shared_ptr<StrikedTypePayoff> >
2210 displacedPayoffs(todaysForwards.size());
2211 for (Size i=0; i<todaysForwards.size(); ++i) {
2212 payoffs[i] = ext::shared_ptr<Payoff>(new
2213 PlainVanillaPayoff(Option::Call, todaysForwards[i]));
2214 //CashOrNothingPayoff(Option::Call, todaysForwards[i], 0.01));
2215 displacedPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
2216 PlainVanillaPayoff(Option::Call, todaysForwards[i]+displacement));
2217 //CashOrNothingPayoff(Option::Call, todaysForwards[i]+displacement, 0.01));
2218 }
2219
2220 for (Size whichProduct=0; whichProduct<2; ++whichProduct)
2221 {
2222 MarketModelPathwiseMultiDeflatedCaplet product1(rateTimes, accruals,
2223 paymentTimes, todaysForwards);
2224
2225 MarketModelPathwiseMultiCaplet product2(rateTimes, accruals,
2226 paymentTimes, todaysForwards);
2227
2228 Clone<MarketModelPathwiseMultiProduct> product;
2229
2230 if (whichProduct == 0)
2231 product = product2;
2232 else
2233 product = product1;
2234
2235
2236 MultiStepOptionlets productDummy(rateTimes, accruals,
2237 paymentTimes, payoffs);
2238
2239
2240
2241 EvolutionDescription evolution = product->evolution();
2242
2243 MarketModelType marketModels[] = {
2244 // CalibratedMM,
2245 // ExponentialCorrelationFlatVolatility,
2246 ExponentialCorrelationAbcdVolatility };
2247
2248 for (Size j=0; j<LENGTH(marketModels); j++)
2249 {
2250
2251 Size testedFactors[] = { 2
2252 //, 4, 8, todaysForwards.size()
2253 };
2254
2255 for (Size m=0; m<LENGTH(testedFactors); ++m)
2256 {
2257 Size factors = testedFactors[m];
2258
2259 MeasureType measures[] = {
2260 MoneyMarket
2261 };
2262
2263 for (Size k=0; k<LENGTH(measures); k++)
2264 {
2265 std::vector<Size> numeraires = makeMeasure(productDummy, measures[k]);
2266
2267 for (Size n=0; n<1; n++)
2268 {
2269 MTBrownianGeneratorFactory generatorFactory(seed_);
2270
2271 bool logNormal = true;
2272 ext::shared_ptr<MarketModel> marketModel =
2273 makeMarketModel(logNormal, evolution, factors,
2274 marketModels[j]);
2275
2276 LogNormalFwdRateEuler evolver(marketModel,
2277 generatorFactory,
2278 numeraires);
2279 SequenceStatisticsInc stats(product->numberOfProducts()*(todaysForwards.size()+1));
2280
2281
2282
2283
2284
2285
2286
2287 Spread forwardBump = 1.0e-6;
2288
2289 std::ostringstream config;
2290 config <<
2291 marketModelTypeToString(marketModels[j]) << ", " <<
2292 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
2293 measureTypeToString(measures[k]) << ", " <<
2294 "MT BGF";
2295 if (printReport_)
2296 BOOST_TEST_MESSAGE(" " << config.str());
2297
2298 Size initialNumeraire = evolver.numeraires().front();
2299 Real initialNumeraireValue =
2300 todaysDiscounts[initialNumeraire];
2301
2302
2303
2304 {
2305
2306 PathwiseAccountingEngine accountingengine(ext::make_shared<LogNormalFwdRateEuler>(evolver), // method relies heavily on LMM Euler
2307 *product,
2308 marketModel, // we need pseudo-roots and displacements
2309 initialNumeraireValue);
2310
2311
2312
2313 accountingengine.multiplePathValues(stats,paths_);
2314 }
2315
2316
2317 std::vector<Real> valuesAndDeltas = stats.mean();
2318 std::vector<Real> errors = stats.errorEstimate();
2319
2320 std::vector<Real> prices(product->numberOfProducts());
2321 std::vector<Real> priceErrors(product->numberOfProducts());
2322
2323 Matrix deltas( product->numberOfProducts(), todaysForwards.size());
2324 Matrix deltasErrors( product->numberOfProducts(), todaysForwards.size());
2325 std::vector<Real> modelPrices(product->numberOfProducts());
2326
2327
2328 for (Size i=0; i < product->numberOfProducts(); ++i)
2329 {
2330 prices[i] = valuesAndDeltas[i];
2331
2332 priceErrors[i] = errors[i];
2333
2334 modelPrices[i] = BlackCalculator(displacedPayoffs[i], todaysForwards[i],
2335 volatilities[i]*sqrt(rateTimes[i]),
2336 todaysDiscounts[i+1]*(rateTimes[i+1]-rateTimes[i])).value();
2337
2338
2339 for (Size j=0; j < todaysForwards.size(); ++j)
2340 {
2341 deltas[i][j] = valuesAndDeltas[(i+1)*product->numberOfProducts()+j];
2342 deltasErrors[i][j] = errors[(i+1)* product->numberOfProducts()+j];
2343
2344 }
2345
2346 }
2347
2348 Matrix modelDeltas(product->numberOfProducts(), todaysForwards.size());
2349
2350
2351 std::vector<DiscountFactor> discPlus(todaysForwards.size()+1, todaysDiscounts[0]);
2352 std::vector<DiscountFactor> discMinus(todaysForwards.size()+1, todaysDiscounts[0]);
2353 std::vector<Rate> fwdPlus(todaysForwards.size());
2354 std::vector<Rate> fwdMinus(todaysForwards.size());
2355
2356
2357 for (Size i=0; i < todaysForwards.size(); ++i)
2358 {
2359 for (Size j=0; j < todaysForwards.size(); ++j)
2360 {
2361 if (i != j)
2362 {
2363 fwdPlus[j] = todaysForwards[j];
2364 fwdMinus[j] = todaysForwards[j];
2365
2366 }
2367 else
2368 {
2369 fwdPlus[j] = todaysForwards[j]+forwardBump;
2370 fwdMinus[j] = todaysForwards[j]-forwardBump;
2371 }
2372
2373 Time tau = rateTimes[j+1]-rateTimes[j];
2374 discPlus[j+1]=discPlus[j]/(1.0+fwdPlus[j]*tau);
2375 discMinus[j+1]=discMinus[j]/(1.0+fwdMinus[j]*tau);
2376
2377 }
2378
2379 for (Size k=0; k < product->numberOfProducts(); ++k)
2380 {
2381 Real tau = rateTimes[k+1] - rateTimes[k];
2382 Real priceUp = BlackCalculator(displacedPayoffs[k], fwdPlus[k],
2383 volatilities[k]*sqrt(rateTimes[k]),
2384 discPlus[k+1]*tau).value();
2385 Real priceDown = BlackCalculator(displacedPayoffs[k], fwdMinus[k],
2386 volatilities[k]*sqrt(rateTimes[k]),
2387 discMinus[k+1]*tau).value();
2388
2389 modelDeltas[k][i] = (priceUp-priceDown)/(2*forwardBump);
2390
2391 }
2392 }
2393
2394
2395 Integer numberErrors =0;
2396
2397 for (Size i=0; i<product->numberOfProducts(); ++i)
2398 {
2399
2400 Real thisPrice = prices[i];
2401 Real thisModelPrice = modelPrices[i];
2402 Real priceErrorInSds = ((thisPrice - thisModelPrice)/priceErrors[i]);
2403
2404 Real errorTheshold = 3.5;
2405
2406 if (fabs(priceErrorInSds) > errorTheshold)
2407 {
2408 BOOST_TEST_MESSAGE("Caplet " << i << " price " << prices[i] << " model price " << modelPrices[i]
2409 << " Standard error: " <<priceErrors[i] << " errors in sds: " << priceErrorInSds);
2410
2411 ++numberErrors;
2412
2413 }
2414
2415 Real threshold =1e-10;
2416
2417 for (Size j =0; j < todaysForwards.size(); ++j)
2418 {
2419 Real delta = deltas[i][j];
2420 Real modelDelta = modelDeltas[i][j];
2421
2422 Real deltaErrorInSds =100;
2423
2424 if (deltasErrors[i][j] > 0.0)
2425 deltaErrorInSds = (( delta - modelDelta )/deltasErrors[i][j]);
2426 else
2427 if (fabs(modelDelta -delta) < threshold) // to cope with zero over zero
2428 deltaErrorInSds =0.0;
2429
2430 if (fabs(deltaErrorInSds) > errorTheshold)
2431 {
2432
2433 BOOST_TEST_MESSAGE("Caplet " << i << " delta " << j << "has value " << deltas[i][j] << " model value " << modelDeltas[i][j]
2434 << " Standard error: " <<deltasErrors[i][j] << " errors in sds: " << deltaErrorInSds);
2435
2436 ++numberErrors;
2437 }
2438
2439
2440 }
2441
2442
2443 }
2444
2445 if (numberErrors >0)
2446 BOOST_FAIL("Pathwise greeks test has " << numberErrors <<"\n");
2447 }
2448 }
2449 }
2450 }
2451 }
2452 }
2453
testPathwiseVegas()2454 void MarketModelTest::testPathwiseVegas()
2455 {
2456
2457 BOOST_TEST_MESSAGE(
2458 "Testing pathwise vegas in a lognormal forward rate market model...");
2459
2460 using namespace market_model_test;
2461
2462 setup();
2463
2464
2465 std::vector<ext::shared_ptr<Payoff> > payoffs(todaysForwards.size());
2466 std::vector<ext::shared_ptr<StrikedTypePayoff> >
2467 displacedPayoffs(todaysForwards.size());
2468 for (Size i=0; i<todaysForwards.size(); ++i) {
2469 payoffs[i] = ext::shared_ptr<Payoff>(new
2470 PlainVanillaPayoff(Option::Call, todaysForwards[i]));
2471 //CashOrNothingPayoff(Option::Call, todaysForwards[i], 0.01));
2472 displacedPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
2473 PlainVanillaPayoff(Option::Call, todaysForwards[i]+displacement));
2474 //CashOrNothingPayoff(Option::Call, todaysForwards[i]+displacement, 0.01));
2475 }
2476
2477
2478 MultiStepOptionlets product(rateTimes, accruals,
2479 paymentTimes, payoffs);
2480
2481 MarketModelPathwiseMultiCaplet caplets(rateTimes, accruals,
2482 paymentTimes, todaysForwards);
2483
2484
2485 MarketModelPathwiseMultiDeflatedCaplet capletsDeflated(rateTimes, accruals,
2486 paymentTimes, todaysForwards);
2487
2488 LMMCurveState cs(rateTimes);
2489 cs.setOnForwardRates(todaysForwards);
2490
2491
2492
2493
2494 EvolutionDescription evolution = product.evolution();
2495 Size steps = evolution.numberOfSteps();
2496 Size numberRates = evolution.numberOfRates();
2497
2498 Real bumpSizeNumericalDifferentiation = 1E-6;
2499 Real vegaBumpSize = 1e-2;
2500 Size pathsToDo =10; // for the numerical differentiation test we are requiring equality on each path so this is actually quite strict
2501 Size pathsToDoSimulation = paths_;
2502 Size bumpIncrement = 1 + evolution.numberOfSteps()/3;
2503 Real numericalBumpSizeForSwaptionPseudo =1E-7;
2504
2505 Real multiplier = 50; // how many times the bump size squared, the numerical differentation is allowed to differ by
2506 // printReport_ = true;
2507 Real maxError =0.0;
2508 Size numberSwaptionPseudoFailures =0;
2509 Size numberCapPseudoFailures = 0;
2510 Size numberCapImpVolFailures = 0;
2511 Size numberCapVolPseudoFailures =0;
2512 Real swaptionPseudoTolerance = 1e-8;
2513 Real impVolTolerance = 1e-5;
2514 Real capStrike = meanForward;
2515 Real initialNumeraireValue =0.95;
2516
2517
2518 MarketModelType marketModels[] =
2519 {
2520 // CalibratedMM,
2521 // ExponentialCorrelationFlatVolatility,
2522 ExponentialCorrelationAbcdVolatility
2523 };
2524 /////////////////////////////////// test derivative of swaption implied vol with respect to pseudo-root elements
2525
2526 for (Size j=0; j<LENGTH(marketModels); j++)
2527 {
2528
2529 Size testedFactors[] = { std::min<Size>(3UL,todaysForwards.size())
2530 // todaysForwards.size()
2531 //, 4, 8,
2532 };
2533
2534
2535
2536
2537 for (Size m=0; m<LENGTH(testedFactors); ++m)
2538 {
2539 Size factors = testedFactors[m];
2540
2541
2542 bool logNormal = true;
2543
2544 ext::shared_ptr<MarketModel> marketModel =
2545 makeMarketModel(logNormal, evolution, factors,
2546 marketModels[j]);
2547
2548 Size startIndex = std::min<Size>(1,evolution.numberOfRates()-2) ;
2549 Size endIndex = evolution.numberOfRates()-1;
2550
2551 SwaptionPseudoDerivative derivative(marketModel,
2552 startIndex,
2553 endIndex);
2554
2555 std::vector<Matrix> pseudoRoots;
2556 for (Size k=0; k < marketModel->numberOfSteps(); ++k)
2557 pseudoRoots.push_back( marketModel->pseudoRoot(k));
2558
2559 // test that the derivative of swaption implied vols to the pseudo-root elements are correct, finite differencing versus analytic value
2560
2561 for (Size step=0; step < evolution.numberOfSteps(); ++ step)
2562 {
2563 for (Size l=0; l < evolution.numberOfRates(); ++l)
2564 for (Size f=0; f < factors; ++f)
2565 {
2566
2567 // change one pseudo root element in the calibration by adding a bump to it
2568
2569 pseudoRoots[step][l][f] += numericalBumpSizeForSwaptionPseudo;
2570
2571
2572 // create new market model with the pseudo root bumped
2573
2574 PseudoRootFacade bumpedUp(pseudoRoots,rateTimes,marketModel->initialRates(),marketModel->displacements());
2575
2576
2577 // compute the implied vol of the swaption with the bumped pseudo roots
2578
2579 Real upImpVol = SwapForwardMappings::swaptionImpliedVolatility(bumpedUp,
2580 startIndex,
2581 endIndex);
2582
2583
2584 // undo the bump
2585
2586 pseudoRoots[step][l][f] -= numericalBumpSizeForSwaptionPseudo;
2587
2588
2589 // bump down
2590
2591 pseudoRoots[step][l][f] -= numericalBumpSizeForSwaptionPseudo;
2592
2593
2594 // create facade for the bumped down pseudo roots
2595
2596 PseudoRootFacade bumpedDown(pseudoRoots,rateTimes,marketModel->initialRates(),marketModel->displacements());
2597
2598 // compute the implied vol of the swaption with the bumped down pseudo roots
2599
2600 Real downImpVol = SwapForwardMappings::swaptionImpliedVolatility(bumpedDown,
2601 startIndex,
2602 endIndex);
2603
2604 // undo bumping
2605
2606 pseudoRoots[step][l][f] += numericalBumpSizeForSwaptionPseudo;
2607
2608 // use symmetric finite differencing to compute the change in the swaptions implied vol for changes in this pseudo-root element
2609
2610 Real volDeriv = (upImpVol-downImpVol)/(2.0*numericalBumpSizeForSwaptionPseudo);
2611
2612 Real modelVal = derivative.volatilityDerivative(step)[l][f];
2613
2614 Real error = volDeriv - modelVal;
2615
2616 if (fabs(error) > swaptionPseudoTolerance)
2617 ++numberSwaptionPseudoFailures;
2618
2619
2620
2621 }
2622
2623 }
2624
2625 if (numberSwaptionPseudoFailures >0)
2626 BOOST_ERROR("swaption pseudo test failed " << numberSwaptionPseudoFailures << " times" );
2627 }
2628 }
2629
2630 /////////////////////////////////////
2631
2632 for (Size j=0; j<LENGTH(marketModels); j++)
2633 {
2634
2635 Size testedFactors[] = { std::min<Size>(3UL,todaysForwards.size())
2636 // todaysForwards.size()
2637 //, 4, 8,
2638 };
2639
2640
2641
2642
2643 for (Size m=0; m<LENGTH(testedFactors); ++m)
2644 {
2645 Size factors = testedFactors[m];
2646
2647
2648 bool logNormal = true;
2649
2650 ext::shared_ptr<MarketModel> marketModel =
2651 makeMarketModel(logNormal, evolution, factors,
2652 marketModels[j]);
2653
2654 for (Size startIndex = 1; startIndex < evolution.numberOfRates()-1; ++startIndex)
2655 for (Size endIndex = startIndex+1; endIndex < evolution.numberOfRates(); ++endIndex)
2656 {
2657
2658 CapPseudoDerivative derivative(marketModel,
2659 capStrike,
2660 startIndex,
2661 endIndex, initialNumeraireValue);
2662
2663 std::vector<Matrix> pseudoRoots;
2664 for (Size k=0; k < marketModel->numberOfSteps(); ++k)
2665 pseudoRoots.push_back( marketModel->pseudoRoot(k));
2666
2667 // test cap price derivatives with respect to pseudo-root elements
2668
2669 for (Size step=0; step < evolution.numberOfSteps(); ++ step)
2670 {
2671 for (Size l=0; l < evolution.numberOfRates(); ++l)
2672 for (Size f=0; f < factors; ++f)
2673 {
2674
2675 // similar to swaption pseudo derivative test but with prices not implied vols
2676
2677 pseudoRoots[step][l][f] += numericalBumpSizeForSwaptionPseudo;
2678
2679 PseudoRootFacade bumpedUp(pseudoRoots,rateTimes,marketModel->initialRates(),marketModel->displacements());
2680
2681 // get total covariances of rates with bumped up pseudo-roots , we really only need the variances
2682 Matrix totalCovUp(bumpedUp.totalCovariance( marketModel->numberOfSteps()-1));
2683
2684
2685 pseudoRoots[step][l][f] -= numericalBumpSizeForSwaptionPseudo;
2686
2687 pseudoRoots[step][l][f] -= numericalBumpSizeForSwaptionPseudo;
2688
2689 PseudoRootFacade bumpedDown(pseudoRoots,rateTimes,marketModel->initialRates(),marketModel->displacements());
2690
2691 // get total covariances of rates with bumped down pseudo-roots , we really only need the variances
2692 Matrix totalCovDown(bumpedDown.totalCovariance( marketModel->numberOfSteps()-1));
2693
2694
2695 pseudoRoots[step][l][f] += numericalBumpSizeForSwaptionPseudo;
2696
2697
2698 // we have to loop through all the caplets underlying the cap to get the price
2699
2700 Real priceDeriv=0.0;
2701 for (Size k=startIndex; k < endIndex; ++k)
2702 {
2703 Real upSd = sqrt(totalCovUp[k][k]);
2704 Real downSd = sqrt(totalCovDown[k][k]);
2705
2706 Real annuity = todaysDiscounts[k+1]* marketModel->evolution().rateTaus()[k];
2707 Real forward = todaysForwards[k];
2708
2709
2710 Real upPrice = blackFormula(Option::Call,
2711 capStrike,
2712 forward,
2713 upSd,
2714 annuity,
2715 marketModel->displacements()[k]);
2716
2717
2718 Real downPrice = blackFormula(Option::Call,
2719 capStrike,
2720 forward,
2721 downSd,
2722 annuity,
2723 marketModel->displacements()[k]);
2724
2725
2726 priceDeriv += (upPrice-downPrice)/(2.0*numericalBumpSizeForSwaptionPseudo);
2727
2728 }
2729
2730 Real modelVal = derivative.priceDerivative(step)[l][f];
2731
2732 Real error = priceDeriv - modelVal;
2733
2734 if (fabs(error) > swaptionPseudoTolerance)
2735 ++numberCapPseudoFailures;
2736
2737
2738
2739 }
2740
2741 }
2742
2743 // test the implied vol of the cap, each underlying caplet has a different implied vol and the cap's is different again
2744
2745 Real impVol = derivative.impliedVolatility();
2746
2747 Matrix totalCov(marketModel->totalCovariance(evolution.numberOfSteps()-1 ) );
2748 Real priceConstVol =0.0;
2749 Real priceVarVol =0.0;
2750
2751 for (Size m= startIndex; m < endIndex; ++m)
2752 {
2753 Real annuity = todaysDiscounts[m+1]* marketModel->evolution().rateTaus()[m];
2754 Real expiry = rateTimes[m];
2755 Real forward = todaysForwards[m];
2756
2757 priceConstVol += blackFormula(Option::Call,
2758 capStrike,
2759 forward,
2760 impVol*sqrt(expiry),
2761 annuity,
2762 marketModel->displacements()[m]);
2763
2764 priceVarVol += blackFormula(Option::Call,
2765 capStrike,
2766 forward,
2767 sqrt(totalCov[m][m]),
2768 annuity,
2769 marketModel->displacements()[m]);
2770
2771 }
2772
2773 if (fabs(priceVarVol - priceConstVol) > impVolTolerance)
2774 ++numberCapImpVolFailures;
2775
2776
2777 }
2778
2779 if (numberCapPseudoFailures >0)
2780 BOOST_ERROR("cap pseudo test failed for prices "
2781 << numberCapPseudoFailures << " times" );
2782
2783 if (numberCapImpVolFailures >0)
2784 BOOST_ERROR("cap pseudo test failed for implied vols "
2785 << numberCapImpVolFailures << " times" );
2786
2787 }
2788
2789 // we have tested the price derivative and the implied vol function, now the derivative of the cap implied vols
2790 // with respect to pseudo-root elements
2791
2792 // since we have already tested the imp vol function we use it here
2793
2794
2795 for (Size m=0; m<LENGTH(testedFactors); ++m)
2796 {
2797 Size factors = testedFactors[m];
2798
2799
2800 bool logNormal = true;
2801
2802 ext::shared_ptr<MarketModel> marketModel =
2803 makeMarketModel(logNormal, evolution, factors,
2804 marketModels[j]);
2805
2806 for (Size startIndex = 1; startIndex < evolution.numberOfRates()-1; ++startIndex)
2807 for (Size endIndex = startIndex+1; endIndex < evolution.numberOfRates(); ++endIndex)
2808 {
2809
2810 CapPseudoDerivative derivative(marketModel,
2811 capStrike,
2812 startIndex,
2813 endIndex,initialNumeraireValue);
2814
2815 std::vector<Matrix> pseudoRoots;
2816 for (Size k=0; k < marketModel->numberOfSteps(); ++k)
2817 pseudoRoots.push_back( marketModel->pseudoRoot(k));
2818
2819
2820 for (Size step=0; step < evolution.numberOfSteps(); ++ step)
2821 {
2822 for (Size l=0; l < evolution.numberOfRates(); ++l)
2823 for (Size f=0; f < factors; ++f)
2824 {
2825 pseudoRoots[step][l][f] += numericalBumpSizeForSwaptionPseudo;
2826
2827 PseudoRootFacade bumpedUp(pseudoRoots,rateTimes,marketModel->initialRates(),marketModel->displacements());
2828
2829 CapPseudoDerivative upDerivative(ext::shared_ptr<MarketModel>(new PseudoRootFacade(bumpedUp)),
2830 capStrike,
2831 startIndex,
2832 endIndex,initialNumeraireValue);
2833
2834 Real volUp = upDerivative.impliedVolatility();
2835
2836
2837
2838
2839 pseudoRoots[step][l][f] -= numericalBumpSizeForSwaptionPseudo;
2840
2841 pseudoRoots[step][l][f] -= numericalBumpSizeForSwaptionPseudo;
2842
2843 PseudoRootFacade bumpedDown(pseudoRoots,rateTimes,marketModel->initialRates(),marketModel->displacements());
2844
2845 CapPseudoDerivative downDerivative(ext::shared_ptr<MarketModel>(new PseudoRootFacade(bumpedDown)),
2846 capStrike,
2847 startIndex,
2848 endIndex,initialNumeraireValue);
2849
2850
2851 Real volDown = downDerivative.impliedVolatility();
2852
2853
2854
2855
2856 pseudoRoots[step][l][f] += numericalBumpSizeForSwaptionPseudo;
2857
2858
2859 Real volDeriv = (volUp-volDown)/(2.0*numericalBumpSizeForSwaptionPseudo);
2860
2861 Real modelVal = derivative.volatilityDerivative(step)[l][f];
2862
2863 Real error = volDeriv - modelVal;
2864
2865 if (fabs(error) > impVolTolerance*10)
2866 ++numberCapVolPseudoFailures;
2867
2868
2869
2870 }
2871
2872
2873
2874 }
2875
2876
2877
2878 }
2879
2880 if (numberCapVolPseudoFailures >0)
2881 BOOST_ERROR("cap pseudo test failed for implied vols "
2882 << numberCapVolPseudoFailures << " times" );
2883
2884 }
2885 }
2886
2887
2888
2889
2890
2891
2892
2893
2894 /////////////////////////////////////
2895
2896 for (Size j=0; j<LENGTH(marketModels); j++)
2897 {
2898
2899 Size testedFactors[] = {
2900 std::min<Size>(1UL,todaysForwards.size())
2901 // todaysForwards.size()
2902 //, 4, 8,
2903
2904 };
2905
2906
2907
2908
2909 for (Size m=0; m<LENGTH(testedFactors); ++m)
2910 {
2911 Size factors = testedFactors[m];
2912 Size factorsToTest = std::min<Size>(2,factors); // doing all possible vegas is combinatorially explosive
2913
2914
2915 MeasureType measures[] = {
2916 MoneyMarket
2917 };
2918
2919 std::vector<Matrix> pseudoBumps;
2920 std::vector<Matrix> pseudoBumpsDown;
2921
2922 for (Size k=0; k < evolution.numberOfRates(); ++k)
2923 {
2924 for (Size f=0; f < factors; ++f)
2925 {
2926 Matrix modelBump(evolution.numberOfRates(), factors,0.0);
2927 modelBump[k][f] =bumpSizeNumericalDifferentiation;
2928 pseudoBumps.push_back(modelBump);
2929 modelBump[k][f] =-bumpSizeNumericalDifferentiation;
2930 pseudoBumpsDown.push_back(modelBump);
2931 }
2932 }
2933
2934 std::vector<std::vector<Matrix> > vegaBumps;
2935
2936 Matrix modelBump(evolution.numberOfRates(), factors,0.0);
2937
2938
2939 for (Size l = 0; l < evolution.numberOfSteps(); ++l)
2940 {
2941 vegaBumps.push_back(std::vector<Matrix>());
2942 for (Size k=0; k < evolution.numberOfRates(); k=k+bumpIncrement)
2943 {
2944 for (Size f=0; f < factorsToTest; ++f)
2945 {
2946
2947 for (Size m=0; m < evolution.numberOfSteps(); ++m)
2948 {
2949 if (l ==m && k >= l)
2950 modelBump[k][f] = vegaBumpSize;
2951
2952 vegaBumps[l].push_back(modelBump);
2953
2954 modelBump[k][f] =0.0;
2955 }
2956 }
2957 }
2958
2959 }
2960
2961
2962
2963
2964
2965 for (Size k=0; k<LENGTH(measures); k++)
2966 {
2967
2968 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
2969
2970 std::vector<RatePseudoRootJacobian> testees;
2971 std::vector<RatePseudoRootJacobianAllElements> testees2;
2972
2973 std::vector<RatePseudoRootJacobianNumerical> testers;
2974 std::vector<RatePseudoRootJacobianNumerical> testersDown;
2975
2976
2977 MTBrownianGeneratorFactory generatorFactory(seed_);
2978
2979 bool logNormal = true;
2980 ext::shared_ptr<MarketModel> marketModel =
2981 makeMarketModel(logNormal, evolution, factors,
2982 marketModels[j]);
2983
2984 for (Size l=0; l < evolution.numberOfSteps(); ++l)
2985 {
2986 const Matrix& pseudoRoot = marketModel->pseudoRoot(l);
2987 testees.push_back(RatePseudoRootJacobian(pseudoRoot,
2988 evolution.firstAliveRate()[l],
2989 numeraires[l],
2990 evolution.rateTaus(),
2991 pseudoBumps,
2992 marketModel->displacements()
2993 ));
2994
2995 testees2.push_back(RatePseudoRootJacobianAllElements(pseudoRoot,
2996 evolution.firstAliveRate()[l],
2997 numeraires[l],
2998 evolution.rateTaus(),
2999 marketModel->displacements()
3000 ));
3001
3002
3003 testers.push_back(RatePseudoRootJacobianNumerical(pseudoRoot,
3004 evolution.firstAliveRate()[l],
3005 numeraires[l],
3006 evolution.rateTaus(),
3007 pseudoBumps,
3008 marketModel->displacements()
3009 ));
3010 testersDown.push_back(RatePseudoRootJacobianNumerical(pseudoRoot,
3011 evolution.firstAliveRate()[l],
3012 numeraires[l],
3013 evolution.rateTaus(),
3014 pseudoBumpsDown,
3015 marketModel->displacements()
3016 ));
3017
3018 }
3019
3020
3021
3022
3023 ext::shared_ptr<BrownianGenerator> generator(generatorFactory.create(factors,
3024 steps));
3025 LogNormalFwdRateEuler evolver(marketModel,
3026 generatorFactory,
3027 numeraires);
3028
3029
3030 std::vector<Real> oldRates(evolution.numberOfRates());
3031 std::vector<Real> newRates(evolution.numberOfRates());
3032 std::vector<Real> gaussians(factors);
3033
3034 std::vector<Size> numberCashFlowsThisStep(product.numberOfProducts());
3035
3036 std::vector<std::vector<MarketModelMultiProduct::CashFlow> > cashFlowsGenerated(product.numberOfProducts());
3037
3038 for (Size i=0; i < product.numberOfProducts(); ++i)
3039 cashFlowsGenerated[i].resize(product.maxNumberOfCashFlowsPerProductPerStep());
3040
3041 Matrix B(pseudoBumps.size(),evolution.numberOfRates());
3042 Matrix B2(pseudoBumps.size(),evolution.numberOfRates());
3043 Matrix B3(pseudoBumps.size(),evolution.numberOfRates());
3044 Matrix B4(pseudoBumps.size(),evolution.numberOfRates());
3045
3046 std::vector<Matrix> globalB;
3047 {
3048 Matrix modelB(evolution.numberOfRates(), factors);
3049 for (Size i=0; i < steps; ++i)
3050 globalB.push_back(modelB);
3051 }
3052
3053 std::vector<Real> oneStepDFs(evolution.numberOfRates()+1);
3054 oneStepDFs[0] = 1.0;
3055
3056
3057 Size numberFailures=0;
3058 Size numberFailures2=0;
3059
3060 for (Size l=0; l < pathsToDo; ++l)
3061 {
3062 evolver.startNewPath();
3063 product.reset();
3064 generator->nextPath();
3065
3066 bool done;
3067 newRates = marketModel->initialRates();
3068 Size currentStep =0;
3069
3070 do
3071 {
3072 oldRates = newRates;
3073
3074
3075 evolver.advanceStep();
3076 done = product.nextTimeStep(evolver.currentState(),
3077 numberCashFlowsThisStep,
3078 cashFlowsGenerated);
3079
3080 newRates = evolver.currentState().forwardRates();
3081
3082 for (Size i=1; i <= evolution.numberOfRates(); ++i)
3083 oneStepDFs[i] = 1.0/(1+oldRates[i-1]*evolution.rateTaus()[i-1]);
3084
3085
3086 generator->nextStep(gaussians);
3087
3088 testees[currentStep].getBumps(oldRates, oneStepDFs, newRates, gaussians, B);
3089 testees2[currentStep].getBumps(oldRates, oneStepDFs, newRates, gaussians, globalB);
3090
3091
3092 testers[currentStep].getBumps(oldRates, oneStepDFs, newRates, gaussians, B2);
3093 testersDown[currentStep].getBumps(oldRates, oneStepDFs, newRates, gaussians, B3);
3094
3095 // now do make out put of allElements class into same form
3096
3097 for (Size i1 =0; i1 < pseudoBumps.size(); ++i1)
3098 {
3099 Size j1=0;
3100
3101 for (; j1 < evolution.firstAliveRate()[i1]; ++j1)
3102 {
3103 B4[i1][j1]=0.0;
3104 }
3105 for (; j1 < numberRates; ++j1)
3106 {
3107 Real sum =0.0;
3108
3109 for (Size k1=evolution.firstAliveRate()[i1]; k1 < numberRates; ++k1)
3110 for (Size f1=0; f1 < factors; ++f1)
3111 sum += pseudoBumps[i1][k1][f1]*globalB[j1][k1][f1];
3112
3113 B4[i1][j1] =sum;
3114
3115 }
3116 }
3117
3118
3119
3120 for (Size j=0; j < B.rows(); ++j)
3121 for (Size k=0; k < B.columns(); ++k)
3122 {
3123 Real analytic = B[j][k]/bumpSizeNumericalDifferentiation;
3124 Real analytic2 = B4[j][k]/bumpSizeNumericalDifferentiation;
3125 Real numerical = (B2[j][k]-B3[j][k])/(2*bumpSizeNumericalDifferentiation);
3126 Real errorSize = (analytic - numerical)/ ( bumpSizeNumericalDifferentiation*bumpSizeNumericalDifferentiation);
3127 Real errorSize2 = (analytic2 - numerical)/ ( bumpSizeNumericalDifferentiation*bumpSizeNumericalDifferentiation);
3128
3129 maxError = std::max(maxError,fabs(errorSize));
3130
3131 if ( fabs( errorSize ) > multiplier )
3132 {
3133 ++numberFailures;
3134 if (printReport_)
3135 BOOST_TEST_MESSAGE("path " << l << " step "
3136 << currentStep << " j " << j
3137 << " k " << k << " B " << B[j][k] << " B2 " << B2[j][k]);
3138
3139 }
3140
3141 if ( fabs( errorSize2 ) > multiplier )
3142 {
3143 ++numberFailures2;
3144 if (printReport_)
3145 BOOST_TEST_MESSAGE("path " << l << " step "
3146 << currentStep << " j " << j
3147 << " k " << k << " B4 " << B4[j][k] << " B2 " << B2[j][k]);
3148
3149 }
3150
3151 }
3152 ++currentStep;
3153 }
3154 while (!done);
3155
3156 }
3157
3158 if (numberFailures >0)
3159 BOOST_FAIL("Pathwise rate pseudoroot jacobian test fails : " << numberFailures <<"\n");
3160
3161
3162 if (numberFailures2 >0)
3163 BOOST_FAIL("Pathwise rate pseudoroot jacobian all elements test fails : " << numberFailures2 <<"\n");
3164 } // end of k loop over measures
3165
3166
3167 // the quick test done now do a simulation test for the vegas for caplets
3168
3169 Size numberDeflatedErrors =0;
3170 Size numberUndeflatedErrors =0;
3171 Real biggestError=0.0;
3172
3173
3174 for (Size deflate =0; deflate <2; ++deflate)
3175 {
3176 Clone<MarketModelPathwiseMultiProduct> productToUse;
3177
3178 if (deflate ==0)
3179 productToUse = caplets;
3180 else
3181 productToUse = capletsDeflated;
3182
3183 for (Size k=0; k<LENGTH(measures); k++)
3184 {
3185
3186 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
3187
3188 MTBrownianGeneratorFactory generatorFactory(seed_);
3189
3190 bool logNormal = true;
3191 ext::shared_ptr<MarketModel> marketModel =
3192 makeMarketModel(logNormal, evolution, factors,
3193 marketModels[j]);
3194
3195 LogNormalFwdRateEuler evolver(marketModel,
3196 generatorFactory,
3197 numeraires);
3198
3199 // SequenceStatistics stats(product.numberOfProducts()*(todaysForwards.size()+1+vegaBumps[0].size()));
3200
3201
3202 std::ostringstream config;
3203 config <<
3204 marketModelTypeToString(marketModels[j]) << ", " <<
3205 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
3206 measureTypeToString(measures[k]) << ", " <<
3207 "MT BGF";
3208 if (printReport_)
3209 BOOST_TEST_MESSAGE(" " << config.str());
3210
3211 Size initialNumeraire = evolver.numeraires().front();
3212 Real initialNumeraireValue =
3213 todaysDiscounts[initialNumeraire];
3214
3215 std::vector<Real> values;
3216
3217 std::vector<Real> errors;
3218
3219 {
3220
3221 PathwiseVegasAccountingEngine accountingengine(ext::make_shared<LogNormalFwdRateEuler>(evolver), // method relies heavily on LMM Euler
3222 productToUse,
3223 marketModel, // we need pseudo-roots and displacements
3224 vegaBumps,
3225 initialNumeraireValue);
3226
3227 accountingengine.multiplePathValues(values,errors,pathsToDoSimulation);
3228 }
3229
3230 // we have computed the vegas now we have to test them against the analytic values
3231
3232 // extract into easier format
3233
3234
3235
3236
3237 Matrix vegasMatrix(caplets.numberOfProducts(), vegaBumps[0].size());
3238 Matrix standardErrors(vegasMatrix);
3239 Matrix deltasMatrix(caplets.numberOfProducts(), numberRates);
3240 Matrix deltasErrors(deltasMatrix);
3241 std::vector<Real> prices(caplets.numberOfProducts());
3242 std::vector<Real> priceErrors(caplets.numberOfProducts());
3243
3244 Size entriesPerProduct = 1+numberRates+vegaBumps[0].size();
3245
3246 for (Size i=0; i < caplets.numberOfProducts(); ++i)
3247 {
3248 prices[i] = values[i*entriesPerProduct];
3249 priceErrors[i] = errors[i*entriesPerProduct];
3250
3251 for (Size j=0; j < vegaBumps[0].size(); ++j)
3252 {
3253 vegasMatrix[i][j] = values[i*entriesPerProduct + numberRates+1 + j];
3254 standardErrors[i][j] = errors[i*entriesPerProduct + numberRates+1 + j];
3255 }
3256 for (Size j=0; j < numberRates; ++j)
3257 {
3258 deltasMatrix[i][j] = values[i*entriesPerProduct +1 + j];
3259 deltasErrors[i][j] = errors[i*entriesPerProduct +1 + j];
3260 }
3261 }
3262
3263
3264
3265 // first get the terminal vols
3266
3267 Matrix totalCovariance(marketModel->totalCovariance(marketModel->numberOfSteps()-1));
3268
3269
3270 std::vector<Real> truePrices(caplets.numberOfProducts());
3271
3272 for (Size r =0; r < truePrices.size(); ++r)
3273 {
3274 truePrices[r] = BlackCalculator(displacedPayoffs[r], todaysForwards[r], sqrt(totalCovariance[r][r]),
3275 todaysDiscounts[r+1]*(rateTimes[r+1]-rateTimes[r])).value();
3276 }
3277
3278
3279 for (Size b =0; b < vegaBumps[0].size(); ++b)
3280 {
3281
3282
3283 std::vector<Real> bumpedPrices(truePrices.size());
3284 std::vector<Real> variances(truePrices.size(),0.0);
3285 std::vector<Real> vegas(truePrices.size());
3286
3287
3288 for (Size step = 0; step < marketModel->numberOfSteps(); ++step)
3289 {
3290 Matrix pseudoRoot( marketModel->pseudoRoot(step));
3291 pseudoRoot += vegaBumps[step][b];
3292
3293 for (Size rate=step; rate<marketModel->numberOfRates(); ++rate)
3294 {
3295 Real variance = 0.0;
3296 for (Size f=0; f < marketModel->numberOfFactors(); ++f)
3297 variance+= pseudoRoot[rate][f]* pseudoRoot[rate][f];
3298
3299 variances[rate]+=variance;
3300 }
3301 }
3302
3303 for (Size r =0; r < truePrices.size(); ++r)
3304 {
3305
3306 bumpedPrices[r] = BlackCalculator(displacedPayoffs[r], todaysForwards[r], sqrt(variances[r]),
3307 todaysDiscounts[r+1]*(rateTimes[r+1]-rateTimes[r])).value();
3308
3309 vegas[r] = bumpedPrices[r] - truePrices[r];
3310
3311 }
3312
3313
3314 for (Size s=0; s < truePrices.size(); ++s)
3315 {
3316 Real mcVega = vegasMatrix[s][b];
3317 Real analyticVega = vegas[s];
3318 Real thisError = mcVega - analyticVega;
3319 Real thisSE = standardErrors[s][b];
3320
3321 if (fabs(thisError) > 0.0)
3322 {
3323 Real errorInSEs = thisError/thisSE;
3324 biggestError = std::max(fabs(errorInSEs),biggestError);
3325
3326 if (fabs(errorInSEs) > 4.5)
3327 {
3328 if (deflate==0)
3329 ++numberUndeflatedErrors;
3330 else
3331 ++numberDeflatedErrors;
3332 }
3333 }
3334
3335 }
3336
3337
3338 }
3339
3340
3341
3342 // for deltas and prices the pathwise vega engine should agree precisely with the pathwiseaccounting engine
3343 // so lets see if it does
3344
3345 Clone<MarketModelPathwiseMultiProduct> productToUse2;
3346
3347 if (deflate ==0)
3348 productToUse2 = caplets;
3349 else
3350 productToUse2 = capletsDeflated;
3351
3352
3353 SequenceStatisticsInc stats(productToUse2->numberOfProducts()*(todaysForwards.size()+1));
3354 {
3355 PathwiseAccountingEngine accountingengine(ext::make_shared<LogNormalFwdRateEuler>(evolver), // method relies heavily on LMM Euler
3356 *productToUse2,
3357 marketModel, // we need pseudo-roots and displacements
3358 initialNumeraireValue);
3359
3360 accountingengine.multiplePathValues(stats,pathsToDoSimulation);
3361 }
3362
3363 std::vector<Real> valuesAndDeltas2 = stats.mean();
3364 std::vector<Real> errors2 = stats.errorEstimate();
3365
3366 std::vector<Real> prices2(productToUse2->numberOfProducts());
3367 std::vector<Real> priceErrors2(productToUse2->numberOfProducts());
3368
3369 Matrix deltas2( productToUse2->numberOfProducts(), todaysForwards.size());
3370 Matrix deltasErrors2( productToUse2->numberOfProducts(), todaysForwards.size());
3371 std::vector<Real> modelPrices2(productToUse2->numberOfProducts());
3372
3373
3374 for (Size i=0; i < productToUse2->numberOfProducts(); ++i)
3375 {
3376 prices2[i] = valuesAndDeltas2[i];
3377 priceErrors2[i] = errors2[i];
3378
3379 for (Size j=0; j < todaysForwards.size(); ++j)
3380 {
3381 deltas2[i][j] = valuesAndDeltas2[(i+1)*productToUse2->numberOfProducts()+j];
3382 deltasErrors2[i][j] = errors2[(i+1)* productToUse2->numberOfProducts()+j];
3383 }
3384 }
3385
3386 for (Size i=0; i < productToUse2->numberOfProducts(); ++i)
3387 {
3388
3389 Real priceDiff = prices2[i] - prices[i];
3390
3391 if (fabs(priceDiff) > 5*priceErrors2[i]) // two sets of standard error
3392 BOOST_FAIL("pathwise accounting engine and pathwise vegas accounting engine not in perfect agreement for price.\n product " << i << ", vega computed price: " << prices[j] << " previous price " << prices2[j] << ", deflate " << deflate << "\n" );
3393
3394 for (Size j=0; j < todaysForwards.size(); ++j)
3395 {
3396 Real error = deltas2[i][j] - deltasMatrix[i][j];
3397 if (fabs(error)> 5* deltasErrors2[i][j] ) // two sets of standard error
3398 BOOST_FAIL("pathwise accounting engine and pathwise vegas accounting engine not in perfect agreement for dealts.\n product " << i << ", rate " << j << " vega computed delta: " << deltasMatrix[i][j] << " previous delta " << deltas2[i][j] << "\n" );
3399 }
3400 }
3401 } // end of k loop over measures
3402 } // end of loop over deflation
3403
3404
3405 if (numberDeflatedErrors+numberUndeflatedErrors >0)
3406 BOOST_FAIL("Model pathwise vega test for caplets fails : " << numberDeflatedErrors <<" deflated errors and " <<numberUndeflatedErrors << " undeflated errors , biggest error in SEs is " << biggestError << "\n");
3407
3408
3409 {
3410 // now do a simulation test for the vegas for caps
3411
3412 std::vector<VolatilityBumpInstrumentJacobian::Cap> caps;
3413
3414 Rate capStrike = todaysForwards[0];
3415
3416 for (Size i=0; i +2 < numberRates; i=i+3)
3417 {
3418 VolatilityBumpInstrumentJacobian::Cap nextCap;
3419 // nextCap.startIndex_ = i;
3420 // nextCap.endIndex_ = i+3;
3421 // nextCap.strike_ = capStrike;
3422 // caps.push_back(nextCap);
3423
3424 // nextCap.startIndex_ = i+1;
3425 // nextCap.endIndex_ = i+3;
3426 // nextCap.strike_ = capStrike;
3427 // caps.push_back(nextCap);
3428
3429 nextCap.startIndex_ = i+2;
3430 nextCap.endIndex_ = i+3;
3431 nextCap.strike_ = capStrike;
3432 caps.push_back(nextCap);
3433
3434 }
3435
3436 std::vector<std::pair<Size,Size> > startsAndEnds(caps.size());
3437
3438 for (Size r=0; r < caps.size(); ++r)
3439 {
3440 startsAndEnds[r].first = caps[r].startIndex_;
3441 startsAndEnds[r].second = caps[r].endIndex_;
3442 }
3443
3444 MarketModelPathwiseMultiDeflatedCap capsDeflated(
3445 rateTimes,
3446 accruals,
3447 paymentTimes,
3448 capStrike,
3449 startsAndEnds);
3450
3451 for (Size k=0; k<LENGTH(measures); k++)
3452 {
3453
3454 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
3455
3456 MTBrownianGeneratorFactory generatorFactory(seed_);
3457 MTBrownianGeneratorFactory generatorFactory2(seed_);
3458
3459 bool logNormal = true;
3460 ext::shared_ptr<MarketModel> marketModel =
3461 makeMarketModel(logNormal, evolution, factors,
3462 marketModels[j]);
3463
3464 LogNormalFwdRateEuler evolver(marketModel,
3465 generatorFactory,
3466 numeraires);
3467
3468 LogNormalFwdRateEuler evolver2(marketModel,
3469 generatorFactory2,
3470 numeraires);
3471
3472 // SequenceStatistics stats(product.numberOfProducts()*(todaysForwards.size()+1+vegaBumps[0].size()));
3473
3474
3475 std::ostringstream config;
3476 config <<
3477 marketModelTypeToString(marketModels[j]) << ", " <<
3478 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
3479 measureTypeToString(measures[k]) << ", " <<
3480 "MT BGF";
3481 if (printReport_)
3482 BOOST_TEST_MESSAGE(" " << config.str());
3483
3484 Size initialNumeraire = evolver.numeraires().front();
3485 Real initialNumeraireValue =
3486 todaysDiscounts[initialNumeraire];
3487
3488 std::vector<Real> values;
3489 std::vector<Real> errors;
3490
3491 std::vector<Real> values2;
3492 std::vector<Real> errors2;
3493
3494
3495 {
3496
3497 PathwiseVegasOuterAccountingEngine accountingengine(ext::make_shared<LogNormalFwdRateEuler>(evolver2), // method relies heavily on LMM Euler
3498 capsDeflated,
3499 marketModel, // we need pseudo-roots and displacements
3500 vegaBumps,
3501 initialNumeraireValue);
3502
3503 accountingengine.multiplePathValues(values2,errors2,pathsToDoSimulation);
3504 }
3505
3506 {
3507
3508 PathwiseVegasAccountingEngine accountingengine(ext::make_shared<LogNormalFwdRateEuler>(evolver), // method relies heavily on LMM Euler
3509 capsDeflated,
3510 marketModel, // we need pseudo-roots and displacements
3511 vegaBumps,
3512 initialNumeraireValue);
3513
3514 accountingengine.multiplePathValues(values,errors,pathsToDoSimulation);
3515 }
3516
3517 // first test to see that the two implementation give the same results
3518
3519 {
3520 Real tol = 1E-8;
3521
3522 Size numberMeanFailures =0;
3523
3524 for (Size i=0; i <values.size(); ++i)
3525 if (fabs(values[i]-values2[i]) > tol)
3526 ++numberMeanFailures;
3527
3528 if (numberMeanFailures >0)
3529 BOOST_FAIL("Comparison of Pathwise vegas accounting engine and PathwiseVegasOuterAccountingEngine yields discrepancies:"
3530 << numberMeanFailures
3531 << " out of "
3532 << values.size() );
3533
3534 }
3535
3536 // we have computed the vegas now we have to test them against the analytic values
3537
3538 // extract into easier format
3539
3540
3541
3542
3543 Matrix vegasMatrix(capsDeflated.numberOfProducts(), vegaBumps[0].size());
3544 Matrix standardErrors(vegasMatrix);
3545 Size entriesPerProduct = 1+numberRates+vegaBumps[0].size();
3546
3547 for (Size i=0; i < capsDeflated.numberOfProducts(); ++i)
3548 for (Size j=0; j < vegaBumps[0].size(); ++j)
3549 {
3550 vegasMatrix[i][j] = values[i*entriesPerProduct + numberRates+1 + j];
3551 standardErrors[i][j] = errors[i*entriesPerProduct + numberRates+1 + j];
3552 }
3553
3554
3555 // first get the terminal vols
3556
3557 Matrix totalCovariance(marketModel->totalCovariance(marketModel->numberOfSteps()-1));
3558
3559 std::vector<Real> trueCapletPrices(numberRates);
3560 ext::shared_ptr<StrikedTypePayoff> dispayoff( new
3561 PlainVanillaPayoff(Option::Call, capStrike+displacement));
3562
3563 for (Size r =0; r < trueCapletPrices.size(); ++r)
3564 trueCapletPrices[r] = BlackCalculator(dispayoff, todaysForwards[r], sqrt(totalCovariance[r][r]),
3565 todaysDiscounts[r+1]*(rateTimes[r+1]-rateTimes[r])).value();
3566
3567 std::vector<Real> trueCapPrices(capsDeflated.numberOfProducts());
3568 std::vector<Real> vegaCaps(capsDeflated.numberOfProducts());
3569
3570
3571 for (Size s=0; s < capsDeflated.numberOfProducts(); ++s)
3572 {
3573
3574 trueCapPrices[s]=0.0;
3575
3576 for (Size t= caps[s].startIndex_; t < caps[s].endIndex_; ++t)
3577 trueCapPrices[s] += trueCapletPrices[t];
3578 }
3579
3580 Size numberErrors =0;
3581
3582
3583 for (Size b =0; b < vegaBumps[0].size(); ++b)
3584 {
3585
3586 std::vector<Real> bumpedCapletPrices(trueCapletPrices.size());
3587 // std::vector<Real> bumpedCapPrices(trueCapPrices.size());
3588
3589 std::vector<Real> variances(trueCapletPrices.size(),0.0);
3590 std::vector<Real> vegasCaplets(trueCapletPrices.size());
3591
3592 for (Size step = 0; step < marketModel->numberOfSteps(); ++step)
3593 {
3594 Matrix pseudoRoot( marketModel->pseudoRoot(step));
3595 pseudoRoot += vegaBumps[step][b];
3596
3597 for (Size rate=step; rate<marketModel->numberOfRates(); ++rate)
3598 {
3599 Real variance = 0.0;
3600 for (Size f=0; f < marketModel->numberOfFactors(); ++f)
3601 variance+= pseudoRoot[rate][f]* pseudoRoot[rate][f];
3602
3603 variances[rate]+=variance;
3604 }
3605 }
3606
3607 for (Size r =0; r < trueCapletPrices.size(); ++r)
3608 {
3609 bumpedCapletPrices[r] = BlackCalculator(dispayoff, todaysForwards[r], sqrt(variances[r]),
3610 todaysDiscounts[r+1]*(rateTimes[r+1]-rateTimes[r])).value();
3611
3612 vegasCaplets[r] = bumpedCapletPrices[r] - trueCapletPrices[r];
3613 }
3614
3615 for (Size s=0; s < capsDeflated.numberOfProducts(); ++s)
3616 {
3617 vegaCaps[s]=0.0;
3618
3619 for (Size t= caps[s].startIndex_; t < caps[s].endIndex_; ++t)
3620 vegaCaps[s] += vegasCaplets[t];
3621 }
3622
3623 for (Size s=0; s < capsDeflated.numberOfProducts(); ++s)
3624 {
3625 Real mcVega = vegasMatrix[s][b];
3626 Real analyticVega = vegaCaps[s];
3627 Real thisError = mcVega - analyticVega;
3628 Real thisSE = standardErrors[s][b];
3629
3630 if (fabs(thisError) > 0.0)
3631 {
3632 Real errorInSEs = fabs(thisError/thisSE);
3633
3634 if (errorInSEs > 4.0)
3635 ++numberErrors;
3636 }
3637
3638 }
3639
3640 }
3641
3642
3643 if (numberErrors >0)
3644 BOOST_FAIL("caps Pathwise vega test fails : " << numberErrors <<"\n");
3645
3646 } // end of k loop over measures
3647 }
3648
3649 }
3650 }
3651
3652 }
3653
testPathwiseMarketVegas()3654 void MarketModelTest::testPathwiseMarketVegas()
3655 {
3656
3657 BOOST_TEST_MESSAGE("Testing pathwise market vegas in a lognormal forward rate market model...");
3658
3659 using namespace market_model_test;
3660
3661 setup();
3662
3663 // specify collection of caps and swaptions and then see if their vegas are correct
3664 // starting by doing a set of co-terminal swaptions
3665 LMMCurveState cs(rateTimes);
3666 cs.setOnForwardRates(todaysForwards);
3667
3668 std::vector<ext::shared_ptr<Payoff> > payoffs(todaysForwards.size());
3669 std::vector<ext::shared_ptr<StrikedTypePayoff> >
3670 displacedPayoffs(todaysForwards.size());
3671 for (Size i=0; i<todaysForwards.size(); ++i) {
3672 payoffs[i] = ext::shared_ptr<Payoff>(new
3673 PlainVanillaPayoff(Option::Call, cs.coterminalSwapRate(i)));
3674
3675 displacedPayoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
3676 PlainVanillaPayoff(Option::Call, cs.coterminalSwapRate(i)+displacement));
3677
3678 }
3679
3680
3681 MultiStepOptionlets dummyProduct(rateTimes, accruals,
3682 paymentTimes, payoffs);
3683
3684 Real bumpSizeNumericalDifferentiation = 1E-6;
3685
3686 MarketModelPathwiseCoterminalSwaptionsDeflated swaptionsDeflated(rateTimes, cs.coterminalSwapRates());
3687 MarketModelPathwiseCoterminalSwaptionsNumericalDeflated swaptionsDeflated2(rateTimes, cs.coterminalSwapRates(),bumpSizeNumericalDifferentiation);
3688
3689
3690 const EvolutionDescription& evolution = dummyProduct.evolution();
3691 Size steps = evolution.numberOfSteps();
3692 Size numberRates = evolution.numberOfRates();
3693
3694
3695 Size pathsToDo =10; // for the numerical differentiation test we are requiring equality on each path so this is actually quite strict
3696 Size pathsToDoSimulation = paths_;
3697
3698 Real multiplier = 50; // how many times the bump size squared, the numerical differentation is allowed to differ by
3699 Real tolerance = 1E-6;
3700
3701 // printReport_ = true;
3702 Real initialNumeraireValue =0.95;
3703
3704 bool allowFactorwiseBumping = true;
3705 std::vector<VolatilityBumpInstrumentJacobian::Cap> caps;
3706
3707 Rate capStrike = todaysForwards[0];
3708
3709 for (Size i=0; i +2 < numberRates; i=i+3)
3710 {
3711 VolatilityBumpInstrumentJacobian::Cap nextCap;
3712 nextCap.startIndex_ = i;
3713 nextCap.endIndex_ = i+3;
3714 nextCap.strike_ = capStrike;
3715 caps.push_back(nextCap);
3716 }
3717 std::vector<std::pair<Size,Size> > startsAndEnds(caps.size());
3718
3719
3720 for (Size j=0; j < caps.size(); ++j)
3721 {
3722 startsAndEnds[j].first = caps[j].startIndex_;
3723 startsAndEnds[j].second = caps[j].endIndex_;
3724
3725
3726 }
3727
3728
3729 MarketModelPathwiseMultiDeflatedCap capsDeflated(
3730 rateTimes,
3731 accruals,
3732 paymentTimes,
3733 capStrike,
3734 startsAndEnds);
3735
3736
3737
3738 std::vector<VolatilityBumpInstrumentJacobian::Swaption> swaptions(numberRates);
3739
3740 for (Size i=0; i < numberRates; ++i)
3741 {
3742 swaptions[i].startIndex_ = i;
3743 swaptions[i].endIndex_ = numberRates;
3744
3745 }
3746
3747
3748
3749 MarketModelType marketModels[] =
3750 {
3751 // CalibratedMM,
3752 // ExponentialCorrelationFlatVolatility,
3753 ExponentialCorrelationAbcdVolatility
3754 };
3755 ///////////////////////////////////
3756 ///////////////////////////////////
3757 // test analytically first, it's faster!
3758
3759 for (Size j=0; j<LENGTH(marketModels); j++)
3760 {
3761
3762 Size testedFactors[] = { std::min<Size>(1UL,todaysForwards.size())
3763 // todaysForwards.size()
3764 //, 4, 8,
3765 };
3766
3767
3768
3769 for (Size m=0; m<LENGTH(testedFactors); ++m)
3770 {
3771 Size factors = testedFactors[m];
3772
3773
3774 bool logNormal = true;
3775
3776 ext::shared_ptr<MarketModel> marketModel =
3777 makeMarketModel(logNormal, evolution, factors,
3778 marketModels[j]);
3779
3780
3781 // we need to work out our bumps
3782
3783 VegaBumpCollection possibleBumps(marketModel,
3784 allowFactorwiseBumping);
3785
3786 OrthogonalizedBumpFinder bumpFinder(possibleBumps,
3787 swaptions,
3788 caps,
3789 multiplier, // if vector length grows by more than this discard
3790 tolerance); // if vector projection before scaling less than this discard
3791 std::vector<std::vector<Matrix> > theBumps;
3792
3793 bumpFinder.GetVegaBumps(theBumps);
3794
3795 // the bumps is now the bumps required to get a one percent implied vol in each instrumnet
3796 // indexed by step, instrument, pseudo-root matrix
3797 // if we dot product with swaption derivatives, we should get a 1% change in imp vol on the diagonal
3798 // and zero off it
3799 {
3800 Matrix swaptionVegasMatrix(swaptionsDeflated.numberOfProducts(), theBumps[0].size());
3801
3802 for (Size i=0; i < swaptionsDeflated.numberOfProducts(); ++i)
3803 {
3804 SwaptionPseudoDerivative thisPseudoDerivative(marketModel,
3805 swaptions[i].startIndex_,
3806 swaptions[i].endIndex_);
3807
3808
3809 for (Size j=0; j < theBumps[0].size(); ++j)
3810 {
3811 swaptionVegasMatrix[i][j] = 0;
3812
3813 for (Size k=0; k < steps; ++k)
3814 for (Size l=0; l < numberRates; ++l)
3815 for (Size m=0; m < factors; ++m)
3816 swaptionVegasMatrix[i][j] += theBumps[k][j][l][m]*thisPseudoDerivative.volatilityDerivative(k)[l][m];
3817 }
3818 }
3819
3820 Size numberDiagonalFailures = 0;
3821 Size offDiagonalFailures=0;
3822
3823 for (Size i=0; i < swaptions.size(); ++i)
3824 {
3825 for (Size j=0; j < theBumps[0].size(); ++j)
3826 {
3827 if (i == j)
3828 {
3829 Real thisError = swaptionVegasMatrix[i][i] - 0.01;
3830
3831 if (fabs(thisError) > 1e-8)
3832 ++numberDiagonalFailures;
3833 }
3834 else
3835 {
3836 Real thisError = swaptionVegasMatrix[i][j];
3837 if (fabs(thisError) > 1e-8)
3838 ++offDiagonalFailures;
3839 }
3840 }
3841 }
3842
3843 if (numberDiagonalFailures + offDiagonalFailures>0 )
3844 BOOST_FAIL("Pathwise market vega analytic test fails for swaptions : " << offDiagonalFailures <<" off diagonal failures \n "
3845 << " and " << numberDiagonalFailures << " on the diagonal." );
3846 }
3847 // now do the caps
3848
3849 Matrix capsVegasMatrix(caps.size(), theBumps[0].size());
3850
3851 for (Size i=0; i < caps.size(); ++i)
3852 {
3853 CapPseudoDerivative thisPseudoDerivative(marketModel,
3854 caps[i].strike_,
3855 caps[i].startIndex_,
3856 caps[i].endIndex_, initialNumeraireValue
3857 );
3858
3859
3860 for (Size j=0; j < theBumps[0].size(); ++j)
3861 {
3862 capsVegasMatrix[i][j] = 0;
3863
3864 for (Size k=0; k < steps; ++k)
3865 for (Size l=0; l < numberRates; ++l)
3866 for (Size m=0; m < factors; ++m)
3867 capsVegasMatrix[i][j] += theBumps[k][j][l][m]*thisPseudoDerivative.volatilityDerivative(k)[l][m];
3868 }
3869 }
3870
3871 Size numberDiagonalFailures = 0;
3872 Size offDiagonalFailures=0;
3873
3874 for (Size i=0; i < caps.size(); ++i)
3875 {
3876 for (Size j=0; j < theBumps[0].size(); ++j)
3877 {
3878 if (i +swaptions.size()== j)
3879 {
3880 Real thisError = capsVegasMatrix[i][j] - 0.01;
3881
3882 if (fabs(thisError) > 1e-8)
3883 ++numberDiagonalFailures;
3884 }
3885 else
3886 {
3887 Real thisError = capsVegasMatrix[i][j];
3888 if (fabs(thisError) > 1e-8)
3889 ++offDiagonalFailures;
3890 }
3891 }
3892 }
3893
3894 if (numberDiagonalFailures + offDiagonalFailures>0 )
3895 BOOST_FAIL("Pathwise market vega analytic test fails for caps : " << offDiagonalFailures <<" off diagonal failures \n "
3896 << " and " << numberDiagonalFailures << " on the diagonal." );
3897
3898
3899 } // end of for (Size m=0; m<LENGTH(testedFactors); ++m)
3900 } // end of for (Size j=0; j<LENGTH(marketModels); j++)
3901 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3902 // test numerically differentiated swaptions against analytically done ones
3903 // we require equality on very path so we don't need many paths
3904
3905
3906 std::vector<Size> numberCashFlowsThisStep1(swaptionsDeflated.numberOfProducts());
3907
3908 std::vector<std::vector<MarketModelPathwiseMultiProduct::CashFlow> > cashFlowsGenerated1(swaptionsDeflated.numberOfProducts());
3909
3910
3911 for (Size i=0; i < swaptionsDeflated.numberOfProducts(); ++i)
3912 {
3913 cashFlowsGenerated1[i].resize(swaptionsDeflated.maxNumberOfCashFlowsPerProductPerStep());
3914 for (Size j=0; j < swaptionsDeflated.maxNumberOfCashFlowsPerProductPerStep(); ++j)
3915 cashFlowsGenerated1[i][j].amount.resize(numberRates+1);
3916 }
3917
3918 std::vector<Size> numberCashFlowsThisStep2(numberCashFlowsThisStep1);
3919 std::vector<std::vector<MarketModelPathwiseMultiProduct::CashFlow> >
3920 cashFlowsGenerated2(cashFlowsGenerated1);
3921
3922
3923
3924 for (Size j=0; j<LENGTH(marketModels); j++)
3925 {
3926
3927 Size testedFactors[] = { std::min<Size>(1UL,todaysForwards.size())
3928 // todaysForwards.size()
3929 //, 4, 8,
3930 };
3931
3932
3933
3934
3935 for (Size m=0; m<LENGTH(testedFactors); ++m)
3936 {
3937 Size factors = testedFactors[m];
3938
3939 MTBrownianGeneratorFactory generatorFactory(seed_);
3940
3941 bool logNormal = true;
3942
3943 ext::shared_ptr<MarketModel> marketModel =
3944 makeMarketModel(logNormal, evolution, factors,
3945 marketModels[j]);
3946
3947 LogNormalFwdRateEuler evolver1(marketModel,
3948 generatorFactory,swaptionsDeflated.suggestedNumeraires()
3949 );
3950
3951 LogNormalFwdRateEuler evolver2(marketModel,
3952 generatorFactory,swaptionsDeflated.suggestedNumeraires()
3953 );
3954
3955 for (Size p=0; p < pathsToDo; ++p)
3956 {
3957 evolver1.startNewPath();
3958 swaptionsDeflated.reset();
3959 evolver2.startNewPath();
3960 swaptionsDeflated2.reset();
3961 Size step =0;
3962
3963 bool done;
3964
3965 do
3966 {
3967 evolver1.advanceStep();
3968 done = swaptionsDeflated.nextTimeStep(evolver1.currentState(),
3969 numberCashFlowsThisStep1,
3970 cashFlowsGenerated1);
3971
3972 evolver2.advanceStep();
3973 bool done2 = swaptionsDeflated2.nextTimeStep(evolver2.currentState(),
3974 numberCashFlowsThisStep2,
3975 cashFlowsGenerated2);
3976
3977 if (done != done2)
3978 BOOST_FAIL("numerical swaptions derivative and swaptions disagree on termination");
3979
3980 for (Size prod = 0; prod < swaptionsDeflated.numberOfProducts(); ++prod)
3981 {
3982 if (numberCashFlowsThisStep1[prod] != numberCashFlowsThisStep2[prod])
3983 BOOST_FAIL("numerical swaptions derivative and swaptions disagree on number of cash flows");
3984
3985 for (Size cf =0; cf < numberCashFlowsThisStep1[prod]; ++cf)
3986 for (Size rate=0; rate<= numberRates; ++rate)
3987 if ( fabs(cashFlowsGenerated1[prod][cf].amount[rate] - cashFlowsGenerated2[prod][cf].amount[rate]) > tolerance )
3988 BOOST_FAIL("numerical swaptions derivative and swaptions disagree on cash flow size. cf = " << cf <<
3989 "step " << step << ", rate " << rate << ", amount1 " << cashFlowsGenerated1[prod][cf].amount[rate]
3990 << " ,amount2 " << cashFlowsGenerated2[prod][cf].amount[rate] << "\n");
3991
3992
3993
3994
3995
3996 }
3997
3998 ++step;
3999
4000
4001 }
4002 while (!done);
4003
4004
4005
4006 }
4007
4008
4009 } // end of for (Size m=0; m<LENGTH(testedFactors); ++m)
4010 } // end of for (Size j=0; j<LENGTH(marketModels); j++)
4011
4012 /////////////////////////////////////
4013
4014 // now time for the full simulation test
4015 // measure vega of each swaption with respect to itself, the other swaptions and the caps
4016 // should get 0.01 and 0 respectively.
4017 for (Size j=0; j<LENGTH(marketModels); j++)
4018 {
4019
4020 Size testedFactors[] = { std::min<Size>(1UL,todaysForwards.size())
4021 // todaysForwards.size()
4022 //, 4, 8,
4023 };
4024
4025
4026
4027
4028 for (Size m=0; m<LENGTH(testedFactors); ++m)
4029 {
4030 Size factors = testedFactors[m];
4031
4032 MTBrownianGeneratorFactory generatorFactory(seed_);
4033
4034 bool logNormal = true;
4035
4036 ext::shared_ptr<MarketModel> marketModel =
4037 makeMarketModel(logNormal, evolution, factors,
4038 marketModels[j]);
4039
4040 LogNormalFwdRateEuler evolver(marketModel,
4041 generatorFactory,swaptionsDeflated.suggestedNumeraires()
4042 );
4043
4044 Size initialNumeraire = evolver.numeraires().front();
4045 Real initialNumeraireValue =
4046 todaysDiscounts[initialNumeraire];
4047
4048
4049 // we need to work out our bumps
4050
4051 VegaBumpCollection possibleBumps(marketModel,
4052 allowFactorwiseBumping);
4053
4054
4055 OrthogonalizedBumpFinder bumpFinder(possibleBumps,
4056 swaptions,
4057 caps,
4058 multiplier, // if vector length grows by more than this discard
4059 tolerance); // if vector projection before scaling less than this discard
4060 std::vector<std::vector<Matrix> > theBumps;
4061
4062 bumpFinder.GetVegaBumps(theBumps);
4063
4064
4065 std::vector<Real> values;
4066
4067 std::vector<Real> errors;
4068
4069 {
4070
4071 PathwiseVegasAccountingEngine
4072 accountingEngine(ext::make_shared<LogNormalFwdRateEuler>(evolver),
4073 swaptionsDeflated,
4074 marketModel,
4075 theBumps,initialNumeraireValue);
4076
4077
4078 accountingEngine.multiplePathValues(values,errors,pathsToDoSimulation);
4079
4080 }
4081
4082 // we now have the simulation vegas, put them in more convenient form
4083
4084
4085 Matrix vegasMatrix(swaptionsDeflated.numberOfProducts(), theBumps[0].size());
4086 Matrix standardErrors(vegasMatrix);
4087 Size entriesPerProduct = 1+numberRates+theBumps[0].size();
4088
4089 for (Size i=0; i < swaptionsDeflated.numberOfProducts(); ++i)
4090 for (Size j=0; j < theBumps[0].size(); ++j)
4091 {
4092 vegasMatrix[i][j] = values[i*entriesPerProduct + numberRates+1+j];
4093 standardErrors[i][j] = errors[i*entriesPerProduct + numberRates+1 +j];
4094 }
4095
4096 // we next get the model vegas for comparison
4097
4098 std::vector<Real> impliedVols_(swaptions.size());
4099
4100 for (Size i=0; i < swaptions.size(); ++i)
4101 impliedVols_[i] = SwapForwardMappings::swaptionImpliedVolatility(*marketModel,
4102 swaptions[i].startIndex_,
4103 swaptions[i].endIndex_);
4104
4105 std::vector<Real> analyticVegas(swaptions.size());
4106 for (Size i=0; i < swaptions.size(); ++i)
4107 {
4108 Real swapRate = cs.coterminalSwapRates()[i];
4109 Real annuity = cs.coterminalSwapAnnuity(0,i)*initialNumeraireValue;
4110 Real expiry = rateTimes[i];
4111 Real sd = impliedVols_[i]*sqrt(expiry);
4112 Real swapDisplacement=0.0;
4113
4114 Real vega = blackFormulaVolDerivative(swapRate,
4115 swapRate,
4116 sd,
4117 expiry,
4118 annuity,
4119 swapDisplacement);
4120
4121 analyticVegas[i] = vega*0.01; // one percent move
4122
4123 }
4124
4125 // diagonal vegas should agree up to standard errors
4126 // off diagonal vegas should be zero
4127
4128 Size numberDiagonalFailures = 0;
4129 Size offDiagonalFailures=0;
4130
4131
4132 for (Size i=0; i < swaptions.size(); ++i)
4133 {
4134 Real thisError = vegasMatrix[i][i] - analyticVegas[i];
4135 Real thisErrorInSds = thisError / (standardErrors[i][i]+1E-6); // silly to penalize for tiny standard error
4136
4137 if (fabs(thisErrorInSds) > 4)
4138 ++numberDiagonalFailures;
4139
4140 }
4141
4142 for (Size i=0; i < swaptions.size(); ++i)
4143 for (Size j=0; j < theBumps[0].size(); ++j)
4144 {
4145 if ( i !=j )
4146 {
4147 Real thisError = vegasMatrix[i][j]; // true value is zero
4148
4149 Real thisErrorInSds = thisError / (standardErrors[i][j]+1E-6);
4150
4151 if (fabs(thisErrorInSds) > 3.5)
4152 ++offDiagonalFailures;
4153 }
4154 }
4155
4156 if (offDiagonalFailures + numberDiagonalFailures >0)
4157 BOOST_FAIL("Pathwise market vega test fails for coterminal swaptions : " << offDiagonalFailures <<" off diagonal failures \n "
4158 << " and " << numberDiagonalFailures << " on the diagonal." );
4159
4160
4161
4162
4163 } // end of for (Size m=0; m<LENGTH(testedFactors); ++m)
4164 } // end of for (Size j=0; j<LENGTH(marketModels); j++)
4165
4166 /////////////////////////////////////
4167 /////////////////////////////////////
4168
4169 // now time for the full simulation test
4170 // measure vega of each caps with respect to itself, the swaptions and the other caps
4171 // should get 0.01, 0 and 0 respectively.
4172 for (Size j=0; j<LENGTH(marketModels); j++)
4173 {
4174
4175 Size testedFactors[] = { std::min<Size>(2UL,todaysForwards.size())
4176 // todaysForwards.size()
4177 //, 4, 8,
4178 };
4179
4180
4181
4182
4183 for (Size m=0; m<LENGTH(testedFactors); ++m)
4184 {
4185 Size factors = testedFactors[m];
4186
4187 MTBrownianGeneratorFactory generatorFactory(seed_);
4188
4189 bool logNormal = true;
4190
4191 ext::shared_ptr<MarketModel> marketModel =
4192 makeMarketModel(logNormal, evolution, factors,
4193 marketModels[j]);
4194
4195 LogNormalFwdRateEuler evolver(marketModel,
4196 generatorFactory,capsDeflated.suggestedNumeraires()
4197 );
4198
4199 Size initialNumeraire = evolver.numeraires().front();
4200 Real initialNumeraireValue =
4201 todaysDiscounts[initialNumeraire];
4202
4203
4204 // we need to work out our bumps
4205
4206 VegaBumpCollection possibleBumps(marketModel,
4207 allowFactorwiseBumping);
4208
4209
4210 OrthogonalizedBumpFinder bumpFinder(possibleBumps,
4211 swaptions,
4212 caps,
4213 multiplier, // if vector length grows by more than this discard
4214 tolerance); // if vector projection before scaling less than this discard
4215 std::vector<std::vector<Matrix> > theBumps;
4216
4217 bumpFinder.GetVegaBumps(theBumps);
4218
4219
4220 std::vector<Real> values;
4221
4222 std::vector<Real> errors;
4223
4224 {
4225
4226 PathwiseVegasAccountingEngine
4227 accountingEngine(ext::make_shared<LogNormalFwdRateEuler>(evolver),
4228 capsDeflated,
4229 marketModel,
4230 theBumps,initialNumeraireValue);
4231
4232
4233 accountingEngine.multiplePathValues(values,errors,pathsToDoSimulation);
4234
4235 }
4236
4237 // we now have the simulation vegas, put them in more convenient form
4238
4239
4240 Matrix vegasMatrix(capsDeflated.numberOfProducts(), theBumps[0].size());
4241 Matrix standardErrors(vegasMatrix);
4242 Size entriesPerProduct = 1+numberRates+theBumps[0].size();
4243
4244
4245 for (Size i=0; i < capsDeflated.numberOfProducts(); ++i)
4246 for (Size j=0; j < theBumps[0].size(); ++j)
4247 {
4248 vegasMatrix[i][j] = values[i*entriesPerProduct +numberRates+j+1];
4249 standardErrors[i][j] = errors[i*entriesPerProduct +numberRates+j+1];
4250 }
4251
4252 // we next get the model vegas for comparison
4253
4254 std::vector<Real> impliedVols_(caps.size());
4255
4256
4257 std::vector<Real> analyticVegas(caps.size());
4258 for (Size i=0; i < caps.size(); ++i)
4259 {
4260
4261 CapPseudoDerivative capPseudo(marketModel,
4262 caps[i].strike_,
4263 caps[i].startIndex_,
4264 caps[i].endIndex_, initialNumeraireValue);
4265
4266 impliedVols_[i] = capPseudo.impliedVolatility();
4267
4268 Real vega=0.0;
4269
4270 for (Size j= caps[i].startIndex_; j< caps[i].endIndex_; ++j)
4271 {
4272
4273 Real forward = cs.forwardRates()[j];
4274 Real annuity = cs.discountRatio(j+1,0)*initialNumeraireValue*accruals[j];
4275 Real expiry = rateTimes[j];
4276 Real sd = impliedVols_[i]*sqrt(expiry);
4277 Real displacement=0.0;
4278
4279 Real capletVega = blackFormulaVolDerivative(caps[i].strike_,forward,
4280 sd,
4281 expiry,
4282 annuity,
4283 displacement);
4284
4285 vega += capletVega;
4286 }
4287
4288
4289
4290 analyticVegas[i] = vega*0.01; // one percent move
4291
4292 }
4293
4294 // diagonal vegas should agree up to standard errors
4295 // off diagonal vegas should be zero
4296
4297 Size numberDiagonalFailures = 0;
4298 Size offDiagonalFailures=0;
4299
4300
4301 for (Size i=0; i < caps.size(); ++i)
4302 {
4303 Real thisError = vegasMatrix[i][i+swaptions.size()] - analyticVegas[i];
4304 Real thisErrorInSds = thisError / (standardErrors[i][i+swaptions.size()]+1E-6); // silly to penalize for tiny standard error
4305
4306 if (fabs(thisErrorInSds) > 4)
4307 {
4308 BOOST_TEST_MESSAGE(" MC cap vega: " <<vegasMatrix[i][i+swaptions.size()] << " Analytic cap vega:" << analyticVegas[i] << " Error in sds:" << thisErrorInSds << "\n");
4309 ++numberDiagonalFailures;
4310 }
4311
4312 }
4313
4314 for (Size i=0; i < caps.size(); ++i)
4315 for (Size j=0; j < theBumps[0].size(); ++j)
4316 {
4317 if ( i+swaptions.size() !=j )
4318 {
4319 Real thisError = vegasMatrix[i][j]; // true value is zero
4320
4321 Real thisErrorInSds = thisError / (standardErrors[i][j]+1E-6);
4322
4323 if (fabs(thisErrorInSds) > 3.5)
4324 ++offDiagonalFailures;
4325 }
4326 }
4327
4328 if (offDiagonalFailures + numberDiagonalFailures >0)
4329 BOOST_FAIL("Pathwise market vega test fails for caps: " << offDiagonalFailures <<" off diagonal failures \n "
4330 << " and " << numberDiagonalFailures << " on the diagonal." );
4331
4332
4333
4334
4335 } // end of for (Size m=0; m<LENGTH(testedFactors); ++m)
4336 } // end of for (Size j=0; j<LENGTH(marketModels); j++)
4337
4338 /////////////////////////////////////
4339
4340
4341
4342
4343 }
4344
4345
4346
4347
4348 //--------------------- Volatility tests ---------------------
4349
testAbcdVolatilityIntegration()4350 void MarketModelTest::testAbcdVolatilityIntegration() {
4351
4352 BOOST_TEST_MESSAGE("Testing Abcd-volatility integration...");
4353
4354 using namespace market_model_test;
4355
4356 setup();
4357
4358 Real a = -0.0597;
4359 Real b = 0.1677;
4360 Real c = 0.5403;
4361 Real d = 0.1710;
4362
4363 const Size N = 10;
4364 const Real precision = 1e-04;
4365
4366 ext::shared_ptr<AbcdFunction> instVol(new AbcdFunction(a,b,c,d));
4367 SegmentIntegral SI(20000);
4368 for (Size i=0; i<N; i++) {
4369 Time T1 = 0.5*(1+i); // expiry of forward 1: after T1 AbcdVol = 0
4370 for (Size k=0; k<N-i; k++) {
4371 Time T2 = 0.5*(1+k); // expiry of forward 2: after T2 AbcdVol = 0
4372 //Integration
4373 for(Size j=0; j<N; j++) {
4374 Real xMin = 0.5*j;
4375 for (Size l=0; l<N-j; l++) {
4376 Real xMax = xMin + 0.5*l;
4377 AbcdSquared abcd2(a,b,c,d,T1,T2);
4378 Real numerical = SI(abcd2,xMin,xMax);
4379 Real analytical = instVol->covariance(xMin,xMax,T1,T2);
4380 if (std::abs(analytical-numerical)>precision) {
4381 BOOST_ERROR(" T1=" << T1 << "," <<
4382 "T2=" << T2 << ",\t\t" <<
4383 "xMin=" << xMin << "," <<
4384 "xMax=" << xMax << ",\t\t" <<
4385 "analytical: " << analytical << ",\t" <<
4386 "numerical: " << numerical);
4387 }
4388 if (T1==T2) {
4389 Real variance = instVol->variance(xMin,xMax,T1);
4390 if (std::abs(analytical-variance)>1e-14) {
4391 BOOST_ERROR(" T1=" << T1 << "," <<
4392 "T2=" << T2 << ",\t\t" <<
4393 "xMin=" << xMin << "," <<
4394 "xMax=" << xMax << ",\t\t" <<
4395 "variance: " << variance << ",\t" <<
4396 "analytical: " << analytical);
4397 }
4398 }
4399 }
4400 }
4401 }
4402 }
4403 }
4404
testAbcdVolatilityCompare()4405 void MarketModelTest::testAbcdVolatilityCompare() {
4406
4407 BOOST_TEST_MESSAGE("Testing different implementations of Abcd-volatility...");
4408
4409 using namespace market_model_test;
4410
4411 setup();
4412
4413 /*
4414 Given the instantaneous volatilities related to forward expiring at
4415 rateTimes[i1] and at rateTimes[i2], the methods:
4416 - LmExtLinearExponentialVolModel::integratedVariance(i1,i2,T)
4417 - Abcd::covariance(T)
4418 return the same result only if T < min(rateTimes[i1],rateTimes[i2]).
4419 */
4420
4421 // Parameters following Rebonato / Parameters following Brigo-Mercurio
4422 // used in Abcd class used in LmExtLinearExponentialVolModel
4423 Real a = 0.0597; // --> d
4424 Real b = 0.1677; // --> a
4425 Real c = 0.5403; // --> b
4426 Real d = 0.1710; // --> c
4427
4428 Size i1; // index of forward 1
4429 Size i2; // index of forward 2
4430
4431 ext::shared_ptr<LmVolatilityModel> lmAbcd(
4432 new LmExtLinearExponentialVolModel(rateTimes,b,c,d,a));
4433 ext::shared_ptr<AbcdFunction> abcd(new AbcdFunction(a,b,c,d));
4434 for (i1=0; i1<rateTimes.size(); i1++ ) {
4435 for (i2=0; i2<rateTimes.size(); i2++ ) {
4436 Time T = 0.;
4437 do {
4438 Real lmCovariance = lmAbcd->integratedVariance(i1,i2,T);
4439 Real abcdCovariance =
4440 abcd->covariance(0,T,rateTimes[i1],rateTimes[i2]);
4441 if(std::abs(lmCovariance-abcdCovariance)>1e-10) {
4442 BOOST_FAIL(" T1=" << rateTimes[i1] << "," <<
4443 "T2=" << rateTimes[i2] << ",\t\t" <<
4444 "xMin=" << 0 << "," <<
4445 "xMax=" << T << ",\t\t" <<
4446 "abcd: " << abcdCovariance << ",\t" <<
4447 "lm: " << lmCovariance);
4448 }
4449 T += 0.5;
4450 } while (T<std::min(rateTimes[i1],rateTimes[i2])) ;
4451 }
4452 }
4453 }
4454
testAbcdVolatilityFit()4455 void MarketModelTest::testAbcdVolatilityFit() {
4456
4457 BOOST_TEST_MESSAGE("Testing Abcd-volatility fit...");
4458
4459 using namespace market_model_test;
4460
4461 setup();
4462
4463 AbcdCalibration instVol(std::vector<Time>(rateTimes.begin(), rateTimes.end()-1), blackVols);
4464 Real a0 = instVol.a();
4465 Real b0 = instVol.b();
4466 Real c0 = instVol.c();
4467 Real d0 = instVol.d();
4468 Real error0 = instVol.error();
4469
4470 instVol.compute();
4471
4472 EndCriteria::Type ec = instVol.endCriteria();
4473 Real a1 = instVol.a();
4474 Real b1 = instVol.b();
4475 Real c1 = instVol.c();
4476 Real d1 = instVol.d();
4477 Real error1 = instVol.error();
4478
4479 if (error1>=error0)
4480 BOOST_FAIL("Parameters:" <<
4481 "\na: " << a0 << " ---> " << a1 <<
4482 "\nb: " << b0 << " ---> " << b1 <<
4483 "\nc: " << c0 << " ---> " << c1 <<
4484 "\nd: " << d0 << " ---> " << d1 <<
4485 "\nerror: " << error0 << " ---> " << error1);
4486
4487 AbcdFunction abcd(a1, b1, c1, d1);
4488 std::vector<Real> k = instVol.k(std::vector<Time>(rateTimes.begin(), rateTimes.end()-1), blackVols);
4489 Real tol = 3.0e-4;
4490 for (Size i=0; i<blackVols.size(); i++) {
4491 if (std::abs(k[i]-1.0)>tol) {
4492 Real modelVol =
4493 abcd.volatility(0.0, rateTimes[i], rateTimes[i]);
4494 BOOST_FAIL("\n EndCriteria = " << ec <<
4495 "\n Fixing Time = " << rateTimes[i] <<
4496 "\n MktVol = " << io::rate(blackVols[i]) <<
4497 "\n ModVol = " << io::rate(modelVol) <<
4498 "\n k = " << k[i] <<
4499 "\n error = " << std::abs(k[i]-1.0) <<
4500 "\n tol = " << tol);
4501 }
4502 }
4503
4504 }
4505
4506 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
testStochVolForwardsAndOptionlets()4507 void MarketModelTest::testStochVolForwardsAndOptionlets() {
4508
4509 BOOST_TEST_MESSAGE(
4510 "Testing exact repricing of "
4511 "forwards and optionlets "
4512 "in a stochastic vol displaced diffusion forward rate market model...");
4513
4514 using namespace market_model_test;
4515
4516 setup();
4517
4518 std::vector<Rate> forwardStrikes(todaysForwards.size());
4519 std::vector<ext::shared_ptr<Payoff> > optionletPayoffs(todaysForwards.size());
4520 /* std::vector<ext::shared_ptr<PlainVanillaPayoff> >
4521 displacedPayoffs(todaysForwards.size()); */
4522 for (Size i=0; i<todaysForwards.size(); ++i)
4523 {
4524 forwardStrikes[i] = todaysForwards[i] + 0.01;
4525 optionletPayoffs[i] = ext::shared_ptr<Payoff>(new
4526 PlainVanillaPayoff(Option::Call, todaysForwards[i]));
4527 /* displacedPayoffs[i] = ext::shared_ptr<PlainVanillaPayoff>(new
4528 PlainVanillaPayoff(Option::Call, todaysForwards[i]+displacement)); */
4529 }
4530
4531 MultiStepForwards forwards(rateTimes, accruals,
4532 paymentTimes, forwardStrikes);
4533 MultiStepOptionlets optionlets(rateTimes, accruals,
4534 paymentTimes, optionletPayoffs);
4535
4536 MultiProductComposite product;
4537 product.add(forwards);
4538 product.add(optionlets);
4539 product.finalize();
4540
4541 EvolutionDescription evolution = product.evolution();
4542
4543 MarketModelType marketModels[] =
4544 {
4545 ExponentialCorrelationFlatVolatility
4546 };
4547
4548 Size firstVolatilityFactor = 2;
4549 Size volatilityFactorStep = 2;
4550
4551 Real meanLevel=1.0;
4552 Real reversionSpeed=1.0;
4553
4554 Real volVar=1;
4555 Real v0=1.0;
4556 Size numberSubSteps=8;
4557 Real w1=0.5;
4558 Real w2=0.5;
4559 Real cutPoint = 1.5;
4560
4561 ext::shared_ptr<MarketModelVolProcess> volProcess(new
4562 SquareRootAndersen(meanLevel,
4563 reversionSpeed,
4564 volVar,
4565 v0,
4566 evolution.evolutionTimes(),
4567 numberSubSteps,
4568 w1,
4569 w2,
4570 cutPoint));
4571
4572 for (Size j=0; j<LENGTH(marketModels); j++)
4573 {
4574
4575 Size testedFactors[] = {1, 2, todaysForwards.size()};
4576 for (Size m=0; m<LENGTH(testedFactors); ++m) {
4577 Size factors = testedFactors[m];
4578
4579 MeasureType measures[] = { MoneyMarket, Terminal };
4580
4581 for (Size k=0; k<LENGTH(measures); k++)
4582 {
4583 std::vector<Size> numeraires = makeMeasure(product, measures[k]);
4584
4585 bool logNormal = true;
4586 ext::shared_ptr<MarketModel> marketModel =
4587 makeMarketModel(logNormal, evolution, factors, marketModels[j]);
4588
4589
4590 for (Size n=0; n<1; n++)
4591 {
4592 MTBrownianGeneratorFactory generatorFactory(seed_);
4593
4594 ext::shared_ptr<MarketModelEvolver> evolver(new SVDDFwdRatePc(marketModel,
4595 generatorFactory,
4596 volProcess,
4597 firstVolatilityFactor,
4598 volatilityFactorStep,
4599 numeraires
4600 ));
4601
4602
4603 std::ostringstream config;
4604 config <<
4605 marketModelTypeToString(marketModels[j]) << ", " <<
4606 factors << (factors>1 ? (factors==todaysForwards.size() ? " (full) factors, " : " factors, ") : " factor,") <<
4607 measureTypeToString(measures[k]) << ", " <<
4608 "SVDDFwdRatePc" << ", " <<
4609 "MT BGF";
4610 if (printReport_)
4611 BOOST_TEST_MESSAGE(" " << config.str());
4612
4613 ext::shared_ptr<SequenceStatisticsInc> stats =
4614 simulate(evolver, product);
4615
4616 std::vector<Real> results = stats->mean();
4617 std::vector<Real> errors = stats->errorEstimate();
4618
4619
4620 // check forwards
4621
4622
4623 for (Size i=0; i < accruals.size(); ++i)
4624 {
4625 Real trueValue = todaysDiscounts[i]- todaysDiscounts[i+1]*(1+ forwardStrikes[i]*accruals[i]);
4626 Real error = results[i] - trueValue;
4627 Real errorSds = error/ errors[i];
4628
4629 if (fabs(errorSds) > 3.5)
4630 BOOST_FAIL("error in sds: " << errorSds << " for forward " << i << " in SV LMM test. True value:" << trueValue << ", actual value: " << results[i] << " , standard error " << errors[i]);
4631
4632
4633
4634 }
4635
4636 for (Size i=0; i < accruals.size(); ++i)
4637 {
4638
4639 Real volCoeff = volatilities[i];
4640 // sqrt(marketModel->totalCovariance(i)[i][i]/evolution.evolutionTimes()[i]);
4641 Real theta = volCoeff*volCoeff*meanLevel;
4642 Real kappa = reversionSpeed;
4643 Real sigma = volCoeff*volVar;
4644 Real rho = 0.0;
4645 Real v1 = v0*volCoeff*volCoeff;
4646
4647
4648
4649
4650 ext::shared_ptr<StrikedTypePayoff> payoff(
4651 new PlainVanillaPayoff(Option::Call, forwardStrikes[i]));
4652
4653
4654
4655 Real trueValue =0.0;
4656 Size evaluations =0;
4657
4658 AnalyticHestonEngine::doCalculation(1.0, // no discounting
4659 1.0 ,// no discounting
4660 todaysForwards[i]+displacement,
4661 todaysForwards[i]+displacement,
4662 rateTimes[i],
4663 kappa,
4664 theta,
4665 sigma,
4666 v1,
4667 rho,
4668 *payoff,
4669 AnalyticHestonEngine::Integration::gaussLaguerre(),
4670 // AnalyticHestonEngine::Integration::gaussLobatto(1e-8, 1e-8),
4671 AnalyticHestonEngine::Gatheral,
4672 0,
4673 trueValue,
4674 evaluations);
4675
4676
4677 trueValue *= accruals[i]*todaysDiscounts[i+1];
4678
4679 // trueValue =
4680 // BlackCalculator(displacedPayoffs[i],
4681 // todaysForwards[i]+displacement,
4682 // volatilities[i]*std::sqrt(rateTimes[i]),
4683 // todaysDiscounts[i+1]*accruals[i]).value();
4684
4685
4686 Real error = results[i+ accruals.size()] - trueValue;
4687 Real errorSds = error/ errors[i];
4688
4689 if (fabs(errorSds) > 4)
4690 BOOST_FAIL("error in sds: " << errorSds << " for caplet " << i << " in SV LMM test. True value:" << trueValue << ", actual value: " << results[i+ accruals.size()] << " , standard error " << errors[i]);
4691
4692
4693
4694
4695 }
4696
4697
4698
4699
4700
4701
4702 }
4703 }
4704 }
4705 }
4706 }
4707
4708
4709
4710 //--------------------- Other tests ---------------------
4711
testDriftCalculator()4712 void MarketModelTest::testDriftCalculator() {
4713
4714 // Test full factor drift equivalence between compute() and
4715 // computeReduced()
4716
4717 BOOST_TEST_MESSAGE("Testing drift calculation...");
4718
4719 using namespace market_model_test;
4720
4721 setup();
4722
4723 Real tolerance = 1.0e-16;
4724 Size factors = todaysForwards.size();
4725 std::vector<Time> evolutionTimes(rateTimes.size()-1);
4726 std::copy(rateTimes.begin(), rateTimes.end()-1, evolutionTimes.begin());
4727 EvolutionDescription evolution(rateTimes,evolutionTimes);
4728 const std::vector<Real>& rateTaus = evolution.rateTaus();
4729 std::vector<Size> numeraires = moneyMarketPlusMeasure(evolution,
4730 measureOffset_);
4731 std::vector<Size> alive = evolution.firstAliveRate();
4732 Size numberOfSteps = evolutionTimes.size();
4733 std::vector<Real> drifts(numberOfSteps), driftsReduced(numberOfSteps);
4734 MarketModelType marketModels[] = {ExponentialCorrelationFlatVolatility,
4735 ExponentialCorrelationAbcdVolatility};
4736 for (Size k=0; k<LENGTH(marketModels); ++k) { // loop over market models
4737 bool logNormal = true;
4738 ext::shared_ptr<MarketModel> marketModel =
4739 makeMarketModel(logNormal, evolution, factors, marketModels[k]);
4740 std::vector<Rate> displacements = marketModel->displacements();
4741 for (Size j=0; j<numberOfSteps; ++j) { // loop over steps
4742 const Matrix& A = marketModel->pseudoRoot(j);
4743 //BOOST_TEST_MESSAGE(io::ordinal(j+1) << " pseudoroot:\n" << A);
4744 Size inf = std::max(0,static_cast<Integer>(alive[j]));
4745 for (Size h=inf; h<numeraires.size(); ++h) { // loop over numeraires
4746 LMMDriftCalculator driftcalculator(A, displacements, rateTaus,
4747 numeraires[h], alive[j]);
4748 driftcalculator.computePlain(todaysForwards, drifts);
4749 driftcalculator.computeReduced(todaysForwards,
4750 driftsReduced);
4751 for (Size i=0; i<drifts.size(); ++i) {
4752 Real error = std::abs(driftsReduced[i]-drifts[i]);
4753 if (error>tolerance)
4754 BOOST_ERROR("MarketModel: " <<
4755 marketModelTypeToString(marketModels[k]) <<
4756 ", " << io::ordinal(j+1) << " step, " <<
4757 ", " << io::ordinal(h+1) << " numeraire, " <<
4758 ", " << io::ordinal(i+1) << " drift, " <<
4759 "\ndrift =" << drifts[i] <<
4760 "\ndriftReduced =" << driftsReduced[i] <<
4761 "\n error =" << error <<
4762 "\n tolerance =" << tolerance);
4763 }
4764 }
4765 }
4766 }
4767 }
4768
testIsInSubset()4769 void MarketModelTest::testIsInSubset() {
4770
4771 // Performance test for isInSubset function (temporary)
4772
4773 BOOST_TEST_MESSAGE("Testing isInSubset function...");
4774
4775 using namespace market_model_test;
4776
4777 setup();
4778
4779 Size dim = 100;
4780 std::vector<Time> set, subset;
4781 for (Size i=0; i<dim; i++) set.push_back(i*1.0);
4782 for (Size i=0; i<dim; i++) subset.push_back(dim+i*1.0);
4783 std::valarray<bool> result = isInSubset(set, subset);
4784 if (printReport_) {
4785 for (Size i=0; i<dim; i++) {
4786 BOOST_TEST_MESSAGE(io::ordinal(i+1) << ":" <<
4787 " set[" << i << "] = " << set[i] <<
4788 ", subset[" << i << "] = " << subset[i] <<
4789 ", result[" << i << "] = " << result[i]);
4790 }
4791 }
4792 }
4793
4794
testAbcdDegenerateCases()4795 void MarketModelTest::testAbcdDegenerateCases() {
4796 BOOST_TEST_MESSAGE("Testing abcd degenerate cases...");
4797
4798 AbcdFunction f1(0.0,0.0,1.0E-15,1.0);
4799 AbcdFunction f2(1.0,0.0,1.0E-50,0.0);
4800
4801 Real cov1 = f1.covariance(0.0,1.0,1.0,1.0);
4802 if (std::fabs(cov1 - 1.0) > 1E-14
4803 || boost::math::isnan(cov1) || boost::math::isinf(cov1))
4804 BOOST_FAIL("(a,b,c,d)=(0,0,0,1): true covariance should be 1.0, "
4805 << "error is " << std::fabs(cov1 - 1.0));
4806
4807 Real cov2 = f2.covariance(0.0,1.0,1.0,1.0);
4808 if (std::fabs(cov2 - 1.0) > 1E-14
4809 || boost::math::isnan(cov2) || boost::math::isinf(cov2))
4810 BOOST_FAIL("(a,b,c,d)=(1,0,0,0): true covariance should be 1.0, "
4811 << "error is " << std::fabs(cov2 - 1.0));
4812 }
4813
testCovariance()4814 void MarketModelTest::testCovariance() {
4815 BOOST_TEST_MESSAGE("Testing market models covariance...");
4816
4817 const Size n = 10;
4818
4819 std::vector<Real> rateTimes;
4820 std::vector<Real> evolTimes1;
4821 std::vector<Real> evolTimes2;
4822 std::vector<Real> evolTimes3;
4823 std::vector<Real> evolTimes4;
4824 std::vector<std::vector<Real> > evolTimes;
4825
4826 for(Size i=1;i<=n;i++) rateTimes.push_back(static_cast<Time>(i));
4827 evolTimes1.push_back(n-1);
4828 for(Size i=1;i<=n-1;i++) evolTimes2.push_back(static_cast<Time>(i));
4829 for(Size i=1;i<=2*n-2;i++) evolTimes3.push_back(0.5*i);
4830 evolTimes4.push_back(0.3);
4831 evolTimes4.push_back(1.3);
4832 evolTimes4.push_back(2.0);
4833 evolTimes4.push_back(4.5);
4834 evolTimes4.push_back(8.2);
4835
4836 evolTimes.push_back(evolTimes1);
4837 evolTimes.push_back(evolTimes2);
4838 evolTimes.push_back(evolTimes3);
4839 evolTimes.push_back(evolTimes4);
4840
4841 std::vector<std::string> evolNames;
4842 evolNames.push_back("one evolution time");
4843 evolNames.push_back("evolution times on rate fixings");
4844 evolNames.push_back("evolution times on rate fixings and midpoints between fixings");
4845 evolNames.push_back("irregular evolution times");
4846
4847 std::vector<Real> ks(n-1,1.0);
4848 std::vector<Real> displ(n-1,0.0);
4849 std::vector<Real> rates(n-1,0.0);
4850 std::vector<Real> vols(n-1,1.0);
4851
4852 Matrix c = exponentialCorrelations(rateTimes,0.5,0.2,1.0,0.0);
4853 ext::shared_ptr<PiecewiseConstantCorrelation> corr(
4854 new TimeHomogeneousForwardCorrelation(c,rateTimes));
4855
4856 std::vector<std::string> modelNames;
4857 modelNames.push_back("FlatVol");
4858 modelNames.push_back("AbcdVol");
4859
4860 for(Size k=0;k<modelNames.size();k++) {
4861 for(Size l=0;l<evolNames.size();l++) {
4862 EvolutionDescription evolution(rateTimes,evolTimes[l]);
4863 ext::shared_ptr<MarketModel> model;
4864 switch(k) {
4865 case 0:
4866 model = ext::shared_ptr<MarketModel>(
4867 new FlatVol(vols,corr,evolution,n-1,rates,displ));
4868 break;
4869 case 1:
4870 model = ext::shared_ptr<MarketModel>(
4871 new AbcdVol(1.0,0.0,1.0E-50,0.0,ks,
4872 corr,evolution,n-1,rates,displ));
4873 break;
4874 default:
4875 BOOST_FAIL("Unknown model " << modelNames[k]);
4876 }
4877 if (model != 0) {
4878 for(Size i=0;i<evolTimes[l].size();i++) {
4879 Matrix cov = model->covariance(i);
4880 Real dt = evolTimes[l][i] - (i>0 ? evolTimes[l][i-1] : 0.0);
4881 for(Size x=0;x<n-1;x++) {
4882 for(Size y=0;y<n-1;y++) {
4883 if(std::min(rateTimes[x],rateTimes[y])>=evolTimes[l][i]
4884 && fabs(cov[x][y]-c[x][y]*dt)>1.0E-14)
4885 BOOST_FAIL("Model " << modelNames[k]
4886 << " with " << evolNames[l]
4887 << ": covariance matrix in step " << i
4888 << ": true value at (" << x << "," << y
4889 << ") is " << c[x][y]*dt
4890 << " actual value is " << cov[x][y]);
4891 }
4892 }
4893 }
4894 }
4895 }
4896 }
4897 }
4898
4899 // --- Call the desired tests
suite(SpeedLevel speed)4900 test_suite* MarketModelTest::suite(SpeedLevel speed) {
4901 test_suite* suite = BOOST_TEST_SUITE("Market-model tests");
4902
4903 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testInverseFloater));
4904
4905 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testPathwiseMarketVegas));
4906 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testPathwiseGreeks));
4907
4908 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testStochVolForwardsAndOptionlets));
4909
4910 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testOneStepForwardsAndOptionlets));
4911 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testOneStepNormalForwardsAndOptionlets));
4912
4913 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testGreeks));
4914
4915 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testAbcdVolatilityIntegration));
4916 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testAbcdVolatilityCompare));
4917 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testAbcdVolatilityFit));
4918
4919 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testPeriodAdapter));
4920
4921 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testDriftCalculator));
4922 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testIsInSubset));
4923
4924 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testAbcdDegenerateCases));
4925 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testCovariance));
4926
4927 if (speed <= Fast) {
4928 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testPathwiseVegas));
4929
4930 using namespace market_model_test;
4931
4932 setup();
4933
4934 MarketModelType marketModels[] = {
4935 ExponentialCorrelationFlatVolatility,
4936 ExponentialCorrelationAbcdVolatility
4937 };
4938
4939 Size testedFactors[] = { 4, 8, todaysForwards.size() };
4940 #define BOOST_PP_LOCAL_MACRO(n) \
4941 suite->add(QUANTLIB_TEST_CASE( \
4942 ext::bind(&MarketModelTest::testCallableSwapAnderson, \
4943 marketModels[n/LENGTH(testedFactors)], \
4944 testedFactors[n%LENGTH(testedFactors)])));
4945
4946 #define BOOST_PP_LOCAL_LIMITS (0, 5)
4947 #include BOOST_PP_LOCAL_ITERATE()
4948 }
4949
4950 if (speed == Slow) {
4951 suite->add(QUANTLIB_TEST_CASE(
4952 &MarketModelTest::testAllMultiStepProducts));
4953 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testCallableSwapNaif));
4954 suite->add(QUANTLIB_TEST_CASE(&MarketModelTest::testCallableSwapLS));
4955 }
4956
4957 return suite;
4958 }
4959