###############################################
#        CASUALTY ACTUARIAL SOCIETY           #
#      OPEN-SOURCE SOFTWARE COMMITTEE         #
###############################################
#       actuar I Companion Script             #
# http://opensourcesoftware.casact.org/actuar #
###############################################

#Use this script as you read along the article
#to quickly execute the examples!


#INTRODUCTION

library(actuar)


#DATA

data(dental)
data(gdental)

dental

gdental


#PROBABILITY LAWS//FUNCTIONS//PROBABILITY DENSITY FUNCTION

dpareto(1.0, shape = 10, scale = 5)


#PROBABILITY LAWS//FUNCTIONS//CUMULATIVE DISTRIBUTION FUNCTION

ppareto(1.0, shape = 10, scale = 5)


#PROBABILITY LAWS//FUNCTIONS//QUANTILE FUNCTION

qpareto(.838, shape = 10, scale = 5)


#PROBABILITY LAWS//FUNCTIONS//RANDOM DATA

r <- rpareto(1000, shape = 10, scale = 5)


#PROBABILITY LAWS//FUNCTIONS//RAW MOMENTS

m1r <- mpareto(1, shape = 10, scale = 5)
m2r <- mpareto(2, shape = 10, scale = 5)
variancer <- m2r - m1r^2
stdevr <- variancer^.5

m1r
stdevr
variancer


#PROBABILITY LAWS//FUNCTIONS//LIMITED MOMENTS

m1l <- levpareto(1.0, shape = 10, scale = 5, order = 1)
m2l <- levpareto(1.0, shape = 10, scale = 5, order = 2)
variancel <- m2l - m1l^2
stdevl <- variancel^.5

m1l
stdevl
variancel

rvl <- matrix(c(m1r,m1l,stdevr,stdevl,variancer,variancel),ncol=2,byrow=TRUE)
colnames(rvl) <- c("Raw","Limited")
rownames(rvl) <- c("Mean","Std Dev","Variance")
rvl <- as.table(rvl)

rvl


#PROBABILITY LAWS//DISTIBUTIONS

showgraphs <- function(fun, par, what = c("d", "p", "m", "lev"), xlim)
{
  dist <- switch(fun,
                 trbeta = "TRANSFORMED BETA DISTRIBUTION",
                 genpareto = "GENERALIZED PARETO DISTRIBUTION",
                 burr = "BURR DISTRIBUTION",
                 invburr = "INVERSE BURR DISTRIBUTION",
                 pareto = "PARETO DISTRIBUTION",
                 invpareto = "INVERSE PARETO DISTRIBUTION",
                 llogis = "LOGLOGISTIC DISTRIBUTION",
                 paralogis = "PARALOGISTIC DISTRIBUTION",
                 invparalogis = "INVERSE PARALOGISTIC DISTRIBUTION",
                 trgamma = "TRANSFORMED GAMMA DISTRIBUTION",
                 invtrgamma = "INVERSE TRANSFORMED GAMMA DISTRIBUTION",
                 invgamma = "INVERSE GAMMA DISTRIBUTION",
                 weibull = "WEIBULL DISTRIBUTION",
                 invweibull = "INVERSE WEIBULL DISTRIBUTION",
                 invexp = "INVERSE EXPONENTIAL DISTRIBUTION",
                 pareto1 = "SINGLE PARAMETER PARETO DISTRIBUTION",
                 lgamma = "LOGGAMMA DISTRIBUTION",
                 genbeta = "GENERALIZED BETA DISTRIBUTION",
                 phtype = "PHASE-TYPE DISTRIBUTION",
                 gamma = "GAMMA DISTRIBUTION",
                 exp = "EXPONENTIAL DISTRIBUTION",
                 chisq = "CHI-SQUARE DISTRIBUTION",
                 lnorm = "LOGNORMAL DISTRIBUTION",
                 invgauss = "INVERSE GAUSSIAN DISTRIBUTION",
                 norm = "NORMAL DISTRIBUTION",
                 beta = "BETA DISTRIBUTION",
                 unif = "UNIFORM DISTRIBUTION")
  
  if (missing(xlim))
  {
    qf <- match.fun(paste("q", fun, sep = ""))
    formals(qf)[names(par)] <- par
    xlim <- c(0, qf(0.999))
  }
  k <- seq.int(4)
  limit <- seq(0, xlim[2], len = 10)
  mfrow = c(ceiling(length(what) / 2), 2)
  op <- par(mfrow = mfrow, oma = c(0, 0, 2, 0))
  
  for (t in what)
  {
    f   <- match.fun(paste(t, fun, sep = ""))
    formals(f)[names(par)] <- par
    
    main <- switch(t,
                   "d" = "Probability Density Function",
                   "p" = "Cumulative Distribution Function",
                   "m" = "Raw Moments",
                   "lev" = "Limited Expected Value Function",
                   "mgf" = "Moment Generating Function")
    
    if (t == "m")
      plot(k, f(k), type = "l", col = 4, lwd = 2, main = main)
    else if (t == "lev")
      plot(limit, f(limit), type = "l", col = 4, lwd = 2, main = main)
    else if (t == "mgf")
      curve(f(x), xlim = c(0, 2), col = 4, lwd = 2, main = main)
    else
      curve(f(x), xlim = xlim, col = 4, lwd = 2, main = main)
    title(main = dist, outer = TRUE)
  }
  par(op)
}


showgraphs("invgamma", list(shape = 6, scale = 10))


#PROBABILITY LAWS//DISTIBUTIONS//TRANSFORMED BETA FAMILY

showgraphs("trbeta", list(shape1 = 3, shape2 = 4, shape3 = 5, scale = 10))

showgraphs("genpareto", list(shape1 = 10, shape2 = 4, scale = 10))

showgraphs("burr", list(shape1 = 3, shape2 = 4, scale = 10))

showgraphs("invburr", list(shape1 = 3, shape2 = 6, scale = 10))

showgraphs("pareto", list(shape = 10, scale = 10))

showgraphs("invpareto", list(shape = 4, scale = 1), what = c("d", "p"))

showgraphs("llogis", list(shape = 6, scale = 10))

showgraphs("paralogis", list(shape = 3, scale = 10))

showgraphs("invparalogis", list(shape = 6, scale = 10))


#PROBABILITY LAWS//DISTRIBUTIONS//TRANSFORMED GAMMA FAMILY

showgraphs("trgamma", list(shape1 = 3, shape2 = 1, scale = 10))

showgraphs("invtrgamma", list(shape1 = 3, shape2 = 2, scale = 10))

showgraphs("invgamma", list(shape = 6, scale = 10))

showgraphs("weibull", list(shape = 1.5, scale = 10))

showgraphs("invweibull", list(shape = 6, scale = 10))

showgraphs("invexp", list(rate = 1), what = c("d", "p"))


#PROBABILITY LAWS//DISTRIBUTIONS//OTHER

showgraphs("pareto1", list(shape = 5, min = 10), xlim = c(0, 50))

showgraphs("lgamma", list(shapelog = 2, ratelog = 5))

showgraphs("genbeta", list(shape1 = 1, shape2 = 2, shape3 = 3, scale = 2))

showgraphs("phtype", list(prob = c(0.5614, 0.4386), rates = matrix(c(-8.64, 0.101, 1.997, -1.095), 2, 2)), what = c("d", "p", "m", "mgf"), xlim = c(0.001, 5))


#PROBABILITY LAWS//DISTRIBUTIONS//INCLUDED IN R

showgraphs("gamma", list(shape = 3, rate = 5), what = c("m", "lev", "mgf"))

showgraphs("chisq", list(df = 3), what = c("m", "lev", "mgf"))

showgraphs("exp", list(rate = 5), what = c("m", "lev", "mgf"))

showgraphs("lnorm", list(meanlog = 1, sdlog = 1), what = c("m", "lev"))

showgraphs("invgauss", list(nu = 1, lambda = 10), what = c("m", "lev", "mgf"), xlim = c(0, 10))

showgraphs("norm", list(mean = 0, sd = 1), what = c("m", "mgf"))

showgraphs("beta", list(shape1 = 1, shape2 = 2), what = c("m", "lev"))

showgraphs("unif", list(min = 0, max = 1), what = c("m", "lev", "mgf"))


#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))

x


#GROUPED DATA//EXTRACTION AND REPLACEMENT

x[, 1]

x[, -1]

x[1:3, ]


x[1, 2] <- 22

x


x[1, c(2, 3)] <- c(22, 19)

x


x[1, 1] <- c(0, 20)

x


x[c(3, 4), 1] <- c(55, 110, 160)

x


#GROUPED DATA//GROUP MEAN AND HISTOGRAMS

mean(x)

hist(x[, -3])


#GROUPED DATA//OGIVE

CDF <- ogive(x)

CDF


knots(CDF)

plot(CDF)

CDF(250)

CDF(knots(CDF))


#EMPIRICAL MOMENTS//KTH EMPIRICAL MOMENT

emm(dental, order = 1:3)

emm(gdental, order = 1:3)


#EMPIRICAL MOMENTS//EMPIRICAL LIMITED EXPECTED VALUE

lev <- elev(dental)

lev


glev <- elev(gdental)

glev


lev(knots(lev))

glev(knots(glev))

plot(lev, type="o", pch = 19)

plot(glev, type="o", pch = 19)


#OPTIMIZATION//MINIMUM DISTANCE ESTIMATION

mde(gdental, pexp, start = list(rate = 1/200), measure = "CvM")

mde(gdental, levexp, start = list(rate = 1/200), measure = "LAS")


#OPTIMIZATION//MINIMUM DISTANCE ESTIMATION//ISSUES

mde(gdental, ppareto, start = list(shape = 3, scale = 600), measure = "CvM")

pparetolog <- function(x, logshape, logscale) ppareto(x,exp(logshape),exp(logscale))


pl <- mde(gdental, pparetolog, start = list(logshape = log(3),logscale = log(600)), measure = "CvM")

pl


exp(pl$estimate)


#COVERAGE MODIFICATIONS

mcovpdf <- coverage(pdf = dgamma, cdf = pgamma, deductible = 1, limit = 10)

mcovcdf <- coverage(cdf = pgamma, deductible = 1, limit = 10)

mcovpdf(5, shape = 5, scale = 3)

mcovcdf(5, shape = 5, scale = 3)


#SOURCES:

#C. Dutang, V. Goulet and M. Pigeon (2008). 
#actuar: An R Package for Actuarial Science. 
#Journal of Statistical Software, vol. 25, no. 7, 1-37. 
#URL http://www.jstatsoft.org/v25/i07

#R Development Core Team (2012). 
#R: A language and environment for statistical computing. 
#R Foundation for Statistical Computing, Vienna, Austria. 
#ISBN 3-900051-07-0, URL http://www.R-project.org/.