C++ Brent::solve方法代码示例

    Real LinearTsrPricer::strikeFromVegaRatio(Real ratio,
                                              Option::Type optionType,
                                              Real referenceStrike) const {

        Real a, b, min, max, k;
        if (optionType == Option::Call) {
            a = swapRateValue_;
            min = referenceStrike;
            b = max = k =
                std::min(smileSection_->maxStrike(), adjustedUpperBound_);
        } else {
            a = min = k =
                std::max(smileSection_->minStrike(), adjustedLowerBound_);
            b = swapRateValue_;
            max = referenceStrike;
        }

        VegaRatioHelper h(&*smileSection_,
                          smileSection_->vega(swapRateValue_) * ratio);
        Brent solver;

        try {
            k = solver.solve(h, 1.0E-5, (a + b) / 2.0, a, b);
        }
        catch (...) {
            // use default value set above
        }

        return std::min(std::max(k, min), max);
    }

    Spread CallableBond::OAS(Real cleanPrice,
                             const Handle<YieldTermStructure>& engineTS,
                             const DayCounter& dayCounter,
                             Compounding compounding,
                             Frequency frequency,
                             Date settlement,
                             Real accuracy,
                             Size maxIterations,
                             Spread guess)
    {
        if (settlement == Date())
            settlement = settlementDate();

        Real dirtyPrice = cleanPrice + accruedAmount(settlement);

        boost::function<Real(Real)> f = NPVSpreadHelper(*this);
        OASHelper obj(f, dirtyPrice);

        Brent solver;
        solver.setMaxEvaluations(maxIterations);

        Real step = 0.001;
        Spread oas=solver.solve(obj, accuracy, guess, step);

        return continuousToConv(oas,
                                *this,
                                engineTS,
                                dayCounter,
                                compounding,
                                frequency);
    }

    Real LinearTsrPricer::strikeFromPrice(Real price, Option::Type optionType,
                                          Real referenceStrike) const {

        Real a, b, min, max, k;
        if (optionType == Option::Call) {
            a = swapRateValue_;
            min = referenceStrike;
            b = max = k =
                std::min(smileSection_->maxStrike(), settings_.upperRateBound_);
        } else {
            a = min = k =
                std::max(smileSection_->minStrike(), settings_.lowerRateBound_);
            b = swapRateValue_;
            max = referenceStrike;
        }

        PriceHelper h(&*smileSection_, optionType, price);
        Brent solver;

        try {
            k = solver.solve(h, 1.0E-5, swapRateValue_, a, b);
        }
        catch (...) {
            // use default value set above
        }

        return std::min(std::max(k, min), max);
    }

    Volatility CalibrationHelper::impliedVolatility(Real targetValue,
                                                    Real accuracy,
                                                    Size maxEvaluations,
                                                    Volatility minVol,
                                                    Volatility maxVol) const {

        ImpliedVolatilityHelper f(*this,targetValue);
        Brent solver;
        solver.setMaxEvaluations(maxEvaluations);
        return solver.solve(f,accuracy,volatility_->value(),minVol,maxVol);
    }

	double Gamma2Distribution::cumulativeInv(const double p, const bool fast, const long bins) {
			if(fast) {
				if(!fastCalculated_ || bins!=fastBins_) {
					fastBins_=bins;
					preCalculateCumulativeInv();
				}
				return fastCumulativeInv(p);
			}
			G2InvHelper hlp(this,p);
			Brent b;
			return b.solve(hlp,INVG2ACC,1.0,0.1);
	}

std::vector<Volatility> OptionletStripper2::spreadsVolImplied() const {

    Brent solver;
    std::vector<Volatility> result(nOptionExpiries_);
    Volatility guess = 0.0001, minSpread = -0.1, maxSpread = 0.1;
    for (Size j=0; j<nOptionExpiries_; ++j) {
        ObjectiveFunction f(stripper1_, caps_[j], atmCapFloorPrices_[j]);
        solver.setMaxEvaluations(maxEvaluations_);
        Volatility root = solver.solve(f, accuracy_, guess,
                                       minSpread, maxSpread);
        result[j] = root;
    }
    return result;
}

 Volatility CallableBond::impliedVolatility(
                           Real targetValue,
                           const Handle<YieldTermStructure>& discountCurve,
                           Real accuracy,
                           Size maxEvaluations,
                           Volatility minVol,
                           Volatility maxVol) const {
     calculate();
     QL_REQUIRE(!isExpired(), "instrument expired");
     Volatility guess = 0.5*(minVol + maxVol);
     blackDiscountCurve_.linkTo(*discountCurve, false);
     ImpliedVolHelper f(*this,targetValue);
     Brent solver;
     solver.setMaxEvaluations(maxEvaluations);
     return solver.solve(f, accuracy, guess, minVol, maxVol);
 }

    Volatility HestonBlackVolSurface::blackVolImpl(Time t, Real strike) const {
        const boost::shared_ptr<HestonProcess> process = hestonModel_->process();

        const DiscountFactor df = process->riskFreeRate()->discount(t, true);
        const DiscountFactor div = process->dividendYield()->discount(t, true);
        const Real spotPrice = process->s0()->value();

        const Real fwd = spotPrice
            * process->dividendYield()->discount(t, true)
            / process->riskFreeRate()->discount(t, true);


        const PlainVanillaPayoff payoff(
            fwd > strike ? Option::Put : Option::Call, strike);

        const Real kappa = hestonModel_->kappa();
        const Real theta = hestonModel_->theta();
        const Real rho   = hestonModel_->rho();
        const Real sigma = hestonModel_->sigma();
        const Real v0    = hestonModel_->v0();

        const AnalyticHestonEngine::ComplexLogFormula cpxLogFormula
            = AnalyticHestonEngine::Gatheral;

        const AnalyticHestonEngine* const hestonEnginePtr = 0;

        Real npv;
        Size evaluations;

        AnalyticHestonEngine::doCalculation(
            df, div, spotPrice, strike, t,
            kappa, theta, sigma, v0, rho,
            payoff, integration_, cpxLogFormula,
            hestonEnginePtr, npv, evaluations);

        if (npv <= 0.0) return std::sqrt(theta);

        Brent solver;
        solver.setMaxEvaluations(10000);
        const Volatility guess = std::sqrt(theta);
        const Real accuracy = std::numeric_limits<Real>::epsilon();

        const boost::function<Real(Real)> f = boost::bind(
            &blackValue, payoff.optionType(), strike, fwd, t, _1, df, npv);

        return solver.solve(f, accuracy, guess, 0.01);
    }

    Real InverseNonCentralChiSquareDistribution::operator()(Real x) const {

        // first find the right side of the interval
        Real upper = guess_;
        Size evaluations = maxEvaluations_;
        while (nonCentralDist_(upper) < x && evaluations > 0) {
            upper*=2.0;
            --evaluations;
        }

        // use a brent solver for the rest
        Brent solver;
        solver.setMaxEvaluations(evaluations);
        return solver.solve(compose(std::bind2nd(std::minus<Real>(),x), 
                                    nonCentralDist_),
                            accuracy_, 0.75*upper, 
                            (evaluations == maxEvaluations_)? 0.0: 0.5*upper,
                            upper);
    }

    Volatility CdsOption::impliedVolatility(
                   Real targetValue,
                   const Handle<YieldTermStructure>& termStructure,
                   const Handle<DefaultProbabilityTermStructure>& probability,
                   Real recoveryRate,
                   Real accuracy,
                   Size maxEvaluations,
                   Volatility minVol,
                   Volatility maxVol) const {
        calculate();
        QL_REQUIRE(!isExpired(), "instrument expired");

        Volatility guess = 0.10;

        ImpliedVolHelper f(*this, probability, recoveryRate,
                           termStructure, targetValue);
        Brent solver;
        solver.setMaxEvaluations(maxEvaluations);
        return solver.solve(f, accuracy, guess, minVol, maxVol);
    }

        Volatility ImpliedVolatilityHelper::calculate(
                                                 const Instrument& instrument,
                                                 const PricingEngine& engine,
                                                 SimpleQuote& volQuote,
                                                 Real targetValue,
                                                 Real accuracy,
                                                 Natural maxEvaluations,
                                                 Volatility minVol,
                                                 Volatility maxVol) {

            instrument.setupArguments(engine.getArguments());
            engine.getArguments()->validate();

            PriceError f(engine, volQuote, targetValue);
            Brent solver;
            solver.setMaxEvaluations(maxEvaluations);
            Volatility guess = (minVol+maxVol)/2.0;
            Volatility result = solver.solve(f, accuracy, guess,
                                             minVol, maxVol);
            return result;
        }

	Real CmsSpreadOption::impliedCorrelation(boost::shared_ptr<CmsPricer> pricer,const Period& pillar, const double price,Size preCalculatedFixings,double preCalculatedPrice) {
		Brent b;
		double res=b.solve(CmsSpreadOptionImplCorrHelper(this,pricer,pillar,price,preCalculatedFixings,preCalculatedPrice),CORRACC,pricer->correlationTermStructure()->correlation(pillar,true),CORRSTEP);
		double res2=atan(res)*2.0/M_PI;
		pricer->correlationTermStructure()->setPillarCorrelation(pillar,res2);
		return res2;
		/*CmsSpreadOptionImplCorrHelper h(boost::shared_ptr<CmsSpreadOption>(this),pricer,pillar,price);
		double x1=-MAXCORR,x2=MAXCORR,x;
		double f1=h(x1);
		double f2=h(x2);
		double f;
		double step=MAXCORR+MAXCORR;
		if(f1*f2>=0.0) QL_FAIL("No sign change for correlations -1,1");
		do {
			x=(x1+x2)/2.0;
			f=h(x);
			if(f1*f<0.0) x2=x; else x1=x;
			step/=step;
			pricer->correlationTermStructure()->setPillarCorrelation(pillar,x);
		} while(step>CORRACC);
		pricer->correlationTermStructure()->setPillarCorrelation(pillar,x);
		return x;*/
	}

    void IterativeBootstrap<Curve>::calculate() const {

        Size n = ts_->instruments_.size();

        // ensure rate helpers are sorted
        std::sort(ts_->instruments_.begin(), ts_->instruments_.end(),
                  detail::BootstrapHelperSorter());

        // check that there is no instruments with the same maturity
        for (Size i=1; i<n; ++i) {
            Date m1 = ts_->instruments_[i-1]->latestDate(),
                 m2 = ts_->instruments_[i]->latestDate();
            QL_REQUIRE(m1 != m2,
                       "two instruments have the same maturity ("<< m1 <<")");
        }

        // check that there is no instruments with invalid quote
        for (Size i=0; i<n; ++i)
            QL_REQUIRE(ts_->instruments_[i]->quote()->isValid(),
                       io::ordinal(i+1) << " instrument (maturity: " <<
                       ts_->instruments_[i]->latestDate() <<
                       ") has an invalid quote");

        // setup instruments
        for (Size i=0; i<n; ++i) {
            // don't try this at home!
            // This call creates instruments, and removes "const".
            // There is a significant interaction with observability.
            ts_->instruments_[i]->setTermStructure(const_cast<Curve*>(ts_));
        }

        // calculate dates and times
        ts_->dates_ = std::vector<Date>(n+1);
        ts_->times_ = std::vector<Time>(n+1);
        ts_->dates_[0] = Traits::initialDate(ts_);
        ts_->times_[0] = ts_->timeFromReference(ts_->dates_[0]);
        for (Size i=0; i<n; ++i) {
            ts_->dates_[i+1] = ts_->instruments_[i]->latestDate();
            ts_->times_[i+1] = ts_->timeFromReference(ts_->dates_[i+1]);
        }

        // set initial guess only if the current curve cannot be used as guess
        if (validCurve_) {
            QL_ENSURE(ts_->data_.size() == n+1,
                      "dimension mismatch: expected " << n+1 <<
                      ", actual " << ts_->data_.size());
        } else {
            ts_->data_ = std::vector<Rate>(n+1);
            ts_->data_[0] = Traits::initialValue(ts_);
            for (Size i=0; i<n; ++i)
                ts_->data_[i+1] = Traits::initialGuess();
        }

        Brent solver;
        Size maxIterations = Traits::maxIterations();

        for (Size iteration=0; ; ++iteration) {
            std::vector<Rate> previousData = ts_->data_;
            // restart from the previous interpolation
            if (validCurve_) {
                ts_->interpolation_ = ts_->interpolator_.interpolate(
                                                      ts_->times_.begin(),
                                                      ts_->times_.end(),
                                                      ts_->data_.begin());
            }
            for (Size i=1; i<n+1; ++i) {

                // calculate guess before extending interpolation
                // to ensure that any extrapolation is performed
                // using the curve bootstrapped so far and no more
                boost::shared_ptr<typename Traits::helper> instrument =
                    ts_->instruments_[i-1];
                Rate guess = 0.0;
                if (validCurve_ || iteration>0) {
                    guess = ts_->data_[i];
                } else if (i==1) {
                    guess = Traits::initialGuess();
                } else {
                    // most traits extrapolate
                    guess = Traits::guess(ts_, ts_->dates_[i]);
                }

                // bracket
                Real min = Traits::minValueAfter(i, ts_->data_);
                Real max = Traits::maxValueAfter(i, ts_->data_);
                if (guess<=min || guess>=max)
                    guess = (min+max)/2.0;

                if (!validCurve_ && iteration == 0) {
                    // extend interpolation a point at a time
                    try {
                        ts_->interpolation_ = ts_->interpolator_.interpolate(
                                                      ts_->times_.begin(),
                                                      ts_->times_.begin()+i+1,
                                                      ts_->data_.begin());
                    } catch(...) {
                        if (!Interpolator::global)
                            throw; // no chance to fix it in a later iteration

                        // otherwise, if the target interpolation is
//.........这里部分代码省略.........

    void InterpolatedYoYOptionletStripper<Interpolator1D>::
    initialize(const boost::shared_ptr<YoYCapFloorTermPriceSurface> &s,
               const boost::shared_ptr<YoYInflationCapFloorEngine> &p,
               const Real slope) const {
        YoYCapFloorTermPriceSurface_ = s;
        p_ = p;
        lag_ = YoYCapFloorTermPriceSurface_->observationLag();
        frequency_ = YoYCapFloorTermPriceSurface_->frequency();
        indexIsInterpolated_ = YoYCapFloorTermPriceSurface_->indexIsInterpolated();
        Natural fixingDays_ = YoYCapFloorTermPriceSurface_->fixingDays();
        Natural settlementDays = 0; // always
        Calendar cal = YoYCapFloorTermPriceSurface_->calendar();
        BusinessDayConvention bdc =
            YoYCapFloorTermPriceSurface_->businessDayConvention();
        DayCounter dc = YoYCapFloorTermPriceSurface_->dayCounter();

        // switch from caps to floors when out of floors
        Rate maxFloor = YoYCapFloorTermPriceSurface_->floorStrikes().back();
        YoYInflationCapFloor::Type useType = YoYInflationCapFloor::Floor;
        Period TPmin = YoYCapFloorTermPriceSurface_->maturities().front();
        // create a "fake index" based on Generic, this should work
        // provided that the lag and frequency are correct
        RelinkableHandle<YoYInflationTermStructure> hYoY(
                                       YoYCapFloorTermPriceSurface_->YoYTS());
        boost::shared_ptr<YoYInflationIndex> anIndex(
                                           new YYGenericCPI(frequency_, false,
                                                            false, lag_,
                                                            Currency(), hYoY));

        // strip each K separatly
        for (Size i=0; i<YoYCapFloorTermPriceSurface_->strikes().size(); i++) {
            Rate K = YoYCapFloorTermPriceSurface_->strikes()[i];
            if (K > maxFloor) useType = YoYInflationCapFloor::Cap;

            // solve for the initial point on the vol curve
            Brent solver;
            Real solverTolerance_ = 1e-7;
             // these are VOLATILITY guesses (always +)
            Real lo = 0.00001, hi = 0.08;
            Real guess = (hi+lo)/2.0;
            Real found;
            Real priceToMatch =
                (useType == YoYInflationCapFloor::Cap ?
                 YoYCapFloorTermPriceSurface_->capPrice(TPmin, K) :
                 YoYCapFloorTermPriceSurface_->floorPrice(TPmin, K));

            try{
                found = solver.solve(
                      ObjectiveFunction(useType, slope, K, lag_, fixingDays_,
                                        anIndex, YoYCapFloorTermPriceSurface_,
                                        p_, priceToMatch),
                      solverTolerance_, guess, lo, hi );
            } catch( std::exception &e) {
                QL_FAIL("failed to find solution here because: " << e.what());
            }

            // ***create helpers***
            Real notional = 10000; // work in bps
            std::vector<boost::shared_ptr<BootstrapHelper<YoYOptionletVolatilitySurface> > > helperInstruments;
            std::vector<boost::shared_ptr<YoYOptionletHelper> > helpers;
            for (Size j = 0; j < YoYCapFloorTermPriceSurface_->maturities().size(); j++){
                Period Tp = YoYCapFloorTermPriceSurface_->maturities()[j];

                Real nextPrice =
                    (useType == YoYInflationCapFloor::Cap ?
                     YoYCapFloorTermPriceSurface_->capPrice(Tp, K) :
                     YoYCapFloorTermPriceSurface_->floorPrice(Tp, K));

                Handle<Quote> quote1(boost::shared_ptr<Quote>(
                                               new SimpleQuote( nextPrice )));
                // helper should be an integer number of periods away,
                // this is enforced by rounding
                Size nT = (Size)floor(s->timeFromReference(s->yoyOptionDateFromTenor(Tp))+0.5);
                helpers.push_back(boost::shared_ptr<YoYOptionletHelper>(
                          new YoYOptionletHelper(quote1, notional, useType,
                                                 lag_,
                                                 dc, cal,
                                                 fixingDays_,
                                                 anIndex, K, nT, p_)));

                boost::shared_ptr<ConstantYoYOptionletVolatility> yoyVolBLACK(
                          new ConstantYoYOptionletVolatility(found, settlementDays,
                                                             cal, bdc, dc,
                                                             lag_, frequency_,
                                                             false,
                                                             // -100% to +300%
                                                             -1.0,3.0));

                helpers[j]->setTermStructure(
                       // gets underlying pointer & removes const
                       const_cast<ConstantYoYOptionletVolatility*>(
                                                          yoyVolBLACK.get()));
                helperInstruments.push_back(helpers[j]);
            }
            // ***bootstrap***
            // this is the artificial vol at zero so that first section works
            Real Tmin = s->timeFromReference(s->yoyOptionDateFromTenor(TPmin));
            Volatility baseYoYVolatility = found - slope * Tmin * found;
            Rate eps = std::max(K, 0.02) / 1000.0;
            Rate minStrike = K-eps;
//.........这里部分代码省略.........

NoArbSabrModel::NoArbSabrModel(const Real expiryTime, const Real forward,
                               const Real alpha, const Real beta, const Real nu,
                               const Real rho)
    : expiryTime_(expiryTime), externalForward_(forward), alpha_(alpha),
      beta_(beta), nu_(nu), rho_(rho), forward_(forward),
      numericalForward_(forward) {

    using namespace ext::placeholders;

    QL_REQUIRE(expiryTime > 0.0 && expiryTime <= detail::NoArbSabrModel::expiryTime_max,
               "expiryTime (" << expiryTime << ") out of bounds");
    QL_REQUIRE(forward > 0.0, "forward (" << forward << ") must be positive");
    QL_REQUIRE(beta >= detail::NoArbSabrModel::beta_min && beta <= detail::NoArbSabrModel::beta_max,
               "beta (" << beta << ") out of bounds");
    Real sigmaI = alpha * std::pow(forward, beta - 1.0);
    QL_REQUIRE(sigmaI >= detail::NoArbSabrModel::sigmaI_min &&
                   sigmaI <= detail::NoArbSabrModel::sigmaI_max,
               "sigmaI = alpha*forward^(beta-1.0) ("
                   << sigmaI << ") out of bounds, alpha=" << alpha
                   << " beta=" << beta << " forward=" << forward);
    QL_REQUIRE(nu >= detail::NoArbSabrModel::nu_min && nu <= detail::NoArbSabrModel::nu_max,
               "nu (" << nu << ") out of bounds");
    QL_REQUIRE(rho >= detail::NoArbSabrModel::rho_min && rho <= detail::NoArbSabrModel::rho_max,
               "rho (" << rho << ") out of bounds");

    // determine a region sufficient for integration in the normal case

    fmin_ = fmax_ = forward_;
    for (Real tmp = p(fmax_);
         tmp > std::max(detail::NoArbSabrModel::i_accuracy / std::max(1.0, fmax_ - fmin_),
                        detail::NoArbSabrModel::density_threshold);
         tmp = p(fmax_)) {
        fmax_ *= 2.0;
    }
    for (Real tmp = p(fmin_);
         tmp > std::max(detail::NoArbSabrModel::i_accuracy / std::max(1.0, fmax_ - fmin_),
                        detail::NoArbSabrModel::density_threshold);
         tmp = p(fmin_)) {
        fmin_ *= 0.5;
    }
    fmin_ = std::max(detail::NoArbSabrModel::strike_min, fmin_);

    QL_REQUIRE(fmax_ > fmin_, "could not find a reasonable integration domain");

    integrator_ =
        ext::make_shared<GaussLobattoIntegral>(
            detail::NoArbSabrModel::i_max_iterations, detail::NoArbSabrModel::i_accuracy);

    detail::D0Interpolator d0(forward_, expiryTime_, alpha_, beta_, nu_, rho_);
    absProb_ = d0();

    try {
        Brent b;
        Real start = std::sqrt(externalForward_ - detail::NoArbSabrModel::strike_min);
        Real tmp =
            b.solve(ext::bind(&NoArbSabrModel::forwardError, this, _1),
                    detail::NoArbSabrModel::forward_accuracy, start,
                    std::min(detail::NoArbSabrModel::forward_search_step, start / 2.0));
        forward_ = tmp * tmp + detail::NoArbSabrModel::strike_min;
    } catch (Error&) {
        // fall back to unadjusted forward
        forward_ = externalForward_;
    }

    Real d = forwardError(std::sqrt(forward_ - detail::NoArbSabrModel::strike_min));
    numericalForward_ = d + externalForward_;
}