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_smmcapletcalibration.hpp"
24#include "utilities.hpp"
25#include <ql/models/marketmodels/correlations/cotswapfromfwdcorrelation.hpp>
26#include <ql/models/marketmodels/correlations/timehomogeneousforwardcorrelation.hpp>
27#include <ql/models/marketmodels/models/piecewiseconstantabcdvariance.hpp>
28#include <ql/models/marketmodels/models/capletcoterminalswaptioncalibration.hpp>
29#include <ql/models/marketmodels/models/cotswaptofwdadapter.hpp>
30#include <ql/models/marketmodels/models/pseudorootfacade.hpp>
31#include <ql/models/marketmodels/products/multistep/multistepcoterminalswaps.hpp>
32#include <ql/models/marketmodels/products/multistep/multistepcoterminalswaptions.hpp>
33#include <ql/models/marketmodels/products/multistep/multistepswap.hpp>
34#include <ql/models/marketmodels/products/multiproductcomposite.hpp>
35#include <ql/models/marketmodels/accountingengine.hpp>
36#include <ql/models/marketmodels/utilities.hpp>
37#include <ql/models/marketmodels/evolvers/lognormalcotswapratepc.hpp>
38#include <ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp>
39#include <ql/models/marketmodels/correlations/expcorrelations.hpp>
40#include <ql/models/marketmodels/models/flatvol.hpp>
41#include <ql/models/marketmodels/models/abcdvol.hpp>
42#include <ql/models/marketmodels/browniangenerators/mtbrowniangenerator.hpp>
43#include <ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp>
44#include <ql/models/marketmodels/swapforwardmappings.hpp>
45#include <ql/models/marketmodels/curvestates/coterminalswapcurvestate.hpp>
46#include <ql/methods/montecarlo/genericlsregression.hpp>
47#include <ql/legacy/libormarketmodels/lmlinexpcorrmodel.hpp>
48#include <ql/legacy/libormarketmodels/lmextlinexpvolmodel.hpp>
49#include <ql/time/schedule.hpp>
50#include <ql/time/calendars/nullcalendar.hpp>
51#include <ql/time/daycounters/simpledaycounter.hpp>
52#include <ql/pricingengines/blackformula.hpp>
53#include <ql/pricingengines/blackcalculator.hpp>
54#include <ql/utilities/dataformatters.hpp>
55#include <ql/math/integrals/segmentintegral.hpp>
56#include <ql/math/statistics/convergencestatistics.hpp>
57#include <ql/math/statistics/sequencestatistics.hpp>
58#include <sstream>
59
60using namespace QuantLib;
61using namespace boost::unit_test_framework;
62
63namespace market_model_smm_caplet_calibration_test {
64
65 Date todaysDate_, startDate_, endDate_;
66 std::vector<Time> rateTimes_;
67 std::vector<Real> accruals_;
68 Calendar calendar_;
69 DayCounter dayCounter_;
70 std::vector<Rate> todaysForwards_, todaysSwaps_;
71 std::vector<Real> coterminalAnnuity_;
72 Size numberOfFactors_;
73 Real alpha_, alphaMax_, alphaMin_;
74 Spread displacement_;
75 std::vector<DiscountFactor> todaysDiscounts_;
76 std::vector<Volatility> swaptionDisplacedVols_, swaptionVols_;
77 std::vector<Volatility> capletDisplacedVols_, capletVols_;
78 Real a_, b_, c_, d_;
79 Real longTermCorrelation_, beta_;
80 Size measureOffset_;
81 unsigned long seed_;
82 Size paths_, trainingPaths_;
83 bool printReport_ = false;
84
85 void setup() {
86
87 // Times
88 calendar_ = NullCalendar();
89 todaysDate_ = Settings::instance().evaluationDate();
90 //startDate = todaysDate + 5*Years;
91 endDate_ = todaysDate_ + 66*Months;
92 Schedule dates(todaysDate_, endDate_, Period(Semiannual),
93 calendar_, Following, Following, DateGeneration::Backward, false);
94 rateTimes_ = std::vector<Time>(dates.size()-1);
95 accruals_ = std::vector<Real>(rateTimes_.size()-1);
96 dayCounter_ = SimpleDayCounter();
97 for (Size i=1; i<dates.size(); ++i)
98 rateTimes_[i-1] = dayCounter_.yearFraction(d1: todaysDate_, d2: dates[i]);
99 for (Size i=1; i<rateTimes_.size(); ++i)
100 accruals_[i-1] = rateTimes_[i] - rateTimes_[i-1];
101
102 // Rates & displacement
103 todaysForwards_ = std::vector<Rate>(accruals_.size());
104 numberOfFactors_ = 3;
105 alpha_ = -0.05;
106 alphaMax_ = 1.0;
107 alphaMin_ = -1.0;
108 displacement_ = 0.0;
109 for (Size i=0; i<todaysForwards_.size(); ++i) {
110 todaysForwards_[i] = 0.03 + 0.0025*i;
111 //todaysForwards_[i] = 0.03;
112 }
113 LMMCurveState curveState_lmm(rateTimes_);
114 curveState_lmm.setOnForwardRates(fwdRates: todaysForwards_);
115 todaysSwaps_ = curveState_lmm.coterminalSwapRates();
116
117 // Discounts
118 todaysDiscounts_ = std::vector<DiscountFactor>(rateTimes_.size());
119 todaysDiscounts_[0] = 0.95;
120 for (Size i=1; i<rateTimes_.size(); ++i)
121 todaysDiscounts_[i] = todaysDiscounts_[i-1] /
122 (1.0+todaysForwards_[i-1]*accruals_[i-1]);
123
124 //// Swaption Volatilities
125 //Volatility mktSwaptionVols[] = {
126 // 0.15541283,
127 // 0.18719678,
128 // 0.20890740,
129 // 0.22318179,
130 // 0.23212717,
131 // 0.23731450,
132 // 0.23988649,
133 // 0.24066384,
134 // 0.24023111,
135 // 0.23900189,
136 // 0.23726699,
137 // 0.23522952,
138 // 0.23303022,
139 // 0.23076564,
140 // 0.22850101,
141 // 0.22627951,
142 // 0.22412881,
143 // 0.22206569,
144 // 0.22009939
145 //};
146
147 //a = -0.0597;
148 //b = 0.1677;
149 //c = 0.5403;
150 //d = 0.1710;
151
152 a_ = 0.0;
153 b_ = 0.17;
154 c_ = 1.0;
155 d_ = 0.10;
156
157 Volatility mktCapletVols[] = {
158 0.1640,
159 0.1740,
160 0.1840,
161 0.1940,
162 0.1840,
163 0.1740,
164 0.1640,
165 0.1540,
166 0.1440,
167 0.1340376439125532
168 };
169
170 //swaptionDisplacedVols = std::vector<Volatility>(todaysSwaps.size());
171 //swaptionVols = std::vector<Volatility>(todaysSwaps.size());
172 //capletDisplacedVols = std::vector<Volatility>(todaysSwaps.size());
173 capletVols_.resize(new_size: todaysSwaps_.size());
174 for (Size i=0; i<todaysSwaps_.size(); i++) {
175 // swaptionDisplacedVols[i] = todaysSwaps[i]*mktSwaptionVols[i]/
176 // (todaysSwaps[i]+displacement);
177 // swaptionVols[i]= mktSwaptionVols[i];
178 // capletDisplacedVols[i] = todaysForwards[i]*mktCapletVols[i]/
179 // (todaysForwards[i]+displacement);
180 capletVols_[i]= mktCapletVols[i];
181 }
182
183 // Cap/Floor Correlation
184 longTermCorrelation_ = 0.5;
185 beta_ = 0.2;
186 measureOffset_ = 5;
187
188 // Monte Carlo
189 seed_ = 42;
190
191#ifdef _DEBUG
192 paths_ = 127;
193 trainingPaths_ = 31;
194#else
195 paths_ = 32767; //262144-1; //; // 2^15-1
196 trainingPaths_ = 8191; // 2^13-1
197#endif
198 }
199
200 enum MarketModelType { ExponentialCorrelationFlatVolatility,
201 ExponentialCorrelationAbcdVolatility/*,
202 CalibratedMM*/
203 };
204
205 enum MeasureType { ProductSuggested, Terminal,
206 MoneyMarket, MoneyMarketPlus };
207
208 enum EvolverType { Ipc, Pc , NormalPc };
209
210}
211
212
213
214void MarketModelSmmCapletCalibrationTest::testFunction() {
215
216 BOOST_TEST_MESSAGE("Testing GHLS caplet calibration "
217 "in a lognormal coterminal swap market model...");
218
219 using namespace market_model_smm_caplet_calibration_test;
220
221 setup();
222
223 Size numberOfRates = todaysForwards_.size();
224
225 EvolutionDescription evolution(rateTimes_);
226 // Size numberOfSteps = evolution.numberOfSteps();
227
228 ext::shared_ptr<PiecewiseConstantCorrelation> fwdCorr(new
229 ExponentialForwardCorrelation(rateTimes_,
230 longTermCorrelation_,
231 beta_));
232
233 ext::shared_ptr<LMMCurveState> cs(new LMMCurveState(rateTimes_));
234 cs->setOnForwardRates(fwdRates: todaysForwards_);
235
236 ext::shared_ptr<PiecewiseConstantCorrelation> corr(new
237 CotSwapFromFwdCorrelation(fwdCorr, *cs, displacement_));
238
239 std::vector<ext::shared_ptr<PiecewiseConstantVariance> >
240 swapVariances(numberOfRates);
241 for (Size i=0; i<numberOfRates; ++i) {
242 swapVariances[i] = ext::shared_ptr<PiecewiseConstantVariance>(new
243 PiecewiseConstantAbcdVariance(a_, b_, c_, d_,
244 i, rateTimes_));
245 }
246
247 // create calibrator
248 std::vector<Real> alpha(numberOfRates, alpha_);
249 bool lowestRoot = true;
250 bool useFullApprox = false;
251 if (printReport_) {
252 BOOST_TEST_MESSAGE("caplet market vols: " << std::fixed <<
253 std::setprecision(4) << io::sequence(capletVols_));
254 BOOST_TEST_MESSAGE("alpha: " << alpha_);
255 BOOST_TEST_MESSAGE("lowestRoot: " << lowestRoot);
256 BOOST_TEST_MESSAGE("useFullApprox: " << useFullApprox);
257 }
258 CTSMMCapletOriginalCalibration calibrator(evolution,
259 corr,
260 swapVariances,
261 capletVols_,
262 cs,
263 displacement_,
264 alpha,
265 lowestRoot,
266 useFullApprox);
267 // calibrate
268 Natural maxIterations = 2;
269 Real capletTolerance = 0.0001;
270 Natural innerMaxIterations = 50;
271 Real innerTolerance = 1e-9;
272 if (printReport_) {
273 BOOST_TEST_MESSAGE("alpha: " << alpha_);
274 BOOST_TEST_MESSAGE("lowestRoot: " << lowestRoot);
275 BOOST_TEST_MESSAGE("useFullApprox: " << useFullApprox);
276 }
277 bool result = calibrator.calibrate(numberOfFactors: numberOfFactors_,
278 maxIterations,
279 tolerance: capletTolerance/10,
280 innerMaxIterations,
281 innerTolerance);
282 if (!result)
283 BOOST_ERROR("calibration failed");
284
285 const std::vector<Matrix>& swapPseudoRoots = calibrator.swapPseudoRoots();
286 ext::shared_ptr<MarketModel> smm(new
287 PseudoRootFacade(swapPseudoRoots,
288 rateTimes_,
289 cs->coterminalSwapRates(),
290 std::vector<Spread>(numberOfRates, displacement_)));
291 CotSwapToFwdAdapter flmm(smm);
292 Matrix capletTotCovariance = flmm.totalCovariance(endIndex: numberOfRates-1);
293
294 std::vector<Volatility> capletVols(numberOfRates);
295 for (Size i=0; i<numberOfRates; ++i) {
296 capletVols[i] = std::sqrt(x: capletTotCovariance[i][i]/rateTimes_[i]);
297 }
298 if (printReport_) {
299 BOOST_TEST_MESSAGE("caplet smm implied vols: " << std::fixed <<
300 std::setprecision(4) << io::sequence(capletVols));
301 BOOST_TEST_MESSAGE("failures: " << calibrator.failures());
302 BOOST_TEST_MESSAGE("deformationSize: " << calibrator.deformationSize());
303 BOOST_TEST_MESSAGE("capletRmsError: " << calibrator.capletRmsError());
304 BOOST_TEST_MESSAGE("capletMaxError: " << calibrator.capletMaxError());
305 BOOST_TEST_MESSAGE("swaptionRmsError: " << calibrator.swaptionRmsError());
306 BOOST_TEST_MESSAGE("swaptionMaxError: " << calibrator.swaptionMaxError());
307 }
308
309 // check perfect swaption fit
310 Real error, swapTolerance = 1e-14;
311 Matrix swapTerminalCovariance(numberOfRates, numberOfRates, 0.0);
312 for (Size i=0; i<numberOfRates; ++i) {
313 Volatility expSwaptionVol = swapVariances[i]->totalVolatility(i);
314 swapTerminalCovariance += swapPseudoRoots[i] * transpose(m: swapPseudoRoots[i]);
315 Volatility swaptionVol = std::sqrt(x: swapTerminalCovariance[i][i]/rateTimes_[i]);
316 error = std::fabs(x: swaptionVol-expSwaptionVol);
317 if (error>swapTolerance)
318 BOOST_ERROR("failed to reproduce " << io::ordinal(i+1) << " swaption vol:"
319 "\n expected: " << io::rate(expSwaptionVol) <<
320 "\n realized: " << io::rate(swaptionVol) <<
321 "\n error: " << error <<
322 "\n tolerance: " << swapTolerance);
323 }
324
325 // check caplet fit
326 for (Size i=0; i<numberOfRates; ++i) {
327 error = std::fabs(x: capletVols[i]-capletVols_[i]);
328 if (error>capletTolerance)
329 BOOST_ERROR("failed to reproduce " << io::ordinal(i+1) << " caplet vol:"
330 "\n expected: " << io::rate(capletVols_[i]) <<
331 "\n realized: " << io::rate(capletVols[i]) <<
332 "\n percentage error: " << error/capletVols_[i] <<
333 "\n error: " << error <<
334 "\n tolerance: " << capletTolerance);
335 }
336}
337
338
339// --- Call the desired tests
340test_suite* MarketModelSmmCapletCalibrationTest::suite() {
341 auto* suite = BOOST_TEST_SUITE("SMM Caplet calibration test");
342
343 suite->add(QUANTLIB_TEST_CASE(&MarketModelSmmCapletCalibrationTest::testFunction));
344
345 return suite;
346}
347

source code of quantlib/test-suite/marketmodel_smmcapletcalibration.cpp