### Replication code for "Feller–Pareto and Related Distributions: Numerical
### Implementation and Actuarial Applications"

###################################################
### code chunk number 1: FP-numeric-comput.Rnw:54-70
###################################################
options(prompt = "R> ", continue = "+", width = 70, useFancyQuotes = FALSE)

## Running this paper requires packages 'rbenchmark', 'qgam'
## 'fitdistrplus', 'SMPracticals' and version 3.0-0 or later
## of actuar.
if (!all(c("actuar", "qgam", "rbenchmark", "fitdistrplus", "SMPracticals")
         %in% installed.packages()))
    stop("missing packages")
if (packageVersion("actuar") < "3.0.0")
    stop("version 3.0-0 or later of actuar needed")

## Load required packages
library("actuar")
library("qgam")
library("rbenchmark")
library("fitdistrplus")


###################################################
### code chunk number 2: FP-numeric-comput.Rnw:429-435
###################################################
rfpareto_gamma <- function(n, min, shape1, shape2, shape3, scale)
{
    x <- rgamma(n, shape3, 1)
    y <- rgamma(n, shape3, 1)
    min + scale * (x/y)^(1/shape2)
}


###################################################
### code chunk number 3: FP-numeric-comput.Rnw:454-459
###################################################
rfpareto_beta <- function(n, min, shape1, shape2, shape3, scale)
{
    x <- rbeta(n, shape1, shape3)
    min + scale * (1/x - 1)^(1/shape2)
}


###################################################
### code chunk number 4: FP-numeric-comput.Rnw:472-477
###################################################
rfpareto_actuar <- function(n, min, shape1, shape2, shape3, scale)
{
    x <- rbeta(n, shape3, shape1)
    min + scale * (1/x - 1)^(-1/shape2)
}


###################################################
### code chunk number 5: FP-numeric-comput.Rnw:838-858
###################################################
pfpareto  <- function(x, min, shape1, shape2, shape3, scale)
{
    y <- ((x - min)/scale)^shape2
    pbeta(y/(1 + y), shape3, shape1, lower = TRUE)
}

pfpareto2  <- function(x, min, shape1, shape2, shape3, scale)
{
    y <- ((x - min)/scale)^shape2
    pbeta(1/(1 + y), shape1, shape3, lower = FALSE)
}

shape1 <- 0.100000
shape3 <- 800
par(mar = c(4, 4, 1, 1))
curve(pfpareto(x, min = 0, shape1, shape2 = 1, shape3, scale = 1),
      from = 1e12, to = 1.01e16, ylab = expression(F(x)),
      lty = 2, col = 2, lwd = 2)
curve(pfpareto2(x, min = 0, shape1, shape2 = 1, shape3, scale = 1),
      add = TRUE, lwd = 2)


###################################################
### code chunk number 6: FP-numeric-comput.Rnw:893-940
###################################################
if (file.exists("cache/benchmark.rgamma.rds")) {
  benchmark.rgamma <- readRDS("cache/benchmark.rgamma.rds")
 } else {
  n <- 1e6
  shape <- runif(2 * n, 0, 100)
  shape1 <- tail(shape, n)
  shape2 <- head(shape, n)
  benchmark.rgamma <- benchmark(rgamma(2 * n, shape),
                                rbeta(n, shape1, shape2))[, 1:4]
  if(!dir.exists("cache/"))
    dir.create("cache/")
  saveRDS(benchmark.rgamma, "cache/benchmark.rgamma.rds")
}
## convert relative times to character mode with equal number of
## decimals for display in table
benchmark.rgamma[, 4] <- format(benchmark.rgamma[, 4], digits = 4)

if (file.exists("cache/benchmark.rfpareto.rds")) {
  benchmark.rfpareto <- readRDS("cache/benchmark.rfpareto.rds")
  } else {
  n <- 1e6
  shape <- runif(2 * n, 0, 100)
  shape1 <- head(shape, n)
  shape2 <- tail(shape, n)
  shape3 <- runif(n, 100, 1000)
  benchmark.rfpareto <-
    benchmark(rfpareto_gamma(n, min = 0, shape1, shape2, shape3, 1),
              rfpareto_beta(n, min = 0, shape1, shape2, shape3, 1),
              rfpareto_actuar(n, min = 0, shape1, shape2, shape3, 1))[, 1:4]
  
  saveRDS(benchmark.rfpareto, "cache/benchmark.rfpareto.rds")
}
## convert relative times to character mode with equal number of
## decimals for display in table
benchmark.rfpareto[, 4] <- format(benchmark.rfpareto[, 4], digits = 4)

## change expressions in first column to references to algorithms
pos <- sapply(c("gamma", "beta", "actuar"), grep, benchmark.rfpareto[, 1])
benchmark.rfpareto[pos, 1] <- paste0("Algorithm \\\\ref{alg:", names(pos), "}")

benchmark.rfpareto


###################################################
### code chunk number 7: FP-numeric-comput.Rnw:1027-1054
###################################################
pfpareto <- function(x, min, shape1, shape2, shape3, scale)
    pbeta(1/(1 + (x/scale)^shape2), shape1, shape3, lower.tail = FALSE);

n <- 1e6
shape1 <- 0.1
shape2 <- 10
shape3 <- 1e6
scale  <- 1

if (!exists("ecdf.rfpareto1"))
{
    if (file.exists("cache/ecdf.rfpareto1.rds"))
        ecdf.rfpareto1 <- readRDS("cache/ecdf.rfpareto1.rds")
    else
    {
        ecdf.rfpareto1 <-
            list(rfpareto_beta = ecdf(rfpareto_beta(n, min = 0, shape1, shape2, shape3, scale)),
                 rfpareto_actuar = ecdf(rfpareto_actuar(n, min = 0, shape1, shape2, shape3, scale)))

        saveRDS(ecdf.rfpareto1, "cache/ecdf.rfpareto1.rds")
    }
}
par(mar = c(4,4,1,1))
curve(pfpareto(x, min = 0, shape1, shape2, shape3, scale), from = 0, to = 1000,
      ylim = c(0, 1), ylab = expression(F[n](x)), lty = 1, lwd = 2)
curve(ecdf.rfpareto1$rfpareto_beta(x), add = TRUE, lty = 1, lwd = 2)
curve(ecdf.rfpareto1$rfpareto_actuar(x), add = TRUE, lty = 2, col = 2, lwd = 2)


###################################################
### code chunk number 8: FP-numeric-comput.Rnw:1059-1085
###################################################
set.seed(564546516)
n <- 1e6
shape1 <- 0.0001
shape2 <- 10
shape3 <- 0.0001
scale <- 1

if (!exists("ecdf.rfpareto2"))
{
    if (file.exists("cache/ecdf.rfpareto2.rds"))
        ecdf.rfpareto2 <- readRDS("cache/ecdf.rfpareto2.rds")
    else
    {
        ecdf.rfpareto2 <-
            list(rfpareto_beta = ecdf(rfpareto_beta(n, min = 0, shape1, shape2, shape3, scale)),
                 rfpareto_gamma = ecdf(rfpareto_gamma(n, min = 0, shape1, shape2, shape3, scale)))

        saveRDS(ecdf.rfpareto2, "cache/ecdf.rfpareto2.rds")
    }
}
par(mar = c(4,4,1,1))
curve(pfpareto(x, min = 0, shape1, shape2, shape3, scale), from = 0, to = 20,
      ylab = expression(F[n](x)), lty = 1, lwd = 2,
      ylim = c(0.498, 0.502))
curve(ecdf.rfpareto2$rfpareto_beta(x), add = TRUE, lty = 2, col = 2, lwd = 2)
curve(ecdf.rfpareto2$rfpareto_gamma(x), add = TRUE, lty = 3, col = 3, lwd = 2)


###################################################
### code chunk number 9: FP-numeric-comput.Rnw:1273-1336
###################################################
levtrbeta.old <- function(limit, shape1, shape2, shape3, scale, order = 1)
{
    u1m <- exp(-log1pexp(shape2 * (log(limit) - log(scale))))
    u <- exp(-log1pexp(-shape2 * (log(limit) - log(scale))))
    tmp <- order/shape2
    a <- shape3 + tmp
    b <- shape1 - tmp

    scale^order * betaint.old(u, a, b, u1m)/(gamma(shape1) * gamma(shape3)) + 
        limit^order * pbeta(u, shape3, shape1, lower.tail = FALSE)
}

betaint.old <- function(x, a, b, x1m)
{
    if (b > 0)
    {
        return(gamma(a) * gamma(b) * pbeta(x, a, b, lower.tail = TRUE))
    }

    r <- floor(-b)

    if (! (a - r - 1 > 0))
        return(NA);

    ap <- a
    bp <- b
    lx <- log(x)
    lx1m  <- log1p(-x)
    x1 <- exp(lx1m - lx)

    ap <- ap - 1
    ## c <- start   #(x^(a - 1) (1 - x)^b) / b
    c <- exp(ap * lx + bp * lx1m)/bp
    sum <- c
    ratio <- 1/bp
    bp <- bp + 1

    i <- 0
    while(i < r)
    {
        tmp <- ap/bp
        c <- c * tmp * x1
        sum <- sum + c;
        ratio <- ratio * tmp
        ap <- ap - 1
        bp <- bp + 1
        i <- i + 1
    }

    -gamma(a + b) * sum + (ratio * ap) *
        exp(lgamma(ap) + lgamma(bp) + pbeta(x, ap, bp, lower.tail = TRUE, log = TRUE))
}

## Generalization (to include 'order') of the function in actuar
elev <- function(x, order = 1)
{
    FUN <- function(limit)
        sapply(limit, function(x, y) mean(pmin(x, y)^order), x = x)
    environment(FUN) <- new.env()
    assign("x", sort(x), envir = environment(FUN))
    assign("order", order, envir = environment(FUN))
    FUN
}


###################################################
### code chunk number 10: FP-numeric-comput.Rnw:1342-1358
###################################################
n <- 1e5
shape1 <- 0.1
shape2 <- 14.5
shape3 <- 1
scale <- 1
order <- 1

E1 <- elev(rfpareto_beta(n, min = 0, shape1 = shape1, shape3 = shape3, shape2 = shape2, scale = scale),
           order = order)
par(mar = c(4,4,1,1))
curve(E1(x), from = 0, to = 25,
      ylim = c(0, levtrbeta.old(25, shape1, shape2, shape3, scale)),
      xlab = "Limit", ylab = "Limited Expected Value",
      lty = 2, col = 2, lwd = 2)
curve(levtrbeta.old(x, shape1, shape2, shape3, scale, order = order),
      add = TRUE, lty = 1, lwd = 2)


###################################################
### code chunk number 11: FP-numeric-comput.Rnw:1363-1373
###################################################
order <- 4

E4 <- elev(rfpareto_beta(n, min = 0, shape1 = shape1, shape3 = shape3, shape2 = shape2, scale = scale),
           order = order)
par(mar = c(4,4,1,1))
curve(E4(x), from = 0, to = 20,
      xlab = "Limit", ylab = "Fourth Limited Moment",
      lty = 2, col = 2, lwd = 2)
curve(levtrbeta.old(x, shape1, shape2, shape3, scale, order = order),
      add = TRUE, lty = 1, lwd = 2)


###################################################
### code chunk number 12: FP-numeric-comput.Rnw:1393-1401
###################################################
order <- 1
par(mar = c(4,4,1,1))
curve(E1(x), from = 0, to = 25,
      ylim = c(0, actuar::levtrbeta(25, shape1, shape2, shape3, scale)),
      xlab = "Limit", ylab = "Limited Expected Value",
      lty = 2, col = 2, lwd = 2)
curve(actuar::levtrbeta(x, shape1, shape2, shape3, scale, order = order),
      add = TRUE, lty = 1, lwd = 2)


###################################################
### code chunk number 13: FP-numeric-comput.Rnw:1406-1413
###################################################
order <- 4
par(mar = c(4,4,1,1))
curve(E4(x), from = 0, to = 20,
      xlab = "Limit", ylab = "Fourth Limited Moment",
      lty = 2, col = 2, lwd = 2)
curve(actuar::levtrbeta(x, shape1, shape2, shape3, scale, order = order),
      add = TRUE, lty = 1, lwd = 2)


###################################################
### code chunk number 14: FP-numeric-comput.Rnw:1426-1428
###################################################
## Cleanup local functions defined above that mask those from actuar
rm("elev", "pfpareto")


###################################################
### code chunk number 15: FP-numeric-comput.Rnw:1447-1449
###################################################
data("danish", package = "SMPracticals")
danish.loss <- as.numeric(danish)


###################################################
### code chunk number 16: FP-numeric-comput.Rnw:1455-1460
###################################################
par(mar = c(4, 4, 1, 1))
## plot(Loss ~ Date, data = danishuni, type = "h", main = "",
##      xlab = "Occurrence date", ylab = "Loss (Million DKK)")
Fn <- ecdf(danish.loss) # empirical CDF
plot(Fn, main = "")


###################################################
### code chunk number 17: FP-numeric-comput.Rnw:1465-1468
###################################################
x <- knots(Fn)             # data points
par(mar = c(4, 4, 1, 1))
plot(log(x), log(1 - Fn(x))) #, main = "Empirical Survival Function (log scale)")


###################################################
### code chunk number 18: FP-numeric-comput.Rnw:1485-1487
###################################################
summary(danish.loss)
quantile(danish.loss, probs = c(0.9, 0.95, 0.99, 1))


###################################################
### code chunk number 19: FP-numeric-comput.Rnw:1527-1532
###################################################
(fit.fp <- fitdist(danish.loss, "fpareto", optim.method = "L-BFGS-B",
                   fix.arg = list(min = 0),
                   start = list(shape1 = 1, shape2 = 1, shape3 = 1,
                                scale = 1),
                   lower = .Machine$double.eps))


###################################################
### code chunk number 20: FP-numeric-comput.Rnw:1543-1553
###################################################
theta1 <- c(16.094, 0.955)
theta2 <- c(1.555, 1.102)
theta <- 0.955 ; phi <- 9.854
pwiw <- function(x)
{
    p1 <- pweibull(x, theta1[1], theta1[2]) / (1 + phi) / 
      pweibull(theta, theta1[1], theta1[2])
    p2 <- (pinvweibull(x, theta2[1], theta2[2]) - pinvweibull(theta, theta2[1], theta2[2])) /
      pinvweibull(theta, theta2[1], theta2[2], lower.tail = FALSE)
    p2 <- p2 * phi/(1 + phi)
    p1 * (x <=  theta) + (1/(1 + phi) + p2) * (x > theta)
}


###################################################
### code chunk number 21: FP-numeric-comput.Rnw:1560-1566
###################################################
par(mar = c(4, 4, 1, 1))
cdfcomp(fit.fp, xlogscale = TRUE, main = "",
        do.points = FALSE, fitcol = 2, fitlty = 2, fitlwd = 2,
        datacol = "black", addlegend = FALSE,
        xlab = "Loss amount (log scale)")
curve(pwiw(x), add = TRUE, col = 3, lty = 3, lwd = 2)


###################################################
### code chunk number 22: FP-numeric-comput.Rnw:1571-1578
###################################################
par(mar = c(4, 4, 1, 1))
cdfcomp(fit.fp, xlogscale = TRUE, main = "",
        do.points = FALSE, fitcol = 2, fitlty = 2, fitlwd = 2,
        datacol = "black", addlegend = FALSE,
        ylim = c(.9, 1), xlim = c(5, 300),
        xlab = "Loss amount (log scale)")
curve(pwiw(x), add = TRUE, col = 3, lty = 3, lwd = 2)


###################################################
### code chunk number 23: FP-numeric-comput.Rnw:1592-1593
###################################################
gofstat(fit.fp)


###################################################
### code chunk number 24: FP-numeric-comput.Rnw:1626-1628
###################################################
par.fp <- fit.fp$estimate
ag <- format(par.fp["shape1"] * par.fp["shape2"], digits = 4)


###################################################
### code chunk number 25: FP-numeric-comput.Rnw:1636-1639
###################################################
par.fp <- fit.fp$estimate
mfpareto(1, min = 0, shape1 = par.fp["shape1"], shape2 = par.fp["shape2"],
         shape3 = par.fp["shape3"], scale = par.fp["scale"])


###################################################
### code chunk number 26: FP-numeric-comput.Rnw:1657-1664
###################################################
EXL.fpareto <- function(deduc, cover)
    levfpareto(deduc + cover, min = 0,
               shape1 = par.fp["shape1"], shape2 = par.fp["shape2"],
               shape3 = par.fp["shape3"], scale = par.fp["scale"]) -
        levfpareto(deduc, min = 0,
                   shape1 = par.fp["shape1"], shape2 = par.fp["shape2"],
                   shape3 = par.fp["shape3"], scale = par.fp["scale"])


###################################################
### code chunk number 27: FP-numeric-comput.Rnw:1669-1677
###################################################
par(mar = c(4, 4, 1, 1))
curve(EXL.fpareto(x, cover = 4), 1, 10,
      ylim = c(0, 1.3), lwd = 2,
      xlab = "Deductible d (million DKK)",
      ylab = "Expected Cost")
curve(EXL.fpareto(x, cover = 3), add = TRUE, col = 2, lty = 2, lwd = 2)
curve(EXL.fpareto(x, cover = 2), add = TRUE, col = 3, lty = 3, lwd = 2)
curve(EXL.fpareto(x, cover = 1), add = TRUE, col = 4, lty = 4, lwd = 2)


###################################################
### code chunk number 28: FP-numeric-comput.Rnw:1716-1720
###################################################
x <- rcomppois(1e6, lambda = 10,
               rfpareto(min = 0, shape1 = par.fp["shape1"],
                        shape2 = par.fp["shape2"], shape3 = par.fp["shape3"],
                        scale = par.fp["scale"]))


###################################################
### code chunk number 29: FP-numeric-comput.Rnw:1731-1732
###################################################
quantile(x, c(0.75, 0.9, 0.95, 0.99))