diff --git a/chainladder/core/correlation.py b/chainladder/core/correlation.py index 73384d8ea..fd1ec829a 100644 --- a/chainladder/core/correlation.py +++ b/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] ) diff --git a/chainladder/core/tests/test_correlation.py b/chainladder/core/tests/test_correlation.py index 3bb9d0b0b..098a01c3f 100644 --- a/chainladder/core/tests/test_correlation.py +++ b/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) diff --git a/conftest.py b/conftest.py index 7818dadaf..a6751be35 100644 --- a/conftest.py +++ b/conftest.py @@ -49,6 +49,10 @@ "data": "xyz", "runs": ["normal_run", "sparse_only_run"], }, + "auto": { + "data": "auto", + "runs": ["normal_run", "sparse_only_run"], + }, }