In actuarial science or finance, we often encounter the problem of fitting distributions to data where the distributions have no closed-form expressions for their densities. These distributions are often infinitely divisible and they happen to be the distributions of the regularly spaced increments of Lévy processes. Beside infinitely divisible distributions, mixture distributions created using a mixing mechanism also provide examples of continuous densities without a closed-form expression. These types of distributions are often encountered in actuarial science. A few examples will be provided as illustrations subsequently.
Likelihood methods might be difficult to implement in such cases, due to the lack of a closed-form expression for the density function. To handle such a situation, we can consider the following approaches:
1) Expectation-maximization (EM) algorithm. Only under special conditions can the EM algorithm be used as it requires some conditional distributions, and these conditional distributions might be difficult to obtain; see McNeil, Frey and Embrechts  (pages 81-85) or McLachlan and Krishnan  .
2) Method of moments. Even though the model density has no closed form, if the model moments can be expressed in closed form, then the method of moments can be used. The main drawback of the method of moments is that estimators thus obtained might not be efficient nor robust for models with three or more parameters as the estimators will depend on a polynomial of degree three or higher, making the methods very sensitive to data which are contaminated; see Küchler and Tappe   for method of moments estimation.
3) The k-L procedure. Even if the density has no closed form, if the model characteristic function has a closed-form expression, then we can select points from the real and imaginary parts of the empirical characteristic function and match them with their model counterparts at the chosen points. This is the k-L procedure as proposed by Feuerverger and McDunnough  (pages 22-24).
4) Indirect inference. These methods are based on simulations and they require two steps. First, we need to choose a proxy model to obtain the estimators which are biased. Second, we remove the bias using simulations. See Garcia, Renault and Veredas  for this method. The proxy models from which the estimators are obtained affect the efficiencies of the estimators. For some models, it is difficult to know which proxy model will generate estimators with high efficiencies.
When implementing these methods for distributions without closed-form densities, there are some drawbacks which motivate us in this paper to extend minimum Hellinger distance methods originally proposed by Beran  to a simulated version (version S) which consists in replacing the model density by a density estimate using a random sample drawn from and minimizing
to obtain the simulated minimum Hellinger distance (SMHD) estimators, where is an empirical density estimate based on the observed data with the property where is the true vector of parameters. This consistency property will imply as ; see section 3 (page 224) of Tamura and Boos  .
Clearly, the new method proposed here will avoid the problem of arbitrariness in the choice of points for the k-L procedure based on characteristic functions. Unlike indirect inference, the proposed method does not need a proxy model. Furthermore, the estimators obtained using the proposed method might be more robust and efficient than method of moments estimators. Besides, the proposed method does not require conditioning, which can be difficult, whereas the EM algorithm does.
It appears that the proposed method, which originally combines simulation with Hellinger distance, adds to the set of statistical techniques that can be useful for financial and actuarial data, yet many of which do not receive much attention in the actuarial literature. SMHD methods depend on being able to draw samples from the parametric family; in general, this is indeed possible. Consequently, SMHD methods also add to the existing literature on simulated inference which is relatively new; see comments by Davidson and MacKinnon  (page 393).
The new method is built on the classical version (version D) of Hellinger distance as proposed by Beran  which consists in minimizing
to obtain the minimum Hellinger distance (MHD) estimators. The MHD estimators have been known to have nice robustness properties with breakdown point greater than 0. Also, they are consistent with, in general, less stringent conditions for consistency than maximum likelihood (ML) estimators. However, more restrictions are placed upon the underlying parametric family for the MHD estimators to attain full efficiency, such assuming having a compact support for example. Despite this drawback, simulation studies often show that the methods perform well across many models. For a literature review of Hellinger distance (HD) methods, see chapters 3 and 10 of the book by Basu, Shioya and Park  . From the literature, it can be seen that HD methods still do not receive proper attention for their use in actuarial science and finance, especially in actuarial science.
In this paper, we introduce a simulated version of HD methods and show that the SMHD estimators are consistent. However, the question of asymptotic normality is still not resolved for the time being. Further work should generate results on asymptotic distributions for the SMHD estimators that shall then be presented in a subsequent paper. In this paper, the methods are presented with fewer technicalities and we relate them with the traditional likelihood methods. In doing so, we wish to encourage practitioners to use these methods for their applied works in their fields. In the next paragraphs, we will consider a few examples for illustrations of the types of distributions without closed-form expressions often encountered in finance and actuarial science where the new simulated method can be particularly useful.
We present here the class of normal mean-variance mixture distributions where the random variable can be represented using equality in distribution as
1) , and are parameters with , , and ;
2) is a nonnegative random variable with an infinitely divisible (ID) distribution;
3) follows a standard normal distribution and is independent of .
The generalized hyperbolic, variance-gamma, and normal-inverse Gaussian distributions belong to this class; see McNeil, Frey and Embrechts  (pages 77-79). By conditioning on first, the moment generating function (mgf) for
can be obtained and given by , where the
moment generating functions of and are given respectively by and . Distributions of the increments observed at regular intervals of a subordinated Brownian motion process belong to this class. It can easily be seen that the density function of depends on the density function of . Consequently, the density function of might not have a closed-form expression in general. Closely related to the variance-gamma distribution is the generalized normal-Laplace (GNL) distribution which is introduced by Reed  and is given in the next example.
A random variable follows a GNL distribution if it can be represented as
1) the parameters are , , , and , with , , , , and ;
2) the random variables and are independent and follow a common gamma distribution with density function ;
3) follows a standard normal distribution, , with being independent of and .
The distribution is infinitely divisible and can display asymmetry and fatter tail than the normal distribution. It will be symmetric if . The vector of parameters is and the mgf for can be obtained using the representation given by Equation (4) and is given by
From the cumulant generating function, the mean and variance are given respectively by
Higher cumulants are
for . (8)
Due to the lack of a closed-form expression for the density function, Reed  (page 477) has proposed using the method of moments and matching the empirical cumulants with the model cumulants to estimate the parameters. He applied the method to data collected on stocks. In the particular case with four parameters, where , moment estimators can be obtained explicitly. However, for the general case with five parameters, the moment equations must be solved numerically. The moment estimators will be discussed in more detail in section 3 and we shall compare their efficiencies with the efficiencies of the SMHD estimators based on simulated samples.
For more on Lévy processes and infinitely divisible distributions used in finance, see chapter 6 of the book by Schoutens  (pages 73-83). For nonnegative infinitely divisible distributions used in actuarial science, see Dufresne and Gerber  , and Luong  . For mixtures of distributions without closed-form density functions, for which the proposed estimators can also be used, see Klugman, Panjer and Willmot  (pages 62-65). We shall consider HD estimation in all those cases.
Assume that we have a random sample of observations and they are independent and identically distributed as the random variable which is continuous with model density given by . The vector of parameters is denoted by . In his seminal paper, Beran  proposes to estimate by the minimum Hellinger distance estimators denoted by which minimize, with respect to , the Hellinger distance between a consistent empirical density estimate and the parametric family with the property pointwise. It leads to minimize the objective function
Beran  also noted that, intuitively, the methods are robust as data are smoothed by a kernel density estimator , and hence the effects of outliers are mitigated. It has been confirmed in various models that the asymptotic break-
down points of the estimators are around and it is well-known that the sam-
ple mean has a breakdown point of 0. See Hogg, McKean and Craig  (pages 594-595), and Maronna, Martin and Yohai  (page 58) for the notions of finite sample and asymptotic breakdown points as measures of robustness of estimators. See Lindsay  for the discussions on robustness and efficiencies of MHD estimators. We also note that, since
and, using the Cauchy-Schwarz inequality, , we find
Moreover, since if and only if
almost everywhere, it implies if and only if almost everywhere.
The objective function is stable and bounded. This might explain why, intuitively, minimizing such an objective function, we obtain estimators that are also stable and therefore robust in some sense.
Kernel density estimators are often used to define . One of the simplest kernel density estimators is the rectangular kernel density estimator which generalizes the usual histogram estimator. In general, kernel density estimators have the form
a) is the bandwidth with the property that and as ;
b) is a density function.
The property specified by a) guarantees the consistency of ; see Corollary 6.4.1 given by Lehmann  (pages 406-408). Subsequently, we implicitly assume that density estimates used with the SMHD method meet the requirements specified by a) and b).
For the rectangular kernel density, the following symmetric density around 0 is chosen with for . The kernel has a compact
support. The density estimate at is then the average of rectangles located within units from . For other kernels and their implementation using the package R, see chapter 10 of the book by Rizzo  (pages 281-318). For Hellinger distance estimation, it is preferable to use a symmetric kernel with a compact support and twice differentiable for meeting the regularity conditions of Theorem 4 as given by Beran  (pages 450-451); also see the discussions by Basu, Shioya and Park  (pages 78-83). In this paper, we only need univariate kernel density estimates but multivariate density estimates based on kernels can also be defined similarly; see Toma  and Scott  .
If has no closed-form expression but random samples can be drawn from the distribution with density , clearly we can use the same type of kernel density estimator, used to define , to estimate . In other words, in order to estimate , we similarly define as being the kernel density estimator based on a random sample of size . Note that as and needs to be reasonably large so that there is little loss of efficiencies due to simulations; we recommend .
Consequently, for the simulated version, we shall minimize the objective function given by
to obtain the SMHD estimators.
For terminology, we shall call the classical version, which is deterministic in terms of , version D, and the simulated version, version S. Since , as given by Equation (13), is not differentiable, a direct simplex search method which is derivative-free is recommended. The R package already has a built-in function for performing the Nelder-Mead simplex method which is a derivative-free method to minimize a function. Also, there is a built-in function to handle density estimates using various kernels. These features will facilitate the implementation of SMHD methods for applied works by practitioners. Furthermore, because the densities and based on a rectangular or triangular kernel are positive only in some finite interval and zero elsewhere, this makes the integration for evaluating Equation (13) easy to handle. A trapezoid quadrature method will suffice to find the SMHD estimators. Note that for the simulated version, we still have
As data are also smoothed, intuitively, these features will again make the simulated version robust.
The paper is organized as follows. In Section 2, we will look into the asymptotic properties of MHD estimators. More precisely, we shall briefly review the asymptotic properties of the classical MHD estimators in Section 2.1 and establish the consistency of SMHD estimators in Section 2.2. Also in Section 2.2, an estimator for the Fisher information matrix is proposed with the use of SMHD estimators. In Section 3, we use a limited simulation study to compare the efficiencies of the SMHD estimators with those of method of moments estimators, using the GNL distribution. Despite being limited, the study seems to show that the SMHD estimators are more efficient than the method of moments estimators. This seems to point to the potential of SMHD methods to generate estimators with good efficiency and further justify their use in actuarial science and finance.
2. Asymptotic Properties
2.1. Asymptotic Properties of the Classical MHD Estimators
MHD estimators can be seen to be consistent in general for version D and version S. In fact, the conditions are even less restrictive than the conditions for maximum likelihood estimators to be consistent. Since we aim for applications, we only consider asymptotic properties under the strict parametric model, i.e., assuming the observations come from the parametric density family , where , and the parameter space is assumed to be compact.
where and are density functions. Note that is a norm in the density functional space and it will respect the triangular inequality.
Tamura and Boos  (page 224) have noted that, if , then , and if for in probability,
it is sufficient for the MHD estimators given by the vector obtained by minimizing Equation (10) to be consistent, i.e., , assuming the parameter space is compact. See Theorem 3.1 by Tamura and Boos  (page 224). Comparing with the regularity conditions for ML estimators as given by Theorem 2.5 of Newey and McFadden  (page 2131), the regularity conditions for MHD estimation do not require that as in likelihood estimation. This makes the MHD estimators consistent in general even with fewer restrictions than ML estimators.
However, for asymptotic normality, they require more stringent conditions to be as efficient as ML estimators. They are found in Theorem 4 given by Beran  (pages 450-451), which is summarized in Theorem 1 below, focusing on the strict parametric model. Beran  (pages 450-451) allows the bandwidth of the kernel to be randomly chosen with , where is a sequence of constants but is a sequence of random variable with . It also requires
a compact support K for both and . Despite these restric-
tions, empirical studies often show that the estimators have high efficiencies in many models without the condition of compact support for the parametric family met. The regularity conditions of Beran’s Theorem 4 when restricted to the strict parametric model are stated using Theorem 1 below. We also require the vector of true parameters to be in , where is compact. Theorem 1 can be viewed as a corollary of Theorem 4 as given by Beran  and the proofs have been given there.
1) The kernel density is symmetric about 0 and has a compact support.
2) The function is twice differentiable and its second derivative is bounded on the compact support.
3) and have a compact support and on .
4) is twice absolutely continuous with its second derivative with respect to being bounded.
5) , , and .
6) There exists a positive constant which might depend on such that is bounded in probability.
Then where is the Fisher information matrix with elements given by
and assumed to exist.
We just give an outline establishing the results of Theorem 1 and focus only on the strict parametric model for applications with the aim that it might help practitioners in the applied fields to follow more easily the arguments needed to develop the new method subsequently.
Note that, beside the rectangular kernel, the triangular kernel with for and the Epanechnikov kernel with for
meet conditions 1 and 2 as required by Theorem 1 and are available in the package R.
For establishing asymptotic normality results for the estimators as indicated by Theorem 1, we can consider a Taylor expansion of the system of equations
around the true vector of parameters . The system of equations implies
with and .
We proceed to perform a Taylor expansion by noting
assuming is differentiable with respect to and
, with , using the compact support assumption for . As a result, we can write that
Therefore, with the regularity conditions met, we will have the representation
where is the remainder term which converges to 0 in probability, which can be re-expressed using the following equality which holds in law,
Using the argument given by Beran  (page 451) allows us to establish the equality in probability,
This can be viewed as a form of generalized delta method to establish equality of the left-hand side and the right-hand side of Equation (23).
Consequently, Equation (22) can be re-expressed, using the equality in distribution, as
as, in general, . Furthermore,
where is the commonly used sample distribution function. This allows the following representation:
For the simulated version, i.e., version S, we can only obtain results for consistency and they will be given in the next section. As for asymptotic normality, we cannot conclude for the time being whether or not conditions of Theorem 7.1 given by Newey and McFadden  (pages 2185-2186) for the asymptotic normality of estimators obtained from a non-smooth function can be met. We hope to have more results on this issue in the future and would present them in a subsequent paper. This does not prevent SMHD estimation from being used as an alternative to methods of moments if the primary interests are in point estimation.
2.2. Asymptotic Properties of the SMHD Estimators
For version S, we minimize
We recommend using the same seed across different values of if possible and the simulated sample size such that at the same rate as . These recommendations conform with other simulated methods of inference such as the method of simulated moments as discussed by Davidson and McKinnon  (page 284) or simulated quasi-likelihood found in Smith  (page S68). The condition of the same seed being used is not necessary for consistency, but it allows to have the value for each fixed each time we want to evaluate ; otherwise, the values might differ slightly due to the fact that simulations are needed to evaluate . With the same seed, behaves like a non-random function with respect to .
with as defined by Equation (28). The following Theorem, which is essentially Theorem (3.1) given by Pakes and Pollard  (page 1038) with the assumption of compactness of the parameter space added, can be used to establish the consistency of SMHD estimators. The proofs of the following Theorem have been given by Pakes and Pollard  (page 1038) using the Euclidean norm. Their proofs are still valid with the norm as defined by Equation (29) and discussed in Section 2.1. It is implicitly assumed that there is no identification problem for the parametric family, i.e., if , then except on a set of measure zero.
1) The parameter space is compact, and .
2) minimizes or equivalently .
3) for each , is bounded in probability, where denotes the norm being used.
Clearly, we have consistency for as and only at .
For the time being, we cannot assert that follows a multivariate normal distribution asymptotically as we cannot verify the regularity conditions of Theorem 7.1 given by Newey and McFadden  (pages 2185-2186) for estimators obtained from a non-smooth objective function. For the simulated unweighted minimum chi-square, Pakes and Pollard  (page 1049) find the as-
ymptotic covariance to be , with being the asymptotic covariance
matrix of the estimators without using simulations. Conforming with other simulated methods which typically give the same type of asymptotic covariance formula, we recommend choosing to minimize the loss of efficiency due to simulations. The matrix , where , can be viewed as a form of benchmark for the approximate asymptotic covariance matrix for if indeed asymptotic normality can be shown. In the absence of a rigorous proof, we have to rely on simulations to evaluate the efficiency of , just as for version D when the support of the distribution is not compact. Further asymptotic results to be obtained in the future will be presented in a subsequent paper.
Since we have estimates for densities, it is natural that we can estimate the Fisher information matrix. Clearly, if the model density has a closed-form expression, then the following matrix
can be used to estimate . Instead of , if it is not available, we can use the kernel density estimate of , and, following a method given by Pakes and Pollard  (page 1043), we can use
with at the rate , where , to estimate , assuming . The vector is a unit vector with 1 in its j-th place and 0 elsewhere. Replacing and
by these estimates will give an estimator for the information matrix. An estimate of the information matrix is useful as the information matrix is related to the Cramer-Rao lower bound.
3. Limited Simulation Study
In this study, we shall compare the efficiencies of the moment estimators for the case with , i.e., the GNL distribution with only four parameters. Reed  (page 477) has given the expressions for the moment estimators using the first six empirical cumulants , with the sample mean . They can be
obtained using central empirical moments as
they follow the same type of relationships which exist between model cumulants and model central moments. Let , and . The following relationships can be found in Stuart and Ord  (pages 90-91) and they are given by
Explicitly, the moments estimators are
, , and . (33)
Reed  (page 477) also notes that method of moments estimators (MM estimators) can take on negative values for positive parameters, and it is not easy to include constraints in method of moments estimation. Also, the use of EM algorithm does not appear to be straightforward for the GNL distribution. SMHD estimation can handle constraints by minimizing the objective function, which is given by Equation (28), with constraints.
A limited simulation study using parameters for the symmetric GNL distribution with four parameters, focusing on parameters in the ranges , , , , has been carried out and the relevant results are summarized in Table 1. The ranges of parameters as indicated are chosen accordingly and conform with the empirical study conducted by Reed  (page 481) using stock data. The simulated sample size for data is and the simulated sample size drawn from the model for SMHD estimation is , hence with . It takes about twenty minutes on a laptop computer to obtain the estimators for samples and, due to the limited computer capacity, we fix samples for each combination of parameters for our study. As we only have access to laptop computers, the scale of the study is limited.
We noticed that the method of moments estimator for
is often negative and we set it equal to zero whenever this is the case, and the comparisons of efficiencies use this version of the method of moments estimator. The density estimate is based on the built-in function of the package R with a rectangular kernel and default bandwidth based on the normal distribution. The overall asymptotic relative efficiency (
with being the commonly used mean square error of the estimator and it is estimated using samples for estimating the expression for ARE and the values of the estimated ARE’s using different sets of parameters are displayed in Table 1.
Despite the scope of the study being limited, it suggests that SMHD estimators perform much better than method of moments estimators overall for the ranges of parameters used in finance. The method of moments estimator for tends to perform better for small values of and deteriorates rapidly as grows larger with even for various parameter values that we tested which lie outside the ranges indicated above and not shown in Table 1. Table 1 is used for illustration and provides a summary of the key findings of the study. Also, in the ranges considered, the method of moments estimator for tends to perform better than its SMHD counterpart, but the overall efficiency of MM estimators still falls behind the overall efficiency of SMHD estimators in general as shown in Table 1. Clearly, more work needs to be done numerically and theoretically, but it shows the potential efficiencies of SMHD methods.
Table 1. Asymptotic relative efficiencies to compare SMHD estimators with MM estimators with .
Note: Tabulated values are estimates of the asymptotic relative efficiencies of the SMHD estimators versus the MM estimators.
Individual ratios of mean square errors for some sets of parameters
, , , ,
, , , ,
, , , ,
As SMHD estimators remain consistent with minimum regularity conditions and despite the lack of results on asymptotic normality, the proposed method appears to be useful for fitting actuarial and financial models using continuous infinitely divisible distributions which arise from Lévy processes or continuous mixture distributions constructed using mixing operations, whenever it is not difficult to simulate from these distributions but the density functions of these distributions have no closed-form expressions. In many models, the proposed method appears to be more efficient than traditional methods such as the method of moments. The proposed method is not difficult to implement but methods based on simulations do not seem to receive much attention in finance and actuarial science. They might be considered as additional robust statistical techniques for analyzing empirical data, especially if point estimation is the main interest.
The helpful comments of an anonymous referee and the kind support of the OJS staff, which led to an improvement in the presentation of the paper, are gratefully acknowledged.
 Küchler, U. and Tappe, S. (2008) Bilateral Gamma Distributions and Processes in Financial Mathematics. Stochastic Processes and their Applications, 118, 261-283.
 Tamura, R.N. and Boos, D.D. (1986) Minimum Hellinger Distance Estimation for Multivariate Location and Covariance. Journal of the American Statistical Association, 81, 223-229.
 Dufresne, F. and Gerber, H.U. (1993) The Probability of Ruin for the Inverse Gaussian and Related Processes. Insurance: Mathematics and Economics, 12, 9-22.
 Luong, A. (2016) Cramér-Von Mises Distance Estimation for Some Positive Infinitely Divisible Parametric Families with Actuarial Applications. Scandinavian Actuarial Journal, 2016, 530-549.
 Newey, W.K. and McFadden, D. (1994) Large Sample Estimation and Hypothesis Testing. In: Engle, R.F. and McFadden, D.L., Eds., Handbook of Econometrics, Volume 4, North Holland, Amsterdam, 2111-2245.