Pearson (1985) gave a four-parameter system of probability density functions and fitted the parameters by the method of moments (MME). Tukey (1960) proposed one-parameter lambda distribution. Tukey’s lambda was generalized, for the purpose of generating random variables for Monte Carlo simulation studies, to the four-parameter generalized lambda distribution (GLD), by Ramberg and Schmeiser (1972) and Ramberg and Schmeiser (1974). Ramberg et al. (1979) developed a four-parameter model together with the necessary tables for fitting a wide variety of curves. Since the early 1970s, the GLD has been applied in many fields of endeavor with continuous probability density functions. The generalized lambda distribution family with four parameters λ₁, λ₂, λ₃, and λ₄, denoted by GLD(λ₁, λ₂, λ₃, λ₄), has the density function,

where Q(y) is the percentile function defined as

where 0 ≤ y ≤ 1. The λ₁ and λ₂ are the location and scale parameter respectively. Further, λ₃ and λ₄ determine the skewness and kurtosis. We note there that not all choices of λ₁, λ₂, λ₃, λ₄ lead to a valid distribution, as described in the following theorem.

Theorem 1

The GLD(λ₁, λ₂, λ₃, λ₄) specifies a valid distribution if and only if

has the same sign for all y in [0, 1], as long as, λ₂, takes that sign. In particular, the GLD(λ₁, λ₂, λ₃, λ₄) specifies a valid distribution if λ₂, λ₃, λ₄all have the same sign.

More details for GLD family refer to Karian and Dudewicz (2000). This paper is organized as follows. In Section 2, we develop the change point detection procedure for the GLD model based on the MIC and provides preliminaries on how to construct a confidence curve and confidence set for the specified level. The simulation results are presented in Section 3. Our new method for the goodness-of-fit of the asymptotic null distribution of S_n is discussed in Section 4. Two real data applications are provided in Section 5. The summary of the results and some discussion are given in Section 6.

2.1. Modified information approach

In this section, we use the modified information approach (MIC) to detect changes in a GLD model. Chen et al. (2006) proposed the MIC as a the modification of the BIC approach by emphasizing the model complexity in the context of change point problems to include the contribution of the location of the change point.

In general, we consider multiple changes in the data set. Vostrikova (1981) proposed the binary segmentation method which could detect multiple structural changes recursively at most one change point at each step. In the first step, this method allows us to scan the whole data set by testing the null hypothesis of no change versus the alternative hypothesis of having one change. Once the first change (if there is any) has been located, the data is divided to two subsequences which are before the change point and after the change point. Then the second step is to repeat the same scanning procedure in Step 1 to these two subsequences respectively by assuming at most one change in each subsequence. Such a process will be repeated until there will be no further subsequences having change points. By doing so, we can find all the possible changes as well as estimate their locations. She also showed that the binary segmentation procedure is consistent. In particular, the multiple change point problem may be viewed as a single problem along with the binary segmentation procedure. Therefore, in this paper, we establish the testing procedure for a single change point detection method, however, it can be easily generalised to multiple change point problem, if needed.

Let X₁, . . . , X_n be a sequence of independent random variables from a GLD model with parameters,

We are interested in testing the changes in all parameters simultaneously. Thus, we would like to test the hypotheses are,

versus,

where 1 ≤ k < n. Then we can define the modified information criteria MIC(n) under H₀ and MIC(k) under H₁ respectively according to Chen et al. (2006) as follows,

where λ̂_i s, i = 1, 2, 3, 4 are MLEs under H₀, and λ^i(1),λ^i(n), i = 1, 2, 3, 4 are MLEs under H₁ which are estimated by R package GLDEX (Su, 2016). Let k be the possible change location in the range of 1 ≤ k < n. Then we accept H₀ if,

which indicates there is no change point in the data set, and we reject H₀ if

which indicates that there exists a change point in the data set. Consequently we can estimate the change point location k̂ by

Further, we define the test statistic S_n based on MIC(n) and MIC(k) to test the null hypothesis of no change versus the alternative hypothesis of one change as follows,

We reject the null hypothesis for a sufficiently large value of S _n. The standardization term 4 log(n) removes the constant term in the difference of MIC(n) and MIC(k).

2.2. Confidence distribution, profile log-likelihood and deviance functions

Confidence distributions (CD) are distribution estimates to be interpreted as distributions of epistemic probabilities. The concept of a CD is similar to a point estimator and it can be referred to as a sample-dependent distribution that can represent confidence intervals of all levels for a parameter of interest. A formal definition of CD can be found in Schweder and Hjort (2002). Furthermore, Schweder and Hjort (2016) systematically studied the theoretical properties of the CD. A detailed analysis of recent developments of CD has been given by Xie and Singh (2013). More applications of the CD can be found in the literature, including bootstrap distributions, p-value functions, normalized likelihood functions, and Bayesian posteriors, among others, Schweder and Hjort (2002), Singhet al. (2005), Singhet al. (2007), Singh and Xie (2012), and Shen et al. (2018).

The CD for change point analysis has been investigated by Cunen et al. (2018) and they construct confidence curves for change locations using the log-likelihood approach. Ratnasingam and Ning (2020) studied the change point detection procedure based on the CD combining with the MIC to construct the confidence set for the change estimate for a skew normal change point model. They also investigated the confidence distribution for detecting and estimating changes in a three-parameter Weibull distribution (Ratnasingam and Ning, 2021). In this paper, we study the CD-based detection procedure along with MIC for a GLD change point model. Next, we describe a procedure to construct a confidence curve for a four-parameter GLD change point model.

Suppose X₁, . . . , X_k be a sequence of independent random variables with the density function f (x,Θ_L) and X_k+1, . . . , X_n coming from the population with the density function f (x,Θ_R). Now the log-likelihood function is defined as

where Θ_L and Θ_R are the parameter space of the pre-change and post-change distributions respectively. By maximizing the log-likelihood function above of a given k, we can obtain the profile log-likelihood function as follows.

where Θ̂_L and Θ̂_R are MLEs of Θ_L and Θ_R respectively. The estimated change point location k̂ corresponds to the maxk(ℓprof(k)). The deviance function is defined as

where x = (x₁, x₂, . . . , x_n). The confidence curve for k based on the deviance function can be obtained through simulation.

where the cc(k, x_obs) < α under the true value of k. By simulation, we compute

for a large number of B of simulated copies of dataset x_obs. For each possible value of k, we simulated data xj*, j = 1, . . . , B from f (x,Θ_L) and f (x,Θ_R) to the left and right side of k, respectively. Furthermore, the change point location is estimated by (2.7). More precisely, our approach depends on the location of the change point. This approach is different from the method used in Cunen et al. (2018). For more details, we refer the readers to Cunen et al. (2018) and Ratnasingam and Ning (2020).

In this section, due to the difficulty in deriving its the analytic properties, we use simulations to investigate the critical values and the performance of the proposed test statistic S_n .

3.1. Critical values

We now describe how to obtain empirical critical values for our test statistic S _n, and for T_n, the test statistics proposed by Ning and Gupta (2009), who considered the classical Bayesian information criterion (BIC) to detect multiple changes in a GLD model, and then compare the performances their performances.

where BIC(n) under H₀ and BIC(k) under H₁ are given by

The major difference between S_n and T_n is that the penalty term in the MIC(k) incorporates the contribution of the location of the change point k associated with the complexity of the model, which is not accounted for by the BIC(k). In order to make a fair power comparison between S_n and T_n, we simulate the critical values for both test statistics under the same null distributions with the same sample sizes for given significance levels.

There are two general approaches available to compute critical values. They are simulation-and bootstrap-based approaches. The simulation-based approach requires the estimation of the test statistics values S_n and T_n under the null hypothesis of a certain number of repetitions, and critical values are equal to the percentiles of the sorted values of the test statistics from the simulations. The second method uses the bootstrap to obtain the asymptotic critical values for the test statistics. In this method, a certain number of samples are drawn from a null distribution with replacement over and over again from a null distribution. These are called bootstrap samples. The critical values for a given significance level correspond to the percentiles of sorted test statistics values.

When using the bootstrap method to obtain simulated critical values of a test statistic, we need to ensure that the bootstrap samples are re-sampled from data under the null distribution. In the simulation-based approach, however, the distribution under the null hypothesis has been determined before re-sampling. Therefore, it is known to satisfy H₀ which can be used to generate a sample. Thus, in simulations, both approaches will obtain similar critical values. However, for real data, it would be an issue for the bootstrap method since we do not know whether the data satisfies H₀ or H₁. Therefore, we can’t perform re-sampling directly on the data. The following strategy will be taken. We first assume the data satisfying H₀, which indicates it should be fitted by a single GLD, GLD₀ = GLD(λ̂₁,λ̂₂,λ̂₃,λ̂₄), where λ̂_i, i = 1, 2, 3, 4 can be obtained by R package GLDEX. Then we generate a random sample based on GLD₀ denoted by x₁, x₂, . . . , x_n. Bootstrap samples are drawn from this generated sample with replacement, denoted by y1(i),y2(i),…,yn(i), i = 1, 2, . . . , B. For each bootstrap sample, we calculate S_n denoted by Sn(i), i = 1, 2, . . . , B. Thus, the p-value can be approximated as follows

where I(·) is the indicator function and Sn(*) is the value of S_n calculated from the original real data.

3.2. Critical values of S_n and T_n

In our simulation study, we set up the null distribution to be GLD(2, 1, 0.19, 0.19) and choose sample sizes n = 50, 60, 80, 100, 150 with significance levels α = 0.01, 0.05, 0.1. We obtain the empirical critical values for S_n (2.1) and T_n (3.1) through the simulations as follows,

• Step 1: We generate data with various sample sizes n = 50, 60, 80, 100, 150 from GLD(2, 1, 0.19, 0.19).

• Step 2: For each generated sample, we calculate S_n and T_n respectively.

• Step 3: We repeat the above steps M = 1,000 times. Then the corresponding percentiles of these S_n and T_n values are the critical values at give significance level α = 0.01, 0.05, and 0.1.

The empirical critical values are provided in Table 1.

We should note here that, for the real data application, we should follow the bootstrap method proposed in Section 3.1 to calculate the p-value since whether the true distribution of the data satisfies H₀ or H₁ is unknown.

In this subsection, we provide results of the simulation study for the coverage probability, confidence sets, and consistency of the change point estimator. Three different sample sizes n = {50, 100, 150} are considered and each with various change point locations. For instance, sample size n = 50, the change point positions are set to k = 10, 15, and 25. The pre-change data are obtained from GLD(2, 1, 0.19, 0.19). One method is considered better than other methods if it produces thinner confident sets and retaining the right coverage for a given test level.

3.3. Power comparison

In this section, we conducted simulations under different scenarios to investigate the performance of test procedures based S_n and T_n in terms of power. We set up the distribution before change to be GLD(2, 1, 0.19, 0.19), and the distribution after the change to be GLD(λ1(n),λ2(n),λ3(n),λ4(n)) with various values listed in Table 2.

Further, various sample sizes n = {50, 100, 150} and different change point locations have been considered under each sample sizes. We consider the change point locations where τ positions {10, 15, 25} for sample size n = 50, {15, 25, 50} for sample size n = 100, and {35, 50, 75} for sample size n = 150. Note that the cases {40, 35, 25} for sample size n = 50, and {85, 75, 50} n = 100, and {155, 100, 75} for sample size n = 150 are fully symmetrical and thus they have the same power, coverage probability, and average size. Our test statistics (2.1) and (3.1) correspond to the MIC and BIC procedures respectively. The critical values at a given significance level are obtained through the simulations proposed in Section 3.2. The simulations results are recorded in Table 2.

Regardless of the method used for the power calculations, the power of the test increases as the sample size becomes larger. Moreover, it appears that the MIC based method gives larger power compared to the BIC based method. This may be due to that MIC method depends on the location of the change point. We also notice that the power of the test increases as the differences between the parameters increases.

3.4. Coverage probability, confidence sets & consistency of the estimator k̂ comparison

First, we examine the coverage probabilities when the method has the exact right coverage for the specified level. We compare our method, the MIC based method, with the log-likelihood based approach proposed in Cunen et al. (2018). The simulation results are listed in Table 3 for when α = {0.50, 0.90, 0.95, 0.99}. The performance of the MIC based method outperforms the log-likelihood based method in all cases. However, both methods provide slightly over coverage for a level α = 0.5 and this is more apparent as the sample size and the difference among the parameters increase. For example, sample size n = 100, the change point location k = 15 with the post-change distribution GLD(2.5, 1.5, 0.69, 0.69) and test level α = 0.90 the coverage level based on MIC method is 0.79 as compared to the log-likelihood based method which gives only 0.74.

Next, we compute the average sizes of confidence sets for MIC and log-likelihood based methods. The results are summarized in Table 4. It can be seen that the MIC based method gives smaller confidence sets compared to the log-likelihood based approach proposed in Cunen et al. (2018). We observe that the size of the confidence set becomes smaller when the difference between the parameters increases. For example, for n = 50, test level α = 0.5, and the post-change distribution GLD(2.5, 1.5, 0.69, 0.69) the average length of the confidence set based on MIC is equal to 7.342, however, the log-likelihood based approach provides slightly large confidence set which average length is 7.556.

We also investigate the consistency of the estimator k̂ of the actual change location k through a numerical study. The difference between the estimated location and the actual change point location is set to δ. Further, we compute the bias and mean squared error (MSE) as well. Table 5 summarizes the simulation results below. It can be concluded from the simulations that the P(|k̂_MIC – k| ≤ δ) is greater than P(|k̂_BIC – k| ≤ δ). However, the difference between these two probabilities decreases as the sample size n increases.

According to Chen et al. (2006), under certain conditions, the asymptotic null distribution of the proposed test statistic S_n is the χ²-distribution with d degrees of freedom, where d is the number of parameters in the null model. In this section, we develop a Goodness-of-fit test to evaluate the asymptotic null distribution of S_n using the energy distance.

Energy distance, as described in Székely (2000), is defined to be the statistical distance between probability distributions. The associated statistics, named energy statistics, are the function of energy distances. The concept is motivated by Newton’s gravitational potential energy which is a function of the distance between two objects. Thus the idea of energy statistics is to consider statistical observations as heavenly bodies governed by a statistical potential energy, which is zero if and only if an underlying statistical null hypothesis is true. Székely and Rizzo (2013) defined the energy distance ℰ(X, Y) between two independent d-dimensional random variables X and Y is computed by

provided E|X|, E|Y| < ∞. Here X′ is and i.i.d copy of X and Y′ is an i.i.d copy of Y. ℰ(X, Y) ≥ 0 and ℰ(X, Y) = 0 if and only if X=dY. There have been numerous studies based on the energy distance. For example, Székely and Rizzo (2005) proposed a test based on the energy distance for multivariate normality. Rizzo (2009), Yang (2012), and Rizzo and Haman (2016) considered one-sample goodness-of-fit tests for Pareto distributions, Univariate stable distributions and asymmetric Laplace distributions based on the energy distance. Székely and Rizzo (2004) constructed an energy-distance-based test for testing equality of distributions under high dimensional settings. Baringhaus and Franz (2004) also proposed a new multivariate two-sample test based on energy distance. Kim et al. (2009) and Matteson and James (2014) studied change point problems incorporating the energy distance. In the univariate case, d = 1, the energy distance (4.1) becomes,

Then the one-sample energy statistic for the goodness-of-fit test based on the energy distance is given by the following definition.

Definition 1

Let X₁, . . . , X_n be a random sample from a univariate population with distribution F and let x₁, . . . , x_n be the observed values of the random variables in the sample. Then the single sample energy statistic for testing the hypotheses H₀ : F = F₀vs H₁ : F ≠ F₀is

where X and X′ are independent and identically distributed with distribution F₀and the expectations are taken with respect to the null distribution F₀.

As we mentioned earlier, according to Chen et al. (2006), the S_n defined in (2.1) follows the χ42 distribution. Thus, our testing hypothesis becomes H0:F=χ42 vs H1:F≠χ42. Therefore, we only need to derive the energy statistics formula for the χ42 distribution.

Theorem 2

Let Y~χ42, then for any fixed x ∈ ℝ

Theorem 3

Let X and X′ be independent identically distributed random variables. Then

where FX-1is the inverse CDF of X, n is the number of equally sized sub-intervals of [0, 1] and y_i is chosen from the i^th sub-interval.

The proof of Theorem 3 is similar to Rizzo and Haman (2016). Thus, details are omitted to conserve space. According to Rizzo (2002), the last term of the (4.2) can be written as,

where y₍₁₎ ≤ y₍₁₎ ≤ ·· · ≤ y_(n) is the ordered sample. Hence, the one sample energy statistic based on Definition 1 can be re-written as follows.

where n^* is the number of equally sized sub-intervals of [0, 1] and y_i is chosen to be in the i^th subinterval. The one sample energy goodness-of-fit test statistic procedure is described below,

• Generate a data x₁, x₂, . . . , x_n from the χ42 distribution.

• Compute the energy statistics of the data x₁, x₂, . . . , x_n using the formula (4).

• Repeat Steps 1 & 2 for 5000 times and obtain Ψn(1),Ψn(2),…,Ψn(5000).

• The critical value can be obtained by finding a 95% quantile of the energy statistics.

• Simulate S_n values using the equation (2.1) and denote them as Sn(1),Sn(2),…,Sn(B).

• Compute the energy statistic of the data Sn(1),Sn(2),…,Sn(B).

• Compare the energy statistic in Step 6 with the critical value found in Step 4, if the critical value exceeds the energy statistic, we conclude that S_n values come from χ42 distribution otherwise not.

We follow the above procedure to conduct the one sample energy goodness-of-fit test statistic. The critical value at 5% significance is 9.5037. The energy statistic for the S_n data is 11.60781 so we reject the null hypothesis at 5% significance level. This suggests that the asymptotic null distribution of S_n does not follow a χ42 distribution. Below we construct the Chi-Square Q-Q plot for S_n test statistic values.

Figure 1 shows that there is a significant deviation from the reference line. This confirms our previous results, that the asymptotic null distribution of S_n does not follow a χ42 distribution. Thus, the MIC statistic derived from the GLD change point model does not comply with the asymptotic properties established in Chen et al. (2006).

In this section, the proposed method is applied to analyze two stock market returns data from the Brazilian and Chilean markets. These data sets were previously used in the literature Ngunkeng and Ning (2014), Ratnasingam and Ning (2020). We assume that changes occur simultaneously across all four parameters. The stock return ratio is obtained through the following transformation.

5.1. Brazilian market return ratio data

In order to identify multiple changes and to create the appropriate confidence sets, we use the proposal method along with binary segmentation procedure. First, MIC(262) = −807.3228. Then the min_1≤k<262 MIC(k) = MIC(87) = −824.2035. Thus the estimated change point locationk̂ = 87, which corresponds to the 88^th location in the data set. The MLEs of the pre-change parameters and post-change parameters are (0.0089, 64.4382, 0.0156, −0.0085) and (0.0010, 32.8107, −0.0777, −0.0389) respectively and the 95% confidence set for the change point estimate is {87, 88, 89}. We then split the data sets into two subsets which are below k(≤ 87) and above k(> 87), then the proposed method is applied recursively in each subset in order to detect all changes in the data sets. This iterative process stops until there are no further changes detected. In particular, all the change points are obtained from our procedure {58, 88, 144, 240, 254}. As opposed to a GLD model used in Ratnasingam and Ning (2020), we identify one additional change point at 253. When compared with Ngunkeng and Ning (2014), we locate an additional change at 58th data point. Figures 2 and 3 show the confidence curves for all change-point estimates and the 95% confidence sets are marked red dashed lines. All change points are graphed in Figure 4.

5.2. Chilean market return ratio data

First, we compute MIC(n) = −1032.118 and min_1≤k<262 MIC(k) = MIC(111) = −1047.216. Thus the change point estimate k̂ is 111. The corresponding change point in the data is 112. Then the MLEs of the pre-change and post-change parameters are (−0.0042, 58.20230.3994, 0.1338) and (0.0012, 48.6707, −0.0726, −0.0061) respectively. Further, the 95% confidence set for the change point is {110, 111, 112, 113, 114}. Next we apply the binary segmentation procedure to detect all the changes in the data set. They are {98, 112, 170, 181}. When compared with Ratnasingam and Ning (2020), we found an additional change point in the Chilean market data and only one change point. Further, we found an additional change point 98 when compared to Ngunkeng and Ning (2014). The confidence curves for all change point estimates and the 95% confidence sets are marked red dashed lines in Figure 6.

In this paper, we propose a change point detection procedure for a GLD model based on the modified information criterion. In order to use as much information about the change point location, the proposed procedure takes into account the effect in terms of model complexity in regards to the location of the change point. We provide confidence sets for the change point location for a specified level α. We also obtain empirical critical values of the test statistics have been found. Simulations conducted with different sample sizes and various change point locations show that our method performs well in terms of a larger power, smaller confidence set, and smaller MSE when compared to other methods. Furthermore, we introduce a new goodness-of-fit test based on energy distance to determine the asymptotic null distribution of S _n. The use and the advantage of the proposed method are illustrated via two stock market data sets.