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

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