A New Goodness-of-Fit Test for Azzalini's Skew-t Distribution Based on the Energy Distance Framework with Applications

1. Introduction

Skewed and heavy-tailed data are prevalent in various applied domains, including econometrics, environmental science, and risk analysis. Income distributions, housing prices, and insurance claims often display asymmetry and excess kurtosis; the tail behavior encodes meaningful extremes, such as financial losses or contaminant spikes, that should not be dismissed as outliers; see, for example, [1,2,3] and Ahmad et al. [4].

A natural starting point for modeling asymmetry is the skew-normal distribution introduced by Azzalini [5] (see also [6]). By augmenting the normal distribution with a skewness parameter, the skew-normal preserves analytical tractability while allowing controlled departures from symmetry. However, because the skew-normal remains light-tailed, it is ill-suited to settings where leptokurtosis is intrinsic to the data-generating process.

To address heavy tails alongside asymmetry, the Skew-t family extends the t distribution by introducing a skewness parameter. An influential approach views Skew-t laws as scale mixtures of skew-normal variables [5,7], a perspective adopted and elaborated by several authors, including Hasan et al. [8]. Related extensions, such as extended skewt (EST) and alternative parameterizations, further enhance the modeling flexibility (see, for example, [9]). Empirically, Skew-t models have seen broad application, from financial time series and risk assessment to environmental monitoring and robust analysis under truncation or censoring [10,11]. Noncentral variants expand the toolkit and have been studied in detail, with their properties and use cases documented in Hasan et al. [8], Hasan [12].

Despite this progress, the assessment of the Skew-t model remains essential. In practice, misspecification of tail weight or asymmetry can distort inference on extremes, dependence, and risk. This motivates rigorous goodness-of-fit (GoF) procedures tailored to Skew-t families—methods that can diagnostically test statistical adequacy against alternatives that differ in tail behavior, asymmetry, or both.

In this paper, we develop an energy-based goodness of fit test for Azzalini’s standard Skew-$t$ distribution proposed in [6] and defined as a random variable with probability density function (pdf) that takes the following form:

$$f_{\nu ,\alpha }\left(x\right)=2g_{\nu }\left(x\right)G_{\nu +1}\left(\alpha x\sqrt{{\displaystyle \frac{\nu +1}{x^2+\nu }}}\right),$$

(1)

where g and G are the regular probability density function (PDF) and the cumulative distribution function (CDF) of the Student’s t distribution. Parameters $\alpha$ and $\nu$ are referred to as skewness and degrees of freedom, respectively. We will denote this density by $St(\nu ,\alpha )$. The special case for $\alpha =0$ reduces to the standard Student’s t distribution with $\nu$ degrees of freedom. The corresponding CDF of the skew t distribution is given by:

$$F_{\nu ,\alpha }\left(x\right)={\int }_{-\infty }^x{f}_{\nu ,\alpha }\left(t\right)dt,$$

(2)

which does not have a closed form. Its first moment is given as

$$E\left(X\right)={\displaystyle \frac{\Gamma \left(\frac{1}{2}(\nu -1)\right)}{\Gamma \left(\frac{\nu }{2}\right)}}\delta \sqrt{{\displaystyle \frac{\nu }{\pi }}},$$

(3)

where $\nu >1$, $\delta ={\displaystyle \frac{\alpha }{\sqrt{1+{\alpha }^2}}}$. It should be noted that the PDF given in equation (Section 1) is referred to as the standard Azzalini’s Skew-t distribution. In the literature, a scale-location extension has been proposed. The reader is referred to [6,13] for more details. Figure 1 provides density estimates of the Skew-t distribution for varying degrees of freedom ($\nu$) when the skewness parameter $\alpha =1$ is held constant.

There is limited research on the goodness-of-fit test for Azzalini’s Skew-t distributions except those based on empirical distribution functions (EDFs) such as Kolmogorov-Sminorv, Cramer-von-Mises, etc., discussed and modified in [14,15]. Recently, [16] proposed a goodness-of-fit test for Azzalini’s Skew-t distribution based on the correlation coefficient (r). Most studies related to Skew-t distributions focus on properties, extensions, generalizations, change point analysis, and applications; see for example, [12,17,18], and [19], among many others. We therefore develop the goodness-of-fit test for the Skew-t distribution using these well-established properties to achieve desirable powers for any given sample size and the chosen parameters $\alpha$ and $\nu$.

In this article, we propose a new one-sample (univariate) energy goodness-of-fit test based on energy statistics proposed by [20] and [21]. In the most recent work, Opperman and Ning [22,23], and [24] proposed goodness-of-fit tests based energy statistics for the Skew-normal, Inverse Gaussian, and Lindley distributions, respectively. For a given sequence of independent random variables of size n and with a cdf G, the test statistic based on energy statistics will reject the null hypothesis that $F=G$ for large values of the test statistic. If the null distribution, F, and the given data come from the same underlying distribution G, then the values of the test statistic are expected to be smaller. Furthermore, there have been numerous studies involving energy statistics such as testing for multivariate normality ([25,26]), testing for equality of distributions ([27,28]), one-sample goodness-of-fit tests ([22,24,29,30]), change point analysis ([31,32,33,34]), among many others.

The energy distance is defined as a statistical distance between the distributions of random vectors which characterizes equality of distributions, see for example [29,35] and [21]. The concept of energy statistics described by Sźekely [21] is based on the notion of Newton’s gravitational potential energy, which is a function of the distance between two bodies. The idea of energy statistics, therefore, 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, see for example [29,36].

Definition 1. 

Sźekely and Rizzo [35] defined the energy distance between distributions of two independent and univariate random samples X and Y with finite expectations as follows

$$\mathcal{E}(\mathbf{X},\mathbf{Y})=2E|X-Y|-E|X-X^{\prime }|-E|Y-Y^{\prime }|\ge 0,$$

(4)

$X\stackrel{d}{=}{X}^{\prime }$, $Y\stackrel{d}{=}{Y}^{\prime }$ and equality holds if and only if X and Y are identically distributed.

1.1. Contribution and Novelty

In this work, we propose a procedure that is more superior for goodness-of-fit test for the Skew-t distribution based on energy distance statistics [35] and [13] definition of the Skew-t distribution. Unlike the proposed method based on energy statistics, many existing methods depend on the distribution function of random variables. Energy statistic-based tests have been shown to be typically more powerful against general alternatives than corresponding tests based on classical statistics (non-energy type) such as Kormogorov-Smirnov, Correlation, Anderson Darling, and Cramer-Von-Mises. In addition, energy statistic-based tests have an invariance property with respect to any distance-preserving transformation of the data set, see for example [29,33,36]. This makes it possible to extend studies involving energy statistics to multivariate and high dimensional settings.

In Section 2, we propose a test procedure based on energy statistics for the goodness-of-fit of Azzalini’s Skew-t distribution and discuss its theoretical properties. We perform various simulations in Section 3 to compare the approach with other existing methods. In Section 4, we apply our method to two real-life datasets. The conclusion is provided in Section 5.

2. Proposed Energy-based Goodness of Fit Test

We propose a one-sample univariate goodness-of-fit test based on the energy statistics proposed by [26,35] for the Skew-t distribution. The null hypothesis is that the data X follow the null distribution $F_0,$ which is a Skew-t distribution, against the alternative that the Skew-t distribution is a poor fit for the data.

Definition 2. 

Let $X_1,...X_n$ be a random sample from a univariate population with distribution F and let $x_1,...x_n$ be the observed values of the random variables in the sample. Then, the one-sample energy statistic goodness-of-fit test for testing the hypothesis $H_0:F=F_0$ vs $H_a:F\ne F_0$ is defined as follows.

$$n{\mathcal{E}}_n=n\left\{{\displaystyle \frac{2}{n}}\sum _{i=1}^{n}E|x_i-X|-E|X-X^{\prime }|-{\displaystyle \frac{1}{n^2}}\sum _{i=1}^n\sum _{j=1}^n|x_i-x_j|\right\},$$

(5)

where X and $X^{\prime }$ are independent and identically distributed variables with distribution $F_0$, and the expectations are taken with respect to the null distribution $F_0$.

The null hypothesis, $F=F_0,$ is rejected for large values of the test statistic $n{\mathcal{E}}_n$. Under the null hypothesis, the limiting distribution of $n{\mathcal{E}}_n$ is a quadratic quantity of the form ${\sum }_{j=1}^{\infty }{\lambda }_j{Z}_j^2$ such that $Z_j,j=1,2,\dots ,$ are i.i.d. standard normal random variables and ${\lambda }_j$ are nonnegative constants that depend on the null distribution. Thus, the goodness-of-fit test can be implemented by finding the constants ${\lambda }_j$. In practice, this could be difficult, and we therefore resort to the use of empirical critical values of $n{\mathcal{E}}_n$ so that $P(n{\mathcal{E}}_n>C_{\alpha })=\alpha .$ This fact is guaranteed since the test based on $n{\mathcal{E}}_n$ is a consistent goodness-of-fit test, see for example, [26] and [25].

The goodness of fit test statistic based on energy statistics is dependent on the derivation of the expected values of $|x-X|$ and $|X-X^{\prime }|$ where X and $X^{\prime }$ are independent and identically distributed random variables from the null distribution $F_0.$

Proposition 1. 

Let $X\sim f\left(x\right)$ where $f\left(x\right)$ is defined in Eq. (2). Then, for any fixed $x\in \mathbb{R},$

$$E|x-X|={\int }_{\mathbb{R}}|x-t|f\left(t\right)dt=2xF\left(x\right)-x+{\displaystyle \frac{\Gamma \left(\frac{1}{2}(\nu -1)\right)}{\Gamma \left(\frac{\nu }{2}\right)}}\delta \sqrt{{\displaystyle \frac{\nu }{\pi }}}-2{\int }_{-\infty }^{x}tf\left(t\right)dt,$$

(6)

where $F(.)$ is the CDF of the Skew-t distribution.

Proof 

(Proof of Proposition 1). Let $X\sim f\left(x\right)$. Then, for every fixed real number x we have

$$\begin{array}{cc}\hfill E|x-X|& ={\int }_{\mathbb{R}}|x-t|f_X\left(t\right)dt\hfill \\ \hfill & ={\int }_{-\infty }^x(x-t)f_X\left(t\right)dt+{\int }_x^{\infty }(t-x)f_X\left(t\right)dt\hfill \\ \hfill & =x{\int }_{-\infty }^x{f}_X\left(t\right)dt-{\int }_{-\infty }^{x}t{f}_X\left(t\right)dt+{\int }_x^{\infty }t{f}_X\left(t\right)dt-x{\int }_x^{\infty }{f}_X\left(t\right)dt\hfill \\ \hfill & =x{F}_X\left(x\right)-{\int }_{-\infty }^{x}t{f}_X\left(t\right)dt+\left({\int }_{-\infty }^{\infty }t{f}_X\left(t\right)dt-{\int }_{-\infty }^{x}t{f}_X\left(t\right)dt\right)-x(1-F_X\left(x\right))\hfill \\ \hfill & =2xF\left(x\right)-x+E\left(X\right)-2{\int }_{-\infty }^{x}t{f}_X\left(t\right)dt\hfill \\ \hfill & =2xF\left(x\right)-x+{\displaystyle \frac{\Gamma \left(\frac{1}{2}(\nu -1)\right)}{\Gamma \left(\frac{\nu }{2}\right)}}\delta \sqrt{{\displaystyle \frac{\nu }{\pi }}}-2{\int }_{-\infty }^{x}t{f}_X\left(t\right)dt,\hfill \end{array}$$

where the integral term can be numerically evaluated in R using the command dst() available in Azzalini sn package.    □

Sometimes the derivation of the second term in Eq. (5) may not be analytically feasible and, therefore, we can use the following approximation as suggested in [22].

Proposition 2.(Quantile-based Approximation) Let X and $X^{\prime }$ be independent and identically distributed random variables with a well-defined cumulative distribution function, $F\left(x\right)$. Given that the quantile or inverse CDF function of X exists, we have the following.

$$E|X-X^{\prime }|={\displaystyle \frac{4}{m}}\sum _{i=1}^m{x}_i{F}^{-1}\left(x_i\right)-{\displaystyle \frac{2}{m}}\sum _{i=1}^m{F}^{-1}\left(x_i\right),$$

(7)

where m is the number of equally sized sub-intervals of $[0,1]$ and $x_i$ is chosen from the $i^{th}$ sub-interval. It is worth noting that the Proposition 2 applies for all distribution functions.

Proof. 

This proof is provided by Opperman and Ning [22].    □

3. Simulations and Results

This section presents a simulation-based assessment of the proposed goodness-of-fit (GoF) test for Azzalini’s Skew-t distribution, built upon the energy distance framework. We evaluate the test’s performance across varying values of the skewness parameter ($\alpha$) and degrees of freedom ($\nu$). Simulations are conducted under both the null hypothesis and multiple alternative distributions, and Type I error and empirical power are computed.

In this simulation study, we calculate two key quantities. The size of the test (type I error) and the empirical power (1-type II error). We used R (R version 4.4.2) and the sn package https://cran.r-project.org/web/packages/sn/refman/sn.html to carry out the simulations. Due to the complexity of repeatedly evaluating the Skew-t CDF and performing numerical integration, all simulations were parallelized using the foreach and doParallel R packages. Simulations were distributed across available cores, significantly accelerating the evaluation of the test under various scenarios.

3.1. The Univariate Energy Test Statistic

The univariate energy test statistic is derived from the expected pairwise distances between observations and a reference Skew-t distribution. Let $F_{\theta }$ be the cumulative distribution function (CDF) of Azzalini’s Skew-t distribution with parameter vector $\theta =(\xi ,\omega ,\alpha ,\nu )$. The test statistic for a sample $\{x_1,\dots ,x_n\}$ is:

$$\begin{array}{cc}\hfill T_n& =n\{{\displaystyle \frac{2}{n}}\sum _{i=1}^n\left[2x_{i}F\left(x_i\right)-x_i+{\displaystyle \frac{\Gamma \left(\frac{1}{2}(\nu -1)\right)}{\Gamma \left(\frac{\nu }{2}\right)}}\delta \sqrt{{\displaystyle \frac{\nu }{\pi }}}-2{\int }_{-\infty }^{x_i}t{f}_X\left(t\right)dt\right]\hfill \\ \hfill & -{\displaystyle \frac{4}{m}}\sum _{i=1}^m{x}_i{F}^{-1}\left(y_i\right)-{\displaystyle \frac{2}{m}}\sum _{i=1}^m{F}^{-1}\left(y_i\right)-{\displaystyle \frac{1}{n^2}}\sum _{i=1}^n\sum _{j=1}^n|x_i-x_j|\}.\hfill \end{array}$$

(8)

This expression is approximated using the following.

(1)

A term involving the integral of $y·f_{ST}(y;\alpha ,\nu )$ up to each data point.

(2)

A quantile-based approximation for the inter-distributional expectation $\mathbb{E}|Y-Y^{\prime }|$.

(3)

A linear-order statistic approximation for the within-sample expectation $\mathbb{E}|x_i-x_j|$.

According to Rizzo [20], the last term of the test statistic $n{\mathcal{E}}_n$ in Eq. (5) can be linearized in order to reduce the computational complexity of the test from $\mathcal{O}\left(n^2\right)$ to $\mathcal{O}(nlogn)$ which is useful during extensive simulations and applications. Let $x_{\left(1\right)},x_{\left(2\right)},\dots ,x_{\left(n\right)}$ be the ordered sample of the random sample $x_1,x_2,\dots ,x_n$. Then the linearization of the double sum of the test $n{\mathcal{E}}_n$ is given as

$$\sum _{i=1}^n\sum _{j=1}^n|x_i-x_j|=2\sum _{k=1}^n((2k-1)-n)x_{\left(k\right)}.$$

(9)

3.2. Critical Value Simulation Under Varying $\alpha$ and $\nu$

In this section, we conduct a simulation study to estimate the critical values for the energy-based goodness-of-fit test under Azzalini’s skew t distribution. The critical value simulations are summarized in Table 1. We selected degrees of freedom $df=5,10,30$ to account for different levels of heavy tail behavior and skewness coefficient $\alpha =1,0,-1$ to cover the three possible scenarios of right skew, symmetry, and left skew. We used 5,000 replicates in this study, with sample sizes $n=50,100,150,200$. Table 2 presents the results of Type I error simulations, illustrating that the empirical values are very close to the theoretical value of 0.05, regardless of the sample size and the parameter values chosen for each sample.

3.3. Type I Error Control

To explore the sensitivity of the test, we vary the skewness parameter $\alpha \in \{-1,0,1\}$ and degrees of freedom $\nu \in \{5,10,30\}$. For each combination of $(\alpha ,\nu )$, we:

(1)

Generate $B=5000$ samples of size $n=50$ from the Skew-t distribution with parameters $(0,1,\alpha ,\nu )$.

(2)

Estimate the parameters $\widehat{\theta }=(\widehat{\xi },\widehat{\omega },\widehat{\alpha },\widehat{\nu })$ using the maximum penalized likelihood method available in the sn package.

(3)

Standardize the data to standard Skew-t random variables and compute the energy goodness-of-fit statistic using the formula in equation (8).

(4)

Determine the empirical ${95}^{\mathrm{th}}$ percentile critical value under the null hypothesis.

The empirical Type I error is then evaluated across different parameter settings to confirm robustness

Table 2 shows that the proposed energy-based test maintained the nominal significance level in all examined configurations. For $n=50$, the empirical Type I error was close to 0.05, although slightly conservative in extreme skewness or low degrees of freedom. As the sample size increased, the test stabilized and consistently aligned with the theoretical level.

Figure 2 shows critical values with varying $\alpha$ for $\nu =5$ and sample size $n=50.$. For each selected level of significance ($LOS=0.01,0.025,0.05,0.10$), the critical values seem to stabilize as the skewness parameter $\alpha$ progresses further from 0. At $\alpha =0$, we see a dramatic increase in the simulated critical value. Note that Azzalini’s skew t distribution reduces to the classic Student t distribution when $\alpha =0$. In this special case, Azzalini’s skewness t might be an overfitting model for the data because the skewness parameter is no longer needed, and Student’s t distribution is a better fit for the data. If $\left|\alpha \right|>2$, the critical value stabilizes. This implies that for a fixed predetermined confidence level, the actual value of the skewness parameter $\alpha$ has little to no effect on the critical value as long as $\left|\alpha \right|>2$; in other words, the data exhibit noticeable skewness.

Figure 2 shows critical values with varying $\nu$ for $\alpha =1.0$ and sample size $n=50.$. When the degrees of freedom $\nu =1$, Azzalini’s skew t distribution reduces to Cauchy’s distribution, and this is why we see higher critical values for the goodness-of-fit test. As the degrees of freedom increase $\nu >10$, critical values for a predetermined significance level stabilize as depicted by the nearly horizontal lines in Figure 2.

3.4. Power Analysis under Various Alternatives

The simulation findings demonstrate the superior statistical power of the proposed energy-based goodness-of-fit,$n{\mathcal{E}}_n$, test statistic compared to several established goodness-of-fit tests, including Anderson-Darling (A-D), Cramér-von Mises (CvM), Watson, Kolmogorov-Smirnov (K-S), and Kuiper tests. The analysis was conducted across various alternative distributions with sample sizes (n) of 50, 100, 150, and 200.

For power studies, samples are drawn from alternative distributions that deviate from Azzalini’s Skew-t distribution family. These include:

• Chi-square: asymmetric, heavy-tailed.
• Standard Gusset’s t: symmetric, heavy-tailed.
• Exponential: lighter tail, positive skew.
• SHASH: The SHASH (Sinh-Arcsinh) distribution that was discussed in [37] is a highly flexible statistical distribution defined by four parameters that separately control the location, scale, skewness, and kurtosis of a variable.
• Generalized t (GT): to assess sensitivity to misspecification.
• Log-normal: asymmetric, heavy-tailed.

For each setting, the test statistic is calculated and compared to the critical value obtained under the corresponding null distribution with ($\alpha$, $\nu$). Empirical power is computed as the proportion of times the null is rejected across $B=5000$ repetitions. Table 3 summarizes the results of the power comparison simulations. For each case, the bold font indicates the highest power achieved at the combination of sample size and alternative distribution displayed in the row heading, and for each of the tests shown in the column heading.

3.5. Superior Performance of the $n{\mathcal{E}}_n$ Test

As illustrated in Table 3, across all scenarios, the energy-based test exhibited the highest power. This dominance is particularly evident in cases of heavy-tailed and skewed distributions.

• For the Log-Normal ($LN(0.5,1)$) distribution, the $n{\mathcal{E}}_n$ test achieved a power of 0.9701 at a sample size of just 50, far exceeding the next best test, Watson test, which had a power of 0.3988.
• Against the GT($\nu =3,\tau =1$) distribution, the $n{\mathcal{E}}_n$ test reached a power of 1.0000 for sample sizes of 150 and 200, indicating perfect detection in the simulation.
• For the Exponential (Exp(1)) distribution, the power of $n{\mathcal{E}}_n$ ranged from 0.8413 to 0.9848, consistently outperforming all other competitors.

The K−S and A−D tests underperformed in almost all scenarios, especially for skewed and heavy-tailed data. The correlation-based test was competitive but slightly less sensitive to departures in kurtosis. The energy test consistently ranked among the top-performing methods in all conditions.

3.6. Effect of Sample Size and Parameters

As anticipated, the power for all tests increased with the sample size. For instance, in the case of the SHASH distribution, the power of the $n{\mathcal{E}}_n$ test rose from 0.5734 for $n=50$ to 0.6749 for $n=200$. It should be noted that for the standard $t_{10}$ distribution, which is a special case of Azzalini’s t distribution when the skewness coefficient $\alpha =0$, all tests showed very low power. The highest power achieved was only 0.2790 by the $n{\mathcal{E}}_n$ test at $n=200$. This suggests that all tests face difficulties in distinguishing the two distributions, although the $n{\mathcal{E}}_n$ test still holds a relative advantage. The proposed energy-based test achieved the highest power under conditions of high skewness ($\alpha =5$) and low degrees of freedom ($\nu =3$). As $\nu$ increased, the power decreased slightly. This decrease in power is likely because the Skew t distribution converges to the Skew-normal family as the degrees of freedom, $\nu \to inf$, making it harder to distinguish Skew-normal from Skew-t distribution as the degrees of freedom increase.

3.7. Summary

In summary, the simulation results provide strong evidence that the energy-based, $n{\mathcal{E}}_n$, statistic is a more powerful goodness-of-fit test than the other methods considered in the study, especially for non-normal distributions with skewness or heavy tails. Overall, the energy-based GoF test for the Skew-t distribution demonstrates:

• Robust control of Type I error
• Superior power in a wide range of alternatives
• Flexibility for skewed, heavy-tailed, and multi-modal data

Its simplicity and effectiveness make it a promising tool for model validation in applied settings. Although the proposed method relies on a Monte Carlo approximation for expectation terms, the computational cost is manageable and can be parallelized. Use of linearized summation and pre-simulated reference distributions further improves efficiency.

4. Real Data Application

To evaluate the practical utility of the proposed test, we applied it to two real-world datasets. These two datasets are: body mass index (BMI) for 102 male Australian athletes (Cook and Weisberg [38]) and the daily rate of returns for Apple stock Macrotrends [39], which provides the historical price of the stock.

Since the underlying distributions of these datasets are not exactly known in advance, we use the bootstrap algorithm to determine whether or not they come from Azzalini’s Skew-t distribution. The bootstrap procedure is performed to approximate the associated p-value of the proposed test as given below.

(1)

Fit the real data $y_1,\dots ,y_n$ with an Azzalini’s Skew-t distribution in Eq. (Section 1), and obtain the maximum likelihood estimates (MLEs) of $\xi ,\omega ,\alpha$ and $\nu$ from the Azzalini’s sn package available in R.

(2)

Use the formula in Eq. (8) to calculate the energy goodness-of-fit statistic of the real data and denote it ${\mathbb{T}}_n^1$.

(3)

Simulate $x_1,\dots ,x_n$, a random sample of size $n,$ from the Azzalini Skew-t distribution with parameters specified as $\widehat{\xi },\widehat{\omega },\widehat{\alpha }$ and $\widehat{\nu }$ that were obtained in step 1.

(4)

Standardize the data to standard Skew-t distribution and compute the energy goodness-of-fit statistic for the simulated data using the formula in Eq. (8) and denote this value as ${\mathbb{T}}_n^{*1}$.

(5)

Repeat this process for B times and obtain B energy goodness-of-fit statistics and denote them by ${\mathbb{T}}_n^{*1},\dots ,{\mathbb{T}}_n^{*B}$.

(6)

The bootstrap p-value is therefore approximated as

$$\widehat{p}={\displaystyle \frac{1}{B}}\sum _{b=1}^{B}I\left({\mathbb{T}}_n^{*b}\ge {\mathbb{T}}_n^1\right),$$

where $I(·)$ is an indicator function that takes the value of one when ${\mathbb{T}}_n^{b*}\ge {\mathbb{T}}_n^1$ and zero otherwise.

In the application, the first dataset represents the body mass index (BMI) of male Australian athletes recorded in Table A1 and is available in the dr package in R, see for example [38]. This dataset was previously analyzed by Maghami and Bahrami [16] with the correlation coefficient as the measure of goodness-of-fit. The dataset exhibits moderate skewness and potential heavy tails, making it an ideal candidate for Skew-t modeling. The second dataset consists of the closing prices of Apple Inc. stock. This data was used as part of [19] analysis for detecting structural changes in the distribution via the MIC-based method. The dataset is obtained from [39], which provides historical stock price information. We limited the data to the range from January 31, 2019, to January 24, 2020 to avoid the heterogeneity introduced by change points. The stock price shows a general upward trend with some noticeable fluctuations. Given that stock prices exhibit time dependence, we transform the raw closing prices, $P_t$, into daily returns, defined as

$$R_t={\displaystyle \frac{P_{t+1}-P_t}{P_t}}.$$

This transformation allows us to analyze relative price changes rather than absolute levels. The resulting data consist of 248 daily return observations, which we suspect follow a Skew-t distribution. The complete transformed data is available in Table A2 in the Appendix.

4.1. Model Fitting

The histograms of the data in Figure 3 and Figure 4 reveal skewness; thus, distributions such as the Skew-t, Skew-normal, Skew-Cauchy, and Normal distributions can be good candidates for modeling these datasets. Parameter estimation is conducted using the maximum likelihood method. The results of the penalized maximum likelihood estimates (PMLE), log-likelihood (LogL), Akaike information criterion (AIC), and Schwarz information criterion (SIC) are reported in Table 4 and Table 5. In both datasets, the Skew-t distribution yielded the lowest AIC values, suggesting a better fit overall. Surprisingly, the skew-normal distribution gave a lower SIC value for the body mass index (BMI) dataset.

In the first dataset (BMI), our proposed test statistic is 1.7109, and its corresponding p-value is 0.5840, supporting the assertion that the BMI data follow the Skew-t distribution. Other tests considered in the study also supported the null hypothesis that body mass index (BMI) data can be modeled using the Skew-t distribution (Table 6). In the second data, the test statistic of the proposed procedure is 3.9635, and the corresponding p-value is given as 0.5473. These results support the fact that the data follow the Azzalini’s Skew-t distribution. Similar results are observed for the empirical distribution function (EDF) tests, as they suggest that the data can be modeled with the Skew-t distribution. These results are summarized in Table 6 and Table 7. Furthermore, density estimates and empirical cumulative functions suggest that the Azzalini Skew-t distribution fits the datasets adequately, as shown in Figure 3 and Figure 3 for the BMI dataset and in Figure 4 and Figure 4 for the Apple stock prices dataset, respectively.

4.2. Discussion

In both applications, the Skew-t distribution fits the data adequately and has the minimum AIC value compared to other distributions. This assertion is also well confirmed by the density and empirical CDF estimates as shown in Figure 3 and Figure 4. Through simulations, 2000 bootstrap samples are drawn from the Skew-t distribution with the specified parameter estimates in Table 4 and Table 5 and the approximate p-values of the proposed test for BMI and Apple stock daily return datasets are obtained as 0.5840 and 0.5473, respectively. We thus fail to reject the null hypothesis and conclude that both data can be modeled with the Skew-t distribution.

This case study highlights the test’s ability to validate Skew-t modeling in real data scenarios, particularly where skewness and heavy tails are present. It further illustrates the utility of energy distance-based approaches in detecting misfit across a range of candidate models.

5. Conclusions

We have developed a new goodness-of-fit test for Azzalini’s Skew-t distribution based on the energy distance framework. This method offers a flexible and powerful alternative to traditional EDF-based and correlation-based GoF tests. Our simulations demonstrate strong control of Type I error and superior power in a wide range of alternatives, particularly those involving skewness, heavy tails, or multi-modal distributions.

The proposed test performs reliably even in small to moderate sample sizes and is computationally feasible through bootstrap calibration and efficient implementation strategies. Real data analysis confirms its practical value in model validation contexts where accurate distributional assumptions are critical.

Future research directions include extensions to multivariate Skew-t distributions and further exploration of analytic approximations to the energy statistic under complex parametric families. In general, the energy-based goodness-of-fit test provides a robust and versatile tool for statistical modeling in both the theoretical and applied domains.