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.
83,731 characters · 16 sections · 56 citation commands
Learning from Forecast Errors: A New Approach to Forecast Combinations
\thispagestyle{empty}
{22pt} \setcounter{page}{1} \setstretch{1.9}
A search for the best forecast combination has been an important on-going research question in economics. CLEMEN1989 pointed out that combining forecasts is \enquote{practical, economical and useful. Many empirical tests have demonstrated the value of composite forecasting. We no longer need to justify that methodology}. However, as demonstrated by DIEBOLD2018, there are still some unresolved issues. Despite the findings based on the theoretical grounds, equal-weighted forecasts have proved surprisingly difficult to beat. Many methodologies that seek for the best forecast combination use equal weights as a benchmark: for instance, DIEBOLD2018 develop \enquote{partially egalitarian Lasso}.
The success of equal weights is partly due to the fact that the forecasters use the same set of public information to make forecasts, hence, they tend to make common mistakes. For example, in the European Central Bank's Survey of Professional forecasters of Euro-area real GDP growth, the forecasters tend to jointly understate or overstate GDP growth. Therefore, we stipulate that the forecast errors include common and idiosyncratic components, which allows the forecast errors to move together due to the common error component. Our paper provides a simple framework to learn from analyzing forecast errors: we separate unique errors from the common errors to improve the accuracy of the combined forecast.
Dating back to GrangerBatesWeights, the well-known expression for the optimal forecast combination weights requires an estimator of inverse covariance (precision) matrix. Graphical models are a powerful tool to estimate precision matrix directly, avoiding the step of obtaining an estimator of covariance matrix to be inverted. Prominent examples of graphical models include Graphical Lasso (GLASSO) and nodewise regression (meinshausen2006). Despite using different strategies for estimating precision matrix, all graphical models assume that the latter is sparse: many entries of precision matrix are zero, which is a necessary condition to consistently estimate inverse covariance. Our paper demonstrates that such assumption contradicts the stylized fact that experts tend to make common mistakes and hence the forecast errors move together through common factors. We show that graphical models fail to recover entries of precision matrix under the factor structure.
This paper overcomes the aforementioned challenge and develops a new precision matrix estimator for the forecast errors under the approximate factor model with unobserved factors. We call our algorithm the Factor Graphical Model. We use a factor model to estimate an idiosyncratic component of the forecast errors, and then apply a Graphical model (Graphical Lasso or nodewise regression) for the estimation of the precision matrix of the idiosyncratic component.
There are a few papers that used graphical models in different contexts to estimate the covariance matrix of the idiosyncratic component when the factors are known and the loadings are assumed to be constant. Brownlees2018JAE estimate a sparse covariance matrix for high-frequency data and construct the realized network for financial data. Brownlees2018EJS develop a power-law partial correlation network based on the Gaussian graphical models. koike2019biased uses the Weighted Graphical Lasso to estimate a sparse covariance matrix of the idiosyncratic component for a factor model with observable factors for high-frequency financial data.
Our paper makes several contributions. First, we allow the forecast errors to be highly correlated due to the common component which is motivated by the stylized fact that the forecasters tend to jointly understate or overstate the predicted series of interest. Second, we develop a high-dimensional precision matrix estimator which combines the benefits of the factor structure and sparsity of the precision matrix of the idiosyncratic component for the forecast combination under the approximate factor model. We prove consistency of forecast combination weights and the Mean Squared Forecast Error (MSFE) estimated using Factor Graphical models. Third, an empirical application to forecasting macroeconomic series in big data environment shows that incorporating the factor structure of the forecast errors into the graphical models improves the performance of a combined forecast over forecast combination using equal weights and graphical models without factors.
The paper is structured as follows: Section 2 reviews Graphical Lasso and nodewise regression. Section 3 studies the approximate factor models for the forecast combination. Section 4 introduces the Factor Graphical Models and discusses the choice of the tuning parameters. Section 5 contains theoretical results and Section 6 validates these results using simulations. Section 7 studies an empirical application for macroeconomic time-series. Section 8 concludes and Section 9 collects the proofs of the theorems.
Notation. For the convenience of the reader, we summarize the notation to be used throughout the paper. Let $\mathcal{S}_p$ denote the set of all $p \times p$ symmetric matrices. For any matrix ${\mathbf C}$, its $(i,j)$-th element is denoted as $c_{ij}$. Given a vector ${\mathbf u}\in \mathbb{R}^d$ and a parameter $a\in \lbrack1,\infty)$, let $\@ifstar{\oldnorm}{\oldnorm*}{{\mathbf u}}_a$ denote $\ell_a$-norm. Given a matrix ${\mathbf U} \in\mathcal{S}_p$, let $\Lambda_{\text{max}}({\mathbf U}) \equiv \Lambda_1({\mathbf U}) \geq \Lambda_2({\mathbf U})\geq \ldots \geq \Lambda_{\text{min}}({\mathbf U}) \equiv \Lambda_p({\mathbf U})$ be the eigenvalues of ${\mathbf U}$. Given a matrix ${\mathbf U} \in \mathbb{R}^{p\times p}$ and parameters $a,b\in \lbrack1,\infty)$, let ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {\mathbf U} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{a,b}\equiv \max_{\@ifstar{\oldnorm}{\oldnorm*}{{\mathbf y}}_a=1}\@ifstar{\oldnorm}{\oldnorm*}{{\mathbf U}{\mathbf y}}_{b}$ denote the induced matrix-operator norm. The special cases are ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {\mathbf U} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1\equiv \max_{1\leq j\leq p}\sum_{i=1}^{p}\@ifstar{\oldabs}{\oldabs*}{u_{i,j}}$ for the $\ell_1/\ell_1$-operator norm; the operator norm ($\ell_2$-matrix norm) ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {\mathbf U} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{2}^{2}\equiv\Lambda_{\text{max}}({\mathbf U}{\mathbf U}')$ is equal to the maximal singular value of ${\mathbf U}$. Finally, $\@ifstar{\oldnorm}{\oldnorm*}{{\mathbf U}}_{\infty}\equiv\max_{i,j}\@ifstar{\oldabs}{\oldabs*}{u_{i,j}}$ denotes the element-wise maximum.
This section briefly reviews a class of models, called graphical models, that search for the estimator of the precision matrix. In graphical models, each vertex represents a random variable, and the graph visualizes the joint distribution of the entire set of random variables. Sparse graphs have a relatively small number of edges.
Suppose we have $p$ competing forecasts of the univariate series $y_t$, $t=1,\ldots,T$. Let ${\mathbf e}_t=(e_{1t},\ldots,e_{pt})' \sim \mathcal{N} (\mathbf{0}, {\bm \Sigma})$ be a $p \times 1$ vector of forecast errors. Assume they follow a Gaussian distribution. The precision matrix ${\bm \Sigma}^{-1}\equiv {\bm \Theta}$ contains information about partial covariances between the variables. For instance, if $\theta_{ij}$, which is the $ij$-th element of the precision matrix, is zero, then the variables $i$ and $j$ are conditionally independent, given the other variables.
Let ${\mathbf W}$ be the estimate of ${\bm \Sigma}$. Given a sample $\{{\mathbf e}_t\}_{t=1}^{T}$, let ${\mathbf S} = (1/T)\sum_{t=1}^{T}({\mathbf e}_t)({\mathbf e}_t)'$ denote the sample covariance matrix, which can be used as a choice for ${\mathbf W}$. Also, let $\widehat{{\mathbf D}}^2\equiv \textup{diag}({\mathbf W})$. We can write down the Gaussian log-likelihood (up to constants) $l({\bm \Theta})=\log\det({\bm \Theta})-\text{trace}({\mathbf W}{\bm \Theta})$. When ${\mathbf W}={\mathbf S}$, the maximum likelihood estimator of ${\bm \Theta}$ is $\widehat{{\bm \Theta}}={\mathbf S}^{-1}$.
In the high-dimensional settings it is necessary to regularize the precision matrix, which means that some edges will be zero. In the following subsections we discuss two most widely used techniques to estimate sparse high-dimensional precision matrices.
The first approach to induce sparsity in the estimation of precision matrix is to add penalty to the maximum likelihood and use the connection between the precision matrix and regression coefficients to maximize the following weighted penalized log-likelihood (Sara2018):
over positive definite symmetric matrices, where $\lambda\geq0$ is a penalty parameter. The subscript $\lambda$ in $\widehat{{\bm \Theta}}_{\lambda}$ means that the solution of the optimization problem in (ref) will depend upon the choice of the tuning parameter. More details on the latter are provided in Subsection 4.1 that describes how to choose the shrinkage intensity in practice. In order to simplify notation, we will omit the subscript.
One of the most popular and fast algorithms to solve the optimization problem in (ref) is called the Graphical Lasso (GLASSO), which was introduced by GLASSO. Define the following partitions of ${\mathbf W}$, ${\mathbf S}$ and ${\bm \Theta}$:
Let ${\bm \beta}\equiv -{\bm \theta}_{12}/\theta_{22}$. The idea of GLASSO is to set ${\mathbf W}= {\mathbf S}+\lambda{\mathbf I}$ in (ref) and combine the gradient of (ref) with the formula for partitioned inverses to obtain the following $\ell_1$-regularized quadratic program
As shown by GLASSO, (ref) can be viewed as a LASSO regression, where the LASSO estimates are functions of the inner products of ${\mathbf W}_{11}$ and $s_{12}$. Hence, (ref) is equivalent to $p$ coupled LASSO problems. Once we obtain $\widehat{{\bm \beta}}$, we can estimate the entries ${\bm \Theta}$ using the formula for partitioned inverses. GLASSO procedure is summarized in Algorithm (ref).
As was shown in GLASSO, the estimator produced by Algorithm (ref) is guaranteed to be positive definite. Furthermore, Sara2018 showed that Algorithm (ref) is guaranteed to converge and produces consistent estimator of precision matrix under certain sparsity conditions.
An alternative approach to induce sparsity in the estimation of precision matrix in equation (ref) is to solve for $\widehat{{\bm \Theta}}$ one column at a time via linear regressions, replacing population moments by their sample counterparts ${\mathbf S}$. When we repeat this procedure for each variable $j=1,\ldots,p$, we will estimate the elements of $\widehat{{\bm \Theta}}$ column by column using $\{{\mathbf e}_t\}_{t=1}^{T}$ via $p$ linear regressions. meinshausen2006 use this approach (which we will refer to as MB) to incorporate sparsity into the estimation of the precision matrix. Instead of running $p$ coupled LASSO problems as in GLASSO, they fit $p$ separate LASSO regressions using each variable (node) as the response and the others as predictors to estimate $\widehat{{\bm \Theta}}$. This method is known as the \enquote{nodewise} regression and it is reviewed below based on Buhlmann2014 and Caner2019.
Let ${\mathbf e}_j$ be a $T \times 1$ vector of observations for the $j$-th regressor, the remaining covariates are collected in a $T \times p$ matrix ${\mathbf E}_{-j}$. For each $j=1,\ldots,p$ we run the following Lasso regressions:
where $\widehat{{\bm \gamma}}_j=\{\widehat{\gamma}_{j,k}; j=1,\ldots,p, k\neq j\}$ is a $(p-1)\times 1$ vector of the estimated regression coefficients that will be used to construct the estimate of the precision matrix, $\widehat{{\bm \Theta}}$. Define
For $j=1,\ldots,p$, define
and write
The approximate inverse is defined as
Similarly to GLASSO, the subscript $\lambda_j$ in $\widehat{{\bm \Theta}}_{\lambda_j}$ means that the estimated ${\bm \Theta}$ will depend upon the choice of the tuning parameter: more details are provided in Subsection 4.1 which discusses how to choose shrinkage intensity in practice. The subscript is omitted to simplify the notation. The procedure to estimate the precision matrix using nodewise regression is summarized in Algorithm (ref).
One of the caveats to keep in mind when using the MB method is that the estimator in (ref) is not self-adjoint. Caner2019 show (see their Lemma A.1) that $\widehat{{\bm \Theta}}$ in (ref) is positive definite with high probability, however, it could still occur that $\widehat{{\bm \Theta}}$ is not positive definite in finite samples. In such cases we use the matrix symmetrization procedure as in fan2018elliptical and then use eigenvalue cleaning as in Callot2017 and Hautsch2012.
The approximate factor models for the forecasts were first considered by StockChan1999. They modeled a panel of ex-ante forecasts of a single time-series as a dynamic factor model and found out that the combined forecasts improved on individual ones when all forecasts have the same information set (up to difference in lags). This result emphasizes the benefit of forecast combination even when the individual forecasts are not based on different information and, therefore, do not broaden the information set used by any one forecaster.
In this paper, we are interested in finding the combination of forecasts which yields the best out-of-sample performance in terms of the mean-squared forecast error. We claim that the forecasters use the same set of public information to make forecasts and hence they tend to make common mistakes. Figure (ref) illustrates this statement: it shows quarterly forecasts of Euro-area real GDP growth produced by the European Central Bank's Survey of Professional Forecasters from 1999Q3 to 2019Q3. As described in DIEBOLD2018, forecasts are solicited for one year ahead of the latest available outcome: e.g., the 2007Q1 survey asked the respondents to forecast the GDP growth over 2006Q3-2007Q3. As evidenced from Figure (ref), forecasters tend to jointly understate or overstate GDP growth, meaning that their forecast errors include common and idiosyncratic parts. Therefore, we can model the tendency of the forecast errors to move together via factor decomposition.
Recall that we have $p$ competing forecasts of the univariate series $y_t$, $t=1,\ldots,T$ and ${\mathbf e}_t=(e_{1t},\ldots,e_{pt})' \sim \mathcal{N} (\mathbf{0}, {\bm \Sigma})$ is a $p \times 1$ vector of forecast errors. Assume that the generating process for the forecast errors follows a $q$-factor model:
where ${\mathbf f}_t=(f_{1t},\ldots, f_{qt})'$ are the common factors of the forecast errors for $p$ models, ${\mathbf B}$ is a $p \times q$ matrix of factor loadings, and ${\bm \varepsilon}_t$ is the idiosyncratic component that cannot be explained by the common factors. Unobservable factors, ${\mathbf f}_{t}$, and loadings, ${\mathbf B}$, are usually estimated by the principal component analysis (PCA), studied in Connor1988,Bai2003,Bai2002,Stock2002. Strict factor structure assumes that the idiosyncratic forecast error terms, ${\bm \varepsilon}_{t}$, are uncorrelated with each other, whereas approximate factor structure allows correlation of the idiosyncratic components (Chamberlain).
We use the following notations: $\operatorname*{\mathbb{E}}{{\bm \varepsilon}_t{\bm \varepsilon}'_t}={\bm \Sigma}_{\varepsilon}$, $\operatorname*{\mathbb{E}}{{\mathbf f}_t{\mathbf f}'_t}={\bm \Sigma}_{f}$, $\operatorname*{\mathbb{E}}{{\mathbf e}_t{\mathbf e}'_t}={\bm \Sigma}={\mathbf B}{\bm \Sigma}_{f}{\mathbf B}'+ {\bm \Sigma}_{\varepsilon}$, and $\operatorname*{\mathbb{E}}{{\bm \varepsilon}_t|{\mathbf f}_{t}}=0$. Let ${\bm \Theta}={\bm \Sigma}^{-1}$, ${\bm \Theta}_{\varepsilon}={\bm \Sigma}_{\varepsilon}^{-1}$ and ${\bm \Theta}_{f}={\bm \Sigma}_{f}^{-1}$ be the precision matrices of forecast errors, idiosyncratic and common components respectively. The objective function to recover factors and loadings from (ref) is:
where (ref) is the assumption necessary for the unique identification of factors. Fixing the value of ${\mathbf B}$, we can project forecast errors ${\mathbf e}_t$ into the space spanned by ${\mathbf B}$: ${\mathbf f}_t=({\mathbf B}'{\mathbf B})^{-1}{\mathbf B}'{\mathbf e}_t={\mathbf B}'{\mathbf e}_t$. When combined with (ref), this yields a concentrated objective function for ${\mathbf B}$:
It is well-known (see Stock2002 among others) that $\widehat{{\mathbf B}}$ estimated from the first $q$ eigenvectors of $\frac{1}{T}\sum_{t=1}^{T}{\mathbf e}_t{\mathbf e}_t'$ is the solution to (ref). Given a sample of the estimated residuals $\{\widehat{{\bm \varepsilon}}_t={\mathbf e}_t-\widehat{{\mathbf B}}\widehat{{\mathbf f}_t}\}_{t=1}^{T}$ and the estimated factors $\{\widehat{{\mathbf f}}_t\}_{t=1}^{T}$, let $\widehat{{\bm \Sigma}}_{\varepsilon} = (1/T)\sum_{t=1}^{T}\widehat{{\bm \varepsilon}}_t\widehat{{\bm \varepsilon}}_t'$ and $\widehat{{\bm \Sigma}}_{f}=(1/T)\sum_{t=1}^{T}\widehat{{\mathbf f}}_t\widehat{{\mathbf f}}_t'$ be the sample counterparts of the covariance matrices.
Moving forward to the forecast combination exercise, suppose we have $p$ competing forecasts, $\widehat{{\mathbf y}}_{t}=(\hat{y}_{1,t},\ldots,\hat{y}_{p,t})'$, of the variable $y_t$, $t=1,\ldots,T$. The forecast combination is defined as follows:
where ${\mathbf w}$ is a $p \times 1$ vector of weights. Define a measure of risk $\text{MSFE}({\mathbf w}, {\bm \Sigma})={\mathbf w}'{\bm \Sigma}{\mathbf w}$. As shown in GrangerBatesWeights, the optimal forecast combination minimizes the variance of the combined forecast error:
where ${\bm \iota}_p$ is a $p\times 1$ vector of ones. The solution to (ref) yields a $p\times 1$ vector of the optimal forecast combination weights:
If the true precision matrix is known, the equation (ref) guarantees to yield the optimal forecast combination. In reality, one has to estimate ${\bm \Theta}$. Hence, the out-of-sample performance of the combined forecast is affected by the estimation error. As pointed out by smith2009simple, when the estimation uncertainty of the weights is taken into account, there is no guarantee that the \enquote{optimal} forecast combination will be better than the equal weights or even improve the individual forecasts. Define $a = {\bm \iota}'_{p}{\bm \Theta}{\bm \iota}_p/p$, and $\widehat{a} = {\bm \iota}'_{p}\widehat{{\bm \Theta}}{\bm \iota}_p/p$. We can write
and
Therefore, in order to control the estimation uncertainty in the MSFE and combination weights, one needs to obtain a consistent estimator of the precision matrix ${\bm \Theta}$. More details are discussed in Subsection 5.2 and Theorems (ref) and (ref).
Since our interest is in constructing weights for the forecast combination, our goal is to estimate a precision matrix of the forecast errors. However, as pointed out by koike2019biased, when common factors are present across the forecast errors, the precision matrix cannot be sparse because all pairs of the forecast errors are partially correlated given other forecast errors through the common factors. To illustrate this point, we generated forecast errors that follow (ref) with $q=2$ and $\varepsilon_t\sim \mathcal{N}(\mathbf{0},{\bm \Sigma}_{\varepsilon})$, where $\sigma_{\varepsilon,ij}=0.4^{\@ifstar{\oldabs}{\oldabs*}{i-j}}$ is the $i,j$-th element of ${\bm \Sigma}_{\varepsilon}$. The vector of factors ${\mathbf f}_{t}$ is drawn from $\mathcal{N}(\mathbf{0},{\mathbf I}_q/10)$, and the entries of the matrix of factor loadings for forecast error $j=1,\ldots,p$, ${\mathbf b}_{j}$, are drawn from $\mathcal{N}(\mathbf{0},{\mathbf I}_q/100)$. The full loading matrix is given by ${\mathbf B}=({\mathbf b}_{1},\ldots,{\mathbf b}_{p})'$. Let $\widehat{q}$ denote the number of factors estimated by the PCA. We set $(T,p)=(1000,50)$ and plot the heatmap and histogram of population partial correlations of forecast errors ${\mathbf e}_{t}$, which are the entries of a precision matrix, in (ref). We now examine the performance of graphical models for estimating partial correlations under the factor structure. (ref) shows the partial correlations estimated by GLASSO that does not take into account factors: due to strict sparsity imposed by graphical models almost all partial correlations are shrunk to zero which degenerates the histogram in (ref). This means that strong sparsity assumption on ${\bm \Theta}$ imposed by classical graphical models (such as GLASSO and nodewise regression from Algorithms (ref)-(ref)) is not realistic under the factor structure.
In order to avoid the aforementioned problem, instead of imposing sparsity assumption on the precision of forecast errors, ${\bm \Theta}$, we require sparsity of the precision matrix of the idiosyncratic errors, ${\bm \Theta}_{\varepsilon}$. The latter is obtained using the estimated residuals after removing the co-movements induced by the factors (see Brownlees2018EJS,Brownlees2018JAE,koike2019biased). Naturally, once we condition on the common components, it is sensible to assume that many remaining partial correlations of ${\bm \varepsilon}_{t}$ will be negligible and thus ${\bm \Theta}_{\varepsilon}$ is sparse.
We use the weighted Graphical Lasso and nodewise regression as shrinkage techniques to estimate the precision matrix of residuals. Once the precision of the low-rank component is obtained, we use the Sherman-Morrison-Woodbury formula to estimate the precision of forecast errors:
To obtain $\widehat{{\bm \Theta}}_{f}=\widehat{{\bm \Sigma}}_{f}^{-1}$, we use $\widehat{{\bm \Sigma}}_{f}=\frac{1}{T}\sum_{t=1}^{T}\widehat{{\mathbf f}}_{t}\widehat{{\mathbf f}}_{t}^{'}$. To get $\widehat{{\bm \Theta}}_{\varepsilon}$, we develop two approaches: the first uses the weighted GLASSO Algorithm (ref), with the initial estimate of the covariance matrix of the idiosyncratic errors calculated as $\widehat{{\bm \Sigma}}_{\varepsilon}=\frac{1}{T}\sum_{t=1}^{T}\widehat{{\bm \varepsilon}}_{t}\widehat{{\bm \varepsilon}}_{t}^{'}$, where $\widehat{{\bm \varepsilon}}_{t}={\mathbf e}_t-\widehat{{\mathbf B}}\widehat{{\mathbf f}}_{t}$. The second uses nodewise regression and applies Algorithm (ref) to $\widehat{{\bm \varepsilon}}_t$. Once we estimate $\widehat{{\bm \Theta}}_{f}$ and $\widehat{{\bm \Theta}}_{\varepsilon}$, we can get $\widehat{{\bm \Theta}}$ using a sample analogue of (ref). We call the proposed procedures Factor Graphical Lasso and Factor nodewise regression and summarize them in Algorithm (ref) and Algorithm (ref) respectively.
Note that Algorithms (ref) and (ref) involve the tuning parameters $\lambda$ and $\lambda_j$, the procedure on how to choose the shrinkage intensity coefficients is described in more detail in Subsection 4.1 that describes how to choose the shrinkage intensity in practice, and Section 5 that establishes sparsity requirements that guarantee convergence of (ref), (ref), (ref), and (ref).
We can use $\widehat{{\bm \Theta}}$ to estimate the forecast combination weights $\widehat{{\mathbf w}}$
where $\widehat{{\bm \Theta}}$ is obtained from Algorithm (ref) or Algorithm (ref). Let us now revisit the motivating example at the beginning of this section: Figures (ref)-(ref) plot the heatmaps and the estimated partial correlations when precision matrix is computed using Factor GLASSO in Algorithm (ref) with $\widehat{q}\in\{1,2,3\}$ statistical factors. The heatmaps and histograms closely resemble population counterparts in Figure (ref), and the result is not very sensitive to over- or under-estimating the number of factors $\widehat{q}$. This demonstrates that using a combination of classical graphical models and factor structure via Factor Graphical Models in Algorithms (ref)-(ref) improves upon the performance of classical graphical models: our approach allows to extract the benefits of modeling common movements in forecast errors, captured by a factor model, and the benefits of using many competing forecasting models that give rise to a high-dimensional precision matrix, captured by a graphical model.
Algorithms (ref)-(ref) require the tuning parameters $\lambda$ (from Algorithm (ref)) and $\lambda_j$ (from Algorithm (ref)) respectively. We now comment on the choices for both tuning parameters.
To motivate the choice of the tuning parameter for GLASSO and Factor GLASSO, we first briefly discuss some of the existing options to motivate our choice of $\lambda$ in (ref) in simulations and the empirical application. Usually $\lambda$ is selected from a grid of values $F_{\lambda}=(\lambda_{\text{min}},\ldots,\lambda_{\text{max}})$ which minimizes the score measuring the goodness-of-fit. Some popular examples include multifold cross-validation (CV), Stability Approach to Regularization Selection (STARS, STARS), and the Extended Bayesian Information Criteria (EBIC, EBIC). Since we are interested in estimating a sparse high-dimensional precision matrix, we need to choose a method for selecting the tuning parameter which is consistent in high-dimensions. MeinshausenCV suggest that CV performs poorly for high-dimensional data, it overfits (STARS), and it does not consistently select models. Zhu_STARS pointed out that the STARS is not computationally efficient. It is consistent under certain conditions, but suffers from the problem of overselection in estimating Gaussian graphical models. In contrast, EBIC is computationally efficient and is considered to be the state-of-the-art technique for choosing the tuning parameter for the undirected graphs. The score measuring the goodness of fit for EBIC can be written as:
where $\eta\in \lbrack0,1\rbrack$, ${\bm \Theta}_{\varepsilon,\lambda}$ is the precision matrix estimated for the tuning parameter $\lambda\in F_{\lambda}$, and the log-likelihood is $l({\bm \Theta}_{\varepsilon,\lambda})=\log\det({\bm \Theta}_{\varepsilon,\lambda})-\text{trace}({\mathbf W}_{\varepsilon}{\bm \Theta}_{\varepsilon})$. For the estimation of graphical models, the degrees of freedom are usually defined as the number of unique non-zero elements in the estimated precision matrix, $\text{df}({\bm \Theta}_{\varepsilon,\lambda})=\sum_{i\leq j} I_{{\bm \Theta}_{\varepsilon,\lambda,i,j}\neq 0}$. Chen_EBIC showed that when $\eta=1$, EBIC is consistent as long as the dimension $p$ does not grow exponentially with the sample size $T$. Hence, in our simulations and the empirical exercise we use EBIC with $\eta=1$ for GLASSO and Factor GLASSO in Algorithms (ref) and (ref).
For Algorithms (ref) and (ref), we follow Caner2019 to choose $\lambda_j$ in (ref) by minimizing the generalized information criterion (GIC). Let $\@ifstar{\oldabs}{\oldabs*}{\widehat{S}_j(\lambda_j)}$ denote the estimated number of nonzero parameters in the vector $\widehat{{\bm \gamma}}_{\varepsilon,j}$:
As pointed out by Caner2019, the GIC selects the true model with probability approaching one both when $p>T$ and when $p\leq T$.
We first introduce some terminology and notations. Let $A\in \mathcal{S}_p$. Define the following set for $j=1,\ldots,p$:
where $d_j(A)$ is the number of edges adjacent to the vertex $j$ (i.e., the degree of vertex $j$), and $d(A)$ measures the maximum vertex degree. Define $S(A)\equiv \bigcup_{j=1}^{p}D_j(A)$ to be the overall off-diagonal sparsity pattern, and $s(A)\equiv \sum_{j=1}^{p}d_j(A)$ is the overall number of edges contained in the graph. Note that $\text{card}(S(A)) \leq s(A)$: when $s(A)=p(p-1)/2$ this would give a fully connected graph.
For the nodewise regression in (ref), denote $D_j\equiv \{ k; \gamma_{j,k}\neq 0 \}$ to be the active set for row ${\bm \gamma}_{j}$, and let $d_j\equiv \@ifstar{\oldabs}{\oldabs*}{D_j}$. Define $\bar{d}\equiv \max_{1\leq j\leq p}d_j$.
We now list the assumptions on the model (ref):
We also impose strong mixing condition. Let $\mathcal{F}_{-\infty}^{0}$ and $\mathcal{F}_{T}^{\infty}$ denote the $\sigma$-algebras that are generated by $\{({\mathbf f}_t,{\bm \varepsilon}_{t}):t\leq 0\}$ and $\{({\mathbf f}_t,{\bm \varepsilon}_{t}):t\geq T\}$ respectively. Define the mixing coefficient
Assumption (ref) divides the eigenvalues into the diverging and bounded ones. This assumption is satisfied by the factor model with pervasive factors, which is stated in Assumption (ref). We say that a factor is pervasive in the sense that it has non-negligible effect on a non-vanishing proportion of individual time-series. Assumptions (ref)-(ref) are crucial for estimating a high-dimensional factor model: they ensure that the space spanned by the principal components in the population level ${\bm \Sigma}$ is close to the space spanned by the columns of the factor loading matrix ${\mathbf B}$. Assumption (ref) is a technical condition which is needed to consistently estimate the factors and loadings.
Let ${\bm \Sigma}={\bm \Gamma}{\bm \Lambda}{\bm \Gamma}^{'}$, where ${\bm \Sigma}$ is the covariance matrix of returns that follow factor structure described in equation (ref). Define $\widehat{{\bm \Sigma}}, \widehat{{\bm \Lambda}}_q,\widehat{{\bm \Gamma}}_q$ to be the estimators of ${\bm \Sigma},{\bm \Lambda},{\bm \Gamma}$. We further let $\widehat{{\bm \Lambda}}_q=\text{diag}(\hat{\lambda}_1,\ldots,\hat{\lambda}_q)$ and $\widehat{{\bm \Gamma}}_q=(\hat{v}_1,\ldots,\hat{v}_q)$ to be constructed by the first $q$ leading empirical eigenvalues and the corresponding eigenvectors of $\widehat{{\bm \Sigma}}$ and $\widehat{{\mathbf B}}\widehat{{\mathbf B}}'=\widehat{{\bm \Gamma}}_q\widehat{{\bm \Lambda}}_q\widehat{{\bm \Gamma}}_{q}^{'}$. Similarly to fan2018elliptical, we require the following bounds on the componentwise maximums of the estimators:
Assumptions (ref)-(ref) are needed in order to ensure that the first $q$ principal components are approximately the same as the columns of the factor loadings. The estimator $\widehat{{\bm \Sigma}}$ can be thought of as any “pilot" estimator that satisfies (ref). For sub-Gaussian distributions, sample covariance matrix, its eigenvectors and eigenvalues satisfy (ref)-(ref).
In addition, the following structural assumptions on the model are imposed:
To study the properties of the combination weights in (ref) and MSFE, we first need to establish the convergence properties of precision matrix produced by Algorithms (ref)-(ref). Let $\omega_{T}\equiv \sqrt{\log p/T} +1/\sqrt{p}$. Also, let $s({\bm \Theta}_{\varepsilon})=\mathcal{O}_P(s_T)$ for some sequence $s_T\in(0,\infty)$ and $d({\bm \Theta}_{\varepsilon})=\mathcal{O}_P(d_T)$ for some sequence $d_T\in(0,\infty)$. The deterministic sequences $s_T$ and $d_T$ will control the sparsity ${\bm \Theta}_{\varepsilon}$ for Factor GLASSO. Note that $d_T$ can be smaller than or equal to $s_T$. The reason why we distinguish between these two sequences is to juxtapose it with the sparsity conditions for the Factor MB, where we will only use the analogue of $d_T$ which was defined as $\bar{d}$ at the beginning of this section.
Let $\varrho_{1T}$ be a sequence of positive-valued random variables such that $\varrho_{1T}^{-1}\omega_{T}\xrightarrow{p}0$ and $\varrho_{1T}d_Ts_T\xrightarrow{\text{p}}0$, with $\lambda \asymp \omega_{T}$ (where $\lambda$ is the tuning parameter for the Factor GLASSO in (ref)). seregina2020FGL show that under the Assumptions (ref)-(ref), (ref)-(ref) and (ref), ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \widehat{{\bm \Theta}}-{\bm \Theta} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{1}=\mathcal{O}_P(\varrho_{1T}d_Ts_T)$ for Factor GLASSO. Furthermore, let $\varrho_{2T}$ be a sequence of positive-valued random variables such that $\varrho_{2T}^{-1}\omega_{T}\xrightarrow{p}0$ and $\varrho_{2T}\bar{d}^{2}\xrightarrow{\text{p}}0$, with $\lambda_j \asymp \omega_{T}$ (where $\lambda_j$ is the tuning parameter for Factor nodewise regression in (ref)). seregina2020sparse shows that under the Assumptions (ref)-(ref), (ref)-(ref), and (ref), we have ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \widehat{{\bm \Theta}}-{\bm \Theta} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1=\mathcal{O}_P(\varrho_{2T}\bar{d}^{2})$. It is interesting to compare the rates for precision matrix obtained by two factor graphical models: if $d_T=s_T$, the rates are similar, whereas if $d_T<s_T$ Factor MB is expected to converge faster. In fact, in high dimensions when $p>T$ and $\omega_{T} \simeq \sqrt{\log p/T}$, Factor MB achieves the minimax rate for this problem (see cai2016convergence for the rate expression).
Having established the convergence rates for precision matrix, we now study the properties of the combination weights and MSFE.
Proofs of Theorems (ref)-(ref) can be found in Section 9. Note that the rates of convergence for MSFE and precision matrix ${\bm \Theta}$ are the same and both are faster than the combination weight rates in Theorem (ref). In contrast to classical graphical models in Algorithms (ref)-(ref), the convergence properties of which were examined by Sara2018 among others, the rates in Theorems (ref)-(ref) depend on the sparsity of ${\bm \Theta}_{\varepsilon}$ rather than of ${\bm \Theta}$. This means that instead of assuming that many partial correlations of forecast errors ${\mathbf e}_{t}$ are negligible, which is not realistic under the factor structure, we impose a milder restriction requiring many partial correlations of ${\bm \varepsilon}_{t}$ to be negligible once the common components have been taken into account. Similarly to the comparison of precision matrix ${\bm \Theta}$ obtained by two graphical models, if $d_T<s_T$ Factor MB is expected to converge faster for combination weights and MSFE. In our simulations the rates of Factor Graphical models are comparable, whereas an empirical application shows that for most macroeconomic series that we studied Factor GLASSO outperforms Factor MB. This suggests that for macroeconomic forecasting using weighted penalized log-likelihood and running $p$ coupled LASSO problems for estimating precision matrix is preferable to fitting $p$ separate LASSO regressions using each variable as the response and the others as predictors.
We divide the simulation results into two subsections. In the first subsection we study the consistency of the Factor GLASSO and Factor MB for estimating precision matrix and the combination weights. In the second subsection we evaluate the out-of-sample forecasting performance of combined forecasts based on the Factor Graphical models from Algorithms (ref)-(ref) in terms of the mean-squared forecast error. We compare the performance of forecast combinations based on the factor models with equal-weighted (EW) forecast combination, forecast combinations using GLASSO and nodewise regression from Algorithms (ref)-(ref). Similarly to the literature on graphical models, all exercises use 100 Monte Carlo simulations.
We consider sparse Gaussian graphical models which may be fully specified by a precision matrix ${\bm \Theta}_0$. Therefore, the random sample is distributed as ${\mathbf e}_t=(e_{1t},\ldots,e_{pt})' \sim \mathcal{N}(0,{\bm \Sigma}_{0})$, where ${\bm \Theta}_0=({\bm \Sigma}_{0})^{-1}$ for $t=1,\ldots, T, \ j=1,\ldots, p$. Let $\widehat{{\bm \Theta}}$ be the precision matrix estimator. We show consistency of the Factor GLASSO (Algorithm (ref)) and Factor MB (Algorithm (ref)), in (i) the operator norm, ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \widehat{{\bm \Theta}}-{\bm \Theta}_{0} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{2}$, (ii) $\ell_{1}$/$\ell_{1}$-matrix norm, ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \widehat{{\bm \Theta}}-{\bm \Theta}_{0} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{1}$, and (iii) in $\ell_1$-vector norm for the combination weights, $\@ifstar{\oldnorm}{\oldnorm*}{\widehat{{\mathbf w}}-{\mathbf w}}_1$, where ${\mathbf w}$ is given by (ref).
The forecast errors are assumed to have the following structure:
where ${\mathbf e}_{t}$ is a $p \times 1$ vector of forecast errors following $\mathcal{N}(\bm{0},{\bm \Sigma})$, ${\mathbf f}_{t}$ is a $q \times 1$ vector of factors, ${\mathbf B}$ is a $p\times q$ matrix of factor loadings, $\phi_f$ is an autoregressive parameter in the factors which is a scalar for simplicity, ${\bm \zeta}_t$ is a $q \times 1$ random vector with each component independently following $\mathcal{N}(0,\sigma^{2}_{\zeta})$, ${\bm \varepsilon}_t$ is a $p \times 1$ random vector following $\mathcal{N}(0,{\bm \Sigma}_{\varepsilon})$, with sparse ${\bm \Theta}_{\varepsilon}$ that has a random graph structure described below. To create ${\mathbf B}$ in (ref) we take the first $q$ columns of an upper triangular matrix from a Cholesky decomposition of the $p \times p$ Toeplitz matrix parameterized by $\rho$: that is, ${\mathbf B} = (b)_{ij}$, where $(b)_{ij}=\rho^{\@ifstar{\oldabs}{\oldabs*}{i-j}}$, $i,j\in \{1,\ldots,p\}$. We set $\rho = 0.2$, $\phi_f = 0.2$ and $\sigma^{2}_{\zeta} = 1$. The specification in (ref) leads to the low-rank plus sparse decomposition of the covariance matrix:
When ${\bm \Sigma}_{\varepsilon}$ has a sparse inverse ${\bm \Theta}_{\varepsilon}$, it leads to the low-rank plus sparse decomposition of the precision matrix ${\bm \Theta}$, such that ${\bm \Theta}$ can be expressed as a function of the low-rank ${\bm \Theta}_{f}$ plus sparse ${\bm \Theta}_{\varepsilon}$.
We consider the following setup: let $p = T^{\delta}$, $\delta = 0.85$, $q = 2(\log(T))^{0.5}$ and $T = \lbrack 2^{\kappa} \rbrack, \ \text{for} \ \kappa=7,7.5,8,\ldots,9.5$. Our setup allows the number of individual forecasts, $p$, and the number of common factors in the forecast errors, $q$, to increase with the sample size, $T$.
A sparse precision matrix of the idiosyncratic components ${\bm \Theta}_{\varepsilon}$ is constructed as follows: we first generate the adjacency matrix using a random graph structure. Define a $p \times p$ adjacency matrix ${\mathbf A}_{\varepsilon}$ which represents the structure of the graph:
where $a_{\varepsilon,ij}$ denotes the $i,j$-th element of the adjacency matrix ${\mathbf A}_{\varepsilon}$. We set $a_{\varepsilon,ij} = a_{\varepsilon,ji}=1, \ \text{for} \ i\neq j$ with probability $\pi$, and $0$ otherwise. Such structure results in $s_T = p(p-1)\pi/2$ edges in the graph. To control sparsity, we set $\pi = 1/(pT^{0.8})$, which makes $s_T = \mathcal{O}(T^{0.05})$. The adjacency matrix has all diagonal elements equal to zero. Hence, to obtain a positive definite precision matrix we apply the procedure described in HUGE: using their notation, ${\bm \Theta}_{\varepsilon}={\mathbf A}_{\varepsilon}\cdot v+{\mathbf I}(\@ifstar{\oldabs}{\oldabs*}{\tau}+0.1+u)$, where $u>0$ is a positive number added to the diagonal of the precision matrix to control the magnitude of partial correlations, $v$ controls the magnitude of partial correlations with $u$, and $\tau$ is the smallest eigenvalue of ${\mathbf A}_{\varepsilon}\cdot v$. In our simulations we use $u=0.1$ and $v=0.3$.
Figures (ref)-(ref) show the averaged (over Monte Carlo simulations) errors of the estimators of the precision matrix ${\bm \Theta}$ and the optimal combination weight versus the sample size $T$ in the logarithmic scale (base 2). The estimate of the precision matrix of the EW forecast combination is obtained using the fact that diagonal covariance and precision matrices imply equal weights. To determine the values of the diagonal elements we use the shrinkage intensity coefficient calculated as the average of the eigenvalues of the sample covariance matrix of the forecast errors (see Ledoit2004). As evidenced by Figures (ref)-(ref), Factor GLASSO and Factor MB demonstrate superior performance over EW and non-factor based models (GLASSO and MB). Furthermore, our method achieves lower estimation error in the combination weights (ref), which leads to lower risk of the combined forecast as shown in (ref). Interestingly, even though the precision matrix estimated using Factor MB has faster convergence rate in ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \cdot \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_2$ and ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \cdot \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_1$ norms as compared to Factor GLASSO, the weights estimated using Factor GLASSO converge faster. Also, note that the precision matrix estimated using the EW method also shows good convergence properties. However, in terms of estimating the combination weight, the performance of EW does not exhibit convergence properties. This is in agreement with previously reported findings (smith2009simple) that equal weights are not theoretically optimal, however, as demonstrated in the next subsection, the EW combination still leads to a relatively good performance in terms of MSFE although the FGM-based combinations outperform it.
We consider the standard forecasting model in the literature (e.g., Stock2002), which uses the factor structure of the high dimensional predictors. Suppose the data is generated from the following data generating process (DGP):
where $y_{t+1}$ is a univariate series of our interest in forecasting, ${\mathbf x}_t$ is an $N \times 1$ vector of regressors (predictors), ${\bm \beta}$ is an $N \times 1$ parameter vector, ${\mathbf g}_{t}$ is an $r \times 1$ vector of factors, ${\bm \Lambda}$ is an $N\times r$ matrix of factor loadings, ${\mathbf v}_t$ is an $N \times 1$ random vector following $\mathcal{N}(0,\sigma^{2}_{v})$, $\phi$ is an autoregressive parameter in the factors which is a scalar for simplicity, ${\bm \xi}_t$ is an $r \times 1$ random vector with each component independently following $\mathcal{N}(0,\sigma^{2}_{\xi})$, $\epsilon_{t+1}$ is a random error following $\mathcal{N}(0,\sigma^{2}_{\epsilon})$, and ${\bm \alpha}$ is an $r\times 1$ parameter vector which is drawn randomly from $\mathcal{N}(1,1)$. We set $\sigma_{\epsilon}=1$. The coefficients $\theta_{s}$ are set according to the rule
as in HANSEN2008. We set $c_1\in\{0,0.75\}$ and $c_2 \in \{0.6, 0.7, 0.8, 0.9\}$. We generate $r$ factors using (ref) with a grid of 10 different AR(1) coefficients $\phi$ equidistant between $0$ and $0.9$. To create ${\bm \Lambda}$ in (ref) we take the first $r$ rows of an upper triangular matrix from a Cholesky decomposition of the $N \times N$ Toeplitz matrix parameterized by $\rho$. We consider a grid of 10 different values of $\rho$ equidistant between $0$ and $0.9$.
One-step ahead forecasts are estimated from the factor-augmented autoregressive (FAR) models of orders $k,l$, denoted as FAR($k,l$):
where the factors $({\hat{g}_{1,t},\ldots,\hat{g}_{k,t}})$ are estimated from equation ((ref)). We consider the FAR models of various orders, with $k=1,\ldots,K$ and $l=1,\ldots,L$. We also consider the models without any lagged $y$ or any factors. Therefore, the total number of forecasting models is $p \equiv (1+K)\times(1+L)$, which includes the forecasting models using naive average or no factors.
The total number of observations is $T$, and the number of observations in the regression period (the train sample) is set to be the first half of the sample, $t=1,\ldots,m \equiv T/2$, to leave the second half of the sample, $t=m+1,\ldots,T$, for the out-of-sample evaluation (the test sample). We roll the estimation window over the test sample of the size $n\equiv T-m$, to update all the estimates in each point of time $t=1,\ldots,m$. Recall that $q$ denotes the number of factors in the forecast errors as in equation ((ref)). We first examine the properties of the combined forecasts based on the Factor Graphical models when $T$ and $p$ vary and compare their performance with the combined forecasts based on the GLASSO, MB and EW forecasts.
We consider a low-dimensional setup to demonstrate the advantage of using FGM even when the number of forecasts, $p$, is small relative to the sample size, $T$: (1) in such scenario EW has an advantage since there are not many models to combine and assigning equal weights should produce satisfactory performance, and (2) non-factor based models have the advantage over the models that estimate factors due to the estimation errors. As a result, this framework with the low-dimensional setup is favorable to EW and non-factor based models. Figure (ref) shows the MSFE for different sample sizes and fixed parameters: we report the results for two values of $c_1\in\{0,0.75\}$. As evidenced from Figure (ref), the models that use the factor structure outperform EW combination and non-factor based counterparts for both values of $c_1$. We see that Factor GLASSO, in general, has lower MSFE than Factor MB. This finding is further supported by our empirical application in Section 7.
In Appendix (ref) we examine the sensitivity of the competing models with respect to variation in the DGP parameters such as number of predictors $N$, values of $c_2$, $\phi$, the strength of factor loadings $\rho$, and the number of factors $q$. We conclude that Factor Graphical Models outperform equally-weighted combinations and the graphical models without factors.
An empirical application to forecasting macroeconomic time series in big data environment highlights the advantage of both Factor Graphical models described in Algorithms (ref)-(ref) in comparison with the existing methods of forecast combination. We use a large monthly frequency macroeconomic database of DataMcCracken, who provide a comprehensive description of the dataset and 128 macroeconomic series. We consider the time period 1960:01-2020:07 with the total number of observations $T=726$, the training sample consists of $m=120$ observations, and the test sample $n\equiv T-m-h+1$, where $h$ is the forecast horizon. We roll the estimation window over the test sample to update all the estimates in each point of time $t=m, \ldots,T-h$. We estimate $h$-step ahead forecasts from FAR($k,l$) which were defined in (ref) with $k=0, 1, \ldots, K=9$, and $l=0, 1,\ldots, L=11$. The total number of forecasting models is $p=120$. The optimal number of factors in the forecast errors (denoted as $q$ in equation (ref)) is chosen using the standard data-driven method that uses the information criterion IC1 described in Bai2002. We note that in the majority of the cases the optimal number of factors was estimated to be equal to 1.
Table (ref) compares the performance of the Factor GLASSO and Factor MB with the competitors for predicting seven representative macroeconomic indicators of the US economy: monthly industrial production (INDPRO), S&P500 composite index (S&P500), Consumer Price Index (CPIAUCSL), real personal consumption (DPCERA3MO86SBEA), M1 money stock (M1SL), civilian unemployment rate (UNRATE), and the effective federal funds rate (FEDFUNDS) using 127 remaining macroeconomic series. Let $\{Y_t\}_{t=1}^{T}$ be the series of interest for forecasting. Similarly to stevanovic2020machine, for INDPROD, S&P500, CPI, Real Personal Consumption and M1 Money Stock we forecast the average growth rate (with logs):
For UNRATE we forecast the average change (without logs):
And for FEDFUNDS we forecast the log of the series:
Table (ref) reports MSFEs of the competing methods with the smallest MSFE in each row in bold font. As evidenced from Table (ref), our methods outperform EW, GLASSO and nodewise regression: accounting for the factor structure results in lower MSFE. Therefore, the FGM framework developed in this paper leads to the superior performance of the combined forecast as compared to EW model even when the models/experts do not contain a lot of unique information. Our empirical application demonstrates that this finding does not originate from the difference in the performance of EW vs graphical models: as evidenced from Table (ref), the performance of GLASSO is worse than that of EW for the FEDFUNDS series, whereas Factor GLASSO outperforms EW. A similar pattern is observed in the performance of nodewise regression for M1 Money Stock. Therefore, the improvement in the combined forecast comes from incorporating the factor structure of the forecast errors into the graphical models. Note that in contrast with EW and non-factor based methods, the performance of Factor GLASSO and Factor MB does not deteriorate significantly when the forecast horizon, $h$, increases. Notice, however, that Factor Graphical Models tend to perform better for $h\geq 2$. In other words, accounting for common factors in forecast errors has greater benefit for longer horizons. Finally, for most series Factor GLASSO outperforms Factor MB, suggesting that for macroeconomic forecasting using weighted penalized log-likelihood and running $p$ coupled LASSO problems for estimating precision matrix is preferable to fitting $p$ separate LASSO regressions using each variable as the response and the others as predictors.
In this paper we overcome the challenge of using graphical models under the factor structure and provide a simple framework that allows practitioners to combine a large number of forecasts when experts tend to make common mistakes. Our new approach to forecast combinations breaks down forecast errors into common and unique parts which improves the accuracy of the combined forecast. The proposed algorithms, Factor Graphical Models, are shown to consistently estimate forecast combination weights and MSFE. Extensive simulations and empirical applications to macroeconomic forecasting in big data environment reveal that FGM outperforms equal-weighted forecasts and combined forecasts produced using graphical models without factors. With the superior performance observed at all forecast horizons, we find that the greater benefit from accounting for the common factors is evidenced at longer horizons.
\cleardoublepage \phantomsection \addcontentsline{toc}{section}{References} {14pt} \cleardoublepage \phantomsection \addcontentsline{toc}{section}{Figures}
\cleardoublepage
\cleardoublepage
\cleardoublepage