key: cord-0221743-p9mrf146
authors: A.H.Nzokem,
title: Fitting Infinitely divisible distribution: Case of Gamma-Variance Model
date: 2021-04-15
journal: nan
DOI: nan
sha: fd5e278d95afbef9998e7e08331040388e732d2d
doc_id: 221743
cord_uid: p9mrf146

The paper examines the Fractional Fourier Transform (FRFT) based technique as a tool for obtaining probability density function and its derivatives, and mainly for fitting stochastic model with the fundamental probabilistic relationships of infinite divisibility. The probability density functions are computed and the distributional proprieties such as leptokurtosis, peakedness, and asymmetry are reviewed for Variance-Gamma (VG) model and Compound Poisson with Normal Compounding model. The first and second derivatives of probability density function of the VG model are also computed in order to build the Fisher information matrix for the Maximum likelihood method. The VG model has been increasingly used as an alternative to the Classical Lognormal Model (CLM) in modelling asset price. The VG model with fives parameters was estimated by the FRFT. The data comes from the daily SPY ETF price data. The Kolmogorov-Smirnov (KS) goodness-of-fit shows that the VG model fits better the empirical cumulative distribution than the CLM. The best VG model comes from the FRFT estimation.

Empirical studies have shown that asset returns are often characterized by peakedness, leptokurtosis and asymmetry. Most of these empirical properties were reviewed in [1] as a set of stylized empirical facts that emerge from price variations in various types of financial markets. These facts provide evidence suggesting the assumptions of the Classical Lognormal Model (CLM) are not consistent with the empirical observations. In order to reduce the theoretical-empirical gap, it has been used in literature two sophisticated stochastic processes with the fundamental probabilistic relationships of infinite divisibility. A natural generalization of the Brownian motion, which is the Lévy process. As a stochastic process, the Lévy process has independent and stationary increments; and as a function, the process is right continuous with left limit. However, the Lévy process model is usually incomplete, as the associated risk-neutral equivalent martingale measure is usually not unique; and this situation leads to many different possible prices for European options [2, 3] . Another way to generalize the Classical Lognormal Model (CLM) is the method of subordination [3] [4] [5] . The subordinated process is obtained by substituting the physical time in the CLM by any independent and stationary increments random process, called the subordinator. If we consider the random process to be a Gamma process, we have a Variance Gamma (VG) model, which is the model the paper will be investigating. The Variance Gamma (VG) model was proposed by [6] as an alternative to the Classical Lognormal Model (CLM). In contrast to the CLM, the VG model does not have an explicit closed-form for probability density function and its derivatives. However, the fundamental probabilistic relationship of infinite indivisibility allows the characteristic function of the VG model and the Fourier of probability density function to have a closed-form function. In the study, the VG model has five parameters: parameters of location (µ), symmetric (δ ), volatility (σ ), and the Gamma parameters of shape (α) and scale (θ ). The VG model density function is proven to be (1.1).

The integral in (1.1) makes it difficult to utilize the density function and its derivatives, and to perform the Maximum likelihood method. However, in the literature, many studies have found a way to circumvent the lack of closed-form by decreasing the number of parameters and making use of approximation function or analytical expression with modified second kind Bessel function. In fact, [7] developed a procedure to approximate (1.1) by Chebyshev Polynomials expansion. [8] and [9] got (1.1) by analytical expression with modified Bessel function of second kind and third kind respectively. [3] got (1.1) through Gauss-Laguerre quadrature approximation with Laguerre polynomial of degree 10. [5] used the Fast Fourier Transform (FFT). The paper will implement the Fractional Fourier Transform (FRFT)-based technique on the Fourier Transform of the VG model function in order to obtain probability density function and its derivatives. Contrary to the Fast Fourier Transform (FFT), the advantage of using the FRFT mentioned in [10] can be found at three levels:(1) both the input function values and the output transform function values are equally spaced; (2) a large fraction of the density function is either zero or smaller than the computer machine epsilon; and (3) only a limited range of the density function are required. The paper is structured as follows, the next section presents the analytical framework, with an overview of the FRFT computing formulas, their applications on three infinitely divisible distributions, and a review of the Maximum likelihood method. The third section will present the Variance Gamma (VG) model and the sample data, before performing the parameter estimations of the VG model and the Kolmogorov-Smirnov (KS) goodness-of-fit test.

The continuous Fourier transform (CFT) of function f (t) and its inverse are defined by:

where i is the imaginary unit. The Fast Fourier Transform (FFT) is commonly used to evaluate the integrals (2.1) and (2.2) . The fundamentally inflexible nature [10] of FFT is the main weakness of the algorithm. The advantages of computing with the FRFT mentioned in [10] can be found at three levels: (1) both the input function values f (x k ) and the output transform values F [ f ](x k ) are equally spaced; (2) a large fraction of f (x k ) are either zero or smaller than the computer machine epsilon; and (3) only a limited range of f (x k ) are required. The FRFT is set up on n-long sequence (x 1 , x 2 , . . . , x n ) by:

For δ = 1 n , we have the Discrete Fourier Transform (DFT). In the general scheme, the parameter δ is not limited to fractions; it can take complex number. It is shown in [11] that G k (x, δ ) is a composition of DFT −1 and DFT .

where y and z are two 2n-long sequences defined by:

As shown in (2.4), the FRFT is efficiently computed through the FFT on 2n-point. We assume

, and β = a n is the step size of the n input values

We have also γ as the step size of the n output values of f (t) and x k = (k − n 2 )γ for 0 ≤ k < n. By choosing the step size β on the input side and the step size γ in the output side, we fix the FRFT parameter δ = β γ 2π . The density function f at x k can be written as (2.7). The proof is provided in Appendix A.1.

In order to perform f (t) function from the Fourier Transform (FT), we assume a = 20, n = 2048, β = γ = a n . For more detail on FRFT, See [10] , [11] .

The FRFT probability density functions are computed for each of the three infinitely divisible distributions: Normal distribution with two parameters, Variance Gamma (VG) Distribution with five parameters, and Compound Poisson with Normal compounding with three parameters. The distributional properties are reviewed and the results also show that the Kurtosis statistics tell very little about the peakedness [12] . For detailed proof of infinite divisibility for each distribution. See [13] , [14] , [15] 2.1.1 Normal distribution.

The normal distribution is an infinitely divisible distribution. It was shown in [13] that the normal distribution N(µ, σ ) is the n th convolution power of a probability measure N( µ n , σ n ), for every n. The density and the Fourier transform functions have an explicit closed form as shown below.

For mean (µ = −2) and standard deviation (σ = 1), Fig 1a displays 

The highest computation error is 1.8738e − 15, which is negligible. The probability density function is symmetric with the Kurtosis statistic 3 for the normal distribution.

2.1.2 Variance Gamma (VG) Distribution.

The Variance Gamma distribution is an infinitely divisible distribution. The density function is (2.11) and does not have a closed form. However, the Fourier transform function has explicit closed form (2.10). See Appendix A.1 for proof

When δ = 0, we have Symmetric Variance Gamma (SVG) Model, which is a special case of the Variance Gamma (VG) Model. It can be shown through Cumulant-generating function [16] that

For Parameter values: µ = −2, δ = 0, σ = 1, α = 1, θ = 1. Fig 2 displays the FRFT estimation of the probability density functions. As shown in Fig 2b, the probability density is left asymmetric and right asymmetric when the parameter (δ ) is negative and positive respectively. For δ = 0, the density function is symmetric, as shown in (2.12). The shape parameter (α) impacts the peakedness and tails of the distribution, as illustrated in Fig 2c and (2.12); heavier is the tails, shorter is the peakedness. θ and σ have the same impact on the distribution. As shown in (2.12), both change only the variance. (c)f (t) and shape parameter (α) 

The density function was shown in [17] to be (2.14). However, the Fourier Transform density has explicit closed form (2.13). See Appendix A.1 for proof.

It can be shown through Cumulant-generating function [16] that

For Parameter values: µ = −0.1, σ = 0.25, λ = 16. The FRFT produces the density functions in Fig 3. As shown in Fig 3b, the probability density is asymmetric, shifts to the left and to the right when the parameter (µ) is negative and positive respectively. For µ = 0, we have a symmetric density function. µ and σ impact the symmetry, peakedness and tails of the distribution, as illustrated in Fig 3c and (2.15); heavier is the tails, shorter is the peakedness. The Poisson parameter (λ ) impacts all the statistics in (2.15). λ affects the asymmetry, peakedness and tails of the distribution simultaneously. (c)f (t) and Volatility parameter (σ ) From a probability density function f (x,V ) with parameter V of size p and the sample data x of size m, we definite the Likelihood Function as shown in (2.16) .

In order to perform the Maximum of the likelihood function (2.16), the quantities dl(x,V ) dV j and

dV k dV j must be computed. We need to compute d f (x,V ) dV j

and d 2 f (x,V ) dV k dV j , which are the first and second order derivative of the density function respectively. Variance Gamma (VG) distribution has five parameters (µ, δ , σ , α, θ ), using the closed form of Fourier transform in (2.10) and the FRFT algorithm, we perform the first derivative functions: 

dV k dV j in (2.19) are critical in computing the Hessian Matrix and the Fisher Information Matrix. For VG distribution, We need fifteen second derivative functions. The results of four functions are shown in Fig 5. The remaining functions are reported in Appendix B.2 and Appendix C.3. 

2 and Appendix C.3. The remaining derivatives have rotational symmetry at the point (−2,0). These features can be checked analytically through the Inverse Fourier Transform in (2.2) and F [ f ] in (2.10). Given the parameters V (µ, δ , σ , α, θ ) and the sample data set x, from the previous development and computations, we have the quantities:

Maximizing L(x,V) in (2.16) with respect to V is equivalent to solving the equation system (2.22) .

The solutions of (2.22) are provided by the Newton-Raphson Iteration Algorithm (2.23). In fact, the Taylor's expansion can be applied to the score function to derive the Newton-Raphson Iteration Algorithm.

For more detail on Maximum likelihood and Newton-Raphson Iteration procedure, see [18] , [19] , [20] , [21] , [22] .

The VG model was introduced by Madan et al [6, 9] . The idea [23] was to model the asset price on bussiness time rather than on calender time. The asset price is modelled on bussiness time (k) as follows:

Where µ, δ ∈ R , σ > 0, α > 0 and θ > 0. µ is the drift of the physical time scale, δ is the drift of the activity time process, and σ is the volatility. {T k } is the activity time process [24] or intrinsic time process [5] and is supposed to be non-negative stationary independent increments. This process is associated with either the flow of new price-sensitive information or the cumulative trading volume process. Y k is the return variable of the stock or index price, we have (3.3) from (3.1) and (3.2) .

The density of Y j in (3.1) and (3.3) has the following expression

The Fourier transform of (3.4) can be shown to be .

The relationship between means (E(Y j )), variance (Var(Y j )), Skewness (Skew(Y j )) and Kurtosis (Kurt(Y j )) indicators and the VG model with five parameters can be shown by Cumulantgenerating function [16] .

By applying the method of Moments, the system of equations has more parameters than equations. The issue can be overcome by imposing [3] E(V ) = 1 or σ = 1. E(V ) = 1 was chosen in [9, 25] and the implication was that the gamma shape and scale parameters should be equal (α = θ ). In the next section, we will assume σ = 1 in order to solve the system of equation (3.7). . . (3.10) . However, the constraints will be relaxed in the Maximum Likelihood method where all five VG parameters will be estimated.

The data comes from the SPY ETF, called SPDR S&P 500 ETF (SPY .11).

One observation will be lost in the computation process. The results of the daily SPY ETF return are shown in Fig 7. The volatility of the SPY ETF is higher in the First quarter of 2020 as The SPY ETF is an inexpensive way for investors to gain diversified exposure to the U.S. stock market; like others stocks and securities, it was unusually volatile in 2020 amid the coronavirus pandemic and massive disruptions in the global economy. 13 daily return observations was identified as outliers and removed from the data set in order to avoid negative impact on the statistics. See Table 1 below. By assuming the independent, identically distribution in the SPY ETF return sample, the method of moments was used to compute statistics in Table 1 for both samples with outliers and without outliers. The outlier effects are greater on Kurtosis, Skewness, and Variance. The sample Skewness (0.4687) is between -0.5 and 0.5, which is the skewness interval for symmetric [5, 26] . 

17.3572 6.6853

The shape of the daily SPY ETF return distribution are estimated from the sample using two non parametric density estimators: Histogram in Fig 8a and Based on the VG model assumption that σ = 1, the equations (3.7)-(3.10) was solved given the statistics provided in Table 1 . The parameter values µ, δ , α, and θ are shown in Table 2 . We have two models, the asymmetric Variance-Gamma model (AVG) with chronometer drift negative (δ = −0.1399); and the symmetric Variance-Gamma model (SVG) without drift. As shown Table 2 , when δ = 0, the drift of the physical time scale µ decreases, but keeps the positive sign. The magnitude of α and θ remains the same. The Gamma shape parameter (α) is less than one. In fact, The larger the sample kurtosis, the smaller the Gamma shape parameter (α) as shown in the relationship (Kurtosis = 3(1 + 1 α )) for SVG model. 

Regarding the estimation by Maximum likelihood method, the previous assumption on the volatility (σ = 1) doest not hold and was dropped. The method has a local property that makes the result depending on the initial value. The estimations by Moments method were used as initial values in the maximization procedure. In table 3, the results are labeled AVG1 for the Asymmetric Variance-Gamma Model and SVG1 for the Symmetric Variance-Gamma Model. Another initial value was chosen with σ = α = θ = 1, and δ = µ = 0. The results are labeled AVG2 and SVG2 respectively for Asymmetric Variance-Gamma Model and Symmetric Variance-Gamma Model. As shown in Appendix D.4, the maximization procedure convergences after 21 iterations For AVG2, as shown in Table 5 . SVG2 also converges, after 15 iterations. See Table 6 in Appendix D.4. It appears that the drift of the physical time scale µ is positive, and the drift of the activity time process δ is negative. The estimation of the Maximum likelihood underestimates the value of the kurtosis.

For the Classical Lognormal Model (CLM), the Method of Moments and the Maximum likelihood yield the same estimation of the mean (μ) and standard deviation (σ ), given a high sample size. The results are displayed in Table 3 . It is shown in (3. 3) that the VG model is a special case of the CLM. It is also known that the parameter values in Table 3 maximized the likelihood function. Therefore, it becomes appropriate to perform the Likelihood Ratio test between the normal distribution and the four other VG models. The Likelihood Ratio estimator (W) is defined as follows.

. P values is the probability that the real values are even more extreme (more in the tail) than our test statistic (Ŵ ). For Likelihood Ratio test review, see section 5.5 in [22] , and [20] , [27] The value of the Log-likelihood (Log(ML)) is displayed in Table 3 for each model. There is a huge discrepancy between the Log-likelihood of the CLM and VG model, which produces P value almost zero. Therefore, the null hypothesis (H 0 ) of normal distribution, is rejected based on the sample data. The results provide evidence that the Variance-Gamma model is an improvement compare to the classical Normal model. 

See Appendix A.1 for (3.15) proof. The cumulative distribution function (F) was computed from the inverse of the Fourier of the cumulative distribution (F [F] ). The results for SVG2 model is shown in Appendix E.5, Fourth column. The same methodology was followed for the remaining VG models. The Kolmogorov-Smirnov (K-S) estimator (D n ) is defined as follows.

where n is the sample size, F n (x) denotes the empirical cumulative distribution of {y 1 , y 2 . . . y n } . The distribution of the Kolmogorov's goodness-of-fit measure D n has been studied extensively in the litterature [28] , [29] , [30] . It was shown that D n distribution is independent of the theoretical distribution (F) under the null hypothesis (H0). More recently in [31] , it was determined numerically the exact distribution of D n . [31] also provides a package KSgeneral [32] in the R software for computing the cumulative distribution of D n when F(x) is continuous and n is a natural number. The cumulative distribution of D n under the null hypothesis was computed and the density function was deduced. The computed density function is shown in Fig 9. Under the null hypothesis (H0), D n has a positively skewed distribution with means (µ = 0.0165) and standard deviation (σ = 5 * 10 −3 ). [33, 34] , d n can be estimated as follows.

The results for SVG2 model is shown in Appendix E.5, Table E Table 4 . As shown in Table 4 , the VG models from the method of moments do not fit the sample data distribution. The P values is less than 5% and the null hypothesis H 0 can not be accepted at some extent. The CLM does not fit the sample data distribution at 5%. Regarding the maximum likelihood method, the SVG2 model has P values = 9.079%, which is high than the classical threshold 5%. Therefore, we can not reject the null hypothesis that SVG2 model fits the empirical distribution. The histogram of the daily SPY ETF return was compared to some models in Table 4 as shown in Fig 10. Compare to the CLM in Fig 10b, 

Proof (2.10) and (3.5 )

Proof (3.6) :

Fourier of the step function

The integral is transformed into convolution of f and u

Using the convolution theorem 

APPENDIX E.5. SVG2 CUMULATIVE FUNCTION (F ) VERSUS EMPIRICAL FUNCTION (F n ), (1) = |F(x j ) −F n (x j−1 )|,(2) = |F(x j ) −F n (x j )|, n = 2755 x j n { j}Fn(x j)F (x j) (1) (2) x j n { j}Fn(x j)F (x j) (1) (2) x j n { j}Fn(x j)F (x j)(1)

Empirical properties of asset returns: stylized facts and statistical issues

Lévy processes in finance: pricing financial derivatives

Option pricing in a dynamic variance-gamma model. Derivatives eJournal

A subordinated stochastic process model with finite variance for speculative prices

Subordinated market index models: A comparison

The variance gamma (vg) model for share market returns

Chebyshev polynomial approximations and characteristic function estimation

The variance gamma process and option pricing

Fitting the variance-gamma model to financial data

A fast method for the numerical evaluation of continuous fourier and laplace transforms

The fractional fourier transform and applications. SIAM review

Kurtosis as peakedness

Probability theory. universitext

An introduction to probability theory and its applications

The origin of infinitely divisible distributions: from de finetti's problem to levy-khintchine formula

The advanced theory of statistics. The advanced theory of statistics

A modified compound poisson process with normal compounding

Sur l'estimateur du maximum de vraisemblance (emv

Maximum likelihood estimation based on newtonraphson iteration for the bivariate random effects model in test accuracy meta-analysis

Applied statistical inference

The variance gamma scaled self-decomposable process in actuarial modelling

The variance gamma (vg) model with long range dependence

Fitting the variance-gamma model: A goodness-of-fit check for emerging markets

statistics for business and economics"(book review)

Unifying political methodology: The likelihood theory of statistical inference

Evaluating kolmogorov's distribution

The kolmogorov-smirnov test for goodness of fit

The efficiency of multiple imputation and maximum likelihood methods for estimating missing values

Computing the kolmogorov-smirnov distribution when the underlying cdf is purely discrete, mixed, or continuous

KSgeneral-package: Computing P-Values of the K-S Test for (Dis)Continuous Null Distribution. R Foundation for Statistical Computing

Rachunek prawdopodobieństwa i statystyka matematyczna w zadaniach

Nig distribution in modelling stock returns with assumption about stochastic volatility

Sis epidemic model: Birth-and-death markov chain approach

Stochastic and Renewal Methods Applied to Epidemic Models

Epidemic dynamics and adaptive vaccination strategy: Renewal equation approach

Age-structured epidemic with adaptive vaccination strategy: Scalar-renewal equation approach