Barnett Zehnwirth documentation and formulation rework (#657)

* Changed basic formuation of PTF in development.ipynb to correspond with implementation. Added note about difference in formulation. Fixed valuation (iota) summation term.

* Added documentation for utils.utility_functions.PTF_formula with example from 2008 paper. Updated reference list

* changes how development and valuation bounds are given to PTF_formula. More changes to BZ documentation.

* reworked pathological BZ example (no valuation coefficients). Added graphs from 2008 paper for 'reasonable model'

* PTF_formula now takes lists. Alphas are the start of every origin bucket. Gammas/iotas are the endpoints of each linear piece. Updated test_barnzehn accordingly

* updated docs to reflect changes to PTF_formula

* More checks for how BZ parameters are passed. Moved PTF_formula to within BZ estimator (called during fit). Changed how gamma is passed: index from 0 and increment grain-agnostic. Docs not updated yet

* updated user_guide/development to reflect PTF_formula integration into BZ

* moved PTF_formula into PatsyFormula

* fix typo

* Revert "fix typo"

This reverts commit df7d310.

* Revert "moved PTF_formula into PatsyFormula"

This reverts commit 840bc4f.

* fix typo and clarify formulation when using patsy directly

---------

Co-authored-by: danielfong-act <danielfong247@gmial.com>

Author: danielfong-act

chainladder/development/barnzehn.py

@@ -9,10 +9,7 @@
 from chainladder.development.glm import TweedieGLM
 from sklearn.linear_model import LinearRegression
 from sklearn.pipeline import Pipeline
-import warnings
-from chainladder.utils.utility_functions import PatsyFormula
-from patsy import ModelDesc
-
+from chainladder.utils.utility_functions import PatsyFormula, PTF_formula
 
 class BarnettZehnwirth(TweedieGLM):
     """ This estimator enables modeling from the Probabilistic Trend Family as
@@ -31,21 +28,38 @@ class BarnettZehnwirth(TweedieGLM):
     response:  str
         Column name for the reponse variable of the GLM.  If ommitted, then the
         first column of the Triangle will be used.
-
+    alpha: list of int
+        List of origin periods denoting the first indices of each group 
+    gamma: list of int
+    iota: list of int
 
     """
 
-    def __init__(self, drop=None,drop_valuation=None,formula='C(origin) + development', response=None):
+    def __init__(self, drop=None,drop_valuation=None,formula=None, response=None, alpha=None, gamma=None, iota=None):
         self.drop = drop
         self.drop_valuation = drop_valuation
-        self.formula = formula
-        self.response = response
 
+        self.response = response
+        if formula and (alpha or gamma or iota):
+            raise ValueError("Model can only be specified by either a formula or some combination of alpha, gamma and iota.")
+        if not (formula or alpha or gamma or iota):
+            raise ValueError("Model must be specified, either a formula or some combination of alpha, gamma and iota.")
+        for Greek in [alpha,gamma,iota]:
+            if Greek:
+                if not ( (type(Greek) is list) and all(type(bound) is int for bound in Greek) ):
+                    raise ValueError("Alpha, gamma and iota must be given as lists of integers, specifying periods.")
+        self.formula = formula
+        self.alpha = alpha
+        self.gamma = gamma
+        self.iota = iota
+        
     def fit(self, X, y=None, sample_weight=None):
         if max(X.shape[:2]) > 1:
             raise ValueError("Only single index/column triangles are supported")
         tri = X.cum_to_incr().log()
         response = X.columns[0] if not self.response else self.response
+        if(not self.formula):
+            self.formula = PTF_formula(self.alpha,self.gamma,self.iota,dgrain=min(tri.development))
         self.model_ = DevelopmentML(Pipeline(steps=[
             ('design_matrix', PatsyFormula(self.formula)),
             ('model', LinearRegression(fit_intercept=False))]),

chainladder/development/tests/test_barnzehn.py

@@ -1,7 +1,6 @@
 import numpy as np
 import chainladder as cl
 import pytest
-from chainladder.utils.utility_functions import PTF_formula
 abc = cl.load_sample('abc')
 
 def test_basic_bz():
@@ -40,13 +39,11 @@ def test_bz_2008():
     exposure=np.array([[2.2], [2.4], [2.2], [2.0], [1.9], [1.6], [1.6], [1.8], [2.2], [2.5], [2.6]])
     abc_adj = abc/exposure
 
-    origin_buckets = [(0,1),(2,2),(3,4),(5,10)]
-    dev_buckets = [(24,36),(36,48),(48,84),(84,108),(108,144)]
-    val_buckets = [(1,8),(8,9),(9,12)]
-    
-    abc_formula = PTF_formula(abc_adj,alpha=origin_buckets,gamma=dev_buckets,iota=val_buckets)
-    
-    model=cl.BarnettZehnwirth(formula=abc_formula, drop=('1982',72)).fit(abc_adj)
+    origin_buckets = [0,2,3,5]
+    dev_buckets = [0,1,2,5,7,10]
+    val_buckets = [0,7,8,11]
+
+    model=cl.BarnettZehnwirth(drop=('1982',72),alpha=origin_buckets,gamma=dev_buckets,iota=val_buckets).fit(abc_adj)
     assert np.all(
         np.around(model.coef_.values,4).flatten()
         == np.array([11.1579,0.1989,0.0703,0.0919,0.1871,-0.3771,-0.4465,-0.3727,-0.3154,0.0432,0.0858,0.1464])

chainladder/utils/utility_functions.py

@@ -769,27 +769,30 @@ def model_diagnostics(model, name=None, groupby=None):
     return concat(triangles, 0)
 
 
-def PTF_formula(tri: Triangle, alpha: ArrayLike = None, gamma: ArrayLike = None, iota: ArrayLike = None):
+def PTF_formula(alpha: list = None, gamma: list = None, iota: list = None,dgrain: int = 12):
     """ Helper formula that builds a patsy formula string for the BarnettZehnwirth 
     estimator.  Each axis's parameters can be grouped together. Groups of origin 
-    parameters (alpha) are set equal, and are specified by a ranges (inclusive). 
+    parameters (alpha) are set equal, and are specified by the first period in each bin. 
     Groups of development (gamma) and valuation (iota) parameters are fit to 
-    separate linear trends, specified as tuples denoting ranges with shared endpoints.
+    separate linear trends, specified a list denoting the endpoints of the linear pieces.
     In other words, development and valuation trends are fit to a piecewise linear model.
     A triangle must be supplied to provide some critical information.
     """
     formula_parts=[]
     if(alpha):
         # The intercept term takes the place of the first alpha
         for ind,a in enumerate(alpha):
-            if(a[0]==0):
+            if(a==0):
                 alpha=alpha[:ind]+alpha[(ind+1):]
-        formula_parts += ['+'.join([f'I({x[0]} <= origin)' for x in alpha])]
-    if(gamma):
-        dgrain = min(tri.development)
-        formula_parts += ['+'.join([f'I((np.minimum({x[1]-dgrain},development) - np.minimum({x[0]-dgrain},development))/{dgrain})' for x in gamma])]
+        formula_parts += ['+'.join([f'I({x} <= origin)' for x in alpha])]
+    if(gamma): 
+        # preprocess gamma to align with grain
+        graingamma = [(i+1)*dgrain for i in gamma]
+        for ind in range(1,len(graingamma)):
+            formula_parts += ['+'.join([f'I((np.minimum({graingamma[ind]},development) - np.minimum({graingamma[ind-1]},development))/{dgrain})'])]
     if(iota):
-        formula_parts += ['+'.join([f'I(np.minimum({x[1]-1},valuation) - np.minimum({x[0]-1},valuation))' for x in iota])]
+        for ind in range(1,len(iota)):
+            formula_parts += ['+'.join([f'I(np.minimum({iota[ind]},valuation) - np.minimum({iota[ind-1]},valuation))'])]
     if(formula_parts):
         return '+'.join(formula_parts)
     return ''

docs/images/plot_ptf_coefficients.png

@@ -87,3 +87,11 @@ @article{shapland2016
 	year = {2016},
 	url = {https://live-casact.pantheonsite.io/sites/default/files/2021-02/04-shapland.pdf}
 }
+
+@article{barnett2008,
+    author = {Barnett, G. and Zehnwirth, B.},
+	title = {Modeling with the {M}ultivariate {P}robabilistic {T}rend {F}amily},
+	journal = {Casualty Actuarial Society E-Forum},
+	year = {2008},
+	url = {https://www.casact.org/sites/default/files/database/forum_08fforum_3barnett_zehnwirth.pdf}
+}