An Automated Approach Towards Sparse Single-Equation Cointegration Modelling

An Automated Approach Towards Sparse Single-Equation Cointegration Modelling

\onehalfspacing

Introduction

In this paper we propose the Single-equation Penalized Error Correction Selector (SPECS) as a tool to perform automated modelling of a potentially large number of time series of unknown order of integration. In many economic applications, datasets will contain possibly (co)integrated time series, which has to be taken into account in the statistical analysis. Traditional approaches include modelling the full system of time series as a vector error correction model (VECM), estimated by methods such as maximum likelihood estimation Johansen1995, or transforming all variables to stationarity before performing further analysis. However, both methods have considerable drawbacks when the dimension of the dataset increases.

While the VECM approach allows for flexible modelling of potentially cointegrated series, these estimators suffer from the curse of dimensionality due to the large number of parameters to estimate. In practice they therefore quickly become difficult to interpret and computationally intractable on even moderately sized datasets. As such, to reliably apply such full-system estimators requires non-trivial a priori choices on the relevance of specific variables to keep the dimension manageable. Moreover, often one only has a single variable of interest, and estimating the parameter-heavy full system is not necessary.

On the other hand, the alternative strategy of prior transformations to stationarity is more easily compatible with single variables of interest and larger dimensions, but requires either a priori knowledge of the order of integration of individual variables, or pre-testing for unit roots, which is prone to errors in particular if the number of variables is large Smeekes2020. Additionally, this approach ignores the presence of cointegration among the variables, which may have detrimental effects on the subsequent analysis.

SPECS is a form of penalized regression designed to sparsely estimate a conditional error correction model (CECM). We demonstrate that SPECS possesses the oracle property as defined in Fan2001; in particular, SPECS simultaneously allows for consistent estimation of the non-zero coefficients and the correct recovery of sparsity patterns in the single-equation model. It therefore provides a fully data-driven way of selecting the relevant variables from a potentially large dataset of (co)integrated time series. Moreover, due to the flexible specification of the single-equation model, SPECS is able to take into account cointegration in the dataset without requiring any form of pre-testing for unit roots or testing for the cointegrating rank, and can thus be applied “as is” to any dataset containing an (unknown) mix of stationary and integrated time series. As a companion to this paper, an R package is made available that implements a fast and easy-to-interpret algorithm for SPECS estimation, and provides immediate access to the dataset used in the empirical application.\footnote{\href{https://github.com/wijler/specs}{\textcolor{blue}{https://github.com/wijler/specs}}}.

Single-equation error correction models are frequently employed in tests for cointegration Engle1987,Phillips1990b,Boswijk1994,Banerjee1998 as well as in forecasting applications EngleY1987,Chou1996, but require a weak exogeneity assumption for asymptotically efficient inference Johansen1992. Weak exogeneity entails the existence of a single cointegrating vector that only appears in the marginal equation for the variable of interest. If this assumption holds, our procedure can be interpreted as an alternative to cointegration testing in the ECM framework Boswijk1994,Palm2010. However, weak exogeneity may not be realistic in large datasets and we provide detailed illustrations of the implications of failure of this assumption and demonstrate that absent of weak exogeneity our procedure consistently estimates a linear combination of the true cointegrating vectors. While this impedes inference on the cointegrating relations, when the main aim of the model is nowcasting or forecasting, our procedure remains theoretically justifiable and provides empirical researchers with a simple and powerful tool for automated analysis of high-dimensional non-stationary datasets. In addition, for modeling a single variable of interest using a large set of potential regressors, SPECS provides a variable selection mechanism, allowing the researcher to discard variables that are irrelevant for this particular analysis. Our simulation results demonstrate strong selective capabilities in both low and high dimensions. Furthermore, a simulated nowcasting application highlights the importance of incorporating cointegration in the data as our proposed estimators obtain higher nowcast accuracies in comparison to a penalized autoregressive distributed lag (ADL) model. This finding is confirmed in an empirical application, where SPECS is employed to nowcast Dutch unemployment rates with the use of a dataset containing Google Trends series.

The use of penalized regression in time series analysis has gained in popularity, with a wide range of variants showing promising performance in applications Smeekes2018b. Recent literature has also seen the development of methods for analyzing high-dimensional (co)integrated time series.

Kock2016 proposes the adaptive lasso to estimate an augmented Dickey-Fuller regression. While this univariate model is inherently different from ours, it provides an insightful demonstration of how the lasso may be used as an alternative to testing for non-stationarity, paralleling our suggestion to consider SPECS as an alternative for cointegration testing under the assumption of weak exogeneity.

For VECM systems, Wilms2016 propose a penalized maximum likelihood approach, with shrinkage performed on the cointegrating vectors, the coefficients regulating the short-run dynamics and the covariance matrix. While their method is shown to obtain forecast gains relative to the traditional Johansen method, no theoretical results are provided. Liao2015 provide an automated method of joint rank selection and parameter estimation with the use of an adaptive penalty and derive oracle properties in a fixed-dimensional framework. Next to this theoretical limitation on its applicability to large datasets, practical implementation is further complicated due to reliance on the eigenvalue decomposition of an asymmetric matrix, which introduces complex values into the corresponding objective function. As noted by Liang2019, this results in a non-standard harmonic function optimization problem. Liang2019 propose joint parameter estimation and rank determination by employing a penalty that makes use of the $QR$-decomposition of the long-run coefficient matrix. This method possesses oracle-like properties under a high-dimensional asymptotic regime, but it requires the availability of an initial OLS estimator, thereby preventing applications on datasets in which the number of variables exceeds, or is close to, the number of available time series observations. Additionally, estimation of the long-run and short-run dynamics is performed sequentially rather than simultaneously, necessitating a two-step procedure.

In a single-equation setting, Lee2018 derive fixed-dimensional oracle properties for the adaptive lasso applied to predictive regressions where the regressors are allowed to be of mixed orders of integration. However, as a consequence of their model formulation in which all variables enter in levels, their estimator appears to be susceptible to spurious regression when the regressors are not cointegrated.

Finally, outside the penalized regression framework, Zhang2018b propose an eigenvalue decomposition to estimate the cointegrating space in the presence of any integer and fractional order of integration of the variables. However, the estimation procedure proposed by Zhang2018b does not perform variable selection, nor does it provide explicit estimates of the transient dynamics in a VECM. Onatski2019 develop a novel inference procedure for the cointegrating rank in high dimensions. Similar to the Johansen procedure, their test is based on the squared canonical correlations, for which they derive the limit spectral distribution under joint asymptotics with the use of arguments from random matrix theory.

Our proposed method provides several contributions to this existing literature. First, our theoretical results are derived in a high-dimensional framework where the number of parameters is allowed to grow with the sample size. This requires non-standard theoretical results on bounds of the smallest eigenvalue of a matrix of (co)integrated regressors, similar to those in Zhang2018b, which are further developed in this paper. Second, unlike many of the penalized regression methods surveyed above, the practical implementation of SPECS is straightforward for large datasets, including cases where the number of parameters is larger than the time dimension. Third, our method completely removes the need for pre-testing for the order of integration or cointegrating rank, and is not sensitive to spurious regression. Fourth, to the best of our knowledge, our paper is the first to explicitly allow for the presence of deterministic components in the theory, a crucial feature for many applications.

The paper is structured as follows. In Section (ref) we discuss the data generating process. Section (ref) describes the SPECS estimator. The main theoretical results of the paper are presented in Section (ref). Section (ref) contains several simulation studies, followed by an empirical application in Section (ref). We conclude in Section (ref). The main proofs and preliminary lemmas needed for them are contained in Appendix (ref), while Appendix (ref) contains results on minimum eigenvalue bounds. Finally, Appendix (ref) contains supplementary material on proofs of preliminary lemmas and additional theorems, as well as further details on the empirical application.

A word on notation. For any an $N$-dimensional vector $\bm x$, $\left\lVert\bm x\right\rVert_p = \left(\sum_{i=1}^N x_i^p\right)^{1/p}$ denotes the $\ell_p$-norm, while for any matrix $\bm D$ with $N$ columns, $\left\lVert\bm D\right\rVert_p = \underset{\bm x \in \mathbb{R}^N}{\text{max}} \frac{\left\lVert\bm D x\right\rVert_p}{\left\lVert\bm x\right\rVert_p}$ is the corresponding induced norm and $\left\lVert\bm D\right\rVert_F$ denotes the Frobenius norm. For an index set $S \subset \{1, \ldots, N\}$, let $\bm x_{S}$ be the vector containing the elements of $\bm x$ corresponding to $S$. Similarly, for a matrix $\bm D$ with $N$ rows, $\bm D_{S}$ is the sub-matrix containing the rows of $\bm D$ indexed by $S$. The orthogonal complement of $\bm D$ is denoted by $\bm D_\perp$, such that $\bm D_\perp^\prime \bm D = \bm{0}$. When $\bm D$ is a square matrix, we denote its $N$ ordered eigenvalues by $\lambda_1(\bm D) \geq \ldots \geq \lambda_N(\bm D)$ and we use $\bm D \succ 0$ to denote that the matrix is positive definite. We use $\bm \iota_N$ to denote a vector of ones of length $N$ and $\bm I_N$ to denote the $N$-dimensional identity matrix. We use $\overset{p}{\to}$ ($\overset{d}{\to}$) to denote convergence in probability (distribution) and $\overset{d}{=}$ denotes equivalence in distribution. Finally, we frequently make use of an arbitrary positive and finite constant $K$ whose value may change throughout the paper, but is always independent of the time and cross-sectional dimensions.

The High-Dimensional Error Correction Model

In this section we first discuss the data generating process for the vector time series along with the assumptions made. Next we transform the multivariate model to a single equation describing our variable of interest.

Data Generating Process

Assume one is interested in modelling a single variable of interest, say $y_t$, based on an $N$-dimensional time series $\bm z_t = (y_t,\bm x_t^\prime)$ observed at $t=1,\ldots,T$. Let $\bm z_t$ be described by

with the stochastic component given by

where $\bm A$ and $\bm B$ are $(N \times r)$-dimensional matrices containing the adjustment rates and cointegrating vectors, respectively. The innovations $\bm \epsilon_t = (\epsilon_{1,t},\bm \epsilon_{2,t}^\prime)^\prime$ satisfy the following assumptions:

This assumption implies that $\bm \epsilon_t$ is an martingale difference sequence with at least (a bit more than) four moments existing. The eigenvalue bounds in the second part place some restrictions on the dependence among the elements of $\bm \epsilon_t$, ruling out for instance a strong common factor affecting all errors. However, a wide range of contemporaneous dependence structures, such as spatial dependence, is still allowed.

The model can be rewritten into a VECM form by substituting (ref) into (ref) to obtain

where $\tau^* = (I-\sum_{j=1}^p\bm \varPhi_j)\bm \tau$. From this representation, it can directly be observed that the presence of a constant in (ref) results in a constant within the cointegrating relationship if $\bm B^\prime \bm \mu \neq \bm{0}$. Furthermore, the linear trend in (ref) appears as a constant in the differenced series and may additionally appear as a trend within the cointegrating vector if $\bm B^\prime \bm \tau \neq \bm{0}$, the latter implying that the equilibrium error $\bm B^\prime \bm z_t$ is a trend stationary process.

The following assumption asserts that the process is (at most) I(1), and the Granger Representation Theorem Johansen1995 can be applied.

Assumption (ref) enables (ref) to be written as a vector moving average (VMA) process

where $\bm C=\bm B_\perp\left(\bm A_\perp^\prime \left(\bm I_N - \sum_{j=1}^p\bm \varPhi_j\right) \bm B_\perp\right)^{-1}\bm A_\perp^\prime$, $\bm s_t = \sum_{s=1}^t \bm \epsilon_s$, $\bm C(L)\bm \epsilon_t$ is a stationary linear process and $\bm z_0$ are initial values. Without loss of generality, we assume henceforth that $\bm z_0 = \bm{0}$.

We need a further restriction on the dependence in the VMA representation in the form of the following assumption, which ensures norm-summability of the coefficients in the Beveridge-Nelson decomposition.

Single-Equation Representation

The number of parameters to estimate in (ref) is at least $2Nr + N^2p$, such that the system quickly grows too large to accurately estimate based on traditional methods. As we assume a single variable $y_t$ is of interest, we therefore instead consider the lighter parameterized single-equation model for $y_t$. To ensure that the variables modelling the variation in $y_t$ remain exogenous, we orthogonalize the errors driving the single-equation model, say $\epsilon_{y,t}$, from the errors driving the marginal equations of the endogenous variables $\bm x_t$. This is achieved by decomposing $\epsilon_{1,t}$ into its best linear prediction based on $\bm \epsilon_{2,t}$ and the corresponding orthogonal prediction error. To this end, partition the covariance matrix of $\bm \epsilon_t$ as

such that we obtain

Define $\pi_0 = \bm \varSigma_{22}^{-1}\bm \sigma_{21}$. Then, writing out (ref) in terms of the observable time series results in the single-equation model

where $\bm \delta^\prime = \left(1,-\bm \pi_0^\prime \right)\bm A\bm B^\prime$, $\bm \pi = (\bm \pi_0^\prime,\ldots,\bm \pi_p^\prime)^\prime$ with $\bm \pi_j^\prime = (1,-\bm \pi_0^\prime)\bm \varPhi_j$ for $j=1,\ldots,p$, $\mu_0 = (1,-\bm \pi_0^\prime)\left(\bm A\bm B^\prime\bm \mu + \bm \tau^*\right)$ and $\tau_0 = (1,-\bm \pi_0^\prime)\bm \tau^*$. Note that $\bm \delta$ is a vector of length $N$, whereas $\bm \pi$ is a vector of length $M = N(p+1)-1$. Additionally, $\bm w_t=(\Delta \bm x_t^\prime,\Delta \bm z_{t-1}^\prime,\ldots,\Delta \bm z_{t-p}^\prime)^\prime$ and $\epsilon_{y,t} = (1 - \bm \pi_0^\prime)\bm \epsilon_t $. Finally, we write the single-equation model in matrix notation as

where $\bm Z_{-1} = (\bm z_0,\ldots,\bm z_{T-1})^\prime$, $\bm W = (\bm w_t,\ldots,\bm w_T)^\prime$, $\bm t = (0,\ldots,T-1)^\prime$, $\bm V=(\bm Z_{-1},\bm W)$, $\bm D = (\bm \iota_T,\bm t)$, $\bm \gamma = (\bm \delta^\prime,\bm \pi^\prime)^\prime$ and $\bm \theta = (\mu_0,\tau_0)^\prime$.

In general, the implied cointegrating vector $\bm \delta$ in the single-equation model for $y_t$ contains a linear combination of the cointegrating vectors in $\bm B$ with their weights being given by $\left(1,-\bm \pi_0^\prime \right)\bm A$. Since the marginal equations of $\bm x_t$ contain information about the cointegrating relationship, efficient estimation within the single-equation model is only attained under an assumption of weak exogeneity. Johansen1992 shows that sufficient conditions for weak exogeneity to hold are (i) normality of $\bm \epsilon_t$, (ii) $\operatorname{\text{rank}}(\bm A \bm B^\prime) = 1$, i.e. there is a single cointegrating $N$-dimensional cointegrating vector $\bm \beta$, and (iii) the vector of adjustment rates takes on the form $\bm \alpha = (\alpha_1,\bm{0}^\prime)^\prime$. However, these conditions are rather restrictive when considering high-dimensional economic datasets that are likely to possess multiple cointegrating relationships and complex covariance structures across the errors. Therefore, we opt to derive our results without assuming weak exogeneity, while acknowledging that direct interpretation of the estimated cointegrating vector will only be valid in the presence of weak exogeneity. Furthermore, we believe that whether the potential loss of asymptotic efficiency in our more parsimonious single-equation model translates to inferior performance in finite samples ultimately remains an empirical question.

As we consider sparse estimation of this single-equation model, let us briefly touch upon the required sparsity. For measuring the sparsity, we work directly in the single-equation representation.\footnote{In absence of weak exogeneity, it may not be directly obvious how we obtain a sparse single-equation model from the VECM. We therefore provide a more detailed discussion of the interpretation of sparsity absent of weak exogeneity in Section (ref). In this section we just take the single-equation model directly as starting point.} Let $S_\delta = \lbrace i \vert \delta_i \neq 0\rbrace$ denote the index set of the non-zero elements in $\bm \delta$, with its cardinality denoted by $ \left\lvertS_\delta\right\rvert$, and let $S_{\pi}$ be defined accordingly for $\bm \pi$. In addition, let $r^*$ denote the dimension of the cointegration space of $\bm z_{S_\delta,t}$, i.e. the number of independent linear stationary combinations of $\bm z_{S_\delta,t}$ (cf. Remark (ref)), and define $s_\delta = \left\lvertS_\delta\right\rvert - r^*$ and $s_\pi = \left\lvertS_\pi\right\rvert+r^*$ as the number of “effective” relevant non-stationary and stationary variables, respectively. Our estimation goal will then be to obtain estimates of $S_\delta$ and $S_{\pi}$, as well as estimate $\bm \delta_{S_\delta}$ and $\bm \pi_{S_\pi}$. To obtain consistency, we need the following assumptions on the amount of sparsity.

Parts (1) and (2) put restrictions on how fast the number of relevant parameters is allowed to grow. The “effective” number of relevant stationary variables ($s_\pi$) is allowed to grow faster than the “effective” number of integrated variables ($s_\delta$), as a result of the collinearity induced by the stochastic trends (cf. Remark (ref)). Part (3) puts an additional restriction on the number of relevant coefficients as a function of the smallest non-zero coefficient. Clearly, if all coefficients are assumed to be fixed, (3) is not binding. In fact, one can allow $\gamma_\min$ to shrink at a rate up to $T^{-1/4}$ before it becomes binding. This assumption may therefore be interpreted as determining the fastest rate at which the population coefficients are allowed to decrease, as a function of $T$, $s_\delta$ and $s_\pi$, to still ensure it can be consistently picked up by our estimation method.

Rotations and Bounds on Eigenvalues

Bounds on eigenvalues play a crucial role in establishing consistency properties of lasso-type penalized regression methods. However, due the mixed integrated nature of our data, where parts of the regressors are stationary, and other parts are only stationary after rotation, the object of our assumptions is not the sample covariance matrix directly, but instead a carefully transformed version. Under Assumptions (ref)-(ref), it is then possible to ensure eigenvalue conditions on the sample covariance matrices. Before we can state the assumption, we must therefore establish some further notation and rotations to be used later.

Let $\bm \gamma = (\bm \delta^\prime, \bm \pi^\prime)^\prime$ and $S_{\gamma}$ its active set. Without loss of generality, we partition the data matrix as $\bm V=(\bm V_{S_\gamma},\bm V_{S_\gamma^c})$, with $\bm V_{S_\gamma} = (\bm Z_{-1,S_\delta},\bm W_{S_\pi})$ representing the time series carrying non-zero coefficients in the population single-equation model, henceforth referred to as the set of relevant variables. In the presence of cointegration, it follows from (ref) that the relevant lagged levels can be written as

where $\bm B_{\perp,S_\delta}$ is an $(\left\lvertS_\delta\right\rvert \times (N-r))$-dimensional matrix containing the rows of $\bm B_\perp$ indexed by $S_\delta$ and $\bm u_{S_\delta,t} = \bm C_{S_\delta}(L)\bm \epsilon_t$. The left null space of $\bm B_{\perp,S_\delta}$, defined as $\bm B^* = \left\lbrace \bm x \in \mathbb{R}^{\left\lvertS_\delta\right\rvert} \vert \ \bm B_{\perp,S_\delta}^\prime \bm x = \bm{0}\right\rbrace$, contains the linear combinations that convert $\bm z_{S_\delta,t}$ to a stationary process. Accordingly, we also refer to this null space as the cointegrating space of $\bm z_{S_\delta,t}$. By construction, $\bm \delta_{S_\delta} \in \bm B^*$, such that this cointegrating space is non-empty whenever $\bm \delta \neq \bm{0}$. In this case, we define $\bm B_{S_\delta}$ as a $(\left\lvertS_\delta\right\rvert \times r^*)$-dimensional basis matrix of $\bm B^*$, with $r^* \leq \left\lvertS_\delta\right\rvert$ representing the dimension of the cointegrating space.\footnote{The matrix $\bm B_{S_\delta}$ is not uniquely defined. However, in most instances, including those contained in the current work, identification of the span of $\bm B_{S_\delta}$ is sufficient.}

Similarly, we define $\bm B_{S_\delta,\perp}$ as a basis matrix of the left null-space of $\bm B_{S_\delta}$, i.e. a $\left(\left\lvertS_\delta\right\rvert\times(\left\lvertS_\delta\right\rvert-r^*)\right)$-dimensional matrix of full column rank with the property that $\bm B_{S_\delta,\perp}^\prime \bm B_{S_\delta} = \bm{0}$. Then, we are able to define a $\bm Q$-transformation that decomposes the reduced system into a stationary and non-stationary contribution as

For the case $\bm \delta = \bm{0}$, we define $\bm Q = \bm I_{\left\lvertS_\pi\right\rvert}$. Post-multiplication of the data matrix by $\bm Q^\prime$ gives

which we refer to as the $\bm Q$-transformed version of $\bm V_{S_\gamma}$. The first $s_\pi = \left\lvertS_\pi\right\rvert + r^*$ columns of (ref), corresponding to $(\bm Z_{-1,S_\delta}\bm B_{S_\delta},\bm W_{S_\pi})$, contain independent stationary linear combinations of the variables that are relevant to $\Delta y_t$ in the single-equation model. The remaining $s_\delta = \left\lvertS_\delta\right\rvert-r^*$ columns, given by $\bm Z_{-1,S_\delta}\bm B_{S_\delta,\perp}$, contain all linearly independent combinations that are integrated of order one.

Finally, we construct a transformed version of the sample covariance matrix based on $\bm V_{S_\gamma}$, which plays a crucial role in the development of our theory. First, to regress out the deterministic components of the observed time series in (ref), we define the matrix $\bm M = \bm I_T - \bm D\left(\bm D^\prime\bm D\right)^{-1}\bm D^\prime$.\footnote{Note that $\bm D$ may vary depending on the deterministic specification of the model; setting $\bm D = (\bm \iota_T,\bm t)$ allows for both a non-zero constant and linear trend, while simply setting $\bm M = \bm I_T$ may be desired (although not required) when it is believed that $\bm \mu=\bm \tau=\bm{0}$.} Then, after rotating by $\bm Q$ and regressing out the deterministic components by $\bm M$, the stationary and non-stationary components are scaled via the matrix $\bm S_T = \operatorname{\text{diag}}(\sqrt{T}\bm I_{s_\pi},\frac{T}{\sqrt{s_\delta}}\bm I_{s_\delta})$. Hence, our transformed sample covariance matrix is defined as

and $\hat{\bm \varSigma}_{22} = \frac{s_\delta}{T^2}\bm B_{S_\delta,\perp}^\prime\bm Z_{-1,S_\delta}^\prime\bm M\bm Z_{-1,S_\delta}\bm B_{S_\delta,\perp}$. We can now state the eigenvalue assumptions.

The first part of Assumption (ref) applies to stationary data and is known to hold when the minimum eigenvalue of the corresponding population covariance matrix is bounded away from zero Medeiros2016. The second part, however, applies to integrated variables and requires arguments that are unique to the non-stationary setting. In particular, we note the necessity of applying a scaling by $\frac{s_\delta}{T^2}$, rather than the usual $\frac{1}{T^2}$ one may expect from the fixed-dimensional literature, cf. Remark (ref). In Appendix (ref), we show several cases under which Assumption (ref) is satisfied.

The Single-Equation Penalized Error Correction Selector

Despite the dimension reduction obtained from moving towards a single-equation representation, regularization remains a necessity in high dimensions. The single-equation model (ref) contains a total of $N(p+2) + 1$ parameters, compared to the $2N(r+1) + N^2p$ parameters in the full-system VECM in (ref), resulting in a substantial reduction in dimensionality. However, the dimension may still grow large when either: (1) the number of potentially relevant variables is large or (ii) when the number of lagged differences required to appropriately model the short-run dynamics is large. Therefore, we consider the use of $\ell_1$-regularization to enable estimation in high dimensions.

The resulting estimator, henceforth referred to as the Single-equation Penalized Error Correction Selector (SPECS), is defined as the minimizer of the following objective function:

where $M = (N+1)p-1$ refers to the number of transformed variables in $\bm w_t$, i.e. the length of $\bm \pi$. We denote the minimizers of (ref) by $\hat{\bm \gamma}$. The group penalty, regulated by $\lambda_G$, serves to promote exclusion of the lagged levels as a group when there is no cointegration present in the data. In this case, the model is effectively estimated in differences and corresponds to a conditional model derived from a vector autoregressive model specified in differences. The individual $\ell_1$-penalties, regulated by $\lambda_I$, serve to enforce sparsity in the coefficient vectors $\bm \delta$ and $\bm \pi$ respectively.

The penalty of each coefficient $\gamma_i$ is weighted by $\omega_i$ to enable simultaneous estimation and selection consistency of the coefficients. Therefore, SPECS resembles a sparse group lasso Simon2013 with adaptive weighting, applied to the conditional error correction model. The weights $\omega_i$ in (ref) are typically derived from an initial estimation procedure such as OLS (if the number of variables is small enough), ridge, or lasso. In particular, let $\hat{\bm \gamma}_I$ denote initial estimates obtained for $\bm \gamma$ using one of the aforementioned methods. The weights can then be constructed as $\omega_i = \left\lvert\hat\gamma_{I,i}\right\rvert^{-k}$ for some $k>0$. As the coefficients of the irrelevant variables tend to zero, this will “blow up” the weights for these coefficients, making them unlikely to be selected in the final estimation. On the other hand, the weights of the relevant coefficients converge to a positive constant leaving them unaffected. This wedge between the weights of relevant and irrelevant coefficients is exactly needed to achieve selection consistency. As demonstrated by Zou2006a, under such assumptions on the weights, the adaptive lasso attains simultaneous selection and estimation consistency, without the necessity for the rather stringent irrepresentable condition in Zhao2006.\footnote{In fact, as the adaptive lasso can be written as a regular lasso on a transformed design matrix, the irrepresentable condition, while still needed, operates on this transformed design matrix and becomes a weighted irrepresentable condition. This condition is then in turn implied by appropriate assumptions on the weights. In this paper we directly take this route rather than going via an irrepresentable condition. Section 7.5 of Buhlmann2011 provides details on the links between these assumptions.} To maintain generality we work with general weights without specifying how they are obtained, and therefore define appropriate assumptions directly on the weights. In Section (ref) we then return to weight construction and propose a feasible way to construct weights that are theoretically shown to satisfy our assumptions.

Part ((ref)) puts an upper bound on the rate at which the weights corresponding to the relevant variables diverge. Part ((ref)) restricts the maximum admissible growth rate of the penalty. Exceeding this rate would in an excess of shrinkage bias that impedes estimation consistency. Finally, part ((ref)) states that the weights of the irrelevant variables -- interacting with the penalty parameter $\lambda_I$ -- grow sufficiently fast in order to guarantee that irrelevant variables are removed from the model with probability converging to one. The required minimum growth rate of the penalty parameter is inversely related to the growth rate of the weights of the irrelevant variables; faster diverging weights require less penalization to identify irrelevant variables.

Two common data-driven ways to select the tuning parameters $\lambda_I$ and $\lambda_G$ are using cross-validation and information criteria. As standard $K$-fold cross-validation does not respect the time order of the data, we instead consider a time series cross-validation (TSCV) scheme as proposed by e.g. Hyndman2018 and Wilms2017, where for different values of $\bm \lambda = (\lambda_I, \lambda_G)^\prime$ the model is estimated on the first part of the sample, and its prediction for the next observation is recorded. The sample is then recursively moved forward towards the end, and the $\bm \lambda$ with the lowest mean squared prediction error is selected. We refer to Smeekes2018a for details on the implementation and a comparison with traditional $K$-fold cross-validation.

While cross-validation works well for prediction chetverikov2016, it tends to generally select fairly low penalty levels and therefore includes many variables. An alternative way to select $\bm \lambda$ is using information criteria, where we find the value of $\bm \lambda$ as

where $\hat{\bm \gamma}(\bm \lambda)$ and $\hat{\bm \theta}$ denote the minimizers of $G_T\left(\bm \gamma,\bm \theta\right)$ in (ref) for a particular value of $\bm \lambda$.\footnote{As explained in Remark (ref), $\hat{\bm \theta}$ does not depend on $\bm \lambda$.} In addition, $\widehat{df}(\bm \lambda)$ is an estimate of the degrees of freedom and $C_T$ is the criterion-specific penalty; for the latter we use the Bayesian Information Criterion Schwarz1978 with $C_T= \ln(T)$.

Zou2007 show that for the (adaptive) lasso the number of non-zero coefficients is an appropriate estimate for the degrees of freedom for model selection using information criteria. For group lasso penalties, estimating the degrees of freedom is more complicated. Yuan2006 propose a heuristic rule, but this requires the least squares estimator which is not available for large $N$. Alternative rules are provided by Breheny2009 and vaiter2012 among others, but none are theoretically valid in our setting. For this reason we propose a simple, heuristic rule where we set $\widehat{df}(\bm \lambda)$ equal to the number of non-zero coefficients. Essentially this means we ignore the strength of the group penalty on the complexity of the model as long as the group is selected, thereby overestimating $df(\bm \lambda)$. As a consequence, we will only choose non-zero values of $\lambda_G$ if they either improve the fit directly or result in setting the whole group to zero without affecting the fit too much. This is an intentional choice, consistent with our theoretical treatment of the group penalty. As discussed in Remark (ref), the group penalty is not necessary and consistency can be achieved even with $\lambda_G = 0$, and can therefore be seen as an optional add-on penalty.

Finally, we note that in practice both methods require the respective objective function to be minimized for a two-dimensional grid of values for $\bm \lambda$. By choosing the lower and upper bounds of the grid carefully, one can ensure that the selected tuning parameters satisfy the assumptions listed in the next subsection. Of course, even though this ensures the theoretical validity of the selection method, its practical performance can still vary considerably. Therefore we investigate the practical performance of BIC and TSCV in the simulations and empirical application respectively.

Theoretical Results

In this section we derive the asymptotic properties of SPECS, describe the construction of the weights and discuss implications for particular model specifications.

Asymptotic Properties

The first result that we pursue is that of selection consistency, i.e. the ability of an estimation procedure to select the correct set of relevant variables with probability converging to one. In fact, Zhao2006 define a stronger property referred to as sign consistency, which additionally requires the procedure to identify the correct signs of the non-zero coefficients with probability converging to one. In the following theorem, we derive sign consistency of SPECS.

Theorem (ref) provides an asymptotic justification for implementing SPECS as a high-dimensional variable selection device. Furthermore, selection consistency is a crucial property when one aims to obtain interpretable solutions or even utilize the estimator as an alternative to classical tests for cointegration. An example of a traditional test for cointegration is the ECM-test by Banerjee1998 which looks at the $t$-ratio of the ordinary least squares coefficient of the lagged dependent variable. Alternatively, Boswijk1994 proposes to test for the joint significance of the least squares coefficients of all lagged variables with a Wald-type test. In our case, one could interpret exclusion of the lagged levels of the dependent variable, or the lagged levels of all variables, as evidence against the presence of cointegration. However, as discussed, an assumption of weak exogeneity is necessary when the aim is a direct interpretation of the estimated cointegration vector. Notwithstanding this caveat, selection consistency offers valuable insights when viewed as a screening mechanism that excludes irrelevant variables even in the absence of weak exogeneity. Moreover, since the set of variables included is strictly smaller than the time series dimension, it is possible to apply a traditional consistent estimator to the selected set of variables Belloni2013. However, ideally SPECS would contain desirable properties that omit the need of a second estimation procedure. For this reason, we establish the simultaneous consistency of the estimated coefficients in the following theorem.

The estimation consistency derived in Theorem (ref) does not place any restrictions on the relative growth rates of $T,N,p$, because it relies solely on high-level assumptions stated in the preceding section. However, when we derive sufficient conditions for the eigenvalue assumptions in Assumption (ref) in Appendix (ref) and provide a feasible method to construct weights that satisfy Assumption (ref) in Section (ref), these restrictions do appear. We refer to Section (ref) for an explicit discussion.

As a corollary to Theorem (ref), it is possible to establish a relationship between the limit distribution of SPECS and the OLS estimator based on the subset of relevant variables.

The oracle results in Corollary (ref), combined with the sign consistency from Theorem (ref), are suggestive of a post-selection inferential procedure. In particular, one may implement a two-step estimation procedure in which SPECS is used to perform variable selection in the first step and a regular OLS regression is performed on the selected variables in the second step. Then, after strengthening part (ref) of Assumption (ref) to $\lambda_I = o\left(\frac{T^{1/2-\xi}}{\sqrt{s_\delta} + \sqrt{s_\pi}}\right)$, Corollary (ref) seems to validate the use of the regular OLS distribution for this two-step estimator, essentially ignoring the variable selection from the first stage. For example, in the case where $\left\lvertS_\gamma\right\rvert$ remains finite, one could use the standard fixed-dimensional results Boswijk1994 to perform inference. However, such a post-selection inferential procedure should be treated with caution, as it is well known that the selection step impacts the sampling properties of the estimator Leeb2005. The convergence results of many selection procedures, SPECS included, hold pointwise only, i.e. the finite-sample distributions do not converge uniformly over the parameter space to their asymptotic distribution. The practical implication is that for certain values in the parameter space, relying on the oracle properties for post-selection test statistics may provide strongly misleading results. While developing a valid post-selection inference procedure to, for example, test for cointegration is certainly of interest, the field of valid post-selection inference is, despite its rapid development, still in its infancy. None of the currently existing methods, such as those considered in Berk2013, Vandegeer2014, Lee2016 or Chernozhukov2018, can easily be adapted to - let alone validated in - our setting. Developing such a method therefore requires a full new theory which is outside the scope of the current paper.

Initial Estimates

In this section, we provide the reader with a directly implementable method to construct weights that satisfy Assumption (ref). As discussed in Section 2.2, we construct the weights as $\omega_i = \left\lvert\hat{\gamma}_{I,i}\right\rvert^{-k}$. For our initial estimator we focus here on the ridge estimator, from which we can derive results for OLS as a special case, and comment on the lasso later on in the section.

Note that the power $k$ gives one the flexibility to adjust how big the wedge between relevant and irrelevant variables is. To illustrate, assume that $\hat{\gamma}_{I,i} = \gamma_i + O_p\left(T^{-a}\right)$ for all $i$. Then, it is clear that $\omega_i = O_p(1)$ when $\gamma_i \neq 0$ and $\omega_i = O_p\left(T^{ka}\right)$ when $\gamma_i = 0$. Therefore, larger values of $k$ will increase the rate at which the weights corresponding to the irrelevant variables diverge. Based on this principle, the availability of a consistent initial estimator allows us to construct weights that satisfy the conditions in Assumption (ref). However, while the idea of adjusting divergence rates through imposing varying values of $k$ seems theoretically attractive, large values of $k$ result in substantial amplification of finite-sample estimation error. As a result, the finite-sample performance of the lasso becomes unstable for large $k$, such that in practice one may want to set the value for $k$ as low as theoretically admissible.

Regardless of the choice of $k$, the basic ingredient for good adaptive weights is a consistent initial estimator. Therefore, we derive the consistency of the ridge estimator. Recall that the ridge estimator is defined as the minimizer of the following objective function:

The properties of the ridge estimator are well-studied in the stationary setting Hastie2008. However, to the best of our knowledge, no explicit results are available in the high-dimensional non-stationary case considered here.

In order to derive consistency of the ridge estimator, we redefine the transformed sample covariance matrix from Section (ref) and the corresponding bound on its minimum eigenvalue. Let $N_\delta = N-r$, $M_\pi = M+r$ and define the new scaling and rotation matrices as $\bm S_R = \operatorname{\text{diag}}\left(\sqrt{T}\bm I_{M_\pi},\frac{T}{\sqrt{N_\delta}}\bm I_{N_\delta}\right)$ and

respectively. The new transformed covariance matrix, based on the full dataset, is given by

with $\hat{\bm \varSigma}_{R,11} = \frac{1}{T}

,$ and $\hat{\bm \varSigma}_{R,22} = \frac{N_\delta}{T^2}\bm B_\perp^\prime\bm Z_{-1}^\prime\bm M\bm Z_{-1}\bm B_\perp$. Then, we extend the minimum eigenvalue bound in Assumption (ref) to (ref) as follows.

We now derive the convergence rate of the ridge estimator under a further restriction on the growth rates of $N,M$. The consistency of the ridge estimator is given in the following theorem.

Similar to Remark (ref), it follows from Theorem (ref) that $\left\lVert\hat{\bm \gamma}_R - \bm \gamma\right\rVert_2 = O_p\left(\frac{N_\delta + \sqrt{M_\pi}}{\sqrt{T}}\right)$. Based on the assumption that $\frac{N_\delta}{T^{1/4}} \to 0$ and $\frac{M_\pi}{\sqrt{T}} \to 0$ in Theorem (ref), it follows directly that $\left\lVert\hat{\bm \gamma}_R - \bm \gamma\right\rVert_2 = o_p(1)$, and therefore ridge can be used to construct weights that satisfy our Assumption (ref). The exact values of $k$ that are needed theoretically vary depending on the number of (total and relevant) variables in the dataset; we return to this issue in Section (ref).

The attentive reader may note that the admissible growth rates of $N_\delta,M_\pi$ in Theorem (ref) are the same as those initially assumed on the subsets of relevant variables, i.e. $s_\delta,s_\pi$, in Theorem (ref). The restriction imposed on the number of stochastic trends, $\frac{N_\delta}{T^{1/4}} \to 0$, corresponds closely to that of Corollary 2.1 in Liang2019, who consider (co)integrated processes as well and roughly require that $\frac{N}{T^{1/4 - \nu}} \to 0$ for some $\nu > 0$. The growth rate of the total number of (implied) stationary variables is restricted to $\frac{M_\pi}{\sqrt{T}} \to 0$. While this may seem limited in comparison to the admissible (near) exponential growth in the stationary setting with i.i.d. Gaussian errors KockCallot2015, we stress that our time series framework is more general, allowing not only for integrated processes, but also substantial dependence in the stationary component. Regarding the latter, our assumptions closely match those in the second row of Table 6 of Medeiros2016 with $\zeta = 1$, where our allowed growth rates are only slightly slower.

Ideally, we would like to allow for faster rates of divergence for the set of the irrelevant variables. A prospective strategy to attain this, would be to implement the lasso as an initial estimator, the consistency of which may be derived with the use of a compatibility condition Buhlmann2011. While desirable, deriving the validity of an appropriate compatibility condition is a considerable task. In addition to the difficulty of showing the theoretical validity of a compatibility condition in the non-stationary setting considered here, the use of a compatibility condition is further complicated by the fact that the stochastic trends asymptotically dominate the variation. More specifically, in order to attain a non-singular limit matrix, a rotation similar to $\bm Q$ is required that separates the stationary and non-stationary components in the full dataset. The standard compatibility condition would have to be adjusted in a non-trivial manner to account for such a rotation. Consequently, we leave the development of a suitable compatibility condition to future research, and instead focus on the ridge estimator under the more stringent growth rates on the number of variables. In the simulations we explore settings beyond these restrictive assumptions, and our adaptive weights continue to function in this case as well. We therefore conjecture that the suitability of the ridge estimator can be extended to a more general setting.

Implications for Particular Model Specifications

To fully appreciate the theoretical results in the preceding section, a detailed understanding of the generality provided by the set of imposed assumptions is helpful. For example, as the results are derived without requiring weak exogeneity, our set of assumptions allows for the presence of stationary variables in the data. However, in the absence of weak exogeneity, model interpretation becomes non-standard and the notion of sparsity carries non-trivial annotations. Therefore, in this section we elaborate on several relevant model specifications to demonstrate the flexibility of the single-equation model and highlight the practical implications of variable selection in such a general framework.

Sparsity and Weak Exogeneity

The benefit of $\ell_1$-regularized estimation stems from its ability to identify sparse parameter structures. However, the concept of sparsity in the conditional models here considered merits additional clarification, as the potential absence of weak exogeneity obscures standard interpretability. Accordingly, in this section we comment on the interplay between weak exogeneity and sparsity and provide several illustrative examples of sparse DGPs. For simplicity of illustration, we assume in this and the following section that $\bm \mu=\bm \tau=\bm{0}$.

In Section (ref) we argue that the coefficients regulating the long-run dynamics in the conditional model are generally derived from linear combinations of the cointegrating vectors in the VECM representation (ref). By decomposing the matrix with adjustment rates as $\bm A = ( \bm \alpha_1, \bm A_2^\prime)^\prime$, we obtain the explicit construction $\bm \delta = \bm B(\bm \alpha_1 - \bm A_2^\prime\bm \varSigma_{\epsilon,22}^{-1}\bm \sigma_{\epsilon,21})$. Hence, it follows that $\delta_i=0$ if the sparsity condition $\bm \beta_i^\prime \left(\bm \alpha_1 - \bm A_2^\prime\bm \varSigma_{\epsilon,22}^{-1}\bm \sigma_{\epsilon,21}\right) = 0$ is satisfied, where $\bm \beta_i$ is the $i$-th row of $\bm B$. While this condition may hold in a variety of non-trivial ways, specific cases of interest that lead to sparsity in $\bm \delta$ can be derived. For example, an integrated variable $x_{i,t}$ that does not cointegrate with any of the variables in the system ($\bm \beta_i = \bm{0}$), will carry a zero coefficient in the derived single-equation long-run equilibrium.

As a more general example, assume that the researcher observes the $N$-dimensional time series $\bm z_t = (\bm z_{1,t}^\prime,\bm z_{2,t}^\prime)^\prime = (y_t,\bm x^\prime_t)^\prime$, from time $t=1,\ldots,T$, where $\bm z_{1,t} = (y_t,\bm x_{1,t}^\prime)^\prime$ is an $N_1$-dimensional time series and $\bm z_{2,t}$ is an $N_2$-dimensional time series. Moreover,

In addition, assume that $\bm \varSigma_\epsilon = \operatorname{\mathbb{E}}\left(\bm \epsilon_t\bm \epsilon_t^\prime\right)$ satisfies Assumption (ref) and can be decomposed as

Then, the quantities appearing in the single-equation model in (ref) take on the form

The definitions in (ref) demonstrate that, under the restriction that the errors driving $\bm z_{1,t}$ and $\bm z_{2,t}$ are uncorrelated, sparsity in the single-equation model arises when (a subset of) $\bm z_{2,t}$ does not Granger-Cause $\bm z_{1,t}$. For example, in the extreme case where $\bm \varPi_{12} = \bm{0}$ and $\bm \varPhi_{12} = \bm{0}$, we have $\bm \delta_2 = \bm{0}$ and $\bm \pi_{j,2} = 0$, respectively. Consequently, then the single-equation model reads as

As an interesting special case, consider the decomposition in (ref) in which $z_{2,t}=\epsilon_{2,t}$ is scalar-valued with $\operatorname{\mathbb{E}}(\epsilon_{2,t}\bm \epsilon_{1,t}) = \bm{0}$. Then, it is straightforward to see that $\bm \pi_{12}=\bm \pi_{21}=\bm{0}$, $\pi_{22} = -1$ and, consequently, $\delta_N = 0$. This finding highlights that stationary variables result in sparsity in $\bm \delta$ only when they are fully exogenous, as said variables may enter the implied cointegrating vector through their correlation structure with the other variables in the system. This further demonstrates the difficulty of direct interpretation of $\bm \delta$ without imposing additional restrictions on the DGP. From a prediction perspective, however, the model's ability to include stationary variables through their correlation structure is clearly a desirable feature.

Finally, we consider a DGP in which $\bm \varSigma_\epsilon$ follows a Toeplitz structure with $\sigma_{\epsilon,ij} = \rho^{\left\lverti-j\right\rvert}$. After partitioning $\bm \varSigma_\epsilon$ as in (ref), we can rewrite

thus showing that $\bm \pi_0 = \bm \varSigma_{\epsilon,22}^{-1}\bm \sigma_{\epsilon,21} = (\rho, 0, \ldots, 0)^\prime$.\footnote{It is straightforward to show that this property carries over to covariance matrices with a block-diagonal Toeplitz structure, with each block $\bm \varSigma_\epsilon^{(k)}$ having the form $\sigma^{(k)}_{i,j}=\rho_{(k)}^{\left\lverti-j\right\rvert}$. The number of non-zero elements in the resulting vector $\bm \pi_0$ will equal the number of blocks in the covariance matrix.} As $\bm \delta^\prime = (1,-\bm \pi_0^\prime)\bm A\bm B^\prime$, this implies that only the long-run equilibria that occur in the equations for $\Delta y_t$ or its cross-sectionally neighbouring variable will be part of the linear combination in the derived the single-equation model. Consequently, any variables in the dataset that are not contained in the equilibria occurring in these equations will induce sparsity in $\bm \delta$.

Mixed Orders of Integration

One of the most prominent benefits of SPECS is the ability to model potentially non-stationary and cointegrated data without the need to adopt a pre-testing procedure with the aim of checking, and potentially correcting, for the order of integration or to decide on the appropriate cointegrating rank of the system. The assumptions under which our theory is developed are compatible with a wide variety of DGPs, including settings where the dataset contains an arbitrary mix of $I(1)$ and $I(0)$ variables. The researcher simply transforms the dataset according to (ref) and SPECS provides consistent estimation of the parameters and identification of the correct implied sparsity pattern. The purpose of this section is to demonstrate this feature by means of some illustrative examples.

The central idea underlying the above feature is that a single-equation model can be derived from any system admitting a finite order VECM representation. In a VECM system containing variables with mixed orders of integration, however, each stationary variable adds an additional trivial cointegrating vector. Such a vector corresponds to a unit vector that equals 1 on the index of the stationary variable. For illustrative purposes, we consider the following general example. Define $\bm z_t = (\bm z_{1,t}^\prime,\bm z_{2,t}^\prime)^\prime$, where $\bm z_{1,t} \sim I(0)$ and $\bm z_{2,t} \sim I(1)$ and possibly cointegrated. Let the dimensions of $\bm z_{1,t}$ and $\bm z_{2,t}$ be $N_1$ and $N_2$ respectively. Then, $\bm z_t$ admits the representation

where $\bm \varPhi(L)$ corresponds to a $p$-dimensional matrix lag polynomial by Assumption (ref) and $\bm \epsilon_t$ satisfies the conditions in Assumption (ref). As long as the design of (ref) conforms to Assumption (ref) and (ref), our main results apply to this setting and both selection and estimation consistency is maintained. For the extreme case in which all variables are integrated of order one, but none are cointegrate, we define $\bm A=\bm B = \bm{0}$. Clearly, it follows that $\bm \delta=\bm{0}$, such that the single-equation model can be seen as a conditional model obtained from a VAR specified in differences. In the other extreme case, when the levels of all variables in the VECM are weakly stationary, decomposition (ref) would simply lead to a VECM in which $-\bm A=\bm B=\bm I_N$, thereby enabling the results in Section (ref) to carry through.\footnote{When all variables are stationary, SPECS can also be shown to consistently estimate the parameters based on the well-documented properties of the adaptive lasso in stationary time series settings, such as those considered in Medeiros2016 and Masini2019.}

Rates of Convergence

We conclude our theoretical analysis with a detailed illustration of the attainable rates of convergence in different asymptotic frameworks. The rates of convergence of $\hat{\bm \gamma}_R$ and $\hat{\bm \gamma}$, as well as the specific construction of the initial weights, are dependent on the growth rates of $N,p,r,\left\lvertS_\delta\right\rvert$ and $\left\lvertS_\pi\right\rvert$. Because of the trade-off between the admissible dimension and the rate of convergence, the choice of the desired asymptotic framework is likely dependent on the specific application. For example, typical macro-economic applications are characterized by short panel datasets which would require a framework in which the cross-sectional dimension grows as fast as theoretically admissible. On the other hand, in applications with a large number of time series observations, such as forecasting based on high-frequency data, the assumption that the number of (potentially) relevant variables grows slow relative to the available time periods seems reasonable. Therefore, to aid interpretation of our results, we provide an overview with different asymptotic frameworks and the corresponding penalty parameters, weight constructions and convergence rates of the initial estimator in Table (ref). The weights for $\delta_i$ and $\pi_j$ are constructed as $\omega_i = \left\lvert\hat{\delta}_{R,i}\right\rvert^{-k_\delta}$ and $\omega_{N+j} = \left\lvert\hat{\pi}_{R,j}\right\rvert^{-k_\pi}$.

The first row of Table (ref) corresponds to the classic fixed-dimensional case. It is reassuring that, similar to the OLS estimator, SPECS obtains $\sqrt{T}$-convergence, with the additional benefit of allowing for consistent recovery of the sparsity pattern. In fact the next three rows highlight that when $N$, $p$ or $r$ diverge, while the number of relevant variables remains fixed, SPECS maintains its $\sqrt{T}$-convergence as long as the penalty weights $k_\delta$ and $k_\pi$ are adjusted appropriately. In the fifth row, we allow the number of relevant stationary variables, i.e. $\left\lvertS_\pi\right\rvert$ to diverge as well. This setting may be preferred when the integrated time series remain persistent after being transformed to stationarity by differencing. We observe that consistency is maintained, although even sharper weights are required and the rate of convergence has reduced to $T^{3/8}$. In the sixth row we additionally allow the number of relevant non-stationary, i.e. $\left\lvertS_\delta\right\rvert$, to increase, whereas the number of cointegrating vectors remains fixed. The increased number of non-zero coefficients corresponding to non-stationary variables reduces the rate of convergence to $T^{1/4}$. Interestingly, in the last row we let the dimension of the cointegrating subspace $r$ grow at the same rate. As illustrated in Remark (ref), this setting naturally occurs when the data is modelled by a non-stationary factor model with idiosyncratic components. In this framework, the number of stochastic trends driving the subset of relevant variables, i.e. $s_\delta$, remains fixed, which positively affects the convergence rate of SPECS.

We consider the theoretical results presented in this section to be of a double nature. On the one hand, it is reassuring that consistent estimation remains feasible in growing dimensions and that suitable weights are available. On the other hand, we acknowledge that the required restrictions on the growth rate of the number of variables seem to caution against application of penalized regression in very high-dimensional settings. However, it is worth noting that the restrictions on $N$ and $p$ largely result from the use of ridge regression as an initial estimator. Indeed, the availability of a novel compatibility condition could justify the use of the lasso as an initial estimator and will allow for generalization of our theoretical results to even higher dimensional asymptotic frameworks. We consider this an interesting avenue for future research.

Simulations

In this section we analyze the selective capabilities and predictive performance of SPECS by means of simulations. We estimate the single-equation model according to the objective function (ref) with the following settings for the penalty rates:

The OLS estimator is only included when feasible according to the dimension of the model to estimate and we additionally include a penalized autoregressive distributed lag model (ADL) with all variables entering in first differences. The latter model can be interpreted as the conditional model one would obtain when ignoring cointegration in the data and specifying a VAR in differences as a model for the full system. The resulting conditional model is the same as the CECM that we consider, but with the built-in restriction $\bm \delta=\bm{0}$.

We estimate the solutions for a grid of penalty values and construct the weights from an initial ridge estimator as proposed in Section (ref). For ADL and SPECS$_1$, we consider 100 possible values for $\lambda_I$ and choose the final model based on the BIC criterion. Alternatively, for SPECS$_2$, the model selection takes place over a two-dimensional grid consisting of 100 values for $\lambda_I$ and 10 possible values for $\lambda_G$, with model selection again being based on the BIC criterion. The weights are defined by $\omega_i = \left\lvert\hat{\gamma}_{R,i}\right\rvert^{-k}$, where $k=2$ for $i \in \lbrace 1,\ldots,N\rbrace$ and $k=1$ for $i \in \lbrace N+1,\ldots,N+M \rbrace$.

We consider three different settings under which we analyze the performance of our estimators; the first setting aims to analyze the effects of dimensionality and weak exogeneity, the second setting explores the effect of the variables' orders of integration and the third setting considers the performance in non-sparse settings. Each setting is described in detail below.

Dimensionality and Weak Exogeneity

In the first part of our simulation study we focus on the effects of dimensionality and weak exogeneity on a (co)integrated dataset. Our simulation DGP takes the form

with $t=1,\ldots, T=100$, $\bm \epsilon_t \sim \mathcal{N}(0,\bm \varSigma)$ and $\sigma_{ij} = 0.8^{|i-j|}$. Furthermore, $\bm \varPhi_1$, the coefficient matrix regulating the short-run dynamics is generated as $0.4 \cdot \bm I_N$, where $N$ varies depending on the specific DGP considered. Based on this DGP, the single-equation model takes on the form

with $\bm \pi_0$ and $\bm \pi_1$ as defined in (ref). We consider a total of four different settings, corresponding to different combinations of (i) dimensionality (low/high) and (ii) weak exogeneity (present/absent). The corresponding parameter settings and implied cointegrating vector $\bm \delta$ are given in Table (ref).

We measure the selective capabilities based on three metrics. The pseudo-power of the models measures the ability to appropriately pick up the presence of cointegration in the underlying DGP. For the OLS procedure we perform the Wald test proposed by Boswijk1994. When the OLS fitting procedure is unfeasible due to the high-dimensionality, we perform the Wald test on the subset of variables included after fitting SPECS$_1$ and refer to this approach as Wald-PS (where PS stands for post-selection). Despite the caveats of oracle-based post-selection inference discussed after Corollary (ref), the inclusion of Wald-PS still offers valuable insights regarding the performance one may expect of such a procedure in light of the aforementioned limitation. SPECS is used as an alternative to this cointegration test by simply checking whether at least one of the lagged levels is included in the model. The percentage of trials in which cointegration is found is then reported as the pseudo-power.

Second, for each trial the Proportion of Correct Selection (PCS) measures the proportion of correctly selected variables, while the Proportion of Incorrect Selection (PICS) describes, as the name may suggest, the proportion of incorrectly selected variables. They are given by

The PCS and PICS are calculated for SPECS$_1$ and SPECS$_2$ and averaged over all trials.

Finally, we consider the predictive performance in a simulated nowcasting application, where we implicitly assume that the information on the latest realization of $\bm x_T$ arrives before the realization of $y_T$. These situations frequently occur in practice, see Giannone2008 and the references therein for an overview as well as the empirical application considered in Section (ref). Due to the construction of the single-equation model, in which contemporaneous values of the conditioning variables contribute to the contemporaneous variation in the dependent variable, our proposed method is particularly well-suited to this application. For any of the considered fitting procedures, the nowcast is given by $\hat{y}_T = \hat{\bm \delta}^\prime \bm z_{T-1} + \hat{\bm \pi}^\prime \Delta \bm x_T + \hat{\bm \phi}^\prime \Delta \bm z_{T-1}$, where by construction $\hat{\bm \delta}=\bm{0}$ in the ADL model. For each method we record the root mean squared nowcast error (RMSNE) relative to the OLS oracle procedure fitted on the relevant variables.

Figure (ref) visually displays the evolution of our performance metrics over a range of values for $a$, representing increasingly faster rates of adjustment towards the long-run equilibrium. The first row of plots shows near-perfect performance of SPECS over all metrics. The pseudo-size is slightly lower than the size of the Wald test when the latter is controlled at 5%, whereas the pseudo-power quickly approaches one. Following expectations, the pseudo-size for SPECS$_2$ is slightly lower as a result of the additional group penalty. Focussing on the selection of variables, we find that for faster adjustment rates, SPECS is able to exactly identify the sparsity pattern with very high frequency, as demonstrated by the PCS approaching 100% and the PICS staying near 0%. Furthermore, the MSNE obtained by our methods is close to the OLS oracle method and is substantially lower than the MSNE obtained by the ADL model for faster adjustment rates, while being almost identical absent of cointegration. The picture remains qualitatively similar when moving away from weak exogeneity while staying in a low-dimensional framework, although the gain in predictive performance over the ADL has decreased somewhat. We postulate that the ADL may benefit from a bias-variance tradeoff, given that the correctly specified single-equation model is sub-optimal in terms of efficiency absent of weak exogeneity compared to a full system estimator. Nonetheless, SPECS is clearly preferred.

The performance in the high-dimensional setting is displayed in rows 3 and 4 of Figure (ref). When the conditioning variables are weakly exogenous with respect to the parameters of interest, the selective capabilities remain strong. The pseudo-power demonstrates the attractive prospect of using our method as an alternative to cointegration testing, especially when taking into consideration that the traditional Wald test is infeasible in the current setting. In addition, the nowcasting performance remains far superior to that of the misspecified ADL. The last row depicts the performance absent of weak exogeneity. In this setting, exact identification of the implied cointegrating vector occurs less frequently, which seems to negatively impact the nowcasting performance. However, the misspecified ADL is still outperformed, despite the deterioration in the selective capabilities of our method.

Mixed Orders of Integration

We next analyze the performance of SPECS on datasets containing variables with mixed orders of integration. The aim of this section is to gain an understanding of the relative performance of SPECS when not all time series are (co)integrated and to compare the performance of SPECS to traditional approaches that rely on pre-testing. The latter goal is attained by adding an additional penalized ADL model to the comparison, namely one in which the data is first corrected for non-stationarity based on a pre-testing procedure in which an Augmented Dickey-Fuller (ADF) test is performed on the individual series. We refer to this procedure as the ADL-ADF model. Based on the general DGP (ref), we distinguish four different cases, corresponding to (i) different orders of the dependent variable ($I(0)$/$I(1)$) and (ii) different degrees of persistence in the stationary variables (low/high). The choice to include varying degrees of persistence is motivated by the conjecture that the performance of the pre-testing procedure incorporated in the ADL-ADF model may deteriorate when the degree of persistence increases, which in turn translates to a decrease in the overall performance of the procedure.

The parameter settings for the varying DGPs, displayed in Table (ref), are chosen such that they allow for a subset of stationary variables in the system. In particular, we first consider a scenario in which the dependent variable itself admits a stationary autoregressive representation in levels. In addition, based on their cross-sectional ordering, the first 15 variables after $y$ are cointegrated based on three cointegrating vectors, the next 10 variables are non-cointegrated random walks, and the last 24 variables all admit a stationary autoregressive structure in levels. The degree of persistence in the stationary variables is regulated by the diagonal matrix $\tilde{\bm B}$ in $\bm B$, with elements $b_{ii} = 1$ in the low persistence case and $b_{ii} \sim U(0,0.2)$ in the high persistence case. It can be seen from the last column in Table (ref), that in line with the stationarity of the dependent variable, the first element in $\bm \delta$ will always be equal to $-1$, whereas an additional five-dimensional cointegrating vector enters the single-equation model for positive values of $a$. For the scenario in which the dependent variable is integrated of order one, the first 15 variables (including $y$) are all cointegrated based on three cointegrating vectors, the next 10 variables are non-cointegrated random walks, whereas the last 15 variables all admit a stationary autoregressive representation. The persistence in the stationary variables is regulated similar to the previous case. Now, however, it is clear from the last column in Table (ref) that $\bm \delta \neq \bm{0}$ only if $a > 0$, such that lagged levels only enter the single-equation when $y$ is cointegrated with its neighbouring variables. We display the performance of the models in Figure (ref).

In the first row of Figure (ref), corresponding to $y \sim I(0)$ and low persistence, SPECS correctly selects the lagged dependent variable in all simulation trials, such that the pseudo-power is always 1. Interestingly, PCS also seems constant around 35%. Upon closer inspection, we find that SPECS chooses an alternative representation of the single-equation model in which the contribution of the non-trivial cointegrating vector seems to be absorbed in the lagged level of the dependent variable. While the resulting model differs from the implied oracle model, which is indeed accurately estimated by the OLS oracle procedure, the model choice seems motivated by a favourable bias-variance trade-off. In line with this conjecture, the nowcast performance of SPECS occasionally exceeds the OLS oracle procedure's where a larger number of parameters is estimated. The standard ADL nowcasts are again inferior, whereas the ADL-ADF model seems to benefit from correct identification of the stationarity of the dependent variable, which is particularly relevant given that the dependent variable itself is a main component in the optimal forecast. However, the nowcast accuracy of SPECS is almost identical to that of the ADL-ADF model, a finding that we interpret as reassuring and confirmatory of our claim that SPECS may be used without any pre-testing procedure. Moreover, the absence of strong persistence in the stationary variables idealizes the results of the ADL-ADF procedure.

In typical macroeconomic applications many time series that are considered as I(0) display much slower mean reversion and, consequently, are more difficult to correctly identify as being stationary.\footnote{For example, the ten time series in the popular Fred-MD dataset which McCracken2016 propose to be I(0), i.e. the series corresponding to a tcode of one, all display strong persistence or near unit root behaviour, with the smallest estimated AR(1) coefficient exceeding 0.86.} Accordingly, in row 2 we display the result for a DGP where the stationary variables display more persistent behaviour. The performance of SPECS remains largely unaffected, whereas the nowcasting performance of the ADL-ADF model deteriorates drastically. We stress the relevance of this result, given that the estimation of ADL models after pre-testing for non-stationarity is fairly common practice. Somewhat surprisingly, the ADL model in differences nowcasts almost as well as SPECS here. Overall, however, the nowcast accuracy of SPECS remains the highest and, equally important, most stable across all specifications.

Continuing the analysis of mixed order datasets, rows 3 and 4 of Figure (ref) display the results for DGPs where the dependent variable is generated as being integrated of order one. The pseudo-power plot clearly reflects that $\bm \delta \neq \bm{0}$ only when $a>0$. Furthermore, while SPECS performs well at removing the irrelevant variables, the relevant variables are not all selected correctly, resulting in somewhat lower values for the PCS metric. Nevertheless, the nowcast performance remains superior to that of the ADL model, especially in the presence of cointegration with fast adjustment rates.

Non-sparse Data Generating Processes

To avoid idealizing the results through a choice of DGPs that suits our estimator, this section considers the performance of the penalized regression estimators in two different non-sparse settings. First, we consider an explicitly constructed VECM that contains many small, but non-zero coefficients. Second, we consider a DGP that contains a non-stationary factor structure on which the single-equation model is likely misspecified.

The non-sparse VECM is generated according to (ref) with $\bm B = \bm I_3 \otimes \tilde{\bm \iota}$, where $\tilde{\bm \iota} = (1,-\bm \iota_4^\prime)^\prime$, and $\bm A = a\bm B$ for $a=0,-0.05,\ldots,-0.5$. Hence, $N=15$ and the total number of parameters to estimate (including a constant and linear trend) is $N(p+2)+1=46$. A major difference with Section (ref) is that we do not generate the covariance matrix of the errors as a Toeplitz-matrix, the latter being a crucial driver of sparsity in the preceding sections. Instead, we implement the procedure detailed in Chang2004, in which we generate a $(N\times N)$ matrix $\bm U$ with $u_{ij} \sim U(0,1)$ to construct the orthonormal matrix $\bm H = \bm U\left(\bm U^\prime\bm U\right)^{-1/2}$, and generate a set of $N$ eigenvalues, $\lambda_1,\ldots,\lambda_N$, where $\lambda_1 = 0.01$, $\lambda_N=1$ and $\lambda_2,\ldots,\lambda_{N-1} \sim U(0.1,1)$ to construct $\bm \varLambda = diag(\lambda_1,\ldots,\lambda_N)$. We then construct the covariance matrix as $\bm \varSigma = \bm H\bm \varLambda\bm H^\prime$. At each simulation trial, we generate a new $\bm \varSigma$ such that the results cannot be attributed to a specific random draw of the covariance matrix. Based on this construction, $\bm \pi_0$, as defined below (ref), and $\bm \delta$ are non-sparse vectors with small elements; even in the setting with the strongest cointegration, i.e. $a=-0.5$, the median magnitude of the coefficients in $\bm \delta$ across all trial is only 0.12. As before, we set $T=100$ and perform 1,000 simulation trials.

The results are displayed in Figure (ref), which contain a number of interesting results. Unsurprisingly, all estimators obtain a substantially lower (pseudo-)power in the current framework. The $\ell_1$-regularized estimators seem more sensitive to this than the traditional Wald estimator considered in Boswijk1994. In line with the weak power, we observe that the PCS for both SPECS$_1$ and SPECS$_2$ is low, with on average only 0.75 out of 15 variables being included in levels.\footnote{The PICS is zero for all $a>0$, simply because the DGP is non-sparse, and is omitted accordingly.} Appropriate inference in the current setting is a difficult task and direct application of SPECS without alteration does not seem to be a feasible strategy. The development of a uniformly valid post-selection inference procedure, such as the desparsified lasso of Vandegeer2014, may alleviate some of these issues. While we consider this an interesting avenue of research, it is outside the scope of the current paper.

While these results may seem discouraging, the results on the nowcast accuracy display a different story. The mean-squared nowcast errors, relative to the OLS oracle procedure, are almost always below one and are similar for the SPECS and penalized ADL estimators. This highlights that the signal of the long-run component is so weak, that the estimation of a misspecified model which ignores cointegration benefits from a favourable bias-variance tradeoff. Therefore, the conclusion remains that SPECS obtains superior predictive performance relative to methods that ignore cointegration when the long-run component provides a strong signal, without sacrificing performance absent of cointegration or in the presence of very weak cointegration.

The second, and final, DGP that we consider contains a non-stationary factor structure and corresponds to setting III in Palm2011. We allow for contemporaneous correlation and dynamic structures in both the error processes driving the “observable” data and the idiosyncratic component in the factor structure. The DGP is given by $\bm z_t = \bm \lambda f_t + \bm \omega_t$, where $\bm z_t$ is a $(50 \times 1)$ time series process, $f_t$ is a single scalar factor and

Furthermore,

where $\bm \epsilon_{1,t} \sim \mathcal{N}(\bm{0},\bm \varSigma)$, with $\bm \varSigma$ again generated as in Chang2004, and $\epsilon_{2,t} \sim \mathcal{N}(0,1)$.

The comparison focuses exclusively on the nowcasting performance for a setting without dynamics ($\bm A_1=\bm B_1=\bm{0}$ and $\alpha_2=\beta_2=0$) and a setting with dynamics ($\alpha_2=\beta_2=0.4$). The construction of $\bm A_1$ and $\bm B_1$ is analogous to Palm2011. We report the RMSNEs of SPECS relative to the ADL in Table (ref). Given that the single-equation model is misspecified in this setup, it is unreasonable to expect SPECS to outperform. Indeed, we observe that the RMSNEs are all very close to one and, while in most cases the ADL model performs slightly better, the difference seems negligible. Hence, the risk of using SPECS to estimate a misspecified model in the sense considered here, does not seem to be higher than the use of the alternative ADL model, whereas the relative merits of SPECS when applied to a wide range of correctly specified models are evident from the first part of the simulations.

Empirical Application

Inspired by Choi2012, we consider nowcasting Dutch unemployment with SPECS based on Google Trends data. Google Trends are time series consisting of normalized indices depicting the volume of search queries entered in Google, originating from a certain geographical area. The Dutch unemployment rates are made available by Statistics Netherlands, an autonomous administrative body focussing on the collection and publication of statistical information. These rates are published on a monthly basis with new releases being made available on the 15th of each new month. This misalignment of publication dates clearly illustrate a practically relevant scenario where improvements upon forward looking predictions of Dutch unemployment rates may be obtained by utilizing contemporaneous Google Trends series.

We collect a novel dataset containing seasonally unadjusted Dutch unemployment rates from the website of Statistics Netherlands\footnote{\href{http://statline.cbs.nl/StatWeb/publication/?VW=T&DM=SLEN&PA=80479eng&LA=EN}{\textcolor{blue}{http://statline.cbs.nl/StatWeb/publication/?VW=T&DM=SLEN&PA=80479eng&LA=EN}}} and a set of manually selected Google Trends time series containing unemployment related search queries, such as “Vacancy", “Resume" and “Unemployment Benefits". The dataset comprises of monthly observations ranging from January 2004 to December 2017. While the full dataset contains 100 unique search queries, a number of these contain zeroes for large sub-periods, indicating insufficient search volumes for those particular series. Consequently, we remove all series that are perfectly correlated over any sub-period consisting of 20% of the total sample.\footnote{The dataset and corresponding R package are available at \href{https://github.com/wijler/specs}{\textcolor{blue}{https://github.com/wijler/specs}}.}

The benchmark model we consider is an ADL model fitted to the differenced data. In detail, let $y_t$ and $\bm x_t$ be the scalar unemployment rate and the vector of Google Trends series observed at time $t$, respectively, and define $\bm z_t = (y_t,\bm x_t^\prime)^\prime$. The benchmark ADL estimator fits

However, this estimator ignores the order of integration of individual time series by differencing the whole dataset, while it is common practice to transform individual series to stationarity based on a preliminary test for unit roots. Hence, similar to Sections (ref) and (ref), we include an additional ADL model where the decision to difference is based on a preliminary ADF test and refer to this method as ADL-ADF.\footnote{We note that none of the time series were found to be integrated of order 2. The outcome of the ADF test is reported for each time series in the online Appendix (ref).} Finally, SPECS estimates

All tuning parameters are obtained by time series cross-validation and we use $k=1.1$ based on a preliminary analysis.\footnote{Comparing the nowcast accuracy for varying $k \in [0,4]$, we found the highest accuracy for $k = 1.1$.} The first nowcast is made by fitting the models on a window containing the first two-thirds of the complete sample, i.e. $t=1,\ldots,T_{c}$ with $T_c = \lceil \frac{2}{3}T \rceil$, based on which the nowcast for $\Delta y_{T_c+1}$ is produced. This procedure is repeated by rolling the window forward by one observation until the end of the sample is reached, producing a total of 54 pseudo out-of-sample nowcasts. Table (ref) reports the MSNE relative to the ADL model for $p = 1,3,6$.

The ADL-ADF estimator does not perform better than the regular ADL model for $p=1,3$, indicating that the potential for errors in pre-testing might lead to unfavourable results. SPECS performs well and is able to obtain smaller mean-squared nowcast errors than the ADL benchmark across almost all specifications, with the combination SPECS$_2$ and $p=1$ being the exception. Moreover, for SPECS$_1$ ($p=3$) and SPECS$_2$ ($p=6$), we find the differences in MSNE to be significant at the 10% level according to the Diebold-Mariano test. The overall (unreported) MSNE is lowest for the SPECS$_1$ estimator based on $p=3$ lagged differences. Given that the addition of lagged levels to the models improves the nowcast performance, the premise of cointegrating relationships between Dutch unemployment rates and Google Trends series seems likely. To further explore the presence of cointegration among our time series we group our variables in five categories; (1) Application Training, (2) General, (3) Job Search, (4) Recruitment Agencies (RA) and (5) Social Security. We narrow down our focus to the nowcasts of models with three lagged difference included, $p=3$, estimated by SPECS$_1$. In Figure (ref) we visually display the share of nowcasts in which the lagged levels of each variable are included in the estimated model. In addition, it depicts the selection stability of those variables, where a green colour indicates that a given variables is included in a given nowcast, and red vice versa. The figure also displays the actual unemployment rates compared to the nowcasted values.

Figure (ref) highlights that only few variables are consistently selected for all nowcasts, although in each category we can distinguish some variables that are included at higher frequencies. The variable whose lagged levels are always selected is “Vakantiebaan", which is a search query for a temporary job during the summer holiday. We postulate that this variable is selected by SPECS to account for seasonality in the Dutch unemployment rates. In an unreported exercise we estimate the model with the addition of a set of eleven unpenalized dummies representing different months of the year. While the variable “Vakantiebaan" is never selected, the mean squared nowcast error increases substantially. Hence, we opt to adhere to our standard model under the caveat that for at least one of the lagged levels included, seasonality effects rather than cointegration seem a more appropriate explanation for its inclusion. Other frequently included variables are queries for vacancies (“uwv.vacatures", 78%), unemployment (“werkloos", 76%) and social benefits (“ww uitkering", 72%), where the stated percentages indicate the proportion of nowcast models in which the respective variables are selected. Furthermore, the last bar represents the frequency in which the lagged level of the Dutch unemployment rate is selected, which occurs for 43 out of 54 nowcasts (80%). The frequent selection of the lagged level of unemployment rates in conjunction with the other lagged levels is indicative of the presence of cointegration among unemployment and Google Trends series. However, we do not attach any structural meaning to the found equilibria based on the difficulty of interpretation when one does not assume the presence of weak exogeneity.

In an attempt to gain insights into the temporal stability of our estimator, we visually display the selection stability in the bottom-left part of Figure (ref). Generally, we see that for the early and later period of the sample very few time series enter the model in levels, whereas for the middle part of the sample the majority of variables are selected. The exact reason for these patterns to occur is unknown and raises questions on the stability of Google trends as informative predictors of Dutch unemployment rates. Standard feasible explanations concern structural instability in the DGP, seasonality effects or data idiosyncrasies. However, there are additional peculiarities specific to the use of Google trends such as normalization, data hubris and search algorithm dynamics, all of which might result in unstable performance Lazer2014. Since the focus of this application is on the relative performance between our estimator and a common benchmark model, rather than on a structural analysis of the relation between Google Trends and unemployment rates, we leave this issue aside as it is outside the scope of the paper. Instead, we focus on the relative empirical performance of our methods, which, notwithstanding the aforementioned caveats, we deem convincingly favourable for SPECS. Finally, on the right of Figure (ref) we display the realized and predicted unemployment rates in levels and differences. Both the penalized ADL model and SPECS seem to follow the actual unemployment rates with reasonable accuracy, with the largest nowcast errors occurring in the first half of 2014. Prior to this period the unemployment rates had been steadily rising in the aftermath of the economic recession, whereas 2014 marks the start of a recovery period. Given that the models are fit on historical data, it is natural that the estimators overestimate the unemployment rate shortly after the start of the economic recovery. Perhaps not entirely coincidental, the start of the period over which the majority of lagged levels are included by SPECS coincides with this recovery period as well, thereby hinting towards structural instability in the DGP as a plausible cause for the observed selection instability.

Conclusion

In this paper, we propose the use of SPECS as an automated approach for sparse single-equation error correction modelling in high-dimensional settings. SPECS is an intuitive estimator that applies penalized regression to a conditional error-correction model. We show that SPECS possesses the oracle property and is able to consistently select the long-run and short-run dynamics in the underlying DGP. These results are derived with the aid of a novel bound on the minimum eigenvalue of the sample covariance matrix containing integrated process, which may be of independent interest. Additionally, in pursuit of suitable weights that aid in the identification of the subset of relevant variables, we derive the consistency of the ridge estimator applied to the same model and demonstrate how ridge regression may be used to construct these weights.

We document favourable finite sample performance of SPECS by means of simulations and an empirical application. The simulation exercise confirms strong selective and predictive capabilities in both low and high dimensions with convincing gains over a benchmark penalized ADL model that ignores cointegration in the dataset. Furthermore, the simulation results demonstrate that the selective capabilities of SPECS remain adequate absent of weak exogeneity and the nowcasting performance remains superior to the benchmark. Finally, we consider an empirical application in which we nowcast the Dutch unemployment rate with the use of Google Trends series. Across all three different dynamic specifications considered, SPECS attains higher nowcast accuracy, thus confirming the findings from our simulation study. As a result, we believe that our proposed estimator, which is easily implemented with readily available tools at a low computational cost, offers a valuable tool for practitioners by enabling automated model estimation on relatively large and potentially non-stationary datasets and, most importantly, allowing to take into account potential (co)integration without requiring pre-testing procedures.

Finally, we highlight several important sources through which the assumptions and asymptotic framework may be generalized further. Sharper and more direct eigenvalue bounds can be utilized to cast SPECS into an even higher-dimensional setting. Similarly, a suitable compatibility condition can be used to validate the lasso as an initial estimator, resulting in improved weights and, again, a less restrictive asymptotic framework. These topics remain subject to our continuing investigation.

\numberwithin{theorem}{section} \numberwithin{lemma}{section} \numberwithin{corollary}{section} \numberwithin{proposition}{section} \numberwithin{table}{section}