| 1 | /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ |
| 2 | |
| 3 | /* |
| 4 | Copyright (C) 2007 Ferdinando Ametrano |
| 5 | Copyright (C) 2007 Marco Bianchetti |
| 6 | Copyright (C) 2007 Cristina Duminuco |
| 7 | Copyright (C) 2007 Mark Joshi |
| 8 | |
| 9 | This file is part of QuantLib, a free-software/open-source library |
| 10 | for financial quantitative analysts and developers - http://quantlib.org/ |
| 11 | |
| 12 | QuantLib is free software: you can redistribute it and/or modify it |
| 13 | under the terms of the QuantLib license. You should have received a |
| 14 | copy of the license along with this program; if not, please email |
| 15 | <quantlib-dev@lists.sf.net>. The license is also available online at |
| 16 | <http://quantlib.org/license.shtml>. |
| 17 | |
| 18 | This program is distributed in the hope that it will be useful, but WITHOUT |
| 19 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS |
| 20 | FOR A PARTICULAR PURPOSE. See the license for more details. |
| 21 | */ |
| 22 | |
| 23 | #include "marketmodel_smmcaplethomocalibration.hpp" |
| 24 | #include "utilities.hpp" |
| 25 | |
| 26 | #include <ql/math/optimization/spherecylinder.hpp> |
| 27 | #include <ql/models/marketmodels/models/fwdtocotswapadapter.hpp> |
| 28 | #include <ql/models/marketmodels/models/fwdperiodadapter.hpp> |
| 29 | #include <ql/models/marketmodels/models/capletcoterminalmaxhomogeneity.hpp> |
| 30 | #include <ql/models/marketmodels/models/alphaformconcrete.hpp> |
| 31 | #include <ql/models/marketmodels/correlations/cotswapfromfwdcorrelation.hpp> |
| 32 | #include <ql/models/marketmodels/correlations/timehomogeneousforwardcorrelation.hpp> |
| 33 | #include <ql/models/marketmodels/models/piecewiseconstantabcdvariance.hpp> |
| 34 | #include <ql/models/marketmodels/models/capletcoterminalswaptioncalibration.hpp> |
| 35 | #include <ql/models/marketmodels/models/cotswaptofwdadapter.hpp> |
| 36 | #include <ql/models/marketmodels/models/pseudorootfacade.hpp> |
| 37 | #include <ql/models/marketmodels/products/multistep/multistepcoterminalswaps.hpp> |
| 38 | #include <ql/models/marketmodels/products/multistep/multistepcoterminalswaptions.hpp> |
| 39 | #include <ql/models/marketmodels/products/multistep/multistepswap.hpp> |
| 40 | #include <ql/models/marketmodels/products/multiproductcomposite.hpp> |
| 41 | #include <ql/models/marketmodels/accountingengine.hpp> |
| 42 | #include <ql/models/marketmodels/utilities.hpp> |
| 43 | #include <ql/models/marketmodels/evolvers/lognormalcotswapratepc.hpp> |
| 44 | #include <ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp> |
| 45 | #include <ql/models/marketmodels/correlations/expcorrelations.hpp> |
| 46 | #include <ql/models/marketmodels/models/flatvol.hpp> |
| 47 | #include <ql/models/marketmodels/models/abcdvol.hpp> |
| 48 | #include <ql/models/marketmodels/browniangenerators/mtbrowniangenerator.hpp> |
| 49 | #include <ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp> |
| 50 | #include <ql/models/marketmodels/swapforwardmappings.hpp> |
| 51 | #include <ql/models/marketmodels/curvestates/coterminalswapcurvestate.hpp> |
| 52 | #include <ql/methods/montecarlo/genericlsregression.hpp> |
| 53 | #include <ql/legacy/libormarketmodels/lmlinexpcorrmodel.hpp> |
| 54 | #include <ql/legacy/libormarketmodels/lmextlinexpvolmodel.hpp> |
| 55 | #include <ql/time/schedule.hpp> |
| 56 | #include <ql/time/calendars/nullcalendar.hpp> |
| 57 | #include <ql/time/daycounters/simpledaycounter.hpp> |
| 58 | #include <ql/pricingengines/blackformula.hpp> |
| 59 | #include <ql/pricingengines/blackcalculator.hpp> |
| 60 | #include <ql/utilities/dataformatters.hpp> |
| 61 | #include <ql/math/integrals/segmentintegral.hpp> |
| 62 | #include <ql/math/statistics/convergencestatistics.hpp> |
| 63 | #include <ql/math/optimization/simplex.hpp> |
| 64 | #include <ql/math/statistics/sequencestatistics.hpp> |
| 65 | #include <ql/models/marketmodels/models/capletcoterminalperiodic.hpp> |
| 66 | #include <ql/models/marketmodels/models/volatilityinterpolationspecifierabcd.hpp> |
| 67 | #include <sstream> |
| 68 | |
| 69 | using namespace QuantLib; |
| 70 | using namespace boost::unit_test_framework; |
| 71 | |
| 72 | using std::fabs; |
| 73 | using std::sqrt; |
| 74 | |
| 75 | namespace market_model_smm_caplet_homo_calibration_test { |
| 76 | |
| 77 | Date todaysDate_, startDate_, endDate_; |
| 78 | std::vector<Time> rateTimes_; |
| 79 | std::vector<Real> accruals_; |
| 80 | Calendar calendar_; |
| 81 | DayCounter dayCounter_; |
| 82 | std::vector<Rate> todaysForwards_, todaysSwaps_; |
| 83 | std::vector<Real> coterminalAnnuity_; |
| 84 | Size numberOfFactors_; |
| 85 | Real alpha_, alphaMax_, alphaMin_; |
| 86 | Spread displacement_; |
| 87 | std::vector<DiscountFactor> todaysDiscounts_; |
| 88 | std::vector<Volatility> swaptionDisplacedVols_, swaptionVols_; |
| 89 | std::vector<Volatility> capletDisplacedVols_, capletVols_; |
| 90 | Real a_, b_, c_, d_; |
| 91 | Real longTermCorrelation_, beta_; |
| 92 | Size measureOffset_; |
| 93 | unsigned long seed_; |
| 94 | Size paths_, trainingPaths_; |
| 95 | bool printReport_ = false; |
| 96 | |
| 97 | void setup() { |
| 98 | |
| 99 | // Times |
| 100 | calendar_ = NullCalendar(); |
| 101 | todaysDate_ = Settings::instance().evaluationDate(); |
| 102 | //startDate = todaysDate + 5*Years; |
| 103 | endDate_ = todaysDate_ + 66*Months; |
| 104 | Schedule dates(todaysDate_, endDate_, Period(Semiannual), |
| 105 | calendar_, Following, Following, DateGeneration::Backward, false); |
| 106 | rateTimes_ = std::vector<Time>(dates.size()-1); |
| 107 | accruals_ = std::vector<Real>(rateTimes_.size()-1); |
| 108 | dayCounter_ = SimpleDayCounter(); |
| 109 | for (Size i=1; i<dates.size(); ++i) |
| 110 | rateTimes_[i-1] = dayCounter_.yearFraction(d1: todaysDate_, d2: dates[i]); |
| 111 | for (Size i=1; i<rateTimes_.size(); ++i) |
| 112 | accruals_[i-1] = rateTimes_[i] - rateTimes_[i-1]; |
| 113 | |
| 114 | // Rates & displacement |
| 115 | todaysForwards_ = std::vector<Rate>(accruals_.size()); |
| 116 | numberOfFactors_ = 3; |
| 117 | alpha_ = 0.0; |
| 118 | alphaMax_ = 1.0; |
| 119 | alphaMin_ = -1.0; |
| 120 | displacement_ = 0.0; |
| 121 | for (Size i=0; i<todaysForwards_.size(); ++i) { |
| 122 | todaysForwards_[i] = 0.03 + 0.0025*i; |
| 123 | // todaysForwards_[i] = 0.03; |
| 124 | } |
| 125 | LMMCurveState curveState_lmm(rateTimes_); |
| 126 | curveState_lmm.setOnForwardRates(fwdRates: todaysForwards_); |
| 127 | todaysSwaps_ = curveState_lmm.coterminalSwapRates(); |
| 128 | |
| 129 | // Discounts |
| 130 | todaysDiscounts_ = std::vector<DiscountFactor>(rateTimes_.size()); |
| 131 | todaysDiscounts_[0] = 0.95; |
| 132 | for (Size i=1; i<rateTimes_.size(); ++i) |
| 133 | todaysDiscounts_[i] = todaysDiscounts_[i-1] / |
| 134 | (1.0+todaysForwards_[i-1]*accruals_[i-1]); |
| 135 | |
| 136 | //// Swaption Volatilities |
| 137 | //Volatility mktSwaptionVols[] = { |
| 138 | // 0.15541283, |
| 139 | // 0.18719678, |
| 140 | // 0.20890740, |
| 141 | // 0.22318179, |
| 142 | // 0.23212717, |
| 143 | // 0.23731450, |
| 144 | // 0.23988649, |
| 145 | // 0.24066384, |
| 146 | // 0.24023111, |
| 147 | // 0.23900189, |
| 148 | // 0.23726699, |
| 149 | // 0.23522952, |
| 150 | // 0.23303022, |
| 151 | // 0.23076564, |
| 152 | // 0.22850101, |
| 153 | // 0.22627951, |
| 154 | // 0.22412881, |
| 155 | // 0.22206569, |
| 156 | // 0.22009939 |
| 157 | //}; |
| 158 | |
| 159 | //a = -0.0597; |
| 160 | //b = 0.1677; |
| 161 | //c = 0.5403; |
| 162 | //d = 0.1710; |
| 163 | |
| 164 | a_ = 0.0; |
| 165 | b_ = 0.17; |
| 166 | c_ = 1.0; |
| 167 | d_ = 0.10; |
| 168 | |
| 169 | Volatility mktCapletVols[] = { |
| 170 | 0.1640, |
| 171 | 0.1740, |
| 172 | 0.1840, |
| 173 | 0.1940, |
| 174 | 0.1840, |
| 175 | 0.1740, |
| 176 | 0.1640, |
| 177 | 0.1540, |
| 178 | 0.1440, |
| 179 | 0.1340376439125532 |
| 180 | }; |
| 181 | |
| 182 | //swaptionDisplacedVols = std::vector<Volatility>(todaysSwaps.size()); |
| 183 | //swaptionVols = std::vector<Volatility>(todaysSwaps.size()); |
| 184 | //capletDisplacedVols = std::vector<Volatility>(todaysSwaps.size()); |
| 185 | capletVols_.resize(new_size: todaysSwaps_.size()); |
| 186 | for (Size i=0; i<todaysSwaps_.size(); i++) { |
| 187 | // swaptionDisplacedVols[i] = todaysSwaps[i]*mktSwaptionVols[i]/ |
| 188 | // (todaysSwaps[i]+displacement); |
| 189 | // swaptionVols[i]= mktSwaptionVols[i]; |
| 190 | // capletDisplacedVols[i] = todaysForwards[i]*mktCapletVols[i]/ |
| 191 | // (todaysForwards[i]+displacement); |
| 192 | capletVols_[i]= mktCapletVols[i]; |
| 193 | } |
| 194 | |
| 195 | // Cap/Floor Correlation |
| 196 | longTermCorrelation_ = 0.5; |
| 197 | beta_ = 0.2; |
| 198 | measureOffset_ = 5; |
| 199 | |
| 200 | // Monte Carlo |
| 201 | seed_ = 42; |
| 202 | |
| 203 | #ifdef _DEBUG |
| 204 | paths_ = 127; |
| 205 | trainingPaths_ = 31; |
| 206 | #else |
| 207 | paths_ = 32767; //262144-1; //; // 2^15-1 |
| 208 | trainingPaths_ = 8191; // 2^13-1 |
| 209 | #endif |
| 210 | } |
| 211 | |
| 212 | enum MarketModelType { ExponentialCorrelationFlatVolatility, |
| 213 | ExponentialCorrelationAbcdVolatility/*, |
| 214 | CalibratedMM*/ |
| 215 | }; |
| 216 | |
| 217 | enum MeasureType { ProductSuggested, Terminal, |
| 218 | MoneyMarket, MoneyMarketPlus }; |
| 219 | |
| 220 | enum EvolverType { Ipc, Pc , NormalPc }; |
| 221 | |
| 222 | } |
| 223 | |
| 224 | |
| 225 | |
| 226 | void MarketModelSmmCapletHomoCalibrationTest::testFunction() { |
| 227 | |
| 228 | BOOST_TEST_MESSAGE("Testing max homogeneity caplet calibration " |
| 229 | "in a lognormal coterminal swap market model..." ); |
| 230 | |
| 231 | using namespace market_model_smm_caplet_homo_calibration_test; |
| 232 | |
| 233 | setup(); |
| 234 | |
| 235 | Size numberOfRates = todaysForwards_.size(); |
| 236 | |
| 237 | EvolutionDescription evolution(rateTimes_); |
| 238 | // Size numberOfSteps = evolution.numberOfSteps(); |
| 239 | |
| 240 | ext::shared_ptr<PiecewiseConstantCorrelation> fwdCorr(new |
| 241 | ExponentialForwardCorrelation(rateTimes_, |
| 242 | longTermCorrelation_, |
| 243 | beta_)); |
| 244 | |
| 245 | ext::shared_ptr<LMMCurveState> cs(new LMMCurveState(rateTimes_)); |
| 246 | cs->setOnForwardRates(fwdRates: todaysForwards_); |
| 247 | |
| 248 | ext::shared_ptr<PiecewiseConstantCorrelation> corr(new |
| 249 | CotSwapFromFwdCorrelation(fwdCorr, *cs, displacement_)); |
| 250 | |
| 251 | std::vector<ext::shared_ptr<PiecewiseConstantVariance> > |
| 252 | swapVariances(numberOfRates); |
| 253 | for (Size i=0; i<numberOfRates; ++i) { |
| 254 | swapVariances[i] = ext::shared_ptr<PiecewiseConstantVariance>(new |
| 255 | PiecewiseConstantAbcdVariance(a_, b_, c_, d_, |
| 256 | i, rateTimes_)); |
| 257 | } |
| 258 | |
| 259 | // create calibrator |
| 260 | Real caplet0Swaption1Priority = 1.0; |
| 261 | if (printReport_) { |
| 262 | BOOST_TEST_MESSAGE("caplet market vols: " << std::fixed << |
| 263 | std::setprecision(4) << io::sequence(capletVols_)); |
| 264 | BOOST_TEST_MESSAGE("caplet0Swapt1Prior: " << caplet0Swaption1Priority); |
| 265 | } |
| 266 | CTSMMCapletMaxHomogeneityCalibration calibrator(evolution, |
| 267 | corr, |
| 268 | swapVariances, |
| 269 | capletVols_, |
| 270 | cs, |
| 271 | displacement_, |
| 272 | caplet0Swaption1Priority); |
| 273 | // calibrate |
| 274 | Natural maxIterations = 10; |
| 275 | Real capletTolerance = 1e-4; // i.e. 1 bp |
| 276 | Natural innerMaxIterations = 100; |
| 277 | Real innerTolerance = 1e-8; |
| 278 | if (printReport_) { |
| 279 | BOOST_TEST_MESSAGE("numberOfFactors: " << numberOfFactors_); |
| 280 | BOOST_TEST_MESSAGE("maxIterations: " << maxIterations); |
| 281 | BOOST_TEST_MESSAGE("capletTolerance: " << io::rate(capletTolerance)); |
| 282 | BOOST_TEST_MESSAGE("innerMaxIterations: " << innerMaxIterations); |
| 283 | BOOST_TEST_MESSAGE("innerTolerance: " << io::rate(innerTolerance)); |
| 284 | } |
| 285 | bool result = calibrator.calibrate(numberOfFactors: numberOfFactors_, |
| 286 | maxIterations, |
| 287 | tolerance: capletTolerance, |
| 288 | innerMaxIterations, |
| 289 | innerTolerance); |
| 290 | if (!result) |
| 291 | BOOST_ERROR("calibration failed" ); |
| 292 | |
| 293 | const std::vector<Matrix>& swapPseudoRoots = calibrator.swapPseudoRoots(); |
| 294 | ext::shared_ptr<MarketModel> smm(new |
| 295 | PseudoRootFacade(swapPseudoRoots, |
| 296 | rateTimes_, |
| 297 | cs->coterminalSwapRates(), |
| 298 | std::vector<Spread>(numberOfRates, displacement_))); |
| 299 | ext::shared_ptr<MarketModel> flmm(new CotSwapToFwdAdapter(smm)); |
| 300 | Matrix capletTotCovariance = flmm->totalCovariance(endIndex: numberOfRates-1); |
| 301 | |
| 302 | std::vector<Volatility> capletVols(numberOfRates); |
| 303 | for (Size i=0; i<numberOfRates; ++i) { |
| 304 | capletVols[i] = std::sqrt(x: capletTotCovariance[i][i]/rateTimes_[i]); |
| 305 | } |
| 306 | if (printReport_) { |
| 307 | BOOST_TEST_MESSAGE("caplet smm implied vols: " << std::fixed << |
| 308 | std::setprecision(4) << io::sequence(capletVols)); |
| 309 | BOOST_TEST_MESSAGE("failures: " << calibrator.failures()); |
| 310 | BOOST_TEST_MESSAGE("deformationSize: " << calibrator.deformationSize()); |
| 311 | BOOST_TEST_MESSAGE("capletRmsError: " << calibrator.capletRmsError()); |
| 312 | BOOST_TEST_MESSAGE("capletMaxError: " << calibrator.capletMaxError()); |
| 313 | BOOST_TEST_MESSAGE("swaptionRmsError: " << calibrator.swaptionRmsError()); |
| 314 | BOOST_TEST_MESSAGE("swaptionMaxError: " << calibrator.swaptionMaxError()); |
| 315 | } |
| 316 | |
| 317 | // check perfect swaption fit |
| 318 | Real error, swapTolerance = 1e-14; |
| 319 | Matrix swapTerminalCovariance(numberOfRates, numberOfRates, 0.0); |
| 320 | for (Size i=0; i<numberOfRates; ++i) { |
| 321 | Volatility expSwaptionVol = swapVariances[i]->totalVolatility(i); |
| 322 | swapTerminalCovariance += swapPseudoRoots[i] * transpose(m: swapPseudoRoots[i]); |
| 323 | Volatility swaptionVol = std::sqrt(x: swapTerminalCovariance[i][i]/rateTimes_[i]); |
| 324 | error = std::fabs(x: swaptionVol-expSwaptionVol); |
| 325 | if (error>swapTolerance) |
| 326 | BOOST_ERROR("failed to reproduce " << io::ordinal(i+1) << " swaption vol:" |
| 327 | "\n expected: " << io::rate(expSwaptionVol) << |
| 328 | "\n realized: " << io::rate(swaptionVol) << |
| 329 | "\n error: " << error << |
| 330 | "\n tolerance: " << swapTolerance); |
| 331 | } |
| 332 | |
| 333 | // check caplet fit |
| 334 | for (Size i=0; i<numberOfRates; ++i) { |
| 335 | error = std::fabs(x: capletVols[i]-capletVols_[i]); |
| 336 | if (error>capletTolerance) |
| 337 | BOOST_ERROR("failed to reproduce " << io::ordinal(i+1) << " caplet vol:" |
| 338 | "\n expected: " << io::rate(capletVols_[i]) << |
| 339 | "\n realized: " << io::rate(capletVols[i]) << |
| 340 | "\n percentage error: " << error/capletVols_[i] << |
| 341 | "\n error: " << error << |
| 342 | "\n tolerance: " << capletTolerance); |
| 343 | } |
| 344 | |
| 345 | Size period =2; |
| 346 | Size offset =0; |
| 347 | std::vector<Spread> adaptedDisplacements; |
| 348 | ext::shared_ptr<MarketModel> adapted(new FwdPeriodAdapter(flmm,period,offset,adaptedDisplacements)); |
| 349 | // FwdToCotSwapAdapter newSwapMM(adapted); |
| 350 | // for (Size i=0; i < newSwapMM.numberOfRates(); ++i) |
| 351 | // BOOST_TEST_MESSAGE("swap MM time dependent vols: "<< i << std::fixed << |
| 352 | // std::setprecision(6) << Array(newSwapMM.timeDependentVolatility(i))); |
| 353 | |
| 354 | |
| 355 | |
| 356 | |
| 357 | |
| 358 | |
| 359 | } |
| 360 | |
| 361 | |
| 362 | |
| 363 | void MarketModelSmmCapletHomoCalibrationTest::testPeriodFunction() |
| 364 | { |
| 365 | |
| 366 | BOOST_TEST_MESSAGE("Testing max homogeneity periodic caplet calibration " |
| 367 | "in a lognormal coterminal swap market model..." ); |
| 368 | |
| 369 | using namespace market_model_smm_caplet_homo_calibration_test; |
| 370 | |
| 371 | setup(); |
| 372 | |
| 373 | Size numberOfRates = todaysForwards_.size(); |
| 374 | Size period=2; |
| 375 | Size offset = numberOfRates % period; |
| 376 | Size numberBigRates = numberOfRates / period; |
| 377 | |
| 378 | EvolutionDescription evolution(rateTimes_); |
| 379 | |
| 380 | std::vector<Time> bigRateTimes(numberBigRates+1); |
| 381 | |
| 382 | for (Size i=0; i <= numberBigRates; ++i) |
| 383 | bigRateTimes[i] = rateTimes_[i*period+offset]; |
| 384 | |
| 385 | ext::shared_ptr<PiecewiseConstantCorrelation> fwdCorr(new |
| 386 | ExponentialForwardCorrelation(rateTimes_, |
| 387 | longTermCorrelation_, |
| 388 | beta_)); |
| 389 | |
| 390 | ext::shared_ptr<LMMCurveState> cs(new LMMCurveState(rateTimes_)); |
| 391 | cs->setOnForwardRates(fwdRates: todaysForwards_); |
| 392 | |
| 393 | ext::shared_ptr<PiecewiseConstantCorrelation> corr(new |
| 394 | CotSwapFromFwdCorrelation(fwdCorr, *cs, displacement_)); |
| 395 | |
| 396 | std::vector<PiecewiseConstantAbcdVariance > |
| 397 | swapVariances; |
| 398 | for (Size i=0; i<numberBigRates; ++i) { |
| 399 | swapVariances.emplace_back(args&: a_, args&: b_, args&: c_, args&: d_, args&: i, args&: bigRateTimes); |
| 400 | } |
| 401 | |
| 402 | VolatilityInterpolationSpecifierabcd varianceInterpolator(period, offset, swapVariances, // these should be associated with the long rates |
| 403 | rateTimes_ // these should be associated with the shorter rates |
| 404 | ); |
| 405 | |
| 406 | |
| 407 | // create calibrator |
| 408 | Real caplet0Swaption1Priority = 1.0; |
| 409 | if (printReport_) { |
| 410 | BOOST_TEST_MESSAGE("caplet market vols: " << std::fixed << |
| 411 | std::setprecision(4) << io::sequence(capletVols_)); |
| 412 | BOOST_TEST_MESSAGE("caplet0Swapt1Prior: " << caplet0Swaption1Priority); |
| 413 | } |
| 414 | |
| 415 | // calibrate |
| 416 | Natural maxUnperiodicIterations = 10; |
| 417 | Real toleranceUnperiodic = 1e-5; // i.e. 1 bp |
| 418 | Natural max1dIterations = 100; |
| 419 | Real tolerance1d = 1e-8; |
| 420 | Size maxPeriodIterations = 30; |
| 421 | Real periodTolerance = 1e-5; |
| 422 | |
| 423 | std::vector<Matrix> swapPseudoRoots; |
| 424 | Real deformationSize; |
| 425 | Real totalSwaptionError; |
| 426 | std::vector<Real> finalScales; //scalings used for matching |
| 427 | Size iterationsDone; // number of period iteratations done |
| 428 | Real errorImprovement; // improvement in error for last iteration |
| 429 | Matrix modelSwaptionVolsMatrix; |
| 430 | |
| 431 | if (printReport_) { |
| 432 | BOOST_TEST_MESSAGE("numberOfFactors: " << numberOfFactors_); |
| 433 | BOOST_TEST_MESSAGE("maxUnperiodicIterations: " << maxUnperiodicIterations); |
| 434 | BOOST_TEST_MESSAGE("toleranceUnperiodic: " << io::rate(toleranceUnperiodic)); |
| 435 | BOOST_TEST_MESSAGE("max1dIterations: " << max1dIterations); |
| 436 | BOOST_TEST_MESSAGE("tolerance1d: " << io::rate(tolerance1d)); |
| 437 | |
| 438 | } |
| 439 | |
| 440 | /*Integer failures =*/ capletSwaptionPeriodicCalibration( |
| 441 | evolution, |
| 442 | corr, |
| 443 | displacedSwapVariances&: varianceInterpolator, |
| 444 | capletVols: capletVols_, |
| 445 | cs, |
| 446 | displacement: displacement_, |
| 447 | caplet0Swaption1Priority, |
| 448 | numberOfFactors: numberOfFactors_, |
| 449 | period, |
| 450 | max1dIterations, |
| 451 | tolerance1d, |
| 452 | maxUnperiodicIterations, |
| 453 | toleranceUnperiodic, |
| 454 | maxPeriodIterations, |
| 455 | periodTolerance, |
| 456 | deformationSize, |
| 457 | totalSwaptionError, // ? |
| 458 | swapCovariancePseudoRoots&: swapPseudoRoots, // the thing we really want the pseudo root for each time step |
| 459 | finalScales, //scalings used for matching |
| 460 | iterationsDone, // number of period iteratations done |
| 461 | errorImprovement, // improvement in error for last iteration |
| 462 | modelSwaptionVolsMatrix // the swaption vols calibrated to at each step of the iteration |
| 463 | ); |
| 464 | |
| 465 | |
| 466 | ext::shared_ptr<MarketModel> smm(new |
| 467 | PseudoRootFacade(swapPseudoRoots, |
| 468 | rateTimes_, |
| 469 | cs->coterminalSwapRates(), |
| 470 | std::vector<Spread>(numberOfRates, displacement_))); |
| 471 | ext::shared_ptr<MarketModel> flmm(new CotSwapToFwdAdapter(smm)); |
| 472 | Matrix capletTotCovariance = flmm->totalCovariance(endIndex: numberOfRates-1); |
| 473 | |
| 474 | |
| 475 | |
| 476 | |
| 477 | std::vector<Volatility> capletVols(numberOfRates); |
| 478 | for (Size i=0; i<numberOfRates; ++i) { |
| 479 | capletVols[i] = std::sqrt(x: capletTotCovariance[i][i]/rateTimes_[i]); |
| 480 | } |
| 481 | |
| 482 | Real error; |
| 483 | Real capletTolerance = 1e-4; // i.e. 1 bp |
| 484 | |
| 485 | // check caplet fit |
| 486 | for (Size i=0; i<numberOfRates; ++i) { |
| 487 | error = std::fabs(x: capletVols[i]-capletVols_[i]); |
| 488 | if (error>capletTolerance) |
| 489 | BOOST_ERROR("failed to reproduce " << io::ordinal(i+1) << " caplet vol:" |
| 490 | "\n expected: " << io::rate(capletVols_[i]) << |
| 491 | "\n realized: " << io::rate(capletVols[i]) << |
| 492 | "\n percentage error: " << error/capletVols_[i] << |
| 493 | "\n error: " << error << |
| 494 | "\n tolerance: " << capletTolerance); |
| 495 | } |
| 496 | |
| 497 | |
| 498 | |
| 499 | std::vector<Spread> adaptedDisplacements(numberBigRates,displacement_); |
| 500 | ext::shared_ptr<MarketModel> adaptedFlmm(new FwdPeriodAdapter(flmm,period,offset,adaptedDisplacements)); |
| 501 | |
| 502 | ext::shared_ptr<MarketModel> adaptedsmm(new FwdToCotSwapAdapter(adaptedFlmm)); |
| 503 | |
| 504 | // check perfect swaption fit |
| 505 | Real swapTolerance = 2e-5; |
| 506 | |
| 507 | Matrix swapTerminalCovariance(adaptedsmm->totalCovariance(endIndex: adaptedsmm->numberOfSteps()-1)); |
| 508 | |
| 509 | for (Size i=0; i<numberBigRates; ++i) { |
| 510 | Volatility expSwaptionVol = swapVariances[i].totalVolatility(i); |
| 511 | // Real cov = swapTerminalCovariance[i][i]; |
| 512 | Time time = adaptedsmm->evolution().rateTimes()[i]; |
| 513 | Volatility swaptionVol = sqrt(x: swapTerminalCovariance[i][i]/time); |
| 514 | |
| 515 | error = std::fabs(x: swaptionVol-expSwaptionVol); |
| 516 | if (error>swapTolerance) |
| 517 | BOOST_ERROR("failed to reproduce " << io::ordinal(i) << " swaption vol:" |
| 518 | "\n expected: " << io::rate(expSwaptionVol) << |
| 519 | "\n realized: " << io::rate(swaptionVol) << |
| 520 | "\n error: " << error << |
| 521 | "\n tolerance: " << swapTolerance); |
| 522 | } |
| 523 | |
| 524 | |
| 525 | |
| 526 | |
| 527 | |
| 528 | |
| 529 | } |
| 530 | |
| 531 | |
| 532 | void MarketModelSmmCapletHomoCalibrationTest::testSphereCylinder() { |
| 533 | |
| 534 | BOOST_TEST_MESSAGE("Testing sphere-cylinder optimization..." ); |
| 535 | |
| 536 | { |
| 537 | Real R =1.0; |
| 538 | Real S =0.5; |
| 539 | Real alpha=1.5; |
| 540 | Real Z1=1.0/sqrt(x: 3.0); |
| 541 | Real Z2=1.0/sqrt(x: 3.0); |
| 542 | Real Z3=1.0/sqrt(x: 3.0); |
| 543 | |
| 544 | SphereCylinderOptimizer optimizer(R, S, alpha, Z1, Z2, Z3); |
| 545 | Size maxIterations=100; |
| 546 | Real tolerance=1e-8; |
| 547 | Real y1, y2, y3; |
| 548 | |
| 549 | optimizer.findClosest(maxIterations, tolerance, y1, y2, y3); |
| 550 | |
| 551 | Real errorTol = 1e-12; |
| 552 | if ( fabs(x: y1-1.0) > errorTol) |
| 553 | BOOST_ERROR("\n failed to reproduce y1=1: " |
| 554 | << y1 << ", " << y2 << ", " << y3); |
| 555 | |
| 556 | if ( fabs(x: y2-0.0) > errorTol) |
| 557 | BOOST_ERROR("\n failed to reproduce y2=0: " |
| 558 | << y1 << ", " << y2 << ", " << y3); |
| 559 | |
| 560 | if ( fabs(x: y3-0.0) > errorTol) |
| 561 | BOOST_ERROR("\n failed to reproduce y3=0: " |
| 562 | << y1 << ", " <<y2 << ", " << y3); |
| 563 | |
| 564 | |
| 565 | optimizer.findByProjection(y1, y2, y3); |
| 566 | |
| 567 | if ( fabs(x: y1-1.0) > errorTol) |
| 568 | BOOST_ERROR("\nfindByProjection failed to reproduce y1=1: " |
| 569 | << y1 << ", " << y2 << ", " << y3); |
| 570 | |
| 571 | if ( fabs(x: y2-0.0) > errorTol) |
| 572 | BOOST_ERROR("\n findByProjection failed to reproduce y2=0: " |
| 573 | << y1 << ", " << y2 << ", " << y3); |
| 574 | |
| 575 | if ( fabs(x: y3-0.0) > errorTol) |
| 576 | BOOST_ERROR("\n findByProjection failed to reproduce y3=0: " |
| 577 | << y1 << ", " <<y2 << ", " << y3); |
| 578 | } |
| 579 | |
| 580 | { |
| 581 | Real R =5.0; |
| 582 | Real S =1.0; |
| 583 | Real alpha=1.0; |
| 584 | Real Z1=1.0; |
| 585 | Real Z2=2.0; |
| 586 | Real Z3=sqrt(x: 20.0); |
| 587 | |
| 588 | SphereCylinderOptimizer optimizer(R, S, alpha, Z1, Z2, Z3); |
| 589 | Size maxIterations=100; |
| 590 | Real tolerance=1e-8; |
| 591 | Real y1,y2,y3; |
| 592 | |
| 593 | optimizer.findClosest(maxIterations, tolerance, y1, y2, y3); |
| 594 | |
| 595 | Real errorTol = 1e-4; |
| 596 | if ( fabs(x: y1-1.03306) > errorTol) |
| 597 | BOOST_ERROR("\n failed to reproduce y1=1.03306: " |
| 598 | << y1 << ", " << y2 << ", " << y3); |
| 599 | |
| 600 | if ( fabs(x: y2-0.999453) > errorTol) |
| 601 | BOOST_ERROR("\n failed to reproduce y2=0.999453: " |
| 602 | << y1 << ", " << y2 << ", " << y3); |
| 603 | |
| 604 | if ( fabs(x: y3-4.78893) > errorTol) |
| 605 | BOOST_ERROR("\n failed to reproduce y3=4.78893: " |
| 606 | << y1 << ", " <<y2 << ", " << y3); |
| 607 | |
| 608 | |
| 609 | optimizer.findByProjection(y1, y2, y3); |
| 610 | |
| 611 | if ( fabs(x: y1-1.0) > errorTol) |
| 612 | BOOST_ERROR("\n findByProjection failed to reproduce y1 =1: " |
| 613 | << y1 << " " << y2 << " " << y3); |
| 614 | |
| 615 | if ( fabs(x: y2-1.0) > errorTol) |
| 616 | BOOST_ERROR("\n findByProjection failed to reproduce y2 =1: " |
| 617 | << y1 << " " << y2 << " " << y3); |
| 618 | |
| 619 | if ( fabs(x: y3-sqrt(x: 23.0)) > errorTol) |
| 620 | BOOST_ERROR("\n findByProjection failed to reproduce y3 =sqrt(23): " |
| 621 | << y1 << " " <<y2 << " " << y3); |
| 622 | |
| 623 | } |
| 624 | } |
| 625 | |
| 626 | |
| 627 | // --- Call the desired tests |
| 628 | test_suite* MarketModelSmmCapletHomoCalibrationTest::suite() { |
| 629 | auto* suite = BOOST_TEST_SUITE("SMM Caplet homogeneous calibration test" ); |
| 630 | |
| 631 | suite->add(QUANTLIB_TEST_CASE( |
| 632 | &MarketModelSmmCapletHomoCalibrationTest::testFunction)); |
| 633 | suite->add(QUANTLIB_TEST_CASE( |
| 634 | &MarketModelSmmCapletHomoCalibrationTest::testPeriodFunction)); |
| 635 | |
| 636 | suite->add(QUANTLIB_TEST_CASE( |
| 637 | &MarketModelSmmCapletHomoCalibrationTest::testSphereCylinder)); |
| 638 | |
| 639 | return suite; |
| 640 | } |
| 641 | |