2. LOCAL-TO-UNITY ASYMPTOTICS FOR ROBUST TESTS

The observations {y_t}_t=1^T come from the data generating process

where {ε_t}_t=1^T is an independent and identically distributed (i.i.d.) sequence with mean 0 and variance σ_ε². The term Γ(L) is the lag polynomial 1 − Γ₁ L − ··· − Γ_p L^p, and we assume that the roots of Γ(z) = 0 all lie outside the unit circle. For convenience we assume that u₀ = 0. We wish to test the unit root hypothesis γ = 0 versus the alternatives γ < 0.

The tests are based on the augmented Dickey and Fuller (1979) representation of the model

where ΔY = (y_p+2 − y_p+1,…,y_T − y_T−1)′, X = (y_p+1,…,y_T−1)′, and Z is the design matrix with row t − p − 1 equal to (1,t,Δy_t−1,…,Δy_t−p). The parameter α is γΓ(1), β is a (p + 2) × 1 parameter vector, and ε = (ε_p+2,…,ε_T)′. For the intercept only model (e.g., μ₂ = 0) the design matrix Z has rows (1,Δy_t−1,…,Δy_t−p).

M-tests are based on M-estimators, which may be characterized as solutions to minimization problems. It will be convenient to adopt the notation that, for any function g(·) mapping

and any n-dimensional column vector x with components x_i, g(x) is the n-dimensional column vector with components g(x_i). The M-estimators

minimize the objective function

where 1 is a vector of ones and g is a function chosen by the researcher. When g(x) = x², estimation is by least squares. Table 1 lists some alternative choices for g that “robustify” the estimators against thick-tailed error distributions. The unit root test rejects the null for small values of

.

Some robust unit root tests

Hasan and Koenker (1997) propose an alternative class of unit root tests based on ranks. The null is rejected for small values of [X′MX]^−1/2X′Mb, where M is the projection matrix I − Z(Z′Z)⁻¹Z′ and b is the vector of regression rank scores

, where

solves the linear programming problem

and φ is a function chosen by the researcher. This construction for b is a way to handle the nuisance parameters β and allows the researcher to consider general functions of the ranked residuals: when Z is a vector of ones and φ(τ) = τ, each element of b is the rank of the corresponding element of ΔY. Table 1 lists some common choices for φ.

In large samples the critical values and power functions of the various tests depend on several nuisance parameters. For the M-tests the nuisance parameters are determined by a scorelike function ψ(x), equal to g′(x) if g is differentiable at x and zero otherwise (we allow g to have finitely many points of nondifferentiability to handle quantile estimators such as least absolute deviations (LAD)). The nuisance parameters are

where

, f is the marginal density of ε₁, and η solves Eψ(ε₁ − η) = 0. When ψ is everywhere differentiable ω takes the more familiar form Eψ′(x − η). The rank tests depend on the same nuisance parameters, with ψ(x − η) replaced by φ(F(x)), where F(x) = Pr[ε₁ ≤ x] . Table 1 lists the g, ψ, and φ functions for some common robust tests, along with the associated formulas for ω.

The tests come in two versions, an e-test and a t-test. The e and t M-tests reject the null for small values of the statistics

where the variables with the [caret ] are consistent estimates of the corresponding nuisance parameters. The e and t rank tests reject for small

Hasan and Koenker (1997) propose only the t version of the test, but for completeness we examine the e version also. Section 4 provides estimates of the nuisance parameters.

Hasan and Koenker (1997) stress that the nuisance parameter λ does not need to be estimated to normalize a rank test. This may be an advantage of the rank tests over the M-tests. In the Monte Carlos in Section 5, the rank tests have more accurate size than the M-tests. Thompson (2004b) argues that this occurs because of estimation error associated with

. Although in large samples the estimation error will have an asymptotically negligible effect on inference, in small samples it could affect test size. In some cases ω is especially poorly estimated—for example, for LAD ω = 2f (η), and the estimate of a density at a point converges at less than a root-T rate.

Because I am interested in inference when γ is close to one, I adopt the local-to-zero reparameterization γ = c/T so that the parameter space is a shrinking neighborhood of zero as the sample size grows. Following Chan and Wei (1987) and Phillips (1987), I take c fixed when making limiting arguments, obtaining asymptotic representations as a function of the local alternative c.

Hasan and Koenker (1997) derive local-to-zero asymptotic representations for the rank tests. In Theorem 1 I provide local-to-zero representations for the M-tests. The theorem applies to objective functions g that are not twice continuously differentiable, such as LAD, Huber's M function, and the quantile regression estimator of Koenker and Bassett (1978). We assume that g is convex, which rules out Student's t objective function. Although the convexity assumption simplifies the proof, alternative methods can be used to show that the theorem holds for Student's t function.

The limiting forms of the test statistics are expressed as functionals of Brownian motion. Define W(·) as standard Brownian motion and define W_c(·) as the Ornstein–Uhlenbeck process

. Define

where D_c(r) = W_c(r) −

in the intercept only model and

in the model with a linear time trend. Under the local-to-null asymptotics, the e-tests M_e and R_e both converge to the limiting distribution represented by the random variable

, and the t-tests M_t and R_t converge in distribution to the random variable

. The limiting forms are

where S is an independent standard normal variable. These results are summarized in the following theorem.

THEOREM 1. Under Assumptions 1 and 2 in the Appendix, M_e and R_e converge weakly to

, and M_t and R_t converge weakly to

.

There are a number of well-known asymptotic results concerning robust inference for a location model (see Rothenberg and Stock, 1997, Appendix). Asymptotic power depends on a single parameter that shifts the limiting distribution of the test statistic. Power is increasing in the parameter, which is maximized under correct specification. The movement of the parameter is asymmetric under incorrect specification, in the sense that with normal data power declines only slightly when using robust methods instead of least squares, but with thick-tailed distributions robust tests lead to significant increases in power.

For the unit root testing problem, large sample power depends on the pair of nuisance parameters ρ and λ. Because we reject for small values of the test statistics, in large samples the power of an e-test is given by

, where q_e(ρ,a) satisfies

and a is the size of the test. Limiting power for the t-test is

, where

. Power obviously increases with λ, which shifts the distribution to the left when c is negative. The term ρ also affects power by changing the shape of the distribution. The results in Section 3.1 suggest that λ dominates the power function relative to ρ.

In general ρ and λ depend on the test and error density. An exception to this is the least squares–based Dickey–Fuller test, where ψ(x) = x and ρ = λ = 1 no matter what the error density. For other tests ρ will differ from 1 and the limiting null representations become a linear combination of a “Dickey–Fuller” term (B_c⁻¹A_c for the e-test and B_c^−1/2A_c for the t-test) and another random variable.

In Table 2 we have calculated ρ and λ by numerical integration for various tests and error densities.¹

Numerical integration was carried out using the Matlab function quad, which implements a quadrature algorithm. All integrals were solved between −40 and 40, so, for example, ω was calculated as

for the Wilcoxon scores tests. Kurtosis is defined to be E [(ε − Eε)⁴]/Var(ε)². We calculate the quantity by numerical integration.

Values of ρ and λ

To get a sense of how the parameters change as we move from normal errors to thicker tailed densities, we give results for error densities of the form f (ε) = A(τ)exp(−|ε|^τ) for τ = 2.0, 1.8, 1.2, and 1.0 and for the thick-tailed Student's t distribution with three degrees of freedom (denoted t3).²

The values of ρ fall as the errors become thicker tailed. For normal errors the parameter is close to 1, except for the LAD and sign rank tests. Under normality the null distributions of the other tests are similar to the null for the least squares–based tests. When the errors have thick tails the null distributions are different from the least squares distributions and from each other.

The term λ increases as the errors become thicker tailed, suggesting that thicker tailed errors lead to increased power. Under normal errors the parameters for the LAD and sign-median rank tests are far below 1, so we expect these tests to have poor power relative to the least squares–based tests when errors are normal. For the rest of the tests λ is close to 1 with normal errors, implying that these tests have power almost as good as least squares when errors are normal and better than least squares when errors are thick tailed. Table 4 confirms this analysis; there the robust tests essentially dominate least squares, even for distributions close to Gaussian.

For the test and errors densities considered in Table 2, λ is maximized by choosing −ψ(x − η) or −φ(F(x)) equal to ∂ log f (x)/dx. For normal errors, λ is largest with the least squares and normal ranks tests. The same result holds for the LAD and sign ranks tests used with double exponential errors and for the t3 test with t3 errors. This property holds for any test and error distribution, which is not surprising because it corresponds to testing under correct specifying of the error distribution. To see that it holds for the rank tests, let f (x) equal

, where A makes the density integrate to one. An argument based on integration by parts can be used to show that ω = σ_ψ². For any

with associated

, we have

Therefore

. The result can be extended to the M-tests by replacing φ(F(x)) with ψ(x − η).

It can be shown, using integration by parts, that ρ = λ⁻¹ when −ψ(x − η) or −φ(F(x)) equals ∂ log f (x)/dx. This may be confirmed in the tables for the least squares and normal ranks tests applied to Gaussian errors. It may also be confirmed for the LAD and sign ranks tests used with double exponential errors and for the t3 test with t3 errors. Thus, under correct specification of the errors only one nuisance parameter appears in the asymptotic representation of the test statistics.

3. OBTAINING CRITICAL VALUES

It is well known that the asymptotic null distributions of the usual Dickey–Fuller tests do not depend on any nuisance parameters. Critical values can be calculated once and tabulated. For robust, non-least squares–based tests, the critical values cannot be tabulated so simply. The asymptotically correct critical values q_e(ρ,a) and q_t(ρ,a) are the ath quantiles of

. For a given robust test ρ must be estimated, so the shape of the asymptotic null distribution is not estimable until we see the data.

There are many ways to accurately approximate the quantiles of

. Shin and So (1999) simulate draws of

for each ρ and use the empirical quantiles as critical values. Seo (1999) and Lucas (1998) also describe simulation-based procedures. Rothenberg and Stock (1997) and Thompson (2004b) use Cornish–Fisher expansions. Hansen (1997) and Juhl and Xiao (2002) simulate to calculate quantiles for similar test statistics.

A computationally convenient approach is to approximate q_e(ρ,a) and q_t(ρ,a) by a third-order polynomial in 1 − ρ. For the

statistic the formula is

where x_ρ = (1, 1 − ρ, (1 − ρ)², (1 − ρ)³)′ and γ_a is a four-dimensional coefficient vector. The formula is defined similarly for the

-statistic. To compute γ_a I simulated N = 5 million draws of

for ρ = ρ₀,…,ρ₁₀₀, where ρ_j = j/100.

as the (aN)th draws of the simulated statistics. For the

test I compute γ_a by regressing

on x_ρj. The minimization problem is

The same procedure is used for

. Table 3 provides the estimated coefficients for a equal to 0.01, 0.025, 0.05, and 0.10.

Estimated coefficients of polynomial approximation to the ath quantiles of the asymptotic null distributions

Local alternative where asymptotic power is approximately 50%

In practice I calculate critical values using a consistent estimate

for ρ. So the e-tests reject the null at level a when M_e or R_e are less than

. Similarly the t-tests reject when M_t or R_t are less than

. The estimation of

and the other nuisance parameters is discussed in Section 4.

Figure 1 plots the simulated quantiles along with the polynomial approximations. The quantiles are nearly linear in ρ, and it is difficult to see a difference between the quantiles and the approximating polynomial. I assess the quality of the approximation in terms of the worst large sample size distortion:

Simulated 1%, 2.5%, 5%, and 10% quantiles (solid lines) plotted along with the polynomial approximations from Table 3 (dotted lines). The approximations are so close that it is difficult to distinguish the solid from the dotted lines. The horizontal axes give values of ρ.

The expression

is an estimate of

. I calculate

by simulating 5 million draws of

, independent of the draws used to estimate γ_α, and calculate

as the percentage of draws below z. The worst size distortion was calculated for the models with and without a trend and for both the

- and

-statistics. The maximum over the four cases was 0.0004942. Thus the polynomial approximations are highly accurate both in terms of approximating the quantiles and minimizing asymptotic size distortions.

3.1. Modifying the Tests to Eliminate the Nuisance Parameter

Herce (1996) and Hasan and Koenker (1997) suggest modifying the test statistics to obtain null distributions that are free of ρ. Recall that

are expressed as linear combinations of the “Dickey–Fuller” terms A₀ /B₀ and

and another random variable. Because the least squares test statistics converge to the Dickey–Fuller terms, we can write down a linear combination of the robust test statistics and of the least squares statistics so that the limiting null distribution does not depend on ρ. Let e_ls and t_ls denote the m_e- and m_t-statistics when estimation is by least squares. The modified tests are based on the relationships

where

is a consistent estimator for ρ. A test that rejects the null when the modified statistics are below the 100a% quantile of S/B₀^1/2 or S will have asymptotic size equal to a. The rank tests may be similarly modified.

Unfortunately, this modification seriously affects the power of the tests. Let Φ and φ denote the distribution function and density function of a standard normal variable. When the errors are i.i.d. normal, the rank test formula for λ is

The same result holds for the M-tests. Thus, under normal errors λ = ρ no matter what the test, so the modified t-test has local power everywhere equal to its size. This may be confirmed in Table 2, where under normality λ equals ρ for every test considered. In Table 4 we also show that the modified e-test also has very little power with normal errors. Lucas (1995b) reports simulation results for M-tests that are consistent with these results. Herce (1996) obtains a stronger result for the LAD estimator—under normal errors the modified LAD t-test has power equal to its size. Although our results apply only to the local power function, in Section 5 we provide Monte Carlo evidence that the actual power function exhibits similar behavior for empirically relevant sample sizes. Robust methods are meant to perform well for plausible misspecifications of the model. By this measure the traditional Dickey–Fuller methods are robust relative to the modified statistics.

Table 4 provides some evidence to compare the modified tests with the unmodified tests that use the polynomial approximation to the critical value. The table lists the local alternative where asymptotic power is approximately ½ in the model that includes a linear time trend (e.g., no restrictions on μ₂). King (1988) and Elliott, Rothenberg, and Stock (1996) make the case for preferring tests that achieve 50% power at alternatives closer to zero. By definition, tests that achieve 50% power at lower values of c have better asymptotic relative efficiency at 50% power (see Stock, 1994, p. 2775).

The modified tests have virtually no power for distributions close to the Gaussian, but sometimes they beat the unmodified tests for the thick-tailed double exponential and t3 errors. Asymptotic power was calculated for local alternatives as small as −100, so a value of −100 in the table means that the test did not achieve 50% power for any alternative greater than −100. Although the modified e-test has power greater than size for Gaussian errors, the modified e-tests considered here do not achieve 50% power for any alternative greater than −100. Both the modified and unmodified e-tests beat the corresponding t-tests for error distributions close to normality, and the t-tests are more powerful for thicker tailed error distributions.

The unmodified robust tests essentially dominate least squares. The least squares–based e- and t-tests of Dickey and Fuller (1979) achieve 50% power at c = −15.75 and −17.5; the unmodified robust tests do slightly worse for normal errors and much better for nonnormal errors. For each robust test the presence of thick-tailed errors leads to greater power relative to normal innovations. There is a larger decline in power under normality for the LAD and sign ranks tests, but even in these two cases the decrease is small relative to the gain in power for thicker tailed distributions.

The normal ranks test is at least as efficient as the least squares test for every error distribution in the table. In the classical location model we can make the stronger statement that the normal ranks test is asymptotically at least as powerful as the least squares t-test no matter what the error distribution. This result is shown by Chernoff and Savage (1958) (also see the discussion in Hettmansperger, 1984, Sec. 2.9). The result may hold for the unit root problem, but I was not able to prove it. At normal errors both the Dickey–Fuller and the normal ranks tests have λ = ρ = 1, so their power functions are identical. I am able to show that for the normal ranks test, λ ≥ 1 and ρ ≥ 1/λ for any error distribution. Power obviously increases in λ, but simulations show that power is not monotonic in ρ. However, I conducted simulations⁴

I carried out the following simulation. Set ρ_i = i/50, λ_j = 50/j, and c_k = −k. For each (ρ, λ, c) = (ρ_i, λ_j, c_k) triple over i = 1,…,50, j = 1,…,50, k = 1,…,50, I simulated 100,000 Monte Carlo draws of

. For each c I computed power for a 5% level test using the polynomial approximations to the critical values given in Table 3.

suggesting that in the model with a time trend, large sample power of the e- and t-tests for λ ≥ 1 and ρ ≥ 1/λ is always at least as high as when λ = ρ = 1. In the model with no time trend the result holds only for the e-test. In that model the t-test has about 1% less power than the Dickey–Fuller test in the range 0.90 ≤ ρ ≤ 0.92 and 1.06 ≤ λ ≤ 1.11. Although I am unable to find an error distribution with these nuisance parameters, one may exist.

Chernoff and Savage (1958) establish that λ ≥ 1 for the normal ranks test; see for example their Theorem 3 or Theorem 2.9.2 of Hettmansperger (1984). To see that ρ ≥ 1/λ, note that σ_ψ² = 1 and assume without loss of generality that σ_ε² = 1. Thus

where f is the density of the errors and φ is a standard normal density. The second equality follows from the change of variables x = Φ⁻¹(F(y)), and the third equality follows from the fact that ∂φ(x)/∂x = −xφ(x). Integration by parts and a change of variables leads to the expressions

Let W = f (ε₁)/φ(Φ⁻¹(F(ε₁))). By Jensen's inequality,

3.2. Lucas's (1995b) Method for Obtaining Critical Values

Lucas (1995b) obtains critical values by computing the distribution of the test statistics from realizations of simulated data. Lucas simulates the data by generating a time series from the data generating process with i.i.d. standard normal innovations v_t and with d_t and γ set to zero. Because the parameter ρ depends on both the choice of estimation method and the error density, this approach will give asymptotically incorrect critical values for the robust tests unless the errors are in fact normal. For the tests based on the t3 estimator, ρ = 0.98 for normal errors and the 5% quantiles of

are −21.8 and −3.42 when the model has a time trend. For t3 errors ρ = 0.71 and the quantiles are −17.5 and −3.01. A simulation shows that when the errors are distributed t3, using the quantiles based on normal errors instead of the correct values results in e- and t-tests with asymptotic sizes of 1.59% and 1.68%, respectively.

⁵

There are 50,000 Monte Carlo replications.

4. CALCULATING THE TESTS

To calculate the test statistics we need consistent estimates for the various nuisance parameters. Likewise, the polynomial approximation to the critical values,

, is evaluated at a consistent estimate

for ρ.

Following Hasan and Koenker (1997), the nuisance parameters are estimated under the null. For a given M-test with associated functions g and ψ, let

denote the vector of residuals

, where

minimizes 1′g(ΔY − Zβ). Theorem 1 implies that under the null,

will converge to the recentered errors ε − η. So we estimate the two variances σ_ε² = Var(ε₁) and σ_ψ² = Var(ψ(ε₁ − η)) with

where

are the sample averages of

. Similarly, the correlation ρ = Corr[ε₁,ψ(ε₁ − η)] is estimated by

For Huber's M-test and the t3-test ω is estimated by

, where ∇ψ(x) equals ψ′(x) if ψ is differentiable at x and zero otherwise. Following Herce (1996), for the LAD test I estimate ω = 2f (η) using the kernel estimator

, where

. The bandwidth h was chosen based on a recommendation in Silverman (1986, Ch. 2). The term

will converge in probability to f (η) because

is converging to ε − η. Following Lucas (1995b), the scale parameter k appearing in Huber's M function is

, where

are the residuals from a preliminary LAD regression.

For the rank tests the nuisance parameters are calculated using the same formulas. I follow Hasan and Koenker (1997) and set the residuals

equal to

, where

minimizes the LAD objective function 1′|ΔY − Zβ|. For the rank tests ψ(x) is replaced by

, where

is an estimate of F. Following Hasan and Koenker (1997), I take

to be the empirical distribution function of

, equal to ({#residuals ≤ x} − 0.5)/(T − p). The term σ_ψ² is known in closed form for the rank tests and does not need to be estimated; σ_ψ² is

for the Wilcoxon ranks test and 1 for the normal and sign ranks tests.

Computing the rank tests requires calculating the vector

, where

solves the linear programming problem in (2). Hasan and Koenker (1997) use an algorithm that calculates

for all s. I instead solve the linear programming problem in expression (2) at the points {τ_i}_i=1⁵⁰ (with τ_i = (i − 1)/49) and use numerical integration. For the Wilcoxon ranks test,

. For the normal ranks test, b is

and for the sign ranks tests it is

. I use the discretized version of the

process because canned algorithms for solving linear programming problems are widely available—for example, I used the Matlab function linprog.

5. MONTE CARLO EVIDENCE

We carried out a Monte Carlo study to evaluate the small sample performance of the various tests. The calculation of the various tests was described in Section 4.

Table 5 lists empirical rejection rates under the null (γ = 0) and for a distant alternative (γ = −0.5). The errors ε_t are i.i.d. draws from either the N(0,1) or Student's t3 distribution. We make the unrealistic assumption that the researcher knows that the errors are i.i.d. and sets the lag length p equal to zero. Results are given for both the intercept only model and for the model with a time trend. Small sample size looks good for all the tests except the modified e. The LAD test exhibits some size distortions. Thompson (2004b) argues that this results from a poor estimate of ω, because

does not converge at a root-T rate.

Empirical rejection rates, i.i.d. errors

Previously we have shown that the modified tests have poor power against the local alternative c. The Monte Carlo results suggest that power is also poor against the distant alternative γ = −0.5. Test power is very close to one for most of the unmodified tests, whereas the modified tests suffer from large power losses. As the asymptotics predict, the modified t-tests have power equal to size for normal errors. The asymptotics also predict that the modified tests will beat the unmodified ones for the thick-tailed Student's t3 distribution, and this does not seem to be the case. One explanation for this is that the local-to-null asymptotic approximations become less accurate for distant alternatives such as γ = −0.5. Thus the Monte Carlo results make a stronger case against the modified tests than the asymptotics.

The calculations in Table 5 assume that the errors are independent. To investigate the impact of serially correlated errors, Table 6 provides evidence about the small sample behavior of the unmodified tests when the errors exhibit serial dependence of unknown form. The data generating processes are

Empirical rejection rates with serially correlated errors

The fundamental innovations ζ_t have initial condition ζ₀ = 0 and were generated from either the N(0,1) or the t3 distribution. The table gives results for the intercept only model and also for the model that includes a time trend.

Bringing the tests to serially correlated data requires an algorithm for selecting the lag length p. The term p is chosen by the Bayesian information criterion of Schwarz (1978) constrained so 2 ≤ p ≤ 6. Because the tests are designed for thick-tailed data, we calculate Schwarz's criterion from a double-exponential likelihood with constant variance.

Size is generally accurate for AR errors, less accurate for MA errors, and worse for GARCH errors. The t-tests exhibit more accurate size than the e-tests. For the rank t-tests, rejection rates under the null range from 0.04 to 0.11. The rank tests in general exhibit better size properties than the M-tests. Thompson (2004b) argues that this occurs because of estimation error associated with

. Computing the rank tests does not require an estimate of

. Although in large samples the estimation error will have an asymptotically negligible effect on inference, in small samples it could affect test size.

As the asymptotic analysis predicted, the robust tests essentially dominate least squares at the alternative γ = −0.1. Except for the LAD and sign ranks tests, the robust tests are about as good as least squares at normal errors and are much more powerful at t3 errors. The asymptotics also suggested that the e-tests would be more powerful than the t-tests for normal errors but less powerful for t3 errors, but in the Monte Carlo the e-tests are at least as powerful as the t-tests for both distributions. Taking the size and power results together, the unmodified Wilcoxon ranks and normal ranks t-tests are nice alternatives to least squares: both rank tests have accurate size, and both do at least as well as least squares over a range of error distributions.