C++ QL_REQUIRE函数代码示例

    std::vector<string> qlRateHelperSelection(
        const std::vector<shared_ptr<QuantLibAddin::RateHelper> >& qlarhs,
        const std::vector<QuantLib::Natural>& priority,
        QuantLib::Natural nImmFutures,
        QuantLib::Natural nSerialFutures,
        QuantLib::Natural frontFuturesRollingDays,
        RateHelper::DepoInclusionCriteria depoInclusionCriteria,
        const std::vector<QuantLib::Natural>& minDist)
    {
        // Checks
        QL_REQUIRE(!qlarhs.empty(), "no instrument given");
        QuantLib::Size nInstruments = qlarhs.size();
        QL_REQUIRE(priority.size()==nInstruments,
                   "priority (" << priority.size() <<
                   ") / instruments (" << nInstruments << ") size mismatch");
        QL_REQUIRE(minDist.size()==nInstruments || minDist.size()==1,
                   "minDist (" << minDist.size() <<
                   ") / instruments (" << nInstruments << ") mismatch");

        // RateHelperItem
        shared_ptr<QuantLibAddin::RateHelper> qlarh;
        shared_ptr<QuantLib::RateHelper> qlrh;
        std::vector<RateHelperItem> rhsAll;
        rhsAll.reserve(nInstruments);
        for (QuantLib::Size i=0; i<nInstruments; ++i) {
            qlarh = qlarhs[i];
            qlarh->getLibraryObject(qlrh);
            string qlarh_id = convert2<string>(
                qlarh->propertyValue("OBJECTID"));
            bool isFutures = dynamic_pointer_cast<FuturesRateHelper>(qlarh);
            bool isImmFutures = false, isSerialFutures = false;
            if (isFutures) {
                isImmFutures = QuantLib::IMM::isIMMdate(qlrh->earliestDate());
                isSerialFutures = !isImmFutures;
            }
            bool isDepo = dynamic_pointer_cast<DepositRateHelper>(qlarh);
            rhsAll.push_back(RateHelperItem(isImmFutures,
                                            isSerialFutures,
                                            isDepo,
                                            qlarh_id,
                                            priority[i],
                                            qlrh->earliestDate(),
                                            qlrh->latestDate(),
                                            minDist.size()==1 ? minDist[0] : minDist[i]));
        }

        // Preliminary sort of RateHelperItems according to
        // their latest date and priority
        std::sort(rhsAll.begin(), rhsAll.end(), RateHelperPrioritySorter());

        // Select input rate helpers according to:
        // expiration, maximum number of allowed Imm and Serial Futures, Depo/Futures priorities
        QuantLib::Natural immFuturesCounter = 0;
        QuantLib::Natural serialFuturesCounter = 0;
        QuantLib::Date evalDate = QuantLib::Settings::instance().evaluationDate();
        std::vector<RateHelperItem> rhs, rhsDepo;

        // Look for the front Futures, if any
        bool thereAreFutures = false;
        QuantLib::Date frontFuturesEarliestDate, frontFuturesLatestDate;
        if (nImmFutures>0 || nSerialFutures>0) {
            QuantLib::Size j=0;
            while (j<nInstruments) {
                if (nImmFutures>0 && rhsAll[j].isImmFutures &&
                        (rhsAll[j].earliestDate-frontFuturesRollingDays >= evalDate)) {
                    thereAreFutures = true;
                    frontFuturesEarliestDate = rhsAll[j].earliestDate;
                    frontFuturesLatestDate = rhsAll[j].latestDate;
                    break;
                }
                if (nSerialFutures>0 && rhsAll[j].isSerialFutures &&
                        (rhsAll[j].earliestDate-frontFuturesRollingDays >= evalDate)) {
                    thereAreFutures = true;
                    frontFuturesEarliestDate = rhsAll[j].earliestDate;
                    frontFuturesLatestDate = rhsAll[j].latestDate;
                    break;
                }
                ++j;
            }
        }

        // If there are NOT Futures, include all Depos
        if (!thereAreFutures)
            depoInclusionCriteria = RateHelper::AllDepos;

        // Start selection
        bool depoAfterFrontFuturesAlreadyIncluded = false;
        for (QuantLib::Size i=0; i<nInstruments; ++i) {
            if (rhsAll[i].earliestDate >= evalDate) {
                if (rhsAll[i].isDepo) {                 // Check Depo conditions
                    switch (depoInclusionCriteria) {
                        case RateHelper::AllDepos:
                       // Include all depos
                            rhs.push_back(rhsAll[i]);
                            break;
                        case RateHelper::DeposBeforeFirstFuturesStartDate:
                        // Include only depos with maturity date before
                        // the front Futures start date
                            if (rhsAll[i].latestDate < frontFuturesEarliestDate)
                                rhs.push_back(rhsAll[i]);
//.........这里部分代码省略.........

	Real CreditDefaultSwap::accrualRebateNPV() const {
		calculate();
		QL_REQUIRE(accrualRebateNPV_ != Null<Real>(),
			       "accrual Rebate NPV not available");
		return accrualRebateNPV_;
	}

    CreditDefaultSwap::CreditDefaultSwap(Protection::Side side,
                                         Real notional,
                                         Rate upfront,
                                         Rate runningSpread,
                                         const Schedule& schedule,
                                         BusinessDayConvention convention,
                                         const DayCounter& dayCounter,
                                         bool settlesAccrual,
                                         bool paysAtDefaultTime,
                                         const Date& protectionStart,
                                         const Date& upfrontDate,
                                         const boost::shared_ptr<Claim>& claim,
										 const DayCounter& lastPeriodDayCounter,
										 const Natural standardCdsStartDelayDays,
										 const Natural standardCdsUpfrontDelayDays)
    : side_(side), notional_(notional), upfront_(upfront),
      runningSpread_(runningSpread), settlesAccrual_(settlesAccrual),
      paysAtDefaultTime_(paysAtDefaultTime), claim_(claim),
      protectionStart_(protectionStart == Null<Date>() ? schedule[0] :
                                                         protectionStart) {

        Date d = upfrontDate == Null<Date>() ? schedule[0] : upfrontDate;

        try {
            if (schedule.rule() ==
                DateGeneration::CDS) { // a standard CDS is identified
                                       // by the schedule rule
                Date evalDate = Settings::instance().evaluationDate();
                if (protectionStart == Null<Date>())
                    protectionStart_ =
                        evalDate + standardCdsStartDelayDays; // if protection
                // start is given it
                // is not set here
                if (upfrontDate == Null<Date>())
                    d = schedule.calendar().advance(
                        evalDate, standardCdsUpfrontDelayDays * Days,
                        schedule.businessDayConvention(),
                        schedule.endOfMonth()); // if upfront date is
                                                // given it is not set
                                                // here
                QL_REQUIRE(protectionStart_ >= schedule[0],
                           "protection must start on or after accrual "
                           "start for standard CDS");
            }
        } catch (...) {
            QL_REQUIRE(protectionStart_ <= schedule[0],
                       "protection can not start after accrual for non "
                       "standard CDS");
        }

        leg_ = FixedRateLeg(schedule)
                   .withNotionals(notional)
                   .withCouponRates(runningSpread, dayCounter)
                   .withPaymentAdjustment(convention)
                   .withLastPeriodDayCounter(lastPeriodDayCounter);

        upfrontPayment_.reset(new SimpleCashFlow(notional*upfront, d));

        QL_REQUIRE(upfrontPayment_->date() >= protectionStart_,
                   "upfront can not be due before contract start");

        if (!claim_)
            claim_ = boost::shared_ptr<Claim>(new FaceValueClaim);
        registerWith(claim_);
    }

 inline Real GeneralStatistics::max() const {
     QL_REQUIRE(samples() > 0, "empty sample set");
     return std::max_element(samples_.begin(),
                             samples_.end())->first;
 }

    void KahaleSmileSection::compute() {

        std::pair<Size,Size> afIdx = ssutils_->arbitragefreeIndices();
        leftIndex_ = afIdx.first;
        rightIndex_ = afIdx.second;

        if(deleteArbitragePoints_) {
            while(leftIndex_>1 || rightIndex_<k_.size()-1) {
                
                ssutils_ = boost::shared_ptr<SmileSectionUtils>(new SmileSectionUtils(*source_,moneynessGrid_,f_));
                std::pair<Size,Size> afIdx = ssutils_->arbitragefreeIndices();

                leftIndex_ = afIdx.first;
                rightIndex_ = afIdx.second;

                QL_REQUIRE(rightIndex_>leftIndex_,
                           "arbitrage free region must at least contain two points (only index is " << leftIndex_ << ")");
            
                if(leftIndex_>1) {
                    moneynessGrid_.erase(moneynessGrid_.begin()+leftIndex_-1);
                    k_.erase(k_.begin()+leftIndex_-1);
                    c_.erase(c_.begin()+leftIndex_-1);
                    leftIndex_--;
                    rightIndex_--; 
                }

                if(rightIndex_<k_.size()-1) {
                    moneynessGrid_.erase(moneynessGrid_.begin()+rightIndex_+1);
                    k_.erase(k_.begin()+rightIndex_+1);
                    c_.erase(c_.begin()+rightIndex_+1);
                    rightIndex_--;
                }
            }
        }

        cFunctions_ = std::vector<boost::shared_ptr<cFunction> >(rightIndex_-leftIndex_+2);

        // extrapolation in the leftmost interval

        Brent brent;
        bool success;
        Real secl = 0.0;

        do {
            success=true;
            try {
                Real k1 = k_[leftIndex_];
                Real c1 = c_[leftIndex_];
                Real c0 = c_[0];
                secl = (c_[leftIndex_]-c_[0]) / (k_[leftIndex_]-k_[0]);
                Real sec = (c_[leftIndex_+1]-c_[leftIndex_]) / (k_[leftIndex_+1]-k_[leftIndex_]);
                Real c1p;
                if(interpolate_) c1p=(secl+sec)/2;
                else {
                    c1p=(blackFormula(Option::Call, k1+gap_, f_, sqrt(source_->variance(k1+gap_)))-
                         blackFormula(Option::Call, k1, f_, sqrt(source_->variance(k1))))/gap_;
                    QL_REQUIRE(secl < c1p && c1p <= 0.0,"dummy"); 
                    // can not extrapolate so throw exception which is caught below
                }
                sHelper1 sh1(k1,c0,c1,c1p);
                Real s = brent.solve(sh1,QL_KAHALE_ACC,0.20,0.00,QL_KAHALE_SMAX); // numerical parameters hardcoded here
                sh1(s);
                boost::shared_ptr<cFunction> cFct1(new cFunction(sh1.f_,s,0.0,sh1.b_));
                cFunctions_[0]=cFct1;
            } catch(...) {
                leftIndex_++;
                success=false;
            }
        } while(!success && leftIndex_ < rightIndex_);

        QL_REQUIRE(leftIndex_ < rightIndex_, "can not extrapolate to left, right index of af region reached (" 
                   << rightIndex_ << ")");

        // interpolation

        Real cp0 = 0.0, cp1 = 0.0;

        if(interpolate_) {

            for(Size i = leftIndex_; i<rightIndex_; i++) {
                Real k0 = k_[i];
                Real k1 = k_[i+1];
                Real c0 = c_[i];
                Real c1 = c_[i+1];
                Real sec = (c_[i+1]-c_[i]) / (k_[i+1]-k_[i]);
                if(i==leftIndex_) cp0 = leftIndex_ > 0 ? (secl + sec) / 2.0 : sec;
                Real secr;
                if(i==rightIndex_-1) secr=0.0;
                else secr = (c_[i+2]-c_[i+1]) / (k_[i+2]-k_[i+1]);
                cp1 = (sec+secr) / 2.0;
                aHelper ah(k0,k1,c0,c1,cp0,cp1);
                Real a;
                try {
                    a = brent.solve(ah,QL_KAHALE_ACC,0.5*(cp1+(1.0+cp0)),cp1+QL_KAHALE_EPS,1.0+cp0-QL_KAHALE_EPS);
                    // numerical parameters hardcoded here
                } catch(...) {
                    // from theory there must exist a zero. if the solver does not find it, it most
                    // probably lies close one of the interval bounds. Just choose the better bound
                    // and relax the accuracy. This does not matter in practice usually.
                    Real la = std::fabs(ah(cp1+QL_KAHALE_EPS));
//.........这里部分代码省略.........

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

        validCurve_ = false;
        Size nInsts = 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<nInsts; ++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<nInsts; ++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<nInsts; ++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_));
        }
        // set initial guess only if the current curve cannot be used as guess
        if (validCurve_)
            QL_ENSURE(ts_->data_.size() == nInsts+1,
                      "dimension mismatch: expected " << nInsts+1 <<
                      ", actual " << ts_->data_.size());
        else {
            ts_->data_ = std::vector<Rate>(nInsts+1);
            ts_->data_[0] = Traits::initialValue(ts_);
        }

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

        LevenbergMarquardt solver(ts_->accuracy_,
                                  ts_->accuracy_,
                                  ts_->accuracy_);
        EndCriteria endCriteria(100, 10, 0.00, ts_->accuracy_, 0.00);
        PositiveConstraint posConstraint;
        NoConstraint noConstraint;
        Constraint& solverConstraint = forcePositive_ ?
            static_cast<Constraint&>(posConstraint) :
            static_cast<Constraint&>(noConstraint);

        // now start the bootstrapping.
        Size iInst = localisation_-1;

        Size dataAdjust = Curve::interpolator_type::dataSizeAdjustment;

        do {
            Size initialDataPt = iInst+1-localisation_+dataAdjust;
            Array startArray(localisation_+1-dataAdjust);
            for (Size j = 0; j < startArray.size()-1; ++j)
                startArray[j] = ts_->data_[initialDataPt+j];

            // here we are extending the interpolation a point at a
            // time... but the local interpolator can make an
            // approximation for the final localisation period.
            // e.g. if the localisation is 2, then the first section
            // of the curve will be solved using the first 2
            // instruments... with the local interpolator making
            // suitable boundary conditions.
            ts_->interpolation_ =
                ts_->interpolator_.localInterpolate(
                                              ts_->times_.begin(),
                                              ts_->times_.begin()+(iInst + 2),
                                              ts_->data_.begin(),
                                              localisation_,
                                              ts_->interpolation_,
                                              nInsts+1);

            if (iInst >= localisation_) {
                startArray[localisation_-dataAdjust] =
                    Traits::guess(iInst, ts_, false, 0); // ?
            } else {
                startArray[localisation_-dataAdjust] = ts_->data_[0];
            }

            PenaltyFunction<Curve> currentCost(
                        ts_,
                        initialDataPt,
//.........这里部分代码省略.........

 inline std::string DayCounter::name() const {
     QL_REQUIRE(impl_, "no implementation provided");
     return impl_->name();
 }

    void FDStepConditionEngine<Scheme>::calculate(
                                            PricingEngine::results* r) const {
        OneAssetOption::results * results =
            dynamic_cast<OneAssetOption::results *>(r);
        setGridLimits();
        initializeInitialCondition();
        initializeOperator();
        initializeBoundaryConditions();
        initializeStepCondition();

        typedef FiniteDifferenceModel<ParallelEvolver<
                    Scheme<TridiagonalOperator> > > model_type;

        typename model_type::operator_type operatorSet;
        typename model_type::array_type arraySet;
        typename model_type::bc_set bcSet;
        typename model_type::condition_type conditionSet;

        prices_ = intrinsicValues_;

        controlPrices_ = intrinsicValues_;
        controlOperator_ = finiteDifferenceOperator_;
        controlBCs_[0] = BCs_[0];
        controlBCs_[1] = BCs_[1];

        operatorSet.push_back(finiteDifferenceOperator_);
        operatorSet.push_back(controlOperator_);

        arraySet.push_back(prices_.values());
        arraySet.push_back(controlPrices_.values());

        bcSet.push_back(BCs_);
        bcSet.push_back(controlBCs_);

        conditionSet.push_back(stepCondition_);
        conditionSet.push_back(boost::shared_ptr<StandardStepCondition>(
                                                   new NullCondition<Array>));

        model_type model(operatorSet, bcSet);

        model.rollback(arraySet, getResidualTime(),
                       0.0, timeSteps_, conditionSet);

        prices_.values() = arraySet[0];
        controlPrices_.values() = arraySet[1];

        boost::shared_ptr<StrikedTypePayoff> striked_payoff =
            boost::dynamic_pointer_cast<StrikedTypePayoff>(payoff_);
        QL_REQUIRE(striked_payoff, "non-striked payoff given");

        Real variance =
            process_->blackVolatility()->blackVariance(
                                     exerciseDate_, striked_payoff->strike());
        DiscountFactor dividendDiscount =
            process_->dividendYield()->discount(exerciseDate_);
        DiscountFactor riskFreeDiscount =
            process_->riskFreeRate()->discount(exerciseDate_);
        Real spot = process_->stateVariable()->value();
        Real forwardPrice = spot * dividendDiscount / riskFreeDiscount;

        BlackCalculator black(striked_payoff, forwardPrice,
                              std::sqrt(variance), riskFreeDiscount);

        results->value = prices_.valueAtCenter()
            - controlPrices_.valueAtCenter()
            + black.value();
        results->delta = prices_.firstDerivativeAtCenter()
            - controlPrices_.firstDerivativeAtCenter()
            + black.delta(spot);
        results->gamma = prices_.secondDerivativeAtCenter()
            - controlPrices_.secondDerivativeAtCenter()
            + black.gamma(spot);
        results->additionalResults["priceCurve"] = prices_;
    }

    void BinomialBarrierEngine<T,D>::calculate() const {

        DayCounter rfdc  = process_->riskFreeRate()->dayCounter();
        DayCounter divdc = process_->dividendYield()->dayCounter();
        DayCounter voldc = process_->blackVolatility()->dayCounter();
        Calendar volcal = process_->blackVolatility()->calendar();

        Real s0 = process_->stateVariable()->value();
        QL_REQUIRE(s0 > 0.0, "negative or null underlying given");
        Volatility v = process_->blackVolatility()->blackVol(
            arguments_.exercise->lastDate(), s0);
        Date maturityDate = arguments_.exercise->lastDate();
        Rate r = process_->riskFreeRate()->zeroRate(maturityDate,
            rfdc, Continuous, NoFrequency);
        Rate q = process_->dividendYield()->zeroRate(maturityDate,
            divdc, Continuous, NoFrequency);
        Date referenceDate = process_->riskFreeRate()->referenceDate();

        // binomial trees with constant coefficient
        Handle<YieldTermStructure> flatRiskFree(
            boost::shared_ptr<YieldTermStructure>(
                new FlatForward(referenceDate, r, rfdc)));
        Handle<YieldTermStructure> flatDividends(
            boost::shared_ptr<YieldTermStructure>(
                new FlatForward(referenceDate, q, divdc)));
        Handle<BlackVolTermStructure> flatVol(
            boost::shared_ptr<BlackVolTermStructure>(
                new BlackConstantVol(referenceDate, volcal, v, voldc)));

        boost::shared_ptr<StrikedTypePayoff> payoff =
            boost::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-striked payoff given");

        Time maturity = rfdc.yearFraction(referenceDate, maturityDate);

        boost::shared_ptr<StochasticProcess1D> bs(
                         new GeneralizedBlackScholesProcess(
                                      process_->stateVariable(),
                                      flatDividends, flatRiskFree, flatVol));

        // correct timesteps to ensure a (local) minimum, using Boyle and Lau
        // approach. See Journal of Derivatives, 1/1994,
        // "Bumping up against the barrier with the binomial method"
        // Note: this approach works only for CoxRossRubinstein lattices, so
        // is disabled if T is not a CoxRossRubinstein or derived from it.
        Size optimum_steps = timeSteps_;
        if (boost::is_base_of<CoxRossRubinstein, T>::value && 
            maxTimeSteps_ > timeSteps_ && s0 > 0 && arguments_.barrier > 0) {
            Real divisor;
            if (s0 > arguments_.barrier)
               divisor = std::pow(std::log(s0 / arguments_.barrier), 2);
            else
               divisor = std::pow(std::log(arguments_.barrier / s0), 2);
            if (!close(divisor,0)) {
                for (Size i=1; i < timeSteps_ ; ++i) {
                    Size optimum = Size(( i*i * v*v * maturity) / divisor);
                    if (timeSteps_ < optimum) {
                        optimum_steps = optimum;
                        break; // found first minimum with iterations>=timesteps
                    }
                }
            }

            if (optimum_steps > maxTimeSteps_) 
               optimum_steps = maxTimeSteps_; // too high, limit
        }

        TimeGrid grid(maturity, optimum_steps);

        boost::shared_ptr<T> tree(new T(bs, maturity, optimum_steps,
                                        payoff->strike()));

        boost::shared_ptr<BlackScholesLattice<T> > lattice(
            new BlackScholesLattice<T>(tree, r, maturity, optimum_steps));

        D option(arguments_, *process_, grid);
        option.initialize(lattice, maturity);

        // Partial derivatives calculated from various points in the
        // binomial tree 
        // (see J.C.Hull, "Options, Futures and other derivatives", 6th edition, pp 397/398)

        // Rollback to third-last step, and get underlying prices (s2) &
        // option values (p2) at this point
        option.rollback(grid[2]);
        Array va2(option.values());
        QL_ENSURE(va2.size() == 3, "Expect 3 nodes in grid at second step");
        Real p2u = va2[2]; // up
        Real p2m = va2[1]; // mid
        Real p2d = va2[0]; // down (low)
        Real s2u = lattice->underlying(2, 2); // up price
        Real s2m = lattice->underlying(2, 1); // middle price
        Real s2d = lattice->underlying(2, 0); // down (low) price

        // calculate gamma by taking the first derivate of the two deltas
        Real delta2u = (p2u - p2m)/(s2u-s2m);
        Real delta2d = (p2m-p2d)/(s2m-s2d);
        Real gamma = (delta2u - delta2d) / ((s2u-s2d)/2);

        // Rollback to second-last step, and get option values (p1) at
//.........这里部分代码省略.........

 void CatBond::arguments::validate() const {
     Bond::arguments::validate();
     QL_REQUIRE(notionalRisk, "null notionalRisk");
 }

    void LinearTsrPricer::initialize(const FloatingRateCoupon &coupon) {

        coupon_ = dynamic_cast<const CmsCoupon *>(&coupon);
        QL_REQUIRE(coupon_, "CMS coupon needed");
        gearing_ = coupon_->gearing();
        spread_ = coupon_->spread();

        fixingDate_ = coupon_->fixingDate();
        paymentDate_ = coupon_->date();
        swapIndex_ = coupon_->swapIndex();

        forwardCurve_ = swapIndex_->forwardingTermStructure();
        if (swapIndex_->exogenousDiscount())
            discountCurve_ = swapIndex_->discountingTermStructure();
        else
            discountCurve_ = forwardCurve_;

        // if no coupon discount curve is given just use the discounting curve
        // from the swap index. for rate calculation this curve cancels out in
        // the computation, so e.g. the discounting swap engine will produce
        // correct results, even if the couponDiscountCurve is not set here.
        // only the price member function in this class will be dependent on the
        // coupon discount curve.

        today_ = QuantLib::Settings::instance().evaluationDate();

        if (paymentDate_ > today_ && !couponDiscountCurve_.empty())
            couponDiscountRatio_ =
                couponDiscountCurve_->discount(paymentDate_) /
                discountCurve_->discount(paymentDate_);
        else
            couponDiscountRatio_ = 1.;

        spreadLegValue_ = spread_ * coupon_->accrualPeriod() *
                          discountCurve_->discount(paymentDate_) *
                          couponDiscountRatio_;

        if (fixingDate_ > today_) {

            swapTenor_ = swapIndex_->tenor();
            swap_ = swapIndex_->underlyingSwap(fixingDate_);

            swapRateValue_ = swap_->fairRate();
            annuity_ = 1.0E4 * std::fabs(swap_->fixedLegBPS());

            ext::shared_ptr<SmileSection> sectionTmp =
                swaptionVolatility()->smileSection(fixingDate_, swapTenor_);

            adjustedLowerBound_ = settings_.lowerRateBound_;
            adjustedUpperBound_ = settings_.upperRateBound_;

            if(sectionTmp->volatilityType() == Normal) {
                // adjust lower bound if it was not set explicitly
                if(settings_.defaultBounds_)
                    adjustedLowerBound_ = std::min(adjustedLowerBound_, -adjustedUpperBound_);
            } else {
                // adjust bounds by section's shift
                adjustedLowerBound_ -= sectionTmp->shift();
                adjustedUpperBound_ -= sectionTmp->shift();
            }

            // if the section does not provide an atm level, we enhance it to
            // have one, no need to exit with an exception ...

            if (sectionTmp->atmLevel() == Null<Real>())
                smileSection_ = ext::make_shared<AtmSmileSection>(
                    sectionTmp, swapRateValue_);
            else
                smileSection_ = sectionTmp;

            // compute linear model's parameters

            Real gx = 0.0, gy = 0.0;
            for (Size i = 0; i < swap_->fixedLeg().size(); i++) {
                ext::shared_ptr<Coupon> c =
                    ext::dynamic_pointer_cast<Coupon>(swap_->fixedLeg()[i]);
                Real yf = c->accrualPeriod();
                Date d = c->date();
                Real pv = yf * discountCurve_->discount(d);
                gx += pv * GsrG(d);
                gy += pv;
            }

            Real gamma = gx / gy;
            Date lastd = swap_->fixedLeg().back()->date();

            a_ = discountCurve_->discount(paymentDate_) *
                 (gamma - GsrG(paymentDate_)) /
                 (discountCurve_->discount(lastd) * GsrG(lastd) +
                  swapRateValue_ * gy * gamma);

            b_ = discountCurve_->discount(paymentDate_) / gy -
                 a_ * swapRateValue_;
        }
    }

    void AnalyticDoubleBarrierBinaryEngine::calculate() const {

        if (arguments_.barrierType == DoubleBarrier::KIKO ||
            arguments_.barrierType == DoubleBarrier::KOKI) {
            boost::shared_ptr<AmericanExercise> ex =
                boost::dynamic_pointer_cast<AmericanExercise>(
                                                   arguments_.exercise);
            QL_REQUIRE(ex, "KIKO/KOKI options must have American exercise");
            QL_REQUIRE(ex->dates()[0] <=
                       process_->blackVolatility()->referenceDate(),
                       "American option with window exercise not handled yet");
        } else {
            boost::shared_ptr<EuropeanExercise> ex =
                boost::dynamic_pointer_cast<EuropeanExercise>(
                                                   arguments_.exercise);
            QL_REQUIRE(ex, "non-European exercise given");
        }
        boost::shared_ptr<CashOrNothingPayoff> payoff =
            boost::dynamic_pointer_cast<CashOrNothingPayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "a cash-or-nothing payoff must be given");

        Real spot = process_->stateVariable()->value();
        QL_REQUIRE(spot > 0.0, "negative or null underlying given");

        Real variance =
            process_->blackVolatility()->blackVariance(
                                             arguments_.exercise->lastDate(),
                                             payoff->strike());
        Real barrier_lo = arguments_.barrier_lo;
        Real barrier_hi = arguments_.barrier_hi;
        DoubleBarrier::Type barrierType = arguments_.barrierType;
        QL_REQUIRE(barrier_lo>0.0,
                   "positive low barrier value required");
        QL_REQUIRE(barrier_hi>0.0,
                   "positive high barrier value required");
        QL_REQUIRE(barrier_lo < barrier_hi,
                   "barrier_lo must be < barrier_hi");
        QL_REQUIRE(barrierType == DoubleBarrier::KnockIn ||
                   barrierType == DoubleBarrier::KnockOut ||
                   barrierType == DoubleBarrier::KIKO ||
                   barrierType == DoubleBarrier::KOKI,
                   "Unsupported barrier type");

        // degenerate cases
        switch (barrierType) {
          case DoubleBarrier::KnockOut:
            if (spot <= barrier_lo || spot >= barrier_hi) {
                // knocked out, no value
                results_.value = 0;
                results_.delta = 0;
                results_.gamma = 0;
                results_.vega = 0;
                results_.rho = 0;
                return;
            }
            break;

          case DoubleBarrier::KnockIn:
            if (spot <= barrier_lo || spot >= barrier_hi) {
                // knocked in - pays
                results_.value = payoff->cashPayoff();
                results_.delta = 0;
                results_.gamma = 0;
                results_.vega = 0;
                results_.rho = 0;
                return;
            }
            break;

          case DoubleBarrier::KIKO:
            if (spot >= barrier_hi) {
                // knocked out, no value
                results_.value = 0;
                results_.delta = 0;
                results_.gamma = 0;
                results_.vega = 0;
                results_.rho = 0;
                return;
            } else if (spot <= barrier_lo) {
                // knocked in, pays
                results_.value = payoff->cashPayoff();
                results_.delta = 0;
                results_.gamma = 0;
                results_.vega = 0;
                results_.rho = 0;
                return;
            }
            break;

          case DoubleBarrier::KOKI:
            if (spot <= barrier_lo) {
                // knocked out, no value
                results_.value = 0;
                results_.delta = 0;
                results_.gamma = 0;
                results_.vega = 0;
                results_.rho = 0;
                return;
            } else if (spot >= barrier_hi) {
                // knocked in, pays
//.........这里部分代码省略.........

    Real DigitalPathPricer::operator()(const Path& path) const {
        Size n = path.length();
        QL_REQUIRE(n>1, "the path cannot be empty");

        Real log_asset_price = std::log(path.front());
        Real x, y;
        Volatility vol;
        TimeGrid timeGrid = path.timeGrid();
        Time dt;
        std::vector<Real> u = sequenceGen_.nextSequence().value;
        Real log_strike = std::log(payoff_->strike());

        Size i;
        switch (payoff_->optionType()) {
          case Option::Call:
            for (i=0; i<n-1; i++) {
                x = std::log(path[i+1]/path[i]);
                // terminal or initial vol?
                vol = diffProcess_->diffusion(timeGrid[i+1],
                                              std::exp(log_asset_price));
                // vol = diffProcess_->diffusion(timeGrid[i+2],
                //                               std::exp(log_asset_price+x));
                dt = timeGrid.dt(i);
                y = log_asset_price +
                    0.5*(x + std::sqrt(x*x-2*vol*vol*dt*std::log((1-u[i]))));
                // cross the strike
                if (y >= log_strike) {
                    if (exercise_->payoffAtExpiry()) {
                        return payoff_->cashPayoff() *
                            discountTS_->discount(path.timeGrid().back());
                    } else {
                        // the discount should be calculated at the exercise
                        // time between path.timeGrid()[i+1] and
                        // path.timeGrid()[i+2]
                        return payoff_->cashPayoff() *
                            discountTS_->discount(path.timeGrid()[i+1]);
                    }
                }
                log_asset_price += x;
            }
            break;
          case Option::Put:
            for (i=0; i<n-1; i++) {
                x = std::log(path[i+1]/path[i]);
                // terminal or initial vol?
                // initial (timeGrid[i+1]) for the time being
                vol = diffProcess_->diffusion(timeGrid[i+1],
                                              std::exp(log_asset_price));
                // vol = diffProcess_->diffusion(timeGrid[i+2],
                //                               std::exp(log_asset_price+x));
                dt = timeGrid.dt(i);
                y = log_asset_price +
                    0.5*(x - std::sqrt(x*x - 2*vol*vol*dt*std::log(u[i])));
                if (y <= log_strike) {
                    if (exercise_->payoffAtExpiry()) {
                        return payoff_->cashPayoff() *
                            discountTS_->discount(path.timeGrid().back());
                    } else {
                        // the discount should be calculated at the exercise
                        // time between path.timeGrid()[i+1] and
                        // path.timeGrid()[i+2]
                        return payoff_->cashPayoff() *
                            discountTS_->discount(path.timeGrid()[i+1]);
                    }
                }
                log_asset_price += x;
            }
            break;
          default:
            QL_FAIL("unknown option type");
        }

        return 0.0;
    }

    Rate YoYInflationIndex::fixing(const Date& fixingDate,
                                   bool /*forecastTodaysFixing*/) const {

        Date today = Settings::instance().evaluationDate();
        Date todayMinusLag = today - availabilityLag_;
        std::pair<Date,Date> lim = inflationPeriod(todayMinusLag, frequency_);
        Date lastFix = lim.first-1;

        Date flatMustForecastOn = lastFix+1;
        Date interpMustForecastOn = lastFix+1 - Period(frequency_);


        if (interpolated() && fixingDate >= interpMustForecastOn) {
            return forecastFixing(fixingDate);
        }

        if (!interpolated() && fixingDate >= flatMustForecastOn) {
            return forecastFixing(fixingDate);
        }

        // four cases with ratio() and interpolated()

        if (ratio()) {

            if(interpolated()){ // IS ratio, IS interpolated

                std::pair<Date,Date> lim = inflationPeriod(fixingDate, frequency_);
                Date fixMinus1Y=NullCalendar().advance(fixingDate, -1*Years, ModifiedFollowing);
                std::pair<Date,Date> limBef = inflationPeriod(fixMinus1Y, frequency_);
                Real dp= lim.second + 1 - lim.first;
                Real dpBef=limBef.second + 1 - limBef.first;
                Real dl = fixingDate-lim.first;
                // potentially does not work on 29th Feb
                Real dlBef = fixMinus1Y - limBef.first;
                // get the four relevant fixings
                // recall that they are stored flat for every day
                Rate limFirstFix =
                IndexManager::instance().getHistory(name())[lim.first];
                QL_REQUIRE(limFirstFix != Null<Rate>(),
                            "Missing " << name() << " fixing for "
                            << lim.first );
                Rate limSecondFix =
                IndexManager::instance().getHistory(name())[lim.second+1];
                QL_REQUIRE(limSecondFix != Null<Rate>(),
                            "Missing " << name() << " fixing for "
                            << lim.second+1 );
                Rate limBefFirstFix =
                IndexManager::instance().getHistory(name())[limBef.first];
                QL_REQUIRE(limBefFirstFix != Null<Rate>(),
                            "Missing " << name() << " fixing for "
                            << limBef.first );
                Rate limBefSecondFix =
                IndexManager::instance().getHistory(name())[limBef.second+1];
                QL_REQUIRE(limBefSecondFix != Null<Rate>(),
                            "Missing " << name() << " fixing for "
                            << limBef.second+1 );

                Real linearNow = limFirstFix + (limSecondFix-limFirstFix)*dl/dp;
                Real linearBef = limBefFirstFix + (limBefSecondFix-limBefFirstFix)*dlBef/dpBef;
                Rate wasYES = linearNow / linearBef - 1.0;

                return wasYES;

            } else {    // IS ratio, NOT interpolated
                Rate pastFixing =
                    IndexManager::instance().getHistory(name())[fixingDate];
                QL_REQUIRE(pastFixing != Null<Rate>(),
                            "Missing " << name() << " fixing for "
                            << fixingDate);
                Date previousDate = fixingDate - 1*Years;
                Rate previousFixing =
                IndexManager::instance().getHistory(name())[previousDate];
                QL_REQUIRE(previousFixing != Null<Rate>(),
                           "Missing " << name() << " fixing for "
                           << previousDate );

                return pastFixing/previousFixing - 1.0;
            }

        } else {  // NOT ratio

            if (interpolated()) { // NOT ratio, IS interpolated

                std::pair<Date,Date> lim = inflationPeriod(fixingDate, frequency_);
                Real dp= lim.second + 1 - lim.first;
                Real dl = fixingDate-lim.first;
                Rate limFirstFix =
                IndexManager::instance().getHistory(name())[lim.first];
                QL_REQUIRE(limFirstFix != Null<Rate>(),
                            "Missing " << name() << " fixing for "
                            << lim.first );
                Rate limSecondFix =
                IndexManager::instance().getHistory(name())[lim.second+1];
                QL_REQUIRE(limSecondFix != Null<Rate>(),
                            "Missing " << name() << " fixing for "
                            << lim.second+1 );
                Real linearNow = limFirstFix + (limSecondFix-limFirstFix)*dl/dp;

                return linearNow;

//.........这里部分代码省略.........

    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<PlainVanillaPayoff> payoff =
            boost::dynamic_pointer_cast<PlainVanillaPayoff>(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<Real>()(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");
        }
    }