actuar : An R Package for Actuarial Science

actuar is a package providing additional Actuarial Science functionality to the R statistical system. The project was launched in 2005 and the package is available on the Comprehensive R Archive Network since February 2006. The current version of the package contains functions for use in the ﬁelds of loss distributions modeling, risk theory (including ruin theory), simulation of compound hierarchical models and credibility theory. This paper presents in detail but with few technical terms the most recent version of the package.


Introduction
actuar is a package providing additional Actuarial Science functionality to the R statistical system (R Development Core Team 2008). Various packages on the Comprehensive R Archive Network (CRAN, http://CRAN.R-project.org/) provide functions useful to actuaries, e.g., copula (Yan 2007), Rmetrics (Wuertz 2007), SuppDists (Wheeler 2008) or the R recommended packages nlme (Pinheiro, Bates, DebRoy, Sarkar, and the R Development Core Team 2007) and survival (Lumley 2008) just to name a few. However, actuar aims to serve as a central location for more specifically actuarial functions and data sets. The project was officially launched in 2005 and is under active development.
The feature set of the package can be split in four main categories: loss distributions modeling, risk theory (including ruin theory), simulation of compound hierarchical models and credibility theory. As much as possible, the developers have tried to keep the "user interface" of the various functions of the package consistent. Moreover, the package follows the general R philosophy of working with model objects.
This paper reviews the various features of version 0.9-7 of actuar. We provide enough actuarial background where needed to make the paper self-contained, but otherwise give numerous references for the reader interested to dive more into the subject.
Future versions of the package can be obtained from CRAN at http://CRAN.R-project. org/package=actuar.

Loss distributions modeling
One important task of actuaries is the modeling of claim amount distributions for ratemaking, loss reserving or other risk evaluation purposes. Package actuar offers many functions for loss distributions modeling. The present section details the following actuar features: 1. introduction of 18 additional probability laws and utility functions to get raw moments, limited moments and the moment generating function; 2. fairly extensive support of grouped data; 3. calculation of the empirical raw and limited moments; 4. minimum distance estimation using three different measures; 5. treatment of coverage modifications (deductibles, limits, inflation, coinsurance).

Probability laws
R already includes functions to compute the probability density function (pdf), the cumulative distribution function (cdf) and the quantile function of a fair number of probability laws, as well as functions to generate variates from these laws. For some root foo , the utility functions are named dfoo , pfoo , qfoo and rfoo , respectively.
The actuar package provides d, p, q and r functions for all the probability laws useful for loss severity modeling found in Appendix A of Klugman, Panjer, and Willmot (2004) and not already present in base R, excluding the inverse Gaussian and log-t but including the loggamma distribution (Hogg and Klugman 1984). Table 1 lists the supported distributions as named in Klugman et al. (2004) along with the root names of the R functions.
In addition to the d, p, q and r functions, the package provides m, lev and mgf functions to compute, respectively, theoretical raw moments theoretical limited moments and the moment generating function when it exists. Every probability law of Table 1 is supported, plus the following ones: beta, exponential, chi-square, gamma, lognormal, normal (no lev), uniform and Weibull of base R and the inverse Gaussian distribution of the package SuppDists (Wheeler 2008). The m and lev functions are especially useful with estimation methods based on the matching of raw or limited moments; see Section 2.4 for their empirical counterparts. The mgf functions come in handy to compute the adjustment coefficient in ruin theory; see Section 3.5.
In addition to the 17 distributions of Table 1, the package provides support for a family of distributions deserving a separate presentation. Phase-type distributions (Neuts 1981) are defined as the distribution of the time until absorption of continuous time, finite state Markov processes with m transient states and one absorbing state. Let be the transition rates (or intensity) matrix of such a process and let (π, π m+1 ) be the initial probability vector. Here, T is an m × m non-singular matrix with t ii < 0 for i = 1, . . . , m and t ij ≥ 0 for i = j, t = −T e and e is a column vector with all components equal to 1. Then the cdf of the time until absorption random variable with parameters π and T is where is the matrix exponential of matrix M .
The exponential, the Erlang (gamma with integer shape parameter) and discrete mixtures thereof are common special cases of phase-type distributions.
The package provides functions {d,p,r,m,mgf}phtype for phase-type distributions. Function pphtype is central to the evaluation of ruin probabilities; see Section 3.6. The core of all the functions presented in this subsection is written in C for speed. The matrix exponential C routine is based on expm() from the package Matrix (Bates and Maechler 2008).

Grouped data
What is commonly referred to in Actuarial Science as grouped data is data represented in an interval-frequency manner. In insurance applications, a grouped data set will typically report that there were n j claims in the interval (c j−1 , c j ], j = 1, . . . , r (with the possibility that c r = ∞). This representation is much more compact than an individual data set -where the value of each claim is known -but it also carries far less information. Now that storage space in computers has almost become a non issue, grouped data has somewhat fallen out of fashion.
Still, grouped data remains in use in some fields of actuarial practice and also of interest in teaching. For this reason, actuar provides facilities to store, manipulate and summarize grouped data. A standard storage method is needed since there are many ways to represent grouped data in the computer: using a list or a matrix, aligning the n j s with the c j−1 s or with the c j s, omitting c 0 or not, etc. Moreover, with appropriate extraction, replacement and summary methods, manipulation of grouped data becomes similar to that of individual data.
First, function grouped.data creates a grouped data object similar to -and inheriting from -a data frame. The input of the function is a vector of group boundaries c 0 , c 1 , . . . , c r and one or more vectors of group frequencies n 1 , . . . , n r . Note that there should be one group boundary more than group frequencies. Furthermore, the function assumes that the intervals are contiguous. For example, the following data Notice how extraction results in a simple vector or matrix if either of the group boundaries or the group frequencies are dropped.
As for replacement operations, the package implements the following.
(i) Replacement of one or more group frequencies: It is not possible to replace the boundaries and the frequencies simultaneously.

Histogram of x[, −3]
x is made simple with a method for the mean function:
The R function hist splits individual data into groups and draws an histogram of the frequency distribution. The package introduces a method for already grouped data. Only the first frequencies column is considered (see Figure 1 for the resulting graph): R has a function ecdf to compute the empirical cdf of an individual data set, where I{A} = 1 if A is true and I{A} = 0 otherwise. The function returns a "function" object to compute the value of F n (x) in any x.
The approximation of the empirical cdf for grouped data is called an ogive (Klugman, Panjer, and Willmot 1998;Hogg and Klugman 1984). It is obtained by joining the known values of F n (x) at group boundaries with straight line segments: The package includes a function ogive that otherwise behaves exactly like ecdf. In particular, methods for functions knots and plot allow, respectively, to obtain the knots c 0 , c 1 , . . . , c r of the ogive and a graph (see Figure 2):

Data sets
This is certainly not the most spectacular feature of actuar, but it remains useful for illustrations and examples: the package includes the individual dental claims and grouped dental claims data of Klugman et al. (2004): R> data("dental") R> dental

Calculation of empirical moments
The package provides two functions useful for estimation based on moments. First, function emm computes the kth empirical moment of a sample, whether in individual or grouped data form: R> emm(dental, order = 1:3) [1] 3.355e+02 2.931e+05 3.729e+08 R> emm(gdental, order = 1:3) [1] 3.533e+02 3.577e+05 6.586e+08 Second, in the same spirit as ecdf and ogive, function elev returns a function to compute the empirical limited expected value -or first limited moment -of a sample for any limit. Again, there are methods for individual and grouped data (see Figure 3 for the graphs): R> lev <-elev(dental) R> lev(knots(lev)) [

Minimum distance estimation
Two methods are widely used by actuaries to fit models to data: maximum likelihood and minimum distance. The first technique applied to individual data is well covered by function fitdistr of the package MASS (Venables and Ripley 2002). The second technique minimizes a chosen distance function between theoretical and empirical distributions. Package actuar provides function mde, very similar in usage and inner working to fitdistr, to fit models according to any of the following three distance minimization methods.
1. The Cramér-von Mises method (CvM) minimizes the squared difference between the theoretical cdf and the empirical cdf or ogive at their knots: for individual data and for grouped data. Here, F (x) is the theoretical cdf of a parametric family, F n (x) is the empirical cdf,F n (x) is the ogive and w 1 ≥ 0, w 2 ≥ 0, . . . are arbitrary weights (defaulting to 1).
2. The modified chi-square method (chi-square) applies to grouped data only and minimizes the squared difference between the expected and observed frequency within each group: where n = r j=1 n j . By default, w j = n −1 j .
3. The layer average severity method (LAS) applies to grouped data only and minimizes the squared difference between the theoretical and empirical limited expected value within each group: where is the empirical limited expected value for grouped data.
The arguments of mde are a data set, a function to compute F (x) or E[X ∧ x], starting values for the optimization procedure and the name of the method to use. The empirical functions are computed with ecdf, ogive or elev.
The expressions below fit an exponential distribution to the grouped dental data set, as per Example 2.21 of Klugman et al. (1998):

Coverage modifications
Let X be the random variable of the actual claim amount for an insurance policy, Y L be the random variable of the amount paid per loss and Y P be the random variable of the amount paid per payment. The terminology for the last two random variables refers to whether or not the insurer knows that a loss occurred. Now, the random variables X, Y L and Y P will differ if any of the following coverage modifications are present for the policy: an ordinary or a franchise deductible, a limit, coinsurance or inflation adjustment (see Klugman et al. 2004, Chapter 5 for precise definitions of these terms). Table 2 summarizes the definitions of Y L and Y P .
Often, one will want to use data Y L 1 , . . . , Y L n (or Y P 1 , . . . , Y P n ) from the random variable Y L (Y P ) to fit a model on the unobservable random variable X. This requires expressing the pdf or cdf of Y L (Y P ) in terms of the pdf or cdf of X. Function coverage of actuar does just that: given a pdf or cdf and any combination of the coverage modifications mentioned above, coverage returns a function object to compute the pdf or cdf of the modified random variable. The function can then be used in modeling like any other dfoo or pfoo function.

Risk theory
Risk theory refers to a body of techniques to model and measure the risk associated with a portfolio of insurance contracts. A first approach consists in modeling the distribution of total claims over a fixed period of time using the classical collective model of risk theory. A second input of interest to the actuary is the evolution of the surplus of the insurance company over many periods of time. In ruin theory, the main quantity of interest is the probability that the surplus becomes negative, in which case technical ruin of the insurance company occurs.
The current version of actuar contains four visible functions related to the above problems: two for the calculation of the aggregate claim amount distribution and two for ruin probability calculations.
We briefly expose the underlying models before we introduce each set of functions.

The collective risk model
Let random variable S represent the aggregate claim amount (or total amount of claims) of a portfolio of independent risks over a fixed period of time, random variable N represent the number of claims (or frequency) in the portfolio over that period, and random variable C j represent the amount of claim j (or severity). Then, we have the random sum where we assume that C 1 , C 2 , . . . are mutually independent and identically distributed random variables each independent of N . The task at hand consists in evaluating numerically the cdf of S, given by where F C (x) = P[C ≤ x] is the common cdf of C 1 , . . . , C n , p n = P[N = n] and F * n C (x) = P[C 1 + · · · + C n ≤ x] is the n-fold convolution of F C (·). If C is discrete on 0, 1, 2, . . . , one has

Discretization of claim amount distributions
Some numerical techniques to compute the aggregate claim amount distribution (see Section 3.3) require a discrete arithmetic claim amount distribution; that is, a distribution defined on 0, h, 2h, . . . for some step (or span, or lag) h. The package provides function discretize to discretize a continuous distribution. (The function can also be used to modify the support of an already discrete distribution, but this requires additional care.) Let F (x) denote the cdf of the distribution to discretize on some interval (a, b) and f x denote the probability mass at x in the discretized distribution. Currently, discretize supports the following four discretization methods.
1. Upper discretization, or forward difference of F (x): The discretized cdf is always above the true cdf.
2. Lower discretization, or backward difference of F (x): The discretized cdf is always under the true cdf.
3. Rounding of the random variable, or the midpoint method: The true cdf passes exactly midway through the steps of the discretized cdf.
4. Unbiased, or local matching of the first moment method: The discretized and the true distributions have the same total probability and expected value on (a, b). Usage of discretize is similar to R's plotting function curve. The cdf to discretize and, for the unbiased method only, the limited expected value function are passed to discretize as expressions in x. The other arguments are the upper and lower bounds of the discretization interval, the step h and the discretization method. For example, upper and unbiased discretizations of a Gamma(2, 1) distribution on (0, 17) with a step of 0.5 are achieved with, respectively, R> fx <-discretize(pgamma(x, 2, 1), method = "upper", + from = 0, to = 17, step = 0.5) R> fx <-discretize(pgamma(x, 2, 1), method = "unbiased", + lev = levgamma(x, 2, 1), from = 0, to = 17, step = 0.5) Function discretize is written in a modular fashion making it simple to add other discretization methods if needed.

Calculation of the aggregate claim amount distribution
Function aggregateDist serves as a unique front end for various methods to compute or approximate the cdf of the aggregate claim amount random variable S. Currently, five methods are supported.
1. Recursive calculation using the algorithm of Panjer (1981). This requires the severity distribution to be discrete arithmetic on 0, 1, 2, . . . , m for some monetary unit and the frequency distribution to be a member of either the (a, b, 0) or (a, b, 1) family of distributions (Klugman et al. 2004). (These families contain the Poisson, binomial, negative binomial and logarithmic distributions and their extensions with an arbitrary mass at x = 0.) The general recursive formula is: with starting value f S (0) = P N (f C (0)), where P N (·) is the probability generating function of N . Probabilities are computed until their sum is arbitrarily close to 1.
The recursions are done in C to dramatically increase speed. One difficulty the programmer is facing is the unknown length of the output. This was solved using a common, simple and fast technique: first allocate an arbitrary amount of memory and double this amount each time the allocated space gets full. (17) and (18). This also requires a discrete severity distribution. However, there is no restriction on the shape of the frequency distribution. The package merely implements the sum (17), the convolutions being computed with R's function convolve, which in turn uses the Fast Fourier Transform. This approach is practical for small problems only, even on today's fast computers.

Exact calculation by numerical convolutions using
3. Normal approximation of the cdf, that is where µ S = E[S] and σ 2 S = VAR [S]. For most realistic models, this approximation is rather crude in the tails of the distribution.

Normal Power II approximation of the cdf, that is
S . The approximation is valid for x > µ S only and performs reasonably well when γ S < 1. See Daykin, Pentikäinen, and Pesonen (1994) for details. 5. Simulation of a random sample from S and approximation of F S (x) by the empirical cdf (8). The simulation itself is done with function simul (see Section 4). This function admits very general hierarchical models for both the frequency and the severity components.
Here also, adding other methods to aggregateDist is simple due to its modular conception.
The arguments of aggregateDist differ depending on the calculation method; see the help page for details. One interesting argument to note is x.scale to specify the monetary unit of the severity distribution. This way, one does not have to mentally do the conversion between the support of 0, 1, 2, . . . assumed by the recursive and convolution methods and the true support of S.
The function returns an object of class "aggregateDist" inheriting from the "function" class. Thus, one can use the object as a function to compute the value of F S (x) in any x.
For illustration purposes, consider the following model: the distribution of S is a compound Poisson with parameter λ = 10 and severity distribution Gamma(2, 1). To obtain an approximation of the cdf of S we first discretize the gamma distribution on (0, 22) with the unbiased method and a step of 0.5, and then use the recursive method in aggregateDist: R> fx <-discretize(pgamma(x, 2, 1), from = 0, to = 22, + step = 0.5, method = "unbiased", lev = levgamma(x, 2, 1)) R> Fs <-aggregateDist("recursive", model. The package defines a few summary methods to extract information from "aggregateDist" objects. First, there are methods of mean and quantile to easily compute the mean and obtain the quantiles of the approximate distribution: Second, the package introduces the generic functions VaR and CTE with methods for objects of class "aggregateDist". The former computes the value-at-risk VaR α such that where α is the confidence level. Thus, the value-at-risk is nothing else than a quantile. As for the method of CTE, it computes the conditional tail expectation Here are examples using object Fs obtained above: To conclude on the subject, Figure 6 shows the cdf of S using five of the many combinations of discretization and calculation method supported by actuar.

The continuous time ruin model
We now turn to the multi-period ruin problem. Let U (t) denote the surplus of an insurance company at time t, c(t) denote premiums collected through time t, and S(t) denote aggregate claims paid through time t. If u is the initial surplus at time t = 0, then a mathematically convenient definition of U (t) is As mentioned previously, technical ruin of the insurance company occurs when the surplus becomes negative. Therefore, the definition of the infinite time probability of ruin is We define some other quantities needed in the sequel. Let N (t) denote the number of claims up to time t ≥ 0 and C j denote the amount of claim j. Then the definition of S(t) is analogous to (16): assuming N (0) = 0 and S(t) = 0 as long as N (t) = 0. Furthermore, let T j denote the time when claim j occurs, such that T 1 < T 2 < T 3 < . . . Then the random variable of the interarrival (or wait) time between claim j − 1 and claim j is defined as W 1 = T 1 and For the rest of this discussion, we make the following assumptions: 1. premiums are collected at a constant rate c, hence c(t) = ct; 2. the sequence {T j } j≥1 forms an ordinary renewal process, with the consequence that random variables W 1 , W 2 , . . . are independent and identically distributed; 3. claim amounts C 1 , C 2 , . . . are independent and identically distributed.

Adjustment coefficient
The quantity known as the adjustment coefficient ρ hardly has any physical interpretation, but it is useful as an approximation to the probability of ruin since we have the inequality The adjustment coefficient is defined as the smallest strictly positive solution (if it exists) of the Lundberg equation where the premium rate c satisfies the positive safety loading constraint E[C − cW ] < 0. If C and W are independent, as in the most common models, then the equation can be rewritten as Function adjCoef of actuar computes the adjustment coefficient ρ from the following arguments: either the two moment generating functions M C (t) and M W (t) (thereby assuming independence) or else function h(t); the premium rate c; the upper bound of the support of M C (t) or any other upper bound for ρ.
For example, if W and C are independent, W ∼ Exponential(2), C ∼ Exponential(1) and the premium rate is c = 2.4 (for a safety loading of 20% using the expected value premium principle), then the adjustment coefficient is R> adjCoef(mgf.claim = mgfexp(x), mgf.wait = mgfexp(x, 2), + premium.rate = 2.4, upper = 1) [1] 0.1667 The function also supports models with proportional or excess-of-loss reinsurance (Centeno 2002). Under the first type of treaty, an insurer pays a proportion α of every loss and the rest is paid by the reinsurer. Then, for fixed α the adjustment coefficient is the solution of Under an excess-of-loss treaty, the primary insurer pays each claim up to a limit L. Again, for fixed L, the adjustment coefficient is the solution of For models with reinsurance, adjCoef returns an object of class "adjCoef" inheriting from the "function" class. One can then use the object to compute the adjustment coefficient for any retention rate α or retention limit L. The package also defines a method of plot for these objects.

Probability of ruin
In this subsection, we always assume that interarrival times and claim amounts are independent.
The main difficulty with the calculation of the infinite time probability of ruin lies in the lack of explicit formulas except for the most simple models. If interarrival times are Exponential(λ) distributed (Poisson claim number process) and claim amounts are Exponential(β) distributed, then If the frequency assumption of this model is defensible, the severity assumption can hardly be used beyond illustration purposes. Fortunately, phase-type distributions have come to the rescue since the early 1990s. Asmussen and Rolski (1991) first show that in the classical Cramér-Lundberg model where interarrival times are Exponential(λ) distributed, if claim amounts are Phase-type(π, T ) distributed, then ψ(u) = 1 − F (u), where F is Phase-type(π + , Q) with with t = −T e as in Section 2.1.
In the more general Sparre Andersen model where interarrival times can have any Phasetype(ν, S) distribution, Asmussen and Rolski (1991) also show that using the same claim severity assumption as above, one still has ψ(u) = 1 − F (u) where F is Phase-type(π + , Q), but with parameters and Q solution of In the above, s = −Se, I n is the n × n identity matrix, ⊗ denotes the usual Kronecker product between two matrices and ⊕ is the Kronecker sum defined as actuar: An R Package for Actuarial Science Function ruin of actuar returns a function object of class "ruin" to compute the probability of ruin for any initial surplus u. In all cases except the exponential/exponential model where (35) is used, the output object calls function pphtype to compute the ruin probabilities.
Some thought went into the interface of ruin. Obviously, all models can be specified using phase-type distributions, but the authors wanted users to have easy access to the most common models involving exponential and Erlang distributions. Hence, one first states the claim amount and interarrival times models with any combination of "exponential", "Erlang" and "phase-type". Then, one passes the parameters of each model using lists with components named after the corresponding parameters of dexp, dgamma and dphtype. If a component "weights" is found in a list, the model is a mixture of exponential or Erlang (mixtures of phase-type are not supported). Every component of the parameter lists is recycled as needed.
The following examples should make the matter clearer.

Simulation of compound hierarchical models
Function simul simulates portfolios of data following compound models of the form (16) where both the frequency and the severity components can have a hierarchical structure. The main characteristic of hierarchical models is to have the probability law at some level in the classification structure be conditional on the outcome in previous levels. For example, consider the following compound hierarchical model: for i = 1, . . . , I, j = 1, . . . , J i , t = 1, . . . , n ij and with The random variables Φ i , Λ ij , Ψ i and Θ ij are generally seen as risk parameters in the actuarial literature. The w ijt s are known weights, or volumes. Using as convention to number the data level 0, the above is a two-level hierarchical model. Goulet and Pouliot (2008) describe in detail the model specification method used in simul.
For the sake of completeness, we briefly outline this method here.
A hierarchical model is completely specified by the number of nodes at each level (I, J 1 , . . . , J I and n 11 , . . . , n IJ , above) and by the probability laws at each level. The number of nodes is passed to simul by means of a named list where each element is a vector of the number of nodes at a given level. Vectors are recycled when the number of nodes is the same throughout a level. Probability models are expressed in a semi-symbolic fashion using an object of mode "expression". Each element of the object must be named -with names matching those of the number of nodes list -and should be a complete call to an existing random number generation function, with the number of variates omitted. Hierarchical models are achieved by replacing one or more parameters of a distribution at a given level by any combination of the names of the levels above. If no mixing is to take place at a level, the model for this level can be NULL.
Function simul also supports usage of weights in models. These usually modify the frequency parameters to take into account the"size"of an entity. Weights are used in simulation wherever the name weights appears in a model.
Hence, function simul has four main arguments: 1) nodes for the number of nodes list; 2) model.freq for the frequency model; 3) model.sev for the severity model; 4) weights for the vector of weights in lexicographic order, that is all weights of entity 1, then all weights of entity 2, and so on.

R> pf
Portfolio of claim amounts The package defines methods for four generic functions to easily access key quantities of the simulated portfolio.
1. By default, the method of aggregate returns the values of aggregate claim amounts S ijt in a regular matrix (subscripts i and j in the rows, subscript t in the columns). The method has a by argument to get statistics for other groupings and a FUN argument to get statistics other than the sum: 3. The method of severity (a generic function introduced by the package) returns the individual claim amounts C ijtu in a matrix similar to those above, but with a number of columns equal to the maximum number of observations per entity, Thus, the original period of observation (subscript t) and the identifier of the severity within the period (subscript u) are lost and each variate now constitutes a "period" of observation. For this reason, the method provides an argument splitcol in case one would like to extract separately the individual claim amounts of one or more periods: In addition, all methods have a classification and a prefix argument. When the first is FALSE, the classification index columns are omitted from the result. The second argument overrides the default column name prefix; see the simul.summaries help page for details.
Function simul was used to simulate the data in Forgues, Goulet, and Lu (2006).

Credibility theory
Credibility models are actuarial tools to distribute premiums fairly among a heterogeneous group of policyholders (henceforth called entities). More generally, they can be seen as prediction methods applicable in any setting where repeated measures are made for subjects with different risk levels.
The credibility theory facilities of actuar consist of the matrix hachemeister containing the famous data set of Hachemeister (1975) and the function cm to fit credibility models.

Hachemeister data set
The data set of Hachemeister (1975) consists of private passenger bodily injury insurance average claim amounts, and the corresponding number of claims, for five U.S. states over 12 quarters between July 1970 and June 1973. The data set is included in the package in the form of a matrix with 5 rows and 25 columns. The first column contains a state index (state), columns 2-13 contain the claim averages (ratio.1, . . . , ratio.12) and columns 14-25 contain the claim numbers (weight.1, . . . , weight.12).

Hierarchical credibility model
The linear model fitting function of R is named lm. Since credibility models are very close in many respects to linear models, and since the credibility model fitting function of actuar borrows much of its interface from lm, we named the credibility function cm.
Function cm acts as a unified interface for all credibility models supported by the package. Currently, these are the unidimensional models of Bühlmann (1969) and Bühlmann and Straub (1970), the hierarchical model of Jewell (1975) (of which the first two are special cases) and the regression model of Hachemeister (1975). The modular design of cm makes it easy to add new models if desired.
This subsection concentrates on usage of cm for hierarchical models.
There are some variations in the formulas of the hierarchical model in the literature. We compute the credibility premiums as given in Bühlmann and Jewell (1987) or Bühlmann and Gisler (2005). We support three types of estimators of the between variance structure parameters: the unbiased estimators of Bühlmann and Gisler (2005) (the default), the slightly different version of Ohlsson (2005) and the iterative pseudo-estimators as found in Goovaerts and Hoogstad (1987) or Goulet (1998). See Belhadj, Goulet, and Ouellet (2008) for further discussion on this topic.
The credibility modeling function assumes that data is available in the format most practical applications would use, namely a rectangular array (matrix or data frame) with entity observations in the rows and with one or more classification index columns (numeric or character). One will recognize the output format of simul and its summary methods.
Then, function cm works much the same as lm. It takes in argument: a formula of the form~terms describing the hierarchical interactions in a data set; the data set containing the variables referenced in the formula; the names of the columns where the ratios and the weights are to be found in the data set. The latter should contain at least two nodes in each level and more than one period of experience for at least one entity. Missing values are represented by NAs. There can be entities with no experience (complete lines of NAs).
In order to give an easily reproducible example, we group states 1 and 3 of the Hachemeister data set into one cohort and states 2, 4 and 5 into another. This shows that data does not have to be sorted by level. The fitted model using the iterative estimators is: R> X <-cbind(cohort = c(1, 2, 1, 2, 2), hachemeister) R> fit <-cm(~cohort + cohort:state, data = X, ratios = ratio.1:ratio.12, + weights = weight.1:weight.12, method = "iterative") R> fit Call: cm(formula =~cohort + cohort:state, data = X, ratios = ratio.1:ratio.12, weights = weight.1:weight.12, method = "iterative") The function returns a fitted model object of class "cm" containing the estimators of the structure parameters. To compute the credibility premiums, one calls a method of predict for this class: The results above differ from those of Goovaerts and Hoogstad (1987) for the same example because the formulas for the credibility premiums are different.

Bühlmann and Bühlmann-Straub models
As mentioned above, the Bühlmann and Bühlmann-Straub models are simply one-level hierarchical models. In this case, the Bühlmann-Gisler and Ohlsson estimators of the between variance parameters are both identical to the usual Bühlmann and Straub (1970) estimator, and the iterative estimator is better known as the Bichsel-Straub estimator.

Regression model of Hachemeister
The regression model of Hachemeister (1975) is a generalization of the Bühlmann-Straub model. If data shows a systematic trend, the latter model will typically under-or overestimate the true premium of an entity. The idea of Hachemeister was to fit to the data a regression model where the parameters are a credibility weighted average of an entity's regression parameters and the group's parameters.
In order to use cm to fit a credibility regression model to a data set, one has to specify a vector or matrix of regressors by means of argument xreg. For example, fitting the model X it = β 0 + β 1 (12 − t) + ε t , t = 1, . . . , 12 to the original data set of Hachemeister is done with R> fit <-cm(~state, hachemeister, xreg = 12:1, ratios = ratio. Computing the credibility premiums requires to give the "future" values of the regressors as in predict.lm, although with a simplified syntax for the one regressor case: R> predict(fit, newdata = 0) [1] 2437 1651 2073 1507 1759

Documentation
In addition to the help pages required by the R packaging system, the package includes vignettes and demonstration scripts; running R> vignette(package = "actuar") and R> demo(package = "actuar") at the R prompt will give the list of each.

Conclusion
The paper presented the facilities of the R package actuar in the fields of loss distribution modeling, risk theory, simulation of compound hierarchical models and credibility theory. We feel this version of the package covers most of the basic needs in these areas. In the future we plan to improve the functions currently available and to start adding more advanced features. For example, future versions of the package should include support for dependence models in risk theory and better handling of regression credibility models.
Obviously, the package left many other fields of Actuarial Science untouched. For this situation to change, we hope that experts in their field will join their efforts to ours and contribute code to the actuar project. The project will continue to grow and to improve by and for the community of developers and users.
Finally, if you use R or actuar for actuarial analysis, please cite the software in publications. Use