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.
89,533 characters · 10 sections · 64 citation commands
Complete Subset Averaging for Quantile Regressions
Quantile regression (QR) has emerged as an essential tool since koenker1978regression (see, e.g.\ koenker_2005). QR estimates the response of conditional quantiles of outcome variables with respect to changes in the covariates. The entire response distribution of outcome variables in economic models provides a broader insight than the classical mean regression. Moreover, in many economic applications, tail quantiles have highly valuable information. See, for example, wage distribution in labor economic applications buchinsky1998dynamics and stock return quantiles (Value-at-Risk) in financial market analysis duffie1997overview. Recently, policymakers have begun to pay attention to the left tail quantiles of GDP growth (Growth-at-Risk) as a measure of downside risks associated with tight financial conditions \citep*{adrian2019vulnerable}. There has also been an increasing interest in climate change, in particular, more frequent and intense extreme weather conditions. A tail quantile is the main object of interest in this analysis \citep*{bhatia2019recent}. Estimation, inference, and prediction of the conditional quantiles are thus important but require a careful econometric analysis due to their nonlinear structure and nonstandard limit theory.
In this paper, we propose a novel prediction method based on complete subset averaging (CSA) for quantile regressions. Following lu2015jackknife, we work on the framework such that all models under consideration are potentially misspecified and that the dimension of regressors goes to infinity as the sample size increases. The CSA method that we propose works as follows. First, pick the numbers of regressors $k$ out of all regressors $K$ available in the data. Then, there exist $K!/(k!(K-k)!)$ complete subsets of size $k$. Second, estimate all the quantile regression models and save all the conditional quantile predictors from each model. Finally, the conditional quantile predictor is constructed as the average of all the quantile predictors estimated in Step 2. Since we average over the complete subsets, the number of models is much larger than the usual model averaging methods selecting the weight of each model. We propose to use an equal weight but select the optimal size of the complete subset $k^*$ based on the leave-one-out cross-validation method.
The CSA approach has a couple of advantages over the existing model averaging method which adopts sophisticated weighting schemes. First, it may produce better forecasts in practice because there is no sampling variance from the weight estimation. This result is already reported both in the forecasting and machine learning literature in the mean regression setup (see, e.g.\ breiman1996bagging, clemen1989combining, stock2004combination, smith2009simple, and \citet*{elliott2013complete}). Second, it does not ask a researcher to choose the initial set of models and the order of each model. In practice, the model averaging methods with different weights usually construct the set of models in an encompassing way and the forecasting performance could depend on the researcher's discretion. Third, CSA averages over a larger number of submodels and one could expect an additional noise reduction from it. However, CSA is possibly more demanding in computation, and we will discuss this issue in detail later.
The contribution of this paper is twofold. First, building upon the theory of lu2015jackknife, we show that the complete subset quantile regression (CSQR) estimator converges the pseudo-true value and satisfies asymptotic normality under mild regularity conditions. The uniform convergence property of CSQR is also provided. Based on these pointwise and uniform limit theories, we prove the asymptotic optimality of $\ensuremath{\widehat}{k}$ in the sense of li1987asymptotic. Second, we implement the CSA method and show that it performs quite well both in simulations and real data sets. Especially, we show that the performance is still satisfactory when we use a fixed number of subsets randomly drawn from the complete subsets when the time budget does not allow estimating the quantile regressions of the whole subsets. We also provide regularity conditions on the choice of the fixed number of subsets. Finally, we provide a theory that compares the performance of equal weighting and optimal weighting in quantile regression. This result justifies our intuition such that optimal weighting forecasts poorly when the number of models increases and extends the existing result in mean regression.
Finally, we summarize related literature. lu2015jackknife and \citet*{elliott2013complete} are closely related to this paper. The former proposes the jackknife model averaging (JMA) method for the quantile prediction problem and derives the nonstandard asymptotic properties of the estimator. Our approach is different from theirs in that we use complete subsets for models to be averaged and that we choose a scalar $\ensuremath{\widehat}{k}$ from the cross-validation method instead of a weighting vector $\ensuremath{\widehat}{w}$. The latter proposes the CSA method in the mean prediction problem and shows by simulation studies that the CSA predictor outperforms alternative methods like bagging, ridge, lasso, and Bayesian model averaging. However, they do not show any optimality result of the estimator. hansen2007least and hansen2012jackknife show the optimality of model averaging based on the Mallows criterion and that of the jackknife model averaging, respectively. ando2014model propose a model averaging method in a high-dimensional setting and show the optimality result. komunjer2013quantile provides a great review on the quantile prediction problem of time-series data. meinshausen2006quantile proposes a quantile prediction method based on random forest. lee2016predictive studies the inference problem of the predictive quantile regression when the regressors are persistent. In the empirical finance literature, \citet*{meligkotsidou2019quantile,meligkotsidou2021out} apply complete subset quantile regression to forecast realized volatility and the risk premium.
The rest of the paper is organized as follows. Section (ref) introduces the model and the CSQR estimator. Section (ref) presents the asymptotic properties of the CSQR estimator and the asymptotic optimality. The Monte Carlo simulation results are reported in Section (ref). Section (ref) investigates two empirical applications and illustrates the advantage of the proposed method. Section (ref) concludes. All the proofs are deferred to the appendix.
We use the following notation. For a matrix $A$, $\left\Vert \cdot \right\Vert $ represents its Frobenius norm $\left\Vert A\right\Vert =\sqrt{tr(AA^{\prime })}$. Let ${\Greekmath 0115} _{\min }\left( A\right) $ and ${\Greekmath 0115} _{\max }\left( A\right) $ denote the smallest and largest eigenvalues of $A$. We use the notation $x_n \approx y_n$ to denote $x_n=y_n+o_p(1)$; and $a_n \ll b_n$ to denote $a_n=o(b_n)$.
In this section, we lay out the model under study and propose the complete subset averaging (CSA) quantile predictor. We also discuss the choice of the subset size based on the cross-validation method.
Consider a random sample $\{(y_{i},x_{i}^{\prime })\}$ for $i=1,\ldots ,n$, where the dimension of $x_{i}$ can be countably infinite. Following lu2015jackknife, we assume that $\{(y_{i},x_{i}^{\prime })\}_{i=1}^n$ is generated from the following linear quantile regression model: for ${\Greekmath 011C} \in (0,1)$,
where ${\Greekmath 0116} _{i}={\Greekmath 0116} _{i}({\Greekmath 011C} ):=\sum_{j=1}^{\infty }{\Greekmath 0112} _{j}x_{ij}$, $ {\Greekmath 0112} _{j}={\Greekmath 0112} _{j}({\Greekmath 011C} )$, ${\Greekmath 0122} _{i}={\Greekmath 0122} _{i}({\Greekmath 011C} ):=y_{i}-Q_{y}({\Greekmath 011C} |x_{i})$, and $Q_y({\Greekmath 011C}|x)$ is the ${\Greekmath 011C}$-th conditional quantile function of $y$ given $x$. Note that we drop ${\Greekmath 011C}$ from each expression for notational simplicity and that ${\Greekmath 0122} _{i}$ satisfies the quantile restriction $P({\Greekmath 0122} _{i}({\Greekmath 011C} )\leq 0|x_{i})={\Greekmath 011C} $. Equivalently, we can also express $Q_{y}\left( {\Greekmath 011C} |x_{i}\right) :=\sum_{j=1}^{\infty }{\Greekmath 0112} _{j}({\Greekmath 011C} )x_{ij}$ as is often done in the quantile regression literature.
We consider a sequence of covariates available, which approximate the above quantile regression model:
where $b_i=b_i({\Greekmath 011C}):={\Greekmath 0116}_i({\Greekmath 011C}) - \sum_{j=1}^{K_n}{\Greekmath 0112} _{j}({\Greekmath 011C} )x_{ij}$ is the approximation error and $K_n$ is the total number of available regressors that may increase as the sample size $n$ increases. Thus, we presume that all models are misspecified in a finite sample as in hansen2007least.\footnote{Using quantile crossings, phillips2015halbert also shows that quantile regresssion is always at the risk of model misspecification unless the parameters are local to constant over ${\Greekmath 011C}$. }
Given $K_n$ regressors, we consider a model composed of $k$ regressors, where $k \in \{1,2,\ldots,K_n\}$. There are $\frac{K_n!}{k!(K_n-k)!}$ different ways to select $k$ regressors out of $K_n$. Therefore, a subset of size $k$ is composed of $M_{(K_n,k)}=\frac{K_n!}{k!(K_n-k)!}$ different elements and a model is defined as a single element of them. We use index $m_{(K_n,k)} \in \{1,2, \ldots ,M_{(K_n,k)}\}$ for each model. For example, consider that we have $K_n=3$ regressors $\{x_{i1}, x_{i2}, x_{i3}\}$ and construct a subset of size $k=2$. Then, we have $M_{(3,2)}=3$ different ways to choose a model as follows: $(x_{i1}, x_{i2}), (x_{i1}, x_{i3}),$ and $(x_{i2}, x_{i3})$. Each model is indexed by $m_{(3,2)}\in\{1,2,3\}$. For succinct notation, we drop all subscripts from $K_n$, $M_{(K_n,k)}$, and $m_{(K_n,k)}$ and denote them as $K$, $M$, and $m$ unless there is any confusion.
We now consider a quantile regression model with regressors in a complete subset. Let model $m$ with a size $k$ be given. For observation $i$, let $x_{i(m,k)}$ be a $k$-dimensional vector of regressors corresponding to model $m$, i.e.\ $x_{i(2,2)}=(x_{i1},x_{i3})$ in the above example. We can construct a linear quantile regression model with regressors $x_{i(m,k)}$:
where $b_{i(m,k)}:={\Greekmath 0116} _{i}-x_{i(m,k)}^{\prime }\Theta _{(m,k)}$ is again the approximation error when we use only $x_{i(m,k)}$ regressors. The model in (ref) is estimated by the standard method in linear quantile regression:
where $\mathbf{\Theta}$ is a parameter space and ${\Greekmath 011A} _{{\Greekmath 011C} }(u):=u({\Greekmath 011C} -1\{u\leq 0\})$ is the check function. Note that the estimator $\ensuremath{\widehat}{\Theta}_{(m,k)}$ is defined for each subset size $k$ and for each model $m$ with $k$ regressors. As noted above, we can think of $M$ different models and corresponding estimators that have $k$ regressors.
We have a few remarks here. First, we use the subscript $(m,k)$ to denote a generic model with $k$ regressors. However, the index set $\{1,\ldots,M_{(K_n,k)}\}$ itself is defined in terms of $k$, which implies that $m$ is also determined by $k$. Recall the original notation $m_{(K_n,k)}$ above. Therefore, model $m \in \{1,\ldots,M_{(K_n,k)}\}$ has the same number of regressors $k$ and we cannot choose $m$ and $k$ in an arbitrary way. Second, we allow that the subset size $k$ goes to infinity as $n$ increases. In other words, there exists a sequence of subset sizes $\{k(n)\}$ that diverges. This setting is natural as the upper bound $K_n$ goes to infinity as $n$ increases. Note that the number of regressors in each model ($k_m$ in their notation) is also allowed to diverge in lu2015jackknife. Both approaches allow more complex models to be averaged as $n$ grows, which is measured by $k$ and $k_m$, respectively. However, lu2015jackknife require controlling the growth rates of $M$ and $\max_m k_m$, separately. The proposed method constructs submodels based on the complete subsets, and $M$ is tightly related to $K$ and $k$. As a result, the regularity condition on the complexity of the models is expressed only in terms of $K_n$ (see Assumption (ref) in Section (ref)).
We finalize this subsection by defining the complete subset averaging (CSA) quantile predictor. Let the size of the complete subset $k$ be given. For each model, we estimate the parameter $\ensuremath{\widehat}{\Theta}_{(m,k)}$ by (ref) and construct the linear index $x_{(m,k)}'\ensuremath{\widehat}{\Theta}_{(m,k)}$. The CSA quantile predictor of $y$ given $x$ is defined as a simple average of those indices over $M$ different models:
The CSA quantile predictor is different from the JMA quantile predictor of lu2015jackknife in two respects. First, we do not select the set of models to be averaged since we average over the complete subsets of size $k$. Second, CSA does not estimate the weights over different models. The idea of averaging over the complete subsets was first introduced by \citet*{elliott2013complete} in the conditional mean prediction setup. Heuristically speaking, since the weights can be seen as additional parameters to be estimated in the model, the equal weight could perform better in a finite sample when the number of models (i.e.\ the dimension of a weight vector) is large.
We propose to choose the subset size $k$ using the leave-one-out cross-validation method. We will show in the next section that the subset size $\ensuremath{\widehat}{k}$ chosen by this method is optimal in the sense that it is asymptotically equivalent to the infeasible optimal choice.
For $k=1,\ldots,K$, we define a cross-validation objective function as follows:
where $\widehat {\Theta}_{i(m,k)}$ is the jackknife estimator for $\widehat { \Theta}_{(m,k)}$, which is estimated by (ref) without using the $ i $-th observation $(x_i,y_i)$, and $\widehat {y}_{i}(k)$ is a corresponding jackknife CSA quantile predictor for the $i$-th outcome variable $y_i$. The prediction error is measured by the check function ${\Greekmath 011A}_{{\Greekmath 011C}}(\cdot)$. Then, we can choose the complete subset size $k$ that minimizes the cross-validation objective function as follows:
After choosing the complete subset size, the CSA quantile predictor is finally defined as
where the plugged-in $\ensuremath{\widehat}{k}$ is chosen by (ref).
We finalize this subsection by adding some remarks on computation. First, we propose to use a fixed number $M_{max}$ of random draws of models when $M$ is too large to implement the method. Since $M=K!/(k!(K-k)!)$, it can be quite large when the model has large potential regressors. The simulation studies in Section (ref) reveal that the CSA quantile predictor still performs well with a feasible size of submodels randomly drawn from the complete subsets. We also provide regularity conditions that assure the asymptotic equivalence between using $M$ and $M_{max}$ in Section (ref). Second, the proposed jackknife method can be immediately extended to the $b$-fold cross-validation method, where $b$ is the partition size of the sample. Algorithm (ref) below summarizes the leave-one-out cross-validation method for choosing $\ensuremath{\widehat}{k}$.
In this section, we investigate the asymptotic properties of the complete subset quantile regression (CSQR) estimator. We first provide the pointwise and uniform convergence results of $\widehat{\Theta}_{(m,k)}$ and $\widehat{\Theta}_{i(m,k)}$, respectively. Then, we show the optimality of CSA in the sense of li1987asymptotic, which implies that $\widehat{k}$ is asymptotically equivalent to the infeasible optimal choice of the subset size.
In addition to the model described in Section (ref), we define some notation for later use. Let $f_{y|x}(\cdot|x)$ be a conditional probability density function for generic random variables $x$ and $y$. Since all models are potentially misspecified in the model averaging literature, we define the pseudo-true parameter value for any given $(m,k)$:
Let ${\Greekmath 0120}_{{\Greekmath 011C}}(c) :={\Greekmath 011C} - 1\{ c \le 0 \}$. For any $(m,k)$ such that $m=1,\ldots,M$ and $k=1,\ldots,K$, we define
and
We need the following regularity conditions.
Conditions (i)--(ii) in Assumption (ref) are the standard i.i.d.\ and the quantile restrictions. Assumption (ref)(iii) requires some finite moment restrictions to achieve the probability bounds of various sample mean objects in the proof. Assumption (ref) allows conditional heteroskedasticity. Note that the eigenvalues of $A_{(m,k)}$ and $B_{(m,k)}$ are bounded and bounded away from zero for a given $(m,k)$. However, these bounds $(\underline{c}_{A(m,k)}, \underline{c}_{B(m,k)}, \overline{c}_{A(m,k)}, \overline{c}_{B(m,k)})$ can converge to zero or diverge to infinity as $n$ increases. The speed of convergence is restricted by Assumption (ref) (iv). These bounded eigenvalue restrictions are commonly imposed in the literature that studies the increasing dimension of parameters (see, e.g.\ portnoy1984asymptotic, portnoy1985asymptotic). Assumptions (ref)--(ref) are standard and similar to those in lu2015jackknife. See the additional remarks therein. Assumption (ref) imposes some regularity conditions on the number of potential regressors $K_n$ and the sequence of the uniform bounds $(\underline{c}_A, \underline{c}_B, \overline{c}_A, \overline{c}_B)$. Different from the regularity condition of JMA in lu2015jackknife, we need not restrict the growth rate of potential models $M$ directly since $M_{(K_n,k)}$ is determined by $K_n$. However, $M_{(K_n,k)}$ increases very quickly at a factorial rate of $K_n$ and we need a stronger restriction on $K_n$. As noted in Assumption (ref)(ii), $K_n$ can increase at most the logarithmic rate of $n$. In the case of JMA, the number of regressors can increase at the polynomial rate if we set $\bar{k}=k_M=M$ in their notation. This is a trade-off in proving the uniform convergence results over a larger index set than that of JMA. We discuss this point in detail below in Theorem (ref). The second part of Assumption (ref)(ii) holds if \b{c}$ _{A}^{3}/\bar{c}_{A}\bar{c}_{B}$ is bounded away from zero or converges to zero at the slower rate than $\log (\log n)/\log n$ when $K$ increases at the rate of $\log n$.
First, we prove the convergence rate and the asymptotic normality of $\ensuremath{\widehat}{\Theta}_{(m,k)}$ when the dimension of parameter $k$ increases.
This theorem provides an asymptotic theory for the quantile regression estimator when the model is misspecified and the number of parameters diverges to infinity as similarly seen in lu2015jackknife. The convergence rate in (i) is a standard result when $k$ diverges as $n$ increases. To show the asymptotic normality with a diverging number of parameters, we also consider an arbitrary linear combination of $\widehat{\Theta}_{(m,k)}$ represented by $C_{(m,k)}$. The difference between two estimators, CSA and JMA, originates from the fact that CSA chooses the total number of the regressors $K_n$ first and the number of complete subset models $M_{(K_n,k)}$ follows automatically for each $k=1,\ldots,K_n$, whereas JSA selects the set of models $M_n$ (in their notation) in advance. Then, the size of regressors $k_m$ in case of JSA is determined by the sequence of models $m=1,\ldots, M_n$ chosen by a researcher. Although there are slight differences in the definition of $c_{A(m,k)}$ and $c_{B(m,k)}$ and their bounds from those in lu2015jackknife, the proof of Theorem (ref) is identical to theirs, so is omitted.
We next turn our attention to the uniform convergence results of $\widehat{\Theta}_{i(m,k)}$ and $\widehat{\Theta}_{(m,k)}$. In addition to its own interest, the uniform convergence rates in the next theorem are required to prove to the asymptotic optimality of $\widehat{k}$.
Since CSA is defined on the index sets of $m$ and $k$, the uniform convergence rates are defined over those sets, $m\in \{1,\ldots, M\}$ and $k\in \{1, \ldots, K\}$. In case of $\widehat{\Theta}_{i(m,k)}$, we need additional uniformity over $i\in\{1,\ldots,n\}$. As a result, the regularity conditions that control the growth rates of $K_n$ and $M_{(K_n,k)}$ are different from those of JMA in Assumption (ref) (ii). As discussed before, since the number of complete subsets increases at the factorial rate of $K_n$, we need a restriction on $K_n$ slightly stronger than that of JMA. We follow the proof strategy in lu2015jackknife which extends the results of rice1984 by using the inequality in shibata1981optimal,shibata1982amendments. To handle the different growth rates, we provide new technical lemmas. The proof of Theorem (ref) as well as these lemmas are provided in the appendix. Finally, the uniform convergence rates are expressed in terms of the sample size $n$ and the total number of regressors $K$ that goes to infinity as $n$ increases.
We next prove the prediction equivalence when we replace $M$ with $M_{max}$. Let $\mathcal{M}_{max}$ be a subset of $\{1,\ldots,M\}$ such that $M_{max}$ elements are randomly drawn. Define $\ensuremath{\widetilde}{y}(k)$ to be the CSA quantile predictor using only $M_{max}$ models:
Let $y^*_k := \lim_{M\rightarrow \infty} M^{-1} \sum_{m=1}^M E\left[x_{(m,k)}' \Theta^*_{(m,k)}\right] < \infty$. We show the validity of $M_{max}$ in the following theorem:
The rate requirement for $M_{max}$ is mild and $M_{max}=O(n^{1/2})$ would work given $K=O(\log n)$. The uniform boundedness assumption on $\Vert x_{(m,k)} \Vert$ is weak and holds easily in most applications. We have some remarks on the uniform convergence assumption of the model average with the pseudo-true parameter $\Theta^*_{(m,k)}$. Let $z_{(m,k)} = x_{(m,k)}'\Theta^*_{(m,k)} - E\left[x_{(m,k)}' \Theta^*_{(m,k)}\right]$. Note that $k$ is discrete and the functional class size over $k$ is small. Thus, it depends on the dependent structure of $z_{(m,k)}$ to hold the uniform law of large numbers. For example, consider the following maximal inequality: for ${\Greekmath 010E}>0$,
where the second line holds from the Markov inequality. Since $K/M =o(1)$, a sufficient condition for the uniform convergence is $\max_{1 \le k \le K}E[\sum_{m=1}^M z_{m,k}]^2/M =O(1)$. If $z_{(m,k)}$ is covariance stationary over $m$ for all $k$, then the sufficient condition becomes the absolute summability condition $\max_{1\le k \le K} \sum_{j=0}^{\infty} \left\vert E[z_{(m,k)} z_{(m+j,k)}] \right\vert < \infty$. See, e.g. fazekas2001general for more general conditions on the partial sums in a different dependent structure.
We now prove the asymptotic optimality of $\widehat{k}$ in the sense of li1987asymptotic. Following lu2015jackknife, we use the final prediction error (FPE, or the out-of-sample quantile prediction error) as a criterion to evaluate the prediction performance:
where $\mathcal{D}_n:=\{(y_i,x_i):i=1,\ldots,n\}$ is a sample. The next theorem shows that $\ensuremath{\widehat}{k}$ is asymptotically equivalent to the infeasible best subset size choice that is defined as a minimizer of $FPE(k)$.
A similar optimality concept has been adopted in the context of the weighted average estimator (e.g.\ hansen2007least, hansen2012jackknife, and lu2015jackknife) and in the context of the IV estimator (e.g.\ donald2001choosing, kuersteiner2010constructing, and lee2018complete). Different from JMA, CSA considers the complete subsets given $(K_n,k)$ and does not require the pre-selection of models to be considered nor the order of models. Thus, the optimality result is also independent of the initial model selection/ordering issue once the total number of regressors is given. The index set $\mathcal{K}$ of CSA is discrete while that of JMA or the jackknife model averaging in hansen2012jackknife is compact. All require the finite moment condition similar to Assumption (A.1) in li1987asymptotic which is assured by Assumption 1 (iii) above. The idea of complete subset averaging has been adopted in the forecasting literature (e.g.\ \citet*{elliott2013complete, elliott2015complete}, \citet*{rapach2010out}). This is the first formal result to show the optimality of the subset size selection.
Finally, we compare the performance of the nonstochastic equal weight with that of the optimal weight. In the mean prediction context, it has been observed that a simple arithmetic mean, i.e. the equal weight, outperforms the estimated optimal weight. This empirical phenomenon is known as the `forecast combination puzzle' and some formal explanations under the mean squared error are provided by smith2009simple, elliott2011averaging, and \citet*{claeskens2016forecast}, to name a few. Heuristically speaking, it happens when the estimation error of the optimal weight is large enough to dominate the efficiency loss caused by the equal weight. We extend this result to the class of smooth expected loss functions. This is crucial in our analysis since the check function ${\Greekmath 011A}_{{\Greekmath 011C}}(\cdot)$ does not give a closed-form solution, which is different from the mean squared error used in the existing literature.
We consider the following simplified framework to focus on the main idea. Let be $\hat{y}_{1}, \ldots, \hat{y}_{M}$ be predictors for $y$ based on $M$ different models. For example, we can think of $\hat{y}_{m} = X'_{(m,k)}\ensuremath{\widehat}{\Theta}_{(m,k)}$ for any given $k$. Let $w$ be an $M$-dimensional weight vector combining the $M$ predictors. We consider only positive weights with $1_M'w=1$, where $1_M$ is an $M$-dimensional unit vector. Let $\hat{y}:=(\hat{y}_1,\ldots, \hat{y}_M)'$ and $e_m:=y - \hat{y}_m$ be the prediction error of $\hat{y}_m$ and $e:=(e_1,\ldots,e_M)'$ be a vector of these prediction errors. We define the prediction error of the combined predictor as $e_c(w):=y-w'\hat{y}=w'(1_M \cdot y - \hat{y})=w'e$. Then, we can define an optimal weight $w^*$ as
where $\Delta^{M-1}$ is the standard $(M-1)$-simplex and $F(w):=E[L(w;e_c)]$ is an expected loss function. For example, the mean squared error in elliott2011averaging can be written in terms of the quadratic loss function: $F(w)=E[e_c^2]=E[w'e e'w]=w'\Sigma w$, where $\Sigma=E[e e']$. The quantile prediction error adopted in this paper can be written in terms of the check function: $F(w) = E[{\Greekmath 011A}_{{\Greekmath 011C}}(e_c)] = E[{\Greekmath 011A}_{{\Greekmath 011C}}(y-w'\hat{y})]$. Let $\bar{w}:=M^{-1}1_M$ be an equal-weight vector and $\hat{w}$ be an estimator for $w^*$ with $\hat{{\Greekmath 0111}} := \hat{w} - w^*$. To illustrate our main point, we further impose that $E[\hat{{\Greekmath 0111}}]=0$ and $\max_m Var(\hat{{\Greekmath 0111}}_m) = \bar{{\Greekmath 011B}}_{{\Greekmath 0111}}^2 > 0$.
We have some remarks. First, it shows that the equal weight $\bar{w}$ may work better than the estimated optimal weight $\hat{w}$ when we average many models, i.e.\ when $M$ is large. Compared to the optimal prediction error $F(w^*)$, the efficiency loss by $\bar{w}$ is bounded by $2^{-1}\bar{{\Greekmath 0115}}_{max}(1+3M^{-1})$, which converges to $2^{-1}\bar{{\Greekmath 0115}}_{max}$ for large enough $M$. On the contrary, the upper bound of the mean efficiency loss by $\hat{w}$ diverges as $M$ increases. We admit that these upper bounds only reflect the worst case scenario. However, it confirms the intuition formally that the equal weight can outperform the estimated optimal weight under the class of smooth expected loss functions. Second, the prediction error of $\bar{w}$ under a quadratic loss function converges to the optimal prediction error as $M$ increases. The same result is also proved in Proposition 1 in elliott2011averaging. Different from his result, it does not require decomposing the prediction error into the common component and the idiosyncratic component. This result is summarized in Corollary (ref) below. Third, to achieve the optimality, the estimation errors of the weight $\hat{w}$, $\{\hat{{\Greekmath 0111}}_m\}$, should vanish fast enough. Let $\bar{{\Greekmath 011B}}_{{\Greekmath 0111}}^2=O(c_n)$. A sufficient condition for the optimality is $c_n=o(M_n^{-1})$. For example, if $c_n$ is a parametric rate, $n^{-1/2}$, then $M_n$ should be bounded by $o(n^{1/2})$. When $M_n=O(\log n)$, for example, this condition is satisfied. However, if $M_n$ increases too fast, then $M_n\bar{{\Greekmath 011B}}_{{\Greekmath 0111}}^2$ will diverge and $\hat{w}$ may work worse than $\bar{w}$.\footnote{We thank an anonymous referee and Co-editor for pointing out this intuition. Also, note that it is one sufficient condition. It is still possible that there exists a different set of conditions that guarantee the optimality.} Fourth, $\left\Vert \triangledown_2 F(w) \right\Vert = (\sum_{m=1}^M {\Greekmath 0115}_m^2 )^{1/2}$, where $\{{\Greekmath 0115}_m\}$ are eigenvalues of $\triangledown_2 F(w)$ since $\triangledown_2 F(w)$ is symmetric. Thus, the uniform bound $C$ exists if $\{{\Greekmath 0115}_m\}$ is absolutely summable, $\sum_{m=1}^{\infty} \left\vert {\Greekmath 0115}_m \right\vert < \infty$. Fifth, if we restrict our attention to the expected check function adopted in this paper, $F(w)$ is twice differentiable if the conditional density $f(y|\hat{y})$ is smooth for all $\hat{y}$. From Theorem 1 in \citet*{angrist2006quantile}, we have
where $Q_{{\Greekmath 011C}}(y|\hat{y})$ is the conditional quantile function of $y$ given $\hat{y}$ and $\bar{{\Greekmath 0121}}_{{\Greekmath 011C}}(\hat{y},w):=\int_0^1 (1-u)\cdot f(u\cdot w'\hat{y} + (1-u) \cdot Q_{{\Greekmath 011C}}(y|\hat{y})|\hat{y})du$. Thus, the smoothness of $F(w)$ is implied by the twice differentiability of $f(y|\hat{y})$. Finally, equation (ref) shows that CSA would not work well if we include many irrelevant models. Similar to the quantile regression specification error in angrist2006quantile, we call $\left(w'\hat{y} - Q_{{\Greekmath 011C}}(y|\hat{y})\right)$ the quantile prediction specification error. If there are many irrelevant models, the optimal weight ${\Greekmath 0121}^*$ would be sparse, i.e.\ many elements of ${\Greekmath 0121}^*$ would be zeros. In such a case, CSA with $\bar{{\Greekmath 0121}}=M^{-1}1_M$ results in a larger quantile prediction specification error given $M$ and $n$. For example, if there is only one relevant regressor and all other coefficients ${\Greekmath 0112}_j({\Greekmath 011C})$ equal zero besides one, the complete subsets will be composed of many irrelevant models. As we will see in the simulations studies in the next section, CSA does not perform well under this situation. Therefore, a pre-screening process is desirable to achieve a satisfactory result of CSA.
In this section, we investigate the finite sample performance of the proposed estimator in simple Monte Carlo experiments. We consider two categories of the simulation designs: (i) all candidate models are misspecified, and (ii) candidate models include the true model.
First, we adopt the following data generating process (DGP):
where $x_{i1}=1$ and $(x_{i2},\ldots, x_{i1000})$ follows a multivariate normal distribution, $N(0,\Sigma)$ with $\Sigma_{jk}={\Greekmath 011A}_x$ if $j\neq k$ and 1 if $j=k$. Therefore, the regressors are possibly dependent on each other, which is a more general feature of the design than the existing literature, see, e.g., hansen2007least and lu2015jackknife. The term ${\Greekmath 0122} _i$ follows $N(0,1)$ independent of $x_{ij}$. The sample is i.i.d. over $i$. The population $R^2:=(Var(y_i)-Var({\Greekmath 0122} _i))/Var(y_i)$ is controlled by ${\Greekmath 0112}$. We consider two sample sizes, $n=50, 150$. The number of potential regressors is set to $K=4\log(n)$, which is 15 and 20, respectively. Note that all candidate models are misspecified since there remain many missing regressors in the sample. We consider various DGPs by combining different $R^2 = \{0.1, \ldots, 0.9\}$, ${\Greekmath 011C} = \{0.1, \ldots, 0.9\}$, and ${\Greekmath 011A}_x = \{0.0, 0.1, 0.2, \ldots, 0.9\}$. We consider 38 different DGPs in total and estimate 74 different quantile models.
We compare the performance of the proposed Complete Subset Averaging estimator (CSA) with the Jackknife Model Averaging estimator (JMA) in lu2015jackknife, the $\ell_1$-penalized quantile regression (L1QR) in \citet*{belloni2011L1}, the bootstrap aggregating methods (BAG) in breiman1996bagging and $\ell_2$-penalized quantile regression. L1QR and L2QR are also called the lasso and the ridge regression in the mean regression setup. The set of models used for JMA is constructed in an encompassing way, e.g.\ $\{x_{i1}\}, \{x_{i1},x_{i2}\},\ldots, \{x_{i1},\ldots, x_{i20}\}$. For CSA, we set the maximum submodels to $M_{max}=100$. Thus, we draw 100 models randomly from the complete subsets of size $k$ if $M=K!/(k!(K-k)!)$ is bigger than 100. Furthermore, we reduce some computational burden by applying 10-fold cross-validation when $n=150$. The tuning parameter of L1QR is chosen by Equation (2.7) in belloni2011L1. The bootstrap size of BAG is set to be 1000. The tuning parameter of L2QR is chosen by 10-fold cross-validation over the set $\{0.01, 0.05, 0.1, 0.5, 1.0\}$ which is constructed after some pre-simulation studies.
To compare the performance, we first compute $\mbox{FPE}(r)$ for each replication $r=1,\ldots,R$ as follows. After estimating the model with $n$ in-sample observations, we generate additional 100 out-of-sample observations. Then, $\mbox{FPE}(r)$ is calculated by
where $\hat{y}_x$ is a predicted value by each method. Then, we construct the following three comparison measures:
where each subscript denotes generic notation for a forecasting method. Note that the loss to CSA ratio provides more direct binary comparison of each method to CSA. We set the total number of replications $R=1000$.
Figures (ref)--(ref) and Tables (ref)--(ref) summarize the simulation results over all designs. Overall, the performance of CSA compared to the alternative is quite satisfactory. We first direct our attention to Figure (ref) and Tables (ref)--(ref). In these simulation designs, we vary $R^2$ over $\{0.1,0.2,\ldots, 0.9\}$ while setting ${\Greekmath 011A}_x=0.9$. We consider two quantiles, ${\Greekmath 011C}=0.1$ and $0.5$, respectively. From the four graphs in Figure (ref), we confirm that CSA outperforms the alternative uniformly over $R^2$'s in terms of FPE in both quantiles. The prediction performance of CSA is better when the sample size is small, $n=50$, and the gap decreases as the sample size increases to $n=150$. At ${\Greekmath 011C}=0.5$, L1QR performs the second when $n=50$ but L2QR does the second when $n=150$. Thus, the performance order next to CSA is not stable. AT ${\Greekmath 011C}=0.1$, BAG performs the second overall but it is deteriorated when $R^2$ is very high, e.g.\ $R^2=0.9$. We also note that the performance of CSA is relatively stable over $R^2$ while that of the alternative increases steeply for larger $R^2$ when $n=50$. The same results are confirmed in Tables (ref)--(ref). CSA shows the highest winning ratios over all designs except ${\Greekmath 011C}=0.5$ and $R^2=0.1$, where that of L2QR is slightly higher. When we conduct the binary comparison (loss to CSA), all methods lose more than 50% to CSA over all designs and more than 80% in some designs. Therefore, we conclude that both the winning ratio and the loss to CSA are more favorable to CSA in this set of simulation designs.
In the next simulation, we study the performance over a wider range of quantiles. We vary the quantile ${\Greekmath 011C}=\{0.1,0.2,\ldots,0.9\}$ while setting $R^2=0.5$ and ${\Greekmath 011A}_x=0.9$. The results are summarized in Figure (ref) and Table (ref). In Figure (ref), CSA outperforms the alternative uniformly over all quantiles in both sample sizes followed by BAG and L2QR. Again, the gap decreases as the sample size increases. It is also interesting that all estimators predict better at the tail distributions and they show the largest prediction errors at the median. The winning ratio and the loss to CSA in Table (ref) are also satisfactory.
Third, we check the performance over different levels of dependency among the predictors. We vary ${\Greekmath 011A}_x=\{0,0.1,0.2,\ldots,0.9\}$ while setting $R^2=0.5$ and ${\Greekmath 011C}=0.5$. Since $(x_{i2},\ldots,x_{i1000})$ are generated from the multivariate normal distribution, they are independent when ${\Greekmath 011A}_x=0$. Figure (ref) reveals an interesting point. CSA performs better than the alternative when there exists any correlation between the predictors, i.e.\ ${\Greekmath 011A}_x > 0$. Recall that most simulation studies in the literature consider independent predictors. As we can see from the empirical applications in the next section, however, the predictors are usually correlated with each other. Therefore, it is promising that CSA performs better when there is any correlation among predictors. \citet*{elliott2013complete} also report in the conditional mean prediction settings that the CSA approach performs better when predictors are correlated with each other. In Table (ref), both the winning ratio and the loss to CSA statistics improve dramatically when ${\Greekmath 011A}_x$ is away from zero, where JMA performs the best.
We next consider the second category of simulation designs, where the candidate models include the true DGP. The new simulations are based on the following model:
where we observe all $K$ predictors in the sample. We consider $K=5,15$ when $n=50$ and $K=10,20$ when $n=150$. Similar to the previous simulations, the population $R^2$ is controlled by ${\Greekmath 0112}$. We set $R^2=0.5$, ${\Greekmath 011C}=0.5$, and ${\Greekmath 011A}_x=0.9$. Instead of varying $R^2$, ${\Greekmath 011C}$, and ${\Greekmath 011A}_X$, we consider three signal structures in this simulation:
Therefore, we consider 12 new DGPs in total.
Tables (ref)-(ref) summarize the simulation results. First of all, we take a look at the loss to CSA ratio in the second column (JMA) in these tables. Note that the loss ratio increases as $K$ increases over all different designs, which is expected by the theoretical results in Theorem (ref). Second, CSA performs worse in the sparse signal models compared to the other two designs. As discussed under equation (ref), this is expected from the theory in Section (ref) since the sparse design generates many subsets with totally irrelevant predictors. Third, it is interesting that JMA does not particularly outperform in this setup, where the candidate models include the true one. Also, note that L1QR does not particularly outperform in the sparse signal model. In fact, L2QR performs well over all three signal designs. Given that L2QR is understudied in the literature, this would be an interesting topic for future research.
In sum, we confirm that CSA shows satisfactory finite sample properties via Monte Carlo simulation studies. Related to the forecast combination puzzle, we observe a similar phenomenon in quantile forecasting and confirm some theoretical predictions developed in Section (ref).
In this section, we investigate the performance of the proposed method with real data sets. Specifically, we revisit two empirical applications in lu2015jackknife: (i) quantile forecast of excess stock returns; and (ii) quantile forecast of wages. Following the simulation studies in Section (ref), we compare the performance of the complete subset averaging (CSA) method to the Jackknife Model Averaging (JMA), the $\ell_1$-penalized quantile regression (L1QR), the bootstrap aggregating method (BAG), and the $\ell_2$-penalized quantile regression (L2QR).
The same data set is composed of monthly observations of the US stock market from January 1950 to December 2005 ($T=672$). The dependent variable is the excess stock return. We use the following twelve regressors: default yield spread, treasury bill rate, net equity expansion, term spread, dividend price ratio, earnings price ratio, long term yield, book-to-market ratio, inflation, return on equity, lagged dependent variable, and smoothed earnings price ratio. See lu2015jackknife and campbell2007predicting for the details of the data set. Note that JMA needs to select the order of important regressors, but we do not need such a selection for CSA, BAG, L2QR. L1QR would select important regressors automatically by the $\ell_1$-penalty.
We forecast the one-period-ahead excess stock returns at 0.5 and 0.05 quantiles using various fixed in-sample sizes, $T_1 = 48, 60, 72, 96, 120, 144,\mbox{ and }180$. The forecast performance is measured by the out-of-sample $R^2$ defined as
where $\widehat {y}_{t+1|t}$ the one-period-ahead ${\Greekmath 011C}$-quantile prediction at time $t$ using the data from the past $T_1$ periods, and $\bar{y}_{t+1|t}$ is the unconditional ${\Greekmath 011C}$-quantile for the same $T_1$ periods. The out-of-sample $R^2$ measures the relative performance of a forecast method compared to the unconditional historical quantile. The higher values of $R^2$ imply better forecasting performance.
Table (ref) summarizes the forecasting results. In addition to $R^2$, we report the ranking of each forecasting method, the mean of $\widehat {k}$, and the median of $\widehat {k}$. The upper panel of Table (ref) reports the results when ${\Greekmath 011C}=0.05$. The $R^2$ of CSA is better than that of JMA uniformly over different sample sizes ($T_1$). The gap between two $R^2$'s is substantial except $T=180$. BAG performs well when $T_1$ is small. The performance of L2QR is not satisfactory over all in-sample sizes. We next turn our attention to the lower panel when ${\Greekmath 011C}=0.5$. Again, CSA performs the best or second best except when $T_1=144$. CSA performs better when $T_1$ is small while JMA does better when $T_1$ is larger. Overall, the gap between $R^2$'s is small when ${\Greekmath 011C}=0.5$. As we have observed from the simulation studies, the performance of the two estimators becomes similar as the sample size increases in both panels. It is also noticeable that the selected $\widehat{k}$ of CSA increases as the sample size increases and that CSA selects relatively large $\widehat{k}$ across all $T_1$ and ${\Greekmath 011C}$. Different from ${\Greekmath 011C}=0.05$, BAG performs poorly when ${\Greekmath 011C}=0.5$. L2QR also shows poor performance.
In sum, the performance of CSA is satisfactory in this forecasting exercise. It is quite stable over different in-sample sizes ($T_1$) and different quantiles in terms of the performance ranking. Among the alternative, BAG and JMA perform well in certain quantiles (0.05 and 05, respectively), but they do poorly when we apply them in different quantiles.
In this subsection we conduct the quantile forecast exercises using the Current Population Survey (CPS) data in 1975. The same data set is also used by lu2015jackknife and hansen2012jackknife for quantile and mean forecast exercises, respectively. The sample size is $n=526$ and we use the logarithm of the average hourly wage as the dependent variable. We use the following ten regressors: professional occupation, years of education, years with current employer, female, service occupation, married, trade, SMSA, services, and clerk occupation.
We split the sample into the estimation sample randomly drawn $n_1$ observations and the evaluation sample of $n-n_1$ observations. The estimation sample size varies $n_1=50,100, 150,$ and $200$ and the random splitting is repeated 200 times for each $n_1$. The out-of-sample $R^2$ is defined as
where $\widehat {y}_s$ is the ${\Greekmath 011C}$-th conditional quantile predictor and $ \bar{y}_s$ is the unconditional ${\Greekmath 011C}$-quantile estimate from the estimation sample. Again, $R^2$ measures the prediction performance relative to the unconditional quantile estimate.
Table (ref) summarizes the exercise results.\footnote{$R^2$s of JMA are different from the numbers reported in Table 5 in lu2015jackknife because they implemented the level of wage as a dependent variable which is supposed to be $\log(wage)$. We use $\log(wage)$ in this empirical illustration.} We confirm that CSA shows good and stable quantile prediction performance. In this application, BAG shows quite a similar performance to CSA. Similar to the stock return application, CSA performs better than BAG when ${\Greekmath 011C}=0.5$ and BAG does when ${\Greekmath 011C}=0.05$. The prediction results of JMA, L1QR, and L2QR are worse than CSA and BAG. The performance gaps are larger when the sample size ($n_1$) is small and they narrow as $n_1$ increases. As predicted by the theory and also confirmed in the stock return application, the selected $\widehat{k}$ increases as $n_1$ increases.
In this paper, we propose a novel conditional quantile prediction method based on complete subset averaging of quantile regressions. We show the asymptotic properties of the estimator when the dimension of regressors diverges to infinity as the sample size increases. The size of the complete subset is chosen by the leave-one-out cross-validation method. We prove that the subset size chosen by this method is optimal in the sense that it is asymptotically equivalent to the infeasible optimal size minimizing the final prediction error. The prediction performance in the simulation studies and empirical applications is satisfactory.
We conclude with two potential extensions of the proposed method. First, we can think of a different approach in choosing the complete subset size. Recently, hirano2019analyzing propose a Laplace cross-validation method, where the tuning parameter of interest is chosen by the pseudo-Bayesian posterior mean, and show that it works better than the standard cross-validation method when the risk function is asymmetric. It would be interesting to check how it performs in the CSA quantile prediction. Second, it will be useful if one can extend the results into the time-series data possibly including persistent regressors (e.g.\ fan2019predictive). We leave them for future research.