A Simple Approximation to the Distribution of the Ridge Regression Estimator

A Simple Approximation to the Distribution of the Ridge Regression Estimator∗ José Luis Montiel Olea

Ryan Strong Amilcar Velez † Haomin Yu

Zhuoheng Xu

arXiv:2608.02539v1 [econ.EM] 3 Aug 2026

August 4, 2026

Abstract We present a simple Gaussian approximation to the finite-sample distribution of the classical ridge regression estimator. Our approximation captures the fact that, in finite samples, the ridge regression estimator trades off bias and variance to reduce estimation and prediction error. Our approximation is based on nonstandard asymptotics where i) we let the estimator’s regularization parameter grow proportionally to the sample size; and ii) we treat the population regression coefficients as local to the reference vector that defines the estimator’s direction of shrinkage. In contrast to other asymptotic approximations in the literature, we allow for general forms of heteroskedasticity and autocorrelation in the data generating process (at the cost of considering a low-dimensional model where the number of covariates is not allowed to grow with the sample size). We use our simple Gaussian approximation to propose two new strategies to select the regularization parameter for the ridge regression estimator. The suggested strategies select the regularization parameter to minimize either average or worst-case excess prediction risk, where risk is computed using our suggested Gaussian approximation.

Keywords: Heteroskedasticity and autocorrelation, nonstandard asymptotics, prediction risk, Ridge regression, regularization parameter selection.

∗

We would like to thank Bruce Hansen, Lihua Lei, Jörg Stoye, and Dacheng Xiu for very helpful comments and suggestions. Montiel Olea gratefully acknowledges financial support by the National Science Foundation Grant SES-2315600. † All of the authors are affiliated to the Department of Economics at Cornell University. Corresponding author: amilcare@cornell.edu

1

1

Introduction

n We have access to a dataset Dn ≡ {(yi , x⊤ i )}i=1 consisting of n observations of a real-valued

outcome variable, yi ∈ R, and a vector of k covariates, xi ∈ Rk . The dataset Dn is assumed to have been generated by the statistical model yi = x⊤ i β + ϵi ,

n {(x⊤ i , ϵi )}i=1 ∼ P.

(1)

We are willing to restrict the data generating processes under consideration by requiring that (1/n)

n X

p xi x⊤ i → Σ,

n √ X d (1/ n) xi ϵi → Nk (0, Ω),

and

i=1

(2)

i=1

where both Σ and Ω are positive definite matrices, and Nk denotes a multivariate Gaussian distribution of dimension k. These assumptions allow for general forms of heteroskedasticity and autocorrelation in the data generating process.1 We are interested in approximating the finite-sample distribution of the ridge regression estimator: βbλn ≡

n X

!−1 xi x⊤ i + λn Ik

i=1

n X

! xi yi + λn β0

.

(3)

i=1

We refer to the nonnegative scalar λn as the regularization parameter and to β0 ∈ Rk as the reference vector. We index the regularization parameter by the sample size to allow for the possibility that it depends on the dataset Dn . Our first result (Theorem 1) shows that there is a sense in which, under (1)-(2), we can approximate the distribution of βbλ − β by the distribution: n

Nk −(λn /n)(Σ + (λn /n)Ik )−1 (β − β0 ), (1/n)(Σ + (λn /n)Ik )−1 Ω(Σ + (λn /n)Ik )−1 . (4) Our simple approximation captures the fact that, in finite samples, the choice of regularization parameter and reference vector affect the bias and variance of the ridge regression estimator. We view the mean and variance in the approximation in (4) as providing simple generalizations of common formulae for the bias and variance of the ridge estimator derived conditional on covariates and under i.i.d. sampling; see, for example, Hansen (2022, Equations 29.8 and 29.9). The Gaussian approximation in (4) can also be viewed as a generalization of the asymptotic distribution of the ridge estimator presented in the seminal work of Knight and Fu (2000), which they derive under conditional homoskedasticity (in 1

We note that (1) and (2) neither require nor imply that the conditional expectation of the outcome given the covariates is linear. In fact, (1) and (2) can be viewed as the typical high-level assumptions used in the analysis of the linear projection model using the terminology in Hansen (2022, Chapter 2, p. 37).

2

our framework, heteroskedasticity and serial autocorrelation are captured in the matrix Ω). In order to derive Theorem 1, we use nonstandard asymptotics.2 As we will explain later, we assume that the true coefficient β in the model (1) (henceforth denoted βn ) depends on the sample size. We treat βn as local to the reference vector β0 that defines the estimator’s direction of shrinkage. More concretely, we assume that: √

n(βn − β0 ) → b,

(5)

for some b ∈ Rk . In addition to (5), our asymptotic analysis accommodates data-dependent regularization parameters that grow proportionally to the sample size and that satisfy p

λn /n → λ ∈ [0, ∞).

(6)

Our second result (Theorem 2) uses the approximation in (4) to provide an approximation to the excess prediction risk of the ridge regression estimator. As we explain later, excess prediction risk is defined as the additional prediction risk relative to an oracle that knows βn . Predictions of the outcome variable based on the ridge regression estimator take the form b aλn (x) = x⊤ βbλn ,

(7)

where βbλn is defined in (3). The exact finite-sample analysis of the prediction risk of ridge regression is challenging even under a stylized homoskedastic, Gaussian regression model. Several recent papers in the statistics literature provide approximations and lower and upper bound for this prediction risk, including Dobriban and Wager (2018), Hastie, Montanari, Rosset, and Tibshirani (2022), Mourtada and Rosasco (2022), Atanasov, Zavatone-Veth, and Pehlevan (2024). We show that it is possible to use our Theorem 1 to provide an approximation to the excess prediction risk of ridge regression as a function of three types of parameters: λ (the probability limit of λn /n); b (the parameter controlling the local-to-β0 regression coefficient βn ); and the variance parameters Σ and Ω. Importantly, the variance parameters can be consistently estimated, but the local parameter b cannot. Our third result (Theorem 3) presents closed-form, data-driven recommendations for the selection of the regularization parameter λn . Our formulae are designed for the special case in which—under the true unknown data generating process—the covariates are asymptotically uncorrelated and have the same variance; that is, Σ = σx2 Ik , where Ik is the identity matrix of dimension k and σx2 is a common variance parameter. Following Hastie et al. (2022), we refer to this setting as one with isotropic features.3 The key step 2 3

See Powell (2017) for an overview of the role of nonstandard asymptotics in modern econometrics. Although we note that in Hastie et al. (2022) Σ equals the identity matrix, whereas we only require

3

in deriving our closed-form formulae is to choose the regularization parameter to optimize the excess prediction risk approximation in Theorem 2. As we discussed before, this approximation depends on the local parameter b and variance parameters σx2 and Ω. While the latter parameters can typically be consistently estimated under mild assumptions, the local parameter b cannot. We use the common decision-theoretic principles of average and worst-case risk, following Ferguson (1967), to account for the unknown parameter b in the approximate excess prediction risk. If π is a prior distribution on b such that Eπ [b⊤ b] < ∞ and Eπ [b⊤ b] > 0, Theorem 3 shows that the π-optimal selection of the regularization parameter λn in ridge regression is b b∗ = n · trace(Ω) . λ n,π 2 ⊤ σ bx Eπ [b b]

(8)

When the localization parameter b admits a known bound B on its norm (∥b∥ ≤ B), the selection of the regularization parameter λn that minimizes approximate worst-case excess risk (and that we term the minimax choice of λn ) is: b trace(Ω) b∗ . λ n,minimax = n · 2 2 σ bx B

(9)

While the formulae in (8)-(9) pertain to the case of isotropic features, the idea of optimizing approximate excess prediction risk (either by using an average or worst-case criterion) is more general. We explain how to operationalize both of these approaches for nonisotropic features; see Sections 3.3.1 and 3.3.2. In both cases, the optimal choice of regularization parameter is the solution of a nonconvex optimization problem over the positive part of the real line. In this case, we suggest choosing λn by evaluating the objective functions on a grid of candidate regularization parameters. This approach is the same as the one used by conventional statistical packages to find and implement the usual cross-validated choice of regularization parameter. Related Literature: The introduction of the ridge estimator to regression analysis is often credited to Hoerl and Kennard (1970), although the estimator is based on the earlier work of Hoerl (1962); see Hoerl (2020). The textbook version of the ridge regression estimator uses a reference vector β0 = 0k×1 ; see, for example, Friedman, Hastie, and Tibshirani (2017, Equation 3.41, p.63) or Hansen (2022, Chapter 29.5). The more general formulation given in (3) has several precedents in the literature; for example, Swindel (1976, Equation 2.1). See also Anup and Maddala (1984). The high-level assumptions in (2) that we use to analyze the distribution of the ridge estimator are common in the econometrics literature. For example, it is well known that Σ = σx2 Ik .

4

(2) can be verified under i.i.d. sampling and standard primitive conditions on the joint distribution of xi and ϵi ; see Hansen (2022, Assumption 7.2 and Theorems 7.1 and 7.2). However, our framework accommodates richer data structures, allowing departures from i.i.d. sampling as well as general forms of heteroskedasticity and autocorrelation in (x⊤ i , ϵi ). The approximation we propose in (4) has not—to the best of our knowledge—appeared before in the literature expressed at this level of generality. The closest reference that we are aware of is the seminal work of Knight and Fu (2000), who present a Gaussian approximation to the distribution of the ridge regression estimator assuming independence between xi and ϵi (and independence across observations). The distribution in (4) coincides with their result when specialized to the case in which Ω is conditionally homoskedastic (that is, Ω = E[ϵ2i ]Σ) and β0 = 0k×1 . We view the mean and variance in the approximation in (4) as providing simple generalizations of common formulae for the bias and variance of the ridge estimator derived conditional on covariates and under i.i.d. sampling; see, for example, Hansen (2022, Equations 29.8 and 29.9) and also Hoerl and Kennard (1970, Section 4). The idea of considering drifting sequences of parameter values—as those in (5) and (6)— in order to improve the distributional approximations provided by standard asymptotic theory has a long history in econometrics. See Powell (2017) for examples and details on the use of nonstandard asymptotics in econometrics. More concretely, the idea of using nonstandard asymptotics to analyze the distribution of the ridge regression estimator also appears in Knight and Fu (2000) (they refer to it as triangular array asymptotics), and more recently in the work of Shen and Xiu (2025), who assess the predictive performance of several machine learning methods in high-dimensional regressions with low signal-to-noise ratios. Drifting sequences of parameter values are also common in the analysis of shrinkage estimators; see, for example, Hansen (2016, Equation 14, p. 119). Local asymptotics are also used in the work of Hirano and Wright (2017) to compare the performance of different forecasting procedures under model uncertainty. There is also a recent literature analyzing the risk of predictions based on the ridge regression estimator, where the outcome variable predictions take the same form as in (7). For example, Hsu, Kakade, and Zhang (2012), Dobriban and Wager (2018), Hastie et al. (2022), Mourtada and Rosasco (2022), Atanasov et al. (2024). All these papers make distributional assumptions that preclude data generating processes that exhibit at least one of the following features: conditional heteroskedasticity, time-series autocorrelation, or a nonlinear conditional expectation function. It is important to mention, however, that the additional assumptions in these papers allow them to consider high-dimensional approximations where the number of covariates can be large relative to the sample size. In contrast, all the analysis in our paper pertains to a low-dimensional model where we treat the number

5

of covariates as fixed, and let the sample size diverge to infinity. Finally, the idea of using an approximation of the risk function to choose regularization parameters has a long history in statistics and econometrics; see, for example, the discussion in Abadie and Kasy (2019, p. 744). To the best of our knowledge the formulae we provide in (8) and (9) are new to the literature. It is important to mention, however, that one of the initial motivations of our paper was to recover an expression analogous to the optimal regularization parameter for ridge regression given by Hastie et al. (2022) in their Corollary 6 to Theorem 6. Their results show that using high-dimensional asymptotics where k/n → γ, the oracle choice of regularization parameter for the ridge estimator (assuming xi and ϵi are independent and that Σ = Ik ) takes the form λ∗ = nE[ϵ2i ]γ/∥β∥2 . Assuming conditional homoskedasticity, our formulae in (8) and (9) become analogous to Hastie et al. (2022)’s formula, but the true unknown ∥βn ∥2 is replaced by either the prior mean of ∥b∥2 or the maximum value of ∥b∥2 . More generally, we note that while ridge regression and lasso-type estimators with cross-validated tuning parameters are routinely used with data that exhibit heteroskedasticity and autocorrelation, it is not entirely clear what are the theoretical guarantees of cross-validation in these environments; see, for example, the discussion in Kock, Pedersen, and Sørensen (2026). Outline: The rest of this paper is organized as follows. Section 2 introduces the notation, framework, and main assumptions. Section 3 presents our main results. Section 4 presents simulation evidence to illustrate and support our main results. Section 5 concludes.

2

Notation and Framework

n An econometrician has access to a dataset Dn ≡ {(yi , x⊤ i )}i=1 consisting of n observations

of a real-valued outcome variable, yi ∈ R, and a vector of k covariates, xi ∈ Rk . The dataset Dn is assumed to have been generated by the statistical model yi = x⊤ i βn + ϵi ,

n {(x⊤ i , ϵi )}i=1 ∼ Pn .

(10)

The parameters of this statistical model are i) the unknown vector of slope coefficients, ⊤ βn ∈ Rk , and ii) the unknown distribution of the tuple ((x⊤ 1 , ϵ1 ), . . . , (xn , ϵn )), which we

denote as Pn . Note that we have indexed both βn and Pn by the sample size n. While this does not make a difference when conducting finite-sample analysis (since, in that case, n is fixed), it will allow us for more flexibility when considering asymptotic approximations to the distribution of different statistics. We implicitly restrict the parameter space of the statistical model in (10) by requiring Pn to satisfy the following high-level assumption:

6

Assumption 1 (High-level assumptions on Pn ). The distribution Pn satisfies 1. (1/n)

⊤ p i=1 xi xi → Σ for a positive definite matrix Σ.

Pn

√ P d 2. (1/ n) ni=1 xi ϵi → Nk (0, Ω) for a positive definite matrix Ω. The high-level Assumption 1 is a convenient way to consider a large class of data generating processes in our analysis. It is well known that Assumption 1 can be verified under i.i.d. sampling and standard primitive conditions on the joint distribution of xi and ϵi ; see Hansen (2022, Assumption 7.2 and Theorems 7.1 and 7.2). However, our framework accommodates richer data structures, allowing departures from i.i.d. sampling as well as general forms of heteroskedasticity and autocorrelation in (x⊤ i , ϵi ). The continuous mapping theorem and Slutsky’s theorem immediately imply the asymptotic normality of the least-squares estimator of βn in the model (10); namely 
√  d n βbOLS − βn → Nk 0, Σ−1 ΩΣ−1 ,

βbOLS ≡

n X

xi x⊤ i

!−1 n X

i=1

xi yi .

(11)

i=1

This slightly generalizes the normal approximation to the distribution of the least-squares estimator derived under i.i.d. sampling; see, for example, Hansen (2022, Theorem 7.3). In our framework, Ω can be the long-run variance of the process {xi ϵi }∞ i=1 , and βn can vary with the sample size.

3

Main Results

3.1

Asymptotic Distribution of the Ridge Estimator

Define the ridge estimator —with regularization parameter λ ∈ R+ and a reference vector β0 ∈ Rk —to be the solution of the following minimization problem: n

k

i=1

j=1

λX 1X 2 (yi − x⊤ β) + (βj − β0j )2 . i n β∈Rk n

(12)

min

The solution to the minimization problem in (12) can be shown to equal

βbλ ≡

n X

!−1 xi x⊤ i + λIk

i=1

n X

! xi yi + λβ0

.

(13)

i=1

It is known that the objective function that defines the ridge estimator intentionally sacrifices training error (since predictors of yi based on βbλ will have a larger training error than

7

predictions based on the least-squares estimator βbOLS ) by penalizing deviations away from the reference vector β0 . We note that in most textbook definitions of the ridge estimator, the reference vector β0 equals a vector of zeros of dimension k × 1; see, for example, Friedman et al. (2017, Equation 3.41, p.63) or Hansen (2022, Chapter 29.5). Since the usual interpretation of (13) is that the ridge estimator “shrinks” the least-squares coefficients, we allow for the possibility that this shrinkage is towards an arbitrary reference vector β0 that could be different from zero. We note that this more general formulation has several precedents in the literature; for example, Swindel (1976, Equation (2.1)). See also Anup and Maddala (1984). As we mentioned in the introduction, the main goal of this paper is to present a useful approximation to the finite-sample distribution of the ridge estimator in (13). We are particularly interested in deriving an approximation that captures the well-known fact that, in finite samples, the ridge estimator presents a bias-variance trade-off. Textbook expressions of the bias and variance of the ridge estimator in finite samples are usually derived assuming i.i.d. sampling and conditional mean independence; i.e., EPn [ϵi |xi ] = EPn [ϵi ] = 0. See, for example, Hansen (2022, Chapter 29.6). This stands in contrast with the generality under which the asymptotic distribution of the least-squares estimator can be derived; see Equation (11) above. Our first result shows that when the regularization parameter λ is negligible relative to the sample size (even if it is data dependent), the standard asymptotic distribution of the ridge estimator will not be useful to capture the bias-variance trade-off that is present in finite samples. More precisely, the following result shows that the asymptotic distribution of the ridge estimator will be asymptotically equivalent to that of the least-squares estimator. In order to state our result (and to allow for the possibility that the regularization parameter λ could be selected in a data-driven manner), we consider a possibly stochastic sequence of regularization parameters {λn }∞ n=1 . Proposition 1 (Approximation to the distribution of βbλn when λn /n is negligible). Suppose that the dataset Dn was generated according to the statistical model (10), and suppose that Pn satisfies Assumption 1. If √ p p λn /n → 0 and (λn /n) n(βn − β0 ) → 0, then

 √   p √  n βbλn − βn − n βbOLS − βn → 0.

Proof. See Appendix A.1.

8

The intuition behind this result follows directly from the expression for βbλn in (13). When n is large and λn /n is close to zero, then βbλn is approximately the same as βbOLS . √ This is why the distribution of n(βbλn − βn ) is approximately the same as the distribution of the least-squares estimator. The proof of Proposition 1 shows that if λn is not negligible relative to the sample size, then the asymptotic distribution of the ridge estimator will differ from that of the leastsquares estimator. Moreover, the asymptotic distribution will capture the bias-variance trade-off that the ridge estimator faces in finite samples. Theorem 1 (Approximation to the distribution of βbλn when λn /n is potentially non-negligible). Suppose that the dataset Dn was generated according to the statistical model (10), and suppose that Pn satisfies Assumption 1. If p

λn /n → λ ∈ [0, ∞) and then

√

n(βn − β0 ) → b ∈ Rk ,


√  d n βbλn − βn → Nk −λ(Σ + λIk )−1 b, (Σ + λIk )−1 Ω(Σ + λIk )−1 .

(14)

(15)

Proof. See Appendix A.2. We view the mean and variance in the limiting distribution in (15) as providing simple generalizations of common formulae for the bias and variance of the ridge estimator derived conditional on covariates and under i.i.d. sampling; see, for example, Hansen (2022, Equations (29.8) and (29.9)). The key assumptions of Theorem 1 are the two conditions in (14). The first condition requires the regularization parameter to be potentially non-negligible relative to the sample size. Note that we allow for the possibility that the regularization parameter is selected in a data-driven way, as long as it has a deterministic probability limit. We note that the idea of using sequences of regularization parameters λn to analyze the asymptotic properties of penalized estimators is not ours and has several precedents in the statistics literature; see, for example, Knight and Fu (2000, Theorems 2 and 3) where different rate conditions on λn are used to analyze the ridge estimator and more general lasso-type estimators. The second condition in (14) can be interpreted as saying that the true regression coefficient βn is “local-to-β0 ”. This means that we are assuming that the reference vector used to compute the ridge estimator in (13) is not too far from the true (and unknown) βn . √ We note that if n(βn −β0 ) is unbounded, then any non-negligible regularization parameter  √ b λn leads to a bias in n βλn − βn that becomes arbitrarily large. The idea of considering drifting sequences of parameter values in order to improve the distributional approximations provided by standard asymptotic theory has a long history 9

in econometrics. See, for example, the local-to-zero asymptotics in the linear instrumental variables model of Staiger and Stock (1997); the local-to-unit-root asymptotics in the study of nearly integrated autoregressive processes in Phillips (1988) and its recent generalization in Dou and Müller (2021); the local-to-identification-failure analysis in the study of nonlinear Generalized Method of Moments models with weak identification in Andrews and Mikusheva (2022); and the growing-number-of-folds asymptotic framework in Velez (2024), which improves distributional approximations for Debiased Machine Learning estimators. See Powell (2017) for more examples and details on the use of nonstandard asymptotics in econometrics. Some recent work more directly related to our set-up is the paper of Shen and Xiu (2025), studying a high-dimensional linear regression model with weak signals. While our work focuses on regression models where the dimension of the covariates is fixed, a version of the weak signals model in Shen and Xiu (2025) can be obtained in the case in which the reference vector equals zero and the true parameter βn drifts towards zero at √ rate 1/ n. Although Theorem 1 follows from elementary asymptotic theory, to our knowledge, we are not aware of any previous work that formally derives the simple approximation in (15) at the level of generality allowed by Assumption 1. Textbook results for the bias and variance of the ridge estimator (usually derived conditional on covariates and under i.i.d. sampling of a linear regression model) can yield a Gaussian distribution similar to (15) in finite samples if the error term in the regression model is also assumed to be Gaussian (conditional on the vector of covariates). But the analysis is more nuanced if the error distribution is non-Gaussian and/or there is serial autocorrelation in the regression residuals. The closest result to (15) that we have found in the literature appears in the seminal paper by Knight and Fu (2000), which provides asymptotic results for lasso-type estimators (including the ridge estimator). On p. 1368 they use their Theorem 4 to derive an approximation to the asymptotic distribution of a ridge estimator with reference vector β0 = 0k×1 . Under i.i.d. sampling and a conditional homoskedasticity assumption (i.e., 2 ⊤ 2 Ω = E[xi x⊤ i ϵi ] = E[xi xi ]E[ϵi ]), they obtain the same formula as in (15) but specialized

to the case in which Ω = E[ϵ2i ]Σ. Thus, our results can be viewed as a generalization of their analysis of the asymptotic distribution of the ridge estimator, but allowing for more general data structures (where heteroskedasticity and serial autocorrelation are permitted and captured by the matrix Ω), and where we also allow for a general reference vector β0 . Knight and Fu (2000) also consider an asymptotic sequence where βn changes with the sample size (they refer to this as triangular array asymptotics). We note that a potentially interesting extension of our results is to derive asymptotic approximations for lasso-type estimators (analogous to part b) of their Theorem 4) under our Assumption 1. These results could then be used to provide concrete recommendations for the choice of regularization

10

parameter. In the remaining part of the paper, we maintain our focus on the ridge estimator and use the asymptotic approximation in Theorem 1 to present approximations to the prediction risk of the ridge estimator, and we use statistical decision theory to provide recommendations on the choice of regularization parameter.

3.2

An Approximation to the Prediction Risk of Ridge Regression

In this section we use Theorem 1 to provide an approximation to the prediction risk of ridge regression. We start by providing a definition of the prediction problem we are interested in, and we then present a formal definition of the risk of predictions based on ridge regression. n Prediction Problem. An econometrician has access to a dataset Dn ≡ {(yi , x⊤ i )}i=1

consisting of n observations of a real-valued outcome variable, yi ∈ R, and a vector of k covariates, xi ∈ Rk . The dataset Dn is assumed to have been generated by the statistical model (10) with parameters (βn , Pn ) that satisfy Assumption 1. A prediction function is a mapping a : Rk → R that maps covariates into a predicted value for the outcome variable. We note that a prediction function is defined for all possible values of covariates, including those that have not been observed in the sample. In order to define a loss function for the prediction problem, we introduce an additional assumption. ∞ Assumption 2. The stochastic process {(x⊤ i , ϵi )}i=1 is strictly stationary in the sense of

Definition 1.3.3 in Brockwell and Davis (2013). The strict stationarity assumption means that—associated to Pn —there exists a probn ability distribution P over Rk+1 such that, for every i = 1, . . . , n, we have (x⊤ i , ϵi ) ∼ P.

P

We will refer to P as the stationary distribution of Pn . We use this stationary distribution to define a loss function for the prediction problem as follows. Suppose we sample a new pair (x⊤ , ϵ) according to the stationary distribution P, and use the slope coefficient βn to construct a new outcome-covariate pair (y, x) = (x⊤ βn + ϵ, x).4 We can then define the loss associated to a prediction function a : Rk → R as L(a, βn , P) ≡ E(y,x)∼(βn ,P) [(y − a(x))2 ]. In a slight abuse of notation, (y, x) ∼ (βn , P) is used to capture the fact that the expectation is taken over outcome-covariate pairs generated using (βn , P) as described above. 4

Note then that in our analysis we do not allow the distribution used to evaluate prediction error to deviate arbitrarily from the distribution that generated the data. As shown by Patil, Du, and Tibshirani (2024), such an assumption can have important implications over the optimal choice of regularization parameter. A framework for analyzing the distributionally robust prediction error for the square-root lasso and related estimators (and for optimally selecting tuning parameters) has been proposed in Montiel Olea, Rush, Velez, and Wiesel (2026), but their framework excludes the ridge regression estimator.

11

Prediction Risk. In a prediction problem, the goal is to construct a data-driven prediction function; that is, we want to use the available training data to choose a prediction function. In a slight abuse of notation, we denote data-driven prediction functions as b a. Formally, we think of b a as a mapping that takes the data set Dn as input and returns a prediction function b a from covariates to outcomes. Because b a is constructed from the training data, it is a random function. The prediction risk is then the expected prediction loss R(â; βn , Pn ) = E(βn ,Pn ) [L(b a; βn , P)],

(16)

where E(βn ,Pn ) means that the expectation is taken over the distribution of Dn , which is parameterized by (βn , Pn ). Approximate Prediction Risk. Consider the prediction function based on the ridge regression estimator b aλn (x) = x⊤ βbλn ,

(17)

where βbλn is defined in (13) with λ replaced by λn . The exact finite-sample analysis of the prediction risk of ridge regression is challenging even under the stylized homoskedastic, Gaussian regression model. However, there are different papers that provide approximations (and lower/upper bounds) to this prediction risk. For example, Dobriban and Wager (2018), Hastie et al. (2022), Mourtada and Rosasco (2022), Atanasov et al. (2024). The following result shows that we can use Theorem 1 to provide an approximation to the excess prediction risk of ridge regression as a function of three types of parameters: λ (the probability limit of λn /n); b (the parameter controlling the local-to-β0 regression coefficient βn ); and the variance parameters Σ and Ω defined in Assumption 1. More precisely, define the approximate excess risk function Rne (λ; b, Σ, Ω) ≡

1 2 ⊤ 1 λ b (Σ + λIk )−1 Σ(Σ + λIk )−1 b + trace (Σ + λIk )−1 Ω(Σ + λIk )−1 Σ . n n (18)

Theorem 2 (Approximation to the excess prediction risk of ridge regression). Suppose that the dataset Dn was generated according to the statistical model (10), that Pn satisfies Assumptions 1 and 2, and that the stationary distribution P satisfies EP [xi ϵi ] = 0k×1 and EP [xi x⊤ i ] = Σ. If p

(i) λn /n → λ ∈ [0, ∞) and

√

n(βn − β0 ) → b ∈ Rk ,

b Ω) b that are consistent for (Σ, Ω), (ii) there exist estimators (Σ, (iii) there exists δ > 0 such that h i sup EPn (Zn⊤ ΣZn )1+δ < ∞, n

12

Zn ≡

√

n(βbλn − βn ),

then b Ω) b + o(β ,P ) (1/n), R(b aλn ; βn , Pn ) − σ 2 = Rne (λ; b, Σ, n n

(19)

where σ 2 ≡ EP [ϵ2 ]. Proof. See Appendix A.3 The difference between the finite sample prediction risk of a given data-driven predictor b a—which we have denoted by R(b an ; βn , Pn )—and the residual variance σ 2 is typically referred to as the excess prediction risk or simply excess risk.5 Theorem 2 says that, in large samples, the excess prediction risk of ridge regression equals—up to a term that is small in probability—by the approximate excess prediction risk function defined in (18), evaluated b Ω). b The result shows that if we fix (Σ, b Ω), b the excess risk is expected to vary at (λ, b, Σ, as a function of i) how large is the penalty parameter relative to the sample size (λ is the probability limit of λn /n) and ii) how close is the reference vector β0 to the true coefficient √ βn (b is the limit of n(βn − β0 )). We note that part ii) of the assumptions of Theorem 2 can easily be verified in lowdimensional environments, where k is fixed and n grows to infinity. However, it is well known that—even if conditional homoskedasticity holds—the consistent estimation of Ω = σϵ2 Σ is a nontrivial problem when k grows proportionally to n; see, for example, Liu, Zheng, and Feng (2020). The estimation of Ω in high-dimensional problems where heteroskedasticity is present is even more challenging; see Cattaneo, Jansson, and Newey (2018) and the references therein. Theorem 2 is conceptually related to an active area of research in the statistics literature providing asymptotic and nonasymptotic approximations to the excess prediction risk of the ridge regression estimator; see, for example, Dobriban and Wager (2018),Hastie et al. (2022), Mourtada and Rosasco (2022) and the references therein. All these papers make distributional assumptions that preclude data generating processes that exhibit conditional heteroskedasticity, time-series autocorrelation, or a nonlinear conditional expectation function. The additional assumptions in these papers, however, allow them to consider highdimensional asymptotics where the number of covariates is allowed to grow relative to the sample size. In contrast, all the analysis in this paper pertains a low-dimensional model where we treat the number of covariates as fixed, and let the sample size diverge to infinity.

3.3

Selection of λ

In this section, we use the approximation to the excess prediction risk of ridge regression given in Theorem 2 to make a concrete recommendation regarding the selection of the 5

See, for example, Mourtada (2022).

13

regularization parameter λn . The key insight is to pick λ—which is the probability limit of b Ω), b which is the leading term in (19). The main λn /n—to optimize the function Re (λ; b, Σ, n

conceptual challenge is that approximation depends on the unknown localization parameter √ b ∈ Rk , which was defined to be the limit of the sequence n(βn − β0 ). Since b cannot be estimated consistently, we suggest handling this localization parameter in a way that we think is entirely analogous to the evaluation of risk functions in statistical decision theory: optimizing either the worst-case or the average risk. Our suggested method provides a b and Ω) b well-defined mapping from the data (through the estimated covariance matrices Σ to values of the regularization parameter. We provide details below. 3.3.1

Choosing λ to minimize the worst-case approximate excess risk

We first consider the problem of minimizing the worst-case approximate excess risk. Let B ⊆ Rk be a user-specified set of potential values for b. For example, one could take B ≡ {b ∈ Rk : ∥b∥ ≤ B}. for some known B > 0, where ∥ · ∥ denotes the standard Euclidean b∗ is a B-minimax selection of the regularization parameter in ridge norm. We say that λ n

b∗ = n · λ b∗ where λ b∗ solves the minimax problem regression if λ n b Ω). b inf sup Rne (λ; b, Σ,

(20)

λ≥0 b∈B

We make two comments about the minimax problem in (20). First, if the set B is unbounded, bn = 0. then the B-minimax choice of the regularization parameter in ridge regression is λ bn under the minimax criterion we This means that in order to obtain a different choice of λ will require a bound on the parameter b ∈ Rk . Second, while the outer optimization problem in (20) could be handled via grid search bn is chosen by means of cross-validation), the inner optimization (as it is done whenever λ problem requires the evaluation of the worst-case approximate excess risk. Since only the b Ω) b in (18) depends on b, we can solve for the first term in the expression of Rne (λ; b, Σ, worst-case risk by solving the following constrained optimization problem: b + λIk )−1 Σ( b Σ b + λIk )−1 b vbworst-case (λ) ≡ sup b⊤ (Σ

subject to b ∈ B.

(21)

b

Algebra shows that when B ≡ {b ∈ Rk : ∥b∥ ≤ B}—for some known B > 0, and with ∥ · ∥ given by the Euclidean norm—the constrained optimization problem in (21) can be solved in closed-form up to a maximum eigenvalue computation (denoted by µmax (·)):   b + λIk )−1 Σ( b Σ b + λIk )−1 . vbworst-case (λ) = B 2 · µmax (Σ

14

Thus, in the special case in which the Euclidean norm of b is bounded by B, the B-minimax bn can be implemented as λ b∗ = n × λ b∗ where λ b∗ solves the problem choice of λ n

   b . (22) b Σ b + λIk )−1 Σ b Σ b + λIk )−1 + trace (Σ b + λIk )−1 Ω( b + λIk )−1 Σ( inf λ2 B 2 µmax (Σ 

λ∈R+

The minimization problem in (22) is a nonlinear (and nonconvex) optimization problem over R+ and, in general, does not have a closed-form solution. Our suggestion is to minimize this function over the same grid of parameter values used in the standard implementation of (leave-one-out) cross-validation for ridge regression provided in standard statistical packages such as glmnet. 3.3.2

Choosing λ to minimize average approximate excess risk

We now consider the problem of choosing the regularization parameter of ridge regression by minimizing the average approximate excess risk. Let π denote a user-specified probability b∗ is a π-optimal selection of the regularization parameter distribution over Rk . We say that λ n

b∗ = n · λ b∗ where λ b∗ solves the problem in ridge regression if λ n b Ω)]. b inf Eb∼π [Rne (λ; b, Σ,

λ≥0

(23)

b Ω) b in (18) depends on Once again, since only the first term in the expression of Rne (λ; b, Σ, b, we can solve for the average excess risk by evaluating the following expectation: i h b Σ b + λIk )−1 b . b + λIk )−1 Σ( vbaverage (λ) ≡ Eb∼π b⊤ (Σ

(24)

This expectation can be evaluated analytically in some cases. For example, suppose b ∼
Nk 0, C 2 Ik , where C is a user-specified hyper-parameter. Algebra shows that   b + λIk )−1 Σ( b Σ b + λIk )−1 . vbaverage (λ) = C 2 · trace (Σ A similar formula can be derived if the distribution over b is elliptical with mean zero and a covariance matrix proportional to the identity matrix; see Appendix B.1. The π-optimal bn in these cases is λ b∗ = n × λ b∗ , where λ b∗ solves the problem choice of λ n

    b + λIk )−1 Σ( b Σ b + λIk )−1 + trace (Σ b + λIk )−1 Ω( b Σ b + λIk )−1 Σ b . (25) inf λ2 C 2 trace (Σ

λ∈R+

Just as for the B-minimax choice of regularization parameter, our suggestion is to minimize this function over the same grid of parameter values used in the standard implementation of (leave-one-out) cross-validation for ridge regression provided in standard statistical

15

packages. 3.3.3

Isotropic Features

We now specialize the derivations in Sections 3.3.1 and 3.3.2 to the case in which Σ— the matrix of second moments of the vector of covariates—is known to be proportional to identity matrix; that is, Σ = σx2 Ik . Hastie et al. (2022) refer to predictions problems of this form as prediction with isotropic features.6 Algebra shows that with isotropic features, both the minimax and average-risk objectives outlined in the previous subsections simplify substantially, allowing explicit characterization of the optimal regularization parameter for ridge regression. The key observation is that under isotropic features, the approximate excess risk in (18) simplifies to Rne (λ; b, σx2 , Ω) =

1 λ2 σx2 1 σx2 ⊤ b b + trace(Ω). n (σx2 + λ)2 n (σx2 + λ)2

(26)

The following theorem presents the optimal selection of the regularization parameter for ridge regression under the minimax and average risk criterion. Theorem 3 (Optimal choice of regularization parameter with isotropic features). Suppose that the assumptions of Theorem 2 hold. Assume that Σ = σx2 Ik , σx2 > 0, and that there b that are consistent for (σx2 , Ω). Then: exist estimators (b σx2 , Ω) (1) If B ≡ {b ∈ Rk

∥b∥ ≤ B}—for some known B > 0, and with ∥ · ∥ given by the

Euclidean norm—the B-minimax selection of the regularization parameter in ridge regression is b trace(Ω) b∗ λ . n,minimax = n · 2 2 σ bx B (2) If Eπ [b⊤ b] < ∞ and Eπ [b⊤ b] > 0, the π-optimal selection of the regularization parameter in ridge regression is b b∗ = n · trace(Ω) . λ n,π σ bx2 Eπ [b⊤ b] Proof. See Appendix A.4 Since the introduction of the ridge regression estimator, the problem of providing a concrete recommendation for its regularization parameter has been analyzed in multiple papers; for example, see Gibbons (1981) and the references therein. The idea of using an approximation of the risk function to choose regularization parameters has a long precedent in statistics and econometrics; see, for example the discussion of Abadie and Kasy (2019) on 6

Hastie et al. (2022) further assume that σx2 = 1.

16

p. 744. To the best of our knowledge the formulae we provide in Theorem 3 are new to the literature. For the formula we provide in (2), the closest reference we were able to find in the literature is the optimal regularization parameter derived in Theorem 2.1 in Dobriban and Wager (2018), p. 254. Their recommendation is based on high-dimensional model where covariates and regression residuals are independent, but covariates need not be isotropic. Assuming a prior on regression coefficients such that Eπ [β] = 0 and Varπ [β] = α2 Ik /k (see their Assumption RRC), Dobriban and Wager (2018) recommend a regularization parameter of the form λ̂ = n ·

k 1 k = 2, n α2 α

which assumes σϵ2 = 1. In the homoskedastic case with isotropic covariates—Ω = σϵ2 σx2 Ik — our formula matches their recommendation since b k k b∗ = n · trace(Ω) = = 2, λ n,π 2 ⊤ ⊤ α σ bx Eπ [b b] Eπ [βn βn ] √ where we have used the fact that the reference vector β0 is zero, that βn = b/ n, and that σϵ2 = 1. However, our result allows for heteroskedasticity and autocorrelation in the data generating process by adjusting the regularization parameter as a function of Ω. The main limitation of our result, as discussed before, is that our theory is only applicable to low-dimensional models. As we discussed after the statement of Theorem 2, it is not entirely clear how to construct a good estimator for Ω in high-dimensional models. While we think that the extension of our main theorems to high-dimensional environments could significantly extend the relevance and applicability of our results, such an extension is outside the scope of this paper (and will likely require different techniques). More generally, it is important to mention again that one of the initial motivations of our paper was to recover an expression analogous to the optimal regularization parameter for the ridge regression estimator given by Hastie et al. (2022) in their Corollary 6 to Theorem 6. Their results show that using high-dimensional asymptotics where k/n → γ, the oracle choice of regularization parameter for the ridge estimator (assuming xi and ϵi are independent and that Σ = Ik ) takes the form λ∗ = nE[ϵ2i ]γ/∥β∥2 . Assuming conditional homoskedasticity, our formulae in (8) and (9) become analogous to Hastie et al. (2022)’s formula, but with the unknown ||βn ||2 replaced by either the prior mean of ∥b∥2 or the maximum value of ∥b∥2 .

4

Simulations

We now examine the extent to which the asymptotic approximations presented in Section 3 capture the finite-sample behavior of the ridge estimator. For this purpose, we present 17

Monte-Carlo simulations to calculate and compare the prediction risk of the following estimators: 1. The ridge regression with the regularization parameter chosen to minimize the worstcase approximate excess-risk criterion in (22). 2. The ridge estimator tuned according to the standard leave-one-out cross-validation; see Section 3 in Patil, Wei, Rinaldo, and Tibshirani (2021). 3. The OLS estimator. Henceforth, we refer to these estimators as minimax ridge, LOO-CV ridge, and OLS, respectively. Recall that Nk denotes a multivariate Gaussian distribution of dimension k.

4.1

DGP-1

We first consider the stylized Gaussian homoskedastic linear regression model yi = x⊤ i βn + ϵi ,

xi ∼ Nk (0, σx2 Ik ),

ϵi ∼ Nk (0, σ 2 ),

b βn = β0 + √ , n

2 with β0 = (0, . . . , 0)⊤ . We assume that ϵi is independent of xi . Let Σ = E[xi x⊤ i ] = σx Ik , 2 2 2 2 and Ω = E[xi x⊤ i ϵi ]. Under homoskedasticity and independence, Ω = σ Σ = σ σx Ik .

Throughout this subsection we set k = 10, σx2 = 1, σ 2 = 1, and b = (1, . . . , 1)⊤ . √ We set B ≡ ∥b∥ = 10. In our simulations, we consider six different sample sizes: n ∈ {500, 1000, 1500, 2000, 2500, 3000}. Each simulation uses 2,000 Monte Carlo repetitions. We use the term minimax ridge to refer to the ridge estimator that uses the regularization parameter given below in (27) and β0 as the reference vector. Since the vector of covariates has isotropic features (Σ = σx2 Ik ), the regularization parameter is chosen as in Theorem 3; that is: b bn,minimax = n trace(Ω) , λ 2 2 σ bx B

(27)

where σ bx2 is estimated from the sample covariance matrix of the covariates, and Ω is consis-

2 tently estimated by an estimator of E[xi x⊤ i ϵi ], obtained by replacing ϵi with OLS residuals

and imposing conditional homoskedasticity. For the LOO-CV (leave-one-out cross-validation) ridge, we use the leave-one-out crossvalidation procedure as in Patil et al. (2021, Section 3).7 The candidate set of regularization parameters λn = nc is based on a nonuniform grid of 500 strictly positive candidate values of c spanning [0.001, 5]. We let RidgeCV choose the minimizer of the leave-one-out criterion over that grid. 7

We set sklearn.linear model.RidgeCV with cv=None and fit intercept=False.

18

n = 500

OLS Minimax ridge OLS limit Minimax ridge limit OLS limit mean Minimax ridge limit mean

n = 1000

OLS Minimax ridge OLS limit Minimax ridge limit OLS limit mean Minimax ridge limit mean

n = 3000

OLS Minimax ridge OLS limit Minimax ridge limit OLS limit mean Minimax ridge limit mean

Density

0.8 0.6 0.4 0.2 0.0

Density

0.8 0.6 0.4 0.2 0.0

Density

0.8 0.6 0.4 0.2 0.0

2

1

0 n( 1

1

n, 1)

2

Figure 1: Distribution of the first coefficient under DGP-1 For each sample size n and m = 1, . . . , 2000, we proceed as follows: (m)

1. Draw an independent sample Dn

(m)

= {(Yi

(m)

, Xi

)}ni=1 from the Gaussian model. (m)

2. Calculate the minimax ridge, LOO-CV ridge, and OLS using Dn . (m) 3. Calculate the prediction error for each estimator βbn ,

h i E(Y,X)∼(βn ,P) (Y − X ⊤ βbn(m) )2 = σ 2 + (βbn(m) − βn )⊤ Σ(βbn(m) − βn ). We then approximate the prediction risk by averaging the prediction errors across repetitions:

2000 h i 1 X E(Y,X)∼(βn ,P) (Y − X ⊤ βbn(m) )2 . 2000 m=1

Figure 1 reports the distribution of

√

n(βb1 −βn,1 ) for the OLS and minimax ridge estima19

2

tor, with the Gaussian limits implied by Theorem 1: N1 (0, σσ2 ) for OLS and N1 (−λb1 /(σx2 + x

2 2

σx λ), (σσ2 +λ) 2 ) for minimax ridge. The approximation is accurate already at n = 500 and x

remains so as n grows. Consistent with our asymptotic approximations, the minimax ridge distribution is shifted away from the origin, but markedly more concentrated than the OLS distribution. In contrast, OLS is correctly centered but presents a large dispersion. We now compare the estimators using the excess prediction risk; that is, R(b aλn ; βn , Pn )− 2 σ . Figure 2 compares the excess prediction risk of the three estimators relative to minimax ridge estimator across different sample sizes. The figure suggests that, in this design, our recommended minimax ridge outperforms both LOO-CV ridge and OLS. Across sample sizes, the risk of LOO-CV ridge is about 20 percent higher than that of minimax ridge. OLS LOO-CV ridge Minimax ridge

Relative excess prediction risk (%)

120 100 80 60 40 20 0 500

1000

1500 2000 Sample size

2500

3000

Figure 2: Relative excess prediction risk under DGP-1 We next evaluate our approximation to the excess prediction risk of ridge estimators presented in Theorem 2. Figure 3 displays a comparison between the Monte Carlo estimate of excess prediction risk R(b aλn ; βn , Pn ) − σ 2 (solid orange line) and our feasible approximate b Ω) b (shaded region) for several values of the scaled regularization paexcess-risk Rne (λ; b, Σ, rameter λn /n ∈ [0.001, 5].8 In all the panels, we scale the vertical axis by the sample size n and report results consistent with (19), showing that the approximate excess-risk calculation tracks the excess-risk curve closely. Figure 3 also presents the approximate excess risk minimizer (vertical dashed line), which is the median (across repetitions) of λn /n based on (27). For all three sample sizes (n = 500, 1,000, and 3,000), the median selected tuning ratio is approximately 0.99, which is close to the theoretical value λn /n = 1. 8 e b Ω) b is data-dependent; therefore, for each For each value of λ, the approximate excess-risk Rn (λ; b, Σ, simulation we have a different value. To capture its range of values, we report the central 95 percent range after dropping the lowest 2.5 percent and highest 2.5 percent of values for each candidate ratio.

20

n times excess risk

11

n = 500

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 0.990

n = 1000

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 0.990

n = 3000

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 0.990

10 9 8 7 6 5

n times excess risk

11 10 9 8 7 6 5

n times excess risk

11 10 9 8 7 6 5 0

1

2

Penalty ratio n/n

3

4

5

Figure 3: Risk curves under DGP-1 bn /n) used by the LOO-CV and We also compare the scaled regularization parameter (λ minimax ridge estimators. Figure 10 in Appendix B.2 confirms that the scaled regularization parameter of the minimax ridge—based on (27)—is close to 1, which is the limit of λn /n, as expected. By contrast, the scaled regularization parameter of the LOO-CV ridge is widely dispersed across the values it can take.

4.2

DGP-2

We next consider a low-dimensional version of the weak-signal design based on Shen and Xiu (2025, Section 3.1). Our asymptotic theory in Section 3 considers models in which the number of covariates is small relative to the sample size. Thus, before considering a high-dimensional version of this design, it is useful to first understand the performance of our approximations in a non-isotropic design with a smaller number of covariates. We note that this design (and also the one presented in the next section) implies a minor departure

21

from Assumption 2 (strict stationarity) by allowing heterogeneous variances. Following Shen and Xiu (2025), we let the outcome be generated by the linear model yi = x⊤ i βn + ϵi , with n = 500 and k = 50. The ridge reference vector is β0 = 0k , so the local √ √ 1/2 1/2 parameter is b = n(βn − β0 ) = nβn . The design matrix is generated as X = Σ1 ZΣ2 , where Z ∈ Rn×k has i.i.d. standard normal entries, Σ1 ∈ Rn×n controls the dependence across observations, and Σ2 ∈ Rk×k controls the dependence across regressors. For Σ1 , we follow Shen and Xiu (2025) and set (Σ1 )ij = 2−|i−j| for 1 ≤ i, j ≤ n. For Σ2 , we also follow Shen and Xiu (2025) and construct the second moments of covariates by orthogonal diagonalization so that Σ2 = QΛQ⊤ . The orthogonal matrix Q is randomly drawn once, the eigenvalues in Λ are drawn once from U [0.5, 1.5], and Σ2 is then held fixed throughout the simulation. We note that, by construction, the matrix Σ2 is not isotropic. In Appendix B.2.2 we describe the distribution of eigenvalues of this matrix and the coefficient vector used in the simulation. Since the diagonal entries of Σ1 are equal to one, we can show that the population prediction covariance is Σ = E[xi x⊤ i ] = Σ2 . The errors are independent of the regressors and are generated with diagonal heteroskedasticity, with the diagonal entries 2 ). Thus, the of the error-variance matrix drawn once from U [0.5, 1.5] and held fixed (σϵ,i P 2 , which is asymptotic variance Ω in Assumption 1 is Ω = σ̄ϵ2 Σ2 , where σ̄ϵ2 ≡ n−1 ni=1 σϵ,i

close to 1 in this design consistent with the strong law of large numbers. Following the coefficient-generation model in Shen and Xiu (2025, Section 3.1), we first draw a preliminary coefficient vector βe = (βe1 . . . , βek ), with independent coordinates satisfying βej = 0 with probability 0.2 and βej ∼ N1 (0, 1.25) with probability 0.8. For each target value of R2 (5%, 20%, and 50%), we rescale this same draw to obtain the corresponding true coefficient vector βn . The resulting βn is then held fixed across the 2, 000 Monte Carlo repetitions. For the minimax criterion, we calibrate the radius B from the same coefficient-generation scheme: for each target R2 , we draw 2,000 coefficient vectors, rescale each draw to the tar√ get R2 , compute n∥βn ∥, and use the 90th percentile as the baseline value of B. The √ √ n normalization is the one used in the local parameter b = nβn . To accommodate the finite-sample effect of estimating a non-negligible number of coefficients, we adjust the approximate excess risk function by multiplying the variance term by n/(n − k). Since this design is not isotropic, minimax ridge uses the general non-isotropic approximate excess-risk criterion implied by Theorem 2. In each Monte Carlo repetition, we draw a new sample 2,(m) b (m) = X (m)⊤ X (m) /(n − k), σ (X (m) , Y (m) ), and we plug in the sample analogues Σ bϵ,df = df (m)2 2,(m) b (m) 1 Pn 2 b (m) = σ bi,OLS , and Ω bϵ,df Σ i=1 u df df . For each target R and each Monte Carlo repetin−k tion, we follow the same steps as in DGP-1. b(m) b∗,(m) , where λ b∗,(m) In each Monte Carlo repetition m, we choose λ n,minimax = (n − k) × λ

22

1.0

R 2 = 5%

OLS simulation Minimax ridge simulation OLS Gaussian approximation Minimax ridge Gaussian approximation Minimax ridge Gaussian approximation mean = -0.386

R 2 = 20%

OLS simulation Minimax ridge simulation OLS Gaussian approximation Minimax ridge Gaussian approximation Minimax ridge Gaussian approximation mean = -0.501

R 2 = 50%

OLS simulation Minimax ridge simulation OLS Gaussian approximation Minimax ridge Gaussian approximation Minimax ridge Gaussian approximation mean = -0.181

Density

0.8 0.6 0.4 0.2 0.0 1.0

Density

0.8 0.6 0.4 0.2 0.0 1.0

Density

0.8 0.6 0.4 0.2 0.0

3

2

1

0 n( 1

n, 1)

1

2

3

Figure 4: Distribution of the first coefficient under DGP-2 minimizes, over a grid of candidate values, the sum of (28) and (29) defined below:   b (m) + λIk )−1 Σ b (m) (Σ b (m) + λIk )−1 , λ2 B 2 µmax (Σ df df df

(28)

h i n b (m) + λIk )−1 Ω b (m) (Σ b (m) + λIk )−1 Σ b (m) . trace (Σ df df df df n−k

(29)

Relative to the criterion in (22), the only modification consists of replacing the plug-in covariance matrices by their degrees-of-freedom adjusted analogues and multiplying the variance component by n/(n − k). After the minimization, the reported ridge tuning pab(m) bn /n. rameter is λ , so all selected tuning ratios are reported on the scale λ n,minimax

For the LOO-CV ridge estimator, we use the same implementation as in DGP-1.9 The number of Monte Carlo repetitions is 2,000 for each target value of R2 . The following figures report our simulation results. Figure 4 reports the simulated dis9

sklearn.linear model.RidgeCV with cv=None and fit intercept=False.

23

tribution of

√

n(βb1 − βn,1 ) for OLS and minimax ridge, together with the Gaussian approx-

imations implied by our theory. The Gaussian curves are computed using the population matrices in the simulation, rather than fitted to the Monte Carlo histograms. Similar to the pattern observed in DGP-1, comparing with the OLS results, the distribution of the minimax ridge estimator is slightly biased but has a smaller variance. 240

OLS LOO-CV ridge mean LOO-CV ridge central 95% range Minimax ridge

Relative excess prediction risk (%)

200 160 120 80 40 0 5

20

R 2 (%)

50

Figure 5: Relative excess prediction risk under DGP-2, LOO-CV and minimax We next compare the three estimators using excess prediction risk. Figure 5 reports the excess prediction risk of OLS and LOO-CV ridge relative to minimax ridge. The shaded region reports the central 95 percent range of the repetition-specific relative excess risk of LOO-CV ridge. The figure shows that OLS is dominated by both ridge estimators, especially when R2 is small. LOO-CV ridge and minimax ridge have very similar excess prediction risk: LOO-CV is slightly above minimax ridge for the lower two values of R2 and slightly below minimax ridge when R2 = 50%. We next evaluate our approximation to the excess prediction risk of ridge estimators presented in Theorem 2 under DGP-2. Figure 6 displays a comparison between the Monte  Carlo estimate of the scaled excess prediction risk, n R(b aλn ; βn , Pn ) − σ 2 (solid orange b Ω) b (shaded region). To line), and our scaled feasible approximate excess risk nRne (λ; b, Σ, capture the range of the feasible approximation (across repetitions), we report the central 95 percent range after dropping the lowest 2.5 percent and highest 2.5 percent of values for each candidate tuning ratio. In this lower-dimensional version of the design, the feasible approximate excess-risk calculation tracks the Monte Carlo risk curve closely, including the location of the low-risk region. Appendix B.2.2 reports additional diagnostics for the eigenvalues of Σ2 , the coefficient vectors, and the selected tuning ratios.

24

n times excess risk

250

R 2 = 5%

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 1.845 Monte Carlo estimate of excess risk minimizer = 1.925

R 2 = 20%

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 0.333 Monte Carlo estimate of excess risk minimizer = 0.333

R 2 = 50%

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 0.095 Monte Carlo estimate of excess risk minimizer = 0.095

200 150 100 50 0

n times excess risk

250 200 150 100 50 0

n times excess risk

250 200 150 100 50

0 0.00

0.25

0.50

0.75

1.00 1.25 Penalty ratio n/n

1.50

1.75

2.00

Figure 6: Risk curves under DGP-2

4.3

DGP-3

We next consider the high-dimensional version of the weak-signal design based on Section 3.1 of Shen and Xiu (2025). Our asymptotic theory focused on models in which the number of covariates is small relative to the sample size. Thus, it is of interest to understand the performance of our approximations in models where the number of covariates are roughly of the same order as the sample size. The setup is exactly the same as the setup in DGP-2 in the previous subsection, except here we set k = 300. The following figures report our simulation results. Figure 7 reports the simulated dis√ tribution of n(βb1 − βn,1 ) for OLS and minimax ridge, together with the Gaussian approximations implied by our theory. The Gaussian curves are computed using the population matrices in the simulation, rather than fitted to the Monte Carlo histograms. Similar to the pattern observed in DGP-1 and DGP-2, comparing with the OLS results, the distribution of the minimax ridge estimator is biased by shrinkage but has a smaller variance.

25

2.0

R 2 = 5%

OLS simulation Minimax ridge simulation OLS Gaussian approximation Minimax ridge Gaussian approximation Minimax ridge Gaussian approximation mean = -0.139

R 2 = 20%

OLS simulation Minimax ridge simulation OLS Gaussian approximation Minimax ridge Gaussian approximation Minimax ridge Gaussian approximation mean = -0.223

R 2 = 50%

OLS simulation Minimax ridge simulation OLS Gaussian approximation Minimax ridge Gaussian approximation Minimax ridge Gaussian approximation mean = -0.241

Density

1.5 1.0 0.5 0.0 2.0

Density

1.5 1.0 0.5 0.0 2.0

Density

1.5 1.0 0.5 0.0

4

2

0 n( 1

n, 1)

2

4

Figure 7: Distribution of the first coefficient under DGP-3 We next compare the three estimators using excess prediction risk. Figure 8 reports the excess prediction risk of OLS and LOO-CV ridge relative to minimax ridge. The shaded region reports the central 95 percent range of the repetition-specific relative excess risk of LOO-CV ridge. The figure shows that OLS is dominated by both ridge estimators, especially when R2 is small. LOO-CV ridge and minimax ridge have comparable excess prediction risk, although LOO-CV ridge has slightly lower mean excess risk than minimax ridge in this high-dimensional design. Thus, in the regime where k/n is large, the approximationbased rule should not be interpreted as uniformly improving on LOO-CV, but it remains competitive with the standard cross-validation benchmark. We next evaluate our approximation to the excess prediction risk of ridge estimators presented in Theorem 2 under DGP-3. Figure 9 displays a comparison between the Monte  Carlo estimate of the scaled excess prediction risk, n R(b aλn ; βn , Pn ) − σ 2 (solid orange b Ω) b (shaded region). To line), and our scaled feasible approximate excess risk nRne (λ; b, Σ, capture the range of the feasible approximation (across repetitions), we report the central 26

OLS LOO-CV ridge mean LOO-CV ridge central 95% range Minimax ridge

3000 2000 Relative excess prediction risk (%)

1000 10 5 0 5 10 15 20 25 30

5

20

R 2 (%)

50

Figure 8: Relative excess prediction risk under DGP-3, LOO-CV and minimax 95 percent range after dropping the lowest 2.5 percent and highest 2.5 percent of values for each candidate tuning ratio. In this high-dimensional version of the design, the feasible approximation is less accurate than in DGP-2, especially away from the low-risk region. This pattern is expected because k/n = 0.6 is far from the low-dimensional asymptotic framework used to derive our approximation. Appendix B.2.3 reports additional diagnostics for the eigenvalues of Σ2 , the coefficient vectors, and the selected tuning ratios.

5

Conclusion

We presented a simple Gaussian approximation to the finite-sample distribution of the ridge regression estimator. Our approximation is based on nonstandard asymptotics where i) we let the estimator’s regularization parameter grow proportionally to the sample size; and ii) we treat the population regression coefficients as local to the reference vector that defines the estimator’s direction of shrinkage. In contrast to other asymptotic approximations in the literature, we allow for general forms of heteroskedasticity and autocorrelation in the data generating process (at the cost of considering a low-dimensional model where the number of covariates is not allowed to grow with the sample size). We used our simple Gaussian approximation to propose two new strategies to select the regularization parameter for the ridge regression estimator. The suggested strategies select the regularization parameter to minimize either average or worst-case excess prediction risk, where risk is computed using our suggested Gaussian approximation.

27

1000

R 2 = 5%

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 11.949 Monte Carlo estimate of excess risk minimizer = 11.710

R 2 = 20%

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 2.481 Monte Carlo estimate of excess risk minimizer = 2.481

R 2 = 50%

Central 95% range of approximate excess risk Monte Carlo estimate of excess risk Approximate excess risk minimizer = 0.652 Monte Carlo estimate of excess risk minimizer = 0.652

n times excess risk

800 600 400 200 0 1000

n times excess risk

800 600 400 200 0

n times excess risk

1000 800 600 400 200 0 0

2

4

6 8 Penalty ratio n/n

10

12

Figure 9: Risk curves under DGP-3

References Abadie, A. and M. Kasy (2019): “Choosing among Regularized Estimators in Empirical Economics: The Risk of Machine Learning,” The Review of Economics and Statistics, 101, 743–762. Andrews, I. and A. Mikusheva (2022): “Optimal decision rules for weak GMM,” Econometrica, 90, 715–748. Anup, A. K. and G. Maddala (1984): “Ridge estimators for distributed lag models,” Communications in Statistics-Theory and Methods, 13, 217–225. Atanasov, A., J. A. Zavatone-Veth, and C. Pehlevan (2024): “Risk and cross validation in ridge regression with correlated samples,” arXiv preprint arXiv:2408.04607.

28

Brockwell, P. J. and R. A. Davis (2013): Time series: theory and methods, Springer Science & Business Media. Cambanis, S., S. Huang, and G. Simons (1981): “On the theory of elliptically contoured distributions,” Journal of Multivariate Analysis, 11, 368–385. Cattaneo, M. D., M. Jansson, and W. K. Newey (2018): “Inference in linear regression models with many covariates and heteroscedasticity,” Journal of the American Statistical Association, 113, 1350–1361. DasGupta, A. (2008): Asymptotic Theory of Statistics and Probability, Springer Verlag. Dobriban, E. and S. Wager (2018): “High-dimensional asymptotics of prediction: Ridge regression and classification,” The Annals of Statistics, 46, 247–279. Dou, L. and U. K. Müller (2021): “Generalized Local-to-Unity Models,” Econometrica, 89, 1825–1854. Ferguson, T. S. (1967): Mathematical Statistics: A Decision Theoretic Approach, New York: Academic Press. Friedman, J., T. Hastie, and R. Tibshirani (2017): The elements of statistical learning: data mining, inference and prediction, vol. 1 of Series in Statistics, New York: Springer, second edition ed. Gibbons, D. G. (1981): “A Simulation Study of Some Ridge Estimators,” Journal of the American Statistical Association, 76, 131–139. Hansen, B. (2022): Econometrics, Princeton University Press. Hansen, B. E. (2016): “Efficient shrinkage in parametric models,” Journal of Econometrics, 190, 115–132. Hastie, T., A. Montanari, S. Rosset, and R. J. Tibshirani (2022): “Surprises in high-dimensional ridgeless least squares interpolation,” The Annals of Statistics, 50, 949– 986. Hirano, K. and J. H. Wright (2017): “Forecasting with model uncertainty: Representations and risk reduction,” Econometrica, 85, 617–643. Hoerl, A. E. (1962): “Application of ridge analysis to regression problems,” Chemical Engineering Progress, 58, 54–59.

29

Hoerl, A. E. and R. W. Kennard (1970): “Ridge regression: Biased estimation for nonorthogonal problems,” Technometrics, 12, 55–67. Hoerl, R. W. (2020): “Ridge regression: a historical context,” Technometrics, 62, 420– 425. Hsu, D., S. M. Kakade, and T. Zhang (2012): “Random design analysis of ridge regression,” in Conference on learning theory, JMLR Workshop and Conference Proceedings, 9–1. Knight, K. and W. Fu (2000): “Asymptotics for lasso-type estimators,” Annals of statistics, 1356–1378. Kock, A. B., R. S. Pedersen, and J. R.-V. Sørensen (2026): “Data-driven tuning parameter selection for high-dimensional vector autoregressions,” Journal of the American Statistical Association, 121, 289–299. Liu, X., S. Zheng, and X. Feng (2020): “Estimation of error variance via ridge regression,” Biometrika, 107, 481–488. Montiel Olea, J. L., C. Rush, A. Velez, and J. Wiesel (2026): “The distributionally robust prediction error of the LASSO and related estimators,” The Annals of Statistics, 54, 1006–1027. Mourtada, J. (2022): “Exact minimax risk for linear least squares, and the lower tail of sample covariance matrices,” The Annals of Statistics, 50, 2157–2178. Mourtada, J. and L. Rosasco (2022): “An elementary analysis of ridge regression with random design,” Comptes Rendus. Mathématique, 360, 1055–1063. Patil, P., J.-H. Du, and R. J. Tibshirani (2024): “Optimal ridge regularization for out-of-distribution prediction,” arXiv preprint arXiv:2404.01233. Patil, P., Y. Wei, A. Rinaldo, and R. Tibshirani (2021): “Uniform consistency of cross-validation estimators for high-dimensional ridge regression,” in International conference on artificial intelligence and statistics, PMLR, 3178–3186. Phillips, P. C. (1988): “Regression theory for near-integrated time series,” Econometrica: Journal of the Econometric Society, 1021–1043. Powell, J. L. (2017): “Identification and Asymptotic Approximations: Three Examples of Progress in Econometric Theory,” Journal of Economic Perspectives, 31, 107–124.

30

Shen, Z. and D. Xiu (2025): “Can Machines Learn Weak Signals?” Working Paper 33421, National Bureau of Economic Research, Cambridge, MA. Staiger, D. and J. H. Stock (1997): “Instrumental Variables Regression with Weak Instruments,” Econometrica, 65, 557–586. Swindel, B. F. (1976): “Good ridge estimators based on prior information,” Communications in Statistics-Theory and Methods, 5, 1065–1075. Velez, A. (2024): “On the asymptotic properties of debiased machine learning estimators,” arXiv preprint arXiv:2411.01864.

31

A

Proofs of Main Results

A.1

Proof of Proposition 1

Proof. Let X be the n × k matrix that contains x⊤ i in its i-th row. Let Y and ϵ be the n × 1 vectors that contain yi and ϵi (respectively) in their i-th row. Define the matrices
b ≡ X ⊤ X/n and A b≡Σ b + (λn /n)Ik . Using this notation, βbλ = A b−1 X ⊤ Y /n + (λn /n)β0 Σ n b −1 (X ⊤ Y )/n. Algebra shows and βbOLS = Σ    √  √ √ b−1 X ⊤ ϵ/ n − (λn /n) n (βn − β0 ) n βbλn − βn = A  √ √  b −1 X ⊤ ϵ/ n . n βbOLS − βn = Σ

(30)

Therefore,  √     √  √ √ b−1 (λn /n) n (βn − β0 ) , b−1 − Σ b −1 (X ⊤ ϵ)/ n − A n βbλn − βn − n βbOLS − βn = A p

b−1 − Σ b −1 is op (1), (ii) Assumption 1 implies which is op (1) because (i) λn /n → 0 implies A √ √ p (X ⊤ ϵ)/ n is Op (1), and (iii) because we are assuming (λn /n) n (βn − β0 ) → 0 in this proposition.

A.2

Proof of Theorem 1 p

p

b ≡ X ⊤ X/n and A b ≡ Σ+(λ b b b Proof. Recall that Σ n /n)Ik . Note that A → Σ+λIk since Σ → Σ √ p p −1 b (λn /n) n (βn − β0 ) → (Σ + λIk )−1 λb. (Assumption 1) and λn /n → λ. Note also that A We conclude by using Equation (30), Assumption 1, and Slutsky’s theorem.

A.3

Proof of Theorem 2

Proof. We establish the result in two steps. We first show that the excess prediction risk Rn (âλn ; βn , Pn ) is approximated by Rne (λ; b, Σ, Ω) up to a remainder error of size o(n−1 ). We then show that Rne (λ; b, Σ, Ω), which depends on population parameters (Σ, Ω), can b Ω), b which uses consistent estimators (Σ, b Ω), b up to error of be approximated by Rne (λ; b, Σ, size o(βn ,Pn ) (1/n). The conclusion follows by combining these two steps and the triangle inequality. Step 1: For a new draw (x⊤ , ϵ) from the stationary distribution P associated to Pn , we use βn to construct a new outcome-covariate pair (y, x) = (x⊤ βn + ϵ, x). The exact finite-sample prediction risk of âλn is:     Rn (âλn ; βn , Pn ) = E(βn ,Pn ) (y − x⊤ βbλn )2 = σ 2 + E(βn ,Pn ) (βbλn − βn )⊤ Σ(βbλn − βn ) .

32

Defining Zn ≡

√ b n(βλn − βn ), this becomes Rn (âλn ; βn , Pn ) = σ 2 +

  1 E(βn ,Pn ) Zn⊤ ΣZn . n

Note that (18) can be rewritten as Rne (λ; b, Σ, Ω) =

1 EZ [Z ⊤ ΣZ], n

where Z ∼ Nk (−(Σ + λIk )−1 λb, (Σ + λIk )−1 Ω(Σ + λIk )−1 ). Therefore, Rn (âλn ; βn , Pn ) − σ 2 − Rne (λ; b, Σ, Ω) =

   1 E(βn ,Pn ) Zn⊤ ΣZn − EZ [Z ⊤ ΣZ] , n

so it suffices to show   E(βn ,Pn ) Zn⊤ ΣZn → EZ [Z ⊤ ΣZ]. d

Under the assumptions of Theorem 2, Slutsky’s theorem and Theorem 1 imply that Zn → Z. Theorem 6.2 in DasGupta (2008) implies that under assumption iii) of Theorem 2 (which   is a sufficient condition for uniform integrability) we have E(βn ,Pn ) Zn⊤ ΣZn → EZ [Z ⊤ ΣZ]. This establishes the oracle approximation Rn (âλn ; βn , Pn ) − σ 2 = Rne (λ; b, Σ, Ω) + o

1 n

.

  b Ω) b − Rne (λ; b, Σ, Ω) = o(β ,P ) (1). We Step 2: It is sufficient to show that n Rne (λ; b, Σ, n n b b claim that this follows from the consistency of the estimators (Σ, Ω) of (Σ, Ω). To see this, note that
nRne (λ; b, Σ, Ω) ≡ λ2 b⊤ (Σ + λIk )−1 Σ(Σ + λIk )−1 b + trace (Σ + λIk )−1 Ω(Σ + λIk )−1 Σ . Then, we conclude by using Slutsky’s theorem and continuous mapping theorem.

A.4

Proof of Theorem 3

Proof. Under isotropic features, Σ = σx2 Ik , the approximation of the excess risk of ridge regression in (18) simplifies to Rne (λ; b, σx2 , Ω) =

1 λ2 σx2 1 σx2 ⊤ b b + trace(Ω). n (σx2 + λ)2 n (σx2 + λ)2

33

(31)

We use the previous closed-form expression to obtain ∂λ Rne (λ; b, σx2 , Ω) =

2σx2 1 (σ 2 (b⊤ b)λ − trace(Ω)) . n (σx2 + λ)3 x

(32)

Note that the previous expression is negative if and only if λ < λ∗ ≡ trace(Ω) . Therefore, σ 2 (b⊤ b) x

Rne (λ; b, σx2 , Ω) as a function of λ is decreasing on [0, λ∗ ) and increasing on [λ∗ , ∞). As a result, we have that λ∗ ∈ arg minλ≥0 Rne (λ; b, σx2 , Ω) for any given (b, σx2 , Ω). We first establish part (1) of Theorem 3. Note that when ∥ · ∥ is the Euclidean norm, the set B ≡ {b ∈ Rk : ∥b∥ ≤ B} = {b ∈ Rk : b⊤ b ≤ B 2 }. Since the coefficient on b⊤ b in the excess risk approximation formula defined by (31) is nonnegative for every λ ≥ 0, the worst-case value of the approximate excess risk is attained when b⊤ b = B 2 . Hence the minimax problem reduces to minimizing the function Qminimax (λ) ≡

1 λ2 σx2 1 σx2 2 B + trace(Ω), n (σx2 + λ)2 n (σx2 + λ)2

λ ≥ 0,

which as a function of λ achieves its minimum at λ∗minimax = trace(Ω) . σ2 B 2 x

We then establish part (2) of Theorem 3. Suppose Eπ [b⊤ b] < ∞. Taking expectation with respect to π gives   1 λ2 σx2 1 σx2 ⊤ Qaverage (λ) := Eπ Rne (λ; b, σx2 , Ω) = E [b b] + trace(Ω), π n (σx2 + λ)2 n (σx2 + λ)2

λ ≥ 0.

If Eπ [b⊤ b] > 0, algebra shows that Qaverage (·) achieves its minimum at λ∗average = σtrace(Ω) 2 E [b⊤ b] . π x

B

Additional Results

B.1

Optimal Choice of λ for Elliptically Contoured Distributions

In this appendix, we extend the result on choosing λ to minimize average approximate excess risk to the more general case in which b follows an elliptically contoured distribution. Following Cambanis, Huang, and Simons (1981), a k-dimensional random vector X is said to follow an elliptical distribution with location parameter µX ∈ Rk , nonnegative definite matrix V ∈ Rk×k , and characteristic generator ψ if the characteristic function of X − µX satisfies ϕX−µX (t) = ψ(t⊤ V t), We denote this by X ∼ ECk (µX , V, ψ).

34

t ∈ Rk .

Now suppose b ∼ ECk (0, Ik , ψ) has finite second moments. Algebra shows that E[bb⊤ ] =

E∥b∥2 Ik . k

Therefore, i h b Σ b + λIk )−1 b b + λIk )−1 Σ( Eb∼π b⊤ (Σ h  i b + λIk )−1 Σ( b Σ b + λIk )−1 bb⊤ = Eb∼π tr (Σ   b Σ b + λIk )−1 Eb∼π [bb⊤ ] b + λIk )−1 Σ( = tr (Σ   2 −1 b b −1 E∥b∥ b = tr (Σ + λIk ) Σ(Σ + λIk ) Ik k  E∥b∥2  b b Σ b + λIk )−1 . tr (Σ + λIk )−1 Σ( = k Hence the average approximate excess risk takes the same form as in Section 3.3.2, with C 2 replaced by E∥b∥2 /k. It follows that the corresponding characterization of the optimal choice of λ extends immediately to the class of elliptical contoured distributions with finite second moment.

B.2

Additional Simulation Results

B.2.1

DGP-1

DGP-1 is the isotropic example. The main additional object of interest is the selected bn /n. Under the population version of the minimax rule in this design, the tuning ratio λ optimal scaled regularization parameter equals c∗ =

trace(Ω) = 1. σx2 B 2

bn /n for LOO-CV ridge Figure 10 reports the distribution of the selected tuning ratio λ and minimax ridge. The feasible minimax rule is tightly centered near the population value c∗ = 1, with dispersion decreasing as the sample size increases. In contrast, the LOO-CV tuning ratio is much more dispersed across Monte Carlo repetitions. This helps explain the risk comparison: in DGP-1, minimax ridge uses a tuning parameter close to the risk-minimizing value, whereas LOO-CV often selects substantially different amounts of shrinkage.

35

17.5

n = 500

LOO-CV ridge Minimax ridge

n = 1000

LOO-CV ridge Minimax ridge

n = 3000

LOO-CV ridge Minimax ridge

15.0 Density

12.5 10.0 7.5 5.0 2.5 0.0 17.5 15.0 Density

12.5 10.0 7.5 5.0 2.5 0.0 17.5 15.0 Density

12.5 10.0 7.5 5.0 2.5 0.0 0.00

0.25

0.50

0.75 1.00 1.25 Selected tuning ratio n/n

1.50

1.75

2.00

bn /n) used by the LOO-CV and minimax ridge Figure 10: Scaled regularization parameter (λ estimators under DGP-1 B.2.2

DGP-2

This appendix records diagnostic information for the covariance matrix and coefficient vectors used in DGP-2. These diagnostics are useful because the design departs from the isotropic benchmark: the prediction covariance is Σ2 , not a scalar multiple of the identity matrix, and the coefficient vector is a fixed draw from the spike-and-slab model described in the main text. First, we verify the equality Σ = E[xi x⊤ i ] = Σ2 used in the simulation section. Let ei ⊤ denote the ith canonical basis vector in Rn . Since the ith row of X is x⊤ i = ei X, we have ⊤ ⊤ E[xi x⊤ i ] = E[X ei ei X] h i 1/2 1/2 1/2 1/2 = Σ2 E Z ⊤ Σ1 ei e⊤ i Σ1 Z Σ2 .

36

For any deterministic n × n matrix A and an n × k matrix Z with i.i.d. standard normal 1/2

1/2

entries, E[Z ⊤ AZ] = trace(A)Ik . Applying this identity with A = Σ1 ei e⊤ i Σ1

gives

  1/2 1/2 1/2 ⊤ 1/2 E[xi x⊤ ] = Σ trace Σ e e Σ Ik Σ2 i i i 2 1 1 1/2

1/2

= Σ2 (e⊤ i Σ1 ei )Ik Σ2 1/2

1/2

= Σ2 (Σ1 )ii Ik Σ2 = Σ2 ,

where the last equality uses (Σ1 )ii = 1. Similarly, because the errors are independent of the 2 , the corresponding variance matrix for regressors, mean zero, and have fixed variances σϵ,i

the score is proportional to Σ2 : n

1X 2 E[xi x⊤ i ϵi ] = n i=1

n

1X 2 σϵ,i n

! Σ2 .

i=1

k = 50 Isotropic benchmark = 1

2.0

Density

1.5 1.0 0.5 0.0

0.6

0.8

1.0 Eigenvalue of 2

1.2

1.4

Figure 11: Eigenvalues of Σ2 under DGP-2 Figure 11 plots the eigenvalues of the fixed Σ2 matrix used in DGP-2. The eigenvalues range from 0.502 to 1.479, with mean 0.981 and max-min spread 0.977. The vertical dashed line marks the isotropic benchmark value one. The figure shows that the design is centered near the isotropic benchmark on average but has meaningful dispersion in eigenvalues. Figure 12 displays the fixed coefficient vector used in DGP-2 after rescaling the same preliminary draw to the three target values of R2 . The labels report ∥βn ∥ for each target value. In this draw, 41 out of the 50 coordinates are nonzero. The figure makes explicit how

37

0.4

R 2 = 5%;

n = 0.2437

R 2 = 20%;

n = 0.5312

R 2 = 50%;

n = 1.0623

n, j

0.2 0.0 0.2 0.4 0.4

n, j

0.2 0.0 0.2 0.4 0.4

n, j

0.2 0.0 0.2 0.4 0

10

20

30 Coefficient index

40

50

Figure 12: Coefficient vectors under DGP-2 increasing the target R2 changes the magnitude of the coefficient vector while preserving the same coefficient pattern. bn /n for LOO-CV ridge Figure 13 reports the distribution of the selected tuning ratio λ and minimax ridge under DGP-2. The minimax rule is more concentrated because it is driven by the approximate-risk calculation and the calibrated radius B, whereas LOO-CV is more dispersed across Monte Carlo repetitions. B.2.3

DGP-3

This appendix records diagnostic information for the covariance matrix and coefficient vectors used in DGP-3. DGP-3 uses the same covariance-generation and coefficient-generation schemes as DGP-2, but increases the number of covariates from k = 50 to k = 300 while keeping n = 500. Thus, the diagnostics below should be read as the high-dimensional counterpart to Appendix B.2.2. The same algebra as in Appendix B.2.2 implies that the population prediction covariance is Σ = E[xi x⊤ i ] = Σ2 , 38

35

R 2 = 5%

LOO-CV ridge Minimax ridge

R 2 = 20%

LOO-CV ridge Minimax ridge

R 2 = 50%

LOO-CV ridge Minimax ridge

30 Density

25 20 15 10 5 0 35 30 Density

25 20 15 10 5 0 35 30 Density

25 20 15 10 5 0 0.0

0.5

1.0

1.5 Selected tuning ratio n/n

2.0

2.5

Figure 13: Selected tuning ratios under DGP-2 because the diagonal entries of Σ1 are equal to one. Since the errors are independent of 2 , the corresponding score variance the regressors, mean zero, and have fixed variances σϵ,i

matrix is

n

1X 2 E[xi x⊤ i ϵi ] = n i=1

n

1X 2 σϵ,i n

! Σ2 .

i=1

Figure 14 plots the eigenvalues of the fixed Σ2 matrix used in DGP-3. The eigenvalues range from 0.502 to 1.499, with mean 1.027 and max-min spread 0.997. As in DGP-2, the design is centered near the isotropic benchmark on average but is clearly non-isotropic. Figure 15 displays the fixed coefficient vector used in DGP-3 after rescaling the same preliminary draw to the three target values of R2 . The labels report ∥βn ∥ for each target value. In this draw, 244 out of the 300 coordinates are nonzero. The figure shows that the signal is dense but individually weak, especially for the lower target value of R2 . bn /n for LOO-CV ridge Figure 16 reports the distribution of the selected tuning ratio λ

39

k = 300 Isotropic benchmark = 1

1.6 1.4

Density

1.2 1.0 0.8 0.6 0.4 0.2 0.0

0.6

0.8

1.0 Eigenvalue of 2

1.2

1.4

R 2 = 5%;

n = 0.2252

n, j

0.15 0.10 0.05 0.00 0.05 0.10 0.15

R 2 = 20%;

n = 0.4908

n, j

0.15 0.10 0.05 0.00 0.05 0.10 0.15 0.15 0.10 0.05 0.00 0.05 0.10 0.15

R 2 = 50%;

n = 0.9815

n, j

Figure 14: Eigenvalues of Σ2 under DGP-3

0

50

100

150 Coefficient index

200

250

300

Figure 15: Coefficient vectors under DGP-3 and minimax ridge under DGP-3. Relative to DGP-2, the LOO-CV choices are substantially more dispersed, especially when R2 is small. This dispersion is consistent with the main-text 40

10

R 2 = 5%

LOO-CV ridge Minimax ridge

R 2 = 20%

LOO-CV ridge Minimax ridge

R 2 = 50%

LOO-CV ridge Minimax ridge

Density

8 6 4 2 0 10

Density

8 6 4 2 0 10

Density

8 6 4 2 0

0

1

2

3 4 5 Selected tuning ratio n/n

6

7

8

Figure 16: Selected tuning ratios under DGP-3 finding that, in the high-dimensional design, LOO-CV and minimax ridge have comparable average risk but can choose very different penalty values in individual samples.

41

Related documents

• Theory of Statistical Estimation : Contributions to Higher Order Asymptotics and Non-Regular Cases / by Masafumi Akahira#1165539
• Estimation Associated with Linear Models: A Revisitation#1149659
• The Distribution of Model Averaging Estimators and an Impossibility Result Regarding Its Estimation#499517
• Elementary Estimates: An Introduction to Nonparametrics.#198473
• Finite-Sample Risk Approximation and Risk-Consistent Tuning for Generalized Ridge Estimation in Nonlinear Models: Controlling Extreme Realizations#141729