Introduction
In the data-rich era, researchers increasingly face two fundamental challenges: variable selection and model specification . The variable selection problem arises from the abundance of available covariates, especially in high-dimensional settings where the number of predictors P P P n n n model specification challenge concerns the choice of an appropriate functional form from an expanding “model zoo.” Parametric models are widely used for their simplicity and interpretability, but they are vulnerable to misspecification. Nonparametric models offer greater flexibility by relaxing functional assumptions, yet they often suffer from the curse of dimensionality , which limits their applicability when the number of covariates is large. Modern machine learning and deep learning methods provide powerful predictive tools, but they typically lack explicit model representations and therefore offer limited interpretability regarding the underlying economic or structural relationships.
A central question in asset pricing is which firm characteristics contain independent information about returns, and in what functional forms they do so. This question is challenging because the literature has documented a “zoo” of characteristics that appear to predict the cross-section of returns, creating both variable-selection and model-specification problems.
In this paper, we study a high-dimensional partially linear semiparametric additive model with interactions, building on LASSO and adaptive group LASSO (Tibshirani 1996, Zou 2006, Yuan & Lin 2006, Wang & Leng 2008). We propose a new specification and estimation framework, Specification LASSO (S-LASSO). Under mild regularity conditions, S-LASSO consistently selects relevant variables and identifies whether they enter linearly, nonlinearly, or through interactions. It therefore addresses variable selection and model specification jointly in high-dimensional environments. In our empirical application, S-LASSO yields a more efficiently specified characteristics-based asset pricing model with better out-of-sample performance.
The method begins with a comprehensive set of candidate components, including linear terms, nonlinear terms, and pairwise interactions. Nonlinear effects are approximated by orthogonal series expansions, interactions by tensor products of basis functions, and the resulting components are selected and estimated through a two-step framework that combines LASSO and adaptive group LASSO.
Our framework is built on a sparse high-dimensional regression structure. Suppose we observe ( Y , P ) (\mathbf{Y}, \mathbf{P}) ( Y , P ) Y \mathbf{Y} Y n × 1 n \times1 n × 1 P \mathbf{P} P n × P n \times P n × P ( X 1 , X 2 , … , X P ) (\mathbf{X}_1, \mathbf{X}_2, \ldots, \mathbf{X}_P) ( X 1 , X 2 , … , X P ) P > n P > n P > n Q \mathbf{Q} Q Q Q Q Q Q Q P P P n → ∞ n \to\infty n → ∞ Q Q Q P \mathbf{P} P
E ( Y i ∣ X 1 i , … , X P i ) = θ + h ( X 1 i , … , X Q i ) , (1) E(Y_i \mid X_{1i},\ldots,X_{Pi})=\theta+h(X_{1i},\ldots,X_{Qi}), \tag*{(1)} E ( Y i ∣ X 1 i , … , X P i ) = θ + h ( X 1 i , … , X Q i ) , (1) where θ \theta θ h ( ⋅ ) h(\cdot) h ( ⋅ ) ϵ i = Y i − E ( Y i ∣ X 1 i , … , X P i ) \epsilon_i=Y_i-E(Y_i\mid X_{1i},\ldots,X_{Pi}) ϵ i = Y i − E ( Y i ∣ X 1 i , … , X P i )
A common starting point is the additive nonparametric model h ( Q i ) = ∑ q = 1 Q f q ( X q i ) h(\mathbf{Q}_i)=\sum_{q=1}^{Q}f_q(X_{qi}) h ( Q i ) = ∑ q = 1 Q f q ( X q i ) f q ( ⋅ ) f_q(\cdot) f q ( ⋅ ) O ( n − 1 / 2 ) O(n^{-1/2}) O ( n − 1/2 )
This motivates the partially linear additive semiparametric model
h ( Q i ) = ∑ l = 1 L β l X i l + ∑ q = L + 1 Q f q ( X i q ) , (2) h(\mathbf{Q}_i)=\sum_{l=1}^{L}\beta_lX_{il}+\sum_{q=L+1}^{Q}f_q(X_{iq}), \tag*{(2)} h ( Q i ) = l = 1 ∑ L β l X i l + q = L + 1 ∑ Q f q ( X i q ) , (2) which separates linear and nonlinear effects. The linear coefficients can then be estimated at the rate O ( n − 1 / 2 ) O(n^{-1/2}) O ( n − 1/2 )
However, these additive models do not account for interactions among covariates. As the following examples show, pairwise interactions are common in economics and finance, so ignoring them can be restrictive.
Example 1.1. In macroeconomics, many production functions incorporate interactions between inputs. For example, the generalized Cobb-Douglas production function can be written as log Y = c + ∑ α = 1 d ∑ β = 1 d c α β log ( X α + X β 2 ) \log Y=c+\sum_{\alpha=1}^{d}\sum_{\beta=1}^{d}c_{\alpha\beta}\log\left(\frac{X_{\alpha}+X_{\beta}}{2}\right) log Y = c + ∑ α = 1 d ∑ β = 1 d c α β log ( 2 X α + X β ) X α X_{\alpha} X α X β X_{\beta} X β c α β c_{\alpha\beta} c α β
Example 1.2. In microeconomics, Deaton & Muellbauer (1980) document a utility model for a household Y Y Y ( X E , X D ) (X_E,X_D) ( X E , X D ) ( X H , X F ) (X_H,X_F) ( X H , X F ) ( X T , X S ) (X_T,X_S) ( X T , X S ) Y = m E D ( X E , X D ) + m H F ( X H , X F ) + m T S ( X T , X S ) + ϵ Y=m_{ED}(X_E,X_D)+m_{HF}(X_H,X_F)+m_{TS}(X_T,X_S)+\epsilon Y = m E D ( X E , X D ) + m H F ( X H , X F ) + m T S ( X T , X S ) + ϵ
Example 1.3. In environmental studies, Dong et al. (2021) examine the effects of C O 2 \mathrm{CO}_2 CO 2
solar irradiance on global sea level rise and specify Y S L = m ( τ , X C O 2 , X S I ) + ϵ Y_{\mathrm{SL}} = m(\tau, X_{\mathrm{CO}_2}, X_{\mathrm{SI}}) + \epsilon Y SL = m ( τ , X CO 2 , X SI ) + ϵ C O 2 \mathrm{CO}_2 CO 2
Example 1.4. In financial studies, Aït-Sahalia et al. (2021) estimate the implied volatility surface using σ ( T , K ) = ∑ j = 0 J ∑ i = 0 L j β i j T i K j \sigma(T, K) = \sum_{j=0}^{J}\sum_{i=0}^{L_j}\beta_{ij}T^iK^j σ ( T , K ) = ∑ j = 0 J ∑ i = 0 L j β ij T i K j K K K T T T
Motivated by these examples, we model the unknown multivariate function as a combination of linear terms, nonlinear terms, and pairwise interactions. Sperlich et al. (2002) emphasize the importance of interactions in additive models and propose marginal integration methods for estimation and testing. Their approach, however, is computationally demanding and restricted to low-dimensional settings.
S-LASSO fills this gap by modeling and selecting pairwise interactions in high-dimensional settings. Pairwise interactions offer a useful compromise between interpretability and flexibility, while the framework can be extended to higher-order interactions if needed. We establish oracle properties for S-LASSO: it consistently selects relevant components and their functional forms, and the estimators of selected components have the same asymptotic distribution as their ordinary least squares counterparts. Monte Carlo experiments confirm S-LASSO’s strong performance even when P > n P > n P > n
Our empirical application studies which firm characteristics explain the cross-section of asset returns, and how they do so. We extend the dataset of Green et al. (2017) and collect 102 firm-level characteristics. Building on this literature, Freyberger et al. (2020) approximate the nonlinear effects of 62 characteristics with B-splines, apply adaptive group LASSO for selection, and introduce interaction terms by multiplying firm size with the remaining characteristics.
We depart from Freyberger et al. (2020) in two ways. First, we allow linear, nonlinear, and interactive effects to coexist. Second, we allow interactions with firm size to take general bivariate forms rather than simple multiplicative ones. Empirically, S-LASSO outperforms the benchmark of Freyberger et al. (2020) both in sample and out of sample, reflecting its ability to separate linear from nonlinear effects and accommodate more flexible size interactions.
Literature Review
Variable selection in high-dimensional settings has been studied more extensively than model specification for unknown multivariate functions. One of the most widely used approaches is LASSO (Least Absolute Shrinkage and Selection Operator). Tibshirani (1996) introduced LASSO to perform regularization and variable selection in linear models:
min α ∑ i = 1 n ( Y i − ∑ j = 1 W α j X i j ) 2 + λ n ∑ j = 1 W ∣ α j ∣ \min_{\alpha} \sum_{i=1}^{n}\left(Y_i-\sum_{j=1}^{W}\alpha_jX_{ij}\right)^2+\lambda_n\sum_{j=1}^{W}|\alpha_j| α min i = 1 ∑ n ( Y i − j = 1 ∑ W α j X ij ) 2 + λ n j = 1 ∑ W ∣ α j ∣ where λ n \lambda_n λ n P ≫ n P \gg n P ≫ n [3 ] and Zou (2006). We revisit this condition in § 4 .
To obtain n \sqrt{n} n
min β ∑ i = 1 n ( Y i − ∑ j = 1 W β j X i j ) 2 + λ n ∑ j = 1 W w ^ j ∣ β j ∣ , \min_{\beta} \sum_{i=1}^{n}\left(Y_i-\sum_{j=1}^{W}\beta_jX_{ij}\right)^2+\lambda_n\sum_{j=1}^{W}\widehat{w}_j|\beta_j|, β min i = 1 ∑ n ( Y i − j = 1 ∑ W β j X ij ) 2 + λ n j = 1 ∑ W w j ∣ β j ∣ , where w ^ j \widehat{w}_j w j w ^ j = ∣ α ^ j ∣ − γ \widehat{w}_j=|\widehat{\alpha}_j|^{-\gamma} w j = ∣ α j ∣ − γ γ > 0 \gamma>0 γ > 0 α ^ j \widehat{\alpha}_j α j
Traditional LASSO methods mainly focus on univariate effects. However, Fan et al. (2016) argue that interactions are also important in many applications. They study Y i = β 0 + ∑ j = 1 d β j X i ( j ) + ∑ j < k β j k X i ( j ) X i ( k ) + ε i Y_i=\beta_0+\sum_{j=1}^{d}\beta_jX_i^{(j)}+\sum_{j<k}\beta_{jk}X_i^{(j)}X_i^{(k)}+\varepsilon_i Y i = β 0 + ∑ j = 1 d β j X i ( j ) + ∑ j < k β j k X i ( j ) X i ( k ) + ε i
Another related method is COSSO , proposed by Lin & Zhang (2006) under the smoothing spline ANOVA framework:
h ( X i ) = b + ∑ j = 1 d f j ( X i ( j ) ) + ∑ j < k f j k ( X i ( j ) , X i ( k ) ) + ⋯ . (3) h(\mathbf{X}_i)=b+\sum_{j=1}^{d}f_j(X_i^{(j)})+\sum_{j<k}f_{jk}(X_i^{(j)},X_i^{(k)})+\cdots. \tag*{(3)} h ( X i ) = b + j = 1 ∑ d f j ( X i ( j ) ) + j < k ∑ f j k ( X i ( j ) , X i ( k ) ) + ⋯ . (3) COSSO is flexible, but it assumes a fixed dimension, does not distinguish linear and nonlinear components, and does not yield an explicit parameterized representation. By contrast, our method accommodates linear terms, performs variable and model selection, and is applicable in high-dimensional settings.
Huang et al. (2010) propose a selection and estimation method for additive nonparametric models based on the group LASSO and adaptive group LASSO. They approximate each unknown function by B-spline bases, f n j ( x ) = ∑ k = 1 m n β j k ϕ k ( x ) f_{nj}(x)=\sum_{k=1}^{m_n}\beta_{jk}\phi_k(x) f nj ( x ) = ∑ k = 1 m n β j k ϕ k ( x )
L n ( μ , β n ) = ∑ i = 1 n ( Y i − μ − ∑ j = 1 P ∑ k = 1 m n β j k ϕ k ( X i j ) ) 2 + λ n ∑ j = 1 P w ^ n j ∥ β n j ∥ 2 , (4) L_n(\mu,\boldsymbol{\beta}_n)=\sum_{i=1}^{n}\left(Y_i-\mu-\sum_{j=1}^{P}\sum_{k=1}^{m_n}\beta_{jk}\phi_k(X_{ij})\right)^2+\lambda_n\sum_{j=1}^{P}\widehat{w}_{nj}\lVert\boldsymbol{\beta}_{nj}\rVert_2, \tag*{(4)} L n ( μ , β n ) = i = 1 ∑ n ( Y i − μ − j = 1 ∑ P k = 1 ∑ m n β j k ϕ k ( X ij ) ) 2 + λ n j = 1 ∑ P w nj ∥ β nj ∥ 2 , (4) where λ n \lambda_n λ n ∥ β n j ∥ 2 \lVert\boldsymbol{\beta}_{nj}\rVert_2 ∥ β nj ∥ 2 β n j = ( β j 1 , … , β j m n ) T \boldsymbol{\beta}_{nj}=(\beta_{j1},\ldots,\beta_{jm_n})^{\mathrm{T}} β nj = ( β j 1 , … , β j m n ) T w ^ n j \widehat{w}_{nj} w nj
The contribution of this paper is threefold. First, we propose a more general and flexible modeling framework that incorporates linear components with faster convergence rates, nonparametric components to capture nonlinear effects, and pairwise interaction terms to allow for richer functional relationships. Second, we develop the S-LASSO method, which simultaneously performs variable selection, model specification, and estimation. We establish its oracle properties in terms of variable selection consistency and show the asymptotic normality of the estimators. Finally, we empirically demonstrate the presence and importance of interactive effects between firm size and other characteristics in characteristics-based asset pricing models.
The rest of the paper is organized as follows. § 2 introduces the model framework. § 3 describes the S-LASSO procedure. § 4 establishes the theoretical results. § 5 reports simulation evidence, and § 6 presents the empirical study. § 7 concludes. All proofs and additional materials are collected in the Appendices.
Model Setup
Partially Linear Additive Semiparametric Model
Suppose that we observe a sample ( Y , P ) (\mathbf{Y},\mathbf{P}) ( Y , P )
Y i = θ + h ( X 1 i , … , X Q i ) + ϵ i , (5) Y_i=\theta+h(X_{1i},\ldots,X_{Qi})+\epsilon_i, \tag*{(5)} Y i = θ + h ( X 1 i , … , X Q i ) + ϵ i , (5) where θ \theta θ ϵ i \epsilon_i ϵ i h ( ⋅ ) h(\cdot) h ( ⋅ ) Y = θ + h ( Q ) + U \mathbf{Y}=\theta+h(\mathbf{Q})+\mathbf{U} Y = θ + h ( Q ) + U θ = θ 1 \theta=\theta\mathbf{1} θ = θ 1 U = ( ϵ 1 , … , ϵ n ) T \mathbf{U}=(\epsilon_1,\ldots,\epsilon_n)^{\mathrm{T}} U = ( ϵ 1 , … , ϵ n ) T h ( Q ) = ( h ( X 11 , … , X Q 1 ) , … , h ( X 1 n , … , X Q n ) ) T h(\mathbf{Q})=(h(X_{11},\ldots,X_{Q1}),\ldots,h(X_{1n},\ldots,X_{Qn}))^{\mathrm{T}} h ( Q ) = ( h ( X 11 , … , X Q 1 ) , … , h ( X 1 n , … , X Q n ) ) T
In this paper, we specify h ( ⋅ ) h(\cdot) h ( ⋅ )
h ( Q ) = ∑ l ∈ L α l X l + ∑ r ∈ R m r ( X r ) + ∑ ( s , s ′ ) ∈ S m s s ′ ( X s , X s ′ ) (6) h(Q) = \sum_{l \in L} \alpha_l X_l + \sum_{r \in R} m_r(X_r) + \sum_{(s,s') \in S} m_{ss'}(X_s, X_{s'}) \tag*{(6)} h ( Q ) = l ∈ L ∑ α l X l + r ∈ R ∑ m r ( X r ) + ( s , s ′ ) ∈ S ∑ m s s ′ ( X s , X s ′ ) (6) where Q Q Q
L = { l : X l has a linear effect on Y } , R = { r : X r has a nonlinear effect on Y } , S = { ( s , s ′ ) : s < s ′ , X s , X s ′ have an interactive effect on Y } . (7) \begin{aligned}
L &= \{l : X_l \text{ has a linear effect on } Y\},\\
R &= \{r : X_r \text{ has a nonlinear effect on } Y\},\\
S &= \{(s,s') : s < s', X_s, X_{s'} \text{ have an interactive effect on } Y\}. \tag*{(7)}
\end{aligned} L R S = { l : X l has a linear effect on Y } , = { r : X r has a nonlinear effect on Y } , = {( s , s ′ ) : s < s ′ , X s , X s ′ have an interactive effect on Y } . (7) The ordering condition s < s ′ s < s' s < s ′
In this setup, m r ( X r ) m_r(X_r) m r ( X r ) m s s ′ ( X s , X s ′ ) m_{ss'}(X_s, X_{s'}) m s s ′ ( X s , X s ′ ) (7) are ∣ L ∣ = L |L| = L ∣ L ∣ = L ∣ R ∣ = R |R| = R ∣ R ∣ = R ∣ S ∣ = S |S| = S ∣ S ∣ = S S S S L C L^C L C R C R^C R C S C S^C S C
The decomposition in Equation (6) mitigates the curse of dimensionality by restricting nonparametric estimation to univariate and bivariate functions while retaining substantial flexibility. Relative to conventional additive models with only univariate components, and relative to Freyberger et al. (2020), which restrict interactions to multiplicative forms involving firm size, our model allows general interactions among covariates and is therefore more flexible. The inclusion of linear components also improves estimation efficiency and reduces computational burden when linear approximations are adequate. Although we focus on pairwise interactions for interpretability and tractability, the framework can be extended to higher-order interactions if needed.
Parameterization
To formalize the method and establish the theoretical results, we introduce notation and define the parameterization. Following Dong et al. (2021), we assume that the j th j^{\text{th}} j th X j \mathbf{X}_j X j n × 1 n \times1 n × 1 X j = ( X j 1 , X j 2 , … , X j n ) T \mathbf{X}_j = (X_{j1}, X_{j2}, \ldots, X_{jn})^{\mathsf{T}} X j = ( X j 1 , X j 2 , … , X j n ) T j = 1 , 2 , … , P j = 1, 2, \ldots, P j = 1 , 2 , … , P Ω n \Omega^n Ω n Ω = R \Omega= \mathbb{R} Ω = R 1
The unknown functions m r : Ω → R m_r:\Omega\to\mathbb{R} m r : Ω → R m a b : Ω 2 → R m_{ab}:\Omega^2 \to\mathbb{R} m ab : Ω 2 → R
m r ( X r ) = ( m r ( X r 1 ) , m r ( X r 2 ) , … , m r ( X r n ) ) T , m a b ( X a , X b ) = ( m a b ( X a 1 , X b 1 ) , m a b ( X a 2 , X b 2 ) , … , m a b ( X a n , X b n ) ) T . (8) \begin{aligned}
\mathbf{m}_r(\mathbf{X}_r) &= \left(m_r(X_{r1}),m_r(X_{r2}),\ldots,m_r(X_{rn})\right)^{\mathsf{T}},\\
\mathbf{m}_{ab}(\mathbf{X}_a,\mathbf{X}_b) &= \left(m_{ab}(X_{a1},X_{b1}),m_{ab}(X_{a2},X_{b2}),\ldots,m_{ab}(X_{an},X_{bn})\right)^{\mathsf{T}}. \tag*{(8)}
\end{aligned} m r ( X r ) m ab ( X a , X b ) = ( m r ( X r 1 ) , m r ( X r 2 ) , … , m r ( X r n ) ) T , = ( m ab ( X a 1 , X b 1 ) , m ab ( X a 2 , X b 2 ) , … , m ab ( X an , X bn ) ) T . (8) We further assume that each m r m_r m r L 2 ( Ω ) L^2(\Omega) L 2 ( Ω ) { ϕ k : k ≥ 0 } \{\phi_k:k\ge0\} { ϕ k : k ≥ 0 } L 2 ( Ω ) = span { ϕ 0 , ϕ 1 , ϕ 2 , … } L^2(\Omega)=\operatorname{span}\{\phi_0,\phi_1,\phi_2,\ldots\} L 2 ( Ω ) = span { ϕ 0 , ϕ 1 , ϕ 2 , … } m r m_r m r m r ( x ) = ∑ k = 0 ∞ b r k ϕ k ( x ) m_r(x)=\sum_{k=0}^{\infty}b_{rk}\phi_k(x) m r ( x ) = ∑ k = 0 ∞ b r k ϕ k ( x ) b r k = ∫ Ω m r ( x ) ϕ k ( x ) d x b_{rk}=\int_{\Omega}m_r(x)\phi_k(x)\,\mathrm{d}x b r k = ∫ Ω m r ( x ) ϕ k ( x ) d x L 2 ( R ) L^2(\mathbb{R}) L 2 ( R ) ϕ i ( x ) = ( π 2 i i ! ) − 1 / 2 H i ( x ) exp ( − x 2 2 ) \phi_i(x)=\left(\sqrt{\pi}2^i i!\right)^{-1/2}H_i(x)\exp\left(-\frac{x^2}{2}\right) ϕ i ( x ) = ( π 2 i i ! ) − 1/2 H i ( x ) exp ( − 2 x 2 ) H i ( x ) H_i(x) H i ( x ) exp ( − x 2 ) \exp(-x^2) exp ( − x 2 )
Introducing a truncation parameter K = K n 2 K=K_n^2 K = K n 2 r r r m r ( X r ) = Φ ( X r ) b r + ξ r \mathbf{m}_r(\mathbf{X}_r)=\Phi(\mathbf{X}_r)\mathbf{b}_r+\boldsymbol{\xi}_r m r ( X r ) = Φ ( X r ) b r + ξ r b r = ( b r 0 , b r 1 , … , b r , K n − 1 ) T \mathbf{b}_r=(b_{r0},b_{r1},\ldots,b_{r,K_n-1})^{\mathsf{T}} b r = ( b r 0 , b r 1 , … , b r , K n − 1 ) T Φ ( X r ) \Phi(\mathbf{X}_r) Φ ( X r )
Φ ( X r ) = [ ϕ 0 ( X r 1 ) ϕ 1 ( X r 1 ) ⋯ ϕ K n − 1 ( X r 1 ) ϕ 0 ( X r 2 ) ϕ 1 ( X r 2 ) ⋯ ϕ K n − 1 ( X r 2 ) ⋮ ⋮ ⋱ ⋮ ϕ 0 ( X r n ) ϕ 1 ( X r n ) ⋯ ϕ K n − 1 ( X r n ) ] . (9) \Phi(\mathbf{X}_r)=
\begin{bmatrix}
\phi_0(X_{r1}) & \phi_1(X_{r1}) & \cdots& \phi_{K_n-1}(X_{r1})\\
\phi_0(X_{r2}) & \phi_1(X_{r2}) & \cdots& \phi_{K_n-1}(X_{r2})\\
\vdots& \vdots& \ddots& \vdots\\
\phi_0(X_{rn}) & \phi_1(X_{rn}) & \cdots& \phi_{K_n-1}(X_{rn})
\end{bmatrix}. \tag*{(9)} Φ ( X r ) = ϕ 0 ( X r 1 ) ϕ 0 ( X r 2 ) ⋮ ϕ 0 ( X r n ) ϕ 1 ( X r 1 ) ϕ 1 ( X r 2 ) ⋮ ϕ 1 ( X r n ) ⋯ ⋯ ⋱ ⋯ ϕ K n − 1 ( X r 1 ) ϕ K n − 1 ( X r 2 ) ⋮ ϕ K n − 1 ( X r n ) . (9) Here, ξ r \boldsymbol{\xi}_r ξ r ξ r = ∑ k = K n ∞ b r k ϕ k ( X r ) \boldsymbol{\xi}_r=\sum_{k=K_n}^{\infty}b_{rk}\phi_k(\mathbf{X}_r) ξ r = ∑ k = K n ∞ b r k ϕ k ( X r )
For the interaction part, we similarly assume that m a b m_{ab} m ab L 2 ( Ω 2 ) L^2(\Omega^2) L 2 ( Ω 2 ) m a b ( x , y ) = ∑ k = 0 ∞ ∑ d = 0 ∞ γ a b k d ϕ k ( x ) ϕ d ( y ) m_{ab}(x,y)=\sum_{k=0}^{\infty}\sum_{d=0}^{\infty}\gamma_{abkd}\phi_k(x)\phi_d(y) m ab ( x , y ) = ∑ k = 0 ∞ ∑ d = 0 ∞ γ abk d ϕ k ( x ) ϕ d ( y ) K = K n K=K_n K = K n
Φ ( X a i ) ⊗ Φ ( X b i ) = ( ϕ 0 ( X a i ) ϕ 0 ( X b i ) , ϕ 0 ( X a i ) ϕ 1 ( X b i ) , … , ϕ 1 ( X a i ) ϕ 0 ( X b i ) , … , ϕ K − 1 ( X a i ) ϕ K − 1 ( X b i ) ) . (10) \begin{aligned}
\Phi(X_{ai})\otimes\Phi(X_{bi})
={}&\bigl(\phi_0(X_{ai})\phi_0(X_{bi}),\phi_0(X_{ai})\phi_1(X_{bi}),\ldots,\\
&\phi_1(X_{ai})\phi_0(X_{bi}),\ldots,\phi_{K-1}(X_{ai})\phi_{K-1}(X_{bi})\bigr). \tag*{(10)}
\end{aligned} Φ ( X ai ) ⊗ Φ ( X bi ) = ( ϕ 0 ( X ai ) ϕ 0 ( X bi ) , ϕ 0 ( X ai ) ϕ 1 ( X bi ) , … , ϕ 1 ( X ai ) ϕ 0 ( X bi ) , … , ϕ K − 1 ( X ai ) ϕ K − 1 ( X bi ) ) . (10) We define the basis matrix as
Φ ( X a ) ⊗ Φ ( X b ) : = [ Φ ( X a 1 ) ⊗ Φ ( X b 1 ) ⋮ Φ ( X a n ) ⊗ Φ ( X b n ) ] . (11) \Phi(\mathbf{X}_a)\otimes\Phi(\mathbf{X}_b):=
\begin{bmatrix}
\Phi(X_{a1})\otimes\Phi(X_{b1})\\
\vdots\\
\Phi(X_{an})\otimes\Phi(X_{bn})
\end{bmatrix}. \tag*{(11)} Φ ( X a ) ⊗ Φ ( X b ) := Φ ( X a 1 ) ⊗ Φ ( X b 1 ) ⋮ Φ ( X an ) ⊗ Φ ( X bn ) . (11) which is the tensor product of the matrices Φ ( X a ) \boldsymbol{\Phi}(\boldsymbol{X}_a) Φ ( X a ) Φ ( X b ) \boldsymbol{\Phi}(\boldsymbol{X}_b) Φ ( X b ) (9) . Thus, the function m a b \boldsymbol{m}_{ab} m ab m a b ( X a , X b ) = Φ ( X a ) ⊗ Φ ( X b ) ⋅ γ a b + ξ a b \boldsymbol{m}_{ab}(\boldsymbol{X}_a,\boldsymbol{X}_b)=\boldsymbol{\Phi}(\boldsymbol{X}_a)\otimes\boldsymbol{\Phi}(\boldsymbol{X}_b)\cdot\boldsymbol{\gamma}_{ab}+\boldsymbol{\xi}_{ab} m ab ( X a , X b ) = Φ ( X a ) ⊗ Φ ( X b ) ⋅ γ ab + ξ ab γ a b = ( γ a b 0 T , γ a b 1 T , ⋯ , γ a b , K n − 1 T ) T \boldsymbol{\gamma}_{ab}=(\boldsymbol{\gamma}_{ab0}^{\mathsf{T}},\boldsymbol{\gamma}_{ab1}^{\mathsf{T}},\cdots,\boldsymbol{\gamma}_{ab,K_n-1}^{\mathsf{T}})^{\mathsf{T}} γ ab = ( γ ab 0 T , γ ab 1 T , ⋯ , γ ab , K n − 1 T ) T γ a b k = ( γ a b k 0 , γ a b k 1 , ⋯ , γ a b k , K n − 1 ) \boldsymbol{\gamma}_{abk}=(\gamma_{abk0},\gamma_{abk1},\cdots,\gamma_{abk,K_n-1}) γ abk = ( γ abk 0 , γ abk 1 , ⋯ , γ abk , K n − 1 ) ξ a b \boldsymbol{\xi}_{ab} ξ ab
Remark. Researchers may use different orthonormal bases when the Hilbert space is equipped with a density function π \pi π L 2 ( R ; π ) = { f : ∫ R f 2 ( x ) π ( x ) d x < ∞ } L^2(\mathbb{R};\pi)=\{f:\int_{\mathbb{R}}f^2(x)\pi(x)\,\mathrm{d}x<\infty\} L 2 ( R ; π ) = { f : ∫ R f 2 ( x ) π ( x ) d x < ∞ } π ( x ) = exp { − x 2 } \pi(x)=\exp\{-x^2\} π ( x ) = exp { − x 2 } { H i , i ≥ 0 } \{H_i,i\ge0\} { H i , i ≥ 0 } ϕ i \phi_i ϕ i H i H_i H i r ∈ R r\in R r ∈ R b r 0 = 0 b_{r0}=0 b r 0 = 0 ( a , b ) ∈ S (a,b)\in S ( a , b ) ∈ S γ a b 00 = γ a b 01 = γ a b 10 = 0 \gamma_{ab00}=\gamma_{ab01}=\gamma_{ab10}=0 γ ab 00 = γ ab 01 = γ ab 10 = 0
Model Approximation and Relevant Covariates
The true model can be approximated as
Y = θ + ∑ l ∈ L α l X l + ∑ r ∈ R Φ ( X r ) b r + ∑ ( s , s ′ ) ∈ S Φ ( X s ) ⊗ Φ ( X s ′ ) γ s s ′ + Ξ + U , (12) Y=\theta+\sum_{l\in L}\alpha_lX_l+\sum_{r\in R}\boldsymbol{\Phi}(\boldsymbol{X}_r)b_r+\sum_{(s,s')\in S}\boldsymbol{\Phi}(\boldsymbol{X}_s)\otimes\boldsymbol{\Phi}(\boldsymbol{X}_{s'})\boldsymbol{\gamma}_{ss'}+\boldsymbol{\Xi}+U, \tag*{(12)} Y = θ + l ∈ L ∑ α l X l + r ∈ R ∑ Φ ( X r ) b r + ( s , s ′ ) ∈ S ∑ Φ ( X s ) ⊗ Φ ( X s ′ ) γ s s ′ + Ξ + U , (12) where Ξ = ∑ r ξ r + ∑ a , b ξ a b \boldsymbol{\Xi}=\sum_r\boldsymbol{\xi}_r+\sum_{a,b}\boldsymbol{\xi}_{ab} Ξ = ∑ r ξ r + ∑ a , b ξ ab K n K_n K n U U U n × 1 n\times1 n × 1 ϵ i \epsilon_i ϵ i α L , b R , γ S \boldsymbol{\alpha}_L,\boldsymbol{b}_R,\boldsymbol{\gamma}_S α L , b R , γ S
α L = ( α l 1 , ⋯ , α l L ) T for l i ∈ L , b R = ( b r 1 T , ⋯ , b r R T ) T for r i ∈ R , \boldsymbol{\alpha}_L=(\alpha_{l_1},\cdots,\alpha_{l_L})^{\mathsf{T}}\quad\text{for }l_i\in L,\qquad\boldsymbol{b}_R=(b_{r_1}^{\mathsf{T}},\cdots,b_{r_R}^{\mathsf{T}})^{\mathsf{T}}\quad\text{for }r_i\in R, α L = ( α l 1 , ⋯ , α l L ) T for l i ∈ L , b R = ( b r 1 T , ⋯ , b r R T ) T for r i ∈ R , γ S = ( γ s 1 , s 1 ′ T , ⋯ , γ s S , s S ′ T ) T for ( s i , s i ′ ) ∈ S . (13) \boldsymbol{\gamma}_S=(\boldsymbol{\gamma}_{s_1,s'_1}^{\mathsf{T}},\cdots,\boldsymbol{\gamma}_{s_S,s'_S}^{\mathsf{T}})^{\mathsf{T}}\quad\text{for }(s_i,s'_i)\in S. \tag*{(13)} γ S = ( γ s 1 , s 1 ′ T , ⋯ , γ s S , s S ′ T ) T for ( s i , s i ′ ) ∈ S . (13) We combine them into the single non-zero coefficient vector β P 1 = ( α L , b R , γ S ) T \boldsymbol{\beta}_{P_1}=(\boldsymbol{\alpha}_L,\boldsymbol{b}_R,\boldsymbol{\gamma}_S)^{\mathsf{T}} β P 1 = ( α L , b R , γ S ) T P 1 = dim ( β P 1 ) P_1=\operatorname{dim}(\boldsymbol{\beta}_{P_1}) P 1 = dim ( β P 1 ) dim ( ⋅ ) \operatorname{dim}(\cdot) dim ( ⋅ )
To facilitate the selection procedure, we define all relevant covariates as Z P 1 = ( Z L , Z R , Z S ) \boldsymbol{Z}_{P_1}=(\boldsymbol{Z}_L,\boldsymbol{Z}_R,\boldsymbol{Z}_S) Z P 1 = ( Z L , Z R , Z S )
Z L = ( X l 1 , ⋯ , X l L ) for l i ∈ L , \boldsymbol{Z}_L=(X_{l_1},\cdots,X_{l_L})\quad\text{for }l_i\in L, Z L = ( X l 1 , ⋯ , X l L ) for l i ∈ L , Z R = ( Φ ( X r 1 ) , ⋯ , Φ ( X r R ) ) for r i ∈ R , \boldsymbol{Z}_R=(\boldsymbol{\Phi}(\boldsymbol{X}_{r_1}),\cdots,\boldsymbol{\Phi}(\boldsymbol{X}_{r_R}))\quad\text{for }r_i\in R, Z R = ( Φ ( X r 1 ) , ⋯ , Φ ( X r R )) for r i ∈ R , Z S = ( Φ ( X s 1 ) ⊗ Φ ( X s 1 ′ ) , ⋯ , Φ ( X s S ) ⊗ Φ ( X s S ′ ) ) for ( s i , s i ′ ) ∈ S . (14) \boldsymbol{Z}_S=\left(\boldsymbol{\Phi}(\boldsymbol{X}_{s_1})\otimes\boldsymbol{\Phi}(\boldsymbol{X}_{s'_1}),\cdots,\boldsymbol{\Phi}(\boldsymbol{X}_{s_S})\otimes\boldsymbol{\Phi}(\boldsymbol{X}_{s'_S})\right)\quad\text{for }(s_i,s'_i)\in S. \tag*{(14)} Z S = ( Φ ( X s 1 ) ⊗ Φ ( X s 1 ′ ) , ⋯ , Φ ( X s S ) ⊗ Φ ( X s S ′ ) ) for ( s i , s i ′ ) ∈ S . (14) Estimation
Based on the model above, our method aims to select the relevant subset Q Q Q P P P
Our procedure builds on the seminal works of LASSO (Tibshirani 1996), adaptive LASSO (Zou 2006), and Group LASSO (Yuan & Lin 2006), and we show that S-LASSO enjoys the oracle properties (Fan & Li 2001, Wang & Leng 2008), which require: (1) Selection consistency : correctly identifies relevant variables with probability approaching 1; (2) Estimation efficiency : yields asymptotically normal estimators for the selected variables, as if the true model were known a priori .
Step 1: Initial Selection
First, we include all possible linear, nonlinear, and pairwise interaction forms of P P P
min α , b , γ ∥ Y − ∑ l = 0 P α l X l − ∑ r = 1 P Φ ( X r ) b r − ∑ a = 1 P − 1 ∑ b > a P Φ ( X a ) ⊗ Φ ( X b ) γ a b ∥ 2 2 + λ n ( ∑ l = 0 P ∣ α l ∣ + ∑ r = 1 P ∥ b r ∥ 1 + ∑ a = 1 P − 1 ∑ b > a P ∣ γ a b ∣ ) (15) \min_{\alpha,b,\gamma} \left\lVert Y-\sum_{l=0}^{P}\alpha_lX_l-\sum_{r=1}^{P}\Phi(X_r)b_r-\sum_{a=1}^{P-1}\sum_{b>a}^{P}\Phi(X_a)\otimes\Phi(X_b)\gamma_{ab}\right\rVert_2^2+\lambda_n\left(\sum_{l=0}^{P}|\alpha_l|+\sum_{r=1}^{P}\lVert b_r\rVert_1+\sum_{a=1}^{P-1}\sum_{b>a}^{P}|\gamma_{ab}|\right) \tag*{(15)} α , b , γ min Y − l = 0 ∑ P α l X l − r = 1 ∑ P Φ ( X r ) b r − a = 1 ∑ P − 1 b > a ∑ P Φ ( X a ) ⊗ Φ ( X b ) γ ab 2 2 + λ n ( l = 0 ∑ P ∣ α l ∣ + r = 1 ∑ P ∥ b r ∥ 1 + a = 1 ∑ P − 1 b > a ∑ P ∣ γ ab ∣ ) (15) where the intercept is included in α 0 \alpha_0 α 0 X 0 = 1 X_0=1 X 0 = 1 ∣ ⋅ ∣ |\cdot| ∣ ⋅ ∣ l 1 l_1 l 1 ∥ ⋅ ∥ 2 \lVert\cdot\rVert_2 ∥ ⋅ ∥ 2 l 2 l_2 l 2 λ n > 0 \lambda_n>0 λ n > 0
Remark. If prior information or conventional specifications are available, we may instead define prior candidate sets L 0 L_0 L 0 R 0 R_0 R 0 S 0 S_0 S 0 (15) becomes
min α , b , γ ∥ Y − ∑ l ∈ L 0 α l X l − ∑ r ∈ R 0 Φ ( X r ) b r − ∑ ( a , b ) ∈ S 0 Φ ( X a ) ⊗ Φ ( X b ) γ a b ∥ 2 2 + λ n ( ∑ l ∈ L 0 ∣ α l ∣ + ∑ r ∈ R 0 ∥ b r ∥ 1 + ∑ ( a , b ) ∈ S 0 ∣ γ a b ∣ ) (16) \min_{\alpha,b,\gamma} \left\lVert Y-\sum_{l\in L_0}\alpha_lX_l-\sum_{r\in R_0}\Phi(X_r)b_r-\sum_{(a,b)\in S_0}\Phi(X_a)\otimes\Phi(X_b)\gamma_{ab}\right\rVert_2^2+\lambda_n\left(\sum_{l\in L_0}|\alpha_l|+\sum_{r\in R_0}\lVert b_r\rVert_1+\sum_{(a,b)\in S_0}|\gamma_{ab}|\right) \tag*{(16)} α , b , γ min Y − l ∈ L 0 ∑ α l X l − r ∈ R 0 ∑ Φ ( X r ) b r − ( a , b ) ∈ S 0 ∑ Φ ( X a ) ⊗ Φ ( X b ) γ ab 2 2 + λ n l ∈ L 0 ∑ ∣ α l ∣ + r ∈ R 0 ∑ ∥ b r ∥ 1 + ( a , b ) ∈ S 0 ∑ ∣ γ ab ∣ (16) Step 2: Group Selection
Following Zou (2006) and Wang & Leng (2008), we next apply adaptive group LASSO using the estimates from Step 1. Specifically, we construct penalty weights and then substitute all basis functions into the adaptive group LASSO criterion:
min α , b , γ ∥ Y − ∑ l = 0 P α l X l − ∑ r = 1 P Φ ( X r ) b r − ∑ a = 1 P − 1 ∑ b > a P Φ ( X a ) ⊗ Φ ( X b ) γ a b ∥ 2 2 + λ ~ n ( ∑ l = 0 P w ^ α l ∣ α l ∣ + ∑ r = 1 P w ^ b r ∥ b r ∥ 2 + ∑ a = 1 P − 1 ∑ b > a P w ^ γ a b ∥ γ a b ∥ 2 ) (17) \min_{\alpha,b,\gamma} \left\|Y-\sum_{l=0}^{P}\alpha_l X_l-\sum_{r=1}^{P}\Phi(X_r)b_r-\sum_{a=1}^{P-1}\sum_{b>a}^{P}\Phi(X_a)\otimes\Phi(X_b)\gamma_{ab}\right\|_2^2+\widetilde{\lambda}_n\left(\sum_{l=0}^{P}\widehat{w}_{\alpha_l}|\alpha_l|+\sum_{r=1}^{P}\widehat{w}_{b_r}\|b_r\|_2+\sum_{a=1}^{P-1}\sum_{b>a}^{P}\widehat{w}_{\gamma_{ab}}\|\gamma_{ab}\|_2\right) \tag*{(17)} α , b , γ min Y − l = 0 ∑ P α l X l − r = 1 ∑ P Φ ( X r ) b r − a = 1 ∑ P − 1 b > a ∑ P Φ ( X a ) ⊗ Φ ( X b ) γ ab 2 2 + λ n ( l = 0 ∑ P w α l ∣ α l ∣ + r = 1 ∑ P w b r ∥ b r ∥ 2 + a = 1 ∑ P − 1 b > a ∑ P w γ ab ∥ γ ab ∥ 2 ) (17) where λ ~ n \widetilde{\lambda}_n λ n
w ^ α l = { ∣ α ~ l ∣ − 1 , if ∣ α ~ l ∣ > 0 , ∞ , if ∣ α ~ l ∣ = 0 , w ^ b r = { N R ∥ b ~ r ∥ 2 − 1 , if ∥ b ~ r ∥ 2 > 0 , ∞ , if ∥ b ~ r ∥ 2 = 0 , \widehat{w}_{\alpha_l}=
\begin{cases}
|\widetilde{\alpha}_l|^{-1}, & \text{if }|\widetilde{\alpha}_l|>0,\\
\infty, & \text{if }|\widetilde{\alpha}_l|=0,
\end{cases}
\qquad
\widehat{w}_{b_r}=
\begin{cases}
\sqrt{N_R}\|\widetilde{b}_r\|_2^{-1}, & \text{if }\|\widetilde{b}_r\|_2>0,\\
\infty, & \text{if }\|\widetilde{b}_r\|_2=0,
\end{cases} w α l = { ∣ α l ∣ − 1 , ∞ , if ∣ α l ∣ > 0 , if ∣ α l ∣ = 0 , w b r = { N R ∥ b r ∥ 2 − 1 , ∞ , if ∥ b r ∥ 2 > 0 , if ∥ b r ∥ 2 = 0 , w ^ γ a b = { N S ∥ γ ~ a b ∥ 2 − 1 , if ∥ γ ~ a b ∥ 2 > 0 , ∞ , if ∥ γ ~ a b ∥ 2 = 0. (18) \widehat{w}_{\gamma_{ab}}=
\begin{cases}
\sqrt{N_S}\|\widetilde{\gamma}_{ab}\|_2^{-1}, & \text{if }\|\widetilde{\gamma}_{ab}\|_2>0,\\
\infty, & \text{if }\|\widetilde{\gamma}_{ab}\|_2=0.
\end{cases} \tag*{(18)} w γ ab = { N S ∥ γ ab ∥ 2 − 1 , ∞ , if ∥ γ ab ∥ 2 > 0 , if ∥ γ ab ∥ 2 = 0. (18) Here, α ~ \widetilde{\alpha} α b ~ \widetilde{b} b γ ~ \widetilde{\gamma} γ N R = K n N_R=K_n N R = K n N S = K n 2 N_S=K_n^2 N S = K n 2
Apart from the relevant variables defined in Equation (14) , we define the irrelevant variables as Z P 2 = ( Z L C , Z R C , Z S C ) Z_{P_2}=(Z_{L^C},Z_{R^C},Z_{S^C}) Z P 2 = ( Z L C , Z R C , Z S C )
Z L C = ( X l 1 , … , X l L C ) for l i ∈ L C , Z_{L^C}=\left(X_{l_1},\ldots,X_{l_{L^C}}\right)\quad\text{for }l_i\in L^C, Z L C = ( X l 1 , … , X l L C ) for l i ∈ L C , Z R C = ( Φ ( X r 1 ) , … , Φ ( X r R C ) ) for r i ∈ R C , Z_{R^C}=\left(\Phi(X_{r_1}),\ldots,\Phi(X_{r_{R^C}})\right)\quad\text{for }r_i\in R^C, Z R C = ( Φ ( X r 1 ) , … , Φ ( X r R C ) ) for r i ∈ R C , Z S C = ( Φ ( X s 1 ) ⊗ Φ ( X s 1 ′ ) , … , Φ ( X s S C ) ⊗ Φ ( X s S C ′ ) ) for ( s i , s i ′ ) ∈ S C (19) Z_{S^C}=\left(\Phi(X_{s_1})\otimes\Phi(X_{s'_1}),\ldots,\Phi(X_{s_{S^C}})\otimes\Phi(X_{s'_{S^C}})\right)\quad\text{for }(s_i,s'_i)\in S^C \tag*{(19)} Z S C = ( Φ ( X s 1 ) ⊗ Φ ( X s 1 ′ ) , … , Φ ( X s S C ) ⊗ Φ ( X s S C ′ ) ) for ( s i , s i ′ ) ∈ S C (19) where L C : = ∣ L C ∣ L^C:=|L^C| L C := ∣ L C ∣ R C : = ∣ R C ∣ R^C:=|R^C| R C := ∣ R C ∣ S C : = ∣ S C ∣ S^C:=|S^C| S C := ∣ S C ∣ n × P 1 n\times P_1 n × P 1 Z P 1 Z_{P_1} Z P 1 β P 1 \beta_{P_1} β P 1 n × P 2 n\times P_2 n × P 2 Z P 2 Z_{P_2} Z P 2 β P 2 \beta_{P_2} β P 2 Z = ( Z P 1 , Z P 2 ) Z=(Z_{P_1},Z_{P_2}) Z = ( Z P 1 , Z P 2 ) β P Z = ( β P 1 T , β P 2 T ) T \beta_{P_Z}=(\beta_{P_1}^{\mathsf{T}},\beta_{P_2}^{\mathsf{T}})^{\mathsf{T}} β P Z = ( β P 1 T , β P 2 T ) T E ( Y ∣ P ) = Z β P Z + Ξ E(Y\mid P)=Z\beta_{P_Z}+\Xi E ( Y ∣ P ) = Z β P Z + Ξ Ξ \Xi Ξ K n K_n K n
The resulting model specification is
h ( Q ^ ) ≈ ∑ l ∈ L ^ β ^ l X l + ∑ r ∈ R ^ m r ( X r ) + ∑ ( s , s ′ ) ∈ S ^ m s s ′ ( X s , X s ′ ) (20) h(\widehat{Q}) \approx\sum_{l\in\widehat{L}} \widehat{\beta}_l X_l + \sum_{r\in\widehat{R}} m_r(X_r) + \sum_{(s,s')\in\widehat{S}} m_{ss'}(X_s,X_{s'}) \tag*{(20)} h ( Q ) ≈ l ∈ L ∑ β l X l + r ∈ R ∑ m r ( X r ) + ( s , s ′ ) ∈ S ∑ m s s ′ ( X s , X s ′ ) (20) where Q ^ \widehat{Q} Q L ^ = { l : ∣ α ^ l ∣ > 0 } \widehat{L}=\{l:|\widehat{\alpha}_l|>0\} L = { l : ∣ α l ∣ > 0 } R ^ = { r : ∥ b ^ r ∥ 2 > 0 } \widehat{R}=\{r:\|\widehat{\mathbf{b}}_r\|_2>0\} R = { r : ∥ b r ∥ 2 > 0 } S ^ = { ( s , s ′ ) : ∥ γ ^ s s ′ ∥ 2 > 0 } \widehat{S}=\{(s,s'):\|\widehat{\boldsymbol{\gamma}}_{ss'}\|_2>0\} S = {( s , s ′ ) : ∥ γ s s ′ ∥ 2 > 0 } α ^ \widehat{\boldsymbol{\alpha}} α b ^ \widehat{\mathbf{b}} b γ ^ \widehat{\boldsymbol{\gamma}} γ L 0 L_0 L 0 R 0 R_0 R 0 S 0 S_0 S 0 L ^ = { l ∈ L 0 : ∣ α ^ l ∣ > 0 } \widehat{L}=\{l\in L_0:|\widehat{\alpha}_l|>0\} L = { l ∈ L 0 : ∣ α l ∣ > 0 } R ^ = { r ∈ R 0 : ∥ b ^ r ∥ 2 > 0 } \widehat{R}=\{r\in R_0:\|\widehat{\mathbf{b}}_r\|_2>0\} R = { r ∈ R 0 : ∥ b r ∥ 2 > 0 } S ^ = { ( s , s ′ ) ∈ S 0 : ∥ γ ^ s s ′ ∥ 2 > 0 } \widehat{S}=\{(s,s')\in S_0:\|\widehat{\boldsymbol{\gamma}}_{ss'}\|_2>0\} S = {( s , s ′ ) ∈ S 0 : ∥ γ s s ′ ∥ 2 > 0 }
Theoretical Results
In this section, we establish the oracle properties of S-LASSO. We begin with the assumptions used in the analysis.
Assumptions
Assumption 1. The samples ( X i T , ϵ i ) (\mathbf{X}_i^{\mathsf{T}},\epsilon_i) ( X i T , ϵ i ) E ( ϵ i ∣ X ) = 0 E(\epsilon_i\mid\mathbf{X})=0 E ( ϵ i ∣ X ) = 0 Var ( ϵ i ∣ X ) = σ 2 \operatorname{Var}(\epsilon_i\mid\mathbf{X})=\sigma^2 Var ( ϵ i ∣ X ) = σ 2 E ( ∥ X i ∥ 2 ) < ∞ E(\|\mathbf{X}_i\|^2)<\infty E ( ∥ X i ∥ 2 ) < ∞ ϵ i \epsilon_i ϵ i z 1 , z 2 > 0 z_1,z_2>0 z 1 , z 2 > 0 P ( ∣ ϵ i ∣ > x ) ≤ exp { − ( x / z 1 ) z 2 } P(|\epsilon_i|>x)\le\exp\{-(x/z_1)^{z_2}\} P ( ∣ ϵ i ∣ > x ) ≤ exp { − ( x / z 1 ) z 2 }
Remark. These conditions are standard moment and tail assumptions in high-dimensional models. They are used to control the stochastic score terms that arise in the selection and estimation analysis.
For any K K K
V ( n ) = 1 n ( Z P 1 , Z P 2 ) T ( Z P 1 , Z P 2 ) = [ V 11 ( n ) V 12 ( n ) V 21 ( n ) V 22 ( n ) ] . (21) \mathbf{V}^{(n)}=\frac{1}{n}(\mathbf{Z}_{P_1},\mathbf{Z}_{P_2})^{\mathsf{T}}(\mathbf{Z}_{P_1},\mathbf{Z}_{P_2})=\begin{bmatrix}V_{11}^{(n)} & V_{12}^{(n)} \\ V_{21}^{(n)} & V_{22}^{(n)}\end{bmatrix}. \tag*{(21)} V ( n ) = n 1 ( Z P 1 , Z P 2 ) T ( Z P 1 , Z P 2 ) = [ V 11 ( n ) V 21 ( n ) V 12 ( n ) V 22 ( n ) ] . (21) Since the Hermite functions are uniformly bounded: sup x ∈ R sup k ≥ 0 ∣ ϕ k ( x ) ∣ ≤ π − 1 / 4 \sup_{x\in\mathbb{R}}\sup_{k\ge0}|\phi_k(x)|\le\pi^{-1/4} sup x ∈ R sup k ≥ 0 ∣ ϕ k ( x ) ∣ ≤ π − 1/4 ϕ k ( X r i ) \phi_k(X_{ri}) ϕ k ( X r i ) ϕ k 1 ( X a i ) ϕ k 2 ( X b i ) \phi_{k_1}(X_{ai})\phi_{k_2}(X_{bi}) ϕ k 1 ( X ai ) ϕ k 2 ( X bi ) V ( n ) \mathbf{V}^{(n)} V ( n ) Z = ( Z P 1 , Z P 2 ) = ( Z 1 , … , Z P 1 , Z P 1 + 1 , … , Z P 1 + P 2 ) = ( z 1 T , … , z n T ) T \mathbf{Z}=(\mathbf{Z}_{P_1},\mathbf{Z}_{P_2})=(Z_1,\ldots,Z_{P_1},Z_{P_1+1},\ldots,Z_{P_1+P_2})=(\mathbf{z}_1^{\mathsf{T}},\ldots,\mathbf{z}_n^{\mathsf{T}})^{\mathsf{T}} Z = ( Z P 1 , Z P 2 ) = ( Z 1 , … , Z P 1 , Z P 1 + 1 , … , Z P 1 + P 2 ) = ( z 1 T , … , z n T ) T 1 n z j T z j ≤ M 1 \frac{1}{n}\mathbf{z}_j^{\mathsf{T}}\mathbf{z}_j\le M_1 n 1 z j T z j ≤ M 1 j = 1 , 2 , … , P 1 + P 2 j=1,2,\ldots,P_1+P_2 j = 1 , 2 , … , P 1 + P 2 M 1 M_1 M 1 [3 ] .
Remark. Assumption 1 can be modified for other settings, for example, by replacing independence with suitable mixing conditions in time series applications, provided that it still implies 1 n Z j T Z j ≤ M 1 \frac{1}{n}Z_j^{\mathsf{T}}Z_j \le M_1 n 1 Z j T Z j ≤ M 1
Assumption 2. There exist constants c 1 c_1 c 1 c 2 c_2 c 2 c 3 c_3 c 3 c 4 c_4 c 4 0 ≤ c 1 < c 2 < 1 0 \le c_1 < c_2 < 1 0 ≤ c 1 < c 2 < 1 c 3 , c 4 > 0 c_3,c_4 > 0 c 3 , c 4 > 0 K n ≤ O ( n c 1 2 ) K_n \le O\left(n^{\frac{c_1}{2}}\right) K n ≤ O ( n 2 c 1 ) P 1 × 1 P_1 \times1 P 1 × 1 β P 1 \boldsymbol{\beta}_{P_1} β P 1 d ( β P 1 ) d(\boldsymbol{\beta}_{P_1}) d ( β P 1 ) ∣ β P 1 ∣ \lvert\boldsymbol{\beta}_{P_1}\rvert ∣ β P 1 ∣ n 1 − c 2 2 d ( β P 1 ) ≥ c 4 n^{\frac{1-c_2}{2}}d(\boldsymbol{\beta}_{P_1}) \ge c_4 n 2 1 − c 2 d ( β P 1 ) ≥ c 4 ρ min ( V 11 ) > c 3 \rho_{\min}(V_{11}) > c_3 ρ m i n ( V 11 ) > c 3 n n n ρ min ( ⋅ ) \rho_{\min}(\cdot) ρ m i n ( ⋅ )
Remark. Condition (a) controls the truncation level, condition (b) imposes a minimum signal strength, and condition (c) excludes severe multicollinearity by requiring the smallest eigenvalue of V 11 ( n ) V_{11}^{(n)} V 11 ( n )
To clarify condition (a), Zhao & Yu [3 ] assume P 1 = O ( n c 1 ) P_1 = O\left(n^{c_1}\right) P 1 = O ( n c 1 ) P 2 = o ( n ( c 2 − c 1 ) k ) P_2 = o\left(n^{(c_2-c_1)k}\right) P 2 = o ( n ( c 2 − c 1 ) k ) k > 0 k > 0 k > 0 E ( ϵ i 2 k ) < ∞ E\left(\epsilon_i^{2k}\right) < \infty E ( ϵ i 2 k ) < ∞ P 1 = O ( K n 2 ) P_1 = O\left(K_n^2\right) P 1 = O ( K n 2 ) P 2 = O ( P 2 ⋅ K n 2 ) P_2 = O\left(P^2 \cdot K_n^2\right) P 2 = O ( P 2 ⋅ K n 2 ) [3 ] ; see Lemma A.4 in the Appendix.
Assumption 3. The sample covariance matrix V ( n ) V^{(n)} V ( n ) [3 ] : there exists a positive P 2 × 1 P_2 \times1 P 2 × 1 η \boldsymbol{\eta} η
∣ V 21 ( n ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ∣ ≤ 1 − η , (22) \left|V_{21}^{(n)}\left(V_{11}^{(n)}\right)^{-1}\operatorname{sign}(\boldsymbol{\beta}_{P_1})\right| \le1-\boldsymbol{\eta}, \tag*{(22)} V 21 ( n ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ≤ 1 − η , (22) where the inequality holds element-wise.
Remark. Assumption 3 is needed for the sign consistency in Theorem 4.1. The weak irrepresentable condition, ∣ V 21 ( n ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ∣ ≤ 1 \left|V_{21}^{(n)}\left(V_{11}^{(n)}\right)^{-1}\operatorname{sign}(\boldsymbol{\beta}_{P_1})\right| \le1 V 21 ( n ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ≤ 1
Assumption 4. The functions m r ( x ) m_r(x) m r ( x ) m a b ( x , y ) m_{ab}(x,y) m ab ( x , y ) r ∈ R r \in\mathbb{R} r ∈ R m r ( x ) m_r(x) m r ( x ) l 1 l_1 l 1 l 1 ≥ 1 l_1 \ge1 l 1 ≥ 1 i = 0 , 1 , … , l 1 i = 0,1,\ldots,l_1 i = 0 , 1 , … , l 1 m r ( l 1 − i ) ( x ) x i ∈ L 2 ( Ω ) m_r^{(l_1-i)}(x)x^i \in L^2(\Omega) m r ( l 1 − i ) ( x ) x i ∈ L 2 ( Ω ) F l 1 ( Ω ) \mathcal{F}^{l_1}(\Omega) F l 1 ( Ω ) ( a , b ) ∈ S (a,b) \in\mathcal{S} ( a , b ) ∈ S f a b x 0 ( y ) = m a b ( x 0 , y ) f_{abx_0}(y)=m_{ab}(x_0,y) f ab x 0 ( y ) = m ab ( x 0 , y ) f a b y 0 ( x ) = m a b ( x , y 0 ) f_{aby_0}(x)=m_{ab}(x,y_0) f ab y 0 ( x ) = m ab ( x , y 0 ) F l 2 ( Ω ) \mathcal{F}^{l_2}(\Omega) F l 2 ( Ω ) x 0 , y 0 x_0,y_0 x 0 , y 0 l 2 ≥ 3 l_2 \ge3 l 2 ≥ 3
Remark. These smoothness conditions are in line with Dong et al. [1 ] and ensure the convergence of the orthogonal series approximation.
Selection Consistency
The following definition formalizes the sign consistency. For brevity, we write β ^ = s β \hat{\beta} =_{s} \beta β ^ = s β sign ( β ^ ) = sign ( β ) \operatorname{sign}(\hat{\beta}) = \operatorname{sign}(\beta) sign ( β ^ ) = sign ( β ) sign ( ⋅ ) \operatorname{sign}(\cdot) sign ( ⋅ )
Definition 4.1. Let β ^ \hat{\beta} β ^ β \beta β β ^ \hat{\beta} β ^ Sign Consistent with β \beta β λ n \lambda_n λ n λ ~ n \tilde{\lambda}_n λ ~ n lim n → ∞ P ( β ^ ( λ n ) = s β ) = 1 \lim_{n \to\infty} P\left(\hat{\beta}(\lambda_n) =_{s} \beta\right) = 1 lim n → ∞ P ( β ^ ( λ n ) = s β ) = 1
Recall P 1 = dim ( β P 1 ) P_1 = \dim(\beta_{P_1}) P 1 = dim ( β P 1 ) P 2 = dim ( β P 2 ) P_2 = \dim(\beta_{P_2}) P 2 = dim ( β P 2 ) P Z = P 1 + P 2 P_Z = P_1 + P_2 P Z = P 1 + P 2 β P Z \beta_{P_Z} β P Z
Theorem 4.1. [Selection Consistency] Under Assumptions 1 1 1 4 4 4 λ n \lambda_n λ n λ n n = o ( n c 2 − c 1 2 ) \frac{\lambda_n}{\sqrt{n}} = o\left(n^{\frac{c_2-c_1}{2}}\right) n λ n = o ( n 2 c 2 − c 1 ) λ n n 1 2 + δ → ∞ \frac{\lambda_n}{n^{\frac{1}{2}+\delta}} \to\infty n 2 1 + δ λ n → ∞ δ > 0 \delta> 0 δ > 0 K n → ∞ K_n \to\infty K n → ∞ n → ∞ n \to\infty n → ∞ lim n → ∞ P ( β ^ P Z ≠ s β P Z ) = 0 \lim_{n \to\infty} P\left(\hat{\beta}_{P_Z} \ne_{s} \beta_{P_Z}\right) = 0 lim n → ∞ P ( β ^ P Z = s β P Z ) = 0
Remark. Theorem 4.1 shows that Step 1 correctly identifies the relevant components and shrinks the irrelevant ones to zero. Because the adaptive weights in Equation (17) assign an infinite penalty to zero first-step estimates, the same zero pattern is preserved in Step 2. Hence, lim n → ∞ P ( β ^ P Z ≠ s β P Z ) = 0 \lim_{n \to\infty} P\left(\hat{\beta}_{P_Z} \ne_{s} \beta_{P_Z}\right) = 0 lim n → ∞ P ( β ^ P Z = s β P Z ) = 0
Remark. In classical LASSO theory [3 ] , the irrepresentable condition is almost necessary and sufficient for exact sign consistency. However, it can be restrictive because all entries in Equation (22) are required to hold simultaneously. This is particularly demanding in asset pricing applications, where firm characteristics are often correlated (Freyberger et al. 2020). When the condition fails, the global sign consistency event { β ^ P Z = s β P Z } \{\hat{\beta}_{P_Z} =_{s} \beta_{P_Z}\} { β ^ P Z = s β P Z } § 9.4 .§ 3 reports an additional diagnostic analysis under empirically calibrated dependence structures.
Asymptotic Normality
We now study the consistency and asymptotic normality of the estimator in Step 2.
Theorem 4.2. [Consistency] Under Assumptions 1 1 1 4 4 4 λ ~ n = o ( n c 2 / 2 ) \tilde{\lambda}_n = o(n^{c_2/2}) λ ~ n = o ( n c 2 /2 ) K n = O ( n ν ) K_n = O(n^\nu) K n = O ( n ν ) 0 < ν < 1 / 4 0 < \nu< 1/4 0 < ν < 1/4 α ^ L \hat{\alpha}_L α ^ L α ^ L → P α L \hat{\alpha}_L \xrightarrow{P} \alpha_L α ^ L P α L m ^ r \hat{m}_r m ^ r L 2 L^2 L 2
consistent, i.e., ∥ m r − m ^ r ∥ L 2 = o P ( 1 ) \lVert m_r-\widehat{m}_r\rVert_{L_2}=o_P(1) ∥ m r − m r ∥ L 2 = o P ( 1 ) r ∈ R r\in R r ∈ R m ^ a b \widehat{m}_{ab} m ab L 2 L_2 L 2 ∥ m a b − m ^ a b ∥ L 2 = o P ( 1 ) \lVert m_{ab}-\widehat{m}_{ab}\rVert_{L_2}=o_P(1) ∥ m ab − m ab ∥ L 2 = o P ( 1 ) ( a , b ) ∈ S (a,b)\in S ( a , b ) ∈ S ∥ m r − m ^ r ∥ L 2 ≤ O P ( n − min { 1 2 − ν , 1 − ν 2 } ) \lVert m_r-\widehat{m}_r\rVert_{L_2}\le O_P\left(n^{-\min\left\{\frac{1}{2}-\nu,\frac{1-\nu}{2}\right\}}\right) ∥ m r − m r ∥ L 2 ≤ O P ( n − m i n { 2 1 − ν , 2 1 − ν } ) ∥ m a b − m ^ a b ∥ L 2 ≤ O P ( n − min { 1 2 − ν , ( l 2 − 1 ) ν 2 } ) \lVert m_{ab}-\widehat{m}_{ab}\rVert_{L_2}\le O_P\left(n^{-\min\left\{\frac{1}{2}-\nu,\frac{(l_2-1)\nu}{2}\right\}}\right) ∥ m ab − m ab ∥ L 2 ≤ O P ( n − m i n { 2 1 − ν , 2 ( l 2 − 1 ) ν } )
Remark. This theorem shows that, once irrelevant variables are removed, Step 2 yields consistent estimators. The L 2 L_2 L 2
We next study the asymptotic distribution of the estimators. Define
T ( x ) = ( α ^ L − α L , ( m ^ r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m ^ a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S ) (23) \mathcal{T}(x)=\begin{pmatrix}
\widehat{\alpha}_L-\alpha_L,\\
\left(\widehat{m}_{r_j}(x_{r_j})-m_{r_j}(x_{r_j})\right)_{1\le j\le R}\\
\left(\widehat{m}_{a_jb_j}(x_{a_j},x_{b_j})-m_{a_jb_j}(x_{a_j},x_{b_j})\right)_{1\le j\le S}
\end{pmatrix} \tag*{(23)} T ( x ) = α L − α L , ( m r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S (23) which collects all parameters and functions at a given x ∈ Ω Q x\in\Omega^Q x ∈ Ω Q L + R + S ≤ 2 Q + Q ( Q − 1 ) 2 L+R+S\le2Q+\frac{Q(Q-1)}{2} L + R + S ≤ 2 Q + 2 Q ( Q − 1 ) T ( x ) \mathcal{T}(x) T ( x )
Theorem 4.3. [n \sqrt{n} n Theorem 4.2 and assume K n = O ( n ν ) K_n=O(n^\nu) K n = O ( n ν ) max { 6 6 l 1 − 5 , 6 6 l 2 − 17 } ≤ ν < 1 / 4 \max\left\{\frac{6}{6l_1-5},\frac{6}{6l_2-17}\right\}\le\nu<1/4 max { 6 l 1 − 5 6 , 6 l 2 − 17 6 } ≤ ν < 1/4
n ( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1 / 2 T ( x ) → d N ( 0 , I L + R + S ) , (24) \sqrt{n}\left(\sigma^2G(x)V_{11}^{-1}G(x)^\mathsf{T}\right)^{-1/2}\mathcal{T}(x)\xrightarrow{d}N(0,I_{L+R+S}), \tag*{(24)} n ( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1/2 T ( x ) d N ( 0 , I L + R + S ) , (24) where G ( x ) = diag { I L , diag { Φ ( x r j ) } 1 ≤ j ≤ R , diag { Φ ( x a j ) ⊗ Φ ( x b j ) } 1 ≤ j ≤ S } G(x)=\operatorname{diag}\left\{I_L,\operatorname{diag}\left\{\Phi(x_{r_j})\right\}_{1\le j\le R},\operatorname{diag}\left\{\Phi(x_{a_j})\otimes\Phi(x_{b_j})\right\}_{1\le j\le S}\right\} G ( x ) = diag { I L , diag { Φ ( x r j ) } 1 ≤ j ≤ R , diag { Φ ( x a j ) ⊗ Φ ( x b j ) } 1 ≤ j ≤ S }
Remark. Theorem 4.3 implies that the linear coefficients in Step 2 are n \sqrt{n} n x ∈ Ω Q x\in\Omega^Q x ∈ Ω Q n \sqrt{n} n
Based on Theorem 4.3, for a specific r 0 r_0 r 0 n ( m ^ r 0 ( x ) − m r 0 ( x ) ) Φ ( x ) T Σ b r 0 − 1 Φ ( x ) → d N ( 0 , σ 2 ) \frac{\sqrt{n}(\widehat{m}_{r_0}(x)-m_{r_0}(x))}{\sqrt{\Phi(x)^\mathsf{T}\Sigma_{b r_0}^{-1}\Phi(x)}}\xrightarrow{d}N(0,\sigma^2) Φ ( x ) T Σ b r 0 − 1 Φ ( x ) n ( m r 0 ( x ) − m r 0 ( x )) d N ( 0 , σ 2 ) Σ b r 0 : = Cov ( Φ ( X r 0 1 ) ∣ z − r 0 , 1 ) \Sigma_{b r_0}:=\operatorname{Cov}(\Phi(X_{r_0 1})\mid z_{-r_0,1}) Σ b r 0 := Cov ( Φ ( X r 0 1 ) ∣ z − r 0 , 1 )
z − r 0 , 1 = ( { X l 1 : l ∈ L } , { Φ ( X r 1 ) : r ≠ r 0 , r ∈ R } , { Φ ( X a 1 ) ⊗ Φ ( X b 1 ) : ( a , b ) ∈ S } ) . (25) z_{-r_0,1}=\left(\{X_{l1}:l\in L\},\{\Phi(X_{r1}):r\ne r_0,r\in R\},\{\Phi(X_{a1})\otimes\Phi(X_{b1}):(a,b)\in S\}\right). \tag*{(25)} z − r 0 , 1 = ( { X l 1 : l ∈ L } , { Φ ( X r 1 ) : r = r 0 , r ∈ R } , { Φ ( X a 1 ) ⊗ Φ ( X b 1 ) : ( a , b ) ∈ S } ) . (25) Similarly, for the interactive function m a 0 , b 0 ( x , y ) m_{a_0,b_0}(x,y) m a 0 , b 0 ( x , y )
n ( m ^ a 0 b 0 ( x , y ) − m a 0 b 0 ( x , y ) ) ( Φ ( x ) ⊗ Φ ( y ) ) T Σ γ a 0 b 0 − 1 ( Φ ( x ) ⊗ Φ ( y ) ) → d N ( 0 , σ 2 ) , (26) \frac{\sqrt{n}\left(\widehat{m}_{a_0b_0}(x,y)-m_{a_0b_0}(x,y)\right)}{\sqrt{(\Phi(x)\otimes\Phi(y))^\mathsf{T}\Sigma_{\gamma_{a_0b_0}}^{-1}(\Phi(x)\otimes\Phi(y))}}\xrightarrow{d}N(0,\sigma^2), \tag*{(26)} ( Φ ( x ) ⊗ Φ ( y ) ) T Σ γ a 0 b 0 − 1 ( Φ ( x ) ⊗ Φ ( y )) n ( m a 0 b 0 ( x , y ) − m a 0 b 0 ( x , y ) ) d N ( 0 , σ 2 ) , (26) where σ 2 = Var ( ϵ i ∣ X ) \sigma^2 = \operatorname{Var}(\epsilon_i \mid\mathbf{X}) σ 2 = Var ( ϵ i ∣ X ) Σ γ a 0 b 0 : = Cov ( Φ ( X a 0 , 1 ) ⊗ Φ ( X b 0 , 1 ) ∣ z − ( a 0 , b 0 ) , 1 ) \boldsymbol{\Sigma}_{\gamma_{a_0b_0}} := \operatorname{Cov}(\Phi(X_{a_0,1}) \otimes\Phi(X_{b_0,1}) \mid z_{-(a_0,b_0),1}) Σ γ a 0 b 0 := Cov ( Φ ( X a 0 , 1 ) ⊗ Φ ( X b 0 , 1 ) ∣ z − ( a 0 , b 0 ) , 1 )
z − ( a 0 , b 0 ) , 1 = ( { X l 1 : l ∈ L } , { Φ ( X r 1 ) : r ∈ R } , { Φ ( X a 1 ) ⊗ Φ ( X b 1 ) : ( a , b ) ≠ ( a 0 , b 0 ) , ( a , b ) ∈ S } ) . (27) z_{-(a_0,b_0),1} = \left(\{X_{l1}: l \in L\}, \{\Phi(X_{r1}): r \in R\}, \{\Phi(X_{a1}) \otimes\Phi(X_{b1}): (a,b) \ne(a_0,b_0), (a,b) \in S\}\right). \tag*{(27)} z − ( a 0 , b 0 ) , 1 = ( { X l 1 : l ∈ L } , { Φ ( X r 1 ) : r ∈ R } , { Φ ( X a 1 ) ⊗ Φ ( X b 1 ) : ( a , b ) = ( a 0 , b 0 ) , ( a , b ) ∈ S } ) . (27) Simulation Study
Data Generating Processes
In this section, we examine the finite-sample performance of S-LASSO using Monte Carlo simulations. We generate the dataset { ( y i , X i T ) : i = 1 , … , n } \{(y_i, X_i^{\mathsf{T}}): i = 1, \ldots, n\} {( y i , X i T ) : i = 1 , … , n }
y i = β x 1 i + m 1 ( x 2 i ) + m 2 ( x 3 i , x 4 i ) + ϵ i , i = 1 , … , n , (28) y_i = \beta x_{1i} + m_1(x_{2i}) + m_2(x_{3i}, x_{4i}) + \epsilon_i,\quad i = 1, \ldots, n, \tag*{(28)} y i = β x 1 i + m 1 ( x 2 i ) + m 2 ( x 3 i , x 4 i ) + ϵ i , i = 1 , … , n , (28) which contains one linear effect, one nonlinear univariate effect, and one interaction effect. The covariates satisfy X i ∼ N ( 0 , Σ ) X_i \sim N(0, \boldsymbol{\Sigma}) X i ∼ N ( 0 , Σ ) N ( 0 , 1 ) N(0,1) N ( 0 , 1 ) ϵ i ∼ N ( 0 , σ 2 ) \epsilon_i \sim N(0,\sigma^2) ϵ i ∼ N ( 0 , σ 2 )
Covariate correlation matrix. The correlation matrix Σ = ( σ i j ) P × P \boldsymbol{\Sigma} = (\sigma_{ij})_{P \times P} Σ = ( σ ij ) P × P σ i j \sigma_{ij} σ ij § 6 . The correlation structure of empirical data is discussed in Appendix § 9.4 .§ 3 .
Two designs for the true functions. We consider two data-generating processes. In DGP 1, we set β = 1 \beta= 1 β = 1 m 1 ( x 2 i ) = 2.24 cos ( x 2 i ) m_1(x_{2i}) = 2.24\cos(x_{2i}) m 1 ( x 2 i ) = 2.24 cos ( x 2 i ) m 2 ( x 3 i , x 4 i ) = 1.43 sin ( x 3 i + x 4 i ) m_2(x_{3i}, x_{4i}) = 1.43\sin(x_{3i} + x_{4i}) m 2 ( x 3 i , x 4 i ) = 1.43 sin ( x 3 i + x 4 i ) β = 1 \beta= 1 β = 1 m 1 ( x 2 i ) = 2.24 cos ( x 2 i ) m_1(x_{2i}) = 2.24\cos(x_{2i}) m 1 ( x 2 i ) = 2.24 cos ( x 2 i ) m 2 ( x 3 i , x 4 i ) = 1.91 sin ( x 3 i x 4 i ) m_2(x_{3i}, x_{4i}) = 1.91\sin(x_{3i}x_{4i}) m 2 ( x 3 i , x 4 i ) = 1.91 sin ( x 3 i x 4 i ) x 3 i x 4 i x_{3i}x_{4i} x 3 i x 4 i
Sample size and error level. We fix the sample size at n = 200 n = 200 n = 200 P = 100 P = 100 P = 100 P = 300 P = 300 P = 300 P < n P < n P < n P > n P > n P > n σ = 1 \sigma= 1 σ = 1
Simulation Results
We compare the proposed S-LASSO with several benchmark models:
LASSO. We apply the standard LASSO to the original linear covariates only.
LN-LASSO. We apply the LASSO to a dictionary containing linear terms and univariate nonlinear basis expansions, but no interaction terms.
FNW. We implement the approach of Freyberger et al. (2020), which includes linear terms, nonlinear univariate effects, and product-type interaction terms taking the form m ( x a x b ) m(x_a x_b) m ( x a x b )
The key benchmark is the FNW method of Freyberger et al. (2020). To provide a balanced comparison, DGP 2 is intentionally specified in a way that is favorable to this method. Since the true functions are smooth, we use Hermite basis expansions for nonlinear terms and tensor-product expansions for interaction terms, with truncation level K = n 1 / 6 K = n^{1/6} K = n 1/6
For all methods and steps, the tuning parameter is selected by the Bayesian Information Criterion (BIC), defined as B I C = log ( M S E ) + d f ⋅ log n n \mathrm{BIC} = \log(\mathrm{MSE}) + df \cdot\frac{\log n}{n} BIC = log ( MSE ) + df ⋅ n l o g n d f df df
Table 1 shows four main patterns. First, S-LASSO recovers the linear and nonlinear effects almost perfectly across both DGPs and both values of P P P P = 300 P = 300 P = 300 P P P
Figure 1 reports the estimated functions under DGP 1 with n = 200 n = 200 n = 200 P = 300 P = 300 P = 300
Empirical Study
In this section, we revisit the question in Freyberger et al. (2020) of how firm characteristics predict stock returns in a nonparametric framework. They model asset returns as additive nonparametric functions of lagged firm characteristics and allow pairwise multiplicative interactions with firm size:
Table 1: Simulation Results. This table compares S-LASSO with three benchmarks: standard linear LASSO, LN-LASSO with linear and univariate nonlinear terms but no interactions, and FNW with univariate nonlinear terms and product-form interactions. ‘LIN.’, ‘NONL.’, and ‘INTE.’ report correct detection rates for the true linear, nonlinear, and interaction effects, respectively. ‘FP’ denotes the false positive rate and ‘MSE’ denotes mean squared error. Entries are averages over 100 replications, with standard deviations in parentheses.
Table 1.
E ( Y ∣ P ) = θ 1 + ∑ r = 1 R m r ( X r ) + ∑ s ≠ s i z e m s , s i z e ( X s ⋅ X s i z e ) , (29) E(Y \mid P) = \theta_{1} + \sum_{r=1}^{R} m_{r}(X_{r}) + \sum_{s \ne\mathrm{size}} m_{s,\mathrm{size}}(X_{s} \cdot X_{\mathrm{size}}), \tag*{(29)} E ( Y ∣ P ) = θ 1 + r = 1 ∑ R m r ( X r ) + s = size ∑ m s , size ( X s ⋅ X size ) , (29) where the unknown interaction functions take inputs of the form X s ⋅ X s i z e = ( X s 1 X s i z e , 1 , ⋯ , X s n X s i z e , n ) T X_{s} \cdot X_{\mathrm{size}} = (X_{s1}X_{\mathrm{size},1}, \cdots, X_{sn}X_{\mathrm{size},n})^{\mathsf{T}} X s ⋅ X size = ( X s 1 X size , 1 , ⋯ , X s n X size , n ) T
We apply S-LASSO to re-examine the effects of firm characteristics on stock returns and their interactions with firm size. Our specification is
E ( Y ∣ P ) = θ 1 + ∑ l = 1 L α l X l + ∑ r = 1 R m r ( X r ) + ∑ s ≠ s i z e m s , s i z e ( X s , X s i z e ) , (30) E(Y \mid P) = \theta_{1} + \sum_{l=1}^{L} \alpha_{l}X_{l} + \sum_{r=1}^{R} m_{r}(X_{r}) + \sum_{s \ne\mathrm{size}} m_{s,\mathrm{size}}(X_{s}, X_{\mathrm{size}}), \tag*{(30)} E ( Y ∣ P ) = θ 1 + l = 1 ∑ L α l X l + r = 1 ∑ R m r ( X r ) + s = size ∑ m s , size ( X s , X size ) , (30) where the notation is consistent with Equation (29) . Relative to Freyberger et al. (2020), our model adds linear terms and allows more flexible interactions with firm size. We also assume that both slope coefficients and characteristic functions are time-invariant within each rolling window.
Figure 1. Estimated and True Functions.
Remark. The interaction terms in S-LASSO can be chosen according to the empirical question. In the main application, we focus on interactions with firm size, motivated by the central role of size in asset pricing [2 ] and by the benchmark of Freyberger et al. (2020). These size-centered interactions improve interpretability and facilitate comparison with the existing nonparametric benchmark. Appendix § 9.4 .§ 2 further reports a full-specification check using all pairwise interactions over consecutive four-year subsamples; firm size is selected in all ten subsamples and is the most frequently connected characteristic, supporting the main empirical specification.
Data
Following Freyberger et al. (2020), our empirical objective is to identify which firm characteristics predict monthly returns and through which functional forms. We use the 102 characteristics listed in Green et al. (2017), all of which can be constructed from CRSP, Compustat, and/or I/B/E/S.
Table 2. Summary of Data. This table summarizes the return data and characteristics data used in this paper. For the return data, the descriptive statistics are calculated across all sample dates and sample firms, where p(J-B) denotes the p p p
Table 2 (PDF p. 20).
We also collect monthly stock returns from CRSP and adjust delisting returns following Hou et al. (2015). The sample spans January 1980 to December 2020. Summary statistics are reported in Table 2. Variable definitions are reported in Appendix § 9.4 .§ 1 . Since only four characteristics have complete coverage (Green et al. 2017), we handle missing data using the method of Freyberger et al. (2025), which applies conditional mean imputation and weighted least squares within a GMM framework.
Specification Results
Selected Characteristics
We apply S-LASSO to Equation (30) to determine which characteristics predict returns and through which functional forms, with tuning parameters selected by BIC. To accommodate potential time-varying parameters and functional forms, we divide the full sample into ten subperiods: 1980–1983, 1984–1987, ⋯ \cdots ⋯
Table 3. Summary of Specified Characteristics by S-LASSO. This table reports how frequently each firm characteristic is selected by S-LASSO across ten four-year subperiods. Panels A and B list characteristics selected frequently as linear and nonlinear effects, respectively. Panel C lists characteristics selected as having interaction effects with firm size, measured by market equity (‘mve’). The tuning parameter is selected by BIC.
Several characteristics, such as ‘cash’ and ‘turn’, are selected repeatedly as linear effects, while others, especially ‘indmom’ and ‘mom1m’, are repeatedly selected as nonlinear effects. Overall, 99 characteristics are selected at least once in the linear component and 87 at least once in the nonlinear component, indicating substantial functional-form heterogeneity. Panel C reports interactions with the market value of equity (‘mve’). In total, 13 characteristics exhibit interaction effects with size, including ‘mom36m’, ‘beta’, ‘mom1m’, and ‘ms’.
These results yield three main findings: (i) Functional-Form Heterogeneity: Firm characteristics predict stock returns through different functional forms. Both linear and nonlinear effects are pervasive, suggesting that restricting all characteristics to enter linearly would miss important nonlinear predictive patterns. (ii) Size-Dependent Effects: We find strong interactions between firm size (‘mve’) and 13 characteristics, including momentum (‘mom36m’ and ‘mom1m’), volatility-related measures (‘std_turn’), and profitability indicators (‘gma’). (iii) Comparison with Existing Studies: Relative to Freyberger et al. (2020), our specification reveals both linear effects and richer size-dependent interactions. In particular, variables such as ‘mom1m’, ‘cash’, and ‘turn’ enter linearly in our models, while characteristics including ‘beta’, ‘gma’, and ‘ms’ show flexible interactions with size beyond the multiplicative form.
Economic Interpretations
To interpret in-sample selection results, we visualize representative estimated functions from the 2016–2019 window.
Representative Nonlinear Effects.
We first examine two selected nonlinear effects: ‘cash’ and ‘egr’. Figure 2 reports their estimated total effects after centering. For ‘cash’, the estimated function displays a nonmonotonic, approximately Z-shaped pattern: firms with very low cash holdings are associated with higher predicted returns, the effect becomes negative for firms with intermediate cash holdings, and it rises again toward the upper tail of the cash distribution. This pattern suggests that cash holdings do not enter the return prediction model through a single monotone channel. One interpretation is that very low cash holdings may proxy for liquidity pressure or financing constraints, whereas very high cash holdings may reflect precautionary savings by firms with greater exposure to aggregate risk, consistent with the mechanism in Palazzo (2012). For ‘egr’, growth in common shareholder equity, the estimated function is hump-shaped: the effect is negative in both tails of the distribution and positive around the middle range. This pattern indicates that common-equity growth does not enter the return prediction model through a simple linear channel. The negative effect for high equity growth is consistent with prior evidence that equity issuance, external financing, and firm expansion are informative for future returns (Pontiff & Woodgate 2008; Cooper et al. 2008).
Figure 2. Estimated nonlinear effects of ‘cash’ and ‘egr’. Each curve shows the estimated total effect after centering. The horizontal axes are standardized.
Momentum Effects and Size Interaction.
In the same subsample, S-LASSO assigns both linear and nonlinear components to ‘mom1m’, ‘mom6m’, ‘mom12m’, and ‘mom36m’, and selects an interaction between ‘mom1m’ and ‘mve’. Figure 3 shows that the selected momentum effects are not linear. The one-month momentum effect is steeper for low recent returns and becomes flatter as recent returns increase. Additionally, the six-month momentum ’mom6m’ is relatively weak. The effect of ‘mom36m’ is negative over most of its range, while its slope also changes with the level of past returns. Thus, momentum is not a homogeneous predictor; its contribution depends on the horizon over which past returns are measured, consistent with evidence that short-, intermediate-, and long-horizon past returns capture different return-predictability effects (Novy-Marx 2012).
Figure 3. Estimated total effects of momentum variables at different horizons. Each curve combines the selected linear and nonlinear components. The horizontal axis is standardized.
Figure 4 plots the estimated interaction surface between one-month momentum and firm size. The surface is non-flat: the interaction effect changes more sharply in the low-‘mve’ region and becomes flatter as ‘mve’ increases. The size-dependent interaction pattern is also broadly consistent with the asset-pricing literature showing that the profitability of momentum-related strategies depends on firm size and related information frictions [2 ] . More broadly, classic studies suggest that return continuation and reversal patterns can vary systematically with firm size. Jegadeesh (1990) documents the existence of short-horizon return reversals and examines serial-correlation patterns across different size groups of stocks. For longer-horizon reversal (‘mom36m’), there is also evidence that the effect is substantially stronger for smaller firms than for larger firms (Chopra et al. (1992)).
Figure 4. Estimated interaction surface between one-month momentum ‘mom1m’ and firm size ‘mve’. Both axes are standardized.
More recent studies likewise show that many anomalies are more pronounced among small firms, including Hou et al. (2020) and Green et al. (2017). Taken together, this evidence motivates our choice to interact characteristics with firm size.
Out-of-sample Performance
We compare S-LASSO with benchmark methods using mean squared prediction error and long-short portfolio performance.
Out-of-sample Mean Squared Error: Using a rolling four-year estimation window, we estimate the model, forecast one-month-ahead returns, and compute the mean squared prediction error (MSE). Repeating this procedure through 2020 allows us to assess predictive stability across different market conditions.
Long–short Portfolios: Following Green et al. (2017), each month we rank stocks by model-predicted returns using only the preceding four years of data, and form an equal-weighted portfolio that goes long the top 20% and short the bottom 20% of
stocks.
Table 4. Out-of-sample Performance. The two portfolio blocks report performance over different evaluation periods, while the MSE block summarizes monthly prediction errors.
Table 4 reports the out-of-sample results. S-LASSO outperforms FNW across the reported prediction-error and portfolio-performance metrics. Its mean squared prediction error falls to 4.03, a 16.7% improvement over FNW’s 4.84, and its forecasts are more stable, with an MSE standard deviation of 16.67 versus 21.72. The corresponding long–short portfolio earns 22.7% per year with a Sharpe ratio of 1.69, compared with 20.8% and 1.56 for FNW, while the maximum drawdown remains nearly identical (40.1% versus 40.3%). Figure 5 plots the log cumulative portfolio returns since 1996 and shows that S-LASSO consistently outperforms the benchmark strategies over the sample period.
Figure 5. Out-of-sample portfolios.
Conclusion
This paper proposes Specification LASSO (S-LASSO), a unified framework for variable selection and model specification in high-dimensional semiparametric models. By combining linear, nonlinear, and pairwise interaction components, S-LASSO selects both relevant covariates and their functional forms. Theoretically, we establish oracle properties, including selection consistency and asymptotic normality. Monte Carlo simulations show that S-LASSO remains robust under correlated covariates and can distinguish linear, nonlinear, and interaction effects in high-dimensional settings. Empirically, we apply S-LASSO to characteristics-based asset pricing. The results show that firm characteristics affect returns through heterogeneous channels, including linear effects, nonlinear effects, and size-dependent interactions. Compared with the benchmark model of Freyberger et al. (2020), S-LASSO improves the out-of-sample predictive performance. Overall, S-LASSO offers an interpretable and practical tool for empirical researchers working with rich covariate spaces. Future work may extend the framework to higher-order interactions and time-varying effects.
Acknowledgments
The authors thank the Editor, the Associate Editor, and the anonymous referee for their constructive comments and suggestions. Chaohua Dong acknowledges the financial support from National Natural Science Foundation of China (project code No. 72473156); Shuyi Ge is supported financially by National Natural Science Foundation of China (project code No. 72673092, 72303119, 72495155); Shaoran Li acknowledges the financial support from National Natural Science Foundation of China (project code No. 72573004, 72203007) and School of Economics, Peking University; Wen Su gratefully acknowledges support from the EPSRC Centre for Doctoral Training in Mathematics of Random Systems (EP/S023925/1). The authors report there are no competing interests to declare.
Supplementary Appendices
Technical Lemmas
This section establishes the fundamental properties of the basis expansions and the resulting design matrix that are used throughout the theoretical analysis. The first lemma records the approximation accuracy of the Hermite expansion, while the second lemma controls the column norms of the expanded design matrix. These two ingredients are used repeatedly below: the approximation bound controls the sieve truncation error, and the column-norm bound ensures that the stochastic score terms in the high-dimensional LASSO analysis are well behaved.
Lemma A.1. Suppose g ( x ) g(x) g ( x ) R \mathbb{R} R m m m m ≥ 1 m \ge1 m ≥ 1 g ( m − l ) ( x ) x l ∈ L 2 ( R ) g^{(m-l)}(x)x^l \in L^2(\mathbb{R}) g ( m − l ) ( x ) x l ∈ L 2 ( R ) l = 0 , 1 , ⋯ , m l = 0,1,\cdots,m l = 0 , 1 , ⋯ , m g ( x ) = ∑ k = 0 ∞ b k ϕ k ( x ) g(x) = \sum_{k=0}^{\infty} b_k\phi_k(x) g ( x ) = ∑ k = 0 ∞ b k ϕ k ( x )
(1) sup x ∣ g ( x ) − ∑ k = 0 K − 1 b k ϕ k ( x ) ∣ = o ( 1 ) K − m − 1 2 − 1 12 \sup_x \left|g(x)-\sum_{k=0}^{K-1}b_k\phi_k(x)\right|=o(1)K^{-\frac{m-1}{2}-\frac{1}{12}} sup x g ( x ) − ∑ k = 0 K − 1 b k ϕ k ( x ) = o ( 1 ) K − 2 m − 1 − 12 1
(2) ∥ g ( x ) − ∑ k = 0 K − 1 b k ϕ k ( x ) ∥ L 2 = o ( 1 ) K − m 2 \left\|g(x)-\sum_{k=0}^{K-1}b_k\phi_k(x)\right\|_{L^2} = o(1)K^{-\frac{m}{2}} g ( x ) − ∑ k = 0 K − 1 b k ϕ k ( x ) L 2 = o ( 1 ) K − 2 m
Proof. See Lemma C.1 in the Supplemental Material of Dong et al. [1 ] .
Lemma A.2. Let Z = ( z i j ) n × ( P 1 + P 2 ) = ( Z 1 , Z 2 , … , Z P 1 , Z P 1 + 1 , … , Z P 1 + P 2 ) Z = (z_{ij})_{n \times(P_1+P_2)} = (Z_1,Z_2,\ldots,Z_{P_1},Z_{P_1+1},\ldots,Z_{P_1+P_2}) Z = ( z ij ) n × ( P 1 + P 2 ) = ( Z 1 , Z 2 , … , Z P 1 , Z P 1 + 1 , … , Z P 1 + P 2 )
1 n X k T X k ≤ M 0 , for any k , (31) \frac{1}{n}X_k^{\mathsf{T}}X_k \leq M_0,\qquad\text{for any } k, \tag*{(31)} n 1 X k T X k ≤ M 0 , for any k , (31) then
1 n Z j T Z j ≤ M 1 , j = 1 , 2 , … , P 1 + P 2 , (32) \frac{1}{n}Z_j^{\mathsf{T}}Z_j \leq M_1,\qquad j=1,2,\ldots,P_1+P_2, \tag*{(32)} n 1 Z j T Z j ≤ M 1 , j = 1 , 2 , … , P 1 + P 2 , (32) where M 1 M_1 M 1
Remark. The condition 1 n X k T X k ≤ M 0 \frac{1}{n}X_k^{\mathsf{T}}X_k \leq M_0 n 1 X k T X k ≤ M 0 Z Z Z
Proof. We consider three cases for z i j z_{ij} z ij Z Z Z
(a) If z i j z_{ij} z ij Z j Z_j Z j
1 n Z j T Z j = 1 n ∑ i = 1 n z i j 2 ≤ M 0 . (33) \frac{1}{n}Z_j^{\mathsf{T}}Z_j = \frac{1}{n}\sum_{i=1}^{n}z_{ij}^{2} \leq M_0. \tag*{(33)} n 1 Z j T Z j = n 1 i = 1 ∑ n z ij 2 ≤ M 0 . (33) (b) If z i j z_{ij} z ij z i j = ϕ k ( X r i ) z_{ij}=\phi_k(X_{ri}) z ij = ϕ k ( X r i ) X r X_r X r k k k ∣ ϕ k ( x ) ∣ ≤ π − 1 / 4 \lvert\phi_k(x)\rvert\leq\pi^{-1/4} ∣ ϕ k ( x )∣ ≤ π − 1/4 z i j 2 ≤ π − 1 / 2 z_{ij}^{2}\leq\pi^{-1/2} z ij 2 ≤ π − 1/2 i i i
1 n Z j T Z j = 1 n ∑ i = 1 n z i j 2 ≤ π − 1 / 2 . (34) \frac{1}{n}Z_j^{\mathsf{T}}Z_j = \frac{1}{n}\sum_{i=1}^{n}z_{ij}^{2} \leq\pi^{-1/2}. \tag*{(34)} n 1 Z j T Z j = n 1 i = 1 ∑ n z ij 2 ≤ π − 1/2 . (34) (c) If z i j z_{ij} z ij z i j = ϕ k 1 ( X a i ) ϕ k 2 ( X b i ) z_{ij}=\phi_{k_1}(X_{ai})\phi_{k_2}(X_{bi}) z ij = ϕ k 1 ( X ai ) ϕ k 2 ( X bi ) ( a , b ) (a,b) ( a , b ) k 1 k_1 k 1 k 2 k_2 k 2 ∣ z i j ∣ ≤ π − 1 / 2 \lvert z_{ij}\rvert\leq\pi^{-1/2} ∣ z ij ∣ ≤ π − 1/2 z i j 2 ≤ π − 1 z_{ij}^{2}\leq\pi^{-1} z ij 2 ≤ π − 1 i i i
1 n Z j T Z j = 1 n ∑ i = 1 n z i j 2 ≤ π − 1 . (35) \frac{1}{n}Z_j^{\mathsf{T}}Z_j = \frac{1}{n}\sum_{i=1}^{n}z_{ij}^{2} \leq\pi^{-1}. \tag*{(35)} n 1 Z j T Z j = n 1 i = 1 ∑ n z ij 2 ≤ π − 1 . (35) Let M 1 = max { π − 1 / 2 , M 0 } M_1=\max\{\pi^{-1/2},M_0\} M 1 = max { π − 1/2 , M 0 } π − 1 ≤ π − 1 / 2 \pi^{-1}\le\pi^{-1/2} π − 1 ≤ π − 1/2 M 1 M_1 M 1
1 n Z j T Z j ≤ M 1 , for all j , (36) \frac{1}{n}Z_j^{\mathsf{T}}Z_j\le M_1,\quad\text{for all }j, \tag*{(36)} n 1 Z j T Z j ≤ M 1 , for all j , (36) which proves the claim.
Lemma A.3. Define W ( n ) = 1 n Z T ( U + Ξ ) \mathbf{W}^{(n)}=\frac{1}{\sqrt{n}}\mathbf{Z}^{\mathsf{T}}(\mathbf{U}+\mathbf{\Xi}) W ( n ) = n 1 Z T ( U + Ξ ) K n → ∞ K_n\to\infty K n → ∞
(a) W ( n ) = W ⋅ { 1 + o P ( 1 ) } \mathbf{W}^{(n)}=\mathbf{W}\cdot\{1+o_P(1)\} W ( n ) = W ⋅ { 1 + o P ( 1 )} W ∼ N ( 0 , σ 2 V ) \mathbf{W}\sim N(0,\sigma^2\mathbf{V}) W ∼ N ( 0 , σ 2 V ) W ( n ) = W ⋅ { 1 + o P ( 1 ) } \mathbf{W}^{(n)}=\mathbf{W}\cdot\{1+o_P(1)\} W ( n ) = W ⋅ { 1 + o P ( 1 )} ϑ ∈ R P Z \boldsymbol{\vartheta}\in\mathbb{R}^{P_Z} ϑ ∈ R P Z ϑ T V ϑ \boldsymbol{\vartheta}^{\mathsf{T}}\mathbf{V}\boldsymbol{\vartheta} ϑ T V ϑ ϑ T W ( n ) = ϑ T W + o P ( 1 ) \boldsymbol{\vartheta}^{\mathsf{T}}\mathbf{W}^{(n)}=\boldsymbol{\vartheta}^{\mathsf{T}}\mathbf{W}+o_P(1) ϑ T W ( n ) = ϑ T W + o P ( 1 )
(b) Let W 1 ( n ) \mathbf{W}_1^{(n)} W 1 ( n ) P 1 P_1 P 1 W ( n ) \mathbf{W}^{(n)} W ( n ) W 1 ( n ) = W 1 ⋅ { 1 + o P ( 1 ) } \mathbf{W}_1^{(n)}=\mathbf{W}_1\cdot\{1+o_P(1)\} W 1 ( n ) = W 1 ⋅ { 1 + o P ( 1 )} ϑ ∈ R P 1 \boldsymbol{\vartheta}\in\mathbb{R}^{P_1} ϑ ∈ R P 1 ϑ T V 11 ϑ \boldsymbol{\vartheta}^{\mathsf{T}}\mathbf{V}_{11}\boldsymbol{\vartheta} ϑ T V 11 ϑ ϑ T W 1 ( n ) = ϑ T W 1 + o P ( 1 ) \boldsymbol{\vartheta}^{\mathsf{T}}\mathbf{W}_1^{(n)}=\boldsymbol{\vartheta}^{\mathsf{T}}\mathbf{W}_1+o_P(1) ϑ T W 1 ( n ) = ϑ T W 1 + o P ( 1 )
Proof. (a) We first separate the stochastic error part from the sieve approximation error. Define
W ~ ( n ) = 1 n Z T U = 1 n ∑ i = 1 n z i ε i . (37) \widetilde{\mathbf{W}}^{(n)}=\frac{1}{\sqrt{n}}\mathbf{Z}^{\mathsf{T}}\mathbf{U}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\mathbf{z}_i\varepsilon_i. \tag*{(37)} W ( n ) = n 1 Z T U = n 1 i = 1 ∑ n z i ε i . (37) For any K n K_n K n Z \mathbf{Z} Z
W ~ ( n ) = W ⋅ { 1 + o P ( 1 ) } , (38) \widetilde{\mathbf{W}}^{(n)}=\mathbf{W}\cdot\{1+o_P(1)\}, \tag*{(38)} W ( n ) = W ⋅ { 1 + o P ( 1 )} , (38) where W ∼ N ( 0 , σ 2 V ) \mathbf{W}\sim N(0,\sigma^2\mathbf{V}) W ∼ N ( 0 , σ 2 V )
It remains to show that replacing U \mathbf{U} U U + Ξ \mathbf{U}+\mathbf{\Xi} U + Ξ W ( n ) \mathbf{W}^{(n)} W ( n ) w j ( n ) w_j^{(n)} w j ( n ) w ~ j ( n ) \widetilde{w}_j^{(n)} w j ( n ) j j j W ( n ) \mathbf{W}^{(n)} W ( n ) W ~ ( n ) \widetilde{\mathbf{W}}^{(n)} W ( n )
w ~ j ( n ) − w j ( n ) = 1 n ∑ i = 1 n z i j ε i − 1 n ∑ i = 1 n z i j ( ε i + Ξ i ) = − 1 n ∑ i = 1 n z i j Ξ i , (39) \widetilde{w}_j^{(n)}-w_j^{(n)}
=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}z_{ij}\varepsilon_i-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}z_{ij}(\varepsilon_i+\Xi_i)
=-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}z_{ij}\Xi_i, \tag*{(39)} w j ( n ) − w j ( n ) = n 1 i = 1 ∑ n z ij ε i − n 1 i = 1 ∑ n z ij ( ε i + Ξ i ) = − n 1 i = 1 ∑ n z ij Ξ i , (39) where Ξ i \Xi_i Ξ i i i i
Ξ i 2 ≤ C Q ( ∑ r ∈ R ξ i r 2 + ∑ ( a , b ) ∈ S ξ i a b 2 ) , i = 1 , 2 , … , n , (40) \Xi_i^2 \le C_Q\left(\sum_{r\in R}\xi_{ir}^2+\sum_{(a,b)\in S}\xi_{iab}^2\right),\quad i=1,2,\ldots,n, \tag*{(40)} Ξ i 2 ≤ C Q r ∈ R ∑ ξ i r 2 + ( a , b ) ∈ S ∑ ξ iab 2 , i = 1 , 2 , … , n , (40) where
ξ i r = m r ( X i r ) − Φ ( X i r ) b r = ∑ k ≥ K b r k ϕ k ( X i r ) , (41) \xi_{ir}=m_r(X_{ir})-\Phi(X_{ir})b_r=\sum_{k\ge K}b_{rk}\phi_k(X_{ir}), \tag*{(41)} ξ i r = m r ( X i r ) − Φ ( X i r ) b r = k ≥ K ∑ b r k ϕ k ( X i r ) , (41) and
ξ i a b = m a b ( X i a , X i b ) − Φ ( X i a ) ⊗ Φ ( X i b ) γ a b . (42) \xi_{iab}=m_{ab}(X_{ia},X_{ib})-\Phi(X_{ia})\otimes\Phi(X_{ib})\gamma_{ab}. \tag*{(42)} ξ iab = m ab ( X ia , X ib ) − Φ ( X ia ) ⊗ Φ ( X ib ) γ ab . (42) By Lemma A.1 , the univariate approximation error satisfies ξ i r = o ( 1 ) K − l 1 − 1 2 − 1 12 \xi_{ir}=o(1)K^{-\frac{l_1-1}{2}-\frac{1}{12}} ξ i r = o ( 1 ) K − 2 l 1 − 1 − 12 1 ξ i r 2 = o ( 1 ) K − l 1 + 5 6 \xi_{ir}^2=o(1)K^{-l_1+\frac{5}{6}} ξ i r 2 = o ( 1 ) K − l 1 + 6 5
m a b ( x , y ) = ∑ k = 0 K − 1 γ a b k ( y ) ϕ k ( x ) + o ( 1 ) K − l 2 − 1 2 − 1 12 = ∑ k = 0 K − 1 ( ∑ d = 0 K − 1 γ a b k d ϕ d ( y ) + o ( 1 ) K − l 2 − 1 2 − 1 12 ) ϕ k ( x ) + o ( 1 ) K − l 2 − 1 2 − 1 12 = Φ ( x ) ⊗ Φ ( y ) γ a b + o ( 1 ) K − l 2 − 1 2 + 11 12 . (43) \begin{aligned}
m_{ab}(x,y)&=\sum_{k=0}^{K-1}\gamma_{abk}(y)\phi_k(x)+o(1)K^{-\frac{l_2-1}{2}-\frac{1}{12}}\\
&=\sum_{k=0}^{K-1}\left(\sum_{d=0}^{K-1}\gamma_{abkd}\phi_d(y)+o(1)K^{-\frac{l_2-1}{2}-\frac{1}{12}}\right)\phi_k(x)+o(1)K^{-\frac{l_2-1}{2}-\frac{1}{12}}\\
&=\Phi(x)\otimes\Phi(y)\gamma_{ab}+o(1)K^{-\frac{l_2-1}{2}+\frac{11}{12}}. \tag*{(43)}
\end{aligned} m ab ( x , y ) = k = 0 ∑ K − 1 γ abk ( y ) ϕ k ( x ) + o ( 1 ) K − 2 l 2 − 1 − 12 1 = k = 0 ∑ K − 1 ( d = 0 ∑ K − 1 γ abk d ϕ d ( y ) + o ( 1 ) K − 2 l 2 − 1 − 12 1 ) ϕ k ( x ) + o ( 1 ) K − 2 l 2 − 1 − 12 1 = Φ ( x ) ⊗ Φ ( y ) γ ab + o ( 1 ) K − 2 l 2 − 1 + 12 11 . (43) Therefore, ξ i a b 2 = o ( 1 ) K − l 2 + 17 6 \xi_{iab}^2=o(1)K^{-l_2+\frac{17}{6}} ξ iab 2 = o ( 1 ) K − l 2 + 6 17
Since Q < ∞ Q<\infty Q < ∞
Ξ i 2 ≤ o ( 1 ) ( K − l 1 + 5 6 + K − l 2 + 17 6 ) = o ( 1 ) , i = 1 , 2 , … , n , (44) \Xi_i^2\le o(1)\left(K^{-l_1+\frac{5}{6}}+K^{-l_2+\frac{17}{6}}\right)=o(1),\quad i=1,2,\ldots,n, \tag*{(44)} Ξ i 2 ≤ o ( 1 ) ( K − l 1 + 6 5 + K − l 2 + 6 17 ) = o ( 1 ) , i = 1 , 2 , … , n , (44) as l 1 ≥ 1 l_1\ge1 l 1 ≥ 1 l 2 ≥ 3 l_2\ge3 l 2 ≥ 3 K n → ∞ K_n\to\infty K n → ∞ Ξ i 2 ≤ a K 2 \Xi_i^2\le a_K^2 Ξ i 2 ≤ a K 2 i i i
a K 2 = o ( 1 ) ( K − l 1 + 5 6 + K − l 2 + 17 6 ) . (45) a_K^2=o(1)\left(K^{-l_1+\frac{5}{6}}+K^{-l_2+\frac{17}{6}}\right). \tag*{(45)} a K 2 = o ( 1 ) ( K − l 1 + 6 5 + K − l 2 + 6 17 ) . (45) We now show that the contribution of Ξ \Xi Ξ ϑ ∈ R P Z \vartheta\in\mathbb{R}^{P_Z} ϑ ∈ R P Z ϑ ⊤ V ϑ \vartheta^\top V\vartheta ϑ ⊤ V ϑ
ϑ ⊤ ( W ( n ) − W ~ ( n ) ) = 1 n ∑ i = 1 n ϑ ⊤ z i Ξ i . (46) \vartheta^\top\left(W^{(n)}-\widetilde{W}^{(n)}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\vartheta^\top z_i\Xi_i. \tag*{(46)} ϑ ⊤ ( W ( n ) − W ( n ) ) = n 1 i = 1 ∑ n ϑ ⊤ z i Ξ i . (46) Under the i.i.d. setting and the column-norm bound in Lemma A.2 ,
Var { ϑ T ( W ( n ) − W ~ ( n ) ) } = Var ( 1 n ∑ i = 1 n ϑ T z i Ξ i ) = Var ( ϑ T z 1 Ξ 1 ) ≤ E [ ( ϑ T z 1 ) 2 Ξ 1 2 ] ≤ a K 2 E [ ( ϑ T z 1 ) 2 ] = a K 2 ϑ T V ϑ = o ( 1 ) . (47) \begin{aligned}
\operatorname{Var}\left\{\vartheta^{\mathsf{T}}\left(\mathbf{W}^{(n)}-\widetilde{\mathbf{W}}^{(n)}\right)\right\}
&=\operatorname{Var}\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\vartheta^{\mathsf{T}}z_i\Xi_i\right)\\
&=\operatorname{Var}\left(\vartheta^{\mathsf{T}}z_1\Xi_1\right)\\
&\le E\left[(\vartheta^{\mathsf{T}}z_1)^2\Xi_1^2\right]\\
&\le a_K^2E\left[(\vartheta^{\mathsf{T}}z_1)^2\right]\\
&=a_K^2\vartheta^{\mathsf{T}}\mathbf{V}\vartheta=o(1). \tag*{(47)}
\end{aligned} Var { ϑ T ( W ( n ) − W ( n ) ) } = Var ( n 1 i = 1 ∑ n ϑ T z i Ξ i ) = Var ( ϑ T z 1 Ξ 1 ) ≤ E [ ( ϑ T z 1 ) 2 Ξ 1 2 ] ≤ a K 2 E [ ( ϑ T z 1 ) 2 ] = a K 2 ϑ T V ϑ = o ( 1 ) . (47) Therefore,
ϑ T ( W ( n ) − W ~ ( n ) ) = o P ( 1 ) . (48) \vartheta^{\mathsf{T}}\left(\mathbf{W}^{(n)}-\widetilde{\mathbf{W}}^{(n)}\right)=o_P(1). \tag*{(48)} ϑ T ( W ( n ) − W ( n ) ) = o P ( 1 ) . (48) Combining this projected negligibility with the projected central limit theorem for W ~ ( n ) \widetilde{\mathbf{W}}^{(n)} W ( n )
W ( n ) = W ⋅ { 1 + o P ( 1 ) } , (49) \mathbf{W}^{(n)}=\mathbf{W}\cdot\{1+o_P(1)\}, \tag*{(49)} W ( n ) = W ⋅ { 1 + o P ( 1 )} , (49) in the sense defined in the lemma.
(b) The vector W 1 ( n ) \mathbf{W}^{(n)}_1 W 1 ( n ) P 1 P_1 P 1 W ( n ) \mathbf{W}^{(n)} W ( n ) V \mathbf{V} V
W 1 ( n ) = W 1 ⋅ { 1 + o P ( 1 ) } , (50) \mathbf{W}^{(n)}_1=\mathbf{W}_1\cdot\{1+o_P(1)\}, \tag*{(50)} W 1 ( n ) = W 1 ⋅ { 1 + o P ( 1 )} , (50) with the same projected interpretation stated in part (b) of the lemma.
Lemma A.4. Under Assumptions 1-2, and provided P ≤ O ( n c P ) P\le O(n^{c_P}) P ≤ O ( n c P ) c P > 0 c_P>0 c P > 0 k > 0 k>0 k > 0 E ( ϵ i 2 k ) < ∞ E(\epsilon_i^{2k})<\infty E ( ϵ i 2 k ) < ∞ P 1 = O ( n c 1 ) P_1=O(n^{c_1}) P 1 = O ( n c 1 ) P 2 = o ( n ( c 2 − c 1 ) k ) P_2=o\left(n^{(c_2-c_1)k}\right) P 2 = o ( n ( c 2 − c 1 ) k )
Proof. The purpose of this lemma is to verify that the dimension of the expanded linear representation satisfies the growth requirements used in the high-dimensional LASSO argument. Condition (a) in Assumption 2 immediately gives P 1 ≤ O ( n c 1 ) P_1\le O(n^{c_1}) P 1 ≤ O ( n c 1 ) O ( K n 2 ) O(K_n^2) O ( K n 2 )
For the inactive part, the largest candidate dictionary is generated by the full set of pairwise interactions. Thus, the number of inactive coefficients is bounded by
P 2 = O ( P 2 ⋅ K n 2 ) . (51) P_2=O(P^2\cdot K_n^2). \tag*{(51)} P 2 = O ( P 2 ⋅ K n 2 ) . (51) :::
Since P ≤ O ( n c P ) P \le O(n^{c_P}) P ≤ O ( n c P ) K n 2 ≤ O ( n c 1 ) K_n^2 \le O(n^{c_1}) K n 2 ≤ O ( n c 1 )
P 2 = O ( P 2 ⋅ K n 2 ) ≤ O ( n 2 c P + c 1 ) . (52) P_2 = O(P^2 \cdot K_n^2) \le O(n^{2c_P+c_1}). \tag*{(52)} P 2 = O ( P 2 ⋅ K n 2 ) ≤ O ( n 2 c P + c 1 ) . (52) Therefore, it is enough to choose k > 0 k > 0 k > 0
2 c P + c 1 < ( c 2 − c 1 ) k . (53) 2c_P + c_1 < (c_2-c_1)k. \tag*{(53)} 2 c P + c 1 < ( c 2 − c 1 ) k . (53) Such a choice is possible because c 2 > c 1 c_2 > c_1 c 2 > c 1 k k k E ( ϵ i 2 k ) < ∞ E(\epsilon_i^{2k}) < \infty E ( ϵ i 2 k ) < ∞ P 2 = o ( n ( c 2 − c 1 ) k ) P_2 = o\left(n^{(c_2-c_1)k}\right) P 2 = o ( n ( c 2 − c 1 ) k ) □ \square □
Lemma A.5. For step 1 estimator β ~ \widetilde{\beta} β we have { β ~ = s β } ⊃ A n ∩ B n \{\widetilde{\beta} =_s \beta\} \supset A_n \cap B_n { β = s β } ⊃ A n ∩ B n where
A n = { ∣ ( V 11 ( n ) ) − 1 W 1 ( n ) ∣ < n ( ∣ β P 1 ∣ − λ n 2 n ∣ ( V 11 ( n ) ) − 1 sign ( β P 1 ) ∣ ) } , B n = { ∣ ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 W 1 ( n ) − W 2 ( n ) ∣ < λ n 2 n ( 1 − ∣ ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ∣ ) } . (54) \begin{aligned}
A_n &= \left\{\left|\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)}\right| < \sqrt{n}\left(\left|\beta_{P_1}\right|-\frac{\lambda_n}{2n}\left|\left(V_{11}^{(n)}\right)^{-1}\operatorname{sign}(\beta_{P_1})\right|\right)\right\},\\
B_n &= \left\{\left|\left(V_{21}^{(n)}\right)\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)}-W_2^{(n)}\right| < \frac{\lambda_n}{2\sqrt{n}}\left(1-\left|\left(V_{21}^{(n)}\right)\left(V_{11}^{(n)}\right)^{-1}\operatorname{sign}(\beta_{P_1})\right|\right)\right\}. \tag*{(54)}
\end{aligned} A n B n = { ( V 11 ( n ) ) − 1 W 1 ( n ) < n ( ∣ β P 1 ∣ − 2 n λ n ( V 11 ( n ) ) − 1 sign ( β P 1 ) ) } , = { ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 W 1 ( n ) − W 2 ( n ) < 2 n λ n ( 1 − ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ) } . (54) All inequalities are understood componentwise.
Proof. Let u ~ = β ~ − β \widetilde{u}=\widetilde{\beta}-\beta u = β − β u ~ \widetilde{u} u
V ( u ~ ) = L ( β ~ ) = ∥ Y − Z β ~ ∥ 2 2 + λ n ∥ β ~ ∥ 1 = ∥ Z β + Ξ + U − Z β ~ ∥ 2 2 + λ n ∣ β ~ ∣ 1 = ∥ Z u ~ − Ξ − U ∥ 2 2 + λ n ∥ u ~ + β ∥ 1 . (55) \begin{aligned}
V(\widetilde{u})=L(\widetilde{\beta})&=\left\|Y-Z\widetilde{\beta}\right\|_2^2+\lambda_n\left\|\widetilde{\beta}\right\|_1\\
&=\left\|Z\beta+\Xi+U-Z\widetilde{\beta}\right\|_2^2+\lambda_n\left|\widetilde{\beta}\right|_1\\
&=\left\|Z\widetilde{u}-\Xi-U\right\|_2^2+\lambda_n\left\|\widetilde{u}+\beta\right\|_1. \tag*{(55)}
\end{aligned} V ( u ) = L ( β ) = Y − Z β 2 2 + λ n β 1 = Z β + Ξ + U − Z β 2 2 + λ n β 1 = ∥ Z u − Ξ − U ∥ 2 2 + λ n ∥ u + β ∥ 1 . (55) Thus, u ~ \widetilde{u} u V ( u ) V(u) V ( u ) V = Ξ + U \mathcal{V}=\Xi+U V = Ξ + U V ( u ) V(u) V ( u )
∥ Z u − Ξ − U ∥ 2 2 = ∥ Z u − V ∥ 2 2 = u T Z T Z u − 2 u T Z T V + V T V . (56) \left\|Zu-\Xi-U\right\|_2^2=\left\|Zu-\mathcal{V}\right\|_2^2=u^{\mathsf{T}}Z^{\mathsf{T}}Zu-2u^{\mathsf{T}}Z^{\mathsf{T}}\mathcal{V}+\mathcal{V}^{\mathsf{T}}\mathcal{V}. \tag*{(56)} ∥ Z u − Ξ − U ∥ 2 2 = ∥ Z u − V ∥ 2 2 = u T Z T Z u − 2 u T Z T V + V T V . (56) Since u T Z T Z u = ( n u ) T V ( n ) ( n u ) u^{\mathsf{T}}Z^{\mathsf{T}}Zu=(\sqrt{n}u)^{\mathsf{T}}V^{(n)}(\sqrt{n}u) u T Z T Z u = ( n u ) T V ( n ) ( n u )
V ( u ) = [ ( n u ) T V ( n ) ( n u ) − 2 u T Z T V ] + λ n ∥ u + β ∥ 1 + V T V = V 1 ( u ) + V 2 ( u ) + V T V . (57) \begin{aligned}
V(u)&=\left[(\sqrt{n}u)^{\mathsf{T}}V^{(n)}(\sqrt{n}u)-2u^{\mathsf{T}}Z^{\mathsf{T}}\mathcal{V}\right]+\lambda_n\|u+\beta\|_1+\mathcal{V}^{\mathsf{T}}\mathcal{V}\\
&=V_1(u)+V_2(u)+\mathcal{V}^{\mathsf{T}}\mathcal{V}. \tag*{(57)}
\end{aligned} V ( u ) = [ ( n u ) T V ( n ) ( n u ) − 2 u T Z T V ] + λ n ∥ u + β ∥ 1 + V T V = V 1 ( u ) + V 2 ( u ) + V T V . (57) Here, V T V \mathbf{V}^{\mathsf{T}}\mathbf{V} V T V u \mathbf{u} u V 1 V_1 V 1
d V 1 d u = 2 n ( V ( n ) n u − W ( n ) ) , (58) \frac{\mathrm{d}V_1}{\mathrm{d}\mathbf{u}}=2\sqrt{n}\left(\mathbf{V}^{(n)}\sqrt{n}\mathbf{u}-\mathbf{W}^{(n)}\right), \tag*{(58)} d u d V 1 = 2 n ( V ( n ) n u − W ( n ) ) , (58) where W ( n ) \mathbf{W}^{(n)} W ( n )
We now construct a candidate solution satisfying the KKT conditions with the correct sign pattern. Partition u \mathbf{u} u W ( n ) \mathbf{W}^{(n)} W ( n )
{ sign ( β P 1 ) u ~ 1 > − ∣ β P 1 ∣ } , (59) \left\{\operatorname{sign}(\beta_{P_1})\widetilde{\mathbf{u}}_1>-\lvert\beta_{P_1}\rvert\right\}, \tag*{(59)} { sign ( β P 1 ) u 1 > − ∣ β P 1 ∣ } , (59) then β ~ P 1 = β P 1 + u ~ 1 \widetilde{\beta}_{P_1}=\beta_{P_1}+\widetilde{\mathbf{u}}_1 β P 1 = β P 1 + u 1 β P 1 \beta_{P_1} β P 1
{ sign ( β P 1 ) u ~ 1 > − ∣ β P 1 ∣ } ⊂ { β ~ P 1 = s β P 1 } . (60) \left\{\operatorname{sign}(\beta_{P_1})\widetilde{\mathbf{u}}_1>-\lvert\beta_{P_1}\rvert\right\}\subset\left\{\widetilde{\beta}_{P_1}=\!_{s}\beta_{P_1}\right\}. \tag*{(60)} { sign ( β P 1 ) u 1 > − ∣ β P 1 ∣ } ⊂ { β P 1 = s β P 1 } . (60) Under this sign pattern and with β ~ P 2 = 0 \widetilde{\beta}_{P_2}=0 β P 2 = 0
V 11 ( n ) n u ~ 1 − W 1 ( n ) = − λ n 2 n sign ( β P 1 ) , ∣ u ~ 1 ∣ < ∣ β P 1 ∣ , − λ n 2 n 1 ≤ V 21 ( n ) n u ~ 1 − W 2 ( n ) ≤ λ n 2 n 1. (61) \begin{aligned}
\mathbf{V}_{11}^{(n)}\sqrt{n}\widetilde{\mathbf{u}}_1-\mathbf{W}_1^{(n)}
&=-\frac{\lambda_n}{2\sqrt{n}}\operatorname{sign}(\beta_{P_1}),\\
\lvert\widetilde{\mathbf{u}}_1\rvert&<\lvert\beta_{P_1}\rvert,\\
-\frac{\lambda_n}{2\sqrt{n}}\mathbf{1}
&\leq\mathbf{V}_{21}^{(n)}\sqrt{n}\widetilde{\mathbf{u}}_1-\mathbf{W}_2^{(n)}
\leq\frac{\lambda_n}{2\sqrt{n}}\mathbf{1}. \tag*{(61)}
\end{aligned} V 11 ( n ) n u 1 − W 1 ( n ) ∣ u 1 ∣ − 2 n λ n 1 = − 2 n λ n sign ( β P 1 ) , < ∣ β P 1 ∣ , ≤ V 21 ( n ) n u 1 − W 2 ( n ) ≤ 2 n λ n 1 . (61) The first equation determines the active perturbation:
n u ~ 1 = ( V 11 ( n ) ) − 1 W 1 ( n ) − λ n 2 n ( V 11 ( n ) ) − 1 sign ( β P 1 ) . (62) \sqrt{n}\widetilde{\mathbf{u}}_1=\left(\mathbf{V}_{11}^{(n)}\right)^{-1}\mathbf{W}_1^{(n)}-\frac{\lambda_n}{2\sqrt{n}}\left(\mathbf{V}_{11}^{(n)}\right)^{-1}\operatorname{sign}(\beta_{P_1}). \tag*{(62)} n u 1 = ( V 11 ( n ) ) − 1 W 1 ( n ) − 2 n λ n ( V 11 ( n ) ) − 1 sign ( β P 1 ) . (62) The event A n A_n A n ∣ u ~ 1 ∣ < ∣ β P 1 ∣ \lvert\widetilde{\mathbf{u}}_1\rvert<\lvert\beta_{P_1}\rvert ∣ u 1 ∣ < ∣ β P 1 ∣
It remains to verify the inactive KKT inequalities. Substituting the expression for n u ~ 1 \sqrt{n}\widetilde{\mathbf{u}}_1 n u 1
V 21 ( n ) n u ~ 1 − W 2 ( n ) = V 21 ( n ) ( V 11 ( n ) ) − 1 W 1 ( n ) − W 2 ( n ) − λ n 2 n V 21 ( n ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) . (63) \begin{aligned}
\mathbf{V}_{21}^{(n)}\sqrt{n}\widetilde{\mathbf{u}}_1-\mathbf{W}_2^{(n)}
&=\mathbf{V}_{21}^{(n)}\left(\mathbf{V}_{11}^{(n)}\right)^{-1}\mathbf{W}_1^{(n)}-\mathbf{W}_2^{(n)}\\
&\quad-\frac{\lambda_n}{2\sqrt{n}}\mathbf{V}_{21}^{(n)}\left(\mathbf{V}_{11}^{(n)}\right)^{-1}\operatorname{sign}(\beta_{P_1}). \tag*{(63)}
\end{aligned} V 21 ( n ) n u 1 − W 2 ( n ) = V 21 ( n ) ( V 11 ( n ) ) − 1 W 1 ( n ) − W 2 ( n ) − 2 n λ n V 21 ( n ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) . (63) By the triangle inequality, the inactive KKT inequalities hold if
∣ ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 W 1 ( n ) − W 2 ( n ) ∣ < λ n 2 n ( 1 − ∣ ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ∣ ) , (64) \left|\left(V_{21}^{(n)}\right)\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)}-W_2^{(n)}\right|<\frac{\lambda_n}{2\sqrt{n}}\left(1-\left|\left(V_{21}^{(n)}\right)\left(V_{11}^{(n)}\right)^{-1}\operatorname{sign}\left(\beta_{P_1}\right)\right|\right), \tag*{(64)} ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 W 1 ( n ) − W 2 ( n ) < 2 n λ n ( 1 − ( V 21 ( n ) ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) ) , (64) which is precisely the event B n B_n B n
Therefore, on A n ∩ B n A_n\cap B_n A n ∩ B n β ~ P 1 = s β P 1 \widetilde{\beta}_{P_1}={}_s\beta_{P_1} β P 1 = s β P 1 β ~ P 2 = 0 \widetilde{\beta}_{P_2}=0 β P 2 = 0
A n ∩ B n ⊂ { β ~ = s β } , (65) A_n\cap B_n\subset\left\{\widetilde{\beta}={}_s\beta\right\}, \tag*{(65)} A n ∩ B n ⊂ { β = s β } , (65) which proves the result.
Proofs
Proof of Theorem 4.1
Proof. The proof proceeds by reducing the selection consistency of the first-step LASSO estimator to the sign-consistency result of Zhao & Yu [3 ] . By Lemma A.5, it is enough to show that the events A n A_n A n B n B_n B n
P ( A n c ) + P ( B n c ) = O ( P Z n k λ n 2 k ) → 0. (66) P\left(A_n^c\right)+P\left(B_n^c\right)=O\left(\frac{P_Zn^k}{\lambda_n^{2k}}\right)\to0. \tag*{(66)} P ( A n c ) + P ( B n c ) = O ( λ n 2 k P Z n k ) → 0. (66) Once this claim is established, Lemma A.5 immediately implies β ~ P Z = s β P Z \widetilde{\beta}_{P_Z}={}_s\beta_{P_Z} β P Z = s β P Z
We now verify the conditions required by Theorem 3 in Zhao & Yu [3 ] . In our notation, these conditions are:
1 n Z j T Z j ≤ M 1 \frac{1}{n}Z_j^{\mathsf T}Z_j\le M_1 n 1 Z j T Z j ≤ M 1 j j j
ρ min ( V 11 ) > c 3 > 0 \rho_{\min}(V_{11})>c_3>0 ρ m i n ( V 11 ) > c 3 > 0
P 1 = O ( n c 1 ) P_1=O(n^{c_1}) P 1 = O ( n c 1 )
n 1 − c 2 2 d ( β P 1 ) ≥ c 4 > 0 n^{\frac{1-c_2}{2}}d(\beta_{P_1})\ge c_4>0 n 2 1 − c 2 d ( β P 1 ) ≥ c 4 > 0
the strong irrepresentable condition holds;
P 2 = o ( n ( c 2 − c 1 ) k ) P_2=o(n^{(c_2-c_1)k}) P 2 = o ( n ( c 2 − c 1 ) k ) k k k E ( ϵ i 2 k ) < ∞ E(\epsilon_i^{2k})<\infty E ( ϵ i 2 k ) < ∞
1 P Z ( λ n n ) 2 k → ∞ . \frac{1}{P_Z}\left(\frac{\lambda_n}{\sqrt{n}}\right)^{2k} \to\infty. P Z 1 ( n λ n ) 2 k → ∞. We check these conditions one by one. Condition (i) follows from Lemma A.2 , which controls the column norms of the expanded design matrix. Conditions (ii) and (iv) are exactly imposed in Assumption 2. Conditions (iii) and (vi) follow from Lemma A.4, which translates the growth of the original covariate dimension and the truncation level into the growth rates of the expanded active and inactive dimensions. Condition (v) is Assumption 3. Finally, Assumption 1 gives finite moments of all orders because the error distribution has an exponential-type tail. Together with the imposed tuning-rate condition on λ n \lambda_n λ n P Z n k / λ n 2 k → 0 P_Z n^k / \lambda_n^{2k} \to0 P Z n k / λ n 2 k → 0
Therefore, by Theorem 3 of Zhao & Yu [3 ] ,
P ( A n c ) + P ( B n c ) = O ( P Z n k λ n 2 k ) → 0. (67) P\left(A_n^c\right) + P\left(B_n^c\right) = O\left(\frac{P_Z n^k}{\lambda_n^{2k}}\right) \to0. \tag*{(67)} P ( A n c ) + P ( B n c ) = O ( λ n 2 k P Z n k ) → 0. (67) Combining this bound with Lemma A.5, we obtain
P ( β ~ P Z ≠ β P Z ) ≤ P ( ( A n ∩ B n ) c ) ≤ P ( A n c ) + P ( B n c ) = O ( P Z n k λ n 2 k ) → 0. (68) P\left(\widetilde{\beta}_{P_Z} \ne\beta_{P_Z}\right) \leq P\left(\left(A_n \cap B_n\right)^c\right) \leq P\left(A_n^c\right) + P\left(B_n^c\right) = O\left(\frac{P_Z n^k}{\lambda_n^{2k}}\right) \to0. \tag*{(68)} P ( β P Z = β P Z ) ≤ P ( ( A n ∩ B n ) c ) ≤ P ( A n c ) + P ( B n c ) = O ( λ n 2 k P Z n k ) → 0. (68) This establishes the sign consistency of the first-step estimator. Since the second-step adaptive group LASSO assigns infinite weights to components estimated as zero in the first step, the same zero pattern is preserved in Step 2. Hence the selected model is sign consistent as stated in Theorem 4.1.
Proof of Theorem 4.2
Proof. (a) We first study the active-set estimator obtained in Step 2. On the sign-consistency event from Step 1, the inactive components are removed by the infinite adaptive weights. Therefore, Step 2 reduces to an adaptive group LASSO problem on the active set.
We define
{ τ α l = w ^ α l sign ( α ^ l ) , l ∈ L , τ b r = w ^ b r b ^ r 1 { ∥ b ^ r ∥ 2 > 0 } ∥ b ^ r ∥ 2 , r ∈ R , τ γ a b = w ^ γ a b γ ^ a b 1 { ∥ γ ^ a b ∥ 2 > 0 } ∥ γ ^ a b ∥ 2 , ( a , b ) ∈ S . (69) \begin{cases}
\tau_{\alpha l} = \widehat{w}_{\alpha l}\operatorname{sign}(\widehat{\alpha}_l), & l \in L,\\
\tau_{b_r} = \widehat{w}_{b_r}\dfrac{\widehat{b}_r\mathbf{1}\{\|\widehat{b}_r\|_2>0\}}{\|\widehat{b}_r\|_2}, & r \in R,\\
\tau_{\gamma_{ab}} = \widehat{w}_{\gamma_{ab}}\dfrac{\widehat{\gamma}_{ab}\mathbf{1}\{\|\widehat{\gamma}_{ab}\|_2>0\}}{\|\widehat{\gamma}_{ab}\|_2}, & (a,b) \in S. \tag*{(69)}
\end{cases} ⎩ ⎨ ⎧ τ α l = w α l sign ( α l ) , τ b r = w b r ∥ b r ∥ 2 b r 1 { ∥ b r ∥ 2 > 0 } , τ γ ab = w γ ab ∥ γ ab ∥ 2 γ ab 1 { ∥ γ ab ∥ 2 > 0 } , l ∈ L , r ∈ R , ( a , b ) ∈ S . (69) which are the sub-gradients of the penalty terms for the active set in Step 2. Under the event C n C_n C n β ^ P 2 = 0 \widehat{\beta}_{P_2}=0 β P 2 = 0 β ^ P 1 \widehat{\beta}_{P_1} β P 1
− 2 Z P 1 T ( Y − Z P 1 β ^ P 1 ) + λ ~ n τ n = 0 , (70) -2Z_{P_1}^{\mathsf{T}}\left(Y-Z_{P_1}\widehat{\beta}_{P_1}\right)+\widetilde{\lambda}_n\tau_n=0, \tag*{(70)} − 2 Z P 1 T ( Y − Z P 1 β P 1 ) + λ n τ n = 0 , (70) where τ n \tau_n τ n Y = Z P 1 β P 1 + U + Ξ Y=Z_{P_1}\beta_{P_1}+U+\Xi Y = Z P 1 β P 1 + U + Ξ
2 ( Z P 1 T Z P 1 ) ( β ^ P 1 − β P 1 ) − 2 Z P 1 T ( U + Ξ ) + λ ~ n τ n = 0 , (71) 2\left(Z_{P_1}^{\mathsf{T}}Z_{P_1}\right)\left(\widehat{\beta}_{P_1}-\beta_{P_1}\right)-2Z_{P_1}^{\mathsf{T}}(U+\Xi)+\widetilde{\lambda}_n\tau_n=0, \tag*{(71)} 2 ( Z P 1 T Z P 1 ) ( β P 1 − β P 1 ) − 2 Z P 1 T ( U + Ξ ) + λ n τ n = 0 , (71) with probability tending to one. Rearranging yields
n ( β ^ P 1 − β P 1 ) = n ( Z P 1 T Z P 1 ) − 1 W 1 ( n ) − 1 2 n ( Z P 1 T Z P 1 ) − 1 λ ~ n τ n = ( V 11 ( n ) ) − 1 W 1 ( n ) − 1 2 n λ ~ n ( V 11 ( n ) ) − 1 τ n , (72) \begin{aligned}
\sqrt{n}\left(\widehat{\beta}_{P_1}-\beta_{P_1}\right)
&=n\left(Z_{P_1}^{\mathsf{T}}Z_{P_1}\right)^{-1}W_1^{(n)}-\frac{1}{2}\sqrt{n}\left(Z_{P_1}^{\mathsf{T}}Z_{P_1}\right)^{-1}\widetilde{\lambda}_n\tau_n\\
&=\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)}-\frac{1}{2\sqrt{n}}\widetilde{\lambda}_n\left(V_{11}^{(n)}\right)^{-1}\tau_n, \tag*{(72)}
\end{aligned} n ( β P 1 − β P 1 ) = n ( Z P 1 T Z P 1 ) − 1 W 1 ( n ) − 2 1 n ( Z P 1 T Z P 1 ) − 1 λ n τ n = ( V 11 ( n ) ) − 1 W 1 ( n ) − 2 n 1 λ n ( V 11 ( n ) ) − 1 τ n , (72) where W 1 ( n ) = n − 1 / 2 Z P 1 T ( U + Ξ ) W_1^{(n)}=n^{-1/2}Z_{P_1}^{\mathsf{T}}(U+\Xi) W 1 ( n ) = n − 1/2 Z P 1 T ( U + Ξ )
We first verify that the active Gram matrix is stable. By independence and the boundedness of the basis functions,
E ∥ V 11 ( n ) − V 11 ∥ 2 2 = 1 n 2 E ∥ Z P 1 T Z P 1 − E ( Z P 1 T Z P 1 ) ∥ 2 2 = 1 n 2 ∑ i = 1 n E ∥ z i P 1 T z i P 1 − E [ z i P 1 T z i P 1 ] ∥ 2 2 ≤ C Q K n 4 n n 2 = O ( K n 4 n ) . (73) \begin{aligned}
E\left\|V_{11}^{(n)}-V_{11}\right\|_2^2
&=\frac{1}{n^2}E\left\|Z_{P_1}^{\mathsf{T}}Z_{P_1}-E\left(Z_{P_1}^{\mathsf{T}}Z_{P_1}\right)\right\|_2^2\\
&=\frac{1}{n^2}\sum_{i=1}^{n}E\left\|z_{iP_1}^{\mathsf{T}}z_{iP_1}-E\left[z_{iP_1}^{\mathsf{T}}z_{iP_1}\right]\right\|_2^2\\
&\leq\frac{CQK_n^4n}{n^2}=O\left(\frac{K_n^4}{n}\right). \tag*{(73)}
\end{aligned} E V 11 ( n ) − V 11 2 2 = n 2 1 E Z P 1 T Z P 1 − E ( Z P 1 T Z P 1 ) 2 2 = n 2 1 i = 1 ∑ n E z i P 1 T z i P 1 − E [ z i P 1 T z i P 1 ] 2 2 ≤ n 2 C Q K n 4 n = O ( n K n 4 ) . (73) This term converges to zero when K n = o ( n 1 / 4 ) K_n=o(n^{1/4}) K n = o ( n 1/4 ) ( V 11 ( n ) ) − 1 \left(V_{11}^{(n)}\right)^{-1} ( V 11 ( n ) ) − 1
It remains to control τ n \tau_n τ n A n ∩ B n A_n\cap B_n A n ∩ B n
− 2 Z P 1 T ( Y − Z P 1 β ~ P 1 ) + λ n sign ( β P 1 ) = 0 , (74) -2Z_{P_1}^{\mathsf{T}}\left(Y-Z_{P_1}\widetilde{\beta}_{P_1}\right)+\lambda_n\operatorname{sign}(\beta_{P_1})=0, \tag*{(74)} − 2 Z P 1 T ( Y − Z P 1 β P 1 ) + λ n sign ( β P 1 ) = 0 , (74) where the sub-gradient equals sign ( β P 1 ) (\beta_{P_1}) ( β P 1 )
Solving this equation gives
n ( β ~ P 1 − β P 1 ) = ( V 11 ( n ) ) − 1 W 1 ( n ) − λ n 2 n ( V 11 ( n ) ) − 1 sign ( β P 1 ) . (75) \sqrt{n}\left(\widetilde{\beta}_{P_1}-\beta_{P_1}\right)=\left(V_{11}^{(n)}\right)^{-1}\mathbf{W}_1^{(n)}-\frac{\lambda_n}{2\sqrt{n}}\left(V_{11}^{(n)}\right)^{-1}\operatorname{sign}(\beta_{P_1}). \tag*{(75)} n ( β P 1 − β P 1 ) = ( V 11 ( n ) ) − 1 W 1 ( n ) − 2 n λ n ( V 11 ( n ) ) − 1 sign ( β P 1 ) . (75) We next show that the first-step active coefficients are uniformly separated from zero. For a vector g \mathbf{g} g
P ( ∣ β ~ P 1 − β P 1 ∣ > g ) ≤ ∑ j = 1 P 1 P ( ∣ β ~ j − β j ∣ > g j ) = ∑ j = 1 P 1 P ( n ∣ β ~ j − β j ∣ > n g j ) = ∑ j = 1 P 1 P ( ∣ ζ j ( n ) − λ n 2 n ϖ j ∣ > n g j ) ≤ ∑ j = 1 P 1 P ( ∣ ζ j ( n ) ∣ > n ( g j + λ n 2 n ϖ j ) ) . (76) \begin{aligned}
P\left(\left|\widetilde{\beta}_{P_1}-\beta_{P_1}\right|>\mathbf{g}\right)
&\leq\sum_{j=1}^{P_1}P\left(\left|\widetilde{\beta}_j-\beta_j\right|>g_j\right)\\
&=\sum_{j=1}^{P_1}P\left(\sqrt{n}\left|\widetilde{\beta}_j-\beta_j\right|>\sqrt{n}g_j\right)\\
&=\sum_{j=1}^{P_1}P\left(\left|\zeta_j^{(n)}-\frac{\lambda_n}{2\sqrt{n}}\varpi_j\right|>\sqrt{n}g_j\right)\\
&\leq\sum_{j=1}^{P_1}P\left(\left|\zeta_j^{(n)}\right|>\sqrt{n}\left(g_j+\frac{\lambda_n}{2n}\varpi_j\right)\right). \tag*{(76)}
\end{aligned} P ( β P 1 − β P 1 > g ) ≤ j = 1 ∑ P 1 P ( β j − β j > g j ) = j = 1 ∑ P 1 P ( n β j − β j > n g j ) = j = 1 ∑ P 1 P ( ζ j ( n ) − 2 n λ n ϖ j > n g j ) ≤ j = 1 ∑ P 1 P ( ζ j ( n ) > n ( g j + 2 n λ n ϖ j ) ) . (76) where β P 1 = ( β 1 , ⋯ , β P 1 ) T \beta_{P_1}=(\beta_1,\cdots,\beta_{P_1})^{\mathrm{T}} β P 1 = ( β 1 , ⋯ , β P 1 ) T ( V 11 ( n ) ) − 1 W 1 ( n ) = ζ ( n ) \left(V_{11}^{(n)}\right)^{-1}\mathbf{W}_1^{(n)}=\zeta^{(n)} ( V 11 ( n ) ) − 1 W 1 ( n ) = ζ ( n ) ( V 11 ( n ) ) − 1 sign ( β P 1 ) = ( ϖ 1 , ⋯ , ϖ P 1 ) T \left(V_{11}^{(n)}\right)^{-1}\operatorname{sign}(\beta_{P_1})=(\varpi_1,\cdots,\varpi_{P_1})^{\mathrm{T}} ( V 11 ( n ) ) − 1 sign ( β P 1 ) = ( ϖ 1 , ⋯ , ϖ P 1 ) T ζ j ( n ) \zeta_j^{(n)} ζ j ( n ) 2 k 2k 2 k
P ( ∣ ζ j ( n ) ∣ > n ( g j + λ n 2 n ϖ j ) ) ≤ O ( [ n ( g j + λ n 2 n ϖ j ) ] − 2 k ) = O ( n − k g j − 2 k + n k λ n 2 k ) . (77) \begin{aligned}
P\left(\left|\zeta_j^{(n)}\right|>\sqrt{n}\left(g_j+\frac{\lambda_n}{2n}\varpi_j\right)\right)
&\leq O\left(\left[\sqrt{n}\left(g_j+\frac{\lambda_n}{2n}\varpi_j\right)\right]^{-2k}\right)\\
&=O\left(n^{-k}g_j^{-2k}+\frac{n^k}{\lambda_n^{2k}}\right). \tag*{(77)}
\end{aligned} P ( ζ j ( n ) > n ( g j + 2 n λ n ϖ j ) ) ≤ O ( [ n ( g j + 2 n λ n ϖ j ) ] − 2 k ) = O ( n − k g j − 2 k + λ n 2 k n k ) . (77) Therefore,
P ( ∣ β ~ P 1 − β P 1 ∣ > g ) ≤ O ( K 2 [ max { g j − 2 k } n k + n k λ n 2 k ] ) . (78) P\left(\left|\widetilde{\beta}_{P_1}-\beta_{P_1}\right|>\mathbf{g}\right)\leq O\left(K^2\left[\frac{\max\{g_j^{-2k}\}}{n^k}+\frac{n^k}{\lambda_n^{2k}}\right]\right). \tag*{(78)} P ( β P 1 − β P 1 > g ) ≤ O ( K 2 [ n k max { g j − 2 k } + λ n 2 k n k ] ) . (78) Choose g 1 = g 2 = ⋯ = g P 1 = c 4 2 n − ( 1 − c 2 ) / 2 ≤ 1 2 d ( β P 1 ) g_1=g_2=\cdots=g_{P_1}=\frac{c_4}{2}n^{-(1-c_2)/2}\leq\frac{1}{2}d(\beta_{P_1}) g 1 = g 2 = ⋯ = g P 1 = 2 c 4 n − ( 1 − c 2 ) /2 ≤ 2 1 d ( β P 1 )
P ( ∣ β ~ P 1 − β P 1 ∣ > g ) ≤ O ( K 2 n − k c 2 ) + O ( K 2 n k λ n 2 k ) = o ( 1 ) , (79) P\left(\left|\widetilde{\beta}_{P_1}-\beta_{P_1}\right|>\mathbf{g}\right)\leq O\left(K^2n^{-kc_2}\right)+O\left(K^2\frac{n^k}{\lambda_n^{2k}}\right)=o(1), \tag*{(79)} P ( β P 1 − β P 1 > g ) ≤ O ( K 2 n − k c 2 ) + O ( K 2 λ n 2 k n k ) = o ( 1 ) , (79) where the last equality follows by choosing k k k λ n \lambda_n λ n d ( β ~ P 1 ) ≥ ( c 4 / 2 ) n − ( 1 − c 2 ) / 2 d(\widetilde{\beta}_{P_1})\geq(c_4/2)n^{-(1-c_2)/2} d ( β P 1 ) ≥ ( c 4 /2 ) n − ( 1 − c 2 ) /2
We now return to the adaptive penalty in Step 2. Since all active first-step estimates are bounded away from zero at the preceding rate, the adaptive weights on the active set are at most of order O P ( n ( 1 − c 2 ) / 2 ) O_P(n^{(1-c_2)/2}) O P ( n ( 1 − c 2 ) /2 )
λ ~ n n ( V 11 ( n ) ) − 1 τ n = 1 P 1 ⋅ O P ( λ ~ n n n ( 1 − c 2 ) / 2 ) = 1 P 1 ⋅ O P ( n − c 2 / 2 λ ~ n ) = 1 P 1 ⋅ o P ( 1 ) . (80) \frac{\widetilde{\lambda}_n}{\sqrt{n}}\left(V_{11}^{(n)}\right)^{-1}\tau_n = \mathbf{1}_{P_1}\cdot O_P\left(\frac{\widetilde{\lambda}_n}{\sqrt{n}}n^{(1-c_2)/2}\right) = \mathbf{1}_{P_1}\cdot O_P\left(n^{-c_2/2}\widetilde{\lambda}_n\right) = \mathbf{1}_{P_1}\cdot o_P(1). \tag*{(80)} n λ n ( V 11 ( n ) ) − 1 τ n = 1 P 1 ⋅ O P ( n λ n n ( 1 − c 2 ) /2 ) = 1 P 1 ⋅ O P ( n − c 2 /2 λ n ) = 1 P 1 ⋅ o P ( 1 ) . (80) where the last equality follows from λ ~ n = o ( n c 2 / 2 ) \widetilde{\lambda}_n=o(n^{c_2/2}) λ n = o ( n c 2 /2 )
α ^ L − α L = O P ( n − 1 / 2 ) , ∥ b ^ r − b r ∥ 2 2 = O P ( K n ) , ∥ γ ^ a b − γ a b ∥ 2 2 = O P ( K 2 n ) . (81) \widehat{\alpha}_L-\alpha_L=O_P(n^{-1/2}),\qquad\left\|\widehat{\mathbf b}_r-\mathbf b_r\right\|_2^2=O_P\left(\frac{K}{n}\right),\qquad\left\|\widehat{\boldsymbol{\gamma}}_{ab}-\boldsymbol{\gamma}_{ab}\right\|_2^2=O_P\left(\frac{K^2}{n}\right). \tag*{(81)} α L − α L = O P ( n − 1/2 ) , b r − b r 2 2 = O P ( n K ) , ∥ γ ab − γ ab ∥ 2 2 = O P ( n K 2 ) . (81) Taking the first L L L β ^ P 1 \widehat{\boldsymbol{\beta}}_{P_1} β P 1 α ^ L − α L = O P ( n − 1 / 2 ) = o P ( 1 ) \widehat{\alpha}_L-\alpha_L=O_P(n^{-1/2})=o_P(1) α L − α L = O P ( n − 1/2 ) = o P ( 1 )
(b) We next prove the L 2 L^2 L 2
m ^ r ( x ) = Φ ( x ) b ^ r = ∑ k = 0 K − 1 b ^ r k ϕ k ( x ) . (82) \widehat{m}_r(x)=\Phi(x)\widehat{\mathbf b}_r=\sum_{k=0}^{K-1}\widehat{b}_{rk}\phi_k(x). \tag*{(82)} m r ( x ) = Φ ( x ) b r = k = 0 ∑ K − 1 b r k ϕ k ( x ) . (82) By orthonormality of the Hermite basis, the squared L 2 L^2 L 2 K K K
∥ m r − m ^ r ∥ L 2 2 = ∑ k = 0 K − 1 ( b r k − b ^ r k ) 2 + ∑ k ≥ K b r k 2 ≤ ∑ k = 0 K − 1 ( b r k − b ^ r k ) 2 + o ( 1 ) K − l 1 , (83) \left\|m_r-\widehat{m}_r\right\|_{L^2}^2=\sum_{k=0}^{K-1}\left(b_{rk}-\widehat{b}_{rk}\right)^2+\sum_{k\geq K}b_{rk}^2
\leq\sum_{k=0}^{K-1}\left(b_{rk}-\widehat{b}_{rk}\right)^2+o(1)K^{-l_1}, \tag*{(83)} ∥ m r − m r ∥ L 2 2 = k = 0 ∑ K − 1 ( b r k − b r k ) 2 + k ≥ K ∑ b r k 2 ≤ k = 0 ∑ K − 1 ( b r k − b r k ) 2 + o ( 1 ) K − l 1 , (83) where the last inequality follows from Lemma A.1. By Equation A-51,
∑ k = 0 K − 1 ( b r k − b ^ r k ) 2 = O P ( K n ) . (84) \sum_{k=0}^{K-1}\left(b_{rk}-\widehat{b}_{rk}\right)^2=O_P\left(\frac{K}{n}\right). \tag*{(84)} k = 0 ∑ K − 1 ( b r k − b r k ) 2 = O P ( n K ) . (84) Hence,
∥ m r − m ^ r ∥ L 2 2 ≤ O P ( K n + K − l 1 ) . (85) \left\|m_r-\widehat{m}_r\right\|_{L^2}^2\leq O_P\left(\frac{K}{n}+K^{-l_1}\right). \tag*{(85)} ∥ m r − m r ∥ L 2 2 ≤ O P ( n K + K − l 1 ) . (85) Since l 1 ≥ 1 l_1\geq1 l 1 ≥ 1 K = o ( n 1 / 4 ) K=o(n^{1/4}) K = o ( n 1/4 )
∥ m r − m ^ r ∥ L 2 2 = o P ( 1 ) . (86) \left\|m_r-\widehat{m}_r\right\|_{L^2}^2=o_P(1). \tag*{(86)} ∥ m r − m r ∥ L 2 2 = o P ( 1 ) . (86) This proves the L 2 L^2 L 2 m ^ r \widehat{m}_r m r
(c) The proof for the interactive components is analogous. The tensor-product Hermite basis gives an orthogonal decomposition of the bivariate L 2 L^2 L 2
∥ m a b − m ~ a b ∥ L 2 2 ≤ O P ( K 2 n ) + o ( 1 ) K − l 2 + 1 . (87) \lVert m_{ab}-\widetilde{m}_{ab}\rVert_{L^2}^2 \le O_P\left(\frac{K^2}{n}\right)+o(1)K^{-l_2+1}. \tag*{(87)} ∥ m ab − m ab ∥ L 2 2 ≤ O P ( n K 2 ) + o ( 1 ) K − l 2 + 1 . (87) The first term vanishes because K = o ( n 1 / 4 ) K=o(n^{1/4}) K = o ( n 1/4 ) K 2 / n = o ( 1 ) K^2/n=o(1) K 2 / n = o ( 1 ) l 2 ≥ 3 l_2\ge3 l 2 ≥ 3 ∥ m a b − m ~ a b ∥ L 2 2 = o P ( 1 ) \lVert m_{ab}-\widetilde{m}_{ab}\rVert_{L^2}^2=o_P(1) ∥ m ab − m ab ∥ L 2 2 = o P ( 1 )
(d) Finally, setting K = O ( n ν ) K=O(n^\nu) K = O ( n ν )
∥ m r − m ~ r ∥ L 2 2 ≤ O P ( n − ( 1 − ν ) + n − l 1 ν ) . (88) \lVert m_r-\widetilde{m}_r\rVert_{L^2}^2 \le O_P\left(n^{-(1-\nu)}+n^{-l_1\nu}\right). \tag*{(88)} ∥ m r − m r ∥ L 2 2 ≤ O P ( n − ( 1 − ν ) + n − l 1 ν ) . (88) For the interactive components,
∥ m a b − m ~ a b ∥ L 2 2 ≤ O P ( n − ( 1 − 2 ν ) + n − l 2 ν + ν ) . (89) \lVert m_{ab}-\widetilde{m}_{ab}\rVert_{L^2}^2 \le O_P\left(n^{-(1-2\nu)}+n^{-l_2\nu+\nu}\right). \tag*{(89)} ∥ m ab − m ab ∥ L 2 2 ≤ O P ( n − ( 1 − 2 ν ) + n − l 2 ν + ν ) . (89) Taking square roots of the two displays yields the rates stated in the theorem.
Proof of Theorem 4.3
Proof. We start from the expansion in Equation A-42:
n ( β ^ P 1 − β P 1 ) = ( V 11 ( n ) ) − 1 W 1 ( n ) − 1 2 n λ ~ n ( V 11 ( n ) ) − 1 τ n . (90) \sqrt{n}\left(\widehat{\beta}_{P_1}-\beta_{P_1}\right)=\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)}-\frac{1}{2\sqrt{n}}\widetilde{\lambda}_n\left(V_{11}^{(n)}\right)^{-1}\tau_n. \tag*{(90)} n ( β P 1 − β P 1 ) = ( V 11 ( n ) ) − 1 W 1 ( n ) − 2 n 1 λ n ( V 11 ( n ) ) − 1 τ n . (90) The proof applies the central limit theorem to the fixed-dimensional projected statistic appearing in the theorem, rather than to the whole growing-dimensional coefficient vector. This is important because the number of basis coefficients may increase with n n n
We first relate the function-value errors to the coefficient errors. For a specific r 0 r_0 r 0
n ( m ~ r 0 ( x r 0 ) − m r 0 ( x r 0 ) ) = Φ ( x r 0 ) n ( b ^ r 0 − b r 0 ) + n ξ r 0 , (91) \sqrt{n}\left(\widetilde{m}_{r_0}(x_{r_0})-m_{r_0}(x_{r_0})\right)=\Phi(x_{r_0})\sqrt{n}\left(\widehat{b}_{r_0}-b_{r_0}\right)+\sqrt{n}\xi_{r_0}, \tag*{(91)} n ( m r 0 ( x r 0 ) − m r 0 ( x r 0 ) ) = Φ ( x r 0 ) n ( b r 0 − b r 0 ) + n ξ r 0 , (91) where ξ r 0 \xi_{r_0} ξ r 0 x r 0 x_{r_0} x r 0 ξ r 0 = o ( 1 ) K − l 1 / 2 + 5 / 12 \xi_{r_0}=o(1)K^{-l_1/2+5/12} ξ r 0 = o ( 1 ) K − l 1 /2 + 5/12
Since K = O ( n ν ) K = O(n^\nu) K = O ( n ν ) ν ≥ 6 / ( 6 l 1 − 5 ) \nu\ge6/(6l_1 - 5) ν ≥ 6/ ( 6 l 1 − 5 )
n ξ r 0 = o ( 1 ) n 1 2 n − l 1 ν 2 + 5 ν 12 = o ( 1 ) n − ( 6 l 1 − 5 ) ν − 6 12 = o ( 1 ) . (92) \sqrt{n}\xi_{r_0} = o(1)n^{\frac{1}{2}}n^{-\frac{l_1\nu}{2}+\frac{5\nu}{12}} = o(1)n^{-\frac{(6l_1-5)\nu-6}{12}} = o(1). \tag*{(92)} n ξ r 0 = o ( 1 ) n 2 1 n − 2 l 1 ν + 12 5 ν = o ( 1 ) n − 12 ( 6 l 1 − 5 ) ν − 6 = o ( 1 ) . (92) Thus, after multiplication by n \sqrt{n} n
Similarly, for a specific ( a 0 , b 0 ) (a_0,b_0) ( a 0 , b 0 )
n ( m ^ a 0 b 0 ( x a 0 , x b 0 ) − m a 0 b 0 ( x a 0 , x b 0 ) ) = ( Φ ( x a 0 ) ⊗ Φ ( x b 0 ) ) n ( γ ^ a 0 b 0 − γ a 0 b 0 ) + n ξ a 0 b 0 , (93) \begin{aligned}
\sqrt{n}\left(\widehat{m}_{a_0b_0}(x_{a_0},x_{b_0})-m_{a_0b_0}(x_{a_0},x_{b_0})\right) \\
&= \left(\Phi(x_{a_0})\otimes\Phi(x_{b_0})\right)\sqrt{n}\left(\widehat{\gamma}_{a_0b_0}-\gamma_{a_0b_0}\right)+\sqrt{n}\xi_{a_0b_0}, \tag*{(93)}
\end{aligned} n ( m a 0 b 0 ( x a 0 , x b 0 ) − m a 0 b 0 ( x a 0 , x b 0 ) ) = ( Φ ( x a 0 ) ⊗ Φ ( x b 0 ) ) n ( γ a 0 b 0 − γ a 0 b 0 ) + n ξ a 0 b 0 , (93) where ξ a 0 b 0 \xi_{a_0b_0} ξ a 0 b 0 ξ a 0 b 0 = o ( 1 ) K − l 2 / 2 + 17 / 12 \xi_{a_0b_0}=o(1)K^{-l_2/2+17/12} ξ a 0 b 0 = o ( 1 ) K − l 2 /2 + 17/12 K = O ( n ν ) K=O(n^\nu) K = O ( n ν ) ν ≥ 6 / ( 6 l 2 − 17 ) \nu\ge6/(6l_2-17) ν ≥ 6/ ( 6 l 2 − 17 )
n ξ a 0 b 0 = o ( 1 ) n 1 2 n − l 2 ν 2 + 17 ν 12 = o ( 1 ) n − ( 6 l 2 − 17 ) ν − 6 12 = o ( 1 ) . (94) \sqrt{n}\xi_{a_0b_0}=o(1)n^{\frac{1}{2}}n^{-\frac{l_2\nu}{2}+\frac{17\nu}{12}}=o(1)n^{-\frac{(6l_2-17)\nu-6}{12}}=o(1). \tag*{(94)} n ξ a 0 b 0 = o ( 1 ) n 2 1 n − 2 l 2 ν + 12 17 ν = o ( 1 ) n − 12 ( 6 l 2 − 17 ) ν − 6 = o ( 1 ) . (94) Hence the bivariate truncation error is also negligible at the n \sqrt{n} n
Collecting the preceding displays for all selected univariate and interactive components gives
n ( m ^ r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R = diag { ( Φ ( x r j ) ) 1 ≤ j ≤ R } n ( β ^ R − β R ) + o P ( 1 ) , n ( m ^ a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S = diag { ( Φ ( x a j ) ⊗ Φ ( x b j ) ) 1 ≤ j ≤ S } n ( γ ^ S − γ S ) + o P ( 1 ) . (95) \begin{aligned}
\sqrt{n}\left(\widehat{m}_{r_j}(x_{r_j})-m_{r_j}(x_{r_j})\right)_{1\leq j\leq R}
&=\operatorname{diag}\left\{\left(\Phi(x_{r_j})\right)_{1\leq j\leq R}\right\}\sqrt{n}\left(\widehat{\boldsymbol{\beta}}_R-\boldsymbol{\beta}_R\right)+o_P(1),\\
\sqrt{n}\left(\widehat{m}_{a_jb_j}(x_{a_j},x_{b_j})-m_{a_jb_j}(x_{a_j},x_{b_j})\right)_{1\leq j\leq S}
&=\operatorname{diag}\left\{\left(\Phi(x_{a_j})\otimes\Phi(x_{b_j})\right)_{1\leq j\leq S}\right\}\sqrt{n}\left(\widehat{\boldsymbol{\gamma}}_S-\boldsymbol{\gamma}_S\right)+o_P(1). \tag*{(95)}
\end{aligned} n ( m r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R n ( m a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S = diag { ( Φ ( x r j ) ) 1 ≤ j ≤ R } n ( β R − β R ) + o P ( 1 ) , = diag { ( Φ ( x a j ) ⊗ Φ ( x b j ) ) 1 ≤ j ≤ S } n ( γ S − γ S ) + o P ( 1 ) . (95) Therefore,
n ( α ^ L − α L ( m ^ r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m ^ a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S ) = G ( x ) n ( β ^ 1 − β 1 ) + o P ( 1 ) , (96) \sqrt{n}
\begin{pmatrix}
\widehat{\alpha}_L-\alpha_L\\
\left(\widehat{m}_{r_j}(x_{r_j})-m_{r_j}(x_{r_j})\right)_{1\leq j\leq R}\\
\left(\widehat{m}_{a_jb_j}(x_{a_j},x_{b_j})-m_{a_jb_j}(x_{a_j},x_{b_j})\right)_{1\leq j\leq S}
\end{pmatrix}
=G(x)\sqrt{n}\left(\widehat{\boldsymbol{\beta}}_1-\boldsymbol{\beta}_1\right)+o_P(1), \tag*{(96)} n α L − α L ( m r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S = G ( x ) n ( β 1 − β 1 ) + o P ( 1 ) , (96) where
G ( x ) = ( I L O O O diag { ( Φ ( x r j ) ) 1 ≤ j ≤ R } O O O diag { ( Φ ( x a j ) ⊗ Φ ( x b j ) ) 1 ≤ j ≤ S } ) . (97) G(x)=\begin{pmatrix}
I_L & O & O \\
O & \operatorname{diag}\left\{\left(\Phi(x_{r_j})\right)_{1\le j\le R}\right\} & O \\
O & O & \operatorname{diag}\left\{\left(\Phi(x_{a_j})\otimes\Phi(x_{b_j})\right)_{1\le j\le S}\right\}
\end{pmatrix}. \tag*{(97)} G ( x ) = I L O O O diag { ( Φ ( x r j ) ) 1 ≤ j ≤ R } O O O diag { ( Φ ( x a j ) ⊗ Φ ( x b j ) ) 1 ≤ j ≤ S } . (97) The matrix G ( x ) G(x) G ( x ) x x x L + R + S L+R+S L + R + S
Substituting Equation A-42 into the last display gives
n ( α ^ L − α L ( m ^ r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m ^ a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S ) = G ( x ) ( V 11 ( n ) ) − 1 W 1 ( n ) − 1 2 λ ~ n n G ( x ) ( V 11 ( n ) ) − 1 τ n + o P ( 1 ) . (98) \sqrt{n}\begin{pmatrix}
\widehat{\alpha}_L-\alpha_L \\
\left(\widehat{m}_{r_j}(x_{r_j})-m_{r_j}(x_{r_j})\right)_{1\le j\le R} \\
\left(\widehat{m}_{a_jb_j}(x_{a_j},x_{b_j})-m_{a_jb_j}(x_{a_j},x_{b_j})\right)_{1\le j\le S}
\end{pmatrix}
=G(x)\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)}-\frac{1}{2}\frac{\widetilde{\lambda}_n}{\sqrt{n}}G(x)\left(V_{11}^{(n)}\right)^{-1}\tau_n+o_P(1). \tag*{(98)} n α L − α L ( m r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S = G ( x ) ( V 11 ( n ) ) − 1 W 1 ( n ) − 2 1 n λ n G ( x ) ( V 11 ( n ) ) − 1 τ n + o P ( 1 ) . (98) We now study the two terms on the right-hand side separately.
First, since L + R + S L+R+S L + R + S G ( x ) ( V 11 ( n ) ) − 1 W 1 ( n ) G(x)\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)} G ( x ) ( V 11 ( n ) ) − 1 W 1 ( n ) W 1 ( n ) W_1^{(n)} W 1 ( n ) Theorem 4.2 gives V 11 ( n ) → P V 11 V_{11}^{(n)}\xrightarrow{P}V_{11} V 11 ( n ) P V 11
( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1 / 2 G ( x ) ( V 11 ( n ) ) − 1 W 1 ( n ) → d N ( 0 , I L + R + S ) . (99) \left(\sigma^2G(x)V_{11}^{-1}G(x)^{\mathsf T}\right)^{-1/2}G(x)\left(V_{11}^{(n)}\right)^{-1}W_1^{(n)}\xrightarrow{d}N\left(0,I_{L+R+S}\right). \tag*{(99)} ( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1/2 G ( x ) ( V 11 ( n ) ) − 1 W 1 ( n ) d N ( 0 , I L + R + S ) . (99) Second, we show that the adaptive penalty term is negligible under the same normalization.
From the proof of Theorem 4.2 , on the active set the sub-gradient vector satisfies ∥ τ n ∥ 2 = O P ( n ( 1 − c 2 ) / 2 ) \|\tau_n\|_2=O_P\left(n^{(1-c_2)/2}\right) ∥ τ n ∥ 2 = O P ( n ( 1 − c 2 ) /2 ) λ ~ n = o ( n c 2 / 2 ) \widetilde{\lambda}_n=o\left(n^{c_2/2}\right) λ n = o ( n c 2 /2 ) λ ~ n n − 1 / 2 ∥ τ n ∥ 2 = o P ( 1 ) \widetilde{\lambda}_n n^{-1/2}\|\tau_n\|_2=o_P(1) λ n n − 1/2 ∥ τ n ∥ 2 = o P ( 1 ) G ( x ) G(x) G ( x ) ( V 11 ( n ) ) − 1 \left(V_{11}^{(n)}\right)^{-1} ( V 11 ( n ) ) − 1
( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1 / 2 λ ~ n n G ( x ) ( V 11 ( n ) ) − 1 τ n = o P ( 1 ) . (100) \left(\sigma^2G(x)V_{11}^{-1}G(x)^{\mathsf T}\right)^{-1/2}\frac{\widetilde{\lambda}_n}{\sqrt{n}}G(x)\left(V_{11}^{(n)}\right)^{-1}\tau_n=o_P(1). \tag*{(100)} ( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1/2 n λ n G ( x ) ( V 11 ( n ) ) − 1 τ n = o P ( 1 ) . (100) Combining the preceding displays and applying Slutsky’s theorem yields
( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1 / 2 n ( α ^ L − α L ( m ^ r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m ^ a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S ) → d N ( 0 , I L + R + S ) , (101) \left(\sigma^{2}G(x)V_{11}^{-1}G(x)^{\mathsf{T}}\right)^{-1/2}\sqrt{n}
\begin{pmatrix}
\hat{\alpha}_{L}-\alpha_{L}\\
\left(\hat{m}_{r_j}(x_{r_j})-m_{r_j}(x_{r_j})\right)_{1\leq j\leq R}\\
\left(\hat{m}_{a_jb_j}(x_{a_j},x_{b_j})-m_{a_jb_j}(x_{a_j},x_{b_j})\right)_{1\leq j\leq S}
\end{pmatrix}
\xrightarrow{d}N\left(0,I_{L+R+S}\right), \tag*{(101)} ( σ 2 G ( x ) V 11 − 1 G ( x ) T ) − 1/2 n α ^ L − α L ( m ^ r j ( x r j ) − m r j ( x r j ) ) 1 ≤ j ≤ R ( m ^ a j b j ( x a j , x b j ) − m a j b j ( x a j , x b j ) ) 1 ≤ j ≤ S d N ( 0 , I L + R + S ) , (101) which proves the result. □ \square □
Subgradient Descent Algorithm for S-LASSO
We first find the classical LASSO solution by minimizing
min α , b , γ ∥ Y − ∑ l = 0 P α l X l − ∑ r = 1 P Φ ( X r ) b r − ∑ a = 1 P − 1 ∑ b > a P Φ ( X a ) ⊗ Φ ( X b ) γ a b ∥ 2 2 + λ n ( ∑ l = 0 P ∣ α l ∣ + ∑ r = 1 P ∣ b r ∣ + ∑ a = 1 P − 1 ∑ b > a P ∣ γ a b ∣ ) . (102) \min_{\alpha,b,\gamma}\left\|\mathbf{Y}-\sum_{l=0}^{P}\alpha_l\mathbf{X}_l-\sum_{r=1}^{P}\Phi(\mathbf{X}_r)b_r-\sum_{a=1}^{P-1}\sum_{b>a}^{P}\Phi(\mathbf{X}_a)\otimes\Phi(\mathbf{X}_b)\gamma_{ab}\right\|_2^2
+\lambda_n\left(\sum_{l=0}^{P}|\alpha_l|+\sum_{r=1}^{P}|b_r|+\sum_{a=1}^{P-1}\sum_{b>a}^{P}|\gamma_{ab}|\right). \tag*{(102)} α , b , γ min Y − l = 0 ∑ P α l X l − r = 1 ∑ P Φ ( X r ) b r − a = 1 ∑ P − 1 b > a ∑ P Φ ( X a ) ⊗ Φ ( X b ) γ ab 2 2 + λ n ( l = 0 ∑ P ∣ α l ∣ + r = 1 ∑ P ∣ b r ∣ + a = 1 ∑ P − 1 b > a ∑ P ∣ γ ab ∣ ) . (102) Then we construct the weights by
w ^ α l = { N L ∣ α ~ l ∣ − 1 , ∣ α ~ l ∣ > 0 , ∞ , ∣ α ~ l ∣ = 0 , w ^ b r = { N R ∥ b ~ r ∥ 2 − 1 , ∥ b ~ r ∥ 2 > 0 , ∞ , ∥ b ~ r ∥ 2 = 0 , \hat{w}_{\alpha_l}=
\begin{cases}
\sqrt{N_L}|\tilde{\alpha}_l|^{-1}, & |\tilde{\alpha}_l|>0,\\
\infty, & |\tilde{\alpha}_l|=0,
\end{cases}
\qquad
\hat{w}_{b_r}=
\begin{cases}
\sqrt{N_R}\|\tilde{b}_r\|_2^{-1}, & \|\tilde{b}_r\|_2>0,\\
\infty, & \|\tilde{b}_r\|_2=0,
\end{cases} w ^ α l = { N L ∣ α ~ l ∣ − 1 , ∞ , ∣ α ~ l ∣ > 0 , ∣ α ~ l ∣ = 0 , w ^ b r = { N R ∥ b ~ r ∥ 2 − 1 , ∞ , ∥ b ~ r ∥ 2 > 0 , ∥ b ~ r ∥ 2 = 0 , w ^ γ a b = { N S ∥ γ ~ a b ∥ 2 − 1 , ∥ γ ~ a b ∥ 2 > 0 , ∞ , ∥ γ ~ a b ∥ 2 = 0. (103) \hat{w}_{\gamma_{ab}}=
\begin{cases}
\sqrt{N_S}\|\tilde{\gamma}_{ab}\|_2^{-1}, & \|\tilde{\gamma}_{ab}\|_2>0,\\
\infty, & \|\tilde{\gamma}_{ab}\|_2=0.
\end{cases} \tag*{(103)} w ^ γ ab = { N S ∥ γ ~ ab ∥ 2 − 1 , ∞ , ∥ γ ~ ab ∥ 2 > 0 , ∥ γ ~ ab ∥ 2 = 0. (103) where α ~ \tilde{\alpha} α ~ b ~ \tilde{b} b ~ γ ~ \tilde{\gamma} γ ~ N L = L N_L=L N L = L N R = R × K N_R=R\times K N R = R × K N S = S × K 2 N_S=S\times K^2 N S = S × K 2
Then the adaptive group LASSO objective function is
L ( α , b , γ ) = ∥ Y − ∑ l = 0 P α l X l − ∑ r = 1 P Φ ( X r ) b r − ∑ a = 1 P − 1 ∑ b > a P Φ ( X a ) ⊗ Φ ( X b ) γ a b ∥ 2 2 + λ ~ n ( ∑ l = 0 P w ^ α l ∣ α l ∣ + ∑ r = 1 P w ^ b r ∥ b r ∥ 2 + ∑ a = 1 P − 1 ∑ b > a P w ^ γ a b ∥ γ a b ∥ 2 ) . (104) L(\alpha,b,\gamma)=\left\|\mathbf{Y}-\sum_{l=0}^{P}\alpha_l\mathbf{X}_l-\sum_{r=1}^{P}\Phi(\mathbf{X}_r)b_r-\sum_{a=1}^{P-1}\sum_{b>a}^{P}\Phi(\mathbf{X}_a)\otimes\Phi(\mathbf{X}_b)\gamma_{ab}\right\|_2^2
+\tilde{\lambda}_n\left(\sum_{l=0}^{P}\hat{w}_{\alpha_l}|\alpha_l|+\sum_{r=1}^{P}\hat{w}_{b_r}\|b_r\|_2+\sum_{a=1}^{P-1}\sum_{b>a}^{P}\hat{w}_{\gamma_{ab}}\|\gamma_{ab}\|_2\right). \tag*{(104)} L ( α , b , γ ) = Y − l = 0 ∑ P α l X l − r = 1 ∑ P Φ ( X r ) b r − a = 1 ∑ P − 1 b > a ∑ P Φ ( X a ) ⊗ Φ ( X b ) γ ab 2 2 + λ ~ n ( l = 0 ∑ P w ^ α l ∣ α l ∣ + r = 1 ∑ P w ^ b r ∥ b r ∥ 2 + a = 1 ∑ P − 1 b > a ∑ P w ^ γ ab ∥ γ ab ∥ 2 ) . (104) One may write the traditional OLS part as Y − Z β Y-Z\beta Y − Z β
L ( β ) = ∥ Y − Z β ∥ 2 2 + λ ~ n ( ∑ l = 0 P w ~ α l ∣ α l ∣ + ∑ r = 1 P w ~ b r ∥ b r ∥ 2 + ∑ a = 1 P − 1 ∑ b > a P w ~ γ a b ∥ γ a b ∥ 2 ) = L O L S ( β ) + L P e n ( β ) . (105) \begin{aligned}
L(\beta) &= \lVert Y-Z\beta\rVert_2^2+\widetilde{\lambda}_n\left(\sum_{l=0}^{P}\widetilde{w}_{\alpha_l}\lvert\alpha_l\rvert+\sum_{r=1}^{P}\widetilde{w}_{b_r}\lVert b_r\rVert_2+\sum_{a=1}^{P-1}\sum_{b>a}^{P}\widetilde{w}_{\gamma_{ab}}\lVert\gamma_{ab}\rVert_2\right)\\
&=L_{\mathrm{OLS}}(\beta)+L_{\mathrm{Pen}}(\beta). \tag*{(105)}
\end{aligned} L ( β ) = ∥ Y − Z β ∥ 2 2 + λ n ( l = 0 ∑ P w α l ∣ α l ∣ + r = 1 ∑ P w b r ∥ b r ∥ 2 + a = 1 ∑ P − 1 b > a ∑ P w γ ab ∥ γ ab ∥ 2 ) = L OLS ( β ) + L Pen ( β ) . (105) The gradient of L O L S L_{\mathrm{OLS}} L OLS − 2 Z T ( Y − Z β ) -2Z^{\mathsf{T}}(Y-Z\beta) − 2 Z T ( Y − Z β ) L P e n L_{\mathrm{Pen}} L Pen
∂ L P e n ∂ α l = λ ~ n w ~ α l sign ( α l ) , l = 0 , 1 ⋯ , P , ∂ L P e n ∂ b r = λ ~ n w ~ b r b r 1 { ∥ b r ∥ 2 > 0 } ∥ b r ∥ 2 , r = 1 , ⋯ , P , ∂ L P e n ∂ γ a b = λ ~ n w ~ γ r γ a b 1 { ∥ γ a b ∥ 2 > 0 } ∥ γ a b ∥ 2 , 1 ≤ a < b ≤ P . (106) \begin{aligned}
\frac{\partial L_{\mathrm{Pen}}}{\partial\alpha_l}&=\widetilde{\lambda}_n\widetilde{w}_{\alpha_l}\operatorname{sign}(\alpha_l),\quad l=0,1\cdots,P,\\
\frac{\partial L_{\mathrm{Pen}}}{\partial b_r}&=\widetilde{\lambda}_n\widetilde{w}_{b_r}\frac{b_r\mathbf{1}_{\{\lVert b_r\rVert_2>0\}}}{\lVert b_r\rVert_2},\quad r=1,\cdots,P,\\
\frac{\partial L_{\mathrm{Pen}}}{\partial\gamma_{ab}}&=\widetilde{\lambda}_n\widetilde{w}_{\gamma_r}\frac{\gamma_{ab}\mathbf{1}_{\{\lVert\gamma_{ab}\rVert_2>0\}}}{\lVert\gamma_{ab}\rVert_2},\quad1\leq a<b\leq P. \tag*{(106)}
\end{aligned} ∂ α l ∂ L Pen ∂ b r ∂ L Pen ∂ γ ab ∂ L Pen = λ n w α l sign ( α l ) , l = 0 , 1 ⋯ , P , = λ n w b r ∥ b r ∥ 2 b r 1 {∥ b r ∥ 2 > 0 } , r = 1 , ⋯ , P , = λ n w γ r ∥ γ ab ∥ 2 γ ab 1 {∥ γ ab ∥ 2 > 0 } , 1 ≤ a < b ≤ P . (106) where we define 0 0 = 0 \frac{0}{0}=0 0 0 = 0 ∂ L P e n ∂ α l = λ ~ n w ~ α l sign ( α l ) = λ ~ n w ~ α l α l 1 { ∣ α l ∣ > 0 } ∣ α l ∣ \frac{\partial L_{\mathrm{Pen}}}{\partial\alpha_l}=\widetilde{\lambda}_n\widetilde{w}_{\alpha_l}\operatorname{sign}(\alpha_l)=\widetilde{\lambda}_n\widetilde{w}_{\alpha_l}\frac{\alpha_l\mathbf{1}_{\{\lvert\alpha_l\rvert>0\}}}{\lvert\alpha_l\rvert} ∂ α l ∂ L Pen = λ n w α l sign ( α l ) = λ n w α l ∣ α l ∣ α l 1 {∣ α l ∣ > 0 } b r b_r b r γ a b \gamma_{ab} γ ab L L L
subg ( β ) = ∂ L ∂ β = − 2 Z T ( Y − Z β ) + ∂ L P e n ∂ β . (107) \operatorname{subg}(\beta)=\frac{\partial L}{\partial\beta}=-2Z^{\mathsf{T}}(Y-Z\beta)+\frac{\partial L_{\mathrm{Pen}}}{\partial\beta}. \tag*{(107)} subg ( β ) = ∂ β ∂ L = − 2 Z T ( Y − Z β ) + ∂ β ∂ L Pen . (107) Provided the initial value β ( 0 ) \beta^{(0)} β ( 0 )
β ( k ) = β ( k − 1 ) − η l ⋅ subg ( β ( k − 1 ) ) , k ≥ 1 , (108) \beta^{(k)}=\beta^{(k-1)}-\eta_l\cdot\operatorname{subg}\left(\beta^{(k-1)}\right),\quad k\geq1, \tag*{(108)} β ( k ) = β ( k − 1 ) − η l ⋅ subg ( β ( k − 1 ) ) , k ≥ 1 , (108) for a given learning rate η l \eta_l η l ∥ β ( k ) − β ( k − 1 ) ∥ 2 \lVert\beta^{(k)}-\beta^{(k-1)}\rVert_2 ∥ β ( k ) − β ( k − 1 ) ∥ 2 β ^ = β ( k ) \widehat{\beta}=\beta^{(k)} β = β ( k )
Additional Tables, Figures and Results
Firm Characteristics
Table 1 (PDF p. 47). Listing of Firm Characteristics
Empirical Results for Full Specification
In the main empirical analysis, we restrict the prior interaction set to pairs involving firm size (‘mve’). This choice is made for two reasons. First, it is computationally motivated. With P = 102 P = 102 P = 102 102 × 101 / 2 = 5,151 102 \times101/2 = 5{,}151 102 × 101/2 = 5 , 151 [2 ] . This makes ‘mve’ a natural and economically meaningful variable to prioritize when introducing interaction terms.
To assess whether this prior-set implementation omits important interaction patterns, we conduct a supplementary full-specification check. Specifically, we split the full sample period from 1980 to 2019 into ten consecutive four-year subsamples, and in each subsample we estimate S-LASSO with the complete set of pairwise interactions. For each characteristic, we then record how often it is identified as participating in at least one selected interaction effect across the ten subsamples. Table 2 summarizes the characteristics most frequently selected in this exercise.
Table 2 (PDF p. 48). Summary of Characteristics Selected for Interaction Effects. This table reports the characteristics most frequently selected as having interaction effects across the ten sub-samples from 1980–1983 to 2016–2019. The reported frequency indicates the number of sub-samples in which the corresponding characteristic is identified as having at least one interaction effect.
Several conclusions emerge from the table. First, ‘mve’ is selected in all ten subsamples, making it the most frequently connected characteristic in the full specification. Second, ‘indmom’ is selected in nine subsamples, and a small group of variables, including ‘baspread’, ‘mom1m’, and ‘chmom’, are also selected repeatedly. Third, beyond these leading variables, the frequencies decline quickly, indicating that strong interaction patterns are concentrated in a relatively small subset of characteristics rather than being spread broadly across the full 102-dimensional characteristic space.
Overall, these results provide ex post support for the prior interaction set used in the main text. They show that firm size is not only economically well motivated, but also empirically the most central variable in the full-specification interaction network. Therefore, focusing on interactions with ‘mve’ in the main analysis captures the most persistent interaction structure while keeping the rolling-window implementation computationally feasible.
Discussion of the Irrepresentable Condition
Assumption 3 is imposed to establish global sign consistency in Theorem 4.1, which may not hold in the case the covariates are highly correlated. The purpose of this part is therefore not to claim that the condition must literally hold in the empirical application, but to assess how restrictive it appears under dependence structures calibrated to our data.
Empirical Distribution of Pairwise Correlations
As a first step, we examine the correlation structure of the 102 firm characteristics used in the empirical analysis. Figure 1 plots the empirical distribution of the 5,151 pairwise correlations, and Table 3 reports the corresponding summary statistics. The distribution is centered close to zero, with mean 0.0079 and median 0.0003. Most pairwise correlations are small in magnitude, although moderate dependence remains in the tails. These summaries do not verify the strong irrepresentable condition, which depends on the full blockwise correlation structure rather than pairwise correlations alone, but they suggest that the empirical design is not dominated by pervasive near-collinearity.
Figure 1 (PDF p. 50). Density of pairwise correlations among the 102 firm characteristics.
Table 3 (PDF p. 50). Summary statistics of pairwise correlations among the 102 firm characteristics
Simulation Diagnostic
We next conduct a simulation diagnostic in the spirit of Zhao & Yu [3 ] . The goal is to quantify how restrictive the strong irrepresentable condition is under dependence structures comparable to those in our data, and to distinguish between its global and entrywise implications.
Specifically, we generate X ∈ R p \mathbf{X} \in\mathbb{R}^{p} X ∈ R p C p × p C_{p \times p} C p × p C p × p = ( c i j ) C_{p \times p} = (c_{ij}) C p × p = ( c ij )
The first q q q X \mathbf{X} X p − q p-q p − q p ∈ { 50 , 100 , 150 , 200 } p \in\{50, 100, 150, 200\} p ∈ { 50 , 100 , 150 , 200 } q ∈ { 5 , 10 , 15 , 20 } q \in\{5, 10, 15, 20\} q ∈ { 5 , 10 , 15 , 20 } ( p , q ) (p,q) ( p , q ) p = 102 p=102 p = 102
We report diagnostics for two designs. The block labeled Original Design evaluates the strong irrepresentable condition for the raw characteristic vector X \mathbf{X} X Expanded Design repeats the same exercise after transforming X \mathbf{X} X Z \mathbf{Z} Z
For each replication and each design matrix, let A \mathcal{A} A A c \mathcal{A}^{c} A c
r = ∣ V A c A V A A − 1 sign ( β A ) ∣ , (109) \mathbf{r} = \left|V_{\mathcal{A}^{c}\mathcal{A}}V_{\mathcal{A}\mathcal{A}}^{-1}\operatorname{sign}(\boldsymbol{\beta}_{\mathcal{A}})\right|, \tag*{(109)} r = V A c A V AA − 1 sign ( β A ) , (109) where the absolute value is taken componentwise. For an inactive component j ∈ A c j \in\mathcal{A}^{c} j ∈ A c r j r_j r j r \mathbf{r} r
The global satisfaction rate is the fraction of replications in which max j ∈ A c r j < 1 \max_{j \in\mathcal{A}^{c}} r_j < 1 max j ∈ A c r j < 1
The entrywise satisfaction rate is the average fraction of inactive components satisfying
r j < 1 r_j < 1 r j < 1
Table 4 (PDF p. 52). Satisfaction rates of the strong irrepresentable condition under empirically calibrated dependence: diagnostics for the original design and the basis-expanded design used in S-LASSO.
Table 4 reports two complementary measures. Panel A gives the global satisfaction rate , namely, the fraction of replications in which the strong irrepresentable condition holds for the entire inactive set. This is a worst-case requirement: a single violating inactive component is enough for the whole design to fail. Panel B gives the entrywise satisfaction rate , defined as the fraction of inactive components that satisfy the corresponding componentwise inequality. This is only a diagnostic summary and should not be interpreted as a substitute for the formal global condition.
Two findings are worth noting. First, the global satisfaction rate declines as q q q
In particular, the global diagnostic is essentially unchanged, while the entrywise diagnostic is often slightly improved after expansion. Hence, the main difficulty appears to come from the original dependence structure itself rather than from the basis expansion.
