[FEAT] Add one-year CDR standard errors to MackChainladder

chainladder/methods/mack.py

@@ -1,6 +1,7 @@
 # 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/.
+import numpy as np
 import pandas as pd
 from chainladder.methods import Chainladder
 
@@ -39,6 +40,12 @@ class MackChainladder(Chainladder):
         The total prediction error by origin period
     total_mack_std_err_:
         The total prediction error across all origin periods
+    cdr_std_err_:
+        The one-year standard error of the claims development result by
+        origin period (Merz and Wüthrich 2008, ``merz2008``)
+    total_cdr_std_err_:
+        The one-year standard error of the claims development result across
+        all origin periods
 
     See Also
     --------
@@ -465,6 +472,197 @@ def _get_total_mack_std_err_(self, obj):
             out = obj.sqrt().values[..., 0, -1]
         return pd.DataFrame(out, index=obj.index, columns=obj.columns)
 
+    @property
+    def cdr_std_err_(self):
+        """
+        One-year standard error of the claims development result by origin
+        period, Result 3.1 (eq. 3.9) of Merz and Wüthrich (2008, ``merz2008``).
+
+        The claims development result is the change in an origin's estimated
+        ultimate once the next diagonal is observed. ``mack_std_err_`` is the
+        counterpart to ultimate. The values follow the ``sigma_interpolation``
+        of :class:`Development`, and missing cells above the latest diagonal
+        give NaN.
+
+        Returns
+        -------
+        Triangle
+            One value per origin, shaped like ``ultimate_``, NaN for the
+            fully developed oldest origin.
+
+        Raises
+        ------
+        ValueError
+            If the triangle is not square, has a tail, or its development
+            factors are not volume-weighted over every link ratio.
+
+        Examples
+        --------
+        The triangle of Merz and Wüthrich (2008), with Mack's extrapolation
+        for the last sigma.
+
+        .. testsetup::
+
+            import chainladder as cl
+
+        .. testcode::
+
+            tr = cl.load_sample('mw2008')
+            dev = cl.Development(sigma_interpolation='mack').fit_transform(tr)
+            model = cl.MackChainladder().fit(dev)
+            print(model.cdr_std_err_)
+
+        .. testoutput::
+
+                          2261
+            2001           NaN
+            2002    566.174395
+            2003   1486.560344
+            2004   3923.098608
+            2005   9722.859763
+            2006  28442.621556
+            2007  20954.286973
+            2008  28119.317963
+            2009  53320.821049
+        """
+        obj = self.ultimate_.set_backend("numpy")
+        obj.values = np.sqrt(self._get_cdr_msep()[0])[..., None]
+        return obj.set_backend(self.X_.array_backend)
+
+    @property
+    def total_cdr_std_err_(self):
+        """
+        One-year standard error of the claims development result summed over
+        all origin periods, Result 3.3 (eq. 3.15) of Merz and Wüthrich (2008,
+        ``merz2008``). It exceeds the quadrature sum of ``cdr_std_err_``
+        because all origins share the estimated factors, and it raises
+        ``ValueError`` where ``cdr_std_err_`` does.
+
+        Returns
+        -------
+        DataFrame
+            One value per (index, column) combination of the fitted Triangle.
+
+        Examples
+        --------
+        On the paper's triangle the one-year error is about three quarters of
+        Mack's error to ultimate.
+
+        .. testsetup::
+
+            import chainladder as cl
+
+        .. testcode::
+
+            tr = cl.load_sample('mw2008')
+            dev = cl.Development(sigma_interpolation='mack').fit_transform(tr)
+            model = cl.MackChainladder().fit(dev)
+            print(model.total_cdr_std_err_)
+            print(model.total_mack_std_err_)
+
+        .. testoutput::
+
+            columns         values
+            (Total,)  81080.546787
+            columns          values
+            (Total,)  108401.387451
+        """
+        obj = self.X_
+        out = np.sqrt(self._get_cdr_msep()[1])
+        return pd.DataFrame(out, index=obj.index, columns=obj.columns)
+
+    def _get_cdr_msep(self):
+        """
+        Mean squared errors of the one-year claims development result.
+
+        Returns the per-origin values, shape (index, columns, origin), and
+        the totals over origins, shape (index, columns), from Results 3.1 and
+        3.3 of ``merz2008`` with I = J = n - 1. Origins i and development
+        periods j are 0-based, as in the paper, and the equation numbers in
+        the comments are the paper's.
+        """
+        from chainladder.development import Development
+
+        X = self.X_.set_backend("numpy")
+        C = X if X.is_cumulative else X.incr_to_cum()
+        a, b, n, m = C.shape
+        msg = "cdr_std_err_ and total_cdr_std_err_ need "
+        upper = np.add.outer(np.arange(n), np.arange(m)) <= n - 1
+        if n != m or not np.array_equal(~np.isnan(X.nan_triangle), upper):
+            raise ValueError(
+                msg + "a square triangle valued at its latest diagonal; this "
+                f"one has {n} origins and {m} development periods."
+            )
+        J = m - 1
+        ldf = X.ldf_.set_backend("numpy").values[..., 0, :]
+        sigma = X.sigma_.set_backend("numpy").values[..., 0, :]
+        if not (
+            np.all(np.nan_to_num(ldf[..., J:], nan=1.0) == 1)
+            and np.all(np.nan_to_num(sigma[..., J:]) == 0)
+        ):
+            raise ValueError(
+                msg + "a pattern without a tail, as Merz and Wüthrich (2008) "
+                "assume; fit without a tail estimator or with TailConstant(1.0)."
+            )
+        average = set(np.asarray(X.average_).ravel())
+        if average != {"volume"}:
+            raise ValueError(
+                msg + "volume-weighted development factors; this pattern uses "
+                f"{sorted(str(v) for v in average)}."
+            )
+        # A default Development on C gives the paper's factors (2.9), and the
+        # fitted w_ and ldf_ must match it.
+        fresh = Development().fit(C)
+        w, w_fresh = np.broadcast_arrays(np.asarray(X.w_), np.asarray(fresh.w_))
+        if not np.array_equal(w, w_fresh, equal_nan=True):
+            raise ValueError(
+                msg + "every link ratio of the triangle; this pattern excludes "
+                "some (drop, drop_high, drop_low, drop_valuation, n_periods)."
+            )
+        f = ldf[..., :J]  # f_j^I, (2.9)
+        if not np.allclose(
+            f, fresh.ldf_.values[..., 0, :], rtol=1e-8, atol=0, equal_nan=True
+        ):
+            raise ValueError(
+                msg + "development factors estimated from this triangle's "
+                "own link ratios; ldf_ differs from them, as after "
+                "Development(fillna=...)."
+            )
+        latest = C.latest_diagonal.values[..., 0]  # C_{i, I-i}
+        diag = latest[..., ::-1][..., :J]  # C_{I-j, j}
+        S1 = C.sum(axis="origin").values[..., 0, :J]  # S_j^{I+1}, (2.10)
+        S = S1 - diag  # S_j^I, (2.9)
+        ult = self.ultimate_.set_backend("numpy").values[..., 0]  # C_{i,J}^I, (2.11)
+        zero = np.zeros((a, b, 1))
+
+        def later(t):
+            """Sum of t over the positions after each position, last axis."""
+            return np.concatenate((np.cumsum(t[..., ::-1], -1)[..., -2::-1], zero), -1)
+
+        with np.errstate(divide="ignore", invalid="ignore"):
+            q = sigma[..., :J] ** 2 / f**2  # sigma_j^2 (3.3) over (f_j^I)^2
+            i = np.arange(1, n)  # origins with a next diagonal
+            d = J - i  # their next development period, I - i
+            qd, Sd, S1d, Ci = q[..., d], S[..., d], S1[..., d], latest[..., i]
+            # Sums over j = I-i+1 .. J-1: the summands of Phi (3.5), and those
+            # of the sum that Delta (3.4) and Lambda (3.13) share.
+            Phi = later(diag * q / S1**2)[..., d]
+            shared = later((diag / S1) ** 2 * q / S)[..., d]
+            Psi = qd / Ci  # (3.6)
+            Delta = qd / Sd + shared  # (3.4)
+            Lam = Ci / S1d * qd / Sd + shared  # Lambda, (3.13)
+            Xi = Phi + qd / S1d  # (3.14)
+            msep = np.full((a, b, n), np.nan)
+            # (3.9), with Gamma = Phi + Psi (3.7)
+            msep[..., 1:] = ult[..., i] ** 2 * (Phi + Psi + Delta)
+            # (3.15): 2 * sum over k > i > 0 of C_{i,J} C_{k,J} (Xi + Lambda)
+            cross = ult[..., i] * later(ult[..., 1:]) * (Xi + Lam)
+            total = msep[..., 1:].sum(-1) + 2 * cross.sum(-1)
+        complete = ~np.isnan(C.values[..., upper]).any(-1)
+        msep[~complete] = np.nan
+        total[~complete] = np.nan
+        return msep, total
+
     @property
     def summary_(self):
         """