Specification procedures for multivariate stable-Paretian laws for independent and for conditionally heteroskedastic data

We consider goodness-of-fit methods for multivariate symmetric and asymmetric stable Paretian random vectors in arbitrary dimension. The methods are based on the empirical characteristic function and are implemented both in the i.i.d. context as well as for innovations in GARCH models. Asymptotic properties of the proposed procedures are discussed, while the finite-sample properties are illustrated by means of an extensive Monte Carlo study. The procedures are also applied to real data from the financial markets.


Introduction
Stable-Paretian (SP) distributions are extremely important from the theoretical point of view as they are closed under convolution, and also they are the only possible limit laws for normalized sums of i.i.d.random variables.In fact, this last feature makes SP laws particularly appealing for financial applications since stock returns and other financial assets often come in the form of sums of a large numbers of independent terms.Moreover the empirical density of such assets is leptokurtic and in many cases skewed, thus making the family of SP laws particularly suited for related applications; see for instance the celebrated stability hypothesis that goes back at least to Mandelbrot (1963) and Fama (1965).
The above findings prompted further research on the stochastic properties and inference of SP laws and on their application potential.The reader is referred to Samorodnitsky and Taqqu (1994), Adler et al. (1998), Uchaikin and Zolotarev (1998), Rachev and Mittnik (2000) and Nolan (2020) for an overview on the stochastic theory, statistical inference and applications of SP laws.
In the aforementioned works, statistical inference and applications are mostly restricted to univariate SP laws.The respective topics for SP random vectors are much less explored.In this connection, certain distributional aspects of SP random vectors are discussed in Press (1972b), while Press (1972a) defines moment-type estimators of the parameters of a multivariate SP distribution utilizing the characteristic function (CF).Maximum-likelihood based estimation is discussed by Ogata (2013), and Nolan (2013), Bayesian methods in Tsionas (2016), Lombardi and Veredas (2009) use indirect estimation methods, whereas Koutrouvelis (1980), Nolan (2013) and Sathe and Upadhye (2020) consider CF-regression methods.As far as testing is concerned, the only available formal method seems to be the CF-based goodness-of-fit test of Meintanis, Ngatchou-Wandji, and Taufer (2015).
In this article we propose CF-based goodness-of-fit procedures for SP random vectors in the elliptically symmetric case.Moreover we also consider the asymmetric case for which to the best of our knowledge goodness-of-fit have not been considered before.The remainder of the article unfolds as follows.In Section 2 we introduce a general goodness-of-fit test for elliptically symmetric SP laws.In Section 3 we particularize this test in terms of computation.Section 4 addresses the problem of estimation of the SP parameters and the study of asymptotics of the proposed test, while in Section 5 we extend the test to a multivariate GARCH model with SP innovations.The results of an extensive simulation study illustrating the finite-sample properties of the method are presented in Section 6.In Section 7 the case of asymmetric SP law is considered, with known as well as unknown characteristic exponent, by means of a different testing procedure that avoids integration over a complicated CF and computation of the corresponding density.The impact of high dimension on the test statistic is also discussed in this section.Applications are given in Section 8, and we conclude in Section 9 with a discussion.Some technical material is deferred to an online supplement along with additional Monte Carlo results.
2 Goodness-of-fit tests Let X be a random vector in general dimension p ≥ 1, with CF ϕ(t) = E(e it ⊤ X ), t ∈ R p , i = √ −1.Here we consider goodness-of-fit tests for the elliptically symmetric SP law.In this connection we note that SP random vectors are parameterized by the triplet (α, δ, Q), where α ∈ (0, 2] denotes the characteristic exponent, and δ ∈ R p and Q ∈ M p , are location vector and dispersion matrix, respectively, with M p being the set of all symmetric positive definite (p × p) matrices.On the basis of i.i.d.copies (X j , j = 1, ..., n) of X we wish to test the null hypothesis where we write X ∼ S α (δ, Q) to denote that X follows a SP law with the indicated parameters.For subsequent use we mention that if X ∼ S α (δ, Q), then it admits the stochastic representation where A is a totally skewed to the right SP random variable with characteristic exponent α/2 (α < 2) and N is a zero-mean Gaussian vector with covariance matrix Q, independent of A; see Samorodnitsky and Taqqu (1994, §2.5).
Our test will make use of the fact that if X ∼ S α (δ, Q), then the CF of X is given by ϕ α (t; δ, Q) = e it ⊤ δ−(t ⊤ Qt) α/2 .The cases α = 2 and α = 1, respectively, correspond to the most prominent members of the SP family, i.e. the Gaussian and Cauchy distributions, while φ α (t) = exp(− t α ), will denote the CF of a spherical SP law, i.e. an SP law with location δ = 0 and dispersion matrix Q set equal to the identity matrix.
As it is already implicit in (1) the parameters (δ, Q) are considered unknown, and hence they will be replaced by estimators ( δ n , Q n ) and the test procedure will be applied on the standardized data Specifically for a non-negative weight function w(•) we propose the test criterion where φ n (t) = n −1 n j=1 exp(it ⊤ Y j ) is the empirical CF computed from (Y j , j = 1, ..., n).Here the characteristic exponent α is considered fixed at value α = α 0 , while the case of unknown α will be considered in Section 7.

Computational aspects
It is already transparent that the test statistic T n,w depends on the weight function w(•), the choice of which we consider in this section.Specifically, from (4) we have by simple algebra where (6)

Using the CF of the Kotz-type distribution
We now discuss the computation of the test statistic figuring in (5).In doing so we will make use of an appropriate weight function w(•) that facilitates explicit representations of the integrals in ( 6).Specifically we choose w(t) = ( t 2 ) ν e −( t 2 ) α 0 /2 as weight function, which for (r, s) = (1, α 0 /2) is proportional to the density of the spherical Kotz-type distribution c( x 2 ) ν e −r( x 2 ) s , with c a normalizing constant; see Nadarajah (2003).With this weight function the integrals in (6) may be derived as special cases of the integral with the cases of interest being 0 < s(= α 0 /2) ≤ 1.In turn this integral can be computed by making use of the CF of the Kotz-type distribution in a form of an absolutely convergent series of the type (see Streit, 1991;Iyengar & Tong, 1989;Kotz & Ostrovskii, 1994) ), which for selected values of s reduces to a finite sum.For more details we refer to Section S1 of the online supplement.
Given I ν,r (•; •), the test statistic figuring in (5) may be written as

Using the inversion theorem
In this section we will use the inversion theorem for CFs in order to compute the integrals defined by (6).Specifically for an absolutely integrable CF ϕ(t), the inversion theorem renders the density f (•) corresponding to ϕ(•) as In this connection, we start from the expression of the test statistic in (5), and adopt the weight function w(t) = e −r t α 0 .This choice amounts to taking ν = 0 in the Kotz-type density, which is the same as if we incorporate the CF φ α0 (•) of the SP law under test in the weight function.
With this weight function, the statistic figuring in (5), say T n,r , becomes where by making use of the inversion theorem in (9), with f α (•) being the density of the spherical SP law with CF φ α (•).
4 On estimation of parameters and limit properties of the test

Estimation of parameters
The parameters δ and Q in (1) are assumed unknown and need to be estimated.
Seeing that reliable procedures for calculation of stable densities are available, we use maximum likelihood estimation.To avoid searching over the space of positive definite matrices, we use the estimators δ n and where L p denotes the space of lower triangular p × p matrices, and, as before, f α (•) denotes the density of a p-dimensional stable distribution with CF φ α0 (•).Initial values for the optimization procedure are obtained using projection estimators of δ and Q as outlined in Nolan (2013, Section 2.3).

Limit null distribution and consistency
We present here the main elements involved in the limit behavior of the test statistic T n,w .In this connection and despite the fact that, as already emphasized, it is computationally convenient to use a weight function that is proportional to the density of a spherical Kotz-type distribution, our limit results apply under a general weight function satisfying certain assumptions and, under given regularity conditions, with arbitrary estimators of the distributional parameters.Specifically, we assume that the weight function satisfies w(t) > 0 (apart from a set of measure zero), w(t) = w(−t) and R p w(t)dt < ∞.
We also suppose that the estimators figuring in the standardization defined by (3) admit the Bahadur-type asymptotic representations Then we may write from (4) It also follows that Ξ n − Ξ n,0 2 P −→ 0, where the approximating process Ξ n,0 (•) admits the i.i.d.representation upon which the central limit theorem applies, and which together with a subsequent application of the continuous mapping theorem entails where Ξ 0 (•) is a zero-mean Gaussian process with covariance kernel, say, K(s, t).In turn the law of T w is that of ∞ j=1 λ j N 2 j , where (N j , j ≥ 1) are i.i.d.standard normal random variables.The covariance kernel K(•, •) of the limit process Ξ 0 (•) enters the distribution of T w via the eigenvalues In this connection the maximum likelihood estimators defined by (12) satisfy certain equivariance/invariance properties (refer to Section S2 of the online supplement), and as a consequence the resulting test statistic is affine invariant.In this case we may set (δ 0 , Q 0 ) equal to the zero vector and identity matrix, respectively, thus rendering the limit null distribution free of true parameter values; see Ebner and Henze (2020) and Meintanis et al. (2015).
Moreover the standing assumptions imply the strong consistency of the new test under fixed alternatives.
Proposition 4.1 Suppose that under the given law of X the estimators of the parameters δ and Q satisfy, ( δn, Qn) a.s. as n → ∞.
Proof Recall from (4) that Now the strong uniform consistency of the empirical CF in bounded intervals (see Csörgő, 1981) entails , an application of Lebesgue's theorem of dominated convergence on ( 14) yields ( 13).
Since w > 0, T w is positive unless T w > 0 unless the CF of X coincides with the CF of a SP law with α = α 0 and (δ, Q) = (δ X , Q X ), and thus by the uniqueness of CFs, the test which rejects the null hypothesis H 0 in (1) for large values of T n,w is consistent against each fixed alternative distribution.The above limit results have been developed in a series of papers, both in the current setting as well as in related settings, and with varying conditions on the weight function and the family of distributions under test; see for instance Henze and Wagner (1997), Gupta et al. (2004), Meintanis et al. (2015), Hadjicosta and Richards (2020a), and Ebner and Henze (2020).In this regard, the solution of the above integral equation, and thus the approximation of the the limit null distribution of T n,w , is extremely complicated, and in fact constitutes a research problem in itself.This line of research has been followed by a few authors.We refer to Matsui and Takemura (2008), Hadjicosta and Richards (2020b), Meintanis et al. (2023) and Ebner and Henze (2023).In these works, the infinite sum distribution of T w is approximated by a corresponding finite sum employing numerically computed eigenvalues and then large-sample critical points for T n,w are found by Monte Carlo.It should be noted that such approximation is specific in several problem-parameters: the distribution under test, the type of estimators of the distributional parameters, and the weight function employed, and thus have to be performed on a case-to-case basis.A different more heuristic approach is moment-matching between the first few moments of T w (computed numerically) and a known distribution, like the gamma distribution or one from a Pearson family of distributions; see Henze (1990Henze ( , 1997) ) and Pfister et al. (2018), while yet another, Satterthwaite-type, approximation is studied by Lindsay et al. (2008).
The validity and usefulness of the above approximation methods notwithstanding, we hereby favor Monte Carlo simulation and bootstrap resampling for the computation of critical points and for test implementation, not only because these are large-sample approximations performed mostly in the univariate setting and thus inappropriate for small sample size n and/or dimension p > 1, but more importantly because, in the case of GARCH-type observations considered in the next section, the finite-sample counterpart of this distribution may even involve true parameter values; see for instance Henze et al. (2019).

The case of the stable-GARCH model
Assume that observations (X j , j = 1, ..., n) arise from a multivariate GARCH model defined by where (ε j , j = 1, ..., n) is a sequence of i.i.d.p-dimensional random variables with mean zero and identity dispersion matrix, and Q j := Q(X j |I j−1 ), with I j denoting the information available at time j, is a symmetric positive definite matrix of dimension (p × p).We wish to test the null hypothesis stated in (1) for the innovations ε j figuring in model (15).Note that in view of (2), Q j may be interpreted as the conditional covariance matrix of the corresponding latent Gaussian vector N .
We also employ initial values in order to start the estimation process.As an estimator of ϑ we use the equation-by-equation (EbE) estimators proposed by Francq and Zakoïan (2016).The reader is referred to Section S3 of the online supplement for more details on the EbE estimator.

Numerical study
We  S8 and S9 of the online supplement.
For comparison, we include results obtained when using the test of Meintanis, Ngatchou-Wandji, and Taufer (2015) using the Gaussian weight function exp(− t 2 ).The test, denoted by M a in the tables, depends on a tuning parameter denoted by a for which we consider the choices a = 4, 6, 10 and a = 15.As pointed out by the authors, the number of operations required to compute their test statistic is in the order of n 2a (for integer-valued a) and becomes time-consuming for larger values of n and a.We therefore use the Monte Carlo approach suggested by the authors to approximate the value of the test statistic (using 1 000 replications for each approximation); see p. 180 of Meintanis et al. (2015).When testing for multivariate normality (i.e.H 0 with α = 2), we consider the test of Henze, Jiménez-Gamero, and Meintanis (2019), denoted by HJM in the tables, with weight function exp(−1.5 t 2 ), which yielded good results in the original paper.
The rejection percentages of the tests are shown in Tables 1-4.All simulations results are based on 1,000 independent Monte Carlo iterations and a significance level of 10% is used throughout.
We first consider the case where we test H 0 with α = 2, i.e. when testing for departures from multivariate normality.Table 1 shows that, when heaviertailed symmetric Paretian alternatives are considered, the newly proposed tests based on T r are more powerful than the existing tests M a of Meintanis et al. (2015), and have power slightly lower but comparable to the test HJM of Henze et al. (2019).In light of the above-mentioned computational complexity of the existing tests, this gain in power makes the new tests attractive for implementation in practice.Moreover, the favorable power is visible for all the considered choices of the tuning parameter r, the choice of which, as opposed to a in M a , has no significant impact on computational complexity.Finally, as is expected, the power of the new tests increase as the sample size is increased.Similar conclusions can also be made in the case of elliptically symmetric t alternatives (see Table S2 of the online supplement).
Considering skew normal alternatives, even more favorable behavior can be observed in the results shown in Table 2. Notice that the tests M a and HJM seem to have very low power against skew normal alternatives, which is not the case with the newly proposed tests.
We now shift our attention to the case of testing H 0 with α = 1.8.Table 3 shows that, compared to the existing tests, the new tests are quite powerful against heavier-tailed alternatives, i.e. alternatives with stability index less than 1.8.Despite the evident dependence of the performance of the new tests on the tuning parameter r, we note that the power is very competitive to the existing tests in most cases, and significantly outperforms the existing tests for most choices of r.For lighter-tailed alternative distributions, the new tests exhibit some under-rejection for certain choices of r.Nevertheless, the problem seems to disappear as the sample size is increased.
Another advantageous feature of tests based on T r is that the high power is also visible in the higher-dimensional setting where p = 6.In fact, in this case the power of the tests seem to increase as the dimension is increased from p = 4 to p = 6.In contrast, the tests M a exhibit a clear decrease in power as the dimension p is increased.
We finally consider testing H 0 with α = 1, i.e. when testing for departures from multivariate Cauchy.Overall, in agreement with previous observations, Table 4 shows that the new tests seem to be competitive in terms of power, with the existing tests having some advantage when the true data-generating distribution has a stability index greater than 1.However, this advantage in power disappears when considering the higher-dimensional case p = 6.Similar behavior can be seen in the case of t and skew Cauchy alternatives (see Tables S6 and S7 in the supplement).

Monte Carlo results for GARCH data
We consider a CCC-GARCH(1, 1) model as defined in ( 16) with κ x = κ q = 1.As parameters, we take A 1 = 0.2I p , B 1 = 0.3I p , and correlation matrix R with all off-diagonal entries set to 0.5.To test whether the innovations ε j were generated by an elliptically symmetric SP law, we use the statistic in (10) applied to the residuals ε j defined in ( 17).Denote the value of the test statistic by T r := T r ( ε 1 , . . ., ε n ).
We use the following bootstrap scheme to determine critical values: 1. Independently generate innovations ε * 1 , . . ., ε * n from the SP law S α0 (0, I p ). 2. Construct a bootstrap sample X * 1 , . . ., X * n using the recursive relation 16) to obtain estimates Q * j , j = 1, . . ., n, and recover the bootstrap residuals ε * j = ( Q * j ) −1/2 X * j , j = 1, . . ., n 4. A bootstrap version of the statistic is given by T Steps 1-4 are repeated many times, say B, to obtain realisations T * r (1), . . ., T * r (B) of the bootstrap statistic.The null hypothesis is rejected at significance level ξ whenever T r exceeds the (1 − ξ)-level empirical quantile of the bootstrap realisations {T * r (b)} B b=1 .In the Monte Carlo simulations that follow, instead of drawing B bootstrap, we employ the warp speed method of Giacomini et al. (2013) which involves drawing only one bootstrap sample for each Monte Carlo iteration.
Table 5 (power against symmetric SP laws) exhibits similar favorable power properties of the newly proposed test as was observed in the i.i.d.case.Notice that the tests all have empirical size close to the nominal level of 10% when using critical values obtained by means of the bootstrap scheme given above.Also see Tables S13 and S14 of the online supplement.
Table S12 (power against skew normal alternatives) shows that the new tests are especially competitive in terms of power when the true innovation distribution is not elliptically symmetric.Finally we mention that although the tests of Meintanis et al. (2015) are competitive when the innovations have an elliptically symmetric Student t-distribution (see Table S11), the slight advantage in power disappears rapidly as the dimension p increases.

Testing asymmetric SP laws
We now shift our focus to the more general case of testing whether multivariate observations from X originate from a SP law, which need not necessarily be elliptically symmetric.In this connection note that the general multivariate SP law depends on a location vector δ ∈ R p and a spectral measure Γ(•) on the unit sphere S p .Accordingly we wish to test the null hypothesis H 0 : X ∼ S α (δ, Γ), for some δ ∈ R p and some spectral measure Γ on S p , where we write X ∼ S α (δ, Γ), when X follows a skew SP law with the indicated parameters.We will be considering the case of H 0 for fixed α = α 0 as well as the case of testing the null hypothesis H 0 with an unspecified α.
In this connection note that if X ∼ S α (δ, Γ), then the CF of X is given by with See, e.g., Nolan et al. (2001).
In testing the null hypothesis H 0 , we consider an entirely different idea for the test statistic first put forward by F. Chen et al. (2022).Specifically, we consider a test statistic formulated as a two-sample test between the original data X n = (X j , j = 1, ..., n) and artificial data X 0n = (X 0j , j = 1, ..., n) generated under the null hypothesis H 0 .More precisely, we propose the test Specification procedures for multivariate stable-Paretian laws criterion where ϕ n (t) = n −1 n j=1 e it ⊤ Xj is the empirical CF of the data at hand, while ϕ 0n (t) = n −1 n j=1 e it ⊤ X0j is the empirical CF computed from the artificial data set X 0n generated under the null hypothesis H 0 with Γ and δ estimated from the original observations X n .
By straightforward computations, we obtain where I w (x) = R p cos(t ⊤ x)w(t)dt, as also defined in (6).Clearly then the numerical approaches of Section 3 are no longer required and the simplicity of this test lies in the fact that in (20) only the computation of I w (•) is needed.Specifically the need for tailor-made weight functions such as those employed in Section 3 is circumvented.Furthermore we no longer need to compute the density of the underlying SP law as in (11).In this connection, suppose that the weight function w(x) figuring in I w (x) above is chosen as the density of any spherical distribution in R p .Then the integral I w (x) gives the CF corresponding to this spherical distribution at the point x.Furthermore it is well known that this CF may be written as Ψ( x 2 ), where Ψ(•) is called the "kernel" of the specific family of spherical distributions.Thus the test statistic figuring in (20) may be written as where the kernel Ψ(•) can be chosen by the practitioner so that the resulting expression in ( 21) is as simple as possible.In this connection, as already clear from the preceding paragraphs, a simple kernel is the kernel of the spherical SP family of distributions with Ψ(ξ) = e −rξ α/2 , r > 0, α ∈ (0, 2].Implementation of the test however relies on estimation of the spectral measure Γ(•) appearing in (18).Motivated by a result of Byczkowski et al. (1993, Theorem 1), we assume that Γ(•) can be approximated by the discrete spectral measure Γ(•) = K k=1 γ k I s k (•), with weights γ k corresponding to mass points s k ∈ S p , k = 1, . . ., K, and I s k (•) being the indicator index.So in order to apply the test, we use the stochastic representation of X 0 ∼ S α0 (δ, Γ) as where (A k , k = 1, ..., K) are i.i.d.(univariate) SP variates following a totally skewed to the right SP law with α = α 0 .In turn this representation is used in order to generate observations (X 0j , j = 1, ..., n) under the null hypothesis H 0 with δ and (γ k , k = 1, ..., K) replaced by appropriate estimates δ and ( γ k , k = 1, ..., K), respectively.The estimates δ and ( γ k , k = 1, ..., K) are obtained as shown in Nolan et al. (2001) and outlined in Section S5 of the online supplement.

Monte Carlo results
We now turn to a simulation study to demonstrate the performance of the test based on (20), say T Ψ for simplicity, in the bivariate case.We specifically consider the following alternative distributions: (A1) the asymmetric SP law S α (δ, Γ) with CF defined in ( 18), where we take δ = 0 and Γ( (A2) spherically symmetric SP distributions, denoted by S α .
The statistic in T Ψ (X n ; X 0n ) in ( 21) is subject to randomness introduced by the artificial data X 0n .To address this randomness in practical implementation, we follow F. Chen et al. (2022) and base our test on the statistic where, for each (r = 1, . . ., m), the set X r 0n is a random sample of observations generated from the SP law S α0 ( δ, Γ), i.e. a random sample satisfying H 0 .
Critical values of the test can be obtained using the following parametric bootstrap scheme: 1. Generate a bootstrap sample X * n from the SP law S α0 ( δ, Γ). 2. Calculate bootstrap estimates δ * and Γ * from X * n and generate m random samples (X * r 0,n , r = 1, . . ., m), from the SP law S α0 ( δ * , Γ * ).6 shows the empirical rejection percentages of tests based on T Ψ with Ψ(ξ) = e −rξ α 0 /2 , where α 0 denotes the hypothesized stability index in H 0 .We write T r , r > 0, for this test.In the case of unknown α in Table 7 we use the same weight function with α 0 = 2.

A bootstrap version of the statistic is T
The left-hand side of Table 6 shows the rejection percentages when observations are sampled from an asymmetric SP distribution, whereas the right-hand side show the results when observations are sampled from a symmetric SP distribution.In all cases, the proposed procedure seems to respect the nominal size of the test, although it seems to be somewhat conservative especially for smaller sample sizes.Nevertheless, the tests have good power against alternatives, which increases with the extent of violation of the null hypothesis.Corresponding results when testing H 0 with α = 1.8 are given in Table S16.

The case of unknown α
In practice, the true value of stability index α will typically be unknown and need to be estimated from data.Suppose the hypothesis of interest is To test this hypothesis, we again use the test based on the statistic in ( 21), but now generate the artificial data X 0n from the SP law S α ( δ, Γ), where α, δ and Γ are projection estimates obtained as outlined in Section S5 of the supplement.
Note that data can be generated from S α ( δ, Γ) using the stocahstic representation in ( 22) with α 0 replaced by α.The bootstrap procedure for obtaining critical values follows similarly to the procedure described in Section 7.1.The empirical rejection percentages (for the bivariate case, p = 2) are shown in Table 7 for several distributions.Note that when observations are  As alternatives, we consider the skew normal (SN ν ), symmetric Laplace and generalized Gaussian (GG ν ) distributions.An observation X from the GG ν distribution (refer to Cadirci et al., 2022) is generated according to X = U V 1/ν , where U is uniform on S p−1 and V ∼ Gamma(p/ν, 2).Table 7 shows that the proposed test procedure has power against non-SP alternatives, and noting the increase in power associated with increasing sample size, the results suggest that the procedure is consistent against non-SP alternatives.

The high-dimensional case
It should be clear that so far, and despite the fact that the new test applies to any dimension p, the underlying setting is not that of high dimension (p > n).In this connection we point out that the extension of goodness-of-fit methods specifically tailored for the classical "small p-large n" regime to high or infinite dimension is not straightforward, and therefore it is beyond the scope of the present article.This is not restricted to our setting alone but applies more generally, and the collections of Goia and Vieu (2016) and Kokoszka et al. (2017) reflect the need for statistical methods specifically tailored for nonclassical settings.If we restrict it to our context, in the latter settings such methods have so far been mostly confined to testing for normality and the interested reader is referred to Bugni et al. (2009), Nieto-Reyes et al. ( 2014), Bárcenas et al. (2017), Kellner and Celisse (2019), Yamada and Himeno (2019), Specification procedures for multivariate stable-Paretian laws Jiang et al. (2019), Górecki et al. (2020), Henze and Jiménez-Gamero (2021), H. Chen and Xia (2023) and W. Chen and Genton (2023).If however the setting is that of regression, with or without conditional heteroskedasticity, the number of parameters rapidly increases with the dimension p.This is one extra reason that corresponding specification methods in high/infinite dimension need to be treated separately, and in this connection the methods of Cuesta-Albertos et al. ( 2019) and Rice et al. (2020) appear to be of the few available for testing the (auto)regression function.
For an illustration of the special circumstances that are brought forward as the dimension grows, consider the test statistic in (10) for α 0 = 2 and without standardization, i.e. replace Y j by X j ∼ S 2 (0, I p ), j = 1, ..., n.Then from (11), and by using the density of the p-variate normal distribution with mean zero and covariance matrix 2I p , we obtain Λ r (x; 2) = π r p 2 e − x 2 /(4r) , and consequently , and thus our test statistic contains sums of terms, each of which is of the order e −2p , as p → ∞. (For simplicity we suppress the terms 4r and 4(r + 1) which occur as denominators in these sums, as they are anyway irrelevant to our argument).
In order to get a feeling of this result, write X j 2 = 2 p k=1 (X jk / √ 2) 2 =: 2S p , where, due to Gaussianity and independence, S p is distributed as chisquared with p degrees of freedom.Thus the expectation of the quantity e − Xj 2 figuring in T n,r coincides with the value of the CF of this chi-squared distribution computed at the point t = 2i.To proceed further, notice that the CF of S p /p at a point t is given by the CF corresponding to S p computed at t/p, and recall that the CF of the chi-squared distribution with p degrees of freedom is given by ϕ Sp (t) = (1 − 2it) −p/2 .Hence, in obvious notation, and reverting back to the CF of S p we get ϕ Sp (t) ≈ e ipt , and thus ϕ Sp (2i) ≈ e −2p , as p → ∞.A similar reasoning applies to X j − X k 2 implying that the expectation of e − Xj −X k 2 (also occurring in T n,r ) being approximated by e −4p , as p → ∞.Hence, by using these approximations, we obtain as p → ∞.Consequently our test statistic degenerates in high dimension, a fact that calls for proper high-dimensional modifications of the test criterion, which is definitely a worthwhile subject for future research.
Nevertheless, we have obtained some initial Monte Carlo results that show a reasonable performance for the test criterion in cases where the dimension is much higher than the maximum p = 6 considered so far.These results may be found in Tables S18 and S19 of the Supplement.

Application to financial data
We consider daily log returns from 4 January 2010 to 30 June 2017 of stocks of two mining companies listed on the London Stock Exchange: Anglo American (AAL) and Rio Tinto (RIO).The complete data set (available from Yahoo! Finance) consists of 1,891 log returns.We model the log returns using the CCC-GARCH(1, 1) model (with intercept) given by where Q j is as in ( 16) with κ x = κ q = 1 and B 1 assumed diagonal.We are interested in determining whether the innovations ε j have a bivariate stable distribution.Since the innovations are unobserved, we apply our test to the residuals ε j = Q −1/2 j (Y j − ω), where the estimates Q j and ω are obtained using EbE estimation as discussed earlier.To obtain critical values of the tests, we apply the bootstrap algorithm of Section 6.2 with B = 1, 000.Table 8 show that, when a stability index of α ∈ {1.75, 1.8} is assumed, the tests based on T r do not reject the null hypothesis that the CCC-GARCH-(1, 1) innovations have a S α distribution (at a 10% level of significance).On the other hand, the null hypothesis of stable innovations is rejected when a stability index of α = {1.7,1.85, 1.9, 2} is assumed.
The correct choice of the innovation distribution has important implications for value-at-risk (VaR) forecasts.For the considered stable distributions, we fit the model in ( 25) to the first 1,000 observations and calculated onestep-ahead 5% and 1% portfolio VaR forecasts for long and short positions for the remaining time period (i.e.891 forecasts for each position).The portfolio is assumed to consist of 50% AAL shares and 50% RIO shares.
Table 8 shows the empirical coverage rates of the forecasted VaR bounds, that is, the proportion of times that the value of the portfolio exceeded the bounds.For the cases where our test supports the null hypothesis, i.e. when α ∈ {1.75, 1.8}, the empirical coverage rates of the value-at-risk bounds are quite close to the nominal rates.In addition, the p-values of Christofferson's (1998) LR cc test (given in brackets in Table 8), indicate that, if the stability index of the innovation distribution is chosen either too high or too low, the true conditional coverage rates of the forecasted VaR bounds are significantly different from the nominal rates.

Conclusion
We have studied goodness-of-fit tests with data involving multivariate SP laws.Our tests include the case of i.i.d.observations as well as the one with observations from GARCH models, and cover both elliptical and skewed distributions.Moreover they refer to hypotheses whereby some parameters are assumed known, as well as to the full composite hypothesis with all parameters estimated from the data at hand.The new procedures are shown to perform well in finite samples and to be competitive against other methods, whenever such methods are available.An application illustrates the usefulness of the new procedure for modeling stock returns and explores the subsequent forecasting implications.
required for the test statistic T n,w can be computed by making use of the CF of the Kotz-type distribution.Specifically for x = 0, we have .
On the other hand, if x = 0, the value of I ν,r (x; s) may be computed from absolutely convergent series which for 1/2 < s < 1 is given by Streit (1991) as , while for 0 < s < 1/2 we have see Kotz and Ostrovskii (1994).
In the special case s = 1 the computation simplifies to , see Iyengar and Tong (1989), and for s = 1/2 the bivariate case was treated by Nadarajah and Kotz (2001) leading to which for ν = 0 simplifies to

S2 Affine invariance
A desirable feature of potential estimators δ n and Q n of the location vector δ and the dispersion matrix Q are the following equivariance/invariance properties for each b ∈ R p and each non-singular (p × p) matrix A. As a consequence the test statistic T n,w := T n,w (X 1 , ..., X n ) satisfies i.e. it is affine invariant.We note that this property is in line with the fact that if X ∼ S α (δ, Q), then AX + b ∼ S α (Aδ + b, AQA ⊤ ), meaning that the SP family of distributions is itself invariant with respect to affine transformations X → AX + b.

S3 EbE estimator for GARCH parameters
We outline the procedure for computing the EbE estimator of ϑ of Francq and Zakoïan (2016).Note in this connection that, under the CCC-GARCH model, n is given by , where fα denotes the density of a symmetric univariate SP law with stability index equal to α, with location zero and dispersion equal to one.Letting D j = diag( q 1,j ( ϑ n )), we then obtain the correlated residuals ε ′ j = D −1 j X j from which R is calculated using maximum likelihood to obtain

S4 Calculation of the stable density
Implementation of our test procedure relies on the evaluation of the density of a p-variate spherical SP law.Efficient evaluation of the density is also needed for maximum likelihood estimation of the parameters.
In our numerical work, we utilize the fact that if X follows a spherical SP law with CF φ α (•), then the density of X can be expressed as where f R (•) is the density of X , the amplitude of X.This reduces the problem of calculating f α to calculating the univariate density f R .Various integral expressions for f R are given in Nolan (2013) and an implementation to numerically evaluate f R is the function damplitude in the R package stable (provided by Robust Analysis Inc., 2016).
To speed up calculations, we pre-calculate f R (u) for each u ∈ {k/(N − k), k = 0, . . ., N − 1}, for some large value of N (= 10, 000 in our simulations).Intermediary points are approximated using cubic spline interpolation and, for u ≥ N , we set f R (u) = 0.

S5 Estimation of the discrete spectral measure
Below we outline the projection estimation procedure of Nolan et al (2001) as implemented in our work.The estimation procedure assumes that the stability index α is known and that the data have been centered.As this is usually not the case, we first estimate the one-dimensional parameters (α j , β j , σ j , δ j ), j = 1, . . ., p, for each of the coordinates of the p-dimensional data set and center the data set using the location estimate δ = (δ 1 , . . ., δ p ).Furthermore, if α is unknown, we estimate it by α = p −1 p j=1 α j .Motivated by a Theorem 1 of Byczkowski et al (1993), we assume that Γ(•) can be approximated by the discrete spectral measure Γ(•) = K k=1 γ k I s k (•) as defined in the main paper, where we take For all estimates calculated in our simulations, we took the number of projections as K = 10.
For each k = 1, . . ., K, the projections X ⊤ 1 s k , . . ., X ⊤ n s k are i.i.d. with a centered, univariate α-stable distribution with skewness and dispersion parameters β k and σ k .These parameters are estimated using maximum likelihood, after which the estimates of the K projections are combined using the method described in Section 2.2 of Nolan et al (2001) to recover an estimate of the spectral measure Γ above.
Table S1 Percentage of rejection of H 0 with α = 2 against stable alternatives.Tests done at the 10% significance level.Due to the computational complexity of the HJM test, it is only included in selected cases to reduce run time.99.9 99.9 99.9 99.8 41.0 36.0 17.5 12.7 11.5 100.0 t4 99.9 99.9 100.0 100.0 63.5 50.Table S3 Percentage of rejection of H 0 with α = 2 against skew normal alternatives.Tests done at the 10% significance level.Due to the computational complexity of the HJM test, it is only included in selected cases to reduce run time.

Table 1
Percentage of rejection of H 0 with α = 2 against stable alternatives.Tests done at the 10% significance level.Additional cases are shown TableS1of the online supplement.

Table 2
Percentage of rejection of H 0 with α = 2 against skew normal alternatives.Tests done at the 10% significance level.Also see TableS3of the online supplement.

Table 3
Percentage of rejection of H 0 with α = 1.8 against stable alternatives.Tests done at the 10% significance level.Also see TableS4of the online supplement.

Table 4
Percentage of rejection of H 0 with α = 1 against stable alternatives.Tests done at the 10% significance level.Also see TableS5of the supplement.

Table 5
Percentage of rejection of H 0 with α = 2 against stable alternatives in the CCC-GARCH(1, 1) case.Also see TableS10of the online supplement.

Table 6
Percentage of rejection of H 0 with α = 1.5 against bivariate stable alternatives.Tests done at the 10% significance level.Also see TableS15of the online supplement.

Table 7
Percentage of rejection of H ′ 0 using tests based on the statistic in (21) with Ψ(ξ) = e −rξ .Tests done at the 10% significance level.Also see TableS17of the online supplement.
sampled from a S 1.8 distribution (symmetric SP law with stability index 1.8), the rejection percentages are close to the nominal level, indicating that the test has reasonable empirical size.

Table 8 p
-values of tests that the CCC-GARCH innovations have a Sα distribution for several choices of α, along with empirical coverage rates of forecasted VaR bounds calculated under H 0 .p-values of the LRcc test are given in brackets.All p-values less than 10% are underlined.

Table S2
Percentage of rejection of H 0 with α = 2 against t alternatives.Tests done at the 10% significance level.Due to the computational complexity of the HJM test, it is only included in selected cases to reduce run time.
TableS5Percentage of rejection of H 0 with α = 1 against stable alternatives.Tests done at the 10% significance level.

Table S6
Percentage of rejection of H 0 with α = 1 against Student t alternatives.Tests done at the 10% significance level.
TableS7Percentage of rejection of H 0 with α = 1 against skew Cauchy alternatives.Tests done at the 10% significance level.

Table S8
Percentage of rejection of H 0 with α = 2 against stable alternatives when using the test in (8) with the Kotz-type weight function.Tests done at the 10% significance level.TableS9Percentage of rejection of H 0 with α = 1 against stable alternatives when using the test in (8) with the Kotz-type weight function.Tests done at the 10% significance level.

Table S10
Percentage of rejection of H 0 with α = 2 against stable alternatives in the CCC-GARCH(1, 1) case.Tests done at the 10% significance level.TableS11Percentage of rejection of H 0 with α = 2 against t alternatives (GARCH model errors).Tests done at the 10% significance level.TableS12Percentage of rejection of H 0 with α = 2 against skew normal alternatives in the CCC-GARCH(1, 1) case.Tests done at the 10% significance level.Due to the computational complexity of the HJM test, it is only included in selected cases to reduce run time.TableS14Percentage of rejection of H 0 with α = 1.5 against stable alternatives (GARCH model errors).Tests done at the 10% significance level.TableS15Percentage of rejection of H 0 with α = 1.5 against stable alternatives.Tests done at the 10% significance level.TableS16Percentage of rejection of H 0 with α = 1.8 against stable alternatives.Tests done at the 10% significance level.TableS17Percentage of rejection of H ′ 0 using tests based on the statistic in (21) with Ψ(ξ) = e −rξ .Tests done at the 10% significance level.TableS19Percentage of rejection of H 0 with α = 1.5 against stable alternatives.For these results, Q and δ were assumed known and not estimated.Tests done at the 10% significance level.