Kernel Ridge-Type Shrinkage Estimators in Partially Linear Regression Models with Correlated Errors

1. Introduction

In many time-series applications, one must balance interpretability and flexibility in modeling. Partially linear regression models achieve this by expressing the response as a parametric linear function of some covariates alongside a nonparametric smooth function of an additional variable. Specifically, we consider the following model setup:

$y_t={\mathbf{x}}_t^T\mathsf{\beta }+f\left(z_t\right)+{\epsilon }_t, t=1,\dots ,n$ (1.1)

where $y_t$ denotes stationary response at time $t$, and its prediction depends on the p-dimensional vector of explanatory variables ${\mathbf{x}}_t=(x_{t1},\dots ,x_{tp})$ with $x_t={(x_1,\dots ,x_n)}^{\text{'}}$ and an extra univariate variable $z_t$. Also, $\mathsf{\beta }\in {\mathbb{R}}^p$ is an unknown p-dimensional parameter vector to be estimated, $f\left(.\right)$ is an unknown smooth function, and ${\epsilon }_i\text{'}s$ are the error terms.

Unlike classical linear models, (1.1) offers the interpretability of $\mathsf{\beta }$ while retaining the capacity to model complex, nonlinear effects through $f(.)$ ([14]; Robinson, 1998). However, two mojor issues often arise:

Multicollinearity among ${\mathbf{x}}_t.$ When the columns of ${\mathbf{x}}_t$ are nearly linearly dependent, ordinary least squares (OLS) or unpenalized nonparametric estimators can have extremely high variances [10].

Autocorrelation or dependence among ${\epsilon }_t$, as is common in time-series or longitudinal settings. We capture this by a first-order autoregressive process:

${\epsilon }_t=\rho {\epsilon }_{t-1}+u_t, \left|\rho \right|<1, u_t~i.i.d. N\left(0,{\sigma }^2\right), t=\mathrm{1,2},\dots ,n$ (1.2)

ensuring that ${\epsilon }_t$ is stationarily dependent over time ([7,15]). Also, to simplify the notational illustration, a matrix–vector form of the model (1.1) can be stated as

$\mathbf{y}=\mathbf{X}\mathsf{\beta }+\mathbf{f}+\mathsf{\epsilon }$ (1.3)

where$\mathbf{y}=\left(y_1,\dots ,y_n\right)$, $\mathbf{X}=({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n)\mathbf{\text{'}}$, $\mathbf{f}=\left(f\left(z_1\right),\dots ,f\left(z_n\right)\right)$ and$\mathsf{\epsilon }=\left({\epsilon }_1,\dots ,{\epsilon }_n\right)$. The main goal is to estimate the unknown parameter vector $\mathsf{\beta }$, the smooth function $f\left(t\right)$ and the mean vector $\mathit{\mu }=\mathbf{X}\mathsf{\beta }+\mathbf{f}$ based on the observations from the data set {${ y}_t, {\mathbf{x}}_t{,z}_t$}. Partially linear models enable easier interpretation of the effect of each variable and owing to the “curse of dimensionality”, these models can be preferred to a completely nonparametric regression model. Specifically, partially linear models are more practical than the classical linear model because they combine both parametric part presented numerically and nonparametric part displayed graphically.

Eventhough partially linear models have specific advantages, these models typically assume uncorrelated errors ([6,8,14]). In real-world time-series or panel data, ignoring (1.2) can degrade efficiency and produce biased inferences ([11,19]). Moreover, multicollinearity among ${\mathbf{x}}_t$​ amplifies these problems, inflating variances to the point that $\mathsf{\beta }$ estimates become unreliable ([13,20]).

To mitigate these issues, we adopt ridge-type kernel smoothing for the parametric component. In essence, a shrinkage penalty is added so that the linear portion $\mathsf{\beta }$ is “shrunk” toward zero to stabilize estimation. Mathematically, for instance, a ridge-type kernel (RTK) estimator of $\mathsf{\beta }$ can be expressed (after suitable data transformation) as in (2.3) which has been extensively studied by [17]. This ridge-inspired framework stems from Stein-type shrinkage theory, which shows that biased estimators can dominate unbiased ones, particularly in moderate-to-high dimensions or under strong correlation ([1,2]). Indeed, Stein’s paradox reveals that for $p\ge 3$ shrinking an estimator can lower the overall mean squared error despite introducing bias, a phenomenon extensively explored by [1,3]).

When errors follow (1.2), we further incorporate transformations to handle the autocovariance structure, yielding generalized ridge-type kernel (GRTK) estimators that optimally exploit the correlation $\mathit{\rho }$. These GRTK methods align with earlier biased estimation approaches for linear or partially linear models with correlated errors ([13,20,21]). By merging the Stein-type shrinkage logic with kernel smoothing, one obtains robust estimators of both $\mathit{\beta }$ and $f\left(.\right)$ that effectively address multicollinearity and autocorrelation in a unified fashion.

In this paper, we formalize and extend such modified ridge-type kernel smoothing estimators for partially linear time-series models. We focus on:

The precise data transformation to accommodate (1.2).

The interplay between shrinkage $k$ and the bandwidth $\lambda$ for kernel smoothing.

Selecting these tuning parameters ($k, \lambda$) using multiple criteria, including GCV, AICc, BIC, and RECP.

This framework is motivated by [3] broader research on shrinkage-based variable selection and penalty methods, which shows that penalized procedures (including ridge, lasso, and Stein-like estimators) can vastly reduce variance in semiparametric regressions—particularly with correlated data or limited samples. Consequently, the main goal is to fill a gap in semiparametric regression by uniting autoregressive error modeling (1.2), handling with multicollinearity and Stein-type shrinkage in kernel-based estimation.

The rest of the paper is organized as follows. Section 2 explains how we fit the autoregressive structure (1.2) and handle collinearity. Section 3 provides the mathematical formulations of our GRTK shrinkage estimators for both $\mathsf{\beta }$ and the smooth function $f\left(.\right)$, including asymptotic properties. Section 4 details the procedures for choosing $k$ and bandwidth $\lambda$. Section 5 illustrates, through extensive simulation and an airline delay time-series dataset, how the proposed methodology outperforms traditional methods under multicollinearity and autocorrelation. Finally, Section 6 concludes with a discussion of key findings and future research directions.

2. Fitting the Model Error Structure

Assume that $y_t=\left(y_1,\dots ,y_n\right)$′ is a realization from a stationary time series described by model (1.1) with autoregressive error terms satisfying the following assumption.

Assumption 2.1: The error ${\mathsf{\epsilon }}_t$ is a stationary dependent sequence with $E\left({\mathsf{\epsilon }}_t\right)=0 , Var\left({\mathsf{\epsilon }}_t\right)={\sigma }_{\epsilon }^2,$and $E{\left|{\mathsf{\epsilon }}_t\right|}^{2+\delta }<\infty$ for some constant$\delta .$ Letting ${\mathsf{\epsilon }}_{\mathit{t}}=\left({\epsilon }_1,\dots ,{\epsilon }_n\right)$ and $Var\left({\mathsf{\epsilon }}_{\mathit{t}}\right)=E\left({\mathit{\epsilon }}_t{\mathit{\epsilon }}_t^{\text{'}}\right)={\sigma }^2\mathit{R}$, where $\mathit{R}$ is a $\left(n\times n\right)$ positive definite symmetric matrix.

Here, Assumption 2.1 is standard in time-series analysis because it allows for modeling of stationary errors with finite variance. Such an assumption is common in many time-series models where the error term’s statistical properties remain constant over time. Note that the autoregressive errors follow an $n-$dimensional multivariate normal distribution with a mean zero and stationary $\left(n\times n\right)$ covariance matrix $\mathsf{\Sigma }$. We may also write this expression in the equivalent form ${\mathsf{\epsilon }}_{\mathit{t}}=\left({\epsilon }_1,\dots ,{\epsilon }_n\right)~N(0,\mathsf{\Sigma })$, where $\mathsf{\Sigma }={\sigma }^2\mathit{R}$ is a $n\times n$ covariance matrix, given by

$\mathit{R}={\displaystyle \frac{1}{1-{\rho }^2}}\left[\begin{array}{ccc}1& \rho & \begin{array}{ccc}{\rho }^2& \dots & {\rho }^{n-1}\end{array}\\ \rho & 1& \begin{array}{ccc}\rho & \dots & {\rho }^{n-2}\end{array}\\ \begin{array}{c}{\rho }^2\\ \begin{array}{c}⋮\\ {\rho }^{n-1}\end{array}\end{array}& \begin{array}{c}\rho \\ \begin{array}{c}⋮\\ {\rho }^{n-2}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}1\\ \begin{array}{c}⋮\\ {\rho }^{n-3}\end{array}\end{array}& \begin{array}{c}\dots \\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}& \begin{array}{c}{\rho }^{n-3}\\ \begin{array}{c}⋮\\ 1\end{array}\end{array}\end{array}\end{array}\right]$ (2.1a)

with the inverse

${\mathit{R}}^{-1}={\displaystyle \frac{1}{1-{\rho }^2}}\left[\begin{array}{ccc}1& -\rho & \begin{array}{ccc} 0& \dots & \begin{array}{cc}0& 0\end{array}\end{array}\\ -\rho & 1+{\rho }^2& \begin{array}{ccc}-\rho & 0& \begin{array}{cc}\dots & 0\end{array}\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ \dots \end{array}& \begin{array}{ccc}\begin{array}{c} ⋮\\ 0\end{array}& \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}⋮\\ -\rho \end{array}& \begin{array}{c}⋮\\ 1\end{array}\end{array}\end{array}\end{array}\right]$ (2.1b)

where $\rho$ denotes the correlation between ${\epsilon }_t$ and ${\epsilon }_{t-1}$, as defined before. It should be emphasized that once we have an estimate $\widehat{\mathsf{\beta }}$ of the parameter vector $\mathsf{\beta }$ and an estimator $\widehat{f}$ of the unknown smooth function $f$, the parameter $\rho$ can be estimated by the residuals from the semiparametric regression $e_t=y_t-\left({\mathbf{x}}_t\widehat{\mathsf{\beta }}+\widehat{f}\left(z_t\right)\right)$, defined as

$\widehat{\rho }={\displaystyle \frac{Cov (e_t,e_{t-1})}{\sqrt{Var\left(e_t\right)}\sqrt{Var\left(e_{t-1}\right)}}}={\displaystyle \frac{\sum _{t=2}^n{e}_t{e}_{t-1}}{\sum _{t=1}^n{e}_t^2}} \left(2.2\right)$There are several methods to estimate $\mathsf{\beta }$ and $f.$ In this work, we adopt the ridge type kernel (RTK) method discussed by [17], generalizes partial kernel method proposed by [14]. Specifically, the RTK estimator $\widehat{\mathsf{\beta }}$ of $\mathsf{\beta }$ is the form

${\widehat{\mathsf{\beta }}}_{RTK}={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}\tilde{\mathbf{y}}$ (2.3)

where$k>0$ is a shrinkage parameter, $\mathbf{I}$ is an $\left(p\times p\right)$ identity matrix, ${\tilde{x}}_t=x_t-\sum _{i=1}^n{W}_{\lambda t}\left(z_i\right)x_i=\left\{\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{X}=\tilde{\mathbf{X}}\right\}$, ${\tilde{y}}_t=y_t-\sum _{i=1}^n{W}_{\lambda t}\left(z_i\right)y_i=\left\{\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}=\tilde{\mathbf{y}}\right\}$, and ${\mathbf{W}}_{\lambda }$ is a kernel smoother matrix of weights satisfying the satisfying Assumption (2.3), given by

$W_{\lambda t}\left(z\right)=K\left({\displaystyle \frac{z_t-z}{\lambda }}\right)/ \sum _{t=1}^{n}K\left({\displaystyle \frac{z_t-z}{\lambda }}\right)={\mathbf{W}}_{\lambda } \left(2.4\right)$where $\lambda$ is a bandwidth parameter to selected, and $K\left(.\right)$ is kernel function defined in Remark 2.2. For example, $K\left(.\right)$ could be Gaussian, Epanechnikov, etc. Then, an estimator $\widehat{f}$ of the nonparametric part $f\left(z\right)$ is obtained by

${\widehat{\mathbf{f}}}_{RTK}=\sum _{i=1}^n{W}_{\lambda t}\left(z\right)\left(y_t-{\mathbf{x}}_t{\widehat{\beta }}_{RTK}\right)={\mathbf{W}}_{\lambda }\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_{RTK}\right)$ (2.5)

Hence, the autoregressive coefficient $\widehat{\rho }$arises from (2.2) and the noise of the AR(1) is obtained by ${\widehat{u}}_t=e_t-\widehat{\rho }{e}_{t-1}$. Note that since the errors in the model (1.1) are serially correlated, the estimators defined in (2.3) and (2.5) is not asymptotically efficient. To improve efficiency, we use kernel ridge type weighted estimation based on transforming data $\tilde{\mathbf{y}}$.

3. Generalized Ridge Type Kernel Estimation

The generalized ridge type kernel (GRTK) estimator of the parametric components in the model (1.1) is constructed by combining the partial kernel method with a suitable transformation to account for $\rho$ in the AR(1) errors. We will then extend this GRTK estimator to incorporate shrinkage estimation, analyze its asymptotic properties, and discuss its statistical characteristics. For simplicity, assume the true parameter vector $\mathsf{\beta }$ and true error covariance matrix $\mathbf{R}$ are given. From the ${E(\epsilon }_t)=0$, one can see that $f\left(z_t\right)=y_t-{\mathbf{x}}_t\mathsf{\beta }, t=1,\dots ,n$. Thus, for a given $\mathsf{\beta }$, the natural estimator of the nonparametric component $f\left(.\right)$ is

$\widehat{\mathbf{f}}=\sum _{t=1}^n{W}_{\lambda t}\left(z\right)\left(y_t-{\mathbf{x}}_t\widehat{\mathsf{\beta }}\right)={\mathbf{W}}_{\lambda }\left(\mathbf{y}-\mathbf{X}\widehat{\mathsf{\beta }}\right) \left(3.1\right)$where ${\mathbf{W}}_{\lambda }$ is a kernel smoother matrix, as defined in (2.4). Since $\mathbf{R}$ defined in Assumption (2.1) is positive definite, there exists a $\left(n\times n\right)-$ dimensional matrix $\mathbf{P}$ such that$\mathbf{P}\mathbf{\text{'}}\mathbf{P}={\mathbf{R}}^{-1}$. We first assume that the matrix $\mathbf{R}$ and, hence, $\mathbf{P}$ are known. Then, we fit the error structure by the following residuals:

$e_t=y_t-{\mathbf{x}}_t\mathsf{\beta }-\widehat{f}={\tilde{y}}_t-{\tilde{\mathbf{x}}}_t\mathsf{\beta },1\le t\le n$ (3.2)

Here, as defined in (2.3), ${\tilde{y}}_t$ and $\tilde{\mathbf{x}}$ are locally centered residuals of $y_t$ and $\mathbf{x}$, respectively. By considering these partial residuals, the GRTK estimator of the vector $\mathsf{\beta }$ is obtained by minimizing the following weighted least squares (WLS) equation.

$WLS\left(\mathsf{\beta }\right)={\left(\mathbf{P}\left({\tilde{y}}_t-{\tilde{\mathbf{x}}}_t\mathsf{\beta }\right)\right)}^{\mathit{\text{'}}}\left(\mathbf{P}\left({\tilde{y}}_t-{\tilde{\mathbf{x}}}_t\mathsf{\beta }\right)\right)+k{\mathsf{\beta }}^{\mathbf{\text{'}}}\mathsf{\beta }$$={‖\mathbf{P}\left(\tilde{\mathbf{y}}-\tilde{\mathbf{X}}\mathsf{\beta }\right)‖}_2^2+k{‖\mathsf{\beta }‖}_2^2$ (3.3)

where ${‖.‖}_2$ denotes the Euclidean norm, $k>0$ is a shrinkage parameter to be selected by a selection method, $\tilde{\mathbf{y}}=\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}$ and $\tilde{\mathbf{X}}=\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{X}$ are the centered residuals, as defined in (2.3). Also, ${\mathbf{W}}_{\lambda }$ is a kernel smoothing matrix based on a bandwidth parameter $\lambda >0$ to be selected.

Algebraically, after some operations, the solution of the equation $WLS\left(\mathsf{\beta }\right)$ yields the following GRTK estimator:

${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathbf{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathbf{R}}^{-1}\tilde{\mathbf{y}}$ (3.4)

Replacing $\mathsf{\beta }$ in the (3.1) with the ${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)$ in (3.4) to produce an GRTK estimator of the nonparametric component in the partially linear model $f\left(.\right)$ given by

${\widehat{\mathbf{f}}}_{GRTK}=\sum _{t=1}^n{W}_{\lambda t}\left(z\right)\left(y_t-{\mathbf{x}}_t{\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right)={\mathbf{W}}_{\lambda }\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right) \left(3.5\right)$This section details the construction of the Generalized Ridge-Type Kernel (GRTK) estimator for the parametric component in model (1.1). This is achieved by combining the partial kernel method with a suitable transformation that accounts for $\rho$ in the AR(1) errors. We will then extend this GRTK estimator to incorporate shrinkage estimation, analyze its asymptotic properties, and discuss its statistical characteristics.

3.1. Shrinkage Estimators with GRTK Estimators

While the GRTK estimator in (3.4) addresses multicollinearity and autocorrelation simultaneously, further variance reduction is possible via Stein-type shrinkage. By combining the full model GRTK estimator with a more parsimonious submodel estimator, we can develop Stein-type shrinkage estimators that may achieve lower risk under certain conditions. This approach follows the framework developed by [3] and extends it to the partially linear time-series context. To develop the shrinkage estimators, we consider two models:

1. Full Model (FM): The complete model with all predictors $\left(p=p_1+p_2\right)$ where $p_1<p$ is the number of significant coefficients and $p_2$ denotes the number of sparse coefficients, estimated using GRTK as in equation (3.4).

2. Submodel (SM): A reduced model with dimension ($p_1<p$), selected via Bayesian information criterion (BIC) to minimize model selection, yielding ${\widehat{\mathsf{\beta }}}_{SM}$.

Let ${\widehat{\mathsf{\beta }}}_{FM}$ be the parametric estimate from the ridge-type kernel smoothing (GRTK) on all predictors as full model ${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)$ from equation (3.4) and let ${\widehat{\mathsf{\beta }}}_{SM}$​ be the corresponding estimate when only a submodel of predictors is retained using BIC. Accordingly, ordinary shrinkage estimator (${\widehat{\mathsf{\beta }}}_S$) can be give as follows:

${\widehat{\mathsf{\beta }}}_S={\widehat{\mathsf{\beta }}}_{FM}-\left[(p_2-2)/\widehat{T}\right]({\widehat{\mathsf{\beta }}}_{FM}-{\widehat{\mathsf{\beta }}}_{SM}),$ (3.6)

where $\widehat{T}$is a “distance measure” given below, and $(p_2-2)/\widehat{T}$​ parallels the classical shrinkage factor. Intuitively, we are pulling ${\widehat{\mathsf{\beta }}}_{FM}$​ toward ${\widehat{\mathsf{\beta }}}_{SM}$​, thus introducing some bias but potentially reducing variance and risk. The distance is given by:

$\widehat{T}=n{\widehat{\sigma }}^{-2}{\widehat{\mathsf{\beta }}}_2^{\text{'}}\left[{\tilde{\mathbf{X}}}_2^{\text{'}}{\tilde{\mathbf{U}}}_1{\tilde{\mathbf{X}}}_2\right]{\widehat{\mathsf{\beta }}}_{2 },$where ${\widehat{\mathsf{\beta }}}_2^{\text{'}}$ is the portion of ${\widehat{\mathsf{\beta }}}_{FM}$ (or ${\widehat{\mathsf{\beta }}}_{SM}$) corresponding to the submodel’s parametric indices. ${\tilde{\mathbf{X}}}_2$ is the part of the design matrix associated with those sparse indices. ${\tilde{\mathbf{U}}}_1={\mathbf{I}}_n-{\tilde{\mathbf{X}}}_1{\left[{\tilde{\mathbf{X}}}_1^{\text{'}}{\tilde{\mathbf{X}}}_1+k{\mathbf{I}}_{p_1}\right]}^{-1}{\tilde{\mathbf{X}}}_1^{\text{'}}$ is a partial residual projection, removing the other block of parametric covariates from the fit. ${\widehat{\sigma }}^2$ is given in (3.17) but note that ${\mathbf{H}}_{\lambda ,k}$ is chosen SM-based partial kernel fit (see [21] for details).

To prevent over-shrinking when $(p_2-2)/\widehat{T}$ is large or negative, we define a positive-part version:

${\widehat{\mathsf{\beta }}}_{PS}={\widehat{\mathsf{\beta }}}_{FM}-{\left[{\displaystyle \frac{p_2-2}{\widehat{T}}}\right]}^{+} \left({\widehat{\mathsf{\beta }}}_{FM}-{\widehat{\mathsf{\beta }}}_{SM}\right),$ (3.7)

where ${\left(x\right)}^{+}=\max \left\{0,x\right\}.$ Thus, if $\left(p_2-2\right)/\widehat{T}$ is negative, we do no shrinkage; if it is positive, we shrink as in (3.6).

Following a partially linear kernel-smoothing approach (GRTK), the nonparametric function $f\left(.\right)$ can be recomputed (or simply computed once) based on these final shrunk parametric estimates using smoothing matrix ${\mathbf{W}}_{\lambda }$:

${\widehat{\mathbf{f}}}_S={\mathbf{S}}_{\lambda }\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_S\right)$, (3.8a)

and

${\widehat{\mathbf{f}}}_{PS}={\mathbf{S}}_{\lambda }\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_{PS}\right)$, (3.8b)

Hence, each final semiparametric estimators become $\left({\widehat{\mathsf{\beta }}}_S,{\widehat{\mathbf{f}}}_s\right), \left({\widehat{\mathsf{\beta }}}_{PS},{\widehat{\mathbf{f}}}_{Ps}\right)$.

3.2. Asymptotic Distribution of the Estimators

This subsection presents the assumptions and theorems necessary to analyze the asymptotic properties of the proposed ${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)$ and ${\widehat{f}}_{GRTK}\left(k\right)$ estimators. Understanding the asymptotic distribution is crucial for statistical inference. In this context, we introduce the following assumptions that can be easily satisfied.

Assumption 3.1. As in setting of the semiparametric regression model, covariates $x_i$ and $z_i$ are related via the following nonparametric regression model.

$x_{ij}=g_j\left(z_i\right)+{\eta }_{ij} ,(1\le i\le n, 1\le j\le p)$ (a)

where ${\eta }_{ij}=({\eta }_{i1},\dots ,{\eta }_{ip})\text{'}$ are real sequence satisfying.

$\underset{n\to \infty }{\lim }{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\eta }_i{\eta }_i^{\text{'}}=\mathit{B} \left(\mathrm{b}\right)$and

${max}_{1\le j\le n}‖\sum _{i=1}^n{W}_{\lambda i}\left(z_i\right){\eta }_i‖=O\left(n^{-1/6}\right)$ (c)

where $\mathbf{B}$ is a $\left(p\times p\right)$ positive definite matrix and $‖.‖$ indicates the Euclidean norm.

Assumption 3.2. The functions $\mathrm{g}\left(.\right)$ and $f_j\left(.\right), 1\le j\le k,$ satisfy a Lipschitz continuous of order 1 on data domain.

Assumption 3.3. The weight functions $W_{\lambda i}$ satisfy these conditions:

${max}_{1\le i\le n}\sum _{j=1}^n{W}_{\lambda i}\left(z_j\right)=O\left(1\right)$${max}_{1\le i,j\le n}{W}_{\lambda i}\left(z_j\right)=O\left(n^{-2/3}\right)$${max}_{1\le j\le n}\sum _{i=1}^n{W}_{\lambda i}\left(z_j\right)I\left(\left|z_i-z_j\right|>O\left(n^{-1/3}\right)\right)=O\left(n^{-1/3}\right),$where $I\left(A\right)$ is the indicator function of a set $A$.

$tr\left(W_{\lambda }^{\text{'}}{W}_{\lambda }\right)=tr\left(W_{\lambda }\right)=O\left({\lambda }^{-1}\right)$, where $tr\left(A\right)$ denotes the trace of a square matrix $A.$Remark 3.1: The ${\eta }_{ij}$ stated in (a) of Assumption 3.1 behaves like as a zero mean uncorrelated sequence of stationary random variables of independent of ${\epsilon }_t$. If the covariates $x_i$ and $z_i$ are the observations of the independent and identically distributed random variables, then $g_j\left(z\right)$ can be considered as $E\left(x_{ij}|z_i\right)=g_j\left(z_i\right)$ $, E({\eta }_i{\eta }_i^{\text{'}})=\mathit{B},\mathrm{a}\mathrm{n}\mathrm{d} {\eta }_{ij}=\left(x_{ij}-g_j\left(z_i\right)\right)$. See [14] for more details. In this case, (b) holds with probability 1 according to the law of large numbers, and (c) is provided by Lemma 1 in [16].

Remark 3.2: If $K\left(x\right)$ is a kernel function, then $K\left(x\right)=K\left({\displaystyle \frac{x}{\lambda }}\right)$ is also a kernel function based on a positive bandwidth parameter $\lambda$, and it satisfies the following properties:

$$\int K\left(x\right)dx=1$$

(a)    ,

$$\begin{gathered} \int xK\left(x\right)dx=0, \\ {\mu }_2=\int x^{2}K\left(x\right)dx<\infty \end{gathered}$$

(b)  

(3.9)

(c) ,

$$K\left(x\right)\ge 0 \mathrm{for} \mathrm{all} x, \mathrm{and} K\left(x\right)=K\left(-x\right)$$

(d)

These properties show that a kernel function needs to be symmetric and continuous probability density function with mean zero and constant variance. Note that the bandwidth parameter $\lambda$ should be chosen optimally by a selection criterion in kernel estimation or smoothing. For example, a large $\lambda$ provides an extremely smooth curve or estimate, while a small $\lambda$ produces a wiggly function curve.

Theorem 3.1: Suppose that the Assumption 2.1 and Assumptions 3.1–3.3 hold, and assume that $\underset{n\to \infty }{\lim }{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\eta }_i{\mathit{R}}^{-1}{\eta }_i^{\text{'}}=\mathit{V}$ is a non-singular covariance matrix with ${\eta }_i={({\eta }_1,\dots ,{\eta }_n)}^{\text{'}}$. Then, as $n\to \infty ,$ we have

$\sqrt{n}\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)-\mathsf{\beta }\right)\stackrel{D}{\to }N\left(0,{\sigma }^2{\mathit{V}}^{-1}\right)$ and $\sqrt{n}(\widehat{\rho }-\rho )\stackrel{D}{\to }N\left(0, n^{-1}\left(1-{\rho }^2\right)\right)$where $\stackrel{D}{\to }$denotes convergence in distribution. See [18] for the proof of the Theorem 3.1.

The asymptotic normality established in Theorem 3.1 allows for valid statistical inference based on the GRTK estimator. Intuitively, this result holds because the transformed data approach effectively accounts for autocorrelation, while the kernel smoothing effectively separates the nonparametric component, allowing the parametric component to be estimated consistently with a regular $\sqrt{\left(n\right)}$ convergence rate.

Theorem 3.2: If the Assumptions 2.1 and 3.1–3.3 hold, then it hold

${max}_{1\le i\le n}\left|{\widehat{f}}_{GRTK}\left(z_i\right)-f\left(z_i\right)\right|=O\left(n^{-1/3} \mathrm{l}\mathrm{o}\mathrm{g}\left(n\right)\right)$a.s.

See [18] for proof of the Theorem 3.2. The asymptotic normality and convergence results established in Theorems 3.1 and 3.2 provide a theoretical foundation for the GRTK estimators. In the following section, we analyze additional statistical properties, such as bias and variance.

This convergence rate for the nonparametric component is optimal in the sense that no estimator can achieve a faster uniform convergence rate without additional structural assumptions. The logarithmic factor arises from the need to establish uniform rather than pointwise convergence.

Note that Theorem 3.2 express that ${\widehat{f}}_{GRTK}$ reaches the optimal strong convergence rate. Also, a consistent estimator of the asymptotic covariance matrix $\mathit{V}$ is required for statistical inference based on ${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)$. The estimate of covariance matrix$\mathit{V}$ can be obtained as

$\widehat{\mathit{V}}={\displaystyle \frac{1}{n}}\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}\right) \left(3.10\right)$where $\widehat{\mathbf{V}}$ is the asymptotic covariance matrix of $\sqrt{n}\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)-\mathsf{\beta }\right)$. Regarding the shrinkage estimator’s asymptotic properties for the currently considered low-dimensional ($p<n$) settings, we outline the key asymptotic results. Let $p_2\ge 3$ be the size of the submodel portion. Suppose ${\beta }_2$ is subject to a local alternatives framework:

${\beta }_{2n}={\displaystyle \frac{1}{\sqrt{n}}}\vartheta ,$for fixed $\vartheta \in {\mathbb{R}}^{p_2},$So that ${\beta }_{2n}\to 0$ as $n\to \infty$. Additionally, let ${\beta }_{1n}$ denotes the complementary block (nonzero coefficients) of dimension $p_1$. Under standard regularity conditions ([4,21]) we can claim that:

Both ${\widehat{\mathsf{\beta }}}_S$ and ${\widehat{\mathsf{\beta }}}_{PS}$ remain consistent for $\mathsf{\beta }$ as a point of consistency

By construction, these estimators introduce shrinkage-based bias in exchange for variance reduction. The bias depends on $\left(p_2-2\right)/\widehat{T}$ or ${\left[{\displaystyle \frac{\left(p_2-2\right)}{\widehat{T}}}\right]}^{+}$ and and generally involves certain noncentral ${\chi }^2$-based expectations.

Regarding the asymptotic quadratic bias (AQDB), asymptotic covariance and asymptotic distributional risk (ADR), similar expansions hold showing how ${\widehat{\mathsf{\beta }}}_S$ and ${\widehat{\mathsf{\beta }}}_{PS}$ incorporate both ridge (FM) plus Stein corrections.

Exact closed-form expressions for the bias, covariance, and risk match those in [21] or [3], with minor notational changes for $\widehat{T}$ and the partial kernel smoothing. Hence, we omit rewriting them here. However, for completeness, detailed bias and risk expansions for these estimators appear in Appendix B. The main takeaway is that ${\widehat{\mathsf{\beta }}}_{PS}$​ typically achieves the smallest asymptotic distributional risk among the four choices (FM, SM, ordinary Stein, positive-part Stein) whenever $\left(p_2-2\right)/\widehat{T}$ would become negative or large, thus showing the advantage of positive-part truncation in a partially linear framework.

3.3. Statistical Properties of the Estimators

Here, we detail the statistical properties of the GRTK estimator, including its bias, variance, and expected value. These properties help to characterize the estimator’s behavior and quality. As defined in (3.4), when $k=0,$ the GRTK estimators of the parametric and nonparametric components reduce to ordinary kernel smoothing (KS) estimators for a partially linear model with correlated error. They can be defined, respectively, as follows:

${\widehat{\mathsf{\beta }}}_{KS}={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{y}}$ and ${\widehat{\mathbf{f}}}_{KS}{=\mathbf{W}}_{\lambda }\left(\mathbf{y}-\mathbf{X}{\widehat{\mathsf{\beta }}}_{KS}\right)$ (3.11)

Using the abbreviation

${G_k=\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}$Expected value, bias and variance of the estimator ${\widehat{\mathsf{\beta }}}_{GRTK}$ can be defined, respectively, as

$E\left[{\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right]=G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}} E\left[\mathsf{\beta }\right]+G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}$$=\mathsf{\beta }-k G_k\mathsf{\beta }+G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}$ (3.12)

$Bias\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right)-\mathsf{\beta }=G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}-k G_k\mathsf{\beta }$ (3.13)

$Var\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right)={{\sigma }^{2}G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}{(\mathbf{I}-{\mathbf{W}}_{\lambda })}^2{\mathit{R}}^{-1}\tilde{\mathbf{X}} G_k$ (3.14)

The implementation details of Equations (3.12) - (3.14) are given in Appendix A1.

Clearly,$E\left[{\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right]\ne \mathsf{\beta }$ for any $>0$ . Hence, the GRTK estimator is a biased estimator. From the expression above it is clear that the expectation of the GRTK estimator vanishes as $k$ tends to infinity

$\underset{k\to \infty }{\lim }{\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}\mathsf{\beta }+{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}=0$Hence, all coefficients of the parametric component are shrunken towards zero as the ridge (or shrinkage) parameter increases.

As shown in the (3.14), although the smoothing method provides the estimates of the components in the model, they do not directly provide an estimate of the variance of the error terms (i.e. ${\sigma }^2$). In a general partially linear model, the estimate of variance ${\sigma }^2$ can be found by the residual sum of squares

$RSS=\left(\mathbf{y}-\widehat{\mathbf{y}}\right)\mathit{\text{'}}\left(\mathbf{y}-\widehat{\mathbf{y}}\right)=\left(\mathbf{y}-\left(\mathbf{X}\widehat{\mathsf{\beta }}+\widehat{\mathbf{f}}\right)\right)\mathit{\text{'}}\left(\mathbf{y}-\left(\mathbf{X}\widehat{\mathsf{\beta }}+\widehat{\mathbf{f}}\right)\right)$$=\left(\mathbf{y}-{\mathbf{H}}_{\lambda ,\mathit{k}}\mathbf{y}\right)\mathit{\text{'}}\left(\mathbf{y}-{\mathbf{H}}_{\mathit{\lambda },k}\mathbf{y}\right)={‖\left({\mathbf{I}-\mathbf{H}}_{\lambda ,k}\right)\mathbf{y}‖}_2^2$ (3.15)

where ${\mathbf{H}}_{\lambda ,k}={\mathbf{W}}_{\lambda }+{(\mathbf{I}-{\mathbf{W}}_{\lambda })}^2\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}$ is a hat matrix which depends on a smoothing parameter $\lambda >0$ and a shrinkage parameter $k$ > 0 for partially linear regression model. Note that the fitted values of the model defined in the (1.1) is obtained by this matrix ${\mathbf{H}}_{\lambda ,k}$ , given by

$\widehat{\mathbf{y}}=\mathbf{X}\widehat{\mathsf{\beta }}+\widehat{\mathbf{f}}=\left[{\mathbf{W}}_{\lambda }+{(\mathbf{I}-{\mathbf{W}}_{\lambda })}^2\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\right]\mathbf{y}={\mathbf{H}}_{\lambda ,k}\mathbf{y}$ (3.16)

where ${G_k=\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}$, as defined in equations (3.12-3.14). The implementation details of Equation (3.16) is given in Appendix A2. Thus, similar to ordinary least squares regression, estimation of the error variance can be stated as

${\widehat{\sigma }}^2={\displaystyle \frac{{‖\left({\mathbf{I}-\mathbf{H}}_{\lambda ,k}\right)\mathbf{y}‖}_2^2}{tr{(\mathrm{I}-{\mathbf{H}}_{\lambda ,k})}^2}}={\displaystyle \frac{RSS}{n-p}} \left(3.17\right)$where $tr{(\mathrm{I}-{\mathbf{H}}_{\lambda ,k})}^2=n-tr\left[2{\mathbf{H}}_{\lambda ,k}-\left({\mathbf{H}}_{\lambda ,k}\right)\mathit{\text{'}}\left({\mathbf{H}}_{\lambda ,k}\right)\right]=n-p$ denotes the residual degrees of freedom (see, for example, [5]). It appears from the denominator of equation (3.17) that the degrees of freedom for $RSS$ can also be expressed as the sample size minus the number of estimated parameters in the model.

3.4. Assessing the Risk and Efficiency

To further evaluate the estimators, we now consider measures of risk and efficiency. As discussed earlier, $Var\left(\widehat{\mathsf{\beta }}\right)$ and $Bias\left(\widehat{\mathsf{\beta }}\right)$ important criteria for assessing estimator quality. This section introduces the Mean Squared Deviation (MSD) and quadratic risk function to compare the performance of different estimators. In general, the ill-effect of $Bias\left(\widehat{\mathsf{\beta }}\right)$ is known as information loss, and this loss can be measured by a function stated as the mean squared deviation (MSD). Note that the MSD is a risk function corresponding to the expected value of the squared error loss for an estimator. In terms of a squared error loss, the $MSD$ can be defined as a matrix consisting of the sum of the variance-covariance matrix and the squared bias. While the MSD matrix provides comprehensive information about estimation quality, comparing matrices directly is challenging. Therefore, we derive a scalar measure by taking the trace of the MSD matrix, which represents the sum of mean squared errors across all coefficient estimates:

$MSD\left(\widehat{\mathsf{\beta }},\mathsf{\beta }\right)=E\left({\left(\widehat{\mathsf{\beta }}-\mathsf{\beta }\right)}^{\text{'}}\left(\widehat{\mathsf{\beta }}-\mathsf{\beta }\right)\right)=Var\left(\widehat{\mathsf{\beta }}\right)+Bias{\left(\widehat{\mathsf{\beta }}\right)}^2\left(3.18\right)$The equation (3.18) gives detailed information about the quality of an estimator. In addition to the $MSD$ matrix, the expected loss, called the scalar-valued version, can also be used to compare different estimators. For convenience, we will work with the scalar valued mean squared deviation.

Definition 3.1: The quadratic risk function of an estimator $\widehat{\mathsf{\beta }}$ of the vector $\mathsf{\beta }$ is defined as the scaler valued version of the mean squared deviation matrix (SMDE), given by

$R\left(\widehat{\mathsf{\beta }},\mathsf{\beta },\mathit{V}\right)=E\left({\left(\widehat{\mathsf{\beta }}-\mathsf{\beta }\right)}^{\text{'}}\mathit{V}\left(\widehat{\mathsf{\beta }}-\mathsf{\beta }\right)\right)=tr\left(\mathit{V}\left(MSD\left(\widehat{\mathsf{\beta }},\mathsf{\beta }\right)\right)\right)=SMSD(\widehat{\mathsf{\beta }},\mathsf{\beta })$where $\mathit{V}$ is a$\left(p\times p\right)$ symmetric and non-negative definite matrix. Based on the above risk function, we can define the following criterion to compare estimators.

Definition 3.2: Let the vectors ${\widehat{\mathsf{\beta }}}_1$ and ${\widehat{\mathsf{\beta }}}_2$ be the two competing estimators of a parameter vector $\mathsf{\beta }$. If the difference of the matrices $MSD$ is non-negative definite, it can be said that the estimator ${\widehat{\mathsf{\beta }}}_2$ is superior to estimator ${\widehat{\mathsf{\beta }}}_1$, and is given by

$∆\left({\widehat{\mathsf{\beta }}}_1,{\widehat{\mathsf{\beta }}}_2\right)=MSD\left({\widehat{\mathsf{\beta }}}_1,\mathsf{\beta }\right)-MSD\left({\widehat{\mathsf{\beta }}}_2,\mathsf{\beta }\right)\ge 0$Accordingly, we get the following result.

Theorem 3.2 ([22]): Let the vectors ${\widehat{\mathsf{\beta }}}_1$ and ${\widehat{\mathsf{\beta }}}_2$ be the two different estimators of a parameter vector $\mathsf{\beta }$. Therefore, the following two expressions are equivalent

$R\left({\widehat{\mathsf{\beta }}}_1,\mathsf{\beta },\mathit{V}\right)-R\left({\widehat{\mathsf{\beta }}}_2,\mathsf{\beta },\mathit{V}\right)\ge 0$ for all non-negative definite matrices $\mathit{V}$ $MSD\left({\widehat{\mathsf{\beta }}}_1,\mathsf{\beta }\right)-MSD\left({\widehat{\mathsf{\beta }}}_2,\mathsf{\beta }\right)$ is a non-negative definite matrix.

The results of Theorem 3.2 denotes that ${\widehat{\mathsf{\beta }}}_2$ has a smaller $MSD\left({\widehat{\mathsf{\beta }}}_2,\mathsf{\beta }\right)$ than the estimator ${\widehat{\mathsf{\beta }}}_1$ if and only if the $R\left({\widehat{\mathsf{\beta }}}_1,\mathsf{\beta },\mathit{V}\right)$ of ${\widehat{\mathsf{\beta }}}_2$averaging over every quadratic risk function is less than that of the estimator ${\widehat{\mathsf{\beta }}}_1$. Thus, the superiority of ${\widehat{\mathsf{\beta }}}_2$ over ${\widehat{\mathsf{\beta }}}_1$ can be observed by comparing the matrices $MSD$.

Substituting equations (3.13) and (3.14) in equation (3.18), the matrix $MSD$ of the proposed estimator ${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)$ is obtained as

$MSD\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right),\mathsf{\beta }\right)=Var\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right)+Bias{\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)\right)}^2$$=G_k\left({\sigma }^2{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}{\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)}^2 {\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}+{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}-k \mathsf{\beta }\right)}^2\right)G_k$Furthermore, considering Definition 3.1, the quadratic risk function for ${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)$ can be stated as follows:

$SMSD\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right),\mathsf{\beta }\right)=R\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right),\mathsf{\beta },\mathit{V}\right)$$=tr\left[\mathit{V}\left(MSD\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right),\mathsf{\beta }\right)\right)\right]=tr\left[MSD\left({\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right),\mathsf{\beta }\right)\right]$$=tr\left[G_k\left({\sigma }^2{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}{\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)}^2 {\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}+{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}-k \mathsf{\beta }\right)}^{\text{'}}\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}-k \mathsf{\beta }\right)\right)G_k\right]$The theoretical properties established in this section demonstrate that GRTK estimators effectively balance bias and variance while accounting for both multicollinearity and error autocorrelation. However, the practical performance of these estimators depends critically on the selection of appropriate tuning parameters ($\lambda ,k$). In the following section, we address this challenge by introducing several parameter selection criteria.

4. Choosing the Penalty Parameters

This section focuses on selecting the bandwidth parameter $\lambda$ and the shrinkage parameter $k$, both of which are crucial components of the generalized weighted least squares equation (3.3). The goal is to determine the optimal values for $\lambda$ and $k$. To achieve this, we employ parameter selection methods. The parameters λ and k are chosen by minimizing specific criteria. The most widely used selection criteria are summarized below:

Generalized cross- validation: The optimal bandwidth $\lambda$ and shrinkage parameter $k$ can be determined by implementing the $GCV$ score function:

$\left(\widehat{\lambda },\widehat{ k}\right)=\underset{\mathit{\lambda }>0,\mathit{k}>0}{\mathrm{argmin}}\mathit{G}\mathit{C}\mathit{V}(\lambda ,k)=\underset{\mathit{\lambda }>0,\mathit{k}>0}{\mathrm{argmin}}\left({\displaystyle \frac{n^{-1}{‖\left({\mathbf{I}-\mathbf{H}}_{\lambda ,k}\right)\mathbf{y}‖}^2}{{ \left[n^{-1}tr\left({\mathbf{I}-\mathbf{H}}_{\lambda ,k}\right)\right]}^2}}\right) \left(4.1\right)$where ${\mathbf{H}}_{\lambda ,k}={\mathbf{W}}_{\lambda }+{(\mathbf{I}-{\mathbf{W}}_{\lambda })}^2\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{G_k=\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1},$ as is defined in the equations (3.15) and (3.16), is the smoother matrix based on the parameters $\lambda$ and $k$.

Improved Akaike information criterion: To eliminate overfitting when the sample size is relatively small, [9] proposed $AI{C}_c$, an improved version of the Classical Akaike Information Criterion

$\mathit{A}\mathit{I}{\mathit{C}}_{\mathit{c}}\left(\lambda ,k\right)=\underset{\lambda >0, k>0}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\mathit{log}\left({\displaystyle \frac{{‖\left({\mathbf{I}-\mathbf{H}}_{\lambda ,k}\right)\mathbf{y}‖}^2}{n-\mathit{tr}\left({\mathbf{H}}_{\lambda ,k}\right)}}\right)+1+{\displaystyle \frac{2\left\{tr\left({\mathbf{H}}_{\lambda ,k}\right)+1\right\}}{n-tr\left({\mathbf{H}}_{\lambda ,k}\right)-2}}\right\} \left(4.2\right)$Bayesian Information Criterion: Schwarz’s criterion, also known as the $BIC$ is another statistical measure for selection of the penalty parameters. The $BIC$ criterion is expressed as

$\mathit{B}\mathit{I}\mathit{C}(\lambda ,k)=\underset{\lambda >0, k>0}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left({\displaystyle \frac{{‖\left({\mathbf{I}-\mathbf{H}}_{\lambda ,k}\right)\mathbf{y}‖}^2}{n-tr\left({\mathbf{H}}_{\lambda ,k}\right)}}\right)+\left({\displaystyle \frac{log\left(n\right)}{n}}\right)tr\left({\mathbf{H}}_{\lambda ,k}\right)\left(4.3\right)$Risk estimation using classical pilots (RECP): The key idea now is to estimate the risk $R(\mathbf{y},\widehat{\mathbf{y}})$ by plugging-in pilot estimates of $\mathbf{y}$ and ${\sigma }^2$ into (4.4), and choose the $\lambda$ and $k$ that minimizes the criterion (see [12])

$R\left(\mathbf{y},\widehat{\mathbf{y}}\right)=\underset{\lambda >0, k>0}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \frac{1}{n-tr\left({\mathbf{H}}_{\lambda ,k}\right)}}\left\{{‖\left({\mathbf{H}}_{\lambda ,k}-\mathbf{I}\right)\mathbf{y}‖}^2+{\sigma }^{2}tr\left({{\mathbf{H}}_{\lambda ,k}}^{\mathit{\text{'}}}{\mathbf{H}}_{\lambda ,k}\right)\right\} \left(4.4\right)$Each of these four criteria offers a valid approach for parameter selection; however, they have different characteristics in practice. GCV generally provides a good balance between bias and variance across different sample sizes. $AI{C}_c$ is generally more suitable for smaller samples, where overfitting is a concern. BIC tends to produce more parsimonious models and is often preferred when the true model is believed to be sparse. RECP often provides good fits for the nonparametric component but can be more sensitive to the correlation structure. The simulation study in the next section provides further guidance on criterion selection under different data scenarios.

5. Numerical Examples

5.1. Simulation Study

This section presents a Monte Carlo simulation study designed to compare the estimation performances of the introduced kernel-type ridge estimator (GRTK) and shrinkage estimators for model (1.1) in the presence of correlated errors. To generate parametric predictors, we use a sparse model with $p=20$. To introduce multicollinearity, covariates are generated with a specific level of collinearity, denoted as $\delta$, the following equation is used:

$x_{ij}=\sqrt{(1-{\delta }^2)}{d}_{ij}+\delta d_{ip}, i=1,\dots ,n, j=1,\dots ,p$ (5.1)

where p is the number of predictors, $d_{ij}$’s generated as $d_{ij}~N\left({\mu }_d=0,{\sigma }_d^2=1\right)$, and $\delta =\left(0.95, 0.99\right)$ denotes the two correlation levels between the predictors of parametric component. In this context, simulated data sets are generated from the following model,

$y_t={\mathbf{x}}_t^{\text{'}}\mathsf{\beta }+f\left(z_t\right)+{\epsilon }_t, t=\mathrm{1,2},\dots ,n,$ (5.2)

as defined in (1.1). In given model (5.2), $\mathsf{\beta }{=\left(-1,1.5, 2, 5, 0_{p_2}\right)}^{\text{'}}$, ${\mathbf{x}}_t^{\text{'}}$ is constructed by equation (5.1); the regression function $f\left(z_t\right)=6z_t{\sin \left(z_t\right)}^2$ and $z_t=5(t-0.5)/n$; the error terms ${\epsilon }_t$are generated using a first order autoregressive process (that is, ${\epsilon }_t=\rho {\epsilon }_{t-1}{+u}_t$) with $\rho =0.5$ and $u_t~NIID\left(0,{\sigma }_u^2=1\right)$. Using data generation procedure given by equations (5.1) and (5.2), we consider the three different sample sizes that are $n=50, 150$ and $300$ to investigate the performance of the introduced estimators for low, medium and large samples, and each we use 1000 repetition for each simulation combination. The simulation results are provided in the following figures and tables. In addition, Figure 1 shows the generated data to clarify the data generation procedure better.

, sparse coefficients with zero and non-zero ones and autocorrelation function of error terms for $n=150, \delta =0.95$. .

From a computational perspective, the GRTK estimation procedure involves several steps: (1) initial parameter estimation to obtain residuals, (2) estimation of the autocorrelation parameter $\rho$, (3) data transformation, (4) parameter selection for ($\lambda ,k$), and (5) final estimation. The most computationally intensive step is typically the selection of optimal ($\lambda ,k$) pairs, requiring a two-dimensional grid search. In our implementation, the GCV-based estimator was computationally most efficient, followed by BIC, AICc, and RECP where RECP involves a pilot variance estimation, and it is more costly than others.

Before presenting the results, as mentioned in the sections above, selection of bandwidth of kernel smoother ($\lambda$) and the shrinkage (or ridge) parameter $\left(k\right)$ has a critical importance due to its effect on the estimation accuracy. Therefore, for different simulation configurations, optimal pairs of $\left(\widehat{\lambda },\widehat{k}\right)$ are presented in Table 1 and Figure 2. Hence, it is possible to see how these parameters are affected by the multicollinearity level between predictors $\left(\delta \right)$ and sample size $\left(n\right)$. When examining Table 1, the chosen bandwidth and shrinkage parameters, along with the specified four criteria, are evident for all possible simulation configurations. From the values, it can be observed that the values of the pair $(\widehat{\lambda }, \widehat{k})$ increase as the correlation level rises, potentially negatively impacting estimation quality in terms of smoothing. On the other hand, the increase in $\widehat{k}$, particularly when multicollinearity is high, is an expected behavior aimed at avoiding indefinability in the variance-covariance matrix. For large sample sizes, the criteria tend to select lower $\widehat{\lambda },$but higher "$\widehat{k}$" values as a general tendency. Similar selection behavior is observed when $\delta =0.99$. We take the behavior of $\left(\widehat{\lambda },\widehat{k}\right)$ obtained in Table 1 and follow similar roadmap for the shrinkage estimators that are obtained rely on the GRTK estimator.

pairs selected by $AI{C}_c$, $BIC$ and $RECP$ critera for different simulation configurations for GRTK.

$$\delta =0.95$$

$$\begin{gathered} \delta =0.99 \\ n \end{gathered}$$

$$GCV$$

$$AI{C}_c$$

$$BIC$$

$$RECP$$

$$GCV$$

$$AI{C}_c$$

$$BIC$$

$$\begin{gathered} RECP \\ \widehat{\lambda } \end{gathered}$$

5.2. Real Data Example: Airline Delay Dataset

In this section, we consider the airline delay dataset to show the performance of the modified GRTK estimators. The Bureau of Transportation Statistics (BTS), a division of the U.S. Department of Transportation (DOT), monitors the punctuality of domestic flights conducted by major air carriers. Commencing in June 2003, BTS initiated the collection of comprehensive information as a daily time series regarding the reasons behind flight delays. Based on this collected data, we obtain the partially linear time-series model with GRTK estimators along with shrinkage estimators and measure the comparative performance according to the four selection criteria for both shrinkage $\left(k\right)$ and bandwidth $\left(\lambda \right)$ parameters. The dataset involves originally 1048576 data points which is almost impossible to process for the analysis. In this paper, we extracted a representative consecutive series from the data with $n=300$ data points. Notice that we do not try to solve a specific data-driven problem but showing the merits of the introduced semiparametric estimators. Therefore dataset with $n=300$is enough to represent both the modelling procedure and its alignment with the simulation settings. To construct the semiparametric time-series model, we consider the following six variables as explanatory variables for the parametric component of the model that are departure time ($DepTime$), CRS departure time ($CRSDepTime$), actual elapsed time ($ActualElapsedTime$), CRS elapsed time ($CRSElapsedTime$), departure delay ($DepDelay$) and distance ($Distance$).The air time of the aircrafts ($AirTime$) is determined as a nonparametric variable for the model. The reason for that can be seen in Figure 5 with the hypothetical curve which can be counted as evidence for the nonlinear relationship between the response variable arrival delays ($ArrDelay$) and the $AirTime$ variable. Accordingly, the model to be estimated is given by:

${ArrDelay}_t={\sum }_{j=1}^{p=6}{x}_{jt}{\beta }_j+f\left({AirTime}_t\right)+{\epsilon }_t, t=1,\dots ,300 \left(5.3\right)$where $x_{jt}$’s denote the six explanatory covariates defined above, ${\beta }_j\text{'}s$ are the corresponding regression coefficients and ${\epsilon }_t$’s are the autocorrelated error terms as shown in (1.2). Autocorrelation parameter is estimated as $\rho =0.01$ for difference of ${ArrDelay}_t$. In Figure 5, autocorrelation functions are provided for both original (non-stationary) and differenced-${ArrDelay}_t$ series.

variable.

Figure 6, each panel includes important information about the data to describe the time-series and its properties. In panel (a), variance inflation factors (VIFs) of the predictors in the parametric component of the model reveal that five out of six variables have VIF values greater than 10, indicating a significant multicollinearity problem that needs to be addressed. In panel (b), the linear relationship between predictors and the response variable $ArrDelay$ can be observed, which is also evident in panel (c), depicting pairs of variables with correlation levels. Finally, in panel (d), the relationship between $ArrDelay$ and the nonparametric covariate $AirTime$ is illustrated with a hypothetical curve.

Table 3. Outcomes of simulations with bias, variance and SMSD scores for ${\widehat{\mathsf{\beta }}}_{PS}.$

Table 3. Outcomes of simulations with bias, variance and SMSD scores for ${\widehat{\mathsf{\beta }}}_{PS}.$

In the simulation study, 3D plots are generated to illustrate the selection of bandwidth ($\lambda$) and shrinkage parameter ($k$) for the four criteria, as shown in Figure 7. According to process in Figure 7, selected pairs of $\left(\widehat{\lambda },\widehat{k}\right)$ for the criteria are as follows: $\left({\widehat{\lambda }}_{GCV}=1.357,{\widehat{k}}_{GCV}=1.289\right)$, $\left({\widehat{\lambda }}_{AI{C}_c}=1.242,{\widehat{k}}_{AI{C}_c}=2\right)$, $\left({\widehat{\lambda }}_{BIC}=1.128,{\widehat{k}}_{BIC}=2\right)$ and $\left({\widehat{\lambda }}_{RECP}=0.10, {\widehat{k}}_{RECP}=0.01\right)$. After obtaining the optimal tuning parameters, the performances of the GRTK estimators are presented in Table 4. This table includes the model variance, MSE of the nonparametric component estimate, and overall variance of the regression coefficient, calculated by $tr\left(Var\right(\widehat{\mathsf{\beta }}\left)\right)$, where $Var\left(\widehat{\mathsf{\beta }}\right)$ is shown in (3.11). The best scores are indicated with bold colors in Table 4. According to the results, consistent with the findings from the simulation study, for $n=300$, we observe closer performances.

based on the four criteria $GCV$, $AI{C}_c$, $BIC$ and $RECP$ for airline delay dataset.

Simultaneously, the $RECP$-based estimator shows good performance, closely followed by the $GCV$-based estimator. $AI{C}_c$ and $BIC$-based estimators also exhibit considerable performances and are close to the other two. Additionally, to assess the statistical significance of the estimated regression coefficients, Student t-test statistics and their p-values are presented in Table 5. Furthermore, the estimated models are tested using the F-test, and p-values are provided. The results indicate that all estimated models based on the four criteria are statistically significant ($p<0.05$). However, ${\widehat{\beta }}_{DepTime}$, ${\widehat{\beta }}_{CRSDepTime}$, and ${\widehat{\beta }}_{Distance}$ do not make a significant contribution to the estimated model, despite the overall significance of the models. In conclusion, it can be stated that all criteria yield meaningful models based on their parameter selection procedures.

Bold color indicates the best scores.

Table 6 presents the final performance measures of the GRTK-based estimators (including their shrinkage variants) for the real airline delay dataset, focusing on both parametric and nonparametric components. These metrics—covering model variance, the MSE of the nonparametric estimate, and the overall variance of the regression coefficients—show that all four selection criteria ($GCV, AICc, BIC, RECP$) yield closely comparable results, mirroring the pattern observed in the simulation study for moderate-to-large sample sizes. In particular, RECP exhibits a slight edge in some measures, yet the other criteria remain competitive. This consistency with the simulation outcomes underscores two key conclusions: (1) the newly introduced GRTK framework and its shrinkage modifications successfully mitigate variance inflation in the presence of multicollinearity; and (2) once an adequately sized dataset is available, different penalty-parameter selection methods converge toward similarly effective solutions in partially linear time-series models.

Figure 8 provides a residual analysis for each of the four GRTK-based models (including shrinkage versions) applied to the airline delay dataset. The top row displays the autocorrelation functions (ACFs) of the residuals, indicating whether significant time dependence remains after fitting. In all four cases, autocorrelation appears modest and tapers off quickly, suggesting that the estimated models have captured most of the serial dependence in the data. The bottom row shows residuals scattered around zero without any pronounced pattern, hinting at an absence of systematic bias or heteroscedasticity. Consequently, the visual diagnostic supports the conclusion that, regardless of the specific selection criterion (GCV, AICc, BIC, RECP), the introduced GRTK estimators adequately address the combined challenges of autocorrelation and multicollinearity in this real-world time-series setting.

Regarding the fitted curves, it presents the estimated curves for the nonparametric component $f\left(z\right)$ under the proposed GRTK framework, incorporating both ordinary Stein shrinkage (S) and positive-part Stein shrinkage (PS), as well as the baseline GRTK. Each panel corresponds to one of the four selection criteria—GCV, AICc, BIC, and RECP—which choose optimal pairs $\left(\widehat{\lambda },\widehat{k}\right)$ to govern bandwidth ($\widehat{\lambda }$) and shrinkage ($\widehat{k}$). Mathematically, all four criteria minimize a penalized objective function (e.g., GCV, AICc, etc.) to balance bias (from a potentially large $\widehat{\lambda }$) and variance (from insufficient shrinkage). As seen in the figure, different choices of $\left(\widehat{\lambda },\widehat{k}\right)$ yield subtle variations in the smoothness and curvature of$f\left(z\right)$. For example, $AI{C}_c$ typically enforces a slightly larger bandwidth to address overfitting risks in moderate samples, resulting in a smoother curve, whereas $BIC$ can pick a smaller $\lambda$ or a larger $k$, leading to more parsimonious fits with sharper inflection points. GCV often strikes a middle ground, and RECP sometimes selects a relatively large bandwidth or smaller shrinkage factor if its pilot-based risk estimation deems this configuration optimal.

Figure 9. Fitted curves obtained based on the four criteria.

Figure 9. Fitted curves obtained based on the four criteria.

Beyond these distinctions, the figure also highlights the effect of adding shrinkage to the baseline GRTK estimates. Visually, the PS-based curves often appear slightly more adaptive in regions where GRTK or S might be overly penalized, reflecting a more balanced bias–variance trade-off. Despite these finer differences, all curves capture the primary nonlinear trends and fall within reasonable bounds, reinforcing the broader conclusion that combining kernel smoothing with either Stein or positive-part Stein shrinkage effectively handles both autocorrelation and multicollinearity in semiparametric time-series models.

6. Conclusions

This paper introduces and analyzes modified ridge-type kernel smoothing estimators tailored for semiparametric time-series models with multicollinearity in the parametric component. By combining generalized ridge-type kernel (GRTK) methodology with shrinkage and positive-part Stein shrinkage versions, we address both near-linear dependencies among regressors and autocorrelation in the error structure. We also explore how to optimally select the bandwidth ($\lambda$) and shrinkage parameter ($k$) using four widely used criteria: $GCV, AI{C}_c, BIC$, and $RECP$.

From the detailed simulation study and the real-world airline delay dataset, the following conclusions can be drawn:

The GCV-based GRTK estimator effectively balances bias and variance for both parametric and nonparametric components. In simulations, it consistently provides stable estimates for linear coefficients $\mathsf{\beta }$ and suitably smooth fits for the nonlinear function $f\left(z\right)$.

All four criteria show instability in small samples, as expected. However, in medium and large samples, the proposed estimators achieve more reliable performance, with reduced bias, variance, and SMSD. The GRTK-based approach effectively manages the autoregressive nature of the errors, highlighting the importance of accounting for correlation in time-series data.

Shrinkage estimators, especially positive-part Stein, excel at mitigating variance inflation and overshrinking under strong multicollinearity. They outperform standard methods when predictors are strongly correlated, particularly in larger samples.

Airline delay data results align with simulations. All selection methods yield comparable models, demonstrating their ability to balance penalty tuning. This advantage becomes most pronounced in moderate and large samples, and GCV often provides a straightforward yet effective method for setting ($\lambda ,k$).

Building on the insights gained here, several cases remain open for further research such as extending to non-stationary time-series data, developing locally adaptive bandwidth selection, incorporating alternative variable selection techniques for high-dimensional data, investigating efficient algorithms for parameter selection.

As a result, the proposed estimators fill an important gap in semiparametric modeling by jointly tackling multicollinearity and autocorrelation, thereby yielding more robust and interpretable results for partially linear models in time-series contexts.

Appendix A1. Derivations of Equations (3.12 - 3.14)

${\widehat{\mathsf{\beta }}}_{GRTK}\left(k\right)={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{y}}$ $={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}\right)$$=G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\left(\mathbf{X}\mathsf{\beta }+\mathbf{f}+\mathsf{\epsilon }\right)\right)$$=G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\tilde{\mathbf{X}}\mathsf{\beta }+\tilde{\mathbf{f}}+\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathsf{\epsilon }\right)$$=G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}\mathsf{\beta }+G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}+G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathsf{\epsilon }$Accordingly, equation (3.12) can be expressed as:

$E\left({\widehat{\mathsf{\beta }}}_{GRTK}\right)=G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}} E\left(\mathsf{\beta }\right)+G_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}$$={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}-k\mathbf{I})\mathsf{\beta }+{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}$$={\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}\left[\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)\mathsf{\beta }-k\mathbf{I}\mathsf{\beta }\right]+{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}$$=\left(\mathbf{I}-k {\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}\right)\mathsf{\beta }+{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}$$=\mathsf{\beta }-k {\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}\mathsf{\beta }+{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{f}}$Equivalently, from (3.9), we obtain the bias as in equation (3.13):

$E\left({\widehat{\mathsf{\beta }}}_R\left(k\right)\right)=E\left({\left({\mathbf{I}}_p+k{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}\right)}^{-1}\right)}^{-1}\widehat{\mathsf{\beta }}\right)=E\left[{\left({\mathbf{I}}_p+k{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}\right)}^{-1}\right)}^{-1}{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\left({\mathbf{I}}_n-{\mathbf{W}}_{\lambda }\right)\mathbf{y}\right]$$={\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}\mathsf{\beta }+{\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{f}}\right)$Also, the variance of an estimator ${\widehat{\mathsf{\beta }}}_R\left(k\right)$ given in (3.14) can be derived as follows:

$Var\left(\left({\widehat{\mathsf{\beta }}}_R\left(k\right)\right)\right)=E\left[\left\{\left({\widehat{\mathsf{\beta }}}_R\left(k\right)\right)-E\left({\widehat{\mathsf{\beta }}}_R\left(k\right)\right)\right\}{\left\{\left({\widehat{\mathsf{\beta }}}_R\left(k\right)\right)-E\left({\widehat{\mathsf{\beta }}}_{\mathrm{R}}\left(k\right)\right)\right\}}^{\text{'}}\right]$$=E{\left\{\begin{array}{c}\left[{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}\mathsf{\beta }+{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{f}}+{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\left({\mathbf{I}}_n-{\mathbf{W}}_{\lambda }\right)\mathsf{\epsilon }\right]\\ -\left[{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}\mathsf{\beta }+{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{f}}\right]\end{array}\right\}}^2$$=E\left\{\left({\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\left({\mathbf{I}}_n-{\mathbf{W}}_{\lambda }\right)\mathsf{\epsilon }\right){\left({\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\left({\mathbf{I}}_n-{\mathbf{W}}_{\lambda }\right)\mathsf{\epsilon }\right)}^{\text{'}}\right\}$$=E\left({\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}{\left({\mathbf{I}}_n-{\mathbf{W}}_{\lambda }\right)}^2\right)\left({\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}\right)E\left({\mathsf{\epsilon }}^2\right)$$={\sigma }^2{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\text{'}}{\left({\mathbf{I}}_n-{\mathbf{W}}_{\lambda }\right)}^2{\left({\tilde{\mathbf{X}}}^{\text{'}}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}\tilde{\mathbf{X}}$As a result, it can be expressed as the following result with abbreviation as in (3.14):

$Var\left(\left({\widehat{\mathsf{\beta }}}_R\left(k\right)\right)\right)={\sigma }^2{\mathbf{G}}_k{\tilde{\mathbf{X}}}^{\text{'}}{\left({\mathbf{I}}_n-{\mathbf{W}}_{\lambda }\right)}^2\tilde{\mathbf{X}}{\mathbf{G}}_k$Hence, derivation has been completed.

Appendix A2. The Proof of the Equation (3.16)

$\widehat{\mathbf{y}}=\mathbf{X}\widehat{\mathsf{\beta }}+\widehat{\mathbf{f}}=\mathbf{X}{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{y}}+{\mathbf{W}}_{\lambda }\left(\mathbf{y}-\mathbf{X}\widehat{\mathsf{\beta }}\right)$$=\mathbf{X}{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}+{\mathbf{W}}_{\lambda }\left(\mathbf{y}-\mathbf{X}{\left({\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\tilde{\mathbf{X}}+k\mathbf{I}\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}\right)$$=\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}+{\mathbf{W}}_{\lambda }\left(\mathbf{y}-\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}\right)$$=\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}+{\mathbf{W}}_{\lambda }\mathbf{y}-{\mathbf{W}}_{\lambda }\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)$$={\mathbf{W}}_{\lambda }\mathbf{y}+\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}-{\mathbf{W}}_{\lambda }\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}$$={\mathbf{W}}_{\lambda }\mathbf{y}+\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{y}$$=\left[{\mathbf{W}}_{\lambda }+\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)\right]\mathbf{y}$Thus, the hat matrix is written as follows:

${\mathbf{H}}_{\mathit{\lambda },\mathit{k}}={\mathbf{W}}_{\lambda }+(\mathbf{I}-{\mathbf{W}}_{\lambda })\mathbf{X}{G}_k{\tilde{\mathbf{X}}}^{\mathbf{\text{'}}}{\mathit{R}}^{-1}\left(\mathbf{I}-{\mathbf{W}}_{\lambda }\right)$ Appendix B . Asymptotic Supplement: Bias, AQDB, and ADR

This appendix provides a treatment of the asymptotic bias, Asymptotic Quadratic Distributional Bias (AQDB), and Asymptotic Distributional Risk (ADR) for the GRTK-based estimators introduced in Section 3. Specifically, we focus on:

The baseline GRTK (FM) estimator ${\widehat{\mathit{\beta }}}_{\mathrm{G}\mathrm{R}\mathrm{T}\mathrm{K}}\left(k\right)={\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}$ from equation (3.4),

The ordinary Stein shrinkage estimator ${\widehat{\mathit{\beta }}}_{\mathrm{S}}$ from equation (3.6), and

The positive-part Stein shrinkage estimator ${\widehat{\mathit{\beta }}}_{\mathrm{P}\mathrm{S}}$ from equation (3.7).

with covariance matrix $\mathbf{R}=\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathit{\epsilon }\right)$ capturing possible autocorrelation (e.g., $\mathrm{A}\mathrm{R}\left(1\right)$). Throughout, we assume $\mathbf{R}$ is known (or replaced with a consistent estimator $\widehat{\mathbf{R}}$, which suffices for the same asymptotic properties). We also assume standard regularity conditions on the kernel smoothing for $f(\cdot )$ (bandwidth $\lambda \to 0$ and smoothness conditions as $n\to \mathrm{\infty }$ ). Under these conditions, all estimators in Section 3 are consistent for $\mathsf{\beta }$, and we can derive their bias, variance, and risk expansions.

Appendix B 1. Baseline GRTK (FM) Estimator

Recall that the full-model GRTK estimator may be written as

${\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}={\left({\tilde{\mathbf{X}}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathrm{\top }}{\mathbf{R}}^{-1}\tilde{\mathbf{y}}$where $\tilde{\mathbf{y}}=({\mathbf{I}}_n-{\mathbf{W}}_{\lambda })\mathbf{y}$ is the partially centered response (subtracting off the estimated nonparametric component $\hat{f}(\cdot )$ ), and $\mathbf{X}$ is the $n\times p$ design matrix for the parametric covariates. Let us Define

$\mathbf{L}={\left({\tilde{\mathbf{X}}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}{\tilde{\mathbf{X}}}^{\mathrm{\top }}{\mathbf{R}}^{-1}$Hence ${\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}=\mathbf{L}\tilde{\mathbf{y}}$. If $\tilde{\mathbf{y}}=\tilde{\mathbf{X}}\mathsf{\beta }+\tilde{\mathsf{\epsilon }}$ (the decomposition after partial kernel smoothing of $f(\cdot )$, then

${\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}=\mathbf{L}\tilde{\mathbf{X}}\mathsf{\beta }+\mathbf{L}\tilde{\mathsf{\epsilon }}$One can show that $\mathbf{L}\tilde{\mathbf{X}}$ is effectively ${\mathbf{I}}_p-k{\left({\tilde{\mathbf{X}}}^{\mathrm{\top }}{\mathbf{R}}^{-1}\tilde{\mathbf{X}}+k{\mathbf{I}}_p\right)}^{-1}$, capturing the ridge shrinkage on the parametric coefficients. Hence,

${\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}=\left[{\mathbf{I}}_{\mathrm{p}}-k{\left({\mathbf{X}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\mathbf{X}+\mathrm{k}{\mathbf{I}}_{\mathrm{p}}\right)}^{-1}\right]\mathsf{\beta }+\mathbf{L}\tilde{\mathsf{\epsilon }}.$Exact Bias. Taking expectation (condition on $\mathbf{X}$ ) and noting $E\left[\tilde{\epsilon }\right]\approx 0$ for large $n$, we get

$\begin{array}{rr}\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\left({\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}\right)& =\mathrm{E}\left[{\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}\right]-\mathsf{\beta }\\ & \approx -k{\left({\tilde{\mathbf{X}}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\tilde{\mathbf{X}}+\mathrm{k}{\mathbf{I}}_{\mathrm{p}}\right)}^{-1}\mathsf{\beta }\end{array}$plus small $o\left(1\right)$ terms from the smoothing residual. Since $‖{\mathbf{X}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\mathbf{X}‖=O\left(n\right)$ in typical regressions and $k$ may vanish with $n$ (or remain bounded), this bias goes to zero as $n\to \mathrm{\infty }$. The asymptotic behavior of this bias term depends crucially on how k scales with n. We can distinguish several cases:

If $k=O\left(1\right)$ (remains constant or bounded as n→∞):Since $‖{\mathbf{X}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\mathbf{X}‖=O\left(n\right)$, we have (${\mathbf{X}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\mathbf{X}+k\mathbf{I}{)}^{-1}=\mathrm{O}(1/\mathrm{n})$, and the bias term becomes $k·O(1/n)·\mathsf{\beta }=O\left({\displaystyle \frac{1}{n}}\right)$, which converges to zero at rate $1/n$.

If $k=O\left(nᵅ\right)$ for some $0<\alpha <1$: The bias term becomes $O(nᵅ⁻¹)$, which still converges to zero when α < 1, but at a slower rate than in case (i).

If $k=O\left(n\right)$: The bias may not vanish asymptotically, potentially compromising consistency.

Therefore, for consistency of the ${\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}$ estimator, we require that $k({\mathbf{X}}^{\mathrm{\text{'}}}{\mathbf{R}}^{-1}\mathbf{X}+k\mathbf{I}{)}^{-1}$→ 0 as $n\to \infty$, which is satisfied when $k=o\left(n\right)$. This asymptotic requirement aligns with our parameter selection methods in Section 4, which implicitly balance the bias-variance tradeoff by choosing appropriate $k$ values. Hence the GRTK estimator is consistent and has $O\left(k\right)$ bias in finite $n$, shrinking $\mathsf{\beta }$ towards $0$. In summary:

$\left({\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}\right)=O\left(k\right)\to 0 \mathrm{as} n\to \mathrm{\infty }$ Appendix B 2. Shrinkage Estimators

Denote by ${\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}$ the submodel estimator, obtained by fitting the same GRTK procedure but only to a subset of $p_1<p$ regressors using BIC criterion. In practice, we set the omitted $(p-{\mathrm{p}}_1)$ coefficients to zero. Then the ordinary Stein-type shrinkage estimator is

${\widehat{\mathsf{\beta }}}_{\mathrm{S}}={\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}+\hat{\theta }\left({\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}-{\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}\right)$where $\hat{\theta }\in \left[\mathrm{0,1}\right]$ is a shrinkage factor determined from the data (discussed in Subsection B3). If $\hat{\theta }=1$, no shrinkage is applied (we use the full model); if $\hat{\theta }=0$, we revert to the submodel; for $0<\hat{\theta }<1$, we form a weighted compromise. The positive-part shrinkage estimator is

${\widehat{\mathsf{\beta }}}_{\mathrm{P}\mathrm{S}}={\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}+max\{0,\hat{\theta }\}\left({\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}-{\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}\right)$so that negative values of $\hat{k}$ (which might arise due to sampling error) are replaced by 0 to prevent over-shrinkage. Regarding the bias of shrinkage estimators,

${\widehat{\mathsf{\beta }}}_{\mathrm{S}}-\mathsf{\beta }=\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}-\mathsf{\beta }\right)+\hat{\theta }\left({\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}-{\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}\right)$note that ${\widehat{\mathit{\beta }}}_{\mathrm{S}\mathrm{M}}$ itself has some bias (particularly if the omitted $(p-p_1)$ coefficients are not truly zero), while ${\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}$ has the $O\left(\mathsf{\theta }\right)$ ridge bias but includes all parameters. Let $\left({\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}\right)={\mathbf{b}}_{\mathrm{S}\mathrm{M}}$ and $\left({\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}\right)={\mathbf{b}}_{\mathrm{F}\mathrm{M}}$. A first-order approximation, assuming $\hat{\theta }$ is not strongly correlated with the random errors in ($\left.{\widehat{\mathsf{\beta }}}_{\mathrm{F}\mathrm{M}}-{\widehat{\mathsf{\beta }}}_{\mathrm{S}\mathrm{M}}\right)$, gives

$\mathrm{E}\left[{\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right]-\mathsf{\beta }\approx [1-\mathrm{E}[\hat{\theta }\left]\right]\left({\mathbf{b}}_{\mathrm{S}\mathrm{M}}\right)+\mathrm{E}\left[\hat{\theta }\right]\left({\mathbf{b}}_{\mathrm{F}\mathrm{M}}\right)$Hence the bias of ${\widehat{\mathsf{\beta }}}_{\mathrm{S}}$ is approximately a linear mixture of submodel bias and fullmodel bias. In large $n$, we typically find ${\mathbf{b}}_{\mathrm{F}\mathrm{M}}=O\left(\theta \right)$ small, and ${\mathbf{b}}_{\mathrm{S}\mathrm{M}}$ is either $O\left(\theta \right)$ (when the submodel is correct) or (if some omitted coefficients are actually nonzero) has a bigger $O\left(1\right)$ or $O\left(n^{-1/2}\right)$ term. Thus, ${\widehat{\mathsf{\beta }}}_{\mathrm{S}}$ has a bias that is smaller than the submodel's in cases where $E\left[\hat{\theta }\right]>0$.

Positive-Part Variant. For ${\widehat{\mathsf{\beta }}}_{\mathrm{PS} }$, the same expansion holds but with $\hat{\theta }$ replaced by $\mathrm{m}\mathrm{a}\mathrm{x}\{\hat{\theta },0\}$. This cannot exceed $E\left[\hat{\theta }\right]$, so

$\left({\widehat{\mathsf{\beta }}}_{\mathrm{P}\mathrm{S}}\right)\le \left({\widehat{\mathsf{\beta }}}_{\mathrm{S}}\right)$elementwise (in a typical risk comparison sense). That is, the positive-part always improves or equals the Stein shrinkage in terms of bias magnitude, since it disallows negative $\hat{\theta }$.

As given in Section 3, the distance measure $T$ quantifies how far ${\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}$ is from the submodel ${\widehat{\mathit{\beta }}}_{\mathrm{S}\mathrm{M}}$. A convenient choice is the squared Mahalanobis distance:

$T={\left({\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}-{\widehat{\mathit{\beta }}}_{\mathrm{S}\mathrm{M}}\right)}^{\mathrm{\top }}{\left[\mathrm{C}\mathrm{o}\mathrm{v}\left({\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}-{\widehat{\mathit{\beta }}}_{\mathrm{S}\mathrm{M}}\right)\right]}^{-1}\left({\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}-{\widehat{\mathit{\beta }}}_{\mathrm{S}\mathrm{M}}\right).$Under the "null" submodel assumption (that the omitted parameters are actually zero), $T$ follows approximately a ${\chi }^2$ distribution with $\nu =(p-p_2)$ degrees of freedom, possibly noncentral if the omitted effects are not truly zero. The main text (see eqn. (3.6)) uses $T$ in the formula for the shrinkage factor:

$\hat{\theta }=1-{\displaystyle \frac{\nu -2}{T}}$where $\nu =(p-p_1)$. If $T<\nu -2,\hat{\theta }$ goes negative; the positive-part approach sets ${\hat{\theta }}_{+}=\mathrm{m}\mathrm{a}\mathrm{x}(0,\hat{\theta })$. Asymptotic Distribution of $T$. Under standard conditions (normal or asymptotically normal errors, large $n$ ):

$T\stackrel{d}{\to }{\chi }_{\nu }^2\left(\mathsf{\Lambda }\right),$where $\mathsf{\Lambda }$ is a noncentrality parameter linked to how large the omitted true coefficients are. In a local-alternatives framework with ${\mathsf{\beta }}_2=\mathsf{\delta }/\sqrt{n}$ for the omitted block, $\mathsf{\Lambda }$ remains finite as $n\to \mathrm{\infty }$. If $\mathsf{\Lambda }=0$ (null case), $E\left[T\right]\approx \nu$ and $\left(T\right)\approx 2\nu$. If $\mathsf{\Lambda }>0,T$ tends to be $>\nu$, so $\hat{\theta }\approx 1$. This adaptivity is what drives the shrinkage phenomenon.

If $T$ is large ( $\gg \nu$ ), we suspect omitted coefficients are nonzero, so $\hat{\theta }\approx 1$ retains the full-model estimate.

If $T$ is near $\nu ,\hat{\theta }$is moderate $(0<\hat{\theta }<1)$, partially shrinking FM toward SM.

If $T<\nu -2$, then $\hat{\theta }$ is negative, which is clipped to 0 by the positive-part rule.

Appendix B 3. Asymptotic Quadratic Bias (AQDB) and Distributional Risk (ADR)

We next assess each estimator's risk and define the Asymptotic Quadratic Distributional Bias (AQDB). For an estimator $\widehat{\mathsf{\beta }}$, define

$\mathrm{M}\mathrm{S}\mathrm{E}\left(\widehat{\mathsf{\beta }}\right)=E\left[\Vert \widehat{\mathsf{\beta }}-\mathsf{\beta }{\Vert }^2\right]=Tr\left[\mathrm{V}\mathrm{a}\mathrm{r}\right(\widehat{\mathsf{\beta }}\left)\right]+\Vert \mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\left(\widehat{\mathsf{\beta }}\right){\Vert }^2$In the local alternatives setting (where the omitted block ${\mathsf{\beta }}_2=\mathsf{\delta }/\sqrt{n}$ ), we often look at $n\mathrm{M}\mathrm{S}\mathrm{E}\left(\widehat{\mathsf{\beta }}\right)$ as $n\to \mathrm{\infty }$, which splits into an asymptotic variance component plus an asymptotic (squared) bias:

$\underset{n\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}} nE\left[\Vert \widehat{\mathit{\beta }}-\mathit{\beta }{\Vert }^2\right]=\mathrm{T}\mathrm{r}\left({\mathsf{\Sigma }}_{\hat{\beta }}\right)+{‖{\mathit{\mu }}_{\hat{\beta }}‖}^2$Here ${\mathsf{\Sigma }}_{\hat{\beta }}=\underset{n\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}} \mathrm{V}\mathrm{a}\mathrm{r}\left(\sqrt{n}\widehat{\mathsf{\beta }}\right)$ and ${\mathit{\mu }}_{\hat{\beta }}=\underset{n\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}} \sqrt{n}\left(E\right[\widehat{\mathsf{\beta }}]-\mathsf{\beta })$. We denote:

$\underset{\mathrm{Asymptotic} \mathrm{Variance} }{\underbrace{\mathrm{T}\mathrm{r}\left({\mathsf{\Sigma }}_{\hat{\beta }}\right)}}+\underset{\mathrm{AQDB} }{\underbrace{{‖{\mathit{\mu }}_{\hat{\beta }}‖}^2}}=\mathrm{A}\mathrm{D}\mathrm{R}\left(\widehat{\mathit{\beta }}\right),$where the Asymptotic Quadratic Distributional Bias (AQDB) is ${‖{\mathit{\mu }}_{\hat{\beta }}‖}^2$ and the Asymptotic Distributional Risk (ADR) is their sum.

Full-Model (FM): has negligible bias (so ${\mathit{\mu }}_{\mathrm{F}\mathrm{M}}=0$ ), but a higher variance from estimating all $p$ coefficients. Hence $\mathrm{A}\mathrm{D}\mathrm{R}\left({\widehat{\mathit{\beta }}}_{\mathrm{F}\mathrm{M}}\right)=\mathrm{T}\mathrm{r}\left({\mathsf{\Sigma }}_{\mathrm{F}\mathrm{M}}\right)$, typically $>\mathrm{T}\mathrm{r}\left({\mathsf{\Sigma }}_{\mathrm{S}\mathrm{M}}\right)$.

Submodel (SM): has lower variance (only $p_1$ parameters) but possibly large bias if ${\mathit{\beta }}_2\ne 0$, giving ${\mathit{\mu }}_{\mathrm{S}\mathrm{M}}\ne 0$ and a big ${‖{\mathit{\mu }}_{\mathrm{S}\mathrm{M}}‖}^2$. Specifically, if the omitted block is $\mathit{\delta }/\sqrt{n}$, then ${\mathit{\mu }}_{\mathrm{S}\mathrm{M}}\approx (0,\mathit{\delta })$ so $\mathrm{A}\mathrm{Q}\mathrm{D}\mathrm{B}\left(\mathrm{S}\mathrm{M}\right)=\Vert \mathit{\delta }{\Vert }^2$.

Stein (S): interpolates. Variance is $<\mathrm{V}\mathrm{a}\mathrm{r}\left(\mathrm{F}\mathrm{M}\right)$ but $>\mathrm{V}\mathrm{a}\mathrm{r}\left(\mathrm{S}\mathrm{M}\right)$. Bias is less than the SM's if omitted effects are actually nonzero, but not zero: $(1-E[\hat{\theta }\left]\right)$ times $\mathsf{\delta }$ roughly, so $\mathrm{A}\mathrm{Q}\mathrm{D}\mathrm{B}\left(\mathrm{S}\right)\approx (1-E[\hat{\theta }]{)}^2\Vert \mathsf{\delta }{\Vert }^2$.

Positive-Part Stein (PS): ensures no negative shrinkage, so $\mathrm{A}\mathrm{Q}\mathrm{D}\mathrm{B}\left(\mathrm{P}\mathrm{S}\right)\approx$ ${\left(1-E\left[{\hat{\theta }}_{+}\right]\right)}^2\Vert \mathsf{\delta }{\Vert }^2$ with $E\left[{\hat{\theta }}_{+}\right]\ge E\left[\hat{\theta }\right]$. Consequently, $\mathrm{A}\mathrm{Q}\mathrm{D}\mathrm{B}\left(\mathrm{P}\mathrm{S}\right)\le \mathrm{A}\mathrm{Q}\mathrm{D}\mathrm{B}\left(\mathrm{S}\right)$. Moreover, $\mathrm{V}\mathrm{a}\mathrm{r}\left(\mathrm{P}\mathrm{S}\right)\le \mathrm{V}\mathrm{a}\mathrm{r}\left(\mathrm{S}\right)$ in typical settings, so ${\widehat{\mathsf{\beta }}}_{\mathrm{P}\mathrm{S}}$ dominates ${\widehat{\mathsf{\beta }}}_{\mathrm{S}}$ in ADR.

Thus, positive-part shrinkage is usually recommended when $p\ge 3$ or under moderate correlation/collinearity, as it yields the lowest or nearly lowest ADR across different parameter regimes. This is consistent with the numerical results in Section 5 of the main text and references such as [21] and [3].