1            /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

2           

3            /*!

4            Copyright (C) 2005, 2006, 2007, 2009 StatPro Italia srl

5           

6            This file is part of QuantLib, a free-software/open-source library

7            for financial quantitative analysts and developers - http://quantlib.org/

8           

9            QuantLib is free software: you can redistribute it and/or modify it

10          under the terms of the QuantLib license.  You should have received a

11          copy of the license along with this program; if not, please email

12          <quantlib-dev@lists.sf.net>. The license is also available online at

13          <http://quantlib.org/license.shtml>.

14         

15          This program is distributed in the hope that it will be useful, but WITHOUT

16          ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS

17          FOR A PARTICULAR PURPOSE.  See the license for more details.

18          */

19         

20          #include <ql/qldefines.hpp>

21          #if !defined(BOOST_ALL_NO_LIB) && defined(BOOST_MSVC)

22          #  include <ql/auto_link.hpp>

23          #endif

24          #include <ql/instruments/vanillaoption.hpp>

25          #include <ql/math/integrals/tanhsinhintegral.hpp>

26          #include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>

27          #include <ql/pricingengines/vanilla/analyticeuropeanvasicekengine.hpp>

28          #include <ql/pricingengines/vanilla/analytichestonengine.hpp>

29          #include <ql/pricingengines/vanilla/baroneadesiwhaleyengine.hpp>

30          #include <ql/pricingengines/vanilla/batesengine.hpp>

31          #include <ql/pricingengines/vanilla/binomialengine.hpp>

32          #include <ql/pricingengines/vanilla/bjerksundstenslandengine.hpp>

33          #include <ql/pricingengines/vanilla/fdblackscholesvanillaengine.hpp>

34          #include <ql/pricingengines/vanilla/integralengine.hpp>

35          #include <ql/pricingengines/vanilla/mcamericanengine.hpp>

36          #include <ql/pricingengines/vanilla/mceuropeanengine.hpp>

37          #include <ql/pricingengines/vanilla/qdfpamericanengine.hpp>

38          #include <ql/time/calendars/target.hpp>

39          #include <ql/utilities/dataformatters.hpp>

40         

41          #include <iostream>

42          #include <iomanip>

43         

44          using namespace QuantLib;

45         

46          int main(int, char* []) {

47         

48              try {

49         

50                  std::cout << std::endl;

51         

52                  // set up dates

53                  Calendar calendar = TARGET();

54                  Date todaysDate(15, May, 1998);

55                  Date settlementDate(17, May, 1998);

56                  Settings::instance().evaluationDate() = todaysDate;

57         

58                  // our options

59                  Option::Type type(Option::Put);

60                  Real underlying = 36;

61                  Real strike = 40;

62                  Spread dividendYield = 0.00;

63                  Rate riskFreeRate = 0.06;

64                  Volatility volatility = 0.20;

65                  Date maturity(17, May, 1999);

66                  DayCounter dayCounter = Actual365Fixed();

67         

68                  std::cout << "Option type = "  << type << std::endl;

69                  std::cout << "Maturity = "        << maturity << std::endl;

70                  std::cout << "Underlying price = "        << underlying << std::endl;

71                  std::cout << "Strike = "                  << strike << std::endl;

72                  std::cout << "Risk-free interest rate = " << io::rate(riskFreeRate)

73                            << std::endl;

74                  std::cout << "Dividend yield = " << io::rate(dividendYield)

75                            << std::endl;

76                  std::cout << "Volatility = " << io::volatility(volatility)

77                            << std::endl;

78                  std::cout << std::endl;

79                  std::string method;

80                  std::cout << std::endl ;

81         

82                  // write column headings

83                  Size widths[] = { 35, 14, 14, 14 };

84                  std::cout << std::setw(widths[0]) << std::left << "Method"

85                            << std::setw(widths[1]) << std::left << "European"

86                            << std::setw(widths[2]) << std::left << "Bermudan"

87                            << std::setw(widths[3]) << std::left << "American"

88                            << std::endl;

89         

90                  std::vector<Date> exerciseDates;

91                  for (Integer i=1; i<=4; i++)

92                      exerciseDates.push_back(settlementDate + 3*i*Months);

93         

94                  auto europeanExercise = ext::make_shared<EuropeanExercise>(maturity);

95         

96                  auto bermudanExercise = ext::make_shared<BermudanExercise>(exerciseDates);

97         

98                  auto americanExercise = ext::make_shared<AmericanExercise>(settlementDate, maturity);

99         

100                Handle<Quote> underlyingH(ext::make_shared<SimpleQuote>(underlying));

101       

102                // bootstrap the yield/dividend/vol curves

103                Handle<YieldTermStructure> flatTermStructure(

104                    ext::make_shared<FlatForward>(settlementDate, riskFreeRate, dayCounter));

105                Handle<YieldTermStructure> flatDividendTS(

106                    ext::make_shared<FlatForward>(settlementDate, dividendYield, dayCounter));

107                Handle<BlackVolTermStructure> flatVolTS(

108                    ext::make_shared<BlackConstantVol>(settlementDate, calendar, volatility,

109                                             dayCounter));

110                auto payoff = ext::make_shared<PlainVanillaPayoff>(type, strike);

111                auto bsmProcess = ext::make_shared<BlackScholesMertonProcess>(

112                        underlyingH, flatDividendTS, flatTermStructure, flatVolTS);

113       

114                // options

115                VanillaOption europeanOption(payoff, europeanExercise);

116                VanillaOption bermudanOption(payoff, bermudanExercise);

117                VanillaOption americanOption(payoff, americanExercise);

118       

119                // Analytic formulas:

120       

121                // Black-Scholes for European

122                method = "Black-Scholes";

123                europeanOption.setPricingEngine(ext::make_shared<AnalyticEuropeanEngine>(bsmProcess));

124                std::cout << std::setw(widths[0]) << std::left << method

125                          << std::fixed

126                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

127                          << std::setw(widths[2]) << std::left << "N/A"

128                          << std::setw(widths[3]) << std::left << "N/A"

129                          << std::endl;

130       

131                //Vasicek rates model for European

132                method = "Black Vasicek Model";

133                Real r0 = riskFreeRate;

134                Real a = 0.3;

135                Real b = 0.3;

136                Real sigma_r = 0.15;

137                Real riskPremium = 0.0;

138                Real correlation = 0.5;

139                auto vasicekProcess = ext::make_shared<Vasicek>(r0, a, b, sigma_r, riskPremium);

140                europeanOption.setPricingEngine(ext::make_shared<AnalyticBlackVasicekEngine>(bsmProcess, vasicekProcess, correlation));

141                std::cout << std::setw(widths[0]) << std::left << method

142                          << std::fixed

143                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

144                          << std::setw(widths[2]) << std::left << "N/A"

145                          << std::setw(widths[3]) << std::left << "N/A"

146                          << std::endl;

147       

148                // semi-analytic Heston for European

149                method = "Heston semi-analytic";

150                auto hestonProcess = ext::make_shared<HestonProcess>(flatTermStructure, flatDividendTS,

151                                      underlyingH, volatility*volatility,

152                                      1.0, volatility*volatility, 0.001, 0.0);

153                auto hestonModel = ext::make_shared<HestonModel>(hestonProcess);

154                europeanOption.setPricingEngine(ext::make_shared<AnalyticHestonEngine>(hestonModel));

155                std::cout << std::setw(widths[0]) << std::left << method

156                          << std::fixed

157                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

158                          << std::setw(widths[2]) << std::left << "N/A"

159                          << std::setw(widths[3]) << std::left << "N/A"

160                          << std::endl;

161       

162                // semi-analytic Bates for European

163                method = "Bates semi-analytic";

164                auto batesProcess = ext::make_shared<BatesProcess>(flatTermStructure, flatDividendTS,

165                                     underlyingH, volatility*volatility,

166                                     1.0, volatility*volatility, 0.001, 0.0,

167                                     1e-14, 1e-14, 1e-14);

168                auto batesModel = ext::make_shared<BatesModel>(batesProcess);

169                europeanOption.setPricingEngine(ext::make_shared<BatesEngine>(batesModel));

170                std::cout << std::setw(widths[0]) << std::left << method

171                          << std::fixed

172                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

173                          << std::setw(widths[2]) << std::left << "N/A"

174                          << std::setw(widths[3]) << std::left << "N/A"

175                          << std::endl;

176       

177                // Barone-Adesi and Whaley approximation for American

178                method = "Barone-Adesi/Whaley";

179                americanOption.setPricingEngine(ext::make_shared<BaroneAdesiWhaleyApproximationEngine>(bsmProcess));

180                std::cout << std::setw(widths[0]) << std::left << method

181                          << std::fixed

182                          << std::setw(widths[1]) << std::left << "N/A"

183                          << std::setw(widths[2]) << std::left << "N/A"

184                          << std::setw(widths[3]) << std::left << americanOption.NPV()

185                          << std::endl;

186       

187                // Bjerksund and Stensland approximation for American

188                method = "Bjerksund/Stensland";

189                americanOption.setPricingEngine(ext::make_shared<BjerksundStenslandApproximationEngine>(bsmProcess));

190                std::cout << std::setw(widths[0]) << std::left << method

191                          << std::fixed

192                          << std::setw(widths[1]) << std::left << "N/A"

193                          << std::setw(widths[2]) << std::left << "N/A"

194                          << std::setw(widths[3]) << std::left << americanOption.NPV()

195                          << std::endl;

196       

197                // QD+ fixed-point engine for American

198                method = "QD+ fixed-point (fast)";

199                americanOption.setPricingEngine(ext::make_shared<QdFpAmericanEngine>

200                                                (bsmProcess, QdFpAmericanEngine::fastScheme()));

201                std::cout << std::setw(widths[0]) << std::left << method

202                          << std::fixed

203                          << std::setw(widths[1]) << std::left << "N/A"

204                          << std::setw(widths[2]) << std::left << "N/A"

205                          << std::setw(widths[3]) << std::left << americanOption.NPV()

206                          << std::endl;

207       

208                method = "QD+ fixed-point (accurate)";

209                americanOption.setPricingEngine(ext::make_shared<QdFpAmericanEngine>

210                                                (bsmProcess, QdFpAmericanEngine::accurateScheme()));

211                std::cout << std::setw(widths[0]) << std::left << method

212                          << std::fixed

213                          << std::setw(widths[1]) << std::left << "N/A"

214                          << std::setw(widths[2]) << std::left << "N/A"

215                          << std::setw(widths[3]) << std::left << americanOption.NPV()

216                          << std::endl;

217       

218                method = "QD+ fixed-point (high precision)";

219                americanOption.setPricingEngine(ext::make_shared<QdFpAmericanEngine>

220                                                (bsmProcess, QdFpAmericanEngine::highPrecisionScheme()));

221                std::cout << std::setw(widths[0]) << std::left << method

222                          << std::fixed

223                          << std::setw(widths[1]) << std::left << "N/A"

224                          << std::setw(widths[2]) << std::left << "N/A"

225                          << std::setw(widths[3]) << std::left << americanOption.NPV()

226                          << std::endl;

227       

228                // Integral

229                method = "Integral";

230                europeanOption.setPricingEngine(ext::make_shared<IntegralEngine>(bsmProcess));

231                std::cout << std::setw(widths[0]) << std::left << method

232                          << std::fixed

233                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

234                          << std::setw(widths[2]) << std::left << "N/A"

235                          << std::setw(widths[3]) << std::left << "N/A"

236                          << std::endl;

237       

238                // Finite differences

239                Size timeSteps = 801;

240                method = "Finite differences";

241                auto fdengine =

242                    ext::make_shared<FdBlackScholesVanillaEngine>(bsmProcess,

243                                                                  timeSteps,

244                                                                  timeSteps-1);

245                europeanOption.setPricingEngine(fdengine);

246                bermudanOption.setPricingEngine(fdengine);

247                americanOption.setPricingEngine(fdengine);

248                std::cout << std::setw(widths[0]) << std::left << method

249                          << std::fixed

250                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

251                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

252                          << std::setw(widths[3]) << std::left << americanOption.NPV()

253                          << std::endl;

254       

255                // Binomial method: Jarrow-Rudd

256                method = "Binomial Jarrow-Rudd";

257                auto jrEngine = ext::make_shared<BinomialVanillaEngine<JarrowRudd>>(bsmProcess, timeSteps);

258                europeanOption.setPricingEngine(jrEngine);

259                bermudanOption.setPricingEngine(jrEngine);

260                americanOption.setPricingEngine(jrEngine);

261                std::cout << std::setw(widths[0]) << std::left << method

262                          << std::fixed

263                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

264                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

265                          << std::setw(widths[3]) << std::left << americanOption.NPV()

266                          << std::endl;

267       

268                // Binomial method: Cox-Ross-Rubinstein

269                method = "Binomial Cox-Ross-Rubinstein";

270               auto crrEngine = ext::make_shared<BinomialVanillaEngine<CoxRossRubinstein>>(bsmProcess, timeSteps);

271                europeanOption.setPricingEngine(crrEngine);

272                bermudanOption.setPricingEngine(crrEngine);

273                americanOption.setPricingEngine(crrEngine);

274                std::cout << std::setw(widths[0]) << std::left << method

275                          << std::fixed

276                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

277                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

278                          << std::setw(widths[3]) << std::left << americanOption.NPV()

279                          << std::endl;

280       

281                // Binomial method: Additive equiprobabilities

282                method = "Additive equiprobabilities";

283                auto aeqpEngine = ext::make_shared<BinomialVanillaEngine<AdditiveEQPBinomialTree>>(bsmProcess, timeSteps);

284                europeanOption.setPricingEngine(aeqpEngine);

285                bermudanOption.setPricingEngine(aeqpEngine);

286                americanOption.setPricingEngine(aeqpEngine);

287                std::cout << std::setw(widths[0]) << std::left << method

288                          << std::fixed

289                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

290                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

291                          << std::setw(widths[3]) << std::left << americanOption.NPV()

292                          << std::endl;

293       

294                // Binomial method: Binomial Trigeorgis

295                method = "Binomial Trigeorgis";

296                auto trigeorgisEngine = ext::make_shared<BinomialVanillaEngine<Trigeorgis>>(bsmProcess, timeSteps);

297                europeanOption.setPricingEngine(trigeorgisEngine);

298                bermudanOption.setPricingEngine(trigeorgisEngine);

299                americanOption.setPricingEngine(trigeorgisEngine);

300                std::cout << std::setw(widths[0]) << std::left << method

301                          << std::fixed

302                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

303                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

304                          << std::setw(widths[3]) << std::left << americanOption.NPV()

305                          << std::endl;

306       

307                // Binomial method: Binomial Tian

308                method = "Binomial Tian";

309                auto tianEngine = ext::make_shared<BinomialVanillaEngine<Tian>>(bsmProcess, timeSteps);

310                europeanOption.setPricingEngine(tianEngine);

311                bermudanOption.setPricingEngine(tianEngine);

312                americanOption.setPricingEngine(tianEngine);

313                std::cout << std::setw(widths[0]) << std::left << method

314                          << std::fixed

315                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

316                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

317                          << std::setw(widths[3]) << std::left << americanOption.NPV()

318                          << std::endl;

319       

320                // Binomial method: Binomial Leisen-Reimer

321                method = "Binomial Leisen-Reimer";

322                auto lrEngine = ext::make_shared<BinomialVanillaEngine<LeisenReimer>>(bsmProcess, timeSteps);

323                europeanOption.setPricingEngine(lrEngine);

324                bermudanOption.setPricingEngine(lrEngine);

325                americanOption.setPricingEngine(lrEngine);

326                std::cout << std::setw(widths[0]) << std::left << method

327                          << std::fixed

328                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

329                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

330                          << std::setw(widths[3]) << std::left << americanOption.NPV()

331                          << std::endl;

332       

333                // Binomial method: Binomial Joshi

334                method = "Binomial Joshi";

335                auto joshiEngine = ext::make_shared<BinomialVanillaEngine<Joshi4>>(bsmProcess, timeSteps);

336                europeanOption.setPricingEngine(joshiEngine);

337                bermudanOption.setPricingEngine(joshiEngine);

338                americanOption.setPricingEngine(joshiEngine);

339                std::cout << std::setw(widths[0]) << std::left << method

340                          << std::fixed

341                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

342                          << std::setw(widths[2]) << std::left << bermudanOption.NPV()

343                          << std::setw(widths[3]) << std::left << americanOption.NPV()

344                          << std::endl;

345       

346                // Monte Carlo Method: MC (crude)

347                timeSteps = 1;

348                method = "MC (crude)";

349                Size mcSeed = 42;

350                auto mcengine1 = MakeMCEuropeanEngine<PseudoRandom>(bsmProcess)

351                    .withSteps(timeSteps)

352                    .withAbsoluteTolerance(0.02)

353                    .withSeed(mcSeed);

354                europeanOption.setPricingEngine(mcengine1);

355                // Real errorEstimate = europeanOption.errorEstimate();

356                std::cout << std::setw(widths[0]) << std::left << method

357                          << std::fixed

358                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

359                          << std::setw(widths[2]) << std::left << "N/A"

360                          << std::setw(widths[3]) << std::left << "N/A"

361                          << std::endl;

362       

363                // Monte Carlo Method: QMC (Sobol)

364                method = "QMC (Sobol)";

365                Size nSamples = 32768;  // 2^15

366       

367                auto mcengine2 = MakeMCEuropeanEngine<LowDiscrepancy>(bsmProcess)

368                    .withSteps(timeSteps)

369                    .withSamples(nSamples);

370                europeanOption.setPricingEngine(mcengine2);

371                std::cout << std::setw(widths[0]) << std::left << method

372                          << std::fixed

373                          << std::setw(widths[1]) << std::left << europeanOption.NPV()

374                          << std::setw(widths[2]) << std::left << "N/A"

375                          << std::setw(widths[3]) << std::left << "N/A"

376                          << std::endl;

377       

378                // Monte Carlo Method: MC (Longstaff Schwartz)

379                method = "MC (Longstaff Schwartz)";

380                auto mcengine3 = MakeMCAmericanEngine<PseudoRandom>(bsmProcess)

381                    .withSteps(100)

382                    .withAntitheticVariate()

383                    .withCalibrationSamples(4096)

384                    .withAbsoluteTolerance(0.02)

385                    .withSeed(mcSeed);

386                americanOption.setPricingEngine(mcengine3);

387                std::cout << std::setw(widths[0]) << std::left << method

388                          << std::fixed

389                          << std::setw(widths[1]) << std::left << "N/A"

390                          << std::setw(widths[2]) << std::left << "N/A"

391                          << std::setw(widths[3]) << std::left << americanOption.NPV()

392                          << std::endl;

393       

394                // End test

395                return 0;

396       

397            } catch (std::exception& e) {

398                std::cerr << e.what() << std::endl;

399                return 1;

400            } catch (...) {

401                std::cerr << "unknown error" << std::endl;

402                return 1;

403            }

404        }