ENH: amend the oaspread algo to higher accuracy, lower efficiency

rateslib/instruments.py

@@ -37,9 +37,9 @@
 
 from rateslib import defaults
 from rateslib.default import NoInput
-from rateslib.calendars import add_tenor, get_calendar, dcf, _get_years_and_months
+from rateslib.calendars import add_tenor, get_calendar, dcf, _get_years_and_months, _DCF1d
 
-from rateslib.curves import Curve, index_left, LineCurve, CompositeCurve, IndexCurve
+from rateslib.curves import Curve, index_left, LineCurve, CompositeCurve, IndexCurve, average_rate
 from rateslib.solver import Solver
 from rateslib.periods import (
     Cashflow,
@@ -1424,6 +1424,7 @@ def fwd_from_repo(
         repo_rate: Union[float, Dual, Dual2],
         convention: Union[str, NoInput] = NoInput(0),
         dirty: bool = False,
+        method: str = "proceeds"
     ):
         """
         Return a forward price implied by a given repo rate.
@@ -1443,6 +1444,8 @@ def fwd_from_repo(
             values.
         dirty : bool, optional
             Whether the input and output price are specified including accrued interest.
+        method : str in {"proceeds", "compounded"}, optional
+            The method for determining the forward price.
 
         Returns
         -------
@@ -1484,9 +1487,17 @@ def fwd_from_repo(
 
         for p_idx in range(settlement_idx, fwd_settlement_idx):
             # deduct accrued coupon from dirty price
-            dcf_ = dcf(self.leg1.periods[p_idx].payment, forward_settlement, convention)
-            accrued_coup = self.leg1.periods[p_idx].cashflow * (1 + dcf_ * repo_rate / 100)
-            total_rtn -= accrued_coup
+            if method.lower() == "proceeds":
+                dcf_ = dcf(self.leg1.periods[p_idx].payment, forward_settlement, convention)
+                accrued_coup = self.leg1.periods[p_idx].cashflow * (1 + dcf_ * repo_rate / 100)
+                total_rtn -= accrued_coup
+            elif method.lower() == "compounded":
+                r_bar, d, _ = average_rate(settlement, forward_settlement, convention, repo_rate)
+                n = (forward_settlement - self.leg1.periods[p_idx].payment).days
+                accrued_coup = self.leg1.periods[p_idx].cashflow * (1 + d * r_bar / 100) ** n
+                total_rtn -= accrued_coup
+            else:
+                raise ValueError("`method` must be in {'proceeds', 'compounded'}.")
 
         forward_price = total_rtn / -self.leg1.notional * 100
         if dirty:
@@ -1850,27 +1861,36 @@ def oaspread(
             self.curves, solver, curves, fx, base, self.leg1.currency
         )
         ad_ = curves[1].ad
-        curves[1]._set_ad_order(2)
-        disc_curve = curves[1].shift(Dual2(0, "z_spread"), composite=False)
-        curves[1]._set_ad_order(0)
         metric = "dirty_price" if dirty else "clean_price"
+
+        curves[1]._set_ad_order(1)
+        disc_curve = curves[1].shift(Dual(0, "z_spread"), composite=False)
         npv_price = self.rate(curves=[curves[0], disc_curve], metric=metric)
 
+        # find a first order approximation of z
+        b = npv_price.gradient("z_spread", 1)[0]
+        c = float(npv_price) - float(price)
+        z_hat = -c / b
+
+        # shift the curve to the first order approximation and fine tune with 2nd order approxim.
+        curves[1]._set_ad_order(2)
+        disc_curve = curves[1].shift(Dual2(z_hat, "z_spread"), composite=False)
+        npv_price = self.rate(curves=[curves[0], disc_curve], metric=metric)
         a, b = (
             0.5 * npv_price.gradient("z_spread", 2)[0][0],
             npv_price.gradient("z_spread", 1)[0],
         )
-        z = _quadratic_equation(a, b, float(npv_price) - float(price))
-        # first z is solved by using 1st and 2nd derivatives to get close to target NPV
+        z_hat2 = _quadratic_equation(a, b, float(npv_price) - float(price))
 
-        # TODO (low) add a tolerance here to continually converge to the solution, via GradDes?
-        disc_curve = curves[1].shift(z, composite=False)
+        # perform one final approximation albeit the additional price calculation slows calc time
+        curves[1]._set_ad_order(0)
+        disc_curve = curves[1].shift(z_hat+z_hat2, composite=False)
         npv_price = self.rate(curves=[curves[0], disc_curve], metric=metric)
-        diff = npv_price - price
-        new_b = b + 2 * a * z
-        z = z - diff / new_b
-        # then a final linear adjustment is made which is usually very small
+        b = b + 2 * a * z_hat2  # forecast the new gradient
+        c = float(npv_price) - float(price)
+        z_hat3 = -c / b
 
+        z = z_hat + z_hat2 + z_hat3
         curves[1]._set_ad_order(ad_)
         return z
 

tests/test_instruments_bonds.py

@@ -747,10 +747,10 @@ def test_fixed_rate_bond_implied_repo_analogue_dirty(self, f_s, f_p):
         assert abs(result - 1.0) < 1e-8
 
     @pytest.mark.parametrize("price, tol", [
-        (112.0, 1e-6),
-        (104.0, 1e-5),
-        (96.0, 1e-3),
-        (91.0, 1e-2)
+        (112.0, 1e-10),
+        (104.0, 1e-10),
+        (96.0, 1e-9),
+        (91.0, 1e-7)
     ])
     def test_oaspread(self, price, tol):
         gilt = FixedRateBond(
@@ -772,6 +772,34 @@ def test_oaspread(self, price, tol):
         result = gilt.rate(curve_z, metric="clean_price")
         assert abs(result - price) < tol
 
+    @pytest.mark.parametrize("price, tol", [
+        (85, 1e-8),
+        (75, 1e-6),
+        (65, 1e-4),
+        (55, 1e-3),
+        (45, 1e-1),
+        (35, 0.20),
+    ])
+    def test_oaspread_low_price(self, price, tol):
+        gilt = FixedRateBond(
+            effective=dt(1998, 12, 7),
+            termination=dt(2015, 12, 7),
+            frequency="S",
+            calendar="ldn",
+            currency="gbp",
+            convention="ActActICMA",
+            ex_div=7,
+            fixed_rate=1.0,
+            notional=-100,
+            settle=0,
+        )
+        curve = Curve({dt(1999, 11, 25): 1.0, dt(2015, 12, 7): 0.85})
+        # result = gilt.npv(curve) = 113.22198344812742
+        result = gilt.oaspread(curve, price=price)
+        curve_z = curve.shift(result, composite=False)
+        result = gilt.rate(curve_z, metric="clean_price")
+        assert abs(result - price) < tol
+
     def test_cashflows_no_curve(self):
         gilt = FixedRateBond(
             effective=dt(2001, 1, 1),