R/DynamicRate.R

In alasdairhunter/HazardLossModel:

#' Function For Calculating Dynamic Rate Parameters For ELTs
#'
#' This function allows you to calculate rate parameters for event loss tables (ELTs).
#' @param elt The data frame or array containing the elt. The first column should contain the annual rate and the second the expected loss.
#' @param mixture_distribution The choice of mixing distribution. Takes values of "Lognormal" or "Bernoulli"
#' @param dispersion The total amount of dispersion for all events, defined as the variance/mean -1
#' @param dynamic_rates If TRUE then the rates include frequency-intensity dependence. 
#' @keywords Dynamic Rate
#' @export
#' @examples
#' DynamicRate()
#Function to fit beta parameters for lognormal/bernoulli clustering models
DynamicRate <- function(elt = NULL, mixture_distribution = NULL, 
                dispersion = 0, extreme_clusters = TRUE) {
  # Order the elt by expected loss - largest to smallest
  elt <- elt[order(elt[, 2], decreasing = T), ]  
  y <- elt[, 2]
  rate <- elt[, 1]  # rates
  n <- length(rate)
  total_rate <- sum(rate)
  # Scaling factor for overdispersion parameter for constant dispersion
  scale_factor <- 1 / total_rate 
  phi <- dispersion * scale_factor
  # Parameters for frequency-intensity independent 
  if (extreme_clusters == FALSE) {
    if (mixture_distribution == "Bernoulli") {
      theta <- 0.5
      beta_1 <- rate * sqrt(phi / (theta * (1 - theta)))
      beta_0 <- rate - (beta_1 * theta)
    } else if (mixture_distribution == "Lognormal") {
      beta_1 <- sqrt(log(1  + phi))
      beta_0 <- log(rate) - 0.5 * log(1 + phi) 
    }
  # parameters for frequency-intensity dependent  
  } else {
    dynamic_phi <- c()
    rate_cumulative <- cumsum(rate)
    
    #Dispersion parameters for the 'dynamic' Bernoulli clustering model
    if (mixture_distribution == "Bernoulli") {
      theta <- 0.5
      strong_events <- rate[rate_cumulative <= sum(rate) / 2]   # Partition the elt into strong and weak events 
      weak_events <- rate[rate_cumulative > sum(rate) / 2]
      phi_strong = rep(dispersion / sum(strong_events), length(strong_events)) 
      beta_1_strong <- strong_events * sqrt(phi_strong[1] / (theta * (1 - theta)))
      dynamic_phi[1:length(strong_events)] <- phi_strong
      sum_beta <- sqrt(dispersion * sum(rate) / (theta * (1 - theta))) - sum(beta_1_strong)
      beta_1_weak <- sum_beta * (weak_events / sum(weak_events))
      rate_dispersion <- sum(sqrt(phi_strong) * rate[1:length(strong_events)])
      dynamic_phi[(length(strong_events) + 1) : n] = ((sqrt(dispersion * sum(rate)) - 
                                                  (rate_dispersion)) / (sum(rate) - sum(strong_events))) ** 2
      beta_1 <- c(beta_1_strong, beta_1_weak)
      beta_0 <- rate - (beta_1 * theta)    
      
      #Dispersion parameters for the 'dynamic' Lognormal clustering model
    } else if (mixture_distribution == "Lognormal") {
      strong_events <- rate[rate_cumulative <= sum(rate)/100 ]
      phi_max <- rep(dispersion / sum(strong_events), length(strong_events))
      dynamic_phi[1:length(strong_events)] <- phi_max
      rate_dispersion <- sum(sqrt(phi_max) * rate[1 : length(strong_events)])
      for (i in (length(strong_events) + 1) : length(rate)) {
        dynamic_phi[i] = ((sqrt(dispersion * rate_cumulative[i]) - (rate_dispersion)) / rate[i]) ** 2
        rate_dispersion <- rate_dispersion + (rate[i] * sqrt(dynamic_phi[i]))
      }
      beta_1 = sqrt(log(1 + dynamic_phi))
      beta_0 = log(rate) - 0.5 * (beta_1 ** 2)
    }
  }
  invisible(list(beta_1 = beta_1, beta_0 = beta_0))
}