From 0f079490d5d7895c748fc2b05af93e0ac1b54739 Mon Sep 17 00:00:00 2001
From: Ethan Kang <ethan.yskang@gmail.com>
Date: Sat, 16 May 2026 01:16:21 -0700
Subject: [PATCH 1/5] docs: add development estimator doctest examples

---
 chainladder/development/barnzehn.py    | 29 ++++++++++
 chainladder/development/glm.py         | 42 +++++++++++++++
 chainladder/development/incremental.py | 41 ++++++++++++++
 chainladder/development/learning.py    | 75 ++++++++++++++++++++++++++
 chainladder/development/outstanding.py | 30 +++++++++++
 chainladder/methods/mack.py            | 31 ++++++++---
 6 files changed, 242 insertions(+), 6 deletions(-)

diff --git a/chainladder/development/barnzehn.py b/chainladder/development/barnzehn.py
index 20ae0b79..d4253ab3 100644
--- a/chainladder/development/barnzehn.py
+++ b/chainladder/development/barnzehn.py
@@ -33,6 +33,35 @@ class BarnettZehnwirth(TweedieGLM):
     gamma: list of int
     iota: list of int
 
+    Examples
+    --------
+    A patsy ``formula`` and the built-in ``alpha`` / ``gamma`` / ``iota`` PTF
+    specification can both fit the same triangle; the leading fitted
+    coefficient differs because the design matrices differ.
+
+    .. testsetup::
+
+        import chainladder as cl
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("abc")
+        m_formula = cl.BarnettZehnwirth(
+            formula="C(origin)+C(development)"
+        ).fit(tri)
+        m_ptf = cl.BarnettZehnwirth(
+            alpha=[0, 5], gamma=[0, 2, 5], iota=[0, 7, 11]
+        ).fit(tri)
+        print(float(np.round(m_formula.coef_.values.flatten()[0], 3)))
+        print(float(np.round(m_ptf.coef_.values.flatten()[0], 3)))
+
+    .. testoutput::
+
+        11.837
+        12.151
+
     """
 
     def __init__(self, drop=None,drop_valuation=None,formula=None, response=None, alpha=None, gamma=None, iota=None):
diff --git a/chainladder/development/glm.py b/chainladder/development/glm.py
index c44bb0d5..a3d50f56 100644
--- a/chainladder/development/glm.py
+++ b/chainladder/development/glm.py
@@ -76,6 +76,48 @@ class TweedieGLM(DevelopmentBase):
     ----------
     model_: sklearn.Pipeline
         A scikit-learn Pipeline of the GLM
+
+    Examples
+    --------
+    ``design_matrix`` controls which patsy terms enter the GLM; dropping
+    ``C(origin)`` changes the first fitted LDF.
+
+    .. testsetup::
+
+        import chainladder as cl
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("genins")
+        m_full = cl.TweedieGLM(
+            power=1, design_matrix="C(development) + C(origin)"
+        ).fit(tri)
+        m_dev = cl.TweedieGLM(power=1, design_matrix="C(development)").fit(tri)
+        print(float(np.round(m_full.ldf_.values[0, 0, 0, 0], 4)))
+        print(float(np.round(m_dev.ldf_.values[0, 0, 0, 0], 4)))
+
+    .. testoutput::
+
+        3.491
+        3.5085
+
+    ``power`` and ``link`` select the Tweedie family; a Normal GLM
+    (``power=0`` with ``link='identity'``) yields a different pattern.
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("genins")
+        m = cl.TweedieGLM(power=0, link="identity").fit(tri)
+        print(float(np.round(m.ldf_.values[0, 0, 0, 0], 2)))
+
+    .. testoutput::
+
+        2.31
+
     """
 
     def __init__(self, design_matrix='C(development) + C(origin)',
diff --git a/chainladder/development/incremental.py b/chainladder/development/incremental.py
index fb0537e0..82266857 100644
--- a/chainladder/development/incremental.py
+++ b/chainladder/development/incremental.py
@@ -69,6 +69,47 @@ class IncrementalAdditive(DevelopmentBase):
         A triangle of full incremental values.
 
 
+    Examples
+    --------
+    Basic fit on ``ia_sample`` with exposure on the latest diagonal.
+
+    .. testsetup::
+
+        import chainladder as cl
+
+    .. testcode::
+
+        tri = cl.load_sample("ia_sample")
+        model = cl.IncrementalAdditive().fit(
+            tri["loss"], sample_weight=tri["exposure"].latest_diagonal
+        )
+        print(model.ldf_.shape)
+
+    .. testoutput::
+
+        (1, 1, 6, 5)
+
+    ``future_trend`` (when non-zero) changes extrapolated incrementals in the
+    lower triangle even when ``trend`` is held at zero; here the summed
+    fitted incrementals increase.
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("ia_sample")
+        loss = tri["loss"]
+        sw = tri["exposure"].latest_diagonal
+        m0 = cl.IncrementalAdditive(trend=0, future_trend=0).fit(loss, sample_weight=sw)
+        m1 = cl.IncrementalAdditive(trend=0, future_trend=0.1).fit(loss, sample_weight=sw)
+        print(float(np.round(np.nansum(m0.incremental_.values), 1)))
+        print(float(np.round(np.nansum(m1.incremental_.values), 1)))
+
+    .. testoutput::
+
+        30988.1
+        33360.1
+
     """
 
     def __init__(
diff --git a/chainladder/development/learning.py b/chainladder/development/learning.py
index 041a9abc..c07a12a4 100644
--- a/chainladder/development/learning.py
+++ b/chainladder/development/learning.py
@@ -50,6 +50,81 @@ class DevelopmentML(DevelopmentBase):
         The estimated loss development patterns.
     cdf_: Triangle
         The estimated cumulative development patterns.
+
+    Examples
+    --------
+    ``fit_incrementals`` toggles whether the pipeline fits on incrementals
+    versus cumulatives, which shifts the implied ``ldf_``.
+
+    .. testsetup::
+
+        import chainladder as cl
+
+    .. testcode::
+
+        import numpy as np
+        from sklearn.linear_model import LinearRegression
+        from sklearn.pipeline import Pipeline
+
+        from chainladder.utils.utility_functions import PatsyFormula
+
+        tri = cl.load_sample("genins")
+        pipe = Pipeline(
+            steps=[
+                ("design_matrix", PatsyFormula("C(development)")),
+                ("model", LinearRegression(fit_intercept=False)),
+            ]
+        )
+        m_incr = cl.DevelopmentML(
+            pipe, y_ml=[tri.columns[0]], fit_incrementals=True
+        ).fit(tri)
+        m_cum = cl.DevelopmentML(
+            pipe, y_ml=[tri.columns[0]], fit_incrementals=False
+        ).fit(tri)
+        print(float(np.round(m_incr.ldf_.values[0, 0, 0, 0], 4)))
+        print(float(np.round(m_cum.ldf_.values[0, 0, 0, 0], 4)))
+
+    .. testoutput::
+
+        3.508
+        3.515
+
+    With ``weighted_step='model'``, ``sample_weight`` is forwarded into the
+    final regressor; squaring the triangle as a crude weight changes the first
+    LDF versus an unweighted fit.
+
+    .. testcode::
+
+        import numpy as np
+        from sklearn.linear_model import LinearRegression
+        from sklearn.pipeline import Pipeline
+
+        from chainladder.utils.utility_functions import PatsyFormula
+
+        tri = cl.load_sample("genins")
+        pipe = Pipeline(
+            steps=[
+                ("design_matrix", PatsyFormula("C(development)")),
+                ("model", LinearRegression(fit_intercept=False)),
+            ]
+        )
+        m0 = cl.DevelopmentML(
+            pipe, y_ml=[tri.columns[0]], fit_incrementals=False
+        ).fit(tri)
+        m1 = cl.DevelopmentML(
+            pipe,
+            y_ml=[tri.columns[0]],
+            fit_incrementals=False,
+            weighted_step="model",
+        ).fit(tri, sample_weight=tri * tri)
+        print(float(np.round(m0.ldf_.values[0, 0, 0, 0], 4)))
+        print(float(np.round(m1.ldf_.values[0, 0, 0, 0], 4)))
+
+    .. testoutput::
+
+        3.515
+        3.4459
+
     """
 
     def __init__(self, estimator_ml=None, y_ml=None, autoregressive=False,
diff --git a/chainladder/development/outstanding.py b/chainladder/development/outstanding.py
index ffaa5528..bc6ddbc5 100644
--- a/chainladder/development/outstanding.py
+++ b/chainladder/development/outstanding.py
@@ -49,6 +49,36 @@ class CaseOutstanding(DevelopmentBase):
         The paid to prior case ratios used for fitting the estimator
     paid_ldf_:
         The selected paid to prior case ratios of the fitted estimator
+
+    Examples
+    --------
+    ``paid_n_periods`` and ``case_n_periods`` control how many recent origin
+    years inform the ``Development`` weights that smooth the paid and case
+    patterns.
+
+    .. testsetup::
+
+        import chainladder as cl
+
+    .. testcode::
+
+        tri = cl.load_sample("usauto")
+        all_years = cl.CaseOutstanding(
+            paid_to_incurred=("paid", "incurred")
+        ).fit(tri)
+        three = cl.CaseOutstanding(
+            paid_to_incurred=("paid", "incurred"),
+            paid_n_periods=3,
+            case_n_periods=3,
+        ).fit(tri)
+        print(round(float(all_years.paid_ldf_.values[0, 0, 0, 0]), 6))
+        print(round(float(three.paid_ldf_.values[0, 0, 0, 0]), 6))
+
+    .. testoutput::
+
+        0.842814
+        0.833138
+
     """
 
     def __init__(
diff --git a/chainladder/methods/mack.py b/chainladder/methods/mack.py
index cfcf71e0..d42d6306 100644
--- a/chainladder/methods/mack.py
+++ b/chainladder/methods/mack.py
@@ -46,17 +46,17 @@ class MackChainladder(Chainladder):
     which combines the deterministic chainladder estimate with Mack's
     stochastic standard error.
 
-    .. testsetup:
+    .. testsetup::
 
         import chainladder as cl
 
-    .. testcode:
+    .. testcode::
 
         tr = cl.load_sample('ukmotor')
         model = cl.MackChainladder().fit(tr)
         print(model.summary_)
 
-    .. testoutput:
+    .. testoutput::
 
                Latest          IBNR      Ultimate  Mack Std Err
         2007  12690.0           NaN  12690.000000           NaN
@@ -71,14 +71,30 @@ class MackChainladder(Chainladder):
     :class:`Chainladder`. Mack's contribution is the stochastic standard error
     in the rightmost column, which can be aggregated across origins.
 
-    .. testcode:
+    .. testcode::
 
         print(model.total_mack_std_err_)
 
-    .. testoutput:
+    .. testoutput::
 
         columns        values
         (Total,)  1424.531543
+
+    Mack's total error depends on how ``ldf_`` and ``sigma_`` were produced.
+    Here the same triangle is pre-smoothed with :class:`Development` using
+    ``average='simple'`` instead of the default volume weights before fitting
+    ``MackChainladder``, which raises the aggregate Mack standard error.
+
+    .. testcode::
+
+        tr = cl.load_sample("ukmotor")
+        tr_simple = cl.Development(average="simple").fit_transform(tr)
+        print(cl.MackChainladder().fit(tr_simple).total_mack_std_err_)
+
+    .. testoutput::
+
+        columns        values
+        (Total,)  1591.603339
     """
 
     def fit(self, X, y=None, sample_weight=None):
@@ -217,7 +233,7 @@ def full_std_err_(self):
             model = cl.MackChainladder().fit(tr)
             print(model.full_std_err_)
 
-        .. testoutput
+        .. testoutput::
 
                         12        24        36        48        60        72   84
             2007  0.047826  0.040745  0.031412  0.010337  0.001431  0.001523  0.0
@@ -268,6 +284,7 @@ def total_process_risk_(self):
         --------
 
         .. testsetup::
+
             import chainladder as cl
 
         .. testcode::
@@ -340,9 +357,11 @@ def mack_std_err_(self):
         error per origin.
 
         .. testsetup::
+
             import chainladder as cl
 
         .. testcode::
+
             tr = cl.load_sample('ukmotor')
             model = cl.MackChainladder().fit(tr)
             print(model.mack_std_err_.iloc[..., -3:, -3:])

From 4a1b11c35bd100a1de63cf62a46a811165aab8c7 Mon Sep 17 00:00:00 2001
From: Ethan Kang <ethan.yskang@gmail.com>
Date: Sat, 16 May 2026 15:54:37 -0700
Subject: [PATCH 2/5] docs: clarify development estimator examples

---
 chainladder/development/barnzehn.py    |  7 ++++---
 chainladder/development/incremental.py | 10 ++++++----
 chainladder/development/learning.py    | 14 ++++++++------
 3 files changed, 18 insertions(+), 13 deletions(-)

diff --git a/chainladder/development/barnzehn.py b/chainladder/development/barnzehn.py
index d4253ab3..e07317fa 100644
--- a/chainladder/development/barnzehn.py
+++ b/chainladder/development/barnzehn.py
@@ -35,9 +35,10 @@ class BarnettZehnwirth(TweedieGLM):
 
     Examples
     --------
-    A patsy ``formula`` and the built-in ``alpha`` / ``gamma`` / ``iota`` PTF
-    specification can both fit the same triangle; the leading fitted
-    coefficient differs because the design matrices differ.
+    Use an explicit patsy ``formula`` when you want direct control over the
+    regression design, or use the PTF ``alpha`` / ``gamma`` / ``iota`` pieces
+    when the model should follow Barnett-Zehnwirth trend groupings. The two
+    specifications fit the same triangle but produce different design matrices.
 
     .. testsetup::
 
diff --git a/chainladder/development/incremental.py b/chainladder/development/incremental.py
index 82266857..7dcf4421 100644
--- a/chainladder/development/incremental.py
+++ b/chainladder/development/incremental.py
@@ -71,7 +71,9 @@ class IncrementalAdditive(DevelopmentBase):
 
     Examples
     --------
-    Basic fit on ``ia_sample`` with exposure on the latest diagonal.
+    Use ``IncrementalAdditive`` when incremental losses should be modeled
+    against an exposure base rather than developed from age-to-age ratios. The
+    ``ia_sample`` data include both loss and exposure columns for this workflow.
 
     .. testsetup::
 
@@ -89,9 +91,9 @@ class IncrementalAdditive(DevelopmentBase):
 
         (1, 1, 6, 5)
 
-    ``future_trend`` (when non-zero) changes extrapolated incrementals in the
-    lower triangle even when ``trend`` is held at zero; here the summed
-    fitted incrementals increase.
+    Apply ``future_trend`` when the completed lower triangle should reflect a
+    prospective trend assumption. Holding historical ``trend`` at zero isolates
+    the effect on the fitted future incrementals.
 
     .. testcode::
 
diff --git a/chainladder/development/learning.py b/chainladder/development/learning.py
index c07a12a4..812075d2 100644
--- a/chainladder/development/learning.py
+++ b/chainladder/development/learning.py
@@ -53,8 +53,10 @@ class DevelopmentML(DevelopmentBase):
 
     Examples
     --------
-    ``fit_incrementals`` toggles whether the pipeline fits on incrementals
-    versus cumulatives, which shifts the implied ``ldf_``.
+    Choose the response scale before fitting a machine-learning development
+    model. Setting ``fit_incrementals=True`` models incremental dollar amounts;
+    setting it to ``False`` models cumulative values before translating the
+    result back into development patterns.
 
     .. testsetup::
 
@@ -89,9 +91,9 @@ class DevelopmentML(DevelopmentBase):
         3.508
         3.515
 
-    With ``weighted_step='model'``, ``sample_weight`` is forwarded into the
-    final regressor; squaring the triangle as a crude weight changes the first
-    LDF versus an unweighted fit.
+    Pass sample weights into the scikit-learn step when observations should not
+    contribute equally to the regression. ``weighted_step='model'`` forwards the
+    weights into the final estimator in the pipeline.
 
     .. testcode::
 
@@ -313,4 +315,4 @@ def transform(self, X):
         X_new.ldf_.valuation_date = pd.to_datetime(options.ULT_VAL)
         X_new._set_slicers()
         X_new.predicted_data_ = predicted_data
-        return X_new
\ No newline at end of file
+        return X_new

From c547ad31f1ba22f083cd09271af3abed5e1ebc9a Mon Sep 17 00:00:00 2001
From: Ethan Kang <ethan.yskang@gmail.com>
Date: Sat, 16 May 2026 16:09:58 -0700
Subject: [PATCH 3/5] docs: expand development example takeaways

---
 chainladder/development/barnzehn.py    | 10 ++++++----
 chainladder/development/incremental.py |  8 +++++---
 chainladder/development/learning.py    |  8 +++++---
 3 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/chainladder/development/barnzehn.py b/chainladder/development/barnzehn.py
index e07317fa..307b6aac 100644
--- a/chainladder/development/barnzehn.py
+++ b/chainladder/development/barnzehn.py
@@ -35,10 +35,12 @@ class BarnettZehnwirth(TweedieGLM):
 
     Examples
     --------
-    Use an explicit patsy ``formula`` when you want direct control over the
-    regression design, or use the PTF ``alpha`` / ``gamma`` / ``iota`` pieces
-    when the model should follow Barnett-Zehnwirth trend groupings. The two
-    specifications fit the same triangle but produce different design matrices.
+    Choose how to describe the regression structure before fitting the model.
+    A patsy ``formula`` is useful when you want direct control over origin and
+    development factors. The PTF ``alpha`` / ``gamma`` / ``iota`` arguments are
+    a shorthand for Barnett-Zehnwirth period, trend, and final-period groupings.
+    Both approaches fit the same triangle here; the different coefficients show
+    that the selected structure changes the design matrix being estimated.
 
     .. testsetup::
 
diff --git a/chainladder/development/incremental.py b/chainladder/development/incremental.py
index 7dcf4421..c3bcf52a 100644
--- a/chainladder/development/incremental.py
+++ b/chainladder/development/incremental.py
@@ -71,9 +71,11 @@ class IncrementalAdditive(DevelopmentBase):
 
     Examples
     --------
-    Use ``IncrementalAdditive`` when incremental losses should be modeled
-    against an exposure base rather than developed from age-to-age ratios. The
-    ``ia_sample`` data include both loss and exposure columns for this workflow.
+    Use ``IncrementalAdditive`` when the reserving assumption is driven by
+    incremental loss per unit of exposure, rather than by observed age-to-age
+    ratios alone. The example passes the latest exposure diagonal as
+    ``sample_weight``; the fitted object still returns development patterns, so
+    it can be used in the same workflow as other development estimators.
 
     .. testsetup::
 
diff --git a/chainladder/development/learning.py b/chainladder/development/learning.py
index 812075d2..ab603d98 100644
--- a/chainladder/development/learning.py
+++ b/chainladder/development/learning.py
@@ -54,9 +54,11 @@ class DevelopmentML(DevelopmentBase):
     Examples
     --------
     Choose the response scale before fitting a machine-learning development
-    model. Setting ``fit_incrementals=True`` models incremental dollar amounts;
-    setting it to ``False`` models cumulative values before translating the
-    result back into development patterns.
+    model. In this class, ``fit_incrementals=True`` means the scikit-learn model
+    is trained on actual incremental dollar amounts. Setting it to ``False``
+    trains on cumulative values instead. Both fitted models are translated back
+    into ``ldf_``, so the comparison is about the training target, not a request
+    to model age-to-age factors directly.
 
     .. testsetup::
 

From 5b0f7c47a0bf6fa801dc7f7cddd40f578f2aff22 Mon Sep 17 00:00:00 2001
From: Ethan Kang <ekang78@gatech.edu>
Date: Sun, 17 May 2026 23:45:41 -0700
Subject: [PATCH 4/5] docs: improve development estimator and Mack docstring
 examples

Address @henrydingliu feedback: reframe all examples as problem-solution
narratives, fix inaccurate fit_incrementals description in DevelopmentML,
add autoregressive example, expand test outputs beyond single LDF values,
and explain weighted_step Pipeline routing mechanism.
---
 chainladder/development/barnzehn.py    | 39 ++++++++---
 chainladder/development/glm.py         | 44 ++++++++----
 chainladder/development/incremental.py | 31 ++++++---
 chainladder/development/learning.py    | 95 +++++++++++++++++++++-----
 chainladder/development/outstanding.py | 29 ++++++--
 chainladder/methods/mack.py            | 60 +++++++++++-----
 6 files changed, 222 insertions(+), 76 deletions(-)

diff --git a/chainladder/development/barnzehn.py b/chainladder/development/barnzehn.py
index 307b6aac..9e162d6a 100644
--- a/chainladder/development/barnzehn.py
+++ b/chainladder/development/barnzehn.py
@@ -35,12 +35,31 @@ class BarnettZehnwirth(TweedieGLM):
 
     Examples
     --------
-    Choose how to describe the regression structure before fitting the model.
-    A patsy ``formula`` is useful when you want direct control over origin and
-    development factors. The PTF ``alpha`` / ``gamma`` / ``iota`` arguments are
-    a shorthand for Barnett-Zehnwirth period, trend, and final-period groupings.
-    Both approaches fit the same triangle here; the different coefficients show
-    that the selected structure changes the design matrix being estimated.
+    Standard ``Development`` assumes loss development ratios are stable across
+    accident years. When a triangle shows a calendar-year effect (e.g., recent
+    diagonals are systematically heavier or lighter across all accident years
+    due to inflation or a shift in case reserving practices), that assumption
+    breaks down. ``BarnettZehnwirth`` addresses this by fitting a log-linear
+    model that decomposes the triangle into separate origin, development, and
+    calendar-year trend components.
+
+    Two interfaces let you specify the trend structure. The ``formula``
+    argument takes any patsy expression directly. The Probabilistic Trend
+    Family (PTF) arguments provide a structured shorthand: ``alpha`` groups
+    accident years that share the same level effect, ``gamma`` defines
+    breakpoints for a piecewise linear development-age trend, and ``iota``
+    defines breakpoints for a piecewise linear calendar-year (diagonal) trend.
+
+    The ``abc`` triangle has 11 accident years (1977-1987) and 11 development
+    ages (12-132 months). Suppose an actuary notices from the triangle (or
+    from external information such as a change in inflation or legal
+    environment) that accident years before 1982 behave differently from those
+    after, that development speed changes at the 36-month and 72-month marks,
+    and that calendar-year trends shift at two points in the diagonal sequence.
+    Those observations translate directly into ``alpha=[0, 5]``,
+    ``gamma=[0, 2, 5]``, and ``iota=[0, 7, 11]``. The first three fitted
+    coefficients differ from the unconstrained ``formula`` model, reflecting
+    the additional structure the actuary has imposed.
 
     .. testsetup::
 
@@ -57,13 +76,13 @@ class BarnettZehnwirth(TweedieGLM):
         m_ptf = cl.BarnettZehnwirth(
             alpha=[0, 5], gamma=[0, 2, 5], iota=[0, 7, 11]
         ).fit(tri)
-        print(float(np.round(m_formula.coef_.values.flatten()[0], 3)))
-        print(float(np.round(m_ptf.coef_.values.flatten()[0], 3)))
+        print(np.round(m_formula.coef_.values.flatten()[:3], 3))
+        print(np.round(m_ptf.coef_.values.flatten()[:3], 3))
 
     .. testoutput::
 
-        11.837
-        12.151
+        [11.837  0.179  0.345]
+        [12.151  0.274 -0.064]
 
     """
 
diff --git a/chainladder/development/glm.py b/chainladder/development/glm.py
index a3d50f56..20c9ca48 100644
--- a/chainladder/development/glm.py
+++ b/chainladder/development/glm.py
@@ -79,8 +79,16 @@ class TweedieGLM(DevelopmentBase):
 
     Examples
     --------
-    ``design_matrix`` controls which patsy terms enter the GLM; dropping
-    ``C(origin)`` changes the first fitted LDF.
+    ``TweedieGLM`` with ``power=1`` and log link implements the
+    Over-Dispersed Poisson (ODP) model, the standard GLM equivalent of
+    volume-weighted chainladder. Its main advantage over column-by-column
+    development is that it fits the entire triangle simultaneously, which
+    allows parameter reduction: replacing categorical dummies (one per
+    accident year or development age) with continuous trend terms lowers the
+    parameter count while often staying close to the traditional result. The
+    default categorical design matrix replicates chainladder; switching to
+    continuous ``development + origin`` terms is a more parsimonious choice
+    when data are sparse.
 
     .. testsetup::
 
@@ -91,32 +99,42 @@ class TweedieGLM(DevelopmentBase):
         import numpy as np
 
         tri = cl.load_sample("genins")
-        m_full = cl.TweedieGLM(
+        m_cat = cl.TweedieGLM(
             power=1, design_matrix="C(development) + C(origin)"
         ).fit(tri)
-        m_dev = cl.TweedieGLM(power=1, design_matrix="C(development)").fit(tri)
-        print(float(np.round(m_full.ldf_.values[0, 0, 0, 0], 4)))
-        print(float(np.round(m_dev.ldf_.values[0, 0, 0, 0], 4)))
+        m_cont = cl.TweedieGLM(
+            power=1, design_matrix="development + origin"
+        ).fit(tri)
+        print(m_cat.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(m_cont.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
 
     .. testoutput::
 
-        3.491
-        3.5085
+        [3.491  1.7474 1.4574 1.1739]
+        [1.9277 1.4465 1.2864 1.2065]
 
-    ``power`` and ``link`` select the Tweedie family; a Normal GLM
-    (``power=0`` with ``link='identity'``) yields a different pattern.
+    The ``power`` parameter selects the variance-mean relationship. For lines
+    where some development periods show no incremental activity, the compound
+    Poisson-Gamma distribution (``1 < power < 2``) is more appropriate than
+    pure Poisson (``power=1``). The ``alpha`` parameter (default 1.0) controls
+    L2 regularization strength: increasing it smooths fitted coefficients on
+    sparse or over-parameterised triangles. If the optimizer raises a
+    ``ConvergenceWarning``, increase ``max_iter`` (default 100).
 
     .. testcode::
 
         import numpy as np
 
         tri = cl.load_sample("genins")
-        m = cl.TweedieGLM(power=0, link="identity").fit(tri)
-        print(float(np.round(m.ldf_.values[0, 0, 0, 0], 2)))
+        m_poisson = cl.TweedieGLM(power=1, link="log").fit(tri)
+        m_tweedie = cl.TweedieGLM(power=1.5, link="log").fit(tri)
+        print(m_poisson.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(m_tweedie.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
 
     .. testoutput::
 
-        2.31
+        [3.491  1.7474 1.4574 1.1739]
+        [3.7032 1.7569 1.4591 1.1776]
 
     """
 
diff --git a/chainladder/development/incremental.py b/chainladder/development/incremental.py
index c3bcf52a..31605c2f 100644
--- a/chainladder/development/incremental.py
+++ b/chainladder/development/incremental.py
@@ -71,11 +71,15 @@ class IncrementalAdditive(DevelopmentBase):
 
     Examples
     --------
-    Use ``IncrementalAdditive`` when the reserving assumption is driven by
-    incremental loss per unit of exposure, rather than by observed age-to-age
-    ratios alone. The example passes the latest exposure diagonal as
-    ``sample_weight``; the fitted object still returns development patterns, so
-    it can be used in the same workflow as other development estimators.
+    Use ``IncrementalAdditive`` when incremental losses per unit of exposure
+    are expected to be more stable across accident years than age-to-age
+    factors. This is common in lines where claim payments are driven by
+    exposure volume: for example, a workers' compensation book where
+    claim costs per payroll dollar are steady even as the cumulative reported
+    losses grow at different rates by accident year. Passing the latest
+    exposure diagonal as ``sample_weight`` normalises each period's
+    incremental by exposure before averaging. The fitted ``ldf_`` slots
+    directly into ``Chainladder`` or any other method.
 
     .. testsetup::
 
@@ -83,19 +87,26 @@ class IncrementalAdditive(DevelopmentBase):
 
     .. testcode::
 
+        import numpy as np
+
         tri = cl.load_sample("ia_sample")
         model = cl.IncrementalAdditive().fit(
             tri["loss"], sample_weight=tri["exposure"].latest_diagonal
         )
-        print(model.ldf_.shape)
+        print(model.ldf_.to_frame(origin_as_datetime=False).iloc[0].values.round(4))
 
     .. testoutput::
 
-        (1, 1, 6, 5)
+        [1.8531 1.3062 1.2332 1.1161 1.0444]
 
-    Apply ``future_trend`` when the completed lower triangle should reflect a
-    prospective trend assumption. Holding historical ``trend`` at zero isolates
-    the effect on the fitted future incrementals.
+    When the reserve will be carried for multiple future years, an actuary
+    may want to apply a prospective trend to the projected incremental
+    payments. Setting ``future_trend=0.1`` grows projected incremental
+    payments by 10% per period relative to the fitted historical pattern,
+    while ``trend=0`` keeps historical development unchanged so only the
+    forward projection is affected. The difference in the sum of all
+    projected incrementals shows how the future trend assumption flows
+    through to the total reserve.
 
     .. testcode::
 
diff --git a/chainladder/development/learning.py b/chainladder/development/learning.py
index ab603d98..29b924d0 100644
--- a/chainladder/development/learning.py
+++ b/chainladder/development/learning.py
@@ -53,12 +53,22 @@ class DevelopmentML(DevelopmentBase):
 
     Examples
     --------
-    Choose the response scale before fitting a machine-learning development
-    model. In this class, ``fit_incrementals=True`` means the scikit-learn model
-    is trained on actual incremental dollar amounts. Setting it to ``False``
-    trains on cumulative values instead. Both fitted models are translated back
-    into ``ldf_``, so the comparison is about the training target, not a request
-    to model age-to-age factors directly.
+    ``DevelopmentML`` bridges scikit-learn and the chainladder workflow:
+    it converts a Triangle to a tabular DataFrame, fits any sklearn-compatible
+    estimator or Pipeline against it, then converts the predictions back into
+    ``ldf_`` patterns usable with tail selection, ``Chainladder``, and other
+    methods.
+
+    ``fit_incrementals`` controls what the estimator is trained to predict.
+    When an actuary believes the period-to-period change in losses is more
+    predictable from development age than the cumulative total is, training
+    on incrementals (``fit_incrementals=True``) is more appropriate. When
+    the cumulative level is the more natural target, use
+    ``fit_incrementals=False``. Both options produce an ``ldf_``, but the
+    fitted pattern differs because the training target changes. (Whether the
+    incremental values represent dollar amounts depends on the estimator:
+    a ``LinearRegression`` with no log transform trains on raw dollar
+    increments, while a log-space transform would fit log-scale values.)
 
     .. testsetup::
 
@@ -85,17 +95,27 @@ class DevelopmentML(DevelopmentBase):
         m_cum = cl.DevelopmentML(
             pipe, y_ml=[tri.columns[0]], fit_incrementals=False
         ).fit(tri)
-        print(float(np.round(m_incr.ldf_.values[0, 0, 0, 0], 4)))
-        print(float(np.round(m_cum.ldf_.values[0, 0, 0, 0], 4)))
+        print(m_incr.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(m_cum.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
 
     .. testoutput::
 
-        3.508
-        3.515
-
-    Pass sample weights into the scikit-learn step when observations should not
-    contribute equally to the regression. ``weighted_step='model'`` forwards the
-    weights into the final estimator in the pipeline.
+        [3.508  1.7435 1.4379 1.1655]
+        [3.515  1.735  1.3993 1.152 ]
+
+    By default the regression treats each cell equally regardless of its loss
+    size. An actuary who prefers development patterns to reflect the
+    experience of larger accident years more heavily (similar to traditional
+    volume-weighted averages) can pass ``sample_weight`` to ``fit``. However,
+    passing ``sample_weight`` alone is not enough: without ``weighted_step``,
+    the weights are computed internally but never forwarded to the sklearn
+    estimator. ``weighted_step`` takes the name of a step in the Pipeline as
+    a string and routes the weights to that step using sklearn's
+    ``step_name__sample_weight`` fit-params convention. In the example below,
+    ``weighted_step='model'`` matches the name given to the
+    ``LinearRegression`` step, so it receives ``sample_weight=tri * tri``
+    (squared cumulative losses) during fitting, giving larger-loss accident
+    years proportionally more influence on the development pattern.
 
     .. testcode::
 
@@ -121,13 +141,52 @@ class DevelopmentML(DevelopmentBase):
             fit_incrementals=False,
             weighted_step="model",
         ).fit(tri, sample_weight=tri * tri)
-        print(float(np.round(m0.ldf_.values[0, 0, 0, 0], 4)))
-        print(float(np.round(m1.ldf_.values[0, 0, 0, 0], 4)))
+        print(m0.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(m1.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+
+    .. testoutput::
+
+        [3.515  1.735  1.3993 1.152 ]
+        [3.4459 1.7749 1.4053 1.1377]
+
+    Loss development has a natural time-series structure: the cumulative loss
+    at one age is directly related to the amount at the prior age. The
+    ``autoregressive`` parameter lets the model re-feed its own predictions
+    as lagged inputs when filling in the lower triangle. Each tuple in the
+    list specifies ``(column_name, lag, source_column)``: the Triangle column
+    ``source_column`` is shifted by ``lag`` periods and added to the tabular
+    DataFrame under ``column_name``. During prediction, each period's output
+    becomes the next period's lagged input automatically. The
+    ``column_name`` must also appear in the Pipeline formula so the estimator
+    can use it as a feature.
+
+    .. testcode::
+
+        import numpy as np
+        from sklearn.linear_model import LinearRegression
+        from sklearn.pipeline import Pipeline
+
+        from chainladder.utils.utility_functions import PatsyFormula
+
+        tri = cl.load_sample("genins")
+        col = tri.columns[0]
+        pipe = Pipeline(
+            steps=[
+                ("design_matrix", PatsyFormula("C(development) + lag_loss")),
+                ("model", LinearRegression(fit_intercept=False)),
+            ]
+        )
+        m = cl.DevelopmentML(
+            pipe,
+            y_ml=[col],
+            autoregressive=[("lag_loss", -1, col)],
+            fit_incrementals=False,
+        ).fit(tri)
+        print(m.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
 
     .. testoutput::
 
-        3.515
-        3.4459
+        [3.4815 1.6246 1.3017 1.0065]
 
     """
 
diff --git a/chainladder/development/outstanding.py b/chainladder/development/outstanding.py
index bc6ddbc5..0a27dee2 100644
--- a/chainladder/development/outstanding.py
+++ b/chainladder/development/outstanding.py
@@ -52,9 +52,22 @@ class CaseOutstanding(DevelopmentBase):
 
     Examples
     --------
-    ``paid_n_periods`` and ``case_n_periods`` control how many recent origin
-    years inform the ``Development`` weights that smooth the paid and case
-    patterns.
+    ``CaseOutstanding`` is appropriate when an actuary tracks both paid and
+    incurred losses and wants development patterns grounded in case reserve
+    movements rather than cumulative-to-cumulative ratios. The method models
+    incremental payments as a percentage of the prior period's case reserve
+    (``paid_ldf_``) and case reserve run-off as a percentage of the prior
+    period's case reserve (``case_ldf_``). The ``paid_to_incurred`` tuple
+    tells the estimator which triangle columns represent the paid and incurred
+    series.
+
+    When case reserving practices have changed recently (e.g., a company
+    tightened case adequacy standards in the last three years), including all
+    historical origin years in the pattern averages would mix old and new
+    reserving regimes. Setting ``paid_n_periods=3`` and ``case_n_periods=3``
+    restricts the averages to the three most recent origin years, reflecting
+    the current reserving philosophy. The shift in the first four ``paid_ldf_``
+    values below shows how much the pattern changes when older data is excluded.
 
     .. testsetup::
 
@@ -62,6 +75,8 @@ class CaseOutstanding(DevelopmentBase):
 
     .. testcode::
 
+        import numpy as np
+
         tri = cl.load_sample("usauto")
         all_years = cl.CaseOutstanding(
             paid_to_incurred=("paid", "incurred")
@@ -71,13 +86,13 @@ class CaseOutstanding(DevelopmentBase):
             paid_n_periods=3,
             case_n_periods=3,
         ).fit(tri)
-        print(round(float(all_years.paid_ldf_.values[0, 0, 0, 0]), 6))
-        print(round(float(three.paid_ldf_.values[0, 0, 0, 0]), 6))
+        print(all_years.paid_ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(three.paid_ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
 
     .. testoutput::
 
-        0.842814
-        0.833138
+        [0.8428 0.71   0.7084 0.6968]
+        [0.8331 0.7008 0.7139 0.7144]
 
     """
 
diff --git a/chainladder/methods/mack.py b/chainladder/methods/mack.py
index d42d6306..f2d9466b 100644
--- a/chainladder/methods/mack.py
+++ b/chainladder/methods/mack.py
@@ -42,9 +42,10 @@ class MackChainladder(Chainladder):
 
     Examples
     --------
-    Fit the Mack chainladder method and inspect the headline summary table,
-    which combines the deterministic chainladder estimate with Mack's
-    stochastic standard error.
+    Use ``MackChainladder`` when the IBNR point estimate alone is not
+    sufficient and a measure of reserve uncertainty is also needed.
+    ``summary_`` shows the deterministic chainladder ultimate alongside
+    Mack's per-origin prediction error.
 
     .. testsetup::
 
@@ -67,9 +68,10 @@ class MackChainladder(Chainladder):
         2012   9650.0   7162.150646  16812.150646    693.166178
         2013   6283.0  14396.919151  20679.919151    901.408385
 
-    The deterministic chainladder ultimates match those of
-    :class:`Chainladder`. Mack's contribution is the stochastic standard error
-    in the rightmost column, which can be aggregated across origins.
+    ``total_mack_std_err_`` aggregates the prediction error across all
+    origins. It exceeds the quadrature sum of the per-origin errors in
+    ``summary_`` because parameter risk is correlated across origins: all
+    origins share the same estimated age-to-age factors.
 
     .. testcode::
 
@@ -80,10 +82,12 @@ class MackChainladder(Chainladder):
         columns        values
         (Total,)  1424.531543
 
-    Mack's total error depends on how ``ldf_`` and ``sigma_`` were produced.
-    Here the same triangle is pre-smoothed with :class:`Development` using
-    ``average='simple'`` instead of the default volume weights before fitting
-    ``MackChainladder``, which raises the aggregate Mack standard error.
+    The Mack standard error is sensitive to how the upstream development
+    factors were estimated. Using simple (unweighted) averaging in
+    :class:`Development` before fitting ``MackChainladder`` gives equal
+    weight to each accident year regardless of size. On a small triangle
+    this raises the aggregate standard error relative to volume weighting,
+    since thinner years contribute more uncertainty.
 
     .. testcode::
 
@@ -114,8 +118,10 @@ def fit(self, X, y=None, sample_weight=None):
 
         Examples
         --------
-        Fitting attaches the ``ultimate_`` and Mack std error attributes to
-        the estimator and returns the estimator itself.
+        After fitting, ``ibnr_`` holds the point estimate per origin and
+        ``mack_std_err_`` holds the prediction error. The ratio of the two
+        gives a coefficient of variation that shows which origins carry the
+        most reserve uncertainty relative to their size.
 
         .. testsetup::
 
@@ -123,12 +129,24 @@ def fit(self, X, y=None, sample_weight=None):
 
         .. testcode::
 
+            import numpy as np
+
             tr = cl.load_sample('ukmotor')
-            cl.MackChainladder().fit(tr)
+            model = cl.MackChainladder().fit(tr)
+            print(model.ibnr_.to_frame(origin_as_datetime=False).round(1))
+            print(np.round(model.mack_std_err_.values[0, 0, :, -1], 1))
 
         .. testoutput::
 
-            MackChainladder()
+                     2261
+            2007      NaN
+            2008    350.9
+            2009   1037.5
+            2010   2044.9
+            2011   3663.4
+            2012   7162.2
+            2013  14396.9
+            [  nan  27.2  36.5 144.5 427.6 693.2 901.4]
         """
         super().fit(X, y, sample_weight)
         if "sigma_" not in self.X_:
@@ -168,8 +186,11 @@ def predict(self, X, sample_weight=None):
 
         Examples
         --------
-        Fit the model and apply it to a Triangle with the same shape, then
-        read the Mack standard error off the resulting Triangle.
+        ``predict`` re-applies the fitted age-to-age factors and sigma
+        estimates to a new triangle without refitting. A common use is
+        sensitivity testing: scale the reported losses by an adverse factor
+        and call ``predict`` to see how the Mack standard error responds,
+        holding the development pattern fixed.
 
         .. testsetup::
 
@@ -179,13 +200,16 @@ def predict(self, X, sample_weight=None):
 
             tr = cl.load_sample('ukmotor')
             model = cl.MackChainladder().fit(tr)
-            predicted = model.predict(tr)
-            print(predicted.total_mack_std_err_)
+            tr_adverse = tr * 1.05
+            print(model.predict(tr).total_mack_std_err_)
+            print(model.predict(tr_adverse).total_mack_std_err_)
 
         .. testoutput::
 
             columns        values
             (Total,)  1424.531543
+            columns        values
+            (Total,)  1475.539173
         """
         X_new = super().predict(X, sample_weight)
         X_new.sigma_ = getattr(X_new, "sigma_", self.X_.sigma_)

From 2f5ee76bf43613d15067800da8ea225b30fb2abd Mon Sep 17 00:00:00 2001
From: Ethan Kang <ethan.yskang@gmail.com>
Date: Mon, 18 May 2026 11:48:34 -0700
Subject: [PATCH 5/5] docs: align development estimator examples with user
 guide

Expand doctest Examples for TweedieGLM, CaseOutstanding,
IncrementalAdditive, DevelopmentML, and Barnett-Zehnwirth to match
guide narratives, reviewer feedback, and stable outputs.

Refs #704

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 chainladder/development/barnzehn.py    |  67 ++++++-----
 chainladder/development/glm.py         | 110 +++++++++++-------
 chainladder/development/incremental.py | 107 ++++++++++++------
 chainladder/development/learning.py    | 147 +++++++++++--------------
 chainladder/development/outstanding.py | 104 ++++++++++-------
 5 files changed, 300 insertions(+), 235 deletions(-)

diff --git a/chainladder/development/barnzehn.py b/chainladder/development/barnzehn.py
index 9e162d6a..7a1580bc 100644
--- a/chainladder/development/barnzehn.py
+++ b/chainladder/development/barnzehn.py
@@ -15,6 +15,12 @@ class BarnettZehnwirth(TweedieGLM):
     """ This estimator enables modeling from the Probabilistic Trend Family as
     described by Barnett and Zehnwirth.
 
+    The model is fit on log-incremental losses and produces multiplicative
+    ``ldf_`` patterns for use with IBNR estimators. Specify the regression
+    structure either with a patsy ``formula`` or with PTF period groupings
+    (``alpha``, ``gamma``, ``iota``) that define origin, trend, and
+    final-period cohorts.
+
     .. versionadded:: 0.8.2
 
     Parameters
@@ -35,31 +41,10 @@ class BarnettZehnwirth(TweedieGLM):
 
     Examples
     --------
-    Standard ``Development`` assumes loss development ratios are stable across
-    accident years. When a triangle shows a calendar-year effect (e.g., recent
-    diagonals are systematically heavier or lighter across all accident years
-    due to inflation or a shift in case reserving practices), that assumption
-    breaks down. ``BarnettZehnwirth`` addresses this by fitting a log-linear
-    model that decomposes the triangle into separate origin, development, and
-    calendar-year trend components.
-
-    Two interfaces let you specify the trend structure. The ``formula``
-    argument takes any patsy expression directly. The Probabilistic Trend
-    Family (PTF) arguments provide a structured shorthand: ``alpha`` groups
-    accident years that share the same level effect, ``gamma`` defines
-    breakpoints for a piecewise linear development-age trend, and ``iota``
-    defines breakpoints for a piecewise linear calendar-year (diagonal) trend.
-
-    The ``abc`` triangle has 11 accident years (1977-1987) and 11 development
-    ages (12-132 months). Suppose an actuary notices from the triangle (or
-    from external information such as a change in inflation or legal
-    environment) that accident years before 1982 behave differently from those
-    after, that development speed changes at the 36-month and 72-month marks,
-    and that calendar-year trends shift at two points in the diagonal sequence.
-    Those observations translate directly into ``alpha=[0, 5]``,
-    ``gamma=[0, 2, 5]``, and ``iota=[0, 7, 11]``. The first three fitted
-    coefficients differ from the unconstrained ``formula`` model, reflecting
-    the additional structure the actuary has imposed.
+    When many accident years are available but you want a smaller number of
+    origin cohorts, specify ``alpha``, ``gamma``, and ``iota`` instead of a
+    separate factor for every year. The fitted design has fewer parameters than
+    a fully saturated origin-by-development formula on the same triangle.
 
     .. testsetup::
 
@@ -67,22 +52,36 @@ class BarnettZehnwirth(TweedieGLM):
 
     .. testcode::
 
-        import numpy as np
-
         tri = cl.load_sample("abc")
-        m_formula = cl.BarnettZehnwirth(
-            formula="C(origin)+C(development)"
-        ).fit(tri)
         m_ptf = cl.BarnettZehnwirth(
             alpha=[0, 5], gamma=[0, 2, 5], iota=[0, 7, 11]
         ).fit(tri)
-        print(np.round(m_formula.coef_.values.flatten()[:3], 3))
-        print(np.round(m_ptf.coef_.values.flatten()[:3], 3))
+        m_full = cl.BarnettZehnwirth(
+            formula="C(origin)+C(development)"
+        ).fit(tri)
+        print(len(m_ptf.coef_.values.flatten()))
+        print(len(m_full.coef_.values.flatten()))
+
+    .. testoutput::
+
+        6
+        21
+
+    Use a patsy ``formula`` when the reserving structure needs explicit terms
+    (for example separate origin and development factors) rather than the PTF
+    cohort shorthand.
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("abc")
+        m = cl.BarnettZehnwirth(formula="C(origin)+C(development)").fit(tri)
+        print(np.round(m.ldf_.values[0, 0, :4, 0], 4))
 
     .. testoutput::
 
-        [11.837  0.179  0.345]
-        [12.151  0.274 -0.064]
+        [2.2854 2.2854 2.2854 2.2854]
 
     """
 
diff --git a/chainladder/development/glm.py b/chainladder/development/glm.py
index 20c9ca48..9db21142 100644
--- a/chainladder/development/glm.py
+++ b/chainladder/development/glm.py
@@ -11,12 +11,21 @@
 
 
 class TweedieGLM(DevelopmentBase):
-    """ This estimator creates development patterns with a GLM using a Tweedie distribution.
-
-    The Tweedie family includes several of the more popular distributions including
-    the normal, ODP poisson, and gamma distributions.  This class is a special case
-    of `DevleopmentML`.  It restricts to just GLM using a TweedieRegressor and
-    provides an R-like formulation of the design matrix.
+    """ GLM reserving with scikit-learn's Tweedie distribution.
+
+    Implements the GLM reserving structure of Taylor and McGuire. The Tweedie
+    family covers normal, ODP Poisson, gamma, and related targets via ``power``
+    and ``link``. Covariates from any triangle axis can enter through a patsy
+    ``design_matrix`` while staying close to traditional chainladder methods when
+    origin and development are coded categorically.
+
+    Triangles are converted to long-format tables internally (as with
+    ``Triangle.to_frame(keepdims=True)``); origin periods are restated as years
+    from the earliest origin for sklearn compatibility, and the response is
+    converted to an incremental basis before fitting. This class is a special
+    case of :class:`~chainladder.DevelopmentML` that uses only
+    :class:`~sklearn.linear_model.TweedieRegressor` behind a
+    :class:`~chainladder.utils.utility_functions.PatsyFormula` step.
 
     .. versionadded:: 0.8.1
 
@@ -29,7 +38,7 @@ class TweedieGLM(DevelopmentBase):
     design_matrix: formula-like
         A patsy formula describing the independent variables, X of the GLM
     response:  str
-        Column name for the reponse variable of the GLM.  If ommitted, then the
+        Column name for the response variable of the GLM. If omitted, then the
         first column of the Triangle will be used.
     power: float, default=1
             The power determines the underlying target distribution according
@@ -79,16 +88,11 @@ class TweedieGLM(DevelopmentBase):
 
     Examples
     --------
-    ``TweedieGLM`` with ``power=1`` and log link implements the
-    Over-Dispersed Poisson (ODP) model, the standard GLM equivalent of
-    volume-weighted chainladder. Its main advantage over column-by-column
-    development is that it fits the entire triangle simultaneously, which
-    allows parameter reduction: replacing categorical dummies (one per
-    accident year or development age) with continuous trend terms lowers the
-    parameter count while often staying close to the traditional result. The
-    default categorical design matrix replicates chainladder; switching to
-    continuous ``development + origin`` terms is a more parsimonious choice
-    when data are sparse.
+    Volume-weighted chainladder development can be replicated with a
+    Poisson-log GLM on incremental paid losses: categorical origin and
+    development in ``design_matrix``, ``power=1``, and ``link='log'``. The
+    resulting ``ldf_`` matches :class:`~chainladder.Development` closely on
+    ``genins``.
 
     .. testsetup::
 
@@ -99,42 +103,68 @@ class TweedieGLM(DevelopmentBase):
         import numpy as np
 
         tri = cl.load_sample("genins")
-        m_cat = cl.TweedieGLM(
-            power=1, design_matrix="C(development) + C(origin)"
-        ).fit(tri)
-        m_cont = cl.TweedieGLM(
-            power=1, design_matrix="development + origin"
+        odp = cl.TweedieGLM(
+            design_matrix="C(development) + C(origin)",
+            power=1,
+            link="log",
         ).fit(tri)
-        print(m_cat.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
-        print(m_cont.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        trad = cl.Development().fit(tri)
+        print(round(float(odp.ldf_.values[0, 0, 0, 0]), 4))
+        print(round(float(trad.ldf_.values[0, 0, 0, 0]), 4))
+        print(np.round(odp.ldf_.values[0, 0, :4, 0], 4))
 
     .. testoutput::
 
-        [3.491  1.7474 1.4574 1.1739]
-        [1.9277 1.4465 1.2864 1.2065]
+        3.491
+        3.4906
+        [3.491 3.491 3.491 3.491]
 
-    The ``power`` parameter selects the variance-mean relationship. For lines
-    where some development periods show no incremental activity, the compound
-    Poisson-Gamma distribution (``1 < power < 2``) is more appropriate than
-    pure Poisson (``power=1``). The ``alpha`` parameter (default 1.0) controls
-    L2 regularization strength: increasing it smooths fitted coefficients on
-    sparse or over-parameterised triangles. If the optimizer raises a
-    ``ConvergenceWarning``, increase ``max_iter`` (default 100).
+    Patsy R-style formulas set ``design_matrix``; continuous ``development``
+    and ``origin`` terms yield a small coefficient table via ``coef_``.
+
+    .. testcode::
+
+        tri = cl.load_sample("genins")
+        glm = cl.TweedieGLM(design_matrix="development + origin").fit(tri)
+        print(len(glm.coef_))
+        print(round(float(glm.coef_.iloc[0, 0]), 6))
+        print(round(float(glm.coef_.iloc[1, 0]), 6))
+
+    .. testoutput::
+
+        3
+        13.516322
+        -0.006251
+
+    On multi-LOB triangles, interaction terms can keep the model parsimonious
+    (10 coefficients here versus 18+ in a full categorical chainladder). The
+    percent difference in ``cdf_`` versus :class:`~chainladder.Development`
+    stays within about 1% at each ultimate lag:
 
     .. testcode::
 
         import numpy as np
 
-        tri = cl.load_sample("genins")
-        m_poisson = cl.TweedieGLM(power=1, link="log").fit(tri)
-        m_tweedie = cl.TweedieGLM(power=1.5, link="log").fit(tri)
-        print(m_poisson.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
-        print(m_tweedie.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        clrd = cl.load_sample("clrd")["CumPaidLoss"].groupby("LOB").sum()
+        clrd = clrd[clrd["LOB"].isin(["ppauto", "comauto"])]
+        dev = cl.TweedieGLM(
+            design_matrix=(
+                "LOB+LOB:C(np.minimum(development, 36))"
+                "+LOB:development+LOB:origin"
+            ),
+            max_iter=1000,
+        ).fit(clrd)
+        trad = cl.Development().fit(clrd)
+        pct = ((dev.cdf_.iloc[..., 0, :] / trad.cdf_) - 1).to_frame().round(3)
+        print(len(dev.coef_))
+        print(np.round(pct.loc["comauto"].values, 3))
+        print(np.round(pct.loc["ppauto"].values, 3))
 
     .. testoutput::
 
-        [3.491  1.7474 1.4574 1.1739]
-        [3.7032 1.7569 1.4591 1.1776]
+        10
+        [ 0.002  0.003 -0.01   0.003  0.011  0.008  0.005 -0.    -0.002]
+        [ 0.006  0.003 -0.     0.001  0.002  0.001  0.001  0.001  0.001]
 
     """
 
diff --git a/chainladder/development/incremental.py b/chainladder/development/incremental.py
index 31605c2f..d1801342 100644
--- a/chainladder/development/incremental.py
+++ b/chainladder/development/incremental.py
@@ -11,15 +11,26 @@
 class IncrementalAdditive(DevelopmentBase):
     """ The Incremental Additive Method.
 
+    This estimator implements the additive method of Schmidt (2006), Section 4.7:
+    expected incremental losses satisfy ``E[Z_{i,k}] = eta_i * gamma_k``, where
+    ``eta_i`` is exposure (``sample_weight``, e.g. premium) for accident year
+    ``i`` and ``gamma_k`` is an incremental loss ratio at development age ``k``
+    that is common to all accident years. The fitted ``zeta_`` estimates those
+    common ``gamma_k``; unobserved incrementals are completed as
+    ``zeta_ * sample_weight``. Dollar ``incremental_`` differ by origin because
+    exposure differs; implied multiplicative ``ldf_`` are derived from the
+    completed incremental triangle and can also differ by origin.
+
     Parameters
     ----------
     trend: float (default=0.0)
-        A multiplicative trend amount used to trend each incremental development
-        period the valuation_date of the Triangle.
+        Implementation extension (not in Schmidt, 2006): multiplicative trend
+        applied to incremental losses before ``zeta_`` is estimated, trending
+        each development period to the triangle valuation date.
     future_trend: float (default=None)
-        The trend to apply to the incremental development periods in the lower
-        half of the completed Triangle. If None, then will be set to the value of
-        the trend parameter.
+        Implementation extension: trend applied when projecting incrementals
+        beyond the valuation date into the lower triangle. If None, uses
+        ``trend``.
     n_periods: integer, optional (default=-1)
         number of origin periods to be used in the ldf average calculation. For
         all origin periods, set n_periods=-1
@@ -54,12 +65,12 @@ class IncrementalAdditive(DevelopmentBase):
         The raw incrementals as a percent of exposure trended to the valuation
         date of the Triangle. Only those used in the fitting.
     zeta_: Triangle
-        The fitted incrementals as a percent of exposure trended to the valuation
-        date of the Triangle.
+        Fitted incremental loss ratios ``gamma_k`` (common across accident years)
+        as a percent of exposure, trended to the valuation date of the Triangle.
     cum_zeta_: Triangle
         The fitted cumulative percent of exposure trended to the valuation date of 
         the Triangle
-    w_: ndarray
+    w_ : ndarray
         The weight used in the zeta fitting
     w_tri_: Triangle
         Triangle of w_
@@ -71,15 +82,10 @@ class IncrementalAdditive(DevelopmentBase):
 
     Examples
     --------
-    Use ``IncrementalAdditive`` when incremental losses per unit of exposure
-    are expected to be more stable across accident years than age-to-age
-    factors. This is common in lines where claim payments are driven by
-    exposure volume: for example, a workers' compensation book where
-    claim costs per payroll dollar are steady even as the cumulative reported
-    losses grow at different rates by accident year. Passing the latest
-    exposure diagonal as ``sample_weight`` normalises each period's
-    incremental by exposure before averaging. The fitted ``ldf_`` slots
-    directly into ``Chainladder`` or any other method.
+    Schmidt (2006), Example F, uses the ``ia_sample`` triangle: cumulative
+    ``loss`` with latest ``exposure`` as ``sample_weight`` (premiums). Fitted
+    ``incremental_`` are dollars by origin and age; ``zeta_`` is one pattern
+    shared across origins; implied ``ldf_`` can still vary by origin.
 
     .. testsetup::
 
@@ -90,40 +96,71 @@ class IncrementalAdditive(DevelopmentBase):
         import numpy as np
 
         tri = cl.load_sample("ia_sample")
-        model = cl.IncrementalAdditive().fit(
+        ia = cl.IncrementalAdditive().fit(
+            tri["loss"], sample_weight=tri["exposure"].latest_diagonal
+        )
+        print(np.round(ia.incremental_.values[0, 0, -1, :], 0))
+        print(np.round(ia.ldf_.values[0, 0, :3, :3], 4))
+
+    .. testoutput::
+
+        [1889. 1811. 1256. 1157.  740.  300.]
+        [[1.8531 1.3062 1.2332]
+         [1.8895 1.3191 1.2336]
+         [1.9233 1.3288 1.2301]]
+
+    A volume-weighted estimate of the common ``gamma_k`` across origins,
+    multiplied by latest exposure, reproduces the fitted incrementals in the
+    lower triangle (here at age 72), as in Schmidt's additive predictors.
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("ia_sample")
+        ia = cl.IncrementalAdditive().fit(
             tri["loss"], sample_weight=tri["exposure"].latest_diagonal
         )
-        print(model.ldf_.to_frame(origin_as_datetime=False).iloc[0].values.round(4))
+        zeta = tri["loss"].cum_to_incr().sum("origin") / tri["exposure"].sum("origin")
+        projected = (
+            zeta.values[0, 0, 0, -1]
+            * tri["exposure"].latest_diagonal.values[0, 0, -1, 0]
+        )
+        fitted = ia.incremental_.values[0, 0, -1, -1]
+        print(np.isclose(projected, fitted))
 
     .. testoutput::
 
-        [1.8531 1.3062 1.2332 1.1161 1.0444]
+        True
 
-    When the reserve will be carried for multiple future years, an actuary
-    may want to apply a prospective trend to the projected incremental
-    payments. Setting ``future_trend=0.1`` grows projected incremental
-    payments by 10% per period relative to the fitted historical pattern,
-    while ``trend=0`` keeps historical development unchanged so only the
-    forward projection is affected. The difference in the sum of all
-    projected incrementals shows how the future trend assumption flows
-    through to the total reserve.
+    The ``trend`` and ``future_trend`` parameters are not part of Schmidt
+    (2006); they are chainladder extensions for trending incrementals before
+    fitting ``zeta_`` and when projecting the lower triangle. The effect is
+    material on projected dollars (not on cumulative link-ratio semantics).
 
     .. testcode::
 
         import numpy as np
 
         tri = cl.load_sample("ia_sample")
-        loss = tri["loss"]
         sw = tri["exposure"].latest_diagonal
-        m0 = cl.IncrementalAdditive(trend=0, future_trend=0).fit(loss, sample_weight=sw)
-        m1 = cl.IncrementalAdditive(trend=0, future_trend=0.1).fit(loss, sample_weight=sw)
-        print(float(np.round(np.nansum(m0.incremental_.values), 1)))
-        print(float(np.round(np.nansum(m1.incremental_.values), 1)))
+        base = cl.IncrementalAdditive().fit(tri["loss"], sample_weight=sw)
+        trended = cl.IncrementalAdditive(trend=0.02, future_trend=0.05).fit(
+            tri["loss"], sample_weight=sw
+        )
+        print(float(np.round(base.incremental_.values[0, 0, -1, -1], 0)))
+        print(float(np.round(trended.incremental_.values[0, 0, -1, -1], 0)))
 
     .. testoutput::
 
-        30988.1
-        33360.1
+        300.0
+        383.0
+
+    References
+    ----------
+    Schmidt, K. (2006). Methods and Models of Loss Reserving Based on Run-Off
+    Triangles: A Unifying Survey. CAS Forum, Fall 2006, Section 4.7 (Additive
+    Method). https://www.casact.org/sites/default/files/database/forum_06fforum_273.pdf
 
     """
 
diff --git a/chainladder/development/learning.py b/chainladder/development/learning.py
index 29b924d0..29f9f1c4 100644
--- a/chainladder/development/learning.py
+++ b/chainladder/development/learning.py
@@ -12,27 +12,32 @@
 
 
 class DevelopmentML(DevelopmentBase):
-    """ A Estimator that interfaces with machine learning (ML) tools that implement
-    the scikit-learn API.
+    """ Interface to scikit-learn estimators for loss development patterns.
 
-    The `DevelopmentML` estimator is used to generate ``ldf_`` patterns from
-    the data.
+    ``DevelopmentML`` lets reserving workflows use any sklearn-compatible
+    regressor (often inside a :class:`~sklearn.pipeline.Pipeline`). It converts
+    a :class:`~chainladder.Triangle` to a tabular design matrix, fits the ML
+    model, predicts through the terminal development age to complete the lower
+    triangle, and expresses the result as ``ldf_`` for tails and IBNR methods.
+    :class:`~chainladder.TweedieGLM` is a special case with
+    :class:`~sklearn.linear_model.TweedieRegressor` as the only ML step.
 
     .. versionadded:: 0.8.1
 
 
     Parameters
     ----------
-    estimator_ml: skearn Estimator
+    estimator_ml: sklearn Estimator
         Any sklearn compatible regression estimator, including Pipelines and
     y_ml: list or str or sklearn_transformer
         The response column(s) for the machine learning algorithm. It must be
         present within the Triangle.
-    autoregressive: tuple, (autoregressive_col_name, lag, source_col_name)
-        The subset of response column(s) to use as lagged features for the
-        Time Series aspects of the model. Predictions from one development period
-        get used as featues in the next development period. Lags should be negative
-        integers.
+    autoregressive: list of tuple
+        Each tuple is ``(feature_name, lag, source_column)``. ``feature_name`` must
+        also appear in the pipeline design matrix. ``DevelopmentML`` fills that
+        column with lagged ``source_column`` values and, when projecting forward,
+        replaces it with the prior development period's prediction. Lags should be
+        negative integers (for example ``-12`` on a monthly triangle is one year).
     weight_step: str
         Step name within estimator_ml that is weighted
     drop: tuple or list of tuples
@@ -53,22 +58,11 @@ class DevelopmentML(DevelopmentBase):
 
     Examples
     --------
-    ``DevelopmentML`` bridges scikit-learn and the chainladder workflow:
-    it converts a Triangle to a tabular DataFrame, fits any sklearn-compatible
-    estimator or Pipeline against it, then converts the predictions back into
-    ``ldf_`` patterns usable with tail selection, ``Chainladder``, and other
-    methods.
-
-    ``fit_incrementals`` controls what the estimator is trained to predict.
-    When an actuary believes the period-to-period change in losses is more
-    predictable from development age than the cumulative total is, training
-    on incrementals (``fit_incrementals=True``) is more appropriate. When
-    the cumulative level is the more natural target, use
-    ``fit_incrementals=False``. Both options produce an ``ldf_``, but the
-    fitted pattern differs because the training target changes. (Whether the
-    incremental values represent dollar amounts depends on the estimator:
-    a ``LinearRegression`` with no log transform trains on raw dollar
-    increments, while a log-space transform would fit log-scale values.)
+    Features from any triangle axis can enter an sklearn
+    :class:`~sklearn.compose.ColumnTransformer` or
+    :class:`~sklearn.pipeline.Pipeline`. On ``clrd`` grouped by line of business,
+    one-hot-encode ``LOB`` and ``development``, pass ``origin`` through, and fit
+    a linear model (the user guide uses ``RandomForestRegressor`` the same way).
 
     .. testsetup::
 
@@ -77,45 +71,37 @@ class DevelopmentML(DevelopmentBase):
     .. testcode::
 
         import numpy as np
+        from sklearn.compose import ColumnTransformer
         from sklearn.linear_model import LinearRegression
         from sklearn.pipeline import Pipeline
+        from sklearn.preprocessing import OneHotEncoder
 
-        from chainladder.utils.utility_functions import PatsyFormula
-
-        tri = cl.load_sample("genins")
-        pipe = Pipeline(
+        clrd = cl.load_sample("clrd").groupby("LOB").sum()["CumPaidLoss"]
+        design_matrix = ColumnTransformer(
+            transformers=[
+                ("dummy", OneHotEncoder(drop="first"), ["LOB", "development"]),
+                ("passthrough", "passthrough", ["origin"]),
+            ]
+        )
+        estimator_ml = Pipeline(
             steps=[
-                ("design_matrix", PatsyFormula("C(development)")),
-                ("model", LinearRegression(fit_intercept=False)),
+                ("design_matrix", design_matrix),
+                ("model", LinearRegression()),
             ]
         )
-        m_incr = cl.DevelopmentML(
-            pipe, y_ml=[tri.columns[0]], fit_incrementals=True
-        ).fit(tri)
-        m_cum = cl.DevelopmentML(
-            pipe, y_ml=[tri.columns[0]], fit_incrementals=False
-        ).fit(tri)
-        print(m_incr.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
-        print(m_cum.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        m = cl.DevelopmentML(estimator_ml=estimator_ml, y_ml="CumPaidLoss").fit(
+            clrd
+        )
+        print(m.ldf_.shape)
+        print(np.round(m.ldf_.values[0, 0, 0, :4], 4))
 
     .. testoutput::
 
-        [3.508  1.7435 1.4379 1.1655]
-        [3.515  1.735  1.3993 1.152 ]
-
-    By default the regression treats each cell equally regardless of its loss
-    size. An actuary who prefers development patterns to reflect the
-    experience of larger accident years more heavily (similar to traditional
-    volume-weighted averages) can pass ``sample_weight`` to ``fit``. However,
-    passing ``sample_weight`` alone is not enough: without ``weighted_step``,
-    the weights are computed internally but never forwarded to the sklearn
-    estimator. ``weighted_step`` takes the name of a step in the Pipeline as
-    a string and routes the weights to that step using sklearn's
-    ``step_name__sample_weight`` fit-params convention. In the example below,
-    ``weighted_step='model'`` matches the name given to the
-    ``LinearRegression`` step, so it receives ``sample_weight=tri * tri``
-    (squared cumulative losses) during fitting, giving larger-loss accident
-    years proportionally more influence on the development pattern.
+        (6, 1, 10, 9)
+        [1.7448 0.9854 0.8117 0.6495]
+
+    ``fit_incrementals`` chooses whether the ML response is built from
+    incremental or cumulative triangle values before ``ldf_`` is derived.
 
     .. testcode::
 
@@ -132,33 +118,23 @@ class DevelopmentML(DevelopmentBase):
                 ("model", LinearRegression(fit_intercept=False)),
             ]
         )
-        m0 = cl.DevelopmentML(
+        m_incr = cl.DevelopmentML(
+            pipe, y_ml=[tri.columns[0]], fit_incrementals=True
+        ).fit(tri)
+        m_cum = cl.DevelopmentML(
             pipe, y_ml=[tri.columns[0]], fit_incrementals=False
         ).fit(tri)
-        m1 = cl.DevelopmentML(
-            pipe,
-            y_ml=[tri.columns[0]],
-            fit_incrementals=False,
-            weighted_step="model",
-        ).fit(tri, sample_weight=tri * tri)
-        print(m0.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
-        print(m1.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(float(np.round(m_incr.ldf_.values[0, 0, 0, 0], 4)))
+        print(float(np.round(m_cum.ldf_.values[0, 0, 0, 0], 4)))
 
     .. testoutput::
 
-        [3.515  1.735  1.3993 1.152 ]
-        [3.4459 1.7749 1.4053 1.1377]
+        3.508
+        3.515
 
-    Loss development has a natural time-series structure: the cumulative loss
-    at one age is directly related to the amount at the prior age. The
-    ``autoregressive`` parameter lets the model re-feed its own predictions
-    as lagged inputs when filling in the lower triangle. Each tuple in the
-    list specifies ``(column_name, lag, source_column)``: the Triangle column
-    ``source_column`` is shifted by ``lag`` periods and added to the tabular
-    DataFrame under ``column_name``. During prediction, each period's output
-    becomes the next period's lagged input automatically. The
-    ``column_name`` must also appear in the Pipeline formula so the estimator
-    can use it as a feature.
+    Autoregressive features use prior development predictions as covariates.
+    The lag column must be named in both ``autoregressive`` and the pipeline
+    (for example in a :class:`~chainladder.PatsyFormula`).
 
     .. testcode::
 
@@ -168,25 +144,28 @@ class DevelopmentML(DevelopmentBase):
 
         from chainladder.utils.utility_functions import PatsyFormula
 
-        tri = cl.load_sample("genins")
+        tri = cl.load_sample("raa")
         col = tri.columns[0]
         pipe = Pipeline(
             steps=[
-                ("design_matrix", PatsyFormula("C(development) + lag_loss")),
+                (
+                    "design_matrix",
+                    PatsyFormula("C(development) + pred_lag"),
+                ),
                 ("model", LinearRegression(fit_intercept=False)),
             ]
         )
         m = cl.DevelopmentML(
             pipe,
-            y_ml=[col],
-            autoregressive=[("lag_loss", -1, col)],
-            fit_incrementals=False,
+            y_ml=col,
+            fit_incrementals=True,
+            autoregressive=[("pred_lag", -12, col)],
         ).fit(tri)
-        print(m.ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(float(np.round(m.ldf_.values[0, 0, 0, 0], 4)))
 
     .. testoutput::
 
-        [3.4815 1.6246 1.3017 1.0065]
+        3.0297
 
     """
 
diff --git a/chainladder/development/outstanding.py b/chainladder/development/outstanding.py
index 0a27dee2..4099eca0 100644
--- a/chainladder/development/outstanding.py
+++ b/chainladder/development/outstanding.py
@@ -11,15 +11,18 @@
 
 
 class CaseOutstanding(DevelopmentBase):
-    """A determinisic method based on outstanding case reserves.
+    """ Deterministic development from prior-lag case reserves.
 
-    The CaseOutstanding method is a deterministic approach that develops
-    patterns of incremental payments as a percent of previous period case
-    reserves as well as patterns for case reserves as a percent of previous
-    period case reserves.  Although the patterns produces by the approach
-    approximate incremental payments and case outstanding, they are converted
-    into comparable multiplicative patterns for usage with the various IBNR
-    methods.
+    Estimates incremental paid amounts and case-reserve runoff as fractions of
+    the prior lag's carried case reserve. Like
+    :class:`~chainladder.MunichAdjustment` and
+    :class:`~chainladder.BerquistSherman`, this is useful when case reserves
+    should inform paid ultimates. A triangle with both paid and incurred columns
+    is required.
+
+    The incremental ``paid_ldf_`` patterns are not multiplicative link ratios;
+    the estimator also builds origin-specific implied multiplicative ``ldf_``
+    so standard IBNR methods can be applied.
 
     .. versionadded:: 0.8.0
 
@@ -33,41 +36,28 @@ class CaseOutstanding(DevelopmentBase):
         all origin periods, set paid_n_periods=-1
     case_n_periods: integer, optional (default=-1)
         number of origin periods to be used in the case pattern averages. For
-        all origin periods, set paid_n_periods=-1
+        all origin periods, set case_n_periods=-1
 
     Attributes
     ----------
     ldf_: Triangle
-        The estimated (multiplicative) loss development patterns.
+        Implied multiplicative loss development patterns (by paid/incurred
+        column); each origin period has its own pattern.
     cdf_: Triangle
         The estimated (multiplicative) cumulative development patterns.
     case_to_prior_case_: Triangle
-        The case to prior case ratios used for fitting the estimator
-    case_ldf_:
-        The selected case to prior case ratios of the fitted estimator
+        Case-to-prior-case incremental ratios by origin (for review).
+    case_ldf_: Triangle
+        Selected case-to-prior-case ratios averaged across origins.
     paid_to_prior_case_: Triangle
-        The paid to prior case ratios used for fitting the estimator
-    paid_ldf_:
-        The selected paid to prior case ratios of the fitted estimator
+        Paid-to-prior-case incremental ratios by origin (for review).
+    paid_ldf_: Triangle
+        Selected paid-to-prior-case ratios averaged across origins.
 
     Examples
     --------
-    ``CaseOutstanding`` is appropriate when an actuary tracks both paid and
-    incurred losses and wants development patterns grounded in case reserve
-    movements rather than cumulative-to-cumulative ratios. The method models
-    incremental payments as a percentage of the prior period's case reserve
-    (``paid_ldf_``) and case reserve run-off as a percentage of the prior
-    period's case reserve (``case_ldf_``). The ``paid_to_incurred`` tuple
-    tells the estimator which triangle columns represent the paid and incurred
-    series.
-
-    When case reserving practices have changed recently (e.g., a company
-    tightened case adequacy standards in the last three years), including all
-    historical origin years in the pattern averages would mix old and new
-    reserving regimes. Setting ``paid_n_periods=3`` and ``case_n_periods=3``
-    restricts the averages to the three most recent origin years, reflecting
-    the current reserving philosophy. The shift in the first four ``paid_ldf_``
-    values below shows how much the pattern changes when older data is excluded.
+    On ``usauto``, incremental paid in 12–24 is about 84% of case outstanding
+    at lag 12 (first entry in ``paid_ldf_`` at development 24–36):
 
     .. testsetup::
 
@@ -78,21 +68,51 @@ class CaseOutstanding(DevelopmentBase):
         import numpy as np
 
         tri = cl.load_sample("usauto")
-        all_years = cl.CaseOutstanding(
+        model = cl.CaseOutstanding(
             paid_to_incurred=("paid", "incurred")
         ).fit(tri)
-        three = cl.CaseOutstanding(
-            paid_to_incurred=("paid", "incurred"),
-            paid_n_periods=3,
-            case_n_periods=3,
-        ).fit(tri)
-        print(all_years.paid_ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
-        print(three.paid_ldf_.to_frame(origin_as_datetime=False).iloc[0].values[:4].round(4))
+        print(np.round(model.paid_ldf_.values[0, 0, 0, :4], 4))
 
     .. testoutput::
 
         [0.8428 0.71   0.7084 0.6968]
-        [0.8331 0.7008 0.7139 0.7144]
+
+    Implied multiplicative ``ldf_`` differ by accident year; the 1998 origin
+    paid pattern is shown below (compare to volume-weighted chainladder).
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("usauto")
+        model = cl.CaseOutstanding(
+            paid_to_incurred=("paid", "incurred")
+        ).fit(tri)
+        print(np.round(model.ldf_["paid"].values[0, 0, 0, :4], 4))
+
+    .. testoutput::
+
+        [1.7925 1.2056 1.0956 1.0457]
+
+    Review origin-level ``paid_to_prior_case_`` and ``case_to_prior_case_``
+    when tuning ``paid_n_periods`` and ``case_n_periods``; fitted selections
+    appear in ``paid_ldf_`` and ``case_ldf_``.
+
+    .. testcode::
+
+        import numpy as np
+
+        tri = cl.load_sample("usauto")
+        model = cl.CaseOutstanding(
+            paid_to_incurred=("paid", "incurred")
+        ).fit(tri)
+        print(np.round(model.case_to_prior_case_.values[0, 0, 0, :4], 4))
+        print(np.round(model.case_ldf_.values[0, 0, 0, :4], 4))
+
+    .. testoutput::
+
+        [0.5378 0.5541 0.5253 0.4981]
+        [0.534  0.5638 0.5296 0.49  ]
 
     """
 
@@ -110,7 +130,7 @@ def fit(self, X, y=None, sample_weight=None):
         Parameters
         ----------
         X : Triangle
-            Set of LDFs to which the munich adjustment will be applied.
+            Triangle with paid and incurred columns for ``paid_to_incurred``.
         y : Ignored
         sample_weight : Ignored