version 1.0-0

DESCRIPTION

@@ -1,19 +1,19 @@
 Package: actuar
 Type: Package
 Title: Actuarial functions
-Version: 0.9-7
-Date: 2008-03-25
-Author: Vincent Goulet, S�bastien Auclair, Christophe Dutang, Tommy Ouellet, Louis-Philippe Pouliot, Mathieu Pigeon
+Version: 1.0-0
+Date: 2008-09-15
+Author: Vincent Goulet, S�bastien Auclair, Christophe Dutang, Xavier Milhaud, Tommy Ouellet, Louis-Philippe Pouliot, Mathieu Pigeon
 Maintainer: Vincent Goulet <[email protected]>
 URL: http:https://www.actuar-project.org
-Description: Collection of functions and data sets related to
- actuarial science applications, mostly loss distributions, risk
- theory (including ruin theory), simulation of compound hierarchical
- models and credibility theory, for the moment.
+Description: Additional actuarial science functionality, mostly in the
+ fields of loss distributions, risk theory (including ruin theory),
+ simulation of compound hierarchical models and credibility theory,
+ for the moment.
 Depends: R (>= 2.6.0)
 License: GPL (>= 2)
 Encoding: latin1
 LazyLoad: yes
 LazyData: yes
 ZipData: yes
-Packaged: Tue Mar 25 15:43:40 2008; vincent
+Packaged: Mon Sep 15 13:52:31 2008; vincent

NAMESPACE

@@ -3,7 +3,7 @@ useDynLib(actuar, .registration = TRUE)
 
 ### Exports
 export(## Credibility theory
- bstraub, cm,
+ cm,
  ## Simulation of compound models
  simul, simpf, severity, unroll,
  ## Risk theory
@@ -82,7 +82,6 @@ S3method(plot, ogive)
 S3method(plot, ruin)
 
 S3method(predict, bstraub)
-S3method(predict, bstraub.old)
 S3method(predict, cm)
 S3method(predict, hache)
 S3method(predict, hierarc)
@@ -97,6 +96,7 @@ S3method(print, summary.aggregateDist)
 S3method(print, summary.cm)
 
 S3method(quantile, aggregateDist)
+S3method(quantile, grouped.data)
 
 S3method(severity, default)
 S3method(severity, portfolio)

R/aggregateDist.R

@@ -38,7 +38,6 @@ aggregateDist <-
  FUN <- npower(moments[1], moments[2], moments[3])
  comment(FUN) <- "Normal Power approximation"
  }
-
  else if (method == "simulation")
  {
  if (missing(nb.simul))
@@ -48,7 +47,6 @@ aggregateDist <-
  FUN <- simS(nb.simul, model.freq = model.freq, model.sev = model.sev)
  comment(FUN) <- "Approximation by simulation"
  }
-
  else
  {
  ## "recursive" and "convolution" cases. Both require a

R/bstraub.R

@@ -1,13 +1,16 @@
 ### ===== actuar: an R package for Actuarial Science =====
 ###
-### Bühlmann-Straub credibility model calculations
+### Bühlmann-Straub credibility model calculations.
+###
+### Computation of the between variance estimators has been moved to
+### external functions bvar.unbiased() and bvar.iterative() to share
+### with hache().
 ###
 ### AUTHORS: Vincent Goulet <[email protected]>,
-### Sébastien Auclair, and Louis-Philippe Pouliot
+### Sébastien Auclair, Louis-Philippe Pouliot
 
 bstraub <- function(ratios, weights, method = c("unbiased", "iterative"),
- tol = sqrt(.Machine$double.eps), maxit = 100,
- echo = FALSE, old.format = TRUE)
+ tol = sqrt(.Machine$double.eps), maxit = 100, echo = FALSE)
 {
  ## If weights are not specified, use equal weights as in
  ## Bühlmann's model.
@@ -24,7 +27,7 @@ bstraub <- function(ratios, weights, method = c("unbiased", "iterative"),
  if (nrow(ratios) < 2)
  stop("there must be more than one node")
  if (!identical(which(is.na(ratios)), which(is.na(weights))))
- stop("missing values are not in the same positions in weights and in ratios")
+ stop("missing values are not in the same positions in 'weights' and in 'ratios'")
  if (all(!weights, na.rm = TRUE))
  stop("no available data to fit model")
 
@@ -44,106 +47,97 @@ bstraub <- function(ratios, weights, method = c("unbiased", "iterative"),
 
  ## Collective weighted average.
  weights.ss <- sum(weights.s)
- ratios.ww <- sum(weights.s * ratios.w) / weights.ss
 
- ## Estimation of s^2.
+ ## Estimation of s^2
  s2 <- sum(weights * (ratios - ratios.w)^2, na.rm = TRUE) / (ntotal - ncontracts)
 
  ## First estimation of a. Always compute the unbiased estimator.
- ac <- weights.ss * (sum(weights.s * (ratios.w - ratios.ww)^2) - (ncontracts - 1) * s2) / (weights.ss^2 - sum(weights.s^2))
+ a <- bvar.unbiased(ratios.w, weights.s, s2, ncontracts)
 
  ## Iterative estimation of a. Compute only if
  ## 1. asked to in argument;
- ## 2. the unbiased estimator is > 0;
- ## 3. weights are not all equal (Bühlmann model).
+ ## 2. weights are not all equal (Bühlmann model).
+ ## 3. the unbiased estimator is > 0;
  method <- match.arg(method)
 
- if (method == "iterative")
+ if (method == "iterative" &&
+ diff(range(weights, na.rm = TRUE)) > .Machine$double.eps^0.5)
  {
- if (ac > 0)
- {
- if (diff(range(weights, na.rm = TRUE)) > .Machine$double.eps^0.5)
- {
- if (echo)
- {
- cat("Iteration\tBetween variance estimator\n")
- exp <- expression(cat(" ", count, "\t\t ", at1 <- at,
- fill = TRUE))
- }
- else
- exp <- expression(at1 <- at)
-
- at <- ac
- count <- 0
- repeat
- {
- eval(exp)
-
- if (maxit < (count <- count + 1))
- {
- warning("maximum number of iterations reached before obtaining convergence")
- break
- }
-
- cred <- 1 / (1 + s2/(weights.s * at))
- ratios.zw <- sum(cred * ratios.w) / sum(cred)
- at <- sum(cred * (ratios.w - ratios.zw)^2) / (ncontracts - 1)
-
- if (abs((at - at1)/at1) < tol)
- break
- }
- }
+ a <-
+ if (a > 0)
+ bvar.iterative(ratios.w, weights.s, s2, ncontracts, start = a,
+ tol = tol, maxit = maxit, echo = echo)
  else
- at <- ac
- }
- else
- at <- 0
- a <- at
- }
- else
- {
- a <- ac
- at <- NULL
+ 0
  }
 
  ## Final credibility factors and estimator of the collective mean.
  if (a > 0)
  {
  cred <- 1 / (1 + s2/(weights.s * a))
- ratios.zw <- sum(cred * ratios.w) / sum(cred)
+ ratios.zw <- drop(crossprod(cred, ratios.w)) / sum(cred)
  }
  else
  {
  cred <- numeric(length(weights.s))
- ratios.zw <- ratios.ww
+ ratios.zw <- drop(crossprod(weights.s, ratios.w)) / sum(weights.s)
  }
 
- if (old.format)
- {
- warning("this output format is deprecated")
- structure(list(individual = ratios.w,
- collective = ratios.zw,
- weights = weights.s,
- s2 = s2,
- unbiased = ac,
- iterative = at,
- cred = cred),
- class = "bstraub.old",
- model = "Buhlmann-Straub")
- }
- else
- structure(list(means = list(ratios.zw, ratios.w),
- weights = list(if (a > 0) sum(cred) else weights.ss, weights.s),
- unbiased = c(ac, s2),
- iterative = if (!is.null(at)) c(at, s2),
- cred = cred,
- nodes = list(nrow(weights))),
- class = "bstraub",
- model = "Buhlmann-Straub")
+ structure(list(means = list(ratios.zw, ratios.w),
+ weights = list(if (a > 0) sum(cred) else weights.ss, weights.s),
+ unbiased = if (method == "unbiased") c(a, s2),
+ iterative = if (method == "iterative") c(a, s2),
+ cred = cred,
+ nodes = list(nrow(weights))),
+ class = "bstraub",
+ model = "Buhlmann-Straub")
 }
 
 predict.bstraub.old <- function(object, ...)
  object$collective + object$cred * (object$individual - object$collective)
 
 predict.bstraub <- function(object, levels = NULL, newdata, ...)
  object$means[[1]] + object$cred * (object$means[[2]] - object$means[[1]])
+
+bvar.unbiased <- function(x, w, within, n)
+{
+ w.s <- sum(w)
+ x.w <- drop(crossprod(w, x)) / w.s
+
+ w.s * (drop(crossprod(w, (x - x.w)^2)) - (n - 1) * within) / (w.s^2 - sum(w^2))
+}
+
+bvar.iterative <- function(x, w, within, n, start,
+ tol = sqrt(.Machine$double.eps), maxit = 100,
+ echo = FALSE)
+{
+ if (echo)
+ {
+ cat("Iteration\tBetween variance estimator\n")
+ exp <- expression(cat(" ", count, "\t\t ", a1 <- a, fill = TRUE))
+ }
+ else
+ exp <- expression(a1 <- a)
+
+ a <- start
+ count <- 0
+
+ repeat
+ {
+ eval(exp)
+
+ if (maxit < (count <- count + 1))
+ {
+ warning("maximum number of iterations reached before obtaining convergence")
+ break
+ }
+
+ cred <- 1 / (1 + within/(w * a))
+ x.z <- drop(crossprod(cred, x)) / sum(cred)
+ a <- drop(crossprod(cred, (x - x.z)^2)) / (n - 1)
+
+ if (abs((a - a1)/a1) < tol)
+ break
+ }
+ a
+}

R/cm.R

@@ -5,7 +5,8 @@
 ### AUTHORS: Louis-Philippe Pouliot, Tommy Ouellet,
 ### Vincent Goulet <[email protected]>.
 
-cm <- function(formula, data, ratios, weights, subset, xreg = NULL,
+cm <- function(formula, data, ratios, weights, subset,
+ regformula = NULL, regdata, adj.intercept = FALSE,
  method = c("Buhlmann-Gisler", "Ohlsson", "iterative"),
  tol = sqrt(.Machine$double.eps), maxit = 100,
  echo = FALSE)
@@ -38,7 +39,7 @@ cm <- function(formula, data, ratios, weights, subset, xreg = NULL,
  ##
  if (any(duplicated(level.numbers)))
  stop("unsupported interactions in 'formula'")
- if (nlevels > 1 && !is.null(xreg))
+ if (nlevels > 1 && !is.null(regformula))
  stop("hierarchical regression models not supported")
  if (missing(ratios) & !missing(weights))
  stop("ratios have to be supplied if weights are")
@@ -107,30 +108,44 @@ cm <- function(formula, data, ratios, weights, subset, xreg = NULL,
  as.matrix(eval(cl, parent.frame())) # weights as matrix
  }
 
- ## Dispatch to appropriate calculation function
+ ## == DISPATCH TO APPROPRIATE CALCULATION FUNCTION ==
+ ##
+ ## Buhlmann-Straub models are handled by bstraub(), regression
+ ## models by hache() and hierarcahical models by hierarc().
  if (nlevels < 2) # one-dimensional model
  {
- if (is.null(xreg)) # Buhlmann-Straub
- {
- ## bstraub() accepts only "unbiased" and "iterative" for
- ## argument 'method'.
- method <- match.arg(method)
- if (method == "Buhlmann-Gisler" || method == "Ohlsson")
- method <- "unbiased"
+ ## One-dimensional models accept only "unbiased" and
+ ## "iterative" for argument 'method'.
+ method <- match.arg(method)
+ if (method == "Buhlmann-Gisler" || method == "Ohlsson")
+ method <- "unbiased"
 
- ## *** The 'old.format = FALSE' argument is necessary in
- ## *** the deprecation phase of the output format of ***
- ## bstraub(). To delete later.
+ if (is.null(regformula)) # Buhlmann-Straub
+ {
  res <- bstraub(ratios, weights, method = method,
- tol = tol, maxit = maxit, echo = echo,
- old.format = FALSE)
+ tol = tol, maxit = maxit, echo = echo)
  }
  else # Hachemeister
- res <- hache(ratios, weights, xreg, tol = tol,
- maxit = maxit, echo = echo)
+ {
+ ## If regression model is actually empty or has only an
+ ## intercept, call bstraub().
+ trf <- terms(regformula)
+ res <-
+ if (length(attr(trf, "factors")) == 0)
+ {
+ warning("empty regression model; fitting with Buhlmann-Straub's model")
+ bstraub(ratios, weights, method = method,
+ tol = tol, maxit = maxit, echo = echo,
+ old.format = FALSE)
+ }
+ else
+ hache(ratios, weights, regformula, regdata,
+ adj.intercept = adj.intercept,
+ method = method, tol = tol,
+ maxit = maxit, echo = echo)
+ }
 
- ## Add quantities not taken into account in calculation
- ## functions to results list.
+ ## Add missing quantities to results.
  res$classification <- levs
  res$ordering <- list(seq_along(levs))
  }