Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 1420618, 9 pages
https://doi.org/10.1155/2020/1420618

Research Article

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

E. Gómez–Déniz¹ and E. Calderín-Ojeda²

¹Department of Quantitative Methods in Economics and Tides Institute, University of Las Palmas de Gran Canaria, Las Palmas de Gran Canaria 35017, Spain
²Centre for Actuarial Studies, Department of Economics, The University of Melbourne, Melbourne, Australia

Correspondence should be addressed to E. Calderín-Ojeda; ecalderin@unimelb.edu.au

Received 24 May 2020; Accepted 17 July 2020; Published 18 September 2020

Academic Editor: Łukasz Jankowski

Copyright © 2020 E. Gómez–Déniz and E. Calderín-Ojeda. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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.

Let and , respectively, be the probability density function (pdf) and the cumulative distribution function (cdf) of a symmetric distribution. A random variable is said to have a skew distribution if its pdf is given by

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 and are replaced in (1) by and , 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

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 Section 2, 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 is a random variable with density (2), , we consider a new random variable . Also, if taking and using a linear transformation, we allow for more flexible location and scale parameters in our model. Thus, we get the pdf given by

where

Here, , , and . Observe that when 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

where is expression (1) with , , , and .

It is straightforward to show that the distribution is unimodal with the mode satisfying the equation:

The mean value and the second order moment of the random variable with density (3), , are given by

Let denote the cdf of . By applying directly result B.21 in Azzalini [19], it is not difficult to observe that this cdf is given by

which is satisfied whenever and . Here, represents Owen’s function (see [24]) given by

If , we deal with the standard normal distribution, and . 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 in the following way:

where is given by (3) and 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, .(2)Generate a random number .(3)If , then set the simulated value from (3) equal to . Otherwise, return to step 1.

By writing as

pdf (3) has mean value equal to and variance provided by

where

Therefore, a regression model can be derived from the LSN distribution. In this model, is the mean of the response variable and can be interpreted as a precision parameter in a way such that, for fixed , the larger the value of , the larger the variance of .

Let us now consider that the random variable is related to a vector of covariates associated with the th observation, where . This vector assigns a weight of observable features is related to through a log link function and is a vector of regressors with for . The regression model takes the form

Observe that log link ensures that falls within the interval .

On the contrary, sometimes it is desirable to deal with a lower value in the range of by shifting this random variable, say with , and therefore a shifted version of the previous distribution is obtained. In this case, it is convenient to work with the pdf given by

where

Figure 1 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 . 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

it is guaranteed that (15) has mean equal to . In this case, for the associated regression model, we will use the link function:

which guarantees that falls within the interval .

It is already known (see Lemma 2 in [25]) that

Then, it is straightforward to see that if , we have

This can be used to easily obtain asymptotic values of the pdf given in (15). In fact, if , we have the following approximation:

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.

Data
Mean
Stand. dev.
Minimum
Maximum

Danish
56.339
534.161
1.00
14239

Norwegian
2217.21
7759.97
500.000
465365

Hospital cost (in log scale)
7.94537
0.75167
6.10504
9.64179

Parameter 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 Appendix A 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 Tables 2 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 , where is the number of free parameters and 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 value of the test statistics, computed using the Monte Carlo method by using 1000 simulations, is shown in brackets. An extremely small 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.

NLL
AIC
CAIC
KS
AD

Danish data

−1.324
3.121
2.075
0.993
3361.486
6730.970
6757.700
0.029
0.001

(0.308)
(0.097)
(0.222)
(0.004)

(0.319)
(0.136)

Norwegian data

657.218
1894.500
1.239
318.836
20966.835
41941.700
41969.100
0.031
0.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

NLL
AIC
CAIC
KS
AD

Pareto

3.620
6.000
937.019
1876.040
1881.300
0.345
0.224

(0.158)

(0.000)
(0.000)

PAT

15.900
647.273
5.307
600.407
1206.330
1222.120
0.846
1.355

(0.628)
(10.897)
(0.083)

(0.000)
(0.000)

LSN
−3.132
7.941
0.277
2.601
590.134
1188.270
1209.330
0.051
0.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 Figure 2. As can be seen, the proposed distribution provides a good fit for the empirical data.

(a)

(b)

(c)

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.

In Table 4, the LSN distribution is fitted to the Wisconsin hospital costs’ dataset. From left to right the parameter estimates, standard errors (SE), and the corresponding values calculated based on the -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.

Variable
Estimate
SE
-Stat.
value

Hospital discharge costs
0.015
0.002
6.37
0.000

Health service area
0.002
0.001
1.74
0.080

Payer
0.012
0.005
2.43
0.014

Size area population
0.015
0.005
3.10
0.001

Number of beds
−0.010
0.008
−1.20
0.228

Income
0.001
0.005
0.31
0.756

Constant
1.607
0.008
200.32
0.000

0.563
0.269
2.09
0.036

−2.974
0.361
−8.21
0.000

0.186
0.014
13.03
0.000

Finally, 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

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 . In this case, the empirical limited expected value function was calculated based on . Obviously, when tends to infinity, and approach and the sample mean, respectively.

Table 5 below exhibits the limited expected value for different values of the policy limit 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 deductible for hospital costs dataset.

Limited expected value

Deductible
Empirical
Pareto
PAT
LSN

6.00
6.00000
6.00000
5.99853
5.99810

6.15
6.14978
6.14346
6.14740
6.14661

6.30
6.29677
6.27480
6.29568
6.29419

6.45
6.44047
6.39529
6.44308
6.44044

6.60
6.58092
6.50605
6.58927
6.58482

6.75
6.71905
6.60806
6.73373
6.72666

6.90
6.85458
6.70218
6.87580
6.86517

7.05
6.98883
6.78918
7.01456
6.99940

7.20
7.11933
6.86973
7.14881
7.12827

7.35
7.24475
6.94442
7.27702
7.25059

7.50
7.36489
7.01380
7.39730
7.36506

7.65
7.47366
7.07833
7.50753
7.47039

7.80
7.57179
7.13843
7.60564
7.56533

7.95
7.65844
7.19450
7.69007
7.64882

8.10
7.73372
7.24686
7.76028
7.72019

8.25
7.79120
7.29582
7.81694
7.77925

8.40
7.83213
7.34167
7.86162
7.82640

8.55
7.86265
7.38465
7.89634
7.86260

8.70
7.88567
7.42498
7.92311
7.88928

8.85
7.90437
7.46287
7.94369
7.90813

9.00
7.91949
7.49850
7.95951
7.92088

9.15
7.93157
7.53205
7.97170
7.92914

9.30
7.94065
7.56367
7.98111
7.93427

9.45
7.94438
7.59349
7.98839
7.93731

9.60
7.94527
7.62165
7.99406
7.93905

9.75
7.94537
7.64826
7.99848
7.94001

9.90
7.94537
7.67342
8.00194
7.94051

10.00
7.94537
7.68945
8.00383
7.94070

(a)

(b)

(c)

Figure 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.

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 dataset danishuni 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 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.

(a)

(b)

(c)

(d)

Figure 4: QQ-plots of randomized quantile residuals for the LSN distribution for the whole dataset danishuni 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.

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.

Appendix

A. PDF of the Alternative Distributions

The pdf of the classical Pareto and PAT distributions used in this work are given by

respectively.

B. Normal Equations

Let us firstly consider the case of the model without covariates and let be a sample obtained from distribution (3). By denoting as the vector of parameters to be estimated, the log-likelihood is proportional to

The required normal equations to compute the maximum likelihood estimates are

Data Availability

The csv data used to support the ﬁndings of this study can be found in the R repository CASdatasets and also at https://instruction.bus.wisc.edu/jfrees/jfreesbooks/Regression%20Modeling/BookWebDec2010/data.html.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

EGD would like to thank the Ministerio de Economía y Competitividad (project ECO2017-85577-P) for partial support of this work (Ministerio de Economía y Competitividad, Spain).

References

A. McNeil, “Estimating the tails of loss severity distribution using extreme value theory,” ASTIN Bulletin, vol. 27, no. 1, pp. 117–137, 1988.

J. Beirlant and J. L. Teugels, “Modeling large claims in non-life insurance,” Insurance: Mathematics and Economics, vol. 11, no. 1, pp. 17–29, 1992.

J. Beirlant, J. Teugels, and P. Vynckier, Practical Analysis of Extreme Values, Leuven University Press, Leuven, Belgium, 1996.

P. Embrechts, S. I. Resnick, and G. Samorodnitsky, “Extreme value theory as a risk management tool,” North American Actuarial Journal, vol. 3, no. 2, pp. 30–41, 1999.

H. Albrecher, J. Beirlant, and J. Teugels, Reinsurance: Actuarial and Statistical Aspects, Wiley, Hoboken, NJ, USA, 2017.

B. Arnold, Pareto Distributions, International Cooperative Publishing House, Silver Spring, MD, USA, 1983.

V. Brazauskas and R. Serfling, “Favorable estimators for fitting Pareto models: a study using goodness-of-fit measures with actual data,” ASTIN Bulletin, vol. 33, no. 2, pp. 365–381, 2003.

M. Rytgaard, “Estimation in the Pareto distribution,” ASTIN Bulletin, vol. 20, no. 2, pp. 201–216, 1990.

A. V. Boyd, “Fitting the truncated Pareto distribution to loss distributions,” Journal of the Staple Inn Actuarial Society, vol. 31, pp. 151–158, 1988.

A. J. Mata, “Pricing excess of loss reinsurance with reinstatements,” ASTIN Bulletin, vol. 30, no. 2, pp. 349–368, 2000.

J. M. Sarabia and F. Prieto, “The Pareto-positive stable distribution: a new descriptive model for city size data,” Physica A: Statistical Mechanics and Its Applications, vol. 388, no. 19, pp. 4179–4191, 2009.

E. Gómez-Déniz and E. Calderín-Ojeda, “A suitable alternative to the Pareto distribution,” Hacettepe Journal of Mathematics and Statistics, vol. 43, no. 5, pp. 843–860, 2014.

E. Gómez-Déniz and E. Calderín-Ojeda, “Modelling insurance data with the Pareto arctan distribution,” ASTIN Bulletin, vol. 45, no. 3, pp. 639–660, 2015.

M. Ghitany, E. Gómez-Déniz, and S. Nadarajah, “A new generalization of the Pareto distribution and its application to insurance data,” Journal of Risk and Financial Management (Special Issue: Extreme Values and Financial Risk), vol. 11, no. 10, pp. 1–14, 2008.

N. Henze, “A probabilistic representation of the skew-normal distribution,” Scandinavian Journal of Statistics, vol. 13, no. 4, pp. 271–275, 1986.

B. Arnold and R. Beaver, “Skewed multivariate models related to hidden truncation and/or selective reporting (with discussion),” Test, vol. 11, no. 1, pp. 1–54, 2002.

A. Azzalini and A. Valle, “The multivariate skew-normal distribution,” Biometrika, vol. 83, no. 4, pp. 715–726, 1996.

A. Azzalini and A. Capitanio, “Statistical applications of the multivariate skew normal distribution,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 61, no. 3, pp. 579–602, 1999.

A. Azzalini, The Skew-Normal and Related Families, Cambridge University Press, Cambridge, UK, 2013.

G. D. Lin and J. Stoyanov, “The logarithmic skew-normal distributions are moment-indeterminate,” Journal of Applied Probability, vol. 46, no. 3, pp. 909–916, 2009.

A. Azzalini, D. Cappello, and S. Kotz, “Log-skew-normal and log-skew-t distributions as model for family income data,” Journal of Income Distribution, vol. 11, pp. 12–20, 2003.

H. Bolfarine, H. W. Gómez, and L. I. Rivas, “The log-bimodal-skew-normal model. A geochemical application,” Journal of Chemometrics, vol. 25, no. 6, pp. 329–332, 2011.

O. Venegas, H. Bolfarine, D. I. Gallardo, A. Vergara-Fernández, and H. W. Gómez, “A note on the log-alpha-skew-normal model with geochemical applications,” Applied Mathematics & Information Sciences, vol. 10, no. 5, pp. 1697–1703, 2016.

D. B. Owen, “Tables for computing bivariate normal probabilities,” The Annals of Mathematical Statistics, vol. 27, no. 4, pp. 1075–1090, 1956.

B. Arnold and G. Lin, “Characterizations of the skew-normal and generalized chi distributions,” Sankhya, vol. 66, pp. 593–606, 2004.

H. Ruskeepaa, Mathematica Navigator. Mathematics, Statistics, and Graphics, Academic Press, Cambridge, MA, USA, Third edition, 2009.

C. Brooks, RATS Handbook to Accompany Introductory Econometrics for Finance, Cambridge University Press, Cambridge, UK, 2009.

C. G. Broyden, “A class of methods for solving nonlinear simultaneous equations,” Mathematics of Computation, vol. 19, no. 92, pp. 577–593, 1965.

C. G. Broyden, “Quasi-Newton methods and their application to function minimisation,” Mathematics of Computation, vol. 21, no. 99, pp. 368–381, 1967.

D. F. Shanno, “Conditioning of quasi-Newton methods for function minimization,” Mathematics of Computation, vol. 24, no. 111, pp. 647–657, 1970.

Y. V. Marchenko and M. G. Genton, “A suite of commands for fitting the skew-normal and skew-t models,” The Stata Journal: Promoting Communications on Statistics and Stata, vol. 10, no. 4, pp. 507–539, 2010.

H. Bozdogan, “Model selection and Akaike’s Information Criterion (AIC): the general theory and its analytical extensions,” Psychometrika, vol. 52, no. 3, pp. 345–370, 1987.