The Workers Compensation Tails

1. Introduction

Workers compensation is a very long-tailed line of insurance, since policies written today may generate payouts for many decades. Thus, the degree of uncertainty surrounding the ultimate loss is significant. The most important input to rational decisions in ratemaking and reserving is historical data, from which the process that generates future losses may be discerned. Unfortunately, loss triangles that extend beyond the cessation of payments in the first accident year (which may occur after 60 development years or more) are hard to come by, which explains the dearth of public knowledge about the trajectory of incremental payments for very mature claims.

This paper attempts to improve the state of information on the payment patterns of workers compensation triangles by analyzing a set of two triangles, the development years of which extend into a maturity range where these triangles start exhibiting incremental payments at zero amounts. The first triangle concerns indemnity benefits and has been made available by SCF (State Compensation Fund) Arizona. The second triangle comprises the medical component of permanent disability (PD) claims and has been provided by SAIF (State Accident Insurance Fund) Oregon.

The SAIF triangle has previously been studied by Sherman and Diss (2006). These authors pose the hypothesis of a bulge in incremental payments at high maturities, as may be caused by “added costs of caring for the elderly.” The analysis presented here does not support this hypothesis.

The following section outlines important properties of the Bayesian loss development model that has been devised for the analysis; a generalized version of this model has been made available in the R software package lossDev. Then, in a subsequent section, the model is applied to the SCF Arizona indemnity triangle (Section 3.1) and the SAIF Oregon triangle (Section 3.2). These two empirical applications of the model offer important insights into the development patterns of indemnity and medical claims payments at high maturities. Section 4 concludes.

2. The model

The statistical model of loss development employed in the analysis is Bayesian and is estimated using MCMC (Markov chain Monte Carlo simulation). The model treats incremental payments as a three-dimensional time series problem, similar to the approach taken by Barnett and Zehnwirth (2000) and Zehnwirth (1994). Specifically, the incremental payments are driven by three time processes, which manifest themselves in exposure growth, development, and the calendar year effect; these processes are illustrated in Figure 1.

Figure 1.Time processes in loss development

In the model, the growth rate that represents the calendar year effect is denoted κ. The rate of exposure growth, η, is net of the calendar year effect. The growth rate δ is the rate of decay in incremental payments, adjusted for the calendar year effect. Incremental payments that have been adjusted for the calendar year effect (and, hence, have been adjusted for inflation) represent consumption of units of (for instance, medical) services. A decline in consumption at the level of the aggregate loss triangle may be due to claimants exiting or remaining claimants decreasing their consumption.

What follows is a description of the main components of the model. After describing the likelihood, there follows a discussion of the modeling of the calendar year effect and of exposure growth. Further, there are details on how to model a structural break in consumption and how to employ Reversible Jump MCMC for estimating the trajectory of the consumption path without and with a structural break. The section concludes with a presentation of the estimation of the probability of payment. Note that both triangles studied here consist of accident year data.

2.1. The Student’s t likelihood

The loss development model fits to logarithmic incremental payments; non-positive incremental payments are treated as missing values. The model uses a Student’s t likelihood to accommodate heavy tails; this t distribution is implemented as a scale mixture of normal distributions (Gelman et al. 2004, 446); further, the degrees of freedom of this t distribution are determined within the model, but are restricted to values greater than 2 to ensure the existence of a finite variance. Finally, this t distribution accommodates heteroskedasticity by allowing the scale parameter to vary with development time; this time variation is estimated using a second-order random walk smoother (Congdon 2005, 141).

The Student’s t likelihood consists of

$$y_{i, j} \sim \text{N}\left(\mu_{i, j}, \omega_{i, j} \cdot \tau_{j}\right) \tag{1.a}$$

$$\mathbf{\mu}_{i, j}=\sum_{k=2}^{i+j} \mathbf{\kappa}_{k}+\sum_{k=1}^{i} \mathbf{\eta}_{k}+S_{j} \tag{1.b}$$

$$\omega_{i, j} \sim \mathrm{Ga}(v / 2, v / 2) \tag{1.c}$$

$$\tau_{j}=1 / \sigma_{j}^{2} \tag{1.d}$$

$$\log \left(\sigma_{j}^{2}\right) \sim \mathbf{N}\left(2 \cdot \log \left(\sigma_{j-1}^{2}\right)-\log \left(\sigma_{j-2}^{2}\right), v\right), j \geq 3 \tag{1.e}$$

$$v \sim \mathrm{Ga}(5,0.5) \tag{1.f}$$

$$\sigma_{j} \sim \mathrm{U}(0,10), j=1,2 \tag{1.g}$$

$$v \sim \chi_{(8)}^{2} \mathrm{I}(2,50) \tag{1.h}$$

where the indexes $i, j$ indicate the row (exposure year) and the column (development year) of the triangle, respectively, and the $\log (\cdot)$ operator refers to the natural logarithm. Equation (1a) states that the log incremental payments $\left(y_{i, j}\right)$ are modeled as a scale mixture of normal distributions. The expected value of a given $\log$ incremental payment $\left(\mu_{i, j}\right)$, as specified in Equation (1b), equals the calendar year effect-adjusted incremental payment of the corresponding column in the first row of the triangle $\left(S_j\right)$, plus the calendar year effects $\left(\kappa_{i+j}\right)$ and the exposure growth $\left(\eta_i\right)$ that has accumulated since the first payment in the first accident year $\left(y_{1,1}\right)$. Note that $\eta_1=0$, since there is no exposure year 0 . Similarly, $\kappa_1=\kappa_2=0$; this is because the index of the calendar year effect refers to the number of the diagonal, and the first calendar year effect is the one for diagonal $i+j=3$ (relative to diagonal $i+j=2$ ). The precision (i.e., inverse of the variance) of the normal distribution is the product of a gamma-distributed random variable $\left(\omega_{i,}\right)$, as defined in Equation (1c), and the inverse of the square of a scale parameter $\left(\sigma_i\right)$ that varies by development year and is defined in Equation (1d). The logarithm of the square of this scale parameter follows a second-order random walk, starting in development year 3, as described in Equation (1e). The innovation variance of this random walk is modeled as an inverse gamma distribution with smoothing parameters 5 and 0.5; see Equation (1f). Equation (1g) presents the prior distributions of the scale parameters of the first two development years. Finally, Equation (1h) defines the prior for the degrees of freedom; following Meyer and Yu (2000), this prior is a χ²₍₈₎ distribution that is restricted to the interval (2,50).

The variable $\left(S_j\right)$ represents (in logarithmic terms) the consumption of (medical or indemnity) services in development year $j$, normalized to the level of exposure of the first accident year. The trajectory of consumption from the first to the final $(j=K)$ development year, $S_1, S_2 \ldots S_K$, is estimated using a linear spline, the number of knots of which are estimated using Reversible Jump MCMC, as is to be discussed in Sections 2.5 and 2.6. The (logarithmic) rate of decay $\left(\delta_j\right)$ is backed out of $S_j$:

$$\delta_{j}=S_{j}-S_{j-1}, j \geq 2 \tag{1.i}$$

Clearly, because there is no development year 0, δ₁ = 0. Solving Equation (1i) for $S_j$ and inserting it into Equation (1b) delivers

$$\mu_{i, j}=S_{1}+\sum_{k=2}^{i+j} \kappa_{k}+\sum_{k=1}^{i} \eta_{k}+\sum_{k=1}^{j} \delta_{k} \tag{1.b′}$$

which states that the expected log incremental payment equals the first expected log payment (located in the first row and first column) (S₁), plus the accumulated calendar year effects $\left(\sum_{k=2}^{i+j} \mathbf{K}_{k}\right)$ the cumulative exposure growth $\left(\sum_{k=1}^{i} \eta_{k}\right),$ and the cumulative decay in consumption $\left(\sum_{k=1}^{j} \mathbf{S}_{k}\right)$

The growth rates κ, η, and δ are transformed from the logarithmic scale to the raw scale (by means of exponentiation) for each individual draw of the MCMC updates. The charts in Section 3 of the decay rate, of the exposure growth rate, and of the calendar year effect all rest on growth rates that have been transformed in such a way.

2.2. The calendar year effect

The calendar year effect (κ) is modeled as a normal distribution around an expert prior (µ_κ), the standard deviation of which has a uniform prior:

$$\kappa_{i+j} \sim \mathrm{N}\left(\mu_{\mathrm{\kappa}, i+j}, \tau_{\mathrm{k}}\right), i+j \geq 3 \tag{2.a}$$

$$\tau_{\mathrm{k}}=1 / \sigma_{\mathrm{k}}^{2} \tag{2.b}$$

$$\sigma_{\mathrm{k}} \sim \mathrm{U}(0,10) \tag{2.c}$$

For indemnity triangles, an appropriate expert prior is the legally stipulated rate of escalation; for Arizona, this rate is zero. For medical triangles, a suitable expert prior is the (logarithmic) rate of inflation of the Medical Care component of the CPI (Consumer Price Index), M-CPI for short. It can be shown that any systematic difference between the logarithmic rate of M-CPI inflation and the (actual but unknown) logarithmic rate of inflation that applies to medical payments on workers compensation claims feeds into the rates of exposure growth and the rates of decay without affecting the estimated log incremental payments. This is because we can rewrite Equation (1b′) as

$$\begin{align} \mu_{i, j}&=S_{1}+\sum_{k=2}^{i+j}\left(\kappa_{k}-c\right)+\sum_{k=1}^{i}\left(\eta_{k}+c\right) +\sum_{k=1}^{j}\left(\delta_{k}+c\right) \end{align} \tag{1.b″}$$

Let us assume that (in logarithmic terms) the M-CPI is systematically below the workers compensation price index by a constant $c>0$; this means that the posterior for the calendar year effect centers on a value that is systematically underestimated. In this case, the constant $c$ feeds into both the rate of decay $\left(\widetilde{\delta}_k=\delta_k+c\right)$ and the rate of exposure growth $\left(\tilde{\eta}_k=\eta_k+c\right)$, thereby leaving the expected log incremental payments $\left(\mu_{i, j}\right)$ unaffected. This arithmetic relation among the three growth rates reduces the discussion about the adequacy of the M-CPI rate of inflation as an expert prior for the workers compensation medical rate of inflation to a matter of interpretation regarding the rate of decay $(\delta)$ and the rate of calendar year effectadjusted exposure growth $(\eta)$.

For the purpose of squaring the triangle and, if applicable, estimating future exposure years, it is necessary to simulate future calendar year effects. Such simulation necessitates an assumption about the future behavior of the expert prior around which the future calendar year effect is normally distributed with the posterior precision τ_κ. For the indemnity triangle, the expert prior of future calendar year effects is simply set to zero. For the medical triangle, the expert prior is the simulated logarithmic rate of M-CPI inflation. This rate of inflation is specified as an Ornstein-Uhlenbeck process, similar to long-term interest rates (see Vasicek 1977). The Ornstein-Uhlenbeck process models the logarithmic rate of inflation as a stochastic process with transitory random shocks, as opposed to the random walk, where all innovations are permanent. As a result, the simulation of future rates of inflation using the Ornstein-Uhlenbeck process implies a return of the inflation rate to its long-term mean.

In discrete time, the Ornstein-Uhlenbeck process can be calibrated by estimating the following first-order autoregressive model (see Dixit and Pindyck 1994):

$$\mu_{\mathrm{\kappa}, {i}+1}=\mathrm{N}\left({a} \cdot \mu_{\mathrm{k}, i}+{b}, \tau_{\varepsilon}\right), 0<{a}<1 \tag{2.d}$$

where a, b, and the precision τ_ε are the parameters to be estimated, and µ_κ,i are the observed (log) rates of M-CPI inflation. The rate of inflation µ_κ,i reverts to the mean b/(1 − a) at the rate −log(a). The innovation variance of the shocks reads (see Dixit and Pindyck 1994)

$$\frac{1}{\tau_{\varepsilon}} \cdot \frac{-2 \log (a)}{(1-a)^{2}} \tag{2.e}$$

The Ornstein-Uhlenbeck process was calibrated using the entire history of the M-CPI rate of inflation up to the year of the final observed diagonal of the SAIF triangle, thus comprising the years 1935 through 2005. The following prior distributions were employed, all of which are highly uninformed:

$$a \sim \operatorname{Beta}(1,1) \tag{2.f}$$

$$b \sim \mathrm{N}\left(0,10^{-6}\right) \tag{2.g}$$

$$\tau_{\varepsilon} \sim \mathrm{Ga}(0.001,0.001) \tag{2.h}$$

After calibrating the Ornstein-Uhlenbeck process to historical data, it is a straightforward procedure to simulate future values for the rate of inflation, which deliver the expert priors for future calendar year effects of the SAIF medical triangle.

2.3. Exposure growth

The rate of exposure growth is modeled as a draw from a normal distribution that accommodates heteroskedasticity by allowing the logarithm of the squared scale parameter to vary with the exposure (e.g., accident) year, similar to the logarithm of the squared scale parameter of the Student’s t distribution. This time variation is estimated using a second-order random walk smoother. The need for a time-varying scale parameter arises from the comparatively high rates of change in volume in the early years of the SCF Arizona and SAIF Oregon triangles. The equations for the rate of exposure growth read

$$\eta_{i} \sim \mathrm{N}\left(\mu_{\eta}, \tau_{\eta}\right) \tag{3.a}$$

$$\mu_{\eta} \sim \text{N}(0,1) \tag{3.b}$$

$$\tau_{\eta, j}=1 / \sigma_{\eta, j}^{2} \tag{3.c}$$

$$\small \log \left(\sigma_{\eta, j}^{2}\right) \sim \mathrm{N}\left(2 \cdot \log \left(\sigma_{\eta, j-1}^{2}\right)-\log \left(\sigma_{\eta, j-2}^{2}\right), v_{\eta}\right), \\ j \geq 4 \tag{3.d}$$

$$v_{\eta} \sim \mathrm{Ga}(5,0.5) \tag{3.e}$$

$$\sigma_{n, j} \sim \mathrm{U}(0,10), j=2,3 . \tag{3.f}$$

Equation (3a) defines the normal distribution, the expected value of which has a highly uninformed prior. The equations for the scale parameter (Equations 3c-f) are equivalent to those contributing to the Student’s t likelihood (Equations 1d-g).

2.4. Structural break

There are two versions of the model, “standard” and “break.” The standard version assumes that the consumption path of services (i.e., the trajectory of exposure-adjusted and calendar year effect-adjusted log incremental payments) is invariant to exposure time (e.g., the accident year), thus resulting in a trajectory S₁, S₂ . . . S_K that applies uniformly to the entire triangle. The break version of the model allows for a (single) structural break in this consumption path, thus differentiating between a pre-structural break trajectory R_1,1, R_1,2 . . . R_1,K and a post-structural break trajectory R_2,1, R_2,2 . . . R_2,L, where L K is the final exposure year of the interval in which the structural break is allowed to occur. Both trajectories are estimated as linear splines using Reversible Jump MCMC, as discussed in Section 2.5.

The model estimates the probability distribution of the change point in the consumption path within a provided interval of exposure years. Such an interval may be set at one year or may comprise multiple years. Within this interval, the consumption paths of the exposure years are mixtures of the pre-break path and the post-break path, except in the extreme case that the model assigns all probability mass to a single year, in which case each year belongs strictly to the pre-break regime or the post-break regime. The change point (C) has a scaled, shifted, and truncated beta prior:

$$\tilde{C} \sim \operatorname{Beta}(2,2) \tag{4.a}$$

$$C=\text { floor }(\tilde{C} \cdot((L+1)-M))+M, \tag{4.b}$$

where M (M ≤ L) is the first year in the interval of exposure years during which the structural break is allowed to occur, L is the last year, and C represents the first year in the new regime; the operator floor(⋅) rounds down to the nearest integer. The density of the beta distribution (Equation 4a) is hump-shaped, thus allocating more probability mass to the center of the interval than to the edges. If the break interval comprises two exposure years, the prior accords each year the same probability of hosting the change point.

In the break version of the model, the consumption path is no longer uniform across the cells of the triangle. Thus, Equation (1b) has to be replaced by

$$\mu_{i, j}=\sum_{k=2}^{i+j} \kappa_{k}+\sum_{k=1}^{i} \eta_{k}+V_{i, j} \tag{4.c}$$

$$V_{i, j}=R_{1, j} \cdot \mathrm{I}(C>\mathrm{i})+R_{2, j} \cdot \mathrm{I}(C \leq \mathrm{i}), \tag{4.d}$$

where I(⋅) is an indicator function that is equal to unity if the statement in the parenthetical is true, and zero otherwise.

2.5. Consumption path

The consumption path is modeled as a linear spline. A linear spline is a piece-wise linear trajectory that changes slopes in locations known as knots. Gimenez et al. (2008) discuss the implementation of a linear spline using MCMC; this spline has a fixed number of knots, but the locations of these knots are endogenous. Because the location of a knot has a posterior distribution instead of just a point mass in a single location, the posterior of a linear spline is not necessarily piece-wise linear but may display a fair degree of smoothness, depending on the number of knots and the posterior uncertainty about their locations. In what follows, a linear spline is implemented where both the number of knots and the locations of these knots are endogenous. The variation in the number of knots implies variation in the dimension of the parameters that need updating. RJMCMC (Reversible Jump MCMC) is an MCMC sampling technique that accommodates such variation in dimension. The implementation of RJMCMC discussed below draws on Lunn, Best, and Whittaker (2008).

2.5.1. Consumption path without structural break

In the version of the model that does not allow for a structural break, the linear spline specification of the consumption path $S_{1 \ldots K}$ reads

$$S_{j}=\beta_{0}+\sum_{g=1}^{\xi} \beta_{g} \cdot \left(j-\theta_{g}\right) \cdot \mathrm{I}\left(j-\theta_{g}>0\right) \tag{5.a}$$

$$\beta_{g} \sim \mathrm{N}(0,0.0001), g=0 \ldots \xi \tag{5.b}$$

$$\vartheta_{1}=1 \tag{5.c}$$

$$\frac{\vartheta_{g}-1}{(b-1)} \sim \operatorname{Beta}(1,2), g=2 \ldots \xi \tag{5.d}$$

$$\theta_{1 . . \xi}=\operatorname{Sort}\left(\vartheta_{1 . .5}\right) \tag{5.e}$$

$$\xi \sim \operatorname{Cat}\left(p_{1} \ldots p_{\lambda}\right), p_{i}=1 / \lambda, i=1 \ldots \lambda \tag{5.f}$$

Equation (5a) describes the linear spline, which consists of an intercept and a series of parameters that indicate the slope of the linear trajectory; the summation index $\xi$ represents the number of knots. In the first development year, the level of consumption $\left(S_1\right)$ equals the intercept $\left(\beta_0\right)$. For $j \leq \theta_2$ (assuming $\xi \geq 2$, otherwise the following equation is unrestricted), $S_j=\beta_0+(j-1) \cdot \beta_1$. The first change in slope (assuming $\xi \geq 2$ ) occurs at knot location $\theta_2$, to the right of which (and up to a possible right-neighboring knot) consumption equals $S_j=\beta_0+(j-1) \cdot \beta_1+\left(j-\theta_2\right) \cdot \beta_2$. To the right of $\theta_3$ (and up to $\theta_4$, assuming $\xi \geq 4$ ), consumption equals $S_j=\beta_0+(j-1) \cdot \beta_1+\left(j-\theta_2\right) \cdot \beta_2+$ $\left(j-\theta_3\right) \cdot \beta_3$, and so on. Equation (5e) sorts the locations of the knots (from smallest to largest development year), thereby allowing for easier interpretation of Equation (5a); the parameter $\vartheta$ indicates the location of a knot before sorting. The prior distributions for the parameters $\beta_{0 \ldots \xi}$ are displayed in Equation (5b).

Equation (5c) sets the location of the first knot to development year 1. Equation (5d) defines the prior distribution for the locations of all other knots as a beta distribution with proper support; the support is defined such that these knots are located in development years between and exclusive of the first development year and the second to last development year that exhibits a positive incremental payment. Equation (5f) affords equal probability to each number of knots (ξ) in the set of integers greater than or equal to unity but less than or equal to λ, where K is the number of columns of the triangle (as mentioned).

2.5.2. Reversible jump Markov chain Monte Carlo implementation

Figure 2 displays the stochastic dependence graph for the parameters updated by means of the RJMCMC algorithm. Solid lines indicate stochastic dependence, whereas dashed lines indicate deterministic dependence. The set {UL} corresponds to the set of indices that pertains to the upper left (UL) portion of the triangle (exclusive of missing or non-positive values, since these values have no bearing on the likelihood). Parameters outside the oval, some of which are shown in Figure 2 to build intuition, are updated (in turn) via more traditional methods, such as slice sampling or Gibbs sampling. The performed RJMCMC update is conditional on the current values of the parameters outside the oval. The nodes with broken outlines need not be included in the RJMCMC update step, since their dimensions are constant; yet, including these nodes greatly aids convergence due to their high correlation with $S_{1 \ldots K}$.

Figure 2.Stochastic dependence graph for RJMCMC

The RJMCMC update process is largely based on Lunn, Best, and Whittaker (2008) and can be organized into four steps: (1) propose a move type, (2) propose a new knot vector, (3) accept or reject the proposed move, and (4) update the remaining parameters inside the oval of Figure 2 via Gibbs sampling. Details of these four steps are presented below.

2.5.2.1. Proposing a move type

At random, a new move type is proposed. Possible moves include birth, death, shuffle, and no change. A birth is the addition of a random number of knots while leaving all current knots in place; a death is the removal of a random number of knots; a shuffle leaves the number of knots unchanged but rearranges a random number of knots; and no change leaves the current number of knots and their positions unchanged.

The type of the proposed move must respect the lower and upper bounds on the number of knots. For instance, if the current number of knots is unity (the lower bound), no death is allowed to occur; similarly, if the current number of knots is equal to λ (the upper bound), no birth is allowed to occur. In order to incorporate these constraints, the probability of proposing a given move type is made a function of these bounds and the current number of knots.

If the proposed move is of type no change, $\xi^{(t+1)}$ and $\theta^{(t+1)}$ are set to $\xi^{(t)}$ and $\theta^{(t)}$, and the process proceeds immediately to step 4. (The superscript (t) refers to the iteration number of the Markov chain.) Otherwise, the number of knots to alter equals unity plus an integer proposed from a Poisson distribution that is truncated so as to respect the mentioned bounds on the number of knots.

2.5.2.2. Proposing a new knot parameter vector

Let the number of knots to alter be denoted $\psi$. If the proposed move is of type birth, then $\xi^{\prime}$ (the proposed value) is set to $\xi^{(t)}+\psi$. Then, for each of the $\psi$ new knots that are to be generated, a new knot is proposed at random from either a uniform distribution or a beta distribution centered on a knot randomly selected from $\theta^{(t)}$, after scaling. If the proposed move is of type death, then $\xi^{\prime}$ is set to $\xi^{(t)}-\psi$, and $\psi$ randomly selected knots are eliminated, in turn, from $\theta^{(t)}$. If the proposed move is of type shuffle, then $\xi^{\prime}$ is set to $\xi^{(t)}$, and $\psi$ randomly selected knots are eliminated from $\theta^{(t)}$. Then, using the surviving $\xi^{(t)}-\psi$ knots as a basis, a birth of $\psi$ knots is performed.

2.5.2.3. Accepting or rejecting the proposed move

In what follows, the acceptance probability ρ = min(ρ′,1) of the newly proposed number of knots (ξ′) and knot vector (θ′) is calculated. As in Lunn, Best, and Whittaker (2008), the acceptance probability is calculated solely based on ξ and θ. This is accomplished by integrating out β, κ, µ_η, and η. Thus, we can write:

$$\begin{array}{l} \rho^{\prime}&=\frac{p\left(y \mid \omega^{(t)}, \tau^{(t)}, \mu_{\mathrm{k}}^{(t)}, \tau_{\mathrm{k}}^{(t)}, \tau_{\eta}^{(t)}, \theta^{\prime}, \xi^{\prime}\right) \times p\left(\xi^{\prime}\right) \times p\left(\theta^{\prime} \mid \xi^{\prime}\right)}{p\left(y \mid \omega^{(t)}, \tau^{(t)}, \mu_{\mathrm{k}}^{(t)}, \tau_{\mathrm{k}}^{(t)}, \tau_{\eta}^{(t)}, \theta^{(t)}, \xi^{(t)}\right)} \\ &\quad \times p\left(\xi^{(t)}\right) \times p\left(\theta^{(t)} \mid \xi^{(t)}\right) \\ &\quad \times \frac{\operatorname{Prop}\left(\left\{\theta^{\prime}, \xi^{\prime}\right\} \rightarrow\left\{\theta^{(t)}, \xi^{(t)}\right\}\right)}{\operatorname{Prop}\left(\left\{\theta^{(t)}, \xi^{(t)}\right\} \rightarrow\left\{\theta^{\prime}, \xi^{\prime}\right\}\right)} \\ \end{array} \tag{5.g}$$

where $\operatorname{Prop}\left(\left\{\theta^{(t)}, \xi^{(t)}\right\} \rightarrow\left\{\theta^{\prime}, \xi^{\prime}\right)\right.$ is the probability of proposing $\left\{\theta^{\prime}, \xi^{\prime}\right\}$, given the current state $\left\{\theta^{(t)}, \xi^{(t)}\right\}$, and $\operatorname{Prop}\left(\left\{\theta^{\prime}, \xi^{\prime}\right\} \rightarrow\left\{\theta^{(t)}, \xi^{(t)}\right\}\right)$ is the probability of proposing the reverse move. Then, with probability $\rho,\left\{\xi^{(t+1)}, \theta^{(t+1)}\right\}$ is set to $\left\{\xi^{\prime}, \theta^{\prime}\right\}$; otherwise, $\left\{\xi^{(t+1)}\right.$, $\left.\theta^{(t+1)}\right\}$ is set to $\left\{\xi^{(t)}, \theta^{(t)}\right\}$.

2.5.2.4. Updating the remaining parameters

Let $\Theta=\left\{\beta_{0 \ldots \xi}, \kappa_{3 \ldots K+1}, \mu_\eta, \eta_{2 \ldots K}\right\}$. All elements of the vector $\Theta$ have a (conditional) normal prior. Furthermore, $y_{i, j} \sim \mathrm{N}\left(\mu_{i, j}, \omega_{i, j} \cdot \tau_j\right)$, where $\mu_{i, j}$ is a linear combination of the elements of $\Theta$. As a result, the (conditional) posterior of $\Theta$ is multivariate normal, which makes the Gibbs sampler available. The algorithm for sampling from $\Theta$ rests on a section of JAGS source code that has been modified to accommodate the variability in the dimension of $\Theta$. The value $\Theta^{(t+1)}$ is drawn directly from the (conditional) posterior of $\Theta$.

2.6. Consumption path with structural break

The break version of the model allows for a structural break in the consumption path across exposure years. A structural break implies different trajectories of log incremental payments (adjusted for exposure growth and calendar year effects) before and after the interval of exposure years in which the structural break is allowed to occur; as discussed, the consumption paths that apply to the individual exposure years within the interval are mixtures of the two the trajectories. Whereas in the standard model this trajectory is denoted $S_{1 \ldots K}$, in the break model, the pre-break trajectory reads $R_{1,1 \ldots K}$ and the post-break trajectory reads $R_{2,1 \ldots K-M+1}$. In accordance to the standard model, the consumption paths are modeled as linear splines:

$$\begin{align} R_{i, j}&=\beta_{i, 0}+\sum_{g=1}^{\xi} \beta_{i, s} \cdot\left(j-\theta_{i, 8}\right) \cdot \mathrm{I}\left(j-\theta_{i, s}>0\right), \\ i&=1,2 \end{align} \tag{5.h}$$

$$\beta_{i, g} \sim \mathrm{N}(0,0.0001), g=0 \ldots \xi_{i}, i=1,2 \tag{5.i}$$

$$\vartheta_{i, 1}=1, i=1,2 \tag{5.j}$$

$$\small{ \frac{\vartheta_{i, s}-1}{\left(b_{i}-1\right)} \sim \operatorname{Beta}\left(\gamma_{i, 1}, \gamma_{i, 2}\right), g=2 \ldots \xi_{i}, i=1,2) \tag{5.k} }$$

$$\vartheta_{i, 1 \ldots \xi_{i}}=\operatorname{Sort}\left(\theta_{i, 1 \ldots \xi_{i}}\right), i=1,2 \tag{5.l}$$

$$\begin{align} \xi_{i} &\sim \operatorname{Cat}\left(p_{1} \ldots p_{\lambda_{i}}\right), \\ p_{j}&=1 / \lambda_{i}, j=1 \ldots \lambda_{i}, i=1,2 .(5 \mathrm{~m}) \end{align} \tag{5.m}$$

All applicable interpretations of variables and parameters are carried forward from Section 2.5.1. Equation (5k) states that the priors for $\vartheta_{1,2 \ldots \xi_1}$ and $\vartheta_{1,2 \ldots \xi_2}$ are beta distributions with proper support. The support is defined in such a way that $\vartheta_{1,2 . . \xi_1}$ are all greater than unity but less than the second to last development year that exhibits a positive incremental payment, and $\vartheta_{2,2, \ldots \xi_2}$ are all greater than unity but less than $K-M$, where $K$ is the number of columns of the triangle (as mentioned) and $M$ is the first year in the interval of exposure years during which the structural break is allowed to occur (as discussed). The parameters for the beta distributions are defined as $\left\{\gamma_{1,1}, \gamma_{1,2}\right\}=\{1,2\}$ and $\left\{\gamma_{2,1}, \gamma_{2,2}\right\}=\{1,1\}$. Equation (5m) affords equal probability to each number of knots $\left(\xi_i\right)$ in the set of integers greater than or equal to unity but less than or equal to $\lambda_i$, where $\lambda_1=$ floor $(K / 2+1)$ and $\lambda_2=\operatorname{floor}((K-M+1) / 2+1)$

2.6.1. Reversible jump Markov chain Monte Carlo implementation

Figure 3 displays the stochastic dependence graph for the parameters updated by means of the RJMCMC algorithm. The meaning of solid and dashed lines is consistent with Figure 2, as is the meaning of the oval and nodes with broken outlines.

Figure 3.Stochastic dependence graph for RJMCMC (break)

The RJMCMC update step for the break case is again largely based on Lunn, Best, and Whittaker (2008) and merely an extension of the algorithm discussed in the context of the standard model. The proposal $\left\{\xi_1^{\prime}, \theta_1^{\prime}\right\}$ is made using the algorithm discussed in Sections 2.5.2.1 and 2.5.2.2, while keeping the current values $\xi_2^{(t)}$ and $\theta_2^{(t)}$ fixed. The acceptance probability for $\left\{\xi_1^{\prime}, \theta_1^{\prime}\right\}$ is calculated as $\rho_1=\min \left(\rho_1^{\prime}, 1\right)$, where

$$\begin{array}{l} &\quad p\left(y \mid \omega^{(i)}, \tau^{(t)}, \mu_{k}^{(t)}, \tau_{k}^{(t)}, \tau_{n}^{(t)}, \theta_{1}^{\prime}, \xi_{1}^{\prime}, \theta_{2}^{(i)}, \xi_{2}^{(t)}\right) \\ p_{1}^{\prime}&=\frac{\times p\left(\xi_{1}^{\prime}\right) \times p\left(\theta_{1}^{\prime} \mid \xi_{1}^{\prime}\right)}{p\left(y \mid \omega^{(t)}, \mathbf{\tau}^{(t)}, \mu_{\mathrm{k}}^{(t)}, \mathbf{\tau}_{k}^{(t)}, \mathbf{\tau}_{n}^{(t)}, \theta_{1}^{\prime}, \xi_{1}^{\prime}, \theta_{2}^{(t)}, \xi_{2}^{(t)}\right)} \\ &\quad \times p\left(\xi_{1}^{(t)}\right) \times p\left(\theta_{1}^{(t)} \mid \xi_{1}^{(t)}\right) \\ &\quad \times \frac{\operatorname{Prop}\left(\left\{\theta_{1}^{\prime}, \xi_{1}^{\prime}\right\} \rightarrow\left\{\theta_{1}^{(t)}, \xi_{1}^{(t)}\right\}\right)}{\operatorname{Prop}\left(\left\{\theta_{1}^{(t)}, \xi_{1}^{(t)}\right\} \rightarrow\left\{\theta_{1}^{\prime}, \xi_{1}^{\prime}\right\}\right)} \\ \end{array} \tag{5.n}$$

The proposal $\left\{\xi_2^{\prime}, \theta_2^{\prime}\right\}$ is also made using the algorithm outlined in Sections 2.5.2.1 and 2.5.2.2, but keeping the updated values $\xi_1^{(t+1)}$ and $\theta_1^{(t+1)}$ fixed. The acceptance probability for $\left\{\xi_2^{\prime}, \theta_2^{\prime}\right\}$ is calculated as $\rho_2=\min \left(\rho_2^{\prime}, 1\right)$, where

$$\begin{array}{l} &\quad p\left(y \mid \omega^{(t)}, \tau^{(t)}, \mu_{\mathrm{k}}^{(t)}, \tau_{k}^{(t)}, \tau_{n}^{(t)}, \theta_{1}^{(t+1)}, \xi_{1}^{(t+1)}, \theta_{2}^{\prime}, \xi_{2}^{\prime}\right) \\ p_{2}^{\prime}&=\frac{\times p\left(\xi_{2}^{\prime}\right) \times p\left(\theta_{2}^{\prime} \mid \xi_{2}^{\prime}\right)}{p\left(y \mid \omega^{(t)}, \tau^{(t)}, \mu_{\mathrm{k}}^{(t)}, \tau_{k}^{(t)}, \mathbf{\tau}_{\eta}^{(t)}, \mathbf{\theta}_{1}^{(t+1)}, \mathbf{\xi}_{1}^{(t+1)}, \mathbf{\theta}_{2}^{(t)}, \mathbf{\xi}_{2}^{(t)}\right)} \\ &\quad \times p\left(\xi_{2}^{(t)}\right) \times p\left(\theta_{2}^{(t)} \mid \xi_{2}^{(t)}\right) \\ &\quad \times \frac{\operatorname{Prop}\left(\left\{\theta_{2}^{\prime}, \xi_{2}^{\prime}\right\} \rightarrow\left\{\theta_{2}^{(t)}, \xi_{2}^{(t)}\right\}\right)}{\operatorname{Prop}\left(\left\{\theta_{2}^{(t)}, \xi_{2}^{(t)}\right\} \rightarrow\left\{\theta_{2}^{\prime}, \xi_{2}^{\prime}\right\}\right)} \\ \end{array} \tag{5.o}$$

The (conditional) posterior for the remaining parameters $\Theta=\left\{\beta_{1,0 \ldots \xi_1}, \beta_{2,0 \ldots \xi_2}, \kappa_{3 \ldots K+1}, \mu_\eta, \eta_{2 \ldots K}\right\}$ is again multivariate normal. Hence, as in Section 2.5.2.4, the value $\Theta^{(t+1)}$ is drawn directly from the conditional posterior of $\Theta$.

2.7. Probability of payment

The discussed model of loss development fits to logarithmic incremental payments. Because the logarithm of zero or less is undefined, incremental payments at amounts of zero or at negative amounts are treated as missing values. Generally, the estimated rate of decay responds to declines in incremental payments as development time progresses; yet, the rate of decay is insensitive to declines to zero because amounts of zero register as missing values on the logarithmic scale. As a result, without accounting for this decline in the probability of payment, the model overestimates the ultimate loss. The solution to this problem is to estimate a trajectory of the probability of payment alongside the model and, draw by draw, multiply the estimated incremental payment by the outcome of a Bernoulli process that models the probability that the payment occurs.

The probability of payment is a function of the number of open claims, which is itself a function of development time and exposure. Because the level of exposure cannot necessarily be assumed as time-invariant, possibly due to a positive trend rate of growth in business volume, the probability of payment may vary along the accident year time axis. At the same time, because all observations of payments at zero amounts that are available for estimating the trajectory of the probability of payment are from high maturities, the estimated trajectory may be appropriate only for the early accident years. For more recent accident years, this trajectory may have to be adjusted, using data that is external to the loss triangle, such as mortality tables or legislative information.

In a given triangle, observations of payments at zero amounts may be sparse, if not nonexistent. Further, although the proportion of payments at zero amounts by column may increase as development time progresses, the total number of observations by column invariably converges to unity before the triangle runs out of reported development years. For this reason, when estimating the Bernoulli process, the prior of the parameter of this Bernoulli distribution, which varies by development year, is determined using a parametric specification. The chosen functional form is a Gompertz curve, which affords the trajectory of the Bernoulli parameter a fair amount of flexibility. The Gompertz curve is restricted to values between zero and unity, thereby leaving two parameters to estimate.

The Bernoulli approach to modeling the probability of payment reads

$$u_{i, j} \sim \operatorname{Bern}\left(p_{j}\right) \tag{6.a}$$

$$p_{j}=1-\exp \left(-\varphi_{1} \cdot \exp \left(-\varphi_{2} \cdot j\right)\right) \tag{6.b}$$

$$\log \left(\varphi_{1}\right) \sim N\left(\log \left(\mu_{\varphi_{1}}\right)-\frac{1}{2 \cdot \tau_{\varphi_{1}}}, \tau_{\varphi_{1}}\right) \tag{6.c}$$

$$\log \left(\varphi_{2}\right) \sim N\left(\log \left(\mu_{\varphi_{, 2}}\right)-\frac{1}{2 \cdot \tau_{\varphi_{2}}}, \tau_{\varphi_{2}}\right) \tag{6.d}$$

$$\tau_{\varphi_{i}} \sim \operatorname{Ga}(10,0.1), i=1,2 \tag{6.e}$$

The expected values $\mu_{\varphi, i=1,2}$ of the prior distributions of the Gompertz parameters in Equations $(6 c, d)$ are chosen based on data exploration prior to estimating the model but, at the same time, impose little restriction on the functional form due to the selected gamma priors. For the SAIF Oregon triangle, these expected values read $\mu_Q=\{1500,0.1355\}$; for the SCF Arizona triangle, it is $\mu_{\varphi}=\{25,000,0.17\}$.

3. Empirical evidence

What follows is an application of the outlined loss development model to the SCF Arizona indemnity triangle and the SAIF Oregon PD medical triangle. The model was estimated in R 2.7.1, using JAGS 1.0.3. The estimation made use of three Markov chains with a burn-in of 2,000,000 iterations and a sample of 2,000,000 iterations, thinned by 5, thus resulting in a sample size of 1,200,000. Each Markov chain was assigned to one processor (on a multi-processor computer) to speed up computation. The estimation results of the SCF indemnity triangle and SAIF medical triangle are presented in Section 3.1 and Section 3.2, respectively.

3.1. Empirical evidence for the SCF Arizona indemnity triangle

The SCF Arizona indemnity triangle consists of accident year data and comprises 66 years of development; the accident years range from 1938 through 2003. The original triangle extends back to accident year 1930, but is only sparsely populated in accident years 1930–1937; hence, these early accident years are excluded from the analysis. The first fully populated diagonal of the SCF triangle falls into calendar year 1988. See Appendix A for a schematic exposition of this triangle.

There were no significant legislative reforms in workers compensation in Arizona during the period 1988 through 2003 (which is the time window covered by fully populated diagonals) that would have given rise to a structural break in the consumption path of indemnity benefits. Hence, no change point specification was included in the model. The results for the SCF Arizona indemnity triangle are summarized in Figures 4 through 16.

Figure 4.SCF indemnity: actual (circles) and estimated (lines) logarithmic incremental payments

Figure 4 presents (in logarithmic terms) the incremental payments, as observed (indicated by circles) and predicted (indicated by lines; each line represents one accident year). In this chart, predicted values are reported only for positive log incremental payments, as indicated by breaks in the lines. Note that no single accident year covers more than 23 development years; most accident years comprise no more than 16 development years; hence, the displayed trajectories of incremental payments each apply only to a limited bracket of (contiguous) development years. Observed and incremental payments line up well, which is in part due to the triangle being one of indemnity payments (as opposed to medical payments); the magnitude of indemnity payments is legally stipulated and highly predictable.

Figure 5 displays (in logarithmic form) the estimated consumption path of indemnity services for the cohort of claimants. By definition, the consumption of indemnity services refers to incremental payments that are adjusted for calendar year effects (and, thus, net of cost-of-living adjustments, where applicable). Another difference to Figure 4 is that the level of exposure is set at the first accident year; the trajectories of all other accident years are simply parallel shifts of the one displayed in the chart, thus adding no information on the development pattern. Initially, the consumption of indemnity services increases (from development year 1 to development year 2), before decreasing monotonically thereafter. This subsequent decline is to be expected, since indemnity services do not lend themselves to utilization increases. At high maturities, indemnity payments typically end with death or, where applicable, are reduced by social security offsets. (No social security offsets apply in the state of Arizona during the studied timeframe.) It is thus not surprising that the rates of decay of the consumption of indemnity services, which are displayed in Figure 6, are correlated with the rate of mortality. The estimated values of these decay rates are provided in Appendix B.

Figure 5.SCF indemnity: trajectory of consumption of indemnity services (level of exposure as of first accident year)

Figure 6.SCF indemnity: trajectory of the rate of decay of consumption

Figure 7 shows the estimated probability of there being an incremental payment in a given development year. As the chart demonstrates, for early development years, the probability of observing a non-zero indemnity payment is essentially unity. This probability starts declining in a meaningful way shortly before development year 50. At the final observed development year (year 66), the probability of observing a non-zero incremental payment is still around 30 percent. Clearly, the decline of this probability is driven by the high development years and, hence, by the early accident years. As a consequence, this probability is greatly influenced by the comparatively low volume of the early accident years and by the mortality rates of the claimants that sustained injuries in those early years. All else being equal, the higher the number of claimants of a given accident year is and the higher the longevity of this cohort of claimants is, the higher is the probability of observing a payment in a given development year. Hence, for the purpose of simulating ultimate losses of more recent (or future) accident years, this trajectory of the probability of payment may need to be shifted to the right.

Figure 7.SCF indemnity: probability of there being a nonzero payment

Figure 8 exhibits the rates of exposure growth. These growth rates are highly variable in the early accident years, possibly due to the comparatively small base or incomplete data. Most notable is a nearly 1,500 percent increase in accident year 1942; this spike in the growth rate follows an accident year with only two recorded payments, one of which is comparatively small (and positive) and the other one is negative. In more recent accident years, the rate of exposure growth is relatively small and smooth. As discussed, the rate of exposure growth was modeled using a normal distribution; the square of the scale parameter of this normal distribution was smoothed (in logarithmic terms) across time using a second-order random walk. Figure 9 displays this scale parameter, alongside 90 percent credible intervals, the bounds of which are subject to the degree of smoothing.

Figure 8.SCF indemnity: rate of exposure growth

Figure 9.SCF indemnity: smoothed scale parameter of normal distribution (exposure growth)

Figure 10 displays the calendar year effect. For the time frame studied here, there was no legally stipulated cost-of-living adjustment of indemnity benefits in Arizona. As a result, the expert prior for the calendar year effect is set to zero; in the chart, this prior is labeled official rate of inflation. Prior to calendar year 1981, the posterior median for the calendar year effect differs from this expert prior only by the sampling error, because the corresponding diagonals are not populated (but instead display missing values). The fully populated diagonals of the time window 1988–2003 exhibit calendar year effects between a negative four percent and a positive ten percent. Past the final observed diagonal, the calendar year effect is simulated from the posterior distribution, as estimated from the preceding set of diagonals.

Figure 10.SCF Indemnity: calendar-year effect

Figure 11 presents the prior and posterior distributions for the degrees of freedom of Student’s t distribution, which is employed in the likelihood that models the logarithmic incremental payments. As mentioned, the prior is a χ²₍₈₎ distribution that is truncated at 2 and 50. The posterior distribution displays a mode of around v = 5 degrees of freedom, which attests to very heavy tails.

Figure 11.SCF Indemnity: prior and posterior distributions of the degrees of freedom of Student’s t distribution

Figure 12 shows the estimates of the scale parameter of Student’s $t$ distribution, $\sigma_j$; the logarithm of the square of this scale parameter was smoothed using a second-order random walk. As suggested by the display of the observed logarithmic payments in Figure 4, the variance of the logarithmic incremental payments increases as the cohort of claimants ages; this is because the number of claimants that contribute to the incremental payments decreases. Not surprisingly, Figure 12 reveals a substantial increase in the scale parameter, which grows to about six times its original value (σ₁) as development time progresses.

Figure 12.SCF Indemnity: smoothed scale parameter of Student’s t distribution

Figure 13 plots the observed logarithmic payments (y_i,j) against the estimated quantiles of the predicated logarithmic incremental payments; this is to gauge the ability of the model to replicate the empirical distribution. At each draw from the posterior, the predicted values (_i,j) were sorted; the chart displays the medians of these sorted values (or, equivalently, estimates of quantiles), which were plotted against the sorted observed values. Most importantly, the chart indicates that the heavy tails of the empirical distributions have been well replicated and that there is no discernible unexplained skewness in the data.

Figure 13.SCF indemnity: sorted observed vs. sorted predicted log incremental payments

Figures 14 through 16 are diagnostic charts that display standardized residuals along the three time axes of the triangle (Figure 14: development year; Figure 15: accident year; Figure 16: calendar year). These standardized residuals are defined as follows (see, for instance, Meyer and Yu 2000):

$$\frac{y_{i, j}-\mu_{i, j}}{\sigma_{j} \cdot \sqrt{\frac{\mathbf{v}}{\mathbf{v}-2}}} \tag{7}$$

Figure 14.SCF indemnity: standardized residuals, by development year

Figure 15.SCF indemnity: standardized residuals, by accident year

Figure 16.SCF indemnity: standardized residuals, by calendar year

The standardized residuals were calculated draw by draw; the charts display the means across the sampled values. The median values of the residuals of a given (development, accident, or calendar) year are indicated by pins with solid pinheads. In a well-calibrated model, these median values show no systematic deviation from zero on the horizontal time axis. Further, with proper adjustment for heteroskedasticity, the variance of these standardized residuals must be uniform on the time axis, which appears to be the case also.

3.2. Empirical evidence for the SAIF Oregon PD medical triangle

The SAIF Oregon accident year triangle comprises 80 years of development of the medical component of PD claims; the accident years run from 1926 through 2005. Sherman and Diss (2006) studied an earlier vintage of this triangle.

Appendix A offers a schematic exposition of the SAIF Oregon triangle. Similar to the SCF Arizona triangle, the SAIF triangle has a large upper left-hand side section of missing values; the first fully populated diagonal of the SAIF triangle falls into calendar year 1985. Due to the sparseness of the data, the first nine accident years (1926–1934) are discarded, thereby reducing the triangle to 71 development years. Starting in development year 42, the triangle is characterized by a high and increasing number of observations with zero payments; past development year 67, all payments are equal to zero.

In 1990, the Oregon Legislature, in special session, made significant changes to the existing Workers’ Compensation Law (SB 1197). This legislative reform may have contributed to a structural break in the consumption path of medical services, potentially by quickening the decay of consumption in the early development years. In order to account for such a potential structural break, the break version of the model is applied. Because the legislative reform occurred in midyear, the Bayesian prior for the location of the change point allocates the same 50 percent probability to the years 1990 and 1991. The empirical results for the SAIF Oregon PD medical triangle are summarized in Figures 17 through 28.

Figure 17.SAIF medical: actual (circles) and estimated (lines) logarithmic incremental payments

Figure 17 displays in logarithmic terms the incremental payments, as observed (indicated by circles) and predicted (indicated by lines; each line represents one accident year). Again, in this chart, predicted values are reported only for positive log incremental payments, as indicated by breaks in the lines. In this triangle, accident years comprise at most 21 observations (except for accident year 1984, which consists of 22 observations). The sets of incremental payments of older accident years tend to start farther to the right and, due to the typically lower volumes of those early years, tend to follow lower trajectories. There is a noticeable increase in the variance of the logarithmic incremental payments as development progresses. Also, as development progresses past development year 20, the lines that represent individual accident years are approximately horizontal; for some accident years, the incremental payments increase slightly over time, whereas for others they decrease somewhat.

Figure 18 presents in logarithmic terms the pre-reform consumption path (gray trajectory) and the post-reform consumption path (black trajectory) of medical services (at the exposure level of the first accident year in the triangle). Consumption of medical services is measured by calendar year effect-adjusted incremental payments. Initially (i.e., from development year 1 to development year 2), consumption increases before decreasing continually in all subsequent development years. Following the 1990 reform, initial exposure-adjusted consumption is slightly higher, but decays faster than before the reform. In the chart, the post-reform consumption path ends with the final development year represented in the time interval of the structural break.

Figure 18.SAIF medical: trajectory of consumption of medical services (level of exposure as of first accident year)

Figure 19 exhibits the pre- and post-reform decay rates of the consumption of medical services. Figure 20 offers a detailed view of these growth rates, starting with the third value. The premise of quantifying a structural break in the decay rate is that the post-reform rate approaches the pre-reform rate as development progresses. This is because, as a cohort of claimants ages, the growth rate of consumption (but not necessarily the level of consumption) can be assumed as increasingly unrelated to the reform and instead increasingly related to mortality. Based on this argument, the trajectory of the post-reform decay rate merges with the pre-reform decay rate where the post-reform consumption path ends.

Figure 19.SAIF triangle: trajectories of the rate of decay of the consumption of medical services, pre-reform and post-reform

Figure 20.SAIF triangle: trajectories of the rate of decay of the consumption of medical services, pre-reform and post-reform (detailed view)

Most interestingly, the decay rate displayed in Figures 19 and 20 stabilizes at a negative 3.93 percent in development year 23. This stationarity after about 22 development years offers an important insight: Generalized to smaller medical triangles (e.g., 30 development years), the estimated decay rate at the final development year may serve as the basis for projecting the consumption path at higher maturities. Appendix B displays the numerical values of the estimated pre-reform and post-reform decay rates.

If, in logarithmic terms, the absolute value of the rate of decay is less than the absolute value of the rate at which claimants leave the cohort, then this implies that the remaining claimants are increasing their consumption. Put differently, although the consumption of medical services decreases for the cohort, the individual claimant’s consumption may increase. In such a situation, an increase in the absolute value of the rate of mortality, as is associated with the aging of a cohort of claimants, implies that the surviving claimants are increasing their consumption at an ever higher rate.

Figure 21 displays the probability of there being a payment in a given development year. The disclaimer made in conjunction with the SCF Arizona indemnity triangle regarding the role of this trajectory for simulated ultimate losses of more recent or future accident years applies here, too.

Figure 21.SAIF triangle: probability of there being a nonzero payment

Figure 22 exhibits the rates of exposure growth. As with the SCF Arizona indemnity triangle, these growth rates are highly variable in the early accident years, but comparatively small and smooth in more recent accident years. Most remarkable is the close to 600 percent increase in exposure in the year 1945, which coincided with World War II veterans rejoining the workforce. Several more but smaller spikes in exposure growth were observed in the following two decades. The displayed rates of exposure growth were affected by the general growth of the economy and by legislative changes. Whereas through 1965, workers compensation benefits in Oregon were paid out of the Industrial Accident Fund, effective January 1, 1966, employers could choose to self-insure, insure with private companies, or stay with the Fund, which is now known as SAIF. Figure 23 shows how the scale parameter of the normal distribution that models exposure growth has declined over time, which is consistent with the smoother and, in absolute value terms, smaller observed rates of growth.

Figure 22.SAIF triangle: rate of exposure growth

Figure 23.SAIF triangle: smoothed scale parameter of normal distribution (exposure growth)

Figure 24 displays the calendar year effect. Prior to calendar year 1984, the posterior median for the calendar year effect differs from the rate of M-CPI inflation only by the sampling error; this is because no incremental payments are available for these diagonals and, hence, no calendar year effect can be estimated. The fully populated diagonals of the time window 1985–2005 exhibit calendar year effects that are several percentage points above or below the rate of M-CPI inflation, which is labeled official rate of inflation in the chart; some of these calendar year effects are negative. Future calendar year effects are simulated using the Ornstein-Uhlenbeck process (which has been calibrated to annual data from 1935 through 2005) and the posterior distribution around this effect. The simulated mean-reverting path of M-CPI inflation is shown by a solid black line; this line increases as it approaches the historical average asymptotically because the M-CPI rate of inflation for 2005 is below the historical average (as recorded through 2005). The credible intervals shown in the chart refer to the rate of M-CPI inflation only, not to the total calendar year effect.

Figure 24.SAIF triangle: calendar year effect

Figure 25 presents the prior and posterior distributions for the degrees of freedom of Student’s t distribution. As with the SCF Arizona indemnity triangle, the low degrees of freedom are indicative of heavy tails.

Figure 25.SAIF triangle: prior and posterior distributions of the degrees of freedom of Student’s t distribution

Figure 26 displays the estimates of the scale parameter of Student’s t distribution. As with the SCF Arizona triangle, the variance of incremental payments increases with development time but, contrary to the SAIF triangle, this variance shows a decline at very high maturities, i.e., before the triangle runs out of payments.

Figure 26.SAIF triangle: smoothed scale parameter of Student’s t distribution

Figure 27 plots the observed logarithmic payments (y_i,j) against the estimated quantiles of the predicated logarithmic incremental payments. Here again, the predicted quantiles line up well with the observed values.

Figure 27.SAIF triangle: sorted observed vs. sorted predicted log incremental payments

Figures 28 through 30 are the diagnostic charts that display the standardized residuals along the three time axes of the triangle. As with the SCF Arizona triangle, the charts indicate that the model is well-specified; most importantly, Figure 28 suggests that the model is capable of accounting for heteroskedasticity on the development time axis.

Figure 28.SAIF triangle: standardized residuals, by development year

Figure 29.SAIF triangle: standardized residuals, by accident year

Figure 30.SAIF triangle: standardized residuals, by calendar year

4. Conclusion

It has been shown that the rate of decay in the consumption of medical services assumes a stationary, negative value after about 20 development years. This implies that, as a cohort of claimants of a given exposure year ages and its number of members shrinks, there is an acceleration of consumption of medical services on the part of the surviving claimants that offsets the acceleration in the rate of mortality, thus resulting in a constant rate of decay in consumption for the cohort. The analysis does not support the hypothesis posited by Sherman and Diss (2006) that there is a bulge in the incremental medical payments at high maturities.

The stationarity of the rate of decay of medical consumption displayed in Appendix B lends itself to straightforward simulation of the development pattern of incremental payments. These incremental payments must not only account for the rate of decay in consumption, but also for price inflation; depending on which of the forces is stronger, incremental payments may increase with development time. Future rates of inflation may be simulated using an Ornstein-Uhlenbeck process that has been calibrated to the rate of M-CPI inflation over the entire history of this price index series. As discussed, any systematic difference between the (logarithmic) rate of M-CPI inflation and the unobserved (logarithmic) rate of inflation for medical services in workers compensation factors into the (logarithmic) rate of decay of consumption of medical services and, hence, is of no concern so long as the relation between these two rates of price inflation can be assumed to be stable.

Finally, it has been shown that legislative reforms, such as the 1990 cost containment reforms implemented in Oregon, may accelerate the decay in the consumption of medical services. Such accelerated decay may spill over into the indemnity triangle, where it manifests itself in shorter claim durations.

Acknowledgment

The author is grateful to SCF Arizona and SAIF Oregon for providing the loss triangles. Thanks also to Jon Evans, Rick Sherman, and Harry Shuford for comments. Chris Laws, Manuel de la Guardia, and Jose Ramos provided research assistance. Chris Laws wrote the C++ Reversible Jump MCMC add-on to JAGS.