Logo
Fully annotated reference manual - version 1.8.12
Loading...
Searching...
No Matches
Functions | Variables
crossassetmodel.cpp File Reference
#include "toplevelfixture.hpp"
#include "utilities.hpp"
#include <boost/test/unit_test.hpp>
#include <boost/test/data/test_case.hpp>
#include <qle/methods/multipathgeneratorbase.hpp>
#include <qle/models/cdsoptionhelper.hpp>
#include <qle/models/cirppconstantfellerparametrization.hpp>
#include <qle/models/commodityschwartzmodel.hpp>
#include <qle/models/commodityschwartzparametrization.hpp>
#include <qle/models/cpicapfloorhelper.hpp>
#include <qle/models/crlgm1fparametrization.hpp>
#include <qle/models/crossassetanalytics.hpp>
#include <qle/models/crossassetanalyticsbase.hpp>
#include <qle/models/crossassetmodel.hpp>
#include <qle/models/crossassetmodelimpliedeqvoltermstructure.hpp>
#include <qle/models/crossassetmodelimpliedfxvoltermstructure.hpp>
#include <qle/models/dkimpliedyoyinflationtermstructure.hpp>
#include <qle/models/dkimpliedzeroinflationtermstructure.hpp>
#include <qle/models/eqbsconstantparametrization.hpp>
#include <qle/models/eqbsparametrization.hpp>
#include <qle/models/eqbspiecewiseconstantparametrization.hpp>
#include <qle/models/fxbsconstantparametrization.hpp>
#include <qle/models/fxbsparametrization.hpp>
#include <qle/models/fxbspiecewiseconstantparametrization.hpp>
#include <qle/models/fxeqoptionhelper.hpp>
#include <qle/models/gaussian1dcrossassetadaptor.hpp>
#include <qle/models/infdkparametrization.hpp>
#include <qle/models/irlgm1fconstantparametrization.hpp>
#include <qle/models/irlgm1fparametrization.hpp>
#include <qle/models/irlgm1fpiecewiseconstanthullwhiteadaptor.hpp>
#include <qle/models/irlgm1fpiecewiseconstantparametrization.hpp>
#include <qle/models/irlgm1fpiecewiselinearparametrization.hpp>
#include <qle/models/jyimpliedzeroinflationtermstructure.hpp>
#include <qle/models/lgm.hpp>
#include <qle/models/lgmimplieddefaulttermstructure.hpp>
#include <qle/models/lgmimpliedyieldtermstructure.hpp>
#include <qle/models/linkablecalibratedmodel.hpp>
#include <qle/models/parametrization.hpp>
#include <qle/models/piecewiseconstanthelper.hpp>
#include <qle/models/pseudoparameter.hpp>
#include <qle/pricingengines/analyticcclgmfxoptionengine.hpp>
#include <qle/pricingengines/analyticdkcpicapfloorengine.hpp>
#include <qle/pricingengines/analyticlgmcdsoptionengine.hpp>
#include <qle/pricingengines/analyticlgmswaptionengine.hpp>
#include <qle/pricingengines/analyticxassetlgmeqoptionengine.hpp>
#include <qle/pricingengines/blackcdsoptionengine.hpp>
#include <qle/pricingengines/crossccyswapengine.hpp>
#include <qle/pricingengines/depositengine.hpp>
#include <qle/pricingengines/discountingcommodityforwardengine.hpp>
#include <qle/pricingengines/discountingcurrencyswapengine.hpp>
#include <qle/pricingengines/discountingequityforwardengine.hpp>
#include <qle/pricingengines/discountingfxforwardengine.hpp>
#include <qle/pricingengines/discountingriskybondengine.hpp>
#include <qle/pricingengines/discountingswapenginemulticurve.hpp>
#include <qle/pricingengines/numericlgmmultilegoptionengine.hpp>
#include <qle/pricingengines/oiccbasisswapengine.hpp>
#include <qle/pricingengines/paymentdiscountingengine.hpp>
#include <qle/processes/commodityschwartzstateprocess.hpp>
#include <qle/processes/crossassetstateprocess.hpp>
#include <qle/processes/irlgm1fstateprocess.hpp>
#include <qle/termstructures/pricecurve.hpp>
#include <ql/currencies/america.hpp>
#include <ql/currencies/europe.hpp>
#include <ql/indexes/ibor/euribor.hpp>
#include <ql/indexes/ibor/gbplibor.hpp>
#include <ql/indexes/ibor/usdlibor.hpp>
#include <ql/indexes/inflation/euhicp.hpp>
#include <ql/indexes/inflation/ukrpi.hpp>
#include <ql/instruments/cpicapfloor.hpp>
#include <ql/instruments/vanillaoption.hpp>
#include <ql/math/optimization/levenbergmarquardt.hpp>
#include <ql/math/randomnumbers/rngtraits.hpp>
#include <ql/math/statistics/incrementalstatistics.hpp>
#include <ql/methods/montecarlo/multipathgenerator.hpp>
#include <ql/methods/montecarlo/pathgenerator.hpp>
#include <ql/models/shortrate/calibrationhelpers/swaptionhelper.hpp>
#include <ql/models/shortrate/onefactormodels/gsr.hpp>
#include <ql/pricingengines/swaption/gaussian1dswaptionengine.hpp>
#include <ql/pricingengines/credit/midpointcdsengine.hpp>
#include <ql/quotes/simplequote.hpp>
#include <ql/termstructures/credit/flathazardrate.hpp>
#include <ql/termstructures/inflation/interpolatedzeroinflationcurve.hpp>
#include <ql/termstructures/inflationtermstructure.hpp>
#include <ql/termstructures/yield/flatforward.hpp>
#include <ql/time/calendars/target.hpp>
#include <ql/time/daycounters/actual360.hpp>
#include <ql/time/daycounters/thirty360.hpp>
#include <boost/make_shared.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/covariance.hpp>
#include <boost/accumulators/statistics/error_of_mean.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/variates/covariate.hpp>

Go to the source code of this file.

Functions

 BOOST_AUTO_TEST_CASE (testBermudanLgm1fGsr)
 
 BOOST_AUTO_TEST_CASE (testBermudanLgmInvariances)
 
 BOOST_AUTO_TEST_CASE (testNonstandardBermudanSwaption)
 
 BOOST_AUTO_TEST_CASE (testLgm1fCalibration)
 
 BOOST_AUTO_TEST_CASE (testCcyLgm3fForeignPayouts)
 
 BOOST_AUTO_TEST_CASE (testLgm5fFxCalibration)
 
 BOOST_AUTO_TEST_CASE (testLgm5fFullCalibration)
 
 BOOST_AUTO_TEST_CASE (testLgm5fMoments)
 
 BOOST_AUTO_TEST_CASE (testLgmGsrEquivalence)
 
 BOOST_AUTO_TEST_CASE (testLgmMcWithShift)
 
 BOOST_AUTO_TEST_CASE (testIrFxCrCirppMartingaleProperty)
 
 BOOST_AUTO_TEST_CASE (testIrFxCrMoments)
 
 BOOST_AUTO_TEST_CASE (testIrFxCrCorrelationRecovery)
 
 BOOST_DATA_TEST_CASE (testIrFxInfCrComMartingaleProperty, bdata::make(infEurFlags) *bdata::make(infGbpFlags) *bdata::make(flatVolsFlags) *bdata::make(driftFreeState), infEurIsDk, infGbpIsDk, flatVols, driftFreeState)
 
 BOOST_DATA_TEST_CASE (testIrFxInfCrComMoments, bdata::make(infEurFlags) *bdata::make(infGbpFlags) *bdata::make(flatVolsFlags) *bdata::make(driftFreeState), infEurIsDk, infGbpIsDk, flatVols, driftFreeState)
 
 BOOST_AUTO_TEST_CASE (testIrFxInfCrEqMartingaleProperty)
 
 BOOST_AUTO_TEST_CASE (testIrFxInfCrEqMoments)
 
 BOOST_AUTO_TEST_CASE (testEqLgm5fPayouts)
 
 BOOST_AUTO_TEST_CASE (testEqLgm5fCalibration)
 
 BOOST_AUTO_TEST_CASE (testEqLgm5fMoments)
 
 BOOST_AUTO_TEST_CASE (testCorrelationRecovery)
 
 BOOST_AUTO_TEST_CASE (testIrFxInfCrCorrelationRecovery)
 
 BOOST_AUTO_TEST_CASE (testIrFxInfCrEqCorrelationRecovery)
 
 BOOST_AUTO_TEST_CASE (testCpiCalibrationByAlpha)
 
 BOOST_AUTO_TEST_CASE (testCpiCalibrationByH)
 
 BOOST_AUTO_TEST_CASE (testCrCalibration)
 

Variables

vector< bool > infEurFlags { true, false }
 
vector< bool > infGbpFlags { true, false }
 
vector< bool > flatVolsFlags { true, false }
 
vector< bool > driftFreeState { true, false }
 

Function Documentation

◆ BOOST_AUTO_TEST_CASE() [1/24]

BOOST_AUTO_TEST_CASE ( testBermudanLgm1fGsr  )

Definition at line 186 of file crossassetmodel.cpp.

186 {
187
188 BOOST_TEST_MESSAGE("Testing consistency of Bermudan swaption pricing in "
189 "LGM 1F and GSR models...");
190
191 BermudanTestData d;
192
193 // we use the Hull White adaptor for the LGM parametrization
194 // which should lead to equal Bermudan swaption prices
195 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm_p = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantHullWhiteAdaptor>(
196 EURCurrency(), d.yts, d.stepTimes_a, d.sigmas_a, d.stepTimes_a, d.kappas_a);
197
198 // fix any T forward measure
199 QuantLib::ext::shared_ptr<Gsr> gsr = QuantLib::ext::make_shared<Gsr>(d.yts, d.stepDates, d.sigmas, d.reversion, 50.0);
200
201 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> lgm = QuantLib::ext::make_shared<LinearGaussMarkovModel>(lgm_p);
202
203 QuantLib::ext::shared_ptr<Gaussian1dModel> lgm_g1d = QuantLib::ext::make_shared<Gaussian1dCrossAssetAdaptor>(lgm);
204
205 QuantLib::ext::shared_ptr<PricingEngine> swaptionEngineGsr =
206 QuantLib::ext::make_shared<Gaussian1dSwaptionEngine>(gsr, 64, 7.0, true, false);
207
208 QuantLib::ext::shared_ptr<PricingEngine> swaptionEngineLgm =
209 QuantLib::ext::make_shared<Gaussian1dSwaptionEngine>(lgm_g1d, 64, 7.0, true, false);
210
211 QuantLib::ext::shared_ptr<PricingEngine> swaptionEngineLgm2 =
212 QuantLib::ext::make_shared<NumericLgmSwaptionEngine>(lgm, 7.0, 16, 7.0, 32);
213
214 d.swaption->setPricingEngine(swaptionEngineGsr);
215 Real npvGsr = d.swaption->NPV();
216 d.swaption->setPricingEngine(swaptionEngineLgm);
217 Real npvLgm = d.swaption->NPV();
218 d.swaption->setPricingEngine(swaptionEngineLgm2);
219 Real npvLgm2 = d.swaption->NPV();
220
221 Real tol = 0.2E-4; // basis point tolerance
222
223 if (std::fabs(npvGsr - npvLgm) > tol)
224 BOOST_ERROR("Failed to verify consistency of Bermudan swaption price "
225 "in IrLgm1f / Gaussian1d adaptor engine ("
226 << npvLgm << ") and Gsr (" << npvGsr << ") models, tolerance is " << tol);
227
228 if (std::fabs(npvGsr - npvLgm2) > tol)
229 BOOST_ERROR("Failed to verify consistency of Bermudan swaption price "
230 "in IrLgm1f / Numeric LGM engine ("
231 << npvLgm2 << ") and Gsr (" << npvGsr << ") models, tolerance is " << tol);
232}

◆ BOOST_AUTO_TEST_CASE() [2/24]

BOOST_AUTO_TEST_CASE ( testBermudanLgmInvariances  )

Definition at line 234 of file crossassetmodel.cpp.

234 {
235
236 BOOST_TEST_MESSAGE("Testing LGM model invariances for Bermudan "
237 "swaption pricing...");
238
239 BermudanTestData d;
240
241 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm_p2 = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantHullWhiteAdaptor>(
242 EURCurrency(), d.yts, d.stepTimes_a, d.sigmas_a, d.stepTimes_a, d.kappas_a);
243
244 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> lgm2 = QuantLib::ext::make_shared<LinearGaussMarkovModel>(lgm_p2);
245
246 QuantLib::ext::shared_ptr<Gaussian1dModel> lgm_g1d2 = QuantLib::ext::make_shared<Gaussian1dCrossAssetAdaptor>(lgm2);
247
248 QuantLib::ext::shared_ptr<PricingEngine> swaptionEngineLgm2 =
249 QuantLib::ext::make_shared<Gaussian1dSwaptionEngine>(lgm_g1d2, 64, 7.0, true, false);
250
251 d.swaption->setPricingEngine(swaptionEngineLgm2);
252 Real npvLgm = d.swaption->NPV();
253
254 lgm_p2->shift() = -5.0;
255 lgm_p2->scaling() = 3.0;
256
257 // parametrizations are not observed, so we have to call update ourselves
258 lgm2->update();
259
260 Real npvLgm2 = d.swaption->NPV();
261
262 Real tol = 1.0E-5;
263
264 if (std::fabs(npvLgm - npvLgm2) > tol)
265 BOOST_ERROR("Failed to verify consistency of Bermudan swaption price "
266 "under LGM model invariances, difference is "
267 << (npvLgm - npvLgm2));
268
269} // testBermudanLgm1fGsr

◆ BOOST_AUTO_TEST_CASE() [3/24]

BOOST_AUTO_TEST_CASE ( testNonstandardBermudanSwaption  )

Definition at line 271 of file crossassetmodel.cpp.

271 {
272
273 BOOST_TEST_MESSAGE("Testing numeric LGM swaption engine for non-standard swaption...");
274
275 BermudanTestData d;
276
277 QuantLib::ext::shared_ptr<NonstandardSwaption> ns_swaption = QuantLib::ext::make_shared<NonstandardSwaption>(*d.swaption);
278
279 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm_p = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantHullWhiteAdaptor>(
280 EURCurrency(), d.yts, d.stepTimes_a, d.sigmas_a, d.stepTimes_a, d.kappas_a);
281
282 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> lgm = QuantLib::ext::make_shared<LinearGaussMarkovModel>(lgm_p);
283
284 QuantLib::ext::shared_ptr<PricingEngine> engine = QuantLib::ext::make_shared<NumericLgmSwaptionEngine>(lgm, 7.0, 16, 7.0, 32);
285 QuantLib::ext::shared_ptr<PricingEngine> ns_engine =
286 QuantLib::ext::make_shared<NumericLgmNonstandardSwaptionEngine>(lgm, 7.0, 16, 7.0, 32);
287
288 d.swaption->setPricingEngine(engine);
289 ns_swaption->setPricingEngine(ns_engine);
290
291 Real npv = d.swaption->NPV();
292 Real ns_npv = d.swaption->NPV();
293
294 Real tol = 1.0E-12;
295 if (std::fabs(npv - ns_npv) >= tol)
296 BOOST_ERROR("Failed to verify consistency of Bermudan swaption price ("
297 << npv << ") and Bermudan nonstandard swaption price (" << ns_npv << "), difference is "
298 << (npv - ns_npv) << ", tolerance is " << tol);
299} // testNonstandardBermudanSwaption

◆ BOOST_AUTO_TEST_CASE() [4/24]

BOOST_AUTO_TEST_CASE ( testLgm1fCalibration  )

Definition at line 301 of file crossassetmodel.cpp.

301 {
302
303 BOOST_TEST_MESSAGE("Testing calibration of LGM 1F model (analytic engine) "
304 "against GSR parameters...");
305
306 // for fixed kappa != 0.0 we calibrate sigma via
307 // the Hull White Adaptor
308
309 SavedSettings backup;
310
311 Date evalDate(12, January, 2015);
312 Settings::instance().evaluationDate() = evalDate;
313 Handle<YieldTermStructure> yts(QuantLib::ext::make_shared<FlatForward>(evalDate, 0.02, Actual365Fixed()));
314 QuantLib::ext::shared_ptr<IborIndex> euribor6m = QuantLib::ext::make_shared<Euribor>(6 * Months, yts);
315
316 // coterminal basket 1y-9y, 2y-8y, ... 9y-1y
317
318 std::vector<QuantLib::ext::shared_ptr<BlackCalibrationHelper> > basket;
319 Real impliedVols[] = { 0.4, 0.39, 0.38, 0.35, 0.35, 0.34, 0.33, 0.32, 0.31 };
320 std::vector<Date> expiryDates;
321
322 for (Size i = 0; i < 9; ++i) {
323 QuantLib::ext::shared_ptr<BlackCalibrationHelper> helper = QuantLib::ext::make_shared<SwaptionHelper>(
324 (i + 1) * Years, (9 - i) * Years, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(impliedVols[i])), euribor6m,
325 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(), yts);
326 basket.push_back(helper);
327 expiryDates.push_back(
328 QuantLib::ext::static_pointer_cast<SwaptionHelper>(helper)->swaption()->exercise()->dates().back());
329 }
330
331 std::vector<Date> stepDates(expiryDates.begin(), expiryDates.end() - 1);
332
333 Array stepTimes_a(stepDates.size());
334 for (Size i = 0; i < stepDates.size(); ++i) {
335 stepTimes_a[i] = yts->timeFromReference(stepDates[i]);
336 }
337
338 Real kappa = 0.05;
339
340 std::vector<Real> gsrInitialSigmas(stepDates.size() + 1, 0.0050);
341 std::vector<Real> lgmInitialSigmas2(stepDates.size() + 1, 0.0050);
342
343 Array lgmInitialSigmas2_a(lgmInitialSigmas2.begin(), lgmInitialSigmas2.end());
344 Array kappas_a(lgmInitialSigmas2_a.size(), kappa);
345
346 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm_p = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantHullWhiteAdaptor>(
347 EURCurrency(), yts, stepTimes_a, lgmInitialSigmas2_a, stepTimes_a, kappas_a);
348
349 // fix any T forward measure
350 QuantLib::ext::shared_ptr<Gsr> gsr = QuantLib::ext::make_shared<Gsr>(yts, stepDates, gsrInitialSigmas, kappa, 50.0);
351
352 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> lgm = QuantLib::ext::make_shared<LinearGaussMarkovModel>(lgm_p);
353
354 QuantLib::ext::shared_ptr<PricingEngine> swaptionEngineGsr =
355 QuantLib::ext::make_shared<Gaussian1dSwaptionEngine>(gsr, 64, 7.0, true, false);
356
357 QuantLib::ext::shared_ptr<PricingEngine> swaptionEngineLgm = QuantLib::ext::make_shared<AnalyticLgmSwaptionEngine>(lgm);
358
359 // calibrate GSR
360
361 LevenbergMarquardt lm(1E-8, 1E-8, 1E-8);
362 EndCriteria ec(1000, 500, 1E-8, 1E-8, 1E-8);
363
364 for (Size i = 0; i < basket.size(); ++i) {
365 basket[i]->setPricingEngine(swaptionEngineGsr);
366 }
367
368 gsr->calibrateVolatilitiesIterative(basket, lm, ec);
369
370 Array gsrSigmas = gsr->volatility();
371
372 // calibrate LGM
373
374 for (Size i = 0; i < basket.size(); ++i) {
375 basket[i]->setPricingEngine(swaptionEngineLgm);
376 }
377
378 lgm->calibrateVolatilitiesIterative(basket, lm, ec);
379
380 Array lgmSigmas = lgm->parametrization()->parameterValues(0);
381
382 Real tol0 = 1E-8;
383 Real tol = 2E-5;
384
385 for (Size i = 0; i < gsrSigmas.size(); ++i) {
386 // check calibration itself, we should match the market prices
387 // rather exactly (note that this tests the lgm calibration,
388 // not the gsr calibration)
389 if (std::fabs(basket[i]->modelValue() - basket[i]->marketValue()) > tol0)
390 BOOST_ERROR("Failed to calibrate to market swaption #"
391 << i << ", market price is " << basket[i]->marketValue() << " while model price is "
392 << basket[i]->modelValue());
393 // compare calibrated model parameters
394 if (std::fabs(gsrSigmas[i] - lgmSigmas[i]) > tol)
395 BOOST_ERROR("Failed to verify LGM's sigma from Hull White adaptor (#"
396 << i << "), which is " << lgmSigmas[i] << " while GSR's sigma is " << gsrSigmas[i] << ")");
397 }
398
399 // calibrate LGM as component of CrossAssetModel
400
401 // create a second set of parametrization ...
402 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm_p21 = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantHullWhiteAdaptor>(
403 USDCurrency(), yts, stepTimes_a, lgmInitialSigmas2_a, stepTimes_a, kappas_a);
404 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm_p22 = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantHullWhiteAdaptor>(
405 EURCurrency(), yts, stepTimes_a, lgmInitialSigmas2_a, stepTimes_a, kappas_a);
406
407 // ... and a fx parametrization ...
408 Array notimes_a(0);
409 Array sigma_a(1, 0.10);
410 QuantLib::ext::shared_ptr<FxBsParametrization> fx_p = QuantLib::ext::make_shared<FxBsPiecewiseConstantParametrization>(
411 EURCurrency(), Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(1.00)), notimes_a, sigma_a);
412
413 // ... and set up an crossasset model with USD as domestic currency ...
414 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
415 parametrizations.push_back(lgm_p21);
416 parametrizations.push_back(lgm_p22);
417 parametrizations.push_back(fx_p);
418 Matrix rho(3, 3, 0.0);
419 rho[0][0] = rho[1][1] = rho[2][2] = 1.0;
420 QuantLib::ext::shared_ptr<CrossAssetModel> xmodel =
421 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, rho, SalvagingAlgorithm::None);
422
423 // .. whose EUR component we calibrate as before and compare the
424 // result against the 1d case and as well check that the USD
425 // component was not touched by the calibration.
426
427 QuantLib::ext::shared_ptr<PricingEngine> swaptionEngineLgm2 = QuantLib::ext::make_shared<AnalyticLgmSwaptionEngine>(xmodel, 1);
428
429 for (Size i = 0; i < basket.size(); ++i) {
430 basket[i]->setPricingEngine(swaptionEngineLgm2);
431 }
432
433 xmodel->calibrateIrLgm1fVolatilitiesIterative(1, basket, lm, ec);
434
435 Array lgmSigmas2eur = xmodel->irlgm1f(1)->parameterValues(0);
436 Array lgmSigmas2usd = xmodel->irlgm1f(0)->parameterValues(0);
437
438 for (Size i = 0; i < gsrSigmas.size(); ++i) {
439 // compare calibrated model parameters against 1d calibration before
440 if (!close_enough(lgmSigmas2eur[i], lgmSigmas[i]))
441 BOOST_ERROR("Failed to verify crossasset LGM1F component calibration "
442 "at parameter #"
443 << i << " against 1d calibration, which is " << lgmSigmas2eur[i] << " while 1d calibration was "
444 << lgmSigmas[i] << ")");
445 // compare usd component against start values (since it was not
446 // calibrated, so should not have changed)
447 if (!close_enough(lgmSigmas2usd[i], lgmInitialSigmas2[i]))
448 BOOST_ERROR("Non calibrated crossasset LGM1F component was changed by "
449 "other's component calibration at #"
450 << i << ", the new value is " << lgmSigmas2usd[i] << " while the initial value was "
451 << lgmInitialSigmas2[i]);
452 }
453
454} // testLgm1fCalibration
QuantExt::close_enough
Filter close_enough(const RandomVariable &x, const RandomVariable &y)
Definition: randomvariable.cpp:838
helper
QuantLib::BootstrapHelper< QuantLib::OptionletVolatilityStructure > helper
Definition: piecewiseatmoptionletcurve.cpp:45
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [5/24]

BOOST_AUTO_TEST_CASE ( testCcyLgm3fForeignPayouts  )

Definition at line 456 of file crossassetmodel.cpp.

456 {
457
458 BOOST_TEST_MESSAGE("Testing pricing of foreign payouts under domestic "
459 "measure in Ccy LGM 3F model...");
460
461 SavedSettings backup;
462
463 Date referenceDate(30, July, 2015);
464
465 Settings::instance().evaluationDate() = referenceDate;
466
467 Handle<YieldTermStructure> eurYts(QuantLib::ext::make_shared<FlatForward>(referenceDate, 0.02, Actual365Fixed()));
468
469 Handle<YieldTermStructure> usdYts(QuantLib::ext::make_shared<FlatForward>(referenceDate, 0.05, Actual365Fixed()));
470
471 // use different grids for the EUR and USD models and the FX volatility
472 // process to test the piecewise numerical integration ...
473
474 std::vector<Date> volstepdatesEur, volstepdatesUsd, volstepdatesFx;
475
476 volstepdatesEur.push_back(Date(15, July, 2016));
477 volstepdatesEur.push_back(Date(15, July, 2017));
478 volstepdatesEur.push_back(Date(15, July, 2018));
479 volstepdatesEur.push_back(Date(15, July, 2019));
480 volstepdatesEur.push_back(Date(15, July, 2020));
481
482 volstepdatesUsd.push_back(Date(13, April, 2016));
483 volstepdatesUsd.push_back(Date(13, September, 2016));
484 volstepdatesUsd.push_back(Date(13, April, 2017));
485 volstepdatesUsd.push_back(Date(13, September, 2017));
486 volstepdatesUsd.push_back(Date(13, April, 2018));
487 volstepdatesUsd.push_back(Date(15, July, 2018)); // shared with EUR
488 volstepdatesUsd.push_back(Date(13, April, 2019));
489 volstepdatesUsd.push_back(Date(13, September, 2019));
490
491 volstepdatesFx.push_back(Date(15, July, 2016)); // shared with EUR
492 volstepdatesFx.push_back(Date(15, October, 2016));
493 volstepdatesFx.push_back(Date(15, May, 2017));
494 volstepdatesFx.push_back(Date(13, September, 2017)); // shared with USD
495 volstepdatesFx.push_back(Date(15, July, 2018)); // shared with EUR and USD
496
497 std::vector<Real> eurVols, usdVols, fxVols;
498
499 for (Size i = 0; i < volstepdatesEur.size() + 1; ++i) {
500 eurVols.push_back(0.0050 + (0.0080 - 0.0050) * std::exp(-0.3 * static_cast<double>(i)));
501 }
502 for (Size i = 0; i < volstepdatesUsd.size() + 1; ++i) {
503 usdVols.push_back(0.0030 + (0.0110 - 0.0030) * std::exp(-0.3 * static_cast<double>(i)));
504 }
505 for (Size i = 0; i < volstepdatesFx.size() + 1; ++i) {
506 fxVols.push_back(0.15 + (0.20 - 0.15) * std::exp(-0.3 * static_cast<double>(i)));
507 }
508
509 Array alphaTimesEur(volstepdatesEur.size()), alphaEur(eurVols.begin(), eurVols.end()), kappaTimesEur(0),
510 kappaEur(1, 0.02);
511 Array alphaTimesUsd(volstepdatesUsd.size()), alphaUsd(usdVols.begin(), usdVols.end()), kappaTimesUsd(0),
512 kappaUsd(1, 0.04);
513 Array fxTimes(volstepdatesFx.size()), fxSigmas(fxVols.begin(), fxVols.end());
514
515 for (Size i = 0; i < alphaTimesEur.size(); ++i) {
516 alphaTimesEur[i] = eurYts->timeFromReference(volstepdatesEur[i]);
517 }
518 for (Size i = 0; i < alphaTimesUsd.size(); ++i) {
519 alphaTimesUsd[i] = eurYts->timeFromReference(volstepdatesUsd[i]);
520 }
521 for (Size i = 0; i < fxTimes.size(); ++i) {
522 fxTimes[i] = eurYts->timeFromReference(volstepdatesFx[i]);
523 }
524
525 QuantLib::ext::shared_ptr<IrLgm1fParametrization> eurLgmParam = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantParametrization>(
526 EURCurrency(), eurYts, alphaTimesEur, alphaEur, kappaTimesEur, kappaEur);
527
528 QuantLib::ext::shared_ptr<IrLgm1fParametrization> usdLgmParam = QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantParametrization>(
529 USDCurrency(), usdYts, alphaTimesUsd, alphaUsd, kappaTimesUsd, kappaUsd);
530
531 // USD per EUR (foreign per domestic)
532 Handle<Quote> usdEurSpotToday(QuantLib::ext::make_shared<SimpleQuote>(0.90));
533
534 QuantLib::ext::shared_ptr<FxBsParametrization> fxUsdEurBsParam =
535 QuantLib::ext::make_shared<FxBsPiecewiseConstantParametrization>(USDCurrency(), usdEurSpotToday, fxTimes, fxSigmas);
536
537 std::vector<QuantLib::ext::shared_ptr<Parametrization> > singleModels;
538 singleModels.push_back(eurLgmParam);
539 singleModels.push_back(usdLgmParam);
540 singleModels.push_back(fxUsdEurBsParam);
541
542 QuantLib::ext::shared_ptr<CrossAssetModel> ccLgm = QuantLib::ext::make_shared<CrossAssetModel>(singleModels);
543
544 Size eurIdx = ccLgm->ccyIndex(EURCurrency());
545 Size usdIdx = ccLgm->ccyIndex(USDCurrency());
546 Size eurUsdIdx = usdIdx - 1;
547
548 ccLgm->setCorrelation(CrossAssetModel::AssetType::IR, eurIdx, CrossAssetModel::AssetType::IR, usdIdx, -0.2);
549 ccLgm->setCorrelation(CrossAssetModel::AssetType::IR, eurIdx, CrossAssetModel::AssetType::FX, eurUsdIdx, 0.8);
550 ccLgm->setCorrelation(CrossAssetModel::AssetType::IR, usdIdx, CrossAssetModel::AssetType::FX, eurUsdIdx, -0.5);
551
552 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> eurLgm = QuantLib::ext::make_shared<LinearGaussMarkovModel>(eurLgmParam);
553 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> usdLgm = QuantLib::ext::make_shared<LinearGaussMarkovModel>(usdLgmParam);
554
555 QuantLib::ext::shared_ptr<StochasticProcess> process = ccLgm->stateProcess();
556 QuantLib::ext::shared_ptr<StochasticProcess> usdProcess = usdLgm->stateProcess();
557
558 // path generation
559
560 Size n = 500000; // number of paths
561 Size seed = 121; // seed
562 // maturity of test payoffs
563 Time T = 5.0;
564 // take large steps, but not only one (since we are testing)
565 Size steps = static_cast<Size>(T * 2.0);
566 TimeGrid grid(T, steps);
567 PseudoRandom::rsg_type sg2 = PseudoRandom::make_sequence_generator(steps, seed);
568
569 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process)) {
570 tmp->resetCache(grid.size() - 1);
571 }
572 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(usdProcess)) {
573 tmp->resetCache(grid.size() - 1);
574 }
575 MultiPathGeneratorMersenneTwister pg(process, grid, seed, false);
576 PathGenerator<PseudoRandom::rsg_type> pg2(usdProcess, grid, sg2, false);
577
578 // test
579 // 1 deterministic USD cashflow under EUR numeraire vs. price on USD curve
580 // 2 zero bond option USD under EUR numeraire vs. USD numeraire
581 // 3 fx option USD-EUR under EUR numeraire vs. analytical price
582
583 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > stat1, stat2a, stat2b, stat3;
584
585 // same for paths2 since shared time grid
586 for (Size j = 0; j < n; ++j) {
587 Sample<MultiPath> path = pg.next();
588 Sample<Path> path2 = pg2.next();
589 Size l = path.value[0].length() - 1;
590 Real fx = std::exp(path.value[2][l]);
591 Real zeur = path.value[0][l];
592 Real zusd = path.value[1][l];
593 Real zusd2 = path2.value[l];
594
595 // 1 USD paid at T deflated with EUR numeraire
596 stat1(1.0 * fx / eurLgm->numeraire(T, zeur));
597
598 // 2 USD zero bond option at T on P(T,T+10) strike 0.5 ...
599 // ... under EUR numeraire ...
600 Real zbOpt = std::max(usdLgm->discountBond(T, T + 10.0, zusd) - 0.5, 0.0);
601 stat2a(zbOpt * fx / eurLgm->numeraire(T, zeur));
602 // ... and under USD numeraire ...
603 Real zbOpt2 = std::max(usdLgm->discountBond(T, T + 10.0, zusd2) - 0.5, 0.0);
604 stat2b(zbOpt2 / usdLgm->numeraire(T, zusd2));
605
606 // 3 USD-EUR fx option @0.9
607 stat3(std::max(fx - 0.9, 0.0) / eurLgm->numeraire(T, zeur));
608 }
609
610 QuantLib::ext::shared_ptr<VanillaOption> fxOption =
611 QuantLib::ext::make_shared<VanillaOption>(QuantLib::ext::make_shared<PlainVanillaPayoff>(Option::Call, 0.9),
612 QuantLib::ext::make_shared<EuropeanExercise>(referenceDate + 5 * 365));
613
614 QuantLib::ext::shared_ptr<AnalyticCcLgmFxOptionEngine> ccLgmFxOptionEngine =
615 QuantLib::ext::make_shared<AnalyticCcLgmFxOptionEngine>(ccLgm, 0);
616
617 ccLgmFxOptionEngine->cache();
618
619 fxOption->setPricingEngine(ccLgmFxOptionEngine);
620
621 Real npv1 = mean(stat1);
622 Real error1 = error_of<tag::mean>(stat1);
623 Real expected1 = usdYts->discount(5.0) * usdEurSpotToday->value();
624 Real npv2a = mean(stat2a);
625 Real error2a = error_of<tag::mean>(stat2a);
626 Real npv2b = mean(stat2b) * usdEurSpotToday->value();
627 ;
628 Real error2b = error_of<tag::mean>(stat2b) * usdEurSpotToday->value();
629 Real npv3 = mean(stat3);
630 Real error3 = error_of<tag::mean>(stat3);
631
632 // accept this relative difference in error estimates
633 Real tolError = 0.2;
634 // accept tolErrEst*errorEstimate as absolute difference
635 Real tolErrEst = 1.0;
636
637 if (std::fabs((error1 - 4E-4) / 4E-4) > tolError)
638 BOOST_ERROR("error estimate deterministic "
639 "cashflow pricing can not be "
640 "reproduced, is "
641 << error1 << ", expected 4E-4, relative tolerance " << tolError);
642 if (std::fabs((error2a - 1E-4) / 1E-4) > tolError)
643 BOOST_ERROR("error estimate zero bond "
644 "option pricing (foreign measure) can "
645 "not be reproduced, is "
646 << error2a << ", expected 1E-4, relative tolerance " << tolError);
647 if (std::fabs((error2b - 7E-5) / 7E-5) > tolError)
648 BOOST_ERROR("error estimate zero bond "
649 "option pricing (domestic measure) can "
650 "not be reproduced, is "
651 << error2b << ", expected 7E-5, relative tolerance " << tolError);
652 if (std::fabs((error3 - 2.7E-4) / 2.7E-4) > tolError)
653 BOOST_ERROR("error estimate fx option pricing can not be reproduced, is "
654 << error3 << ", expected 2.7E-4, relative tolerance " << tolError);
655
656 if (std::fabs(npv1 - expected1) > tolErrEst * error1)
657 BOOST_ERROR("can no reproduce deterministic cashflow pricing, is "
658 << npv1 << ", expected " << expected1 << ", tolerance " << tolErrEst << "*" << error1);
659
660 if (std::fabs(npv2a - npv2b) > tolErrEst * std::sqrt(error2a * error2a + error2b * error2b))
661 BOOST_ERROR("can no reproduce zero bond option pricing, domestic "
662 "measure result is "
663 << npv2a << ", foreign measure result is " << npv2b << ", tolerance " << tolErrEst << "*"
664 << std::sqrt(error2a * error2a + error2b * error2b));
665
666 if (std::fabs(npv3 - fxOption->NPV()) > tolErrEst * error3)
667 BOOST_ERROR("can no reproduce fx option pricing, monte carlo result is "
668 << npv3 << ", analytical pricing result is " << fxOption->NPV() << ", tolerance is " << tolErrEst
669 << "*" << error3);
670
671} // testCcyLgm3fForeignPayouts
QuantExt::MultiPathGeneratorMersenneTwister
Instantiation of MultiPathGenerator with standard PseudoRandom traits.
Definition: multipathgeneratorbase.hpp:59
steps
std::vector< Size > steps
Definition: commodityschwartzmodel.cpp:109
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [6/24]

BOOST_AUTO_TEST_CASE ( testLgm5fFxCalibration  )

Definition at line 967 of file crossassetmodel.cpp.

967 {
968
969 BOOST_TEST_MESSAGE("Testing fx calibration in Ccy LGM 5F model...");
970
971 Lgm5fTestData d;
972
973 // we test the 5f model against the 3f model eur-gbp
974 std::vector<QuantLib::ext::shared_ptr<Parametrization> > singleModelsProjected;
975 singleModelsProjected.push_back(d.eurLgm_p);
976 singleModelsProjected.push_back(d.gbpLgm_p);
977 singleModelsProjected.push_back(d.fxGbp_p);
978
979 Matrix cProjected(3, 3);
980 for (Size i = 0, ii = 0; i < 5; ++i) {
981 if (i != 0 && i != 3) {
982 for (Size j = 0, jj = 0; j < 5; ++j) {
983 if (j != 0 && j != 3)
984 cProjected[ii][jj++] = d.c[i][j];
985 }
986 ++ii;
987 }
988 }
989
990 QuantLib::ext::shared_ptr<CrossAssetModel> ccLgmProjected =
991 QuantLib::ext::make_shared<CrossAssetModel>(singleModelsProjected, cProjected, SalvagingAlgorithm::None);
992
993 QuantLib::ext::shared_ptr<AnalyticCcLgmFxOptionEngine> ccLgmFxOptionEngineUsd =
994 QuantLib::ext::make_shared<AnalyticCcLgmFxOptionEngine>(d.ccLgmExact, 0);
995
996 QuantLib::ext::shared_ptr<AnalyticCcLgmFxOptionEngine> ccLgmFxOptionEngineGbp =
997 QuantLib::ext::make_shared<AnalyticCcLgmFxOptionEngine>(d.ccLgmExact, 1);
998
999 QuantLib::ext::shared_ptr<AnalyticCcLgmFxOptionEngine> ccLgmProjectedFxOptionEngineGbp =
1000 QuantLib::ext::make_shared<AnalyticCcLgmFxOptionEngine>(ccLgmProjected, 0);
1001
1002 ccLgmFxOptionEngineUsd->cache();
1003 ccLgmFxOptionEngineGbp->cache();
1004 ccLgmProjectedFxOptionEngineGbp->cache();
1005
1006 // while the initial fx vol starts at 0.2 for usd and 0.15 for gbp
1007 // we calibrate to helpers with 0.15 and 0.2 target implied vol
1008 std::vector<QuantLib::ext::shared_ptr<BlackCalibrationHelper> > helpersUsd, helpersGbp;
1009 for (Size i = 0; i <= d.volstepdatesFx.size(); ++i) {
1010 QuantLib::ext::shared_ptr<BlackCalibrationHelper> tmpUsd = QuantLib::ext::make_shared<FxEqOptionHelper>(
1011 i < d.volstepdatesFx.size() ? d.volstepdatesFx[i] : d.volstepdatesFx.back() + 365, 0.90, d.fxEurUsd,
1012 Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.15)), d.ccLgmExact->irlgm1f(0)->termStructure(),
1013 d.ccLgmExact->irlgm1f(1)->termStructure());
1014 QuantLib::ext::shared_ptr<BlackCalibrationHelper> tmpGbp = QuantLib::ext::make_shared<FxEqOptionHelper>(
1015 i < d.volstepdatesFx.size() ? d.volstepdatesFx[i] : d.volstepdatesFx.back() + 365, 1.35, d.fxEurGbp,
1016 Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.20)), d.ccLgmExact->irlgm1f(0)->termStructure(),
1017 d.ccLgmExact->irlgm1f(2)->termStructure());
1018 tmpUsd->setPricingEngine(ccLgmFxOptionEngineUsd);
1019 tmpGbp->setPricingEngine(ccLgmFxOptionEngineGbp);
1020 helpersUsd.push_back(tmpUsd);
1021 helpersGbp.push_back(tmpGbp);
1022 }
1023
1024 LevenbergMarquardt lm(1E-8, 1E-8, 1E-8);
1025 EndCriteria ec(1000, 500, 1E-8, 1E-8, 1E-8);
1026
1027 // calibrate USD-EUR FX volatility
1028 d.ccLgmExact->calibrateBsVolatilitiesIterative(CrossAssetModel::AssetType::FX, 0, helpersUsd, lm, ec);
1029 // calibrate GBP-EUR FX volatility
1030 d.ccLgmExact->calibrateBsVolatilitiesIterative(CrossAssetModel::AssetType::FX, 1, helpersGbp, lm, ec);
1031
1032 Real tol = 1E-6;
1033 for (Size i = 0; i < helpersUsd.size(); ++i) {
1034 Real market = helpersUsd[i]->marketValue();
1035 Real model = helpersUsd[i]->modelValue();
1036 Real calibratedVol = d.ccLgmExact->fxbs(0)->parameterValues(0)[i];
1037 if (std::fabs(market - model) > tol) {
1038 BOOST_ERROR("calibration for fx option helper #" << i << " (USD) failed, market premium is " << market
1039 << " while model premium is " << model);
1040 }
1041 // the stochastic rates produce some noise, but do not have a huge
1042 // impact on the effective volatility, so we check that they are
1043 // in line with a cached example (note that the analytic fx option
1044 // pricing engine was checked against MC in another test case)
1045 if (std::fabs(calibratedVol - 0.143) > 0.01) {
1046 BOOST_ERROR("calibrated fx volatility #" << i << " (USD) seems off, expected to be 0.143 +- 0.01, but is "
1047 << calibratedVol);
1048 }
1049 }
1050 for (Size i = 0; i < helpersGbp.size(); ++i) {
1051 Real market = helpersGbp[i]->marketValue();
1052 Real model = helpersGbp[i]->modelValue();
1053 Real calibratedVol = d.ccLgmExact->fxbs(1)->parameterValues(0)[i];
1054 if (std::fabs(market - model) > tol)
1055 BOOST_ERROR("calibration for fx option helper #" << i << " (GBP) failed, market premium is " << market
1056 << " while model premium is " << model);
1057 // see above
1058 if (std::fabs(calibratedVol - 0.193) > 0.01)
1059 BOOST_ERROR("calibrated fx volatility #" << i << " (USD) seems off, expected to be 0.193 +- 0.01, but is "
1060 << calibratedVol);
1061 }
1062
1063 // calibrate the projected model
1064
1065 for (Size i = 0; i < helpersGbp.size(); ++i) {
1066 helpersGbp[i]->setPricingEngine(ccLgmProjectedFxOptionEngineGbp);
1067 }
1068
1069 ccLgmProjected->calibrateBsVolatilitiesIterative(CrossAssetModel::AssetType::FX, 0, helpersGbp, lm, ec);
1070
1071 for (Size i = 0; i < helpersGbp.size(); ++i) {
1072 Real fullModelVol = d.ccLgmExact->fxbs(1)->parameterValues(0)[i];
1073 Real projectedModelVol = ccLgmProjected->fxbs(0)->parameterValues(0)[i];
1074 if (std::fabs(fullModelVol - projectedModelVol) > tol)
1075 BOOST_ERROR("calibrated fx volatility of full model @"
1076 << i << " (" << fullModelVol << ") is inconsistent with that of the projected model ("
1077 << projectedModelVol << ")");
1078 }
1079
1080} // testLgm5fFxCalibration

◆ BOOST_AUTO_TEST_CASE() [7/24]

BOOST_AUTO_TEST_CASE ( testLgm5fFullCalibration  )

Definition at line 1082 of file crossassetmodel.cpp.

1082 {
1083
1084 BOOST_TEST_MESSAGE("Testing full calibration of Ccy LGM 5F model...");
1085
1086 Lgm5fTestData d;
1087
1088 // calibration baskets
1089
1090 std::vector<QuantLib::ext::shared_ptr<BlackCalibrationHelper> > basketEur, basketUsd, basketGbp, basketEurUsd, basketEurGbp;
1091
1092 QuantLib::ext::shared_ptr<IborIndex> euribor6m = QuantLib::ext::make_shared<Euribor>(6 * Months, d.eurYts);
1093 QuantLib::ext::shared_ptr<IborIndex> usdLibor3m = QuantLib::ext::make_shared<USDLibor>(3 * Months, d.usdYts);
1094 QuantLib::ext::shared_ptr<IborIndex> gbpLibor3m = QuantLib::ext::make_shared<GBPLibor>(3 * Months, d.gbpYts);
1095
1096 for (Size i = 0; i <= d.volstepdates.size(); ++i) {
1097 Date tmp = i < d.volstepdates.size() ? d.volstepdates[i] : d.volstepdates.back() + 365;
1098 // EUR: atm+200bp, 150bp normal vol
1099 basketEur.push_back(QuantLib::ext::shared_ptr<SwaptionHelper>(new SwaptionHelper(
1100 tmp, 10 * Years, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.015)), euribor6m, 1 * Years, Thirty360(Thirty360::BondBasis),
1101 Actual360(), d.eurYts, BlackCalibrationHelper::RelativePriceError, 0.04, 1.0, Normal)));
1102 // USD: atm, 20%, lognormal vol
1103 basketUsd.push_back(QuantLib::ext::shared_ptr<SwaptionHelper>(
1104 new SwaptionHelper(tmp, 10 * Years, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.30)), usdLibor3m,
1105 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(), d.usdYts,
1106 BlackCalibrationHelper::RelativePriceError, Null<Real>(), 1.0, ShiftedLognormal, 0.0)));
1107 // GBP: atm-200bp, 10%, shifted lognormal vol with shift = 2%
1108 basketGbp.push_back(QuantLib::ext::shared_ptr<SwaptionHelper>(new SwaptionHelper(
1109 tmp, 10 * Years, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.30)), gbpLibor3m, 1 * Years, Thirty360(Thirty360::BondBasis),
1110 Actual360(), d.usdYts, BlackCalibrationHelper::RelativePriceError, 0.02, 1.0, ShiftedLognormal, 0.02)));
1111 }
1112
1113 for (Size i = 0; i < d.volstepdatesFx.size(); ++i) {
1114 Date tmp = i < d.volstepdatesFx.size() ? d.volstepdatesFx[i] : d.volstepdatesFx.back() + 365;
1115 // EUR-USD: atm, 30% (lognormal) vol
1116 basketEurUsd.push_back(QuantLib::ext::make_shared<FxEqOptionHelper>(
1117 tmp, Null<Real>(), d.fxEurUsd, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.20)), d.eurYts, d.usdYts,
1118 BlackCalibrationHelper::RelativePriceError));
1119 // EUR-GBP: atm, 10% (lognormal) vol
1120 basketEurGbp.push_back(QuantLib::ext::make_shared<FxEqOptionHelper>(
1121 tmp, Null<Real>(), d.fxEurGbp, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.20)), d.eurYts, d.gbpYts,
1122 BlackCalibrationHelper::RelativePriceError));
1123 }
1124
1125 // pricing engines
1126
1127 QuantLib::ext::shared_ptr<PricingEngine> eurSwEng = QuantLib::ext::make_shared<AnalyticLgmSwaptionEngine>(d.ccLgmExact, 0);
1128 QuantLib::ext::shared_ptr<PricingEngine> usdSwEng = QuantLib::ext::make_shared<AnalyticLgmSwaptionEngine>(d.ccLgmExact, 1);
1129 QuantLib::ext::shared_ptr<PricingEngine> gbpSwEng = QuantLib::ext::make_shared<AnalyticLgmSwaptionEngine>(d.ccLgmExact, 2);
1130
1131 QuantLib::ext::shared_ptr<AnalyticCcLgmFxOptionEngine> eurUsdFxoEng =
1132 QuantLib::ext::make_shared<AnalyticCcLgmFxOptionEngine>(d.ccLgmExact, 0);
1133 QuantLib::ext::shared_ptr<AnalyticCcLgmFxOptionEngine> eurGbpFxoEng =
1134 QuantLib::ext::make_shared<AnalyticCcLgmFxOptionEngine>(d.ccLgmExact, 1);
1135
1136 eurUsdFxoEng->cache();
1137 eurGbpFxoEng->cache();
1138
1139 // assign engines to calibration instruments
1140 for (Size i = 0; i < basketEur.size(); ++i) {
1141 basketEur[i]->setPricingEngine(eurSwEng);
1142 }
1143 for (Size i = 0; i < basketUsd.size(); ++i) {
1144 basketUsd[i]->setPricingEngine(usdSwEng);
1145 }
1146 for (Size i = 0; i < basketGbp.size(); ++i) {
1147 basketGbp[i]->setPricingEngine(gbpSwEng);
1148 }
1149 for (Size i = 0; i < basketEurUsd.size(); ++i) {
1150 basketEurUsd[i]->setPricingEngine(eurUsdFxoEng);
1151 }
1152 for (Size i = 0; i < basketEurGbp.size(); ++i) {
1153 basketEurGbp[i]->setPricingEngine(eurGbpFxoEng);
1154 }
1155
1156 // calibrate the model
1157
1158 LevenbergMarquardt lm(1E-8, 1E-8, 1E-8);
1159 EndCriteria ec(1000, 500, 1E-8, 1E-8, 1E-8);
1160
1161 d.ccLgmExact->calibrateIrLgm1fVolatilitiesIterative(0, basketEur, lm, ec);
1162 d.ccLgmExact->calibrateIrLgm1fVolatilitiesIterative(1, basketUsd, lm, ec);
1163 d.ccLgmExact->calibrateIrLgm1fVolatilitiesIterative(2, basketGbp, lm, ec);
1164
1165 d.ccLgmExact->calibrateBsVolatilitiesIterative(CrossAssetModel::AssetType::FX, 0, basketEurUsd, lm, ec);
1166 d.ccLgmExact->calibrateBsVolatilitiesIterative(CrossAssetModel::AssetType::FX, 1, basketEurGbp, lm, ec);
1167
1168 // check the results
1169
1170 Real tol = 1E-6;
1171
1172 for (Size i = 0; i < basketEur.size(); ++i) {
1173 Real model = basketEur[i]->modelValue();
1174 Real market = basketEur[i]->marketValue();
1175 if (std::abs((model - market) / market) > tol)
1176 BOOST_ERROR("calibration failed for instrument #"
1177 << i << " in EUR basket, model value is " << model << " market value is " << market
1178 << " relative error " << std::abs((model - market) / market) << " tolerance " << tol);
1179 }
1180 for (Size i = 0; i < basketUsd.size(); ++i) {
1181 Real model = basketUsd[i]->modelValue();
1182 Real market = basketUsd[i]->marketValue();
1183 if (std::abs((model - market) / market) > tol)
1184 BOOST_ERROR("calibration failed for instrument #"
1185 << i << " in USD basket, model value is " << model << " market value is " << market
1186 << " relative error " << std::abs((model - market) / market) << " tolerance " << tol);
1187 }
1188 for (Size i = 0; i < basketGbp.size(); ++i) {
1189 Real model = basketGbp[i]->modelValue();
1190 Real market = basketGbp[i]->marketValue();
1191 if (std::abs((model - market) / market) > tol)
1192 BOOST_ERROR("calibration failed for instrument #"
1193 << i << " in GBP basket, model value is " << model << " market value is " << market
1194 << " relative error " << std::abs((model - market) / market) << " tolerance " << tol);
1195 }
1196 for (Size i = 0; i < basketEurUsd.size(); ++i) {
1197 Real model = basketEurUsd[i]->modelValue();
1198 Real market = basketEurUsd[i]->marketValue();
1199 if (std::abs((model - market) / market) > tol)
1200 BOOST_ERROR("calibration failed for instrument #"
1201 << i << " in EUR-USD basket, model value is " << model << " market value is " << market
1202 << " relative error " << std::abs((model - market) / market) << " tolerance " << tol);
1203 }
1204 for (Size i = 0; i < basketEurUsd.size(); ++i) {
1205 Real model = basketEurGbp[i]->modelValue();
1206 Real market = basketEurGbp[i]->marketValue();
1207 if (std::abs((model - market) / market) > tol)
1208 BOOST_ERROR("calibration failed for instrument #"
1209 << i << " in EUR-GBP basket, model value is " << model << " market value is " << market
1210 << " relative error " << std::abs((model - market) / market) << " tolerance " << tol);
1211 }
1212}

◆ BOOST_AUTO_TEST_CASE() [8/24]

BOOST_AUTO_TEST_CASE ( testLgm5fMoments  )

Definition at line 1214 of file crossassetmodel.cpp.

1214 {
1215
1216 BOOST_TEST_MESSAGE("Testing analytic moments vs. Euler and exact discretization "
1217 "in Ccy LGM 5F model...");
1218
1219 Lgm5fTestData d;
1220
1221 QuantLib::ext::shared_ptr<StochasticProcess> p_exact = d.ccLgmExact->stateProcess();
1222 QuantLib::ext::shared_ptr<StochasticProcess> p_euler = d.ccLgmEuler->stateProcess();
1223
1224 Real T = 10.0; // horizon at which we compare the moments
1225 Size steps = static_cast<Size>(T * 10.0); // number of simulation steps
1226 Size paths = 25000; // number of paths
1227
1228 Array e_an = p_exact->expectation(0.0, p_exact->initialValues(), T);
1229 Matrix v_an = p_exact->covariance(0.0, p_exact->initialValues(), T);
1230
1231 TimeGrid grid(T, steps);
1232
1233 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_euler)) {
1234 tmp->resetCache(grid.size() - 1);
1235 }
1236 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_exact)) {
1237 tmp->resetCache(grid.size() - 1);
1238 }
1239 MultiPathGeneratorSobolBrownianBridge pgen(p_euler, grid);
1240 MultiPathGeneratorSobolBrownianBridge pgen2(p_exact, grid);
1241
1242 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > e_eu[5], e_eu2[5];
1243 accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > v_eu[5][5], v_eu2[5][5];
1244
1245 for (Size i = 0; i < paths; ++i) {
1246 Sample<MultiPath> path = pgen.next();
1247 Sample<MultiPath> path2 = pgen2.next();
1248 for (Size ii = 0; ii < 5; ++ii) {
1249 Real cii = path.value[ii].back();
1250 Real cii2 = path2.value[ii].back();
1251 e_eu[ii](cii);
1252 e_eu2[ii](cii2);
1253 for (Size jj = 0; jj <= ii; ++jj) {
1254 Real cjj = path.value[jj].back();
1255 v_eu[ii][jj](cii, covariate1 = cjj);
1256 Real cjj2 = path2.value[jj].back();
1257 v_eu2[ii][jj](cii2, covariate1 = cjj2);
1258 }
1259 }
1260 }
1261
1262 Real errTolLd[] = { 0.2E-4, 0.2E-4, 0.2E-4, 10.0E-4, 10.0E-4 };
1263
1264 for (Size i = 0; i < 5; ++i) {
1265 // check expectation against analytical calculation (Euler)
1266 if (std::fabs(mean(e_eu[i]) - e_an[i]) > errTolLd[i]) {
1267 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
1268 << ") is inconsistent with numerical value (Euler "
1269 "discretization, "
1270 << mean(e_eu[i]) << "), error is "
1271 << e_an[i] - mean(e_eu[i]) << " tolerance is "
1272 << errTolLd[i]);
1273 }
1274 // check expectation against analytical calculation (exact disc)
1275 if (std::fabs(mean(e_eu2[i]) - e_an[i]) > errTolLd[i]) {
1276 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
1277 << ") is inconsistent with numerical value (Exact "
1278 "discretization, "
1279 << mean(e_eu2[i]) << "), error is "
1280 << e_an[i] - mean(e_eu2[i]) << " tolerance is "
1281 << errTolLd[i]);
1282 }
1283 }
1284
1285 // we have to deal with different natures of volatility
1286 // for ir (normal) and fx (ln) so different error
1287 // tolerances apply
1288 Real tollNormal = 0.1E-4; // ir-ir
1289 Real tolMixed = 0.25E-4; // ir-fx
1290 Real tolLn = 8.0E-4; // fx-fx
1291 Real tol; // set below
1292
1293 for (Size i = 0; i < 5; ++i) {
1294 for (Size j = 0; j <= i; ++j) {
1295 if (i < 3) {
1296 tol = tollNormal;
1297 } else {
1298 if (j < 3) {
1299 tol = tolMixed;
1300 } else {
1301 tol = tolLn;
1302 }
1303 }
1304 if (std::fabs(boost::accumulators::covariance(v_eu[i][j]) - v_an[i][j]) > tol) {
1305 BOOST_ERROR("analytical covariance at ("
1306 << i << "," << j << ") (" << v_an[i][j]
1307 << ") is inconsistent with numerical "
1308 "value (Euler discretization, "
1309 << boost::accumulators::covariance(v_eu[i][j]) << "), error is "
1310 << v_an[i][j] - boost::accumulators::covariance(v_eu[i][j]) << " tolerance is " << tol);
1311 }
1312 if (std::fabs(boost::accumulators::covariance(v_eu2[i][j]) - v_an[i][j]) > tol) {
1313 BOOST_ERROR("analytical covariance at ("
1314 << i << "," << j << ") (" << v_an[i][j]
1315 << ") is inconsistent with numerical "
1316 "value (Exact discretization, "
1317 << boost::accumulators::covariance(v_eu2[i][j]) << "), error is "
1318 << v_an[i][j] - boost::accumulators::covariance(v_eu2[i][j]) << " tolerance is " << tol);
1319 }
1320 }
1321 }
1322
1323} // testLgm5fMoments
QuantExt::MultiPathGeneratorSobolBrownianBridge
Instantiation using SobolBrownianGenerator from models/marketmodels/browniangenerators.
Definition: multipathgeneratorbase.hpp:161
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [9/24]

BOOST_AUTO_TEST_CASE ( testLgmGsrEquivalence  )

Definition at line 1325 of file crossassetmodel.cpp.

1325 {
1326
1327 BOOST_TEST_MESSAGE("Testing equivalence of GSR and LGM models...");
1328
1329 SavedSettings backup;
1330
1331 Date evalDate(12, January, 2015);
1332 Settings::instance().evaluationDate() = evalDate;
1333 Handle<YieldTermStructure> yts(QuantLib::ext::make_shared<FlatForward>(evalDate, 0.02, Actual365Fixed()));
1334
1335 Real T[] = { 10.0, 20.0, 50.0, 100.0 };
1336 Real sigma[] = { 0.0050, 0.01, 0.02 };
1337 Real kappa[] = { -0.02, -0.01, 0.0, 0.03, 0.07 };
1338
1339 for (Size i = 0; i < LENGTH(T); ++i) {
1340 for (Size j = 0; j < LENGTH(sigma); ++j) {
1341 for (Size k = 0; k < LENGTH(kappa); ++k) {
1342
1343 std::vector<Date> stepDates;
1344 std::vector<Real> sigmas(1, sigma[j]);
1345
1346 QuantLib::ext::shared_ptr<Gsr> gsr = QuantLib::ext::make_shared<Gsr>(yts, stepDates, sigmas, kappa[k], T[i]);
1347
1348 Array stepTimes_a(0);
1349 Array sigmas_a(1, sigma[j]);
1350 Array kappas_a(1, kappa[k]);
1351
1352 // for shift = -H(T) we change the LGM measure to the T forward
1353 // measure effectively
1354 Real shift = close_enough(kappa[k], 0.0) ? -T[i] : -(1.0 - std::exp(-kappa[k] * T[i])) / kappa[k];
1355 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm_p =
1356 QuantLib::ext::make_shared<IrLgm1fPiecewiseConstantHullWhiteAdaptor>(EURCurrency(), yts, stepTimes_a,
1357 sigmas_a, stepTimes_a, kappas_a);
1358 lgm_p->shift() = shift;
1359
1360 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> lgm = QuantLib::ext::make_shared<LinearGaussMarkovModel>(lgm_p);
1361
1362 QuantLib::ext::shared_ptr<StochasticProcess1D> gsr_process = gsr->stateProcess();
1363 QuantLib::ext::shared_ptr<StochasticProcess1D> lgm_process =
1364 QuantLib::ext::dynamic_pointer_cast<StochasticProcess1D>(lgm->stateProcess());
1365
1366 Size N = 10000; // number of paths
1367 Size seed = 123456;
1368 Size steps = 1; // one large step
1369 Real T2 = T[i] - 5.0; // we check a distribution at this time
1370
1371 TimeGrid grid(T2, steps);
1372
1373 PseudoRandom::rsg_type sg = PseudoRandom::make_sequence_generator(steps * 1, seed);
1374 PathGenerator<PseudoRandom::rsg_type> pgen_gsr(gsr_process, grid, sg, false);
1375 PathGenerator<PseudoRandom::rsg_type> pgen_lgm(lgm_process, grid, sg, false);
1376
1377 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean>, tag::variance> > stat_lgm, stat_gsr;
1378
1379 Real tol = 1.0E-12;
1380 for (Size ii = 0; ii < N; ++ii) {
1381 Sample<Path> path_lgm = pgen_lgm.next();
1382 Sample<Path> path_gsr = pgen_gsr.next();
1383 Real yGsr = (path_gsr.value.back() - gsr_process->expectation(0.0, 0.0, T2)) /
1384 gsr_process->stdDeviation(0.0, 0.0, T2);
1385 Real xLgm = path_lgm.value.back();
1386 Real gsrRate = -std::log(gsr->zerobond(T2 + 1.0, T2, yGsr));
1387 // it's nice to have uniform interfaces in all models ...
1388 Real lgmRate = -std::log(lgm->discountBond(T2, T2 + 1.0, xLgm));
1389 stat_gsr(gsrRate);
1390 stat_lgm(lgmRate);
1391 // check pathwise identity
1392 if (std::fabs(gsrRate - lgmRate) >= tol) {
1393 BOOST_ERROR("lgm rate (" << lgmRate << ") deviates from gsr rate (" << gsrRate << ") on path #"
1394 << i);
1395 }
1396 }
1397
1398 // effectively we are checking a pathwise identity
1399 // here as well, but the statistics seems to better
1400 // summarize a possible problem, so we output differences
1401 // in the mean as well
1402 if (std::fabs(mean(stat_gsr) - mean(stat_lgm)) > tol ||
1403 std::fabs(boost::accumulators::variance(stat_gsr) - boost::accumulators::variance(stat_lgm)) >
1404 tol) {
1405 BOOST_ERROR("failed to verify LGM-GSR equivalence, "
1406 "(mean,variance) of zero rate is ("
1407 << mean(stat_gsr) << "," << boost::accumulators::variance(stat_gsr) << ") for GSR, ("
1408 << mean(stat_lgm) << "," << boost::accumulators::variance(stat_lgm)
1409 << ") for LGM, for T=" << T[i] << ", sigma=" << sigma[j] << ", kappa=" << kappa[k]
1410 << ", shift=" << shift);
1411 }
1412 }
1413 }
1414 }
1415
1416} // testLgmGsrEquivalence
LENGTH
#define LENGTH(a)
Definition: utilities.hpp:27
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [10/24]

BOOST_AUTO_TEST_CASE ( testLgmMcWithShift  )

Definition at line 1418 of file crossassetmodel.cpp.

1418 {
1419 BOOST_TEST_MESSAGE("Testing LGM1F Monte Carlo simulation with shifted H...");
1420
1421 // cashflow time
1422 Real T = 50.0;
1423
1424 // shift horizons
1425 Real T_shift[] = { 0.0, 10.0, 20.0, 30.0, 40.0, 50.0 };
1426
1427 // tolerances for error of mean
1428 Real eom_tol[] = { 0.17, 0.05, 0.02, 0.01, 0.005, 1.0E-12 };
1429
1430 Handle<YieldTermStructure> yts(QuantLib::ext::make_shared<FlatForward>(0, NullCalendar(), 0.02, Actual365Fixed()));
1431
1432 QuantLib::ext::shared_ptr<IrLgm1fParametrization> lgm =
1433 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(EURCurrency(), yts, 0.01, 0.01);
1434 QuantLib::ext::shared_ptr<StochasticProcess> p = QuantLib::ext::make_shared<IrLgm1fStateProcess>(lgm);
1435
1436 QuantLib::ext::shared_ptr<LinearGaussMarkovModel> model = QuantLib::ext::make_shared<LinearGaussMarkovModel>(lgm);
1437
1438 Size steps = 1;
1439 Size paths = 10000;
1440 Size seed = 42;
1441 TimeGrid grid(T, steps);
1442
1443 MultiPathGeneratorMersenneTwister pgen(p, grid, seed, true);
1444
1445 for (Size ii = 0; ii < LENGTH(T_shift); ++ii) {
1446
1447 lgm->shift() = -(1.0 - exp(-0.01 * T_shift[ii])) / 0.01;
1448
1449 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > e_eu;
1450
1451 for (Size i = 0; i < paths; ++i) {
1452 Sample<MultiPath> path = pgen.next();
1453 Sample<MultiPath> path_a = pgen.next();
1454 e_eu(1.0 / model->numeraire(T, path.value[0].back()));
1455 e_eu(1.0 / model->numeraire(T, path_a.value[0].back()));
1456 }
1457
1458 Real discount = yts->discount(T);
1459
1460 if (error_of<tag::mean>(e_eu) / discount > eom_tol[ii]) {
1461 BOOST_ERROR("estimated error of mean for shifted mc simulation with shift "
1462 << T_shift[ii] << " can not be verified (" << error_of<tag::mean>(e_eu) / discount
1463 << "), tolerance is 1E-8");
1464 }
1465
1466 if (std::fabs(mean(e_eu) / discount - 1.0) > eom_tol[ii]) {
1467 BOOST_ERROR("estimated error for shifted mc simulation with shift " << T_shift[ii]
1468 << " can not "
1469 "be verified ("
1470 << mean(e_eu) / discount - 1.0
1471 << "), tolerance is 1E-8");
1472 }
1473 }
1474
1475} // testLgmMcWithShift
QuantExt::exp
CompiledFormula exp(CompiledFormula x)
Definition: compiledformula.cpp:151
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [11/24]

BOOST_AUTO_TEST_CASE ( testIrFxCrCirppMartingaleProperty  )

Definition at line 1477 of file crossassetmodel.cpp.

1477 {
1478
1479 BOOST_TEST_MESSAGE("Testing martingale property in ir-fx-cr(lgm)-cf(cir++) model for "
1480 "Euler and exact discretizations...");
1481
1482 IrFxCrModelTestData d(false);
1483 IrFxCrModelTestData d_cirpp(true);
1484
1485 QuantLib::ext::shared_ptr<StochasticProcess> process1 = d.modelExact->stateProcess();
1486 QuantLib::ext::shared_ptr<StochasticProcess> process2 = d_cirpp.modelEuler->stateProcess();
1487
1488 Size n = 10000; // number of paths
1489 Size seed = 18; // rng seed
1490 Time T = 2.0; // maturity of payoff
1491 Time T2 = 20.0; // zerobond maturity
1492 Size steps = static_cast<Size>(T * 24); // number of steps taken (euler / Brigo-Alfonsi)
1493
1494 LowDiscrepancy::rsg_type sg1 = LowDiscrepancy::make_sequence_generator(process1->factors() * 1, seed);
1495 LowDiscrepancy::rsg_type sg2 = LowDiscrepancy::make_sequence_generator(process2->factors() * steps, seed);
1496
1497 TimeGrid grid1(T, 1);
1498 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process1)) {
1499 tmp->resetCache(grid1.size() - 1);
1500 }
1501 MultiPathGenerator<LowDiscrepancy::rsg_type> pg1(process1, grid1, sg1, false);
1502 TimeGrid grid2(T, steps);
1503 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process2)) {
1504 tmp->resetCache(grid2.size() - 1);
1505 }
1506 MultiPathGenerator<LowDiscrepancy::rsg_type> pg2(process2, grid2, sg2, false);
1507
1508 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > eurzb1, usdzb1, gbpzb1, n1eur1, n2usd1,
1509 n3gbp1;
1510 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > eurzb2, usdzb2, gbpzb2, n1eur2, n2usd2,
1511 n3gbp2, n1cir2, n2cir2, n3cir2;
1512
1513 for (Size j = 0; j < n; ++j) {
1514 Sample<MultiPath> path1 = pg1.next();
1515 Sample<MultiPath> path2 = pg2.next();
1516 Size l1 = path1.value[0].length() - 1;
1517 Size l2 = path2.value[0].length() - 1;
1518 Real zeur1 = path1.value[0][l1];
1519 Real zusd1 = path1.value[1][l1];
1520 Real zgbp1 = path1.value[2][l1];
1521 Real fxusd1 = std::exp(path1.value[3][l1]);
1522 Real fxgbp1 = std::exp(path1.value[4][l1]);
1523 Real crzn11 = path1.value[5][l1];
1524 Real cryn11 = path1.value[6][l1];
1525 Real crzn21 = path1.value[7][l1];
1526 Real cryn21 = path1.value[8][l1];
1527 Real crzn31 = path1.value[9][l1];
1528 Real cryn31 = path1.value[10][l1];
1529 Real zeur2 = path2.value[0][l2];
1530 Real zusd2 = path2.value[1][l2];
1531 Real zgbp2 = path2.value[2][l2];
1532 Real fxusd2 = std::exp(path2.value[3][l2]);
1533 Real fxgbp2 = std::exp(path2.value[4][l2]);
1534 Real crzn12 = path2.value[5][l2]; // CR LGM
1535 Real cryn12 = path2.value[6][l2];
1536 Real crzn22 = path2.value[7][l2];
1537 Real cryn22 = path2.value[8][l2];
1538 Real crzn32 = path2.value[9][l2];
1539 Real cryn32 = path2.value[10][l2];
1540 Real ciry12 = path2.value[11][l2]; // CR CIR++
1541 Real cirn12 = path2.value[12][l2];
1542 Real ciry22 = path2.value[13][l2];
1543 Real cirn22 = path2.value[14][l2];
1544 Real ciry32 = path2.value[15][l2];
1545 Real cirn32 = path2.value[16][l2];
1546
1547 // EUR zerobond
1548 eurzb1(d.modelExact->discountBond(0, T, T2, zeur1) / d.modelExact->numeraire(0, T, zeur1));
1549 // USD zerobond
1550 usdzb1(d.modelExact->discountBond(1, T, T2, zusd1) * fxusd1 / d.modelExact->numeraire(0, T, zeur1));
1551 // GBP zerobond
1552 gbpzb1(d.modelExact->discountBond(2, T, T2, zgbp1) * fxgbp1 / d.modelExact->numeraire(0, T, zeur1));
1553 // EUR defaultable zerobond for name 1
1554 std::pair<Real, Real> sn11 = d.modelExact->crlgm1fS(0, 0, T, T2, crzn11, cryn11);
1555 n1eur1(sn11.first * sn11.second * d.modelExact->discountBond(0, T, T2, zeur1) / d.modelExact->numeraire(0, T, zeur1));
1556 // USD defaultable zerobond for name 2
1557 std::pair<Real, Real> sn21 = d.modelExact->crlgm1fS(1, 1, T, T2, crzn21, cryn21);
1558 n2usd1(sn21.first * sn21.second * d.modelExact->discountBond(1, T, T2, zusd1) * fxusd1 /
1559 d.modelExact->numeraire(0, T, zeur1));
1560 // GBP defaultable zerobond for name 3
1561 std::pair<Real, Real> sn31 = d.modelExact->crlgm1fS(2, 2, T, T2, crzn31, cryn31);
1562 n3gbp1(sn31.first * sn31.second * d.modelExact->discountBond(2, T, T2, zgbp1) * fxgbp1 /
1563 d.modelExact->numeraire(0, T, zeur1));
1564
1565 // EUR zerobond
1566 eurzb2(d_cirpp.modelEuler->discountBond(0, T, T2, zeur2) / d_cirpp.modelEuler->numeraire(0, T, zeur2));
1567 // USD zerobond
1568 usdzb2(d_cirpp.modelEuler->discountBond(1, T, T2, zusd2) * fxusd2 / d_cirpp.modelEuler->numeraire(0, T, zeur2));
1569 // GBP zerobond
1570 gbpzb2(d_cirpp.modelEuler->discountBond(2, T, T2, zgbp2) * fxgbp2 / d_cirpp.modelEuler->numeraire(0, T, zeur2));
1571 // EUR defaultable zerobond for name 1
1572 std::pair<Real, Real> sn12 = d_cirpp.modelEuler->crlgm1fS(0, 0, T, T2, crzn12, cryn12);
1573 n1eur2(sn12.first * sn12.second * d_cirpp.modelEuler->discountBond(0, T, T2, zeur2) / d_cirpp.modelEuler->numeraire(0, T, zeur2));
1574 // USD defaultable zerobond for name 2
1575 std::pair<Real, Real> sn22 = d_cirpp.modelEuler->crlgm1fS(1, 1, T, T2, crzn22, cryn22);
1576 n2usd2(sn22.first * sn22.second * d_cirpp.modelEuler->discountBond(1, T, T2, zusd2) * fxusd2 /
1577 d_cirpp.modelEuler->numeraire(0, T, zeur2));
1578 // GBP defaultable zerobond for name 3
1579 std::pair<Real, Real> sn32 = d_cirpp.modelEuler->crlgm1fS(2, 2, T, T2, crzn32, cryn32);
1580 n3gbp2(sn32.first * sn32.second * d_cirpp.modelEuler->discountBond(2, T, T2, zgbp2) * fxgbp2 /
1581 d_cirpp.modelEuler->numeraire(0, T, zeur2));
1582 // EUR defaultable zerobond for name 1 (CIR++)
1583 std::pair<Real, Real> sc12 = d_cirpp.modelEuler->crcirppS(3, T, T2, ciry12, cirn12);
1584 n1cir2(sc12.first * sc12.second * d_cirpp.modelEuler->discountBond(0, T, T2, zeur2) / d_cirpp.modelEuler->numeraire(0, T, zeur2));
1585 // USD defaultable zerobond for name 2 (CIR++)
1586 std::pair<Real, Real> sc22 = d_cirpp.modelEuler->crcirppS(4, T, T2, ciry22, cirn22);
1587 n2cir2(sc22.first * sc22.second * d_cirpp.modelEuler->discountBond(1, T, T2, zusd2) * fxusd2 /
1588 d_cirpp.modelEuler->numeraire(0, T, zeur2));
1589 // GBP defaultable zerobond for name 3 (CIR++)
1590 std::pair<Real, Real> sc32 = d_cirpp.modelEuler->crcirppS(5, T, T2, ciry32, cirn32);
1591 n3cir2(sc32.first * sc32.second * d_cirpp.modelEuler->discountBond(2, T, T2, zgbp2) * fxgbp2 /
1592 d_cirpp.modelEuler->numeraire(0, T, zeur2));
1593 }
1594
1595 BOOST_TEST_MESSAGE("EXACT:");
1596 BOOST_TEST_MESSAGE("EUR zb = " << mean(eurzb1) << " +- " << error_of<tag::mean>(eurzb1) << " vs analytical "
1597 << d.eurYts->discount(T2));
1598 BOOST_TEST_MESSAGE("USD zb = " << mean(usdzb1) << " +- " << error_of<tag::mean>(usdzb1) << " vs analytical "
1599 << d.usdYts->discount(T2) * d.fxEurUsd->value());
1600 BOOST_TEST_MESSAGE("GBP zb = " << mean(gbpzb1) << " +- " << error_of<tag::mean>(gbpzb1) << " vs analytical "
1601 << d.gbpYts->discount(T2) * d.fxEurGbp->value());
1602 BOOST_TEST_MESSAGE("N1 zb EUR = " << mean(n1eur1) << " +- " << error_of<tag::mean>(n1eur1) << " vs analytical "
1603 << d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2));
1604 BOOST_TEST_MESSAGE("N2 zb USD = " << mean(n2usd1) << " +- " << error_of<tag::mean>(n2usd1) << " vs analytical "
1605 << d.fxEurUsd->value() * d.usdYts->discount(T2) *
1606 d.n2Ts->survivalProbability(T2));
1607 BOOST_TEST_MESSAGE("N3 zb GBP = " << mean(n3gbp1) << " +- " << error_of<tag::mean>(n3gbp1) << " vs analytical "
1608 << d.fxEurGbp->value() * d.gbpYts->discount(T2) *
1609 d.n3Ts->survivalProbability(T2));
1610
1611 BOOST_TEST_MESSAGE("\nEULER:");
1612 BOOST_TEST_MESSAGE("EUR zb = " << mean(eurzb2) << " +- " << error_of<tag::mean>(eurzb2) << " vs analytical "
1613 << d_cirpp.eurYts->discount(T2));
1614 BOOST_TEST_MESSAGE("USD zb = " << mean(usdzb2) << " +- " << error_of<tag::mean>(usdzb2) << " vs analytical "
1615 << d_cirpp.usdYts->discount(T2) * d_cirpp.fxEurUsd->value());
1616 BOOST_TEST_MESSAGE("GBP zb = " << mean(gbpzb2) << " +- " << error_of<tag::mean>(gbpzb2) << " vs analytical "
1617 << d_cirpp.gbpYts->discount(T2) * d_cirpp.fxEurGbp->value());
1618 BOOST_TEST_MESSAGE("N1 zb EUR = " << mean(n1eur2) << " +- " << error_of<tag::mean>(n1eur2) << " vs analytical "
1619 << d_cirpp.eurYts->discount(T2) * d_cirpp.n1Ts->survivalProbability(T2));
1620 BOOST_TEST_MESSAGE("N2 zb USD = " << mean(n2usd2) << " +- " << error_of<tag::mean>(n2usd2) << " vs analytical "
1621 << d_cirpp.fxEurUsd->value() * d_cirpp.usdYts->discount(T2) *
1622 d_cirpp.n2Ts->survivalProbability(T2));
1623 BOOST_TEST_MESSAGE("N3 zb GBP = " << mean(n3gbp2) << " +- " << error_of<tag::mean>(n3gbp2) << " vs analytical "
1624 << d_cirpp.fxEurGbp->value() * d_cirpp.gbpYts->discount(T2) *
1625 d_cirpp.n3Ts->survivalProbability(T2));
1626 BOOST_TEST_MESSAGE("N1 zb EUR = " << mean(n1cir2) << " +- " << error_of<tag::mean>(n1cir2) << " vs analytical "
1627 << d_cirpp.eurYts->discount(T2) * d_cirpp.n1Ts->survivalProbability(T2));
1628 BOOST_TEST_MESSAGE("N2 zb USD = " << mean(n2cir2) << " +- " << error_of<tag::mean>(n2cir2) << " vs analytical "
1629 << d_cirpp.fxEurUsd->value() * d_cirpp.usdYts->discount(T2) *
1630 d_cirpp.n2Ts->survivalProbability(T2));
1631 BOOST_TEST_MESSAGE("N3 zb GBP = " << mean(n3cir2) << " +- " << error_of<tag::mean>(n3cir2) << " vs analytical "
1632 << d_cirpp.fxEurGbp->value() * d_cirpp.gbpYts->discount(T2) *
1633 d_cirpp.n3Ts->survivalProbability(T2));
1634
1635 Real tol1 = 2.0E-4; // EXACT
1636 Real tol2 = 12.0E-4; // EULER
1637
1638 Real ev = d.eurYts->discount(T2);
1639 if (std::abs(mean(eurzb1) - ev) > tol1)
1640 BOOST_TEST_ERROR("Martingale test failed for eurzb (exact discr.), expected " << ev << ", got " << mean(eurzb1)
1641 << ", tolerance " << tol1);
1642 ev = d.usdYts->discount(T2) * d.fxEurUsd->value();
1643 if (std::abs(mean(usdzb1) - ev) > tol1)
1644 BOOST_TEST_ERROR("Martingale test failed for usdzb (exact discr.), expected " << ev << ", got " << mean(usdzb1)
1645 << ", tolerance " << tol1);
1646 ev = d.gbpYts->discount(T2) * d.fxEurGbp->value();
1647 if (std::abs(mean(gbpzb1) - ev) > tol1)
1648 BOOST_TEST_ERROR("Martingale test failed for gbpzb (exact discr.), expected " << ev << ", got " << mean(gbpzb1)
1649 << ", tolerance " << tol1);
1650 ev = d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2);
1651 if (std::abs(mean(n1eur1) - ev) > tol1)
1652 BOOST_TEST_ERROR("Martingale test failed for n1eur (exact discr.), expected " << ev << ", got " << mean(n1eur1)
1653 << ", tolerance " << tol1);
1654 ev = d.fxEurUsd->value() * d.usdYts->discount(T2) * d.n2Ts->survivalProbability(T2);
1655 if (std::abs(mean(n2usd1) - ev) > tol1)
1656 BOOST_TEST_ERROR("Martingale test failed for n2usd (exact discr.), expected " << ev << ", got " << mean(n2usd1)
1657 << ", tolerance " << tol1);
1658 ev = d.fxEurGbp->value() * d.gbpYts->discount(T2) * d.n3Ts->survivalProbability(T2);
1659 if (std::abs(mean(n3gbp1) - ev) > tol1)
1660 BOOST_TEST_ERROR("Martingale test failed for n3gbp (exact discr.), expected " << ev << ", got " << mean(n3gbp1)
1661 << ", tolerance " << tol1);
1662
1663 ev = d_cirpp.eurYts->discount(T2);
1664 if (std::abs(mean(eurzb2) - ev) > tol2)
1665 BOOST_TEST_ERROR("Martingale test failed for eurzb (Euler discr.), expected " << ev << ", got " << mean(eurzb2)
1666 << ", tolerance " << tol2);
1667 ev = d_cirpp.usdYts->discount(T2) * d_cirpp.fxEurUsd->value();
1668 if (std::abs(mean(usdzb2) - ev) > tol2)
1669 BOOST_TEST_ERROR("Martingale test failed for usdzb (Euler discr.), expected "
1670 << ev << ", got " << mean(usdzb2) << ", tolerance " << tol2 * error_of<tag::mean>(usdzb2));
1671 ev = d_cirpp.gbpYts->discount(T2) * d_cirpp.fxEurGbp->value();
1672 if (std::abs(mean(gbpzb2) - ev) > tol2)
1673 BOOST_TEST_ERROR("Martingale test failed for gbpzb (Euler discr.), expected " << ev << ", got " << mean(gbpzb2)
1674 << ", tolerance " << tol2);
1675
1676 // CR LGM
1677 ev = d_cirpp.eurYts->discount(T2) * d_cirpp.n1Ts->survivalProbability(T2);
1678 if (std::abs(mean(n1eur2) - ev) > tol2)
1679 BOOST_TEST_ERROR("Martingale test failed for n1eur (Euler discr.), expected " << ev << ", got " << mean(n1eur2)
1680 << ", tolerance " << tol2);
1681 ev = d_cirpp.fxEurUsd->value() * d_cirpp.usdYts->discount(T2) * d_cirpp.n2Ts->survivalProbability(T2);
1682 if (std::abs(mean(n2usd2) - ev) > tol2)
1683 BOOST_TEST_ERROR("Martingale test failed for n2usd (Euler discr.), expected " << ev << ", got " << mean(n2usd2)
1684 << ", tolerance " << tol2);
1685 ev = d_cirpp.fxEurGbp->value() * d_cirpp.gbpYts->discount(T2) * d_cirpp.n3Ts->survivalProbability(T2);
1686 if (std::abs(mean(n3gbp2) - ev) > tol2)
1687 BOOST_TEST_ERROR("Martingale test failed for n3gbp (Euler discr.), expected " << ev << ", got " << mean(n3gbp2)
1688 << ", tolerance " << tol2);
1689
1690 // CR CIR
1691 ev = d_cirpp.eurYts->discount(T2) * d_cirpp.n1Ts->survivalProbability(T2);
1692 if (std::abs(mean(n1cir2) - ev) > tol2)
1693 BOOST_TEST_ERROR("Martingale test failed for n1cir (Euler discr.), expected " << ev << ", got " << mean(n1cir2)
1694 << ", tolerance " << tol2);
1695 ev = d_cirpp.fxEurUsd->value() * d_cirpp.usdYts->discount(T2) * d_cirpp.n2Ts->survivalProbability(T2);
1696 if (std::abs(mean(n2cir2) - ev) > tol2)
1697 BOOST_TEST_ERROR("Martingale test failed for n2cir2 (Euler discr.), expected " << ev << ", got " << mean(n2cir2)
1698 << ", tolerance " << tol2);
1699 ev = d_cirpp.fxEurGbp->value() * d_cirpp.gbpYts->discount(T2) * d_cirpp.n3Ts->survivalProbability(T2);
1700 if (std::abs(mean(n3cir2) - ev) > tol2)
1701 BOOST_TEST_ERROR("Martingale test failed for n3cir2 (Euler discr.), expected " << ev << ", got " << mean(n3cir2)
1702 << ", tolerance " << tol2);
1703
1704} // testIrFxCrMartingaleProperty

◆ BOOST_AUTO_TEST_CASE() [12/24]

BOOST_AUTO_TEST_CASE ( testIrFxCrMoments  )

Definition at line 1706 of file crossassetmodel.cpp.

1706 {
1707
1708 BOOST_TEST_MESSAGE("Testing analytic moments vs. Euler and exact discretization "
1709 "in ir-fx-cr model...");
1710
1711 IrFxCrModelTestData d;
1712
1713 QuantLib::ext::shared_ptr<StochasticProcess> p_exact = d.modelExact->stateProcess();
1714 QuantLib::ext::shared_ptr<StochasticProcess> p_euler = d.modelEuler->stateProcess();
1715
1716 Real T = 2.0; // horizon at which we compare the moments
1717 Size steps = static_cast<Size>(T * 10); // number of simulation steps (Euler and exact)
1718 Size paths = 50000; // number of paths
1719
1720 Array e_an = p_exact->expectation(0.0, p_exact->initialValues(), T);
1721 Matrix v_an = p_exact->covariance(0.0, p_exact->initialValues(), T);
1722
1723 Size seed = 18;
1724 TimeGrid grid(T, steps);
1725
1726 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_exact)) {
1727 tmp->resetCache(grid.size() - 1);
1728 }
1729 MultiPathGeneratorSobolBrownianBridge pgen(p_euler, grid, SobolBrownianGenerator::Diagonal, seed,
1730 SobolRsg::JoeKuoD7);
1731 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_euler)) {
1732 tmp->resetCache(grid.size() - 1);
1733 }
1734 MultiPathGeneratorSobolBrownianBridge pgen2(p_exact, grid, SobolBrownianGenerator::Diagonal, seed,
1735 SobolRsg::JoeKuoD7);
1736
1737 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > e_eu[11], e_eu2[11];
1738 accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > v_eu[11][11], v_eu2[11][11];
1739
1740 for (Size i = 0; i < paths; ++i) {
1741 Sample<MultiPath> path = pgen.next();
1742 Sample<MultiPath> path2 = pgen2.next();
1743 for (Size ii = 0; ii < 11; ++ii) {
1744 Real cii = path.value[ii].back();
1745 Real cii2 = path2.value[ii].back();
1746 e_eu[ii](cii);
1747 e_eu2[ii](cii2);
1748 for (Size jj = 0; jj <= ii; ++jj) {
1749 Real cjj = path.value[jj].back();
1750 v_eu[ii][jj](cii, covariate1 = cjj);
1751 Real cjj2 = path2.value[jj].back();
1752 v_eu2[ii][jj](cii2, covariate1 = cjj2);
1753 }
1754 }
1755 }
1756
1757 for (Size i = 0; i < 11; ++i) {
1758 BOOST_TEST_MESSAGE("E_" << i << " " << e_an[i] << " " << mean(e_eu[i]) << " " << mean(e_eu2[i]));
1759 }
1760 BOOST_TEST_MESSAGE("==================");
1761
1762 BOOST_TEST_MESSAGE("one step analytical");
1763 for (Size i = 0; i < 11; ++i) {
1764 std::ostringstream tmp;
1765 for (Size j = 0; j <= i; ++j) {
1766 tmp << v_an[i][j] << " ";
1767 }
1768 BOOST_TEST_MESSAGE(tmp.str());
1769 }
1770 BOOST_TEST_MESSAGE("==================");
1771
1772 BOOST_TEST_MESSAGE("euler");
1773 for (Size i = 0; i < 11; ++i) {
1774 std::ostringstream tmp;
1775 for (Size j = 0; j <= i; ++j) {
1776 tmp << boost::accumulators::covariance(v_eu[i][j]) << " ";
1777 }
1778 BOOST_TEST_MESSAGE(tmp.str());
1779 }
1780 BOOST_TEST_MESSAGE("==================");
1781
1782 BOOST_TEST_MESSAGE("exact");
1783 for (Size i = 0; i < 11; ++i) {
1784 std::ostringstream tmp;
1785 for (Size j = 0; j <= i; ++j) {
1786 tmp << boost::accumulators::covariance(v_eu2[i][j]) << " ";
1787 }
1788 BOOST_TEST_MESSAGE(tmp.str());
1789 }
1790 BOOST_TEST_MESSAGE("==================");
1791
1792 Real errTolLd[] = { 0.5E-4, 0.5E-4, 0.5E-4, 10.0E-4, 10.0E-4, 0.7E-4, 0.7E-4, 0.7E-4, 0.7E-4, 0.7E-4, 0.7E-4 };
1793
1794 for (Size i = 0; i < 11; ++i) {
1795 // check expectation against analytical calculation (Euler)
1796 if (std::fabs(mean(e_eu[i]) - e_an[i]) > errTolLd[i]) {
1797 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
1798 << ") is inconsistent with numerical value (Euler "
1799 "discretization, "
1800 << mean(e_eu[i]) << "), error is "
1801 << e_an[i] - mean(e_eu[i]) << " tolerance is "
1802 << errTolLd[i]);
1803 }
1804 // check expectation against analytical calculation (exact disc)
1805 if (std::fabs(mean(e_eu2[i]) - e_an[i]) > errTolLd[i]) {
1806 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
1807 << ") is inconsistent with numerical value (exact "
1808 "discretization, "
1809 << mean(e_eu2[i]) << "), error is "
1810 << e_an[i] - mean(e_eu2[i]) << " tolerance is "
1811 << errTolLd[i]);
1812 }
1813 }
1814
1815 // this is a bit rough compared to the more differentiated test
1816 // of the IR-FX model ...
1817 Real tol = 10.0E-4;
1818
1819 for (Size i = 0; i < 11; ++i) {
1820 for (Size j = 0; j <= i; ++j) {
1821 if (std::fabs(boost::accumulators::covariance(v_eu[i][j]) - v_an[i][j]) > tol) {
1822 BOOST_ERROR("analytical covariance at ("
1823 << i << "," << j << ") (" << v_an[i][j]
1824 << ") is inconsistent with numerical "
1825 "value (Euler discretization, "
1826 << boost::accumulators::covariance(v_eu[i][j]) << "), error is "
1827 << v_an[i][j] - boost::accumulators::covariance(v_eu[i][j]) << " tolerance is " << tol);
1828 }
1829 if (std::fabs(boost::accumulators::covariance(v_eu2[i][j]) - v_an[i][j]) > tol) {
1830 BOOST_ERROR("analytical covariance at ("
1831 << i << "," << j << ") (" << v_an[i][j]
1832 << ") is inconsistent with numerical "
1833 "value (exact discretization, "
1834 << boost::accumulators::covariance(v_eu2[i][j]) << "), error is "
1835 << v_an[i][j] - boost::accumulators::covariance(v_eu2[i][j]) << " tolerance is " << tol);
1836 }
1837 }
1838 }
1839
1840} // testIrFxCrMoments
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [13/24]

BOOST_AUTO_TEST_CASE ( testIrFxCrCorrelationRecovery  )

Definition at line 1842 of file crossassetmodel.cpp.

1842 {
1843
1844 BOOST_TEST_MESSAGE("Test if random correlation input is recovered for "
1845 "small dt in ir-fx-cr model...");
1846
1847 class PseudoCurrency : public Currency {
1848 public:
1849 PseudoCurrency(const Size id) {
1850 std::ostringstream ln, sn;
1851 ln << "Dummy " << id;
1852 sn << "DUM " << id;
1853 data_ = QuantLib::ext::make_shared<Data>(ln.str(), sn.str(), id, sn.str(), "", 100, Rounding(), "%3% %1$.2f");
1854 }
1855 };
1856
1857 Real dt = 1.0E-6;
1858 Real tol = 1.0E-7;
1859
1860 // for ir-fx this fully specifies the correlation matrix
1861 // for new asset classes add other possible combinations as well
1862 Size currencies[] = { 1, 2, 3, 4, 5, 10, 20 };
1863 Size creditnames[] = { 0, 1, 5, 10 };
1864
1865 MersenneTwisterUniformRng mt(42);
1866
1867 Handle<YieldTermStructure> yts(QuantLib::ext::make_shared<FlatForward>(0, NullCalendar(), 0.01, Actual365Fixed()));
1868
1869 Handle<DefaultProbabilityTermStructure> hts(
1870 QuantLib::ext::make_shared<FlatHazardRate>(0, NullCalendar(), 0.01, Actual365Fixed()));
1871
1872 Handle<Quote> fxspot(QuantLib::ext::make_shared<SimpleQuote>(1.00));
1873
1874 Array notimes(0);
1875 Array fxsigma(1, 0.10);
1876
1877 for (Size ii = 0; ii < LENGTH(currencies); ++ii) {
1878 for (Size jj = 0; jj < LENGTH(creditnames); ++jj) {
1879
1880 std::vector<Currency> pseudoCcy;
1881 for (Size i = 0; i < currencies[ii]; ++i) {
1882 PseudoCurrency tmp(i);
1883 pseudoCcy.push_back(tmp);
1884 }
1885
1886 Size dim = 2 * currencies[ii] - 1 + creditnames[jj];
1887
1888 // generate random correlation matrix
1889 Matrix b(dim, dim);
1890 Size maxTries = 100;
1891 bool valid = true;
1892 do {
1893 Matrix a(dim, dim);
1894 for (Size i = 0; i < dim; ++i) {
1895 for (Size j = 0; j <= i; ++j) {
1896 a[i][j] = a[j][i] = mt.nextReal() - 0.5;
1897 }
1898 }
1899 b = a * transpose(a);
1900 for (Size i = 0; i < dim; ++i) {
1901 if (b[i][i] < 1E-5)
1902 valid = false;
1903 }
1904 } while (!valid && --maxTries > 0);
1905
1906 if (maxTries == 0) {
1907 BOOST_ERROR("could no generate random matrix");
1908 return;
1909 }
1910
1911 Matrix c(dim, dim);
1912 for (Size i = 0; i < dim; ++i) {
1913 for (Size j = 0; j <= i; ++j) {
1914 c[i][j] = c[j][i] = b[i][j] / std::sqrt(b[i][i] * b[j][j]);
1915 }
1916 }
1917
1918 // set up model
1919
1920 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
1921
1922 // IR
1923 for (Size i = 0; i < currencies[ii]; ++i) {
1924 parametrizations.push_back(
1925 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(pseudoCcy[i], yts, 0.01, 0.01));
1926 }
1927 // FX
1928 for (Size i = 0; i < currencies[ii] - 1; ++i) {
1929 parametrizations.push_back(QuantLib::ext::make_shared<FxBsPiecewiseConstantParametrization>(
1930 pseudoCcy[i + 1], fxspot, notimes, fxsigma));
1931 }
1932 // CR
1933 for (Size i = 0; i < creditnames[jj]; ++i) {
1934 parametrizations.push_back(
1935 QuantLib::ext::make_shared<CrLgm1fConstantParametrization>(pseudoCcy[0], hts, 0.01, 0.01));
1936 }
1937
1938 // get QuantLib::Error: negative eigenvalue(s) (-3.649315e-16) with SalvagingAlgorithm::None
1939 QuantLib::ext::shared_ptr<CrossAssetModel> modelExact =
1940 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, c, SalvagingAlgorithm::Spectral,
1941 IrModel::Measure::LGM, CrossAssetModel::Discretization::Exact);
1942 QuantLib::ext::shared_ptr<CrossAssetModel> modelEuler =
1943 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, c, SalvagingAlgorithm::Spectral,
1944 IrModel::Measure::LGM, CrossAssetModel::Discretization::Euler);
1945
1946 QuantLib::ext::shared_ptr<StochasticProcess> peuler = modelEuler->stateProcess();
1947 QuantLib::ext::shared_ptr<StochasticProcess> pexact = modelExact->stateProcess();
1948
1949 Matrix c1 = peuler->covariance(dt, peuler->initialValues(), dt);
1950 Matrix c2 = pexact->covariance(0.0, peuler->initialValues(), dt);
1951
1952 Matrix r1(dim, dim), r2(dim, dim);
1953
1954 for (Size i = 0; i < dim; ++i) {
1955 for (Size j = 0; j <= i; ++j) {
1956 // there are two state variables per credit name,
1957 Size subi = i < 2 * currencies[ii] - 1 ? 1 : 2;
1958 Size subj = j < 2 * currencies[ii] - 1 ? 1 : 2;
1959 for (Size k1 = 0; k1 < subi; ++k1) {
1960 for (Size k2 = 0; k2 < subj; ++k2) {
1961 Size i0 = i < 2 * currencies[ii] - 1
1962 ? i
1963 : 2 * currencies[ii] - 1 + 2 * (i - (2 * currencies[ii] - 1)) + k1;
1964 Size j0 = j < 2 * currencies[ii] - 1
1965 ? j
1966 : 2 * currencies[ii] - 1 + 2 * (j - (2 * currencies[ii] - 1)) + k2;
1967 r1[i][j] = r1[j][i] = c1[i0][j0] / std::sqrt(c1[i0][i0] * c1[j0][j0]);
1968 r2[i][j] = r2[j][i] = c2[i0][j0] / std::sqrt(c2[i0][i0] * c2[j0][j0]);
1969 if (std::fabs(r1[i][j] - c[i][j]) > tol) {
1970 BOOST_ERROR("failed to recover correlation matrix from "
1971 "Euler state process (i,j)=("
1972 << i << "," << j << "), (i0,j0)=(" << i0 << "," << j0
1973 << "), input correlation is " << c[i][j] << ", output is " << r1[i][j]
1974 << ", difference " << (c[i][j] - r1[i][j]) << ", tolerance " << tol
1975 << " test configuration is " << currencies[ii] << " currencies and "
1976 << creditnames[jj] << " credit names");
1977 }
1978 if (k1 == k2) {
1979 if (std::fabs(r2[i][j] - c[i][j]) > tol) {
1980 BOOST_ERROR("failed to recover correlation matrix "
1981 "from "
1982 "exact state process (i,j)=("
1983 << i << "," << j << "), (i0,j0)=(" << i0 << "," << j0
1984 << "), input correlation is " << c[i][j] << ", output is " << r2[i][j]
1985 << ", difference " << (c[i][j] - r2[i][j]) << ", tolerance " << tol
1986 << " test configuration is " << currencies[ii] << " currencies and "
1987 << creditnames[jj] << " credit names");
1988 }
1989 }
1990 }
1991 }
1992 }
1993 }
1994 } // for creditnames
1995 } // for currencies
1996
1997} // testIrFxCrCorrelationRecovery

◆ BOOST_DATA_TEST_CASE() [1/2]

BOOST_DATA_TEST_CASE ( testIrFxInfCrComMartingaleProperty  ,
bdata::make(infEurFlags) *bdata::make(infGbpFlags) *bdata::make(flatVolsFlags) *bdata::make(driftFreeState)  ,
infEurIsDk  ,
infGbpIsDk  ,
flatVols  ,
driftFreeState   
)

Definition at line 2286 of file crossassetmodel.cpp.

2288 {
2289
2290 BOOST_TEST_MESSAGE("Testing martingale property in ir-fx-inf-cr-com model for Euler and exact discretizations...");
2291 BOOST_TEST_MESSAGE("EUR inflation model is: " << (infEurIsDk ? "DK" : "JY"));
2292 BOOST_TEST_MESSAGE("GBP inflation model is: " << (infGbpIsDk ? "DK" : "JY"));
2293 BOOST_TEST_MESSAGE("Using " << (flatVols ? "" : "non-") << "flat volatilities.");
2294
2295 IrFxInfCrComModelTestData d(infEurIsDk, infGbpIsDk, flatVols, driftFreeState);
2296
2297 BOOST_TEST_MESSAGE("get exact state process");
2298 QuantLib::ext::shared_ptr<StochasticProcess> process1 = d.modelExact->stateProcess();
2299 BOOST_TEST_MESSAGE("get Euler state process");
2300 QuantLib::ext::shared_ptr<StochasticProcess> process2 = d.modelEuler->stateProcess();
2301
2302 Size n = 5000; // number of paths
2303 Size seed = 18; // rng seed
2304 Time T = 2.0; // maturity of payoff
2305 Time T2 = 20.0; // zerobond maturity
2306 Size steps = static_cast<Size>(T * 40); // number of steps taken (euler)
2307
2308 BOOST_TEST_MESSAGE("build sequence generators");
2309 LowDiscrepancy::rsg_type sg1 = LowDiscrepancy::make_sequence_generator(process1->factors() * 1, seed);
2310 LowDiscrepancy::rsg_type sg2 = LowDiscrepancy::make_sequence_generator(process2->factors() * steps, seed);
2311
2312 BOOST_TEST_MESSAGE("build multi path generator");
2313 TimeGrid grid1(T, 1);
2314 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process1)) {
2315 tmp->resetCache(grid1.size() - 1);
2316 }
2317 MultiPathGenerator<LowDiscrepancy::rsg_type> pg1(process1, grid1, sg1, false);
2318 TimeGrid grid2(T, steps);
2319 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process2)) {
2320 tmp->resetCache(grid2.size() - 1);
2321 }
2322 MultiPathGenerator<LowDiscrepancy::rsg_type> pg2(process2, grid2, sg2, false);
2323
2324 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > eurzb1, usdzb1, gbpzb1, infeur1, infgbp1,
2325 n1eur1, commodityA_1, commodityB_1;
2326 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > eurzb2, usdzb2, gbpzb2, infeur2, infgbp2,
2327 n1eur2, commodityA_2, commodityB_2;
2328
2329 for (Size j = 0; j < n; ++j) {
2330 Sample<MultiPath> path2 = pg2.next();
2331 Size l2 = path2.value[0].length() - 1;
2332 Sample<MultiPath> path1 = pg1.next();
2333 Size l1 = path1.value[0].length() - 1;
2334 Real zeur1 = path1.value[0][l1];
2335 Real zusd1 = path1.value[1][l1];
2336 Real zgbp1 = path1.value[2][l1];
2337 Real fxusd1 = std::exp(path1.value[3][l1]);
2338 Real fxgbp1 = std::exp(path1.value[4][l1]);
2339 Real infeurz1 = path1.value[5][l1];
2340 Real infeury1 = path1.value[6][l1];
2341 Real infgbpz1 = path1.value[7][l1];
2342 Real infgbpy1 = path1.value[8][l1];
2343 Real crzn11 = path1.value[9][l1];
2344 Real cryn11 = path1.value[10][l1];
2345 Real coma1 = path1.value[11][l1];
2346 Real comb1 = path1.value[12][l1];
2347 Real zeur2 = path2.value[0][l2];
2348 Real zusd2 = path2.value[1][l2];
2349 Real zgbp2 = path2.value[2][l2];
2350 Real fxusd2 = std::exp(path2.value[3][l2]);
2351 Real fxgbp2 = std::exp(path2.value[4][l2]);
2352 Real infeurz2 = path2.value[5][l2];
2353 Real infeury2 = path2.value[6][l2];
2354 Real infgbpz2 = path2.value[7][l2];
2355 Real infgbpy2 = path2.value[8][l2];
2356 Real crzn12 = path2.value[9][l2];
2357 Real cryn12 = path2.value[10][l2];
2358 Real coma2 = path2.value[11][l2];
2359 Real comb2 = path2.value[12][l2];
2360 // EUR zerobond
2361 eurzb1(d.modelExact->discountBond(0, T, T2, zeur1) / d.modelExact->numeraire(0, T, zeur1));
2362 // USD zerobond
2363 usdzb1(d.modelExact->discountBond(1, T, T2, zusd1) * fxusd1 / d.modelExact->numeraire(0, T, zeur1));
2364 // GBP zerobond
2365 gbpzb1(d.modelExact->discountBond(2, T, T2, zgbp1) * fxgbp1 / d.modelExact->numeraire(0, T, zeur1));
2366 // EUR CPI indexed bond
2367 bool indexIsInterpolated = true;
2368 if (infEurIsDk) {
2369 std::pair<Real, Real> sinfeur1 = d.modelExact->infdkI(0, T, T2, infeurz1, infeury1);
2370 infeur1(sinfeur1.first * sinfeur1.second * d.modelExact->discountBond(0, T, T2, zeur1) /
2371 d.modelExact->numeraire(0, T, zeur1));
2372 } else {
2373 infeur1(exp(infeury1) * inflationGrowth(d.modelExact, 0, T, T2, zeur1, infeurz1, indexIsInterpolated) *
2374 d.modelExact->discountBond(0, T, T2, zeur1) / d.modelExact->numeraire(0, T, zeur1));
2375 }
2376 // GBP CPI indexed bond
2377 if (infGbpIsDk) {
2378 std::pair<Real, Real> sinfgbp1 = d.modelExact->infdkI(1, T, T2, infgbpz1, infgbpy1);
2379 infgbp1(sinfgbp1.first * sinfgbp1.second * d.modelExact->discountBond(2, T, T2, zgbp1) * fxgbp1 /
2380 d.modelExact->numeraire(0, T, zeur1));
2381 } else {
2382 infgbp1(exp(infgbpy1) * inflationGrowth(d.modelExact, 1, T, T2, zgbp1, infgbpz1, indexIsInterpolated) *
2383 d.modelExact->discountBond(2, T, T2, zgbp1) * fxgbp1 / d.modelExact->numeraire(0, T, zeur1));
2384 }
2385 // EUR defaultable zerobond
2386 std::pair<Real, Real> sn11 = d.modelExact->crlgm1fS(0, 0, T, T2, crzn11, cryn11);
2387 n1eur1(sn11.first * sn11.second * d.modelExact->discountBond(0, T, T2, zeur1) / d.modelExact->numeraire(0, T, zeur1));
2388
2389 // EUR zerobond
2390 eurzb2(d.modelExact->discountBond(0, T, T2, zeur2) / d.modelExact->numeraire(0, T, zeur2));
2391 // USD zerobond
2392 usdzb2(d.modelExact->discountBond(1, T, T2, zusd2) * fxusd2 / d.modelExact->numeraire(0, T, zeur2));
2393 // GBP zerobond
2394 gbpzb2(d.modelExact->discountBond(2, T, T2, zgbp2) * fxgbp2 / d.modelExact->numeraire(0, T, zeur2));
2395 // EUR CPI indexed bond
2396 if (infEurIsDk) {
2397 std::pair<Real, Real> sinfeur2 = d.modelExact->infdkI(0, T, T2, infeurz2, infeury2);
2398 infeur2(sinfeur2.first * sinfeur2.second * d.modelExact->discountBond(0, T, T2, zeur2) /
2399 d.modelExact->numeraire(0, T, zeur2));
2400 } else {
2401 infeur2(exp(infeury2) * inflationGrowth(d.modelExact, 0, T, T2, zeur2, infeurz2, indexIsInterpolated) *
2402 d.modelExact->discountBond(0, T, T2, zeur2) / d.modelExact->numeraire(0, T, zeur2));
2403 }
2404 // GBP CPI indexed bond
2405 if (infGbpIsDk) {
2406 std::pair<Real, Real> sinfgbp2 = d.modelExact->infdkI(1, T, T2, infgbpz2, infgbpy2);
2407 infgbp2(sinfgbp2.first * sinfgbp2.second * d.modelExact->discountBond(2, T, T2, zgbp2) * fxgbp2 /
2408 d.modelExact->numeraire(0, T, zeur2));
2409 } else {
2410 infgbp2(exp(infgbpy2) * inflationGrowth(d.modelExact, 1, T, T2, zgbp2, infgbpz2, indexIsInterpolated) *
2411 d.modelExact->discountBond(2, T, T2, zgbp2) * fxgbp2 / d.modelExact->numeraire(0, T, zeur2));
2412 }
2413 // EUR defaultable zerobond
2414 std::pair<Real, Real> sn12 = d.modelExact->crlgm1fS(0, 0, T, T2, crzn12, cryn12);
2415 n1eur2(sn12.first * sn12.second * d.modelExact->discountBond(0, T, T2, zeur2) / d.modelExact->numeraire(0, T, zeur2));
2416 // commodity forward prices
2417 commodityA_1(d.comModelA->forwardPrice(T, T2, Array(1, coma1)));
2418 commodityB_1(d.comModelB->forwardPrice(T, T2, Array(1, comb1)));
2419 commodityA_2(d.comModelA->forwardPrice(T, T2, Array(1, coma2)));
2420 commodityB_2(d.comModelB->forwardPrice(T, T2, Array(1, comb2)));
2421 }
2422
2423 BOOST_TEST_MESSAGE("EXACT:");
2424 BOOST_TEST_MESSAGE("EUR zb = " << mean(eurzb1) << " +- " << error_of<tag::mean>(eurzb1) << " vs analytical "
2425 << d.eurYts->discount(T2));
2426 BOOST_TEST_MESSAGE("USD zb = " << mean(usdzb1) << " +- " << error_of<tag::mean>(usdzb1) << " vs analytical "
2427 << d.usdYts->discount(T2) * d.fxEurUsd->value());
2428 BOOST_TEST_MESSAGE("GBP zb = " << mean(gbpzb1) << " +- " << error_of<tag::mean>(gbpzb1) << " vs analytical "
2429 << d.gbpYts->discount(T2) * d.fxEurGbp->value());
2430 BOOST_TEST_MESSAGE("IDX zb EUR = " << mean(infeur1) << " +- " << error_of<tag::mean>(infeur1) << " vs analytical "
2431 << d.eurYts->discount(T2) *
2432 std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2));
2433 BOOST_TEST_MESSAGE("IDX zb GBP = " << mean(infgbp1) << " +- " << error_of<tag::mean>(infgbp1) << " vs analytical "
2434 << d.gbpYts->discount(T2) *
2435 std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) *
2436 d.fxEurGbp->value());
2437 BOOST_TEST_MESSAGE("N1 zb EUR = " << mean(n1eur1) << " +- " << error_of<tag::mean>(n1eur1) << " vs analytical "
2438 << d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2));
2439
2440 BOOST_TEST_MESSAGE("\nEULER:");
2441 BOOST_TEST_MESSAGE("EUR zb = " << mean(eurzb2) << " +- " << error_of<tag::mean>(eurzb2) << " vs analytical "
2442 << d.eurYts->discount(T2));
2443 BOOST_TEST_MESSAGE("USD zb = " << mean(usdzb2) << " +- " << error_of<tag::mean>(usdzb2) << " vs analytical "
2444 << d.usdYts->discount(T2) * d.fxEurUsd->value());
2445 BOOST_TEST_MESSAGE("GBP zb = " << mean(gbpzb2) << " +- " << error_of<tag::mean>(gbpzb2) << " vs analytical "
2446 << d.gbpYts->discount(T2) * d.fxEurGbp->value());
2447 BOOST_TEST_MESSAGE("IDX zb EUR = " << mean(infeur2) << " +- " << error_of<tag::mean>(infeur2) << " vs analytical "
2448 << d.eurYts->discount(T2) *
2449 std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2));
2450 BOOST_TEST_MESSAGE("IDX zb GBP = " << mean(infgbp2) << " +- " << error_of<tag::mean>(infgbp2) << " vs analytical "
2451 << d.gbpYts->discount(T2) *
2452 std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) *
2453 d.fxEurGbp->value());
2454 BOOST_TEST_MESSAGE("N1 zb EUR = " << mean(n1eur2) << " +- " << error_of<tag::mean>(n1eur2) << " vs analytical "
2455 << d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2));
2456
2457 // a bit higher than for plain zero bond , since we look at indexed zero
2458 // bonds, too
2459 Real tol1 = 5.0E-4; // EXACT
2460 Real tol2 = 14.0E-4; // EULER
2461
2462 Real ev = d.eurYts->discount(T2);
2463 if (std::abs(mean(eurzb1) - ev) > tol1)
2464 BOOST_TEST_ERROR("Martingale test failed for eurzb (exact discr.),"
2465 "expected "
2466 << ev << ", got " << mean(eurzb1) << ", tolerance " << tol1);
2467 ev = d.usdYts->discount(T2) * d.fxEurUsd->value();
2468 if (std::abs(mean(usdzb1) - ev) > tol1)
2469 BOOST_TEST_ERROR("Martingale test failed for eurzb (exact discr.),"
2470 "expected "
2471 << ev << ", got " << mean(usdzb1) << ", tolerance " << tol1);
2472 ev = d.gbpYts->discount(T2) * d.fxEurGbp->value();
2473 if (std::abs(mean(gbpzb1) - ev) > tol1)
2474 BOOST_TEST_ERROR("Martingale test failed for eurzb (exact discr.),"
2475 "expected "
2476 << ev << ", got " << mean(gbpzb1) << ", tolerance " << tol1);
2477 ev = d.eurYts->discount(T2) * std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2);
2478 if (std::abs(mean(infeur1) - ev) > tol1)
2479 BOOST_TEST_ERROR("Martingale test failed for idx eurzb (exact discr.),"
2480 "expected "
2481 << ev << ", got " << mean(infeur1) << ", tolerance " << tol1);
2482 ev = d.gbpYts->discount(T2) * std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) * d.fxEurGbp->value();
2483 if (std::abs(mean(infgbp1) - ev) > tol1)
2484 BOOST_TEST_ERROR("Martingale test failed for idx gbpzb (exact discr.),"
2485 "expected "
2486 << ev << ", got " << mean(infgbp1) << ", tolerance " << tol1);
2487 ev = d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2);
2488 if (std::abs(mean(n1eur1) - ev) > tol1)
2489 BOOST_TEST_ERROR("Martingale test failed for def eurzb (exact discr.),"
2490 "expected "
2491 << ev << ", got " << mean(n1eur1) << ", tolerance " << tol1);
2492
2493 ev = d.eurYts->discount(T2);
2494 if (std::abs(mean(eurzb2) - ev) > tol2)
2495 BOOST_TEST_ERROR("Martingale test failed for eurzb (Euler discr.),"
2496 "expected "
2497 << ev << ", got " << mean(eurzb2) << ", tolerance " << tol2);
2498 ev = d.usdYts->discount(T2) * d.fxEurUsd->value();
2499 if (std::abs(mean(usdzb2) - ev) > tol2)
2500 BOOST_TEST_ERROR("Martingale test failed for usdzb (Euler discr.),"
2501 "expected "
2502 << ev << ", got " << mean(usdzb2) << ", tolerance " << tol2 * error_of<tag::mean>(usdzb2));
2503 ev = d.gbpYts->discount(T2) * d.fxEurGbp->value();
2504 if (std::abs(mean(gbpzb2) - ev) > tol2)
2505 BOOST_TEST_ERROR("Martingale test failed for gbpzb (Euler discr.),"
2506 "expected "
2507 << ev << ", got " << mean(gbpzb2) << ", tolerance " << tol2);
2508 ev = d.eurYts->discount(T2) * std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2);
2509 if (std::abs(mean(infeur2) - ev) > tol2)
2510 BOOST_TEST_ERROR("Martingale test failed for idx eurzb (Euler discr.),"
2511 "expected "
2512 << ev << ", got " << mean(infeur2) << ", tolerance " << tol1);
2513 ev = d.gbpYts->discount(T2) * std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) * d.fxEurGbp->value();
2514 if (std::abs(mean(infgbp2) - ev) > tol2)
2515 BOOST_TEST_ERROR("Martingale test failed for idx gbpzb (Euler discr.),"
2516 "expected "
2517 << ev << ", got " << mean(infgbp2) << ", tolerance " << tol1);
2518 ev = d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2);
2519 if (std::abs(mean(n1eur2) - ev) > tol2)
2520 BOOST_TEST_ERROR("Martingale test failed for def eurzb (Euler discr.),"
2521 "expected "
2522 << ev << ", got " << mean(n1eur1) << ", tolerance " << tol1);
2523
2524 // commodity A forward prices
2525 Real tol;
2526 ev = d.comParametrizationA->priceCurve()->price(T2);
2527 tol = error_of<tag::mean>(commodityA_1);
2528 if (std::abs(mean(commodityA_1) - ev) > tol)
2529 BOOST_TEST_ERROR("Martingale test failed for commodity A (exact discr.),"
2530 "expected " << ev << ", got " << mean(commodityA_1) << " +- " << tol);
2531 tol = error_of<tag::mean>(commodityA_2);
2532 if (std::abs(mean(commodityA_2) - ev) > tol)
2533 BOOST_TEST_ERROR("Martingale test failed for commodity A (Euler discr.),"
2534 "expected " << ev << ", got " << mean(commodityA_2) << " +- " << tol);
2535
2536 // commodity B forward prices
2537 ev = d.comParametrizationB->priceCurve()->price(T2);
2538 tol = error_of<tag::mean>(commodityB_1);
2539 if (std::abs(mean(commodityB_1) - ev) > tol)
2540 BOOST_TEST_ERROR("Martingale test failed for commodity B (exact discr.),"
2541 "expected " << ev << ", got " << mean(commodityB_1) << " +- " << tol);
2542 tol = error_of<tag::mean>(commodityB_2);
2543 if (std::abs(mean(commodityB_2) - ev) > tol)
2544 BOOST_TEST_ERROR("Martingale test failed for commodity B (Euler discr.),"
2545 "expected " << ev << ", got " << mean(commodityB_2) << " +- " << tol);
2546
2547
2548} // testIrFxInfCrMartingaleProperty
QuantExt::inflationGrowth
Real inflationGrowth(const QuantLib::ext::shared_ptr< CrossAssetModel > &model, Size index, Time S, Time T, Real irState, Real rrState, bool indexIsInterpolated)
Definition: jyimpliedzeroinflationtermstructure.cpp:61
driftFreeState
vector< bool > driftFreeState
Definition: crossassetmodel.cpp:2284
+ Here is the call graph for this function:

◆ BOOST_DATA_TEST_CASE() [2/2]

BOOST_DATA_TEST_CASE ( testIrFxInfCrComMoments  ,
bdata::make(infEurFlags) *bdata::make(infGbpFlags) *bdata::make(flatVolsFlags) *bdata::make(driftFreeState)  ,
infEurIsDk  ,
infGbpIsDk  ,
flatVols  ,
driftFreeState   
)

Definition at line 2550 of file crossassetmodel.cpp.

2552 {
2553
2554 BOOST_TEST_MESSAGE("Testing analytic moments vs. Euler and exact discretization in ir-fx-inf-cr-com model...");
2555 BOOST_TEST_MESSAGE("EUR inflation model is: " << (infEurIsDk ? "DK" : "JY"));
2556 BOOST_TEST_MESSAGE("GBP inflation model is: " << (infGbpIsDk ? "DK" : "JY"));
2557 BOOST_TEST_MESSAGE("Using " << (flatVols ? "" : "non-") << "flat volatilities.");
2558
2559 IrFxInfCrComModelTestData d(infEurIsDk, infGbpIsDk, flatVols, driftFreeState);
2560
2561 Size n = d.modelExact->dimension();
2562
2563 QuantLib::ext::shared_ptr<StochasticProcess> p_exact = d.modelExact->stateProcess();
2564 QuantLib::ext::shared_ptr<StochasticProcess> p_euler = d.modelExact->stateProcess();
2565
2566 Real T = 2.0; // horizon at which we compare the moments
2567 Size steps = static_cast<Size>(T * 10); // number of simulation steps (Euler and exact)
2568 Size paths = 10000; // number of paths
2569
2570 Array e_an = p_exact->expectation(0.0, p_exact->initialValues(), T);
2571 Matrix v_an = p_exact->covariance(0.0, p_exact->initialValues(), T);
2572
2573 Size seed = 18;
2574 TimeGrid grid(T, steps);
2575
2576 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_euler)) {
2577 tmp->resetCache(grid.size() - 1);
2578 }
2579 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_exact)) {
2580 tmp->resetCache(grid.size() - 1);
2581 }
2582 MultiPathGeneratorSobolBrownianBridge pgen(p_euler, grid, SobolBrownianGenerator::Diagonal, seed,
2583 SobolRsg::JoeKuoD7);
2584 MultiPathGeneratorSobolBrownianBridge pgen2(p_exact, grid, SobolBrownianGenerator::Diagonal, seed,
2585 SobolRsg::JoeKuoD7);
2586
2587 vector<accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > > e_eu(n);
2588 vector<accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > > e_eu2(n);
2589 vector<vector<accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > > > v_eu(n);
2590 vector<vector<accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > > > v_eu2(n);
2591 for (Size i = 0; i < n; ++i) {
2592 v_eu[i] = vector<accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > >(n);
2593 v_eu2[i] = vector<accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > >(n);
2594 }
2595
2596 for (Size i = 0; i < paths; ++i) {
2597 Sample<MultiPath> path = pgen.next();
2598 Sample<MultiPath> path2 = pgen2.next();
2599 for (Size ii = 0; ii < n; ++ii) {
2600 Real cii = path.value[ii].back();
2601 Real cii2 = path2.value[ii].back();
2602 e_eu[ii](cii);
2603 e_eu2[ii](cii2);
2604 for (Size jj = 0; jj <= ii; ++jj) {
2605 Real cjj = path.value[jj].back();
2606 v_eu[ii][jj](cii, covariate1 = cjj);
2607 Real cjj2 = path2.value[jj].back();
2608 v_eu2[ii][jj](cii2, covariate1 = cjj2);
2609 }
2610 }
2611 }
2612
2613 for (Size i = 0; i < n; ++i) {
2614 BOOST_TEST_MESSAGE("E_" << i << " " << std::fixed << std::setprecision(12) <<
2615 e_an[i] << " " << mean(e_eu[i]) << " " << mean(e_eu2[i]));
2616 }
2617 BOOST_TEST_MESSAGE("==================");
2618
2619 BOOST_TEST_MESSAGE("one step analytical");
2620 for (Size i = 0; i < n; ++i) {
2621 std::ostringstream tmp;
2622 for (Size j = 0; j <= i; ++j) {
2623 tmp << v_an[i][j] << " ";
2624 }
2625 BOOST_TEST_MESSAGE(tmp.str());
2626 }
2627 BOOST_TEST_MESSAGE("==================");
2628
2629 BOOST_TEST_MESSAGE("euler");
2630 for (Size i = 0; i < n; ++i) {
2631 std::ostringstream tmp;
2632 for (Size j = 0; j <= i; ++j) {
2633 tmp << boost::accumulators::covariance(v_eu[i][j]) << " ";
2634 }
2635 BOOST_TEST_MESSAGE(tmp.str());
2636 }
2637 BOOST_TEST_MESSAGE("==================");
2638
2639 BOOST_TEST_MESSAGE("exact");
2640 for (Size i = 0; i < n; ++i) {
2641 std::ostringstream tmp;
2642 for (Size j = 0; j <= i; ++j) {
2643 tmp << boost::accumulators::covariance(v_eu2[i][j]) << " ";
2644 }
2645 BOOST_TEST_MESSAGE(tmp.str());
2646 }
2647 BOOST_TEST_MESSAGE("==================");
2648
2649 Real errTolLd[] = { 0.5E-4, 0.5E-4, 0.5E-4, 10.0E-4, 10.0E-4, 1E-4, 1E-4, 1E-4, 1E-4, 1E-4, 1E-4, 1E-4, 1E-4 };
2650
2651 for (Size i = 0; i < n; ++i) {
2652 // check expectation against analytical calculation (Euler)
2653 if (std::fabs(mean(e_eu[i]) - e_an[i]) > errTolLd[i]) {
2654 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
2655 << ") is inconsistent with numerical value (Euler "
2656 "discretization, "
2657 << mean(e_eu[i]) << "), error is "
2658 << e_an[i] - mean(e_eu[i]) << " tolerance is "
2659 << errTolLd[i]);
2660 }
2661 // check expectation against analytical calculation (exact disc)
2662 if (std::fabs(mean(e_eu2[i]) - e_an[i]) > errTolLd[i]) {
2663 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
2664 << ") is inconsistent with numerical value (exact "
2665 "discretization, "
2666 << mean(e_eu2[i]) << "), error is "
2667 << e_an[i] - mean(e_eu2[i]) << " tolerance is "
2668 << errTolLd[i]);
2669 }
2670 }
2671
2672 // as above, this is a bit rough compared to the more differentiated
2673 // test of the IR-FX model ...
2674 Real tol = 10.0E-4;
2675
2676 for (Size i = 0; i < n; ++i) {
2677 for (Size j = 0; j <= i; ++j) {
2678 if (std::fabs(boost::accumulators::covariance(v_eu[i][j]) - v_an[i][j]) > tol) {
2679 BOOST_ERROR("analytical covariance at ("
2680 << i << "," << j << ") (" << v_an[i][j]
2681 << ") is inconsistent with numerical "
2682 "value (Euler discretization, "
2683 << boost::accumulators::covariance(v_eu[i][j]) << "), error is "
2684 << v_an[i][j] - boost::accumulators::covariance(v_eu[i][j]) << " tolerance is " << tol);
2685 }
2686 if (std::fabs(boost::accumulators::covariance(v_eu2[i][j]) - v_an[i][j]) > tol) {
2687 BOOST_ERROR("analytical covariance at ("
2688 << i << "," << j << ") (" << v_an[i][j]
2689 << ") is inconsistent with numerical "
2690 "value (exact discretization, "
2691 << boost::accumulators::covariance(v_eu2[i][j]) << "), error is "
2692 << v_an[i][j] - boost::accumulators::covariance(v_eu2[i][j]) << " tolerance is " << tol);
2693 }
2694 }
2695 }
2696
2697} // testIrFxInfCrMoments
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [14/24]

BOOST_AUTO_TEST_CASE ( testIrFxInfCrEqMartingaleProperty  )

Definition at line 2910 of file crossassetmodel.cpp.

2910 {
2911
2912 BOOST_TEST_MESSAGE("Testing martingale property in ir-fx-inf-cr-eq model for "
2913 "Euler and exact discretizations...");
2914
2915 IrFxInfCrEqModelTestData d;
2916
2917 QuantLib::ext::shared_ptr<StochasticProcess> process1 = d.modelExact->stateProcess();
2918 QuantLib::ext::shared_ptr<StochasticProcess> process2 = d.modelEuler->stateProcess();
2919
2920 Size n = 50000; // number of paths
2921 Size seed = 18; // rng seed
2922 Time T = 2.0; // maturity of payoff
2923 Time T2 = 20.0; // zerobond maturity
2924 Size steps = static_cast<Size>(T * 24); // number of steps taken (euler)
2925
2926 // this can be made more accurate by using LowDiscrepancy instead
2927 // of PseudoRandom, but we use an error estimator for the check
2928 LowDiscrepancy::rsg_type sg1 = LowDiscrepancy::make_sequence_generator(process1->factors() * 1, seed);
2929 LowDiscrepancy::rsg_type sg2 = LowDiscrepancy::make_sequence_generator(process2->factors() * steps, seed);
2930
2931 TimeGrid grid1(T, 1);
2932 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process1)) {
2933 tmp->resetCache(grid1.size() - 1);
2934 }
2935 MultiPathGenerator<LowDiscrepancy::rsg_type> pg1(process1, grid1, sg1, false);
2936 TimeGrid grid2(T, steps);
2937 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process2)) {
2938 tmp->resetCache(grid2.size() - 1);
2939 }
2940 MultiPathGenerator<LowDiscrepancy::rsg_type> pg2(process2, grid2, sg2, false);
2941
2942 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > eurzb1, usdzb1, gbpzb1, infeur1, infgbp1,
2943 n1eur1, eqsp1, eqlh1;
2944 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > eurzb2, usdzb2, gbpzb2, infeur2, infgbp2,
2945 n1eur2, eqsp2, eqlh2;
2946
2947 for (Size j = 0; j < n; ++j) {
2948 Sample<MultiPath> path1 = pg1.next();
2949 Sample<MultiPath> path2 = pg2.next();
2950 Size l1 = path1.value[0].length() - 1;
2951 Size l2 = path2.value[0].length() - 1;
2952 Real zeur1 = path1.value[0][l1];
2953 Real zusd1 = path1.value[1][l1];
2954 Real zgbp1 = path1.value[2][l1];
2955 Real fxusd1 = std::exp(path1.value[3][l1]);
2956 Real fxgbp1 = std::exp(path1.value[4][l1]);
2957 Real infeurz1 = path1.value[5][l1];
2958 Real infeury1 = path1.value[6][l1];
2959 Real infgbpz1 = path1.value[7][l1];
2960 Real infgbpy1 = path1.value[8][l1];
2961 Real crzn11 = path1.value[9][l1];
2962 Real cryn11 = path1.value[10][l1];
2963 Real eq11 = path1.value[11][l1];
2964 Real eq21 = path1.value[12][l1];
2965 Real zeur2 = path2.value[0][l2];
2966 Real zusd2 = path2.value[1][l2];
2967 Real zgbp2 = path2.value[2][l2];
2968 Real fxusd2 = std::exp(path2.value[3][l2]);
2969 Real fxgbp2 = std::exp(path2.value[4][l2]);
2970 Real infeurz2 = path2.value[5][l2];
2971 Real infeury2 = path2.value[6][l2];
2972 Real infgbpz2 = path2.value[7][l2];
2973 Real infgbpy2 = path2.value[8][l2];
2974 Real crzn12 = path2.value[9][l2];
2975 Real cryn12 = path2.value[10][l2];
2976 Real eq12 = path2.value[11][l2];
2977 Real eq22 = path2.value[12][l2];
2978
2979 // EUR zerobond
2980 eurzb1(d.modelExact->discountBond(0, T, T2, zeur1) / d.modelExact->numeraire(0, T, zeur1));
2981 // USD zerobond
2982 usdzb1(d.modelExact->discountBond(1, T, T2, zusd1) * fxusd1 / d.modelExact->numeraire(0, T, zeur1));
2983 // GBP zerobond
2984 gbpzb1(d.modelExact->discountBond(2, T, T2, zgbp1) * fxgbp1 / d.modelExact->numeraire(0, T, zeur1));
2985 // EUR CPI indexed bond
2986 std::pair<Real, Real> sinfeur1 = d.modelExact->infdkI(0, T, T2, infeurz1, infeury1);
2987 infeur1(sinfeur1.first * sinfeur1.second * d.modelExact->discountBond(0, T, T2, zeur1) /
2988 d.modelExact->numeraire(0, T, zeur1));
2989 // GBP CPI indexed bond
2990 std::pair<Real, Real> sinfgbp1 = d.modelExact->infdkI(1, T, T2, infgbpz1, infgbpy1);
2991 infgbp1(sinfgbp1.first * sinfgbp1.second * d.modelExact->discountBond(2, T, T2, zgbp1) * fxgbp1 /
2992 d.modelExact->numeraire(0, T, zeur1));
2993 // EUR defaultable zerobond
2994 std::pair<Real, Real> sn11 = d.modelExact->crlgm1fS(0, 0, T, T2, crzn11, cryn11);
2995 n1eur1(sn11.first * sn11.second * d.modelExact->discountBond(0, T, T2, zeur1) / d.modelExact->numeraire(0, T, zeur1));
2996 // EQ
2997 eqsp1(std::exp(eq11) * fxusd1 / d.modelExact->numeraire(0, T, zeur1));
2998 eqlh1(std::exp(eq21) / d.modelExact->numeraire(0, T, zeur1));
2999
3000 // EUR zerobond
3001 eurzb2(d.modelExact->discountBond(0, T, T2, zeur2) / d.modelExact->numeraire(0, T, zeur2));
3002 // USD zerobond
3003 usdzb2(d.modelExact->discountBond(1, T, T2, zusd2) * fxusd2 / d.modelExact->numeraire(0, T, zeur2));
3004 // GBP zerobond
3005 gbpzb2(d.modelExact->discountBond(2, T, T2, zgbp2) * fxgbp2 / d.modelExact->numeraire(0, T, zeur2));
3006 // EUR CPI indexed bond
3007 std::pair<Real, Real> sinfeur2 = d.modelExact->infdkI(0, T, T2, infeurz2, infeury2);
3008 infeur2(sinfeur2.first * sinfeur2.second * d.modelExact->discountBond(0, T, T2, zeur2) /
3009 d.modelExact->numeraire(0, T, zeur2));
3010 // GBP CPI indexed bond
3011 std::pair<Real, Real> sinfgbp2 = d.modelExact->infdkI(1, T, T2, infgbpz2, infgbpy2);
3012 infgbp2(sinfgbp2.first * sinfgbp2.second * d.modelExact->discountBond(2, T, T2, zgbp2) * fxgbp2 /
3013 d.modelExact->numeraire(0, T, zeur2));
3014 // EUR defaultable zerobond
3015 std::pair<Real, Real> sn12 = d.modelExact->crlgm1fS(0, 0, T, T2, crzn12, cryn12);
3016 n1eur2(sn12.first * sn12.second * d.modelExact->discountBond(0, T, T2, zeur2) / d.modelExact->numeraire(0, T, zeur2));
3017 // EQ
3018 eqsp2(std::exp(eq12) * fxusd2 / d.modelExact->numeraire(0, T, zeur2));
3019 eqlh2(std::exp(eq22) / d.modelExact->numeraire(0, T, zeur2));
3020 }
3021
3022 BOOST_TEST_MESSAGE("EXACT:");
3023 BOOST_TEST_MESSAGE("EUR zb = " << mean(eurzb1) << " +- " << error_of<tag::mean>(eurzb1) << " vs analytical "
3024 << d.eurYts->discount(T2));
3025 BOOST_TEST_MESSAGE("USD zb = " << mean(usdzb1) << " +- " << error_of<tag::mean>(usdzb1) << " vs analytical "
3026 << d.usdYts->discount(T2) * d.fxEurUsd->value());
3027 BOOST_TEST_MESSAGE("GBP zb = " << mean(gbpzb1) << " +- " << error_of<tag::mean>(gbpzb1) << " vs analytical "
3028 << d.gbpYts->discount(T2) * d.fxEurGbp->value());
3029 BOOST_TEST_MESSAGE("IDX zb EUR = " << mean(infeur1) << " +- " << error_of<tag::mean>(infeur1) << " vs analytical "
3030 << d.eurYts->discount(T2) *
3031 std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2));
3032 BOOST_TEST_MESSAGE("IDX zb GBP = " << mean(infgbp1) << " +- " << error_of<tag::mean>(infgbp1) << " vs analytical "
3033 << d.gbpYts->discount(T2) *
3034 std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) *
3035 d.fxEurGbp->value());
3036 BOOST_TEST_MESSAGE("N1 zb EUR = " << mean(n1eur1) << " +- " << error_of<tag::mean>(n1eur1) << " vs analytical "
3037 << d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2));
3038 BOOST_TEST_MESSAGE("EQSP USD = " << mean(eqsp1) << " +- " << error_of<tag::mean>(eqsp1) << " vs analytical "
3039 << d.spSpotToday->value() * d.eqDivSp->discount(T) * d.fxEurUsd->value());
3040 BOOST_TEST_MESSAGE("EQLH EUR = " << mean(eqlh1) << " +- " << error_of<tag::mean>(eqlh1) << " vs analytical "
3041 << d.lhSpotToday->value() * d.eqDivLh->discount(T));
3042
3043 BOOST_TEST_MESSAGE("\nEULER:");
3044 BOOST_TEST_MESSAGE("EUR zb = " << mean(eurzb2) << " +- " << error_of<tag::mean>(eurzb2) << " vs analytical "
3045 << d.eurYts->discount(T2));
3046 BOOST_TEST_MESSAGE("USD zb = " << mean(usdzb2) << " +- " << error_of<tag::mean>(usdzb2) << " vs analytical "
3047 << d.usdYts->discount(T2) * d.fxEurUsd->value());
3048 BOOST_TEST_MESSAGE("GBP zb = " << mean(gbpzb2) << " +- " << error_of<tag::mean>(gbpzb2) << " vs analytical "
3049 << d.gbpYts->discount(T2) * d.fxEurGbp->value());
3050 BOOST_TEST_MESSAGE("IDX zb EUR = " << mean(infeur2) << " +- " << error_of<tag::mean>(infeur2) << " vs analytical "
3051 << d.eurYts->discount(T2) *
3052 std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2));
3053 BOOST_TEST_MESSAGE("IDX zb GBP = " << mean(infgbp2) << " +- " << error_of<tag::mean>(infgbp2) << " vs analytical "
3054 << d.gbpYts->discount(T2) *
3055 std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) *
3056 d.fxEurGbp->value());
3057 BOOST_TEST_MESSAGE("N1 zb EUR = " << mean(n1eur2) << " +- " << error_of<tag::mean>(n1eur2) << " vs analytical "
3058 << d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2));
3059 BOOST_TEST_MESSAGE("EQSP USD = " << mean(eqsp2) << " +- " << error_of<tag::mean>(eqsp2) << " vs analytical "
3060 << d.spSpotToday->value() * d.eqDivSp->discount(T) * d.fxEurUsd->value());
3061 BOOST_TEST_MESSAGE("EQLH EUR = " << mean(eqlh2) << " +- " << error_of<tag::mean>(eqlh2) << " vs analytical "
3062 << d.lhSpotToday->value() * d.eqDivLh->discount(T));
3063
3064 // a bit higher than for plain zero bond , since we look at indexed zero
3065 // bonds, too
3066 Real tol1 = 3.0E-4, tol1r = 0.001; // EXACT
3067 Real tol2 = 14.0E-4, tol2r = 0.01; // EULER
3068
3069 Real ev = d.eurYts->discount(T2);
3070 if (std::abs(mean(eurzb1) - ev) > tol1)
3071 BOOST_TEST_ERROR("Martingale test failed for eurzb (exact discr.),"
3072 "expected "
3073 << ev << ", got " << mean(eurzb1) << ", tolerance " << tol1);
3074 ev = d.usdYts->discount(T2) * d.fxEurUsd->value();
3075 if (std::abs(mean(usdzb1) - ev) > tol1)
3076 BOOST_TEST_ERROR("Martingale test failed for eurzb (exact discr.),"
3077 "expected "
3078 << ev << ", got " << mean(usdzb1) << ", tolerance " << tol1);
3079 ev = d.gbpYts->discount(T2) * d.fxEurGbp->value();
3080 if (std::abs(mean(gbpzb1) - ev) > tol1)
3081 BOOST_TEST_ERROR("Martingale test failed for eurzb (exact discr.),"
3082 "expected "
3083 << ev << ", got " << mean(gbpzb1) << ", tolerance " << tol1);
3084 ev = d.eurYts->discount(T2) * std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2);
3085 if (std::abs(mean(infeur1) - ev) > tol1)
3086 BOOST_TEST_ERROR("Martingale test failed for idx eurzb (exact discr.),"
3087 "expected "
3088 << ev << ", got " << mean(infeur1) << ", tolerance " << tol1);
3089 ev = d.gbpYts->discount(T2) * std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) * d.fxEurGbp->value();
3090 if (std::abs(mean(infgbp1) - ev) > tol1)
3091 BOOST_TEST_ERROR("Martingale test failed for idx gbpzb (exact discr.),"
3092 "expected "
3093 << ev << ", got " << mean(infgbp1) << ", tolerance " << tol1);
3094 ev = d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2);
3095 if (std::abs(mean(n1eur1) - ev) > tol1)
3096 BOOST_TEST_ERROR("Martingale test failed for def eurzb (exact discr.),"
3097 "expected "
3098 << ev << ", got " << mean(n1eur1) << ", tolerance " << tol1);
3099 ev = d.spSpotToday->value() * d.eqDivSp->discount(T) * d.fxEurUsd->value();
3100 if (std::abs(mean(eqsp1) - ev) / ev > tol1r)
3101 BOOST_TEST_ERROR("Martingale test failed for eq sp (exact discr.),"
3102 "expected "
3103 << ev << ", got " << mean(eqsp1) << ", rel tolerance " << tol1r);
3104 ev = d.lhSpotToday->value() * d.eqDivLh->discount(T);
3105 if (std::abs(mean(eqlh1) - ev) / ev > tol1r)
3106 BOOST_TEST_ERROR("Martingale test failed for eq lh (exact discr.),"
3107 "expected "
3108 << ev << ", got " << mean(eqlh1) << ", rel tolerance " << tol1r);
3109
3110 ev = d.eurYts->discount(T2);
3111 if (std::abs(mean(eurzb2) - ev) > tol2)
3112 BOOST_TEST_ERROR("Martingale test failed for eurzb (Euler discr.),"
3113 "expected "
3114 << ev << ", got " << mean(eurzb2) << ", tolerance " << tol2);
3115 ev = d.usdYts->discount(T2) * d.fxEurUsd->value();
3116 if (std::abs(mean(usdzb2) - ev) > tol2)
3117 BOOST_TEST_ERROR("Martingale test failed for usdzb (Euler discr.),"
3118 "expected "
3119 << ev << ", got " << mean(usdzb2) << ", tolerance " << tol2 * error_of<tag::mean>(usdzb2));
3120 ev = d.gbpYts->discount(T2) * d.fxEurGbp->value();
3121 if (std::abs(mean(gbpzb2) - ev) > tol2)
3122 BOOST_TEST_ERROR("Martingale test failed for gbpzb (Euler discr.),"
3123 "expected "
3124 << ev << ", got " << mean(gbpzb2) << ", tolerance " << tol2);
3125 ev = d.eurYts->discount(T2) * std::pow(1.0 + d.infEurTs->zeroRate(T2 - d.infLag), T2);
3126 if (std::abs(mean(infeur2) - ev) > tol2)
3127 BOOST_TEST_ERROR("Martingale test failed for idx eurzb (Euler discr.),"
3128 "expected "
3129 << ev << ", got " << mean(infeur2) << ", tolerance " << tol2);
3130 ev = d.gbpYts->discount(T2) * std::pow(1.0 + d.infGbpTs->zeroRate(T2 - d.infLag), T2) * d.fxEurGbp->value();
3131 if (std::abs(mean(infgbp2) - ev) > tol2)
3132 BOOST_TEST_ERROR("Martingale test failed for idx gbpzb (Euler discr.),"
3133 "expected "
3134 << ev << ", got " << mean(infgbp2) << ", tolerance " << tol2);
3135 ev = d.eurYts->discount(T2) * d.n1Ts->survivalProbability(T2);
3136 if (std::abs(mean(n1eur2) - ev) > tol2)
3137 BOOST_TEST_ERROR("Martingale test failed for def eurzb (Euler discr.),"
3138 "expected "
3139 << ev << ", got " << mean(n1eur1) << ", tolerance " << tol2);
3140 ev = d.spSpotToday->value() * d.eqDivSp->discount(T) * d.fxEurUsd->value();
3141 if (std::abs(mean(eqsp2) - ev) / ev > tol2r)
3142 BOOST_TEST_ERROR("Martingale test failed for eq sp (Euler discr.),"
3143 "expected "
3144 << ev << ", got " << mean(eqsp2) << ", rel tolerance " << tol2r);
3145 ev = d.lhSpotToday->value() * d.eqDivLh->discount(T);
3146 if (std::abs(mean(eqlh2) - ev) / ev > tol2r)
3147 BOOST_TEST_ERROR("Martingale test failed for eq lh (exact discr.),"
3148 "expected "
3149 << ev << ", got " << mean(eqlh2) << ", rel tolerance " << tol2r);
3150
3151} // testIrFxInfCrEqMartingaleProperty

◆ BOOST_AUTO_TEST_CASE() [15/24]

BOOST_AUTO_TEST_CASE ( testIrFxInfCrEqMoments  )

Definition at line 3153 of file crossassetmodel.cpp.

3153 {
3154
3155 BOOST_TEST_MESSAGE("Testing analytic moments vs. Euler and exact discretization "
3156 "in ir-fx-inf-cr-eq model...");
3157
3158 IrFxInfCrEqModelTestData d;
3159
3160 const Size n = 13; // d.modelExact->dimension();
3161
3162 QuantLib::ext::shared_ptr<StochasticProcess> p_exact = d.modelExact->stateProcess();
3163 QuantLib::ext::shared_ptr<StochasticProcess> p_euler = d.modelExact->stateProcess();
3164
3165 Real T = 2.0; // horizon at which we compare the moments
3166 Size steps = static_cast<Size>(T * 10); // number of simulation steps (Euler and exact)
3167 Size paths = 60000; // number of paths
3168
3169 Array e_an = p_exact->expectation(0.0, p_exact->initialValues(), T);
3170 Matrix v_an = p_exact->covariance(0.0, p_exact->initialValues(), T);
3171
3172 Size seed = 18;
3173 TimeGrid grid(T, steps);
3174
3175 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_euler)) {
3176 tmp->resetCache(grid.size() - 1);
3177 }
3178 MultiPathGeneratorSobolBrownianBridge pgen(p_euler, grid, SobolBrownianGenerator::Diagonal, seed,
3179 SobolRsg::JoeKuoD7);
3180 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_exact)) {
3181 tmp->resetCache(grid.size() - 1);
3182 }
3183 MultiPathGeneratorSobolBrownianBridge pgen2(p_exact, grid, SobolBrownianGenerator::Diagonal, seed,
3184 SobolRsg::JoeKuoD7);
3185
3186 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > e_eu[n], e_eu2[n];
3187 accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > v_eu[n][n], v_eu2[n][n];
3188
3189 for (Size i = 0; i < paths; ++i) {
3190 Sample<MultiPath> path = pgen.next();
3191 Sample<MultiPath> path2 = pgen2.next();
3192 for (Size ii = 0; ii < n; ++ii) {
3193 Real cii = path.value[ii].back();
3194 Real cii2 = path2.value[ii].back();
3195 e_eu[ii](cii);
3196 e_eu2[ii](cii2);
3197 for (Size jj = 0; jj <= ii; ++jj) {
3198 Real cjj = path.value[jj].back();
3199 v_eu[ii][jj](cii, covariate1 = cjj);
3200 Real cjj2 = path2.value[jj].back();
3201 v_eu2[ii][jj](cii2, covariate1 = cjj2);
3202 }
3203 }
3204 }
3205
3206 for (Size i = 0; i < n; ++i) {
3207 BOOST_TEST_MESSAGE("E_" << i << " " << e_an[i] << " " << mean(e_eu[i]) << " " << mean(e_eu2[i]));
3208 }
3209 BOOST_TEST_MESSAGE("==================");
3210
3211 BOOST_TEST_MESSAGE("one step analytical");
3212 for (Size i = 0; i < n; ++i) {
3213 std::ostringstream tmp;
3214 for (Size j = 0; j <= i; ++j) {
3215 tmp << v_an[i][j] << " ";
3216 }
3217 BOOST_TEST_MESSAGE(tmp.str());
3218 }
3219 BOOST_TEST_MESSAGE("==================");
3220
3221 BOOST_TEST_MESSAGE("euler");
3222 for (Size i = 0; i < n; ++i) {
3223 std::ostringstream tmp;
3224 for (Size j = 0; j <= i; ++j) {
3225 tmp << boost::accumulators::covariance(v_eu[i][j]) << " ";
3226 }
3227 BOOST_TEST_MESSAGE(tmp.str());
3228 }
3229 BOOST_TEST_MESSAGE("==================");
3230
3231 BOOST_TEST_MESSAGE("exact");
3232 for (Size i = 0; i < n; ++i) {
3233 std::ostringstream tmp;
3234 for (Size j = 0; j <= i; ++j) {
3235 tmp << boost::accumulators::covariance(v_eu2[i][j]) << " ";
3236 }
3237 BOOST_TEST_MESSAGE(tmp.str());
3238 }
3239 BOOST_TEST_MESSAGE("==================");
3240
3241 Real errTolLd[] = { 0.5E-4, 0.5E-4, 0.5E-4, 10.0E-4, 10.0E-4, 0.7E-4, 0.7E-4,
3242 0.7E-4, 0.7E-4, 0.7E-4, 0.7E-4, 10.0E-4, 10.0E-4 };
3243 Real errTolLdEuler[] = { 0.5E-4, 0.5E-4, 0.5E-4, 10.0E-4, 10.0E-4, 0.7E-4, 0.7E-4,
3244 0.7E-4, 0.7E-4, 0.7E-4, 0.7E-4, 40.0E-4, 40.0E-4 };
3245
3246 for (Size i = 0; i < n; ++i) {
3247 // check expectation against analytical calculation (Euler)
3248 if (std::fabs(mean(e_eu[i]) - e_an[i]) > errTolLdEuler[i]) {
3249 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
3250 << ") is inconsistent with numerical value (Euler "
3251 "discretization, "
3252 << mean(e_eu[i]) << "), error is "
3253 << e_an[i] - mean(e_eu[i]) << " tolerance is "
3254 << errTolLdEuler[i]);
3255 }
3256 // check expectation against analytical calculation (exact disc)
3257 if (std::fabs(mean(e_eu2[i]) - e_an[i]) > errTolLd[i]) {
3258 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
3259 << ") is inconsistent with numerical value (exact "
3260 "discretization, "
3261 << mean(e_eu2[i]) << "), error is "
3262 << e_an[i] - mean(e_eu2[i]) << " tolerance is "
3263 << errTolLd[i]);
3264 }
3265 }
3266
3267 // as above, this is a bit rough compared to the more differentiated
3268 // test of the IR-FX model ...
3269 Real tol = 10.0E-4, tolEuler = 65.0E-4;
3270
3271 for (Size i = 0; i < n; ++i) {
3272 for (Size j = 0; j <= i; ++j) {
3273 if (std::fabs(boost::accumulators::covariance(v_eu[i][j]) - v_an[i][j]) > tolEuler) {
3274 BOOST_ERROR("analytical covariance at ("
3275 << i << "," << j << ") (" << v_an[i][j]
3276 << ") is inconsistent with numerical "
3277 "value (Euler discretization, "
3278 << boost::accumulators::covariance(v_eu[i][j]) << "), error is "
3279 << v_an[i][j] - boost::accumulators::covariance(v_eu[i][j]) << " tolerance is "
3280 << tolEuler);
3281 }
3282 if (std::fabs(boost::accumulators::covariance(v_eu2[i][j]) - v_an[i][j]) > tol) {
3283 BOOST_ERROR("analytical covariance at ("
3284 << i << "," << j << ") (" << v_an[i][j]
3285 << ") is inconsistent with numerical "
3286 "value (exact discretization, "
3287 << boost::accumulators::covariance(v_eu2[i][j]) << "), error is "
3288 << v_an[i][j] - boost::accumulators::covariance(v_eu2[i][j]) << " tolerance is " << tol);
3289 }
3290 }
3291 }
3292
3293} // testIrFxInfCrEqMoments
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [16/24]

BOOST_AUTO_TEST_CASE ( testEqLgm5fPayouts  )

Definition at line 3459 of file crossassetmodel.cpp.

3459 {
3460
3461 BOOST_TEST_MESSAGE("Testing pricing of equity payouts under domestic "
3462 "measure in CrossAsset LGM model...");
3463
3464 IrFxEqModelTestData d;
3465 Settings::instance().evaluationDate() = d.referenceDate;
3466
3467 QuantLib::ext::shared_ptr<StochasticProcess> process = d.ccLgmExact->stateProcess();
3468 QuantLib::ext::shared_ptr<StochasticProcess> process2 = d.ccLgmEuler->stateProcess();
3469
3470 // path generation
3471
3472 Size n = 500000; // number of paths
3473 Size seed = 121; // seed
3474 // maturity of test payoffs
3475 Time T = 5.0;
3476 // take large steps, but not only one (since we are testing)
3477 Size steps = static_cast<Size>(T * 2.0);
3478 Size steps_euler = static_cast<Size>(T * 52.0);
3479 TimeGrid grid(T, steps);
3480 TimeGrid grid_euler(T, steps_euler);
3481 PseudoRandom::rsg_type sg2 = PseudoRandom::make_sequence_generator(steps, seed);
3482
3483 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process)) {
3484 tmp->resetCache(grid.size() - 1);
3485 }
3486 MultiPathGeneratorMersenneTwister pg(process, grid, seed, false);
3487 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process2)) {
3488 tmp->resetCache(grid_euler.size() - 1);
3489 }
3490 MultiPathGeneratorMersenneTwister pg2(process2, grid_euler, seed, false);
3491
3492 // test
3493 // 1 LH (EUR) forward under CrossAsset numeraire vs. analytic pricing engine
3494 // 2 SP (USD) forward (converted to EUR) under CrossAsset numeraire vs. analytic pricing engine
3495 // 3 LH (EUR) EQ option under CrossAsset numeraire vs. analytic pricing engine
3496 // 4 SP (USD) EQ option under CrossAsset numeraire vs. analytic pricing engine
3497
3498 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > stat1, stat2, stat3a, stat3b, stat4a, stat4b;
3499
3500 Real strikeLh = 12.7;
3501 Real strikeSp = 2150.0;
3502
3503 // same for paths2 since shared time grid
3504 for (Size j = 0; j < n; ++j) {
3505 Sample<MultiPath> path = pg.next();
3506 // Sample<MultiPath> path = pg2.next();
3507 Size l = path.value[0].length() - 1;
3508 Real eurusdfx = std::exp(path.value[2][l]);
3509 Real zeur = path.value[0][l];
3510 Real eqsp = std::exp(path.value[3][l]);
3511 Real eqlh = std::exp(path.value[4][l]);
3512 Real ccnum = d.ccLgmExact->numeraire(0, T, zeur);
3513
3514 // 1 LH forward settled at T deflated with numeraire
3515 Real lhFwd = eqlh - strikeLh;
3516 stat1(lhFwd / ccnum);
3517
3518 // 2 SP forward settled at T (converted to base) deflated with numeraire
3519 Real spFwd = eurusdfx * (eqsp - strikeSp);
3520 stat2(spFwd / ccnum);
3521
3522 // 3 LH option exercise at T deflated with numeraire
3523 Real lhCall = std::max(lhFwd, 0.0);
3524 Real lhPut = std::max(-1.0 * lhFwd, 0.0);
3525 stat3a(lhCall / ccnum);
3526 stat3b(lhPut / ccnum);
3527
3528 // 4 SP option exercise at T (converted to base) deflated with numeraire
3529 Real spCall = std::max(spFwd, 0.0);
3530 Real spPut = std::max(-1.0 * spFwd, 0.0);
3531 stat4a(spCall / ccnum);
3532 stat4b(spPut / ccnum);
3533 }
3534
3535 Date tradeMaturity = d.referenceDate + 5 * 365;
3536
3537 QuantLib::ext::shared_ptr<EquityForward> lhFwdTrade =
3538 QuantLib::ext::make_shared<EquityForward>("LH", EURCurrency(), Position::Long, 1, tradeMaturity, strikeLh);
3539 QuantLib::ext::shared_ptr<EquityForward> spFwdTrade =
3540 QuantLib::ext::make_shared<EquityForward>("SP", USDCurrency(), Position::Long, 1, tradeMaturity, strikeSp);
3541
3542 QuantLib::ext::shared_ptr<VanillaOption> lhCall =
3543 QuantLib::ext::make_shared<VanillaOption>(QuantLib::ext::make_shared<PlainVanillaPayoff>(Option::Call, strikeLh),
3544 QuantLib::ext::make_shared<EuropeanExercise>(d.referenceDate + 5 * 365));
3545 QuantLib::ext::shared_ptr<VanillaOption> lhPut =
3546 QuantLib::ext::make_shared<VanillaOption>(QuantLib::ext::make_shared<PlainVanillaPayoff>(Option::Put, strikeLh),
3547 QuantLib::ext::make_shared<EuropeanExercise>(d.referenceDate + 5 * 365));
3548 QuantLib::ext::shared_ptr<VanillaOption> spCall =
3549 QuantLib::ext::make_shared<VanillaOption>(QuantLib::ext::make_shared<PlainVanillaPayoff>(Option::Call, strikeSp),
3550 QuantLib::ext::make_shared<EuropeanExercise>(d.referenceDate + 5 * 365));
3551 QuantLib::ext::shared_ptr<VanillaOption> spPut =
3552 QuantLib::ext::make_shared<VanillaOption>(QuantLib::ext::make_shared<PlainVanillaPayoff>(Option::Put, strikeSp),
3553 QuantLib::ext::make_shared<EuropeanExercise>(d.referenceDate + 5 * 365));
3554
3555 QuantLib::ext::shared_ptr<DiscountingEquityForwardEngine> lhFwdEngine =
3556 QuantLib::ext::make_shared<DiscountingEquityForwardEngine>(d.eurYts, d.eqDivLh, d.lhSpotToday, d.eurYts);
3557 QuantLib::ext::shared_ptr<DiscountingEquityForwardEngine> spFwdEngine =
3558 QuantLib::ext::make_shared<DiscountingEquityForwardEngine>(d.usdYts, d.eqDivSp, d.spSpotToday, d.usdYts);
3559
3560 lhFwdTrade->setPricingEngine(lhFwdEngine);
3561 spFwdTrade->setPricingEngine(spFwdEngine);
3562
3563 QuantLib::ext::shared_ptr<AnalyticXAssetLgmEquityOptionEngine> spEqOptionEngine =
3564 QuantLib::ext::make_shared<AnalyticXAssetLgmEquityOptionEngine>(
3565 d.ccLgmExact, d.eqSpIdx, d.ccLgmExact->ccyIndex(d.ccLgmExact->eqbs(d.eqSpIdx)->currency()));
3566 QuantLib::ext::shared_ptr<AnalyticXAssetLgmEquityOptionEngine> lhEqOptionEngine =
3567 QuantLib::ext::make_shared<AnalyticXAssetLgmEquityOptionEngine>(
3568 d.ccLgmExact, d.eqLhIdx, d.ccLgmExact->ccyIndex(d.ccLgmExact->eqbs(d.eqLhIdx)->currency()));
3569
3570 lhCall->setPricingEngine(lhEqOptionEngine);
3571 lhPut->setPricingEngine(lhEqOptionEngine);
3572 spCall->setPricingEngine(spEqOptionEngine);
3573 spPut->setPricingEngine(spEqOptionEngine);
3574
3575 Real npv1 = mean(stat1);
3576 Real error1 = error_of<tag::mean>(stat1);
3577 Real expected1 = lhFwdTrade->NPV();
3578
3579 Real npv2 = mean(stat2);
3580 Real error2 = error_of<tag::mean>(stat2);
3581 Real expected2 = d.usdEurSpotToday->value() * spFwdTrade->NPV();
3582
3583 Real npv3a = mean(stat3a);
3584 Real error3a = error_of<tag::mean>(stat3a);
3585 Real expected3a = lhCall->NPV();
3586 Real npv3b = mean(stat3b);
3587 Real error3b = error_of<tag::mean>(stat3b);
3588 Real expected3b = lhPut->NPV();
3589
3590 Real npv4a = mean(stat4a);
3591 Real error4a = error_of<tag::mean>(stat4a);
3592 Real expected4a = d.usdEurSpotToday->value() * spCall->NPV();
3593 Real npv4b = mean(stat4b);
3594 Real error4b = error_of<tag::mean>(stat4b);
3595 Real expected4b = d.usdEurSpotToday->value() * spPut->NPV();
3596
3597 Real tolErrEst = 1.5; // allow absolute diffs to be within 1.5 standard errors
3598 // TODO : IS THIS STRICT ENOUGH?
3599 BOOST_CHECK_LE(std::fabs(npv1 - expected1), tolErrEst * error1);
3600 BOOST_CHECK_LE(std::fabs(npv2 - expected2), tolErrEst * error2);
3601 BOOST_CHECK_LE(std::fabs(npv3a - expected3a), tolErrEst * error3a);
3602 BOOST_CHECK_LE(std::fabs(npv3b - expected3b), tolErrEst * error3b);
3603 BOOST_CHECK_LE(std::fabs(npv4a - expected4a), tolErrEst * error4a);
3604 BOOST_CHECK_LE(std::fabs(npv4b - expected4b), tolErrEst * error4b);
3605
3606} // testEqLgm5fPayouts
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [17/24]

BOOST_AUTO_TEST_CASE ( testEqLgm5fCalibration  )

Definition at line 3608 of file crossassetmodel.cpp.

3608 {
3609
3610 BOOST_TEST_MESSAGE("Testing EQ calibration of IR-FX-EQ LGM 5F model...");
3611
3612 IrFxEqModelTestData d;
3613 Settings::instance().evaluationDate() = d.referenceDate;
3614
3615 // calibration baskets
3616 std::vector<QuantLib::ext::shared_ptr<BlackCalibrationHelper> > basketSp, basketLh;
3617
3618 for (Size i = 0; i < d.volstepdatesEqSp.size(); ++i) {
3619 Date tmp = i < d.volstepdatesEqSp.size() ? d.volstepdatesEqSp[i] : d.volstepdatesEqSp.back() + 365;
3620 basketSp.push_back(QuantLib::ext::make_shared<FxEqOptionHelper>(
3621 tmp, Null<Real>(), d.spSpotToday, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.20)), d.usdYts, d.eqDivSp,
3622 BlackCalibrationHelper::RelativePriceError));
3623 }
3624 for (Size i = 0; i < d.volstepdatesEqLh.size(); ++i) {
3625 Date tmp = i < d.volstepdatesEqLh.size() ? d.volstepdatesEqLh[i] : d.volstepdatesEqLh.back() + 365;
3626 basketLh.push_back(QuantLib::ext::make_shared<FxEqOptionHelper>(
3627 tmp, Null<Real>(), d.lhSpotToday, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(0.20)), d.eurYts, d.eqDivLh,
3628 BlackCalibrationHelper::RelativePriceError));
3629 }
3630
3631 // pricing engines
3632 QuantLib::ext::shared_ptr<AnalyticXAssetLgmEquityOptionEngine> spEqOptionEngine =
3633 QuantLib::ext::make_shared<AnalyticXAssetLgmEquityOptionEngine>(
3634 d.ccLgmExact, d.eqSpIdx, d.ccLgmExact->ccyIndex(d.ccLgmExact->eqbs(d.eqSpIdx)->currency()));
3635 QuantLib::ext::shared_ptr<AnalyticXAssetLgmEquityOptionEngine> lhEqOptionEngine =
3636 QuantLib::ext::make_shared<AnalyticXAssetLgmEquityOptionEngine>(
3637 d.ccLgmExact, d.eqLhIdx, d.ccLgmExact->ccyIndex(d.ccLgmExact->eqbs(d.eqLhIdx)->currency()));
3638
3639 // assign engines to calibration instruments
3640 for (Size i = 0; i < basketSp.size(); ++i) {
3641 basketSp[i]->setPricingEngine(spEqOptionEngine);
3642 }
3643 for (Size i = 0; i < basketLh.size(); ++i) {
3644 basketLh[i]->setPricingEngine(lhEqOptionEngine);
3645 }
3646
3647 // calibrate the model
3648 LevenbergMarquardt lm(1E-8, 1E-8, 1E-8);
3649 EndCriteria ec(1000, 500, 1E-8, 1E-8, 1E-8);
3650
3651 d.ccLgmExact->calibrateBsVolatilitiesIterative(CrossAssetModel::AssetType::EQ, d.eqSpIdx, basketSp, lm, ec);
3652 d.ccLgmExact->calibrateBsVolatilitiesIterative(CrossAssetModel::AssetType::EQ, d.eqLhIdx, basketLh, lm, ec);
3653
3654 // check the results
3655 Real tol = 1E-6;
3656
3657 for (Size i = 0; i < basketSp.size(); ++i) {
3658 Real model = basketSp[i]->modelValue();
3659 Real market = basketSp[i]->marketValue();
3660 if (std::abs((model - market) / market) > tol)
3661 BOOST_ERROR("calibration failed for instrument #"
3662 << i << " in SP basket, model value is " << model << " market value is " << market
3663 << " relative error " << std::abs((model - market) / market) << " tolerance " << tol);
3664 }
3665 for (Size i = 0; i < basketLh.size(); ++i) {
3666 Real model = basketLh[i]->modelValue();
3667 Real market = basketLh[i]->marketValue();
3668 if (std::abs((model - market) / market) > tol)
3669 BOOST_ERROR("calibration failed for instrument #"
3670 << i << " in LH basket, model value is " << model << " market value is " << market
3671 << " relative error " << std::abs((model - market) / market) << " tolerance " << tol);
3672 }
3673} // testEqLgm5fCalibration

◆ BOOST_AUTO_TEST_CASE() [18/24]

BOOST_AUTO_TEST_CASE ( testEqLgm5fMoments  )

Definition at line 3675 of file crossassetmodel.cpp.

3675 {
3676
3677 // TODO : REVIEW TOLERANCES
3678 BOOST_TEST_MESSAGE("Testing analytic moments vs. Euler and exact discretization "
3679 "in IR-FX-EQ LGM 5F model...");
3680
3681 IrFxEqModelTestData d;
3682 Settings::instance().evaluationDate() = d.referenceDate;
3683
3684 QuantLib::ext::shared_ptr<StochasticProcess> p_exact = d.ccLgmExact->stateProcess();
3685 QuantLib::ext::shared_ptr<StochasticProcess> p_euler = d.ccLgmEuler->stateProcess();
3686
3687 Real T = 2.0; // horizon at which we compare the moments
3688 Size steps_euler = static_cast<Size>(T * 50.0); // number of simulation steps
3689 Size steps_exact = 1;
3690 Size paths = 25000; // number of paths
3691
3692 Array e_an = p_exact->expectation(0.0, p_exact->initialValues(), T);
3693 Matrix v_an = p_exact->covariance(0.0, p_exact->initialValues(), T);
3694 Matrix v_an_eu = p_euler->covariance(0.0, p_euler->initialValues(), T);
3695
3696 TimeGrid grid_euler(T, steps_euler);
3697 TimeGrid grid_exact(T, steps_exact);
3698
3699 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_euler)) {
3700 tmp->resetCache(grid_euler.size() - 1);
3701 }
3702 MultiPathGeneratorSobolBrownianBridge pgen(p_euler, grid_euler);
3703
3704 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_exact)) {
3705 tmp->resetCache(grid_exact.size() - 1);
3706 }
3707 MultiPathGeneratorSobolBrownianBridge pgen2(p_exact, grid_exact);
3708
3709 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > e_eu[5], e_eu2[5];
3710 accumulator_set<double, stats<tag::covariance<double, tag::covariate1> > > v_eu[5][5], v_eu2[5][5];
3711
3712 for (Size i = 0; i < paths; ++i) {
3713 Sample<MultiPath> path = pgen.next();
3714 Sample<MultiPath> path2 = pgen2.next();
3715 for (Size ii = 0; ii < 5; ++ii) {
3716 Real cii = path.value[ii].back();
3717 Real cii2 = path2.value[ii].back();
3718 e_eu[ii](cii);
3719 e_eu2[ii](cii2);
3720 for (Size jj = 0; jj <= ii; ++jj) {
3721 Real cjj = path.value[jj].back();
3722 v_eu[ii][jj](cii, covariate1 = cjj);
3723 Real cjj2 = path2.value[jj].back();
3724 v_eu2[ii][jj](cii2, covariate1 = cjj2);
3725 }
3726 }
3727 }
3728
3729 Real errTol[] = { 0.2E-4, 0.2E-4, 10.0E-4, 10.0E-4, 10.0E-4 };
3730
3731 for (Size i = 0; i < 5; ++i) {
3732 // check expectation against analytical calculation (Euler)
3733 if (std::fabs(mean(e_eu[i]) - e_an[i]) > errTol[i]) {
3734 BOOST_ERROR("analytical expectation for component #"
3735 << i << " (" << e_an[i]
3736 << ") is inconsistent with numerical value (Euler "
3737 "discretization, "
3738 << mean(e_eu[i]) << "), error is " << e_an[i] - mean(e_eu[i]) << " tolerance is " << errTol[i]);
3739 }
3740 // check expectation against analytical calculation (exact disc)
3741 if (std::fabs(mean(e_eu2[i]) - e_an[i]) > errTol[i]) {
3742 BOOST_ERROR("analytical expectation for component #" << i << " (" << e_an[i]
3743 << ") is inconsistent with numerical value (Exact "
3744 "discretization, "
3745 << mean(e_eu2[i]) << "), error is "
3746 << e_an[i] - mean(e_eu2[i]) << " tolerance is "
3747 << errTol[i]);
3748 }
3749 }
3750
3751 // we have to deal with different natures of volatility
3752 // for ir (normal) and fx (ln) so different error
3753 // tolerances apply
3754 Real tollNormal = 0.1E-4; // ir-ir
3755 Real tolMixed = 0.25E-4; // ir-fx, ir-eq
3756 Real tolLn = 8.0E-4; // fx-fx, fx-eq
3757 Real tolEq = 12.0E-4; // eq-eq (to account for higher eq vols)
3758 Real tol; // set below
3759
3760 for (Size i = 0; i < 5; ++i) {
3761 for (Size j = 0; j <= i; ++j) {
3762 if (i < 2) {
3763 tol = tollNormal;
3764 } else if ((i >= 3) && (j >= 3)) {
3765 tol = tolEq;
3766 } else {
3767 if (j < 2) {
3768 tol = tolMixed;
3769 } else {
3770 tol = tolLn;
3771 }
3772 }
3773 if (std::fabs(boost::accumulators::covariance(v_eu[i][j]) - v_an[i][j]) > tol) {
3774 BOOST_ERROR("analytical covariance at ("
3775 << i << "," << j << ") (" << v_an[i][j]
3776 << ") is inconsistent with numerical "
3777 "value (Euler discretization, "
3778 << boost::accumulators::covariance(v_eu[i][j]) << "), error is "
3779 << v_an[i][j] - boost::accumulators::covariance(v_eu[i][j]) << " tolerance is " << tol);
3780 }
3781 if (std::fabs(boost::accumulators::covariance(v_eu2[i][j]) - v_an[i][j]) > tol) {
3782 BOOST_ERROR("analytical covariance at ("
3783 << i << "," << j << ") (" << v_an[i][j]
3784 << ") is inconsistent with numerical "
3785 "value (Exact discretization, "
3786 << boost::accumulators::covariance(v_eu2[i][j]) << "), error is "
3787 << v_an[i][j] - boost::accumulators::covariance(v_eu2[i][j]) << " tolerance is " << tol);
3788 }
3789 }
3790 }
3791
3792 BOOST_TEST_MESSAGE("Testing correlation matrix recovery in presence of equity simulation");
3793
3794 // reset caching so that we can retrieve further info from the processes
3795 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_euler)) {
3796 tmp->resetCache(0);
3797 }
3798 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(p_exact)) {
3799 tmp->resetCache(grid_euler.size() - 1);
3800 }
3801
3802
3803 Matrix corr_input = d.ccLgmExact->correlation();
3804 BOOST_CHECK(corr_input.rows() == corr_input.columns());
3805 Size dim = corr_input.rows();
3806 BOOST_CHECK(corr_input.rows() == 5);
3807 Matrix r1(dim, dim), r2(dim, dim);
3808 Real dt = 1.0E-6;
3809 Real tol_corr = 1.0E-7;
3810 Matrix v_an_dt = p_exact->covariance(0.0, p_exact->initialValues(), dt);
3811 Matrix v_an_eu_dt = p_euler->covariance(0.0, p_euler->initialValues(), dt);
3812 BOOST_CHECK(v_an_dt.rows() == v_an_eu_dt.rows());
3813 BOOST_CHECK(v_an_dt.columns() == v_an_eu_dt.columns());
3814 BOOST_CHECK(corr_input.rows() == v_an_dt.rows());
3815 BOOST_CHECK(corr_input.columns() == corr_input.columns());
3816
3817 for (Size i = 0; i < dim; ++i) {
3818 for (Size j = 0; j <= i; ++j) {
3819 r1[i][j] = r1[j][i] = v_an_dt[i][j] / std::sqrt(v_an_dt[i][i] * v_an_dt[j][j]);
3820 r2[i][j] = r2[j][i] = v_an_eu_dt[i][j] / std::sqrt(v_an_eu_dt[i][i] * v_an_eu_dt[j][j]);
3821 BOOST_CHECK_MESSAGE(std::fabs(r1[i][j] - corr_input[i][j]) < tol_corr,
3822 "failed to recover correlation matrix from "
3823 "exact state process (i,j)=("
3824 << i << "," << j << "), input correlation is " << corr_input[i][j] << ", output is "
3825 << r1[i][j] << ", difference " << (corr_input[i][j] - r1[i][j]) << ", tolerance "
3826 << tol_corr);
3827 BOOST_CHECK_MESSAGE(std::fabs(r2[i][j] - corr_input[i][j]) < tol_corr,
3828 "failed to recover correlation matrix from "
3829 "Euler state process (i,j)=("
3830 << i << "," << j << "), input correlation is " << corr_input[i][j] << ", output is "
3831 << r2[i][j] << ", difference " << (corr_input[i][j] - r2[i][j]) << ", tolerance "
3832 << tol_corr);
3833 }
3834 }
3835
3836 for (Size i = 0; i < 5; ++i) {
3837 Real meu = mean(e_eu[i]);
3838 Real s_meu = error_of<tag::mean>(e_eu[i]);
3839 BOOST_TEST_MESSAGE(i << ";EULER;" << e_an[i] << ";" << meu << ";" << s_meu);
3840 Real meu2 = mean(e_eu2[i]);
3841 Real s_meu2 = error_of<tag::mean>(e_eu2[i]);
3842 BOOST_TEST_MESSAGE(i << ";EXACT;" << e_an[i] << ";" << meu2 << ";" << s_meu2);
3843 }
3844 for (Size i = 0; i < 5; ++i) {
3845 for (Size j = 0; j <= i; ++j) {
3846 Real cov = boost::accumulators::covariance(v_eu[i][j]);
3847 BOOST_TEST_MESSAGE(i << ";" << j << ";"
3848 << "EULER;" << v_an[i][j] << ";" << cov);
3849 Real cov2 = boost::accumulators::covariance(v_eu2[i][j]);
3850 BOOST_TEST_MESSAGE(i << ";" << j << ";"
3851 << "EXACT;" << v_an[i][j] << ";" << cov2);
3852 }
3853 }
3854} // testEqLgm5fMoments
+ Here is the call graph for this function:

◆ BOOST_AUTO_TEST_CASE() [19/24]

BOOST_AUTO_TEST_CASE ( testCorrelationRecovery  )

Definition at line 3856 of file crossassetmodel.cpp.

3856 {
3857
3858 BOOST_TEST_MESSAGE("Test if random correlation input is recovered for "
3859 "small dt in Ccy LGM model...");
3860
3861 class PseudoCurrency : public Currency {
3862 public:
3863 PseudoCurrency(const Size id) {
3864 std::ostringstream ln, sn;
3865 ln << "Dummy " << id;
3866 sn << "DUM " << id;
3867 data_ = QuantLib::ext::make_shared<Data>(ln.str(), sn.str(), id, sn.str(), "", 100, Rounding(), "%3% %1$.2f");
3868 }
3869 };
3870
3871 Real dt = 1.0E-6;
3872 Real tol = 1.0E-7;
3873
3874 // for ir-fx this fully specifies the correlation matrix
3875 // for new asset classes add other possible combinations as well
3876 Size currencies[] = { 1, 2, 3, 4, 5, 10, 20 };
3877
3878 MersenneTwisterUniformRng mt(42);
3879
3880 Handle<YieldTermStructure> yts(QuantLib::ext::make_shared<FlatForward>(0, NullCalendar(), 0.01, Actual365Fixed()));
3881
3882 Handle<Quote> fxspot(QuantLib::ext::make_shared<SimpleQuote>(1.00));
3883
3884 Array notimes(0);
3885 Array fxsigma(1, 0.10);
3886
3887 for (Size ii = 0; ii < LENGTH(currencies); ++ii) {
3888
3889 std::vector<Currency> pseudoCcy;
3890 for (Size i = 0; i < currencies[ii]; ++i) {
3891 PseudoCurrency tmp(i);
3892 pseudoCcy.push_back(tmp);
3893 }
3894
3895 Size dim = 2 * currencies[ii] - 1;
3896
3897 // generate random correlation matrix
3898 Matrix b(dim, dim);
3899 Size maxTries = 100;
3900 bool valid = true;
3901 do {
3902 Matrix a(dim, dim);
3903 for (Size i = 0; i < dim; ++i) {
3904 for (Size j = 0; j <= i; ++j) {
3905 a[i][j] = a[j][i] = mt.nextReal() - 0.5;
3906 }
3907 }
3908 b = a * transpose(a);
3909 for (Size i = 0; i < dim; ++i) {
3910 if (b[i][i] < 1E-5)
3911 valid = false;
3912 }
3913 } while (!valid && --maxTries > 0);
3914
3915 if (maxTries == 0) {
3916 BOOST_ERROR("could no generate random matrix");
3917 return;
3918 }
3919
3920 Matrix c(dim, dim);
3921 for (Size i = 0; i < dim; ++i) {
3922 for (Size j = 0; j <= i; ++j) {
3923 c[i][j] = c[j][i] = b[i][j] / std::sqrt(b[i][i] * b[j][j]);
3924 }
3925 }
3926
3927 // set up model
3928
3929 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
3930
3931 // IR
3932 for (Size i = 0; i < currencies[ii]; ++i) {
3933 parametrizations.push_back(
3934 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(pseudoCcy[i], yts, 0.01, 0.01));
3935 }
3936 // FX
3937 for (Size i = 0; i < currencies[ii] - 1; ++i) {
3938 parametrizations.push_back(
3939 QuantLib::ext::make_shared<FxBsPiecewiseConstantParametrization>(pseudoCcy[i + 1], fxspot, notimes, fxsigma));
3940 }
3941
3942 QuantLib::ext::shared_ptr<CrossAssetModel> modelExact =
3943 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, c, SalvagingAlgorithm::None, IrModel::Measure::LGM,
3944 CrossAssetModel::Discretization::Exact);
3945 QuantLib::ext::shared_ptr<CrossAssetModel> modelEuler =
3946 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, c, SalvagingAlgorithm::None, IrModel::Measure::LGM,
3947 CrossAssetModel::Discretization::Euler);
3948
3949 QuantLib::ext::shared_ptr<StochasticProcess> peuler = modelExact->stateProcess();
3950 QuantLib::ext::shared_ptr<StochasticProcess> pexact = modelEuler->stateProcess();
3951
3952 Matrix c1 = peuler->covariance(dt, peuler->initialValues(), dt);
3953 Matrix c2 = pexact->covariance(0.0, peuler->initialValues(), dt);
3954
3955 Matrix r1(dim, dim), r2(dim, dim);
3956
3957 for (Size i = 0; i < dim; ++i) {
3958 for (Size j = 0; j <= i; ++j) {
3959 r1[i][j] = r1[j][i] = c1[i][j] / std::sqrt(c1[i][i] * c1[j][j]);
3960 r2[i][j] = r2[j][i] = c2[i][j] / std::sqrt(c2[i][i] * c2[j][j]);
3961 if (std::fabs(r1[i][j] - c[i][j]) > tol) {
3962 BOOST_ERROR("failed to recover correlation matrix from "
3963 "Euler state process (i,j)=("
3964 << i << "," << j << "), input correlation is " << c[i][j] << ", output is " << r1[i][j]
3965 << ", difference " << (c[i][j] - r1[i][j]) << ", tolerance " << tol);
3966 }
3967 if (std::fabs(r2[i][j] - c[i][j]) > tol) {
3968 BOOST_ERROR("failed to recover correlation matrix from "
3969 "exact state process (i,j)=("
3970 << i << "," << j << "), input correlation is " << c[i][j] << ", output is " << r2[i][j]
3971 << ", difference " << (c[i][j] - r2[i][j]) << ", tolerance " << tol);
3972 }
3973 }
3974 }
3975
3976 } // for currencies
3977
3978} // test correlation recovery

◆ BOOST_AUTO_TEST_CASE() [20/24]

BOOST_AUTO_TEST_CASE ( testIrFxInfCrCorrelationRecovery  )

Definition at line 3980 of file crossassetmodel.cpp.

3980 {
3981
3982 BOOST_TEST_MESSAGE("Test if random correlation input is recovered for "
3983 "small dt in ir-fx-inf-cr model...");
3984
3985 SavedSettings backup;
3986 Settings::instance().evaluationDate() = Date(30, July, 2015);
3987
3988 class PseudoCurrency : public Currency {
3989 public:
3990 PseudoCurrency(const Size id) {
3991 std::ostringstream ln, sn;
3992 ln << "Dummy " << id;
3993 sn << "DUM " << id;
3994 data_ = QuantLib::ext::make_shared<Data>(ln.str(), sn.str(), id, sn.str(), "", 100, Rounding(), "%3% %1$.2f");
3995 }
3996 };
3997
3998 Real dt = 1.0E-6;
3999 Real tol = 1.0E-7;
4000
4001 // for ir-fx this fully specifies the correlation matrix
4002 // for new asset classes add other possible combinations as well
4003 Size currencies[] = { 1, 2, 3, 4, 5, 10, 20 };
4004 Size cpiindexes[] = { 0, 1, 10 };
4005 Size creditnames[] = { 0, 1, 5 };
4006
4007 MersenneTwisterUniformRng mt(42);
4008
4009 Handle<YieldTermStructure> yts(QuantLib::ext::make_shared<FlatForward>(0, NullCalendar(), 0.01, Actual365Fixed()));
4010
4011 std::vector<Date> infDates;
4012 std::vector<Real> infRates;
4013 infDates.push_back(Date(30, April, 2015));
4014 infDates.push_back(Date(30, July, 2015));
4015 infRates.push_back(0.01);
4016 infRates.push_back(0.01);
4017 Handle<ZeroInflationTermStructure> its(
4018 QuantLib::ext::make_shared<ZeroInflationCurve>(Settings::instance().evaluationDate(), NullCalendar(), Actual365Fixed(),
4019 3 * Months, Monthly, infDates, infRates));
4020
4021 Handle<DefaultProbabilityTermStructure> hts(
4022 QuantLib::ext::make_shared<FlatHazardRate>(0, NullCalendar(), 0.01, Actual365Fixed()));
4023
4024 Handle<Quote> fxspot(QuantLib::ext::make_shared<SimpleQuote>(1.00));
4025
4026 Array notimes(0);
4027 Array fxsigma(1, 0.10);
4028
4029 for (Size ii = 0; ii < LENGTH(currencies); ++ii) {
4030 for (Size kk = 0; kk < LENGTH(cpiindexes); ++kk) {
4031 for (Size jj = 0; jj < LENGTH(creditnames); ++jj) {
4032
4033 std::vector<Currency> pseudoCcy;
4034 for (Size i = 0; i < currencies[ii]; ++i) {
4035 PseudoCurrency tmp(i);
4036 pseudoCcy.push_back(tmp);
4037 }
4038
4039 Size dim = 2 * currencies[ii] - 1 + cpiindexes[kk] + creditnames[jj];
4040
4041 // generate random correlation matrix
4042 Matrix b(dim, dim);
4043 Size maxTries = 100;
4044 bool valid = true;
4045 do {
4046 Matrix a(dim, dim);
4047 for (Size i = 0; i < dim; ++i) {
4048 for (Size j = 0; j <= i; ++j) {
4049 a[i][j] = a[j][i] = mt.nextReal() - 0.5;
4050 }
4051 }
4052 b = a * transpose(a);
4053 for (Size i = 0; i < dim; ++i) {
4054 if (b[i][i] < 1E-5)
4055 valid = false;
4056 }
4057 } while (!valid && --maxTries > 0);
4058
4059 if (maxTries == 0) {
4060 BOOST_ERROR("could no generate random matrix");
4061 return;
4062 }
4063
4064 Matrix c(dim, dim);
4065 for (Size i = 0; i < dim; ++i) {
4066 for (Size j = 0; j <= i; ++j) {
4067 c[i][j] = c[j][i] = b[i][j] / std::sqrt(b[i][i] * b[j][j]);
4068 }
4069 }
4070
4071 // set up model
4072
4073 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
4074
4075 // IR
4076 for (Size i = 0; i < currencies[ii]; ++i) {
4077 parametrizations.push_back(
4078 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(pseudoCcy[i], yts, 0.01, 0.01));
4079 }
4080 // FX
4081 for (Size i = 0; i < currencies[ii] - 1; ++i) {
4082 parametrizations.push_back(QuantLib::ext::make_shared<FxBsPiecewiseConstantParametrization>(
4083 pseudoCcy[i + 1], fxspot, notimes, fxsigma));
4084 }
4085 // INF
4086 for (Size i = 0; i < cpiindexes[kk]; ++i) {
4087 parametrizations.push_back(
4088 QuantLib::ext::make_shared<InfDkConstantParametrization>(pseudoCcy[0], its, 0.01, 0.01));
4089 }
4090 // CR
4091 for (Size i = 0; i < creditnames[jj]; ++i) {
4092 parametrizations.push_back(
4093 QuantLib::ext::make_shared<CrLgm1fConstantParametrization>(pseudoCcy[0], hts, 0.01, 0.01));
4094 }
4095
4096 // get QuantLib::Error: negative eigenvalue(s) (-3.649315e-16) with SalvagingAlgorithm::None
4097 QuantLib::ext::shared_ptr<CrossAssetModel> modelEuler =
4098 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, c, SalvagingAlgorithm::Spectral,
4099 IrModel::Measure::LGM, CrossAssetModel::Discretization::Euler);
4100 QuantLib::ext::shared_ptr<CrossAssetModel> modelExact =
4101 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, c, SalvagingAlgorithm::Spectral,
4102 IrModel::Measure::LGM, CrossAssetModel::Discretization::Exact);
4103
4104 QuantLib::ext::shared_ptr<StochasticProcess> peuler = modelEuler->stateProcess();
4105 QuantLib::ext::shared_ptr<StochasticProcess> pexact = modelExact->stateProcess();
4106
4107 Matrix c1 = peuler->covariance(dt, peuler->initialValues(), dt);
4108 Matrix c2 = pexact->covariance(0.0, peuler->initialValues(), dt);
4109
4110 Matrix r1(dim, dim), r2(dim, dim);
4111
4112 for (Size i = 0; i < dim; ++i) {
4113 for (Size j = 0; j <= i; ++j) {
4114 // there are two state variables per credit name,
4115 // and per inflation index
4116 Size subi = i < 2 * currencies[ii] - 1 ? 1 : 2;
4117 Size subj = j < 2 * currencies[ii] - 1 ? 1 : 2;
4118 for (Size k1 = 0; k1 < subi; ++k1) {
4119 for (Size k2 = 0; k2 < subj; ++k2) {
4120 Size i0 = i < 2 * currencies[ii] - 1
4121 ? i
4122 : 2 * currencies[ii] - 1 + 2 * (i - (2 * currencies[ii] - 1)) + k1;
4123 Size j0 = j < 2 * currencies[ii] - 1
4124 ? j
4125 : 2 * currencies[ii] - 1 + 2 * (j - (2 * currencies[ii] - 1)) + k2;
4126 r1[i][j] = r1[j][i] = c1[i0][j0] / std::sqrt(c1[i0][i0] * c1[j0][j0]);
4127 r2[i][j] = r2[j][i] = c2[i0][j0] / std::sqrt(c2[i0][i0] * c2[j0][j0]);
4128 if (std::fabs(r1[i][j] - c[i][j]) > tol) {
4129 BOOST_ERROR("failed to recover correlation matrix "
4130 "from "
4131 "Euler state process (i,j)=("
4132 << i << "," << j << "), (i0,j0)=(" << i0 << "," << j0
4133 << "), input correlation is " << c[i][j] << ", output is " << r1[i][j]
4134 << ", difference " << (c[i][j] - r1[i][j]) << ", tolerance " << tol
4135 << " test configuration is " << currencies[ii] << " currencies and "
4136 << cpiindexes[kk] << " cpi indexes and " << creditnames[jj]
4137 << " credit names");
4138 }
4139 if (k1 == k2) {
4140 if (std::fabs(r2[i][j] - c[i][j]) > tol) {
4141 BOOST_ERROR("failed to recover correlation "
4142 "matrix "
4143 "from "
4144 "exact state process (i,j)=("
4145 << i << "," << j << "), (i0,j0)=(" << i0 << "," << j0
4146 << "), input correlation is " << c[i][j] << ", output is "
4147 << r2[i][j] << ", difference " << (c[i][j] - r2[i][j])
4148 << ", tolerance " << tol << " test configuration is "
4149 << currencies[ii] << " currencies and " << cpiindexes[kk]
4150 << " cpi indexes and " << creditnames[jj] << " credit names");
4151 }
4152 }
4153 }
4154 }
4155 }
4156 }
4157 } // for creditnames
4158 } // for cpiindexes
4159 } // for currencies
4160} // testIrFxInfCrCorrelationRecovery

◆ BOOST_AUTO_TEST_CASE() [21/24]

BOOST_AUTO_TEST_CASE ( testIrFxInfCrEqCorrelationRecovery  )

Definition at line 4162 of file crossassetmodel.cpp.

4162 {
4163
4164 BOOST_TEST_MESSAGE("Test if random correlation input is recovered for "
4165 "small dt in ir-fx-inf-cr-eq model...");
4166
4167 SavedSettings backup;
4168 Settings::instance().evaluationDate() = Date(30, July, 2015);
4169
4170 class PseudoCurrency : public Currency {
4171 public:
4172 PseudoCurrency(const Size id) {
4173 std::ostringstream ln, sn;
4174 ln << "Dummy " << id;
4175 sn << "DUM " << id;
4176 data_ = QuantLib::ext::make_shared<Data>(ln.str(), sn.str(), id, sn.str(), "", 100, Rounding(), "%3% %1$.2f");
4177 }
4178 };
4179
4180 Real dt = 1.0E-6;
4181 Real tol = 1.0E-7;
4182
4183 // for ir-fx this fully specifies the correlation matrix
4184 // for new asset classes add other possible combinations as well
4185 Size currencies[] = { 1, 2, 3, 4, 5 };
4186 Size cpiindexes[] = { 0, 1, 10 };
4187 Size creditnames[] = { 0, 1, 5 };
4188 Size eqs[] = { 0, 1, 5 };
4189
4190 MersenneTwisterUniformRng mt(42);
4191
4192 Handle<YieldTermStructure> yts(QuantLib::ext::make_shared<FlatForward>(0, NullCalendar(), 0.01, Actual365Fixed()));
4193
4194 std::vector<Date> infDates;
4195 std::vector<Real> infRates;
4196 infDates.push_back(Date(30, April, 2015));
4197 infDates.push_back(Date(30, July, 2015));
4198 infRates.push_back(0.01);
4199 infRates.push_back(0.01);
4200 Handle<ZeroInflationTermStructure> its(
4201 QuantLib::ext::make_shared<ZeroInflationCurve>(Settings::instance().evaluationDate(), NullCalendar(), Actual365Fixed(),
4202 3 * Months, Monthly, infDates, infRates));
4203
4204 Handle<DefaultProbabilityTermStructure> hts(
4205 QuantLib::ext::make_shared<FlatHazardRate>(0, NullCalendar(), 0.01, Actual365Fixed()));
4206
4207 Handle<Quote> fxspot(QuantLib::ext::make_shared<SimpleQuote>(1.00)), eqspot(QuantLib::ext::make_shared<SimpleQuote>(1.00));
4208
4209 Array notimes(0);
4210 Array fxsigma(1, 0.10), eqsigma(1, 0.10);
4211
4212 for (Size ii = 0; ii < LENGTH(currencies); ++ii) {
4213 for (Size kk = 0; kk < LENGTH(cpiindexes); ++kk) {
4214 for (Size jj = 0; jj < LENGTH(creditnames); ++jj) {
4215 for (Size ee = 0; ee < LENGTH(eqs); ++ee) {
4216
4217 std::vector<Currency> pseudoCcy;
4218 for (Size i = 0; i < currencies[ii]; ++i) {
4219 PseudoCurrency tmp(i);
4220 pseudoCcy.push_back(tmp);
4221 }
4222
4223 Size dim = 2 * currencies[ii] - 1 + cpiindexes[kk] + creditnames[jj] + eqs[ee];
4224
4225 // generate random correlation matrix
4226 Matrix b(dim, dim);
4227 Size maxTries = 100;
4228 bool valid = true;
4229 do {
4230 Matrix a(dim, dim);
4231 for (Size i = 0; i < dim; ++i) {
4232 for (Size j = 0; j <= i; ++j) {
4233 a[i][j] = a[j][i] = mt.nextReal() - 0.5;
4234 }
4235 }
4236 b = a * transpose(a);
4237 for (Size i = 0; i < dim; ++i) {
4238 if (b[i][i] < 1E-5)
4239 valid = false;
4240 }
4241 } while (!valid && --maxTries > 0);
4242
4243 if (maxTries == 0) {
4244 BOOST_ERROR("could no generate random matrix");
4245 return;
4246 }
4247
4248 Matrix c(dim, dim);
4249 for (Size i = 0; i < dim; ++i) {
4250 for (Size j = 0; j <= i; ++j) {
4251 c[i][j] = c[j][i] = b[i][j] / std::sqrt(b[i][i] * b[j][j]);
4252 }
4253 }
4254
4255 // set up model
4256
4257 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
4258
4259 // IR
4260 for (Size i = 0; i < currencies[ii]; ++i) {
4261 parametrizations.push_back(
4262 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(pseudoCcy[i], yts, 0.01, 0.01));
4263 }
4264 // FX
4265 for (Size i = 0; i < currencies[ii] - 1; ++i) {
4266 parametrizations.push_back(QuantLib::ext::make_shared<FxBsPiecewiseConstantParametrization>(
4267 pseudoCcy[i + 1], fxspot, notimes, fxsigma));
4268 }
4269 // INF
4270 for (Size i = 0; i < cpiindexes[kk]; ++i) {
4271 parametrizations.push_back(
4272 QuantLib::ext::make_shared<InfDkConstantParametrization>(pseudoCcy[0], its, 0.01, 0.01));
4273 }
4274 // CR
4275 for (Size i = 0; i < creditnames[jj]; ++i) {
4276 parametrizations.push_back(
4277 QuantLib::ext::make_shared<CrLgm1fConstantParametrization>(pseudoCcy[0], hts, 0.01, 0.01));
4278 }
4279 // EQ
4280 for (Size i = 0; i < eqs[ee]; ++i) {
4281 parametrizations.push_back(QuantLib::ext::make_shared<EqBsPiecewiseConstantParametrization>(
4282 pseudoCcy[0], "dummy", eqspot, fxspot, notimes, eqsigma, yts, yts));
4283 }
4284
4285 // get QuantLib::Error: negative eigenvalue(s) (-3.649315e-16) with SalvagingAlgorithm::None
4286 QuantLib::ext::shared_ptr<CrossAssetModel> modelEuler = QuantLib::ext::make_shared<CrossAssetModel>(
4287 parametrizations, c, SalvagingAlgorithm::Spectral, IrModel::Measure::LGM,
4288 CrossAssetModel::Discretization::Euler);
4289 QuantLib::ext::shared_ptr<CrossAssetModel> modelExact = QuantLib::ext::make_shared<CrossAssetModel>(
4290 parametrizations, c, SalvagingAlgorithm::Spectral, IrModel::Measure::LGM,
4291 CrossAssetModel::Discretization::Exact);
4292
4293 QuantLib::ext::shared_ptr<StochasticProcess> peuler = modelEuler->stateProcess();
4294 QuantLib::ext::shared_ptr<StochasticProcess> pexact = modelExact->stateProcess();
4295
4296 Matrix c1 = peuler->covariance(dt, peuler->initialValues(), dt);
4297 Matrix c2 = pexact->covariance(0.0, peuler->initialValues(), dt);
4298
4299 Matrix r1(dim, dim), r2(dim, dim);
4300
4301 Size sizeIrFx = 2 * currencies[ii] - 1;
4302
4303 for (Size i = 0; i < dim; ++i) {
4304 for (Size j = 0; j <= i; ++j) {
4305 // there are two state variables per credit name,
4306 // and per inflation index
4307 Size subi = i < sizeIrFx || i >= sizeIrFx + cpiindexes[kk] + creditnames[jj] ? 1 : 2;
4308 Size subj = j < sizeIrFx || i >= sizeIrFx + cpiindexes[kk] + creditnames[jj] ? 1 : 2;
4309 for (Size k1 = 0; k1 < subi; ++k1) {
4310 for (Size k2 = 0; k2 < subj; ++k2) {
4311 Size i0 = i < sizeIrFx
4312 ? i
4313 : (i < sizeIrFx + cpiindexes[kk] + creditnames[jj]
4314 ? sizeIrFx + 2 * (i - sizeIrFx) + k1
4315 : sizeIrFx + 2 * cpiindexes[kk] + 2 * creditnames[jj] +
4316 (i - sizeIrFx - cpiindexes[kk] - creditnames[jj]));
4317 Size j0 = j < sizeIrFx
4318 ? j
4319 : (j < sizeIrFx + cpiindexes[kk] + creditnames[jj]
4320 ? sizeIrFx + 2 * (j - sizeIrFx) + k2
4321 : sizeIrFx + 2 * cpiindexes[kk] + 2 * creditnames[jj] +
4322 (j - sizeIrFx - cpiindexes[kk] - creditnames[jj]));
4323 r1[i][j] = r1[j][i] = c1[i0][j0] / std::sqrt(c1[i0][i0] * c1[j0][j0]);
4324 r2[i][j] = r2[j][i] = c2[i0][j0] / std::sqrt(c2[i0][i0] * c2[j0][j0]);
4325 if (std::fabs(r1[i][j] - c[i][j]) > tol) {
4326 BOOST_ERROR("failed to recover correlation matrix "
4327 "from "
4328 "Euler state process (i,j)=("
4329 << i << "," << j << "), (i0,j0)=(" << i0 << "," << j0
4330 << "), input correlation is " << c[i][j] << ", output is "
4331 << r1[i][j] << ", difference " << (c[i][j] - r1[i][j])
4332 << ", tolerance " << tol << " test configuration is "
4333 << currencies[ii] << " currencies and " << cpiindexes[kk]
4334 << " cpi indexes and " << creditnames[jj] << " credit names"
4335 << " and " << eqs[ee] << " equities");
4336 }
4337 if (k1 == k2) {
4338 if (std::fabs(r2[i][j] - c[i][j]) > tol) {
4339 BOOST_ERROR("failed to recover correlation "
4340 "matrix "
4341 "from "
4342 "exact state process (i,j)=("
4343 << i << "," << j << "), (i0,j0)=(" << i0 << "," << j0
4344 << "), input correlation is " << c[i][j] << ", output is "
4345 << r2[i][j] << ", difference " << (c[i][j] - r2[i][j])
4346 << ", tolerance " << tol << " test configuration is "
4347 << currencies[ii] << " currencies and " << cpiindexes[kk]
4348 << " cpi indexes and " << creditnames[jj] << " credit names"
4349 << " and " << eqs[ee] << " equities");
4350 }
4351 }
4352 }
4353 }
4354 }
4355 }
4356 } // for equity
4357 } // for creditnames
4358 } // for cpiindexes
4359 } // for currencies
4360} // testIrFxInfCrEqCorrelationRecovery

◆ BOOST_AUTO_TEST_CASE() [22/24]

BOOST_AUTO_TEST_CASE ( testCpiCalibrationByAlpha  )

Definition at line 4362 of file crossassetmodel.cpp.

4362 {
4363
4364 BOOST_TEST_MESSAGE("Testing calibration to ZC CPI Floors (using alpha) and repricing via MC...");
4365
4366 // set up IR-INF model, calibrate to given premiums and check
4367 // the result with a MC simulation
4368
4369 SavedSettings backup;
4370 Date refDate(30, July, 2015);
4371 Settings::instance().evaluationDate() = Date(30, July, 2015);
4372
4373 // IR
4374 Handle<YieldTermStructure> eurYts(QuantLib::ext::make_shared<FlatForward>(refDate, 0.01, Actual365Fixed()));
4375 QuantLib::ext::shared_ptr<Parametrization> ireur_p =
4376 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(EURCurrency(), eurYts, 0.01, 0.01);
4377
4378 // INF
4379 Real baseCPI = 100.0;
4380 std::vector<Date> infDates;
4381 std::vector<Real> infRates;
4382 infDates.push_back(Date(30, April, 2015));
4383 infDates.push_back(Date(30, July, 2015));
4384 infRates.push_back(0.0075);
4385 infRates.push_back(0.0075);
4386 Handle<ZeroInflationTermStructure> infEurTs(QuantLib::ext::make_shared<ZeroInflationCurve>(
4387 refDate, TARGET(), Actual365Fixed(), 3 * Months, Monthly, infDates, infRates));
4388 infEurTs->enableExtrapolation();
4389 Handle<ZeroInflationIndex> infIndex(QuantLib::ext::make_shared<EUHICPXT>(false, infEurTs));
4390
4391 infIndex->addFixing(Date(1, April, 2015), 100);
4392
4393 Real premium[] = { 0.0044, 0.0085, 0.0127, 0.0160, 0.0186 };
4394
4395 std::vector<QuantLib::ext::shared_ptr<BlackCalibrationHelper> > cpiHelpers;
4396 Array volStepTimes(4), noTimes(0);
4397 Array infVols(5, 0.01), infRev(1, 1.5); // !!
4398
4399 Time T;
4400 for (Size i = 1; i <= 5; ++i) {
4401 Date maturity = refDate + i * Years;
4402 QuantLib::ext::shared_ptr<CpiCapFloorHelper> h(new CpiCapFloorHelper(Option::Put, baseCPI, maturity, TARGET(),
4403 ModifiedFollowing, TARGET(), ModifiedFollowing,
4404 0.01, infIndex, 3 * Months, premium[i - 1]));
4405 Real t = inflationYearFraction(Monthly, false, Actual365Fixed(), infEurTs->baseDate(),
4406 h->instrument()->fixingDate());
4407 cpiHelpers.push_back(h);
4408 if (i <= 4)
4409 volStepTimes[i - 1] = t;
4410 T = t;
4411 }
4412
4413 QuantLib::ext::shared_ptr<InfDkPiecewiseConstantParametrization> infeur_p =
4414 QuantLib::ext::make_shared<InfDkPiecewiseConstantParametrization>(EURCurrency(), infEurTs, volStepTimes, infVols,
4415 noTimes, infRev);
4416
4417 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
4418 parametrizations.push_back(ireur_p);
4419 parametrizations.push_back(infeur_p);
4420
4421 QuantLib::ext::shared_ptr<CrossAssetModel> model =
4422 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, Matrix(), SalvagingAlgorithm::None);
4423
4424 model->setCorrelation(CrossAssetModel::AssetType::IR, 0, CrossAssetModel::AssetType::INF, 0, 0.33);
4425
4426 // pricing engine
4427 QuantLib::ext::shared_ptr<AnalyticDkCpiCapFloorEngine> engine =
4428 QuantLib::ext::make_shared<AnalyticDkCpiCapFloorEngine>(model, 0, baseCPI);
4429
4430 for (Size i = 0; i < cpiHelpers.size(); ++i) {
4431 cpiHelpers[i]->setPricingEngine(engine);
4432 }
4433
4434 // calibration
4435 LevenbergMarquardt lm;
4436 EndCriteria ec(1000, 500, 1E-8, 1E-8, 1E-8);
4437 model->calibrateInfDkVolatilitiesIterative(0, cpiHelpers, lm, ec);
4438
4439 for (Size i = 0; i < cpiHelpers.size(); ++i) {
4440 BOOST_TEST_MESSAGE("i=" << i << " modelvol=" << model->infdk(0)->parameterValues(0)[i] << " market="
4441 << cpiHelpers[i]->marketValue() << " model=" << cpiHelpers[i]->modelValue()
4442 << " diff=" << (cpiHelpers[i]->marketValue() - cpiHelpers[i]->modelValue()));
4443 }
4444
4445 // reprice last ZC floor with Monte Carlo
4446 Size n = 50000; // number of paths
4447 Size seed = 18; // rng seed
4448 Size steps = 1; // number of discretization steps
4449
4450 QuantLib::ext::shared_ptr<StochasticProcess> process = model->stateProcess();
4451 LowDiscrepancy::rsg_type sg = LowDiscrepancy::make_sequence_generator(process->factors() * steps, seed);
4452 TimeGrid grid(T, steps);
4453 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process)) {
4454 tmp->resetCache(grid.size() - 1);
4455 }
4456 MultiPathGenerator<LowDiscrepancy::rsg_type> pg(process, grid, sg, false);
4457
4458 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > floor;
4459
4460 Real K = std::pow(1.0 + 0.01, T);
4461
4462 for (Size i = 0; i < n; ++i) {
4463 Sample<MultiPath> path = pg.next();
4464 Size l = path.value[0].length() - 1;
4465 Real irz = path.value[0][l];
4466 Real infz = path.value[1][l];
4467 Real infy = path.value[2][l];
4468 Real I = model->infdkI(0, T, T, infz, infy).first;
4469 floor(std::max(-(I - K), 0.0) / model->numeraire(0, T, irz));
4470 }
4471
4472 BOOST_TEST_MESSAGE("mc floor 5y = " << mean(floor) << " +- ");
4473
4474 // check model calibration
4475 Real tol = 1.0E-12;
4476 for (Size i = 0; i < cpiHelpers.size(); ++i) {
4477 if (std::abs(cpiHelpers[i]->modelValue() - cpiHelpers[i]->marketValue()) > tol) {
4478 BOOST_ERROR("Model calibration for ZC CPI Floor #"
4479 << i << " failed, market premium is " << cpiHelpers[i]->marketValue() << ", model value is "
4480 << cpiHelpers[i]->modelValue() << ", difference is "
4481 << (cpiHelpers[i]->marketValue() - cpiHelpers[i]->modelValue()) << ", tolerance is " << tol);
4482 }
4483 }
4484 // check repricing with MC
4485 tol = 1.0E-5;
4486 Real mcPrice = mean(floor);
4487 if (std::abs(mcPrice - cpiHelpers[4]->modelValue()) > tol) {
4488 BOOST_ERROR("Failed to reprice 5y ZC CPI Floor with MC ("
4489 << mcPrice << "), analytical model price is " << cpiHelpers[4]->modelValue() << ", difference is "
4490 << mcPrice - cpiHelpers[4]->modelValue() << ", tolerance is " << tol);
4491 }
4492} // testCpiCalibrationByAlpha

◆ BOOST_AUTO_TEST_CASE() [23/24]

BOOST_AUTO_TEST_CASE ( testCpiCalibrationByH  )

Definition at line 4494 of file crossassetmodel.cpp.

4494 {
4495
4496 BOOST_TEST_MESSAGE("Testing calibration to ZC CPI Floors (using H) and repricing via MC...");
4497
4498 // set up IR-INF model, calibrate to given premiums and check
4499 // the result with a MC simulation
4500
4501 SavedSettings backup;
4502 Date refDate(30, July, 2015);
4503 Settings::instance().evaluationDate() = Date(30, July, 2015);
4504
4505 // IR
4506 Handle<YieldTermStructure> eurYts(QuantLib::ext::make_shared<FlatForward>(refDate, 0.01, Actual365Fixed()));
4507 QuantLib::ext::shared_ptr<Parametrization> ireur_p =
4508 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(EURCurrency(), eurYts, 0.01, 0.01);
4509
4510 // INF
4511 Real baseCPI = 100.0;
4512 std::vector<Date> infDates;
4513 std::vector<Real> infRates;
4514 infDates.push_back(Date(30, April, 2015));
4515 infDates.push_back(Date(30, July, 2015));
4516 infRates.push_back(0.0075);
4517 infRates.push_back(0.0075);
4518 Handle<ZeroInflationTermStructure> infEurTs(QuantLib::ext::make_shared<ZeroInflationCurve>(
4519 refDate, TARGET(), Actual365Fixed(), 3 * Months, Monthly, infDates, infRates));
4520 infEurTs->enableExtrapolation();
4521 Handle<ZeroInflationIndex> infIndex(QuantLib::ext::make_shared<EUHICPXT>(false, infEurTs));
4522 infIndex->addFixing(Date(1, April, 2015), 100);
4523 Size nMat = 14;
4524 Real premium[] = { 0.000555, 0.000813, 0.000928, 0.00127, 0.001616, 0.0019, 0.0023,
4525 0.0026, 0.0029, 0.0032, 0.0032, 0.0033, 0.0038, 0.0067 };
4526 Period maturity[] = { 1 * Years, 2 * Years, 3 * Years, 4 * Years, 5 * Years, 6 * Years, 7 * Years,
4527 8 * Years, 9 * Years, 10 * Years, 12 * Years, 15 * Years, 20 * Years, 30 * Years };
4528
4529 std::vector<QuantLib::ext::shared_ptr<BlackCalibrationHelper> > cpiHelpers;
4530 Array volStepTimes(13), noTimes(0);
4531 Array infVols(14, 0.0030), infRev(14, 1.0); // init vol and rev !!
4532 Real strike = 0.00; // strike !!
4533
4534 Time T = Null<Time>();
4535 for (Size i = 1; i <= nMat; ++i) {
4536 Date mat = refDate + maturity[i - 1];
4537 QuantLib::ext::shared_ptr<CpiCapFloorHelper> h(new CpiCapFloorHelper(Option::Put, baseCPI, mat, TARGET(),
4538 ModifiedFollowing, TARGET(), ModifiedFollowing,
4539 strike, infIndex, 3 * Months, premium[i - 1]));
4540 Real t = inflationYearFraction(Monthly, false, Actual365Fixed(), infEurTs->baseDate(),
4541 h->instrument()->fixingDate());
4542 cpiHelpers.push_back(h);
4543 if (i <= nMat - 1)
4544 volStepTimes[i - 1] = t;
4545 T = t;
4546 }
4547
4548 QuantLib::ext::shared_ptr<InfDkPiecewiseLinearParametrization> infeur_p =
4549 QuantLib::ext::make_shared<InfDkPiecewiseLinearParametrization>(EURCurrency(), infEurTs, volStepTimes, infVols,
4550 volStepTimes, infRev);
4551
4552 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
4553 parametrizations.push_back(ireur_p);
4554 parametrizations.push_back(infeur_p);
4555
4556 QuantLib::ext::shared_ptr<CrossAssetModel> model =
4557 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, Matrix(), SalvagingAlgorithm::None);
4558
4559 model->setCorrelation(CrossAssetModel::AssetType::IR, 0, CrossAssetModel::AssetType::INF, 0, 0.33);
4560
4561 // pricing engine
4562 QuantLib::ext::shared_ptr<AnalyticDkCpiCapFloorEngine> engine =
4563 QuantLib::ext::make_shared<AnalyticDkCpiCapFloorEngine>(model, 0, baseCPI);
4564
4565 for (Size i = 0; i < cpiHelpers.size(); ++i) {
4566 cpiHelpers[i]->setPricingEngine(engine);
4567 }
4568
4569 // calibration
4570 LevenbergMarquardt lm;
4571 EndCriteria ec(1000, 500, 1E-8, 1E-8, 1E-8);
4572 // model->calibrateInfDkVolatilitiesIterative(0, cpiHelpers, lm, ec);
4573 model->calibrateInfDkReversionsIterative(0, cpiHelpers, lm, ec);
4574
4575 for (Size i = 0; i < cpiHelpers.size(); ++i) {
4576 BOOST_TEST_MESSAGE("i=" << i << " modelvol=" << model->infdk(0)->parameterValues(0)[i] << " modelrev="
4577 << model->infdk(0)->parameterValues(1)[i] << " market=" << cpiHelpers[i]->marketValue()
4578 << " model=" << cpiHelpers[i]->modelValue()
4579 << " diff=" << (cpiHelpers[i]->marketValue() - cpiHelpers[i]->modelValue()));
4580 }
4581
4582 // reprice last ZC floor with Monte Carlo
4583 Size n = 100000; // number of paths
4584 Size seed = 18; // rng seed
4585 Size steps = 1; // number of discretization steps
4586
4587 QuantLib::ext::shared_ptr<StochasticProcess> process = model->stateProcess();
4588 LowDiscrepancy::rsg_type sg = LowDiscrepancy::make_sequence_generator(process->factors() * steps, seed);
4589 TimeGrid grid(T, steps);
4590 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process)) {
4591 tmp->resetCache(grid.size() - 1);
4592 }
4593 MultiPathGenerator<LowDiscrepancy::rsg_type> pg(process, grid, sg, false);
4594
4595 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > floor;
4596
4597 Real K = std::pow(1.0 + strike, T);
4598
4599 for (Size i = 0; i < n; ++i) {
4600 Sample<MultiPath> path = pg.next();
4601 Size l = path.value[0].length() - 1;
4602 Real irz = path.value[0][l];
4603 Real infz = path.value[1][l];
4604 Real infy = path.value[2][l];
4605 Real I = model->infdkI(0, T, T, infz, infy).first;
4606 floor(std::max(-(I - K), 0.0) / model->numeraire(0, T, irz));
4607 }
4608
4609 BOOST_TEST_MESSAGE("mc (low disc) floor last = " << mean(floor) << " +- " << error_of<tag::mean>(floor));
4610
4611 // check model calibration
4612 Real tol = 1.0E-12;
4613 for (Size i = 0; i < cpiHelpers.size(); ++i) {
4614 if (std::abs(cpiHelpers[i]->modelValue() - cpiHelpers[i]->marketValue()) > tol) {
4615 BOOST_ERROR("Model calibration for ZC CPI Floor #"
4616 << i << " failed, market premium is " << cpiHelpers[i]->marketValue() << ", model value is "
4617 << cpiHelpers[i]->modelValue() << ", difference is "
4618 << (cpiHelpers[i]->marketValue() - cpiHelpers[i]->modelValue()) << ", tolerance is " << tol);
4619 }
4620 }
4621 // check repricing with MC
4622 tol = 2.0E-4;
4623 Real mcPrice = mean(floor);
4624 if (std::abs(mcPrice - cpiHelpers[nMat - 1]->modelValue()) > tol) {
4625 BOOST_ERROR("Failed to reprice last ZC CPI Floor with MC ("
4626 << mcPrice << "), analytical model price is " << cpiHelpers[nMat - 1]->modelValue()
4627 << ", difference is " << mcPrice - cpiHelpers[nMat - 1]->modelValue() << ", tolerance is " << tol);
4628 }
4629} // testCpiCalibrationByH

◆ BOOST_AUTO_TEST_CASE() [24/24]

BOOST_AUTO_TEST_CASE ( testCrCalibration  )

Definition at line 4631 of file crossassetmodel.cpp.

4631 {
4632
4633 BOOST_TEST_MESSAGE("Testing calibration to CDS Options and repricing via MC...");
4634
4635 // set up IR-CR model, calibrate to given CDS Options and check
4636 // the result with a MC simulation
4637
4638 SavedSettings backup;
4639 Date refDate(30, July, 2015);
4640 Settings::instance().evaluationDate() = Date(30, July, 2015);
4641
4642 // IR (zero vol)
4643 Handle<YieldTermStructure> eurYts(QuantLib::ext::make_shared<FlatForward>(refDate, 0.01, Actual365Fixed()));
4644 QuantLib::ext::shared_ptr<Parametrization> ireur_p =
4645 QuantLib::ext::make_shared<IrLgm1fConstantParametrization>(EURCurrency(), eurYts, 0.00, 0.01);
4646
4647 // CR
4648 Handle<DefaultProbabilityTermStructure> prob(QuantLib::ext::make_shared<FlatHazardRate>(refDate, 0.01, Actual365Fixed()));
4649
4650 Size nMat = 10;
4651 Real impliedVols[] = { 0.10, 0.12, 0.14, 0.16, 0.18, 0.2, 0.21, 0.215, 0.22, 0.225 };
4652 Period maturity[] = { 1 * Years, 2 * Years, 3 * Years, 4 * Years, 5 * Years,
4653 6 * Years, 7 * Years, 8 * Years, 9 * Years, 10 * Years };
4654
4655 std::vector<QuantLib::ext::shared_ptr<BlackCalibrationHelper> > cdsoHelpers;
4656 Array volStepTimes(nMat - 1), noTimes(0);
4657 Array crVols(nMat, 0.0030), crRev(nMat, 0.01); // init vol and rev
4658
4659 Time T = Null<Time>();
4660 Date lastMat;
4661 for (Size i = 1; i <= nMat; ++i) {
4662 Date mat = refDate + maturity[i - 1];
4663 Schedule schedule(mat, mat + 10 * Years, 3 * Months, TARGET(), Following, Following, DateGeneration::CDS,
4664 false);
4665 // ensure that cds starts after option expiry
4666 if (schedule.date(0) < mat)
4667 schedule = Schedule(schedule.date(1), mat + 10 * Years, 3 * Months, TARGET(), Following, Following,
4668 DateGeneration::CDS, false);
4669 QL_REQUIRE(schedule.date(0) >= mat,
4670 "CDS start (" << schedule.date(0) << ") should be on or after option maturity (" << mat << ")");
4671 QuantLib::ext::shared_ptr<CdsOptionHelper> h(
4672 new CdsOptionHelper(mat, Handle<Quote>(QuantLib::ext::make_shared<SimpleQuote>(impliedVols[i - 1])),
4673 Protection::Buyer, schedule, Following, Actual360(), prob, 0.4, eurYts));
4674 Real t = eurYts->timeFromReference(mat);
4675 cdsoHelpers.push_back(h);
4676 if (i <= nMat - 1)
4677 volStepTimes[i - 1] = t;
4678 T = t;
4679 lastMat = mat;
4680 }
4681
4682 QuantLib::ext::shared_ptr<CrLgm1fPiecewiseConstantParametrization> creur_p =
4683 QuantLib::ext::make_shared<CrLgm1fPiecewiseConstantParametrization>(EURCurrency(), prob, volStepTimes, crVols,
4684 volStepTimes, crRev);
4685
4686 std::vector<QuantLib::ext::shared_ptr<Parametrization> > parametrizations;
4687 parametrizations.push_back(ireur_p);
4688 parametrizations.push_back(creur_p);
4689
4690 QuantLib::ext::shared_ptr<CrossAssetModel> model =
4691 QuantLib::ext::make_shared<CrossAssetModel>(parametrizations, Matrix(), SalvagingAlgorithm::None);
4692
4693 model->setCorrelation(CrossAssetModel::AssetType::IR, 0, CrossAssetModel::AssetType::CR, 0, 0.33);
4694
4695 // pricing engine
4696 QuantLib::ext::shared_ptr<AnalyticLgmCdsOptionEngine> engine =
4697 QuantLib::ext::make_shared<AnalyticLgmCdsOptionEngine>(model, 0, 0, 0.4);
4698
4699 for (Size i = 0; i < cdsoHelpers.size(); ++i) {
4700 cdsoHelpers[i]->setPricingEngine(engine);
4701 }
4702
4703 // calibration
4704 LevenbergMarquardt lm;
4705 EndCriteria ec(1000, 500, 1E-8, 1E-8, 1E-8);
4706 model->calibrateCrLgm1fVolatilitiesIterative(0, cdsoHelpers, lm, ec);
4707
4708 for (Size i = 0; i < cdsoHelpers.size(); ++i) {
4709 BOOST_TEST_MESSAGE("i=" << i << " modelvol=" << model->crlgm1f(0)->parameterValues(0)[i]
4710 << " modelrev=" << model->crlgm1f(0)->parameterValues(1)[i] << " market="
4711 << cdsoHelpers[i]->marketValue() << " model=" << cdsoHelpers[i]->modelValue()
4712 << " diff=" << (cdsoHelpers[i]->marketValue() - cdsoHelpers[i]->modelValue()));
4713 }
4714
4715 // check model calibration
4716 Real tol = 1.0E-12;
4717 for (Size i = 0; i < cdsoHelpers.size(); ++i) {
4718 if (std::abs(cdsoHelpers[i]->modelValue() - cdsoHelpers[i]->marketValue()) > tol) {
4719 BOOST_ERROR("Model calibration for CDSO #"
4720 << i << " failed, market premium is " << cdsoHelpers[i]->marketValue() << ", model value is "
4721 << cdsoHelpers[i]->modelValue() << ", difference is "
4722 << (cdsoHelpers[i]->marketValue() - cdsoHelpers[i]->modelValue()) << ", tolerance is " << tol);
4723 }
4724 }
4725
4726 Real lastModelValue = cdsoHelpers[nMat - 1]->modelValue();
4727
4728 // reprice last CDSO with Monte Carlo
4729 // note that the IR vol is zero (same as assumption in CDSO analytic engine)
4730 Size n = 10000; // number of paths
4731 Size seed = 18; // rng seed
4732 Size steps = 1; // number of discretization steps
4733
4734 QuantLib::ext::shared_ptr<StochasticProcess> process = model->stateProcess();
4735 LowDiscrepancy::rsg_type sg = LowDiscrepancy::make_sequence_generator(process->factors() * steps, seed);
4736 TimeGrid grid(T, steps);
4737 if (auto tmp = QuantLib::ext::dynamic_pointer_cast<CrossAssetStateProcess>(process)) {
4738 tmp->resetCache(grid.size() - 1);
4739 }
4740 MultiPathGenerator<LowDiscrepancy::rsg_type> pg(process, grid, sg, false);
4741
4742 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean> > > cdso;
4743
4744 QuantLib::ext::shared_ptr<CreditDefaultSwap> underlying =
4745 QuantLib::ext::static_pointer_cast<CdsOptionHelper>(cdsoHelpers.back())->underlying();
4746 Real K = underlying->fairSpreadClean();
4747 BOOST_TEST_MESSAGE("Last CDSO fair spread is " << K);
4748
4749 Settings::instance().evaluationDate() = lastMat;
4750 QuantLib::ext::shared_ptr<LgmImpliedDefaultTermStructure> probMc =
4751 QuantLib::ext::make_shared<LgmImpliedDefaultTermStructure>(model, 0, 0);
4752 QuantLib::ext::shared_ptr<LgmImpliedYieldTermStructure> ytsMc =
4753 QuantLib::ext::make_shared<LgmImpliedYieldTermStructure>(model->lgm(0));
4754 QuantLib::ext::shared_ptr<QuantExt::MidPointCdsEngine> dynamicEngine = QuantLib::ext::make_shared<QuantExt::MidPointCdsEngine>(
4755 Handle<DefaultProbabilityTermStructure>(probMc), 0.4, Handle<YieldTermStructure>(ytsMc));
4756 underlying->setPricingEngine(dynamicEngine);
4757
4758 for (Size i = 0; i < n; ++i) {
4759 Sample<MultiPath> path = pg.next();
4760 Size l = path.value[0].length() - 1;
4761 Real irz = path.value[0][l];
4762 Real crz = path.value[1][l];
4763 Real cry = path.value[2][l];
4764 probMc->move(lastMat, crz, cry);
4765 ytsMc->move(lastMat, irz);
4766 Real surv = model->crlgm1fS(0, 0, T, T, crz, cry).first;
4767 Real npv = surv * std::max(underlying->NPV(), 0.0) / model->numeraire(0, T, irz);
4768 cdso(npv);
4769 }
4770
4771 BOOST_TEST_MESSAGE("mc (low disc) cdso last = " << mean(cdso) << " +- " << error_of<tag::mean>(cdso));
4772
4773 // check repricing with MC
4774 tol = 3.0E-4;
4775 Real mcPrice = mean(cdso);
4776 if (std::abs(mcPrice - lastModelValue) > tol) {
4777 BOOST_ERROR("Failed to reprice last CDSO with MC (" << mcPrice << "), analytical model price is "
4778 << lastModelValue << ", difference is "
4779 << mcPrice - lastModelValue << ", tolerance is " << tol);
4780 }
4781} // testCrCalibration

Variable Documentation

◆ infEurFlags

vector<bool> infEurFlags { true, false }

Definition at line 2281 of file crossassetmodel.cpp.

◆ infGbpFlags

vector<bool> infGbpFlags { true, false }

Definition at line 2282 of file crossassetmodel.cpp.

◆ flatVolsFlags

vector<bool> flatVolsFlags { true, false }

Definition at line 2283 of file crossassetmodel.cpp.

◆ driftFreeState

vector<bool> driftFreeState { true, false }

Definition at line 2284 of file crossassetmodel.cpp.

Generated by Doxygen 1.9.5