7 Appendix

The appendix presents the R code used in Sections 4 and 5. The code can be copied and pasted into an R session. At the time of writing R version 4.0.2 (2020-06-22) was used, with brms 2.13.5 and rstan 2.21.2.

library(rstan)
library(brms)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

7.1 R code from section 4

The “GenIns” triangle is part of the ChainLadder package. The triangle is transformed into a long table format, with premiums, incremental paid and loss ratio columns added:

library(ChainLadder)
library(data.table)
data(GenIns)
lossDat <- data.table(
  AY = rep(1991:2000, 10),
  dev = sort(rep(1:10, 10)),
  premium = rep(10000000+400000*0:9, 10),
  cum_loss = as.vector(GenIns),
  incr_loss = as.vector(cum2incr(GenIns))
)[order(AY, dev),
  `:=`(cum_lr = cum_loss/premium,
       incr_lr = incr_loss/premium)]

The next code chunk shows how the loss emergence pattern is modelled using differential equations in Stan. The Stan code is stored as a character string in R, and later passed on into brm.

myFuns <- "
real[] ode_lossemergence(real t, real [] y,  real [] theta, 
                         real [] x_r, int[] x_i){
  real dydt[3];
  real ke = theta[1];
  real dr = theta[2];
  real kp1 = theta[3];
  real kp2 = theta[4];
  
  dydt[1] = pow(ke, dr) * pow(t, dr - 1) * exp(-t * ke)/tgamma(dr)
          - (kp1 + kp2) * y[1];
  dydt[2] = kp2 * (y[1] - y[2]);
  dydt[3] = (kp1 *  y[1] + kp2 * y[2]);
 
  return dydt;
}
real int_lossemergence(real t, real ke, real dr, 
                        real kp1, real kp2){
  real y0[3]; real y[1, 3]; real theta[4];
  
  y0[1] = 0; y0[2] = 0; y0[3] = 0; 
  
  theta[1] = ke; 
  theta[2] = dr;
  theta[3] = kp1;
  theta[4] = kp2;
  
  y = integrate_ode_rk45(ode_lossemergence, 
                        y0, 0, rep_array(t, 1), theta,
                        rep_array(0.0, 0), rep_array(1, 1),
                        0.0001, 0.0001, 500); // tolerances, steps
  return (y[1, 3]);
}
real lossemergence(real t, real devfreq, real ke, real dr, 
                   real kp1, real kp2){
    real out = int_lossemergence(t, ke, dr, kp1, kp2);
    if(t > devfreq){ // paid greater dev period 1
    // incremental paid
     out = out - int_lossemergence(t - devfreq, ke, dr, kp1, kp2);
    }
    return(out);
}
"

The following code defines the hierarchical structure using the formula interface in brms:

frml <- bf(incr_lr ~ eta,
           nlf(eta ~ log(ELR * lossemergence(dev, 1.0, ke, dr, kp1, kp2))),
           nlf(ke ~ exp(oke * 0.5)),
           nlf(dr ~ 1 + 0.1 * exp(odr * 0.5)),
           nlf(kp1 ~ 0.5 * exp(okp1 * 0.5)),
           nlf(kp2 ~ 0.1 * exp(okp2 * 0.5)),
           ELR ~ 1 + (1 | AY), 
           oke  ~ 1 + (1 | AY), odr ~ 1 + (1 | AY), 
           okp1 ~ 1 + (1 | AY), okp2 ~ 1 + (1 | AY),
           nl = TRUE)

7.1.1 Multilevel effects with narrow priors

Model run with narrow priors for the multilevel effects:

mypriors <- c(prior(inv_gamma(4, 2), nlpar = "ELR", lb=0),
              prior(normal(0, 1), nlpar = "oke"),
              prior(normal(0, 1), nlpar = "odr"),
              prior(normal(0, 1), nlpar = "okp1"),
              prior(normal(0, 1), nlpar = "okp2"),
              prior(student_t(10, 0, 0.1), class = "sd", nlpar = "ELR"),
              prior(student_t(10, 0, 0.1), class = "sd", nlpar = "oke"),
              prior(student_t(10, 0, 0.1), class = "sd", nlpar = "odr"),
              prior(student_t(10, 0, 0.1), class = "sd", nlpar = "okp1"),
              prior(student_t(10, 0, 0.1), class = "sd", nlpar = "okp2"),
              prior(student_t(10, 0, 1), class = "sigma"))
fit_loss <- brm(frml, prior = mypriors,
                data = lossDat, family = lognormal(), seed = 12345,
                stanvars = stanvar(scode = myFuns, block = "functions"),
                file="models/section_4/GenInsIncModelLog")
fit_loss
#>  Family: lognormal 
#>   Links: mu = identity; sigma = identity 
#> Formula: incr_lr ~ eta 
#>          eta ~ log(ELR * lossemergence(dev, 1, ke, dr, kp1, kp2))
#>          ke ~ exp(oke * 0.5)
#>          dr ~ 1 + 0.1 * exp(odr * 0.5)
#>          kp1 ~ 0.5 * exp(okp1 * 0.5)
#>          kp2 ~ 0.1 * exp(okp2 * 0.5)
#>          ELR ~ 1 + (1 | AY)
#>          oke ~ 1 + (1 | AY)
#>          odr ~ 1 + (1 | AY)
#>          okp1 ~ 1 + (1 | AY)
#>          okp2 ~ 1 + (1 | AY)
#>    Data: lossDat (Number of observations: 55) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Group-Level Effects: 
#> ~AY (Number of levels: 10) 
#>                    Estimate Est.Error l-95% CI
#> sd(ELR_Intercept)      0.04      0.03     0.00
#> sd(oke_Intercept)      0.08      0.06     0.00
#> sd(odr_Intercept)      0.09      0.07     0.00
#> sd(okp1_Intercept)     0.08      0.06     0.00
#> sd(okp2_Intercept)     0.09      0.07     0.00
#>                    u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(ELR_Intercept)      0.11 1.00     1784     1690
#> sd(oke_Intercept)      0.24 1.00     3352     2226
#> sd(odr_Intercept)      0.26 1.00     4252     2019
#> sd(okp1_Intercept)     0.23 1.00     3068     2082
#> sd(okp2_Intercept)     0.26 1.00     3572     1614
#> 
#> Population-Level Effects: 
#>                Estimate Est.Error l-95% CI u-95% CI
#> ELR_Intercept      0.49      0.03     0.43     0.56
#> oke_Intercept     -0.90      0.54    -1.92     0.13
#> odr_Intercept      0.49      1.04    -1.62     2.47
#> okp1_Intercept    -0.40      0.58    -1.40     0.82
#> okp2_Intercept     0.01      0.97    -1.87     2.00
#>                Rhat Bulk_ESS Tail_ESS
#> ELR_Intercept  1.00     4459     3290
#> oke_Intercept  1.00     3769     2887
#> odr_Intercept  1.00     6278     3278
#> okp1_Intercept 1.00     3746     2407
#> okp2_Intercept 1.00     5210     3193
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat
#> sigma     0.37      0.04     0.30     0.45 1.00
#>       Bulk_ESS Tail_ESS
#> sigma     6018     2946
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Population-level posterior parameters on original scale:

x <- posterior_samples(fit_loss, "^b")
mySummary <- function(x){
  c(Estimate = mean(x), Est.Error = sd(x),
    `l-95% CI` = as.numeric(quantile(x, probs = 0.025)),
    `u-95% CI` = as.numeric(quantile(x, probs = 0.975)))
}
rbind(
  ELR = mySummary(x[, 'b_ELR_Intercept']),
  ke = mySummary(exp(x[, 'b_oke_Intercept'] * 0.5)),
  dr = mySummary(1 + 0.1 * exp(x[, 'b_odr_Intercept'] * 0.5)),
  kp1 = mySummary(0.5 * exp(x[, 'b_okp1_Intercept'] * 0.5)),
  kp2 = mySummary(0.1 * exp(x[, 'b_okp2_Intercept'] * 0.5))
  )
#>     Estimate Est.Error l-95% CI u-95% CI
#> ELR   0.4912   0.03330  0.42824   0.5595
#> ke    0.6595   0.17935  0.38199   1.0696
#> dr    1.1454   0.07655  1.04457   1.3441
#> kp1   0.4281   0.13434  0.24851   0.7521
#> kp2   0.1134   0.05934  0.03931   0.2723

7.1.2 Multilevel effects with wider priors

Model run with wider priors for the multilevel effects:

mypriors2 <- c(prior(inv_gamma(4, 2), nlpar = "ELR", lb=0),
              prior(normal(0, 1), nlpar = "oke"),
              prior(normal(0, 1), nlpar = "odr"),
              prior(normal(0, 1), nlpar = "okp1"),
              prior(normal(0, 1), nlpar = "okp2"),
              prior(student_t(10, 0, 1), class = "sd", nlpar = "ELR"),
              prior(student_t(10, 0, 1), class = "sd", nlpar = "oke"),
              prior(student_t(10, 0, 1), class = "sd", nlpar = "odr"),
              prior(student_t(10, 0, 1), class = "sd", nlpar = "okp1"),
              prior(student_t(10, 0, 1), class = "sd", nlpar = "okp2"),
              prior(student_t(10, 0, 1), class = "sigma"))
fit_loss2 <- brm(frml, prior = mypriors2,
                data = lossDat, family = lognormal(), seed = 12345,
                control = list(adapt_delta = 0.9, max_treedepth=15),
                stanvars = stanvar(scode = myFuns, block = "functions"),
                file="models/section_4/GenInsIncModelLog2")
fit_loss2
#>  Family: lognormal 
#>   Links: mu = identity; sigma = identity 
#> Formula: incr_lr ~ eta 
#>          eta ~ log(ELR * lossemergence(dev, 1, ke, dr, kp1, kp2))
#>          ke ~ exp(oke * 0.5)
#>          dr ~ 1 + 0.1 * exp(odr * 0.5)
#>          kp1 ~ 0.5 * exp(okp1 * 0.5)
#>          kp2 ~ 0.1 * exp(okp2 * 0.5)
#>          ELR ~ 1 + (1 | AY)
#>          oke ~ 1 + (1 | AY)
#>          odr ~ 1 + (1 | AY)
#>          okp1 ~ 1 + (1 | AY)
#>          okp2 ~ 1 + (1 | AY)
#>    Data: lossDat (Number of observations: 55) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Group-Level Effects: 
#> ~AY (Number of levels: 10) 
#>                    Estimate Est.Error l-95% CI
#> sd(ELR_Intercept)      0.05      0.04     0.00
#> sd(oke_Intercept)      0.24      0.20     0.01
#> sd(odr_Intercept)      0.67      0.51     0.03
#> sd(okp1_Intercept)     0.24      0.20     0.01
#> sd(okp2_Intercept)     0.81      0.61     0.04
#>                    u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(ELR_Intercept)      0.13 1.00     1764     1934
#> sd(oke_Intercept)      0.73 1.00     2545     1999
#> sd(odr_Intercept)      1.93 1.00     2642     1788
#> sd(okp1_Intercept)     0.75 1.00     2603     2620
#> sd(okp2_Intercept)     2.26 1.00     3638     2070
#> 
#> Population-Level Effects: 
#>                Estimate Est.Error l-95% CI u-95% CI
#> ELR_Intercept      0.49      0.04     0.42     0.57
#> oke_Intercept     -0.88      0.56    -1.94     0.18
#> odr_Intercept      0.29      1.00    -1.70     2.20
#> okp1_Intercept    -0.44      0.60    -1.50     0.84
#> okp2_Intercept     0.00      0.98    -1.92     2.00
#>                Rhat Bulk_ESS Tail_ESS
#> ELR_Intercept  1.00     3294     2689
#> oke_Intercept  1.00     2795     2949
#> odr_Intercept  1.00     5223     2804
#> okp1_Intercept 1.00     2771     3067
#> okp2_Intercept 1.00     3941     2730
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat
#> sigma     0.37      0.04     0.30     0.45 1.00
#>       Bulk_ESS Tail_ESS
#> sigma     5560     2997
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Population-level posterior parameters on original scale:

x <- posterior_samples(fit_loss2, "^b")
rbind(
  ELR = mySummary(x[, 'b_ELR_Intercept']),
  ke = mySummary(exp(x[, 'b_oke_Intercept'] * 0.5)),
  dr = mySummary(1 + 0.1 * exp(x[, 'b_odr_Intercept'] * 0.5)),
  kp1 = mySummary(0.5 * exp(x[, 'b_okp1_Intercept'] * 0.5)),
  kp2 = mySummary(0.1 * exp(x[, 'b_okp2_Intercept'] * 0.5))
  )
#>     Estimate Est.Error l-95% CI u-95% CI
#> ELR   0.4921   0.03795  0.42340   0.5741
#> ke    0.6698   0.19095  0.37871   1.0928
#> dr    1.1310   0.06812  1.04282   1.3004
#> kp1   0.4213   0.13688  0.23628   0.7600
#> kp2   0.1130   0.05958  0.03821   0.2717

Note that in order to predict the models, the user defined Stan functions have to be exported to R via:

expose_functions(fit_loss, vectorize = TRUE)

7.2 R code from case study in section 5

7.2.1 Data

The data used for the case study is a subset of the wkcomp data set from the raw R package (Fannin (2018)).

library(raw)
data(wkcomp)
library(data.table)
library(tidyverse)
# Convert to tibble, rename cols, add calendar year and loss ratio columns
wkcomp <- wkcomp %>% 
  as_tibble() %>% 
  rename(accident_year = AccidentYear, dev_year = Lag, 
         entity_id = GroupCode) %>% 
  mutate(cal_year = accident_year + dev_year - 1,
         paid_loss_ratio = CumulativePaid/DirectEP,
         os_loss_ratio = (CumulativeIncurred - CumulativePaid)/DirectEP)

# Add incremental paid loss ratio column
wkcomp <- wkcomp %>% 
  group_by(entity_id, accident_year) %>% 
  arrange(dev_year) %>% 
  mutate(incr_paid_loss_ratio = paid_loss_ratio - 
           shift(paid_loss_ratio, n=1, fill=0, 
                 type="lag")) %>% 
  ungroup() %>% 
  arrange(entity_id, accident_year, dev_year)

# Stack paid and os into one column + define train and test
wkcomp2 <- wkcomp %>% 
  transmute(
    entity_id, accident_year, dev_year, cal_year, 
    premium = DirectEP, delta = 1, deltaf = "paid",
    loss_ratio_train = ifelse(cal_year < max(accident_year), 
                              incr_paid_loss_ratio, 
                              NA),
    loss_ratio_test = ifelse(cal_year >= max(accident_year), 
                             incr_paid_loss_ratio, 
                             NA),
    loss_amount_train = ifelse(cal_year < max(accident_year), 
                               CumulativePaid, 
                               NA),
    loss_amount_test = ifelse(cal_year >= max(accident_year), 
                              CumulativePaid, 
                              NA)
  ) %>% 
  bind_rows(
    wkcomp %>% 
      transmute(
        entity_id, accident_year, dev_year, cal_year, 
        premium = DirectEP, delta = 0, deltaf = "os",
        loss_ratio_train = ifelse(cal_year < max(accident_year), 
                                  os_loss_ratio, 
                                  NA),
        loss_ratio_test = ifelse(cal_year >= max(accident_year), 
                                 os_loss_ratio, 
                                 NA),
        loss_amount_train = ifelse(cal_year < max(accident_year), 
                                   CumulativeIncurred - CumulativePaid, 
                                   NA),
        loss_amount_test = ifelse(cal_year >= max(accident_year), 
                                  CumulativeIncurred - CumulativePaid, 
                                  NA)
      )
  )

Filter for company “337”:

dat337 <- wkcomp2 %>% filter(entity_id ==337)

7.2.2 Model 1

myFunsCumPaid <- "
real paid(real t, real ker, real kp, real RLR, real RRF){
 return(
  RLR*RRF/(ker - kp) * (ker *(1 - exp(-kp*t)) - 
  kp*(1 - exp(-ker*t)))
 );
}
real os(real t, real ker, real kp, real RLR){
 return(
  (RLR*ker/(ker - kp) * (exp(-kp*t) - exp(-ker*t)))
 );
}
real claimsprocess(real t, real ker, real kp, 
                   real RLR, real RRF, real delta){
    real out; 
    out = os(t, ker, kp, RLR) * (1 - delta) + 
          paid(t, ker, kp, RLR, RRF) * delta;
    
    return(out);
}
"
frml1 <- bf(loss_amount_train ~ premium * claimsprocess(dev_year, ker, kp,
                                                        RLR, RRF, delta),
            nlf(ker ~ 3 * exp(oker * 0.1)),
            nlf(kp ~ 1 * exp(okp * 0.1)),
            nlf(RLR ~ 0.7 * exp(oRLR * 0.2)),
            nlf(RRF ~ 0.8 * exp(oRRF * 0.1)),
            oRLR ~ 1 + (1 | ID | accident_year),
            oRRF ~ 1 + (1 | ID | accident_year),
            oker ~ 1,  okp ~ 1, 
            sigma ~ 0 + deltaf,
            nl = TRUE)
mypriors1 <- c(prior(normal(0, 1), nlpar = "oRLR"),
               prior(normal(0, 1), nlpar = "oRRF"),
               prior(normal(0, 1), nlpar = "oker"),
               prior(normal(0, 1), nlpar = "okp"),
               prior(student_t(1, 0, 1000), class = "b", 
                     coef="deltafpaid", dpar= "sigma"),
               prior(student_t(1, 0, 1000), class = "b", 
                     coef="deltafos", dpar= "sigma"),
               prior(student_t(10, 0, 0.2), class = "sd", nlpar = "oRLR"),
               prior(student_t(10, 0, 0.1), class = "sd", nlpar = "oRRF"),
               prior(lkj(1), class="cor"))
m1fit <- brm(frml1, data = dat337[!is.na(loss_ratio_train)], 
             family = gaussian(),
             prior = mypriors1,
             stanvars = stanvar(scode = myFunsCumPaid, block = "functions"),
             control = list(adapt_delta = 0.99, max_treedepth=15),
             file="models/section_5/CaseStudy_Model_1",
             seed = 123, iter = 2000, chains = 4)
m1fit
#>  Family: gaussian 
#>   Links: mu = identity; sigma = log 
#> Formula: loss_amount_train ~ premium * claimsprocess(dev_year, ker, kp, RLR, RRF, delta) 
#>          ker ~ 3 * exp(oker * 0.1)
#>          kp ~ 1 * exp(okp * 0.1)
#>          RLR ~ 0.7 * exp(oRLR * 0.2)
#>          RRF ~ 0.8 * exp(oRRF * 0.1)
#>          oRLR ~ 1 + (1 | ID | accident_year)
#>          oRRF ~ 1 + (1 | ID | accident_year)
#>          oker ~ 1
#>          okp ~ 1
#>          sigma ~ 0 + deltaf
#>    Data: dat337[!is.na(loss_ratio_train)] (Number of observations: 90) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Group-Level Effects: 
#> ~accident_year (Number of levels: 9) 
#>                                    Estimate Est.Error
#> sd(oRLR_Intercept)                     0.64      0.15
#> sd(oRRF_Intercept)                     0.74      0.22
#> cor(oRLR_Intercept,oRRF_Intercept)     0.54      0.26
#>                                    l-95% CI u-95% CI
#> sd(oRLR_Intercept)                     0.40     1.00
#> sd(oRRF_Intercept)                     0.33     1.22
#> cor(oRLR_Intercept,oRRF_Intercept)    -0.01     0.97
#>                                    Rhat Bulk_ESS
#> sd(oRLR_Intercept)                 1.00     1779
#> sd(oRRF_Intercept)                 1.00     1763
#> cor(oRLR_Intercept,oRRF_Intercept) 1.00     1300
#>                                    Tail_ESS
#> sd(oRLR_Intercept)                     2221
#> sd(oRRF_Intercept)                     1536
#> cor(oRLR_Intercept,oRRF_Intercept)     1532
#> 
#> Population-Level Effects: 
#>                  Estimate Est.Error l-95% CI u-95% CI
#> oRLR_Intercept       1.44      0.26     0.92     1.95
#> oRRF_Intercept      -1.36      0.42    -2.17    -0.50
#> oker_Intercept      -5.40      0.81    -6.87    -3.71
#> okp_Intercept       -8.49      0.36    -9.17    -7.79
#> sigma_deltafos       8.25      0.16     7.96     8.57
#> sigma_deltafpaid     6.66      0.20     6.31     7.08
#>                  Rhat Bulk_ESS Tail_ESS
#> oRLR_Intercept   1.00     1137     2014
#> oRRF_Intercept   1.00     2518     2684
#> oker_Intercept   1.00     1821     2368
#> okp_Intercept    1.00     2491     2826
#> sigma_deltafos   1.00     1485     2306
#> sigma_deltafpaid 1.00     1825     2109
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

7.2.3 Model 2

myFunsIncrPaid <- "
real paid(real t, real ker, real kp, real RLR, real RRF){
 return(
  RLR*RRF/(ker - kp) * (ker *(1 - exp(-kp*t)) - 
  kp*(1 - exp(-ker*t)))
 );
}
real os(real t, real ker, real kp, real RLR){
 return(
  (RLR*ker/(ker - kp) * (exp(-kp*t) - exp(-ker*t)))
 );
}
real claimsprocess(real t, real devfreq, real ker, real kp, 
                   real RLR, real RRF, real delta){
    real out; 
    out = os(t, ker, kp, RLR) * (1 - delta) + 
          paid(t, ker, kp, RLR, RRF) * delta;
    
    if( (delta > 0) && (t > devfreq) ){ // paid greater dev period 1
    // incremental paid
     out = out - paid(t - devfreq, ker, kp, RLR, RRF)*delta;
    }
    return(out);
}
"
frml2 <- bf(loss_ratio_train ~ eta,
            nlf(eta ~ log(claimsprocess(dev_year, 1.0, ker, kp,
                                        RLR, RRF, delta))),
            nlf(ker ~ 3 * exp(oker * 0.1)),
            nlf(kp ~ 1 * exp(okp * 0.1)),
            nlf(RLR ~ 0.7 * exp(oRLR * 0.2)),
            nlf(RRF ~ 0.8 * exp(oRRF * 0.1)),
            oRLR ~ 1 + (1 | ID | accident_year) + (1 | dev_year),
            oRRF ~ 1 + (1 | ID | accident_year) + (1 | dev_year),
            oker ~ 1  + (1 | accident_year) + (1 | dev_year), 
            okp ~ 1 + (1 | accident_year) + (1 | dev_year), 
            sigma ~ 0 + deltaf, nl = TRUE)
mypriors2 <- c(prior(normal(0, 1), nlpar = "oRLR"),
               prior(normal(0, 1), nlpar = "oRRF"),
               prior(normal(0, 1), nlpar = "oker"),
               prior(normal(0, 1), nlpar = "okp"),
               prior(normal(log(0.2), 0.2), class = "b", 
                     coef="deltafpaid", dpar= "sigma"),
               prior(normal(log(0.2), 0.2), class = "b", 
                     coef="deltafos", dpar= "sigma"),
               prior(student_t(10, 0, 0.3), class = "sd", nlpar = "oker"),
               prior(student_t(10, 0, 0.3), class = "sd", nlpar = "okp"),
               prior(student_t(10, 0, 0.7), class = "sd", nlpar = "oRLR"),
               prior(student_t(10, 0, 0.5), class = "sd", nlpar = "oRRF"),
               prior(lkj(1), class="cor"))
m2fit <- brm(frml2, data = dat337[!is.na(loss_ratio_train)], 
             family = brmsfamily("lognormal", link_sigma = "log"),
             prior = mypriors2,
             stanvars = stanvar(scode = myFunsIncrPaid, block = "functions"),
             control = list(adapt_delta = 0.99, max_treedepth=15),
             file="models/section_5/CaseStudy_Model_2",
             seed = 123, iter = 2000, chains = 4)
m2fit
#>  Family: lognormal 
#>   Links: mu = identity; sigma = log 
#> Formula: loss_ratio_train ~ eta 
#>          eta ~ log(claimsprocess(dev_year, 1, ker, kp, RLR, RRF, delta))
#>          ker ~ 3 * exp(oker * 0.1)
#>          kp ~ 1 * exp(okp * 0.1)
#>          RLR ~ 0.7 * exp(oRLR * 0.2)
#>          RRF ~ 0.8 * exp(oRRF * 0.1)
#>          oRLR ~ 1 + (1 | ID | accident_year) + (1 | dev_year)
#>          oRRF ~ 1 + (1 | ID | accident_year) + (1 | dev_year)
#>          oker ~ 1 + (1 | accident_year) + (1 | dev_year)
#>          okp ~ 1 + (1 | accident_year) + (1 | dev_year)
#>          sigma ~ 0 + deltaf
#>    Data: dat337[!is.na(loss_ratio_train)] (Number of observations: 90) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Group-Level Effects: 
#> ~accident_year (Number of levels: 9) 
#>                                    Estimate Est.Error
#> sd(oRLR_Intercept)                     0.69      0.31
#> sd(oRRF_Intercept)                     0.73      0.47
#> sd(oker_Intercept)                     0.26      0.22
#> sd(okp_Intercept)                      2.23      1.35
#> cor(oRLR_Intercept,oRRF_Intercept)     0.37      0.47
#>                                    l-95% CI u-95% CI
#> sd(oRLR_Intercept)                     0.15     1.37
#> sd(oRRF_Intercept)                     0.04     1.79
#> sd(oker_Intercept)                     0.01     0.80
#> sd(okp_Intercept)                      0.57     5.44
#> cor(oRLR_Intercept,oRRF_Intercept)    -0.69     0.98
#>                                    Rhat Bulk_ESS
#> sd(oRLR_Intercept)                 1.00     1032
#> sd(oRRF_Intercept)                 1.00     1042
#> sd(oker_Intercept)                 1.00     3430
#> sd(okp_Intercept)                  1.00      470
#> cor(oRLR_Intercept,oRRF_Intercept) 1.00     2682
#>                                    Tail_ESS
#> sd(oRLR_Intercept)                     1089
#> sd(oRRF_Intercept)                     1527
#> sd(oker_Intercept)                     1638
#> sd(okp_Intercept)                      1269
#> cor(oRLR_Intercept,oRRF_Intercept)     2614
#> 
#> ~dev_year (Number of levels: 9) 
#>                    Estimate Est.Error l-95% CI
#> sd(oRLR_Intercept)     0.35      0.26     0.01
#> sd(oRRF_Intercept)     1.11      0.49     0.13
#> sd(oker_Intercept)     0.27      0.23     0.01
#> sd(okp_Intercept)      0.45      0.43     0.02
#>                    u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(oRLR_Intercept)     0.98 1.00     1072     1669
#> sd(oRRF_Intercept)     2.12 1.00      949      921
#> sd(oker_Intercept)     0.86 1.00     3410     1546
#> sd(okp_Intercept)      1.32 1.01      552      903
#> 
#> Population-Level Effects: 
#>                  Estimate Est.Error l-95% CI u-95% CI
#> oRLR_Intercept       1.54      0.44     0.70     2.47
#> oRRF_Intercept      -1.45      0.69    -2.78    -0.07
#> oker_Intercept      -1.31      1.08    -3.40     0.84
#> okp_Intercept       -5.07      1.75    -7.71    -1.50
#> sigma_deltafos      -1.78      0.16    -2.13    -1.49
#> sigma_deltafpaid    -1.89      0.16    -2.20    -1.57
#>                  Rhat Bulk_ESS Tail_ESS
#> oRLR_Intercept   1.00     1167     1897
#> oRRF_Intercept   1.00     1689     2552
#> oker_Intercept   1.00     4150     2685
#> okp_Intercept    1.00      499     1205
#> sigma_deltafos   1.00     1109     1860
#> sigma_deltafpaid 1.00     1198     2632
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

7.2.4 Model 3

CycleIndex <- data.table(
  accident_year = 1988:1997,
  RLM = c(1, 1.18, 1.22, 1.05, 1, 0.87, 0.94, 1.34, 1.64, 2.14)^0.6,
  RRM = c(1, 1.05, 1.05, 1.01, 1.0, 0.95, 0.99, 1.1, 1.25, 1.35)^0.6
)
setkey(dat337, accident_year)
setkey(CycleIndex, accident_year)
dat337 <- CycleIndex[dat337]
frml3 <- bf(loss_ratio_train ~ eta,
            nlf(eta ~ log(claimsprocess(dev_year, 1.0, ker, kp,
                                        RLR, RRF, delta))),
            nlf(ker ~ 3 * exp(oker * 0.1)),
            nlf(kp ~ 1 * exp(okp * 0.1)),
            nlf(RLR ~ 0.7 * exp(oRLR * 0.2) * (RLM^lambda1)),
            nlf(RRF ~ 0.8 * exp(oRRF * 0.1) * (RRM^lambda2)),
            oRLR ~ 1 + (1 | ID | accident_year) + (1 | dev_year),
            oRRF ~ 1 + (1 | ID | accident_year) + (1 | dev_year),
            lambda1 ~  1,
            lambda2 ~ 1,
            oker ~ 1  + (1 | accident_year) + (1 | dev_year), 
            okp ~ 1 + (1 | accident_year) + (1 | dev_year), 
            sigma ~ 0 + deltaf, nl = TRUE)
mypriors3 <- c(prior(normal(0, 1), nlpar = "oRLR"),
               prior(normal(0, 1), nlpar = "oRRF"),
               prior(normal(0, 1), nlpar = "oker"),
               prior(normal(0, 1), nlpar = "okp"),
               prior(normal(log(0.2), 0.2), class = "b", 
                     coef="deltafpaid", dpar= "sigma"),
               prior(normal(log(0.2), 0.2), class = "b", 
                     coef="deltafos", dpar= "sigma"),
               prior(student_t(10, 0, 0.3), class = "sd", nlpar = "oker"),
               prior(student_t(10, 0, 0.3), class = "sd", nlpar = "okp"),
               prior(student_t(10, 0, 0.7), class = "sd", nlpar = "oRLR"),
               prior(student_t(10, 0, 0.5), class = "sd", nlpar = "oRRF"),
               prior(normal(1, 0.25), nlpar = "lambda1"), 
               prior(normal(1, 0.25), nlpar = "lambda2"), 
               prior(lkj(1), class="cor"))
m3fit <- brm(frml3, data = dat337[!is.na(loss_ratio_train)], 
             family = brmsfamily("lognormal", link_sigma = "log"),
             prior = mypriors3,
             stanvars = stanvar(scode = myFunsIncrPaid, block = "functions"),
             control = list(adapt_delta = 0.99, max_treedepth=15),
             file="models/section_5/CaseStudy_Model_3",
             seed = 123, iter = 2000, chains = 4)
m3fit
#>  Family: lognormal 
#>   Links: mu = identity; sigma = log 
#> Formula: loss_ratio_train ~ eta 
#>          eta ~ log(claimsprocess(dev_year, 1, ker, kp, RLR, RRF, delta))
#>          ker ~ 3 * exp(oker * 0.1)
#>          kp ~ 1 * exp(okp * 0.1)
#>          RLR ~ 0.7 * exp(oRLR * 0.2) * (RLM^lambda1)
#>          RRF ~ 0.8 * exp(oRRF * 0.1) * (RRM^lambda2)
#>          oRLR ~ 1 + (1 | ID | accident_year) + (1 | dev_year)
#>          oRRF ~ 1 + (1 | ID | accident_year) + (1 | dev_year)
#>          lambda1 ~ 1
#>          lambda2 ~ 1
#>          oker ~ 1 + (1 | accident_year) + (1 | dev_year)
#>          okp ~ 1 + (1 | accident_year) + (1 | dev_year)
#>          sigma ~ 0 + deltaf
#>    Data: dat337[!is.na(loss_ratio_train)] (Number of observations: 90) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Group-Level Effects: 
#> ~accident_year (Number of levels: 9) 
#>                                    Estimate Est.Error
#> sd(oRLR_Intercept)                     0.29      0.21
#> sd(oRRF_Intercept)                     0.80      0.40
#> sd(oker_Intercept)                     0.26      0.21
#> sd(okp_Intercept)                      1.69      1.28
#> cor(oRLR_Intercept,oRRF_Intercept)     0.17      0.52
#>                                    l-95% CI u-95% CI
#> sd(oRLR_Intercept)                     0.01     0.80
#> sd(oRRF_Intercept)                     0.08     1.62
#> sd(oker_Intercept)                     0.01     0.78
#> sd(okp_Intercept)                      0.40     5.02
#> cor(oRLR_Intercept,oRRF_Intercept)    -0.85     0.96
#>                                    Rhat Bulk_ESS
#> sd(oRLR_Intercept)                 1.00     1435
#> sd(oRRF_Intercept)                 1.00     1383
#> sd(oker_Intercept)                 1.00     4237
#> sd(okp_Intercept)                  1.00      392
#> cor(oRLR_Intercept,oRRF_Intercept) 1.00     1747
#>                                    Tail_ESS
#> sd(oRLR_Intercept)                     1803
#> sd(oRRF_Intercept)                     1811
#> sd(oker_Intercept)                     2390
#> sd(okp_Intercept)                      1217
#> cor(oRLR_Intercept,oRRF_Intercept)     2039
#> 
#> ~dev_year (Number of levels: 9) 
#>                    Estimate Est.Error l-95% CI
#> sd(oRLR_Intercept)     0.37      0.27     0.02
#> sd(oRRF_Intercept)     1.16      0.52     0.17
#> sd(oker_Intercept)     0.27      0.23     0.01
#> sd(okp_Intercept)      0.48      0.53     0.02
#>                    u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(oRLR_Intercept)     1.02 1.00     1269     1897
#> sd(oRRF_Intercept)     2.24 1.00     1020      905
#> sd(oker_Intercept)     0.83 1.00     4717     2299
#> sd(okp_Intercept)      1.61 1.00      924      867
#> 
#> Population-Level Effects: 
#>                   Estimate Est.Error l-95% CI u-95% CI
#> oRLR_Intercept        1.33      0.42     0.58     2.25
#> oRRF_Intercept       -1.58      0.70    -2.91    -0.14
#> lambda1_Intercept     1.06      0.21     0.65     1.46
#> lambda2_Intercept     1.02      0.25     0.53     1.50
#> oker_Intercept       -1.35      1.09    -3.46     0.75
#> okp_Intercept        -5.78      1.69    -7.95    -1.85
#> sigma_deltafos       -1.80      0.15    -2.11    -1.52
#> sigma_deltafpaid     -1.91      0.16    -2.22    -1.60
#>                   Rhat Bulk_ESS Tail_ESS
#> oRLR_Intercept    1.00     1434     2624
#> oRRF_Intercept    1.00     2272     2852
#> lambda1_Intercept 1.00     6838     3540
#> lambda2_Intercept 1.00     8468     2802
#> oker_Intercept    1.00     4601     2525
#> okp_Intercept     1.01      450     1590
#> sigma_deltafos    1.00     1567     2758
#> sigma_deltafpaid  1.00     1535     2518
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Note, in order to predict from the models, the user defined Stan functions have to be exported to R via:

expose_functions(m1, vectorize = TRUE)
expose_functions(m2, vectorize = TRUE)
expose_functions(m3, vectorize = TRUE)

7.3 Session info

utils:::print.sessionInfo(session_info, local=FALSE)
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04 LTS
#> 
#> Matrix products: default
#> BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-openmp/libopenblasp-r0.3.8.so
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets 
#> [6] methods   base     
#> 
#> other attached packages:
#>  [1] cowplot_1.0.0        ggridges_0.5.2      
#>  [3] raw_0.1.6            MASS_7.3-51.6       
#>  [5] knitr_1.29           modelr_0.1.8        
#>  [7] forcats_0.5.0        stringr_1.4.0       
#>  [9] dplyr_1.0.0          purrr_0.3.4         
#> [11] readr_1.3.1          tidyr_1.1.0         
#> [13] tibble_3.0.3         tidyverse_1.3.0     
#> [15] tidybayes_2.1.1      latticeExtra_0.6-29 
#> [17] lattice_0.20-41      ChainLadder_0.2.11  
#> [19] data.table_1.13.0    bayesplot_1.7.2     
#> [21] brms_2.13.5          Rcpp_1.0.5          
#> [23] rstan_2.21.2         ggplot2_3.3.2       
#> [25] StanHeaders_2.21.0-6 deSolve_1.28

References

Fannin, Brian A. 2018. Raw: R Actuarial Workshops. https://CRAN.R-project.org/package=raw.