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.
82,010 characters · 12 sections · 66 citation commands
Combining Forecasts under Structural Breaks Using Graphical LASSO
\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 ECB SPF 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. Precision matrix represents a network of interacting entities, such as corporations or genes. When the data is Gaussian, the sparsity in the precision matrix encodes the conditional independence graph - two variables are conditionally independent given the rest if and only if the entry corresponding to these variables in the precision matrix is equal to zero. Graphical models are a powerful tool to directly estimate precision matrix, avoiding the step of obtaining an estimator of covariance matrix to be inverted. Prominent examples of graphical models include GL (GLASSO) and nodewise regression (meinshausen2006). Despite using different strategies for estimating precision matrix, all graphical models assume that precision matrix 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. fgl show that graphical models fail to recover the entries of a nonsparse precision matrix under the factor structure and propose FGL that combines the benefits of graphical models and factor models.
At the same time, the network of experts changes over time, that is, the relationships between forecasts produced by different experts or models can change either smoothly or abruptly (e.g., as a response to an unexpected policy shock, or in the times of economic downturns). Such changes give rise to different regimes and it is important to account for changes in optimal forecast combination weights induced by structural breaks. This paper augments fgl and develops a unified framework to generalize network inference in the presence of structural breaks. As a first extension, we model structural changes in factor loadings. As a second extension, we model structural changes in the precision matrix of the idiosyncratic component after removing common factors. We estimate regime-dependent precision matrix for forecast combination using both pre- and post-break data when forecast errors are driven by common factors. We call the proposed algorithm Regime-Dependent Factor Graphical LASSO (RD-FGL) and develop its scalable implementation using the Alternating Direction Method of Multipliers (ADMM).
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, to tackle changing relationships between forecasts produced by different experts or models as a response to unstable environments, we develop a unified framework to generalize network inference in the presence of structural breaks. We propose RD-FGL that models structural changes in factor loadings and idiosyncratic precision matrix. We develop scalable implementation of RD-FGL using ADMM to estimate regime-dependent forecast combination weights. Third, an empirical application to forecasting macroeconomic series using the data of the ECB SPF shows that incorporating (i) factor structure in the forecast errors together with (ii) sparsity in the precision matrix of the idiosyncratic components and (iii) regime-dependent combination weights improves the performance of a combined forecast over forecast combinations using equal weights.
We emphasize that in this paper our goal is to develop a framework for forecast combinations that incorporates structural breaks which have already occurred in the past. We neither consider the possibility of breaks over the forecast horizon (as explored in pesaran2006hierarchical), nor study the case of the out-of-sample breaks. These scenarios are interesting extensions of this paper, however they lie outside the scope of this paper. We also emphasize that we take the individual forecasts to be combined as given and do not discuss how the forecasts are generated.
The paper is structured as follows. Section 2 studies the approximate factor model for the forecast errors. Section 3 reviews FGL and contains theoretical results on the consistency of the FGL estimator for forecast combinations. Section 4 introduces Regime-Dependent graphical model and discusses its implementation using ADMM. Section 5 validates theoretical results using simulations. Section 6 studies an empirical application for macroeconomic time-series forecasting. Section 7 concludes.
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_{ij}}$ 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}}_{\text{max}}\equiv\max_{i,j}\@ifstar{\oldabs}{\oldabs*}{u_{ij}}$ denotes the element-wise maximum. For two sequences $a_{T,p}$ and $b_{T,p}$, we denote $a_{T,p}\asymp b_{T,p}$ if there exist constants $c_1, c_2 >0$ such that $c_1 a_{T,p} \leq b_{T,p} \leq c_2 a_{T,p}$.
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 ECB SPF 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.
Suppose we have $p$ competing forecasts of the univariate series $y_t$, $t=1,\ldots,T$ and $\widetilde{{\mathbf e}}_t=(\widetilde{e}_{1t},\ldots,\widetilde{e}_{pt})' \sim \mathcal{N} ({\mathbf m}, {\bm \Sigma})$ is a $p \times 1$ vector of forecast errors. Note that we allow bias in the forecasts. In fact, Figure (ref) demonstrates that the individual forecasts are indeed biased. Assume that the generating process for the forecast errors follows a $q$-factor model: $\widetilde{{\mathbf e}}_t = {\mathbf m} + {\mathbf B}{\mathbf f}_{t}+{\bm \varepsilon}_{t}$, 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. Define demeaned forecast errors as $ {\mathbf e}_t \equiv \widetilde{{\mathbf e}}_t - {\mathbf m}$ such that:
Unobservable factors, ${\mathbf f}_{t}$, and loadings, ${\mathbf B}$, are usually estimated by the principal component analysis (PCA), studied in 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}$, and $\operatorname*{\mathbb{E}}{{\mathbf e}_t{\mathbf e}'_t}={\bm \Sigma}={\mathbf B}{\bm \Sigma}_{f}{\mathbf B}'+ {\bm \Sigma}_{\varepsilon}$. 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 the constraint is 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 $\widehat{y}_{t}^{c}={\mathbf w}'\widehat{{\mathbf y}}_{t}$, 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 MSFE of the combined forecast error:
where ${\bm \iota}_p$ is a $p\times 1$ vector of ones.\footnote{As noted in timmerman_handbook_of_forecasting, global mean-variance optimization problem in finance is similar to the forecast combination problem. Combined forecasts can be viewed as the portfolio and the source of risk reflects incomplete information about the target variable and model misspecification possibly due to non-stationarities in the underlying data generating process.} 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}$. 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 $\@ifstar{\oldabs}{\oldabs*}{\frac{\text{MSFE}(\widehat{{\mathbf w}},\widehat{{\bm \Sigma}})}{\text{MSFE}({\mathbf w},{\bm \Sigma})} -1} = \@ifstar{\oldabs}{\oldabs*}{ \frac{\hat{a}^{-1} }{a^{-1}}-1}=\frac{\@ifstar{\oldabs}{\oldabs*}{a-\hat{a}}}{\@ifstar{\oldabs}{\oldabs*}{\hat{a}}}$ and $\@ifstar{\oldnorm}{\oldnorm*}{\widehat{{\mathbf w}}-{\mathbf w}}_1 = \Big[(a\widehat{{\bm \Theta}}{\bm \iota}_p) - (a{\bm \Theta}{\bm \iota}_p) + (a{\bm \Theta}{\bm \iota}_p) - (\hat{a}{\bm \Theta}{\bm \iota}_p) \Big]/p \cdot (\hat{a}a)$. 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 Section 3 and Theorem (ref).
Since our interest is in constructing weights for the forecast combination, our goal is to estimate a precision matrix of the forecast errors. This brings us to consider a family of graphical models, which have evolved from the connection between partial correlations and the entries of an adjacency matrix. The adjacency matrix has zero or one in its entries, with a zero entry indicating that two variables are independent conditional on the rest. The adjacency matrix is sometimes referred to as a “graph". 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. Weighted Graphical Lasso (GL) procedure (GLASSO) described in Supplemental Appendix (ref) is a representative member of graphical models family.
Before estimating precision matrix of forecast errors, ${\bm \Theta}$, we first obtain estimates of factors, $\widehat{{\mathbf f}}_t$, and factor loadings, $\widehat{{\mathbf B}}$, using PCA. Second, we obtain $\widehat{{\bm \Sigma}}_{f}=\frac{1}{T}\sum_{t=1}^{T}\widehat{{\mathbf f}}_{t}\widehat{{\mathbf f}}_{t}^{'}$, $\widehat{{\bm \Theta}}_{f}=\widehat{{\bm \Sigma}}_{f}^{-1}$, $\widehat{{\bm \varepsilon}}_t = {\mathbf e}_t-\widehat{{\mathbf B}}\widehat{{\mathbf f}_t}$, and $\widehat{{\bm \Sigma}}_{\varepsilon}=\frac{1}{T}\sum_{t=1}^{T}\widehat{{\bm \varepsilon}}_{t}\widehat{{\bm \varepsilon}}_{t}^{'}$. Third, we note that when common factors are present across the forecast errors, the precision of forecast errors, ${\bm \Theta}$, cannot be sparse because all pairs of the forecast errors are partially correlated given other forecast errors through the common factors. Therefore, instead of imposing sparsity assumption on ${\bm \Theta}$ we require sparsity of the precision matrix of the idiosyncratic errors, ${\bm \Theta}_{\varepsilon}$.
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. Hence, it is estimated with the Weighted Graphical Lasso penalty:
initialized with ${\mathbf W}_{\varepsilon}=\widehat{{\bm \Sigma}}_{\varepsilon}+\tau{\mathbf I}$, where $\widehat{\gamma}_{\varepsilon,ii}$ is the $(i,i)$-th element of $\widehat{{\bm \Gamma}}_{\varepsilon}^2\equiv \textup{diag}({\mathbf W}_{\varepsilon})$. The subscript $\tau$ in $\widehat{{\bm \Theta}}_{\varepsilon,\tau}$ means that the solution of the optimization problem in (ref) will depend upon the choice of the tuning parameter $\tau$. In order to simplify notation, we will omit the subscript $\tau$.
Finally we put all estimates together using the Sherman-Morrison-Woodbury formula to estimate the precision of forecast errors:
The aforementioned procedure introduced in fgl is called Factor Graphical LASSO (FGL).\footnote{Instead of decomposing covariance matrix into low-rank and idiosyncratic components, $\operatorname*{\mathbb{E}}{{\mathbf e}_t{\mathbf e}'_t}={\bm \Sigma}={\mathbf B}{\bm \Sigma}_{f}{\mathbf B}'+ {\bm \Sigma}_{\varepsilon}$, shi2020l2 regularize ${\bm \Sigma}$ directly with the factor structure (more broadly, latent group structure) in mind. FGL approach is different since it focuses on estimating precision matrix directly which is the main focus of the theoretical analysis in fgl. } It is summarized in Supplemental Appendix (ref), where we also discuss the choice of the tuning parameter $\tau$ in (ref).
We can use $\widehat{{\bm \Theta}}$ to estimate the forecast combination weights $ \widehat{{\mathbf w}}={\bm \Theta}{\bm \iota}_p/{\bm \iota}_p'\widehat{{\bm \Theta}}{\bm \iota}_p$. This 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.
Let us now examine the asymptotic properties of FGL. 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. 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 $ \alpha(T)=\sup_{A\in \mathcal{F}_{-\infty}^{0}, B \in \mathcal{F}_{T}^{\infty}}\@ifstar{\oldabs}{\oldabs*}{\Pr{A}\Pr{B}-\Pr{AB}}$.
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. Part (ii) of Assumption (ref) is needed for consistent estimation of the optimal forecast combination weights. 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. Assumptions (ref)(i), (ref), and (ref) are standard assumptions and are used in fan2013POET.
Let ${\bm \Lambda}_q = \text{diag}(\lambda_{1},\ldots,\lambda_q)$ be a matrix of $q$ leading eigenvalues of ${\bm \Sigma}$, and ${\mathbf V}_q = ({\mathbf v}_1,\ldots,{\mathbf v}_q)$ is a $p\times q$ matrix of their corresponding leading eigenvectors. Define $\widehat{{\bm \Sigma}}, \widehat{{\bm \Lambda}}_q,\widehat{{\mathbf V}}_q$ to be the estimators of ${\bm \Sigma},{\bm \Lambda}_q,{\mathbf V}_q$. We further let $\widehat{{\bm \Lambda}}_q=\text{diag}(\hat{\lambda}_1,\ldots,\hat{\lambda}_q)$ and $\widehat{{\mathbf V}}_q=(\hat{{\mathbf v}}_1,\ldots,\hat{{\mathbf 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{{\mathbf V}}_q\widehat{{\bm \Lambda}}_q\widehat{{\mathbf V}}_{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 assumption on the model is imposed:
Note that Assumptions (ref)-(ref) and (ref) are standard assumptions and are used in fan2018elliptical.
To study the properties of the combination weights and MSFE, we first need to establish the convergence properties of precision matrix produced by Algorithm (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 FGL. Note that $d_T$ can be smaller than or equal to $s_T$.
Let $\varrho_{T}$ be a sequence of positive-valued random variables such that $\varrho_{T}^{-1}\omega_{T}\xrightarrow{P}0$ and $\varrho_{T}d_Ts_T\xrightarrow{P}0$, with $\tau \asymp \omega_{T}$ (where $\tau$ is the tuning parameter for the FGL in (ref)). fgl 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_{T}d_Ts_T) = o_P(1)$ and ${\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}_{2}=\mathcal{O}_P(\varrho_{T}s_T) = o_P(1)$ for FGL.
Having established the convergence rates for precision matrix, we now study the properties of the combination weights and the resulted MSFE.
The proof of Theorem (ref) can be found in Supplementary Appendix (ref). Note that the rates of convergence for MSFE and precision matrix ${\bm \Theta}$, which was derived in fgl, are the same and both are faster than the combination weight rates. In contrast to the classical graphical model in Algorithm (ref), the convergence properties of which were examined by Sara2018 among others, the rates in Theorem (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.
There are two streams of literature that study time-varying networks. The first one models dynamics in the precision matrix locally. Zhou2010 develop a nonparametric method for estimating time-varying graphical structure for multivariate Gaussian distributions using an $\ell_1$-penalized log-likelihood. They find out that if the covariances change smoothly over time, the covariance matrix can be estimated well in terms of predictive risk even in high-dimensional problems. Kolar2018 introduce nonparanormal graphical models that allow to model high-dimensional heavy-tailed systems and the evolution of their network structure. They show that the estimator consistently estimates the latent inverse Pearson correlation matrix. The second stream of literature allows the network to vary with time by introducing two different frequencies. tvgl study time-varying Graphical LASSO with smoothing evolutionary penalty.
We augment the framework in Section 3 to account for regime switching by modeling the change in precision matrix due to $N$ structural breaks. Define $n_j \equiv t_j - t_{j-1}$ to be the sample between the $j$-th and $(j-1)$-th break points, where $j=1,\ldots, N+1$, $\sum_{j=1}^{N+1}n_j = T$, $t_0 = 0$, $N\leq T$.
Macroeconomic and financial datasets typically span a long time period, hence, the assumption of time-invariant factor loadings is restrictive. As a first extension to FGL, we model structural changes in factor loadings using a framework similar to su2017time. For now assume a single known break $N=1$ which occurs at $T_1$.\footnote{The possible presence of break over the forecast horizon, as explored in pesaran2006hierarchical, would result in common forecast errors of all forecasting models, which can be captured by using the factor model. This would be another scenario of common forecast errors as motivated by Figure (ref).} Write equation (ref) as:
To estimate $\{{\mathbf b}_{i}\}_{i=1}^{p}$ and $\{{\mathbf f}_t\}_{t=1}^{T}$, we can consider the following weighted least squares problem:
subject to certain identification restrictions to be specified later on. Here, $K_{\gamma t} = \gamma \mathds{1}\left[t \leq T_1\right] + \mathds{1}\left[t > T_1\right]$ is a discrete kernel as in li2013categorical with $\gamma \in [0,1]$. Since more recent information is usually more relevant to forecasting, such kernel-weight estimator gives weight 1 to post-break observations and weight $\gamma$ to pre-break observations.\footnote{Adjusting the parameter to control the degree to which pre-break data are discounted by the model is useful for applied researchers. Bayesian framework with hierarchical priors provides an alternative way to combine the sample information contained in the objective function with prior information about the values of the model parameters and the relations among them (see pastor2001hierarchical,pesaran2006hierarchical).}
Define the $T \times p$ matrices ${\mathbf E}(\gamma) = \Big({\mathbf e}_1(\gamma),\ldots, {\mathbf e}_p(\gamma) \Big)$, $\mathcal{E}(\gamma) = \Big({\bm \varepsilon}_1(\gamma),\ldots, {\bm \varepsilon}_p(\gamma) \Big)$, where\\ ${\mathbf e}_{i}(\gamma) = \Big(K_{\gamma 1}^{1/2}e_{i1},\ldots, K_{\gamma T}^{1/2}e_{iT}\Big)'$ and ${\bm \varepsilon}_{i}(\gamma) = \Big(K_{\gamma 1}^{1/2}\varepsilon_{i1},\ldots, K_{\gamma T}^{1/2}\varepsilon_{iT}\Big)'$. Also, let\\ ${\mathbf F}(\gamma) = \Big(K_{\gamma 1}^{1/2}{\mathbf f}_{1},\ldots, K_{\gamma T}^{1/2}{\mathbf f}_{T}\Big)'$ be a $T \times q$ matrix collecting factors. In matrix notation, the transformed model in (ref) can be written as ${\mathbf E}(\gamma) = {\mathbf F}(\gamma){\mathbf B}' + \mathcal{E}(\gamma)$, where ${\mathbf B}= ({\mathbf b}_{1},\ldots, {\mathbf b}_{p})'$ is a $p \times q$ matrix of factor loadings.
As shown in su2017time for the continuous kernel, the minimization problem in (ref) reduces to:
The problem in (ref) is the conventional PCA problem. The estimated factor matrix $\widehat{{\mathbf F}}(\gamma) = \Big(K_{\gamma 1}^{1/2}\widehat{{\mathbf f}}_{1},\ldots, K_{\gamma T}^{1/2}\widehat{{\mathbf f}}_{T}\Big)'$ is $\sqrt{T}$ times eigenvectors corresponding to the $q$ largest eigenvalues of ${\mathbf E}(\gamma){\mathbf E}'(\gamma)$, arranged in descending order, and $\widehat{{\mathbf B}}'(\gamma) = (\widehat{{\mathbf F}}(\gamma)\widehat{{\mathbf F}}'(\gamma))^{-1}\widehat{{\mathbf F}}'(\gamma){\mathbf E}(\gamma) = \widehat{{\mathbf F}}'(\gamma) {\mathbf E}(\gamma)/T $ are the estimators of the corresponding time-varying factor loadings, where $\widehat{{\mathbf B}}(\gamma) = (\widehat{{\mathbf b}}_{1}(\gamma),\ldots, \widehat{{\mathbf b}}_{p}(\gamma))'$ is $p \times q$.
Since the estimator $\widehat{{\mathbf F}}(\gamma)$ is only consistent up to a rotation, we use a two-stage estimation procedure to obtain a consistent estimator (su2017time). Based on the consistent estimators of ${\mathbf b}_{i}$'s obtained from the first stage, consistent estimators of ${\mathbf f}_{t}(\gamma)$ can be obtained by considering the following least squares problem $\widehat{{\mathbf f}}_t(\gamma) = \arg\!\min_{{\mathbf f}_t} \sum_{i=1}^{p} \Big[e_{it} - \widehat{{\mathbf b}}^{'}_{i}(\gamma){\mathbf f}_{t} \Big]^2$ which yields the solution $\widehat{{\mathbf f}}_t(\gamma) = \Big( \sum_{i=1}^{p} \widehat{{\mathbf b}}_{i}(\gamma) \widehat{{\mathbf b}}^{'}_{i}(\gamma) \Big)^{-1}\Big( \sum_{i=1}^{p} \widehat{{\mathbf b}}_{i}(\gamma) e_{it} \Big)$.
As in su2017time, we assume that $\operatorname*{\mathbb{E}}{{\mathbf f}_t{\mathbf f}'_t}$ is homogeneous over $t$. This assumption is not restrictive, since if $\operatorname*{\mathbb{E}}{{\mathbf f}_t{\mathbf f}'_t} = {\bm \Sigma}_{f,t}$, we can rewrite the common component as ${\mathbf b}'_{i}{\mathbf f}_t = \Big( {\bm \Sigma}_{f}^{-1/2}{\bm \Sigma}_{f,t}^{1/2}{\mathbf b}_{i} \Big)' {\bm \Sigma}_{f}^{1/2} {\bm \Sigma}_{f,t}^{-1/2}{\mathbf f}_{t} = {\mathbf b}^{*'}_{i}{\mathbf f}^{*'}_{t}$, where ${\mathbf b}^{*}_{i} = {\bm \Sigma}_{f}^{-1/2}{\bm \Sigma}_{f,t}^{1/2}{\mathbf b}_{i} $, and ${\mathbf f}^{*}_{t} = {\bm \Sigma}_{f}^{1/2} {\bm \Sigma}_{f,t}^{-1/2}{\mathbf f}_{t}$ satisfies $\operatorname*{\mathbb{E}}{{\mathbf f}^{*}_{t}{\mathbf f}^{*'}_{t}} ={\bm \Sigma}_{f} $ for each $t$.
To choose the optimal tuning parameter $\gamma$ in (ref), we use the cross-validation and solve the following minimization problem:
where $\widehat{{\mathbf b}}^{(-s)}_{i}(\gamma)$ and $\widehat{{\mathbf f}}^{(-s)}_{s}(\gamma)$ are estimated by leaving the $s$-th time series observation out of the PCA procedure.
As a second extension to FGL, we model structural changes in the precision matrix of the idiosyncratic component. Let ${\bm \Sigma}_{\varepsilon,j}$ and ${\bm \Sigma}_j$ be covariance matrices of idiosyncratic part and forecast errors in regime $j$. Define the corresponding precision matrices to be ${\bm \Theta}_{\varepsilon,j}\equiv{\bm \Sigma}_{\varepsilon,j}^{-1}$ and ${\bm \Theta}_j \equiv {\bm \Sigma}_{j}^{-1}$. Similarly to the previous subsection, without loss of generality we assume ${\bm \Sigma}_{f_j}={\bm \Sigma}_f$ for all regimes $j$.
Let $\widehat{{\bm \Sigma}}_{\varepsilon,j}=\frac{1}{n_j}\sum_{k=1}^{n_j}\widehat{{\bm \varepsilon}}_{j,k}\widehat{{\bm \varepsilon}}_{j,k}'$. To model dynamics in $\{{\bm \Theta}_{\varepsilon,j}\}_{j=1}^{N+1}$ we use the following optimization problem:
where the penalty for the off-diagonal (od) elements is $\@ifstar{\oldnorm}{\oldnorm*}{{\bm \Theta}_{\varepsilon,j}}_{\text{od},1} = \sum_{l\neq q}\widehat{\gamma}_{\varepsilon,ll,j}\widehat{\gamma}_{\varepsilon,qq,j}\@ifstar{\oldabs}{\oldabs*}{\theta_{\varepsilon,lq,j}}$, $\widehat{\gamma}_{\varepsilon,ll,j}$ is the $(l,l)$-th element of $\widehat{{\bm \Gamma}}_{\varepsilon,j}^2\equiv \textup{diag}(\widehat{{\bm \Sigma}}_{\varepsilon,j})$ and $\theta_{\varepsilon,lq,j}$ is the $lq$-th element of matrix ${\bm \Theta}_{\varepsilon,j}$. Figure (ref) visualizes dynamics of the precision matrix.
The optimization problem in (ref) has two tuning parameters: $\alpha$, which determines the sparsity level of the network, and $\beta$, which controls the strength of resemblance between two neighboring precision estimators.\footnote{pesaran2006hierarchical use Bayesian framework to handle parameter instability using hierarchical priors. As pointed out in their paper, “intuition for the use of hierarchical priors comes from the shrinkage literature, since one can think of the parameters within the individual regimes as being shrunk towards a set of the so-called hyperparameters that characterize the “top” layer of the hierarchy" (p. 1059).} In simulations and the empirical application we use the following procedure for tuning $\alpha$ and $\beta$: first, we set a grid of values $(\alpha,\beta) \in \{0,0.25,0.5,1,10,30\}$. Second, we use the first 2/3 of the training data to estimate forecast combination weights and jointly tune $\alpha$ and $\beta$ in the remaining 1/3 to yield the smallest value of the objective function, which is chosen to be either ${\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$-loss of precision matrix for simulations in Subsection 5.1, or MSFE for simulations in Subsection 5.2 and the empirical application. Note that when $\beta=0$, the optimization in (ref) reduces to estimating ${\bm \Theta}_{\varepsilon,i}$ using Algorithm (ref) in each regime separately. Naturally, this incorporates the case when the structural break is strong and only the post-break data is used for producing forecast combination weights. When $\beta$ is large, there are weak structural breaks in ${\bm \Theta}_{\varepsilon,j}$, and ${\bm \Theta}_{\varepsilon,j}$'s are estimated by using the data across different regimes. Section 6 provides more discussion on this in the context of our empirical application.
The smoothing function $\psi(\cdot)$ in (ref) can be LASSO ($\psi=\sum_{l,q}\@ifstar{\oldabs}{\oldabs*}{\cdot}$), Group LASSO ($\psi=\sum_{q}\@ifstar{\oldnorm}{\oldnorm*}{\cdot_q}_2$), or Ridge ($\psi=\sum_{l,q}(\cdot_{lq})^{2}$). LASSO penalty encourages small changes in the precision matrix over time: when the $lq$-th element changes at two consecutive times, the penalty forces the rest of the elements of the precision to remain the same. Group LASSO penalty allows the entire graph to restructure at some time points. This penalty is useful for anomaly detection, since it can identify structural changes in the network structure. Ridge penalty allows the network to change smoothly over time. This penalty is less strict than the LASSO penalty: instead of encouraging the graphs to be exactly the same, it allows smooth transitions. In our empirical application we use Ridge penalty to accommodate smooth transitions of precision over time.
To estimate (ref) we use the ADMM algorithm described in details in Supplementary Appendix (ref). Once ${\bm \Theta}_{\varepsilon,i}$ is estimated, we combine estimated factors, loadings and precision matrix of the idiosyncratic components using Sherman-Morrison-Woodbury formula to estimate the final precision matrix of forecast errors and use it to compute optimal forecast combination weights. We call the aforementioned procedure RD-FGL and summarize it in Algorithm (ref).
We develop a scalable implementation of (ref) for the RD-FGL in Algorithm (ref) through ADMM, which is extensively discussed in Supplementary Appendix (ref). ADMM is a distributed convex optimization approach (BoydProximalAlgorithms) that allows us to split the optimization problem in (ref) into a series of subproblems. As pointed out in tvgl, the scalability of ADMM comes from the improved runtime: to estimate a $p \times p$ matrix, the cost per iteration of ADMM is $\mathcal{O}(p^3)$ (which is the cost of an eigendecomposition of the ${\bm \Theta}$ step in Supplemental Appendix (ref)). In contrast, the runtime of general interior-point methods is $\mathcal{O}(p^6)$ (mohan2014node).
The previous two subsections assumed that the number and location of breaks are known. We now relax these assumptions. First, assume that the number of breaks in factor loadings, $N_B$, and the number of breaks in idiosyncratic precision, $N_{{\bm \Theta}}$, are known and $N_B=N_{{\bm \Theta}}=1$, but their locations are unknown and might differ from each other.
To estimate the location of the break in factor loadings, we adapt the procedure in bai2020breakinloadings. For a given break point in loadings, $T_{1}$, define the sum of squared residuals (SSR) as in (ref):
where $K_{\gamma t} = \gamma \mathds{1}\left[t \leq T_{1}\right] + \mathds{1}\left[t > T_{1}\right]$ is a discrete kernel. The estimated break date is given by $\widehat{T}_{1}=\arg\!\min_{1\leq T_{1}\leq T-1} \text{SSR}(T_{1})$.
To estimate the location of the break in ${\bm \Theta}_{\varepsilon}$ we use the procedure similar to bai2010breakinprecision. Define $t_{1}$ to be a break point in ${\bm \Theta}_{\varepsilon}$. Recall, $n_j = t_j - t_{j-1}$, where $j=1,2$. Note that the number of observations in each regime depends on $t_1$: $n_1 = t_1 - t_{0}$ and $n_2 = t_2 - t_{1}$. For a given break point in idiosyncratic precision $t_1$, define the following objective function as in (ref):
The estimated break date is given by $\widehat{t}_{1}=\arg\!\min_{1\leq t_{1}\leq T-1} L(t_{1})$.\footnote{As noted by bai1998estimating, it is difficult to detect the break near the end of the sample. It is reasonable to expect that if break magnitude is reduced, the procedures in (ref) and (ref) might not detect the break or detect it with a delay. smith2021break point out that in such situations only few observations from the current regime are available to estimate the model parameters, leading to volatile and inaccurate forecasts. To address slow detection of breaks, smith2021break exploit information in the cross-section to detect breaks more rapidly in real time. To address small sample size problem, smith2021break adopt a Bayesian approach that uses economically motivated priors to shrink the parameters towards sensible values that rule out economically implausible values. As we discussed in Footnote 4, the Bayesian framework is alternative to using kernel-weighted observations.}
When the number of breaks in either loadings is known and greater than 1 ($N_B>1$) and/or $N_{{\bm \Theta}}>1$, we can use the one-at-a-time approach as in bai2010breakinprecision: the objective functions are identical to (ref) and (ref). The breaks are estimated sequentially. Once the first break is obtained, we split the sample at the estimated break point, resulting in two subsamples. A single break point in each subsample is estimated, but only one that achieves the smallest objective function ((ref) or (ref)) is retained. If the number of breaks is equal to two, the procedure is stopped. Otherwise, we continue splitting into subsamples until all breaks are estimated.
If the number of breaks is unknown, we proceed as suggested in bai2010breakinprecision: in the aforementioned one-at-a-time approach apply the test for existence of break point (bai2003multiplebreak) to each subsample before estimating a break point.
We divide the simulation results into two subsections. In the first subsection we study the consistency of the FGL and RD-FGL for estimating precision matrix and the combination weights. In the second subsection we evaluate the out-of-sample forecasting performance of combined forecasts in terms of MSFE. We compare the performance of forecast combinations based on the factor models in Algorithms (ref), (ref) with equal-weighted (EW) forecast combination\footnote{As pointed out by the referee, EW arises when the forecast errors follow a factor structure (one factor, homogeneous idiosyncratic variance). It can be viewed as one of “factor-based" methods.}, and combinations that use GL without factor structure (Algorithm (ref)). We examine the performance of RD-FGL for different specifications of the smoothing function $\psi(\cdot)$ as described in Subsection 4.2. LASSO penalty is denoted as $\ell_1$, Group LASSO as $\ell_g$, and Ridge as $\ell_2$. Similarly to the literature on graphical models, all exercises use 100 Monte Carlo simulations. We present simulation results with a structural break in both ${\mathbf B}$ and ${\bm \Theta}_{\varepsilon}$. The results without a break, with break only in ${\bm \Theta}_{\varepsilon}$, and with multiple breaks can be found in Supplemental Appendix (ref).
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, \ i=1,\ldots, p$. Let $\widehat{{\bm \Theta}}$ be the precision matrix estimator. We show consistency of the FGL 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}$, and (ii) 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)_{lm}$, where $(b)_{lm}=\rho^{\@ifstar{\oldabs}{\oldabs*}{l-m}}$, $l,m\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 $\operatorname*{\mathbb{E}}{{\mathbf e}_t{\mathbf e}'_t}={\bm \Sigma}={\mathbf B}{\bm \Sigma}_{f}{\mathbf B}'+ {\bm \Sigma}_{\varepsilon}$. 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 with $a_{\varepsilon,lm}$ being the $l,m$-th element of the adjacency matrix ${\mathbf A}_{\varepsilon}$. We set $a_{\varepsilon,lm} = a_{\varepsilon,ml}=1, \ \text{for} \ l\neq m$ 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 = 500/(pT^{0.8})$, which makes $s_T = \mathcal{O}(T^{0.05})$. The adjacency matrix has all diagonal elements equal to zero. To generate a sparse symmetric positive-definite precision matrix we use Scikit-Learn datasets package in Python (scikit). To control the magnitude of partial correlations, the value of the smallest coefficient is set to 0.1 and the value of the largest coefficient is set to 0.3.
To incorporate structural breaks in ${\bm \Theta}_{\varepsilon}$ and factor loadings ${\mathbf B}$, we proceed as follows. We fix a single break point in the middle of the sample size, $T/2$: in the precision matrix of the idiosyncratic errors before the break, referred to as ${\bm \Theta}_{\varepsilon,1}$, the value of the largest coefficient is set to 0.4; whereas in the precision matrix of the idiosyncratic errors after the break, ${\bm \Theta}_{\varepsilon,2}$, the value of the largest coefficient is set to 0.6. As a consequence, even though both matrices are still sparse, ${\bm \Theta}_{\varepsilon,2}$ has larger partial correlations. We use ${\bm \Theta}_{\varepsilon,1}$ and ${\bm \Theta}_{\varepsilon,2}$ to generate ${\bm \varepsilon}_t$ in (ref). For the structural break in factor loadings (which is assumed to happen at the same time as the structural change in ${\bm \Theta}_{\varepsilon}$), before the break we set $\rho_1=0.2$ in the Toeplitz matrix used to generate ${\mathbf B}$ (i.e., ${\mathbf B} = (b)_{lm}$, where $(b)_{lm}=\rho^{\@ifstar{\oldabs}{\oldabs*}{l-m}}$), and after the break we set $\rho_2=0.6$.
Figure (ref) shows 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).
Figure (ref) examines the performance when there are breaks in both ${\bm \Theta}_{\varepsilon}$ and ${\mathbf B}$: accounting for the break significantly reduces the estimation error of precision matrix and combination weights. We report the results for the case when $\gamma$ is estimated using cross-validation ($\gamma = \hat{\gamma}$) (as discussed in Section 4). Supplemental Appendix (ref) presents the results for the case when the break is only in ${\bm \Theta}_{\varepsilon}$.
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 $M \times 1$ vector of regressors (predictors), ${\bm \alpha}$ is an $M \times 1$ parameter vector, ${\mathbf g}_{t}$ is an $r \times 1$ vector of factors, ${\bm \Lambda}$ is an $M\times r$ matrix of factor loadings, ${\mathbf v}_t$ is an $M \times 1$ random vector following $\mathcal{N}(0,\sigma^{2}_{v}{\mathbf I}_M)$, $\phi$ is an autoregressive parameter in the factors which is a scalar for simplicity, ${\bm \xi}_t$ is an $M \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 $\theta_{s}=(1+s)^{c_1}c_{2}^{s}$ as in HANSEN2008. We set $c_1=0.75$. We set $M=100$ and generate $r=5$ factors. To create ${\bm \Lambda}$ in (ref) we take the first $r$ rows of an upper triangular matrix from a Cholesky decomposition of the $M \times M$ Toeplitz matrix parameterized by $\rho = 0.9$. The ranking of competing models was not very sensitive to varying values of $\phi$, $\rho$, $c_2$, and $r$.
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. We set $K=2$ and $L=7$.
The total number of observations is $m$. The period for training the models is set to be $m_1 = m/2$ -- this is used to train competing FAR models in (ref). The remaining part of the sample, $m_2 = m-m_1$ is split as follows: the estimation window for training competing models (that is, EW, GL, FGL, and RD-FGL) is set to be of size $= m_2/2$. We roll the estimation window over the the test sample of the size $m_2/2$ to update all the estimates in each point of time. Recall that $q$ denotes the number of factors in the forecast errors as in equation ((ref)).
To incorporate structural break we proceed as follows. The period for training the models is set to be $m_1 = m/3$ -- this is used to train competing FAR models in (ref). The remaining part of the sample, $m_2 = m-m_1$ is split as follows: the estimation window for training competing models is set to be of size $= m_2/2$. We roll the estimation window over the test sample of the size $m_2/2$. The break point is fixed at 1/2 of the first estimation window. Before the break, when generating $\theta_{s}$ we set $c_2 = 0.3$, and after the break $c_2 = 0.9$. All other parameters stay unchanged. Notice that the break in $c_2$ can propagate into both a break in precision matrix and factor loadings.
Similarly to the previous subsection, we include different specifications of the smoothing function $\psi(\cdot)$. Figure (ref) shows the performance of all models including RD-FGL with $\gamma$ estimated using cross-validation: similarly to the conclusions in the previous subsection, accounting for the break significantly reduces MSFE of the combined forecast.
We use quarterly forecasts on the expected rates of inflation, real GDP growth and unemployment rate in the Euro area published by the \href{https://www.ecb.europa.eu/stats/ecb_surveys/survey_of_professional_forecasters/html/index.en.html}{ECB}. The raw data records 119 forecasters in total, but the panel is highly unbalanced with many missing values due to entry and exit in the long span. We follow Shi_ForecastCombinations to obtain most qualified forecasters: first, we filter out irregular respondents if they missed more than 45% of the observations; second, we use a random forest imputation algorithm (MissForest-R,MissForest-Paper) to interpolate the remaining missing values. We consider the forecasts of three main economic indicators: (1) Real GDP growth defined as the year-on-year (YoY) percentage change of real GDP, based on standardized European System of National and Regional Accounts (ESA) definition. The time period under consideration is 1999:Q3-2023Q1 (which yields the total number of observations equal to 95), the final number of forecasters is $p=59$, and the prediction horizon is 2-quarters ahead. (2) Inflation which is defined as the YoY percentage change of the Harmonised Index of Consumer Prices (HICP) published by Eurostat. The time period under consideration is 1999:Q4-2023Q1 (which yields the total number of observations equal to 95), the final number of forecasters is $p=56$, and the prediction horizon is 2-quarters ahead. (3) Unemployment rate which refers to Eurostat's definition and it is calculated as percentage of the labor force. The time period under consideration is 1999:Q3-2022Q4 (which yields the total number of observations equal to 94), the final number of forecasters is $p=45$, and the prediction horizon is 2-quarters ahead.
We consider four choices of the training sample: $R\in\{20,30,40,50\}$, the estimation window is rolled over the test sample to update the estimates in each point of time. 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. In the majority of the cases the optimal number of factors was estimated to be equal to 1. To explore the benefits of using FGL and RD-FGL for forecast error quantification, we consider several alternative estimators of covariance/precision matrix of the idiosyncratic component in (ref): (1) linear shrinkage estimator of covariance developed by Ledoit2004 further referred to as Factor LW (FLW); (2) nonlinear shrinkage estimator of covariance by ledoit2017nonlinear (Factor NLW or FNLW); (3) POET (fan2013POET); (4) constrained $\ell_{1}$-minimization for inverse matrix estimator, CLIME (cai2011constrained) (Factor CLIME or FCLIME); (5) nodewise regression developed by meinshausen2006 (Factor MB or FMB). To examine the benefits of imposing sparsity on ${\bm \Theta}_{\varepsilon}$ we also include the factor model without sparsity assumption on the idiosyncratic error precision matrix (referred to as Not Sparse) -- this corresponds to imposing $\tau=0$ in (ref). To examine the benefits of using factor structure, we include several counterparts of the aforementioned models that directly estimate precision of the forecast errors without estimating factors and loadings: GL, LW, NLW, CLIME, and MB. For RD-FGL, similarly to the simulations, we include different specifications of the smoothing function $\psi(\cdot)$ and report the results for the case when the break parameter in factor loadings is estimated using cross-valiation ($\gamma=\hat{\gamma}$).
Our benchmark is the simple average with equal weights on all forecasters (referred to as EW). Going back to the discussion in Section 4 regarding setting $\beta=0$ in equation (ref): as we pointed out, this corresponds to using only post-break sample for estimation which is suboptimal since the value of $\beta$ is already chosen optimally from the grid that includes $\beta = 0$ to minimize the MSFE. Hence, by construction, RD-FGL is superior to using only post-break data.
For RD-FGL the number of breaks for loadings and precision is estimated using the test for existence of break point (bai2003multiplebreak): using their sequential procedure we search for up to three breaks and set the trimming parameter to 10% of the total number of observations, and the significance level at 5%. The location of the break points for each series is estimated using the one-at-a-time approach described in Subsection 4.3.
Table (ref) compares the performance of FGL and RD-FGL with the competitors for predicting three macroeconomic indicators for Euro-area using a combination of ECB SPF forecasts. It reports the ratios of MSFE of each method to the MSFE of the EW combined forecast. Using the Model Confidence Set (MCS) of hansenMCS, we identify the set of superior models (SSM) for each series and horizon at 90% confidence level. Once SSM is identified, we rank these models according to the relative sample loss of the $i$-th model relative to the average across models in SSM, and report the ranking in Table (ref).
There are three main findings that we learn from analyzing Table (ref): (1) for all series factor-based models outperform non-factor ones. This means that incorporating the factor structure in the forecast errors improves forecasting performance. (2) for all series the Not Sparse model provides one of the worst performances. This means that the factor structure per se is not sufficient to achieve performance gains over EW, hence, it is necessary to impose sparsity on the precision matrix of the idiosyncratic components. (3) For real GDP growth and inflation series RD-FGL is always included in the SSM. For the unemployment rate, FGL outperforms RD-FGL. This result is supported by the behavior observed in the actual series: real GDP growth and inflation exhibit strong breaks following the global financial crisis and Covid pandemic, however this is not the case for the unemployment rate series that did not have strong breaks throughout the whole sample period.
In this paper we develop a unified framework to generalize network inference under a factor structure in the presence of structural breaks. We overcome the challenge of using graphical models under the factor structure and provide a simple approach that allows practitioners to combine a large number of forecasts when experts tend to make common mistakes. Using pre- and post-break data, our new approach to forecast combinations breaks down forecast errors into common and unique parts which improves the accuracy of the combined forecast. We allow the structural breaks to affect factor loadings and idiosyncratic precision matrix. For the ease of practical use we develop a scalable optimization procedure for RD-FGL, based on the ADMM. The empirical application to forecasting macroeconomic series using the data of the ECB Survey of Professional Forecasters shows that incorporating (i) factor structure in the forecast errors together with (ii) sparsity in the precision matrix of the idiosyncratic components and (iii) regime-dependent combination weights improves the performance of a combined forecast.
\cleardoublepage
\phantomsection
\addcontentsline{toc}{section}{References}
{14pt}
\cleardoublepage