### ### TITLE: actuar: An R Package for Actuarial Science ### AUTHORS: Christophe Dutang, Vincent Goulet*, Mathieu Pigeon ### * Corresponding author ### AFFILIATION: École d'actuariat, Université Laval, Québec, Canada ### ADDRESS: Université Laval, École d'actuariat, Pavillon ### Alexandre-Vachon, 1045, avenue de la Médecine, local ### 1620, Québec (QC) G1V~0A6, Canada ### PHONE: +418 656 5736 ### FAX: +418 656 7790 ### EMAIL: vincent.goulet@act.ulaval.ca ### library("actuar") ### ### LOSS DISTRIBUTIONS MODELING ### ## Grouped data x <- grouped.data(Group = c(0, 25, 50, 100, 150, 250, 500), Line.1 = c(30, 31, 57, 42, 65, 84), Line.2 = c(26, 33, 31, 19, 16, 11)) class(x) x x[, 1] # group boundaries x[, -1] # group frequencies x[1:3,] # first 3 groups x[1, 2] <- 22; x # frequency replacement x[1, c(2, 3)] <- c(22, 19); x # frequency replacement x[1, 1] <- c(0, 20); x # boundary replacement x[c(3, 4), 1] <- c(55, 110, 160); x # boudaries replacement mean(x) # mean of grouped data hist(x[, -3]) # histogram Fnt <- ogive(x) # ogive knots(Fnt) # group boundaries Fnt(knots(Fnt)) # ogive at group boundaries plot(Fnt) # plot of the ogive ## Data sets data("dental"); dental data("gdental"); gdental ## Calculation of empirical moments emm(dental, order = 1:3) # first three moments emm(gdental, order = 1:3) # idem for grouped data lev <- elev(dental) # empirical LEV for ungrouped data lev(knots(lev)) # ELEV at data points plot(lev, type = "o", pch = 19) # plot of the ELEV function lev <- elev(gdental) # empirical LEV for grouped data lev(knots(lev)) # ELEV at data points plot(lev, type = "o", pch = 19) # plot of the ELEV function ## Minimum distance estimation mde(gdental, pexp, start = list(rate = 1/200), measure = "CvM") mde(gdental, pexp, start = list(rate = 1/200), measure = "chi-square") mde(gdental, levexp, start = list(rate = 1/200), measure = "LAS") ## Coverage modifications f <- coverage(pdf = dgamma, cdf = pgamma, deductible = 1, limit = 10) f(0, shape = 5, rate = 1) f(5, shape = 5, rate = 1) f(9, shape = 5, rate = 1) f(12, shape = 5, rate = 1) g <- coverage(dgamma, pgamma) # simpler case g ### ### RISK THEORY ### ## Discretization of claim amount distributions fu <- discretize(plnorm(x), method = "upper", from = 0, to = 5) fl <- discretize(plnorm(x), method = "lower", from = 0, to = 5) fr <- discretize(plnorm(x), method = "rounding", from = 0, to = 5) fb <- discretize(plnorm(x), method = "unbiased", from = 0, to = 5, lev = levlnorm(x)) op0 <- par(mfrow = c(2, 2), mar = c(5, 4, 4, 2)) curve(plnorm(x), from = 0, to = 5, lwd = 2, main = "Upper", ylab = "F(x)") op <- par(col = "blue") plot(stepfun(0:4, diffinv(fu)), pch = 20, add = TRUE) par(op) curve(plnorm(x), from = 0, to = 5, lwd = 2, main = "Lower", ylab = "F(x)") op <- par(col = "blue") plot(stepfun(0:5, diffinv(fl)), pch = 20, add = TRUE) par(op) curve(plnorm(x), from = 0, to = 5, lwd = 2, main = "Rounding", ylab = "F(x)") op <- par(col = "blue") plot(stepfun(0:4, diffinv(fr)), pch = 20, add = TRUE) par(op) curve(plnorm(x), from = 0, to = 5, lwd = 2, main = "Unbiased", ylab = "F(x)") op <- par(col = "blue") plot(stepfun(0:5, diffinv(fb)), pch = 20, add = TRUE) par(c(op0, op)) ## Calculation of the aggregate claim amount distribution fx <- discretize(pgamma(x, 2, 1), from = 0, to = 22, step = 0.5, method = "unbiased", lev = levgamma(x, 2, 1)) Fs <- aggregateDist("recursive", model.freq = "poisson", model.sev = fx, lambda = 10, x.scale = 0.5) summary(Fs) # summary method knots(Fs) # support of Fs.b (knots) plot(Fs, do.points = FALSE, verticals = TRUE, xlim = c(0, 60)) # plot mean(Fs) # empirical mean quantile(Fs) # quantiles quantile(Fs, 0.999) # quantile VaR(Fs) # value-at-risk CTE(Fs) # conditional tail expectation fx.u <- discretize(pgamma(x, 2, 1), from = 0, to = 22, step = 0.5, method = "upper") Fs.u <- aggregateDist("recursive", model.freq = "poisson", model.sev = fx.u, lambda = 10, x.scale = 0.5) fx.l <- discretize(pgamma(x, 2, 1), from = 0, to = 22, step = 0.5, method = "lower") Fs.l <- aggregateDist("recursive", model.freq = "poisson", model.sev = fx.l, lambda = 10, x.scale = 0.5) Fs.n <- aggregateDist("normal", moments = c(20, 60)) Fs.s <- aggregateDist("simulation", model.freq = expression(y = rpois(10)), model.sev = expression(y = rgamma(2, 1)), nb.simul = 10000) op0 <- par(col = "black") plot(Fs, do.points = FALSE, verticals = TRUE, xlim = c(0, 60), sub = "") par(col = "blue") plot(Fs.u, do.points = FALSE, verticals = TRUE, add = TRUE, sub = "") par(col = "red") plot(Fs.l, do.points = FALSE, verticals = TRUE, add = TRUE, sub = "") par(col = "green") plot(Fs.s, do.points = FALSE, verticals = TRUE, add = TRUE, sub = "") par(col = "magenta") plot(Fs.n, add = TRUE, sub = "") legend(30, 0.4, c("recursive + unbiased", "recursive + upper", "recursive + lower", "simulation", "normal approximation"), col = c("black", "blue", "red", "green", "magenta"), lty = 1, text.col = "black") par(op0) ## Adjustment coefficient adjCoef(mgf.claim = mgfexp(x), mgf.wait = mgfexp(x, 2), premium.rate = 2.4, upper = 1) # no reinsurance mgfx <- function(x, y) mgfexp(x * y) p <- function(x) 2.6 * x - 0.2 rho <- adjCoef(mgfx, mgfexp(x, 2), premium = p, upper = 1, reins = "prop", from = 0, to = 1) # with reinsurance rho(c(0.75, 0.8, 0.9, 1)) # evaluation for some limits plot(rho) # plot ## Probability of ruin psi <- ruin(claims = "e", par.claims = list(rate = 5), wait = "e", par.wait = list(rate = 3)) psi psi(0:10) ruin(claims = "e", par.claims = list(rate = c(3, 7), weights = 0.5), wait = "e", par.wait = list(rate = 3)) prob <- c(0.5614, 0.4386) rates <- matrix(c(-8.64, 0.101, 1.997, -1.095), 2, 2) ruin(claims = "p", par.claims = list(prob = prob, rates = rates), wait = "e", par.wait = list(rate = c(5, 1), weights = c(0.4, 0.6))) psi <- ruin(claims = "p", par.claims = list(prob = prob, rates = rates), wait = "e", par.wait = list(rate = c(5, 1), weights = c(0.4, 0.6))) plot(psi, from = 0, to = 50) ### ### SIMULATION OF COMPOUND HIERARCHICAL MODELS ### nodes <- list(cohort = 2, entity = c(4, 3), year = c(4, 4, 4, 4, 5, 5, 5)) mf <- expression(cohort = rexp(2), entity = rgamma(cohort, 1), year = rpois(weights * entity)) ms <- expression(cohort = rnorm(2, sqrt(0.1)), entity = rnorm(cohort, 1), year = rlnorm(entity, 1)) wijt <- runif(31, 0.5, 2.5) pf <- simul(nodes = nodes, model.freq = mf, model.sev = ms, weights = wijt) pf # default print method aggregate(pf) # aggregate claim amounts aggregate(pf, by = c("cohort", "year"), FUN = mean) frequency(pf) # number of claims frequency(pf, by = "cohort") severity(pf) # individual claim amounts severity(pf, splitcol = 1) weights(pf) # weights ### ### CREDIBILITY THEORY ### ## Hierarchical credibility model X <- cbind(cohort = c(1, 2, 1, 2, 2), hachemeister) fit <- cm(~cohort + cohort:state, data = X, ratios = ratio.1:ratio.12, weights = weight.1:weight.12, method = "iterative") fit predict(fit) # credibility premiums summary(fit) # summary of results summary(fit, levels = "cohort") # summary for cohorts only predict(fit, levels = "cohort") # cohorts cred. premiums only ## Regression model of Hachemeister fit <- cm(~state, hachemeister, xreg = 12:1, ratios = ratio.1:ratio.12, weights = weight.1:weight.12) fit predict(fit, newdata = 0) # future value of regressor