The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
75,569 characters
\thispagestyle{empty}
\begin{center}
{\huge Time-varying Forecast Combination for High-Dimensional Data\footnote[1]{We thank Xu Cheng, Francis X. Diebold, Atsushi Inoue, Frank Schorfheide, Nese Yildiz and seminar participants at the
University of Pennsylvania, University of Rochester, the 2019 Southern Economic Society meeting and the 12th Econometric Society World Congress for their useful comments and discussions. Any remaining errors are
solely ours. }}\bigskip
\bigskip
$
\begin{array}{c}
\text{{\large Bin Chen}} \\
\text{University of Rochester} \\
\text{ {\large Kenwin Maung}} \\
\text{University of Rochester}
\end{array}
$\bigskip
\bigskip
\bigskip
\end{center}
\textit{Abstract: } In this paper, we propose a new nonparametric estimator of time-varying forecast combination weights. When the number of individual forecasts is small, we study the asymptotic properties of the local linear estimator. When the number of candidate forecasts exceeds or diverges with the sample size, we consider penalized local linear estimation with the group SCAD penalty. We show that the estimator exhibits the oracle property and correctly selects relevant forecasts with probability approaching one. Simulations indicate that the proposed estimators outperform existing combination schemes when structural changes exist. Two empirical studies on inflation forecasting and equity premium prediction highlight the merits of our approach relative to other popular methods. \bigskip \bigskip
\noindent \textit{JEL Classifications}: \textit{C12, C14, C22}
\bigskip
\noindent \textit{Key words: }Cross validation, Forecast combination, High
dimension, Local linear estimation, SCAD, Sparsity.\bigskip
\verb||
\verb||
\pagebreak \setcounter{page}{1}
\section{Introduction}
Multiple forecasts of the same variable are often available to decision makers. As pointed out in the seminal paper by \cite{bates1969combination}, combinations of individual forecasts can outperform individual forecasts as economic systems are highly complex and even the most sophisticated model is likely to be misspecified. It is unlikely that a single model will dominate uniformly. Even if such a best model exists, it is very difficult to identify it in practice since many forecasts might have similar predictive accuracy.
Structural breaks in predictive relationships pose additional challenges when generating out-of-sample forecasts. Individual forecasts may vary with structural changes caused by changes in preferences, institutional evolution or technological progress, among other reasons. It is likely that the combination of forecasts from models with different degrees of adaptability would average out individual effects and outperform forecasts from one specific model. Forecast combination can thus be viewed as a strategy against potential structural changes, in the spirit of portfolio hedging, by offering diversification gains. \citep[see e.g.][]{aiolfi2006persistence, timmermann2006forecast, elliott2005optimal}. Understandably, this ability to deal with both model uncertainty and structural changes in forecasting has motivated many authors to apply forecast combination in various fields, ranging from macroeconomics \citep{elliott2005optimal, stock2004combination} to empirical asset pricing \citep{lin2018forecasting, rapach2010out}, with many reporting significant performance gains over prevailing methods.
Given that the relative performance of different forecasts is likely to change over time, it is natural to consider forecast combination with time-varying weights. Time-varying forecast combination was first proposed by \cite{bates1969combination}, who developed several adaptive estimation schemes for time-varying weights based on exponential discounting or rolling estimation. In the regression context, \cite{diebold1987structural} generalize these schemes to select combination weights that minimize the weighted average of forecast errors. \cite{deutsch1994combination} and \cite{elliott2005optimal} consider a parametric interpretation of the time-varying combination weights by allowing them to be driven by smooth transitions or switching. In an empirical study, \cite{lin2018forecasting} consider the iterated mean combination and the iterated weighted combination, which improve on existing forecast combination schemes by combining them with the historical sample mean forecast. These methods rely on either rolling estimation with a fixed window size or impose certain parametric functional form assumption on the combination weights, which may be restrictive. Therefore, it is desirable to develop an alternative time-varying combination scheme which can hedge against structural changes of an unknown form.
Recently, nonparametric time-varying parameter models have proven to be a reliable tool in identifying the smoothly-varying coefficient functions. Furthermore, it has been shown to adequately capture the evolutionary behavior of economic relationships through various empirical applications. Such models were first introduced by \cite{robinson1989nonparametric, robinson1991time} and further studied by \cite{cai2007trending}, \cite{chen2012testing}, \cite{kristensen2012non}, \cite{zhang2012inference}, \cite{dahlhaus2019towards}, \cite{hongsunwang} among many others. One advantage of the nonparametric
time-varying parameter model is that little restriction is imposed on the functional forms of coefficients, apart from the regularity condition that they evolve smoothly over time. Motivated by this flexibility, we will adopt this framework in estimating the time-varying combination weights.
This paper develops two time-varying forecast combination schemes. When the number of forecasts is small, we consider a new nonparametric estimator for the combination weights and study its asymptotic properties. Our framework is general enough to accommodate some degree of non-stationarity, in particular that of local stationarity, and thus we can avoid taking a hard stance on the time series behavior of forecasts. As such, our estimator can be viewed as a generalization of the classical \cite{granger1984improved} regression estimator. To implement our nonparametric combination scheme, we consider a cross-validation (CV) bandwidth selection method and show that the selected bandwidth converges to the theoretical optimal bandwidth, which minimizes the integrated mean squared combined forecast errors (IMSCFE). Bandwidth selection here is analogous to the optimal selection of window size for classical rolling regression \citep[see for e.g.][]{hongsunwang}.
When the number of potential forecasts is allowed to be at the same order as or even larger than the sample size, we consider a two-stage penalized local linear procedure similar to that of \cite{li2015model} in studying varying coefficient models. Model selection has been an increasingly important topic in econometrics and statistics in the past twenty years, and various penalized likelihood or least-square methods have been studied to handle model selection for high-dimensional data. Examples of
commonly-used penalization schemes include the Lasso
\citep{tibshirani1996regression}, smoothly clipped absolute deviation (SCAD) \citep{fan2001variable}, group Lasso \citep{yuan2006model}, adaptive Lasso \citep{zou2006adaptive}, and the minimax concave penalty (MCP) \citep{zhang2010nearly}. In a related context, model selection for functional coefficient models under the $i.i.d.$ assumption are considered in \cite{wang2009shrinkage}, \cite{wei2011variable} and \cite{li2015model}, and model selection for
high-dimsional linear time series models are studied in \cite{kock2015oracle}, \cite{han2020high} and \cite{diebold2019machine}. We contribute to this growing literature of high-dimensional model selection in time series econometrics by investigating the asymptotic properties of our two-stage estimator and showing that it possesses the oracle property.
Our proposed approach has a number of appealing features. First, the forecast combination weights are modeled as some nonparametric function of time. Unlike \cite{deutsch1994combination} and \cite{elliott2005optimal}, we do not impose parametric assumptions on the functional form of time variation. Second, the finite-dimensional case and the high-dimensional cases are studied in a unified framework. The two-step scheme shrinks the weights of irrelevant forecasts to 0 and provides a practical tool to reduce dimensionality. Third, unlike the majority of the literature on forecast combinations, we investigate the asymptotic properties of our estimator, and establish results on both estimation and selection consistency.
The rest of the paper is organized as follows. In Section 2, we introduce the framework of our nonparametric time-varying combination scheme and develop the estimator when the number of forecasts is small. Section 3 derives the asymptotic properties of the estimator and the CV-selected bandwidth. Section 4 introduces the two-stage penalized estimator for combination weights in high-dimensional forecast combination. Section 5 discusses its implementation and Section 6 studies the oracle property of the estimator. In Section 7, a simulation study is conducted to assess the reliability of the low- and high-dimensional estimators in finite samples. Two empirical examples on inflation forecasting and equity premium prediction are used to illustrate the merits of our approaches in Section 8. Main mathematical proofs are collected in the appendix, while the proofs of some technical results and additional simulations are contained in an online appendix \citep{onlineappendix}.
\section{Nonparametric Forecast Combination}
\label{nonpar}
Assume that a decision maker is interested in predicting some univariate
series $y_{t+1}$, conditional on $I_{t},$ the information available at time $
t$, which consists a set of individual forecasts $f_{t}=
\left(f_{t1,}f_{t2},...,f_{td}\right) ^{\top}$ in addition to current and
past values of $y$, i.e. $I_{t}=\left( y_{s},f_{s}\right)_{s=1}^{t}.$ The
vast majority of studies in the forecasting literature considers a linear
forecast combination model:
\begin{equation*}
y_{t+1}=\omega _{0}+f_{t}^{\top }\omega _{1}+\varepsilon _{t+1},\qquad
t=1,...,T,
\end{equation*}
where $\omega _{0}$ is an intercept and $\omega _{1}$ is a $d\times 1$ vector,
whose component $\omega _{1i}$ can be viewed as the weight assigned to the $
i^{th}$ forecast, where $i=1,...,d.$ Note that the sum of $\omega _{1i}$
need not be unity.\footnote{
Alternatively, a transformed model, which imposes the constraint that $
\sum_{i=1}^{d}\omega _{1i}=1$ and $\omega _{0}=0$ can be considered. Namely,
$y_{t+1}-f_{td}=\sum_{i=1}^{d-1}\omega _{1i}\left( f_{ti}-f_{td}\right)
+\varepsilon _{t+1}.$} Individual forecasting models can be viewed as
local approximations to the true data generating process (DGP) and their forecast
ability is likely to change over time due to the prevalence of structural
changes. Therefore, we consider forecast combination with time-varying
weights:
\begin{equation*}
y_{t+1}=\omega _{0t}+f_{t}^{\top }\omega _{1t}+\varepsilon _{t+1},\qquad
t=1,...,T,
\end{equation*}
where $\left( \omega _{0t},\omega _{1t}^{\top }\right) ^{\top }$ are adapted
to the current information set $I_{t}$.
We opt to estimate the model without imposing any parametric functional form on the time variation of combination weights\footnote{There are at least three ways to estimate the time-varying weights \citep{elliott2005optimal, timmermann2006forecast}. The first method is based on a rolling window estimation with some fixed window length $c$, where $c$ is often selected arbitrarily in empirical studies. The second method assumes the form of a time-varying parameter model, where the combination weights
are assumed to follow a multivariate unit root process. The third method assumes that
weights are driven by switching \citep{elliott2005optimal} or smooth transitions \citep{deutsch1994combination} with some observed or latent state
variable. }. Specifically, we adopt the following framework of the nonparametric time-varying parameter model:
\begin{equation}
y_{t+1}=\omega _{0}\left( t/T\right) +f_{t}^{\top }\omega _{1}\left(
t/T\right) +\varepsilon_{t+1},\qquad t=1,...,T, \label{yt1}
\end{equation}
where $\omega_{0}:[0,1]\rightarrow \mathbb{R}^{1}$ and $\omega_{1}:[0,1] \rightarrow \mathbb{R}^{d}$ are smooth functions of the standardized time $t/T$ over [0,1]. This model was introduced by \cite{robinson1989nonparametric, robinson1991time} and has been studied extensively. The specification that $\omega _{jt} \equiv \omega_j(t/T)$, for $j=0,1$, are functions of the ratio $t/T$ rather than time $t$ itself is a common scaling scheme in the literature which guarantees that the amount of local information increases suitably with the sample size. One difficulty regarding estimation of a predictive model is that no symmetric data is available when making a forecast at any time point $t$. In other words, we are unable to use data from $t+1$ onwards when producing forecasts at time $t$. In the context of nonparametric regression, this is essentially the boundary problem. Note that although local linear smoothing can enhance the convergence rate of the asymptotic bias in the boundary, the asymptotic variance at a boundary point is inevitably larger because we have fewer observations contributing to the estimator on a smaller data interval. To further reduce the variance, we adopt the reflection method following Hall and Wehrly (1991) and Chen and Hong (2012). Specifically, we reflect the data at each data point $t$ and obtain pseudodata $(y_{s+1},f_{s}^{\top})=(y_{2t-s+1},f_{2t-s}^{\top })$ for $t+1\leq s\leq t+\lfloor Th\rfloor,$ where $\lfloor Th\rfloor $ denotes the integer part of $Th$ and $h$ is the bandwidth used in estimation. We use the synthesized data (the union of the original data and pseudodata) to estimate $\beta _{t}=( \omega _{0t},\omega _{1t}^{^{\top }}) ^{\top }$ for each $t$ via local linear estimation.
Let $z_{st}=\left( 1,\frac{s-t}{T}\right) ^{\top }$ and $k_{st}=h^{-1}k \left( \frac{s-t}{Th}\right) ,$ where the kernel $k :[ -1,1] \rightarrow \mathbb{R}^{+}$ is a prespecified symmetric probability density. Examples of $k(\cdot )$ include the uniform, Epanechnikov and quartic kernels. As discussed in \cite{hongsunwang}, $\lfloor Th\rfloor $ is analogous to the window length of a rolling window regression. The local linear parameter estimator for combination weights at time $t$ is obtained by minimizing the local sum of squared residuals:
\begin{equation}
\underset{\gamma \in \mathbb{R}^{2(d+1)}}{\min }T^{-1}\sum_{s=t-\lfloor
Th\rfloor ,s\neq t}^{t+\lfloor Th\rfloor }k_{st}\left[ y_{s+1}-\alpha
_{0}^{\top }x_{s}-\alpha _{1}^{\top }\left( \frac{s-t}{T}\right) x_{s}\right]
^{2}=T^{-1}\sum_{s=t-\lfloor Th\rfloor , s\neq t }^{t+\lfloor Th\rfloor
}k_{st}(y_{s+1}-\gamma ^{\top }q_{st})^{2}, \label{LLeast}
\end{equation}
where $\gamma =(\alpha _{0}^{\top },\alpha _{1}^{\top })^{\top }$ is a $
2(d+1)\times 1$ vector, $\alpha _{j}$ is a $(d+1)\times 1$ coefficient
vector for $(\frac{s-t}{T})^{j}x_{s},$ $j=0,1,$ $q_{st}=z_{st}\otimes x_{s}$
is a $2(d+1)\times 1$ vector, and $\otimes $ is the Kronecker product.
Minimizing \eqref{LLeast} with respect to $\gamma_t$ yields the local linear
estimate of $\beta(t/T)$,
\begin{equation}
\hat{\beta}_t = \hat{\beta}(t/T) = (e_1^\top \otimes I_{(d+1)}) \hat{\gamma}_t, \label{betasol}
\end{equation}
where $e_1 = (1,0)^\top$, $I_{(d+1)}$ is a $(d+1) \times (d+1)$ identity
matrix, and
\begin{equation}
\hat{\gamma}_t = \bigg(\sum_{ s=t-\floor*{Th} , s \neq t}^{t+
\floor*{Th}} k_{st} q_{st} q_{st}^\top \bigg)^{-1} \sum_{ s=t-
\floor*{Th} , s \neq t}^{t+\floor*{Th}} k_{st} q_{st} y_{s+1}.
\label{alter}
\end{equation}
The estimator takes the form of a leave-one-out local
linear estimator considered in \citet{chen2012testing} because of the
predictive design of the regression, which refers to the fact that we do not observe $y_{t+1}$ at time $t$ and hence cannot use it in the estimation.
\section{Asymptotic Properties}
\label{sect.asym}
In this section, we consider the asymptotic properties of the estimator when the number of forecasts is relatively small ($d \ll T$) and thus no regularization is required. To begin, we impose the following regularity conditions.
\bigskip
\noindent \textbf{Assumption A.1 (Mixing condition):} The process $\{R_t\}= (\varepsilon_{t+1}, X_t^\top)^\top$ is a $\beta$-mixing process with mixing coefficient ${\beta^*(j)}$ satisfying $\sum_{l=1}^\infty l^4 \beta^*(l)^{\delta/(1+\delta)} <
\infty$ for some $\delta > 0 $.
\bigskip
\noindent \textbf{Assumption A.2 (Moment conditions):} The following moment
conditions are satisfied: \setlist{nolistsep}
\begin{enumerate}[label=(\roman*), noitemsep]
\item $\sup_{1 \leq t \leq T} E \| R_t \|^{4+\delta} < \infty$ for some $
\delta >0$,
\item $\beta_t = \beta(t/T)$ is a smooth function such that its second order
derivative is continuous in $[0,1]$,
\item $M(t/T) = E(X_t X_t^\top)$, $\sigma^2(t/T) = E(\varepsilon_t^2)$ and $
V(t/T) = E(X_t X_t^\top \varepsilon_{t+1}^2)$, where $M(\tau)$, $
\sigma^2(\tau)$ and $V(\tau)$ are Lipschitz continuous for all $\tau \in
[0,1]$, and $M(\tau)$ is positive definite.
\end{enumerate}
\bigskip
\noindent \textbf{Assumption A.3 (Forecast error):} Let $\{ \varepsilon_t\}$
be a martingale difference sequence (m.d.s.). In particular, $
E(\varepsilon_{t+1}|\mathcal{I}_t) = 0$, where $\mathcal{I}_t = \{X_t^\top,
X_{t-1}^\top, \ldots, \varepsilon_{t}, \varepsilon_{t-1}, \ldots \}$.
\bigskip
\noindent \textbf{Assumption A.4 (Kernel):} $k: [-1,1] \rightarrow \mathbb{R}
^+$ is a symmetric bounded probability density function. Further, assume $
k(0) \geq k(u)$ for all $u \in [-1,1]$, and $\int k^2(u) du < \infty$.
\bigskip
Assumption A.1 limits the temporal dependence in $\{R_{t}\}$ under a $\beta $-mixing structure. Assumption A.2 imposes common smoothness restrictions on the functions of interest \citep[see][]{cai2007trending, orbe2005nonparametric, robinson1989nonparametric}, and requires slightly more than four moments of the data. More importantly, unlike \citet{cai2007trending} and \citet{chen2012testing}, we allow for time-varying moments which means that the data does not have to be stationary. This is highly relevant because the stationarity of many macroeconomic variables that are of forecasting interest, such as inflation, is still subject to debate. Notably, these conditions are sufficiently general to accommodate nonlinear
locally stationary processes as defined in \citet{dahlhaus2019towards} and \citet{vogt2012nonparametric}, which include time-varying parameter autoregressive processes. Assumption A.3 allows for conditional heteroscedasticity of an unknown form but rules out potential serial correlation in the forecast errors. This condition is reasonable if we expect forecasters to have included all the information known to them at time $t$ in making their forecasts. When more forecasts are used in the combination, the condition is more likely to hold. The m.d.s. assumption greatly simplifies the analysis with high-dimensional data, but we also consider relaxing the assumption to allow for serial correlation below. Lastly, A.4 is a standard assumption for kernel regressions. We note that commonly used second-order kernels, such as the Epanechnikov, uniform and quartic kernels, satisfy this condition. Furthermore, A.4 implies $\int_{-1}^{1}k(u)du=1$, $\int_{-1}^{1}uk(u)du=0$, and $\int_{-1}^{1}u^{2}k(u)du<\infty $.
We now state the asymptotic properties of $\hat{\beta}(\tau)$, which can be viewed as an extension of Theorem 4 of \citet{cai2007trending} for forecast combination with the reflection method.
\begin{proposition}
\thlabel{consis} If assumptions A.1-4 hold, and $h=O(T^{-1/5})$, then for
all $\tau \in \lbrack 0,1]$, we have
\begin{equation}
\sqrt{Th}\bigg[\hat{\beta}(\tau )-\beta (\tau )-\frac{h^{2}\beta ^{^{\prime
\prime }}(\tau )\mu _{2}}{2}+o_{p}(h^{2})\bigg]\rightarrow ^{d}N(0,2\nu
_{0}M^{-1}(\tau )V(\tau )M^{-1}(\tau )).
\end{equation}
where $\nu _{0}=\int_{-1}^{1}k^{2}(u)du$, $\mu _{2}=\int u^{2}k(u)du$, and $
V(\tau)$ is defined in A.2.
\end{proposition}
\thref{consis} shows that $\hat{\beta}(\tau )$ is a consistent estimator of $\beta (\tau)$ and the asymptotic bias depends on the curvature of $\beta(\tau).$ Numerical analysis with commonly used kernels shows that the asymptotic variance here is much smaller than what we might expect if we had not used the data reflection. For example, with the Epanechinikov kernel, the reflection method can reduce the asymptotic variance by more than 70\%.
As \citet{diebold1988serial} points out, regression-based methods of forecast combination might lead to serially correlated errors. Therefore, we relax the m.d.s. assumption and consider the following alternative.
\noindent \textbf{Assumption A.3* (Serial correlation):} Let {$\varepsilon_t$
} satisfy: (i) $E( \varepsilon_{t+1}|X_{t}) =0$, and (ii) $\Gamma
_{j}(t/T)=Cov(X_{t}\varepsilon_{t+1},X_{t+j}\varepsilon _{t+j+1}) $, where $
\Gamma _{j}(\tau)$ is Lipschitz continuous for all $\tau \in [0,1]$.
\begin{proposition}
\thlabel{consis2} If assumptions A.1,A.2,A.3* and A.4 hold, and $
h=O(T^{-1/5})$, then for all $\tau \in \lbrack 0,1]$, we have
\begin{equation}
\sqrt{Th}\bigg[\hat{\beta}(\tau )-\beta (\tau )-\frac{h^{2}\beta ^{^{\prime
\prime }}(\tau )\mu _{2}}{2}+o_{p}(h^{2})\bigg]\rightarrow ^{d}N(0,2\nu
_{0}M^{-1}(\tau )\Omega (\tau )M^{-1}(\tau )). \label{betanorm}
\end{equation}
where $\Omega (\tau )=\sum_{j=-\infty }^{\infty }\Gamma _{j}(\tau )$ .
\end{proposition}
\thref{consis2} shows that asymptotic normality continues to hold with potential serial correlation, while assumptions A.1 and A.2 guarantee the existence of the long-run variance $\Omega(\tau)$ for each $\tau \in [0,1]$.
Next, we study the optimal choice of the bandwidth, $h$, under the m.d.s. assumption (A.3). The choice of the bandwidth $h$ is generally believed to be more important than the choice of the kernel function $k(\cdot )$ in estimation. A small $h$ tends to reduce the bias in $\hat{\beta}
_{t} $ at the expense of variance, and vice versa with large $h$. Hence, we
opt for a data-driven method to select the bandwidth by optimizing some
metric of forecast errors. From \thref{consis}, we obtain the mean squared
combined forecast errors (MSCFE) as
\begin{align}
MSCFE_{t}(h)& =E[(y_{t+1}-X_{t}^{\top }\beta _{t})^{2}] \notag \\
& =E[\varepsilon _{t+1}^{2}]+E[(\hat{\beta}_{t}-\beta _{t})^{\top
}X_{t}X_{t}^{\top }(\hat{\beta}_{t}-\beta _{t})].
\end{align}
Here, we eliminate the cross product term since $\hat{\beta}_{t} $ only uses information up to time $t$. Subsequently, define the integrated MSCFE as
\begin{align}
IMSCFE(h)& =\int_{0}^{1}MSCFE_{t}(h)d(t/T) \notag \label{IMSCFE} \\
& =\int_{0}^{1}\sigma^2 (\tau )d\tau +\int_{0}^{1}\operatorname{Tr}{\bigg[M(\tau)\bigg\{
\frac{h^4\mu_2^2}{4} \beta^{''}(\tau)\beta^{''}(\tau)^\top + \frac{2\nu_0
V_{\beta}(\tau)}{Th}\bigg\}\bigg]}d\tau ,
\end{align}
where $V_{\beta }(\tau )\equiv M^{-1}(\tau )V(\tau )M^{-1}(\tau )$ and label the second term in \eqref{IMSCFE} as $IMSCFE(h)_{L}$. To obtain the optimal bandwidth, we minimize the $IMSCFE(h)$ or equivalently $IMSCFE(h)_{L}$ with respect to $h$, which yields
\begin{equation}
h^{opt}=T^{-\frac{1}{5}}\bigg(\frac{2\nu _{0}\int \operatorname{Tr}\lbrack {V(\tau
)M^{-1}(\tau )]d\tau }}{\mu _{2}^{2}\int \operatorname{Tr}[{M(\tau) \beta^{''}(\tau)
\beta^{''}(\tau)^\top}]d\tau }\bigg)^{\frac{1}{5}},
\end{equation}
and hence the optimal convergence rate of the IMSCFE is of the order $O(T^{-4/5})$.
In practice, we can use a leave-one-out cross-validation (CV) to select the bandwidth. Specifically, a data-driven choice of $h$ is obtained by solving the following problem,
\begin{equation}
\label{CVselect}
\hat{h}_{CV}=\underset{c_{1}T^{-1/5}\leq h\leq c_{2}T^{-1/5}}{\operatorname{arg \ min}}CV(h),
\end{equation}
where $CV(h)=T^{-1}\sum_{s=1}^{T}(y_{s+1}-X_{s}^{\top }\hat{\beta}_{s})^{2}$
, and $c_{1}$ and $c_{2}$ are suitable constants. To formalize the
optimality of the bandwidth selected by CV, we require
stronger moment assumptions.
\noindent \textbf{Assumption A.5 (CV moment condition):} Assume $\sup_{1\leq
t\leq T}E\Vert R_{t}\Vert ^{12}<\infty $.
Assumption A.5 is imposed to facilitate technical derivation. \cite{hardle1985optimal} and \cite{xia2002asymptotic} impose similar moment conditions.
\begin{theorem}
\thlabel{CVunif} Suppose assumptions A.1-5 are satisfied. As $T\rightarrow
\infty $
\begin{equation}
\hat{h}_{CV}/h^{opt}\rightarrow^{p} 1.
\end{equation}
\end{theorem}
The estimated bandwidth derived from minimizing CV is asymptotically optimal in the sense that it minimizes the IMSCFE. Unlike plug-in methods which are directly based on the theoretical optimal bandwidth $h^{opt}$, the CV method does not require the preliminary estimation of asymptotic bias or variance. It can be implemented automatically and we expect that it would have reasonable finite sample performance.
\section{High-Dimensional Forecast Combination}
\label{highdimforecast}
When the dimension of forecasts is large, the local linear estimation does not work well. Moreover in practice, some individual forecasts may or may not be important or relevant. Therefore, we combine the local linear estimation with regularization for model selection and estimation of the combination weights in a high-dimensional context.
Consider
\begin{equation*}
y_{t+1}=\omega _{0t}+f_{t}^{\top }\omega _{1t}+\varepsilon _{t+1},\qquad
t=1,...,T,
\end{equation*}
where the $f_{t}=\left( f_{t1,}f_{t2},...,f_{tp_{T}}\right) ^{\top }$\text{
is }$p_{T}\times 1$ and $p_{T}$, the number of candidate forecasts, can be larger than the sample size $T.$ We assume that there exists $d\ll T$ and $1\leq d<p_{T}$\text{ such that }$
\omega _{1t,j}\neq 0$ for $1\leq j\leq d$ and $\omega _{1t,j}=0$
for $d<j\leq p_{T}.$ In other words, there are $d$ relevant forecasts. Moreover, the dimension of important individual forecasts $d$ may diverge with $T$. A natural
estimator of forecast weights would be the local linear estimator with the
Lasso penalty:
\begin{equation}
\tilde{\gamma}_{t}=\underset{\gamma _{t}\in \mathbb{R}^{2(p_{T}+1)}}{\arg
\min }T^{-1}\sum_{s=t-\floor* {Th},s\neq t}^{t+\floor* {Th}}k_{st}
\left[ y_{s+1}-\alpha _{0t}^{\top }X_{s}-\alpha _{1t}^{\top }\left( \frac{s-t
}{T}\right) X_{s}\right] ^{2}+\lambda _{1}\left\vert \alpha _{0t}\right\vert
+\lambda _{2}\left\vert h\alpha _{1t}\right\vert , \label{prelim}
\end{equation}
where $\lambda _{1}$ and $\lambda _{2}$ are two tuning parameters, and $
\gamma _{t}=(\alpha _{0t}^{\top },\alpha _{1t}^{\top })^{\top }$.
Then, the local linear estimator for $\beta _{t}$ is given by
\begin{equation}
\tilde{\beta}_{t}=(e_{1}^{\top }\otimes I_{(p_{T}+1)})\tilde{\gamma}_{t}.
\label{betat}
\end{equation}
Define $\tilde{\Upsilon}=\left( \tilde{\gamma}_{1},\tilde{\gamma}_{2},...,
\tilde{\gamma}_{T}\right) ^{\top }$ and $\tilde{B}=\left( \tilde{\beta}_{1},
\tilde{\beta}_{2},...,\tilde{\beta}_{T}\right) ^{\top }$. It is easy to
verify that
\begin{align}
\tilde{\Upsilon}=& \underset{\Upsilon \in \mathbb{R}^{T\times 2(p_{T}+1)}}{
\arg \min }T^{-1}\sum_{t=1}^{T}\sum_{s=t-\floor* {Th},s\neq t}^{t+\lfloor
Th\rfloor }k_{st}\left[ y_{s+1}-\alpha _{0t}^{\top }X_{s}-\alpha _{1t}^{\top
}\left( \frac{s-t}{T}\right) X_{s}\right] ^{2} \notag \\
& +\lambda _{1}\left\vert \alpha _{0t}\right\vert +\lambda _{2}\left\vert
h\alpha _{1t}\right\vert . \label{fslasso}
\end{align}
The Lasso-based local linear estimator is estimation consistent (see \thref{fsprop}) but requires very strong assumptions for selection consistency. Instead, following \cite{li2015model}, we minimize \eqref{fslasso} to obtain preliminary estimates for use in a second-stage penalized optimization with the (group) SCAD penalty proposed by \cite{fan2001variable}. In particular, we use the initial
estimates from $\tilde{B}$ to solve
\begin{eqnarray}
\hat{\Upsilon}^{h} &=&\underset{\Upsilon \in \mathbb{R}^{T\times 2(p_{T}+1)}}
{\arg \min }T^{-1}\sum_{t=1}^{T}\sum_{s=t-\floor* {Th},s\neq t}^{t+\floor* {Th}}k_{st}\left[ y_{s+1}-\alpha _{0t}^{\top }X_{s}-\alpha _{1t}^{\top
}\left( \frac{s-t}{T}\right) X_{s}\right] ^{2} \label{llscad} \\
&&+\sum_{j=1}^{p_{T}+1}p_{\lambda _{3}}^{^{\prime }}\left( \left\Vert \tilde{
B}_{j}\right\Vert \right) \left\Vert \alpha _{0,j}\right\Vert
+\sum_{j=1}^{p_{T}+1}p_{\lambda _{4}}^{^{\prime }}\left( \tilde{D}
_{j}\right) \left\Vert h\alpha _{1,j}\right\Vert , \notag
\end{eqnarray}
where $\lambda _{3}$ and $\lambda _{4}$ are two tuning parameters, $\alpha
_{i}=\left( \alpha _{i1},\alpha _{i2},...,\alpha _{iT}\right) ^{\top }$ for $
i=0,1,$ and $\alpha _{i,j}$ is the $j^{th}$ column of $\alpha _{i},$ $\tilde{
B}_{j}$ is the $j^{th}$ column of the first-stage Lasso-based local linear
estimator $\tilde{B},$ and
\begin{equation*}
\tilde{D}_{j}=\left\{ \sum_{t=1}^{T}\left[ \tilde{\beta}_{t,j}-\frac{1}{T}
\sum_{s=1}^{T}\tilde{\beta}_{s,j}\right] ^{2}\right\} ^{1/2},
\end{equation*}
which measures the smoothness of the LASSO-based local linear estimator $
\tilde{B}.$ Moreover, $p_{\lambda }^{'}\left( \cdot \right) $ is the
derivative of the SCAD penalty function with regularization parameter $
\lambda $ defined by
\begin{equation*}
p_{\lambda }^{' }\left( x\right) =\lambda \left[ 1\left( x\leq \lambda
\right) +\frac{\left( a\lambda -x\right) _{+}}{\left( a-1\right) \lambda }
1\left( x>\lambda \right) \right]
\end{equation*}
and $a=3.7$ as suggested in \cite{fan2001variable}. Instead of the SCAD
penalty itself, we use a local linear approximation of the penalty to
overcome difficulties due to non-convexity of the SCAD penalty
\citep{zou2008one, fan2014strong}. Then the local linear estimator for $
\beta _{t}$ with the group SCAD penalty is
\begin{equation}
\hat{\beta}_{t}^{h}=(e_{1}^{\top }\otimes I_{(p_{T}+1)})\hat{\gamma}_{t}^{h},
\label{betah}
\end{equation}
where $\hat{\gamma}_{t}^{h}$ is the $t^{th}$ row of $\hat{\Upsilon}^{h}.$
\section{Computational Algorithm}
\label{notation}
We approach the estimation of our forecast weights with a two-stage strategy\footnote{Throughout the process, we standardize our data even though the forecasts and the variable of interest are expected to share the same scale because we find that it helps with the stability of the algorithm and is standard
practice in the Lasso literature. This involves centering the data and dividing by its standard deviation using the whole sample.}. In the first
stage, we solve \eqref{prelim} to obtain preliminary coefficient estimates. Subsequently, we use these to initialize the group coordinate descent algorithm (Yuan and Lin, 2006; Wei et al., 2011) in order to solve the penalized regression with the group SCAD penalties. The first-stage is a standard problem that can be solved efficiently for each $t$ by accessible statistical programs, such as {\fontfamily{lmtt}\selectfont glmnet} in {\fontfamily{lmtt}\selectfont R}, so we focus our attention to the second-stage problem with the group SCAD penalties.
For convenience, we rewrite the first term of \eqref{llscad} in matrix
notation as
\begin{equation}
\mathcal{L}^\diamond(\alpha_0, \alpha_1) = T^{-1} \bigg(\overline{Y}
-\sum_{j=1}^{p_{T}+1}\Xi _{j}\alpha _{0,j}-\sum_{j=1}^{p_{T}+1}\Xi
_{j+p_{T}+1}\alpha _{1,j}\bigg)^{\top }\overline{K}\bigg(\overline{Y}
-\sum_{j=1}^{p_{T}+1}\Xi _{j}\alpha _{0,j}-\sum_{j=1}^{p_{T}+1}\Xi
_{j+p_{T}+1}\alpha _{1,j}\bigg) \label{llmatrix}
\end{equation}
where $\alpha_{i,j}$ is as previously, the $T \times 1$ $j^{th}$ column of $
\alpha_i$ for $i=0,1$. $\overline{Y}$ is a $2\floor*{Th}T\times 1$ vector
obtained by stacking $Y_{t}$ for $t=1,\ldots ,T$ which is in turn a $2\floor*
{Th} \times 1$ vector obtained by stacking $y_{s+1}$ for $s$ from $t-\floor*{
Th}$ to $t+\floor*{Th}$ excluding $t$. $\overline{K}$ is a $2\floor*{Th}
T\times 2\floor*{Th}T$ block-diagonalization of $\{K_{t}\}_{t}^{T}$, where $
K_t$ is a diagonal matrix with diagonal elements corresponding to $k_{t-
\floor*{Th},t},\ldots,k_{t+\floor*{Th},t}$ excluding $k_{t,t}$. $\Xi_{i}$ is
a selection matrix such that
\begin{equation}
\underset{2\floor*{Th}T \times T}{\Xi_{i}}=
\begin{bmatrix}
e_{1}^{\top }\otimes Q_{1}e_{i,2(p_{T}+1)} \\
\vdots \\
e_{T}^{\top }\otimes Q_{T}e_{i,2(p_{T}+1)}
\end{bmatrix}
\label{Xi}
\end{equation}
where $Q_t$ is obtained by vertically stacking $(X_s^\top, X_s^\top(\frac{s-t
}{T}))$ for $s$ from $t-\floor*{Th}$ to $t+\floor*{Th}$ excluding $t$, and $
e_t$ and $e_{i,2(p_T+1)}$ are $T \times 1$ and $2(p_T+1) \times 1$ unit
vectors with unity in the $t^{th}$ and $i^{th}$ coordinates respectively.
To use the group coordinate descent algorithm as in \cite{wei2011variable}
and \cite{yuan2006model}, we need to orthogonalize the matrix $\Xi
_{i}^{\top } \overline{K} \Xi _{i}$. This can be achieved by
post-multiplying $\Xi_i$ with the inverse of the Cholesky decomposition of $
\Xi _{i}^{\top } \overline{K} \Xi _{i}$ (let it be $A_i$), so that $(\Xi_i
A_{i})^{\top } \overline{K} (\Xi_i A_{i})=I$. Hence, we proceed with the
assumption that $\Xi _{i}$ has been orthogonalized.
Then, it can be shown that
\begin{equation*}
\alpha _{0,i}=\bigg(1-\frac{\tau _{i}}{\Vert \tilde{S}_{i}\Vert }\bigg)_{+}
\tilde{S}_{i}\quad \text{and}\quad \alpha _{1,i}=\bigg(1-\frac{h\tau
_{i}^{\ast }}{\Vert \tilde{S}_{i}^{\ast }\Vert }\bigg)_{+}\tilde{S}
_{i}^{\ast }
\end{equation*}
where $\tau _{i}=p_{\lambda _{3}}^{' }\left( \left\Vert \tilde{B}
_{i}\right\Vert \right) ,\tau _{i}^{\ast }=p_{\lambda 4}^{' }\left(
\tilde{D}_{i}\right)$, and $\tilde{S}_{i}=\Xi_{i}^{\top } \overline{K}(
\overline{Y}-\sum_{j \neq i} \Xi_{j} \alpha_{0j}-\sum_{j=1}^{p_{T}+1}
\Xi_{j+p_{T}+1} \alpha_{1j})$ and $\tilde{S}_{i}^{\ast
}=\Xi_{i+p_{T}+1}^{\top }\overline{K}(\overline{Y}-\sum_{j=1}^{p_{T}+1}
\Xi_{j}\alpha _{0,j}-\sum_{j\neq i} \Xi_{j+p_{T}+1}\alpha _{1,j})$. These
relationships allow us to iteratively compute the parameters through the
following algorithm.
\smallskip
Step 1. Initialize with estimates from the first-stage lasso.
In other words, $\alpha_{0,i}^{(0)} = \tilde{\alpha}_{0,i}$ or $
\alpha_{1,i}^{(0)} = \tilde{\alpha}_{1,i}$. Define $r_0^{(0)} = r_1^{(0)} = \overline{Y}$.
\smallskip
Step 2. Construct $\tilde{S}_i^{(k+1)} = \Xi_i^\top \overline{K
} r_0^{(k)} + \alpha_{0,i}^{(k)}$ or $\tilde{S}_i^{*(k+1)} = \Xi_{i+d+1}^\top
\overline{K} r_1^{(k)} + \alpha_{1,i}^{(k)}$.
\smallskip
Step 3. Use the relation $\alpha_{0,i}^{(k+1)} = (1 - \frac{
\tau_i^{(k)}}{\| \tilde{S}_i^{(k)} \|})_+ \tilde{S}^{(k+1)}_i$ or $
\alpha_{1,i}^{(k+1)} = (1 - \frac{h\tau_i^{*(k)}}{\| \tilde{S}_i^{*(k)} \|}
)_+ \tilde{S}_i^{*(k+1)}$ to update, and use these to form $\tau_i^{(k+1)}$
and $\tau_i^{*(k+1)}$.
\smallskip
Step 4. Construct a new $r$ with $r_0^{(k+1)} = r_0^{(k)} -
\Xi_i^\top \overline{K} ( \alpha_{0,i}^{(k+1)} - \alpha_{0,i}^{(k)})$ or $
r_1^{(k+1)} = r_1^{(k)} - \Xi_{i+d+1}^\top \overline{K} ( \alpha_{1,i}^{(k+1)} -
\alpha_{1,i}^{(k)})$.
\smallskip
Step 5. Repeat for all the coefficient vectors until a
reasonable tolerance is achieved. We use $1\times10^{-3}$ in our simulations
and applications.
\smallskip
Step 6. Recover the original parameters by applying the
reverse transformation for orthogonalization and standardization.
To implement the local linear estimation with the group SCAD penalty, we need to choose tuning parameters $\lambda _{j}$ and $h$. First, for the preliminary estimates, the tuning parameters $\lambda_{1}$ and $\lambda _{2}$ are obtained via K-fold CV. Given the potentially high computational costs involved in our two-step procedure, CV is attractive because it is readily accessible as the default option of many Lasso-type algorithms in statistical programs\footnote{For example, {\fontfamily{lmtt}\selectfont glmnet} and {\fontfamily{lmtt}\selectfont lars} in {\fontfamily{lmtt}\selectfont R}}. Furthermore, in a high-dimensional setting, \cite{homrighausen2017risk} have shown that the CV estimate for Lasso is risk consistent for the oracle tuning parameter. This result might be stronger than what we require in section \ref{high.asymp} since we do not expect the preliminary estimator to be selection consistent for our asymptotic results.
Following \cite{li2015model}, we set the bandwidth as $h=C[\log(p_{T}+1)/T]^{1/5}$ to minimize the computational burden from additional tuning. Simulation studies show that the preliminary
estimation is not very sensitive to bandwidth selection. For the second-stage estimation, the tuning parameters $\lambda_{3}$ and $\lambda_{4}$ are selected via a modified version of BIC:
\begin{equation*}
BIC=\log (SSR)+C_{T}\left( l\log \floor*{Th}/\floor*{Th}\right) ,
\end{equation*}
where $SSR=T^{-1}\sum_{t=1}^{T}(y_{t+1}-x_{t}^{\top }\hat{\beta}
_{t}^{h})^{2},$ $l$ is the number of significant forecasts (i.e. maximum possible value of $l$ is $d$) obtained given a pair of candidate tuning parameters, and $\lfloor Th\rfloor $ is the effective sample size for estimating time-varying parameters. $C_{T}=\log p_{T}$ is selected to guarantee the consistency of BIC in a high-dimensional regression, following \cite{wang2009shrinkagetuning}. Here, BIC is a popular tuning approach for variable selection problems with SCAD penalties. For example, \cite{cai2015functional} use the BIC for tuning the SCAD penalty in penalized functional coefficient models, although we remark that the tuning parameters we seek, $\lambda_3$ and $\lambda_4$, are less restrictive than theirs\footnote{In particular, we require $\lambda_3 \propto \lambda_4 = o(T^{1/2})$ in HD.3(ii), while $\lambda = o(T^{1/10})$ in \cite{cai2015functional}}.
\section{Asymptotic Analysis with High Dimension}
\label{high.asymp}
To study the asymptotic properties of $\tilde{\beta}_{t}$ and $\hat{\beta}
_{t}^{h},$ we impose the following additional assumptions using the notation established in the previous section.
\bigskip
\noindent \textbf{Assumption HD.1 (Moment conditions):} \setlist{nolistsep}
\begin{enumerate}[label=(\roman*), noitemsep]
\item Let $X_t^{o}$ contain the first $d$ relevant forecasts. Define $M^{o}(t/T) = E(X_t^{o} X_t^{o^\top})$, $\sigma^2(t/T) = E(\varepsilon_t^2)$ and $
V^{o}(t/T) = E(X_t^{o} X_t^{o^\top} \varepsilon_{t+1}^2)$, where $M^{o}(\tau)$, $\sigma^2(\tau)$ and $V^{o}(\tau)$ are Lipschitz continuous for all $\tau \in
[0,1]$, and $M^{o}(\tau)$ is positive definite.
\item Assume that uniformly in $t$, $\max_{1\leq j\leq
(p_{T}+1)}E[|\varepsilon _{t+1}X_{tj}|^{\psi }]<\infty $ for some sufficiently large $\psi >2+
\frac{\delta _{1}}{1-\delta _{2}}+\delta $, $\delta >0,$ where $\delta _{1}$
and $\delta _{2}$ are defined in assumption HD.3 below, and $X_{sj}$ refers to the $
j^{th}$ element of $X_{s}$.
\end{enumerate}
\bigskip
\noindent \textbf{Assumption HD.2 (Restricted eigenvalues):} Define the set $
S$ as
\begin{align*}
S=\bigg\{ &v=(v_{11},\ldots,v_{1p_T+1},v_{21},\ldots,v_{2p_T+1})^\top \in
\mathbb{R}^{2(p_T+1)}: \|v\|=1, \\
&\sum_{m=1}^{p_T+1}(|v_{1m}|+|v_{2m}|)\leq 2(1+\delta)\sum_{m=1}^{d}
(|v_{1m}|+|v_{2m}|) \bigg\},
\end{align*}
for some $\delta>0$. Then there exists positive constants $
0<\rho_1\leq\rho_2<\infty$ with probability approaching one such that,
\begin{equation*}
\rho_1 \leq T^{-1} \inf_{\tau \in [0.1]} \inf_{v \in S} v^\top
Q_\tau^{h\top} K_\tau Q_\tau^h v \leq T^{-1} \sup_{\tau \in [0.1]} \sup_{v
\in S} v^\top Q_\tau^{h\top} K_\tau Q_\tau^h v \leq \rho_2
\end{equation*}
where $Q^h_\tau= Q_\tau H$, $H=diag\{I_{p_T+1\times p_T+1}, \mathbf{h}
^{-1}\} $, and $\mathbf{h}^{-1}$ is a $p_T+1\times p_T+1$ diagonal matrix
with diagonal elements $1/h$.\footnote{
Effectively, $Q_t^h$ is the matrix constructed by vertically stacking $
(X_s^\top,X_s^\top(\frac{t-s}{Th}))$ from $s = t-\floor*{Th}$ to $t+\floor*{
Th}$ excluding $t$.} $Q_\tau$ and $K_\tau$ are defined in the discussion of
\eqref{llmatrix} in Section \ref{notation}.
\noindent \textbf{Assumption HD.3 (Rates and tuning parameters):}
\setlist{nolistsep}
\begin{enumerate}[label=(\roman*), noitemsep]
\item Let $p_{T}=c_{1}T^{\delta _{1}}$, $h=c_{2}T^{-\delta _{2}}$, $\lambda
_{1}\varpropto \lambda _{2}$, where $0\leq \delta _{1}<\infty $, $0<\delta
_{2}<1.$ The bandwidth and the tuning parameter $\lambda _{1}$ satisfy $
d h^{2}\lambda _{1}^{-1}+d h^{-2}\lambda_{1}^{2}+ d\lambda _{1}^{1/2}+\left( \log h^{-1}/Th\right) ^{1/2}\lambda _{1}^{-1}\rightarrow 0$
.
\item Let $dh^2 \varpropto (Th)^{-1/2}$, $\lambda _{3}\varpropto \lambda
_{4}$, $\lambda_3 = o(T^{1/2})$, and $h^{-1/2}[(\log h^{-1})^{1/2}+
d^{1/2} + \lambda_1 h^{1/2}\sqrt{Td}]\lambda _{3}^{-1} \rightarrow 0.$
\item With probability approaching one, there exists a positive constant $
b_\diamond$ such that
\begin{equation*}
\min_{1 \leq j \leq d} \|B_j\| \geq b_\diamond T^{1/2},\text{ and } \min_{1
\leq j \leq d_1} D_{j} \geq b_\diamond T^{1/2}.
\end{equation*}
\end{enumerate}
\bigskip
Assumption HD.1(i) is the high-dimensional counterpart of (iii) in assumption 2. HD.1(ii) is a moment condition similar to that in \cite{li2015model}, while HD.2 is a generalization of the restricted eigenvalue conditions in \cite{bickel2009simultaneous}. The regularity conditions in HD.3(i) allows the number of forecasts to increase at a polynomial rate and imposes restrictions on the penalty parameters, $\lambda _{1}$ and $\lambda _{2}$ and the bandwidth, $h$, which are required for showing \thref{fsprop}. On the other hand, HD.3(ii) allows $\lambda _{3}\rightarrow \infty $, albeit at a slower rate than $\sqrt{T}$, and is used to show \thref{ssoracle}. Finally, HD.3(iii) requires the coefficients on relevant forecasts to be bounded away from 0, which is used to show the oracle property in \thref{oracleprop}.
We first establish the asymptotic properties of the first-stage estimator $
\tilde{\beta}_{t}.$
\begin{proposition}
\thlabel{fsprop} If assumptions A.1,A.2(i)-(ii),A.3, A.4, HD.1-2, and HD.3(i) hold, we have
\begin{equation}
\max_{t}\Vert \tilde{\beta}_{t}-\beta _{t}\Vert \rightarrow^{p} 0
\end{equation}
as $T\rightarrow \infty .$
\end{proposition}
\thref{fsprop} shows that the local linear estimator with Lasso penalty is estimation consistent but it is not variable selection consistent in the absence of strong "irrepresentability" conditions for the Lasso in a high-dimensional setting \citep[see for e.g.][]{zhang2010nearly}. Therefore, we only use $\tilde{\beta}_{t}$ as the first-stage estimator to figure out the initial weights for use in the group SCAD penalty. To study the selection consistency of
the proposed two-stage method, we define $S=\left\{ j_{1},...,j_{d^{\ast
}}\right\} $ as the index set of an arbitrary model with a total of $0\leq d^{\ast }\leq p_{T}
$ non-zero coefficients (i.e.$X_{tj_{1}},...,X_{tj_{d^{\ast }}}).$ Then we
use $S_{0}=\left\{ 1,...,d\right\} $ to denote the index set of the true model and $\hat{S}
=\{ j:\Vert \hat{B}_{j}^{h}\Vert >0\} $ to represent
the model selected by the two-stage procedure, where $\hat{B}_{j}^{h}=( \hat{\beta}_{1,j}^{h},...,
\hat{\beta}_{T,j}^{h}) ^{\top }.$
\begin{theorem}
\thlabel{ssoracle} Assume assumptions A.1, A.2(i)-(ii), A.3, A.4, and HD.1-3
hold. We have
\begin{equation} \label{select}
P\left( \hat{S}=S_{0}\right) \rightarrow 1
\end{equation}
and
\begin{equation} \label{est}
\max_{t} \|\hat{\beta}_t^{h} - \beta_t\| \rightarrow^{p} 0
\end{equation}
as T$\rightarrow \infty .$
\end{theorem}
\thref{ssoracle} shows that the two-stage method is not only estimation consistent but can also consistently select all relevant individual forecasts. Next, we establish the oracle property. Let $\hat{\beta_t}^{o,h}$ and $\beta_t^{o}$ represent the first $d$ non-zero forecast weights in $\hat{\beta_t}^h$ and $\beta_t$ respectively. The following theorem states that the two-stage group SCAD estimator is asymptotically normal.
\begin{theorem}
\thlabel{oracleprop} Assume A.1, A.2(i)-(ii), A.3, A.4, and HD.1-3 hold. Then, for all $\tau \in [0,1]$,
\begin{equation}
\sqrt{Th} A_T \Omega^{o^{-1/2}}(\tau)\bigg\{ \hat{\beta}^{o,h}(\tau) - \beta^{o}(\tau) - \frac{h^2}{2} \mu_2 \beta^{o^{''}}(\tau) + o_p(h^2) \bigg\} \rightarrow^{d} N(0,G)
\end{equation}
as $T\rightarrow \infty$, where $\Omega^o(\tau) = 2 \nu_0 M^{o^{-1}}(\tau)V^{o}(\tau)M^{o^{-1}}(\tau)$, and $A_T$ is an arbitrary $q \times d$ matrix\footnote{Similar to \cite{fan2004nonconcave}, we consider the asymptotic normality of arbitrary linear combinations of $\hat{\beta}^{o,h}(\tau)$ by pre-multiplying $A_T$ because the dimensions depend on $T$ and might diverge as $T$ goes to infinity.} such that $A_T A_T^\top \rightarrow G$ for a given finite $q$
\end{theorem}
\thref{oracleprop} is related to the oracle property because it is derived by showing that the two-stage SCAD estimator is asymptotically equivalent to the oracle estimator\footnote{The estimator obtained from optimizing the likelihood while having \textit{a priori} knowledge of the relevant estimators and excluding the irrelevant ones.}, which implies that both estimators share the same asymptotic behavior. By further establishing the asymptotic normality of the oracle estimator, we can thus translate the property to our estimator and adopt similar tools for statistical inference. We note that in the special case where $d$ does not depend on $T$, the oracle estimator will be identical to the low-dimensional estimator studied in section \ref{sect.asym}, and a direct application of \thref{consis} is possible.
\section{Monte Carlo Simulation}
In this section, we contrast the out-of-sample forecasting performance of
the nonparametric estimator with that of common forecast combination
techniques, which include both static and time-varying approaches. We first look at the case with only two forecasts. Then we study a high-dimensional setting to evaluate the
finite-sample oracle properties of the proposed estimator with the group
SCAD penalty.
\subsection{Forecast combination with low-dimensional data}
\label{lowsim}
To start off, we consider the following time-varying coefficient model:
\begin{equation*}
y_{t+1} = \omega_{0t} + \omega_{1t} f_{1,t} + \beta_{2t} f_{2,t} + u_{t+1}
\end{equation*}
where $\omega_{0t} = \exp(-3 + 2.5 \tau)$, $\omega_{1t} = 0.5(1.5
\tau - 0.8)^3 + 0.5$, and $\omega_{2t} = 0.2\sin(4 \tau) + 0.4$
, for $\tau = t/T$, while the forecasts evolve according to:
\begin{align*}
&f_{1,t} = 0.5 + 0.8y_t + e_{1,t} \\
&f_{2,t} = 0.5 + 0.3\sin\bigg(2\tau + 0.25\bigg)y_t + e_{2,t}.
\end{align*}
Lastly, $u_{t+1}$, $e_{1,t}$ and $e_{2,t}$ are normally distributed with
mean 0 and unit variance.
We study this model specification because it can be endowed with an economic
interpretation consistent with the time-varying common factor framework for
forecast combinations in \cite{elliott2005optimal}. Specifically, it can be
shown that this model reduces to a time-varying parameter AR(1) process with
heteroskedastic errors, and hence in this case, the factor is observable and
completely captured by the first lag of $y$.
The experiment is conducted with three samples: $T \in \{200,300,500\}$ and
an out-of-sample period of 50 time points in excess of \textit{T}. For each
case, we require a holdout or burn-in period of $2 \times T$ for the process
to stabilize and to provide points for bandwidth CV.
In addition, we compare the proposed nonparametric estimator to alternative
approaches that can be classified as either static or adaptive. For the
static case, we use $T$ points to estimate the weights, while we employ an
expanding window ($T + k$) for the adaptive estimators as $k$ increases. The
adaptive estimators suffer from the same issue as that of the nonparametric
estimator in that $y_{t+1}$ is not available at time $t$ for estimation. Hence, we estimate $\hat{\beta}_{t-1}$ and use it to approximate $\hat{\beta}_t$ to form a forecast of $y_{t+1}$. This reflects the practice of professional forecasters and is justified by our assumption that $\beta_t$ is smooth.
Finally, we evaluate the performance of all the models by calculating the
average squared combined forecast error (ASCFE), which is defined as the sum
of the squared deviations of the computed $\hat{y_t}$ from $y_t$ obtained
from the out-of-sample period, i.e. $ASCFE = 1/50\sum_{t=T+2}^{T+51}(y_t-
\hat{y}_t)^2$. The experiment is repeated 500 times to obtain the mean and
standard deviations of the ASCFE.
\subsubsection{Competing forecast combination methods}
\label{competing} \smallskip \noindent \textbf{Nonparametric estimation}
\smallskip
We consider the local linear estimator with data reflection as established in section \ref{sect.asym}. For bandwidth selection, we adopt the CV method introduced in \eqref{CVselect}. For ease of presentation, we label the nonparametric data
reflection method "NPRf". In all subsequent simulations and empirical
applications, we use the Epanechnikov kernel, i.e. $K(u)=0.75(1-u^{2})_{+}$.
\bigskip \noindent \textbf{Bates and Granger (1969)} \smallskip
A commonly used time-varying combination scheme, originally
suggested by \cite{bates1969combination}, is an adaptive updating method which
assigns a higher weight to forecasts that perform comparatively well in the
recent past. In particular, we consider an expanding window version of their
estimator such that the combination weight of forecast $1$ is given by
\begin{equation*}
\omega _{1t}^{BG}=\frac{\hat{e}_{1,t}^{-1}}{\hat{e}_{1,t}^{-1}+\hat{e}_{2,t}^{-1}}
\end{equation*}
where $\hat{e}_{i,t}=\frac{1}{t}\sum_{l=1}^{t}(y_{l+1}-f_{i,l})^{2}$. This is
defined analogously for $f_{2}$. Note that there is no intercept in the
model, and the weights sum to 1. We term this method "BG".
\bigskip \noindent \textbf{Least squares regression} \smallskip
We consider three static combination schemes due to \cite{granger1984improved}. They are:
\begin{align*}
& y_{t+1}=\omega _{0}+\omega _{1}f_{1,t}+\omega _{2}f_{2,t}+z_{t}, \\
& y_{t+1}=\omega _{1}f_{1,t}+\omega _{2}f_{2,t}+z_{t}, \\
& y_{t+1}=\omega _{1}f_{1,t}+\omega _{2}f_{2,t}+z_{t},\text{where }\omega
_{1}+\omega _{2}=1.
\end{align*}
These models are respectively labeled as "GRregconst" for regression with a
constant, "GRreg" for regression without the intercept term, and "GRconstr"
for constrained regression. It is known that despite biased
forecasts, "GRregconst" performs favorably in terms of MSCFE because $\omega _{0}$ is able to capture the bias \cite[see][]{timmermann2006forecast}. Since we have introduced time-variation in our
experiment, we also consider adaptive versions of the models above, whose
weights are estimated ex-ante. These models are named "TVGRregconst",
"TVGRreg", and "TVGRregconstr" accordingly.
\bigskip \noindent \textbf{Equal weights} \smallskip
Lastly, we include the simple strategy of assigning equal weights
to the forecasts, which in this case is expressed by: $\hat{y}_{t+1} =
0.5f_{1,t} + 0.5f_{2,t}$. This is intended to assess whether the proposed
estimators suffer from the "forecast combination puzzle", which refers to
the commonly observed empirical fact that weights derived from a simple
arithmetic mean often outperform theoretical optimal weights based on
sophisticated estimation. We label this "EQ".
\subsubsection{Simulation results}
Table \ref{tab:table1} reports the mean and standard deviation of the ASCFE.
We can see that our estimator outperforms alternative methods by achieving
both the lowest ASCFE, and the smallest variance, at all sample sizes.
\begin{table}[tbp]
\centering
\caption{Simulation results with low-dimensional data (2 forecasts).}
\setlength{\tabcolsep}{20pt}
\begin{threeparttable}
\begin{tabular}{lccc}
\toprule & \multicolumn{3}{c}{Sample size (T)} \\
Estimation & 200 & 300 & 500 \\
\midrule \textit{Adaptive} & & & \\
NPRf & 1.06 & 1.06 & 1.03 \\
& (0.22) & (0.21) & (0.21) \\
BG & 1.19 & 1.22 & 1.20 \\
& (0.24) & (0.23) & (0.24) \\
TVGRregconst & 1.08 & 1.09 & 1.07 \\
& (0.22) & (0.21) & (0.22) \\
TVGRreg & 1.13 & 1.14 & 1.12 \\
& (0.23) & (0.22) & (0.23) \\
TVGRregconstr & 1.17 & 1.19 & 1.17 \\
& (0.23) & (0.22) & (0.23) \\
\textit{Static} & & & \\
GRregconst & 1.09 & 1.11 & 1.08 \\
& (0.24) & (0.22) & (0.22) \\
GRreg & 1.14 & 1.15 & 1.13 \\
& (0.23) & (0.22) & (0.23) \\
GRregconstr & 1.18 & 1.20 & 1.18 \\
& (0.24) & (0.23) & (0.23) \\
EQ & 1.29 & 1.34 & 1.34 \\
& (0.26) & (0.25) & (0.26) \\
\midrule Best & NPRf & NPRf & NPRf \\
\bottomrule & & &
\end{tabular}
\begin{tablenotes}[flushleft]
\small
\item Notes: (1) Mean and standard deviation in parentheses of ASCFE from 500 iterations. (2) NPRf: nonparametric estimator with data reflection, BG: \cite{bates1969combination} adaptive etimator, GRregconst: OLS regression with intercept, GRreg: OLS regression without intercept, GRregconstr: OLS regression with sum of coefficients constrained to unity, TV: time-varying weights, EQ: equal weights.
\end{tablenotes}
\end{threeparttable}
\label{tab:table1}
\end{table}
In addition, two observations are salient. First, adaptive estimators
perform better than their static counterparts, which is expected because the
true DGP contains many time-varying parameters. In
addition, within the two categories, approaches that allow for an
intercept (TVGRregconst and GRregconst) achieve relatively low ASCFEs which
is likely attributable to their ability in capturing the (time-varying) bias
in the true model. By this logic, it is reasonable that NPRf performs well because it
accommodates both time-varying biases and weights. The lackluster
performance of EQ is also consistent with this reasoning as it permits
neither.
For robustness, we look at 3 more cases in the online appendix. The cases are
modifications of the original DGP: (1) bias is removed
and time-varying weights replaced with constant weights; (2) a constant bias
is permitted and time-varying weights replaced with constant weights; and
(3) bias is removed but weights remain time-varying. The key message from
this exercise is that NPRf performs the best when forecast weights are
indeed changing over time. Regression-based methods appear to perform
slightly better otherwise, although NPRf is a close runner-up.
\subsection{Forecast combination with high-dimensional data}
\label{highsim}
Now we consider the performance of NPRf with the group SCAD penalty in
situations where the number of forecasts may be larger than the effective
sample size. The DGP of section \ref{lowsim} is extended by
considering $J$ additional forecasts $\{f_{j,t}\}_{j=1}^{J}$, for each $t$,
where $J\in \{10,50,100\}$. However, these forecasts are redundant in that $
\beta _{jt}=0$ for all $j=1,\ldots ,J$ and all $t$. The redundant forecasts
are generated from a joint normal distribution with mean \textbf{0} and
variance-covariance matrix $\Sigma $ such that $cov(f_{j,t},f_{j^{'},t})=2\exp (-|j-j^{^{' }}|)$ and $j,j^{^{'}}=1,\ldots ,J$. To
minimize the computational burden, we consider an out-of-sample period of 10
points for the sample sizes 50,100, and 150, and with 200 Monte Carlo
simulations.
As mentioned in section \ref{notation}, we choose the respective penalty parameters via K-fold CV in the first-stage and BIC in the second. The bandwidth is selected using a rule as in \cite{li2015model}: $h=(\log
(J+3)/T)^{0.2}$. The initial values used in the
group coordinate descent algorithm are obtained from a first-stage Lasso.
Since NPRf is a local linear estimator, the effective number of regressors
(EReg) is thus $2 \times (J + 3)$ where the addition of 3 refers to the
intercept and the two relevant forecasts. For large $J$, EReg often exceeds
the effective sample size, especially when the bandwidth is small, and the
solution to the weighted least squares problem implicit in NPRf is no longer
unique\footnote{Given that EReg is often larger than the sample size, we do not consider the alternative methods from section \ref{competing} as they would no longer produce reliable estimates}. Hence the goal of this analysis is to assess whether the
penalization scheme can accurately select the relevant forecasts and
estimate their weights.
\begin{table}
\centering
\caption{Simulation results for high-dimensional data.}
\begin{threeparttable}
\begin{tabular}{rrrr}
\toprule
\multicolumn{1}{l}{Sample size} & \multicolumn{3}{c}{\textbf{group SCAD}} \\
\cmidrule{2-4} & \multicolumn{1}{c}{ASCFE} & \multicolumn{1}{c}{\% correct for OOS period} & \multicolumn{1}{c}{Relevant included for OOS period} \\
\midrule
\multicolumn{1}{l}{10 extra forecasts} & & & \\
50 & 1.55 & 0.81 & 1.00 \\
100 & 1.34 & 0.91 & 1.00 \\
150 & 1.34 & 0.96 & 1.00 \\
\multicolumn{1}{l}{50 extra forecasts} & & & \\
50 & 1.53 & 0.73 & 0.98 \\
100 & 1.40 & 0.90 & 1.00 \\
150 & 1.25 & 0.91 & 1.00 \\
\multicolumn{1}{l}{100 extra forecasts} & & & \\
50 & 1.62 & 0.70 & 0.96 \\
100 & 1.34 & 0.80 & 1.00 \\
150 & 1.24 & 0.81 & 1.00 \\
& & & \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\small
\item Notes: (1) Mean of ASCFE and share of iterations that selected exactly or included the relevant forecasts for the whole out-of-sample (OOS) period of 10 time points. (2) Results are from 200 iterations.
\end{tablenotes}
\end{threeparttable}
\label{tab:table2}
\end{table}
Table \ref{tab:table2} reports the mean ASCFE for the group SCAD strategy,
given different combinations of $T$ and $J$. An immediate observation is
that the ASCFE is comparable to that obtained in the low-dimension case in
table \ref{tab:table1}. Given a fixed $J$, we observe that the ASCFE falls
and the share of iterations that accurately captures the relevant forecasts
increases with larger samples. This provides some evidence that the
estimator is consistent in estimating and selecting important
forecasts.
\section{Empirical Applications}
\label{empirical}
We illustrate the use of our nonparametric estimator by considering two
applications of combined forecasting in macroeconomics and finance.
Specifically, we review the results on forecasting inflation in \cite
{ang2007macro} by extending their analysis to include recent data.
Subsequently, we follow \cite{rapach2010out} to examine the predictability
of equity returns with forecast combinations using 13 predictors. In the
latter context, since the effective number of variables for the local linear
estimator is two times that of the original, it may still be large relative
to the effective sample size even if it does not exceed it. We show that
applying the group SCAD strategy in this case yields two benefits over
conventional methods: estimation is more precise in terms of smaller errors,
and important equity return predictors can be identified. Predictor selection is of independent interest in the return
predictability literature.
To facilitate statistical comparison between forecast strategies, we conduct the
\cite{diebold1995comparing} test (DM test) and the 'Reality Check' (RC) test introduced by \cite
{white2000reality}. The latter can be used to compare forecasts generated from
both nested and non-nested models, which is especially relevant in the
high-dimensional context because of variable selection.
\subsection{Forecasting inflation}
\cite{ang2007macro} studied the performance of various strategies in
forecasting inflation in the US from 1985 to 2002. These approaches can be
classified into four broad categories: ARIMA-type time series models,
Phillips curve-implied forecasts, term structure models, and survey-based
measures. Two major results from their investigation are particularly
striking. First, they find that the median of survey forecasts, in
particular, the Livingstone survey and the Survey of Professional
Forecasters (SPF)\footnote{
See \cite{ang2007macro} for a detailed description of the surveys. We do not consider the Livingstone survey in our application because it is a biannual survey.}, consistently outperform models from the other three categories. Furthermore, they show that combining forecasts across the different categories using various parsimonious methods, such as least squares regression and equal weights, do not generally yield more accurate forecasts compared to using only survey information. Our goal is thus twofold: we are interested in assessing whether this claim remains valid given an updated evaluation period, and to ascertain the potential gains from using the nonparametric estimator as
opposed to existing methods in combining inflation forecasts.
To achieve this, we use the basic forecasting models considered by \cite{ang2007macro} from each of the four categories: ARMA(1,1), PC1, TS1, and the SPF for a survey-based forecast\footnote{Naming convention and models follow that of \cite{ang2007macro}}. We consider three CPI-based measures of inflation: CPI for all urban consumers (PUNEW), CPI less housing and shelter (PUXHS), and CPI less
food and energy (PUXX). We generate inflation forecasts for a period of 1981Q3 to 2018Q2, and use the last 3 years for our out-of-sample evaluation. We refer interested readers to the online appendix for a full description of the set-up.
For all measures of inflation, we look at the same combination
techniques as detailed in section \ref{competing}. In addition, we consider individual predictions made by the SPF and the ARMA(1,1) model on quarterly inflation which Ang et al. (2007) regarded as a benchmark. We note that the methods TVGRreg and EQ respectively correspond to the "OLS" and "Mean" methods employed in their paper. Table \ref{tab:table4} reports the ASCFEs, the best methods, and the p-values for one-sided DM and RC tests.
\begin{table}[tbp]
\centering
\caption{Forecast combination results for inflation.}
\begin{threeparttable}
\begin{tabular}{lccc}
\toprule
Estimation & PUNEW & PUXHS & PUXX \\
\midrule
\textit{Adaptive} & & & \\
NPRf & 0.182 & 0.451 & 0.050 \\
BG & 0.457 & 0.679 & 0.117 \\
TVGRregconst & 0.186 & 0.792 & 0.048 \\
TVGRreg & 0.475 & 0.512 & 0.108 \\
TVGRregconstr & 0.446 & 0.791 & 0.098 \\
ARMA(1,1) & 0.206 & 1.897 & 0.061 \\
\textit{Static} & & & \\
GRregconst & 0.186 & 0.929 & 0.048 \\
GRreg & 0.477 & 0.518 & 0.110 \\
GRregconstr & 0.445 & 0.803 & 0.099 \\
EQ & 0.485 & 0.702 & 0.127 \\
SPF & 0.467 & 0.842 & 0.441 \\
\midrule
Best & NPRf & NPRf & TVGRregconst \\
\midrule
& & & \\
DM test (p-value) & & & \\
NPRf & & & \\
$<$ EQ & 0.142 & 0.314 & 0.084 \\
$<$ SPF & 0.063 & 0.231 & 0.014 \\
& & & \\
RC test (p-value) & & & \\
NPRf & & & \\
$<$ EQ & 0.031 & 0.179 & 0.033 \\
$<$ SPF & 0.025 & 0.057 & 0.003 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\small
\item Notes: (1) ASCFEs for inflation forecasting in the top panel. (2) NPRf: nonparametric estimator with data reflection, BG: \cite{bates1969combination} adaptive etimator, GRregconst: OLS regression with intercept, GRreg: OLS regression without intercept, GRregconstr: OLS regression with sum of coefficients constrained to unity, TV: time-varying weights, EQ: equal weights, SPF: Survey of Professional Forecasters, PUNEW: CPI (all urban) inflation, PUXHS: CPI (less housing and shelter) inflation, PUXX: CPI (less food and energy) inflation. (3) DM test: \protect\cite{diebold1995comparing} test. RC test:
'Reality Check' test \citep{white2000reality}. (4) "$x < y$" indicates a test of null hypothesis of equal predictive ability between $x$ and $y$, with the one-sided alternative of superior predictive ability of $x$ over $y$.
\end{tablenotes}
\end{threeparttable}
\label{tab:table4}
\end{table}
Several comments are in order here. First, the nonparametric method, NPRf, achieves the best performance in two out of three inflation indices, but we note that the various OLS methods are close contenders. For PUXX, both TVGRregconst and GRregconst achieve the lowest ASCFE, but NPRf is a close runner-up. Furthermore, there is statistical evidence that NPRf improves over equal weighting and simple SPF-implied forecasts for PUNEW and PUXX. This finding contrasts with \cite{ang2007macro}, who found that even with forecast combinations, the weight on the SPF forecast dominated all other contributions. Here, we find that other types of forecast can improve inflation forecasting above and beyond using just survey estimates.
It has to be emphasized that an additional merit of employing NPRf in this situation is that we do not have to take a strong and often contentious stand on the stationarity of inflation, since given our formulation of the estimator, it accommodates some extent of non-stationarity including locally stationary processes. Regardless, the evidence presented thus far appear to favor the use of forecast combination for inflation forecasting, which is a notable departure from the conclusions arrived at in \cite{ang2007macro}.
\subsection{Predictability of stock returns}
\label{stock}
Many popular macroeconomic and financial variables do not posses
out-of-sample predictive power in forecasting stock returns even if they may exhibit good in-sample performance \citep{bossaerts1999implementing, welch2008comprehensive}. However, models that do successfully deliver
statistically significant and economically relevant out-of-sample
forecasting gains are often approaches that accommodate model uncertainty and parameter instability \citep{rapach2013forecasting}. This is not particularly surprising given that forecasters do not possess \textit{a priori} information on the "best" models. Furthermore, there is strong evidence that many predictive relationships of stock returns are unstable over time \citep{chen2012testing}. Forecast combination fits well into this framework because it helps to attenuate uncertainty in individual forecasts and can accommodate time variation. In fact, \cite{rapach2010out} reported significant out-of-sample predictive gains in forecasting US stock returns compared to the historical average using combination strategies that overlap with those that we have introduced in section \ref{competing}. Here, we investigate whether our nonparameteric estimator can achieve similar favorable out-of-sample performace.
To do so, we use updated data (till 2018) from \cite{welch2008comprehensive} and construct the predictors and excess return in a similar fashion. The dependent variable is the stock return derived from the S\&P 500 index and is defined as $\Delta P_{t+1}=[\log (P_{t+1}+D_{t})-\log(P_{t})]-R_{t}$, where $P_{t}$ is the index value, $D_{t}$ is the dividends paid on the index, and $R_{t}$ is the 3-month Treasury bill rate. We use 13 predictors following Rapach et al. (2010), whose full description can be found on the online appendix. Our sample period is from 1947Q2 to 2018Q3, and we use the last 5 years, 2013Q4 to 2018Q3, for out-of-sample evaluation.
As usual, we compare the performance of static and adaptive estimators as introduced in section \ref{competing}. In addition, we look at the
performance of the prevailing historical average defined as $\bar{\Delta} P_{t+1,T-1}=1/(T-1)\sum_{t=1}^{T-1}\Delta P_{t+1}$, which is essentially a moving average with an expanding window, as this was the benchmark in \cite{rapach2010out}, and is supposedly difficult to beat \citep{welch2008comprehensive}.
Table \ref{tab:table5} presents the ASCFEs scaled by a multiplication of 1000. In general, the adaptive methods seem to perform better than the static OLS-based approaches which is consistent with the notion of parameter instability. However, we observe that the NPRf does not perform favorably in this scenario, as it beats neither the historical average nor the equal weighting scheme. One possible reason for this could be due to the high number of predictors ($13 \times 2$) relative to the sample size. We check this intuition by considering the group SCAD (gSCAD) penalized version of NPRf.
Investigating the value of gSCAD for variable selection is particularly relevant given that the number of predictors in the return predictability literature has been on the rise. Furthermore, predictor selection \textit{per se} is of strong independent interest in asset pricing. To further evaluate our penalization method, we compare it with the partially egalitarian approaches introduced by \cite{diebold2019machine}, which seeks to select important variables in a first-stage Lasso and subsequently shrink their weights to the arithmetic mean with a modified Ridge estimator. Their method capitalizes on the frequently reported empirical finding that simple averages of forecasts tend to outperform sophisticated estimation strategies.
\begin{table}[htbp]
\centering
\caption{Forecast combination results for stock returns.}
\begin{threeparttable}
\begin{tabular}{lclc}
\toprule
Estimation & & \multicolumn{2}{c}{ASCFE} \\
\midrule
\textit{Adaptive} & & & \\
NPRf & & \multicolumn{2}{c}{3.183} \\
BG & & \multicolumn{2}{c}{2.879} \\
TVGRregconst & & \multicolumn{2}{c}{4.992} \\
TVGRreg & & \multicolumn{2}{c}{4.663} \\
TVGRregconstr & & \multicolumn{2}{c}{2.991} \\
ARMA(1,1) & & \multicolumn{2}{c}{2.807} \\
Historical average & & \multicolumn{2}{c}{3.034} \\
\textit{Static} & & & \\
GRregconst & & \multicolumn{2}{c}{12.506} \\
GRreg & & \multicolumn{2}{c}{9.969} \\
GRregconstr & & \multicolumn{2}{c}{3.125} \\
EQ & & \multicolumn{2}{c}{2.888} \\
& & \multicolumn{2}{c}{} \\
\midrule
Penalized estimation & ASCFE & \multicolumn{2}{c}{Selected predictors} \\
\cmidrule{3-4} & & For all periods & Sometimes selected \\
\midrule
gSCAD & 2.650 & \multicolumn{1}{p{10em}}{E/P, SVAR, NTIS, \newline{}TBL, I/K} & \multicolumn{1}{p{10em}}{B/M} \\
peLasso & 3.349 & TBL, I/K & \\
& & & \\
\midrule
Best & gSCAD & \multicolumn{2}{c}{} \\
& & & \\
\midrule
DM test (p-value) & & & \\
gSCAD & & & \\
$<$ Hist. avg. & & \multicolumn{2}{c}{0.065} \\
$<$ EQ & & \multicolumn{2}{c}{0.106} \\
$<$ peLasso & & \multicolumn{2}{c}{0.027} \\
& & & \\
RC test (p-value) & & & \\
gSCAD & & & \\
$<$ Hist. avg. & & \multicolumn{2}{c}{0.129} \\
$<$ EQ & & \multicolumn{2}{c}{0.232} \\
$<$ peLasso & & \multicolumn{2}{c}{0.093} \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\small
\item Notes: (1) ASCFEs multiplied by 1000. (2) NPRf: nonparametric estimator with data reflection, BG: \cite{bates1969combination} adaptive etimator, GRregconst: OLS regression with intercept, GRreg: OLS regression without intercept, GRregconstr: OLS regression with sum of coefficients constrained to unity, TV: time-varying weights, EQ: equal weights, gSCAD: NPRf with group SCAD penalties, peLasso: partially egalitarian Lasso \citep{diebold2019machine}. (3) DM test: \protect\cite{diebold1995comparing} test. RC test:
'Reality Check' test \citep{white2000reality}. (4) "$x < y$" indicates a test of null hypothesis of equal predictive ability between $x$ and $y$, with the one-sided alternative of superior predictive ability of $x$ over $y$. See online appendix for predictor description.
\end{tablenotes}
\end{threeparttable}
\label{tab:table5}
\end{table}
The results in table \ref{tab:table5} show that gSCAD is able to achieve the lowest ASCFE and provide some statistical evidence that gSCAD possesses superior predictive ability over the historical average and the competing regularized estimator, peLasso. Although it appears to perform better than EQ, the statistical evidence is borderline. This result is not entirely unexpected given the short out-of-sample evaluation period and arguably unsatisfactory finite-sample properties of tests for out-of-sample predictive ability \citep{elliott2005optimal, chao2001out}. On another note, it is interesting to see that both gSCAD and peLasso agree on the relevance of the contributions from macro-level variables, the treasury rate (TBL) and the investment-to-capital ratio (I/K), in predicting equity returns for the considered sample period. Although gSCAD suggests that at least 4 other predictors are also important. Capturing the predictors that were excluded by peLasso might have contributed to a better performance.
\section{Concluding Remarks}
In this paper, we have proposed a new approach to forecast combination with time-varying weights by means of a nonparametric local linear estimator with data reflection. The estimator can be adapted to situations where the number of forecasts is large, possibly larger than the sample size, through the use of sparsity and penalization techniques. Theoretically, we have shown that the nonparametric estimator with the group SCAD penalty is consistent in both selecting relevant forecasts and in estimating the weights. Simulations have shown that the nonparametric estimator performs well in situations of parameter instability, and the relative success of our estimator in empirical applications suggests that predictive instability is the norm rather than the exception. Our proposed nonparametric estimator may also be generalized to forecast combinations with volatility and density forecasts, and it may be interesting to explore these avenues in future research given their popularity in applications.
\bibliographystyle{elsarticle-harv}
\bibliography{Chen_Maung}
\processdelayedfloats
\makeatletter
\efloat@restorefloats \makeatother