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