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.
124,386 characters · 13 sections · 61 citation commands
Optimal Bias-Correction and Valid Inference in High-Dimensional Ridge Regression: A Closed-Form Solution
\setstretch{1.3}
\setcounter{page}{2}
\abovedisplayskip=0.1pt \belowdisplayskip=0.1pt Regularization theory was one of the first signs of the existence of intelligent inference.{---Vladimir N. Vapnik (p.9, vapnik2013nature)}
Ridge regression, or more formally $\ell_2$-regularisation estimation, is a fundamental tool in econometrics, statistics, and machine learning with applications in many fields of science, technology, engineering, mathematics, medicine, social sciences, and humanities. The idea of $\ell_2$-regularisation appeared in the early 1940s for the stability of inverse problems; see tikhonov1943stability. It was first introduced to data analysis by hoerl1959optimum and later formulated in hoerl1970aridge,hoerl1970bridge for providing a robust solution to some of the persistent challenges encountered in traditional linear regression techniques; see hoerl1985ridge for a nice review. Emerging as a fundamental technique in predictive modelling, ridge regression addresses issues such as multicollinearity and overfitting, which commonly afflict predictive models dealing with high-dimensional data. Since its inception, ridge regression's practical adoption persists due to its superior performance over the least-squares estimator in various scenarios, evident in applications across neuroscience, chemistry, biology, and economics; see leonard2023large, zahrt2019prediction, otwinowski2014inferring, giannone2021economic, and abadie2019choosing, among others, underscoring its empirical effectiveness. From a shrinkage perspective, the ridge estimator also dominates the least-squares solutions in the sense that its mean-squared errors (MSEs) can be smaller, which provides a reasonable explanation on the empirical effectiveness of ridge estimators. See theobald1974generalizations, athey2019machine, hastie2020ridge, hansen2022econometrics, and a comprehensive introduction of ridge regression in van2023lecture.
The ridge estimator offers a closed-form expression that simplifies both theoretical and empirical analyses. It aligns with the dense modelling techniques of giannone2021economic, which acknowledge the potential significance of all explanatory variables for prediction. Empirical studies, such as those in giannone2021economic, indicate that dense models generally tend to outperform the sparse ones in out-of-sample economic prediction. Similarly, abadie2019choosing find that the ridge estimators dominate the Lasso and the pre-testing estimators in terms of risks when the effects of different predictors on the dependent variable are “smoothly distributed”. These results suggest that ridge estimators indeed constitute a crucial tool in statistical modelling and economic forecasting, especially in the big data era.
However, as highlighted in Section 2.8 of athey2019machine, constructing valid confidence intervals remains a challenge for many regularised methods, including ridge regression, even in asymptotic settings. This long-standing challenge in ridge-type regression involves at least two critical aspects within a linear regression framework: (a) performing hypothesis testing on specific linear combinations of the regression coefficients using the ridge estimators; and (b) deriving confidence or prediction intervals based on the ridge estimators in empirical applications. The primary reason for these challenges is that the inherent bias of ridge estimators poses significant challenges, compromising both statistical efficiency and scalability across various applications. To date, the feasibility of conducting statistical inferences and hypothesis testing on ridge-type estimators without imposing additional structural constraints remains largely unexplored in the literature. This complexity arises from the ridge estimator's intrinsic bias, which complicates direct statistical inferences despite its elegant closed-form expression. As a result, research addressing the inference challenges of ridge regression in high-dimensional settings is limited, leading to its widespread application across disciplines without comprehensive theoretical investigations.
To the best of our knowledge, there are only a few works in the literature concerning the bias and inference of ridge estimators under different scenarios. Assuming a sparse structure, shao2012estimation proposed a threshold ridge regression method and proved its consistency. The method therein actually estimates the projected coefficient vector rather than the true one in the linear model. dobriban2018high derived the limit of the predictive risk of ridge regression and regularised discriminant analysis in a dense random effects model. buhlmann2013statistical proposed to use the Lasso to correct the bias of ridge estimators. However, the estimation depends on the existence of an initial estimator which must be sufficiently accurate. zhang2022ridge adopted a similar approach as that in shao2012estimation and proposed a threshold ridge regression and a bootstrap method to make inferences. However, the method therein still estimates the projected coefficient vector rather than the true one, and the remaining bias may not be asymptotically negligible in general.
In this paper, we introduce a systematic approach tailored specifically for mitigating the bias in ridge-type estimators for high-dimensional linear regression models. Leveraging the closed-form expression of the ridge estimators, the bias term can also be established in an analytic form. Although the bias term involves the true parameters which are unknown in practice, we found that replacing the true parameters with the ridge estimators turns out to be an effective way to mitigate the bias. Therefore, the proposed method employs an iterative bias-correction strategy, and the bias can be reduced substantially when the number of iterations is sufficiently large. Notably, it achieves complete bias correction if the covariate dimension $p$ is smaller than the sample size $n$, and can reduce considerable bias if $p$ surpasses $n$. We show that our bias-correction method is an optimal one in the sense that the bias can be completely corrected (asymptotically) when $p< n$ with a sufficient number of the proposed bias-correction iterations, and the remaining bias of the de-biased ridge estimator when $p>n$ is unattainable through any linear transformations of the response vector.
To further combat the remaining bias in the de-biased ridge estimators when $p>n$, we introduce a novel ridge-screening (RS) method for selecting covariates prior to applying our bias correction procedure. The RS approach constructs a restricted model that inherently encompasses the true model as a subset. This is based on the assumption that only a subset of the covariates holds significance in the linear model. Specifically, we postulate that the number of significant covariates should be less than the sample size, which can appropriately diverge relative to the dimension $p$ and sample size $n$. Crucially, we demonstrate that under certain mild conditions, the true model is inherently nested within the selected one, establishing RS as a novel variable selection approach that offers additional value beyond bias correction. Leveraging this restricted model, our bias-correction procedure can be further applied to the restricted model and effectively rectifies the bias in the resulting ridge estimators
The derivation of the proposed methodology and theory is mainly based on a fixed design, which is similar to the setting in hansen2022modern, but the results throughout this paper remain valid for random regressors by conditioning on the design matrix. Thus, this paper rigorously establishes asymptotic properties and provides valid inferences of the proposed bias-corrected estimators for both $p< n$ and $p>n$ under some relaxed and intuitive conditions. Furthermore, we delve into the bias-variance trade-off of our de-biased ridge estimators, examining its relationship with the number of iterations in the bias-correction procedure. To validate our methodology, we provide empirical evidence using both simulated and real-world data. The prediction intervals constructed by the proposed method are indeed satisfactory for forecasting the U.S. macroeconomic series using factor-augmented regression. Moreover, the RS method can further enhance the coverage rates of these prediction intervals. These empirical validations highlight the practical efficacy and adaptability of the proposed approach across a wide range of high-dimensional regression settings.
The contributions of this paper are multi-fold. First, the proposed bias-correction method is simple and easy to implement. In fact, the proposed approach is a systematic procedure and the resulting de-biased estimator has a closed-form expression. Second, the optimality of the proposed bias-correction procedure consists of two aspects: (a) it achieves complete bias correction when the covariate dimension $p$ is smaller than the sample size $n$; and (b) the remaining bias of the de-biased ridge estimator in the scenario when $p>n$ is unattainable through any linear transformations of the response vector. These results distinguish our work from the existing ones that only part of the bias can be corrected in most of the aforementioned literature. Third, we also propose a novel variable selection procedure, namely the ridge-screening (RS) method, which screens out some insignificant variables based on the de-biased ridge estimator due to its optimality. Fourth, we establish certain asymptotic properties of the de-biased ridge estimators in both scenarios when $p< n$ and $p>n$, and it is shown that the de-biased ridge estimators are asymptotically normal, which provides valid inference methods for the ridge estimators. Fifth, we develop a procedure to construct confidence and prediction intervals in ridge regressions using the proposed de-biased estimators and associated inference methods. Finally, we establish the bias-variance trade-off of the proposed approach both theoretically and through validation with simulated data. It's important to note that, unlike the scenario described in Section 2.8 of athey2019machine where many inference approaches for regularised machine learning methods compromise predictive performance, our proposed procedure focuses on correcting the bias and rendering the estimators suitable for statistical inferences, without adversely affecting their predictive performance.
We highlight that the asymptotic framework adopted in this paper is slightly different from traditional approaches in the literature. Typically, in conventional frameworks, the asymptotic properties are established when the dimension and/or the sample size are approaching infinity. However, in our study, most of the asymptotic results are derived under a different scenario that the number of iterations in the bias-correction procedure tends to infinity for any given configuration of the dimension $p$ and the sample size $n$. This configuration can be sufficiently large to encompass the framework of big data analysis. One of the primary motivations behind this choice is to demonstrate the validity of the proposed procedure by explicitly providing the exact bias and covariance terms of the de-biased estimators in this paper. This approach is also reasonable because it reflects common scenarios encountered in empirical data analysis, where datasets often have fixed dimensions and sample sizes. Without this setting, the dimensions of the bias and covariance terms of the de-biased estimators would expand to infinity if we considered $n$ and $p$ approaching infinity, making them challenging to formulate and describe theoretically. Moreover, our asymptotic results remain valid as $n$ and $p$ approach infinity. This holds true so long as we initially allow the number of iterations to increase towards infinity at a moderate rate. In this asymptotic manner, we can symbolically retain the forms of the bias and covariance terms for a growing configuration of $(p, n)$. Therefore, it's important to emphasise that our chosen asymptotic framework does not undermine the validity of the proposed bias-correction procedures.
The rest of the paper is organised as follows: Section (ref) introduces the ridge estimation and its bias-correction procedure in scenarios when $p< n$ and $p>n$, where a ridge-screening method is also introduced for variable selection. Section (ref) presents the inference methodologies for the proposed de-biased estimators. Section (ref) examines the finite-sample performance pf the proposed method via Monte-Carlo simulations. Section (ref) provides an empirical application of the proposed method, and Section (ref) offers some concluding insights. All proofs and derivations for the asymptotic results are available in an online Supplementary Material.
{\bf Notation:} We use the following notation. For a $p\times 1$ vector ${\mathbf u}=(u_1,..., u_p)'$, $\|{\mathbf u}\|_2=\sqrt{\sum_{i=1}^p|u_i|^2}$ is the $\ell_1$-norm and $\|{\mathbf u}\|_\infty=\max_{1\leq i\leq p}|u_i|$ is the $\ell_\infty$-norm. ${\mathbf I}_p$ denotes the $p\times p$ identity matrix. For a matrix ${\mathbf H}$, its operator norm is $\|{\mathbf H} \|_2=\sqrt{\lambda_{\max} ({\mathbf H}' {\mathbf H} ) }$, where $\lambda_{\max} (\cdot) $ denotes the largest eigenvalue of a matrix, and $\|{\mathbf H}\|_{\min}$ is the square root of the minimum non-zero eigenvalue of ${\mathbf H}{\mathbf H}'$. $|{\mathbf H}|$ denotes the absolute value of ${\mathbf H}$ elementwisely. The superscript ${'}$ denotes the transpose of a vector or matrix. We also use the notation $a\asymp b$ to denote $a=O(b)$ and $b=O(a)$ or $a$ and $b$ are of the same order.
Let $\{({\mathbf x}_1,y_1),...,({\mathbf x}_n,y_n)\}$ be a given sample of centred observable data. We consider the problem of estimating a $p$-dimensional vector $\boldsymbol{\beta}$ from the following linear model:
where $y_i$ is a scalar response variable, ${\mathbf x}_i$ is a $p$-dimensional covariate vector, and ${\varepsilon}_i$ is a random error term with mean zero and finite variance. Similar to the setting in hansen2022modern and Ch. 29 of hansen2022econometrics, we treat the $n\times p$ design matrix ${\mathbf X}=({\mathbf x}_1,...,{\mathbf x}_n)'$ consisting of $p$ covariates as a fixed one. But the estimation results throughout the paper remain valid for random regressors by conditioning on the design matrix data ${\mathbf X}$. Note that Model ((ref)) can be expressed in vector form
where ${\mathbf y}=(y_1,...,y_n)'$ is an $n$-dimensional response vector, and $\mbox{\boldmath$\varepsilon$}=({\varepsilon}_1,...,{\varepsilon}_n)'$ is an $n$-dimensional vector of noises with $E(\mbox{\boldmath$\varepsilon$})={\bf 0}$ and $\mbox{Cov}(\mbox{\boldmath$\varepsilon$})=\boldsymbol{\Sigma}_{\varepsilon}$, where $\boldsymbol{\Sigma}_{\varepsilon}$ is a diagonal matrix with positive and bounded diagonal elements. If $p< n$ and ${\mathbf X}'{\mathbf X}$ is an invertible matrix, the least-squares estimator for $\boldsymbol{\beta}$ is
This least-squares estimator $\widehat\boldsymbol{\beta}_{lse}$ is only well-defined if $({\mathbf X}'{\mathbf X})^{-1}$ exists. In a high-dimensional setting, if the columns of the design matrix ${\mathbf X}$ are exactly linearly dependent, for example, this is obviously true when $p>n$, this collinearity among the columns implies that ${\mathbf X}'{\mathbf X}$ is singular, rendering $\widehat\boldsymbol{\beta}_{lse}$ an invalid estimator. To make the least-squares estimator in ((ref)) a well-defined quantity, we modify the definition in ((ref)) as
where $({\mathbf X}'{\mathbf X})^{+}$ is the Moore-Penrose generalised inverse. It is not hard to see that the estimator in ((ref)) reduces to the one in ((ref)) if $p<n$ and ${\mathbf X}'{\mathbf X}$ is invertible. Therefore, we will denote the estimator in ((ref)) as the least-squares solution throughout this article.
The ridge regression estimator was first proposed by hoerl1959optimum; see the review article of hoerl1985ridge. It essentially comprises of an ad-hoc fix to resolve the singularity issue of ${\mathbf X}'{\mathbf X}$ in the presence of many covariates. Suppose $({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)$ is invertible for a given $\lambda>0$, the ridge estimator is defined as
which simply replaces ${\mathbf X}'{\mathbf X}$ by ${\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p$ with a tuning parameter $\lambda>0$ in the least-squares estimator of ((ref)). From a regression point of view, the ridge estimator can be obtained in the following way. Let ${\mathbf y}^*=({\mathbf y}',{\bf 0})'$ and ${\mathbf X}^*=({\mathbf X}',\sqrt{\lambda}{\mathbf I}_p)'$ be the augmented data, the ridge estimator is a solution to the following optimisation problem:
From the expression of the ridge estimator for a given $\lambda>0$, it is not hard to see that \[\widehat\boldsymbol{\beta}(\lambda)\rightarrow\widehat\boldsymbol{\beta}_{lse},\,\, \text{as}\,\, \lambda\rightarrow 0,\] and \[\lambda\widehat\boldsymbol{\beta}(\lambda)\rightarrow{\mathbf X}'{\mathbf y},\,\,\text{as}\,\, \lambda\rightarrow \infty,\] which is the componentwise regression estimator if each covariate is standardised. Therefore, a large $\lambda$ would reduce the variance of the estimator, but may increase the bias. In fact, the bias increases as $\lambda$ grows.
In this paper, we only focus on a given $\lambda>0$ and investigate the bias-correction and inference issues for $\widehat\boldsymbol{\beta}(\lambda)$. For the purpose of comparisons with the proposed approach, we first specify the initial bias of the ridge estimator of ((ref)) in the following theorem.
From Theorem 1, the bias of the ridge estimator depends on the unknown true parameter $\boldsymbol{\beta}$. Therefore, it is fundamentally challenging to make any statistical inference on the ridge estimator or to construct any confidence or prediction intervals involving the ridge estimator. Consequently, it is important to seek an effective way to correct or reduce the bias of a ridge estimator.
In the next section, we will tackle this issue by proposing an iterative procedure to reduce the bias of the ridge estimator.
We divide the discussion into two key parts depending on whether ${\mathbf X}'{\mathbf X}$ is invertible or singular. For simplicity and in line with the setting in wang2016high, we assume that ${\mathbf X}'{\mathbf X}$ is invertible when $p< n$ and singular when $p>n$ throughout this article. An alternative framework is to consider the scenarios that $p/n\in (0,1)$ and $p/n\in(1,\infty)$, a common setting in random matrix theory, which is also helpful in establishing the asymptotic results if we further allow $n,p\rightarrow\infty$ later. This assumption is well-justified as we consider a fixed design matrix ${\mathbf X}$ in this paper and ${\mathbf X}'{\mathbf X}$ naturally maintains its invertibility when $p< n$. Furthermore, in extreme cases where highly correlated variables are present within ${\mathbf X}$ (in a random sense), we may implement specific transformations based on prior knowledge or statistical methods such as the hierarchical clustering or the $k$-means algorithm to mitigate these correlations prior to conducting ridge regression; see the discussion in Section 4.1.2 of fan2008sure. Consequently, we only focus on the bias-correction issue in this paper and rule out the case when some covariates in ${\mathbf X}$ are highly correlated.
The ridge estimator in ((ref)) can be written as
The bias-correction procedure is based on this expression and ((ref)) in Theorem (ref). The rationale for the procedure is as follows. Since $\boldsymbol{\beta}$ in the bias term of ((ref)) is unknown, we first replace it by $\widehat\boldsymbol{\beta}(\lambda)$ defined in ((ref)) and construct a first-step de-biased ridge estimator as
Plugging $\widehat\boldsymbol{\beta}(\lambda)$ from ((ref)) into ((ref)), we obtain:
Consequently, the bias term of $\widehat\boldsymbol{\beta}_{c,1}(\lambda)$ is
It is then apparent that the $\ell_2$-norm of the bias term ${\mathbf b}_{\lambda,1}$ produced by $\widehat\boldsymbol{\beta}_{c,1}(\lambda)$ is smaller than that of the initial bias ${\mathbf b}_{\lambda,0}$ in Theorem (ref) under some mild conditions, implying that the bias ${\mathbf b}_{\lambda,0}$ has been partially corrected by $\widehat\boldsymbol{\beta}_{c,1}(\lambda)$. To see this, we conduct a singular-value decomposition on ${\mathbf X}$ or a spectral decomposition on ${\mathbf X}'{\mathbf X}$, and the effectiveness of the bias-correction approach depends on two observations: (a) the eigenvalues of $\lambda({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-1}$ are positive and strictly less than one; and (b) the eigenvalues of $\lambda^2({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-2}$ becomes smaller in the de-biased estimator of ((ref)) compared to those of $\lambda({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-1}$ in ${\mathbf b}_{\lambda,0}$.
Following the first correction, we replace the unknown vector $\boldsymbol{\beta}$ in ((ref)) by the ridge estimator $\widehat\boldsymbol{\beta}(\lambda)$ again to obtain a second-step de-biased estimator:
By a similar argument, the bias term of $\widehat\boldsymbol{\beta}_{c,2}(\lambda)$ is
where the eigenvalues of $\lambda^3({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-3}$ are even smaller compared to those of $\lambda^2({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-2}$ in ${\mathbf b}_{\lambda,1}$. Consequently, we can show that the $\ell_2$-norm of the bias ${\mathbf b}_{\lambda,2}$ is smaller than that of ${\mathbf b}_{\lambda,1}$ under the same framework as before. We repeat this procedure and denote
as a de-biased estimator at the $k$-th step, where we define $\widehat\boldsymbol{\beta}_{c,0}(\lambda)=\widehat\boldsymbol{\beta}(\lambda)$ for $k=0$. To characterize the effect of the bias correction, we make an assumption on the singular-value decomposition (SVD) of ${\mathbf X}$ first.
Assumption (ref) is intuitive for a fixed design ${\mathbf X}$ with a large and fixed configuration of $(p,n)$. For example, the eigenvalues of ${\mathbf X}'{\mathbf X}$ are of order $n$ if the entries of ${\mathbf X}$ are independent copies of a random variable with zero mean, unit variance, and finite fourth moment if $p/n\in (0,1)$; see the Bai-Yin's law in bai1993limit. Consequently, the eigenvalues of $\lambda({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-1}$ are strictly less than one, which ensures that our iterative procedure can substantially reduce the bias. In fact, we have the following theorem for the de-biased ridge estimator in ((ref)).
Theorem (ref) implies that we can completely correct the bias term incurred by the ridge estimator $\widehat\boldsymbol{\beta}(\lambda)$ if $p<n$ so long as we conduct a sufficient number of iterations. This result is particularly useful for empirical data analysis when $n$ and $p$ are given.
To investigate the performance of the de-biased estimator $\widehat\boldsymbol{\beta}_{c,k}(\lambda)$ when $p>n$ under which ${\mathbf X}'{\mathbf X}$ is singular, we first perform a singular-value-decomposition on ${\mathbf X}$. Suppose the true rank of ${\mathbf X}$ is $\text{rank}({\mathbf X})=p^*\leq \min(p,n)=n$, where we can simply set $p^*=n$ since we deal with a given and fixed configuration of $(p,n)$ without including highly correlated covariates. However, the results in Theorem (ref) still hold for $p^*<n$. By abuse of notation, there exist semi-orthogonal matrices ${\mathbf V}_1\in R^{n\times p^*}$ and ${\mathbf U}_1\in R^{p\times p^*}$, and a diagonal matrix ${\mathbf D}_1=\mbox{diag}(d_1,...,d_{p^*})$ with $d_1\geq ...\geq d_{p^*}>0$ such that
Since ${\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p$ is symmetric and invertible, there also exists an orthogonal complement matrix ${\mathbf U}_2\in R^{p\times (p-p^*)}$ of ${\mathbf U}_1$ such that
We formulate the above description in Assumption (ref) below.
Then, we have the following theorem.
A key insight from Theorem (ref) is that there is a remaining bias term ${\mathbf U}_2{\mathbf U}_2'\boldsymbol{\beta}$ that cannot be corrected by the proposed method. This is understandable since the projection of $\boldsymbol{\beta}$ on the singular directions ${\mathbf U}_2$ is not captured in Model ((ref)) according to the singular-value-decomposition of ${\mathbf X}$ in ((ref)). If the vector $\boldsymbol{\beta}$ belongs to the space spanned by the columns of ${\mathbf U}_1$, there will be no need to correct the bias since ${\mathbf U}_2'\boldsymbol{\beta}={\bf 0}$. Otherwise, it is challenging to make such a correction for the bias in Theorem (ref). As a matter of fact, we have the following theorem regarding the bias term in Theorem (ref).
We focus on linear transformations of the data in Theorem (ref) because the least-squares estimator and the ridge estimator are all linear combinations of the data ${\mathbf y}$. Theorem (ref) indicates that the remaining bias term ${\mathbf U}_2{\mathbf U}_2'\boldsymbol{\beta}$ is unattainable through any linear transformation of the data ${\mathbf y}$, showing that the proposed approach does its best to correct the bias. The results in Theorems (ref)--(ref) also indicate that our bias-correction method is an optimal one for any given $(p,n)$.
Simulation results (e.g., Figure (ref)) in Section (ref) suggest that the uncorrectable bias ${\mathbf U}_2{\mathbf U}_2'\boldsymbol{\beta}$ can be significant even when the number of nonzero elements in $\boldsymbol{\beta}$ is small. This underscores that any method that only corrects part of the bias may lead to poor inference performance.
To address further the uncorrectable bias identified in Theorem (ref) and to ensure valid statistical inferences when $p>n$, we impose some additional structures on the parameters or covariates. Without these additional structures, as demonstrated in Theorem (ref), the bias-correction becomes challenging. In the subsequent analysis, we use $C$ or $c$ to denote a generic constant, the specific value of which may vary across different contexts.
In this section, we propose a ridge-screening approach to select the significant variables in Model ((ref)). It is feasible for both scenarios when $p<n$ and $p>n$. Therefore, it can also be treated as a new variable selection approach especially useful in a high-dimensional setting, which is of independent interest to statisticians, econometricians, and data scientists. We only focus on the scenario when $p>n$ and ${\mathbf X}{\mathbf X}'$ is singular. Note that the reason for the remaining bias term in Theorem (ref) that cannot be corrected is that some projection directions of the coefficients $\boldsymbol{\beta}$ are not captured in Model ((ref)). For example, this is the case when some parameters in $\boldsymbol{\beta}$ are redundant if certain covariates in ${\mathbf X}$ are strongly correlated. Therefore, it is reasonable to assume that some covariates in ${\mathbf X}$ are not useful in the linear regression ((ref)), and hence, they can be dropped before establishing valid ridge estimators. To embrace the sparsity assumption in high-dimensional data analysis, we assume the true parameter vector $\boldsymbol{\beta}$ belongs to the following submodel class
with cardinality $|\mathcal{M}_0|=s^*< \min(p,n)$, where $\kappa$ is a small positive number. It is obvious that we require $s^*< p^*$, which is the rank of ${\mathbf X}$ defined in Assumption (ref). This cardinality assumption is natural and it reflects the idea that many coefficient parameters are relatively small among the $p$-dimensional vector $\boldsymbol{\beta}$. We treat them as zero elements only for ease of exploitation. According to the results in Theorem (ref) and Theorem (ref), we have that
as $k\rightarrow\infty$ for a given $\lambda>0$.
It is intuitive to expect that the components of $\widehat\boldsymbol{\beta}_{c,k}(\lambda)$ corresponding to positions in the submodel $\mathcal{M}_0$ will be greater than those at positions in $\mathcal{M}_0^c$, the complement set of $\mathcal{M}_0$. This is due to the following reason: Assuming the number of nonzero elements in $\boldsymbol{\beta}$ is finite or relatively small compared to $p$, then the $\ell_2$-norm of the $p$-dimensional dense vector ${\mathbf U}_2{\mathbf U}_2'\boldsymbol{\beta}$ is also of finite or relatively small order. Consequently, the magnitude of each projected coordinate in ${\mathbf U}_2{\mathbf U}_2'\boldsymbol{\beta}$ is of a smaller order compared to the nonzero elements in $\boldsymbol{\beta}$ on average. Therefore, we propose a Ridge-Screening (RS) method that selects the submodel class
where $\widehat\boldsymbol{\beta}_{c,k}(\lambda^*)=(\widehat\beta_{c,k,1}(\lambda^*),...,\widehat\beta_{c,k,p}(\lambda^*))'$ and we use a different penalty $\lambda^*>0$ to distinguish it from the one utilised in the subsequent step. We clarify the notation used here. For $p>n$, $s^*$ is the number of nonzero elements in $\boldsymbol{\beta}$, $p^*$ denotes the rank of ${\mathbf X}$, and $n^*$ is the number of covariates selected by the RS method. Obviously, we require that $s^*<n^*\leq p^*\leq\min(p,n)=n$ for $p>n$. We can also simply set $p^*=n$ since we deal with a fixed design.
The ranking method to derive the submodel in ((ref)) is similar to the approach presented in Section 2.2 of fan2008sure. However, the method in fan2008sure is based on marginal correlations between the response and the features, while ((ref)) is a utilisation of a de-biased estimator, which is more proximate to the true parameter than the conventional ridge estimator in Eq. (5) of fan2008sure. Furthermore, our methodology is not constrained by a specific choice of $\lambda^*$, and the optimal one can be chosen by an information criterion or a cross-validation method as discussed at the end of Section (ref).
In practice, we may select $n^*< \min(p,n)=n$ in ((ref)) through cross-validation along with $\lambda^*$ because the actual design matrix with $n^*$ columns of covariates is no longer singular in such cases. Additionally, we can show that the model with $n^*$ covariates asymptotically encompasses the submodel $\mathcal{M}_0$ in ((ref)) under some mild conditions, i.e., the true submodel $\mathcal{M}_0$ in ((ref)) is nested within the one with $n^*$ covariates.
Next, we restrict the design matrix to the submodel class $\mathcal{M}_k(\lambda^*)$ and denote the restricted design as ${\mathbf X}_{\mathcal{M}_k}\in R^{n \times n^*}$. The new ridge estimator becomes
where $n^*< n$, and $\lambda$ can be different from the $\lambda^*$ used in the ridge screening approach in ((ref)). When $\lambda^*$ and $n^*$ are optimally selected using the method outlined at the end of Section (ref) below, we can set $\lambda=\lambda^*$ because it is optimal for the restricted ridge estimator with a given set of $n^*$ significant variables. We observe that the eigenvalues of ${\mathbf X}_{\mathcal{M}_k}'{\mathbf X}_{\mathcal{M}_k}$ are strictly greater than zero for a given configuration of $(p,n)$. This essentially reduces the scenario to the case when $p<n$ in Section (ref). Consequently, we can employ the bias-correction procedure detailed in Section (ref) over $l$ iterations to obtain a de-biased estimator:
To establish the asymptotic properties of the RS method and the de-biased estimator for the restricted one in ((ref)), we make a few intuitive assumptions below.
Assumptions (ref) and (ref) are natural conditions about the orders of the singular values of ${\mathbf X}\in R^{n\times p}$ and ${\mathbf X}_{\mathcal{M}_k}$, and the penalty parameter $\lambda^*$ (or $\lambda$) is comparable to the magnitude of $d_j^2$ (or $d_{\mathcal{M}_k,j}^2$). Similar to the illustrations for Assumption (ref), the magnitude in Assumptions (ref) and (ref) can be easily verified if the entries of ${\mathbf X}$ are independent copies of a random variable with zero mean, unit variance, and finite fourth moment under the setting $p/n\in(1,\infty)$; see bai1993limit. In fact, the orders specified in Assumptions (ref) and (ref) are only employed to establish the validity of the ridge-screening method in Theorem (ref) below, and they are not the only ones capable of achieving this, so long as the rates can be properly controlled in the proof of Theorem (ref). For any fixed $(p,n)$, which can be large, the efficacy of the bias-correction method and the inference methods in Section (ref) below remain valid so long as the nonzero singular values and the penalty parameters are strictly greater than zero. Assumption (ref) indicates that the minimum nonzero element in $\boldsymbol{\beta}$ cannot be too small, and the norm of the $i$-th projected coordinate $({\mathbf U}_2{\mathbf U}_2'\boldsymbol{\beta})_i$ is bounded by the magnitude of its original coordinate $|\beta_i|$ (up to a small constant). This is actually a reasonable and intuitive assumption. For instance, in Assumption (ref), if $\|\boldsymbol{\beta}\|_2\leq C\sqrt{s^*}$ (at most) because there are only $s^*$ nonzero elements in $\boldsymbol{\beta}$, then $\|{\mathbf U}_1{\mathbf U}_1'\boldsymbol{\beta}\|_2=O_p(\sqrt{s^*})$, but it is a $p$-dimensional vector, meaning that the magnitude of each $({\mathbf U}_1{\mathbf U}_1'\boldsymbol{\beta})_i$ is of order $\sqrt{s^*/p}$ on average. We may postulate that $\sqrt{s^*/p}=o(n^{-\tau})$ in such a case and Assumption (ref) is even slightly weaker than this situation as we can allow that the projected coordinate is of the same order as that of the original one. Assumption (ref) is a general condition that includes Gaussian distributions as a special case.
We have the following theorem regarding the de-biased estimator in ((ref)) after the ridge-screening.
To conclude this subsection, we provide a summary of the pseudo-code for the proposed bias-correction procedures in Algorithm (ref) and Algorithm (ref), corresponding to the scenarios of $p < n$ and $p > n$, respectively. For the convergence criterion in Algorithm (ref), we can choose a small threshold $\eta>0$, say $\eta=10^{-2}$, and the convergence of the algorithm is determined by checking whether the following inequality hold:
It is important to note that the convergence is guaranteed because the non-random sequence ${\mathbf A}_k:=\sum_{j=1}^k\lambda^j({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-j}$ converges, as demonstrated in the proof of Theorem (ref) in the Appendix. Additionally, we emphasise that the number of iterations is typically not large, as only a logarithmic order of the dimension is required to ensure convergence, as discussed in Remark (ref). Furthermore, the computation is efficient since we only need to perform the SVD on ${\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p$ once, as we do for the classical ridge estimator. After this single SVD, the iterative process can be applied to the diagonal singular-value matrix in the expansion.
Finally, we briefly discuss the methods for selecting the unknown parameters in Algorithms (ref) and (ref) of the proposed approaches. Firstly, the sole unknown parameter in Algorithm (ref) is $\lambda$. Throughout this paper, we assume $\lambda$ to be given, as our focus is mainly on the bias-correction issue for ridge estimators. Empirically, one can employ the widely-accepted information criteria (e.g., AIC or BIC) or utilise cross-validation methods to determine an appropriate $\lambda$. For a comprehensive understanding, readers are referred to Section 1.8 of van2023lecture. Secondly, in the ridge-screening method of Algorithm (ref) when $p>n$, the parameters $(\lambda^*,n^*)$ are unknown. Here, one can employ the aforementioned information criteria or cross-validation methods, as discussed in Sections 1.8.1-1.8.3 of van2023lecture, to simultaneously determine $(\lambda^*,n^*)$. This can be achieved by selecting a candidate for $\lambda^*$ and then determining the optimal $n^*$ such that the pair $(\lambda^*,n^*)$ minimises the information criterion or yields the best prediction performance on test sets using the restricted ridge estimator $\widehat\boldsymbol{\beta}_{\mathcal{M}_k}(\lambda^*)$. The final choice of $(\lambda^*,n^*)$ can be determined by evaluating the performance across each grid point of the penalty parameters. After identifying the significant variables through the RS method, the parameter $\lambda$ used in bias-correction is equal to $\lambda^*$, as $\lambda^*$ represents the optimal choice for the associated $n^*$. Given that these methods are well-established in the literature, we omit further details here to save space.
In this section, we briefly introduce the inference method for the de-biased ridge estimators in Sections (ref) and (ref). The following theorem can be derived immediately based on the proofs of Theorems (ref) and (ref).
Hence, by correcting the bias, we can proceed to make statistical inferences and construct confidence intervals using distributions in ((ref)) and ((ref)) without the need to identify the sparse structure in Model ((ref)). This is possible because all biases can be approximated by the data, and the covariance terms in the limiting distributions are of full rank, and they can be estimated from the available data.
In this section, we briefly illustrate the application of the proposed bias-correction method in time series ridge regression. In particular, we consider the following Autoregressive (AR) model:
where ${\mathbf x}_t=(y_{t-1},...,y_{t-p})'$ is the covariate vector consisting of the $p$ lagged variables of $y_t$, and $\boldsymbol{\beta}=(\beta_1,...,\beta_p)'$ is the associated parameter vector. Note that $E(y_t) = 0$ for simplicity. We restrict our consideration to cases where $p<n$. It's uncommon to encounter situations where the number of lagged regressors in an AR model exceeds the sample size. As a matter of fact, it is rare to see that the order of an AR model exceeds 20 for non-seasonal time series in empirical applications. See, for instance, Woodward2012.
Let $\widetilde{\mathbf y}=(y_{p+1},...,y_n)'$, $\widetilde{\mathbf X}=({\mathbf x}_{p+1},...,{\mathbf x}_{n})'$, and $\widetilde\mbox{\boldmath$\varepsilon$}=({\varepsilon}_{p+1},...,{\varepsilon}_{p})'$, then, Model ((ref)) can be expressed as
where we assume $\mbox{Cov}(\widetilde\mbox{\boldmath$\varepsilon$})=\sigma^2{\mathbf I}_{n-p}$ for simplicity. For a proper $\lambda>0$, the de-biased estimator in ((ref)) can be written as \[\widehat\boldsymbol{\beta}_{c,k}(\lambda) =\boldsymbol{\beta}-\lambda^{k+1}(\widetilde{\mathbf X}'\widetilde{\mathbf X}+\lambda{\mathbf I}_p)^{-(k+1)}\boldsymbol{\beta}+\sum_{j=0}^k\lambda^j(\widetilde{\mathbf X}'\widetilde{\mathbf X}+\lambda{\mathbf I}_p)^{-(j+1)}\widetilde{\mathbf X}'\widetilde\mbox{\boldmath$\varepsilon$}.\] For weakly stationary time series sequence $\{y_t\}$, the ergodic theorem guarantees that \[\frac{1}{n}\widetilde{\mathbf X}'\widetilde{\mathbf X}=\frac{1}{n}\sum_{t=p+1}^n{\mathbf x}_t{\mathbf x}_{t}'\rightarrow_p E({\mathbf x}_t{\mathbf x}_t'),\,\,\text{as}\,\, n\rightarrow\infty.\] Suppose the covariance matrix $E({\mathbf x}_t{\mathbf x}_t')$ admits a similar spectral decomposition as ${\mathbf X}'{\mathbf X}/n$ in Assumption (ref) above, and the the magnitude of the penalty $\lambda$ satisfies $\lambda\asymp n$ as that in Assumption (ref), we can show that \[\sum_{j=0}^k\lambda^{j+1}(\widetilde{\mathbf X}'\widetilde{\mathbf X}+\lambda{\mathbf I}_p)^{-(j+1)}=O_p(1),\,\,\text{and}\quad\frac{1}{\lambda}\widetilde{\mathbf X}'\widetilde\mbox{\boldmath$\varepsilon$}=O_p(\frac{1}{\sqrt{n}}).\] For large $n$, it follows that
and
Equation ((ref)) suggests that the bias-correction procedure remains applicable to time series regression with weakly dependent regressors. Additionally, Equation ((ref)) confirms that the asymptotic distribution of Theorem (ref)(i) continues to hold for time series data as well.
In this section, we explore methods for constructing confidence and prediction intervals based on ridge regression. For simplicity, we assume $\boldsymbol{\Sigma}_{\varepsilon}=\sigma^2{\mathbf I}_n$ to facilitate the illustration.
We begin by outlining the construction of confidence intervals for the mean response $E(y|{\mathbf x}_0)={\mathbf x}_0'\boldsymbol{\beta}$ with a given covariate ${\mathbf x}_0$ when $p<n$. Note that $\widehat y_0={\mathbf x}_0'\widehat\boldsymbol{\beta}(\lambda)$ on a ridge regression, then, by ((ref)) in Remark (ref)(ii) with a sufficiently large $k$, we have
where \[\boldsymbol{\Sigma}_{1,k}(\lambda)=\sigma^2\sum_{j=0}^k \lambda^j({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-(j+1)}{\mathbf X}'{\mathbf X}\sum_{j=0}^k \lambda^j({\mathbf X}'{\mathbf X}+\lambda{\mathbf I}_p)^{-(j+1)}.\] In practice, we may replace $\sigma$ in $\boldsymbol{\Sigma}_{1,k}(\lambda)$ by $\widehat\sigma$, which is estimated from the data as that in Remark (ref)(iv), and we denote the estimated variance in ((ref)) by $\widehat\boldsymbol{\Sigma}_{1,k}(\lambda)$. Denote $z_{\alpha}$ the critical value of a standard normal distribution such that its tail probability is $\alpha$, then, the $(1-\alpha)$-confidence interval for ${\mathbf x}_0'\boldsymbol{\beta}$ is $[L_1,U_1]$, where
and
Next, we discuss the way to construct prediction intervals for a future value $y_{n+1}$ with a new covariate ${\mathbf x}_{n+1}$. According to Ch. 2 of montgomery2012introduction, by a similar argument as that for the confidence intervals, we can construct the $(1-\alpha)$-prediction intervals of the future value $y_{n+1}$ as $[L_2,U_2]$, where
and
The values presented in ((ref))-((ref)) are directly estimated from the data. Consequently, we are equipped to construct valid confidence and prediction intervals in ridge regression via employing the proposed bias-correction techniques. Analogously, confidence and prediction intervals can also be established for models utilising de-biased ridge estimators post ridge-screening when $p>n$. For instance, the $(1-\alpha)$-prediction interval of the future value $y_{n+1}$ produced by the restricted ridge estimator is $[L_3,U_3]$ with
and
where ${\mathbf x}_{0,k}$ is the $n^*$-dimensional sub-vector of the new covariate ${\mathbf x}_0$ selected by the ridge-screening method.
The bias-variance trade-off is a fundamental concept in machine learning and statistics. It refers to the delicate balance between two sources of error in a predictive model. It is a way of analysing a learning algorithm's expected errors. Therefore, it would be interesting if we can derive the bias-variance trade-off of our de-biased estimators as we increase the number of iterations.
In this section, we study the MSE of the de-biased ridge estimator $\widehat\boldsymbol{\beta}_{c,k}(\lambda)$ in the $k$-th iteration for $p<n$ with a given $\lambda>0$. The analysis of the scenario with $p>n$ is similar, and hence, we only focus on the case of $p<n$. The MSE of the $k$-th de-biased ridge estimator is defined as
For simplicity, we assume $\boldsymbol{\Sigma}_{\varepsilon}=\sigma^2{\mathbf I}_n$. By Assumption (ref) and Theorem (ref)(i),
where ${\mathbf D}_{k}(\lambda)=\mbox{diag}(d_{k,1}(\lambda),...,d_{k,p}(\lambda))$ with \[d_{k,i}(\lambda)=\frac{1}{d_i^2}\left[1-(\frac{\lambda}{\lambda+d_i^2})^{k+1}\right]^2.\] By a similar argument, we can show that
where $\boldsymbol{\Lambda}_{k}(\lambda)=\mbox{diag}(\gamma_{k,1}(\lambda),...,\gamma_{k,p}(\lambda))$ with \[\gamma_{k,i}(\lambda)=[\frac{\lambda}{\lambda+d_i^2}]^{2(k+1)}.\] As the number of iterations $k$ increases, we can readily see that the bias term in ((ref)) decreases while the variance term in ((ref)) increases. Therefore, it is interesting to see whether there is a bias-variance trade-off in the proposed bias-correction procedure.
In the following theorem, we provide some sufficient conditions under which we have a bias-variance trade-off for the proposed method.
From Theorem 2 in theobald1974generalizations, the MSE of the classical ridge estimator can be smaller than that of the least-squares estimator for certain properly chosen $\lambda>0$. In our framework, we observe a similar paradigm concerning the number of iterations $k>0$, as described in the following corollary
We conduct Monte-Carlo experiments to demonstrate the performance of the proposed procedure by considering two scenarios where $p<n$ and $p>n$. Our goal is to study the effect of the bias-correction procedure and the approximation of the asymptotic normalities established in Section (ref).\\
{\bf Example 1.} In this example, we employ the data-generating process in ((ref)) for different settings of $(p,n)$ with $p<n$. For each configuration of $(p,n)$, we set the seed number in R software to 1234 and generate an $n\times p$ matrix ${\mathbf H}$, where its elements are drawn independently from a Uniform distribution $U(-2,2)$. We then perform a singular-value decomposition on ${\mathbf H}$ and obtain its left and right singular vectors ${\mathbf M}\in R^{n\times p}$ and ${\mathbf N}\in R^{p\times p}$, then we let ${\mathbf X}={\mathbf M}{\mathbf N}'$. The first $p/2$ elements of $\boldsymbol{\beta}$ are generated from $U(-2,-1)$ independently, and the remaining $p/2$ elements are from $U(1,2)$. In each replication, the noise $\mbox{\boldmath$\varepsilon$}$ is generated from multivariate normal distribution $N({\bf 0},{\mathbf I}_n)$. We consider the scenarios where $(p,n)=(50,100)$, $(50,400)$, $(100,200)$, and $(100,500)$. The choices of $\lambda$ are $\lambda=0.05n$, $0.1n$, $0.3n$, and $0.5n$. 1000 replications are used in each setting throughout the experiments.
We report the estimation errors of the ridge estimators and the de-biased ones in Table (ref). For each configuration of $(p,n,\lambda)$, we define the empirical mean squared errors (MSEs) as
and
where $\widehat\boldsymbol{\beta}^{(j)}(\lambda)$ and $\widehat\boldsymbol{\beta}_{c,k}^{(j)}(\lambda)$ are the ridge estimators and the de-biased one with $k$ iterations in the $j$-th replication, respectively. The standard errors $\widehat\sigma_0$ and $\widehat\sigma_{k}$ are estimated via
Thus, we evaluate the estimations of the standard errors by the ridge estimators and the de-biased ones.
For each setting in Table (ref), we see that the MSE first decreases as the number of iterations increases, but it increases if we use too many iterations to correct the bias terms. This is understandable because it is in line with the bias-variance trade-off in the machine learning literature. See, for example, hastie2009elements. Specifically, we plot the MSEs for $(p,n)=(50,100)$ and $(100,200)$ with $\lambda=0.05n$ in Figure (ref), where we can clearly see that the MSEs can achieve a minimum point as we increase the number of iterations, or equivalently, the model complexity, and the MSE will increase if we try to completely correct the biases. Finally, the MSE becomes stable, which is consistent with our asymptotic theory. Furthermore, comparing $\widehat{\sigma}$ and $\widehat{\sigma}_{100}$, we see that the standard error estimated using the debiased ridge estimator is closer to the true one (unity), showing the effectiveness of our bias-correction procedure.
Next, we study the consistency of the de-biased estimators using the proposed method. For simplicity, we investigate the average estimation errors (AEEs) of the ridge and de-biased estimators. These AEE are defined, respectively, as
and
based on the 1000 simulations. These measures, where we use $\ell_2$-norm to quantify the empirical biases. The AEEs are reported in Table (ref). From the table, we can see a decreasing pattern for each $(p,n)$ as the number of iterations increases, implying that our de-biased estimators are consistent to the true parameters for sufficiently large $k$.
Finally, we investigate the performance of statistical inference using the de-biased estimators given in Theorem (ref) for $p<n$. For simplicity, we only consider the case when $(p,n)=(100,200)$ and $\lambda=0.3n$, and can produce similar results for other settings. Let ${\mathbf e}_i$ be the $i$-th standard unit vector of $R^p$, \[\boldsymbol{\theta}_1=(0.8,-1,0.5,{\bf 0}_{p-3}')',\,\,\text{and}\,\,\boldsymbol{\theta}_2=(-1,0.5,0.8,{\bf 0}_{p-3}')',\] where ${\bf 0}_s$ is an $s$-dimensional vector of zeros. Figure (ref) presents the histograms of $\sqrt{n}{\mathbf e}_1'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$, $\sqrt{n}{\mathbf e}_2'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$, $\sqrt{n}\boldsymbol{\theta}_1'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$, and $\sqrt{n}\boldsymbol{\theta}_2'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$ obtained from 1000 experiments. From Figure (ref), we see that the biases of the empirical means in all the histogram plots are significantly large and diverge from zero, implying that the traditional ridge estimators are not appropriate for making statistical inference as they overlook the bias terms. In contrast, we plot the empirical histograms of the de-biased estimators in Figure (ref) under the same setting. From Figure (ref), we see that the finite sample performance is quite satisfactory, and the bias effect of the de-biased estimators is small. In addition, the curve of the normal distribution in Theorem (ref)(i) is added to the corresponding histograms, where the standard error term $\sigma$ is estimated from the data as that in Remark (ref)(iv) and Table (ref). From these curves, we can further confirm that the proposed inference method is valid for the de-biased ridge estimators, indicating that one can make valid inferences using the ridge estimators with our bias-correction method in practice.\\
{\bf Example 2.} In this example, we consider the scenario of $p>n$. The seed number used in this example is the same as that in Example 1. Each row of the design matrix ${\mathbf X}$ is generated independently from multivariate standard normal distribution $N({\bf 0},{\mathbf I}_p)$. For each configuration of $(p,n)$, the first $5$ elements in $\boldsymbol{\beta}$ are generated from $U(-5,-2)$, the $6$th to the $10$th elements are from $U(2,5)$, and the remaining ones are zeros. In other words, we consider a sparse vector $\boldsymbol{\beta}$ with $p^*=10$ in ((ref)). The dimensions and the sample sizes are $(p,n)=(150,120), (150,140), (220,180)$, and $(220,200)$, and we consider $\lambda=0.1n$, $0.3n$, and $0.8n$ in each replication for any given $(p,n)$. For simplicity, we choose $n^*=40$ in the ridge screening of ((ref)) with $k=100$ iterations. $1000$ replications are used for each configuration of $(p,n)$ throughout the experiments.
In Table (ref), we report the empirical MSEs of the ridge estimators, the de-biased estimators before applying the ridge-screening method, and those after the ridge-screening with $k=100$ iterations in the first stage of bias-corrections. Specifically, the upper panel in Table (ref) for each $\lambda$ presents the MSEs of the ridge estimators and the de-biased estimators with 100 iterations before applying the ridge-screening approach. We see that the MSEs of the de-biased estimators become quite stable as we increase the number of iterations. This result is consistent with our findings that there exists some minor bias that cannot be corrected from the data. In addition, the estimated standard errors are not close to unity, which is the true value in the experiments. In the lower panel of Table (ref) for each $\lambda$, we apply the ridge-screening approach with $k=100$ iterations and sort out the largest $n^*=40$ components of $\widehat\boldsymbol{\beta}_{c,100}(\lambda)$ to obtain the restricted ridge estimators. We further apply the bias-correction approach, and the de-biased estimators are obtained in another $100$ iterations. From the lower panels of Table (ref) for each $\lambda$, we see that the bias of the de-biased estimators decreases sharply to a relatively stable value for each configuration of $(p,n)$, which is understandable since the selected model consists of more parameters than those in the true one. In addition, the estimation errors are significantly smaller than those without the ridge-screening approach, and the estimated standard errors after bias correction are closer to the true one than those of the methods without using the ridge-screening and the bias-correction. We also note that the estimated standard errors are still larger than the true one (unity) when $\lambda=0.8n$ after the ridge-screening, but they are much closer to the true ones than those produced by the original ridge estimators. In addition, the result can be optimised by choosing more appropriate tuning parameter $\lambda$ (e.g. $\lambda=0.1n$ or $0.3n$) such that the estimated standard errors are close to one. Furthermore, we also note that there is also a bias-variance trade-off in the MSEs of the de-biased estimators with or without the ridge-screening. Figure (ref) in Sec. (ref) of the Supplement plots the the empirical MSEs of the de-biased ridge estimators for different number of iterations for $\lambda=0.5n$ and $(p,n)=(150,120)$, where Part (a) plots the MSEs of the de-baised ridge estimators before applying the ridge-screening method, and Part (b) provides the MSEs of the de-biased ridge estimators after the ridge-screening. From the plots, it is clear that there is an optimal number of iteration that minimises the MSE, which is in line with our asymptotic results in Theorem (ref).
We also study the performance of the proposed ridge-screening method in Table (ref) of the Supplement, where the empirical probability (EP) is calculated by
where $\mathcal{M}_{k}^i(\lambda^*)$ is the recovery in the $i$-th experiment and $|\mathcal{M}_0\cap \mathcal{M}_{k}^i(\lambda^*)|$ is the cardinality of the set $\mathcal{M}_0\cap \mathcal{M}_{k}^i(\lambda^*)$. From Table (ref), the proposed method provides satisfactory performance in variable selections. The comparison results with other variable selection approaches, such as the LASSO, are not provided in this experiment due to the empirical probability of correct recoveries by the proposed RS method being $100\%$ in all settings of Table (ref).
Similarly to the experiments in Example 1, we present the empirical average estimation errors (AEEs) in Table (ref). The settings of $(p,n)$, $\lambda^*$, $\lambda$, $k$, and $l$ are the same as those in Table (ref). From Table (ref), we see that the average estimation errors before ridge-screening are also quite stable as we increase the number of iterations, which is similar to the findings in Table (ref). After we apply the ridge-screening approach to the de-biased estimators in the upper panel of Table (ref) with $k=100$ for each $\lambda^*$, the bias terms (in the lower panel for each $\lambda$) decrease substantially from their previous values, and the biases also decrease to stable values as we increase the number of iterations in the second-stage bias correction. Notably, the bias terms are all significantly smaller than those in the upper panel without the ridge-screening procedure. This finding is in agreement with our theoretical results in Theorem (ref) and Theorem (ref).
Finally, we study the performance of the inference method in Theorem (ref)(ii). We only consider the case when $(p,n)=(220,200)$, $\lambda^*=\lambda=0.1n$, and $n^*=40$, and we can produce similar results for other settings. Define \[\boldsymbol{\gamma}_1=(0.8,-1,0.5,{\bf 0}_{p-3}')',\quad\boldsymbol{\gamma}_2=(0,-1,0.5,0.8,{\bf 0}_{p-4}')',\] \[\boldsymbol{\gamma}_3=({\bf 0}_3',-0.9,0.4,-0.8,{\bf 0}_{p-6}')',\,\,\text{and}\,\,\boldsymbol{\gamma}_4=({\bf 0}_4',0.5,0.7,-0.8,{\bf 0}_{p-7}')'.\] We plot the empirical histograms of $\sqrt{n}\boldsymbol{\gamma}_1'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$, $\sqrt{n}\boldsymbol{\gamma}_2'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$, $\sqrt{n}\boldsymbol{\gamma}_3'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$, and $\sqrt{n}\boldsymbol{\gamma}_4'(\widehat\boldsymbol{\beta}(\lambda)-\boldsymbol{\beta})$ in Figure (ref), from which it is observed that most of the estimators are not centred around zero. Subsequently, we apply the bias-correction procedure to the original ridge estimators, and the empirical histograms of the bias-corrected estimators are presented in Figure (ref) as compared to Figure (ref). From Figure (ref), it can be seen that while some estimators become more centred around zero (e.g., (a) and (b)) due to the bias-correction procedure, the validity of inference without ridge-screening cannot be assured, as (c) and (d) still exhibit significant biases.
Furthermore, we apply the ridge-screening approach to the de-biased estimators depicted in Figure (ref), and subsequently obtain restricted ridge estimators along with their de-biased counterparts. Figure (ref) showcases the empirical histograms of: (a) $\sqrt{n}\boldsymbol{\gamma}_1'(\widehat\boldsymbol{\beta}_{\mathcal{M}_k,l}^*(\lambda)-\boldsymbol{\beta})$; (b) $\sqrt{n}\boldsymbol{\gamma}_2'(\widehat\boldsymbol{\beta}_{\mathcal{M}_k,l}^*(\lambda)-\boldsymbol{\beta})$; (c) $\sqrt{n}\boldsymbol{\gamma}_3'(\widehat\boldsymbol{\beta}_{\mathcal{M}_k,l}^*(\lambda)-\boldsymbol{\beta})$; and (d) $\sqrt{n}\boldsymbol{\gamma}_4'(\widehat\boldsymbol{\beta}_{\mathcal{M}_k,l}^*(\lambda)-\boldsymbol{\beta})$, where $\widehat\boldsymbol{\beta}_{\mathcal{M}_k,l}^*(\lambda)$ is a $p$-dimensional vector with $\widehat\boldsymbol{\beta}_{\mathcal{M}_k,l}(\lambda)$ as its sub-vector and the remaining elements are zero. Similar to the curves in Figure (ref), the density curve in Figure (ref) is plotted based on the limiting distribution specified in Theorem (ref)(ii). From Figure (ref), it is evident that the de-biased restricted ridge estimators predominantly cluster around zero. The density curves for the corresponding estimators are quite satisfactory, suggesting the validity of the proposed inference approach.
In this section, we apply the proposed method to forecast the U.S. macroeconomic series and study the out-of-sample prediction intervals using ridge regression. We consider the widely used macroeconomic variables studied by Stock2002, mccracken2016fred, giannone2021economic, and gao2022modeling,gao2023supervised, among many others. The data are obtained from the FRED-MD database which is maintained by St. Louis Federal Reserve Bank. See \url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}. There are 127 variables in the online data set, and we remove 4 of them because of missing values therein. Consequently, we consider the 123 macro variables spanning from July 1962 to December 2019 as all the series have no missing values during this period. The resulting data set is the same as that in gao2023supervised. The detailed variables and transformation codes to ensure the stationarity of each macro variable are provided in Table IA.I of gao2023supervised. The sample size is $n = 690$.
Similarly to Stock2002b, bai2006, and cheng2015forecasting, we adopt the following factor-augmented regression model:
where ${\mathbf X}_t$ consists of 122 macroeconomic variables, $y_{t}$ is the remaining target one, and ${\mathbf f}_t=(f_{1,t},...,f_{r,t})'$ is an $r$-dimensional latent factor process which can be estimated by applying Principal Component Analysis (PCA) to ${\mathbf X}_t$. We set $q=10$ and $r=60$, and apply the ridge regression to estimate the parameters $\alpha_i$ and $\gamma_j$, for $1\leq i\leq q$ and $1\leq j\leq r$, with estimated factors. Hence, the number of covariates in the ridge regression is $p=q+r=70$. We focus on the prediction of the consumer price index: all (CPI-All), which is an important index related to the inflation. See also Stock2002 and gao2022modeling for similar studies.
First, we set $\lambda\in\{0.05n,0.1n,0.2n,...,1.5n\}$, which contains 16 candidates for the penalty parameter. We split the sample into two subsamples, where the first one consists of the first $80\%$ of the data for modelling and the second subsample of the remaining $20\%$ of the data for out-of-sample prediction. We adopt the rolling-window scheme as that in gao2023supervised, that is, we train the factors and the forecasting coefficients using the first subsample to predict the next target data point. Then we repeat the above procedure after moving the next available observation of predictors and target CPI-All variable from the second subsample to the first one to obtain the next prediction. This rolling-window scheme is terminated when there is no more observation to compute the forecasting error. We choose the optimal penalty in the ridge estimator such that the out-of-sample mean-squared forecast errors (MSFEs) achieve a minimum value. Thus, this setting is in line with the advocated models combining regularised estimation and data-driven selection of regularisation parameters in abadie2019choosing. In this experiment, we find that the optimal penalty $\widehat\lambda=0.8n$ for $h=1$ and $\widehat\lambda=0.7n$ for $h=2$.
Next, for each configuration of $(\widehat\lambda,h)=(0.8n,1)$ and $(0.7n,2)$, we execute Algorithm (ref). We observe that the bias-correction procedure terminates within 3 to 5 iterations when setting $\eta=O(10^{-4})$ in equation ((ref)). For simplicity, we opt for $k=10$ in the bias-correction, a choice deemed sufficiently large. Following the guidelines provided in Section (ref), we compute pointwise prediction intervals for the data points in the testing subsample. The standard errors for each prediction are determined using the method described in Remark (ref)(iv). We present the $95\%$ pointwise prediction intervals for both 1-step ahead and 2-step ahead forecasts in Figure (ref) and Figure (ref), respectively. It is evident that the majority of the true observations fall within these prediction intervals. Additionally, we evaluate the coverage rates of the prediction intervals across the 138 points in the testing subsample, and find a consistent coverage rate of $92.03\%$ for both the 1-step and 2-step ahead forecasts, indicating that the $95\%$ prediction intervals perform effectively in this empirical example.
Finally, we investigate the effectiveness of the ridge-screening (RS) approach in out-of-sample predictions. We take the 1-step ahead prediction as an example. For each $\lambda^*\in\{0.05n,0.1n,0.2n,...,1.5n\}$, which is the same as our previous candidate set, we adopt $k=10$ and vary $n^*$ from $10$ to $70$ in ((ref)), that is, we select the largest $10$ to $70$ variables based on their magnitudes from the set $\{|\widehat\boldsymbol{\beta}_{c,10,1}(\lambda^*)|,...,|\widehat\boldsymbol{\beta}_{c,10,70}(\lambda^*)|\}$ and compute their corresponding out-of-sample forecast errors. By minimizing the out-of-sample MSFEs, we determined that the optimal parameters are $(\widehat\lambda^*,n^*)=(0.5n,31)$. This implies that the submodel, consisting of covariates corresponding to the largest 31 elements in $|\widehat\boldsymbol{\beta}_{c,10}(0.5n)|$, yields the smallest out-of-sample MSFEs. Consequently, we chose $\widehat\lambda=\widehat\lambda^*=0.5n$ for the subsequent bias-correction procedure in the restricted ridge estimators, conducting it with $l=10$ iterations. Figure (ref) displays the $95\%$ out-of-sample prediction intervals based on the selected submodel and the corresponding restricted ridge estimators for $h=1$. Interestingly, we observe that some observations not covered by the prediction intervals in Figure (ref) are encompassed by the new intervals in Figure (ref). The coverage rate of these $95\%$ prediction intervals stands at $93.5\%$, surpassing the $92.03\%$ achieved by the previous approach without ridge-screening. This suggests that the ridge-screening method enhances the performance of the prediction intervals.
A fundamental question in ridge regression revolves around the possibility of correcting the bias term using available data, alongside ensuring the validity of statistical inferences without affecting the predictive performance of the original ridge estimator. This paper tackled these long-standing challenges by introducing a straightforward and readily implementable iterative procedure, along with a ridge-screening method. In cases where the dimension $p$ is smaller than the sample size $n$, our procedure effectively corrects the bias term under some mild assumptions. However, when $p>n$ and the projection matrix becomes singular, we propose a ridge-screening approach to eliminate coefficients that are relatively insignificant compared to others. We then specifically focus on the restricted model, which retains a sufficient number of parameters to encapsulate most of the information regarding the response vector. Our bias-correction procedure can be further applied to the restricted model, enabling the correction of the bias. Moreover, the ridge-screening method offers a novel approach to variable selection, which is of independent interest beyond bias correction. We derive the limiting distributions of the de-biased estimators, to facilitate making statistical inferences. Simulated and real data examples are used to corroborate the effectiveness and validity of our proposed methodology. The proposed inference solution for ridge regression serves as an illustrative example addressing the inference challenges of regularised machine learning methods outlined in Section 2.8 of athey2019machine without adversely affecting the predictive performance.
\singlespacing \let\oldbibliography \oldbibliography {3pt}
\setstretch{1.8}
\setcounter{page}{1}