What Determines EP Curve Shape?

Abstract

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.

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 $\prod _i{p_i}^{n_i}$ . It is easily seen that for the optimum solution we have $\frac{p_i}{p_j}=\frac{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\left[r\right]\equiv \int f^r\mathit{d\mu },r>0,$

• Absolute Moment (AM), $\mathit{AM}\left[r\right]\equiv \int {\left|,f\right|}^r\mathit{d\mu },r>0,$

• Central Moment (CM), $\mathit{CM}\left[r\right]\equiv \int {\left(,f-m_f\right)}^r\mathit{d\mu },r>0,$

• Absolute Central Moment (ACM), $\mathit{ACM}\left[r\right]\equiv \int {\left|,f-m_f\right|}^r\mathit{d\mu },r>0,$

• Skewness (S), $S\equiv \frac{\mathit{CM}\left[3\right]}{\mathit{CM}{\left[2\right]}^{\frac{3}{2}}},$

• Kurtosis (K), $K\equiv \frac{\mathit{CM}\left[4\right]}{\mathit{CM}{\left[2\right]}^2},$

• Shape Factor (SF), $\mathit{SF}\equiv \frac{K}{S^2}=\frac{\mathit{CM}\left[4\right]\ast \mathit{CM}\left[2\right]}{\mathit{CM}{\left[3\right]}^2}$ .

We can prove by Hölder inequality (https://en.wikipedia.org/wiki/Hölder's_inequality) that.

$\mathit{SF}\ge 1$ :

$$\mid \int {\left(,f-m_f\right)}^3\mathit{d\mu }\mid \le \int {\left|,f-m_f\right|}^3\mathit{d\mu }=\int {\left|,f-m_f\right|}^2{\left|,f-m_f\right|}^1\mathit{d\mu }$$ E1

$$\le {\left(\int {\left|,f-m_f\right|}^4\mathit{d\mu }\right)}^{\frac{1}{2}}{\left(\int {\left|,f-m_f\right|}^2\mathit{d\mu }\right)}^{\frac{1}{2}}.$$ E2

A better inequality $K\ge S^2+1$ is proved in [5, 6, 7]. But by Hölder inequality we can also know that $\frac{\mathit{ACM}\left[4\right]\ast \mathit{ACM}\left[2\right]}{\mathit{ACM}{\left[3\right]}^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 $\mathit{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:

$$\frac{\mathit{AM}\left[4\right]\ast \mathit{AM}\left[2\right]}{\mathit{AM}{\left[3\right]}^2}\ge 1\text{and}\frac{\mathit{AM}\left[3\right]\ast \mathit{AM}\left[1\right]}{\mathit{AM}{\left[2\right]}^2}\ge 1.$$

From absolute central moment define:

$$\mathit{SF}1\equiv \frac{\mathit{ACM}\left[4\right]\ast \mathit{ACM}\left[2\right]}{\mathit{ACM}{\left[3\right]}^2}\ge 1\text{and}\mathit{SF}2\equiv \frac{\mathit{ACM}\left[3\right]\ast \mathit{ACM}\left[1\right]}{\mathit{ACM}{\left[2\right]}^2}\ge 1.$$

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

$$\frac{M\left[4\right]\ast M\left[2\right]}{M{\left[3\right]}^2}\ge 1\text{and}\frac{M\left[3\right]\ast M\left[1\right]}{M{\left[2\right]}^2}\ge 1.$$

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

$$\mathit{ACM}\left[r\right]\le \mathit{ACM}{\left[s\right]}^{\frac{r}{s}},\mathit{\text{where}}0<r<s,$$

$\mathit{SF}3\left[r\right]\equiv \frac{\mathit{ACM}\left[r\right]}{\mathit{ACM}{\left[1\right]}^r}\le 1,\mathit{\text{where}}0<r<1$ and $\mathit{SF}3\left[s\right]\equiv \frac{\mathit{ACM}\left[s\right]}{\mathit{ACM}{\left[1\right]}^s}\ge 1,\mathit{\text{where}}s>1.$

Similar definition from absolute moment:

$$\mathit{AM}\left[r\right]\le \mathit{AM}{\left[s\right]}^{\frac{r}{s}},\mathit{\text{where}}0<r<s,$$

$\mathit{SF}4\left[r\right]\equiv \frac{\mathit{AM}\left[r\right]}{\mathit{AM}{\left[1\right]}^r}\le 1,\mathit{\text{where}}0<r<1$ and $\mathit{SF}4\left[s\right]\equiv \frac{\mathit{AM}\left[s\right]}{\mathit{AM}{\left[1\right]}^s}\ge 1,\mathit{\text{where}}s>1.$

Checking against $\mathit{\text{NormalDistribution}}\left[\mu ,\sigma \right]$ , we see their minimum based on absolute moment: $\frac{\mathit{AM}\left[4\right]\ast \mathit{AM}\left[2\right]}{\mathit{AM}{\left[3\right]}^2}$ , $\frac{\mathit{AM}\left[3\right]\ast \mathit{AM}\left[1\right]}{\mathit{AM}{\left[2\right]}^2}$ , and $\mathit{SF}4\left[2\right]=\frac{\mathit{AM}\left[2\right]}{\mathit{AM}{\left[1\right]}^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 $\frac{M\left[4\right]\ast M\left[2\right]}{M{\left[3\right]}^2}$ and $\frac{M\left[3\right]\ast M\left[1\right]}{M{\left[2\right]}^2}$ is when the numerical calculation error with extreme parameters arrive at negative kurtosis, then the calculated $\mathit{SF}$ are meaningless (An example of $\mathit{\text{BesselK}}$ function inaccuracy brings about negative kurtosis for generalized hyperbolic distribution can be found in [8]).

Even though both $\mathit{SF}$ and $\mathit{SF}1$ are invariant under linear transformation of the distribution, and both measure the convexity, $\mathit{SF}\ge \mathit{SF}1$ can additionally measure the asymmetry, combining these two into one quantity. Since most distributions in reinsurance are not symmetric, $\mathit{SF}$ is preferred over $\mathit{SF}1$ . That only $\mathit{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\in \left(0,\infty \right),n>0$ , which is $\mathit{\text{WeibullDistribution}}\left[n,1\right]$ or $\mathit{\text{GammaDistribution}}\left[1,1,n,0\right]$ , where only $\mathit{SF}$ has a nontrivial interior global minimum.

Figure 1.

$\mathit{SF}\mathit{SF}1\mathit{SF}2\mathit{SF}3\left[2\right]$ 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 $\mathit{SF}$ in favor of skewness and kurtosis alone is provided by studying the simple example power distribution family with PDF $\frac{n+1}{n}{x}^{\frac{1}{n}},x\in \left[0,1\right],n<-1\parallel n>0$ (or $\mathit{\text{BetaDistribution}}\left[1/n+1,1\right]$ ). This distribution family has the largest value of skewness and kurtosis, and at the same time the smallest shape factor $\mathit{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 $\mathit{SF}$ takes the finite and distribution family specific value of 1.125. The shape factor $\mathit{SF}$ thus makes meaning out of the meaningless infinities.

2.2 Alternative way of defining shape factor for symmetric distribution

For symmetric distribution, $\mathit{CM}\left[3\right]=0$ , our $\mathit{SF}$ will be indiscriminately infinity. We can now employ $\mathit{SF}1$ in place of $\mathit{SF}$ . Other measures from $\mathit{ACM}$ such as $\mathit{SF}2$ and $\mathit{SF}3$ may also be candidates. From the $\mathit{SF}3$ plot Figure 2 of $\mathit{\text{NormalDistribution}}\left[\mu ,\sigma \right]$ we see that ${\displaystyle \underset{0<r<1}{min}}\mathrm{SF}3\left[r\right]=0.919824$ . The lower the value of $\mathit{SF}3\left[2\right]$ , the higher the ${\displaystyle \underset{0<r<1}{min}}\mathrm{SF}3\left[r\right]$ . We can use either $\mathit{SF}3\left[2\right]$ or ${\displaystyle \underset{0<r<1}{min}}\mathrm{SF}3\left[r\right]$ as a shape factor for symmetric distribution to describe the convexity of the $\mathit{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: ${\displaystyle \underset{n>0}{max}} {\displaystyle \underset{0<r<1}{min}}\mathrm{SF}3\left[r\right]=0.942085$ , higher than the Normal distribution family.

Figure 2.

$\mathit{SF}3$ plot of Normal distribution. The horizontal axis is the order $r$ of the power and the vertical axis is $\mathit{SF}3\left[r\right]$ .

When all $\mathit{SF}$ , $\mathit{SF}1$ , and $\mathit{SF}2$ are available, however, we will prefer $\mathit{SF}$ to $\mathit{SF}1$ and $\mathit{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 $\mathit{SF}$ contours are almost lines.

Figure 3.

$\mathit{SF}1$ contour plot of Beta distribution. The horizontal axis is $\alpha$ and the vertical axis is $\beta$ .

Figure 4.

$\mathit{SF}2$ contour plot of Beta distribution. The horizontal axis is $\alpha$ and the vertical axis is $\beta$ .

Figure 5.

$\mathit{SF}$ contour plot of Beta distribution. The horizontal axis is $\alpha$ and the vertical axis is $\beta$ .

2.3 Lower bound of $\mathit{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 $\mathit{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 $\mathit{SF}$	Location of the Min
$\text{FrechetDistribution}\left[\mathrm{\alpha },\mathrm{\beta },\mathrm{\mu }\right]$	2.9555	α → 7.9305
$\text{ExtremeValueDistribution}\left[\mathrm{\alpha },\mathrm{\beta }\right]$	4.15843	any α, β
$\text{MaxStableDistribution}\left[\mu ,\sigma ,\xi \right]$	1.91227	ξ → -1.55970090120176
$\text{InverseGaussianDistribution}\left[\mu ,\lambda ,\theta \right]$	1.5	λ/μ → 0
$\text{SkewNormalDistribution}\left[\mu ,\sigma ,\alpha \right]$	3.90603	α → ∞
$\text{ExpGammaDistribution}\left[\kappa ,\theta ,\mu \right]$	2.25	κ → 0
$\text{BirnbaumSaundersDistribution}\left[\alpha ,\lambda \right]$	1.63481	α → ∞
$\text{MeixnerDistribution}\left[a,b,m,d\right]$	1.5	d → 0,b → ±π

Table 2.

Lower bound of $\mathit{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 $\mathit{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 $\mathit{SF}$ 4.95061 when α is 2.07764, from its proposition 2.3, is thus also not eligible for NATH.

The global lower bound of $\mathit{SF}$ for parametric distribution can be used to filter out those distributions whose values are larger than the losses data $\mathit{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 $\mathit{min}\mathit{SF}=1.0$ , so we need an in-detail study of why it cannot fit NATH.

When matching a $\mathit{\text{BetaDistribution}}\left[\alpha ,\beta \right]$ 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 $\mathit{\text{BetaDistribution}}\left[\alpha ,\beta \right]$ :

$$\alpha \to \frac{\mathrm{cv}-\mathrm{std}-{\mathrm{cv}}^2\mathrm{std}}{{\mathrm{cv}}^3},\beta \to \frac{{\mathrm{cv}}^2-2\text{cvstd}-{\mathrm{cv}}^3\mathrm{std}+{\mathrm{std}}^2+{\mathrm{cv}}^2{\mathrm{std}}^2}{{\mathrm{cv}}^3\mathrm{std}}.$$

Since α > 0, we must have:

$$\mathrm{std}<\frac{\mathrm{cv}}{1+{\mathrm{cv}}^2},$$

or

$$\frac{1-\sqrt{1-4{\mathrm{std}}^2}}{2\mathrm{std}}<\mathrm{cv}<\frac{1+\sqrt{1-4{\mathrm{std}}^2}}{2\mathrm{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 $\mathit{SF}$ takes minimum for a given cv.

The overall observation is that when cv < 1, $\mathit{SF}$ approaches infinity in the middle value of std, and decreases when deviating from it. When cv > 1, $\mathit{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 $\mathit{SF}$ must be attained either at the global extreme where the derivative of $\mathit{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 $\frac{\mathrm{cv}}{1+{\mathrm{cv}}^2}$ with value:

$$\underset{0<\mathit{std}<\frac{\mathit{cv}}{1+{\mathit{cv}}^2}}{\mathit{min}}\mathit{SF}=\frac{1-{\mathit{cv}}^2+{\mathit{cv}}^4}{{\left(-1+{\mathit{cv}}^2\right)}^2},\mathit{\text{when}}\mathit{cv}<0.707107\parallel \mathit{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:

$$\underset{0<\mathit{std}<\frac{\mathit{cv}}{1+{\mathit{cv}}^2}}{\mathit{min}}\mathit{SF}=1.5+\frac{0.75}{{\mathit{cv}}^2},\mathit{\text{when}}0.707107\le \mathit{cv}\le 1.024766.$$ E4

When CV is between 1.024766 and 2.48239, the minimum $\mathit{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 $\mathit{SF}$ value.

3.2 Shape factor range for given skewness

By solving Beta distribution parameters α and β through skewness $\mathit{sk}$ and kurtosis $\mathit{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 $\beta$ parameter. The horizontal axis is the skewness and the vertical axis is the kurtosis.

For a fixed skewness, $\alpha$ is monotonic increasing with respect to kurtosis; on the other hand, $\beta$ 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 $\beta$ parameter and $\alpha$ parameter vs. kurtosis for a given skewness 5.99378.

Solving for that singular point we get the permissible kurtosis upper bound $3+\frac{3}{2}{\mathit{sk}}^2$ , and solve for $\beta =0$ get the permissible kurtosis lower bound $1+{\mathit{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 $\alpha$ and $\beta$ , letting $\beta \to \infty$ , and then eliminating $\alpha$ 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 $\alpha$ by skewness and $\beta$ , substituting the real solution into shape factor, and then take the limit for $\beta \to \infty$ .

All three methods get the same upper bound of $\mathit{SF}=\frac{3}{2}+\frac{3}{{\mathit{sk}}^2}.$

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

$$1+\frac{1}{S^2}\le \mathit{SF}\le 1.5+\frac{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 $\mathit{\text{KumaraswamyDistribution}}\left[\alpha ,\beta \right]$ [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 $\beta \to \infty$ . For a given skewness or a given kurtosis, there exists a maximum allowable α that is arrived when $\beta \to \infty$ . 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 $\beta \to \infty$ , 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 $\alpha$ parameter and the vertical axis is the $\beta$ 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:

$$\underset{\beta \to \infty }{\text{limit}}S=\frac{2\text{Gamma}{\left[\frac{1}{\alpha }\right]}^3-6\alpha \text{Gamma}\left[\frac{1}{\alpha }\right]\text{Gamma}\left[\frac{2}{\alpha }\right]+3{\alpha }^2\text{Gamma}\left[\frac{3}{\alpha }\right]}{{\left(\alpha \left(-\alpha \text{Gamma}{\left[1+\frac{1}{\alpha }\right]}^2+2\text{Gamma}\left[\frac{2}{\alpha }\right]\right)\right)}^{3/2}},$$ E6

$$\underset{\beta \to \infty }{\text{limit}}K=\frac{-3\text{Gamma}\left[\frac{1}{\alpha }\right]\left(\text{Gamma}{\left[\frac{1}{\alpha }\right]}^3-4\alpha \text{Gamma}\left[\frac{1}{\alpha }\right]\text{Gamma}\left[\frac{2}{\alpha }\right]+4{\alpha }^2\text{Gamma}\left[\frac{3}{\alpha }\right]\right)+{\alpha }^4\text{Gamma}\left[\frac{4+\alpha }{\alpha }\right]}{\text{Gamma}{\left[\frac{1}{\alpha }\right]}^4-4\alpha \text{Gamma}{\left[\frac{1}{\alpha }\right]}^2\text{Gamma}\left[\frac{2}{\alpha }\right]+{\alpha }^4\text{Gamma}{\left[\frac{2+\alpha }{\alpha }\right]}^2},$$ E7

$$\underset{\beta \to \infty }{\text{limit}}\frac{K}{S^2}=\frac{{\alpha }^3{\left(-\alpha \text{Gamma}{\left[1+\frac{1}{\alpha }\right]}^2+2\text{Gamma}\left[\frac{2}{\alpha }\right]\right)}^3\left(-3\text{Gamma}\left[\frac{1}{\alpha }\right]\left(\text{Gamma}{\left[\frac{1}{\alpha }\right]}^3-4\alpha \text{Gamma}\left[\frac{1}{\alpha }\right]\text{Gamma}\left[\frac{2}{\alpha }\right]+4{\alpha }^2\text{Gamma}\left[\frac{3}{\alpha }\right]\right)+{\alpha }^4\text{Gamma}\left[\frac{4+\alpha }{\alpha }\right]\right)}{{\left(2\text{Gamma}{\left[\frac{1}{\alpha }\right]}^3-6\alpha \text{Gamma}\left[\frac{1}{\alpha }\right]\text{Gamma}\left[\frac{2}{\alpha }\right]+3{\alpha }^2\text{Gamma}\left[\frac{3}{\alpha }\right]\right)}^2\left(\text{Gamma}{\left[\frac{1}{\alpha }\right]}^4-4\alpha \text{Gamma}{\left[\frac{1}{\alpha }\right]}^2\text{Gamma}\left[\frac{2}{\alpha }\right]+{\alpha }^4\text{Gamma}{\left[\frac{2+\alpha }{\alpha }\right]}^2\right)}.$$ E8

Its plot Figure 16 has two branches, the dividing point is $\alpha \to 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 $\mathrm{\alpha }=1.80143\ast {10}^{-9},\mathrm{\beta }=0.247044$ . The conditional minimum of the shape factor when skewness = 5.99378 is about 1.04753 when $\mathrm{\alpha }={10}^{-10.5},\mathrm{\beta }=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 $\mathit{GB}1$ distribution [15], since $\mathit{\text{KumaraswamyDistribution}}\left[\alpha ,\beta \right]\approx \mathit{\text{GeneralizedBetaDistributionI}}\left[1,\beta ,\alpha ,1\right],$ and the maximum shape factor plot in Figure 18 is lower than that of LogNormalDistribution, the upper bound of $\mathit{GB}1$ . The following section will study a sibling distribution to $\mathit{GB}1$ , fitted as good as $\mathit{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 $\mathit{GB}2$ ([15]), or $\mathit{\text{BetaPrimeDistribution}}\left[p=\alpha ,q=\beta ,\alpha =1,\beta =1\right]$ ([11]): $\mathit{\text{TransformedDistribution}}\left[\frac{x}{1-x},x\approx \mathit{\text{BetaDistribution}}\left[\alpha ,\beta \right]\right]\approx \mathit{\text{BetaPrimeDistribution}}\left[\alpha ,\beta \right]$ . The minimum shape factor of Beta Distribution is 1, but that of the transformed is 1.5:

$$\mathit{\text{NMinimize}}\left[\left\{\frac{3\left(-3+\beta \right)\left(2{\left(-1+\beta \right)}^2+{\alpha }^2\left(5+\beta \right)+\alpha \left(-1+\beta \right)\left(5+\beta \right)\right)}{4\left(-4+\beta \right){\left(-1+2\alpha +\beta \right)}^2},\alpha >0,\beta >4\right\},\left\{\alpha ,\beta \right\}\right]=\left\{1.5000000239052607,\left\{\alpha \to 0.,\beta \to 6.274769836372949\times {10}^7\right\}\right\}.$$

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

$$\mathit{\text{FindMinimum}}\left[\left\{\frac{\mathit{\text{Kurtosis}}\left[\mathit{\text{BetaPrimeDistribution}}\left[p,q,\alpha ,\beta \right]\right]}{\mathit{\text{Skewness}}{\left[\mathit{\text{BetaPrimeDistribution}}\left[p,q,\alpha ,\beta \right]\right]}^2}/.\alpha \to \frac{1}{x}/.q\to 4x+ y/.x\to {10}^z,1.>p>0.,y>1.,-4.<z<-1.\right\},\left\{\left\{p,6.384125235007732\times {10}^{-10}\right\},\left\{y,1.0032844709998097\right\},\left\{z,-2.157370895027263\right\}\right\},\mathit{\text{MaxIterations}}\to 5000\right]=\left\{1.1250258984236121, ,\left\{p\to 2.083731454230264\times {10}^{-8}y\to \mathrm{42:816363091057056} z\to -2.6498169598310573\right\}\right\}.$$

This is the same value as the minimum shape factor for $\mathit{\text{GammaDistribution}}\left[\alpha ,\beta ,\gamma ,\mu \right]$ (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 $\mathit{GB}2$ shape factor change tendency with respect to α, Figure 19, and shape factor change tendency with respect to p, Figure 20.

Figure 19.

$\mathit{GB}2$ distribution shape factor vs. α for fixed p = 10^-3.312 = 0.000487528.

Figure 20.

$\mathit{GB}2$ distribution shape factor vs. p for fixed α = 10^2.6498169598310573 = 446.495.

The $\mathit{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.

$\mathit{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.

$\mathit{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 $\mathit{GB}2$ fitting problem.

5.1 Minimum shape factor for $\mathit{GB}\mathbf{2}$

The skewness and kurtosis matching problem for $\mathit{GB}2$ is very sensitive to the initial parameter ranges given. A study of the minimum shape factor of $\mathit{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 $\mathit{\text{NMinimize}}$ , minimizing the log shape factor instead can overcome this difficulty. The plot is in Figure 22.

Figure 22.

The numerical minimum $\mathit{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 $\mathit{GB}2$ is a very smooth curve. The fitted formula of $\mathit{GB}2$ $\mathit{min}\mathit{SF}$ for given $p$ by Mathematica’s machine learning function $\mathit{\text{FindFormula}}$ is Eq. (9).

$$min\frac{K}{S^2}=1.1593871374775397+1.4702458297305288\ast {0.5148499158800361}^{\frac{1}{p^{0.3215433282777008}}}$$ E9

As a test, for NATH the log shape factor is $\mathit{Log}\left[1.83408\right]=0.60654412$ , the solution of Eq. (9) for $p$ with NATH SF is $p=0.608342$ ; the minimum log shape factor of $\mathit{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 $\mathit{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=\frac{\lambda }{\alpha },q=\frac{4+\nu }{\alpha }$ , the minimization of shape factor for $\mathit{GB}2$ , for given $\lambda =\mathit{p\alpha }$ , is easier to perform. The plot is included in Figure 23 as well.

Figure 23.

The numerical minimum $\mathit{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,\alpha >0.5,\mathit{p\alpha }<0.5.$ This is confirmed by $\mathit{GB}2$ fit practice. The best fit by $\mathit{GB}2$ for NATH is at $w=-0.329075005,p=0.468732$ , with about 5% error from input TVaR. The discontinuity of fitted $\mathit{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: $\alpha \leftrightarrow \alpha ,\beta \leftrightarrow \theta ,\gamma \leftrightarrow \beta ,\mu \leftrightarrow a.$

For generalized gamma distribution $\mathit{\text{GammaDistribution}}\left[\alpha ,\beta ,\gamma ,\mu \right]$ , the shape factor depends only on $\alpha$ and $\gamma$ . It seems the smaller the $\alpha$ , and the bigger the $\gamma$ , the smaller the $\frac{K}{S^2}$ . When $\alpha =3.318512677036329\times {10}^{-12},\gamma =8811.572418686921,\frac{K}{S^2}=1.125$ , close to the global minimum 1 of K/S^2.

So there arises the question: the generalized gamma and $\mathit{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 $\mathit{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], ${\gamma }_4>1.738{\gamma }_3^2-0.3544{\gamma }_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 $\mathit{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), $\mathit{\text{HyperbolicDistribution}}\left[\lambda ,\alpha ,\beta ,\delta ,\mu \right]$ , 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 $\mathit{\text{BesselK}}\left[\lambda ,a\right]$ , or $K_{\lambda }\left(a\right)$ 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.

$$\underset{\alpha ,\beta ,\delta ,\mu }{\mathit{min}}\mathit{SF}=1.5+0.75\lambda ,\mathit{\text{when}}\lambda \ge 0,$$ E10

$$\underset{\alpha ,\beta ,\delta ,\mu }{\mathit{min}}\mathit{SF}=1+\frac{1}{2+\lambda },\mathit{\text{when}}-0.6\le \lambda \le 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 $\mathit{\text{FindFormula}}$ .

$$\underset{\alpha ,\beta ,\delta ,\mu }{\mathit{min}}\mathit{SF}=1.1130471668735116-1.6512030619809768\lambda -1.6137376956833365{\lambda }^2-1.1485038172210114{\lambda }^3-0.5421785615853132{\lambda }^4-0.17094578834265223{\lambda }^5-0.03603744794749387{\lambda }^6-0.005000441043297472{\lambda }^7-0.0004372189547557593{\lambda }^8-0.000021791071048963054{\lambda }^9-4.711954312790356\times {10}^{-7}{\lambda }^{10}$$ E12

$$\underset{\alpha ,\beta ,\delta ,\mu }{\mathit{min}}\mathit{SF}=2.2104215691249425-0.6522131009473879\ast {1.6355318649123258}^{\lambda }+\frac{0.018965779149540653}{{\lambda }^3}-0.1051542360603726\lambda$$ 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 $\mathit{\text{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 $\mathit{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. 5–8, 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 $\mathit{GB}2$ case. The optimization function $\mathit{\text{NMinimize}}$ and $\mathit{\text{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 C²oDE 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 $B1$ 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.

Acknowledgments

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.