From afbc8d132c9bc2a4e04df382ff75194918a1eb38 Mon Sep 17 00:00:00 2001 From: klausspanderen Date: Thu, 13 May 2010 22:39:57 +0000 Subject: [PATCH 001/410] fixed equal_with bug git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17297 8618b1d8-e22c-0410-b026-96a5dba3e089 --- ql/methods/montecarlo/lsmbasissystem.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ql/methods/montecarlo/lsmbasissystem.cpp b/ql/methods/montecarlo/lsmbasissystem.cpp index 8d27c5cfc0..352354f60a 100644 --- a/ql/methods/montecarlo/lsmbasissystem.cpp +++ b/ql/methods/montecarlo/lsmbasissystem.cpp @@ -152,7 +152,8 @@ namespace QuantLib { // find duplicates for (k=0; k(10*v[k]*QL_EPSILON), + bind(equal_within( + std::fabs(50*v[k]*QL_EPSILON)), v[k], _1) ) == v.begin() + k) { From 3fa07915f3aca455dd0a147cd26bf98e226e3aff Mon Sep 17 00:00:00 2001 From: klausspanderen Date: Tue, 18 May 2010 22:02:05 +0000 Subject: [PATCH 002/410] improved LSM basis system (thanks to Kakhkhor Abdijalilov) git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17298 8618b1d8-e22c-0410-b026-96a5dba3e089 --- ql/methods/montecarlo/lsmbasissystem.cpp | 245 +++++++++-------------- ql/methods/montecarlo/lsmbasissystem.hpp | 14 +- 2 files changed, 106 insertions(+), 153 deletions(-) diff --git a/ql/methods/montecarlo/lsmbasissystem.cpp b/ql/methods/montecarlo/lsmbasissystem.cpp index 352354f60a..64c6b0650f 100644 --- a/ql/methods/montecarlo/lsmbasissystem.cpp +++ b/ql/methods/montecarlo/lsmbasissystem.cpp @@ -2,7 +2,8 @@ /* Copyright (C) 2006 Klaus Spanderen - + Copyright (C) 2010 Kakhkhor Abdijalilov   + This file is part of QuantLib, a free-software/open-source library for financial quantitative analysts and developers - http://quantlib.org/ @@ -20,194 +21,148 @@ /*! \file lsmbasissystem.cpp \brief utility classes for longstaff schwartz early exercise Monte Carlo */ +// lsmbasissystem.hpp -#include -#include #include -#include -#include -#include -#include -#include +#include -using boost::bind; +#include +#include +#include namespace QuantLib { - namespace { + // makes typing a little easier + typedef std::vector > VF_R; + typedef std::vector > VF_A; + typedef std::vector > VV; + Real (GaussianOrthogonalPolynomial::*ptr_w)(Size, Real) const = + &GaussianOrthogonalPolynomial::weightedValue; + + // pow(x, order) class MonomialFct : public std::unary_function { public: - MonomialFct(Size order) : order_(order) {} - + MonomialFct(Size order): order_(order) {} inline Real operator()(const Real x) const { Real ret = 1.0; - for (Size i=0; i Real] functors + to create [Array -> Real] functor */ + class MultiDimFct : public std::unary_function { + public: + MultiDimFct(const VF_R b): b_(b) { + QL_REQUIRE(b_.size()>0, "zero size basis"); + } + inline Real operator()(const Array& a) const { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(b_.size()==a.size(), "wrong argument size"); + #endif + Real ret = b_[0].operator()(a[0]); + for(Size i=1; i > - LsmBasisSystem::pathBasisSystem(Size order, PolynomType polynomType) { + // build order N+1 tuples from order N tuples + VV next_order_tuples(const VV& v) { + const Size order = std::accumulate(v[0].begin(), v[0].end(), 0u); + const Size dim = v[0].size(); + + check_tuples(v, dim, order); + + // the set of unique tuples + std::set > tuples; + std::vector x; + for(Size i=0; i > ret; + VF_R LsmBasisSystem::pathBasisSystem(Size order, PolynomType polyType) { + VF_R ret(order+1); for (Size i=0; i<=order; ++i) { - switch (polynomType) { + switch (polyType) { case Monomial: - ret.push_back(MonomialFct(i)); + ret[i] = MonomialFct(i); break; case Laguerre: - ret.push_back( - bind(&GaussianOrthogonalPolynomial::weightedValue, - GaussLaguerrePolynomial(), i, _1)); + ret[i] = boost::bind(ptr_w, GaussLaguerrePolynomial(), i, _1); break; case Hermite: - ret.push_back( - bind(&GaussianOrthogonalPolynomial::weightedValue, - GaussHermitePolynomial(), i, _1)); + ret[i] = boost::bind(ptr_w, GaussHermitePolynomial(), i, _1); break; case Hyperbolic: - ret.push_back( - bind(&GaussianOrthogonalPolynomial::weightedValue, - GaussHyperbolicPolynomial(), i, _1)); + ret[i] = boost::bind(ptr_w, GaussHyperbolicPolynomial(), i, _1); break; case Legendre: - ret.push_back( - bind(&GaussianOrthogonalPolynomial::weightedValue, - GaussLegendrePolynomial(), i, _1)); + ret[i] = boost::bind(ptr_w, GaussLegendrePolynomial(), i, _1); break; case Chebyshev: - ret.push_back( - bind(&GaussianOrthogonalPolynomial::weightedValue, - GaussChebyshevPolynomial(), i, _1)); + ret[i] = boost::bind(ptr_w, GaussChebyshevPolynomial(), i, _1); break; case Chebyshev2nd: - ret.push_back( - bind(&GaussianOrthogonalPolynomial::weightedValue, - GaussChebyshev2ndPolynomial(), i, _1)); + ret[i] = boost::bind(ptr_w,GaussChebyshev2ndPolynomial(),i, _1); break; default: QL_FAIL("unknown regression type"); } } - return ret; } - - std::vector > - LsmBasisSystem::multiPathBasisSystem(Size dim, Size order, - PolynomType polynomType) { - - const std::vector > b - = pathBasisSystem(order, polynomType); - - std::vector > ret; - ret.push_back(bind(constant(1.0), - bind(f_workaround, _1, 0, dim))); - for (Size i=1; i<=order; ++i) { - const std::vector > a - = w(dim, i, polynomType, b); - - for (std::vector >::const_iterator - iter = a.begin(); iter < a.end(); ++iter) { - ret.push_back(*iter); - } - } - - // remove-o-zap: now remove redundant functions. - // usually we do have a lot of them due to the construction schema. - // We use a more "hands on" method here. - std::deque rm(ret.size(), true); - - Array x(dim), v(ret.size()); - MersenneTwisterUniformRng rng(1234UL); - - for (Size i=0; i<10; ++i) { - Size k; - - // calculate random x vector - for (k=0; k( - std::fabs(50*v[k]*QL_EPSILON)), - v[k], _1) ) - == v.begin() + k) { - - // don't remove this item, it's unique! - rm[k] = false; - } - } - } - - std::vector >::iterator - iter = ret.begin(); - - for (Size i=0; i < rm.size(); ++i) { - if (rm[i]) { - iter = ret.erase(iter); - } - else { - ++iter; + VF_A LsmBasisSystem::multiPathBasisSystem(Size dim, Size order, + PolynomType polyType) { + QL_REQUIRE(dim>0, "zero dimension"); + // get single factor basis + VF_R pathBasis = pathBasisSystem(order, polyType); + VF_A ret; + // 0-th order term + VF_R term(dim, pathBasis[0]); + ret.push_back(MultiDimFct(term)); + // start with all 0 tuple + VV tuples(1, std::vector(dim)); + // add multi-factor terms + for(Size i=1; i<=order; ++i) { + tuples = next_order_tuples(tuples); + // now we have all tuples of order i + // for each tuple add the corresponding term + for(Size j=0; j > - LsmBasisSystem::w(Size dim, Size order, PolynomType polynomType, - const std::vector > & b) { - - std::vector > ret; - - for (Size i=order; i>=1; --i) { - const std::vector > left - = w(dim, order-i, polynomType, b); - - for (Size j=0; j a - = bind(b[i], bind(f_workaround, _1, j, dim)); - - if (i == order) { - ret.push_back(a); - } - else { - // add linear combinations - for (Size j=0; j > - pathBasisSystem(Size order, PolynomType polynomType); - - static std::vector > - multiPathBasisSystem(Size dim, Size order, - PolynomType polynomType); + pathBasisSystem(Size order, PolynomType polyType); - private: static std::vector > - w(Size dim, Size order, PolynomType polynomType, - const std::vector > & basis); + multiPathBasisSystem(Size dim, Size order, PolynomType polyType); }; + } #endif From d62419755d2d9540ba51a9e630c9e5da30088e09 Mon Sep 17 00:00:00 2001 From: klausspanderen Date: Sat, 29 May 2010 10:03:35 +0000 Subject: [PATCH 003/410] changed order of member variable initialization git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17305 8618b1d8-e22c-0410-b026-96a5dba3e089 --- ql/experimental/fx/blackdeltacalculator.cpp | 3 ++- ql/experimental/fx/deltavolquote.cpp | 2 +- ql/experimental/fx/deltavolquote.hpp | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/ql/experimental/fx/blackdeltacalculator.cpp b/ql/experimental/fx/blackdeltacalculator.cpp index 59122fd7b7..8e37c139cd 100644 --- a/ql/experimental/fx/blackdeltacalculator.cpp +++ b/ql/experimental/fx/blackdeltacalculator.cpp @@ -28,8 +28,9 @@ namespace QuantLib { DiscountFactor dDiscount, // domestic discount DiscountFactor fDiscount, // foreign discount Real stdDev): - dt_(dt), spot_(spot), stdDev_(stdDev), ot_(ot), + dt_(dt), ot_(ot), dDiscount_(dDiscount), fDiscount_(fDiscount), + stdDev_(stdDev), spot_(spot), forward_(spot*fDiscount/dDiscount), phi_(Integer(ot)) { QL_REQUIRE(spot_>0.0, diff --git a/ql/experimental/fx/deltavolquote.cpp b/ql/experimental/fx/deltavolquote.cpp index 1892911378..bde5b69590 100644 --- a/ql/experimental/fx/deltavolquote.cpp +++ b/ql/experimental/fx/deltavolquote.cpp @@ -36,7 +36,7 @@ namespace QuantLib { DeltaType deltaType, Time maturity, AtmType atmType) - : vol_(vol), deltaType_(deltaType), atmType_(atmType), maturity_(maturity) { + : vol_(vol), deltaType_(deltaType), maturity_(maturity), atmType_(atmType) { registerWith(vol_); diff --git a/ql/experimental/fx/deltavolquote.hpp b/ql/experimental/fx/deltavolquote.hpp index 29927075c5..f37a1fca27 100644 --- a/ql/experimental/fx/deltavolquote.hpp +++ b/ql/experimental/fx/deltavolquote.hpp @@ -80,8 +80,8 @@ namespace QuantLib { Real delta_; Handle vol_; DeltaType deltaType_; - AtmType atmType_; Time maturity_; + AtmType atmType_; }; } From 81afbe57166f0574bd9c8a10a275c20ac9f7b659 Mon Sep 17 00:00:00 2001 From: lballabio Date: Mon, 7 Jun 2010 12:55:38 +0000 Subject: [PATCH 004/410] Cleaned up a bit doc generation. The basepath is no longer passed from the top-level makefile. The online docs are generated in a folder named after the library version. git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17306 8618b1d8-e22c-0410-b026-96a5dba3e089 --- Docs/Makefile.am | 16 ++++++++++------ Makefile.am | 11 +++++------ 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/Docs/Makefile.am b/Docs/Makefile.am index 88d3f6e160..aa246cdbbd 100644 --- a/Docs/Makefile.am +++ b/Docs/Makefile.am @@ -14,6 +14,8 @@ DOXYGEN_CUSTOM := quantlib.css print.css \ DOXYGEN_INPUT := $(shell find ${top_srcdir} -name *.hpp) \ $(wildcard pages/*.docs) +BASEPATH = @abs_top_srcdir@ + %.gz: % rm -f $@ gzip --best $< @@ -56,16 +58,17 @@ docs-online: .time-stamp-online -e "s/SHOW_FILES = YES/SHOW_FILES = NO/" \ -e "s/INCLUDE_GRAPH = YES/INCLUDE_GRAPH = NO/" \ -e "s/SEARCHENGINE = NO/SEARCHENGINE = YES/" \ - -e "s/HTML_OUTPUT = html/HTML_OUTPUT = html-online/" \ + -e "s/HTML_OUTPUT = html/HTML_OUTPUT = reference-$(VERSION)/" \ -e "s|ql_basepath|${BASEPATH}/|" \ -e "s|ql_version|$(VERSION)|" \ quantlib.doxy > .quantlib.doxy $(DOXYGEN) .quantlib.doxy rm -f .quantlib.doxy - rm -f html-online/*.md5 - cp -f images/*.jpg html-online - cp -f images/*.png html-online - cp -f images/*.ico html-online + rm -f reference-$(VERSION)/*.md5 + rm -f reference-$(VERSION)/*.map + cp -f images/*.jpg reference-$(VERSION) + cp -f images/*.png reference-$(VERSION) + cp -f images/*.ico reference-$(VERSION) touch .time-stamp-online docs-man: .time-stamp-man @@ -83,7 +86,8 @@ docs-man-dist: docs-man docs-clean: rm -Rf man - rm -Rf html-online html + rm -Rf html + rm -Rf reference-* rm -f .time-stamp* diff --git a/Makefile.am b/Makefile.am index 9414999524..f2ef2cd1f5 100644 --- a/Makefile.am +++ b/Makefile.am @@ -41,17 +41,16 @@ benchmark: $(MAKE) -C test-suite benchmark .PHONY: docs docs-man docs-online docs-dist docs-clean -BASEPATH = @abs_top_srcdir@ docs: - BASEPATH=$(BASEPATH) $(MAKE) -C Docs docs-all + $(MAKE) -C Docs docs-all html-local: - BASEPATH=$(BASEPATH) $(MAKE) -C Docs docs-html + $(MAKE) -C Docs docs-html docs-man: - BASEPATH=$(BASEPATH) $(MAKE) -C Docs docs-man + $(MAKE) -C Docs docs-man docs-online: - BASEPATH=$(BASEPATH) $(MAKE) -C Docs docs-online + $(MAKE) -C Docs docs-online docs-dist: - BASEPATH=$(BASEPATH) $(MAKE) -C Docs docs-dist + $(MAKE) -C Docs docs-dist docs-clean: $(MAKE) -C Docs docs-clean From 72eb0cbc6b5aa66e25e0f899ffcd21e788643a89 Mon Sep 17 00:00:00 2001 From: nando Date: Tue, 8 Jun 2010 08:58:49 +0000 Subject: [PATCH 005/410] git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17309 8618b1d8-e22c-0410-b026-96a5dba3e089 --- ql/instruments/makevanillaswap.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ql/instruments/makevanillaswap.cpp b/ql/instruments/makevanillaswap.cpp index b00e20d0ed..f228300cdf 100644 --- a/ql/instruments/makevanillaswap.cpp +++ b/ql/instruments/makevanillaswap.cpp @@ -89,7 +89,7 @@ namespace QuantLib { floatTenor_, floatCalendar_, floatConvention_, floatTerminationDateConvention_, - floatRule_ , floatEndOfMonth_, + floatRule_, floatEndOfMonth_, floatFirstDate_, floatNextToLastDate_); Rate usedFixedRate = fixedRate_; From 19bf416c294514b8fc9bcb1af2f42f1cf8758f4c Mon Sep 17 00:00:00 2001 From: nando Date: Tue, 8 Jun 2010 09:01:00 +0000 Subject: [PATCH 006/410] minor changes git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17310 8618b1d8-e22c-0410-b026-96a5dba3e089 --- ql/math/solvers1d/bisection.hpp | 8 ++++---- ql/math/solvers1d/falseposition.hpp | 16 +++++++--------- ql/math/solvers1d/newton.hpp | 2 +- ql/math/solvers1d/newtonsafe.hpp | 2 +- 4 files changed, 13 insertions(+), 15 deletions(-) diff --git a/ql/math/solvers1d/bisection.hpp b/ql/math/solvers1d/bisection.hpp index 9ae5db1534..08610eb852 100644 --- a/ql/math/solvers1d/bisection.hpp +++ b/ql/math/solvers1d/bisection.hpp @@ -35,7 +35,8 @@ namespace QuantLib { class Bisection : public Solver1D { public: template - Real solveImpl(const F& f, Real xAccuracy) const { + Real solveImpl(const F& f, + Real xAccuracy) const { /* The implementation of the algorithm was inspired by Press, Teukolsky, Vetterling, and Flannery, @@ -58,12 +59,11 @@ namespace QuantLib { dx /= 2.0; xMid=root_+dx; fMid=f(xMid); - evaluationNumber_++; + ++evaluationNumber_; if (fMid <= 0.0) root_=xMid; - if (std::fabs(dx) < xAccuracy || fMid == 0.0) { + if (std::fabs(dx) < xAccuracy || fMid == 0.0) return root_; - } } QL_FAIL("maximum number of function evaluations (" << maxEvaluations_ << ") exceeded"); diff --git a/ql/math/solvers1d/falseposition.hpp b/ql/math/solvers1d/falseposition.hpp index cd894af83f..edea9fd749 100644 --- a/ql/math/solvers1d/falseposition.hpp +++ b/ql/math/solvers1d/falseposition.hpp @@ -35,7 +35,8 @@ namespace QuantLib { class FalsePosition : public Solver1D { public: template - Real solveImpl(const F& f, Real xAccuracy) const { + Real solveImpl(const F& f, + Real xAccuracy) const { /* The implementation of the algorithm was inspired by Press, Teukolsky, Vetterling, and Flannery, @@ -43,8 +44,7 @@ namespace QuantLib { Cambridge University Press */ - Real fl, fh, xl, xh, dx, del, froot; - + Real fl, fh, xl, xh; // Identify the limits so that xl corresponds to the low side if (fxMin_ < 0.0) { xl=xMin_; @@ -57,13 +57,13 @@ namespace QuantLib { xh=xMin_; fh = fxMin_; } - dx=xh-xl; + Real del, froot; while (evaluationNumber_<=maxEvaluations_) { // Increment with respect to latest value - root_=xl+dx*fl/(fl-fh); + root_=xl+(xh-xl)*fl/(fl-fh); froot=f(root_); - evaluationNumber_++; + ++evaluationNumber_; if (froot < 0.0) { // Replace appropriate limit del=xl-root_; xl=root_; @@ -73,11 +73,9 @@ namespace QuantLib { xh=root_; fh=froot; } - dx=xh-xl; // Convergence criterion - if (std::fabs(del) < xAccuracy || froot == 0.0) { + if (std::fabs(del) < xAccuracy || froot == 0.0) return root_; - } } QL_FAIL("maximum number of function evaluations (" diff --git a/ql/math/solvers1d/newton.hpp b/ql/math/solvers1d/newton.hpp index 660e99e63f..175a5db22a 100644 --- a/ql/math/solvers1d/newton.hpp +++ b/ql/math/solvers1d/newton.hpp @@ -67,7 +67,7 @@ namespace QuantLib { return root_; froot = f(root_); dfroot = f.derivative(root_); - evaluationNumber_++; + ++evaluationNumber_; } QL_FAIL("maximum number of function evaluations (" diff --git a/ql/math/solvers1d/newtonsafe.hpp b/ql/math/solvers1d/newtonsafe.hpp index 271f3a6c63..14e20c563b 100644 --- a/ql/math/solvers1d/newtonsafe.hpp +++ b/ql/math/solvers1d/newtonsafe.hpp @@ -91,7 +91,7 @@ namespace QuantLib { return root_; froot = f(root_); dfroot = f.derivative(root_); - evaluationNumber_++; + ++evaluationNumber_; if (froot < 0.0) xl=root_; else From 1cc5d6d760250d0ac12681b233ad6ccc122a4e93 Mon Sep 17 00:00:00 2001 From: nando Date: Wed, 9 Jun 2010 09:28:08 +0000 Subject: [PATCH 007/410] using close(x, y) instead of x==y git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17314 8618b1d8-e22c-0410-b026-96a5dba3e089 --- ql/math/solver1d.hpp | 15 +++++++----- ql/math/solvers1d/bisection.hpp | 9 ++++---- ql/math/solvers1d/brent.hpp | 10 ++++---- ql/math/solvers1d/falseposition.hpp | 27 +++++++++++----------- ql/math/solvers1d/newton.hpp | 8 +++---- ql/math/solvers1d/newtonsafe.hpp | 22 +++++++++--------- ql/math/solvers1d/ridder.hpp | 36 ++++++++++++++--------------- ql/math/solvers1d/secant.hpp | 32 ++++++++++++------------- 8 files changed, 80 insertions(+), 79 deletions(-) diff --git a/ql/math/solver1d.hpp b/ql/math/solver1d.hpp index 0d8a3ae801..61173087b4 100644 --- a/ql/math/solver1d.hpp +++ b/ql/math/solver1d.hpp @@ -24,6 +24,7 @@ #ifndef quantlib_solver1d_hpp #define quantlib_solver1d_hpp +#include #include #include #include @@ -98,9 +99,9 @@ namespace QuantLib { fxMax_ = f(root_); // monotonically crescent bias, as in optionValue(volatility) - if (fxMax_ == 0.0) + if (close(fxMax_,0.0)) return root_; - else if (fxMax_ > 0.0) { + else if (close(fxMax_, 0.0)) { xMin_ = enforceBounds_(root_ - step); fxMin_ = f(xMin_); xMax_ = root_; @@ -114,8 +115,10 @@ namespace QuantLib { evaluationNumber_ = 2; while (evaluationNumber_ <= maxEvaluations_) { if (fxMin_*fxMax_ <= 0.0) { - if (fxMin_ == 0.0) return xMin_; - if (fxMax_ == 0.0) return xMax_; + if (close(fxMin_, 0.0)) + return xMin_; + if (close(fxMax_, 0.0)) + return xMax_; root_ = (xMax_+xMin_)/2.0; return this->impl().solveImpl(f, accuracy); } @@ -183,11 +186,11 @@ namespace QuantLib { << ") > enforced hi bound (" << upperBound_ << ")"); fxMin_ = f(xMin_); - if (fxMin_ == 0.0) + if (close(fxMin_, 0.0)) return xMin_; fxMax_ = f(xMax_); - if (fxMax_ == 0.0) + if (close(fxMax_, 0.0)) return xMax_; evaluationNumber_ = 2; diff --git a/ql/math/solvers1d/bisection.hpp b/ql/math/solvers1d/bisection.hpp index 08610eb852..843b2d99e9 100644 --- a/ql/math/solvers1d/bisection.hpp +++ b/ql/math/solvers1d/bisection.hpp @@ -57,12 +57,12 @@ namespace QuantLib { while (evaluationNumber_<=maxEvaluations_) { dx /= 2.0; - xMid=root_+dx; - fMid=f(xMid); + xMid = root_+dx; + fMid = f(xMid); ++evaluationNumber_; if (fMid <= 0.0) - root_=xMid; - if (std::fabs(dx) < xAccuracy || fMid == 0.0) + root_ = xMid; + if (std::fabs(dx) < xAccuracy || (close(fMid, 0.0))) return root_; } QL_FAIL("maximum number of function evaluations (" @@ -72,5 +72,4 @@ namespace QuantLib { } - #endif diff --git a/ql/math/solvers1d/brent.hpp b/ql/math/solvers1d/brent.hpp index 8c70aa914d..f81c1587c3 100644 --- a/ql/math/solvers1d/brent.hpp +++ b/ql/math/solvers1d/brent.hpp @@ -35,7 +35,8 @@ namespace QuantLib { class Brent : public Solver1D { public: template - Real solveImpl(const F& f, Real xAccuracy) const { + Real solveImpl(const F& f, + Real xAccuracy) const { /* The implementation of the algorithm was inspired by Press, Teukolsky, Vetterling, and Flannery, @@ -70,14 +71,14 @@ namespace QuantLib { // Convergence check xAcc1=2.0*QL_EPSILON*std::fabs(root_)+0.5*xAccuracy; xMid=(xMax_-root_)/2.0; - if (std::fabs(xMid) <= xAcc1 || froot == 0.0) + if (std::fabs(xMid) <= xAcc1 || (close(froot, 0.0))) return root_; if (std::fabs(e) >= xAcc1 && std::fabs(fxMin_) > std::fabs(froot)) { // Attempt inverse quadratic interpolation s=froot/fxMin_; - if (xMin_ == xMax_) { + if (close(xMin_,xMax_)) { p=2.0*xMid*s; q=1.0-s; } else { @@ -109,7 +110,7 @@ namespace QuantLib { else root_ += sign(xAcc1,xMid); froot=f(root_); - evaluationNumber_++; + ++evaluationNumber_; } QL_FAIL("maximum number of function evaluations (" << maxEvaluations_ << ") exceeded"); @@ -122,5 +123,4 @@ namespace QuantLib { } - #endif diff --git a/ql/math/solvers1d/falseposition.hpp b/ql/math/solvers1d/falseposition.hpp index edea9fd749..c5004f2cb7 100644 --- a/ql/math/solvers1d/falseposition.hpp +++ b/ql/math/solvers1d/falseposition.hpp @@ -47,34 +47,34 @@ namespace QuantLib { Real fl, fh, xl, xh; // Identify the limits so that xl corresponds to the low side if (fxMin_ < 0.0) { - xl=xMin_; + xl = xMin_; fl = fxMin_; - xh=xMax_; + xh = xMax_; fh = fxMax_; } else { - xl=xMax_; + xl = xMax_; fl = fxMax_; - xh=xMin_; + xh = xMin_; fh = fxMin_; } Real del, froot; while (evaluationNumber_<=maxEvaluations_) { // Increment with respect to latest value - root_=xl+(xh-xl)*fl/(fl-fh); - froot=f(root_); + root_ = xl+(xh-xl)*fl/(fl-fh); + froot = f(root_); ++evaluationNumber_; if (froot < 0.0) { // Replace appropriate limit - del=xl-root_; - xl=root_; - fl=froot; + del = xl-root_; + xl = root_; + fl = froot; } else { - del=xh-root_; - xh=root_; - fh=froot; + del = xh-root_; + xh = root_; + fh = froot; } // Convergence criterion - if (std::fabs(del) < xAccuracy || froot == 0.0) + if (std::fabs(del) < xAccuracy || (close(froot, 0.0))) return root_; } @@ -85,5 +85,4 @@ namespace QuantLib { } - #endif diff --git a/ql/math/solvers1d/newton.hpp b/ql/math/solvers1d/newton.hpp index 175a5db22a..40922666cf 100644 --- a/ql/math/solvers1d/newton.hpp +++ b/ql/math/solvers1d/newton.hpp @@ -38,7 +38,8 @@ namespace QuantLib { class Newton : public Solver1D { public: template - Real solveImpl(const F& f, Real xAccuracy) const { + Real solveImpl(const F& f, + Real xAccuracy) const { /* The implementation of the algorithm was inspired by Press, Teukolsky, Vetterling, and Flannery, @@ -52,10 +53,10 @@ namespace QuantLib { dfroot = f.derivative(root_); QL_REQUIRE(dfroot != Null(), "Newton requires function's derivative"); - evaluationNumber_++; + ++evaluationNumber_; while (evaluationNumber_<=maxEvaluations_) { - dx=froot/dfroot; + dx = froot/dfroot; root_ -= dx; // jumped out of brackets, switch to NewtonSafe if ((xMin_-root_)*(root_-xMax_) < 0.0) { @@ -77,5 +78,4 @@ namespace QuantLib { } - #endif diff --git a/ql/math/solvers1d/newtonsafe.hpp b/ql/math/solvers1d/newtonsafe.hpp index 14e20c563b..8ee3c9f916 100644 --- a/ql/math/solvers1d/newtonsafe.hpp +++ b/ql/math/solvers1d/newtonsafe.hpp @@ -38,7 +38,8 @@ namespace QuantLib { class NewtonSafe : public Solver1D { public: template - Real solveImpl(const F& f, Real xAccuracy) const { + Real solveImpl(const F& f, + Real xAccuracy) const { /* The implementation of the algorithm was inspired by Press, Teukolsky, Vetterling, and Flannery, @@ -51,26 +52,26 @@ namespace QuantLib { // Orient the search so that f(xl) < 0 if (fxMin_ < 0.0) { - xl=xMin_; - xh=xMax_; + xl = xMin_; + xh = xMax_; } else { - xh=xMin_; - xl=xMax_; + xh = xMin_; + xl = xMax_; } // the "stepsize before last" - dxold=xMax_-xMin_; + dxold = xMax_-xMin_; // it was dxold=std::fabs(xMax_-xMin_); in Numerical Recipes // here (xMax_-xMin_ > 0) is verified in the constructor // and the last step - dx=dxold; + dx = dxold; froot = f(root_); dfroot = f.derivative(root_); QL_REQUIRE(dfroot != Null(), "NewtonSafe requires function's derivative"); - evaluationNumber_++; + ++evaluationNumber_; while (evaluationNumber_<=maxEvaluations_) { // Bisect if (out of range || not decreasing fast enough) @@ -82,8 +83,8 @@ namespace QuantLib { dx = (xh-xl)/2.0; root_=xl+dx; } else { - dxold=dx; - dx=froot/dfroot; + dxold = dx; + dx = froot/dfroot; root_ -= dx; } // Convergence criterion @@ -105,5 +106,4 @@ namespace QuantLib { } - #endif diff --git a/ql/math/solvers1d/ridder.hpp b/ql/math/solvers1d/ridder.hpp index 9be100c7d4..538f75dd79 100644 --- a/ql/math/solvers1d/ridder.hpp +++ b/ql/math/solvers1d/ridder.hpp @@ -35,7 +35,8 @@ namespace QuantLib { class Ridder : public Solver1D { public: template - Real solveImpl(const F& f, Real xAcc) const { + Real solveImpl(const F& f, + Real xAcc) const { /* The implementation of the algorithm was inspired by Press, Teukolsky, Vetterling, and Flannery, @@ -54,12 +55,12 @@ namespace QuantLib { root_ = QL_MIN_REAL; while (evaluationNumber_<=maxEvaluations_) { - xMid=0.5*(xMin_+xMax_); + xMid = 0.5*(xMin_+xMax_); // First of two function evaluations per iteraton - fxMid=f(xMid); - evaluationNumber_++; + fxMid = f(xMid); + ++evaluationNumber_; s = std::sqrt(fxMid*fxMid-fxMin_*fxMax_); - if (s == 0.0) + if (close(s, 0.0)) return root_; // Updating formula nextRoot = xMid + (xMid - xMin_) * @@ -67,25 +68,25 @@ namespace QuantLib { if (std::fabs(nextRoot-root_) <= xAccuracy) return root_; - root_=nextRoot; + root_ = nextRoot; // Second of two function evaluations per iteration - froot=f(root_); - evaluationNumber_++; - if (froot == 0.0) + froot = f(root_); + ++evaluationNumber_; + if (close(froot, 0.0)) return root_; // Bookkeeping to keep the root bracketed on next iteration if (sign(fxMid,froot) != fxMid) { - xMin_=xMid; - fxMin_=fxMid; - xMax_=root_; - fxMax_=froot; + xMin_ = xMid; + fxMin_ = fxMid; + xMax_ = root_; + fxMax_ = froot; } else if (sign(fxMin_,froot) != fxMin_) { - xMax_=root_; - fxMax_=froot; + xMax_ = root_; + fxMax_ = froot; } else if (sign(fxMax_,froot) != fxMax_) { - xMin_=root_; - fxMin_=froot; + xMin_ = root_; + fxMin_ = froot; } else { QL_FAIL("never get here."); } @@ -104,5 +105,4 @@ namespace QuantLib { } - #endif diff --git a/ql/math/solvers1d/secant.hpp b/ql/math/solvers1d/secant.hpp index d65a32cef1..06d5e26d24 100644 --- a/ql/math/solvers1d/secant.hpp +++ b/ql/math/solvers1d/secant.hpp @@ -35,7 +35,8 @@ namespace QuantLib { class Secant : public Solver1D { public: template - Real solveImpl(const F& f, Real xAccuracy) const { + Real solveImpl(const F& f, + Real xAccuracy) const { /* The implementation of the algorithm was inspired by Press, Teukolsky, Vetterling, and Flannery, @@ -48,24 +49,24 @@ namespace QuantLib { // Pick the bound with the smaller function value // as the most recent guess if (std::fabs(fxMin_) < std::fabs(fxMax_)) { - root_=xMin_; - froot=fxMin_; - xl=xMax_; - fl=fxMax_; + root_ = xMin_; + froot = fxMin_; + xl = xMax_; + fl = fxMax_; } else { - root_=xMax_; - froot=fxMax_; - xl=xMin_; - fl=fxMin_; + root_ = xMax_; + froot = fxMax_; + xl = xMin_; + fl = fxMin_; } while (evaluationNumber_<=maxEvaluations_) { - dx=(xl-root_)*froot/(froot-fl); - xl=root_; - fl=froot; + dx = (xl-root_)*froot/(froot-fl); + xl = root_; + fl = froot; root_ += dx; - froot=f(root_); - evaluationNumber_++; - if (std::fabs(dx) < xAccuracy || froot == 0.0) + froot = f(root_); + ++evaluationNumber_; + if (std::fabs(dx) < xAccuracy || (close(froot, 0.0))) return root_; } QL_FAIL("maximum number of function evaluations (" @@ -75,5 +76,4 @@ namespace QuantLib { } - #endif From 8dcfe6338e964c02086237bd9d63691800418e09 Mon Sep 17 00:00:00 2001 From: klausspanderen Date: Sat, 26 Jun 2010 10:46:03 +0000 Subject: [PATCH 008/410] added analytic pricing engine for the piecewise constant time dependent heston model (incl. calibration model) git-svn-id: https://quantlib.svn.sourceforge.net/svnroot/quantlib/trunk/QuantLib@17315 8618b1d8-e22c-0410-b026-96a5dba3e089 --- ql/models/equity/Makefile.am | 6 +- .../piecewisetimedependenthestonmodel.cpp | 66 ++++ .../piecewisetimedependenthestonmodel.hpp | 85 +++++ ql/pricingengines/vanilla/Makefile.am | 2 + .../vanilla/analyticptdhestonengine.cpp | 213 +++++++++++++ .../vanilla/analyticptdhestonengine.hpp | 81 +++++ test-suite/hestonmodel.cpp | 290 ++++++++++++++---- test-suite/hestonmodel.hpp | 4 +- 8 files changed, 679 insertions(+), 68 deletions(-) create mode 100644 ql/models/equity/piecewisetimedependenthestonmodel.cpp create mode 100644 ql/models/equity/piecewisetimedependenthestonmodel.hpp create mode 100644 ql/pricingengines/vanilla/analyticptdhestonengine.cpp create mode 100644 ql/pricingengines/vanilla/analyticptdhestonengine.hpp diff --git a/ql/models/equity/Makefile.am b/ql/models/equity/Makefile.am index a24ab77646..13b2230948 100644 --- a/ql/models/equity/Makefile.am +++ b/ql/models/equity/Makefile.am @@ -7,13 +7,15 @@ this_include_HEADERS = \ batesmodel.hpp \ gjrgarchmodel.hpp \ hestonmodel.hpp \ - hestonmodelhelper.hpp + hestonmodelhelper.hpp \ + piecewisetimedependenthestonmodel.hpp libEquityModels_la_SOURCES = \ batesmodel.cpp \ gjrgarchmodel.cpp \ hestonmodel.cpp \ - hestonmodelhelper.cpp + hestonmodelhelper.cpp \ + piecewisetimedependenthestonmodel.cpp noinst_LTLIBRARIES = libEquityModels.la diff --git a/ql/models/equity/piecewisetimedependenthestonmodel.cpp b/ql/models/equity/piecewisetimedependenthestonmodel.cpp new file mode 100644 index 0000000000..9305744151 --- /dev/null +++ b/ql/models/equity/piecewisetimedependenthestonmodel.cpp @@ -0,0 +1,66 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2010 Klaus Spanderen + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +#include +#include + +namespace QuantLib { + + PiecewiseTimeDependentHestonModel::PiecewiseTimeDependentHestonModel( + const Handle& riskFreeRate, + const Handle& dividendYield, + const Handle& s0, + Real v0, + const Parameter& theta, + const Parameter& kappa, + const Parameter& sigma, + const Parameter& rho, + const TimeGrid& timeGrid) + : CalibratedModel(5), + s0_ (s0), + riskFreeRate_ (riskFreeRate), + dividendYield_(dividendYield), + timeGrid_ (timeGrid) { + + arguments_[0] = theta; + arguments_[1] = kappa; + arguments_[2] = sigma; + arguments_[3] = rho; + arguments_[4] = ConstantParameter(v0, PositiveConstraint()); + + registerWith(s0); + registerWith(riskFreeRate); + registerWith(dividendYield); + } + + const TimeGrid& PiecewiseTimeDependentHestonModel::timeGrid() const { + return timeGrid_; + } + + const Handle& + PiecewiseTimeDependentHestonModel::dividendYield() const { + return dividendYield_; + } + + const Handle& + PiecewiseTimeDependentHestonModel::riskFreeRate() const { + return riskFreeRate_; + } +} + diff --git a/ql/models/equity/piecewisetimedependenthestonmodel.hpp b/ql/models/equity/piecewisetimedependenthestonmodel.hpp new file mode 100644 index 0000000000..f65d73c9f5 --- /dev/null +++ b/ql/models/equity/piecewisetimedependenthestonmodel.hpp @@ -0,0 +1,85 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2010 Klaus Spanderen + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +/*! \file piecewisetimedependenthestonmodel.hpp + \brief piecewise constant time dependent Heston-model +*/ + +#ifndef quantlib_piecewise_time_dependent_heston_model_hpp +#define quantlib_piecewise_time_dependent_heston_model_hpp + +#include +#include + +namespace QuantLib { + + //! Piecewise time dependent Heston model + /*! References: + + Heston, Steven L., 1993. A Closed-Form Solution for Options + with Stochastic Volatility with Applications to Bond and + Currency Options. The review of Financial Studies, Volume 6, + Issue 2, 327-343. + + A. Elices, Models with time-dependent parameters using + transform methods: application to Heston’s model, + http://arxiv.org/pdf/0708.2020 + */ + class PiecewiseTimeDependentHestonModel : public CalibratedModel { + public: + PiecewiseTimeDependentHestonModel( + const Handle& riskFreeRate, + const Handle& dividendYield, + const Handle& s0, + Real v0, + const Parameter& theta, + const Parameter& kappa, + const Parameter& sigma, + const Parameter& rho, + const TimeGrid& timeGrid); + + // variance mean version level + Real theta(Time t) const { return arguments_[0](t); } + // variance mean reversion speed + Real kappa(Time t) const { return arguments_[1](t); } + // volatility of the volatility + Real sigma(Time t) const { return arguments_[2](t); } + // correlation + Real rho(Time t) const { return arguments_[3](t); } + // spot variance + Real v0() const { return arguments_[4](0.0); } + // spot + Real s0() const { return s0_->value(); } + + + const TimeGrid& timeGrid() const; + const Handle& dividendYield() const; + const Handle& riskFreeRate() const; + + protected: + const Handle s0_; + const Handle riskFreeRate_; + const Handle dividendYield_; + const TimeGrid timeGrid_; + }; +} + + +#endif + diff --git a/ql/pricingengines/vanilla/Makefile.am b/ql/pricingengines/vanilla/Makefile.am index e0c5ca1770..7511eacfec 100644 --- a/ql/pricingengines/vanilla/Makefile.am +++ b/ql/pricingengines/vanilla/Makefile.am @@ -11,6 +11,7 @@ this_include_HEADERS = \ analyticgjrgarchengine.hpp \ analytichestonengine.hpp \ analytichestonhullwhiteengine.hpp \ + analyticptdhestonengine.hpp \ baroneadesiwhaleyengine.hpp \ batesengine.hpp \ binomialengine.hpp \ @@ -47,6 +48,7 @@ libVanillaEngines_la_SOURCES = \ analyticgjrgarchengine.cpp \ analytichestonengine.cpp \ analytichestonhullwhiteengine.cpp \ + analyticptdhestonengine.cpp \ baroneadesiwhaleyengine.cpp \ batesengine.cpp \ bjerksundstenslandengine.cpp \ diff --git a/ql/pricingengines/vanilla/analyticptdhestonengine.cpp b/ql/pricingengines/vanilla/analyticptdhestonengine.cpp new file mode 100644 index 0000000000..65d129e5d0 --- /dev/null +++ b/ql/pricingengines/vanilla/analyticptdhestonengine.cpp @@ -0,0 +1,213 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2010 Klaus Spanderen + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +/*! \file analyticptdhestonengine.cpp + \brief analytic piecewise time dependent Heston-model engine +*/ + +#include +#include +#include + + +namespace QuantLib { + // helper class for integration + class AnalyticPTDHestonEngine::Fj_Helper + : public std::unary_function { + + public: + Fj_Helper( + const Handle& model, + const AnalyticPTDHestonEngine* const engine, + Time term, Real strike, Size j); + + Real operator()(Real phi) const; + + private: + const Size j_; + const Time term_; + const Real v0_, x_, sx_; + + const AnalyticPTDHestonEngine* const engine_; + + std::vector r_, q_; + const boost::shared_ptr qTS_; + const Handle model_; + + const TimeGrid timeGrid_; + }; + + AnalyticPTDHestonEngine::Fj_Helper::Fj_Helper( + const Handle& model, + const AnalyticPTDHestonEngine* const engine, + Time term, Real strike, Size j) + : j_(j), + term_(term), + + v0_(model->v0()), + x_ (std::log(model->s0())), + sx_(std::log(strike)), + engine_(engine), + r_(model->timeGrid().size()-1), + q_(model->timeGrid().size()-1), + model_(model), + timeGrid_(model->timeGrid()){ + + for (Size i=0; i riskFreeRate() + ->forwardRate(begin, end, Continuous, NoFrequency).rate(); + q_[i] = model->dividendYield() + ->forwardRate(begin, end, Continuous, NoFrequency).rate(); + } + + QL_REQUIRE(term_ < model_->timeGrid().back(), "maturity is too large"); + } + + Real AnalyticPTDHestonEngine::Fj_Helper::operator()(Real phi) const { + + // avoid numeric overflow for phi->0. + // todo: use l'Hospital's rule use to get lim_{phi->0} + phi = std::max(Real(std::numeric_limits::epsilon()), phi); + + std::complex D = 0.0; + std::complex C = 0.0; + + for (Size i=timeGrid_.size(); i > 0; --i) { + const Time begin = timeGrid_[i-1]; + if (begin < term_) { + const Time end = std::min(term_, timeGrid_[i]); + const Time tau = end-begin; + const Time t = 0.5*(end+begin); + + const Real rho = model_->rho(t); + const Real sigma = model_->sigma(t); + const Real kappa = model_->kappa(t); + const Real theta = model_->theta(t); + + const Real sigma2 = sigma*sigma; + const Real t0 = kappa - ((j_== 1)? rho*sigma : 0); + const Real rpsig = rho*sigma*phi; + + const std::complex t1 = t0+std::complex(0, -rpsig); + const std::complex d = std::sqrt(t1*t1 - sigma2*phi + *std::complex(-phi, (j_== 1)? 1 : -1)); + const std::complex g = (t1-d)/(t1+d); + const std::complex gt + = (t1-d - D*sigma2)/(t1+d - D*sigma2); + + D = (t1+d)/sigma2*(g-gt*std::exp(-d*tau)) + /(1.0-gt*std::exp(-d*tau)); + + const std::complex lng + = std::log((1.0 - gt*std::exp(-d*tau))/(1.0 - gt)); + + C =(kappa*theta)/sigma2*((t1-d)*tau-2.0*lng) + + std::complex(0.0, phi*(r_[i-1]-q_[i-1])*tau) + C; + } + } + return std::exp(v0_*D+C+std::complex(0.0, phi*(x_ - sx_))).imag() + /phi; + } + + AnalyticPTDHestonEngine::AnalyticPTDHestonEngine( + const boost::shared_ptr& model, + Size integrationOrder) + : GenericModelEngine(model), + integration_(new AnalyticHestonEngine::Integration( + AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder))) { + } + + AnalyticPTDHestonEngine::AnalyticPTDHestonEngine( + const boost::shared_ptr& model, + Real relTolerance, Size maxEvaluations) + : GenericModelEngine(model), + integration_(new AnalyticHestonEngine::Integration( + AnalyticHestonEngine::Integration::gaussLobatto( + relTolerance, Null(), maxEvaluations))) { + } + + void AnalyticPTDHestonEngine::calculate() const { + // this is an european option pricer + QL_REQUIRE(arguments_.exercise->type() == Exercise::European, + "not an European option"); + + // plain vanilla + boost::shared_ptr payoff = + boost::dynamic_pointer_cast(arguments_.payoff); + QL_REQUIRE(payoff, "non-striked payoff given"); + + const Real v0 = model_->v0(); + const Real spotPrice = model_->s0(); + QL_REQUIRE(spotPrice > 0.0, "negative or null underlying given"); + + const Real strike = payoff->strike(); + const Real term + = model_->riskFreeRate()->dayCounter().yearFraction( + model_->riskFreeRate()->referenceDate(), + arguments_.exercise->lastDate()); + const Real riskFreeDiscount = model_->riskFreeRate()->discount( + arguments_.exercise->lastDate()); + const Real dividendDiscount = model_->dividendYield()->discount( + arguments_.exercise->lastDate()); + + //average values + const TimeGrid& timeGrid = model_->timeGrid(); + const Size n = timeGrid.size()-1; + Real kappaAvg = 0.0, thetaAvg = 0.0, sigmaAvg=0.0, rhoAvg = 0.0; + + for (Size i=1; i <= n; ++i) { + const Time t = 0.5*(timeGrid[i-1] + timeGrid[i]); + kappaAvg += model_->kappa(t); + thetaAvg += model_->theta(t); + sigmaAvg += model_->sigma(t); + rhoAvg += model_->rho(t); + } + kappaAvg/=n; thetaAvg/=n; sigmaAvg/=n; rhoAvg/=n; + + const Real c_inf = std::min(10.0, std::max(0.0001, + std::sqrt(1.0-square()(rhoAvg))/sigmaAvg)) + *(v0 + kappaAvg*thetaAvg*term); + + const Real p1 = integration_->calculate(c_inf, + Fj_Helper(model_, this, term, strike, 1))/M_PI; + + const Real p2 = integration_->calculate(c_inf, + Fj_Helper(model_, this, term, strike, 2))/M_PI; + + switch (payoff->optionType()) + { + case Option::Call: + results_.value = spotPrice*dividendDiscount*(p1+0.5) + - strike*riskFreeDiscount*(p2+0.5); + break; + case Option::Put: + results_.value = spotPrice*dividendDiscount*(p1-0.5) + - strike*riskFreeDiscount*(p2-0.5); + break; + default: + QL_FAIL("unknown option type"); + } + } +} diff --git a/ql/pricingengines/vanilla/analyticptdhestonengine.hpp b/ql/pricingengines/vanilla/analyticptdhestonengine.hpp new file mode 100644 index 0000000000..1809c1f056 --- /dev/null +++ b/ql/pricingengines/vanilla/analyticptdhestonengine.hpp @@ -0,0 +1,81 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2010 Klaus Spanderen + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +/*! \file analyticptdhestonengine.hpp + \brief analytic piecewise time dependent Heston-model engine +*/ + +#ifndef quantlib_analytic_piecewise_time_dependent_heston_engine_hpp +#define quantlib_analytic_piecewise_time_dependent_heston_engine_hpp + +#include +#include +#include + + +namespace QuantLib { + + //! analytic piecewise constant time dependent Heston-model engine + + /*! References: + + Heston, Steven L., 1993. A Closed-Form Solution for Options + with Stochastic Volatility with Applications to Bond and + Currency Options. The review of Financial Studies, Volume 6, + Issue 2, 327-343. + + J. Gatheral, The Volatility Surface: A Practitioner's Guide, + Wiley Finance + + A. Elices, Models with time-dependent parameters using + transform methods: application to Heston’s model, + http://arxiv.org/pdf/0708.2020 + + \ingroup vanillaengines + */ + class AnalyticPTDHestonEngine + : public GenericModelEngine { + public: + // Simple to use constructor: Using adaptive + // Gauss-Lobatto integration and Gatheral's version of complex log. + // Be aware: using a too large number for maxEvaluations might result + // in a stack overflow as the Lobatto integration is a recursive + // algorithm. + AnalyticPTDHestonEngine( + const boost::shared_ptr& model, + Real relTolerance, Size maxEvaluations); + + // Constructor using Laguerre integration + // and Gatheral's version of complex log. + AnalyticPTDHestonEngine( + const boost::shared_ptr& model, + Size integrationOrder = 144); + + void calculate() const; + + private: + class Fj_Helper; + + const boost::shared_ptr integration_; + }; +} + +#endif diff --git a/test-suite/hestonmodel.cpp b/test-suite/hestonmodel.cpp index 7010867a5b..04d0c8d574 100644 --- a/test-suite/hestonmodel.cpp +++ b/test-suite/hestonmodel.cpp @@ -1,7 +1,7 @@ /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ /* - Copyright (C) 2005, 2007, 2009 Klaus Spanderen + Copyright (C) 2005, 2007, 2009, 2010 Klaus Spanderen This file is part of QuantLib, a free-software/open-source library for financial quantitative analysts and developers - http://quantlib.org/ @@ -24,11 +24,13 @@ #include #include #include +#include #include #include #include #include #include +#include #include #include #include @@ -49,6 +51,90 @@ using namespace QuantLib; using namespace boost::unit_test_framework; +namespace { + + struct CalibrationMarketData { + Handle s0; + Handle riskFreeTS, dividendYield; + std::vector > options; + }; + + CalibrationMarketData getDAXCalibrationMarketData() { + /* this example is taken from A. Sepp + Pricing European-Style Options under Jump Diffusion Processes + with Stochstic Volatility: Applications of Fourier Transform + http://math.ut.ee/~spartak/papers/stochjumpvols.pdf + */ + + Date settlementDate(Settings::instance().evaluationDate()); + + DayCounter dayCounter = Actual365Fixed(); + Calendar calendar = TARGET(); + + Integer t[] = { 13, 41, 75, 165, 256, 345, 524, 703 }; + Rate r[] = { 0.0357,0.0349,0.0341,0.0355,0.0359,0.0368,0.0386,0.0401 }; + + std::vector dates; + std::vector rates; + dates.push_back(settlementDate); + rates.push_back(0.0357); + Size i; + for (i = 0; i < 8; ++i) { + dates.push_back(settlementDate + t[i]); + rates.push_back(r[i]); + } + // FLOATING_POINT_EXCEPTION + Handle riskFreeTS( + boost::shared_ptr( + new ZeroCurve(dates, rates, dayCounter))); + + Handle dividendYield( + flatRate(settlementDate, 0.0, dayCounter)); + + Volatility v[] = + { 0.6625,0.4875,0.4204,0.3667,0.3431,0.3267,0.3121,0.3121, + 0.6007,0.4543,0.3967,0.3511,0.3279,0.3154,0.2984,0.2921, + 0.5084,0.4221,0.3718,0.3327,0.3155,0.3027,0.2919,0.2889, + 0.4541,0.3869,0.3492,0.3149,0.2963,0.2926,0.2819,0.2800, + 0.4060,0.3607,0.3330,0.2999,0.2887,0.2811,0.2751,0.2775, + 0.3726,0.3396,0.3108,0.2781,0.2788,0.2722,0.2661,0.2686, + 0.3550,0.3277,0.3012,0.2781,0.2781,0.2661,0.2661,0.2681, + 0.3428,0.3209,0.2958,0.2740,0.2688,0.2627,0.2580,0.2620, + 0.3302,0.3062,0.2799,0.2631,0.2573,0.2533,0.2504,0.2544, + 0.3343,0.2959,0.2705,0.2540,0.2504,0.2464,0.2448,0.2462, + 0.3460,0.2845,0.2624,0.2463,0.2425,0.2385,0.2373,0.2422, + 0.3857,0.2860,0.2578,0.2399,0.2357,0.2327,0.2312,0.2351, + 0.3976,0.2860,0.2607,0.2356,0.2297,0.2268,0.2241,0.2320 }; + + Handle s0(boost::shared_ptr(new SimpleQuote(4468.17))); + Real strike[] = { 3400,3600,3800,4000,4200,4400, + 4500,4600,4800,5000,5200,5400,5600 }; + + std::vector > options; + + for (Size s = 0; s < 13; ++s) { + for (Size m = 0; m < 8; ++m) { + Handle vol(boost::shared_ptr( + new SimpleQuote(v[s*8+m]))); + + Period maturity((int)((t[m]+3)/7.), Weeks); // round to weeks + options.push_back(boost::shared_ptr( + new HestonModelHelper(maturity, calendar, + s0->value(), strike[s], vol, + riskFreeTS, dividendYield, + CalibrationHelper::ImpliedVolError))); + } + } + + CalibrationMarketData marketData + ={ s0, riskFreeTS, dividendYield, options }; + + return marketData; + } + +} + + void HestonModelTest::testBlackCalibration() { BOOST_MESSAGE( "Testing Heston model calibration using a flat volatility surface..."); @@ -148,11 +234,6 @@ void HestonModelTest::testBlackCalibration() { void HestonModelTest::testDAXCalibration() { - /* this example is taken from A. Sepp - Pricing European-Style Options under Jump Diffusion Processes - with Stochstic Volatility: Applications of Fourier Transform - http://math.ut.ee/~spartak/papers/stochjumpvols.pdf - */ BOOST_MESSAGE( "Testing Heston model calibration using DAX volatility data..."); @@ -162,63 +243,14 @@ void HestonModelTest::testDAXCalibration() { Date settlementDate(5, July, 2002); Settings::instance().evaluationDate() = settlementDate; - DayCounter dayCounter = Actual365Fixed(); - Calendar calendar = TARGET(); - - Integer t[] = { 13, 41, 75, 165, 256, 345, 524, 703 }; - Rate r[] = { 0.0357,0.0349,0.0341,0.0355,0.0359,0.0368,0.0386,0.0401 }; - - std::vector dates; - std::vector rates; - dates.push_back(settlementDate); - rates.push_back(0.0357); - Size i; - for (i = 0; i < 8; ++i) { - dates.push_back(settlementDate + t[i]); - rates.push_back(r[i]); - } - // FLOATING_POINT_EXCEPTION - Handle riskFreeTS( - boost::shared_ptr( - new ZeroCurve(dates, rates, dayCounter))); - - Handle dividendTS( - flatRate(settlementDate, 0.0, dayCounter)); - - Volatility v[] = - { 0.6625,0.4875,0.4204,0.3667,0.3431,0.3267,0.3121,0.3121, - 0.6007,0.4543,0.3967,0.3511,0.3279,0.3154,0.2984,0.2921, - 0.5084,0.4221,0.3718,0.3327,0.3155,0.3027,0.2919,0.2889, - 0.4541,0.3869,0.3492,0.3149,0.2963,0.2926,0.2819,0.2800, - 0.4060,0.3607,0.3330,0.2999,0.2887,0.2811,0.2751,0.2775, - 0.3726,0.3396,0.3108,0.2781,0.2788,0.2722,0.2661,0.2686, - 0.3550,0.3277,0.3012,0.2781,0.2781,0.2661,0.2661,0.2681, - 0.3428,0.3209,0.2958,0.2740,0.2688,0.2627,0.2580,0.2620, - 0.3302,0.3062,0.2799,0.2631,0.2573,0.2533,0.2504,0.2544, - 0.3343,0.2959,0.2705,0.2540,0.2504,0.2464,0.2448,0.2462, - 0.3460,0.2845,0.2624,0.2463,0.2425,0.2385,0.2373,0.2422, - 0.3857,0.2860,0.2578,0.2399,0.2357,0.2327,0.2312,0.2351, - 0.3976,0.2860,0.2607,0.2356,0.2297,0.2268,0.2241,0.2320 }; - - Handle s0(boost::shared_ptr(new SimpleQuote(4468.17))); - Real strike[] = { 3400,3600,3800,4000,4200,4400, - 4500,4600,4800,5000,5200,5400,5600 }; - - std::vector > options; + CalibrationMarketData marketData = getDAXCalibrationMarketData(); + + const Handle riskFreeTS = marketData.riskFreeTS; + const Handle dividendTS = marketData.dividendYield; + const Handle s0 = marketData.s0; - for (Size s = 0; s < 13; ++s) { - for (Size m = 0; m < 8; ++m) { - Handle vol(boost::shared_ptr( - new SimpleQuote(v[s*8+m]))); - - Period maturity((int)((t[m]+3)/7.), Weeks); // round to weeks - options.push_back(boost::shared_ptr( - new HestonModelHelper(maturity, calendar, - s0->value(), strike[s], vol, - riskFreeTS, dividendTS, - CalibrationHelper::ImpliedVolError))); - } - } + const std::vector > options + = marketData.options; const Real v0=0.1; const Real kappa=1.0; @@ -234,14 +266,14 @@ void HestonModelTest::testDAXCalibration() { boost::shared_ptr engine( new AnalyticHestonEngine(model, 64)); - for (i = 0; i < options.size(); ++i) + for (Size i = 0; i < options.size(); ++i) options[i]->setPricingEngine(engine); LevenbergMarquardt om(1e-8, 1e-8, 1e-8); model->calibrate(options, om, EndCriteria(400, 40, 1.0e-8, 1.0e-8, 1.0e-8)); Real sse = 0; - for (i = 0; i < 13*8; ++i) { + for (Size i = 0; i < 13*8; ++i) { const Real diff = options[i]->calibrationError()*100.0; sse += diff*diff; } @@ -961,6 +993,131 @@ void HestonModelTest::testMultipleStrikesEngine() { } + +void HestonModelTest::testAnalyticPiecewiseTimeDependent() { + BOOST_MESSAGE("Testing analytic piecewise time dependent Heston prices..."); + + SavedSettings backup; + + Date settlementDate(27, December, 2004); + Settings::instance().evaluationDate() = settlementDate; + DayCounter dayCounter = ActualActual(); + Date exerciseDate(28, March, 2005); + + boost::shared_ptr payoff( + new PlainVanillaPayoff(Option::Call, 1.0)); + boost::shared_ptr exercise(new EuropeanExercise(exerciseDate)); + + std::vector dates; + dates.push_back(settlementDate); dates.push_back(Date(01, January, 2007)); + std::vector irates; + irates.push_back(0.0); irates.push_back(0.2); + Handle riskFreeTS( + boost::shared_ptr( + new ZeroCurve(dates, irates, dayCounter))); + + std::vector qrates; + qrates.push_back(0.0); qrates.push_back(0.3); + Handle dividendTS( + boost::shared_ptr( + new ZeroCurve(dates, qrates, dayCounter))); + + + const Real v0 = 0.1; + Handle s0(boost::shared_ptr(new SimpleQuote(1.0))); + + ConstantParameter theta(0.09, PositiveConstraint()); + ConstantParameter kappa(3.16, PositiveConstraint()); + ConstantParameter sigma(4.40, PositiveConstraint()); + ConstantParameter rho (-0.8, BoundaryConstraint(-1.0, 1.0)); + + boost::shared_ptr model( + new PiecewiseTimeDependentHestonModel(riskFreeTS, dividendTS, + s0, v0, theta, kappa, + sigma, rho, TimeGrid(20.0, 2))); + + VanillaOption option(payoff, exercise); + option.setPricingEngine(boost::shared_ptr( + new AnalyticPTDHestonEngine(model))); + + const Real calculated = option.NPV(); + boost::shared_ptr hestonProcess( + new HestonProcess(riskFreeTS, dividendTS, s0, v0, + kappa(0.0), theta(0.0), sigma(0.0), rho(0.0))); + boost::shared_ptr hestonModel(new HestonModel(hestonProcess)); + option.setPricingEngine(boost::shared_ptr( + new AnalyticHestonEngine(hestonModel))); + + const Real expected = option.NPV(); + + if (std::fabs(calculated-expected) > 1e-12) { + BOOST_FAIL("failed to reproduce heston prices " + << "\n calculated: " << calculated + << "\n expected: " << expected); + } +} + +void HestonModelTest::testDAXCalibrationOfTimeDependentModel() { + BOOST_MESSAGE( + "Testing Time dependent Heston model calibration ..."); + + SavedSettings backup; + + Date settlementDate(5, July, 2002); + Settings::instance().evaluationDate() = settlementDate; + + CalibrationMarketData marketData = getDAXCalibrationMarketData(); + + const Handle riskFreeTS = marketData.riskFreeTS; + const Handle dividendTS = marketData.dividendYield; + const Handle s0 = marketData.s0; + + const std::vector > options + = marketData.options; + + std::vector