Open access peer-reviewed chapter

What Determines EP Curve Shape?

Written By

Frank Xuyan Wang

Submitted: 12 September 2018 Reviewed: 01 December 2018 Published: 07 January 2019

DOI: 10.5772/intechopen.82832

From the Edited Volume

Applied Mathematics

Edited by Bruno Carpentieri

Chapter metrics overview

1,329 Chapter Downloads

View Full Metrics


Propose use kurtosis divided by skewness squared as shape factor, and use the global or conditional minimum/maximum of this shape factor for selecting and differentiating distribution families. Semi-empirical formulas for that lower/upper bound are calculated for various distribution families, with the aid of Computer Algebra System, for fitting hard to match distributions. Previous studies show high CV distribution is hard to fit and simulate, this study extends that conclusion to cases with low CV but still hard to match EP curves, characterized by having shape factors close to 1. The maximal likelihood approach of distribution fit can tell us which distribution family is better suited for an empirical distribution, but the shape factor range information can tell us why a distribution cannot fit well, or is not suitable at all. So the shape factor, in a sense, determines the EP curve shape.


  • Skewness
  • kurtosis
  • TVaR
  • shape factor
  • reinsurance
  • computer algebra system
  • Beta distribution
  • Kumaraswamy distribution
  • asymptotic expansion
  • GB2 distribution
  • numerical optimization
  • generalized hyperbolic distribution

1. Introduction

In reinsurance industry, losses for a contract are simulated and represented by the losses cumulative distribution function (CDF), survival, or quantile functions. The plots of these functions are called the EP curves with the following terminology [1]: for a given annual or aggregated loss, the probability of seeing annual loss exceeding that loss is the exceeding-probability (EP) or aggregate-exceeding-probability (AEP). The average of all annual losses exceeding that given loss is the AEP tail value at risk, called the AEP TVaR, or simply TVaR. The EP curve is represented by a table consisting of pairs of probability and loss. It is desirable to fit a parametric distribution to this table for a more succinct representation and more reasonable interpolations for values not in the table. Then which distribution family to use and what characteristics of the data are needed or determine the distribution are the questions to answer.

The (scaled) Beta distribution is widely used in reinsurance for fit loss or loss ratio, perhaps because the Beta distribution has only two parameters and very simple formulas for mean and standard deviation using these parameters, whose inverse function also has simple formulas, so that the two statistics of mean and standard deviation can be used to easily determine the parameters.

For about 85% of the perils, this approach works well, in the sense that the TVaR of the fitted distribution for quantile of interest, such as the 0.96, 0.99 or 0.996 quantile TVaR which is needed for pricing and risk monitoring, is close to a few percent of the original data TVaR. The remaining 15% perils, such as the North American Tornado Hail (NATH), Australia Wind Storm (AUWS), Hawaii Wind Storm (HIWS), and Mexico Earthquake (MXEQ), can have more than 10% deviations.

The maximum likelihood estimation method is a way to find alternative fitting distributions [2, 3]. Instead of finding approximations of the smoothed empirical distribution, we optimize an objective function whose optimum solution gives us the candidate distribution form. Suppose the annual losses x i occurred n i times in our observation; to find a probability function that gives probabilities p i for these losses, we just maximize the objective function i p i n i . It is easily seen that for the optimum solution we have p i p j = n i n j : the relative occurring frequency is maintained in the probability function. In the objective function, if we replace the p i by a power function of p i , the conclusion still holds, but not if we use a logarithm or exponential function.

While the maximum likelihood approach works well for many perils and identifies a few best fitted distribution families (Mathematica has more than 200 distribution families that can be used for extensive searches), it did not work for the NATH peril. The NATH has {Mean, StandardDeviation, Skewness, Kurtosis, 0.99TVaR} = {7418611.10904006, 9517336.93024634, 5.99378199789956, 65.8901734355745, 68867612.8345741}.

This is not contradictory to the maximum likelihood principal, since in any implementation, only known forms of the probability density function (PDF) and as-small-as-possible numbers of parameters can be used. To overcome this limitation, we need to look into the particularity of those distributions and come up with or select more suitable function forms for the PDF or CDF. In [4] it is found that a high coefficient of variation (CV) distribution is hard to fit or simulate. But the NATH has a small CV of 1.28. The skewness and kurtosis alone also not differentiate them from other distributions.

Trial and error found the empirical rule that these hard distributions have small values of kurtosis divided by skewness squared, Table 1. This finding prompted us for the study of the property of kurtosis/skewness^2 (K/S^2), henceforth will be called the shape factor (SF).

Peril CV Shape factor
NATH 1.283 1.834
AUWS 5.711 1.260
HIWS 4.678 1.238
MXEQ 3.930 1.878

Table 1.

Numerical characteristics of a few hard to fit and simulate perils.

Numerical optimization or solution will be our primary tool for this SF study. Analytical deduction, symbolic algebra, and symbolic limit from computer algebra system (CAS) Mathematica will be another major tool, as well as Mathematica’s plot functions. Those plots can help reveal the patterns or tendencies of functions. The found pattern can in turn aid in taking special directional/constraint limit or substitutions in CAS to get the analytical formula for SF bound when it is possible.

The overall lower bound we find of SF is presented in Section 2, through the triple analytical, graphical, and numerical methods. Followed by in-detail studies of SF of various selected distribution families, which are either widely used in practice, such as the Beta distribution in Section 3 and the generalized Gamma distribution in Section 6, or is most simple to simulate, such as the Kumaraswamy distribution in Section 4. The most inclusive distribution, BetaPrime distribution, is in Section 5, for which we do not get an analytical formula, so the empirical formula for SF lower bound is provided. Some distributions that have wide matching capabilities, but for the NATH may have fitted distribution facing numerical difficulties, such as the Fleishman distribution, whose fit has non-monotonically increasing polynomial form and hence is hard to solve for inverse CDF, are only briefly mentioned in Section 7. The top distribution found through maximum likelihood fit, the generalized hyperbolic distribution (GH), even with the most complex PDF, has unexpectedly simple and beautiful analytical formulas for SF lower bound; the results are in the final Section 8. All our studies will focus on SF bound deductions and applications.


2. Lower bound of the shape factor

For a random variable f with mean m f , the following characteristics are defined:

  • Moment (M), M r f r , r > 0 ,

  • Absolute Moment (AM), AM r f r , r > 0 ,

  • Central Moment (CM), CM r f m f r , r > 0 ,

  • Absolute Central Moment (ACM), ACM r f m f r , r > 0 ,

  • Skewness (S), S CM 3 CM 2 3 2 ,

  • Kurtosis (K), K CM 4 CM 2 2 ,

  • Shape Factor (SF), SF K S 2 = CM 4 CM 2 CM 3 2 .

We can prove by Hölder inequality (ölder's_inequality) that.

SF 1 :

f m f 3 f m f 3 = f m f 2 f m f 1 E1
f m f 4 1 2 f m f 2 1 2 . E2

A better inequality K S 2 + 1 is proved in [5, 6, 7]. But by Hölder inequality we can also know that ACM 4 ACM 2 ACM 3 2 = 1 iff f is constant: if f is not constant, the shape factor must be larger than the lower bound 1.

The contribution to SF > 1 plausibly comes from two parts: Eq. (1) due to symmetry, the more symmetric the distribution, the larger the contribution to SF, or conversely, the smaller the SF, the more asymmetric the distribution; and Eq. (2) due to ACM convexity or steepness, the steeper the PDF, the smaller the SF.

This property of the shape factor identified our exceptional perils as possessing very steep and asymmetric PDF whose SF are small.

2.1 Are there better definitions of shape factor?

To measure the steepness or the convexity, we can get similar inequality to Eq. (2) by Hölder inequality for absolute moment:

AM 4 AM 2 AM 3 2 1 and AM 3 AM 1 AM 2 2 1 .

From absolute central moment define:

SF 1 ACM 4 ACM 2 ACM 3 2 1 and SF 2 ACM 3 ACM 1 ACM 2 2 1 .

For nonnegative random variables such as the reinsurance contract loss distribution, use the following inequality for moment:

M 4 M 2 M 3 2 1 and M 3 M 1 M 2 2 1 .

From another application of Hölder inequality, we get yet other measures of convexity from absolute central moment:

ACM r ACM s r s , where 0 < r < s ,

SF 3 r ACM r ACM 1 r 1 , where 0 < r < 1 and SF 3 s ACM s ACM 1 s 1 , where s > 1 .

Similar definition from absolute moment:

AM r AM s r s , where 0 < r < s ,

SF 4 r AM r AM 1 r 1 , where 0 < r < 1 and SF 4 s AM s AM 1 s 1 , where s > 1 .

Checking against NormalDistribution μ σ , we see their minimum based on absolute moment: AM 4 AM 2 AM 3 2 , AM 3 AM 1 AM 2 2 , and SF 4 2 = AM 2 AM 1 2 , are all 1, but that by absolute central moment are not: min SF1 = 1.1781, min SF2 = 1.27324, min SF3 [2]=1.5708. Moreover, the convex index SF1, SF2, and SF3 out of absolute central moment are shift invariant besides the scale transformation invariant of the random variable, so they are preferred over the ones based on absolute moment.

The only case in favor of M 4 M 2 M 3 2 and M 3 M 1 M 2 2 is when the numerical calculation error with extreme parameters arrive at negative kurtosis, then the calculated SF are meaningless (An example of BesselK function inaccuracy brings about negative kurtosis for generalized hyperbolic distribution can be found in [8]).

Even though both SF and SF 1 are invariant under linear transformation of the distribution, and both measure the convexity, SF SF 1 can additionally measure the asymmetry, combining these two into one quantity. Since most distributions in reinsurance are not symmetric, SF is preferred over SF 1 . That only SF measured both asymmetry and convexity, while the others cannot, can be seen from Figure 1, for the case of exponential distribution family with PDF e x n n x 1 + n , x 0 , n > 0 , which is WeibullDistribution n 1 or GammaDistribution 1 1 n 0 , where only SF has a nontrivial interior global minimum.

Figure 1.

SF SF 1 SF 2 SF 3 2 plot of exponential distribution e x n n x 1 + n . The horizontal axis is the order of the exponential and the vertical axis is shape factors values.

An intuitive reason for why using shape factor SF in favor of skewness and kurtosis alone is provided by studying the simple example power distribution family with PDF n + 1 n x 1 n , x 0 1 , n < 1 n > 0 (or BetaDistribution 1 / n + 1 1 ). This distribution family has the largest value of skewness and kurtosis, and at the same time the smallest shape factor SF when n turns to −1, where the PDF is the steepest, but the skewness and kurtosis take the indistinguishable value of infinity. In comparison, the shape factor SF takes the finite and distribution family specific value of 1.125. The shape factor SF thus makes meaning out of the meaningless infinities.

2.2 Alternative way of defining shape factor for symmetric distribution

For symmetric distribution, CM 3 = 0 , our SF will be indiscriminately infinity. We can now employ SF 1 in place of SF . Other measures from ACM such as SF 2 and SF 3 may also be candidates. From the SF 3 plot Figure 2 of NormalDistribution μ σ we see that min 0 < r < 1 SF 3 r = 0.919824 . The lower the value of SF 3 2 , the higher the min 0 < r < 1 SF 3 r . We can use either SF 3 2 or min 0 < r < 1 SF 3 r as a shape factor for symmetric distribution to describe the convexity of the ACM curve. The second measure has the merit of independence to the power order r , by additional efforts of numerical minimization. For our power distribution family, the maximum of the minimum is: max n > 0   min 0 < r < 1 SF 3 r = 0.942085 , higher than the Normal distribution family.

Figure 2.

SF 3 plot of Normal distribution. The horizontal axis is the order r of the power and the vertical axis is SF 3 r .

When all SF , SF 1 , and SF 2 are available, however, we will prefer SF to SF 1 and SF 2 since its dependency on parameters show simpler patterns than the other two; this can be shown from their contour plots for Beta distribution Figures 3, 4, 5, where SF contours are almost lines.

Figure 3.

SF 1 contour plot of Beta distribution. The horizontal axis is α and the vertical axis is β .

Figure 4.

SF 2 contour plot of Beta distribution. The horizontal axis is α and the vertical axis is β .

Figure 5.

SF contour plot of Beta distribution. The horizontal axis is α and the vertical axis is β .

2.3 Lower bound of SF for well-known distributions

Using numerical optimization [9, 10], for most of the top-fitted distributions from the maximum log likelihood approach, we get the minimum SF values, with distribution definition in [11] whose naming and parameterization for probability distributions will be used throughout this chapter, in Table 2.

Distribution Min SF Location of the Min
FrechetDistribution α β μ 2.9555 α → 7.9305
ExtremeValueDistribution α β 4.15843 any α, β
MaxStableDistribution μ σ ξ 1.91227 ξ → -1.55970090120176
InverseGaussianDistribution μ λ θ 1.5 λ/μ → 0
SkewNormalDistribution μ σ α 3.90603 α → ∞
ExpGammaDistribution κ θ μ 2.25 κ → 0
BirnbaumSaundersDistribution α λ 1.63481 α → ∞
MeixnerDistribution a b m d 1.5 d → 0,b → ±π

Table 2.

Lower bound of SF for some well-known distributions.

From this table, we know that most of the distributions are not able to describe NATH since NATH has SF 1.834. More involved numerical integration and optimization also eliminated the Beckmann Distribution [12], with admissible SF range 3.63–8.16, being the top four-parameter-distribution in another distribution fit case study that has SF 4.58.

The Alpha-Skew-Normal Distribution from [13] has minimum SF 4.95061 when α is 2.07764, from its proposition 2.3, is thus also not eligible for NATH.

The global lower bound of SF for parametric distribution can be used to filter out those distributions whose values are larger than the losses data SF , so that we can focus on distributions that do not violate the bound. In the following sections we will study typical distribution SF bound, beginning with the Beta distribution.


3. Beta distribution

Regardless of the fact that multitude distribution types have been used for the frequency and severity distribution of individual contract losses, the aggregated portfolio losses for the majority of perils can be fitted by a compound Poisson distribution with Beta distribution as the severity, somehow an attest of its prevalence. Beta distribution has min SF = 1.0 , so we need an in-detail study of why it cannot fit NATH.

When matching a BetaDistribution α β for skewness 5.99378 and kurtosis 65.8902, we must have β < 0. When matching a Beta distribution for CV(=std/mean, the standard deviation divided by the mean) 1.2829 and either skewness 5.99378 or kurtosis 65.8902, we must have either both α < 0 and β < 0 or at least one of α or β less than 0. Since CV, skewness, and kurtosis are scale invariant, so no scaled Beta distribution can at the same time match any two of the three statistics CV, skewness, and kurtosis.

3.1 Minimum shape factor for given CV

Using Mathematica, we can solve the parameter α and β by cv and std for BetaDistribution α β :

α cv std cv 2 std cv 3 , β cv 2 2 cvstd cv 3 std + std 2 + cv 2 std 2 cv 3 std .

Since α > 0, we must have:

std < cv 1 + cv 2 ,


1 1 4 std 2 2 std < cv < 1 + 1 4 std 2 2 std .

We also know std must be between 0 and 0.5 for these solutions to exist. By computer-aided exploration through contour plot, we can find the location of the std where SF takes minimum for a given cv.

The overall observation is that when cv < 1, SF approaches infinity in the middle value of std, and decreases when deviating from it. When cv > 1, SF approaches its minimum in the middle value of std and increases when departing. Together with the fact that std has an allowable upper bound of cv/(1 + cv^2) and lower bound of 0, the minimum of SF must be attained either at the global extreme where the derivative of SF with respect to std is zero or at the two boundaries when cv > 1, and attained at the two boundaries when cv < 1.

Using Mathematica to take the derivative of the shape factor with respect to std to find the std where shape factor attained extreme values, and solving it for the intersection with std upper and lower bound, we know the minimal shape factor for Beta distribution for a given CV when CV is below 0.707107 or above 2.48239 (intersecting std upper bound) is attained at std upper bound cv 1 + cv 2 with value:

min 0 < std < cv 1 + cv 2 SF = 1 cv 2 + cv 4 1 + cv 2 2 , when cv < 0.707107 cv > 2.48239 . E3

When CV is between 0.707107 and 1.024766 (intersecting std lower bound) the minimal shape factor is attained at std lower bound 0 with value:

min 0 < std < cv 1 + cv 2 SF = 1.5 + 0.75 cv 2 , when 0.707107 cv 1.024766 . E4

When CV is between 1.024766 and 2.48239, the minimum SF is attained at std that is the zero derivative points of the shape factor. The piecewise curve plot of the minimum SF for given CV is in Figure 6. The formula for the central piece, minshape, is given in Figure 7 which is too complex for manual derivation without the aid of computer algebra system.

Figure 6.

Plot of Beta distribution min shape factor for given cv.

Figure 7.

Formula for minshape obtained using Mathematica.

From the curve we know when CV = 1.28, the minimal shape factor is 1.88, larger than 1.83 of NATH. In the best effort to match the input, we may elect to relax CV, for example, to 1.3, then the minimum shape factor is 1.85. With the constraint of a given CV, the minimum shape factor of the Beta Distribution may be significantly larger than its global minimum 1, so that it cannot attain to the wanted SF value.

3.2 Shape factor range for given skewness

By solving Beta distribution parameters α and β through skewness sk and kurtosis kt , and examining the contour plot of β, we can see the allowable region is bound by two parabolas, Figure 8.

Figure 8.

Contour plot of Beta distribution β parameter. The horizontal axis is the skewness and the vertical axis is the kurtosis.

For a fixed skewness, α is monotonic increasing with respect to kurtosis; on the other hand, β has a singular point in some kurtosis, below that kurtosis is positive and monotonic increasing(in the region where α is positive), Figure 9.

Figure 9.

Plots of Beta distribution β parameter and α parameter vs. kurtosis for a given skewness 5.99378.

Solving for that singular point we get the permissible kurtosis upper bound 3 + 3 2 sk 2 , and solve for β = 0 get the permissible kurtosis lower bound 1 + sk 2 .

Observe that the upper bound is when β turns to infinity, we can also get a simpler derivation of the upper bound by representing skewness and shape factor in α and β , letting β , and then eliminating α to get shape factor as a function of skewness (Mathematica cannot solve equation for skewness which includes square root expression, we get around that by solving equation for the square of skewness, and then abandoning the negative solution introduced by this square).

A third way of more tedious calculation is through solving α by skewness and β , substituting the real solution into shape factor, and then take the limit for β .

All three methods get the same upper bound of SF = 3 2 + 3 sk 2 .

So for Beta distribution, the allowable region of skewness and kurtosis is bound below by kurtosis = skewness^2 + 1 where β 0 , and above by kurtosis = 3 + 1.5*skewness^2 where β :

1 + 1 S 2 SF 1.5 + 3 S 2 . E5

For the given skewness of 5.99378 of NATH, the maximum allowable kurtosis is 56.88813, less than the wanted 65.8902. So NATH cannot be fitted by any affine transformation of Beta distribution, certifying NATH as a trying case for distribution fitting. We will use it to test many of the well-known distributions in later sections. We also see surprisingly that unlike many of the other distribution families whose shape factors are too high, the Beta distributions have the shape factor range too low, or too close to 1. This suggests us to search for distributions with shape factors ranges in between.


4. Kumaraswamy distribution

Using the same approach as in the Beta distribution, we first study the skewness and kurtosis tendency of KumaraswamyDistribution α β [14], since the latter tested to be a better choice in our experiment and is also the easiest for simulation, Figures 10, 11, 12; and then study the SF bound for given skewness.

Figure 10.

Contour plot of Kumaraswamy distribution skewness.

Figure 11.

Contour plot of Kumaraswamy distribution kurtosis.

Figure 12.

Contour plot of Kumaraswamy distribution shape factor.

From these plots, we see an overall rough tendency of the skewness, kurtosis and shape factor. For a given α, the shape factor converges to a finite limit when β . For a given skewness or a given kurtosis, there exists a maximum allowable α that is arrived when β . In the parameters space of (α,β), for a given α, the kurtosis is increasing with respect to β in the top left portion where the skewness is positive, and decreasing in the right bottom portion where the skewness is negative. And in the parameters space of (α,β), for a given α, the shape factor is decreasing with respect to β in the top left portion where the skewness is positive, and increasing in the right bottom portion where the skewness is negative. But we will see later that the tendencies are more delicate than the monotonicity shown through visual observation.

Combining the tendency of shape factor and the contour plot for given skewness, kurtosis, and shape factor as in Figure 13, we may guess that for a given positive skewness, when α turn to its upper limit and β turn to infinity, the shape factor will converge to its minimum. We use Mathematica to calculate the asymptotic expansion of the Gamma function and the quotient of Gamma function at infinity for orders up to 4 or 2, take the symbolic limit for β , to get these boundary values, Figures 14 and 15.

Figure 13.

Contour plot of Kumaraswamy distribution skewness, kurtosis, and shape factor for given values 5.99, 65.89, and 1.83. The horizontal axis is the α parameter and the vertical axis is the β parameter.

Figure 14.

Derivation of Kumaraswamy distribution skewness upper bound for given α.

Figure 15.

Derivation of Kumaraswamy distribution kurtosis upper bound for given α and shape factor boundary value for given α when β → ∞.

We thus have a simple formula for boundary value of Kumaraswamy distribution shape factor:

limit β S = 2 Gamma 1 α 3 6 α Gamma 1 α Gamma 2 α + 3 α 2 Gamma 3 α α α Gamma 1 + 1 α 2 + 2 Gamma 2 α 3 / 2 , E6
limit β K = 3 Gamma 1 α Gamma 1 α 3 4 α Gamma 1 α Gamma 2 α + 4 α 2 Gamma 3 α + α 4 Gamma 4 + α α Gamma 1 α 4 4 α Gamma 1 α 2 Gamma 2 α + α 4 Gamma 2 + α α 2 , E7
limit β K S 2 = α 3 α Gamma 1 + 1 α 2 + 2 Gamma 2 α 3 3 Gamma 1 α Gamma 1 α 3 4 α Gamma 1 α Gamma 2 α + 4 α 2 Gamma 3 α + α 4 Gamma 4 + α α 2 Gamma 1 α 3 6 α Gamma 1 α Gamma 2 α + 3 α 2 Gamma 3 α 2 Gamma 1 α 4 4 α Gamma 1 α 2 Gamma 2 α + α 4 Gamma 2 + α α 2 . E8

Its plot Figure 16 has two branches, the dividing point is α 3.602349425719043 where the skewness is zero, and below it is mainly the positive skewness region while above it is the negative skewness region.

Figure 16.

Plot of Kumaraswamy distribution shape factor boundary value for given α when β → ∞.

The minimum value at the left branch of Figure 16 is 1.91227 and arrived at α = 0.641149. When α > 1000 the numerical value for that boundary can be negative and is thus unreliable. The value 1.91227 is not the global minimum of the shape factor: for α = 0.641149 the shape factor plot Figure 17 with respect to β decreases first, at the point 10.6095 arriving at the minimum value of 1.80935, and increasing after the point 10.6095.

Figure 17.

Plot of Kumaraswamy distribution shape factor for given α = 0.641149, β in the range 0.3–1 and 1–300.

In principle, the extreme value of the shape factor for a given skewness will arrive either at the upper boundary where β → ∞ or at the lower boundary where α → 0, or at some middle point where the contour plot of the skewness and the contour plot of the kurtosis will be tangent to each other. The Mathematica contour plot does not work for a very small α, but by numerical minimization we know the global minimum of the Kumaraswamy distribution shape factor is 1.03709 when α = 1.80143 10 9 , β = 0.247044 . The conditional minimum of the shape factor when skewness = 5.99378 is about 1.04753 when α = 10 10.5 , β = 0.149286 through list calculation; this is higher than 1 + 1/S^2 = 1 + 1/5.99378199789956^2 = 1.02784, the lower boundary of Beta distribution.

The Mathematica contour plot works for large α, and we see the shape factor is increasing along the contour of skewness, which attains its maximum when β → ∞. For example, for NATH skewness 5.99378199789956, the maximum shape factor is 1.97131, arriving at α = 0.5239510562868946. The maximum shape factor of Kumaraswamy distribution for given skewness is in Figure 18, which is algebraically represented by the parametric curve of Eq. (6) and Eq. (8).

Figure 18.

Plot of Eq. (6)(8) and plot of Kumaraswamy distribution maximum shape factor for given skewness.

So the permissible shape factor range of the Kumaraswamy distribution still spans the lower end of the whole allowable range of (1,∞), but higher than that of the Beta distribution. Affine transformed Kumaraswamy distribution can fit all the first four moments of NATH, with the fitted distribution TVaR close to NATH TVaR in the error range of 5–6%, while the best effort affine transformed Beta distribution is in the error range of 9–10%.

To further improve the fit, we need additional freedom in parameters, such as the GB 1 distribution [15], since KumaraswamyDistribution α β GeneralizedBetaDistributionI 1 β α 1 , and the maximum shape factor plot in Figure 18 is lower than that of LogNormalDistribution, the upper bound of GB 1 . The following section will study a sibling distribution to GB 1 , fitted as good as GB 1 , but is more widely known.


5. BetaPrime distribution

Beta distribution and Kumaraswamy distribution are a few exceptions which have analytical formulas for the shape factor bounds; for other distributions to be studied, numerical optimization and empirical plot or formula will be the only feasible approach.

Transformation of Beta Distribution by x/(1-x) is the GB 2 ([15]), or BetaPrimeDistribution p = α q = β α = 1 β = 1 ([11]): TransformedDistribution x 1 x x BetaDistribution α β BetaPrimeDistribution α β . The minimum shape factor of Beta Distribution is 1, but that of the transformed is 1.5:

NMinimize 3 3 + β 2 1 + β 2 + α 2 5 + β + α 1 + β 5 + β 4 4 + β 1 + 2 α + β 2 α > 0 β > 4 α β = 1.5000000239052607 α 0 . β 6.274769836372949 × 10 7 .

Empirically, the larger the third parameter α, the smaller the minimum shape factor. The smallest shape factor we get of the BetaPrimeDistribution is 1.125, when α = 446.49537:

FindMinimum Kurtosis BetaPrimeDistribution p q α β Skewness BetaPrimeDistribution p q α β 2 / . α 1 x / . q 4 x +   y / . x 10 z 1 . > p > 0 . y > 1 . 4 . < z < 1 . p 6.384125235007732 × 10 10 y 1.0032844709998097 z 2.157370895027263 MaxIterations 5000 = 1.1250258984236121   p 2.083731454230264 × 10 8 y 42:816363091057056   z 2.6498169598310573 .

This is the same value as the minimum shape factor for GammaDistribution α β γ μ (in Section 6). When α > 10,000, the Gamma function involved will not calculate or will calculate incorrectly.

With the transformation of p-> 10^w, α-> 10^-z, q-> 4*10^z + y, we can study the GB 2 shape factor change tendency with respect to α, Figure 19, and shape factor change tendency with respect to p, Figure 20.

Figure 19.

GB 2 distribution shape factor vs. α for fixed p = 10^-3.312 = 0.000487528.

Figure 20.

GB 2 distribution shape factor vs. p for fixed α = 10^2.6498169598310573 = 446.495.

The GB 2 shape factor is mainly determined by α and p, only slightly changing with respect to q when q is smaller than 5. The change with respect to α and p is similar, having two peaks, or three peaks if we regard the two sides of the infinity as two branches since that border is not easily crossable for searching or optimization algorithms.

GB 2 shape factor’s dependency with p and α, or w and z through transformation p = 10^w, α = 10^-z, is mostly unaffected by q except for right-most values of z. They are μ-shaped (Figure 21), this is different from Hyperbolic Distribution (in Section 8), whose shape factor dependency with λ is V-shaped. We guess V-shaped curves have unique global minimums, but μ-shaped curves will show bifurcation behavior: the converged solution in optimization will be very different when the initial point or interval is slightly different.

Figure 21.

GB 2 distribution skewness kurtosis and shape factor vs. α or z, vs. p or w plots for fixed y = q-4/α.

The knowledge that the shape factor curve attained extreme values in −3.3,-1.25 and 1 with respect to z, and attained extreme values in −2.65, −1.11 and 1 with respect to w, can be used to set the initial interval, the paramount factor determining the quality of the numerical optimization solution, for solving the GB 2 fitting problem.

5.1 Minimum shape factor for GB 2

The skewness and kurtosis matching problem for GB 2 is very sensitive to the initial parameter ranges given. A study of the minimum shape factor of GB 2 with respect to each parameter will give us permissible ranges for those parameters. Direct work with shape factor encounters problems from Mathematica’s numerical optimization function NMinimize , minimizing the log shape factor instead can overcome this difficulty. The plot is in Figure 22.

Figure 22.

The numerical minimum GB 2 shape factor for given p in horizontal axis.

In the range (0.0001, 5.0) of p , the numerical minimum shape factor plot of GB 2 is a very smooth curve. The fitted formula of GB 2 min SF for given p by Mathematica’s machine learning function FindFormula is Eq. (9).

min K S 2 = 1.1593871374775397 + 1.4702458297305288 0.5148499158800361 1 p 0.3215433282777008 E9

As a test, for NATH the log shape factor is Log 1.83408 = 0.60654412 , the solution of Eq. (9) for p with NATH SF is p = 0.608342 ; the minimum log shape factor of GB 2 for this p is 0.60603997, only 0.08% smaller than input.

From the contour plot Figure 20 we know for given α, the shape factor of GB 2 has two singular points with p or 10 w . The minimization for given α needs to carry out in each of the three regions cut by these two singular points. The plot is in Figure 23. With a new parameterization, p = λ α , q = 4 + ν α , the minimization of shape factor for GB 2 , for given λ = , is easier to perform. The plot is included in Figure 23 as well.

Figure 23.

The numerical minimum GB 2 shape factor for given α or given pα in horizontal axis.

Figures 22 and 23 show that the permissible parameters for NATH are p < 0.63 , α > 0.5 , < 0.5 . This is confirmed by GB 2 fit practice. The best fit by GB 2 for NATH is at w = 0.329075005 , p = 0.468732 , with about 5% error from input TVaR. The discontinuity of fitted GB 2 TVaR with respect to parameter change is also observed, this w value is such a critical point.


6. Generalized gamma distribution

The generalized gamma distribution in Mathematica is the Amoroso distribution [16], with the parameter correspondence: α α , β θ , γ β , μ a .

For generalized gamma distribution GammaDistribution α β γ μ , the shape factor depends only on α and γ . It seems the smaller the α , and the bigger the γ , the smaller the K S 2 . When α = 3.318512677036329 × 10 12 , γ = 8811.572418686921 , K S 2 = 1.125 , close to the global minimum 1 of K/S^2.

So there arises the question: the generalized gamma and GB 2 can match smaller shape factors than Hyperbolic Distribution (Section 8), why they cannot fit as good as the latter for NATH with shape factor 1.83409?

One explanation is that the numerical solution for GB 2 or generalized Gamma distribution is trapped in the shape factor curve right branch by the combined constraints of skewness and kurtosis, which is not the branch that can attain 1.125, unlike the generalized hyperbolic distribution whose shape factor has a global minimum in λ = 0.


7. Fleishman distribution

We guess 1.5 is the lower bound of shape factor for most unbounded parametric distribution families. For example, for Fleishman distribution, by the empirical formula from [5], γ 4 > 1.738 γ 3 2 0.3544 γ 3 + 1.978 , the minimum shape factor is 1.72213, larger than 1.5.

The lower bound of shape factor from unbounded distributions seems, in general, to be higher than bound distributions’. Outside of the latter’s upper bound and near the former’s lower bound, for a SF value slightly larger than 1.5, in practice, most parametric distributions have difficult matching both the kurtosis and skewness: the comparatively best one is selected for study in the next section.


8. Hyperbolic distribution

Taking a sequence of numerical minimization of the shape factor, for various values of fixed λ, we get the empirical minimum shape factor curve for generalized hyperbolic distribution (GH), HyperbolicDistribution λ α β δ μ , in Figure 24.

Figure 24.

V-shape of the numerical minimum GH shape factor for given λ in horizontal axis.

We observed that when λ > −0.6, the minimum shape factor is attained when α^2-β^2-> 0 and β-> 0, that is, it is attained by a skew hyperbolic t distribution [17, 18, 19]. When looking at the plot of shape factor with respect to λ, we feel that it must have some simple formula. So we utilize Mathematica symbolic calculation to expand the shape factor with asymptotic expansion for BesselK λ a , or K λ a in [20], with respect to α^2-β^2 and then take the symbolic limit, Figure 25.

Figure 25.

Derivation of the GH shape factor limit when λ > −2.

The semi-empirical formula for the minimum shape factor in this range thus obtained is very simple, Eq. (10–11), which has the global minimum of 1.5 when λ turns to 0.

min α , β , δ , μ SF = 1.5 + 0.75 λ , when λ 0 , E10
min α , β , δ , μ SF = 1 + 1 2 + λ , when 0.6 λ 0 E11

When λ ≤ −0.65, however, the minimum shape factor is not attained when α^2-β^2-> 0. When λ is in the interval [−9,−0.65], the attainable smallest shape factor is between 3.15 and 1.74, with an empirical 10th order polynomial formula Eq. (12), or less accurately a mixed exponential and power function Eq. (13), found through the Mathematica FindFormula .

min α , β , δ , μ SF = 1.1130471668735116 1.6512030619809768 λ 1.6137376956833365 λ 2 1.1485038172210114 λ 3 0.5421785615853132 λ 4 0.17094578834265223 λ 5 0.03603744794749387 λ 6 0.005000441043297472 λ 7 0.0004372189547557593 λ 8 0.000021791071048963054 λ 9 4.711954312790356 × 10 7 λ 10 E12
min α , β , δ , μ SF = 2.2104215691249425 0.6522131009473879 1.6355318649123258 λ + 0.018965779149540653 λ 3 0.1051542360603726 λ E13

So for each given K/S^2 value, there exists a permissible interval of λ, whose lower bound is calculated via Eq. (11–12) and upper bound is calculated via Eq. (10). When λ changes inside this interval, we noticed that the 0.99 TVaR of the first four moments matched generalized hyperbolic distribution will increase with respect to λ. If the lower bound still has 0.99 TVaR bigger than the input TVaR, then it is not possible to fit with moments matched HyperbolicDistribution . The opposite statement is also valid for the interval upper bound.

With this knowledge, the NATH permissible λ interval is [−0.8439,0.4454], and the left end point still have 0.99 TVaR larger than the input TVaR, but now only by 4.05%, better than the 5% error of GB 2 .


9. Conclusion and discussions

We proposed using the ratio of kurtosis by squared skewness as the best candidate for shape factor that can characterize the distribution asymmetry, as well as the PDF steepness. The closer this factor to 1, the more asymmetry and the steeper the PDF. The asymptotic approximation and symbolic limit is used to calculate the boundary of this factor for various distributions: the Beta, the Kumaraswamy, and the Hyperbolic Distributions, for example. This range information of the shape factor, with the surprisingly simple formulas in the three above examples (Eqs. 58, 10, 11), can be used to select or eliminate candidate distributions for fitting. The plot of the shape factor together with plot for skewness and kurtosis can aid in setting the initial value or parameter intervals when fitting distribution to data by numerical optimization, which usually would not work well without this information.

The idea of the shape factor and the careful study of each distribution for this shape factor is the preliminary for the numerical optimization that finally finds the best fit. The information provided by shape factor plot is rough but the numerical optimization’s dependence on initial value or intervals is delicate, exemplified by GB 2 case. The optimization function NMinimize and FindMinimum in Mathematica sometimes can only find a local optimum at best. As shown in [21, 22], the DyHF and the CMODE algorithms are the two best no-adjustment-needed global optimization algorithms. Now that the C2oDE algorithm is better than these two [23], it would be desirable to see how it works on the GB2 fit problem. With a foolproof universally applicable global optimization algorithm, the ado with shape factor and their boundaries will no longer be needed, or be used merely as some validations; but before that time, the hard earned knowledge about shape factor through CAS is still indispensable. This is a good topic for subsequent research.

Our shape factor idea is only a small step ahead of the skewness-kurtosis plot of Pearson [6] and McDonald et al [15, 24, 25, 26]. Or we just made the idea implicitly in their plot explicit. But with this clearly defined form, anyone can readily start calculating it for any interested candidate distribution.

Our formula Eq. (5) is not new, since Beta distribution has the same range of skewness, kurtosis, and shape factor as the scaled Beta distribution, the B 1 distribution in [15]. Our presentation is an example of how our method can be used to easily arrive at those formulas. Theoretically equivalent expressions are not equivalent in application. With data distributions usually not having small skewness, Eq. (5) says that the Beta distribution has a shape factor roughly in the range of (1, 1.5), this not only reveals an intrinsic property of Beta distribution, but is also more easily applicable in practice than the skewness-kurtosis plot.

The residual error of all the distributions tested so far indicates that the power function or simple exponential function PDF is not enough to provide the additional freedom of shifting for the EP curve on the condition of matched first four moments. Other forms such as mixtures, combinations, or transformations of distributions may need to be considered. A previous study indicated the following transformations are good candidates [4, 27, 28, 29, 30, 31, 32]: EWGU, KGG, EG, EWED, LIG, THT. Further research will be done along these lines.



This research is supported by Validus Research Inc. The author thanks his colleagues and former colleagues for helpful suggestions and feedbacks, the editor for revision suggestions, and Claire Wang for help with grammar and language corrections.

Conflict of interest

The author declares no conflict of interest.


  1. 1. Sharma K. Natural Catastrophe Modeling for Pricing in Insurance [Thesis]. Tartu: University of Tartu; 2014
  2. 2. Currie ID. Maximum likelihood estimation and mathematica. Applied Statistics. 1995;44(3):379-394
  3. 3. Wang F. Dfittool for Mathematica [Internet]. 2016. Available from: [Accessed: 2018-11-09]
  4. 4. Wang FX. An inequality for reinsurance contract annual loss standard deviation and its application. In: Salman A, Razzaq MGA, editors. Accounting from a Cross-Cultural Perspective. IntechOpen; 2018. pp. 73-89. DOI: 10.5772/intechopen.76265 Available from:
  5. 5. Ferenci T. The Use of Fleishman Distribution in the Empirical Investigation of Statistical Tests [Thesis]. Budapest University of Technology and Economics; 2012. Available from: [Accessed: 2018-05-02]
  6. 6. Pearson K IX. Mathematical contributions to the theory of evolution.—XIX. Second supplement to a memoir on skew variation. Published 1 January 1916. DOI: 10.1098/rsta.1916.0009. Available from: [Accessed: 2018-04-23]
  7. 7. Klaassen CAJ, Mokveld PJ, van Es B. Squared Skewness minus kurtosis bounded by 186/125 for unimodal distributions. Statistics & Probability Letters. 2000;50(2):131-135
  8. 8. Wang F. Problem with BesselK Function [Internet]. 2018. Available from: [Accessed: 2018-11-12]
  9. 9. Wolfram Mathematica Tutorial Collection. Constrained Optimization. Wolfram Research, Inc; 2008
  10. 10. Loehle C. Global optimization using mathematica: A test of software tools. Mathematica in Education and Research. 2006:139-152
  11. 11. Marichev O, Trott M. The Ultimate Univariate Probability Distribution Explorer [Internet]. 2013. Available from: [Accessed: 2018-06-06]
  12. 12. Hill RJ, Frehlich RG. Probability distribution of irradiance for the onset of strong scintillation. Journal of the Optical Society of America. A. 1997;14(7):1530-1540. DOI: 10.1364/JOSAA.14.001530
  13. 13. Olivero DE. Alpha-skew-normal distribution. Proyecciones Journal of Mathematics. 2010;29(3):224-240. Available from: [Accessed: 2018-05-01]
  14. 14. de Pascoa MAR, Ortega EMM, Cordeiro GM. The Kumaraswamy generalized gamma distribution with application in survival analysis. Statistical Methodology. 2011;8(5):411-433. DOI: 10.1016/j.stamet.2011.04.001 Available from:
  15. 15. McDonald JB, Sorensen J, Turley PA. Skewness and kurtosis properties of income distribution models. LIS working paper series, No. 569, 2011. Review of Income and Wealth. 2011. DOI: 10.1111/j.1475-4991.2011.00478.x Available from:
  16. 16. Crooks GE. The Amoroso Distribution [Internet]. 2010. Available from: [Accessed 2018-05-03]
  17. 17. Aas K, Haff IH. The generalized hyperbolic skew Student’s t-distribution. Journal of Financial Econometrics. 2006;4(2):275-309. DOI: 10.1093/jjfinec/nbj006
  18. 18. Hu W, Kercheval A. Risk management with generalized hyperbolic distributions. In: Locke P, editor. Proceedings of the Fourth IASTED International Conference on Financial Engineering and Applications (FEA '07). Anaheim, CA, USA: ACTA Press; 2007. pp. 19-24
  19. 19. Scott DJ, Würtz D, Dong C, Tran TT. Moments of the generalized hyperbolic distribution. Computational Statistics. 2011;26(3):459-476. DOI: 10.1007/s00180-010-0219-z
  20. 20. Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York: Dover; 1972
  21. 21. Wang FX. Relay Optimization Method [Internet]. May 2014. Available from: [Accessed 2018-07-28]
  22. 22. Wang FX. Design index-based hedging: Bundled loss property and hybrid genetic algorithm. In: Tan Y, Shi Y, Buarque F, Gelbukh A, Das S, Engelbrecht A, editors. Advances in Swarm and Computational Intelligence. ICSI 2015. Lecture Notes in Computer Science, vol 9140. Cham: Springer; 2015. pp. 266-275. DOI: 10.1007/978-3-319-20466-6_29
  23. 23. Wang BC, Li HX, Li JP, Wang Y. Composite differential evolution for constrained evolutionary optimization. IEEE Transactions on Systems, Man, and Cybernetics: Systems. . DOI: 10.1109/TSMC.2018.2807785
  24. 24. Vargo E, Pasupathy R, Leemis LM. Moment-ratio diagrams for Univariate distributions. Journal of Quality Technology. 2010;42(3):1-11
  25. 25. Celikoglu A, Tirnakli U. Skewness and kurtosis analysis for non-Gaussian distributions. Physica A: Statistical Mechanics and its Applications. 2018;499:325-334. DOI: 10.1016/j.physa.2018.02.035
  26. 26. Huerlimann W. Normal variance-mean mixtures (I) an inequality between skewness and kurtosis. Advances in Inequalities and Applications. 2014;2014(2)
  27. 27. Cordeiro GM, Ortega EMM, Ramires TG. A new generalized Weibull family of distributions: Mathematical properties and applications. Journal of Statistical Distributions and Applications. 2015;2(13). DOI: 10.1186/s40488-015-0036-6
  28. 28. Barreto-Souza W, Santos AHS, Cordeiro GM. The Beta generalized exponential distribution. Journal of Statistical Computation and Simulation. 2010;80:159-172
  29. 29. Lemonte AJ, Cordeiro GM. The exponentiated generalized inverse Gaussian distribution. Statistics & Probability Letters. 2011;81(4):506-517. DOI: 10.1016/j.spl.2010.12.016
  30. 30. Alzaghal A, Famoye F, Lee C. Exponentiated T-X family of distributions with some applications. International Journal of Statistics and Probability. 2013;2(3). DOI: 10.5539/ijsp.v2n3p31
  31. 31. Okorie IE, Akpanta AC, Ohakwe J, Chikezie DC, Shiraishi H. The modified power function distribution. Cogent Mathematics. 2017;4(1). DOI: 10.1080/23311835.2017.1319592
  32. 32. Borzadaran GR, Borzadaran HAM. Log-concavity property for some well-known distributions. Surveys in Mathematics and its Applications. 2011;6:203-219

Written By

Frank Xuyan Wang

Submitted: 12 September 2018 Reviewed: 01 December 2018 Published: 07 January 2019