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.
135,992 characters · 11 sections · 0 citation commands
Detecting Sparse Cointegration
\def\spacingset#1{ {#1}}
\pretitle{\flushleft} \posttitle
\preauthor{\flushleft} \postauthor
{\it Keywords:} Cointegration, High Dimensional Data, Adaptive lasso, Unit Roots.
JEL: C32, C52.
\spacingset{1.5}
Understanding long-run equilibrium relationships among economic and financial variables is a central part of economic modelling. These relationships often result in time series moving together due to common trends, a phenomenon known as cointegration. Cointegration analysis, which identifies if nonstationary variables share common stochastic trends, has typically been performed in low-dimensional settings. Classical methods for detecting cointegration like the Engle-Granger two-step procedure (Engle and Granger (1987)) or Johansen's maximum likelihood approach (Johansen (1991)) are well-suited for environments with few variables. However, with the advent of high-dimensional data where the number of series of interest to an investigation may be large, traditional methods become inadequate for identifying cointegrating relationships or testing for their presence. In high-dimensional settings, standard estimation techniques like least squares also become unstable, even when computationally feasible. Moreover, the lack of oracle knowledge makes it challenging to identify a small subset of cointegration-inducing variables from a large pool of candidates.
The goal of this paper is to propose a simple method to test for the presence of cointegration among a large pool of I(1) candidate series, particularly in cases where cointegration, if present, is assumed to be sparse and one wishes to identify potential cointegration-inducing covariates. The environment is that of a single equation cointegration setting where a target series of interest potentially cointegrates with a small number of series from a large pool of candidates. The key insight of our approach is that we can detect cointegration by examining the properties of residuals from a sparse regularized regression. Consider the problem of determining whether the stock price of a specific constituent of a stock market index, such as the FTSE100 or S&P500, cointegrates with other stocks in the index. This question holds important implications for portfolio diversification, as cointegrated stocks exhibit stable long-run relationships that mitigate risks. Similarly, it may be of interest to determine whether stocks within similar industries share a common stochastic trend. Such analyses require cointegration methods capable of handling a large number of covariate candidates to detect meaningful relationships. The sparse nature of cointegration that we operate under reflects the idea that most variables in high-dimensional economic and financial data will be irrelevant to any given cointegrating relationship. For example, while a stock's price may cointegrate with a handful of other stocks in the market (e.g., stocks in related industries), it is unlikely to share such a relationship with all 99 other constituents in an index like the FTSE100. Identifying a relevant subset of predictors is therefore valuable for interpretation and subsequent analysis.
A primary challenge in this high-dimensional sparse context is the selection of a manageable set of candidate covariates that may form a cointegrating relationship. To address this, our first step employs the adaptive lasso (least absolute shrinkage and selection operator). This method is chosen for its ability to produce parsimonious models while delivering consistent estimates of cointegrating relationships. As we show theoretically and empirically, the adaptive lasso can achieve model selection consistency under favorable conditions or the weaker sure-screening property, and yields super-consistent estimators for non-zero slope coefficients, though these ideal properties may not hold universally when regressors themselves are cointegrated.
The second critical challenge lies in robustly determining whether the residuals from the relationship estimated in the first step are stationary (I(0)), indicating cointegration, or nonstationary (I(1)), implying a spurious regression. Although one may be inclined to invoke existing techniques such as ADF type unit root tests (Engle and Granger (1987), Engle and Yoo (1987), Phillips and Ouliaris (1990)) or KPSS type stationarity tests (Kwiatowski et al. (1992), Shin (1994)), both of these face limitations that tend to amplify in high dimensional contexts (e.g., non-standard limiting distributions that depend on the number of fitted covariates requiring model specific tabulations). In addition and as we show in the sequel, under no-cointegration, the adaptive lasso may spuriously select a non-sparse set of predictors (possibly all), making traditional testing based approaches conceptually difficult to design and implement. In this paper we depart from these testing based methods and propose an information-theoretic approach that avoids reliance on asymptotic distribution-based inferences. Crucially, we demonstrate that this residual-based inference remains valid even when variable selection is imperfect as it will typically be the case when the right-hand side covariates are cointegrated themselves. An important additional advantage of such an approach is its robustness to phenomena such as endogeneity and serial correlation while also being immune to the the number of fitted covariates. The viewing of inferences about stationarity and non-stationarity as a model selection problem has been explored in the context of unit-root detection in Phillips and Ploberger (1996), Phillips (2008) among others and in the context of vector error correction models in Gonzalo and Pitarakis (1998).
This work adds to the growing research on high-dimensional estimation and inference in nonstationary settings. Incorporating modern high-dimensional statistical methods into time series analysis is particularly important for understanding economic data. However, the combination of high dimensionality and nonstationarity presents unique challenges. High dimensionality can lead to misleading results in nonstationary environments (e.g. spurious cointegration), as highlighted in Gonzalo and Pitarakis (1999, 2021) and Onatski and Wang (2018). Recent research has developed theoretical tools to address these challenges. Kock (2016) studied the properties of the adaptive lasso in autoregressive models with unit roots, while Koo et al. (2020) and Lee, Shi and Gao (2022) established precise limiting distributions for lasso-based estimators in high-dimensional predictive regressions with unit-root covariates. More recently, Bykhovskaya and Gorin (2022, 2024) and Onatski and Wang (2018, 2021) have developed high-dimensional extensions of Johansen-type maximum eigenvalue and trace tests for no cointegration in large VAR systems, focusing on system-wide inference for the cointegrating rank. Closer to our single-equation perspective, Smeekes and Wijler (2021) proposed an automated sparse cointegration modelling approach based on penalised error-correction models, and Mei and Shi (2024) analysed lasso's properties in high-dimensional predictive regressions with I(1) predictors and a focus on forecasting.
Building on this literature, our paper develops a simple, practical method for detecting sparse cointegrating relationships in a single-equation regression with a potentially large pool of $p$ I(1) covariates, combining adaptive-lasso–based selection of the active set with an information-criterion–based decision on whether the resulting disturbance is I(0) or I(1). Our theoretical analysis characterizes both the ideal behavior of the adaptive lasso under identifiable conditions and its more robust properties for residual-based inference. Our results are developed in a fixed $p$ asymptotic framework, in line with the oracle property literature on adaptive lasso, while our Monte-Carlo experiments explore finite sample behaviour for moderate values of $p$ across various designs.
The remainder of the paper is structured as follows. Section 2 introduces the theoretical framework. Section 3 presents the two-step approach to detecting cointegration and derives its theoretical properties. Section 4 discusses practical implementation issues, including penalty selection and finite-sample adjustments. Section 5 reports Monte Carlo results, and Section 6 concludes. A supplementary appendix accompanying this paper contains a selection of additional simulation based outcomes.
We consider a single-equation cointegration setting where the target variable \( y_t \) may cointegrate with a subset of predictors drawn from a large pool of \( p \) I(1) series. The operating model is given by
where \( y_t \in \mathbb{R} \) is the target variable, \( \bm{x}_t = (x_{1t}, x_{2t}, \ldots, x_{pt})' \in \mathbb{R}^p \) is a \( p \)-dimensional vector of covariates, \( \bm{\beta} = (\beta_1, \beta_2, \ldots, \beta_p)' \in \mathbb{R}^p \) is the vector of unknown slope coefficients, \( \beta_0 \in \mathbb{R} \) is the intercept term, and \( z_t \in \mathbb{R} \) represents the deviation from the equilibrium relationship.
The covariates \( \bm{x}_t \) are modeled as I(1) processes, possibly correlated across dimensions and with serially correlated disturbances
The deviation term \( z_t \) is modelled as
Letting $\bm \eta_{t}=(u_{t},v_{1t},\ldots,v_{pt})'$, we model these $p+1$ disturbance series as
where $\bm C(L)=\sum_{i=0}^{\infty}C_{i}e_{t-i}$ for $\sum_{i=0}^{\infty}i |\bm C_{i}|<\infty$, $\bm C_{0}=I$ and ${\bm e}_{t} \sim i.i.d.(0,\bm \Sigma_{e})$ with $E\|\bm e\|^{2+\delta}<\infty$ for some positive $\delta$. These assumptions are standard in this literature and essentially ensure that an FCLT holds for the ${\bm \eta}_{t}$ sequence (see Phillips and Durlauf (1986)). Under cointegration, the existence of a long-run equilibrium relationship implies that \( z_t \) is stationary with $|\rho|<1$, even though both \( y_t \) and the components of \( \bm{x}_t \) are individually I(1) processes. If $\rho=1$ instead, ((ref)) is viewed as a spurious regression. By “spurious regression” under $\rho=1$ we mean that there is no vector $\bm b \in \mathbb{R}^p$ such that $y_t - \bm x_t' \bm b$ is $I(0)$, i.e.\ $y_t$ is not cointegrated with any linear combination of the regressors.
To formalise the notion of sparse cointegration, let
denote the set of active covariates inducing cointegration, with \( |S| = s \ll p \). The remaining covariates with \( \beta_j = 0 \) are irrelevant for the cointegrating relationship. Note that throughout this paper we denote the true parameter vector by $\bm\beta^0$ and write $\beta_j^0$ for its $j$th component. When there is no ambiguity, we may drop the superscript and write simply $\bm\beta$. The notion of sparse cointegration posits that the target variable \( y_t \) is cointegrated with only a small number of covariates among the large pool of \( p \) candidates. Formally, this is expressed by assuming that the cardinality of the active set \( S \) is much smaller than \( p \), i.e., \( s \ll p \). This assumption reflects the realistic scenario where most covariates are irrelevant for the equilibrium relationship, as often encountered in applications with high-dimensional economic and financial data. Sparse cointegration implies that \( \bm{\beta} \) is a sparse vector, where
In line with the above description of the notion of sparse cointegration, we let the $s \times 1$ vector ${\bm \beta}_{S}$ collect the parameters whose covariates induce cointegration, when the latter is present. Similarly we let ${\bm x}_{{S},t}$ denote the $s-vector$ of active covariates associated with ${\bm \beta}_{S}$. The parameters associated with the variables that do not actively enter ((ref)) are in turn collected in the $(p-s)$ vector ${\bm \beta}_{{S}^{c}}$ while ${\bm x}_{{S}^{c},t}$ collects the $(p-s)$ inactive series.
Our goal is to assess whether ((ref)) is truly a cointegrating relationship by analyzing the stationarity of a sample counterpart of $z_{t}$ when $p$ is large and only a small unknown number of these series induce cointegration, if the latter is truly present. If sparse cointegration is present for instance we expect the adaptive lasso based residuals to behave like an I(0) process which our approach will be designed to detect. Similarly, if the regression is spurious with $\rho=1$ we expect these residuals to behave like an I(1) process.
Throughout this paper we adopt a fixed-$p$ asymptotic framework with $n \to \infty$. This choice is in line with the oracle-property literature on adaptive lasso (e.g., Zou (2006), Kock (2016)) and allows us to obtain transparent conditions under which our two-stage procedure behaves as intended. We do not claim that fixed-$p$ asymptotics are the only or universally preferable framework in high-dimensional settings. Recent contributions also focused on establishing increasing-dimension asymptotics for lasso-type estimators in nonstationary environments. In Smeekes and Wijler (2021) and Mei and Shi (2024) for instance, the authors developed complementary high-dimensional asymptotic theories in which $p$ is allowed to grow with $n$, typically under different objectives (e.g., automated ECM specification or predictive regression based inferences). Our notion of “high dimensionality” here refers instead to the large {\it model space} generated by $p$ potential regressors. Even for moderate values of $p$, the number of possible models grows as $2^{p}$, rendering exhaustive or stepwise classical selection methods impractical. Moreover, it is well known that ordinary least squares becomes unstable when $p/n$ is large, even if $p < n$, and that estimating a large number of slope coefficients by least squares typically leads to residuals that accumulate substantial estimation error and are therefore noisy objects for deciding about cointegration. In this sense, our fixed-$p$ asymptotic analysis should be viewed as a tractable baseline that captures the sparse-cointegration structure and is complemented by Monte Carlo experiments that explore finite-sample performance in moderate-$p$ settings that are empirically relevant.
We initially motivate and describe the use of the adaptive lasso as a method to perform model-selection and estimation in a single shot within the cointegrating regression model in ((ref)). The adaptive lasso estimator of the model parameters is given by
where $\lambda_{n}$ is a regularization parameter (penalty term) and the $w_{j}$'s are the adaptive weights associated with the $\beta_{j}'s$. Specifically, given an initial estimator $\widetilde{\beta}_{j}$ (e.g., OLS) these weights are set as $w_{j}=1/|\widetilde{\beta}_{j}|^{\gamma}$ with $\gamma>0$. The intercept $\beta_0$ is not penalised in (ref).
The adaptive lasso is particularly well-suited to the cointegration setting considered in this paper, as it addresses several limitations of the standard lasso. In high-dimensional regression problems, the standard lasso suffers from shortcomings that are particularly detrimental to our goal of consistently estimating the residuals of a sparse cointegrating regression. Unless strong assumptions are imposed (e.g., irrepresentable condition), it does not satisfy the oracle property, which guarantees asymptotic consistency and correct model selection. The adaptive lasso aims to resolve these issues by introducing data-dependent weights that penalize small coefficients more heavily, with the goal of allowing for consistent estimation of large coefficients while effectively shrinking irrelevant predictors to zero. For instance, if $\widetilde{\beta}_{j}$ is near zero, the associated $w_{j}$ will be large, which can help penalize parameters associated with variables that are unlikely to be active.
Unlike the standard lasso, which imposes the same penalty on all coefficients regardless of their magnitude, the adaptive lasso adjusts the penalty weights based on an initial estimate of the coefficients. By applying smaller penalties to coefficients with larger initial estimates, the adaptive lasso allows these coefficients to converge more accurately to their true values, while continuing to shrink the irrelevant coefficients to zero. This feature makes the adaptive lasso particularly suitable for constructing residuals $\hat{z}_{t}$ that mimic the true $z_{t}'s$.
Before proceeding further it is also important to discuss the explicit inclusion of an intercept term and the use of raw (unstandardised) I(1) series in (ref). In typical lasso applications with stationary (I(0)) data, it is common practice to standardize predictors to have zero mean and unit variance. This ensures that the $\ell_{1}$ penalty is applied on a comparable scale across coefficients. For nonstationary I(1) series, however, standard “standardisation” (subtracting a sample mean and dividing by a sample standard deviation) is problematic. Their sample mean and variance are not well-defined in the usual sense, and such a transformation would also complicate the interpretation of cointegrating coefficients, which are defined in terms of the original scales of the variables. Moreover, in a cointegration setting like ours, the intercept $\beta_{0}$ is part of the long-run equilibrium relationship and is not merely a nuisance parameter. Our model in (ref) includes $\beta_{0}$ and estimates it as part of the adaptive lasso procedure outlined in (ref). Crucially, the penalty term $\lambda_{n}\sum_{j=1}^{p}w_{j}|\beta_{j}|$ applies only to the slope coefficients $\beta_{j}$ and not to the intercept. The Karush–Kuhn–Tucker (KKT) conditions ensure that the residuals sum to zero, analogous to ordinary least squares, and by the Frisch–Waugh–Lovell theorem (see Yamada (2017)) this is equivalent to running the adaptive lasso on demeaned variables while keeping the intercept unpenalised.
While the standard lasso is sensitive to the scale of predictors, the adaptive lasso mitigates this issue via its data-based weights, which automatically adjust penalties based on the apparent importance of each variable. The adaptive weights $w_{j}$ adjust the penalty for each coefficient. If an initial estimate $\widetilde{\beta}_{j}$ is large (which may reflect a strong underlying effect and/or the scale of $x_{jt}$), its corresponding weight will be small, leading to less shrinkage for that coefficient; conversely, a small $\widetilde{\beta}_{j}$ leads to a large weight and thus stronger shrinkage. These adaptive data-dependent weights therefore absorb much of the variation in the magnitudes and scales of the I(1) predictors and mitigate the need for prior standardization purely for the purpose of equitable penalization. In what follows, the theoretical properties of the adaptive lasso, including model selection consistency and $n$-consistency of the estimated active slope parameters, are derived in the context of such I(1) regressors without prior standardization.
Locating the optima of a function such as ((ref)) is a convex optimisation problem which guarantees that any local minimum is also a global minimum. The structure of the program however is challenging due to the inclusion of the $\ell_{1}$ norm penalty in the objective function. Note for instance that although the $\ell_{1}$ norm is convex, it is non-differentiable at points where $\beta_{j}=0$. This creates difficulties for standard optimisation techniques such as gradient descent which relies on smoothness. Instead, coordinate descent type of algorithms which avoid the need to compute a full gradient at non-differentiable points are typically considered.
Given $\hat{\bm \beta}^{AL}$ estimated from (ref), we define
as the active set of predictors selected by the adaptive lasso, and in the sequel refer to the covariates associated with this estimated active set as $\bm{x}_{\hat{S},t}$. We next use this estimated active set $\hat{S}$ to form the residuals of interest. These can be obtained using the adaptive lasso estimates directly
noting that the vast majority of the components of $\hat{\bm \beta}^{AL}$ are set to zero by the effect of the $\ell_{1}$ penalization. Under sparse cointegration with $|\rho|<1$, these residuals behave as I(0) processes, while under $\rho=1$ they behave as I(1) processes. Thus for our second-stage procedure, they serve as valid proxies for $z_t$ in terms of their order of integration properties
Given the residual sequence $\hat{z}_{t}$ constructed as above, our next goal is to assess whether or not it contains a unit root in its autoregressive representation. For this purpose consider the following two competing models
Model ${\cal M}_{0}$ imposes I(1) behaviour by fitting an AR($k$) process in first differences to the residuals. In contrast, model ${\cal M}_{1}$ augments this specification with a lagged level $\hat{z}_{t-1}$ and allows the $\hat{z}_{t}$'s to be stationary when $\phi<0$. Because the $\hat{z}_{t}$'s are centered by construction, the inclusion of an intercept term is not crucial for our asymptotic arguments. The order $k$ of the autoregressions is directly linked to the behaviour of the error process $u_{t}$ driving the equilibrium errors as formulated in (ref).
We now view the objective of testing for cointegration as a model selection problem between models ${\cal M}_{0}$ and ${\cal M}_{1}$. Selecting model ${\cal M}_{0}$ indicates absence of cointegration (i.e. (ref) is a spurious regression). In contrast, support for model ${\cal M}_{1}$ implies that the $\hat{z}_{t}$'s behave like an I(0) process so that $y_{t}$ and the set of covariates selected by the adaptive lasso are cointegrated. Letting $\overline{\sigma}^{2}_{0}$ and $\overline{\sigma}^{2}_{1}$ denote the residual variances from ((ref)) and ((ref)), selection between the two models is made using the criteria
where $c_{n}$ is a penalty term and $k\in \{0,1,\ldots,k_{max}\}$. The proposed model-selection based approach is based on comparing $IC_{n,0}(\hat{z})$ with $IC_{n,1}(\hat{z})$ and leads to choosing ${\cal M}_{0}$ if $IC_{n,0}(\hat{z})\leq IC_{n,1}(\hat{z})$ and to choosing ${\cal M}_{1}$ if $IC_{n,0}(\hat{z})>IC_{n,1}(\hat{z})$. Note that this model selection based approach bears strong resemblance with a likelihood ratio type test since the requirement $IC_{n,0}(\hat{z})>IC_{n,1}(\hat{z})$ for {\it rejecting} ${\cal M}_{0}$ reduces to $n \ln \overline{\sigma}^{2}_{0}/\overline{\sigma}^{2}_{1}>c_{n}$. Here the penalty term $c_{n}$ plays a similar role to the critical values used to form the rejection region of likelihood ratio type statistics, although we emphasise that our analysis proceeds entirely via information criteria rather than the null distribution of a test statistic. It is indeed useful to observe that the above decision based inferences essentially rely on the sign of $\Delta IC_{n}(\hat{z}) \coloneqq IC_{n,0}(\hat{z})-IC_{n,1}(\hat{z})$, $\Delta IC_{n}(\hat{z}) \gtrless 0$, rather than its actual distribution. \\
In summary, our two-stage procedure for detecting sparse cointegration operates as follows:
\paragraph{Remark 1.} In (ref) and (ref) the order $k$ of the fitted autoregressions is meant to capture the potential presence of serial correlation in the $u_{t}$'s driving the equilibrium errors. In classical unit-root and cointegration testing, such augmentation is fundamental for achieving nuisance-parameter-free asymptotics for test statistics. In our model-selection context, however, the asymptotic ability of the proposed procedure to distinguish between ${\cal M}_{0}$ and ${\cal M}_{1}$ hinges on the relative orders of the residual sums of squares under I(0) and I(1) specifications, and not on the detailed limiting distribution of any test statistic. Heuristically, if $c_{n}\rightarrow \infty$, correct decisions will typically be ensured provided that $n \ln (\overline{\sigma}^{2}_{0}/\overline{\sigma}^{2}_{1})$ is $O_{p}(1)$ under the appropriate model, regardless of whether its limit is characterised by nuisance parameters induced by serial correlation and/or endogeneity. Nevertheless, the inclusion of a sufficiently rich lag structure will favourably influence finite-sample properties.
\@startsection{subsubsection}{3}{\z@} {1.5ex \@plus 1ex \@minus .2ex} {0.8ex \@plus .2ex} {\normalfont} {Asymptotic Properties of the Adaptive Lasso (Stage 1): Consistency and Support Recovery}
We first analyze the adaptive lasso's behaviour in ideal conditions to gain intuition on its characteristics within our non-stationary environment. Our analysis requires distinguishing across alternative structures on the cointegration properties of $\bm{x}_{t}$ itself, beyond the generic framework introduced in Section 2 and which remained silent about the presence or absence of cointegration within the ${\bm x}_{t}'s$ themselves.
Our propositions below are understood to hold under the stochastic framework described in Section 2 and more specifically (ref)-(ref). For the purpose of establishing the selection and estimation properties of the adaptive lasso under $|\rho|<1$ and $\rho=1$ (Propositions 1 and 2) we will need to operate under additional conditions on the cointegration properties of the ${\bm x}_{t}'s$ themselves, akin to the single-equation environments considered in the classical Engle-Granger type of analyses or the predictive regression literature. Crucially however, our core methodological contribution on the actual detection of sparse cointegration via the proposed model selection approach based on (ref)-(ref) (Proposition 3) in Stage 2 will be seen to remain valid regardless of whether the right hand side components are cointegrated or not (including the possibility of cointegration within the active predictors themselves). Differently put, the environment under which our Proposition 3 holds can be much weaker than the restrictions we impose in Propositions 1 and 2. \\
Proposition 1. Suppose that $y_t$ and $\bm{x}_t$ satisfy the framework of Section 2 with $|\rho|<1$. Let $S = \{j:\beta^0_j\neq 0\}$ denote the active index set, and let $\widehat{\bm\beta}^{AL}$ be the adaptive lasso estimator based on the full regression of $y_t$ on $(\bm x_t)_{t=1}^n$ with weights $w_j = |\widetilde{\beta}^{\,ols}_j|^{-\gamma}$, $\gamma>0$, where $\widetilde{\bm\beta}^{\,ols}$ denotes the full-sample OLS estimator.
Proposition 1 collects a set of basic properties of the adaptive lasso in the cointegrating regression with $|\rho|<1$. We note for instance that on $S$ the adaptive lasso behaves like a version of the OLS estimator and has the same super-consistency as in the unpenalised cointegrating regression. This result holds even if the covariates in ${\bm x}_{S^{c},t}$ are themselves cointegrated, but requires ruling out any form of cointegration involving the components of ${\bm x}_{S,t}$ as in the classical cointegration literature. This is to ensure that the support estimated by the adaptive lasso contains at least all of the truly active covariates. Part a(ii) of Proposition 1 formalises a “sure screening” property. With probability tending to one, no truly active regressor is dropped from the selected model $\hat S$. If cointegration within the ${\bm x}_{t}'s$ is entirely ruled out however and an additional condition holds on the penalty sequence, then full model selection consistency is achieveable as stated in part (b)(i) of the Proposition.
\paragraph{Remark 2.} The intuition behind the model selection consistency and its weaker sure-screening counterpart of the adaptive lasso is best understood from the formulation of the KKT conditions associated with the objective function in ((ref)). Let \[ S_{n,j} \;\coloneqq\; \frac{1}{n}\sum_{t=1}^n x_{j,t}\,\widehat z_t \] denote the normalised score component for regressor $j$, where $\widehat z_t$ are the adaptive lasso residuals. At any adaptive lasso solution the KKT conditions require that for each coordinate $j$ \[ |S_{n,j}| \;\le\; \frac{\lambda_n}{2 n} w_j \quad\text{whenever }\widehat\beta^{AL}_j = 0. \] Intuitively a regressor can be excluded from the selected model only if its sample covariance with the residuals is dominated by the adaptive penalty. When right hand side cointegration is excluded as in part (b) of Proposition 1, we have for $j \in S$, $S_{n,j}=O_{p}(1)$ with its specific continuous asymptotic distribution containing no atom at zero. At the same time $w_{j}=O_{p}(1)$, and $\lambda_{n}w_{j}/n\rightarrow 0$ so that the KKT condition fails and $P(\hat{\beta}^{AL}_{j}=0)\rightarrow 0$. This is the sure screening property whereby active coefficients cannot be set to zero. For the inactives with $j \in S^{c}$, we also have $S_{n,j}=O_{p}(1)$ but $w_{j}=O_{p}(n^{\gamma})$ and the requirement $\lambda_{n} n^{\gamma-1}\rightarrow \infty$ comes into play. Thus the KKT condition holds and all inactives are dropped, yielding full model selection consistency as stated in part (b) of Proposition 1. \\
Before proceeding with the stage 2 inferences we also document the adaptive lasso's behaviour in a spurious regression setting (i.e., under $\rho=1$). This is summarized in Proposition 2 below restricted to an environment where the right hand side regressors are assumed not to be cointegrated. \\
Proposition 2. Suppose that $y_t$ and $\bm{x}_t$ satisfy the framework of Section 2 with $\rho=1$ in ((ref)), so that $z_{t}$ is a random walk and ((ref)) is a spurious regression in the sense that $y_{t}$ is not cointegrated with any linear combination of ${\bm x}_{t}$. Let $S = \{j:\beta^0_j\neq 0\}$ denote the active index set with $|S|=s$ fixed, and let $\widehat{\bm\beta}^{AL}$ be the adaptive lasso estimator based on the full regression of $y_t$ on $(\bm x_t)_{t=1}^n$ with weights $w_j = |\widetilde{\beta}^{\,ols}_j|^{-\gamma}$, $\gamma>0$, where $\widetilde{\bm\beta}^{\,ols}$ denotes the full-sample OLS estimator. Assume that every non-zero linear combination of the ${\bm x}_{t}'s$ is I(1) i.e. there is no cointegration within the active block, within the inactive block, or between the two blocks. Suppose further that the penalty sequence $\lambda_n$ satisfies
Then, with $p$ fixed and $n \rightarrow \infty$, we have (i) $\hat{\bm \beta}^{AL}=O_{p}(1)$, and (ii) for each coordinate $j \in \{1,2,\ldots,p\}$, $P(\hat{\bm \beta}^{AL}_{j}=0)\rightarrow 0$ , and hence with $p$ fixed $P(\hat{S}=\{1,2,\ldots,p\})\rightarrow 1$, $\hat{S} \coloneqq \{ j: \hat{\beta}_{j}^{AL}\neq 0 \}$ } \\
Proposition 2 establishes that the adaptive lasso asymptotically selects all regressors in the spurious regression with a full rank I(1) regressor system (no right hand side cointegration). The penalty cannot shrink any coefficient exactly to zero and the adaptive lasso essentially behaves much like OLS in terms of support selection.
\paragraph{Remark 3.} Our result in Proposition 2 has ruled out any type of cointegration within the ${\bm x}_{t}'s$ in the right hand side. We do not characterize the limiting support of the adaptive lasso in such instances. {\it Internal} cointegration among inactive regressors as considered in Proposition 1(a) for instance would create stationary directions within the design matrix, and we may expect that the adaptive lasso will eliminate redundant predictors even under spurious regression ($\rho=1$) to minimize the weighted penalty. Consequently, while the selected set $\hat{S}$ is likely to remain large, it need not converge to the full set $\{1,\dots,p\}$. Our Proposition 3 below establishes that the resulting residuals $\hat{z}_t$ consistently inherit the true order of integration of $z_t$ (behaving as $I(0)$ under cointegration and $I(1)$ under spurious regression) regardless of this internal regressor structure.
\@startsection{subsubsection}{3}{\z@} {1.5ex \@plus 1ex \@minus .2ex} {0.8ex \@plus .2ex} {\normalfont} {ADF regression based model selection and sparse cointegration detection (Stage 2)}
We now focus on detecting whether the estimated residuals $\{\hat z_t\}$ are consistent with an $I(0)$ or an $I(1)$ process under $|\rho|<1$ and $\rho=1$ respectively. Recall from (ref)–(ref) that our model–selection rule compares two ADF–type specifications for $\hat z_t$: a unit–root–restricted model ${\cal M}_0$ and an unrestricted “stationary” model ${\cal M}_1$, with corresponding information criteria $IC_{n,0}(\hat{z})$ and $IC_{n,1}(\hat{z})$. We select ${\cal M}_0$ if $IC_{n,0}(\hat{z})<IC_{n,1}(\hat{z})$ and ${\cal M}_1$ otherwise. With $\Delta IC_n(\hat z) \coloneqq IC_{n,0}(\hat z) - IC_{n,1}(\hat z)$, $\Delta IC_n(\hat z)<0$ corresponds to selecting the restricted ADF model ${\cal M}_0$ while $\Delta IC_n(\hat z)>0$ corresponds to selecting the unrestricted model ${\cal M}_1$. We interpret $\Delta IC_n(\hat z)>0$ as evidence in favour of sparse cointegration and $\Delta IC_n(\hat z)<0$ as evidence in favour of a spurious regression. The following result shows that this selection rule is asymptotically consistent under very mild conditions on the time–series behaviour of $\{\hat z_t\}$. \\
Proposition 3. (Residual-Based Cointegration Detection). Suppose that $y_{t}$ and ${\bm x}_{t}$ satisfy the framework of Section 2 with $|\rho|\le 1$. Let $\hat z_t$ be the Stage 1 residuals constructed from the adaptive lasso estimator in ((ref)). Assume that the penalty sequence $c_n$ used in (ref)–(ref) satisfies $c_n\to\infty$ and $c_n/n\to 0$. Then, regardless of the cointegration structure among the regressors ${\bm x}_t$, as $n\to\infty$, (i) If $|\rho|<1$, then $P(\Delta IC_n(\hat{z})>0) \to 1$, so the information criterion selects the unrestricted ADF model ${\cal M}_{1}$ and we correctly conclude in favour of sparse cointegration. (ii) If $\rho=1$, then $P(\Delta IC_n(\hat{z})<0) \to 1$ so the information criterion selects the unit-root restricted model ${\cal M}_{0}$ and we correctly conclude that the regression is spurious.\\
Proposition 3 establishes that the residual–based procedure asymptotically selects the stationary ADF model when the true errors are stationary and the unit–root–restricted ADF model when the true errors follow a random walk. Crucial to the validity of Proposition 3 is the requirement that the $\hat{z}_{t}$'s are prediction consistent under cointegration, and under spurious regression they preserve the unit-root behaviour of $z_{t}$ (see Lemma C1 in the appendix).
A criterion such as the BIC with $c_n=\ln n$ satisfies both $c_n\to\infty$ and $c_n/n\to 0$ and is therefore model–selection consistent for the problem of discriminating between ${\cal M}_0$ and ${\cal M}_1$ on the basis of the behaviour of $\hat z_t$. In practice, however, different choices of $c_n$ entail different finite–sample trade–offs between detecting I(0) behaviour and guarding against spurious rejections under I(1).
\paragraph{Remark 4.} Proposition 3 imposes no restrictions on the internal cointegration of the regressors $\bm x_t$. The logic relies on the properties of the subspace spanned by $\bm x_t$ rather than the uniqueness of the parameter vector estimated by the adaptive lasso. Under cointegration ($|\rho|<1$), the adaptive lasso maintains prediction consistency even if $\bm x_t$ is rank deficient, identifying a linear combination that eliminates the trend so that $\hat z_t$ is $I(0)$. Conversely, under spurious regression ($\rho=1$), $y_t$ contains a stochastic trend asymptotically independent of $\bm x_t$. Consequently, no linear combination of the regressors can eliminate this trend, ensuring that $\hat z_t$ remains $I(1)$ regardless of the rank or endogeneity of $\bm x_t$.
In the context of Proposition 3, it is also important to emphasise that our main result is driven by the relative orders of magnitude of the residual sums of squares under ${\cal M}_0$ and ${\cal M}_1$, rather than by the detailed limiting distribution of any unit–root test statistic. The requirement \[ IC_{n,0}(\hat{z}) > IC_{n,1}(\hat{z}) \quad \Longleftrightarrow \quad n \ln \Bigl(\frac{\overline{\sigma}^2_{0}}{\overline{\sigma}^2_{1}}\Bigr) > c_n \] resembles a likelihood–ratio type rejection rule, but our analysis relies only on the fact that $n \ln (\overline{\sigma}^2_{0}/\overline{\sigma}^2_{1})$ is $O_p(1)$ under the relevant model. No nuisance–parameter–free limit distribution is needed. In particular, the generated–regressor aspect of $\hat z_t$ does not affect the asymptotic selection rule. Replacing $z_t$ by $\hat z_t$ perturbs both residual variances $\overline{\sigma}^2_{0}$ and $\overline{\sigma}^2_{1}$ by $o_p(1)$, so that the sign of $\Delta IC_{n}(\hat{z})$ is unchanged asymptotically. This, however, does not mean that any $c_n$ satisfying $c_n\to\infty$ and $c_n/n\to 0$ will provide accurate decision frequencies in finite samples. From Proposition 2, for instance, we can infer that the adaptive lasso will select a large number of I(1) predictors when the true model is spurious (no cointegration). The resulting residuals $\hat z_t$, obtained using a large number of estimated slope parameters may require substantial sample sizes in order to mimic an I(1) process accurately, and the finite–sample behaviour of $n \ln (\overline{\sigma}^2_{0}/\overline{\sigma}^2_{1})$ can be sensitive to over–selection in the first stage.
Penalty term in the adaptive lasso estimation (stage 1)
The implementation of the adaptive lasso optimisation programme in (ref) requires choosing the penalty parameter $\lambda_{n}$ and setting the adaptive weights $w_{j}$. As is common in the model selection literature, the theoretical requirements on $\lambda_{n}$ leave a multitude of valid choices in practice. In finite samples, $\lambda_n$ is treated as a tuning parameter and selected in a data-driven way. In the context of unit-root autoregressions, Kock (2016) proposes to select $\lambda_{n}$ via a BIC criterion. This entails estimating the model via the adaptive lasso across a grid of $\lambda_{n}$ values and picking the one that minimises BIC. We adopt the same principle. Letting \[ RSS(\lambda_{n}) = \sum_{t=1}^{n}\bigl(y_{t}-\hat{\beta}_{0}^{AL}(\lambda_n)-{\bm x}_{t}'\hat{\bm \beta}^{AL}(\lambda_n)\bigr)^{2} \] denote the residual sum of squares obtained using a given $\lambda_{n}$, we define
where $k(\lambda_{n})$ is the number of nonzero coefficients in $\hat{\bm \beta}^{AL}(\lambda_{n})$. The data-based penalty used in our adaptive lasso estimation is then
where $\Lambda$ is a finite grid of candidate values. In practice, it is convenient to base $\Lambda$ on the usual lasso path: one can first compute the largest penalty $\lambda_{\max}$ for which all slope coefficients are shrunk to zero, and then construct a decreasing grid from a small fraction of $\lambda_{\max}$ up to $\lambda_{\max}$ on a logarithmic scale. In our Monte Carlo experiments we implemented this idea using a grid of 20 logarithmically spaced values spanning a wide range of penalties. Although the BIC-selected $\hat\lambda_n^{BIC}$ is random, the theoretical conditions on $\lambda_n$ ensure that, asymptotically, the admissible range of penalties is broad, and our simulations indicate that BIC tuning produces values of $\hat\lambda_n^{BIC}$ consistent with the behaviour predicted by Propositions 1 and 2.
Adaptive lasso weights (stage 1)
Regarding the adaptive weights $w_j$, we follow the standard practice in the adaptive lasso literature (Zou (2006)) and base them on a preliminary estimator $\widetilde{\beta}_j$ of the slope coefficients in (ref), typically the OLS estimator. Concretely, we set $w_{j}=1/|\widetilde{\beta}_{j}^{ols}|^{\gamma}$, $\gamma>0$. Intuitively, coefficients with larger preliminary estimates receive smaller penalties, while coefficients whose preliminary estimates are close to zero are penalised more heavily. This encourages the adaptive lasso to retain truly active regressors and to shrink irrelevant ones exactly to zero. The exponent $\gamma$ controls the strength of the adaptation. In applied work it is most commonly set to $\gamma=1$ or $\gamma=2$. A larger value of $\gamma$ induces more aggressive penalisation of coefficients with small preliminary estimates, allowing for a slower growth rate of $\lambda_n$ to achieve model-selection consistency; a smaller $\gamma$ implies milder adaptation and thus requires a faster growth rate of $\lambda_n$.
Optimisation algorithm for the adaptive lasso problem (stage 1)
The optimisation problem in (ref) is convex, and in practice we compute $(\hat{\beta}_0^{AL},\hat{\bm \beta}^{AL})$ using a standard coordinate-descent implementation of the lasso. Concretely, after computing the initial OLS-based weights $w_{j}$, we rescale the regressors by forming $\widetilde{X}=X \,\mathrm{diag}(w^{-1})$ and fit a plain lasso to $(\widetilde{X},y)$ over the grid $\Lambda$. The resulting lasso estimates, say $\hat{\bm \beta}^{plain}(\lambda)$, are then mapped back to obtain the adaptive lasso counterparts $\hat{\bm\beta}^{AL}(\lambda) = \mathrm{diag}(w)\,\hat{\bm \beta}^{plain}(\lambda)$. This reparametrisation is standard and allows us to exploit existing, numerically stable lasso solvers without having to code a bespoke adaptive lasso routine.
Treatment of $\hat{z}$ in the detection of sparse cointegration (stage 2)
Our two-step cointegration detection procedure hinges on examining the adaptive-lasso residual series $\{\hat{z}_{t}\}$ via the ADF-type regressions in (ref)–(ref). Proposition 2 shows that under no cointegration (spurious regression) the adaptive lasso tends to select nearly all $p$ predictors as active in large samples. Hence for magnitudes of $\rho$ at or near 1, the residual sequence $\hat{z}_{t}$ is likely formed from a highly over-parameterised fit, with many small coefficients. As a result, its sampling variability may deviate from that of a genuine I(1) process in finite samples, possibly yielding a behaviour that wrongly mimics that of an I(0) process. This may make it harder for the BIC–type criterion in Stage 2 to correctly favour the unit–root–restricted ADF model. For this reason, in our empirical implementation we apply a simple capping rule based on which the original residuals are re-estimated. This is purely a finite-sample device as none of our theoretical results need to invoke such a cap under $n\rightarrow \infty$ and $p$ fixed. The specifics of the proposed mechanism are motivated and discussed in Section 5 below.
Penalty sequence and lag length in the $\hat{z}$ based model-selection procedure (stage 2)
A classical BIC type penalty $c_{n}=\ln n$ in (ref)–(ref) may be too weak unless an unrealistically large sample size is available. Once $\hat{z}_{t}$ is overfitted, adding the single “lagged-level” regressor $\hat{z}_{t-1}$ in ${\cal M}_{1}$ can give such a sharp drop in residual variance that BIC's extra penalty term $(\ln n)/n$ may not be enough to compensate for it. As a result, the procedure implemented with a BIC penalty may frequently produce spurious findings of stationarity even when the true $z_{t}$ is I(1). In our empirical implementation we specify the penalty as $c_{n}=c \ln n$, $c>0$, which clearly satisfies the requirements in Proposition 3 while also allowing us to fine-tune outcomes by experimenting across multiple magnitudes of $c$.
In addition to the above considerations, our two-step cointegration detection procedure requires fitting ADF-type regressions as in (ref)–(ref), and this raises the issue of lag length determination. In line with most of the literature on lag-length selection in ADF-type regressions, we estimate $k$ via a BIC criterion implemented on the fitted model in (ref) using a given upper bound $k_{max}$. Note that in the present context the accuracy of $\hat{k}_{BIC}$ is not fundamental for achieving model selection consistency when comparing models ${\cal M}_{0}$ and ${\cal M}_{1}$. This is because the selection consistency result of Proposition 3 does not require serially uncorrelated in (ref)–(ref); rather, it depends on the relative orders of the residual variances under I(0) and I(1) specifications, which are robust to mild serial correlation captured by a finite lag length.
The Monte Carlo study is designed to complement the asymptotic results in Propositions 1–3 and to document their finite–sample behaviour under a range of DGPs. In particular, we aim to address three questions: (i) How well does the adaptive lasso estimate $S$ and the corresponding slope coefficients when the conditions of Proposition 1 hold, and how does its performance deteriorate when the regressors exhibit more complicated cointegration patterns (Scenarios labelled as C0–C3 in the DGP design below)? (ii) In the spurious–regression case $\rho=1$ of Proposition 2, to what extent does the adaptive lasso behave as predicted, and how sensitive is this behaviour to cointegration among the regressors? (iii) For Proposition 3, does the residual sequence $\{\hat z_t\}$ produced by Stage 1 behave as an $I(0)$ or $I(1)$ process in accordance with the true error $z_t$, is $\hat z_t$ a good proxy for $z_t$, and how accurately does the Stage 2 ADF–based model–selection rule discriminate between $|\rho|<1$ and $\rho=1$?
To answer these questions, we collect three groups of finite–sample metrics from each Monte Carlo design.
\paragraph{Model–selection performance (stage 1).} For each replication we record the selected set $\widehat S$, and summarise performance across replications through the following metrics. Selection error rates: the false positive rate (FPR) on $S^c$, the false negative rate (FNR) on $S$, and the false discovery rate (FDR) summarising the proportion of spurious inclusions among selected predictors. Correct set: the proportion of replications in which the adaptive lasso recovers the true active set exactly, $\widehat S = S$. Selected model size: the average cardinality $|\widehat S|$, which reveals whether the method selects parsimonious models under $|\rho|<1$ and inflates towards the full support in the spurious case $\rho=1$, as suggested by Proposition 2.
\paragraph{Residual behaviour and approximation of $z_t$'s behaviour with $\hat z_t$ (stage 2).} To support the residual–based arguments underlying Proposition 3, we examine the behaviour of the true errors $\{z_t\}$ and the adaptive–Lasso residuals $\{\hat z_t\}$ side by side. For each replication we compute: the mean squared difference $n^{-1}\sum_{t=1}^n (\hat z_t - z_t)^2$, the sample correlation between $\{\hat z_t\}$ and $\{z_t\}$, and the ratio of sample variances $\widehat{\mathrm{Var}}(\hat z_t)/\widehat{\mathrm{Var}}(z_t)$. These quantities quantify how well $\hat z_t$ proxies the latent error $z_t$ and thus how severe the generated–regressor effect is in practice.
\paragraph{Residual–based classification between $|\rho|<1$ and $\rho=1$ (stage 2).} Finally, we assess the performance of the Stage 2 ADF regression-based model–selection procedure of Proposition 3. For each replication we compute the information criteria $IC_{n,0}(\cdot)$ and $IC_{n,1}(\cdot)$ for the unit–root–restricted and unrestricted ADF specifications, respectively, both when they are fitted to the true errors $z_t$ (oracle benchmark) and when they are fitted to the estimated residuals $\hat z_t$. We then record the classification frequencies based on $\hat z_t$ i.e. the proportions of replications in which the rule based on $\Delta IC_n(\hat z)=IC_{n,0}(\hat z)-IC_{n,1}(\hat z)$ selects the stationary specification (interpreted as $|\rho|<1$) or the unit–root–restricted specification (interpreted as $\rho=1$), and the oracle classification frequencies i.e. the analogous frequencies when the same decision rule is applied to the unobserved true errors $z_t$, using $\Delta IC_n(z)$.
We consider the single–equation regression
with $p\in\{10,50,100\}$ and a sparse coefficient vector on $S=\{1,\dots,5\}$, so that $s\coloneqq |S|=5$ and $|S^c|=p-5$. The five active covariates are taken as the first five of the $x_{j,t}$'s, and we experiment with two signal strengths on the active set, namely $\beta_{S}^{(\text{weak})} = (0.25,\,0.25,\,0.25,\,0.25,\,0.25)'$ and $\beta_{S}^{(\text{strong})} = (1,\,0.5,\,1.5,\,0.8,\,1)'$. The disturbance $z_t$ follows an AR(1) process
with $\rho\in\{0,0.5,0.8, 0.9, 1\}$ and where $e_{t}$ is a stationary random shock (denoted $u_t$ in Section 2); we relabel it here for notational convenience. When $|\rho|<1$, $z_t$ is stationary and (ref) is a cointegrating regression (sparse cointegration). When $\rho=1$, $z_t$ is a random walk and (ref) is a spurious regression.
Next, we specify the dynamics of the regressors $x_{j,t}$. Although our key objective is to illustrate our theoretical results in Propositions 1-3 we also wish to document what happens when the restrictions related to the potential presence of cointegration among the right hand side covariates are violated. For this purpose we design the dynamics of the $x_{j,t}$'s in a way that allows for flexible cointegration patterns on the right-hand side. We introduce $K$ common I(1) factors $s_t=(s_{1,t},\dots,s_{K,t})'$:
and a stationary idiosyncratic component $u_t=(u_{1,t},\dots,u_{p,t})'$. Each regressor is then defined as
where $g:\{1,\dots,p\}\to\{1,\dots,K\}$ is a group–assignment function that controls the cointegration structure. If two regressors $x_{j,t}$ and $x_{k,t}$ share the same group index $g(j)=g(k)$, their difference $x_{j,t}-x_{k,t}=u_{j,t}-u_{k,t}$ is I(0) and they are cointegrated. If they do not share any factor, no such exact cointegration arises. In all designs we take $K=p$ and construct different scenarios by changing the assignment $g(\cdot)$. This ensures that any desired pattern of cointegration (or absence) can be induced by suitable grouping.
To capture both cross–sectional correlation among shocks and potential endogeneity between regressors and the error term $z_t$, we stack all innovations into the single vector \[ \eta_t :=
\in\mathbb{R}^{K+p+1}, \] and assume
This specification induces: (i) correlation across the I(1) factors via $\text{cov}(\xi_{k,t},\xi_{\ell,t})$, (ii) correlation across the stationary components via $\text{cov}(u_{j,t},u_{k,t})$, and (iii) endogeneity through nonzero covariances between the shocks to $z_t$ and the shocks to the regressors, $\text{cov}(e_t,\xi_{k,t})$ and $\text{cov}(e_t,u_{j,t})$. All processes are initiated at zero, $s_0=0$, $u_0=0$ and $z_0=0$. Despite its simplicity, the specification (ref)–(ref) captures the three key phenomena we wish to study: I(1) regressors with possibly complicated cointegration patterns, serial correlation in the error, and endogeneity between $\{x_{j,t}\}$ and $z_t$. \\
Right hand side cointegration scenarios: The function $g(\cdot)$ in (ref) allows us to generate several distinct designs for the cointegration structure of the right–hand side covariates (recall that the active set is $S=\{1,\dots,5\}$ and the inactive set is $S^c=\{6,\dots,p\}$).
As noted in Section 4, the adaptive lasso's tendency to overfit under spurious regression ($\rho=1$) typically yields residuals that appear stationary in finite samples, particularly when $p$ is large. While Proposition 3 guarantees consistency as $n \to \infty$, the sample size required for $\hat{z}_t$ to mimic an $I(1)$ process in high-dimensional settings can be substantial. To mitigate this finite-sample pathology while adhering to the sparse-cointegration perspective, we introduce a parsimony cap on the number of selected coefficients:$$s_{\mathrm{cap}} = \min\bigl\{p,\, \max\bigl\{0,\ \mathrm{round}(S_{\mathrm{base}} + F_{\log p}\,\ln p)\bigr\} \bigr\},$$where $S_{\mathrm{base}}$ and $F_{\log p}$ are user-specified constants (e.g., $S_{\mathrm{base}}=4$, $F_{\log p}=1.5$). In practice, we implement this by retaining only the $s_{\mathrm{cap}}$ coefficients with the largest absolute values from the first-stage estimation. This procedure acts as a finite-sample stabilization device akin to cardinality-constrained regression (e.g., Blumenthal and Davies1; Bertsimas et al. 2). Crucially, with $p$ fixed and $s_{\mathrm{cap}} \ge s$, the cap has no impact on the asymptotic validity of our results. Under cointegration ($|\rho|<1$), the capped support continues to contain the true active set with probability approaching one. Under spurious regression ($\rho=1$), trimming the support reduces overfitting without eliminating the unit-root component of the residuals, thereby ensuring the robustness of the Stage 2 decision rule.
\paragraph{Remark 5.} Our Stage 2 procedure relies on $\hat z_t$ inheriting the correct order of integration from $z_t$, rather than converging to it in mean-square. Under $|\rho|<1$, Proposition 1 ensures $\hat z_t$ is a consistent $I(0)$ proxy. In the spurious case ($\rho=1$), however, the uncapped adaptive lasso absorbs substantial random-walk variation from $y_t$ (Proposition 2). The capping mechanism addresses this by limiting the variance reduction achievable through spurious correlations, preserving the $I(1)$ signal required for detection without requiring $\hat z_t$ to be a precise reconstruction of $z_t$. In the empirical analysis that follows, we focus on the capped implementation; comparative results for the uncapped version are provided in the Supplementary Appendix.
\@startsection{subsubsection}{3}{\z@} {1.5ex \@plus 1ex \@minus .2ex} {0.8ex \@plus .2ex} {\normalfont} {Adaptive lasso estimation and model selection consistency (stage1)}
{Favourable designs C0 and C1.} Scenarios C0 (no cointegration among regressors) and C1 (cointegration among inactives only) satisfy the conditions of Proposition 1. Monte Carlo results confirm the expected behaviour when $|\rho|<1$. For $p=10$ and $\rho\in\{0,0.5\}$, the adaptive lasso performs nearly perfectly. At $n=500$, error rates are negligible (e.g., in C0 with $\rho=0$, $\Pr(\widehat S=S)\approx 0.98$). As $n$ increases to $1000$, both FPR and FNR essentially vanish. As dimensionality increases to $p=50$ and $p=100$, selection becomes more difficult but the pattern holds. While the exact support recovery probability $\Pr(\widehat S=S)$ declines (e.g., to $\approx 0.35$ for $p=100$), the FNR remains extremely low. Crucially, the active covariates are almost always included, ensuring reliable Stage 2 residuals. Scenario C1 behaves similarly to C0, with slightly higher difficulty due to cointegration among inactives, but again, truly active regressors are rarely dropped. As persistence increases ($\rho=0.8, 0.9$), the procedure begins to reflect the influence of near-unit-root errors, with slightly rising FPRs and model sizes. This is partly driven by the capping mechanism, which mitigates overfitting as $\rho \to 1$.
{Scenario C2: cointegration among actives.} Scenario C2 violates the identification assumption of Proposition 1. As expected, dropping active regressors becomes common. For instance, with $\rho=0$ and $p=10$, FNR rises to $\approx 0.14$, and $\Pr(\widehat S=S)$ drops to $\approx 0.30$. In high dimensions ($p=50, 100$), $\Pr(\widehat S=S)$ is effectively zero. The adaptive lasso typically selects a sparse model including only a subset of the cointegrated actives, reflecting the multiplicity of sparse representations available when regressors share stochastic trends. As $\rho$ increases, performance further deteriorates, with FDR exceeding $0.6$–$0.7$ near the unit root boundary.
{Scenario C3: cointegration within and across $S$ and $S^c$.} Scenario C3 allows cointegration both within the active set and across the active and inactive blocks. The Monte Carlo results confirm that (capped) adaptive lasso struggles in this environment even when $\rho$ is small. For $\rho=0$, $p=10$, $n=500$ we observe large error rates (FPR $\approx 0.21$, FNR $\approx 0.21$) and a Correct–Set probability near zero.
{Behaviour as $\rho\to 1$ and Proposition 2.} Table 1 illustrates the overfitting predicted by Proposition 2. Under $\rho=1$, the selected model size saturates near the imposed upper bound $s_{cap}$ (e.g., 7 for $p=50$). Table 2 further corroborates this by reporting average selection probabilities. Under cointegration ($\rho \le 0.5$), the probability of selecting actives $\overline{\pi}_{A} \approx 1$ while for inactives $\overline{\pi}_{I} \approx 0$. Under spurious regression ($\rho=1$), $\overline{\pi}_{A}$ drops and $\overline{\pi}_{I}$ increases sharply (e.g., to $\approx 0.47$), confirming the tendency to select spurious predictors.
\paragraph{Remark 6.} The counterpart to Table 1 under no capping imposed on the adapted lasso is Table S1 in the supplementary appendix. We note that key differences only kick in as $\rho$ approaches the unit-root boundary (as expected). Under $p=100$ and $\rho=0.9$ for instance, the capped lasso hovers around $\hat{S}\approx 8$ while without capping $\hat{S}\approx 24$. Such overfitting almost certainly distorts the stochastic properties of $\hat{z}_{t}$ in finite samples as illustrated in Table S4 which presents decision frequencies for sparse cointegration detection (see Supplementary Appendix).
\@startsection{subsubsection}{3}{\z@} {1.5ex \@plus 1ex \@minus .2ex} {0.8ex \@plus .2ex} {\normalfont} {Residual analysis: comparing $\hat{z}$ with $z$}
Table 3 reports Monte Carlo diagnostics for the Stage 2 residuals $\hat z_t$ under the capped implementation, for all scenarios $C0$--$C3$.
Low and moderate persistence ($\rho=0, 0.5$). When $z_t$ is stationary, $\hat z_t$ tracks the true error extremely well. Correlations are near unity ($0.94$--$1.00$) and variance ratios are close to one. This confirms that when Stage 1 consistently identifies the active regressors, Stage 2 delivers residuals that converge to the true error process.
Role of cointegration in $x_t$. Comparing scenarios $C0$--$C3$ shows that the internal cointegration structure of $x_t$ has only a second-order impact on residual diagnostics. Differences in MSE and correlation are marginal and shrink with $n$. This supports Proposition 3. The limiting behaviour of Stage 2 residuals is driven by $z_t$ and is largely invariant to regressor cointegration.
High persistence and Unit Roots ($\rho \ge 0.8$). As $\rho \to 1$, approximating $z_t$ becomes difficult. For $\rho \in \{0.8, 0.9\}$, the variance ratio drops below one, indicating "shrinkage" of the residuals. Crucially, the capping mechanism mitigates this: for $p=100, \rho=0.9$, the capped version preserves $\approx 70\%$ of the variance versus $50\%$ for the uncapped version. When $z_t$ is exactly $I(1)$ ($\rho=1$), MSEs become large and variance ratios drop, as expected. However, the capped procedure retains a non-trivial fraction of the unit-root variation, preventing the over-aggressive removal of the trend observed in the uncapped case.
\@startsection{subsubsection}{3}{\z@} {1.5ex \@plus 1ex \@minus .2ex} {0.8ex \@plus .2ex} {\normalfont} {Sparse cointegration detection (stage2)}
Table 4 reports the performance of the Stage 2 ADF model-selection rule.
Interior stationary regime ($\rho \le 0.8$). For $|\rho|<1$ and away from the boundary, the procedure performs nearly perfectly. Across all scenarios and dimensions, both the oracle (using true $z_t$) and residual-based rules select the stationary model with probability $\approx 1$. This confirms that Stage 2 is robust to the presence or absence of regressor cointegration.
Boundary case ($\rho=0.9$). At $n=500$, even the oracle criterion tends to misclassify $\rho=0.9$ as a unit root ($\approx 80\%$ of replications). Interestingly, the residual-based rule (with capping) is more favourable to stationarity, correcting this bias in higher dimensions. At $n=1000$, both oracle and residual-based decisions align with the asymptotic prediction of stationarity.
Spurious regression regime ($\rho=1$). This is central for Proposition 2. The oracle correctly selects the unit-root model with probability $\approx 1$. The residual-based performance depends on $p$ and capping. For low dimension ($p=10$), the rule performs well ($\approx 97$--$99\%$ accuracy). For high dimension ($p=50, 100$), the parsimony cap is essential. The capped rule selects the unit-root model in $80$--$89\%$ of replications (compared to far lower rates for the uncapped version, see Supplementary Appendix). Although strictly lower than the oracle, the cap effectively restores the correct qualitative behaviour, allowing the information criterion to recognise the absence of cointegration in the vast majority of samples.
Robustness. The results in Table 4 show that the cointegration structure of the regressors ($C0$--$C3$) has virtually no impact on the decision frequencies. The procedure correctly identifies spurious regressions and cointegrating relationships alike, provided Stage 1 overfitting is controlled via capping.
We have proposed a two-step procedure for detecting sparse cointegration in high-dimensional settings. The first stage employs the adaptive lasso to identify a subset of covariates and estimate the long-run relationship. We showed that this estimator achieves model selection consistency in classical settings and retains a sure-screening property even when the regressors exhibit internal cointegration.
The second stage reframes the unit-root testing problem as a model selection procedure. By comparing a unit-root-restricted ADF specification against an unrestricted one via an information criterion, we avoid the non-standard, nuisance-parameter-dependent limiting distributions associated with residual-based tests. Our asymptotic analysis confirms that this decision rule is consistent provided the first-stage residuals inherit the correct order of integration.
A crucial theoretical insight is the adaptive lasso's tendency to over-select predictors under spurious regression ($\rho=1$). This overfitting yields residuals that appear artificially stationary in finite samples. To address this, we introduced a finite-sample parsimony cap and a robust penalty term. While not required for the asymptotic theory, these adjustments are important for preventing spurious findings of cointegration when $p$ is large. Monte Carlo experiments confirm that, away from the unit-root boundary, the procedure mimics an oracle observing the true errors. It maintains high power and controlled size across diverse correlation structures. In spurious regression scenarios, the capped implementation effectively limits overfitting, preserving the unit-root signal in the residuals where the uncapped version fails.
Future work may extend this framework to multiple cointegrating relationships, develop a complete post-selection inference theory, and explore the properties of residual-based decision rules under increasing-dimension asymptotics where $p$ grows with $n$.
{\bf Notation.} (i) For a vector $\bm v\in \mathds{R}^{p}$, the sign vector is $sgn(\bm v)=(sgn(v_{1}),\ldots,sgn(v_{p}))'$, where for any $c\in \mathds{R}$, $sgn(c)=1$ if $c>0$, $sgn(c)=0$ if $c=0$, and $sgn(c)=-1$ if $c<0$. (ii) The subdifferential set of $\|\cdot\|_{1}$ at ${\bm v}$ is $\partial \|\bm v\|_{1}\coloneqq (\partial |v_{1}|,\ldots,\partial |v_{p}|)'$, where $\partial |v_{i}|=\{1\}$ if $v_{i}>0$, $\partial |v_{i}|=[-1,1]$ if $v_{i}=0$, and $\partial |v_{i}|=\{-1\}$ if $v_{i}<0$. Thus $\partial|v_{i}|=sgn(v_{i})$ if $v_{i}\neq 0$ and $\partial|v_{i}|=[-1,1]$ if $v_{i}=0$.\\
{Unit root / cointegration limits.} We work under the operating assumptions stated in Section 2. Under these assumptions the following large–sample results are standard in the unit–root and cointegration literature (see, e.g., Phillips and Durlauf (1986)):
for some random stochastically bounded matrices $\bm A,\bm B$. Note that $\bm B$ is bounded only when $z_t$ is $I(0)$. If $z_t$ is $I(1)$, the term in (ref) is $O_p(n)$. Also, the limiting matrix ${\bm A}$ may be rank deficient if the ${\bm x}_{t}'s$ are cointegrated themselves. Partitioning ${\bm x}_{t}$ as $({\bm x}_{S,t}',{\bm x}_{S^{c},t}')'$ we refer to the corresponding blocks of ${\bm A}$ as ${\bm A}_{S.S}$, ${\bm A}_{S.S^{c}}$, ${\bm A}_{S^{c}.S}$, and ${\bm A}_{S^{c}.S^{c}}$. The presence of cointegration solely within the inactive regressors would imply a rank deficient ${\bm A}_{S^{c}.S^{c}}$ while ${\bm A}_{S.S}\succ 0$.
Our proofs of Propositions 1-2 distinguish between two alternative designs for the p I(1) regressors. We consider two cases: (i) there is no right-hand side cointegration within ${\bm x}_{t}$ itself, that is, $\forall a\neq 0: a'{\bm x}_{t} \sim I(1)$. Equivalently $\sum_{t}{\bm x}_{t}{\bm x}_{t}'/n^{2}\Rightarrow \bm A\succ 0$. We will refer to this regime as $\textbf{X-NC}$. (ii) Only the inactive right-hand side predictors inside $S^{c}$ may be cointegrated, but in a way that does not involve the actives (i.e., no cointegration within ${\bm x}_{S,t}$, no cointegration between ${\bm x}_{S,t}$ and ${\bm x}_{S^{c},t}$, but ${\bm x}_{S^{c},t}$ may have an internal cointegrating structure). We will refer to this regime that implies a rank deficient ${\bm A}_{S^{c}.S^{c}}$ as X-IC. Note also that these restrictions will not be needed for establishing the validity of Proposition 3. Finally, and solely for notational simplicity, all our proofs omit the fitted intercept. \\
{Preliminary Lemmas.} We next introduce a series of intermediate results. Lemmas A1 and A2 relate to the asymptotic properties of the initial OLS estimators and their implications for the behaviour of the adaptive lasso weights. Some of these intermediate results follow from existing results in the unit-root and cointegration literature (e.g., Phillips and Durlauf (1986)). We therefore restrict their proofs to their implications for the asymptotic behaviour of the adaptive lasso weights $w_{j}=|\widetilde{\beta}_{j}^{ols}|^{-\gamma}$, which are invoked in the proofs of Propositions 1-2. Lemma A3 in turn recalls the asymptotic properties of the score terms $\sum_{t} x_{j,t}z_{t}$. \\
Lemma A1 (Weights under Cointegration $|\rho| < 1$). Suppose $y_{t}$ and ${\bm x}_{t}$ satisfy the framework of Section 2 with $|\rho|<1$ (stationary errors).
Proof of Lemma A1 Part (1) is a standard result for OLS in cointegrating regressions (super-consistency). We focus on Part (2). Under (X-NC), for $j \in S^{c}$, standard asymptotics imply $n(\widetilde{\beta}_{j}^{ols}-\beta_{j}^{0}) \Rightarrow L_j$. Since $\beta_j^0 = 0$, $n\widetilde{\beta}_{j}^{ols} \Rightarrow L_j$. The limit $L_j$ is a functional of Brownian motions possessing a continuous distribution with $P(L_j=0)=0$. To show $w_j = O_p(n^\gamma)$, observe that $w_j n^{-\gamma} = |n \widetilde{\beta}_j^{ols}|^{-\gamma}$. We must show this quantity is stochastically bounded. For any $\epsilon > 0$, we rely on the continuous mapping theorem and the property $P(L_j=0)=0$. Choose $\delta > 0$ such that $P(|L_j| < \delta) < \epsilon/2$. By weak convergence, there exists $N$ such that for all $n \geq N$, $P(|n \widetilde{\beta}_j^{ols}| < \delta) < \epsilon$. Let $M = \delta^{-\gamma}$. Then: \[ P(w_j n^{-\gamma} > M) = P(|n \widetilde{\beta}_j^{ols}|^{-\gamma} > \delta^{-\gamma}) = P(|n \widetilde{\beta}_j^{ols}| < \delta) < \epsilon. \] Thus, $w_j n^{-\gamma} = O_p(1)$, or $w_j = O_p(n^\gamma)$. \qed \\
Lemma A2 (OLS under Spurious Regression $\rho=1$). Suppose the conditions of Lemma A1(1) hold but the errors are $I(1)$ (i.e., $\rho=1$, so $y_t$ and ${\bm x}_t$ are not cointegrated). Then: (1) $\widetilde{\beta}_{j}^{ols}=O_{p}(1)$ and $\widetilde{\beta}_{j}^{ols}\Rightarrow \xi_{j}$, where $\xi_{j}$ is a non-degenerate random variable, $P(\xi_{j}=0)=0$. (2) $w_{j}=|\widetilde{\beta}_{j}^{ols}|^{-\gamma}=O_{p}(1)$. \\
\noindentProof of Lemma A2. Part (1) is the standard spurious regression result (Phillips, 1986), where the OLS estimator converges to a functional of independent Brownian motions rather than to 0. Part (2) follows from Part (1) and the continuous mapping theorem, as the limit $\xi_j$ has no mass at zero (Phillips, 1986). \qed\\
Lemma A3 (score asymptotics). (1) Under the conditions of Lemma A1(1), for each $j \in S$, the normalized score $S_{n,j} \coloneqq \frac{1}{n}\sum_{t=1}^n x_{j,t}z_{t} \Rightarrow Z_{j}$, where $Z_{j}$ is non-degenerate continuous random variable, $P(Z_{j}=0)=0$. (2) Under $\rho=1$ and (X-NC) assumed to hold, the score term $S_{n,j}=O_{p}(n)$ and $S_{n,j}/n$ has an $O_{p}(1)$ non-degenerate limit with no atom at 0 for each $j$. \\
\noindentProof of Lemma A3. Both parts follow straightforwardly from well-known results in the unit-root and cointegration literature (e.g., Phillips and Durlauf (1986)). \qed \\
Proof of Proposition 1. We proceed under the conditions collected in the statement of Proposition 1. \\
Part a(i) - ($n$-consistency of the adaptive lasso estimator under (X-IC) and $|\rho|<1$): Consider the adaptive lasso objective under the local parameterization ${\bm \beta}={\bm \beta}^{0}+n^{-1} {\bm u}$ for a generic vector ${\bm u}\in \mathbb{R}^{p}$. We let ${\bm u}=({\bm u}_{S}',{\bm u}_{S^{c}}')'$. Using $y_{t}={\bm x}_{t}'{\bm\beta}^{0}+z_{t}$ (recall that for notational simplicity our derivations ignore the intercept) and standard algebra, we write the objective function as
and let $\hat{\bm u}^{AL} = n (\hat{\bm \beta}^{AL} - \bm \beta^0)$ be the minimizer of $\mathcal{L}_n({\bm u})$. Expanding the square in the first term gives
It is now also convenient to substract $\mathcal{L}_n(\bm 0)$ from $\mathcal{L}_n({\bm u})$, consider
and focus on the optimisation of $V_{n}(\bm u)$. By standard FCLT results, $\bm A_n \Rightarrow \bm A$ and $\bm B_n \Rightarrow \bm B$. The assumptions for part (a) of Proposition 1 state that $\bm x_{S,t}$ is not internally cointegrated and not cointegrated with $\bm x_{S^c,t}$. However, $\bm x_{S^c,t}$ may be internally cointegrated. While the full matrix $\bm A$ may therefore be singular due to internal cointegration within the inactive set $S^c$, the assumption of no cross-cointegration ensures that the stochastic trends in the active set $S$ are asymptotically linearly independent of those in $S^c$. This implies that the Schur complement of the active block is positive definite almost surely, i.e., $\bm A_{SS \cdot S^c} \coloneqq \bm A_{SS} - \bm A_{S S^c} (\bm A_{S^c S^c})^{-} \bm A_{S^c S} \succ 0$. Note that the absence of cross-cointegration between active and inactive predictors is key to ensure this property. Mathematically, it implies that the null space of ${\bm A}$ is entirely contained within the subspace spanned by the inactive parameters. This guarantees that the concentrated objective function remains strictly convex with respect to ${\bm u}_{S}$.
We next formalize the asymptotic behaviour of ${V}_{n}({\bm u})$ by considering the two components in (ref) separately.
Bounding the Quadratic Term: we initially consider the behaviour of the quadratic term $Q_n(\bm u)$ with respect to $\bm u_S$. For any fixed $\bm u_{S^c}$, the quadratic form is minimized with respect to $\bm u_S$ when the gradient is zero. Due to the positive definiteness of $\bm A_{SS \cdot S^c}$, the profile quadratic function (concentrated on $\bm u_S$) grows quadratically: \[ \min_{\bm u_{S^c}} \left( \bm u' \bm A \bm u \right) = \bm u_S' (\bm A_{SS \cdot S^c}) \bm u_S \ge c \|\bm u_S\|^2, \] for some $c > 0$. Intuitively, this ensures that regardless of the behavior of the inactive estimates, the active parameter estimates cannot diverge without pushing the objective function to infinity. \\
\noindentBounding the Penalty Term: We decompose the penalty component by active and inactive sets.
Combining the bounds, for sufficiently large $n$ with high probability the objective function satisfies \[ V_n(\bm u) \ge \bm u_S' \bm A_{SS \cdot S^c} \bm u_S - 2 \bm u' \bm B - o_p(1) \|\bm u_S\|. \] Since $\bm A_{SS \cdot S^c} \succ 0$ and $\bm B = O_p(1)$, the quadratic growth in $\bm u_S$ dominates both $-2\bm u' \bm B$ and the vanishing penalty term. As $\|\bm u_S\| \to \infty$, $V_n(\bm u) \to \infty$ regardless of $\bm u_{S^c}$. By convexity, the minimizer $\widehat{\bm u}_S^{AL} = n(\widehat{\bm\beta}^{AL}_S - \bm\beta^0_S)$ must therefore be stochastically bounded: $n(\widehat{\bm\beta}^{AL}_S - \bm\beta^0_S) = O_p(1).$ \\
Part a(ii) - Sure-Screening Property: The proof of the sure-screening property relies directly on the convergence rate established in Part (a)(i). It holds under X-IC (where inactives may be cointegrated) because, as established in Part (a)(i), the lack of cross-cointegration between actives and inactives ensures the active estimator $\widehat{\bm\beta}^{AL}_S$ remains consistent even if the inactive estimator is not unique. We aim to show that for every truly active variable $j \in S$, the probability of it being excluded from the model converges to zero, i.e., $P(\widehat{\beta}^{AL}_j = 0) \to 0$ as $n \to \infty$. From Part (a)(i), we established the super-consistency of the estimator on the active set which implies $\widehat{\beta}^{AL}_j \xrightarrow{p} \beta^0_j, \quad \forall j \in S.$ By the beta-min condition assumed in the proposition, the true coefficients are bounded away from zero: $\min_{j \in S} |\beta^0_j| \ge \underline{b} > 0.$ Consider the event where an active variable is dropped incorrectly, $\{ \widehat{\beta}^{AL}_j = 0 \}$. On this event, the estimation error must be at least as large as the signal \[ \left| \widehat{\beta}^{AL}_j - \beta^0_j \right| = \left| 0 - \beta^0_j \right| = |\beta^0_j| \ge \underline{b} \] which implies $\{ \widehat{\beta}^{AL}_j = 0 \} \subseteq \{ |\widehat{\beta}^{AL}_j - \beta^0_j| \ge \underline{b} \}.$ Taking probabilities \[ P(\widehat{\beta}^{AL}_j = 0) \le P\left( \left| \widehat{\beta}^{AL}_j - \beta^0_j \right| \ge \underline{b} \right). \]
Since $\widehat{\beta}^{AL}_j$ is consistent for $\beta^0_j$, the probability of the estimation error exceeding any fixed constant $\underline{b} > 0$ converges to zero. Thus, $P(\widehat{\beta}^{AL}_j = 0) \to 0$ for each $j \in S$. Finally, since $|S|=s$ fixed, we apply the union bound to the event of missing any active variable \[ P(S \not\subseteq \widehat{S}) = P\left( \bigcup_{j \in S} \{ \widehat{\beta}^{AL}_j = 0 \} \right) \le \sum_{j \in S} P(\widehat{\beta}^{AL}_j = 0) \to 0. \] Consequently, $P(S \subseteq \widehat{S}) \to 1$. \\
Part b(i) - $n$-consistency of the adaptive lasso estimator under (X-NC) and full model selection consistency: We operate under (X-NC), where there is no cointegration among any components of $\bm x_t$. This implies that the limit matrix $\bm A$ is positive definite almost surely. We use the same objective function expansion $V_n(\bm u)$ as in Part (a)(i) but unlike in Part (a)(i) where we relied on a profile argument, here the quadratic form $\bm u' \bm A_n \bm u$ provides a strictly convex lower bound with respect to the entire vector $\bm u$. We focus on bounding the penalty term to ensure that it does not destabilize the quadratic minimum.
Combining these, for sufficiently large $n$ and with high probability \[ V_n(\bm u) \ge \bm u' \bm A_n \bm u - 2 \bm u' \bm B_n - o_p(1)\|\bm u\|. \]
Let $\lambda_{\min}(\bm A)$ denote the minimum eigenvalue of the limit matrix. Since $\bm A \succ 0$, we have $\lambda_{\min}(\bm A) > c$ for some $c > 0$. Thus, the quadratic term grows as $c\|\bm u\|^2$, which asymptotically dominates the linear signal term ($2\bm u' \bm B_n$) and the vanishing penalty term. Because the objective function diverges to $+\infty$ as $\|\bm u\| \to \infty$ and convex, the minimizer $\widehat{\bm u}$ must lie within a stochastically bounded set. This implies $\widehat{\bm u} = n(\widehat{\bm\beta}^{AL} - \bm\beta^0) = O_p(1)$.\\
From Part (a)(ii) (Sure Screening), we know that $P(S \subseteq \widehat{S}) \to 1$. To prove full consistency, it suffices to show that no inactive variables are selected asymptotically, i.e., $P(j \in \widehat{S}) \to 0$ for all $j \in S^c$. Consider an inactive variable $j \in S^c$ (where $\beta^0_j = 0$). Suppose, for the sake of contradiction, that this variable is selected ($\widehat{\beta}^{AL}_j \neq 0$). The KKT stationarity condition must hold:
We analyze the asymptotic order of magnitude of both sides.\\
Right-Hand Side (Penalty): Under X-NC, Lemma A1(2) applies. The initial OLS estimator satisfies $n\widetilde{\beta}_j^{ols} = O_p(1)$ (non-degenerate limit), which implies that the weights explode at rate $n^\gamma$ and $w_j = n^\gamma O_p(1).$ Substituting this into the penalty term \[ \text{RHS} = \frac{\lambda_n}{n} O_p(n^\gamma) = O_p(\lambda_n n^{\gamma-1}). \] By the assumption $\lambda_n n^{\gamma-1} \to \infty$, the RHS diverges to infinity in probability.\\
Left-Hand Side (Score): Using the score decomposition:
Combining these terms, the Left-Hand Side is stochastically bounded: $\text{LHS} = O_p(1).$ Thus the KKT equality (ref) implies an equality between a bounded term ($O_p(1)$) and a diverging term ($O_p(\lambda_n n^{\gamma-1}) \to \infty$). \[ P\left( |\mathcal{S}_{n,j}| = \frac{\lambda_n}{n} w_j \right) \le P\left( O_p(1) \ge \text{diverging} \right) \to 0. \] Therefore, the probability that an inactive variable is selected approaches zero. Since $S^c$ is a fixed finite set, $P(\widehat{S} \cap S^c = \emptyset) \to 1$. Combined with the sure screening property ($P(S \subseteq \widehat{S}) \to 1$), we conclude: $P(\widehat{S} = S) \to 1.$ \qed \\
Proof of Proposition 2. We analyze the asymptotic behavior of the Adaptive lasso objective function scaled by $n^{-2}$. As we ignore the inclusion of the intercept, the estimator is defined as:
Since $\rho=1$, $y_t$ is an $I(1)$ process. By assumption, $\bm x_t$ contains $p$ non-cointegrated $I(1)$ processes. We recall that the standard FCLT implies the following joint weak convergence as $n \to \infty$:
Part (i): Boundedness ($\widehat{\bm \beta}^{AL} = O_p(1)$). Consider the rescaled objective function $\widetilde{Q}_n(\bm \beta) \coloneqq n^{-2} \sum_{t=1}^n (y_t - \bm x_t' \bm \beta)^2 + n^{-2} \lambda_n \sum_{j=1}^p w_j |\beta_j|$. Expanding the quadratic term \[ \widetilde{Q}_n(\bm \beta) = \frac{1}{n^2}\sum_{t=1}^n y_t^2 - 2 \bm \beta' \left( \frac{1}{n^2}\sum_{t=1}^n \bm x_t y_t \right) + \bm \beta' \left( \frac{1}{n^2}\sum_{t=1}^n \bm x_t \bm x_t' \right) \bm \beta + \frac{\lambda_n}{n^2} \sum_{j=1}^p w_j |\beta_j|. \] We analyze the penalty term. From Lemma A2, under spurious regression, the OLS estimator $\widetilde{\beta}_j^{ols} = O_p(1)$ converges to a non-degenerate random variable. Thus, the weights $w_j = |\widetilde{\beta}_j^{ols}|^{-\gamma} = O_p(1)$. Given the assumption $\lambda_n / n \to 0$, it follows that $\lambda_n / n^2 \to 0$. Therefore for any fixed ${\bm \beta}$ $$\text{Pen}_n(\bm \beta) = \left( \frac{\lambda_n}{n^2} \right) \sum_{j=1}^p O_p(1) |\beta_j| \xrightarrow{p} 0.$$ Consequently, the objective function $\widetilde{Q}_{n}(\bm \beta)$ converges pointwise in distribution to the stochastic limit process $\widetilde{Q}(\bm \beta)$ \[ \widetilde{Q}_n(\bm \beta) \Rightarrow \widetilde{Q}(\bm \beta) \coloneqq c - 2 \bm \beta' \bm b + \bm \beta' \bm A \bm \beta. \] Since $\bm A \succ 0$ almost surely, the limit function $\widetilde{Q}(\bm \beta)$ is strictly convex with a unique minimizer $\bm \beta^* = \bm A^{-1} \bm b$. By convexity and the pointwise convergence to $\widetilde{Q}(\bm \beta)$, the Argmax Continuous Mapping Theorem implies that the minimizer of $\widetilde{Q}_n(\bm \beta)$ converges in distribution to the minimizer of $\widetilde{Q}$ \[ \widehat{\bm \beta}^{AL} \Rightarrow \bm \beta^*. \] Since $\bm \beta^*$ is a well-defined random variable (a functional of Brownian motions), it is stochastically bounded. Thus, $\widehat{\bm \beta}^{AL} = O_p(1)$.
Part (ii): Selection Failure ($P(\widehat{\beta}^{AL}_j = 0) \to 0$). We establish that the adaptive lasso estimator asymptotically collapses to the unpenalized OLS estimator, which is non-zero with probability 1. The unpenalized OLS estimator is given by $\widetilde{{\bm \beta}}^{ols} = (\sum \bm x_t \bm x_t')^{-1} \sum \bm x_t y_t$. Using the limits in (ref) and (ref), $\widetilde{\bm \beta}^{ols} \Rightarrow \bm A^{-1} \bm b$. The KKT optimality condition for the adaptive lasso estimator $\widehat{\bm \beta}^{AL}$ is: \[ 2 \sum_{t=1}^n \bm x_t (y_t - \bm x_t' \widehat{\bm \beta}^{AL}) = \lambda_n \bm W \bm s, \] where $\bm W = \text{diag}(w_1, \dots, w_p)$ and $\bm s$ is the subgradient vector with elements $s_j \in \text{sgn}(\widehat{\beta}_j^{AL})$. Rearranging terms \[ 2 \sum_{t=1}^n \bm x_t y_t - 2 \left( \sum_{t=1}^n \bm x_t \bm x_t' \right) \widehat{\bm \beta}^{AL} = \lambda_n \bm W \bm s. \] Given the positive definiteness of the limit matrix ${\bm A}$, for sufficiently large $n$ we can pre-multiply by the inverse of the Gram matrix $(\sum \bm x_t \bm x_t')^{-1}$ \[ 2 \widetilde{\bm \beta}^{ols} - 2 \widehat{\bm \beta}^{AL} = \lambda_n \left( \sum_{t=1}^n \bm x_t \bm x_t' \right)^{-1} \bm W \bm s. \] Solving for the difference between the estimators gives
We analyze the order of the right-hand side of (ref): $( n^{-2} \sum \bm x_t \bm x_t' )^{-1} \Rightarrow \bm A^{-1} = O_p(1)$, $\bm W = O_p(1)$ (as established in Lemma A3), $\bm s$ is bounded by 1 in element-wise magnitude. The scalar factor is $\lambda_n / n^2$. Since $\lambda_n / n \to 0$, it follows that $\lambda_n / n^2 \to 0$. \\
Consequently, the difference vanishes asymptotically: \[ \widehat{\bm \beta}^{AL} - \widetilde{\bm \beta}^{ols} = O_p(1) \cdot o(1) \cdot O_p(1) = o_p(1). \]
Next, let $\xi_j$ denote the $j$-th element of the limit random vector $\bm A^{-1} \bm b$. The distribution of $\xi_j$ is a ratio of functionals of Brownian motions and is absolutely continuous. Therefore, $P(\xi_j = 0) = 0$. Since $\widehat{\beta}^{AL}_j \xrightarrow{p} \xi_j$ and $P(\xi_j = 0) = 0$, it follows that the probability of the estimator being exactly zero converges to zero: $P(\widehat{\beta}^{AL}_j = 0) \to 0$. Thus the sparse penalty is asymptotically insufficient to force any coefficient to zero against the spurious trend. Taking the intersection over the fixed number of predictors $p$, we conclude $P(\widehat{S} = \{1, \dots, p\}) \to 1$. \qed \\
Remark A1. It is instructive to contrast the failure in Proposition 2 with the success in Proposition 1. In the cointegrated setting (Proposition 1), model selection consistency is achieved provided $\lambda_n n^{\gamma-1} \to \infty$. This condition ensures the penalty term dominates the score term, which is $O_p(n)$ (driven by $I(1) \cdot I(0)$ correlations), because the weights explode at rate $n^\gamma$. In the spurious regression setting (Proposition 2), the context changes fundamentally in two ways. First, the score term becomes $O_p(n^2)$ (driven by the spurious correlation of independent $I(1)$ processes). Second, the initial OLS estimator is inconsistent, converging to a random variable rather than zero. Consequently, the adaptive weights do not explode but remain $O_p(1)$. As a result, satisfying the condition $\lambda_n n^{\gamma-1} \to \infty$ is insufficient to suppress the noise. To force the inactive coefficients to zero in this regime, the penalty parameter would need to grow faster than the spurious trend itself (i.e., $\lambda_n > O_p(n^2)$), a rate that is incompatible with the standard tuning assumption $\lambda_n = o(n)$. \\
We now proceed with the Proof of Proposition 3. A key feature of the proposed model selection based approach for detecting sparse cointegration is its robustness to the “internal” cointegration properties of the ${\bm x}_{t}'s$. We initially introduce an additional lemma (Lemma A4) designed to establish the fact that under both $|\rho|<1$ and $\rho=1$ the adaptive lasso residuals $\hat{z}_{t}$ share the same I(0)/I(1)'ness properties as the $z_{t}'s$ themselves. \\
Lemma A4 (Universal Residual Behavior). Let $\hat{z}_t$ be the adaptive lasso residuals. Regardless of internal cointegration in $\bm x_t$: (i) If $|\rho|<1$ (cointegration): $\frac{1}{n} \sum \hat{z}_t^2 = O_p(1)$. (ii) If $\rho=1$ (spurious): $\frac{1}{n^2} \sum \hat{z}_t^2 \Rightarrow V > 0$. \\
Proof of Lemma A4. The adaptive lasso estimator $\widehat{\bm \beta}^{AL}$ is defined as the global minimizer of the penalized objective function. Consequently, the value of the objective function at the estimator must be less than or equal to its value at any other point, including the true parameter vector $\bm \beta^0$ $$\sum_{t=1}^n \hat{z}_t^2 + \lambda_n \sum_{j=1}^p w_j |\widehat{\beta}_j^{AL}| \leq \sum_{t=1}^n (y_t - \bm x_t' \bm \beta^0)^2 + \lambda_n \sum_{j=1}^p w_j |\beta_j^0|.$$ Dropping the non-negative penalty term on the LHS, we obtain the upper bound for the RSS \[ \sum_{t=1}^n \hat{z}_t^2 \leq \sum_{t=1}^n (y_t - \bm x_t' \bm \beta^0)^2 + \lambda_n \sum w_j |\beta_j^0|. \]
We next analyze the order of the RHS terms:
Part (i): Under cointegration ($|\rho|<1$), $y_t - \bm x_t' \bm \beta^0 = z_t$, which is a stationary $I(0)$ process. The sum of squares of a stationary process satisfies $\sum z_t^2 = O_p(n)$. For the penalty, recall that for active parameters ($j \in S$), the weights converge to constants $w_j = O_p(1)$. Since $\beta_j^0$ is fixed and $\lambda_n/n \to 0$, the penalty term is $O_p(\lambda_n) = o_p(n)$. Thus, the total upper bound is $O_p(n) + o_p(n) = O_p(n)$. Dividing by $n$ gives $$\frac{1}{n} \sum_{t=1}^n \hat{z}_t^2 \leq O_p(1).$$ This bound holds universally because even if $\bm \beta^0$ is not unique (due to rank deficiency), there exists at least one valid cointegrating vector for which this inequality holds.
Part (ii): Under spurious regression ($\rho=1$), $y_t$ is an $I(1)$ process. By the Functional Central Limit Theorem (FCLT), $n^{-1/2} y_{\lfloor nr \rfloor} \Rightarrow B_y(r)$ and also $n^{-1/2} \bm x_{\lfloor nr \rfloor} \Rightarrow \bm B_x(r)$, where $B_y$ is a scalar Brownian motion and $\bm B_x$ is a vector Brownian motion. From Proposition 2, we know the estimator is stochastically bounded, $\widehat{\bm \beta}^{AL} = O_p(1)$, and converges in distribution to some random vector $\bm \xi$. The residuals satisfy \[ \frac{\hat{z}_{\lfloor nr \rfloor}}{\sqrt{n}} = \frac{y_{\lfloor nr \rfloor}}{\sqrt{n}} - \left(\frac{\bm x_{\lfloor nr \rfloor}}{\sqrt{n}}\right)' \widehat{\bm \beta}^{AL} \Rightarrow B_y(r) - \bm B_x(r)' \bm \xi \coloneqq Q(r) \] and the limit process $Q(r)$ is non-zero with probability 1. By the Continuous Mapping Theorem, the normalized sum of squared residuals converges to the integral of this non-zero process \[ \frac{1}{n^2} \sum_{t=1}^n \hat{z}_t^2 = \frac{1}{n} \sum_{t=1}^n \left( \frac{\hat{z}_t}{\sqrt{n}} \right)^2 \Rightarrow \int_0^1 Q(r)^2 dr > 0. \] Thus, the residual variance scales with $n^2$, confirming that $\hat{z}_t$ remains $I(1)$. \qed \\
Proof of Proposition 3. Let $\overline{\sigma}^{2}_{0} = n^{-1} \sum_{t} (\Delta \hat{z}_t - \hat{\bm \beta}_{0}'\boldsymbol{w}_t)^2$ and $\overline{\sigma}^{2}_{1} = n^{-1} \sum_{t} (\Delta \hat{z}_t - \hat{\phi}\hat{z}_{t-1} - \hat{\bm \beta}_{1}'\boldsymbol{w}_t)^2$ denote the residual variance estimators under the restricted (${\cal M}_0$, $\phi=0$) and unrestricted (${\cal M}_1$, $\phi \neq 0$) specifications, respectively. The difference in information criteria is defined as
Part (i) - Cointegration ($|\rho| < 1$). By Lemma A4(i), the first-stage residuals satisfy $n^{-1}\sum \hat{z}_t^2 = O_p(1)$. They behave asymptotically as the stationary error process $z_t$. Consider the auxiliary ADF regression. Under the alternative hypothesis of cointegration, the true coefficient on the lagged level is $\phi = \rho - 1 < 0$. The unrestricted model ${\cal M}_1$ produces a consistent estimator $\hat{\phi} \xrightarrow{p} \phi \neq 0$. Consequently, the estimated error variance converges to the variance of the innovations: $\overline{\sigma}^{2}_{1} \xrightarrow{p} \sigma_{\epsilon}^2$. The restricted model ${\cal M}_0$ imposes $\phi=0$, which is a misspecification. The estimator minimizes the variance of the differences $\Delta \hat{z}_t$ without utilizing the mean-reverting level component $\hat{z}_{t-1}$. Thus, $\overline{\sigma}^{2}_{0} \xrightarrow{p} \sigma_{\Delta z}^2 > \sigma_{\epsilon}^2$. The ratio of the variance estimators converges to a constant strictly greater than 1 \[ \frac{\overline{\sigma}^{2}_{0}}{\overline{\sigma}^{2}_{1}} \xrightarrow{p} \frac{\sigma_{\Delta z}^2}{\sigma_{\epsilon}^2} \coloneqq \kappa > 1. \] Substituting into (ref): \[ \Delta IC_n(\hat{z}) \xrightarrow{p} \ln(\kappa) - 0 > 0. \] Thus, the probability of selecting the unrestricted model approaches 1.
Part (ii): Spurious Regression ($\rho = 1$). By Lemma A4(ii), the first-stage residuals are $I(1)$ such that $n^{-2}\sum \hat{z}_t^2$ converges to a positive random variable. The dependent variable in the ADF regressions is $\Delta \hat{z}_t$, which is $I(0)$. Under the null hypothesis of a unit root, the true parameter is $\phi = 0$. The unrestricted model ${\cal M}_1$ estimates this zero parameter. The t-statistic for $\phi=0$, denoted as $t_{\phi}$, corresponds to a residual-based unit root statistic. From standard least-squares algebra \[ n\overline{\sigma}^{2}_{0} - n\overline{\sigma}^{2}_{1} = \frac{n\overline{\sigma}^{2}_{1}}{n-k-1} t_{\phi}^2 \approx \overline{\sigma}^{2}_{1} t_{\phi}^2. \] From conventional unit-root statistics, we have that $\overline{\sigma}^{2}_{1}=O_{p}(1)$ (the variance of the dependent variable $\Delta \hat{z}$ is stable) and we also know that $t_{\phi}^{2}=O_{p}(1)$ (see e.g. Phillips and Ouliaris (1990)). Rearranging for the ratio of variances we can therefore write \[ \frac{\overline{\sigma}^{2}_{0}}{\overline{\sigma}^{2}_{1}} = 1 + \frac{\overline{\sigma}^{2}_{0} - \overline{\sigma}^{2}_{1}}{\overline{\sigma}^{2}_{1}} = 1 + \frac{1}{n} \left( \frac{ RSS_0 - RSS_1 }{ \overline{\sigma}^{2}_{1} } \right) = 1 + \frac{O_p(1)}{n}. \] Using the Taylor expansion $\ln(1+x) = x + O(x^2)$ for small $x$: \[ \ln\left( \frac{\overline{\sigma}^{2}_{0}}{\overline{\sigma}^{2}_{1}} \right) = \frac{O_p(1)}{n} = O_p(n^{-1}). \] Now consider the penalty term in (ref): \[ \Delta IC_n(\hat{z}) = O_p(n^{-1}) - \frac{c_n}{n} = \frac{1}{n} \left( \underbrace{O_p(1)}_{\text{Likelihood Ratio}} - \underbrace{c_n}_{\text{Penalty}} \right). \] Since $c_n \to \infty$ while the likelihood ratio term is $O_p(1)$, the term inside the parentheses becomes negative with probability approaching 1. Although the factor $1/n$ drives the magnitude of $\Delta IC_n$ to zero, it is strictly positive and does not alter the sign. Consequently, $\Delta IC_n(\hat{z}) < 0$ asymptotically, leading to the correct selection of the restricted (unit root) model.\qed