# On the Usefulness of the Logarithmic Skew Normal Distribution for Describing Claims Size Data

Authors: E. Gómez–Déniz; E. Calderín-Ojeda
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, the three-parameter skew lognormal distribution is proposed to model actuarial data concerning losses. This distribution yields a satisfactory fit to empirical data in the whole range of the empirical distribution as compared to other distributions used in the actuarial statistics literature. To the best of our knowledge, this distribution has not been used in insurance context and it might be suitable for computing reinsurance premiums in situations where the right tail of the empirical distribution plays an important role. Furthermore, a regression model can be simply derived to explain the response variable as a function of a set of explanatory variables.


---

## Body

## 1. Introduction

The modelling of large claims is a topic of relevant importance in general insurance and reinsurance, particularly in the field of mathematical risk theory, see, for instance, McNeil [1]; Beirlant and Teugels [2]; Beirlant et al. [3]; and Embrechts et al. [4]. For a comprehensive study about reinsurance, see the recent book of Albrecher et al. [5]. It is our interest to find simple statistical distributions appropriate for modelling both, smaller and medium-size losses with a high frequency and large losses with a low frequency. This is an issue of particular interest in the context of reinsurance and premium calculation principles. In this sense, the classical Pareto distribution [6–8] has been traditionally considered as a suitable claims’ size distribution in relation to rating problems. Concerning this, the single parameter Pareto distribution not only has nice statistical properties but also provides a good description of the random behaviour of large losses (e.g., the right tail of the distribution). Particularly, when calculating deductibles and excess-of-loss levels of reinsurance, the simple Pareto distribution has been proved to be convenient (see, for example, [9, 10]). As an alternative to the classical Pareto distribution, other models have been recently introduced in the actuarial literature by Sarabia et al. [11–13] and Ghitany et al. [14].On the contrary, there exist many situations where the empirical data show slight or marked asymmetry. This is frequently the case, for example, with actuarial and financial data that also have heavy tails reflecting the existence of extreme values. These two features imply that the data cannot be adequately modelled by the Gaussian or normal distribution.Letg and G, respectively, be the probability density function (pdf) and the cumulative distribution function (cdf) of a symmetric distribution. A random variable Z is said to have a skew distribution if its pdf is given by(1)fZz=2gzGγz,−∞<z<∞,γ∈ℝ.This family of distributions has been widely studied as an extension of the normal distribution via a shape parameter,γ, that accounts for the skewness of the model. When g⋅ and G⋅ are replaced in (1) by ϕz and Φz, i.e., the standard normal density and distribution function, respectively, the resulting model is called the skew normal distribution.In this paper, special attention is paid to the generalized skew normal density provided in Henze [15] and also studied by Arnold and Beaver [16]. Its pdf is given by(2)fZz=ϕzΦγ0+γ1zΦγ0/1+γ12,−∞<z<∞,γ0,γ1∈ℝ.For multivariate extensions, see for instance Azzalini and Valle [17], Azzalini and Capitanio [18], and Arnold and Beaver [16]. For an exhaustive and comprehensive study of the skew-normal distribution, see the recent book of Azzalini [19].This paper is organized as follows. In Section2, the proposed distribution is studied and some interesting properties are given. Numerical applications are provided in Section 3, and conclusions are drawn in Section 4.


## 2. Modelling the Size of Losses

The basic skew lognormal distribution has been studied by Lin and Stoyanov [20] (see also [19]; Chap. 2, [21–23]). Nevertheless, in this paper, we will pay special attention to the distribution arising from exponentiation of (2) in the following sense: if X is a random variable with density (2), X∼fZ, we consider a new random variable Y∼expX. Also, if taking γ0=γ1=λ and using a linear transformation, we allow for more flexible location and scale parameters in our model. Thus, we get the pdf given by(3)fμ,σ,λx=ϕrμ,σxσxΦ1+rμ,σxλΦλ0,x>0,where(4)rμ,σx=logx−μσ,λ0=λ1+λ2.Here,λ∈ℝ, μ∈ℝ, and σ>0. Observe that when λ=0 expression (3) reduces to the classical lognormal distribution. It can be seen that the parameter λ regulates the shape of the distribution. Furthermore, Lin and Stoyanov [20] established that the distribution has heavy tails and therefore it is a suitable distribution to be incorporated to the wide catalogue of heavy tails (see [5]). Additionally Lin and Stoyanov [20] provided a stochastic representation of this distribution and determined that the distribution can be obtained as the product of two independent random variables, one of which is lognormal and the other one a logarithmic half-normal. Simple calculations show that density (3) can be also written as(5)fμ,σ,λx=exp1/21+2rμ,σxΦλ0fZz,where fZz is expression (1) with z=1+rμ,σx, γ=λ, gz=ϕz, and Gz=Φz.It is straightforward to show that the distribution is unimodal with the mode satisfying the equation:(6)λϕ1+rμ,σzΦ1+rμ,σz=rμ,σz+σ.The mean value and the second order moment of the random variableX with density (3), X∼LSNμ,σ,λ, are given by(7)EX=Φ1+σλ0Φλ0expμ+σ22,EX2=Φ1+2σλ0Φλ0exp2μ+σ2.LetFμ,σ,λx=PrX≤x denote the cdf of X∼LSNμ,σ,λ. By applying directly result B.21 in Azzalini [19], it is not difficult to observe that this cdf is given by(8)Fμ,σ,λx=Φrμ,σx+1Φλ0Trμ,σx,λ0rμ,σx+Tλ0,rμ,σxλ0−Trμ,σx,1+rμ,σxλrμ,σx−Tλ0,λ2+1+λ2rμ,σxλ,which is satisfied whenever λ≠0 and x>0. Here, Tx,a represents Owen’s function (see [24]) given by(9)Tx,a=12π∫0a11+t2exp−12x21+t2dt,a∈ℝ+.Ifλ=0, we deal with the standard normal distribution, and Tx,0=0. The latter function can be easily computed by using the Mathematica software package.The Acceptance-Rejection method of simulation can be used to generate random variates from (3). We begin by simulating a value from a lognormal distribution with parameters μ and σ. Then, having chosen an alternative random variable that follows a lognormal probability distribution to simulate from, we define a constant c in the following way:(10)c=maxxfμ,σ,λxgx,where fμ,σ,λx is given by (3) and gx is the standard lognormal density. The algorithm for simulating a random variate from the LSN distribution is as follows:(1)
Generate a random variate from the lognormal distribution,x1.(2)
Generate a random numberu1.(3)
Ifu1≤Φ1+rμ,σx1λ, then set the simulated value from (3) equal to x1. Otherwise, return to step 1.By writingμ as(11)μ=−σ22+logθΦλ0Φ1+σλ0,pdf (3) has mean value equal to θ>0 and variance provided by(12)varX=ψ1θψ2σ,λ,where(13)ψ1θ=θ2,ψ2σ,λ=Φ1+2σλ0Φλ02exp−σ2Φ1+σλ02−1.Therefore, a regression model can be derived from the LSN distribution. In this model,θ is the mean of the response variable and ψ2σ,λ can be interpreted as a precision parameter in a way such that, for fixed θ, the larger the value of ψ2σ,λ, the larger the variance of X.Let us now consider that the random variableXi is related to a vector of k covariates yi=y1i,…,yki′ associated with the ith observation, where i=1,2…,n. This vector assigns a weight of observable features is related to θi through a log link function and β=β1,…,βk′ is a vector of regressors with βj∈ℝ for j=1,…,k. The regression model takes the form(14)Xi∼LSNθi,σ,λ,logθi=yi′β.Observe that log link ensures thatθi falls within the interval 0,∞.On the contrary, sometimes it is desirable to deal with a lower value in the range ofX by shifting this random variable, say X−α with α>0, and therefore a shifted version of the previous distribution is obtained. In this case, it is convenient to work with the pdf given by(15)fμ,σ,αx=ϕr˜μ,σ,αxσx−αΦ1+r˜μ,σ,αxλΦλ0,x>α,where(16)r˜μ,σ,αx=logx−α−μσ.Figure1 shows pdf (15) as compared to the Pareto and PAT (Pareto ArcTan distribution) (this generalization of the Pareto distribution was proposed by Gómez-Déniz and Calderín-Ojeda [13]), all of them with the same mean, say (22), and a threshold of α=1. It is discernible that the tail of the LSN density seems larger than those ones of the other models.Figure 1
Graph of the Pareto (thin line), PAT (dashed line), and LSN (thick line) with the same mean given by (22).By taking(17)μ=−σ22+logθ−αΦλ0Φ1+σλ0,it is guaranteed that (15) has mean equal to θ>α. In this case, for the associated regression model, we will use the link function:(18)θi=α+expyi′β,which guarantees that θi falls within the interval α,∞.It is already known (see Lemma 2 in [25]) that(19)limz⟶∞−lnΦ−zz2=12.Then, it is straightforward to see that ifλ<0, we have(20)limx⟶∞−lnΦ1+rμ,σ,αxλ1+rμ,σ,αx2=12.This can be used to easily obtain asymptotic values of the pdf given in (15). In fact, if λ<0, we have the following approximation:(21)fx≈ϕrμ,σ,αxσx−αΦλ0exp−121+rμ,σ,αxλ2,x>α,x⟶∞.


## 3. Empirical Illustrations

In this section, the versatility of the proposed skew lognormal model, as compared with the Pareto distribution and the Pareto ArcTan Distribution [13], is tested using three datasets. The first dataset is the danishuni that can be found in the R package CASdatasets collected at Copenhagen Reinsurance and comprises 2167 fire losses over the period 1980 to 1990, adjusted for inflation to reflect 1985 values and are expressed in millions of Danish Krone. The second dataset is norfire comprises 9181 fire losses over the period 1972 to 1992 from an unknown Norwegian company. A priority of 500 thousands of Norwegian Krone (NKR) (if this amount is exceeded, the reinsurer becomes liable to pay) was applied to obtain this dataset. The third dataset considers hospital costs in the state of Wisconsin. It can be found in the personal web page of Professor E. Frees. It examines the impact of several predictors on hospital charges obtained from the Office of Health Care Information, Wisconsin’s Department of Health and Human Services, in the year 1989. The response variable is the logarithm of total hospital charges per number of discharges, and the explanatory variables considered are the health service area, hospital discharge costs, the type of payer, size of the area population, number of hospital beds, and average income within the area. For comparison reasons some alternative distributions have been considered in this paper, their pdf’s are provided in Appendix A and B. Some descriptive statistics of the three variables of interest are shown in Table 1.Table 1
Descriptive statistics of the dependent variable.
DataMeanStand. dev.MinimumMaximumDanish56.339534.1611.0014239Norwegian2217.217759.97500.000465365Hospital cost (in log scale)7.945370.751676.105049.64179Parameter estimation for all the models considered in this paper has been carried out by the method of maximum likelihood using Mathematica v.12.0® and also confirmed by using WinRATS v.7.0. Both moments and the maximum likelihood methods are suitable to estimate the vector of parameters of the distribution via sample observations, as shown in AppendixA and B. Codes are available from the authors upon request. For details about software, see Ruskeepaa [26] and Brooks [27].Since closed expressions are not available for the maximum likelihood estimates and the computation of the global maximum of the logarithm function of the likelihood is not guaranteed, thus it is advisable to use several seed points as the starting value. It is also prudent to use different optimization methods (Newton–Raphson, Broyden–Fletcher–Goldfarb–Sanno, BGGS) (this is an iterative method for solving unconstrained nonlinear optimization problems which can be seen in Broyden [28, 29] and Shanno [30]) that ensure that the same solution is obtained from any of these methods. The standard errors of the estimates can be calculated by inverting the Hessian matrix. In this sense, both Mathematica and WinRats have at least two methods to reach it. The first is to retrieve them from the Cholesky factors (this package is available on the web upon request). The second, faster, is to obtain them by finite differentiation. Furthermore, WinRats package also offers the possibility to directly compute the maximum of the log-likelihood providing the elements of the Fisher information matrix. In fact, for the examples considered in this section, these two packages were used to quickly compute the maximum likelihood estimates. Commands for fitting the skew normal and log-skew normal distributions are also available in stata (see [31]).Parameter estimates and their standard errors (in brackets), negative of the maximum of the log likelihood function, and AIC for the three datasets considered are shown in Tables2 and 3. Model assessment is derived through the following information criteria. Negative log likelihood (NLL) is calculated by taking the negative of the value of the log-likelihood evaluated at the ML estimates; Akaike information criterion (AIC) is computed as twice the NLL and evaluated at the ML estimates, plus twice the number of estimated parameters; Consistent Akaike Information Criteria (CAIC), a corrected version of the AIC is proposed by Bozdogan [32] to overcome the tendency of the AIC overestimating the complexity of the underlying model as it lacks the asymptotic property of consistency. In order to calculate the CAIC, a correction factor based on the sample size is used to compensate for the overestimating nature of AIC. The CAIC is defined as twice the NLL plus k1+lnn, where k is the number of free parameters and n refers to the sample size. Note that a model with a lower statistics value is preferred to one with a higher value. All these results are shown in Table 2. Furthermore, we also include the Kolmogorov–Smirnov test (KS) and the Anderson–Darling test (AD) to express the fit of the model to the data in terms of the distribution function. For these statistics, smaller values indicate a better fit of the model to the data. Note that they not only provide a way to measure the fit in terms of distribution functions but also allow us to perform hypothesis testing for model selection purposes. The p value of the test statistics, computed using the Monte Carlo method by using 1000 simulations, is shown in brackets. An extremely small p value may lead to a confident rejection of the null hypothesis that the data comes from the proposed model. It can be seen that LSN distribution is not rejected at the usual significance levels for both tests.Table 2
Parameter estimates and their standard errors (in brackets), negative of the maximum of the log likelihood function, and AIC without covariates.
λ^θ^σ^α^NLLAICCAICKSADDanish data−1.3243.1212.0750.9933361.4866730.9706757.7000.0290.001(0.308)(0.097)(0.222)(0.004)(0.319)(0.136)Norwegian data657.2181894.5001.239318.83620966.83541941.70041969.1000.0310.001(18.987)(43.876)(0.027)(9.600)(0.148)(0.049)Table 3
Parameter estimates and standard errors (in brackets) of the homogeneous model together with some measures of model selection for different distributions considered for Wisconsin hospital costs’ dataset.
Modelλ^θ^σ^α^NLLAICCAICKSADPareto3.6206.000937.0191876.0401881.3000.3450.224(0.158)(0.000)(0.000)PAT15.900647.2735.307600.4071206.3301222.1200.8461.355(0.628)(10.897)(0.083)(0.000)(0.000)LSN−3.1327.9410.2772.601590.1341188.2701209.3300.0510.004(0.338)(0.031)(0.007)(0.123)(0.492)(0.260)Graphs of histogram of the data and superimposed fitted densities are given in Figure2. As can be seen, the proposed distribution provides a good fit for the empirical data.Figure 2
Smooth histogram of the empirical Danish, Norwegian insurance data, and hospital cost (in log scale) as compared to fitted models. Pareto (dashed line), PAT (thin line), and LSN (thick line). (a) Danish insurance loss data. (b) Norwegian insurance loss data. (c) Hospital costs’ data.
(a)(b)(c)In Table4, the LSN distribution is fitted to the Wisconsin hospital costs’ dataset. From left to right the parameter estimates, standard errors (SE), and the corresponding p values calculated based on the t-Wald statistic are displayed. Besides, AIC and BIC values for each model are provided in the last two rows of the table. The estimates of the regressors associated with the explanatory variables, number of beds and income, are not significant at the usual nominal level.Table 4
Parameter estimates for Wisconsin hospital costs dataset including covariates.
VariableEstimateSEt-Stat.p valueHospital discharge costs0.0150.0026.370.000Health service area0.0020.0011.740.080Payer0.0120.0052.430.014Size area population0.0150.0053.100.001Number of beds−0.0100.008−1.200.228Income0.0010.0050.310.756Constant1.6070.008200.320.000α0.5630.2692.090.036λ−2.9740.361−8.210.000σ0.1860.01413.030.000NLL=564.619AIC=1149.240CAIC=1201.890Finally, in order to choose a model that provides an acceptable description of the loss process for the Danish and Norwegian datasets, we should verify that the population limited expected values, computed numerically by(22)Lx=EminX,x=∫0xydFy+xF¯x,are close to the empirical ones. As is well known, (22) is the expected amount per claim retained by the insured on a policy with a fixed amount deductible of x. In this case, the empirical limited expected value function was calculated based on Enx=1/n∑i=1nminxi,x. Obviously, when x tends to infinity, Lx and Enx approach EX and the sample mean, respectively.Table5 below exhibits the limited expected value for different values of the policy limit x considered for the hospital costs’ dataset. It is observed that the values obtained from the LSN distribution adheres closely to the observed empirical limited expected values obtained from the Pareto and PAT distributions. Similarly, in Figure 3, empirical and fitted limited expected value function for this dataset and also for the Danish and Norwegian data are shown. As can be observed, the LSN distribution stays closer to the empirical limited expected values for the three datasets considered.Table 5
:Limited expected value for the different distributions considered and different values of the fixed amount deductiblex for hospital costs dataset.
Limited expected valueDeductibleEmpiricalParetoPATLSN6.006.000006.000005.998535.998106.156.149786.143466.147406.146616.306.296776.274806.295686.294196.456.440476.395296.443086.440446.606.580926.506056.589276.584826.756.719056.608066.733736.726666.906.854586.702186.875806.865177.056.988836.789187.014566.999407.207.119336.869737.148817.128277.357.244756.944427.277027.250597.507.364897.013807.397307.365067.657.473667.078337.507537.470397.807.571797.138437.605647.565337.957.658447.194507.690077.648828.107.733727.246867.760287.720198.257.791207.295827.816947.779258.407.832137.341677.861627.826408.557.862657.384657.896347.862608.707.885677.424987.923117.889288.857.904377.462877.943697.908139.007.919497.498507.959517.920889.157.931577.532057.971707.929149.307.940657.563677.981117.934279.457.944387.593497.988397.937319.607.945277.621657.994067.939059.757.945377.648267.998487.940019.907.945377.673428.001947.9405110.007.945377.689458.003837.94070Figure 3
Empirical (thick line) and fitted limited expected values: Pareto distribution (dashed line), PAT distribution (dotted line), and LSN (thin line) for the (a) Danish data, (b) Norwegian data, and (c) hospital cost data.
(a)(b)(c)


### 3.1. Out-of-Sample Validation of LSN Distribution

We are interested now in showing the power of the LSN distribution to predict the number of claims out-of-sample via quantile-quantile plots (QQ-plots) of randomized quantile residuals. For this reason, the datasetdanishuni is randomly partitioned into two disjoint subsets of different sizes. The first subset (dataset A) is used for fitting the models, whereas the second subset (dataset B) is used for graphing the quantile-quantile plots of the randomized quantile residuals to check for normality. The residuals for the out-of-sample datasets are exhibited in Figure 4 for different in-sample sizes, i.e., size of dataset A. For comparison purposes, we have we have plotted the residuals for the whole dataset (top-left graph) with sample size 2167. A perfect alignment with the 45° line implies the residuals are normally distributed. It is observable that the residuals for LSN distribution adheres reasonably well to the line throughout the residual distribution. However, for all different in-sample sizes considered, the model tends to overestimate the lower quantiles and underestimate the upper quantiles. In general, the predictive power of the model is acceptable.Figure 4
QQ-plots of randomized quantile residuals for the LSN distribution for the whole datasetdanishuni and different in-sample sizes. (a) Whole dataset size-2167. (b) In-sample size-1500. (c) In-sample size-1000. (d) In-sample size-750.
(a)(b)(c)(d)


## 3.1. Out-of-Sample Validation of LSN Distribution

We are interested now in showing the power of the LSN distribution to predict the number of claims out-of-sample via quantile-quantile plots (QQ-plots) of randomized quantile residuals. For this reason, the datasetdanishuni is randomly partitioned into two disjoint subsets of different sizes. The first subset (dataset A) is used for fitting the models, whereas the second subset (dataset B) is used for graphing the quantile-quantile plots of the randomized quantile residuals to check for normality. The residuals for the out-of-sample datasets are exhibited in Figure 4 for different in-sample sizes, i.e., size of dataset A. For comparison purposes, we have we have plotted the residuals for the whole dataset (top-left graph) with sample size 2167. A perfect alignment with the 45° line implies the residuals are normally distributed. It is observable that the residuals for LSN distribution adheres reasonably well to the line throughout the residual distribution. However, for all different in-sample sizes considered, the model tends to overestimate the lower quantiles and underestimate the upper quantiles. In general, the predictive power of the model is acceptable.Figure 4
QQ-plots of randomized quantile residuals for the LSN distribution for the whole datasetdanishuni and different in-sample sizes. (a) Whole dataset size-2167. (b) In-sample size-1500. (c) In-sample size-1000. (d) In-sample size-750.
(a)(b)(c)(d)


## 4. Conclusions and Extensions

In this paper, the use of the skew lognormal distribution has been proposed as a suitable model for describing claims’ size empirical data. To the best of our knowledge, this probabilistic family has not been previously considered in the actuarial literature. The advantages of this distribution are twofold. On the one hand, it includes heavy tails which makes it an interesting model to describe severity data and; on the other hand, the distribution can be rewritten to allow for the incorporation of predictors to explain a response variable.Finally, it is worthy to point out that the distribution introduced in this paper could be extended to other areas of risk theory. For example, the LSN distribution could be employed as the secondary distribution (the distribution of the claims size) in the compound collective risk model. In particular, this distribution could be used as an approximation of the distribution of the total claims’ size instead of the traditionally considered normal distribution. Furthermore, calculation of reinsurance premiums based on this distribution is an issue that deserves to be studied.



---

*Source: 1420618-2020-09-18.xml*