Log-moment estimators for the generalized Linnik and Mittag-Leffler distributions with applications to financial modeling

We propose formal estimation procedures for the parameters of the generalized, three-parameter Linnik $gL(\alpha,\mu, \delta)$ and Mittag-Leffler $gML(\alpha,\mu, \delta)$ distributions. The estimators are derived from the moments of the log-transformed random variables, and are shown to be asymptotically unbiased. The estimation algorithms are computationally efficient and the proposed procedures are tested using the daily S\&P 500 and Dow Jones index data. The results show that the standard two-parameter Linnik and Mittag-Leffler models are not flexible enough to accurately model the current stock market data.

The main goal of this paper is to estimate the parameters of the heavy-tailed three-parameter generalized Linnik family of one-dimensional distributions, gL(α, δ, µ), with the characteristic function with µ > 0, δ > 0, and 0 < α ≤ 2. Another objective of this paper is to estimate parameters of the heavy-tailed three-parameter generalized Mittag-Leffler distribution gM L(α, δ, µ) (see, e.g., Laskin (2003)) with the Laplace transform with the corresponding density function is the generalized Mittag-Leffler function (see, e.g., Cahoy and Polito (2013)), with (η) r = η(η + 1) . . . (η + r − 1), η = 0, representing the classical Pochhammer symbol. We emphasize that estimation procedure for gM L(α, δ, µ = 1) was developed in Cahoy (2013). Note that if α = 1 and the data support is R + , we obtain the gamma distribution. When α = δ = 1 and and data support is R + , then we obtain the exponential distribution. It can be shown that gM L(α, δ, µ) is a mixture of generalized gamma densities with the strictly α + -stable density as the mixing distribution. With δ = 1, we have the usual Mittag-Leffler distribution (see. e.g., Pillai (1990)) which can be interpreted as a mixture of Weibull densities.
Finally, we also compare the efficiency of the above three-parameter models with the existing models (see e.g., Kozubowski (1999Kozubowski ( , 2001) using stock market S&P 500 and Dow Jones index data.
The paper is organized as follows: In Section 2, we provide structural representations of the generalized Linnik gL(α, δ, µ) and the generalized Mittag-Leffler gM L(α, δ, µ) random variables. In Section 3, we derive the method-of-moments estimators based on the log-transformed data. In Section 4, we test the algorithms using synthetic data. Section 5 shows the analyses of the S&P 500 and Dow Jones data. We conclude in Section 6 with a discussion of the key points of this work and possible future extensions of our study.

Mixture representations and moments
In this section we provide representations for random variables with generalized Mittag-Leffler and Linnik distributions employing the standard Lévy α-stable random variables and review related results for completeness.
The proof is straightforward: The proof can also be found, for example, in Pakes (1998). Recall that the α + -stable random variable can be conveniently generated using the classical Kanter (1975) where U is uniformly distributed in [0, π], and E is a standard exponential random variable (with rate/scale parameter one) independent of U. The q-th fractional moment of X can be easily derived from the above result, and is given below.
Theorem 2. Let 0 < α < 1. Then The proof follows directly from the standard moment formulas See also Cahoy and Polito (2013).

Generalized Linnik distributions on the entire real line
Theorem 3. Let 0 < α ≤ 2, and S α be a random variable with a symmetric αstable distribution with characteristic function exp(−|t| α ), and U be an independent gamma distributed random variable with density (2.1). Then the random variable has the gL(α, δ, µ) distribution.
The proof follows from the proof of Theorem 1. Note that Devroye (1990) had the proof for δ = 1/δ , µ = 1. Apparently, the case α = 1 is essentially different in both families.
The symmetric α-stable random variable S α can be generated using the standard Chambers-Mallows-Stuck (1976) formula where U 2 is uniformly distributed on [−π/2, π/2], and E is independent of U 2 and exponentially distributed with parameter one. An expression for the q-th fractional moment of Y is derived below.
Proof. Note that Using the q-th fractional moment of the symmetric stable random variable S α (see, Bening et al., 2004) .
3 Parameter estimation via the logarithmic moments 3.1 Generalized 3-parameter Mittag-Leffler distribution gM L(α, δ, µ) Following Cahoy, Uchaikin and Woyczynski (2010), we apply the log transformation to the random variable X given in (2.2) as where X = ln(X), U = ln(U ), and S = ln(S). For reproducibility, we can recall the first four log-moments of S from Zolotarev (1986), and Cahoy, Uchaikin and Woyczynski (2010): where C 0.5772 is the Euler's constant. It is straightforward to show the probability density of U as Using the polygamma function of order k, , the first four logmoments of U are Using the above moments, the estimating equations are as follows (see, also Cahoy and Polito (2013) where they were mentioned without showing the elementary (although tedious) algebra of moments): where ζ(·) is the Riemann Zeta function. Finally, using the estimatorsμ 3 andσ 2 X , we can solve the above equations for the variance and the third central moment, perhaps using a numerical software to obtain the estimatorsδ andα. Pluggingα andδ into the mean equation above, we obtain the following estimator of the parameter µ: 3.2 2-parameter Mittag-Leffler distribution gM L(α, 1, µ) We start by emphasizing that this two-parameter version is different from what had been studied in Cahoy, Uchaikin and Woyczynski (2010), which is gM L(α, 1, µ = λ −α ) in section 1, and from Cahoy (2013), which is gM L(α, δ, µ = 1). If δ = 1 then ψ(1) = −C, ψ (1) (1) = π 2 /6, and ψ (2) (1) = −2ζ (3). In addition, (3), and From the first two moments we obtain the following closed-form expressions of the estimators of α and µ: Note that these estimators are always non-negative as required and are asymptotically unbiased as shown in Proposition 2 below. where Proof. The proof directly follows from the asymptotic normality of sample moments and the multivariate delta method, where with the variance-covariance matrix and T is the gradient matrix. The above results can be used to approximate the (1 − ν)100% confidence intervals for α and µ.

Generalized 3-parameter Linnik distribution gL(α, δ, µ)
Applying the log transformation to the absolute value of the generalized Linnik random variable Y given in (2.5), we get an expression where S α = ln(|S α |). The first four integer-order log-moments of S α (see, Cahoy (2012)) are as follows: The moments above yield the same mean µ Y and the centered third order moment µ 3 as in the previous subsection. The variance then can be calculated to be Now the estimation approach employed in the previous subsection for the generalized Mittag-Leffler distribution gM L(α, δ, µ) can also be applied in the present case. The only difference here lies in the formula for the variance being used in the minimization process.

2-parameter Linnik distribution gL(α, 1, µ)
We start by emphasizing that this two-parameter version is different from what had been studied in Cahoy (2012), which is gL(α, 1, µ = λ −α ) in the first un-numbered equation in section 1. If δ = 1 then Moreover, we obtain the following closed-form expressions of the estimators of α and µ:α where Proof. The proof directly follows from Proposition 2 above where T , and the components of the variance-covariance matrix Σ are given in the beginning of this subsection. The above results can be used to approximate the (1 − ν)100% confidence intervals for α and µ.

Testing our estimation procedures on simulated data
In this section we will test the performance of our estimators using simulated data. Furthermore, to quantify the performance errors we will calculate the mean bias,

Generalized Mittag-Leffler distribution gM L(α, δ, µ)
For reproducibility, we used the optim function in R to minimize (σ 2 X −σ 2 X ) 2 +(µ 3 − µ 3 ) 2 with respect to α and δ using the initial value (0.1, 1). Note that expressions for σ 2 X and µ 3 are in Section 3.1. We emphasize that other built-in functions in R like the polygamma are used as well in the calculation process. However, the stable random variables are generated following Kanters formula (2.3) and C-M-S formula (2.6) due to their elegance.
The point estimates ofα andδ are then plugged in the point estimatorμ. From Table 1, the bias ofμ is around 19% when n = 10 3 and is around 6% when n = 10 4 . The CV fluctuates around 7.6% when n = 10 4 . Generally, Table 1 indicated positive results for the proposed method. Table 1: The mean bias and CV of the proposed estimators for the gM L(α, δ, µ) family using three different values of α, δ = 0.5, and µ = 1, for sample sizes n = 10 2 , 10 3 ,

Generalized Linnik distribution gL(α, δ, µ)
In this subsection we are providing results from testing our estimation procedures for the parameters in the gL(α, δ, µ) family. The approach is similar to the one we used for the generalized Mittag-Leffler distributions. The initial value pair used is (α 0 , δ 0 ) = (1, 1). We also calculated the same statistics for comparison. From Table  2, the bias went down to as little as 2.4% and went as high as 9.8% when n = 10 4 . The CV ranges from 3.2% to 12.4%. Note that the results for n = 100 suggest larger samples are needed or better optimization procedure (like the L-BFGS-B method in R). Also, the estimator for µ seems to get worse as the true α value approaches two. Overall, Table 2 provided favorable results for the proposed method especially for large samples. Note that in practice one can use bootstrap to quantify the variability of these estimators. In the entire analyses, we generated 1000 bootstrap samples to calculate the point and the 95% interval estimates of the parameters. We also used the boundary corrected kernel density estimate of the evmix package of R to compare the fits of gM L(α, δ, µ) and gM L(α, 1, µ) whenever possible. R scripts are available from the authors upon request.

Standard and Poor's (S&P) 500 index
It was originally called the "Composite Index" when it was first introduced as a stock market index in 1923. Three years later, the Composite Index expanded to 90 stocks, and then in 1957 -to its current 500, and renamed S&P 500 Index. It was the first index to be published daily. It contains 500 of the largest stocks in the United States. It is a benchmark for gauging the overall health of the large American companies, and the U.S. stock market in general. More than $7.8 trillion is benchmarked to the index (Source: Investopedia).
5.1.1 Comparison between gM L(α, δ, µ) and gM L(α, 1, µ) distributions We fitted the gM L(α, δ, µ) to the absolute values of the negative adjusted closing log returns (n = 9005) from the S&P 500 data. Table 3 clearly indicates that α is favored to be less than one and δ to be larger than unity, which suggests that a twoparameter Mittag-Leffler model is not adequate for this data. This observation is reinforced by the two-parameter estimates from the same table. In particular,α > 1 despite the relatively large sample size. The estimates of µ are however similar. Table 3. Parameter estimates for gM L(α, δ, µ) and gM L(α, 1, µ) models applied to (S&P) 500 data. To examine the model fit, we simulated data (sample size 2n = 18, 010) from the estimated model above. Specifically, we superimposed the boundary corrected kernel density estimates of the simulated data on the histogram of the observed data. Figure 1 below shows the good fit of the proposed model to the daily negative adjusted closing S&P 500 log returns. The graph demonstrates the advantage of flexibility of the three-parameter model as opposed to the two-parameter gM L(α, 1, µ) distribution in capturing the peak near the origin. Withα > 1, plotting the fit of gM L(α, 1,μ) is meaningless and computationally impossible.

Comparison between gL(α, δ, µ) and gL(α, 1, µ) distributions
We also analyzed the entire log adjusted closing returns (n = 17, 025). From the the estimates in Table 4, the estimates for δ favor values larger than one, which implies that the daily S&P 500 log returns (adjusted closing) are not adequately described by the two-parameter Linnik distribution (δ = 1). Note that we are not able to get an interval estimate for δ as the optim function gives the same value as a root for every bootstrap sample. The table also indicates that α is likely to be less than two andμ is way larger than the estimate obtained in the preceding section.
The two-parameter estimate of α exceeds two indicating a bad fit of the model to the data. But the estimates of µ from both two-and three-parameter models are comparable.  Figure 2 below confirms the good fit of the gL(α, δ, µ) family (using 2n = 33798 simulated observations) to the log adjusted closing returns. It also reveals that the flexibility of the proposed three-parameter gL(α, δ, µ) permits better capturing of the peak of the data at the origin. Note that the algorithm used in the calculation was not able to generate a comparable fit asα is way larger than two (upper bound). The analysis here is similar to the one we carried out in the previous subsection and deals with the absolute values of the negative adjusted closing log returns(n = 4, 359) from the Dow Jones index. The 95% CI for µ is between 154 and 173. The estimates of α strongly favor values less than unity. The point and interval estimates of δ indicate that δ > 1, which implies superiority of the generalized gM L(α, δ, µ) distribution for the absolute values of Dow daily Jones log returns over the gM L(α, 1, µ) distribution. Observe that the gM L(α, 1, µ) fit provides similar estimates for µ but not for α, and that its kernel density estimate is missing asα exceeds one. Again, as in the S&P 500 case discussed above, we constructed the graphs (using 2n = 8, 718 simulated observations) to investigate the model adequacy. The smoothed density of the gM L(0.983, 1.211, 162.596) is in Figure 3 . It basically confirms what we have already observed above, that is, the three-parameter model provides more flexibility in capturing the peak of the 'cupping' near the origin than the two-parameter Mittag-Leffler distributions.

5.2.2
Comparison between gL(α, δ, µ) and gL(α, 1, µ) distributions We also applied the generalized Linnik distribution to the whole adjusted closings of the daily Dow Jones log returns (n = 17, 025). Looking at the estimate of δ, δ > 1, it is clear that the daily Dow Jones log returns (using the adjusted closing) cannot be adequately described by the two-parameter Linnik distribution (δ = 1).
Moreover, the table below shows α to be likely less than 1.9.  Figure 4 shows the fits of the gL(α, δ, µ) and gL(α, δ = 1, µ) models. Notice that the algorithm was able to produce a comparable fit in this case even if α > 2. Furthermore, it validates the previous observation that the two-parameter gL(α, δ = 1, µ) is not adequate for the description of the daily Dow Jones adjusted closing log returns especially in capturing the peak at the origin.

Concluding remarks
The article proposes formal statistical inference procedures for the heavy-tailed three-parameter generalized Linnik and three-parameter generalized Mittag-Leffler families of distributions. The models provide considerable flexibility in modeling stationary discrete-time processes. The consistency and unbiasedness of the point estimators were computationally tested and seemed to be acceptable. Furthermore, the structural representations and the random number generation algorithms were provided for convenience. The paper provides guidance to how to distinguish different subcases of these models that exist in the literature.
The heavy-tailed three-parameter generalized Linnik and generalized Mittag-Leffler models present evidence that the adjusted S&P 500 and Dow Jones log returns  Figure 4: The histogram (using 150 bins) of the observed data and the kernel density plots (bandwidth = 0.001) of the simulated data (red dashed for 2-parameter gL(α, δ = 1, µ) and solid blue for 3-parameter gL(α, δ, µ) ) using the obtained estimates.
can obey these probabilistic laws. The comparison of the proposed three-parameter models with the two-parameter models clearly demonstrated inadequacy of the latter ones especially when one considers approximations around the origin in modeling daily log returns of stock market data (see also Kozubowski, 2001). Improvements of these procedures using robust, Bayesian approaches or more efficient algorithms, and the derivation of the trivariate or joint asymptotic distribution of the three point estimators would be worth exploring in the future.