Unit Root CADF Testing with R

This paper describes CADFtest, an R package for testing for the presence of a unit root in a time series using the covariate-augmented Dickey-Fuller (CADF) test proposed in Hansen (1995). The procedures presented here are user friendly, allow fully automatic model specification, and allow computation of the asymptotic p values of the test.


Introduction and statistical background
Testing for a unit root is a frequent problem in macroeconomic modeling and forecasting.Although many tests have been developed so far, the practitioner's workhorse in this field is still probably the Augmented Dickey-Fuller (ADF) test (Dickey andFuller 1979, 1981;Said and Dickey 1984), which is known to have good size but low power under many conditions. 1 However, the literature on unit root testing has largely proceeded in a univariate framework, with notable exceptions being represented by the panel unit root tests (see Choi 2006b, for a recent survey).In fact, reality is hardly univariate at all and, although univariate representations of multivariate time series exist (see e.g., Zellner and Palm 1974), nevertheless reckoning explicitly the multivariate nature of most economic time series can in principle lead to better testing procedures.In a seminal paper Hansen (1995b) suggests that, when testing for a unit root, a viable way to exploit the information embodied in related series and increase power is to consider stationary covariates in an otherwise conventional ADF framework.Unless the variable of interest is independent of the stationary covariates considered in the analysis, by the Neyman-Pearson lemma the most powerful test makes use of the information embodied in the stationary covariates themselves (Hansen 1995b(Hansen , p. 1152)).As a consequence, not considering the multivariate dimension of the problem leads to a loss in the power of the test.Using covariates also allows to some extent to couple unit root testing and economic theory, because economic theory can be used as a guideline to chose the appropriate covariates to be included in the analysis (see e.g., Amara and Papell 2006;Elliott and Pesavento 2006), although other approaches can be used as well (Lee and Tsong 2009).
Formally, Hansen (1995b) considers a univariate time series y t composed of a deterministic 1 However, large negative MA components are well known to have adverse effects on the size of the ADF test (see e.g., Schwert 1989).An interesting summary of many Monte Carlo results on different unit root tests can be found in Stock (1994).Haldrup and Jansson (2006) is a recent survey of the methods proposed to improve the size and the power of unit roots tests.and a stochastic component such that a(L)∆s t = δs t−1 + v t (2) where d t is the deterministic term (usually a constant or a constant and a linear trend), a(L) := (1 − a 1 L − a 2 L 2 − . . .− a p L p ) is a polynomial in the lag operator L, ∆x t is a vector of stationary covariates, µ ∆x := E(∆x), b(L) := (b q 2 L −q 2 + . . .+ b q 1 L q 1 ) is a polynomial where both leads and lags are allowed for.Furthermore, the long-run (zero-frequency) covariance matrix Ω := ∞ k=−∞ E t t−k , with t := (v t , e t ) can be defined, from which the long-run squared correlation between v t and e t , ρ 2 , can be derived.When ∆x t has no explicative power on the long-run movement of v t , then ρ 2 ≈ 1.On the contrary, when ∆x t explains nearly all the zero-frequency variability of v t , then ρ 2 ≈ 0. The case ρ 2 = 0 is ruled out for technical reasons (see Hansen 1995bHansen , p. 1151)).It should be noted that this restriction excludes the possibility that the variable y t be cointegrated with the cumulated stationary covariate(s).
As with the ADF test, Hansen's CADF test is based on different models, according to the different deterministic kernels that the investigator may wish to consider.For example, the model with constant and linear trend is Similarly to the conventional ADF test, the CADF test is based on the t-statistic for δ, t(δ), with the null hypothesis being that a unit root is present, i.e.H 0 : δ = 0, against the one-sided alternative H 1 : δ < 0. Hansen (1995b) refers to the test statistic as the CADF(p, q 1 , q 2 ) statistic.Hansen (1995bHansen ( , p. 1154) ) proves that under the unit-root null, if conventional weak dependence and moment restrictions hold, t(δ) in the model without deterministic terms is such that where ⇒ denotes weak convergence, W is a standard Wiener process, and N (0, 1) is a standard normal independent of W .It is interesting to note that (5) is the distribution of a weighted sum of a Dickey-Fuller and a standard normal random variable.If a model with constant or a model with constant and trend are considered, the standard Wiener process W in (5) has to be replaced by a demeaned or a detrended Wiener process, respectively.Note that the asymptotic distribution of the test statistic depends on the nuisance parameter ρ 2 but, provided ρ 2 is given, it can be simulated using standard techniques (see e.g., Hatanaka 1996).Hansen (1995bHansen ( , p. 1155) ) provides the asymptotic critical values of the test, while here we offer a practical way to compute its p values.
Hansen's CADF test is firmly rooted in the ADF tradition and for this reason it can be more familiar and attractive to practitioners than other tests, although Elliott and Jansson (2003) show that Hansen's CADF test is not the point optimal test in general.Feasible point optimal tests in the presence of deterministic components are developed in Elliott and Jansson (2003).
A quite different approach is used in this case, the feasible tests being based on VAR models.
However, Monte Carlo simulations reported in Elliott and Jansson (2003) suggest that power gains with respect to Hansen's CADF test can be obtained at the cost of slightly worse size performances.
In this paper we present the R (R Development Core Team 2009) package CADFtest that allows users to perform Hansen's CADF unit root test easily.The main function CADFtest() returns a CADFtest class object that not only contains the test statistic, but also its asymptotic p value and many other useful details.In fact, the class CADFtest inherits from the class htest,2 so that no special print() method is needed.However, dedicated summary() and plot() methods have been developed in order to allow the user to analyze the test results more in detail.A specialized update() method is also available that ease testing using different options.
The remainder of the paper is structured as follows: Section 2 discusses the way Hansen's CADF test has been implemented in the function CADFtest(), and some applications are illustrated in Section 3. In Section 4, the method to compute the asymptotic p values is illustrated in detail along with the use of the function CADFpvalues().Section 5 offers some comparisons with other existing R packages performing the ADF test.A summary is offered in Section 6.

Implementation and use of the function CADFtest()
In principle, carrying out Hansen's CADF test is no more complicated than carrying out an ordinary ADF test.What makes things more complicated is the presence of the nuisance parameter ρ 2 in the asymptotic distribution (5).In fact, a consistent estimate of ρ 2 has to be derived in order to choose the correct asymptotic critical value and/or to compute the correct asymptotic p value of the test.The problem is solved into two steps.First, êt and vt are derived; then, their long-run covariance matrix Ω is estimated using a HAC covariance estimator (see e.g., Andrews 1991;Zeileis 2004Zeileis , 2006;;Kleiber and Zeileis 2008).
Once the kind of model (no constant, with constant, with constant and trend) has been chosen, using CADFtest() the investigator can either set the polynomial orders p, q 2 and q 1 to fixed values, or can decide the maximum value for each and let the procedure to select and estimate the model according to different information criteria.
In order to estimate ρ 2 it is necessary to estimate e t and v t first.For example, if the model with constant and trend (4) is used, then e t and v t are estimated as where " " denote parameters estimated by ordinary least squares and ∆x is the sample average of ∆x t .Once êt and vt have been computed, a kernel-based HAC covariance estimator (Andrews 1991) is used to estimate Ω and hence ρ 2 .In order to estimate ρ 2 in a rather flexible way, in CADFtest() the kernHAC() function included in the sandwich package (Zeileis 2004(Zeileis , 2006) ) is used.This allows the investigator to chose the kernel to be applied, the bandwidth, and if and how prewhitening should be performed.Differently from Hansen (1995b) where a Parzen kernel without prewhitening is used, the default choice in CADFtest() is a quadratic spectral kernel with VAR(1) prewhitening.The bandwidth is adaptively chosen using the method proposed in Andrews (1991), but the user is free to change any of these default choices.
The usage of the function is extremely simple: CADFtest(model, X = NULL, type = c("trend", "drift", "none"), data = list(), max.lag.y= 1, min.lag.X = 0, max.lag.X = 0, dname = NULL, criterion = c("none", "BIC", "AIC", "HQC", "MAIC"), ...) The minimal required input is CADFtest(y), where y can be either a vector or a time series.However, if no stationary covariate is specified, an ordinary ADF test is performed.In fact, the ordinary ADF test can be carried out with R using other existing packages such as fUnitRoots (Wuertz 2009), tseries (Trapletti 2009) and urca (Pfaff 2008).In this respect there would be no need to add one further package.However, given that the ADF test can be seen as a particular case of the more general CADF test, it seems logical to leave the user the opportunity to carry out both tests in the same framework, using the same conventions and allowing for the computation of the test p values.Furthermore, as will be shown in section 3, the interface to CADFtest() is very flexible and intuitive, and the results are easy to read: this can make carrying out conventional ADF tests using CADFtest() very appealing.As far as the computation of the p values of the ADF test is concerned, CADFtest() exploits the facilities offered by punitroots() implemented in the package urca (Pfaff 2008) that uses the method proposed in MacKinnon (1994MacKinnon ( , 1996)).In principle, it would have been possible (and easy) to use the function CADFpvalues() described in detail in section 4, but given that MacKinnon (1996) describes a method to compute approximate finite sample, rather than asymptotic, p values, it seems fair to refer directly to a function that implements this procedure.However, note that MacKinnon (1996) derives the finite sample p values for Gaussian data: in non-Gaussian settings the finite sample p values are not necessarily more accurate than the asymptotic ones.
All the arguments, with the exception of min.lag.X and max.lag.X that are relevant only when a CADF test is carried out, work irrespective of the test being ADF or CADF.However, if a proper CADF test has to be performed, at least a stationary covariate must be passed to the procedure.The covariates are passed in a very simple way, using a formula (model) statement.For example, suppose we want to test the variable y using x 1 and x 2 as stationary covariates: if the other default options are accepted, then it is sufficient to specify CADFtest(y ~x1 + x2).Note that the formula that is passed as argument to the CADFtest() function is not the complete model to be used, but it just indicates which variable has to be tested for a unit root (y) and which are to be used as stationary covariates in the test (x1 and x2).A formula statement can be used also to specify an ordinary ADF test by typing CADFtest(y 1), where the term "1" does not imply that a model with constant will be used, but it simply means that no stationary covariate is passed to the procedure (the deterministic kernel is always defined by the argument type described below).
Other arguments are used to specify the deterministic kernel to be used in the model (type), the lead and lag orders (max.lag.y,min.lag.X, max.lag.X), and if the model has to be fixed or selected using a criterion such as "AIC" (Akaike 1973), "BIC" (Schwarz 1978), "HQC" (Hannan and Quinn 1979) or "MAIC" (Ng and Perron 2001).Indeed, given that a number of competing models with potentially many regressors have to be compared, information criteria offer a handy solution.Furthermore, Hall (1994) shows that when the data are generated by an ARIMA(p 0 , 1, 0) process, then the distribution of the ADF test statistic under the null converges asymptotically to the correct distribution when the number of lags in the empirical model is determined by using either the AIC, the BIC, or the HQC.On the other hand, Ng and Perron (2001) argue that standard information criteria should be modified to take into account the fact that the series are I(1) under the null and propose their modified AIC (MAIC) that should be more robust in the presence of negative moving-average errors.Ng and Perron's MAIC is computed by CADFtest() in the OLS-detrended version suggested by Perron and Qu (2007).However, notice that although the MAIC should in principle work well also in the CADF framework, its usefulness has been proved only in the simpler ADF context.
When no stationary covariate is passed to the procedure, then lag selection is obviously limited to the lags of ∆y t , but when a proper CADF test is performed, then model selection implies the joint determination of the lags of the differenced dependent variable and the leads and lags of the stationary covariates as well.If criterion = "none" (the default choice) is specified, no automatic model selection is performed and the lag orders are fixed to the values passed to the procedure.In particular, max.lag.ycorresponds to p, the lag order of a(L) in ( 4), and it is set to 1 by default: it can be set equal to 0 or to a positive integer.min.lag.X corresponds to q 2 , the maximum lead in b(L) in ( 4), and it is equal to 0 by default: if modified, it must be set equal to a negative integer (a negative lag is a lead).max.lag.X correspond to q 1 , the maximum lag in b(L) in ( 4), and the default choice is 0: if modified, it must be set equal to a positive integer.If criterion is different from "none", then all the models with lags polynomials up to the specified orders (of both the y and the covariates) are estimated and the final model to be used is selected on the basis of the chosen criterion.The deterministic components to be used in the model are specified using the conventions utilized in the R package urca (Pfaff 2008).The default value (type = "trend") implies that the model with constant and trend (4) is used.If type = "drift" or type = "none" is specified, then the model with constant or the model without deterministic components are utilized.data is the data set to be used and dname is the name of data: in general there is no need to change dname, given that it is automatically computed by the function itself, unless one wants to indicate for example that a specific data set has been used.Further arguments can be passed to the procedure to control the parameters to be used in the HAC covariance estimation.These further arguments can be passed using the conventions valid for the command kernHAC() (see the package sandwich : Zeileis 2004: Zeileis , 2006)).If Hansen's results have to be replicated, then kernel = "Parzen" and prewhite = FALSE have to be specified, otherwise a quadratic spectral kernel with VAR(1) prewhitening is used by default.
The function CADFtest() returns an object of class CADFtest containing the test statistic (statistic), the p value of the test (p.value), the lag structure of the selected model (max.lag.y,min.lag.X, max.lag.X), the value of the information criteria (AIC, BIC, HQC, MAIC), the estimated value of ρ 2 (parameter), and the full estimated model (est.model).Other returned information concern the nature of the test (either CADF or ADF) stored in method, the name of data used (data.name), the value of δ under the null (null.value), the description of the alternative (alternative) and the estimated value of δ (estimate).
A summary of the test can be obtained just by using a print() command.Given that the class CADFtest inherits from the class htest, the print() command produces the standard R output of the htest class.However, the summary() command is also allowed that returns a more detailed account of the test results.For greater flexibility, print() can be applied to a CADFtestsummary object (produced by summary()) to control further printing options.For example, the number of significant digits can be controlled by digits, while significance stars can be avoided by setting signif.stars= FALSE.

Some examples of application
We provide here some simple examples of application of the function CADFtest().Data are taken from the R package urca (Pfaff 2008) and refer to the extended Nelson and Plosser (1982) data set used in Schotman and Van Dijk (1991).These are the same data used in Hansen (1995b), so that we will be able to replicate some of the empirical applications proposed there.
First, we load the data and the required package CADFtest: all the following examples assume that both have been loaded.R> data("npext", package = "urca") R> library("CADFtest") A complete description of the data can be retrieved simply by typing ?npext in R.
We first replicate the analysis carried out in Hansen (1995bHansen ( , p. 1165) by testing for the presence of a unit root in the log per capita US real GNP using a standard ADF test with constant, trend and three lags: The p value of the test is stored in ADFt$p.valueand it is easily accessible: As already mentioned, the finite sample p value is computed using punitroots() implemented in package urca (Pfaff 2008).In principle, it would have been possible also to compute the asymptotic p value by using the function CADFpvalues to be described in detail in the next section by invoking R> CADFpvalues(ADFt$statistic, type = "trend", rho2 = 1) [1] 0.07589502 When a standard Dickey-Fuller test is performed, CADFtest() acts as an interface to existing commands.For example, in the case above equation ( 4) is estimated using the package dynlm (Zeileis 2009) and the test p value is computed using punitroots().
Even if all the results are readily accessible, a summary of the test can be obtained just by typing R> print(ADFt) ADF test data: npext$gnpperca ADF(3) = -3.2606,p-value = 0.08082 alternative hypothesis: true delta is less than 0 sample estimates: delta -0.2014652 The function correctly warns the user that a conventional ADF test has been performed and reports the main results along with the number of lags used in the test.
If we want to obtain a more detailed summary that includes the details of the estimated model, we can just type R> summary(ADFt) The model output uses the same conventions utilized in the package dynlm (Zeileis 2009): trnd is the deterministic linear trend, L(y, 1) stands for y t−1 and L(d(y), i) represents ∆y t−i .Note that the p value of the lagged dependent variable refers to the unit root null and is therefore consistent with the test p value.Note also that, differently from the conventional usage, the F statistic here pertains to the joint significance of the stationary regressors, so that under the null it has the standard F distribution (Sims, Stock, and Watson 1990).If a simple DF test is performed, then the F -statistic is not computed and a NA value is returned.If more control on the output summary is desired, then it is possible to store the summary results in an object (of class CADFtestsummary) and print it using print() with the desired options (for example, digits = 3, signif.stars= FALSE).
Further details about the test can be gathered from the estimated residuals and from residuals plots.The function residuals() can be used to extract the estimated residuals as in the following line:

R> res.ADFt <-residuals(ADFt)
The extracted residuals can then be used in any desired analysis.Fairly informative plots can also be easily produced by a plot() command as in

R> plot(ADFt)
Four plots are produced by default, as in Figure 1.In particular, the standardized residuals, the estimated residuals density along with the test for normality proposed in Jarque and Bera (1980), the estimated residuals autocorrelation function (ACF) and partial autocorrelation function (PACF) are plotted.However, any combination of these plots can be produced as well.For example, if the residuals density is not needed, then it is sufficient to specify R> plot(ADFt, plots = c(1, 3, 4)) to produce a visualization as in Figure 2.
In order to show other useful features of the CADFtest() command, we carry out now a few data transformations: Data are now interpreted as annual time series starting in 1860.The sample ends in 1988 (this is easy to verify by invoking the tsp() function).Given that unemploy is the log of the unemployment rate, while we need the unemployment rate, the series in levels used by Hansen (1995b) is computed.The time series in levels are stored in L, while D stores the first differences of the original variables, that will be used as stationary covariates in the CADF tests.S contains all the series over a common sample that starts in 1909, as in Hansen (1995b).
The ADF test could have been also performed by invoking R> ADFt <-CADFtest(L.gnpperca~1, data = S, max.lag.y= 3) Since no stationary covariate is explicitly indicated in the model, the test is performed as an ordinary ADF test (with constant and trend and three lags, as before).
Automatic lag selection can be achieved by using the criterion argument in the CADFtest() command.For example, in the following we fix the maximum lag order p of ∆y t−p to p = 4 and let the final model to be selected by the BIC, possibly highlighting that the data are from the extended Nelson and Plosser (1982)  or by updating the object ADFt that contains the results of a previous test: R> ADFt2 <-update(ADFt, change = list("max.lag.y= 4", "criterion = 'BIC'", + "dname = 'Extended Nelson-Plosser data'")) Of course, when automatic lag selection is enabled, all the models are estimated on the same common sample.
Let's now turn back again to our ADF(3) test performed by using CADFtest(): R> ADFt <-CADFtest(L.gnpperca~1, data = S, max.lag.y= 3) Now, suppose that we want to run Hansen's CADF test on the log-real GNP per capita by using the first difference of the unemployment rate as a stationary covariate.The test is carried out with constant and trend and allowing 3 lags for the (differences of the) dependent variable and 0 lags for the covariate.For the results to be fully consistent with Hansen (1995b, Table 8, column 2, p. 1166), kernel = "Parzen" and prewhite = FALSE have to be specified.Given the last ADF was carried out using 3 lags, we can perform the CADF test just by calling R> CADFt <-update(ADFt, change = list("+ D.unemrate", "kernel = 'Parzen'", + "prewhite = FALSE")) R> print(CADFt) CADF test data: L.gnpperca ~D.unemrate CADF(3,0,0) = -3.413,rho2 = 0.064, p-value = 0.001729 alternative hypothesis: true delta is less than 0 sample estimates: delta -0.08720302 Differently from Hansen (1995b), we not only verify that the test is significant at the asymptotic 1% level, but we can also give a precise assessment of the test p value.

p values computation and the function CADFpvalues()
The possibility of computing the p values of a test greatly increases the chances that the test is effectively used by practitioners.This is a fortiori true when the test procedure requires the use of non-standard tables available only in few specialized papers.There are even instances where computation of the p values is necessary for further investigations, as is the case for some panel unit root tests (see e.g., Maddala and Wu 1999;Choi 2001Choi , 2006a;;Costantini, Lupi, and Popp 2007).The R function CADFpvalues() presented here allows the computation of asymptotic p values of the CADF test proposed in Hansen (1995b).CADFpvalues() is used within the CADFtest() function to compute the p values of the test along with the other test results already discussed.However, CADFpvalues() can also be used separately from the main testing procedure.
The method used to compute the p values has been originally proposed in Costantini et al. (2007) and is based on a response surface approach similar to that proposed in MacKinnon (1994MacKinnon ( , 1996) ) for the p values of the ADF test (classical references on the estimation and use of response surfaces are, among others Hendry 1984;Ericsson 1986).Differently from what happens with reference to the Dickey-Fuller distribution which is free from nuisance parameters, the asymptotic distribution of the CADF test statistic depends on the nuisance parameter 0 < ρ 2 ≤ 1, so that the asymptotic distribution (5) has to be simulated over a grid of values for ρ 2 .When ρ 2 = 1 the distribution coincides with the ordinary Dickey-Fuller distribution.
In order to obtain fairly good approximations, here a grid of 40 values ρ 2 ∈ {0.025, 0.050, 0.0725, . . ., 1} is considered.For each of the three models (without deterministic components, with constant, with constant and linear trend), 100,000 replications3 have been used for each value of ρ 2 .The Wiener functionals have been simulated using a standard approach (see e.g., Hatanaka 1996, p. 67) with T = 5, 000 (for the "no constant", "constant" and "constant plus trend" case, standard, demeaned and detrended Wiener processes have been used, respectively).On the basis of the simulated values, for each value of ρ 2 1,005 asymptotic quantiles q ρ have been derived corresponding to the probabilities p = (0.00025, 0.00050, 0.00075, . .., 0.001, 0.002, . .., 0.998, 0.999, 0.99925, 0.99950, 0.99975).As a result, we obtained a 1, 005 × 40 matrix of estimated quantiles.Along the rows of the matrix it is possible to read how a given quantile varies with ρ 2 .Indeed, the estimated quantiles vary very smoothly with ρ 2 (see Costantini et al. 2007).
For each row of the quantile matrix the model is estimated and the β's are saved in a 1, 005 × 4 table.The tables of estimated coefficients for the "no constant", "constant" and "constant plus trend" case, respectively are used by the function CADFpvalues() in order to compute the asymptotic p values for any value of 0 < ρ 2 ≤ 1 for the relevant model.4 The way the computation of the p values proceeds in CADFpvalues() is essentially the following: 1.The relevant table of parameters is read, depending on the specific model used ("no constant", "constant" or "constant plus trend").
2. For any desired value ρ 2 0 of ρ 2 , the estimated parameters are used to compute for all the 1,005 probability values p the fitted quantiles q ρ 0 (p) as 3. The approach suggested in MacKinnon (1994MacKinnon ( , 1996) ) can now be used on q ρ 0 to derive the p value.First, given the value t(δ) of the test statistic, it is necessary to find the fitted quantile q ρ 0 that is closest to t(δ) and the corresponding probability p.

The regression
where Φ −1 (p) is the inverse of the cumulative standard normal distribution is estimated locally on an interval of p centered on p.In CADFpvalues() local interpolation takes place using 11 values centered on p.
5. The p value associated with the estimated test statistic t(δ) is finally obtained from Φ γ0 + γ1 t(δ) + γ2 t(δ) While computation is rather involved, from the user's viewpoint the usage of the function is extremely simple: CADFpvalues(t0, rho2 = 0.5, type = "trend") where t0 is the value of the test statistic t(δ), rho2 is the estimated value of ρ 2 , and type assumes the values "trend" (the default), "drift" or "none" as discussed above when a model with constant plus trend, with constant or without constant is considered.
If desired, CADFpvalues() can be used also to compute the asymptotic p values of the ordinary ADF test, as shown above in Section 3. In fact, it is sufficient to set rho2 = 1 to obtain the asymptotic p values of the Dickey-Fuller distribution.For example R> CADFpvalues(-0.44,type = "drift", rho2 = 1) [1] 0.9018844 computes a p value that can be compared directly with the values reported for example in Table 4.2 in Banerjee, Dolado, Galbraith, and Hendry (1993).

Other R implementations of the ADF test
While CADFtest implements Hansen's covariate augmented Dickey-Fuller test and includes the ADF test as a special case, other R packages can perform the ADF test.However, we believe that the function CADFtest() has a more flexible and convenient interface than other existing functions have.
urca (Pfaff 2008) is a leading R package for the analysis of integrated and cointegrated time series that includes the ur.df() function for the ADF test.However, this command cannot deal with missing values and does not have a data argument.Therefore, in order to carry out the same ADF test with three lags we did before, we need to specify all the data details manually, leaving out the first 49 observations for which we have missing values: R> library("urca") R> adf.urca <-ur.df(npext$gnpperca[-(1:49)],type = "trend", lags = 3) Apart from the call being less flexible, the results are not as easy to read as are those that can be obtained from CADFtest().Only the critical values, taken from Dickey and Fuller (1981) and Hamilton (1994), are available to judge the significance of the test, while the test p value is not reported.In fact, if a summary() is performed, the p value associated to the coefficient of the lagged dependent variable is computed using the t distribution that is obviously incorrect under the null.Automatic lag selection is possible on the basis of the AIC and BIC criteria (not the HQ or the MAIC) and a plot() method is available that, when applied to adf.urca, produces a result similar to Figure 2.
As with ur.df(), adf.test() in package tseries (Trapletti 2009) does not allow for missing values.The typical call would be R> library("tseries") R> adf.tseries <-adf.test(npext$gnpperca[-(1:49)], k = 3) The output is easily interpretable, but no summary() or plot() methods are offered.Furthermore, the model is restricted to the case with constant and trend only, and the p value is computed using a simplified procedure based on the interpolation of the values reported in Banerjee et al. (1993, Table 4.2, p. 103).Finally, adf.test() does not offer automatic lag selection options.
fUnitRoots (Wuertz 2009) is another important R package performing unit root tests.In particular, the function unitrootTest() performs the ADF test and computes the relevant finite sample p values using the approach developed in MacKinnon (1996).However, no automatic lag selection is performed.The object passed to the procedure must be a vector or a time series, but no data argument is used, so that the usual call would be something like R> library("fUnitRoots") R> adf.fUnitRoots <-unitrootTest(npext$gnpperca, lags = 3, type = "ct") The function can cope with leading missing values, in the sense that the model is correctly estimated.However, in our example the initial 49 missing values are considered in the length of the series, so that the finite sample p value is incorrect.In order to obtain correct p values, no missing values should be present so that in our case the series should again be adjusted manually.

Summary
This paper presents the R package CADFtest that allows unit root testing using the Covariate Augmented Dickey-Fuller (CADF) test originally proposed in Hansen (1995b).
Differently from the already available routines written in GAUSS and in MATLAB (Hansen 1995a), the present functions are easy to use, do not require the user to modify the programs, and allow the computation of the asymptotic p values of the tests.Beside being extremely useful in general, p values computation opens to the possibility of using the CADF tests in unit root combination tests, for example in the context of macro panels (see e.g., Maddala and Wu 1999;Choi 2001Choi , 2006a;;Costantini et al. 2007).When used to perform conventional ADF tests, CADFtest also encompasses the main features of other existing R packages, with a more flexible and intuitive interface.
CADFtest can be downloaded from the Comprehensive R Archive Network (CRAN) at http: //CRAN.r-project.org/package=CADFtest.