Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
97,620 characters · 18 sections · 135 citation commands
bootUR: An R Package for Bootstrap Unit Root Tests
In this paper, we introduce the \pkg{bootUR} package bootUR for \proglang{R} R, which implements several bootstrap tests for unit roots. Unit root testing is an essential part of any statistical analysis of time series. Given the crucial role of unit root testing in time series analysis, surprisingly few \proglang{R} packages exist that allow for easy and comprehensive unit root testing. The \pkg{bootUR} package aims to fill this gap by offering three major contributions to existing \proglang{R} packages. First, it offers a comprehensive, easy-to-use and reliable set of unit root tests not found as generally in other packages. Second, it offers accurate $p$ values based on bootstrap methods. Third, its functions are not only directly applicable to single time series, but also to datasets consisting of a potentially large set of time series. With these contributions the \pkg{bootUR} package provides practitioners with a single source to fill their unit root testing needs.
Proper handling of unit roots is of paramount importance before commencing any form of analysis on the time series of interest. The by far most important use of unit root tests is therefore as a pre-test to determine whether differencing of the series is needed to eliminate the trend and render the time series stationary. Ignoring unit roots, or stochastic trends, essentially invalidates any subsequent statistical analysis: the stochastic trend, and associated non-decaying dependence of the present on the far past of a series, yields standard inference inapplicable. Probably the most famous consequence of ignoring unit roots is the `spurious regression phenomenon', where one finds seemingly important relations (high $R^2$s and highly significant $t$ statistics) between unrelated time series with stochastic trends. These results have a long history and are well-established and extensively documented in the time series literature. A reader new to unit roots may, for instance, consult Enders08 for a classical textbook treatment of this spurious regression phenomenon as well as the more general problems associated with unit roots.
Currently, unit root tests are scattered across several packages in the \proglang{R} environment for statistical computing and graphics, making it difficult for a practitioner to find and apply an appropriate and reliable test. The most popular unit root test is the classical augmented-Dickey Fuller (ADF) test DickeyFuller79,dickey1981likelihood. Implementations of the ADF test are incorporated in various packages, in particular \pkg{CADFtest} lupi2009unit, \pkg{fUnitRoots} fUnitRoots, \pkg{tseries} tseries, and \pkg{urca} urcabook.\footnote{The \pkg{mleur} package mleurwp also implements the ADF test, but links to \pkg{urca} for this purpose. The package \pkg{uroot} uroot used to have the ADF test implemented but it is no longer supported in the package's current version, hence disregarded from the overview.}
As we will argue in the next section, most `standard' unit root tests, such as the ones implemented in these packages, require practitioners to choose which test to use together with its various model specifications, that may have a major impact on the performance of the unit root tests. As its first major contribution, the \pkg{bootUR} package instead implements the user-friendly union of rejections principle HLT09, HLT12, SmeekesTaylor12 that relieves the user from the burden of having to choose the right specification and performs this task automatically.
With the exception of the HEGY seasonal unit root test in the \pkg{uroot} package uroot, current \proglang{R} implementations of unit root tests rely on asymptotic inference when returning critical values or $p$ values for the unit root test.\footnote{Another exception is the repository \pkg{URT} URT, available on GitHub, which includes bootstrap unit root tests. In the remainder, we only focus on packages that are currently maintained on the Comprehensive \proglang{R} Archive Network (CRAN).} However, unit root tests are very sensitive to size distortions in smaller samples due to, for example, neglected serial correlation Schwert89. Size distortions due to features such as time-varying volatility even persist asymptotically Cavaliere05. Bootstrap unit root tests have therefore become a commonly used alternative to asymptotic inference. The bootstrap approximates the exact distribution of the unit root test statistic by repeatedly drawing new samples from the original sample, thereby capturing the features of the time series of interest that affect the distribution of the test, which ensures that bootstrap unit root tests have accurate size properties under very general conditions. The addition of a comprehensive set of bootstrap tests constitutes the second major contribution of our package. The `boot' in \pkg{bootUR} stands for bootstrap to emphasize that most unit roots tests we provide rely on various bootstrap methods for constructing $p$ values.
Finally, most datasets contain multiple, sometimes even many, time series to be tested for unit roots, often leading practitioners to apply unit root tests to each time series separately. Such a practice does not only suffer from multiple testing issues, rejecting several tests by chance alone, but also disregards similarities between individual time series which, if exploited, could increase the often limited power of the individual tests. Although some packages provide joint unit root tests for multivariate or panel data (\pkg{pdR}, pdR; \pkg{plm}, plmbook, plmJSS),\footnote{The packages \pkg{PANICr} PANICr and \pkg{punitroots} punitroots also provide panel unit root tests, but the former has been removed from CRAN and the latter is only available on \proglang{R-Forge} theussl2009collaborative.} such tests may increase power but do not allow one to determine the properties of individual series. For this goal one would need tests accounting for multiple testing, but proper implementations of multiple testing corrections are currently lacking for unit root tests. Therefore, the third major contribution of \pkg{bootUR} is to implement easy tools for applying unit root tests to multivariate time series with automatic multiple testing control.
With these contributions, the \pkg{bootUR} package provides a unified framework for easy and comprehensive unit root testing based on the following philosophy. 1) for novice users, the tests should be easy to implement with sensible default options; 2) those default options should lead to reliable and accurate unit root tests, applicable in general situations; 3) expert users, familiar with the unit root literature, should be able to easily tweak and adjust the tests to their desired settings; 4) all tests should be easily scalable to large datasets without additional effort by the user, thereby providing `automatic' functionality.
To accomplish our philosophy, the package has a simple structure, yet it offers users a wide variety of unit root tests. In particular, unit root tests can directly be performed on single time series or multiple time series. To this end, we deliberately created separate functions that serve these purposes: the \code{adf()} function provides the traditional `textbook' ADF test DickeyFuller79 based on asymptotic $p$ values, while the functions \code{boot_adf()} and \code{boot_union()} can be used for bootstrap-based inference on single time series. The function \code{boot_ur()} is the main function for multiple time series without multiple-testing control, \code{boot_sqt()} and \code{boot_fdr()} are set up for multiple time series with multiple-testing control, and the \code{boot_panel()} function offers a panel unit root test. For each unit root test, the bootstrap method can be chosen by the end-user. To this end, all functions make use of the universal argument \code{bootstrap}. Via suitable warning and error messages, user-friendly advise is provided on the (non)-applicability of certain bootstrap methods in certain situations. Model specifications (such as deterministic components, lag length selection, detrending methods) are under the user's full control, with the option to have them implemented automatically according the union of unit root tests principle which ensures reliable tests across potentially heterogeneous series. Each function contains many options whose syntax is shared across the package, thereby facilitating usability and control by the end-user.
We have also added several functions, based around the core functions above, that aid in the practical implementation of the unit root tests. Most importantly, the function \code{order_integration()} provides an automatic way to determine the order of integration of each series in a dataset, based on a sequence of one of the aforementioned unit root tests. As it also directly outputs the correctly differenced time series that remove all stochastic trends, it provides the user with the option to conduct the entire unit root pre-analysis with a single command. Additionally, we provide several functions that easily allow the user to assess and visualize properties of the data and outcomes of the tests.
The package is available from the Comprehensive \proglang{R} Archive Network (CRAN) at \url{https://CRAN.R-project.org/package=bootUR}. In addition, the latest (development) version is available on GitHub at \url{https://github.com/smeekes/bootUR}. The core of the package is written in \proglang{C++}, where we make use of the packages \pkg{Rcpp} eddelbuettel2011rcpp,Rcpp13,Rcpp17 and \pkg{RcppArmadillo} eddelbuettel2014rcpparmadillo to facilitate seamless integration with \proglang{R}. In addition we use the \pkg{RcppParallel} RcppParallel package, which provides parallel functionality on all platforms to ensure scalability to large datasets. Version 0.5.0 of the \pkg{bootUR} package and version 4.2.1 of \proglang{R} were used in this paper.
Adhering to the four points of our philosophy not only requires thoughts on how to implement the tests and design the API, but it also requires a careful choice of the appropriate statistical methods. We therefore first consider the problem from a statistical point of view in Section (ref), where we discuss the unit root test for single time series and multiple time series, and in Section (ref), where we discuss the bootstrap methods. We then continue with the package's implementation in Section (ref). Section (ref) uses three empirical applications to compare \pkg{bootUR}'s unit root functions to implementations in other \proglang{R} packages and illustrate its usefulness for practitioners. Section (ref) concludes. Appendix (ref) contains a simulation study of the finite sample performance of the methods included in the package.
We first discuss unit root tests for individual time series (Section (ref)), followed by testing multiple series for unit roots (Section (ref)). In our discussion, paralleling SmeekesWijler20, we do not focus on theory, but on the issues that arise for practitioners when implementing these tests on their time series. For a more extensive and theoretical overview of unit root testing, we refer the interested reader to Choi15.
Consider the case where we have $T$ observations from a time series $y_t$ ($t=1,\ldots,T$) generated according to the data generating process (DGP)
where $d_t$ are deterministic functions of time. In particular, three cases are commonly considered: $d_t=0$ (no deterministic components), $d_t=1$ (intercept only), and $d_t=(1,t)^\top$ (intercept and linear trend). The error process $u_t$ is allowed to be serially correlated and heteroskedastic. The presence of serial correlation in $u_t$ has to be accounted for in inference. Typically, $u_t$ is modelled as an invertible infinite order linear process, for instance as
where $\varepsilon_t$ is typically assumed to be a martingale difference sequence. This linearity motivates the use of adding lagged differences of the time series to account for the serial dependence, as in the classical ADF test DickeyFuller79. However, paparoditis2018the show that ADF-type approaches are valid under much more general forms of dependence in $u_t$.
We focus on testing whether or not $y_t$ contains a unit root, that is on testing
in Equation (ref). Under the null hypothesis of a unit root, $y_t$ contains a stochastic trend, and equivalently $y_t$ is being said to be integrated of order 1 ($I(1)$), while the alternative postulates that $y_t$ is integrated of order 0 ($I(0)$), which is generally taken as synonymous to $y_t$ being stationary. Here `integrated of order $d$' means that $y_t$ should be differenced $d$ times to achieve a process that does not contain a stochastic trend anymore.\footnote{Although stationary is generally used as synonym for $I(0)$, an $I(0)$ process can still be non-stationary, for instance through a shift in the variance. Despite this distinction, we follow tradition and use `$I(0)$' and `stationary' interchangeably.}
To test the null hypothesis of a unit root, the classical ADF test DickeyFuller79,dickey1981likelihood remains the pre-dominant choice in practice. For this reason it also forms the backbone of the \pkg{bootUR} package. However, even in its most basic form, practitioners are required to make several non-trivial choices that affect its performance. Table (ref) summarizes these choices and indicates how the various \proglang{R} packages address each of them. In this section, we first discuss the ADF test and the choices that need to be made, before discussing the union of unit root tests principle proposed by HLT09, HLT12 which alleviates many of the concerns.
\paragraph{ADF test} Let $\Delta$ be the difference operator defined as $\Delta y_t:= y_{t} - y_{t-1}$. If no deterministic components are present, the ADF regression is given by
where the lagged differences of $y_t$ are added to the regression to capture the serial correlation present in $u_t$. Testing the null of a unit root then boils down to testing the significance of the parameter $\gamma$ in Equation (ref).
If the time series $y_t$ is suspected to have deterministic components as well, testing becomes more complicated. The traditional one-step procedure adds the relevant deterministic components directly in Equation (ref). However, the one-step detrending does not lend itself to be extended to other forms of detrending such as the popular generalized least squares (GLS) or quasi-differencing (QD) framework ERS96 discussed later. Most contemporaneous literature ChangPark02,HLT09 therefore considers a two-step approach. Here, a first stage regression is run of $y_t$ on the deterministic components $d_t$, and in the second stage the ADF regression
is run on the residuals of the first stage regression, $y_t^d = y_t - \hat{\beta}^\top d_t$, commonly referred to as the detrended time series. When estimated by ordinary least squares (OLS), the one-step and two-step detrending are generally used interchangeably and will give very similar results.\footnote{Although they are not equal, and may even have different limit behavior, HLT09.} For the traditional \code{adf()} test, \pkg{bootUR} offers the user the choice to implement one-step as an alternative to the default two-step detrending. This option is mainly provided to allow replicability of results obtained through other software or textbook examples. The two-step procedure is implemented for all bootstrap unit root tests in the \pkg{bootUR} package, including the \code{boot_adf()} function for single time series.
ERS96 proposed the DF-GLS test, in which $\hat{\beta}$ is obtained by a GLS regression where the (near-)unit root in $y_t$ is first removed by quasi-differencing; the regression is then performed by OLS for $y_t - \left(1 - \frac{c}{T} \right) y_{t-1}$ on $d_t - \left(1 - \frac{c}{T} \right) d_{t-1}$, where $c$ is a parameter that determines how close to differencing the GLS step is; ERS96 recommend that $c=7$ for the case $d_t = 1$ and $c=13.5$ for the case $d_t = (1,t)^\top$ to yield tests with good power properties. The DF-GLS, or alternatively, DF-QD test,\footnote{To avoid confusion with a `proper' GLS estimation that also takes into account higher-order serial dependence and heteroskedasticty, we follow HLT09 and refer to this test as the quasi-differenced (QD) test rather than GLS.} is often considered to be more powerful than the ADF test, although as shown by MuellerElliott03 and HLT09 inter alia, the QD test is only more powerful if the initial condition, that is the deviation of the start of the time series from equilibrium, is small. When the initial condition is large, the standard OLS-detrended ADF test is considerably more powerful. Among other \proglang{R} packages, only \pkg{urca} provides a version of the DF-QD test, see Table (ref).
While both options are implemented in the function \code{boot_adf()} for varying choice of $d_t$, the issue of including deterministic components presents the practitioner with two difficult choices: which deterministic components to include, and how to perform the detrending. These choices can have a major impact on the performance. If too few deterministic components are included, deterministic trends are detected as stochastic trends, and the test becomes inconsistent. On the other hand, adding too many deterministic components reduces the power of the test considerably, and should also be avoided HLT09. Yet, all unit root tests in current \proglang{R} packages ask the user to make these choices. Similarly, the initial condition is unobservable, such that the user has to choose which detrending method to use.
The \pkg{bootUR} package is, to the best of our knowledge, the first \proglang{R} package which does not force the user to make these choices, but instead offers the function \code{boot_union()} for a data-driven alternative via the union of rejections principle introduced by HLT09, HLT12. Before discussing this in detail, we first turn to the third choice a user has to make: selecting the lag length $p$ in Equation (ref).
The lag length choice concerns a trade off between size distortions incurred from including too few lags to capture all serial correlation, and power loss incurred from including too many lags. Although theory (and some \proglang{R} packages such as \pkg{tseries}) generally assume $p$ to be a deterministic function of the sample size, in practice data-driven selection allows for a more precise choice of $p$ that can achieve a better trade off between size and power. \pkg{bootUR} offers a set of information criteria for automatic data-driven lag length selection. Information criteria trade off model fit (through the residual sum of squares) and overfitting (through a penalty on the number of parameters). The lag length is estimated as
where $\hat{\sigma}^2_k:= (T-p_{\max})^{-1}\sum_{t=p_{\max}+1}^{T}(\hat{\varepsilon}_{k,t}^{d})^2$ with $\hat{\varepsilon}_{k,t}^{d}$ the OLS residuals from the ADF regression with lag length $k$ in Equation (ref), and $C_T$ is a penalty function that differs according to the information criterion used. We consider two penalties: one corresponding to the Akaike information criterion (AIC; $C_T = 2$) and the other to the Bayesian information criterion (BIC; $C_T = \text{ln}T$).
Next to the original criteria, \pkg{bootUR} also implements their modified variants proposed by NgPerron01. These modifications are specifically motivated for lag length selection in the ADF regression. They are given by
where $\xi_T(k):= (\hat{\sigma}^2_k)^{-1} \hat{\gamma}^2\sum_{t=p_{\max}+1}^{T}(y_{t-1}^d)^2$. The lag length is then estimated as in Equation (ref), with $IC(k)$ replaced by $MIC(k)$. The modified AIC (MAIC) is obtained by taking $C_T=2$, the modified BIC (MBIC), by taking $C_T=\text{ln}T$. NgPerron01 show that the MICs yield large size improvements over the ICs for the purpose of unit root testing. perron2007simple recommend to always use the MICs with the OLS rather than QD-detrended data (even if the unit root test itself makes use of QD detrending) since this improves the test's power properties; \pkg{bootUR} follows this recommendation. In addition, there are various seemingly minor aspects of how the lag selection is implemented that influence its performance, such as how many observations are used to calculate the residual sum of squares. NgPerron05 provide a detailed study and guidelines for these choices; \pkg{bootUR} implements the scheme they recommend as optimal.
CPST15 find that heteroskedasticity affects the performance of information criteria, leading to less accurate choices of $p$ and consequent power loss of the unit root tests. They propose rescaled information criteria, where the time series $y_t$ is rescaled with a nonparametric estimate of its (time-varying) standard deviation, thereby eliminating the heteroskedastictiy. The information criterion is then applied to this rescaled series. These rescaled ICs are generally more powerful in the presence of heteroskedasticity, yet very similar to the original ones without. \pkg{bootUR} therefore performs the rescaling by default (with the option not to consider it) since it is a safe choice and relieves the user of the burden to check whether heteroskedasticity is present.
\paragraph{Union of rejections test} As mentioned above, choosing the right deterministic components to include and the right detrending method to use, is crucial to obtain tests with good power properties. HLT09, HLT12 take a different approach based on a very simple principle. Roughly speaking, for both specification issues, we have one powerful test and one not powerful test. A logical step would therefore be to perform both tests and reject whenever one of them rejects the null hypothesis - the logic being that the one rejecting is then the powerful one. With two tests performed simultaneously, one must control for multiple testing and adjust the tests with a Bonferroni-type adjustment to control size at the desired level. HLT09 introduced this union of rejections idea for the two specification issues separately, while HLT12 combined the two approaches to consider a union of four tests -- intercept only or intercept with trend in combination with OLS or QD detrending -- that guards against both uncertainty over the trend and the initial condition.
This makes the union test a safe option for quick or automatic unit root testing where careful manual specification is not viable, and makes it therefore very suitable for \pkg{bootUR}'s philosophy that the default option provides a reliable and accurate test, for which no in depth knowledge is needed about either the data or the applicability of various unit root tests. Moreover, it scales easily to large datasets with many series, where careful manual considerations about these specifications are not possible regardless of the expertise of the user.
The \pkg{bootUR} package implements the bootstrap version of the union test developed by SmeekesTaylor12, which uses the bootstrap both for determining the appropriate size correction and for obtaining the test's $p$ values. The test statistic takes the form
where $ADF$ and $QD$ are the ADF and QD detrended tests, and superscript $\mu$ and $\tau$ respectively indicating whether the series are demeaned or detrended. The critical values $c_{\cdot}^{\cdot} (\alpha)$ are determined in a preliminary bootstrap step as the individual level $\alpha$ critical values of the four tests; weighting with their inverse is needed to bring the four tests on the same scale. The variable $s$ is a scaling factor to which the statistics are scaled. Any $s < 0$ suffices to preserve the left-tail rejection region; in \pkg{bootUR} we scale to -1. This bootstrap union test is made available through the function \code{boot_union()}. Finally, note that this union-based approach still requires one to select the lag lengths in each of the four ADF regressions. To this end, any of the four information criteria, AIC, BIC, MAIC and MBIC can be used.
Of course, various other unit root tests such as the phillips1988testing (PP) test, the KPSS kwiatkowski1992testing stationarity test, the HEGY test hylleberg1990seasonal for seasonal unit roots, or the covariate-augmented Dickey-Fuller test of hansen1995rethinking exist. While these tests undoubtedly have a lot of value, we intentionally do not implement them in \pkg{bootUR} to avoid overloading the user with choices in tests and options. Instead, we focus on the ADF test which is by far the most popular in practice. Practitioners interested in applying these other unit root tests are recommended to consult the packages \pkg{fUnitRoots}, \pkg{stats}, \pkg{tseries} or \pkg{urca} for the PP test, and \pkg{fUnitRoots}, \pkg{tseries} or \pkg{urca} for the KPSS test. We also want to highlight the \pkg{uroot} package which focuses on the HEGY seasonal unit root test, and the \pkg{CADFtest} package which offers the covariate-augmented ADF test, which exploits correlation with known stationary covariates to improve power. These two packages packages are entirely devoted to these specific unit root tests and already provide excellent functionality for the practitioner in need of these specific tests.
Practitioners often make use of several time series in their analysis, and typically need to test all for unit roots. While performing a unit root test for each series separately is normal practice for a small number of time series, this becomes more complicated if the number of series is large. First, performing many unit root tests simultaneously suffers from multiple testing issues as the probability of incorrect classifications increases with the number of performed tests. Second, we would like to exploit the similarity between different time series to improve the power of the unit root tests, in particular if the time dimension is relatively small.
In the \pkg{bootUR} package, we consider three different ways to approach the testing problem with multiple time series. First, the simplest option of ignoring the test multiplicity issue by just performing unit root tests separately for each series. To this end, the function \code{boot_ur()} from \pkg{bootUR} can be used. Second, we consider the traditional approach of panel unit root tests, where on pools the information in all series to obtain a more powerful test. The function \code{boot_panel()} offers such a test. Third, we can consider individual tests but then with appropriate control of multiple testing error rates. \pkg{bootUR} considers two such tests, namely \code{boot_sqt()} and \code{boot_fdr()}.
Despite the large literature on this topic, software implementations for multiple unit roots are mostly lacking. While there is some support for panel unit root testing as discussed hereafter, methods to control multiple testing in the context of unit root testing are, to the best of our knowledge, not available. While several general purpose multiple testing packages exist, using these in a proper way with unit root tests requires considerable effort and expertise from the user. For instance, some standard corrections may be overly conservative, such as the Bonferroni correction, or only applicable under specific conditions on the dependence, such as BenjaminiHochberg95's BenjaminiHochberg95 method to control the false discovery rate. As argued by for instance RSW08, bootstrap methods for controlling multiple testing allow for general forms of dependence and avoid being too conservative. However, such bootstrap methods need to be integrated with the unit root testing, which is the approach taken in \pkg{bootUR}.
Throughout this section, we use the following notation. Consider $N$ time series for which one would like to test the presence of a unit root. We denote their respective individual unit root test statistics by $UR_i, \ 1\leq i \leq N$. Typically these would correspond to one of the tests discussed in Section (ref). Without loss of generality, we assume that rejections occur for small values of the test statistic.
\paragraph{Panel unit root tests} Panel unit root tests view the multiple time series as a coherent panel dataset, and exploit the similarity between such time series to pool the information in them and achieve more powerful tests. They have a long tradition in econometrics, see e.g., breitung2008unit or Choi15 for reviews. A typical panel unit root test has the null hypothesis that all series have a unit root. Rejection of this null hypothesis is then typically interpreted as evidence that a `significant proportion' of the series is stationary. However, how large that proportion is, or which series are stationary is not revealed by the test. Pesaran12 therefore suggests to use panel unit root tests as an initial screening tool for analyzing multiple series; if the panel unit root test rejects the null, this indicates that the individual series need to examined further; if not, treating the full dataset as $I(1)$ may be a reasonable choice.
We implement the bootstrap Group-Mean ($GM$) test of PSU11
in the function \code{boot_panel()} which is based on averaging the unit root test statistics $UR_i$ ($1\leq i \leq N$) of the $N$ individual time series. This test is valid under very general forms of dependence within the dataset, yet does not require modelling it. Panel unit roots tests are scarcely available for \proglang{R} users. Currently, only two packages with panel unit root tests, namely \pkg{plm} and \pkg{pdR}, are being maintained. The package \pkg{plm} was the first to offer panel unit root tests and provides the tests introduced in maddala1999comparative, choi2001unit, levin2002unit, im2003testing. However, none of them allow for cross-sectional dependence (see kleiber2011panel for a discussion). The package \pkg{pdR} offers the panel unit root test of chang2002nonlinear and a panel version of the seasonal unit root test in hylleberg1990seasonal.
\paragraph{Multiple testing} Multiple testing methods allow practitioners to determine the order of integration for each series in their dataset while controlling the overall probability of making false rejections. First note that in order to properly rank and compare tests for different series, the individual test statistics should have the same marginal distributions. Then, the ranking
corresponds to a ranking from `most significant' to `least significant', when the $i$-th order statistic of $UR_1, \ldots, UR_N$ is denoted by $UR_{(i)}$. To ensure the comparability of these statistics, nuisance parameters need to be eliminated from the distribution of the test statistics. \pkg{bootUR} does this automatically for the union test by scaling all test statistics in Equation (ref) towards $s=-1$; if the user chooses to set specifications manually, it is up to the user to choose them such that any nuisance parameters are eliminated.
The goal is to find an appropriate cut-off point $R$ such that the null of a unit root is rejected for all statistics less than or equal to $UR_{(R)}$, while it is not rejected for all statistics larger. How this threshold is determined, depends on how one controls for multiple testing. \pkg{bootUR} implements two ways to do this: the sequential testing procedure of Smeekes15, which also encompasses the Step-M method of RomanoWolf05 to control the familywise error rate, and the false discovery rate (FDR) controlling approach of RSW08, MoonPerron12.
\paragraph{Sequential quantile test} Smeekes15 proposes the Bootstrap Sequential Quantile Test (BSQT) for multiple unit root testing, that acts as an intermediate between panel unit root testing and full multiple testing control. The method proceeds by sequentially testing groups of time series for unit roots, where the user decides the group sizes. At step 1, we test whether the first $p_1$ series are stationary. Here `first' does not refer to the order in the dataset (which is arbitrary), but to the most significant tests as found via Equation (ref). If the null hypothesis that all $p_1$ units have a unit root cannot be rejected, the test stops. If we do observe a rejection, we move on to the second group where we test if the first $p_2$ are stationary. However, as we already concluded that the first $p_1$ units are stationary, in this second step the actual test is whether the next $p_2 - p_1$ units are stationary as well. We continue this testing procedure until no rejection is observed anymore or we tested all series in the dataset. The BSQT can be performed by using the function \code{boot_sqt()}.
More formally, let $p_1, \ldots, p_K$ be the number of series to be tested as stationary in each of the steps $k=1,\ldots,K$. In the sequential step $k$ we then test
As the first test should have as $H_0$ that all units are $I(0)$, $p_0=0$ by default. Furthermore, $p_K=N$ to complete the testing procedure. The number of steps $K$ and the intermediate numbers $p_1, \ldots, p_{K-1}$ can be chosen by the practitioner. Instead of thinking in terms of $p_k$ series, it may be easier to think in terms of quantiles $q_k$, and set $p_k = [q_k N]$. A practitioner may for instance think “I want to split my series in 10 equally-sized groups.” In that case the practitioner simply sets $q_k = 0.1 k$.
The choice of $\{p_k\}$ does require input and consideration from the user. It can be chosen based on the nature of the dataset and the desired level of precision of the practitioner. Smeekes15 shows that if $p_k$ units are found to be $I(0)$, the probability that the true number of stationary series lies outside the interval $[p_{k-1}, p_{k+1}]$ is at most the chosen significance level of the test. Finding that $p_k$ series are $I(0)$ should therefore be interpreted as finding that the number of $I(0)$ series is in the interval $[p_{k-1}, p_{k+1}]$. In the end, if $p_2, \ldots, p_{K-1}$ are chosen sensibly and not spaced too far apart, the series that lie in the `uncertain interval' are likely those series which are `just about' significant, and correspond to time series with a $\rho$ parameter very close to 1. The practical consequences of incorrect classification of these series are typically small, as their behavior makes them fit reasonably well in both classes of $I(1)$ and $I(0)$ series.
One special case worth mentioning -- set as the default in \code{boot_sqt()} -- is when we set $p_k = k$, such that each series gets tested sequentially. Smeekes15 shows that in this case the BSQT method coincides with the popular Step-M method of RomanoWolf05 to control the familywise error rate (FWE). The FWE is defined as the probability of making at least one false rejection, and is typically controlled via the Bonferroni or holm1979simple approach.
\paragraph{FDR-controlled test} The false discovery rate (FDR), originally proposed by BenjaminiHochberg95, is defined as $FDR = \mathbb{E} \left[\frac{F}{R} \mathbbm{1} (R > 0)\right]$, where $R$ denote the total number of rejections, and $F$ the number of false rejections. It is more appropriate for larger $N$ than the FWE, as it aims to control the proportion of false rejections to the total, rather than the probability of a single false rejection. RSW08fdr develop a bootstrap method to control the FDR, and show that unlike the classical way to control FDR, the bootstrap is appropriate under very general forms of dependence between series. MoonPerron12 applied this method to unit root testing, and it is their method that is implemented in the \code{boot_fdr()} function of the \pkg{bootUR} package. FDR control can be combined with any unit root test specification considered in Section (ref).
To decide on whether to use BSQT or FDR control, relative sample sizes can be considered. The Monte Carlo comparison of Smeekes15 reveals that the FDR-controlling test is somewhat more accurate when the sample size $T$ is at least of equal magnitude as the number of time series $N$, whereas the BSQT method is clearly preferable when $T$ is much smaller than $N$, since the FDR-controlling test then suffers from a lack of power. While Smeekes15 investigates the performance of (various) bootstrap unit root tests based on the standard ADF test with intercept only, we conduct an additional Monte Carlo simulation experiment in Appendix (ref) that compares the asymptotic \code{adf()} test to the bootstrap tests\code{boot_ur()}, \code{boot_sqt()} and \code{boot_fdr()} based on the union of rejections test, which is the default in the package \pkg{bootUR}.
We mainly rely on bootstrap methods to obtain critical values and/or $p$ values for all of the unit root tests discussed in Section (ref). In the \pkg{bootUR} package, six bootstrap methods are implemented: the sieve bootstrap (SB), moving block bootstrap (MBB), sieve wild bootstrap (SWB), dependent wild bootstrap (DWB), block wild bootsrap (BWB) and autoregressive wild bootstrap (AWB). Their properties are summarized in Table (ref), and discussed more extensively below. As immediately apparent from Table (ref), any `off-the-shelf' time series bootstrap method may be used to counteract size distortions arising from neglected serial correlation Schwert89; whereas a wild bootstrap method is needed to deal with general forms of heteroskedasticity CavaliereTaylor08, CavaliereTaylor09ER. General forms of cross-sectional dependence can be captured by any bootstrap method apart from the sieve ones.
Next to correcting the size of unit root tests, bootstrap methods have other advantages. First, the bootstrap offers an automatic $p$ value. Second, the bootstrap directly allows for implementation of multiple testing techniques such as those discussed above. Third, it guards against misspecification and uncertainty regarding the lag length selection in the ADF. As \pkg{bootUR} re-selects the lag lengths within the bootstrap replications, it automatically takes effects of lag selection into account. This, coupled with the fact that the bootstrap captures any dependence missed by the lagged differences in the ADF regression, adds another layer of protection to the tests.
The sieve bootstrap (SB) has been extensively considered in the context of unit root testing; see among others Psaradakis01, chang2003sieve, PapPol05, PSU08 and smeekes2013detrending. It estimates the dependence as an autoregressive (AR) process, resamples the residuals of the AR fit, and then re-applies the AR model recursively to place the dependence back into the bootstrap sample. \pkg{bootUR} determines the required order of the AR model by the order of the ADF model, combining these in a single step as they should conceptually coincide.
While it is able to capture general forms of serial dependence KPP11, it is mostly suited for tests on single time series. SmeekesUrbain14 show that it is not suited to capture general forms of cross-sectional dependence, making it invalid for joint or multiple testing. The \pkg{bootUR} package therefore advises to only use it for unit root testing of a single series or on multivariate series without multiple testing control, throwing a warning to alert the user otherwise. When still applied multivariately, users should also realize that each time series is required to be observed over the same periods, which we refer to as balanced datasets. This often forces practitioners to delete observations for series that have been observed for a longer period. The reason for this limitation is that resampling step of the sieve bootstrap would reshuffle the missing values, creating bootstrap sample with `holes' in it.
The moving block bootstrap (MBB) is another traditional bootstrap method that has not only been used for univariate unit root testing in paparoditis2003residual, but also for multivariate unit root testing in MoonPerron12 and Smeekes15, as well as for panel unit root testing in PSU11. It works by dividing the data in overlapping blocks of data and resampling those blocks to create bootstrap series by laying them end-to-end. The blocks are taken in the time dimension and encompass all series. The block length $\ell$ is set automatically by \pkg{bootUR} as a function of the sample size, following a rule proposed by PSU11 that they showed to perform well in many different circumstances. However, it is easily adjusted by the user to experiment with different lengths and assess the sensitivity of the results for varying block lengths.
The MBB can accommodate any form of serial dependence as long as it `fits' into an adequately sized block, which is a wide class. Unlike the SB, the MBB can also handle general forms of dependence between series, including but not limited to common factor structures. From a practical point of view an attractive feature is that it can be applied without requiring one to model the serial and/or cross-sectional dependence. PSU11 show its validity for mixed $I(1)/I(0)$ panel datasets under such general forms of dependence. The MBB still has, however, two disadvantages: it cannot handle unbalanced datasets and is sensitive to unconditional heteroskedasticity. To handle both issues, users should switch to one of the wild bootstrap methods available in \pkg{bootUR}.
The wild, or multiplier, bootstrap Mammen93, DavidsonFlachaire08 is known to be robust against general forms of heteroskedasticity, however it cannot handle serial dependence. Nonetheless, if combined with a sieve bootstrap, we get the best of both worlds. That is, by replacing the resampling step applied to the residuals of the AR model with a multiplication by independent and identically distributed (iid) random variables with mean zero and variance one, we obtain the sieve wild bootstrap (SWB). CavaliereTaylor09ER, cavaliere2009heteroskedastic and SmeekesTaylor12 among others apply this sieve wild bootstrap for bootstrap unit root testing. The method is perfectly suited to individual unit root testing, but due the AR estimation, suffers from the same inability to capture complex dependence across series as explained by SmeekesUrbain14 for the SB. For the generation of the iid random variables, \pkg{bootUR} follows the unit root papers cited above and uses the normal distribution.
The three remaining bootstrap methods implemented in the package are all wild bootstrap methods adjusted to deal with dependence. However unlike the SWB, here the multiplicative random variables themselves are adjusted to be dependent over time. This setup allows these bootstrap methods to capture complex serial and cross-series dependence structures as well as (unconditional) heteroskedasticity. In addition, no resampling takes place for the dependent wild bootstrap, such that missing values `stay in their place' which makes the method applicable to unbalanced datasets.
The three wild bootstrap methods only differ in how the multiplier variables are made time-dependent. The dependent wild bootstrap method (DWB), originally introduced by Shao10, draws random variables from a $T$-dimensional $N(0, \Sigma)$ distribution, where the elements in $\Sigma$ decrease with the distance between them. Shao10 proposes to use a kernel function to achieve this, along with a bandwidth $\ell$ which ensures that variables more than $\ell$ time points apart are independent. This way $\ell$ has a similar interpretation as the block length in the MBB. RhoShao19 and SmeekesUrbain14RM study the DWB for unit root testing, the latter focusing on multivariate settings.
We consider two more variations. The block wild bootstrap (BWB) shao2011bootstrap, zhang2014bootstrapping is a direct alternative to the MBB, where for each block of size $\ell$, we use the same multiplier variable, and the variables are independent between blocks. The autoregressive wild bootstrap (AWB) SmeekesUrbain14RM, FSU20 generates the multiplier variables as a first-order autoregressive process. Unlike the BWB and DWB who have a block length $\ell$ tuning parameter, the tuning parameter of the AWB is the first-order AR parameter. To be able to use the same tuning parameter $\ell$, we use the conversion formula proposed by SmeekesUrbain14RM and FSU20 that writes the AR parameter as a function of $\ell$, though \pkg{bootUR} also allows to set the AR parameter directly. The default setting for $\ell$ in \pkg{bootUR} uses the same rule as for the MBB, which was also tested for the three wild bootstrap methods by SmeekesUrbain14RM. They also provide theoretical results on the validity of these methods under general forms of dependence and heteroskedasticity.
For completeness, in Algorithm (ref) we present the six bootstrap methods and their role in the general bootstrap algorithm. Note that the outcome of the bootstrap algorithm is a collection of bootstrap unit root test statistics $UR_i^b$ for the series $i=1,\ldots, N$ and bootstrap replications $b=1\ldots,B$. How these are then used depends on the multiple testing approach taken. For instance, if we ignore multiple testing, we simply calculate the bootstrap $p$ values
For the BSQT and FDR tests more involved processing is needed; for details we refer to Smeekes15 and RSW08fdr respectively.
The \pkg{bootUR} package has a simple structure with 13 user-accessible functions. Section (ref) presents three functions to check if the data are suitable to be bootstrapped. Sections (ref) and (ref) introduce the seven core functions for unit root testing on respectively individual and multiple time series. Section (ref) presents three useful functions for determining the order of integration of each series in a particular dataset.
The package's functions will now be presented together with examples of their specific use. To this end, we make use of the dataset \code{"MacroTS"} which contains a collection of 20 macroeconomic time series taken from Eurostat and comes with the package. A complete description of the data can be obtained by typing \code{?MacroTS} in \proglang{R}. The following examples require that the \pkg{bootUR} package and the data have been loaded. As random number generation is required to draw bootstrap samples, we first set the seed of the random number generator to obtain reproducible results:
The package \pkg{bootUR} offers three functions, namely \code{check_missing_insample_values()}, \code{find_nonmissing_subsample()} and \code{plot_missing_values()}, to check if a particular dataset is suitable to be bootstrapped. While the bootstrap tests do not work with missing data, unbalanced datasets are generally allowed (see Table (ref)).
The function \code{check_missing_insample_values()} checks if a particular dataset contains missing values. Its usage is extremely simple, as it only requires the data as input,
which can be a vector, matrix, data frame or in time series format (e.g., \code{ts}, \code{zoo} or \code{xts}). It returns an $N$-dimensional Boolean vector which indicates for each series whether missing values are present (\code{TRUE}) or not (\code{FALSE}).
If a dataset contains series with different starting and end points, the bootstrap methods SWB, DWB, BWB and AWB can still be used. The function \code{find_nonmissing_subsample()} lets users check the start and end points of each series as follows:
The output slot \code{range} returns a $(2\times N)$-matrix displaying the first and last non-missing value for each series, the logical slot \code{all_equal} provides a quick check to see if all time series have the same non-missing indices (\code{TRUE}) or not (\code{FALSE}).
Finally, to display missingness in the dataset, we can use
which displays present cell values in green, missing values at the start or end (`Unbalanced NAs') in purple and internal missing values in red, see Figure (ref). Only the latter are problematic for the wild bootstrap methods, while the purple values also need to be avoided for the resampling-based bootstraps.
\pkg{bootUR} has one function, \code{adf()}, that implements a traditional asymptotic ADF test and two functions to perform a bootstrap unit root test on a single series: \code{boot_adf()} for a standard bootstrap ADF test and \code{boot_union()} for a bootstrap union test.
We start by discussing the standard \code{adf()} function followed by its bootstrap analogue \code{boot_adf()} with the many options users can tweak. As \pkg{bootUR} shares its syntax across the various functions, the majority of function arguments remains identical across \pkg{bootUR}'s functions, which facilitates usability and control by the end-user. In the remainder, we therefore only highlight the differences compared to the \code{adf()} and/or \code{boot_adf()} function.
\paragraph{Asymptotic ADF test} To perform a standard ADF test with asymptotic $p$ value, \pkg{bootUR} offers the \code{adf()} function:
The minimum required input is \code{adf(data)}, where the time series \code{data} can be a vector or a time series object. The optional argument \code{data_name} allows users to name their data for use in the printed output. The remaining set of arguments relates to the ADF regression. The deterministic components can be tweaked via the argument \code{deterministics} with options \code{"none", "intercept", "trend"}, lag length selection in the ADF regression is controlled via the arguments \code{min_lag} and \code{max_lag} that respectively set the minimum and maximum lag length. For the latter, we use the sample size-based default $ \lfloor 12\cdot(T/100)^{1/4} \rfloor$. The information criterion can be selected via the argument \code{criterion} (with options \code{"AIC", "BIC", "MAIC", "MBIC"}) and the argument \code{criterion_scale} lets practitioners choose to use the rescaled information criteria of CPST15. To overwrite data-driven lag selection with a pre-specified lag length, users can simply put both \code{min_lag} and \code{max_lag} equal to the desired lag length. Finally, the argument \code{two_step} allows users to specify whether OLS detrending should be performed following the one-step or two-step procedure, with the latter being the default.
We apply the asymptotic variant of the ADF test with both one-step and two-step OLS detrending to Dutch Gross Domestic Product (GDP). An intercept and linear time trend are added as deterministic components:
The function returns a `\code{bootUR}' class object that inherits from the class `\code{htest}'. The outcome of the unit root test can be easily read from the console via the \code{print()} function\footnote{We slightly adjust the standard \code{R} output from the `\code{htest}' class to make it better fit for our unit root tests. Practitioners interested in retrieving the standard \code{R} output from the `\code{htest}' class, can use \code{stats:::print.htest(adf_out1)} instead.}:
The tests differ slightly numerically but both give the same evidence regarding the presence or absence of a unit root. Return values such as sample estimate $\hat{\gamma}$, test statistic and $p$ value can be easily retrieved via:
Next to the standard slots in `\code{htest}', the object contains the additional slots \code{details} and \code{specifications}. The output slot \code{details} contains an overview of the sample estimate $\hat{\gamma}$ in the slot \code{`individual estimates'}, the test statistic in the slot \code{`individual statistics'}, the $p$ value in \code{`individual p-values'}, and the selected number of lagged differences in the ADF regression in the slot \code{`selected lags`}. The output slot \code{specifications} contains the specifications (deterministics, minimum and maximum lag order, criterion for lag selection, with or without re-scaling and de-trending procedure; one-step or two-step) of the unit root test that is performed.
\paragraph{Bootstrap ADF test} To perform a standard ADF bootstrap unit root test on a single series, the \code{boot_adf()} function can be used. The function is structured as follows:
Many arguments used for \code{adf()} re-appear in the \code{boot_adf()} function and play the same role. We focus the discussion on the new arguments. The arguments related to the bootstrap specifications include the bootstrap method (\code{bootstrap}) and number of bootstrap replications (\code{B}). If a user chooses the bootstrap method \code{"MBB"}, \code{"DBB"} \code{"BWB"} or \code{"AWB"}, the desired block length can be controlled via the argument \code{block_length}. By default, we use \code{block_length =} $\lfloor1.75 \cdot T^{1/3} \rfloor$, as recommended in PSU11. While for the first three, this argument concerns the genuine block length, for the latter, the block length is transformed into an autoregressive parameter \code{ar_AWB} via the formula $0.01^{(1/\texttt{l})}$ as in SmeekesUrbain14RM; this can be overwritten by setting \code{ar_AWB} directly. The following set of arguments relates to the ADF regression and most arguments have a one-to-one correspondence with those in the \code{adf()} function. For \code{boot_adf()}, the user can additionally choose the type of detrending via \code{detrend}, with options \code{"OLS"} (default) and \code{"QD"}.
The argument \code{show_progress} provides live progress updates on the bootstrap. The latter is particularly useful for large values of the argument \code{B}. Finally, the option \code{do_parallel = TRUE} allows the bootstrap to be executed in parallel via \pkg{RcppParallel}; the argument \code{cores} allows users to specify how many cores should be used for the parallel loops. By default, all but two cores are used.
We illustrate the bootstrap ADF test on Dutch GDP, with the sieve bootstrap (\code{bootstrap = "SB"}) as used by PSU08 and smeekes2013detrending. We use the same specification for the deterministic components as for the \code{adf()} test above and detrending is done via OLS.
As for \code{adf()}, the function returns an object of class `\code{bootUR}' and `\code{htest}'. The bootstrap ADF test returns the same value of the sample estimate and test statistic as for the standard ADF test with two-step OLS detrending as it should, only the $p$ value differs. Also in this case, the test indicates no rejection of the unit root.
\paragraph{Union of rejections test} To perform a bootstrap union unit root test on a single series, the \code{boot_union()} function can be used. It shares all its arguments with \code{boot_adf()} except for \code{deterministics} and \code{detrend} which are omitted since \code{boot_union()} performs tests for all the possible options, then combines the outcomes of the four unit root tests as in Equation (ref) to produce a single $p$ value. By combining these tests, the sample estimate is not available and hence set to \code{NA} in the output. The function introduces one new argument \code{union_quantile}, which allows the user to choose the $\alpha$ in Equation (ref) used to construct the critical values that serve as weights in the union. If the user has a particular significance level in mind when performing the test, $\alpha$ should be set to this. If not, the default 0.05 is used.
The bootstrap union test for Dutch GDP with the sieve wild bootstrap as proposed by SmeekesTaylor12 can be obtained via
Finally, for the union test, the output slot \code{details} contains an overview of the estimates, test statistics, $p$ values and selected lag orders of the four tests, namely the ADF unit root tests with intercept only or intercept with trend in combination with OLS or GD trending.
Below, we discuss the various approaches \pkg{bootUR} offers to approach the testing problem with multiple series.
\paragraph{Separate unit root tests} To perform individual ADF tests on multiple time series simultaneously without multiple testing control, the function \code{boot_ur()} can be used:
Compared to the syntax of \code{boot_adf()}, it has three additional arguments. The argument \code{union} controls whether a bootstrap union test is used (\code{TRUE}) or not (\code{FALSE}). If \code{union = TRUE} (default), the arguments \code{deterministics} and \code{detrend} are ignored, and a warning message is returned if the user would have provided specifications for these anyway. If set to \code{FALSE}, the deterministic components and detrending methods can be specified as for the \code{boot_adf()} function. The argument \code{union_quantile} is the same as for \code{boot_union()}. Finally, the optional argument \code{level} offers the user the possibility to give a significance level for the tests to be performed. Doing so adds a slot in the output object which contains a Boolean vector indicating for each series if the unit root null hypothesis is rejected or not. This may be convenient for automated post-processing if the unit root tests are performed as an initial step in an analysis. If a significance level is given, \code{union_quantile} is overwritten with this value.
Furthermore, since the bootstrap is performed for all series simultaneously, the bootstrap methods \code{"SB"} or \code{"MBB"}, that cannot handle unbalanced datasets, should not be used. If the user were to specify these anyway, the function will revert to splitting the bootstrap up and performing it separately for each time series. A warning message is then returned to alert the user. If a vector (or univariate time series) instead of a matrix (multivariate time series) is given for \code{data}, a single unit root test is performed; in this case the function acts as an alternative to \code{boot_adf()} and \code{boot_union()}.\footnote{In fact, internally these functions call \code{boot_ur()}.}
We illustrate the function's usage by performing individual ADF tests with the \verb|"MBB"| bootstrap on the first five series of the unbalanced dataset \code{MacroTS}, which correspond to the real GDP in Belgium, Germany, France, the Netherlands and the United Kingdom respectively.
An object of class `\code{bootUR}' and `\code{mult_htest}', for multiple hypothesis testing, is returned. We constructed the class `\code{mult_htest}' as an extension of `\code{htest}' to make it suitable for multiple hypothesis testing. In addition to the return values as for the class `\code{htest}', it has additional output slots \code{series.names}, which contains the names of the individual time series, and \code{rejections}, which optionally contains the aforementioned vector with rejections in case this is requested.
Looking at the results, none of the time series is stationary at the 5% significance level, as can be seen from the printed $p$ values. Information on, for instance, the test statistics and $p$ values can also be accessed directly via:
\paragraph{Panel unit root test} To perform a panel unit root test, the function \code{boot_panel()} can be used. It shares its syntax with \code{boot_ur()}. Unlike for the latter, usage of the \code{"MBB"} or \code{"SB"} bootstrap methods for a panel unit root test on unbalanced datasets will result in an error-- not a warning --since the unbalancedness cannot be reverted. Therefore, users should switch to one of the wild bootstrap methods or balance their dataset. Besides, sieve bootstrap methods can be used, but they are not suited to capture general forms of dependence across units (see Table (ref)). The code therefore warns users against their usage.
We illustrate the usage of the panel unit root test on the five GDP time series with the \verb|"DWB"| bootstrap of Shao10 and RhoShao19:
Since the null is not rejected at any reasonable level, treating all five GDP series as $I(1)$ is reasonable.
\paragraph{Sequential quantile test} To perform the BSQT for multiple unit root testing, the function \code{boot_sqt()} should be used. It has one additional argument compared to the \code{boot_panel()} function, namely \code{steps} which sets the group sizes. These can either be set in units or in quantiles. To split the series in, for instance, $K$ equally sized groups, use \code{steps = 0:K / K}. By the convention of Smeekes15, the first entry of the vector should be equal to zero, while the second entry indicates the end of the first group, and so on. If the initial zero value or the final value ($N$ or 1 for quantiles) are accidentally omitted, the function automatically adds them back. The default \code{steps = 0:NCOL(data)} corresponds to the Step-M method of RomanoWolf05. Regarding the bootstrap methods, the same warning and error messaging as for the \code{boot_panel()} apply. To set the significance level of the individual tests performed sequentially, the argument \code{SQT_level} can be set.
We illustrate the BSQT on the five GDP series with the \code{"AWB"} (default) bootstrap method of SmeekesUrbain14RM and FSU20:
The details on the number of series $p_k$ to be tested as stationary in step $k$ under $H_0$ and $H_1$ (first two columns) is printed to the console, together with, the test-statistic and $p$ value (last two columns) for each of the sequential steps until no rejection occurs. The latter information is also accessible through the output slot \code{details\$SQT}, information on the (non) rejection of the unit root null for each of the series separately can be accessed via the slot \code{rejections}.
\paragraph{FDR-controlled test} To perform a multiple unit root test by controlling the FDR, the function \code{boot_fdr()} should be used. Its arguments are the same as for the other multivariate unit root tests, with the addition of the argument \code{FDR_level} which sets the FDR level. We illustrate it here with the \code{"BWB"} bootstrap method of shao2011bootstrap and SmeekesUrbain14RM:
The procedure developed by RSW08, MoonPerron12 does not provide $p$ values, therefore critical values are returned instead here.
Finally, \pkg{bootUR} offers three useful functions for determining the order of integration of each series in dataset: \code{order_integration()}, \code{diff_mult()} and \code{plot_order_integration()}. The main function is \code{order_integration()} which applies the `Pantula principle' Pantula89 to determine the order of integration of each series
The argument \verb|max_order| sets the maximum order of integration that should be considered for each series. Generally the default of two should generally suffice, with series of order three or higher only very rarely occurring in practice. The user can choose the unit root test through the argument \code{method} depending on whether a single (\code{"adf"}, \code{"boot_adf"} or \code{"boot_union"}) or a multiple time series (\code{"boot_ur"}, \code{"boot_sqt"} or \code{"boot_fdr"}) is considered. All arguments used in these functions can also be passed on via \code{order_integration()}. As setting the order of integration automatically requires “yes/no” decisions about the presence of unit roots, a significance level has to be set via the argument \code{level}. For \code{"boot_sqt"} or \code{"boot_fdr"} this is automatically transformed to \code{SQT_level} and \code{FDR_level}, respectively.
The Pantula principle works as follows. It starts by setting $d = \verb|max_order| - 1$ and testing for a unit root on the $\Delta^{d} y_t$ series. The series for which the unit root null cannot be rejected are classified as $I(d+1)$ and subsequently removed from the dataset. In the next step, $d = d-1$ and the remaining series are tested and classified accordingly. Under the default \verb|max_order = 2|, this second round involves testing the series in levels and classifying them as either $I(1)$ (if the unit root null is not rejected) or $I(0)$ (if the null is rejected).
The function returns an object of class `\code{bootUR}' and `\code{order_integration}' with two accessible elements. The slot \code{diff_data} contains a matrix whose columns are $\Delta^{d_i}y_{i,t}$ with $d_i$ indicating the order of integration of the $i^{th}$ series $(i = 1, \ldots, N)$. This matrix is generated by the user-accessible function \code{diff_mult(data, d)}, where \code{data} is the original dataset and \code{d} is an $N$-dimensional vector indicating each series' order of integration. It contains the same number of rows as the original dataset (since the default setting \code{keep_NAs = TRUE} in \code{diff_mult()} is used), thereby indicating lost observations as missing. It can be tweaked if a practitioner directly makes use of this function. Note that the object of the input series (e.g., matrix, data frame or time series object) is preserved for the differenced series. The output slot \code{order_int} makes the vector containing the found orders of integration available to the end-user.
Finally, if the argument \code{plot_orders} in the function \code{order_integration()} is set to \code{TRUE}, a plot is provided which displays each series' order of integration. To this end, it uses the function \code{plot_order_integration()} which takes the object that \code{order_integration()} outputs as its input. This function is also made accessible if the end-user wishes to further adjust the display of the variable names, legend and colours through its optional arguments \code{show_names}, \code{show_legend}, \code{names_size}, \code{legend_size} and \code{cols}.
We illustrate the methods on three datasets: the income data used in Smeekes15, the \verb|MacroTS| dataset which comes with the package, and the FRED-QD dataset, which is widely used for macro-economic analysis.
\paragraph{Replication Smeekes15} We replicate the empirical application of Smeekes15 based on income data from the Panel Study of Income Dynamics (PSID) used in pesaran2007simple. The income dataset represents a balanced micro-panel with $T=22$ time periods and $N=181$ units split into three subsamples consisting of college graduates (CLG, $N=58$), high school graduates (HSG, $N=87$), and high school dropouts (HSD, $N=36)$.
We use the functions \code{boot_ur()}, \code{boot_sqt()} and \code{boot_fdr()} with the MBB bootstrap to determine the proportion of stationary units in the PSID dataset with the same specifications as in Smeekes15: the bootstrap replications $B$ are set to 4999, the block size to five, and the lags in the ADF regression are selected using the rescaled MAIC criterion with a maximal lag length of three.
Table (ref) presents the rejection results of a unit root. \code{boot_ur()} and \code{boot_fdr()} indicate that around 30% of the units are stationary. As discussed in Smeekes15, \code{boot_fdr()} returns only few units as stationary and thereby likely suffers from low power due to the small sample size. Finally, in line with Smeekes15, the proportion of stationary units are fairly similar across subsamples.
\paragraph{MacroTS} The \verb|MacroTS| dataset contains $N=20$ macro-economic time series collected from Eurostat (\href{https://ec.europa.eu/eurostat/data/database}{https://ec.europa.eu/eurostat/data/database}) and is included in the package. Quarterly observations from 1995-2019 ($T=100$) are available on GDP, consumption, inflation and unemployment for Belgium, Germany, France, the Netherlands and the United Kingdom. The dataset is unbalanced, see Figure (ref).
\paragraph{FRED-QD} This is a quarterly version of the monthly Federal Reserve Economic Data (FRED) database introduced in McCrackenNg16. It contains $N=248$ macro-economic time series and was imported into \proglang{R} using the commands
This paper uses the data from 1959 Quarter 2 to 2019 Quarter 4 ($T=244$) to avoid possible structural breaks due to the COVID-19 pandemic in 2020. If a researcher wishes to import the up-to-date version of the dataset, \code{2020-06.csv} should be changed to \code{current.csv}. As can be seen from Figure (ref), the dataset contains one internal NA, since the third observation of variable 188 (\verb|UMCSENTx|: Consumer Expectations) is missing while the second observation is not. \pkg{bootUR} cannot handle internal missing values but this can be easily fixed by setting the second observation to \verb|NA|, which results in the first three observations of this variable being `unbalanced NAs' that can be handled by \pkg{bootUR}. The resulting dataset then contains 38 macro-economic indicators with missing values at the start of the sample. Finally, note that all FRED-QD series have been classified into $I(0), I(1), I(2)$ by the transformation codes provided in mccracken2020fred. However, the authors themselves indicate several discrepancies between these codes and the outcome of unit root tests. We therefore use the transformation codes as a benchmark for the classifications obtained through the unit root tests but do not necessarily consider the classification closest to theirs to be the best.
Since some of the macro-economic series are likely to be $I(2)$, we use the \code{order_integration()} function (with its defaults) to implement the Pantula principle. All unit root tests in the \pkg{bootUR} are performed with their default settings, which means that union tests are performed with the \verb|AWB| bootstrap method, and lag length selection is done via the re-scaled MAIC. Throughout this section, a significance level of 5% is used. For \code{boot_sqt()}, the default (i.e., Step-M method) is reported as well as results for evenly spaced 0.1 quantiles (\code{steps = 0:10/10}, for \code{MacroTS}), and 0.05 quantiles (\code{steps = 0:20/20}, for FRED-QD).
We compare \pkg{bootUR}'s unit root tests to the \proglang{R} packages reported in Table (ref). We hereby use the following specifications: For the function \code{CADFtest()} (package \pkg{CADFtest}), we perform ADF regressions with intercept and trend (\verb|type = "trend"|), and lag length selection with MAIC (\verb|criterion = "MAIC"|) thereby considering a maximum of $\lfloor 12\cdot(T/100)^{1/4} \rfloor$ lags, set via the argument \verb|max.lag.y|. These lag length specifications correspond to the defaults used in \pkg{bootUR}. For \code{unitrootTest()} (package \pkg{fUnitRoots}), we perform ADF regressions with intercept and trend (\verb|type = "ct"|). By default, one lagged difference is included. For \code{adf.test()} (package \pkg{tseries}), we use its default settings which implies ADF-regressions with intercept and trend and the number of lags fixed to $\lfloor (T-1)^{1/3}\rfloor$, a deterministic function of the sample size. For \code{ur.df()} (package \pkg{urca}), we use ADF-regressions with intercept and trend (\verb|type="trend"|), lag length selection via AIC (\verb|selectlags = "AIC"|), thereby considering a maximum of $\lfloor 12\cdot(T/100)^{1/4} \rfloor$ lags, set via the argument \verb|lags|. Finally, for \code{ur.ers()} (package \pkg{urca}), we use an intercept and trend for detrending (\verb|model| = "trend"). By default, four lagged differences are included in the ADF-regression. Unlike the other packages, \pkg{urca} only comes with critical values to judge the significance of the unit root test, the $p$ value is not reported, see Table (ref). As discussed in lupi2009unit, the $p$ value reported under \code{summary()} is computed using the $t$ distribution, which is incorrect under the unit root null. Finally, only the packages \pkg{CADFtest} and \pkg{fUnitRoots} can handle missing values, for the other packages, we removed missing values prior to performing the unit root tests.
Before applying the various unit root tests to the two datasets, we perform the \code{boot_panel()} (with default settings) to all series taken in first differences, and to all series in levels. Table (ref) reports the $p$ values of the panel unit root tests. For both datasets, the panel unit root tests on the series in first differences indicates that the unit root null is rejected, thereby indicating that a `significant proportion' of the series is stationary in first differences (hence not $I(2)$). The panel unit root test on the series in levels indicates non-rejection of the unit root null.
To shed further light on the order of integration for each of the individual series, the bootstrap unit root tests are applied and compared to the implementations from other \proglang{R} packages. Figure (ref) presents the obtained orders of integration on the \verb|MacroTS| dataset, Figures (ref) and (ref) on the FRED-QD dataset. Note that on such a large dataset, performing the bootstrap tests -- even twice in \code{order_integration()} if \code{d=2} -- can become computationally intensive. Table (ref) gives the computation times relative to our \code{adf()} function.\footnote{As recorded on a MacBook Pro (macOS 11.6.7) with Quad-Core Intel Core i5 2,4 GHz processor.} While the \code{adf()} function took just over a second to run, it took roughly between 7 and 8 minutes to perform each of the bootstrap tests. There is therefore a price to pay for bootstrap accuracy. A simple way to save computation time, while still preserving reasonable accuracy, is to reduce the number of bootstrap replications, by altering the argument \code{B}.
Globally speaking, most unit root tests agree upon a series' classification into $I(0), I(1), I(2)$, which is comforting. Still, several interesting remarks can be made. First, the results of \code{boot_ur()} are fairly similar to \code{adf()}, \code{boot_sqt()} and \code{boot_fdr()} but it classifies a considerable amount of series as $I(0)$ instead of $I(1)$ on the FRED-QD dataset. This illustrates that ignoring multiple testing can quickly lead to a considerable number of misclassifications on such large datasets. Second, among the \code{boot_sqt()} procedures, the default Step-M method tends to classify more series as $I(2)$ than the other procedure. On the smaller \verb|MacroTS| dataset the two versions of the \code{boot_sqt()} show more agreement than on the larger FRED-QD dataset. Third, \code{boot_fdr()} tends to classify more series as $I(1)$ than the other tests. For a more elaborate discussion of this tendency, we refer the interested reader to SmeekesWijler20. Fourth, among the unit root tests from the other \proglang{R} packages, \code{CADFtest()} produces most similar results to \pkg{bootUR}. The function \code{unitrootTest()} detects far less series as $I(2)$.
This paper presents the \proglang{R} package \pkg{bootUR} that provides a unified framework for bootstrap unit root testing on single and multiple time series. To this end, the package builds upon the popular ADF test with a union of rejections principle. Unlike existing packages on unit root tests, \pkg{bootUR} (i) provides a large collection of easy-to-use, fully-controllable and reliable unit root tests, including the union of rejections test which is set as default to enable quick, automatic unit root testing, (ii) ensures accurate inference through bootstrap methods with easy-to-read output (including $p$ values), (iii) allows for testing the presence of unit roots in datasets containing many time series by relying on fast \proglang{C++} implementations.
We thank the editors and referees for their thorough review and highly appreciate their constructive comments which substantially improved the quality of the package and paper. The first author was financially supported by the Netherlands Organization for Scientific Research (NWO) under grant number 452-17-010, the second author by the European Union's Horizon 2020 research and innovation programme under the Marie Sk\lodowska-Curie grant agreement No 832671. We gratefully acknowledge the comments and checks provided by Robert Adamek, Rui Jorge Almeida, Nalan Ba\c{s}t\"urk, Caterina Schiavoni and Eti\"enne Wijler on earlier versions of the package. All remaining errors are our own.