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) 2006 François du Vignaud
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19 */
20
21 #include "swapforwardmappings.hpp"
22 #include "utilities.hpp"
23 #include <ql/models/marketmodels/swapforwardmappings.hpp>
24 #include <ql/models/marketmodels/correlations/timehomogeneousforwardcorrelation.hpp>
25 #include <ql/models/marketmodels/curvestates/lmmcurvestate.hpp>
26 #include <ql/models/marketmodels/evolutiondescription.hpp>
27 #include <ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp>
28 #include <ql/models/marketmodels/models/flatvol.hpp>
29 #include <ql/models/marketmodels/correlations/expcorrelations.hpp>
30 #include <ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp>
31 #include <ql/models/marketmodels/products/multistep/multistepcoterminalswaptions.hpp>
32 #include <ql/models/marketmodels/accountingengine.hpp>
33 #include <ql/models/marketmodels/models/cotswaptofwdadapter.hpp>
34 #include <ql/models/marketmodels/curvestates/coterminalswapcurvestate.hpp>
35 #include <ql/time/schedule.hpp>
36 #include <ql/time/daycounters/simpledaycounter.hpp>
37 #include <ql/math/statistics/sequencestatistics.hpp>
38 #include <ql/pricingengines/blackcalculator.hpp>
39
40 #include <ql/models/marketmodels/products/multistep/multistepswaption.hpp>
41
42 #if defined(BOOST_MSVC)
43 #include <float.h>
44 //namespace { unsigned int u = _controlfp(_EM_INEXACT, _MCW_EM); }
45 #endif
46
47 using namespace QuantLib;
48 using namespace boost::unit_test_framework;
49
50 using std::fabs;
51 using std::sqrt;
52
53 #define BEGIN(x) (x+0)
54 #define END(x) (x+LENGTH(x))
55
56 namespace {
57
58 class MarketModelData{
59 public:
60 MarketModelData();
rateTimes()61 const std::vector<Time>& rateTimes(){return rateTimes_;}
forwards()62 const std::vector<Rate>& forwards(){return forwards_;}
volatilities()63 const std::vector<Volatility>& volatilities(){return volatilities_;}
displacements()64 const std::vector<Rate>& displacements(){return displacements_;}
discountFactors()65 const std::vector<DiscountFactor>& discountFactors(){return discountFactors_;}
nbRates() const66 Size nbRates() const { return nbRates_; }
67
68 private:
69 std::vector<Time> rateTimes_, accruals_;
70 std::vector<Rate> forwards_;
71 std::vector<Spread> displacements_;
72 std::vector<Volatility> volatilities_;
73 std::vector<DiscountFactor> discountFactors_;
74 Size nbRates_;
75 };
76
MarketModelData()77 MarketModelData::MarketModelData(){
78 // Times
79 Calendar calendar = NullCalendar();
80 Date todaysDate = Settings::instance().evaluationDate();
81 Date endDate = todaysDate + 9*Years; // change back
82 Schedule dates(todaysDate, endDate, Period(Semiannual),
83 calendar, Following, Following, DateGeneration::Backward, false);
84 nbRates_ = dates.size()-2;
85 rateTimes_ = std::vector<Time>(nbRates_+1);
86 //paymentTimes_ = std::vector<Time>(rateTimes_.size()-1);
87 accruals_ = std::vector<Time>(nbRates_);
88 DayCounter dayCounter = SimpleDayCounter();
89 for (Size i=1; i<nbRates_+2; ++i)
90 rateTimes_[i-1] = dayCounter.yearFraction(todaysDate, dates[i]);
91
92 displacements_ = std::vector<Rate>(nbRates_, .0);
93
94 forwards_ = std::vector<Rate>(nbRates_);
95 discountFactors_ = std::vector<Rate>(nbRates_+1);
96 discountFactors_[0] = 1.0; // .95; fdv1-> WHY ???????
97 for (Size i=0; i<nbRates_; ++i){
98 forwards_[i] = 0.03 + 0.0010*i;
99 accruals_[i] = rateTimes_[i+1] - rateTimes_[i];
100 discountFactors_[i+1] = discountFactors_[i]
101 /(1+forwards_[i]*accruals_[i]);
102 }
103 Volatility mktVols[] = {0.15541283,
104 0.18719678,
105 0.20890740,
106 0.22318179,
107 0.23212717,
108 0.23731450,
109 0.23988649,
110 0.24066384,
111 0.24023111,
112 0.23900189,
113 0.23726699,
114 0.23522952,
115 0.23303022,
116 0.23076564,
117 0.22850101,
118 0.22627951,
119 0.22412881,
120 0.22206569,
121 0.22009939
122 /*
123 0.2,
124 0.2,
125 0.2,
126 0.2,
127 0.2,
128 0.2,
129 0.2,
130 0.2,
131 0.2,
132 0.2,
133 0.2,
134 0.2,
135 0.2,
136 0.2,
137 0.2,
138 0.2,
139 0.2,
140 0.2,
141 0.2,
142 0.2
143 */
144
145 };
146 volatilities_ = std::vector<Volatility>(nbRates_);
147 for (Size i = 0; i < volatilities_.size(); ++i)
148 volatilities_[i] = mktVols[i];//.0;
149 }
150
151 ext::shared_ptr<SequenceStatisticsInc>
simulate(const std::vector<Real> & todaysDiscounts,const ext::shared_ptr<MarketModelEvolver> & evolver,const MarketModelMultiProduct & product)152 simulate(const std::vector<Real>& todaysDiscounts,
153 const ext::shared_ptr<MarketModelEvolver>& evolver,
154 const MarketModelMultiProduct& product) {
155 Size paths_;
156 #ifdef _DEBUG
157 paths_ = 127;// //
158 #else
159 paths_ = 32767; //262144-1; // //; // 2^15-1
160 #endif
161
162 Size initialNumeraire = evolver->numeraires().front();
163 Real initialNumeraireValue = todaysDiscounts[initialNumeraire];
164
165 AccountingEngine engine(evolver, product, initialNumeraireValue);
166 ext::shared_ptr<SequenceStatisticsInc> stats(new
167 SequenceStatisticsInc(product.numberOfProducts()));
168 engine.multiplePathValues(*stats, paths_);
169 return stats;
170 }
171
makeMultiStepCoterminalSwaptions(const std::vector<Time> & rateTimes,Real strike)172 MultiStepCoterminalSwaptions makeMultiStepCoterminalSwaptions(
173 const std::vector<Time>& rateTimes, Real strike ){
174 std::vector<Time> paymentTimes(rateTimes.begin(), rateTimes.end()-1);
175 std::vector<ext::shared_ptr<StrikedTypePayoff> > payoffs(paymentTimes.size());
176 for (Size i = 0; i < payoffs.size(); ++i){
177 payoffs[i] = ext::shared_ptr<StrikedTypePayoff>(new
178 PlainVanillaPayoff(Option::Call, strike));
179 }
180 return MultiStepCoterminalSwaptions (rateTimes,
181 paymentTimes, payoffs);
182
183 }
184
185 }
186
187
testForwardSwapJacobians()188 void SwapForwardMappingsTest::testForwardSwapJacobians()
189 {
190 {
191 BOOST_TEST_MESSAGE("Testing forward-rate coinitial-swap Jacobian...");
192 MarketModelData marketData;
193 const std::vector<Time>& rateTimes = marketData.rateTimes();
194 const std::vector<Rate>& forwards = marketData.forwards();
195 const Size nbRates = marketData.nbRates();
196 LMMCurveState lmmCurveState(rateTimes);
197 lmmCurveState.setOnForwardRates(forwards);
198
199 Real bumpSize = 1e-8;
200
201 std::vector<Rate> bumpedForwards(forwards);
202
203 Matrix coinitialJacobian(nbRates,nbRates);
204
205 for (Size i=0; i < nbRates; ++i)
206 for (Size j=0; j < nbRates; ++j)
207 {
208 bumpedForwards = forwards;
209 bumpedForwards[j]+= bumpSize;
210 lmmCurveState.setOnForwardRates(bumpedForwards);
211 Real upRate = lmmCurveState.cmSwapRate(0,i+1);
212 bumpedForwards[j]-= 2.0*bumpSize;
213 lmmCurveState.setOnForwardRates(bumpedForwards);
214 Real downRate = lmmCurveState.cmSwapRate(0,i+1);
215 Real deriv = (upRate-downRate)/(2.0*bumpSize);
216 coinitialJacobian[i][j] = deriv;
217
218 }
219
220 Matrix modelJacobian(SwapForwardMappings::coinitialSwapForwardJacobian(lmmCurveState));
221
222 Real errorTolerance = 1e-5;
223
224
225 for (Size i=0; i < nbRates; ++i)
226 for (Size j=0; j < nbRates; ++j)
227 if( fabs(modelJacobian[i][j]-coinitialJacobian[i][j]) > errorTolerance)
228 {
229 BOOST_TEST_MESSAGE("rate " << i
230 << ", sensitivity " << j
231 << ", formula value " << modelJacobian[i][j]
232 << " bumping value " << coinitialJacobian[i][j]
233 << "\n");
234
235 BOOST_ERROR("test failed");
236 }
237 }
238
239 {
240
241 BOOST_TEST_MESSAGE("Testing forward-rate constant-maturity swap Jacobian...");
242 MarketModelData marketData;
243 const std::vector<Time>& rateTimes = marketData.rateTimes();
244 const std::vector<Rate>& forwards = marketData.forwards();
245 const Size nbRates = marketData.nbRates();
246 LMMCurveState lmmCurveState(rateTimes);
247 lmmCurveState.setOnForwardRates(forwards);
248
249 Real bumpSize = 1e-8;
250
251 for( Size spanningForwards = 1; spanningForwards < nbRates; ++spanningForwards)
252 {
253
254 std::vector<Rate> bumpedForwards(forwards);
255
256 Matrix cmsJacobian(nbRates,nbRates);
257
258 for (Size i=0; i < nbRates; ++i)
259 for (Size j=0; j < nbRates; ++j)
260 {
261 bumpedForwards = forwards;
262 bumpedForwards[j]+= bumpSize;
263 lmmCurveState.setOnForwardRates(bumpedForwards);
264 Real upRate = lmmCurveState.cmSwapRate(i,spanningForwards);
265 bumpedForwards[j]-= 2.0*bumpSize;
266 lmmCurveState.setOnForwardRates(bumpedForwards);
267 Real downRate = lmmCurveState.cmSwapRate(i,spanningForwards);
268 Real deriv = (upRate-downRate)/(2.0*bumpSize);
269 cmsJacobian[i][j] = deriv;
270
271 }
272
273 Matrix modelJacobian(SwapForwardMappings::cmSwapForwardJacobian(lmmCurveState, spanningForwards));
274
275 Real errorTolerance = 1e-5;
276
277
278 for (Size i=0; i < nbRates; ++i)
279 for (Size j=0; j < nbRates; ++j)
280 if( fabs(modelJacobian[i][j]-cmsJacobian[i][j]) > errorTolerance)
281 {
282 BOOST_TEST_MESSAGE(
283 "rate " << i
284 << ", sensitivity " << j
285 << ", formula value " << modelJacobian[i][j]
286 << " bumping value " << cmsJacobian[i][j]
287 << "\n");
288
289 BOOST_ERROR("test failed");
290
291 }
292 }
293
294 }
295 }
296
297
testForwardCoterminalMappings()298 void SwapForwardMappingsTest::testForwardCoterminalMappings() {
299
300 BOOST_TEST_MESSAGE("Testing forward-rate coterminal-swap mappings...");
301 MarketModelData marketData;
302 const std::vector<Time>& rateTimes = marketData.rateTimes();
303 const std::vector<Rate>& forwards = marketData.forwards();
304 const Size nbRates = marketData.nbRates();
305 LMMCurveState lmmCurveState(rateTimes);
306 lmmCurveState.setOnForwardRates(forwards);
307
308 const Real longTermCorr=0.5;
309 const Real beta = .2;
310 Real strike = .03;
311 MultiStepCoterminalSwaptions product
312 = makeMultiStepCoterminalSwaptions(rateTimes, strike);
313
314 const EvolutionDescription& evolution = product.evolution();
315 const Size numberOfFactors = nbRates;
316 Spread displacement = marketData.displacements().front();
317 Matrix jacobian =
318 SwapForwardMappings::coterminalSwapZedMatrix(
319 lmmCurveState, displacement);
320
321 Matrix correlations = exponentialCorrelations(evolution.rateTimes(),
322 longTermCorr,
323 beta);
324 ext::shared_ptr<PiecewiseConstantCorrelation> corr(new
325 TimeHomogeneousForwardCorrelation(correlations,
326 rateTimes));
327 ext::shared_ptr<MarketModel> smmMarketModel(new
328 FlatVol(marketData.volatilities(),
329 corr,
330 evolution,
331 numberOfFactors,
332 lmmCurveState.coterminalSwapRates(),
333 marketData.displacements()));
334
335 ext::shared_ptr<MarketModel>
336 lmmMarketModel(new CotSwapToFwdAdapter(smmMarketModel));
337
338 SobolBrownianGeneratorFactory generatorFactory(SobolBrownianGenerator::Diagonal);
339 std::vector<Size> numeraires(nbRates,
340 nbRates);
341 ext::shared_ptr<MarketModelEvolver> evolver(new LogNormalFwdRatePc
342 (lmmMarketModel, generatorFactory, numeraires));
343
344 ext::shared_ptr<SequenceStatisticsInc> stats =
345 simulate(marketData.discountFactors(), evolver, product);
346 std::vector<Real> results = stats->mean();
347 std::vector<Real> errors = stats->errorEstimate();
348
349 const std::vector<DiscountFactor>& todaysDiscounts = marketData.discountFactors();
350 const std::vector<Rate>& todaysCoterminalSwapRates = lmmCurveState.coterminalSwapRates();
351 for (Size i=0; i<nbRates; ++i) {
352 const Matrix& cotSwapsCovariance = smmMarketModel->totalCovariance(i);
353 //Matrix cotSwapsCovariance= jacobian * forwardsCovariance * transpose(jacobian);
354 //Time expiry = rateTimes[i];
355 ext::shared_ptr<PlainVanillaPayoff> payoff(
356 new PlainVanillaPayoff(Option::Call, strike+displacement));
357 //const std::vector<Time>& taus = lmmCurveState.rateTaus();
358 Real expectedSwaption = BlackCalculator(payoff,
359 todaysCoterminalSwapRates[i]+displacement,
360 std::sqrt(cotSwapsCovariance[i][i]),
361 lmmCurveState.coterminalSwapAnnuity(i,i) *
362 todaysDiscounts[i]).value();
363 if (fabs(expectedSwaption-results[i]) > 0.0001)
364 BOOST_ERROR(
365 "expected\t" << expectedSwaption <<
366 "\tLMM\t" << results[i]
367 << "\tstdev:\t" << errors[i] <<
368 "\t" <<std::fabs(results[i]- expectedSwaption)/errors[i]);
369 }
370 }
371
testSwaptionImpliedVolatility()372 void SwapForwardMappingsTest::testSwaptionImpliedVolatility()
373 {
374
375 BOOST_TEST_MESSAGE("Testing implied swaption vol in LMM using HW approximation...");
376 MarketModelData marketData;
377 const std::vector<Time>& rateTimes = marketData.rateTimes();
378 const std::vector<Rate>& forwards = marketData.forwards();
379 const Size nbRates = marketData.nbRates();
380 LMMCurveState lmmCurveState(rateTimes);
381 lmmCurveState.setOnForwardRates(forwards);
382
383 const Real longTermCorr=0.5;
384 const Real beta = .2;
385 Real strike = .03;
386
387 for (Size startIndex = 1; startIndex+2 < nbRates; startIndex = startIndex+5)
388 {
389
390 Size endIndex = nbRates-2;
391
392 ext::shared_ptr<StrikedTypePayoff> payoff(new
393 PlainVanillaPayoff(Option::Call, strike));
394 MultiStepSwaption product(rateTimes, startIndex, endIndex,payoff );
395
396 const EvolutionDescription& evolution = product.evolution();
397 const Size numberOfFactors = nbRates;
398 Spread displacement = marketData.displacements().front();
399 Matrix jacobian =
400 SwapForwardMappings::coterminalSwapZedMatrix(
401 lmmCurveState, displacement);
402
403 Matrix correlations = exponentialCorrelations(evolution.rateTimes(),
404 longTermCorr,
405 beta);
406 ext::shared_ptr<PiecewiseConstantCorrelation> corr(new
407 TimeHomogeneousForwardCorrelation(correlations,
408 rateTimes));
409 ext::shared_ptr<MarketModel> lmmMarketModel(new
410 FlatVol(marketData.volatilities(),
411 corr,
412 evolution,
413 numberOfFactors,
414 lmmCurveState.forwardRates(),
415 marketData.displacements()));
416
417
418 SobolBrownianGeneratorFactory generatorFactory(SobolBrownianGenerator::Diagonal);
419 std::vector<Size> numeraires(nbRates,
420 nbRates);
421 ext::shared_ptr<MarketModelEvolver> evolver(new LogNormalFwdRatePc
422 (lmmMarketModel, generatorFactory, numeraires));
423
424 ext::shared_ptr<SequenceStatisticsInc> stats =
425 simulate(marketData.discountFactors(), evolver, product);
426 std::vector<Real> results = stats->mean();
427 std::vector<Real> errors = stats->errorEstimate();
428
429
430 Real estimatedImpliedVol = SwapForwardMappings::swaptionImpliedVolatility(*lmmMarketModel,startIndex,endIndex);
431
432 Real swapRate = lmmCurveState.cmSwapRate(startIndex,endIndex-startIndex);
433 Real swapAnnuity = lmmCurveState.cmSwapAnnuity(startIndex,startIndex,endIndex-startIndex)*marketData.discountFactors()[startIndex];
434
435 ext::shared_ptr<PlainVanillaPayoff> payoffDis( new PlainVanillaPayoff(Option::Call, strike+displacement));
436
437 Real expectedSwaption = BlackCalculator(payoffDis,
438 swapRate+displacement, estimatedImpliedVol *sqrt(rateTimes[startIndex]),
439 swapAnnuity).value();
440
441 Real error = expectedSwaption - results[0];
442 Real errorInSds = error/errors[0];
443 if (fabs(errorInSds) > 3.5 )
444 BOOST_ERROR(
445 "expected\t" << expectedSwaption <<
446 "\tLMM\t" << results[0]
447 << "\tstdev:\t" << errors[0] <<
448 "\t" <<errorInSds);
449 }
450
451 }
452
453
454
suite()455 test_suite* SwapForwardMappingsTest::suite() {
456 test_suite* suite = BOOST_TEST_SUITE("swap-forward mappings tests");
457
458
459 suite->add(QUANTLIB_TEST_CASE(
460 &SwapForwardMappingsTest::testSwaptionImpliedVolatility));
461
462 suite->add(QUANTLIB_TEST_CASE(
463 &SwapForwardMappingsTest::testForwardSwapJacobians));
464 // #if !defined(QL_NO_UBLAS_SUPPORT)
465 // suite->add(QUANTLIB_TEST_CASE(
466 // &SwapForwardMappingsTest::testForwardCoterminalMappings));
467 // #endif
468 return suite;
469 }
470
471