feature: Estimate convergence factor alpha

Estimate convergence factor alpha. The parameter determines the convergence speed of the yield curve towards the Ultimate Forward Rate (UFR). The parameter is estimated by finding the smallest value such that the difference between forward rate at convergence maturity and UFR is smaller than 1bps.

Closes PR #4

.github/workflows/python-package.yml

@@ -76,6 +76,5 @@ jobs:
         TWINE_USERNAME: __token__
         TWINE_PASSWORD: ${{ secrets.PYPI_TOKEN }}
       run: |
-        twine upload -r testpypi dist/*
-
-    # twine upload dist/*
+        twine upload dist/*
+        # For testing only: twine upload -r testpypi dist/*

README.md

@@ -4,7 +4,8 @@ This Python package implements the Smith-Wilson yield curve fitting algorithm. I
 <br /><br />
 
 ## How to use the package
-1. To use the Smith-Wilson fitting algorithm, first import the Python package and specify the inputs. In the example below the inputs are zero-coupon rates with annual frequency up until year 25. The UFR is 2.9% and the convergence parameter alpha is 0.128562. The `terms` list defines the list of maturities, in this case `[1.0, 2.0, 3.0, ..., 25.0]`
+1. Install the package with `pip install smithwilson`
+2. To use the Smith-Wilson fitting algorithm, first import the Python package and specify the inputs. In the example below the inputs are zero-coupon rates with annual frequency up until year 25. The UFR is 2.9% and the convergence parameter alpha is 0.128562. The `terms` list defines the list of maturities, in this case `[1.0, 2.0, 3.0, ..., 25.0]`
     ```py
     import smithwilson as sw
 
@@ -22,7 +23,7 @@ This Python package implements the Smith-Wilson yield curve fitting algorithm. I
 
     ```
 
-1. Specify the targeted output maturities. This is the set of terms you want to get rates fitted by Smith-Wilson.
+3. Specify the targeted output maturities. This is the set of terms you want to get rates fitted by Smith-Wilson.
    Expand the set of rates beyond the Last Liquid Point (e.g. extrapolate to 150 years with annual frequency):
    ```py
    # Extrapolate to 150 years
@@ -41,14 +42,17 @@ This Python package implements the Smith-Wilson yield curve fitting algorithm. I
    terms_target = [0.25, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0]
    ```
 
-1. Call the Smiwth-Wilson fitting algorithm. This returns the rates as numpy vector with each element corresponding to the maturity in `terms_target`
+4. Call the Smiwth-Wilson fitting algorithm. This returns the rates as numpy vector with each element corresponding to the maturity in `terms_target`
    ```py
    # Calculate fitted rates based on actual observations and two parametes alpha & UFR
    fitted_rates = sw.fit_smithwilson_rates(rates_obs=rates, t_obs=terms,
-                                           t_target=terms_target, alpha=alpha, ufr=ufr)
+                                           t_target=terms_target, ufr=ufr,
+                                           alpha=alpha)  # Optional
    ```
 
-1. To display the results and/or processing them it can be useful to turn them into a table, here using the pandas library:
+   The convergence parameter alpha is optional and will be estimated if not provided. The parameter determines the convergence speed of the yield curve towards the Ultimate Forward Rate (UFR). The parameter is estimated by finding the smallest value such that the difference between forward rate at convergence maturity and UFR is smaller than 1bps.
+
+5. To display the results and/or processing them it can be useful to turn them into a table, here using the pandas library:
    ```py
    # Ensure pandas package is imported
    import pandas as pd
@@ -102,6 +106,3 @@ In the last case, `t` can be any maturity vector, i.e. with additional maturitie
 
 [EIOPA (2018). Technical documentation of the methodology to derive EIOPA’srisk-free interest rate term structures](https://eiopa.europa.eu/Publications/Standards/Technical%20Documentation%20(31%20Jan%202018).pdf); p.37-46
 <br /><br />
-
-## Author
-[Dejan Simic](https://www.linkedin.com/in/dejsimic/)

setup.py

@@ -5,9 +5,8 @@
 
 setuptools.setup(
     name="smithwilson",
-    version="0.1.0",
+    version="0.2.0",
     author="Dejan Simic",
-    # author_email="dejan.simic",
     description=
     "Implementation of the Smith-Wilson yield curve fitting algorithm in Python for interpolations and extrapolations of zero-coupon bond rates",
     long_description=long_description,
@@ -23,5 +22,5 @@
         "Operating System :: OS Independent",
     ],
     python_requires='>=3.7',
-    install_requires=["numpy>=1.21.5"],
+    install_requires=["numpy>=1.21.5", "scipy>=1.7.0"],
 )

smithwilson/__init__.py

@@ -1,2 +1 @@
-
-from smithwilson.core import *
+from smithwilson.core import calculate_prices, fit_convergence_parameter, fit_smithwilson_rates, ufr_discount_factor, fit_parameters, wilson_function

smithwilson/core.py

@@ -1,6 +1,7 @@
 from math import log
 import numpy as np
-from typing import Union, List
+from scipy import optimize
+from typing import Union, List, Optional
 
 
 def calculate_prices(rates: Union[np.ndarray, List[float]], t: Union[np.ndarray, List[float]]) -> np.ndarray:
@@ -49,7 +50,9 @@ def ufr_discount_factor(ufr: float, t: Union[np.ndarray, List[float]]) -> np.nda
     return np.exp(-ufr * t)
 
 
-def wilson_function(t1: Union[np.ndarray, List[float]], t2: Union[np.ndarray, List[float]], alpha: float, ufr: float) -> np.ndarray:
+def wilson_function(t1: Union[np.ndarray, List[float]],
+                    t2: Union[np.ndarray, List[float]],
+                    alpha: float, ufr: float) -> np.ndarray:
     """Calculate matrix of Wilson functions
 
     The Smith-Wilson method requires the calculation of a series of Wilson
@@ -94,7 +97,9 @@ def wilson_function(t1: Union[np.ndarray, List[float]], t2: Union[np.ndarray, Li
     return W
 
 
-def fit_parameters(rates: Union[np.ndarray, List[float]], t: Union[np.ndarray, List[float]], alpha: float, ufr: float) -> np.ndarray:
+def fit_parameters(rates: Union[np.ndarray, List[float]],
+                   t: Union[np.ndarray, List[float]],
+                   alpha: float, ufr: float) -> np.ndarray:
     """Calculate Smith-Wilson parameter vector ζ
 
     Given the Wilson-matrix, vector of discount factors and prices,
@@ -129,8 +134,10 @@ def fit_parameters(rates: Union[np.ndarray, List[float]], t: Union[np.ndarray, L
     return zeta
 
 
-def fit_smithwilson_rates(rates_obs: Union[np.ndarray, List[float]], t_obs: Union[np.ndarray, List[float]],
-                          t_target: Union[np.ndarray, List[float]], alpha: float, ufr: float) -> np.ndarray:
+def fit_smithwilson_rates(rates_obs: Union[np.ndarray, List[float]],
+                          t_obs: Union[np.ndarray, List[float]],
+                          t_target: Union[np.ndarray, List[float]],
+                          ufr: float, alpha: Optional[float] = None) -> np.ndarray:
     """Calculate zero-coupon yields with Smith-Wilson method based on observed rates.
 
     This function expects the rates and initial maturity vector to be
@@ -159,8 +166,9 @@ def fit_smithwilson_rates(rates_obs: Union[np.ndarray, List[float]], t_obs: Unio
         rates_obs: Initially observed zero-coupon rates vector before LLP of length n
         t_obs: Initially observed time to maturity vector (in years) of length n
         t_target: New targeted maturity vector (in years) with interpolated/extrapolated terms
-        alpha: Convergence speed parameter
         ufr: Ultimate Forward Rate (annualized/annual compounding)
+        alpha: (optional) Convergence speed parameter. If not provided estimated using
+            the `fit_convergence_parameter()` function
 
     Returns:
         Vector of zero-coupon rates with Smith-Wilson interpolated or extrapolated rates
@@ -173,6 +181,9 @@ def fit_smithwilson_rates(rates_obs: Union[np.ndarray, List[float]], t_obs: Unio
     t_obs = np.array(t_obs).reshape((-1, 1))
     t_target = np.array(t_target).reshape((-1, 1))
 
+    if alpha is None:
+        alpha = fit_convergence_parameter(rates_obs=rates_obs, t_obs=t_obs, ufr=ufr)
+
     zeta = fit_parameters(rates=rates_obs, t=t_obs, alpha=alpha, ufr=ufr)
     ufr_disc = ufr_discount_factor(ufr=ufr, t=t_target)
     W = wilson_function(t1=t_target, t2=t_obs, alpha=alpha, ufr=ufr)
@@ -183,3 +194,54 @@ def fit_smithwilson_rates(rates_obs: Union[np.ndarray, List[float]], t_obs: Unio
 
     # Transform price vector to zero-coupon rate vector (1/P)^(1/t) - 1
     return np.power(1 / P, 1 / t_target) - 1
+
+
+def fit_convergence_parameter(rates_obs: Union[np.ndarray, List[float]],
+                              t_obs: Union[np.ndarray, List[float]],
+                              ufr: float) -> float:
+    """Fit Smith-Wilson convergence factor (alpha).
+
+    Args:
+        rates_obs: Initially observed zero-coupon rates vector before LLP of length n
+        t_obs: Initially observed time to maturity vector (in years) of length n
+        ufr: Ultimate Forward Rate (annualized/annual compounding)
+
+    Returns:
+        Convergence parameter alpha
+    """
+
+    # Last liquid point (LLP)
+    llp = np.max(t_obs)
+
+    # Maturity at which forward curve is supposed to converge to ultimate forward rate (UFR)
+    # See: https://www.eiopa.europa.eu/sites/default/files/risk_free_interest_rate/12092019-technical_documentation.pdf (chapter 7.D., p. 39)
+    convergence_t = max(llp + 40, 60)
+
+    # Optimization function calculating the difference between UFR and forward rate at convergence point
+    def forward_difference(alpha: float):
+        # Fit yield curve
+        rates = fit_smithwilson_rates(rates_obs=rates_obs,     # Input rates to be fitted
+                                      t_obs=t_obs,             # Maturities of these rates
+                                      t_target=[convergence_t, convergence_t + 1],     # Maturity at which curve is supposed to converge to UFR
+                                      alpha=alpha,             # Optimization parameter
+                                      ufr=ufr)                 # Ultimate forward rate
+
+        # Calculate the forward rate at convergence maturity - this is an approximation since
+        # according to the documentation the minimization should be based on the forward intensity, not forward rate
+        forward_rate = (1 + rates[1])**(convergence_t + 1) / (1 + rates[0])**(convergence_t) - 1
+
+        # Absolute difference needs to be smaller than 1 bps
+        return -abs(forward_rate - ufr) + 1 / 10_000
+
+    # Minimize alpha w.r.t. forward difference criterion
+    root = optimize.minimize(lambda alpha: alpha, x0=0.15, method='SLSQP', bounds=[[0.05, 1.0]],
+                             constraints=[{
+                                 'type': 'ineq',
+                                 'fun': forward_difference
+                             }],
+                             options={
+                                 'ftol': 1e-6,
+                                 'disp': True
+                             })
+
+    return float(root.x)

smithwilson/tests/test_core.py

@@ -145,9 +145,10 @@ def test_fit_parameters(self):
 
     def test_fit_smithwilson_rates_actual(self):
         """Test estimation of yield curve fitted with the Smith-Wilson algorithm.
-           This example uses an actual example from EIOPA. Deviations must be less than 1bps (0.01%).
-           Source: https://eiopa.europa.eu/Publications/Standards/EIOPA_RFR_20190531.zip
-                   EIOPA_RFR_20190531_Term_Structures.xlsx; Tab: RFR_spot_no_VA; Switzerland
+
+        This example uses an actual example from EIOPA. Deviations must be less than 1bps (0.01%).
+        Source: https://eiopa.europa.eu/Publications/Standards/EIOPA_RFR_20190531.zip
+                EIOPA_RFR_20190531_Term_Structures.xlsx; Tab: RFR_spot_no_VA; Switzerland
         """
 
         # Input
@@ -186,11 +187,53 @@ def test_fit_smithwilson_rates_actual(self):
         np.testing.assert_almost_equal(actual, expected, decimal=4, err_msg="Fitted rates not matching")
 
 
-    def test_fit_smithwilson_rates_random(self):
-        """Test estimation of yield curve fitted with the Smith-Wilson algorithm
-           This test uses random data points.
+    def test_fit_smithwilson_rates_incl_convergence(self):
+        """Test estimation of yield curve without known convergence factor alpha.
+
+        This example uses an actual example from EIOPA. Deviations must be less than 1bps (0.01%).
+        Source: https://eiopa.europa.eu/Publications/Standards/EIOPA_RFR_20190531.zip
+                EIOPA_RFR_20190531_Term_Structures.xlsx; Tab: RFR_spot_no_VA; Switzerland
         """
 
+        # Input
+        r = np.array([-0.00803, -0.00814, -0.00778, -0.00725, -0.00652,
+                      -0.00565, -0.0048, -0.00391, -0.00313, -0.00214,
+                      -0.0014, -0.00067, -0.00008, 0.00051, 0.00108,
+                      0.00157, 0.00197, 0.00228, 0.0025, 0.00264,
+                      0.00271, 0.00274, 0.0028, 0.00291, 0.00309]).reshape((-1, 1))
+        t = np.array([float(y + 1) for y in range(len(r))]).reshape((-1, 1)) # 1.0, 2.0, ..., 25.0
+        ufr = 0.029
+        alpha = 0.128562
+
+        t_target = np.array([float(y + 1) for y in range(65)]).reshape((-1, 1))
+
+        # Expected Output
+        expected = np.array([-0.00803, -0.00814, -0.00778, -0.00725, -0.00652,
+                             -0.00565, -0.0048, -0.00391, -0.00313, -0.00214,
+                             -0.0014, -0.00067, -0.00008, 0.00051, 0.00108,
+                             0.00157, 0.00197, 0.00228, 0.0025, 0.00264,
+                             0.00271, 0.00274, 0.0028, 0.00291, 0.00309,
+                             0.00337, 0.00372, 0.00412, 0.00455, 0.00501,
+                             0.00548, 0.00596, 0.00644, 0.00692, 0.00739,
+                             0.00786, 0.00831, 0.00876, 0.00919, 0.00961,
+                             0.01002, 0.01042, 0.01081, 0.01118, 0.01154,
+                             0.01189, 0.01223, 0.01255, 0.01287, 0.01318,
+                             0.01347, 0.01376, 0.01403, 0.0143, 0.01456,
+                             0.01481, 0.01505, 0.01528, 0.01551, 0.01573,
+                             0.01594, 0.01615, 0.01635, 0.01655, 0.01673]).reshape((-1, 1))
+
+        # Actual Output
+        actual = sw.fit_smithwilson_rates(rates_obs=r, t_obs=t, t_target=t_target, ufr=ufr)
+
+        # Assert - Precision of 4 decimal points equals deviatino of less than 1bps
+        self.assertEqual(type(actual), type(expected), "Returned types not matching")
+        self.assertTupleEqual(actual.shape, expected.shape, "Shapes not matching")
+        np.testing.assert_almost_equal(actual, expected, decimal=4, err_msg="Fitted rates not matching")
+
+
+    def test_fit_smithwilson_rates_random(self):
+        """Test estimation of yield curve fitted with the Smith-Wilson algorithm using random data points."""
+
         # Input
         r = np.array([0.02, 0.025, -0.033, 0.01, 0.0008]).reshape((-1, 1))
         t = np.array([0.25, 1.0, 5.0, 20.0, 25.0]).reshape((-1, 1))
@@ -209,4 +252,33 @@ def test_fit_smithwilson_rates_random(self):
         # Assert
         self.assertEqual(type(actual), type(expected), "Returned types not matching")
         self.assertTupleEqual(actual.shape, expected.shape, "Shapes not matching")
-        np.testing.assert_almost_equal(actual, expected, decimal=8, err_msg="Fitted rates not matching")
+        np.testing.assert_almost_equal(actual, expected, decimal=8, err_msg="Fitted rates not matching")
+
+
+
+    def test_fit_alpha(self):
+        """Test estimation of convergence factor alpha.
+
+        This example uses an actual example from EIOPA. Deviations must be less than 0.001.
+        Source: https://eiopa.europa.eu/Publications/Standards/EIOPA_RFR_20190531.zip
+                EIOPA_RFR_20190531_Term_Structures.xlsx; Tab: RFR_spot_no_VA; Switzerland
+        """
+
+        # Input
+        r = np.array([-0.00803, -0.00814, -0.00778, -0.00725, -0.00652,
+                      -0.00565, -0.0048, -0.00391, -0.00313, -0.00214,
+                      -0.0014, -0.00067, -0.00008, 0.00051, 0.00108,
+                      0.00157, 0.00197, 0.00228, 0.0025, 0.00264,
+                      0.00271, 0.00274, 0.0028, 0.00291, 0.00309]).reshape((-1, 1))
+        t = np.array([float(y + 1) for y in range(len(r))]).reshape((-1, 1)) # 1.0, 2.0, ..., 25.0
+        ufr = 0.029
+
+        # Expected Output
+        alpha_expected = 0.128562
+
+        # Actual Output
+        alpha_actual = sw.fit_convergence_parameter(rates_obs=r, t_obs=t, ufr=ufr)
+
+        # Assert - Precision of 4 decimal points equals deviatino of less than 1bps
+        self.assertEqual(type(alpha_actual), type(alpha_expected), "Returned types not matching")
+        self.assertAlmostEqual(alpha_actual, alpha_expected, msg="Alpha not matching", delta=0.001)