1.1in1.1in0.55in0.8in

# High-Dimensional Knockoffs Inference for Time Series Data††thanks:  Chien-
Ming Chi is Assistant Research Fellow, Institute of Statistical Science,
Academia Sinica, Taiwan (E-mail: xbbchi@stat.sinica.edu.tw). Yingying Fan is
Centennial Chair in Business Administration and Professor, Data Sciences and
Operations Department, Marshall School of Business, University of Southern
California, Los Angeles, CA 90089 (E-mail: fanyingy@marshall.usc.edu). Ching-
Kang Ing is Distinguished Professor, Institute of Statistics, National Tsing
Hua University, Hsinchu 30013, Taiwan (E-mail: cking@stat.nthu.edu.tw). Jinchi
Lv is Kenneth King Stonier Chair in Business Administration and Professor,
Data Sciences and Operations Department, Marshall School of Business,
University of Southern California, Los Angeles, CA 90089 (E-mail:
jinchilv@marshall.usc.edu). This work was supported by NSF Grants DMS-1953356
and EF-2125142, and by Grants 111-2118-M-001-012-MY2 and 109-2118-M-007-MY3
from the National Science and Technology Council, Taiwan.

Chien-Ming Chi1, Yingying Fan2, Ching-Kang Ing3 and Jinchi Lv2
Academia Sinica1, University of Southern California2
and National Tsing Hua University3

(May 18, 2023)

###### Abstract

The model-X knockoffs framework provides a flexible tool for achieving finite-
sample false discovery rate (FDR) control in variable selection in arbitrary
dimensions without assuming any dependence structure of the response on
covariates. It also completely bypasses the use of conventional p-values,
making it especially appealing in high-dimensional nonlinear models. Existing
works have focused on the setting of independent and identically distributed
observations. Yet time series data is prevalent in practical applications in
various fields such as economics and social sciences. This motivates the study
of model-X knockoffs inference for time series data. In this paper, we make
some initial attempt to establish the theoretical and methodological
foundation for the model-X knockoffs inference for time series data. We
suggest the method of time series knockoffs inference (TSKI) by exploiting the
ideas of subsampling and e-values to address the difficulty caused by the
serial dependence. We also generalize the robust knockoffs inference in [5] to
the time series setting and relax the assumption of known covariate
distribution required by model-X knockoffs, because such an assumption is
overly stringent for time series data. We establish sufficient conditions
under which TSKI achieves the asymptotic FDR control. Our technical analysis
reveals the effects of serial dependence and unknown covariate distribution on
the FDR control. We conduct power analysis of TSKI using the Lasso coefficient
difference knockoff statistic under linear time series models. The finite-
sample performance of TSKI is illustrated with several simulation examples and
an economic inflation study.

Running title: TSKI

Key words: Model-X knockoffs; Time series; High dimensionality; FDR control;
Power analysis; Sparsity; Interpretable forecasting; E-values

##  1 Introduction

Selecting important covariates contributing to a response variable from a set
of candidate covariates is an important problem arising in various scientific
fields including social and medical sciences. When the number of candidate
variables is large, controlling the false discovery rate (FDR) of the variable
selection procedure is an effective and popular way to control the error rate.
Most existing works achieve this goal by building procedures based on the
p-values constructed for assessing the importance of individual variables;
see, e.g., the seminal works of [7, 8]. Yet, the high dimensionality in
covariates and the possibly complicated nonlinear model structure that are
frequently encountered in the big data era make many conventional ways of
p-value calculations inapplicable or even completely fail [20]. To overcome
such difficulties, the framework of knockoffs inference was proposed in [4,
11] to achieve the goal of exact FDR control in variable selection in finite
samples completely bypassing the use of the conventional p-values in high-
dimensional regression models. It allows arbitrary dependence structure of the
response on covariates and arbitrary dimensionality of covariates at the cost
of assuming the known joint covariate distribution. While various efforts have
been made to ease the implementation and relax the technical assumptions on
knockoffs inference [27, 19, 22, 3, 41], such endeavors are almost always in
the setting of independent observations without serial dependence. As a
result, it is unclear whether these existing works are applicable to the time
series setting where serial dependence is present. Our work tries to fill the
gap by developing a theoretical and methodological foundation of model-X
knockoffs inference [11] in the time series setting.

We now briefly review the model-X knockoffs to better lay out the foundation.
Let $\boldsymbol{y}$ be an $n$-dimensional response vector and
$\boldsymbol{X}$ an $n\times p$ matrix of covariates across $n$ independent
observations. The knockoffs inference aims at selecting relevant covariates
according to Definition 2 (see Section 2.1) while keeping the FDR under
control. To this end, the model-X knockoffs framework [11] constructs an
$n\times p$ matrix
$\widetilde{\boldsymbol{X}}=(\widetilde{\boldsymbol{x}}_{1},\cdots,\widetilde{\boldsymbol{x}}_{p})$
with $\widetilde{\boldsymbol{x}}_{j}\in\mathbb{R}^{n}$, $j=1,\cdots,p$, using
the joint distribution of $\boldsymbol{X}$ (assumed to be known) such that

$\widetilde{\boldsymbol{X}}\perp\\!\\!\\!\perp\boldsymbol{y}|\boldsymbol{X}\
\text{ and }\
(\boldsymbol{X},\widetilde{\boldsymbol{X}})_{\textnormal{swap}(S)}\overset{d}{=}(\boldsymbol{X},\widetilde{\boldsymbol{X}})$
(1)

for each $S\subset\\{1,\cdots,p\\}$, where swap($S$) denotes the swapping
operation meaning that for each $j\in S$, columns $j$ and $j+p$ are swapped,
and $\overset{d}{=}$ stands for equal in distribution. Since
$\widetilde{\boldsymbol{X}}$ is constructed conditionally independent of
$\boldsymbol{y}$, it holds that
$(\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{y})_{\textnormal{swap}(S)}\overset{d}{=}(\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{y})$
for each $S$ containing only irrelevant covariates. At a high level, the
model-X knockoffs inference creates “fake” covariates
$\widetilde{\boldsymbol{X}}$ which perfectly mimic the “behavior” of the
original covariates. By using these “fake” covariates as control, the
importance of original covariates can be inferred. The fake covariates
$\widetilde{\boldsymbol{x}}_{j}$’s are referred to as model-X knockoff
variables (knockoffs for short) in [11].

Two important assumptions in [11] are that 1) the observations in
$(\boldsymbol{X},\boldsymbol{y})\in\mathbb{R}^{n\times(p+1)}$ across the rows
are independent and identically distributed (i.i.d.) and 2) the joint
distribution of $\boldsymbol{X}$ is known for generating the knockoff matrix
$\widetilde{\boldsymbol{X}}$. The first assumption ensures that the knockoff
variables can be constructed in a rowwise fashion independent of other rows.
The second assumption makes it possible to construct ideal knockoffs
satisfying the two critical conditions in (1). However, for time series data
with serial dependence, it is unclear whether the rowwise knockoff variable
construction still gives us nearly valid knockoffs inference and what the
precise effect of serial dependence on FDR control is. Also, in time series
applications when the covariates are lagged variables with stationarity
assumption, assuming known covariate distribution directly gives us the set of
important variables, invalidating the problem of variable selection. In this
paper, we will address theses challenges for time series applications by
relaxing the two aforementioned assumptions.

Regarding the assumption of i.i.d. rows, we overcome the difficulties caused
by serial dependence with the idea of subsampling. Meanwhile, to address the
challenges related to known covariate distribution, we generalize the robust
knockoffs inference in [5] to the time series setting, where the robust
knockoffs inference [5] is a recent innovative work that allows for
approximate covariate distribution, as opposed to known covariate
distribution, in non-time series applications. We will name the knockoff
variables generated using approximate covariate distribution as approximate
knockoffs to ease the presentation. Specifically, our analyses reveal that
with approximate knockoffs generated in a rowwise fashion ignoring the serial
dependence, the FDR inflation has an upper bound depending on the
Kullback–Leibler (KL) divergence between the distributions of data matrices
corresponding to approximate and exact model-X knockoff variables (see Theorem
2). There is generally no guarantee that such KL divergence could
asymptotically vanish in the existence of serial dependence, indicating that
the corresponding FDR is generally uncontrolled. To address such difficulty,
we apply the subsampling technique to obtain asymptotically independent rows
[40] under the assumption that the time series are $\beta$-mixing. In contrast
to most of the existing literature (e.g., [40]), our analysis of
asymptotically independent rows takes into account the response variables and
(approximate) knockoff covariates, creating additional technical challenges.
We then apply the knockoffs inference method to the subsampled data. The
asymptotically independent rows in the subsampled data make it possible for
the KL divergence to asymptotically vanish. However, a major disadvantage of
subsampling is the inefficient use of data. We propose to apply the knockoffs
inference to each subsample of data, and then aggregate their results via the
e-value method [37, 32]. The resulting knockoffs inference method with
subsampling and e-value aggregation is presented in Section 2 and named as the
time series knockoffs inference (TSKI). We provide an explicit
characterization on how the serial dependence and the accuracy of the
approximate knockoffs affect the FDR control. Particularly, the TSKI procedure
can achieve asymptotic FDR control when the subsampling is done in an
appropriate fashion.

It is well-known that the FDR and power are two sides of the same coin. We
then turn to the power analysis of TSKI in Section 3. Assuming linear time
series models and regularity conditions, we show that TSKI without
subsampling, which is the e-value version of the original knockoff plus
inference [5, 4], enjoys the asymptotic sure screening property [18] under a
signal strength condition that is no stronger than the beta-min condition. For
TSKI with subsampling, our results reveal an interesting phenomenon that with
asymptotic probability one, the set of selected variables is either empty or
enjoys the sure screening property, when the tuning parameters are chosen
appropriately. The event of no discovery (i.e., the empty selected set)
usually occurs when most individual knockoff filters for different subsamples
select an abundantly large number of false positives so that the resulting
e-values for all variables become too small to be selected. We further show
that if the FDR for each individual knockoff filter can be controlled, then an
upper bound on the number of false positives for each knockoff filter can be
established. Consequently, the event of no discoveries can be ruled out, and
the set of selected variables by TSKI enjoys asymptotic power one. These
results are summarized formally in Theorem 3.

We test the empirical performance of TSKI on both simulated and real data. Our
simulation results in Section 4 demonstrate that the TSKI with subsampling
controls stably the time series FDR for some popular $\beta$-mixing processes.
The selection power is generally satisfactory with sufficient sample size.
Further, in Section 5 we apply the TSKI to study the temporal relations
between the inflation and other macroeconomic time series from the U.S.
economy during the recent ten years.

###  1.1 Notation

To facilitate the technical presentation, let us first introduce some
necessary notation. Let $(\Omega,\mathcal{F},\mathbb{P})$ and $\mathcal{R}$ be
the underlying probability space and the Borel $\sigma$-algebra on the real
line $\mathbb{R}$, respectively. We will use the boldface for random vectors
and matrices, and the tilde for the knockoff variables. The distribution of a
random mapping $\boldsymbol{X}$ is denoted as $\mu_{\boldsymbol{X}}$. For any
vector $\boldsymbol{x}$, define $\boldsymbol{x}_{-j}$ as $\boldsymbol{x}$
without the $j$th coordinate. We also use parentheses for matrix
concatenation. For a matrix $\boldsymbol{X}$ and an index subset $S$,
$\boldsymbol{X}(S)$ represents a submatrix of $\boldsymbol{X}$ containing only
columns with indices in $S$. For
$\vec{x}\coloneqq(x_{1},\cdots,x_{n})^{{}^{T}}\in\mathbb{R}^{n}$, we define
$\left\lVert\vec{x}\right\rVert_{k}\coloneqq(\sum_{i=1}^{n}|x_{i}|^{k})^{1/k}$,
$\left\lVert\vec{x}\right\rVert_{\infty}\coloneqq\max_{1\leq i\leq n}|x_{i}|$,
and
$\left\lVert\vec{x}\right\rVert_{0}\coloneqq\sum_{i=1}^{n}\boldsymbol{1}_{\\{x_{i}\neq
0\\}}$ with $\boldsymbol{1}_{\\{\cdot\\}}$ being the indicator function. Note
that we use the vector notation $\vec{x}$ to denote vectors in the Euclidean
space. The total variation (TV) norm for any measures $\mu_{1}$ and $\mu_{2}$
on $(\Omega,\mathcal{F})$ is defined as
$\left\lVert\mu_{1}-\mu_{2}\right\rVert_{TV}\coloneqq
2\sup_{\mathcal{D}\in\mathcal{F}}|\mu_{1}(\mathcal{D})-\mu_{2}(\mathcal{D})|$.
For any real sequences $\\{a_{n}\\}$ and $\\{b_{n}\\}$, $a_{n}=O(b_{n})$ means
$\limsup_{n\to\infty}\left|\frac{a_{n}}{b_{n}}\right|<\infty$, and
$a_{n}=o(b_{n})$ means
$\limsup_{n\to\infty}\left|\frac{a_{n}}{b_{n}}\right|=0$. We use $\\#S$ to
denote the cardinality of a given set $S$. Moreover, a transition kernel is
defined as a map
$p:(\mathbb{R}^{k_{1}},\mathcal{R}^{k_{2}})\longrightarrow[0,1]$ for some
positive integers $k_{1}$ and $k_{2}$ satisfying that (i) for each
$\mathcal{D}\in\mathcal{R}^{k_{2}}$, $p(\cdot,\mathcal{D})$ is a measurable
function and (ii) for each $x\in\mathbb{R}^{k_{1}}$, $p(x,\cdot)$ is a
probability measure.

##  2 Robust time series knockoffs inference with TSKI

Given an observed sample $\\{Y_{i},\boldsymbol{x}_{i}\\}_{i=1}^{n}$ that may
involve stationary time series variables, we are interested in selecting
relevant features in the covariate vector $\boldsymbol{x}_{1}$ with respect to
response $Y_{1}$ while controlling the false discovery rate (FDR), where the
definition of the null feature is given in Definition 1 below.

######  Definition 1.

(Null covariate) Consider the response $Y$ and the covariate vector
$\boldsymbol{x}=(X_{1},\cdots,X_{p})^{{}^{T}}$. Covariate $X_{j}$ with
$j\in\\{1,\cdots,p\\}$ is said to be null if and only if $X_{j}$ is
independent of response $Y$ conditional on
$(X_{1},\cdots,X_{j-1},X_{j+1},\cdots,X_{p})^{T}$.

We denote the set of null features in $\boldsymbol{x}_{i}$ with respect to
response $Y_{i}$ as $\mathcal{H}_{0}\subset\\{1,\cdots,p\\}$ for each $i$,
where $\mathcal{H}_{0}$ is independent of $i$ because
$(Y_{i},\boldsymbol{x}_{i})$’s are assumed to have the same distribution. For
a set of selected important features $\widehat{S}\subset\\{1,\cdots,p\\}$
returned by some algorithm, we are interested in controlling the false
discovery rate (FDR) defined as

$\textnormal{FDR}\coloneqq\mathbb{E}\left(\frac{\\#(\widehat{S}\cap\mathcal{H}_{0})}{(\\#\widehat{S})\vee
1}\right).$

In this paper, we will develop an inference procedure for controlling the FDR
for time series data with serial dependency. Our method builds on the recent
work of robust knockoffs inference [5] proposed for the non-time series
setting. We will introduce some key definitions motivated from [5] that will
be used throughout the presentation next.

###  2.1 Robust knockoffs inference framework

Throughout this paper, we will consider the knockoff generator given in
Definition 2 below.

######  Definition 2.

A knockoff generator
$\kappa:\mathbb{R}^{p}\times\mathcal{R}^{p}\longmapsto\mathbb{R}^{p}$ is a
transition kernel such that $\kappa(\vec{z},\cdot)$ is a probability measure
for each $\vec{z}\in\mathbb{R}^{p}$ and $\kappa(\cdot,\mathcal{A})$ is a
measurable function for each $\mathcal{A}\in\mathcal{R}^{p}$, where
$\mathcal{R}$ is the smallest $\sigma$-algebra of $\mathbb{R}$.

Given a knockoff generator and a random vector $\boldsymbol{x}$, we sample a
knockoff random vector $\widetilde{\boldsymbol{x}}$ of $\boldsymbol{x}$ from
the distribution $\kappa(\boldsymbol{x},\cdot)$. From now on, we let
$\widetilde{\boldsymbol{x}}_{i}$ be the knockoff vector of
$\boldsymbol{x}_{i}$ generated from the knockoff generator $\kappa$. Denote by
$\boldsymbol{X}=(\boldsymbol{x}_{1},\cdots,\boldsymbol{x}_{n})^{{}^{T}}$,
$\widetilde{\boldsymbol{X}}=(\widetilde{\boldsymbol{x}}_{1},\cdots,\widetilde{\boldsymbol{x}}_{n})^{{}^{T}}$,
and $\boldsymbol{Y}=(Y_{1},\cdots,Y_{n})^{{}^{T}}$. To facilitate our
technical analysis, we will introduce some regularity conditions.

######  Condition 1.

Assume that the density function of the distribution of
$(\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{Y})$ exists and
$(Y_{i},\boldsymbol{x}_{i})$’s are identically distributed. In addition, the
supports of $(\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{Y})$ and
$[\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{Y}]_{\textnormal{swap}(\\{j\\})}$
are the same for each $j\in\\{1,\cdots,p\\}$, where
$[\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{Y}]_{\textnormal{swap}(S)}$
is obtained by swapping the $j$th and $(j+p)$th columns of
$(\boldsymbol{X},\widetilde{\boldsymbol{X}},\boldsymbol{Y})$ for each $j\in
S\subset\\{1,\cdots,p\\}$.

######  Condition 2.

The knockoff generator $\kappa(\cdot,\cdot)$ is given before we observe
$\\{\boldsymbol{x}_{i},Y_{i}\\}_{i=1}^{n}$.

Condition 1 is a basic regularity condition. Condition 2 may be relaxed if we
apply the sample splitting to obtain an asymptotically independent training
subsample for training the knockoff generator. With Condition 2, $Y_{i}$ is
independent of $\widetilde{\boldsymbol{x}}_{i}$ conditional on
$\boldsymbol{x}_{i}$, and
$(Y_{i},\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i})$’s are identically
distributed if $(Y_{i},\boldsymbol{x}_{i})$’s are identically distributed. In
addition, without Condition 2,
$(\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i})$’s may not be i.i.d. even
if $\boldsymbol{x}_{i}$’s are i.i.d. themselves.

Let $(Y,\boldsymbol{x},\widetilde{\boldsymbol{x}})$ with
$\widetilde{\boldsymbol{x}}=(\widetilde{X}_{1},\cdots,\widetilde{X}_{p})^{{}^{T}}$
be an independent random vector identically distributed as
$(Y_{1},\boldsymbol{x}_{1},\widetilde{\boldsymbol{x}}_{1})$. A good knockoff
generator is expected to generate $\widetilde{\boldsymbol{x}}$ with
distribution as close to the model-X knockoffs as possible, where the model-X
knockoffs random vector
$\widetilde{\boldsymbol{x}}^{\textnormal{X}}\in\mathbb{R}^{p}$ of
$\boldsymbol{x}$ is such that 1) for each $S\subset\\{1,\cdots,p\\}$,
$[\boldsymbol{x},\widetilde{\boldsymbol{x}}^{\textnormal{X}}]_{\textnormal{swap}(S)}\overset{d}{=}(\boldsymbol{x},\widetilde{\boldsymbol{x}}^{\textnormal{X}})$
with $\overset{d}{=}$ standing for equal in distribution and 2)
$\widetilde{\boldsymbol{x}}^{\textnormal{X}}$ is independent of response $Y$
conditional on $\boldsymbol{x}$. However, generating the exact model-X
knockoffs requires the knowledge on the true underlying distribution of
$\boldsymbol{x}$, which is an impractical assumption for time series data. To
understand this, note that for time series data, covariates $\boldsymbol{x}$
may be lags of the response, and thus assuming a known distribution of
$\boldsymbol{x}$ gives us the serial distribution of the time series that
leads directly to known variable importance in such a scenario.

To overcome such difficulty, we will allow certain deviation of the
distribution of $\widetilde{\boldsymbol{x}}$ from that of the model-X
knockoffs. In other words, $\widetilde{\boldsymbol{x}}$ generated from
$\kappa(\boldsymbol{x},\cdot)$ is only a random vector of approximate knockoff
variables. The use of approximate knockoffs instead of the exact model-X
knockoffs may incur a potential FDR inflation. To control the FDR inflation,
we require that our knockoff generator $\kappa$ satisfies Condition 3 below,
which is the pairwise exchangeability [5]. Specifically, in addition to the
knockoff generator $\kappa(\cdot,\cdot)$, let us assume that we have a
coordinatewise knockoff generator
$\kappa_{j}:\mathbb{R}^{p-1}\times\mathcal{R}^{p-1}\longmapsto\mathbb{R}$ for
each $j\in\\{1,\cdots,p\\}$ that generates $\widetilde{X}_{j}^{\dagger}$ given
$\boldsymbol{x}_{-j}$, where the distribution of $\widetilde{X}_{j}^{\dagger}$
conditional on $\boldsymbol{x}_{-j}$ is an approximation to the true
underlying distribution of $X_{j}$ conditional on $\boldsymbol{x}_{-j}$.
Condition 3 is an assumption on all $p+1$ knockoff generators
$\kappa(\cdot,\cdot),\kappa_{1}(\cdot,\cdot),\cdots,\kappa_{p}(\cdot,\cdot)$.

######  Condition 3.

For each $1\leq j\leq p$, 1) if
$\widetilde{\boldsymbol{z}}=(\widetilde{Z}_{1},\cdots,\widetilde{Z}_{p})^{T}$
is sampled from the conditional distribution
$\kappa((X_{1},\cdots,X_{j-1},\widetilde{X}_{j}^{\dagger},X_{j+1},\cdots,X_{p}),\cdot)$,
then
$(\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j})$
and
$(\widetilde{Z}_{j},\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j})$
have the same distribution with $\widetilde{X}_{j}^{\dagger}$ sampled from
$\kappa_{j}(\boldsymbol{x}_{-j},\cdot)$, and 2) the density function of the
distribution of
$(\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j})$
exists.

To avoid potential confusion, throughout the paper we distinguish the knockoff
generator $\kappa(\cdot,\cdot)$ from $\kappa_{j}(\cdot,\cdot)$’s by referring
to the former as “the knockoff generator” and the latter as “the
coordinatewise knockoff generators.” Also, the dagger tilde notation denotes a
coordinatewise knockoff variable generated from $\kappa_{j}(\cdot,\cdot)$’s,
while the tilde notation represents a knockoff variable produced by the
knockoff generator $\kappa(\cdot,\cdot)$.

We would like to point out that for practical implementation, only
$\kappa(\cdot,\cdot)$ is needed for generating the knockoff variables, and
$\kappa_{j}(\cdot,\cdot)$’s are not needed. Nevertheless, from a theoretical
viewpoint, the existence of some $\kappa_{j}(\cdot,\cdot)$’s such that our
knockoff generator satisfies Condition 3 is crucial for controlling the FDR
inflation. Some examples of the knockoff generator that satisfy Condition 3
with respect to some coordinatewise knockoff generators are given in [5]. When
Condition 3 is satisfied, the FDR inflation can be measured by the
Kullback–Leibler (KL) divergence between the distribution of
$\widetilde{X}_{j}^{\dagger}$ conditional on $\boldsymbol{x}_{-j}$ and that of
$X_{j}$ conditional on $\boldsymbol{x}_{-j}$, which will be summarized
formally in Theorems 1–2 to be presented in the later sections.

###  2.2 Robust TSKI

Our robust time series knockoffs inference (TSKI) procedure applies the robust
knockoffs inference [5] to multiple subsamples of the given time series
observations and then establishes an ensemble inference using the e-values
[37]. Specifically, we consider $q+1$ subsamples and their index sets
$H_{k}=\\{k+s(q+1):s=0,1,\cdots,\lfloor\frac{n-k}{q+1}\rfloor\\}$ for each
$k\in\\{1,\cdots,q+1\\}$. To simplify the technical presentation, let
$\\{V_{i},\boldsymbol{u}_{i},\widetilde{\boldsymbol{u}}_{i}\\}_{i=1}^{N}$ be a
generic subsample of
$\\{Y_{i},\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i}\\}_{i=1}^{n}$.
Denote by $\boldsymbol{v}\coloneqq(V_{1},\cdots,V_{N})^{{}^{T}}$,
$\boldsymbol{U}\coloneqq(\boldsymbol{u}_{1},\cdots,\boldsymbol{u}_{N})^{{}^{T}}$,
and
$\widetilde{\boldsymbol{U}}\coloneqq(\widetilde{\boldsymbol{u}}_{1},\cdots,\widetilde{\boldsymbol{u}}_{N})^{{}^{T}}$.

Now the knockoff statistics $W_{j}$’s based on $\boldsymbol{v}$ and the
augmented design $(\boldsymbol{U},\widetilde{\boldsymbol{U}})$ are calculated
such that $W_{j}$ satisfies the sign-flip property [4, 11]: for each
$S\subset\\{1,\cdots,p\\}$ and each $1\leq j\leq p$,

$W_{j}(\boldsymbol{v},[\boldsymbol{U},\widetilde{\boldsymbol{U}}]_{\textnormal{swap}(S)})=\begin{cases}-W_{j}(\boldsymbol{v},\boldsymbol{U},\widetilde{\boldsymbol{U}})\quad\textnormal{
if }j\in S,\\\
W_{j}(\boldsymbol{v},\boldsymbol{U},\widetilde{\boldsymbol{U}})\quad\textnormal{
otherwise.}\end{cases}$ (2)

In fact, the condition can be weakened from “for each
$S\subset\\{1,\cdots,p\\}$” to “for each $S\subset\mathcal{H}_{0}$ with
$\\#S=1$,” but we choose not to pursue such a direction for technical
simplicity. Examples 1–2 below are two important examples of the knockoff
statistics that satisfy the sign-flip property. The random forests model in
Example 2 can be replaced with other learning models such as the deep learning
model.

######  Example 1 (Lasso coefficient difference (LCD)).

For a given sample
$(\boldsymbol{v},\boldsymbol{U},\widetilde{\boldsymbol{U}})$ and a tuning
parameter $\lambda\geq 0$, we define

$W_{j}=W_{j}(\boldsymbol{v},\boldsymbol{U},\widetilde{\boldsymbol{U}})=|\widehat{\beta}_{j}|-|\widehat{\beta}_{j+p}|,$

where
$(\widehat{\beta}_{1},\cdots,\widehat{\beta}_{2p})^{T}=\arg\min_{\boldsymbol{\beta}\in\mathbb{R}^{2p}}\left\\{n^{-1}\sum_{i=1}^{n}(V_{i}-(\boldsymbol{u}_{i}^{{}^{T}},\widetilde{\boldsymbol{u}}_{i}^{{}^{T}})\boldsymbol{\beta})^{2}+\lambda\sum_{j=1}^{2p}|\beta_{j}|\right\\}$
with $\boldsymbol{\beta}=(\beta_{1},\cdots,\beta_{2p})^{T}$.

######  Example 2 (Random forests mean decrease accuracy (MDA)).

For a given sample
$(\boldsymbol{v},\boldsymbol{U},\widetilde{\boldsymbol{U}})$, we define

$W_{j}=W_{j}(\boldsymbol{v},\boldsymbol{U},\widetilde{\boldsymbol{U}})=N^{-1}\sum_{i=1}^{N}\big{\\{}[V_{i}-\widehat{m}(\boldsymbol{u}_{i}^{(j)},\widetilde{\boldsymbol{u}}_{i})]^{2}-[V_{i}-\widehat{m}(\boldsymbol{u}_{i},\widetilde{\boldsymbol{u}}_{i}^{(j)})]^{2}\big{\\}}$

for each $j\in\\{1,\cdots,p\\}$, where
$(\boldsymbol{u}_{i}^{(j)},\widetilde{\boldsymbol{u}}_{i}^{(j)})=[\boldsymbol{u}_{i},\widetilde{\boldsymbol{u}}_{i}]_{\textnormal{swap}(\\{j\\})}$
and $\widehat{m}:\mathbb{R}^{2p}\longmapsto\mathbb{R}$ is the random forests
regression function trained by regressing $V_{i}$’s on
$(\boldsymbol{u}_{i},\widetilde{\boldsymbol{u}}_{i})$’s.

1 Let $0<\tau_{1}<1$ be a constant and $0<\tau^{*}<1$ the target FDR level.

2For each $k\in\\{1,\cdots,q+1\\}$, calculate the knockoff statistics
$W_{1}^{k},\cdots,W_{p}^{k}$ satisfying (2) using sample
$\\{\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i},Y_{i}\\}_{i\in H_{k}}$.

3Calculate the e-value statistics $e_{j}=(q+1)^{-1}\sum_{k=1}^{q+1}e_{j}^{k}$,
where111$\min\emptyset$ and $\max\emptyset$ are defined as $\infty$ and $0$,
respectively.
$\begin{split}e_{j}^{k}&=\frac{p\times\boldsymbol{1}_{\\{W_{j}^{k}\geq
T^{k}\\}}}{1+\sum_{s=1}^{p}\boldsymbol{1}_{\\{W_{s}^{k}\leq-T^{k}\\}}},\\\
T^{k}&=\min\left\\{t\in\mathcal{W}_{+}^{k}:\frac{1+\\#\\{j:W_{j}^{k}\leq-t\\}}{\\#\\{j:W_{j}^{k}\geq
t\\}\vee 1}\leq\tau_{1}\right\\},\end{split}$ (3) and
$\mathcal{W}_{+}^{k}=\\{|W_{s}^{k}|:|W_{s}^{k}|>0\\}$ for each
$k\in\\{1,\cdots,q+1\\}$.

4Let $\widehat{S}=\\{j:e_{j}\geq p(\tau^{*}\times\widehat{k})^{-1}\\}$ with
$\widehat{k}=\max\\{k:e_{(k)}\geq p(\tau^{*}\times k)^{{-1}}\\}$, where
$e_{(j)}$’s are the ordered statistics of $e_{j}$’s such that
$e_{(1)}\geq\cdots\geq e_{(p)}$.

5

Algorithm 1 Robust time series knockoffs inference (TSKI) via e-values

Algorithm 1 above provides a way to combine the sets of variables selected by
multiple individual knockoff filters via the e-value approach. The e-BH
procedure is needed here since naively taking the intersection or union over
multiple selected sets by the knockoff filters (i.e., $\\{j:W_{j}^{k}\geq
T^{k}\\}$ with $k=1,\cdots,q+1$) would not be guaranteed to control the FDR.
Algorithm 1 can be viewed as a sophisticated majority voting procedure. In
particular, if $\widehat{S}\not=\emptyset$, then $\widehat{S}$ includes the
intersection of variable sets selected by individual knockoff filters as a
subset and $\widehat{S}$ is also a subset of the union of all variable sets
selected by individual knockoff filters; that is,

$\cap_{k=1}^{q+1}\\{j:W_{j}^{k}\geq
T^{k}\\}\subset\widehat{S}\subset\cup_{k=1}^{q+1}\\{j:W_{j}^{k}\geq T^{k}\\}.$

Meanwhile, the tuning parameter $\tau_{1}$ should be set to be smaller than
$\tau^{*}$ in order to reduce the selection power loss compared to the
underlying knockoffs procedure. See Theorem 3 in Section 3 for details.

######  Condition 4 ($h$-step $\beta$-mixing with exponential decay).

Assume that the process $\\{\boldsymbol{x}_{t}\\}$ is a $p$-dimensional
stationary Markov chain with a transition kernel
$p:\mathbb{R}^{p}\times\mathcal{R}^{p}\longmapsto\mathbb{R}$ and a stationary
distribution $\pi$. There exist a positive integer $h$, a measurable function
$V:\mathbb{R}^{p}\longrightarrow[0,\infty)$, and some constants $0\leq\rho<1$
and $0<C_{0}<\infty$ such that for each $\vec{x}\in\mathbb{R}^{p}$,

$\left\lVert p^{h}(\vec{x},\cdot)-\pi(\cdot)\right\rVert_{TV}\leq
C\rho^{h}V(\vec{x}),$

where $C>0$ is some constant with $C_{0}\geq
C\int_{\mathbb{R}^{p}}V(\vec{x})d\pi(\vec{x})$ and
$\left\lVert\cdot\right\rVert_{TV}$ denotes the total variation (TV) norm
associated with measures. Moreover, for each $\vec{x}\in\mathbb{R}^{p}$,
$p(\vec{x},\cdot)$ is absolutely continuous with respect to the Lebesgue
measure.

Let
$\\{\boldsymbol{x}_{i}^{\pi},\widetilde{\boldsymbol{x}}_{i}^{\pi},Y_{i}^{\pi}\\}_{i=1}^{n}$
be a sequence of i.i.d. random vectors such that
$(\boldsymbol{x}_{1}^{\pi},\widetilde{\boldsymbol{x}}_{1}^{\pi},Y_{1}^{\pi})$
and $(\boldsymbol{x}_{1},\widetilde{\boldsymbol{x}}_{1},Y_{1})$ have the same
distribution,
$\mathcal{X}_{k}=\\{\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i},Y_{i}\\}_{i\in
H_{k}}$, and
$\mathcal{X}_{k}^{\pi}=\\{\boldsymbol{x}_{i}^{\pi},\widetilde{\boldsymbol{x}}_{i}^{\pi},Y_{i}^{\pi}\\}_{i\in
H_{k}}$ for each $k\in\\{1,\cdots,q+1\\}$. Let $f_{\boldsymbol{z}}(\cdot)$ be
the density function of random vector $\boldsymbol{z}$. The FDR inflation due
to the knockoff variable approximation can be measured and controlled by the
KL divergence given in Theorem 1 below.

######  Theorem 1.

Let $\widehat{S}$ be the set of selected variables by the TSKI with Algorithm
1. Then under Condition 1 and the assumption of positive $T^{k}$’s in (3), we
have

$\begin{split}\textnormal{FDR}&\leq\inf_{\varepsilon>0}\Big{[}\tau^{*}e^{\varepsilon}+\sum_{k=1}^{q+1}\mathbb{P}(\max_{1\leq
j\leq
p}\widehat{\textnormal{KL}}_{j}^{k\pi}>\varepsilon)\Big{]}+\sum_{k=1}^{q+1}\sup_{\mathcal{D}\in\mathcal{R}^{\\#H_{k}\times(2p+1)}}|\mathbb{P}(\mathcal{X}_{k}\in\mathcal{D})-\mathbb{P}(\mathcal{X}_{k}^{\pi}\in\mathcal{D})|,\end{split}$
(4)

where $0<\tau^{*}<1$ is the target FDR level and for each $1\leq k\leq q+1$
and $1\leq j\leq p$,

$\begin{split}\widehat{\textnormal{KL}}_{j}^{k\pi}&=\sum_{i\in
H_{k}}\log{\left(\frac{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(X_{ij}^{\pi},\widetilde{X}_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi},\widetilde{\boldsymbol{x}}_{-ij}^{\pi})}{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(\widetilde{X}_{ij}^{\pi},X_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi},\widetilde{\boldsymbol{x}}_{-ij}^{\pi})}\right)}.\end{split}$
(5)

The proof of Theorem 1 is provided in Section A.1 of the Supplementary
Material. Recall that $(\boldsymbol{x},\widetilde{\boldsymbol{x}},Y)$ and
$(\boldsymbol{x}_{1},\widetilde{\boldsymbol{x}}_{1},Y_{1})$ have the same
distribution. The KL divergence given in (5) that appears on the right-hand
side of (4) can be simplified further below provided that some additional
conditions are satisfied.

######  Corollary 1.

Assume that all the conditions of Theorem 1 hold. If further Condition 3 is
satisfied, then (4) holds with

$\begin{split}\widehat{\textnormal{KL}}_{j}^{k\pi}&=\sum_{i\in
H_{k}}\log{\left(\frac{f_{X_{j},\boldsymbol{x}_{-j}}(X_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi})f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(\widetilde{X}_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi})}{f_{X_{j},\boldsymbol{x}_{-j}}(\widetilde{X}_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi})f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(X_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi})}\right)}\\\
&=\sum_{i\in
H_{k}}\log{\left(\frac{f_{X_{j}|\boldsymbol{x}_{-j}}(X_{ij}^{\pi}|\boldsymbol{x}_{-ij}^{\pi})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(\widetilde{X}_{ij}^{\pi}|\boldsymbol{x}_{-ij}^{\pi})}{f_{X_{j}|\boldsymbol{x}_{-j}}(\widetilde{X}_{ij}^{\pi}|\boldsymbol{x}_{-ij}^{\pi})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(X_{ij}^{\pi}|\boldsymbol{x}_{-ij}^{\pi})}\right)},\end{split}$

where $\widetilde{X}_{j}^{\dagger}$’s are given in Condition 3 and
$f_{\boldsymbol{z}_{1}|\boldsymbol{z}_{2}}(\boldsymbol{z}_{1}|\boldsymbol{z}_{2})=f_{\boldsymbol{z}_{1},\boldsymbol{z}_{2}}(\boldsymbol{z}_{1},\boldsymbol{z}_{2})[f_{\boldsymbol{z}_{2}}(\boldsymbol{z}_{2})]^{-1}$
which is the conditional probability density function of $\boldsymbol{z}_{1}$
given $\boldsymbol{z}_{2}$.

######  Corollary 2.

Assume that all the conditions of Theorem 1 hold. If further Condition 2 is
satisfied, $\\{\boldsymbol{x}_{i}\\}_{i\geq 1}$ satisfies Condition 4 with
$q$-step and constants $C_{0}>0$ and $0\leq\rho<1$, and $Y_{i}$ is
$\boldsymbol{x}_{i+1}$-measurable, then (4) holds with

$\sum_{k=1}^{q+1}\sup_{\mathcal{D}\in\mathcal{R}^{\\#H_{k}\times(2p+1)}}|\mathbb{P}(\mathcal{X}_{k}\in\mathcal{D})-\mathbb{P}(\mathcal{X}_{k}^{\pi}\in\mathcal{D})|\leq
C_{0}\times\rho^{q}\times n.$ (6)

Moreover, when $(Y_{i},\boldsymbol{x}_{i})$’s are i.i.d., (6) holds with
$\rho=0$.

To make our results more concrete, let us consider an example with decreasing
KL divergence. Assume that $\\{\boldsymbol{x}_{t}\\}$ follows a stationary
linear Gaussian process as in Example 3 in Section 2.3.1 with zero mean and
covariance matrix
$\Theta=\mathbb{E}(\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{{}^{T}})$, and the
knockoff generator is such that $\kappa(\vec{z},\cdot)$ follows a Gaussian
distribution with mean $(\boldsymbol{I}_{p}-D\widehat{\Theta})\vec{z}$ and
variance $2D-D\widehat{\Theta}D$ for each $\vec{z}\in\mathbb{R}^{p}$, where
$\boldsymbol{I}_{p}$ is the $p$-dimensional identity matrix,
$\widehat{\Theta}$ is the estimated covariance matrix, and $D$ is a diagonal
matrix with nonnegative entries such that $2D-D\widehat{\Theta}D$ is positive
semidefinite. Then it has been shown in Lemma 5 of [5] that when the data
consists of i.i.d. observations without serial dependency, it holds that for
each $\varepsilon>0$,

$\lim_{n\rightarrow\infty}\sum_{k=1}^{q+1}\mathbb{P}(\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}^{k\pi}>\varepsilon)=0$ (7)

as long as Condition 3 is satisfied, $p\gg q$, and

$\lim_{n\rightarrow\infty}\left\\{\max_{1\leq j\leq
p}\Theta_{jj}^{-\frac{1}{2}}\left\lVert\Theta^{-\frac{1}{2}}(\widehat{\Theta}_{j}-\Theta_{j})\right\rVert_{2}\times\sqrt{n\log{p}}\right\\}=0,$

where $\Theta_{jj}$ denotes the $j$th diagonal entry of $\Theta$ and
$\Theta_{j}$ represents the $j$th column of $\Theta$. We omit the dependence
of the parameters on sample size $n$ here for simplicity. More examples on
asymptotically vanishing KL divergence for non-time series data can be found
in the same paper above. Similar results can be proved for our applications of
time series data using the proof techniques in [5] by replacing the
concentration inequalities for i.i.d. observations with those for
$\beta$-mixing data; we omit the details here because the proof modification
is straightforward.

###  2.3 Stationary processes satisfying Condition 4

####  2.3.1 Linear ARX process

Example 3 below is a benchmark model for describing the behavior of
macroeconomic variables.

######  Example 3 (Autoregressive models with exogenous variables (ARX)).

For each $t$, we define

$Y_{t}=\sum_{j=1}^{k_{1}}\alpha_{j}Y_{t-j}+\sum_{l=1}^{k_{2}}\sum_{j=1}^{k_{3l}}\beta_{j}^{(l)}H_{t-j+1}^{(l)}+\varepsilon_{t}$

with
$H_{t}^{(l)}=\epsilon_{t}^{(l)}+\sum_{j=1}^{k_{3l}}b_{j}^{(l)}H_{t-j}^{(l)}$,
where for some positive constants $L_{0},L_{1}$, and $L_{2}$, it satisfies
that $1-\sum_{j=1}^{k_{3l}}b_{j}^{(l)}z^{j}\not=0$ and
$1-\sum_{j=1}^{k_{1}}\alpha_{j}z^{j}\not=0$ for each $|z|\leq 1+L_{0}$ and
each $l\in\\{1,\cdots,k_{2}\\}$;
$\sum_{l=1}^{k_{2}}\sum_{j=1}^{k_{3l}}|\beta_{j}^{(l)}|<L_{1}$; and
$k_{1},k_{31},\cdots,k_{3k_{2}}<L_{2}$. Here, $k_{1}$, $k_{2}$, and
$k_{31},\cdots,k_{3k_{2}}$ are all positive integers. Moreover,
$(\epsilon_{t}^{(1)},\cdots$, $\epsilon_{t}^{(k_{2})},\varepsilon_{t})$’s are
$(k_{2}+1)$-dimensional i.i.d. Gaussian random vectors with zero mean and
positive definite covariance matrix $\Sigma_{0}$ that satisfy
$\mathbb{E}(\varepsilon_{t}\epsilon_{t}^{(l)})=0$ for each $l$.

It is seen that the covariate vector with respect to response $Y_{t}$ is
$\boldsymbol{x}_{t}=(Y_{t-1},\cdots,Y_{t-k_{1}},\boldsymbol{h}_{t}^{{}^{T}})^{{}^{T}}$
with
$\boldsymbol{h}_{t}=(H_{t}^{(1)},\cdots,H_{t-k_{31}+1}^{(1)},\cdots,H_{t}^{(k_{2})},\cdots,H_{t-k_{3k_{2}}+1}^{(k_{2})})^{{}^{T}}$.
It has been shown (e.g., [2]) that the stationary process
$\\{\boldsymbol{x}_{t}\\}$ in Example 3 with fixed dimensionality satisfies
Condition 4 with $h$-step and some constants $\rho,C_{0}$ for each positive
integer $h$. When the dimensionality of $\boldsymbol{x}_{t}$ increases,
certain growth conditions such as (8) below on the value of $h$ and the
dimensionality of $\boldsymbol{x}_{t}$ are needed for Condition 4 to hold.
Hence, we make the dependence of the stationary processes on $h$ explicit:
hereafter $\\{\boldsymbol{x}_{t}^{(h)}\\}$ denotes a $p_{h}$-dimensional
stationary process satisfying Condition 4 for $h=1,2,\cdots$. Proposition 1
below provides some sufficient conditions on the growth rate of $p_{h}$ as a
function of $h$.

######  Proposition 1.

Let $\\{\boldsymbol{x}_{t}^{(h)}\\}$ be a sequence of $p_{h}$-dimensional
linear process in Example 3 with constant $L_{i}$’s and uniformly positive
definite $\Sigma_{0}^{(h)}$’s. Assume that for some constant $C_{2}>0$ and
sufficiently small $s_{2}>0$,

$\sup_{h>0}\\{p_{h}\exp{(-s_{2}h)}\\}\leq C_{2}.$ (8)

Then for some constants $0\leq\rho<1$ and $C_{0}>0$ and all large $h$,
$\\{\boldsymbol{x}_{t}^{(h)}\\}$ satisfies Condition 4 with $h$-step.

Let arbitrary constants $K_{0}>0$ and $\delta>0$ be given. In view of Theorem
1 and Proposition 1, if we perform subsampling such that
$q=q_{n}=\lceil(\log{n})^{1+\delta}\rceil$ in the TSKI algorithm (Algorithm 1)
and assume Example 3 with $p=p_{q_{n}}=O(n^{K_{0}})$, then the $\beta$-mixing
convergence rate required by Theorem 1 (Condition 4 with $q_{n}$-step and
constants $C_{0},\rho$) is satisfied by the subsampled data for all large $n$
according to Proposition 1.

####  2.3.2 Various nonlinear AR-type processes

Other than Example 3, many time-homogeneous Markov chains also satisfy
Condition 4. To name a few, [35, 2] showed that with some additional mild
regularity conditions, $\\{(Y_{t-1},\cdots,Y_{t-k_{1}})\\}$ of Example 4 below
satisfies Condition 4 for all $h>0$ with some constants $C_{0}$ and
$0\leq\rho<1$.

######  Example 4 (Nonlinear AR [35, 2]).

Let a measurable function $G:\mathbb{R}^{k_{1}}\longrightarrow\mathbb{R}$ be
given such that $\sup_{\vec{z}\in\mathbb{R}^{k_{1}}}|G(\vec{z})|<\infty$, and
$\\{\varepsilon_{t}\\}$ a sequence of i.i.d. model errors. For each $t$, we
define

$Y_{t}=G(Y_{t-1},\cdots,Y_{t-k_{1}})+\varepsilon_{t}.$

The self-exciting threshold autoregressive models (SETAR) [36, 23] also
satisfies Condition 4 for all $h>0$ according to [2]. For more examples such
as the exponential AR models, see [31, 2] and the references therein.

####  2.3.3  ARCH-type process

######  Example 5 (AR($k_{1}$)-X-ARCH($k_{3}$) [14, 30]).

Let $\varepsilon_{t}$’s be i.i.d. random variables with zero mean and
$\mathbb{E}\varepsilon_{1}^{2}=1$, and
$\\{\boldsymbol{h}_{t}\coloneqq(H_{t,1},\cdots,H_{t,k_{2}})\\}$ a sequence of
$k_{2}$-dimensional time series covariates that is independent of
$\varepsilon_{t}$’s. Let measurable functions
$G_{1}:\mathbb{R}^{k_{3}}\rightarrow(0,\infty)$,
$G_{2}:\mathbb{R}^{k_{2}}\rightarrow\mathbb{R}$, and
$\gamma_{j}:\mathbb{R}^{k_{1}}\rightarrow\mathbb{R}$ where $1\leq j\leq k_{1}$
be given such that
$\sum_{j=1}^{k_{1}}\sup_{\vec{z}\in\mathbb{R}^{k_{1}}}|\gamma_{j}(\vec{z})|<1$.
The functional-coefficient ARX-ARCH model [12] is given by

$Y_{t}=\sum_{j=1}^{k_{1}}\gamma_{j}(Y_{t-1},\cdots,Y_{t-k_{1}})Y_{t-j}+G_{2}(\boldsymbol{h}_{t})+\sigma_{t}\varepsilon_{t}$

with
$\sigma_{t}=G_{1}(\sigma_{t-1}\varepsilon_{t-1},\cdots,\sigma_{t-k_{3}}\varepsilon_{t-k_{3}})$.

Example 5 above is an ARCH model with the mean function consisting of an AR
component and exogenous covariates. The covariate vector is
$\boldsymbol{x}_{t}=(Y_{t-1},\cdots,Y_{t-k_{1}-k_{3}},\boldsymbol{h}_{t},\cdots,\boldsymbol{h}_{t-k_{3}})^{{}^{T}}$
for response $Y_{t}$. In particular, Model 3 in Section 4 is an example of the
ARX-ARCH model above when $\gamma_{j}$’s are constants and the model error
follows a standard ARCH process [17].

The AR component in Example 5 above can be a general functional-coefficient
autoregressive model [12], and the ARCH component can take the form of smooth
transition ARCH model [26]. With $G_{2}=0$ and additional mild regularity
conditions, the Markov chain $\\{Y_{t},\cdots,Y_{t-k_{1}-k_{3}+1}\\}_{t}$
satisfies Condition 4 for each $h>0$ with constants $C_{0}>0$ and
$0\leq\rho<1$ according to [14, 30]. In the presence of exogenous covariates
$\boldsymbol{h}_{t}$, additional regularity conditions on
$\boldsymbol{h}_{t}$’s are needed for Condition 4 to hold for the Markov
chain; such study is beyond the scope of the current paper and is left for
future investigation.

One challenge in variable selection problem with time series data is that
there does not always exist an obvious definition of the covariate vector.
Taking Example 5 for instance, the existence of the ARCH component requires us
to take into account both the mean function and variance function when
selecting the set of non-null variables in the broad sense according to
Definition 1. To better understand this, let us consider an ARX-ARCH(1) model
with a standard ARCH component such that
$\sigma_{t}=\sqrt{0.1+0.9(\sigma_{t-1}\varepsilon_{t-1})^{2}}$, and write

$\sigma_{t-1}\varepsilon_{t-1}=Y_{t-1}-\sum_{j=1}^{k_{1}}\gamma_{j}(Y_{t-2},\cdots,Y_{t-1-k_{1}})Y_{t-1-j}-G_{2}(\boldsymbol{h}_{t-1}).$
(9)

In this example, in addition to variables
$(Y_{t-1},\cdots,Y_{t-k_{1}},\boldsymbol{h}_{t})$ which affect the mean
regression function, we should also take into account the lagged covariates in
the ARCH component $(Y_{t-1},\cdots,Y_{t-k_{1}-1},\boldsymbol{h}_{t-1})$ when
conducting variable selection in the broad sense according to Definition 1.
That is, one sensible choice of the covariate vector is
$\boldsymbol{x}_{t}=(Y_{t-1},\cdots,Y_{t-k_{1}-1},\boldsymbol{h}_{t},\boldsymbol{h}_{t-1})^{{}^{T}}$.
Omitting variables in the variance function (i.e., the ARCH component) and
defining the covariate vector as
$(Y_{t-1},\cdots,Y_{t-k_{1}},\boldsymbol{h}_{t})$ may give us a nonsparse set
of non-null variables according to Definition 1. Nevertheless, the actual
variable selection performance of the TSKI depends on the specific choice of
the knockoff statistics, as will be shown in the later sections. If the
knockoff statistics are constructed based on the mean regression function
alone (e.g., the LCD and MDA discussed earlier), then the corresponding TSKI
cannot be expected to have power in selecting variables that affect only the
variance function. In this sense, our results in such scenario should be
interpreted as for selecting the important variables contributing to the mean
regression function alone.

######  Remark 1.

We have left out the GARCH-type process [9] in our discussion because it can
be challenging to formulate meaningful covariates for variable selection
purpose in such a setting. Note that a GARCH-type process can be represented
as an ARCH-type process with infinite order. Thus, if the covariates vector is
not well-formulated such that some active covariates are not included, the
resulting regression model may no longer be sparse, rendering the FDR control
problem invalid. We shall leave the variable selection problem for GARCH-type
process in future work.

###  2.4 Robust TSKI without subsampling

In this section, we consider a special case of the TSKI with Algorithm 1 when
we set $q=0$, that is, no subsampling is used. We provide a full description
of the corresponding algorithm in Algorithm 2 below since the related method
and theory have standalone value to the time series literature. Two major
contributions here are 1) extending the robust knockoffs inference [5] to the
e-value setting and 2) further extending the results to time series
applications. By similar analysis as in Theorem 2 below, we can show that the
robust knockoffs inference [5] (without using the e-values) can also be
extended to time series applications, but the details are omitted here for
simplicity. We emphasize that our results (11)–(12) in Theorem 2 below assume
neither i.i.d. observations nor the pairwise exchangibility Condition 3.

1 Let $0<\tau_{1}<1$ be a constant and $0<\tau^{*}<1$ the target FDR level.

2Calculate the knockoff statistics $W_{1},\cdots,W_{p}$ satisfying (2) with
the full sample
$\\{Y_{i},\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i}\\}_{i=1}^{n}$.

3Let $\mathcal{W}_{+}=\\{|W_{s}|:|W_{s}|>0\\}$. Calculate the e-value
statistics $e_{j}$’s such that222$\min\emptyset$ and $\max\emptyset$ are
defined as $\infty$ and $0$, respectively.
$\begin{split}e_{j}=\frac{p\times\boldsymbol{1}_{\\{W_{j}\geq
T\\}}}{1+\sum_{s=1}^{p}\boldsymbol{1}_{\\{W_{s}\leq-T\\}}},\ \
T=\min\left\\{t\in\mathcal{W}_{+}:\frac{1+\\#\\{j:W_{j}\leq-t\\}}{\\#\\{j:W_{j}\geq
t\\}\vee 1}\leq\tau_{1}\right\\}.\end{split}$ (10)

4Let $\widehat{S}=\\{j:e_{j}\geq p(\tau^{*}\times\widehat{k})^{-1}\\}$ with
$\widehat{k}=\max\\{k:e_{(k)}\geq p(\tau^{*}\times k)^{{-1}}\\}$, where
$e_{(j)}$’s are the ordered statistics of $e_{j}$’s such that
$e_{(1)}\geq\cdots\geq e_{(p)}$.

Algorithm 2 Robust time series knockoffs inference (TSKI) via e-values without
subsampling

Let $\boldsymbol{X}_{-j}$ be the submatrix of $\boldsymbol{X}$ with the $j$th
column removed, and $\boldsymbol{X}_{j}$ and $\widetilde{\boldsymbol{X}}_{j}$
the $j$th columns of $\boldsymbol{X}$ and $\widetilde{\boldsymbol{X}}$,
respectively. Recall that $(Y,\boldsymbol{x},\widetilde{\boldsymbol{x}})$ is
an independent copy of
$(Y_{1},\boldsymbol{x}_{1},\widetilde{\boldsymbol{x}}_{1})$ and
$\widetilde{X}_{j}^{\dagger}$ is given in Condition 3 by the $j$th
coordinatewise knockoff generator.

######  Theorem 2.

Let $\widehat{S}$ be the set of selected variables by the TSKI with Algorithm
2 and $0<\tau^{*}<1$ the target FDR level. Assume that Condition 1 holds and
$T$ in (10) is positive. Then we have

$\begin{split}\textnormal{FDR}&\leq\inf_{\varepsilon>0}\Big{[}\tau^{*}\times
e^{\varepsilon}+\mathbb{P}(\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}>\varepsilon)\Big{]},\end{split}$ (11)

where for each $1\leq j\leq p$,

$\widehat{\textnormal{KL}}_{j}=\log{\left(\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}\right)}.$
(12)

If we further assume that Condition 2 is satisfied and $\boldsymbol{X}_{j}$ is
independent of $\boldsymbol{Y}$ conditional on $\boldsymbol{X}_{-j}$ for each
$j\in\mathcal{H}_{0}$, then we have

$\widehat{\textnormal{KL}}_{j}=\log{\left(\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j}}(\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j})}{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j}}(\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j})}\right)}.$
(13)

Moreover, if
$(\boldsymbol{x},\widetilde{\boldsymbol{x}}),(\boldsymbol{x}_{1},\widetilde{\boldsymbol{x}}_{1}),\cdots,(\boldsymbol{x}_{n},\widetilde{\boldsymbol{x}}_{n})$
are further assumed to be i.i.d., then we have

$\begin{split}\widehat{\textnormal{KL}}_{j}&=\sum_{i=1}^{n}\log{\left(\frac{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(X_{ij},\widetilde{X}_{ij},\boldsymbol{x}_{-ij},\widetilde{\boldsymbol{x}}_{-ij})}{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(\widetilde{X}_{ij},X_{ij},\boldsymbol{x}_{-ij},\widetilde{\boldsymbol{x}}_{-ij})}\right)}.\end{split}$
(14)

If further Condition 3 is satisfied, then we have

$\begin{split}\widehat{\textnormal{KL}}_{j}&=\sum_{i=1}^{n}\log{\left(\frac{f_{X_{j},\boldsymbol{x}_{-j}}(X_{ij},\boldsymbol{x}_{-ij})f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(\widetilde{X}_{ij},\boldsymbol{x}_{-ij})}{f_{X_{j},\boldsymbol{x}_{-j}}(\widetilde{X}_{ij},\boldsymbol{x}_{-ij})f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(X_{ij},\boldsymbol{x}_{-ij})}\right)}\\\
&=\sum_{i=1}^{n}\log{\left(\frac{f_{X_{j}|\boldsymbol{x}_{-j}}(X_{ij}|\boldsymbol{x}_{-ij})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(\widetilde{X}_{ij}|\boldsymbol{x}_{-ij})}{f_{X_{j}|\boldsymbol{x}_{-j}}(\widetilde{X}_{ij}|\boldsymbol{x}_{-ij})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(X_{ij}|\boldsymbol{x}_{-ij})}\right)},\end{split}$
(15)

where
$f_{\boldsymbol{z}_{1}|\boldsymbol{z}_{2}}(\boldsymbol{z}_{1}|\boldsymbol{z}_{2})$
denotes the conditional probability density function of $\boldsymbol{z}_{1}$
given $\boldsymbol{z}_{2}$.

The proof of Theorem 2 follows mainly those in [5, 32, 37] and is presented in
Section A.2 of the Supplementary Material. Comparing (15) to (14), we see that
the KL divergences become invariant to $\widetilde{\boldsymbol{x}}_{-j}$
thanks to the additional assumption Condition 3. The simplified form in (15)
is important in proving the asymptotic FDR control as in (7). In addition,
Condition 3 allows for deviation of the conditional distribution of
$\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}$ from the true underlying
conditional distribution of $X_{j}|\boldsymbol{x}_{-j}$, making the procedure
more practically applicable, as verified in examples given in [5].

##  3 Selection power analysis of TSKI under linear models

In this section, we will analyze the selection power of TSKI with Algorithms
1–2 in time series data applications. Since the selection power of any
procedure would depend on the signal strength in the data, we showcase the
power analysis using linear time series models where the signal strength can
be measured conveniently by regression coefficients. Correspondingly, we
consider the LCD knockoff statistic which is used popularly in linear models
for variable selection.

Given the augmented data
$\\{Y_{i},\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i}\\}_{i=1}^{n}$ with
both original covariate vectors $\boldsymbol{x_{i}}$’s and their knockoff
copies $\widetilde{\boldsymbol{x}}_{i}$’s, the Lasso estimate of the
regression coefficient vector with tuning parameter $\lambda_{n}\geq 0$ is
given by

$(\widehat{\beta}_{1},\cdots,\widehat{\beta}_{2p})^{T}=\arg\min_{\vec{\beta}\in\mathbb{R}^{2p}}\left\\{n^{-1}\sum_{i=1}^{n}(Y_{i}-(\boldsymbol{x}_{i}^{{}^{T}},\widetilde{\boldsymbol{x}}_{i}^{{}^{T}})\vec{\beta})^{2}+\lambda_{n}\sum_{j=1}^{2p}|\beta_{j}|\right\\},$

where $\vec{\beta}=(\beta_{1},\cdots\beta_{2p})^{T}$. We also denote by
$\vec{\beta}^{*}:=(\beta_{1}^{*},\cdots,\beta_{2p}^{*})^{T}$ and
$S^{*}=\\{j:|\beta_{j}^{*}|>0\\}$ the true coefficient vector and the set of
relevant features, respectively, in the linear regression model

$Y_{i}=(\boldsymbol{x}_{i}^{T},\widetilde{\boldsymbol{x}}_{i}^{T})\vec{\beta}^{*}+\varepsilon_{i},$

where $\varepsilon_{i}$ is the model error. It is seen that $\beta_{j}^{*}=0$
for all $j=p+1,\cdots,2p$ because they correspond to the knockoff variables.

For the time series setting when
$(Y_{i},\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i})$ exhibits serial
dependency across $i$, the Lasso estimation error has been investigated
thoroughly under various model settings and serial dependence structures; see,
for example, [25, 39, 6, 29, 38, 1]. Specifically, the linear AR and ARCH
model with sub-Weibull model errors has been considered in [38], the AR
distributed lag models with GARCH errors are discussed in [1, 29], and the
misspecified AR models and AR models with weak sparsity have been studied in
[1]. Rather than reproducing these existing results, we summarize the generic
condition needed for the TSKI power analysis in Condition 5 below. We remark
that Condition 5 have been proved to hold in the aforementioned works with
appropriately chosen $\lambda_{n}$ and their corresponding regularity
conditions.

######  Condition 5.

For some constant $c_{0}>0$ and sequence $k_{3n}>0$ with
$\lim_{n\rightarrow\infty}k_{3n}=0$, it holds that
$\mathbb{P}\left(\sum_{j=1}^{2p}|\widehat{\beta}_{j}-\beta_{j}^{*}|\leq
c_{0}(\\#S^{*})\lambda_{n}\right)\geq 1-k_{3n}$.

######  Condition 6.

There exists some sequence $k_{1n}>0$ with $k_{1n}q^{-1}\rightarrow\infty$ as
$n\rightarrow\infty$ such that $\min_{j\in
S^{*}}|\beta_{j}^{*}|>k_{1n}\lambda_{n}$.

######  Condition 7.

For some constant $c_{1}\in(0,1)$ and sequence $\\{k_{2n}\\}$ with
$\lim_{n\rightarrow\infty}k_{2n}=0$, it holds that
$2(\tau_{1}\\#S^{*})^{-1}<c_{1}$, $\mathbb{P}(\\#\\{j:W_{j}\geq T\\}\geq
c_{1}(\\#S^{*}))\geq 1-k_{2n}$ for Algorithm 2, and
$\mathbb{P}(\\#\\{j:W_{j}^{k}\geq T^{k}\\}\geq c_{1}(\\#S^{*}))\geq 1-k_{2n}$
with $k\in\\{1,\cdots,q+1\\}$ for Algorithm 1.

Recall that $q+1$ is the number of subsamples in Algorithm 1, $W_{j}^{k}$’s
and $T^{k}$’s are given in Algorithm 1, and $W_{j}$’s and $T$’s are given in
Algorithm 2. Also, we set $q=O((\log{n})^{1+\delta})$ for some small
$\delta>0$ in our applications. Thus, when $\\#S^{*}\gg(\log{n})^{1+\delta}$,
Condition 6 is weaker than the popularly imposed beta-min assumption which
requires that $\min_{j\in S^{*}}|\beta_{j}^{*}|>(\\#S^{*})\times\lambda_{n}$
for ensuring the Lasso model selection consistency.

Condition 7 is a high-level assumption that can be shown to hold when the
following Lasso estimation error under the $L_{2}$-norm is imposed
additionally

$\mathbb{P}\left(\sqrt{\sum_{j=1}^{2p}(\widehat{\beta}_{j}-\beta_{j}^{*})^{2}}\leq
c_{0}\sqrt{\\#S^{*}}\lambda_{n}\right)\geq 1-k_{2n}.$ (16)

See [19] for details on deriving Condition 7 under Condition 5 and (16). Note
that the error bounds in Condition 5 and (16) are the standard Lasso error
bounds and have been derived in, for example, [25, 39, 6, 29, 38, 1], in time
series settings, where $k_{1n},k_{2n},k_{3n}$ decrease generally with sample
size $n$ at some polynomial rates. For technical simplicity in the power
analysis, we assume that there are no ties in the magnitude of nonzero
knockoff statistics and the Lasso solutions, and that $\max_{1\leq k\leq
q+1}\\{T,T^{k}\\}<\infty$ in (3) and (10) almost surely.

######  Theorem 3.

Assume that $S^{*}\subset\\{1,\cdots,p\\}$ with $\\#S^{*}>0$.

  * 1)

Let $\widehat{S}$ be returned by the TSKI with Algorithm 2 with constant
tuning parameters $0<\tau_{1}\leq\tau_{*}<1$. Then under Conditions 5–7, it
holds that for all large $n$,

$\mathbb{P}\left(\frac{\\#(S^{*}\cap\widehat{S})}{\\#S^{*}}\geq
1-4c_{0}k_{1n}^{-1}\right)\geq 1-k_{2n}-k_{3n}.$ (17)

  * 2)

Let $\widehat{S}$ be returned by the TSKI with Algorithm 1 with
$\tau_{1}\in(0,1)$ and $\tau_{*}\in(0,1)$. Assume that Conditions 6–7 are
satisfied and Condition 5 holds for the Lasso estimates applied to each
subsample in $H_{k}$ in Algorithm 1. Then it holds that for all large $n$,

$\mathbb{P}\left(\\{\widehat{S}=\emptyset\\}\cup\left\\{\frac{\\#(S^{*}\cap\widehat{S})}{\\#S^{*}}\geq
1-4c_{0}(1+q)k_{1n}^{-1}\right\\}\right)\geq 1-(q+1)\times(k_{2n}+k_{3n}).$
(18)

If we further assume that $\tau_{1}=\tau^{*}\times
K^{-1}\times(1-4(q+1)c_{0}k_{1n}^{-1})$, then we have that for each $K>1$ and
all large $n$,

$\mathbb{E}\left(\frac{\\#(S^{*}\cap\widehat{S})}{\\#S^{*}}\right)\geq\left(1-\frac{(q+1)(\tau_{1}+\theta_{\varepsilon})K}{K-1}-(q+1)\times(k_{2n}+k_{3n})\right)\times
k_{4n},$ (19)

where $\\{k_{4n}\\}$ is some increasing sequence with
$\lim_{n\rightarrow\infty}k_{4n}=1$ and

$\begin{split}\theta_{\varepsilon}=\inf\left\\{\theta\geq 0:\max_{1\leq k\leq
q+1}\mathbb{E}\left(\frac{\\#(\\{j:W_{j}^{k}\geq
T^{k}\\}\cap(S^{*})^{c})}{\\#\\{j:W_{j}^{k}\geq T^{k}\\}\vee
1}\right)\leq\tau_{1}+\theta\right\\}.\end{split}$ (20)

The proof of Theorem 3 is provided in Section A.3 of the Supplementary
Material. As discussed previously, $k_{2n}$ and $k_{3n}$ converge generally to
zero at some polynomial rates of $n$. If we set $q=O((\log{n})^{1+\delta})$
for some $\delta>0$ as discussed earlier, then the second term in the
parentheses on the right-hand side of (18) is asymptotically negligible. In
addition, we need $(\tau_{1}+\theta_{\varepsilon})q$ to be asymptotically
negligible for making the probability of $\\{\widehat{S}=\emptyset\\}$ vanish
asymptotically, which further ensures the high selection power, in view of
(18)–(19). The intuition is that when $(\tau_{1}+\theta_{\varepsilon})q$ is
sufficiently small, the false discovery proportion for each knockoff filter is
well controlled, yielding an upper bound on
$\sum_{s=1}^{p}\mathbf{1}_{\\{W_{s}^{k}\leq-T^{k}\\}}$ for each
$k=1,\cdots,q+1$. This further entails that the $j$th e-value statistic given
in Algorithm 1 is sufficiently large for each $j\in S^{*}$ so that it can be
selected by the e-value procedure and hence $\widehat{S}\not=\emptyset$. In
implementation, in light of the definition of $\tau_{1}$, we can make
$\tau_{1}q$ asymptotically negligible by choosing $K\gg
q=O((\log{n})^{1+\delta})$.

##  4 Simulation studies

We now investigate the finite-sample performance of the TSKI procedure with
Algorithm 1 and Algorithm 2. Two selection methods considered in this section
are TSKI with LCD and TSKI with random forests mean decrease accuracy (MDA),
which are abbreviated as TSKI-LCD and TSKI-MDA, respectively. The theoretical
foundation of the TSKI inference with $q=0$ and $q>0$ for time series
application has been established in Theorem 2 and Theorem 1, respectively.

###  4.1 Simulation settings

We begin with introducing the data generating processes for three simulation
examples. Specifically, we will consider below three popular time series
models: the autoregressive model with exogenous variables (ARX), the self-
exciting threshold autoregressive model with exogenous variables (SETARX)
[36], and the autoregressive conditional heteroskedasticity model with
exogenous variables (ARX-ARCH) [17]. In Models 1–3 below, $\beta$ denotes the
magnitude of the autoregressive coefficient, $H_{t,j}$’s represent the time
series covariates, and $\varepsilon_{t}$’s are model errors.

######  Model 1 (ARX model).

For each integer $t$, we define

$Y_{t}=\sum_{j=1}^{2}(-0.5)^{j-1}\beta Y_{t-j}+\sum_{j=1}^{15}0.6\times
H_{t,j}+\varepsilon_{t}.$

######  Model 2 (SETARX model).

For each integer $t$, we define

$Y_{t}=\begin{cases}\sum_{j=1}^{2}(-0.5)^{j-1}\beta
Y_{t-j}+\sum_{j=1}^{15}0.6\times H_{t,j}+\varepsilon_{t}\quad\textnormal{ if
}Y_{t-d}>r,\\\ \sum_{j=1}^{2}-(-0.5)^{j-1}\beta
Y_{t-j}+\sum_{j=1}^{15}0.6\times H_{t,j}+\varepsilon_{t}\quad\textnormal{
otherwise,}\end{cases}$

where we choose $r=0.5$ and $d=1$ as the threshold value and the threshold
lag, respectively.

######  Model 3 (AR-X-ARCH model).

For each $t$, we define

$Y_{t}=\sum_{j=1}^{2}(-0.5)^{j-1}\beta Y_{t-j}+\sum_{j=1}^{15}0.6\times
H_{t,j}+\sigma_{t}\varepsilon_{t}$

with $\sigma_{t}^{2}=0.1+0.9(\sigma_{t-1}\varepsilon_{t-1})^{2}$.

The set of model errors $\\{\varepsilon_{t}\\}$ is an independent sequence of
i.i.d. zero-mean Gaussian random variables with
$\mathbb{E}\varepsilon_{1}^{2}=1$, and the magnitude of the autoregressive
coefficient is such that $\beta\in\\{0.3,0.7\\}$. The time series covariates
are given by $H_{t,j}=0.2\times H_{t-1,j}+\epsilon_{t,j}$ with
$j\in\\{1,\cdots,50\\}$, where $(\epsilon_{t,1},\cdots,\epsilon_{t,50})$’s are
i.i.d. Gaussian random vectors with zero mean and
$\mathbb{E}(\epsilon_{t,k}\epsilon_{t,l})=0.2^{|k-l|}$ for all $k,l$. We
formulate the coveriate vector with respect to response $Y_{t}$ as
$\boldsymbol{x}_{t}=(Y_{t-1},\cdots,Y_{t-20},\boldsymbol{h}_{t},\boldsymbol{h}_{t-1},\boldsymbol{h}_{t-2},\boldsymbol{h}_{t-3},\boldsymbol{h}_{t-4})$,
where $\boldsymbol{h}_{t}=(H_{t,1},\cdots,H_{t,50})$ and $p=270$.

In Models 1 and 3, the mean functions both depend linearly on the covariates,
while in Model 2 the mean function is a piecewise linear function of the
covariates. For Models 1–2, the relevant index set according to Definition 1
is $S_{0}=\\{1,2,21,\cdots,35\\}$, giving rise to the set of null variables
$\mathcal{H}_{0}=\\{3,\cdots,20,36,\cdots,270\\}$. For Model 3, the relevant
set is $S_{\textnormal{arch}}=\\{1,2,3,21,\cdots,35,71,\cdots 85\\}$ and the
null set is
$\mathcal{H}_{\textnormal{arch}}=\\{4,\cdots,20,36,\cdots,70,86,\cdots,270\\}$,
while $S_{0}$ and $\mathcal{H}_{0}$ defined before are the sets of active and
null covariates, respectively, in the mean regression function. Although
$S_{0}$ and $\mathcal{H}_{0}$ differ from $S_{\textnormal{arch}}$ and
$\mathcal{H}_{\textnormal{arch}}$, respectively, in Model 3, we examine the
empirical power and FDR of the TSKI with respect to $\mathcal{S}_{0}$ and
$\mathcal{H}_{0}$ for two reasons: 1) this is an interesting problem in time
series inference and 2) random forests and Lasso are both algorithms designed
for fitting the mean regression and thus are not expected to detect variables
that affect only the variance function.

For implementation, the target FDR level is set as $\tau^{*}=0.2$, and the R
packages glmnet and randomForest are used for calculating the Lasso estimates
and the random forests MDA, respectively. We generate the approximate knockoff
variables using the idea of the second-order approximation where we pretend
that $\boldsymbol{x}_{t}$ had a multivariate Gaussian distribution, and then
use equation (3) of [19] to simulate the knockoff variables from a conditional
Gaussian distribution with the precision matrix estimated using the method
developed in [21]. Note that such knockoff variables have the first two
moments approximately match the those of the original covariates and thus,
this procedure is named as the second-order approximation. The TSKI Algorithm
1 with subsampling parameter $q\in\\{0,1,2\\}$,
$\tau_{1}\in\\{0.066,0.1,0.2\\}$, and $\tau^{*}=0.2$ is considered in our
simulation. The R code for the simulation experiments is available in the
online Supplementary Material.

Let $\widehat{S}$ be the set of selected features returned by the TSKI-LCD or
TSKI-MDA. Then the empirical versions of the FDR and power are given by

$\begin{split}\textnormal{FDP}\coloneqq\frac{\\#(\widehat{S}\cap\mathcal{H}_{0})}{\max\\{\\#\widehat{S},1\\}}\
\text{ and }\ \textnormal{Sample power}&\coloneqq\frac{\\#(\widehat{S}\cap
S_{0})}{\\#S_{0}},\end{split}$

respectively, where the FDP stands for the false discovery proportion. The
values of $(\beta,q,n)$ are included in both Tables 1–2 and we choose
$\tau_{1}=\tau^{*}/(q+1)$.

###  4.2 Empirical performance of TSKI

$n/\beta$ | Method | $q$ | FDR | Power
---|---|---|---|---
$200$/0.7 | TSKI-LCD | 0 | 0.237 | 0.992
TSKI-LCD | 1 | 0.108 | 0.529
TSKI-LCD | 2 | 0.013 | 0.016
TSKI-MDA | 0 | 0.273 | 0.391
TSKI-MDA | 1 | 0.026 | 0.021
TSKI-MDA | 2 | 0 | 0
$500$/0.7 | TSKI-LCD | 0 | 0.188 | 1
TSKI-LCD | 1 | 0.142 | 0.999
TSKI-LCD | 2 | 0.121 | 0.999
TSKI-MDA | 0 | 0.208 | 0.729
TSKI-MDA | 1 | 0.068 | 0.164
TSKI-MDA | 2 | 0.009 | 0.011

$n/\beta$ | Method | $q$ | FDR | Power
---|---|---|---|---
$200$/0.7 | TSKI-LCD | 0 | 0.203 | 0.985
TSKI-LCD | 1 | 0.122 | 0.755
TSKI-LCD | 2 | 0.027 | 0.057
TSKI-MDA | 0 | 0.220 | 0.292
TSKI-MDA | 1 | 0.044 | 0.021
TSKI-MDA | 2 | 0 | 0
$500$/0.7 | TSKI-LCD | 0 | 0.181 | 0.999
TSKI-LCD | 1 | 0.166 | 0.996
TSKI-LCD | 2 | 0.141 | 0.996
TSKI-MDA | 0 | 0.142 | 0.476
TSKI-MDA | 1 | 0.037 | 0.076
TSKI-MDA | 2 | 0.006 | 0.005

Table 1: The simulation results on the empirical FDR and power for the TSKI
with $\tau_{1}=\tau^{*}/(q+1)$ and $\tau^{*}=0.2$ under Model 1 (left panel)
and Model 3 (right panel) in Section 4.1.

$n/\beta$ | Method | $q$ | FDR | Power
---|---|---|---|---
$200$/0.3 | TSKI-LCD | 0 | 0.168 | 0.874
TSKI-LCD | 1 | 0.076 | 0.411
TSKI-LCD | 2 | 0.005 | 0.004
TSKI-MDA | 0 | 0.171 | 0.584
TSKI-MDA | 1 | 0.047 | 0.056
TSKI-MDA | 2 | 0 | 0
$500$/0.3 | TSKI-LCD | 0 | 0.164 | 0.886
TSKI-LCD | 1 | 0.104 | 0.882
TSKI-LCD | 2 | 0.096 | 0.835
TSKI-MDA | 0 | 0.142 | 0.893
TSKI-MDA | 1 | 0.112 | 0.762
TSKI-MDA | 2 | 0.011 | 0.036

$n/\beta$ | Method | $q$ | FDR | Power
---|---|---|---|---
$200$/0.7 | TSKI-LCD | 0 | 0.150 | 0.685
TSKI-LCD | 1 | 0.043 | 0.082
TSKI-LCD | 2 | 0 | 0
TSKI-MDA | 0 | 0.183 | 0.443
TSKI-MDA | 1 | 0.021 | 0.026
TSKI-MDA | 2 | 0 | 0
$500$/0.7 | TSKI-LCD | 0 | 0.165 | 0.914
TSKI-LCD | 1 | 0.105 | 0.884
TSKI-LCD | 2 | 0.041 | 0.194
TSKI-MDA | 0 | 0.165 | 0.752
TSKI-MDA | 1 | 0.114 | 0.602
TSKI-MDA | 2 | 0.015 | 0.041

Table 2: The simulation results on the empirical FDR and power for the TSKI
with $\tau_{1}=\tau^{*}/(q+1)$ and $\tau^{*}=0.2$ under Model 2 in Section
4.1.

For all simulation experiments reported in Tables 1–2, both TSKI-LCD and TSKI-
MDA with $q>0$ control the FDR in finite samples at the target value of
$\tau^{*}=0.2$, but at the cost of lower selection power compared to the case
of $q=0$. In contrast, TSKI with $q=0$ (i.e., no subsampling) has FDR
exceeding the target FDR level in Model 1 when the sample size is small (i.e.,
$n=200$). When $q=2$ and the sample size is small, TSKI becomes overly
conservative with both low FDR and low power in most cases, with the MDA based
method suffering more severely.

In all but one settings (Table 2, $(n,\beta)=(500,0.3)$, the LCD based method
outperforms the MDA based method in terms of power. Three plausible
explanations for this are 1) the mean regression functions in all three models
are not too far away from linear functions, 2) MDA is a nonparameteric method
which requires larger sample size to perform well compared to LCD, and 3) MDA
suffers from the potential bias issue due to correlations between covariates
[34, 13].

Overall, our simulation experiments show that the time series FDR is
controlled stably by the TSKI procedure with subsampling parameter $q>0$ in
finite samples, whereas the selection power depends on the sample size, the
subsampling parameter $q\geq 0$, the signal strength $\beta$, and the choice
of the feature importance measure. Regarding the subsampling parameter $q$, we
remark that by the construction of (3) in Algorithm 1, TSKI may have decent
asymptotic power only when the number of relevant features is no less than
$\tau_{1}^{-1}=(q+1)/\tau^{*}$. To see the intuition, let us consider a case
when the number of relevant features is less than $\tau_{1}^{-1}$ and their
knockoff statistics are the few largest (positive) statistics asymptotically.
Then even if the other knockoff statistics are all zero, we have
$T^{k}=\infty$ in (3), and hence the knockoff filter screens out all features.
It is worth mentioning that a similar but more restrictive requirement is
assumed in Condition 7 for the asymptotic power.

##  5 Real data application

Identifying important time series economic variables that can affect the
inflation has been an active research problem with a long history [24, 15],
due to the importance of inflation. In this section, we will analyze the
temporal relations between the U.S. inflation and other major time series
covariates of the U.S. economy from May 2013 ($t=1$) to January 2023 ($t=117$)
using monthly data in this period. From Figure 1, we see that the inflation
series consists of two nonstationary points at January 2015 and April 2020,
corresponding to events A and B, respectively. Event A corresponds to a sharp
gasoline index decline, whereas event B is due to the onset of the COVID-19
pandemic. These events cannot be predicted based on the available information
before their occurrences. In addition, the U.S. economy has been experiencing
economic recovery in the post-COVID-19 era, and there seem to be effects of
Russia–Ukraine war, as revealed by the clear structural change during time
period C indicated in Figure 1. In light of these versatile time varying
patterns of the inflation series, we facilitate our analysis by considering
five-year rolling windows, where the ends of the rolling windows start from
April 2018 to January 2023 (a total of $58$ rolling windows).

Figure 1: The U.S. inflation from May 2013 to January 2023.

Figure 2:  The left panel displays the averages of “having any selections”
indicators over 100 repetitions, where the indicator at each rolling window is
one if and only if any covariates are selected, and the $x$-axis indicates the
ending time of each rolling window. The right panel shows the analogous
results, but with the indicator being one at each covariate index if that
covariate is selected at any rolling windows. The first 127 covariates are
current time covariates, and the $128$th to $254$th covariates are one month
lag covariates in the AR(2) model. Covariates measuring similar economic
values are clustered closer (see [28] for detailed definitions of these
covariates). The selection method here is the TSKI-LCD without subsampling
($q=0$).

Figure 3:  These two figures are analogous to Figure 2 but with $q=1$ for the
TSKI-LCD procedure.

(a) The red curve is ACOGNO

(b) The red curve is EXCAUSx

(c) The red curve is CLAIMSx

Figure 4: The black curves in the three panels are the inflation series at
time $t+1$. The red curve in panel (a) is the number of new orders for
consumer goods at time $t$, the red curve in panel (b) indicates the
U.S./Canada exchange rate at time $t$, and the red curve in panel (c) is the
U.S. initial claims for unemployment benefits at time $t$. All curves here are
standardized and adjusted for visual comparison, and hence the values of these
time series are not reported on the $y$-axis.

Specifically, we apply the TSKI-LCD with $q\in\\{0,1\\}$ to each five-year
rolling window to investigate the temporal relations between the inflation and
other time series lags. These monthly economic time series, including numerous
types of consumer price indices, unemployment rates, and housing prices, can
be obtained from the FRED-MD database [28]333The website URL:
https://research.stlouisfed.org/econ/mccracken/fred-databases/. and the U.S.
Bureau of Labor Statistics. These time series have been adjusted according to
the instruction of the FRED-MD database, which can be found in the work of
[28] or the codes provided in the paper. In particular, the inflation at time
$t$ is defined as the adjusted consumer price index for all goods

$\textnormal{Inflation}_{t}\coloneqq\left(\frac{\textnormal{CPI}_{t}-\textnormal{CPI}_{t-1}}{\textnormal{CPI}_{t-1}}\times
100\right)\%,$

where CPIt is the consumer price index for all goods at time (month) $t$. Each
time series has a FRED-MD code. For example, CPIAUCSL is the FRED-MD code of
the inflation series. To reduce randomness resulting from the use of the
knockoffs construction, we repeat the inference procedure $100$ times, and
report the average results in Figure 2 with $q=0$, where the left panel
displays whether any significant covaraites have been found at each rolling
window period, and the right panel shows the selection frequency over $58$
rolling window periods of the 127 time series covariates and their one month
lags. That is, we model the one month ahead inflation series
$\text{Inflation}_{t+1}$ (response) using an AR(2) model with $127$ time
series covariates at current time $t$ and $127$ one month lags of these time
series covariates. So the covariate dimensionality is $p=254$.

As can be seen in Figure 2, the TSKI-LCD with $q=0$ identifies some active
windows around the B and C periods, and the selection frequency concentrates
on a sparse set of covariates. The majority of the selected variables are
covariates at current time $t$. The top $10$ most frequently selected
covariates by the TSKI-LCD with $q=0$ are employment related series (HWI,
CLAIMSx), consumption related indices (ACOGNO, CPIAUCSL), housing related
series (PERMITS), U.S. bond yields (GS5, GS10), stock market indices
(S.P.div.yield, S.P.500), and exchange rates (EXCAUSx) at their current time
$t$, where their FRED-MD codes are given in the parentheses.

The simulation results in Section 4 suggest that the choice of $q=1$ has
better FDR control especially when the sample size is limited. Motivated by
our simulation results, we next apply the TSKI-LCD with subasmpling $q=1$ in
Figure 3 with the expectation of better FDR control. The results of Figure 3
are more conservative in comparison to those in Figure 2. Despite being
conservative, these new results also suggest some active windows in the B and
C periods, and the selected variables are also mostly covariates at current
time $t$. In addition, most covariates selected by the TSKI-LCD with $q=1$
belong to the set selected by the TSKI-LCD with $q=0$. In particular, the top
10 selected covariates are CPIAUCSL, CLAIMSx, PERMITS, AMDMUOx, GS10, ACOGNO,
EXCAUSx, EXUSUKx, INDPRO, IPFUELS, where only PERMITS is one month lag
covariate at $t-1$ in the AR(2) model. Among them, EXUSUKx, AMDMUOx, INDPRO,
and IPFUELS are new in comparison to the list of selected set when $q=0$,
where the first two are the U.S./U.K. exchange rate and another type of
consumption related index (the number of unfilled orders for durable goods),
respectively, and the last two are industrial production indices that are
related to consumption price indices. The difference in the selected sets of
covariates are attributed to the fact that some economic covariates are
designed to track similar economic factors and tend to be highly correlated.

The recent literature [33] suggests that the inflation foresting is a
difficult task in the sense that AR models with additional time series
covariates rarely outperform simple AR models with only inflation lags. In
other words, conditional on the lagged inflation series, additional covariates
do not carry strong signals in the inflation forecasting. The fact that the
TSKI-LCD with $q=1$ selects only a few time series covariates indeed supports
such argument. It is also interesting to notice that the stock market indices
are not considered as important covariates by the TSKI-LCD with subsampling
(i.e., $q=1$), but are selected as important covariates when $q=0$. In
particular, the selection frequencies of S&P dividend yields and S&P 500 (both
at lag $t-1$) are $100\%$ and $79\%$, respectively, in Figure 2, while only
$5\%$ and $1\%$, respectively, in Figure 3, suggesting that stock market
indices could be spurious findings.

The selection results by the TSKI-LCD motivate us to further investigate the
dependency of inflation on a few time series covariates. In Figure 4, we
further plot three selected series, namely ACOGNO, EXCAUSx, and CLAIMSx, that
are among the top 10 lists both when $q=0$ and $q=1$. These three selected
covariates are all at their current time $t$ in the AR(2) model. ACOGNO is the
number of new orders for consumer goods which is an important consumption
index, EXCAUSx is the exchange rate from U.S. dollar to Canadian dollar, and
CLAIMSx is the initial claim for unemployment benefits. We also include the
inflation series at time $t+1$ in the same period (the black curve in each
panel). For better visual comparison, all curves in Figure 4 have been
standardized and adjusted.

A visual inspection of Figure 4 shows that the impacts of the COVID-19
pandemic in April 2020 on the U.S. economy are stronger than those of gasoline
price shock in January 2015. This potentially explains why our empirical
findings of significant covariates concentrate mostly around periods B and C.
From Figure 4, we see that the gasoline price shock in January 2015 affects
the consumption index series ACOGNO more mildly compared to the impact of
COVID-19 in April 2020. Also, although there is some variation in the exchange
rate EXCAUSx after January 2015, it is unclear whether such variation was
caused by the gasoline price shock. In contrast, most time series of the U.S.
economy responded clearly to the pandemic to an unigonorable degree.
Particularly, the exchange rate and the number of initial claims dropped in
March 2020, suggesting that these covariates are one-month ahead indicators
for the inflation drop in April 2020. In summary, we have applied the newly
suggested tool of TSKI to study the U.S. economy. Our empirical results
illustrate the potential of the TSKI on obtaining more instructive findings in
real data applications.

## References

  * Adamek et al. [2022] Adamek, R., S. Smeekes, and I. Wilms (2022).  Lasso inference for high-dimensional time series.  Journal of Econometrics, to appear.
  * An and Huang [1996] An, H. and F. Huang (1996).  The geometrical ergodicity of nonlinear autoregressive models.  Statistica Sinica 6, 943–956.
  * Bai et al. [2021] Bai, X., J. Ren, Y. Fan, and F. Sun (2021).  KIMI: knockoff inference for motif identification from molecular sequences with controlled false discovery rate.  Bioinformatics 37, 759–766.
  * Barber and Candès [2015] Barber, R. F. and E. J. Candès (2015).  Controlling the false discovery rate via knockoffs.  The Annals of Statistics 43, 2055–2085.
  * Barber et al. [2020] Barber, R. F., E. J. Candès, and R. J. Samworth (2020).  Robust inference with knockoffs.  The Annals of Statistics 48, 1409–1431.
  * Basu and Michailidis [2015] Basu, S. and G. Michailidis (2015).  Regularized estimation in sparse high-dimensional time series models.  The Annals of Statistics 43, 1535–1567.
  * Benjamini and Hochberg [1995] Benjamini, Y. and Y. Hochberg (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal Statistical Society Series B 57, 289–300.
  * Benjamini and Yekutieli [2001] Benjamini, Y. and D. Yekutieli (2001).  The control of the false discovery rate in multiple testing under dependency.  The Annals of statistics 29, 1165–1188.
  * Bollerslev [1986] Bollerslev, T. (1986).  Generalized autoregressive conditional heteroskedasticity.  Journal of Econometrics 31, 307–327.
  * Brockwell et al. [1991] Brockwell, P. J., R. A. Davis, and S. E. Fienberg (1991).  Time Series: Theory and Methods.  Springer Science & Business Media.
  * Candès et al. [2018] Candès, E., Y. Fan, L. Janson, and J. Lv (2018).  Panning for gold: ‘model-X’ knockoffs for high dimensional controlled variable selection.  Journal of the Royal Statistical Society Series B 80, 551–577.
  * Chen and Tsay [1993] Chen, R. and R. S. Tsay (1993).  Functional-coefficient autoregressive models.  Journal of the American Statistical Association 88(421), 298–308.
  * Chi et al. [2022] Chi, C.-M., Y. Fan, and J. Lv (2022).  FACT: high-dimensional random forests inference.  arXiv preprint arXiv:2207.01678.
  * Cline and Pu [2004] Cline, D. B. H. and H.-M. Pu (2004).  Stability and the Lyapounov exponent of threshold AR-ARCH models.  The Annals of Applied Probability 14(4), 1920–1949.
  * Crump et al. [2022] Crump, R. K., S. Eusepi, M. Giannoni, and A. Şahin (2022).  The unemployment-inflation trade-off revisited: the Phillips curve in COVID times.  Technical report, National Bureau of Economic Research.
  * Durrett [2019] Durrett, R. (2019).  Probability: Theory and Examples, Volume 49.  Cambridge University Press.
  * Engle [1982] Engle, R. F. (1982).  Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation.  Econometrica 50, 987–1007.
  * Fan and Lv [2008] Fan, J. and J. Lv (2008).  Sure independence screening for ultrahigh dimensional feature space (with discussion).  Journal of the Royal Statistical Society Series B 70, 849–911.
  * Fan et al. [2020] Fan, Y., E. Demirkaya, G. Li, and J. Lv (2020).  RANK: large-scale inference with graphical nonlinear knockoffs.  Journal of the American Statistical Association 115, 362–379.
  * Fan et al. [2019] Fan, Y., E. Demirkaya, and J. Lv (2019).  Nonuniformity of p-values can occur early in diverging dimensions.  Journal of Machine Learning Research 20(77), 1–33.
  * Fan and Lv [2016] Fan, Y. and J. Lv (2016).  Innovated scalable efficient estimation in ultra-large Gaussian graphical models.  The Annals of Statistics 44(5), 2098–2126.
  * Fan et al. [2020] Fan, Y., J. Lv, M. Sharifvaghefi, and Y. Uematsu (2020).  IPAD: stable interpretable forecasting with knockoffs inference.  Journal of the American Statistical Association 115, 1822–1834.
  * Hansen [2011] Hansen, B. E. (2011).  Threshold autoregression in economics.  Statistics and Its Interface 4, 123–127.
  * King et al. [1995] King, R. G., J. H. Stock, and M. W. Watson (1995).  Temporal instability of the unemployment-inflation relationship.  Economic Perspectives 19(3), 2–13.
  * Kock and Callot [2015] Kock, A. B. and L. Callot (2015).  Oracle inequalities for high dimensional vector autoregressions.  Journal of Econometrics 186(2), 325–344.
  * Lanne and Saikkonen [2005] Lanne, M. and P. Saikkonen (2005).  Non-linear GARCH models for highly persistent volatility.  The Econometrics Journal 8(2), 251–276.
  * Lu et al. [2018] Lu, Y., Y. Fan, J. Lv, and W. S. Noble (2018).  DeepPINK: reproducible feature selection in deep neural networks.  Advances in Neural Information Processing Systems (NeurIPS 2018).
  * McCracken and Ng [2016] McCracken, M. W. and S. Ng (2016).  FRED-MD: a monthly database for macroeconomic research.  Journal of Business & Economic Statistics 34, 574–589.
  * Medeiros and Mendes [2016] Medeiros, M. C. and E. F. Mendes (2016).  L1-regularization of high-dimensional time-series models with non-Gaussian and heteroskedastic errors.  Journal of Econometrics 191, 255–271.
  * Meitz and Saikkonen [2010] Meitz, M. and P. Saikkonen (2010).  A note on the geometric ergodicity of a nonlinear AR-ARCH model.  Statistics & probability letters 80(7-8), 631–638.
  * Ozaki [1985] Ozaki, T. (1985).  2 non-linear time series models and dynamical systems.  Handbook of Statistics 5, 25–83.
  * Ren and Barber [2022] Ren, Z. and R. F. Barber (2022).  Derandomized knockoffs: leveraging e-values for false discovery rate control.  arXiv preprint arXiv:2205.15461.
  * Stock and Watson [2007] Stock, J. H. and M. W. Watson (2007).  Why has U.S. inflation become harder to forecast?  Journal of Money, Credit and Banking 39, 3–33.
  * Strobl et al. [2007] Strobl, C., A.-L. Boulesteix, A. Zeileis, and T. Hothorn (2007).  Bias in random forest variable importance measures: illustrations, sources and a solution.  BMC Bioinformatics 8(1), 1–21.
  * Tjøstheim [1990] Tjøstheim, D. (1990).  Non-linear time series and Markov chains.  Advances in Applied Probability 22, 587–611.
  * Tong and Lim [1980] Tong, H. and K. S. Lim (1980).  Threshold autoregression, limit cycles and cyclical data.  Journal of the Royal Statistical Society Series B 42, 245–292.
  * Wang and Ramdas [2020] Wang, R. and A. Ramdas (2020).  False discovery rate control with e-values.  arXiv preprint arXiv:2009.02824.
  * Wong et al. [2020] Wong, K. C., Z. Li, and A. Tewari (2020).  Lasso guarantees for $\beta$-mixing heavy-tailed time series.  The Annals of Statistics 48, 1124–1142.
  * Wu and Wu [2016] Wu, W.-B. and Y. N. Wu (2016).  Performance bounds for parameter estimates of high-dimensional linear models with correlated errors.  Electronic Journal of Statistics 10, 352–379.
  * Yu [1994] Yu, B. (1994).  Rates of convergence for empirical processes of stationary mixing sequences.  The Annals of Probability, 94–116.
  * Zhu et al. [2021] Zhu, Z., Y. Fan, Y. Kong, J. Lv, and F. Sun (2021).  DeepLINK: deep learning inference using knockoffs with applications to genomics.  Proceedings of the National Academy of Sciences of the United States of America 118, e2104683118.

Supplementary Material to “High-Dimensional Knockoffs Inference for Time
Series Data”

Chien-Ming Chi, Yingying Fan, Ching-Kang Ing and Jinchi Lv

This Supplementary Material contains the proofs of all main results and
technical lemmas, and some additional technical details. All the notation is
the same as defined in the main body of the paper.

##  Appendix A Proofs of Theorems 1–3, Corollaries 1–2, and Proposition 1

###  A.1 Proof of Theorem 1

For simplicity, in this proof we use the notation
$\boldsymbol{U},\widetilde{\boldsymbol{U}},\boldsymbol{V}$ to denote a generic
subsample of
$\\{\boldsymbol{x}_{i},\widetilde{\boldsymbol{x}}_{i},Y_{i}\\}_{i=1}^{n}$ or
their independent and identically distributed (i.i.d.) counterparts
$\\{\boldsymbol{x}_{i}^{\pi},\widetilde{\boldsymbol{x}}_{i}^{\pi},Y_{i}^{\pi}\\}_{i=1}^{n}$,
where their exact meaning will be made explicit whenever confusion is
possible. Let us define

$\widehat{\textnormal{KL}}_{j}^{k}\equiv\log{\left(\frac{f_{\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V}}(\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V})}{f_{\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V}}(\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V})}\right)},$
(A.1)

where $\boldsymbol{U}=(\boldsymbol{x}_{i},i\in H_{k})^{{}^{T}}$,
$\widetilde{\boldsymbol{U}}=(\widetilde{\boldsymbol{x}}_{i},i\in
H_{k})^{{}^{T}}$, and $\boldsymbol{V}=(Y_{i},i\in H_{k})^{{}^{T}}$. We can use
(A.7) in the proof of Theorem 2 in Section A.2 to conclude that for each
$k\in\\{1,\cdots,q+1\\}$,

$\mathbb{E}(\sum_{j\in\mathcal{H}_{0}}e_{j}^{k}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}^{k}\leq\varepsilon\\}})\leq
p\times e^{\varepsilon}.$

Then it holds that

$\mathbb{E}(\sum_{j\in\mathcal{H}_{0}}e_{j}^{(\varepsilon)})\leq p\times
e^{\varepsilon},$

where

$e_{j}^{(\varepsilon)}=(q+1)^{-1}\sum_{k=1}^{q+1}e_{j}^{k}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}^{k}\leq\varepsilon\\}}.$

Denote by $\widehat{S}_{\varepsilon}$ be the set of selected features when
applying the e-BH method to $e_{j}^{(\varepsilon)}$’s at the target FDR level
$\tau^{*}$. Then using similar arguments to those for (A.8) in Section A.2, we
can show that

$\mathbb{E}\left(\frac{\\#(\widehat{S}_{\varepsilon}\cap\mathcal{H}_{0})}{(\\#\widehat{S}_{\varepsilon})\vee
1}\right)\leq\tau^{*}\times e^{\varepsilon}.$

Thus, it follows that

$\begin{split}&\mathbb{E}\left(\frac{\\#(\widehat{S}\cap\mathcal{H}_{0})}{(\\#\widehat{S})\vee
1}\right)\leq\mathbb{E}\left(\frac{\\#(\widehat{S}_{\varepsilon}\cap\mathcal{H}_{0})}{(\\#\widehat{S}_{\varepsilon})\vee
1}\times\boldsymbol{1}_{\\{\widehat{S}=\widehat{S}_{\varepsilon}\\}}+\boldsymbol{1}_{\\{\widehat{S}\not=\widehat{S}_{\varepsilon}\\}}\right)\\\
&\leq\tau^{*}\times e^{\varepsilon}+\sum_{k=1}^{q+1}\mathbb{P}(\max_{1\leq
j\leq p}\widehat{\textnormal{KL}}_{j}^{k}>\varepsilon)\end{split}$ (A.2)

in view of
$\\{\widehat{S}\not=\widehat{S}_{\varepsilon}\\}\subset\\{\max_{1\leq k\leq
q+1}\max_{1\leq j\leq p}\widehat{\textnormal{KL}}_{j}^{k}>\varepsilon\\}$.
Since (LABEL:theorem1.eq.4) holds for each $\varepsilon>0$, we can further
obtain that

$\begin{split}\mathbb{E}\left(\frac{\\#(\widehat{S}\cap\mathcal{H}_{0})}{(\\#\widehat{S})\vee
1}\right)\leq&\inf_{\varepsilon>0}\left\\{\tau^{*}e^{\varepsilon}+\sum_{k=1}^{q+1}\mathbb{P}(\max_{1\leq
j\leq p}\widehat{\textnormal{KL}}_{j}^{k\pi}>\varepsilon)\right\\}\\\
&+\sum_{k=1}^{q+1}\Big{(}\mathbb{P}(\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}^{k}>\varepsilon)-\mathbb{P}(\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}^{k\pi}>\varepsilon)\Big{)}.\end{split}$ (A.3)

It remains to bound the second term on the right-hand side of (A.3) above.
Recall that $\widehat{\textnormal{KL}}_{j}^{k\pi}$ is defined analogously to
(A.1) but based on i.i.d. sample
$\\{\boldsymbol{x}_{i}^{\pi},\widetilde{\boldsymbol{x}}_{i}^{\pi},Y_{i}^{\pi}\\}_{i=1}^{n}$.
With $\boldsymbol{U}=(\boldsymbol{x}_{i}^{\pi},i\in H_{k})^{{}^{T}}$,
$\widetilde{\boldsymbol{U}}=(\widetilde{\boldsymbol{x}}_{i}^{\pi},i\in
H_{k})^{{}^{T}}$, and $\boldsymbol{V}=(Y_{i}^{\pi},i\in H_{k})^{{}^{T}}$, we
can deduce that

$\begin{split}\widehat{\textnormal{KL}}_{j}^{k\pi}&=\sum_{i\in
H_{k}}\log{\left(\frac{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(X_{ij}^{\pi},\widetilde{X}_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi},\widetilde{\boldsymbol{x}}_{-ij}^{\pi})}{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(\widetilde{X}_{ij}^{\pi},X_{ij}^{\pi},\boldsymbol{x}_{-ij}^{\pi},\widetilde{\boldsymbol{x}}_{-ij}^{\pi})}\right)}\\\
&=\log{\left(\frac{f_{\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j}}(\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j})}{f_{\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j}}(\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j})}\right)}\\\
&=\log{\left(\frac{f_{\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V}}(\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V})}{f_{\boldsymbol{U}_{j},\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V}}(\widetilde{\boldsymbol{U}}_{j},\boldsymbol{U}_{j},\boldsymbol{U}_{-j},\widetilde{\boldsymbol{U}}_{-j},\boldsymbol{V})}\right)},\end{split}$
(A.4)

where the third equality above follows from similar analysis to that of (13)
in Theorem 1. The conditional column independence required by (13) holds for
each $j\in\mathcal{H}_{0}$ because (A.4) involves i.i.d. samples.

We can now see that $\widehat{\textnormal{KL}}_{j}^{k}$ is well-defined thanks
to Condition 1. Moreover, by the fact that the supports of
$(\boldsymbol{x},\widetilde{\boldsymbol{x}})$ and
$[\boldsymbol{x},\widetilde{\boldsymbol{x}}]_{\textnormal{swap}(\\{j\\})}$ are
the same (as guaranteed by Condition 1) and the definition of
$(\boldsymbol{x}_{i}^{\pi},\widetilde{\boldsymbol{x}}_{i}^{\pi})$’s,
$\widehat{\textnormal{KL}}_{j}^{k\pi}$ is well-defined. Hence, in light of
(A.1) and (A.4), there exists some measurable function
$g:\mathbb{R}^{\\#H_{k}\times(2p+1)}\longmapsto\mathbb{R}$ such that
$\widehat{\textnormal{KL}}_{j}^{k}=g(\mathcal{X}_{k})$ and
$\widehat{\textnormal{KL}}_{j}^{k\pi}=g(\mathcal{X}_{k}^{\pi})$ for each
$\varepsilon\geq 0$, which entails that there exists some
$\mathcal{D}\in\mathcal{R}^{\\#H_{k}\times(2p+1)}$ such that

$\begin{split}\\{\mathcal{X}_{k}\in\mathcal{D}\\}&=\\{\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}^{k}>\varepsilon\\},\\\
\\{\mathcal{X}_{k}^{\pi}\in\mathcal{D}\\}&=\\{\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}^{k\pi}>\varepsilon\\}.\end{split}$ (A.5)

With the aid of (A.5), it holds that

$\begin{split}&|\mathbb{P}(\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}^{k}>\varepsilon)-\mathbb{P}(\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}^{k\pi}>\varepsilon)|\\\
&\leq\sup_{\mathcal{D}\in\mathcal{R}^{\\#H_{k}\times(2p+1)}}|\mathbb{P}(\mathcal{X}_{k}\in\mathcal{D})-\mathbb{P}(\mathcal{X}_{k}^{\pi}\in\mathcal{D})|.\end{split}$
(A.6)

Therefore, from (A.3)–(A.4) and (LABEL:theorem1.eq.5) we can obtain the
desired conclusion, which completes the proof of Theorem 1.

###  A.2 Proof of Theorem 2

Let us first make a useful claim that with $\widehat{\textnormal{KL}}_{j}$’s
given in (12), it holds that for each $\varepsilon>0$,

$\sum_{j\in\mathcal{H}_{0}}\mathbb{E}(e_{j}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}})\leq
p\times e^{\varepsilon}.$ (A.7)

Then we consider an application of the e-BH method [37] to
$e_{j}^{(\varepsilon)}\coloneqq
e_{j}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}}$
with the target FDR level $\tau^{*}$, yielding a set of selected features
$\widehat{S}_{\varepsilon}\subset\\{1,\cdots,p\\}$ defined as

$\widehat{S}_{\varepsilon}=\\{j:e_{j}^{(\varepsilon)}\geq
p(\tau^{*}\times\widehat{k}_{\varepsilon})^{-1}\\}$

with $\widehat{k}_{\varepsilon}\equiv\max\\{k:e_{(k)}^{(\varepsilon)}\geq
p(\tau^{*}\times k)^{{-1}}\\}$. Here, $e_{(j)}^{(\varepsilon)}$’s are the
ordered statistics of $e_{j}^{(\varepsilon)}$’s such that
$e_{(1)}^{(\varepsilon)}\geq\cdots\geq e_{(p)}^{(\varepsilon)}$. It is easy to
see that $\\#\widehat{S}_{\varepsilon}=\widehat{k}_{\varepsilon}$.

In view of the definition of $\widehat{k}_{\varepsilon}$, we can deduce that

$\begin{split}\mathbb{E}\left(\frac{\\#(\widehat{S}_{\varepsilon}\cap\mathcal{H}_{0})}{(\\#\widehat{S}_{\varepsilon})\vee
1}\right)&=\mathbb{E}\left(\frac{\sum_{j\in\mathcal{H}_{0}}\boldsymbol{1}_{\\{j\in\widehat{S}_{\varepsilon}\\}}}{\widehat{k}_{\varepsilon}\vee
1}\right)\\\
&\leq\mathbb{E}\left(\frac{\sum_{j\in\mathcal{H}_{0}}\boldsymbol{1}_{\\{j\in\widehat{S}_{\varepsilon}\\}}\times\tau^{*}\times
e_{j}^{(\varepsilon)}}{p}\right)\\\
&\leq\tau^{*}p^{-1}\times\mathbb{E}\left(\sum_{j\in\mathcal{H}_{0}}e_{j}^{(\varepsilon)}\right)\\\
&\leq\tau^{*}e^{\varepsilon},\end{split}$ (A.8)

where the last inequality above is from (A.7). Then combining (A.8) and
$\\{\widehat{S}\not=\widehat{S}_{\varepsilon}\\}\subset\cup_{j=1}^{p}\\{\widehat{\textnormal{KL}}_{j}>\varepsilon\\}$
leads to

$\begin{split}\mathbb{E}\left(\frac{\\#(\widehat{S}\cap\mathcal{H}_{0})}{(\\#\widehat{S})\vee
1}\right)&\leq\mathbb{E}\left(\boldsymbol{1}_{\\{\widehat{S}=\widehat{S}_{\varepsilon}\\}}\times\frac{\\#(\widehat{S}_{\varepsilon}\cap\mathcal{H}_{0})}{(\\#\widehat{S}_{\varepsilon})\vee
1}+\boldsymbol{1}_{\\{\widehat{S}\not=\widehat{S}_{\varepsilon}\\}}\right)\\\
&\leq\tau^{*}e^{\varepsilon}+\mathbb{P}(\max_{1\leq j\leq
p}\widehat{\textnormal{KL}}_{j}>\varepsilon)\end{split}$

for each $\varepsilon\geq 0$. This concludes the proof for the desired result
(11). We will provide the proofs of (A.7) and (13)–(15), separately.

Proof of (A.7). Let us define

$T_{j}\equiv\min\left\\{t\in\mathcal{W}_{+}^{\dagger}:\frac{1+\\#\\{s:W_{s}^{\dagger}\leq-t\\}}{\\#\\{s:W_{s}^{\dagger}\geq
t\\}\vee 1}\leq\tau_{1}\right\\},$

where $W_{k}^{\dagger}=W_{k}$ if $k\not=j$ and $W_{j}^{\dagger}=|W_{j}|$,
$\mathcal{W}_{+}^{\dagger}=\\{|W_{s}^{\dagger}|:|W_{s}^{\dagger}|>0\\}$, and
$\min\emptyset$ is defined as infinity. We further define
$\boldsymbol{X}_{j}^{(0)}$ and $\boldsymbol{X}_{j}^{(1)}$ such that
$\boldsymbol{X}_{j}^{(0)}=\boldsymbol{X}_{j}$ and
$\boldsymbol{X}_{j}^{(1)}=\widetilde{\boldsymbol{X}}_{j}$ if $W_{j}\geq 0$,
and $\boldsymbol{X}_{j}^{(1)}=\boldsymbol{X}_{j}$ and
$\boldsymbol{X}_{j}^{(0)}=\widetilde{\boldsymbol{X}}_{j}$ if $W_{j}<0$.

For each $\varepsilon\geq 0$, we can deduce that

$\begin{split}&\sum_{j\in\mathcal{H}_{0}}\mathbb{E}\left(\frac{\boldsymbol{1}_{\\{W_{j}\geq
T\\}}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}}}{1+\sum_{s=1}^{p}\boldsymbol{1}_{\\{W_{s}\leq-T\\}}}\right)=\sum_{j\in\mathcal{H}_{0}}\mathbb{E}\left(\frac{\boldsymbol{1}_{\\{W_{j}\geq
T_{j}\\}}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}}}{1+\sum_{s=1}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{j}\\}}}\right)\\\
&=\sum_{j\in\mathcal{H}_{0}}\mathbb{E}\left(\frac{\boldsymbol{1}_{\\{W_{j}\geq
T_{j}\\}}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}}}{1+\sum_{s=1,s\not=j}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{j}\\}}}\right)\\\
&=\sum_{j\in\mathcal{H}_{0}}\mathbb{E}\left(\frac{\boldsymbol{1}_{\\{W_{j}>0\\}}\times\boldsymbol{1}_{\\{|W_{j}|\geq
T_{j}\\}}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}}}{1+\sum_{s=1,s\not=j}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{j}\\}}}\right)\\\
&=\sum_{j\in\mathcal{H}_{0}}\mathbb{E}\left(\frac{\mathbb{P}(W_{j}>0,\widehat{\textnormal{KL}}_{j}\leq\varepsilon\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\times\boldsymbol{1}_{\\{|W_{j}|\geq
T_{j}\\}}}{1+\sum_{s=1,s\not=j}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{j}\\}}}\right),\end{split}$ (A.9)

where the first three equalities above hold because when
$\boldsymbol{1}_{\\{W_{j}\geq T\\}}=1$, we have $W_{j}>0$, $T=T_{j}$, and
$\boldsymbol{1}_{\\{W_{j}\leq-T_{j}\\}}=0$, and the last equality above holds
since $|W_{j}|$, $T_{j}$, and $W_{1},\cdots,W_{j-1}$, $W_{j+1},\cdots,W_{p}$
are functions of
$(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})$
due to the sign-flip property (2).

From the definitions of $\boldsymbol{X}_{j}^{(0)}$,
$\boldsymbol{X}_{j}^{(1)}$, and $\widehat{\textnormal{KL}}_{j}$, we can obtain
that

$\begin{split}&\mathbb{P}(W_{j}>0,\widehat{\textnormal{KL}}_{j}\leq\varepsilon\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=\mathbb{P}(W_{j}>0,\widehat{\textnormal{KL}}_{j}^{(01)}\leq\varepsilon\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=\mathbb{P}(W_{j}>0\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}^{(01)}\leq\varepsilon\\}},\end{split}$
(A.10)

where

$\widehat{\textnormal{KL}}_{j}^{(01)}\equiv\log{\left(\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}\right)}.$

Furthermore, we will show that it holds almost surely that

$\begin{split}&\mathbb{P}(W_{j}>0\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&\leq\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}\\\
&\qquad\times\mathbb{P}(W_{j}<0\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=e^{\widehat{\textnormal{KL}}_{j}^{(01)}}\times\mathbb{P}(W_{j}<0\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}),\end{split}$
(A.11)

where Condition 1 is assumed to avoid division by zero on the right-hand side
(RHS) of the second inequality above. The proof of (LABEL:robust.eq.2) is
deferred to right after the proof of (A.7).

By (LABEL:robust.eq.18)–(LABEL:robust.eq.2), it holds that

$\begin{split}&\textnormal{RHS of \eqref{robust.eq.15}}\\\
&\leq\sum_{j\in\mathcal{H}_{0}}\mathbb{E}\left(\frac{e^{\widehat{\textnormal{KL}}_{j}^{(01)}}\mathbb{P}(W_{j}<0\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\,\boldsymbol{1}_{\\{|W_{j}|\geq
T_{j}\\}}\,\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}^{(01)}\leq\varepsilon\\}}}{1+\sum_{s=1,s\not=j}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{j}\\}}}\right)\\\ &\leq
e^{\varepsilon}\,\sum_{j=1}^{p}\mathbb{E}\left(\frac{\boldsymbol{1}_{\\{W_{j}\leq-
T_{j}\\}}}{1+\sum_{s=1,s\not=j}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{j}\\}}}\right)\\\ &\leq
e^{\varepsilon}\,\mathbb{E}\left(\frac{\sum_{j=1}^{p}\boldsymbol{1}_{\\{W_{j}\leq-
T_{j}\\}}}{1+\sum_{s=1,s\not=j}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{s}\\}}}\right)\\\ &\leq
e^{\varepsilon}\,\mathbb{E}\left(\frac{\sum_{j=1}^{p}\boldsymbol{1}_{\\{W_{j}\leq-
T_{j}\\}}}{1\vee\left(\sum_{s=1}^{p}\boldsymbol{1}_{\\{W_{s}\leq-
T_{s}\\}}\right)}\right)\\\ &\leq e^{\varepsilon},\end{split}$ (A.12)

where RHS is short for the right-hand side and the third inequality above
follows from Lemma 6 in [5]. Hence, by resorting to (LABEL:robust.eq.15),
(LABEL:robust.eq.23), and the fact that

$\sum_{j\in\mathcal{H}_{0}}\mathbb{E}(e_{j}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}})=p\times\sum_{j\in\mathcal{H}_{0}}\mathbb{E}\left(\frac{\boldsymbol{1}_{\\{W_{j}\geq
T\\}}\times\boldsymbol{1}_{\\{\widehat{\textnormal{KL}}_{j}\leq\varepsilon\\}}}{1+\sum_{s=1}^{p}\boldsymbol{1}_{\\{W_{s}\leq-T\\}}}\right),$

we can establish (A.7).

Proof of (LABEL:robust.eq.2). Denote by
$F_{>0}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})$
and
$F_{<0}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})$
the versions of
$\mathbb{P}(W_{j}>0|\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})$
and
$\mathbb{P}(W_{j}<0|\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})$,
respectively. We will show that functions
$F_{>0}:\mathbb{R}^{n(2p+1)}\longmapsto\mathbb{R}$ and
$F_{<0}:\mathbb{R}^{n(2p+1)}\longmapsto\mathbb{R}$ satisfiy that

$\begin{split}F_{>0}(\vec{z})&=\frac{\boldsymbol{1}_{\\{w_{j}(\vec{z})>0\\}}\times
f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\vec{z})}{f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\vec{z})},\\\
F_{<0}(\vec{z})&=\frac{\boldsymbol{1}_{\\{w_{j}(\vec{z}_{\textnormal{swap}})<0\\}}\times
f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\vec{z}_{\textnormal{swap}})}{f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\vec{z})},\end{split}$
(A.13)

respectively, where
$\vec{z}=(\vec{z}_{1},\vec{z}_{2},\vec{z}_{3},\vec{z}_{4},\vec{z}_{5})$,
$\vec{z}_{\textnormal{swap}}=(\vec{z}_{2},\vec{z}_{1},\vec{z}_{3},\vec{z}_{4},\vec{z}_{5})$
with $\vec{z}_{1}\in\mathbb{R}^{n}$, $\vec{z}_{2}\in\mathbb{R}^{n}$,
$\vec{z}_{3}\in\mathbb{R}^{n(p-1)}$, $\vec{z}_{4}\in\mathbb{R}^{n(p-1)}$,
$\vec{z}_{5}\in\mathbb{R}^{n}$, and
$w_{j}:\mathbb{R}^{n(2p+1)}\longmapsto\mathbb{R}$ denotes the knockoff
statistic function of $W_{j}$. From the definitions of
$\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)}$, and $w_{j}(\cdot)$ along
with the sign-flip property (2), we have that almost surely,

$\begin{split}w_{j}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})&\geq
0,\\\
w_{j}(\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})&<0.\end{split}$
(A.14)

Then an application of (A.13)–(A.14) and the fact that the probability density
function is nonnegative yields that

$\begin{split}&F_{>0}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\leq\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}{f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})},\\\
&F_{<0}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})=\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}{f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})},\end{split}$

which entail that

$\begin{split}&\mathbb{P}(W_{j}>0|\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=F_{>0}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&\leq\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}{f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})},\\\
&\mathbb{P}(W_{j}<0|\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=F_{<0}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=\frac{f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}{f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})}.\end{split}$
(A.15)

Hence, a combination of (LABEL:robust.eq.4) and Condition 1 (which ensures
that the denominator is nonzero) establishes the result in
(LABEL:robust.eq.2).

It remains to prove (A.13). To this end, observe that for any Borel sets
$A_{1}\in\mathcal{R}^{n}$, $A_{2}\in\mathcal{R}^{n}$,
$A_{3}\in\mathcal{R}^{n(p-1)}$, $A_{4}\in\mathcal{R}^{n(p-1)}$, and
$A_{5}\in\mathcal{R}^{n}$, it holds that

$\begin{split}&\int_{\vec{z}\in A_{1}\times\cdots\times
A_{5}}F_{>0}(\vec{z})f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\vec{z})d\vec{z}\\\
&=\mathbb{P}(W_{j}>0,\boldsymbol{X}_{j}\in
A_{1},\widetilde{\boldsymbol{X}}_{j}\in A_{2},\boldsymbol{X}_{-j}\in
A_{3},\widetilde{\boldsymbol{X}}_{-j}\in A_{4},\boldsymbol{Y}\in A_{5})\\\
&=\mathbb{P}(W_{j}>0,\boldsymbol{X}_{j}^{(0)}\in
A_{1},\boldsymbol{X}_{j}^{(1)}\in A_{2},\boldsymbol{X}_{-j}\in
A_{3},\widetilde{\boldsymbol{X}}_{-j}\in A_{4},\boldsymbol{Y}\in A_{5})\\\
&=\int_{\\{\boldsymbol{X}_{j}^{(0)}\in A_{1},\boldsymbol{X}_{j}^{(1)}\in
A_{2},\boldsymbol{X}_{-j}\in A_{3},\widetilde{\boldsymbol{X}}_{-j}\in
A_{4},\boldsymbol{Y}\in A_{5}\\}}\mathbb{P}(W_{j}>0\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})d\mathbb{P},\end{split}$
(A.16)

where the second equality above holds by the definitions of
$\boldsymbol{X}_{j}^{(0)}$ and $\boldsymbol{X}_{j}^{(1)}$. Similarly, for any
Borel sets $A_{1}\in\mathcal{R}^{n}$, $A_{2}\in\mathcal{R}^{n}$,
$A_{3}\in\mathcal{R}^{n(p-1)}$, $A_{4}\in\mathcal{R}^{n(p-1)}$, and
$A_{5}\in\mathcal{R}^{n}$, we have

$\begin{split}&\int_{\vec{z}\in A_{1}\times\cdots\times
A_{5}}F_{<0}(\vec{z})f_{\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\vec{z})d\vec{z}\\\
&=\mathbb{P}(W_{j}<0,\widetilde{\boldsymbol{X}}_{j}\in
A_{1},\boldsymbol{X}_{j}\in A_{2},\boldsymbol{X}_{-j}\in
A_{3},\widetilde{\boldsymbol{X}}_{-j}\in A_{4},\boldsymbol{Y}\in A_{5})\\\
&=\mathbb{P}(W_{j}<0,\boldsymbol{X}_{j}^{(0)}\in
A_{1},\boldsymbol{X}_{j}^{(1)}\in A_{2},\boldsymbol{X}_{-j}\in
A_{3},\widetilde{\boldsymbol{X}}_{-j}\in A_{4},\boldsymbol{Y}\in A_{5})\\\
&=\int_{\\{\boldsymbol{X}_{j}^{(0)}\in A_{1},\boldsymbol{X}_{j}^{(1)}\in
A_{2},\boldsymbol{X}_{-j}\in A_{3},\widetilde{\boldsymbol{X}}_{-j}\in
A_{4},\boldsymbol{Y}\in A_{5}\\}}\mathbb{P}(W_{j}<0\big{|}\
\boldsymbol{X}_{j}^{(0)},\boldsymbol{X}_{j}^{(1)},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})d\mathbb{P},\end{split}$
(A.17)

where the second equality above holds by the definitions of
$\boldsymbol{X}_{j}^{(0)}$ and $\boldsymbol{X}_{j}^{(1)}$. With the aid of
(LABEL:robust.eq.12)–(LABEL:robust.eq.25), we can resort to the $\pi-\lambda$
Theorem [16] and the definition of the conditional expectation to obtain
(A.13). Thus, we have established (LABEL:robust.eq.2).

Proof of (13). We now aim to show (13) under Condition 2 and the assumption
that $\boldsymbol{X}_{j}$ is independent of $\boldsymbol{Y}$ conditional on
$\boldsymbol{X}_{-j}$ for each $j\in\mathcal{H}_{0}$. First, in light of
Condition 2 and Definition 2 of the knockoff generator, we see that
$\boldsymbol{Y}$ is independent of $\widetilde{\boldsymbol{X}}$ conditional on
$\boldsymbol{X}$. Next, we can deduce that

$\begin{split}&f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\vec{z}_{1},\vec{z}_{2},\vec{z}_{3},\vec{z}_{4},\vec{z}_{5})\\\
&=f_{\boldsymbol{X}}(\vec{z}_{1},\vec{z}_{3})\times
f_{\widetilde{\boldsymbol{X}}_{j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}|\boldsymbol{X}}(\vec{z}_{2},\vec{z}_{4},\vec{z}_{5}|\vec{z}_{1},\vec{z}_{3})\\\
&=f_{\boldsymbol{X}}(\vec{z}_{1},\vec{z}_{3})\times
f_{\boldsymbol{Y}|\boldsymbol{X}}(\vec{z}_{5}|\vec{z}_{1},\vec{z}_{3})\times
f_{\widetilde{\boldsymbol{X}}|\boldsymbol{X}}(\vec{z}_{2},\vec{z}_{4}|\vec{z}_{1},\vec{z}_{3})\\\
&=f_{\boldsymbol{X},\widetilde{\boldsymbol{X}}}(\vec{z}_{1},\vec{z}_{3},\vec{z}_{2},\vec{z}_{4})\times
f_{\boldsymbol{Y}|\boldsymbol{X}}(\vec{z}_{5}|\vec{z}_{1},\vec{z}_{3})\\\
&=f_{\boldsymbol{X},\widetilde{\boldsymbol{X}}}(\vec{z}_{1},\vec{z}_{3},\vec{z}_{2},\vec{z}_{4})\times
f_{\boldsymbol{Y}|\boldsymbol{X}_{-j}}(\vec{z}_{5}|\vec{z}_{3}),\end{split}$
(A.18)

where the second equality above follows from the conditional independence
property of the knockoffs and the last equality above holds because of the
columnwise conditional independence assumption on the null features.

From (LABEL:robust.eq.1), it holds that

$\begin{split}&f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j}}(\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j})\times
f_{\boldsymbol{Y}|\boldsymbol{X}_{-j}}(\boldsymbol{Y}|\boldsymbol{X}_{-j})\end{split}$

and

$\begin{split}&f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y}}(\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j},\boldsymbol{Y})\\\
&=f_{\boldsymbol{X}_{j},\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j}}(\widetilde{\boldsymbol{X}}_{j},\boldsymbol{X}_{j},\boldsymbol{X}_{-j},\widetilde{\boldsymbol{X}}_{-j})\times
f_{\boldsymbol{Y}|\boldsymbol{X}_{-j}}(\boldsymbol{Y}|\boldsymbol{X}_{-j}),\end{split}$

which establish (13).

Proofs of (14) and (15). The proof of (14) is straightforward using the
additional assumption of i.i.d. observations and hence, is omitted here for
simplicity. We now focus on proving (15). Fixing a feature index $j$, let us
consider a random vector
$(\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}})$
such that $\widetilde{X}_{j}^{\dagger}$ is generated by the $j$th
coordinatewise knockoff generator $\kappa_{j}(\boldsymbol{x}_{-j},)$ given
$\boldsymbol{x}_{-j}$ and that
$\widetilde{\boldsymbol{z}}=(\widetilde{Z}_{j},\widetilde{\boldsymbol{z}}_{-j})$
is a knockoff vector of $(\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j})$
generated from the knockoff filter
$\kappa((\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}),)$, where
$\kappa_{j}$ and $\kappa$ are as given in Condition 3. In view of Condition 3,
we see that
$(\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j})$
and
$(\widetilde{Z}_{j},\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j})$
have the same distribution and the corresponding density functions exist,
which entail that for each
$(z_{1},z_{2},\vec{z}_{3},\vec{z}_{4})\in\mathbb{R}^{2p}$,

$f_{\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j}}(z_{1},z_{2},\vec{z}_{3},\vec{z}_{4})=f_{\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j}}(z_{2},z_{1},\vec{z}_{3},\vec{z}_{4}).$
(A.19)

Next, it follows from the definition of
$(\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j})$
that

$\begin{split}f_{\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j}}(z_{1},z_{2},\vec{z}_{3},\vec{z}_{4})&=f_{\boldsymbol{x}_{-j}}(\vec{z}_{3})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(z_{1}|\vec{z}_{3})f_{\widetilde{Z}_{j},\widetilde{\boldsymbol{z}}_{-j}|\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(z_{2},\vec{z}_{4}|z_{1},\vec{z}_{3})\\\
&=f_{\boldsymbol{x}_{-j}}(\vec{z}_{3})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(z_{1}|\vec{z}_{3})f_{\widetilde{X}_{j},\widetilde{\boldsymbol{x}}_{-j}|X_{j},\boldsymbol{x}_{-j}}(z_{2},\vec{z}_{4}|z_{1},\vec{z}_{3}),\end{split}$
(A.20)

where
$f_{\widetilde{Z}_{j},\widetilde{\boldsymbol{z}}_{-j}|\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(z_{2},\vec{z}_{4}|z_{1},\vec{z}_{3})=f_{\widetilde{X}_{j},\widetilde{\boldsymbol{x}}_{-j}|X_{j},\boldsymbol{x}_{-j}}(z_{2},\vec{z}_{4}|z_{1},\vec{z}_{3})$
because a knockoff generator outputs random vectors with the same distribution
if the input values are the same due to Definition 2. Similarly, it holds that

$\begin{split}f_{\widetilde{X}_{j}^{\dagger},\widetilde{Z}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{z}}_{-j}}(z_{2},z_{1},\vec{z}_{3},\vec{z}_{4})&=f_{\boldsymbol{x}_{-j}}(\vec{z}_{3})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(z_{2}|\vec{z}_{3})f_{\widetilde{Z}_{j},\widetilde{\boldsymbol{z}}_{-j}|\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(z_{1},\vec{z}_{4}|z_{2},\vec{z}_{3})\\\
&=f_{\boldsymbol{x}_{-j}}(\vec{z}_{3})f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(z_{2}|\vec{z}_{3})f_{\widetilde{X}_{j},\widetilde{\boldsymbol{x}}_{-j}|X_{j},\boldsymbol{x}_{-j}}(z_{1},\vec{z}_{4}|z_{2},\vec{z}_{3}).\end{split}$
(A.21)

From (A.19)–(A.21), we can deduce that

$\begin{split}f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(z_{1}|\vec{z}_{3})f_{\widetilde{X}_{j},\widetilde{\boldsymbol{x}}_{-j}|X_{j},\boldsymbol{x}_{-j}}(z_{2},\vec{z}_{4}|z_{1},\vec{z}_{3})=f_{\widetilde{X}_{j}^{\dagger}|\boldsymbol{x}_{-j}}(z_{2}|\vec{z}_{3})f_{\widetilde{X}_{j},\widetilde{\boldsymbol{x}}_{-j}|X_{j},\boldsymbol{x}_{-j}}(z_{1},\vec{z}_{4}|z_{2},\vec{z}_{3}),\end{split}$

which results in

$\begin{split}\frac{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(z_{1},z_{2},\vec{z}_{3},\vec{z}_{4})}{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(z_{2},z_{1},\vec{z}_{3},\vec{z}_{4})}=\frac{f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(z_{2},\vec{z}_{3})\times
f_{X_{j},\boldsymbol{x}_{-j}}(z_{1},\vec{z}_{3})}{f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(z_{1},\vec{z}_{3})\times
f_{X_{j},\boldsymbol{x}_{-j}}(z_{2},\vec{z}_{3})}.\end{split}$ (A.22)

Setting
$(z_{1},z_{2},\vec{z}_{3},\vec{z}_{4})=(X_{ij},\widetilde{X}_{ij},\boldsymbol{x}_{-ij},\widetilde{\boldsymbol{x}}_{-ij})$
in (A.22) above, we can obtain that

$\begin{split}&\sum_{i=1}^{n}\log{\left(\frac{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(X_{ij},\widetilde{X}_{ij},\boldsymbol{x}_{-ij},\widetilde{\boldsymbol{x}}_{-ij})}{f_{X_{j},\widetilde{X}_{j},\boldsymbol{x}_{-j},\widetilde{\boldsymbol{x}}_{-j}}(\widetilde{X}_{ij},X_{ij},\boldsymbol{x}_{-ij},\widetilde{\boldsymbol{x}}_{-ij})}\right)}\\\
&=\sum_{i=1}^{n}\log{\left(\frac{f_{X_{j},\boldsymbol{x}_{-j}}(X_{ij},\boldsymbol{x}_{-ij})f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(\widetilde{X}_{ij},\boldsymbol{x}_{-ij})}{f_{X_{j},\boldsymbol{x}_{-j}}(\widetilde{X}_{ij},\boldsymbol{x}_{-ij})f_{\widetilde{X}_{j}^{\dagger},\boldsymbol{x}_{-j}}(X_{ij},\boldsymbol{x}_{-ij})}\right)},\end{split}$

which establishes (15). This concludes the proof of Theorem 2.

###  A.3 Proof of Theorem 3

Proof of (17). We begin with proving the first assertion (17) in Theorem 3.
The notation of $W_{j}$ and $T$ is the same as in Algorithm 2. By the
definition of $T$ and the assumption that $T<\infty$, we can show that

$1+\\#\\{j:W_{j}\leq-T\\}\leq\tau_{1}\times\\#\\{j:W_{j}\geq T\\}.$

Moreover, we have that

$e_{j}=\frac{p}{\\#\\{s:W_{s}\leq-T\\}+1}>0$

if and only if $j\in\\{s:W_{s}\geq T\\}$. Thus, for nonzero $e_{j}$’s or
equivalently any $j\in\\{s:W_{s}\geq T\\}$, it holds that

$\begin{split}e_{j}&\geq p\times(\tau_{1}\times\\#\\{j:W_{j}\geq T\\})^{-1}\\\
&\geq p\times(\tau^{*}\times\\#\\{j:W_{j}\geq T\\})^{-1},\end{split}$

where the second inequality above is from the assumption that
$\tau_{1}\leq\tau^{*}$. Using the above result and the construction of the
e-BH procedure, we can show that

$\\{j:W_{j}\geq T\\}\subset\widehat{S}.$ (A.23)

In light of (A.23), we need only to prove that conditional on event
$\\{\sum_{j=1}^{2p}|\widehat{\beta}_{j}-\beta_{j}^{*}|\leq
c_{0}(\\#S^{*})\lambda_{n}\\}\cap\\{\\#\\{j:W_{j}\geq T\\}\geq
c_{1}(\\#S^{*})\\}$, it holds that

$\begin{split}\frac{\\#(S^{*}\cap\\{j:W_{j}\geq T\\})}{\\#S^{*}}\geq
1-(1+\phi)c_{0}k_{1n}^{-1}.\end{split}$ (A.24)

Then the first assertion of Theorem 3 follows from Conditions 5–6.

We now proceed with establishing (A.24). Assume without loss of generality
that

$|W_{1}|\geq\cdots\geq|W_{p}|.$

Let $j^{*}\in\\{1,\cdots,p\\}$ be given such that $j^{*}\in\\{s:|W_{s}|=T\\}$.
Such $j^{*}$ always exists because of the assumption that $T<\infty$. Then it
follows that

$-T<W_{j^{*}+1}\leq 0$

by the definition of $T$ (because otherwise $T$ would take a smaller value
than $|W_{j*}|$) and the assumption that there are no ties in
$\\{|W_{j}|:|W_{j}|>0\\}$. We will analyze two cases separately, where the
first case considers $W_{j^{*}+1}=0$ and the second case considers
$-T<W_{j^{*}+1}<0$.

Let us consider the first case of $W_{j^{*}+1}=0$. Denote by
$\widetilde{q}=\phi c_{0}k_{1n}^{-1}$ with $\phi>0$ and $\phi^{2}-\phi-1=0$.
We will discuss the scenarios of
$\\#\\{j:W_{j}<0\\}\leq\widetilde{q}(\\#S^{*})$ and
$\\#\\{j:W_{j}<0\\}>\widetilde{q}(\\#S^{*})$ separately, where the former case
will be examined here and the latter one will be left to a later part. For the
scenario of $\\#\\{j:W_{j}<0\\}\leq\widetilde{q}(\\#S^{*})$, some simple
calculations together with $W_{j*+1}=0$ give that

$\begin{split}\\#(\\{j:W_{j}\geq T\\}\cap S^{*})&=\\#(\\{j:|W_{j}|>0\\}\cap
S^{*})-\\#(\\{j:W_{j}<0\\}\cap S^{*})\\\ &\geq\\#(\\{j:|W_{j}|>0\\}\cap
S^{*})-\widetilde{q}(\\#S^{*}).\end{split}$ (A.25)

We will deal with the term $\\#(\\{j:|W_{j}|>0\\}\cap S^{*})$ on the RHS of
(A.25) below. On the event
$\\{\sum_{j=1}^{2p}|\widehat{\beta}_{j}-\beta_{j}^{*}|\leq
c_{0}(\\#S^{*})\lambda_{n}\\}$, we can deduce that

$\begin{split}c_{0}\lambda_{n}(\\#S^{*})&\geq\sum_{j\in\widehat{S}_{1}\cap
S^{*}}|\beta_{j}^{*}|\\\ &\geq\\#(\widehat{S}_{1}\cap S^{*})\times(\min_{j\in
S^{*}}|\beta_{j}^{*}|),\end{split}$

where $\widehat{S}_{1}=\\{j:|\widehat{\beta}_{j}|=0\\}$. Such result and
Condition 6 entail that

$c_{0}(\\#S^{*})k_{1n}^{-1}\geq\\#(\widehat{S}_{1}\cap S^{*}).$

Hence, it follows from the assumption that there are no ties in
$\\{|\widehat{\beta}_{j}|:|\widehat{\beta}_{j}|>0\\}$ that

$\begin{split}\\#(\\{j:|W_{j}|>0\\}\cap S^{*})&=\\#((\widehat{S}_{1})^{c}\cap
S^{*})\\\ &\geq(1-c_{0}k_{1n}^{-1})\times(\\#S^{*}).\end{split}$ (A.26)

Then combining (A.25)–(A.26), we can obtain that conditional on event
$\\{\sum_{j=1}^{2p}|\widehat{\beta}_{j}-\beta_{j}^{*}|\leq
c_{0}(\\#S^{*})\lambda_{n}\\}$,

$\begin{split}\frac{\\#(\\{j:W_{j}\geq T\\}\cap S^{*})}{\\#S^{*}}&\geq
1-c_{0}k_{1n}^{-1}-\widetilde{q}\\\ &=1-(1+\phi)c_{0}k_{1n}^{-1},\end{split}$
(A.27)

which establishes (A.24). Moreover, observe that the second scenario of
$\\#\\{j:W_{j}<0\\}>\widetilde{q}(\\#S^{*})$ when $W_{j^{*}+1}=0$ implies that

$\\#\\{j:W_{j}\leq-T\\}=\\#\\{j:W_{j}<0\\}>\widetilde{q}(\\#S^{*})>\widetilde{q}(\\#S^{*}),$

which reduces to the same form as in (A.28) below. Thus, the proof provided
below can be applied here to conclude the proof for the first case
$W_{j^{*}+1}=0$.

We now consider the case of $-T<W_{j^{*}+1}<0$. From the definition of $T$ and
the assumption that there are no ties in $\\{|W_{j}|:|W_{j}|>0\\}$, it holds
on event $\\{\\#\\{j:W_{j}\geq T\\}\geq c_{1}(\\#S^{*})\\}$ that

$\begin{split}\\#\\{j:W_{j}\leq-T\\}+2&\geq\tau_{1}\times\\#\\{j:W_{j}\geq
T\\}\\\ &\geq\tau_{1}c_{1}(\\#S^{*}).\end{split}$ (A.28)

Meanwhile, on the event
$\\{\sum_{j=1}^{2p}|\widehat{\beta}_{j}-\beta_{j}^{*}|\leq
c_{0}(\\#S^{*})\lambda_{n}\\}$, since $\beta_{j+p}^{*}=0$ for all $j>0$ and
$|\widehat{\beta}_{j+p}|\geq T+|\widehat{\beta}_{j}|$ for all
$j\in\\{s:W_{s}\leq-T\\}$, we have that

$\begin{split}c_{0}\lambda_{n}(\\#S^{*})&\geq\sum_{j:W_{j}\leq-T}|\widehat{\beta}_{j+p}-\beta_{j+p}^{*}|\\\
&=\sum_{j:W_{j}\leq-T}|\widehat{\beta}_{j+p}|\\\
&\geq\\#\\{j:W_{j}\leq-T\\}\times T.\end{split}$ (A.29)

Then from (A.28)–(A.29) and Conditions 6–7, we can obtain that conditional on
event $\\{\\#\\{j:W_{j}\geq T\\}\geq
c_{1}(\\#S^{*})\\}\cap\\{\\#\\{j:W_{j}\geq T\\}\geq c_{1}(\\#S^{*})\\}$,

$\begin{split}T\leq\frac{c_{0}\lambda_{n}(\\#S^{*})}{c_{1}\tau_{1}(\\#S^{*})-2}\leq
k_{1n}\lambda_{n}\phi^{-1}\end{split}$ (A.30)

for all large $n$.

Further, conditional on event
$\\{\sum_{j=1}^{2p}|\widehat{\beta}_{j}-\beta_{j}^{*}|\leq
c_{0}(\\#S^{*})\lambda_{n}\\}\cap\\{\\#\\{j:W_{j}\geq T\\}\geq
c_{1}(\\#S^{*})\\}$, we can deduce that for all large $n$,

$\begin{split}c_{0}\lambda_{n}(\\#S^{*})&\geq\sum_{j\in
S^{*}\cap(\\{j:W_{j}\geq
T\\})^{c}}(|\widehat{\beta}_{j}-\beta_{j}^{*}|+|\widehat{\beta}_{j+p}|)\\\
&\geq\sum_{j\in S^{*}\cap(\\{j:W_{j}\geq
T\\})^{c}}(|\widehat{\beta}_{j}-\beta_{j}^{*}|+|\widehat{\beta}_{j}|-T)\\\
&\geq\sum_{j\in S^{*}\cap(\\{j:W_{j}\geq T\\})^{c}}(|\beta_{j}^{*}|-T)\\\
&\geq\\#(S^{*}\cap(\\{j:W_{j}\geq T\\})^{c})\times
k_{1n}(1-\phi^{-1})\lambda_{n},\end{split}$ (A.31)

where the second inequality above is from the fact that
$|\widehat{\beta}_{j+p}|>|\widehat{\beta}_{j}|-T$ for $j$ in $\\{j:W_{j}\geq
T\\}^{c}$, the third inequality above is due to the triangle inequality, and
the last inequality above results from Condition 6 and (A.30).

In light of (A.31), it holds on event
$\\{\sum_{j=1}^{2p}|\widehat{\beta}_{j}-\beta_{j}^{*}|\leq
c_{0}(\\#S^{*})\lambda_{n}\\}\cap\\{\\#\\{j:W_{j}\geq T\\}\geq
c_{1}(\\#S^{*})\\}$ that for all large $n$,

$\begin{split}\frac{\\#(S^{*}\cap\\{j:W_{j}\geq
T\\})}{\\#S^{*}}&=1-\frac{\\#(S^{*}\cap(\\{j:W_{j}\geq
T\\})^{c})}{\\#S^{*}}\\\ &\geq 1-\frac{c_{0}}{k_{1n}(1-\phi^{-1})}\\\
&=1-(1+\phi)c_{0}k_{1n}^{-1},\end{split}$ (A.32)

where the second equality above follows from the definition of $\phi$. This
establishes (A.24). Thus, combining the above results concludes the proof for
the first assertion of Theorem 3.

Proof of (19). We now aim to prove the third assertion (19) of Theorem 3, and
will defer the proof of (18) to the end. Let the knockoff thresholds
$T^{k}$’s, knockoff statistics $W_{j}^{k}$’s, statistics $e_{j}^{k}$’s, and
e-values $e_{j}$’s be given as in Algorithm 1. Let us first outline the proof
idea for (19) and we will build our proof on the first assertion of Theorem 3.
Using the inclusion-exclusion principle, we will show that
$\cap_{k=1}^{q+1}\\{j:e_{j}^{k}>0\\}$ includes most features in $S^{*}$ with
high probability. We then prove that each $e_{j}$ with $j$ in
$\cap_{k=1}^{q+1}\\{j:e_{j}^{k}>0\\}$ is sufficiently large to be selected by
the e-BH procedure, that is,
$\cap_{k=1}^{q+1}\\{j:e_{j}^{k}>0\\}\subset\widehat{S}$. Combining all these
results will complete the proof of (19). We will provide the full details of
the proof next.

First, recall that $\widehat{S}=\\{j:e_{j}\geq
p(\tau^{*}\times\widehat{k})^{-1}\\}$ with $\widehat{k}=\max\\{k:e_{(k)}\geq
p(\tau^{*}\times k)^{{-1}}\\}$, where $e_{(j)}$’s are the ordered statistics
of $e_{j}$’s such that $e_{(1)}\geq\cdots\geq e_{(p)}$. Let $K>1$ be the
constant specified in Theorem 3. Let us consider two events given by

$\cap_{k=1}^{{q+1}}\\{\\#\\{j:W_{j}^{k}\geq T^{k}\\}\leq K(\\#S^{*})\\}$
(A.33)

and

$\cap_{k=1}^{q+1}\left\\{\frac{\\#(S^{*}\cap\\{j:W_{j}^{k}\geq
T^{k}\\})}{\\#S^{*}}\geq 1-(1+\phi)c_{0}k_{1n}^{-1}\right\\}.$ (A.34)

We will show that conditional on the two events in (A.33) and (A.34) above, it
holds that

$e_{(\\#\mathcal{S})}=\min_{j\in\mathcal{S}}e_{j}\geq
p(\tau^{*}(\\#\mathcal{S}))^{-1},$ (A.35)

where

$\mathcal{S}\coloneqq\cap_{k=1}^{q+1}\\{j:W_{j}^{k}\geq T^{k}\\}.$ (A.36)

Then it follows from (A.35) and the definition of $\widehat{k}$ that
$\widehat{k}\geq\\#\mathcal{S}$ and $\mathcal{S}\subset\widehat{S}$. Such
results along with an application of the inclusion–exclusion principle entail