Nonlinear Boosting with Multiple Testing in High-Dimensional Generalised Linear Models with Binary Responses

Nonlinear Boosting with Multiple Testing in High-Dimensional Generalised Linear Models with Binary Responses

Introduction

Statistical models with binary response variables and a large number of potential covariates (regressors) are widely used in economics and finance, applications including, inter alia, default or firm-exit prediction, credit scoring, prediction of the probability of recession or deflation, and analyses of treatment adoption or policy compliance. In such settings, the researcher typically encounters a large number of candidate covariates (possibly even larger than the number of observations available on them), strong cross-sectional dependence, and only a small number of truly relevant covariates that contribute to the response. As a consequence, the problem of variable selection in high-dimensional nonlinear settings is of particular interest and importance. By identifying the subset of relevant covariates, variable selection can improve estimation accuracy and lead to models which are parsimoniously parameterised and easily interpretable.

Existing approaches to high-dimensional model selection are dominated by two classes of methods. The first consists of penalised-likelihood procedures, such as the least absolute shrinkage and selection operator (LASSO) tibshirani1996, the elastic net zou2005regularization and other related approaches fan2001scan,park2007lassologit, to name a few. These methods are computationally attractive and often achieve high true positive rates (sensitivity). However, they are known to have non-negligible false positive rates (fallout) in the presence of pseudo-signals, that is, covariates which do not belong to the data-generating process (DGP) but are correlated with signal covariates that do. Such pseudo-signal effects are extremely common in applications involving macroeconomic and financial time series, many of which are partially driven by a small number of latent common factors. The second class of variable selection approaches comprises sure independent screening procedures fan2008sure,FanSong2010 and multiple-testing approaches such as the one covariate at a time multiple testing (OCMT) approach chudik2018ocmt. While these methods offer explicit error control, they rely on marginal or weakly conditional statistics and, therefore, inherit similar vulnerabilities with respect to pseudo-signal selection.

Building on the boosting with multiple testing (BMT) framework introduced for linear regression models in kapetanios2026, the present paper develops a nonlinear BMT procedure appropriate for high-dimensional generalised linear models (GLMs) with binary responses. The procedure combines two key features. First, selection proceeds in a forward stagewise manner analogous to that associated with boosting BuhlmannHothorn2007, adding at most one covariate per iteration based on a conditional likelihood-based test statistic. Second, at each stage, candidate covariates are filtered using a family-wise multiple-testing threshold similar to that employed in OCMT, which serves both as a screening device and as a stopping rule.

A distinctive feature of the nonlinear BMT method is that it does not rely on sparsity-inducing penalties or marginal screening. Instead, identification and selection are driven by conditional likelihood comparisons that adapt to the current model. This makes the procedure particularly well suited to environments with factor-driven dependence, where many covariates are individually correlated with the response but only a small subset carries independent explanatory power. The procedure is intentionally conservative, aiming to exclude pseudo-signals and noise variables even when covariates are highly collinear. This is partly achieved by selecting a single additional covariate at each stage (based on a marginal significance test), even though many candidate covariates may (individually) have statistically significant coefficients. At each subsequent stage, the remaining covariates are re-assessed conditionally on the expanded model, and again at most one additional covariate is selected.

The main theoretical contribution is to establish finite-stage selection results for nonlinear BMT that parallel, but do not rely on, linear-model arguments. We adopt a general binary-response GLM setting with a smooth link function and allow the number of candidate covariates to grow with the sample size --- while the number of signals remains fixed and finite. We establish three types of results. First, we consider the behaviour of stagewise likelihood-based test statistics along the BMT path in an approximating model that contains all true signals and a finite number of pseudo-signals. Our results establish stochastic control of false selections and justify the stopping rule on which BMT relies. Second, we introduce explicit stagewise dominance conditions under which signals are selected before any pseudo-signal or noise covariates. This yields exact model recovery in both single-signal and multiple-signal cases, in the sense that BMT asymptotically selects all signals and no other covariate with probability tending to one. Third, conditional on correct model selection, we show that the post-BMT maximum-likelihood estimator (MLE) of the GLM parameters is asymptotically equivalent to the infeasible oracle estimator that knows the signals in advance. In other words, the coefficients on the signal variables are estimated as accurately as when the correct sparse model is known.

In addition to theoretical results, we provide extensive simulation evidence that demonstrates the good performance of the nonlinear BMT procedure in finite samples and its superiority over some competing methods. Specifically, in a binary logistic setting with varying degrees of sparsity and collinearity among the covariates, BMT delivers consistently strong performance across a wide range of dimensionalities, combining highly accurate selection of signals, model sizes that closely track the true sparsity level of the DGP and low parameter estimation error. Moreover, BMT outperforms OCMT and LASSO, which tend to be less accurate in selecting signals and yield overparameterised models and higher estimation errors.

We illustrate the practical use of nonlinear BMT by considering a predictive model for the probability that the U.S. inflation rate will, on average, be higher than 2.5% over a 12-month horizon. Using a large set of variables taken from the FRED-MD database of Mccracken2021, BMT selects predictors that yield highly accurate out-of-sample identification of high-inflation periods. This is especially true when the predictors and their estimated coefficients are updated in a recursive manner. While models selected using OCMT and LASSO achieve marginally better in-sample fit, they rely on approximately seven and ten times as many predictors as BMT, respectively, and exhibit poorer out-of-sample performance.

Although the paper focuses on binary-response GLMs, the proposed methodology is considerably more general. Specifically, Appendix (ref) develops a BMT framework for high-dimensional GLMs with responses drawn from an one-parameter exponential family, thereby encompassing a wide range of general models used in applied work. We show that the central theoretical properties of the procedure, including uniform control of the stagewise test statistics, exact model recovery and oracle post-selection inference, continue to hold in this broader setting under suitable conditions. To maintain a clear exposition, the binary-response case is treated in the main text, while the more general theory is relegated to the appendix.

The remainder of the paper is organised as follows. Section (ref) introduces the general GLM framework, defines the stagewise submodels, and formalises the nonlinear BMT procedure. Section (ref) contains the high-level assumptions used in the analysis. Section (ref) provides approximating-model results governing the behaviour of the test statistics and the stopping rule used in the BMT algorithm. Section (ref) contains exact-recovery results for the single-signal and multiple-signal cases, establishes the oracle properties of the post-selection parameter estimator and verifies the required conditions for binary logit and probit models. Section (ref) provides primitive sufficient conditions for a key dominance condition required for exact model recovery and oracle inference, and verifies one of these for binary logit and probit models. Section (ref) provides simulation evidence on the finite-sample properties of nonlinear BMT and competing procedures. Section (ref) discusses an empirical illustration that highlights the usefulness of BMT in a predictive model for the probability of U.S. inflation exceeding a given threshold value over a specified time horizon. Section (ref) summarises and concludes. Proofs and technical lemmas are collected in Appendix (ref). The general exponential-family extension is given in Appendix (ref), Appendix (ref) compares BMT with binary LASSO, and Appendix (ref) gives primitive sufficient conditions for the uniform Wald approximation. Additional simulation results are available from the authors on request.

Setting and nonlinear BMT

This section describes our GLM setup, which follows the classical formulation of mccullagh1989glm. The single-index conditional mean below is maintained for the exact-recovery results. We nevertheless use quasi-likelihood notation, both to permit robust inference and because the stagewise submodels are generally misspecified when they omit active covariates. Their population targets are therefore pseudo-true parameters in the sense of white1982misspec,white94. The stagewise submodels and the BMT algorithm are defined below.

GLM and inference

The parametric model of interest is the single-index GLM

where $y_{t}$ is a response variable taking values in $\{0,1\}$, $\boldsymbol{x}_{t}=(x_{1,t},\dots,x_{p,t})^{\prime}\in\mathbb{R}^{p}$ is a covariate vector whose dimension $p$ may exceed the number of available observations $T$, $\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{p})^{\prime} \in\mathbb{R}^{p}$ is an unknown parameter, and $G$ is a known smooth and monotone (inverse) link function defined on $\mathbb{R}$.\footnote{Without loss of generality, an intercept is not included explicitly in the index $\boldsymbol{x}_{t}^{\prime}\boldsymbol{\beta}$. The covariate corresponding to an intercept is typically treated as an \textquotedblleft always-in\textquotedblright\ control.} In high-dimensional settings where $p$ is large, it is typically assumed that only a small number of covariates contribute to the conditional mean of the response, which amounts to the parameter $\boldsymbol{\beta}$ and its (unknown) true value, denoted $\boldsymbol{\beta}_{0}=(\beta_{0,1},\ldots,\beta_{0,p})^{\prime}\in \mathbb{R}^{p}$, being sparse.

In what follows, $F_{p}=\{1,2,\ldots,p\}$ denotes the set of candidate covariates. Let $S^{(0)}\subseteq F_p$ be a fixed, possibly empty, set of always-in controls. The true active set relevant for selection is the set of active non-always-in covariates,

It is maintained throughout that $S_{0}$ is non-empty with cardinality $k$ that does not depend on $T$. By contrast, the cardinality $p=p_{T}$ of $F_{p}$ is allowed to grow with $T$. A candidate covariate $x_{j,t}$ is said to be a signal variable if $j\in S_{0}$, a pseudo-signal (or proxy) variable if $j\notin S^{(0)}\cup S_{0}$ and $\mathrm{Cov}(x_{j,t},x_{i,t})\neq0$ for some $i\in S_{0}$, and a noise variable if $j\notin S^{(0)}\cup S_{0}$ and $\mathrm{Cov}(x_{j,t},x_{i,t})=0$ for all $i\in S_{0}$.

Inference on the parameter $\boldsymbol{\beta}$ in ((ref)) is typically likelihood-based. The following standard $M$-estimation discussion applies to a fixed-dimensional fitted model, including each admissible stagewise submodel because $k_{\max}$ is fixed. To fix notation, let $l_{t}(\eta_{t},y_{t})$ denote a postulated log-density of $y_{t}$ given $\boldsymbol{x}_{t}$ evaluated at index $\eta_{t}=\eta_{t}(\boldsymbol{\beta })=\boldsymbol{x}_{t}^{\prime}\boldsymbol{\beta}$ and response $y_{t}$; put $\dot{l}_{t}(\eta,y)=\partial l_{t}(\eta,y)/\partial\eta$ and $\ddot{l} _{t}(\eta,y)=\partial^{2}l_{t}(\eta,y)/\partial\eta^{2}$. A (quasi-) MLE $\widehat{\boldsymbol{\beta}}=(\widehat{\beta}_{1},\ldots,\widehat{\beta} _{p})^{\prime}$ is obtained as a maximiser with respect to $\boldsymbol{\beta }$ of the average (quasi-) log-likelihood

For such a fixed-dimensional model, if $l_{t}(\eta_{t},y_{t})$ is correctly specified (which amounts to correct specification of the response probability $\mathbb{P}(y_{t}=1|\boldsymbol{x}_{t})$), the maximiser of the empirical criterion ((ref)) is $\sqrt{T}$-consistent for the corresponding true parameter and asymptotically normal. Under quasi-likelihood or other criterion misspecification, the same conclusion holds for the pseudo-true parameter, denoted $\boldsymbol\beta^*$, under the usual $M$-estimation conditions white1982misspec,white94. When the (quasi-) score sequence is serially uncorrelated, a consistent estimate of the asymptotic covariance matrix of $\sqrt T(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast})$ is the one-period sandwich matrix

where $\boldsymbol{s}_{t}(\boldsymbol{\beta})=\partial l_{t}(\eta_{t} ,y_{t})/\partial\boldsymbol{\beta}=\boldsymbol{x}_{t}\dot{l}_{t}(\eta _{t},y_{t})$ and $\boldsymbol{H}_{t}(\boldsymbol{\beta})=\partial^{2} l_{t}(\eta_{t},y_{t})/\partial\boldsymbol{\beta}\partial\boldsymbol{\beta }^{\prime}=\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{\prime}\ddot{l}_{t}(\eta _{t},y_{t})$. Robustness with respect to dependence in the (quasi-) scores $\boldsymbol{s}_{t}(\boldsymbol{\cdot})$ across observations may be achieved by replacing the second factor in ((ref)) with $T^{-1}\sum _{t=1}^{T}\sum_{i=1}^{T}K(|t-i|/\varpi)\boldsymbol{s}_{t} (\widehat{\boldsymbol{\beta}})\boldsymbol{s}_{i}(\widehat{\boldsymbol{\beta} })^{\prime}$, where $K:\mathbb{R}\rightarrow\lbrack-1,1]$ is a suitable kernel function and $\varpi=\varpi_{T}>0$ is a bandwidth parameter that is an increasing function of $T$. Under the standard GLM assumptions of correct specification and conditionally independent responses, this asymptotic covariance matrix reduces to the inverse of the average observed information $-T^{-1}\sum_{t=1}^{T}\boldsymbol{H}_{t}(\widehat{\boldsymbol{\beta}})$. The covariance matrix of $\widehat{\boldsymbol{\beta}}$ itself is estimated by $T^{-1}$ times the relevant asymptotic covariance matrix; stagewise standard errors below use this $T^{-1/2}$ scaling.

Stagewise submodels and pseudo-true parameters

Let $S\subseteq F_{p}$ denote a conditioning set of selected covariates and $\boldsymbol{x}_{S,t}$ be the subvector of $\boldsymbol{x}_{t}$ indexed by $S$. For a candidate $j\notin S$, the associated augmented submodel is defined as the GLM with parameter $\boldsymbol{\vartheta}_{j} =(\boldsymbol{\gamma}_{j}^{\prime},\theta_{j})^{\prime}\in\mathbb{R}^{|S|+1}$ and index

Letting $L_{T,j}(\boldsymbol{\vartheta}_{j};S)=T^{-1}\sum_{t=1}^{T}l_{t} (\eta_{j,t}(\boldsymbol{\vartheta}_{j};S),y_{t})$ be the average (quasi-) log-likelihood for the augmented submodel, denote a (possibly local) maximiser of it by

and write $\widehat{\theta}_{j}(S)$ for the component of $\widehat{\boldsymbol{\vartheta}}_{j}(S)$ corresponding to $x_{j,t}$.

For theoretical comparisons, we also require the pseudo-true parameter $\boldsymbol{\vartheta}_{j}^{\ast}(S)$, defined as the maximiser of the expected (quasi-) log-likelihood, i.e.,

where the expectation is taken under the true conditional distribution of $y_{t}$ given $\boldsymbol{x}_{t}$; we write $\theta_{j}^{\ast}(S)$ for the component of $\boldsymbol{\vartheta}_{j}^{\ast}(S)$ corresponding to $x_{j,t} $. Under correct specification and when $S\supseteq S^{\dagger}=S^{(0)}\cup S_{0}$, $\theta_{j}^{\ast}(S)=0$ for every $j\notin S^{\dagger}$; under misspecification, $\theta_{j}^{\ast}(S)$ may be nonzero and captures the population \textquotedblleft proxy\textquotedblright\ effect of $x_{j,t}$ after conditioning on $\boldsymbol{x}_{S,t}$.

For subsequent theoretical results, the maximiser in (ref) is assumed to be unique for every admissible pair $(S,j)$. If the population criterion has multiple maximisers, the results should be read under a fixed deterministic tie-breaking rule and the additional requirement that all tied maximisers yield the same candidate coefficient and the same relevant curvature and variance quantities. This convention makes $\theta_j^*(S)$ and the population noncentrality defined below single-valued.

Stagewise Wald statistics and score implementation

For a given $S$ and $j\notin S$, define the (absolute) stagewise Wald statistic for testing $H_{0}:\theta_{j}=0$ in the augmented submodel (ref) as

Throughout the theoretical analysis, $\widehat{\mathrm{se}}_{j}(S)$ is based on an estimator that is uniformly consistent for $T^{-1/2}\{V_{\theta\theta,j}^{\ast}(S)\}^{1/2}$, where $V_{\theta\theta,j}^{\ast}(S)$ is obtained from the long-run sandwich matrix in Assumption (ref). The inverse observed information may be used when the information equality holds for the particular augmented submodel, and the one-period sandwich may be used when its stagewise scores are serially uncorrelated. Otherwise, a uniformly consistent long-run or heteroskedasticity and autocorrelation consistent (HAC) variance estimator is required. The uniform rate needed when population noncentralities diverge is part of Assumption A3.4; pointwise relative consistency alone is not sufficient. Subject to this requirement, the ordering arguments do not depend on the particular implementation of the standard error.

In addition to (ref), we will use a Lagrange multiplier (LM)/score statistic that is convenient for probability bounds. Let $\widehat{\boldsymbol{\vartheta}}_{0}(S)$ denote the restricted (quasi-) MLE obtained subject to $\theta_j=0$, let $\widehat\eta_t(S)=\eta_{j,t}(\widehat{\boldsymbol{\vartheta}}_{0}(S);S)$, and put $\widehat p_t(S)=G\{\widehat\eta_t(S)\}$. A convenient residualised score form is available when the restricted Bernoulli submodel is correctly specified. In that case, define \[ \widehat a_t(S) = \frac{G'\{\widehat\eta_t(S)\}} {\widehat p_t(S)\{1-\widehat p_t(S)\}}, \qquad \widehat w_t(S) = \frac{G'\{\widehat\eta_t(S)\}^{2}} {\widehat p_t(S)\{1-\widehat p_t(S)\}}. \] Let $\widetilde{x}_{j,t}(S)=x_{j,t}-\boldsymbol{x}_{S,t}^{\prime} \widehat{\boldsymbol{\pi}}_{j}(S)$ be the residual from the weighted least-squares projection of $x_{j,t}$ on $\boldsymbol{x}_{S,t}$ using weights $\widehat w_t(S)$. The corresponding normalised score statistic is

For the logit link, $\widehat a_t(S)=1$ and $\widehat w_t(S)=\widehat p_t(S)\{1-\widehat p_t(S)\}$. For a non-canonical link when the restricted submodel is misspecified, the LM statistic should instead be constructed directly from the restricted score, the Schur complement of the actual restricted sample Hessian, and the relevant sandwich or long-run score covariance; the simple weights above need not equal the curvature weights. The formal theory below is Wald-based, so no global score representation is required.

The formal results below are stated for the Wald statistic $W_{T,j}(S)$. If the score statistic in (ref) is used instead, the same arguments apply after imposing the analogue of Assumption (ref).4 with $\widetilde W_{T,j}(S)$ replacing $W_{T,j}(S)$. Appendix (ref) records the standard local Wald--score equivalence for population-null and local alternatives; no global equivalence is required for exact-recovery results.

Nonlinear BMT algorithm

To give a formal description of the BMT procedure, fix an initial set (possibly empty) of \textquotedblleft always-in\textquotedblright\ controls $S^{(0)}$ and a maximum number of stages $k_{\max}$ (with $k_{\max}\geq k$). For $\ell\geq0$, given the current set $S^{(\ell)}$, the nonlinear BMT algorithm comprises four steps.

Theoretical results in subsequent sections establish: (i) high-probability stochastic bounds for the test statistics along the selection path (approximating-model results in Theorems (ref) and (ref)); (ii) sufficient conditions ensuring that the stagewise maximiser is always a signal until $S_{0}$ is exhausted (exact recovery results in Theorems (ref) and (ref)), leading to oracle post-selection inference (Theorem (ref)).

Assumptions

This section collects the assumptions used in the theoretical analysis. They are organised to mirror the logical roles of the main theorems. Assumptions (ref)--(ref) are sufficient for the approximating-model results, while Assumptions (ref)--(ref) introduce additional structure required for exact model recovery and oracle inference. Henceforth, limits in order symbols are taken as $T\rightarrow\infty$, unless stated otherwise.

For later use, let $K=k_{\max}+1$, $d_0=s_0+k_{\max}+1$ and $\ell_T=\log(p\vee T)$, and define the admissible class $\mathcal A_T=\{(S,j):S^{(0)}\subseteq S\subseteq F_p, |S\setminus S^{(0)}|\le k_{\max},\ j\notin S\}$.

For an admissible pair $m=(S,j)\in\mathcal A_T$, write $l_{m,t}(\boldsymbol\vartheta)= l_t\{\eta_{j,t}(\boldsymbol\vartheta;S),y_t\}$. Let $F_{1,m,t}$, $F_{2,m,t}$ and $F_{3,m,t}$ be coordinatewise envelopes, over $\boldsymbol\vartheta\in\mathcal B_{j,S}$ with the compact sets defined in A3.2, for the score, negative Hessian and third derivative of $l_{m,t}$, respectively. Let $F_{\Omega,0,m,t}$ envelope the coordinates of the score outer product and let $F_{\Omega,1,m,t}$ envelope their first parameter derivatives. For $\nu>0$, write $\|Z\|_{\psi_\nu}=\inf\{c>0:\mathbb E\exp[(|Z|/c)^\nu]\le2\}$.

For a given conditioning set $S$, define the population normalised conditional strength, or population Wald noncentrality, of covariate $x_{j,t}$ in the augmented submodel with index (ref) as

The factor $\sqrt{T}$ is essential: $\mathrm{NC}_{j}^{\ast}(S)$ is the population analogue of the absolute Wald statistic, not the unscaled pseudo-true coefficient. When the particular augmented submodel is correctly specified and its observations are conditionally independent, $V_{\theta\theta,j}^{\ast}(S)$ reduces to the corresponding diagonal element of the inverse information matrix.

\paragraph{Remark.} Because BMT adds one covariate at each stage, the exact-recovery induction also holds under the weaker requirement that, for every admissible $S$, the strongest remaining signal satisfies \(\max_{j\in S_0\setminus S}\mathrm{NC}_j^*(S)\ge \max_{m\notin S_0,\,m\notin S}\mathrm{NC}_m^*(S)+d_T\) and \(\max_{j\in S_0\setminus S}\mathrm{NC}_j^*(S)\ge c_T+d_T\), with the same condition $\mathcal E_T/d_T=o_p(1)$. Assumption (ref) uses the minimum over remaining signals and is therefore stronger: it ensures that every remaining signal dominates every non-signal.

Assumption (ref) replaces pointwise ergodicity by dependence conditions that support maximal inequalities. Under A2-PM, the displayed summability condition is satisfied, for example, if $\sup_T\alpha_T(h)\le C h^{-a}$ with $a>q(q+\delta)/(2\delta)$. The fixed values of $s_0$ and $k_{\max}$ imply that every augmented GLM has dimension at most $d_0=s_0+k_{\max}+1$. Moreover, \[ |\mathcal A_T| \le p\sum_{\ell=0}^{k_{\max}}\binom p\ell \le Cp^{k_{\max}+1}=Cp^K. \] Thus the high-dimensional cost is a union bound over polynomially many fixed-dimensional likelihood problems, rather than control over all sparse directions.

Assumption (ref) supplies envelopes for the score, Hessian, score outer product, its first parameter derivative and the third likelihood derivative. For logit links, the scalar score and curvature weights are bounded. Thus, A2-ET follows, for example, from uniformly sub-Gaussian covariates, while A2-PM follows from sufficiently high covariate moments. For probit links, the Mills-ratio terms are unbounded away from compact index sets; a compact-index condition or correspondingly stronger covariate tails are therefore required.

Assumption (ref) remains stated at the level used by the recovery proofs. Appendix (ref) shows that, for information standard errors under the information equality or one-period sandwich standard errors with serially uncorrelated scores, A3.4 follows from A1--A3.3, the envelopes in A2 and the mild parameter-set regularity stated there. Under A2-ET, if $\ell_T^3/T\to0$, then $\mathcal E_T=O_p(\sqrt{\ell_T})$. Under A2-PM, the appendix gives $a_{T,p}^{\mathrm{PM}}=p^{K/q}$ and requires $p^{K/q}=o(\sqrt T)$. The high-level formulation is retained because HAC standard errors add a bandwidth-dependent long-run-variance estimation error. This approximation includes nuisance re-estimation, Hessian approximation and standard-error estimation; pointwise relative consistency of $\widehat{\mathrm{se}}_j(S)$ is not sufficient when population noncentralities diverge.

Assumption (ref), which is only required for exact model recovery, formalises stagewise separation with the correct stochastic margin. The population dominance gap $d_T$ must dominate the maximal sample approximation error $\mathcal E_T$, and the remaining signals must also clear the testing threshold. A fixed positive population-coefficient gap is sufficient in common fixed-signal settings because the Wald noncentrality in (ref) contains the factor $\sqrt T$, but the theorem is stated in the more general rate form needed for high-dimensional testing.

Assumption (ref) ensures that false positives are asymptotically eliminated. It requires the threshold to dominate the maximal uniform stochastic approximation error. Bonferroni or Gaussian-max thresholds must, therefore, be chosen with sufficient slack for (ref) to hold; the boundary choice $c_T=\sqrt{2\log p}$ need not be sufficient without such slack.

Approximating-model results

This section contains results relating to the behaviour of the test statistics and the stopping rule that form the basis of the nonlinear BMT procedure. The results are formulated at the level of the approximating model, that is, a model whose covariates include all signals as well as one or more of the pseudo-signals. No dominance or exact-recovery assumptions are imposed and the conditioning set may be misspecified.

The first result establishes uniform stochastic control of the stagewise Wald test statistics along the BMT path.

\paragraph{Comments.}

Theorem (ref) provides a uniform approximation result for the nonlinear BMT statistics along the entire selection path. Notably, no assumptions are made regarding the relative magnitude of the population strengths $\mathrm{NC}_{j}^{\ast}(S)$ across covariates. The theorem, therefore, applies to arbitrary approximating models, including cases in which strong pseudo-signal effects are present.

Parts (a) and (b) establish that, under either exponential-tail or polynomial-moment conditions, the stagewise Wald statistics concentrate uniformly around their population counterparts. Part (c) shows that the multiple-testing threshold eliminates covariates whose population conditional strength is zero, uniformly across stages.

The results in Theorem (ref) form the basis for the stopping and no-false-inclusion properties of the nonlinear BMT procedure, to which we turn our attention to next. Specifically, we consider two key properties: (i) covariates with zero population conditional strength are not selected (uniformly along the path); (ii) the procedure stops once the approximating model is exhausted, in the sense that no remaining covariate has population conditional strength large enough to cross the multiple-testing threshold.

Before stating the relevant theorem, recall that, for a conditioning set $S$ and a candidate $j\notin S$, $\mathrm{NC}_{j}^{\ast}(S)$ denotes the population normalised conditional strength of $x_{j,t}$ (see (ref)). Define the population-null set at $S$ as

and the maximal remaining population strength at stage $\ell$ as

\paragraph{Comments.}

Part (i) in both regimes is the formal \textquotedblleft no-null\textquotedblright\ property: variables whose population conditional strength is zero are not selected, uniformly over all stages $\ell\leq k_{\max}$. The proof is a direct union-bound argument over $\ell$ combined with Theorem (ref)(c), exploiting that $k_{\max}$ is fixed.

Part (ii) formalises the stopping logic in an approximating-model sense. The condition $M_{\ell^{\dagger}}^{\ast}=o_p(c_{T})$ means that by stage $\ell^{\dagger}$ the remaining population conditional strengths are asymptotically too small to cross the multiple-testing threshold. Since Theorem (ref) controls the maximal deviation of $W_{T,j}(S^{(\ell)})$ from $\mathrm{NC}_{j}^{\ast}(S^{(\ell)})$ uniformly in stage, it follows that the maximal sample statistic falls below $c_{T}$ at (or before) $\ell^{\dagger}$ with probability tending to one.

Exact recovery and oracle inference

This section provides results relating to exact recovery of single or multiple signals by the nonlinear BMT algorithm. The results rely on suitable dominance conditions and yield statements about the recovery of the true model, that is the model containing all signals and no pseudo-signals or noise variables. Conditional on such results, the oracle properties of the post-selection MLE of the GLM parameters can be established.

Single signal

Recall that, at any stage $\ell\geq0$, the BMT procedure compares the stagewise statistics $\{W_{T,j}(S^{(\ell)}):j\notin S^{(\ell)}\}$ and selects the largest one among those exceeding the threshold. Exact recovery in the single-signal case requires that the population conditional strength of the signal dominates that of all other covariates uniformly in the conditioning set. This is formalised in the following dominance condition.

\setcounter{assumption}{0}

Under the restriction in Assumption (ref), the conditioning set before the signal is selected is $S^{(0)}$. The slightly more general notation is retained to match the stagewise formulation used below.

We have the following exact-recovery result in this setting.

\paragraph{Proof outline.}

We sketch the key steps in the proof here; a complete proof is given in Appendix (ref).

Multiple signals

Next, we extend the single-signal recovery result to the general fixed-$k$ case. The key additional difficulty is that, after some signals are selected, the remaining signal variables must continue to dominate all pseudo-signal and noise covariates conditionally on the currently selected set. This is captured by the stagewise dominance condition of Assumption (ref).

The following theorem is the nonlinear analogue of the multiple-signal separation and exact recovery theorem established in linear BMT analysis kapetanios2026.

\paragraph{Comments.}

Part (a) is the central stagewise recovery statement: the algorithm selects a signal at each of the first $k$ stages, hence it necessarily includes all signals after $k$ selections. The nontrivial requirement is the uniform dominance gap in (ref), which must dominate the maximal sample approximation error and must be large enough for the remaining signals to cross the threshold at every intermediate stage.

Part (b) adds a standard \textquotedblleft no-residual-signal\textquotedblright\ condition after $S_{0}$ is included. Under correct specification of the single-index GLM, the condition is satisfied with equality to zero: once $S_{0}$ is included, any covariate $m\notin S_{0}$ has zero population conditional strength in the augmented model. The weaker $o(c_T)$ condition is sufficient for stopping because the threshold dominates the sample approximation error.

\paragraph{Proof outline.}

A complete proof is given in Appendix (ref); here, we summarise the main argument.

Oracle estimation and inference

We now establish the oracle properties of the post-nonlinear-BMT estimator: conditional on exact model recovery, as established in Theorem (ref), the post-selection MLE of the parameters of the GLM behaves asymptotically (to first order) as if the true model that contains only signal variables were known in advance.

Letting $\widehat{S}$ denote the set selected by the nonlinear BMT procedure, write \(S^\dagger=S^{(0)}\cup S_0.\) Define $\widehat{\boldsymbol\beta}_{\rm post}\in\mathbb R^p$ by fitting the (quasi-) MLE on $\widehat S$ and setting all coordinates outside $\widehat S$ equal to zero. Write $\widehat{\boldsymbol\beta}_{{\rm post},S^\dagger}$ for its fixed-dimensional $S^\dagger$ subvector. Let $\boldsymbol\beta_{0,S^\dagger}$ denote the true parameter vector restricted to the oracle support $S^\dagger$, including the always-in controls, and write $\bar d=|S^\dagger|$. The infeasible oracle estimator is computed using the same maximisation and deterministic tie-breaking rule as the post-BMT refit, so the two estimators coincide whenever $\widehat S=S^\dagger$.

\setcounter{assumption}{5}

We have the following result.

\paragraph{Proof outline.}

On the event $\{\widehat S=S^\dagger\}$, which occurs with probability tending to one, the fixed-dimensional post-selection subvector equals the oracle estimator: $\widehat{\boldsymbol\beta}_{{\rm post},S^\dagger} =\widehat{\boldsymbol\beta}_{\rm or}$. Assumption (ref) gives the fixed-dimensional asymptotic linear representation and central limit theorem for the oracle estimator. The high-probability equality between the post-BMT and oracle estimators then transfers the oracle limit distribution to the post-BMT estimator. A complete proof is given in Appendix (ref).

Logit and probit models

We verify that the conditions of Theorems (ref) --(ref) are satisfied for logit and probit models, where the (inverse) link function is the logistic distribution function $\Lambda(z)=1/(1+e^{-z} )$ and the standard normal distribution function $\Phi$, respectively.

\paragraph{Comment.}

For logit, the scalar score and curvature weights are bounded, so the derivative-envelope conditions in Assumption (ref) follow from the corresponding covariate tails. For probit, compactness of the index controls the Mills-ratio terms. Conditions A3.2--A3.3 retain their population identification and nonsingularity roles, while Appendix (ref) provides the uniform Wald approximation. The vanishing of $\mathrm{NC}_{m}^{\ast}(S)$ after inclusion of $S^{\dagger}$ follows from correct specification of the conditional mean.

Together with Theorems (ref)--(ref), the results in Lemma (ref) show that the nonlinear BMT procedure results in consistent model selection and oracle inference for high-dimensional binary-response models under transparent and verifiable conditions.

Primitive conditions for stagewise dominance

In this section, we provide primitive sufficient conditions for the stagewise dominance conditions ((ref)) and ((ref)) required for the exact-recovery results of Theorems (ref) and (ref). We consider two distinct routes relying on assumptions that serve different empirical and conceptual purposes. The local (small-index) regime (Assumptions (ref)--(ref) below) is most appropriate when the single-index remains uniformly moderate, so that the nonlinear GLM admits a stable linear approximation --- in this case dominance reduces to a partial-correlation gap. The global regime is designed for settings with strong dependence and potential proxy covariates, where dominance must be preserved after conditioning on previously selected variables. It imposes stagewise conditional proxy weakness through Assumption (ref). Assumption (ref) gives a primitive moment-based sufficient condition, whereas Assumption (ref) is a partial-projection diagnostic that is useful when additional model-specific restrictions link it to Assumption (ref).

The local conditions and the global-regime condition are independent alternatives. Either Assumptions (ref)--(ref), or Assumption (ref) verified through an appropriate primitive route, is sufficient for stagewise dominance. Theorems (ref) and (ref) require only one of these alternatives. Recalling the definition of the population normalised conditional strength $\mathrm{NC}_{j}^{\ast}(S)$ of covariate $x_{j,t}$, given in (ref), we note that a dominance condition such as Assumption (ref) requires that, for any admissible conditioning set $S$ with $|S\cap S_{0}|\leq k-1$, the remaining signals dominate all non-signals in $\mathrm{NC}_{j}^{\ast}(S)$ by a uniform gap.

Local regime linearisation

The local regime exploits the fact that, when the index $\boldsymbol{x} _{t}^{\prime}\boldsymbol{\beta}_{0}$ is uniformly small, the (inverse) link $G$ can be linearised and the GLM behaves like a weighted linear model. We make the following three assumptions.

\setcounter{assumption}{0}

The next lemma provides the link between the local regime and stagewise dominance.

\paragraph{Comment.}

The local regime reduces dominance to a linear gap in partial correlations with the index. In binary models, it is most plausible when the response probability stays away from $0$ and $1$ and the true coefficients are small enough that $\boldsymbol{x}_{t}^{\prime} \boldsymbol{\beta}_{0}$ remains uniformly moderate. Assumption (ref) is therefore either a fixed small-index condition for bounded designs or a local triangular-array condition. The rate requirement $r_T=o(\underline\kappa_T)$ makes explicit that the linearisation error must be small relative to the partial-correlation gap.

Global regime proxy weakness

Next, we formalise the global regime condition (Assumption (ref)), then give a moment-based sufficient condition and a complementary partial-projection diagnostic.

\setcounter{assumption}{0}

Stagewise dominance follows from ((ref)) by virtue of the next result.

The first primitive form of the global regime condition controls proxies through a conditional score-moment restriction. For an admissible $S$, let $\boldsymbol\gamma^*(S)$ maximise $\mathbb E[l_t(\boldsymbol x_{S,t}'\boldsymbol\gamma,y_t)]$, put $\eta_{S,t}^*=\boldsymbol x_{S,t}'\boldsymbol\gamma^*(S)$, $\xi_t(S)=\dot l_t(\eta_{S,t}^*,y_t)$, and $w_t(S)=-\mathbb E[\ddot l_t(\eta_{S,t}^*,y_t)\mid\boldsymbol x_t]$. Let $\widetilde x_{j,t}(S)$ be the residual defined by $\mathbb E[w_t(S)\boldsymbol x_{S,t}\widetilde x_{j,t}(S)]=\boldsymbol0$. The relevant weighted Gram matrices are assumed nonsingular.

\setcounter{assumption}{0}

The following lemma gives conditions under which ((ref)) is sufficient for ((ref)) to hold.

The second primitive form of the global regime condition controls proxies through a weighted partial-projection restriction. For $R=S_0\setminus S$, let $\widetilde{\boldsymbol x}_{R,t}(S)$ collect the curvature-weighted residuals of the remaining signals after projection on $\boldsymbol x_{S,t}$, and let $\mathcal P^w_{R\mid S}\widetilde x_{m,t}(S)$ denote the fitted value from the population weighted linear projection of $\widetilde x_{m,t}(S)$ on $\widetilde{\boldsymbol x}_{R,t}(S)$.

The bound in (ref) is useful only when it can be shown to deliver the moment attenuation in (ref); by itself it does not control the sign or size of the score moments.

Logit and probit links

Finally, we verify the key analytic condition used in Lemma (ref) for logit and probit links.

In logit and probit models, Assumption (ref), or Assumption (ref) when it is verified to imply the moment attenuation in (ref), implies Assumption (ref), which in turn implies Assumption (ref) (Lemma (ref); see also Remark (ref)). Alternatively, the local small-index conditions ((ref)--(ref)) yield dominance directly (Lemma (ref)). Hence, stagewise dominance is verifiable via multiple primitive routes, depending on whether the setting is closer to a local-linear or global proxy-controlled regime.

\paragraph{Comparison with binary LASSO.} Appendix (ref) compares the stagewise population ordering used by BMT with the generalised irrepresentability condition for binary LASSO. The comparison is not one-sided. Theorem (ref) gives a Gaussian proxy design in which BMT recovers the true support although binary LASSO violates generalised irrepresentability, whereas Theorem (ref) gives a design in which binary LASSO is sign-consistent but an inactive covariate has the largest first-stage BMT statistic. Thus, neither condition implies the other without further restrictions.

Theorem (ref) gives a qualified positive relation. In a local binary model, a strict generalised-irrepresentability margin implies one-signal-at-a-time BMT dominance when conditioning on previously selected signals does not exhaust that margin in the direction of the remaining active score. This directional-stability condition allows correlated active covariates and mixed-sign projection coefficients. The result concerns the strongest remaining signal, which is sufficient for a procedure that adds one covariate at each stage; the stronger all-remaining-signals formulation in Assumption (ref) additionally requires the active-score balance condition stated in Appendix (ref).

Monte Carlo simulations

The purpose of this section is to assess the finite-sample properties of the BMT selection procedure via Monte Carlo simulations. We begin by describing the experimental design and evaluation metrics, and proceed to discuss the results of the simulations.\footnote{Simulations were run on a Linux server (AMD EPYC 7302, 64 threads, 503 GiB RAM, Ubuntu 24.04.3 LTS) using MATLAB.}

Experimental design

The DGP is the binary-response GLM with linear single-index structure given by

where $y_{t}\in\{0,1\}$ are conditionally independent Bernoulli variables, $\boldsymbol{x}_{t}=(x_{1,t},\ldots,x_{p,t})^{\prime}$ is a covariate vector and $T\in\{150,300\}$. The number of covariates is $p\in\{100,200\}$ when $T=150$ and $p\in\{200,400\}$ when $T=300$. The intercept $\alpha\in\{0,1,2\}$ controls the number of observed $0$ and $1$ responses in artificial samples from (ref) --- the larger the value of $\alpha$ is, the more unbalanced the set of binary data is likely to be. For all pairs $(T,p)$ under consideration, $\boldsymbol{\beta}_{0}=2\cdot(\boldsymbol{1}_{k}^{\prime },\boldsymbol{0}_{p-k}^{\prime})^{\prime}$, with $k\in\{1,4\}$, where $\boldsymbol{1}_{n}$ and $\mathbf{0}_{n}$ denote $n\times1$ vectors of all ones and zeros, respectively. The DGP is high-dimensional and sparse, with $k$ signal variables.

Signals $x_{1,t},\ldots,x_{k,t}$ and pseudo-signals $x_{k+1,t},\ldots ,x_{2k,t}$ are generated as

while the remaining noise variables $x_{2k+1,t},\ldots,x_{p,t}$ are obtained as

In (ref)--(ref), $\varepsilon_{j,t}$ are stochastic components that represent idiosyncratic variations across the covariates, while $g_{t}$ and $f_{t}$ are stochastic components that capture local and global common factors, respectively. Note that the local factor $g_{t}$ affects only the signals and pseudo-signals, whereas the global factor $f_{t}$ influences all $p$ covariates. These components are simulated as stationary first-order autoregressive processes with zero mean, unit variance and a common autocovariance structure. Specifically,

where the innovations $\{e_{j,t}^{\varepsilon}\}$, $\{e_{t}^{f}\}$ and $\{e_{t}^{g}\}$ are collections of independent standard normal random variables and are independent of each other. We set $\rho=0.6$.

The inclusion of local and common factor components in the covariates introduces collinearity among the covariates, as is typical in real-world high-dimensional datasets where correlation among signals and pseudo-signals arises as a result of shared latent structures. Additionally, such components reflect the co-movement among covariates that is often observed in big-data settings due to the influence of underlying systemic forces. Examples of such forces include macroeconomic trends, technological innovations, regulatory changes, demographic shifts, geopolitical developments, and sector-wide shocks such as fluctuations in commodity prices, global supply chain disruptions, or monetary policy shifts. Modelling $g_{t}$, $f_{t}$, and $\varepsilon_{j,t}$ as autoregressive processes allows for serial correlation and persistence that are typical of economic time series.

The role of collinearity is assessed by considering a range of values for the parameters $\nu_{g}$ and $\nu_{f}$. Defining the variance inflation factor associated with a signal or pseudo-signal $x_{j,t}$ as $\mathrm{VIF} _{j}=1/(1-\mathrm{R}_{j}^{2})$, where $\mathrm{R}_{j}^{2}$ is the coefficient of determination for the population regression of $x_{j,t}$ on all other $p-1$ candidate covariates, it is easy to show that $\mathrm{VIF}_{j}=1+k(\nu _{f}^{2}+\nu_{g}^{2})\equiv\mathrm{VIF}$ for all $1\leq j\leq2k$. Hence, we set \[ \nu_{f}=\left[ \frac{\omega(\mathrm{VIF}-1)}{k}\right] ^{1/2},\quad\nu _{g}=\left[ \frac{(1-\omega)(\mathrm{VIF}-1)}{k}\right] ^{1/2}, \] for some $\omega\in(0,1)$. When $\mathrm{VIF}=1$, signals and pseudo-signals are mutually uncorrelated and the factors $g_{t}$ and $f_{t}$ do not affect any of the covariates ($\nu_{g}=\nu_{f}=0$). When $\mathrm{VIF}>1$, $100\omega\%$ of the collinearity among covariates, as quantified by $\mathrm{VIF}$, is attributable to the global factor $f_{t}$, while the remaining $100(1-\omega)\%$ arises from the local factor $g_{t}$. We consider $(\mathrm{VIF},\omega)\in\{(1,0),(2,0.25),(2,0.75),(3,0.25),(3,0.75)\}$ in the experiments.

The BMT procedure is implemented as described in Section (ref), using a correctly specified likelihood and family-wise threshold $c_{T}=\Phi^{-1}\{1-0.05/(2aT)\}$, with $a=1$ in the first stage ($\ell=1$) and $a=2$ in subsequent stages ($\ell\geq2$). For comparison, we also include in the experiments the OCMT and LASSO procedures (adapted to the GLM setting). The OCMT implementation uses the same threshold $c_{T}$ as BMT, while the regularisation parameter for the LASSO penalty is chosen by 10-fold cross-validation.

Performance metrics

The properties of covariate-selection procedures are assessed using a variety of metrics. These include: the true positive rate $TPR=TP/(TP+FN)$; the false positive rate $FPR=FP/(FP+TN)$; the true discovery rate $TDR=TP/(TP+FP)$; the false discovery rate $FDR=FP/(FP+TP)$; the $F1$ score: \( F1=\frac{2\cdot TP}{(2\cdot TP)+FP+FN}; \) the Matthews correlation coefficient

\(MCC=\frac{(TP\cdot TN)-(FP\cdot FN)}{[(TP+FP)(TP+FN)(TN+FP)(TN+FN)]^{1/2}}.\) The components of these performance measures are defined as

where $\widehat{D}_{j}=1$ if covariate $x_{j,t}$ is selected and $\widehat{D}_{j}=0$ otherwise, $\beta_{0,j}$ is the true value of the coefficient on $x_{j,t}$ in (ref), and $\mathbb{I}(\cdot)$ is the indicator function.

The $TPR$ metric (sensitivity) captures the share of signals that are correctly selected, whereas $TDR$ (precision) measures the share of selected covariates that are signals. Conversely, $FPR$ is the share of pseudo-signal and noise variables that are erroneously selected, whereas $FDR$ measures the share of selected covariates that are not signals. Large values of $TPR$ and $TDR$ are, evidently, indicative of accurate selection, as are small values of $FPR$ and $FDR$. The $F1$ measure is the harmonic mean of $TPR$ and $TDR$; its lowest value 0 is obtained when no signal variables are selected, while the maximum value 1 is associated with perfectly accurate signal selection. The $MCC$ is, arguably, the most comprehensive accuracy measure, since it synthesises information about all four possible outcomes associated with the quantities in (ref)--(ref); its extreme values of $-1$ and $1$ indicate complete failure and success, respectively, in selecting the signals, while a zero value represents no better than random selection.\footnote{Baldi2020 and Chico2020, among others, provide useful discussions on the relative merits of these accuracy metrics.}

In addition to the above set of classification-based performance measures, we also consider: the number of covariates $(\widehat{k})$ selected (referred to henceforth as model size); the high-dimensional Bayesian information criterion \( BIC=-2\widehat{L}_{T}+\widehat{k}\ln(pT), \) where $\widehat{L}_{T}$ is the maximised log-likelihood associated with the selected $\widehat{k}$-covariate model Chen2012; (iii) the root mean squared error of the MLE, computed as \[ RMSE=\left( N^{-1}\sum_{i=1}^{N}\left\Vert \widehat{\boldsymbol{\beta}} ^{(i)}-\boldsymbol{\beta}_{0}\right\Vert ^{2}\right) ^{1/2}, \] where $\widehat{\boldsymbol{\beta}}^{(i)}$ is the $p$-dimensional vector with $\widehat{k}$ entries equal to the estimated coefficients in the $\widehat{k} $-covariate model selected in the $i$th out of $N$ Monte Carlo replications and all other entries equal to zero. As in the case of the estimation error, the value of each performance criterion for a given design point is computed as an average over $N=1000$ replications.

Simulation results

Our discussion here focuses primarily on the $MCC$, the $F1$ score, the model size, the $BIC$ for the selected models, and the $RMSE$ of the estimated coefficients. For each of the two values of $k$ (number of signals) under consideration, Tables (ref)$-$ (ref) contain summary statistics for the performance measures over all 20 design points (combinations of $p$, $T$ and $(\mathrm{VIF},\omega)$) with $\alpha=0$. The statistics reported include the median (MED), the interquartile range (IQR), the minimum (MIN) and maximum (MAX) values, and the average ranking across the three selection methods (RANK).

It is clear from the results in Table (ref) that BMT outperforms the two competing methods on the $MCC$ metric, ranking first across all experiments on average (mean rank is $1.0$). It attains a median $MCC$ of 0.99, with only a narrow interquartile range in the case of multiple signals. By contrast, the second-ranked OCMT has median $MCC$ of only $0.53$ when $k=1$, rising to $0.76$ when $k=4$. For both values of $k$, the interquartile range for OCMT is substantial. The $MCC$ values for LASSO exhibit low dispersion, by comparison, but they never exceed $0.42$.

The superior performance of BMT is also evident in the summary statistics for the $F1$ score reported in Table (ref). Once again, BMT ranks first in all experiments, with median $F1$ values of $0.99$ for both $k=1$ and $k=4$, and low interquartile range. OCMT does better than LASSO in terms of the median $F1$ value, although this is at the cost of considerable variation.

Another attractive feature of BMT is that selection accuracy is not achieved by sacrificing model parsimony. This can be seen in Table (ref), which contains summary statistics for the number of selected covariates. When there is a single signal, the median model size for BMT is $1.02$ with interquartile range of only $0.01$. OCMT, on the other hand, has median model size $6.30$ and very large dispersion, the maximum number of covariates selected being $325.76$, on average. LASSO also selects many additional covariates beyond the single signal, although no more than $17.48$, on average. Qualitatively similar results are obtained when multiple signals are present in the DGP, LASSO again selecting very large models (a median of $31.11$ covariates), as opposed to the median model size of $7.06$ and $3.97$ of OCMT and BMT, respectively.

Turning to the complexity-penalised measure of fit, Table (ref) records summary statistics for the $BIC$. BMT delivers models with the smallest median $BIC$, regardless of the number of signals present. OCMT is, once again, second best, while the sets of covariates selected by LASSO yield the largest median $BIC$ values.

Summary statistics for the $RMSE$ of the estimated coefficients in the models selected by each of the three procedures are given in Table (ref). Once again, BMT performs best, with a median $RMSE$ of $0.33$ and $0.93$ when $k=1$ and $k=4$, respectively. The corresponding values for OCMT and LASSO are more than twice as large. Although OCMT occasionally attains almost as low an $RMSE$ as that of BMT, there are cases where estimation errors are excessively large (the maximum $RMSE$ values being $1553.86$ for $k=1$ and $166.62$ for $k=4$). In spite of its tendency to select models with a large number of covariates, LASSO does only slight worse than OCMT when $k=1$ and outperforms OCMT when $k=4$.

Similar conclusions emerge from the experiments with $\alpha=1$ and $\alpha=2$. To meet the journal's space constraints, the corresponding summary tables and the design-level results for all nine performance metrics are omitted from the online appendix; they are available from the authors on request. Across these experiments, BMT continues to deliver smaller models, lower parameter-estimation error and higher overall selection accuracy than OCMT and LASSO. It also has a computational advantage over LASSO, as shown in Table (ref), and is faster than OCMT in the very sparse single-signal design.

Empirical illustration

In 2012, the U.S. Federal Reserve (Fed) announced the implementation of an explicit inflation targeting policy, setting the inflation target at 2%, as measured by the annual growth rate of the personal consumption expenditures (PCE) price index. This monetary policy strategy aims to maintain well-anchored inflation expectations and a moderate interest rate environment, thereby contributing to a well-functioning economy. To improve the monitoring of price pressures in the economy, the Federal Reserve Bank of St. Louis has developed a factor-augmented ordered probit model, which is used to forecast the probability that the average annual PCE inflation rate over the next 12 months will exceed 2.5% Jackson2015. These probability forecasts, known as the Price Pressures Measure (PPM), are updated and released on a monthly basis.\footnote{The PPM is available at \url{https://fred.stlouisfed.org/series/STLPPM}.}

Our aim here is to illustrate the practical usefulness of the BMT procedure in the context of forecasting the probability that U.S. inflation will, on average, be higher than the 2.5% PPM value over a 12-month horizon. More specifically, letting $\pi_{t}$ denote the (year-to-year) inflation rate in month $t$, we are interested in the probability of the event $\{\bar{\pi }_{t+12}>2.5\}$, where $\bar{\pi}_{t+12}=(1/12)\sum_{i=1}^{12}\pi_{t+i}$ is the inflation rate averaged over the following 12 months.

Our predictive model for the probabilities associated with the PPM is the binary logit model

where $y_{t+12}=\mathbb{I}(\bar{\pi}_{t+12}>2.5)$ and $\boldsymbol{x}_{t}$ is a $p$-dimensional vector of potential predictors. The latter are taken from the FRED-MD database of macroeconomic and financial variables.\footnote{A complete description of this database can be found in Mccracken2021. The data can be downloaded from \url{https://www.stlouisfed.org/research/economists/mccracken/fred-databases}.} After eliminating variables with missing values and duplicates of consumer price indices, our dataset comprises monthly observations on the PCE inflation rate and on $p=102$ predictors, spanning the period January 1961 to June 2025. The predictors belong to the following general economic groups: output and income (Group 1); labour market (Group 2); housing (Group 3); consumption, orders, and inventories (Group 4); money and credit (Group 5); interest and exchange rates (Group 6); prices (Group 7); stock market (Group 8). When necessary, the transformations suggested by Mccracken2021 are used to ensure approximate stationarity of the data.

To benchmark the performance of BMT, we also implement the OCMT and LASSO variable-selection procedures considered in Section (ref). The predictive ability of the three methods is evaluated using a pseudo-out-of-sample forecasting exercise based on a $85/15$ sample split. In other words, $85\%$ of available observations (January 1961 to December 2015) are used for parameter estimation and variable selection (the training sample), while the remaining $15\%$ (January 2016 to June 2025) are used for forecast evaluation (the evaluation sample). It is worth noting that the evaluation sample includes two policy-relevant episodes: the COVID-19 pandemic in 2020 and the global inflation surge in 2021--2023. This makes the out-of-sample forecasting exercise more challenging and interesting.

As evaluation metrics, we use the $MCC$, the area under the receiver operating characteristic curve ($AUROC$) Bamber1975 and the square root of the quadratic probability score ($SQPS$). The components of the $MCC$ are obtained in this case as

where $\widehat{y}_{t+12}\in(0,1)$ is the fitted (predicted) conditional mean of $y_{t+12}$ from (ref), $\tau\in(0,1)$ is a prespecified threshold value, and $\mathcal{T}$ is the set of all $t$ for which $\widehat{y}_{t+12}$ is computed. We adopt the conventional choice $\tau=0.5$, although other threshold values may also be used. The $SQPS$ and $AUROC$ measures are defined as \[ SQPS=\left( \left\vert \mathcal{T}\right\vert ^{-1}\sum_{t\in\mathcal{T} }\left( \widehat{y}_{t+12}-y_{t+12}\right) ^{2}\right) ^{1/2}, \] and \[ AUROC=\left( \left\vert \mathcal{T}_{0}\right\vert \left\vert \mathcal{T} _{1}\right\vert \right) ^{-1}\sum_{t\in\mathcal{T}_{0}}\sum_{h\in \mathcal{T}_{1}}\mathbb{I}\left( \widehat{y}_{h+12}-\widehat{y} _{t+12}>0\right) , \] where $\mathcal{T}_{i}=\{t\in\mathcal{T}:$ $y_{t+12}=i\}$, $i\in\{0,1\}$. Values of $SQPS$ and $AUROC$ near $0$ and $1$, respectively, are an indication of accurate (in-sample or out-of-sample) identification of high-inflation regimes, that is, periods associated with $\bar{\pi}_{t+12}>2.5$.

In-sample results

In-sample results for the three covariate-selection methods are summarised in Table (ref).\footnote{The BMT, OCMT and LASSO procedures are implemented as in the simulations (see Section (ref)). To account for serial dependence, standard errors are obtained from a kernel HAC covariance estimator with the Bartlett kernel and data-dependent bandwidth determined by the method of Newey1994.} The table reports three measures of predictive performance, namely $MCC$, $AUROC$ and $SQPS$, the Bayesian Information Criterion ($BIC$), and the number of selected predictors $(\widehat{k})$.

All three methods achieve perfect or near-perfect performance with respect to the $AUROC$ metric. This outcome is perhaps unsurprising for OCMT and LASSO, which select 33 and 49 predictors, respectively. By contrast, BMT attains a comparable level of discriminatory accuracy using only 5 predictors.\footnote{These include industrial production measures, capacity utilization, wages and personal consumption expenditures price index.} The resulting gains in parsimony are reflected in the $BIC$, which favours BMT and suggests that it provides a more efficient representation of the predictive information contained in the data

Out-of-sample results

The out-of-sample performance of the variable-selection methods is assessed in the context of two different forecasting schemes. In the first, the estimated parameters of the logit model are updated recursively when generating forecasts (by extending the training sample by one month), while the set of predictors determined by each method in the initial training sample remains the same. By contrast, in the second forecasting scheme, both the selected variables and their estimated coefficients are updated recursively after each out-of-sample forecast. Although this is computationally more expensive, it represents a more realistic approach to out-of-sample evaluation.

Results are reported in Table (ref). BMT consistently outperforms the OCMT and LASSO methods on all three evaluation metrics, irrespective of the forecasting scheme used. The ability of BMT to select predictors that deliver accurate out-of-sample identification of high-inflation regimes is particularly impressive in the second setup, with $AUROC$ around $0.92$.

Figures (ref) and (ref) offer some insight into the differences in the out-of-sample performance of the three variable-selection methods under each of the two forecasting schemes used. The heatmaps in these figures show the number of selected predictors within each of the eight economic groups in our dataset at each point in the evaluation period; white indicates that no variables are selected, while dark blue/red indicates that up to 15/25 variables are selected from a given economic group. As can be seen in Figure (ref), LASSO selects a stable subset of predictors across the evaluation period, but this subset tends to be large. OCMT also selects a large number of predictors across the evaluation period, although not necessarily the same ones. BMT, on the other hand, selects parsimonious models with a stable set of predictors, resulting in lower model and parameter uncertainties compared to OCMT and LASSO. Needless to say, the set of predictors chosen by each selection procedure remains unaltered over the evaluation period in Figure (ref) due to the design of the forecasting scheme.

Based on the results obtained by the three methods considered here, most selected predictors for price pressures in the U.S. fall into the economic groups 1 (output and income), 2 (labour market), 6 (interest/exchange rates) and 7 (prices). These align with traditional specifications of the Phillips curve, which often involve labour market indicators, such as the unemployment rate or unit labour costs, as well as economic activity measures, such as the output gap Stock1999,Faust2013. They also align with popular forecasting models for inflation in which factors such as financial variables Forni2003 and commodity prices Chen2014 have been found to have significant predictive power.\newline

Conclusion

This paper has considered the problem of variable selection in binary-response GLMs when the number of possible covariates is large (and potentially exceeds the number of available observations). Extending the approach of kapetanios2026, we propose a nonlinear BMT procedure to select covariates, combining forward stagewise variable addition (boosting) with a family-wise multiple testing stopping rule. A distinctive feature of the procedure is that covariates are introduced sequentially, with at most one covariate added at each stage based on its conditional statistical significance. Following each addition, all remaining candidate covariates are re-evaluated relative to the expanded model before any further selection takes place. This cautious stagewise strategy helps limit the inclusion of pseudo-signals that are correlated with true signals but do not contribute independently to the conditional mean of the response variable.

We have shown that, under the stated conditions, the nonlinear BMT recovers the true model by selecting all the signals and no other covariates with probability tending to one. What is more, the post-BMT (quasi-) MLE of the parameters of the model enjoys an oracle property, in the sense that it has the same asymptotic properties as the (infeasible) oracle estimator that knows the signal variables in advance.

An extensive simulation study has demonstrated the good finite-sample performance of BMT, and its superiority over OCMT and LASSO, in sparse high-dimensional logit models. Across a variety of designs with varying collinearity patterns and numbers of covariates and signals, BMT was found to deliver parsimonious specifications, highly accurate covariate selection and low parameter estimation errors. Combined with the relatively low computational cost of BMT (relative to competing procedures such as OCMT and LASSO) and the attractive asymptotic properties of the procedure, these findings show that BMT provides a useful and powerful tool for variable selection.

Finally, an empirical application has illustrated the practical usefulness of the BMT procedure in the context of forecasting the probability that U.S. inflation will, on average, be higher than 2.5% over a 12-month horizon. We have found that BMT delivers a relatively small set of predictors that can identify accurately out-of-sample periods in which average inflation is high.