SDD : An R Package for Serial Dependence Diagrams

Detecting and measuring lag-dependencies is very important in time-series analysis. This study is commonly carried out by focusing on the linear lag-dependencies via the well-known autocorrelogram. However, in practice, there are many situations in which the autocorrelogram fails because of the nonlinear structure of the serial dependence. To cope with this problem, in this paper the R package SDD is introduced. Among the available approaches to analyze the lag-dependencies in an omnibus way, the SDD package considers the autodependogram and some of its variants. The autodependogram, deﬁned by computing the classical Pearson χ 2 -statistic at various lags, is a graphical device recently proposed in the literature to analyze lag-dependencies. The concept of reproducibility probability, and several density-based measures of divergence, are considered to deﬁne the variants of the autodependogram. An application to daily returns of the Swiss Market Index is also presented to exemplify the use of the package.


Introduction
Investigating the temporal dependence structure is of fundamental importance in time-series analysis. The autocorrelogram, which measures the strength of linear dependencies (autocorrelations) as a function of the time lags, has been one of the primary tools for exploring and testing serial dependence for many decades. However, in practice, there are many situations in which the autocorrelogram fails because of the nonlinear structure of the serial dependence. A typical example is given by the GARCH model, whose components are uncorrelated but dependent (see Diks 2009, p. 6253, for another example). Despite the fast advancement in nonlinear time series models, there are few tools which can explore the complex dependence structures in nonlinear time series, like the autocorrelogram does for the linear ones. Starting from these considerations, several diagrams have been recently proposed which are very similar, in aspect and intent, to the autocorrelogram but they are widely applicable to both linear and nonlinear time series (see Anderson and Vahid 2005, Bagnato, Punzo, and Nicolis 2012, Bagnato, De Capitani, and Punzo 2014, 2013a, and Zhou 2012. In this paper we present the R (R Core Team 2014) package SDD (Bagnato, De Capitani, Mazza, and Punzo 2015) which allows for the calculation and display of some of the dependence measures discussed above. Note that, although R is well-provided with functions (such as acf()) and packages (such as tseries, Trapletti and Hornik 2013), to analyze time series (see also McLeod, Yu, and Krougly 2007, Hyndman and Khandakar 2008, Gasparrini 2011, Cowpertwait and Metcalfe 2009, and Hyndman 2014 for an overview), it does not offer graphical tools of the type described above.
The paper is organized as follows. Sections 2, 3, and 4, review the theoretical foundations of the dependence measures implemented in the SDD package. The relevance of these diagrams is shown, via a real data set, in Section 5, and conclusions are finally given in Section 6.

The autodependogram
Let {X t } t∈N represent a strictly stationary and ergodic stochastic process. Moreover, let (X 1 , . . . , X n ) be an observed time series of length n from {X t } t∈N . To study the generic dependence of lag r, r < n, let us consider the n r = n−r couples {(X i , X i+r )} nr i=1 . Bagnato and Punzo (2010) propose to group these couples in a k×k contingency table, with marginal sets of adjacent intervals C . In detail, assume that p (r) u+ = P X t ∈ C (r) u > 0 for all u, p (r) +v = P X t ∈ D (r) v > 0 for all v, and P X t ∈ k u=1 C (r) u = P X t ∈ k v=1 D (r) v = 1. Moreover, define p (r) uv = P (X t , X t+r ) ∈ C (r) . The statistical hypothesis of interest is with θ +v . This hypothesis can be tested using the Pearson χ 2 statistic where p uv , respectively. In particular, it results that p  it is shown by simulations that, also in the serial context, the large sample null distribution of δ r is well-approximated by the χ 2 with (k − 1) 2 degrees of freedom. This fact allows to test the null hypothesis of independence for lag r using δ r as test statistic: denoting with χ 2 [η;q] the q-quantile of the χ 2 distribution with η degrees of freedom, the null hypothesis is rejected at level α if δ r > χ 2 [(k−1) 2 ;1−α] .
In conformity with the autocorrelogram (ACF), the diagram obtained by plotting δ r as a function of the time lags r, r = 1, . . . , l, is called autodependogram (ADF) by . The autodependogram may be applied to time series with missing data and can be used for model diagnostic checking of nonlinear models (Bagnato and Punzo 2013).
To completely specify the test statistic δ r it is necessary to define the partitions {C u } k u=1 and {D v } k v=1 . Following , the so-called equifrequency interval, which assigns equal frequencies to each interval C u and to each interval D v , will be adopted. Using the equifrequency intervals, only the value of k must be selected. As a default in the SDD package, k is chosen such that k = min {k s , k p } with k s = n l 5 1 2 and k p = 2 11 10 where · denotes the floor function while z 1−α stands for the (1 − α)-quantile of the standard normal distribution (see  to find out more on the motivation of this rule and for simulation results confirming its validity). Note that, based on (3), k is a function of n.
In particular, if n → ∞, then also k → ∞ (although with a lower rate of divergence) and the limiting null distribution of δ r , conveniently standardized, will be a standard normal (Morris 1975). As suggested by the discussion in Mann and Wald (1942), this asymptotic normality can be heuristically justified by the fact that the χ 2 distribution tends to the normal when the degrees of freedom diverge. However, for low values of k, the normal approximation of the χ 2 distribution is poor and, since rule (3) provides small values of k even when n is large (as an example, for n = 100 and n = 1000 we obtain k = 4 and k = 7, respectively), the SSD package always uses the χ 2 distribution.
Being k fixed and equal in each lag r, the same level-α critical value χ 2 [(k−1) 2 ;1−α] can be used for each of the l tests of lag-independence. This allows a horizontal line at height χ 2 [(k−1) 2 ;1−α] to be added to the autodependogram. Analogously with reference lines used in autocorrelograms, this line (hereafter level-α critical line or, simply, critical line) demarcates acceptance and rejection regions for each lag.

Some considerations
As documented in Bagnato et al. (2014), although the autodependogram is very similar in aspect to the autocorrelogram, the latter graphically represents the strength of the linear lag-dependencies through the autocorrelations, while the former displays the evidence of the presence of linear/nonlinear lag-dependencies. To clarify this distinction, it is sufficient to observe the different behavior of these diagrams when n diverges. In this case, the level-α critical lines of the autocorrelogram collapse to 0 while its bars tend toward the true values of the autocorrelations; on the other hand, the level-α critical line of the autodependogram remains fixed and the bars diverge with n under the alternative, and converge to zero under the null. Another important difference is that, the autodependogram is able to capture various and general serial dependence structures while the autocorrelogram points out only linear relationships. This feature of the autodependogram motivates the adjective omnibus used in . Obviously, this generality is paid for in terms of power, and hence in the resulting descriptive ability under certain dependence structures. For example, the autodependogram will be less informative with respect to the autocorrelogram when the lagdependencies are linear, or with respect to a representation using rank correlation statistics when the lag-dependencies are monotonic.

Normalizations
Lags may be ranked according to the evidence of their dependence, as measured by the autodependogram bars. Nevertheless, being the not normalized χ 2 statistic, these bars do not reveal clearly the strength of this evidence and do not allow comparisons between different time series. A better alternative, proposed in Bagnato et al. (2012, Section 6), consists in substituting δ r with the Cramer contingency coefficient (C-ADF) Because the denominator in (4) is equal to the maximum value that the numerator can assume, this statistic can be considered as a "normalized dependence measure" assuming values between 0 (independence) and 1 (maximum dependence). It is interesting to note that, through the application of the transformation in (4) to the critical value χ 2 [(k−1) 2 ;1−α] , the (different) critical values for the tests based on ν r are obtained. These values allow to superimpose on the diagram of ν r , r = 1, . . . , l, a level-α critical line which, in this case, will not be horizontal but increasing in r. However, as it is well-known, ν r attains values near to its maximum sporadically (see Bagnato et al. 2012, Section 6).
A further normalized diagram representing the evidence of the presence of dependence can be obtained substituting δ r with where p r = 1−F (k−1) 2 δ r denotes the p value associated with δ r and F p is the cumulative χ 2 distribution with p degrees of freedom. Also in this case a level-α critical line at height 1 − α can be added to the resulting diagram. Nevertheless, p r unbalances the two decision regions by mapping the rejection one only on (1 − α, 1]. To avoid the last problem, it is possible to introduce a particular monotonic decreasing transformation p * r = g ( p r ), on [0, 1], such that g (α) = 1/2. In this way the acceptance and the rejection regions are mapped on [0, 1/2] and (1/2, 1], respectively. For example: The transformation in (6) preserves the autodependogram's ability to detect autodependencies.

The reproducibility probability autodependogram
A particularly appealing alternative diagram, proposed by Bagnato et al. (2014), stems from the reproducibility probability (RP; Goodman 1992) representation of the χ 2 -test of independence.
Under local alternatives of dependence of lag r, Bagnato et al. (2014) show that the large sample distribution of δ r is generally well-approximated by a noncentral χ 2 with (k − 1) 2 degrees of freedom and noncentrality parameter This result allows the definition of the approximated power function where F η,q denotes the cumulative noncentral χ 2 distribution with η degrees of freedom and noncentrality parameter q. The function π α,n (·) associates the probability of rejecting the null hypothesis of lag-independence to each possible value of δ r . It is strictly increasing and continuous since F η,q (x) is strictly decreasing in q for each x > 0 and η ∈ N (see Johnson, Kotz, and Balakrishnan 1995, p. 444). Then, the χ 2 -test of lag-independence is strictly unbiased and it can be re-defined using the RP-testing technique, according to Martini (2008).
To understand this technique (see also De Capitani 2013, for a simple introduction) it is useful to represent the test by the so called critical function: Representation (8) highlights that the test is a Bernoulli random variable and so the random nature of the statistical test results.
Given n, denote with δ * r the true but unknown value of the noncentrality parameter. The value assumed by the power function in correspondence to δ * r , that is π α,n (δ * r ), is the true power of the test and it coincides with the only unknown parameter of the Bernoulli random variable Ψ (r) α . Thus, the randomness of the test is completely described by π α,n (δ * r ) and this makes the true power the key element to evaluate the reliability of the test result. The true power is also known as reproducibility probability (RP; Goodman 1992) since it is the probability of obtaining a rejection of the null hypothesis in subsequent and identical replications of the test.
From a practical point of view, the RP is unknown since δ * r is unknown too. The most natural way to obtain an RP estimator is to plug δ r (which is an estimator of δ r ) into the power function obtaining π r = π α,n δ r . A further possible solution considers the estimator which satisfies P δ • r ≤ δ * r = 1/2. Martini (2008) shows that π r , differently from π r , allows the replication of the test (8) with the following critical function: In practice: "when the probability to reject the null hypothesis is estimated to be greater than that to accept it, then the null is rejected" (Martini 2008, Remark 2). The usefulness of the RP-based decision rule is also due to the fact that the RP estimate provides a rational indicator of the reproducibility of the test outcome. It is worth noting that, the equivalence between the RP-based decision rule and the classical one holds for all the values of α and n since their effect is implicit in the definition of the power function (Martini 2008 andDe Capitani andMartini 2011).
The RP-autodependogram (RP-ADF) of Bagnato et al. (2014) is obtained by substituting δ r with π r = π α,nr δ • r , r = 1, . . . , l, where π r can be interpreted as reproducibility probability of the level-α independence test at lag r. This additional information is very precious since it clearly quantifies the evidence of the presence of lag dependence. Moreover, for the values of α commonly used, the RPautodependogram bars are almost normalized since the values of π r range in [α, 1]. Thanks to the RP-testing result in (10), this diagram is also endowed with a level-α critical line, at height 1/2, which balances the acceptance and the rejection regions by mapping them on [α, 1/2] and (1/2, 1], respectively. Note that the sample size n r and the level α do not influence the height of the critical line but they have a great impact on the bars of the RPautodependogram. In more detail: the greater the value of α, the higher the bars of the RP-autodependogram (ceteris paribus); the greater the sample sizes n r , the higher the bars of the RP-autodependogram (ceteris paribus). Note that this behavior is in agreement with the influence of n r and α on the power of the test (Bagnato et al. 2014).

The divergence-based autopedendograms
In addition to the requirements introduced at the beginning of Section 2, suppose that X 1 has an absolutely continuous density g with support IR. Moreover, assume that (X 1 , X 1+r ) has an absolutely continuous joint density f r with support IR 2 . For a motivation about these additional assumptions, see Bagnato et al. (2013a). The presence of dependence for lag r can be checked by testing the statistical hypothesis To evaluate the discrepancy between f r (x, y) and g (x) g (y), with the aim to obtain a test statistic for the testing problem (12), several divergence functionals can be considered (see, e.g., Diks 2009). By analogy with Bagnato, De Capitani, and Punzo (2013b), all the functionals considered in SDD have the form where D is a real-valued function. An example of (13) is the generalized entropy of Tsallis It is easy to note that: ∆ (r) 1/2 coincides with the Hellinger metric, ∆ 1 is the Kullback-Leibler (KL) divergence, while ∆ (r) 2 can be interpreted as the "continuous counterpart" of the Pearson χ 2 statistic adopted in Section 2. A further intuitive dependence measure is the L 1 -distance: It can be interpreted as the "continuous counterpart" of the well-known Mortara dependence index (Mortara 1922). Similarly, in line with Rosenblatt (1975) and Skaug and Tjøstheim (1993), the following functionals can be respectively introduced: Naturally, all the aforementioned functionals are sensitive to departures from independence and, consequently, the testing problem (12) can be solved using an estimator ∆ (r) of ∆ (r) as test statistic with the null hypothesis rejected for large values of ∆ (r) .
Once chosen ∆, the corresponding autodependogram, hereafter simply denoted as ∆-autodependogram (∆-ADF), is obtained by substituting p r in (6) with the p value q r related to ∆ (r) : The SDD package implements the ∆-ADF based on the eight functionals ∆ 1 seems to be, among the eight aforementioned functionals, the best performer. Several methodologies are proposed in the literature to implement ∆ (r) ; the differences among the various approaches stem from the way: (i) the densities f r and g are estimated, (ii) the integral in (13) is computed, (iii) the p values q r , r = 1, . . . , l, are obtained. Among the alternatives considered by Bagnato et al. (2013b), the Gaussian kernel density estimator, an approximated numerical integration (similarly to Granger, Maasoumi, and Racine 2004), and the permutation approach, appear to be the best solutions for (i), (ii) and (iii), respectively. Details on these settings are given below.

Gaussian kernel density estimator
The Gaussian kernel (GK) estimator for the univariate density g is where K h (x; X i ) = (2πh 2 ) −1/2 exp − 1 2 h −1 (x − X i ) 2 and h > 0 is the bandwidth. Simi-larly, the GK (product) estimator for the bivariate density f r is where K h is defined as in (19).
To apply the GK density estimators (19) and (20), a value for the bandwidth h needs to be chosen. With this aim, in the literature, several data dependent procedures have been proposed. Examples are the Silverman's rule of thumb and the likelihood cross-validation method (see Silverman 1986, p. 45 and Section 3.4.4, respectively). The bandwidth obtained with the latter procedure, denoted with h LCV , is particularly useful in this context because, as observed in Granger et al. (2004, p. 654), it produces optimal density estimators according to the KL-criterion. In detail, working on the estimator in (19), h LCV is obtained by minimizing the cross-validation statistic is the estimated density constructed from all the data points except X i . The obtained value of h LCV is then used also in (20).
However, it is well-known that a bandwidth which is optimal for estimation is usually suboptimal for testing. In particular, although h LCV suffices to establish consistency of the test statistics, this choice could not be optimal in terms of the power of the resulting tests. As observed by Anderson, Hall, and Titterington (1994), in testing procedures a relative oversmoothing may be appropriate for some dependence functionals (test statistics). Nevertheless, the simulation results of Bagnato et al. (2013b) highlight that when the GK density estimator is adopted to define the estimator of ∆ (r) , the use of h LCV is appropriate and, then, no oversmoothing is applied in SDD.

Estimation of the dependence functional
Simulation results in Bagnato et al. (2013b) indicate that a good solution for the estimation of ∆ (r) is where f r and g are the GK estimates of f and g, and where, following the default setting of the sm package (Bowman and Azzalini 2014), with a = x (n) − x (1) /4 and x (1) (x (n) ) denoting the minimum (maximum) observed value. The valuesỹ j are defined in the same way. Note that (21) is an approximation of based on the 100 × 100 grid of equally spaced values {(x i ,ỹ j ) : i, j = 1, . . . , 100}.

Computing p values
Among the various proposals to compute the p value q r , the simulation studies of Bagnato et al. (2013b,a) suggest that a permutation approach represents a good compromise between simplicity and performance. It exploits the fact that, conditionally on the observed data x 1 , . . . , x n , each of the possible n! permutations is equally likely under the assumption of serial independence (see Diks 2009). In detail, let ∆ (r,0) denote the value assumed by ∆ (r) for the observed data. Analogously, let ∆ (r,b) be the dependence functional estimate obtained from a random permutation of the original data, with b = 1, . . . , B. Under the assumption of serial independence, ∆ (r,1) , . . . , ∆ (r,B) are equally likely and the p value can be defined as in Diks and Panchenko (2007): where L is defined as follows. Let Z = # ∆ (r,s) : ∆ (r,s) = ∆ (r,0) ; s = 0, 1, . . . , B ≥ 1 denote the number of ties plus one. If Z = 1 then L = 1, while if Z > 1 then L is drawn from the discrete uniform distribution on {1, . . . , Z}. The above procedure for computing the p value leads to a randomized test having an exact level α if the null is rejected whenever q r ≤ α and 0 < α = c/(B + 1) < 1 for some integer c.

Testing independence on more than one lag simultaneously
The SDD package also allows to test independence on a set of lags, say R, specified by the user. Such a procedure is typically used to test the more general hypothesis of serial independence (cf. Diks 2009). Among the available procedures (see Bagnato et al. 2013b, for some examples), the SDD package handles the "multiple-lag testing problem" almost surely ∀ r ∈ R in two simple ways.
A first way, the most common, consists in building a "Portmanteau test" (P-test); when the ∆-ADF is considered, the test statistic is Q R = r∈R ∆ (r) . In this case, recalling (23), the p value is computed as where Q (s) R = r∈R ∆ (r,s) . If the user adopts δ r instead of ∆ (r) (similar reasoning holds for the C-ADF), the asymptotic distribution of Q R = r∈R δ r , used to compute the p value, is a χ 2 with (k − 1) 2 card (R) degrees of freedom (see Punzo 2010, 2012, for details), where card (R) denotes the cardinality of R.
A second simpler way consists in using a "Simultaneous test" (S-test); it is based on the adjustment of the p values q r , r ∈ R, in order to take a decision about the simultaneous independence for the lags in R. A lot of methods exist for adjustment (see, e.g., Wright 1992); the package SDD implements all the methods allowed by the function p.adjust() of the stats package. As an example, the famous Bonferroni adjustment leads to the following decision rule: "do not reject H R 0 if all the single tests, for lags belonging to R, do not reject at level α/card (R)".
As can be inferred from the simulations of Bagnato et al. (2013b), the P-test seems to be a better choice. A further alternative, whose performance is comparable to the P-test, would be the use of the multiple-lag procedure described in Bagnato et al. (2013b). However, because the computational burden required is substantial, this technique is not considered here.

Package description and illustrative example
Package SDD is developed in an object-oriented design, using the standard S3 paradigm. Its main function, ADF(), allows for the computation and the plot of all the different types of serial dependence diagrams illustrated in the paper and it returns an object of class 'SDD'. Its arguments, along with their description, are listed in Table 1.
To illustrate the use of the package, the SMI dataset included in the SDD package and already analyzed in Bagnato et al. (2014) and Bagnato and Punzo (2013) is considered. Data consist of n = 660 daily returns of the Swiss Market Index spanning the period from 2009-08-12 to 2012-03-06 (the share prices used to compute the daily returns are downloadable from http://finance.yahoo.com/). To begin the analysis, data are loaded by R> library("SDD") R> data("SMI", package = "SDD") In the financial context, a first glance to the serial dependence can be obtained by considering the autocorrelogram of the squared raw series. The corresponding plot is obtained via the command R> res1 <-ADF(SMI^2, dtype = "ACF", main = "") It produces the plot in Figure 1(a), which shows linear dependencies, on the squared raw data, for the majority of the considered lags. The function acf(), of the stats package, is used and the autocorrelogram is displayed without the first bar referred to lag zero. The presence of serial dependence on the squared raw data is also corroborated by the autocorrelation-based tests, whose results are printed via the print() method dtype Autodependence function to be computed: "ACF" for the ACF; "ADF" (default) for the ADF; "CADF" for the C-ADF; "RPADF" for the RP-ADF; "DeltaADF" for the ∆-ADFs.

lag.max
Maximum lag l considered.
alpha Significance level α for the critical line.

num.clas
Value of k in the contingency table used for the χ 2 -based ADFs. If not specified, it is determined internally using the rule of thumb (3).
bandwidth Bandwidth h in the GK density estimator. If not specified, h LCV is used.

B
Number of permutations B used to compute p values for the ∆-ADF.
fres Alternative external function specifying the resampling method from the raw series; it is used to compute p values.
fdenest Alternative external function to perform estimation of g (x) and f r (x, y).
plot If TRUE (default), the diagram specified by dtype (and eventually by delta) is plotted.

R
The set R of Section 4.4.
Instead of considering directional (specific) types of lag-dependencies, we can act in an omnibus way by displaying the autodependogram variants. The following command R> res2 <-ADF(SMI, main = "") produces the ADF in Figure 1( 2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28 Also in this case, serial independence is rejected at any significance level by the Portmanteau test and at levels greater than 0.038 by the simultaneous test. The RP-ADF and the ∆ 1 -ADF, shown respectively in Figure 1(c) and Figure 1(d), are obtained with the commands R> res3 <-ADF(SMI, dtype = "RPADF", main = "") R> set.seed(1) R> ADF(SMI, dtype = "DeltaADF", main = "") In particular, the RP-ADF allows for additional interpretation; by focusing on the first two bars π 1 and π 2 , using the following command R> res3$res$vbar[1:2] [1] 0.8890191 0.8085375 we can conclude that the null hypothesis of lag-1 (lag-2) independence is rejected and we estimate that this hypothesis will be rejected with probability 0.88901 (0.80853) if the same test is performed on a different time series of the same length from the same generating process. On the other hand, the ∆ 1 -ADF in Figure 1(d) detects more lag-dependencies. From this point of view, it seems more similar to the ACF on the squared residuals. As stated in Bagnato et al. (2013a), it should be due to the higher power of the ∆ 1 -ADF among the implemented autodependograms.

Conclusions
In this paper we have presented the SDD package for the R environment. It provides several diagrams to analyze linear/nonlinear lag-dependencies for time series. Differently from other implemented and existing diagrams, like the autocorrelogram, this analysis is carried out in an omnibus way, without focusing on specific (directional) forms of lag-dependence. In other words, we can detect lag-dependencies that are more general than, for example, the linear or monotonic ones. The package also includes an illustrative data set on daily returns of the Swiss Market Index (SMI). Through the application of variants of the autodependogram to this data set, we have presented the usefulness of the SDD package for detecting and evaluating the lag-dependencies present in both raw data and residuals of a fitted model. We believe the availability of such diagnostic checks will be appreciated by other R users as well. Extensions to either "partial" or "cross" variants of the considered diagrams will be implemented in future versions of the package.