Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
97,315 characters · 21 sections · 33 citation commands
Boosted p-Values for High-Dimensional Vector Autoregression
\doublespacing
JEL Classification: C18, C32
Keywords: Least-squares boosting, linear regression, vector autoregression, p-values
\doublespacing
Since its introduction as a macroeconometric framework in sims1980var, 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 tibshirani1996lasso or imposing other penalties are flourishing. The use of the lasso or its variants can be found in fujita2007,lozano2009grplasso,basu2015lasso,medeiros2016lasso,nicholson2020grplasso,wong2020lasso, among others. The Dantzig selector in han2015directVAR and the tensor method in wang2021tensor also directly address the estimation problem in high-dimensional VARs. Methods such as those in chen2012srrr,uematsu2019sofar 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 friedman2001gbm 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 lutz2006boosting. 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 (ref)). 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 wasserman2009hdselection,meinshausen2009pval,zhangzhang2014JRSSBlowdimension,vandeGeer2014AOSoptimalinterval, 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 (ref).
A related line of research is the recent development in post-selection inference for the lasso, where works in lockhart2014lassopval,lee2016lassopval,tibshiranietal2016lassopval 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, lee2016lassopval,tibshiranietal2016lassopval 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 freundetal2017boosting 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 mccracken2016data 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 \url{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.
We begin by introducing the standard VAR model and some notations. Consider a $d\times1$ vector time series $y_t$ with $t=1,\cdots,T$. A $p$th-order VAR (VAR($p$)) is defined as
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 (ref) in a matrix format. Let $'$ denote matrix/vector transpose. Define
In matrix form, (ref) can be written as
The multivariate regression format is given by
where
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 lutz2006boosting. Alternatively, one can use the “row-boosting" procedure in lutz2006boosting 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.
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). (ref) becomes
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 (ref) so that they are grouped by variables
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 (ref) can be written as
where $\bm{\phi}_g = [\phi_{(1)},\cdots,\phi_{(j)},\cdots, \phi_{(d)}]'$. $\phi_{(j)}'$ is a $p \times d$ coefficient matrix, each row of which comes from the $j$th row of $\phi_1',\cdots,\phi_p'$.
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$
where $ \hat{\mathbf{R}}^{(k-1)} = \hat{\mathbf{R}}^{(0)} = \mathbf{Y} $ if $ k=1 $.
If we define the $p \times d$ matrix
the estimate in (ref) with a learning rate $\nu$, can be updated as
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 \ge 1$, the algorithm becomes
We call this algorithm LS-Boost1. The matrix $\mathbf{B}^{(k)}$ in (ref) is the common hat matrix, and the sum of its diagonal elements gives the degree of freedom up to boosting step $k$.
An approach with a sparser solution is to, instead of using all lags of variable $\mathbf{X}_{(j)}$ in (ref) 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}'$ be the $s$th row of $\phi_{(j)}'$, the regression equation for the “row-boosting" algorithm can be written as
where $\phi_{(j)s}'$ is a $1 \times d$ row vector, which is the $s$th column of $\phi_{(s)}$. Define
The update equation for the coefficient estimate at step $k$ is given by
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 \ge 1$, the algorithm becomes
We call this algorithm LS-Boost2. Note that, although (ref) may look the same as (ref), their values are different since (ref) is based on $\mathbf{X}_{(j_k)s_k}$.
A few remarks are in order.
We discuss the asymptotic result for the two boosting methods introduced in (ref).
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 (ref), rewrite (ref) in a recursive form to have
Following the definitions in (ref), (ref) can be further written as
where we define
(ref) 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 (ref). 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 (ref) to rewrite
It can be shown that terms such as $\tilde{\mathbf{A}}_{j}^{(k)} \mathbf{X}_{(1)} \phi_{(1)}'$ 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.
(ref) are standard in VAR modeling. In (ref), 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 (ref) 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 $\text{vec}$ be the matrix vectorization operator and $\otimes$ be the Kronecker product.
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 (ref)). (ref) gives a familiar form of the asymptotic result for a VAR that can be found in standard textbooks such as hamilton1994time,lutkepohl2005time. 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 (ref), 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} $.
To gain more insights into (ref), 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
Under (ref), all cross product terms can be shown to converge to a constant as $T \rightarrow \infty$, and we have
where $\mathbf{X}_{(j_1),t}'$ is the $t$th column of $\mathbf{X}_{(j_1)}'$, $u_t$ is the $t$th row of $\mathbf{u}$, and
Let the limit of $ \nu \left(\frac{\mathbf{X}_{(j_1)}'\mathbf{X}_{(j_1)}}{T}\right)^{-1}$ be $c$. Vectorizing (ref) gives
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$. (ref) 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 (ref). It can be shown that $E(u_{t_i}|\mathbf{X}_{(j_1),t}') = 0$, because $u_{t_i}$ is the error at time $t$ while $\mathbf{X}_{(j_1),t}'$ includes $p$ lags before time $t$.
where we replace $c$ with $\nu (\frac{\mathbf{X}_{(j_1)}'\mathbf{X}_{(j_1)}}{T})$, $E(\mathbf{X}_{(j_1),t}' \mathbf{X}_{(j_1),t})$ with $\mathbf{X}_{(j_1)}' \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.
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 (ref). 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 (ref) becomes
Since $\mathbf{H}$ is idempotent with eigenvalues $0$ and $1$, Theorem 4.3.1 in hornjohnson1985matrixanalysis 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}'$, 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 buhlmannandyu2003JASAL2boosting 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 (ref) 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 (ref) and will also converge to the same limit. In (ref), we provide a different perspective on the convergence property of the LS-Boost estimator when $k$ is large.
Next we discuss the asymptotic results for the estimator described in (ref), and it will be similar to (ref).
By the definition of $\hat{\beta}_{(j)s}^{(k)\prime}$ in (ref), we have
where
Hence, similar to (ref), we have
The derivation of the asymptotic distribution result largely follows the proof of (ref), and we give the result as a corollary in the following.
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$.
While (ref) and (ref) 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 (ref), 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 (ref) 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 freundetal2017boosting. The following result is non-asymptotic, but it can characterize the behavior of the estimator as $k \rightarrow \infty$.
Consider the VAR model in (ref). Let $\lambda_{\text{pmin}}(\mathbf{X}_g'\mathbf{X}_g)$ be the smallest non-zero eigenvalue of $\mathbf{X}_g'\mathbf{X}_g$ and define the linear convergence rate coefficient
freundetal2017boosting 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 (ref).
See the supplement for the proof. (ref) is similar to Theorem 2.1 in freundetal2017boosting. 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 (ref), $\hat{\bm{\phi}}_g^{(k)}$ in (ref) is defined as $\hat{\bm{\phi}}_g^{(k)} = \left[\hat{\phi}_{(1)}^{(k)}, \cdots, \hat{\phi}_{(d)}^{(k)}\right]'$.
Notice that results in (ref) require very few assumptions to hold. They are derived for the fixed-$T$ case, and moment conditions in (ref) 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}}$. (ref) 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, (ref) 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 lutz2006boosting obtains a result for the consistency of LS-Boost prediction when $T$ and $k$ go to infinity simultaneously, whereas (ref) 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 (ref), we can immediately deduce the behavior of $\hat{\bm{\phi}}_g^{(k)}$ when $k \rightarrow \infty$ with a fixed $T$, and (ref) 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.
We add some additional discussions on the estimator and its p-value.
Incremental hypothesis testing. It is mentioned in (ref) 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 (ref). (ref) implies that when $k=1$, as $T \rightarrow \infty$, $\text{vec}(\hat{\phi}_{(j)}^{(k)\prime})$ converges to $c_{(j_1)}^{(1)}$ in (ref) 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 (ref). 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 (ref), $\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
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 (ref) 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 (ref) is not about the parameter in the population regression equation. Thus, when rejecting the null in (ref) 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
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 (ref) is obtained from the asymptotic distribution in (ref). 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. (ref) 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 hastie2015slsparcity and Figure 1 in tibshiranietal2016lassopval.
Since $\phi_{\text{(j),LS},rq}^{(k)}$ in (ref) 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 (ref) 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 \textit{p}-value. An incorrect application of the \textit{p}-value is to use c.v. first to select a model then apply the \textit{p}-values to further pare down parameters. It is incorrect because the \textit{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 \textit{p}-values so that the \textit{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 (ref)). 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.
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.
Let us consider the following stationary bivariate VAR selected from equation (4.2.1) in lutkepohl2005time.
We simulate a sample of $500$ observations and report in (ref) 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 hurvichandtsai1989bcAIC for the two boosting methods.
In (ref), 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.
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. (ref) plots the average s.e. of $100$ A.I.C.-selected LS-Boost1 models for the $8$ parameters in (ref) 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. (ref) 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 (ref), the bivariate VAR has a unique LS solution and (ref) predicts the LS-Boost1 estimator converges to this unique LS solution.
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.
We use a VAR(2) model with no intercept in our simulation. Consider the model
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 (ref) for more details).
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 (ref) in a VAR(1) format (see equation 10.1.11 in hamilton1994time)
where
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 (ref) 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 (ref).
(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.
We mentioned in (ref) 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 (ref) 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 $\texttt{family = "mgaussian"}$. For the $j$th column in $\mathbf{X}$ in (ref), 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 uematsu2019sofar and the sparse reduced-rank regression (SRRR) in chen2012srrr, 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 (ref), this is a valid way to use the proposed p-value.
(ref) report the simulation results for the three models described in (ref) 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 (ref), 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 (ref) and its performance in (ref) 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 (ref) 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.
In this section, we apply the proposed p-value to the VAR modeling of the monthly macroeconomic data set provided in mccracken2016data. 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.
(ref) 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 (ref) in the following.
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.
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. (ref) give an illustration of the VAR(1) coefficients. Variables appear in (ref) are described in (ref).
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 (ref) 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 (ref) because most of the corresponding coefficients are positive; the sector will be completely filled if we combine the links in (ref).
The links reveals various properties in the data. Consider the GS1 sector in (ref) 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 mccracken2016data. 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 (\texttt{FEDFUNDS}). Using the $758$ observations in the data, we add \texttt{FEDFUNDS} back to \texttt{AAAFFM} and compute its correlation with \texttt{GS1}, and the result is $0.899$, indicating these two time series are positively correlated. Hence, it is important to interpret the result between \texttt{GS1} and \texttt{AAAFFM} in (ref) as the negative relationship between the short-term interest rate and corporate bond yield \textit{spread}. Similar explanations applies to other links between \texttt{GS1} and other bond yields or long-term interest rates in that figure.
(ref) 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 (ref) 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 (ref).
(ref) 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.
(ref) 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.
(ref) 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 (ref). Our empirical example has $104$ equations. We select the equation for civilian unemployment rate (UNRATE) as an example.
(ref) 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 \texttt{GS1} is also small, but climbed above the $5\%$ cutoff after around step $50$, but is below $5\%$ at step $17$. The \textit{p}-values for \texttt{BAAFFM} and \texttt{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.
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.
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.
\spacing{1.45}
\setcounter{page}{1} \spacing{1.42}