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.
75,569 characters · 15 sections · 95 citation commands
\thispagestyle{empty}
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.
JEL Classifications: C12, C14, C22
Key words: Cross validation, Forecast combination, High dimension, Local linear estimation, SCAD, Sparsity.
\verb||
\verb||
\setcounter{page}{1}
Multiple forecasts of the same variable are often available to decision makers. As pointed out in the seminal paper by 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. 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 elliott2005optimal, stock2004combination to empirical asset pricing 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 bates1969combination, who developed several adaptive estimation schemes for time-varying weights based on exponential discounting or rolling estimation. In the regression context, diebold1987structural generalize these schemes to select combination weights that minimize the weighted average of forecast errors. deutsch1994combination and 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, 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 robinson1989nonparametric, robinson1991time and further studied by cai2007trending, chen2012testing, kristensen2012non, zhang2012inference, dahlhaus2019towards, 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 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 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 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 tibshirani1996regression, smoothly clipped absolute deviation (SCAD) fan2001variable, group Lasso yuan2006model, adaptive Lasso zou2006adaptive, and the minimax concave penalty (MCP) zhang2010nearly. In a related context, model selection for functional coefficient models under the $i.i.d.$ assumption are considered in wang2009shrinkage, wei2011variable and li2015model, and model selection for high-dimsional linear time series models are studied in kock2015oracle, han2020high and 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 deutsch1994combination and 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 onlineappendix.
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:
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:
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 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 elliott2005optimal or smooth transitions deutsch1994combination with some observed or latent state variable. }. Specifically, we adopt the following framework of the nonparametric time-varying parameter model:
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 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 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:
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 (ref) with respect to $\gamma_t$ yields the local linear estimate of $\beta(t/T)$,
where $e_1 = (1,0)^\top$, $I_{(d+1)}$ is a $(d+1) \times (d+1)$ identity matrix, and
The estimator takes the form of a leave-one-out local linear estimator considered in 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.
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.
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 $.
Assumption A.2 (Moment conditions): The following moment conditions are satisfied: \setlist{nolistsep}
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 \}$.
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$.
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 cai2007trending, orbe2005nonparametric, robinson1989nonparametric, and requires slightly more than four moments of the data. More importantly, unlike cai2007trending and 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 dahlhaus2019towards and 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 cai2007trending for forecast combination with the reflection method.
\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 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.
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]$.
\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
Here, we eliminate the cross product term since $\hat{\beta}_{t} $ only uses information up to time $t$. Subsequently, define the integrated MSCFE as
where $V_{\beta }(\tau )\equiv M^{-1}(\tau )V(\tau )M^{-1}(\tau )$ and label the second term in (ref) 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
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,
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.
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. hardle1985optimal and xia2002asymptotic impose similar moment conditions.
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.
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
where the $f_{t}=\left( f_{t1,}f_{t2},...,f_{tp_{T}}\right) ^{\top }$ 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}$ 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:
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
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
The Lasso-based local linear estimator is estimation consistent (see \thref{fsprop}) but requires very strong assumptions for selection consistency. Instead, following li2015model, we minimize (ref) to obtain preliminary estimates for use in a second-stage penalized optimization with the (group) SCAD penalty proposed by fan2001variable. In particular, we use the initial estimates from $\tilde{B}$ to solve
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
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
and $a=3.7$ as suggested in 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 zou2008one, fan2014strong. Then the local linear estimator for $ \beta _{t}$ with the group SCAD penalty is
where $\hat{\gamma}_{t}^{h}$ is the $t^{th}$ row of $\hat{\Upsilon}^{h}.$
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 (ref) 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 (ref) in matrix notation as
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
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 wei2011variable and 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
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.
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}$.
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)}$.
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)}$.
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)})$.
Step 5. Repeat for all the coefficient vectors until a reasonable tolerance is achieved. We use $1\times10^{-3}$ in our simulations and applications.
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, 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) since we do not expect the preliminary estimator to be selection consistent for our asymptotic results.
Following 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:
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 wang2009shrinkagetuning. Here, BIC is a popular tuning approach for variable selection problems with SCAD penalties. For example, 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 cai2015functional}.
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.
Assumption HD.1 (Moment conditions): \setlist{nolistsep}
Assumption HD.2 (Restricted eigenvalues): Define the set $ S$ as
for some $\delta>0$. Then there exists positive constants $ 0<\rho_1\leq\rho_2<\infty$ with probability approaching one such that,
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 (ref) in Section (ref).
Assumption HD.3 (Rates and tuning parameters): \setlist{nolistsep}
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 li2015model, while HD.2 is a generalization of the restricted eigenvalue conditions in 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}.$
\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 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 }.$
\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.
\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 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), and a direct application of \thref{consis} is possible.
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.
To start off, we consider the following time-varying coefficient model:
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:
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 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 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.
Nonparametric estimation
We consider the local linear estimator with data reflection as established in section (ref). For bandwidth selection, we adopt the CV method introduced in (ref). 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})_{+}$.
Bates and Granger (1969)
A commonly used time-varying combination scheme, originally suggested by 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
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".
Least squares regression
We consider three static combination schemes due to granger1984improved. They are:
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 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.
Equal weights
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".
Table (ref) 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.
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.
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) 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 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), 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 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) 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.
Table (ref) 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). 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.
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 ang2007macro by extending their analysis to include recent data. Subsequently, we follow 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 diebold1995comparing test (DM test) and the 'Reality Check' (RC) test introduced by 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.
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 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 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 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). 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) reports the ASCFEs, the best methods, and the p-values for one-sided DM and RC tests.
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 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 ang2007macro.
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 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 rapach2013forecasting. This is not particularly surprising given that forecasters do not possess a priori information on the "best" models. Furthermore, there is strong evidence that many predictive relationships of stock returns are unstable over time chen2012testing. Forecast combination fits well into this framework because it helps to attenuate uncertainty in individual forecasts and can accommodate time variation. In fact, 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). 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 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). 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 rapach2010out, and is supposedly difficult to beat welch2008comprehensive.
Table (ref) 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 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 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.
The results in table (ref) 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 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.
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.
\processdelayedfloats
\makeatletter \efloat@restorefloats \makeatother