Reconciling to Mack papers

chainladder/development/tests/test_development.py

@@ -814,46 +814,4 @@ def test_pipeline(clrd):
         dev2 = pipe.fit(X=clrd)
     assert np.array_equal(
         dev1.w_v2_.values, dev2.named_steps.drop_hilo.w_v2_.values, True
-    )
-
-def test_sigma1(atol):
-    data = {
-        "valuation": [
-            1981,1982,1983,1984,1985,
-            1982,1983,1984,1985,
-            1983,1984,1985,
-            1984,1985,
-            1985,
-        ],
-        "origin": [
-            1981,1981,1981,1981,1981,
-            1982,1982,1982,1982,
-            1983,1983,1983,
-            1984,1984,
-            1985,
-        ],
-        "values": [
-            1000,1385,1700,1905,2000,
-            1000,1395,1700,1895,
-            1000,1405,1700,
-            1000,1415,
-            1000,
-        ],
-    }
-    tri = cl.Triangle(
-        data,
-        origin="origin",
-        development="valuation",
-        columns=["values"],
-        cumulative=True,
-    )
-    assert np.allclose(
-        cl.Development(average='simple').fit(tri).sigma_.values,
-        cl.Development(average='volume').fit(tri).sigma_.values,
-        atol = atol
-    )        
-    assert np.allclose(
-        cl.Development(average='simple').fit(tri).sigma_.values,
-        cl.Development(average='regression').fit(tri).sigma_.values,
-        atol = atol
-    )      
+    )

chainladder/methods/tests/test_mack.py

@@ -27,4 +27,46 @@ def test_multi_triangle_mack(clrd,atol):
     mack = cl.MackChainladder().fit(tri)
     for i in range(len(tri.index)):
         for j in range(len(tri.columns)):
-            assert np.all(abs(mack.full_std_err_.iloc[i,j].values-cl.MackChainladder().fit(tri.iloc[i,j]).full_std_err_.values) < atol)
+            assert np.all(abs(mack.full_std_err_.iloc[i,j].values-cl.MackChainladder().fit(tri.iloc[i,j]).full_std_err_.values) < atol)
+
+def test_mack1997_hardcode():
+    """
+    Reconciles key MackChainladder statistics to values provided in the paper
+    """
+    #sourced from Table 1, p365 of Mack(1997)
+    ldf_se = [2.24,.517,.122,.051,.042,.023,.015,.012]
+    sigma = [1337,988.5,440.1,207.0,164.2,74.6,35.49,16.89]
+    #sourced from Table 2, p366 of Mack(1997)
+    ibnr_se = [0,61,140,319,596,1038,1298,1806,2182]
+
+    tri = cl.load_sample("mortgage")
+    dev = cl.Development(sigma_interpolation = 'mack').fit_transform(tri)
+    model = cl.MackChainladder().fit(dev)
+    ldf_rhs = dev.std_err_.values.flatten()
+    assert np.allclose(ldf_se[0],ldf_rhs[0],atol=0.01)
+    assert np.allclose(ldf_se[1:],ldf_rhs[1:],atol=0.001)
+    sigma_rhs = dev.sigma_.values.flatten()
+    assert np.allclose(sigma[0],sigma_rhs[0],atol=1)
+    assert np.allclose(sigma[1:],sigma_rhs[1:],atol=0.1)
+    ibnr_rhs = model.summary_.values[0,0,:,-1]/1000
+    assert np.allclose(ibnr_se,np.nan_to_num(ibnr_rhs,nan=0),atol=1,equal_nan=True)
+
+def test_mack1994_hardcode(raa):
+    """
+    Reconciles key MackChainladder statistics to values provided in the paper
+    """
+    #sourced from top table on p130 of Mack(1994)
+    sigma_sq = [27883,1109,691,61.2,119,40.8,1.34,7.88]
+    #sourced from bottom table on p130 of Mack(1994)
+    ibnr_se = [206,623,747,1469,2002,2209,5358,6333,24566]
+
+    dev = cl.Development(sigma_interpolation = 'mack').fit_transform(raa)
+    model = cl.MackChainladder().fit(dev)
+    sigma_rhs = dev.sigma_.values.flatten() ** 2
+    assert np.allclose(sigma_sq[:3],sigma_rhs[:3],atol=1)
+    assert np.allclose(sigma_sq[3:4],sigma_rhs[3:4],atol=0.1)
+    assert np.allclose(sigma_sq[4:5],sigma_rhs[4:5],atol=1)
+    assert np.allclose(sigma_sq[5:6],sigma_rhs[5:6],atol=0.1)
+    assert np.allclose(sigma_sq[6:8],sigma_rhs[6:8],atol=0.01)
+    ibnr_rhs = model.summary_.values[0,0,:,-1]
+    assert np.allclose(ibnr_se,ibnr_rhs[1:],atol=1)

chainladder/utils/weighted_regression.py

@@ -1,33 +1,94 @@
 # This Source Code Form is subject to the terms of the Mozilla Public
 # License, v. 2.0. If a copy of the MPL was not distributed with this
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
+from __future__ import annotations
+
 import numpy as np
 from chainladder.utils.sparse import sp
 from sklearn.base import BaseEstimator
 import warnings
 
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:  # pragma: no cover
+    from types import ModuleType
+    from typing import Literal
+    from chainladder.core.typing import BackendArray
 
 class WeightedRegression(BaseEstimator):
-    """Helper class that fits a system of regression equations
+    """
+    Helper class that fits a system of regression equations
     as a closed-form solution.  This greatly speeds up
     the implementation of the Mack stochastic properties.
+
+    Parameters
+    ----------
+    axis: integer (default = 2)
+        the axis along with the perform the regression; 
+        axis of 2 is along the origin periods; 
+        axis of 3 is along the development periods;
+    thru_orig: bool (default = False)
+        whether the regression is forced to go through the origin
+    xp: ModuleType (default = numpy)
+        array module for calculations
+
+    Attributes
+    ----------
+    slope_: Triangle
+        coefficients of the regression
+    sigma_: Triangle
+        standard deviation of the Triangle values according to Mack (1997), not be confused with sigma from the regression
+    std_err_: DatetimeIndex
+        standard error of estimator``slope_``
     """
 
-    def __init__(self, axis=None, thru_orig=False, xp=None, average=None):
+    def __init__(
+            self, 
+            axis: int = 2, 
+            thru_orig: bool = False, 
+            xp: ModuleType = np, 
+    ):
         self.axis = axis
         self.thru_orig = thru_orig
         self.xp = xp
-        self.average = average
 
     def infer_x_w(self):
+        """
+        Creates dummy X and/or w for the regression when none are given
+        """
         xp = self.xp
         if self.w is None:
             self.w = xp.ones(self.y.shape)
         if self.x is None:
             self.x = xp.cumsum(xp.ones(self.y.shape), self.axis)
         return self
 
-    def fit(self, X, y=None, sample_weight=None, average=None):
+    def fit(
+            self, 
+            X:BackendArray, 
+            y:BackendArray|None=None, 
+            sample_weight:BackendArray|None=None, 
+            average: Literal["volume", "simple", "regression", "geometric"] | None = None
+    ):
+        """
+        Fit the model with X.
+
+        Parameters
+        ----------
+        X : Array
+            independent variable for the regression
+        y : Array or None (default = None)
+            dependent variable for the regression
+        sample_weight : Array or None (default = None)
+            which (x,y) pairs should be used in the regression
+        average: literal (or list of literals), or None, optional (default = None)
+            type of averaging to use; dictates the weights used in the regression
+
+        Returns
+        -------
+        self : object
+            Returns the instance itself.
+        """
         self.x = X
         self.y = y
         self.w = sample_weight
@@ -44,6 +105,10 @@ def fit(self, X, y=None, sample_weight=None, average=None):
         return self
 
     def _fit_OLS_thru_orig(self):
+        """
+        Given a set of w, x, y, and an axis, this Function returns OLS slope
+        and other statistics, while forcing an intercept of 0
+        """
         from chainladder.utils.utility_functions import num_to_nan
 
         x = self.x
@@ -106,19 +171,14 @@ def _fit_OLS_thru_orig(self):
         mse_denom = xp.nansum((y * 0 + 1) * (xp.nan_to_num(w) != 0), axis) - 1
         mse_denom = num_to_nan(mse_denom)
         mse = wss_residual / mse_denom
-
+        sigma = xp.sqrt(mse)
         std_err = xp.sqrt(mse / denominator)
-        sigma = std_err * xp.sqrt(mse_denom + 1)
-
-        coef = coef[..., None]
-        sigma = sigma[..., None]
-        std_err = std_err[..., None]
 
         self._w_reg = w
 
-        self.slope_ = coef
-        self.sigma_ = sigma
-        self.std_err_ = std_err
+        self.slope_ = coef[..., None]
+        self.sigma_ = sigma[..., None]
+        self.std_err_ = std_err[..., None]
 
         return self