Package 'actuar'

Actuarial Functions and Heavy Tailed Distributions

Description

Functions and data sets for actuarial science: modeling of loss distributions; risk theory and ruin theory; simulation of compound models, discrete mixtures and compound hierarchical models; credibility theory. Support for many additional probability distributions to model insurance loss size and frequency: 23 continuous heavy tailed distributions; the Poisson-inverse Gaussian discrete distribution; zero-truncated and zero-modified extensions of the standard discrete distributions. Support for phase-type distributions commonly used to compute ruin probabilities. Main reference: <doi:10.18637/jss.v025.i07>. Implementation of the Feller-Pareto family of distributions: <doi:10.18637/jss.v103.i06>.

Details

actuar provides additional actuarial science functionality and support for heavy tailed distributions to the R statistical system.

The current feature set of the package can be split into five main categories.

1. Additional probability distributions: 23 continuous heavy tailed distributions from the Feller-Pareto and Transformed Gamma families, the loggamma, the Gumbel, the inverse Gaussian and the generalized beta; phase-type distributions; the Poisson-inverse Gaussian discrete distribution; zero-truncated and zero-modified extensions of the standard discrete distributions; computation of raw moments, limited moments and the moment generating function (when it exists) of continuous distributions. See the “distributions” package vignette for details.

2. Loss distributions modeling: extensive support of grouped data; functions to compute empirical raw and limited moments; support for minimum distance estimation using three different measures; treatment of coverage modifications (deductibles, limits, inflation, coinsurance). See the “modeling” and “coverage” package vignettes for details.

3. Risk and ruin theory: discretization of the claim amount distribution; calculation of the aggregate claim amount distribution; calculation of the adjustment coefficient; calculation of the probability of ruin, including using phase-type distributions. See the “risk” package vignette for details.

4. Simulation of discrete mixtures, compound models (including the compound Poisson), and compound hierarchical models. See the “simulation” package vignette for details.

5. Credibility theory: function cm fits hierarchical (including Bühlmann, Bühlmann-Straub), regression and linear Bayes credibility models. See the “credibility” package vignette for details.

Author(s)

Christophe Dutang, Vincent Goulet, Mathieu Pigeon and many other contributors; use packageDescription("actuar") for the complete list.

Maintainer: Vincent Goulet.

References

Dutang, C., Goulet, V. and Pigeon, M. (2008). actuar: An R Package for Actuarial Science. Journal of Statistical Software, 25(7), 1–37. doi:10.18637/jss.v025.i07.

Dutang, C., Goulet, V., Langevin, N. (2022). Feller-Pareto and Related Distributions: Numerical Implementation and Actuarial Applications. Journal of Statistical Software, 103(6), 1–22. doi:10.18637/jss.v103.i06.

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

See Also

For probability distributions support functions, use as starting points: FellerPareto, TransformedGamma, Loggamma, Gumbel, InverseGaussian, PhaseType, PoissonInverseGaussian and, e.g., ZeroTruncatedPoisson, ZeroModifiedPoisson.

For loss modeling support functions: grouped.data, ogive, emm, elev, mde, coverage.

For risk and ruin theory functions: discretize, aggregateDist, adjCoef, ruin.

For credibility theory functions and datasets: cm, hachemeister.

Examples

## The package comes with extensive demonstration scripts;
## use the following command to obtain the list.
## Not run: demo(package = "actuar")

---

Adjustment Coefficient

Description

Compute the adjustment coefficient in ruin theory, or return a function to compute the adjustment coefficient for various reinsurance retentions.

Usage

adjCoef(mgf.claim, mgf.wait = mgfexp, premium.rate, upper.bound,
        h, reinsurance = c("none", "proportional", "excess-of-loss"),
        from, to, n = 101)

## S3 method for class 'adjCoef'
plot(x, xlab = "x", ylab = "R(x)",
     main = "Adjustment Coefficient", sub = comment(x),
     type = "l", add = FALSE, ...)

Arguments

mgf.claim	an expression written as a function of x or of x and y, or alternatively the name of a function, giving the moment generating function (mgf) of the claim severity distribution.
mgf.wait	an expression written as a function of x, or alternatively the name of a function, giving the mgf of the claims interarrival time distribution. Defaults to an exponential distribution with mean 1.
premium.rate	if reinsurance = "none", a numeric value of the premium rate; otherwise, an expression written as a function of y, or alternatively the name of a function, giving the premium rate function.
upper.bound	numeric; an upper bound for the coefficient, usually the upper bound of the support of the claim severity mgf.
h	an expression written as a function of x or of x and y, or alternatively the name of a function, giving function $h$ in the Lundberg equation (see below); ignored if mgf.claim is provided.
reinsurance	the type of reinsurance for the portfolio; can be abbreviated.
from, to	the range over which the adjustment coefficient will be calculated.
n	integer; the number of values at which to evaluate the adjustment coefficient.
x	an object of class "adjCoef".
xlab, ylab	label of the x and y axes, respectively.
main	main title.
sub	subtitle, defaulting to the type of reinsurance.
type	1-character string giving the type of plot desired; see plot for details.
add	logical; if TRUE add to already existing plot.
...	further graphical parameters accepted by plot or lines.

Details

In the typical case reinsurance = "none", the coefficient of determination is the smallest (strictly) positive root of the Lundberg equation

$h(x) = E[e^{x B - x c W}] = 1$

on $[0, m)$, where $m =$ upper.bound, $B$ is the claim severity random variable, $W$ is the claim interarrival (or wait) time random variable and $c =$ premium.rate. The premium rate must satisfy the positive safety loading constraint $E[B - c W] < 0$.

With reinsurance = "proportional", the equation becomes

$h(x, y) = E[e^{x y B - x c(y) W}] = 1,$

where $y$ is the retention rate and $c(y)$ is the premium rate function.

With reinsurance = "excess-of-loss", the equation becomes

$h(x, y) = E[e^{x \min(B, y) - x c(y) W}] = 1,$

where $y$ is the retention limit and $c(y)$ is the premium rate function.

One can use argument h as an alternative way to provide function $h(x)$ or $h(x, y)$. This is necessary in cases where random variables $B$ and $W$ are not independent.

The root of $h(x) = 1$ is found by minimizing $(h(x) - 1)^2$.

Value

If reinsurance = "none", a numeric vector of length one. Otherwise, a function of class "adjCoef" inheriting from the "function" class.

Author(s)

Christophe Dutang, Vincent Goulet [email protected]

References

Bowers, N. J. J., Gerber, H. U., Hickman, J., Jones, D. and Nesbitt, C. (1986), Actuarial Mathematics, Society of Actuaries.

Centeno, M. d. L. (2002), Measuring the effects of reinsurance by the adjustment coefficient in the Sparre-Anderson model, Insurance: Mathematics and Economics 30, 37–49.

Gerber, H. U. (1979), An Introduction to Mathematical Risk Theory, Huebner Foundation.

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2008), Loss Models, From Data to Decisions, Third Edition, Wiley.

Examples

## Basic example: no reinsurance, exponential claim severity and wait
## times, premium rate computed with expected value principle and
## safety loading of 20%.
adjCoef(mgfexp, premium = 1.2, upper = 1)

## Same thing, giving function h.
h <- function(x) 1/((1 - x) * (1 + 1.2 * x))
adjCoef(h = h, upper = 1)

## Example 11.4 of Klugman et al. (2008)
mgfx <- function(x) 0.6 * exp(x) + 0.4 * exp(2 * x)
adjCoef(mgfx(x), mgfexp(x, 4), prem = 7, upper = 0.3182)

## Proportional reinsurance, same assumptions as above, reinsurer's
## safety loading of 30%.
mgfx <- function(x, y) mgfexp(x * y)
p <- function(x) 1.3 * x - 0.1
h <- function(x, a) 1/((1 - a * x) * (1 + x * p(a)))
R1 <- adjCoef(mgfx, premium = p, upper = 1, reins = "proportional",
              from = 0, to = 1, n = 11)
R2 <- adjCoef(h = h, upper = 1, reins = "p",
             from = 0, to = 1, n = 101)
R1(seq(0, 1, length = 10))	# evaluation for various retention rates
R2(seq(0, 1, length = 10))	# same
plot(R1)		        # graphical representation
plot(R2, col = "green", add = TRUE) # smoother function

## Excess-of-loss reinsurance
p <- function(x) 1.3 * levgamma(x, 2, 2) - 0.1
mgfx <- function(x, l)
    mgfgamma(x, 2, 2) * pgamma(l, 2, 2 - x) +
    exp(x * l) * pgamma(l, 2, 2, lower = FALSE)
h <- function(x, l) mgfx(x, l) * mgfexp(-x * p(l))
R1 <- adjCoef(mgfx, upper = 1, premium = p, reins = "excess-of-loss",
             from = 0, to = 10, n = 11)
R2 <- adjCoef(h = h, upper = 1, reins = "e",
             from = 0, to = 10, n = 101)
plot(R1)
plot(R2, col = "green", add = TRUE)

---

Aggregate Claim Amount Distribution

Description

Compute the aggregate claim amount cumulative distribution function of a portfolio over a period using one of five methods.

Usage

aggregateDist(method = c("recursive", "convolution", "normal",
                         "npower", "simulation"),
              model.freq = NULL, model.sev = NULL, p0 = NULL,
              x.scale = 1, convolve = 0, moments, nb.simul, ...,
              tol = 1e-06, maxit = 500, echo = FALSE)

## S3 method for class 'aggregateDist'
print(x, ...)

## S3 method for class 'aggregateDist'
plot(x, xlim, ylab = expression(F[S](x)),
     main = "Aggregate Claim Amount Distribution",
     sub = comment(x), ...)

## S3 method for class 'aggregateDist'
summary(object, ...)

## S3 method for class 'aggregateDist'
mean(x, ...)

## S3 method for class 'aggregateDist'
diff(x, ...)

Arguments

method	method to be used
model.freq	for "recursive" method: a character string giving the name of a distribution in the $(a, b, 0)$ or $(a, b, 1)$ families of distributions. For "convolution" method: a vector of claim number probabilities. For "simulation" method: a frequency simulation model (see rcomphierarc for details) or NULL. Ignored with normal and npower methods.
model.sev	for "recursive" and "convolution" methods: a vector of claim amount probabilities. For "simulation" method: a severity simulation model (see rcomphierarc for details) or NULL. Ignored with normal and npower methods.
p0	arbitrary probability at zero for the frequency distribution. Creates a zero-modified or zero-truncated distribution if not NULL. Used only with "recursive" method.
x.scale	value of an amount of 1 in the severity model (monetary unit). Used only with "recursive" and "convolution" methods.
convolve	number of times to convolve the resulting distribution with itself. Used only with "recursive" method.
moments	vector of the true moments of the aggregate claim amount distribution; required only by the "normal" or "npower" methods.
nb.simul	number of simulations for the "simulation" method.
...	parameters of the frequency distribution for the "recursive" method; further arguments to be passed to or from other methods otherwise.
tol	the resulting cumulative distribution in the "recursive" method will get less than tol away from 1.
maxit	maximum number of recursions in the "recursive" method.
echo	logical; echo the recursions to screen in the "recursive" method.
x, object	an object of class "aggregateDist".
xlim	numeric of length 2; the $x$ limits of the plot.
ylab	label of the y axis.
main	main title.
sub	subtitle, defaulting to the calculation method.

Details

aggregateDist returns a function to compute the cumulative distribution function (cdf) of the aggregate claim amount distribution in any point.

The "recursive" method computes the cdf using the Panjer algorithm; the "convolution" method using convolutions; the "normal" method using a normal approximation; the "npower" method using the Normal Power 2 approximation; the "simulation" method using simulations. More details follow.

Value

A function of class "aggregateDist", inheriting from the "function" class when using normal and Normal Power approximations and additionally inheriting from the "ecdf" and "stepfun" classes when other methods are used.

There are methods available to summarize (summary), represent (print), plot (plot), compute quantiles (quantile) and compute the mean (mean) of "aggregateDist" objects.

For the diff method: a numeric vector of probabilities corresponding to the probability mass function evaluated at the knots of the distribution.

Recursive method

The frequency distribution must be a member of the $(a, b, 0)$ or $(a, b, 1)$ families of discrete distributions.

To use a distribution from the $(a, b, 0)$ family, model.freq must be one of "binomial", "geometric", "negative binomial" or "poisson", and p0 must be NULL.

To use a zero-truncated distribution from the $(a, b, 1)$ family, model.freq may be one of the strings above together with p0 = 0. As a shortcut, model.freq may also be one of "zero-truncated binomial", "zero-truncated geometric", "zero-truncated negative binomial", "zero-truncated poisson" or "logarithmic", and p0 is then ignored (with a warning if non NULL).

(Note: since the logarithmic distribution is always zero-truncated. model.freq = "logarithmic" may be used with either p0 = NULL or p0 = 0.)

To use a zero-modified distribution from the $(a, b, 1)$ family, model.freq may be one of standard frequency distributions mentioned above with p0 set to some probability that the distribution takes the value $0$. It is equivalent, but more explicit, to set model.freq to one of "zero-modified binomial", "zero-modified geometric", "zero-modified negative binomial", "zero-modified poisson" or "zero-modified logarithmic".

The parameters of the frequency distribution must be specified using names identical to the arguments of the appropriate function dbinom, dgeom, dnbinom, dpois or dlogarithmic. In the latter case, do take note that the parametrization of dlogarithmic is different from Appendix B of Klugman et al. (2012).

If the length of p0 is greater than one, only the first element is used, with a warning.

model.sev is a vector of the (discretized) claim amount distribution $X$; the first element must be $f_X(0) = \Pr[X = 0]$.

The recursion will fail to start if the expected number of claims is too large. One may divide the appropriate parameter of the frequency distribution by $2^n$ and convolve the resulting distribution $n =$ convolve times.

Failure to obtain a cumulative distribution function less than tol away from 1 within maxit iterations is often due to too coarse a discretization of the severity distribution.

Convolution method

The cumulative distribution function (cdf) $F_S(x)$ of the aggregate claim amount of a portfolio in the collective risk model is

$F_S(x) = \sum_{n = 0}^{\infty} F_X^{*n}(x) p_n,$

for $x = 0, 1, \dots$; $p_n = \Pr[N = n]$ is the frequency probability mass function and $F_X^{*n}(x)$ is the cdf of the $n$th convolution of the (discrete) claim amount random variable.

model.freq is vector $p_n$ of the number of claims probabilities; the first element must be $\Pr[N = 0]$.

model.sev is vector $f_X(x)$ of the (discretized) claim amount distribution; the first element must be $f_X(0)$.

Normal and Normal Power 2 methods

The Normal approximation of a cumulative distribution function (cdf) $F(x)$ with mean $\mu$ and standard deviation $\sigma$ is

$F(x) \approx \Phi\left( \frac{x - \mu}{\sigma} \right).$

The Normal Power 2 approximation of a cumulative distribution function (cdf) $F(x)$ with mean $\mu$, standard deviation $\sigma$ and skewness $\gamma$ is

$F(x) \approx \Phi \left( -\frac{3}{\gamma} + \sqrt{\frac{9}{\gamma^2} + 1 + \frac{6}{\gamma} \frac{x - \mu}{\sigma}} \right).$

This formula is valid only for the right-hand tail of the distribution and skewness should not exceed unity.

Simulation method

This methods returns the empirical distribution function of a sample of size nb.simul of the aggregate claim amount distribution specified by model.freq and model.sev. rcomphierarc is used for the simulation of claim amounts, hence both the frequency and severity models can be mixtures of distributions.

Author(s)

Vincent Goulet [email protected] and Louis-Philippe Pouliot

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

Daykin, C.D., Pentikäinen, T. and Pesonen, M. (1994), Practical Risk Theory for Actuaries, Chapman & Hall.

See Also

discretize to discretize a severity distribution; mean.aggregateDist to compute the mean of the distribution; quantile.aggregateDist to compute the quantiles or the Value-at-Risk; CTE.aggregateDist to compute the Conditional Tail Expectation (or Tail Value-at-Risk); rcomphierarc.

Examples

## Convolution method (example 9.5 of Klugman et al. (2012))
fx <- c(0, 0.15, 0.2, 0.25, 0.125, 0.075,
        0.05, 0.05, 0.05, 0.025, 0.025)
pn <- c(0.05, 0.1, 0.15, 0.2, 0.25, 0.15, 0.06, 0.03, 0.01)
Fs <- aggregateDist("convolution", model.freq = pn,
                    model.sev = fx, x.scale = 25)
summary(Fs)
c(Fs(0), diff(Fs(25 * 0:21))) # probability mass function
plot(Fs)

## Recursive method (example 9.10 of Klugman et al. (2012))
fx <- c(0, crossprod(c(2, 1)/3,
                     matrix(c(0.6, 0.7, 0.4, 0, 0, 0.3), 2, 3)))
Fs <- aggregateDist("recursive", model.freq = "poisson",
                    model.sev = fx, lambda = 3)
plot(Fs)
Fs(knots(Fs))		      # cdf evaluated at its knots
diff(Fs)                      # probability mass function

## Recursive method (high frequency)
fx <- c(0, 0.15, 0.2, 0.25, 0.125, 0.075,
        0.05, 0.05, 0.05, 0.025, 0.025)
## Not run: Fs <- aggregateDist("recursive", model.freq = "poisson",
                    model.sev = fx, lambda = 1000)
## End(Not run)
Fs <- aggregateDist("recursive", model.freq = "poisson",
                    model.sev = fx, lambda = 250, convolve = 2, maxit = 1500)
plot(Fs)

## Recursive method (zero-modified distribution; example 9.11 of
## Klugman et al. (2012))
Fn <- aggregateDist("recursive", model.freq = "binomial",
                    model.sev = c(0.3, 0.5, 0.2), x.scale = 50,
                    p0 = 0.4, size = 3, prob = 0.3)
diff(Fn)

## Equivalent but more explicit call
aggregateDist("recursive", model.freq = "zero-modified binomial",
              model.sev = c(0.3, 0.5, 0.2), x.scale = 50,
              p0 = 0.4, size = 3, prob = 0.3)

## Recursive method (zero-truncated distribution). Using 'fx' above
## would mean that both Pr[N = 0] = 0 and Pr[X = 0] = 0, therefore
## Pr[S = 0] = 0 and recursions would not start.
fx <- discretize(pexp(x, 1), from = 0, to = 100, method = "upper")
fx[1L] # non zero
aggregateDist("recursive", model.freq = "zero-truncated poisson",
              model.sev = fx, lambda = 3, x.scale = 25, echo=TRUE)

## Normal Power approximation
Fs <- aggregateDist("npower", moments = c(200, 200, 0.5))
Fs(210)

## Simulation method
model.freq <- expression(data = rpois(3))
model.sev <- expression(data = rgamma(100, 2))
Fs <- aggregateDist("simulation", nb.simul = 1000,
                    model.freq, model.sev)
mean(Fs)
plot(Fs)

## Evaluation of ruin probabilities using Beekman's formula with
## Exponential(1) claim severity, Poisson(1) frequency and premium rate
## c = 1.2.
fx <- discretize(pexp(x, 1), from = 0, to = 100, method = "lower")
phi0 <- 0.2/1.2
Fs <- aggregateDist(method = "recursive", model.freq = "geometric",
                    model.sev = fx, prob = phi0)
1 - Fs(400)			# approximate ruin probability
u <- 0:100
plot(u, 1 - Fs(u), type = "l", main = "Ruin probability")

---

The “Beta Integral”

Description

The “beta integral” is just a multiple of the non regularized incomplete beta function. This function provides an R interface to the C level routine. It is not exported by the package.

Usage

betaint(x, a, b)

Arguments

x	vector of quantiles.
a, b	parameters. See Details for admissible values.

Details

Function betaint computes the “beta integral”

$B(a, b; x) = \Gamma(a + b) \int_0^x t^{a-1} (1-t)^{b-1} dt$

for $a > 0$, $b \neq -1, -2, \ldots$ and $0 < x < 1$. (Here $\Gamma(\alpha)$ is the function implemented by R's gamma() and defined in its help.) When $b > 0$,

$B(a, b; x) = \Gamma(a) \Gamma(b) I_x(a, b),$

where $I_x(a, b)$ is pbeta(x, a, b). When $b < 0$, $b \neq -1, -2, \ldots$, and $a > 1 + [-b]$,

$\begin{array}{rcl} B(a, b; x) &=& \displaystyle -\Gamma(a + b) \left[ \frac{x^{a-1} (1-x)^b}{b} + \frac{(a-1) x^{a-2} (1-x)^{b+1}}{b (b+1)} \right. \\ & & \displaystyle\left. + \cdots + \frac{(a-1) \cdots (a-r) x^{a-r-1} (1-x)^{b+r}}{b (b+1) \cdots (b+r)} \right] \\ & & \displaystyle + \frac{(a-1) \cdots (a-r-1)}{b (b+1) \cdots (b+r)} \Gamma(a-r-1) \\ & & \times \Gamma(b+r+1) I_x(a-r-1, b+r+1), \end{array}$

where $r = [-b]$.

This function is used (at the C level) to compute the limited expected value for distributions of the transformed beta family; see, for example, levtrbeta.

Value

The value of the integral.

Invalid arguments will result in return value NaN, with a warning.

Note

The need for this function in the package is well explained in the introduction of Appendix A of Klugman et al. (2012). See also chapter 6 and 15 of Abramowitz and Stegun (1972) for definitions and relations to the hypergeometric series.

Author(s)

Vincent Goulet [email protected]

References

Abramowitz, M. and Stegun, I. A. (1972), Handbook of Mathematical Functions, Dover.

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

Examples

x <- 0.3
a <- 7

## case with b > 0
b <- 2
actuar:::betaint(x, a, b)
gamma(a) * gamma(b) * pbeta(x, a, b)    # same

## case with b < 0
b <- -2.2
r <- floor(-b)        # r = 2
actuar:::betaint(x, a, b)

## "manual" calculation
s <- (x^(a-1) * (1-x)^b)/b +
    ((a-1) * x^(a-2) * (1-x)^(b+1))/(b * (b+1)) +
    ((a-1) * (a-2) * x^(a-3) * (1-x)^(b+2))/(b * (b+1) * (b+2))
-gamma(a+b) * s +
    (a-1)*(a-2)*(a-3) * gamma(a-r-1)/(b*(b+1)*(b+2)) *
    gamma(b+r+1)*pbeta(x, a-r-1, b+r+1)

---

Raw and Limited Moments of the Beta Distribution

Description

Raw moments and limited moments for the (central) Beta distribution with parameters shape1 and shape2.

Usage

mbeta(order, shape1, shape2)
levbeta(limit, shape1, shape2, order = 1)

Arguments

order	order of the moment.
limit	limit of the loss variable.
shape1, shape2	positive parameters of the Beta distribution.

Details

The $k$th raw moment of the random variable $X$ is $E[X^k]$ and the $k$th limited moment at some limit $d$ is $E[\min(X, d)^k]$, $k > -\alpha$.

The noncentral beta distribution is not supported.

Value

mbeta gives the $k$th raw moment and levbeta gives the $k$th moment of the limited loss variable.

Invalid arguments will result in return value NaN, with a warning.

Author(s)

Vincent Goulet [email protected] and Mathieu Pigeon

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

See Also

Beta for details on the beta distribution and functions [dpqr]beta.

Examples

mbeta(2, 3, 4) - mbeta(1, 3, 4)^2
levbeta(10, 3, 4, order = 2)

---

The Burr Distribution

Description

Density function, distribution function, quantile function, random generation, raw moments and limited moments for the Burr distribution with parameters shape1, shape2 and scale.

Usage

dburr(x, shape1, shape2, rate = 1, scale = 1/rate,
      log = FALSE)
pburr(q, shape1, shape2, rate = 1, scale = 1/rate,
      lower.tail = TRUE, log.p = FALSE)
qburr(p, shape1, shape2, rate = 1, scale = 1/rate,
      lower.tail = TRUE, log.p = FALSE)
rburr(n, shape1, shape2, rate = 1, scale = 1/rate)
mburr(order, shape1, shape2, rate = 1, scale = 1/rate)
levburr(limit, shape1, shape2, rate = 1, scale = 1/rate,
        order = 1)

Arguments

x, q	vector of quantiles.
p	vector of probabilities.
n	number of observations. If length(n) > 1, the length is taken to be the number required.
shape1, shape2, scale	parameters. Must be strictly positive.
rate	an alternative way to specify the scale.
log, log.p	logical; if TRUE, probabilities/densities $p$ are returned as $\log(p)$.
lower.tail	logical; if TRUE (default), probabilities are $P[X \le x]$, otherwise, $P[X > x]$.
order	order of the moment.
limit	limit of the loss variable.

Details

The Burr distribution with parameters shape1 $= \alpha$, shape2 $= \gamma$ and scale $= \theta$ has density:

$f(x) = \frac{\alpha \gamma (x/\theta)^\gamma}{ x [1 + (x/\theta)^\gamma]^{\alpha + 1}}$

for $x > 0$, $\alpha > 0$, $\gamma > 0$ and $\theta > 0$.

The Burr is the distribution of the random variable

$\theta \left(\frac{X}{1 - X}\right)^{1/\gamma},$

where $X$ has a beta distribution with parameters $1$ and $\alpha$.

The Burr distribution has the following special cases:

• A Loglogistic distribution when shape1 == 1;

• A Paralogistic distribution when shape2 == shape1;

• A Pareto distribution when shape2 == 1.

The $k$th raw moment of the random variable $X$ is $E[X^k]$, $-\gamma < k < \alpha\gamma$.

The $k$th limited moment at some limit $d$ is $E[\min(X, d)^k]$, $k > -\gamma$ and $\alpha - k/\gamma$ not a negative integer.

Value

dburr gives the density, pburr gives the distribution function, qburr gives the quantile function, rburr generates random deviates, mburr gives the $k$th raw moment, and levburr gives the $k$th moment of the limited loss variable.

Invalid arguments will result in return value NaN, with a warning.

Note

levburr computes the limited expected value using betaint.

Distribution also known as the Burr Type XII or Singh-Maddala distribution. See also Kleiber and Kotz (2003) for alternative names and parametrizations.

The "distributions" package vignette provides the interrelations between the continuous size distributions in actuar and the complete formulas underlying the above functions.

Author(s)

Vincent Goulet [email protected] and Mathieu Pigeon

References

Kleiber, C. and Kotz, S. (2003), Statistical Size Distributions in Economics and Actuarial Sciences, Wiley.

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

See Also

dpareto4 for an equivalent distribution with a location parameter.

Examples

exp(dburr(1, 2, 3, log = TRUE))
p <- (1:10)/10
pburr(qburr(p, 2, 3, 2), 2, 3, 2)

## variance
mburr(2, 2, 3, 1) - mburr(1, 2, 3, 1) ^ 2

## case with shape1 - order/shape2 > 0
levburr(10, 2, 3, 1, order = 2)

## case with shape1 - order/shape2 < 0
levburr(10, 1.5, 0.5, 1, order = 2)

---

Moments and Moment Generating Function of the (non-central) Chi-Squared Distribution

Description

Raw moments, limited moments and moment generating function for the chi-squared ($\chi^2$) distribution with df degrees of freedom and optional non-centrality parameter ncp.

Usage

mchisq(order, df, ncp = 0)
levchisq(limit, df, ncp = 0, order = 1)
mgfchisq(t, df, ncp = 0, log= FALSE)

Arguments

order	order of the moment.
limit	limit of the loss variable.
df	degrees of freedom (non-negative, but can be non-integer).
ncp	non-centrality parameter (non-negative).
t	numeric vector.
log	logical; if TRUE, the cumulant generating function is returned.

Details

The $k$th raw moment of the random variable $X$ is $E[X^k]$, the $k$th limited moment at some limit $d$ is $E[\min(X, d)]$ and the moment generating function is $E[e^{tX}]$.

Only integer moments are supported for the non central Chi-square distribution (ncp > 0).

The limited expected value is supported for the centered Chi-square distribution (ncp = 0).

Value

mchisq gives the $k$th raw moment, levchisq gives the $k$th moment of the limited loss variable, and mgfchisq gives the moment generating function in t.

Invalid arguments will result in return value NaN, with a warning.

Author(s)

Christophe Dutang, Vincent Goulet [email protected]

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

Johnson, N. L. and Kotz, S. (1970), Continuous Univariate Distributions, Volume 1, Wiley.

See Also

Chisquare

Examples

mchisq(2, 3, 4)
levchisq(10, 3, order = 2)
mgfchisq(0.25, 3, 2)

---

Credibility Models

Description

Fit the following credibility models: Bühlmann, Bühlmann-Straub, hierarchical, regression (Hachemeister) or linear Bayes.

Usage

cm(formula, data, ratios, weights, subset,
   regformula = NULL, regdata, adj.intercept = FALSE,
   method = c("Buhlmann-Gisler", "Ohlsson", "iterative"),
   likelihood, ...,
   tol = sqrt(.Machine$double.eps), maxit = 100, echo = FALSE)

## S3 method for class 'cm'
print(x, ...)

## S3 method for class 'cm'
predict(object, levels = NULL, newdata, ...)

## S3 method for class 'cm'
summary(object, levels = NULL, newdata, ...)

## S3 method for class 'summary.cm'
print(x, ...)

Arguments

formula	character string "bayes" or an object of class "formula": a symbolic description of the model to be fit. The details of model specification are given below.
data	a matrix or a data frame containing the portfolio structure, the ratios or claim amounts and their associated weights, if any.
ratios	expression indicating the columns of data containing the ratios or claim amounts.
weights	expression indicating the columns of data containing the weights associated with ratios.
subset	an optional logical expression indicating a subset of observations to be used in the modeling process. All observations are included by default.
regformula	an object of class "formula": symbolic description of the regression component (see lm for details). No left hand side is needed in the formula; if present it is ignored. If NULL, no regression is done on the data.
regdata	an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the regression model.
adj.intercept	if TRUE, the intercept of the regression model is located at the barycenter of the regressor instead of the origin.
method	estimation method for the variance components of the model; see Details.
likelihood	a character string giving the name of the likelihood function in one of the supported linear Bayes cases; see Details.
tol	tolerance level for the stopping criteria for iterative estimation method.
maxit	maximum number of iterations in iterative estimation method.
echo	logical; whether to echo the iterative procedure or not.
x, object	an object of class "cm".
levels	character vector indicating the levels to predict or to include in the summary; if NULL all levels are included.
newdata	data frame containing the variables used to predict credibility regression models.
...	parameters of the prior distribution for cm; additional attributes to attach to the result for the predict and summary methods; further arguments to format for the print.summary method; unused for the print method.

Details

cm is the unified front end for credibility models fitting. The function supports hierarchical models with any number of levels (with Bühlmann and Bühlmann-Straub models as special cases) and the regression model of Hachemeister. Usage of cm is similar to lm for these cases. cm can also fit linear Bayes models, in which case usage is much simplified; see the section on linear Bayes below.

When not "bayes", the formula argument symbolically describes the structure of the portfolio in the form $~ terms$. Each term is an interaction between risk factors contributing to the total variance of the portfolio data. Terms are separated by + operators and interactions within each term by :. For a portfolio divided first into sectors, then units and finally contracts, formula would be ~ sector + sector:unit + sector:unit:contract, where sector, unit and contract are column names in data. In general, the formula should be of the form ~ a + a:b + a:b:c + a:b:c:d + ....

If argument regformula is not NULL, the regression model of Hachemeister is fit to the data. The response is usually time. By default, the intercept of the model is located at time origin. If argument adj.intercept is TRUE, the intercept is moved to the (collective) barycenter of time, by orthogonalization of the design matrix. Note that the regression coefficients may be difficult to interpret in this case.

Arguments ratios, weights and subset are used like arguments select, select and subset, respectively, of function subset.

Data does not have to be sorted by level. Nodes with no data (complete lines of NA except for the portfolio structure) are allowed, with the restriction mentioned above.

The print methods use the option deparse.cutoff to control the printing of the call to cm.

Value

Function cm computes the structure parameters estimators of the model specified in formula. The value returned is an object of class cm.

An object of class "cm" is a list with at least the following components:

means	a list containing, for each level, the vector of linearly sufficient statistics.
weights	a list containing, for each level, the vector of total weights.
unbiased	a vector containing the unbiased variance components estimators, or NULL.
iterative	a vector containing the iterative variance components estimators, or NULL.
cred	for multi-level hierarchical models: a list containing, the vector of credibility factors for each level. For one-level models: an array or vector of credibility factors.
nodes	a list containing, for each level, the vector of the number of nodes in the level.
classification	the columns of data containing the portfolio classification structure.
ordering	a list containing, for each level, the affiliation of a node to the node of the level above.

Regression fits have in addition the following components:

adj.models	a list containing, for each node, the credibility adjusted regression model as obtained with lm.fit or lm.wfit.
transition	if adj.intercept is TRUE, a transition matrix from the basis of the orthogonal design matrix to the basis of the original design matrix.
terms	the terms object used.

The method of predict for objects of class "cm" computes the credibility premiums for the nodes of every level included in argument levels (all by default). Result is a list the same length as levels or the number of levels in formula, or an atomic vector for one-level models.

Hierarchical models

The credibility premium at one level is a convex combination between the linearly sufficient statistic of a node and the credibility premium of the level above. (For the first level, the complement of credibility is given to the collective premium.) The linearly sufficient statistic of a node is the credibility weighted average of the data of the node, except at the last level, where natural weights are used. The credibility factor of node $i$ is equal to

$\frac{w_i}{w_i + a/b},$

where $w_i$ is the weight of the node used in the linearly sufficient statistic, $a$ is the average within node variance and $b$ is the average between node variance.

Regression models

The credibility premium of node $i$ is equal to

$y^\prime b_i^a,$

where $y$ is a matrix created from newdata and $b_i^a$ is the vector of credibility adjusted regression coefficients of node $i$. The latter is given by

$b_i^a = Z_i b_i + (I - Z_I) m,$

where $b_i$ is the vector of regression coefficients based on data of node $i$ only, $m$ is the vector of collective regression coefficients, $Z_i$ is the credibility matrix and $I$ is the identity matrix. The credibility matrix of node $i$ is equal to

$A^{-1} (A + s^2 S_i),$

where $S_i$ is the unscaled regression covariance matrix of the node, $s^2$ is the average within node variance and $A$ is the within node covariance matrix.

If the intercept is positioned at the barycenter of time, matrices $S_i$ and $A$ (and hence $Z_i$) are diagonal. This amounts to use Bühlmann-Straub models for each regression coefficient.

Argument newdata provides the “future” value of the regressors for prediction purposes. It should be given as specified in predict.lm.

Variance components estimation

For hierarchical models, two sets of estimators of the variance components (other than the within node variance) are available: unbiased estimators and iterative estimators.

Unbiased estimators are based on sums of squares of the form

$B_i = \sum_j w_{ij} (X_{ij} - \bar{X}_i)^2 - (J - 1) a$

and constants of the form

$c_i = w_i - \sum_j \frac{w_{ij}^2}{w_i},$

where $X_{ij}$ is the linearly sufficient statistic of level $(ij)$; $\bar{X_{i}}$ is the weighted average of the latter using weights $w_{ij}$; $w_i = \sum_j w_{ij}$; $J$ is the effective number of nodes at level $(ij)$; $a$ is the within variance of this level. Weights $w_{ij}$ are the natural weights at the lowest level, the sum of the natural weights the next level and the sum of the credibility factors for all upper levels.

The Bühlmann-Gisler estimators (method = "Buhlmann-Gisler") are given by

$b = \frac{1}{I} \sum_i \max \left( \frac{B_i}{c_i}, 0 \right),$

that is the average of the per node variance estimators truncated at 0.

The Ohlsson estimators (method = "Ohlsson") are given by

$b = \frac{\sum_i B_i}{\sum_i c_i},$

that is the weighted average of the per node variance estimators without any truncation. Note that negative estimates will be truncated to zero for credibility factor calculations.

In the Bühlmann-Straub model, these estimators are equivalent.

Iterative estimators method = "iterative" are pseudo-estimators of the form

$b = \frac{1}{d} \sum_i w_i (X_i - \bar{X})^2,$

where $X_i$ is the linearly sufficient statistic of one level, $\bar{X}$ is the linearly sufficient statistic of the level above and $d$ is the effective number of nodes at one level minus the effective number of nodes of the level above. The Ohlsson estimators are used as starting values.

For regression models, with the intercept at time origin, only iterative estimators are available. If method is different from "iterative", a warning is issued. With the intercept at the barycenter of time, the choice of estimators is the same as in the Bühlmann-Straub model.

Linear Bayes

When formula is "bayes", the function computes pure Bayesian premiums for the following combinations of distributions where they are linear credibility premiums:

• $X|\Theta = \theta \sim \mathrm{Poisson}(\theta)$ and $\Theta \sim \mathrm{Gamma}(\alpha, \lambda)$;

• $X|\Theta = \theta \sim \mathrm{Exponential}(\theta)$ and $\Theta \sim \mathrm{Gamma}(\alpha, \lambda)$;

• $X|\Theta = \theta \sim \mathrm{Gamma}(\tau, \theta)$ and $\Theta \sim \mathrm{Gamma}(\alpha, \lambda)$;

• $X|\Theta = \theta \sim \mathrm{Normal}(\theta, \sigma_2^2)$ and $\Theta \sim \mathrm{Normal}(\mu, \sigma_1^2)$;

• $X|\Theta = \theta \sim \mathrm{Bernoulli}(\theta)$ and $\Theta \sim \mathrm{Beta}(a, b)$;

• $X|\Theta = \theta \sim \mathrm{Binomial}(\nu, \theta)$ and $\Theta \sim \mathrm{Beta}(a, b)$;

• $X|\Theta = \theta \sim \mathrm{Geometric}(\theta)$ and $\Theta \sim \mathrm{Beta}(a, b)$.

• $X|\Theta = \theta \sim \mathrm{Negative~Binomial}(r, \theta)$ and $\Theta \sim \mathrm{Beta}(a, b)$.

The following combination is also supported: $X|\Theta = \theta \sim \mathrm{Single~Parameter~Pareto}(\theta)$ and $\Theta \sim \mathrm{Gamma}(\alpha, \lambda)$. In this case, the Bayesian estimator not of the risk premium, but rather of parameter $\theta$ is linear with a “credibility” factor that is not restricted to $(0, 1)$.

Argument likelihood identifies the distribution of $X|\Theta = \theta$ as one of "poisson", "exponential", "gamma", "normal", "bernoulli", "binomial", "geometric", "negative binomial" or "pareto".

The parameters of the distributions of $X|\Theta = \theta$ (when needed) and $\Theta$ are set in ... using the argument names (and default values) of dgamma, dnorm, dbeta, dbinom, dnbinom or dpareto1, as appropriate. For the Gamma/Gamma case, use shape.lik for the shape parameter $\tau$ of the Gamma likelihood. For the Normal/Normal case, use sd.lik for the standard error $\sigma_2$ of the Normal likelihood.

Data for the linear Bayes case may be a matrix or data frame as usual; an atomic vector to fit the model to a single contract; missing or NULL to fit the prior model. Arguments ratios, weights and subset are ignored.

Author(s)

Vincent Goulet [email protected], Xavier Milhaud, Tommy Ouellet, Louis-Philippe Pouliot

References

Bühlmann, H. and Gisler, A. (2005), A Course in Credibility Theory and its Applications, Springer.

Belhadj, H., Goulet, V. and Ouellet, T. (2009), On parameter estimation in hierarchical credibility, Astin Bulletin 39.

Goulet, V. (1998), Principles and application of credibility theory, Journal of Actuarial Practice 6, ISSN 1064-6647.

Goovaerts, M. J. and Hoogstad, W. J. (1987), Credibility Theory, Surveys of Actuarial Studies, No. 4, Nationale-Nederlanden N.V.

See Also

subset, formula, lm, predict.lm.

Examples

data(hachemeister)

## Buhlmann-Straub model
fit <- cm(~state, hachemeister,
          ratios = ratio.1:ratio.12, weights = weight.1:weight.12)
fit				# print method
predict(fit)			# credibility premiums
summary(fit)			# more details

## Two-level hierarchical model. Notice that data does not have
## to be sorted by level
X <- data.frame(unit = c("A", "B", "A", "B", "B"), hachemeister)
fit <- cm(~unit + unit:state, X, ratio.1:ratio.12, weight.1:weight.12)
predict(fit)
predict(fit, levels = "unit")	# unit credibility premiums only
summary(fit)
summary(fit, levels = "unit")	# unit summaries only

## Regression model with intercept at time origin
fit <- cm(~state, hachemeister,
          regformula = ~time, regdata = data.frame(time = 12:1),
          ratios = ratio.1:ratio.12, weights = weight.1:weight.12)
fit
predict(fit, newdata = data.frame(time = 0))
summary(fit, newdata = data.frame(time = 0))

## Same regression model, with intercept at barycenter of time
fit <- cm(~state, hachemeister, adj.intercept = TRUE,
          regformula = ~time, regdata = data.frame(time = 12:1),
          ratios = ratio.1:ratio.12, weights = weight.1:weight.12)
fit
predict(fit, newdata = data.frame(time = 0))
summary(fit, newdata = data.frame(time = 0))

## Poisson/Gamma pure Bayesian model
fit <- cm("bayes", data = c(5, 3, 0, 1, 1),
          likelihood = "poisson", shape = 3, rate = 3)
fit
predict(fit)
summary(fit)

## Normal/Normal pure Bayesian model
cm("bayes", data = c(5, 3, 0, 1, 1),
   likelihood = "normal", sd.lik = 2,
   mean = 2, sd = 1)

---

Density and Cumulative Distribution Function for Modified Data

Description

Compute probability density function or cumulative distribution function of the payment per payment or payment per loss random variable under any combination of the following coverage modifications: deductible, limit, coinsurance, inflation.

Usage

coverage(pdf, cdf, deductible = 0, franchise = FALSE,
         limit = Inf, coinsurance = 1, inflation = 0,
         per.loss = FALSE)

Arguments

pdf, cdf	function object or character string naming a function to compute, respectively, the probability density function and cumulative distribution function of a probability law.
deductible	a unique positive numeric value.
franchise	logical; TRUE for a franchise deductible, FALSE (default) for an ordinary deductible.
limit	a unique positive numeric value larger than deductible.
coinsurance	a unique value between 0 and 1; the proportion of coinsurance.
inflation	a unique value between 0 and 1; the rate of inflation.
per.loss	logical; TRUE for the per loss distribution, FALSE (default) for the per payment distribution.

Details

coverage returns a function to compute the probability density function (pdf) or the cumulative distribution function (cdf) of the distribution of losses under coverage modifications. The pdf and cdf of unmodified losses are pdf and cdf, respectively.

If pdf is specified, the pdf is returned; if pdf is missing or NULL, the cdf is returned. Note that cdf is needed if there is a deductible or a limit.

Value

An object of mode "function" with the same arguments as pdf or cdf, except "lower.tail", "log.p" and "log", which are not supported.

Note

Setting arguments of the function returned by coverage using formals may very well not work as expected.

Author(s)

Vincent Goulet [email protected] and Mathieu Pigeon

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

See Also

vignette("coverage") for the exact definitions of the per payment and per loss random variables under an ordinary or franchise deductible.

Examples

## Default case: pdf of the per payment random variable with
## an ordinary deductible
coverage(dgamma, pgamma, deductible = 1)

## Add a limit
f <- coverage(dgamma, pgamma, deductible = 1, limit = 7)
f <- coverage("dgamma", "pgamma", deductible = 1, limit = 7) # same
f(0, shape = 3, rate = 1)
f(2, shape = 3, rate = 1)
f(6, shape = 3, rate = 1)
f(8, shape = 3, rate = 1)
curve(dgamma(x, 3, 1), xlim = c(0, 10), ylim = c(0, 0.3))    # original
curve(f(x, 3, 1), xlim = c(0.01, 5.99), col = 4, add = TRUE) # modified
points(6, f(6, 3, 1), pch = 21, bg = 4)

## Cumulative distribution function
F <- coverage(cdf = pgamma, deductible = 1, limit = 7)
F(0, shape = 3, rate = 1)
F(2, shape = 3, rate = 1)
F(6, shape = 3, rate = 1)
F(8, shape = 3, rate = 1)
curve(pgamma(x, 3, 1), xlim = c(0, 10), ylim = c(0, 1))    # original
curve(F(x, 3, 1), xlim = c(0, 5.99), col = 4, add = TRUE)  # modified
curve(F(x, 3, 1), xlim = c(6, 10), col = 4, add = TRUE)    # modified

## With no deductible, all distributions below are identical
coverage(dweibull, pweibull, limit = 5)
coverage(dweibull, pweibull, per.loss = TRUE, limit = 5)
coverage(dweibull, pweibull, franchise = TRUE, limit = 5)
coverage(dweibull, pweibull, per.loss = TRUE, franchise = TRUE,
         limit = 5)

## Coinsurance alone; only case that does not require the cdf
coverage(dgamma, coinsurance = 0.8)

---

Conditional Tail Expectation

Description

Conditional Tail Expectation, also called Tail Value-at-Risk.

TVaR is an alias for CTE.

Usage

CTE(x, ...)

## S3 method for class 'aggregateDist'
CTE(x, conf.level = c(0.9, 0.95, 0.99),
         names = TRUE, ...)

TVaR(x, ...)

Arguments

x	an R object.
conf.level	numeric vector of probabilities with values in $[0, 1)$.
names	logical; if true, the result has a names attribute. Set to FALSE for speedup with many probs.
...	further arguments passed to or from other methods.

Details

The Conditional Tail Expectation (or Tail Value-at-Risk) measures the average of losses above the Value at Risk for some given confidence level, that is $E[X|X > \mathrm{VaR}(X)]$ where $X$ is the loss random variable.

CTE is a generic function with, currently, only a method for objects of class "aggregateDist".

For the recursive, convolution and simulation methods of aggregateDist, the CTE is computed from the definition using the empirical cdf.

For the normal approximation method, an explicit formula exists:

$\mu + \frac{\sigma}{(1 - \alpha) \sqrt{2 \pi}} e^{-\mathrm{VaR}(X)^2/2},$

where $\mu$ is the mean, $\sigma$ the standard deviation and $\alpha$ the confidence level.

For the Normal Power approximation, the explicit formula given in Castañer et al. (2013) is

$\mu + \frac{\sigma}{(1 - \alpha) \sqrt{2 \pi}} e^{-\mathrm{VaR}(X)^2/2} \left( 1 + \frac{\gamma}{6} \mathrm{VaR}(X) \right),$

where, as above, $\mu$ is the mean, $\sigma$ the standard deviation, $\alpha$ the confidence level and $\gamma$ is the skewness.

Value

A numeric vector, named if names is TRUE.

Author(s)

Vincent Goulet [email protected] and Tommy Ouellet

References

Castañer, A. and Claramunt, M.M. and Mármol, M. (2013), Tail value at risk. An analysis with the Normal-Power approximation. In Statistical and Soft Computing Approaches in Insurance Problems, pp. 87-112. Nova Science Publishers, 2013. ISBN 978-1-62618-506-7.

See Also

aggregateDist; VaR

Examples

model.freq <- expression(data = rpois(7))
model.sev <- expression(data = rnorm(9, 2))
Fs <- aggregateDist("simulation", model.freq, model.sev, nb.simul = 1000)
CTE(Fs)

---

Individual Dental Claims Data Set

Description

Basic dental claims on a policy with a deductible of 50.

Usage

dental

Format

A vector containing 10 observations

Source

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (1998), Loss Models, From Data to Decisions, Wiley.

---

Discretization of a Continuous Distribution

Description

Compute a discrete probability mass function from a continuous cumulative distribution function (cdf) with various methods.

discretise is an alias for discretize.

Usage

discretize(cdf, from, to, step = 1,
           method = c("upper", "lower", "rounding", "unbiased"),
           lev, by = step, xlim = NULL)

discretise(cdf, from, to, step = 1,
           method = c("upper", "lower", "rounding", "unbiased"),
           lev, by = step, xlim = NULL)

Arguments

cdf	an expression written as a function of x, or alternatively the name of a function, giving the cdf to discretize.
from, to	the range over which the function will be discretized.
step	numeric; the discretization step (or span, or lag).
method	discretization method to use.
lev	an expression written as a function of x, or alternatively the name of a function, to compute the limited expected value of the distribution corresponding to cdf. Used only with the "unbiased" method.
by	an alias for step.
xlim	numeric of length 2; if specified, it serves as default for c(from, to).

Details

Usage is similar to curve.

discretize returns the probability mass function (pmf) of the random variable obtained by discretization of the cdf specified in cdf.

Let $F(x)$ denote the cdf, $E[\min(X, x)]$ the limited expected value at $x$, $h$ the step, $p_x$ the probability mass at $x$ in the discretized distribution and set $a =$ from and $b =$ to.

Method "upper" is the forward difference of the cdf $F$:

$p_x = F(x + h) - F(x)$

for $x = a, a + h, \dots, b - step$.

Method "lower" is the backward difference of the cdf $F$:

$p_x = F(x) - F(x - h)$

for $x = a + h, \dots, b$ and $p_a = F(a)$.

Method "rounding" has the true cdf pass through the midpoints of the intervals $[x - h/2, x + h/2)$:

$p_x = F(x + h/2) - F(x - h/2)$

for $x = a + h, \dots, b - step$ and $p_a = F(a + h/2)$. The function assumes the cdf is continuous. Any adjusment necessary for discrete distributions can be done via cdf.

Method "unbiased" matches the first moment of the discretized and the true distributions. The probabilities are as follows:

$p_a = \frac{E[\min(X, a)] - E[\min(X, a + h)]}{h} + 1 - F(a)$

$p_x = \frac{2 E[\min(X, x)] - E[\min(X, x - h)] - E[\min(X, x + h)]}{h}, \quad a < x < b$

$p_b = \frac{E[\min(X, b)] - E[\min(X, b - h)]}{h} - 1 + F(b),$

Value

A numeric vector of probabilities suitable for use in aggregateDist.

Author(s)

Vincent Goulet [email protected]

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

See Also

aggregateDist

Examples

x <- seq(0, 5, 0.5)

op <- par(mfrow = c(1, 1), col = "black")

## Upper and lower discretization
fu <- discretize(pgamma(x, 1), method = "upper",
                 from = 0, to = 5, step = 0.5)
fl <- discretize(pgamma(x, 1), method = "lower",
                 from = 0, to = 5, step = 0.5)
curve(pgamma(x, 1), xlim = c(0, 5))
par(col = "blue")
plot(stepfun(head(x, -1), diffinv(fu)), pch = 19, add = TRUE)
par(col = "green")
plot(stepfun(x, diffinv(fl)), pch = 19, add = TRUE)
par(col = "black")

## Rounding (or midpoint) discretization
fr <- discretize(pgamma(x, 1), method = "rounding",
                 from = 0, to = 5, step = 0.5)
curve(pgamma(x, 1), xlim = c(0, 5))
par(col = "blue")
plot(stepfun(head(x, -1), diffinv(fr)), pch = 19, add = TRUE)
par(col = "black")

## First moment matching
fb <- discretize(pgamma(x, 1), method = "unbiased",
                 lev = levgamma(x, 1), from = 0, to = 5, step = 0.5)
curve(pgamma(x, 1), xlim = c(0, 5))
par(col = "blue")
plot(stepfun(x, diffinv(fb)), pch = 19, add = TRUE)

par(op)

---

Empirical Limited Expected Value

Description

Compute the empirical limited expected value for individual or grouped data.

Usage

elev(x, ...)

## Default S3 method:
elev(x, ...)

## S3 method for class 'grouped.data'
elev(x, ...)

## S3 method for class 'elev'
print(x, digits = getOption("digits") - 2, ...)

## S3 method for class 'elev'
summary(object, ...)

## S3 method for class 'elev'
knots(Fn, ...)

## S3 method for class 'elev'
plot(x, ..., main = NULL, xlab = "x", ylab = "Empirical LEV")

Arguments

x	a vector or an object of class "grouped.data" (in which case only the first column of frequencies is used); for the methods, an object of class "elev", typically.
digits	number of significant digits to use, see print.
Fn, object	an R object inheriting from "ogive".
main	main title.
xlab, ylab	labels of x and y axis.
...	arguments to be passed to subsequent methods.

Details

The limited expected value (LEV) at $u$ of a random variable $X$ is $E[X \wedge u] = E[\min(X, u)]$. For individual data $x_1, \dots, x_n$, the empirical LEV $E_n[X \wedge u]$ is thus

$E_n[X \wedge u] = \frac{1}{n} \left( \sum_{x_j < u} x_j + \sum_{x_j \geq u} u \right).$

Methods of elev exist for individual data or for grouped data created with grouped.data. The formula in this case is too long to show here. See the reference for details.

Value

For elev, a function of class "elev", inheriting from the "function" class.

Author(s)

Vincent Goulet [email protected] and Mathieu Pigeon

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (1998), Loss Models, From Data to Decisions, Wiley.

See Also

grouped.data to create grouped data objects; stepfun for related documentation (even though the empirical LEV is not a step function).

Examples

data(gdental)
lev <- elev(gdental)
lev
summary(lev)
knots(lev)            # the group boundaries

lev(knots(lev))       # empirical lev at boundaries
lev(c(80, 200, 2000)) # and at other limits

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

---

Empirical Moments

Description

Raw empirical moments for individual and grouped data.

Usage

emm(x, order = 1, ...)

## Default S3 method:
emm(x, order = 1, ...)

## S3 method for class 'grouped.data'
emm(x, order = 1, ...)

Arguments

x	a vector or matrix of individual data, or an object of class "grouped data".
order	order of the moment. Must be positive.
...	further arguments passed to or from other methods.

Details

Arguments ... are passed to colMeans; na.rm = TRUE may be useful for individual data with missing values.

For individual data, the $k$th empirical moment is $\sum_{j = 1}^n x_j^k$.

For grouped data with group boundaries $c_0, c_1, \dots, c_r$ and group frequencies $n_1, \dots, n_r$, the $k$th empirical moment is

$\frac{1}{n} \sum_{j = 1}^r \frac{n_j (c_j^{k + 1} - c_{j - 1}^{k + 1})}{ (k + 1) (c_j - c_{j - 1})},$

where $n = \sum_{j = 1}^r n_j$.

Value

A named vector or matrix of moments.

Author(s)

Vincent Goulet [email protected] and Mathieu Pigeon

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (1998), Loss Models, From Data to Decisions, Wiley.

See Also

mean and mean.grouped.data for simpler access to the first moment.

Examples

## Individual data
data(dental)
emm(dental, order = 1:3)

## Grouped data
data(gdental)
emm(gdental)
x <- grouped.data(cj = gdental[, 1],
                  nj1 = sample(1:100, nrow(gdental)),
                  nj2 = sample(1:100, nrow(gdental)))
emm(x) # same as mean(x)

---

Moments and Moment Generating Function of the Exponential Distribution

Description

Raw moments, limited moments and moment generating function for the exponential distribution with rate rate (i.e., mean 1/rate).

Usage

mexp(order, rate = 1)
levexp(limit, rate = 1, order = 1)
mgfexp(t, rate = 1, log = FALSE)

Arguments

order	order of the moment.
limit	limit of the loss variable.
rate	vector of rates.
t	numeric vector.
log	logical; if TRUE, the cumulant generating function is returned.

Details

The $k$th raw moment of the random variable $X$ is $E[X^k]$, the $k$th limited moment at some limit $d$ is $E[\min(X, d)^k]$ and the moment generating function is $E[e^{tX}]$, $k > -1$.

Value

mexp gives the $k$th raw moment, levexp gives the $k$th moment of the limited loss variable, and mgfexp gives the moment generating function in t.

Invalid arguments will result in return value NaN, with a warning.

Author(s)

Vincent Goulet [email protected], Christophe Dutang and Mathieu Pigeon.

References

Johnson, N. L. and Kotz, S. (1970), Continuous Univariate Distributions, Volume 1, Wiley.

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2012), Loss Models, From Data to Decisions, Fourth Edition, Wiley.

See Also

Exponential

Examples

mexp(2, 3) - mexp(1, 3)^2
levexp(10, 3, order = 2)
mgfexp(1,2)

---

Extract or Replace Parts of a Grouped Data Object

Description

Extract or replace subsets of grouped data objects.

Usage

## S3 method for class 'grouped.data'
x[i, j]
## S3 replacement method for class 'grouped.data'
x[i, j] <- value

Arguments

x	an object of class grouped.data.
i, j	elements to extract or replace. i, j are numeric or character or, for [ only, empty. Numeric values are coerced to integer as if by as.integer. For replacement by [, a logical matrix is allowed, but not replacement in the group boundaries and group frequencies simultaneously.
value	a suitable replacement value.

Details

Objects of class "grouped.data" can mostly be indexed like data frames, with the following restrictions:

1. For [, the extracted object must keep a group boundaries column and at least one group frequencies column to remain of class "grouped.data";

2. For [<-, it is not possible to replace group boundaries and group frequencies simultaneously;

3. When replacing group boundaries, length(value) == length(i) + 1.

x[, 1] will return the plain vector of group boundaries.

Replacement of non adjacent group boundaries is not possible for obvious reasons.

Otherwise, extraction and replacement should work just like for data frames.

Value

For [ an object of class "grouped.data", a data frame or a vector.

For [<- an object of class "grouped.data".

Note

Currently [[, [[<-, $ and $<- are not specifically supported, but should work as usual on group frequency columns.

Author(s)

Vincent Goulet [email protected]

See Also

[.data.frame for extraction and replacement methods of data frames, grouped.data to create grouped data objects.

Examples

data(gdental)

(x <- gdental[1])         # select column 1
class(x)                  # no longer a grouped.data object
class(gdental[2])         # same
gdental[, 1]              # group boundaries
gdental[, 2]              # group frequencies

gdental[1:4,]             # a subset
gdental[c(1, 3, 5),]      # avoid this

gdental[1:2, 1] <- c(0, 30, 60) # modified boundaries
gdental[, 2] <- 10              # modified frequencies
## Not run: gdental[1, ] <- 2   # not allowed

---