The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
102,920 characters
Nonlinear Boosting with Multiple Testing in High-Dimensional Generalised Linear Models with Binary Responses
\title{Nonlinear Boosting with Multiple Testing in High-Dimensional Generalised
Linear Models with Binary Responses}
\author{\textbf{Charisios Grivas}\thanks{Department of Mathematics, Aalborg
University, Denmark (\href{[email removed]}{[email removed]}).}\\Aalborg University
\and \textbf{George Kapetanios}\thanks{Department of Banking and Finance, King's
College London, U.K.
(\href{[email removed]}{[email removed]}). }\\King's College London
\and \textbf{Zacharias Psaradakis}\thanks{Birkbeck Business School, Birkbeck,
University of London, U.K.
(\href{[email removed]}{[email removed]}).}\\Birkbeck, University of London
\and \textbf{Vasilis Sarafidis}\thanks{Brunel Business School, Brunel University of
London, U.K.
(\href{[email removed]}{[email removed]}).}\\Brunel University London
\and \textbf{Mari\'{a}n V\'{a}vra}\thanks{Research Department, National Bank of
Slovakia and Institute of Forecasting, Slovak Academy of Sciences, Slovak Republic (\href{[email removed]}{[email removed]}
).}\\National Bank of Slovakia\\Slovak Academy of Sciences
\and \textbf{Alexia Ventouri}\thanks{Department of Banking and Finance, King's
College London, U.K.
(\href{[email removed]}{[email removed]}). }\\King's College London}
\date{\today}
\maketitle
\clearpage\newpage
\begin{abstract}
\noindent This paper proposes a nonlinear boosting with multiple testing (BMT)
approach to variable selection in high-dimensional generalised linear models
with binary responses. At each stage of the BMT procedure, the model is
updated by adding only the most significant covariate, conditional on those already selected in previous stages, while taking into account the multiple testing nature of the
problem. It is shown that, under the stated conditions, the BMT procedure selects
all covariates whose true coefficients are nonzero, and no other covariates,
with probability tending to one. Furthermore, the procedure enjoys an
oracle property, in the sense that the post-BMT maximum likelihood estimator
of the parameters of the model is asymptotically equivalent to an oracle
estimator that knows the correct sparse model in advance. Monte Carlo
experiments demonstrate that BMT outperforms competing methods, delivering
high covariate-selection accuracy and low parameter estimation error. An
empirical example illustrates that BMT delivers a predictive model for the
probability that U.S. inflation exceeds a given threshold over a 12-month
horizon which has very good out-of-sample performance.\bigskip\ \newline
\noindent\textit{JEL classification}: C12; C13; C25; C52; C55.\newline
\noindent\textit{Key words}: Binary response; Boosting; Generalised linear
models; High dimensionality; Multiple testing; Oracle estimator; Variable selection.
\end{abstract}
\clearpage\newpage
\section{Introduction}
\label{sec: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)
\citep{tibshirani1996}, the elastic net \citep{zou2005regularization} and
other related approaches \citep[e.g.,][]{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 \citep[e.g.,][]{fan2008sure,FanSong2010} and multiple-testing
approaches such as the one covariate at a time multiple testing (OCMT)
approach \citep{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 \emph{boosting with multiple testing} (BMT) framework
introduced for linear regression models in \citet{kapetanios2026}, the present
paper develops a \emph{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
\citep[e.g.,][]{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 \cite{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{app:addGLM} 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{sec:glm_bmt}
introduces the general GLM framework, defines the stagewise submodels, and
formalises the nonlinear BMT procedure. Section~\ref{sec:assumptions} contains
the high-level assumptions used in the analysis. Section~\ref{sec:approx}
provides approximating-model results governing the behaviour of the test
statistics and the stopping rule used in the BMT algorithm.
Section~\ref{sec:exact}~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{sec:primitive}
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{sec:Monte_Carlo}
provides simulation evidence on the finite-sample properties of nonlinear BMT
and competing procedures. Section~\ref{sec:Example} 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{sec:conclusion} summarises and concludes.
Proofs and technical lemmas are collected in Appendix~\ref{app:prelim}. The general exponential-family extension is given in Appendix~\ref{app:addGLM}, Appendix~\ref{app:binary_lasso} compares BMT with binary LASSO, and Appendix~\ref{app:primitive_wald} gives primitive sufficient conditions for the uniform Wald approximation. Additional simulation results are available from the authors on request.
\section{Setting and nonlinear BMT}
\label{sec:glm_bmt}
This section describes our GLM setup, which follows the classical formulation
of \citet{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
\citet{white1982misspec,white94}. The stagewise submodels and the BMT algorithm
are defined below.
\subsection{GLM and inference}
\label{sec:glm_model}
The parametric model of interest is the single-index GLM
\begin{equation}
\mathbb{E}(y_{t}|\boldsymbol{x}_{t})=G(\boldsymbol{x}_{t}^{\prime
}\boldsymbol{\beta}),\quad t=1,2,\dots,T, \label{eq:glm_single_index}
\end{equation}
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 \emph{true active set} relevant for selection is the
set of active non-always-in covariates,
\begin{equation}
S_{0}=\{j\in F_{p}\setminus S^{(0)}:\beta_{0,j}\neq0\},\qquad
\left\vert S_{0}\right\vert =k.
\label{eq:true_set}
\end{equation}
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
\emph{signal} variable if $j\in S_{0}$, a \emph{pseudo-signal} (or
\emph{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
\emph{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{eq:glm_single_index})
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
\begin{equation}
L_{T}(\boldsymbol{\beta})=T^{-1}\sum_{t=1}^{T}l_{t}(\eta_{t},y_{t}).
\label{eq:criterion}
\end{equation}
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{eq:criterion}) 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 \citep{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
\begin{equation}
\left( T^{-1}\sum_{t=1}^{T}\boldsymbol{H}_{t}(\widehat{\boldsymbol{\beta}
})\right) ^{-1}\left( T^{-1}\sum_{t=1}^{T}\boldsymbol{s}_{t}
(\widehat{\boldsymbol{\beta}})\boldsymbol{s}_{t}(\widehat{\boldsymbol{\beta}
})^{\prime}\right) \left( T^{-1}\sum_{t=1}^{T}\boldsymbol{H}_{t}
(\widehat{\boldsymbol{\beta}})\right) ^{-1},\label{eq:cov_est}
\end{equation}
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{eq:cov_est}) 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.
\subsection{Stagewise submodels and pseudo-true parameters}
\label{sec:pseudotrue}
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 \emph{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
\begin{equation}
\eta_{j,t}(\boldsymbol{\vartheta}_{j};S)=\boldsymbol{x}_{S,t}^{\prime
}\boldsymbol{\gamma}_{j}+\theta_{j}\,x_{j,t}. \label{eq:aug_index}
\end{equation}
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
\begin{equation}
\widehat{\boldsymbol{\vartheta}}_{j}(S)\in\arg\max_{\boldsymbol{\vartheta}
_{j}}L_{T,j}(\boldsymbol{\vartheta}_{j};S), \label{eq:mle_aug}
\end{equation}
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 \emph{pseudo-true} parameter
$\boldsymbol{\vartheta}_{j}^{\ast}(S)$, defined as the maximiser of the
expected (quasi-) log-likelihood, i.e.,
\begin{equation}
\boldsymbol{\vartheta}_{j}^{\ast}(S)\in\arg\max_{\boldsymbol{\vartheta}_{j}
}\;\mathbb{E}\bigl[l_{t}(\eta_{j,t}(\boldsymbol{\vartheta}_{j};S),y_{t}
)\bigr], \label{eq:pseudotrue}
\end{equation}
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 \eqref{eq:pseudotrue} 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.
\subsection{Stagewise Wald statistics and score implementation}
\label{sec:wald}
For a given $S$ and $j\notin S$, define the (absolute) \emph{stagewise Wald
statistic} for testing $H_{0}:\theta_{j}=0$ in the augmented submodel
\eqref{eq:aug_index} as
\begin{equation}
W_{T,j}(S)=\frac{\bigl|\widehat{\theta}_{j}(S)\bigr|}{\widehat{\mathrm{se}
}_{j}(S)}. \label{eq:wald_stat}
\end{equation}
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{ass:A3}. 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 \eqref{eq:wald_stat}, 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
\begin{equation}
\widetilde{W}_{T,j}(S)
=
\frac{\left|T^{-1/2}\sum_{t=1}^{T}\widetilde{x}_{j,t}(S)
\widehat a_t(S)\{y_t-\widehat p_t(S)\}\right|}
{\left[T^{-1}\sum_{t=1}^{T}\widetilde{x}_{j,t}^{2}(S)
\widehat w_t(S)\right]^{1/2}}.
\label{eq:score_form}
\end{equation}
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 \eqref{eq:score_form} is used instead, the same
arguments apply after imposing the analogue of Assumption~\ref{ass:A3}.4 with
$\widetilde W_{T,j}(S)$ replacing $W_{T,j}(S)$. Appendix~\ref{app:wald_score}
records the standard local Wald--score equivalence for population-null and
local alternatives; no global equivalence is required for exact-recovery
results.
\subsection{Nonlinear BMT algorithm}
\label{sec:algo}
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.
\begin{enumerate}
\item \emph{Restricted fit}\textbf{.} Estimate the parameters of a GLM with
covariates indexed by $S^{(\ell)}$ and obtain the fitted objects of interest
(residual-like terms and weights).
\item \emph{Candidate screening}\textbf{.} For each $j\notin S^{(\ell)}$,
compute the stagewise statistic $W_{T,j}\bigl(S^{(\ell)}\bigr)$. A score
statistic $\widetilde{W}_{T,j}\bigl(S^{(\ell)}\bigr)$ may be used instead
when the uniform approximation condition is imposed for that statistic.
\item \emph{Multiple-testing filter}\textbf{.} For a family-wise threshold
$c_{T}>0$, form the passing set
\begin{equation}
\mathcal{J}^{(\ell)}=\Bigl\{j\notin S^{(\ell)}:W_{T,j}\bigl(S^{(\ell
)}\bigr)\geq c_{T}\Bigr\}. \label{eq:pass_set}
\end{equation}
\item \emph{Stagewise selection}\textbf{.} If $\mathcal{J}^{(\ell
)}=\varnothing$, stop and return $\widehat{S}=S^{(\ell)}$. Otherwise, select
\[
j^{(\ell)}\in\arg\max_{j\in\mathcal{J}^{(\ell)}}W_{T,j}\bigl(S^{(\ell
)}\bigr),\qquad S^{(\ell+1)}=S^{(\ell)}\cup\{j^{(\ell)}\},
\]
and continue until the algorithm stops or $k_{\max}$ successful selections
have been made.
\end{enumerate}
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{thm:approx_stats} and
\ref{thm:approx_stop}); (ii) sufficient conditions ensuring that the stagewise
maximiser is always a signal until $S_{0}$ is exhausted (exact recovery
results in Theorems~\ref{thm:k1_recovery} and \ref{thm:fixedk_recovery}),
leading to oracle post-selection inference (Theorem~\ref{thm:oracle}).
\section{Assumptions}
\label{sec: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{ass:A1}--\ref{ass:A3} are sufficient for the approximating-model results,
while Assumptions \ref{ass:A4}--\ref{ass:A5} 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.
\begin{assumption} \label{ass:A1}\qquad
\begin{enumerate}
\item[A1.1] The triangular array
$\{(y_t,\boldsymbol x_t')\}_{t=1}^T$ is row-wise strictly stationary and
strongly mixing. Let $\alpha_T(h)$ denote its strong-mixing coefficient,
uniformly over rows. Under A2-ET, there are constants
$C_\alpha,c_\alpha>0$ such that
$\sup_T\alpha_T(h)\le C_\alpha\exp(-c_\alpha h)$ for every $h\ge1$.
Under A2-PM, with $q$ and $\delta$ as defined there,
\[
\sum_{h=1}^{\infty}h^{q/2-1}
\sup_T\alpha_T(h)^{\delta/(q+\delta)}<\infty.
\]
\item[A1.2] The number of covariates $|F_p|=p=p_T$ may diverge with $T$,
while the number of non-always-in signals $|S_0|=k$ and the maximum number of
successful selections $k_{\max}$ are fixed, positive and finite. The number
of always-in covariates $s_0=|S^{(0)}|$ is fixed and finite and may be zero.
The always-in set and the signal set are disjoint:
$S^{(0)}\cap S_0=\varnothing$.\medskip
\end{enumerate}
\end{assumption}
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\}$.
\begin{assumption}
\label{ass:A2} One of the following sets of conditions is satisfied:
\begin{enumerate}
\item[A2-ET] Uniformly over admissible $m$ and $t$,
\[
\|F_{1,m,t}\|_{\psi_2}
+\|F_{2,m,t}\|_{\psi_1}
+\|F_{\Omega,0,m,t}\|_{\psi_1}
+\|F_{3,m,t}\|_{\psi_{2/3}}
+\|F_{\Omega,1,m,t}\|_{\psi_{2/3}}
\le C_F<\infty.
\]
\item[A2-PM] There exist $q>\max\{4,d_0\}$ and $\delta>0$ such that,
uniformly over admissible $m$ and $t$,
\[
\mathbb E\!\left(
F_{1,m,t}^{q+\delta}+F_{2,m,t}^{q+\delta}
+F_{\Omega,0,m,t}^{q+\delta}+F_{3,m,t}^{q+\delta}
+F_{\Omega,1,m,t}^{q+\delta}
\right)\le C_F<\infty.
\]
\end{enumerate}
\end{assumption}
\begin{assumption}
\label{ass:A3}\qquad
\begin{enumerate}
\item[A3.1] The (inverse) link function
$G:\mathbb{R}\rightarrow(0,1)$ is strictly increasing and three-times
continuously differentiable, with derivatives that are bounded on compact
subsets of $\mathbb{R}$.
\item[A3.2] For each admissible pair $(S,j)$, there is a compact parameter
set $\mathcal B_{j,S}\subset\mathbb R^{|S|+1}$ containing
$\boldsymbol\vartheta_j^*(S)$ in its interior. Uniformly over all stagewise
submodels with $|S\setminus S^{(0)}|\le k_{\max}$ and all $j\notin S$, the
norms of the pseudo-true parameters and the diameters of the sets
$\mathcal B_{j,S}$ are bounded, and the distances of the pseudo-true
parameters from the boundaries of $\mathcal B_{j,S}$ are bounded away from
zero. Writing
$Q_{j,S}(\boldsymbol\vartheta)=\mathbb E
l_{m,t}(\boldsymbol\vartheta)$, the maximisers are uniformly well separated:
for every $\epsilon>0$ there is $c(\epsilon)>0$ such that
\[
\inf_{(S,j)\in\mathcal A_T}
\left[
Q_{j,S}\{\boldsymbol\vartheta_j^*(S)\}
-
\sup_{\substack{\boldsymbol\vartheta\in\mathcal B_{j,S}:\\
\|\boldsymbol\vartheta-\boldsymbol\vartheta_j^*(S)\|\ge\epsilon}}
Q_{j,S}(\boldsymbol\vartheta)
\right]\ge c(\epsilon).
\]
\item[A3.3] For any admissible
conditioning set $S$ with $|S\setminus S^{(0)}|\leq k_{\max}$ and any
$j\notin S$, the population information matrix
\[
\boldsymbol{J}_{j}^{\ast}(S)
= -\mathbb{E}\left[
\frac{\partial^{2}l_{t}(\eta_{j,t}(\boldsymbol{\vartheta}_{j}^{\ast}(S);S),y_t)}
{\partial\boldsymbol{\vartheta}_{j}\partial\boldsymbol{\vartheta}_{j}^{\prime}}
\right]
\]
is finite and nonsingular, with smallest eigenvalue uniformly bounded away
from zero. The long-run score covariance matrix
\[
\boldsymbol{\Omega}_{j}^{\ast}(S)
=
\sum_{h=-\infty}^{\infty}
\mathbb{E}\!\left[
\boldsymbol{s}_{j,t}(\boldsymbol{\vartheta}_{j}^{\ast}(S);S)
\boldsymbol{s}_{j,t-h}(\boldsymbol{\vartheta}_{j}^{\ast}(S);S)^{\prime}
\right]
\]
exists and is nonsingular uniformly, where
\(
\boldsymbol{s}_{j,t}(\boldsymbol{\vartheta}_{j};S)
=
\frac{\partial l_t(\eta_{j,t}(\boldsymbol{\vartheta}_{j};S),y_t)}
{\partial\boldsymbol{\vartheta}_{j}}.
\)
The sandwich matrix
\(
\boldsymbol{V}_{j}^{\ast}(S)
=
\boldsymbol{J}_{j}^{\ast}(S)^{-1}
\boldsymbol{\Omega}_{j}^{\ast}(S)
\boldsymbol{J}_{j}^{\ast}(S)^{-1}
\)
is finite, and its candidate-coefficient diagonal element
$V_{\theta\theta,j}^{\ast}(S)$ is uniformly bounded away from zero and infinity.
If the stagewise scores are serially uncorrelated, the displayed long-run matrix
$\boldsymbol{\Omega}_{j}^{\ast}(S)$ reduces to the usual one-period covariance.
\item[A3.4] (High-level uniform Wald approximation.) For the admissible
class $\mathcal A_T$ defined above, define
\[
\mathcal E_T
=
\max_{(S,j)\in\mathcal A_T}
\left|W_{T,j}(S)-\mathrm{NC}_{j}^{\ast}(S)\right|.
\]
The following uniform approximation rates hold. Under Assumption A2-ET,
\(\mathcal E_T=O_p(\sqrt{\ell_T}).\)
Under Assumption A2-PM, there exists a deterministic maximal-deviation rate
$a_{T,p}^{\mathrm{PM}}$ such that
\(\mathcal E_T=O_p(a_{T,p}^{\mathrm{PM}}).\)
Primitive sufficient conditions and explicit rates are given in
Appendix~\ref{app:primitive_wald}.
\end{enumerate}
\end{assumption}
For a given conditioning set $S$, define the population \emph{normalised
conditional strength}, or population Wald noncentrality, of covariate
$x_{j,t}$ in the augmented submodel with index \eqref{eq:aug_index} as
\begin{equation}
\mathrm{NC}_{j}^{\ast}(S)
=
\frac{\sqrt{T}\,|\theta_{j}^{\ast}(S)|}
{\{V_{\theta\theta,j}^{\ast}(S)\}^{1/2}}.
\label{eq:nc_def_prim}
\end{equation}
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.
\begin{assumption}
\label{ass:A4}
There exists a deterministic sequence $d_T>0$ such that
\begin{equation}
\frac{\mathcal E_T}{d_T}=o_p(1),
\label{eq:dom_error_margin}
\end{equation}
and, for any conditioning set $S$ satisfying
$S^{(0)}\subseteq S\subseteq S^{(0)}\cup S_0$ and $|S\cap S_0|<k$,
\begin{equation}
\min_{j\in S_{0}\setminus S}\mathrm{NC}_{j}^{\ast}(S)
\ge
\max_{\substack{m\notin S_{0}\\m\notin S}}
\mathrm{NC}_{m}^{\ast}(S)+d_T.
\label{eq:stg_dom}
\end{equation}
In addition, every remaining signal is above the testing threshold by the same
population margin:
\begin{equation}
\min_{j\in S_{0}\setminus S}\mathrm{NC}_{j}^{\ast}(S)
\ge c_T+d_T.
\label{eq:stg_threshold}
\end{equation}
\end{assumption}
\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{ass:A4}
uses the minimum over remaining signals and is therefore stronger: it ensures
that every remaining signal dominates every non-signal.
\begin{assumption}
\label{ass:A5}
The threshold sequence $\{c_T\}$ used in the nonlinear BMT procedure satisfies
$c_T\rightarrow\infty$ and $c_T=o(\sqrt T)$. Moreover, there exists a constant
$\bar{\eta}_{0}\in(0,1)$ such that
\begin{equation}
\Pr\!\left(\mathcal E_T\le (1-\bar{\eta}_{0})c_T\right)\rightarrow1,\quad\mathrm{as}\text{ }T\rightarrow\infty.
\label{eq:threshold_calibration}
\end{equation}
\end{assumption}
Assumption \ref{ass:A1} 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{ass:A2} 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{ass:A3} remains stated at the level used by the recovery
proofs. Appendix~\ref{app:primitive_wald} 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{ass:A4}, 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 \eqref{eq:nc_def_prim}
contains the factor $\sqrt T$, but the theorem is stated in the more general
rate form needed for high-dimensional testing.
Assumption \ref{ass:A5} 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 \eqref{eq:threshold_calibration} to hold; the boundary choice
$c_T=\sqrt{2\log p}$ need not be sufficient without such slack.
\section{Approximating-model results}
\label{sec:approx}
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 \emph{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.
\begin{theorem}
\label{thm:approx_stats} Suppose Assumptions~\ref{ass:A1}--\ref{ass:A3} hold.
Let $\{S^{(\ell)}\}_{\ell\geq0}$ denote any sequence of conditioning sets that
can be generated by the nonlinear BMT procedure, with $|S^{(\ell)}\setminus
S^{(0)}|\leq k_{\max}$ for all $\ell$. Then, as $T\rightarrow\infty$, the
following statements hold. \newline\noindent\textbf{(a) Exponential-tail
regime.} If Assumption A2-ET holds, then
\begin{equation}
\max_{\ell\leq k_{\max}}\max_{j\notin S^{(\ell)}}\bigl|W_{T,j}(S^{(\ell
)})-\mathrm{NC}_{j}^{\ast}(S^{(\ell)})\bigr|=O_{p}\!\left(\sqrt{\ell_T}\right).\label{eq:thm1_ET}
\end{equation}
\noindent\textbf{(b) Polynomial-moment regime.} If Assumption A2-PM holds,
then
\begin{equation}
\max_{\ell\leq k_{\max}}\max_{j\notin S^{(\ell)}}\bigl|W_{T,j}(S^{(\ell
)})-\mathrm{NC}_{j}^{\ast}(S^{(\ell)})\bigr|=O_{p}(a_{T,p}^{\mathrm{PM}
}).\label{eq:thm1_PM}
\end{equation}
\noindent\textbf{(c) Implication for population-null covariates.} If, in
addition, Assumption~\ref{ass:A5} holds, then under either tail regime
\begin{equation}
\Pr\!\left( \exists\,\ell\leq k_{\max},\ \exists\,j\notin S^{(\ell
)}:\mathrm{NC}_{j}^{\ast}(S^{(\ell)})=0\ \text{and}\ W_{T,j}(S^{(\ell)})\geq
c_{T}\right) \rightarrow0.\label{eq:thm1_null}
\end{equation}
\end{theorem}
\paragraph{Comments.}
\label{sec:thm1_discussion}
Theorem~\ref{thm:approx_stats} 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.
\bigskip
The results in Theorem~\ref{thm:approx_stats} 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 \emph{zero population conditional strength}
are not selected (uniformly along the path); (ii) the procedure stops once the
approximating model is \emph{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 \eqref{eq:nc_def_prim}). Define the
\emph{population-null} set at $S$ as
\begin{equation}
\mathcal{N}^{\ast}(S)=\{j\notin S:\ \mathrm{NC}_{j}^{\ast}(S)=0\},
\label{eq:nullset_def}
\end{equation}
and the maximal remaining population strength at stage $\ell$ as
\begin{equation}
M_{\ell}^{\ast}=\max_{j\notin S^{(\ell)}}\mathrm{NC}_{j}^{\ast}\bigl(S^{(\ell
)}\bigr). \label{eq:Mell_def}
\end{equation}
\begin{theorem}
\label{thm:approx_stop}
Suppose Assumptions~\ref{ass:A1}--\ref{ass:A3} and \ref{ass:A5} hold. Let
$\widehat L\leq k_{\max}$ denote the number of successful selections made
before the nonlinear BMT procedure stops, with $\widehat L=k_{\max}$ if the
prespecified stage cap is reached, and let $\widehat S=S^{(\widehat L)}$.
For notational convenience, extend the path after stopping by setting
$S^{(\ell)}=\widehat S$ for $\ell\geq\widehat L$. Then, as
$T\rightarrow\infty$:
\begin{enumerate}
\item[\textup{(i)}] (\emph{No selection of population-null covariates,
uniformly in stage})
\begin{equation}
\Pr\!\left(
\exists\,\ell\leq k_{\max},\ \exists\,j\in
\mathcal N^*\!\left(S^{(\ell)}\right)
\text{ such that BMT selects }j\text{ at stage }\ell
\right)\rightarrow0.
\label{eq:thm2_nonull}
\end{equation}
\item[\textup{(ii)}] (\emph{Stopping under approximating-model exhaustion}) If
there exists a possibly random stage $\ell^\dagger\le k_{\max}$ such that
\begin{equation}
M_{\ell^\dagger}^*=o_p(c_T),
\label{eq:exhaust_general}
\end{equation}
then
\begin{equation}
\Pr(\widehat L\le \ell^\dagger)\rightarrow1.
\label{eq:thm2_stop}
\end{equation}
\end{enumerate}
Part~\textup{(i)} may equivalently be written along the realised selected
path as
\begin{equation}
\Pr\!\left(
\exists\,\ell<\widehat L:
j^{(\ell)}\in\mathcal N^*\!\left(S^{(\ell)}\right)
\right)\rightarrow0.
\label{eq:thm2_path_nonull}
\end{equation}
\end{theorem}
\paragraph{Comments.}
\label{sec:thm2_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{thm:approx_stats}(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{thm:approx_stats} 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 \emph{sample} statistic falls below $c_{T}$
at (or before) $\ell^{\dagger}$ with probability tending to one.
\section{Exact recovery and oracle inference}
\label{sec:exact}
This section provides results relating to \emph{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 \emph{oracle properties}
of the post-selection MLE of the GLM parameters can be established.
\subsection{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}
\begin{assumption}
\label{ass:A41}
There exists a deterministic sequence $d_{1T}>0$ such that
\begin{equation}
\frac{\mathcal E_T}{d_{1T}}=o_p(1),
\label{eq:dom_k1_error}
\end{equation}
and, for every admissible conditioning set $S$ satisfying
$S^{(0)}\subseteq S\subseteq S^{(0)}\cup\{j_0\}$ and $j_0\notin S$,
\begin{equation}
\mathrm{NC}_{j_{0}}^{\ast}(S)
\ge
\max_{\substack{m\neq j_0\\m\notin S}}
\mathrm{NC}_{m}^{\ast}(S)+d_{1T},
\label{eq:dom_k1}
\end{equation}
and
\begin{equation}
\mathrm{NC}_{j_{0}}^{\ast}(S)
\ge c_T+d_{1T}.
\label{eq:dom_k1_threshold}
\end{equation}
\end{assumption}
Under the restriction in Assumption~\ref{ass:A41}, 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.
\begin{theorem}
\label{thm:k1_recovery} Suppose
Assumptions~\ref{ass:A1}--\ref{ass:A3}, \ref{ass:A41} and \ref{ass:A5} hold,
$k=|S_0|=1$, with $S_0=\{j_0\}$, and $k_{\max}\ge1$ in the nonlinear BMT
procedure. Then, as $T\rightarrow\infty$:
\medskip\noindent\newline\textbf{(a)} With probability
tending to one, the first-stage passing set is nonempty and the first selected
covariate is the signal:
\begin{equation}
\Pr\!\left(j^{(0)}=j_0\right)\rightarrow1.
\label{eq:thm3_first}
\end{equation}
\medskip\noindent\textbf{(b)} Consequently,
\begin{equation}
\Pr\!\left(j_0\in\widehat S\right)\rightarrow1.
\label{eq:thm3_include}
\end{equation}
\medskip\noindent\textbf{(c)} If, in addition,
\begin{equation}
\max_{\substack{m\neq j_0\\m\notin S^{(0)}}}
\mathrm{NC}_m^*\!\left(S^{(0)}\cup\{j_0\}\right)=o(c_T),
\label{eq:k1_exhaustion}
\end{equation}
then
\begin{equation}
\Pr\!\left(\widehat S=S^{(0)}\cup\{j_0\}\right)\rightarrow1.
\label{eq:thm3_exact}
\end{equation}
\end{theorem}
\paragraph{Proof outline.}
\label{sec:thm3_proofoutline}
We sketch the key steps in the proof here; a complete proof is given in
Appendix~\ref{app:thm3}.
\begin{enumerate}
\item \emph{Population separation with the correct margin}. By
\eqref{eq:dom_k1}, for any admissible $S$ not containing the signal,
$\mathrm{NC}_{j_0}^*(S)$ exceeds every competing population noncentrality by
at least $d_{1T}$.
\item \emph{Uniform sample approximation}. By definition of $\mathcal E_T$,
\[
W_{T,j_0}(S)\ge \mathrm{NC}_{j_0}^*(S)-\mathcal E_T,
\qquad
W_{T,m}(S)\le \mathrm{NC}_{m}^*(S)+\mathcal E_T.
\]
Since $\mathcal E_T/d_{1T}=o_p(1)$, the sample ordering inherits the
population ordering with probability tending to one.
\item \emph{Threshold crossing}. The additional condition
\eqref{eq:dom_k1_threshold} ensures that the signal is not only the largest
sample statistic but also exceeds the multiple-testing threshold.
\item \emph{Stopping after selecting the signal}. Condition
\eqref{eq:k1_exhaustion} ensures that, after the signal is included, all
remaining population noncentralities are asymptotically below the threshold.
The stopping theorem then rules out further selections.
\end{enumerate}
\subsection{Multiple signals}
\label{sec:theorem4}
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{ass:A4}.
The following theorem is the nonlinear analogue of the multiple-signal
separation and exact recovery theorem established in linear BMT analysis \citep[cf.][]{kapetanios2026}.
\begin{theorem}
\label{thm:fixedk_recovery} Suppose
Assumptions~\ref{ass:A1}--\ref{ass:A5} hold and $k_{\max}\geq k\geq1$ in the
nonlinear BMT procedure. Then, the following statements are true as
$T\rightarrow\infty$. \medskip\noindent\newline\textbf{(a) No false selections
before exhausting }$S_{0}$\textbf{:} With probability tending to one, at each
stage $\ell=0,1,\dots,k-1$,
\begin{equation}
j^{(\ell)}\in S_{0}\setminus S^{(\ell)}, \label{eq:thm4_true_each_stage}
\end{equation}
i.e.,\ the covariate selected at stage $\ell$ is a previously unselected
signal. In particular,
\begin{equation}
\Pr\!\left( S^{(k)}=S^{(0)}\cup S_{0}\right) \ \rightarrow\ 1,
\label{eq:thm4_after_k}
\end{equation}
where $S^{(k)}$ denotes the set after $k$ successful selections.
\medskip\noindent\newline\textbf{(b) Exact model recovery and stopping:}
Assume additionally that
\begin{equation}
\max_{\substack{m\notin S_{0}\\m\notin S^{(0)}}}
\mathrm{NC}_{m}^{\ast}\bigl(S^{(0)}\cup S_{0}\bigr)=o(c_T),
\label{eq:thm4_null_after_true}
\end{equation}
i.e., once all signals are included, the remaining population conditional
strengths are asymptotically below the testing threshold. Then,
\begin{equation}
\Pr\!\left( \widehat{S}=S^{(0)}\cup S_{0}\right) \rightarrow\ 1,
\label{eq:thm4_exact}
\end{equation}
i.e.,\ nonlinear BMT selects all and only the true signals and stops after $k$
stages (up to the always-in controls).
\end{theorem}
\paragraph{Comments.}
\label{sec:thm4_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 \eqref{eq:stg_dom},
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.}
\label{sec:thm4_proofoutline}
A complete proof is given in Appendix~\ref{app:thm4}; here, we summarise the
main argument.
\begin{enumerate}
\item \emph{Induction over stages}. Suppose that after $\ell<k$ stages the
selected non-always-in covariates are all signals, so that
$S^{(\ell)}\subseteq S^{(0)}\cup S_0$ and
$|S^{(\ell)}\cap S_0|=\ell$.
\item \emph{Population dominance at stage $\ell$}. Assumption~\ref{ass:A4}
gives a population gap of at least $d_T$ between every remaining signal and
every remaining non-signal, and it also ensures that the remaining signals
cross the threshold by at least $d_T$.
\item \emph{Uniform sample approximation}. Since $\mathcal E_T/d_T=o_p(1)$,
the sample Wald statistics preserve the population ordering with probability
tending to one. Hence a non-signal cannot maximise the stagewise statistic at
any stage before $S_0$ is exhausted, and the selected covariate must be a new
signal.
\item \emph{Stopping after exhaustion}. Once $S^{(k)}=S^{(0)}\cup S_0$,
condition \eqref{eq:thm4_null_after_true} and Theorem~\ref{thm:approx_stop}
imply that all remaining sample statistics are below $c_T$ with probability
tending to one. The algorithm therefore stops without false inclusions.
\end{enumerate}
\subsection{Oracle estimation and inference}
\label{sec:oracle}
We now establish the oracle properties of the post-nonlinear-BMT
estimator: conditional on exact model recovery, as established in
Theorem~\ref{thm:fixedk_recovery}, 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}
\begin{assumption}
\label{ass:A6}
The oracle criterion over $S^\dagger$ has a unique pseudo-true value
$\boldsymbol\beta_{0,S^\dagger}$. Its score
\[
\boldsymbol\psi_{0,t}
=
\frac{\partial}{\partial\boldsymbol\beta_{S^\dagger}}
l_t(\boldsymbol x_{S^\dagger,t}'\boldsymbol\beta_{0,S^\dagger},y_t)
\]
satisfies
\[
T^{-1/2}\sum_{t=1}^T\boldsymbol\psi_{0,t}
\Rightarrow N(\boldsymbol{0},\boldsymbol\Omega_0),
\]
where $\boldsymbol\Omega_0$ is the relevant long-run score covariance matrix ($\Rightarrow$ signifies convergence in distribution as
$T\rightarrow\infty$).
The oracle Hessian
\[
\boldsymbol J_0
=
-\mathbb E\left[
\frac{\partial^2}{\partial\boldsymbol\beta_{S^\dagger}\partial
\boldsymbol\beta_{S^\dagger}'}
l_t(\boldsymbol x_{S^\dagger,t}'\boldsymbol\beta_{0,S^\dagger},y_t)
\right]
\]
is nonsingular, and the oracle estimator $\widehat{\boldsymbol\beta}_{\rm or}$
satisfies the asymptotic linear expansion
\[
\sqrt T(\widehat{\boldsymbol\beta}_{\rm or}-\boldsymbol\beta_{0,S^\dagger})
=
\boldsymbol J_0^{-1}T^{-1/2}\sum_{t=1}^T\boldsymbol\psi_{0,t}+o_p(1).
\]
\end{assumption}
We have the following result.
\begin{theorem}
\label{thm:oracle} Suppose
Assumption~\ref{ass:A6} holds and
\(\Pr(\widehat S=S^\dagger)\rightarrow1\) as $T\rightarrow\infty$.
Then,
\begin{equation}
\widehat{\boldsymbol\beta}_{{\rm post},S^\dagger}-\boldsymbol\beta_{0,S^\dagger}
=o_p(1),
\label{eq:oracle_consistency}
\end{equation}
and
\begin{equation}
\sqrt{T}\bigl(\widehat{\boldsymbol\beta}_{{\rm post},S^\dagger}-
\boldsymbol\beta_{0,S^\dagger}\bigr)
\Rightarrow
N\!\left(\boldsymbol{0},\boldsymbol J_0^{-1}\boldsymbol\Omega_0\boldsymbol J_0^{-1}\right).
\label{eq:oracle_an}
\end{equation}
Consequently, the post-BMT estimator has the same first-order limiting
distribution as the infeasible oracle estimator that knows $S^\dagger$ in advance.
\end{theorem}
\paragraph{Proof outline.}
\label{sec:thm5_proofsketch}
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{ass:A6} 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{app:thm5}.
\subsection{Logit and probit models}
\label{sec:logit_probit}
We verify that the conditions of Theorems~\ref{thm:approx_stats}
--\ref{thm:oracle} 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.
\begin{lemma}
\label{lem:logit_probit} Let $G=\Lambda$ or $G=\Phi$.
Suppose Assumptions~\ref{ass:A1}--\ref{ass:A2} and conditions A3.2--A3.3
hold. Suppose also that the pseudo-true parameters of all admissible stagewise
submodels lie in a common compact set, that their associated indices lie in a
fixed compact interval, and that the additional parameter-set, growth and
standard-error conditions of Theorem~\ref{thm:primitive_wald} hold. Then A3.4
holds with the rates stated in that theorem, and hence
Assumption~\ref{ass:A3} holds. If, in addition, the postulated conditional
mean is correctly specified and
$S\supseteq S^{\dagger}=S^{(0)}\cup S_0$, then
$\theta_m^*(S)=0$ and hence $\mathrm{NC}_m^*(S)=0$ for every candidate
$m\notin S$.
\end{lemma}
\paragraph{Comment.}
For logit, the scalar score and curvature weights are bounded, so the
derivative-envelope conditions in Assumption~\ref{ass:A2} 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{app:primitive_wald}
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. \medskip
Together with Theorems~\ref{thm:approx_stats}--\ref{thm:oracle}, the
results in Lemma~\ref{lem:logit_probit} 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.
\section{Primitive conditions for stagewise dominance}
\label{sec:primitive}
In this section, we provide primitive sufficient conditions for the stagewise
dominance conditions (\ref{eq:stg_dom}) and (\ref{eq:dom_k1}) required for the
exact-recovery results of Theorems~\ref{thm:k1_recovery} and
\ref{thm:fixedk_recovery}. We consider two distinct routes relying on
assumptions that serve different empirical and conceptual purposes. \emph{The
local (small-index) regime} (Assumptions \ref{ass:L1}--\ref{ass:L3} 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 \emph{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
\emph{conditional proxy weakness} through Assumption~\ref{ass:GM2}.
Assumption~\ref{ass:G1} gives a primitive moment-based sufficient condition,
whereas Assumption~\ref{ass:G2} is a partial-projection diagnostic that is
useful when additional model-specific restrictions link it to
Assumption~\ref{ass:G1}.
The local conditions and the global-regime condition are independent
alternatives. Either Assumptions~\ref{ass:L1}--\ref{ass:L3}, or
Assumption~\ref{ass:GM2} verified through an appropriate primitive route, is
sufficient for stagewise dominance. Theorems~\ref{thm:k1_recovery} and
\ref{thm:fixedk_recovery} 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 \eqref{eq:nc_def_prim}, we note that a dominance
condition such as Assumption \ref{ass:A4} 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.
\subsection{Local regime linearisation}
\label{sec:local_regime}
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}
\begin{assumption}
\label{ass:L1} There are a constant $\delta_0>0$ and a deterministic sequence
$0<\delta_T\leq\delta_0$ such that
\[
\sup_{1\leq t\leq T}|\boldsymbol{x}_{t}^{\prime}\boldsymbol{\beta}_{0}|
\leq\delta_T
\]
with probability tending to one. The sequence $\delta_T$ may be fixed at a
sufficiently small value or may converge to zero in a local triangular array.
\end{assumption}
\begin{assumption}
\label{ass:L2} On $[-\delta_0,\delta_0]$, $G$
satisfies $0<\underline{g}\leq G^{\prime}(z)\leq\bar{g}<\infty$ and
$|G^{\prime\prime}(z)|\leq M<\infty$ for some constants $\underline{g}$,
$\bar{g}$ and $M$.
\end{assumption}
\begin{assumption}
[Scaled local partial-correlation gap]\label{ass:L3} For every admissible
conditioning set $S^{(0)}\subseteq S\subseteq S^{(0)}\cup S_{0}$ with
$|S\cap S_{0}|\leq k-1$, let $r_{j,t}(S)$ denote the residual from the
least-squares projection of $x_{j,t}$ on $\boldsymbol{x}_{S,t}$ and put
\(M_j(S)=\left|\mathbb E[r_{j,t}(S)\boldsymbol x_t'\boldsymbol\beta_0]\right|.\)
There exist constants $0<C_L\leq C_U<\infty$ and a deterministic sequence
$\underline{\kappa}_T>0$ such that
\begin{equation}
C_L\min_{j\in S_0\setminus S}M_j(S)
\geq
C_U\max_{\substack{m\notin S_0\\m\notin S}}M_m(S)
+\underline{\kappa}_T.
\label{eq:local_gap}
\end{equation}
\end{assumption}
The next lemma provides the link between the local regime and stagewise dominance.
\begin{lemma}
[Local regime implies dominance]\label{lem:local_implies_A4} Under Assumptions
A3.3 and \ref{ass:L1}--\ref{ass:L3}, suppose that the uniform local expansion
takes the form
\begin{equation}
\sup_{\substack{S^{(0)}\subseteq S\subseteq S^{(0)}\cup S_0\\j\notin S}}
\left|
\frac{|\theta_j^*(S)|}{\{V_{\theta\theta,j}^*(S)\}^{1/2}}
-C_j(S)M_j(S)
\right|
\leq r_T,
\qquad
C_L\leq C_j(S)\leq C_U,
\label{eq:local_uniform_expansion}
\end{equation}
where $r_T$ is deterministic and
$r_T=o(\underline\kappa_T)$. If
$c_T+\mathcal E_T=o_p(\sqrt T\,\underline\kappa_T)$, then
Assumption~\ref{ass:A4} holds with
$d_T=\sqrt T\,\underline\kappa_T/4$.
\end{lemma}
\paragraph{Comment.}
The local regime reduces dominance to a \emph{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{ass:L1}
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.
\subsection{Global regime proxy weakness}
\label{sec:global_regime}
Next, we formalise the global regime condition (Assumption \ref{ass:GM2}),
then give a moment-based sufficient condition and a complementary
partial-projection diagnostic.
\setcounter{assumption}{0}
\begin{assumption}
\label{ass:GM2} There exists $\rho\in(0,1)$
such that, for any admissible $S^{(0)}\subseteq S\subseteq S^{(0)}\cup S_0$
with $|S\cap S_0|\le k-1$,
\begin{equation}
\max_{\substack{m\notin S_0\\m\notin S}}
\frac{|\theta_m^*(S)|}{\{V_{\theta\theta,m}^*(S)\}^{1/2}}
\le
\rho
\min_{j\in S_0\setminus S}
\frac{|\theta_j^*(S)|}{\{V_{\theta\theta,j}^*(S)\}^{1/2}}.
\label{eq:GM2_theta}
\end{equation}
In addition, the remaining signals are strong enough in Wald-noncentrality
scale:
\begin{equation}
b_T
:=
\inf_{S}
\min_{j\in S_0\setminus S}\mathrm{NC}_j^*(S)
\quad\text{satisfies}\quad
\frac{c_T+\mathcal E_T}{b_T}=o_p(1),
\label{eq:GM2_signal_strength}
\end{equation}
where the infimum is over the same admissible conditioning sets.
\end{assumption}
Stagewise dominance follows from (\ref{eq:GM2_theta}) by virtue of the next result.
\begin{lemma}
[Global regime implies dominance]\label{lem:GM2_implies_A4} Under
Assumptions~A3.3 and \ref{ass:GM2}, the stagewise dominance condition of
Assumption~\ref{ass:A4} is satisfied with
$d_T=(1-\rho)b_T/2$.
\end{lemma}
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}
\begin{assumption}
\label{ass:G1} There exists $\rho
\in(0,1)$ such that, for all admissible $S$,
\begin{equation}
\max_{\substack{m\notin S_{0}\\m\notin S}}\Bigl|\mathbb{E}\bigl[\widetilde{x}_{m,t}(S)\,\xi
_{t}(S)\bigr]\Bigr|\ \leq\ \rho\,\min_{j\in S_{0}\setminus S}\Bigl|\mathbb{E}
\bigl[\widetilde{x}_{j,t}(S)\,\xi_{t}(S)\bigr]\Bigr|. \label{eq:G1}
\end{equation}
\end{assumption}
The following lemma gives conditions under which (\ref{eq:G1}) is sufficient
for (\ref{eq:GM2_theta}) to hold.
\begin{lemma}
\label{lem:G1_implies_GM2} Suppose
Assumption~\ref{ass:A3} holds and there exist constants
$0<C_L\le C_U<\infty$ such that, uniformly over admissible $S$ and $j\notin S$,
\[
C_L\left|\mathbb E[\widetilde x_{j,t}(S)\xi_t(S)]\right|
\le
\frac{|\theta_j^*(S)|}{\{V_{\theta\theta,j}^*(S)\}^{1/2}}
\le
C_U\left|\mathbb E[\widetilde x_{j,t}(S)\xi_t(S)]\right|.
\]
If Assumption~\ref{ass:G1} holds with a constant $\rho$ satisfying
$\rho C_U/C_L<1$, and if the signal-strength condition
\eqref{eq:GM2_signal_strength} holds, then Assumption~\ref{ass:GM2} holds.
\end{lemma}
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)$.
\begin{assumption}
\label{ass:G2} The curvature weights are positive for the designs to which
this condition is applied. There exists $\rho\in(0,1)$ such that, for every
admissible $S$,
\begin{equation}
\sup_{\substack{m\notin S_{0}\\m\notin S}}
\frac{\mathbb E\!\left[w_t(S)
\{\mathcal P^w_{R\mid S}\widetilde x_{m,t}(S)\}^{2}\right]}
{\mathbb E\!\left[w_t(S)\widetilde x_{m,t}^{2}(S)\right]}
\leq\rho,
\qquad R=S_0\setminus S.
\label{eq:G2}
\end{equation}
Thus the weighted partial $R^2$ of an inactive candidate with the remaining
signals, after conditioning on the current model, is bounded away from one.
\end{assumption}
The bound in \eqref{eq:G2} is useful only when it can be shown to deliver the
moment attenuation in \eqref{eq:G1}; by itself it does not control the sign or
size of the score moments.
\begin{remark}
\label{rem:G2_diagnostic}
Assumption~\ref{ass:G2} is therefore a design diagnostic rather than an
independent sufficient condition for stagewise dominance. When additional
model-specific restrictions imply \eqref{eq:G1} with a constant satisfying
the restriction in Lemma~\ref{lem:G1_implies_GM2}, that lemma, together with
\eqref{eq:GM2_signal_strength}, yields Assumption~\ref{ass:GM2}. No such
implication is claimed from Assumption~\ref{ass:G2} alone.
\end{remark}
\subsection{Logit and probit links}
\label{sec:link_verification_GM2}
Finally, we verify the key analytic condition used in
Lemma~\ref{lem:G1_implies_GM2} for logit and probit links.
\begin{lemma}
\label{lem:logit_probit_mapping} For the logit and probit links, suppose
Assumption~\ref{ass:A3} holds and the relevant indices lie in a compact set.
Suppose further that the profiled population score for the candidate
coefficient has Schur-complement curvature uniformly bounded away from zero
and infinity on the line segment between $0$ and $\theta_j^*(S)$, uniformly
over admissible $S$ and $j$. Then, there exist constants
$0<C_L\le C_U<\infty$ such that, uniformly over admissible $S$ and $j$,
\[
C_L\left|\mathbb{E}[\widetilde{x}_{j,t}(S)\xi_t(S)]\right|
\le
\frac{|\theta_j^*(S)|}{\{V_{\theta\theta,j}^*(S)\}^{1/2}}
\le
C_U\left|\mathbb{E}[\widetilde{x}_{j,t}(S)\xi_t(S)]\right|.
\]
\end{lemma}
In logit and probit models, Assumption~\ref{ass:G1}, or Assumption~\ref{ass:G2}
when it is verified to imply the moment attenuation in \eqref{eq:G1}, implies
Assumption \ref{ass:GM2}, which in turn implies Assumption \ref{ass:A4}
(Lemma~\ref{lem:G1_implies_GM2}; see also
Remark~\ref{rem:G2_diagnostic}). Alternatively,
the local small-index conditions
(\ref{ass:L1}--\ref{ass:L3}) yield dominance directly
(Lemma~\ref{lem:local_implies_A4}). 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{app:binary_lasso} compares the stagewise population ordering used
by BMT with the generalised irrepresentability condition for binary LASSO. The
comparison is not one-sided. Theorem~\ref{thm:app_lasso_bmt_success} gives a
Gaussian proxy design in which BMT recovers the true support although binary
LASSO violates generalised irrepresentability, whereas
Theorem~\ref{thm:app_lasso_lasso_success} 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{thm:app_lasso_gi_bmt} 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{ass:A4} additionally requires the active-score balance
condition stated in Appendix~\ref{app:binary_lasso}.
\section{Monte Carlo simulations}
\label{sec:Monte_Carlo}
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.}
\subsection{Experimental design}
\label{sec:Monte_Carlo_design}
The DGP is the binary-response GLM with linear single-index structure given
by
\begin{equation}
\mathbb{E}(y_{t}|\boldsymbol{x}_{t})=\Lambda(\alpha+\boldsymbol{x}_{t}
^{\prime}\boldsymbol{\beta}_{0}),\quad t=1,2,\ldots,T,\label{eq_dgp1}
\end{equation}
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 \eqref{eq_dgp1} --- 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
\begin{equation}
x_{j,t}=\frac{\varepsilon_{j,t}+\nu_{g}g_{t}+\nu_{f}f_{t}}{(1+\nu_{g}^{2}
+\nu_{f}^{2})^{1/2}},\quad j=1,\ldots,2k, \label{eq_dgp2}
\end{equation}
while the remaining noise variables $x_{2k+1,t},\ldots,x_{p,t}$ are obtained
as
\begin{equation}
x_{j,t}=\frac{\varepsilon_{j-1,t}+\varepsilon_{j,t}+\nu_{f}f_{t}}{(2+\nu
_{f}^{2})^{1/2}},\quad j=2k+1,\ldots,p. \label{eq_dgp3}
\end{equation}
In \eqref{eq_dgp2}--\eqref{eq_dgp3}, $\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,
\begin{align}
\varepsilon_{j,t} & =\rho\varepsilon_{j,t-1}+(1-\rho^{2})^{1/2}
e_{j,t}^{\varepsilon},\label{eq_dgp41}\\
f_{t} & =\rho f_{t-1}+(1-\rho^{2})^{1/2}e_{t}^{f},\label{eq_dgp42}\\
g_{t} & =\rho g_{t-1}+(1-\rho^{2})^{1/2}e_{t}^{g}, \label{eq_dgp43}
\end{align}
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{sec:algo}, 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.
\subsection{Performance metrics}
\label{sec:diagnostics}
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
\begin{align}
TP & =\sum_{j=1}^{p}\mathbb{I}(\widehat{D}_{j}=1,\,\beta_{0,j}\neq0),\quad
FP=\sum_{j=1}^{p}\mathbb{I}(\widehat{D}_{j}=1,\,\beta_{0,j}=0),\label{eq_tp}\\
TN & =\sum_{j=1}^{p}\mathbb{I}(\widehat{D}_{j}=0,\,\beta_{0,j}=0),\quad
FN=\sum_{j=1}^{p}\mathbb{I}(\widehat{D}_{j}=0,\,\beta_{0,j}\neq0),
\label{eq_tn}
\end{align}
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 \eqref{eq_dgp1}, 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 \eqref{eq_tp}--\eqref{eq_tn}; 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{\citet{Baldi2020} and \citet{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 \citep{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.
\subsection{Simulation results}
\label{sec:MCresults}
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{tab_summary_mcc_alpha0}$-$
\ref{tab_summary_rmse_alpha0} 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).
\begin{table}[h]
\caption{Summary Statistics for Matthews Correlation Coefficient ($MCC$):
$\alpha=0$}
\label{tab_summary_mcc_alpha0}
\centering
\begin{tabular}
[c]{ccccccc}\hline\hline
Signals & Method & MED & IQR & MIN & MAX & RANK\\
\hline
$k=1$ & BMT & 0.99 & 0.00 & 0.99 & 1.00 & 1.00 \\
& OCMT & 0.53 & 0.65 & 0.02 & 0.99 & 2.30 \\
& LASSO & 0.33 & 0.05 & 0.29 & 0.39 & 2.70 \\
$k=4$ & BMT & 0.99 & 0.08 & 0.84 & 1.00 & 1.00 \\
& OCMT & 0.76 & 0.40 & 0.13 & 1.00 & 2.10 \\
& LASSO & 0.35 & 0.07 & 0.29 & 0.42 & 2.90 \\
\hline\hline
\end{tabular}
\begin{minipage}{11cm}\vspace{0.25cm}
\footnotesize{Notes: $MCC$ ranges between $-1$ (perfect misclassification) and $1$ (perfect classification), with 0 indicating random classification.}
\end{minipage}\end{table}
It is clear from the results in Table~\ref{tab_summary_mcc_alpha0} 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$.
\begin{table}[h]
\caption{Summary Statistics for $F1$ Score: $\alpha=0$}
\label{tab_summary_f1_alpha0}
\centering
\begin{tabular}
[c]{ccccccc}\hline\hline
Signals & Method & MED & IQR & MIN & MAX & RANK\\
\hline
$k=1$ & BMT & 0.99 & 0.00 & 0.99 & 1.00 & 1.00 \\
& OCMT & 0.45 & 0.70 & 0.01 & 0.99 & 2.30 \\
& LASSO & 0.22 & 0.05 & 0.17 & 0.29 & 2.70 \\
$k=4$ & BMT & 0.99 & 0.09 & 0.83 & 1.00 & 1.00 \\
& OCMT & 0.74 & 0.46 & 0.06 & 1.00 & 2.10 \\
& LASSO & 0.25 & 0.09 & 0.18 & 0.35 & 2.90 \\
\hline\hline
\end{tabular}
\begin{minipage}{11cm}\vspace{0.25cm}
\footnotesize{Notes: $F1$ ranges between 0 (no signal selected) and 1 (all signals selected).}
\end{minipage}\end{table}
The superior performance of BMT is also evident in the summary statistics for
the $F1$ score reported in Table~\ref{tab_summary_f1_alpha0}. 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.
\begin{table}[h]
\caption{Summary Statistics for Model Size: $\alpha=0$}
\label{tab_summary_khat_alpha0}
\centering
\begin{tabular}
[c]{ccccccc}\hline\hline
Signals & Method & MED & IQR & MIN & MAX & RANK\\\hline
$k=1$ & BMT & 1.02 & 0.01 & 1.01 & 1.03 & 1.00 \\
& OCMT & 6.30 & 52.93 & 1.02 & 325.76 & 2.40 \\
& LASSO & 12.87 & 4.00 & 9.36 & 17.48 & 2.60 \\
$k=4$ & BMT & 3.97 & 0.58 & 2.94 & 4.02 & 1.10 \\
& OCMT & 7.06 & 15.30 & 2.94 & 160.72 & 2.05 \\
& LASSO & 31.11 & 13.06 & 20.45 & 46.13 & 2.85 \\
\hline\hline
\end{tabular}
\begin{minipage}{11cm}\vspace{0.25cm}
\footnotesize{Notes: Model size is the number of selected covariates.}
\end{minipage}
\end{table}
Another attractive feature of BMT is that selection accuracy is not achieved
by sacrificing model parsimony. This can be seen in
Table~\ref{tab_summary_khat_alpha0}, 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.
\begin{table}[h]
\caption{Summary Statistics for Bayesian Information Criterion ($BIC$):
$\alpha=0$}
\label{tab_summary_bic_alpha0}
\centering
\begin{tabular}
[c]{ccccccc}\hline\hline
Signals & Method & MED & IQR & MIN & MAX & RANK\\\hline
$k=1$ & BMT & 215.94 & 140.00 & 145.90 & 287.93 & 1.00 \\
& OCMT & 295.24 & 441.51 & 146.08 & 3524.60 & 2.40 \\
& LASSO & 328.57 & 193.81 & 217.01 & 459.13 & 2.60 \\
$k=4$ & BMT & 146.91 & 75.18 & 101.77 & 216.45 & 1.00 \\
& OCMT & 213.18 & 132.92 & 115.70 & 1852.99 & 2.10 \\
& LASSO & 417.32 & 248.08 & 247.74 & 678.28 & 2.90 \\
\hline\hline
\end{tabular}
\end{table}
Turning to the complexity-penalised measure of fit,
Table~\ref{tab_summary_bic_alpha0} 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{tab_summary_rmse_alpha0}. 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$.
\begin{table}[h]
\caption{Summary Statistics for Root Mean Squared Error ($RMSE$): $\alpha=0$}
\label{tab_summary_rmse_alpha0}
\centering
\begin{tabular}
[c]{ccccccc}\hline\hline
Signals & Method & MED & IQR & MIN & MAX & RANK\\\hline
$k=1$ & BMT & 0.33 & 0.12 & 0.27 & 0.41 & 1.15 \\
& OCMT & 1.63 & 74.23 & 0.26 & 1553.86 & 2.40 \\
& LASSO & 0.77 & 0.18 & 0.61 & 0.99 & 2.45 \\
$k=4$ & BMT & 0.93 & 0.30 & 0.69 & 1.20 & 1.15 \\
& OCMT & 2.79 & 13.80 & 0.68 & 166.62 & 2.45 \\
& LASSO & 1.90 & 0.30 & 1.66 & 2.24 & 2.40 \\
\hline\hline
\end{tabular}
\end{table}
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{tab_summary_time_alpha0}, and is faster than OCMT in the very sparse
single-signal design.
\begin{table}[h]
\caption{Computation Time (in hours): $\alpha=0$}
\label{tab_summary_time_alpha0}
\centering
\begin{tabular}
[c]{ccccc}\hline\hline
\multicolumn{1}{c}{} & \multicolumn{2}{c}{$T=150$} &
\multicolumn{2}{c}{$T=300$}\\\cline{2-5}
Method & $k=1$ & $k=4$ & $k=1$ & $k=4$\\\hline
BMT & 0.02 & 0.07 & 0.09 & 0.28\\
OCMT & 0.03 & 0.04 & 0.28 & 0.27\\
LASSO & 1.62 & 0.65 & 3.73 & 1.67\\\hline\hline
\end{tabular}
\end{table}
\section{Empirical illustration}
\label{sec:Example}
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\% \citep{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
\begin{equation}
\mathbb{E}(y_{t+12}|\boldsymbol{x}_{t})=\Lambda(\alpha+\boldsymbol{x}
_{t}^{\prime}\boldsymbol{\beta}),\quad t\geq1, \label{eq_ppm}
\end{equation}
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 \cite{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 \cite{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{sec:Monte_Carlo}. 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$) \citep{Bamber1975} and the square root of the
quadratic probability score ($SQPS$). The components of the $MCC$ are obtained
in this case as
\begin{align*}
TP & =\sum_{t\in\mathcal{T}}\mathbb{I}(\widehat{y}_{t+12}>\tau
,y_{t+12}=1),\quad FP=\sum_{t\in\mathcal{T}}\mathbb{I}(\widehat{y}_{t+12}
>\tau,y_{t+12}=0),\\
TN & =\sum_{t\in\mathcal{T}}\mathbb{I}(\widehat{y}_{t+12}\leq\tau
,y_{t+12}=0),\quad FN=\sum_{t\in\mathcal{T}}\mathbb{I}(\widehat{y}_{t+12}
\leq\tau,\,y_{t+12}=1),
\end{align*}
where $\widehat{y}_{t+12}\in(0,1)$ is the fitted (predicted) conditional mean
of $y_{t+12}$ from \eqref{eq_ppm}, $\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$.
\subsection{In-sample results}
In-sample results for the three covariate-selection methods are summarised in
Table~\ref{tab_fit}.\footnote{The BMT, OCMT and LASSO procedures are
implemented as in the simulations (see Section~\ref{sec:Monte_Carlo_design}). 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
\cite{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})$.
\begin{table}[h]
\caption{In-Sample Performance}
\centering
\begin{tabular}
[c]{c|ccc}\hline\hline
& BMT & OCMT & LASSO\\
\hline
$MCC$ & 0.803 & 0.939 & 0.994 \\
$AUROC$ & 0.974 & 0.996 & 1.000 \\
$SQPS$ & 0.255 & 0.159 & 0.070 \\
$BIC$ & 318.0 & 470.6 & 575.2 \\
$\widehat{k}$ & 5 & 33 & 49\\
\hline\hline
\end{tabular}
\label{tab_fit}
\end{table}
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
\subsection{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{tab_oos}. 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$.
\begin{table}[h]
\caption{Out-of-Sample Performance}
\label{tab_oos}
\centering
\begin{tabular}
[c]{c|ccc|ccc}\hline\hline
& \multicolumn{3}{c|}{Scheme 1: Fixed Predictors,} & \multicolumn{3}{c}{Scheme
2: Updated Predictors,}\\
& \multicolumn{3}{c|}{Updated Coefficients} & \multicolumn{3}{c}{Updated
Coefficients}\\\hline
& BMT & OCMT & LASSO & BMT & OCMT & LASSO\\\hline
$MCC$ & 0.524 & 0.433 & 0.059 & 0.704 & 0.479 & 0.347 \\
$AUROC$ & 0.863 & 0.640 & 0.281 & 0.917 & 0.695 & 0.616 \\
$SQPS$ & 0.419 & 0.497 & 0.631 & 0.333 & 0.474 & 0.549 \\
\hline\hline
\end{tabular}
\end{table}
\begin{figure}[h]
\caption{Selected Predictors Across Economic Groups Under Forecasting Scheme 1
(Fixed Predictors and Updated Coefficients)}
\label{fig_heat1}
\centering \includegraphics[width=\textwidth]{fig_heat1.png}\end{figure}
\begin{figure}[h]
\caption{Selected Predictors Across Economic Groups Under Forecasting Scheme 2
(Updated Predictors and Coefficients)}
\label{fig_heat2}
\centering
\includegraphics[width=\textwidth]{fig_heat2.png}\end{figure}
Figures~\ref{fig_heat1} and \ref{fig_heat2} 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{fig_heat2}, 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{fig_heat1} 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 \citep[e.g.,][]{Stock1999,Faust2013}. They also align with popular
forecasting models for inflation in which factors such as financial variables
\citep{Forni2003} and commodity prices \citep{Chen2014} have been found to
have significant predictive power.\newline
\section{Conclusion}
\label{sec: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 \citet{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.
\bigskip
\bibliographystyle{plainnat}
\bibliography{literature}
\clearpage