[BRK] Count the factor pairs present in DevelopmentCorrelation

chainladder/core/correlation.py

@@ -38,12 +38,13 @@ class DevelopmentCorrelation:
     ----------
     t_critical: DataFrame
         Boolean value for whether correlation is too high at significance level
-        ``p_critical``.
+        ``p_critical``. False where ``t_expectation`` is NaN.
     t_expectation: DataFrame
-        Values representing the Spearman rank correlation
-    t_variance: float
-        Variance measure of Spearman rank correlation
-    confidence_interval: tuple
+        Values representing the Spearman rank correlation. NaN for a triangle
+        with no column of at least two factor pairs.
+    t_variance: DataFrame
+        Variance measure of Spearman rank correlation, one value per triangle.
+    confidence_interval: tuple of DataFrame
         Range around zero that contains ``t_expectation`` with probability of
         about ``1 - p_critical`` when the development factors are uncorrelated.
         ``t_critical`` is True when ``t_expectation`` falls outside it.
@@ -76,8 +77,8 @@ class DevelopmentCorrelation:
         raa = cl.load_sample('raa')
         dc = cl.DevelopmentCorrelation(raa, p_critical=0.5)
         print(round(float(dc.t_expectation.iloc[0, 0]), 4))
-        print(round(float(dc.confidence_interval[0]), 4))
-        print(round(float(dc.confidence_interval[1]), 4))
+        print(round(float(dc.confidence_interval[0].iloc[0, 0]), 4))
+        print(round(float(dc.confidence_interval[1].iloc[0, 0]), 4))
         print(bool(dc.t_critical.iloc[0, 0]))
 
     .. testoutput::
@@ -100,8 +101,8 @@ class DevelopmentCorrelation:
     .. testcode::
 
         tight = cl.DevelopmentCorrelation(raa, p_critical=0.99)
-        print(round(float(tight.confidence_interval[0]), 4))
-        print(round(float(tight.confidence_interval[1]), 4))
+        print(round(float(tight.confidence_interval[0].iloc[0, 0]), 4))
+        print(round(float(tight.confidence_interval[1].iloc[0, 0]), 4))
         print(bool(tight.t_critical.iloc[0, 0]))
 
     .. testoutput::
@@ -121,70 +122,79 @@ def __init__(self, triangle, p_critical: float = 0.5):
             triangle = triangle.set_backend("numpy")
         xp = triangle.get_array_module()
 
-        m1 = triangle.link_ratio
-
-        # Rank link ratios by development period, assigning a score of 1 for the lowest
-        m1_val = xp.apply_along_axis(
-            func1d=rankdata, axis=2, arr=m1.values, nan_policy="omit"
-        ) * (m1.values * 0 + 1)
-
-        # Remove the last element from each column, and then rank again
-        m2 = triangle[triangle.valuation < triangle.valuation_date].link_ratio
-        m2.values = xp.apply_along_axis(
-            func1d=rankdata, axis=2, arr=m2.values, nan_policy="omit"
-        ) * (m2.values * 0 + 1)
-
-        m1 = m2.copy()
-
-        # remove the first column from m1 since it is not used in the comparison to m2
-        m1.values = m1_val[..., : m2.shape[2], 1:]
-
-        # Apply Spearman Rank Correlation formula
-        # numerator is the one in formula G4 of Mack (1994), "Measuring the
-        # Variability of Chain Ladder Reserve Estimates"
-        numerator = ((m1 - m2) ** 2).sum("origin")
-
-        # remove last column because it was not part of the comparison with m2
-        numerator.values = numerator.values[..., :-1]
-        numerator.ddims = numerator.ddims[:-1]
+        link_ratio = triangle.link_ratio.values
+
+        # T_k compares each factor F_ik of column k with F_i,k-1 of the same
+        # origin. A complete triangle has I - k such pairs in column k. The
+        # pairs present are counted per triangle, so a missing factor or an
+        # origin without data lowers the count.
+        current = link_ratio[..., 1:]
+        previous = link_ratio[..., :-1]
+        paired = ~xp.isnan(current) & ~xp.isnan(previous)
+        n_pairs = paired.sum(axis=2)
+
+        # Rank both columns among the paired factors only, 1 for the lowest,
+        # and sum the squared rank differences (numerator of formula G4)
+        squared_diff = xp.zeros(n_pairs.shape)
+        if paired.size:
+            ranks = [
+                xp.apply_along_axis(
+                    func1d=rankdata,
+                    axis=2,
+                    arr=xp.where(paired, column, xp.nan),
+                    nan_policy="omit",
+                )
+                for column in (current, previous)
+            ]
+            ranks = [xp.where(paired, rank, xp.nan) for rank in ranks]
+            squared_diff = xp.nansum((ranks[0] - ranks[1]) ** 2, axis=2)
+
+        # Formula G4 for each column with at least two pairs, weighted by one
+        # less than its number of pairs, as Mack proposes
+        usable = n_pairs >= 2
+        weight = xp.where(usable, n_pairs - 1, 0)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            t_k = xp.where(
+                usable, 1 - 6 * squared_diff / (n_pairs**3 - n_pairs), xp.nan
+            )
 
-        # n_dev_periods is the number of development periods in the triangle
-        # ("I" in Mack (1994)).
-        n_dev_periods = len(triangle.development)
+        # Drop trailing columns that are unusable in every triangle, such as
+        # the last column of a complete triangle, which has one pair
+        has_weight = xp.any(weight > 0, axis=(0, 1))
+        n_columns = int(xp.flatnonzero(has_weight)[-1]) + 1 if has_weight.any() else 0
+        t_k, weight = t_k[..., :n_columns], weight[..., :n_columns]
 
         # k values are the column indexes for which we are calculating T_k
-        k = xp.array(range(2, 2 + numerator.shape[3]))
-
-        # denominator is the one in formula G4 of Mack (1994)
-        denominator = ((n_dev_periods - k) ** 3 - n_dev_periods + k)[None, None, None]
-
-        # complete formula G4, results in array of each T_k value
-        self.t = 1 - 6 * xp.nan_to_num(numerator.values) / denominator
-
-        # per Mack, weight is one less than the number of pairs for each T_k
-        weight = (n_dev_periods - k - 1)[None, None, None]
-
-        # Calculate big T, the weighted average of the T_k values
-        t_expectation = (
-            xp.sum(xp.nan_to_num(weight * self.t), axis=3) / xp.sum(weight, axis=3)
-        )[..., None]
+        k = xp.array(range(2, 2 + n_columns))
+
+        # Big T is the weighted average of the T_k values. Its variance is one
+        # over the total weight, which on a complete triangle is formula G6,
+        # 2 / ((I - 2)(I - 3)). A triangle without a usable column gets NaN
+        # for both.
+        total_weight = xp.sum(weight, axis=2)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            t_expectation = xp.where(
+                total_weight > 0,
+                xp.sum(xp.nan_to_num(weight * t_k), axis=2) / total_weight,
+                xp.nan,
+            )
+            t_variance = xp.where(total_weight > 0, 1 / total_weight, xp.nan)
 
         idx = triangle.index.set_index(triangle.key_labels).index
 
-        # variance is result of formula G6
-        self.t_variance = 2 / ((n_dev_periods - 2) * (n_dev_periods - 3))
-
-        # array of t values
-        self.t = pd.DataFrame(self.t[0, 0, ...], columns=k, index=["T_k"])
+        # array of t values for the first triangle
+        self.t = pd.DataFrame(t_k[0, 0][None], columns=k, index=["T_k"])
 
-        # array of weights
-        self.weights = pd.DataFrame(weight[0, 0, ...], columns=k, index=["I-k-1"])
+        # array of weights for the first triangle
+        self.weights = pd.DataFrame(weight[0, 0][None], columns=k, index=["I-k-1"])
 
         # final big T
         self.t_expectation = pd.DataFrame(
-            t_expectation[..., 0, 0], columns=triangle.columns, index=idx
+            t_expectation, columns=triangle.columns, index=idx
         )
 
+        self.t_variance = pd.DataFrame(t_variance, columns=triangle.columns, index=idx)
+
         # table of Spearman's rank coefficients Tk, can be used to verify consistency
         # with paper
         self.corr = pd.concat([self.t, self.weights])
@@ -197,7 +207,8 @@ def __init__(self, triangle, p_critical: float = 0.5):
             norm.ppf(0.5 + (1 - p_critical) / 2) * xp.sqrt(self.t_variance),
         )
 
-        # if T lies outside this range, we reject the null hypothesis
+        # if T lies outside this range, we reject the null hypothesis; a NaN T
+        # compares False, so a triangle without a usable column is not flagged
         self.t_critical = (self.t_expectation < self.confidence_interval[0]) | (
             self.t_expectation > self.confidence_interval[1]
         )

chainladder/core/tests/test_correlation.py

@@ -2,6 +2,7 @@
 
 from typing import TYPE_CHECKING
 
+import numpy as np
 import pytest
 
 from scipy.stats import norm
@@ -93,12 +94,110 @@ def test_dev_corr_matches_mack(raa: Triangle) -> None:
     assert dc.t.iloc[0].tolist() == pytest.approx(t_k)
     assert dc.weights.iloc[0].tolist() == [7, 6, 5, 4, 3, 2, 1]
     assert dc.t_expectation.iloc[0, 0] == pytest.approx(0.070, abs=5e-4)
-    assert dc.t_variance == 1 / 28
-    assert dc.confidence_interval[0] == pytest.approx(-0.127, abs=5e-4)
-    assert dc.confidence_interval[1] == pytest.approx(0.127, abs=5e-4)
+    assert dc.t_variance.iloc[0, 0] == 1 / 28
+    assert dc.confidence_interval[0].iloc[0, 0] == pytest.approx(-0.127, abs=5e-4)
+    assert dc.confidence_interval[1].iloc[0, 0] == pytest.approx(0.127, abs=5e-4)
     assert not dc.t_critical.iloc[0, 0]
 
 
+def test_dev_corr_complete_triangle_counts(genins: Triangle) -> None:
+    """
+    On a complete triangle the pair counts are Mack's I - k.
+
+    genins has ten origins and ten development periods, so the weights are
+    I - k - 1 = 7, ..., 1 and Var(T) = 2 / ((I - 2)(I - 3)) = 1/28, as in
+    formula G6 of Mack (1994).
+    """
+    dc = genins.development_correlation(p_critical=0.5)
+    assert dc.weights.iloc[0].tolist() == [7, 6, 5, 4, 3, 2, 1]
+    assert dc.t_variance.iloc[0, 0] == 1 / 28
+
+
+def test_dev_corr_no_data(raa: Triangle) -> None:
+    """
+    A triangle without link ratios gets NaN and is not flagged.
+
+    GH #1577: with no pairs to rank there is no T_k to average, so T and its
+    variance are NaN, and t_critical is False, as ValuationCorrelation does
+    for empty diagonals.
+    """
+    dc = (raa * np.nan).development_correlation(p_critical=0.5)
+    assert np.isnan(dc.t_expectation.iloc[0, 0])
+    assert np.isnan(dc.t_variance.iloc[0, 0])
+    assert not dc.t_critical.iloc[0, 0]
+
+
+def test_dev_corr_empty_member(auto: Triangle) -> None:
+    """
+    An empty triangle among several leaves the others unchanged.
+
+    GH #1577: the auto sample has no incurred losses for CommercialAuto. That
+    triangle gets NaN and is not flagged, and each of the other three gets
+    the result it gets on its own.
+    """
+    dc = auto.development_correlation(p_critical=0.5)
+    assert np.isnan(dc.t_expectation.loc["CommercialAuto", "incurred"])
+    assert not dc.t_critical.loc["CommercialAuto", "incurred"]
+    others = [
+        ("CommercialAuto", "paid"),
+        ("PersonalAuto", "incurred"),
+        ("PersonalAuto", "paid"),
+    ]
+    for lob, column in others:
+        alone = auto.loc[lob][column].development_correlation(p_critical=0.5)
+        expected = alone.t_expectation.iloc[0, 0]
+        assert dc.t_expectation.loc[lob, column] == pytest.approx(expected)
+        assert dc.t_variance.loc[lob, column] == alone.t_variance.iloc[0, 0]
+
+
+def test_dev_corr_counts_origins_with_data(raa: Triangle) -> None:
+    """
+    I is the number of origins with data, not of development columns.
+
+    GH #1577: the last six origins of raa, kept on the full ten-column grid,
+    give the same T and variance as the same data cut to a 6 x 6 triangle,
+    where Var(T) = 2 / ((6 - 2)(6 - 3)) = 1/6.
+    """
+    recent = raa[raa.origin >= "1985"]
+    square = recent[recent.development <= 72]
+    on_grid = recent.development_correlation(p_critical=0.5)
+    cut = square.development_correlation(p_critical=0.5)
+    expected = cut.t_expectation.iloc[0, 0]
+    assert on_grid.t_expectation.iloc[0, 0] == pytest.approx(expected)
+    assert on_grid.t_variance.iloc[0, 0] == pytest.approx(1 / 6)
+    assert cut.t_variance.iloc[0, 0] == pytest.approx(1 / 6)
+
+
+@pytest.mark.parametrize("size", [2, 3])
+def test_dev_corr_too_small(raa: Triangle, size: int) -> None:
+    """
+    Triangles with fewer than four origins give NaN instead of raising.
+
+    GH #1577: T_k needs at least two pairs, which a triangle first has at
+    I = 4. A 3 x 3 triangle raised AttributeError and a 2 x 2 triangle
+    IndexError.
+    """
+    small = raa[raa.origin >= str(1991 - size)]
+    small = small[small.development <= 12 * size]
+    dc = small.development_correlation(p_critical=0.5)
+    assert np.isnan(dc.t_expectation.iloc[0, 0])
+    assert not dc.t_critical.iloc[0, 0]
+
+
+def test_dev_corr_smallest_triangle(raa: Triangle) -> None:
+    """
+    A 4 x 4 triangle has one T_k, from two pairs, with Var(T) = 1.
+
+    Formula G6 gives 2 / ((4 - 2)(4 - 3)) = 1, and k = 2 is the only column
+    with two pairs.
+    """
+    small = raa[raa.origin >= "1987"]
+    small = small[small.development <= 48]
+    dc = small.development_correlation(p_critical=0.5)
+    assert dc.weights.iloc[0].tolist() == [1]
+    assert dc.t_variance.iloc[0, 0] == 1
+
+
 def test_validate_critical(raa: Triangle) -> None:
     with pytest.raises(ValueError):
         raa.valuation_correlation(p_critical=1.5, total=True)

conftest.py

@@ -49,6 +49,10 @@
         "data": "xyz",
         "runs": ["normal_run", "sparse_only_run"],
     },
+    "auto": {
+        "data": "auto",
+        "runs": ["normal_run", "sparse_only_run"],
+    },
 }