brms_freq_sev: Bayesian Frequency Severity Model with BRMS

In ChrisWaller26/bayesact: The Bayesian Actuarial Package

Bayesian Frequency Severity Model with BRMS

Description

This function allows you to create a multivariate Bayesian frequency-severity model in which the frequency parameter is adjusted by the survival function of the severity model at the deductible/excess/attachment point.

This allows for both linear and non-linear functions of the variates for both the frequency and severity model components.

Usage

brms_freq_sev(
  freq_formula = NULL,
  sev_formula = NULL,
  freq_family = NULL,
  sev_family = NULL,
  freq_data = NULL,
  sev_data = NULL,
  prior = NULL,
  ded_name = "ded",
  freq_adj_fun = NULL,
  stanvars = NULL,
  mle = FALSE,
  chains = 2,
  iter = 1000,
  warmup = 250,
  refresh = 100,
  sample_prior = "no",
  ded_adj_min = 0,
  backend = "rstan",
  adapt_delta = 0.8,
  max_treedepth = 8,
  control = NULL,
  seed = sample.int(.Machine$integer.max, 1),
  save_pars = NULL,
  ...
)

Arguments

freq_formula	BRMS Formula; Linear/Non-linear formula for frequency model
sev_formula	BRMS Formula; Linear/Non-linear formula for severity model
freq_family	Family; Family for frequency model
sev_family	Family; Family for severity model
freq_data	Data Frame; The data required for the frequency model
sev_data	Data Frame; The data required for the severity model
prior	BRMS Prior; The set of priors for both the frequency and severity models
ded_name	Character; The column name for the deductible/excess/attachment point in the frequency data
freq_adj_fun	Character; The Stan function used to adjust the mean frequency parameter. If NULL, the survival function of the severity model at the deductible will be used.
mle	Boolean; If TRUE, the optimize function is used to create parameter point estimates via Maximum-Likelihood Estimation
ded_adj_min	Numeric; The minimum value the deductible adjustment can be. This can help when some deductibles are very high.
...	Additional accepted BRMS fit parameters

Value

BRMS Fit

Examples

#### Simulate Frequency Data ####

options(stringsAsFactors = FALSE,
        mc.cores = parallel::detectCores())

#' Assuming one rating factor: region.

set.seed(123456)

# Region Names

regions = c("EMEA", "USC")

# Number of frequency samples

freq_n = 5e3

# Defines a function for lambda

freq_lambda = exp(c(EMEA = 0.5, USC = 1))

# Generate samples for ground-up frequency data

freq_data =
  data.frame(
    pol_id =  seq(freq_n),
    ded = runif(freq_n, 1e3, 5e3),
    lim = runif(freq_n, 50e3, 100e3),
    region = sample(regions, freq_n, replace = T)
  ) %>%
  mutate(
    freq_lambda = freq_lambda[region],
    claimcount_fgu =
      rpois(freq_n, freq_lambda)
  )

#### Simulate severity Data ####

mu_vec = c(EMEA = 8, USC = 9)
sigma_vec = exp(c(EMEA = 0, USC = 0.4))

sev_data =
  data.frame(
    ded = rep(freq_data$ded,
              freq_data$claimcount_fgu),
    lim = rep(freq_data$lim,
              freq_data$claimcount_fgu),
    region = rep(freq_data$region,
                 freq_data$claimcount_fgu)
  ) %>%
  mutate(
    loss_uncapped =
      unlist(
        lapply(
          seq(freq_n),
          function(i){

            rlnorm(freq_data$claimcount_fgu[i],
                   mu_vec[freq_data$region[i]],
                   sigma_vec[freq_data$region[i]]
            )

          }
        )
      )
  ) %>%
  mutate(
    pol_id = rep(seq(freq_n), freq_data$claimcount_fgu)
  ) %>%
  filter(
    loss_uncapped > ded
  ) %>%
  mutate(
    claim_id = row_number(),
    lim_exceed = as.integer(loss_uncapped >= lim),
    loss = pmin(loss_uncapped, lim)
  )

# Frequency data filtered for losses below the deductible

freq_data_net =
  freq_data %>%
  left_join(
    sev_data %>%
      group_by(
        pol_id
      ) %>%
      summarise(
        claimcount = n()
      ) %>%
      ungroup(),
    by = "pol_id"
  ) %>%
  mutate(
    claimcount = coalesce(claimcount, 0)
  )

#### Run Model ####

mv_model_fit =
  brms_freq_sev(

    freq_formula =
      bf(claimcount ~ 1 + region),

    sev_formula =
      bf(loss | trunc(lb = ded) + cens(lim_exceed) ~
           1 + region,
         sigma ~ 1 + region
      ),

    freq_family = poisson(),
    sev_family = lognormal(),

    freq_data = freq_data_net,
    sev_data = sev_data,

    prior = c(prior(normal(0, 1),
                     class = Intercept,
                     resp = claimcount),

               prior(normal(0, 1),
                     class = b,
                     resp = claimcount),

               prior(normal(8, 1),
                     class = Intercept,
                     resp = loss),

               prior(lognormal(0, 1),
                     class = Intercept,
                     dpar = sigma,
                     resp = loss),

               prior(normal(0, 1),
                     class = b,
                     dpar = sigma,
                     resp = loss)
    ),

    ded_name = "ded",

    backend = "cmdstanr",

    chains = 4,
    iter = 1000,
    warmup = 250,
    refresh = 50,
    adapt_delta = 0.999,
    max_treedepth = 15
  )

#### Results ####

model_post_samples =
  posterior_samples(
    mv_model_fit
  ) %>%
  transmute(
    s1_emea = b_loss_s1_Intercept,
    s1_usc  = b_loss_s1_Intercept +
      b_loss_s1_regionUSC,

    sigma_emea = exp(b_sigma_loss_Intercept),
    sigma_usc  = exp(b_sigma_loss_Intercept +
                       b_sigma_loss_regionUSC),

    f1_emea = b_claimcount_f1_Intercept,
    f1_usc  = b_claimcount_f1_Intercept +
      b_claimcount_f1_regionUSC
  )

model_output =
  model_post_samples %>%
  sapply(
    function(x) c(lower = quantile(x, 0.025),
                  mean  = mean(x),
                  upper = quantile(x, 0.975))
  ) %>%
  as.data.frame() %>%
  bind_rows(
    data.frame(
      s1_emea = 8,
      s1_usc  = 9,

      sigma_emea = exp(0),
      sigma_usc  = exp(0.4),

      f1_emea = 0.5,
      f1_usc  = 1
    )
  )

ChrisWaller26/bayesact documentation built on April 3, 2022, 10:26 p.m.