The BDS test of independence

Abstract

In this article, we describe and implement the Brock, Dechert, and Scheinkman (1987, Working paper) test of independence of the elements of a time series.

1 Introduction

The question of whether time-series data are independently distributed has a long history in econometrics. In Stata, runtest (see [R] runtest) is a nonparametric test of the hypothesis that the observations of the series occur in a random order by counting how many successive observations lie above or below a threshold such as the median of the series. Swed and Eisenhart (1943) provide exact critical values for this test. ¹ The interest in testing for independence in time series was reinvigorated in the early 1980s by the nonlinear dynamics and chaos literature with the goal of distinguishing deterministic systems from random systems. Brock, Dechert, and Scheinkman (1987) and Brock et al. (1996) introduced a test known as the Brock, Dechert, and Scheinkman (BDS) test for detecting dependence in time-series data based on the correlation dimension of the process. Although the popularity of the BDS test was initially derived from its link with deterministic chaos, it has proved to be useful as a residual diagnostic test because it can detect deviations from dependence in time-series data quite reliably. Furthermore, its asymptotic properties as a residual diagnostic are well understood. Extensive Monte Carlo results have proved the BDS test useful in relatively small samples (Brock, Hsieh, and LeBaron 1991).

The idea behind the BDS is fairly simple. If the data are independently and identically distributed (i.i.d.), then for a given distance ε, the probability that the difference between pairs of m-dimensional points is less than or equal to ε will be a constant. This probability is denoted C_m (ε) and is known as the correlation integral (Grassberger and Procaccia 1983). The intuition of the BDS test is that with i.i.d. data, C_m (ε) will simply be the product of the individual pairs, so C_m (ε) = C ₁(ε) ^m . However, the sample analogues of these quantities will not satisfy this condition exactly. The BDS test provides a formal statistical evaluation of the significance of this divergence.

Although the intuition of the BDS test is straightforward, calculating the BDS statistic is not easy, and it is even more challenging to perform the computations with the speed necessary to make bootstrap resampling feasible. To complement earlier implementations of the test, LeBaron (1997) describes a fast algorithm and provides C source code to implement the test. The initial component of LeBaron’s approach is based on an efficient sorting algorithm that, for a time series of length T, requires O(T log T) operations and hence pays dividends for large values of T . The main disadvantage of this fast algorithm is that it is not completely transparent.

By contrast, the approach to compute the BDS test described in this article offers transparency of operation because it does not require sorting of the data to estimate the correlation integral. Furthermore, the implementation in Mata makes the procedure accessible on any machine running Stata. This desire to facilitate ease of use does imply a speed penalty. Consequently, our variant is intrinsically O(T ²), which makes it considerably slower than LeBaron’s algorithm when applied to a lengthy time series. However, the absence of overhead costs makes the approach efficient for smaller datasets.

In addition to describing the command that implements the test, we also outline an elegant simplification for the computation of the variance of the BDS statistic.

2 The BDS test

The correlation integral introduced by Grassberger and Procaccia (1983) is a method for measuring the frequency with which temporal patterns are repeated in data. Given observations of a time series x_t, t = 1,…T, the m-history of the time series is

$x_t^m=\left(x_t,x_{t+1},x_{t+2},\dots ,x_{t+m-1}\right)$

A time series expressed in this form is described as being embedded of dimension m, with m = 1 representing the case where only the original time series is considered. The sample correlation integral for a given ε > 0 and embedding dimension m = 1 is

${\widehat{C}}_1=\frac{2}{n(n-1)}{\displaystyle \sum _{j=1}^{n-1}}\hspace{0.17em}{\displaystyle \sum _{k=j+1}^n}I\left(\left|x_j-x_k\right|\le \epsilon \right)$

1

where n = (T − m + 1) and I(S) is the indicator function taking the value 1 if S is a true statement and 0 otherwise. For ease of notation, the dependence of this quantity on ε is suppressed. Of course, if ε is chosen so that all pairs satisfy the condition, then ${\widehat{C}}_1=1$ ; if ε is chosen so that no pairs satisfy the condition, then ${\widehat{C}}_1=0$ . Consequently, the correlation integral has the interpretation of measuring spatial correlation. Often in practice, ε is set in terms of standard deviations of the data. ²

Strictly speaking, if there are T data points, the computation of ${\widehat{C}}_1$ with m = 1 can use all T elements of the sample to construct the (x_j, x_k ) pairs required in (1). However, if higher-order embedding dimensions are to be computed, it is important to ensure that the correlation integrals ${\widehat{C}}_1,\hspace{0.17em}.\hspace{0.17em}.\hspace{0.17em}.,{\widehat{C}}_m$ are each calculated on the same number of (x_j, x_k ) pairs. Consequently, the observed sample is always truncated to x ₁ ,…, x_n with n = (T − M + 1), where M is the maximum order of the embedding dimension.

The correlation integral ${\widehat{C}}_m$ is then given by

${\widehat{C}}_1=\frac{2}{n(n-1)}{\displaystyle \sum _{j=1}^{n-1}}\hspace{0.17em}{\displaystyle \sum _{k=j+1}^n}I\left(\left|x_j-x_k\right|\le \epsilon \right)$

Therefore, the mth-order correlation integral computes the joint probability

$\mathrm{Pr}\left(|x_j-x_k|<\epsilon ,|x_{j+1}-x_{k+1}|<\epsilon ,\dots ,|x_{j+m-1}-x_{k+m-1}|<\epsilon \right)$

For a given value of ε, define I_j,k = I(|x_j − x_k| ≤ ε). Then, ${\widehat{C}}_1$ and ${\widehat{C}}_m$ have simplified expressions

${\widehat{C}}_1=\frac{2}{n(n-1)}{\displaystyle \sum _{j=1}^{n-1}}\hspace{0.17em}{\displaystyle \sum _{k=j+1}^n}\hspace{0.17em}{I}_{j,k},\hspace{1em}{\widehat{C}}_m=\frac{2}{n(n-1)}{\displaystyle \sum _{j=1}^{n-1}}\hspace{0.17em}{\displaystyle \sum _{k=j+1}^n}\hspace{0.17em}{\displaystyle \prod _{r=0}^{m-1}}{I}_{j+r,k+r}$

2

in which occurrences of ε have been suppressed for representational simplicity.

The BDS test proceeds by noting that under the assumption of i.i.d. data, this probability will simply be the product of the individual probabilities for each pair if the observations are independent. Under this null hypothesis, it follows from (3) that

$\begin{array}{l}\text{E}\left({\widehat{C}}_1\right)=\frac{2}{n(n-1)}{\displaystyle \sum _{j=1}^{n-1}\hspace{0.17em}{\displaystyle \sum _{k=j+1}^n\text{E}}}\left(I_{j,k}\right)=\frac{2}{n(n-1)}{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{k=j+1}^n{C}_1}}=C_1\\ \text{E}\left({\widehat{C}}_m\right)=\frac{2}{n(n-1)}{\displaystyle \sum _{j=1}^{n-1}\hspace{0.17em}{\displaystyle \sum _{k=j+1}^n\hspace{0.17em}{\displaystyle \prod _{r=0}^{m-1}\text{E}}}}\left(I_{j+r,k+r}\right)=C_1^m\end{array}$

The BDS test provides a formal basis for judging the size of this error. Given a value for ε, Brock et al. (1996) defined the BDS statistic as

$W_{n,m}(\epsilon )=\sqrt{n}\left(\frac{{\widehat{C}}_m-{\widehat{C}}_1^m}{{\sigma }_{n,m}}\right)$

where the standard deviation σ_n,m is computed from the variance

${\sigma }_{n,m}^2=4\left\{{\beta }^m+2{\displaystyle \sum _{j=1}^{m-1}{\beta }^{m-j}}{\alpha }^{2j}+{(m-1)}^2{\alpha }^{2m}-m^2\beta {\alpha }^{2m-2}\right\}$

3

in which α and β are defined as

$\alpha =\frac{1}{n^2}{\displaystyle \sum _{j=1}^n\hspace{0.17em}{\displaystyle \sum _{k=1}^n{I}_{j,k}}},\hspace{1em}\beta =\frac{1}{n^3}{\displaystyle \sum _{i=1}^n\hspace{0.17em}{\displaystyle \sum _{j=1}^n\hspace{0.17em}{\displaystyle \sum _{k=1}^n{I}_{i,j}}}}{I}_{j,k}$

4

Brock et al. (1996) demonstrate that under the null hypothesis of i.i.d. data, this statistic converges in distribution to the standard normal

$W_{n,m}(\epsilon )\stackrel{d}{\to }N(0,1)$

for a given ε.

3 Simplifying the computation of the variance

It is clear from (3) that the expression for ${\sigma }_{n,m}^2$ can be regarded as a polynomial of degree m in β with coefficients that are functions of α; that is,

${\sigma }_{n,m}^2(\beta )=4\left\{{\beta }^m+2{\displaystyle \sum _{j=1}^{m-1}{\beta }^{m-j}}{\alpha }^{2j}+{(m-1)}^2{\alpha }^{2m}-m^2\beta {\alpha }^{2m-2}\right\}$

5

A straightforward calculation indicates that ${\sigma }_{n,m}^2\left({\alpha }^2\right)=0$ . The remainder theorem now asserts that (β − α ²) is a factor of (5), which when factored into the expression gives

${\sigma }_{n,m}^2(\beta )=\{\begin{array}{ll}4{\left(\beta -{\alpha }^2\right)}^2\hfill & m=2\hfill \\ 4\left(\beta -{\alpha }^2\right)g(\beta )\hfill & m\ge 3\hfill \end{array}$

where g(β) is the polynomial of degree (m − 1) with expression

$g(\beta )={\displaystyle \sum _{j=1}^{m-1}(2j-1)}{\beta }^{m-j}{\alpha }^{2(j-1)}-{(m-1)}^2{\alpha }^{2(m-1)}$

The simplest way to continue the analysis of ${\sigma }_{n,m}^2(\beta )$ is to replace (2j − 1) with the algebraically identical expression j ² − (j − 1)². With this substitution in place,

$g(\beta )={\displaystyle \sum _{j=1}^{m-1}{j}^2}{\beta }^{m-j}{\alpha }^{2(j-1)}-{\displaystyle \sum _{j=1}^{m-1}{(j-1)}^2}{\beta }^{m-j}{\alpha }^{2(j-1)}-{(m-1)}^2{\alpha }^{2(m-1)}$

The second summation is reindexed by replacing j − 1 with j to obtain

$\begin{array}{c}g(\beta )={\displaystyle \sum _{j=1}^{m-1}{j}^2}{\beta }^{m-j}{\alpha }^{2(j-1)}-{\displaystyle \sum _{j=1}^{m-2}{j}^2}{\beta }^{m-j-1}{\alpha }^{2j}-{(m-1)}^2{\alpha }^{2(m-1)}\\ ={\displaystyle \sum _{j=1}^{m-1}{j}^2}{\beta }^{m-j}{\alpha }^{2(j-1)}-{\displaystyle \sum _{j=1}^{m-1}{j}^2}{\beta }^{m-j-1}{\alpha }^{2j}\end{array}$

Both summations are now combined to give

$\begin{array}{c}g(\beta )={\displaystyle \sum _{j=1}^{m-1}{j}^2}\left\{{\beta }^{m-j}{\alpha }^{2(j-1)}-{\beta }^{m-j-1}{\alpha }^{2j}\right\}\\ =\left(\beta -{\alpha }^2\right){\displaystyle \sum _{j=1}^{m-1}{j}^2}{\beta }^{m-j-1}{\alpha }^{2(j-1)}\end{array}$

6

The case m = 2 can be incorporated into (6) to derive the final result

${\sigma }_{n,m}^2=4{\left(\beta -{\alpha }^2\right)}^2{\displaystyle \sum _{j=1}^{m-1}{j}^2}{\beta }^{m-j-1}{\alpha }^{2(j-1)},\hspace{1em}m\ge 2$

4 The bds command

The bds command calculates the BDS test as described in section 5.

4.1 Syntax

Before using the bds command and other similar Stata time-series commands, one must tsset or xtset the data so that the variable of interest is defined as a proper time series. The command syntax is

bds varname [ if ] [ in ] [ , m( integer ) eps( numlist )]

Note that varname may not contain gaps within the specified sample. varname can contain time-series operators. The command can be applied to one unit of a panel.

4.2 Options

The command supports the following options:

m( integer ) specifies the maximum embedding dimension, which evaluates sequences of length 2,…, m. The default is m(3).

eps( numlist ) provides a list of values of ε to be used in the comparison with the differences of pairs’ values. The test evaluates the number of absolute differences that are less than the epsilon value. If numlist values are not provided, six defaults are considered: the 70th percentile of the distribution of absolute differences, followed by 0.5, 1.0, 1.5, 2.0, and 2.5 times the standard deviation of the data. If values are provided, they must be given in ascending order.

4.3 Stored results

bds stores the following in r():

5 The algorithm

Given ε > 0, the quantity I_j,k has value 1 if |x_j − x_k| ≤ ε and 0 otherwise. The first step in the computation of the BDS statistic is to create a lookup table, which in Stata can be visualized as an upper triangular matrix of one-byte integers with T rows and T columns. The (j, k)th entry of this matrix is 1 if |x_j − x_k| ≤ ε and 0 otherwise. ³ The (k, k)th diagonal entry of this lookup matrix compares |x_k − x_k| with ε and therefore is 1 for all values of k.

5.1 First-order correlation integral

The sample estimate of the first-order correlation integral is the expected value of I_j,k taken over all values of j and k satisfying j < k ≤ n, or equivalently, the fraction of all pairs (x_j, x_k ) for which |x_j − x_k| ≤ ε. Put simply, this is the total count of all the nonzero entries in the upper triangle of the lookup matrix, excluding the main diagonal divided by the total number of distinct pairings, namely, n(n − 1)/2.

One way to think of this total count is as a sum of row counts r ₁ , r ₂ ,…, r_n , in which r_j denotes the sum of the jth row of the lookup matrix, excluding the contribution from the main diagonal. The alternative view is to think of the total count as a sum of column counts c ₁ , c ₂ ,…, c_n , in which c_j denotes the sum of the jth column of the lookup matrix, excluding the contribution from the main diagonal. In the former, r_n = 0, and in the latter, c ₁ = 0. Both strategies will give the same total count, but each will be composed of different partial counts. The explicit expressions for r_j and c_j are, respectively,

$r_j={\displaystyle \sum _{k=j+1}^n{I}_{j,k}},\hspace{1em}{c}_j={\displaystyle \sum _{k=1}^{j-1}{I}_{k,j}}$

The total number of counts, say, C, is therefore

${\displaystyle \sum _{j=1}^n{r}_j}=C={\displaystyle \sum _{j=1}^n{c}_j}\hspace{1em}\to \hspace{1em}2C={\displaystyle \sum _{j=1}^n\left(r_j+c_j\right)}$

and the estimated value of the first-order correlation integral is therefore 2C/n(n − 1).

Value of α

The usual decomposition of a double summation into diagonal and off-diagonal contributions allows α in the first expression in (4) to have representation

$\alpha =\frac{1}{n^2}\left({\displaystyle \sum _{j=1}^n{I}_{j,j}}+2{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{k=j+1}^n{I}_{j,k}}}\right)=\frac{n+2C}{n^2}$

Value of β

The calculation of the value of β begins by noting that (4) may be usefully rewritten as

$\beta =\frac{1}{n^3}{\displaystyle \sum _{j=1}^n\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^n{I}_{i,j}}}{I}_{j,k}\right)}$

This representation indicates that the value of β is constructed from a count of all triples (x_i, x_j, x_k ) for which I_i,jI_j,k = 1, including the possibilities i = j, k = j, or both. Suppose that this counting process has reached x_j . Then, the previous analysis has demonstrated that there are m = r_j + c_j + 1 data points within ε of x_j , including x_j itself. Clearly, the value of m will depend on j, but to maintain representational clarity, we suppress this dependence. Furthermore, suppose that these m data (which include x_j ) are x_{p 1} ,…, x_{p m} , ignoring again the dependence of p ₁ ,…, p_m on j.

Importantly, x_{p 1} ,…, x_{p m} all share the property that I_{j,p 1} = I_{p 1 ,j} · · · = I_{p m ,j} = I_{j,p m} = 1 for a fixed value of j. Consequently,

$r_j+c_j+1={\displaystyle \sum _{k=1}^m{I}_{j,p_k}}\hspace{1em}\to \hspace{1em}{\left({\displaystyle \sum _{r=1}^m{I}_{j,p_r}}\right)}^2={\left(r_j+c_j\right)}^2+2\left(r_j+c_j\right)+1$

However, by construction,

${\left({\displaystyle \sum _{r=1}^m{I}_{j,p_r}}\right)}^2={\displaystyle \sum _{r=1}^m{\displaystyle \sum _{s=1}^m{I}_{p_r,j}}}{I}_{j,p_s}={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^n{I}_{i,j}}}{I}_{j,k}$

In conclusion,

$\begin{array}{l}\beta =\frac{1}{n^3}{\displaystyle \sum _{j=1}^n\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^n{I}_{i,j}}}{I}_{j,k}\right)}=\frac{1}{n^3}\left\{{\displaystyle \sum _{j=1}^n{\left(r_j+c_j\right)}^2}+2{\displaystyle \sum _{j=1}^n\left(r_j+c_j\right)}+{\displaystyle \sum _{j=1}^{n}1}\right\}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\frac{1}{n^3}\left\{{\displaystyle \sum _{j=1}^n{\left(r_j+c_j\right)}^2}+2C+n\right\}\end{array}$

5.2 Higher-order correlation integrals

The sample estimate of the correlation integral at embedding dimension m ≤ M is the expected value of all quantities of type

${\displaystyle \prod _{p=0}^{m-1}{I}_{j+p,k+p}}$

7

taken over the n(n − 1)/2 possible values of (j, k) for which 1 ≤ j < k ≤ n, or equivalently, the fraction of all pairs (x_j, x_k ) for which expression (7) is 1, that is, I_{j + p,k + p} = 1 for all values of p satisfying 0 ≤ p < m. Given a pair (x_j, x_k ), the contributions made by that pair to the correlation integrals at each depth of embedding may be calculated simultaneously by noting that

${\displaystyle \prod _{p=0}^{m-1}{I}_{j+p,k+p}}=\left({\displaystyle \prod _{p=0}^{m-2}{I}_{j+p,k+p}}\right)\times I_{j+m-1,k+m-1}$

8

This means that the computational penalty involved in the calculation of higher-order correlation integrals is small, provided the contributions from each pair, say, (x_j, x_k ), at all requested levels of embedding are done simultaneously.

The matrix representation of the lookup table allows (8) to be computed efficiently. Because (k +p)−(j +p) = k −j = r is independent of p (and consequently the depth of embedding), the value of I_{j + p,k + p} is the entry in the matrix lookup table at row (j + p) and column (j + p) + r. In overview, the contributions to the estimated correlation integral at embedding depth m can be visualized as the sum of the products of m consecutive elements of the super diagonals of the matrix lookup table taken over the entire upper triangle (that is, excluding the main diagonal) for values of m from m = 2 to m = M. The estimate of the correlation integral is this sum divided by n(n − 1)/2.

6 Empirical applications

6.1 Sunspots

Sunspots are regions on the surface of the sun with magnetic field strengths thousands of times stronger than the earth’s magnetic field. They appear as dark spots on the surface of the sun and typically last for several days. The data, y_t , which are the annual averages of daily sunspot numbers from 1700 to 2017, were compiled by the Solar Influences Data Analysis Center in Belgium. The series is plotted in figure 1. This series not only is the longest directly observed index of solar activity but also has interesting time-series properties.

Figure 1.

Plot of average annual sunspots series from 1700 to 2017

The sunspot numbers have been of interest to climatologists who hypothesized a link between sunspot activity and climate change. White and Liu (2008) provide evidence that the solar cycle may be the trigger for El Niño and La Niña episodes from 1900–2005, suggesting that higher solar activity implies weaker and less frequent El Niño events. This view is controversial, and no generally accepted statistical link has been established. The so-called Maunder Minimum ⁴ between 1645 and 1715 was a period in which sunspots were scarce and the winters harsh, strongly suggesting a link between solar activity and climate change. The current view is that there has been no significant long-term upward trend in solar activity since 1700. This implies that rising global temperatures since the industrial revolution cannot be attributed to increased solar activity.

A simple application of the BDS test to the annual sunspot data confirms the strong rejection of the null hypothesis of i.i.d. data. The value of ε is set to one standard deviation of the data, and the maximum embedding dimension is 4.

In fact, the sunspot data have given rise to many attempts to use nonlinear modeling to capture the key features of the time series. In particular, many different self-exciting threshold autoregressive models have been proposed in the literature and applied to the sunspot data. Consider the following model, which is similar to those of Tong and Lim (1980) and Battaglia and Orfei (2005):

$y_t=\{\begin{array}{ll}{\alpha }_0+{\alpha }_1{y}_{t-1}+{\alpha }_2{y}_{t-2}+v_{1t}\hfill & \text{ if }{y}_{t-3}\le k\hfill \\ {\beta }_0+{\beta }_1{y}_{t-1}+{\beta }_2{y}_{t-2}+v_{2t}\hfill & \text{ if }{y}_{t-3}>k\hfill \end{array}$

9

Fitting the model using the Stata command threshold (see [TS] threshold) yields

The optimal threshold value returned by the search is k = 60.7, which is slightly larger than the values reported by Tong and Lim (1980) and Battaglia and Orfei (2005), although they use annual data over a shorter sample period in their work and include a much more complex dynamic structure. The autoregressive coefficients in both regimes have a somewhat similar pattern in terms of sign, although the difference in their sizes is slightly more marked in regime 2. These results suggest a rather complex dynamic pattern, although it is possible to conjecture that the large positive first-order autocorrelation coefficient is offset by a much smaller negative second-order autocorrelation coefficient, which implies that the process persists a little longer in regime 2.

The residuals from the threshold regression are plotted in figure 2. The first impression is that there is less structure in the residuals than in the original series, but the real question to be addressed is whether the residuals are now i.i.d. Application of the BDS, again with a maximum embedding dimension of 4 specified as one of the options, yields the following results:

Figure 2.

Residuals obtained from fitting the model in (9) using the annual sunspot data

For choices of ε up to two standard deviations, the test statistic is significant and the null hypothesis of i.i.d. is strongly rejected. This result is a warning that the apparent lack of structure obtained from a visual impression can be misleading. Interestingly, at 2.5 standard deviations, the null hypothesis cannot be rejected. It may be that this choice of ε is simply too large for these data.

6.2 U.S. equity returns

Consider the example given in Hurn et al. (2020), in which real U.S. equity returns are to be forecast using an autoregressive [AR(1)] model based on the assumption of normality. The model is

$r_t={\varphi }_0+{\varphi }_1{r}_{t-1}+v_t,\hspace{1em}{v}_t~N\left(0,{\sigma }_v^2\right)$

If the observed values r_t are indeed generated correctly according to this simple model, then the transformed quantity

$u_t=\Phi \left(\frac{v_t}{{\sigma }_v}\right),\hspace{1em}t=1,2,\dots ,T$

takes values in the unit interval [0, 1] because of the fundamental property of Φ(·), which is the cumulative distribution function (CDF) of the standard normal distribution. This transformation is known as the probability integral transform (Diebold, Gunther, and Tay 1998). Furthermore, Rosenblatt (1952) demonstrates that if the model is correctly specified, the probability integral transform, u_t , will be independent and uniformly distributed on the unit interval. The null hypothesis that the model is correctly specified can then be tested in terms of the BDS test applied to u_t .

Using monthly observations for the period January 1871 to September 2016 on real equity returns, we find estimation of the model gives the following results. ⁵

The probability integral transform is given by

$u_t=\Phi \left(\frac{{\widehat{v}}_t}{{\widehat{\sigma }}_v}\right)$

in which ${\widehat{\sigma }}_v$ is the standard error of the regression. The BDS for independence applied to u_t with embedding dimension m = 4 gives

The null hypothesis is rejected, and the conclusion is that the AR(1) model of equity returns is misspecified because u_t is not i.i.d. A histogram of the transformed time series, u_t , given in figure 3 suggests that the distribution of u_t is also not uniform. The interior peak of the distribution of u_t and also the peak at zero suggest that equity returns, r_t , are not consistent with the specification of an AR(1) model with normally distributed errors.

Figure 3.

Probability integral transform applied to the forecast errors of the AR(1) model of United States equity returns, January 1871 to September 2016

An interesting feature of the results is that the BDS statistic is much smaller for the choice of ε based on the fraction of pairs rather than on the standard deviations. This may be suggestive that the fraction of pairs method for choosing ε is more robust to the distribution of the underlying series.

7 Conclusion

In this article, we introduced the bds command, which computes the BDS test of the null hypothesis that the time series to which it is applied consists of i.i.d. observations. The algorithm is adapted from LeBaron’s C codes and produces the same output as these codes. In addition to providing the command, we also provided an elegant simplification of the computation of the variance of the statistic.

The fact that the procedure does not use compiled code but is written in Mata to run on all machines on which Stata is loaded means that the routine is computationally less efficient when the sample size is large. Consequently, this routine is not suitable for large-scale simulation studies. Two empirical examples demonstrated how the routine is implemented in practice.

8 Programs and supplemental materials

Supplemental Material, sj-zip-1-stj-10.1177_1536867X211025796 - The BDS test of independence

Supplemental Material, sj-zip-1-stj-10.1177_1536867X211025796 for The BDS test of independence by Christopher F. Baum, Stan Hurn and Kenneth Lindsay in The Stata Journal