
==== Front
J Appl Stat
J Appl Stat
Journal of Applied Statistics
0266-4763
1360-0532
Taylor & Francis

2315450
10.1080/02664763.2024.2315450
Version of Record
Articles
Research Article
Testing nonlinearity of heavy-tailed time series
JOURNAL OF APPLIED STATISTICS
J. G. De Gooijer
De Gooijer Jan G.
Amsterdam School of Economics, University of Amsterdam, Amsterdam, The Netherlands
CONTACT Jan G. De Gooijer j.degooijer@contact.uva.nl Amsterdam School of Economics, University of Amsterdam, PO Box 15867, Amsterdam 1001 NJ, The Netherlands
11 2 2024
2024
11 2 2024
51 13 26722689
6 5 2023
14 1 2024
Nova techset5 2 2024
Converted to JATS 1.2 by Nova Techset5 2 2024
© 2024 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.
2024
The Author(s)
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives License (http://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited, and is not altered, transformed, or built upon in any way. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.

A test statistic for nonlinearity of a given heavy-tailed time series process is constructed, based on the sub-sample stability of Gini-based sample autocorrelations. The finite-sample performance of the proposed test is evaluated in a Monte Carlo study and compared to a similar test based on the sub-sample stability of a heavy-tailed analogue of the conventional sample autocorrelation function. In terms of size and power properties, the quality of our test outperforms a nonlinearity test for heavy-tailed time series processes proposed by [S.I. Resnick and E. Van den Berg, A test for nonlinearity of time series with infinite variance, Extremes 3 (2000), pp. 145–172.]. A nonlinear Pareto-type autoregressive process and a nonlinear Pareto-type moving average process are used as alternative specifications when comparing the power of the proposed test statistic. The efficacy of the test is illustrated via the analysis of a heavy-tailed actuarial data set and two time series of Ethernet traffic.

Keywords

Gini-based autocorrelation
heavy tails
nonlinear Pareto-type models
sub-sample stability
nonlinearity tests
AMS CLASSIFICATIONS

62F03
62M10
62P05
==== Body
pmc1. Introduction

It is not uncommon in practice to encounter time series which exhibit non-standard features such as nonlinearity and heavy tails.1 For instance, many data sets in areas such as hydrology [19], internet traffic [29], and finance [21] appear to be compatible with the assumption of heavy-tailed marginal distributions; see, e.g. Castillo et al. [7] and the examples therein. Additionally, the book by Embrechts et al. [13] contains many examples involving heavy-tailed actuarial data sets; see also Beirlant et al. [4] and Section 7.1 of this paper.

Various types of heavy-tailed distributions have been proposed such as the Pareto, Student-t, log-normal, Weibull, and Cauchy distribution. Many of those distributions only possess a finite number of moments. This requires special consideration when establishing limit properties of estimation methods, nonlinearity test statistics, and diagnostic tests in the case of nonlinear time series with heavy tails. For instance, under the assumptions of finite-variance errors, asymptotic theory of the well-known Pearson sample autocorrelation function (ACF) may not hold in the presence of heavy-tailed nonlinear time series models unless moment conditions are satisfied. It appears, that for many of these models this statistic converges (in probability) to a non-degenerate random variable (r.v.); see, e.g. Davis and Resnick [12].

Motivated by the above consideration and in the absence of a proper nonlinearity tests of time series with infinite variance, Resnick and Van den Berg [32] proposed a test statistic based on the difference between two sub-sample ACFs (Section 4). They established the limiting null distribution of the test statistic. The null hypothesis of linearity is a moving average process driven by Pareto distributed innovations. The finite-but-large sample performance of the proposed test is assessed by three independent realizations of length 10,000 for a simple bilinear data generating process with Pareto distributed innovations. The test statistic has not been studied for small and moderate sample sizes using more elaborate nonlinear heavy-tailed processes.

In this paper, we develop and study a nonlinearity test statistic for heavy-tailed time series based on sub-sample stability of the Gini autocorrelation function (G-ACF). Clearly, this idea parallels the nonlinearity test statistic for heavy-tailed time series briefly discussed above. However, Carcea and Serfling [6] and Shelef and Schechtman [37] showed that G-ACFs are useful in measuring association between time series from heavy-tailed distributions; they can find a balance between robustness and efficiency [cf. 35, 39]. Moreover, G-ACFs require only first-order moment assumptions. Another difference from the existing definition of the ACF is that the G-ACF is not symmetric in its time series variables in general. Since there is an asymmetric inertia in many economic time series this particular feature of the G-ACF can be useful in finding a preliminary nonlinear model specification; see Carcea [5] and Shelef [36] for further discussions and applications. Another objective of the current paper is to compare the finite-sample performance of our test with the Resnick–Van den Berg test statistic within the context of analysing heavy-tailed Pareto-type nonlinear time series processes.

The rest of the paper is organized as follows. Section 2 briefly introduces Gini-based autocovariances and autocorrelations. Both population and sample quantities will be presented. In Section 3, we first discuss the ACF-based test for nonlinearity. Next, we present a similar test statistic based on sub-sample Gini-based ACFs. In Section 5, we introduce two Pareto-type nonlinear time series models. These models are used in Section 6 as alternative specifications when investigating the power of the two test statistics. Section 7 provides two empirical applications of the nonlinearity test statistics. Finally, we end the paper in Section 8 with discussing some possible extensions.

2. Gini autocovariances and autocorrelations

2.1. Population

Let {Yt,t∈Z} be a strictly stationary discrete time series process with finite first-order moment and a continuous distribution function F(Yt). The first-order moment assumption allows for developing concepts and methods in situations when the second-order moment assumption is not satisfied (E(Yt2)=∞) as in the case of heavy-tailed distributions and data. However, care must be taken in model selection using the traditional sample ACF as an estimator of the lag ℓ, ℓ∈N. In fact, in a series of papers, Resnick and his co-authors showed how dangerous it is to rely on sample ACFs as a model fitting tool in the heavy-tailed case; see, e.g. Cohen et al. [10], Resnick [30], and Resnick et al. [31].

The Pearson (P) ACF ρ(P)(ℓ)=Cov(Yt,Yt−ℓ)/Var(Yt) is defined under second-order moment assumptions. It measures the degree of linearity in the relationship between Yt and Yt−ℓ. As an alternative and due to its potential application in infinite variance time series, the G-ACF has been suggested; see, e.g. Carcea and Serfling [6] and Shelef and Schechtman [37]. For lag ℓ, two notions of this statistic are defined as follows (1) ρ(G1)(ℓ)=γ(G1)(ℓ)/γ(G)(0)andρ(G2)(ℓ)=γ(G2)(ℓ)/γ(G)(0),

where, paralleling the Gini covariance for a bivariate random vector with joint distribution function [35], the autocovariances for lag ℓ are defined by

and

where E[F(Yt)]=1/2 since F(Yt) is continuous and symmetric about zero, and with γ(G)(0)≡γ(G1)(0)=γ(G2)(0). Both functions ρ(G1)(ℓ) and ρ(G2)(ℓ) take values in the interval [−1,1].

The functions γ(G1)(ℓ) and γ(G2)(ℓ) are defined under merely first-order assumptions. The first quantity measures the degree of monotonicity in representing Yt as a function of F(Yt−ℓ), while γ(G2)(ℓ) measures the degree of monotonicity in representing F(Yt) as a function of Yt−ℓ. Observe that these functions measure asymmetry in a data generating process via forward and backward dependencies.

It is easy to see that Cov(Yt,F(Yt−ℓ))≠Cov(Yt−ℓ,F(Yt)) in general. If the distribution of the random vector (Yt,F(Yt+ℓ)) is exchangeable for all t∈Z and lag ℓ∈N up to a linear transformation then Cov(Yt,F(Yt−ℓ))=Cov(Yt−ℓ,F(Yt)). As an example, consider the following stationary autoregressive process of order one (AR(1)): Yt=ϕYt−1+ϵt with innovations {ϵt}. Then, γ(G1)(ℓ)=Cov(ϕYt−1+ϵt,F(Yt−ℓ))=ϕℓγ(G1)(0). But γ(G2)(ℓ)=Cov(Yt−ℓ,F(ϕYt−1+ϵt)) is not necessarily equal to γ(G1)(ℓ). In this case ρ(G1)(ℓ)=ϕℓ=ρ(P)(ℓ), which indicates that the estimator for ρ(G1)(ℓ) (to be defined in Section 2.2) can be used to measure the autocorrelation when the second moment does not exist. Since γ(G1)(ℓ)=γ(G2)(−ℓ) and γ(G2)(ℓ)=γ(G1)(−ℓ), we focus only on an estimator of ρ(G1)(ℓ) in the rest of the paper, and replace the superscript G1 by G.

2.2. Sample

Let {Yt}t=1n represent a set of observations generated by a strictly stationary process {Yt,t∈Z} under first-order assumptions. Following Shelef and Schechtman [37], and parallel to the classical centered sample ACF, an estimator of ρ(G)(ℓ) is given by (2) ρˆ(G)(ℓ)=γˆ(G)(ℓ)γˆ(G)(0)=∑t=1n−ℓ(Yt+ℓ−Y¯((ℓ+1):n))(R(Yt)−R¯(Y(1:(n−ℓ)))/(n−ℓ−1)∑t=1n(Yt−Y¯(1:n))(R(Yt)−R¯(Y(1:n))/(n−1),

where Y¯((ℓ+1):n)=∑t=ℓ+1nYt/(n−ℓ), R(Yt) denotes the rank of Yt (divided by the sample size), and R¯(Y(i:j))=∑t=ijR(Yt)/(j−i+1). Note that the ranks used in the numerator of (2) are based on the set of observations {(Y1+ℓ,Y1),(Y2+ℓ,Y2),…,(Yn,Yn−ℓ)}, whereas the ranks in its denominator are based on all observations.

As an alternative of γˆ(G)(0), Carcea and Serfling [6] proposed the following centered estimator of γ(G)(0): (3) γ~(G)(0)=2n2∑t=1n(2t−n)Y(t:n).

For the Gini autocovariance estimator at lag ℓ>0, these authors first introduced the sample distribution function Fˆ(⋅) as an unbiased estimator of F(Yt), i.e. Fˆ(y)=n−1∑t=1nI(Yt≤y),−∞<y<∞,

where I(⋅) is the indicator function, i.e. I(z)=1 if z>1 and I(z)=0 if z≤0. Next, parallel to the definition of the standard Gini covariance function for a pair of random variables (X,Y), Carcea and Serfling [6] define the following estimator of γ(G)(ℓ): (4) γ~(G)(ℓ)=2n−ℓ∑t=1n−ℓ(2Fˆ(Y(t:n−ℓ))−1)Y(t:n−ℓ),1,2,ℓ≥1.

Here Y(t:n−ℓ),1,2 is the tth ordered (i.e. sorted in size) first component value that is concomitant to the tth ordered second component value, relative to the (n−ℓ) lag ℓ bivariate pairs of observations in the set {(Y1+ℓ,Y1),(Y2+ℓ,Y2),…,(Yn,Yn−ℓ)}. The difference between γˆ(G)(ℓ) and γ~(G)(ℓ) (ℓ=0,1,2,…) goes to zero as n→∞; see [Appendix 1][37]. Asymptotic normality of γ~(G)(0) follows from Sandström [33, Thm. 1] using Gâteaux differentiation as a first order approximation to its estimation error; see Kattumannil [17] and Kattumannil et al. [18].

3. Testing for nonlinearity

3.1. Null hypothesis

Many time series models for the infinite variance case are based on symmetric α-stable r.v.'s, 0<α<2. In this section, we assume that a sequence {Zt,t∈Z} of i.i.d. r.v.'s (i.e. the innovation process) follows a Pareto type III distribution with location parameter zero and we write Zt∼i.i.d.P(σ,α). Its survival function is of the form F¯Z(z)=P(Z>z)=[1+(z/σ)α]−1 (z≥0), where the parameters σ>0 and α>0 correspond to scale and tail index respectively. Arnold [2] discusses various properties of the Type III P(σ,α) distribution. The δth moment of Zt can be expressed as E(Ztδ)=σδ(π/α)csc⁡(πδ/α), δ∈(−α,α). It is well known that E(Ztδ)=∞ for δ≥ α, i.e. Var(Zt)=∞ for 0<α≤2.

The null hypothesis of interest is the following discrete time moving average (MA) process (5) H0:Yt=∑j=0∞ψjZt−j,(ψ0=1),

where {ψj,j∈N} is a sequence of constants such that ∑j|ψj|ξ<∞ for some ξ<min(1,α).

3.2. Uncentered sub-sample ACFs

Clearly, when Var(Yt)=∞ there is no need for centering the data at Y¯≡Y¯(1:n)=∑t=1nYt/n. Then one can use the following heavy-tailed analogue of the sample Pearson-based ACF: (6) ρˆ(ℓ)=γˆ(ℓ)/γˆ(0),ℓ∈{0,1,…,m},

where γˆ(ℓ)=∑t=1n−ℓYtYt+ℓ and m≪n is a positive integer. Davis and Resnick [12] have shown that for heavy-tailed linear MA (∞) processes ρˆ(ℓ)→Pρ(ℓ),ℓ≥0,

where (7) ρ(ℓ)=∑j=0∞ψjψj+ℓ∑j=0∞ψj2.

Observe that ρ(ℓ) in (7) is not the usual theoretical Pearson-based ACF which does not exist in this case. Consequently, if the model is nonlinear, ρˆ(ℓ) does not converge in probability to a constant but instead converges weakly to a nondegenerate r.v. So deviations from the linearity assumption for a heavy-tailed model which is not MA (∞) can be reflected in large differences between uncentered sample ACFs, computed for two non-overlapping sub-samples {Yt1}t1=1N and {Yt2}t2=N+12N of {Yt}t=1n of equal length N, say. The resulting uncentered sample ACFs are given by (8) ρˆ(1)(ℓ)=∑t1=1N−ℓYt1Yt1+ℓ∑t1=1NYt12andρˆ(2)(ℓ)=∑t2=N+12N−ℓYt2Yt2+ℓ∑t2=N+12NYt22,

where N=n/2 when n is even, and N=(n−1)/2 if n is odd. Analogue to (8), we define the uncentered sub-sample Gini-based sample ACFs as: (9) ρˆ(1)(G)(ℓ)=∑t1=1N−ℓYt1+ℓR(Yt1)∑t1=1NYt1R(Yt1).,ρˆ(2)(G)(ℓ)=∑t2=N+12N−ℓYt2+ℓR(Yt2)∑t2=N+12NYt2R(Yt2)..

These statistics are natural estimators of E[YtF(Yt−ℓ)]/E[YtF(Yt)].

4. Test statistics

In this section, we first introduce in Theorem 4.1 a well-established test statistic for nonlinearity of heavy-tailed time series processes using ρˆ(u)(ℓ), u = 1, 2. Next, we discuss a similar test statistic based on ρˆ(u)(G)(ℓ).

4.1. A test based on sub-sample stability of ρˆ(ℓ)

First, we introduce four independent and stable r.v.'s, to be used in Theorem 4.1. In particular, let Sj,α, Sj,α′, j≥1, denote two i.i.d. α-stable r.v.'s with characteristic function E(eiuSα)=exp⁡{−Γ(1−α)cos⁡(πα/2)|u|α}, α≠1. Also, let S0,α/2, S0,α/2′ be two i.i.d. positive stable r.v.'s with index α/2 and characteristic function E(eiuS0,α/2)=exp⁡{−Γ(1−α/2)cos⁡(πα/4)|u|α/2(1−i⋅sign(u)tan⁡(πα/4))}; Zolotarev [43].2 Then, it can be shown [32] that the following theorem holds.

Theorem 4.1 Suppose that {Yt,t∈Z} follows the {MA}(∞) process defined in (5) with tail index α∈(0,2). Furthermore, if 0<α<1, let ρˆ(1)(ℓ) and ρˆ(2)(ℓ) denote the uncentered sub-sample ACFs defined in (8). Then, for a fixed value m>0, and as n→∞, (10) (nlog⁡n)1/α⋁ℓ=1m|ρˆ(1)(ℓ)−ρˆ(2)(ℓ)|⟶D⋁ℓ=1m|Y(ℓ)|,

where (11) Y(ℓ)=D∑j=1∞(ρ(ℓ+j)+ρ(ℓ−j)−2ρ(j)ρ(ℓ))(Sj,αS0,α/2−Sj,α′S0,α/2′),

and where the notation =D denotes equality in distribution.

The first term (n/log⁡n)1/α in (10) is a scaling constant given by Davis and Resnick [12] for (stable) Pareto-type distributions of |Zt|. Theorem 4.1 remains true if ρˆ1(ℓ) and ρˆ2(ℓ) are replaced by their mean-corrected versions. For practical purpose, Y(ℓ) in (11) can be replaced by the following estimator (12) Yˆ(ℓ)=∑j=1r(n)(ρˆ(ℓ+j)+ρˆ(ℓ−j)−2ρˆ(j)ρˆ(ℓ))(Sj,αˆS0,αˆ/2−Sj,αˆ′S0,αˆ/2′),ℓ∈{0,1,…,m},

where ρˆ(ℓ) is defined by (6), r(n)=o((log⁡n)n(2ξ/α)−1), and where αˆ is a consistent estimator of the tail index α. Now as an immediate result of Theorem 4.1 and (12), the Resnick–Van den Berg test statistic for linearity is given by (13) T(m)=(nlog⁡n)1/αˆ⋁ℓ=1m|ρˆ(1)(ℓ)−ρˆ(2)(ℓ)|,m>0.

For applications, it is quite common to use the well-known Hill maximum likelihood estimator as an estimator of α; Embrechts et al. [cf. 13]. Another recommendation is to replace r(n) in (12) by r(n)=m+⌊(n/2)0.499⌋ if αˆ<4/3, and r(n)=m+⌊(n/2)(2/αˆ)−1⌋ otherwise. Finally, for using the asymptotic distribution of T(m) it is recommended to choose the number of simulated limit vectors (Yˆ(1),…,Yˆ(m)) large, say 10,000.

4.2. A test based on sub-sample stability of ρˆ(G)(ℓ)

Similar to (13), the proposed Gini-based test statistic for linearity of heavy-tailed time series is given by (14) T(G)(m)=(nlog⁡n)1/αˆ⋁ℓ=1m|ρˆ(1)(G)(ℓ)−ρˆ(2)(G)(ℓ)|,m>0.

Finding the asymptotic distribution of T(G)(m), however, fails because limit theory for ρˆ(G)(ℓ) for MA(q) (q≥1) models with Pareto innovations is only available for the case q = 0; see Section 2.2.

As a logical solution to this problem, we employ the well-known moving block bootstrap procedure. For this procedure, the relevant simulation steps are as follows: (i) Generate a real-valued stationary time series process {Yt}t=1n of the form (5); (ii) Divide the series into b consecutive observations: Bi={Yi,Yi+1,…,Yi+b−1}, where i=1,…,n−b+1. In each bootstrap replication, resample (with replacement) k blocks B1∗,…,Bk∗ where k≈(n/b)∈N. Next, join these k blocks together side by side to construct the bootstrapped sample. Calculate the bootstrap sub-sample estimators ρˆ(u)(G,∗)(ℓ) (u=1,2;ℓ=1,…,m) and, for fixed values of αˆ and m, the associated test statistic T(G,∗)(m)=(n/log⁡n)1/αˆ⋁ℓ=1m|ρˆ(1)(G,∗)(ℓ)−ρˆ(2)(G,∗)(ℓ)|; (iii) Repeat step ii) a large number of times, M (typically M≥500), in order to obtain the set of bootstrapped test statistics {T(j)(G,∗)(m)}j=1M; iv) Using the sampling quantiles of {T(j)(G,∗)(m)}j=1M, we can determine critical values or we can approximate a p-value for T(G)(m) by calculating (#(T(j)(G,∗)(m)>T(G)(m))+1)/(M+1). In our study we set M = 500.

The above procedure can be easily modified to determine critical values of the test statistic (13) under the null hypothesis (5). In fact, to put the finite-sample comparison of T(m) and T(G)(m) on an equal footing, we use the block bootstrap procedure for both test statistics. The optimal block-length b is selected by cross-validation using the bootstrap procedure proposed by Hall et al. [16].3

Lahiri [22] introduced an empirical example where the moving block bootstrap surpasses all other versions of block bootstrap methods in terms of mean squared error (MSE). Radovanov and Marcikí [27] compared four different block bootstrap methods: non-overlapping block bootstrap, overlapping block bootstrap (moving block bootstrap), stationary block bootstrap and subsampling. Overall, the moving block bootstrap procedure shows the lowest MSE.

5. Models

Several classes of heavy-tailed time series processes have been proposed in the literature, including various modifications of linear time series models. But often these processes do not have closed-form densities for their marginal distributions. By contrast, in Sections 5.1 and 5.2, we present two nonlinear time series models each with a Pareto Type III marginal distribution and with known survival function. The first model has an AR(1) minification structure. i.e. Yt=Kmin(Yt−1,Et), where K>1 is an appropriate constant. The innovations {Et} are usually positive-valued and assumed to be independent of Yt−1. Explicit marginal distributions have been obtained when {Et} follows an exponential distribution or a Weibull distribution; see, e.g. Han [15] for a brief review. Other examples include the logistic processes by Arnold [2], Sim [38] and the semi-Pareto process by Pillai [26]. The second nonlinear model presented in Section 5.2, is a nonlinear Pareto moving average (MA) as a direct modification of the process defined in Section 5.1.

5.1. Nonlinear Pareto AR(1) process

A time series process {Yt,t∈Z} is called a first-order Pareto Type III autoregressive (PAR(1)) process if it has the form (15) Yt=min(p−1/αYt−1,11−WtZt),

where 1/0 is interpreted as +∞, and where {Wt} denotes a sequence of i.i.d. Bernoulli(p) r.v.'s, which is independent of Yt and Zt. If the initial distribution is Y0∼P(σ,α) then the process is strictly stationary.

Model (15) was proposed by Yeh [41] and by Yeh et al. [42]. The motivation behind the PAR(1) model was three-fold: partly as an alternative to the assumption of marginal Gaussianity which underlies the Gaussian linear stationary time series model, partly as a model for correlated positive r.v.'s with P(σ,α) marginal distributions; but chiefly as a simple point-process with which to analyse non-Poisson series of events. The economic interpretation of the PAR(1) model is explained by an example in Yeh [41].

The sample path of Yt exhibit sudden large peaks followed by sudden drops to lower levels; cf. Ferreira [14]. Explicit expressions for ρ(G)(ℓ) are available only in the case ℓ=1; Carcea and Serfling [6]. In general, expressions ρ(G)(ℓ) for ℓ>1 seem intractable for PAR(1) processes. Note that the PAR(1) process has three parameters to estimate: σ, α, and p. Initial estimates of σ and α can be obtained by equating suitable sample moments to their expectations. Then a variable metric numerical method can be used to estimate all three parameters simultaneously through a maximum likelihood function; see, e.g. Balakrishnan [3, Ch. 5] and Yeh et al. [42].

5.2. Nonlinear Pareto MA(1) process

One problem the PAR(1) model has is that it generates sample paths in which large values are followed by runs of decreasing values, with the runs having geometrically distributed lengths. The large values arise when Zt is included (i.e.  Wt=0) while the falling values stem from the deterministic part of (15) (i.e.  Wt=1). This behavioral characteristic is likely to limit the scope for practical application of the PAR(1) model. To overcome this deficiency, we propose a so-called first-order Pareto moving average (PMA(1)) process. Specifically, the Pareto Type III PMA(1) process is defined by (16) Yt=min(p−1/αZt−1,11−WtZt).

So in the first term on the right-hand side of (15), Yt−1 has been replaced by Zt−1. The mean and variance of {Yt,t∈Z} are easy to derive from the survival function F¯Yt(y)=1/(1+(y/σ)α). However, it is an open problem to obtain an explicit expression for the lag-one G-ACF for 1<α<2. Figure 1 shows two simulated sample paths of the PMA(1) process. Figure 1. Simulated sample paths of a PMA(1) process with Zt∼i.i.d.P(1,1.5) and n = 200: (a) p = 0.3; (b) p = 0.5.

6. Simulation study

We conducted three simulation experiments to evaluate the performance of the T(m) and T(G)(m) test statistics for sample sizes n=200, 400, and 1000. In all experiments, the innovations {Zt} are simulated as i.i.d. r.v.'s following a P(1,α) distribution with α=1.1,1.5, and 2.1. The ACFs of the partitioned (sub-sample) data in (13) and (14) are computed over m=5, 10 and 20 lags. For the models considered under the alternative hypothesis, the first 100 observations were discarded to avoid initialization effects.

6.1. Size

First, we consider the empirical size of each nonlinearity test statistic. For the linear model under H0, we selected an MA(1) process with parameters ψ1=0.5, and a zero-order MA (i.e. ψ1=0) process. Some test results are reported in Table 1. When α=1.1,1.5, m = 20, and n = 200, 400, the performance of T(m) is poor. In contrast, the empirical rejection frequencies of the test statistic T(G)(m) are close to the 5% nominal significance level for all combinations of α, m and n. Overall, the best performance of T(G)(m) seems to occur at the maximum lag m = 10. Finding the best value of m for the T(m) test statistic is more problematic. It appears that the asymptotic distribution of T(m) as given by Theorem 4.1 is sensitive to the value of m, and hence an overall recommendation is hard to give. Table 1. Empirical size (in %) of the T(m) and T(G)(m) test statistics at the 5% nominal significance level and computed from 1000 independent realizations of a linear MA(1) model of the form (5) with ψ1=0 and ψ1=0.5.

 	 	ψ1=0	ψ1=0.5	
 	 	α=1.1	α=1.5	α=2.1	α=1.1	α=1.5	α=2.1	
n	m	T	T(G)	T	T(G)	T	T(G)	T	T(G)	T	T(G)	T	T(G)	
200	5	3.7	4.2	2.5	5.2	2.2	4.7	4.2	5.1	4.8	5.6	4.6	4.6	
 	10	3.8	4.9	2.1	4.2	2.6	5.0	5.1	5.1	4.1	4.6	6.7	5.3	
 	20	1.6	3.5	1.7	4.1	1.4	5.9	5.1	3.7	2.8	4.1	4.2	5.1	
400	5	2.6	3.6	1.9	4.7	2.1	4.5	4.1	4.3	5.4	5.2	6.1	3.3	
 	10	3.1	3.6	3.3	4.5	2.4	5.6	3.9	4.7	4.5	5.9	5.3	4.6	
 	20	3.3	3.1	2.7	5.0	2.0	5.1	4.1	4.3	3.7	4.4	3.9	4.8	

6.2. Power

One way to evaluate the empirical power in relation to the size of a given test is to study the so-called size-power tradeoff curve as suggested by Davidson and MacKinnon [11]. Figure 2 shows tradeoff curves for the T(m) and T(G)(m) test statistics for n = 400, 1000, m = 10, 15, 20, α=1.5, and p = 0.33. For each curve, the horizontal axis shows the size computed for the data generating process that satisfies the null hypothesis (5). The vertical axis shows the empirical power for the PAR(1) and the PMA(1) processes, respectively. The upper right-hand corner of each graph corresponds to a critical value of zero. Both size and power are one at this point. The lower left-hand corner corresponds to a very large critical value, so large that the test statistic under study will never exceed it. Both size and power are 0 at this point. A test for which size always equals power has a size-power tradeoff curve by the 45° line (dotted line). Figure 2. Size-power tradeoff curves.

From Figure 2, we see that the T(G)(m) test statistic is systematically more powerful than T(m) for a given size of test. For m = 20, the size-power tradeoff curves of T(G)(m) are slightly more powerful than for m = 10 and m = 15, with curves further from the 45° line. The overall picture is completely different for T(m). Indeed, for the relevant size range (0, 0.15] the size of the test statistic T(m) exceeds its power for both nonlinear models with curves below the 45° line. A test statistic with this type of size-power tradeoff function is called a ‘biased test’. It is clearly undesirable to use a test which is more likely to reject the H0 when its true than when its false. In fact, biased tests are not of very much use in practice; see, e.g. Kendall and Stuart [20, Ch. 23] for a discussion on this topic.

As expected, for n = 1000 the size-power curves of T(G)(m) indicate that the test is slightly more powerful than for n = 400 for all values of m. On the other hand, the bias in T(m) does not seem to decrease as the sample size increases. Numerous other experiments, for which results are not shown, produced similar ordering of T(m) and T(G)(m). The above observations remain valid when the uncentered G-ACF is replaced by its centered counterpart.

6.3. P-values

Another way to explore the performance of both test statistics is to establish their statistical significance via p-value computation. Table 2 shows averages of observed p-values, i.e. the fraction of times that ∨ℓ=1m|Yˆ(ℓ)| was greater than the value of a particular test statistic. Hence the smaller the p-value, the stronger the evidence to reject the linear model (5). We see that there is a clear distinction between the results when n = 200 and n = 400 for all values of α. High p-values for the test statistic T(m) in the third column indicate that the process under study is linear when n = 200. By contrast, when n = 400 all p-values of this test statistic with the PAR(1) model indicate that the H0 of linearity should be rejected at a 5% nominal significance level. The p-values for T(G)(m) with the PMA(1) model as an alternative indicate that the null hypothesis should be rejected. So, it seems that T(m) is sensitive to parameter values α and sample size n. Similar test results were obtained for other values of the Bernoulli probability parameter p. In particular, we noted that for p↑1 there is increasing dependence in the simulated time series resulting in lower p-values. Indeed, for parameter values p near to 0 the PMA(1) model reduces to Yt=Zt with probability 1 while for p = 1 the PAR(1) model becomes Yt=Yt−1 with probability 1. Table 2. Averages of p-values for the T(m) and T(G)(m) test statistics.

 	 	n = 200	n = 400	
 	 	PAR(1)	PMA(1)	PAR(1)	PMA(1)	
α	m	T	T(G)	T	T(G)	T	T(G)	T	T(G)	
1.1	5	0.041	0.015	0.038	0.018	0.017	0.004	0.015	0.005	
 	10	0.045	0.016	0.043	0.020	0.018	0.004	0.018	0.005	
 	20	0.059	0.025	0.062	0.027	0.023	0.006	0.020	0.007	
1.5	5	0.056	0.012	0.057	0.012	0.023	0.005	0.057	0.012	
 	10	0.061	0.014	0.060	0.014	0.024	0.006	0.060	0.014	
 	20	0.076	0.020	0.075	0.019	0.028	0.007	0.075	0.019	
2.1	5	0.055	0.020	0.054	0.019	0.022	0.008	0.054	0.019	
 	10	0.060	0.022	0.059	0.022	0.023	0.009	0.059	0.022	
 	20	0.078	0.030	0.075	0.029	0.027	0.011	0.075	0.029	
Note: Embolded entries show values larger than a nominal rejection level of 5%.

7. Two applications

7.1. Fire insurance losses

In the actuarial context, Pareto distributed individual loss r.v.'s are utilized to model extreme loss and risky types of insurance coverage. In applications involving extreme value theory, the Danish fire data set has been the subject of many studies; see, e.g. McNeil [24] and Resnick [28]. The data set consists of inflation-adjusted losses to buildings, losses to contents, and losses to profits expressed in millions of Danish Kroner (1985 prices) and covering the period 1980–1993.4 Although the full data set consists of 2,493 observations, we confine our attention to strictly positive components, resulting in n = 517 observations. Figure 3 shows plots of the three series. Some summary statistics are also presented. The heavy-tailed nature of each series is clearly suggested by the high kurtosis measure as well as by the extreme maximum observed values. Figure 3. Danish fire data; n = 517.

Table 3 shows observed p-values of the test statistics when m = 5, 10, 20 with three fixed values of the tail index α. Several features are noteworthy. First, for the variables ‘Buildings’ and ‘Contents’ all observed p-values indicate that the null hypothesis (5) should not be rejected by T(m) and T(G)(m) at the 5% nominal significance level. However, the results are markedly different for the variable ‘Profits’ where the p-values of T(G)(m) are below the 5% nominal significance level regardless of the value of α while the p-values of T(m) yield no ground for rejecting H0. Furthermore, we see that for both test statistics the reported p-values are not sensitive to the value of α. So, for the variable ‘Profits’ the p-values of the test statistic T(G)(m) provide strong evidence against the null hypothesis that the data under study are generated by a linear MA model with heavy-tailed Pareto-type innovations. On the other hand, for the variables ‘Buildings’ and ‘Contents’ there is very little evidence arguing against the null hypothesis (5). Table 3. Fire insurance losses: P-values for the test statistics T(m) and T(G)(m).

 	 	α=1.1	α=1.5	α=2.1	
Series	m	T	T(G)	T	T(G)	T	T(G)	
Buildings	5	0.267	0.069	0.271	0.072	0.243	0.076	
 	10	0.272	0.077	0.250	0.079	0.270	0.065	
 	20	0.257	0.087	0.250	0.066	0.272	0.077	
Contents	5	0.976	0.118	0.978	0.122	0.976	0.125	
 	10	0.977	0.116	0.979	0,124	0.973	0.123	
 	20	0.972	0.121	0.976	0.119	0.972	0.122	
Profits	5	0.112	0.031	0.119	0.024	0.106	0.025	
 	10	0.116	0.037	0.112	0.034	0.092	0.041	
 	20	0.111	0.040	0.093	0.040	0.096	0.038	

7.2. Ethernet traffic

As a second application of the nonlinearity test statistics, we consider two time series of byte-rate variation in Ethernet traffic data. Many empirical studies have shown that Ethernet traffic appears to possess the properties of a self-similar process [23,40], and have subsequently led to insights into the burstiness properties of individual transport protocol connections that constitute teletraffic data; see Paxon and Floyd [25]. For instance, the distribution of the byte-rate is heavy-tailed and close to a Pareto distribution; see, e.g.Resnick [29] and the references therein.

Here, we focus on the so-called BellCore Ethernet traffic measurements described and analysed by Willinger et al. [40].5 Figure 4 shows the byte-rate time series for fixed time intervals of 1 (n = 3,142) and 10 (n=314) seconds. Some summary statistics are also presented. The distribution function of both series is right-skewed. The extreme maximum values indicate the heavy-tailed nature of the data. Figure 4. Time series plot of arrival rate aggregated over 1 second and 10 second intervals.

Table 4 contains the p-values of the T(m) and T(G)(m) test statistics. For the variable ‘1 sec’ the p-values for T(G)(m) provide strong evidence to reject the null hypothesis (5) with values smaller than the 5% nominal significance level. By contrast, the p-values of T(m) indicate that there is no numerical evidence that the process generating this variable is nonlinear. It is interesting to compare this results with test results obtained for the aggregated variable ‘10 secs’. Here all observed p-values for T(G)(m) are above 0.05, indicating that the underlying process is linear. In comparison, some p-values of the test statistic T(m) are smaller than 0.05 and some are slightly larger. Thus, in this case there is not enough evidence for the presence of nonlinear effects in the data. Table 4. Ethernet traffic data: P-values for the test statistics T(m) and T(G)(m).

 	 	α=1.1	α=1.5	α=2.1	
Series	m	T	T(G)	T	T(G)	T	T(G)	
1 sec	5	0.291	0.028	0.262	0.029	0.277	0.030	
 	10	0.289	0.039	0.285	0.048	0.248	0.036	
 	20	0.249	0.036	0.286	0.048	0.271	0.037	
10 secs	5	0.047	0.247	0.052	0.256	0.052	0.231	
 	10	0.061	0.233	0.049	0.265	0.051	0.253	
 	20	0.046	0.254	0.052	0.254	0.071	0.246	

The above observations on the variable ‘1 sec’ using T(G)(m) are supported by previous time series studies of Ethernet traffic data. For instance, Resnick [29] noted that empirical graphical techniques provided an indication of nonlinearity in a data set of inter arrival times of ISDN D-channel packets. Moreover, he expressed the need for ‘tools’ to cope with nonlinearities. Another study is by Chandra et al. [8]. They concluded that byte-rate data are nonlinear and can be adequately modeled by a so-called self-exciting threshold AR process. Our test results for the variable ‘10 secs’ are somewhat mixed and should be interpreted cautiously since aggregation over a fixed time interval may well destroy nonlinear features in the data.

Evidently, one of the clearest messages conveyed by the p-values reported in Tables 3 and 4 is that the overall conclusion about the null hypothesis may change depending on the use of T(m) or T(G)(m). Thus, one should not rely completely on the performance of one test statistic for assessing the appropriateness of time series processes with heavy-tailed innovations. However, only adopting T(m) may result in type I or type II errors. Recall that the difference between T(m) and T(G)(m) stems from the way the sub-sample ACFs are calculated. A test based on ρ(G)(ℓ) seems to detect a more intricate dependence structure (see, e.g.[37]) than a test based on ρ(ℓ), and hence is certainly worth applying in practice.

8. Discussion

The G-ACF is recommended in case there are indications that the process is heavy-tailed with Pareto-type nonlinearities. More precisely, the proposed Gini-based test statistic T(G)(m) can aid in detecting nonlinear features in heavy-tailed data. Simulation experiments showed that the size and power of the test T(G)(m) are good in comparison with the test statistic T(m). In fact, the size distortion of T(m) is quite severe and the test also seems biased for the sample sizes under consideration. Also, we have seen that the test statistic T(G)(m) provides evidence of nonlinear features in two empirical data sets. Of course, fitting PAR and/or PMA models to the data would be an interesting next step. At this point, however, we would like to stress that the empirical illustrations are only the beginning of a further study, not an in-depth analysis.

Several possible extensions to this paper can be pursued. For instance, the simulated series under the alternative hypothesis were driven by PAR(1) and PMA(1) processes. It is straightforward, to generalize these two processes to higher order schemes and bring them together to form, for instance, a nonlinear PARMA(1,1) process. However, it is not easy to calculate the survival function for these processes and their corresponding autocorrelation structure. Of course, other forms of heavy-tailed nonlinear processes are also worthy to consider.

Finally, recall that the Gini autocorrelations are generally not symmetric in the variables Yt and F(Yt−ℓ). For the Gini correlation, measuring the association between two r.v.'s X and Y, Sang et al. [34] proposed a symmetric measure based on a joint rank. It is straightforward to extend this measure to time series variables. Using this extension, it is worthwhile to compare the performance of T(G)(m) with a test statistic based on sub-sample symmetric Gini autocorrelations. Alternatively, the arithmetic mean (γ(G1)(ℓ)+γ(G2)(ℓ))/2 and the geometric mean (γ(G1)(ℓ)γ(G2)(ℓ)) can be used as symmetric alternatives of the traditional G-ACF.

Acknowledgments

We thank Amit Shelef for helpful comments on an earlier draft of the paper. Also, we express our sincere thanks to an Associate Editor and two anonymous referees for their helpful comments and suggestions.

Notes

Disclosure statement

The author reports no conflict of interest regarding this paper.

Data and source code availability statement

The R software source code for the simulations and the analysed data are available at the author's website: https://www.jandegooijer.nl.

1 Let X, Xi, i≥1, be independent non-negative random variables with common distribution function FX(⋅). Let F¯ denote the tail distribution of X, i.e. FX(x)=P(X>x). Then F is said to be heavy-tailed if, for all ϵ>0, E(eϵX)=∞, or equivalently, if for all ϵ>0, eϵXP(X>x)→∞, as x→∞. Basically, F is heavy-tailed if its tail decreases more slowly than exponentially. The above definition was taken from Chen [9, Definition. 1.1.1]; see, also, Resnick [29], Resnick and Van den Berg [32] and Embrechts et al. [13].

2 The independent stable r.v.'s S0,α/2,S1,α,…,Sj,α can be computed by a method due to Ament and O'Neil [1]; see https://gitlab.com/s.ament/qastable for MATLAB routines.

3 The HHL bootstrap algorithm is included in the R–blocklength package.

4 The data set is included in the R-package fitdistrplus.

5 The data is often referred to as ‘BCAug89’ traces and is available at the Internet Traffic Archive site https://ita.ee.lbl.gov/html/contrib/BC.html
==== Refs
References

1 S. Ament and M. O'Neil, Accurate and efficient numerical calculation of stable densities via optimized quadrature and asymptotics, Stat. Comput. 28 (2018), pp. 171–185.
2 B.C. Arnold, Pareto Distributions, Int. Cooperative Publishing House, Fairland, MD, 1983.
3 N. Balakrishnan, Non-Gaussian Autoregressive-Type Time Series, Springer-Verlag, Berlin, 2021.
4 J. Beirlant, Y. Goegebeur, J. Segers, and J. Teugels, Statistics of Extremes: Theory and Applications, Wiley, New York, 2004.
5 M.D. Carcea, Gini autocovariance function used for time series with heavy-tail distributions, WIREs Comput. Stat. 10 (2018), p. e1428.
6 M.D. Carcea and R. Serfling, A Gini autocovariance function for time series modelling, J. Time Ser. Anal. 36 (2015), pp. 817–838.
7 E. Castillo, A.S. Hadi, N. Balakrishnan, and J.M. Sarabia, Extreme Value and Related Models with Application in Engineering and Science, Wiley, New York, 2005.
8 K. Chandra, C. You, G. Olowoyeye, and C. Thompson, Non-linear time series models of Ethernet traffic, in Proc. 24th Conference on Local Computer Networks, LCN99, 1999, pp. 164–171.
9 B. Chen, Heavy tails: Asymptotics, algorithms, applications, Ph.D. dissertation, Technical University Eindhoven, The Netherlands.
10 J. Cohen, S.I. Resnick, and G. Samorodnitsky, Sample correlations of infinite variance time series models: An empirical and theoretical study, J. Appl. Math. Stoch. Anal. 11 (1988), pp. 255–282.
11 R. Davidson and J.G. MacKinnon, Graphical methods for investigating the size and power of hypothesis tests, Manchester School Econ. Soc. Stud. Univ. Manchester 66 (1998), pp. 1–26. Available at https://russell-davidson.arts.mcgill.ca/articles/gmn-euro.pdf.
12 R.A. Davis and S.I. Resnick, Limit theory for the sample covariance and correlation functions of moving averages, Ann. Statist. 14 (1986), pp. 533–558.
13 P. Embrechts, C. Klüppelberg, and T. Mikosch, Modelling Extremal Events for Insurance and Finance, Springer-Verlag, Berlin, 1997.
14 M. Ferreira, On the extremal behavior of a Pareto process; An alternative for ARMAX modeling, Kybernetika 48 (2012), pp. 31–49.
15 L. Han, W.J. Braun, and J. Loeppky, Random coefficient minification processes, Stat. Pap. 61 (2020), pp. 1741–1762.
16 P. Hall, J.L. Horowitz, and B.-Y. Jing, On blocking rules for the bootstrap with dependent data, Biometrika 82 (1995), pp. 561–574.
17 S.K. Kattumannil, A class of estimators: A unifying tool towards the estimation of Gini index and its variants, Technical Report RM 707, Dept. of Statistics and Probability, Michigan State University, 2014.
18 S.K. Kattumannil, N. Sreelakshmi, and N. Balakrishnan, Non-parametric inference for Gini covariance and its variants, Sankhyā A-84 (2022), pp. 790–807.
19 R.W. Katz, M.B Parlange, and P. Naveau, Statistics of extremes in hydrology, Adv. Water Resour. 25 (2002), pp. 1287–1304.
20 M.G. Kendall and A. Stuart, The Advanced Theory of Statistics, Vol. II , 4th ed., Charles Griffin, London, 1979.
21 G.G. Koedijk, M.M.A. Schafgans, and C.G. de Vries, The tail index of exchange rate returns, J. Int. Econ. 29 (1990), pp. 93–108.
22 S.N. Lahiri, Resampling Methods for Dependent Data, Springer–Verlag, New York, 2003.
23 W.E. Leland, M.S. Taqqu, W. Willinger, and D.V. Wilson, On the self-similarity nature of Ethernet traffic (extended version), IEEE/ACM Trans. Netw. 2 (1994), pp. 1–15.
24 A.J. McNeil, Estimating the tails of loss severity distributions using extreme value theory, ASTIN Bull. 27 (1997), pp. 117–137.
25 V. Paxon and S. Floyd, Wide-area traffic: The failure of Poisson modeling, IEEE/ACM Trans. Netw. 3 (1996), pp. 226–244.
26 R.N. Pillai, Semi-Pareto processes, J. Appl. Probab. 28 (1991), pp. 461–465.
27 B. Radovanov and A. Marcikí, A comparison of four different block bootstrap methods, Croat. Oper. Res. Rev. 5 (2014), pp. 189–202.
28 S.I. Resnick, Discussion of the Danish data on large fire insurance losses, ASTIN Bull. 27 (1997), pp. 139–151.
29 S.I. Resnick, Heavy tail modeling and teletraffic data, Ann. Statist. 25 (1997), pp. 1805–1869.
30 S.I. Resnick, Why non-linearities can ruin the heavy-tailed modelers day, in A Practical Guide to Heavy Tails: Statistical Techniques for Analyzing Heavy Tailed Distributions, R. Adler, R. Feldman and M. Taqqu, eds., Birkhäuser Verlag, Boston, 1997, pp. 219–240.
31 S.I. Resnick, G. Samorodnitsky, and F. Xue, How misleading can sample ACFs of stable MA's be? (Very!), Ann. Appl. Prob. 9 (1999), pp. 797–817.
32 S.I. Resnick and E. Van den Berg, A test for nonlinearity of time series with infinite variance, Extremes 3 (2000), pp. 145–172.
33 A. Sandström, Asymptotic normality of linear functions of concomitants of order statistics, Metrika 34 (1987), pp. 129–142.
34 Y. Sang, X. Dang, and H. Sang, Symmetric Gini covariance and correlation, Canad. J. Statist. 44 (2016), pp. 323–342.
35 E. Schechtman and S. Yitzhaki, A measure of association based on Gini mean difference, Comm. Statist. Theory Methods 16 (1987), pp. 207–231.
36 A. Shelef, A Gini-based unit root test, Comput. Stat. Data Anal. 100 (2016), pp. 763–772.
37 A. Shelef and E. Schechtman, A Gini-based time series analysis and test for reversibility, Statist. Papers 60 (2019), pp. 687–716.
38 C.H. Sim, First-order autoregressive logistic processes, J. Appl. Probab. 30 (1993), pp. 467–470.
39 C. Vanderford, On symmetric and computationally efficient Gini correlations in several bivariate distributions, Ph.D. diss., University of Mississippi, 2022.
40 W. Willinger, M. Taqqu, W. Leland, and D. Wilson, Self-similarity in high-speed packet traffic analysis and modeling of ethernet traffic measurements, Statist. Sci. 10 (1995), pp. 67–85.
41 H.-C. Yeh, Pareto processes, Ph.D. diss, University of California Riverside, 1983.
42 H.-C. Yeh, B.C. Arnold, and C.A. Robertson, Pareto processes, J. Appl. Probab. 25 (1988), pp. 291–301.
43 V.M. Zolotarev, One-Dimensional Stable Distributions, Translations of Mathematical Monographs, Vol. 65, American Mathematical Society, Providence, RI, 1986.
