Feller-Pareto and Related Distributions: Numerical Implementation and Actuarial Applications

Actuaries model insurance claim amounts using heavy tailed probability distributions. They routinely need to evaluate quantities related to these distributions such as quantiles in the far right tail, moments or limited moments. Furthermore, actuaries often resort to simulation to solve otherwise untractable risk evaluation problems. The paper discusses our implementation of support functions for the Feller-Pareto distribution for the R package actuar . The Feller-Pareto defines a large family of heavy tailed distributions encompassing the transformed beta family and many variants of the Pareto distribution.


Introduction
Actuaries are experts in evaluating the likelihood and financial consequences of future events. A pivotal part of their work is the modeling of the size and frequency of insurance claims. With probability models for the claims process in hand, actuaries can compute insurance premiums, determine the amount a company has to set aside in its reserve to cope with future events, evaluate the risk the company will not be able to meet its obligations, or run simulations to compare business strategies or solve otherwise untractable problems.
Needless to say, the probability distribution is the cornerstone of actuarial and statistical modeling. In this data science golden age, one needs to be able to compute or simulate various quantities for a large array of probability distributions. Actuarial applications, for example, require accurate computation of probabilities or quantiles in the far right tail of distributions; easy computation of raw and limited moments; fast and reliable simulation of random variates (see, e.g., McNeil, Frey, and Embrechts 2015). The R statistical system (R Core Team 2022a) provides in its base distribution a number of functions to compute the probability density function (PDF) or the probability mass function (PMF), the cumulative distribution function (CDF) and the quantile function, as well as functions to generate variates for the most common discrete and continuous distributions of statistics.
The number of probability distributions implemented in base R has remained remarkably stable over time. The R Core Team has elected to delegate support for additional distributions to extension packages developed and maintained by the R community. The Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org is the central repository of R packages. In their extensive CRAN Task View on probability distributions, Dutang and Kiener (2022) list packages that support probability distributions not found in base R (263 packages as of March 07, 2022).
First released on CRAN in 2005, actuar (Dutang, Goulet, and Pigeon 2008), available at https://CRAN.R-project.org/package=actuar, provides specialized functions for actuarial applications, but also support functions for a large number of continuous size distributions useful for loss severity modeling; for phase-type distributions used in computation of ruin probabilities; for zero-truncated and zero-modified extensions of the discrete distributions commonly used in loss frequency modeling; for the heavy tailed Poisson-inverse Gaussian discrete distribution. The package also introduces support functions to compute raw moments, limited moments and the moment generating function (when it exists) of continuous distributions. These functions all prove useful in loss modeling and risk evaluation.
We introduce with this paper our latest addition to the range of heavy tailed probability models supported by actuar: the Feller-Pareto distribution and related Pareto distributions with a location parameter. To the best our knowledge, the Feller-Pareto distribution has not yet been extensively used for insurance loss modeling. However, Brazauskas (2002) proposes its use to model the size of insurance claims and Cummins, Dionne, McDonald, and Pritchett (1990) relies on a special case to estimate the tail distribution and reinsurance premiums for fire claims.
The paper is organized as follows. For the reader not familiar with the Feller-Pareto family of distributions, we first summarize how it extends the transformed beta family in Section 2. Next, we discuss the theoretical and numerical challenges faced during the numerical implementation in Section 3. In particular in Section 3.2, we also compare various algorithms to generate random variates from the Feller-Pareto and related distributions. In Section 4, we provide a small case study based on real world data. The workhorses behind the R support functions are implemented in C for unrivaled speed and efficiency. We close the paper by briefly explaining how R package developers can now take advantage of these routines through actuar's new package API in Section 5.
For the purpose of this paper, a heavy tailed distribution is a positive distribution without a finite moment generating function. Other definitions exist, such as the domain of attraction in extreme value theory (Embrechts, Klüppelberg, and Mikosch 1997).

Feller-Pareto distribution
The Pareto distribution is a well known probability distribution originally used by economists to model income in a society. Due to its heavy tail, the distribution also plays an important role in actuarial science to model loss cost, notably in reinsurance. In recent years, many authors proposed increasingly flexible variants of the Pareto that eventually led to the distri-bution now generally known as the Feller-Pareto; see Arnold (2015); Kleiber and Kotz (2003) for exhaustive surveys. Arnold (2015) proposes the Feller-Pareto distribution as the distribution of the random variable
One will note the presence of the location parameter µ. In actuarial applications, this proves useful for the modeling of claim amounts above a certain threshold such as, for example, for reinsurance contracts stipulating that the reinsurance company is responsible for the total amount of losses above a certain limit.
It is well known that the distribution of the ratio U/V in (1) is a beta distribution of the second kind (see Johnson, Kotz, and Balakrishnan 1994a). Therefore, the CDF of X is given in terms of the (incomplete) beta function by where I(a, b; x) is the (regularized) incomplete beta function -or incomplete beta ratiorepresenting the CDF of the beta distribution (henceforth denoted B(a, b)). We define as the incomplete beta function, and as the beta function (Abramowitz and Stegun 1972).
The PDF corresponding to the CDF F (x) above is We can rewrite this equation as follows: which highlights a link between the Feller-Pareto distribution and the beta (of the first kind). Indeed, it is also well known that if U ∼ G(τ, 1), V ∼ G(α, 1) and U and V are independent, then B = V /(V + U ) ∼ B(α, τ ). We can rewrite (1) as thereby defining the Feller-Pareto as a transformation of a beta distribution. It is actually under this form that Feller (1968) and Johnson, Kotz, and Balakrishnan (1994b); Johnson et al. (1994a) introduce the Feller-Pareto.
One final remark that will prove useful in the sequel. Another well known result about beta distributions is that if B has a B(α, τ ) distribution, then the reflected variableB = 1 − B has a B(τ, α) distribution. Therefore, by a simple manipulation of (4), we can equivalently define the Feller-Pareto as the distribution of the random variable whereB ∼ B(τ, α).

Moments
The computation of moments plays a central role in actuarial applications. The first raw moment -or expected value -is tightly linked to the concept of insurance premium, whereas the second central moment -or variance -measures the risk associated with this premium.
A general formula for the (raw) moment of order k is easy to derive for the Feller-Pareto distribution with µ = 0 using definition (1) as a ratio of independent gamma distributions. Indeed, in that case, As Arnold (2015) and Klugman, Panjer, and Willmot (2019) show, we readily obtain For the more general case µ ̸ = 0, one derives the moment of order k using the binomial expansion The resulting general equation for the raw moment of order k of the Feller-Pareto is: Note, however, that due to constraints on the validity of the binomial expansion, the above result is only valid for nonnegative integral values of k < αγ (Arnold 2015).

Limited moments
Limited moments play an important role in actuarial science, especially the limited expected value as the (pure) premium for an insurance contract subject to a maximal insurable loss as is common in homeowners insurance (Klugman et al. 2019, Chapter 3).
Let a ∧ b = min(a, b) denote the minimum of a and b. We define the limited moment of order k of a random variable X as The limited expected value is the special case with k = 1.
Like raw moments, the limited moments of the Feller-Pareto are simpler to derive in the case µ = 0. Klugman et al. (2019, p. 463) already provide for x ≥ 0, k > −τ γ, and with, as in Section 2.1, u = y/(1 + y) and y = ((x − µ)/θ) γ . (We will further discuss in Section 3.4 how Γ(a)Γ(b)I(a, b; x) may also be valid for b < 0.) Using again binomial expansion, we can similarly generalize the above formula for the limited moment of order k when µ ̸ = 0: for x ≥ µ and, this time, nonnegative integral values of k.

Generation of random variates
The definitions of Section 2.1 provide us with three straightforward algorithms to generate random variates from the Feller-Pareto distribution.
The first algorithm relies on definition (1) of the Feller-Pareto as a ratio of Gamma random variables.
The second and third algorithms are very similar in that they both use transformations of a beta of the first kind, one using definition (4), the other using definition (5). We introduce the third algorithm here since it was actually the one used in actuar -with µ = 0 -to simulate variates from the transformed beta distribution.

Historical remarks
The later part of the twentieth century saw a lot of simultaneous independent work on probability distributions to model income and loss cost. On the one side, economists proposed distributions to model income -such as the Singh-Maddala or Dagum distributions -until McDonald (1984) introduced a four-parameter generalized beta of the second kind distribution with PDF This is exactly a Feller-Pareto with µ = 0, b = θ, b = γ, p = τ and q = α. See Atkinson and Bourguignon (2015) for a thorough discussion of income distributions and Cummins et al. (1990, Table 2) for the interpretation as a mixture.
On the other side, as reported by Kleiber and Kotz (2003), actuaries also proposed various distributions to model insurance loss cost -such as the Fisk or Lomax distributionsuntil Venter (1983) proposed a distribution with a PDF identical to (10) and called it the transformed beta distribution.

Special cases
One will find many variants of the Pareto distribution in the literature, most of which can be considered special cases of the Feller-Pareto. We review some of them below. For the reader interested to know more, we highly recommend the extensive work of Arnold (2015) on the Pareto family of distributions.
As a first case of the Feller-Pareto distribution, let us consider the aptly named Pareto I, or Single Parameter Pareto (Klugman et al. 2019). We obtain this distribution by setting the location parameter µ equal to the scale parameter θ, and by setting the shape parameters γ and τ equal to 1. The resulting CDF is Even though this distribution appears to have two parameters, its name stems from the fact that only α is an actual parameter for estimation purposes, θ being a known threshold.
Next, the Pareto II distribution is obtained by setting γ = τ = 1 in the Feller-Pareto. The CDF is then When µ = 0, we obtain what is generally simply called the Pareto distribution.
Then we have the Pareto III distribution, for which we can find many different definitions in the literature. The definition of Arnold (2015) is in line with the Pareto II, but with a different shape parameter. The CDF is which is a Feller-Pareto with shape parameters α = τ = 1. Johnson et al. (1994b,a), on the other hand, use the original definition of Pareto himself However, as explained in Kleiber and Kotz (2003), parameter estimation is hard with this definition of the Pareto III. For this reason, we will keep the above Arnold (2015) definition.
Finally, the Pareto IV is a Feller-Pareto where only the shape parameter τ is set equal to 1. The CDF is The Feller-Pareto generalizes other common continuous extreme value distributions found in Appendix A of Klugman et al. (2019): the case µ = 0 is the transformed beta distribution, as mentioned above; the case µ = 0, τ = 1 is the Burr distribution type 12; the case µ = 0, τ = 1, α = 1 is the loglogistic distribution; the case µ = 0, α = 1 is the inverse Burr or Burr type 3 distribution. Figure 1 shows the relationships between the Feller-Pareto and the distributions of the transformed beta family of Klugman et al. (2019). One may connect the diagram of Figure 1 to the impressive diagram of relationships between univariate distributions of Leemis and McQueston (2008), starting with the Pareto or generalized Pareto distributions. (Cummins et al. 1990) also propose a diagram for the relationships between the members of generalized beta of the second kind family.

R implementation
As already mentioned in the introduction of this paper, R includes functions to compute the PDF, the CDF and the quantile function of standard probability laws, as well as functions to generate variates from these laws. For some root foo (e.g., pois for the Poisson distribution or gamma for the gamma distribution), the support functions are named dfoo, pfoo, qfoo and rfoo, respectively.
Since version 0.9-4, released in 2007, 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 et al. (2019) and not already present in base R, excluding the log-t but including the Loggamma distribution of Hogg and Klugman (1984). Version 2.0-0 of the package also added support for the Inverse Gaussian distribution. The package vignette "distributions" keeps the up-to-date list of supported distributions and root names of the R functions.
In addition to the d, p, q and r functions, actuar provides m, lev and mgf functions to compute, respectively, theoretical raw moments m k , theoretical limited moments E (X ∧ x) k and the moment generating function Parameter Argument when the latter exists. The moment generating function of the Feller-Pareto distribution does not exist since the density (3) does not decrease at an exponential rate.
With version 3.0-0, released alongside this paper, actuar adds support for the Feller-Pareto, Pareto IV, Pareto III and Pareto II distributions. The Pareto I and the standard Pareto distribution without a location parameter were already supported by actuar. The root of the support functions are, respectively, fpareto, pareto4, pareto3 and pareto2. Therefore, for the Feller-Pareto, the package now provides the functions dfpareto, pfpareto, qfpareto, rfpareto, mfpareto and levfpareto. Table 1 shows the correspondence between the parameters of the Feller-Pareto as defined in Section 2.1 and the arguments of the R functions.
In the sequel, we only discuss the implementation of the Feller-Pareto functions since the Pareto IV, III and II functions are merely special cases.

Computational challenges
This section discusses numerical and performance issues we encountered while implementing in R some of the functions of the fpareto family.

Cumulative distribution function
We first look at the implementation of pfpareto to compute the CDF (2) of the Feller-Pareto distributions. This requires to evaluate the incomplete beta ratio I (a, b; x). Since this is just the CDF of the beta distribution, we have at our disposal the function pbeta in base R to carry the computation efficiently.
Actually, it proves not so simple. In our accuracy tests, we observed that when y = [(x−µ)/θ] γ gets "large", then u = y/(1 + y) becomes numerically equal to 1 and, accordingly, pbeta returns 1. This can represent a significant jump in the CDF of the Feller-Pareto; see the dashed line in Figure 2. In actuarial applications, events in the right tail of a heavy tailed probability distribution are the most critical ones for an insurance company: although very unlikely, they represent a huge financial exposure.
We needed a fix for the computation of the probabilities -or, equivalently, of the quantilesin the far right tail of the distribution. It came from the symmetry property of the incomplete beta ratio. As stated in Abramowitz and Stegun (1972, Section 6.6.3): Therefore, we can also write the CDF (2) of the Feller-Pareto distribution as Now, the value 1 − u = 1/(1 + y) underflows to zero much later than y/(1 + y) "overflows" to one for large values of y. Using the approach in (16) coupled with the standard argument lower.tail = FALSE of pbeta to compute 1 − I(a, b; x), we are able to compute probabilities accurately much farther in the right tail of the distribution; see the solid line in Figure 2.
The previous discussion also holds for very small values of y -although to a lesser extend since the set of representable floating point numbers is much more dense near zero (Goldberg 1991). Since any accuracy improvement in the right tail of the distribution should not be canceled by a deterioration in the left tail (near µ), our implementation of pfpareto actually computes the CDF with (2) when u = y/(1 + y) ≤ 0.5, and with (16) otherwise.

Generation of random variates
The other main issue we faced in our R implementation was with the rfpareto function to generate random variates from the Feller-Pareto distribution. Section 2.4 outlined the three simulation algorithms at our disposal.
When comparing simulation algorithms, we are mostly interested in their relative speed and ability to reproduce the underlying distribution, especially for limit cases. We use functions rgamma and rbeta of base R to generate gamma and beta variates, respectively. Therefore, we leave aside considerations of randomness and efficiency since we are building on well
established algorithms and routines to generate variates from these distributions (Cheng 1978;Ahrens and Dieter 1982).
One will first note that Algorithm 1 requires two gamma variates to generate a single Feller-Pareto variate, whereas Algorithms 2 and 3 require only one beta variate. This is an advantage for the beta-based algorithms, but not in a 2 : 1 ratio. Indeed, rgamma takes less time to generate a variate than rbeta. As Table 2 shows, generating two gamma variates is actually only about 30% slower that generating a single beta variate. These performance measurements carry over to the generation of Feller-Pareto variates, as shown in Table 3.
(We obtained benchmarking times with the very convenient function benchmark of package rbenchmark (Kusnierczyk 2012).) The other issue we investigated was performance of the algorithms in the main limit cases. First, when α → 0 and τ → ∞, variates from a B(τ, α) tend to 1. Due to the negative power in its second step, Algorithm 3 tends, in such a situation, to yield indeterminate values -"not a number", or NaN in R terminology. Algorithm 2 does not exhibit this behavior. Figure 3 (left graphic) shows the effect on the empirical CDF of samples generated with Algorithms 2 and 3 for a small value of α and a large value of τ . Algorithm 2 is clearly preferable.
The indeterminate value issue of Algorithm 3 also arises in Algorithm 1 when α → τ → 0 and  Table 3: Relative time required to simulate 10 6 Feller-Pareto variates with three different algorithms (on average over the number of replications). Location parameter µ = 0; scale parameter θ = 1; shape parameters α and γ uniformly distributed on (0, 100); shape parameter τ uniformly distributed on (100, 1000). both gamma variates may be numerically equal to 0. Figure 3 (right graphic) compares the theoretical CDF with typical empirical CDF of samples generated with Algorithms 1 and 2 (notice on the y-axis how this graphic is zoomed in). In our various tests, the empirical CDF with Algorithm 1 was almost consistently farther from the true CDF than the one with Algorithm 2.
As a result of the above observations, we recommend Algorithm 2 to generate Feller-Pareto random variates. This is the the algorithm we implemented in actuar. Accordingly, we also changed the algorithm used in the simulation of transformed beta variates from Algorithm 3 to Algorithm 2.

Moments
Numerical computation of raw moments does not pose significant difficulties. It is worth noting, though, that the function mfpareto computes the ratio of gamma functions in (6) and (7) as a ratio of beta functions since This way, we take advantage of the optimizations and careful handling of limit cases of routine beta of the C API of R (R Core Team 2022b).
Furthermore, in order to reduce the number of operations, we actually compute the moment of order k in the case µ ̸ = 0 with the following expression that is algebraically equivalent to (7):

Limited moments
Before we turn to the computation of limited moments of the Feller-Pareto family of distributions, we must first study an alternative definition of the (regularized) incomplete beta function that is valid for negative values of the parameters. Klugman et al. (2019) introduce the alternative definition to extend the range of admissible values for limited expected value functions.
As seen in Abramowitz and Stegun (1972, Section 6.6), we have the following relation for the integral on the right hand side of (17): is the Gauss hypergeometric series. With the above definition, the incomplete beta function also admits negative, non integral values for parameters a and b. Of more interest here is the case where b < 0, b ̸ = −1, −2, . . . and a > 1 + ⌊−b⌋. Integration by parts of (18) yields where r = ⌊−b⌋. Following the terminology of actuar, we will call (18) the beta integral. (It is also known in the literature as the normalized incomplete beta function.) The package includes a C routine betaint -and a non exported R interface of the same name -that is a rather straightforward C implementation of (19). Its sole usage in the package is to evaluate the limited moments of distributions of the Feller-Pareto family.
Given the above, we rewrite the general expression (9) of the limited moment of order k as or, equivalently, (τ, α; u)).
Note that for the expression above to be valid, α − j/γ, j = 0, 1, . . . , k must not be a negative integer. This is the expression that we implemented in the C routine levfpareto.
For the case µ = 0 -distributions of the transformed beta family -we use (8) rewritten as still with α − k/γ not a negative integer.
It is worth noting that the C routines behind limited moments functions carry a lot of computations in the log scale to delay numerical overflow and, therefore, to gain accuracy in the tails. For example, levfpareto and related routines compute u = y 1 + y = 1 1 + y −1 , y = as u = e log u , with The routines make use of function log1pexp of the C API of R to compute log(1 + e x ) accurately, notably for large x.
Upon analyzing the accuracy of the levtrbeta function in a previous version of actuar, we discovered discrepancies in the tail of the limited expected value function which in practice may lead to erroneous premium computation; see Section 4.2. Figure 4 shows the problem whether the beta integral (18) is called with a positive (left graphic) or negative (right graphic) value of its second parameter b. The cause is exactly the same as in Section 3.1: u gets numerically equal to 1 when y gets large. The cure is also the same: compute any instance of I(a, b; u) as 1 − I(b, a; 1 − u) using 1 − u = (1 + y) −1 when u > 0.5. This also led us to add an argument to the betaint routine in order to pass the numerically accurate value of 1 − x to the incomplete beta function in (19) and, simultaneously, to all instances of this term in the sum.
The limited moments functions for all distributions of the Feller-Pareto family and the beta integral function incorporate the above improvements in version 3.0-0 of actuar. Figure 5 shows that their behavior in the tail is now consistent with the theoretical expressions.

Case study
We illustrate in this section how the Feller-Pareto distribution -and the functions we provide in actuar -can prove useful to model insurance loss costs, evaluate premiums and simulate data.
To this end, we use a dataset well known in actuarial circles as the Danish dataset. Collected by Copenhagen Reinsurance, it consists of aggregate claim amounts for 2,492 fire losses over . The dataset has been studied by, e.g., Embrechts et al. (1997);McNeil (1997); Resnick (1997); Davison (2003). The R package SMPracticals (Davison 2019) provides the Danish dataset in the time series object danish. Since the dates of the events are irrelevant in our analysis, we use the data as a simple atomic vector.
R> data("danish", package = "SMPracticals") R> danish.loss <-as.numeric(danish) As Figure 6 shows, the Danish loss amounts exhibit a heavy, Pareto type, tail. McNeil (1997); Resnick (1997) propose alternative depictions of this behavior using a QQ-plot against the exponential distribution, a mean excess plot and a Hill plot. Summary statistics and right tail quantiles are also quite revealing.  McNeil 1997) fitted three distributions on the loss amounts: a Pareto I, a shifted lognormal and a generalized Pareto, the latter being a Pareto II distribution with a different parametrization widely used in extreme value theory (see, e.g., Embrechts et al. 1997). None of these distributions catch the steep curvature of the cumulative distribution function, although (Mc-Neil 1997) shows the generalized Pareto successfully catches the tail behavior for losses over DKK 20,000,000.

Loss modeling
The first step in our case study will be to fit the Feller-Pareto distribution to the data using maximum likelihood estimation (MLE). Our main R tool will be function fitdist from package fitdistrplus (Delignette-Muller and Dutang 2015), which wraps and extends the well known function fitdistr from package MASS (Venables and Ripley 2002).
Since the Danish dataset contains positive loss amounts, we have a natural value of µ = 0 for the location parameter of the Feller-Pareto. We could try to estimate this parameter along with the shape parameters and the scale parameter, but this proves numerically cumbersome. Brazauskas (2002) provides explicit formulas for the Fisher information matrix when the location parameter is set to 0. However, for the current example, we simply rely on the numerical values provided by the R optimization routines.
We use the bounded BFGS optimization algorithm to restrict the parameters to positive values (Bonnans, Gilbert, Lemaréchal, and Sagastizábal 2006). For numerical reasons, we actually restrict the shape and scale parameters to [ϵ, ∞), where ϵ = 2 −52 is the machine epsilon in double precision.

Fitting of the distribution ' fpareto ' by maximum likelihood
In addition to the analysis of the goodness-of-fit plots, we can also check statistical measures or criteria (Klugman et al. 2019, Chapter 15)  Our intent here is merely to illustrate how the numerical tools we provide may prove useful to fit Feller-Pareto distributions to loss data. Still, the above statistics show that the model fares very well against the Weibull-Inverse Weibull mixture determined as best among 256 composite models by Grün and Miljkovic (2019) for the same dataset.

Expected losses
Our Feller-Pareto model has αγ = 1.297. By (6), its expected value exists. We may therefore easily compute the expected losses (in million DKK) for the Danish dataset using function mfpareto.  After McNeil (1997), we also compute the expected losses for an excess-of-loss reinsurance contract where the reinsurer is responsible for losses between d and d + c, where d is the deductible value and c, the cover. Let X be the original loss amount and X re , the reinsurance cost. We have X re = min(max(X − d, 0), c) = min(X, d + c) − min(X, d).
The expected loss cost for the reinsurer is E [X re ] = E [min(X, d + c)]−E [min(X, d)]. We used function levfpareto to compute these expected values and to draw Figure 8, which displays the expected cost E [X re ] as a function of the deductible d for different values of covers c. One can see the drastic effect of the deductible and the cover on expected costs.

Simulation of aggregate claim amounts
For pricing and risk evaluation purposes, actuaries usually need to take into account not only the severity of losses, but also their frequency. 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 X j represent the amount of claim j (or severity). Then, we have the random sum where we assume that X 1 , X 2 , . . . are mutually independent and identically distributed random variables each independent of N . When N has a Poisson distribution, we say that the distribution of S is a compound Poisson, see, e.g., Klugman et al. (2019).
But for a few simple cases, finding the distribution of the random variable S is a difficult task. However, one may turn to simulation to get good estimates of quantities related to this distribution.
Let us simulate a large sample from a compound Poisson distribution where the severity of losses is distributed according to our Feller-Pareto. This is straightforward using function rcomppois from package actuar. (See Goulet and Pouliot (2008) for an explanation of the syntax to specify simulation models.) R> x <-rcomppois(1e6, lambda = 10, rfpareto(min = 0, + shape1 = par.fp["shape1"], shape2 = par.fp["shape2"], + shape3 = par.fp["shape3"], scale = par.fp["scale"])) We may now, for example, get fairly good approximations of the upper tail quantiles of the distribution of S. Since the sample size of x is one million, we can reasonably assess quantiles up to probability level p = 0.999999 = 1 − 10 −6 based on the non-parametric estimator provided by base R function quantile. Otherwise, we will need either to have a higher number of simulations or to use the extreme value theory, e.g., McNeil (1997).

Accessing the C routines
The actual workhorses behind the R functions of actuar presented in this paper are C routines. Package developers who need one or many of these routines in their C code can take advantage of them through the new package API of actuar. The header file inst/include/actuarAPI.h in the package installation directory contains the declarations of the more than 200 routines part of the API.
We direct the interested reader to the "distributions" package vignette for further instructions to access the C routines through the API. The companion package expint (Goulet 2022) ships with a complete test package implementing these instructions. See the vignette of the latter package for more information.

Conclusion
This paper introduced the implementation of the Feller-Pareto distribution and related Pareto distributions with a location parameter in version 3.0-0 of the R package actuar. The Feller-Pareto defines a large family of distributions encompassing the transformed beta family and many variants of the Pareto distribution.
A large proportion of the numerical challenges discussed in the paper are related to the fact that, for most distributions of the Feller-Pareto family, we need to evaluate the incomplete beta function for a value that can jump to 1 if not carefully computed to avoid numerical overflows. By applying a simple correction to many functions present in previous versions of actuar, we were able to fix long standing discrepancies in CDF and limited moment functions, as well as to provide robust new functions for the Feller-Pareto distribution. We also improved the simulation algorithm for all distributions of the Feller-Pareto family.
A case study based on a famous heavy tailed insurance dataset showed how one can use the various support functions for the Feller-Pareto distribution in modeling, graphics or simulation.