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.
95,290 characters · 19 sections · 81 citation commands
Determination of the effective cointegration rank in high-dimensional time-series predictive regressions
\affil[1]{ School of Economics, Zhejiang University} \affil[2]{ Center for Data Science, Zhejiang University} \affil[3]{ Booth School of Business, University of Chicago}
The availability of large-scale or vast time-series data in recent years brings new challenges and opportunities to time series modeling. Analysis of high-dimensional (HD) time series has emerged as one of the important and active research areas in statistics, economics, finance, and engineering, among other scientific fields. For example, returns of a large number of assets form a HD time series and play an important role in asset pricing, portfolio allocation, and risk management. Environmental studies often employ HD time series consisting of a large number of pollution indexes collected from many monitoring stations over time. In many applications, data often exhibit characteristics of unit-root nonstationarity. For instance, the series of quarterly gross domestic products, total exports, and total imports of an economy tend to contain unit roots. In theory, the vector autoregressive moving-average (VARMA) models can be used to analyze such data, but they often encounter the difficulties of cointegration testing, overparametrization, and lack of identifiability. See, for example, johansen2002small,tiao1989model,lutkepohl2006new,tsay2014multivariate, and the references therein. To overcome these difficulties, dimension reduction or structural regularization becomes a necessity, and various methods have been developed in the literature including the regularized estimation method for HD VAR models in lin2017regularized and the factor modeling by stock2005implications,bai2002determining,forni2005generalized,pena2006nonstationary,lam2011estimation,lam2012factor,gao2019structural,gao2021modeling_JASA,gao2021modeling_IJF,gao2022divide, among others. However, most of the studies mentioned above focus on stationary processes and are not applicable to unit-root nonstationary series. The only exceptions are bai2004estimating,pena2006nonstationary,gao2021modeling_IJF. On the other hand, the unit-root nonstationarity is commonly seen in many empirical applications and the complexity of the dynamical dependence in such data requires further investigation.
It is well known that cointegration is often used to account for common trends and to avoid non-invertibility induced by over-differencing unit-root time series. See engle1987co,johansen1988statistical,johansen1991estimation,tsay2014multivariate, and the references therein. In practice, the cointegration rank of a given vector time series is unknown, and many approaches have been proposed to estimate the rank; see, for example, engle1987co,johansen1988statistical,johansen1991estimation,saikkonen2000testing,aznar2002selecting. However, these methods are rarely applied to HD time series due to their poor finite-sample performance, as discussed in johansen2002small. Yet there are many real applications that involve HD time series. For example, banerjee2014forecasting emphasized the importance of testing for no cross-sectional cointegration in panel cointegration analysis, and the cross-sectional dimension of modern macroeconomic panel can easily be as large as several hundreds. Recently, there are some studies on identifying the cointegration rank of unit-root time series from a factor modeling perspective. See pena2006nonstationary for the case of fixed dimensions and bai2004estimating,zhang2019identifying,gao2021modeling_IJF for HD time series. However, the situation changes in the case of growing dimension because the estimated cointegration rank usually grows as the dimension increases and the cointegration relationships are often hard to interpret when there are many cointegrating vectors.
This paper marks a further development in estimating the cointegration rank of HDUR time series from a predictive perspective. To avoid employing a large number of cointegrating vectors given by a high-dimensional method, we estimate the {\it effective cointegration rank} in a predictive framework. Specifically, suppose our goal is to predict the future values of a HD stationary time series $\by\in \bbR^{p}$ using $\bx$ as predictors. It is well known that only the cointegrated series have potential predictive power for the stationary process $\by$. If the number of cointegrating vectors is large, the stacked variables obtained by cointegrating vectors form a HD stationary time series, and can be used as potential predictors. But not all cointegrated series have predictive power for $\by$, and we define the {\it effective cointegration rank} as the effective dimension of the stacked variables that have predictive power for $\by$. The resulting effective rank can be much smaller than the cointegration rank of $\bx$.
The proposed method consists of a two-step estimation procedure. First, we postulate that the HDUR time series follows a factor model as that specified in bai2004estimating, where the common factors capture the nonstationary common trends of all the components, and the idiosyncratic term is a stationary process. We apply the Principal Component Analysis (PCA) to estimate the common stochastic trends and their associated loading matrix, and the orthogonal complement of the loading matrix consists of the cointegrating vectors. Second, we put together all stationary series obtained by the cointegrating vectors of the first step to form a set of predictors, and perform a reduced-rank regression between the $\by_t$ series of interest and the predictors. To further explain the variability of the data, we also include some lagged variables of $\by_t$ in the regression and assume their coefficient matrices are of low-dimensional structures. We propose two procedures to estimate all the coefficient matrices depending on whether the autoregressive (AR) matrices are sparse or of low-rank. When the AR coefficient matrices are sparse, we apply the nuclear norm penalty to the regression coefficient matrix of the stationary predictors obtained from the first step, and the LASSO penalty to the coefficients of the lagged variables. When both the AR matrices and the coefficient matrix of the predictors are of low-rank, we propose an integrative reduced-rank approach to estimate all unknown parameters. Two iterative, alternating procedures are proposed to estimate all unknown coefficients under the two aforementioned scenarios. Theoretical properties of the estimators are established under the assumption that the dimensions $p$ and $N$ and the sample size $T \to \infty$. Both simulated and real examples are used to illustrate the proposed procedure. The empirical application suggests that the 13 macroeconomic variables from welch2008comprehensive provide satisfactory performance as predictors in forecasting the returns of 79 stocks in the S&P 500 index.
The idea of using predictive regression to estimate the cointegrating vector can be found in, for example, koo2020high. However, the method of koo2020high only identifies one cointegrating vector in predicting another univariate time series, whereas the proposed method not only recovers the total cointegration rank, but also identifies the effective cointegration rank in predicting a large panel of time series. In addition, the proposed estimation method is different from theirs as we use PCA, reduced-rank, and LASSO techniques to achieve our goals while they focus mainly on the use of LASSO regularization. Note also that our framework is established for data with time series dependence structure and we use a combination of reduced-rank and sparsity techniques in the estimation procedure, which is different from most of the methods discussed in reinsel2022multivariate that focus on the reduced-rank techniques for {\it i.i.d.} observations. The only exception is the work of lin2017regularized in studying regularized estimation for multi-block stationary VAR models using the tools and techniques developed in negahban2011estimation and agarwal2012noisy. Furthermore, none of the work mentioned above deals with HDUR time series data.
This paper makes multiple contributions. First, the cointegration problem has been a central issue in modeling HDUR time series, but the lack of clear interpretations of a large number of cointegrating vectors renders the existing HD methods less appealing. We define the effective cointegration rank to select the most significant cointegration relationships from a predictive point of view using reduced-rank method. This method often produces a small number of significant cointegrating vectors which are easier to interpret in general. Second, our predictive regression model consists of both nonstationary and stationary variables as predictors and has a wide range of applications including the prediction of stock returns using macroeconomic series in finance and the prediction of PM$_{2.5}$ values using other air pollution and meteorological indexes in environmental studies. Third, the proposed approach combines the advantages of using two regularization methods, reduced-rank and LASSO, to reduce the dimension of a large system, and the asymptotic results derived suggest that properties of both methods continue to hold when they are used simultaneously in a regression model with serially dependent data. This is a theoretical contribution.
This paper is organized as follows. We introduce the proposed model, estimation methodology, and the modeling procedure in Section (ref). Section (ref) is devoted to theoretical properties of the proposed model and its associated estimates, and Section (ref) presents some simulation results to demonstrate the performance of the proposed method in finite samples. In Section (ref), we apply the proposed method to the prediction of stock returns using some commonly used macroeconomic predictors. Section (ref) provides some discussions and concluding remarks. All technical proofs of the theorems are relegated to an online supplement.
Let $\by_t = (y_{1t},y_{2t},\cdots,y_{pt})'$ be an observable $p$-dimensional stationary time series, and $\bx_t = (x_{1t},x_{2t},\cdots,x_{Nt})'$ an observable $N$-dimensional $I(1)$ process. We consider the following predictive regression model:
where $\mat{W}$ is a $p\times N$ coefficient matrix associated with the $I(1)$ process $\bx_t$, and $\mat{\Phi}_i$ is the $p\times p$ coefficient matrix of $\by_{t-i}$, for $1\leq i\leq d$, and $\vect{e}_t\sim$ WN$(0,\bSigma_e)$ is a white noise error term with mean zero and a nonsingular covariance $\bSigma_e$. Our goal is to estimate $\bW$ and $\bPhi_i$ based on a given sample, and to forecast future values of $\by_t$. For simplicity, all variables are set to zero if the time index is not positive. Also, Model ((ref)) can be extended to multi-step ahead predictions for $\by_{t+h}$ with $h > 0$.
In Model ((ref)), $\bx_t$ is nonstationary but all other variables are stationary so that it only makes sense if some variables in $\bx_t$ are cointegrated, otherwise, $\bW$ would essentially be a zero matrix because the correlation between a stationary process and a unit-root nonstationary one is zero in general. If we blindly apply the Least Squares (LS) method to estimate the model, the number of parameters to be estimated is large, and the resulting estimator $\wh \bW$ would be hard to interpret as we do not know whether all or only a few rows in $\wh\bW$ are the estimated cointegrating vectors. In theory, if all the cointegrating vectors of $\bx_t$ are known, then the resulting linear combinations of the $I(1)$ variables are stationary and can be useful predictors in Model ((ref)). However, not all cointegration relationships are helpful in predicting $\by_{t+h}$ in general, especially when the dimension $N$ of $\bx_t$ is large.
In view of the above discussion, we modify Model ((ref)) as follows. First, similarly to the setting in bai2004estimating, we assume that $\bx_t$ admits a latent factor structure:
where $\bff_t=(f_{1t},f_{2t},\ldots,f_{rt})'$ is an $r$-dimensional factor process that constitutes the common stochastic trends of $\bx_t$, that is,
where $\vect{u}_t$ is an $r$-dimensional zero-mean stationary process that drives $\bff_t$. The idiosyncratic term $\bve_t$ in ((ref)) is assumed to be a stationary process independent of the common factors $\bff_t$. Therefore, the cointegration rank of $\bx_t$ is $N-r$. For ease in model identification, we assume that $\bB$ is an orthonormal matrix such that $\bB'\bB=\bI_r$; see also bai2002determining,fan2013large for details.
Let $\bB_c\in\bbR^{N\times(N-r)}$ be an orthogonal complement matrix of $\bB$ such that $\bB_c'\bB_c=\bI_{N-r}$ and $\bB_c'\bB={\bf 0}$. It follows from Model ((ref)) that the columns of $\bB_c$ can be treated as a set of cointegrating vectors of $\bx_t$ because $\bB_c'\bx_t=\bB_c'\bve_t$ is stationary. Letting
we define $\bW=\bA\bB_c'$ and rewrite Model ((ref)) as follows:
where $\bz_t$ is now a stationary process defined in ((ref)). Similarly to the identifiability issue in factor models, $\bA$ and $\bB_c$ are not uniquely defined. Nonetheless, the product, $\bW=\bA\bB_c'$, is uniquely defined. Therefore, we split Model ((ref)) into ((ref))--((ref)), and our goal is to estimate the factor loading matrix $\bB$ or equivalently the cointegrating vector matrix $\bB_c$, the coefficient matrices $\bA$ and $\bPhi_i$, for $1\leq i\leq d$.
Although $\bB$, $\bB_c$ and $\bA$ are not uniquely defined due to the identification issue, the linear spaces spanned by the columns of $\bB$ and $\bB_c$, denoted as $\mathcal{M}(\bB)$ and $\mathcal{M}(\bB_c)$ respectively, are uniquely defined. For any specific choice of $\bB_c$, $\bA$ can also be uniquely determined. Therefore, when we mention the estimation or consistency of the loading matrix $\bB$ or $\bB_c$ in the sequel, we always refer to their column spaces to avoid any confusion. The estimation of $\bA$ is also based on a given and fixed $\bB_c$ so that the procedure is valid.
We consider two approaches to estimating the effective cointegration rank, or equivalently, the reduced-rank of the coefficient matrix $\bA$, and the AR coefficients $\bPhi_i$'s for high-dimensional cases under different assumptions. The first approach is based on imposing a reduced-rank structure on the matrix $\bA$ and some sparsity assumptions on the AR coefficient matrices. The second approach requires that all predictors, including the lagged variables, have their own low-rank coefficient matrices.
In this section, we introduce a Reduced-Rank and Sparse Regression approach (RRSRA) to estimating the coefficient matrices $\bA$ and $\bPhi_i$ for observed data $\{\bx_1,...,\bx_T\}$ and $\{\by_1,...,\by_T\}$. Note that the dimensions of $\bA\in \bbR^{p\times (N-r)}$ and $\bPhi_i\in\bbR^{p\times p}$ can be very large under the assumption that the number of common stochastic trends $r$ is finite as $p,N\rightarrow\infty$. Even if $\{\bz_1,...,\bz_T\}$ were given, the traditional LS method would lead to overfitting because there are many parameters to estimate. Therefore, some structure regularization must be imposed on the coefficient matrices. For simplicity, we assume the matrix $\bA$ is singular and has a reduced-rank form with $r_{\bA}=\text{rank}(\bA)\ll\min(p,N-r)$, and the AR coefficient matrices $\bPhi_i$'s are sparse in the sense that only a small number of elements in each matrix are nonzero, for $1\leq i\leq d$.
Assume that the number of common stochastic trends $r$ in Model ((ref)) and the order $d\geq 1$ in Model ((ref)) are known. Their selections will be discussed below. Note that $\bz_t$ is unobservable in Model (ref) and needs to be estimated from the data $\bx_t$. We briefly introduce the proposed two-step estimation procedure. First, similarly to that in bai2004estimating, we estimate the factor loading matrix $\bB$ by solving the following optimization problem:
where $\mat{X}=[\bx_1,\bx_2,\ldots,\bx_T]$ and $\mat{F}=[\bff_1,\bff_2,\ldots,\bff_T]$ are the stacked matrices across the time horizon. It is not hard to show that the optimization method in ((ref)) is equivalent to Principal Component estimation, and the columns of $\wh\bB$ are just the $r$ standardized eigenvectors of $\bX\bX'$ associated with the $r$ largest eigenvalues. Therefore, we choose $\wh\bB_c$ such that its columns are the $N-r$ standardized eigenvectors associated with the $N-r$ smallest eigenvalues of $\bX\bX'$. Then, we define $\wh\bz_t=\wh\bB_c'\bx_t$, which serves as a proxy of $\bz_t$ and will be used as predictors in the second step of estimation.
Next, we introduce a method to estimate the coefficient matrices $\bA$ and $\bPhi_i$, for $1\leq i\leq d$. To begin, define $\bPhi=[\bPhi_1,...,\bPhi_d]\in \bbR^{p\times dp}$ and $\bP_{t-1}=(\by_{t-1}',...,\by_{t-d}')'$. For any given penalty parameters $\lambda_{\bA}>0$ and $\lambda_{\bPhi}>0$, we solve the following optimization problem:
where the data are set to ${\bf 0}$ if the subscript $t\leq 0$. For reduced-rank regression, we refer the readers to the new monograph by reinsel2022multivariate. In particular, its Chapters 9 to 12 discuss some recent developments in reduced-rank regressions under high-dimensional settings, including the use of nuclear-norm penalty in ((ref)). Similar ideas can also be found in negahban2011estimation,chen2013reduced, among others. However, most of the methods considered in the aforementioned literature only deal with i.i.d. data, while we consider serially dependent data in this paper both theoretically and empirically.
It is generally not easy to obtain the true global solutions to the optimization problem in (ref) because the objective function in the bracket of ((ref)) involves different types of penalties. Therefore, we formulate an iterative procedure to obtain an approximate set of numerical solutions to ((ref)) in Algorithm (ref). Specifically, for a fixed $\bA$, we can estimate $\bPhi$ via a standard LASSO procedure, and there are several methods and software packages available to obtain sparse solutions. See, for example, hastie2015statistical. When $\bPhi$ is fixed, the estimation of $\bA$ is an instance of a {\it semidefinite program}. See vandenberghe1996semidefinite,ji2009accelerated. Since the objective function is convex, it is also biconvex in both sets of parameters. If the estimates in all iterations lie within a small ball around the true parameters, the convergence of the estimates to a stationary point is guaranteed. Because the function is convex, the estimates also achieve a global minimum. See, for example, tseng2001convergence and burai2013necessary. The theoretical results in Section (ref) below are developed for the optimal solutions $\wh\bA$ and $\wh\bPhi$. The simulation results in Section (ref) suggest that the initial values in Algorithm (ref) have little impact on the asymptotic behavior of the estimates.
Next we turn to the interpretation of the low-rank structure of the matrix $\bA$ in Model ((ref)). From negahban2011estimation, we see that the estimation of the rank of $\bA$ is equivalent to an optimal selection of the penalty parameter $\lambda_{\bA}$. A similar argument applies to the sparsity of $\bPhi$ and the choices of $\lambda_{\bPhi}$. See also, hastie2015statistical. Suppose the true rank of $\bA$ is equal to $k_0\ll\min\{p,N-r\}$, we may decompose $\bA$ as $\bA=\bC\bR'$ with $\bC\in \bbR^{p\times k_0}$ and $\bR\in \bbR^{(N-r)\times k_0}$. Therefore, $\bR'\bz_t$ is a $k_0$-dimensional stationary random vector and has some predictive power for the future values of $\by_t$. Note that $\bR'\bz_t=\bR'\bB_c'\bx_t$, implying that $\bR'\bB_c'$ is the reduced-rank matrix consisting of $k_0$ significant cointegrating vectors that play important roles in predicting the future values of $\by_t$. In other words, the cointegration rank $N-r$ can be reduced to a smaller number $k_0$ which is useful in prediction and is also easier to interpret. We call $k_0$ the {\em effective cointegration rank} of such a prediction application.
In this section, we introduce an integrative reduced-rank approach (IRRA) to estimating all the coefficient matrices of Model ((ref)). The approach is similar to the setting in Chapter 10 of reinsel2022multivariate for i.i.d. observations, but we focus on Model ((ref)) with time-series dependence. Specifically, under Models ((ref))--((ref)), the IRRA assumes that each set of predictors has its own low-rank coefficient matrix, that is, in addition to the assumption in Section 2.2.1 that $r_{\bA}\ll \min(p,N-r)$, we also assume that $0\leq r_i=\text{rank}(\bPhi_i)\ll p$, for $1\leq i\leq d$, when both $p$ and $N$ are large. This approach bridges the reduced-rank and the sparse models in the sense that the coefficient matrix $\bPhi_i$ is fully sparse with all entries being zero if $r_i=0$. On the other hand, the groupwise low-rank structure in IRRA is more flexible and different from a globally low-rank structure for $\bPhi$ defined in Section 2.2.1. The low-rankness of $\bPhi_i$'s does not necessarily imply that $\bPhi$ is of low rank, while a low-rank matrix $\bPhi$ implies that each $\bPhi_i$ is of low-rank, which cannot exceed that of $\bPhi$. Under this assumption, we consider the following convex optimization problem:
where $\lambda_{\bA}$ and $\lambda_i$ are the penalty parameters associated with $\bA$ and $\bPhi_i$, respectively. Similarly to the setting of reinsel2022multivariate, we may rewrite $\lambda_i$ as $\lambda_i=\lambda_{\bPhi}w_i$ for a global penalty $\lambda_{\bPhi}$ and some prescribed constant $w_i$, for $1\leq i\leq d$. It is clear that $\lambda_i$ is a tuning parameter controlling the amount of regularization applied to $\bPhi_i$. If $w_i=1$, and hence $\lambda_1=...=\lambda_d$, all penalty parameters of $\bPhi_i$ are the same. A simple choice is to take \[w_i=\sigma_1(\bY)\{\sqrt{p}+\sqrt{rank(\bY)}\}/T,\,\,i=1,...,d,\] so that we only have a single parameter $\lambda_{\bPhi}$ to control the regularization of the coefficient $\bPhi_i$, for $1\leq i\leq d$.
Because the objective function in ((ref)) is convex, there are several feasible algorithms available to solve the optimization problem therein. For example, following the recipe in boyd2011distributed, li2019integrative proposed an Alternating Direction Method of Multipliers (ADMM) algorithm to fit a model similar to that in ((ref)) with reduced-rank structures. However, the ADMM algorithm is relatively more involved as it alternates between a primal step and a dual step. In this paper, we propose an easy-to-implement iterative procedure to estimate all the coefficient matrices with reduced-rank structures, which is similar to the block coordinate descent method in tseng2001convergence. The detailed procedure is outlined in Algorithm (ref) below. Since the objective function is convex, by the argument in tseng2001convergence, the convergence of the estimators via Algorithm (ref) to a stationary point is guaranteed. On the other hand, from boyd2004convex, we know that the conjugate of a conjugate function of a convex one is itself, by Theorem 2 of burai2013necessary, the stationary point obtained by Algorithm (ref) is a global minimum. Simulation results in Section 4.2 suggest that the estimators obtained by Algorithm (ref) are comparable to those obtained by the ADMM method, while the former is much easier to implement than the latter in practice. Similarly to the argument used at the end of Section 2.2.1, the cointegration rank has been reduced to a much smaller and effective one due to the reduced-rank structure of $\bA$. We omit the details to save space.
The estimation of $\wh\bB$ and its orthogonal complement $\wh\bB_c$ in the prior sections is based on a given $r$, which is unknown in practice. There are several methods available in the literature to determine the number of unit-root factors in Equation ((ref)). See, for example, the information criterion in bai2004estimating, the Canonical Correlation Analysis (CCA) method in pena2006nonstationary, the autocorrelation-based method in zhang2019identifying and its modified version in gao2021modeling_IJF, among others.
In this paper, we adopt the auto-correlation based method of gao2021modeling_IJF. Specifically, let $\wh\bXi=(\wh\bxi_1,...,\wh\bxi_N):=[\wh\bB,\wh\bB_c]$ be the matrix containing all the eigenvectors of $\bX\bX'$ and $\wh f_{j,t}=\wh\bxi_j'\bx_t$ be the $j$-th principal component, for $1\leq j\leq N$. For some prescribed integer $\bar{k}>0$, define
where $\wh\rho_j(k)$ is the lag-$k$ sample autocorrelation function (ACF) of the principal component $\wh f_{j,t}$, for $1\leq j\leq N$. If $\wh f_{j,t}$ is stationary, then under some mild conditions, $\wh\rho_j(k)$ decays to zero exponentially as $k$ increases, and $\lim_{\bar{k}\rightarrow\infty} S_{j}(\bar{k})<\infty$ as $T\rightarrow\infty$. If $\wh f_{j,t}$ is unit-root nonstationary, then $\wh\rho_j(k)\rightarrow 1$, and $\lim_{\bar{k}\rightarrow\infty} S_{j}(\bar{k})=\infty$ as $T\rightarrow\infty$. Therefore, we start with $j=1$. If the average of the absolute sample ACFs $S_{j}(\bar{k})/\bar{k}\geq \delta_0$ for some constant $0< \delta_0<1$, then $\wh f_{j,t}$ has a unit root and we increase $j$ by $1$ to repeat the detecting process. This detecting process is continued until $S_{j}(\bar{k})/\bar{k}< \delta_0$ or $j=N$. If $S_{j}(\bar{k})/\bar{k}\geq \delta_0$ for all $j$, then $\wh r=N$; otherwise, we denote $\wh r=j-1$.
In this section, we briefly introduce a way to choose the tuning parameters $\lambda_{\bA}$ and $\lambda_{\bPhi}$, and the order $d$ in ((ref)). We only consider the procedure introduced in Section 2.2.1 since the one in Section 2.2.2 is similar. We first fix the order $d$ and consider the subsamples $\{\by_1,...,\by_{T_1+j}\}$ and $\{\bx_{1},...,\bx_{T_1+j-1}\}$, for $0\leq j\leq T-T_1-1$ and $T_1<T$. We then adopt a rolling-window-based method to select $\lambda_{\bA}$ and $\lambda_{\bPhi}$ from a forecasting perspective. Specifically, we prescribe two candidate intervals $[a_1,a_2]$ and $[b_1,b_2]$ with $a_2>a_1>0$ and $b_2>b_1>0$, and choose $(\lambda_{\bA},\lambda_{\bPhi})$ from $[a_1,a_2]\times [b_1,b_2]$ via a grid-search approach. For any pair $(\lambda_{\bA},\lambda_{\bPhi})\in [a_1,a_2]\times[b_1,b_2]$ and each $0\leq j\leq T-T_1-1$, we first estimate the loading matrix and obtain the stationary process $\{\wh\bz_1,...,\wh\bz_{T_1+j-1}\}$ based on the sample $\{\bx_1,...,\bx_{T_1+j-1}\}$, and apply the iterative procedure in Algorithm (ref) to obtain the estimators for all the coefficients based on the subsample $\{\by_1,...,\by_{T_1+j}\}$. We can then obtain the predicted value $\wh\by_{T_1+j+1}$ for $\by_{T_1+j+1}$. We repeat the above procedure for $0\leq j\leq T-T_1-1$ and obtain all the forecasts $\{\wh\by_{T_1+j+1},...,\wh\by_{T}\}$. Define the average of forecast errors as
Note that the forecast errors defined in (ref) also depend on the value of $d$, which itself is unknown in practice. We may prescribe an integer $\bar{d}>0$ and search the optimal one over $0\leq \wh d\leq \bar{d}$ such that the forecast error is minimized. Consequently, the optimal tuning parameters are chosen as
In practice, for simplicity, $\bar{d}$ is often chosen as a small integer provided that the series under study is not seasonal. This choice can also be justified theoretically, because the marginal model of a $p$-dimensional VAR($d$) process is ARMA($pd, p(d-1)$) the order of which can be sufficiently high when $p$ is large; see, for instance, Chapter 2 of tsay2014multivariate. In this paper, we choose $\bar{d}=3$ and the proposed model and procedure work sufficiently well in the real data analysis.
In this section, we investigate some theoretical properties of the coefficient estimates $\wh{\bB}$, $\wh{\mat{A}}$, and $\wh{\mat{\Phi}}$ under the condition that $p,N,T\rightarrow\infty$. We start with some assumptions and postpone proofs of all theorems to an online supplement.
With the identification condition $\bB'\bB=\bI_r$, the processes $\bff_t$ and $\bu_t$ have an additional strength of $\sqrt{N}$. For the stationary process $\bu_t$ in ((ref)), define a normalized process \[ \vect{S}_{T}^r (\vect{t}) = (S_{T}^1 (t_1),\ldots,S_{T}^r (t_r))' = \left(\frac{1}{\sqrt{NT}}\sum_{s=1}^{[T t_1]}u_{1s}, \ldots,\frac{1}{\sqrt{NT}}\sum_{s=1}^{[T t_r]}u_{rs}\right)', \] where $\vect{t}=(t_1,t_2,\ldots,t_r)'$ is a constant vector with $0\leq t_1\leq\cdots\leq t_r \leq 1$.
Assumption (ref) is standard for dependent random processes. For a theoretical justification of the mixing conditions for VAR models, see gao2019banded. Assumption (ref) implies that all moment conditions for the idiosyncratic terms in bai2004estimating are satisfied. Assumption (ref) is used to characterize the limiting behavior of the unit-root factors. Similar assumptions are used in bai2004estimating, zhang2019identifying, and gao2021modeling_IJF, among others. Assumptions (ref)-(ref) imply that all conditions for the common factors and the idiosyncratic terms in bai2004estimating hold. Assumption (ref) is used to control the sample covariance between the common factors and the idiosyncratic terms. The rate in Assumption (ref) is not strong and can be established under the setting of stock1987asymptotic, where we can assume the factors and idiosyncratic terms have similar structure as those in (2.4) therein. Assumption (ref) is the standard stationarity condition for a VAR process.
Turn to the convergence of the estimated loading matrix and its orthogonal complement. Note that the loading matrix $\bB$ is not uniquely defined due to the identification issue, only the linear space spanned by its columns, denoted by $\mathcal{M}(\bB)$, or the matrix product $\bB\bB'$ is uniquely defined. We state the convergence of the estimated loading matrix and its orthogonal complements in the following theorem.
The following theorem establishes the convergence of the estimated number of common stochastic trends.
Next, turn to the convergence of the estimated regression coefficients obtained in Section (ref). To control the errors between the estimated coefficients and the true ones, we introduce a Restricted Strong Convexity (RSC) condition which is often used in high-dimensional regularized estimation problems. See agarwal2012noisy and wainwright2019high for details. For any given $\lambda_{\bA},\lambda_{\bPhi}>0$, and a matrix $\bDelta\in \bbR^{p\times (N-r+dp)}=[\bDelta_1,\bDelta_2]$ with $\bDelta_1\in\bbR^{p\times (N-r)}$ and $\bDelta_2\in\bbR^{p\times dp}$, we use a weighted combination to define an associated norm as follows:
The restricted strong convexity condition under our setting is defined below.
When $\tau_T=0$, the RSC condition in Definition (ref) is called a locally strong convexity condition. See wainwright2019high. Denote $\bDelta=[\bDelta_{\bA},\bDelta_{\bPhi}]$ with $\bDelta_{\bA}=\wh\bA-\bA$ and $\bDelta_{\bPhi}=\wh\bPhi-\bPhi$. We now establish the convergence rates of the estimated coefficient matrices below.
To establish properties of the estimated coefficients using the IRRA of Section 2.2.2, we first introduce a restricted set that is constructed by a projection of any matrix onto a subspace generated by another one of the same shape. Specifically, for any $m\times n$ matrix $\mat{\Theta}$, we perform a singular value decomposition (SVD) $\mat{\Theta} = \mat{U}\mat{D}\mat{V}'$ with a partition as follows,
where $\mat{U}_k\in\mathbb{R}^{m\times k}$ and $\mat{V}_k\in\mathbb{R}^{n\times k}$ are the sub-matrices consisting of the left and right singular vectors associated with the $k$ largest singular values of $\mat{\Theta}$, respectively, and $\mat{U}_{k,c} \in\mathbb{R}^{m\times (m-k)}$ and $\mat{V}_{k,c} \in\mathbb{R}^{n\times (n-k)}$ are the remaining ones. Similarly to negahban2011estimation, we define two subspaces as follows,
For any matrix $\mat{M}\in\mathbb{R}^{m\times n}$, we decompose it as $\mat{M}=\mat{M}_1+\mat{M}_2$, where
Because $\mat{M}_2\in \mathcal{S}_\mat{\Theta}^\perp (k)$, we use $\Pi_{\mathcal{S}_\mat{\Theta}^\perp (k)}(\mat{M}) = \mat{M}_2$ to denote the projection of matrix $\mat{M}$ onto the subspace $\mathcal{S}_\mat{\Theta}^\perp (k)$.
Turn to the estimated coefficients using the IRRA of Section 2.2.2. By an abuse of notation, we define $\bDelta=[\bDelta_{\bA},\bDelta_{\bPhi}]=[\bDelta_{\bA},\bDelta_{\bPhi_1},...,\bDelta_{\bPhi_d}]$ with $\bDelta_{\bA}=\wh\bA-\bA\in \mathbb{R}^{p \times (N-r)}$ and $\bDelta_{\bPhi_i}=\wh\bPhi_i-\bPhi_i\in\mathbb{R}^{p \times p}$, and hence $\bDelta_{\bPhi}\in\mathbb{R}^{p\times dp}$. We decompose $\bDelta_{\bA}$ as $\bDelta_{\bA} = \bDelta_{\bA,1} + \bDelta_{\bA,2}$ and $\bDelta_{\bPhi_i}$ as $\bDelta_{\bPhi_i} = \bDelta_{\bPhi_i,1} + \bDelta_{\bPhi_i,2}$, for $1\leq i\leq d$. It follows that $\bDelta_{\bA,2} = \Pi_{\mathcal{S}_{\bA}^\perp(r_{\bA})}(\bDelta_{\bA})$, and $\bDelta_{\bPhi_i,2} = \Pi_{\mathcal{S}_{\bPhi_i}^\perp(r_i)}(\bDelta_{\bPhi_i})$, for $1\leq i\leq d$. We define a restricted set $\mathcal{C}$ as
We make an additional assumption below.
Note that Assumption 6 is a locally restricted strong convexity condition by setting $\tau_T=0$ in Definition (ref). Similar assumptions are also considered in Chapter 10 of reinsel2022multivariate for {\it i.i.d.} data. We next state the convergence of the estimated coefficients based on the IRRA of Section 2.2.2.
In this section, we evaluate the finite-sample performance of the proposed methodologies under the scenarios when both $p$ and $N$ are increasing from small to large. Though the dimensions of $\bB$ and $\wh\bB$ are not necessarily the same, as estimation error in $r$ may occur, the discrepancy measure adopted in Theorem (ref) remains valid. To simplify the presentation and without loss of generality, we set $d=1$ in (ref), and similar results can also be obtained for other choices of finite $d$.
We follow the data generating process in (ref) and (ref) and consider a three-factor model, where the factors are $I(1)$ processes generated by (ref). {We further multiply the factors by $\sqrt{N}$ because we will use orthonormal loading matrices below and the strength of general loadings is imposed on the factors in line with the assumptions and identification conditions.} While the number of factors $r=3$ is fixed, we set $p=20, 40, 60$, and $N=20, 40, 60$, respectively, and in each configuration of $(p,N)$, we set the sample size $T=400, 800, 1200$ to illustrate the proposed method and to exam certain theoretical properties of the estimators. In order to obtain reproducible results, we initialize a random generator in the NumPy package in Python by setting the seed to $1024$, and this seed is used throughout the simulation.
To begin, we need to obtain the coefficient matrices of the model. We start with generating the loading matrix $\bB$ and its corresponding orthogonal complement $\bB_c$. As $[\wh\bB,\wh\bB_c]$ is an $N\times N$ full-rank orthonormal matrix, we first randomly generate an $N\times N$ orthogonal matrix, and divide its columns in such a way that the submatrix with the first $r$ columns is chosen as $\bB$ and the remaining columns form naturally the $\bB_c$ matrix. For the low-rank matrix $\mat{A}$, we first randomly generate two orthonormal matrices $\mat{U}\in\bbR^{p\times p}$ and $\mat{V}\in \mathbb{R}^{(N-r)\times (N-r)}$, and a $p\times (N-r)$ rectangular diagonal matrix $\mat{D}$ with only five positive entries on the upper left of the diagonal while all the other entries are set to zero. The positive diagonal entries in $\bD$ are drawn independently from a uniform distribution on the interval of $[0.1,1)$ so that all the five elements are strictly greater than $0$. The matrix $\mat{A}$ with rank $r_{\bA}=5$ is then chosen as $\mat{A} = \mat{U}\mat{D}\bV'$. Next, for the sparse matrix $\bPhi$, we first create a sparse matrix $\bPhi_1$ with only $20$ randomly located non-zero entries each of which is drawn uniformly on the intervals $(-1,-0.1] \cup [0.1,1)$. In order to guarantee the stationarity of $\by_t$ in (ref), we use the normalized matrix $\bPhi = 0.9 \times \bPhi_1/\norm{\bPhi_1}_2$ as the autoregressive coefficient matrix, which implies that Assumption (ref) holds.
For each configuration of $(p,N,T)$, with the coefficient matrices $\bB, \bB_c, \bA$ and $\bPhi$ chosen by the aforementioned methods, we generate $\bx_t, \bz_t$ and $\by_t$ according to Models (ref), (ref) and (ref), respectively. To obtain stable results, we use $500$ replications for each $(p,N,T)$ configuration and set $\bve_t\sim N(\vect{0},\bI_N)$, $\vect{u}_t \sim N(\vect{0},\bI_r)$, and $\vect{e}_t \sim N(\vect{0},\bI_p)$ in each realization.
We first study the performance of ((ref)) in estimating the number of factors. Because the data generating process $\bx_t$ of the previous section is independent of the dimension $p$, we only illustrate the proposed method for the case of $p=20$, and similar results can also be obtained for other cases. Table (ref) reports the empirical probabilities of $P(\wh{r} = r)$ based on $500$ repetitions for each $(N,T)$ configuration when $p=20$, where we use the method described in Section (ref) with $\bar{k}=10$ and $\delta_0=0.3$. From Table (ref), we see that the auto-correlation based method can successfully recover the number of common stochastic trends. This is understandable because all the factors used in the simulation are strong ones. Similar results can also be found in bai2002determining and lam2012factor.
Next, we consider the estimation accuracy of the loading matrix $\bB$, which is measured by $\norm{\bB\bB' - \wh{\bB}\wh{\bB}'}_2$ over 500 replications. For the same reason mentioned before, we only show the results for the case of $p=20$. Boxplots of the discrepancies are shown in Figure (ref), from which we see that for each $N$, the discrepancy between the estimated loading matrix and the true one decreases as the sample size $T$ increases. This result is in agreement with our theorems. Furthermore, we also evaluate the estimation errors of the extracted factors. For each $(N,T)$ configuration, we define the the root-mean-squared-error (RMSE) of the estimated factors as
which quantifies the accuracy in recovering the common stochastic trends. Figure (ref) shows the results via boxplots using 500 replications. From Figure (ref), we see clearly that, the recovery errors of the common factors decrease as the sample size $T$ increases, which is consistent with the theoretical results in Theorem (ref).
We then study the estimation accuracy of the low-rank matrix $\bA$ and the sparse matrix $\mat{\Phi}$ using the procedure in Algorithm (ref). For simplicity, we set the tuning parameters $\lambda_{\bA} = \sqrt{(p+N)/T}$ and $\lambda_{\bPhi} = \sqrt{\log(p)/T}$, which are just taken from the rates discussed in Remark (ref)(ii) by setting $C_*=1$ and this choice is good enough to produce satisfactory performance in the simulation. In practice, we may choose an optimal $C_*$ from an interval using grid search. Due to the identification issue as that for $\bB$, we also use $\norm{\bA\bA' - \wh{\bA}\wh{\bA}'}_2$ to evaluate the discrepancy between $\wh\bA$ and $\bA$. Because there is no identification issues with $\bPhi$ and the estimated $\wh\bPhi$, we use $\norm{\mat{\Phi} - \wh{\mat{\Phi}}}_2$ to measure the estimation accuracy of the autoregressive coefficients. Boxplots of the estimation errors for $\wh\bA$ and $\wh\bPhi$ are presented in the Figures (ref) and (ref), respectively. As expected from Theorem (ref), in each case of $(p,N)$, the estimation errors of $\bA$ and $\bPhi$ both decrease as the sample size $T$ increases, which is also consistent with our theoretical properties.
Finally, we consider the estimation errors of the estimated explanatory variables and the true ones in Model ((ref)). Similarly to that in ((ref)), we define the RMSE for the regression model ((ref)) as
which is similar to the in-sample errors of a regression model. Figure (ref) displays boxplots of the RMSEs in ((ref)). From the plot, we see that the patterns of the boxplots are similar to those obtained before. For each given $(p,N)$, the RMSEs decrease as the sample size $T$ increases, illustrating the efficacy of the proposed method. Overall, the simulation results indicate that the proposed procedure works well in recovering the estimated coefficients.
In this example, we investigate the performance of IRRA of Section (ref). First, we generate the data $\bx_t$ using the same method as that of Section (ref). Second, unlike the sparse autoregressive matrices in Example 1, we generate two low-rank matrices $\bPhi_1 \in \bbR^{p\times p}$ and $\bPhi_2\in \bbR^{p\times p}$ under the context of the IRRA. Without loss of generality, we generate the those low-rank matrices in the same way as that of $\bA$, and set $\mathit{rank}(\bPhi_1)= \mathit{rank}(\bPhi_2)=3$. Third, the process $\by_t$ is then generated according to (ref) with the coefficients given above, where we choose $d=2$.
Similarly to the procedure in Example 1, we apply ((ref)) and Algorithm (ref) to estimate the number of factors and the coefficients, respectively. Since the performance of the auto-correlation based method is shown in Example 1, we omit the details here. Figures (ref), (ref) and (ref) show the discrepancies between the estimated coefficients and the true ones using Algorithm (ref). From these boxplots, we see that, for each configuration of $(p, N)$, all three coefficient estimates $\wh\bA, \wh\bPhi_1$ and $\wh\bPhi_2$ converge to the true ones as $T \to \infty$, which is consistent with our theory. For comparison, we also test the ADMM algorithm of li2019integrative, and find that the results of ADMM are quite close to those of the Algorithm (ref) in the sense that the distance $\norm{\wh{\bTheta}^{(\text{Ite})}-\wh{\bTheta}^{(\text{ADMM})}}_2/\norm{\wh{\bTheta}^{(\text{Ite})}}_2$ is less than $10\%$ in most cases, where $\wh\bTheta=\wh\bA, \wh\bPhi_1$ or $\wh\bPhi_2$, and $\wh{\bTheta}^{(\text{Ite})}$ and $\wh{\bTheta}^{(\text{ADMM})}$ are the coefficient matrices estimated by Algorithm (ref) and ADMM, respectively. Therefore, we omit the results obtained by the ADMM algorithm to save space.
In this section, we apply the proposed method to predicting monthly stock returns. welch2008comprehensive examined the predictability of some macroeconomic variables to the equity premium, and concluded that the performance of the predictions, both in-sample and out-of-sample, is poor and unstable. Using the same set of predictors, koo2020high exploited the cointegration relationship of the predictors, and showed that LASSO can improve the predictability of the macroeconomic variables in forecasting the equity premium of S&P 500 index. We use an extended data set and conduct a forecasting experiment using the proposed method. Note that we predict the stock returns of a cross-section, instead of the equity premium of an individual stock or index.
Consider the monthly returns of selected stocks in the S&P 500 index. Using the constituents of the index in January 2011 and the data structures in the CRSP (Center for Research in Security Prices) database, we select 79 stocks, which have no missing values during the time span from January 1960 to December 2019, as our sample. Therefore, we have 720 monthly observations of 79 $I(0)$ processes. We also collect the monthly macroeconomic variables in welch2008comprehensive as the predictors. An updated version of the data can be downloaded from Prof. Amit Goyal's personal website (\url{https://sites.google.com/view/agoyal145}), where we choose 13 macroeconomic predictors as the $I(1)$ processes $\bx_t$ from December 1959 to November 2019. Therefore, we have $N=13$, $p=79$ and $T=720$ in this illustration.
Table (ref) presents some descriptive statistics of the macro predictors, including their first-order sample autocorrelation coefficients $\rho(1)$ over the entire sample period. As shown in the table, nine predictors have a first-order sample autocorrelation coefficient higher than 0.95, but four variables (inflation, long-term yield, corporate bond returns, and stock variance) show little persistence. Therefore, we see that most of the variables are highly persistent and can be used as the $I(1)$ predictors in our model.
For the purpose of evaluating the forecasting performance of the proposed method, we adopt the out-of-sample $R^2$ measure commonly used in the literature regarding the prediction of stock returns, see gu2020empirical. At each time point, we define the out-of-sample $R^2$ as \[ R^2_{\mathrm{OOS}}(t) = 1 - \dfrac{\norm{\by_{t} - (\wh{\bA} \wh{\bz}_{t-1} + \wh{\bPhi} \bP_{t-1})}_2^2}{\norm{\by_{t}}_2^2}, \] which differs from that in gu2020empirical by being a function of the time index $t$, which denotes the forecasting origin.
Our empirical analysis works as follows. First, we divide the time span of 60 years into two periods. The first period is the initial estimation period from January 1960 to December 2010, and the second one is the testing period from January 2011 to December 2019. Specifically, we conduct the empirical test according to the following procedure, which is similar to that commonly used in the asset pricing literature. At the beginning, we use the data from January 1960 to December 2010 to estimate the coefficient matrices $\bA$ and $\mat{\Phi}$ of the model (ref), then predict the returns of the 79 stocks of January 2011 and calculate the out-of-sample $R^2$ of the prediction. We then add the returns of January 2011 to the estimation period, and refit the model (ref) with data from January 1960 to January 2011 to obtain updated coefficient matrices. The updated model is used to predict the returns of February 2011 with the model ((ref)) and to calculate the out-of-sample $R^2$ again. We repeat this estimation-prediction process by adding one-month returns to the estimation period in each iteration until November 2019, which enables us to predict the returns for December 2019. In addition, we choose the tuning parameters $\lambda_{\bA}$ and $\lambda_{\bPhi}$ based on the procedure described in Section (ref) but letting $\lambda_{\bA}$ and $\lambda_{\bPhi}$ be proportional to $\sqrt{(p+N)/T}$ and $\sqrt{\log(p)/T}$, respectively. See Remark (ref)(ii).
For comparison, we consider some alternative models commonly seen in the literature as benchmarks. The first benchmark is the naive VAR($d$) model, that is, \[ \by_t = \bPsi \bP_{t-1}+\be_t, \ t=1,2,\ldots,T. \] For each series of $\by_t$ and the corresponding row of $\bPsi$, we can treat the above equation as a simple regression problem with $dp$ explanatory variables and $T$ observations. Therefore, we may use the LS method to estimate each row of $\bPsi$, then put them together to construct an estimator of $\bPsi$. We use the model in koo2020high as another benchmark, where the tuning parameter $\lambda$ is fixed to $\frac{\log(p)}{10\sqrt{T}}$. Since the original model in koo2020high is developed for predicting a scalar time series, we apply their model with $13$ macroeconomic predictors described above to predict each stock separately, then stack the predictions together to calculate the out-of-sample $R^2$ for all $79$ stocks. Our final benchmark is the random walk model, in which we predict the returns of the next period using returns of the current period. For a more comprehensive comparison, we conduct the experiment for $d=1,2$ and $3$, respectively, for the naive VAR, the RRSRA, and the IRRA models.
We evaluate the empirical performance of the RRSRA in this section. To begin, Figure (ref) shows a time plot of the estimated number of common trends by the method described in Section (ref) in the testing period. The figure shows that, except for the first three months of 2011 that may be affected by some economic crisis, the estimated number of common trends within $\bx_t$ is four, which is fairly stable over the entire test period. Thus, we have nine cointegrating vectors to produce the stationary process $\wh{\bz}_t$ as a proxy of macroeconomic predictors. Before analyzing the forecasting performance of these estimated $\wh{\bz}_t$ variables, we take a look at the number of parameters to be estimated in the two coefficient matrices $\wh\bA$ and $\wh\bPhi$. Suppose that $\wh{r}=4$, then there are $9$ cointegrating vectors and hence, the matrix $\wh{\bA}$ has $79\times 9=711$ entries, and the matrix $\wh{\mat{\Phi}}$ has $79\times 79=6241$ entries to be estimated, both of which are relatively large. Therefore, we expect that the dimensions of the two matrices can further be reduced to low-rank or sparse ones, which is commonly assumed in the literature to avoid over-fitting and to produce better forecasting performance. For this reason, we expect that the tuning parameters $\lambda_{\mat{A}}$ and $\lambda_{\mat{\Phi}}$ in our framework should be relatively large to guarantee that the dimensions can be reduced.
Figure (ref) shows the estimated rank of $\wh{\mat{A}}$ and the estimated number of non-zero entries of $\wh{\mat{\Phi}}$ for the proposed model with $d=1$ at each prediction time point. The average rank of $\wh{\mat{A}}$ is $1.97$ and the average number of non-zero entries of $\wh{\bPhi}$ is $5.88$. Except for the first three months in 2011, the estimated rank of $\wh\bA$ is $2$ in the estimation period. In addition, $\wh{\bPhi}$ has at most $13$ non-zero entries at all time points, which is extremely small compared to $6241$ of the total number of entries. Overall, the proposed method provides an effective way to reduce the number of parameters and the dimension of the coefficient matrices.
Next we show the forecasting results in detail by following the method described in Section (ref) and the forecasting procedure mentioned above to evaluate the performance of different models. Table (ref) reports the overall comparisons of our proposed method against the three benchmarks mentioned before, in terms of $R^2_{\text{OOS}}$. From Panels A and C of Table (ref), we see that the proposed RRSRA substantially outperforms the naive VAR model (denoted by VAR($d$)), the method of koo2020high (denoted by Koo), and the random walk model (denoted by RW). The mean of out-of-sample $R^2$ of our method with $d=1$ is $0.91\%$ and the result is nearly the same for $d=2$ or $3$. We also note that our results are slightly better than those in gu2020empirical, where the highest monthly out-of-sample $R^2$ for all stocks is 0.40% among all machine learning methods considered in their paper, and it is $0.70\%$ for the top $1,000$ stocks and $0.47\%$ for the bottom $1,000$ stocks by market values. In particular, the VAR model performs relatively poorly, as the out-of-sample $R^2$s with different lags all assume negative values with large magnitudes, which are $-22.54\%$, $-41.08\%$ and $-87.81\%$, respectively. One possible reason is that the number of parameters to be estimated in VAR models is significantly large and this often leads to severe over-fitting, which in turn produces high variations in out-of-sample forecasting. When the lag order $d$ increases, the number of parameters also increases, so the performance would further deteriorate. For the model in koo2020high, the results in Table (ref) imply that the Koo method has limited predictive power when forecasting the returns of individual stocks. Finally, we see that the random walk model performs the worst. In summary, the proposed model has marked advantages in prediction over the three benchmark models considered.
To explore the in-sample goodness of fit, we apply all entertained models except the random walk to the entire data set, and calculate the in-sample $R^2$ at each time point. The results over the time are shown in Table (ref). As expected, the VAR model, which has many more degrees of freedom than the others, produces the highest in-sample $R^2$. Our model and the one of koo2020high provide a robust in-sample fit, and the result by the Koo method is only slightly worse than those of the proposed models.
To check whether our model outperforms the benchmarks uniformly over the in-sample and the out-of-sample periods, we plot the in-sample $R^2$s and out-of-sample ones of the RRSRA$(1)$ model in Figures (ref) and (ref), respectively, where the time index is on the horizontal axis. For a better illustration, the points in Figure (ref) are the in-sample $R^2$s based on the data of each year from 1960 to 2019. In both figures, we also plot the VAR(1) as a benchmark. An additional plot of the random walk model is also included in Figure (ref) as another benchmark. Because the results produced by the Koo method in koo2020high are very close to ours, they are omitted. From Figure (ref), we see that our method fits the data relatively poorly compared to the VAR model in most years according to the in-sample $R^2$. This is understandable since the VAR model fits the data via the LS method to minimize the squared distance between the fitted values and the true ones, while our method adopts regularization, which often introduces some in-sample biases in order to provide more stable predictions in out-of-samples. Furthermore, Figure (ref) shows that our method produces more robust predictions than the VAR and the random walk model, and outperforms them over most of the time points based on the out-of-sample $R^2$. This illustrates the predictive advantages of using the proposed method.
In this section, we evaluate the predictive performance of IRRA models described in Section (ref). The procedure of estimating the factor model (ref) and obtaining $\bz_t$ are exactly the same as those in Section (ref). With the estimated $\wh\bz_t$, we fit the data via (ref) using both iterative method and ADMM method, and evaluate its predictive performance. We expect that the results of IRRA are close to those of RRSRA, because the tuning parameter $\lambda_{\bPhi}$ selected by a grid search is relatively large in both models. When the tuning parameter $\lambda_{\bPhi}$ in both methods tends to infinity, the estimated coefficients obtained by the two algorithms tend to be the same.
We first examine the estimated coefficient matrices. Figure (ref) shows the rank of $\wh\bA$ and $\wh\bPhi$ estimated by Algorithm (ref) in the case of $d=1$ at each time point. We find that the rank of $\wh\bA$ is reduced to $1$ or $2$ over the entire time horizon, which is the same as that of the RRSRA, implying that the efficient cointegration rank is low in this particular application. In addition, the rank of $\wh\bPhi$ is also $1$ or $2$ over time.
Panel B of Table (ref) shows the predictive performance of (ref), where the coefficients are estimated by both Algorithm (ref) and the ADMM method with $d=1,2,3$, respectively. All six results are positive and close to each other. They are also close to, but a little worse than, those of RRSRA. One possible reason is that both IRRA and RRSRA are constrained regressions but the latter one produces sparse solutions and reduces the model complexity more substantially compared to the low-rank structures. Similarly to the conclusion of the RRSRA procedure, the $R^2_{\text{OOS}}$ results of (ref) also outperform those of the benchmarks.
Finally, we see that the performance of IRRA is close to that of RRSRA not only with respect to the out-of-sample $R^2$, but also with respect to the in-sample $R^2$; see Tables (ref) and (ref). Panel B of Table (ref) shows the in-sample $R^2$ results for IRRA, with all six estimation settings. Once again, we find that the results are close to those in Panel A, but IRRA fits the data slightly better, which may be a consequence of higher degrees of freedom in IRRA.
Finding proper cointegration relationships is an important topic in Econometrics and Statistics, yet the interpretation of cointegrating structures might become complicated if the dimension of the system under study is high. This paper introduced the concept of {\it effective cointegration rank} and considered a new method to identify the important cointegration relationships among a high-dimensional $I(1)$ series from a predictive perspective. In a nutshell, the effective cointegration rank is the number of cointegrating relationships that can produce useful predictors in a given forecasting application. The proposed method consists of a two-step estimation procedure, where we first use the Principal Component Analysis to estimate the common stochastic trends of the $I(1)$ series and to identify all possible cointegrating vectors. We then employ all stationary series obtained via the cointegrating vectors and some lagged values of dependent variables to form predictors of the second-step estimation. A reduced-rank regression technique is applied to the co-integrated predictors and the dimension of relevant cointegrating vectors is defined as the effective cointegration rank. We also applied the LASSO penalty or reduced rank constraints to the coefficients of the lagged variables in the second step, and an iterative procedure is proposed to estimate the unknown coefficients.
Our proposed method has a wide range of applications in many scientific areas, including Economics, Finance, and Environmental studies, because it is common in these areas to use nonstationary variables or factors to predict stationary series in empirical applications. We applied the proposed method to the problem of predicting cross-sectional stock returns, and illustrated clearly the predictive advantages of the proposed procedure over some commonly used benchmarks available in the literature.
\printbibliography