1	/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2	
3	/*!
4	 Copyright (C) 2002, 2003 Sadruddin Rejeb
5	 Copyright (C) 2004 Ferdinando Ametrano
6	 Copyright (C) 2005, 2006, 2007 StatPro Italia srl
7	
8	 This file is part of QuantLib, a free-software/open-source library
9	 for financial quantitative analysts and developers - http://quantlib.org/
10	
11	 QuantLib is free software: you can redistribute it and/or modify it
12	 under the terms of the QuantLib license.  You should have received a
13	 copy of the license along with this program; if not, please email
14	 <quantlib-dev@lists.sf.net>. The license is also available online at
15	 <http://quantlib.org/license.shtml>.
16	
17	 This program is distributed in the hope that it will be useful, but WITHOUT
18	 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19	 FOR A PARTICULAR PURPOSE.  See the license for more details.
20	*/
21	
22	#include <ql/qldefines.hpp>
23	#ifdef BOOST_MSVC
24	#  include <ql/auto_link.hpp>
25	#endif
26	#include <ql/instruments/swaption.hpp>
27	#include <ql/pricingengines/swap/discountingswapengine.hpp>
28	#include <ql/pricingengines/swaption/treeswaptionengine.hpp>
29	#include <ql/pricingengines/swaption/jamshidianswaptionengine.hpp>
30	#include <ql/pricingengines/swaption/g2swaptionengine.hpp>
31	#include <ql/pricingengines/swaption/fdhullwhiteswaptionengine.hpp>
32	#include <ql/pricingengines/swaption/fdg2swaptionengine.hpp>
33	#include <ql/models/shortrate/calibrationhelpers/swaptionhelper.hpp>
34	#include <ql/models/shortrate/onefactormodels/blackkarasinski.hpp>
35	#include <ql/math/optimization/levenbergmarquardt.hpp>
36	#include <ql/indexes/ibor/euribor.hpp>
37	#include <ql/cashflows/coupon.hpp>
38	#include <ql/quotes/simplequote.hpp>
39	#include <ql/termstructures/yield/flatforward.hpp>
40	#include <ql/time/calendars/target.hpp>
41	#include <ql/time/daycounters/thirty360.hpp>
42	#include <ql/utilities/dataformatters.hpp>
43	
44	#include <boost/timer.hpp>
45	#include <iostream>
46	#include <iomanip>
47	
48	using namespace QuantLib;
49	
50	#if defined(QL_ENABLE_SESSIONS)
51	namespace QuantLib {
52	
53	    Integer sessionId() { return 0; }
54	
55	}
56	#endif
57	
58	
59	//Number of swaptions to be calibrated to...
60	
61	Size numRows = 5;
62	Size numCols = 5;
63	
64	Integer swapLenghts[] = {
65	      1,     2,     3,     4,     5};
66	Volatility swaptionVols[] = {
67	  0.1490, 0.1340, 0.1228, 0.1189, 0.1148,
68	  0.1290, 0.1201, 0.1146, 0.1108, 0.1040,
69	  0.1149, 0.1112, 0.1070, 0.1010, 0.0957,
70	  0.1047, 0.1021, 0.0980, 0.0951, 0.1270,
71	  0.1000, 0.0950, 0.0900, 0.1230, 0.1160};
72	
73	void calibrateModel(
74	          const boost::shared_ptr<ShortRateModel>& model,
75	          const std::vector<boost::shared_ptr<CalibrationHelper> >& helpers) {
76	
77	    LevenbergMarquardt om;
78	    model->calibrate(helpers, om,
79	                     EndCriteria(400, 100, 1.0e-8, 1.0e-8, 1.0e-8));
80	
81	    // Output the implied Black volatilities
82	    for (Size i=0; i<numRows; i++) {
83	        Size j = numCols - i -1; // 1x5, 2x4, 3x3, 4x2, 5x1
84	        Size k = i*numCols + j;
85	        Real npv = helpers[i]->modelValue();
86	        Volatility implied = helpers[i]->impliedVolatility(npv, 1e-4,
87	                                                           1000, 0.05, 0.50);
88	        Volatility diff = implied - swaptionVols[k];
89	
90	        std::cout << i+1 << "x" << swapLenghts[j]
91	                  << std::setprecision(5) << std::noshowpos
92	                  << ": model " << std::setw(7) << io::volatility(implied)
93	                  << ", market " << std::setw(7)
94	                  << io::volatility(swaptionVols[k])
95	                  << " (" << std::setw(7) << std::showpos
96	                  << io::volatility(diff) << std::noshowpos << ")\n";
97	    }
98	}
99	
100	int main(int, char* []) {
101	
102	    try {
103	
104	        boost::timer timer;
105	        std::cout << std::endl;
106	
107	        Date todaysDate(15, February, 2002);
108	        Calendar calendar = TARGET();
109	        Date settlementDate(19, February, 2002);
110	        Settings::instance().evaluationDate() = todaysDate;
111	
112	        // flat yield term structure impling 1x5 swap at 5%
113	        boost::shared_ptr<Quote> flatRate(new SimpleQuote(0.04875825));
114	        Handle<YieldTermStructure> rhTermStructure(
115	            boost::shared_ptr<FlatForward>(
116	                      new FlatForward(settlementDate, Handle<Quote>(flatRate),
117	                                      Actual365Fixed())));
118	
119	        // Define the ATM/OTM/ITM swaps
120	        Frequency fixedLegFrequency = Annual;
121	        BusinessDayConvention fixedLegConvention = Unadjusted;
122	        BusinessDayConvention floatingLegConvention = ModifiedFollowing;
123	        DayCounter fixedLegDayCounter = Thirty360(Thirty360::European);
124	        Frequency floatingLegFrequency = Semiannual;
125	        VanillaSwap::Type type = VanillaSwap::Payer;
126	        Rate dummyFixedRate = 0.03;
127	        boost::shared_ptr<IborIndex> indexSixMonths(new
128	            Euribor6M(rhTermStructure));
129	
130	        Date startDate = calendar.advance(settlementDate,1,Years,
131	                                          floatingLegConvention);
132	        Date maturity = calendar.advance(startDate,5,Years,
133	                                         floatingLegConvention);
134	        Schedule fixedSchedule(startDate,maturity,Period(fixedLegFrequency),
135	                               calendar,fixedLegConvention,fixedLegConvention,
136	                               DateGeneration::Forward,false);
137	        Schedule floatSchedule(startDate,maturity,Period(floatingLegFrequency),
138	                               calendar,floatingLegConvention,floatingLegConvention,
139	                               DateGeneration::Forward,false);
140	
141	        boost::shared_ptr<VanillaSwap> swap(new VanillaSwap(
142	            type, 1000.0,
143	            fixedSchedule, dummyFixedRate, fixedLegDayCounter,
144	            floatSchedule, indexSixMonths, 0.0,
145	            indexSixMonths->dayCounter()));
146	        swap->setPricingEngine(boost::shared_ptr<PricingEngine>(
147	                                 new DiscountingSwapEngine(rhTermStructure)));
148	        Rate fixedATMRate = swap->fairRate();
149	        Rate fixedOTMRate = fixedATMRate * 1.2;
150	        Rate fixedITMRate = fixedATMRate * 0.8;
151	
152	        boost::shared_ptr<VanillaSwap> atmSwap(new VanillaSwap(
153	            type, 1000.0,
154	            fixedSchedule, fixedATMRate, fixedLegDayCounter,
155	            floatSchedule, indexSixMonths, 0.0,
156	            indexSixMonths->dayCounter()));
157	        boost::shared_ptr<VanillaSwap> otmSwap(new VanillaSwap(
158	            type, 1000.0,
159	            fixedSchedule, fixedOTMRate, fixedLegDayCounter,
160	            floatSchedule, indexSixMonths, 0.0,
161	            indexSixMonths->dayCounter()));
162	        boost::shared_ptr<VanillaSwap> itmSwap(new VanillaSwap(
163	            type, 1000.0,
164	            fixedSchedule, fixedITMRate, fixedLegDayCounter,
165	            floatSchedule, indexSixMonths, 0.0,
166	            indexSixMonths->dayCounter()));
167	
168	        // defining the swaptions to be used in model calibration
169	        std::vector<Period> swaptionMaturities;
170	        swaptionMaturities.push_back(Period(1, Years));
171	        swaptionMaturities.push_back(Period(2, Years));
172	        swaptionMaturities.push_back(Period(3, Years));
173	        swaptionMaturities.push_back(Period(4, Years));
174	        swaptionMaturities.push_back(Period(5, Years));
175	
176	        std::vector<boost::shared_ptr<CalibrationHelper> > swaptions;
177	
178	        // List of times that have to be included in the timegrid
179	        std::list<Time> times;
180	
181	        Size i;
182	        for (i=0; i<numRows; i++) {
183	            Size j = numCols - i -1; // 1x5, 2x4, 3x3, 4x2, 5x1
184	            Size k = i*numCols + j;
185	            boost::shared_ptr<Quote> vol(new SimpleQuote(swaptionVols[k]));
186	            swaptions.push_back(boost::shared_ptr<CalibrationHelper>(new
187	                SwaptionHelper(swaptionMaturities[i],
188	                               Period(swapLenghts[j], Years),
189	                               Handle<Quote>(vol),
190	                               indexSixMonths,
191	                               indexSixMonths->tenor(),
192	                               indexSixMonths->dayCounter(),
193	                               indexSixMonths->dayCounter(),
194	                               rhTermStructure)));
195	            swaptions.back()->addTimesTo(times);
196	        }
197	
198	        // Building time-grid
199	        TimeGrid grid(times.begin(), times.end(), 30);
200	
201	
202	        // defining the models
203	        boost::shared_ptr<G2> modelG2(new G2(rhTermStructure));
204	        boost::shared_ptr<HullWhite> modelHW(new HullWhite(rhTermStructure));
205	        boost::shared_ptr<HullWhite> modelHW2(new HullWhite(rhTermStructure));
206	        boost::shared_ptr<BlackKarasinski> modelBK(
207	                                        new BlackKarasinski(rhTermStructure));
208	
209	
210	        // model calibrations
211	
212	        std::cout << "G2 (analytic formulae) calibration" << std::endl;
213	        for (i=0; i<swaptions.size(); i++)
214	            swaptions[i]->setPricingEngine(boost::shared_ptr<PricingEngine>(
215	                new G2SwaptionEngine(modelG2, 6.0, 16)));
216	
217	        calibrateModel(modelG2, swaptions);
218	        std::cout << "calibrated to:\n"
219	                  << "a     = " << modelG2->params()[0] << ", "
220	                  << "sigma = " << modelG2->params()[1] << "\n"
221	                  << "b     = " << modelG2->params()[2] << ", "
222	                  << "eta   = " << modelG2->params()[3] << "\n"
223	                  << "rho   = " << modelG2->params()[4]
224	                  << std::endl << std::endl;
225	
226	
227	
228	        std::cout << "Hull-White (analytic formulae) calibration" << std::endl;
229	        for (i=0; i<swaptions.size(); i++)
230	            swaptions[i]->setPricingEngine(boost::shared_ptr<PricingEngine>(
231	                new JamshidianSwaptionEngine(modelHW)));
232	
233	        calibrateModel(modelHW, swaptions);
234	        std::cout << "calibrated to:\n"
235	                  << "a = " << modelHW->params()[0] << ", "
236	                  << "sigma = " << modelHW->params()[1]
237	                  << std::endl << std::endl;
238	
239	        std::cout << "Hull-White (numerical) calibration" << std::endl;
240	        for (i=0; i<swaptions.size(); i++)
241	            swaptions[i]->setPricingEngine(boost::shared_ptr<PricingEngine>(
242	                                     new TreeSwaptionEngine(modelHW2, grid)));
243	
244	        calibrateModel(modelHW2, swaptions);
245	        std::cout << "calibrated to:\n"
246	                  << "a = " << modelHW2->params()[0] << ", "
247	                  << "sigma = " << modelHW2->params()[1]
248	                  << std::endl << std::endl;
249	
250	        std::cout << "Black-Karasinski (numerical) calibration" << std::endl;
251	        for (i=0; i<swaptions.size(); i++)
252	            swaptions[i]->setPricingEngine(boost::shared_ptr<PricingEngine>(
253	                                      new TreeSwaptionEngine(modelBK, grid)));
254	
255	        calibrateModel(modelBK, swaptions);
256	        std::cout << "calibrated to:\n"
257	                  << "a = " << modelBK->params()[0] << ", "
258	                  << "sigma = " << modelBK->params()[1]
259	                  << std::endl << std::endl;
260	
261	
262	        // ATM Bermudan swaption pricing
263	
264	        std::cout << "Payer bermudan swaption "
265	                  << "struck at " << io::rate(fixedATMRate)
266	                  << " (ATM)" << std::endl;
267	
268	        std::vector<Date> bermudanDates;
269	        const std::vector<boost::shared_ptr<CashFlow> >& leg =
270	            swap->fixedLeg();
271	        for (i=0; i<leg.size(); i++) {
272	            boost::shared_ptr<Coupon> coupon =
273	                boost::dynamic_pointer_cast<Coupon>(leg[i]);
274	            bermudanDates.push_back(coupon->accrualStartDate());
275	        }
276	
277	        boost::shared_ptr<Exercise> bermudanExercise(
278	                                         new BermudanExercise(bermudanDates));
279	
280	        Swaption bermudanSwaption(atmSwap, bermudanExercise);
281	
282	        // Do the pricing for each model
283	
284	        // G2 price the European swaption here, it should switch to bermudan
285	        bermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
286	            new TreeSwaptionEngine(modelG2, 50)));
287	        std::cout << "G2 (tree):      " << bermudanSwaption.NPV() << std::endl;
288	        bermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
289	            new FdG2SwaptionEngine(modelG2)));
290	        std::cout << "G2 (fdm) :      " << bermudanSwaption.NPV() << std::endl;
291	
292	        bermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
293	            new TreeSwaptionEngine(modelHW, 50)));
294	        std::cout << "HW (tree):      " << bermudanSwaption.NPV() << std::endl;
295	        bermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
296	            new FdHullWhiteSwaptionEngine(modelHW)));
297	        std::cout << "HW (fdm) :      " << bermudanSwaption.NPV() << std::endl;
298	
299	        bermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
300	            new TreeSwaptionEngine(modelHW2, 50)));
301	        std::cout << "HW (num, tree): " << bermudanSwaption.NPV() << std::endl;
302	        bermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
303	            new FdHullWhiteSwaptionEngine(modelHW2)));
304	        std::cout << "HW (num, fdm) : " << bermudanSwaption.NPV() << std::endl;
305	
306	        bermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
307	            new TreeSwaptionEngine(modelBK, 50)));
308	        std::cout << "BK:             " << bermudanSwaption.NPV() << std::endl;
309	
310	
311	        // OTM Bermudan swaption pricing
312	
313	        std::cout << "Payer bermudan swaption "
314	                  << "struck at " << io::rate(fixedOTMRate)
315	                  << " (OTM)" << std::endl;
316	
317	        Swaption otmBermudanSwaption(otmSwap,bermudanExercise);
318	
319	        // Do the pricing for each model
320	        otmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
321	            new TreeSwaptionEngine(modelG2, 50)));
322	        std::cout << "G2 (tree):       " << otmBermudanSwaption.NPV()
323	                  << std::endl;
324	        otmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
325	            new FdG2SwaptionEngine(modelG2)));
326	        std::cout << "G2 (fdm) :       " << otmBermudanSwaption.NPV()
327	                  << std::endl;
328	
329	        otmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
330	            new TreeSwaptionEngine(modelHW, 50)));
331	        std::cout << "HW (tree):       " << otmBermudanSwaption.NPV()
332	                  << std::endl;
333	        otmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
334	            new FdHullWhiteSwaptionEngine(modelHW)));
335	        std::cout << "HW (fdm) :       " << otmBermudanSwaption.NPV()
336	                  << std::endl;
337	
338	        otmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
339	            new TreeSwaptionEngine(modelHW2, 50)));
340	        std::cout << "HW (num, tree):  " << otmBermudanSwaption.NPV()
341	                  << std::endl;
342	        otmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
343	            new FdHullWhiteSwaptionEngine(modelHW2)));
344	        std::cout << "HW (num, fdm):   " << otmBermudanSwaption.NPV()
345	                  << std::endl;
346	
347	        otmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
348	            new TreeSwaptionEngine(modelBK, 50)));
349	        std::cout << "BK:              " << otmBermudanSwaption.NPV()
350	                  << std::endl;
351	
352	
353	        // ITM Bermudan swaption pricing
354	
355	        std::cout << "Payer bermudan swaption "
356	                  << "struck at " << io::rate(fixedITMRate)
357	                  << " (ITM)" << std::endl;
358	
359	        Swaption itmBermudanSwaption(itmSwap,bermudanExercise);
360	
361	        // Do the pricing for each model
362	        itmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
363	            new TreeSwaptionEngine(modelG2, 50)));
364	        std::cout << "G2 (tree):       " << itmBermudanSwaption.NPV()
365	                  << std::endl;
366	        itmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
367	            new FdG2SwaptionEngine(modelG2)));
368	        std::cout << "G2 (fdm) :       " << itmBermudanSwaption.NPV()
369	                  << std::endl;
370	
371	        itmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
372	            new TreeSwaptionEngine(modelHW, 50)));
373	        std::cout << "HW (tree):       " << itmBermudanSwaption.NPV()
374	                  << std::endl;
375	        itmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
376	            new FdHullWhiteSwaptionEngine(modelHW)));
377	        std::cout << "HW (fdm) :       " << itmBermudanSwaption.NPV()
378	                  << std::endl;
379	
380	        itmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
381	            new TreeSwaptionEngine(modelHW2, 50)));
382	        std::cout << "HW (num, tree):  " << itmBermudanSwaption.NPV()
383	                  << std::endl;
384	        itmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
385	            new FdHullWhiteSwaptionEngine(modelHW2)));
386	        std::cout << "HW (num, fdm) :  " << itmBermudanSwaption.NPV()
387	                  << std::endl;
388	
389	        itmBermudanSwaption.setPricingEngine(boost::shared_ptr<PricingEngine>(
390	            new TreeSwaptionEngine(modelBK, 50)));
391	        std::cout << "BK:              " << itmBermudanSwaption.NPV()
392	                  << std::endl;
393	
394	        double seconds = timer.elapsed();
395	        Integer hours = int(seconds/3600);
396	        seconds -= hours * 3600;
397	        Integer minutes = int(seconds/60);
398	        seconds -= minutes * 60;
399	        std::cout << " \nRun completed in ";
400	        if (hours > 0)
401	            std::cout << hours << " h ";
402	        if (hours > 0 || minutes > 0)
403	            std::cout << minutes << " m ";
404	        std::cout << std::fixed << std::setprecision(0)
405	                  << seconds << " s\n" << std::endl;
406	
407	        return 0;
408	    } catch (std::exception& e) {
409	        std::cerr << e.what() << std::endl;
410	        return 1;
411	    } catch (...) {
412	        std::cerr << "unknown error" << std::endl;
413	        return 1;
414	    }
415	}