# Loss Developments via Growth Curves and Stan

Last week I posted a biological example of fitting a non-linear growth curve with Stan/RStan. Today, I want to apply a similar approach to insurance data using ideas by David Clark [1] and James Guszcza [2].

Instead of predicting the growth of dugongs (sea cows), I would like to predict the growth of cumulative insurance loss payments over time, originated from different origin years. Loss payments of younger accident years are just like a new generation of dugongs, they will be small in size initially, grow as they get older, until the losses are fully settled.

Here is my example data set:

```
library(ChainLadder)
round(GenIns/1000)
```

```
## dev
## origin 1 2 3 4 5 6 7 8 9 10
## 1 358 1125 1735 2218 2746 3320 3466 3606 3834 3901
## 2 352 1236 2170 3353 3799 4120 4648 4914 5339 NA
## 3 291 1292 2219 3235 3986 4133 4629 4909 NA NA
## 4 311 1419 2195 3757 4030 4382 4588 NA NA NA
## 5 443 1136 2128 2898 3403 3873 NA NA NA NA
## 6 396 1333 2181 2986 3692 NA NA NA NA NA
## 7 441 1288 2420 3483 NA NA NA NA NA NA
## 8 359 1421 2864 NA NA NA NA NA NA NA
## 9 377 1363 NA NA NA NA NA NA NA NA
## 10 344 NA NA NA NA NA NA NA NA NA
```

```
plot(GenIns/1000, lattice=TRUE,
layout=c(5,2), pch=19, cex=0.5,
main="Cumulative development by origin year")
```

Following the articles cited above I will assume that the growth can be explained by a Weibull curve, with two parameters \(\theta\) (scale) and \(\omega\) (shape): \[ G(dev|\omega, \theta) = 1 - \exp\left(-\left(\frac{dev}{\theta}\right)^\omega\right) \]

Inspired by the classical over-dispersed Poisson GLM in actuarial science, Guszcza [2] assumes a power variance structure for the process error as well. For the purpose of this post I will assume that claims cost follow a Normal distribution, to make a comparison with a the classical least square regression easier. With a prior estimate of the ultimate claims cost the cumulative loss can be modelled as:

\[ \begin{aligned} CL_{AY, dev} & \sim \mathcal{N}(\mu_{dev}, \sigma^2_{dev})\\ \mu_{dev} & = Ult \cdot G(dev|\omega, \theta)\\ \sigma_{dev} & = \sigma \sqrt{\mu_{dev}} \end{aligned} \]

Perhaps, the above formula suggests a hierarchical modelling approach, with different ultimate loss costs by accident year. Indeed, this is the focus of [2] and I will endeavour to reproduce the results with Stan in my next post, but today I will stick to a pooled model that assumes a constant ultimate loss across all accident years, i.e. \(Ult_{AY} = Ult \;\forall\; AY\).

To prepare my analysis I read in the data as a long data frame instead of the matrix structure above. Additionally, I compose another data frame that I will use later to predict payments two years beyond the first 10 years. Furthermore, to align the output with [2] I relabelled the development periods from years to months, so that year 1 becomes month 6, year 2 becomes month 18, etc. The accident years run from 1991 to 2000, while the variable `origin`

maps those years from 1 to 10.

```
# Read data as long table
url <- "https://raw.githubusercontent.com/mages/diesunddas/master/Data/ClarkTriangle.csv"
dat <- read.csv(url)
dat$origin <- dat$AY-min(dat$AY)+1
dat <- dat[order(dat$dev),]
# Add future dev years
nyears <- 12
newdat <- data.frame(
origin=rep(1:10, each=nyears),
AY=rep(sort(unique(dat$AY)), each=nyears),
dev=rep(seq(from=6, to=nyears*12-6, by=12), 10)
)
newdat <- merge(dat, newdat, all=TRUE)
newdat <- newdat[order(newdat$dev),]
```

To get a reference point for Stan I start with a non-linear least square regression:

```
## Pooled regression
start.vals <- c(ult = 5000, omega = 1.4, theta = 45)
library(nlme)
w0 <- gnls(cum ~ ult*(1 - exp(-(dev/theta)^omega)),
weights = varPower(fixed=.5),
data=dat, start = start.vals)
w0
```

```
## Generalized nonlinear least squares fit
## Model: cum ~ ult * (1 - exp(-(dev/theta)^omega))
## Data: dat
## Log-likelihood: -389.072
##
## Coefficients:
## ult omega theta
## 4642.641548 1.349402 40.054344
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 0.5
## Degrees of freedom: 55 total; 52 residual
## Residual standard error: 6.711733
```

The output above doesn’t look unreasonable, apart from the accident years 1991, 1992 and 1998. The output of `gnls`

gives me an opportunity to provide my prior distributions with good starting points. I will use an inverse Gamma as a prior for \(\sigma\), constrained Normals for the parameters \(\theta, \omega\) and \(Ult\) as well. If you have better ideas, then please get in touch.

The Stan code below is very similar to last week. Again, I am interested here in the posterior distributions, hence I add a block to generate quantities from those. Note, Stan comes with a build-in function for the cumulative Weibull distribution `weibull_cdf`

.

```
stanmodel <-
"data {
int<lower=0> N; // total number of rows
real cum[N]; // cumulative paid
real dev[N]; // development period
// Variables for predictions
int<lower=0> new_N; // number of rows of prediction data set
real new_dev[new_N]; // development periods to predict
}
parameters {
real<lower=0> theta; // scale parameter
real<lower=0> omega; // shape parameter
real<lower=0> ult; // ultimate loss
real<lower=0> tau; // model precession
}
transformed parameters {
real sigma; // process error
sigma = 1 / sqrt(tau);
}
model {
real mu[N];
real disp_sigma[N];
// Priors
theta ~ normal(46, 10); // scale parameter
omega ~ normal(1.4, 1); // shape parameter
tau ~ gamma(.0001, .0001);
ult ~ normal(5000, 1000);
// Hyperparameters: Modelled parameters
for (i in 1:N){
mu[i] <- ult * weibull_cdf(dev[i], omega, theta);
disp_sigma[i] = sigma * sqrt(mu[i]);
}
// Likelihood: Modelled data
cum ~ normal(mu, disp_sigma);
}
generated quantities{
real Y_mean[new_N];
real Y_pred[new_N];
for(i in 1:new_N){
// Posterior parameter distribution of the mean
Y_mean[i] = ult * weibull_cdf(new_dev[i], omega, theta);
// Posterior predictive distribution
Y_pred[i] = normal_rng(Y_mean[i], sigma*sqrt(Y_mean[i]));
}
}"
cat(stanmodel, file="Pooled_Model_for_Loss_Reserving.stan")
```

I store the Stan code in a separate text file, which I read into R to compile the model. The compilation takes a little time. The sampling itself is done in a few seconds.

Let’s take a look at the output:

```
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
pooledDso <- stan_model(file = "Pooled_Model_for_Loss_Reserving.stan",
model_name = "GrowthCurve")
fit <- sampling(pooledDso, iter=5500, warmup=500,
thin=2, cores=2, chains=4, seed=1037025002,
data = list(N = nrow(dat),
cum=dat$cum,
dev=dat$dev,
new_N=nrow(newdat),
new_dev=newdat$dev))
```

```
print(fit, pars=c("ult", "omega", "theta", "sigma"),
probs=c(0.5, 0.75, 0.95))
```

```
## Inference for Stan model: GrowthCurve.
## 4 chains, each with iter=5500; warmup=500; thin=2;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
## mean se_mean sd 50% 75% 95% n_eff Rhat
## ult 4758.39 3.10 242.73 4739.08 4900.08 5186.92 6146 1
## omega 1.35 0.00 0.07 1.35 1.40 1.47 6743 1
## theta 41.74 0.05 3.51 41.45 43.75 47.96 5928 1
## sigma 6.83 0.01 0.68 6.77 7.25 8.04 8027 1
##
## Samples were drawn using NUTS(diag_e) at Sun Aug 16 08:58:26 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
```

`rstan::traceplot(fit, pars=c("ult", "omega", "theta", "sigma"))`

```
stan_dens(fit, pars=c("ult", "omega", "theta", "sigma"),
fill="skyblue")
```

The parameters are a little different to the output of `gnls`

, but well within the standard error. From the plots I notice that the samples for the ultimate loss as well as for \(\theta\) are a little skewed to the right. Well, assuming cumulative losses to follow a Normal distribution was a bit a of stretch to start with. Still, the samples seem to converge.

Finally, I can plot the 90% credible intervals of the posterior distributions.

```
Y_mean <- rstan::extract(fit, "Y_mean")
Y_mean_cred <- apply(Y_mean$Y_mean, 2, quantile, c(0.05, 0.95))
Y_mean_mean <- apply(Y_mean$Y_mean, 2, mean)
Y_pred <- rstan::extract(fit, "Y_pred")
Y_pred_cred <- apply(Y_pred$Y_pred, 2, quantile, c(0.05, 0.95))
Y_pred_mean <- apply(Y_pred$Y_pred, 2, mean)
plot(dat$cum ~ dat$dev, col=(dat$AY+1),
xlab="dev", ylab="cum",
xlim=range(newdat$dev),
main="Non-linear Growth Curve")
lines(newdat$dev, Y_mean_mean)
points(newdat$dev, Y_pred_mean, pch=19)
lines(newdat$dev, Y_mean_cred[1,], col=4)
lines(newdat$dev, Y_mean_cred[2,], col=4)
lines(newdat$dev, Y_pred_cred[1,], col=2)
lines(newdat$dev, Y_pred_cred[2,], col=2)
legend(x=70, y=1500, bty="n", lwd=2, lty=c(NA, NA, 1, 1,1),
legend=c("observation", "prediction", "mean prediction",
"90% mean cred. interval", "90% pred. cred. interval"),
col=c(1,1,1,4,2), pch=c(1, 19, NA, NA, NA))
```

The 90% prediction credible interval captures most of the data and although this model might not be suitable for reserving individual accident years, it could provide an initial starting point for further investigations. Additionally, thanks to the Bayesian model I have access to the full distribution, not just point estimators and standard errors.

My next post will continue with this data set and the ideas of James Guszcza by allowing the ultimate loss cost to vary by accident year, treating it as a random effect. Here is a teaser of what the output will look like:

### References

[1] David R. Clark. *LDF Curve-Fitting and Stochastic Reserving: A Maximum Likelihood Approach.* Casualty Actuarial Society, 2003. CAS Fall Forum.

[2] Guszcza, J. *Hierarchical Growth Curve Models for Loss Reserving*, 2008, CAS Fall Forum,
pp. 146–173.

### Session Info

```
R version 3.2.2 (2015-08-14)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.1 (El Capitan)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] rstan_2.8.0 ggplot2_1.0.1 Rcpp_0.12.1
[4] nlme_3.1-122 lattice_0.20-33 ChainLadder_0.2.3
loaded via a namespace (and not attached):
[1] nloptr_1.0.4 plyr_1.8.3 tools_3.2.2
[4] digest_0.6.8 lme4_1.1-10 statmod_1.4.21
[7] gtable_0.1.2 mgcv_1.8-8 Matrix_1.2-2
[10] parallel_3.2.2 biglm_0.9-1 SparseM_1.7
[13] proto_0.3-10 gridExtra_2.0.0 coda_0.18-1
[16] stringr_1.0.0 MatrixModels_0.4-1 stats4_3.2.2
[19] lmtest_0.9-34 grid_3.2.2 nnet_7.3-11
[22] inline_0.3.14 tweedie_2.2.1 cplm_0.7-4
[25] minqa_1.2.4 reshape2_1.4.1 car_2.1-0
[28] actuar_1.1-10 magrittr_1.5 codetools_0.2-14
[31] scales_0.3.0 MASS_7.3-44 splines_3.2.2
[34] systemfit_1.1-18 rsconnect_0.3.79 pbkrtest_0.4-2
[37] colorspace_1.2-6 labeling_0.3 quantreg_5.19
[40] sandwich_2.3-4 stringi_1.0-1 munsell_0.4.2
[43] zoo_1.7-12
```