diff --git a/ql/cashflows/cashflows.cpp b/ql/cashflows/cashflows.cpp index 46148d491f3..093ad45ba07 100644 --- a/ql/cashflows/cashflows.cpp +++ b/ql/cashflows/cashflows.cpp @@ -618,6 +618,9 @@ namespace QuantLib { Time t = 0.0; Date lastDate = npvDate; const DayCounter& dc = y.dayCounter(); + bool stc = (y.compounding() == SimpleThenCompounded); + DiscountFactor discount = 1.0; + bool firstCoupon = stc; for (const auto& i : leg) { if (i->hasOccurred(settlementDate, includeSettlementDateFlows)) continue; @@ -627,8 +630,22 @@ namespace QuantLib { c = 0.0; } - t += getStepwiseDiscountTime(i, dc, npvDate, lastDate); - DiscountFactor B = y.discountFactor(t); + Time dt = getStepwiseDiscountTime(i, dc, npvDate, lastDate); + t += dt; + DiscountFactor B; + if (stc) { + DiscountFactor b; + if (firstCoupon) { + b = 1.0 / (1.0 + y.rate() * dt); + firstCoupon = false; + } else { + b = y.discountFactor(dt); + } + discount *= b; + B = discount; + } else { + B = y.discountFactor(t); + } P += c * B; dPdy += t * c * B; @@ -660,6 +677,10 @@ namespace QuantLib { Natural N = y.frequency(); Date lastDate = npvDate; const DayCounter& dc = y.dayCounter(); + bool stc = (y.compounding() == SimpleThenCompounded); + DiscountFactor discount = 1.0; + Real alpha = 0.0; + bool firstCoupon = stc; for (const auto& i : leg) { if (i->hasOccurred(settlementDate, includeSettlementDateFlows)) continue; @@ -669,8 +690,24 @@ namespace QuantLib { c = 0.0; } - t += getStepwiseDiscountTime(i, dc, npvDate, lastDate); - DiscountFactor B = y.discountFactor(t); + Time dt = getStepwiseDiscountTime(i, dc, npvDate, lastDate); + t += dt; + DiscountFactor B; + if (stc) { + DiscountFactor b; + if (firstCoupon) { + b = 1.0 / (1.0 + r * dt); + alpha = dt * b; + firstCoupon = false; + } else { + b = y.discountFactor(dt); + alpha += dt / (1.0 + r / N); + } + discount *= b; + B = discount; + } else { + B = y.discountFactor(t); + } P += c * B; switch (y.compounding()) { case Simple: @@ -683,10 +720,7 @@ namespace QuantLib { dPdy -= c * B * t; break; case SimpleThenCompounded: - if (t<=1.0/N) - dPdy -= c * B*B * t; - else - dPdy -= c * t * B/(1+r/N); + dPdy -= c * B * alpha; break; case CompoundedThenSimple: if (t>1.0/N) @@ -833,6 +867,7 @@ namespace QuantLib { DiscountFactor discount = 1.0; Date lastDate = npvDate; const DayCounter& dc = y.dayCounter(); + bool firstCoupon = (y.compounding() == SimpleThenCompounded); for (const auto& i : leg) { if (i->hasOccurred(settlementDate, includeSettlementDateFlows)) continue; @@ -842,7 +877,14 @@ namespace QuantLib { amount = 0.0; } - DiscountFactor b = y.discountFactor(getStepwiseDiscountTime(i, dc, npvDate, lastDate)); + Time dt = getStepwiseDiscountTime(i, dc, npvDate, lastDate); + DiscountFactor b; + if (firstCoupon) { + b = 1.0 / (1.0 + y.rate() * dt); + firstCoupon = false; + } else { + b = y.discountFactor(dt); + } discount *= b; lastDate = i->date(); @@ -992,6 +1034,10 @@ namespace QuantLib { Rate r = y.rate(); Natural N = y.frequency(); Date lastDate = npvDate; + bool stc = (y.compounding() == SimpleThenCompounded); + DiscountFactor discount = 1.0; + Real alpha = 0.0, beta = 0.0; + bool firstCoupon = stc; for (const auto& i : leg) { if (i->hasOccurred(settlementDate, includeSettlementDateFlows)) continue; @@ -1001,8 +1047,26 @@ namespace QuantLib { c = 0.0; } - t += getStepwiseDiscountTime(i, dc, npvDate, lastDate); - DiscountFactor B = y.discountFactor(t); + Time dt = getStepwiseDiscountTime(i, dc, npvDate, lastDate); + t += dt; + DiscountFactor B; + if (stc) { + DiscountFactor b; + if (firstCoupon) { + b = 1.0 / (1.0 + r * dt); + alpha = dt * b; + beta = -dt * dt * b * b; + firstCoupon = false; + } else { + b = y.discountFactor(dt); + alpha += dt / (1.0 + r / N); + beta -= dt / (N * (1.0 + r / N) * (1.0 + r / N)); + } + discount *= b; + B = discount; + } else { + B = y.discountFactor(t); + } P += c * B; switch (y.compounding()) { case Simple: @@ -1015,10 +1079,7 @@ namespace QuantLib { d2Pdy2 += c * B*t*t; break; case SimpleThenCompounded: - if (t<=1.0/N) - d2Pdy2 += c * 2.0*B*B*B*t*t; - else - d2Pdy2 += c * B*t*(N*t+1)/(N*(1+r/N)*(1+r/N)); + d2Pdy2 += c * B * (alpha*alpha - beta); break; case CompoundedThenSimple: if (t>1.0/N) diff --git a/test-suite/bonds.cpp b/test-suite/bonds.cpp index 88f1956e7c3..0c38b3da60b 100644 --- a/test-suite/bonds.cpp +++ b/test-suite/bonds.cpp @@ -1892,6 +1892,60 @@ BOOST_AUTO_TEST_CASE(testFixingConvention) { BOOST_CHECK_EQUAL(couponF->fixingDate(), expectedFollowing); } +BOOST_AUTO_TEST_CASE(testSimpleThenCompoundedLongFirstPeriod) { + + BOOST_TEST_MESSAGE("Testing SimpleThenCompounded with a long first coupon period..."); + + // A semiannual bond with an 8-month first coupon period. The ICMA + // specification requires simple discounting for the first period + // regardless of its length; the time-magnitude threshold in + // InterestRate::compoundFactor() produces the wrong result in this case. + + Date settlement(15, January, 2023); + Settings::instance().evaluationDate() = settlement; + + // Backward schedule with an explicit first coupon date produces a long + // (8-month) first stub: Jan 15 -> Sep 15 -> Mar 15 -> ... -> Mar 15, 2026 + Schedule sch(settlement, + Date(15, March, 2026), + Period(Semiannual), + NullCalendar(), + Unadjusted, Unadjusted, + DateGeneration::Backward, + false, + Date(15, September, 2023)); + + DayCounter dc = Thirty360(Thirty360::BondBasis); + FixedRateBond bond(0, 100.0, sch, std::vector(1, 0.06), dc); + + InterestRate yield(0.05, dc, SimpleThenCompounded, Semiannual); + + // Expected dirty price computed from the ICMA formula: + // B1 = 1/(1 + 0.05*8/12) = 30/31 (simple for the long first period) + // B_k = B1 * (1/1.025)^(k-1) (compound for subsequent periods) + // cash flows: 4, 3, 3, 3, 3, 103 + Real expectedPrice = 102.8931; + Real tolerance = 1.0e-4; + + Real calculated = BondFunctions::dirtyPrice(bond, yield, settlement); + if (std::fabs(calculated - expectedPrice) > tolerance) + BOOST_ERROR("failed to reproduce dirty price for long-first-period bond" + << std::fixed << std::setprecision(6) + << "\n calculated: " << calculated + << "\n expected: " << expectedPrice + << "\n tolerance: " << tolerance); + + // Yield roundtrip: the implied yield must recover the input. + Bond::Price price{calculated, Bond::Price::Dirty}; + Rate impliedYield = BondFunctions::yield(bond, price, + dc, SimpleThenCompounded, Semiannual, + settlement, 1.0e-10, 100, 0.05); + if (std::fabs(impliedYield - yield.rate()) > 1.0e-7) + BOOST_ERROR("yield roundtrip failed for long-first-period bond" + << "\n recovered: " << impliedYield + << "\n expected: " << yield.rate()); +} + BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE_END()