2020
2121#include < ql/experimental/credit/randomdefaultmodel.hpp>
2222#include < ql/math/solvers1d/brent.hpp>
23+ #include < ql/math/solvers1d/bisection.hpp>
2324#include < utility>
2425
2526using namespace std ;
@@ -34,7 +35,8 @@ namespace QuantLib {
3435 Root (Handle<DefaultProbabilityTermStructure> dts, Real pd)
3536 : dts_(std::move(dts)), pd_(pd) {}
3637 Real operator ()(Real t) const {
37- QL_REQUIRE (t >= 0.0 , " t < 0" );
38+ QL_REQUIRE (t >= 0.0 , " GaussianRandomDefaultModel: internal error, t < 0 ("
39+ << t << " ) during root searching." );
3840 return dts_->defaultProbability (t, true ) - pd_;
3941 }
4042 private:
@@ -72,9 +74,21 @@ namespace QuantLib {
7274 Real p = CumulativeNormalDistribution ()(y);
7375
7476 if (dts->defaultProbability (tmax) < p)
75- pool_->setTime (name, tmax+1 );
76- else
77- pool_->setTime (name, Brent ().solve (Root (dts,p),accuracy_,0 ,1 ));
77+ pool_->setTime (name, tmax + 1 );
78+ else {
79+ // we know there is a zero of f(t) = dts->defaultProbability(t) - p in [0, tmax]
80+ try {
81+ // try bracketing the root and find it with Brent
82+ Brent brent;
83+ brent.setLowerBound (0.0 );
84+ brent.setUpperBound (tmax);
85+ pool_->setTime (name, brent.solve (Root (dts, p), accuracy_, tmax / 2.0 , 1.0 ));
86+ } catch (...) {
87+ // if Brent fails, use Bisection, this is guaranteed to find the root
88+ pool_->setTime (
89+ name, Bisection ().solve (Root (dts, p), accuracy_, tmax / 2.0 , 0.0 , tmax));
90+ }
91+ }
7892 }
7993 }
8094
0 commit comments