# Boosted p-Values for High-Dimensional Vector Autoregression

Xiao Huang Department of Economics, Finance, and Quantitative Analysis, Coles
College of Business, Kennesaw State University, GA 30144, USA. Email:
<EMAIL_ADDRESS>
###### Abstract

Assessing the statistical significance of parameter estimates is an important
step in high-dimensional vector autoregression modeling. Using the least-
squares boosting method, we compute the p-value for each selected parameter at
every boosting step in a linear model. The p-values are asymptotically valid
and also adapt to the iterative nature of the boosting procedure. Our
simulation experiment shows that the p-values can keep false positive rate
under control in high-dimensional vector autoregressions. In an application
with more than $100$ macroeconomic time series, we further show that the
p-values can not only select a sparser model with good prediction performance
but also help control model stability. A companion R package boostvar is
developed.

JEL Classification: C18, C32

Keywords: Least-squares boosting, linear regression, vector autoregression,
p-values

##  1 Introduction

Since its introduction as a macroeconometric framework in Sims (1980), vector
autoregression (VAR) has become one of the most valuable econometric tools for
describing the joint dynamics of multiple time series. Its application can
also be found in many other disciplines such as biology and neuronscience.
With increasingly available data, there have been many developments in
estimating high-dimensional VARs. Ideas of using the lasso in Tibshirani
(1996) or imposing other penalties are flourishing. The use of the lasso or
its variants can be found in Fujita et al. (2007); Lozano et al. (2009); Basu
and Michailidis (2015); Medeiros and Mendes (2016); Nicholson et al. (2020);
Wong et al. (2020), among others. The Dantzig selector in Han et al. (2015)
and the tensor method in Wang et al. (2021) also directly address the
estimation problem in high-dimensional VARs. Methods such as those in Chen and
Huang (2012); Uematsu et al. (2019) are designed for multivariate linear
regression, and they can also be used for estimation. A data analyst now has
many methods for estimation in high-dimensional VARs.

One unanswered question in this active area of research is how to assign
p-values to parameter estimates in a high-dimensional VAR. In a standard
linear regression, we provide both parameter estimates and their p-values, and
the same should be expected in high-dimensional VARs. Without a p-value, it is
not possible to assess the statistical uncertainty associated with an
estimate. If the research purpose is model interpretation, obtaining a valid
p-value is perhaps as important as having a parameter estimate. If the goal is
prediction, the model is typically selected via cross-validation (c.v.).
However, cross-validation itself guarantees no statistical significance, and
knowing the p-values can help improve the prediction for certain data. Given
the abundant evidence in the literature on cross-section regression that tools
like the lasso can have high false positive rate, it is reasonable to
conjecture that many of the non-zero parameter estimates produced by the lasso
(and some other methods such as boosting) can be statistically insignificant.
A p-value can certainly help reveal better the underlying structure in the
data.

Since each estimation method generates its own sparsity, it is clear that we
need different ways of computing the standard error (s.e.) to reflect the
unique uncertainty embedded in each method. We choose to work with the least-
squares boosting (LS-Boost) method in this paper. LS-Boost is introduced in
Friedman (2001) as an example of gradient boosting algorithm. It can easily
handle a regression model with high dimension. Its application in high-
dimensional VAR estimation is discussed in Lutz and Bühlmann (2006). The focus
of this paper is not to develop another algorithm to generate a new sparsity
in high-dimensional VARs. Instead, we study the following question: If LS-
Boost is used in a high-dimensional VAR for estimation, can we derive a
standard error and compute the p-value for each selected parameter at every
boosting step? We give an affirmative answer to this question.

When the number of boosting step is large, the LS-Boost estimator converges to
the LS estimator. We provide a simulation example to show our p-value can
converge to the p-value of the LS estimator in a bivariate VAR. When the VAR
dimension is large, there is no unique LS solution to a VAR; the LS-Boost
estimator converges to one of the LS solutions as boosting step increases. For
LS-Boost, our method can compute the p-value for every estimate at a given
boosting step. Just like the lasso and many other penalized estimator, the LS-
Boost estimator is biased at a given boosting step. In this case, the t-test
uses a biased parameter estimator to test the null hypothesis of whether the
cumulative incremental changes in a parameter is statistically different from
$0$ at every boosting step (see more detailed discussions in Section 5). Let
$\beta_{j}$ be a scalar parameter for the variable $X_{j}$. The classical
fixed hypothesis is $H_{0}{:}\;\beta_{j}=0$. One is supposed to use an
unbiased estimator $\hat{\beta}_{j}$ for hypothesis testing. In this regard,
our null hypothesis is not equivalent to $H_{0}{:}\;\beta_{j}=0$. This,
however, should not be the ground on which to discount the usefulness of our
p-value. Consider a case of a cross-section regression with $1,000$ variables
and $100$ observations. To test $H_{0}{:}\;\beta_{j}=0$, we need a reliable
$\hat{\beta}_{j}$ and its s.e. But to obtain a $\hat{\beta}_{j}$ in a high-
dimensional regression, we will inevitably resort to some adaptive procedures
such as the lasso. As soon as one uses the lasso for even an initial estimate,
the lasso selection and estimation risk starts to contribute the overall
estimation risk of $\hat{\beta}_{j}$, and these additional risks must be
accounted for when testing the null, an intrinsically hard problem. In
addition, no matter which penalized method is used in high-dimensional model,
one will likely end up with a biased estimator for $\beta_{j}$ and use it in
fixed hypothesis testing. As a result, hypothesis testing in high dimensional
models is challenging. Various approaches have been proposed for fixed
hypotheses testing in high-dimensional cross-section regressions, including
Wasserman and Roeder (2009); Meinshausen et al. (2009); Zhang and Zhang
(2014); van de Geer et al. (2014), among others. These approaches typically
make some assumptions on the sparsity and the variables of the model. These
papers also have one thing in common: they all target the unknown population
parameter $\beta_{j}$ and try to test the hypothesis $H_{0}{:}\;\beta_{j}=0$.

We take a different approach in this paper. In high dimensional linear models,
since it is practically impossible to obtain an unbiased estimator for
$\beta_{j}$, can we estimate $\beta_{j}$ partially, in an incremental way, and
perform a t-test on these partial estimates? The new perspective this paper
explores is we can use a sequence of biased LS-Boost estimates to continuously
check the statistical significance of a coefficient up to each boosting step.
Put differently, if the numerator of a t-statistic is biased, as long as we
can obtain a valid s.e. for the biased estimator and correctly quantify the
statistical uncertainty of the numerator, the t-statistic will still have an
asymptotic normal distribution. This argument will be made more precise in the
discussion in Section 5.

A related line of research is the recent development in post-selection
inference for the lasso, where works in Lockhart et al. (2014); Lee et al.
(2016); Tibshirani et al. (2016) focus on how to assign a valid p-value to the
lasso estimator by adaptively incorporating the uncertainty in variable
selection. Our work differs from theirs in several aspects: (i) this paper
studies the issue of assigning p-values in a VAR for LS-Boost while their work
is focused on the lasso (and several other sequential procedures) in a cross-
section regression; (ii) more importantly, Lee et al. (2016); Tibshirani et
al. (2016) assume that data are normal and use polyhedral conditioning sets to
characterize the lasso selection event and to address the post-selection
inference in the lasso; whereas we make no assumption of data normality and
derive the estimator’s variance directly from its closed-form solution; (iii)
we also provide a tractable asymptotic distribution for each parameter
estimate at every boosting step. Our method of computing the p-values in VAR
nests cross-section regression as a special case. In this special case, we use
three common data sets to demonstrate the computation of p-values in a cross-
section regression. The discussion can be found in the online supplement.

This paper makes the following contributions to the literature on high-
dimensional VARs. First, we provide an asymptotic distribution result for
every selected parameter at each boosting step for a stationary VAR. We do not
make any sparsity assumption on the parameter; nor do we impose any assumption
on the relative growth rate between sample size and the dimension of the
model. The asymptotic variance allows us to construct a t-test easily. Second,
we extend the computational bound results for boosting in Freund et al. (2017)
to VAR. When combined with the asymptotic result, this allows us to
characterize the behavior of LS-Boost estimator when both the boosting step
and sample size go to infinity. Third, our simulation result shows that using
a p-value can help reduce the false positive rate (FPR) and improve the F
score of a model. Our application using the macroeconomic data in McCracken
and Ng (2016) further demonstrates that a p-value-adjusted model can give good
prediction performance and good model stability. With the help of p-value, we
can remove more than $95\%$ of the nonzero LS-Boost estimates while still
improving the prediction in this application. An R package boostvar that
implements our method can be found at https://github.com/xhuang20/boostvar.

The rest of the paper is organized as follows. Section 2 discusses two LS-
Boost algorithms in VAR. Section 3 derives the standard error and asymptotic
distribution for the LS-Boost estimator. Section 4 presents the computation
bounds for boosting in a VAR model. Simulation results are discussed in
Section 5, and Section 6 includes an application. Section 7 concludes. The
online supplement contains all proofs, additional discussions and tables.

##  2 Two LS-Boost algorithms in VAR

We begin by introducing the standard VAR model and some notations. Consider a
$d\times 1$ vector time series $y_{t}$ with $t=1,\cdots,T$. A $p$th-order VAR
(VAR($p$)) is defined as

$y_{t}=\phi_{1}y_{t-1}+\cdots+\phi_{p}y_{t-p}+u_{t},$ (1)

where $\phi_{1},\cdots,\phi_{p}$ are $d\times d$ matrices of coefficients and
$u_{t}$ is a $d\times 1$ vector of error terms. We omit the intercept in
formulating the model by assuming data are demeaned. The total number of
parameters is $p\times d^{2}$, which can easily exceed the number of
observations $T$ and make the LS method infeasible.

It is often convenient to rewrite eq. 1 in a matrix format. Let ′ denote
matrix/vector transpose. Define

$\underset{T\times d}{\mathbf{Y}}=\begin{bmatrix}y_{1}^{\prime}\\\ \vdots\\\
y_{T}^{\prime}\end{bmatrix},\underset{T\times
d}{\mathbf{Y}_{-1}}=\begin{bmatrix}y_{0}^{\prime}\\\ \vdots\\\
y_{T-1}^{\prime}\end{bmatrix},\cdots,\underset{T\times
d}{\mathbf{Y}_{-p}}=\begin{bmatrix}y_{1-p}^{\prime}\\\ \vdots\\\
y_{T-p}^{\prime}\end{bmatrix},\underset{T\times
d}{\mathbf{u}}=\begin{bmatrix}u_{1}^{\prime}\\\ \vdots\\\
u_{T}^{\prime}\end{bmatrix},\text{ and }\underset{pd\times
d}{\bm{\phi}}=\begin{bmatrix}\phi_{1}^{\prime}\\\ \vdots\\\
\phi_{p}^{\prime}\end{bmatrix}.$

In matrix form, eq. 1 can be written as

$\mathbf{Y}=\mathbf{Y}_{-1}\phi_{1}^{\prime}+\cdots+\mathbf{Y}_{-p}\phi_{p}^{\prime}+\mathbf{u}.$
(2)

The multivariate regression format is given by

$\mathbf{Y}=\mathbf{X}\bm{\phi}+\mathbf{u},$ (3)

where

$\mathbf{X}=\left[\mathbf{Y}_{-1},\cdots,\mathbf{Y}_{-p}\right]\text{ and
}\bm{\phi}=[\phi_{1},\cdots,\phi_{p}]^{\prime}.$ (4)

The basic idea of the LS-Boost is to recursively select a column of the
covariates in $\mathbf{X}$ that gives the best LS fit of $\mathbf{Y}$ in the
first boosting step or the current residuals in later boosting steps. There is
no unique way to implement the LS-Boost algorithm in a VAR. At each boosting
step, one can select a column in $\mathbf{X}$ that gives the best fit for a
selected column in $\mathbf{Y}$ or the residual matrix and update a single
element in $\bm{\phi}$, which is the componentwise linear least squares
procedure described in Lutz and Bühlmann (2006). Alternatively, one can use
the “row-boosting” procedure in Lutz and Bühlmann (2006) to select a column in
$\mathbf{X}$ to gives the best fit for the entire $\mathbf{Y}$ or the residual
matrix and update a row of entries in $\bm{\phi}$. In the following, we
describe a slightly more general procedure that nests the “row-boosting” as a
special case.

###  2.1 A Group (by variable) LS-Boost Algorithm

To motivate the group LS-boost procedure, let us consider a VAR(2) model with
three variables: unemployment rate (UNR), consumer price index (CPI), and
inflation rate (INF). Equation 3 becomes

$[\underbrace{\text{UNR},\text{CPI},\text{INF}}_{T\times
3}]=[\underbrace{\text{UNR}_{-1},\text{CPI}_{-1},\text{INF}_{-1},\text{UNR}_{-2},\text{CPI}_{-2},\text{INF}_{-2}}_{T\times
6}]\bm{\phi}+\mathbf{u}.$ (5)

Our first algorithm works by selecting a variable with all its lags that gives
the best fit at each boosting step. Hence, if the variable CPI is found to
give the best fit for $\mathbf{Y}$, the matrix
$[\text{CPI}_{-1},\text{CPI}_{-2}]$ will be the selected variable matrix for
that step. Since we select each variable with all its lags as a group, it will
be helpful to rearrange columns in eq. 4 so that they are grouped by variables

$\mathbf{X}_{g}=[\mathbf{X}_{(1)},\cdots,\mathbf{X}_{(j)},\cdots,\mathbf{X}_{(d)}]_{T\times
pd},$ (6)

where $X_{(j)}$ is a $T\times p$ matrix that collects $p$ columns of lags in
variable $j$ such as $[\text{CPI}_{-1},\text{CPI}_{-2}]$. And the VAR($p$)
process in eq. 1 can be written as

$\displaystyle\mathbf{Y}$
$\displaystyle=\mathbf{X}_{(1)}\phi_{(1)}^{\prime}+\cdots+\mathbf{X}_{(j)}\phi_{(j)}^{\prime}+\cdots+\mathbf{X}_{(d)}\phi_{(d)}^{\prime}+\mathbf{u}$
$\displaystyle=\mathbf{X}_{g}\bm{\phi}_{g}+\mathbf{u},$ (7)

where
$\bm{\phi}_{g}=[\phi_{(1)},\cdots,\phi_{(j)},\cdots,\phi_{(d)}]^{\prime}$.
$\phi_{(j)}^{\prime}$ is a $p\times d$ coefficient matrix, each row of which
comes from the $j$th row of $\phi_{1}^{\prime},\cdots,\phi_{p}^{\prime}$.

Let $k$ denote the iteration step and $k=0,1,\cdots,k_{\text{stop}}$, and
$k_{\text{stop}}$ is a prespecified stopping number in boosting iteration. Let
the superscript $(k)$ denote a quantity associated with step $k$ so that, for
example, $\hat{\mathbf{R}}^{(k)}$ is the $T\times d$ residual matrix generated
at step $k$. Because we select all lags of a variable at each step, the LS
regression at step $k$ will take the following form for variable $j$

$\underset{T\times
d}{\vphantom{\phi_{(j)}^{\prime}}\hat{\mathbf{R}}^{(k-1)}}=\underset{T\times
p}{\vphantom{\phi_{(j)}^{\prime}}\mathbf{X}_{(j)}}\underset{p\times
d}{\phi_{(j)}^{\prime}}+\text{error term},$ (8)

where $\hat{\mathbf{R}}^{(k-1)}=\hat{\mathbf{R}}^{(0)}=\mathbf{Y}$ if $k=1$.

If we define the $p\times d$ matrix

$\mathbf{A}_{j}=(\mathbf{X}_{(j)}^{\prime}\mathbf{X}_{(j)})^{-1}\mathbf{X}_{(j)}^{\prime},$
(9)

the estimate in eq. 8 with a learning rate $\nu$, can be updated as

$\hat{\phi}_{(j)}^{(k)\prime}=\hat{\phi}_{(j)}^{(k-1)\prime}+\nu\mathbf{A}_{j}\hat{\mathbf{R}}^{(k-1)}.$
(10)

Initialize
$\hat{\phi}_{(1)}^{(0)}=\cdots=\hat{\phi}_{(d)}^{(0)}=\mathbf{0}_{d\times p}$.
Let $\lVert\cdot\rVert_{2}$ denote the matrix and vector Euclidean norm. For
each boosting iteration $k\geq 1$, the algorithm becomes

  1. Step 1. Select the variable $\mathbf{X}_{(j_{k})}$ with

$j_{k}\in\operatorname*{argmin}_{1\leq j\leq
d}\frac{1}{Td}\lVert\hat{\mathbf{R}}^{(k-1)}-\mathbf{X}_{(j)}\beta_{(j)}^{(k)\prime}\rVert_{2}^{2}\text{
with }\hat{\beta}_{(j)}^{(k)\prime}=\mathbf{A}_{j}\hat{\mathbf{R}}^{(k-1)}.$
(11)

  2. Step 2. Update $\hat{\phi}_{(j)}^{(k)}$ and $\hat{\mathbf{R}}^{(k)}$.

$\displaystyle\hat{\phi}_{(j)}^{(k)\prime}$
$\displaystyle=\begin{cases}\hat{\phi}_{(j_{k})}^{(k-1)\prime}+\nu\hat{\beta}_{(j_{k})}^{(k)\prime}&\text{if
}j=j_{k},\\\ \hat{\phi}_{(j)}^{(k-1)\prime}&\text{if }j\neq j_{k},\end{cases}$
(12) $\displaystyle\hat{\mathbf{R}}^{(k)}$
$\displaystyle=\hat{\mathbf{R}}^{(k-1)}-\nu\mathbf{X}_{(j_{k})}\hat{\beta}_{(j_{k})}^{(k)\prime}.$
(13)

Update several other quantities.

$\displaystyle\mathbf{A}_{j}^{(k)}$
$\displaystyle=\begin{cases}\mathbf{A}_{j_{k}}&\text{if }j=j_{k},\\\
0&\text{if }j\neq j_{k},\end{cases}$ (14) $\displaystyle\mathbf{H}^{(k)}$
$\displaystyle=\mathbf{X}_{(j_{k})}(\mathbf{X}_{(j_{k})}^{\prime}\mathbf{X}_{(j_{k})})^{-1}\mathbf{X}_{(j_{k})}^{\prime},$
(15) $\displaystyle\mathbf{B}^{(k)}$
$\displaystyle=\mathbf{I}_{T}-(\mathbf{I}_{T}-\nu\mathbf{H}^{(k)})\cdots(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)}),$
(16) $\displaystyle\hat{\mathbf{F}}^{(k)}$
$\displaystyle=\hat{\mathbf{F}}^{(1)}+\cdots+\hat{\mathbf{F}}^{(k-1)}+\nu\mathbf{H}^{(k)}\mathbf{R}^{(k-1)}=\mathbf{B}^{(k)}\mathbf{Y},$
(17) $\displaystyle\hat{\mathbf{u}}^{(k)}$
$\displaystyle=\mathbf{Y}-\hat{\mathbf{F}}^{(k)}=(\mathbf{I}_{T}-\mathbf{B}^{(k)})\mathbf{Y}.$
(18)

  3. Step 3. Repeat Steps 1 and 2 many times until step $k_{\text{stop}}$.

We call this algorithm LS-Boost1. The matrix $\mathbf{B}^{(k)}$ in eq. 16 is
the common hat matrix, and the sum of its diagonal elements gives the degree
of freedom up to boosting step $k$.

###  2.2 A sparser LS-Boost Algorithm

An approach with a sparser solution is to, instead of using all lags of
variable $\mathbf{X}_{(j)}$ in eq. 8 at each boosting step $k$, select only a
single column of $\mathbf{X}_{(j)}$ to fit $\hat{\mathbf{R}}^{k-1}$. There is
no guarantee that a recent lag will enter the model earlier than a distant
lag. For example, the column $\text{CPI}_{-2}$ can enter the model earlier
than $\text{CPI}_{-1}$ does. Let $\mathbf{X}_{(j)s}$ be the $s$th lag of the
$j$th variable (the $s$th column in $\mathbf{X}_{(j)}$) and
$\phi_{(j)s}^{\prime}$ be the $s$th row of $\phi_{(j)}^{\prime}$, the
regression equation for the “row-boosting” algorithm can be written as

$\underset{T\times
d}{\vphantom{\phi_{(j)}^{\prime}}\hat{\mathbf{R}}^{(k-1)}}=\underset{T\times
1}{\vphantom{\phi_{(j)}^{\prime}}\mathbf{X}_{(j)s}}\underset{1\times
d}{\phi_{(j)s}^{\prime}}+\text{error term},$ (19)

where $\phi_{(j)s}^{\prime}$ is a $1\times d$ row vector, which is the $s$th
column of $\phi_{(s)}$. Define

$\underset{1\times
T}{\vphantom{\mathbf{X}_{(j)s}^{\prime}}\mathbf{A}_{js}}=(\mathbf{X}_{(j)s}^{\prime}\mathbf{X}_{(j)s})^{-1}\mathbf{X}_{(j)s}^{\prime}.$
(20)

The update equation for the coefficient estimate at step $k$ is given by

$\hat{\phi}_{(j)s}^{(k)\prime}=\hat{\phi}_{(j)s}^{(k-1)\prime}+\nu\mathbf{A}_{js}\hat{\mathbf{R}}^{(k-1)}.$
(21)

Initialize $\hat{\phi}_{(j)s}^{(0)}=\mathbf{0}_{d\times 1}\forall
j=1,\cdots,d$ and $s=1,\cdots,p$. For each boosting step $k\geq 1$, the
algorithm becomes

  1. Step 1. Select the variable $\mathbf{X}_{(j_{k})s_{k}}$ with

$j_{k},s_{k}\in\operatorname*{argmin}_{\begin{subarray}{c}1\leq j\leq d\\\
1\leq s\leq
p\end{subarray}}\frac{1}{Td}\lVert\hat{\mathbf{R}}^{(k-1)}-\mathbf{X}_{(j)s}\beta_{(j)s}^{(k)\prime}\rVert_{2}^{2}\text{
with }\hat{\beta}_{(j)s}^{(k)\prime}=\mathbf{A}_{js}\hat{\mathbf{R}}^{(k-1)}.$
(22)

  2. Step 2. Update $\hat{\phi}_{(j)s}^{(k)}$ and $\mathbf{R}^{(k)}$.

$\displaystyle\hat{\phi}_{(j)s}^{(k)\prime}$
$\displaystyle=\begin{cases}\hat{\phi}_{(j_{k})s_{k}}^{(k-1)\prime}+\nu\hat{\beta}_{(j_{k})s_{k}}^{(k)\prime}&\text{if
}j=j_{k}\text{ and }s=s_{k},\\\
\hat{\phi}_{(j)s}^{(k-1)\prime}&\text{otherwise},\end{cases}$ (23)
$\displaystyle\hat{\mathbf{R}}^{(k)}$
$\displaystyle=\hat{\mathbf{R}}^{(k-1)}-\nu\mathbf{X}_{(j_{k})s_{k}}\hat{\beta}_{(j_{k})s_{k}}^{(k)\prime}.$
(24)

Update several other quantities.

$\displaystyle\mathbf{A}_{js}^{(k)}$
$\displaystyle=\begin{cases}\mathbf{A}_{j_{k}s_{k}}&\text{if }j=j_{k}\text{
and }s=s_{k}\\\ 0&\text{otherwise},\end{cases}$ (25)
$\displaystyle\mathbf{H}^{(k)}$
$\displaystyle=\mathbf{X}_{(j_{k})s_{k}}(\mathbf{X}_{(j_{k})s_{k}}^{\prime}\mathbf{X}_{(j_{k})s_{k}})^{-1}\mathbf{X}_{(j_{k})s_{k}}^{\prime},$
(26) $\displaystyle\mathbf{B}^{(k)}$
$\displaystyle=\mathbf{I}_{T}-(\mathbf{I}_{T}-\nu\mathbf{H}^{(k)})\cdots(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)}),$
(27) $\displaystyle\hat{\mathbf{F}}^{(k)}$
$\displaystyle=\hat{\mathbf{F}}^{(1)}+\cdots+\hat{\mathbf{F}}^{(k-1)}+\nu\mathbf{H}^{(k)}\mathbf{R}^{(k-1)}=\mathbf{B}^{(k)}\mathbf{Y},$
(28) $\displaystyle\hat{\mathbf{u}}^{(k)}$
$\displaystyle=\mathbf{Y}-\hat{\mathbf{F}}^{(k)}=(\mathbf{I}_{T}-\mathbf{B}^{(k)})\mathbf{Y}.$
(29)

  3. Step 3. Repeat Steps 1 and 2 many times until step $k_{\text{stop}}$.

We call this algorithm LS-Boost2. Note that, although eqs. 27, 28 and 29 may
look the same as eqs. 16, 17 and 18, their values are different since eq. 26
is based on $\mathbf{X}_{(j_{k})s_{k}}$.

A few remarks are in order.

######  Remark 1.

Since the dimension of $\phi_{(j)}$ in eq. 8 is $d\times p$, the total number
of nonzero parameters will increase by $dp$ if a new variable is added. When
the dimension of $d$ is large, $dp$ can be large. However, since we have $d$
equations, on average, each equation will add $p$ more parameters, and $p$ is
usually a small number such as $2$ or $3$. Thus the model size for each of the
$d$ equation will not increase significantly when using LS-Boost1.

######  Remark 2.

LS-Boost2 is more flexible than LS-Boost1 and will generate a sparser model in
general. LS-Boost1 preserves the time series structure throughout variable
selection and estimation, and it can work better for certain data.

######  Remark 3.

The most flexible boosting algorithm is to update the parameters one at a
time. While the previous two algorithms will add either $dp$ or $d$ parameters
for a selected variable, we can choose to update a single parameter at each
boosting iteration step, which will require us to take many sweeps in the
data. This approach, called the componentwise procedure in Lutz and Bühlmann
(2006), completely ignores the time series structure of the data. Its
implementation will be similar to that of LS-Boost in a cross-section
regression. We focus on the use of LS-Boost1 and LS-Boost2 in this paper.

##  3 The asymptotic standard error and the p-value

We discuss the asymptotic result for the two boosting methods introduced in
Section 2.

###  3.1 Asymptotic results for LS-Boost1 when $k$ is fixed

We first discuss the computation of the standard error for LS-Boost1 at each
boosting step. Using the definition of $\hat{\beta}_{(j)}^{(k)\prime}$ in eq.
11, rewrite eq. 10 in a recursive form to have

$\displaystyle\hat{\phi}_{(j)}^{(k)\prime}$
$\displaystyle=\hat{\phi}_{(j)}^{(k-1)\prime}+\nu\hat{\beta}_{(j)}^{(k)\prime}=\hat{\phi}_{(j)}^{(k-2)\prime}+\nu\hat{\beta}_{(j)}^{(k-1)\prime}+\nu\hat{\beta}_{(j)}^{(k)\prime}=\cdots$
$\displaystyle=\hat{\phi}_{(j)}^{(0)\prime}+(\nu\hat{\beta}_{(j)}^{(1)\prime}+\cdots+\nu\hat{\beta}_{(j)}^{(k)\prime})$
$\displaystyle=\nu\hat{\beta}_{(j)}^{(1)\prime}+\cdots+\nu\hat{\beta}_{(j)}^{(k)\prime}\qquad\text{
since }\hat{\phi}_{(j)}^{(0)\prime}=0.$ (30)

Following the definitions in eqs. 14, 15, 16, 17 and 18, section 3.1 can be
further written as

$\displaystyle\hat{\phi}_{(j)}^{(k)\prime}$
$\displaystyle=\nu\mathbf{A}_{j}^{(1)}\mathbf{Y}+\nu\mathbf{A}_{j}^{(2)}(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)})\mathbf{Y}+\nu\mathbf{A}_{j}^{(3)}(\mathbf{I}_{T}-\nu\mathbf{H}^{(2)})(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)})\mathbf{Y}$
$\displaystyle\quad+\cdots+\nu\mathbf{A}_{j}^{(k)}(\mathbf{I}_{T}-\nu\mathbf{H}^{(k-1)})\cdots(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)})\mathbf{Y}$
$\displaystyle=\tilde{\mathbf{A}}_{j}^{(k)}\mathbf{Y},$ (31)

where we define

$\tilde{\mathbf{A}}_{j}^{(k)}=\left[\nu\mathbf{A}_{j}^{(1)}+\nu\mathbf{A}_{j}^{(2)}(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)})+\cdots+\nu\mathbf{A}_{j}^{(k)}(\mathbf{I}_{T}-\nu\mathbf{H}^{(k-1)})\cdots(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)})\right].$
(32)

Section 3.1 is a matrix form and it holds for all $j=1,\cdots,d$ and any
iteration step $k$, regardless of whether the $j$th variable (and its lags) is
selected at boosting step $k$. Recall the definition of $\mathbf{A}_{j}^{(k)}$
in eq. 14. If the $j$th variable is never selected up to the iteration step
$k$, we have $\mathbf{A}_{j}^{(1)}=\cdots=\mathbf{A}_{j}^{(k)}=0$ and
$\hat{\phi}_{(j)}^{(k)\prime}=0$; if the $j$th variable is selected in step
$k-1$ but not $k$, we have $\mathbf{A}_{j}^{(k-1)}\neq 0$ and
$\mathbf{A}_{j}^{(k)}=0$.

For LS-Boost1, we can use $\mathbf{X}_{g}$ in eq. 6 to rewrite

$\displaystyle\hat{\phi}_{(j)}^{(k)\prime}$
$\displaystyle=\tilde{\mathbf{A}}_{j}^{(k)}\left[\mathbf{X}_{(1)}\phi_{(1)}^{\prime}+\cdots+\mathbf{X}_{(d)}\phi_{(d)}^{\prime}+\mathbf{u}\right]$
$\displaystyle=\tilde{\mathbf{A}}_{j}^{(k)}\mathbf{X}_{(1)}\phi_{(1)}^{\prime}+\cdots+\tilde{\mathbf{A}}_{j}^{(k)}\mathbf{X}_{(d)}\phi_{(d)}^{\prime}+\tilde{\mathbf{A}}_{j}^{(k)}\mathbf{u}$
(33)

It can be shown that terms such as
$\tilde{\mathbf{A}}_{j}^{(k)}\mathbf{X}_{(1)}\phi_{(1)}^{\prime}$ will
converge to a constant as $T\rightarrow\infty$, and we only need to consider
the term $\tilde{\mathbf{A}}_{j}^{(k)}\mathbf{u}$ when computing the variance
of the elements in $\hat{\phi}_{(j)}^{(k)\prime}$. In fact, we will derive the
asymptotic distribution of $\hat{\phi}_{(j)}^{(k)\prime}$, which shows the
standard error is asymptotically valid.

We make the following assumptions.

######  Assumption 1.

$y_{t}$ is a VAR($p$) process defined in eq. 1 and is strictly stationary.

######  Assumption 2.

The errors $\left\\{u_{t}\right\\}$ are i.i.d. with $E(u_{t})=0$ and
$\text{Var}(u_{t})=\Omega$. The errors have bounded fourth moment,
$E(u_{i_{1}t}u_{i_{2}t}u_{i_{3}t}u_{i_{4}t})<\infty$ for all
$i_{1},i_{2},i_{3},i_{4}=1,\cdots,d$.

Assumptions 1 and 2 are standard in VAR modeling. In Assumption 1, we
implicitly assume the lag order $p$ is fixed. In practice, $p$ is usually a
small number and can be easily tuned or selected by an information criterion.
The bounded fourth moment condition in Assumption 2 implies that the process
$\\{y_{t}\\}$ itself has a bounded fourth moment, a required condition to
apply the central limit theory for a vector martingale difference sequence. If
one is only interested in applying the LS-Boost algorithm, the i.i.d.
assumption can be relaxed; we only need uncorrelated errors in a VAR($p$)
model and the errors may be heteroskedastic or conditional heteroskedastic.
The i.i.d. assumption are imposed for the convenience of deriving the standard
errors.

Let vec be the matrix vectorization operator and $\otimes$ be the Kronecker
product.

######  Theorem 1.

Under Assumptions 1 and 2, as $T\rightarrow\infty$, for the selected variable
$j$ at boosting step $k\geq 1$,

$\sqrt{T}(\text{vec}(\hat{\phi}_{(j)}^{(k)\prime})-c_{(j)}^{(k)})\rightarrow
N\left(\mathbf{0},\Omega^{(k)}\otimes\mathbf{Q}_{(j)}^{(k)}\right),$ (34)

where $c_{(j)}^{(k)}$, $\Omega^{(k)}$, and $\mathbf{Q}_{(j)}^{(k)}$ are
constants that are unique to the $j$th variable at boosting step $k$.

The explicit expressions of both $c_{(j)}^{(k)}$ and $\mathbf{Q}_{(j)}^{(k)}$
are given in the proof in the online supplement. All quantities such as
$c_{(j)}^{(k)}$, $\Omega^{(k)}$, and $\mathbf{Q}_{(j)}^{(k)}$ are constants
for a given $k$ and they do not depend on the sample data (see the proof in
LABEL:eq:phi_term1, LABEL:eq:phi_term1_proof, LABEL:eq:phi_jk5,
LABEL:eq:phi_jk5_cjk_term and LABEL:eq:Omega_and_Q). Theorem 1 gives a
familiar form of the asymptotic result for a VAR that can be found in standard
textbooks such as Hamilton (1994); Lütkepohl (2005). What is new here are the
expressions for $c_{(j)}^{(k)}$, $\Omega^{(k)}$, and $\mathbf{Q}_{(j)}^{(k)}$.
These constants adapt to the iterative nature of the boosting algorithm, and
$\mathbf{Q}_{(j)}^{(k)}$ keeps accumulating the statistical uncertainty as the
boosting procedure progresses. The diagonal elements of
$\hat{\Omega}^{(k)}\otimes\hat{\mathbf{Q}}_{(j)}^{(k)}$ will allow us to
construct the t statistic and compute the asymptotically valid p-value. In
writing eq. 34, we assume the jth variable is selected as least once up to
step k; otherwise, $\phi_{(j)}^{(k)}$ is not estimated and its value is simply
$\mathbf{0}$.

######  Remark 4.

It is clear from Theorem 1 that, for a given $k$, the boosting estimate
$\hat{\phi}_{(j)}^{(k)}$ is biased and inconsistent since $c_{(j)}^{(k)}\neq
0$. In high-dimensional modeling, working with a biased estimate in practice
is quite common for many adaptive methods such as the lasso and LS-boosting.
The new approach this paper adopts is to use the standard error in eq. 34 to
construct the t-statistic for a biased estimate. Since obtaining a full and
unbiased LS solution is practically impossible, working with a biased estimate
for hypothesis testing provides an alternative. Consequently, the hypothesis
is not for the population parameter
$H_{0}{:}\;\text{vec}(\phi_{(j)}^{\prime})=\boldmath{0}$. Instead, we test the
statistical significance of a LS boosting estimate up to every step $k$. See
Section 5 for further discussions.

######  Remark 5.

In estimation, $\hat{\Omega}^{(k)}$ can be obtained from the cross product of
the residuals ($\hat{\mathbf{u}}^{(k)}$). When the dimension $d$ is large,
this estimate becomes singular and one has to introduce sparsity in
$\Omega^{(k)}$ to get a feasible estimate. There is a large body of research
on estimating high-dimensional variance covariance matrix. Notice that we only
need to estimate the diagonal elements of
$\hat{\Omega}^{(k)}\otimes\hat{\mathbf{Q}}_{(j)}^{(k)}$ for standard errors,
but not the entire matrix. This greatly reduces the number of parameters and
makes our procedure feasible for high-dimensional VARs.

To gain more insights into Theorem 1, let us consider the case at boosting
step $k=1$ when the $j_{1}$th variable is selected ($j=j_{1}$). When $k=1$, we
have

$\displaystyle\hat{\phi}_{(j_{1})}^{(1)\prime}$
$\displaystyle=\nu\mathbf{A}_{j_{1}}^{(1)}\mathbf{Y}$
$\displaystyle=\nu(\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})})^{-1}\mathbf{X}_{(j_{1})}^{\prime}\mathbf{Y}$
$\displaystyle=\nu(\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})})^{-1}\mathbf{X}_{(j_{1})}^{\prime}(\mathbf{X}_{(1)}\phi_{(1)}^{\prime}+\cdots+\mathbf{X}_{(j_{1})}\phi_{(j_{1})}^{\prime}+\cdots+\mathbf{X}_{(d)}\phi_{(d)}^{\prime}+\mathbf{u})$
$\displaystyle=\nu\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}}{T}\right)^{-1}\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(1)}}{T}\right)\phi_{(1)}^{\prime}+\cdots+\nu\phi_{(j_{1})}^{\prime}+\cdots$
$\displaystyle\quad+\nu\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}}{T}\right)^{-1}\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(d)}}{T}\right)\phi_{(d)}^{\prime}+\nu(\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})})^{-1}\mathbf{X}_{(j_{1})}^{\prime}\mathbf{u}.$

Under Assumptions 1 and 2, all cross product terms can be shown to converge to
a constant as $T\rightarrow\infty$, and we have

$\hat{\phi}_{(j_{1})}^{(1)\prime}-c_{(j_{1})}^{(1)}=\nu\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}}{T}\right)^{-1}\frac{1}{T}\sum_{t=1}^{T}\mathbf{X}_{(j_{1}),t}^{\prime}u_{t},$
(35)

where $\mathbf{X}_{(j_{1}),t}^{\prime}$ is the $t$th column of
$\mathbf{X}_{(j_{1})}^{\prime}$, $u_{t}$ is the $t$th row of $\mathbf{u}$, and

$\displaystyle c_{(j_{1})}^{(1)}$
$\displaystyle=\lim\limits_{T\rightarrow\infty}\nu\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}}{T}\right)^{-1}\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(1)}}{T}\right)\phi_{(1)}^{\prime}+\cdots+\nu\phi_{(j_{1})}^{\prime}+\cdots$
$\displaystyle\quad+\nu\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}}{T}\right)^{-1}\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(d)}}{T}\right)\phi_{(d)}^{\prime}.$
(36)

Let the limit of
$\nu\left(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}}{T}\right)^{-1}$
be $c$. Vectorizing eq. 35 gives

$\begin{bmatrix}\hat{\phi}_{(j_{1}),1}^{(1)\prime}-c_{(j_{1}),1}^{(1)}\\\
\vdots\\\ \hat{\phi}_{(j_{1}),i}^{(1)\prime}-c_{(j_{1}),i}^{(1)}\\\ \vdots\\\
\hat{\phi}_{(j_{1}),d}^{(1)\prime}-c_{(j_{1}),d}^{(1)}\end{bmatrix}=\begin{bmatrix}\frac{1}{T}\sum_{t=1}^{T}c\mathbf{X}_{(j_{1}),t}^{\prime}u_{t_{1}}\\\
\vdots\\\
\frac{1}{T}\sum_{t=1}^{T}c\mathbf{X}_{(j_{1}),t}^{\prime}u_{t_{i}}\\\
\vdots\\\
\frac{1}{T}\sum_{t=1}^{T}c\mathbf{X}_{(j_{1}),t}^{\prime}u_{t_{d}}\end{bmatrix}_{dp\times
1},$ (37)

where $\hat{\phi}_{(j_{1}),i}^{(1)\prime}$ and $c_{(j_{1}),i}^{(1)}$ are the
$i$th column of $\hat{\phi}_{(j_{1})}^{(1)\prime}$ and $c_{(j_{1})}^{(1)}$,
respectively; $u_{t_{i}}$ is the $i$th element of the $1\times d$ row vector
$u_{t}$. Equation 37 gives a convenient form to derive the asymptotic
distribution. To compute the (asymptotic) standard error, consider the $i$th
$p\times 1$ vector $\hat{\phi}_{(j_{1}),i}^{(1)\prime}-c_{(j_{1}),i}^{(1)}$ on
the l.h.s. of eq. 37. It can be shown that
$E(u_{t_{i}}|\mathbf{X}_{(j_{1}),t}^{\prime})=0$, because $u_{t_{i}}$ is the
error at time $t$ while $\mathbf{X}_{(j_{1}),t}^{\prime}$ includes $p$ lags
before time $t$.

$\displaystyle\text{Var}(\hat{\phi}_{(j_{1}),i}^{(1)\prime}-c_{(j_{1}),i}^{(1)})$
$\displaystyle=\text{Var}\left(\frac{1}{T}\sum_{t=1}^{T}c\mathbf{X}_{(j_{1}),t}^{\prime}u_{t_{i}}\right)$
$\displaystyle=\frac{1}{T}\text{Var}(c\mathbf{X}_{(j_{1}),t}^{\prime}u_{t_{i}})\qquad\qquad\text{
because $u_{t}$ is i.i.d.}$
$\displaystyle=\frac{1}{T}E(c\mathbf{X}_{(j_{1}),t}^{\prime}u_{t_{i}}^{2}\mathbf{X}_{(j_{1}),t}c^{\prime})\qquad\text{because
$E(u_{t_{i}}|\mathbf{X}_{(j_{1}),t}^{\prime})=0$}$
$\displaystyle=\frac{1}{T}cE(\mathbf{X}_{(j_{1}),t}^{\prime}\mathbf{X}_{(j_{1}),t})c^{\prime}\sigma_{i}^{2},$

where we replace $c$ with
$\nu(\frac{\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}}{T})$,
$E(\mathbf{X}_{(j_{1}),t}^{\prime}\mathbf{X}_{(j_{1}),t})$ with
$\mathbf{X}_{(j_{1})}^{\prime}\mathbf{X}_{(j_{1})}/T$ in estimation.
$\sigma_{i}^{2}$ is the $i$th diagonal element of $\Omega^{(1)}$ and can be
estimated based on the residuals at step $k$.

When $k>1$, the computation for the standard error becomes more involved but
essentially follows the same steps.

######  Remark 6.

From the definition of $c_{(j_{1})}^{(1)}$ in eq. 36, it is clear that the
parameter estimator $\hat{\phi}_{(j_{1})}^{(1)\prime}$ in eq. 35 is biased,
which is expected for an iterative procedure like the LS-Boost. The key result
is, despite of the bias, eq. 35 allows us to compute the asymptotic variance
of $\hat{\phi}_{(j_{1})}^{(1)\prime}$ at step $k=1$, and this asymptotic
variance characterizes the statistical uncertainty in estimating
$\hat{\phi}_{(j_{1})}^{(1)\prime}$.

Next, we use a special case to illustrate the behavior of
$\hat{\phi}_{(j_{1})}^{(1)\prime}$ when $k$ is large. Consider the expression
for $\hat{\phi}_{(j)}^{k\prime}$ in section 3.1. Assume variable $j$ is always
chosen throughout all $k$ steps so that
$\mathbf{A}_{j}^{(1)}=\cdots=\mathbf{A}_{j}^{(k)}=\mathbf{A}_{j}$ and
$\mathbf{H}^{(1)}=\cdots=\mathbf{H}^{(k)}=\mathbf{H}$, and section 3.1 becomes

$\displaystyle\hat{\phi}_{(j)}^{(k)\prime}$
$\displaystyle=\nu\mathbf{A}_{j}\mathbf{Y}+\nu\mathbf{A}_{j}(\mathbf{I}_{T}-\nu\mathbf{H})\mathbf{Y}+\nu\mathbf{A}_{j}(\mathbf{I}_{T}-\nu\mathbf{H})(\mathbf{I}_{T}-\nu\mathbf{H})\mathbf{Y}$
$\displaystyle\quad+\cdots+\nu\mathbf{A}_{j}(\mathbf{I}_{T}-\nu\mathbf{H})\cdots(\mathbf{I}_{T}-\nu\mathbf{H})\mathbf{Y}$
$\displaystyle=\nu\mathbf{A}_{j}\left[(\mathbf{I}_{T}-\nu\mathbf{H})^{0}+(\mathbf{I}_{T}-\nu\mathbf{H})^{1}+\cdots+(\mathbf{I}_{T}-\nu\mathbf{H})^{(k-1)}\right]\mathbf{Y}$
$\displaystyle\rightarrow\nu\mathbf{A}_{j}\left[\mathbf{I}_{T}-(\mathbf{I}_{T}-\nu\mathbf{H})\right]^{-1}\mathbf{Y}\quad\text{
as }k\rightarrow\infty$
$\displaystyle=\mathbf{A}_{j}\mathbf{Y}=(\mathbf{X}_{(j)}^{\prime}\mathbf{X}_{(j)}^{-1}\mathbf{X}_{(j)}^{\prime}\mathbf{Y}.$
(38)

Since $\mathbf{H}$ is idempotent with eigenvalues $0$ and $1$, Theorem 4.3.1
in Horn and Johnson (1985) to show the eigenvalues of
$\mathbf{I}_{T}-\nu\mathbf{H}$ is between $0$ and $1$. Diagonalization gives
$(\mathbf{I}_{T}-\nu\mathbf{H})^{k}=\mathbf{U}\mathbf{D}_{k}\mathbf{U}^{\prime}$,
where $\mathbf{U}$ is orthonormal and $\mathbf{D}_{k}$ is diagonal with all
eigenvalues of $(\mathbf{I}_{T}-\nu\mathbf{H})^{k}\ $, which will allows us to
derive a bias expression similar to that in Proposition 3 and Theorem 1 of
Bühlmann and Yu (2003) for the case of multivariate regression. Since
$\hat{\phi}_{(j)}^{(k)\prime}$ converges to the LS estimator, it is unbiased
as $k\rightarrow\infty$ and will have the same variance as the LS estimator.

In practice, we will have many variables and the boosting algorithm will not
update the same variable at every step so that some (possibly many) of the
$\mathbf{A}_{j}^{(k)}$ in section 3.1 are zero, But it can be roughly seen
that, as long as the boosting step $k$ is large enough, coefficient matrix for
variable $j$ will get many updates. These incremental updates constitute a
subsequence of the geometric matrix series in section 3.1 and will also
converge to the same limit. In Section 4, we provide a different perspective
on the convergence property of the LS-Boost estimator when $k$ is large.

###  3.2 Asymptotic results for LS-Boost2 when $k$ is fixed

Next we discuss the asymptotic results for the estimator described in Section
2.2, and it will be similar to Theorem 1.

By the definition of $\hat{\beta}_{(j)s}^{(k)\prime}$ in eq. 22, we have

$\displaystyle\hat{\phi}_{(j)s}^{(k)\prime}$
$\displaystyle=\hat{\phi}_{(j)s}^{(k-1)\prime}+\nu\hat{\beta}_{(j)s}^{(k)\prime}$
$\displaystyle=\nu\hat{\beta}_{(j)s}^{(1)\prime}+\cdots+\nu\hat{\beta}_{(j)s}^{(k)\prime}$
$\displaystyle=\tilde{\mathbf{A}}_{js}^{(k)}\mathbf{Y},\text{ similar to
}\lx@cref{creftype~refnum}{eq:phi jk2}$ (39)

where

$\tilde{\mathbf{A}}_{js}^{(k)}=\left[\nu\mathbf{A}_{js}^{(1)}+\nu\mathbf{A}_{js}^{(2)}(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)})+\cdots+\nu\mathbf{A}_{js}^{(k)}(\mathbf{I}_{T}-\nu\mathbf{H}^{(k-1)})\cdots(\mathbf{I}_{T}-\nu\mathbf{H}^{(1)})\right].$
(40)

Hence, similar to section 3.1, we have

$\displaystyle\hat{\phi}_{(j)s}^{(k)\prime}$
$\displaystyle=\tilde{\mathbf{A}}_{js}^{(k)}\mathbf{X}_{(1)s}\phi_{(1)s}^{\prime}+\cdots+\tilde{\mathbf{A}}_{js}^{(k)}\mathbf{X}_{(d)s}\phi_{(d)s}^{\prime}+\tilde{\mathbf{A}}_{js}^{(k)}\mathbf{u}.$
(41)

The derivation of the asymptotic distribution result largely follows the proof
of Theorem 1, and we give the result as a corollary in the following.

######  Corollary 1.

Under Assumptions 1 and 2, as $T\rightarrow\infty$, for every selected
variable $j$ and its lag $s$ at every boosting step $k\geq 1$,

$\sqrt{T}(\hat{\phi}_{(j)s}^{(k)}-c_{(j)s}^{(k)})\rightarrow
N\left(\mathbf{0},\Omega^{(k)}\cdot\mathbf{Q}_{(j)s}^{(k)}\right),$ (42)

where $c_{(j)s}^{(k)}$, $\Omega^{(k)}$, and $\mathbf{Q}_{(j)s}^{(k)}$ are
constant that are unique to the $j$th variable and its $s$th lag at boosting
step $k$.

See the supplement for the proof. Note that $\mathbf{Q}_{(j)s}^{(k)}$ is a
scalar in the above corollary. Again we assume the parameter $\phi_{(j)s}$
gets updated at least once between steps $1$ and $k$.

##  4 Convergence of LS-Boost1 estimator for fixed-T

While Theorem 1 and Corollary 1 characterize the behavior of the estimator
when $k$ is fixed and $T\rightarrow\infty$, it will also be useful to study
its behavior when $k\rightarrow\infty$ and $T$ is fixed. At the end of Section
3.1, we briefly discuss the general case of a VAR with many variables and
possibly multiple lags. The LS-boost procedure will likely select a different
variable at each step, and the expression for $\tilde{\mathbf{A}}_{(j)}^{(k)}$
in eq. 32 becomes more complicated. To see the limit of the estimator when
$k\rightarrow\infty$ and $T$ is fixed , we instead derive a computation bound
result for the VAR LS-Boost estimator, similar to the cross-section result in
Freund et al. (2017). The following result is non-asymptotic, but it can
characterize the behavior of the estimator as $k\rightarrow\infty$.

Consider the VAR model in eq. 7. Let
$\lambda_{\text{pmin}}(\mathbf{X}_{g}^{\prime}\mathbf{X}_{g})$ be the smallest
non-zero eigenvalue of $\mathbf{X}_{g}^{\prime}\mathbf{X}_{g}$ and define the
linear convergence rate coefficient

$\gamma=1-\frac{\nu(2-\nu)\lambda_{\text{pmin}}(\mathbf{X}_{g}^{\prime}\mathbf{X}_{g})}{4d}.$
(43)

Freund et al. (2017) shows that $0.75\leq\gamma<1$ when columns of
$\mathbf{X}_{g}$ are normalized with unit $\ell_{2}$ norm. To handle time
series correlation in a VAR, we further assume the data $\mathbf{X}_{(j)}$ are
re-scaled so that its inner product is an identity matrix. This normalization
is only used for the convenience of deriving the theoretical results, and it
is not needed for the actual implementation of the boosting algorithms. More
details are provided in the proof of the following theorem. Let the LS
solution at boosting step $k$ be $\bm{\phi}^{(k)}_{g,\text{LS}}$ and let its
estimator be $\hat{\bm{\phi}}^{(k)}_{g,\text{LS}}$. Also let
$\hat{\bm{\phi}}_{g,\text{LS}}$ be the LS estimator for the VAR model in eq.
7.

######  Theorem 2.

For the LS-Boost1 algorithm described in Section 2.1, under Assumption 1 with
i.i.d. errors, the estimated coefficient bound is given by

$\lVert\hat{\bm{\phi}}_{g}^{(k)}-\hat{\bm{\phi}}^{(k)}_{g,\text{LS}}\rVert_{2}\leq\frac{\lVert\mathbf{X}_{g}\hat{\bm{\phi}}_{\text{g,LS}}\rVert_{2}^{2}}{\lambda_{\text{pmin}}(\mathbf{X}_{g}^{\prime}\mathbf{X}_{g})}\gamma^{k/2}.$
(44)

The prediction bound is given by

$\lVert\mathbf{X}_{g}\hat{\bm{\phi}}_{g}^{(k)}-\mathbf{X}_{g}\hat{\bm{\phi}}_{g,\text{LS}}\rVert_{2}\leq\lVert\bm{\mathrm{X}}\hat{\bm{\phi}}_{\text{g,LS}}\rVert_{2}\gamma^{k/2}.$
(45)

See the supplement for the proof. Theorem 2 is similar to Theorem 2.1 in
Freund et al. (2017). Our results differ from theirs in two aspects: our
results are derived for a VAR and we need to take care of additional
dependence in time series. Similar to the definition of $\bm{\phi}_{g}$ below
eq. 7, $\hat{\bm{\phi}}_{g}^{(k)}$ in eq. 44 is defined as
$\hat{\bm{\phi}}_{g}^{(k)}=\left[\hat{\phi}_{(1)}^{(k)},\cdots,\hat{\phi}_{(d)}^{(k)}\right]^{\prime}$.

Notice that results in eqs. 44 and 45 require very few assumptions to hold.
They are derived for the fixed-$T$ case, and moment conditions in Assumption 2
are not needed here.

Normally, we will simply write $\bm{\phi}^{(k)}_{g}$ instead of
$\bm{\phi}^{(k)}_{g,\text{LS}}$ for the unknown parameter. The subscript “LS”
in $\bm{\phi}^{(k)}_{g,\text{LS}}$ is added to emphasize that, at boosting
step $k$ and corresponding to the specific sparsity at step $k$, there exist
an (or multiple) unknown LS solution $\bm{\phi}^{(k)}_{g,\text{LS}}$. Equation
44 states that $\hat{\bm{\phi}}_{g}^{(k)}$ converges to the closest (in
$\ell_{2}$ norm) LS estimate . Notice that $\bm{\phi}^{(k)}_{g,\text{LS}}$
changes as $k$ changes. Depending on $k$, $\bm{\phi}^{(k)}_{g,\text{LS}}$ may
or may not be unique. When $\bm{\phi}^{(k)}_{g,\text{LS}}$ is non-unique,
$\hat{\bm{\phi}}^{(k)}_{g,\text{LS}}$ estimates one of the LS solutions at
step $k$ and $\hat{\bm{\phi}}_{g}^{(k)}$ converges to one of the many
$\bm{\phi}^{(k)}_{g,\text{LS}}$s. To see why the limit
$\bm{\phi}^{(k)}_{g,\text{LS}}$ may change with $k$, consider the following
two scenarios: When $k$ is small, LS-Boost1 generates very few nonzeros in the
coefficient matrices and a LS solution for these nonzeros is unique for the
model at step $k$; when $k$ is large, LS-Boost1 may have generated a large
number of nonzeros in the coefficient matrices and a LS solution for these
nonzeros exists but is non-unique. In the case of prediction, Equation 45
directly shows that the LS-Boost1 prediction converges to the LS prediction at
a linear rate. Note that, by using an $\alpha$-mixing condition for the data
and several other assumptions, Theorem 2 in Lutz and Bühlmann (2006) obtains a
result for the consistency of LS-Boost prediction when $T$ and $k$ go to
infinity simultaneously, whereas Theorems 1 and 2 in this paper permits the
computation of the s.e. and can characterize precisely the convergence of LS-
Boost prediction as a function of the boosting step $k$.

From Theorem 2, we can immediately deduce the behavior of
$\hat{\bm{\phi}}_{g}^{(k)}$ when $k\rightarrow\infty$ with a fixed $T$, and
Theorem 1 gives the result when $T\rightarrow\infty$ with a fixed $k$. What
would happen if both $k$ and $T\rightarrow\infty$? The result will depend on
the relative speed at which $dp$, $T$ and $k\rightarrow\infty$. Loosely
speaking, if $T>dp$, the estimator from LS-Boost1 will converge to the unique
LS solution and to one of the non-unique LS solutions if $T<dp$ as
$k\rightarrow\infty$. The final result will also be influenced by the
interplay of all these factors and the actual sparsity of the coefficient
matrices.

Our discussion in this section focuses on LS-Boost1 with the understanding the
LS-Boost2 is a special case of LS-Boost1.

##  5 Additional discussions on the p-value

We add some additional discussions on the estimator and its p-value.

Incremental hypothesis testing. It is mentioned in Section 1 that the
incremental hypothesis studied in this paper is not equivalent to the
classical hypothesis for the population regression parameters. Consider the
$k=1$ example discussed in Section 3.1. Theorem 1 implies that when $k=1$, as
$T\rightarrow\infty$, $\text{vec}(\hat{\phi}_{(j)}^{(k)\prime})$ converges to
$c_{(j_{1})}^{(1)}$ in eq. 36 and $c_{(j_{1})}^{(1)}$ can be viewed as the
vectorized version of $\phi_{(j_{1}),\text{LS}}^{(1)}$, the $j_{1}$th block of
$\bm{\phi}_{\text{g,LS}}^{(1)}$ in Theorem 2. Hence, for a fixed boosting step
$k$, the true unknown value of $\hat{\phi}_{(j)}^{(k)\prime}$ is
$\phi_{\text{(j),LS}}^{(k)}$, which is different from the LS solution to eq.
7, $\phi_{(j),\text{LS}}$. At each boosting step, the null hypothesis of the
t-test for each element of the coefficient matrix of the $j$th variable is

$H_{0}{:}\;\phi_{\text{(j),LS},rq}^{(k)}=0,\text{ for }r=1,\cdots,d\text{ and
}q=1,\cdots,p,$ (46)

where $\phi_{\text{(j),LS},rq}^{(k)}$ is the $rq$th element of the matrix
$\phi_{\text{(j),LS}}^{(k)}$, whose transpose,
$\phi_{\text{(j),LS}}^{(k)\prime}$, is the $j$th $p\times d$ block of
$\bm{\phi}_{\text{g,LS}}^{(k)}$, and the value of $\phi_{\text{(j),LS}}^{(k)}$
may change with $k$. The LS-Boost1 estimator $\hat{\phi}_{(j)}^{(k)}$ is an
estimator for $\phi_{\text{(j),LS}}^{(k)}$ and it reflects the cumulative
incremental changes in the $j$th coefficient matrix. Hence, the t-test and its
p-value repeatedly test the null in eq. 46 for each nonzero entry in
$\hat{\phi}_{(j)}^{(k)}$ at every boosting step. If any element in
$\hat{\phi}_{(j)}^{(k)}$ is $0$, we do not have a boosting estimate yet and
cannot perform a t-test for the corresponding element in
$\phi_{\text{(j),LS}}^{(k)}$.

Interpretation of an incremental p-value. Note that
$\phi_{\text{(j),LS},rq}^{(k)}\neq\phi_{\text{(j),LS},rq}$, i.e., the null in
eq. 46 is not about the parameter in the population regression equation. Thus,
when rejecting the null in eq. 46 at $5\%$ significance level, we conclude
that $\phi_{\text{(j),LS},rq}^{(k)}$ is statistically significant but with no
implication on the significance of $\phi_{\text{(j),LS},rq}$. Alternatively,
since the t-statistic is

$t\text{-statistic}=\frac{\hat{\phi}_{(j),rq}^{(k)}-0}{\sqrt{\widehat{\text{var}}(\hat{\phi}_{(j),rq}^{(k)})}},$
(47)

we conclude that, given the model up to step $k$, the cumulative incremental
change, $\hat{\phi}_{(j),rq}^{(k)}$, is statistically different from $0$ when
the null is rejected. The denominator of eq. 47 is obtained from the
asymptotic distribution in Theorem 1. If a parameter estimate is $0.1$ and is
significant at step $k$, when it becomes $0.2$ at step $k+1$, is it also
significant at step $k+1$? Not necessary. Figure 4 shows some examples where
estimates are significant when $k$ is small but become insignificant when $k$
is large. Through repeated hypothesis testing we can track the significance of
a parameter as it evolves along the path of the boosting algorithm. Our s.e.
and p-value are computed by assuming the model up to step $k$ is given and
they do not quantify the statistical uncertainty associated with model
selection at step $k$. Explicitly incorporating the model selection
uncertainty can sometimes lead to wide confidence intervals, see examples for
cross-section regression in Figure 6.12 in Hastie et al. (2015) and Figure 1
in Tibshirani et al. (2016).

Since $\phi_{\text{(j),LS},rq}^{(k)}$ in eq. 47 changes with $k$, it is clear
that our null hypothesis also changes at each boosting step when LS-Boost adds
a new variable in the current step. Hence, the null in our test is not as
fixed as the classical null of $H_{0}{:}\;\phi_{\text{(j)},rq}=0$, but they
are similar: In classical hypothesis testing, given all the variables and a
linear model, we test if $\phi_{\text{(j)},rq}=0$ without questioning where
does the model (or all the variables) come from; in the case of LS-Boost, we
test, given the model at a boosting step, if
$\phi_{\text{(j),LS},rq}^{(k)}=0$.

The superscript $(k)$ makes the null hypothesis in eq. 46 a little unusual.
But it is no more unusual than some of the common practices in statistics. We
draw an analogy between our incremental hypothesis testing and the use of many
shrinkage methods in linear regression. Let $\phi_{(j),rq}$ be a true unknown
parameter in a linear regression, say, VAR. A shrinkage method such as the
lasso or LS-Boost suggests that, even when the LS estimate is available, it is
better not to use it in prediction because a LS solution only minimizes the
error on the training data but not the error on the test data (test error);
instead, use c.v. or AIC to screen the solution path spawn by the algorithm to
get a biased estimate of $\phi_{(j),rq}$, which typically gives better
prediction result and can minimize the test error. It would seem unusual not
to use a LS estimate, but it reflects the practice of bias-variance trade-off
in statistics. Similarly in hypothesis testing, when testing the significance
of a coefficient, we propose not to use the full LS solution — even when such
a LS solution is unique and/or computable. Instead, use the sequential s.e.
and p-value to inspect when a fraction of $\phi_{(j),rq}$ is significant and
when it is not. Hence, it is proper to use a superscript so that it becomes
$\phi_{(j),rq}^{(k)}$. The purpose of sequential hypothesis testing is similar
to the use of c.v. in model selection for the lasso and many other shrinkage
methods: to uncover better solutions in spite of their bias.

Proper use of the p-value. We give two examples of proper use of the proposed
p-value. Example $1$: Given the LS-Boost solution paths of many parameters, an
analyst can choose any boosting step $k$ and use the p-value to inspect the
statistical significance of a selected parameter with a possible Bonferroni
correction. Confidence intervals can be constructed based on the associated
standard errors. Example 2: Given the LS-Boost solution paths of many
parameters, the analyst can use the p-values to remove non-significant
parameters at each boosting step and apply c.v. or validation procedure for
model selection afterwards. Since the p-value is used before c.v., the
statistical uncertainty associated with c.v. will not affect the validity of
the p-value. An incorrect application of the p-value is to use c.v. first to
select a model then apply the p-values to further pare down parameters. It is
incorrect because the p-values are computed at the same time when the boosting
algorithm generates the solution path; using c.v. incurs additional
uncertainty that are not reflected in the p-values so that the p-values cannot
be used unless such additional uncertainty is small.

Trade-off between the false positive rate and the false negative rate. Our
simulation results suggest that we can use a $p$-value to keep the FPR of a
model under control at the cost of increasing the FNR (FPR and FNR are defined
in Section 6.2.1). A data analyst has to decide whether FPR or FNR is of the
most concern. Typically, it is FPR. An equally important topic is proper
control of the false discovery rate (FDR). In a time series model like VAR,
the p-values are dependent not only across variables at each boosting step but
also across boosting steps for each variable. It will be very useful to design
a method to incorporate the standard Benjamini-Hochberg procedure for FDR
control in high-dimensional VARs. This topic is beyond the scope of the
current paper.

A choice between prediction and model interpretation. When a p-value is
supplied, one generally expects it will improve the model. The specific nature
of such improvement depends on the goal of research. Removing a variable based
on its p-value does not necessarily improve prediction. We provide some mixed
results in the simulation exercise. Hence, it is understood that the goal
leans more towards better model interpretation when a p-value is applied. Our
application, though, does provide an example where a p-value-adjusted boosting
model improves both prediction and model interpretation over an unadjusted
boosting model.

##  6 Simulation results

We conduct Monte Carlo simulation in this section to study the use of the
proposed p-values in high-dimensional VAR modeling. We first consider a small
model and compare the boosted p-values to those obtained from the LS method,
and then move onto a larger scale simulation study with high-dimensional VARs.

###  6.1 A simple example of bivariate VAR

Let us consider the following stationary bivariate VAR selected from equation
(4.2.1) in Lütkepohl (2005).

$y_{t}=\begin{bmatrix}0.02\\\ 0.03\end{bmatrix}+\begin{bmatrix}0.5&0.1\\\
0.4&0.5\end{bmatrix}y_{t-1}+\begin{bmatrix}0&0\\\
0.25&0\end{bmatrix}y_{t-2}+u_{t}\text{ and }u_{t}\sim
N\left(\mathbf{0},\begin{bmatrix}0.09&0\\\ 0&0.04\end{bmatrix}\right).$ (48)

We simulate a sample of $500$ observations and report in Table 1 the estimates
and the p-values from the LS method along with those of LS-Boost1 and LS-
Boost2. The learning rate $\nu$ is chosen to be $0.1$. The selected stopping
step at each replication is determined by the bias-corrected Akaike
information criterion (AIC) in Hurvich and Tsai (1989) for the two boosting
methods.

Table 1: Average estimates and p-values for the bivariate VAR model

| $\hat{\phi}_{1}$ | $\hat{\phi}_{2}$ |  p-value for $\hat{\phi}_{1}$ |  p-value for $\hat{\phi}_{2}$
---|---|---|---|---
LS-Boost1 | 0.494 | 0.095 | -0.001 | -0.006 | 0.000 | 0.255 | 0.499 | 0.553
| 0.403 | 0.476 | 0.258 | 0.009 | 0.000 | 0.000 | 0.000 | 0.513
LS-Boost2 | 0.470 | 0.095 | 0.022 | -0.006 | 0.000 | 0.089 | 0.328 | 0.620
| 0.419 | 0.478 | 0.234 | 0.010 | 0.000 | 0.000 | 0.000 | 0.465
LS | 0.495 | 0.092 | 0.002 | -0.005 | 0.000 | 0.288 | 0.499 | 0.554
| 0.403 | 0.490 | 0.247 | 0.005 | 0.000 | 0.000 | 0.000 | 0.521

  * Notes: This table reports estimation output with a sample size of $500$. All numbers are averages over $100$ replications. Results for the intercept are skipped in this table.

In Table 1, the LS-Boost estimates are very close to those of the LS. More
interestingly, the statistical decision based on LS-Boost p-values gives the
identical results compared to the LS if we use $5\%$ as the critical value.
Compared to the result of LS-Boost2, the result of LS-Boost1 is slightly
closer to that of the LS, perhaps due to the fact that, on average, columns
$\mathbf{X}_{g}$ in LS-Boost1 is updated more often since LS-Boost1 selects
all lags of a single variable.

((a)) LS-Boost1 s.e. and LS s.e.  ((b)) Path of LS-Boost1 p-values and the LS
p-value

Figure 1: Figure 1(a) plots the average s.e. of the $8$ coefficient estimators
in eq. 48 of AIC-selected LS-Boost1 models in $100$ replications. The height
of the white bar outlined with a black line is the average s.e. of the
corresponding LS in $100$ replications. Figure 1(b) plots the entire path of
the mean (over $100$ replications) and the median p-value of LS-Boost1. The
black solid and dashed lines are the mean and median p-values of the LS in
$100$ replications.

It is also interesting to observe the boosting p-value can be very close to
that of the LS. For example, the two p-values of LS-Boost1, $0.499$ and
$0.553$, are close to $0.499$ and $0.554$ of the LS. Figure 1(a) plots the
average s.e. of $100$ A.I.C.-selected LS-Boost1 models for the $8$ parameters
in eq. 48 vs. those of the LS. The s.e. of A.I.C.-selected models are close to
those of the LS, in spite of the fact A.I.C. typically stops the boosting
iteration very early in the simulation. Figure 1(b) offers a more complete
picture by plotting the average of LS-Boost1 s.e. for the entire $500$ steps
for the parameter $\phi_{1,12}$. It demonstrates that our p-value can converge
to that of the LS when $k$ is large. This is not a surprise since the sample
size ($500$) is large relative to the number of parameters in eq. 48, the
bivariate VAR has a unique LS solution and Theorem 2 predicts the LS-Boost1
estimator converges to this unique LS solution.

###  6.2 Simulation results for high-dimensional VARs

In this section, we consider high-dimensional VARs and focus on using the
$5\%$ cutoff to remove variables and study its impact on model performance.
This is arguably the most familiar application of a p-value.

####  6.2.1 Simulation setup

We use a VAR(2) model with no intercept in our simulation. Consider the model

$y_{t}=\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+u_{t},$ (49)

where we make the following specifications:

(i) Both $\phi_{1}$ and $\phi_{2}$ are $d\times d$ matrices and are generated
according to a sparsity parameter $s$. Three combinations of sample size
$(T)$, model dimension $(d)$, and coefficient sparsity $(s)$ are considered
(see Table 2 for more details).

Table 2: Sample size, model dimension, and sparsity model type | $T$ | $d$ | $s$ | No. of nonzeros | No. of parameters | sparsity rate
---|---|---|---|---|---|---
|  |  |  | $d\times s\times 2$ | $d\times d\times 2$ |
model type 1 | 200 | 50 | 5 | 500 | 5,000 | 10%
model type 2 | 50 | 100 | 5 | 1,000 | 20,000 | 5%
model type 3 | 100 | 200 | 10 | 4,000 | 80,000 | 5%

For $\phi_{1}$, we randomly select $s$ columns to be nonzero and fill them
with uniform random variables on $[-0.5,0.5]$. The same procedure is used to
generate $\phi_{2}$. To guarantee stationarity of the simulated process, we
need to impose some restriction on the maximum eigenvalue of the coefficient
matrices. Rewrite the VAR(2) in eq. 49 in a VAR(1) format (see equation
10.1.11 in Hamilton (1994))

$\mathcal{Y}_{t}=\mathbf{F}\mathcal{Y}_{t-1}+\mathbf{v},$ (50)

where

$\mathbf{F}=\begin{bmatrix}\phi_{1}&\phi_{2}\\\ \mathbf{I}_{d}&\mathbf{0}\\\
\end{bmatrix}.$ (51)

Given $\phi_{1}$ and $\phi_{2}$, we recursively shrink their value until the
maximum eigenvalue of $\mathbf{F}$ is less than $1$, and this will make
$y_{t}$ in eq. 49 covariance-stationary.

(iii) Create a Toeplitz matrix $\tilde{\Omega}$, whose $ij$th element is
$\rho^{|i-j|}$ with $\rho$ chosen to be $0.5$. For a given signal-to-noise
(SNR) level, choose a value for $\sigma^{2}$ so that it meets the desired SNR
level, i.e.,
$\text{SNR}=\frac{\lambda_{\text{max}}(\mathbf{F})}{\sigma^{2}\lambda_{\text{max}}(\tilde{\Omega})}$.
Let the error variance-covariance matrix be $\Omega=\sigma^{2}\tilde{\Omega}$.
Draw the i.i.d. $d\times 1$ error vector $u_{t}$ from $N(\mathbf{0},\Omega)$.
We consider three SNR levels in our study: $0.5$, $1.0$, and $3.0$.

(iv) With $\phi_{1}$, $\phi_{2}$, and $\Omega$, we can simulate the process
according to eq. 49.

(v) Run various methods (discussed below), including LS-Boost1 and LS-Boost2,
on the simulated data. For each method, we estimate the model based on a
training data set, choose the tuning parameter on a validation set, and
evaluate the model performance on a test set. The number of observations in
both validation and test sets is $200$. All model metrics reported are
averages of $100$ replications. Alternatively, one can use an information
criterion to determine the tuning parameter. Since there are quite a few
information criteria to choose from, using a validation set provides a simple
and uniform approach for model selection across different methods.

(vi) For each method, we report the following metrics: Mean squared error
(MSE) for estimated coefficients:
$\text{MSE}=\lVert\hat{\bm{\phi}}-\bm{\phi}\rVert_{2}^{2}/(p\times d^{2})$;
Mean squared prediction error (MSPE) on the test set:
$\text{MSPE}=\lVert\mathbf{X}(\hat{\bm{\phi}}-\bm{\phi})\rVert_{2}^{2}/(T_{\text{test}}\times
d)$, where $T_{\text{test}}=200$; false positive rate (FPR):
$\text{FPR}=\text{FP}/(\text{TN}+\text{FP})$ and false negative rate (FNR):
$\text{FNR}=\text{FN}/(\text{TP}+\text{FN})$, where TN, FP, TP, and FN are
number of true zeros, false nonzeros, true nonzeros, and false nonzeros;
F-score: $\text{F-score}=2\text{TP}/(2\text{TP}+\text{FP}+\text{FN})$. The
F-score measures the accuracy of support recovery, and it balances FPR and
FNR. Finally, we also report the model size by counting the number of nonzeros
in the selected model from validation for each method. All evaluation metrics
are averages over $100$ replications.

####  6.2.2 Methods under consideration

We mentioned in Section 1 that there is now a long list of statistical methods
available for estimation in high-dimensional VARs. Our goal is not an
extensive comparison of different methods but the study of estimation and use
of the proposed p-value in LS-Boost. Hence, we select only a few of the other
methods for demonstration purposes.

The first method to consider is the lasso. Various papers have used the lasso
in high-dimensional VARs, and the actual implementation of this technique
varies; variants such as the adaptive lasso or the group lasso can also be
used. One common approach is to apply the lasso to each column of $\mathbf{Y}$
in eq. 3 separately. However, in a VAR setup, it is important to do estimation
jointly to incorporate possible correlation across the $d$ variables. We
instead use the multi-response approach in the R package glmnet by setting the
option family = "mgaussian". For the $j$th column in $\mathbf{X}$ in eq. 3,
the algorithm uses a group lasso penalty to penalize the entire $j$th row of
$\bm{\phi}$. Hence, the modeling strategy is: a column of $\mathbf{X}$ will be
included in the model for all $d$ response variables or excluded for all $d$
response variables. This is exactly the same modeling strategy LS-Boost2
adopts. We note that a row in $\bm{\phi}$ is a column in either $\phi_{1}$ or
$\phi_{2}$. Given that we simulate $\phi_{1}$ and $\phi_{2}$ by randomly
assign $s$ nonzero columns, this way of applying the lasso to the simulated
data will give the lasso a slight edge over LS-Boost1.

We also consider two other recent methods that are suitable for estimation in
high-dimensional VARs: the sparse orthogonal factor regression (SOFAR) in
Uematsu et al. (2019) and the sparse reduced-rank regression (SRRR) in Chen
and Huang (2012), both of which are implemented in the R package rrpack. Both
method requires a prespecified rank parameter. We tune this nrank parameter on
the grid of $\\{1,2,3,4,5\\}$. These two methods also uses the lasso to
penalize certain aspects of the coefficient matrix $\bm{\phi}$ to generate
sparsity, and we use $100$ for the length of the lasso tuning sequence. All
other model control parameters are the default values provided by the rrpack
package. Both methods are very sophisticated and flexible in application. For
example, the sofar function in the package has more then $10$ model control
parameters. It is almost impossible to experiment with all model parameter
configurations. We use the default values for simplicity purposes.

Since we use $5$ rank parameters and a sequence of $100$ penalty parameters in
both the SOFAR and the SRRR, we set the length of tuning sequence in the lasso
and LS-Boost to $500$. In both LS-Boost1 and LS-Boost2, we use the learning
rate $0.1$, and train the model up to $500$ iterations and send these
solutions to a validation set for model selection and a test set for model
evaluation.

Finally, after getting the boosting solution paths and the p-values, we use
$5\%$ as the p-value cutoff to remove nonzeros in the boosting estimates and
call these estimates LS-Boost1p and LS-Boost2p, respectively. Apply LS-Boost1p
and LS-Boost2p to the validation set for model selection and to the test set
for model evaluation. As it is discussed in Section 5, this is a valid way to
use the proposed p-value.

####  6.2.3 The simulation results and the evaluation metrics

Tables 3, 4 and 5 report the simulation results for the three models described
in Table 2 at three SNR levels. Our major conclusions are the proposed p-value
can help LS-Boost improve FPR and F score, and it can also help LS-Boost
significantly reduce model size without compromising MSE or MSPE too much.

FPR, FNR and F score. From these tables, it is clear that all regularized
methods under consideration in high dimensional VAR can have very high FPR,
implying that their Type I error rate would be too high to make these model
useful for interpretation. After applying the p-value to remove many nonzeros,
the FPR for LS-Boost1p and LS-Boost2p is roughly around $5\%$, bringing the
FPR under control. A downside is the FNR increases as the FPR decreases.
Inspecting the F score will be helpful. Overall, the p-value-adjusted models
have a higher F score compared to that of LS-Boost1 and LS-Boost2 and the
other two methods.

Model size, MSE and MSPE. The MSE and MSPE of all methods are more or less on
the same scale. When the model dimension is very large in Table 5, the lasso
and LS-Boost start to outperform the SOFAR and SRRR by a small margin.
Notably, LS-Boost1p and LS-Boost2p usually deliver similar MSE and MSPE with
far fewer parameters, as can be seen from the Model size column in all tables.
If we compare the MSPE of LS-Boost1 and LS-Boost1p, LS-Boost1p can improve
MSPE in Table 3 and its performance in Table 5 is a little worse. Overall,
based on these tables, we cannot conclude that applying p-values to the
original LS-Boost models will always improve MSE or MSPE. However, it should
also be noted that estimators such as LS-Boost1p and LS-Boost2p are lot more
effective in using nonzero parameters. Consider the MSPE in Table 5 when
$\text{SNR}=3.0$, LS-Boost2 uses $12,832$ parameters to deliver an MSPE of
$0.0823$ and the lasso uses $28,050$ parameters to obtain an MSPE of $0.0813$.
But it takes only $4,501$ parameters for LS-Boost2p to yield a similar MSPE of
$0.0862$. Similar pattern can be observed throughout all tables and also in
our application example.

Table 3: Evaluation metrics for model type $1$ with $T=200,d=50,s=5$

Estimator | MSE | MSPE | FPR | FNR | F score | Model size
---|---|---|---|---|---|---
| $\text{SNR}=0.5$
SOFAR | 0.0039 | 0.3907 | 0.262 | 0.380 | 0.320 | 1487.7
SRRR | 0.0043 | 0.4064 | 0.985 | 0.004 | 0.183 | 4929.5
Lasso | 0.0040 | 0.3939 | 0.423 | 0.274 | 0.264 | 2266.0
LS-Boost1 | 0.0043 | 0.4167 | 0.586 | 0.000 | 0.282 | 3135.0
LS-Boost2 | 0.0041 | 0.3975 | 0.321 | 0.324 | 0.305 | 1780.5
LS-Boost1p | 0.0042 | 0.4081 | 0.064 | 0.575 | 0.424 | 501.9
LS-Boost2p | 0.0041 | 0.3981 | 0.049 | 0.585 | 0.448 | 427.3
| $\text{SNR}=1.0$
SOFAR | 0.0039 | 0.1953 | 0.259 | 0.382 | 0.320 | 1476.2
SRRR | 0.0040 | 0.1958 | 0.849 | 0.075 | 0.194 | 4284.0
Lasso | 0.0040 | 0.1970 | 0.423 | 0.274 | 0.264 | 2266.0
LS-Boost1 | 0.0043 | 0.2084 | 0.586 | 0.000 | 0.282 | 3135.0
LS-Boost2 | 0.0041 | 0.1988 | 0.321 | 0.324 | 0.305 | 1780.5
LS-Boost1p | 0.0042 | 0.2041 | 0.064 | 0.575 | 0.424 | 501.9
LS-Boost2p | 0.0041 | 0.1990 | 0.049 | 0.585 | 0.448 | 427.3
| $\text{SNR}=3.0$
SOFAR | 0.0039 | 0.0651 | 0.243 | 0.391 | 0.332 | 1396.9
SRRR | 0.0038 | 0.0637 | 0.225 | 0.404 | 0.335 | 1309.0
Lasso | 0.0040 | 0.0657 | 0.423 | 0.282 | 0.261 | 2263.0
LS-Boost1 | 0.0043 | 0.0695 | 0.584 | 0.000 | 0.282 | 3126.0
LS-Boost2 | 0.0041 | 0.0663 | 0.315 | 0.333 | 0.304 | 1753.0
LS-Boost1p | 0.0042 | 0.0680 | 0.065 | 0.572 | 0.426 | 505.3
LS-Boost2p | 0.0041 | 0.0664 | 0.050 | 0.582 | 0.448 | 434.1

  * Notes: All evaluation metrics reported in the table are averages over $100$ replications. In model type 1, the true nonzeros are $500$ and the total number of parameters is $5,000$. See Section 6.2.1 for the computation of MSE, MSPE, FPR, FNR, and F score. Model size is the total number of nonzeros in the selected model for each method.

Table 4: Evaluation metrics for model type $2$ with $T=50,d=100,s=5$

Estimator | MSE | MSPE | FPR | FNR | F score | Model size
---|---|---|---|---|---|---
| $\text{SNR}=0.5$
SOFAR | 0.0028 | 0.8827 | 0.315 | 0.372 | 0.169 | 6609.9
SRRR | 0.0025 | 0.8428 | 0.495 | 0.250 | 0.135 | 10156.0
Lasso | 0.0026 | 0.8553 | 0.287 | 0.398 | 0.173 | 6046.0
LS-Boost1 | 0.0031 | 0.9559 | 0.313 | 0.000 | 0.260 | 6952.0
LS-Boost2 | 0.0028 | 0.8802 | 0.183 | 0.446 | 0.230 | 4035.0
LS-Boost1p | 0.0031 | 0.9523 | 0.059 | 0.621 | 0.304 | 1502.1
LS-Boost2p | 0.0029 | 0.9000 | 0.040 | 0.666 | 0.320 | 1097.0
| $\text{SNR}=1.0$
SOFAR | 0.0028 | 0.4414 | 0.323 | 0.364 | 0.168 | 6777.4
SRRR | 0.0025 | 0.4145 | 0.426 | 0.280 | 0.147 | 8812.0
Lasso | 0.0026 | 0.4277 | 0.286 | 0.398 | 0.173 | 6043.0
LS-Boost1 | 0.0031 | 0.4780 | 0.313 | 0.000 | 0.260 | 6952.0
LS-Boost2 | 0.0028 | 0.4401 | 0.183 | 0.446 | 0.230 | 4035.0
LS-Boost1p | 0.0031 | 0.4761 | 0.059 | 0.621 | 0.304 | 1502.1
LS-Boost2p | 0.0029 | 0.4500 | 0.040 | 0.666 | 0.320 | 1097.0
| $\text{SNR}=3.0$
SOFAR | 0.0029 | 0.1466 | 0.334 | 0.364 | 0.163 | 6981.1
SRRR | 0.0023 | 0.1329 | 0.195 | 0.387 | 0.235 | 4313.0
Lasso | 0.0026 | 0.1415 | 0.290 | 0.411 | 0.168 | 6106.0
LS-Boost1 | 0.0031 | 0.1578 | 0.332 | 0.000 | 0.250 | 7300.0
LS-Boost2 | 0.0028 | 0.1458 | 0.182 | 0.450 | 0.228 | 4009.0
LS-Boost1p | 0.0031 | 0.1579 | 0.058 | 0.623 | 0.305 | 1478.1
LS-Boost2p | 0.0029 | 0.1490 | 0.040 | 0.667 | 0.320 | 1092.4

  * Notes: All evaluation metrics reported in the table are averages over $100$ replications. In model type 2, the true nonzeros are $1,000$ and the total number of parameters is $20,000$. See Section 6.2.1 for the computation of MSE, MSPE, FPR, FNR, and F score. Model size is the total number of nonzeros in the selected model for each method.

Table 5: Evaluation metrics for model type $3$ with $T=100,d=200,s=10$

Estimator | MSE | MSPE | FPR | FNR | F score | Model size
---|---|---|---|---|---|---
| $\text{SNR}=0.5$
SOFAR | 0.0015 | 0.5919 | 0.242 | 0.423 | 0.190 | 20678.7
SRRR | 0.0015 | 0.5865 | 0.546 | 0.224 | 0.128 | 44628.0
Lasso | 0.0012 | 0.4889 | 0.331 | 0.349 | 0.164 | 27786.0
LS-Boost1 | 0.0013 | 0.5476 | 0.254 | 0.000 | 0.307 | 23292.0
LS-Boost2 | 0.0012 | 0.4951 | 0.136 | 0.446 | 0.279 | 12590.0
LS-Boost1p | 0.0014 | 0.5532 | 0.063 | 0.668 | 0.263 | 6137.0
LS-Boost2p | 0.0013 | 0.5181 | 0.042 | 0.701 | 0.286 | 4404.7
| $\text{SNR}=1.0$
SOFAR | 0.0015 | 0.2960 | 0.243 | 0.422 | 0.189 | 20778.0
SRRR | 0.0014 | 0.2911 | 0.503 | 0.244 | 0.134 | 41238.0
Lasso | 0.0012 | 0.2444 | 0.331 | 0.349 | 0.164 | 27786.0
LS-Boost1 | 0.0013 | 0.2738 | 0.254 | 0.000 | 0.307 | 23292.0
LS-Boost2 | 0.0012 | 0.2476 | 0.136 | 0.446 | 0.279 | 12590.0
LS-Boost1p | 0.0014 | 0.2766 | 0.063 | 0.668 | 0.263 | 6137.0
LS-Boost2p | 0.0013 | 0.2591 | 0.042 | 0.701 | 0.286 | 4404.7
| $\text{SNR}=3.0$
SOFAR | 0.0015 | 0.0984 | 0.256 | 0.420 | 0.183 | 21777.4
SRRR | 0.0014 | 0.0950 | 0.334 | 0.340 | 0.165 | 28030.0
Lasso | 0.0012 | 0.0813 | 0.335 | 0.344 | 0.164 | 28050.0
LS-Boost1 | 0.0013 | 0.0911 | 0.254 | 0.000 | 0.310 | 23340.0
LS-Boost2 | 0.0012 | 0.0823 | 0.139 | 0.440 | 0.278 | 12832.0
LS-Boost1p | 0.0014 | 0.0919 | 0.064 | 0.668 | 0.263 | 6159.0
LS-Boost2p | 0.0013 | 0.0862 | 0.043 | 0.700 | 0.284 | 4501.2

  * Notes: All evaluation metrics reported in the table are averages over $100$ replications. In model type 3, the true nonzeros are $4,000$ and the total number of parameters is $80,000$. See Section 6.2.1 for the computation of MSE, MSPE, FPR, FNR, and F score. Model size is the total number of nonzeros in the selected model for each method.

##  7 Application

In this section, we apply the proposed p-value to the VAR modeling of the
monthly macroeconomic data set provided in McCracken and Ng (2016). This
dataset is frequently updated and well-maintained, and we use version 2022-02
that records $128$ economic variables from $1959/01$ to $2022/01$. After
removing missing observations and performing the suggested transformation such
as difference and twice difference as described by the TCODE variable in the
paper, we end up with $104$ variables and $758$ observations. Out of these
$758$ observations, $50\%$, $25\%$ and $25\%$ are used for training,
validation and test, respectively. An intercept is included in all model
fitting.

Exactly the same tuning process used in the simulation section is applied to
each method in this application except that for SOFAR and SRRR. We increase
the penalty parameter sequence for these two methods from $100$ to $500$ while
still tuning over five rank values, giving these two methods a large grid for
the search of a good solution. SOFAR sometimes does not identify any model on
the validation set and we will use the model selected by AIC instead. Both
SOFAR and SRRR have many model parameters to tune in the R package rrpack. We
will not make further effort to tune these parameters in their R code while
keeping in mind that their predication performance in our reported table can
be likely improved if more tuning is implemented. We focus on the proposed
p-value and its impact on the LS-Boost procedure.

Table 6: MSPE and model size in application

| SOFAR | SRRR | Lasso | LS-Boost1 | LS-Boost2 | LS-Boost1p | LS-Boost2p
---|---|---|---|---|---|---|---
| MSPE
VAR(1) | 3.399 | 0.176 | 270.718 | 391.494 | 391.494 | 0.157 | 0.157
VAR(2) | 6.750 | 0.190 | 182.781 | 255.170 | 354.143 | 0.159 | 0.157
VAR(3) | 3.247 | 0.211 | 105.392 | 230.217 | 200.226 | 0.168 | 0.158
VAR(4) | 3.788 | 0.653 | 87.728 | 187.724 | 213.099 | 0.171 | 0.159
| Model size (No. of nonzero parameter estimates)
VAR(1) | 528 | 1144 | 6240 | 7072 | 7072 | 411 | 411
VAR(2) | 968 | 1352 | 8216 | 14352 | 10296 | 425 | 410
VAR(3) | 32448 | 1560 | 6656 | 15288 | 7488 | 519 | 411
VAR(4) | 43264 | 2184 | 7904 | 20384 | 10608 | 618 | 410

  * Notes: The sample size for the training, validation and test sets are $378,189$ and $189$. The MSPE is computed for the test set. Intercept is excluded when computing the model size of each method. The total number of parameters, excluding the intercept, for VAR(1), VAR(2), VAR(3), and VAR(4) are $10816,21632,32448$, and $43264$, respectively.

Table 6 reports the main results of the application. Two model evaluation
metrics are included: MSPE and model size. The MSPE is computed based on the
test set and model size is the number of nonzero parameter estimates for the
model selected from the validation data set for each method. The MSPE shows
that both Lasso and LS-Boost fail to produce a result comparable to that of
SOFAR or SRRR. This should not be construed as a failure of the lasso method
in general. Our earlier discussion points out that there is no unique way to
implement the lasso in VAR. Other forms of the lasso can also be implemented.
We discuss a few aspects of Table 6 in the following.

###  7.1 Model sparsity

One reason that LS-Boost fails to deliver any meaningful result is probably
due to the noise in the data. Noise and correlation make it very difficult for
the boosting algorithm to focus on the most essential subset of variables;
instead, the algorithm keeps adding new variables to the model. The model size
quickly becomes out of control with the number of nonzeros running well into
thousands or tens of thousands. For example, in a VAR(4) model, LS-Boost1
gives a model with $2,0384$ nonzero parameter estimates. Given the VAR system
have $104$ equations, it implies, on average, each equation has $196$ nonzero
parameters on the r.h.s. The possible high FPR (and possibly high FDR)
translates to a very poor performance on the test set. Using the p-value
salvages the boosting method. By simply removing a coefficient estimate with a
p-value larger than $5\%$, we significantly reduces model size. To appreciate
the reduction in model size, let us consider the VAR(4) model. A VAR(4) has a
total of $43,264(=104\times 104\times 4)$ parameters, excluding the intercept.
The algorithm LS-Boost2p gives only $410$ nonzeros, achieving a sparsity of
about $0.9\%$, reducing the average number of variables per equation from
$196$ to about $4$, making the selected model highly interpretable.

###  7.2 Model stability

A surprising benefit is p-values can help control model stability. Use the LS-
Boost2p as an example. As the model changes from VAR(1) to VAR(4), the model
size barely changes, indicating that the $5\%$ critical value removes most of
the nonzeros. In fact, for LS-Boost2p, all the nonzeros appear in the
coefficient matrix for $\mathbf{Y}_{t-1}$ and all coefficients in the second,
third, and fourth lags are zero. Figures 2(a) and 2(b) give an illustration of
the VAR(1) coefficients. Variables appear in Figures 2(a) and 2(b) are
described in Table 7.

Table 7: Variable definitions in Figures 2(a), 2(b) and 3

Variable | Group No. | Description
---|---|---
AAAFFM | 6 I.R.E.R. | Moody’s Aaa corporate bond minus federal funds rate
BAAFFM | 6 I.R.E.R. | Moody’s Baa corporate bond minus federal funds rate
GS1 | 6 I.R.E.R. | 1-year treasury rate
T10YFMM | 6 I.R.E.R. | 10-year treasury constant maturity minus federal funds rate
T1YFMM | 6 I.R.E.R. | 1-year treasury constant maturity minus federal funds rate
T5YFMM | 6 I.R.E.R. | 5-year treasury constant maturity minus federal funds rate
TB3SMFMM | 6 I.R.E.R. | 3-month treasury constant maturity minus federal funds rate
TB6SMFMM | 6 I.R.E.R. | 6-month treasury constant maturity minus federal funds rate

  * Notes: All variables are from Group $6$ described in McCracken and Ng (2016). I.R.E.R. is the abbreviation for the group name Interest Rate and Exchange Rates.

((a)) Negative links between $\mathbf{Y}_{t}$ and $\mathbf{Y}_{t-1}$ ((b))
Positive links between $\mathbf{Y}_{t}$ and $\mathbf{Y}_{t-1}$

Figure 2: Figures 2(a) and 2(b) describe the negative and positive links in
the VAR(1) coefficient matrix from LS-Boost2p, respectively. In each
subfigure, the right half circle includes variables in $\mathbf{Y}_{t-1}$ with
a nonzero (negative or positive) coefficient on $\mathbf{Y}_{t}$; the left
half circle includes variables in $\mathbf{Y}_{t}$ that the $6$ variables
affect. Green and red links indicate negative and positive coefficient
estimates, respectively. There are a total of $411$ links, equivalent to $411$
nonzeros reported in Table 6 for LS-Boost2p.

Out of the $104$ variables, $6$ of them are consistently selected by the
boosting procedure. It is important to recall that, for LS-Boost2p, all
nonzeros concentrate in the coefficient for $\mathbf{Y}_{t-1}$. Hence, these
$6$ variables are the only (statistically significant) source that drives
$\mathbf{Y}_{t}$, not only in VAR(1) but also in higher-order VARs (at least
up to oder $4$). Consider Figure 2(a) that plots only the negative
coefficients in the VAR(1) model based on LS-Boost2p. The green color for the
links indicates a negative coefficient and all links start from the right half
circle to the left half circle with an arrow pointing to the left, mimicking
the usual setup in a regression where $\mathbf{Y}_{t-1}$ appears on the r.h.s.
and $\mathbf{Y}_{t}$ on the l.h.s. The width of a link is proportional to the
absolute value of the corresponding coefficient estimate, and a wider link
indicates a negative coefficient with a larger absolute value. The colored
sectors have different width. A wider sector on the right half circle suggests
a variable has a larger overall influence (sum of absolute values of the
coefficients). For example, AAAFFM has the widest sector due to its overall
large coefficients. The sector for AAAFFM is only partially filled in Figure
2(a) because most of the corresponding coefficients are positive; the sector
will be completely filled if we combine the links in Figures 2(a) and 2(b).

The links reveals various properties in the data. Consider the GS1 sector in
Figure 2(a) as an example. GS1 is the 1-year treasury rate, a key variable for
short-term interest rate. There are four relatively wide arrows leaving the
sector GS1 and pointing at four variables in $\mathbf{Y}_{t}$. These four
negative links may puzzle a researcher at a first look, but they are largely
due to the way variables are defined in McCracken and Ng (2016). For example,
we see GS1 has a relatively large, negative impact on AAAFFM, the corporate
bond yield spread. We would expect that bond yield will increase as the short-
term interest rate increases. A negative link is counter-intuitive. Notice
that AAAFFM is defined as corporate bond yield minus the federal funds rate
(FEDFUNDS). Using the $758$ observations in the data, we add FEDFUNDS back to
AAAFFM and compute its correlation with GS1, and the result is $0.899$,
indicating these two time series are positively correlated. Hence, it is
important to interpret the result between GS1 and AAAFFM in Figure 2(a) as the
negative relationship between the short-term interest rate and corporate bond
yield spread. Similar explanations applies to other links between GS1 and
other bond yields or long-term interest rates in that figure.

Figure 2(b) describes all positive links between the $6$ variables and
$\mathbf{Y}_{t}$. We observe that most of the $411$ coefficient estimates from
LS-Boost2p are positive. The $6$ variables have relatively large
autoregressive coefficients, but there are also many positive links to other
variables. In a VAR(1) model with $104$ variables, a dense coefficient matrix
will have $10,816(=104\times 104)$ links. Hence Figure 2(b) is already a very
sparse representation of the coefficient matrix. To further help interpret the
result, we remove all coefficients that are $<0.1$ and re-draw the links in
Figure 3.

Figure 3: Large positive links with coefficient estimates $\geq 0.1$ in Figure
2(b).

Figure 3 is now even more interpretable. It clearly reveals that a risk
premium such as AAAFFM, BAAFFM, and T10YFFM has a highly persistent pattern
with its own lag, and their impact on other variables is also obvious. The
dynamic feedback mechanism is asymmetric in a VAR. For example, T10YFFM has a
positive effect on T1YFFM but not vice versa. Although we selectively choose
positive coefficients $\geq 0.1$, one should also take into consideration the
scale of the data when measuring the impact of a lag variable.

Figures 2(a), 2(b) and 3 illustrate many important economic relationships
uncovered from high-dimensional VAR analysis. The p-value approach
significantly shrinks the model size and keep the model stable in high
dimension. It provides a valuable tool for empirical researcher to investigate
large data sets.

###  7.3 Sequential testing

Section 5 points out the sequential nature of our testing procedure. We
further illustrate it with a full plot of the $500$ p-values for each of the
$6$ variables shown in Figures 2(a), 2(b) and 3. Our empirical example has
$104$ equations. We select the equation for civilian unemployment rate
(UNRATE) as an example.

Figure 4: Figure 4 plots the entire path of the p-values for the $6$ variables
in the equation for Civilians Unemployment Rate (UNRATE). Use the right scale
to read the p-values for BAAFFM and TB6SMFFM. The validation step chooses to
stop at step $17$, marked by a vertical dotted line.

Figure 4 plots the path of p-values for all $6$ variables. The validation step
selects the model at step $17$. Three variables, AAAFFM, T10YFFM and TB3SMFFM
have very small p-values throughout the $500$ boosting steps. The p-values for
GS1 is also small, but climbed above the $5\%$ cutoff after around step $50$,
but is below $5\%$ at step $17$. The p-values for BAAFFM and TB6SMFFM,
however, become very large in later steps and can exceed $50\%$. Had the
validation procedure stopped at, say, step $300$, both variables would have
been excluded from the model.

No Bonferroni correction is used at each step when excluding variables. The
Bonferroni correction will give very conservative results and will likely
yield an even sparser model. Whether one should use such adjustment is still
an open question in practice.

##  8 Conclusions

This paper proposes a p-value for the LS-Boost algorithm in high-dimensional
VARs that adapts to its iterative nature. Our p-value can be used for
incremental hypothesis testing that tests the statistical significance of
every selected parameter at each boosting step. We derive the asymptotic
distribution for the estimator for a given boosting step and also discuss its
convergence when the boosting step is large. Simulation results indicate a
p-value-adjusted model can improve the FPR and F score of a model produced by
LS-Boost. The application further reveals that the proposed p-value can help
control both model size and stability, making model interpretation relatively
easy in high dimensions.

LS-Boost is a classical regression algorithm and our results complement its
wide use in practice. The proposed p-value is conceptually simple and
technically practical, and our R package boostvar also makes it
computationally easy. In addition, when both sample size and boosting step are
large, we show the p-value can converge to that of LS estimator in both VAR
and cross-section regression (see tables in the online supplement), providing
an implicit computational guarantee for the p-value.

It will be interesting to design a procedure to properly control the FDR based
on the sequential p-values in a VAR. In addition, more applications in both
cross-section regression and VAR are needed to further study the property of
the p-value. We leave these topics for future work.

## Acknowledgments

The author thanks the Department of Economics, Finance, and Quantitative
Analysis for financial support and the Office of Research at Kennesaw State
University for computation support.

## References

  * Basu and Michailidis (2015) Basu, S. and Michailidis, G. (2015). Regularized estimation in sparse high-dimensional time series models. The Annals of Statistics, 43 (4), 1535 – 1567.
  * Bühlmann and Yu (2003) Bühlmann, P. and Yu, B. (2003). Boosting with the l2 loss. Journal of the American Statistical Association, 98 (462), 324–339.
  * Chen and Huang (2012) Chen, L. and Huang, J. Z. (2012). Sparse reduced-rank regression for simultaneous dimension reduction and variable selection. Journal of the American Statistical Association, 107 (500), 1533–1545.
  * Freund et al. (2017) Freund, R. M., Grigas, P. and Mazumder, R. (2017). A new perspective on boosting in linear regression via subgradient optimization and relatives. The Annals of Statistics, 45 (6), 2328–2364.
  * Friedman (2001) Friedman, J. H. (2001). Greedy function approximation: A gradient boosting machine. The Annals of Statistics, 29 (5), 1189–1232.
  * Fujita et al. (2007) Fujita, A., Sato, J. R., Garay-Malpartida, H. M., Yamaguchi, R., Miyano, S., Sogayar, M. C. and Ferreira, C. E. (2007). Modeling gene expression regulatory networks with the sparse vector autoregressive model. BMC Systems Biology, 1 (1), 1–11.
  * Hamilton (1994) Hamilton, J. (1994). Time Series Analysis. Princeton University Press.
  * Han et al. (2015) Han, F., Lu, H. and Liu, H. (2015). A direct estimation of high dimensional stationary vector autoregressions. Journal of Machine Learning Research, 16 (97), 3115–3150.
  * Hastie et al. (2015) Hastie, T., Tibshirani, R. and Wainwright, M. (2015). Statistical Learning with Sparcity: The Lasso and Generalizations. Chapman and Hall/CRC.
  * Horn and Johnson (1985) Horn, R. A. and Johnson, C. R. (1985). Matrix Analysis. Cambridge University Press.
  * Hurvich and Tsai (1989) Hurvich, C. M. and Tsai, C.-L. (1989). Regression and time series model selection in small samples. Biometrika, 76 (2), 297–307.
  * Lee et al. (2016) Lee, J. D., Sun, D. L., Sun, Y. and Taylor, J. E. (2016). Exact post-selection inference, with application to the lasso. The Annals of Statistics, 44 (3), 907 – 927.
  * Lockhart et al. (2014) Lockhart, R., Taylor, J., Tibshirani, R. J. and Tibshirani, R. (2014). A significance test for the lasso. The Annals of Statistics, 42 (2), 413 – 468.
  * Lozano et al. (2009) Lozano, A. C., Abe, N., Liu, Y. and Rosset, S. (2009). Grouped graphical granger modeling for gene expression regulatory networks discovery. Bioinformatics, 25 (12), 110–118.
  * Lütkepohl (2005) Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer Berlin, Heidelberg.
  * Lutz and Bühlmann (2006) Lutz, R. W. and Bühlmann, P. (2006). Boosting for high-multivariate responses in high-dimensional linear regression. Statistica Sinica, 16 (2), 471–494.
  * McCracken and Ng (2016) McCracken, M. W. and Ng, S. (2016). Fred-md: A monthly database for macroeconomic research. Journal of Business & Economic Statistics, 34 (4), 574–589.
  * Medeiros and Mendes (2016) Medeiros, M. C. and Mendes, E. F. (2016). L1-regularization of high-dimensional time-series models with non-gaussian and heteroskedastic errors. Journal of Econometrics, 191 (1), 255–271.
  * Meinshausen et al. (2009) Meinshausen, N., Meier, L. and Bühlmann, P. (2009). p-values for high-dimensional regression. Journal of the American Statistical Association, 104 (488), 1671–1681.
  * Nicholson et al. (2020) Nicholson, W. B., Wilms, I., Bien, J. and Matteson, D. S. (2020). High dimensional forecasting via interpretable vector autoregression. Journal of Machine Learning Research, 21 (166), 1–52.
  * Sims (1980) Sims, C. A. (1980). Macroeconomics and reality. Econometrica, 48 (1), 1–48.
  * Tibshirani (1996) Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 58 (1), 267–288.
  * Tibshirani et al. (2016) Tibshirani, R. J., Taylor, J., Lockhart, R. and Tibshirani, R. (2016). Exact post-selection inference for sequential regression procedures. Journal of the American Statistical Association, 111 (514), 600–620.
  * Uematsu et al. (2019) Uematsu, Y., Fan, Y., Chen, K., Lv, J. and Lin, W. (2019). Sofar: Large-scale association network learning. IEEE Transactions on Information Theory, 65 (8), 4924–4939.
  * van de Geer et al. (2014) van de Geer, S., Bühlmann, P., Ritov, Y. and Dezeure, R. (2014). On asymptotically optimal confidence regions and tests for high-dimensional models. The Annals of Statistics, 42 (3), 1166 – 1202.
  * Wang et al. (2021) Wang, D., Zheng, Y., Lian, H. and Li, G. (2021). High-dimensional vector autoregressive time series modeling via tensor decomposition. Journal of the American Statistical Association, 0 (0), 1–19.
  * Wasserman and Roeder (2009) Wasserman, L. and Roeder, K. (2009). High-dimensional variable selection. The Annals of Statistics, 37 (5A), 2178 – 2201.
  * Wong et al. (2020) Wong, K. C., Li, Z. and Tewari, A. (2020). Lasso guarantees for $\beta$-mixing heavy-tailed time series. The Annals of Statistics, 48 (2), 1124 – 1142.
  * Zhang and Zhang (2014) Zhang, C.-H. and Zhang, S. S. (2014). Confidence intervals for low dimensional parameters in high dimensional linear models. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 76 (1), 217–242.