<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/correlated-log-normal-chain-ladder-model/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://magesblog.com/post/correlated-lognormal-chain-ladder-model/"> R on mages' blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>On 23 November Glenn Meyers gave a fascinating talk about <a rel="nofollow" href="https://github.com/mages/diesunddas/raw/master/BayesianMixer/Bayesian%20Loss%20Reserving%20Glenn%20Meyers.pdf" target="_blank"><em>The Bayesian Revolution in Stochastic Loss Reserving</em></a> at the <a rel="nofollow" href="https://www.meetup.com//Bayesian-Mixer-London/events/243775592/?showComments=true" target="_blank">10th Bayesian Mixer Meetup</a> in London.</p>
<div>
<img alt="Glenn Meyers speaking at the Bayesian Mixer" src="https://i1.wp.com/magesblog.com/img/magesblog/IMG_1782.jpg?w=456&ssl=1" data-recalc-dims="1">
<p><em>Glenn Meyers speaking at the Bayesian Mixer</em></p>
</div>
<p>Glenn worked for many years as a research actuary at <a rel="nofollow" href="https://www.verisk.com/" target="_blank">Verisk/ ISO</a>, he helped to set up the <a rel="nofollow" href="http://www.casact.org/research/index.cfm?fa=loss_reserves_data" target="_blank"><em>CAS Loss Reserve Database</em></a> and published a monograph on <a rel="nofollow" href="https://www.casact.org/pubs/monographs/papers/01-Meyers.PDF" target="_blank"><em>Stochastic loss reserving using Bayesian MCMC models</em></a>.</p>
<p><a rel="nofollow" href="https://github.com/mages/diesunddas/raw/master/BayesianMixer/Bayesian%20Loss%20Reserving%20Glenn%20Meyers.pdf" target="_blank"><img alt="Slides" src="https://i0.wp.com/magesblog.com/img/magesblog/BayesianLossReservingGlennMeyers.png?w=456&ssl=1" data-recalc-dims="1"></a></p>
<p>In this blog post I will go through the <em>Correlated Log-normal Chain-Ladder Model</em> from his presentation. It is discussed in more detailed in his monograph.</p>
<p>Glenn kindly <a rel="nofollow" href="https://raw.githubusercontent.com/mages/diesunddas/master/BayesianMixer/CCL%20for%20London%20Bayesian%20Meetup.R" target="_blank">shared his code as well</a>, which I have used as a basis for this post.</p>
<div>
<h1>Getting the data</h1>
<p>The CAS Loss Reserve Database is an excellent data source to test reserving models. It is hosted on the <a rel="nofollow" href="http://www.casact.org/research/index.cfm?fa=loss_reserves_data" target="_blank">CAS website</a> and contains historical regulatory filings of US insurance companies.</p>
<p>The following code allows me to download the data and extract the information for one company, here company 353, which was the example company Glenn used as well.</p>
<pre><code>library(data.table)
CASdata &lt;- fread(&quot;http://www.casact.org/research/reserve_data/comauto_pos.csv&quot;)

createLossData &lt;- function(CASdata, company_code){
  
  compData &lt;- CASdata[GRCODE==company_code,
                      c(&quot;EarnedPremNet_C&quot;, &quot;AccidentYear&quot;, &quot;DevelopmentLag&quot;, 
                        &quot;IncurLoss_C&quot;, &quot;CumPaidLoss_C&quot;, &quot;BulkLoss_C&quot;)]
  # rename column names
  setnames(compData, names(compData),
           c(&quot;premium&quot;, &quot;accident_year&quot;, &quot;dev&quot;,
             &quot;incurred_loss&quot;, &quot;paid_loss&quot;, &quot;bulk_loss&quot;))
  compData &lt;- compData[, `:=`(
    origin = accident_year - min(accident_year) + 1,
    # origin period starting from 1, id to identify first origin period
    origin1id = ifelse(accident_year == min(accident_year), 0,  1),
    # Reported incurred loss
    # set losses &lt; 0 to 1, as we will take the log later
    reported_incurred_loss=pmax(incurred_loss - bulk_loss, 1))] 
  # add calendar period and sort data by dev and then origin
  compData &lt;- compData[, cal := origin + dev - 1][order(dev, origin)]
  compData &lt;- compData[, `:=`(
    reported_incurred_loss_train = ifelse(cal &lt;= max(origin),
                                     reported_incurred_loss,
                                     NA),
    reported_incurred_loss_test = ifelse(cal > max(origin),
                                         reported_incurred_loss,
                                         NA))]
  traintest &lt;- rbindlist(list(compData[cal &lt;= max(origin)],
                              compData[cal > max(origin)]))
  return(traintest)
}

lossData &lt;- createLossData(CASdata, company_code = 353)</code></pre>
<p>The data shows the historical annual developments of incurred claims for accident years 1988 to 1997. I separated the data into a training and test data set. The training data will have incurred claims reported up to calendar year 1997, while data from 1998 to 2006 will be used to test the model.</p>
<p>Let’s take a look at the data for <em>Commercial Auto</em> of company 353.</p>
<pre><code>devPlot &lt;- function(x, data, company_code, ...){
  library(lattice)
  key &lt;- list(rep=FALSE, columns=2,
              lines=list(col=c(&quot;#00526D&quot;, &quot;#00526D&quot;), 
                         type=c(&quot;p&quot;, &quot;p&quot;), 
                         pch=c(19, 1)),
              text=list(lab=c(&quot;Observation&quot;, &quot;Hold out observation&quot;)))
  xyplot(x, data, t=c(&quot;b&quot;, &quot;b&quot;), pch=c(19, 1), 
         col=c(&quot;#00526D&quot;, &quot;#00526D&quot;),
         layout=c(5,2),
         par.settings = list(strip.background=list(col=&quot;#CBDDE6&quot;)),
         par.strip.text = list(font = 2),
         key = key, 
         as.table=TRUE,
         scales=list(alternating=1),
         ylab = &quot;Reported incurred loss ($k)&quot;,
         xlab=&quot;Development year&quot;,
         main=&quot;Reported incurred loss development by accident year&quot;,
         sub=paste(&quot;Data source: CAS Loss Reserving Database, Comm. Auto, Comp.&quot;,
                   company_code), ...)
}

devPlot(reported_incurred_loss_train + reported_incurred_loss_test 
        ~ dev | factor(accident_year), 
        data=lossData, company_code = 353)</code></pre>
<p><img src="https://i0.wp.com/magesblog.com/post/2017-12-02-correlated-lognormal-chain-ladder-model_files/figure-html/devPlot-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>The data appears to be reasonably well behaved. Although the years converge to different levels of ultimate incurred claims, the shape of the curves look similar for most years. The accident years 1988 and 1993 show a distinctly different pattern.</p>
</div>
<div>
<h1>The correlated log-normal chain-ladder model</h1>
<p>The correlated log-normal chain-ladder model combines ideas of the classical chain-ladder method, the growth curve model and Bayesian inference.</p>
<p><span>\[<br>
\begin{align}<br>
C_{i,j} &amp; \sim \log\mathcal{N}(\mu_{i,j}, \sigma^2_j)\\<br>
L_i &amp; = P_i \cdot \ell\\<br>
\mu_{1,j} &amp; = \log(L_1) + \beta_j \\<br>
\mu_{i,j} &amp; = \log(L_i) + \alpha_i + \beta_j + \rho \cdot \left(\log(C_{i-1,j}) – \mu_{i-1,j}\right) \\<br>
&amp; \mbox{for } i = 2, \dots, M \mbox{ and } j = 1, \dots, N+1-i \\<br>
&amp; \mbox{and } \alpha_1 := 0,\, \beta_N := 0\\<br>
\sigma^2_j &amp; = \sum_{k=j}^N a_k \\<br>
&amp; \mbox{with priors}\\<br>
\alpha_i &amp; \sim \mathcal{N}(0, 10) \\<br>
\beta_i &amp; \sim \mathcal{N}(0, 10) \\<br>
\log(\ell) &amp; \sim \mathcal{N}(0, 1) \\<br>
\rho &amp; \sim \mbox{Beta}(2, 2)\\<br>
a_k &amp; \sim \mbox{Uniform}(0,1)\\<br>
&amp;\mbox{and}\\<br>
C_{i,j} &amp;:=\mbox{cumulative incurred claims of origin } i\mbox{, development }j \\<br>
P_{i} &amp; := \mbox{premium of origin } i\\<br>
\ell &amp; := \mbox{'population' loss ratio across all origin periods}\\<br>
L_i &amp; := \mbox{initial expected ultimate loss for origin } i\\<br>
M &amp; := \mbox{number of origin periods} \\<br>
N &amp; := \mbox{number of development periods}<br>
\end{align}<br>
\]</span></p>
<p>Given that we assume the losses follow a log-normal distribution we have:</p>
<p><span>\[<br>
\begin{align}<br>
\mbox{Median}(C_{i,j}) &amp; = L_i \cdot \exp(\mu_{i,j})\\<br>
\mathbb{E}(C_{i,j}) &amp; = L_i \cdot \exp\left(\mu_{i,j} + \frac{1}{2} \sigma_j^2 \right)<br>
\end{align}<br>
\]</span> Similar to the chain-ladder method we have one parameter for each development period <span>\(\beta_j\)</span> and one parameter for each origin period <span>\(\alpha_i\)</span>.</p>
<p>But unlike the Mack-chain-ladder model, Glenn allows for correlation between origin periods <span>\(i\)</span> and <span>\(i-1\)</span> in his model: <span>\(\rho \cdot \left(\log(C_{i-1,j}) – \mu_{i-1,j}\right)\)</span>.</p>
<p>Interestingly as well, Glenn assumes that the volatility shrinks as claims mature over development periods, modelling <span>\(\sigma_j^2\)</span> as the sum of uniform variances.</p>
<p>This model is certainly not short on parameters. I wonder, if it is over-parametrised? Does it provide good in sample prediction, but perform poorly out of sample?</p>
</div>
<div>
<h1>Stan model</h1>
<p>So let’s put this model into Stan. Again, I follow Glenn’s implementation, but I have changed some of the variable names and added some lines in the generated quantities block to predict incurred claims for the test data period.</p>
<pre><code>library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())</code></pre>
<pre><code>data{
  int &lt;lower=1> len_data; // number of rows with data
  int&lt;lower=1> len_pred; // number of rows to predict
  // indicator of first origin period
  int&lt;lower=0, upper=1> origin1id[len_data + len_pred]; 
  real logprem[len_data + len_pred];
  real logloss[len_data];
  int&lt;lower=1> origin[len_data + len_pred]; // origin period
  int&lt;lower=1> dev[len_data + len_pred]; // development period
}
transformed data{
  int n_origin = max(origin);
  int n_dev = max(dev);
  int len_total = len_data + len_pred;
}
parameters{
  real r_alpha[n_origin - 1];
  real r_beta[n_dev - 1];
  real log_elr;
  real&lt;lower=0, upper=100000> a_ig[n_dev];
  real&lt;lower=0, upper=1> r_rho;
  real logloss_pred[len_pred];
}
transformed parameters{
  real alpha[n_origin];
  real beta[n_dev];
  real sig2[n_dev];
  real sig[n_dev];
  real mu[len_data];
  real mu_pred[len_pred];
  real &lt;lower=-1, upper=1> rho;
  
  rho = -2*r_rho + 1;
  
  for (i in 1:(n_dev - 1)){
    beta[i] = r_beta[i];
  }
  beta[n_dev] = 0;
  alpha[1] = 0;
  
  for (i in 2:n_origin){
    alpha[i] = r_alpha[i-1];
  }
  
  // Create ascending set of sig2
  sig2[n_dev] = gamma_cdf(1/a_ig[n_dev],1,1); // map into [0,1]
  for (i in 1:(n_dev-1)){ 
    sig2[n_dev - i] = sig2[n_dev + 1 - i] + gamma_cdf(1/a_ig[i],1,1);
  }
  for (i in 1:n_dev){ 
    sig[i] = sqrt(sig2[i]);
  }
  // first origin and dev period (top left corner of triangle)
  mu[1] = logprem[1] + log_elr + beta[dev[1]];
 
  for (i in 2:len_data){
    mu[i] = logprem[i] + log_elr + alpha[origin[i]] +  beta[dev[i]] +
      rho*(logloss[i-1] - mu[i-1]) * origin1id[i];
  }
  
  mu_pred[1] =  logprem[(len_data) + 1] + alpha[origin[len_data + 1]] + 
  log_elr + beta[dev[len_data + 1]] +
      rho*(logloss[len_data] - mu[len_data]) * origin1id[len_data + 1];
      
  for (i in 2:len_pred){
    mu_pred[i] = logprem[len_data + i] + alpha[origin[len_data + i]] + 
    log_elr + beta[dev[len_data + i]] +
      rho*(logloss_pred[i-1] - mu_pred[i-1]) * origin1id[len_data + i];
  }
}
model {
  log_elr ~ normal(0, 1);
  r_alpha ~ normal(0, sqrt(10/1.0));
  r_beta ~ normal(0, sqrt(10/1.0));
  a_ig ~ inv_gamma(1,1); // inverse gamma for numerical resaons
  r_rho ~ beta(2,2);
  // model where we have data
  for (i in 1:(len_data)){
    logloss[i] ~ normal(mu[i], sig[dev[i]]);
  }
  // model where data is missing, the prediction period
  for (i in 1:(len_pred)){
    logloss_pred[i] ~ normal(mu_pred[i], sig[dev[len_data + i]]);
  }
}
generated quantities{
  vector[len_data] log_lik;
  vector[len_total] ppc_loss;
  
  for (i in 1:len_data){ 
    log_lik[i] = normal_lpdf(logloss[i] | mu[i], sig[dev[i]]);
  }
  // simulate posterior predicted losses
  for (i in 1:len_data){ 
    ppc_loss[i] = exp(normal_rng(mu[i], sig[dev[i]]));
  }
  for (i in 1:len_pred){ 
    ppc_loss[len_data + i] = exp(normal_rng(mu_pred[i], sig[dev[len_data + i]]));
  }
}</code></pre>
<p>The next code block prepares the data as an input for Stan.</p>
<pre><code>createStanDataList &lt;- function(lossData){
  with(lossData, {
    train_idx &lt;- !is.na(reported_incurred_loss_train)
    list(
      len_data = sum(train_idx),
      len_pred = sum(!train_idx),
      logprem = log(premium),
      logloss = log(reported_incurred_loss_train[train_idx]),
      origin = origin,
      dev = dev,
      origin1id = origin1id)
  })
}
stan_data &lt;- createStanDataList(lossData)</code></pre>
</div>
<div>
<h1>Model run for company 353</h1>
<p>Finally, we can run our Stan model for company 353.</p>
<pre><code>fitCCL353 &lt;- sampling(CCLmodel,
               data=stan_data,
               seed = 1234, iter = 4000,
               control=list(adapt_delta = 0.99,
                            max_treedepth = 10))</code></pre>
<pre><code>## Warning: There were 4824 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded</code></pre>
<pre><code>## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low</code></pre>
<pre><code>## Warning: Examine the pairs() plot to diagnose sampling problems</code></pre>
<p>Stan comes back with a couple of warnings and messages to investigate the model in more detail.</p>
<p>I get pretty much the same output as Glenn (see slides 19, 20):</p>
<pre><code>print(fitCCL353, 
      pars=c(&quot;alpha&quot;, &quot;beta&quot;, &quot;rho&quot;, &quot;log_elr&quot;, &quot;sig&quot;), 
      probs=c(0.025, 0.975))</code></pre>
<pre><code>## Inference for Stan model: 05311aae980b7925c9bb0a91cd33e3a7.
## 4 chains, each with iter=4000; warmup=2000; thin=1; 
## post-warmup draws per chain=2000, total post-warmup draws=8000.
## 
##            mean se_mean   sd  2.5% 97.5% n_eff Rhat
## alpha[1]   0.00    0.00 0.00  0.00  0.00  8000  NaN
## alpha[2]  -0.26    0.00 0.02 -0.29 -0.23  2770 1.00
## alpha[3]   0.11    0.00 0.02  0.06  0.15  1279 1.00
## alpha[4]   0.21    0.00 0.03  0.16  0.26  1314 1.00
## alpha[5]   0.01    0.00 0.03 -0.06  0.07   699 1.00
## alpha[6]  -0.06    0.00 0.04 -0.13  0.03  1694 1.00
## alpha[7]   0.45    0.00 0.06  0.33  0.56  1432 1.00
## alpha[8]   0.03    0.00 0.08 -0.13  0.20  1790 1.00
## alpha[9]   0.15    0.00 0.14 -0.14  0.43  1341 1.00
## alpha[10]  0.18    0.01 0.29 -0.36  0.77  1300 1.01
## beta[1]   -0.59    0.00 0.11 -0.83 -0.38  4196 1.00
## beta[2]   -0.19    0.00 0.07 -0.32 -0.05  2987 1.00
## beta[3]   -0.10    0.00 0.05 -0.20  0.00  1821 1.00
## beta[4]   -0.03    0.00 0.04 -0.10  0.05  1352 1.00
## beta[5]   -0.01    0.00 0.03 -0.07  0.05  1174 1.00
## beta[6]    0.00    0.00 0.03 -0.06  0.06  1169 1.00
## beta[7]    0.00    0.00 0.03 -0.05  0.06  1010 1.00
## beta[8]    0.01    0.00 0.02 -0.05  0.06  1257 1.00
## beta[9]    0.00    0.00 0.02 -0.05  0.05  1662 1.00
## beta[10]   0.00    0.00 0.00  0.00  0.00  8000  NaN
## rho        0.17    0.00 0.21 -0.27  0.55  2523 1.00
## log_elr   -0.39    0.00 0.01 -0.43 -0.37  1207 1.00
## sig[1]     0.26    0.00 0.09  0.14  0.51  3428 1.00
## sig[2]     0.15    0.00 0.04  0.09  0.26  3884 1.00
## sig[3]     0.09    0.00 0.03  0.05  0.16  1856 1.00
## sig[4]     0.06    0.00 0.02  0.03  0.11   467 1.01
## sig[5]     0.04    0.00 0.02  0.02  0.09   328 1.01
## sig[6]     0.04    0.00 0.02  0.02  0.07   290 1.01
## sig[7]     0.03    0.00 0.01  0.01  0.06   291 1.01
## sig[8]     0.02    0.00 0.01  0.01  0.05   289 1.01
## sig[9]     0.02    0.00 0.01  0.01  0.04   294 1.01
## sig[10]    0.01    0.00 0.01  0.00  0.03   330 1.01
## 
## Samples were drawn using NUTS(diag_e) at Sun Dec  3 21:26:56 2017.
## 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).</code></pre>
<p>Some of the parameters don’t appear significant, such as <span>\(\beta_5\)</span> to <span>\(\beta_9\)</span>, but then again they shouldn’t have done any harm here either, as they will have had little impact.</p>
<p>However, I am most interested in comparing the posterior predictive distribution with the training and test data. How does the model perform out of sample?</p>
<p>The next code chunk extracts the 95% credible interval from the posterior predictive simulated losses and plots the output.</p>
<pre><code>createPlotData &lt;- function(stanfit, data, probs=c(0.25, 0.975)){
  
  ppc_loss &lt;- as.matrix(extract(stanfit, &quot;ppc_loss&quot;)$ppc_loss)
  ppc_loss_summary &lt;- cbind(
    mean=apply(ppc_loss, 2, mean),
    t(apply(ppc_loss, 2, quantile, probs=probs))
  )
  colnames(ppc_loss_summary) = c(
    &quot;Y_pred_mean&quot;,
    paste0(&quot;Y_pred_cred&quot;, gsub(&quot;\\.&quot;, &quot;&quot;, probs[1])),
    paste0(&quot;Y_pred_cred&quot;, gsub(&quot;\\.&quot;, &quot;&quot;, probs[2])))
  
  return(cbind(ppc_loss_summary, data))
}</code></pre>
<pre><code>plotDevBananas &lt;- function(x, data, company_code,
                           xlab=&quot;Development year&quot;, 
                           ylab=&quot;Reported incurred loss ($k)&quot;, 
                           main=&quot;Correlated Log-normal Chain Ladder Model&quot;,
                           ...){
  key &lt;- list(
    rep=FALSE, 
    lines=list(col=c(&quot;#00526D&quot;, &quot;#00526D&quot;, &quot;purple&quot;), 
               type=c(&quot;p&quot;, &quot;p&quot;, &quot;l&quot;), 
               pch=c(19, 1, NA)),
    text=list(lab=c(&quot;Observation&quot;, &quot;Hold out observation&quot;, &quot;Mean estimate&quot;)),
    rectangles = list(col=adjustcolor(&quot;yellow&quot;, alpha.f=0.5), border=&quot;grey&quot;),
    text=list(lab=&quot;95% Prediction credible interval&quot;))
  
  xyplot(x, data=data,as.table=TRUE,xlab=xlab, ylab=ylab, main=main,
         sub=paste(&quot;Data source: CAS Loss Reserving Database,&quot;,
                   &quot;Comm. Auto, Comp.&quot;,
                   company_code),
         scales=list(alternating=1), layout=c(5,2), key=key,
         par.settings = list(strip.background=list(col=&quot;#CBDDE6&quot;)),
         par.strip.text = list(font = 2),
         panel=function(x, y){
           n &lt;- length(x)
           divisor &lt;- 5
           cn &lt;- c(1:(n/divisor))
           upper &lt;- y[cn+n/divisor*0]
           lower &lt;- y[cn+n/divisor*1]
           x &lt;- x[cn]
           panel.polygon(c(x, rev(x)), c(upper, rev(lower)),
                         col = adjustcolor(&quot;yellow&quot;, alpha.f = 0.5), 
                         border = &quot;grey&quot;)
           panel.lines(x, y[cn+n/divisor*2], col=&quot;purple&quot;)
           panel.points(x, y[cn+n/divisor*4], lwd=1, col=&quot;#00526D&quot;)
           panel.points(x, y[cn+n/divisor*3], lwd=1, pch=19, col=&quot;#00526D&quot;)
         }, ...)
}</code></pre>
<pre><code>plotData &lt;- createPlotData(fitCCL353, data=lossData) 

plotDevBananas(Y_pred_cred025 + Y_pred_cred0975 + Y_pred_mean + 
  reported_incurred_loss_train + 
  reported_incurred_loss_test ~ dev | factor(accident_year),
  data=plotData, company_code = 353)</code></pre>
<p><img src="https://i2.wp.com/magesblog.com/post/2017-12-02-correlated-lognormal-chain-ladder-model_files/figure-html/plotModelOuput-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>This looks pretty good!</p>
<p>We can see for the older years that the credible interval shrinks as the development years progress. The skewness of the log-normal distribution is clearly visible for the most recent accident years. In all cases the hold out observations are within the 95% prediction credible interval, often very close to the mean, apart from the years 1994 and 1995, where we observe a more unusual pattern in the data.</p>
<p>Finally, here is the distribution of the ‘population’ loss ratio <span>\(\ell\)</span> across all accident years and correlation parameter <span>\(\rho\)</span>.</p>
<pre><code>summary(elr &lt;- exp(extract(fitCCL353, &quot;log_elr&quot;)$log_elr))</code></pre>
<pre><code>##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6033  0.6695  0.6740  0.6739  0.6786  0.7380</code></pre>
<pre><code>summary(rho &lt;- extract(fitCCL353, &quot;rho&quot;)$rho)</code></pre>
<pre><code>##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.86328  0.03551  0.18278  0.17077  0.31764  0.87749</code></pre>
<p>The distribution of <span>\(\rho\)</span> is very interesting. Many of the traditional reserving methods, including the Mack-chain-ladder model assume that the accident years are independent. However, for this data set it is unlikely to be valid. Indeed, the probability of a positive correlation is 80%.</p>
<pre><code>library(MASS); library(latex2exp)
truehist(rho, col=&quot;skyblue&quot;, xlab=TeX(&quot;$\\\\rho$&quot;)); 
abline(v=0, col=2)</code></pre>
<p><img src="https://i2.wp.com/magesblog.com/post/2017-12-02-correlated-lognormal-chain-ladder-model_files/figure-html/rhoplot-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<div>
<h3>Company 833</h3>
<p>Let’s look at another company. Here is the data for company 833.</p>
<pre><code>lossData &lt;- createLossData(CASdata, company_code = 833)
devPlot(reported_incurred_loss_train + reported_incurred_loss_test 
        ~ dev | factor(accident_year), 
        data=lossData, company_code = 833)</code></pre>
<p><img src="https://i2.wp.com/magesblog.com/post/2017-12-02-correlated-lognormal-chain-ladder-model_files/figure-html/comp833data-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>This data is not as well behaved as for company 353. The shape of the curves do vary quite a bit from one accident year to the next. Let’s find out how well our model predicts the hold out sample in this case.</p>
<pre><code>stan_data &lt;- createStanDataList(lossData)</code></pre>
<pre><code>fitCCL833 &lt;- sampling(CCLmodel,
                      data=stan_data,
                      seed = 1234, iter = 4000,
                      control=list(adapt_delta = 0.99,
                                   max_treedepth = 10))</code></pre>
<pre><code>## Warning: There were 28 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup</code></pre>
<pre><code>## Warning: There were 1443 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded</code></pre>
<pre><code>## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low</code></pre>
<pre><code>## Warning: Examine the pairs() plot to diagnose sampling problems</code></pre>
<p>Stan gives messages back, which I shall ignore for the time being and instead look at my ‘banana’ plot again.</p>
<pre><code>plotData &lt;- createPlotData(fitCCL833, lossData) 
plotDevBananas(Y_pred_cred025 + Y_pred_cred0975 + Y_pred_mean + 
                 reported_incurred_loss_train + 
                 reported_incurred_loss_test ~ dev | factor(accident_year),
               data=plotData, company_code = 833)</code></pre>
<p><img src="https://i2.wp.com/magesblog.com/post/2017-12-02-correlated-lognormal-chain-ladder-model_files/figure-html/comp833plot-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>Interesting! Not too bad either.</p>
</div>
<div>
<h3>Company 25275</h3>
<p>Here is another company, whose data show some unusual patterns.</p>
<pre><code>## Other companies
lossData &lt;- createLossData(CASdata, company_code = 25275)
devPlot(reported_incurred_loss_train + reported_incurred_loss_test 
        ~ dev | factor(accident_year), 
        data=lossData, company_code = 25275)</code></pre>
<p><img src="https://i1.wp.com/magesblog.com/post/2017-12-02-correlated-lognormal-chain-ladder-model_files/figure-html/comp25275Data-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>Let’s fit this data with the same model as well.</p>
<pre><code>stan_data &lt;- createStanDataList(lossData)</code></pre>
<pre><code>fitCCL25275 &lt;-  sampling(CCLmodel,
                         data=stan_data,
                         seed = 1234, iter = 4000,
                         control=list(adapt_delta = 0.99,
                                      max_treedepth = 10))</code></pre>
<pre><code>## Warning: There were 132 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup</code></pre>
<pre><code>## Warning: There were 1738 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded</code></pre>
<pre><code>## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low</code></pre>
<pre><code>## Warning: Examine the pairs() plot to diagnose sampling problems</code></pre>
<pre><code>plotData &lt;- createPlotData(fitCCL25275, lossData) 

plotDevBananas(Y_pred_cred025 + Y_pred_cred0975 + Y_pred_mean + 
                 reported_incurred_loss_train + 
                 reported_incurred_loss_test ~ dev | factor(accident_year),
               data=plotData, company_code = 25275,
               ylim=c(0, 500))</code></pre>
<p><img src="https://i0.wp.com/magesblog.com/post/2017-12-02-correlated-lognormal-chain-ladder-model_files/figure-html/comp25275plot-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>Well, what do you think about that?</p>
</div>
</div>
<div>
<h1>References</h1>
<p><a rel="nofollow" href="https://www.casact.org/pubs/monographs/papers/01-Meyers.PDF" target="_blank"><em>Stochastic Loss Reserving Using Bayesian McMc Models</em></a>. Glenn Meyers. CAS monograph series. Number 1. <strong>Casualty Actuarial Society</strong>. 2015</p>
</div>
<div>
<h1>Session Info</h1>
<pre><code>sessionInfo()</code></pre>
<pre><code>## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.1
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## 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] methods   stats     graphics  grDevices utils     datasets  base     
## 
## other attached packages:
## [1] latex2exp_0.4.0      MASS_7.3-47          rstan_2.16.2        
## [4] StanHeaders_2.16.0-1 ggplot2_2.2.1        lattice_0.20-35     
## [7] data.table_1.10.4-3 
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.14     knitr_1.17       magrittr_1.5     munsell_0.4.3   
##  [5] colorspace_1.3-2 rlang_0.1.4      stringr_1.2.0    plyr_1.8.4      
##  [9] tools_3.4.2      parallel_3.4.2   grid_3.4.2       gtable_0.2.0    
## [13] htmltools_0.3.6  yaml_2.1.14      lazyeval_0.2.1   rprojroot_1.2   
## [17] digest_0.6.12    tibble_1.3.4     bookdown_0.5     gridExtra_2.3   
## [21] codetools_0.2-15 inline_0.3.14    evaluate_0.10.1  rmarkdown_1.8   
## [25] blogdown_0.3     stringi_1.1.6    compiler_3.4.2   scales_0.5.0    
## [29] backports_1.1.1  stats4_3.4.2</code></pre>
</div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://magesblog.com/post/correlated-lognormal-chain-ladder-model/"> R on mages' blog</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/iurkui7nJNA" height="1"></a></div>