C++ DayCounter类代码示例

void DayCounterTest::testThirty360_BondBasis() {

    BOOST_MESSAGE("Testing thirty/360 day counter (Bond Basis)...");

    // http://www.isda.org/c_and_a/docs/30-360-2006ISDADefs.xls
    // Source: 2006 ISDA Definitions, Sec. 4.16 (f)
    // 30/360 (or Bond Basis)

    DayCounter dayCounter = Thirty360(Thirty360::BondBasis);
    std::vector<Date> testStartDates;
    std::vector<Date> testEndDates;
    Time calculated;

    // ISDA - Example 1: End dates do not involve the last day of February
    testStartDates.push_back(Date(20, August, 2006)); testEndDates.push_back(Date(20, February, 2007));
    testStartDates.push_back(Date(20, February, 2007)); testEndDates.push_back(Date(20, August, 2007));
    testStartDates.push_back(Date(20, August, 2007)); testEndDates.push_back(Date(20, February, 2008));
    testStartDates.push_back(Date(20, February, 2008)); testEndDates.push_back(Date(20, August, 2008));
    testStartDates.push_back(Date(20, August, 2008)); testEndDates.push_back(Date(20, February, 2009));
    testStartDates.push_back(Date(20, February, 2009)); testEndDates.push_back(Date(20, August, 2009));

    // ISDA - Example 2: End dates include some end-February dates
    testStartDates.push_back(Date(31, August, 2006)); testEndDates.push_back(Date(28, February, 2007));
    testStartDates.push_back(Date(28, February, 2007)); testEndDates.push_back(Date(31, August, 2007));
    testStartDates.push_back(Date(31, August, 2007)); testEndDates.push_back(Date(29, February, 2008));
    testStartDates.push_back(Date(29, February, 2008)); testEndDates.push_back(Date(31, August, 2008));
    testStartDates.push_back(Date(31, August, 2008)); testEndDates.push_back(Date(28, February, 2009));
    testStartDates.push_back(Date(28, February, 2009)); testEndDates.push_back(Date(31, August, 2009));

    //// ISDA - Example 3: Miscellaneous calculations
    testStartDates.push_back(Date(31, January, 2006)); testEndDates.push_back(Date(28, February, 2006));
    testStartDates.push_back(Date(30, January, 2006)); testEndDates.push_back(Date(28, February, 2006));
    testStartDates.push_back(Date(28, February, 2006)); testEndDates.push_back(Date(3, March, 2006));
    testStartDates.push_back(Date(14, February, 2006)); testEndDates.push_back(Date(28, February, 2006));
    testStartDates.push_back(Date(30, September, 2006)); testEndDates.push_back(Date(31, October, 2006));
    testStartDates.push_back(Date(31, October, 2006)); testEndDates.push_back(Date(28, November, 2006));
    testStartDates.push_back(Date(31, August, 2007)); testEndDates.push_back(Date(28, February, 2008));
    testStartDates.push_back(Date(28, February, 2008)); testEndDates.push_back(Date(28, August, 2008));
    testStartDates.push_back(Date(28, February, 2008)); testEndDates.push_back(Date(30, August, 2008));
    testStartDates.push_back(Date(28, February, 2008)); testEndDates.push_back(Date(31, August, 2008));
    testStartDates.push_back(Date(26, February, 2007)); testEndDates.push_back(Date(28, February, 2008));
    testStartDates.push_back(Date(26, February, 2007)); testEndDates.push_back(Date(29, February, 2008));
    testStartDates.push_back(Date(29, February, 2008)); testEndDates.push_back(Date(28, February, 2009));
    testStartDates.push_back(Date(28, February, 2008)); testEndDates.push_back(Date(30, March, 2008));
    testStartDates.push_back(Date(28, February, 2008)); testEndDates.push_back(Date(31, March, 2008));

    int expected[] = { 180, 180, 180, 180, 180, 180,
                       178, 183, 179, 182, 178, 183,
                        28,  28,   5,  14,  30,  28,
                       178, 180, 182, 183, 362, 363,
                       359,  32,  33};

    for (Size i = 0; i < testStartDates.size(); i++) {
        calculated = dayCounter.dayCount(testStartDates[i], testEndDates[i]);
        if (calculated != expected[i]) {
                BOOST_ERROR("from " << testStartDates[i]
                            << " to " << testEndDates[i] << ":\n"
                            << "    calculated: " << calculated << "\n"
                            << "    expected:   " << expected[i-1]);
        }
    }
}

    void AnalyticDividendEuropeanEngine::calculate() const {

        QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
                   "not an European option");

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

        Date settlementDate = process_->riskFreeRate()->referenceDate();
        Real riskless = 0.0;
        Size i;
        for (i=0; i<arguments_.cashFlow.size(); i++)
            if (arguments_.cashFlow[i]->date() >= settlementDate)
                riskless += arguments_.cashFlow[i]->amount() *
                    process_->riskFreeRate()->discount(
                                              arguments_.cashFlow[i]->date());

        Real spot = process_->stateVariable()->value() - riskless;
        QL_REQUIRE(spot > 0.0,
                   "negative or null underlying after subtracting dividends");

        DiscountFactor dividendDiscount =
            process_->dividendYield()->discount(
                                             arguments_.exercise->lastDate());
        DiscountFactor riskFreeDiscount =
            process_->riskFreeRate()->discount(arguments_.exercise->lastDate());
        Real forwardPrice = spot * dividendDiscount / riskFreeDiscount;

        Real variance =
            process_->blackVolatility()->blackVariance(
                                              arguments_.exercise->lastDate(),
                                              payoff->strike());

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

        results_.value = black.value();
        results_.delta = black.delta(spot);
        results_.gamma = black.gamma(spot);

        DayCounter rfdc  = process_->riskFreeRate()->dayCounter();
        DayCounter voldc = process_->blackVolatility()->dayCounter();
        Time t = voldc.yearFraction(
                                 process_->blackVolatility()->referenceDate(),
                                 arguments_.exercise->lastDate());
        results_.vega = black.vega(t);

        Real delta_theta = 0.0, delta_rho = 0.0;
        for (i = 0; i < arguments_.cashFlow.size(); i++) {
            Date d = arguments_.cashFlow[i]->date();
            if (d >= settlementDate) {
                delta_theta -= arguments_.cashFlow[i]->amount() *
                  process_->riskFreeRate()->zeroRate(d,rfdc,Continuous,Annual)*
                  process_->riskFreeRate()->discount(d);
                Time t = process_->time(d);
                delta_rho += arguments_.cashFlow[i]->amount() * t *
                             process_->riskFreeRate()->discount(t);
            }
        }
        t = process_->time(arguments_.exercise->lastDate());
        try {
            results_.theta = black.theta(spot, t) +
                             delta_theta * black.delta(spot);
        } catch (Error&) {
            results_.theta = Null<Real>();
        }

        results_.rho = black.rho(t) +
                       delta_rho * black.delta(spot);
    }

 inline bool operator==(const DayCounter& d1, const DayCounter& d2) {
     return (d1.empty() && d2.empty())
         || (!d1.empty() && !d2.empty() && d1.name() == d2.name());
 }

    void BinomialVanillaEngine<T>::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<PlainVanillaPayoff> payoff =
            boost::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-plain payoff given");

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

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

        TimeGrid grid(maturity, timeSteps_);

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

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

        DiscretizedVanillaOption option(arguments_, *process_, grid);

        option.initialize(lattice, maturity);

        // Partial derivatives calculated from various points in the
        // binomial tree (Odegaard)

        // Rollback to third-last step, and get underlying price (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 p2h = va2[2]; // high-price
        Real s2 = lattice->underlying(2, 2); // high price

        // Rollback to second-last step, and get option value (p1) at
        // this point
        option.rollback(grid[1]);
        Array va(option.values());
        QL_ENSURE(va.size() == 2, "Expect 2 nodes in grid at first step");
        Real p1 = va[1];

        // Finally, rollback to t=0
        option.rollback(0.0);
        Real p0 = option.presentValue();
        Real s1 = lattice->underlying(1, 1);

        // Calculate partial derivatives
        Real delta0 = (p1-p0)/(s1-s0);   // dp/ds
        Real delta1 = (p2h-p1)/(s2-s1);  // dp/ds

        // Store results
        results_.value = p0;
        results_.delta = delta0;
        results_.gamma = 2.0*(delta1-delta0)/(s2-s0);    //d(delta)/ds
        results_.theta = blackScholesTheta(process_,
                                           results_.value,
                                           results_.delta,
                                           results_.gamma);
    }

    void ContinuousArithmeticAsianLevyEngine::calculate() const {
        QL_REQUIRE(arguments_.averageType == Average::Arithmetic,
                   "not an Arithmetic average option");
        QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
                   "not an European Option");

        DayCounter rfdc  = process_->riskFreeRate()->dayCounter();
        DayCounter divdc = process_->dividendYield()->dayCounter();
        DayCounter voldc = process_->blackVolatility()->dayCounter();
        Real spot = process_->stateVariable()->value();

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

        // original time to maturity
        Date maturity = arguments_.exercise->lastDate();
        Time T = rfdc.yearFraction(startDate_,
                                   arguments_.exercise->lastDate());
        // remaining time to maturity
        Time T2 = rfdc.yearFraction(process_->riskFreeRate()->referenceDate(),
                                    arguments_.exercise->lastDate());

        Real strike = payoff->strike();

        Volatility volatility =
            process_->blackVolatility()->blackVol(maturity, strike);

        CumulativeNormalDistribution N;

        Rate riskFreeRate = process_->riskFreeRate()->
            zeroRate(maturity, rfdc, Continuous, NoFrequency);
        Rate dividendYield = process_->dividendYield()->
            zeroRate(maturity, divdc, Continuous, NoFrequency);
        Real b = riskFreeRate - dividendYield;
        QL_REQUIRE(b != 0.0, "null cost of carry not allowed by Levy engine");

        Real Se = (spot/(T*b))*(exp((b-riskFreeRate)*T2)-exp(-riskFreeRate*T2));

        Real X;
        if (T2 < T) {
            QL_REQUIRE(!currentAverage_.empty() && currentAverage_->isValid(),
                       "current average required");
            X = strike - ((T-T2)/T)*currentAverage_->value();
        } else {
            X = strike;
        }

        Real M = (2*spot*spot/(b+volatility*volatility)) *
            (((exp((2*b+volatility*volatility)*T2)-1)
              / (2*b+volatility*volatility))-((exp(b*T2)-1)/b));

        Real D = M/(T*T);

        Real V = log(D)-2*(riskFreeRate*T2+log(Se));

        Real d1 = (1/sqrt(V))*((log(D)/2)-log(X));
        Real d2 = d1-sqrt(V);

        if(payoff->optionType()==Option::Call)
            results_.value = Se*N(d1) - X*exp(-riskFreeRate*T2)*N(d2);
        else
            results_.value = Se*N(d1) - X*exp(-riskFreeRate*T2)*N(d2)
                             - Se + X*exp(-riskFreeRate*T2);
    }

void MargrabeOptionTest::testGreeks() {

    BOOST_MESSAGE("Testing analytic European exchange option greeks...");

    SavedSettings backup;

    std::map<std::string,Real> calculated, expected, tolerance;
    tolerance["delta1"]  = 1.0e-5;
    tolerance["delta2"]  = 1.0e-5;
    tolerance["gamma1"]  = 1.0e-5;
    tolerance["gamma2"]  = 1.0e-5;
    tolerance["theta"]   = 1.0e-5;
    tolerance["rho"]     = 1.0e-5;

    Option::Type types[] = {Option::Call};
    Real underlyings1[]  = { 22.0 };
    Real underlyings2[]  = { 20.0 };
    Rate qRates1[]       = { 0.06, 0.16, 0.04 };
    Rate qRates2[]       = { 0.04, 0.14, 0.02 };
    Rate rRates[]        = { 0.1, 0.2, 0.08 };
    Time residualTimes[] = { 0.1, 0.5 };
    Volatility vols1[]   = { 0.20 };
    Volatility vols2[]   = { 0.15, 0.20, 0.25};

    DayCounter dc = Actual360();
    Date today = Date::todaysDate();
    Settings::instance().evaluationDate() = today;

    boost::shared_ptr<SimpleQuote> spot1(new SimpleQuote(0.0));
    boost::shared_ptr<SimpleQuote> spot2(new SimpleQuote(0.0));

    boost::shared_ptr<SimpleQuote> qRate1(new SimpleQuote(0.0));
    boost::shared_ptr<YieldTermStructure> qTS1 = flatRate(qRate1, dc);
    boost::shared_ptr<SimpleQuote> qRate2(new SimpleQuote(0.0));
    boost::shared_ptr<YieldTermStructure> qTS2 = flatRate(qRate2, dc);

    boost::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
    boost::shared_ptr<YieldTermStructure> rTS = flatRate(rRate, dc);

    boost::shared_ptr<SimpleQuote> vol1(new SimpleQuote(0.0));
    boost::shared_ptr<BlackVolTermStructure> volTS1 = flatVol(vol1, dc);
    boost::shared_ptr<SimpleQuote> vol2(new SimpleQuote(0.0));
    boost::shared_ptr<BlackVolTermStructure> volTS2 = flatVol(vol2, dc);

    for (Size i=0; i<LENGTH(types); i++) {
        for (Size k=0; k<LENGTH(residualTimes); k++) {
          Date exDate = today + timeToDays(residualTimes[k]);
          boost::shared_ptr<Exercise> exercise(new EuropeanExercise(exDate));

          // option to check
          boost::shared_ptr<BlackScholesMertonProcess> stochProcess1(new
              BlackScholesMertonProcess(Handle<Quote>(spot1),
                                      Handle<YieldTermStructure>(qTS1),
                                      Handle<YieldTermStructure>(rTS),
                                      Handle<BlackVolTermStructure>(volTS1)));

          boost::shared_ptr<BlackScholesMertonProcess> stochProcess2(new
              BlackScholesMertonProcess(Handle<Quote>(spot2),
                                      Handle<YieldTermStructure>(qTS2),
                                      Handle<YieldTermStructure>(rTS),
                                      Handle<BlackVolTermStructure>(volTS2)));

          std::vector<boost::shared_ptr<StochasticProcess1D> > procs;
          procs.push_back(stochProcess1);
          procs.push_back(stochProcess2);

          //The correlation -0.5 can be different real between -1 and 1 for more tests
          Real correlation = -0.5;
          Matrix correlationMatrix(2,2, correlation);
          for (Integer j=0; j < 2; j++) {
              correlationMatrix[j][j] = 1.0;

          boost::shared_ptr<PricingEngine> engine(
                             new AnalyticEuropeanMargrabeEngine(stochProcess1,
                                                                stochProcess2,
                                                                correlation));

          //The quantities of S1 and S2 can be different from 1 & 1 for more tests
          MargrabeOption margrabeOption(1, 1, exercise);

          // analytic engine
          margrabeOption.setPricingEngine(engine);

              for (Size l=0; l<LENGTH(underlyings1); l++) {
                for (Size m=0; m<LENGTH(qRates1); m++) {
                  for (Size n=0; n<LENGTH(rRates); n++) {
                    for (Size p=0; p<LENGTH(vols1); p++) {
                      Real u1 = underlyings1[l],
                           u2 = underlyings2[l],
                           u;
                      Rate q1 = qRates1[m],
                           q2 = qRates2[m],
                           r = rRates[n];
                      Volatility v1 = vols1[p],
                                 v2 = vols2[p];

                      spot1 ->setValue(u1);
                      spot2 ->setValue(u2);
                      qRate1->setValue(q1);
                      qRate2->setValue(q2);
//.........这里部分代码省略.........

void AnalyticKouEuropeanEngine::calculate() const
{

	QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
		"not an European option");

	boost::shared_ptr<StrikedTypePayoff> payoff =
		boost::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);

	QL_REQUIRE(payoff, "non-striked payoff given");

	/*
	*calculate to be implemented
	*/
	Option::Type type = payoff->optionType();
	Real eta1 = process_->posJumpMean();
	QL_REQUIRE(eta1 > 1, "invalid eta1 input!!");
	Real eta2 = process_->negJumpMean();
	QL_REQUIRE(eta2 > 0, "invalid eta2 input!!");
	Real lambda = process_->jumpIntensity();
	Real p = process_->posProbability();

	//here assume the three term structure have the same daycounter
	DayCounter daycounter = process_->riskFreeRate()->dayCounter();
	Real sigma = process_->blackVolatility()->blackVol(arguments_.exercise->lastDate(), payoff->strike());
	Real variance = sigma*sigma;
	// the paper ignore the case when dividends are payed out
	Real effectiveR =
		process_->riskFreeRate()->zeroRate(arguments_.exercise->lastDate(), daycounter, Continuous)
		- process_->dividendYield()->zeroRate(arguments_.exercise->lastDate(), daycounter, Continuous);
	Real S = process_->stateVariable()->value();
	QL_REQUIRE(S > 0.0, "negative or null underlying given");
	Real K = payoff->strike();
	Time T = daycounter.yearFraction(process_->riskFreeRate()->referenceDate(),
		arguments_.exercise->lastDate());

	/*
	*auxillary variables
	*/
	Real zeta = p*eta1 / (eta1 - 1) + (1 - p)*eta2 / (eta2 + 1) - 1;
	Real adjustEta1 = eta1 - 1;
	Real adjustEta2 = eta2 + 1;
	Real adjustLambda = process_->jumpIntensity()*(zeta + 1);
	Real effectiveP = p / (1 + zeta)*eta1 / (eta1 - 1);

	//price of call option
	Real price =
		S*process_->dividendYield()->discount(arguments_.exercise->lastDate())
		*KouHelper::Gamma(effectiveR + 0.5*sigma*sigma - lambda*zeta, sigma, adjustLambda, effectiveP, adjustEta1, adjustEta2, log(K / S), T, tolerance_) -
		K*process_->riskFreeRate()->discount(arguments_.exercise->lastDate())
		*KouHelper::Gamma(effectiveR - 0.5*sigma*sigma - lambda*zeta, sigma, lambda, p, eta1, eta2, log(K / S), T, tolerance_);

	if (type == Option::Type::Call)
	{
		results_.value = price;
	}

	else results_.value = price - S*process_->dividendYield()->discount(arguments_.exercise->lastDate())
		+ K*process_->riskFreeRate()->discount(arguments_.exercise->lastDate());

	//more results if needed
}

void SwingOptionTest::testExtOUJumpSwingOption() {

    BOOST_TEST_MESSAGE("Testing simple swing option pricing for Kluge model...");

    SavedSettings backup;

    Date settlementDate = Date::todaysDate();
    Settings::instance().evaluationDate() = settlementDate;
    DayCounter dayCounter = ActualActual();
    Date maturityDate = settlementDate + Period(12, Months);

    boost::shared_ptr<StrikedTypePayoff> payoff(
                                     new PlainVanillaPayoff(Option::Put, 30));

    std::vector<Date> exerciseDates(1, settlementDate+Period(1, Months));
    while (exerciseDates.back() < maturityDate) {
        exerciseDates.push_back(exerciseDates.back()+Period(1, Months));
    }
    boost::shared_ptr<SwingExercise> swingExercise(
                                            new SwingExercise(exerciseDates));

    std::vector<Time> exerciseTimes(exerciseDates.size());
    for (Size i=0; i < exerciseTimes.size(); ++i) {
        exerciseTimes[i]
                 = dayCounter.yearFraction(settlementDate, exerciseDates[i]);
    }

    TimeGrid grid(exerciseTimes.begin(), exerciseTimes.end(), 60);
    std::vector<Size> exerciseIndex(exerciseDates.size());
    for (Size i=0; i < exerciseIndex.size(); ++i) {
        exerciseIndex[i] = grid.closestIndex(exerciseTimes[i]);
    }

    boost::shared_ptr<ExtOUWithJumpsProcess> jumpProcess = createKlugeProcess();

    const Rate irRate = 0.1;
    boost::shared_ptr<YieldTermStructure> rTS(
                                flatRate(settlementDate, irRate, dayCounter));

    boost::shared_ptr<PricingEngine> swingEngine(
                new FdSimpleExtOUJumpSwingEngine(jumpProcess, rTS, 25, 50, 25));

    boost::shared_ptr<PricingEngine> vanillaEngine(
                new FdExtOUJumpVanillaEngine(jumpProcess, rTS, 25, 50, 25));

    VanillaOption bermudanOption(payoff, swingExercise);
    bermudanOption.setPricingEngine(vanillaEngine);
    const Real bermudanOptionPrices = bermudanOption.NPV();

    const Size nrTrails = 16000;
    typedef PseudoRandom::rsg_type rsg_type;
    typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
    rsg_type rsg = PseudoRandom::make_sequence_generator(
                    jumpProcess->factors()*(grid.size()-1), BigNatural(421));

    MultiPathGenerator<rsg_type> generator(jumpProcess, grid, rsg, false);

    for (Size i=0; i < exerciseDates.size(); ++i) {
        const Size exerciseRights = i+1;

        VanillaSwingOption swingOption(payoff, swingExercise,
                                       exerciseRights, exerciseRights);
        swingOption.setPricingEngine(swingEngine);
        const Real swingOptionPrice = swingOption.NPV();

        const Real upperBound = exerciseRights*bermudanOptionPrices;

        if (swingOptionPrice - upperBound > 2e-2) {
            BOOST_ERROR("Failed to reproduce upper bounds"
                        << "\n    upper Bound: " << upperBound
                        << "\n    Price:       " << swingOptionPrice);
        }

        Real lowerBound = 0.0;
        for (Size j=exerciseDates.size()-i-1; j < exerciseDates.size(); ++j) {
            VanillaOption europeanOption(payoff, boost::shared_ptr<Exercise>(
                                     new EuropeanExercise(exerciseDates[j])));
            europeanOption.setPricingEngine(
                boost::shared_ptr<PricingEngine>(vanillaEngine));
            lowerBound += europeanOption.NPV();
        }

        if (lowerBound - swingOptionPrice > 2e-2) {
            BOOST_ERROR("Failed to reproduce lower bounds"
                       << "\n    lower Bound: " << lowerBound
                       << "\n    Price:       " << swingOptionPrice);
        }

        // use MC plus perfect forecast to find an upper bound
        GeneralStatistics npv;
        for (Size n=0; n < nrTrails; ++n) {
            sample_type path = generator.next();

            std::vector<Real> exerciseValues(exerciseTimes.size());
            for (Size k=0; k < exerciseTimes.size(); ++k) {
                const Real x = path.value[0][exerciseIndex[k]];
                const Real y = path.value[1][exerciseIndex[k]];
                const Real s = std::exp(x+y);

                exerciseValues[k] =(*payoff)(s)*rTS->discount(exerciseDates[k]);
//.........这里部分代码省略.........

void DefaultProbabilityCurveTest::testDefaultProbability() {

    BOOST_TEST_MESSAGE("Testing default-probability structure...");

    Real hazardRate = 0.0100;
    Handle<Quote> hazardRateQuote = Handle<Quote>(
                boost::shared_ptr<Quote>(new SimpleQuote(hazardRate)));
    DayCounter dayCounter = Actual360();
    Calendar calendar = TARGET();
    Size n = 20;

    double tolerance = 1.0e-10;
    Date today = Settings::instance().evaluationDate();
    Date startDate = today;
    Date endDate = startDate;

    FlatHazardRate flatHazardRate(startDate, hazardRateQuote, dayCounter);

    for(Size i=0; i<n; i++){
        startDate = endDate;
        endDate = calendar.advance(endDate, 1, Years);

        Probability pStart = flatHazardRate.defaultProbability(startDate);
        Probability pEnd = flatHazardRate.defaultProbability(endDate);

        Probability pBetweenComputed =
            flatHazardRate.defaultProbability(startDate, endDate);

        Probability pBetween = pEnd - pStart;

        if (std::fabs(pBetween - pBetweenComputed) > tolerance)
            BOOST_ERROR(
                "Failed to reproduce probability(d1, d2) "
                << "for default probability structure\n"
                << std::setprecision(12)
                << "    calculated probability: " << pBetweenComputed << "\n"
                << "    expected probability:   " << pBetween);

        Time t2 = dayCounter.yearFraction(today, endDate);
        Probability timeProbability = flatHazardRate.defaultProbability(t2);
        Probability dateProbability =
            flatHazardRate.defaultProbability(endDate);

        if (std::fabs(timeProbability - dateProbability) > tolerance)
            BOOST_ERROR(
                "single-time probability and single-date probability do not match\n"
                << std::setprecision(10)
                << "    time probability: " << timeProbability << "\n"
                << "    date probability: " << dateProbability);

        Time t1 = dayCounter.yearFraction(today, startDate);
        timeProbability = flatHazardRate.defaultProbability(t1, t2);
        dateProbability = flatHazardRate.defaultProbability(startDate, endDate);

        if (std::fabs(timeProbability - dateProbability) > tolerance)
            BOOST_ERROR(
                "double-time probability and double-date probability do not match\n"
                << std::setprecision(10)
                << "    time probability: " << timeProbability << "\n"
                << "    date probability: " << dateProbability);
    }
}

void BatesModelTest::testAnalyticVsBlack() {

    BOOST_MESSAGE("Testing analytic Bates engine against Black formula...");

    SavedSettings backup;

    Date settlementDate = Date::todaysDate();
    Settings::instance().evaluationDate() = settlementDate;

    DayCounter dayCounter = ActualActual();
    Date exerciseDate = settlementDate + 6*Months;

    boost::shared_ptr<StrikedTypePayoff> payoff(
                                     new PlainVanillaPayoff(Option::Put, 30));
    boost::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));

    Handle<YieldTermStructure> riskFreeTS(flatRate(0.1, dayCounter));
    Handle<YieldTermStructure> dividendTS(flatRate(0.04, dayCounter));
    Handle<Quote> s0(boost::shared_ptr<Quote>(new SimpleQuote(32.0)));

    Real yearFraction = dayCounter.yearFraction(settlementDate, exerciseDate);
    Real forwardPrice = s0->value()*std::exp((0.1-0.04)*yearFraction);
    Real expected = blackFormula(payoff->optionType(), payoff->strike(),
        forwardPrice, std::sqrt(0.05*yearFraction)) *
                                            std::exp(-0.1*yearFraction);
    const Real v0 = 0.05;
    const Real kappa = 5.0;
    const Real theta = 0.05;
    const Real sigma = 1.0e-4;
    const Real rho = 0.0;
    const Real lambda = 0.0001;
    const Real nu = 0.0; 
    const Real delta = 0.0001;

    VanillaOption option(payoff, exercise);

    boost::shared_ptr<BatesProcess> process(
        new BatesProcess(riskFreeTS, dividendTS, s0, v0, 
                         kappa, theta, sigma, rho, lambda, nu, delta));

    boost::shared_ptr<PricingEngine> engine(new BatesEngine(
        boost::shared_ptr<BatesModel>(new BatesModel(process)), 64));

    option.setPricingEngine(engine);
    Real calculated = option.NPV();

    Real tolerance = 2.0e-7;
    Real error = std::fabs(calculated - expected);
    if (error > tolerance) {
        BOOST_ERROR("failed to reproduce Black price with BatesEngine"
                    << QL_FIXED
                    << "\n    calculated: " << calculated
                    << "\n    expected:   " << expected
                    << QL_SCIENTIFIC
                    << "\n    error:      " << error);
    }

    engine = boost::shared_ptr<PricingEngine>(new BatesDetJumpEngine(
        boost::shared_ptr<BatesDetJumpModel>(
            new BatesDetJumpModel( process, 1.0, 0.0001)), 64));

    option.setPricingEngine(engine);
    calculated = option.NPV();

    error = std::fabs(calculated - expected);
    if (error > tolerance) {
        BOOST_ERROR("failed to reproduce Black price with " \
                    "BatesDetJumpEngine"
                    << QL_FIXED
                    << "\n    calculated: " << calculated
                    << "\n    expected:   " << expected
                    << QL_SCIENTIFIC
                    << "\n    error:      " << error);
    }

    engine = boost::shared_ptr<PricingEngine>(new BatesDoubleExpEngine(
        boost::shared_ptr<BatesDoubleExpModel>(
            new BatesDoubleExpModel(process, 0.0001, 0.0001, 0.0001)), 64));

    option.setPricingEngine(engine);
    calculated = option.NPV();

    error = std::fabs(calculated - expected);
    if (error > tolerance) {
        BOOST_ERROR("failed to reproduce Black price with BatesDoubleExpEngine"
                    << QL_FIXED
                    << "\n    calculated: " << calculated
                    << "\n    expected:   " << expected
                    << QL_SCIENTIFIC
                    << "\n    error:      " << error);
    }

    engine = boost::shared_ptr<PricingEngine>(new BatesDoubleExpDetJumpEngine(
        boost::shared_ptr<BatesDoubleExpDetJumpModel>(
            new BatesDoubleExpDetJumpModel(
                process, 0.0001, 0.0001, 0.0001, 0.5, 1.0, 0.0001)), 64));

    option.setPricingEngine(engine);
    calculated = option.NPV();

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

    void FdHullWhiteSwaptionEngine::calculate() const {
        const Handle<YieldTermStructure> ts = model_->termStructure();

        // 1. Layout
        const boost::shared_ptr<FdmLinearOpLayout> layout(
            new FdmLinearOpLayout(std::vector<Size>(1u, xGrid_)));

        // 2. Mesher
        const DayCounter dc = ts->dayCounter();
        const Date referenceDate = ts->referenceDate();
        const Time maturity = dc.yearFraction(referenceDate,
                                              arguments_.exercise->lastDate());


        const boost::shared_ptr<OrnsteinUhlenbeckProcess> process(
            new OrnsteinUhlenbeckProcess(model_->a(), model_->sigma()));

        const boost::shared_ptr<Fdm1dMesher> shortRateMesher(
            new FdmSimpleProcess1dMesher(xGrid_, process, maturity,1,invEps_));

        const boost::shared_ptr<FdmMesher> mesher(
            new FdmMesherComposite(layout,
                std::vector<boost::shared_ptr<Fdm1dMesher> >(
                    1, shortRateMesher)));

        // 3. Inner Value Calculator
        const std::vector<Date>& exerciseDates = arguments_.exercise->dates();
        std::map<Time, Date> t2d;

        for (Size i=0; i < exerciseDates.size(); ++i) {
            const Time t = dc.yearFraction(referenceDate, exerciseDates[i]);
            QL_REQUIRE(t >= 0, "exercise dates must not contain past date");

            t2d[t] = exerciseDates[i];
        }

        const Handle<YieldTermStructure> disTs = model_->termStructure();
        const Handle<YieldTermStructure> fwdTs
            = arguments_.swap->iborIndex()->forwardingTermStructure();

        QL_REQUIRE(fwdTs->dayCounter() == disTs->dayCounter(),
                "day counter of forward and discount curve must match");
        QL_REQUIRE(fwdTs->referenceDate() == disTs->referenceDate(),
                "reference date of forward and discount curve must match");

        const boost::shared_ptr<HullWhite> fwdModel(
            new HullWhite(fwdTs, model_->a(), model_->sigma()));

        const boost::shared_ptr<FdmInnerValueCalculator> calculator(
             new FdmAffineModelSwapInnerValue<HullWhite>(
                 model_.currentLink(), fwdModel,
                 arguments_.swap, t2d, mesher, 0));

        // 4. Step conditions
        const boost::shared_ptr<FdmStepConditionComposite> conditions =
             FdmStepConditionComposite::vanillaComposite(
                 DividendSchedule(), arguments_.exercise,
                 mesher, calculator, referenceDate, dc);

        // 5. Boundary conditions
        const std::vector<boost::shared_ptr<FdmDirichletBoundary> > boundaries;

        // 6. Solver
        FdmSolverDesc solverDesc = { mesher, boundaries, conditions,
                                     calculator, maturity,
                                     tGrid_, dampingSteps_ };

        const boost::scoped_ptr<FdmHullWhiteSolver> solver(
            new FdmHullWhiteSolver(model_, solverDesc, schemeDesc_));

        results_.value = solver->valueAt(0.0);
    }

    void AnalyticCliquetEngine::calculate() const {

        QL_REQUIRE(arguments_.accruedCoupon == Null<Real>() &&
                   arguments_.lastFixing == Null<Real>(),
                   "this engine cannot price options already started");
        QL_REQUIRE(arguments_.localCap == Null<Real>() &&
                   arguments_.localFloor == Null<Real>() &&
                   arguments_.globalCap == Null<Real>() &&
                   arguments_.globalFloor == Null<Real>(),
                   "this engine cannot price capped/floored options");

        QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
                   "not an European option");

        boost::shared_ptr<PercentageStrikePayoff> moneyness =
            boost::dynamic_pointer_cast<PercentageStrikePayoff>(
                                                           arguments_.payoff);
        QL_REQUIRE(moneyness, "wrong payoff given");

        std::vector<Date> resetDates = arguments_.resetDates;
        resetDates.push_back(arguments_.exercise->lastDate());

        Real underlying = process_->stateVariable()->value();
        QL_REQUIRE(underlying > 0.0, "negative or null underlying");
        Real strike = underlying * moneyness->strike();
        boost::shared_ptr<StrikedTypePayoff> payoff(
                      new PlainVanillaPayoff(moneyness->optionType(),strike));

        results_.value = 0.0;
        results_.delta = results_.gamma = 0.0;
        results_.theta = 0.0;
        results_.rho = results_.dividendRho = 0.0;
        results_.vega = 0.0;

        for (Size i = 1; i < resetDates.size(); i++) {

            Real weight =
                process_->dividendYield()->discount(resetDates[i-1]);
            DiscountFactor discount =
                process_->riskFreeRate()->discount(resetDates[i]) /
                process_->riskFreeRate()->discount(resetDates[i-1]);
            DiscountFactor qDiscount =
                process_->dividendYield()->discount(resetDates[i]) /
                process_->dividendYield()->discount(resetDates[i-1]);
            Real forward = underlying*qDiscount/discount;
            Real variance =
                process_->blackVolatility()->blackForwardVariance(
                                        resetDates[i-1],resetDates[i],strike);

            BlackCalculator black(payoff, forward, std::sqrt(variance), discount);

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

            results_.value += weight * black.value();
            results_.delta += weight * (black.delta(underlying) +
                                        moneyness->strike() * discount *
                                        black.beta());
            results_.gamma += 0.0;
            results_.theta += process_->dividendYield()->forwardRate(
                resetDates[i-1], resetDates[i], rfdc, Continuous, NoFrequency) *
                weight * black.value();

            Time dt = rfdc.yearFraction(resetDates[i-1],resetDates[i]);
            results_.rho += weight * black.rho(dt);

            Time t = divdc.yearFraction(
                                    process_->dividendYield()->referenceDate(),
                                    resetDates[i-1]);
            dt = divdc.yearFraction(resetDates[i-1],resetDates[i]);
            results_.dividendRho += weight * (black.dividendRho(dt) -
                                              t * black.value());

            dt = voldc.yearFraction(resetDates[i-1], resetDates[i]);
            results_.vega += weight * black.vega(dt);
        }

    }

    void BjerksundStenslandApproximationEngine::calculate() const {

        QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
                   "not an American Option");

        boost::shared_ptr<AmericanExercise> ex =
            boost::dynamic_pointer_cast<AmericanExercise>(arguments_.exercise);
        QL_REQUIRE(ex, "non-American exercise given");
        QL_REQUIRE(!ex->payoffAtExpiry(),
                   "payoff at expiry not handled");

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

        Real variance =
            process_->blackVolatility()->blackVariance(ex->lastDate(),
                                                       payoff->strike());
        DiscountFactor dividendDiscount =
            process_->dividendYield()->discount(ex->lastDate());
        DiscountFactor riskFreeDiscount =
            process_->riskFreeRate()->discount(ex->lastDate());
        Real spot = process_->stateVariable()->value();
        QL_REQUIRE(spot > 0.0, "negative or null underlying given");
        Real strike = payoff->strike();

        if (payoff->optionType()==Option::Put) {
            // use put-call simmetry
            std::swap(spot, strike);
            std::swap(riskFreeDiscount, dividendDiscount);
            payoff = boost::shared_ptr<PlainVanillaPayoff>(
                                new PlainVanillaPayoff(Option::Call, strike));
        }

        if (dividendDiscount>=1.0) {
            // early exercise is never optimal - use Black formula
            Real forwardPrice = spot * dividendDiscount / riskFreeDiscount;
            BlackCalculator black(payoff, forwardPrice, std::sqrt(variance),
                                  riskFreeDiscount);

            results_.value        = black.value();
            results_.delta        = black.delta(spot);
            results_.deltaForward = black.deltaForward();
            results_.elasticity   = black.elasticity(spot);
            results_.gamma        = black.gamma(spot);

            DayCounter rfdc  = process_->riskFreeRate()->dayCounter();
            DayCounter divdc = process_->dividendYield()->dayCounter();
            DayCounter voldc = process_->blackVolatility()->dayCounter();
            Time t =
                rfdc.yearFraction(process_->riskFreeRate()->referenceDate(),
                                  arguments_.exercise->lastDate());
            results_.rho = black.rho(t);

            t = divdc.yearFraction(process_->dividendYield()->referenceDate(),
                                   arguments_.exercise->lastDate());
            results_.dividendRho = black.dividendRho(t);

            t = voldc.yearFraction(process_->blackVolatility()->referenceDate(),
                                   arguments_.exercise->lastDate());
            results_.vega        = black.vega(t);
            results_.theta       = black.theta(spot, t);
            results_.thetaPerDay = black.thetaPerDay(spot, t);

            results_.strikeSensitivity  = black.strikeSensitivity();
            results_.itmCashProbability = black.itmCashProbability();
        } else {
            // early exercise can be optimal - use approximation
            results_.value = americanCallApproximation(spot,
                                                       strike,
                                                       riskFreeDiscount,
                                                       dividendDiscount,
                                                       variance);
        }
    }

void qlHazardRateFactory::calculate()
{
    try {

        TiXmlDocument document;
        document.Parse(instInfo_.c_str(), 0, TIXML_ENCODING_UTF8);
        TiXmlNode* rootNode = document.FirstChild("defaultCurveCalculation");

        Settings::instance().evaluationDate() =
            strToDate(rootNode->FirstChildElement("evaluationTime")->GetText());

        DayCounter daycounter = Actual365Fixed();

        TiXmlElement* pricingParametorNode = rootNode->FirstChildElement("pricingParameter");
        TiXmlElement* referenceCurveNode = pricingParametorNode->FirstChildElement("referenceCurves");

        Calendar calendar = SouthKorea();

        Date refereceDate = Settings::instance().evaluationDate();

        // must be a business day
        refereceDate = calendar.adjust(refereceDate);

        // dummy curve
        boost::shared_ptr<Quote> flatRate(new SimpleQuote(0.01));

        //Handle<YieldTermStructure> tsCurve(
        //      boost::shared_ptr<FlatForward>(
        //              new FlatForward(todaysDate, Handle<Quote>(flatRate),
        //                              Actual365Fixed())));


        TiXmlNode* pricingParametorNode = document.FirstChild("pricingParameter");

        TiXmlElement* referenceCurveNode = pricingParametorNode->FirstChildElement("referenceCurves");

        std::vector<boost::shared_ptr<HazardRateStructure>> hazardRateCurves;

        TiXmlElement* tsNode = referenceCurveNode->FirstChildElement("termstructure");

        if(tsNode)
        {
            for(tsNode; tsNode; tsNode = tsNode->NextSiblingElement("termstructure"))
            {
                hazardRateCurves.push_back(this->hazardRateTS(tsNode));
            }
        }

        TiXmlNode* yieldNode = document.FirstChild("yieldNode");

        qlYieldTermStructureFactory yieldTsFactory = qlYieldTermStructureFactory();
        this->tsCurve_ = yieldTsFactory.yieldTSHandle(yieldNode);

        // outputTenor Build
        std::vector<Period> outputTenors;
        Size outputNum = outputTenors.size();

        qlTimeUnitFactory tuFactory = qlTimeUnitFactory();
        TiXmlElement* outputTenorNode = pricingParametorNode->FirstChildElement("outputCurveTemplate");

        TiXmlElement* tenorNode = outputTenorNode->FirstChildElement("tenorList");

        if(tenorNode)
        {
            for(tenorNode; tenorNode; tenorNode = tenorNode->NextSiblingElement("tenor"))
            {
                outputTenors.push_back(tuFactory.timeUnit(tenorNode));
            }
        }

        //compositeCalculation

        boost::shared_ptr<DefaultProbabilityTermStructure> copulaHazardTS =
            boost::shared_ptr<DefaultProbabilityTermStructure>(new CopulaHazardRate<ProbabilityOfAtLeastNEvents,Linear>(
                        dates,hazardRateCurves,0.9,1,50));



        //result


        for (Size i=0 ; i<outputNum ; ++i)
        {
            dates_.push_back(calendar.advance(refereceDate,outputTenors[i]));
            dateTimes_.push_back(daycounter.yearFraction(refereceDate,dates_[i]));
            data_.push_back(copulaHazardTS->defaultProbability(dateTimes_[i]));

        }

        // makeXml

        this->buildResultXml();

        //this->buildResultXml(inst_->additionalResults());

    } catch (std::exception& e) {
        std::cerr << e.what() << std::endl;
    } catch (...) {
        std::cerr << "unknown error" << std::endl;
    }
//.........这里部分代码省略.........

    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
//.........这里部分代码省略.........